text:
# Superconductivity in the Orbital Degenerate Model for Heavy Fermion Systems

## I Introduction

Superconductivity in strongly correlated electron systems has been one of central issues in the research field of condensed matter physics, since the pioneering discovery of superconductivity in CeCu2Si2 . The subsequent discovery of high temperature superconductivity in cuprates has accelerated further investigations of this subject, leading to unambiguous identification of the unconventional nature of $`d_{x^2y^2}`$-wave pairing in cuprates. Recently, new heavy fermion compounds CeTIn5 (T=Rh, Ir, and Co) have been discovered . Among them, CeRhIn5 exhibits an antiferromagnetic (AFM) transition at a Néel temperature $`T_\mathrm{N}`$=3.8K and becomes superconducting only under hydrostatic pressure larger than 15 kbar. On the other hand, both CeIrIn5 and CeCoIn5 are superconducting at ambient pressure with transition temperatures $`T_\mathrm{c}`$=0.4K and 2.3K, respectively. Especially, the $`T_\mathrm{c}`$ of CeCoIn5 is the highest among those yet observed for heavy fermion superconductors .

Reflecting the fact that CeTIn5 has the HoCoGa5-type tetragonal crystal structure, quasi two-dimensional Fermi surfaces have been observed in de Haas-van Alphen experiments for the compounds, consistent with the band-structure calculation results . Concerning the superconducting state, nuclear relaxation rate in CeTIn5 exhibits $`T^3`$ behavior below $`T_\mathrm{c}`$  and thermal conductivity in CeCoIn5 is found to include a component with four-fold symmetry , strongly suggesting $`d_{x^2y^2}`$-wave pairing in the superconducting phase of CeTIn5. Furthermore, it has been shown that in the alloy system CeRh1-xIrxIn5, the superconducting phase appears in the neighborhood of the AFM phase . A natural consequence of these experimental results is that superconductivity in CeTIn5 compounds is induced by AFM spin fluctuations, similarly to high-$`T_\mathrm{c}`$ superconductivity.

In spite of the similarities mentioned above, here we should emphasize several differences between heavy fermion superconductors such as CeTIn5 and high-$`T_\mathrm{c}`$ cuprates. First of all, relevant electrons in high-$`T_\mathrm{c}`$ superconductors are almost itinerant $`3d`$-electrons, while for heavy fermion superconductors, relevant ones are $`f`$-electrons, which almost localized and their dispersion is mainly determined by hybridization with conduction electrons. Another important difference is concerned with the electronic states relevant to low-energy physics. For cuprates, it is widely recognized that a single-band model for holes in $`d_{x^2y^2}`$ orbitals on the square lattice is a good starting point. On the other hand, for heavy fermion superconductors, several Fermi surfaces are observed in general. Such complex electronic states may be eventually traced back to the orbital degeneracy and relatively weak crystalline electric field (CEF) effect of $`f`$-electrons compared with the $`3d`$-electrons in the CuO2 plane. It implies that construction of a realistic microscopic model is not an easy work for heavy fermion systems. Therefore, up to now, most of theoretical investigations for superconductivity in heavy fermion systems have been restricted in a phenomenological level.

In this paper, we discuss the effect of orbital fluctuations on superconductivity based on a microscopic theory applied to a microscopic model. In the next section, we introduce an orbital degenerate model obtained including important characters of CeTIn5. Then, in order to study the superconducting transition induced by spin and/or orbital fluctuations, we develop a strong-coupling theory using fluctuation exchange (FLEX) approximation  in which spin and orbital fluctuations as well as the single-particle spectrum are determined self-consistently. Finally, we discuss experimental results for CeTIn5 in the light of the present theory.

## II Model Hamiltonian

In order to introduce a minimal model for $`f`$-electron systems, we start from the local basis for a Ce3+ ion. Among 14-fold degenerate $`4f`$-electronic states, due to the effect of strong spin-orbit coupling, only the $`j`$=5/2 sextet effectively contributes to the low-energy excitations ($`j`$ is total angular momentum). Under the cubic CEF, this sextet is further split into $`\mathrm{\Gamma }_7`$ doublet and $`\mathrm{\Gamma }_8`$ quartet, where the corresponding eigen-states are given by $`|\mathrm{\Gamma }_{7\pm }`$= $`\sqrt{1/6}|\pm 5/2`$$``$$`\sqrt{5/6}|3/2`$, $`|\mathrm{\Gamma }_{8\pm }^{(1)}`$= $`\sqrt{5/6}|\pm 5/2`$+$`\sqrt{1/6}|3/2`$, and $`|\mathrm{\Gamma }_{8\pm }^{(2)}`$=$`|\pm 1/2`$. Here $`+`$ and $``$ in the subscripts denote “pseudo-spin” up and down in each Kramers doublet, respectively. For simplicity, we consider only $`|\mathrm{\Gamma }_8`$ states, further split into two Kramers doublets as $`|\mathrm{\Gamma }_{8\pm }^{(1)}`$ and $`|\mathrm{\Gamma }_{8\pm }^{(2)}`$ under the additional tetragonal CEF. Note that $`|\mathrm{\Gamma }_{8\pm }^{(1)}`$ and $`|\mathrm{\Gamma }_{8\pm }^{(2)}`$ belong to $`\mathrm{\Gamma }_7`$ and $`\mathrm{\Gamma }_6`$ irreducible representations, respectively, in the tetragonal system. Validity of this assumption for CeTIn5 will be discussed later in this paper, but here we stress that the Hamiltonian constructed from $`\mathrm{\Gamma }_8`$ quartet is the simplest model including essential physics of interplay between spin and orbital fluctuations.

In order to consider itinerant features of $`4f`$-electrons, we take into account nearest-neighbor hopping of $`f`$-electrons by the tight-binding method . It should be noted that the matrix elements of the nearest-neighbor hopping depend on not only the orbital but also the hopping direction, since the forms of wavefunctions of $`|\mathrm{\Gamma }_{8\pm }^{(1)}`$ and $`|\mathrm{\Gamma }_{8\pm }^{(2)}`$ states are different from each other. We can determine the hopping matrix elements by assuming that the overlap integral through the $`\sigma `$-bond ($`ff\sigma `$) is dominant. Since CeTIn5 has a tetragonal crystal structure and quasi two-dimensional Fermi surfaces have been experimentally observed , it is natural to consider the two-dimensional square lattice composed of Ce3+ ions. Although the hybridization of $`f`$ electrons with In $`5p`$ electronic states may be the main source of kinetic energy for $`f`$ electrons, in the present scheme such an effect is considered as the effective hopping amplitude for $`f`$ quasi-particles, after the $`p`$-electron degrees of freedom are integrated out.

By further adding the on-site Coulomb interaction terms among $`f`$-electrons, the effective Hamiltonian with orbital degeneracy is given by

$`H`$ $`=`$ $`{\displaystyle \underset{\mathrm{𝐢𝐚}\tau \tau ^{}\sigma }{}}t_{\tau \tau ^{}}^𝐚f_{𝐢\tau \sigma }^{}f_{𝐢+𝐚\tau ^{}\sigma }\epsilon {\displaystyle \underset{𝐢}{}}(n_{𝐢1\sigma }n_{𝐢2\sigma })/2`$ (1)
$`+`$ $`U{\displaystyle \underset{𝐢\tau }{}}n_{𝐢\tau }n_{𝐢\tau }+U^{}{\displaystyle \underset{𝐢\sigma \sigma ^{}}{}}n_{𝐢1\sigma }n_{𝐢2\sigma ^{}},`$ (2)
where $`f_{𝐢\tau \sigma }`$ is the annihilation operator for an $`f`$-electron with pseudo-spin $`\sigma `$ in the $`\tau `$-orbital state $`\mathrm{\Gamma }_8^{(\tau )}`$ at site $`𝐢`$, $`𝐚`$ is the vector connecting nearest-neighbor sites, and $`n_{𝐢\tau \sigma }`$=$`f_{𝐢\tau \sigma }^{}f_{𝐢\tau \sigma }`$. The first term represents the nearest-neighbor hopping of $`f`$-electrons with the amplitude $`t_{\tau \tau ^{}}^𝐚`$ between $`\tau `$ and $`\tau ^{}`$ orbitals along the $`𝐚`$-direction, given by $`t_{11}^𝐱`$=$`\sqrt{3}t_{12}^𝐱`$= $`\sqrt{3}t_{21}^𝐱`$=$`3t_{22}^𝐱`$=1 for $`𝐚`$=$`𝐱`$ and $`t_{11}^𝐲`$=$`\sqrt{3}t_{12}^𝐲`$= $`\sqrt{3}t_{21}^𝐲`$=$`3t_{22}^𝐲`$=1 for $`𝐚`$=$`𝐲`$, respectively, in energy units where $`t_{11}^𝐱`$=1. The second term denotes the tetragonal CEF, leading to an energy splitting $`\epsilon `$ between the two orbitals. In the third and fourth terms, $`U`$ and $`U^{}`$ are the intra- and inter-orbital Coulomb interactions, respectively. Due to the rotational invariance in the orbital space for the interaction part, $`U^{}`$ should be equal to $`U`$, when we ignore the Hund’s rule coupling, since it is irrelevant in the large-$`U`$ limit for the quarter-filling case with one $`f`$-electron per site. Thus, in this paper, we restrict ourselves to the case of $`U`$=$`U^{}`$. Note also that in the quarter-filling case, the present model is virtually reduced to the half-filled single-orbital Hubbard model in the limit of $`\epsilon `$=$`\mathrm{}`$.

## III FLEX approximation

In our previous work, we have developed a weak-coupling theory for superconductivity based on the same orbital degenerate model described above, using the static spin and orbital fluctuations obtained within the random phase approximation (RPA) . Various superconducting phase have been found around varieties of ordered phases, whose boundaries are determined by the RPA instability. In order to develop more sophisticated theory, we should include effects of (1) mode-mode coupling and (2) quasi-particle damping, neglected in the previous work. Regarding (1), the RPA does not incorporate effects of mode-mode coupling between fluctuations, but the mode-mode coupling modifies significantly the temperature and frequency dependences of spin and orbital fluctuations. Concerning (2), within the weak-coupling theory for superconductivity, damping of quasi-particles by the scattering due to spin and orbital fluctuations is not taken into account, which leads to the suppression of superconductivity. Therefore, the weak-coupling theory generally overestimates the region of superconducting phase.

In the present paper, we apply the fluctuation exchange (FLEX) approximation  to the orbital degenerate model discussed in the preceding section. The FLEX approximation has the following two features: (1) This is a kind of the mode-mode coupling theory, where spin and orbital fluctuations and the spectra of $`f`$-electrons are self-consistently determined through the fluctuation exchange self-energy. (2) It provides the Dyson-Gorkov equation where the normal and anomalous self-energies are obtained on an equal footing. Here we emphasize that the FLEX approximation has been successful to understand consistently the normal and superconducting states of high-$`T_\mathrm{c}`$ cuprates, in particular with reasonable estimation of $`T_\mathrm{c}`$ . Therefore, it is interesting to apply the FLEX approximation to the orbital degenerate model to understand its properties concerning superconductivity.

In the doubly degenerate case, the Green’s functions for $`f`$-electrons form a 2$`\times `$2 matrix and they follow the Dyson-Gorkov equations

$$\widehat{G}(k)=\widehat{G}^{(0)}(k)+\widehat{G}^{(0)}(k)\widehat{\mathrm{\Sigma }}^{(1)}(k)\widehat{G}(k),$$
(3)
where $`\widehat{G}^{(0)}(k)`$ is a matrix for the non-interacting Green’s function. Within the FLEX approximation, components of the self-energy matrix $`\mathrm{\Sigma }_{ml}^{(1)}(k)`$ are given by

$$\mathrm{\Sigma }_{ml}^{(1)}(k)=\frac{T}{N_0}\underset{q}{}\underset{\mu \nu }{}V_{\mu m,\nu l}^{\mathrm{eff}}(q)G_{\mu \nu }(kq),$$
(4)
with

