Domain percolation in a quenched ferromagnetic spinor condensate

We show that the easy-axis (EA) magnetic domains formed in a quenched ferromagnetic spinor condensate are described by percolation theory. We introduce a generalized spin rotation to vary the proportion of positive and negative EA domains, allowing us to explore domain percolation. Using simulations we investigate the finite-size scaling behaviour to extract the correlation length critical exponent and the transition point. We analyse the sensitivity of our results to the early-time dynamics of the system, the quadratic Zeeman energy, and the threshold condition used to define the positive (percolating) domains.

In this work we focus on the geometrical-statistical properties of the spin domains that form in a spinor condensate prepared in an unstable initial state using a generalized quench. These spin domains emerge in the easy-axis (EA) phase of a ferromagnetic spin-1 condensate, and prefer to have their magnetization either aligned (positive) or anti-aligned (negative) with the external magnetic field. The domains initially grow randomly, seeded by quantum or thermal noise. We analyse these domains in terms of percolation theory, canonically formulated to describe the behaviour of connected clusters in a random graph. It is useful to briefly introduce the central idea of percolation theory (e.g. see [25]). The so called site percolation problem involves a lattice in which sites are randomly and independently occupied with a probability p. The central question of percolation theory is for a given p, what is the probability that a path (of occupied sites) exists extending across the lattice? Such a path is also called a percolating cluster. For an infinite lattice there is a critical value p c (the percolation threshold) such that for p < p c a percolating cluster never occurs, while for p > p c it always occurs. It is found that a "geometrical" phase transition occurs at p c , such that the geometric properties of the clusters near p c are characterized by universal critical exponents [26].
The mapping of the continuous domains that occur in the EA phase of a spinor condensate onto a site percolation problem is one of the issues we discuss in this paper. Our basic approach is to consider positively magnetized regions as "occupied" sites, and to then analyse whether there exists a contiguous positive domain that spans the system. In order to explore the percolation behaviour it is necessary to vary the effective p value, i.e. relative proportion of the system occupied by positive domains. We demonstrate that this can be done using a generalized spin rotation to modify the initial state magnetization. The particular rotations we introduce are also chosen to reduce heating that will occur as the EA domains form.
Takeuchi et al. showed that the domains forming in the immiscibility transition of a binary condensate are described by percolation theory [27]. The EA ferromagnetic phase of a spinor condensate and immiscible phase of a binary condensate also have been revealed to have similar behaviour in phase ordering dynamics (c.f. [23,[28][29][30][31]). Spinor systems have some potential advantages for exploring phase transition dynamics, such as ease of preparing initial states and techniques for directly probing the spin degrees of freedom.
We now briefly outline the paper. In Sec. 2 we introduce the basic formalism for describing the dynamics of a spin-1 condensate. We present the initialization procedure we have developed for the quench dynamics to produce systems in which the relative proportion of positive domains can be controllably varied. We also introduce the numerical scheme we use to simulate the dynamics of a uniform quasi-two-dimensional (quasi-2D) system. Section 3 contains the main results of the paper. We define an effective occupation probability p in terms of the conserved system magnetization, and measure percolation by identifying domains that span or wrap around the periodic simulation grid. We study the percolation behaviour of the system using an ensemble of simulations to evaluate the probability that percolating domains occur. We perform a finite-size scaling analysis to demonstrate how the percolation probabilities change as the system size increases. This allows us to accurately extract the correlation length critical exponent ν, and extrapolate to the infinite system percolation threshold p c . We show that this threshold is independent of the percolation measure (i.e. spanning or wrapping domains), but is sensitive to the time after the quench, the quadratic Zeeman energy, and the condition used to identify the positive domains. In Sec. 4 we conclude and discuss the outlook for this work. Since our main results of the paper are presented for a uniform quasi-2D system we also briefly touch on the feasibility for observation in experiments and present a result for experimentally realistic parameters.

Hamiltonian and evolution
We study a homogeneous quasi-2D ferromagnetic spin-1 condensate. Such a system can be realized with 87 Rb atoms prepared in the ground state F = 1 hyperfine manifold and confined in an optical trapping potential. The system is then described by the spinor field ψ = (ψ 1 , ψ 0 , ψ −1 ), where ψ m is the component of the system in the m Zeeman sublevel.
The dynamics of this system is described by the spin-1 Gross-Pitaevskii equation (GPE) where q is the quadratic Zeeman energy shift, g n is the positive coupling constant for the density dependent interaction, and g s is the negative (i.e. ferromagnetic) coupling constant for the spin dependent interaction [32,33]. We have also introduced as the total (areal) density and the ν-component of spin density, where f ≡ (f x , f y , f z ) are the spin-1 matrices, and F = (F x , F y , F z ) is the spin density vector. The effect of the linear Zeeman shift is trivially removed by moving to a frame rotating at the Larmor frequency, and is neglected in Eq. (1). Evolution according to the GPE (1) conserves energy and total particle number. As the system is axially symmetric in spin space (invariant to spin rotations about z), the z-magnetization is also a constant of motion. The magnetization controls the effective proportion of positive and negative domains in the EA phase, and in the next subsection we consider how to alter M z in a way suitable for exploring the percolation transition.

Initial state preparation and quench
The spin domain formation dynamics we consider here can be implemented as a quench in which a single parameter of the system is changed and the system transitions towards a new equilibrium state. The initial state we consider is the uniform (zero-momentum) spinor where n 0 is the condensate density, and ϕ is an angle parameterizing the state. This miscible two-component state is the ground state for a spin-1 condensate with anti-ferromagnetic interactions (g s > 0) at q < 0. The quench of interest is to suddenly change g s to a negative value, whereby the two components are immiscible and EA domains will form. In principle the sign of g s can be changed using optical [34] or microwave-induced Feshbach resonances [35]. However, as these techniques have not yet been demonstrated in spinor experiments, we propose a different scheme to equivalently initialize the dynamics in a ferromagnetic condensate. While (5) is not an equilibrium state for a ferromagnetic condensate, it can be prepared starting from a polar state ψ P ∼ [0, 1, 0] T , which is the unmagnetized ground state for large positive q, by driving the atomic internal states using two subsequent electromagnetic pulses. The first pulse is the spin-1 rotation e −i π 2 fy , which produces the intermediate state g. see [36,37]). The second pulse is the pseudo-spin half rotation e −iϕσy of variable angle ϕ, performed on the m = ±1 levels, where σ y is the y-Pauli spin matrix. The result of both pulses is then the desired initial state (5). The second rotation could be driven as a two-photon transition between the m = ±1 states using microwave radiation detuned from the intermediate |F = 2, m = 0 state. We emphasize that the atomic physics toolbox of coherent manipulations presents various ways to quickly and reliably engineer (5), also see [38].
We briefly comment on the motivation for initializing the system to this particular initial state. First, this state has a magnetization of M z = N sin 2ϕ controlled by the angle of the spin rotation, where N = d 2 x n is the total number of particles. Second, particularly in the regime of small rotations (|ϕ| π/2) of interest here, this state is close to the state T , which undergoes less heating ‡ during the formation of EA domains (see discussion in Ref. [30]), and thus produces cleaner domains.

