Nematic ordering dynamics of an anti-ferromagnetic spin-1 condensate

TL;DR

This paper investigates the formation and evolution of nematic order in a quasi-2D anti-ferromagnetic spin-1 condensate after a phase quench, highlighting the role of half-quantum vortices in the coarsening process.

Contribution

It introduces a new analysis of spin-nematic and superfluid order dynamics, emphasizing the role of half-quantum vortices in the coarsening process of the condensate.

Findings

Dynamical scaling observed in late-time order growth

Both types of order exhibit diffusive growth laws

Half-quantum vortices drive the annihilation process

Abstract

We consider the formation of order in a quasi-two-dimensional (quasi-2D) anti-ferromagnetic spin-1 condensate quenched from an easy-axis (EA) to an easy-plane (EP) nematic phase. We define the relevant order parameter to quantify the spin-nematic degrees of freedom and study the evolution of the spin-nematic and superfluid order during the coarsening dynamics using numerical simulations. We observe dynamical scaling in the late time dynamics with both types of order extending across the system with a diffusive growth law. We identify half-quantum vortices (HQVs) as the relevant topological defects of the ordering dynamics, and demonstrate that the growth of both types of order is determined by the mutual annihilation of these vortices.

Figures

Equations

$$\mathbf{\psi }\equiv ({\psi }_1,{\psi }_0,{\psi }_{-1}{)}^T,$$

$${\mathcal{H}}_{\mathrm{int}}=\frac{g_n}{2}{n}^2+\frac{g_s}{2}{\left|\mathbf{F}\right|}^2.$$

$${\mathcal{H}}_{\mathrm{QZ}}=q{\mathbf{\psi }}^{†}{f}_z^2\mathbf{\psi },$$

$$\vec{\psi }=e^{i\theta }(\vec{u}+i\vec{v}),$$

$$\mathbf{F}=-i{\vec{\psi }}^{\ast }\times \vec{\psi }=2\vec{u}\times \vec{v},$$

$${\mathcal{N}}_{ab}$$

$$=n{\delta }_{ab}-\frac{1}{2}({\vec{\psi }}^{\ast }\otimes \vec{\psi }+\vec{\psi }\otimes {\vec{\psi }}^{\ast }),$$

$$=n{\delta }_{ab}-(\vec{u}\otimes \vec{u}+\vec{v}\otimes \vec{v}).$$

$$\alpha =\vec{\psi }\cdot \vec{\psi }={\psi }_0^2-2{\psi }_1{\psi }_{-1},$$

$${\vec{\psi }}_{\mathrm{P}}=e^{i\theta }\vec{u},$$

$$\theta \to \theta +\pi ,\mathrm{and}\vec{u}\to -\vec{u}.$$

$$|\mathbf{F}{|}^2+{\mathcal{A}}^2=n^2,$$

$$Q$$

$$\displaystyle=\left(\begin{array}[]{cc}Q_{xx}&Q_{xy}\\ Q_{xy}&-Q_{xx}\end{array}\right),$$

$${\mathcal{A}}_{\perp }=|{\alpha }_{\perp }|,$$

$${\alpha }_{\perp }\equiv -2{\psi }_1{\psi }_{-1}.$$

$$i\hslash \frac{\partial \mathbf{\psi }}{\partial t}=\left(-\frac{{\hslash }^2{\nabla }^2}{2M}+q{f}_z^2+g_{n}n+g_s\mathbf{F}\cdot \mathbf{f}\right)\mathbf{\psi }.$$

$$\displaystyle\bm{\psi}(\mathbf{x},t=0)=\sqrt{n_{c}}\left(\begin{array}[]{c}0\\ 1\\ 0\end{array}\right)+\bm{\delta}(\mathbf{x}),$$

$$⟨O^2(t)⟩=\frac{1}{l^2}\int d^2\mathbf{x}{O}^2(\mathbf{x},t).$$

$$R_{m{m}^{\prime }}(t)=⟨\delta n_m\delta n_{m^{\prime }}⟩,$$

$$⟨{\mathcal{A}}_{\perp }^2{⟩}_{\mathrm{uc}}\approx 4⟨n_1⟩⟨n_{-1}⟩.$$

$$R_{\Theta }(t)\equiv ⟨{\psi }_{-1}{\psi }_1{\psi }_0^{\ast }{\psi }_0^{\ast }⟩.$$

$$⟨|F_{\perp }{|}^2⟩$$

$$⟨{\mathcal{A}}^2⟩$$

$$G_{\varphi }(\mathbf{r},t)$$

$$G_{\theta }(\mathbf{r},t)$$

$$\displaystyle\bm{\psi}_{\mathrm{EP}}=\sqrt{\frac{n_{c}}{2}}e^{i\theta}\left(\begin{array}[]{c}-e^{-i\phi}\\ 0\\ e^{i\phi}\end{array}\right),$$

$$G_{\varphi }^{\mathrm{EP}}(\mathbf{r})$$

$$G_{\theta }^{\mathrm{EP}}(\mathbf{r})$$

$$G_{\varphi }(r,t)=\int d{\Omega }_r\int \frac{d^2{\mathbf{x}}^{\prime }}{l^2}\frac{2}{n_c^2}⟨\mathrm{Tr}\{Q({\mathbf{x}}^{\prime })Q({\mathbf{x}}^{\prime }+\mathbf{r})\}{⟩}_t,$$

$$H_{\varphi }(r)=G_{\varphi }(r/L_{\varphi }(t),t).$$

$$H_{\theta }(r)=G_{\theta }(r/L_{\theta }(t),t).$$

$$S_{\varphi }(\mathbf{k},t)=\int d^2\mathbf{r}{G}_{\varphi }(\mathbf{r},t)e^{i\mathbf{k}\cdot \mathbf{r}}.$$

Full text

I Introduction

Spin-1 condensates Stenger et al. (1998); Ho (1998); Ohmi and Machida (1998) with anti-ferromagnetic interactions prefer to order into spin-nematic phases Kawaguchi and Ueda (2012). Such phases have a vanishing average spin-density, and are instead characterized by the nematic tensor $\mathcal{N}_{ab}=\tfrac{1}{2}\langle f_{a}f_{b}+f_{b}f_{a}\rangle$, where $f_{a\in\{x,y,z\}}$ are the spin matrices. The ground states of this system have an axially symmetric nematic tensor (uniaxial nematic) with a preferred axis (but not direction) characterized by a director $\vec{u}$ in spin-space (i.e. $\vec{u}$ and $-\vec{u}$ are equivalent). Recently experimental evidence was presented for spin-nematic order in an anti-ferromagnetic condensate Zibold et al. (2016).

