text:
# The differential sum rule for the relaxation rate in the cuprates

## Abstract

Motivated by recent experiments by Basov et al, we study the differential sum rule for the effective scattering rate $`1/\tau (\omega )`$. We show that in a dirty BCS superconductor, the area under $`1/\tau (\omega )`$ does not change between the normal and the superconducting states. For magnetically mediated pairing, similar result holds between $`T<T_c`$ and $`TT_c`$, while in the pseudogap phase, $`1/\tau (\omega )`$ is just suppressed compared to $`1/\tau (\omega )`$ in the normal state. We argue that this violation of the differential sum rule is due to the fact that low-energy fermions in the pseudogap phase are still incoherent.

The analysis of the optical sum rules in condensed matter systems is a valuable tool that helps one to understand the key physics and relevant energy scales in the problem tom . The focus of this report are recent experimental results basov  for the effective relaxation rate $`\tau ^1(\omega )=(4\pi /\omega _{pl}^2)Re[1/\sigma (\omega )]`$ where $`\sigma (\omega )=\sigma _1(\omega )+i\sigma _2(\omega )`$ is the optical conductivity, $`\omega _{pl}^2=4\pi ne^2/m`$ is the plasma frequency, and $`n`$ is the density of particles. The data analysis for optimally doped $`YBa_2Cu_3O_{6.95}`$ basov\_1  and $`Tl_2Ba_2CuO_{6+x}`$ puchkov  and underdoped $`YBa_2Cu_4O_8`$ and $`Bi_2Sr_2CaCu_2O_{8+x}`$ basov  revealed an approximate differential sum rule for $`\tau ^1(\omega )`$ between $`TT_c`$ and $`T<T_c`$: although $`𝑑\omega \tau ^1(\omega )`$ does not converge, it changes very little when the system enters into the superconducting state. This differential sum rule, however, is not satisfied between the normal and the pseudogap phases – $`1/\tau (\omega )`$ in the pseudogap phase appears to be just suppressed.

Quite generally, exact sum rules are related to conservation laws. The $`f`$sum rule for the optical conductivity states that at a given density of particles, the total absorbing power of the solid characterized by $`\sigma _1`$ does not depend on the details of the interactions and is determined only by the total number of particles in the system pi\_noz . The total absorption power is given by $`_0^{\mathrm{}}\sigma _1(\omega )`$ Applying the Kubo formula that relates $`\sigma (\omega )`$ with the full retarded current-current correlator $`\mathrm{\Pi }(\omega )`$: $`\sigma (\omega )=(\omega _{pl}^2/4\pi )\mathrm{\Pi }(\omega )/(i\omega +\delta )`$, separating the frequency integral into the integral over infinitesimally small $`\omega `$ and the rest and using the Kramers-Kronig relation for $`\mathrm{\Pi }(\omega )1`$ that vanishes at the highest frequencies, we obtain $`_0^{\mathrm{}}\sigma _1(\omega )=\omega _{pl}^2/8`$, independent of the form of $`\mathrm{\Pi }(\omega )`$.

Is there an analogous sum rule for $`1/\tau (\omega )`$? Using $`1/\tau (\omega )=Im[\omega ^2/\mathrm{\Pi }(\omega )]/\omega `$ and applying the Kramers-Kronig relation, we find

$$_0^{\mathrm{}}\frac{d\omega }{\tau (\omega )}=\frac{\pi }{2}\left[\mathrm{}\frac{\omega ^2}{\mathrm{\Pi }(\omega )}_{\omega 0}+C\right]=\frac{\pi }{2}C$$
(1)
The constant $`C`$ again has to be chosen such that $`\omega ^2/\mathrm{\Pi }(\omega )+C`$ vanishes at $`\omega \mathrm{}`$. However, $`C`$ turns out to be infinite as at high frequencies $`\mathrm{\Pi }(\omega )1`$, and $`\omega ^2/\mathrm{\Pi }(\omega )`$ diverges. This divergence implies that there is no conservation law associated with the relaxation rate and hence no sum rule for $`1/\tau (\omega )`$.

Still, one can try to see whether one can get useful information by comparing $`1/\tau (\omega )`$ for two different system parameters, e.g. temperatures, which do not affect the system behavior at high frequencies. Suppose momentarily that $`1/\tau (\omega ,T_1)1/\tau (\omega ,T_2)`$ vanishes at high frequencies. The Kramers-Kronig relation is then fully satisfied and using (1) we immediately find that $`𝑑\omega (\tau ^1(\omega ,T_1)\tau ^1(\omega ,T_2))=0`$, i.e., the area under $`1/\tau (\omega )`$ does not change with $`T`$. The conservation of the total area under $`1/\tau (\omega )`$ would create a valuable tool to study the evolution of the spectral weight between, e.g., the normal and the superconducting states. This new differential sum rule, however, is not associated with a conservation law and therefore is not guaranteed to be satisfied – only explicit calculations can determine whether or not the temperature dependence in $`\mathrm{\Pi }(\omega ,T)`$ is weak enough to ensure the convergence of $`1/\tau (\omega ,T_1)1/\tau (\omega ,T_2)`$.

In this communication we study under which conditions the differential sum rule for $`1/\tau (\omega )`$ actually holds, and at which frequencies it is exhausted. We consider the magnetic scenario for the pairing in the cuprates, and argue that although the differential sum rule is approximately satisfied for all $`T`$, it is exhausted at energies comparable to the pairing gap only if there is a strong feedback effect from the pairing on the fermionic propagator – otherwise, the sum rule is exhausted at much larger frequencies. We associate the first regime with $`T<T_c`$, and the second one with the pseudogap phase.

To put our analysis of the spin mediated pairing into perspective we first analyze the situation in a dirty BSC $`s`$wave superconductor at $`T=0`$, when the pairing causes a strong feedback on the fermionic propagator, and in the toy model where there is no feedback from the pairing on the fermionic self-energy.

The theory of a dirty superconductor is well developed ag ; anderson . In the normal state, the inelastic scattering by impurities yields a retarded fermionic self-energy $`\mathrm{\Sigma }(\omega )=i/2\tau `$, In a superconducting state, this self-energy is modified due to a feedback from superconductivity and takes the form $`\mathrm{\Sigma }(\omega )=(i/2\tau )\omega /\sqrt{\omega ^2\mathrm{\Delta }^2}`$ where $`\mathrm{\Delta }`$ is the superconducting gap ag . Substituting these forms into the current-current polarization bubble and performing the momentum integration we obtain

$`\mathrm{\Pi }(\omega )={\displaystyle _0^{\mathrm{}}}𝑑\mathrm{\Omega }{\displaystyle \frac{1}{\left(\sqrt{\mathrm{\Omega }_+^2\mathrm{\Delta }^2}+\sqrt{\mathrm{\Omega }_{}^2\mathrm{\Delta }^2}+i/\tau \right)}}`$
$`\times {\displaystyle \frac{\sqrt{\mathrm{\Omega }_+^2\mathrm{\Delta }^2}\sqrt{\mathrm{\Omega }_{}^2\mathrm{\Delta }^2}\mathrm{\Delta }^2\mathrm{\Omega }_+\mathrm{\Omega }_{}}{\sqrt{\mathrm{\Omega }_+^2\mathrm{\Delta }^2}\sqrt{\mathrm{\Omega }_{}^2\mathrm{\Delta }^2}}}`$ (2)
where $`\mathrm{\Omega }_\pm =\mathrm{\Omega }\pm \omega /2`$. In the normal state, this reduces to a conventional Drude form $`\mathrm{\Pi }(\omega )=\omega /(\omega +i/\tau )`$. In the superconducting state, the frequency integral in (2) can be evaluated analytically in the clean limit $`\mathrm{\Delta }\tau 1`$. After lengthy but straightforward calculations we found that both $`\sigma _1(\omega )`$ and $`1/\tau (\omega )`$ vanish at $`\omega <2\mathrm{\Delta }`$, while at larger frequencies

$$\sigma _1(\omega )=\frac{\omega _{pl}^2}{4\pi \tau \omega ^2}E\left(\sqrt{1\frac{4\mathrm{\Delta }^2}{\omega ^2}}\right),\frac{1}{\tau (\omega )}=\frac{4\pi \sigma _1(\omega )\omega ^2}{\omega _{pl}^2}$$
(3)
where $`E(x)`$ is the complete elliptic integral of the second kind morr . At $`\omega =2\mathrm{\Delta }+0`$, $`E=\pi /2`$ and both $`\sigma _1(\omega )`$ and $`1/\tau (\omega )`$ jump to finite values. At high frequencies, $`E(x1)1`$, $`\sigma _1(\omega )`$ vanishes as $`\omega ^2`$, and $`1/\tau (\omega )`$ approaches the normal state result $`\tau (\omega )=\tau `$. To the same order, we also have $`\mathrm{\Pi }(0)=1\pi /(8\mathrm{\Delta }\tau )`$. We checked analytically that the $`f`$sum rule $`(8/\omega _{pl}^2)_{+0}^{\mathrm{}}𝑑\omega \sigma _1(\omega )=1\mathrm{\Pi }(0)`$ is indeed satisfied.

Expanding $`E(x)`$ near $`x=1`$, we find that at high frequencies $`\tau ^1(\omega )\tau ^1(2\mathrm{\Delta }^2/\omega ^2\tau )(\mathrm{log}(2\omega /\mathrm{\Delta })0.5)`$, i.e. $`𝑑\omega (1/\tau (\omega )1/\tau )`$ converges. The convergence implies that for a dirty BCS superconductor, the differential sum rule for $`1/\tau (\omega )`$ is an exact one.

The plots of $`\sigma _1(\omega )`$ and $`1/\tau (\omega )`$ are presented in Fig. 1 together with the results for $`I_\sigma (\omega )=(8/\omega _{pl}^2)_0^\omega 𝑑x\sigma _1(x)`$ and $`I_\tau (\omega )=(\tau /2\mathrm{\Delta })_0^\omega 𝑑x(1/\tau (x)1/\tau )`$. It follows from (3) that the sum rules for $`\sigma _1(\omega )`$ and $`1/\tau (\omega )`$ should be exhausted at frequencies of the order of $`\mathrm{\Delta }`$. Fig. 1 shows that this is indeed the case.

We next consider the behavior of $`\sigma _1(\omega )`$ and $`1/\tau (\omega )`$ in the toy model in which we assume that the pairing does not change the fermionic self-energy. Having the cuprates in mind, we also assume, that the normal state is not a Fermi liquid, i.e., fermionic self-energy at low frequencies behaves as $`\mathrm{\Sigma }(\omega )=(i\omega )^\alpha \overline{\omega }^{1\alpha }`$ with $`\alpha <1`$,

For the frequency dependent self-energy $`\mathrm{\Sigma }(\omega )`$, the current-current correlator $`\mathrm{\Pi }(\omega )`$ is still given by Eq. (2), but now $`\mathrm{\Omega }`$ is substituted by $`\mathrm{\Omega }+\mathrm{\Sigma }(\mathrm{\Omega })`$. For definiteness, we present the results for $`\alpha =1/2`$, which is the normal state quantum-critical exponent in the spin-fermion theory abanov01 , but the results are qualitatively the same for all $`\alpha `$ including the marginal Fermi liquid limit $`\alpha 0`$ av . For $`\alpha =1/2`$, the fermionic density of states $`N(\omega )=Im[\stackrel{~}{\mathrm{\Sigma }}(\omega )/(\mathrm{\Delta }^2\stackrel{~}{\mathrm{\Sigma }}^2(\omega ))^{1/2}]`$ has a maximum at $`\omega =\stackrel{~}{\mathrm{\Delta }}\mathrm{\Delta }^2/\overline{\omega }`$, that could be identified with the pseudogap. Still, $`N(\omega )`$ remains finite at all frequencies. The assumption that the self-energy does not change is internally consistent if indeed the anomalous normal state behavior extends to frequencies which well exceed the pairing gap, i.e., when $`\overline{\omega }\stackrel{~}{\mathrm{\Delta }}`$.