$`V_{\mu m,\nu l}^{\mathrm{eff}}(q)=[{\displaystyle \frac{3}{2}}\widehat{U}^\mathrm{s}\widehat{\chi }^\mathrm{s}(q)\widehat{U}^\mathrm{s}+{\displaystyle \frac{1}{2}}\widehat{U}^\mathrm{o}\widehat{\chi }^\mathrm{o}(q)\widehat{U}^\mathrm{o}`$ (5)
$`{\displaystyle \frac{1}{4}}(\widehat{U}^\mathrm{s}+\widehat{U}^\mathrm{o})\widehat{\overline{\chi }}(q)(\widehat{U}^\mathrm{s}+\widehat{U}^\mathrm{o})+{\displaystyle \frac{3}{2}}\widehat{U}^\mathrm{s}{\displaystyle \frac{1}{2}}\widehat{U}^\mathrm{o}]_{\mu m,\nu l},`$ (6)
where the first and second terms of $`V_{\mu m,\nu l}^{\mathrm{eff}}(q)`$ give contributions of the spin and orbital fluctuations, respectively, to the self-energy.

In Eq. (4), $`\widehat{\chi }^\mathrm{s}(q)`$ and $`\widehat{\chi }^\mathrm{o}(q)`$ are the 4$`\times `$4 matrices of the spin and orbital fluctuations, given by

$`\widehat{\chi }^\mathrm{s}(q)=[\widehat{1}\widehat{U}^\mathrm{s}\widehat{\overline{\chi }}(q)]^1\widehat{\overline{\chi }}(q),`$ (7)
$`\widehat{\chi }^\mathrm{o}(q)=[\widehat{1}+\widehat{U}^\mathrm{o}\widehat{\overline{\chi }}(q)]^1\widehat{\overline{\chi }}(q),`$ (8)
where the matrix element of the irreducible susceptibility $`\overline{\chi }_{ij,st}(q)`$ is given by

$$\overline{\chi }_{ij,st}(q)=\frac{T}{N_0}\underset{k}{}G_{si}(k+q)G_{jt}(k),$$
(9)
For the susceptibility matrices, the labels of row and column appear in the order 11, 22, 12, and 21, pairs of orbital indices 1 and 2. In these expressions, $`T`$ is a temperature, $`N_0`$ is number of unit cells, abbreviations $`k`$=$`(𝐤,i\omega _n)`$ and $`q`$=$`(𝐪,i\mathrm{\Omega }_n)`$ ($`\omega _n`$=$`(2n+1)\pi T`$, $`\mathrm{\Omega }_n`$=$`2n\pi T`$) are used. The interaction matrices $`\widehat{U}^\mathrm{s}`$ and $`\widehat{U}^\mathrm{o}`$ are given by

$`\widehat{U}^\mathrm{s}=\left[\begin{array}{cccc}U& 0& 0& 0\\ 0& U& 0& 0\\ 0& 0& U^{}& 0\\ 0& 0& 0& U^{}\end{array}\right],\widehat{U}^\mathrm{o}=\left[\begin{array}{cccc}U& 2U^{}& 0& 0\\ 2U^{}& U& 0& 0\\ 0& 0& U^{}& 0\\ 0& 0& 0& U^{}\end{array}\right].`$ (18)
By solving the self-consistent equations, the spin and orbital fluctuations and the Green’s functions for $`f`$-electrons are determined simultaneously.

In order to discuss superconductivity, it is necessary to calculate anomalous self-energy. The matrix element of the anomalous self-energy $`\mathrm{\Sigma }_{ml}^{(2)}(k)`$ is obtained by the functional derivative of the thermodynamic potential with respect to the anomalous Green’s function as follows

$$\mathrm{\Sigma }_{ml}^{(2)}(k)=\frac{T}{N_0}\underset{q}{}\underset{\mu \nu }{}V_{\mu m,l\nu }^\xi (q)F_{\mu \nu }(kq),$$
(19)
with

$$\widehat{F}(k)=\widehat{G}(k)\widehat{\mathrm{\Sigma }}^{(2)}(k)\widehat{G}(k),$$
(20)
where the matrix elements of the effective pairing interactions for spin-singlet and spin-triplet channels are given, respectively, by

$`V_{\mu m,l\nu }^\mathrm{S}(q)=[{\displaystyle \frac{3}{2}}\widehat{U}^\mathrm{s}\widehat{\chi }^\mathrm{s}(q)\widehat{U}^\mathrm{s}+{\displaystyle \frac{1}{2}}\widehat{U}^\mathrm{o}\widehat{\chi }^\mathrm{o}(q)\widehat{U}^\mathrm{o}`$ (21)
$`{\displaystyle \frac{1}{2}}(\widehat{U}^\mathrm{s}+\widehat{U}^\mathrm{o})]_{\mu m,l\nu },`$ (22)
$`V_{\mu m,l\nu }^\mathrm{T}(q)=[{\displaystyle \frac{1}{2}}\widehat{U}^\mathrm{s}\widehat{\chi }^\mathrm{s}(q)\widehat{U}^\mathrm{s}+{\displaystyle \frac{1}{2}}\widehat{U}^\mathrm{o}\widehat{\chi }^\mathrm{o}(q)\widehat{U}^\mathrm{o}`$ (23)
$`{\displaystyle \frac{1}{2}}(\widehat{U}^\mathrm{s}+\widehat{U}^\mathrm{o})]_{\mu m,l\nu }.`$ (24)
The $`T_\mathrm{c}`$ is obtained by the temperature at which the maximum eigenvalue of Eq. (8) becomes unity. As already pointed out in the previous works based on the RPA , one can see from these effective pairing interactions that developments of both spin and orbital fluctuations have a destructive interference for the singlet channel, while they are constructive for the triplet one.

## IV Calculated results

The FLEX calculation is numerically carried out for each value of $`\epsilon `$ at fixed parameter values of $`U`$=$`U^{}`$=4 and $`n`$=1 corresponding to one $`f`$-electron density per site. All summations involved in the above self-consistent equations are performed using the fast Fourier transformation algorithm for the $`𝐤`$-space with $`32\times 32`$ meshes in the first Brillouin zone and for Matsubara frequency sum with energy cut-off five times larger than the relevant band width. In Fig. 1, $`𝐪`$ dependences of the principal components of $`\widehat{\chi }^\mathrm{s}(𝐪,0)`$ and $`\widehat{\chi }^\mathrm{o}(𝐪,0)`$ are shown for a fixed temperature $`T`$=0.02 for different level splitting, $`\epsilon `$: The upper panel for $`\epsilon `$=0 and the lower panel for $`\epsilon `$=2. For $`\epsilon `$=0, corresponding to the orbitally degenerate case, the AFM spin fluctuation in the $`\tau `$=1 orbital is enhanced, but not sufficiently developed to induce $`d_{x^2y^2}`$-wave superconductivity. With increasing the orbital splitting energy to $`\epsilon `$=2, the AFM spin fluctuation for the $`\tau `$=1 orbital further develops, and orbital fluctuations are completely suppressed, in comparison with the developed AFM spin fluctuation.

In Fig. 2, the phase diagram obtained within the FLEX approximation for the orbital degenerate system are shown, where the solid and open circles describe the superconducting and AFM transitions, respectively. The dotted curve is a schematic phase boundary expected between the two phases. From Fig. 2, we can see that (1) the spin-singlet superconducting phase with $`B_{1\mathrm{g}}`$-symmetry appears next to the AFM phase and (2) $`T_\mathrm{c}`$ is enhanced with increasing the orbital splitting energy $`\epsilon `$. From these observations, we can conclude that the superconducting phase is induced by the development of the AFM spin fluctuations for the $`\tau `$=1 orbital with increasing the orbital splitting energy $`\epsilon `$.

## V Summary and discussion

Let us now discuss the experimental results for CeTIn5 in comparison with the present theoretical results. First we should note that in the actual tetragonal crystal, the four-fold degenerate $`\mathrm{\Gamma }_8`$ states in the cubic notation split into two Kramers doublets, $`\mathrm{\Gamma }_7`$ and $`\mathrm{\Gamma }_6`$, as we mentioned in the Sec. 2. Therefore, in the tetragonal system, the $`j`$=5/2 states split into two $`\mathrm{\Gamma }_7`$ and one $`\mathrm{\Gamma }_6`$. Analyses of experimental data of magnetic susceptibility of CeTIn5 by using the CEF theory seem to be consistent more with the level scheme where two $`\mathrm{\Gamma }_7`$ are lower than the $`\mathrm{\Gamma }_6`$ . The energy splitting between the two $`\mathrm{\Gamma }_7`$ is estimated as 68K for CeRhIn5, 61K for CeIrIn5, and 151K for CeCoIn5. According to the present analysis, higher $`T_\mathrm{c}`$ is obtained for larger $`\epsilon `$, consistent with the tendency in $`T_\mathrm{c}`$ of the two superconducting materials: $`T_\mathrm{c}`$=0.4K for CeIrIn5 and $`T_\mathrm{c}`$=2.3K for CeCoIn5.

However, the antiferromagnetically ordered CeRhIn5 with $`T_\mathrm{N}`$=3.8K has an intermediate value for $`\epsilon `$. One possible scenario to understand the discrepancy is to consider difference of quasi-two dimensionality, as we have pointed out in . To make more direct and quantitative comparison with experimental results, however, one should be aware of an assumption for the model used in the present study, which may be called as $`\mathrm{\Gamma }_8`$ model. Namely, this $`\mathrm{\Gamma }_8`$ model assumes that $`\mathrm{\Gamma }_7`$ and $`\mathrm{\Gamma }_6`$ are the lower two Kramers doublets, which may be different from the level scheme obtained from the experiments. Thus, it will be an interesting future problem to elucidate the role of orbital splitting for the more realistic two $`\mathrm{\Gamma }_7`$ model, which includes not only the nearest neighbor hopping but also the next nearest neighbor one and thus, may reproduce the realistic electronic states better than the $`\mathrm{\Gamma }_8`$ model .

In order to construct a whole story for CeTIn5 compounds, it may be required to use the even more realistic $`f`$-$`p`$ model including the $`f`$-$`p`$ hybridization explicitly. In the $`f`$-$`p`$ model, we can have two different energy scales: The large energy scale is the band width for the conduction electron and/or the Coulomb interaction, while the small energy scale corresponds to the energy splitting of the CEF levels discussed here. The present study indicates the possibility that the small energy scale of the CEF level splitting plays a key role as the controlling parameter of quantum phase transitions.

In summary, based on the effective microscopic model with orbital degeneracy for $`f`$-electron systems, we have proposed that the orbital splitting energy is the parameter controlling the change from the paramagnetic to the AFM phase with the $`d_{x^2y^2}`$-wave superconducting phase in between. Actually, the self-consistent FLEX approximation applied to the orbital degenerate model shows that the $`d_{x^2y^2}`$-wave superconducting phase is induced by increasing the orbital splitting energy.

## Acknowledgement

The authors would like to thank T. Maehira for many helpful discussions. T.H. and K.U. are supported by the Grant-in-Aid for Scientific Research from Japan Society for the Promotion of Science.

## References


format: text

text:
# Stability and Decay Rates of Non-Isotropic Attractive Bose-Einstein Condensates

## I Introduction

Experimental Bose-Einstein condensation (BEC) with attractive interactions was first realized in ultra cold vapors of $`^7`$Li atoms Bradley1 , opening a new field in the study of macroscopic quantum phenomena. Such attractive condensates are known to be metastable in spatially localized systems, provided that the number of condensed particles is below a critical value $`𝒩_c`$ Bradley2 . More recently, Feshbach resonances in BEC of $`^{85}`$Rb atoms were used to investigate the stability and dynamics of condensates with two-body interactions going from repulsive to attractive values PRLRoberts .

Experimental atomic traps generally use a harmonic and slightly asymmetric potential. Thus, for most of the condensates produced so far, the geometry is nearly spherical. However, extremely asymmetric traps have been recently employed in experimental investigations of cigar-like expcig1 ; expcig2 ; expcig3  or pancake-like expcrepe  condensates.