The concept of nematic order is typically discussed in the context of liquid crystals, where the order is associated with the orientation of long molecules. Indeed, many beautiful studies of phase transition dynamics and coarsening have been performed in liquid crystal systems (e.g. see Orihara and Ishibashi (1986); Chuang et al. (1991a, b); Yurke et al. (1992); Nagaya et al. (1992); Blundell and Bray (1992); Chuang et al. (1993)). A sudden change in conditions (e.g. temperature or pressure) is used to quench this system from an isotropic phase (unoriented molecules) into the nematic phase, and the formation of order and defect dynamics can be observed optically.

In this paper we develop a theory for the ordering dynamics (coarsening) of an anti-ferromagnetic spin-1 condensate. There has been considerable theoretical work on the coarsening dynamics of ferromagnetic spin-1 condensates Mukerjee et al. (2007); Kudo and Kawaguchi (2013, 2015); Williamson and Blakie (2016a), however this area is largely unexplored in the anti-ferromagnetic system. Our interest is in the symmetry breaking phase transition from an EA phase (with $\vec{u}$ along the direction set by the external magnetic field) at positive quadratic Zeeman energy $q$, to an EP phase ($\vec{u}$ transverse to the external field) at negative $q$ (see Fig. 1(d) and Refs. Bookjans et al. (2011); Vinit et al. (2013); Phuc et al. (2013); Jiang et al. (2014); Vinit and Raman (2017)). We consider a quench between these phases implemented by a sudden change in $q$, e.g. using microwave dressing (see Bookjans et al. (2011); Vinit et al. (2013); Jiang et al. (2014)). Upon entering this new phase, the system breaks the continuous axial symmetry of the initial state by developing transverse spin-nematic domains. Here our interest lies in characterizing the dynamics of the phase transition, with a particular emphasize on the late-time coarsening dynamics. That is, to understand the universal aspects of how small domains created after the quench anneal together to bring the system towards an ordered equilibrium state. To undertake this study we first discuss how the nematic order is characterized in a spinor condensate, and develop an appropriate order parameter for the EP phase. Using numerical simulations we study how the EP order forms in the system. We demonstrate that the late-time coarsening behaviour exhibits dynamical scaling with a diffusive domain growth law of $L(t)\sim[t/\ln(t)]^{1/2}$, where $L$ is the size of the ordered domains and $t$ is the time after the quench. We separately consider the superfluid order and show that it grows with an identical law to the spin-nematic order, in contrast to recent results for the ferromagnetic spin-1 system Bourges and Blakie (2017). The order parameter growth is determined by the dynamics of HQVs in the system, and we verify that the number of these vortices scales as $L(t)^{-2}$, i.e. that coarsening proceeds by vortex anti-vortex pairs mutually annihilating. Recent experiments have demonstrated that it is possible to measure HQVs in anti-ferromagnetic spin-1 condensates Seo et al. (2015, 2016) due to their ferromagnetic cores Lovegrove et al. (2012). Thus, measuring the HQV distribution as a function of time after the quench could be a practical method for experiments to quantify the coarsening of this system. Alternatively, it may be possible to directly image Carusotto and Mueller (2004); Higbie et al. (2005) or probe Baillie and Blakie (2016) nematic properties of the condensate.

We note that the symmetries and defects of the EP phase are similar to those of a (two-component) binary condensate in the miscible regime. Indeed, work by Karl et al. Karl et al. (2013) discussed the role of the equivalent vortices in the ordering dynamics of a two-component system, although that work focused on understanding the emergence of power-law behavior in various momentum correlation functions, and relating these to turbulence cascades.

The outline of the paper is as follows. In Sec. II we introduce the basic formalism for spin-1 condensates and consider how to quantify spin-nematic order. We discuss the EA to EP quench and introduce the relevant order parameter for this phase transition. In Sec. III we start by introducing the quasi-2D system, the equation of motion and simulation technique we use to study the quench dynamics. We present results for the evolution of various local densities and correlation functions that illustrate the early time dynamics of the quench, and show the emergence of EP order. We then focus on the late time dynamics of the system and characterise the phase ordering dynamics. To do this we introduce correlation functions for the spin-nematic and superfluid order. We evaluate these using an ensemble of large-scale simulations and demonstrate correlation function collapse (dynamic scaling) and extract the relevant growth laws. Finally, we examine the role of HQVs and show that the average distance between vortices characterises the growth of order. Then we conclude in Sec. IV.

II Formalism

II.1 Spin-1 Anti-ferromagnetic condensate

A spin-1 condensate is described by the spinor field

[TABLE]

where the three components describe the condensate amplitude in the spin levels $m=1,0,-1$, respectively. The short-ranged interactions between atoms are described by the rotationally invariant Hamiltonian density

[TABLE]

The first term describes the density dependent interactions, with coupling constant $g_{n}$, where $n\equiv\bm{\psi}^{\dagger}\bm{\psi}$ is the total density. The second term describes the spin-dependent interactions $g_{s}|\bm{F}|^{2}$, with coupling constant $g_{s}$, where $\bm{F}\equiv\bm{\psi}^{\dagger}\bm{f}\bm{\psi}$ is the spin density and $\bm{f}\equiv(f_{x},f_{y},f_{z})$ are the spin-1 matrices. For the case $g_{s}>0$, known as anti-ferromagnetic interactions, the condensate prefers to minimise the spin-density to reduce the interaction energy. In addition to interactions, the (uniform) quadratic Zeeman shift

[TABLE]

also plays a role in determining the preferred spin-ordering of the condensate. The quadratic Zeeman energy $q$ can be controlled using the magnetic bias field, it can also be varied by using microwave dressing (e.g. see Gerbier et al. (2006); Leslie et al. (2009)).

II.2 Nematic order

