text:
# References

DARK ENERGY, DARK MATTER AND
THE CHAPLYGIN GAS
R.Colistete Jr.111e-mail: coliste@ccr.jussieu.fr,J. C. Fabris222e-mail: fabris@cce.ufes.br, S.V.B. Gonçalves333e-mail: sergio@cce.ufes.br and P.E. de Souza444e-mail: patricia.ilus@bol.com.br
Departamento de Física, Universidade Federal do Espírito Santo, CEP29060-900, Vitória, Espírito Santo, Brazil

## Abstract

The possibility that the dark energy may be described by the Chaplygin gas is discussed. Some observational constraints are established. These observational constraints indicate that a unified model for dark energy and dark matter through the employement of the Chaplygin gas is favored.
PACS number(s): 98.80.Bp, 98.65.Dx

The possible existence of dark matter and dark energy in the Universe constitutes one of the most important problems in physics today. One calls dark matter the mysterious component of galaxies and clusters of galaxies that manifestates itself only through gravitational effects, affecting only the dynamics of those object. Dark matter appears in collapsed objects. The existence of dark energy, on the other hand, has been speculated due to the observation of the Universe as a whole. In particular, the position of the first accoustic peak in the spectrum of the anisotropy of the cosmic microwave background radiation and the data from supernova type Ia suggest that the dynamics of the Universe is such that the most important matter component is non-clustered, exhibiting negative pressure. The current observations favor $`\mathrm{\Omega }_{dm}0.3`$ and $`\mathrm{\Omega }_{de}0.7`$ for the proportion of dark matter and dark energy.

The most natural candidate to represent dark energy is a cosmological constant. However, it is necessary a fine tunning of $`120`$ orders of magnitude in order to obtain agreement with observations. Another popular candidate today is quintessence, a self-interacting scalar field. But the quintessence program suffers from fine tuning of microphysical paremeters. In this work, we discuss the possibility that the dark energy may be represented by the Chaplygin gas, which is characterized by the equation of state

$$p=\frac{A}{\rho }.$$
(1)
The Chaplygin gas has an interesting motivation connected with string theory. In fact, if we consider a $`d`$-brane configuration in the $`d+2`$ Nambu-Goto action, the employement of the light-cone parametrization leads to the action of a newtonian fluid with the equation of state (1), whose symmetries are the same as the Poincaré group. Hence, the relativistic character of the action is somehow hidden in the equation of state (1). For a review of the properties of the Chaplygin gas see reference .

Using the relativistic equation of conservation for a fluid in a FRW background, we obtain that the density of the Chaplygin gas depends on the scale factor as

$$\rho =\sqrt{A+\frac{B}{a^6}}.$$
(2)
Hence, initially the Chaplygin gas behaves as a pressurelless fluid and assumes asymptotically a behaviour typical of a cosmological constant. In this sense, the Chaplygin gas may be an alternative for the description of dark energy. Remark that, even if the pressure associated with the Chaplygin gas is negative, the sound velocity is real. Hence, no instability problem occurs as it happens with other fluids with negative pressure.

In reference  the density perturbations in a Universe dominated by the Chaplygin gas was investigate. In that work, the fact that the Chaplygin gas may be represented by a newtonian fluid was explored. The density contrast behaves as

$$\delta =t^{1/6}[C_1J_\nu (\mathrm{\Sigma }^2nt^{7/3})+C_2J_\nu (\mathrm{\Sigma }^2nt^{7/3})],$$
(3)
where $`n`$ is the wavenumber of the perturbations, $`\nu =5/14`$ and $`\mathrm{\Sigma }^2=\frac{54}{49}\frac{A\pi G}{a_0^2}`$, $`a_0`$ is a constant and $`C_1`$ and $`C_2`$ are constants that depend on $`n`$. A simple analysis of this solution indicates that initially the perturbations behave as in the case of the pressurelless fluid, presenting asymptotically oscillatory behaviour with decreasing amplitude, approaching a zero value typical of a cosmological constant. Even if this result was obtained in a newtonian context, a numerical analysis of the corresponding relativistic equations exhibit the same behaviour.

In reference , the mass power spectrum of the perturbations in a Universe dominated by the Chaplygin gas was computed. The matter content was assumed to be composed of radiation, dark matter and dark energy, the latter represented by the Chaplygin gas. Due to the complexity of perturbed equations, a numerical analysis was performed. The initial spectrum, at the moment of decoupling of radiation and matter, was supposed to be the Harrisson-Zeldovich scale invariant spectrum. Taken at constant time, this implies to impose an initial spectral index $`n_s=5`$. The system evolves and the power spectrum is computed for the present time. The final results were obtained for a Universe dominated uniquely by a pressurelless fluid (baryonic model), for a cosmological constant model and for a mixed of dark matter and Chaplygin gas, with different values for the sound velocity of the Chaplygin gas, defined as $`v_s^2=\overline{A}=A/\rho _{c0}^2`$, where $`\rho _{c0}`$ is the Chaplygin gas density today. When $`\overline{A}=1`$ the Chaplygin gas becomes identical to a cosmological constant. It was verified that the Chaplygin gas models interpolate the baryonic model and the cosmological constant model as $`\overline{A}`$ varies from zero to unity. The results for the power spectrum $`P(n)=n^{3/2}\delta _n`$ in terms of $`n/n_0`$, where $`n_0`$ is reference scale corresponding to $`100Mpc`$, are displayed in figure $`1`$ for the baryonic model, cosmological constant model and for the Chaplygin gas model with different values for $`\overline{A}`$. The spectral indice for perturbations of some hundreds of megaparsecs up to the Hubble radius is $`n_s4.7`$ for the baryonic model, $`n_s4.2`$ for the cosmological constant model, $`n_s4.5`$ for the Chaplygin gas with $`\overline{A}=0.5`$.

Another important test to verify if the Chaplygin gas model may represent dark energy is the comparison with the supernova type Ia data. In order to do so, we evaluate the luminosity distance in the Chaplygin gas cosmological model. In such a model, the luminosity distance, for a flat Universe, reads

$$D_L=(1+z)_0^z\frac{dz^{}}{\sqrt{\mathrm{\Omega }_{m0}(1+z^{})^3+\mathrm{\Omega }_{c0}\sqrt{\overline{A}+(1\overline{A})(1+z^{})^6}}},$$
(4)
$`z`$ being the redshift. From this expression the following relation between the apparent magnitude $`m`$ and the absolute magnitude $`M`$ is obtained:

$$mM=5\mathrm{log}\{(1+z)_0^z\frac{dz^{}}{\sqrt{\mathrm{\Omega }_{m0}(1+z^{})^3+\mathrm{\Omega }_{c0}\sqrt{\overline{A}+(1\overline{A})(1+z^{})^6}}}\}.$$
(5)
In order to compare with the supernova data, we compute the quantity

$$\mu _0=5\mathrm{log}\left(\frac{D_L}{Mpc}\right)+25.$$
(6)
The quality of the fitting is characterized by the parameter

$$\chi ^2=\underset{i}{}\frac{[\mu _{0,i}^o\mu _{0,i}^t]^2}{\sigma _{\mu _0,i}^2+\sigma _{mz,i}^2}.$$
(7)
In this expression, $`\mu _{0,i}^o`$ is the measured value, $`\mu _{0,i}^t`$ is the value calculated through the model described above, $`\sigma _{\mu _0,i}^2`$ is the measuremente error, $`\sigma _{mz,i}^2`$ is the dispersion in the distance modulus due to the dispersion in galaxy redshift due to perculiar velocities. This quantity we will taken as

$$\sigma _{mz}=\frac{\mathrm{log}D_L}{z}\sigma _z,$$
(8)
where, following , $`\sigma _z=200km/s`$. We evaluate, in fact, $`\chi _\nu ^2`$, the estimated errors for degree of freedom.

The best fitting was obtained for the case where $`\mathrm{\Omega }_{m0}=0`$, $`\mathrm{\Omega }_{c0}=1`$, $`\overline{A}0.847`$, $`H_062.1`$, with $`T_015.1Gy`$ for the age of Universe. In this case, $`\chi _\nu ^20.760`$. As a comparison, if we substitute the Chaplygin gas by a cosmological constant ($`\overline{A}=1`$), the best fitting is given by $`\mathrm{\Omega }_{m0}0.244`$, $`\mathrm{\Omega }_\mathrm{\Lambda }0.756`$, $`H_061.8`$, and the age of the Universe is $`T_016.5Gy`$. In this model $`\chi _\nu ^20.774`$. In figure $`2`$ we display this fitting.

These results deserve some comments. First, if we consider a cosmological constant model, the best fitting is in quite good agreement with other observational data, mainly those from $`CMB`$ and the dynamics of clusters of galaxies. However, the Chaplygin gas model indicates no need of dark matter, and the fitting quality is slightly better than the cosmological model. This is interesting since there are claims that the Chaplygin gas may unify dark matter with dark energy: at small scales, the Chaplygin gas may cluster, playing the role of dark matter, remaining a smooth component at large scale .

The results presented above indicate that the Chaplygin gas may be taken seriously as a candidate to describe dark energy and dark matter. However, more analysis are needed in order to verify to which extent it can be a better alternative with respect to the cosmological constant model combined with dark energy. A crucial test is the analysis of the spectrum of anisotropy of the cosmic microwave background radiation.


format: text

text:
# Calculating two- and three-body decays with FeynArts and FormCalc

## 1 Introduction

The analytical calculation of Feynman diagrams and their subsequent numerical evaluation represent the standard technique for calculating cross sections of scattering processes and widths of particle decays in High Energy Physics. Unfortunately, the number of diagrams to be computed grows rapidly with the number of external particle legs, internal particle loops, and exchanged particles. This makes an automatic generation and computation of Feynman diagrams a highly desirable feature, not only to reduce the time needed for the calculations, but also to eliminate possible sources of errors.