Various physical processes compete to determine the lifetime of attractive condensates. The processes considered in this paper are macroscopic quantum tunneling (MQT) Stoof1 ; Ueda , inelastic two and three body collisions (ICO) Shi ; Dodd ; Kagan  and thermally induced collapse (TIC) Stoof1 ; Sackett . The MQT and TIC contributions have been evaluated in the literature using a variational Gaussian approximation to the condensate wave function. However, this approximation is known to be in substantial quantitative error – e.g. as high as $`17\%`$ for $`𝒩_c`$ Ueda ; Ruprecht ; PRLiso  – when compared to the exact solution of the Gross-Pitaevskii (G-P) equation.

In the nearly spherical isotropic case, both the elliptic (stable) and the hyperbolic (unstable) exact stationary solutions of the G-P equation were obtained numerically by Newton’s method in PRLiso . These solution branches were shown to meet at $`𝒩_c`$ through a generic Hamiltonian Saddle Node (HSN) bifurcation. While the Gaussian approximation presents an analogous HSN bifurcation, the amplitudes of its associated scaling laws were found to be in substantial ($`14\%`$) error. A method for computing the unstable branch in the isotropic case via a shooting method was outlined in Gammal1999 , but generalizing this procedure to higher dimensions would be inefficient, and impossible in non-rectangular domains. The decay rates for the processes of MQT, ICO and TIC were also computed, in the spherical case, from the numerical G-P solutions. They were shown to obey universal scaling laws. Experimentally significant quantitative differences were found between the exact rates and those based on the Gaussian approximation PRLiso .

In the extreme anisotropic cases, the variational Gaussian approximation has been computed and compared to the G-P solution on the elliptic (stable) branch GammalPRA64 ; GammalPreprint . This has allowed a more reliable determination of the critical value $`𝒩_c`$ than can be obtained by the Gaussian approximation GammalPreprint . However, a faithful determination of the lifetimes needs the computation of the hyperbolic (unstable) branch PRLiso , which has not yet been performed in the anisotropic case.

The main purpose of the present paper is to show that it is possible to compute the full HSN bifurcation diagram, and the corresponding lifetimes, in extreme anisotropic cases. We will do so by studying a cigar-like and a pancake-like condensate, and will obtain their MQT, ICO and TIC decay rates. While we have concentrated, for simplicity, on these two axisymmetric cases, the new numerical methods developed in this work are capable of solving the general anisotropic problem.

The paper is organized as follows. In section II we present the model considered throughout this work. After defining our working form of the G-P equation, we explain the methods that we used to obtain the stationary states and their linearized stability. Section III is devoted to the numerical determination of the bifurcation diagram and stability of the stationary states. Isotropic and non-isotropic cases are compared and the dynamics is discussed in terms of the HSN bifurcation. In section IV we define and compute the decay rates. Isotropic and non-isotropic rates are discussed, their similarity is analyzed in terms of the spontaneous isotropization of condensates. Finally, section V is our conclusion. Details of our numerical methods are given in the Appendix.

## II Presentation of the Model

### II.1 Gross-Pitaevskii Equation

At low enough temperatures, neglecting the thermal and quantum fluctuations, a Bose condensate can be represented by a complex wave function $`\mathrm{\Psi }(𝐱,t)`$ that obeys the dynamics of the G-P equation Gross0 ; Pitaevskii1 . Specifically, we consider a condensate of $`𝒩`$ particles of mass $`m`$ and (negative) effective scattering length $`\stackrel{~}{a}`$ in a confining harmonic potential $`V(\stackrel{~}{𝐱})=m(\stackrel{~}{\omega }_x^2\stackrel{~}{x}^2+\stackrel{~}{\omega }_y^2\stackrel{~}{y}^2+\stackrel{~}{\omega }_z^2\stackrel{~}{z}^2)/2`$ where $`\stackrel{~}{𝐱}=(\stackrel{~}{x},\stackrel{~}{y},\stackrel{~}{z})`$ is the position vector.

These variables can be rescaled with respect to any reference frequency $`\widehat{\omega }`$ by using the natural quantum harmonic oscillator units of time $`\tau _0=1/\widehat{\omega }`$ and length $`L_0=\sqrt{\mathrm{}/m\widehat{\omega }}`$. In terms of the non-dimensional variables $`t=\stackrel{~}{t}/\tau _0`$, $`𝐱=\stackrel{~}{𝐱}/L_0`$, $`a=4\pi \stackrel{~}{a}/L_0`$, $`\omega _x=\stackrel{~}{\omega }_x/\widehat{\omega }`$, $`\omega _y=\stackrel{~}{\omega }_y/\widehat{\omega }`$ and $`\omega _z=\stackrel{~}{\omega }_z/\widehat{\omega }`$, the condensate is described by the action

$$𝒜=𝑑t\left\{d^3x\frac{i}{2}\left(\overline{\mathrm{\Psi }}\frac{\mathrm{\Psi }}{t}\mathrm{\Psi }\frac{\overline{\mathrm{\Psi }}}{t}\right)\right\},$$
(1)
with

$$=\mu 𝒩$$
(2)
where