To quantify the spin-order it is useful to introduce the Cartesian representation of the spinor field $\vec{\psi}\equiv(\psi_{x},\psi_{y},\psi_{z})$, where $\psi_{x}=(\psi_{-1}-\psi_{1})/\sqrt{2}$, $\psi_{y}=-i(\psi_{1}+\psi_{-1})/\sqrt{2}$, and $\psi_{z}=\psi_{0}$. We will give results in both the cartesian $\vec{\psi}$ and spherical [$\bm{\psi}$, see Eq. (1)] bases as needed.

A general spinor can be decomposed in the form

[TABLE]

where $\theta$ is the global phase and $\{\vec{u},\vec{v}\}$ are mutually orthogonal real vectors satisfying $|\vec{u}|^{2}+|\vec{v}|^{2}=n$, and $|\vec{u}|\geq|\vec{v}|$ (also see Ohmi and Machida (1998); Mueller (2004); Yukawa and Ueda (2012); Zibold et al. (2016)). For a spin-1 spinor, the local spin information described by the spin density vector is

[TABLE]

and the symmetric nematic (or quadrupolar) tensor density is

[TABLE]

The nematic tensor describes the anisotropy of the spin fluctuations, and in general has the symmetries of an ellipsoid. This is revealed by diagonalizing $\mathcal{N}$, giving $\{\vec{u},\vec{v},\mathbf{F}\}$ as the eigenvectors with respective eigenvalues $\lambda_{u}=\frac{1}{2}(n-\mathcal{A})$, $\lambda_{v}=\frac{1}{2}(n+\mathcal{A})$ and $\lambda_{F}=n$. Here $\mathcal{A}=2|\vec{u}|^{2}-n\geq 0$ is the alignment parameter Zibold et al. (2016), which characterizes the relative fluctuations of magnetization along the directions orthogonal to $\mathbf{F}$. The alignment is related to the spin-singlet amplitude111Note our definition differs by a constant factor from Kawaguchi and Ueda (2012).

[TABLE]

as $\mathcal{A}=|\alpha|$. It is conventional to take the eigenvector associated with the smallest eigenvalue of $\mathcal{N}$ as the nematic director, i.e. the vector $\vec{u}$. We can use the eigenvectors and eigenvalues to represent the nematic tensor density as an ellipsoid [see Fig. 1(a)]. We also note that $\lambda_{u}=|\vec{v}|^{2}$ and $\lambda_{v}=|\vec{u}|^{2}$, so that the extent of the ellipsoid along the $\vec{u}$ direction is the squared length of $\vec{v}$, and the extent of the ellipsoid along the $\vec{v}$ direction is the squared length of $\vec{u}$.

Two limiting states are of interest. First, the fully magnetized ferromagnetic state with $|\mathbf{F}|=n$, where $|\vec{u}|=|\vec{v}|=\sqrt{n/2}$, and $\mathcal{A}=0$. Second, and of primary concern in our work, is the fully polar (or spin-nematic) state which has the form

[TABLE]

with $|\vec{u}|=\sqrt{n}$, $\mathcal{A}=n$ and $\mathbf{F}=0$ [see Figs. 1(b) and (c)]. The spin properties of this state are completely characterized by the director $\vec{u}$, and the state is invariant under the transformation

[TABLE]

For general spin-1 states the relation

[TABLE]

holds, so that $\mathcal{A}$ can be used to characterize how close a state is to the limiting cases of ferromagnetic ($\mathcal{A}=0$) or polar ($\mathcal{A}=n$) order.

II.3 Order parameter for the EA to EP phase transition

Here we are concerned with an anti-ferromagnetic condensate in which a quench is performed by a sudden change in the quadratic Zeeman energy from a positive value to a negative value222The $z$-magnetization $M_{z}\equiv\int d^{2}\mathbf{x}\,F_{z}$ of the system is conserved, and here we focus on the case $M_{z}=0$ where the transition occurs at $q=0$. crossing a quantum phase transition between two different ground states [see Fig. 1(d)]. For both cases the ground state is fully polar $\vec{\psi}=e^{i\theta}\vec{u}$. For $q>0$ the director ($\vec{u}$) is along the $z$ axis [EA phase, see Fig. 1(c)]. For $q<0$ the director lies in the $xy$-plane [EP phase, see Fig. 1(b)]. Thus the EP phase breaks the axial symmetry (invariance to spin rotations about $z$) of the Hamiltonian. This type of quench in an anti-ferromagnetic spinor condensate of 23Na atoms has been performed in a number of experiments Bookjans et al. (2011); Vinit et al. (2013); Jiang et al. (2014); Vinit and Raman (2017); Kang et al. (2017), however the EP nematic order was not directly probed in these studies (c.f. Zibold et al. (2016)). We also note that other phase transitions can be considered in this system, e.g. Witkowska et al. Witkowska et al. (2014) considered a $q$ quench for an anti-ferromagnetic condensate with a non-zero (conserved) $z$-magnetization, where a transition to a phase separated state occurs.

We would like to obtain an order parameter that can distinguish between these two states, notably the order parameter should be zero in the EA phase and non-zero in the EP phase. To do this we note that in the EA phase the nematic tensor is isotropic in the $xy$-plane [see Fig. 1(c)], while in the EP phase the nematic tensor is anisotropic in the $xy$-plane [see Fig. 1(b)]. To quantify the EP nematic order, and taking motivation from nematic liquid crystals de Gennes and Prost (1995), we use a traceless symmetric tensor to quantify order in this system. Particular to the EA to EP phase transition we use the planar tensor:

[TABLE]