For more than ten years, a well-tested and easy-to-use Feynman diagram generator has existed in the form of the Mathematica package FeynArts . More recently, a fast and equally reliable computer algebra program for the evaluation of the Dirac and color traces inherent in the generated diagrams has been implemented within Form (FormCalc ) . Together with LoopTools  (an adapted version of the Fortran program FF  for the numerical evaluation of tensor loop integrals), the FeynArts and FormCalc program package allows for an automatic computation of tree- and loop-level $`22`$ and $`23`$ scattering processes within the Standard Model and the Minimal Supersymmetric Standard Model. Since the phase space for each particular scattering process depends on the number of in- and out-going external particles, it is best evaluated directly within Fortran with specific kinematical routines. So far, only $`22`$ and $`23`$ scattering processes have been implemented within FeynArts and FormCalc . The aim of this Paper to extend their applicability also to two- and three-body decays of massive particles and to resonant particle production. We provide four new kinematical routines and one adapted routine, which allow for the computation of the total two- and three-body decay widths and of the energy distribution of one decay product for three-body decays.

The remainder of this Paper is organized as follows: In Sec. 2 we collect the well-known kinematical formulæ for two- and three-body decays. Their implementation in the new and modified Fortran routines is described in Sec. 3. As a test of the new routines, we calculate in Sec. 4 three standard examples for scalar, fermion, and vector boson decays. A short summary is given in Sec. 5.

## 2 Kinematics for two- and three-body decays

The decay width of a heavy particle with mass $`m_1`$,

$$\mathrm{d}\mathrm{\Gamma }_i=\frac{1}{F}||^2\mathrm{d}\mathrm{\Phi }_i,$$
(1)
depends on the flux factor $`F`$, i.e. the number of decaying particles per unit volume, the squared invariant scattering amplitude $`||^2`$, and the phase space volume $`\mathrm{d}\mathrm{\Phi }_i`$ for $`i`$ final-state particles. In the rest frame of the decaying particle, $`k_1=(m_1,0,0,0)`$, the flux factor is simply

$$F=2m_1.$$
(2)
The four-momenta of two decay products with masses $`m_2,m_3`$,

$`k_2`$ $`=`$ $`(E_2,|k_2|\mathrm{sin}\theta ,\mathrm{\hspace{0.17em}0},|k_2|\mathrm{cos}\theta ),`$
$`k_3`$ $`=`$ $`(E_3,|k_3|\mathrm{sin}\theta ,\mathrm{\hspace{0.17em}0},|k_3|\mathrm{cos}\theta ),`$ (3)
are uniquely defined by their energies $`E_i^2=|k_i|^2+m_i^2`$, squared three-momenta

$$|k_2|^2=|k_3|^2=\frac{(m_1^2m_2^2+m_3^2)^2}{4m_1^2}m_3^2,$$
(4)
and by the center-of-mass scattering angle $`\theta `$. The phase space volume for two-body decays is then

$`\mathrm{d}\mathrm{\Phi }_2`$ $`=`$ $`{\displaystyle \frac{\mathrm{d}^3k_2}{(2\pi )^32E_2}}{\displaystyle \frac{\mathrm{d}^3k_3}{(2\pi )^32E_3}}(2\pi )^4\delta ^{(4)}(k_1k_2k_3)`$ (5)
$`=`$ $`{\displaystyle \frac{1}{(2\pi )^2}}{\displaystyle \frac{|k_2|}{4m_1}}\mathrm{d}\mathrm{cos}\theta \mathrm{d}\varphi ,`$
where the dependence on the azimuthal angle $`\varphi `$ integrates to $`2\pi `$. The total two-body decay width is obtained by integrating over the polar angle $`\theta `$. \[Note that in the rest frame of the decaying particle, the distribution in the polar angle $`\theta `$ depends on the arbitrary definition of a reference axis.\]

For three-particle decays, the phase space volume is

$`\mathrm{d}\mathrm{\Phi }_3`$ $`=`$ $`{\displaystyle \frac{\mathrm{d}^3k_2}{(2\pi )^32E_2}}{\displaystyle \frac{\mathrm{d}^3k_3}{(2\pi )^32E_3}}{\displaystyle \frac{\mathrm{d}^3k_4}{(2\pi )^32E_4}}(2\pi )^4\delta ^{(4)}(k_1k_2k_3k_4)`$ (6)
$`=`$ $`{\displaystyle \frac{1}{(2\pi )^4}}{\displaystyle \frac{1}{8}}\mathrm{d}E_2\mathrm{d}\varphi _2\mathrm{d}E_4\mathrm{d}\mathrm{cos}\theta _4.`$
Here, the $`\delta `$-function has been used to eliminate the integration over $`k_3`$ and $`\theta _2`$, and the trivial integration over the azimuthal angle $`\varphi _4`$ has already been performed. For three-body decays, we will not only consider total decay widths, but also distributions in the decay energy $`E_4`$.

While the standard version of FeynArts and FormCalc diverges for $`22`$ and $`23`$ cross sections, which proceed through an intermediate $`s`$-channel resonance, our modified version makes it possible to calculate also these resonance cross sections using the formula

$$\sigma _{abRX}(\sqrt{s})=\frac{4\pi \mathrm{\Gamma }(Rab)\mathrm{\Gamma }(RX)}{(sm_R^2)^2+m_R^2(\mathrm{\Gamma }_{\mathrm{tot}}^R)^2},$$
(7)
where $`s`$ is the squared center-of-mass energy of the initial state particles $`a`$ and $`b`$, $`m_R`$ is the mass of the resonance $`R`$, and the decay widths are calculated as described above.

## 3 Description of the new and modified Fortran routines

The phase space integrals described in the previous Section have been implemented in FeynArts and FormCalc with two new Fortran routines for two- and three-body decays,

* 1to2.F,
* 1to3.F,

and two corresponding common blocks,

* 1to2.h,
* 1to3.h.

The new Fortran routines have been adapted from the existing routines for $`22`$ and $`23`$ particle scattering by modifying variable and center-of-mass energy definitions, the flux factor, and the initial-state four-momenta. The Fortran routine

* num.F,

which handles the averaging and the summation over initial and final state spins, has been modified to accomodate single-particle initial states. For this purpose, the routine num.F checks if the new preprocessor variable DECAY has been defined in 1to2.F or 1to3.F or not.111The new Fortran routines can be obtained from klasen@mail.desy.de.

## 4 Three examples: $`hb\overline{b}`$, $`\mu e\overline{\nu }_e\nu _\mu `$, and $`Z\nu _e\overline{\nu }_e`$

In order to demonstrate that the new kinematical routines work properly within the extended package of FeynArts and FormCalc and are straightforward to use, we compute in this Section three standard decay widths. Note that the FormCalc routine process.h has to be adapted to every considered scattering process. In particular, the correct phase space generator has to be included at the end of process.h. However, the calling structure of the main program run.F has not been changed. For more details we refer the reader to the FeynArts  and FormCalc  manuals.

### 4.1 $`hb\overline{b}`$

The isotropic, scalar two-body decay $`hb\overline{b}`$ represents the dominant Higgs decay mode for a Standard Model Higgs boson with mass below $`m_H=140`$ GeV. Using the default parameter settings of FeynArts and FormCalc  and the Standard Model initialization file sm\_ini.F, we obtain the decay widths listed in Table 1.

Identical results are obtained with the tree-level formula

$$\mathrm{\Gamma }(Hb\overline{b})=\frac{\alpha N_Cm_Hm_b^2}{8\mathrm{sin}^2\theta _Wm_W^2}\left(1\frac{4m_b^2}{m_H^2}\right)^{3/2}$$
(8)
where $`1/\alpha =137`$ is the inverse of the electromagnetic fine structure constant, $`\mathrm{sin}^2\theta _W=1m_W^2/m_Z^2`$ is the electroweak mixing angle, $`N_C=3`$ accounts for the color degree of freedom, $`m_H`$, $`m_Z=91.1882`$ GeV, and $`m_W=80.419`$ GeV are the masses of the Higgs, $`Z`$-, and $`W`$-bosons, and $`m_b=4.7`$ GeV is the mass of the bottom quark. We have checked that the cross section differential in the scattering angle $`\theta `$ is indeed constant.

### 4.2 $`\mu e\overline{\nu }_e\nu _\mu `$

In the Standard Model, the three-body decay of muons with mass $`m_\mu =105.658`$ MeV into electrons and neutrinos proceeds through a virtual $`W`$-boson and accounts for almost 100% of the muonic decay width. With the default parameter settings of FeynArts and FormCalc , we obtain

$$\mathrm{\Gamma }(\mu e\overline{\nu }_e\nu _\mu )=2.8\times 10^{10}\mathrm{eV}$$
(9)
which agrees with the result obtained from the lowest order formula

$$\mathrm{\Gamma }(\mu e\overline{\nu }_e\nu _\mu )=\frac{G^2m_\mu ^5}{192\pi ^3},$$
(10)
where we have averaged over the initial-state spin, Fermi’s constant $`G`$ has been computed from $`G=\pi \alpha /(\sqrt{2}\mathrm{sin}\theta _W^2m_W^2)`$, and the electron mass has been neglected. Furthermore, the calculated muon life time $`\tau =\mathrm{}/\mathrm{\Gamma }=2.3\times 10^6`$ s is in good agreement with the measured value of $`2.2\times 10^6`$ s. As mentioned above, for three-body decays we consider also the energy distribution of one of the decay products, in this case of the observed electron

$$\frac{\mathrm{d}\mathrm{\Gamma }}{\mathrm{d}E_e}(\mu e\overline{\nu }_e\nu _\mu )=\frac{G^2m_\mu ^2}{12\pi ^3}E_e^2\left(3\frac{4E_e}{m_\mu }\right).$$
(11)
As can be seen in Fig. 1, the numerical evaluation of Eq.

(11) agrees very well with the result of FeynArts and FormCalc .

### 4.3 $`Z\nu _e\overline{\nu }_e`$

In our final example, we consider the invisible branching ratio of the $`Z`$-boson in the Standard Model, which represents an important boundary condition for the number of neutrinos realized in nature. The tree-level formula

$$\mathrm{\Gamma }(Z\nu _e\overline{\nu }_e)=\frac{\alpha m_Z}{24\mathrm{sin}^2\theta _W\mathrm{cos}^2\theta _W}$$
(12)
and our modified program package FeynArts and FormCalc unanimously lead to the numerical result

$$\mathrm{\Gamma }(Z\nu _e\overline{\nu }_e)=160\mathrm{MeV}.$$
(13)
For three generations of neutrinos and employing the measured total width of the $`Z`$-boson, $`\mathrm{\Gamma }_{\mathrm{tot}}=2.49`$ GeV, this translates into an invisible branching ratio of