$`𝒩`$ $`=`$ $`{\displaystyle d^3x|\mathrm{\Psi }|^2}`$ (3)
$``$ $`=`$ $`{\displaystyle d^3x\left[\frac{1}{2}|\mathrm{\Psi }|^2+V(𝐱)|\mathrm{\Psi }|^2+\frac{a}{2}|\mathrm{\Psi }|^4\right]}`$ (4)
$`V(𝐱)`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\omega _x^2x^2+\omega _y^2y^2+\omega _z^2z^2)`$ (5)
The Euler-Lagrange equation corresponding to $`𝒜`$ is our working form of the Gross-Pitaevskii equation:

$`i{\displaystyle \frac{\mathrm{\Psi }}{t}}`$ $`=`$ $`{\displaystyle \frac{\delta }{\delta \overline{\mathrm{\Psi }}}}`$
$`=`$ $`\left[{\displaystyle \frac{1}{2}}^2V(𝐱)\left(a|\mathrm{\Psi }|^2\mu \right)\right]\mathrm{\Psi }.`$
Our goal is to numerically determine the stable and unstable stationary states of (II.1) and the eigenvalues of (II.1) linearized about these stationary states. We will carry out this calculation for various values of a cylindrical potential defined by $`\omega _r\omega _x=\omega _y`$ and $`\omega _z`$: the isotropic case $`\omega _r=\omega _z`$, a cigar case $`\omega _r/5=\omega _z`$ and a pancake case $`\omega _r=\omega _z/5`$. We will then use these results to calculate the condensate decay rates and compare these decay rates to those produced by the Gaussian approximation.

### II.2 Stationary States

Stationary states of (II.1) corresponding to minima of $``$ at a given value of $`𝒩`$ can be obtained by integrating to relaxation the diffusion equation

$`{\displaystyle \frac{\mathrm{\Psi }}{t}}`$ $`=`$ $`{\displaystyle \frac{\delta }{\delta \overline{\mathrm{\Psi }}}}=\left[{\displaystyle \frac{1}{2}}^2V(𝐱)\left(a|\mathrm{\Psi }|^2\mu \right)\right]\mathrm{\Psi }`$ (7)
$`{\displaystyle \frac{𝒩}{t}}`$ $`=`$ $`0,`$ (8)
using initial data $`\mathrm{\Psi }(t=0)`$ with a total number of particles $`𝒩`$. The condition (8) fixes the value of the Lagrange multiplier $`\mu `$ during the relaxation. This relaxation method yields both the solution $`\mathrm{\Psi }`$ and the Lagrange multiplier $`\mu `$. It is equivalent to that used in Ruprecht  and can only reach the stable stationary solutions of (7). Unstable stationary solutions to (II.1) and (7) are obtained by a Newton branch following method detailed in the Appendix.

Note that the Lagrange multiplier $`\mu `$ can only affect the solutions of (II.1) through a homogeneous rotating phase factor $`e^{i\mu t}`$, in contrast to its particle number conservation effect on equation (7). However, every stationary solution to (II.1) is indexed by the unique $`\mu `$-value that makes it time-independent, as shown in figure 1.

### II.3 Linearized Stability

We now turn our attention to computing the linear stability of the G-P equation about a stationary state. We first write (II.1) in the abbreviated form:

$$i\frac{\mathrm{\Psi }}{t}=L\mathrm{\Psi }+W(\mathrm{\Psi })$$
(9)
where

$`L\mathrm{\Psi }`$ $``$ $`{\displaystyle \frac{1}{2}}^2\mathrm{\Psi }`$ (10)
$`W(\mathrm{\Psi })`$ $``$ $`\left[V(𝐱)a|\mathrm{\Psi }|^2+\mu \right]\mathrm{\Psi }`$ (11)
The stationary states of (9) satisfy:

$$0=L\mathrm{\Psi }+W(\mathrm{\Psi })$$
(12)
Without loss of generality, $`\mathrm{\Psi }`$ can be chosen to be real. Our objective is to calculate the eigenpairs $`(\lambda ,\psi )`$ of the operator that results from linearizing (9) about a stationary state $`\mathrm{\Psi }`$. (We use $`\mathrm{\Psi }`$ to designate solutions to the nonlinear problem (12) and $`\psi `$ to designate eigenvectors, which are solutions to the linear problem to be defined below.) In order to correctly formulate the linear stability problem, it is necessary to first decompose $`\psi =\psi ^R+i\psi ^I`$. We write the linearized evolution equation

$$\frac{\psi }{t}=i(L+DW(\mathrm{\Psi }))(\psi ^R+i\psi ^I)$$
(13)
where $`DW(\mathrm{\Psi })`$ is the Fréchet derivative, or Jacobian, of $`W`$ evaluated at $`\mathrm{\Psi }`$. $`DW(\mathrm{\Psi })`$ acts on $`\psi `$ via:

$$DW\psi =DW^R\psi ^R+iDW^I\psi ^I$$
(14)
where we have omitted the functional dependence of $`DW`$, $`DW^R`$, and $`DW^I`$ on $`\mathrm{\Psi }`$, and where

$`DW^R`$ $``$ $`\mu V(𝐱)3a\mathrm{\Psi }^2`$ (15a)
$`DW^I`$ $``$ $`\mu V(𝐱)a\mathrm{\Psi }^2.`$ (15b)
Equation (13) is then written in matrix form as

$$\frac{}{t}\left(\begin{array}{c}\psi ^R\\ \psi ^I\end{array}\right)=\left[\begin{array}{cc}0& (L+DW^I)\\ L+DW^R& 0\end{array}\right]\left(\begin{array}{c}\psi ^R\\ \psi ^I\end{array}\right)$$
(16)
The eigenmodes $`(\lambda ,\psi ^R,\psi ^I)`$ satisfy

$$\lambda \left(\begin{array}{c}\psi ^R\\ \psi ^I\end{array}\right)=\left[\begin{array}{cc}0& (L+DW^I)\\ L+DW^R& 0\end{array}\right]\left(\begin{array}{c}\psi ^R\\ \psi ^I\end{array}\right)$$
(17)
Note that this eigensystem is usually presented in the literature (see, for example, reference Griffin ) in terms of the variables

$$(\omega ^\mathrm{B},\psi ,\psi ^{})(i\lambda ,\psi ^R+i\psi ^I,\psi ^Ri\psi ^I)$$
(18)
as the equivalent Bogoliubov–de Gennes coupled equations

$$\omega ^\mathrm{B}\left(\begin{array}{c}\psi \\ \psi ^{}\end{array}\right)=\left[\begin{array}{cc}L+DW^\mathrm{B}& a\mathrm{\Psi }^2\\ a\mathrm{\Psi }^2& (LDW^\mathrm{B})\end{array}\right]\left(\begin{array}{c}\psi \\ \psi ^{}\end{array}\right)$$
(19)
where

$$DW^\mathrm{B}\mu V(𝐱)2a\mathrm{\Psi }^2$$
(20)
In the following, we will work with matrix formulation (17) because it avoids a potential notational inconsistency of (19) arising from the fact that $`\psi `$ and $`\psi ^{}`$ are complex conjugates only when $`\omega ^\mathrm{B}`$ is imaginary.

It is more convenient to work with the square of the matrix in (17):

$`\lambda ^2\left(\begin{array}{c}\psi ^R\\ \psi ^I\end{array}\right)=`$ (23)
$`\left[\begin{array}{cc}\left(L+DW^I\right)\left(L+DW^R\right)& 0\\ 0& \left(L+DW^R\right)\left(L+DW^I\right)\end{array}\right]\left(\begin{array}{c}\psi ^R\\ \psi ^I\end{array}\right)`$ (28)
Because (23) is block diagonal, it can be separated into the two problems:

$`\lambda ^2\psi ^R=(L+DW^I)(L+DW^R)\psi ^R`$ (29a)
$`\lambda ^2\psi ^I=(L+DW^R)(L+DW^I)\psi ^I`$ (29b)
Problems (29a) and (29b) are closely related, Since the operators $`L`$, $`DW^I`$, and $`DW^R`$ are all self-adjoint under the standard Euclidean inner product, the operators in (29a) and (29b) are adjoint to each other. If $`\psi ^R`$ is an eigenvector of (29a) with non-zero eigenvalue $`\lambda ^2`$, then $`(L+DW^R)\psi ^R`$ is an eigenvector of (29b) with the same eigenvalue. (Similarly, if $`(\lambda ,\psi ^R,\psi ^I)`$ is an eigenmode of (17), then $`(\lambda ,\psi ^R,\psi ^I)`$ is also an eigenmode of (17).) Thus, we solve only (29a). The eigenvalues $`\lambda ^2`$ of (29a)-(29b) must be either complex conjugate pairs or real. We find them to be real and (almost all) negative, perturbed only slightly from the eigenvalues of $`L^2`$. The eigenvalues $`\lambda `$ of (17) are therefore found to be either pure imaginary or pure real, with most imaginary.

Problems (17) and (29a)-(29b) have neutral eigenmodes which reflect the physical invariances of the problem. Since $`DW^I\mathrm{\Psi }=W(\mathrm{\Psi })`$, then the stationary state $`\mathrm{\Psi }`$ is a neutral mode of $`L+DW^I`$ and hence of problem (29b). This neutral mode is the phase mode of (II.1), since its existence is a consequence of the invariance of solutions $`\mathrm{\Psi }`$ to (12) under multiplication by any complex number on the unit circle. The corresponding eigenmode of problem (29a) is $`d\mathrm{\Psi }/d\mu `$. This neutral mode can be understood as a consequence of differentiating (II.1) with respect to $`\mu `$:

$$0=\frac{d}{d\mu }[(L+W)\mathrm{\Psi }]=(L+DW^R)\frac{d\mathrm{\Psi }}{d\mu }+\mathrm{\Psi }$$
(30)
Thus,

$$(L+DW^I)(L+DW^R)\frac{d\mathrm{\Psi }}{d\mu }=(L+DW^I)\mathrm{\Psi }=0.$$
(31)
In terms of the original problem (17), the phase mode $`(\lambda ,\psi ^R,\psi ^I)=(0,0,\mathrm{\Psi })`$ is a neutral eigenvector, while (30) shows that $`(\lambda ,\psi ^R,\psi ^I)=(0,d\mathrm{\Psi }/d\mu ,0)`$ is a neutral generalized eigenvector, the two modes forming a Jordan block for (17).

In practice, we fix $`\mu `$ to calculate the stationary states $`\mathrm{\Psi }`$ and the eigenvalues. The operators of (29a)-(29b) depend on $`\mu `$ both explicitly and through $`\mathrm{\Psi }`$. For $`\mu `$ above a critical value $`\mu _c`$ all eigenvalues $`\lambda `$ are imaginary, i.e. $`\mathrm{\Psi }`$ is an elliptic stationary state of (9). As $`\mu `$ crosses $`\mu _c`$, we will see that one imaginary pair fuses at zero, and becomes real, with one positive and one negative value of $`\lambda `$ for $`\mu <\mu _c`$. Stationary states for $`\mu <\mu _c`$ are thus hyperbolic in the directions corresponding to these eigenvalues.

## III Bifurcation and Stability of Condensates

In this section we will find the stationary solutions and study the stability of isotropic ($`\omega _r=\omega _z`$), cigar-like ($`\omega _r/5=\omega _z`$) and pancake-like ($`\omega _r=\omega _z/5`$) condensates. These results were obtained by solving equation (12) for the stationary states and (17) or (29a) for the corresponding bifurcating eigenvalues. The system is discretized using pseudo-spectral methods in a spherical domain for the isotropic case and in a periodic Cartesian domain for the non-isotropic cases. We use Newton’s method to calculate the branches of stationary states. The bifurcating eigenvalue is found in the isotropic case by diagonalizing the matrix corresponding to (17). In the non-isotropic case, we use instead the iterative inverse Arnoldi method, which requires only actions of the operator in (29a). The BiCGSTAB variant of the conjugate gradient method is used to solve the linear systems required by both Newton’s method and the inverse Arnoldi method. The numerical methods we use are described in greater detail in the Appendix.

### III.1 Isotropic Condensate

In order to compare our results with the existing experiments on quasi-isotropic condensates, we will use the following physical constants, corresponding to 7Li atoms in a radial trap: $`m=1.16\times 10^{26}\mathrm{kg}`$, $`\stackrel{~}{a}=27.3a_0`$ (with $`a_0`$ the Bohr radius) and $`\widehat{\omega }=(\stackrel{~}{\omega }_x\stackrel{~}{\omega }_y\stackrel{~}{\omega _z})^{1/3}=908.41\mathrm{s}^1`$. These values yield $`a=5.74\times 10^3`$. With these parameters, the mean-field approximation (II.1) is expected to be very reliable. Note that we ignore the contributions of non-condensed atoms. They interact with the condensate only through a nearly constant background density term, inducing no significant change in the dynamics of the system Hourbiers .

The values of the energy functional $``$ and the (smallest absolute value) square eigenvalue $`\lambda ^2`$ versus particle number $`𝒩`$ are shown as solid lines on figure 2 (top and bottom, respectively). The eigenvalues are imaginary on the metastable elliptic lower branch ($`\lambda ^2<0`$) and real on the unstable hyperbolic upper branch ($`\lambda ^2>0`$). Using (2) on stationary solutions we obtain $`d/d𝒩=\mu `$. Thus $`\mu `$ is the slope of $``$ and the lower branches $`_{}`$, $`\lambda _{}^2`$ (respectively upper branches $`_+`$, $`\lambda _+^2`$) are scanned for $`\mu >\mu _c`$ (respectively $`\mu <\mu _c`$). The point $`\mu =\mu _c`$ determines the maximum number of particles $`𝒩=𝒩_c`$ for which stationary solutions exist. We have checked that all the other pairs of eigenvalues are imaginary on both branches (data not shown).

The dashed curves on figure 2 are derived from the Gaussian variational approximation which will be defined in section III.3 for the general anisotropic case. In the present isotropic case, this approximation can be solved in closed form, yielding the expressions PRLiso :

$`𝒩(\mu )={\displaystyle \frac{4\sqrt{2\pi ^3}\left(8\mu +3\sqrt{7+4\mu ^2}\right)}{7|\stackrel{~}{a}|\left(2\mu +\sqrt{7+4\mu ^2}\right)^{3/2}}},`$ (32)
$`=𝒩(\mu )\left(\mu +3\sqrt{7+4\mu ^2}\right)/7.`$ (33)
The number of particles $`𝒩`$ is maximal at $`𝒩_c^G=8\sqrt{2\pi ^3}/|5^{5/4}a|`$ for $`\mu =\mu _c^G=1/2\sqrt{5}`$. The eigenvalues can also be obtained in closed form from the linearized equations of motion PRLiso :

$$\lambda ^2(\mu )=8\mu ^24\mu \sqrt{7+4\mu ^2}+2$$
(34)
By inspection of figure 2 it is apparent that both the solution of the GP equation and the Gaussian variational approximation share the same qualitative behavior, with quantitative discrepancies. Figure 3 shows the physical origin of the quantitative errors in the Gaussian approximation. It is apparent that the exact solution is well approximated by a Gaussian only for small $`𝒩`$ on the stable (elliptic) branch.

### III.2 Hamiltonian Saddle Node Normal Form

The qualitative behavior displayed on figure 2 by the physical quantities $``$ and $`\lambda ^2`$ near the critical value $`𝒩=𝒩_c`$ is the generic signature of a HSN bifurcation defined, at lowest order, by the normal form GH ; Elphick

$$\ddot{q}=\delta \beta q^2,$$
(35)
where $`\delta =(1𝒩/𝒩_c)`$ is the bifurcation parameter, $`\beta `$ is a dimensionless constant, and $`q`$ is the coordinate describing the state of the system in the direction of the phase space that becomes unstable. Indeed, introducing the additional non-dimensional quantities $``$ and $`\gamma `$ to define the appropriate energy

$$=_0+\frac{1}{2}\dot{q}^2\delta q+\frac{1}{3}\beta q^3\gamma \delta ,$$
(36)
it is straightforward to derive from (35) that, close to the critical point $`\delta =0`$, the universal scaling laws are given by

$`_\pm `$ $`=`$ $`_c_l\delta \pm _\mathrm{\Delta }\delta ^{3/2},`$ (37)
$`\lambda _\pm ^2`$ $`=`$ $`\pm \lambda _\mathrm{\Delta }^2\delta ^{1/2},`$ (38)
where $`_c=_0`$, $`_l=\gamma `$, $`_\mathrm{\Delta }=2/3\sqrt{\beta }`$ and $`\lambda _\mathrm{\Delta }^2=2\sqrt{\beta }`$. Note that these relations can be inverted to obtain the parameters in (35) from the critical data. For the Gaussian approximation the critical amplitudes can be computed from equations (32) and (33). One finds

$`_c`$ $`=`$ $`{\displaystyle \frac{4\sqrt{2\pi ^3}}{5^{3/4}|a|}}`$ (39)
$`_\mathrm{\Delta }`$ $`=`$ $`{\displaystyle \frac{64\sqrt{\pi ^3}}{5^{9/4}|a|}}`$ (40)
$`\lambda _\mathrm{\Delta }^2`$ $`=`$ $`4\sqrt{10}.`$ (41)
For the exact solutions, we obtain the critical amplitudes $`_\mathrm{\Delta }=1340`$ and $`\lambda _\mathrm{\Delta }^2=14.68`$ by performing fits on the numerical data. Comparing both results, we find that the Gaussian approximation captures the bifurcation qualitatively, but with quantitative errors of $`17\%`$ for $`𝒩_c`$ Ruprecht , $`24\%`$ for $`_\mathrm{\Delta }`$ and $`14\%`$ for $`\lambda _\mathrm{\Delta }^2`$ in the isotropic case PRLiso .

The phase portrait of the normal form is shown on Figure 4. When $`\delta =(1𝒩/𝒩_c)>0`$, equation (35) admits two fixed points $`Q_\pm =\sqrt{\delta /\beta }`$, as shown in Fig. 4A. Thus, a hyperbolic stationary state and an elliptic stationary state coexist. The phase space is separated into two regions by a separatrix which is a homoclinic orbit linking the hyperbolic stationary state to itself. Trajectories inside the orbit remain bounded near the elliptic fixed point. If the condensate is taken beyond the separatrix by a perturbation (e.g. thermal excitations or quantum tunneling, see below section IV), it will fall into unbounded (hyperbolic) trajectories and collapse. As $`𝒩`$ is increased, the hyperbolic and elliptic stationary states approach one another, (Fig. 4B) and the homoclinic orbit inside which orbits are bounded is reduced. The two stationary states join at $`𝒩=𝒩_c`$ (Fig. 4C), at which the HSN occurs. No stationary state exists for $`𝒩>𝒩_c`$.

### III.3 Non-isotropic Condensates

We now briefly present the main expressions obtained from a Gaussian variational analysis of the G-P equation with a cylindrical potential trap. Some of these results have been previously obtained by other authors GammalPRA64 ; GammalPreprint ; PerezGarciaAsym1 ; PerezGarciaAsym2 ; AbdullaevPRA63 ; AdhikariPRE65 . We therefore restrict our discussion to the equations that will be used in our analysis of the condensate lifetimes.

The trial function is a Gaussian solution to the linear $`(a=0)`$ Schrödinger equation in which we incorporate eight variational parameters in order to take into account the anisotropy of the system. The form of the ansatz is given by:

$`\mathrm{\Psi }(x,y,z,t)=\left(A_r(t)+iA_i(t)\right)`$
$`\times \mathrm{exp}\{({\displaystyle \frac{1}{X(t)^2}}+i\varphi _X(t)){\displaystyle \frac{x^2}{2}}`$
$`\left({\displaystyle \frac{1}{Y(t)^2}}+i\varphi _Y(t)\right){\displaystyle \frac{y^2}{2}}`$
$`({\displaystyle \frac{1}{Z(t)^2}}+i\varphi _Z(t)){\displaystyle \frac{z^2}{2}}\}`$ (42)
where the real parameters $`\{A_r,A_i\}`$, $`\{\varphi _X,\varphi _Y,\varphi _Z\}`$ and $`\{X,Y,Z\}`$ are related to the amplitude, the phase and the width of the Gaussian profile respectively. The Euler-Lagrange equations associated with the trial function (III.3) and the action defined in equation (1) can be reduced to the following system of second-order differential equations:

$`{\displaystyle \frac{d^2X}{dt^2}}`$ $`=`$ $`\omega _x^2X{\displaystyle \frac{\nu }{X^2YZ}}+{\displaystyle \frac{1}{X^3}}`$
$`{\displaystyle \frac{d^2Y}{dt^2}}`$ $`=`$ $`\omega _y^2Y{\displaystyle \frac{\nu }{XY^2Z}}+{\displaystyle \frac{1}{Y^3}}`$
$`{\displaystyle \frac{d^2Z}{dt^2}}`$ $`=`$ $`\omega _z^2Z{\displaystyle \frac{\nu }{Z^2XY}}+{\displaystyle \frac{1}{Z^3}}`$ (43)
where

$$\nu =\sqrt{\frac{2}{\pi }\frac{m(\omega _x\omega _y\omega _z)^{1/3}}{\mathrm{}}}|\stackrel{~}{a}|𝒩.$$
(44)
The evolution of the condensate is better understood by drawing an analogy between its width and the motion of a particle with coordinates $`(X,Y,Z)`$ moving in the potential

$`U(X,Y,Z)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(\omega _x^2X^2+\omega _y^2Y^2+\omega _z^2Z^2\right){\displaystyle \frac{\nu }{XYZ}}`$ (45)
$`+`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{X^2}}+{\displaystyle \frac{1}{Y^2}}+{\displaystyle \frac{1}{Z^2}}\right).`$
Indeed, defining $`P_x=dX/dt`$, $`P_y=dY/dt`$, $`P_z=dZ/dt`$ and the Hamiltonian