Simulations
It is useful to introduce q 0 = 2|g s |n 0 as characteristic spin energy, and use it to define the spin time § t s ≡ /q 0 , and the spin healing length ξ s ≡ / √ q 0 M . Simulations are performed using the spin-1 GPE with weak noise (representing the vacuum fluctuations) added to seed the dynamic instabilities. The noise is added to the polar condensate state as described in Ref. [39] and then the spin rotations described in the previous subsection are applied to this state to prepare the initial condition. The simulation is performed on a 2D square grid of spatial dimension L × L covered by a grid of N L × N L equally spaced points, with periodic boundary conditions. The GPE is evolved using the fourth order symplectic technique described in Ref. [40].

Percolating species and effective p value
The initial condition (5) is unstable for a condensate with ferromagnetic interactions. For q < 0 the ground state of the system is an EA ferromagnetic state, i.e. a phase in which the ground state condensate maximally aligns or anti-aligns along the spin z axis (i.e. the system prefers to have F z ≈ ±n 0 ). For a spatially extended system the initial condition will ‡ Atoms m = 0 liberate energy of −q when they spin mix into the m = ±1 levels to form EA domains. § This is the characteristic timescale over which spin can evolve when driven by the spin-dependent interaction. evolve into the EA phase by producing small domains, which can be labelled as positive and negative by the sign of F z . Our interest here is to characterise the properties of these domains. To map the EA system onto a percolation problem we will consider positive domains to be "occupied". In practice due to domain walls, and heating from the quench, the system is not perfectly polarized. We take the positive domains to be specified by the set of points where ≥ 0 is a constant that defines the threshold condition for a point to be included in a positive domain. We denote the total area of σ + as A + , i.e., the total area of positive domains. The ratio A + /L 2 quantifies the probability that a randomly selected point in the system is in a positive domain, and would be a useful choice for the occupation probability p used in standard percolation theory. However, A + is not a constant of motion and changes during evolution. We instead choose to use as an effective p value, defined in terms of the constant of motion M z . The value of M z and hence p is set by the initial state (5) using the spin rotation angle ϕ. For sharp domain walls and perfectly polarized condensate domains (i.e. everywhere F z = ±n 0 ), then this definition is equivalent to A + /L 2 .