Evaluating $`\mathrm{\Pi }(\omega )`$ and substituting it into $`\sigma _1(\omega )`$ and $`1/\tau (\omega )`$ we found that in the normal state, $`\sigma _{1,n}(\omega )(\omega \overline{\omega })^{1/2}`$ at $`\omega \overline{\omega }`$ and $`\sigma _{1,n}(\overline{\omega }/\omega ^3)^{1/2}`$ at $`\omega \overline{\omega }`$, while $`1/\tau _n(\omega )(\omega \overline{\omega })^{1/2}`$ in both limits. For $`\mathrm{\Delta }0`$, we found that for $`\omega \stackrel{~}{\mathrm{\Delta }}`$, $`\sigma _1(\omega )(\omega \overline{\omega })^{1/2}(\omega /\stackrel{~}{\mathrm{\Delta }})^{5/2}`$ and $`1/\tau (\omega )(\omega \overline{\omega })^{1/2}(\omega /\stackrel{~}{\mathrm{\Delta }})^{7/2}`$. We see that both $`\sigma _1(\omega )`$ and $`1/\tau (\omega )`$ are reduced compared to their normal state values but are still finite. At $`\omega \stackrel{~}{\mathrm{\Delta }}`$, both $`\sigma _1(\omega )`$ and $`1/\tau (\omega )`$ become comparable to the normal state results. At larger $`\stackrel{~}{\mathrm{\Delta }}\omega \overline{\omega }`$, $`\sigma _1(\omega )=\sigma _{1,n}(\omega )1.992(\omega _{pl}^2/4\pi )\stackrel{~}{\mathrm{\Delta }}(\omega ^3\overline{\omega })^{1/2}`$ and $`1/\tau (\omega )=1/\tau _n(\omega )3.51\stackrel{~}{\mathrm{\Delta }}(\overline{\omega }/\omega )^{1/2}`$. Finally, at $`\omega \overline{\omega }`$, $`\sigma _1(\omega )\sigma _{1,n}(\omega )\omega ^{7/2}\mathrm{log}\omega `$ and $`1/\tau (\omega )1/\tau _n(\omega )\omega ^{3/2}\mathrm{log}\omega `$.

We see that $`\sigma _1(\omega )`$ converges to its normal state value at frequencies of order $`\stackrel{~}{\mathrm{\Delta }}`$, as in a dirty BCS superconductor; the sum rule for $`\sigma _1(\omega )`$ is then exhausted at $`\omega \stackrel{~}{\mathrm{\Delta }}`$. This behavior is illustrated in Fig. 2a where we present the results of our numerical calculations – $`I_\sigma (\omega )`$ converges to $`I_\sigma (\mathrm{})=1\mathrm{\Pi }(0)`$ ($`0.67`$ for our choice of $`\overline{\omega }=2\mathrm{\Delta }`$) already at $`\mathrm{\Omega }\stackrel{~}{\mathrm{\Delta }}`$. On the other hand, $`\tau ^1(\omega )\tau _n^1(\omega )`$ scales as $`\omega ^{1/2}`$ between $`\omega \stackrel{~}{\mathrm{\Delta }}`$ and $`\omega \overline{\omega }`$ such that at these frequencies $`I_\tau (\omega )=𝑑\omega (\tau ^1(\omega )\tau _n^1(\omega ))`$ does not converge. Furthermore, at these frequencies, $`1/\tau (\omega )`$ is still smaller than $`1/\tau _n(\omega )`$. This result holds for all $`\alpha <1`$ as one can straightforwardly verify. Only at $`\omega >\overline{\omega }`$, $`\tau ^1(\omega )`$ finally becomes larger than $`\tau _n^1(\omega )`$, and $`I_\tau (\omega )`$ converges. The convergence implies that the differential sum rule for $`1/\tau (\omega )`$ is again exactly satisfied, however it is exhausted only at frequencies that well exceed the pseudogap. We present the numerical results for $`1/\tau (\omega )`$ and $`I_\tau (\omega )`$ in Fig. 2b.

We now present the results for $`\sigma _1(\omega )`$ and $`1/\tau (\omega )`$ for spin-fluctuation mediated $`d`$wave pairing. We obtained these results by solving a set of coupled Eliashberg equations for the spin-fermion model that describes the spin fluctuation exchange at low energies acs . We will demonstrate that at low $`T`$, the behavior of the conductivity and the relaxation rate resembles that in a dirty BSC superconductor, while immediately below the pairing instability the system behavior is similar to that in the toy model for the pseudogap.

The spin-fermion model is characterized by a single dimensionless coupling constant $`\lambda `$ and a single overall energy $`\overline{\omega }`$ that scales with the effective spin-fermion interaction abanov01 . We will also use a characteristic energy scale for the spin fluctuations $`\omega _{sf}=\overline{\omega }/4\lambda ^2`$. A fit to the NMR, ARPES and neutron experiments yields $`\lambda 12`$ near the optimal doping abanov01 . We refer the readers to Ref abanov01  for the discussion of the applicability of the model to the cuprates and the justification of the Eliashberg approach at strong spin-fermion coupling despite the formal absence of the Migdal theorem.

We begin with the normal state. In Fig. 3a we present our results for $`1/\tau (\omega )`$ and $`I_\tau (\omega )`$ at various $`T`$. For definiteness we used $`\lambda =2`$. We see that $`I_\tau (\omega )`$ diverges at high frequencies, i.e., the differential sum rule is not satisfied. We checked analytically that this is caused by the $`1/\omega `$ behavior of the integrand in $`I_\tau (\omega )`$.

For completeness, in Fig. 3b we present the results for $`\sigma _1(\omega )`$ and $`I_\sigma (\omega )`$ at various $`T`$. We see that $`I_\sigma (\omega )`$ flattens up at $`\omega 10\omega _{sf}`$, but its value is still by about $`30\%`$ smaller than it should be for $`\omega =\mathrm{}`$. The full sum rule is exhausted only at unrealistically large $`\omega 10^3\omega _{sf}`$ (Ref. abanov01 ) where the low-energy theory is clearly inapplicable. The weak convergence of $`I_\sigma (\omega )`$ is related to the fact that over a wide frequency range $`\sigma _1(\omega )`$ is inversely proportional to $`\omega `$, and $`I_\sigma (\omega )`$ increases as $`\mathrm{log}\omega `$. The practical application of this result is that the plasma frequency $`\omega _{pl}`$ obtained by integrating the measured $`\sigma _1(\omega )`$ up to, say, $`1eV`$ should be about $`30\%`$ less than the actual plasma frequency. This is consistent with the fact that to match the measured $`\sigma _1(\omega 0)`$ with the theoretical $`\sigma _1(\omega 0)`$ obtained using $`\mathrm{\Sigma }(\omega )`$ extracted from the ARPES data, one needs the plasma frequency of about $`20000cm^1`$ abanov01 ; millis , while $`\omega _{pl}`$ extracted from the “partial” sum rule puchkov  is around $`16000cm^1`$.

We next consider what happens below the pairing instability temperature $`T_{ins}0.2\overline{\omega }`$ acf . Earlier Schmalian and the two of us have found acs  that at $`TT_{ins}`$, the the fermionic self-energy remains large at the smallest $`\omega `$ and gradually evolves from its normal state value. This gradual behavior is qualitatively different from a dirty BCS superconductor, as in the latter the quasiparticle spectral function and the density of states instantly drop to zero at frequencies below $`\mathrm{\Delta }`$ due to a feedback from the pairing ag ; anderson . We found in acs  that the behavior similar to that in a BSC superconductor occurs only below $`T_c<T_{ins}`$, and the difference between $`T_{ins}`$ and $`T_c`$ increases with increasing $`\lambda `$. We conjectured that at $`T_c<T<T_{ins}`$, fluctuations destroy coherent superconductivity, i.e., the system is in the pseudogap regime.

In Fig. 4 we present the results for $`\sigma _1(\omega )`$ and $`1/\tau (\omega )`$ for two different $`\lambda `$ and for three different temperatures: $`TT_c`$, $`T_c<T<T_{ins}`$ and $`T=T_{ins}`$ where $`\sigma _1(\omega )`$ and $`1/\tau (\omega )`$ are the same as in the normal state. We clearly see that between $`T_c<T<T_{ins}`$, $`1/\tau (\omega )`$ is nearly homogeneously suppressed, while at smaller $`T`$ it shifts from the lowest frequencies to $`\omega 2\mathrm{\Delta }`$ ($`0.8\overline{\omega }`$) and develops an overshoot at $`\omega 2\mathrm{\Delta }`$. The magnitude of the overshoot depends on the coupling and is larger at larger $`\lambda `$. We compared our results with the fermionic self-energy acs , and found that the temperature where $`1/\tau (\omega )`$ develops an overshoot, agrees with the onset temperature for the feedback effects on the fermionic propagator.

For quantitative analysis of the sum rule we looked analytically into the high frequency behavior of $`1/\tau (\omega )`$. We found that the temperature variation of $`1/\tau (\omega )`$ at high frequencies is still strong enough to prevent the differential sum rule for $`1/\tau `$. On the other hand, we see from the insert in Fig. 4b that this variation can be safely neglected if the frequency integration is restricted to $`\omega 34\overline{\omega }10\mathrm{\Delta }`$ (for optimally doped $`YBCO`$, $`10\mathrm{\Delta }2500cm^1`$). In Fig. 4b we show the results for the differential sum rule between $`TT_c`$ and $`T_c<T<T_{ins}`$. We see that the spectral weight is almost exactly conserved if the upper limit of the frequency integral is chosen close to $`34\overline{\omega }`$ . The conservation of the spectral weight between $`TT_c`$ and $`TT_c`$, and the apparent loss of the spectral weight between $`T_c`$ and $`T_{ins}`$ are the main results of recent experimental analysis basov . We see that these results are fully reproduced in our analysis.

We also see from Fig. 4 that at small $`\omega `$, $`\sigma _1(\omega )`$ keeps increasing between $`T_{ins}`$ and $`T_c`$. This indicates that the development of the pseudogap does not give rise to a suppression of the conductivity at the lowest frequencies. The latter is only reduced below $`T_c`$. The sensitivity of $`\sigma _1(\omega =0)`$ to $`T_c`$ rather than to the pseudogap temperature is consistent with the data nicole . The low frequency behavior of $`\sigma _1(\omega =0)`$ below $`T_c`$ is not captured in our theory as it is predominantly determined by impurities lee . Finally, we found both analytically and numerically that at large $`\omega `$, $`\sigma _1(\omega )`$ also turns out to be sensitive to $`T_c`$ (see the inset in Fig. 4a). This also agrees with the data vdm .

To conclude, in this paper we considered the differential sum rule for the effective scattering rate $`1/\tau (\omega )`$ (the difference between the area under $`1/\tau (\omega )`$ for two different temperatures). We argued that for spin-fluctuation mediated pairing, this sum rule is generally not an exact one, but is rather well satisfied below $`T_c`$ and is exhausted at frequencies compared to the pairing gap, $`\mathrm{\Delta }`$. We identified this behavior with the strong feedback from the pairing on the fermionic self-energy. We found that in the pseudogap region, where feedback effects are small, the differential sum rule is (approximately) exhausted only at much larger energies, comparable to the upper limit of non-Fermi-liquid behavior in the normal state, while at $`\omega =O(\mathrm{\Delta })`$, $`1/\tau (\omega )`$ is nearly homogeneously suppressed compared to the normal state. We argued that this behavior as well as the behavior of $`\sigma _1(\omega )`$ are consistent with the experimental data on the cuprates.