$`H(P_x,P_y,P_z,X,Y,Z)={\displaystyle \frac{1}{2}}\left(P_x^2+P_y^2+P_z^2\right)+U(X,Y,Z),`$
we find that equations (III.3) transform into Hamilton’s equations of motion.

If we consider now a potential trap (5) with cylindrical symmetry ($`\omega _r\omega _x=\omega _y`$) equations (III.3) can be simplified by using $`X(t)=Y(t)`$. We thus find that (III.3) yields two fixed points $`(X_+,Z_+)`$ and $`(X_{},Z_{})`$ which describe the stationary solutions for $`\mathrm{\Psi }(x,y,z,t)`$. These obey

$`0`$ $`=`$ $`\omega _r^2{\displaystyle \frac{4\mu }{7X_\pm ^2}}{\displaystyle \frac{5}{7X_\pm ^4}}+{\displaystyle \frac{2}{7}}\left({\displaystyle \frac{1}{X_\pm ^2}}+{\displaystyle \frac{1}{Z_\pm ^2}}\right){\displaystyle \frac{1}{X_\pm ^2}}`$
$`0`$ $`=`$ $`\omega _z^2{\displaystyle \frac{4\mu }{7Z_\pm ^2}}{\displaystyle \frac{5}{7Z_\pm ^4}}+{\displaystyle \frac{4}{7}}\left({\displaystyle \frac{1}{X_\pm ^2Z_\pm ^2}}\right),`$ (46)
where the chemical potential $`\mu `$ is related to the total number of particles through

$`𝒩`$ $`=`$ $`{\displaystyle \frac{2L_0}{7|\stackrel{~}{a}|}}\sqrt{2\pi }X^2Z\left({\displaystyle \frac{1}{X^2}}+{\displaystyle \frac{1}{2Z^2}}\mu \right).`$ (47)
The fixed points correspond to a metastable center $`(X_+,Z_+)`$ and to an unstable saddle point $`(X_{},Z_{})`$, respectively. They are analogous to the $`Q_+`$ and $`Q_{}`$ points appearing in the phase portraits on Figure 4. The solutions to (III.3) can be computed numerically, together with the linearized variational equations evaluated at every stationary point.

Figures 5 and 6 show $``$ and $`\lambda ^2`$ for the cigar and pancake cases, respectively. The solid lines present the values obtained by discretizing and solving numerically the original differential equations (12) and (29a), using the methods described in the Appendix. The dashed lines were computed using the Gaussian approximation described above. Both the isotropic and non-isotropic cases display saddle-node bifurcations. This is to be expected, since the saddle-node bifurcation is the generic way in which stable and unstable branches meet GH .

It is apparent from figures 2, 5, and 6 that the exact critical number of particles $`𝒩_c^E`$ is smaller than the Gaussian value $`𝒩_c^G`$ for all three geometries Ueda ; Ruprecht ; GammalPRA64 ; GammalPreprint . Table 1 compares the different critical $`𝒩`$ values obtained.

In order to compare properly the HSN bifurcations obtained for the three aspect ratios studied, we can rescale the intensity of the potential to obtain the same $`𝒩_c^E`$ for all cases. In general, any confining harmonic potential with frequencies $`\omega _r`$ and $`\omega _z`$ that produces a critical number of particles $`𝒩_c`$ can be rescaled by a factor

$$c=\left(\frac{𝒩_c}{𝒩_{c}^{}{}_{}{}^{}}\right)^2,$$
(48)
to obtain the a new potential with frequencies $`\omega _r^{}=c\omega _r`$ and $`\omega _z^{}=c\omega _z`$, which will have the critical number of particles $`𝒩_c^{}`$. The remaining physical quantities for the new potential are obtained through the following transformations:

$`\mathrm{\Psi }^{}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Psi }}{c^{1/4}}}`$ (49)
$`𝒩^{}`$ $`=`$ $`{\displaystyle \frac{𝒩}{\sqrt{c}}}`$ (50)
$`^{}`$ $`=`$ $`{\displaystyle \frac{}{\sqrt{c}}}`$ (51)
$`\lambda ^{}`$ $`=`$ $`\lambda .`$ (52)
We choose arbitrarily to rescale the potential intensity so that all $`𝒩_c^E`$ are equal to that for the isotropic case $`𝒩_c^{\mathrm{iso}}`$. Table 2 shows the value of the rescaling factors $`c_{\mathrm{cig}}`$ and $`c_{\mathrm{pan}}`$ for the cigar and pancake cases respectively, as obtained from equation (48) using $`𝒩_c^{}=𝒩_c^{\mathrm{iso}}`$. The last two columns of this table show the critical amplitudes obtained for the rescaled $``$ and $`\lambda ^2`$ curves. These were obtained by fitting the HSN asymptotic forms given in relations (37) and (38) to the rescaled data.

## IV Lifetime of Condensates

In this section we first find expressions for the TIC, MQT and ICO decay rates. Using the numerical data presented in the previous section, we then compute these decay rates for the $`\omega _r=\omega _z`$ (isotropic), $`\omega _r/5=\omega _z`$ (cigar) and $`\omega _r=\omega _z/5`$ (pancake) cases. Finally, we compare the results obtained for these three potential geometries by studying the spontaneous isotropization of the condensates.

### IV.1 Definition and Computations of Decay Rates

The TIC (thermally induced collapse) decay rate $`\mathrm{\Gamma }_T`$ is estimated using the formula Gardiner

$$\frac{\mathrm{\Gamma }_T}{\widehat{\omega }}=\frac{|\lambda _+|}{2\pi }\mathrm{exp}\left[\frac{\mathrm{}\widehat{\omega }}{k_BT}\left(_+_{}\right)\right]$$
(53)
where $`\mathrm{}\widehat{\omega }(_+_{})`$ is the (dimensionalized) height of the nucleation energy barrier (with $`\widehat{\omega }`$ the reference frequency introduced in Section II.1), $`T`$ is the temperature of the condensate and $`k_B`$ is the Boltzmann constant. Note that the prefactor characterizes the typical decay time which is controlled by the slowest part of the nucleation dynamics: the top-of-the-barrier saddle point eigenvalue $`\lambda _+`$ and not $`\lambda _{}`$ as used in Stoof1 . However, near the bifurcation both eigenvalues scale in the same way and the behavior of $`\mathrm{\Gamma }_T`$ can be obtained directly from the universal saddle-node scaling laws (37) and (38). Thus the exponential factor and the prefactor vanish respectively as $`\delta ^{3/2}`$ and $`\delta ^{1/4}`$.

We estimate the MQT (macroscopic quantum tunneling) decay rate using an instanton technique that takes into account the semi classical trajectory giving the dominant contribution to the quantum action path integral Stoof1 ; Ueda . We approximate this so-called bounce trajectory by the solution of the equation of motion

$$\frac{d^2q(t)}{dt^2}=\frac{d\stackrel{~}{\mathrm{\Phi }}(q)}{dq}$$
(54)
starting and ending at the fixed point $`q_f`$ of the phase space where $`(q_f)=_{}`$. The Euclidean potential $`\stackrel{~}{\mathrm{\Phi }}(q)`$ is defined so that $`\stackrel{~}{\mathrm{\Phi }}(q)`$ reconstructs the Hamiltonian dynamics in the region scanned by the bounce trajectory (see Figure 7). We represent it by a fourth-order polynomial of the form

$$\stackrel{~}{\mathrm{\Phi }}(q)=\alpha _0+\alpha _2q^2+\alpha _3q^3+\alpha _4q^4$$
(55)
coefficients $`\alpha _0`$, $`\alpha _2`$, $`\alpha _3`$ and $`\alpha _4`$ chosen such that

$`\stackrel{~}{\mathrm{\Phi }}(0)`$ $`=`$ $`_+`$ (56a)
$`\stackrel{~}{\mathrm{\Phi }}(q_f)`$ $`=`$ $`_{}`$ (56b)
$`_q^2\stackrel{~}{\mathrm{\Phi }}(0)`$ $`=`$ $`\lambda _+(𝒩)`$ (56c)
$`_q^2\stackrel{~}{\mathrm{\Phi }}(q_f)`$ $`=`$ $`\lambda _{}(𝒩)`$ (56d)
We thus obtain a semi-analytic polynomial expression for $`\stackrel{~}{\mathrm{\Phi }}(q)`$ where the coefficients are determined through the numerical values presented in figures 2, 5, and 6.

Once $`\stackrel{~}{\mathrm{\Phi }}(q)`$ and the bounce point $`q_b`$ (defined through the relation $`\stackrel{~}{\mathrm{\Phi }}(q_b)=\stackrel{~}{\mathrm{\Phi }}(q_f)`$) are known, the MQT rate is estimated as

$$\frac{\mathrm{\Gamma }_Q}{\widehat{\omega }}=\sqrt{\frac{|\lambda _{}|v_0^2}{4\pi }}\mathrm{exp}\left[\frac{4}{\sqrt{2}}_{q_f}^{q_b}\sqrt{\stackrel{~}{\mathrm{\Phi }}(q)\stackrel{~}{\mathrm{\Phi }}(q_f)}𝑑q\right],$$
(57)
where $`v_0`$ is defined by the asymptotic form of the bounce trajectory $`q(t)`$ as it approaches $`q_f`$ Stoof1 , given by $`q(\tau )q_f+(v_0/|\lambda _{}|)\mathrm{exp}[|\lambda _{}\tau |]`$.

