Collective modes of vortex lattices in two-component Bose-Einstein condensates under synthetic gauge fields

We study collective modes of vortex lattices in two-component Bose-Einstein condensates subject to synthetic magnetic fields in mutually parallel or antiparallel directions. By means of the Bogoliubov theory with the lowest-Landau-level approximation, we numerically calculate the excitation spectra for a rich variety of vortex lattices that appear commonly for parallel and antiparallel synthetic fields. We find that in all of these cases, there appear two distinct modes with linear and quadratic dispersion relations at low energies, which exhibit anisotropy reflecting the symmetry of each lattice structure. Remarkably, the low-energy spectra for the two types of fields are found to be related to each other by simple rescaling when vortices in different components overlap owing to an intercomponent attraction. These results are consistent with an effective field theory analysis. However, the rescaling relations break down for interlaced vortex lattices appearing with an intercomponent repulsion, indicating a nontrivial effect of an intercomponent vortex displacement beyond the effective field theory. We also find that high-energy parts of the excitation bands exhibit line or point nodes as a consequence of a fractional translation symmetry present in some of the lattice structures.


Introduction
Formation of quantized vortices under rotation is a hallmark of superfluidity. When quantized vortices proliferate under rapid rotation, they organize into a regular lattice owing to their mutual repulsion. The resulting triangular vortex lattice structure was originally predicted by Abrikosov [1] for type-II superconductors in a magnetic field, and observed in superconducting materials [2], superfluid 4 He [3,4], and Bose-Einstein condensates (BEC) [5,6,7] and fermionic superfluids [8] of ultracold atoms. In ultracold atomic gases, in particular, the rotation frequency can be tuned over a wide range, and the equilibrium and dynamical properties of vortex lattices can be investigated in considerable detail [9,10,11]. Rotation can be viewed as the standard way to induce a synthetic gauge field for neutral atoms since the Hamiltonian in the rotating frame of reference is equivalent to that of charged particles in a uniform magnetic field. Notably, experimental techniques for producing synthetic gauge fields via optical dressing of atoms have also been developed over the past decade [12,13], and the very first experiment of these techniques successfully created around 10 vortices in a BEC without rotating the gas [14].
A BEC under rotation (or in a synthetic magnetic field) undergoes different regimes with increasing the rotation frequency Ω [11]. When a BEC rotates slowly, the size of the vortex core is much smaller than the intervortex separation. In this regime, the spatial variation of the BEC density, ∇|Ψ|, can be ignored, and the Thomas-Fermi (TF) approximation is applicable [15]. This regime is called the mean-field TF regime. With increasing Ω, the intervortex separation decreases and eventually becomes comparable with the vortex core size. Then the BEC flattens to an effectively two-dimensional (2D) system, and the interaction energy per particle becomes small compared with the kinetic energy per particle. It is thus reasonable to assume that atoms reside in the lowest-Landau-level (LLL) manifold for the motion in the 2D plane and to perform the mean-field calculation in this manifold [16,17]. This regime is called the mean-field LLL regime [10]. As Ω is further increased, the mean-field description breaks down, and the system is expected to enter a highly correlated regime. In particular, in a regime where the number of vortices N v becomes comparable with the number of atoms N , it has been predicted that the vortex lattice melts and a variety of quantum Hall states appear at integral and fractional values of the filling factor ν := N/N v [10,18,19].
A vortex lattice supports an elliptically polarized oscillatory mode, which was predicted by Tkachenko [20,21,22] and observed in superfluid 4 He [23]. While Tkachenko's original work predicted a linear dispersion relation for an incompressible fluid, a number of theoretical studies have been done to take into account a finite compressibility of the fluid [24,25,26,27,28]. It has been shown that the compressibility leads to hybridization with sound waves and qualitatively changes the dispersion relation into a quadratic form for small wave vectors. Collective modes of a vortex lattice have been observed over a wide range of rotation frequencies in a harmonically trapped BEC [29]. Theoretical analyses of the observed modes have been conducted with the hydrodynamic theory [30,31,32] and Gross-Pitaevskii (GP) mean-field calculations [33,34]. For a uniform BEC in the mean-field LLL regime, the dispersion relation of the Tkachenko mode can analytically be obtained within the Bogoliubov theory, and it is found to take a quadratic form [35,36,37]. Effective field theory for the Tkachenko mode has been developed in Refs. [38,39].
The properties of vortex latices can further be enriched in multicomponent BECs, such as those made up of different hyperfine spin states of identical atoms. For two-component BECs under rotation, GP mean-field calculations have shown that several different types of vortex lattices appear as the ratio of the intercomponent coupling g ↑↓ to the intracomponent one g > 0 is varied (see Fig. 1) [40,41,42]. Among them, interlaced square vortex lattices [ Fig. 1(d)], which are unique to these systems, have been observed experimentally [43]. Furthermore, optical dressing techniques can produce a variety of (possibly non-Abelian) gauge fields in multicomponent gases [12,13,44,45]. In particular, mutually antiparallel magnetic fields have been induced in two-component BECs, leading to the observation of the spin Hall effect [46]. If the antiparallel fields are made even higher, such systems are expected to show a rich phase diagram consisting of vortex lattices and (fractional) quantum spin Hall states [47,48,49]. Notably, it has been shown within the GP theory that BECs in antiparallel magnetic fields exhibit the same vortex lattice phase diagram as BECs in parallel magnetic fields [49] (see also Sec. 2.1). It is thus interesting to ask whether and how the difference between the two types of systems arises in other properties such as collective modes. In this context, it is worth noting that the two types of systems show markedly different phase diagrams in the quantum Hall regime [49,50,51,52,53], which has been interpreted in light of pseudopotentials and entanglement formation [53].
In this paper, we study collective modes of vortex lattices in two-component BECs in parallel and antiparallel synthetic magnetic fields in the mean-field LLL regime. On the basis of the Bogoliubov theory in the LLL basis, we numerically calculate excitation spectra for all the vortex-lattice structures shown in Fig. 1. We find that in all the cases, there appear two modes with linear and quadratic dispersion relations at low energies, which show anisotropy reflecting the symmetry of each lattice structure. Remarkably, the low-energy spectra for the two types of synthetic fields are related to each other by simple rescaling in the case of overlapping vortex lattices [ Fig. 1(a)] that appear for an intercomponent attraction. These results are consistent with an effective field theory analysis for low energies, which is a generalization of Ref. [38] aided with a symmetry consideration of the elastic energy of a vortex lattice. However, the rescaling relations are found to be violated for interlaced vortex lattices [ Fig. 1(b)-(e)] that appears for an intercomponent repulsion, presumably due to a nontrivial effect of a vortex displacement between the components beyond the effective field theory. We also find some interesting features of the excitation bands at high energies, such as line and point nodes, which arise from "fractional" translation symmetries or special structures of the Bogoliubov Hamiltonian matrix.
Here we comment on some related studies. Keçeli and Oktel [54] have studied collective excitation spectra in two-component BECs in parallel fields by means of the hydrodynamic theory, and predicted the appearance of two low-energy modes with linear and quadratic dispersion relations similar to ours. Our calculation is based on the Bogoliubov theory, provides unbiased results for weak interactions, and also contains information on the higher-energy part of the spectra. Furthermore, in the effective field theory analysis, we point out a term missing in Ref. [54], which is responsible for the anisotropy of the quadratic dispersion relation for interlaced triangular lattices [Fig.

1(b)]
. We also note that Woo et al. [55] have numerically investigated excitation spectra in rotating two-component BECs in a harmonic trap, and have identified a variety of excitations such as Tkachenko modes and surface waves.
The rest of this paper is organized as follows. In Sec. 2, we introduce the systems that we study in this paper, and formulate the problem by means of the Bogoliubov theory in the LLL basis. We then present our numerical results of Bogoliubov excitation spectra. In Sec. 3, we use an effective field theory to derive analytical formulae of lowenergy excitation spectra. In particular, we find the remarkable rescaling relations between the spectra for the two types of synthetic magnetic fields. In Sec. 4, we analyze the anisotropy of low-energy excitation spectra using the numerical data, and discuss its consistency with the effective field theory. In Sec. 5, we summarize the main results and discuss the outlook for future studies. In Appendix A, we describe the derivation of the matrix elements of the interaction used in Sec. 2. In Appendix B, we give precise definitions of the fractional translation operators used in Sec. 2. In Appendix C, we discuss some features of the Bogoliubov excitation spectra at high-symmetry points (found in Sec. 2) by using the data of the Bogoliubov Hamiltonian matrix. In Appendix D, we present a symmetry consideration of the elastic energy of vortex lattices, which is used in Sec. 3.

