Multifaceted phase ordering kinetics of an antiferromagnetic spin-1 condensate

We study phase domain coarsening in the long time limit after a quench of magnetic field in a quasi one-dimensional spin-1 antiferromagnetic condensate. We observe that the growth of correlation length obeys scaling laws predicted by the two different models of phase ordering kinetics, namely the binary mixture and vector field. We derive regimes of clear realization for both of them. We demonstrate appearance of atypical scaling laws, which emerge in intermediate regions.

www.nature.com/scientificreports/ Halperin classification. All this occurs depending on the system's parameters. The scaling exponent z d = 4 has not yet been reported in the PhOK studies with quantum systems.
Using parameters for which we observe a change of scaling exponent passing from the short to long time scales, the characteristic length L(t) can be even a subject to multi-scaling behaviour 30,38,39 . The reason for such a multifaceted PhOK is that the spin-1 antiferromagnetic condensate effectively behaves as a binary mixture or vector fields model which is assigned to the H or B model, respectively. Therefore, the PhOK of the system can exhibit features characteristic for the H or B models independently, or both of them simultaneously. A study of the interplay between the models can be made. This is interesting because to date studies of the PhOK attributed system's behaviour mostly to a particular scaling law characterised by a particular model. A natural questions arises whether distinct models, and hence their physical mechanisms, are mutually compatible, or under what conditions they co-occur. Here, we characterize and classify the appearance of various scaling exponents. We calculate borders for the limit cases in which the B and H models can be realized in their forms. In the region around borders, both models compete, leading to various scaling laws in which scaling exponents smoothly change among the two limiting cases.

Model and methods
A spinor Bose-Einstein condensate of N sodium atoms is considered 40,41 . The system is represented by the vector � ψ = (ψ 1 , ψ 0 , ψ −1 ) T , whose components describe atoms in the corresponding Zeeman levels numerated by the magnetic numberm F = 0, ±1 . We assume a ring-shaped quasi-one-dimensional geometry with periodic boundary conditions [42][43][44] , where transverse degrees of freedom are confined in a strong potential with frequency ω ⊥ . The Hamiltonian of the system is where m is the atomic mass, n = n m F = ψ † m F ψ m F is the local atom density and F = (ψ † f x ψ, ψ † f y ψ, ψ † f z ψ) is the spin density with the spin-1 matrices f x,y,z . The spin-independent and spin-dependent interaction coefficients, c 0 and c 2 , are both positive for sodium atoms. Namely, they are c 0 = 2 ω ⊥ (2a 2 + a 0 )/3 and c 2 = 2 ω ⊥ (a 2 − a 0 )/3 , where a S is the s-wave scattering length for pairs of colliding atoms with total spin S 40 . The term q c 2 ρ is the quadratic Zeeman energy, where the dimensionless parameter q can be controlled using magnetic field or the microwave dressing 40 , and ρ is the mean density of the system. The Hamiltonian conserves the total atom number N = m F N m F and the magnetization M = N 1 − N −1 . The ground state of the system in the thermodynamic limit 45 is presented in Fig. 1.
In the present study we consider the non-negative magnetization and choose the 2C state as an initial state. Next, the quadratic Zeeman shift is set to a fixed value q > q c , and the evolution starts. The reason to consider quench from 2C towards 2C + ρ 0 , ρ + + ρ 0 phases is that the symmetry breaking occurs in this direction. We describe dynamics of the system on the mean-field level by solving the time-dependent Gross-Pitaevskii (GP) equations Figure 1. (i) The phase diagram of antiferromagnetic BECs 41,45 in the thermodynamic limit (when the spin healing length ξ s = / √ 2m c 2 ρ is much smaller than the linear system size L, i.e., ξ s ≪ L ) and positive magnetization. Three homogeneous phases (marked by colors) can be revealed: 2C phase (green) in which the components m F = ±1 coexist; ρ 0 (red) and ρ + (yellow) phases in which atoms occupy the m F = 0 and m F = +1 Zeeman state, respectively. The ground state of the system is (ii) the 2C phase when q < q 1 = (M/N) 2 /2 (solid white line), (iii) the phase separated into 2C and ρ 0 domains for q ∈ (q 1 , q 2 ),and (iv) the phase separated into ρ + and ρ 0 domains when q ≫ q 2 . Here q 2 is a ' contractual' crossover between the ρ 0 + 2C and ρ + + ρ 0 phases. In 45 the value of q 2 = 1/2 was derived in the thermodynamic limit omitting any excitation. In fact, q 2 depends on the system size and, e.g., thermal excitation as was pointed out in the recent experimental work 41 . The 2C phase remains a local energy minimum up to q c = 1 − 1 − (M/N) 2 (solid black line) 46 . The region between the white and black solid lines is where the two phases 2C and ρ 0 + 2C are stable. The vertical thick lines in (iii) and (iv) illustrate domain walls. The blue arrow in (i) shows an example of the considered quench for M/N = 0.7 and q = 0.5. www.nature.com/scientificreports/ where n m F = |ψ m F | 2 and δ m F ,0 is the Kronecker delta function, see e.g. 47 . To obtain the initial state for an arbitrary chosen value of M, all atoms are prepared in the polar ground state � ψ pgs = (0, ψ pgs (x), 0) T . This state is then subject to double rotations: (i) a spin-1 rotation e if y π/2 , which produces the intermediate state (1, 0, −1) T , and (ii) a rotation e −iσ y θ through angle θ that is performed on the m F = ±1 levels around the y-Pauli matrix σ y . The above procedure leads to � ψ M = ψ pgs √ 2 (sin θ + cos θ, 0, sin θ − cos θ) T , and the desired state for a given M is constructed when 2θ = arcsin M N . The state for arbitrary M can be also prepared experimentally by applying the two subsequent electromagnetic pulses 48 . In our calculations, stochastic white noise with variance = 1 2 particle momentum mode is added to all Zeeman components of the initial ψ M to seed the formation of symmetrybreaking phase domains. The calculations are made for an ensemble of N r = 300 realizations, the number of grid points N = 2 10 , 2 11 on the box size L large enough to avoid finite size effects, i.e., L = 7 × 10 2 , 10 3 , 5 × 10 3 , 10 4 µ m, which corresponds to ρ = 14.3, 10, 2, 1 µm −1 , respectively. We set the number of atoms to N = 10 4 .