Application of percolation analysis
In Figs. 1(a)-(c) we present examples of the F z density at a time of t = 50 t s after a quench with p = 0.5 (i.e. M z = 0) for three different realizations for initial noise. In the simulations the initial magnetization is uniform (F z = 0 for p = 0.5 in Fig. 1), but dynamic instabilities lead to the formation of EA domains. These domains become sufficiently well formed by around t ∼ 10 t s that after this time we can analyse their properties. To perform this analysis we apply the threshold condition (6) at every computational grid point to construct a binary image: points satisfying the threshold condition are assigned a value of 1 and all other points are assigned a value of 0. This allows us to calculate the total domain area as where N + is the total number of points with a value of 1, and ∆x = L/N L is the grid point spacing. We then apply the Hoshen-Kopelman algorithm [41] to the binary image to enumerate distinct positive domains (i.e. sets of occupied points in connected clusters). Our main analysis is to perform a quench simulation and then at some final time t quantify if any domain percolates in the following ways: Spanning domain: A domain which touches two opposite sides of the system (i.e. left to right, or top to bottom).
Wrapping domain: A spanning domain whose spanning ends overlap when periodic boundary conditions are imposed. We emphasize that the definition of a percolating domain is not uniquely defined, and other definitions could be used. Our results will show that while these two choices are different where N span (N wrap ) is the number of trajectories where at least one spanning (wrapping) domain was found. Since all wrapping domains are spanning domains, we have that P wrap ≤ P span . In Fig. 2(a) we plot the spanning and wrapping probabilities as a function of p, calculated using N traj = 590 simulation trajectories for each p value. These probabilities [fits shown in Fig. 2(a)] to determine the effective percolation threshold p eff c (L) and the width of the percolation transition ∆(L).
We can perform similar simulations and analysis to that in Fig. 2(a) but for simulations of different systems sizes L. From the fits to these results we can then determine the finite size scaling of p eff c (L) and ∆(L). The scaling result for the transition width allows us to extract the ν critical exponent. In standard percolation theory this critical exponent describes the divergence of the correlation length ξ as the percolation transition is approached: ξ ∼ |p − p c | −ν . The scaling relationship (9) usually furnishes a more accurate value for ν than would be obtained by measuring the correlation length divergence [42]. In Fig. 2(b) we plot ∆(L) −1 versus L for systems of 11 different sizes from L = 50ξ s up to 565.6ξ s (at fixed grid point spacing). The results for ∆(L) −1 demonstrate that for both spanning and wrapping percolation measures the transition sharpens (i.e. ∆ narrows) as the system size increases. The best power-law fits to both sets of results are shown as straight lines in Fig. 2(b), and the inverse slopes of both lines are the same within error bars giving ν = 1.34 ± 0.04 [c.f. Eq. (9)]. This result is in agreement with the value of 4/3 expected for standard percolation in 2D.
We can now examine the percolation threshold p c , i.e. the p value in the infinite system where the probability for finding a percolation domain abruptly changes from zero to unity. Finite size scaling predicts that the difference between the finite system percolation threshold and the infinite threshold p c scales as Using the value of ν extracted above, we plot p eff c (L) versus L −1/ν in Fig. 2(c). The results for both percolation measures are well fitted by straight lines. These lines are of different slopes, but extrapolate to the same y-intercept corresponding to the infinite system limit (i.e. L −1/ν → 0) of p c = 0.512 ± 0.002. This shows that while in a finite size system the various percolation measures are distinct, in the infinite system limit the percolation transition is sudden and insensitive to the particular definition we use for a percolating domain.