Bogoliubov analysis of excitation spectra
In this section, we introduce the systems that we study in this paper, and formulate the problem by means of the Bogoliubov theory in the LLL basis. Our formulation is closely related to those in Refs. [35,36,37]. In particular, the LLL magnetic Bloch states [37,56,57], which have a periodic pattern of zeros, play a crucial role here. We then present our numerical results of Bogoliubov excitation spectra and discuss their low-and high-energy characteristics.

Systems
We consider a system of a 2D pseudospin-1 2 Bose gas having two hyperfine spin states (labeled by α =↑, ↓). The spin-α component is subject to a synthetic magnetic field B α in the z direction. In the case of a gas rotating with an angular frequency Ω, parallel fields B ↑ = B ↓ = 2M Ω/q are induced in the two components in the rotating frame of reference, where M and q are the mass and the fictitious charge, respectively, of a neutral atom. An optical dressing technique of Ref. [46], in contrast, can be used to produce antiparallel fields B ↑ = −B ↓ . In the second-quantized form, the Hamiltonian of the system is given by where r = (x, y) is the coordinate on the 2D plane, p = −i (∂ x , ∂ y ) is the momentum, andψ α (r) is the bosonic field operator for the spin-α component satisfying the commutation relations [ψ α (r),ψ † β (r )] = δ αβ δ (2) (r − r ) and [ψ α (r),ψ β (r )] = [ψ † α (r),ψ † β (r )] = 0. The gauge field for the spin-α component is given by where we assume B > 0 and ↑ = ↓ = 1 ( ↑ = − ↓ = 1) for parallel (antiparallel) fields. For a 2D system of area A, the number of magnetic flux quanta piercing each component (or the number of vortices) is given by N v = A/(2π 2 ), where = /qB is the magnetic length. The total number of atoms is given by N = N ↑ + N ↓ , where N α is the number of spin-α bosons.
In the Hamiltonian (1), we assume contact interactions between atoms. For a gas tightly confined by a harmonic potential with frequency ω z in the z direction, the effective coupling constants in the 2D plane are given by g αα = a α 8π 3 ω z /M and g ↑↓ = g ↓↑ = a ↑↓ 8π 3 ω z /M , ‡ where a α and a ↑↓ are the s-wave scattering lengths between like and unlike bosons, respectively, in the 3D space. For simplicity, we set g ↑↑ = g ↓↓ ≡ g > 0 and N ↑ = N ↓ in the following. We further assume that the synthetic magnetic fields B α are sufficiently high or the interactions are sufficiently weak so that the energy scales of the interaction per atom, |g αβ |n, are much smaller than the Landaulevel spacing ω c := qB/M , where n := N ↑ /A = N ↓ /A is the density of atoms in each component. In this situation, it is legitimate to employ the LLL approximation in which the Hilbert space is restricted to the lowest Landau level [10,16,17].
When the filling factor ν ≡ N/N v is sufficiently high (ν 1), the system is well described by the GP mean-field theory. In this theory, the GP energy functional E[ψ ↑ , ψ ↓ ] is introduced by replacing the field operatorψ α (r) by the condensate wave function ψ α (r) in the Hamiltonian (1); then, the functional is minimized under the conditions d 2 r|ψ α | 2 = N α (α =↑, ↓) to determine the ground-state wave functions {ψ α (r)}. Using the LLL wave functions which have periodic patterns of zeros (and are equivalent to the LLL magnetic Bloch states described in Sec. 2.2), Mueller and Ho [40] have obtained a rich ground-state phase diagram for the parallel-field case, which consists of five different vortex-lattice structures as shown in the upper panels of Fig. 1. Notably, the GP energy functionals for the parallel-and antiparallel-field cases are related to each other as [49]. This implies that within the GP theory, the ground-state wave function of one case can be obtained from that of the other through the complex conjugation of the spin-↓ component. § Therefore, BECs in antiparallel fields also exhibit a rich variety of vortex-lattice structures as shown in Fig. 1 in the same way as BECs in parallel fields.