In the same way as was done for the TIC, universal scaling laws can be derived close to criticality from (35), (37) and (38). We find that the exponential factor in (57) follows the same scaling as $`\sqrt{|_+_{}|}dq`$. It therefore vanishes as $`\sqrt{\delta ^{3/2}}\delta ^{1/2}=\delta ^{5/4}`$. The asymptotic form of $`q(t)`$ shows that $`dq`$ follows the same law as $`v_0/|\lambda _{}|`$. Thus $`v_0\delta ^{3/4}`$, and the prefactor vanishes as $`\sqrt{\delta ^{1/4}}\delta ^{3/4}=\delta ^{7/8}`$. Note that these universal scaling laws agree with those already derived in the Gaussian case in Ueda .

The TIC (53) and MQT (57) decay rates obtained for the exact and Gaussian stationary states are shown in Figure 8. To validate these results we checked that the Gaussian TIC decay rates computed in Sackett  are found when we (incorrectly at a finite distance from criticality) replace $`\lambda _+`$ by $`\lambda _{}`$ in equation (53) (data not shown). We also checked that our Gaussian MQT decay rate agrees with the one previously computed in Ueda .

The ICO (inelastic two and three body collision) atomic decay rates are evaluated using the formula $`d𝒩/dt=f_C(𝒩)`$ with

$$f_C(𝒩)=K|\mathrm{\Psi }|^4d^3𝐱+L|\mathrm{\Psi }|^6d^3𝐱,$$
(58)
where $`K=3.8\times 10^4[s^1]`$ and $`L=2.6\times 10^7[s^1]`$ as in Shi ; Dodd . In order to compare the particle decay rate (58) to the condensate collective decay rates obtained for the TIC and MQT, we compute the condensate ICO half-life using

$$\tau _{1/2}(𝒩)=_{𝒩/2}^𝒩\frac{dn}{f_C(n)}$$
(59)
and plot $`\tau _{1/2}^1`$ on figure 8.

Figures 8 and 9 compare the condensate decay rates for the isotropic and the cigar and pancake potentials, rescaled by $`c_{\mathrm{cig}}`$ and $`c_{\mathrm{pan}}`$ as described in Section III.3. We note that the three aspect ratios generate very similar results after rescaling. The relative magnitudes of the different decay rates – TIC, MQT, and ICO – are the same for the three cases. At $`T1[\mathrm{nK}]`$ the MQT effect becomes important compared to the ICO decay in a region very close to $`𝒩_c^E`$ ($`\delta 8\times 10^3`$). This was shown in Ueda  using Gaussian computations but evaluating them with the exact maximal number of condensed particles $`𝒩_c^E`$. Figure 9 shows that even for temperatures as low as $`2[\mathrm{nK}]`$, the TIC decay rate exceeds the MQT rate except in a region extremely close to $`𝒩_c`$ ($`\delta <5\times 10^3`$), where the condensates will live less than $`10^1[s]`$. Thus, in the experimental case of 7Li atoms, the relevant effects are ICO and TIC, with the crossover that is shown on Figure 9.

### IV.2 Spontaneous Isotropization of Condensates

The decay rates of the isotropic and non-isotropic cases shown in figure 8 are quite similar, despite the fact that $`\omega _z`$ and $`\omega _r`$ differ by a factor of five. We have investigated this question by examining the wave functions $`\mathrm{\Psi }`$ for the pancake and cigar cases. These wave functions are peaked at the origin, as shown on figure 3. Their characteristic length scales in the axial and radial directions, $`\mathrm{}_z`$ and $`\mathrm{}_r`$, can be measured by computing the ratios of the value of $`\mathrm{\Psi }`$ to its curvature at the origin.

More specifically, we define

$`\mathrm{}_z^2`$ $``$ $`\mathrm{\Psi }\left({\displaystyle \frac{^2\mathrm{\Psi }}{z^2}}\right)^1(r=0,z=0)`$ (60a)
$`\mathrm{}_r^2`$ $``$ $`\mathrm{\Psi }\left({\displaystyle \frac{1}{r}}{\displaystyle \frac{}{r}}r{\displaystyle \frac{\mathrm{\Psi }}{r}}\right)^1(r=0,z=0)`$ (60b)
We then obtain the ellipticity of the wavefunction as the ratio $`\mathrm{}`$ of these length scales: $`\mathrm{}=\mathrm{}_r/\mathrm{}_z`$ for the cigar and $`\mathrm{}=\mathrm{}_z/\mathrm{}_r`$ for the pancake. These ellipticity ratios are shown in figure 10 as a function of $`\mu `$. For large $`\mu `$, i.e. away from the saddle-node bifurcation along the stable branch, $`\mathrm{}`$ decreases rapidly away from one, indicating that the wave function is highly non-isotropic. At the saddle-node bifurcation, $`\mathrm{}=0.89`$ at $`\mu =0.38`$ for the cigar and $`\mathrm{}=0.80`$ at $`\mu =0.31`$ for the pancake. As $`\mu `$ is decreased, i.e. as we leave the saddle-node bifurcation along the unstable branch, $`\mathrm{}`$ approaches one as the wave function becomes more nearly spherically symmetric. This trend is present both in the numerical solution to the G-P equation and in the Gaussian approximation, as can be seen on figure 10. Since the decay rates result from the scaling behavior near the saddle-node bifurcation, where the condensate is fairly isotropic, it follows that the decay rates are similar for the cigar, pancake, and spherically symmetric geometries, as we have shown on Figures 8 and 9.

The spontaneous isotropization of condensates when $`\mu `$ is decreased can be understood by the following phenomenological reasoning. When $`\mu `$ grows, the balance of terms in the right hand side of equation (II.1) changes. For small $`\mu `$, it is dominated by the the isotropic $`^2`$ and the anisotropic $`V(x)`$ terms. But for large $`\mu `$, the wavefunction $`\mathrm{\Psi }`$ is strongly peaked and the $`^2`$ and nonlinear terms, both isotropic, become dominant.

Figure 10 also provides a test of our numerical spatial resolution. By computing the ellipticity $`\mathrm{}`$ for different Fourier truncation levels, we show that $`\mathrm{}`$ changes with the resolution for low $`\mu `$, especially for the pancake case, where we used fewer Fourier modes than in the other calculations. Note however that our decay rate calculations only use results near the saddle-node bifurcation, where $`\mathrm{}`$ varies by less than 0.1% when different truncation levels are used. This indicates that $`\mathrm{\Psi }`$ was adequately resolved in the region of interest.

As $`\mu `$ is decreased, the wave functions become more highly peaked for both our numerical results and for the Gaussian approximation (see figures 3 and 11). This is the main reason for the declining accuracy. To continue the computations further, the size of the periodic box should be reduced along with $`\mu `$. We believe that, with adequate resolution, all of the exact wave functions would become spherically symmetric as $`\mu `$ decreases, as do the Gaussian approximations.

## V Conclusion

We have demonstrated that it is possible to numerically compute the stationary states, the bifurcating eigenvalues and the lifetime of anisotropic attractive Bose Einstein condensates.

The Gaussian mean field approximation was found to have significant quantitative errors for all the different confining potential geometries that were studied, when compared with numerical solutions to the G-P equation.

Spontaneous isotropization of the metastable condensate was found to occur as the critical number of particles is approached, yielding a life-time that depends weakly on the anisotropy of the confining potential.

Direct methods – Gaussian elimination and diagonalization – were used in treating the spherically symmetric case, of size $`128`$, but are far too costly for the three-dimensional case, of size $`10^6`$. In fact, since we only calculated axisymmetric stationary states and eigenvectors with an additional midplane symmetry, an intermediate two-dimensional axisymmetric cylindrical representation could have been implemented, of size 5000, permitting the use of direct methods. Our purpose, however, has been to construct and explore numerical methods appropriate for a general non-isotropic case.

The methods used to compute stationary states and bifurcating eigenvalues for the non-isotropic cases where essentially analogous. Each consists of a powerful and rapid outer iteration: Newton’s method for the stationary states and the inverse Arnoldi method for the eigenvalues. The large linear systems that need to be inverted within each method are solved by the same inner biconjugate gradient iteration – BiCGSTAB – and constitutes the main numerical difficulty. Its convergence is greatly improved by an inverse Laplacian preconditioning which is empirically tuned by adjusting the pseudo-timestep $`\sigma `$ in Newton’s method or the shift $`s`$ in Arnoldi’s method.

Our results and implementation have demonstrated that all these numerical techniques can be successfully combined to calculate the stationary states and eigenvectors for the G-P equation in a confining potential with an arbitrary three-dimensional geometry.

Acknowledgments: This work was supported by ECOS-CONICYT program no. C01E08 and by NSF grant DMR-0094569. Computations were performed at the Institut du Développement et des Ressources en Informatique Scientifique (IDRIS) of the CNRS. *

## Appendix A Numerical Methods

### A.1 Spatial Discretization

The operators $`L`$ and $`W`$ defined in (10) and (11) are spatially discretized using the pseudospectral method Got-Ors . For the isotropic case, the spherically symmetric $`\mathrm{\Psi }(r,t)`$ is expanded as a series of even Chebyshev polynomials $`T_{2n}(r/R)`$, on which the boundary condition $`\mathrm{\Psi }(R,t)=0`$ is imposed. The domain is taken to be $`0rR=4`$ and the resolution used is $`N_R=128`$. For the non-isotropic cases, we use a three-dimensional periodic Cartesian domain and $`\mathrm{\Psi }`$ is expanded as a three-dimensional trigonometric (Fourier) series. The cigar case is solved in a periodic domain of size $`(L_x,L_y,L_z)=(5.39,5.39,12.04)`$ in units of $`L_0`$, using $`(N_x,N_y,N_z)=(96,96,96)`$ gridpoints or trigonometric modes (with a $`2/3`$ dealiasing rule,) so the total number of gridpoints or trigonometric functions is as high as $`N_{3D}=10^6`$. (The more poorly resolved pancake case was calculated using $`(L_x,L_y,L_z)=(12.04,12.04,5.39)`$ and $`(N_x,N_y,N_z)=(48,48,96)`$.) The harmonic potential (5) is approximated by a periodic potential by writing $`x=\mathrm{arcsin}(\mathrm{sin}(x))`$ and Taylor-expanding the $`\mathrm{arcsin}`$ function. This leads to a Fourier series for the potential, which is truncated according to the resolution used.

Pseudospectral methods require performing over $`\mathrm{\Psi }`$, at every iteration, a Chebyshev transform in the isotropic case or a Fourier transform in the non-isotropic case. These operations consume a time proportional to $`N_R\mathrm{log}N_R`$ or $`N_{3D}\mathrm{log}(N_{3D})`$, respectively. Actions and inversions of the Laplacian $`L`$ are carried out on the Chebyshev or Fourier representations of $`\mathrm{\Psi }`$, while actions of the multiplicative operator $`W`$ are carried out on its grid representations. The time required by these operations scales approximately linearly in $`N_R`$ or $`N_{3D}`$.

### A.2 Stationary States

As stated in section II.2, the stationary states of (II.1) that correspond to minima of $``$ at a given value of $`𝒩`$ can be obtained by integrating to relaxation the diffusion equation

$$\frac{\mathrm{\Psi }}{t}=L\mathrm{\Psi }+W(\mathrm{\Psi })$$
(61)
where the initial data $`\mathrm{\Psi }(t=0)`$ has a total number of particles $`𝒩`$ and the value of the Lagrange multiplier $`\mu `$ is fixed during the relaxation by the condition $`𝒩/t=0`$.

To integrate (61) a mixed implicit-explicit first-order time-stepping scheme is used:

$$\mathrm{\Psi }(t+\sigma )=(I\sigma L)^1(I+\sigma W)\mathrm{\Psi }(t)$$
(62)
where $`I`$ is the identity operator. The Helmholtz operator $`(I\sigma L)^1`$ is easily inverted in the Chebyshev or Fourier representation. The motivation for integrating $`L`$ implicitly is to avoid the extremely small timesteps that would otherwise be necessitated by the wide range of eigenvalues of the Laplacian.