Scientific Reports
As illustrated in Fig. 2 in the example for the m F = 0 component, the growth of phase domains and their coarsening is observed. The evolution of the density and phase can be divided into three stages: (i) creation of domains seed followed by spin domain formation around t = 10τ , (ii) early dynamics characterized by fast reduction of the number of domains, and (iii) further dynamics leading to domains merging at the longest time scale.
The initial exponential growth followed by the phase domains formation is well understood 46 so far relying on unstable modes. Here, in turn, PhOK in the long time limit driven by phase domains merging is described. The description relies on the evolution of the correlation length determined from the first-order correlation function g is a solution of the Gross-Pitaevskii equations (GPEs) for the m F = 0 Zeeman component and N = � dx|ψ 0 | 2 � is the normalization factor. The computed correlation length l h has such a property that g 2 ). Our initial condition, although trivialized, determines the following. No atoms in the m F = 0 state implies that its l 1/2 = 0 at t = 0 . Since atoms entirely occupied m F = ±1 states, their correlation length is of the order of the system size. The situation reverses dramatically, when the system is abruptly quenched. Therefore, the analysis of the superfluid order from the correlation length of the m F = 0 spin component seems suitable for the PhOK investigation in our system. The units chosen for the characteristic length and time in the spin-1 system are ξ s and τ = /(c 2 ρ) , respectively. Their values depends on the system density ρ . When the grow of correlation length is re-scaled by the two units, ξ s and τ , the scaling laws for different values of ρ overlap each other, see example in Fig. 5. The dimensionless quantities l 1/2 /ξ s and t/τ allows analysing the scale-invariant character of scaling laws.