Lowest-Landau-level magnetic Bloch states
To describe the excitation properties of a vortex lattice, it is important to choose the basis consistent with the periodicity of the lattice. Following Refs. [37,56,57], we utilize the LLL magnetic Bloch states for this purpose. Let a 1 and a 2 be the primitive vectors of a vortex lattice as shown in Fig. 1(f). These vectors satisfy which implies the presence of one vortex in each component per unit cell. The reciprocal primitive vectors are then given by ‡ These are obtained by multiplying the coupling constants g (3D) α = 4π 2 a α /M and g (3D) ↑↓ = 4π 2 a ↑↓ /M for the 3D contact interactions by the factor M ω z /(2π ). This factor arises from the restriction to the ground state of the confinement potential in the z direction. § A similar situation arises for the ferromagnetic and antiferromagnetic Heisenberg models on a bipartite lattice, whose classical Hamiltonians are related to each other through the spin inversion S j → −S j on one of the two sublattices.
which satisfy a i · b j = 2πδ ij (i, j = 1, 2). Using the pseudomomentum for a spin-α particle we introduce the magnetic translation operator as T α (s) = e −iKα·s/ [58]. We note that the pseudomomentum K α = (K α,x , K α,y ) satisfies the commutation relation [K α,x , K α,y ] = −i α 2 / 2 . Starting from the most localized symmetric LLL wave function c 0 (r) = e −r 2 /4 2 / √ 2π 2 , we construct a set of LLL wave functions by multiplying two translation operators as commute with each other since one magnetic flux quantum pierces each unit cell as seen in Eq. (3); this property justifies the application of Bloch's theorem. By superposing c mα (r), we can construct the LLL magnetic Bloch state as [56] Ψ kα (r) = 1 with the normalization factor This state is an eigenstate of T α (a j ) with an eigenvalue e −ik·a j . The LLL magnetic Bloch state Ψ kα (r) represents a vortex lattice with a periodic pattern of zeros for any value of the wave vector k. Indeed, by rewriting Eq. (6) as and comparing it with the complex conjugate of the Perelomov overcompleteness equation m (−1) m 1 +m 2 c mα (r) = 0 [59], we find that Ψ kα (r) has zeros at [57] When one describes a triangular vortex lattice of a scalar BEC using a LLL magnetic Bloch state, the choice of the wave vector k is arbitrary once the primitive vectors a 1 and a 2 are set appropriately. This is because a change in k only leads to a translation of zeros as seen in Eq. (8). The vortex lattices of two-component BECs in Fig. 1 can Mueller and Ho [40] instead use the Jacobi theta function to express a vortex-lattice wave function. Such an expression is obtained by performing the Poisson resummation in Eq. (6) for N v → ∞. also be described by the LLL magnetic Bloch states Ψ qα,α (r) (α =↑, ↓); however, the wave vectors q ↑ and q ↓ have to be chosen in a way consistent with the displacement u 1 a 1 + u 2 a 2 between the components [see Fig. 1