This relaxation method is equivalent to that used in Ruprecht  and can only reach the stable stationary solutions of (61). In order to also capture unstable stationary solutions Seydel  we implemented a Newton branch-following algorithm PRLiso ; TuckBark . We search for fixed points of (62), a condition strictly equivalent to the stationarity of (II.1):

$`0=B\mathrm{\Psi }(t)`$ $``$ $`\mathrm{\Psi }(t+\sigma )\mathrm{\Psi }(t)`$ (63)
$`=`$ $`(I\sigma L)^1(I+\sigma W)\mathrm{\Psi }(t)\mathrm{\Psi }(t)`$
$`=`$ $`\left[(I\sigma L)^1(I+\sigma W)I\right]\mathrm{\Psi }(t)`$
$`=`$ $`(I\sigma L)^1\left[(I+\sigma W)(I\sigma L)\right]\mathrm{\Psi }(t)`$
$`=`$ $`(I\sigma L)^1[(\sigma (L+W)]\mathrm{\Psi }(t).`$
Solutions to (63) are found using Newton’s method. We begin with an initial estimate $`\mathrm{\Psi }`$, in our case the solution at a neighboring value of $`\mu `$. Newton’s method calls for approximating the nonlinear operator $`B`$ whose roots are sought by its linearization $`B_\mathrm{\Psi }`$ about $`\mathrm{\Psi }`$. We seek a decrement $`\psi `$ such that $`\mathrm{\Psi }\psi `$ solves this linearized equation:

$`0`$ $`=`$ $`B(\mathrm{\Psi }\psi )B(\mathrm{\Psi })B_\mathrm{\Psi }\psi `$
$`B_\mathrm{\Psi }\psi `$ $`=`$ $`B(\mathrm{\Psi }).`$ (64)
$`\mathrm{\Psi }`$ is then replaced by $`\mathrm{\Psi }\psi `$ and equation (64) solved again for a further decrement. The process is iterated until $`B(\mathrm{\Psi })`$ or $`\psi `$ is sufficiently small. In our case, equation (64) takes the form

$$(I\sigma L)^1\sigma (L+DW)\psi =(I\sigma L)^1\sigma (L+W)\mathrm{\Psi }.$$
(65)
We will explain how we solve the large linear problem (65) in section A.3.

The role of $`\sigma `$ is formally that of the timestep in (62), but in (65), its value can be taken to be arbitrarily large. For $`\sigma \mathrm{}`$, (65) becomes:

$$L^1(L+DW)\psi =L^1(L+W)\mathrm{\Psi }.$$
(66)
For the spherically symmetric case, the linear system (66) is of size $`N_R=128`$ and can be solved by standard Gaussian elimination. The boundary condition $`\psi (r=R)=0`$ is imposed by modifying the operator $`L^1`$ or $`(I\sigma L)^1`$, as it is in the time-stepping algorithm (62).

To compute the full branch of solutions as a function of $`\mu `$, we begin from a stable state of (61) at a small value of $`𝒩`$ obtained by time-integration. Each stationary state is computed in $`35`$ Newton iterations.

### A.3 Conjugate Gradient Solution of Linear Systems

For the periodic Cartesian case, the linear system (66) of size $`N_{3D}=10^6`$ is too large to be stored or inverted directly: the operation count for Gaussian elimination would be of the order of $`N_{3D}^3`$. Instead, we use BiCGSTAB Vandervorst , a variant of the well-known conjugate gradient method, developed for linear systems which are not symmetric definite. Such methods are matrix-free, meaning that they require only the right-hand-side of (66), and a subroutine which acts with the linear operator of the left-hand-side. A solution to the linear system is constructed as a carefully chosen linear combination of powers of the linear operator acting on the right-hand-side.

For a periodic Cartesian geometry, conjugate gradient methods are particularly economical, since operator actions are all accomplished in a time proportional to $`N_{3D}`$. However, conjugate gradient methods for nonsymmetric definite systems may converge slowly (requiring a large number of evaluations of the linear operator) or even not at all. This happens when the operator is poorly conditioned, i.e. roughly when it has a wide range of eigenvalues. One must then precondition the linear system, i.e. multiply both sides of the system by a matrix which improves its conditioning and accelerates convergence. Since for operators such as $`L+DW`$, the wide range of eigenvalues is due primarily to those of $`L`$, we expect $`L^1`$ to be an effective preconditioner. From (65), it can be seen that $`\sigma `$ allows us to interpolate between linear operators $`\sigma (L+DW)`$ and $`L^1(L+DW)`$. We vary $`\sigma `$ empirically to optimize the convergence of BiCGSTAB. A few hundred BiCGSTAB iterations are usually required to solve the linear system.

A further advantage of iterative inversion methods is that they can produce a (non-unique) solution even when the linear operator is singular. This is the case for our operators, which have the neutral modes described in section II.3, as well as other neutral modes related to symmetries and the Fourier representation. The preconditioner, however, is inverted exactly. If (66) is used, the constant Fourier mode is treated separately which allows us to construct an invertible version of the singular operator $`L`$.

### A.4 Eigenvalue Problem

We now describe our numerical method for calculating the linear stability of the stationary states. For the spherically symmetric case, the eigenvalues of (17) are computed by constructing and diagonalizing the corresponding matrix for each converged stationary solution. The results reported were generated with a Mathematica code running on a workstation. With the values $`R=4`$, $`N_R=128`$, the first two eigenvalues of the harmonic oscillator are obtained with a precision better than $`0.05\%`$.

For the three-dimensional case, it is again not possible to construct and diagonalize the matrix of size $`N_{3D}=10^6`$ directly: the operation count for diagonalization is also of the order $`N_{3D}^3`$. Instead, we calculate only eigenvalues of interest, using a variant of the iterative inverse power method. The inverse power method calculates the eigenvalues of a matrix $`M`$ closest to a value $`s`$ by means of the sequence defined by:

$$(MsI)\psi _{j+1}=\psi _j$$
(67)
The sequence $`\{\psi _j\}`$ converges rapidly to the eigenvector whose eigenvalue is nearest to $`s`$, with the eigenvalue $`\lambda `$ of $`M`$ estimated by $`(\lambda s)^1\psi _{j+1},\psi _j/\psi _j,\psi _j`$.

In order to calculate complex or multiple eigenvalues, and to obtain more precise eigenvalues and error estimates, we use the sequence generated by (67) to implement the more general Arnoldi or Krylov method Arnoldi . Instead of retaining only the last two members of the sequence, the last $`K`$ members (typically 4 or 6) are orthonormalized and then assembled into the $`K\times K`$ matrix $`H_{jk}\psi _j,(MsI)^1\psi _k`$. The eigenvalues of $`H`$ provide estimates of up to $`K`$ of the eigenvalues $`(\lambda s)^1`$ of $`(MsI)^1`$.

In our implementation of the Arnoldi method for (29a), we seek the eigenvalues $`\lambda ^2`$ of the matrix $`(L+DW^I)(L+DW^R)`$. Rather than solving (67), we solve the equivalent preconditioned problem TuckBark

$$L^2[(L+DW^I)(L+DW^R)sI]\psi _{j+1}=L^2\psi _j$$
(68)
by using BiCGSTAB. From (68) we obtain a sequence of vectors containing an increasing proportion of the desired eigenvectors, but since our solution of (68) is not exact, we then construct $`H`$ by multiplication rather than inversion via $`H_{jk}\psi _j,M\psi _k`$. We can then estimate the eigenvalues $`\lambda ^2`$ by those of $`H`$. Although the formal role of $`s`$ is that of a shift which focuses the inverse iteration on the eigenvalues being sought, here we also use it empirically to improve the convergence of BiCGSTAB.

The inverse Arnoldi method requires between 3 and 10 iterations to converge, each of which requires several hundred BiCGSTAB iterations in order to solve its associated linear system.


text:
# Memory Effects in Long-Time Lindblad Motion

## Abstract

It is shown that the Lindblad equation accounts for memory effects. That is to say, Lindblad operators can be constructed in a natural manner such that a memory term appears in the asymtotic $`\left(time\mathrm{}\right)`$ region; at the same time the expectation values depend on the initial state. Furthermore a procedure to extend the Lindblad equation to an equation of motion for an ideal Bose ’gas’ of ’particles’,i.e. systems with non-trivial internal structure,is described.Initially in some quantum state this collection of ’particles’ will asymptotically turn into an equilibrium ensemble whose probability distribution is determined by the Lindblad operators building the dissipative part of the equation of motion.

A)The Lindblad generalisation  of Schroedinger motion to dissipative motion hinges, apart from technical assumptions, only on very general physical notions:

* The abelian group of unitary Schroedinger motion generated by the Hamiltonian generalizes to a set of abelian semigroups characterized by a collection of operators $`V_J`$,the Lindblad operators. The semigroup structure - the non-existence of time-reversed motion - accounts for absorption.
* Complete positivity, interpreted physically, means that positive motions in Hilbert space of states (system 1) can be extended to a positive motion in the product space resulting by adjoining a second Hilbert space (system2), a construction leading to entanglement of both systems.

It should be noted that in the derivation of the Lindblad equation concepts used in the derivation of master equations - for instance the decomposition of the space of states in a product $`{}_{}{}^{}system^{}^{}bath^{}`$ \- do not play any role, at no place. Nor is there any conceptual relation to open systems.The relation of master-equations for open systems and Lindblad equations has been clarified to some extent      and shown to be controlled by relative scales.
We take up this observation and consider a system whose degrees of freedom interact with scale-dependent Hamiltonians and look for stationary states evolving from given initial states,i.e. we construct maps ($`\rho `$ is the density operator of our system)

$$\tau \left(V\right):\varrho |_{t=0}\varrho |_{t=\mathrm{}}$$
(1)
and discuss their dependence on the Lindblad operators $`V_J`$ which together with the Hamiltonian are supposed to differ for different scales - time scales,energy scales etc.
We present an explicit construction of asymptotic stationary states which will be seen to contain memory terms.

B)In this section we consider the case of only one Lindblad operator $`V`$ and write the Lindblad equation of motion

$$\dot{B}=i[H,B]+V^+BV\frac{1}{2}[V^+V,B]_+$$
(2)
where $`H`$ is the Hamiltonian and $`B`$ an observable.***V is assumed to be invertible; B,H are bounded operators acting in a separable Hilbert space for which the Lindblad equation has been proved.This fact allows us to use interchangeably the notions ’operator’ and ’matrix’ and treat the question of dimensions - finite or infinite - in a rather cavalier way.
Using the polar decomposition (U is a unitary operator)

$$V=U\sqrt{V^+V}$$
(3)
we rewrite (2) as ( note that the assumption of unitarity of $`U`$ excludes zero-modes in $`V^+V`$ )

$`{\displaystyle \frac{1}{\sqrt{V^+V}}}\dot{B}{\displaystyle \frac{1}{\sqrt{V^+V}}}=`$
$`i{\displaystyle \frac{1}{\sqrt{V^+V}}}[H,B]{\displaystyle \frac{1}{\sqrt{V^+V}}}+U^+BU`$
$`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{\sqrt{V^+V}}}B\sqrt{V^+V}+\sqrt{V^+V}B{\displaystyle \frac{1}{\sqrt{V^+V}}}\right)`$
The observation (see below) that

$$W:=\frac{1}{V^+V}$$
(5)
is a (non-normalized) probability distribution leads us to the physically plausible assumption

$$W=W(H,\mathrm{})$$
(6)
where the dots indicate further observables commuting with H. Tracing the equation of motion we immediately see that the trace of the rhs of the equation of motion vanishes identically and hence

$$tr(\dot{B}W)=0$$
(7)
or

$$tr(BW)=const$$
(8)
( $`\dot{W}=0`$ since $`W`$ depends only on conserved quantities). Needless to say we tacitly assume W to be traceclass.

In  we have demonstrated the following asymptotic form for B

* Irreducible V

$$B|_{t=\mathrm{}}=b(\mathrm{})𝕀$$
(9)
where $`𝕀`$ is the unit operator in $``$.
* Reducible V, i.e.

