Collective magnetism in arrays of spinor Bose–Einstein condensates

We elucidate dipolar magnetic interaction effects in spinor condensates at a low-background field. In particular, we show that arrays of ferromagnetic spinor Bose–Einstein condensates show a rich phase structure that is strongly influenced by the competition of shape anisotropy and higher-order magnetostatic contributions both depending on the condensate's form.


Introduction
Dipolar magnetism has triggered curiosity and applications from the compass to data storage. However, on the microscopic scale, the long-range dipolar interaction is in general overwhelmed by Heisenberg exchange interaction. Only a few frustrated condensed matter systems, such as spin-ice [1,2], reveal a glimpse of the wealth of phenomena dipolar systems offer. These perspectives have triggered significant efforts in artificially created dipolar nano-materials [3]. Spinor Bose-Einstein condensates (BEC) promise a new route into this area of physics [4][5][6] encouraged by recent experiments on quantum quench dynamics [7,8]. The problem is that the interparticle dipolar interaction has to dominate over interactions with an external field, which requires an ultralow magnetic field. Meystre and co-workers [9,10] suggested a solution in the form of a lattice of mini-BECs with a collectively enhanced magnetic moment of the condensates at each site. However, the dipolar physics in realistic experimental systems goes far beyond what could be anticipated within the point-dipole assumption made in these pioneering investigations.
In this paper, we show that competition of two additional geometry-dependent energy scales adds to the richness of the system. The first originates from the higher-order multipolar contributions to the magnetic moment of a homogeneously magnetized spheroid [11,12]. The higher-order terms have shorter range than the dipolar one and decay with the distance as 1/| r | 2l+1 , where l is the rank of the expansion. Nevertheless, these contributions are not negligible if the distance between the interacting bodies is comparable with their largest dimension [11]. The second is the magnetostatic self-energy. The separation of magnetic poles in each mini-BEC causes a self-induced demagnetizing field H d ↑↓ M leading to the so-called shape anisotropy E d causing a geometry-dependent preference of a certain magnetization direction.
We analytically calculate the self-energy and multipole moments of individual minicondensates and introduce the obtained data as input parameters into Monte-Carlo (MC) simulations. As a result of the numerical modeling, we derive thermodynamically stable magnetic configurations. We find that the shape anisotropy of a ferromagnetic mini-condensate dominates over the inter-condensate interactions if the long axis 2c of a spheroid is > √ 2 times larger than its short axis 2a. For that reason, the ground states of realistic BEC arrays differ from those predicted in [10,14] and show non-trivial structures, such as spin spirals. Furthermore, we predict the ground state of a two-dimensional (2D) square array of spin-domain condensates [13].

Magnetostatic interactions in arrays of ferromagnetic mini-spinor condensates
We first discuss the energy scales and their geometry dependence. In analogy to [10,14], we consider a one-dimensional (1D) or 2D array of ferromagnetic mini-spinor condensates at negligible external magnetic field. The term ferromagnetic refers to the local contact interaction between atoms mediated by spin-dependent s-wave collisions [15,16], which give rise to a mean-field energy minimum for aligned spins. An example is a 87 Rb spinor condensate in the F = 1 hyperfine state, where typical peak densities of 5 × 10 14 cm −3 give rise to a ferromagnetic interaction [17][18][19] of the order of h × 20 Hz per atom with Planck's constant h. This interaction sets the upper limit for dynamics and poses constraints on the allowable external magnetic field as well as the maximum multipolar interaction between sites. For example, in 87 Rb spinor condensates, a magnetic field of just 1 µG corresponds to a Zeeman energy of h × 0.7 Hz per atom; hence, any relevant dipolar interaction energies would have to be of the order of a few Hz per atom.
To estimate the inter-site dipolar interaction E dd , we approximate each condensate via a single macroscopic dipole moment M = g F m F µ B N located at its center, where N is the number of atoms in the condensate. For a typical peak density of 5 × 10 14 cm −3 , we get E dd = h × 8 Hz × R z R 2 ⊥ /a 3 per atom, where R z and R ⊥ are the Thomas-Fermi radii of the condensate and a is the lattice constant. The calculated energy corresponds to the interaction between two sites. As each site in a lattice has η next neighbors, E dd per atom in this case is at least η times larger. In a 2D lattice of prolate condensates, the macroscopic dipole approximation can be expected to remain valid up to R z ≈ a. With four next neighbors η = 4, to obtain an interaction energy of the order of h × 1 Hz per atom, it is necessary to have R ⊥ a/5. Similarly, for a 1D lattice of oblate condensates, R ⊥ ≈ a and R z a/16.
The above argument demonstrates that a situation where inter-site dipolar interactions are relevant is likely to involve condensates with large aspect ratios whose longer dimension is comparable with the lattice spacing. This geometry causes higher-order multipolar inter-well interaction as well as the magnetostatic self-interaction in each well to play an important role in the system.

