Tuning the structural and dynamical properties of a dipolar Bose-Einstein condensate: Ripples and instability islands

It is now well established that the stability of aligned dipolar Bose gases can be tuned by varying the aspect ratio of the external harmonic confinement. This paper extends this idea and demonstrates that a Gaussian barrier along the strong confinement direction can be employed to tune both the structural properties and the dynamical stability of an oblate dipolar Bose gas aligned along the strong confinement direction. In particular, our theoretical mean-field analysis predicts the existence of instability islands immersed in otherwise stable regions of the phase diagram. Dynamical studies indicate that these instability islands, which can be probed experimentally with present-day technology, are associated with the going soft of a Bogoliubov--de Gennes excitation frequency with radial breathing mode character. Furthermore, we find dynamically stable ground state densities with ripple-like oscillations along the radial direction. These structured ground states exist in the vicinity of a dynamical radial roton-like instability.


I. INTRODUCTION
Dipolar Bose Einstein condensates (BECs) such as magnetic Cr condensates are characterized by angledependent long-range interactions [1,2,3]. Usually, the magnetic dipole-dipole interactions compete with the isotropic short-range s-wave contact interactions, which dictate the behaviors of alkali atom BECs [4] as well as of two-component Fermi gases such as 40 K and 6 Li [5,6]. In Cr, the s-wave scattering length a s can be tuned to zero through the application of an external field in the vicinity of a Fano-Feshbach resonance, thus allowing for the realization of a pure dipolar Bose gas [7,8,9]. Although seemingly simple, dipolar Bose gases in which the dipoles are aligned along a particular laboratory axis have been shown to exhibit a variety of intriguing and unique features such as so-called roton instabilities [10,11] and structured density profiles [12,13,14]. Importantly, although the strength of the magnetic dipole-dipole interactions is fixed by the electronic structure of the atoms, it can be tuned effectively by varying the aspect ratio of the external confining potential [15,16]. This effective tunability of the dipolar interactions underlies a series of studies of dipolar gases and can be readily understood intuitively. Assuming the dipoles are aligned along the z-axis, an effectively one-dimensional confinement that predominantly allows for motion along the z-axis leads to an attraction between the dipoles and thus collapse. An effectively two-dimensional confinement that predominantly allows for motion in the xy-plane, in contrast, leads to a repulsion between the dipoles and thus stabilization.
This paper considers a pure dipolar Bose gas, aligned along the z-axis, in an oblate trapping geometry with cylindrical symmetry characterized by the aspect ratio λ. In addition, a repulsive Gaussian barrier with fixed width b and variable height A is added along the z direction (centered at z = 0). When the chemical potential µ is much larger than the barrier height A, the Gaussian potential serves as a small perturbation. For µ/A ≪ 1, in contrast, the fairly narrow Gaussian barrier significantly modifies the system behavior and approximately splits the dipolar condensate into two "sheets" or "pancakes". We determine the phase and stability diagram of the dipolar gas as a function of the mean-field strength D and the barrier height A for various aspect ratios λ.
Our key results are as follows: (i) We find an instability island immersed in otherwise stable regions of the D versus A phase diagram. Our analysis of the Bogoliubov-de Gennes excitation frequencies and corresponding eigenmodes indicates that these unstable islands arise due to a breathing mode-like instability. (ii) We find dynamically stable density profiles of Gaussian shape with superimposed ripple-like structure. For certain parameter combinations near the instability line, these ripple-like structures exist for vanishing as well as finite Gaussian barrier. As the mean-field strength increases, one of the Bogoliubov-de Gennes frequencies with radial roton type character becomes soft, inducing a dynamical instablility.
Experimentally, the instability island predicted by our calculations can be probed straightforwardly by following two different types of trajectories in the D versus A phase diagram: (i) The first experiment probes the instability by adiabatically increasing the barrier height from A 1 to A 1c or by adiabatically decreasing the barrier height from A 2 to A 2c for fixed aspect ratio λ, fixed mean-field strength D and fixed barrier width b (see trajectories 1 and 2 in Fig. 1). The experimental signature of the collapse would be the onset of significant atom loss at A 1c and A 2c , respectively. (ii) The second experiment by-passes the instability island by moving from A 1 to A 2 and back to A 1 by simultaneously varying A and D while keeping λ and b fixed (trajectories 3 and 4 in Fig. 1 trace out one of many possible paths). In the case of Cr, varying the mean-field strength D, which is determined by the number of particles N and the square of the dipole mo- The remainder of this paper is organized as follows. Section II introduces the system under study and the mean-field formalism. Section III presents our results and Sec. IV concludes.