Generality of results
The analysis of the last subsection showed that the formation of percolating EA domains exhibits a well-defined transition in the infinite system, consistent with the standard 2D percolation transition. However, in order to better understand the generality of this result we investigate how the percolation behaviour changes with time  Fig. 3(c), which is seen to dip down and briefly oscillate at around t ≈ 15 t s , and then slowly increase back towards its initial value as time progresses. Spin mixing is more significant at small negative q values where the m = 0 level is energetically accessible [we discuss the role of q further in Sec. 3.3.2].
The tendency of the positive domain to percolate is related to the total domain area, and hence N + . The reduction in N + occurring when the domains first form, means that a higher initial N + population is needed to see percolation, or equivalently a higher magnetization (i.e. p value). For this reason p c is higher early on and decreases with time as the N + population increases.
It is also useful at this point to return to the definition of the p value given in Eq. (3.1). An alternative definition for the occupation probability p is the fraction of the system covered by positive domains, i.e.
We have explicitly given this quantity a t dependence to indicate that unlike Eq. (3.1) this quantity is not a constant of motion. An example of the evolution of A + is given in Fig. 3(c), showing that it rapidly grows from zero as the domains initially form, and then more slowly at later times as N + increases. We can adapt the analysis of our percolation results by plotting the probabilities at each time t against p (t) [c.f. Fig. 2(a)], where p (t) is the average over the N traj trajectories of the value of p at time t. Performing fits and finite size scaling analysis (as described in Sec. 3.2) we can then extract a percolation threshold p c and the critical exponent ν. The results of this analysis, applied to the same trajectories analyzed above in terms of p is also shown in Fig. 3(a) and (b). We see that the percolation threshold p c is now essentially constant at a value of p c ≈ 0.49. The critical exponent values are almost identical to the previous analysis.

Dependence on quadratic Zeeman energy
Our results thus far have been presented for the case of q = −0.3q 0 . The quadratic Zeeman energy shifts the m = ±1 Zeeman sublevels relative to the m = 0 sublevel. For small negative values of the quadratic Zeeman energy (−q 0 < q < 0) interactions are able to drive the spin-mixing of atoms back into the m = 0 level in the initial unstable dynamics. We observed this as a dominant source of the variation in p c with time in Fig. 3. For larger negative values (q ≤ −q 0 ) spin-mixing is energetically unfavourable. Indeed, an analysis of the quasi-particle excitations on the initial state (for ϕ = 0) is shown in Figs. 4(a) and (b) for the cases q = −0.3q 0 and −q 0 , respectively. There are three excitation branches with dispersion relations µ (k), which we label as µ = {0, 1, 2} (see [43]), with analytically known properties for the ϕ = 0 case [1]. For q = −0.3q 0 the magnon branches  Fig. 3(a)] and for q = −q 0 (circles). (d) Evolution of the relative population of m = 1 atoms (i.e. N + /N ) and the relative positive domain area (i.e. A + /L 2 ) for a simulation with L = 400 ξ s and p = 0.515. Inset: same results plotted with different range to show the initial (t < 10 t s ) behaviour. Results analyzed using the density threshold condition = 0.1.
magnons", which have amplitude in m = 0 sublevel. When these excitations grow on top of the condensate (noting the condensate is mainly in the m = ±1 sublevels) it leads to the formation of transverse magnetization F ⊥ = (F x , F y ). Significantly, the growth of these modes corresponds to spin-mixing of atomic population into the m = 0 sublevel. The 1 excitation branch consists of "axial magnons", which have amplitude in the m = ±1 sublevels. As these modes grow they cause the population in the m = ±1 to spatially separate into positive and negative domains, i.e. these modes directly initiate the immiscibility dynamics leading to the formation of the EA F z domains. For q = −q 0 [ Fig. 4(b)] only the axial magnons are dynamically unstable, and in this case there will be no spin-mixing.
An example contrasting the quenches to q = −0.3q 0 and −q 0 is shown in Fig. 4(d). We see that the N + population for the −q 0 quench is constant (c.f. −0.3q 0 , where spin-mixing causes N + to dip by about 20%), and as a result the domain area A + is seen to saturate faster, and vary more slowly in time, than the q = −0.3q 0 case. As a result the percolation analysis (in terms of p) shows that p c for the q = −q 0 case has a much weaker time dependence than for the q = −0.3q 0 case. Deeper q quenches will not lead to any additional changes in the dynamics, as the m = 0 level [already absent from our initial condition (5)], will remain frozen out from the dynamics. Indeed, in this regime the system behaves as a two-component system.

