# Nucleation of Spontaneous vortices in trapped Fermi gases undergoing a BCS-BEC crossover

###### Abstract

We study the spontaneous formation of vortices during the superfluid condensation in a trapped fermionic gas subjected to a rapid thermal quench via evaporative cooling. Our work is based on the numerical solution of the time dependent crossover Ginzburg-Landau equation coupled to the heat diffusion equation. We quantify the evolution of condensate density and vortex length as a function of a crossover phase parameter from BCS to BEC. The more interesting phenomena occur somewhat nearer to the BEC regime and should be experimentally observable; during the propagation of the cold front, the increase in condensate density leads to the formation of supercurrents towards the center of the condensate as well as possible condensate volume oscillations.

## Introduction

Atomic gases undergoing a superfluid condensation are very clean and simple systems to address non-equilibrium quantum dynamics of superconductors. As an example of important non-equilibrium behavior, the present paper focuses on the way in which a rapid quench from the normal phase leads to the formation of a random tangle of vortices. In addition to cosmology Zurek (), the character and dynamics of these defects is of great interest to a number of different physics subdisciplines ranging from atomic phy sics, weiler+n08 (); berloff+pra02 () to condensed matter bunkov (); ruutu () and high temperature superconductivity Ong2 (); carmi (); maniv (). To this end, theoretical investigations of rapid quenches of atomic Bose gases have been undertaken weiler+n08 () and appear consistent with experiments. In this Rapid Communication we make predictions for future observations on trapped Fermi gases throughout the entire crossover from BCS to BEC, with special emphasis on the unitary regime.

These fermionic systems are in many ways more suitable systems to study the dynamics of quantum gases. The Pauli principle suppresses three body collisions making it possible to study the strong interaction regime. For the Bose counterparts, collisions are frequent and lead to condensate decoherence. A related advantage of the Fermi gases is the opportunity to explore tunable interaction strengths which, in turn, affects the dynamics. This tunability is associated with the crossover between the BCS (where the inter-fermionic attraction is weak) and the Bose-Einstein condensed (BEC) limits (where the attraction is strong). This crossover is entirely accessible through application of a magnetic field, and the exploitation of so-called Feshbach resonances Levinreview (); Zwierlein (). Our work is based on a numerical simulation of of the time-dependent Ginzburg-Landau equation (TDGLE) aranson+prl99 (); aranson+rmp02 () in the presence of thermal (white) noise . This TDGLE equation, which in many ways is one of the most fundamental equations in condensed matter physics, represents a differential equation characterizing the dynamics and spatial dependence of the pairing gap parameter. Microscopically, one can demonstrate that Randeria (); Maly (); Levinreview () the TDGLE coincides with the energy conserving Gross-Pitaevski description in the strong attraction or BEC limit, where the dissipation is minimal Similarly in the weak attraction limit the TDGLE leads to the well known diffusive dynamics of BCS theory.

It should be stressed that TDGL theory refers to the gap or order parameter, and not to the underlying fermions. In this way the number of fermions, or the Fermi-Fermi interaction parameters do not enter, as they have been effectively integrated out. In contrast to static GL theory, which is a consequence of Gorkov theory, there is no fully rigorous derivation of this dynamics; alternative dynamical schemes such as time dependent Bogoliubov deGennes (TDBDG) theory BulgacScience (), which have similar antecedents in Gorkov theory, are also not rigorous.

By contrast to this TDBDG approach, our studies are in the highly non-equilibrium regime. We stress that there are essentially no alternative effective simulation techniques for addressing the condensate evolution, associated with rapid quenches. Monte Carlo and first principle quantum dynamics are limited to small system sizes, short times and typically zero temperature; moreover, they deal with the dynamics of single particles, not the collective dynamics of the order parameter (in the strong coupling limit) which we consider here.

## Theory

We model the temperature equilibration in a physical way as a non-uniform process, associated with evaporative cooling. The effects of the BCS-BEC crossover enter most prominently via the change in the complex time relaxation rate of the TDGLE Randeria (); Maly (); Levinreview (). As a function of the continuous crossover from BCS to BEC we see that the number and lifetime of spontaneous vortices increases, so that the steady state condensate is slower to form. Among the most interesting observations in this paper pertains to the near unitary regime (or mid-point between BCS and BEC) where the interactions are strongest and the TDGLE- based approach is the most appropriate.