Scaling hypothesis and various scaling laws: numerical results
The scaling hypothesis states that during the PhOK process a single length scale L(t) is expected which increases in time as L(t) ∼ t 1/z d , where 1/z d is the scaling exponent. The whole variation of the correlation function (1) that occurs during PhOK becomes independent of time, when it is scaled by L(t). We verify the hypothesis in our system and confirm that the correlation function g (1) N (x, t) is one parameter-dependent after the proper re-scaling, i.e., g N (x, t) = f (x/L(t)) . However, the value of the scaling exponent 1/z d strongly depends on the system parameters.
N (x, t) (upper panels) and f (x/L(t)) (bottom panels) with the properly chosen the L(t) scaling. Our numerical calculations show that the scaling exponent is 1/3 and 3/2 (or equivalently the dynamical exponent z d = 3 and z d = 3/2 ) for M = 0 and for a nonzero magnetization value when q is relatively large, respectively, see an example in Fig. 3c,d,g,h. The two scaling exponents have already been reported for spin-1 Bose-Einstein condensates with ferromagnetic interactions characterized by negative c 2 17,21 . These results can be explained on the basis of the hydrodynamic (binary mixture) model, as discussed by us in Appendix 5. In the case of macroscopic magnetization and a relatively low value of q the best matching to numerical data is observed when using scaling exponents predicted by the vector model with L(t) ∼ (t/ ln(t)) 1/4 , see an example in Fig. 3b,f. The temporal scaling L(t) ∼ t 1/2 can also be observed, see panels (a) and (e).
In the next two subsections, we discuss in details the particular scaling laws characteristic for the H and B models, and the extent of their occurrences. Our numerical findings are summarized in the last subsection in the form of the phase diagram for the spin-1 antiferromagnetic BECs.
In the limit of zero magnetization. Let us start the analysis with the zero magnetization case for which the scaling laws associated with the binary mixture or H model in the Hohenberg and Halperin classification of dynamic critical phenomena can be observed.
In the model the fluid flow contributes to the transport of the order parameter 1 , in general. This is why hydrodynamical processes like inertial or viscous growth along with diffusion mechanism are included. Each mechanism dominates at a different stage of the domains formation, and eventually one wins at the longest time scale. The H model predicts the following scaling laws: ∼ t 1/3 when the diffusive transport of the order parameter is dominant, ∼ t 2/3 when the inertia of fluids are important, and ∼ t 1 if the viscous process wins 1,49 .
Our numerical GPEs results demonstrate that the diffusive ∼ t 1/3 and inertial hydrodynamic ∼ t 2/3 scaling laws are revealed in the long time limit, while the scaling law ∼ t 1 from a pure viscous effect is absent. The exact derivation of scaling exponents from GPEs is difficult, but the estimated analysis can be performed-a similar one was already done on the hydrodynamic equations for the binary mixture 1 . We discuss this in Appendix 5. In general, the estimation of scaling laws for the H model requires introduction of a non-zero surface tension. In one-dimensional system described by GPEs, it equals zero. The scaling analysis we performed show that the interaction coefficient c 0 plays the role of surface tension and in a result, he scaling analysis of the hydrodynamic form of GPEs gives the same scaling laws as the typical H model, see Appendix 5. Consequently, the appearance of the two scaling exponents 1/3 and 2/3 in the systems described by the GPEs can be justified.
The transition from the diffusive to the intertial hydrodynamic scaling laws is clearly visible in Fig. 4b. We observe that the average domain wall width is comparable with the width of the phase domain itself at the transition point where the scaling exponent changes. If the size of the phase domain is larger than the size of the average domain wall, then the diffusion transport defines the physics of the system and the scaling exponent. We use that reasoning to estimate the transition point between z −1 d = 1/3 and z −1 d = 2/3 scaling laws. To proceed, let us assume that the fractional size of phase domain over the entire system is given by the fractional volume occupied by the ground state phase ρ 0 , composed of spin domains, i.e., x 0 (q) = 1 − √ q 1 √ q . This formula can be established from the analysis of equilibrium conditions for the coexistence of 2C and ρ 0 phases 46 . The width of The relation x 0 (q th ) ≈ ξ 2C /ξ s or equivalently 1 − √ q 1 √ q th ≈ √ 1 − q th /q th , gives the desired condition for the transition point between the two different scaling laws. This point is noted here as q th . Whenever q > q th , the diffusive transport governs domains coarsening. The resulting estimate for q th is shown in Fig. 4b by the vertical solid line, and in Fig. 7 by the dashed line. At this juncture, it is useful to connect the appearance of domains of phase (PhOK) and also domains in the real space, to the occupations of atoms in particular Zeeman states. The system we consider is three component in general. When M decreases at relatively low q the m F = 0 component becomes macroscopically occupied at longer time scales, while the remaining two components m F = ±1 are occupied marginally, both to almost the same extent. For the purposes of the dynamics one can expect that the properties of atoms remaining in the m F = ±1 components become identical. The system starts then to behave like a binary mixture composed of two species of atoms: these in the m F = 0 state and those in m F = ±1 treated together as the second species. Similarly, in the large q limit occupation of the m F = −1 component becomes marginal while the remaining two m F = 0, 1 are of the only importance. In both regimes of parameters, the binary mixture description of PhOK could become suitable for our system, and therefore, the resulting scaling exponents could appear.