$$\mathrm{BR}_{\mathrm{inv}}=3\frac{\mathrm{\Gamma }(Z\nu _e\overline{\nu }_e)}{\mathrm{\Gamma }_{\mathrm{tot}}}=19.3\%,$$
(14)
which is close to the experimental value of 20.0%.

## 5 Summary

In this Paper, we have described the implementation of two- and three-body decays and resonance cross sections in the computer algebra package FeynArts and FormCalc . We have provided four new Fortran routines, which generate the phase space for two- and three-particle decays, and we have modified one existing routine in order to account for correct spin averages of one-particle initial states. While the new routines have been obtained in a straight-forward way and are easy to use, they will hopefully serve a large variety of purposes and a wide community in High Energy Physics to obtain interesting new results for tree- and loop-level particle decays or resonance cross sections.

## Acknowledgments

This work has been supported by the Deutsche Forschungsgemeinschaft through Grant No. KL 1266/1-3.


text:
# Field-Dependent Critical Current in Type-II Superconducting Strips: Combined Effect of Bulk Pinning and Geometrical Edge Barrier

## Abstract

Recent theoretical and experimental research on low-bulk-pinning superconducting strips has revealed striking dome-like magnetic-field distributions due to geometrical edge barriers. The observed magnetic-flux profiles differ strongly from those in strips in which bulk pinning is dominant. In this paper we theoretically describe the current and field distributions of a superconducting strip under the combined influence of both a geometrical edge barrier and bulk pinning at the strip’s critical current $`I_c`$, where a longitudinal voltage first appears. We calculate $`I_c`$ and find its dependence upon a perpendicular applied magnetic field $`H_a`$. The behavior is governed by a parameter $`p`$, defined as the ratio of the bulk-pinning critical current $`I_p`$ to the geometrical-barrier critical current $`I_{s0}`$. We find that when $`p>2/\pi `$ and $`I_p`$ is field-independent, $`I_c`$ vs $`H_a`$ exhibits a plateau for small $`H_a`$, followed by the dependence $`I_cI_pH_a^1`$ in higher magnetic fields.

The combination of a geometrical edge barrier and bulk pinning recently has been shown to strongly affect the properties of low-dimensional superconductors (thin films, single crystals, and tapes with high demagnetizing factors) placed in either a perpendicular magnetic field  or a transport-current-carrying state . While most experimental studies of the field dependence of the critical current $`I_c`$ are being interpreted solely on the basis of bulk-pinning theory (see for example ), a number of works  have shown that a geometrical edge barrier (or surface barrier) may strongly affect $`I_c`$. In this paper we study the combined effect of a geometrical edge barrier and bulk pinning upon the magnetic field dependence of $`I_c`$ for type-II superconducting strips. We shall show how the dependence of $`I_c`$ upon $`H_a`$ is controlled by the parameter $`p=I_p/I_{s0}`$, where $`I_p`$ is the bulk-pinning critical current in the absence of a geometrical edge barrier, and $`I_{s0}`$ is the geometrical-barrier critical current in the absence of bulk pinning.

We consider a superconducting strip of thickness $`d`$ ($`|y|<d/2`$) and width $`2W`$ ($`|x|<W`$) centered on the $`z`$ axis. We assume that $`d`$ is less than the London penetration depth $`\lambda `$ and that $`W`$ is much larger than the two-dimensional screening length $`\mathrm{\Lambda }=2\lambda ^2/d`$. The strip is subjected to a perpendicular applied magnetic field $`𝐇_a=(0,H_a,0)`$, and it carries a total current $`I`$ in the $`z`$ direction described by a spatially dependent sheet current density $`𝐊(x)=𝐉d=[0,0,K_z(x)]`$. We wish to determine the current-density and magnetic-field distributions at the critical current at which a steady-state flux-flow voltage appears along the length of the strip. For a strip containing no magnetic flux, $`K_z(x)`$ is the sum of two contributions,

$$K_{Iz}(x)=\frac{I}{\pi \sqrt{W^2x^2}},$$
(1)
the Meissner-state current density generated by the applied current $`I`$, and

$$K_{az}(x)=\frac{2H_ax}{\sqrt{W^2x^2}},$$
(2)
the Meissner-state current density induced by the applied field $`H_a`$. The divergences in Eqs. (1) and (2) at $`|x|=W`$ are cut off when $`x`$ is within $`\mathrm{\Lambda }`$ of the edge.

To account for the edge barrier, we assume that vortices nucleate and enter the superconductor when $`K_z`$ at either sample edge reaches the value $`K_s=j_sd`$ at which the barrier is overcome. For an ideal edge, $`j_s`$ is equal to the Ginzburg-Landau depairing current density $`j_{GL}`$ , but for an extremely defected edge, $`j_s`$ may become negligibly small. When $`H_a=0`$, the sheet current at both edges is approximately $`K_{Iz}I/\pi \sqrt{2W\mathrm{\Lambda }}`$, such that the edge-barrier critical current in zero external magnetic field is $`I_{s0}\pi K_s\sqrt{2W\mathrm{\Lambda }}`$. When $`H_a>0`$, the net sheet current at $`x=W`$ is approximately $`K_z(I+2\pi H_aW)/(\pi \sqrt{2W\mathrm{\Lambda }})`$, and the edge-barrier critical current becomes $`I_s(H_a)/I_{s0}=1h`$ for small $`h`$, where $`h=H_a/(I_{s0}/2\pi W).`$ This result is essentially the same as that found in Ref. for the critical current in low applied magnetic fields for bulk-pinning-free strips.

We next account for bulk pinning, characterized via a bulk-pinning critical sheet current density, $`K_p=j_pd`$, such that the critical current in the absence of an edge barrier is $`I_p=2K_pW`$. We first consider the case of relatively weak bulk pinning when $`I_p<(2/\pi )I_{s0}`$, i.e., $`p<2/\pi .`$ In low fields ($`0<H_a<H_d`$, region I of Fig. 1), vortices nucleate on the right-hand side at $`x=W`$ when $`I`$ slightly exceeds $`I_s(H)`$. As long as $`K_z(x)=K_{Iz}(x)+K_{az}(x)`$ exceeds $`K_p`$, these vortices are driven entirely across the strip, traveling rapidly (speed $`v`$ governed solely by the force-balance equation $`[K_z(x)K_p]\varphi _0=\eta vd`$, where $`\eta `$ is the viscous drag coefficient), and annihilating with their images on the opposite side of the strip. The critical current is then $`I_c(H_a,p)=I_s(H_a)`$, and the normalized critical current is

$$i_c(h,p)=I_c(H_a,p)/I_{s0}=1h.$$
(3)
However, $`K_{min}`$, the minimum value of $`K_z(x)`$, decreases with increasing $`H_a`$ and reaches $`K_p`$ at $`I=I_s(H_a)`$ when $`H_a=H_d`$, where

$$h_d=H_d/(I_{s0}/2\pi W)=\frac{1}{2}[1(\pi p/2)^2].$$
(4)
When $`p<2/\pi `$ and $`h>h_d`$ or when $`p>2/\pi `$, i.e., for $`h`$ and $`p`$ outside region I of Fig. 1, nucleated vortices stop inside the strip, and dome-like flux distributions occur at the critical current. To determine the critical current with domes present, one must first calculate the vortex (and antivortex) density $`n(x)=\mu _0H_y(x,0)/\varphi _0`$ inside the dome, where $`K_z(x)=K_p`$, and the sheet current density $`K_z(x)`$ outside, where $`n(x)=0`$. We have obtained these mathematically by using the Cauchy integral inversion method to invert the Biot-Savart law. However, we shall use the method of complex fields to give a physical interpretation of the mathematical results.

Following Ref., we express the two-dimensional field distribution as an analytic function $`H(\zeta )=H_y(x,y)+iH_x(x,y)`$ of the complex variable $`\zeta =x+iy`$, such that the Biot-Savart law becomes

$$H(\zeta )=H_a+\frac{1}{2\pi }_W^W\frac{K_z(u)}{\zeta u}𝑑u.$$
(5)
The inversion procedure yields $`H(\zeta )`$ of the following form at the critical current $`I_c(H_a)`$:

$$H(\zeta )=\frac{(\zeta a)^{1/2}(\zeta b)^{1/2}}{(\zeta ^2W^2)^{1/2}}[H_a+\frac{K_p}{2\pi }Q(a,b,\zeta )],$$
(6)
where

$`Q(a,b,\zeta )`$ $`={\displaystyle _a^b}{\displaystyle \frac{\sqrt{W^2u^2}}{(\zeta u)\sqrt{(ua)(bu)}}}𝑑u`$ (8)
$`={\displaystyle \frac{2(W+a)}{\sqrt{(Wa)(W+b)}}}[{\displaystyle \frac{(W\zeta )}{(\zeta a)}}𝚷({\displaystyle \frac{(\zeta +W)(ba)}{(\zeta a)(W+b)}},q)`$ (9)
$`+𝚷({\displaystyle \frac{ba}{W+b}},q)],`$ (10)
and

$$q^2=\frac{2W(ba)}{(Wa)(W+b)}.$$
(11)
In Eq. (6), the term proportional to $`H_a`$ is simply the complex field describing the Meissner-state response to the applied field $`H_a`$ of two parallel superconducting strips ($`W<x<a`$ and $`b<x<W`$). The term proportional to $`K_p`$ is the complex field describing the image-current response of the two strips to currents $`K_pdu`$ summed over the region $`a<u<b`$. We have evaluated the integral in Eq. (8) in terms of complete elliptic integrals of the third kind $`𝚷(n,k)`$, where $`n`$ is called either the characteristic or the parameter, and $`k`$ is called the modulus. Equation (10) can be used to evaluate $`H_y(x,0)=\mathrm{Re}H(x)`$ and $`K_z(x)=2\mathrm{I}\mathrm{m}H(x+iϵ)`$.