For notational purposes, it is convenient to introduce a parameter that tunes the equation from BCS () to BEC (). Physically this parameter can be viewed as representing an applied magnetic field in conjunction with a Feshbach resonance. This tunes the strength of interaction between the fundamental fermionic particles (see e.g. Randeria (); Maly (); Levinreview () for a rigorous derivation).

The TDGLE is given by

(1) |

where is the order parameter and is uniformly distributed thermal noise with fluctuation temperature .

The coefficient takes into account the evaporative cooling of the system and the trapping potential, i.e. , where is the trapping potential which is determined by the trap size. Eq. (1) is written in dimensionless units, with length scale , time scale , and the order parameter scale . Temperatures are measured in units of . For (fermionic) Li these values are estimated as: m, s, and concentration m. The dimensionless fluctuation temperature for a condensate at nK is in the interval , depending on the condensed atoms.

The trap is assumed to be spherically-symmetric. Thus we calculate the temperature profile due to evaporative cooling by solving the spherically-symmetric heat diffusion equation

(2) |

with an initial uniform temperature and final temperature that is fixed at the boundary . The heat diffusion constant is renormalized by ; typical values for Li: .

All simulation runs start with random initial conditions in the normal phase () at temperatures larger than the critical temperature . The evaporative cooling is initiated at time on the surface of our spherical trap with radius . The cold front surface propagates towards the center and quenches the atomic gas below the critical temperature. The dimensionless diffusion constant was typically set to and the fluctuation strength . This mechanism of condensate nucleation is to be contrasted with that described in Ref. weiler+n08 () where the cooling occurred uniformly in space, thus corresponding to an infinite diffusion constant. Another crucial difference with this previous work is the equilibration mechanism: the time evolution of the condensate used in weiler+n08 () was based on the Gross-Pitaevskii model where the dissipation was introduced by truncation of high-frequency modes. By contrast, our TDGLE model includes a complex relaxation rate which is determined by the crossover physics through the phase factor . Experiments on these cold gases Levinreview () cannot probe very deeply into the BCS or BEC regimes, but are confined to the so-called “unitary” midpoint. Nevertheless, at unitarity, pairs are reasonably long lived Maly () so that we can associate with the physically accessible regime.

## Numerical Results

Our numerical calculations were done for volumes discretized in up to grid points, averaged over up to different initial conditions and time evolved for different crossover phase values. The timestep in dimensionless units was chosen to be and the total simulation time up to at -values close to the BEC limit. Our quasi-spectral split step method to solve the TDGLE which uses fast Fourier transforms, is much more stable than the traditional finite difference method. As a result of employing modern graphics processing units, our systems are an order of magnitude larger than in recent work on the bosonic BEC weiler+n08 (). This allows us to simulate more realistic physical situations CUDA ().

We begin with an illustration of the time evolved condensation process for the BCS limit () and a near-BEC situation (). Throughout this paper we avoid the strict fermion-based BEC limit () since without dissipation the condensate does not form and nucleation of vortices is completely suppressed. By introducing a complex relaxation rate in the TDGLE, we avoid having to include the interactions between condensed and non-condensed pairs which do not enter as naturally, in the fermionic- BEC limit. These were essential for addressing the bosonic counterpart experiments weiler+n08 (); bradley+pra08 ().

Plotted in Fig. 1 are the isosurfaces for constant condensate magnitude. We also show cross-sections in the -plane through the center of the system, indicating the phase of the order parameter . These latter plots appear below the counterpart isosurface pictures and contain information about spontaneous vortices. The isosurfaces for the condensate magnitude are color-coded according to fixed condensate density values , where we used . The three different panels from left to right correspond to three different times close to the typical time at which the condensate forms, i.e. the time when the condensate volume reaches half of its saturation value.

One can see from Fig. 1, that most spontaneous vortices are connected to the surface of the condensate. Their time evolution is entirely accessible in our simulations suppl () and one can follow their decay in the bulk and at the condensate surface. In the phase cross-section plots, these vortices are even more visible appearing as topological defects (phase singularities). Although vortices are most plentiful at the surface of the condensate, we find their decay occurs mostly in the bulk.