It is our pleasure to thank D. N. Basov, G. Blumberg, N. Bontemps, C. Homes, D. van der Marel, A. Millis, M. Norman, D. Pines, and J. Tu for useful conversations. We are also thankful to D. N. Basov, N. Bontemps, and J. Tu for sharing the unpublished results with us. The research was supported by NSF DMR-9979749 (A. Ch.) and by DR Project 200153 at Los Alamos National Laboratory (Ar. A.).


format: text

text:
# Periodic orbits for three particles with finite angular momentum

## Introduction

Recently C. Moore  applied the theory of braids to find a classification of periodic orbits for the motion of several particles moving in a plane under the action of pairwise attractive power-law forces, obtaining numerical results for the case of 2 and 3 strands. In particular, for three particles of equal mass, he found a solution where these particles travel on a common orbit in the shape of a figure eight. This orbit was found also independently by A. Chenciner and R. Montgomery , who have proved its existence for the case of zero total angular momentum and inverse square forces. Moore’s numerical result that this orbit is stable has also been confirmed by computations of C Simó . Previously, one of the simplest known periodic orbits without collisions for this problem  was Lagrange’s classic solution  where the three particles form the vertices of an equilateral triangle with each particle moving on an elliptic orbit with focus in the common center of mass. In this paper we show that for attractive power-law forces of the form $`1/r^{p+1}`$ with $`2<p0`$, the circular Lagrange orbit is the limit of a family of orbits which are periodic in a uniformly rotating frame of reference. For small deviations, we find analytically that this orbit is a distorted ellipse which, as the angular momentum decreases, resembles a distorted figure eight shape, but it does not approach the recently discovered zero-angular-momentum figure eight orbit which has a different winding number. For this figure eight orbit, which exists for $`2<p`$, we obtain an analytic approximation, and we prove by perturbation theory that it can be extended to finite angular momentum. We also investigate the evolution of these two orbits for a finite range of total angular momentum by a novel iterative method, and consider briefly the case of such periodic orbits for n equal mass particles, and for three particles of different masses.

## Computational method for periodic orbits

Our method is based on expanding in a Fourier series the Cartesian coordinates for the periodic orbit on which all the three equal mass particles travel, and then finding the extremun of the action for this system with respect to the Fourier coefficients in a frame rotating with a fixed frequency $`\mathrm{\Omega }`$. We obtain analytic approximations by perturbation theory, and numerical solutions by a novel iterative method which converges rapidly. Setting $`\theta =\omega t`$, where $`t`$ is the time, $`\omega =2\pi /\tau `$ is the frequency, and $`\tau `$ is the period of the orbit), we have

$$x(\theta )=\underset{k=1}{\overset{k=n}{}}a_s(k_o)sin(k_o\theta )+a_c(k_e)cos(k_e\theta )$$
(1)
and

$$y(\theta )=\underset{k=1}{\overset{k=n}{}}b_s(k_e)sin(k_e\theta )+b_c(k_o)cos(k_o\theta ),$$
(2)
where $`k_o=2k1`$ and $`k_e=2k`$, excluding multiples of the integer $`3`$. In practice, the number of terms $`n`$ in the Fourier sums, Eqs. 1 and 2, depends on the desired accuracy of the solution. The Cartesian coordinates of each of the three particles is then given by $`x_i=x(\theta _i)`$ and $`y_i=y(\theta _i)`$, for $`i=1,2,3`$, and the invariance of the action under the cyclic transformation $`1231`$ is satisfied by the requirement $`\theta _1=\theta `$, $`\theta _2=\theta +2\pi /3`$ and $`\theta _3=\theta +4\pi /3`$. These conditions are satisfied by the equations which determine the extremun of the action in a frame rotating with frequency $`\mathrm{\Omega }`$ in the plane of the orbit  which takes the form

$$A=_0^{2\pi }𝑑\theta [K(\theta )P(\theta )+\frac{1}{2}\mathrm{\Omega }^2I(\theta )+\mathrm{\Omega }L(\theta )],$$
(3)
where $`K(\theta )`$ is the kinetic energy, $`P(\theta )`$ the potential energy, $`I(\theta )`$ the moment of inertia, and $`L(\theta )`$ the angular momentum in the rotating frame. These are given by the following expressions:

$$K(\theta )=\frac{1}{2}\omega ^2\underset{i=1}{\overset{i=3}{}}(\frac{dx_i}{d\theta })^2+(\frac{dy_i}{d\theta })^2,$$
(4)
$$P(\theta )=(\frac{1}{p})(\frac{1}{r_{12}^p}+\frac{1}{r_{13}^p}+\frac{1}{r_{23}^p}),$$
(5)
$$I(\theta )=\underset{i=1}{\overset{i=3}{}}x_i^2+y_i^2,$$
(6)
and

$$L(\theta )=\omega \underset{i=1}{\overset{i=3}{}}x_i\frac{dy_i}{d\theta }y_i\frac{dx_i}{d\theta },$$
(7)
where $`r_{ij}=\sqrt{(x_ix_j)^2+(y_iy_j)^2}`$ is the distance between particles $`i`$ and $`j`$, and $`p`$ is the power of the force law,

The extrema of the action, Eq. 3, are obtained by finding the Fourier coefficients for which the partial derivatives of the action with respect to these coefficients vanish,

$$\frac{A}{a_s(k_o)}=u(k_0)a_s(k_o)f_s(k_o)v(k_0)b_c(k_o)=0,$$
(8)
$$\frac{A}{b_c(k_o)}=u(k_0)b_c(k_o)g_c(k_o)v(k_0)a_s(k_o)=0,$$
(9)
$$\frac{A}{b_s(k_e)}=u(k_0)b_s(k_e)g_s(k_e)+v(k_0)a_c(k_e)=0,$$
(10)
$$\frac{A}{a_c(k_e)}=u(k_0)a_c(k_e)f_c(k_e)+v(k_0)b_s(k_e)=0,$$
(11)
where $`u(k)=k^2\omega ^2+\mathrm{\Omega }^2`$ and $`v(k)=2k\omega \mathrm{\Omega }`$. Here the functions

$$f_s(k_o)=\frac{1}{\pi }_0^{2\pi }𝑑\theta \frac{P}{x_1}sin(k_o\theta ),$$
(12)
$$f_c(k_e)=\frac{1}{\pi }_0^{2\pi }𝑑\theta \frac{P}{x_1}cos(k_e\theta ),$$
(13)
$$g_s(k_e)=\frac{1}{\pi }_0^{2\pi }𝑑\theta \frac{P}{y_1}sin(k_e\theta ),$$
(14)
$$g_c(k_o)=\frac{1}{\pi }_0^{2\pi }𝑑\theta \frac{P}{y_1}cos(k_o\theta ).$$
(15)
depend also on all the Fourier coefficients.

We obtain numerical solutions of these equations by starting with the Fourier coefficients to some conjectured approximate period orbit, e.g. according to Moore’s classification , and then evaluating the partial derivatives of the action, Eqs. 8-11. If the extrema of the action is at a maximum (minimum), an improved orbit is then obtained by changing each coefficient in proportion to plus(minus) the corresponding partial derivative of the action. Thus, starting with some initial value for the coefficient $`a_s(k_o)`$, we obtain a new value $`a_s(k_o)`$ by the relation

$$a_s^{}(k_o)=a_s(k_o)+\frac{\delta s}{u(k_o)}\frac{A}{a_s(k_o)}$$
(16)
with similar expressions for the other Fourier coefficients where $`\delta s`$ is a positive (negative) parameter if the action is a maximum (minimum). This procedure is then iterated until the partial derivatives of the action are sufficiently small to any desired accuracy. Our method differs in an essential way from the standard method of steepest descent by the important factor $`1/u(k)`$ which is essential for the convergence of the iterations with moderate values of $`\delta s`$ (order $`.1`$).

To obtain analytic approximations of Eqs. 8-11, it is convenient to determine the frequencies $`\omega `$ and $`\mathrm{\Omega }`$ by fixing the values of $`a_s(1)`$ and $`b_c(1)`$, which give the scale of the orbit and the magnitude of the total angular momentum in the rest frame, $`L_t=L+\mathrm{\Omega }I`$, where $`L`$ is the mean angular momentum in the rotating frame, Eq. 7,

$$L=3\omega \underset{k=1}{\overset{k=n}{}}k_oa_s(k_o)b_c(k_o)+k_ea_c(k_e)b_s(k_e),$$
(17)
and $`I`$ is the mean moment of inertia, Eq. 6,

$$I=\frac{3}{2}\underset{k=1}{\overset{k=n}{}}a_s(k_o)^2+a_c(k_e)^2+b_s(k_e)^2+b_c(k_o)^2.$$
(18)
We have

$$\omega =\frac{1}{2}(\sqrt{\alpha }+\sqrt{\beta })$$
(19)
and

$$\mathrm{\Omega }=\frac{1}{2}(\sqrt{\alpha }\sqrt{\beta }),$$
(20)
where

$$\alpha =\frac{f_s(1)g_c(1)}{(1b_c(1))}$$
(21)
and

$$\beta =\frac{f_s(1)+g_c(1)}{(1+b_c(1))}$$
(22)
are obtained from Eqs. 8 and 9 for $`k_o=1`$ and $`a_s(1)=1`$. These frequencies scale with the size $`s`$ of the orbit according to the generalized Kepler relations that $`\omega ^2s^{p+2}`$ and $`\mathrm{\Omega }^2s^{p+2}`$ are constants.

The generalization of Lagrange’s solution  is obtained by setting the Fourier coefficients $`a_c(k_e)`$ and $`b_s(k_e)`$ equal to zero, which according to Eqs. 11 and 10 implies that $`f_c(k_e)`$ and $`g_s(k_e)`$ also vanish. Small deviations from Lagrange’s circular solution can be obtained by setting $`a_s(1)=1`$ and $`b_c(1)=1\delta `$, and expanding Eqs. 8 and 9 in a power series $`\delta `$. In the limit of vanishing $`\delta `$

$$\alpha =\frac{1}{3^{p/4}}\sqrt{p/2}$$
(23)
which shows that these solutions exist only for $`p0`$. In particular, for the special case that $`p=2`$ corresponding to harmonic forces, $`\omega =\sqrt{3}`$, $`\mathrm{\Omega }=0`$, and all the Fourier coefficients $`a_s(k_0)`$ and $`b_c(k_0)`$ vanish for $`k_0>1`$.

## Numerical solutions and analytic approximations

Two numerical solutions of the extended Lagrange problem are shown in Figs. 1 and 2, where segments of the orbit occupied by each particle during a third of a period are indicated by a dark line, a dotted line and a dashed line respectively, and ten successive locations on these segments are marked at equal intervals of time to illustrate the relative position and variable speed of the three particles moving on this orbit. In this solution the angular momentum of the orbit decreases as $`b_c(1)`$ decreases, and the approximate elliptical orbit starts to get pinched in the middle. This is shown in Fig.1 for $`p=.05`$ at the special value $`b_c(1)=.049`$ where the orbit also develops cusps at the two ends. At these values of the parameters each particle comes momentarily to rest at these points because the attractive forces due to the other two particles balance the centrifugal force in the rotating frame.

Further decrease of the value of $`b_c(1)`$, corresponding to decreasing angular momentum, requires that the y-component of velocity at these ends change direction, which is possible only if an additional lobe appears at each end. A further decrease of $`b_c(1)`$ leads to an additional lobe at the center, as shown in Fig. 2 for $`p=.05`$ and $`b_c(1)=.0171`$.

The recently discovered  \-  zero angular momentum periodic orbit in the shape of a figure eight is a solution of Eqs. 8 and 10 with the Fourier coefficients $`a_c(k_e)`$ and $`b_c(k_o)`$ set equal to zero. According to Eqs. 11 and 9, this implies that $`f_c(k_e)=g_c(k_o)=0`$. For the special case that $`p=2`$, it can be shown that the moment of inertia, Eq. 6, is a constant. Consequently the shape of this orbit can be determined from purely geometrical considerations by imposing the condition that the Fourier coefficients of an expansion of this moment and that of the angular momentum vanish. Keeping terms in the sums in Eqs. 1 and 2 up to n=4, we obtain