(f)]. One useful choice is
Here, we displace the spin-↑ component by 1 2 (u 1 a 1 + u 2 a 2 ) and the spin-↓ component by − 1 2 (u 1 a 1 + u 2 a 2 ) instead of displacing only one of the components. This is useful for avoiding zeros of the normalization factor ζ(k) at some high-symmetry points in the first Brillouin zone [56]. ¶

Representation of the Hamiltonian
Using the magnetic Bloch states (6), we expand the field operator asψ α (r) = k Ψ kα (r)b kα , where k runs over the first Brillouin zone, and b kα is a bosonic annihilation operator satisfying [b kα , b † k α ] = δ kk δ αα . Substituting this expansion into the Hamiltonian, we obtain where ω c /2 is the LLL single-particle zero-point energy andN α = k b † kα b kα is the number operator for the spin-α component. The interaction matrix element V αβ (k 1 , k 2 , k 3 , k 4 ) is given by As described in Appendix A, this matrix element is calculated to be Here, δ P kk := G δ k,k +G is the periodic Kronecker's delta, where G runs over the reciprocal lattice vectors. In the case of parallel fields, the function S αβ (k 1 , k 2 , k 3 ) does not depend on α or β, and is given by ¶ If we set q ↑ = (b 1 − b 2 )/2 and q ↓ = 0 for square lattices, for example, we have ζ(q ↑ ) = 0 and Eq. (6) is not well-defined unless we factor out a nonanalytic dependence around the point of our concern [56].
In the case of antiparallel fields, S αβ (k 1 , k 2 , k 3 ) depends on α and β, and is given in terms of S(k 1 , k 2 , k 3 ) defined above by

Bogoliubov approximation
At high filling factors, the condensate is only weakly depleted and we can apply the Bogoliubov approximation [60,35,36,37]. + Provided that the condensation occurs at wave vector q α in the spin-α component, it is useful to introducẽ and retaining terms up to the second order inb kα andb † kα (k = 0), we obtain the following Bogoliubov Hamiltonian: Here, the matrix M(k) is given by + In the thermodynamic limit, however, this approximation is not valid since the fraction of quantum depletion diverges as 1 [35,37]. The Bogoliubov theory is still applicable since N v is at most of the order of 100 in experiments of ultracold atomic gases [7].
To diagonalize the Bogoliubov Hamiltonian (18), we perform the Bogoliubov transformation  where W (k) is a paraunitary matrix satisfying Equation (22) ensures the invariance of the bosonic commutation relation. If the matrix W (k) is chosen to satisfy the Bogoliubov Hamiltonian is diagonalized as By multiplying Eq. (23) from the left by W (k)τ 3 and using Eq. (22), one finds Therefore, the excitation energies E i (k) (i = 1, 2) can be obtained as the right eigenvalues of τ 3 M(k).

Numerical results
We use the formulation described above to numerically calculate the Bogoliubov excitation spectrum {E i (k)} in the following way. For a given wave vector k, we calculate the matrix M(k) in Eq. (19) by using Eqs. (12), (13), and (15). We note that each of the functions ζ(·) andζ(·) used in Eq. (13) involves an infinite sum but only with respect to two integer variables [see Eqs. (7) and (14)], which can numerically be taken with high accuracy. We then calculate the right eigenvalues of τ 3 M(k) to obtain {E i (k)}. Figure 2 presents the obtained energy spectra for all the lattice structures in Fig. 1 and for both the parallel-and antiparallel-field cases. In all the cases, we find that there appear two modes with linear and quadratic dispersion relations at low energies around the Γ point. Furthermore, we find anisotropy of the coefficients of these dispersion relations. For example, such anisotropy can clearly be seen along the path M 1 → Γ → R for (c) rhombic, (d) square, and (e) rectangular lattices. We discuss such anisotropy in detail in later sections. Apart from the low-energy features, the spectra in Fig. 2 also exhibit unique structures of band touching at some high-symmetry points or along lines in the Brillouin zone. In particular, the spectra for (c) rhombic, (d) square, and (e) rectangular lattices in parallel fields exhibit line nodes, whose locations in the Brillouin zones are shown in Fig.  2(f). This can be understood as a consequence of a "fractional" translation symmetry * [61,62]. Namely, in these cases, the system is invariant under the product T (P) of the translation by a 3 /2 and the spin reversal ↑↔↓, where a 3 := a 1 +a 2 . Since T (P) commutes with the Bogoliubov Hamiltonian and T (P) 2 gives the translation by a 3 , the Bloch states at k can be chosen to be the eigenstates of T (P) with T (P) |w ± k = ±e −ik·a 3 /2 |w ± k . For k → k + b i (i = 1, 2), the two eigenstates must switch places, indicating the occurrence of an odd number of degeneracies. In Fig. 2(f), we can indeed confirm that starting from any point other than the line nodes, the degeneracy occurs once or three times for the above changes of k. The emergence of point nodes at the M 1 and M 2 points for the same lattices (c), (d), and (e) for antiparallel fields can be understood by considering the symmetry under the product T (AP) of the time reversal and the translation by a 3 /2. Since T (AP) 2 is equal to the translation by a 3 , we have (T (AP) ) 2 = e −ik·a 3 in the subspace with the wave vector k. The Kramers degeneracy thus occurs at time-reversal-invariant momenta with e −ik·a 3 = 1, which is the case for k = b 1 /2 and b 2 /2 (M 1 and M 2 points). In Appendix C, we further discuss some other features of the spectra at high-symmetry points, such as the coincidence of the excitation energies between the two types of fields at the M 1 and M 2 points in Fig. 2(c), (d), and (e) by using the numerical data of the Bogoliubov Hamiltonian matrix M(k).