II. MEAN-FIELD DESCRIPTION OF DIPOLAR BECS
Our description of the aligned dipolar Bose gas with vanishing s-wave interactions is based on the timedependent Gross-Pitaevskii equation [17,18,19] i where m denotes the mass of the dipoles, d the strength of the dipole moment, and ϑ the angle between the relative distance vector r − r ′ and the z-axis. Throughout, we employ cylindrical coordinates r = (ρ, ϕ, z). The external confining potential V ext (ρ, z), is characterized by the transverse and longitudinal angular trapping frequencies ω ρ and ω z , respectively, and the height A and width b of the repulsive Gaussian barrier. Throughout, we fix b, i.e., we set b = 0.2a z , where a z = /(mω z ). The angular trapping frequencies define the aspect ratio λ, Throughout, we consider oblate systems with λ = 6 − 9.
If the Gaussian barrier is sufficiently large, the system consists-to lowest order-of two staggered pancakes; in this case, we define an effective aspect ratio λ eff , which is obtained by analyzing the energy spectrum for large A and vanishing interactions, i.e., for D = 0 (see below). To quantify the strength of the mean-field term in Eq. (1), we define the dimensionless quantity D, where E ρ = ω ρ and a ρ = /(mω ρ ). Note that the quantity D defined in Eq. (4) differs from the definition used in our previous work [20] but agrees with that used in Ref. [12].
We consider stationary and dynamical properties of the dipolar Bose gas. Our stationary calculations determine the energetically lowest-lying gas-like state of the dipolar system. We write Ψ( r, t) = exp(−iµt/ )ψ 0 (ρ, ϕ, z), where ψ 0 denotes the stationary ground state solution and µ the corresponding chemical potential. Since the system is cylindrically symmetric, ψ 0 is a function of ρ and z only and not of ϕ. We solve the resulting twodimensional Schrödinger equation by propagating a given initial state in imaginary time. In selected cases, we additionally employ a basis set expansion type approach. We have checked that the results of these two distinctly different approaches, whose implementation follows that outlined in Ref. [12,20,21], agree within our numerical uncertainties.
In addition to the ground state properties, we determine the Bogoliubov-de Gennes eigenfrequencies and eigenmodes, which provide insights into the dynamical stability of the system [4,22]. To this end, we solve the Bogoliubov-de Gennes eigenequations, which are derived within linear response theory assuming that the ground state ψ 0 is perturbed by δψ( r, t), iteratively using the Arnoldi method [23]. Here, u( r) and v( r) are the so-called Bogoliubov-de Gennes functions and ω is the Bogoliubov-de Gennes eigenfrequency. Following the approach introduced in Ref. [23], we solve the Bogoliubov-de Gennes eigenequations for ω 2 and for the so-called Bogoliubov-de Gennes eigenmode f ( r), where f ( r) = u( r) + v( r). Because of the cylindrical symmetry of the system, the ϕ dependence can be separated off, i.e., f ( r) =f (ρ, z)h(ϕ) with h(ϕ) = exp(ikϕ)/ √ 2π. This implies that the Bogoliubov-de Gennes excitation frequencies ω can, for a given azimuthal quantum number k, be obtained by solving a two-dimensional eigenvalue problem. Typically, we consider excitations with up to k = 4. For the discussions in Sec. III, it is important to recall that a negative ω 2 , and thus a purely imaginary ω, signals a dynamically unstable ground state ψ 0 ( r). Correspondingly, an analysis of the Bogoliubovde Gennes eigenmodes f ( r) near the dynamical instability point provides insights into the decay mechanism; in particular, the correction to the time-dependent density is, to lowest order in the perturbation δψ, given by 2 cos(ωt)ψ 0 ( r)f ( r).
To elucidate the energy and length scales of V ext (ρ, z), we consider the non-interacting system (i.e., we set D = 0). In this case, the linear Schrödinger equation separates and the solutions are readily obtained. In the (ρ, ϕ)-plane, we have a two-dimensional harmonic oscillator with eigenenergies (2n ρ + |k| + 1)E ρ , where n ρ and k denote the principal and azimuthal quantum numbers, respectively (n ρ = 0, 1, · · · and k = 0, ±1, · · · ). In the z-direction, the energy spectrum depends on two parameters, the axial frequency ω z and the barrier height A (recall, the barrier width b is kept constant throughout this paper). Figure 2 shows the lowest six eigenenergies, obtained by diagonalizing the Hamiltonian matrix expressed in terms of harmonic oscillator functions, as a function of A for b = 0.2a z . For A = 0, the eigenenergies follow the harmonic oscillator pattern (n z + 1/2)E z , where n z = 0, 1, · · · and E z = ω z . For large A, in contrast, the eigenenergies appear in roughly equally spaced pairs [ Fig. 2(a) shows three of these pairs]. We use the energy spacing between the lowest two pairs to define an effective angular frequency ω z,eff that characterizes the left and the right well of the double well potential. Extrapolating to the large barrier height limit, we find ω z,eff ≈ 2ω z . Correspondingly, we define an effective harmonic oscillator length a z,eff through a z,eff = /(mω z,eff ), leading to a z,eff ≈ a z / √ 2. While this analysis is approximate-our definition of a z,eff , e.g., assumes a harmonic potential around the minimum of the left well and of the right well-, it provides insights into how the energy and length scales associated with the zdegree of freedom vary with increasing barrier height.
Lastly, we point out that the Gaussian barrier leads to an energetically low-lying tunneling splitting mode for comparatively large A and D = 0. For A = 16E z , e.g., we find that the density at z = 0 is about 25 times smaller than the peak density, suggesting that the left well and the right well are, to a good approximation, decoupled. Correspondingly, as shown in Fig. 2(b), the lowest frequency decreases from 1ω z for A = 0 to about 0.044ω z for A = 16E z . This implies that the tunneling splitting mode for D = 0 is, in the large A limit and for λ = 6 − 9, energetically lower-lying than the lowest mode along the ρ-direction which has a frequency of ω ρ . Our discussion of the excitation spectrum for vanishing D presented here serves as a guide to understanding the Bogoliubov-de Gennes excitation spectra for finite D (see Sec. III).