$$b_s(4)=\frac{b}{d}\frac{b^3}{2d^4}(46+97b^2+49b^4),$$
(24)
$$a_s(5)=\frac{b^2}{d}\frac{b^2}{2d^4}(25+49b^2+22b^4),$$
(25)
$$a_s(7)=\frac{b^2}{2d^2},$$
(26)
$$b_s(8)=\frac{b^3}{2d^2},$$
(27)
where $`b=b_s(2)`$ and $`d=5+7b^2`$. The magnitude of $`b`$ is determined by setting $`\omega =\sqrt{f_s(1)}`$, Eq. 19, and solving Eq. 10 for $`k_e=2`$ which takes the form $`b=h(b)`$, where

$$h(b)=\frac{g_s(2)}{4f_s(1)}.$$
(28)
We find that $`h(b)`$ is a positive monotonically decreasing function for $`b>0`$ which diverges as $`b`$ approaches zero, and vanishes as $`b`$ increases indefinitely. Therefore, there exists only a single solution which we find at $`b=.3793`$ corresponding to $`\omega =1.448`$. This result is in good agreement with our numerical solution of Eqs. 8 and 10, which gives $`b=.3776`$ and $`\omega =1.452`$. The resulting figure eight orbit obtained from the Fourier coefficients given by Eqs. 24 \- 28 is shown in Fig. 3, which is graphically indistinguishable from the orbit obtained from the numerical solution of Eqs. 8 and 10. Surprisingly, these analytic coefficients also give a good approximation to the zero-angular-momentum figure eight orbit for $`p<2`$, because it turns out that the moment of inertia is approximately constant in this case . As $`p`$ decreases we find that $`b`$ decreases and vanishes when $`p=2`$, because in this limit of a harmonic force law the figure eight orbit becomes a linear collision orbit.

The zero-angular-momentum orbit provides the basis for a perturbation calculation of the corresponding orbit for small values of the angular momentum. Expanding the coefficients $`a_c(k_e)`$ and $`b_c(k_o)`$ in powers of $`b_c(1)=b`$, these coefficients together with the unperturbed coefficients $`a_s(k_o)`$ and $`b_s(k_e)`$ determine the angular momentum of the orbit. For example, for $`p=1`$ we find that $`\mathrm{\Omega }=1.647b`$, and substituting this value of $`\mathrm{\Omega }`$ and the lowest order Fourier coefficients in Eqs. 17 and 18 we obtain $`L_t=.5276b`$. Solutions for larger values of $`b`$ are obtained numerically by iterating these equations, fixing $`a_s(1)`$ and starting with some arbitrary value for $`b_s(2)`$, while all other Fourier coefficients are initially set equal to zero. As $`b`$ increases the figure eight orbit becomes more asymmetric with respect to the y-axis, as expected from the increasing action of the centrifugal and Coriolis forces in the rotating frame of reference, and at a critical value of the angular momentum the y-component of the velocity of the particle at the end of the orbit in a direction opposite to the rotation vanishes giving rise to a cusp. A further increase of the angular momentum then requires that this component of velocity change sign at this point, which give rise to an additional loop in the orbit, as shown in Fig. 4.

Finite angular momentum solutions can also be obtained by rotating the figure eight orbit about an axis in the plane of the original orbit leading to a three dimensional orbit as shown in Fig. 5

Solutions of our equations provide initial conditions for a numerical integration of the equations of motion which can then be applied to study the stability of these periodic orbits with respect to small perturbations. A preliminary investigation indicates that some of these solutions are stable over a restricted range of angular momentum and values of $`p`$ . Our method can readily be extended to investigate similar periodic motions for an arbitrary number of equal mass particles moving on the same orbit, many of which have been discovered recently . For example, any odd number $`n_o`$ of particles can travel on a figure eight orbit similar to that shown in Fig. 3. We find that this orbit, Fig. 6, approaches a unique limit as $`n_o\mathrm{}`$, but the frequency increases logarithmically with $`n_o`$.

In contrast, for the extension of Lagrange’s circular orbit each additional particle leads to a different orbital curve. Finally we have also applied our method to study the case of three particles with different masses which must move, however, along different orbital curves. An example is shown in Fig. 7 which is an extension of an orbit discovered by Moore  for the case of equal mass particles.

## Acknowledgements

I would like to thank Richard Montgomery for stimulating my interest in this problem and for many useful discussions, and Alan Chenciner and Charles Simó for helpul comments.


text:
# True and fake Lax pairs: how to distinguish them

## 1 Introduction

In this paper, we use the gauge-invariant description of zero-curvature representations of evolution equations, which was developed in , to solve the problem of how to distinguish the fake Lax pairs, introduced in , from the true Lax pairs. This approach is algorithmic and exploits only the $`x`$-part of a studied Lax pair. The main difference between the fake Lax pairs and the true ones turns out to be clearly seen from the structure of their *cyclic bases* (we use the terminology of  and explain it step by step in what follows).

Note that a more general gauge-invariant description, applicable to zero-curvature representations of any non-overdetermined systems of partial differential equations, was introduced in , in a very abstract form, and later the explicit expressions for the objects involved were provided in . The concept of cyclic bases of zero-curvature representations, however, did not appear in .

In Sections 2 and 3, respectively, we study the true and fake Lax pairs considered in . In Section 4, we discuss some further applications of the cyclic basis technique.

## 2 True Lax pairs

Let us study a true Lax pair first. We take the true Lax pair of the modified KdV equation $`u_t+u_{xxx}6u^2u_x=0`$ from :

$$\psi _{xx}=\left(u^2+\mathrm{i}k\frac{u_x}{u}k^2\right)\psi +\frac{u_x}{u}\psi _x,$$
(1)
$$\psi _t=\left(\mathrm{i}k\frac{u_{xx}}{u}2k^2\frac{u_x}{u}\right)\psi +\left(\frac{u_{xx}}{u}+2u^2+2\mathrm{i}k\frac{u_x}{u}+4k^2\right)\psi _x,$$
(2)
where $`k`$ is an essential (spectral) parameter. Introducing the vector $`\mathrm{\Phi }`$,

$$\mathrm{\Phi }=\left(\begin{array}{c}\varphi _1\\ \varphi _2\end{array}\right),\varphi _1=\psi ,\varphi _2=\psi _x,$$
(3)
we rewrite the Lax pair (1)–(2) as

$$\mathrm{\Phi }_x=X\mathrm{\Phi },\mathrm{\Phi }_t=T\mathrm{\Phi },$$
(4)
where

$$X=\left(\begin{array}{cc}0& 1\\ u^2+\mathrm{i}k\frac{u_x}{u}k^2& \frac{u_x}{u}\end{array}\right)$$
(5)
(the explicit form of the matrix $`T`$ will not be used in what follows). The compatibility condition

$$D_tX=D_xT[X,T]$$
(6)
of the over-determined linear system (4) is the *zero-curvature representation* of the modified KdV equation ($`D_t`$ and $`D_x`$ stand for the total derivatives, the square brackets denote the commutator).

For a generic matrix $`X(x,u,u_x,\mathrm{},u_{x\mathrm{}x})`$ and a generic evolution equation $`u_t=h(x,u,u_x,\mathrm{},u_{x\mathrm{}x})`$, we can rewrite the zero-curvature representation (6) in its *characteristic form* :

$$hC_u=H,$$
(7)
where $``$ is defined by $`M=D_xM[X,M]`$ for any matrix $`M`$, the matrix $`H`$ is determined by $`X`$ and $`T`$, and the matrix $`C_u`$ is the *characteristic* of $`X`$ with respect to $`u`$, defined as

$$C_u=\frac{X}{u}\left(\frac{X}{u_x}\right)+^2\left(\frac{X}{u_{xx}}\right)^3\left(\frac{X}{u_{xxx}}\right)+\mathrm{}.$$
(8)
Under the *gauge transformation*

$$X^{}=SXS^1+\left(D_xS\right)S^1,T^{}=STS^1+\left(D_tS\right)S^1,$$
(9)
generated by the transformation $`\mathrm{\Phi }^{}=S\mathrm{\Phi }`$ with any matrix $`S`$ ($`detS0`$), the corresponding transformation of (7) is tensor: $`C_u^{}=SC_uS^1`$, $`^{}=SS^1`$, $`H^{}=SHS^1`$. The *cyclic basis* is the sequence of linearly independent matrices $`C_u,C_u,\mathrm{},^{n1}C_u`$, where $`n`$ is maximal. In the *closure equation* of the cyclic basis,

$$^nC_u=a_0C_u+\mathrm{}+a_{n1}^{n1}C_u,$$
(10)
the coefficients $`a_0,\mathrm{},a_{n1}`$ (and, of course, the order $`n`$) are invariants with respect to the transformation (9).

Returning to the particular matrix $`X`$ (5), we compute $`C_u`$ and $`^iC_u`$, $`i=1,2,3`$, and find that $`C_u`$, $`C_u`$ and $`^2C_u`$ are linearly independent, whereas

$$^3C_u=4k^2\frac{u_x}{u}C_u+4\left(u^2k^2\right)C_u+\frac{u_x}{u}^2C_u.$$
(11)
We see that, in the case of the true Lax pair of the modified KdV equation, the cyclic basis is three-dimensional, with the closure equation given by (11). It is very important that some coefficients of the closure equation contain a parameter: from one hand, this means that the parameter $`k`$ cannot be removed (‘gauged out’) from $`X`$ by (9), since the coefficients of (11) are invariants, and, from the other hand, this has an essential influence on the structure of the class of all evolution equations which admit the linear problem (4) with this particular $`X`$, as we will see now. Let us find all those evolution equations $`u_t=h(x,u,u_x,\mathrm{},u_{x\mathrm{}x})`$ which are represented by (4), where $`X`$ is given by (5) and $`T(x,u,u_x,\mathrm{},u_{x\mathrm{}x})`$ is an arbitrary traceless $`2\times 2`$ matrix. Decomposing $`H`$ in (7) over the cyclic basis as $`H=pC_u+qC_u+r^2C_u`$ and using (11) for $`^3C_u`$, we find

$$q=\left(D_x+\frac{u_x}{u}\right)r,p=\left(D_x^2+D_x\frac{u_x}{u}+4\left(k^2u^2\right)\right)r$$
(12)
and

$$h=\left(A+4k^2B\right)r,$$
(13)
where

$$A=D_x^2u^1D_xu4D_xu^2,B=u^1D_xu.$$
(14)
Taking into account that $`r`$ can contain only finite-order derivatives of $`u`$ and that $`h`$ must be independent of $`k`$, we find from (13) that $`r`$ is a polynomial in $`k`$ and that the general expression for $`h`$ is $`h=R^n0`$, $`n=1,2,3,\mathrm{}`$, where

$$R=AB^1=D_x^24u^24u_xD_x^1u$$
(15)
is the recursion operator and $`D_x^10`$ must be interpreted as any constant. We have obtained a *discrete class* of evolution equations, namely, the integrable hierarchy of the modified KdV equation, whereas a zero-curvature representation without any parameter in the closure equation of its cyclic basis always leads to a *continual class* of evolution equations .

## 3 Fake Lax pairs

Now, let us study the fake Lax pair of the modified KdV equation, found in . Again, we need only the $`x`$-part of this Lax pair,

$$\psi _{xx}=\left(u^2+\mathrm{i}k\frac{u_x}{u}+k^2\right)\psi +\left(\frac{u_x}{u}2\mathrm{i}k\right)\psi _x,$$
(16)
which is remarkably similar to (1). Using (3), we rewrite (16) as $`\mathrm{\Phi }_x=X\mathrm{\Phi }`$, where

