Thermal quenches in the stochastic Gross-Pitaevskii equation: Morphology of the vortex network

Three-dimensional complex field theories with continuous symmetry breaking are used to describe phase transitions in a host of physical systems including superfluid ⁴He [1–3], type-II superconductors [3,4], nematic liquid crystals [5], magnetic materials [6], weakly interacting bosons [7], and the early universe [8,9]. Such phase transitions lead to the formation of topological defects. In recent years, an impressive theoretical, numerical and experimental effort has been devoted to the measurement of the density of defects after slow quenches through second-order phase transitions [10]. However, very little is known about the size distribution, geometric properties and spatial organisation of the defects inherited from fast and slow quenches. Such questions have been addressed in spin models, and mostly in two dimensions [11,12].

Of particular interest are linear defects, be them vortices, disclinations, or cosmic strings supported by 3d field theories with global U(1) symmetry that capture, in different limits, the systems mentioned in the first paragraph. In this Letter we study the dynamics engendered by the stochastic Gross-Pitaevskii equation for weakly interacting bosons [13] from the point of view of the vortex tangle. We show that the number density of vortex lengths is the one of fully packed models at high T and we fully characterise its time-dependence after sudden quenches. (Relativistic or underdamped extensions yield similar results apart from short-time differences [14].)

The stochastic dynamics of the space-time–dependent complex scalar field $\psi({\mathbf x},t)$ are ruled by the Langevin equation

$$\begin{equation} (2 i \mu - \gamma_{\text{L}}) \dot{\psi} = - \nabla^2 \psi + g(|\psi|^2 - \rho) \psi - \sqrt{\gamma_{\text{L}} T} \xi. \end{equation} \tag{ 1 }$$

The complex noise ξ has Gaussian statistics with zero mean, $\langle \xi_a({\mathbf x},t) \rangle =0$, and $\langle \xi_a({\mathbf x},t) \xi_b({\mathbf x}', t') \rangle = \delta_{ab} \delta({\mathbf x}-{\mathbf x}') \delta(t-t')$ for $a,b=1,2$ the real and imaginary components. g and ρ are real parameters in the Mexican-hat potential energy with degenerate minima at $|\psi|^2= \rho$, the zero temperature equilibrium density. $\gamma_{\text{L}}$ is a friction coefficient. In the following we measure the order parameter in units of $\sqrt{\rho}$, and space and time in units of the T = 0 mean-field equilibrium correlation length, $x_0 \equiv 1/\sqrt{g \rho}$, and correlation time, $t_0 \equiv \gamma_{\text{L}} / x_0^2$. The chemical potential μ, temperature T, and Langevin noise ξ are measured in units of $t_0 / x_0^2$, $x_0^{-2}$, and $\sqrt{\rho / t_0}$, respectively (we set the Boltzmann constant to one hereafter). We place the fields in a three-spatial-dimensional lattice and we solve the stochastic equation numerically with the lowest-order Runge-Kutta method, periodic boundary conditions (PBC), and parameters $g = 1$, $\mu = 0.5$ and $\gamma_{\text{L}} = 1$. We show data obtained using discretization steps $\Delta x = 1$ and $\Delta t = 0.01$, averaged over at least 100 noise realisations. We checked that the equilibrium and dynamic results are not modified using smaller $\Delta x$ and $\Delta t$ when studied as functions of $T/T_{\text{c}}$ (the critical temperature $T_{\text{c}}$ itself depends on $\Delta x$).

Vortices are centred at points where the field vanishes and its phase changes by $2\pi n$, with n a non-vanishing integer, along a closed contour around them. A vortex element threads a square plaquette when the phase of the field on its four vertices changes by $\pm 2 \pi$. When more than one vortex element enter (and exit) a unit cell, we use the stochastic scheme to recombine them. The length l of a vortex is the number of elements along it times $\Delta x$. (We will discuss the effects of other reconnection rules in [14], see also [15,16].) For PBC the lines are closed and can wrap around the system.