III. RESULTS
This section presents the energetics, the density profiles, and selected Bogoliubov-de Gennes excitation spectra and eigenmodes for trapped dipolar Bose gases. In particular, we detail calculations that predict an instability island (see Fig. 1) and density profiles with ripple-like oscillations.
To understand how the instability island for λ = 8 emerges, Figs. 3(a)-(c) show the D versus A phase diagram for fixed b and three different λ values, i.e., for λ = 6, 8 and 9. Compared to Fig. 1, the phase diagram for λ = 8 [see Fig. 3(b)] shows an extended (D, A) parameter regime. As in Fig. 1, mechanically stable regions of the phase diagram are separated by solid lines from mechanically unstable regions of the phase diagram. Dynamically stable regions are separated from dynamically unstable regions by symbols; filled circles, squares and diamonds indicate that the dynamical instability is triggered by a k = 0, 2 and 3 mode, respectively. The mechanical and dynamical instability points that separate the "upper" and "lower" regions of the phase diagram differ by up to about 30 %. Along these instabil-  Fig. 1. The dynamically stable regions are separated from the dynamically unstable regions by symbols: Filled circles, squares and diamonds indicate that the dynamical instability is triggered by a k = 0, 2 and 3 mode, respectively. Dashed lines separate stationary ground state densities with maximum at ρ = 0 (labeled "S0") from those with maximum at ρ > 0 (labeled "S>0"). In certain regions of the phase diagram, the ground state densities possess ripplelike oscillations (the corresponding phase boundaries are not indicated here; see text and Fig. 10).
ity lines, the decay is triggered by modes with vanishing or non-vanishing azimuthal quantum number k and has previously been discussed for A = 0 in Ref. [12] and for A = 12E z in Ref. [20]. In the vicinity of the instability island, which is centered at A ≈ 3.5E z and D ≈ 13 [see elliptically-shaped solid lines in Fig. 3(b) and Fig. 1], the k = 0 Bogoliubov-de Gennes mode becomes soft first; in this region, the dynamically unstable region is only slightly larger than the mechanically unstable region.
Inspection of Figs. 3(a)-(c) shows that the instability island for λ = 8 emerges from the break-up of the unstable region of the phase diagram for λ = 6 into two pieces [see Fig. 3(a)]. For λ = 7 (not shown), we find a somewhat larger mechanically unstable island than for as for λ = 8. Thus, as λ increases the instability island shrinks, and it is absent for λ = 9. Figure 3(b) shows that the instability island is surrounded by a dynamically stable region, suggesting that the experimental sequence discussed in the introduction in the context of Fig. 1 is feasible with present-day technology. Specifically, for a Cr condensate with vanishing s-wave scattering length, confined by a trap with ω ρ = 2π × 100 Hz and ω z = 2π × 800 Hz, point A 1 in Fig. 1 corresponds to N ≈ 7200 atoms, while the dynamical instability points for A = 3.5E z in Fig. 3(b) correspond to (from bottom to top) N ≈ 5900, 9300 and 14200. In determining the D versus A phase diagrams shown in Figs. 3(a)-(c), our ground state calculations employed a grid with widths 4 and 1E z in D and A, respectively. In a second set of calculations, a finer grid in D was employed to determine the mechanical instability point more accurately and to analyze the dynamical stability. We cannot rule out the existence of additional instability islands with smaller widths than our numerical resolution. Figure 4 shows the energy per particle E tot /N (solid line) and the chemical potential µ (dashed line) for λ = 8 and b = 0.2a z : Fig. 4(a) shows the energetics for A = 3.5E z (corresponding to 28E ρ ) as a function of D [i.e., along a vertical cut in Fig. 3(b)], while Fig. 4(b) shows the energetics for D = 13 as a function of A [i.e., along a horizontal cut in Fig. 3(b)]. For fixed D [see Fig. 4(b)], the energy and chemical potential increase monotonically till the mechanical instability point is reached at A ≈ 2.85E z . At this point, the chemical potential µ is about two times smaller than the barrier height. Beyond the mechanically unstable region, i.e., for A 4.4E z , both E tot /N and µ continue to increase monotonically. The energy per particle and the chemical potential behave similarly as a function of D for fixed A [see Fig. 4(a)]. Interestingly, the presence of the instability island manifests itself in the form of a "gap" in the energy and the chemical potential but does not notably change the slope of either of these quantities (i.e., E tot /N and µ appear to be changing smoothly across the instability island). We find a similar behavior for a subset of the Bogoliubov-de Gennes eigenmodes (see below).
In addition to the energetics, we analyze the structural properties of dipolar Bose gases in a double-well potential. Dashed lines in Figs. 3(a)-(c) separate the parameter region where the ground state density |ψ 0 | 2 has its maximum at ρ = 0 (labeled S 0 ) from that where |ψ 0 | 2 has its maximum at ρ > 0 (labeled S >0 ). The transition from the S 0 type to the S >0 type density is smooth and the dashed lines are based on an analysis of the integrated density n(ρ), n(ρ) = 2π |ψ 0 (ρ, z)| 2 dz. We denote the ρ value at which n(ρ) is maximal by ρ max and define the transition from S 0 to S >0 type densities, i.e., the dashed lines in Fig. 3, by the condition n(0)/n(ρ max ) = 0.98. Figures 5(a) and (b) show examplary integrated densities n(ρ) of types S 0 and S >0 , respectively. The latter density profile has, for a system with vanishing barrier, been previously termed "red blood cell" [12]; this name is motivated by the fact that the system's isodensity plot resembles the shape of a red blood cell. For all parame- ter combinations investigated, the dynamical instability is triggered by a k = 0 mode if the ground state density is of S 0 type and by a finite k mode if the ground state density is of S >0 type. In particular, Fig. 3(b) shows that the density profiles in the vicinity of the instability island have a simple Gaussian shape and that the dynamical instability in this regime is associated with a k = 0 mode. To shed light on the dynamics in the vicinity of the instability island, Fig. 6 shows the Bogoliubov-de Gennes excitation spectrum as a function of D for k = 0, A = 3.5E z , λ = 8 and b = 0.2a z . The lowest excitation frequency becomes purely imaginary at D ≈ 10, signaling the dynamical instability. The frequency "recovers" at somewhat larger D values (D ≈ 16). We find that the real parts of the eigenfrequencies with non-vanishing azimuthal quantum number k remain finite in the vicinity of the instability island, implying that the dynamics in this parameter space of the phase diagram is dominated by k = 0 modes. Our analysis of the Bogoliubov-de Gennes eigenmodes (see Figs. 7 and 8 as well as the discussion below) shows that the eigenmodesf corresponding to the Bogoliubov-de Gennes excitation frequencies shown by solid lines in Fig. 6 have nodal lines that are to a good approximation independent of z while those corresponding to the Bogoliubov-de Gennes excitation frequencies shown by dashed lines have, among other nodal lines, a nodal line given by z = 0. The eigenfrequencies shown by solid lines are affected by the instability (i.e., the slope of these frequencies changes near the dynamical instability), while those shown by dashed lines appear to change smoothly across the instability. This suggests that these two classes of eigenfrequencies are approximately decoupled. We find that the finite k frequencies also change as though they are uneffected by the instability island. For larger D values (D 24.5), however, modes with finite azimuthal quantum number become relevant in determining the system's dynamical stability for a fairly wide range of A values. In fact, as indicated in Fig. 3(b), the instability for comparatively large D values and A = 0 to A ≈ 10E z is triggered by k = 2 and 3 modes (i.e., for these A values the energetically lowest-lying k = 2 or 3 frequency becomes purely imaginary before the lowest k = 0 frequency does). Figure 9 shows the real part of the energetically lowestlying k = 0 Bogoliubov-de Gennes eigenfrequency as a function of D for λ = 8, b = 0.2a z and various A values. The frequency shown in Fig. 9(c) for A = 3.5E z has just been discussed in the context of Fig. 6. Inspection of Fig. 9 shows that the vanishing of the real part of the energetically lowest-lying breathing mode type frequency for A = 3.5E z is accompanied by a "dip" in the corresponding frequencies for smaller and larger A [see Figs. 9(b) and (d)] around D ≈ 10 and 15, respectively. For these A values, however, the mode "recovers" as D increases and the system is mechanically and dynamically stable up to comparatively large D values. We now discuss the Bogoliubov-de Gennes eigenmodes shown in Fig. 7 and in Fig. 8 in more detail. The k = 0 eigenmodes shown in Fig. 8 show one nodal line that is given by z = 0 and a second nodal line that depends on ρ and z. The former nodal line reflects the symmetry of the confining geometry and corresponds to oscillations between the left well and the right well, while the latter is of secondary importance since the amplitude off is small along this nodal line. We classify the eigenfrequencies associated with eigenmodes of the type shown in Fig. 8 as having tunneling splitting mode character. The k = 0 eigenmodes shown in Fig. 7, in contrast, are characterized by a single nodal line which depends to first order only on ρ and not on z. We classify the eigenfrequencies associated with eigenmodes of the type shown in Fig. 7 as having breathing mode character. An analysis of the eigenmodes corresponding to the energetically lowest-lying k = 0 frequency for other D values but the same A, λ and b shows that the nodal line is located at ρ = a ρ for D = 0 and bends slightly as D increases while remaining located around ρ ≈ a ρ [see Fig. 7(a)-(b)]. Past the instability, the nodal line off is again located at ρ ≈ a ρ [see Fig. 7(c)] and moves to larger ρ values with increasing D [see Fig. 7(d)]. Eventually, a new nodal line moves in from ρ = 0. Near the dynamical instability at D ≈ 36, the eigenmode is characterized by two approximately equally spaced nodal lines with approximately constant ρ (not shown; qualitatively, the eigenmode is similar to those shown in Fig. 11 but without the third nodal line). For large D [D ≈ 36, see Fig. 9(c)], the instability thus acquires some similarities with a radial roton instability (see below).
Lastly, we show that cylindrically symmetric systems with modest λ support, in addition to density profiles of simple Gaussian shape and of red blood cell type shape, density profiles with ripple-like oscillations. Figures 10(a)-(c) show examplary integrated density profiles n(ρ) with ripple-like oscillations for λ = 9 and three different (D, A) parameter combinations. In general, we find density profiles with ripple-like oscillations in a relatively small region near the dynamical instability points for λ = 9 but not for λ = 8 and 6 [the parameter combinations that support ripple-like density profiles fall into the S 0 region of the phase diagram shown in Fig. 3(c); they are not indicated explicitly in Fig. 3(c)]. Ripplelike structures have very recently been predicted to exist in dipolar systems with non-vanishing s-wave scattering length consisting of pancake-shaped stacks, where each pancake is characterized by a large aspect ratio λ, e.g., λ = 340 [24]. Ripple-like oscillations have also been predicted to exist for dipolar systems with vortices [14]. Here, we find ripple-like oscillations in the dynamically stable region of the phase diagram of pure dipolar gases confined by a cylindrically-symmetric trapping geometry with and without Gaussian barrier and with modest aspect ratios λ.
The integrated density profiles shown in Fig. 10(a)-(c) for A = 0 to 12E z possess density oscillations with characteristic length scale of the order of a ρ . Figures 11(a)-(c) show the corresponding eigenmodesf of the energetically lowest-lying k = 0 Bogoliubov-de Gennes eigenfrequency. The Gaussian barrier modifies the amplitude of the eigenmode near z = 0 but otherwise leaves the overall structure of the eigenmodes unaffected. All three eigenmodes possess roughly equally spaced nodal lines that are to a good approximation independent of z, and which can, roughly speaking, be characterized-just as the ripple-like oscillations-by ρ node ≈ na ρ , where n = 1, 2 and 3. Using λ = 9, it can be seen that ρ node is determined by ρ node ≈ πa ρ / √ λ or ρ node ≈ πa z [10,12]. The latter expressions reflect that the condensate develops three-dimensional character, i.e., that it "gets chopped up" into smaller pieces along the ρ-direction. The dynamical instability that arises when D is increased somewhat compared to the D values chosen in Figs. 10 and 11 [see also Fig. 3(c)] is thus identified as a radial roton-like instability [10,12]. Intuitively, one might expect that the length scale that characterizes the ripple-like oscillations and the nodal lines of the Bogoliubov-de Gennes eigenmodes would depend on the barrier height A or, equivalently, the effective aspect ratio λ eff . Figures 10 and 11, however, suggest that the system behavior is determined by λ instead. We note, though, that the change of λ eff as a function of A is fairly small, implying that further studies are needed to determine the relevant length scale unambiguously.