For $`p>2/\pi `$ and small values of $`h`$, i.e., for $`h`$ and $`p`$ in region IV of Fig. 1, the vortex distribution at the critical current can be described as a double dome, consisting of a vortex dome adjacent to an antivortex dome (see inset). Just above the critical current, vortices nucleate at $`x=W(x^{}=1)`$, where $`K_z(W\mathrm{\Lambda })=K_s`$, move rapidly to the left through an otherwise vortex-free region ($`b<x<W`$), and then move slowly to the left through a vortex-filled region (the vortex dome), $`a_+<x<b`$. Antivortices nucleate at $`x=W`$, where $`K_z(W+\mathrm{\Lambda })=K_s`$, move rapidly to the right through an otherwise vortex-free region ($`W<x<a`$), and then move slowly to the right through an antivortex-filled region (the antivortex dome), $`a<x<b_{}`$. Vortices and antivortices annihilate where the two domes meet at $`x=b_{}=a_+`$.

Two equations must be solved simultaneously for $`a`$ and $`b`$ at the critical current $`I_c(H_a)`$ for known values of $`h`$ and $`p`$ in region IV of Fig. 1. One condition is that $`K_z(W\mathrm{\Lambda })=K_s`$, which yields from Eqs. (6) and (10)

$`\sqrt{(1a^{})(1b^{})}h+(1+a^{})\sqrt{{\displaystyle \frac{1b^{}}{1+b^{}}}}𝚷({\displaystyle \frac{b^{}a^{}}{1+b^{}}},q)p`$ (12)
$`=1,`$ (13)
where we use the normalized quantities $`a^{}=a/W`$ and $`b^{}=b/W`$. The other condition, that $`K_z(W+\mathrm{\Lambda })=K_s`$, yields

$`\sqrt{(1+a^{})(1+b^{})}h+(1b^{})\sqrt{{\displaystyle \frac{1+a^{}}{1a^{}}}}𝚷({\displaystyle \frac{b^{}a^{}}{1a^{}}},q)p`$ (14)
$`=1.`$ (15)
Expansion of Eq. (5) for large $`\zeta `$ yields $`H(\zeta )=H_a+I/2\pi \zeta +\mathrm{O}(1/\zeta ^2)`$. Expanding Eqs. (6) and (8) in powers of $`\zeta `$ and making use of Eq. (13), we obtain the normalized critical current, $`i_c=I_c/I_{s0}`$:

$`i_c(h,p)={\displaystyle \frac{a^{}+b^{}}{2\sqrt{(1a^{})(1b^{})}}}`$ (16)
$`+{\displaystyle \frac{1}{2}}\sqrt{{\displaystyle \frac{1+b^{}}{1a^{}}}}[(1a^{})𝐄(q)+(1+a^{})𝐊(q)]p,`$ (17)
where $`a^{}`$ and $`b^{}`$ are determined from Eqs. (13) and (15) for the desired values of $`h`$ and $`p`$.

The double-dome vortex-antivortex distribution (region IV of Fig. 1) occurs at $`i_c`$ only for $`h`$ in the range $`0<h<h_1`$ for $`p>2/\pi `$. Here, $`h_1(p)`$ is the lowest value of $`h`$ that makes $`H_y(x,0)>0`$ in the region $`a<x<b`$. Thus, one of the equations determining $`h_1`$ is $`H_a+(K_p/2\pi )Q(a,b,a+ϵ)=0`$ \[see Eq. (6)\], which yields

$`h+{\displaystyle \frac{1}{(b^{}a^{})\sqrt{(1a^{})(1+b^{})}}}`$ (18)
$`[(1a^{})(1+b^{})𝐄(q)(1+a^{})(1b^{})𝐊(q)`$ (19)
$`(b^{}a^{})(1b^{})𝚷({\displaystyle \frac{b^{}a^{}}{1a^{}}},q)]p=0.`$ (20)
$`h_1`$ can be determined for a given value of $`p`$ as the value of $`h`$ when Eqs. (13), (15), and (20) are simultaneously solved for $`a^{},b^{},`$ and $`h`$.

For known values of $`h`$ and $`p`$ in region III of Fig. 1, the left and right boundaries of the vortex dome $`a^{}`$ and $`b^{}`$ are determined by simultaneously solving Eqs. (13) and (20); Eq. (20) also gives the condition that $`dK_z(x)/dx=0`$ at $`x=a`$. Once $`a^{}`$ and $`b^{}`$ are found, Eq. (17) again can be used to calculate the critical current. Just above the critical current, vortices nucleate at $`x=W`$, move rapidly to the left through the vortex-free region $`b<x<W`$, move slowly to the left through the vortex dome $`a<x<b`$, escape from the dome, and finally move rapidly to the left through the vortex-free region $`W<x<a`$.

For increasing values of $`h`$, the left boundary of the vortex-filled region moves closer to the left edge of the strip; $`a`$ becomes equal to $`W+\mathrm{\Lambda }`$ when $`h=h_2`$, which can be determined for a given value of $`p`$ as the value of $`h`$ when Eqs. (13) and (20) are numerically solved for $`b^{}`$ and $`h`$, taking $`a^{}=1+\mathrm{\Lambda }/W`$. (For Figs. 1 and 2, $`\mathrm{\Lambda }/W=0.01`$ was assumed.) For $`h>h_2`$, the field and current distributions and $`i_c`$ can be calculated with good accuracy by simply setting $`a=W`$ in Eq. (6). The complex field in region II of Fig. 1 is then

$$H_{II}(\zeta )=\frac{(\zeta b)^{1/2}}{(\zeta W)^{1/2}}[H_a+\frac{K_p}{2\pi }Q(b,\zeta )],$$
(21)
where

$`Q(b,\zeta )={\displaystyle _a^b}{\displaystyle \frac{\sqrt{Wu}}{(\zeta u)\sqrt{bu}}}𝑑u`$ (23)
$`=2\mathrm{sinh}^1\sqrt{{\displaystyle \frac{W+b}{Wb}}}2\mathrm{sinh}^1\sqrt{{\displaystyle \frac{(W+b)(\zeta W)}{(Wb)(\zeta +W)}}}.`$ (24)
The condition $`K_z(W\mathrm{\Lambda })=K_s`$, which determines $`b^{}=b/W`$ at $`I_c`$ for $`h`$ and $`p`$ in region II of Fig. 1, becomes

$$\sqrt{2(1b^{})}[h+p\mathrm{sinh}^1\sqrt{(1+b^{})/(1b^{})}]=1$$
(25)
instead of Eq. (13), and the normalized critical current can be expressed as

$$i_c(h,p)=\frac{1}{4}\sqrt{2(1b^{})}+\frac{p}{2}\sqrt{2(1+b^{})}$$
(26)
instead of Eq. (17).

The reduced critical current $`i_c`$ as a function of $`h`$, calculated from Eq. (17) or (26), is shown in Fig. 2 for various values of $`p`$. For $`p=0`$ and $`h>1/2`$, we find $`i_c=1/4h`$, as obtained in  for pin-free strips with an edge barrier. In the opposite limit, $`p1`$, we obtain $`i_cp`$, as expected for bulk-pinning-dominated behavior. For $`hp`$, we see from Eq. (25) that $`b^{}`$ approaches 1, and Eq. (26) yields $`i_cp+1/4h.`$ A generic feature of Fig. 2 is the plateau in $`i_c`$ vs $`h`$ to the left of the dashed curve in region IV; actually, $`i_c`$ increases by a few percent as $`h`$ increases from 0 to $`h_1`$. This increase is due to a significant decrease in the width $`ba`$ over which $`K_z(x)`$ is restricted to $`K_p`$. This effect is partially compensated by a change in shape of $`K_z(x)`$ in the vortex-free zones \[e.g., $`dK_z(x)/dx=0`$ at $`x=a`$ at $`h=h_1`$\]. Field-dependent critical current densities $`j_c(H_a)`$ have been found experimentally in Refs., but the behavior was interpreted solely in terms of bulk pinning.

In agreement with earlier work, our results show that the critical current of a strip is not a simple superposition of currents $`I_s`$ and $`I_p`$, as was suggested in Refs. . Only in the limit $`hp`$ is it possible to express the critical current as $`I_c(H_a)=I_p+I_s(H_a)`$.

We have assumed here that $`K_p`$ is a constant. Because of the nonlocal current-field relation \[Eq. (5)\], it would be a challenging task to find $`I_c(H_a)`$ when $`K_p`$ depends upon the local magnetic field $`H_y(x,0)`$.

In summary, we have solved for the field dependence of the critical current density in a superconducting strip accounting for both bulk pinning and a geometrical edge barrier, and we have developed a procedure for finding the magnetic-field and current-density distributions inside the strip at the critical current. In the presence of a strong edge barrier, we have found strong field dependencies of the critical current. Such effects should be taken into account when interpreting experimental critical currents in low and moderate magnetic fields.

We thank Y. Mawatari for stimulating discussions. This work was funded by the Basic Research Foundation of Russia through Grant No. 01-02-16593, the Education Ministry (Grant No. E00-3.4-331) and the Science Ministry (Project 107-1(00)-P) of the Russian Federation, by Iowa State University of Science and Technology under Contract No. W-7405-ENG-82 with the U.S. Department of Energy, and by the National Academy of Sciences under the Collaboration in Basic Science and Engineering Program supported by Contract No. INT-0002341 from the National Science Foundation.


text:
# Renormalization and destruction of 1/𝛾² tori in the standard nontwist map

## I Introduction

In this paper we consider the standard nontwist map (SNM) $`M`$, as introduced in Ref. diego2, :

$`x_{n+1}`$ $`=`$ $`x_n+a\left(1y_{n+1}^2\right)`$
$`y_{n+1}`$ $`=`$ $`y_nb\mathrm{sin}\left(2\pi x_n\right),`$ (1)
where $`(x,y)𝕋\times `$, $`a(0,1)`$, and $`b(\mathrm{},\mathrm{})`$. The map $`M`$ is area-preserving and violates the twist condition

$$\frac{x_{i+1}(x_i,y_i)}{y_i}0(x_i,y_i),$$
(2)
along a curve in phase space, which has been recently called the nonmonotone curve.petrisor1  Traditionally, most studies of area-preserving maps have dealt with the twist case, but in recent years more and more research has been focused on the nontwist case.