Magnetostatic self-energy of a mini-spinor condensate
The magnetostatic self-energy density E d of a homogeneously magnetized body is described by the so-called demagnetizing factor N : where M s is the saturation magnetization and µ 0 M 2 s the magnetostatic energy of an infinite continuous magnet. The demagnetizing factor and, hence, E d is smallest when all magnetic moments compensate each other and the total magnetic charge is zero. E d increases whenever uncompensated magnetic poles (charges) are created in a material or at the boundaries. The preference of a certain magnetization orientation in order to decrease E d is known as the shape anisotropy E d . It corresponds to the difference between the self-energy of the most unfavorable E max d and most favorable E min d states. The demagnetizing factors for spheroids are well known [20]. Generally, to minimize the uncompensated magnetic charges, the magnetization of an individual, freestanding ellipsoid is collinear with its longest axis.
The shape anisotropy E d of a condensate can be approximated by treating it as a rotationally symmetric ellipsoid of volume V with homogeneous magnetization m = g F m F µ B n and n = N /V . For 87 Rb in the F = 1 state, this results in a shape anisotropy of Hz per atom, where D z and D ⊥ are the demagnetizing factors for axial and perpendicular magnetization, respectively. For a prolate ellipsoid with aspect ratio 5:1, D z − D ⊥ ≈ −0.31, for an oblate ellipsoid with aspect ratio 1: The multipole moments are the coefficients of a series expansion of a potential (typically of the form 1/ r − r ) due to a charge distribution. The main purpose of a multipole expansion is the approximation of the potential outside a finite charged body, which then can be used to calculate the stray field or the interaction energy with a second body exposed to this field. The stray field of a homogeneously magnetized body is also due to the surface or volume magnetic poles. In spherical coordinates, the moments have the form where R is the regular normalized spherical harmonic function. The analytical solutions for multipole moments Q lm of general ellipsoids with semi-axis (a, b and c) and arbitrary magnetization direction have been derived in [12]. According to these calculations, a singlespin-domain ellipsoid possesses odd multipolar contributions (l = 1, 3, 5, . . .), while a two spindomain ellipsoid even (l = 2, 4, 6, . . .) multipolar contributions only. The set of involved ranks of the multipole expansion strongly depends on the aspect ratio c/a of the spheroidal drop and on the orientation of its magnetization. A perfect sphere possesses a dipolar moment only, while a strongly prolate or an oblate spheroid has non-negligible higher-order contributions up to the rank l = 7. The higher-order contributions to the dipolar moments are visualized in figures 1(a) and (b) by means of the real-spherical harmonics corresponding to equipotential surfaces and the color reflects the sign of the potential. A dipole in this representation consists of two potential lobes of different signs only. An ellipsoidal drop, magnetized along its longer axis, possesses octopolar (l = 3) contributions appearing like additional 'collars'. If the magnetization is rotated in arbitrary direction, the shape of equipotential surfaces deforms as the non-negligible azimuthal contributions m appear. Hence, the interaction energy between two BECs, E i j , which can be calculated using equation (37) of [11] for non-intersecting charge distributions, depends on the shape of a BEC and magnetization orientation. Generally, the lowest interaction energy between one-domain condensates corresponds to the 'head-to-tail' magnetization of the neighboring spheroids magnetized along their shorter axis.