$$X=\left(\begin{array}{cc}0& 1\\ u^2+\mathrm{i}k\frac{u_x}{u}+k^2& \frac{u_x}{u}2\mathrm{i}k\end{array}\right).$$
(17)
For this $`X`$, we compute the characteristic $`C_u`$ (8), and then find that

$$C_u=0.$$
(18)
Thus, in this ‘fake’ case, the cyclic basis is one-dimensional, and the closure equation (18) contains no parameters. As a result of this, any evolution equation of the form $`u_t=D_xp(x,u,u_x,\mathrm{},u_{x\mathrm{}x})`$, where $`p`$ is arbitrary, admits a zero-curvature representation (6) with $`X`$ (17), as we can show by decomposing $`H`$ in (7) over the cyclic basis as $`H=pC_u`$. (Note that the class of represented equations may be even wider if we take into account the *singular basis* , but we will not consider this possibility here.) No discrete hierarchy and no recursion operator appeared in this case: we obtained a continual class of evolution equations, as should be expected from the structure of the cyclic basis.

In conclusion, let us study the most general fake Lax pair from . We rewrite its $`x`$-part (see (31a) in ) in the matrix form $`\mathrm{\Phi }_x=X\mathrm{\Phi }`$ with

$$X=\left(\begin{array}{cc}0& 1\\ \lambda f^2+\eta \mu f\eta ^2+\eta \frac{D_xf}{f}& \mu f2\eta +\frac{D_xf}{f}\end{array}\right),$$
(19)
where $`\lambda `$, $`\eta `$ and $`\mu `$ are arbitrary parameters and $`f(t,x,u,u_x,\mathrm{},u_{x\mathrm{}x})`$ is any function. Then we compute $`C_u`$ and $`C_u`$ for this $`X`$ (19), and find that

$$C_u=\frac{D_x\left(Ef\right)}{Ef}C_u,$$
(20)
where $`E`$ is the Euler operator, $`Ef=f/uD_x\left(f/u_x\right)+D_x^2\left(f/u_{xx}\right)\mathrm{}`$. We see that the cyclic basis for this fake Lax pair is one-dimensional again. The structure of the closure equation (20) suggests that the class of evolution equations, represented by (6) with $`X`$ (19), is continual: it contains at least all those evolution equations, for which $`f`$ is a conserved density. Indeed, considering $`\left(1\times 1\right)`$-dimensional matrices $`X`$ and $`T`$, for which $`[X,T]=0`$, and choosing $`X=f`$, we see that (6) is nothing but a conservation law with a conserved density $`f`$, whereas the closure equation for this $`X=f`$ is exactly the same as for $`X`$ (19), namely (20). Moreover, since the structure of the closure equation is invariant under the gauge transformations, the transformation (9) of $`X`$ (19) with

$$S=\mathrm{e}^{\eta x}\left(\begin{array}{cc}1& 0\\ \frac{\eta }{f}& \frac{1}{f}\end{array}\right),X^{}=\left(\begin{array}{cc}0& 1\\ \lambda & \mu \end{array}\right)f$$
(21)
reveals the origin of (20).

## 4 Discussion

We summarize that the difference between true and fake Lax pairs, considered in , is clearly seen from the structure of their cyclic bases. Below we give some comments on further applications of the cyclic basis technique.

Since the dimensions and the coefficients of cyclic bases of zero-curvature representations are gauge invariants, it is very convenient to use these invariants to verify whether a parameter in a given Lax pair is an essential (‘spectral’) parameter or this parameter can be removed (‘gauged out’) by a gauge transformation, and whether two given Lax pairs can (or cannot) be related to each other by a gauge transformation .

If a given Lax pair contains no essential parameter, there is always a continual class of evolution equations associated with the $`x`$-part of this linear problem, and most of those equations are non-integrable, of course . However, continual classes can also appear for linear problems with essential parameters, and this indicates that two dependent variables in such a spectral problem can be merged into one new dependent variable by a gauge transformation .

Last but not least, the cyclic basis technique is very useful to obtain hierarchies and recursion operators for given spectral problems .


text:
# Recent results from heavy ion collisions

## 1 Introduction

The purpose of the current heavy ion programs at CERN (Switzerland) and Brookhaven National Laboratory (USA) is to probe strongly interacting matter under extreme conditions, i.e. at high densities and temperatures. The central subject of these studies is the transition from the quark-gluon plasma to hadronic matter. In the early phases of ultra-relativistic heavy ion collisions, when a hot and dense region is formed in the center of the reaction, there is copious production of up, down, and strange quarks. Transverse expansion is driven by multiple scattering among the incoming and produced particles. As the medium expands and cools, the quarks/gluons combine to form the hadrons that are eventually observed.

In this paper, we will summarize the recent experimental results on transverse momentum distributions and particle ratios. The related physics issues are chemical equilibrium and collective expansion. While the former has to do with inelastic collisions, the later is dominated by the elastic cross section. We will focus at the relatively low transverse momentum region, $`p_t`$ $`2`$ GeV/c, where the bulk production occurs. Experimental results are from refs. \[1-9\] and for results on high momentum transfer, global measurements, and correlations readers are referred to ref. .

## 2 Baryonic Physics: mid-rapidity $`\overline{p}/p`$ ratios and inclusive net-proton distributions

All RHIC experiments have consistently measured the ratio of $`\overline{p}/p`$ at mid-rapidity. It is found that at mid-rapidity, the ratio is not unity and a dramatic increase in the ratio is observed from SPS to RHIC, meaning that the system created at $`\sqrt{s_{_{\mathrm{NN}}}}`$ = 130 GeV is not yet net-baryon free.

One of the unsolved problems in high energy nuclear collisions is the baryon transfer that occurs at the early stage of the collision. The later dynamic evolution of the system is largely determined at this moment since the available energy is fixed at this stage of the collision. On the theoretical side, a novel mechanism, the non-perturbative gluon junction, was proposed to address this problem. From the STAR measured $`\overline{p}/p`$ ratios and the STAR preliminary results of anti-proton yields, we have extracted the net-proton yields at mid-rapidity. The values are shown in Fig. 1 as a function of collision centrality. No hyperon decay corrections have been applied as that requires the knowledge of the yields of hyperons.

One can see that the number of net-protons increases as the number of the negatively charged particles increases. This implies that more and more protons (baryons) are transfered from the incoming nuclei to mid-rapidity as the overlap of the two nuclei increases. More data are needed to study the role of the junction mechanism in heavy ion collisions. In this respect, the rapidity distributions of the net-protons as a function of collision centrality will be of crucial importance.

## 3 Chemical Freeze-out: particle ratios

Recently, much theoretical effort has been devoted to the analysis of particle production within the framework of statistical models. These approaches are applied to the results of both elementary collisions ($`e+e^{}`$, $`p+p`$) and heavy ion collisions ($`Au+Au`$ and $`Pb+Pb`$). Many features of the data imply that a large degree of chemical equilibration may be reached both at AGS and SPS energies. The three most important results are: (i) at high energy collisions the chemical freeze-out (inelastic collisions cease) occurs at about 160-180 MeV and it is ‘universal’ to both elementary and heavy ion collisions; (ii) the kinetic freeze-out (elastic scatterings cease) occurs at a lower temperature $`120140`$ MeV; (iii) the compilation of freeze-out parameters in heavy ion collisions in the energy range from 1 - 200 A$``$GeV shows that a constant energy per particle $`E/N1`$ GeV can reproduce the behavior in the temperature-potential ($`T_{ch}\mu __B`$) plane.

The systematics for chemical parameters is shown in the phase plot in Fig. 2 where the dashed lines represent the boundary between interactions involving hadronic and partonic degrees of freedom. The dotted lines represent the results of ref. . The lattice prediction on $`T_c`$ ($`160`$ MeV) is shown at $`\mu _B=0`$ and the baryon density for neutron stars is indicated at $`T_{ch}=0`$.

## 4 Kinetic Freeze-out: transverse momentum distributions

The measured transverse momentum distributions have been fitted by the exponential function $`f=Aexp(m_t/T),`$ where $`T`$ is the slope parameter and $`A`$ is the normalization constant. The magnitude of the slope parameter provides information on temperature (random motion in local rest frame) and collective transverse flow. Fig. 3(a) shows the measured particle slope parameters from Pb+Pb central collisions at SPS ($`\sqrt{s_{_{\mathrm{NN}}}}`$ = 8.8 GeV and $`\sqrt{s_{_{\mathrm{NN}}}}`$ = 17.2 GeV) energies. Recent results are:

* Particle distributions from 40 A$``$ GeV ($`\sqrt{s_{_{\mathrm{NN}}}}`$ = 8.8 GeV) collisions were reported by NA45 and NA49. The slope parameters are similar to the ones observed at 158A$``$ GeV, i.e., they follow the established systematic trend;
* The systematic of the transverse momentum distributions for $`J/\psi `$ was reported by the NA50 experiment. It is interesting to note that the slope parameter of $`J/\psi `$ is similar to those for $`\varphi `$ and $`\mathrm{\Omega }`$;
* For central collisions the $`\varphi `$ slope parameter of NA49 is about 300 MeV whereas that of the NA50 is about 240 MeV. NA49 and NA50 reconstructed $`\varphi `$ mesons via $`K^+K^{}`$ and $`\mu ^+\mu ^{}`$ channels, respectively. As discussed in refs. , part of the difference may be caused by the final state interaction of the decay kaons. In fact, if one studies the centrality dependence of the $`\varphi `$ slope parameter, one finds that at peripheral collisions, both experimental results agree with each other and the value is close to that from $`p+p`$ collisions.

Preliminary results of the slope parameters from Au+Au central collisions at RHIC ($`\sqrt{s_{_{\mathrm{NN}}}}`$ = 130 GeV) are shown in Fig. 3(b). It is obvious that the mass dependence of the slope parameter is stronger than that from collisions at SPS energies.

As one can see from Fig. 3(a), the slope parameters appear to fall into two groups: group (I) is flat as a function of particle mass, whereas the slopes in group (II) increase strongly with particle mass. At the RHIC energy, the slope parameter systematic of $`\pi ,K,`$ and $`p`$ shows an even stronger dependence on the particle mass. The strong energy dependence of the slope parameter might be the result of the larger pressure gradient at the RHIC energy. With a set of reasonable initial/freeze-out conditions and equation of state, the stronger transverse expansion at the RHIC energy was, in deed, predicted by hydrodynamic calculations. In addition, the results are consistent with the large value of the event anisotropy parameter $`v_2`$ at the RHIC energy.

Within the hadronic gas, the interaction cross section for particles like $`\varphi ,\mathrm{\Omega },`$ and $`J/\psi `$ are smaller than that of $`\pi ,K,`$ and $`p`$. Therefore the interactions between them and the rest of the system are weak, leading to the flat band behavior in Fig. 3(a). On the other hand, the slope parameter of these weakly interacting particles may reflect some characteristics of the system at hadronization. Then it should be sensitive to the strength of the color field. Under this assumption, the fact that the weakly interacting particles show a flat slope parameter as a function of their mass would indicate that the $`flow`$ develops at a later stage of the collision. Should the collective flow develop at the partonic level, one would expect to see a mass dependence of the slope parameters for all particles.