Applications of nontwist maps occur in many fields, for example, the study of magnetic field lines in toroidal plasma devices (see e.g. Refs. horton1, ; balescu, ), in celestial mechanics,kyner  fluid dynamicsdiego2  and atomic physics.chandre2  It has been showndullin ; vander  that nontwist regions appear generically in area-preserving maps that have a tripling bifurcation of an elliptic fixed point. In addition to these applications, the map is quite interesting from a mathematical standpoint because many important theorems in the theory of area-preserving maps assume the validity of the twist condition, e.g. the KAM theorem and the Poincare-Birkhoff theorem. The SNM can serve as a model for the development of new proofs. Up to now, only a few mathematical results exist for nontwist maps (see e.g. Ref. delshams, ; franks, ; petrisor1, ; simo, ).

We continue the work of del-Castillo-Negrete, Greene and Morrison,diego4 ; diego5  who studied the breakup of the shearless invariant torus with winding number $`1/\gamma `$, where $`\gamma =(1+\sqrt{5})/2`$ is the golden mean. We present the analysis of the breakup of the shearless invariant torus with winding number (in continued fraction representation)

$$\omega =[0,2,1,1,\mathrm{}]=1/\gamma ^2.$$
(3)
Because this winding number is a noble number (its continued fraction expansion ends with $`[1,1,1,\mathrm{}]`$), the behavior of the residues of the approximating periodic orbits is expected to be the same as in the $`1/\gamma `$ case, i.e. we should find the same fixed point of the renormalization group operator with the same unstable eigenvalues that were found in Ref. diego5, . But, the form of the renormalization group operator, which is defined later in Sec. IV, is different from the $`1/\gamma `$ case. Also, the region of parameter space we study is different. Additionally, since the periods of approximating periodic orbits are bigger than those for the $`1/\gamma `$ case, the present work serves as a test case for improved numerical techniques described later in Sec. III.1.

A different approach, which yields rough parameter values for the breakup of invariant tori, was used by Shinohara and Aizawa in Ref. shin1, , who showed that a shearless invariant torus crosses the $`x`$-axis at two points,note1  $`x_A=a/21/4`$ and $`x_B=a/2+1/4`$.

For a given $`(a,b)`$, a point on the shearless torus, $`(a/2+1/4,0)`$, is iterated many times (we used $`10^6`$). If the $`y`$ value stays below a threshold (we used $`|y|<0.52`$), it is assumed that the shearless curve exists and the point is plotted. Figure 1 depicts our duplication of their procedure.

We see that the critical point for the $`1/\gamma ^2`$ shearless curve (indicated by $``$) lies on the boundary in Fig. 1. Thus, the boundary points of Fig. 1 represent the critical function for the SNM. This is a generalization of the definition of the critical function for the standard twist map (see e.g. Ref. stark, ), which has only one parameter e.g. $`k`$. The critical function in the twist case is then defined as $`k_c(\omega )`$. Here, we have two parameters, but the shearless invariant torus of a given winding number $`\omega `$ exists only for parameter values belonging to a curve $`(a,b(a;\omega ))`$ in the parameter space. Thus, we can define the critical function by the critical points on each of those curves by $`(a_c,b(a_c;\omega ))`$. By finding the critical points for many other winding numbers (both nobles and non-nobles), we hope to find a more accurate critical function curve than the one shown in Fig. 1.

In Sec. II, we review some basic properties of the SNM. The detailed breakup of the shearless invariant torus with winding number $`1/\gamma ^2`$ is presented in Sec. III, which also contains a discussion of the numerical procedures involved. In Sec. IV, we interpret the results within the framework of the renormalization group. Section V a summary and some directions of future research.

## II Review of residue criterion and standard nontwist map

In this section, we give a brief review of some basic concepts of the theory of area-preserving maps in the context of the SNM. For a more in-depth discussion the reader is referred to Ref. diego4,  and references therein.

### II.1 Periodic orbits and residue criterion

Since the pioneering work of Greene,greene2  periodic orbits have proven to be very useful for studying the breakup of invariant tori in area-preserving maps. Below are some standard definitions.

An orbit of an area-preserving map $`M`$ is a sequence of points $`\left\{(x_i,y_i)\right\}`$ such that $`M(x_i,y_i)=(x_{i+1},y_{i+1})`$. The winding number $`\omega `$ of an orbit is defined as $`\omega =lim_i\mathrm{}(x_i/i)`$ when the limit exists. Here the $`x`$-coordinate is “lifted” from $`𝕋`$ to $``$. A periodic orbit of period $`n`$ is a sequence of $`n`$ points $`\left\{(x_i,y_i)\right\}_{i=1}^n`$, such that $`M^n(x_i,y_i)=(x_i+m,y_i)`$ for all $`i=1,\mathrm{},n`$, and $`m`$ is an integer. Periodic orbits have rational winding numbers $`\omega =m/n`$. An invariant torus is an orbit with irrational winding number that covers densely a one-dimensional set in phase space. Of particular importance are the invariant tori that wind around the $`x`$-domain because, in two-dimensional maps, they act as transport barriers.

The linear stability of a periodic orbit is determined by the value of its residue,greene2  $`R`$, which is defined as $`R:=\left[2\text{Tr}(L)\right]/4`$. Here, $`L`$ is the map $`M^n`$ linearized about the periodic orbit of interest and Tr denotes the trace. If $`0<R<1`$, the orbit is stable or elliptic; if $`R<0`$ or $`R>1`$, it is unstable or hyperbolic; in the degenerate cases $`R=0`$ and $`R=1`$, it is parabolic.

Periodic orbits can be used to systematically approximate invariant tori.greene2  The method is based on the observation that given a sequence of rational numbers $`\left\{m_i/n_i\right\}`$ whose limit is $`\omega `$, the sequence of periodic orbits with winding numbers $`\left\{m_i/n_i\right\}`$ approaches the invariant torus with winding number $`\omega `$ in phase space. It is important to find the “best” possible sequence, i.e. the sequence that converges to $`\omega `$ the fastest. The elements of the best possible sequence (see e.g. Ref. khinchin, ) are the convergents that are obtained from successive truncations of the continued fraction expansion of $`\omega `$.

The residue criteriongreene2  can be stated as follows: Consider an invariant torus with winding number $`\omega `$. Let $`\{m_i/n_i\}`$ be the sequence of convergents approximating $`\omega `$, and $`R_i`$ the residues of their corresponding periodic orbits.

1. If $`lim_i\mathrm{}|R_i|=0`$, the invariant torus exists.
2. If $`lim_i\mathrm{}|R_i|=\mathrm{}`$, the invariant torus is destroyed.
3. At the boundary in parameter space between those two limits, the invariant torus is at the threshold of destruction and the residues either converge to a constant, non-zero value, or there are convergent subsequences.

This criterion is based on the idea that the destruction of an invariant torus is caused by the de-stabilization of nearby periodic orbits. The residue criterion has been used successfully in many cases to predict with high precision the threshold for the destruction of invariant tori. Several theorems have been proved that lend mathematical support to the criterion.falc ; mackay2

The numerical search for periodic orbits is difficult because, in principle, it is a two-dimensional root finding problem. However, the task is considerably simplified for reversible maps,greene2 ; devogel  which are maps that can be factored as $`M=I_1I_0`$, where $`I_{0,1}`$ are involution maps that satisfy $`I_1^2=I_0^2=1`$. The sets of fixed points of the involution maps, $`\mathrm{\Gamma }_{0,1}=\left\{(x,y)\right|I_{0,1}(x,y)=(x,y)\}`$, are one-dimensional sets, called symmetry lines of the map. Once we know $`\mathrm{\Gamma }_{0,1}`$, the search for periodic orbits is reduced to a one-dimensional root finding problem, as explained below in Sec. III.1.1.

### II.2 Standard nontwist map

The SNM is reversible. The symmetry lines $`\mathrm{\Gamma }_0`$, composed of fixed points of $`I_0`$ are $`s_1=\left\{(x,y)\right|x=0\}`$ and $`s_2=\left\{(x,y)\right|x=1/2\}`$. The symmetry lines $`\mathrm{\Gamma }_1`$, composed of fixed points of $`I_1`$ are $`s_3=\left\{(x,y)\right|x=a\left(1y^2\right)/2\}`$ and $`s_4=\left\{(x,y)\right|x=a\left(1y^2\right)/2+1/2\}`$.

A major difference between the standard nontwist map and twist maps is that there are two periodic orbits, if they exist, with the same winding number on each symmetry line. This can be seen easily in the integrable case. For $`b=0`$, the $`m/n`$ periodic orbits on the $`s_1`$ symmetry line are located at

$$(x,y)=(0,\pm \sqrt{1(m/n)/a}).$$
(4)
We will call the orbit with the bigger (smaller) $`y`$-coordinate the up (down) periodic orbit.

The SNM is also invariant with respect to the transformation

$$T(x,y)=(x+1/2,y).$$
(5)
The coordinates of the up and down periodic orbits on the symmetry lines $`s_i`$, denoted by $`(x_{ui},y_{ui})`$ and $`(x_{di},y_{di})`$ respectively, are related by this symmetry as follows:

$`(x_{d2},y_{d2})=T\left((x_{u1},y_{u1})\right),(x_{u2},y_{u2})=T\left((x_{d1},y_{d1})\right),`$
$`(x_{d4},y_{d4})=T\left((x_{u3},y_{u3})\right),(x_{u4},y_{u4})=T\left((x_{d3},y_{d3})\right).`$ (6)
Therefore, it is actually enough to compute periodic orbits on $`s_1`$ and $`s_3`$, since the orbits along the other symmetry lines can be obtained from (6).

### II.3 Periodic orbit collisions and bifurcation curves

Periodic orbits in the SNM can undergo a particular form of bifurcation that occurs when the up and down periodic orbits of the same winding number meet on the symmetry line. These collisions were detected numerically in Refs. stix, ; howard,  and diego4, . Further studies of this bifurcation can be found in Refs. petrisor1, ; simo, .

From (4) it follows that, for a given $`a`$, only periodic orbits with $`m/n<a`$ exist at $`b=0`$. As the value of $`b`$ increases, the up and down orbits approach each other and at the bifurcation value, they collide and annihilate each other. For higher values of $`b`$, both orbits no longer exist. Figure 2 illustrates the behavior of periodic orbits as we increase $`b`$ from $`b=0`$. Here the $`y`$-coordinates of the $`m/n=3/8`$ periodic orbits on $`s_1`$ is shown as a function of $`b`$ for the fixed value of $`a=0.4`$.