Effective field theory for low-energy excitation spectra
We have seen in the preceding section that vortex lattices of two-component BECs exhibit two excitation modes with linear and quadratic dispersion relations at low energies. Here we derive such low-energy dispersion relations by using an effective field theory. Specifically, we apply the formalism for a scalar BEC developed by Watanabe and Murayama [38] to the present two-component case. This approach is equivalent to the hydrodynamic theory applied by Keçeli and Oktel [54] to two-component BECs in parallel fields. However, we point out that an important term is missing in the elastic energy of vortex lattices used in Ref. [54]. This term is crucial for explaining the anisotropy of the quadratic dispersion relation for interlaced triangular lattices. Furthermore, we derive remarkable "rescaling" relations between the spectra for the two types of synthetic fields; these relations are confirmed for overlapping triangular lattices in Sec. 4.

Effective Lagrangian for phase variables
The Lagrangian density of the two-component BECs corresponding to the Hamiltonian (1) is given by [60] where ψ α (r, t) is the bosonic field for the spin-α component. To describe the low-energy properties of the BECs, it is useful to decompose the field as ψ α = √ n α exp(−iθ α ), where n α (r, t) and θ α (r, t) are the density and phase variables, respectively. Substituting this into Eq. (26) and keeping only the leading terms in the derivative expansion, we obtain where is an effective chemical potential for the spin-α component. Introducing n ± := n ↑ ± n ↓ and g ± := g ± g ↑↓ , we can rewrite Eq. (27) as By integrating out n ± (r, t), we obtain the effective Lagrangian for the phase variables {θ α (r, t)} as

Coupling between vortex displacement and phase variables
In the presence of vortices, the phase variables {θ α (r, t)} involve singularities. It is thus useful to decompose θ α into regular and singular parts as θ α = θ reg,α + θ sing,α . We also introduce the displacement of a vortex from the equilibrium position as u α (r, t) = r − X α (r, t), where r is the equilibrium position of the vortex and X α is the position at time t. The derivatives of the singular part θ sing,α of the phase can be related to the displacement u α as [38,39] where ij is an antisymmetric tensor with xy = − yx = +1. The effective chemical potential in Eq. (28) can then be expressed in terms of {θ reg,α , u α } as One should also note that the displacement u α (r, t) leads to a change in the elastic energy d 2 r E el (u α , ∂ i u α ). Here, the form of the elastic energy density E el depends on the type of a lattice as discussed in the next section and Appendix D. The effective Lagrangian in terms of {θ reg,α , u α } is then obtained as The ground state of H − µ 0 (N ↑ + N ↓ ) is given by θ reg = µ 0 t/ and u α = 0. To discuss the low-energy properties, it is therefore useful to introduce ϕ α = µ 0 t/ − θ reg,α and expand the Lagrangian (31) in terms of {ϕ α , u α }. Keeping only the quadratic terms in these variables, we obtain where ϕ ± := ϕ 1 ± ϕ 2 . Because {u α } have the mass term −u 2 α , one can expect that they can safely be integrated out in the discussion of low-energy dynamics. To do so, it is useful to derive the Euler-Lagrange equations for {u α }: Similar relations are also found in hydrodynamic theory [24,25,26,27,28,30,31,32,54]. Introducing u ± := u ↑ ± u ↓ and using the cyclotron frequency ω c = qB/M and the magnetic length = /qB, these relations can be rewritten as The second term on the left-hand side can be ignored in the LLL approximation (|g αβ |n ω c ). The relations (34) indicate that the vortex displacements u ± and the phases ϕ ± are coupled in opposite manners between the parallel-and antiparallelfield cases. Namely, the symmetric u + (antisymmetric u − ) is coupled to the symmetric ϕ + (antisymmetric ϕ − ) in parallel fields, while they are coupled in a crossed manner in antiparallel fields.
Substituting Eq. (34) into the Lagrangian density (32), we obtain the Lagrangian density in terms of {ϕ ± }, which can be used to determine the excitation spectrum. For this purpose, we need to determine the form of the elastic energy density E el , which is done next.