Equation (1) takes the system to Gibbs-Boltzmann equilibrium and captures a second-order thermodynamic phase transition at $T_{\text{c}}= 2.26 \pm 0.02$ for the parameters chosen. The critical temperature can be estimated from the analysis of the average over the noise and/or initial conditions of the long-time limit of

$$\begin{equation} \overline \psi \equiv \frac{1}{L^3} \ \left| \sum\limits_{i} \psi_{i} \right|, \end{equation} \tag{ 2 }$$

where the sum runs over all sites on the 3d lattice. This quantity acts as an order parameter vanishing above $T_{\text{c}}$ and taking a non-zero value below $T_{\text{c}}$, see fig. 1. As in the study of usual second-order phase transitions, we computed the Binder parameter and the linear susceptibility associated to it, as well as the specific heat and other quantities. Consistently, all these measurements yield the value of $T_{\text{c}}$ reported above. Moreover, the critical exponents estimated numerically [14] are very close to the ones expected from the -expansion and Monte Carlo simulations of the 3d XY model [17,18].

Zoom In Zoom Out Reset image size

A typical equilibrium state at very high temperature is such that two vortex elements enter and exist each unit cubic cell. The number density of vortex lengths normalised by the sample volume, $N(l)/L^3$, is shown in fig. 2 for various L. The data crossover from an algebraic decay with exponent 5/2, as for Gaussian random walks [19], to a weaker decay with exponent 1, the result for fully-packed oriented loop models with the two-in two-out constraint [20]. Due to the Gaussian statistics short loops feel the size of the simulation box at scales ${\mathcal O}(L^2)$, as confirmed by the collapse of lN vs. $l/L^2$ In short,

