Symmetry breaking in binary Bose-Einstein condensates in the presence of an inhomogeneous artificial gauge field

We study a two component Bose-Einstein condensate in the presence of an inhomogeneous artificial gauge field. In response to this field, the condensate forms a localised vortex lattice structure that leads to a non-trivial symmetry breaking in the phase separated regime. The underlying physical mechanism can be understood by considering the energy landscape and we present a simplified model that is capable of reproducing the main features of the phase separation transition. The intuition gained by numerically solving this simplified model is then corroborated using the analytical Thomas-Fermi model.


I. INTRODUCTION
Ultracold gases of neutral atoms have, in the past two decades, evolved into highly controllable systems that allow one to study and simulate numerous fundamental quantum mechanical effects [1][2][3]. One of the reasons for this is the large experimental toolbox for tuning almost all of the terms of their Hamiltonians using static or timedependent external fields. This includes using spin-orbit coupling [4] or artificial gauge fields [5,6] to affect the kinetic part, optical lattices [2], density-dependent gauge potential [7] or painted potentials [8] to adjust the external trapping terms, or Feshbach resonances to control the non-linear interaction terms [9], to name just a few. Additionally, systems with different symmetries can be created using multicomponent setups [10,11].
The simplest multi-component system is a binary Bose-Einstein condensate (BEC) made either from bosonic atoms in two different hyperfine states [12], two different isotopes [13], or two different elements [14][15][16][17]. These systems show intriguing physics related to interpenetrating superfluidity [18,19] and in particular possess a de-mixing phase transition [13,[20][21][22]. The latter is mostly determined by the interplay between the different interaction energies and, in free space, occurs when the square of the inter-component interaction strength exceeds the product of the two intra-component interaction strengths. It is worth noting that for certain atomic condensate settings these three interaction strengths can in principle be tuned independently [13,23]. In nonhomogeneous systems the point of the separation transition can be shifted, as the effects of the kinetic energy have to be taken into account [24]. Other terms that can be present in the Hamiltonian, e.g. accounting for Rabi coupling [25], spin-orbit coupling [26], or rotations [7,27], are known to have an influence on the phase separation threshold as well. * sahar.hejazi@oist.jp In this work we are interested in the phase separation process in a two-component system in the presence of angular momentum. However, contrary to previously considered situations [19,28], we will investigate systems where the rotational energy is not homogeneously distributed over the whole condensate. While such a situation can in principle be realised experimentally by locally creating vortices through phase imprinting [29,30], this technique usually leads to non-equilibrium situations as the condensate has to adjust its density to accommodate the imprinted phase distribution. To avoid excitations, such as phonon modes, which can have non-negligible effects on the phase-separation transition [31], we instead consider spatially inhomogeneous artificial gauge fields that only induce rotation in certain areas of the condensate [32][33][34].
Condensates in harmonic traps have been shown to respond to homogeneous rotation with the formation of triangular vortex lattices [35], whereas in different external potentials different geometrical arrangements of the vortices are possible [36][37][38][39]. Furthermore, condensates that encompass low-density regions, either due to local potential maxima or in the phase separated limit of a multicomponent system, can support the so called hidden or ghost vortices located in these regions [27,28,[40][41][42][43]. While the response to localised rotation through a gauge field has already been explored for single-component condensates [32][33][34], the effect on the phase separation transition in two-component systems has not yet been discussed.
In order to clearly isolate the effects of localised rotation, we consider systems with as many symmetries as possible: both condensates are made from atoms of the same species, both have the same number of particles and both have identical intra-component interaction strengths. We also restrict ourselves to a fundamental two-dimensional dynamics and a rotationally isotropic trapping geometry. Without rotation the separation transition in such a system leads to a straight phase boundary that cuts through the center of the trapping potential and whose direction is due to spontaneous symmetry breaking. In the presence of strong, global external rotation, this is no longer the case and the phase separation dynamics becomes highly complex and breaks all spatial symmetries by forming unordered serpentine vortex sheets [27].
To create a situation which lies in between the nonrotating and globally rotating settings we consider a gauge field that originates from an evanescent optical field above the surface of a prism, close to which a twocomponent condensate is trapped. The short-range exponential decay of the evanescent field in the direction perpendicular to the prism surface then results in an artificial magnetic field with a pronounced maximum at some distance from the surface. While in the miscible regime this produces a localised vortex distribution in the direction parallel to the surface that is in principal consistent with the symmetric splitting of the two components, we show that the interplay between the kinetic and the interaction energy leads to additional symmetry breaking that is not purely determined by minimising the length of the phase boundary.
This manuscript is organised as follows. In Section II, we describe our model for a two-component Bose-Einstein condensate in the presence of a nonhomogeneous artificial gauge field originating from the evanescent field created at the surface of a dielectric prism. In Section III we show how this artificial gauge field affects the miscible and immiscible regimes and in Section IV we study, through a simplified model, the physical mechanisms behind the symmetry breaking observed in the immiscible regime. The numerical results obtained are then corroborated using an analytical Thomas-Fermi model in Section V, and finally we conclude in Section VI.

II. MODEL
We consider a two-component Bose-Einstein condensate of neutral alkali atoms that is tightly confined in one spatial direction, such that it can be effectively treated using a two-dimensional description. In particular, we choose a harmonic trapping potential of the form , where the frequency ω is the same in both directions, so that the trap is symmetric in the x-z plane. Furthermore we assume that the atoms in both components have the same mass M , which can be achieved by trapping a single species and condensing the atoms in two different internal states. The entire system is located just above the surface of a dielectric prism with refractive index n, so that the atoms can interact with the evanescent field, see Fig. 1.
Within the mean field approach, the two-component BEC can be described by a set of coupled Gross- Pitaevskii equations (GPEs) of the form [44] i where g ll = 4N π 2 a l M and g lm = 4N π 2 a lm M are the intraand inter-component scattering strengths, respectively, with l, m = {1, 2} and l = m. As usual, a l is the s-wave scattering length between atoms of the same component and a lm for atoms of different components. The condensate wave-function is described by φ l and is normalized as |φ l | 2 dx dz = 1. The vector field A represents the gauge potential, which stems from the evanescent field emanating from the prism surface.
To describe the gauge field we consider a laser field with a wave-vector k and frequency ω L , chosen to be close to the resonance of the atomic transition. This field propagates inside the prism at an angle θ with respect to its surface. When this angle is larger than the critical angle, θ 0 = arcsin( 1 n ), the beam undergoes total internal reflection and an evanescent field is created at the surface of the prism. The electric field, E(x, z, t), propagates in x-z plane with an amplitude E 0 and decays from the surface in the positive z direction with a penetration depth d = (k 0 n 2 sin 2 θ − 1) −1 . It takes the form where t TE (θ) = 2n cos θ n cos θ + i n 2 sin 2 θ − 1 −1 corresponds to the transmission coefficient, and the running phase is given by ϕ(x) = xk 0 n sin θ [33]. with the same colour coding, as plotted in [33].
The interaction between the evanescent field and the atoms in the condensate occurs via dipole coupling, d · E(x, z), where d is the electric dipole moment of the atoms. Without loss of generality we assume it to be the same for both components. In the rotating wave approximation this then leads to a dressed state of the form [33] |χ where Φ(x, z) = arctan |κ(x,z)| ∆ , κ(x, z) = d · E(x, z)/ and ∆ = ω L − ω A is the detuning of the laser light from the atomic resonance frequency, ω A , which we assume again to be the same for both components. Assuming that the atoms move slowly enough to adiabatically follow this spatially inhomogeneous eigenstate, they pick up a geometrical Berry phase which can be written as the appearance of a vector potential A = i χ|∇χ , which has the explicit form An artificial magnetic field can then be calculated from the vector potential via B = ∇ × A as [33] B |∆| . From Eqs. (4) and (5) one can directly see that, since the evanescent field decays with increasing distance above the surface of the prism, the gauge field and the B-field will have to be inhomogeneous as well. For given sets of laser parameters this is shown in Fig. 2(a) for the A-field and in Fig. 2(b) for the B-field. In particular one can see from these plots that the artificial magnetic field has a maximum at finite distance away from the surface. The position of this maximum strongly depends on the value of s, while a change in the angle of the incident beam mostly affects the amplitude of the magnetic field. The atoms in a condensate trapped within the evanescent field will therefore experience effects corresponding to the presence of a spatially inhomogeneous B field [32,33].

III. EFFECTS OF THE INHOMOGENEOUS MAGNETIC FIELD
To explore the effects of the inhomogeneous artificial gauge potential, we numerically solve the two coupled Gross-Pitaevskii Equations (1) by using a standard FFT/split-operator method [45]. From here onward we work in harmonic oscillator units, that is x → x/a 0 , z → z/a 0 and t → tω, with a 0 = /M ω. We choose equal intra-component coupling coefficient, g 11 = g 22 = g and inter-component coupling is given by g 12 = g 21 ≡ αg. For stability reasons we only consider repulsive interactions. Thus, the condensate is in the miscible regime for 0 < α 1, and in the phase separated regime for α 1.
In Fig. 3 we show examples of ground state density profiles within the miscible and immiscible regimes. One can immediately note that the vortices only appear in a localised area, which corresponds to the region where the B-field is largest [32,33]. In the miscible regime and for the parameters chosen in Fig. 3 they form a single line along the maximum of the B-field, with each component having an offset with respect to the other such that they minimise the interaction energy. However, for less localised B-fields they can also arrange in a localised triangular lattices that converges to the full Abrikosov geometry for global fields [33].
While in the immiscible regime the condensate components separate, as expected, it is immediately clear from Fig. 3(b) that the separation can not be driven by the minimisation of the interaction energy alone. Naively one could expect that the separation boundary would be a straight line along the z-direction at x = 0, which would lead to minimising the interaction energy and the kinetic energy stemming from the boundary, while ensuring that both condensates have the same amount of vortices and energy. However, this would not necessarily minimise the overall energy of the system, as additional kinetic energy is associated with the vortices. In fact, Figure 3(b) shows that only one of the components carries visible vortices and that, even though the external parameters are the same for both situations, the number of vortices is not the same in the miscible and the immiscible regimes. This clearly indicates that some of the vorticity in the system is hidden in so-called ghost vortices, which are located in the low density areas at the phase boundary, so that the large rotational energy required to rotate high densities is avoided. To understand the interplay between the interaction and rotational energy in more detail, we will in the following explore a toy model of the inhomogeneous B-field that captures all relevant processes.

IV. TOY MODEL
Since the main characteristic of the inhomogeneous magnetic field is the existence of a localised maximum (see Fig. 2(b)), we will in the following consider the limit where the B-field is tightly localized in space. This can be achieved by assuming a step-like gauge-potential given by where Θ(z) is the Heaviside step function, z 0 is the shift of the Heaviside function in the z direction, and A 0 is the strength of the artificial gauge potential (see Fig. 1(b)). The field is constant in the x-direction. This form of the gauge potential catches the physical parameters that are related to the evanescent electric field created at the surface of the prism in a physically realistic and clean way: A 0 accounts for all the experimental parameters that characterize the strength of the realistic gauge potential (see Eq. 4) and z 0 accounts for the shift in real space due to the s parameter. In all our simulations below we chose A 0 such that its maximum value has the same order of magnitude as the realistic model. Typical ground state density distributions in the miscible and phase-separated regimes of the two-component system are shown in Fig. 4. In the miscible regime, α < 1, the localised B-field leads to a single-line of vortices, see Fig. 4(a). Due to the repulsive interactions between components, vortices within each component arrange themselves with an offset with respect to their counterparts, effectively filling the low density vortex cores of the other component. This is very similar to the realistic setting considered above, see Fig. 3(a).
For the phase separation regime, α 1, we show in Fig. 4(b-d) the density distribution for three different values of z 0 . In panel (c) the B-field is located exactly at the center of the BEC (z 0 = 0), and one can see that this leads to a separation of the two components into two clouds with essentially mirror-symmetric density profiles. The phase boundary is exactly along the line of the finite B field and corresponds to the shortest line possible. No vortices are visible and all vorticity is carried by ghostvortices located in the low density area between the two components [27,28,38]. This solution clearly minimizes the interaction and the kinetic energy of the system and is reminiscent of the standard phase separation in twocomponent systems without vorticity. However, the direction of symmetry breaking is now determined by the B-field and not chosen spontaneously.
For z 0 = 0 this simple picture breaks down and the additional kinetic energy in the system plays a crucial role in how the phase separation occurs. In Fig. 4(b) and (d) we show the density profiles for α = 1.5 and the B-field located at z 0 = −2.5 and z 0 = 2.5, respectively. One can immediately see that the two components separate in a non-symmetric way, which strongly depends on the position of the B-field, and that one component still possesses vortices, while the other does not. In particular one can see that the phase boundary is only partly along the line of the B field, before turning to be more aligned along the x-direction. The part along the z-direction increases in length with increasing z 0 , becomes the length of the whole condensate at z 0 = 0, and then decreases again almost symmetrically for z 0 > 0. While this symmetry breaking behaviour seems unusual at first sight, it can be intuitively understood by realising that the system is still trying to reduce the rotational energy by creating ghost vortices in the phase boundary region. Yet, when the B-field line does not cross the condensate symmetrically, separating the components fully along this line would lead to one component having a significantly smaller area available compared to the other. As the interaction energy is non-linear, this would lead to a significant increase in the overall energy, which is unfavourable. Thus, the system uses the B-field line partly to minimise the rotational energies, but then minimzes the interaction energies by departing from it. The phase separation is therefore a careful balance between the minimisation of the interaction and the rotational energies. It is important to realise that the situations for values of z 0 that are symmetric around zero are not fully identical (see Figs. 4(b) and (c)), as the A-field breaks the system's symmetry along z. Again, this is qualitatively the same behaviour that is also found in the realistic model shown in Fig. 3(b).
The intuition for the immiscible regime developed above can be supported by looking at the kinetic and interaction energy of each component given by One can see from Fig. 5 that the kinetic energy of the vortex-carrying component grows initially much faster as the B-field line moves through the condensate, compared to that of the second component. However, with increasing values of z 0 , the length of the phase separation border along the B-field line grows, leading to more and more vortices turning into ghost vortices. The kinetic energy therefore decreases again, until the same value is reached for both components when z 0 = 0, i.e. when all vortices have become ghost vortices. The same process then repeats in the second component, which starts carrying the vortices once z 0 > 0. The graph of the interaction energy as a function of z 0 (see Fig. 5(b)) shows that the component carrying the vortices has generally a lower non-linear energy than the one that carries no angular momentum. This is due to the additional centrifugal forces in the vortex carrying component, which allow the system to achieve a lower density. Again, these variations in energy go to zero when all vortices have been turned into ghost vortices at z 0 = 0 and the role of the two components flips subsequently. It is worth noting that when the B-field passes z 0 = 0 a jump in the interaction energy can be seen as there exist a sudden point when the last visible vortex has been turned into a ghost vortex. Again, the asymmetries present in the kinetic and interaction energies are due to the effect of the A-field, which increases the total energy of the system as it increasingly envelopes the entirety of the two component BEC.

V. SYMMETRY BREAKING
While the part of the phase separation line along the B-field line is set by external parameters, the remaining question is about the position of the turning point and the direction of the break away from it. Intuitively it should be as short as possible, which for rotationally isotropic geometries should lead to a break at a right angle. One can see from Fig. 4(b)-(d) that this is approximately the case and below we confirm this intuition by determining the break-off point (x 0 , z 0 ) by energy minimization using the analytical Thomas-Fermi (TF) wavefunction, φ TF l , obtained from solving where µ TF l is the chemical potential of each component and we are again using harmonic oscillator units. This approximation is valid when the kinetic energy terms of the Hamiltonian can be neglected as they are much smaller than the non-linear ones (see Fig. 4).
However, as the kinetic energy clearly plays an important role in the phase separation, we take its effect into account by fixing the phase separation line along the maximum of the magnetic field in the x-direction up to a value of x 0 . We then approximate the rest of the phase separation border by a straight line along the z-direction, so that both parts have a sharp π/2 angle between them (see inset of Fig. 6). These conditions are encoded in the limits of integration for all integrals which depend on the position of the vertical part of the phase boundary, i.e. x 0 . For simplicity we also fix the TF radius to R TF = 2 g/π in all limits of integration, so that no extra functional dependence on x 0 or z 0 appears. This allows us to avoid coupled transcendental equations for the chemical potentials of both components.
Within this model, we then minimize the total energy of the two-component system as a function of the position of the break away from the phase separation along the B-field line at x 0 (see [46] for details). The results from this analytical approach are shown in Fig. 6(a) and one can immediately see that they very closely match the ones found from numerically solving the full twocomponent GPEs. This indicates that the straight line boundary along the x-direction connected at a right angle to the first part of the boundary along z 0 provides the lowest energy solution for the system to phase separate and forces it to spontaneously break the symmetry. We also show in Fig. 6(b) that the length of the phase boundary exceeds the diameter of the condensate whenever the position of the B-field breaks the symmetry of the system. Finally, it is worth noting that while the assumption of a right angle connection between the two parts of the phase boundary is a good assumption in the toy model of the gauge field, it only approximately holds in the realistic model discussed in Sec. III. This is due to the B-field being spread out over a larger range in  6. (a) Turning points (x0, z0) of the asymmetric phase separation (see Fig. 4 for density plots) comparing the analytical Thomas-Fermi solution given by the solid-black line with the GPE numerical simulations, green dots (see [46] for details) for α = 1.5. This plot shows the first quadrant of the 2D representation of the BEC, with (0, 0) corresponding to the center of the harmonic trap. The dashed-red line identifies the area delimited by the Thomas-Fermi radius found at α = 0. (b) Length of the phase boundary, L b , within the TF approximation as a function of the position of the maxima of the magnetic field z0. The dashed-red line indicates the shortest path for the phase boundary, i.e., 2RT F . the z-direction and the vortices being discreet along the B-field line.

VI. CONCLUSIONS
In this work we have investigated the ground state of a two-component Bose-Einstein condensate in the presence of a inhomogeneous artificial gauge potential. This situation appears when the gauge field is created by an evanescent field at the surface of a dielectric material, close to which the BEC is trapped, and it is therefore experimentally realistic. While in the miscible regime the systems responds to the angular momentum imposed by the gauge field in an expected manner, in the phase separated regime a careful balance between the need to minimise the interaction and the kinetic energy leads to a phase separation that spontaneously breaks the symmetry in unusual ways. In particular, the phase separation border is no longer just a straight line that crosses the system symmetrically as would be the case in the absence of the gauge field. Using a toy model we have carefully explored the mechanism behind this symmetry breaking two-component state and clearly described the importance of the kinetic energy in the phase separation process. To confirm our numerical results, we have also presented an energy minimization calculation using the analytical TF solution that allows to determine the position where the phase boundary turns away from being a straight line.
Using inhomogeneous gauge potentials to induce rotation locally into condensates holds the potential to be a valuable way to engineer and study interesting super-fluid dynamics. These can range from the above study on phase separation in multi-component condensates to creating well-defined initial states to study quantum turbulence [47]. The fact that such systems are experimentally possible using today's technology makes this an exciting and promising direction of research. This project was supported by Okinawa Institute of Science and Technology Graduate University. JP also acknowledges the JSPS KAKENHI Grant Number 20K14417.