Elastic energy
Since the elastic energy is invariant under a uniform change in u + (r, t) (i.e., translation of the lattices), E el should be a function of ∂ i u + (i = x, y) and u − to the leading order in the derivative expansion. We therefore introduce the form To express E (+) el , it is useful to introduce w 0 := ∂ x u x + + ∂ y u y + , w 1 := ∂ x u x + − ∂ y u y + , w 2 := ∂ y u x + + ∂ x u y + .
In the LLL regime, the vortex density stays constant, and therefore w 0 = 0; this can also be confirmed by using Eq. (34). From a symmetry consideration (see Appendix D), each term in Eq. (35) can be expressed as where n := N ↑ /A = N ↓ /A is the average number density of each component. For each of the vortex lattices in Fig. 1(a)-(e), the dimensionless elastic constants Keçeli and Oktel [54] have considered an elastic energy consisting of E (±) el above, but have not included E (+−) el . Therefore, in their work, the symmetric and antisymmetric displacements u ± were decoupled from each other in collective modes. In our analysis, E (+−) el is found to be allowed by symmetry for interlaced triangular lattices. As we see below, this part crucially changes the low-energy spectrum, and explains the anisotropy of the spectrum for the concerned lattices.
We note that within the mean-field theory, the elastic energy density E el should take the same form [Eqs. (35) and (37)] for the parallel-and antiparallel-field cases because of the exact correspondence of the GP energy functionals between the two cases [49]. The dimensionless elastic constants are also expected to take the same values between the two cases. However, as we will see in Sec. 4, the elastic constants estimated from the numerical results of the energy spectra are different between the two cases. We discuss this puzzling issue in Sec. 4.2.

Excitation spectrum
The Lagrangian density in terms of ϕ ± is obtained by substituting Eq. (34) into Eq. (32) and using the above E el . After performing the Fourier transformation ϕ ± (r, t) = k dω 2π √ A e i(k·r−ωt) ϕ ± (k, ω), we obtain the action where (40) is the inverse of Green's function in Fourier space with In Eq. (40) [and Eqs. (42), (43), and (45) below], the upper and lower of the double signs correspond to the parallel-and antiparallel-field cases, respectively. The excitation spectrum corresponds to the poles of the Green's function, and can thus be obtained by solving the equation det[iG(k, ω) −1 ] = 0. Since Γ − (k) Γ(k) Γ + (k) for k 1, we obtain the low-energy dispersion relations as Using Eq. (41) and the fact that Γ − (k) is isotropic when F 1 = 0 [see Eq. (38)], we obtain the following explicit expressions with C 4 := F 1 2 /4D 1 . We thus find the emergence of quadratic and linear dispersion relations whose anisotropy reflects the symmetry of each lattice structure.
To discuss the anisotropy further, we parametrize the wave vector in terms of polar coordinates as k = k(cos θ, sin θ) (k 1) and introduce the dimensionless functions Using the dispersion relations (43) obtained from the effective field theory, these functions are calculated as In this result [and also in Eqs. (42) and (43)], the dependence on the type of synthetic fields occurs only in the coefficients g ± /g. This observation leads to the following remarkable relations: where the superscripts P and AP refer to the parallel-and antiparallel-field cases, respectively. Namely, the functions {f P/AP i (θ)} for the two types of synthetic fields are related to each other by simple rescaling. While these rescaling relations are expected for all the lattice structures within the effective theory, we see in the next section that the relations hold only for overlapping triangular lattices and break down for the other lattices.

Anisotropy of low-energy excitation spectra
We have seen in Sec. 2.5 that the Bogoliubov excitation spectrum exhibits linear and quadratic dispersion relations at low energies with significant anisotropy in some cases.
In this section, we analyze this anisotropy further by calculating the dimensionless functions {f i (θ)} defined in Eq. (44) for the cases shown in Fig. 2. We compare the numerical results with the analytical expressions (43) obtained by the effective field theory. We also examine whether the numerical results satisfy the rescaling relations (46) derived by the effective field theory.

Overlapping triangular lattices
For (a) overlapping triangular lattices, by using Eqs. (38) and (45), the analytic expressions of {f P/AP i (θ)} for parallel (P) and antiparallel (AP) fields are obtained as Notably, these functions show no dependence on θ in the effective field theory.
In numerical calculations, we obtain {f P/AP i (θ)} from the data of the Bogoliubov excitation spectra along a circular path k = k(cos θ, sin θ) with sufficiently small k and arbitrary θ ∈ [0, 2π). Figure 3(a) presents numerical results for g ↑↓ /g = −0.2. We find that the functions {f P/AP i (θ)} stay constant to a good accuracy consistent with the analytical expressions (47). The figure also shows the rescaled functions [the leftand right-hand sides of Eq. (46)], clearly demonstrating the rescaling relations (46). The dimensionless elastic constants C and D thus take the same values for the two types of fields and are plotted as functions of g ↑↓ /g in Fig. 4(a). Both constants are linear functions of g ↑↓ /g, which is consistent with the fact that the elastic energy is a linear function of g ↑↓ /g for a fixed vortex-lattice structure (see also Fig. 4 of Ref. [54]). Thus the numerical results are consistent with the effective field theory in the case of overlapping triangular lattices.  (θ) (right) for parallel (P; gray) and antiparallel (AP; pink) fields for the same cases as in Fig. 2(a)-(e). These functions are calculated from the Bogoliubov excitation spectra {E i (k)} along a circular path k = k(cos θ, sin θ) with k = 0.001a/ 2 and θ ∈ [0, 2π). The rescaled functions [left-and right-hand sides of Eq. (46)] are also shown (black and red), confirming the rescaling relations (47) through the overlap of the curves for (a) overlapping triangular lattices.