It is clear that in the BCS limit the condensate forms quickly occupying the entire trap volume in a relatively short time period. At much longer time scales (than shown here), the condensate will appear uniform. The near-BEC limit exhibits an interesting contrast: Here the condensate expands slowly, and vortices are much more plentiful, decaying even more gradually. An interesting feature is shown in the middle panel of the last row of Fig. 1. We observe in the course of condensate formation, concentric rings in the phase cross-section interrupted by trapped topological defects. This corresponds to the generation of supercurrents from the surface of the condensate towards its center. The phenomenon of transient supercurrent generation can be understood as follows: In the BEC limit the TDGLE can be written in the hydrodynamic form landau_hf () with the condensate density satisfying the continuity equation . During the quench, changes from zero to its maximal value leading to the formation of supercurrents.

Figure 2 presents a plot of the time evolution of the condensate volume , as well as the vortex length for different values of the crossover phase . Both quantities are normalized to the asymptotic final condensate volume. Here the condensate volume is calculated from the number of grid points with condensate density larger than . Assuming a spherical shape of the condensate (reflecting the trapping potential) we establish the diameter of the condensate; we then calculate the number of grid points with inside this sphere to estimate the vortex length. In the BCS limit the dissipation is maximum which is responsible for the rapid formation of the condensate. Increasing as one approaches the near-BEC gradually delays the formation of the condensate leading to increase of the maximum total vortex length.

The vortices show a roughly -decay of the total vortex length berloff+pra02 (); Zurek (). In Fig. 2b the time dependence of the total vortex length is plotted from the BCS to the near-BEC limit. The inset shows some of these same curves on a log-log scale and compares with a linear curve (dashed) representing a decay. In two dimensions, bulk vortices annihilate in pairs, the relaxation should be proportional to the square of the vortex concentration, i.e. the relaxation follows a power-law vinen (). Similarly, a is expected for the decay of vortex lines in 3D for the BCS case and near the BEC limit Kobayashi (). The deviations from power law are likely caused by finite size effects.

The condensate volume shows a well pronounced overshoot with subsequent oscillations in time. This appears around unitarity, albeit closer to the BEC limit. This overshoot is not generic, since it disappears in the BCS limit, presumably as a consequence of the large dissipation. In the near-BEC limit, where the dissipation is strongly suppressed, the condensate takes longer to form than the time associated with a thermal quench. We quantify this “condensate breathing” by studying the height of this condensate volume overshoot in the inset of Fig. 3. Plotted here is the maximum of the volume of the condensate as a function of . Also indicated is the maximum vortex length within the condensate. It can be seen that the maximum of the condensate overshoot occurs at . The appearance of the overshoot is also characterized by the time scale , where the condensate reaches its maximum, which diverges below (since there is no overshoot) and near the BEC limit (because the time of the occurrence of the overshoot diverges for ). A minimum of is observed at , i.e., the overshoot happens fastest. Fig. 3 indicates the typical time scale of the condensate formation and the nearly identical time scale , which is the time at which the vortex length is maximum. Both these times generally depend only weakly on , except near . Also the typical decay time of the vortices appears to be nearly independent of below (see Fig. 2b). Using the experimental time scales for Li, we can conclude from Figs. 2 & 3, that the condensate forms within a few seconds and vortices annihilate on the same time scale.

The overshoot and subsequent time oscillations of the condensate may be related to excitation of sound waves produced in a spherical trap by the quench. To further clarify this phenomenon we solved the spherically-symmetric TDGLE and obtained similar behavior as shown in the inset to Fig. 2a; the main difference being a delay in the condensate nucleation oneD (). In the process we verified that the period of these oscillations is close to the time needed for sound waves to traverse the condensate. Indeed in our units the speed of sound (which emerges from the BEC limit of the TDGLE) is equal to unity in the long-wavelength limit, while the oscillation period obtained from Fig. 2a is approximately time units. This is close to the time needed for the sound wave to cross our system with diameter length units. Using the value for and for Li, the speed of the second sound is on the order m/s, comparable to the mean velocity of the atoms at nK.