where $\mathcal{N}_{2\times 2}$ is the $xy$-submatrix of $\mathcal{N}$, and $I_{2}$ is the identity matrix. Evaluating this expression we find that $Q_{xx}=\mathrm{Re}\{\psi_{1}^{*}\psi_{-1}\}$ and $Q_{xy}=\mathrm{Im}\{\psi_{1}^{*}\psi_{-1}\}$, i.e. it depends on the relative phase coherence between the $\psi_{1}$ and $\psi_{-1}$ components of the system. While $Q$ is traceless by construction, $\mathrm{Tr}(Q^{2})=0$ only when the spin fluctuations are isotropic in the $xy$-plane. The EP phase is thus revealed by $\mathrm{Tr}(Q^{2})$ becoming non-zero, thus demonstrating how $Q$ serves as an order parameter. We can write the eigenvalues of $Q$ as $\{-\frac{1}{2}\mathcal{A}_{\perp},\frac{1}{2}\mathcal{A}_{\perp}\}$, where we have defined a “transverse alignment” parameter333$\mathcal{A}_{\perp}$ is sensitive to anisotropy of ${Q}$, but does not completely distinguish between polar and ferromagnetic states as does $\mathcal{A}$. E.g., the fully ferromagnetic state with $\mathbf{F}=n\hat{\mathbf{x}}$ has $\mathcal{A}=0$, but $\mathcal{A}_{\perp}=\frac{1}{2}n$ (c.f. the pure EP polar state with $\vec{u}=\sqrt{n}\hat{\mathbf{x}}$ for which $\mathcal{A}=\mathcal{A}_{\perp}=n$). As is apparent from Fig. 1(a) a state with $\mathbf{F}\neq 0$ also has an anisotropic nematic tensor, just to a lesser extent than a polar state.

[TABLE]

and have introduced [c.f. Eq. (9)]

[TABLE]

Using this result gives $\mathrm{Tr}(Q^{2})=\frac{1}{2}\mathcal{A}_{\perp}^{2}$. In Appendix A we present an alternative formulation of the planar tensor $Q$ and order parameter results.

III Results

III.1 Quasi-two-dimensional quench

In order to explore the quench dynamics we focus on a quasi-2D system. In this regime the extent of the condensate in one direction (which we take to be $z$) is less than the spin healing length, so spin motion is effectively frozen out in this direction. This regime has been realized in experiments by applying a tight optical trap in this direction (e.g. see Sadler et al. (2006); Seo et al. (2015)). Additionally, we neglect any transverse confinement and take the condensate to be homogeneous in the plane. The dynamics of this system is described by the spin-1 Gross-Pitaevskii equation (GPE)

[TABLE]

Note we have neglected the linear Zeeman shift which can be removed from the equation of motion by transforming to a rotating frame.

To numerically solve this equation we represent each component of the spinor field $\bm{\psi}$ on a 2D square region of dimensions $l\times l$ covered by an $N\times N$ grid of equally spaced points. Taking periodic boundary conditions for the solution we evaluate spatial derivatives in the kinetic energy term of Eq. (19) with spectral accuracy using fast Fourier transforms. To evolve the GPE in time we use the second order symplectic method presented in Ref. Symes et al. (2016).

The initial condition for the simulations is a uniform EA ground state (in the spherical basis)

[TABLE]

where ${n_{c}}$ is the condensate (areal) density and $\bm{\delta}$ is a small noise field added to seed the growth of unstable modes following the quench. The late-time results are insensitive to the form of white spatial noise we add to the initial state as long as the noise is weak ($|\bm{\delta}|^{2}\ll n_{c}$). We choose to add noise according to the truncated Wigner prescription Blakie et al. (2008), which is consistent with the quantum vacuum noise on the initial state (see Williamson and Blakie (2016b) for details). We introduce the characteristic spin energy $q_{0}\equiv 2g_{s}n_{c}$, and associated spin healing length $\xi_{s}=\hbar/\sqrt{Mq_{0}}$ and spin time $t_{s}=\hbar/q_{0}$ as convenient units.

III.2 Early-time dynamics: development of local order

Immediately following the quench the initial EA state is unstable and begins to evolve towards the new phase. Aspects of these early time dynamics, and the emergence of local EP order can be revealed by studying the behaviour of the spin and alignment densities. Since some of these densities (e.g. $F_{z}$) can be locally negative, we quantify the development of a particular density of interest $O$ by spatially averaging $O^{2}$, i.e. we evaluate

[TABLE]

We present results for a variety of densities of interest in Fig. 2(a). These results show that immediately following the quench the EA state becomes dynamically unstable to magnon excitations which grow exponentially and cause the system to develop transverse magnetization [i.e. $\mathbf{F}_{\perp}=(F_{x},F_{y})$]. The precise nature of the instability and the wavevectors of the unstable modes depends upon the value of $q$, and aspects of this have already been explored in experiments Bookjans et al. (2011); Vinit et al. (2013); Vinit and Raman (2017); Kang et al. (2017). The axial magnetization ($F_{z}$) similarly experiences exponential growth. The general behavior of spin density growth we observe is similar for quenched condensates with ferromagnetic interactions (e.g. see Saito et al. (2007); Lamacraft (2007); Leslie et al. (2009); Barnett et al. (2011)). Noting that the average $z$-magnetization of the initial state is zero (and conserved), the quantity $\langle F_{z}^{2}(t)\rangle$ corresponds to the fluctuations in magnetization studied in recent experiments Kang et al. (2017).

More direct insight into the change in nematic order is provided by the alignment densities $\{\mathcal{A},\mathcal{A}_{\perp}\}$ discussed in Sec. II.2. The initial EA state is fully aligned (i.e. $\mathcal{A}=n_{c}$), but this dips down in the early dynamics as the magnetization develops [as required by the relation (12)]. As the alignment is restored for $t\gtrsim 20\,t_{s}$ it is of a different character, consistent with EP order emerging. We see this by evaluating the transverse alignment $\mathcal{A}_{\perp}$ order which is initially negligible, but then grows and is seen to saturate towards the value of $\mathcal{A}$.

Various in situ measurements of correlations between components of the density have been performed in spinor condensate experiments (e.g. see Vinit et al. (2013); Guzman et al. (2011); Seo et al. (2016); Kang et al. (2017)). Most relevant to our system are the measurements of Vinit et al. Vinit et al. (2013) of the time evolution of the local pair correlation function following the EA to EP quench of a quasi-one-dimensional anti-ferromagnetic condensate. The correlation functions measured were444Here and for the remainder of this subsection all expectations will be taken to be spatially averaged as in Eq. (24).

[TABLE]