Interlaced lattices
We have performed similar analyses for interlaced lattices as shown in Fig. 3(b)-(e). The functions {f P/AP i (θ)} displayed in the figure show anisotropy except in the right panels for (b) interlaced triangular and (d) square lattices. These behaviors are consistent with the analytical results in Eqs. (38) and (45). Indeed, we can fit the numerical data perfectly using Eq. (45) if we determine C P/AP i (i = 1, 2, 3, 4) and D P/AP i (i = 1, 2, 3) separately for parallel or antiparallel fields. Figure 4 plots the determined constants {C i } and {D i }. We note that the constant C 4 , which is newly introduced in this work and originates from the coupling between the symmetric and antisymmetric vortex displacements u ± , is indeed finite for (b) interlaced triangular lattices.
However, the rescaling relations (46) derived from the effective field theory do not hold in Fig. 3(b)-(e). It can also be seen in different values of the constants {C P/AP i } and {D P/AP i } between the parallel-and antiparallel-field cases in Fig. 4(b)-(e). The difference between the two cases tends to increase with increasing the ratio g ↑↓ /g > 0. Furthermore, the constants for (d) square lattices show nonlinear dependences on g ↑↓ /g, which is inconsistent with the expected linear dependences for a fixed vortex-lattice structure (see Fig. 6 of Ref. [54]). These results cannot be explained within our effective field theory.
As discussed in the last paragraph of Sec. 3.3, the elastic constants should take the same values between the parallel-and antiparallel-field cases because of the exact correspondence of the GP energy functionals between the two cases [49]. Therefore, a possible insufficiency of our effective field theory would reside in how the elastic constants are related to the coefficients in the dispersion relations. We infer that the derivative expansions used in the derivation of the effective Lagrangian should be improved for interlaced vortex lattices which have a finite displacement between the components.

Summary and outlook
We have studied collective excitation modes of vortex lattices in two-component BECs subject to synthetic magnetic fields in parallel or antiparallel directions. Our motivation for studying the two types of synthetic fields stems from the fact that they lead to the same mean-field ground-state phase diagram [49] consisting of a variety of vortexlattice phases [40,41]-it is interesting to investigate what similarities and differences occur in collective modes. Our analyses are based on a microscopic calculation using the Bogoliubov theory and an analytical calculation using a low-energy effective field theory. We have found that there appear two distinct modes with linear and quadratic dispersion relations at low energies for all the lattice structures and for both types of synthetic fields. These dispersion relations show anisotropy that reflects the symmetry of each lattice structure. In particular, we have pointed out that the anisotropy of the quadratic dispersion relation for interlaced triangular lattices can be explained by a term in the elastic energy that mixes the symmetric and antisymmetric vortex displacements-such a term was missing in a previous study [54]. We have also found that the low-energy spectra for the two types of synthetic fields are related by simple rescaling in the case of overlapping triangular lattices that appears for intercomponent attraction (−1 < g ↑↓ /g < 0). However, contrary to the effective field theory prediction, such relations are found to break down for interlaced vortex lattices, which appear for intercomponent repulsion (g ↑↓ /g > 0) and involve a vortex displacement between the components. This indicates a nontrivial effect of an intercomponent vortex displacement on excitation properties that cannot be captured by the effective field theory developed in this paper. We have also found that the spectra exhibit unique structures of band touching at some high-symmetry points or along lines in the Brillouin zone. We have discussed their physical origins on the basis of fractional translation symmetries and the numerical data of the Bogoliubov Hamiltonian matrix.
The information of Bogoliubov excitation spectra studied in this work can be used to calculate the quantum correction to the ground-state energy due to zero-point fluctuations [see Eq. (24)], where the correction is expected to be enhanced as the filling factor ν is reduced. Despite the exact equivalence of the mean-field ground states between the parallel-and antiparallel-field cases [49], we have found quantitatively different Bogoliubov excitation spectra for the two cases as shown in Fig. 2. It is thus interesting to investigate how the quantum corrections affect the rich vortex lattice phase diagrams in the two cases. The present work would be a step toward understanding how the systems evolve from equivalent phase diagrams in the mean-field regime to markedly different phase diagrams in the quantum Hall regime [49,50,51,52,53] as the filling factor is lowered.