Based on these numerical observations, the notion of a bifurcation curve in parameter space was defined in Ref. diego4, . The $`m/n`$-bifurcation curve $`b=\mathrm{\Phi }_{m/n}(a)`$ is the set of $`(a,b)`$ values for which the $`m/n`$ up and down periodic orbits are at the point of collision. The main property of this curve is that for $`(a,b)`$ values below $`b=\mathrm{\Phi }_{m/n}(a)`$, the $`r/s`$ periodic orbits, with $`r/s<m/n`$ exist. Thus, $`m/n`$ is the maximum winding number for parameter values along the $`m/n`$-bifurcation curve.

The idea of approximating invariant tori with irrational winding numbers by periodic orbits is used to define the bifurcation curve for an invariant torus as follows:diego4  The $`\omega `$-bifurcation curve $`b=\mathrm{\Phi }_\omega (a)`$ for an irrational $`\omega `$ is the set of $`(a,b)`$ values such that $`b=\mathrm{\Phi }_\omega (a)=lim_i\mathrm{}\mathrm{\Phi }_{m_i/n_i}(a)`$, where $`\mathrm{\Phi }_{m_i/n_i}(a)`$ is the $`m_i/n_i`$-bifurcation curve and $`\{m_i/n_i\}`$ are the convergents of $`\omega `$. For $`(a,b)`$ points along the $`\omega `$-bifurcation curve the invariant torus with irrational winding number $`\omega `$ is the curve of maximum winding number and is called shearless. Figure 3 depicts the bifurcation curves for several convergents of $`1/\gamma ^2`$. This figure also makes it plausible that the limit in the above definition exists.

## III Breakup of Torus with $`\omega =1/\gamma ^2`$

In this section, we present the analysis of the breakup of the shearless invariant torus in the standard nontwist map with winding number $`\omega =1/\gamma ^2`$. Tables 1 and 2 list the convergents used for these calculations. For more details see Ref. wurm, .

### III.1 Numerical methods

The computational steps necessary to find the critical point and the residue behavior of the approximating periodic orbits are as follows:

1. Find a good approximation to the $`1/\gamma ^2`$-bifurcation curve in $`(a,b)`$-space using the bifurcation curves for its convergents.
2. Along this bifurcation curve, find the up and down periodic orbits on the symmetry line $`s_1`$ that approximate the invariant torus, and compute their residues.
3. Locate the $`(a,b)`$ point along the curve at which the residues exhibit critical behavior.
4. Find the residues of the periodic orbits at criticality along the remaining symmetry lines.
5. Find the eigenvalues of the unstable eigenmodes of the renormalization group operator. The details of how to do this depend crucially on the type of critical scaling behavior that is exhibited by the residues.

#### III.1.1 Searching for periodic orbits

Periodic orbits on the symmetry lines can be be computed relatively easily for reversible maps using the following result:diego4  If $`(x,y)\mathrm{\Gamma }_{0,1}`$ then $`M^n(x,y)=(x,y)`$ if and only if $`M^{n/2}(x,y)\mathrm{\Gamma }_{0,1}`$ (for $`n`$ even) or $`M^{(n\pm 1)/2}(x,y)\mathrm{\Gamma }_{1,0}`$ (for $`n`$ odd). Thus, for example, periodic orbits with odd period $`n`$ on the $`s_1`$ symmetry line can be obtained by looking for points on $`s_1`$ that are mapped to $`s_3`$ or $`s_4`$ after $`(n+1)/2`$ iterations. This can be implemented as a one-dimensional root finding problem by considering the zeros of the function $`F(y)=\mathrm{sin}\left[2\pi \left(\widehat{x}a\left(1\widehat{y}^2\right)/2\right)\right]`$, where $`(\widehat{x},\widehat{y}):=M^{(n+1)/2}(0,y)`$. The sine function is used to eliminate the difference between $`s_3`$ and $`s_4`$. Similar ideas can be applied to find other orbits.

#### III.1.2 Finding $`m/n`$-bifurcation curves

Recall that the bifurcation curve $`\mathrm{\Phi }_{m/n}(a)`$ of a periodic orbit of winding number $`m/n`$ was defined in Sec. II.3 to be the set of points $`(a,b)`$, at which the up and down periodic orbits of winding number $`m/n`$ collide along the $`s_1`$ symmetry line. Thus, at a given value of $`a`$, the function $`F(y)`$ has two roots for $`b<\mathrm{\Phi }_{m/n}(a)`$, no roots (locally) for $`b>\mathrm{\Phi }_{m/n}(a)`$ and a single root, which is also an extremum, for $`b=\mathrm{\Phi }_{m/n}(a)`$. We thus search for the zero of the extremum of $`F(y)`$ as $`b`$ is varied.

To find the whole (or large portions) of a bifurcation curve, we use the monotonic nature of the curve (see Fig. 3) as follows: Given a point $`(a_1,b_1)`$ on the bifurcation curve i.e. $`b_1=\mathrm{\Phi }_{m/n}(a_1)`$, we increase $`a`$ by a fixed amount to $`a_2=a_1+a_{\text{step}}`$. We then start at the point $`(a_2,b_1)`$ and increase $`b`$ until we reach $`b_2=\mathrm{\Phi }_{m/n}(a_2)`$. To make sure that we are finding the correct bifurcation curve, we start searching $`(a,b)`$-space at $`(a,b)=(m/n,0)`$. Even then, the steps in $`a`$ cannot be taken to be too large. Experience has shown that steps in $`a`$ of $`1\times 10^5`$ or $`1\times 10^6`$ are safe. This method is unfortunately very slow because the part of the curve at small $`b`$ values is very steep and the interesting (near critical) part of the curve is far away from the $`b=0`$ limit.

We managed to drastically improve the speed of these calculations by using the following ideas:

1. Numerical evidence strongly suggests that a bifurcation curve is smooth and monotonically increasing, although it is not proved mathematically.note2  So we use linear extrapolation from two previous points to find the new value of $`b`$ around which to search for the bifurcation point. It was found that any higher order extrapolation did not improve the algorithm further.
2. To find bifurcation curves for periodic orbits with very large periods (e.g. of the order of several million) the following procedure is used: Starting at the bifurcation curve of a smaller period, we increase $`b`$ until the bifurcation curve of the higher period is reached. The advantage of this procedure is that we do not need to do the extremely time consuming calculations of the bifurcation curves for very high period orbits starting at $`b=0`$, but rather we can search for them near the region of interest.

#### III.1.3 Finding $`1/\gamma ^2`$-bifurcation curve

Recall that the $`1/\gamma ^2`$-bifurcation curve was defined as the limit of $`m_i/n_i`$ bifurcation curves, where $`m_i/n_i`$ are convergents of $`1/\gamma ^2`$. It was numerically observed that close to criticality, this limit is approached in accordance with the following scaling relation:diego4

$$\mathrm{\Phi }_{[n+1]}(a)=\mathrm{\Phi }_{1/\gamma ^2}(a)+B_n(a)\nu _1^{n/12},$$
(7)
where the $`\mathrm{\Phi }_{[n]}(a)`$ denotes the bifurcation curve of the periodic orbit with winding number $`[n]=F_n/F_{n+2}`$, $`\nu _1`$ is a number to be determined later, and $`B_n(a)`$ is a period-twelve function, i.e. $`B_{n+12}(a)=B_n(a)`$.

If Eq. (7) holds, it follows that for fixed $`a`$

$$\mathrm{\Phi }_{1/\gamma ^2}=\underset{n\mathrm{}}{lim}\frac{\mathrm{\Phi }_{[n+1]}\mathrm{\Phi }_{[n+12]}\mathrm{\Phi }_{[n]}\mathrm{\Phi }_{[n+13]}}{\left(\mathrm{\Phi }_{[n+1]}\mathrm{\Phi }_{[n]}\right)\left(\mathrm{\Phi }_{[n+13]}\mathrm{\Phi }_{[n+12]}\right)}.$$
(8)
We obtained the $`1/\gamma ^2`$-bifurcation curve using $`n=19`$ in Eq. (8), i.e. using the bifurcation curves for $`[32]`$, $`[31]`$, $`[20]`$ and $`[19]`$ (see Tables 1 and 2).

Now one can justify a posteriori the use of Eq. (7). Solving Eq. (7) with $`a=a_c`$ for $`\nu _1`$ yields:

$$\nu _1=\underset{n\mathrm{}}{lim}\left(\frac{\mathrm{\Phi }_{[n+13]}\left(a_c\right)b_c}{\mathrm{\Phi }_{[n+1]}\left(a_c\right)b_c}\right),$$
(9)
and

$$B_n(a_c)=\left(\mathrm{\Phi }_{[n+1]}\left(a_c\right)b_c\right)\nu _1^{n/12},$$
(10)
where $`(a_c,b_c)`$ is the critical point for breakup of the shearless $`1/\gamma ^2`$ invariant torus i.e. $`b_c=\mathrm{\Phi }_{1/\gamma ^2}(a_c)`$. We found that $`\nu _1^{1/12}=2.678`$. Some numerical evidence for the periodicity of $`B_n(a_c)`$ is given in Table 3.

### III.2 Results

In this subsection, we present the results of our computations.

#### III.2.1 Residue behavior at criticality

We computed bifurcation curves up to $`[32]=3524578/9227465`$ and found the critical points along them, i.e. the parameter values along those curves for which the residues of approximating periodic orbits neither converge to zero nor diverge to infinity.

Figure 4 shows the critical residue behavior of the up and down periodic orbits on the symmetry line $`s_1`$ along several different bifurcation curves. For lower period bifurcation curves, the residues first show signs of a six-cycle (to be discussed later in greater details), but then converge to $`\left|R_i\right|0.25`$. This is because the invariant torus we are studying is not quite shearless. Thus we see the same behavior of the residues as in the case of a twist map. As we proceed to higher period bifurcation curves, the behavior of the residues of the approximating periodic orbits found along the $`s_1`$ symmetry line resembles more and more a six-cycle. A renormalization group interpretation of these results is given in Sec. IV.

Finally, we found the critical point $`(a_c,b_c)`$ along the $`1/\gamma ^2`$-bifurcation curve to be the following:

$$a_c=0.425160543,b_c=0.9244636470355.$$
(11)
At the critical parameter values $`(a_c,b_c)`$, the residues of the down periodic orbits on $`s_1`$, which are equal to the residues of the up periodic orbits on $`s_2`$ because of the symmetry of the map (see Eq. (6)), converge to the six-cyclenote3  $`\{C_1,C_2,C_3,C_4,C_5,C_6\}`$, where

$$\begin{array}{cc}C_1=0.609\pm 0.005,\hfill & C_2=1.288\pm 0.002,\hfill \\ C_3=2.593\pm 0.005,\hfill & C_4=1.584\pm 0.008,\hfill \\ C_5=2.336\pm 0.006,\hfill & C_6=2.593\pm 0.005.\hfill \end{array}$$
(12)
The six-cycle can clearly be seen in Fig. 5 which shows the residues of the up and down periodic orbits at the critical point along the $`s_1`$ symmetry line.

The residues of the up periodic orbits on $`s_1`$ (and of the down periodic orbits on $`s_2`$) converge to the six-cycle $`\{D_1,D_2,D_3,D_4,D_5,D_6\}`$, where

$$\begin{array}{cc}D_1=1.584\pm 0.008,\hfill & D_2=1.288\pm 0.002,\hfill \\ D_3=2.630\pm 0.006,\hfill & D_4=0.609\pm 0.005,\hfill \\ D_5=2.336\pm 0.006,\hfill & D_6=2.630\pm 0.006.\hfill \end{array}$$
(13)
The residue convergence for other symmetry lines is shown in Table 4 where we denote by $`R_{u_i}`$ and $`R_{d_i}`$ the residues of the up and down periodic orbit on the symmetry line $`s_i`$.

Note that the six-cycle $`\{D_i\}`$ of $`R_{u_1}`$ and $`R_{d_2}`$ (respectively, the six-cycle $`\{C_i\}`$ of $`R_{u_2}`$ and $`R_{d_1}`$) is observed to be the same as that of $`R_{u_3}`$ and $`R_{d_4}`$ (respectively, $`R_{u_4}`$ and $`R_{d_3}`$) except it is shifted. The two six-cycles are related because of the symmetry of the map as follows: $`D_1=C_4`$, $`D_2=C_2`$, $`D_4=C_1`$, $`D_5=C_5`$, $`C_3=C_6`$, and $`D_3=D_6`$. It was numerically observed that $`C_6D_6`$, and therefore $`C_3D_3`$. Using these relations we see that there are only five independent residues which we take to be $`C_1`$,$`C_2`$,$`C_3`$,$`C_4`$, and $`C_5`$.

We compared the values of the residues at three different points along the $`1/\gamma ^2`$-bifurcation curve, one point below criticality, one at criticality, and one above criticality: $`(a_{},b_{})=(0.425160540,0.9244636195728)`$, $`(a_c,b_c)=(0.425160543,0.9244636470355)`$ and $`(a_+,b_+)=(0.425160545,0.9244636653440)`$, respectively. The numerical results for the $`C_i`$ are listed in Table 5.

We see that each element of the six-cycle tends to zero for $`(a_{},b_{})`$, to infinity for $`(a_+,b_+)`$, while it tends to the critical value at $`(a_c,b_c)`$.

Figure 6 clearly illustrates this behavior.

A comparison with the results of Ref. diego4,  shows that, within numerical accuracy, we found the same values for the residues of the six-cycle, but the sequence has shifted by two: $`C_1=H_3`$, $`C_2=H_4`$, $`C_3=H_5`$, $`C_4=H_6`$, $`C_5=H_1`$ and $`C_6=H_2`$, where $`H_i`$ denote the residues for the $`1/\gamma `$ case found in Ref. diego4, . A similar shift by two occurs for other symmetry lines.

#### III.2.2 Spatial scaling at criticality

As expected, the shearless curve exhibits scale invariance at criticality, which can be demonstrated explicitly by using symmetry line coordinatesdiego4  $`(\widehat{x},\widehat{y})`$ defined by $`\widehat{x}=xa(1y^2)/2`$ and $`\widehat{y}=yy_s`$. In these coordinates, the $`s_3`$ symmetry line becomes a straight line that intersects the shearless curve at the origin. We find that, in symmetry line coordinates, the shearless $`1/\gamma ^2`$ invariant torus at criticality remains invariant under a scale change $`(x,y)(\alpha ^{12}x,\beta ^{12}y)`$. This property is illustrated in Fig. 7.

As described in Ref. diego5, , we can find $`y_s`$ using

$$y_s=\underset{i\mathrm{}}{lim}\frac{y_{[2i+1]}y_{[2i+11]}y_{[2i1]}y_{[2i+13]}}{\left(y_{[2i+1]}y_{[2i1]}\right)\left(y_{[2i+13]}y_{[2i+11]}\right)}0.47253494777,$$
(14)
where $`y_{[n]}`$ denotes the $`y`$-coordinate of the periodic orbit $`[n]`$ along the $`s_3`$ symmetry line. To obtain the quoted value of $`y_s`$ we used $`i=10`$. We then obtained $`\alpha `$ and $`\beta `$ as follows:note4

$$\alpha =\underset{n\mathrm{}}{lim}\left|\frac{\widehat{x}_{[2i+1]}}{\widehat{x}_{[2i+13]}}\right|^{1/12}1.61759$$
(15)
and

$$\beta =\underset{n\mathrm{}}{lim}\left|\frac{\widehat{y}_{[2i+1]}}{\widehat{y}_{[2i+13]}}\right|^{1/12}1.65702,$$
(16)
where $`(\widehat{x}_{[n]},\widehat{y}_{[n]})`$ are symmetry line coordinates of the point of the periodic orbit $`[n]`$ that is the closest to the origin. Within numerical accuracy, these values are the same as in Ref. diego5, .

Further numerical analysis shows that periodic orbits also exhibit scaling behavior locally near the $`s_3`$ symmetry line. Figure 8 shows points of the periodic orbit $`[n]=[21]`$ (in symmetry line coordinates) and points of the periodic orbit $`[33]`$ with the $`x`$ and $`y`$ coordinates rescaled by $`\alpha ^{12}`$ and $`\beta ^{12}`$ respectively. The result suggests that periodic orbits remain invariant under a simultaneous spatial rescaling and shifting of the winding number by twelve from $`[n]`$ to $`[n+12]`$.

### III.3 Numerical accuracy

We conclude this section with comments about the numerical accuracy of the results.

1. Points on the $`m/n`$-bifurcation curves were found with an accuracy ranging between $`10^{12}`$ and $`10^{15}`$, where the larger value corresponds to larger periods. We obtain this measure of accuracy from the condition

$$F(y)=0,F^{}(y)=0\text{ and }F^{\prime \prime }(y)0,$$
(17)
as explained previously. The numbers quoted above are the values of $`F(y)`$ obtained at the numerically found minima in $`y`$.
2. Periodic orbits along the different symmetry lines around the critical point where found with an accuracy ranging between $`10^{15}`$ and $`10^{17}`$. Here, the criterion is the difference between the winding number $`m/n`$ of the periodic orbit of interest, and the winding number of the orbit that results when starting at the numerically found location of the periodic orbit on the respective symmetry line, and then iterating the map.
3. A criterion for the accuracy of the scaled bifurcation curve $`\mathrm{\Phi }_{1/\gamma ^2}`$ is harder to find, since the location of the actual curve is unknown. An upper bound on the error, though, should be the distance between the $`\mathrm{\Phi }_{1/\gamma ^2}`$ and $`\mathrm{\Phi }_{[37]}`$, since the latter definitely lies on the other side of $`\mathrm{\Phi }_{1/\gamma ^2}`$. This error was found to be approximately $`2\times 10^{13}`$.
4. A criterion for the accuracy of the critical point in parameter space, $`(a_c,b_c)`$ is even harder to define, since we cannot actually find the residues of existing orbits of all  periods, which is required to check if the six-cycle at that point continues ad infinitum. We believe that the value for $`a_c`$ is accurate up to $`1\times 10^9`$.
5. The uncertainties for the critical residues, quoted above in Eq. (12), were computed from the variation in numerical values of residues at $`a_c`$ for the three or four highest period orbits found (e.g. for $`C_1`$, using residues of $`[13]`$, $`[19]`$, $`[25]`$ and $`[31]`$).
6. If we evaluate the residues of the up and down periodic orbits along the $`s_1`$ and $`s_3`$ symmetry lines, then the residues on the other two symmetry lines can be constructed using symmetry arguments (see Sec. II.1). But, as a check of the numerical procedures, we independently evaluated the residues on all the four symmetry lines and confirmed the symmetry arguments.

## IV Renormalization group interpretation

In this section, we interpret the above results within the renormalization group framework. The analysis follows Refs. diego5,  closely, since, as expected, the residue behavior exhibits a six-cycle at criticality. But, because of the different winding number (i.e. different from $`1/\gamma `$) the renormalization group operator will have a different form.

Renormalization ideas have been used fruitfully in area-preserving maps and Hamiltonian flows. (See e.g. Refs. shenker, ; mackay1, ; greene4, ; escande, ; chandre3, ; abad, ; koch,  and references therein.) In contrast to mathematical KAM theory, which proves the existence of dense sets of invariant tori in regions of phase space, the renormalization group approach addresses the problem of the destruction of an invariant torus with a specific winding number under strong perturbation.

The following renormalization approach (see e.g. Refs. greene6, ; diego5, ) is based on the residue criterion (Sec. II.1). To study the breakup of an invariant torus of winding number $`\omega `$, we loosely represent the map $`M`$ as:

$$M=(R_1,R_2,R_3,\mathrm{}),$$
(18)
where the $`\{R_i\}`$ are the residues of the periodic orbits with winding numbers $`\{m_i/n_i\}`$, the convergents of $`\omega `$. For example, the integrable map will be represented by $`(0,0,0,\mathrm{})`$ because all the orbits are parabolic in that case. The key idea is to construct an operator $``$ that explores the infinite tail of (18) by mapping a map given by (18) to another map, $`(M)`$, represented by $`\left(M\right)=(\widehat{R}_1,\widehat{R}_2,\mathrm{},\widehat{R}_i,\mathrm{})`$, where $`\widehat{R}_i=R_{i+1}`$. This operation can be interpreted as a time renormalization since periodic orbits with large periods are transformed into periodic orbits with smaller periods, which amounts to a rescaling of time.