Figure 4 shows the bombarding energy dependence of the freeze-out temperature $`T_{fo}`$ and the average collective velocity $`\beta _t`$. These parameters were extracted by hydrodynamic motivated model fitting. It is interesting to observe that both $`T_{fo}`$ and $`\beta _t`$ saturate at a beam energy of about 10 A$``$GeV. The saturation temperature is about 100 - 120 MeV, very close to the mass of the lightest meson. The steep rise of the $`T_{fo}`$ and $`\beta _t`$ up to about 5 A$``$GeV incident energy indicates that at low energy collisions the thermal energy essentially goes into kinetic degrees of freedom. The saturation at $``$10 A$``$GeV shows that particle generation becomes important. As proposed in, for a pure hadronic scenario there may be a limiting temperature $`T_c140`$ MeV in high-energy collisions, although the underlying physics for both, the transition from partonic to hadronic degrees of freedom and the transition from interacting hadrons to free-streaming is not clear at the moment. By coupling the limiting temperature idea to a hydrodynamic model calculation, Stöcker et al. successfully predicted the energy dependence of the freeze-out temperature.

At the RHIC energy, the collective velocity parameter seems to be larger than that from collisions at AGS/SPS energies. This can already be seen in Fig. 4 where the increase from pion to proton at $`\sqrt{s_{_{\mathrm{NN}}}}`$ = 130 GeV (Fig. 3(b)) is much faster than at $`\sqrt{s_{_{\mathrm{NN}}}}`$ = 17 GeV (Fig. 3(a)). Is this a manifestation of van Hove’s picture? On the other hand, compared to results from lower energy collisions, the temperature parameters seem to be lower. Is this the consequence of hydrodynamic expansion from a higher initial density fireball? An energy scan between $`\sqrt{s_{_{\mathrm{NN}}}}`$ = 20 - 130 GeV will be extremely important in order to study this evolution in more detail.

## 5 Summary

In summary, the most interesting results from collisions at RHIC are that the system is indeed approaching net-baryon free and the transverse expansion is found to be stronger than that from collisions at AGS/SPS energies.

In order to understand the trend of the collective velocity, an energy scan between $`\sqrt{s_{_{\mathrm{NN}}}}`$ = 20 - 200 GeV, is important. In addition, systematic studies on the anisotropy parameter $`v_2`$ and the transverse momentum distributions of $`\varphi ,\mathrm{\Omega },`$ and $`J/\psi `$ are necessary as they will help in determining whether the collectivity is developed at the partonic stage.

—–

We are grateful for many enlightening discussions with Drs. P. Braun-Munzinger, M. Gyulassy, H.G. Ritter, K. Schweda, E.V. Shuryak, F.Q. Wang, and X.N. Wang. This work has been supported by the U.S. Department of Energy under Contract No. DE-AC03-76SF00098.


text:
# 1 Introduction

## 1 Introduction

Inflation probably took place in the early universe. If the scale of (the latest) inflation was of the order of the electroweak scale, the subsequent reheating would provide ideal conditions for electroweak baryogenesis.

The original scenario of Kuzmin, Rubakov, and Shaposhnikov  was based on a brilliant idea that all three Sakharov’s conditions  necessary for a successful baryogenesis were satisfied, at least qualitatively, in the early universe due to the properties of the Standard Model at finite temperature. First, the baryon number was violated by sphaleron transitions. Second, the universe was out of thermal equilibrium during a phase transition. Finally, CP was broken, because it is not an exact symmetry in the Standard Model.

Unfortunately, this scenario, in its original form, cannot account for the observed baryon-to-photon ratio $`\eta 10^{10}`$. First problem has to do with the fact that the electroweak phase transition is not sufficiently strongly first-order unless the Higgs mass is below about 45 GeV, which is ruled out by experiment. Second, the CP violation from the Cabibbo-Kobayashi-Maskawa matrix is way too small for baryogenesis because its contribution to $`\eta `$ is suppressed by high powers of the Yukawa couplings.

At the end of inflation, however, the universe was emphatically out of equilibrium. Could the baryon number violating processes take place in that non-thermal environment? The answer is yes; this was clearly demonstrated by a consensus of numerical  and analytical  arguments. Finally, the requisite CP violation may come from time-dependent solutions for the scalar zero-modes in the background . Altogether, electroweak baryogenesis can be very efficient at the end of inflation.

## 2 Electroweak baryogenesis at preheating

The idea  of using for baryogenesis the non-equilibrium nature of preheating , along with the new sources of CP violation that preheating offers , is very appealing because it is simple and requires only a minor “upgrade” of the Standard Model. The inflaton has to be added to the Standard Model, but this is a widely accepted piece of new physics that is desirable for other reasons. The only uncomfortable features of the scenario proposed independently in Refs.  and  are the relatively low scale of inflation and the need to couple the inflaton to the Higgs sector. In combination with the constraint that inflation must last for more than 60 e-folds, this forces the inflaton potential to be unnaturally flat.

Of course, if the universe has undergone more than one period of inflation, the latest, electroweak-scale inflation could have lasted for only a small number of e-folds. In this case, there is no problem with naturalness.

Low-scale inflation models have been constructed . One of the side benefits of lowering the inflation scale is avoiding the gravitino over-production constraints. There are other particle-physics motivations for using the TeV scale, which is associated with supersymmetry breaking in a class of models .

Although the simplest models of hybrid inflation do not yield parametric resonance , there is no fundamental reason why a weak-scale inflation could not be followed by resonant preheating . The latter possibility allows for efficient electroweak baryogenesis at preheating.

### 2.1 Baryon number non-conservation during preheating

Numerical simulations in 1+1 dimensions  demonstrate that the baryon number violating processes can go unsuppressed in a non-thermal plasma at preheating. One way to think about this is in terms of “effective temperature” . Although plasma is far from equilibrium, the Chern-Simons number changes at the same rate as in thermal plasma at temperature $`T_{\mathrm{eff}}`$. Since sphalerons, being extended objects, are mainly affected by the long-wavelength modes, it is the infra-red part of the spectrum that determines $`T_{\mathrm{eff}}`$. At preheating, the long-wavelength modes are overpopulated as compared to the thermal ensemble. Hence, the effective temperature $`T_{\mathrm{eff}}`$ is higher than the energy density. When thermal equilibrium is restored, the final temperature $`T_R<100\mathrm{GeV}<T_{\mathrm{eff}}`$. At that point the rate of baryon number violating processes is negligible. This behavior, favorable to baryogenesis, is shown in Figs. 1,2.

One can try to describe this behavior semi-analytically by isolating certain gauge degrees of freedom and studying their evolution at preheating , when the vector boson mass, $`M_W`$ changes with time due to the oscillations of the inflaton coupled to the Higgs sector. Let us introduce an ansatz  for the gauge potential

$$gA_\mu =(\frac{\tau _a}{2i})A_\mu ^a.$$
(1)
Our spatially-homogeneous ansatz is:

$$gA_0=0;gA_i=(\frac{\tau _i}{2i})\varphi (t)$$
(2)
in which the group index is tied to the spatial index. By the conventional rules of charge conjugation and parity for the gauge potential, $`\varphi `$ is C even, P odd, CP odd.

It is important to note that this ansatz does not correspond to a non-vanishing VEV for an EW field. Gauge invariance alone is enough to ensure that there can be no expectation value coupling the space-time indices to group indices.

One readily calculates the EW electric and magnetic fields:

$$gE_iG_{0i}=(\frac{\tau _i}{2i})\dot{\varphi }(t);gB_i\frac{1}{2}ϵ_{ijk}G_{jk}=(\frac{\tau _i}{2i})\varphi ^2.$$
(3)
Then one calculates the density $`W`$ of the Chern-Simons number as:

$$W=(\frac{1}{8\pi ^2})\varphi ^3.$$
(4)
It is straightforward to check that $`\dot{W}`$ is the topological charge density $`Q`$, related to B+L violation through the anomaly equation.

With the assumption of a given Higgs VEV, the equations of motion for the gauge potential are:

$$[D^\mu ,G_{\mu \nu }]+M_W^2(t)(A_\nu +(_\nu U)U^1)=0.$$
(5)
Here the unitary matrix $`U`$ represents the Goldstone (phase) part of the Higgs field. The mass term will be assumed to have the form:

$$M_W^2(t)=m^2(1+ϵ\mathrm{cos}(\omega t))$$
(6)
where $`m`$ is the value of $`M_W`$ with no oscillations.

This ansatz has some important properties. First, it carries some non-zero Chern-Simons density. Second, the equation of motion for the field $`\varphi `$,

$$\ddot{\varphi }+2\varphi ^3+(1+ϵ\mathrm{cos}rt)\varphi =0,$$
(7)
has resonantly growing solutions. Finally, the ansatz can evolve for a long time without breaking spatial homogeneity .

The Chern-Simons density grows as the third power of $`\varphi `$.

### 2.2 CP violation at preheating

CP and time-reversal symmetries at preheating may be broken by the time-dependent classical motions of scalar fields in the background . This is similar to the source of CP violation used by Affleck and Dine . CP violation of this kind is poorly constrained by experiment. Several viable scenarios for electroweak baryogenesis at preheating were presented in Ref. . For example, a modified spontaneous baryogenesis a la Cohen, Kaplan, and Nelson  becomes very efficient at preheating. Their original scenario used the variation of the Higgs field inside a wall of a bubble formed in a first-order phase transition. A similar effect can occur at preheating uniformly in space, on the horizon scales . One can get the desired baryon asymmetry in a Standard Model supplemented by an additional Higgs doublet and an inflaton sector . The difference with the scenario proposed by Cohen, Kaplan and Nelson  is that in our case CP violation occurs homogeneously in space, as opposed to in a bubble wall. In addition, the final prediction for baryon asymmetry in the CKN scenario was very far from the equilibrium value because the sphaleron rate was slow on the time scales associated with the growth of the bubble. In our case, the Higgs parameters change slowly in time while the baryon number non-conservation is rapid. This allows a slow adiabatic adjustment of the baryon number to that which minimizes the free energy.

## 3 Other possibilities for baryogenesis at preheating

Non-thermal production of heavy particles at the end of inflation  could result in baryogenesis , or a leptogenesis  via the decay of the heavy particles.

Alternatively, in models with low-energy supersymmetry, weak-scale inflation can lead to formation of the Affleck-Dine condensate . In general, this condensate does not remain homogeneous and can break up  into SUSY Q-balls . This affords a number of interesting possibilities for generating baryons and dark matter simultaneously .

If there were complex scalar fields not associated with supersymmetry, their fragmentation in a similar fashion could produce stable Q-balls that could make up the dark matter. Interestingly enough, this form of dark matter can be strongly self-interacting  and can, therefore, resolve the problems with non-interacting cold dark matter on small scales .

## 4 Conclusion

A highly non-thermal state of the universe at the end of inflation provides a fertile ground for baryogenesis. In addition to a dramatic departure from thermal equilibrium, preheating creates conditions for baryon number non-conservation and activates sources of CP violation that are largely unconstrained by the experimental data. Several modest modifications of the Standard Model allow for efficient electroweak baryogenesis in the wake of inflation.


text:
# Vortices in atomic-molecular Bose-Einstein condensates

## I Introduction

Recent experimental observation of the photoassociation of Bose condensed atoms via the stimulated Raman process and the production of ultracold molecules  is, undoubtedly, the first step towards the creation of a coherently coupled atomic-molecular condensate. Such a system would represent a rich ground for exploring a variety of the nonlinear atom-optics effects, analogous to those known for the processes of parametric up and down conversion of light waves . From this point of view, the excited topological states of the macroscopic wavefunction– vortices– in the atomic-molecular Bose-Einstein condensates are of great interest as a direct counterpart of the parametric vortices in nonlinear optics . On the other hand, the study of such topological states is important for understanding the properties of composite superfluids, such as the recently proposed molecular Bose-Einstein condensate (BEC) coupled to a cloud of Fermi atoms , or the composite atomic-molecular BEC (AMBEC) .

From the experimental point of view, it is important to not only establish the conditions for the formation of topological states, but also identify the states which are dynamically stable, i.e. are not destroyed by exponentially growing collective excitations of the condensate. One of the proposed methods for the observation of vortices involves imaging of the expanding condensate released from a trap, with a vortex core growing during the expansion. The atomic-molecular condensate represents an interesting experimental challenge as the excited states may be formed in a dynamically stable, self-confined droplet , after the hybrid condensate is released from a trap.