IV. SUMMARY
We have investigated oblate dipolar Bose gases with vanishing s-wave scattering length in a double well po- tential as a function of the barrier height A, the meanfield strength D and the aspect ratio λ. Our stationary and dynamical mean-field calculations add two new aspects to the already extensive list of intriguing behaviors of dipolar Bose gases: (i) We find an instability island immersed in a mechanically and dynamically stable region of the phase diagram. This instability island can, as has been outlined in Secs. I and III, be probed with presentday technology and arises due to the going soft of a radial breathing mode-like frequency. (ii) We find structured ground state densities with ripple-like oscillations in the dynamically stable region of the phase diagram for moderate aspect ratios and cylindrically-symmetric confining geometries. These ripple-like oscillations exist, for appropriately chosen parameter combinations, for vanishing and non-vanishing Gaussian barriers; an increase of the mean-field strength eventually induces a radial roton-like instability. In the future, it will be interesting to investigate if the ground state densities with ripple-like oscillations can be probed experimentally using time-of-flight expansion techniques. Theoretically, this question can be addressed by determining the real time dynamics after turning off the confining potential and by comparing the results with those for simple Gaussian density profiles. We emphasize that the emergence of the instability island and of density profiles with ripple-like oscillations is unique to dipolar gases-these features are a direct signature of the long-range, anisotropic dipole-diple interaction and absent for s-wave interacting Bose gases.
Current experiments are directed at loading molecular gases with appreciable electric dipole moment [25,26] into one-dimensional lattices. The premise of these experiments is to exploit the effectively two-dimensional geometry of each lattice site as a stabilization mechanism. In this context, several theoretical studies [24,27,28] have investigated stacks of effectively two-dimensional dipolar Bose gases. Common to these studies is the assumption of a very tight confinement along the z-axis, i.e., a large aspect ratio for each individual lattice site, and a comparatively large spacing between individual lattice sites. In those studies, roton and phonon instabilities emerge. The studies presented here can be viewed as an alternative approach to tackling the multi-site optical lattice system. Using a single mean-field equation and working in a regime where the extent of the cloud in one well is approximately equal to the separation between clouds, the full system dynamics is investigated for moderate aspect ratios. Future studies will investigate the system behavior as a function of the separation between sites, using a single mean-field equation and relatively small barrier heights as well as a set of coupled mean-field equations; these future studies will facilitate direct comparisons between the approach pursued here and Refs. [24,27,28]. Support by the NSF through grant PHY-0855332 is gratefully acknowledged.