where $\delta n_{m}(\mathbf{x},t)=n_{m}(\mathbf{x},t)-\langle n_{m}\rangle$ is the $m$-component density fluctuation operator, with $n_{m}=|\psi_{m}|^{2}$ and $\langle n_{m}\rangle$ being the mean density of this component. We have evaluated the same correlation functions measured in experiments (c.f. Fig. 3 of Ref. Vinit et al. (2013)) and present the results in Fig. 2(b). We find similar qualitative behavior to their results, however note that their measurements were for a shallow quench (to $q\approx-0.02q_{0}$) and with appreciable thermal effects. These same types of local density measurements could be used to evaluate the alignment densities. Indeed, noting that $\langle\mathcal{A}_{\perp}^{2}\rangle=4\langle n_{1}n_{-1}\rangle$ [see Eqs. (17) and (18)] , taking $n_{1}$ and $n_{-1}$ as uncorrelated, we can make the estimate

[TABLE]

For the uniform system $\langle n_{m}\rangle=N_{m}/l^{2}$, and thus $\langle\mathcal{A}_{\perp}^{2}\rangle_{\mathrm{uc}}$ is determined by the component populations $N_{m}=\int d^{2}\mathbf{x}\,n_{m}$, which are readily measured in experiments. As can be seen from Fig. 2(a) the uncorrelated approximation tends to overestimate the EP order ($\langle\mathcal{A}_{\perp}^{2}\rangle$) once it develops ($t\gtrsim 20\,t_{s}$). Noting that $\langle n_{1}n_{-1}\rangle=\langle n_{1}\rangle\langle n_{-1}\rangle+R_{1,-1}$, this overestimate of Eq. (26) is due to the negative value $R_{1,-1}$ takes for $t\gtrsim 20\,t_{s}$ [Fig. 2(b)]. Evidence for $R_{1,-1}$ becoming negative was also found in experiments at late times Vinit et al. (2013).

Finally we examine the system evolution to quantify the local “phase locking” of the $m=\pm 1$ components relative to the $m=0$ component. This was recently observed in experiments by applying a spin rotation to the system and measuring the resulting magnetic fluctuations Zibold et al. (2016). In our simulations we can directly access this from the local (spatially averaged) correlation function

[TABLE]

Taking $\psi_{m}=\sqrt{n_{m}}e^{i\theta_{m}}$, we see that $R_{\Theta}\sim e^{i(\theta_{1}+\theta_{-1}-2\theta_{0})}$, which conventionally defines the relative phase $\Theta\equiv\theta_{1}+\theta_{-1}-2\theta_{0}$. To understand the physical relevance of this correlation function, we note that the transverse spin density squared and the alignment density squared are

[TABLE]

respectively. Thus varying the real part of $R_{\Theta}$ the system can enhance or reduce the spin density, while having the opposite effect on the alignment [also see Eq. (12)]. Anti-ferromagnetic systems prefer $\Theta=\pi$ to reduce the spin-density. The behaviour of $R_{\Theta}$ is shown in Fig. 2(c), noting that we have normalized $R_{\Theta}$ by the average densities of each component [using $\langle n_{-1}\rangle\approx\langle n_{1}\rangle$, also see inset to Fig. 2(c)] so that the magnitude measures the concentration of $\Theta$. These results show that after the early dynamics settles down ($t\gtrsim 25t_{s}$) the function $R_{\Theta}$ approaches a negative real value, i.e. $\Theta\to\pi$. The $m=0$ component is unoccupied in the EP ground state, but maintains a small population [see inset to Fig. 2(c)] at late times due to heating from the quench. The $m=0$ component of the system is noisy (consistent with a thermalized gas, e.g. see Williamson and Blakie (2016b)) and the amplitude of the $R_{\Theta}$ correlation function is significantly reduced by these fluctuations. However, our results show that there is still a tendency for the spin-dependent interactions to lock the relative phase of the $m=\pm 1$ components relative to the $m=0$ component.

III.3 Late-time Universal coarsening dynamics

In addition to considering the emergence of local spin-nematic order we wish to examine the spatial dependence of the textures (domains) that develop and how these evolve in time. In Fig. 3 we visualize the system order in a region of a simulation soon after local order is established [Fig. 3(a)] and at a later time [Fig. 3(b)]. This visualization is performed by decomposing the spinor field at each simulation point according to Eq. (4) to obtain $\vec{u}(\mathbf{x})$ and $\theta(\mathbf{x})$. The results in Fig. 3 demonstrate that the spin-nematic and superfluid (i.e. global phase $\theta$) order tends to extend over larger length scales as time passes, showing that the system is coarsening toward an EP state with (quasi)-long range order.

To quantify the spatial dependence of the ordering we introduce the correlation functions

[TABLE]

for the spin-nematic and superfluid orders, respectively, evaluated at time $t$ after the quench. See Appendix B for more details about how these correlation functions relate to the atomic field operators.

To illustrate the use of these correlation functions, we consider the EP ground state spinor

[TABLE]

where the angle $\phi$ is associated with the spin-nematic order [i.e. $\vec{u}\propto\cos\phi\,\hat{\mathbf{x}}+\sin\phi\,\hat{\mathbf{y}}$], and the global phase $\theta$ is associated with the superfluid order. Taking $\theta$ and $\phi$ to be spatially dependent random variables, we use $\bm{\psi}_{\mathrm{EP}}$ to evaluate the correlation functions (30) and (31), yielding

[TABLE]

In practice we compute the spin-nematic order parameter correlation function as

[TABLE]

which includes averaging to improve the statistics of our results: $\langle\,\rangle_{t}$ denotes an average over trajectories (simulations with different seeding noise). The integral $\int d\Omega_{r}$ is an angular average in 2D position space (utilizing the isotropy of the system) and $l^{-2}\int d^{2}\mathbf{x}^{\prime}$ denotes spatial averaging. The convolutions are efficiently computed using fast Fourier transforms. We also apply these additional averaging steps when computing the $G_{\theta}(r,t)$ correlation function.

Results for the evolution of $G_{\phi}(r,t)$ are shown in Fig. 4(a). As time increases the correlation function is seen to decay more slowly, indicating that the in-plane spin-nematic order is extending over larger distances. We can investigate if the growth of this order exhibits dynamic scaling whereby the nematic domains are statistically self-similar at different times, up to an overall length scale that grows with time. This property often holds in the late time (when the domain sizes are much larger than microscopic length scales of the system) phase ordering dynamics of systems Bray (1994). To verify dynamic scaling we demonstrate that the correlation function collapses to a universal (time-independent) function under time-dependent rescaling of space, i.e. by showing that with an appropriate choice of $L_{\phi}(t)$ we have