$$\begin{eqnarray} N(l)/L^3 \simeq \left\{ \begin{array}{l@{\qquad}l} l^{-5/2}, & l \ll L^2,\\ l^{-1}\ L^{-3}, & l \gg L^2, \end{array} \right. \end{eqnarray} \tag{ 3 }$$

at $T \gg T_{\text{c}}$. Gaussian statistics were found in a Z₃ model for the initial condition of cosmological models [21] and in random optical fields [22,23] observed in laser speckels and modelled with random wave superpositions. The crossover to length-scales longer than L² was not discussed in these works.

Zoom In Zoom Out Reset image size

The detailed geometric analysis of line ensembles has been very successful in, e.g., polymer science [19]. At a critical point, be it thermodynamic or geometric, the critical objects are fractal. Their Hausdorff dimension D relates the linear length along the loop, l ($\gg l_0$ a microscopic length-scale), and the radius of the smallest sphere that contains the loop, R, as $l \simeq R^{D}$. In the thermodynamic limit,

$$\begin{equation} N(l)/L^3 \simeq l^{-\alpha_{\text{L}}}\ e^{- l m_{\text{L}}}, \end{equation} \tag{ 4 }$$

with $m_{\text{L}}$ the line tension and $\alpha_{\text{L}}$ the "Fisher exponent". (This form should be corrected to capture finite-size corrections.) At criticality $m_{\text{L}} \simeq |T-T_{\text{L}}|^{\beta_{\text{L}}}$ with $\beta_{\text{L}}$ another characteristic exponent. D and $\alpha_{\text{L}}$ are linked by $\alpha_{\text{L}} = 1+d/D$.

In equilibrium, the vortices undergo a geometric transition between a "localised" phase with finite loops only, and an "extended" phase in which the length of a finite fraction of lines diverges in the thermodynamic limit. The geometric threshold falls within the thermodynamically ordered phase, at $T_{\text{L}} \simeq 0.98 T_{\text{c}}$, where m_L vanishes. In fig. 3 we display $N(l)/L^3$ at $T_{\text{L}}$. An algebraic data fit yields $\alpha_{\text{L}} = 2.17(2)$ and the dashed line uses $\alpha_{\text{L}} = 2.17$, accordingly $D\simeq 2.56$, a value to be associated with self-seeking random walks. Following the arguments in [24] the value of $\alpha_{\text{L}}$ should coincide with the one for interacting random walks in 3d studied in [25], indeed, a very close numerical value was measured in this work, $\alpha_{\text{L}}=2.18(3)$. In the inset, the variation of $m_{\text{L}}$ approaching $T_{\text{L}}$ suggests $\beta_{\text{L}} \simeq 1.7$.

Zoom In Zoom Out Reset image size

We confirmed $T_{\text{L}} <T_{\text{c}}$ with the analysis of other observables. We defined the size of a loop as the maximal side of the parallelepiped covering the vortex in the three Cartesian directions, i.e., the length that the string would have after smoothing out all small scale irregularities. We then counted the loops with larger size than the system size. Besides, we counted the non-contractible loops, i.e., those with non-vanishing winding number in at least one of the three spatial directions. The two numbers behave as order parameters for the geometric transition. Figure 4 displays the temperature dependence of the averaged $N_{\text{non-cont}}$ together with $T_{\text{L}}$ and $T_{\text{c}}$ shown with vertical dotted lines.

Zoom In Zoom Out Reset image size

We conclude that the geometric threshold differs from the thermodynamic critical point. The same fact was found in the 3d XY model [15] and the O(2) field theory [16]. These findings refute claims of the coincidence of the two transitions in the context of cosmology [26,27], in the field of superfluidity and type-II superconductivity [28,29], and in general [30]. See [24] for the field theory of the geometric transition.

We now turn to the dynamics after instantaneous quenches. A system prepared in equilibrium at $T\gg T_{\text{c}}$ is subsequently evolved at T = 0. For a sufficiently large L the spatial averaged initial field is practically zero and the distribution of vortex lengths is very close to the one in fig. 2. After a transient, the system enters a dynamic scaling regime in which the space-time correlation, $C(r,t) \equiv \langle \psi^*({\mathbf x}, t) \psi({\mathbf x}', t) \rangle$ with $r \equiv |{\mathbf x} -{\mathbf x}'|$, or the dynamic structure factor, scale as $C(r,t) \simeq\ f(r/\xi_{\text{d}}(t))$, or $S(k,t) \simeq \xi^d_{\text{d}}(t) \Phi(k \xi_{\text{d}}(t))$, with a growing length $\xi_{\text{d}}(t) = \gamma_{\text{v}} t^{1/2}$, see fig. 5. This time-dependence was predicted in [31] for continuous spin models with non-conserved order parameter. Consequently, the vortex density decays as $t^{-1}$, see the inset. (We defined $\rho_{\text{vortex}}$ as the number of plaquettes with non-vanishing flux divided by their total number. Its "value" depends on $\Delta x$ but its time-dependence does not.) Finite-size effects, once analysed via $\xi_{\text{d}}/L$ against $t/L^2$ yield a first saturation to $\xi_{\text{d}} \simeq 0.7 \, L$ at $t\simeq 0.5L^2$, before the final approach to $\xi_{\text{d}} = L$ at $t \simeq L^2$. Vortices inhibit long-range order in the entire space but, as soon as they disappear at the end of the first plateau, the system rapidly enters the state with large phase correlation. This is the reason why $\xi_{\text{d}}$ increases suddenly at $t\simeq 0.5L^2$ and quickly approaches the system size L.

Zoom In Zoom Out Reset image size

We focus now on the time-evolving number of rings with length l. We first describe the numerical data and we later give functional forms for the various regimes.

Early regime, fig. 6: At t = 1, N is still very close to the initial one. Soon after, a short-length regime with positive slope starts to develop. To the right of the maximum, the algebraic decay is no longer given by the Gaussian 5/2 but by a value that is numerically indistinguishable from the one of the critical threshold, $\alpha_{\text{L}}= 2.17$. The time needed to reach this algebraic behaviour is $t_{\text{p}} \simeq 7$ for L = 100. The precise determination of the variation of $t_{\text{p}}$ with L is a hard task [32,33] and goes beyond the scope of this letter. The weight of the longest length-scales remains the time-independent $l^{-1}$. The crossover between the two power laws occurs at a shorter l than for the initial state because the weight of the shorter loops has decreased by a factor of time.

Zoom In Zoom Out Reset image size

Dynamic scaling regime, fig. 7: Three scales can be identified for L = 100 (above) and L = 500 (below). At short lengths the curves are well represented by a linear dependence (dotted lines) until a maximum at a scale of the order of the growing correlation length $\xi_{\text{d}}$. Beyond it a regime with weight $l^{-2.17}$ is seen in both panels (data fits yield $\alpha_L = 2.17(9)$ and $\alpha_L=2.17(5)$). For the smaller system size, it is progressively erased and at t = 30 only the final $l^{-1}$ decay remains. For the larger system size, the curves are parallelly translated down. The last $l^{-1}$ regime is also seen. In the inset to the lower panel the scaling t²N vs. $l/t^{1/2}$ describes the first two regimes quite accurately.

Zoom In Zoom Out Reset image size

Late epochs, fig. 8: At very long times just a few non-contractible long vortices remain, until they break and finally disappear, letting the system phase order. The number of loops has a peak at $l\simeq L$. In the inset, the average number of loops, $N_{\text{loop}}$, and non-contractible loops, $N_{\text{non-cont}}$, as functions of time. In equilibrium these quantities vanish and, accordingly, the dynamic values approach zero.

Zoom In Zoom Out Reset image size

A simple argument along the lines of the one used in [11,12] to deal with the statistics of clusters in spin models leads to an analytic —though approximate— expression for N.

The first remark is that the double algebraic decay of N at t = 0 or $t_{\text{p}}$ is well approximated by

$$\begin{equation} N(l, t) \simeq l^{-\alpha_{\text{L}}} \ [c_1^{n}(t) + c_2^n(t) l^{(\alpha_{\text{L}}-1) n}]^{1/n}. \end{equation} \tag{ 5 }$$

A sharp crossover between the two power laws is obtained for large n (n = 6 is sufficient, see fig. 6). The crossover takes place at $l^*(t) \simeq [c_1(t)/c_2(t)]^{1/(\alpha_{\text{L}}-1)}$. If $\alpha_{\text{L}}=5/2$, as for the initial state, $c_1(0) \propto L^3$ and $c_2(0)$ a finite constant ensure $l^*(0) \simeq L^2$ and the limit forms in eq. (3). At $t_{\text{p}}$ the weight of the loops with $l > l^*(t_{\text{p}})$ is not modified with respect to the initial one, and $c_2(t_{\text{p}}) = c_2$. The total number of loops diminishes in time but remains $O(L^3)$ until the very late epochs, see the inset in the upper panel in fig. 7. For N as in (5), to leading order in L, and ignoring constants, $N_{\text{loop}}(t) = \int \text{d}l \, N(l,t) \simeq c_1(t) \int_{l_0}^{l^*}\text{d}l \, l^{-\alpha_{\text{L}}} + c_2(t) \int_{l^*}^{L^3}\text{d}l\, l^{-1} \simeq c_1(t) l_0^{1-\alpha_{\text{L}}}$ that scales as L³ if $c_1(t) =L^3 \overline c_1(t_{\text{p}})$. According to fig. 6, $l^*(t_{\text{p}}) < l^*(0)$.

We assume that after the transient $t_{\text{p}}$, the length of each vortex is reduced at the same rate as the dynamic correlation length $\xi_d$ grows

$$\begin{equation} l(t, l_{\text{p}}) \simeq \gamma_{\text{v}} \sqrt{t_{\text{a}} - t}, \end{equation} \tag{ 6 }$$

with $\gamma_{\text{v}}$ a parameter, $t_{\text{a}} = t_{\text{p}} + l_{\text{p}}^2 / \gamma_{\text{v}}^2$ the annihilation time at which $l(t_{\text{a}}) =0$, and $l_{\text{p}}$ the length of the vortex ring at $t_{\text{p}}$. We suppose that the vortices are sufficiently long and far apart that they evolve independently of each other. Neglecting the fact that they break up and disappear in the course of evolution we use $N(l,t) \simeq \int \text{d}l_{\text{p}} \, N(l_{\text{p}}, t_{\text{p}}) \, \delta (l - l(t, l_{\text{p}}))$ to obtain

$$\begin{equation} N(l,t) \simeq \frac{l\ N(\sqrt{\gamma_{\text{v}}^2 (t-t_{\text{p}}) + l^2}, t_{\text{p}})}{\left[\gamma_{\text{v}}^2 (t-t_{\text{p}}) + l^2\right]^{1/2}} \; . \end{equation} \tag{ 7 }$$

Let us first focus on long length scales. From the numerator one estimates a crossover at a dynamic length $l^*(t)$, that is advected towards smaller scales as time evolves,

$$\begin{eqnarray} &&{l^*}(t) \approx \sqrt{(c_1(t_{\text{p}})/c_2)^{1/(\alpha_{\text{L}}-1)}- \gamma_{\text{v}}^2 (t-t_{\text{p}})}, \end{eqnarray} \tag{ 8 }$$

as observed in the numerical data. We recover

$$\begin{equation} N(l,t) \simeq c_2 \ l^{-1}, \qquad l \gg l^*(t), \end{equation} \tag{ 9 }$$

independently of time. Instead, for $l < l^*(t)$ we need to correct (7) to take into account the annihilation of vortices with short length that implies $N_{\text{loop}}(t) \simeq t^{-\zeta}$, see the inset in the upper panel in fig. 7, and the consequent reduction of the averaged length size $\langle l \rangle \simeq t^{-\zeta+1/2}$. We enforce this scaling heuristically, by simply multiplying eq. (7) by $(\gamma_{\text{v}}^2 t)^{-\zeta}$. Proceeding in this way, and taking $t\gg t_{\text{p}}$,

$$\begin{equation} (\gamma_{\text{v}}^2 t)^{\zeta+\alpha_{\text{L}}/2}\ N(l,t)\simeq \frac{ c_1(t_{\text{p}})~l/(\gamma_{\text{v}}^2 t)^{1/2} }{\left[ 1+ l^2/(\gamma_{\text{v}}^2 t) \right]^{(1+\alpha_{\text{L}})/2}}. \end{equation} \tag{ 10 }$$

Finally, this regime can also be split in two:

$$\begin{eqnarray} (\gamma_{\text{v}}^2 t)^{\zeta+\frac{\alpha_{\text{L}}}{2}} \ \frac{N(l,t)}{L^3} \simeq\left\{ \begin{array}{ll} [l/(\gamma_{\text{v}}^2 t)^{1/2}], & l \ll \xi_{\text{d}}(t), \\ {[l/(\gamma_{\text{v}}^2 t)^{1/2}]}^{-\alpha_{\text{L}}}, & l \gg \xi_{\text{d}}(t), \end{array} \right. \end{eqnarray} \tag{ 11 }$$

For $\alpha_{\text{L}} \simeq 2.17$ and $\zeta \simeq 1.1$, the exponent on the left-hand side is close to 2, the value used in the inset in fig. 7 (lower panel).

In conclusion, we presented a full characterisation of the vortex network in and out of equilibrium. It should be possible to observe these static and dynamic vortex statistics in experimental systems, by, for example, examining the 3d vortex configurations in the turbulent state of a trapped Bose-Einstein condensate shortly after its expansion [34].

We thank I. Carusotto, J. T. Chalker, P. Comaron, F. Larcher, M. Picco and N. Proukakis for useful discussions. This research was supported in part by NSF under Grant No. PHY11-25915, KAKENHI (Grant No. 26870295), Global COE Program "the Physical Sciences Frontier", the Photon Frontier Network Program, MEXT, Japan, Grant-in-Aid for Scientific Research on Innovative Areas "Fluctuation & Structure" (Grant No. 26103519) from the Ministry of Education, Culture, Sports, Science, and Technology of Japan, the JSPS Core-to-Core program "Non-Equilibrium Dynamics of Soft-Matter and Information.", and the IRSES European Project "SoftActive". LFC is a member of the Institut Universitaire de France.