Assuming further experimental progress towards the creation of the molecular and hybrid atomic-molecular BECs, in this paper we analyse the dynamics of the parametrically coupled atomic and molecular condensates, in the framework of a two-component mean-field model  that takes into account all types of the mean-field atomic and molecular interactions . For the first time to our knowledge, we study the structure, stability, and dynamics of topological states in AMBEC, and, in particular, we reveal that the dynamically stable two-component vortex solitons can be formed in both trapped and untrapped AMBEC clouds. Unlike the conventional vortices predicted and already observed in atomic BEC, vortex solitons can exist as stable self-trapped states even without a trapping potential, supported by an effectively attractive interaction between the parametrically coupled atomic and molecular components.

Another crucial problem is the dynamics of AMBEC and its topological states in the presence of dissipation, i.e. loss of atoms and molecules from the corresponding fractions of the hybrid BEC. The formation of BEC in the presence of losses due to either a coherent Raman photoassociation process or Feshbach resonance has been analyzed previously , however the dissipative dynamics of vortices, earlier studied for a single-component atomic BEC , has never been discussed in the context of hybrid condensates.

Since, so far, no conclusive experimental data are available that would allow us to estimate the rates of the condensate losses rigorously, in this work we take them into account by including the phenomenological dissipative terms, and then considering the dissipation rates for which the topological states can still be formed. We analyze the dynamics of the vortex solitons in the presence of dissipation and reveal that they can exhibit a nontrivial decay scenario, with the decay rates much smaller than the characteristic loss rates. Since the physical principles underlying the existence and stability of these topological states are rather generic, we expect that the similar structures can be found in more realistic models of the atomic-molecular BEC dynamics, and eventually observed in experiment.

## II Model

We study the dynamics of the hybrid atom-molecular condenstates, produced by a coherent stimulated Raman photoassociation process in a strongly anisotropic trap with a tight confinement along one spatial dimension. We model the parametrically interacting condensates by a system of two coupled Gross-Pitaevskii (GP) equations for the macroscopic wave functions of the atomic ($`\mathrm{\Psi }_a`$) and molecular ($`\mathrm{\Psi }_m`$) components (see details and references, e.g., in Ref. ):

$`i{\displaystyle \frac{\mathrm{\Psi }_a}{t}}+{\displaystyle \frac{1}{2}}\mathrm{\Delta }\mathrm{\Psi }_aV_a\mathrm{\Psi }_a\chi \mathrm{\Psi }_a^{}\mathrm{\Psi }_me^{i\delta t}`$ (1)
$`(U_{aa}|\mathrm{\Psi }_a|^2+U_{am}|\mathrm{\Psi }_m|^2)\mathrm{\Psi }_a=iR_a,`$ (2)
$`i{\displaystyle \frac{\mathrm{\Psi }_m}{t}}+{\displaystyle \frac{1}{4}}\mathrm{\Delta }\mathrm{\Psi }_mV_m\mathrm{\Psi }_m{\displaystyle \frac{1}{2}}\chi \mathrm{\Psi }_a^2e^{i\delta t}`$ (3)
$`(U_{mm}|\mathrm{\Psi }_m|^2+U_{am}|\mathrm{\Psi }_a|^2)\mathrm{\Psi }_m=iR_m,`$ (4)
where the wave functions, time, and spatial coordinates are measured in the units of $`x_0^{3/2}`$, $`\omega _{}^1`$, and $`x_0(\mathrm{}/m\omega _{})^{1/2}`$, respectively. Here $`\omega _{}`$ is the characteristic frequency in the plane $`(x,y)`$ of the weak confinement of the trap for the atomic BEC fraction. Consequently, $`V_a=(1/2)r^2+(\lambda ^2/2)z^2`$ and $`V_m=r^2+\lambda ^2z^2`$, with $`r^2=x^2+y^2`$, and $`\lambda ^2=\omega _z/\omega _{}1`$. In these notations, the interaction strengths $`U_{ij}`$ are measured in the units of $`(\mathrm{}\omega _{})^1x_0^3`$, and the Raman-induced coupling, $`\chi `$, in the units of $`(\mathrm{}\omega _{})^1x_0^{3/2}`$. The terms $`R_{a,m}`$ on the right hand side are the phenomenological loss terms which will be discussed and specified below in Sec. 4.2.

Due to significantly dissimilar scales of spatial confinement in the different directions, and the quasi-two-dimensonal geometry of the trap, the form of the condensate wavefunction in the $`z`$ direction is completely determined by the structure of the harmonic trap. We therefore assume that the wavefunction can be factorized in the following way: $`\mathrm{\Psi }_{a,m}=\varphi _{a,m}(x,y,t)\mathrm{\Phi }_{a,m}(z)`$, where $`\mathrm{\Phi }_{a,m}(z)`$ satisfy the one-dimentional harmonic oscillator equations:

$$\frac{1}{2}\frac{^2\mathrm{\Phi }_{a,m}}{z^2}v_{a,m}(z)\mathrm{\Phi }_{a,m}+ϵ_{a,m}\mathrm{\Phi }_a=0,$$
(5)
where $`v_a=(1/2)\lambda ^2z^2`$ and $`v_m=2\lambda ^2z^2`$. The ground state solutions of these equations correspond to $`ϵ_a=\lambda /2`$ and $`ϵ_m=\lambda `$, so that $`\mathrm{\Phi }_a=C_a\mathrm{exp}(\lambda z/2)`$ and $`\mathrm{\Phi }_m=C_m\mathrm{exp}(\lambda z)`$. The normalization conditions $`\mathrm{\Phi }_{a,m}^2=1`$ yield $`C_a=(\lambda /\pi )^{1/4}`$ and $`C_m=(2\lambda /\pi )^{1/4}`$.

Substituting the factorized wavefunctions into the model equations, multiplying the first and second equations by $`\mathrm{\Phi }_a`$ and $`\mathrm{\Phi }_m`$, respectively, integrating to eliminate the $`z`$ dependence, and introducing $`\psi _{a,m}(x,y,t)=\varphi _{a,m}(x,y,t)\mathrm{exp}(i\lambda t/2)`$, we obtain the final evolution equations for the wavefunctions of a “pancake”-shape AMBEC:

$`i{\displaystyle \frac{\psi _a}{t}}+{\displaystyle \frac{1}{2}}\mathrm{\Delta }_{}\psi _aV_a(r)\psi _a\stackrel{~}{\chi }\psi _a^{}\psi _me^{i\delta t}`$ (6)
$`\left(\stackrel{~}{U}_{aa}|\psi _a|^2+\stackrel{~}{U}_{am}|\psi _m|^2\right)\psi _a=i\stackrel{~}{R}_a,`$ (7)
$`i{\displaystyle \frac{\psi _m}{t}}+{\displaystyle \frac{1}{4}}\mathrm{\Delta }_{}\psi _mV_m(r)\psi _m{\displaystyle \frac{1}{2}}\stackrel{~}{\chi }\psi _a^2e^{i\delta t}`$ (8)
$`(\stackrel{~}{U}_{mm}|\psi _m|^2+\stackrel{~}{U}_{am}|\psi _a|^2)\psi _m=i\stackrel{~}{R}_m,`$ (9)
where $`\mathrm{\Delta }_{}=^2/x^2+^2/y^2`$, $`V_a(r)=(1/2)r^2`$, and $`V_m(r)=r^2`$. Here we introduce the renormalized interaction strengths: $`\stackrel{~}{U}_{aa}=U_{aa}(\lambda /2\pi )^{1/2}`$, $`\stackrel{~}{U}_{mm}=U_{mm}(\lambda /\pi )^{1/2}`$, $`\stackrel{~}{U}_{am}=U_{am}(2\lambda /3\pi )^{1/2}`$, $`\stackrel{~}{\chi }=\chi (\lambda /2\pi )^{1/4}`$ and, correspondingly, the renormalized loss terms $`\stackrel{~}{R}_{a,m}`$. For the sake of clarity, we omit the tilde from the equations throughout the following text. One can see that the renormalization factors absorbed by the interaction coefficients as a result of the reduction to the 2D geometry are, in fact, close to $`1`$ for typical $`\lambda 10`$.

For definiteness, we consider the case of $`{}_{}{}^{87}Rb`$, for which the production of ultracold molecules were recently demonstrated experimentally , and assume that the Raman detuning parameter vanishes, $`\delta =0`$. Using the recently reported data for the intra- and inter-species scattering lengths , we take $`U_{aa}=0.062`$; $`U_{am}=0.084(\pm 0.07)`$, and $`\chi =1.09`$. The large uncertainty in the value of $`U_{am}`$ corresponds to the s-wave scattering length, $`a_{am}=180a_0\pm 150a_0`$, recently measured experimentally . Since the molecular BEC has not been experimentally realized yet, the exact value of $`U_{mm}`$ is not known, although it can be estimated to be of the order of the other interaction strengths (see details in Ref. ) and, therefore, it cannot be ignored in the model (6).

It has also been previously established that the qualitative behaviour of the system depends on the sign of the determinant $`\mathrm{\Delta }UU_{aa}U_{mm}U_{am}^2`$ , which defines the relative strength of the cubic interactions. Parameter $`U_{mm}`$ can therefore be used to investigate the system in two distinct regimes of the net attractive $`\mathrm{\Delta }U>0`$ and net repulsive $`\mathrm{\Delta }U<0`$ cubic interactions . At the low particle numbers, i.e. for the low densities of the condensate wavefunctions, and in the parameter region where $`\mathrm{\Delta }U0`$, the interactions are dominated by the parametric, i.e. quadratic, coupling ($`\chi `$) which can be effectively attractive or repulsive, depending on the relative phase of the molecular and atomic components . The structure and dynamical properties of all topological structures in the AMBEC is controlled by the complex interplay of the quadratic and cubic inter- and intra-component mean-field interactions.

## III Vortex solitons

Spatially localized ground-state solutions of the model (6) at $`R_a=R_m=0`$, with and without the trapping potential $`V_j(r)`$, have been analyzed earlier . Here, we are interested in the excited states of the model (6) that correspond to the topological vortex states of the hybrid two-component condensate. In order to find the vortex states, we consider the model equations (6) without the dissipative terms (i.e. $`R_a=R_m=0`$), and look for the radially symmetric solutions with a nontrivial phase in the form:

$`\psi _a(r,\theta ,t)=\varphi _a(r)e^{ik_a\theta i\mu _at},`$ (10)
$`\psi _m(r,\theta ,t)=\varphi _m(r)e^{ik_m\theta i\mu _mt},`$ (11)
where $`(r,\theta )`$ are the polar coordinates in the plane of the trap, and $`\mu _{a,m}`$ are the chemical potentials of the corresponding AMBEC fractions. Due to the parametric coupling between the atomic and molecular matter waves, the existence of stationary states of the AMBEC requires the phase-matching conditions, $`\mu _m=2\mu _a2\mu `$, and $`k_m=2k_a2n`$, so that the integer $`n`$ defines the topological charge of the AMBEC vortex soliton via its atomic component.

The spatially localized, radially symmetric solutions of the model (6) in the form (10) can be found by a standard numerical relaxation technique. Different types of such topological states are summarized in Fig. 1 for the case $`n=1`$ (single-charged vortices). The upper solid curve in the middle part of Fig. 1 \[that includes the point (d)\] corresponds to stable vortices in a trapping potential; such vortices are supported by the effectively repulsive cubic interaction since $`\mathrm{\Delta }U<0`$. These topological states are somewhat similar to the vortices of a single-component atomic BEC, and they expand when released from the trap. At the low particle numbers \[see the inset in Fig. 1, dashed curve with the point (c)\], these vortex solutions become unstable due to the competition between the effectively attractive quadratic and repulsive cubic interaction. Thus, the repulsive nature of the inter-species cubic interaction manifests itself only for large densities of the components $`\varphi _a`$ and $`\varphi _m`$.

