Relaxation dynamics of half-quantum vortices in a two-dimensional two-component Bose-Einstein condensate

We study the relaxation dynamics of quantum turbulence in a two-component Bose-Einstein condensate containing half-quantum vortices. We find a temporal scaling regime for the number of vortices and the correlation lengths that at early times is strongly dependent on the relative strength of the inter-species interaction. At later times we find that the scaling becomes universal, independent of the inter-species interaction, and approaches that numerically observed in a scalar Bose-Einstein condensate.

Introduction. -Since the realization of superfluidity, quantum turbulence (QT) has been studied in systems ranging from superfluid liquid Helium [1,2] to quasi-particle condensates in solid-state systems [3]. Due to their unprecedented experimental accessibility, QT in Bose-Einstein condensates (BECs) in dilute, ultracold atomic gases have attracted considerable theoretical [4][5][6][7][8][9] and experimental [10][11][12][13][14][15] interest in both 2D and 3D configurations. In a scalar BEC, the QT state is made up of a large number of vortices with quantised circulation. The collective behaviour of the vortices plays a key role in the hydrodynamics, recovering features of classical turbulence that can exhibit the characteristic Kolmogorov power-law spectrum [16].
In contrast to the scalar superfluids, multicomponent and spinor BECs are described by multicomponent order parameters and allow for a wider range of topological defects, which give rise to novel dynamics [17][18][19][20]. Consequently, there has been increasing interest in the properties of QT and non-equilibrium dynamics in such systems [21][22][23][24][25]. The simplest non-scalar topological excitation appears in a two-component BEC, described by two complex fields, as the appearance of a phase singularity in only one component. When the atomic mass and mean density of the components are equal, such vortices are often referred to as half-quantum vortices (HQVs), due to their similarities with vortices carrying half a quantum of superfluid circulation in superfluid 3 He [26,27] and spin-1 BECs [28,29]. The study of QT in BECs can be separated into two distinct categories: 1) forced turbulence where a statistically stationary state is established; 2) decaying turbulence where a non-equilibrium initial condition, typically involving vortices, relaxes towards equilibrium. Here, we numerically investigate the spatial and temporal properties of the relaxation dynamics of a non-equilibrium initial state in a two-dimensional two-component system containing HQVs. Using a pseudospin interpretation, we compute the temporal scaling of the correlation functions associated with the spin-and mass-superfluid ordering. We relate these to the vortex decay rate and analyse how this depends on the intra-component interaction strength.
We contrast our observations for this system with similar simulations that have been performed for scalar BECs and reported in [30][31][32].
The two-component BEC. -We consider an untrapped two-component BEC described by the Gross-Pitaevskii (GP) mean-field theory subject to periodic boundary conditions. The dynamics of the condensate is described by the two coupled GP equations where ψ j is the condensate wavefunction and m j (j = 1, 2) is the atomic mass for the jth component. The strength of inter-and intra-component interactions are described by g j and g 12 , respectively. We consider a condensate where m 1 = m 2 = m, as is the case, e.g., when the two components are different hyperfine states of the same atomic species, and also assume g 1 = g 2 = g. The key parameter is then the ratio of intra-to inter-species interactions which in experiment could be tuned using magnetic [33] or microwave-induced [34] Feshbach resonances. Here we consider 0 < γ < 1, such that all interactions are repulsive, while keeping the condensate stable against separation of the components. The vortex states of the two-component BEC may be understood as follows: We write the two-component wavefunction as the vector (ψ 1 , ψ 2 ) T . Taking θ j = Arg(ψ j ), this may be decomposed as where Gradients in Φ can then be interpreted in terms of pseudospin currents, while gradients in Θ may be associated with a total, superfluid mass current. Now consider a vortex state consisting of a phase singularity in ψ 1 , around which θ 1 winds by 2π while θ 2 remains unchanged, such that where φ is the azimuthal angle around the vortex. The vortex is thus equivalently described by a π change in Θ (and a simultaneous π change in Φ) along a closed path encircling the vortex. Since Θ can be associated with a total mass current in the two components together, these vortex states are often referred to as HQVs and we adopt this language from here on. However, the two-component vortices are topologically distinct from HQVs in the A and polar phases of superfluid 3 He [26,27] and in the uniaxial nematic phase of spin-1 BECs [28,29]. A pseudospin picture also allows us to understand the size of HQV cores in terms of an energetic hierarchy of length scales arising from the inter-and intra-component interactions. These length scales are associated, respectively, with variations of the total superfluid density and of the density difference between the components. We thus define the density and spin healing lengths as [35] where n 0 is the number density of each component in a uniform system. Since a HQV consists of a phase singularity in only one condensate component, the remaining component is free to fill the vortex core. This can be interpreted as a variation of the pseuodspin z-component, whose size is determined by the spin healing length. When ξ s ξ d , the vortex core can thus expand, lowering the total energy. Therefore, γ directly determines the sizes of the vortex cores in the system. A similar energetic hierarchy of length scales leads to dramatic defect-core deformations in spinor BECs [36], including splitting of singly quantised vortices into HQVs [29,37].
Numerical method. -To study the dynamics of vortices in a turbulent regime we numerically evolve the time-dependent two-component Gross-Pitaevskii equations using a split-step algorithm [38]. We write eq. (1) in terms of the dimensionless variables:r = r/a s ,t = t/τ , g = 2mg/ 2 andψ j = a s ψ j , where a s is the lattice spacing and τ = 2ma 2 s / is the lattice time. The resulting equations then become where we have dropped the tildes for notational convenience. Our simulations were performed on a periodic domain of non-dimensional area L 2 with side length L = N s where N 2 s is the number of grid points. We solve eq. (7) on a grid of 1024 2 points with a s = 1. Motivated by similar work in a scalar BEC [32], we take N = 3.2 × 10 9 atoms per component and dimensionless g = L 2 /4N . The non-dimensional density healing length is thus fixed at N s /(gN ) 1/2 = 2. We now explore the role of the inter-component interaction by varying γ within the range 0 < γ < 1.
The initial condition for the GP evolution is constructed as a grid of vortex positions containing 48 2 vortices in each component with the grids of each component offset in both x and y to avoid overlapping positions. We then add a small, random displacement to each position to create an irregular distribution of vortices. This facilitates the development of an initially chaotic and subsequently turbulent vortex evolution during the relaxation dynamics. The phase of each component is subsequently constructed as an alternating 2π winding around each vortex position using the method described in ref. [5] that also accounts for the periodic boundary conditions. An initial short period of imaginary-time propagation, keeping the phase profile fixed, allows the vortex cores to form. The resulting HQVs consist of a density depletion in one component at the position of the phase singularity, which is then filled with atoms of the other component, as illustrated in fig. 1. From this initial state, the system is evolved according to eq. (7). HQVs with opposite circulation but with the phase singularity in the same component may annihilate which leads to a decay of the total vortex number. Results. -We first investigate the effect of γ on the relaxation dynamics of HQVs. Fig. 2(a)-(c) shows the density of the ψ 1 component for γ = 0.1, 0.6, 0.8. HQVs with a phase singularity in this component are readily apparent by the corresponding density depletion, and have a core size that grows with increasing γ. For γ ≥ 0.6, high density peaks also become noticeable and correspond to the positions of HQVs with phase singularity in ψ 2 . This can be understood from the healing lengths, eq. (6). For small γ, ξ s ∼ ξ d . As γ increases, the spin healing length also increases. Consequently, the cores of the HQVs fill with atoms from the other component as the resulting lowering of the kinetic energy offsets the cost in interaction energy. This causes the vortex cores to expand to a size similar to the spin healing length, as borne out by our simulations.
To track the vortex positions, we evaluate the pseudovorticity [39,40] where is the mass current of component j = 1, 2. The pseudovorticity remains regular and non-zero within the core of each vortex, and relaxes to zero away from the vortex singularity (at length scales exceeding the spin-healing length, ξ s ), as shown in fig. 2(d)-(f). The sign of the pseudo-vorticity also determines the charge of the vortex. The pseudo-vorticity shows the vortex positions particularly sharply for small γ, where the vortex cores are small.
We now investigate the spatial properties of our turbulent system. We split the kinetic energy, E kin = E v + E q , into classical (E v ), and quantum-pressure (E q ) contributions. These are given by where n j = |ψ j | 2 for j = 1, 2. The energy spectra for these contributions involve the Fourier transforms of the generalised velocities for the incompressible (i), compressible (c), and quantum pressure (q) parts [21], defined as The incompressible and compressible components of the velocity field are recovered from a Helmholtz decomposition into a divergence free, incompressible part ∇ · v i = 0, and an irrotational, compressible part ∇ × v c = 0. The kinetic energy spectrum can then be calculated by integrating the corresponding Fourier transforms over the full k-space angle for wave number k = |k|. The total kinetic energy is given by integrating over all k and summing over the different contributions: E kin = δ dkE δ (k) for δ = (i, c, q). The occupation numbers corresponding to the different energy contributions are then Fig . 3 shows the occupation number for each energy contribution along with the total occupation number n(k) for the case of γ = 0.6 at a time t = 2 × 10 5 τ . The total single-particle spectrum obeys the predicted scaling n(k) ∼ k −4 in the infrared (IR) region and n(k) ∼ k −2 in the ultraviolet (UV) seen for some turbulent, 2D, scalar BEC systems [31,32,42]. Decomposing the kinetic energy into its constituent parts, we see that the incompressible contribution dominates in the IR and is responsible for the change in scaling to k −4 in this region. This incompressible contribution is associated with the vortices in the system [24]. At large k, the spectrum is dominated by the compressible and quantum pressure contributions exhibiting the weak-wave-turbulence scaling k −2 . This scaling of the energy is qualitatively insensitive to variations in γ.
Next, we consider the time-dependent properties of the turbulent dynamics. For this purpose, we will concentrate on the correlation functions for the spin and mass parts of the pseudospinor order parameter. For a homogeneous  n(k) n q (k) n i (k) n c (k) k 2 k 4 Fig. 3: Occupation numbers for different fractions of the system with γ = 0.6 at t = 5 × 10 4 ξ 2 d : single particle spectrum for quantum pressure (purple diamonds), incompressible (red diamonds) and compressible (blue diamonds) contributions. The total occupation number (black diamonds) for the single particle spectrum is obtained by summing the corresponding fractions from each condensate component. The single particle spectrum obeys a k −2 scaling (dotted line) in the ultra-violet and a k −4 scaling (dashed line) in the infrared regions. turbulent system these are defined, respectively, as [43] G Φ (r, t) = 2 where · denotes ensemble averaging. Here, the matrix where Q xx = Re{ψ * 1 ψ 2 } and Q xy = Im{ψ * 1 ψ 2 }, is associated with spin ordering in the system, while α = −2ψ 1 ψ 2 is an alignment parameter. Exploiting the fact that our turbulent system is homogeneous, we can replace ensemble averages with spatial averages. The spin correlation function is then equivalently defined as [43] where dΩ r denotes angular integration. We perform the same averaging for the superfluid correlation function. In fig. 4(a) we plot the results for the spin correlation function at different times for γ = 0. 6 the system. We verify the same behaviour for the masscorrelation function. From these correlation functions we obtain the correlation length, L δ (t) for δ = {Φ, Θ}, which we take as the distance at which the corresponding correlation function decays to a quarter of its value at r = 0: The correlation functions are said to exhibit dynamical scaling when their form at different times remains self similar. This means that they collapse to a universal, time-independent function when scaled by the correlation lengths, i.e. H δ (r) = G δ (r/L δ (t), t). The inset in fig. 4 shows this collapse of the spin correlation function in our system, indicating that G Φ (r, t) does indeed exhibit dynamical scaling. We again verify the same behaviour for G Θ (r, t). Grid Random Random w/ noise t 1/2 t 1/5 Fig. 5: Mean vortex distance in a scalar BEC for three different initial conditions using the same parameters as in ref. [32]. The t 1/5 early-time as well as the t 1/2 late-time scaling regimes are recovered. Fig. 4(b) shows both correlation lengths L Φ,Θ (t) as a function of time for γ = 0.3, γ = 0.6 and γ = 0.8. After the initial evolution the temporal scaling of the correlation lengths becomes universal for all values of γ. However, the effect of γ is apparent in the early time evolution where a larger γ leads to a faster growth of the correlation lengths. This is indicative of a difference in the decay rate of the vortices in the early-time dynamics.
We can investigate this behaviour by considering the total number of vortices in the system as a function of time. We extract the mean distance between vortices as where N vort is the total number of vortices in the system. As a point of reference, in a scalar BEC initially containing a large number of vortices, d ∼ t β [32], where β characterises the vortex annihilation rate. In particular, after some (possibly short) period of evolution, a β = 1/5 scaling is observed. For comparison, we have indicated this theoretically expected scaling in fig. 4(b) for our two-component BEC. At late times, a β = 1/2 scaling appears in the scalar BEC, whose onset is delayed if the initial vortex distribution is highly clustered [32]. In fig. 5 we reproduce this late-time scaling using the parameters of ref. [32] for an initial grid of elementary vortices analogous to our two-component initial state, as well as for a random vortex distribution with and without noise added to the energy spectrum. In all cases we recover both the t 1/5 scaling after initial evolution and the t 1/2 late-time scaling, indicating that this behavior is robust and qualitatively insensitive to details of the initial condition.
Motivated by this previous work, we perform a similar analysis to establish how these results extend to a two-component system with HQVs and how the vortex annihilation rate depends on γ. We focus on the early vortex evolution, where fig. 4(b) suggests that the γ- dependence is significant. Fig. 6 shows N vort as a function of time for three different values of γ. For γ = 0.7 and 0.9, a new scaling regime emerges at early times (2.5 × 10 2 ξ 2 d t 2.5 × 10 3 ξ 2 d ), where N vort (t) decays as t −1 (γ = 0.7) and t −1.5 (γ = 0.9). For t 2.5 × 10 3 ξ 2 d , N vort (t) approaches a universal t −2/5 scaling corresponding to d ∼ t 1/5 , similar to the scalar BEC also shown. These results imply a better agreement with the theoretical t 1/5 scaling than indicated from the correlation lengths [ fig. 4(b)]. This suggests that although their growth is driven by vortex annihilation, the length scales L Φ,Θ (t) are not fully equivalent to d (t). The region of interest in fig. 6 only extends up to t = 5×10 4 ξ 2 d and we therefore expect a universal transition to d ∼ t 1/2 at times extending beyond the time interval of our simulations.
Previous work has demonstrated that, for a sufficiently high γ 0.6, a dipole consisting of HQVs with opposite phase winding in the same component will shrink in size as the vortices move toward one another and annihilate [17]. We therefore attribute the different scaling regime at early times, when the mean inter-vortex separation is small, to this behaviour. This is further supported by the fact that we do not see such scaling for γ 0.6, where such rapid annihilation rate is not prevalent. Within that range of values for γ, the vortex dynamics begins to recover the behavior observed in a scalar BEC.
We can model the vortex decay rate by a kinetic-like equation of the form where η > 1. The dependence of N vort on the right-hand side of the equation indicates that the decay rate is a function of the number of vortices that are involved in facilitating the annihilation. Using this simple model, we can derive temporal scaling of the total vortex number as [44] N vort ∼ t −2/z , where z = −2(1 − η). We note that an exponent of z = 2 corresponds to a two-body collision process whereas z = 5 corresponds to three-body collisions [32]. In fig. 7, we quantify the γ dependence of the early-time scaling by considering the exponent z in the region 2.5 × 10 2 ξ 2 d < t < 2.5 × 10 3 ξ 2 d . We see a rapid decrease of the exponent after γ > 0.6, when the more rapid annihilation becomes prevalent. The observed decrease in the value of z with γ in our simulations signals an additional interaction effect not present in the scalar system.
Conclusions. -We have investigated spatial and temporal aspects of a decaying turbulent two-component BEC containing HQVs. The occupation-number spectrum is found to show a scaling behaviour consistent with similar results for a scalar BEC across a wide range of values of the inter-component interaction strength.
However, we find that a new interaction-dependent scaling regime appears in the temporal properties of the massand spin-correlation functions, as well as the mean intervortex separation. For large values of the relative intercomponent interaction strength, γ 0.6, these exhibit a γ-dependent scaling that is markedly different from the universal behavior, which conforms to that of a scalar BEC at a similar stage of time evolution. Modelling the total vortex number using a simple kinetic equation, we have found that this early-time decay rate for high γ cannot be explained by simple two-or three-body collisions. The observed enhanced vortex decay rate at early times for large γ may be due to the role played by an additional inter-vortex force that arises between vortices in the same component. The results suggest that this force is shortrange and appears in addition to the well-known 1/R intervortex force. This latter force appears to dominate once the vortex density drops significantly following the rapid vortex annihilations occurring at early times. * * * The results presented were obtained using the High Performance Computing Cluster supported by the Research and Specialist Computing Support service at the University of East Anglia.