[TABLE]

Results showing the collapse are presented in Fig. 4(b), where we have taken $L_{\phi}(t)$ to be the correlation length defined by the distance over which the correlation function decays to $\frac{1}{4}$ of its local value, i.e. $G_{\phi}(L_{\phi},t)=\frac{1}{4}G_{\phi}(0,t)$. The collapse is reasonably good except at short length scales ($r\ll L_{\phi}$) where the correlation function sharpens as $t$ increases.

The length scale $L_{\phi}(t)$ is not unique and can be multiplied by a constant and still yield correlation function collapse. However, as chosen $L_{\phi}(t)$ gives a reasonable characterization of the domain size555Domain size cannot be uniquely defined because the in-plane nematic order varies continuously. in the ordering EP system. From considering the evolution of $L_{\phi}(t)$ we can extract the dynamic critical exponent $z_{\phi}$ as $L_{\phi}(t)\sim t^{1/z_{\phi}}$, providing a key characterization of dynamic universality class of the system. In Fig. 4(e) we show the time evolution of $L_{\phi}(t)$ on a log-log graph and find that at late times ($t\gtrsim 10^{3}\,t_{s}$) this grows as $L_{\phi}(t)\sim[t/\ln(t/t_{0})]^{1/2}$, i.e. with a dynamic critical exponent of $z_{\phi}=2$ and logarithmic corrections. We find that the growth law exhibits a slight bulge (i.e. above the asymptotic growth law) extending from early times up until times of the order $10^{3}t_{s}$. We find that this correlates with the time period over which the magnetic fluctuations evolve appreciably in the system [see Fig. 4(f)], suggesting that the decay of magnetic fluctuations may set an important time scale for the system entering into the late-time coarsening regime (also see Kang et al. (2017)).

The $L_{\phi}(t)\sim[t/\ln(t/t_{0})]^{1/2}$ growth law we obtain here is the same form of growth known from the dissipative 2D XY model Yurke et al. (1993); Bray et al. (2000) (also see Kulczykowski and Matuszewski (2017)), and was established in early work considering the coarsening dynamics of smectic liquid crystal films Pargellis et al. (1992) (also see Nam et al. (2011); Singh et al. (2012); Singh and Singh (2013)). Singh et al. Singh et al. (2012) have predicted an analytic form of $H_{\phi}$ for nematic liquid crystals, which they have favourably compared to the results of Monte Carlo simulations using of a spin-nematic liquid crystal model Blundell and Bray (1992). We however, find that this result is not a good fit to the $H_{\phi}$ we obtain.

We can also consider the superfluid scaling in this system, with examples of the evolving $G_{\theta}$ correlation function shown in Fig. 4(c). We verify dynamic scaling in a similar way to the spin-nematic order by finding a length scale $L_{\theta}(t)$ such that we have correlation function collapse:

[TABLE]

Results showing this collapse are presented in Fig. 4(d), where again we have taken $L_{\theta}(t)$ to be the distance over which the correlation function decays to $\frac{1}{4}$ of the its local value. These results also reveal that the late-time superfluid correlation function $G_{\theta}$ has a similar shape to the spin-nematic correlation function $G_{\phi}$. By definition both correlation functions have the same local value, i.e. $G_{\theta}(0)=G_{\phi}(0)=\langle\mathcal{A}_{\perp}^{2}\rangle/n_{c}^{2}$. However, in general the superfluid correlation function decays more slowly and has a slightly longer characteristic length than the spin-nematic correlation function [e.g. see inset in Fig. 4(c)].

In Fig. 4(e) see that $L_{\theta}$ grows in a similar way to $L_{\phi}$, consistent with the same dynamical critical exponent, i.e. $z_{\theta}\approx z_{\phi}\approx 2$ (to within $\log$-corrections). Thus we find that the superfluid and spin-nematic order grow together in this system. This is different to recent results for the ordering of an EA ferromagnetic phase of a spin-1 condensate, which found that the superfluid order grows significantly slower than the spin order Bourges and Blakie (2017) (also see Williamson and Blakie (2017)).

It is conventional to also analyze the structure factors associated with the order parameter correlation function. The structure factor for spin-nematic order is defined as

[TABLE]

The structure factors also collapse with dynamic scaling according to

[TABLE]

where $\hat{h}$ is the Fourier transform of $H_{\phi}$ (39). Results for the $S_{\phi}$ structure factor are shown in Fig. 5. For $k$ vectors in the range $L_{\phi}^{-1}<k\ll\xi_{s}^{-1}$ (i.e. length scales between the microscopic healing length and the domain size) the structure factor exhibits a power law decay that is approximately of the form $k^{-3}$. This differs from the generalized Porod law result of $k^{-4}$ decay expected in 2D spin models Bray and Puri (1991); Singh and Singh (2013). The $k^{-3}$ decay law is also found for the first order structure factors (single-particle momentum spectra) in studies of binary condensates in relevant regimes Karl et al. (2013), and is analyzed in terms of turbulence scaling.

We can similarly define a superfluid structure factor $S_{\theta}$ from $G_{\theta}$. This structure factor has a similar collapse and power-law decay to what we have presented for $S_{\phi}(k)$.

III.4 Topological defects

It is of interest to consider HQVs, which are the topological defects supported by the EP order parameter. To illustrate the properties of HQVs we first consider a single HQV located at the origin. Away from the core the wave function is approximately of the form

[TABLE]

where we have set $\theta\to q_{\theta}\varphi$ and $\phi\to q_{\phi}\varphi$ in Eq. (35), $\varphi$ is the azimuthal angle about the core, and $\{q_{\theta},q_{\phi}\}$ are the winding numbers. In Eq. (49) we have also introduced the component windings

[TABLE]

where $\sigma_{m}$ denotes the net phase winding in the $m$-th component of the field. The $\sigma_{m}$ must be integer for the field to be single valued. The cases $\sigma_{1}=\pm 1$ (with $\sigma_{-1}=0$) and $\sigma_{-1}=\pm 1$ (with $\sigma_{1}=0$) define the four HQVs, corresponding to $q_{\theta}=\pm\frac{1}{2}$, $q_{\phi}=\pm\frac{1}{2}$, i.e. vortices with half-quantized values of the windings in $\theta$ and $\phi$ (see Fig. 6).