Competition between the magnetostatic self-energy and the magnetostatic interactions
As the discussion above shows, the shape anisotropy and the magnetostatic interaction between individual BECs require different magnetization configurations. To compare E d and E i j , one needs a reliable normalization. For that purpose, we note that the shape anisotropy of an oblate spheroid of volume V with the aspect ratio c/a → 0 is E d = µ 0 V M 2 s /2. This value is identical to the magnetostatic energy of an infinite 2D array of touching spheres (see figure 2(a)). Indeed, as a perfect sphere corresponds to a dipole posed into its center, the maximal energy for such an array is [21]. In the following, all energy densities will be expressed in units of µ 0 M 2 s /2. In figure 1(c), the interaction energy between a pair of BECs is compared with their shape anisotropy for different aspect ratios of the spheroidal drop. For the sake of comparison, the distance between two spheroids was fixed to r i j = d, while semi-axis a r i j , b ⊥ r i j (in-plane) and c Oz could be varied between d/2 and zero. The situation c = a = b = d/2 corresponds to two touching spheres, c = b > a to two oblate, while a = b < c to two prolate spheroids (see figure 2(a)). The described geometry corresponds to a typical optical 1D chain of oblate or a 2D array of prolate BECs (see figure 2(b)).
The shape anisotropy of a perfect sphere is zero, because E d is independent of the magnetization orientation. The interaction energy between two such touching spheres placed on the x-axis and magnetized in z-direction, in contrast, is maximal and equal to the interaction between point dipoles. The situation reverses for the limiting cases of pancake-or cigar-like spheroids (see figures 1(c) and (d) respectively). Here, E d becomes maximal, while E dd vanishes. As pointed out above, all non-spherical BECs possess higher-order magnetostatic contributions. Consideration of the higher-order terms changes the strength of E dd as a function of the aspect ratio (dashed lines in figures 1(c) and (d)). This change seems to be small in figure 1. However, the quantitative analysis demonstrates that the octopolar contribution enlarges the interaction energy by up to 40% with respect to that of pure dipole-dipole coupling for an oblate spheroid and diminishes it by up to 20% for a prolate one. The functions E d (a/d) and E Q1m+Q3m dd (a/d) cross at a ≈ 0.26d for oblate and a ≈ 0.32d for prolate spheroids. It means that for all aspect ratios larger than c/a(oblate) ≈ 1.9 and c/a(prolate) ≈ 1.6, the shape  figure 1(f), magnetization can take any orientation in the zx-plane parallel to the bond, while in figure 1(e) the zy-plane perpendicular to the bond is analyzed. The angles β 1 and β 2 refer to the polar angles of magnetization of the first and the second spinors, respectively. The absolute minimum belongs to the configuration shown in the inset to figure 1(e), which coincides with the minimum of the self-energy. However, for a 5c, the local minimum of figure 1(f) has comparable energy. The insets show corresponding magnetization configurations. Similarly to the previous case, the higher-order contributions do not change the position of energy minima, but make them much deeper.