Macroscopic magnetization.
In the case of macroscopic magnetization M, we can talk about the regime of parameters for which occupations of all three Zeeman components are significant after the quench in the long time limit. Also then the variation of the scaled correlation length l 1/2 /ξ s versus scaled time t/τ is expected, as shown in Fig. 5 for various system densities ρ . This time, an increase of the correlation length does not exhibit the scaling attached to the binary mixture models. Instead, we observe that l 1/2 reveals growth typical for vector field models 1 . Those models predict L(t) ∼ t 1/2 (model A) and L(t) ∼ (t/lnt) 1/4 (model B) for non conserved and conserved order parameters, respectively.
Better understanding of the evolution of the order parameter in the case of considered system can be provided with definition of n 0 = N 0 /N-the fractional number of atoms in the m F = 0 Zeeman component, where Fig. 5b. Notice that for early times n 0 is susceptible to large fluctuations.  www.nature.com/scientificreports/ In the range ∼ 10τ . . . 40τ , the order parameter is found to be evidently non conserved. The correlation length scaling exponent is 1/2 there, see Fig. 5c. However, for the longest time scale which we are interested in, the change of n 0 is less significant although its variance is far from being equal to zero. Despite that fact, the value of the scaling exponent changes to 1/4 with logarithmic correction. The two scaling laws can also be observed by properly re-scaling the correlation function g N (x, t) as demonstrated in Fig. 3. Note, however, that the universal character can have only the scaling revealed in the long time limit, that is z −1 d = 1/4 . The temporal 1/2 scaling is mentioned by us only to demonstrate the multiscaling behaviour of PhOK characteristic for our system in early times. The similar observation and conclusion was also made for ferromagnetic spin-1 Bose-Einstein condensates 15 .
The character of the system smoothly transforms to the binary mixture when the value of q increases. The occupation of the m F = −1 Zeeman component diminishes to zero and the remaining two, m F = 0, 1 are of the only importance. The change of the scaling exponent to 1/3 typical for the H model is observed by us independently of M.
Classification of scaling laws. To summarize the results at this stage: we show the variation of scaling laws in time for q = 0.5 , 0.75, 1.0, 1.25 in Fig. 6 depending on the magnetization value. Our observations confirm the presence of mechanisms which are competitive with each other and they lead finally to different behaviour in the long times limit. Therefore, several values of the growth exponents can emerge during the later evolution, namely z d → 3/2 , z d → 3 or even to z d → 4 with logarithmic correction. The change of scaling exponent while varying magnetization is typical for our system and appears in a wide range of the value of q.
Careful analysis of the numerical results indicates an important role of occupations of particular Zeeman states. We conclude that occupations of atoms in particular Zeeman components give an intuitive picture of which of the scaling law will be revealed in the long dynamics of our system. The vector (B) model fits the best when occupations of all the three Zeeman components are significant. It takes the place for macroscopic magnetization for values of q a bit larger than the critical one, q c . In turn, a high concentration of atoms can stimulate stronger interactions, the liquid character and thus the hydrodynamic description would be more accurate. It is when magnetization is small and in the large q limit, q ≫ q c . In this regime, the hydrodynamic (H) model seems to match the best.
The border between both models can be deduced by matching when the number of atoms in the m F = −1 component vanishes, i.e., N −1 (M, q)/N → 0 . This is illustrated in Fig. 7. The vector field model is realized below the gray solid border line, while the binary mixture model applies above it. The change between the two regions is smooth, and so z d does not perfectly reflect the particular model all around the border line.

