Anomalous phase ordering of a quenched ferromagnetic superfluid

Coarsening dynamics, the canonical theory of phase ordering following a quench across a symmetry breaking phase transition, is thought to be driven by the annihilation of topological defects. Here we show that this understanding is incomplete. We simulate the dynamics of an isolated spin-1 condensate quenched into the easy-plane ferromagnetic phase and find that the mutual annihilation of spin vortices does not take the system to the equilibrium state. A nonequilibrium background of long wavelength spin waves remain at the Berezinskii-Kosterlitz-Thouless temperature, an order of magnitude hotter than the equilibrium temperature. The coarsening continues through a second much slower scale invariant process with a length scale that grows with time as $t^{1/3}$. This second regime of coarsening is associated with a turbulent cascade that transports spin wave energy from low to high wavevectors, bringing about the the eventual equilibrium state. Because the relevant spin-waves are noninteracting, the cascade occurs through a dynamic coupling to other degrees of freedom of the system. Strongly coupling the system to a reservoir destroys the second regime of coarsening, allowing the system to thermalise following the annihilation of vortices. Our results, which provide a connection between the currently disparate fields of phase ordering dynamics and wave turbulence, offer a paradigm for understanding phase ordering dynamics beyond the current framework built around topological defects.


Introduction
Quenching a system across a continuous phase transition from a high to low symmetry phase causes the system to spontaneously break symmetry. Immediately after the quench causally disconnected regions of the system will break symmetry independently, resulting in the formation of domains with independent order parameter orientation. The subsequent growth of these domains toward the global equilibrium state is known as coarsening dynamics. Although the microscopic details of coarsening are usually extremely complicated, at a macroscopic level a much simpler scaling regime can emerge for large average domain size L. Spatial correlations of the order parameter at different times t then collapse onto a single curve when rescaled by L, and the domains grow as L ∼ t 1/η with the scaling exponent η determined by the dynamic universality class [1]. Such universal dynamics has been explored in a vast variety of systems, ranging from the early universe [2] to superfluid formation [3] to opinion spreading in sociology [4]. When the quench produces topological defects, the decay of these defects has long been thought to provide a unifying framework for understanding the coarsening [1].
In this work we study the easy-plane ferromagnetic ordering of a homogeneous quasi-2D spin-1 condensate after all vortices have annihilated. Remarkably, we find that the annihilation of vortices does not take the system to the equilibrium state. Instead, a nonequilibrium background of spin waves remain at the Berezinskii-Kosterlitz-Thouless (BKT) temperature, an order of magnitude hotter than the eventual equilibrium temperature. The coarsening then continues via spin wave energy transport from low to high wavevectors that displays features of a novel turbulent cascade, relevant to the emerging area of spin turbulence (e.g. see [13,[29][30][31][32]). We argue that the nonlinear transverse spin wave dynamics arises from a dynamic coupling to interacting axial spin degrees of freedom. Order parameter correlations show dynamic scale invariance during the spin wave coarsening, with a length scale that grows as t 1/3 . This scaling is distinct from that during the vortex driven coarsening, showing that there are two renormalisation group fixed points affecting the phase ordering of this system. Strongly coupling the system to a reservoir of energy and particles destroys the second scaling regime, allowing the system to thermalise following the annihilation of vortices. Our results give new insights into the phase ordering dynamics of isolated systems and provide a potential profitable connection between phase ordering and wave turbulence.

Background
A spin-1 condensate can be described by three interacting classical fields ψ m for condensates in the three spin components with spin projections m = −1, 0, 1. The quasi-2D Hamiltonian density within a uniform trap [11] is [33][34][35][36], where M is the atom mass, n = m |ψ m | 2 is the areal density, g n is the quasi-2D density interaction strength and H s encompasses the spin dependent terms, The first term in H s is the spin interaction energy, with spin density F = mm ψ * m f mm ψ m /n 2 for spin-1 matrices ( f x , f y , f z ) ≡ f, and quasi-2D spin interaction strength g s . The sign of g s determines whether the interactions are ferromagnetic (g s < 0), which occurs in 87 Rb [37], or antiferromagnetic (g s > 0), which occurs in 23 Na [35]. Here we consider the ferromagnetic case. The second term in H s is a quadratic Zeeman splitting of the spin components, which can be induced using either DC magnetic fields or AC microwave stark shifts [15,38]. A linear Zeeman term pnF z can also be included, but conservation of nF z means this term does not affect the system dynamics and can be removed via the unitary transformation e −ipmt/ħ h ψ m → ψ m . The quasi-2D regime is obtained from a 3D system by tightly confining the system in one direction and integrating over the resulting spatial profile along that direction [5,36].
The relative strength of the two terms in H s produces a rich phase diagram, from which a variety of quenches can be explored. The zero temperature mean field phase diagram for ferromagnetic interactions and with q > 0 is shown in Figure 1(a). A quantum critical point at q = q 0 ≡ 2|g s |n 0 (n 0 is the mean condensate density) separates the unmagnetised polar phase (all atoms in the m = 0 condensate) from the easy-plane ferromagnetic phase with spin order parameter F ⊥ ≡ (F x , F y ) (for quantization along F z ). The order parameter manifold of F ⊥ is SO(2) with transverse spin vortices as topological defects. These vortices consist of a positive or negative phase winding of the transverse spin angle θ (tan θ = F y /F x ), and can only decay via the mutual annihilation of two vortices of opposite sign. Vortices with negative phase winding are also termed antivortices. The energy scale q 0 defines a time scale t s ≡ ħ h/q 0 and the spin healing length ξ s ≡ ħ h/ M q 0 .

Quench dynamics: anomalous phase ordering
We simulate the condensate dynamics following an instantaneous quench of the quadratic Zeeman energy from deep in the polar phase to q = 0.3q 0 in the easy-plane ferromagnetic phase, see Fig. 1(a),(b). Symmetry breaking and the production of transverse spin vortices following Figure 1: (a) A spin-1 ferromagnetic condensate is unmagnetised ("polar") for q > q 0 and magnetises in the transverse plane ("easy-plane") for 0 < q < q 0 . The point q = q 0 is a quantum critical point (QCP). We explore the ordering of transverse spin following a quench from q q 0 to q = 0.3q 0 . (b) Coarsening of transverse spin domains [colormap shown in (a)]. This is associated with collisions between transverse spin vortices (red triangles) and antivortices (black circles), which can then mutually annihilate, resulting in a growing intervortex spacing L v (red bar). A second length scale L sw giving the thermal wavelength of spin waves grows much more slowly (black bar), such that the transverse spin remains out of equilibrium long after all transverse spin vortices have annihilated. The central time axis quantifies the first stage of vortex driven coarsening, the time of last vortex annihilation, and the subsequent stage of spin wave thermalisation. (c) Spatial correlations of transverse spin at different times, showing that long after all vortices have annihilated the correlations still decay more rapidly than the equilibrium prediction (upper black dashed line). such a quench have been observed in experiments with 87 Rb [5]. Conservative dynamics of our system is simulated by numerically integrating the three coupled Gross-Pitaevskii equations (GPEs) obtained from Eq. (1) [14], Further numerical details are described in Appendix A.1. A homogeneous system can be realised in experiments using a flat bottomed trap [11,39]. We quantify order in the system by spatial correlations of F ⊥ , where angular brackets denote an ensemble average (see Appendix A.2). Figure 1(c) shows the evolving correlation function Eq. (4). For times 10 2 t s t 10 3 t s , the growth of order is scale invariant and driven by the mutual annihilation of transverse spin vortices of opposite sign, see Fig. 1(b), with correlations decaying to zero at a length scale on the order of the intervortex spacing. This vortex driven coarsening has been described in previous work [19,20]. We find that all vortices have annihilated by a time t ≈ 2.8 × 10 3 t s , after which correlations extend to the boundary. The correlations can then be compared to the equilibrium (thermalised) prediction [36,40], Here T BKT = πK/2k B is the BKT temperature associated with the unbinding of transverse spin vortices [41], with K = ħ h 2 n 0 (1 − q/q 0 )/2M the spin wave stiffness and k B Boltzmann's constant. The equilibrium temperature T eq of our microcanonical system is calculated by equipartitioning the energy liberated by the quench amongst all collective modes of the system [20], which gives ν ≈ 0.011. This equilibrium prediction is shown in Fig. 1(c). Surprisingly, even after very long simulation times t = 10 5 t s , correlations of transverse spin only agree with the equilibrium prediction for length scales r 5ξ s . For larger length scales the correlations decay more rapidly. This absence of equilibrium following the annihilation of topological defects is not predicted by the current theory of coarsening dynamics [1].

Spin wave energy transport driving phase ordering
To identify the origin of the unexpectedly slow ordering displayed in Fig. 1(c) we look at the distribution of energy in the gradient of the transverse spin angle ∇θ (this vector field is proportional to currents of F z magnetization [42]). We firstly perform a Helmholtz decomposition The first contribution v i , known as the incompressible field, arises from vortex excitations while the second contribution v c , known as the compressible field, arises from transverse spin wave excitations. The spectral energies of the incompressible and compressible fields are given by, whereṽ µ (k) = l −1 d 2 r v µ (r)e −ik·r is the Fourier transform of v µ (r) and angular brackets denote an ensemble average (see Appendix A.2). The evolving spectral energies ε µ (k, t) are shown in Fig. 2(a),(b). The incompressible spectral energy, Fig. 2(a), shows a k −2 decay when vortices are present, in agreement with the infrared (ξ s k < 1) scaling of a distribution of quantum vortices [43,44]. Once all vortices have annihilated the spectral energy drops abruptly. In comparison, the compressible spectral energy, Fig. 2(b), shows nonequilibrium features across the duration of the simulation. The initial condition of our simulation results in a flat high energy distribution ε c (k, 0) ≈ 200k B T eq . For times t 10 3 t s , the compressible spectral energy shows three approximate regimes, We have introduced the evolving wavevectors k lw (t) and k eq (t) to signify the boundaries between the three regimes of ε c (k, t). The spectral energy ε lw (k, t) is the long wavelength portion of ε c (k, t), with energy per mode ε lw (k, t) ≈ 10k B T eq ≈ k B T BKT /2 in the wavevector window The evolving compressible field spectral energy ε c (k, t) (solid lines) shows three regions: a persistent high temperature long wavelength region with a temperature approximately equal to T BKT ; a steep region with an approximate k −4 scaling; and a short wavelength thermal region. The spectral energy of F z excitations ε F z (k, t) (dashed lines) closely follows ε c (k, t) for times t 400t s . The interacting F z fluctuations can mediate the thermalisation of the transverse spin waves. (c) The combined spectral energy of transverse and axial spin waves, E sw (t), is decomposed into a low wavevector portion E low (t), which decreases in time, and a high wavevector portion E high (t), which increases in time, consistent with a cascade of energy from low to high wavevectors. The total spin wave energy E sw (t) also decreases in time indicating energy flow away from spin waves.
considered. This nonequilibrium temperature, being approximately at the BKT temperature, corresponds to the typical energy of a single transverse spin vortex [41,45], and may be a remnant of interactions between spin waves and vortices during the vortex driven coarsen-ing. For k > k lw ε c (k, t) decays steeply as k −α with an exponent α ≈ 4 until the equilibrium distribution ε c (k, t) = k B T eq /2 is reached at a wavevector k eq (t). The structure of ε c (k, t) is suggestive of a turbulent cascade, with a high temperature long wavelength energy source cascading to a short wavelength thermal field. We provide further evidence of this shortly. With no vortices present, the persistent nonequilibrium features of ε c (k, t) must be responsible for the anomalously slow ordering that we observe in Fig. 1(c).
The observed dynamics of epsil on c (k, t) necessarily involves nonlinear interactions. The field conjugate to the transverse spin phase θ , i.e. the generator of rotations of θ , is nF z , leading to dynamic coupling between transverse spin waves and the axial spin waves of F z [20,36]. The interacting axial spin waves can therefore mediate transverse spin wave interactions. (Transverse spin waves can also interact via nematic N zz = |psi −1 | 2 + |ψ 1 | 2 fluctuations [36], however we find N zz fluctuations are thermalised after the vortices have annihilated and are therefore relatively small.) Expanding the system Hamiltonian to quadratic order in F z and N zz fluctuations [36] gives the spectral energy of axial spin fluctuations, whereF z (k) ≡ l −1 d 2 r F z (r)e ik·r and angular brackets denote an ensemble average (see Appendix A.2). For times t 400t s the spectral energy ε F z (k, t) closely follows ε c (k, t), see Fig. 2(b), indicating that the dynamics of the two spectra are coupled and in equilibrium with each other. The nonlinear interactions of axial spin waves can allow the redistribution of energy in ε F z (k, t) and then dynamic coupling to transverse spin waves can actuate the same effect in ε c (k, t).
To provide further evidence for the presence of an energy cascade in Fig. 2(b) we decompose the total spin wave energy into a low wavevector portion and a high wavevector portion where we choose k mid = 0.5ξ −1 s . In Fig. 2(c) we plot the energy changes ∆E(t) ≡ E(t) − E(10 5 t s ) of these three quantities for times after all vortices have annihilated. The energy E low decreases in time while E high increases, consistent with an energy cascade from k < k mid to k ≥ k mid . There is also a net decrease in the total spin wave energy E sw , showing that energy is also lost from the spin wave excitations, either to other quadratic excitations [36,[46][47][48] or to excitations beyond quadratic order. In principle, one could solve for the dynamics of these additional excitations to obtain effective spin wave dynamics that would transport energy from low to high wavevectors. Figure 2(b) shows that the spin wave energy transport is associated with an approximate k −4 scaling of ε c (k, t) and ε F z (k, t). (Note the spectral energies in most studies of 2D turbulence include a k phase space factor so that the k −4 scaling observed here would normally be described as k −3 scaling.) There are currently no predictions for such a cascade within weak wave turbulence theory [49,50]. To confirm that the energy transport shown in Figs. 2(b),(c) is a turbulent cascade would require showing that the energy transport is local in wavevector space. The length scale L sw (t) grows as t 1/3 for t > 10 3 t s , much slower than the t/ ln(t/t s ) growth of average intervortex spacing L v (t). The largest thermarlised wavelength extracted from ε c (k, t) is 2πk eq (t) −1 , which follows the growth of L sw (t).

A second regime of scale invariance
The robust shape of ε c (k, t) for times t 10 3 t s (see Fig. 2(b)) suggests a regime of scale invariance driven by spin waves, beyond the scale invariant coarsening dynamics driven by vortex annihilation. To explore this we consider the late time dynamics of correlations of transverse spin, Eq. (4), which in a scale invariant regime will evolve as [51,52], for some universal function f and growing length scale L(t). The r −ν correction factor ensures G(r, ∞) ∼ r −ν , consistent with equilibrium. Since ν ≈ 0.011 1, the correction is only significant when G(r) is close to ordered.
The evolving correlation function for times after all vortices have annihilated is shown in the inset to Fig. 3(a). The correlations exhibit a short wavelength ordered portion that grows slowly in time and a nonequilibrium long wavelength portion. The correlation functions collapse onto a single curve after rescaling according to Eq. (12), see Fig. 3(a). We define the rescaling factor L sw (t) by G(L sw , t) = 0.8G(0, t), which follows the boundary between the ordered portion of the correlation function and the nonequilibrium portion. This length scale is governed by spin waves and grows as a power law L sw ∼ t 1/3 for times t 10 3 t s , i.e. times after all vortices have annihilated, see Fig. 3(b). The length scale 2πk eq (t) −1 , where k eq (t) is introduced in Eq. (7) and defined more precisely in Appendix A.3, follows the growth of L sw (t). For comparison, the scale invariance during vortex driven coarsening is associated with the more rapidly growing average intervortex spacing L v (t) (defined in Appendix A.4) [19]. The nonequilibrium portions of the correlation functions in Fig. 3(a) clearly exhibit an additional algebraic decay G(r) ∼ r −0.21−ν ≈ r −0.22 . The value of the decay exponent corresponds to a temperature of T ≈ 0.9T BKT , see Eq. (5), and is consistent with the nonequilbrium temperature of ε lw (k, t) from Eq. (7).

Comparison with open system dynamics
Our analysis so far has considered isolated, energy conserving dynamics. It is of interest to compare our results with open system quench dynamics, where the condensate is coupled to a reservoir of energy and particles. Using a stochastic Gross-Pitaevskii theory (see Appendix A.5), we model a spin-1 condensate strongly coupled to a reservoir with fixed temperature and chemical potential, which we choose such that the equilibriated energy and particle number matches those of the conservative dynamics. We then simulate the same quench as for the isolated system dynamics. Figure 4 Fig. 4(a) are in stark contrast to the results in Fig. 1(c) for the isolated system. Indeed, differences in the two cases are apparent from the evolving spin domains, Fig. 1(b) and Fig. 4(b), with the open system being more ordered in the spaces between vortices. For the large reservoir coupling strength we have used here, spin waves in the open system are able to rapidly thermalise directly with the reservoir rather than via interactions with other spin waves. However, we emphasise that microscopically derived reservoir coupling strengths are much smaller than the value we use here [53], and therefore the isolated system dynamics are a realistic approximation to experiments.
The growing length scales L v and L sw for the open system dynamics, defined as for the conservative dynamics, are shown in Fig. 4(c). The growth of L v in the open system is very similar to the isolated system growth (denoted by L v,ISO in this figure). In the open system, however, there is no second growing length scale, and L sw follows the growth of L v .
The decay of transverse spin correlations for open system dynamics following quenches to different values of q show good agreement with Eq. (5) once all vortices have annihilated, see  Fig. 1(a)]. Spin vortices and antivortices are marked by red triangles and black circles respectively. Comparing with Fig. 1(b), it is clear that the transverse spin is more ordered in the space between vortices in the open system dynamics. (c) The growth of intervortex spacing L v follows the isolated system growth L v,ISO . The length scale L sw follows the growth of intervortex spacing L v , indicating the absence of a second coarsening process. Inset: Algebraic decay exponent ν fit for different q quenches (dots) obtained from single trajectory simulations by fitting to the transverse spin correlation function for 2ξ s ≤ r ≤ 100ξ s and averaging the result across times 5 × 10 3 t s ≤ t ≤ 10 4 t s . The error bars give the standard error of this mean. For q > 0.1q 0 the fitted exponents agree well with the equilibrium prediction from Eq. (5) (solid line). Fig. 4(c) inset. (The temperature and chemical potential for these quenches have been adjusted as a function of q; see Appendix A.5.) The small deviation at the smallest q value might be caused by axial spin fluctuations, which become stronger as q → 0 due to a diminishing energy gap. Indeed, we expect that the physics will be modified in the limit q → 0, since the ground state manifold changes from SO(2) × U(1) to SO(3), resulting in changes in collective mode excitations [47,48] and vortex topology [21].

Conclusion
We have shown that vortex driven coarsening of an isolated easy-plane ferromagnetic spin-1 condensate does not take the system to equilibrium. Instead, a second regime of scale invariant coarsening associated with transport of spin wave energy scales more slowly as t 1/3 . Strongly coupling the system to a reservoir of energy and particles destroys this second coarsening process and equilibrium is reached after the vortex driven coarsening.
The presence of two dynamic scaling regimes in the isolated dynamics shows that there are two renormalisation group fixed points affecting the phase ordering. The first, associated with vortices, has been ascribed to the model E dynamic universality class [19]. The second slower scaling, L sw ∼ t 1/3 , matches that of the scalar model B dynamic universality class [1,54], however the order parameter does not: the scalar model B universality class describes a one component conserved order parameter, whereas the order parameter in our system has two components and is not conserved. There is, however, a second important field in our system that does satisfy the properties of the scalar model B universality class: the conserved nF z field. It could be possible that the dynamics of the nF z field belongs to the scalar model B dynamic universality class, even though this is not the order parameter of the system, and that the dynamic coupling between nF z and θ leads to model B scaling emerging in the correlations of transverse spin. However, model B is also a dissipative model whereas we have shown that strong dissipation destroys the second scaling regime. Hence the second scaling regime might alternatively belong to a currently unidentified dynamic universality class unique to isolated systems. We have shown that this second scaling regime displays features of a spin wave energy cascade, thus identifying a potential connection between the fields of wave turbulence and phase ordering dynamics.
The nonequilibrium background of spin waves that remain after the vortices have annihilated is at a temperature very close to T BKT . These spin waves may have thermalised with the vortex field during vortex driven coarsening, either via scattering off of vortices or via spin wave production after vortex annihilation (see [55]). The absence of interactions once all vortices have annihilated would then leave behind high temperature spin waves, reminiscent of photons decoupling from matter in the early universe to produce the cosmic microwave background. This intriguing process may be ubiquitous in phase ordering systems involving topological defects interacting with collective mode excitations.

Funding information
We acknowledge support from the Marsden Fund of the Royal Society of New Zealand. LW acknowledges support from the University of Otago Postgraduate Publishing Bursary.

A.1 GPE simulations
The GPE simulations are conducted on a 2D square grid with dimensions l × l = 400ξ s ×400ξ s covered by an N × N = 512 × 512 grid of equally spaced points. In experiments in 87 Rb, g n /|g s | ∼ 100 [37]. We use a more modest ratio g n /|g s | = 10, which is sufficient to suppress density fluctuations at the energy scale we are interested in. The mean condensate density is taken to be n 0 = 10 4 ξ −2 s . We evolve our system using a recently developed fourth order symplectic integrator [56] to ensure that energy, atom number and nF z magnetization are conserved effectively. We find that total energy and atom number are conserved to within a factor of 10 −9 across the full simulation time. The total axial magnetization d 2 r n(r)F z (r) remains below 10 −6 n 0 l 2 . We use a time step of 0.02t s for each integration step. The kinetic energy time step is evaluated spectrally using fast Fourier transforms, and we employ periodic boundary conditions. Our initial state is the polar state (ψ 1 , ψ 0 , ψ −1 ) = n 0 (0, 1, 0) + δ, where δ is noise added to Bogoliubov modes on top of the ground state at q = ∞, as in [20], which seeds the symmetry breaking evolution. Noise added this way corresponds to adding on average half a particle per mode according to the truncated Wigner prescription [57]. We then evolve our system using Eq. (3) at a quadratic Zeeman energy q = 0.3q 0 , so that the quench is effectively instantaneous at t = 0.

A.2 Ensemble averaging
Correlations Eq. (4) and spectral energies Eq. (6) and Eq. (8) are computed using an ensemble average of the form,ḡ where u = r, k and g(u) denotes the result of a single simulation trajectory. In the GPE simulations, the ensemble average is over 30 simulation trajectories conducted with independent initial noise. In the open system simulations, the ensemble average is over 10 simulation trajectories. We also averageḡ(u) over azimuthal angles of the coordinate u, such that g(u) →ḡ(u) for u ≡ |u|. Correlation functions are additionally averaged over space, i.e. we replace F ⊥ (0) · F ⊥ (r) by F ⊥ (r ) · F ⊥ (r + r) in Eq. (4) and average over the spatial coordinate r .

A.3 Definition of k eq
The wavevector k eq introduced in Eq. (7) is obtained as follows. We firstly skew the spectrum ε c (k) by multiplying by k. We then define k eq as the position of the local minimum that appears in kε c (k) at the start of the equilibrium portion of the spectral energy. To improve resolution, we firstly interpolate the numerical values for kε c around its minimum and then find the minimum point of the more highly resolved interpolated data.

A.4 Vortex detection and averaging
We detect vortices by evaluating the phase winding of the transverse spin angle around plaquettes of our simulation grid. The average intervortex spacings in Fig. 3(b) and Fig. 4 is the vortex number for a single simulation trajectory at time t and angular brackets denote an ensemble average over the 30 simulation trajectories for the GPE results and the 10 simulation trajectories for the open system results. The results for L v in Fig. 1(b) are for the single trajectory displayed.

A.5 Open system simulations
To model open system evolution we couple our condensate to a reservoir of energy and particles with fixed temperature T and chemical potential µ. The dynamics is simulated using the simple growth stochastic Gross-Pitaevskii equations (SGPEs) [53,[57][58][59], Here L m [ψ m ] = − ħ h 2 ∇ 2 2M + qm 2 + g n n ψ m + g s n m F · f mm ψ m (15) is the conservative evolution operator from Eq. (3), µ is the chemical potential and γ is a dimensionless damping. The precise value chosen for γ will not affect equilibrium properties, but will affect the rate that equilibrium is approached. The term dW (r, t) is Gaussian distributed complex noise with delta correlations, The SGPEs Eq. (14) take the form of Langevin equations. The temperature (as a function of q) is chosen to be that obtained by equiparitioning the energy liberated by the quench amongst the 3N 2 numerical modes, as was done in the calculation of T eq for the microcanonical system. The energy liberated is the energy of the polar state evaluated at the final quadratic Zeemen energy q [20]. The temperature is then, The chemical potential (as a function of q) is chosen to be that of a zero temperature spin-1 condensate in the easy-plane phase, which is obtained by solving L m ψ m = µψ m . This gives [14], µ = g n n 0 + g s n 0 + q 2 .
These choices of T and µ give steady state energy and atom number within 1% of the conservative GPE results. We use γ = 10 −2 , which gives |dW | q 0 |ψ m |d t, and therefore reservoir scattering events occur within the time scale of spin interactions. The microscopically derived value for γ will be considerably smaller than this [53], resulting in the GPE dynamics overwhelming the reservoir interactions. We evolve Eq. (14) using an interaction picture fourth order Runge-Kutta integrator with periodic boundary conditions and kinetic energy evaluated to spectral accuracy. The noise is added in a single step following the Runge-Kutta integration of the (1 − iγ)(L m [ψ m ] − µψ m ) term [60]. Numerical parameters and initial condition sampling are the same as for the GPE simulations.