$$V=\underset{\alpha }{}V_\alpha $$
(10)
where the $`V_\alpha `$ are matrices in orthogonal subspaces $`^\alpha `$ of $``$, yield

$$B|_{t=\mathrm{}}=\underset{\alpha }{}b^\alpha (\mathrm{})𝕀^\alpha .$$
(11)

In the following we consider only the irreducible case and derive

$$B|_{t=\mathrm{}}=\frac{tr(B|_{t=0}W)}{tr(W)}𝕀$$
(12)
The expectation value of the asymptotic configuration then is

$$=tr(B|_{t=0}\varrho _0)=\frac{tr(B|_{t=0}W)}{tr(W)}$$
(13)
for all states $`\varrho _0`$, i.e. the expectation value is independent of the initial state, no memory effects are present.We see that

$$P_W=\frac{W}{tr(W)}$$
(14)
is a normalized probability distribution derived from the Lindblad operator $`V`$.Translating this result into the Schroedinger picture we derive that any initial state $`\varrho _0`$ tends to $`P_W`$ for $`t\mathrm{}`$, i.e.

$$\tau \left(V\right):\varrho _0\varrho _0|_{t=\mathrm{}}=P_W$$
(15)
for all initial states $`\varrho _0`$.

We now turn to the question of memory effects. To show that they can be incorporated we extend the Lindblad equation, without changing its formal content, to an equation of motion for quantum subsystems separated, e.g. by scales, from the system built up by these subsystems. As an example we could take a molecule: the subsystems are spanned by the states corresponding to the inner degrees of freedom of the atoms composing the molecule, the system - the molecule - is built up by the atomic states of outer shells.
To realize this construct we endow the input matrices V and H with a direct product structure and, to simplify matters, choose the ansaetze

$$V=\left(\stackrel{~}{V}_{ik}𝕀\right)\sqrt{n}$$
(16)
$$H=\left(\stackrel{~}{H}_{i,k}\right)$$
(17)
$$\stackrel{~}{V}_{ik},\stackrel{~}{H}_{ik}$$
(18)
where $`𝕀`$ is the $`n\times n`$ unit matrix and $``$ is a $`n\times n`$ matrix sub-Hamiltonian, identical for all sites $`(i,k)`$, i.e.$`\stackrel{~}{V}`$ and $`\stackrel{~}{H}`$ are matrices with $`n\times n`$ matrix valued entries indexed by $`(i,k)`$. Our ansatz for V guarantees that,in the terms of our example, the Lindblad operator leaves the inner degrees of freedom unaffected. The observable B is written as a matrix of $`n\times n`$ matrices $`B_{ik}`$

$$B=\left(B_{ik}\right)$$
(19)
The probability distribution is then

$$W=(\left(\stackrel{~}{V}^+\stackrel{~}{V}\right)^1𝕀)=:(\stackrel{~}{W}_{ik}𝕀),$$
(20)
V has the polar decomposition

$$V=U\sqrt{V^+V}=\left(\stackrel{~}{U}𝕀\right)\sqrt{\stackrel{~}{V}^+\stackrel{~}{V}𝕀}.$$
(21)
Tracing the equation of motion with respect to the indices $`(i,k)`$ then yields instead of (7)

$$\stackrel{~}{Tr}\left(\dot{B}W\right)=\stackrel{~}{Tr}\left(\stackrel{~}{W}\stackrel{~}{H}B\right)\stackrel{~}{Tr}\left(\stackrel{~}{W}B\stackrel{~}{H}\right)$$
(22)
which leads to

$$\stackrel{~}{Tr}\left(BW\right)=C+_0^t\left(\stackrel{~}{Tr}\left(\stackrel{~}{W}\stackrel{~}{H}B(t)\right)\stackrel{~}{Tr}\left(\stackrel{~}{W}B(t)\stackrel{~}{H}\right)\right)𝑑t$$
(23)
It has to be stated that taking the trace of this $`n\times n`$ matrix we obtain a vanishing result

$$Tr\left(\dot{B}W\right)=tr_n\stackrel{~}{Tr}\left(\dot{B}W\right)=0$$
(24)
in accordance with equation (7) .This is because we have assumed

$$0=[W,H]=[\stackrel{~}{W},\stackrel{~}{H}]$$
(25)
and, thus

$$[\stackrel{~}{W},\stackrel{~}{H}]=0.$$
(26)
Following the derivation given in  we find for the asymtotic configuration

$$B|_{t=\mathrm{}}=b(\mathrm{})\left(\delta _{ik}\right)$$
(27)
where $`b(\mathrm{})`$ is now a $`n\times n`$ matrix which is read off

$`b\left(\mathrm{}\right)=`$
$`{\displaystyle \frac{1}{\left(\stackrel{~}{Tr}\stackrel{~}{W}\right)}}\times `$
$`\left(\stackrel{~}{Tr}\left(B|_{t=0}W\right)+{\displaystyle _0^{\mathrm{}}}\left(\stackrel{~}{Tr}\left(\stackrel{~}{W}\stackrel{~}{H}B(t)\right)\stackrel{~}{Tr}\left(\stackrel{~}{W}B(t)\stackrel{~}{H}\right)\right)𝑑t\right).`$
We note the explicit appearance of a memory term.Calculating the expectation value of $`B|_{t=\mathrm{}}`$ in some state $`\varrho _0`$ written as

$$\varrho _{(0)}=\left(\varrho _{ik}^{(0)}\right)$$
(29)
where the $`\varrho _{ik}^{(0)}`$ are $`n\times n`$ matrices, we find

$$=\underset{i}{}tr_n\left(\varrho _{ii}^{(0)}b(\mathrm{})\right)$$
(30)
and observe that now the asymptotic expectation value does depend on the initial state in concordance with the appearance of a memory term.

So we have seen that a simple and intuitively clear generalisation of the Lindblad equation to an equation for dynamical degrees of freedom of subsystems leads to memory effects; the asymptotic subsystem variables given in equations (23) should be interpreted as the new dynamical subsystem variables obtained from an asymptotic averaging procedure over those degrees of freedom of the total system living on ’lower’ scales:the equations (23) are clearly seen as an elimination procedure for ’environment’ variables separated into a statistical average and a memory term.

C)We now turn to the case of more than one, say N, Lindblad operators. We take $`N`$ finite with the provision of eventually letting $`N\mathrm{}`$ as certain physical models might require.The equation of motion then reads

$$\dot{B}=i[H,B]+\underset{J}{}\left(V_J^+BV_J\frac{1}{2}[V_J^+V_J,B]_+\right)$$
(31)
We rewrite this equation as an equation operating in a direct sum of identical spaces

$$_N=\underset{1N}{}$$
(32)
and define

$$B_N:=B𝕀_N$$
(33)
$$H_N:=H𝕀_N$$
(34)
$$V_N:=\sqrt{N}\left(V_J\delta _{JK}\right)$$
(35)
to arrive at

$$\dot{B}_N=i[H_N,B_N]+V_N^+B_NV_N\frac{1}{2}[V_N^+V_N,B_N]_+.$$
(36)
The polar decompositions

$$V_J=U_J\sqrt{V_J^+V_J}$$
(37)
lead to the polar decomposition

$$V_N=U_N\sqrt{V_N^+V_N}$$
(38)
where

$$U_N=\left(U_J\delta _{JK}\right)$$
(39)
is unitary. Employing the same procedure as above we have

$`{\displaystyle \frac{1}{\sqrt{V_N^+V_N}}}\dot{B}_N{\displaystyle \frac{1}{\sqrt{V_N^+V_N}}}=`$
$`i{\displaystyle \frac{1}{\sqrt{V_N^+V_N}}}[H_N,B_N]{\displaystyle \frac{1}{\sqrt{V_N^+V_N}}}+U_N^+BU_N`$
$`{\displaystyle \frac{1}{2}}\left(\sqrt{V_N^+V_N}B_N{\displaystyle \frac{1}{\sqrt{V_N^+V_N}}}+{\displaystyle \frac{1}{\sqrt{V_N^+V_N}}}B_N\sqrt{V_N^+V_N}\right)`$
. (40)
Assuming either independence of $`W_J`$ on J ( $`U_J`$ does depend on $`J`$ in general) or, alternatively, $`V_J`$ positive and in analogy to (6)

$$W_J:=\frac{1}{V_J^+V_J}=W_J(H,\mathrm{})$$
(41)
and taking the total trace ($`tr_N`$ pertains to the matrix indices of the $`N\times N`$ matrices introduced above, $`tr_{}`$ to the operators on $``$) we find

$$tr_{}tr_N\left(\dot{B}_NW_N\right)=0$$
(42)
and thus

$$tr_{}\left(\dot{B}\underset{1N}{}W_J\right)=0$$
(43)
and

$$=\frac{tr_{}\left(B|_{t=0}W\right)}{tr_{}\left(W\right)}$$
(44)
$$W=\underset{J}{}W_J$$
(45)
To illustrate this result we take $`V_J`$ positive and assume the following specific ansatz

$$W_J=\left(V^+V\right)^J$$
(46)
$`J=1,2,\mathrm{}`$ and put

$$V=\mathrm{exp}\left(\frac{\beta H}{2}\right)$$
(47)
The expectation value for B reaches asymptotically

$$=\frac{1}{\overline{N}}tr_{}\left(\frac{B|_{t=0}}{\mathrm{exp}\left(\beta H\right)1}\right)$$
(48)
which is simply the expectation value of the ’particle’ observable B in an ideal Bose ’gas’ of, in the average, $`\overline{N}`$ ’particles’ at inverse temperature $`\beta `$; ’particle’ is just a more intuitive name for the physical object dubbed ’system’ up to now.
This interpretation deserves further clarifications. To this end we reformulate equation (Memory Effects in Long-Time Lindblad Motion) as an equation in Fock space $`_F`$, aiming at the Bose nature of the ideal gas to be introduced.We define

$$_F:=\underset{J}{}^J$$
(49)
and

$$V_F:=\underset{J}{}V^J$$
(50)
so that

$$W_F=V_F^+V_F=\underset{J}{}\left(V^+V\right)^J.$$
(51)
In the product space $`^J`$ we select as physically relevant states symmetric states which we take as superpositions of symmetric ’system’ product states - we introduce many-’particle’ boson states .The observable $`B`$ is extended to a symmetrically operating operator

$$B_F=\underset{J}{}\left(\mathrm{}𝕀B\mathrm{}\right)$$
(52)
where B stands consecutively on all positions of the J-fold product.
The expectation value at $`t=\mathrm{}`$ is, in strict analogy to (13) and (44)

$$<B_F|_{t=\mathrm{}}>=\frac{tr__F\left(B_F|_{t=0}W_F\right)}{tr__F\left(W_F\right)}$$
(53)
where the trace is now to be calculated with a symmetric product basis in $`^J`$ for all J. Computing this trace one encounters disconnected terms (matrix elements now pertain to $``$)

$$\underset{i_1,\mathrm{},i_L}{}<i_1|BV^{J_1}|i_1><i_2|V^{J_2}|i_2>\mathrm{}<i_L|V^{J_k}|i_L>$$
(54)
$$\underset{l}{}J_l=J$$
(55)
All these terms sum up to the same common factor in the nominator and denominator - the connected cluster theorem - so that we reproduce (48) with W$`{}_{J}{}^{}=(V^+V)^J`$. We conclude that any symmetric quantum many-particle state composed of whatever complex quantum systems - a Bose many-particle state - is transported by Lindblad motion into an equilibrium ensemble with a probability distribution

$$W_{Bose}=\frac{1}{W_{}1}$$
(56)
where

$$W_{}=\frac{1}{V^+V}$$
(57)
is an operator acting in the space of the system’s states. This derivation is a first step towards a dissipative quantum field theory: the case of free fields, although we never explicitly introduced this concept.I shall return to the extension to more complicated cases in a forthcoming publication.

Acknowledgment:
It is a pleasure to thank Jan-Michael Rost, Max Planck Institut fuer die Physik Komplexer Systeme, Dresden and Axel Schenzle,Sektion Physik, Universitaet Muenchen for the stimulating ambiente during my stay at these instutions where this work has been done.