Much of our theoretical understanding of HQV dynamics has come from studies of miscible two-component condensates Öhberg and Santos (2002); Eto et al. (2011); Kasamatsu et al. (2016), which also support HQVs (also see Ji et al. (2008)). Notably, Eto et al. Eto et al. (2011) have shown that the interaction potential between two HQVs separated by a distance $R$ (for $R\gg\xi_{s}$) is of the form

[TABLE]

where

[TABLE]

with $(q_{\theta}^{(1)},q_{\phi}^{(1)})$ and $(q_{\theta}^{(2)},q_{\phi}^{(2)})$, [or $(\sigma_{1}^{(1)},\sigma_{-1}^{(1)})$ and $(\sigma_{1}^{(2)},\sigma_{-1}^{(2)})$] being the sets of winding numbers specifying HQV 1 and 2, respectively. For the case where both HQVs have winding in the same component (i.e. both having $|\sigma_{1}|=1$ or $|\sigma_{-1}|=1$ ) then $|\kappa|=\frac{1}{2}$ and the interaction is of the same form as that for $U(1)$ vortices in a scalar condensate. When the vortices occur in different components $\kappa=0$ and there is no long ranged interaction. However, a short ranged repulsive interaction is predicted, extending over a length scale comparable to the vortex core size Ji et al. (2008); Eto et al. (2011); Shirley et al. (2014). Two HQVs with opposite circulation in the same component (e.g. a HQV with $\sigma_{1}=1$ and a HQV with $\sigma_{1}=-1$) can collide and annihilate, as has been recently observed in experiments Seo et al. (2016).

Coarsening dynamics can be viewed in terms of the dynamics of topological defects of the order parameter which are generated in the early stages of the quench dynamics. The windings associated with these defects disrupt the order, and as they mutually annihilate order is able to extend over larger length scales. We show the locations of HQVs in Fig. 3, which reveals a qualitative relationship between the domain sizes and the vortex locations. To quantify the role of defects we detect the number of vortices in our simulations during the evolution. In practice we count the total number of integer phase windings in the $m=\pm 1$ components. In the early time dynamics not all vortices detected are HQVs, but we find that only HQVs persist at late times $(t\gtrsim 100\,t_{s}$). In Fig. 7 we show the averaged total number of vortices $N_{\mathrm{vort}}$ as a function of time. The number of vortices decreases as the coarsening progresses. We can compare these results to the characteristic length scales discussed in Sec. III.3. Crudely, if the characteristic length scale is taken to be the distance between vortices then we would expect

[TABLE]

We have added these results for the characteristic length to Fig. 7 verifying that the relationship (53) holds.

As we noted above a pair of $\sigma_{1}=1$ and $\sigma_{1}=-1$ HQVs (or a $\sigma_{-1}=1$ and $\sigma_{-1}=-1$ pair) evolve similarly to a vortex anti-vortex pair in a scalar condensate, and have the potential to mutually annihilate. In such a case each component vortex experiences a Magnus force which causes the pair to move with uniform velocity in a direction perpendicular to the line joining them. Such motion, without some other source of dissipation, does not lead to the vortices meeting and annihilating. This is in contrast to oppositely charged polar core spin vortices (the topological defects of the easy-plane ferromagnetic phase) that accelerate towards each other and annihilate Turner (2009); Williamson and Blakie (2016c). We expect that in our system dissipation will arise from the interaction between the vortices and the sound waves (spin waves) excited by the quench. However, recent results on HQVs suggest an additional dissipative mechanism even in the absence of spin waves: GPE simulations of a quiet binary condensate (without excitation) have found that such a pair of HQVs move together and annihilate (see Sec. IV of Kasamatsu et al. (2016)). This effect was observed to be dependent on the interaction parameter regime, only occurring for $\gamma>0.5$, where $\gamma$ is the ratio of the inter- to intra-species interaction in the binary condensate. In the spin-1 system666This mapping is made by neglecting the $\psi_{0}$ component in the spin-1 GPE. this parameter relates to the interaction parameters as $\gamma\approx(g_{n}-g_{s})/(g_{n}+g_{s})$. Since our main simulations presented are for $\gamma=0.5$, where this additional dissipative effect is expected to be negligible, it is of interest to see if our coarsening dynamics changes for a larger value of $\gamma$. To explore this issue we have conducted quench simulations for $g_{s}=g_{n}/12$ ($\gamma\approx 0.85$). The results for these simulations are roughly comparable to our main results in Fig. 4 (which are for $g_{s}=g_{n}/3$), and do not indicate that the coarsening proceeds at a faster rate. Nevertheless a better understanding of HQV dynamics, particularly in the spin-1 system at finite $q$ values, would be a valuable direction for future research. Also, a more detailed study of the dynamics and correlations between HQVs during the coarsening will be needed to illuminate the microscopic processes that are important in the system evolution (c.f. Schole et al. (2012)).

IV Conclusion

In this paper we have presented a theory for quantifying order formation in an anti-ferromagnetic spin-1 condensate. We have used this to study the dynamics of a quasi-2D system quenched into the EP spin-nematic phase. This topic has been of growing interest with a number of experimental developments motivating this work. This includes studies of correlations and spatial ordering in a quasi-one-dimensional system Bookjans et al. (2011); Vinit et al. (2013), and evolution of magnetic fluctuations and HQV formation in a quasi-2D system Kang et al. (2017). A key issue has been identifying appropriate observables to quantify spin-nematic order. This issue has been explored by Zibold et al. Zibold et al. (2016) who developed a novel measurement scheme to demonstrate spin-nematic order in the single mode regime Zibold et al. (2016). We motivate and define order parameters for the system to quantify the spin-nematic and superfluid order, and in doing this we have connected our formalism to quantities that have already been measured in experiments.

We have also studied the universal coarsening regime emerges at late-times after the quench. We evaluate the evolution of the order parameter correlation functions by averaging over an ensemble of large-scale numerical simulations and show that both types of order exhibit dynamic scaling, with a characteristic length scale that grows as $L\sim[t/\log(t/t_{0})]^{1/2}$. Our results also show that the coarsening is determined by the mutual annihilation of HQVs produced in the early stages of the quench. In experiments it may be difficult to directly measure the order parameter correlation function, whereas the average distance between HQVs (which can be directly imaged Seo et al. (2015, 2016); Kang et al. (2017)) will be a more convenient method to measure a characteristic length scale of order in the system.

