Evolution and decay of an Alice ring in a spinor Bose-Einstein condensate

.


I. INTRODUCTION
Owing to their persistence through phase transitions and dissipation, topological defects are important in multiple areas of physics such as in cosmology [1] and condensed-matter physics [2].The Alice string is a topological line defect, a vortex line, that was initially studied in the context of field theories [3,4].It has an unusual feature of changing the sign of a monopole charge that encircles it [5][6][7][8].It is possible for an Alice string to appear also as a closed loop, an Alice ring (AR), as a result of a continuous deformation of a monopole defect [9,10].
Bose-Einstein condensates (BECs) with an internal spin degree of freedom provide an ideal platform for studying topological defects due to their versatile order parameter space and symmetries [11].The multicomponent order parameter fields of spinor BECs are able to host a wide variety of topological defects, such as coreless vortices [11], Dirac monopoles [12,13], skyrmions [14,15], and knots [16,17].In particular, the order parameter space of the polar phase of the BEC with hyperfine spin F = 1 supports vortices, quantum knots, and point-like monopole defects [18,19], analogous to the 't Hooft-Polyakov monopole [20,21].We refer to these topological monopole defects as isolated monopoles since they are free of any other attached defects, unlike the Dirac monopole.
It is energetically favorable for the isolated monopole to decay into a closed vortex ring that conserves the topological monopole charge [22].Consequently, such a vortex has a winding number of 1/2 for the scalar phase, allowed by the polar order parameter ζ p = e iϕ d that is invariant under the simultaneous substitutions of d → − d and ϕ → ϕ + π , where d ∈ R 3 is the unit-length director vector and ϕ ∈ [0, π ) is the scalar phase.Along a path encircling once the half-quantum vortex, the director continuously changes its direction, and the scalar phase winds π , rendering the vortex ring analogous to the AR considered in particle physics [23].
The long-time evolution of an isolated monopole in 87 Rb spin-1 BEC has previously been studied in Refs.[24,25], though in a quadrupole magnetic field or without any external magnetic fields.In both cases, the condensate eventually experiences a polar-to-ferromagnetic phase transition caused by dynamical instabilities driving the spinor degrees of freedom towards the ferromagnetic ground-state phase of 87 Rb.This change of magnetic phase inevitably leads to the decay of any produced ARs since the ferromagnetic order parameter does not support such a topological defect.In addition, even though the produced AR is observed in the previous literature [26], its evolution is not studied in detail.In Ref. [27], the stability and dynamics of an AR was studied in the case of 23 Na, the ground-state phase of which is polar.Their work was also conducted in the absence of magnetic fields, which prevents examining the effects of the Zeeman shifts on the decay dynamics.In addition, stationary-state AR solutions were used as initial conditions, which led to two solution branches with different ring radii.Thus, it remains open whether a homogeneous constant magnetic field can be used to improve the stability of ARs in polar BECs.Long-lived ARs would be of great interest thanks to their peculiar topological features.
In this paper, we present the evolution and decay dynamics of ARs in purely optically trapped 87 Rb spin-1 BECs in the presence of a homogeneous magnetic field.Since the total magnetization of the BEC is conserved over the time scale of the simulations, the minimization of energy to find the effective ground-state magnetic phase does not depend on the linear Zeeman energy, and the external-magneticfield-induced quadratic Zeeman energy shifts the ground-state phase of the 87 Rb condensate to polar and hence alters the decay dynamics of the polar phase.Instead of a stationary-state AR solution, we produce the AR by simulating the creation and decay of the initial isolated monopole.With a purely polar initial state and perfectly centered monopole, we are able to verify the existence of the AR after an evolution of over 160 ms.Interestingly, near the end of this period we find a second AR to emerge from the boundary of the condensate.We also study in detail the effects of the quadratic Zeeman shift, added noise, and the location of the monopole on the dynamics and decay of the AR which emerges as the decay product of the initial monopole.We find that by creating the initial monopole near the condensate boundary, we are able to change the winding direction of the emerging AR, compared with the centered case.This allows us to create two ARs, with opposite winding directions, which eventually annihilate each other.Experimental nonidealities, such as noise and inaccuracy in the monopole positioning, are observed to shorten the lifetime of the AR, though the ring is still clearly noticeable beyond 80 ms after its creation.These long-lived ARs suggest a possibility to theoretically and experimentally demonstrate their topological charge negation property on monopoles that travel through them [3].
Furthermore, we pay attention to the observation in Ref. [27] that the instability of the AR may be caused by the error arising from the finite spatial discretization.It has been shown in the case of the scalar Gross-Pitaevskii equation [28,29] that more isotropic discretization grids than the widely used evenly spaced Cartesian grids can improve the accuracy and the reliability of the numerical solution.Therefore, we utilize a spatial discretization consisting of a truncated octahedral tiling in order to decrease the AR decay caused by the discretization error.The used numerical method is also verified to conserve the total magnetization to ensure the polar ground-state phase of the BEC in the presence of a strong enough homogeneous magnetic field.
The rest of this paper is organized as follows: In Sec.II, we introduce the mean-field theory of spin-1 BECs and the characteristics of the polar order parameter.Section III describes our methods for the AR creation and numerical temporal integration.The results of the simulations are presented in Sec.IV.Finally, Sec.V summarizes and concludes this article.