Equation (
A.1) can then be rewritten as Therefore, the interaction matrix element can be expressed by Eq. (12) with Let us focus on the case of parallel fields ( ↑ = ↓ = +1). In this case, the function S αβ (k 1 , k 2 , k 3 ) does not depend on α or β, and therefore we drop the subscript α, β. Using 4 2F (r n 1 , r n 2 , r n 3 ) = − j r 2 n j + (r n 2 · r n 3 + r n 1 · r n 3 ) − i(r n 2 × r n 3 + r n 1 × r n 3 ) z , we find S(k 1 , k 2 , k 3 ) = n (−1) n 1 n 2 exp −r 2 n /4 2 + ik 3 · r n × ζ k 1 + (r n × e z + ir n )/4 2 ζ k 2 + (r n × e z + ir n )/4 2 , where the sums over n 1 and n 2 are rewritten in terms of ζ(k) in Eq. (7), and the remaining dummy variable n 3 is replaced by n. We can further rewrite this by exploiting the following property of ζ(k) for s ∈ Z 2 : By setting n = 2s + p with s ∈ Z 2 and p ∈ {0, 1} 2 , Eq. (A.2) can be rewritten as In the case of antiparallel fields, S ↑↑ (k 1 , k 2 , k 3 ) is given by S(k 1 , k 2 , k 3 ) above. The other S αβ (k 1 , k 2 , k 3 )'s can be obtained by using the relation Ψ k↓ (r) = Ψ * −k↑ (r), leading to the result in Eq. (15).
Using T (a j )T = e −[K·a j ,K·a 3 /2]/ 2T T (a j ) = −T T (a j ) (j = 1, 2), one can confirm that Ψ k+q,↑ (r) defined in this way has the expected momentum: Furthermore, by operatingT on Ψ k+q,↑ (r), we havẽ Equations (B.3) and (B.4) indicate that the operatorT has the role of interchanging Ψ k−q,↓ (r) and Ψ k+q,↑ (r) with the multiplication of the same phase factor e −ik·a 3 /2 , which is a useful feature of the present basis. In this representation, one can show where the bar on α or β indicates the spin reversal ↑↔↓ and we use the invariance of the interaction g αβ δ (2) (r − r ) under the translation and the spin reversal (g αβ = gᾱβ). For a single particle, we define the fractional translation as the wave function changes byT in Eqs. (B.3) and (B.4) followed by the spin reversal σ x . For many particles, the fractional translation operator T (P) can be expressed in the second-quantized form as Using Eq. (B.5), one can confirm that the Bogoliubov Hamiltonian (18) is invariant under T (P) . The ground state |GS is obtained as the vacuum annihilated by the Bogolon annihilation operators γ k,j (j = 1, 2) in Eq. (21). The single-particle excitations γ † k,j |GS (j = 1, 2) can be used for the Bloch states |w ± k in the argument of Sec. 2.5.

Appendix B.2. Case of antiparallel fields
In the case of antiparallel fields ( ↑ = − ↓ = 1), we again modify the basis slightly from the magnetic Bloch states introduced in Sec. 2.2. While we use the same magnetic Bloch states as in Sec. 2.2 for the spin-↓ component, we define Ψ k+q,↑ (r) for the spin-↑ component viaT Appendix C. Bogoliubov Hamiltonian matrix at high-symmetry points In Sec. 2.5, we have discussed the origins of point and line nodes in the Bogoliubov excitation spectra in Fig. 2(c), (d), and (e) from the viewpoint of fractional translational symmetries. In Fig. 2, we further notice the following interesting features of the spectra at high-symmetry points: (i) coincidence of the excitation energies between the two types of fields at the M 1 and M 2 points for (c) rhombic, (d) square, and (e) rectangular lattices, (ii) the point node at the K 1 point for (a) overlapping triangular lattices in antiparallel fields, and (iii) the point node at the K 2 point for (b) interlaced triangular lattices in antiparallel fields. We have not succeeded in explaining these features from a symmetry viewpoint. Here, we instead discuss their origins on the basis of the numerical data of the Bogoliubov Hamiltonian matrix M(k). where the upper and lower of the double signs correspond to the parallel-and antiparallel-field cases, respectively, and "0" indicates elements whose numerical values vanish with high accuracy. The structure of the matrix indicates that the spin-↑ and ↓ components are completely decoupled at this wave vector. This remarkable decoupling may be understood pictorially as in Fig. C1; here, only the spin-↑ vortices show displacements, and the forces acting on each spin-↓ vortex from the surrounding spin-↑ vortices cancel out owing to the staggered nature of the displacements. Once the two components are decoupled in this way, they independently exhibit collective modes with identical spectra irrespective of the direction of the synthetic field. This explains the two-fold degeneracy of eigenenergies and the coincidence of those energies between the parallel-and antiparallel-field cases. Similar structures of the matrix M(k) are also seen at the This matrix consists of two independent blocks-a block corresponding to a spin-↑ particle and a spin-↓ hole and a block corresponding to a spin-↓ particle and a spin-↑ hole.
Since the two blocks have identical matrix elements, they show identical eigenenergies, which leads to the two-fold degeneracy at the K 1 point. In this matrix, there is no coupling between a particle and a hole or between spin-↑ and ↓ particles. Thus, the spin-↑ and ↓ particles independently exhibit excitation modes, which leads to the two-fold degeneracy at the K 2 point. In contrast, there is a coupling between the spin-↑ and ↓ holes, which leads to excitations with non-degenerate negative eigenenergies. Such excitations appear with non-degenerate positive eigenenergies at the K 1 point, which is seen in Fig. 2(b). Understanding how the matrix structures in Eqs. (C.2) and (C.3) emerge is still elusive.

Appendix D. Symmetry consideration of the elastic energy
Here we consider the elastic energy density E el (u α , ∂ i u α ) of the vortex lattices of twocomponent BECs shown in Fig. 1, and discuss how the symmetry constrains it into the form of Eqs. (35), (37), and (38). We start from the quadratic forms of w := (w 1 , w 2 ) t and u − : where C, D, and F are real 2 × 2 matrices, and C and D can be assumed to be symmetric. We assume that the vortex lattices have the symmetry under the coordinate Under this transformation, while u − is transformed by the same matrix Λ, w is, in general, transformed by a different matrixΛ. In order for the elastic energy to be invariant under this transformation, the following equations must be satisfied: and we finally obtain Eqs. (37) and (38).