Thick solid curve in Fig. 1 (and the insert) describes a family of the vortices of the single-component purely molecular condensate wavefunction $`\psi _m`$ \[see the case (a)\], whereas a dotted curve above it corresponds to another family of hybrid vortices that bifurcates from the purely molecular states. This second family of the two-component vortices is supported by an effectively repulsive quadratic nonlinearity allowing such states to exist for low particle numbers but only in the presence of a trapping potential \[cf. an example in (b)\].

An important feature of the model (6) is the existence of untrapped self-localized vortex solitons for $`\mu <0`$. These structures are supported by an effectively attractive quadratic interaction between the atomic and molecular components. As the parameter $`\mu `$ decreases, the solution becomes broader due to the presence of effectively repulsive cubic nonlinearity, with a saturation amplitude of the wavefunction corresponding to the amplitude of the spatially homogeneous (free-space) condensate wavefunction \[see the cases (e) and (f) in Fig. 1\].

In this paper, we determine the vortex stability by a direct numerical propagation, but a more rigorous stability analysis will be presented elsewhere. We find that the broad vortex solitons \[shown, e.g., in Fig. 1(f)\] are stable, whereas narrow vortex solitons \[e.g. those shown in Fig. 1(e)\] are unstable, in agreement with the earlier analysis of the parametric vortices in optics . Interestingly, for the unstable vortices we observe two types of the instability-induced dynamics, i.e. the break-up into two or three filaments, for larger and smaller particle numbers, respectively. This dynamical behaviour is a clear signature of the transition between the regime of effectively cubic interaction in which the vortex soliton decays into $`2|n|`$ filaments , and the regime of effectively quadratic interaction when the vortex soliton breaks up into $`2|n|+1`$ filaments, as usually observed for quadratic parametric beams .

## IV Dissipation-induced vortex dynamics

### A Physical estimates of the losses

In the formation of ultracold molecules from an atomic BEC, the processes of the incoherent production of molecules, decays of the atomic excited state, and inelastic collisions between atoms and molecules, all contribute to the the effective loss of particles from the corresponding fractions of the BEC, resulting in the dissipative dynamics of the localized AMBEC states.

To take into account all possible losses in the hybrid condensate, we introduce the following linear and nonlinear (i.e., density-dependent) dissipative terms into the GP model (6),

$$\begin{array}{c}R_a=\nu _a\psi _a\gamma _a|\psi _a|^2\psi _a\gamma _{am}|\psi _m|^2\psi _a,\hfill \\ R_m=\nu _m\psi _m\gamma _m|\psi _m|^2\psi _m\gamma _{am}|\psi _a|^2\psi _m,\hfill \end{array}$$
(12)
where the dimensionless coefficients, $`\nu _j`$ and $`\gamma _j`$, stand for the decay rates of the atomic excited state ($`\nu _a`$) and the excited molecules into high bound states of the atomic BEC or “hot” (non-BEC) atoms ($`\nu _m`$), the rates of inelastic collisions and decay of the atomic states ($`\gamma _a`$), and inelastic molecule-molecule collisions and dissociation into high bound state or continuum state of the molecular BEC ($`\gamma _m`$), and atom-molecule inelastic collisions resulting in the production of “hot” atoms and molecules ($`\gamma _{am}`$).

The linear loss terms in Eq. (12) are similar to those discussed earlier in Ref. , however the latter work did not take into account the losses due to the rotationally or vibrationally inelastic atom-molecule collisions, as well as the molecule-molecule collision losses. System (6) with the loss terms (12) is formally equivalent to the mathematical description of the formation of AMBEC through Feshbach resonances , except for the terms responsible for the linear losses from the excited atomic state ($`\nu _i`$) which are specific for the AMBEC production via the coherent Raman photoassociation process. We stress that no conclusive experimental data for the losses that would enable the coherent AMBEC production are available at the moment, so some reasonable estimates need to be made.

In Eqs. (12), the coefficients for the collisional losses, $`\gamma _j`$, are measured the units of $`\mathrm{\Gamma }=x_0^3\omega _{}`$, and for the typical values of $`x_05`$ $`\mu \mathrm{m}`$, and $`\omega _{}10^2`$ $`\mathrm{Hz}`$, $`\mathrm{\Gamma }10^8`$ $`\mathrm{cm}^3/\mathrm{s}`$. In what follows, we investigate the effect of relatively “weak”, $`\gamma _j=10^3`$ and “strong”, $`\gamma _j=10^1`$ losses, which nevertheless corresponds to physically “strong” inelastic collisions rates of $`\mathrm{\Gamma }=10^{11}10^9`$ $`\mathrm{cm}^3/\mathrm{s}`$. For simplicity, we assume that all the inelastic loss rates in the model are equal to each other.

As was estimated in Ref. , the linear loss rates that still enable the coherent dynamics of the AMBEC are of the order of $`10^2`$ $`\mathrm{s}^1`$, which corresponds to our dimensionless values of $`\nu _j1`$. Assuming this to be the case, we note that the term responsible for the “weakest” nonlinear losses in our model still an order of magnitude larger than the corresponding linear loss, i.e. $`\mathrm{\Gamma }n_j10^3`$ $`\mathrm{s}^1`$ for the typical condensate density of $`n_j10^{14}`$ $`\mathrm{cm}^3`$. This means that, without loss of generality, we can assume that $`\gamma _j=\nu _j`$ in our model, and bear in mind that the effect of the linear losses is only significant when the particle numbers (and condensate densities) are small.

The other (strong) loss mechanism for the molecular fraction, which is absent in our effectively two-mode consideration, is the coupling of a molecule to the other, excited states of the atomic BEC due to the closeness of the excited molecular level to the BEC atom-pair level. This loss channel has been recently shown to be responsible for the enhancement of the condensate loss rates . One can argue, however, that on the time scales large compared with the corresponding rates of decay into the higher-order modes, the lowest-order nonlinear excited modes of the hybrid AMBEC that we consider here become decoupled from the rest of the bound state spectrum, and all losses associated with such coupling can be effectively incorporated into the phenomenological terms proportional to $`\nu _m`$.

### B Lossy dynamics of the vortex soliton

In the presence of losses, the total particle number in both components of AMBEC decreases. First of all, when the molecular component is small, the effective equations (6) can be decoupled and the dissipative dynamics of two components can be described separately. In particular, neglecting the molecular component, and averaging over the spatial structure of the atomic BEC wavefunction, from Eqs. (6),(12) we obtain the simple equation for the density $`n_a=|\psi _a|^2`$ of the atomic component

$$i\frac{dn_a}{dt}=2\nu _an_a2\gamma _an_a^2,$$
(13)
which has the analytic solution,

$$n_a(t)=\frac{\nu _ae^{2\nu _at}}{(C\gamma _ae^{2\nu _at})},$$
(14)
where $`C=\gamma _a+\nu _a/n_a(0)`$. When the effect of linear losses can be neglected, the decay of the atomic component follows a power law.

We have conducted the numerical simulations of the dissipation-induced dynamics of the vortex soliton varying the loss terms. In Fig. 3, we show the decay of the maximum amplitude of the atomic component $`n_a^{\mathrm{max}}=|\psi _a|_{\mathrm{max}}^2`$ of the two-component vortex soliton in the presence of the losses of a different strength. For “strong” and “medim” losses ($`\gamma =\nu =0.1`$ and $`0.01`$ respectively), the decay follows approximately the analytical dependence (14) and, therefore, the modal (spatial) structure of the vortex state is not very important.

For weaker losses, we reveal novel features in the vortex dynamics associated with the interplay between the instability-induced vortex fragmentation and the decay of the vortex amplitude due to losses. As an example, in Fig. 4(top) we show the propagation of the initially stable, untrapped, vortex soliton \[shown as (f) in Fig. 1\] in the presence of weak dissipation, $`\gamma _i=\nu _i=0.001`$. While this vortex is absolutely stable in the lossless regime (see Fig. 2, top), its number of particles decreases and the corresponding vortex solution moves along the dashed curve in Fig. 1 from the point (f) to the point (e), entering the lower part of the branch corresponding to unstable vortices. Then, the vortex fragmentation into two or three filaments occurs, as shown in Fig. 4(top). This leads us to conclude that the fragmentation to components occurs when the particle number has dropped far enough to place the vortex in the region of instability for stationary solutions and, overall, the vortex parameters adiabatically follow the family of untrapped localized solutions.

For comparison, in Fig. 4(bottom) we show the vortex dynamics for stronger losses, when no fragmentation is observed at all. In this case, the vortex soliton does not follow a family of localized solutions adiabatically. Moreover, it appears that the dissipation leads to a switching from the fundamental vortex state to one of the higher-order vortex modes. The additional “rings” in the density distribution observed in Fig. 4 (bottom) therefore seem to be a manifestation of such a higher-order mode.

### C Evolution of a vortex released from a trap

Although different possibilities to create the topological states in experiment may exist, it is probably safe to assume that the presence of the confining trap is paramount. The question therefore is whether a dynamically stable topological state of AMBEC created in a trap can be transformed into a dynamically stable vortex in a self-confined AMBEC droplet, after the hybrid condensate is released from the trap. We investigate such a possiblity by conducting additional numerical studies.

A trapped vortex, corresponding to the stable family shown in the main part of Fig. 1 \[the branch that includes the point (d)\] is released from the trap at $`t=0`$. We find that for weak or no dissipation this vortex structure diffracts rapidly, as shown in Fig. 5(top). However, when the losses become large enough, the released vortex state gets self-trapped, as shown in Fig. 5(bottom).

In order to verify the effect of self-trapping and get a deeper insight into our results, we use the parametric plot of the system Hamiltonian $`H(t)`$ vs. the total particle number $`N(t)`$, which both evolve in time in the presence of losses. For weak losses, we reveal that the vortex dynamics follows the “trapped” family of the vortex states (dashed line in Fig. 6). The family of the trapped vortex states in this plot corresponds to the initial conditions of a trapped vortex stationary solution with the trapping potential turned off. This vortex state diffracts when it is released from the trap due to the dominant repulsive cubic nonlinearity. However, for stronger losses we observe a novel effect when the vortex, being initially released from the trap, does not rapidly diffract but, instead, it switches to a self-trapped vortex soliton, as shown by a dotted curve in Fig. 6.

## V Conclusions

In the framework of a two-component mean-field model of hybrid coherently coupled condensates, we have analyzed the structure, stability, and dynamics of the topological states in the atomic-molecular BEC created via the induced two-color Raman photoassociation process. We have predicted the existence of novel two-component vortex solitons supported by the coherent parametric coupling between the atomic and molecular BEC fractions, and demonstrated that such topological states can exist even in the absence of a trapping potential.

Since the losses are known to be very important in the formation of the hybrid atomic-molecular condensates, we have analyzed the effect of weak and strong losses of a different physical origin on the dynamics of the vortex solitons, the type of problems never addressed in other fields such as nonlinear optics. In particular, for the first time to our knowledge, we have described an interplay between the nonlinearity-induced vortex fragmentation and dissipation-induced vortex decay. Moreover, we have revealed that the losses may be very useful to identify the self-trapped vortex states in the hybrid condensates in experiment, after releasing the atomic-molecular condensate from a trap.

Although our analysis is valid for a quasi-two-dimensional model, the generalization to a spherically symmetric case seems straightforward. Additionally, in spite of the limited validity of the two-component model, we believe that the similar types of topological states can be found in more realistic models of hybrid condensates, as well as in other types of coherently coupled mixtures of the multi-component condensates such as spinor Bose-Einstein condensates.

The work has been supported by the Planning and Performance Fund of the Institute of Advanced Studies at the Australian National University.