A. Mean-field theory
In the mean-field approximation of a dilute spin-1 Bose-Einstein condensate, the order parameter may be expressed as (r, t ) = √ n(r, t )ζ (r, t ), where n is the particle density of the condensate and ζ denotes a three-component complexvalued spinor that satisfies ζ † ζ = 1.The components ζ m are indexed by the magnetic quantum number m ∈ {+1, 0, −1}.The temporal evolution of the order parameter at adequately low temperatures is described by the Gross-Pitaevskii meanfield equation where F = (F x , F y , F z ) is a vector of the standard dimensionless spin-1 matrices, F = † F is the local average spin, is the three-body loss term, c 0 and c 2 are effective atom-atom interaction strengths, and ĥ0 is the single-particle Hamiltonian given by ĥ0 (r, t where m is the atom mass, g F is the Landé g factor, μ B is the Bohr magneton, B is the external magnetic field, and q is the strength of the quadratic Zeeman shift.The harmonic optical trapping potential is given by V = m(ω r x 2 + ω r y 2 + ω z z 2 )/2, where ω r and ω z are the radial and axial trapping frequencies, respectively.The spin-1 matrices F α (α = x, y, z) are given in the zquantized basis, and hence the spinor ζ is presented in the eigenbasis of F z .In the analysis of topological defects, it is convenient to express the spinor also in the Cartesian basis, in which, the components can be obtained by the unitary transformation ⎛ ⎜ ⎝ and the spinor and the order parameter are denoted by ζ and = √ n ζ , respectively.

B. Polar order parameter manifold
For spin-1 BECs in the absence of external magnetic fields, two different magnetic phases can be identified, ferromagnetic and polar, which are determined by the local average spin magnitudes | F | = 1 and | F | = 0, respectively.In the polar phase, the spinor can be generally expressed in the eigenbasis 023104-2 of F z as where ϕ is the scalar phase and {d α } are real-valued numbers satisfying d 2 x + d 2 y + d 2 z = 1.This representation of the spinor ζ p in the eigenbasis of F z can be transformed using Eq. ( 3) such that the Cartesian-basis spinor becomes ζ p (r, t ) = e iϕ(r,t ) d(r, t ), where d = (d x , d y , d z ) is a realvalued unit vector.The polar spinor ζ p is invariant under SO(2) rotations about the director d, and it also possesses a discrete Z 2 spin-gauge symmetry of e iϕ d = e i(ϕ+π ) (− d), which allows an apparent scalar-phase discontinuity of π in a continuous order parameter field.The resulting polar order parameter space is M p = [U (1) × S 2 ]/Z 2 , the first and second homotopy groups of which are isomorphic to the additive group of integers π 1 (M p ) ∼ = π 2 (M p ) ∼ = Z.This permits the order parameter field to host both monopoles and vortices as topological defects, with full and fractional-integer charges, respectively.

C. Alice ring
An isolated monopole in the polar phase of a spin-1 BEC, with the topology defined by the director d, is schematically illustrated in Fig. 1(a).Such a topological defect causes the particle density to vanish at its core.This has a high energy cost and hence promotes the defect unstable against decay into a closed vortex loop, which preserves the topology of the initial monopole far from its core.Such a vortex loop may have a lower energy cost of filling its core with finite particle density in the ferromagnetic phase than the corresponding cost of forcing the density to vanish at the core of the vortex.A schematic illustration of such a closed ferromagnetic-core vortex loop, an Alice ring, is shown in Figs.1(b) and 1(c).
The figure illustrates how the initial monopole topology of the director field d is preserved far from the torus-shaped vortex, and hence on any closed poloidal path about the vortex core, such as L in Fig. 1(b), the director d undergoes a full reversal of direction, quantifying the topological charge of the AR as 1/2.Additionally, on any such path similar to L, the scalar phase ϕ of the polar order parameter also experiences a winding of π , rendering the ring a half-quantum vortex.Such a path must always pass through a point where the director d suddenly flips its sign, and the scalar phase ϕ changes by π radians, although the exact location of this is not uniquely determinable since it is not physical and depends on the choice of gauge.In Fig. 1 the gauge has been chosen such that this apparent discontinuity is located in the middle of the AR.
In the situation described above, the order parameter holds a rotational symmetry about the vertical ring axis, as presented in Fig. 1(c).
If a half-quantum vortex is line-like, i.e., not wrapped into a loop as illustrated in Fig. 2, it may be considered to be an Alice string.Any path about this line-like defect, such as L in Fig. 2, holds identical properties to those of L described above.In Fig. 2, the gauge has been chosen such that the

A. Monopole creation
The creation of the initial monopole defect follows the method originally developed for Dirac monopoles in Ref. [12] and is extended to isolated monopoles in Ref. [18], where the internal degrees of freedom of the order parameter 023104-3 adiabatically follow the changes in the external magnetic field.In summary, the condensate is initially in the polar-phase ground state = √ ne iφ (0, 1, 0) , in an external magnetic field of the form B(r, t ) = B b (t ) + B q (r, t ), where B b (t ) = B b (t )ẑ is a homogeneous bias field and B q (r) = b q (x x + yŷ − 2zẑ) is a quadrupole magnetic field.Initially, the bias field strength is B b = 1.0 G and the quadrupole gradient b q is linearly increased from zero to 4.3 G/cm in 10 ms.The bias field is subsequently decreased to 45 mG in 10 ms, placing the point at which the magnetic field vanishes 52 µm above the center of the condensate and keeping the director d approximately aligned with ẑ.The magnetic field zero is then brought to the center of the condensate by reducing the bias field linearly to zero at the rate dB b /dt = −0.25 G/s for 180 ms, causing the director to adiabatically follow the magnetic field and ideally result in a texture d = (x, y, −2z)/ x 2 + y 2 + (2z) 2 .Note that the parameters we employ here are experimentally demonstrated in Refs.[18,19].Immediately after the monopole creation process, we extinguish the quadrupole field in 200 µs and simultaneously ramp up the bias field to 1.2 G in 500 µs.Subsequently, the condensate is allowed to evolve in the presence of the harmonic optical trap and the homogeneous bias field B b .
The atom number in the simulated condensate is N = 2.5 × 10 5 and the radial and axial optical trapping frequencies are ω r = 2π × 120 Hz and ω z = 2π × 160 Hz.We set the three-body loss term to = h × 2.9 × 10 −30 cm 6 /s and the atom-atom interaction strengths to c 0 = (4π h2 /3m)(2a 2 + a 0 ) and c 2 = (4π h2 /3m)(a 2 − a 0 ), where a 0 = 101.8× a B and a 2 = 100.4× a B are the scattering lengths for 87 Rb and a B is the Bohr radius.Depending on the scenario we study, the strength of the quadratic Zeeman shift is set to either zero, q = 0, or to q = (g F μ B ) 2 / E hf , where E hf ≈ (6.8 GHz) × 2π h is the 87 Rb hyperfine energy splitting and g F = −1/2.

B. Numerical methods
In the numerical integrator, we scale the time, energy, and spatial distances into the units of hω r , 1/ω r , and a r = √ h/mω r , respectively, and numerically solve the resulting dimensionless Gross-Pitaevskii equation (GPE) using discrete exterior calculus (DEC) [30].We briefly introduce our utilization of DEC to the GPE in the Appendix, whereas a detailed description of this matter is given in Ref. [29].In our DEC method, the spatial domain is discretized using a pair of meshes: a primal mesh and its dual mesh.We choose the body-centered cubic (BCC) grid as the primal mesh and its Voronoi diagram, the truncated octahedral grid, as the dual mesh, and present the order parameter on the latter.As a point of comparison, we run the simulations also using cubic grids for both meshes.The ferromagnetic-phase ground state is solved using the imaginary-time evolution method in combination with DEC, without any external magnetic fields.Subsequently, the ground state is forced to the polar phase by replacing the spinor with ζ = (0, 1, 0) .
During the evolution, the order parameter does not persist purely in the polar phase, though it is still possible to find a unique director and phase also in the mixed phase, which have identical properties to those in the pure polar phase, as required by the continuity of the order parameter.We extract ϕ and d by first decomposing the spinor in the Cartesian basis into two real-component vectors as ζ = u + iv.Then, we find a phase ϕ , such that for the gauge transformation e iϕ (u . For a purely polar order parameter, we have b = 0, d = d/|d|, and ϕ = ϕ , and hence, we use these quantities for the scalar phase and the director also in the case of mixed phases [31].
The numerical simulations are carried out using a cubic computational domain with the side length of L = 24 × a r and spatial discretizations with 112 3 × 12 and 256 3 degrees of freedom for the truncated octahedral mesh and for the cubic mesh, respectively.

A. Decay dynamics and lifetime
First, we examine the evolution of the AR using the physical parameters presented in Table I and I.
field is extinguished.The decay of the initial monopole into the AR is presented in Fig. 3, which shows that at T = 0, the director is well aligned with the quadrupole field and a torusshaped domain with elevated local average spin magnitude (| F | = 0.8) has already emerged around the center of the BEC.For the director and scalar phase, Fig. 3 focuses only on the x = 0 plane, but the structure is essentially symmetric about the z axis.After T = 0, the radius of the ring and | F | start to increase, and at T = 4 ms, the vortex core has been filled with almost pure ferromagnetic phase (| F | > 0.95).The ferromagnetic filling is energetically favorable over forcing the particle density to zero at the vortex core since for 87 Rb the spin healing length ξ s is significantly greater than the density healing length ξ n (ξ s ≈ 14.7 × ξ n ). Figure 3 also clearly illustrates that on a closed poloidal path about the core, the scalar phase ϕ winds and the director d rotates by π radians.Throughout the rest of our analysis, we utilize the ferromagnetic core, in combination with the scalar phase and director, to verify the existence of the AR.
First, we examine the long-time evolution in these ideal circumstances where the initial state 0 = √ n g (0, 1, 0) , where n g is the density of the ground state that is free of noise, and the monopole is placed exactly at the center of the BEC. Figure 4 shows the ferromagnetic core of the AR vortex at different points in time, which reveals that during the evolution, the AR travels back and forth along the z axis and that its radius varies in time.During the first 90 ms of the evolution, the AR stays within a harmonic oscillator length a r from the origin and the radius of the AR increases from 0.8 × a r to 1.5 × a r .Despite the changes in the position and size, the ferromagnetic core retains its ring-like form for over 160 ms, after which the ring breaks at a single point and the vortex extends into a line-like defect, with both ends at the condensate boundary still possessing the half quantum of phase winding for the full length of the vortex, i.e., the Alice ring becomes a curvy Alice string, illustrated in Fig. 5.The figure also shows how one end of the Alice string starts to align with the z axis.
The decay of the AR can be observed also in the evolution of the total magnetization M = F dr, visualized in Fig. 6, where the x and y components of M start to oscillate around the zero value approximately at the same time the AR decays.The figure also verifies the conservation of the z magnetization in the numerical method, which is important to avoid the transition to the ferromagnetic phase.
Interestingly, just before the decay of the initial AR into an Alice string, another AR emerges from the boundary of the condensate as visualized in Fig. 7.When compared to the original AR, the second ring has an opposite topological charge, visible from the rotation direction of the director d on  a poloidal path about the vortex.Despite this, the two ARs have identical scalar-phase winding directions, which keeps them apart and therefore tends to protect them from annihilating each other.After the original AR decays into the previously mentioned Alice string, it subsequently merges with this newly produced ferromagnetic domain near the condensate boundary and causes them together to break into arbitrary sections with no clearly identifiable structure.At this point of time, we consider the half-quantum vortex to be fully decayed.
Next, we simulate the AR in conditions similar to their experimental realization [26].We perturb the initial state with noise as 0 = √ n g (0, 1, 0) + 0.1 × (X, Y, Z ) , where X, Y, Z ∈ C with real and imaginary parts being random numbers according to the standard normal distribution.The noise amplitude of 0.1 is chosen empirically in an attempt to maximize the matching between the simulated results and the experimental observations in Ref. [26].The noise populates dynamically unstable modes which draw more population exponentially in time, and hence the lifetime of the AR can be expected to decrease only logarithmically with respect to the initial noise amplitude.After this, the order parameter is normalized to match the atom number used in the experiments [26] and to prevent its fluctuation with the random seed.Furthermore, we create the initial monopole randomly inside the origin-centered ball of radius a r .We observe that these factors prompt arbitrary ferromagnetic-phase domains to emerge from the condensate boundary, which causes the AR to decay sooner than in the ideal case considered above.The displacement of the initial monopole is also observed to affect the orientation of the AR and lead to a bend of the ring axis away from the z axis as shown in Fig. 8.Given the employed three-body loss rate , a characteristic number for 87 Rb spin-1 BECs, the particle number Ñ (t ) = n(r, t )dr essentially vanishes approximately at T = 3 s, and causes Ñ to decrease only by 5% of the initial number N r .Random noise with a magnitude of 10% of the order parameter amplitude is added to the initial state and the monopole is placed randomly in the origin-centered ball of radius a r by the creation ramp.See Table I for parameter values.
during the lifetime of the AR.The simulation was run also without the three-body loss term, but no significant difference in the AR evolution or lifetime was observed.We also study the effects of the quadratic Zeeman effect and the spatial discretization on the lifetime and decay dynamics of the AR.The AR lifetime is not found to be significantly longer when the strength of the quadratic Zeeman shift is set to zero, hence the decay is unlikely to be caused by the quadratic Zeeman effect.On the other hand, the evenly spaced Cartesian grid, commonly used with the time-splitting spectral methods [32], is found to cause faster AR decay such that, for example, in the ideal conditions, a lifetime of only 100 ms is observed.

B. Winding direction of the scalar phase
Let us examine the noise-free case.By comparing the initial state of the AR in Fig. 3 for T = 4 ms to the decayed states in Figs. 5 and 7 for T > 160 ms, we observe that the half-quantum winding of the scalar phase changes its direction during the evolution.This spontaneous flip in the winding direction takes place approximately at the half-lifetime point T ≈ 82 ms and is illustrated in Fig. 9.The figure shows how at first, the winding is such that the scalar phase increases along the positive z direction when traveled through the ring.Next, the AR shrinks significantly, after which it grows back to its original size, but with a winding such that the scalar phase increases in the negative z direction when traveled through the AR.Furthermore, closer examination of the dynamics of the topological defect during the flip reveals that, in between the different scalar phase winding directions, the AR contracts back into a monopole, presented in Fig. 10.These observa- tions motivate us to further investigate the factors determining the winding direction and the cause for its flip.A small displacement, not exceeding a r from the origin, of the initial monopole is not noticed to affect the winding direction.It seems to rather be determined by the noise perturbation, or more precisely, by the spin polarization of the initial state: how much the local spin magnitude | F | differs from zero seems to be correlated with the winding direction of the emerging half-quantum vortex.Therefore, in the experimental conditions, the AR possesses a winding opposite to the ideal case already from the beginning and does not experience the direction change during the evolution.From these observations it appears that in the ideal conditions, the spin polarization of the AR-supporting order parameter increases during the evolution, and once it reaches a certain threshold, the winding direction flips.The winding direction can also be identified from the integrated m = 0 particle density on the y = 0 plane as shown in Fig. 11.The different winding directions cause either a V-shaped or a -shaped gap to the density profile, which we will subsequently call positive and negative windings, respectively.
Even though the aforementioned small displacement of the initial monopole does not noticeably affect the winding, the case is different with more substantial distances away from the origin.In particular, if the monopole is displaced along the z axis such that the AR is produced near the condensate boundary, the winding may be opposite to that of the wellcentered case.This happens, for example, if the preferred winding direction is such that it would cause the AR to travel upwards, but the AR is created near the top of the condensate so that it does not have room to travel up.This causes the AR to have its winding to favor downwards movement.Note that also for the well-centered AR, the direction of movement in Fig. 4 agrees with the winding directions shown in Fig. 9.

C. Annihilation
In Sec.IV A, we mentioned the emergence of the second AR with the opposite topological charge before the decay of the original AR.However, during their coexistence, we observe that the ARs arising from the ideal conditions have identical winding directions of the scalar phase (Fig. 7).By utilizing the observations obtained in Sec.IV B, we are also able to produce two coexisting ARs with both opposite charges and opposite scalar-phase windings.This is achieved by creating the initial monopole close to the top boundary of the cloud, which forces its decay product AR to flip its winding with respect to the z = 0 case.The second AR still similarly emerges near z = 0, holding its original winding.We use such noise to perturb the initial condition that is known to cause a -shaped m = 0 density profile when the monopole is placed at the origin, i.e., negative scalar-phase winding.The initial monopole is then placed near the top boundary of the condensate cloud at z = 2.42 × a r , which causes the produced AR to have positive winding instead.Again, in the decay dynamics, a second AR emerges from the side boundary but now, contrary to the original AR, the second one possesses the negative winding that an initially centered monopole would also have.
The opposite winding directions of the two ARs cause them to approach each other until their ferromagnetic cores eventually merge together and the opposite scalar phase windings and director-determined topological charges cancel each other, which subsequently releases the merged ferromagnetic cores and allows them to return to the polar phase.This intriguing annihilation process of the AR and anti-AR is presented in Fig. 12, where at T = 10 ms, we observe the scalar-phase winding to be zero on a poloidal path about the point where the ferromagnetic cores have merged.The phase change can also be observed as a drop in the evolution of the expectation value of the local average spin magnitude | F | = n| F |dr/N presented in Fig. 13.The noise perturbation in the initial condition and the extreme placement of the monopole cause the decay to happen a lot quicker than in our other simulations, already at around T = 10 ms.

V. CONCLUSION
We have numerically studied the long-time evolution and decay dynamics of Alice rings in 87 Rb spin-1 BECs.The evolution was studied in the presence of an external homogeneous magnetic field using a natural-crystal-structure spatial discretization in our numerical solver to maximize the lifetime of the Alice ring.The numerical simulations show that with an initially unperturbed ground state, an Alice ring created at the center of the condensate has a lifetime of over 160 ms.In conditions closer to those of experimental setups, the Alice ring decays after a 90-ms evolution.A detailed inspection of the decay dynamics also revealed the emergence of a second Alice ring into the condensate from its boundary.The director field of the emerging Alice ring has an opposite monopole winding from that of the original ring.In this case, the two Alice rings have an identical winding direction of the scalar phase, and the original ring decays by breaking into an Alice string instead of annihilating with the coexisting anti-Alice ring.The scalar-phase winding of the centered Alice ring is observed to be correlated with the magnetic polarization of the order parameter that supports the initial monopole defect.The winding is opposite that of the centered case if the Alice ring is created near the condensate boundary such that its natural traveling direction would otherwise be blocked.This causes the anti-Alice ring to emerge with an opposite scalar-phase winding, which enables the annihilation of the two Alice rings inside the condensate.In future studies, a more thorough investigation of the cooperative action of noise and three-body losses to the dynamics and lifetime of the Alice ring may be conducted.Most interestingly, the long-lived Alice rings may be utilized to demonstrate their intriguing property of negating the charge of a monopole that passes through them.

ACKNOWLEDGMENTS
We acknowledge financial support from the Academy of Finland through its Centre of Excellence in Quantum Technology (Grant No. 336810) and computational resources from CSC -IT Center for Science, Finland.We thank David Hall, Tuomas Ollikainen, and Toni Annala for useful discussions.

APPENDIX: DISCRETE EXTERIOR CALCULUS
Discrete exterior calculus provides a discrete extension of differential geometry and exterior algebra, and in this description of its utilization, we assume some prior knowledge with its concepts as defined in Ref. [30].A more detailed description of the utilization of discrete exterior calculus to general wave propagation problems is given in Ref. [33] and particularly to the Gross-Pitaevskii equation in Refs.[28,29].
We begin by regarding the order parameter as a complex vector-valued discrete differential 0-form ψ, which allows us to write the Laplacian in the single-particle Hamiltonian [Eq.( 2 is a matrix representing the discrete exterior derivative that operates on discrete differential 2-forms, and 2 and 3 are matrices representing the discrete Hodge stars that operate on discrete differential 2-forms and 3-forms, respectively.
We discretize the spatial domain using a pair of meshes consisting of a primal mesh and a dual mesh.The bodycentered cubic grid is used as the primal mesh and its Voronoi diagram, the truncated octahedral grid, as the dual mesh.Because of the choice for the dual mesh, its edges are always orthogonal to the faces of the primal mesh, which renders the discrete Hodge stars diagonal matrices.The discrete differential 0-form ψ corresponds to the nodes and the discrete differential 1-form u = d 2 ψ to the edges of the dual mesh.The discrete codifferential − 3 d 2 −1 2 maps the discrete differential 1-form u first to primal mesh faces as a discrete differential 2-form, and then through a 3-form on primal mesh bodies back to the dual nodes, which finalizes the full discretization of the Laplacian.

FIG. 1 .
FIG. 1. Schematic illustrations of (a) a monopole and (b and c) an Alice ring.In all panels, the greenish ellipsoid represents the polar-phase BEC cloud and the red domain is the ferromagnetic core of the Alice ring.In (a) and (b), the cloud is cropped away for y > 0 and the director field d is visualized with black arrows.In (a), the black dot indicates the location of the zero density.In (b), a closed poloidal path L about the vortex core is colored with the scalar phase ϕ according to the given color bar.The director d is also shown at evenly spaced points on the path L and colored with the scalar phase ϕ.In (c), the BEC cloud is cropped away for z > 0.

FIG. 2 .
FIG. 2. Schematic illustration of an Alice string.The greenish ellipsoid represents the polar-phase BEC cloud and the red domain is the ferromagnetic core of the Alice string.The cloud is cut in half along a horizontal axis and the director field d is visualized on the intersection with black arrows.A closed poloidal path L about the vortex core is colored with the scalar phase ϕ according to the provided color bar.The director d is also shown at evenly spaced points on the path L .

4 FIG. 3 .
FIG. 3. Decay of a well-centered initial monopole into an Alice ring in the ideal noise-free case.Cross section (x = 0 plane) of the BEC order parameter with significantly nonvanishing particle density for (a and b) T = 0 ms and (c) T = 4 ms, with the director d and the scalar phase ϕ of the polar phase represented by the black arrows and the background color map, respectively.In (b), the monopole is shown in a zoomed-in and rotated (about the y axis) view with a three-dimensional (3D) visualization of the domain of the elevated local average spin magnitude (| F | = 0.8) as a bright red torus (color not denoting ϕ).The semitransparency of the red surface denotes that the domain is not purely ferromagnetic.In (c), the view angle is as in (b) for T = 4 ms with the fully opaque bright red torus-shaped 3D domain presenting the almost pure ferromagnetic phase (| F | 0.95).In (a) the field of view is 14 × 10 a 2 r .The field of view of (b) and (c) is denoted by the black-border isosceles trapezoid in (a).The white lines in the scalar phase color maps denote the instantaneous jumps of π and the simultaneous reversals of d.The physical parameters used in the simulations are presented in TableI.

FIG. 4 .
FIG. 4. Ferromagnetic cores (| F | > 0.95) of the Alice ring shown by different colors for different evolution times T after the monopole creation in the case of a well-centered monopole and no initial noise.The 3D toruses are shown directly along the y axis and hence their ring-like shapes are not fully visible.See Table I for parameter values.

FIG. 5 .
FIG. 5. Alice string at T = 180 ms, as a product of an Alice ring breaking off at a single point.In (a), the director d and the scalar phase ϕ are presented on the y = 0 plane at regions of significant particle density, with black arrows and background color map, respectively, and the almost pure ferromagnetic phase (| F | > 0.95), i.e., the core of the curvy Alice string, is represented by the bright red 3D domain (color not denoting ϕ).In (b), the same ferromagneticphase core of the string as in (a) is colored by the scalar phase ϕ and shown from a different view angle.The field of view in (a) is 12 × 9 a 2 r and in (b) is 6 × 6 a 2 r .The simulation was conducted in the ideal conditions.See Table I for parameter values.

023104- 6 FIG. 8 .
FIG. 8. Long-lived Alice ring in experimental-like conditions at T = 84 ms.The director d (black arrows) and the scalar phase ϕ (background color) are shown in the x = y plane and the red domain denotes the almost pure ferromagnetic phase (| F | > 0.95).The field of view is 8 × 8 a 2r .Random noise with a magnitude of 10% of the order parameter amplitude is added to the initial state and the monopole is placed randomly in the origin-centered ball of radius a r by the creation ramp.See TableIfor parameter values.

FIG. 9 .
FIG. 9. Spontaneous flip in the winding direction during the time interval from T = 81 ms to 84 ms, as indicated in each panel.The scalar phase ϕ on the x = 0 plane is represented by the background color and ferromagnetic phase by the red domain.The field of view in all of the panels is 6 × 6 a 2 r .See Table I for parameter values.

FIG. 10 .
FIG. 10.Monopole as a result of Alice ring contraction at T = 82.4ms.The director d and the scalar phase ϕ on the x = 0 plane are represented by the black arrows and the background color, respectively.The field of view is 4 × 4 a 2 r .SeeTable I for parameter values.

023104- 8 FIG. 12 .
FIG.12.Annihilation of the Alice ring and anti-Alice ring.The ferromagnetic-phase cores (red domains), the director d (black arrows), and the scalar phase ϕ (background color) are shown for evolution times from T = 6 ms to T = 14 ms.In all of the panels, the field of view is 14 × 10 a 2 r .The initial state is perturbed with noise, 10% of the order parameter amplitude, and the initial monopole is placed near the top boundary of the cloud, at z = 2.42 × a r .See TableIfor parameter values.

FIG. 13 .
FIG. 13.Expectation value of the local average spin magnitude as a function of the evolution time T corresponding to Fig. 12.The initial state is perturbed with noise, 10% of the order parameter amplitude, and the initial monopole is placed near the top boundary of the cloud, at z = 2.42 × a r .The dashed red line marks the time instant when the Alice ring and anti-Alice ring annihilate.See Table I for parameter values.

TABLE I .
Physical parameters used in the numerical simulations during the Alice ring evolution.