Having developed and applied formalism for non-ferromagnetic ordering in a spin-1 system we open the door to other studies of ordering in spinor systems. This includes the rich array of spin order that emerge in higher spin systems (e.g. see Kawaguchi and Ueda (2012)).

Acknowledgments

The authors acknowledge support from the Marsden Fund of the Royal Society of New Zealand. PBB thanks Y. Kawaguchi for feedback on an early draft of the formalism, acknowledges useful discussions with L. Williamson, and thanks B.T. Wong for support of this research.

Appendix A Planar treatment of spin-nematic order

We can formulate our order parameters by considering the cartesian spinor field projected onto the plane:

[TABLE]

Recalling $\psi_{x}=\frac{1}{\sqrt{2}}(\psi_{-1}-\psi_{1})$, $\psi_{y}=-\frac{i}{\sqrt{2}}(\psi_{1}+\psi_{-1})$, we see that the planar treatment only depends on the $\{\psi_{1},\psi_{-1}\}$ spherical components of the spinor.

We now proceed to develop a mathematical description of the spin properties of the planar-spin system analogously to the three-dimensional treatment developed in Sec. II.2. We can decompose the planar spinor into two real planar vectors

[TABLE]

which are orthogonal and satisfy the normalization condition

[TABLE]

where $n_{\perp}=\vec{\psi}_{\perp}^{*}\cdot\vec{\psi}_{\perp}=n_{1}+n_{-1}$. We choose $\vec{u}_{\perp}$ to be the effective planar director and take it to be the longest vector, i.e. $|\vec{u}_{\perp}|^{2}\geq\frac{1}{2}n_{\perp}\geq|\vec{v}_{\perp}|^{2}$. We emphasize that the vectors $\{\vec{u}_{\perp},\vec{v}_{\perp}\}$ are not in general the projected versions of the three-dimensional vectors in Eq. (4) (e.g. projection of $\{\vec{u},\vec{v}\}$ does not preserve their orthogonality).

Because our vectors are 2D we can only obtain a $z$-component of the cross product, which yields the usual $F_{z}=|\psi_{1}|^{2}-|\psi_{-1}|^{2}$ magnetization density, i.e.

[TABLE]

The $m=0$ component projected out of the spinor prohibits us from quantifying the transverse magnetization. The singlet-amplitude to the planar system is defined as

[TABLE]

and we have the relation [c.f. Eq. (12)]

[TABLE]

We can construct a symmetric traceless tensor [i.e. the one introduced in Eq. (13)] as

[TABLE]

As noted in Sec. II.3 the elements of $Q$ in spherical spinor components are $Q_{xx}=\mathrm{Re}\{\psi_{1}^{*}\psi_{-1}\}=-Q_{yy}$ and $Q_{xy}=\mathrm{Im}\{\psi_{1}^{*}\psi_{-1}\}$, with $\mathrm{det}(Q)=-n_{1}n_{-1}$.

By inspection of Eq. (61) we see $\{\vec{u}_{\perp},\vec{v}_{\perp}\}$ are eigenvectors of $Q$ with eigenvalues $\lambda_{u}=\frac{1}{2}{n_{\perp}}-|\vec{u}_{\perp}|^{2}$ and $\lambda_{v}=\frac{1}{2}{n_{\perp}}-|\vec{v}_{\perp}|^{2}$, respectively. Given our convention to choose $\vec{u}_{\perp}$ as the longer vector we have that $\lambda_{u}$ is negative (i.e. the director corresponds to the lowest eigenvalue). Because the matrix is traceless the eigenvalues are given by $\pm\sqrt{-\mathrm{det}(Q)}$, i.e. $\lambda_{u}=-\sqrt{n_{1}n_{-1}}$ and $\lambda_{v}=\sqrt{n_{1}n_{-1}}$. The trace of $Q^{2}$ is then just the sum of the eigenvalues squared, and recalling the transverse alignment $\mathcal{A}_{\perp}=|\alpha_{\perp}|=\sqrt{2n_{1}n_{-1}}$, we obtain

[TABLE]

We also note that $Q$ can be written in the form

[TABLE]

where we have introduced $\varphi\equiv\frac{1}{2}\mathrm{Arg}(\psi_{1}^{*}\psi_{-1})$, i.e. $\psi_{1}^{*}\psi_{-1}=\frac{1}{2}\mathcal{A}_{\perp}e^{2i\varphi}$. Note that this has eigenvalues and eigenvectors:

[TABLE]

where the hats emphasize that these are unit vectors. We observe that the relative phase of the $\psi_{1}$ and $\psi_{-1}$ components directly determines the orientation $\varphi$ of the planar director $\vec{u}_{\perp}$. Note that this result is general for any spin-1 spinor, however for the particular case of the EP ground state (35) we have $\varphi\to\phi$, $\mathcal{A}_{\perp}\to n_{c}$.

Appendix B Correlation functions

Using the results of the previous section we can provide an alternative motivation for the correlation functions used in the paper. Firstly, we will consider the orientation of the director at two different points in space. For a spin model this might be characterized by a correlation function of the form

[TABLE]

where the inner product is squared to account for $\vec{u}$ and $-\vec{u}$ being the same. In terms of the fields our relevant quantity is the complex density $\Phi\equiv\psi_{1}^{*}\psi_{-1}=\frac{1}{2}\mathcal{A}_{\perp}e^{2i\varphi}$. Correlating this at two points in space we have

[TABLE]

which is identical to $G_{\phi}$ as defined in (30) if we normalize by a factor of $4/n_{c}^{2}$ .

From Eqs. (55) and (58) we see that the superfluid phase $\theta_{\perp}$ is related to the singlet-amplitude as

[TABLE]

where we can take $\theta_{\perp}=\frac{1}{2}\mathrm{Arg}(\psi_{1}\psi_{-1})$. Thus to correlate this superfluid order at two points we can consider the pairing-like field $\alpha_{\perp}$ at those two locations, i.e.

[TABLE]

Normalizing by a factor of $n_{c}^{-2}$ gives $G_{\theta}$ [Eq. (31)].