There has been recent excitement over numerical approaches for addresses vortices in Fermionic superfluids, particularly with the implementation BulgacScience () of a time dependent Bogoliubov deGennes numerical scheme. Unlike TDGLE, this approach addresses the fermionic degrees of freedom in conjunction with the fermionic gap parameter. Generally it has been applied to . Importantly, it does not incorporate the fact that the superconducting order parameter need not be the same as fermionic excitation gap. Related to this last observation is the fact that these schemes do not accomodate non-condensed bosons which are expected to enter into the dynamics away from the strict BCS regime. As a consequence, our understanding of vortices in fermionic superfluids will require the exploration of a number of different numerical approaches, including the TDGLE scheme we address here.

We thank Kara Lamb, Matt Davis, Chih-Chun Chien and Nate Gemelke for useful discussions. This work was supported by the by the U.S. DOE, Office of Basic Energy Sciences, Division of Materials Science and Engineering, under Contract DEAC02-06CH11357, and by NSF-MRSEC Grant No. 0820054 (KL).

## References

- (1) N.D. Antunes, L.M.A. Bettencourt, W.H. Zurek, Phys. Rev. Lett. 82, 2824 (1999).
- (2) C. N. Weiler T W. Neely, D. R. Scherer, A. S. Bradley, M. J. Davis, and B. P. Anderson, Nature 455, 948 (2008).
- (3) N.G. Berloff and B.V. Svistunov, Phys. Rev. A 66, 013603 (2002).
- (4) C. Bäuerle, Yu. M. Bunkov, S. N. Fisher, H. Godfrin, and G. R. Pickett, Nature 382, 332 (1996).
- (5) V. M. H. Ruutu, V. B. Eltsov, A. J. Gill, T. W. B. Kibble, M. Krusius, Yu. G. Makhlin, B. Placais, G. E. Volovik, and W. Xu, Nature 382, 334 (1996).
- (6) Y. Wang, Z. Xu, T. Kakeshita, S. Uchida and N. P. Ong, Phys. Rev. B. 64, 224519 (2001).
- (7) R. Carmi, E. Polturak, G. Koren, and A. Auerbach, Nature 404, 853 (2000).
- (8) A.Maniv, E. Polturak and G. Koren, Phys. Rev. Lett. 91, 197001 (2003).
- (9) Q.J. Chen, J. Stajic, S, Tan and K Levin, Phys. Rep. 412, 1 (2005).
- (10) M. Zwierlein, C. Schunck, A. Schirotzek, W. Ketterle, Nature 442, 52 (2006).
- (11) I. S. Aranson, N. B. Kopnin, and V. M. Vinokur, Phys. Rev. Lett. 83, 2600 (1999); Phys. Rev. B 63, 184501 (2001)
- (12) I. Aranson and L. Kramer, Rev. Mod. Phys. 74, 99 (2002).
- (13) C. A. R. Sa de Melo, M. Randeria, and J. R. Engelbrecht, Phys. Rev. Lett. 71, 3202 (1993)
- (14) J. Maly, B. Jankó, and K. Levin, Physica C 321, 113 (1999).
- (15) We used a cluster of Nvidia Tesla C2050/2070 cards) that increased the speed of our simulations by up to a factor 100 compared to traditional CPU-based computations.
- (16) A. Bulgac, Y-L Luo, P. Magierski, K. J. Roche, and Y. Yu, Science 332, 1288 (2011)
- (17) A.S. Bradley, C.W. Gardiner, and M.J. Davis, Phys. Rev. A 77, 033616 (2008).
- (18) Supplementary information and movies can be found at PROLA or http://mti.msd.anl.gov/highlights/bcsbec .
- (19) T. Frisch, Y. Pomeau and S. Rica, Phys. Rev. Lett. 69, 1644 (1992).
- (20) W. F. Vinen, Proc. R. Soc. London A 242, 493 (1957).
- (21) M. Kobayashi and M. Tsubota, J. Low Temp. Phys, 145, 209 (2006).
- (22) The spherically-symmetric TDGLE cannot capture the full dynamics and does not admit vortex solutions because of the reduced dimensionality. Correspondingly the noise amplitude is reduced and becomes a function of the distance from the center of the condensate.