Discussion and conclusion
In this paper, we have explored the superfluid phase-ordering dynamics of an antiferromagnetic spin-1 condensate quenched from the antiferromagnetic state to a state where domains of atoms with different spin projections separate. We have found that the growth of domains is scale-invariant with various dynamic critical exponent z d that are typical for the B and H models. We classified the various scaling exponents due to the system parameters. www.nature.com/scientificreports/ When the occupation of the m F = −1 Zeeman state is significant in the final state, we observe that the scaling exponents typical for vectors fields model appear ( z d = 4 with logarithmic correction). It is when the system is quenched from the initial 2C phase to the final ρ 0 + 2C phase excluding the region where M → 0 and q < q 2 and to the final ρ + + ρ − phase for macroscopic magnetization and q < 1 . At short times we observe also that temporally the correlation length l 1/2 exhibits the scaling L(t) ∼ t 1/2 . These all properties are characteristic for multiscaling of the PhOK. On the other hand, when the occupation of the m F = −1 component is marginal in the final state the scaling laws typical for the H model are realized, namely z d = 3/2 when q and magnetization are small, and z d = 3 in the large q limit for any M. It is when the system is quenched from the initial 2C phase to the ρ 0 + ρ + phase including the region of the ρ 0 + 2C phase where M → 0 and q < q 2 . While the appearance of scaling laws typical for the H model can be justified from the hydrodynamic form of GPEs, the occurrence of scaling typical for B model is not evident. The latter provides an interesting direction for future work.
To conclude: the antiferromagnetic spin-1 condensate captures the universal two-model feature of the PhOK in the long time limit. The system parameters (q, M/N) set the physics, and determine the dynamical scaling exponent corresponding to the model H or B. Switching between the models by changing the initial state or the q parameter is allowed.

A Scaling analysis from spin-1 hydrodynamic equations
In the following, we present the hydrodynamic forms of the GPEs. The wave functions in Eq. (3) describing atoms in particular components are replaced by ψ m F = √ n m F e i θ m F -according to Madelung transformation. After some algebra, we obtain corresponding equations for the densities n m F , namely and the corresponding Navier-Stokes equations for velocities v m F = m ∇θ m F w h e r e θ = θ 1 + θ −1 − 2θ 0 , a n d U 1 = c 0 n + qc 2 ρ + c 2 (n 1 − n −1 + n 0 ) , U 0 = c 0 n + c 2 (n 1 + n −1 ) , U −1 = c 0 n + q c 2 ρ + c 2 (n −1 − n 1 + n 0 ).
We will perform here the scaling analysis corresponding to the one presented in 1 for the hydrodynamic (B) model. The chemical potential can be identified as μ = U 0 , see details in 46 , and therefore the derivative of chemical potential as φ∇μ → ∇U 0 . Taking this as an assumption one can perform the scaling analysis for m F = 0 as in below. 1 Scaling by diffusion mechanism: In the diffusive regime it is assumed that the rate of change of the domain size dL dt is associated with the chemical potential gradient ∇U 0 . The scaling of U 0 is linked to the scaling of the local density U 0 ∼ c 0 n ∼ c 0 N/L , therefore one has: Note, the role of the surface tension σ typically introduced in analysis of binary mixture 1 is played here by the interaction coefficient c 0 .
2 Scaling by inertial growth: If the inertial terms become important, the scaling law is given by the relation (v 0 ∇)v 0 ∼ ∇U 0 . Then, one recovers the same scaling as for the H model, namely (4) n m F + v m F ∇n m F = − n m F ∇v m F + 2 c 2 (|m F | − 2) n 0 √ n 1 n −1 sin(θ)