Monte-Carlo simulations and discussion
To analyze the consequences of multipolar contributions and the shape anisotropy on ground states of BEC arrays, MC calculations have been performed. The Hamiltonian of magnetostatic interactions in spherical coordinates reads where Q A l A m A and Q B l B m B are the moments of multipoles A and B expressed in spherical harmonics [11] and T l A l B m A m B ( R AB ) is the geometric interaction tensor depending on the interparticle distance vector R AB between multipoles on sites A and B The distance dependence is given by the complex conjugate of the irregular normalized spherical harmonic function and K x is the easy-plane shape anisotropy due to the magnetostatic self-energy (it is a standard procedure to replace the intra-body self-energy E d by an additional anisotropy term [22]). The technical details of the MC technique are given in [11]. At each MC step the magnetization of a randomly chosen mini-condensate has been rotated by a random angle. The multipole moment of the BEC with rotated magnetization has been recalculated on the basis of formulae given in [12] and the new energy calculated using equation (2). The new configuration has been accepted or rejected in the framework of slow tempered annealing with a single spin update based on the Metropolis algorithm [11]. First, we discuss the thermodynamically stable states of linear chains of pancake-like BECs. In contrast to [9,14], the Hamiltonian we use includes an additional anisotropy term +K S 2 x corresponding to the shape anisotropy described in figure 1 and higher-order multipoles. For K = 0 (spherical BECs), an absolute energy minimum is a ferromagnetic, head-to-tail alignment of magnetization along the x-axis in agreement with [9]. For K → ∞, corresponding to c/a → ∞, a perfect, mutually antiparallel alignment of neighboring moments perpendicular to the chain has been found. For values of K E dd /ζ (with ζ -the Riemann zeta-function accounting for the dipolar sum in an infinite chain), however, the stable state of a finite chain is a modulated antiferromagnetic spiral depicted in figure 3(a). The easy-plane S z and S y components of magnetization show an antiferromagnetic modulation, while the magnetization along the chain axis S x demonstrates a vanishingly small signal of same periodicity. As follows from the magnetization distribution in the zy-plane (side view), the full 2π rotation is achieved on the scale of approximately 50 a for K ≈ 1.7E dd corresponding to the spheroids with the modest aspect ratio of 1:10. The shape anisotropy tries to keep the magnetization in the zy-plane, while the magnetostatic coupling favors the ferromagnetic alignment along the chain axis. The modulated spiral is a compromise solution allowing to minimize both energy contributions. The modulation period increases with increasing anisotropy and with increasing chain length and becomes infinitely large for K 5E dd .
The internal energy of such a spiral found in our simulations is negligibly larger than that of a perfect antiferromagnetic alignment E 0.01E dd . This result has also been confirmed by the analytical evaluation of the sum of equation (2) for the equilibrium distribution of magnetization S y (r ) = sin(δπr )sin(qπr ), where δπr is the angle made by the magnetization with respect to the zy plane. For a given aspect ratio of a mini-BEC, a specific set of δ and q exists for which the energy of the modulated spiral (q = 0) is close to that of the antiferromagnetic energy minimum (q = 0). The stability of the spiral state can be explained by the configurational entropy. While there are only two perfect antiferromagnetic configurations, at least 4πρ modulated states of identical energy exist (ρ is the density of states). The entropy ln(4π)ρk B T decreases the free energy favoring the non-trivial spiral configuration found in numerical simulations. The higher-order multipolar contributions strengthen the shape anisotropy and lead to an increase of the modulation period. The spiral is stable in the frame of the chain, but slowly rotates in the global coordinate frame because of the rotational symmetry. Therefore, the time-averaged magnetization is close to zero. The MC solution for the ground state of a 2D square lattice of pure dipoles (spherical BECs) is identical to the solution a = b of [14]. Inclusion of the shape anisotropy for spheroids with c > 2a transforms this pattern into the checkerboard 2 × 2 out-of-plane configuration for any multipole rank. The mostly non-trivial stable configuration has been found for BECs with spin domains. The corresponding configuration is visualized in figure 3(b). The minimum of the free energy is given by the configuration of figure 2(e) for individual condensates. The macroscopic configuration of the entire array consists of a checkerboard pattern of these microscopic spin-domains with two-in and two-out m F . Occasionally, the array of BECs acquire a domain wall dividing the sample into two phase shifted macroscopic domains. The magnetization of BECs inside of the wall rotates from a two-in to the two-out alignment ( figure 3(b)).
The results of calculations described above demonstrate that the multipolar magnetic interactions and the shape anisotropy strongly influence the ground states of optical BEC lattices and lead to non-trivial, stable configurations. This opens a large area for future investigations potentially including artificial spin-ice geometries with magnetic monopole excitations in coldatom systems.