The residue criterion can now be rephrased in this framework:

1. If $`lim_n\mathrm{}^n(M)=(0,0,0,\mathrm{})`$, the invariant torus exists.
2. If $`lim_n\mathrm{}^n(M)=(\pm \mathrm{},\pm \mathrm{},\mathrm{})`$, the invariant torus is destroyed.
3. If $`lim_n\mathrm{}^n(M)`$ is a map for which the residues have finite, nonzero values, i.e. a map that is invariant under the action of $`^m`$ (a fixed point of $`^m`$) for some $`m>0`$, the invariant torus is at the threshold of destruction. Possible scenarios are the convergence of the residues to a fixed value or to a convergent subsequence.

There are two kinds of fixed points: simple fixed points and critical fixed points. In the case of area-preserving maps, we come to the following interpretation. A simple fixed point is an integrable map (all the residues are zero), and its basin of attraction contains all the maps for which the invariant torus exists. A critical fixed point is a map for which the invariant torus under consideration is at criticality. All the maps in its basin of attraction exhibit the same universal behavior at the critical breakup.

### IV.1 Renormalization group operator

Following the discussion in Refs. mackay1,  and diego5, , we use pairs of commuting maps because they provide a simple way to define the renormalization operators for invariant tori.

A pair of commuting maps is an ordered pair of maps, $`(U,T)`$, such that $`UT=TU`$. An orbit of a point $`(x,y)`$ generated by $`(U,T)`$ is the set of points $`\{U^mT^n(x,y)\}`$, where $`m`$ and $`n`$ are integers. A periodic orbit of period $`m/n`$ is an orbit for which $`U^mT^n(x_i,y_i)=(x_i,y_i)`$.

For the breakup of the invariant torus with winding number $`\omega =1/\gamma ^2`$, we define the renormalization group operator by

$$\left(\begin{array}{c}U\\ T\end{array}\right):=B\left(\begin{array}{c}U^1T^1\\ UT^2\end{array}\right)B^1.$$
(19)
As for the case of $`1/\gamma `$,diego5  this operator contains both time and space renormalization as follows:

The space renormalization is represented by the operator $`B`$, which rescales the $`(x,y)`$ coordinates, i.e. $`(x,y)B(x,y)`$ where

$$B=\left(\begin{array}{cc}r& 0\\ 0& s\end{array}\right).$$
(20)
At the critical fixed point studied in this paper, we see that $`r=\alpha `$ and $`s=\beta `$ given by Eqs. (15)-(16).

The time renormalization is, again, accomplished by the specific combination of the commuting maps. If $`(x,y)`$ is a periodic orbit of $`(U,T)`$ with winding number $`F_i/F_{i+2}`$, then $`B(x,y)`$ is a periodic orbit of $`(\widehat{U},\widehat{T})=(U,T)`$ with winding number $`F_{i1}/F_{i+1}`$, as can be verified as follows:

$`\widehat{U}^{F_{i1}}\widehat{T}^{F_{i+1}}B(x,y)`$ $`=`$ $`B(UT)^{F_{i1}}\left(UT^2\right)^{F_{i+1}}(x,y)`$
$`=`$ $`BU^{F_{i1}+F_{i+1}}T^{F_{i1}+2F_{i+1}}(x,y)`$
$`=`$ $`BU^{F_i}T^{F_{i+2}}(x,y)`$
$`=`$ $`B(x,y).`$
By induction, an orbit with winding number $`F_i/F_{i+2}`$ under $`(U,T)`$ is transformed into an orbit of $`^n(U,T)`$ with winding number $`F_{in}/F_{i+2n}`$.

### IV.2 Simple periodic orbit of $``$

We can find the integrable period-two orbit $`(U_\pm ,T_\pm )`$ of the renormalization operator (19) by requiring that $`(U_\pm ,T_\pm )=(U_{},T_{})`$. This two-cycle is given by the following pairs of maps:

$$U_\pm \left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{c}x\gamma ^2y^2/\gamma ^2\\ y\end{array}\right),T_\pm \left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{c}x+1\pm y^2\\ y\end{array}\right),$$
(21)
where the rescaling of the coordinates is given by

$$B=\left(\begin{array}{cc}\gamma & 0\\ 0& \pm \gamma \end{array}\right).$$
(22)
Using the definition $`U_\pm ^mT_\pm ^n(x,y)=(x,y)`$ of the periodic orbits of period $`m/n`$, we get the rotation number as a function of $`y`$:

$`\omega _\pm (y)`$ $`=`$ $`{\displaystyle \frac{1\pm y^2}{\gamma ^2y^2/\gamma ^2}}={\displaystyle \frac{1}{\gamma ^2}}\left(1\pm y^2\right)\left(1\pm {\displaystyle \frac{y^2}{\gamma ^2}}\right)^1`$
$``$ $`{\displaystyle \frac{1}{\gamma ^2}}\left[1\pm \left(1{\displaystyle \frac{1}{\gamma ^4}}\right)y^2+\mathrm{}\right].`$
Thus we see that the map $`(U_{},T_{})`$ is locally equivalent, under a change of coordinates, to the SNM with parameters $`(a,b)=(1/\gamma ^2,0)`$.

### IV.3 Critical periodic orbit of $``$

The next step is to analyze the critical periodic orbit of $``$. Consider the nontwist map

$$𝒞=(C_1,,C_2,,C_3,,C_4,,C_5,,C_6,,C_1,,C_2,\mathrm{}),$$
(23)
where the $`C_i`$ are the elements of the six-cycle computed earlier, and the “-” denote the periodic orbits that do not exist (see Table 2). By construction, this map is a period-12 orbit of the renormalization group operator (a fixed point of $`^{12}`$), i.e.

$$^{12}𝒞=𝒞.$$
(24)
In Sec. III.2, we found that the residues of the periodic orbits approximating the $`1/\gamma ^2`$-shearless curve in the standard nontwist map exhibit convergence to the six-cycle $`\left\{C_i\right\}`$. Assuming that we can fine-tune the results for $`(a_c,b_c)`$, we expect that $`lim_n\mathrm{}^nM(a_c,b_c)=𝒞`$.

If we are studying the breakup of the $`1/\gamma ^2`$-shearless curve for parameter values along the bifurcation curve for one of the low-period convergents, then we start near the stable manifold of the critical periodic orbit of $``$. But, under the action of $``$, we are pushed along an unstable direction. Thus, we see parts of the six-cycle of residues (see Fig. 4), but the limiting residue behavior is observed to be $`lim_i\mathrm{}\left|R_i\right|0.25`$, which is characteristic for the critical fixed point of twist maps (see e.g. Ref. mackay1, ). In renormalization group language, this means that part of the unstable manifold of the critical nontwist fixed point (maps for which $`(a,b)`$ is below $`(a_c,b_c)`$) is in the basin of attraction of the critical twist fixed point.

### IV.4 Eigenvalues

As shown in Ref. diego5, , it is possible to use our numerical data to draw further conclusions about the renormalization group operator $``$, in particular to compute its unstable eigenvalues. The main difficulty in computing these eigenvalues is that the space of maps is infinite-dimensional whereas the $`(a,b)`$ parameter space has only two dimensions. But the fact that we can find an isolated point $`(a_c,b_c)`$ in parameter space at which the map is at criticality means that the dimension of the unstable manifold is two. The map $`M`$ at $`(a_c,b_c)`$ is the intersection point of the two-parameter family of maps with the stable manifold of the fixed point, i.e. values $`a_c`$ and $`b_c`$ describe the location of the critical fixed point of $`^{12}`$ in its unstable manifold.

Below, we first compute the eigenvalues characterizing the approach to $`(a_c,b_c)`$ in the $`(a,b)`$ parameter space using the numerical results from above. As shown in Ref. diego5, , based on the type of numerical data obtained, the two eigenvalues can be found by

$$\nu _1=\underset{n\mathrm{}}{lim}\left(\frac{\mathrm{\Phi }_{[n+12]}\left(a_c\right)b_c}{\mathrm{\Phi }_{[n]}\left(a_c\right)b_c}\right)$$
(25)
and

$$\nu _2=\underset{n\mathrm{}}{lim}\left(\frac{a_{c[2n+12]}a_c}{a_{c[2n]}a_c}\right).$$
(26)
The last step is to relate the values $`\nu _i`$ to the unstable eigenvalues $`\delta _i`$ of the renormalization group operator $``$. The key idea is to study the behavior of the residues of the periodic orbits approximating $`\omega `$ at the $`(a,b)`$ values used in the computation of $`\nu _i`$. For details see Ref. diego5, .

Denoting the unstable eigenvalues of $``$ by $`\delta _1`$ and $`\delta _2`$, we conclude that $`\delta _i=\left(1/\nu _i\right)^{1/12}`$. We find numerical values of

$$\delta _12.678,\delta _21.583.$$
(27)
Comparing this to the values found in Ref. diego5,  shows that within (assumed) numerical uncertainty these values are the same as those for $`1/\gamma `$, as predicted. The $`a_c`$ values used to determine $`\delta _2`$ were $`a_{c[26]}`$ and $`a_{c[14]}`$, which explains the larger discrepancy. Work is under way to improve this result.

## V Conclusion

We have shown through numerical simulations that the critical residue values at the breakup of the $`1/\gamma ^2`$-shearless curve in the standard nontwist map coincide with those of the $`1/\gamma `$-shearless curve. In addition, the critical scaling parameters and the unstable eigenvalues of the renormalization group operator were found to be the same for both cases. The main differences are the location of the respective critical point in parameter space and the detailed form of the renormalization group operator in terms of commuting maps pairs.

Future work includes the search for the breakup values of more winding numbers to map out the details of the critical function depicted in Fig. 1. In addition, new fixed points of the renormalization group operator might be obtained by studying the breakup of shearless curves with non-noble winding numbers.

###### Acknowledgements.

The authors would like to thank John Greene and Diego del-Castillo-Negrete for many helpful discussions. This research was supported in part by U.S. Department of Energy Contract No. DE-FG01-96ER-54346 and by an appointment of A. Wurm to the U.S. Department of Energy Fusion Energy Postdoctoral Research Program administered by the Oak Ridge Institute for Science and Education.