Sensitivity to the density threshold condition
Our analysis thus far has been based on defining positive domains via the density threshold condition (6) with = 0.1. We now consider the effect of changing . In Fig. 5(a) we show the F z density of a domain produced from a simulation with p = 0.5. Figures 5(b) and (c) show the binary image and largest domain constructed from that result using = 0.05 and = 0.1, respectively. Most noticeably the extent of the largest domain reduces with increasing . This is because often the domains are connected by tenuous regions with |F z | n 0 , which are sensitive to the threshold condition. In

Outlook and Conclusions
In this paper we have investigated the percolation of EA magnetic domains forming in spinor condensate. This could be investigated using a quench of the spin-dependent interaction, however we have instead proposed a scheme using a quadratic Zeeman quench and a genrealized spin rotation applied to a ferromagnetic spin-1 condensate, which will be feasible to implement in current experiments. By varying the conserved magnetization M z of the initial state, and hence the proportion of positive and negative EA domains, this system is able to explore the percolation transition. Using an ensemble of simulations of a quasi-2D system we have quantified the probability of percolation occurring as M z is varied, and for systems of various sizes. From these results we use finite-size scaling to extract the correlation length critical exponent, obtaining a value consistent with standard 2D percolation. We also use finite size scaling to extrapolate to the infinite system percolation threshold.
We have explored various aspects of the system dynamics, including the role of spinmixing in the early time evolution of the percolating clusters, showing that the percolation threshold is time dependent. We showed that this effect can be mitigated by directly measuring the positive domain area to define an occupation probability p , or by quenching to deeper values of the quadratic Zeeman energy, where spin-mixing is suppressed. At late times (t > 10 2 t s ), not considered here, the system will begin to phase order, and the typical size of domains will grow as l d (t) ∼ t 2/3 [21,39]. In this regime the system should exhibit the phenomenon of "phase-ordering percolation" [27], whereby the effective size of the system diminishes in time as L/l d (t). This would be an interesting direction for future investigation.
Our results indicate that it may be viable to study percolation in experiments. One Figure 6. Example of domains in a trapped system using 87 Rb parameters. (a) Normalized magnetization column densityF z ≡ dzF z (x) /n peak 2D , where n peak 2D is the peak (areal) column density of the initial condensate. (b) The binary image constructed usingF z > , with = 0.1 as the threshold condition. Positive domains are shaded white, except for the largest domain which is shaded red. Simulation for p = 0.5 case, with q = −0.5q 0 . Condensate of 6 × 10 6 atoms in flat trap: V (x) = V ρ ( x 2 + y 2 ) + 1 2 M ω 2 z z 2 , where V ρ (ρ) = 1 2 V 0 {tanh[(ρ − R)/w ρ ] + 1}, with ω z /2π = 10 3 s −1 , V 0 /h = 4 × 10 3 s −1 , R = 114 µm and w ρ = 7.6 µm. Using the peak initial condensate density n peak we have q 0 /h = 14.2 Hz, ξ s = 2.9 µm and t s = 19.9 ms. Initial preparation: a ground state condensate for potential V (x) is produced with all atoms in m = 0 state. This state has n peak 2D = 1.75 × 10 14 m −2 and n peak = 1.97 × 10 20 m −3 . Vacuum Wigner noise is added to planar momentum modes k ρ = (k x , k y ) on the condensate, restricted to 2 k 2 ρ /2M < ω z . This noise is projected onto the spatial region inside the trap, i.e. x 2 + y 2 < R. The spin rotations described in Sec. 2.2 are applied to this condensate with noise to produce the initial condition for the simulation. important reason is that the domains form in the early time dynamics [i.e. t ∼ O(10 t s )], which is a time scale easily accessible to experiments with 87 Rb (ferromagnetic) spinor condensates, which have a small value of |g s | and hence t s is large. Also it is feasible to manipulate the quadratic Zeeman energy on such time-scales with negligible heating (e.g. see [44]). A full study of the experimental system is outside the scope of this paper, nevertheless it is useful to explore the feasibility of observing EA domains in a realistic scenario. We indicate a result from such a calculation in Fig. 6, showing the domains formed and a domain that vertically spans the system. This result is constructed from the column density of a 87 Rb condensate in a flat-bottomed trap. Importantly this result shows that it is possible to get a reasonable number of domains forming, such that percolation properties will be nontrivial. It would also be interesting to consider a highly oblate harmonic trap. However, as the density decreases as we move towards the edge of such a trap the spin healing length, and hence the typical domain size, will also increase. Another area for future exploration is the role of finite temperature, whereby the initial state will have thermally occupied excitations. This will cause the domains to form more rapidly, but could also influence the nature of the domains that form.