Keywords: multicomponent Bose-Einstein condensates, synthetic gauge fields,
vortex lattices, quantum entanglement

# Vortex lattices in binary Bose-Einstein condensates: Collective modes,
quantum fluctuations, and intercomponent entanglement

Takumi Yoshino1, Shunsuke Furukawa2, Masahito Ueda1,3,4 1Department of
Physics, University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan
2Department of Physics, Keio University, 3-14-1 Hiyoshi, Kohoku-ku, Yokohama
223-8522, Japan
3RIKEN Center for Emergent Matter Science (CEMS), Wako, Saitama 351-0198,
Japan
4Institute for Physics of Intelligence, University of Tokyo, 7-3-1 Hongo,
Bunkyo-ku, Tokyo 113-0033, Japan<EMAIL_ADDRESS>and
<EMAIL_ADDRESS>
###### Abstract

We study binary Bose-Einstein condensates subject to synthetic magnetic fields
in mutually parallel or antiparallel directions. Within the mean-field theory,
the two types of fields have been shown to give the same vortex-lattice phase
diagram. We develop an improved effective field theory to study properties of
collective modes and ground-state intercomponent entanglement. Here, we point
out the importance of introducing renormalized coupling constants for coarse-
grained densities. We show that the low-energy excitation spectra for the two
kinds of fields are related to each other by suitable rescaling using the
renormalized constants. By calculating the entanglement entropy, we find that
for an intercomponent repulsion (attraction), the two components are more
strongly entangled in the case of parallel (antiparallel) fields, in
qualitative agreement with recent studies for a quantum (spin) Hall regime. We
also find that the entanglement spectrum exhibits an anomalous square-root
dispersion relation, which leads to a subleading logarithmic term in the
entanglement entropy. All of these are confirmed by numerical calculations
based on the Bogoliubov theory with the lowest-Landau-level approximation.
Finally, we investigate the effects of quantum fluctuations on the phase
diagrams by calculating the correction to the ground-state energy due to zero-
point fluctuations in the Bogoliubov theory. We find that the boundaries
between rhombic-, square-, and rectangular-lattice phases shift appreciably
with a decrease in the filling factor.

##  1 Introduction

Engineering synthetic gauge fields and observing their physical effects in
ultracold atomic gases have been a subject of great interest in recent years
[1, 2, 3, 4, 5]. While atomic gases are charge neutral, effective gauge fields
can be induced by rotating gases [6, 7] or optically dressing atoms [8].
Atomic Bose-Einstein condensates (BECs) in synthetic magnetic fields have
close analogy with type-II superconductors in magnetic fields. Real or
synthetic magnetic fields induce quantized vortices in these systems; when
vortices proliferate in high fields, they form a regular lattice pattern owing
to their repulsion, as originally predicted by Abrikosov [9]. The resulting
triangular vortex lattice structure has been observed experimentally in
rapidly rotating BECs [10, 11, 12]. A vortex lattice exhibits elliptically
polarized oscillations known as the Tkachenko mode [13, 14, 15, 16], which has
also been observed in a trapped BEC [12, 17]. In a uniform system, the
Tkachenko mode has a quadratic dispersion relation at low frequencies [18, 19,
20, 21, 22], and is understood as a Nambu-Goldstone mode associated with
spontaneously broken U(1) symmetry and magnetic translation symmetries [23].
For sufficiently high synthetic magnetic fields, atoms are expected to be in
the lowest Landau level. A key parameter in this regime is the filling factor
$\nu:=N/N_{\mathrm{v}}$, where $N$ is the number of atoms and $N_{\mathrm{v}}$
is the number of flux quanta piercing the system. While the Gross-Pitaevskii
(GP) mean-field theory is applicable for $\nu\gg 1$ [24], quantum fluctuations
become significant as $\nu$ is lowered [18, 20, 21, 22, 25]. Theory predicts
that when $\nu$ is below a critical value $\nu_{\mathrm{c}}$, the vortex
lattice melts and incompressible quantum Hall states appear at various integer
and fractional values of $\nu$ [6, 26, 27, 28]. Estimates of
$\nu_{\mathrm{c}}$ vary between $\nu_{\mathrm{c}}\simeq 2-6$ from exact
diagonalization [27, 29, 30] and $\nu_{\mathrm{c}}\simeq 5-14$ from a
Lindemann criterion [18, 20, 22, 27].

For binary (or pseudospin-$\frac{1}{2}$) BECs, which are populated in two
hyperfine spin states of the same atomic species, a richer variety of
synthetic gauge fields have been realized, such as a uniform magnetic field by
rotation [31], and spin-orbit couplings [32, 33, 34] and pseudospin-dependent
antiparallel magnetic fields [35] by optical dressing techniques. For binary
BECs under rotation, GP mean-field calculations have revealed that five
vortex-lattice phases appear as the ratio of the intercomponent coupling
$g_{\uparrow\downarrow}$ to the intracomponent one $g>0$ is varied (see Fig.
1) [36, 37, 38, 39]. Square vortex lattices [Fig. 1(d)] have indeed been
observed experimentally [31]. Meanwhile, a spin Hall effect due to pseudospin-
dependent Lorentz forces has been observed in binary BECs in antiparallel
magnetic fields [35]. For high antiparallel fields, theory predicts a rich
phase diagram that consists of vortex lattices and (fractional) quantum spin
Hall states [40, 41, 42]. Remarkably, within the GP mean-field theory, binary
BECs in antiparallel magnetic fields exhibit the same vortex phase diagram as
those in parallel magnetic fields (i.e., under rotation) [42]. This is because
the GP energy functionals as well as the time-independent GP equations for the
two systems are related to each other through the complex conjugation of the
spin-$\downarrow$ condensate wave function. It is interesting to further
investigate similarities and differences between the two systems. As the
systems in parallel and antiparallel magnetic fields are closely related with
bilayer quantum Hall systems [43] and quantum spin Hall systems [44],
respectively, a comparative study of these systems can make a link between the
two research fields. Such studies have been conducted on collective modes of
vortex lattices [45, 46] and phase diagrams in a quantum (spin) Hall regime
[42, 47, 48, 49, 50, 51, 52].

Figure 1:  Upper panels: vortex-lattice structures in the ground state of the
binary BECs under synthetic magnetic fields described by the Lagrangian in Eq.
(1) [36, 37, 38, 39] (see also Refs. [53, 54, 55, 56, 57, 58, 59]). Within the
GP mean-field theory, the systems in parallel and antiparallel fields exhibit
the same phase diagrams [42]. Five different structures appear as the
interaction ratio $g_{\uparrow\downarrow}/g$ is varied (see Fig. 10): (a)
overlapping triangular lattices, (b) interlaced triangular lattices, (c)
rhombic lattices, (d) square lattices, and (e) rectangular lattices. Black
(grey) circles show the vortex positions in the spin-$\uparrow$ ($\downarrow$)
component. As shown in (f), each lattice structure is specified by the
primitive vectors ${\mathbf{a}}_{1}=(a,0)$ and
${\mathbf{a}}_{2}=b(\cos\theta,\sin\theta)$ as well as the vortex displacement
$u_{1}{\mathbf{a}}_{1}+u_{2}{\mathbf{a}}_{2}$ in one component relative to the
other. The area of the unit cell is given by
$({\mathbf{a}}_{1}\times{\mathbf{a}}_{2})_{z}=ab\sin\theta=2\pi\ell^{2}$. The
angle $\theta$ and the aspect ratio $b/a$ vary continuously in the rhombic-
and rectangular-lattice phases, respectively, as shown later in Fig. 11. Lower
panels: the Brillouin zone for each lattice structure shown above. The vectors
${\mathbf{b}}_{1}$ and ${\mathbf{b}}_{2}$ show the reciprocal primitive
vectors, while uppercase letters show high-symmetry points. The excitation
spectra and the single-particle ES presented in Figs. 4 and 7 are calculated
along dotted arrows.

Collective excitation spectra can be different between the parallel- and
antiparallel-field cases as there is no correspondence between the two cases
for the time-dependent GP equations. Keçeli and Oktel [45] have calculated
excitation spectra in binary BECs under rotation by using a hydrodynamic
theory. In a previous work [46], we have conducted a comparative study of
excitation modes in the parallel- and antiparallel-field cases by means of the
Bogoliubov theory and an effective field theory. All the calculations in these
works are based on the lowest-Landau-level (LLL) approximation. For both types
of fields and for all the lattice structures in Fig. 1, it is found that there
appear two distinct modes with linear and quadratic dispersion relations at
low energies, which exhibit anisotropy reflecting the symmetry of each lattice
structure. Furthermore, for overlapping triangular lattices [Fig. 1(a)], the
low-energy spectra for the two types of fields are found to be related to each
other by simple rescaling [46]. This indicates a nontrivial correspondence
between the two types of fields in excitation properties. However, such
rescaling relations do not hold for other lattices, which appears inconsistent
with the effective field theory prediction. A more refined description of low-
energy modes has remained an open issue.

Numerical studies on the quantum (spin) Hall regime with $\nu=O(1)$ have
revealed that binary Bose gases in parallel and antiparallel synthetic
magnetic fields exhibit markedly different phase diagrams [42, 47, 48, 49, 50,
51, 52]. For parallel fields, the product states of a pair of nearly
uncorrelated quantum Hall states are robust against an intercomponent
attraction $g_{\uparrow\downarrow}<0$ and persist even when
$|g_{\uparrow\downarrow}|$ is close to $g$ [52]. Meanwhile, a variety of spin-
singlet quantum Hall states with high intercomponent entanglement emerge for
$g_{\uparrow\downarrow}\approx g$ [47, 48, 49, 50, 51]. For antiparallel
fields, (fractional) quantum spin Hall states approximated by products of
quantum Hall states with opposite chiralities are robust against an
intercomponent repulsion $g_{\uparrow\downarrow}>0$ [42]. The phase diagrams
for the two types of fields thus exhibit opposite behavior in view of
intercomponent entanglement. An interpretation of these results has been given
in light of pseudopotential representations of interactions [52]. As there is
no intercomponent entanglement in the GP mean-field theory valid for $\nu\gg
1$, it is interesting to investigate how the intercomponent entanglement
arises in the two systems as $\nu$ is lowered from the mean-field regime.

In this paper, we present a detailed comparative study of vortex lattices of
binary BECs in parallel and antiparallel fields concerning ground-state and
excitation properties. We first formulate an improved effective field theory
for such vortex lattices, and derive some properties of collective modes and
ground-state intercomponent entanglement. Here, a major improvement is the
introduction of renormalized coupling constants for coarse-grained
densities.111 In the formulation of Ref. [46], the necessity of
renormalization was not transparent. This was because we integrated out the
density variables, before coarse graining, to obtain an effective Lagrangian
for phase variables. In this paper, we keep both the phase and density
variables in the Lagrangian, and perform coarse graining of both the variables
at the same time. The renormalization of density-density interactions can
naturally be understood in this formulation. This formulation also allows us
to move easily to the operator formalism, in which ground-state and excitation
properties can be studied in an algebraic manner.  We show that the low-energy
excitation spectra for the two types of fields are related to each other by
suitable rescaling using the renormalized constants. Namely, the rescaling
relations proposed in our previous work [46] must be modified using the
renormalized constants. By calculating the entanglement entropy (EE), we find
that for an intercomponent repulsion (attraction), the two components are more
strongly entangled in the case of parallel (antiparallel) fields, in
qualitative agreement with recent exact diagonalization results for a quantum
(spin) Hall regime [42, 52]. As a by-product, we also find that the
entanglement spectrum (ES) exhibits an anomalous square-root dispersion
relation, and that the EE exhibits a volume-law scaling with a subleading
logarithmic correction. This anomalous feature of the ES is associated with
the emergence of a long-range interaction in terms of the density in the
entanglement Hamiltonian. All these predictions are confirmed by numerical
calculations based on the Bogoliubov theory with the LLL approximation.
Finally, we investigate the effects of quantum fluctuations on the phase
diagrams by calculating the Lee-Huang-Yang correction, which is a correction
to the ground-state energy due to zero-point fluctuations in the Bogoliubov
theory [60]. We find that the boundaries between rhombic-, square-, and
rectangular-lattice phases shift appreciably with a decrease in $\nu$. Here,
the shift occurs more significantly for parallel fields.

Let us comment on the relation to another recent work of our own [61]
concerning intercomponent ES and EE in binary BECs in $d$ spatial dimensions
in the absence of synthetic gauge fields. Here we employ effective field
theory to show that the ES exhibits a gapless square-root dispersion relation
in the presence of an intercomponent tunneling (a Rabi coupling) and a gapped
dispersion relation in its absence (see also Refs. [62, 63] for related
results in two coupled Tomonaga-Luttinger liquids). In the present work, in
contrast, the ES exhibits a square-root dispersion relation in the absence of
an intercomponent tunneling in both the parallel- and antiparallel-field
cases. This qualitative distinction is related to the fact that binary BECs in
parallel and antiparallel fields have a higher density of low-energy
excitations and thus experience larger quantum fluctuations than those without
synthetic gauge fields. References [61, 63] and the present paper demonstrate
that a variety of long-range interactions can be emulated in a subsystem of
multicomponent BECs that have only short-range interactions. We also note that
the field-theoretical methods for investigating intercomponent entanglement in
the present paper are closely analogous to those in Refs. [61, 62, 63].

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 low-energy effective
field theory to derive some properties of excitation spectra, intercomponent
entanglement, and correlation functions. In Sec. 3, we briefly review the
Bogoliubov theory with the LLL approximation, which has been adapted to the
present problem in Ref. [46], and then give the expressions of the
intercomponent ES and EE. In Sec. 4, we present numerical results based on the
Bogoliubov theory. We confirm the field-theoretical predictions on the
excitation spectra, the intercomponent entanglement, and the fraction of
quantum depletion. Furthermore, we investigate the effects of quantum
fluctuations on the ground-state phase diagrams. In Sec. 5, we present a
summary of the present study and an outlook for future studies. In Appendices,
we describe some technical details of Secs. 2.

##  2 Effective field theory

Effective field theory for a vortex lattice in a scalar BEC has been developed
in Refs. [18, 23, 64, 65]. In our previous work [46], we applied the
formulation by Watanabe and Murayama [23] to binary BECs in parallel and
antiparallel fields. For parallel fields, this approach was essentially
equivalent to the hydrodynamic theory by Keçeli and Oktel [45]. Here, we
formulate an improved effective field theory for binary BECs in parallel and
antiparallel fields by introducing renormalized coupling constants
$\bar{g}_{\alpha\beta}$ for coarse-grained densities. In doing so, we keep
both the phase and density variables in the Lagrangian, rather than
integrating out the density variables as in Ref. [46]. We then move to the
operator formalism, and study excitation spectra, intercomponent entanglement,
and correlation functions in an algebraic manner (see Refs. [61, 62, 63] for
analogous calculations in different systems).

###  2.1 Systems

We consider a system of 2D binary (pseudospin-$\frac{1}{2}$) BECs having two
hyperfine spin states (labeled by $\alpha=\uparrow,\downarrow$) and subject to
synthetic magnetic fields $B_{\uparrow}$ and $B_{\downarrow}$ in mutually
parallel or antiparallel directions. The Lagrangian density of the system is
given by

$\mathcal{L}=\sum_{\alpha=\uparrow,\downarrow}\left[\frac{i\hbar}{2}(\psi_{\alpha}^{*}\dot{\psi}_{\alpha}-\dot{\psi}_{\alpha}^{*}\psi_{\alpha})-\frac{1}{2M}|(-i\hbar\nabla-q\mathbf{A}_{\alpha})\psi_{\alpha}|^{2}\right]-\sum_{\alpha,\beta=\uparrow,\downarrow}\frac{g_{\alpha\beta}}{2}|\psi_{\alpha}|^{2}|\psi_{\beta}|^{2},$
(1)

where $\psi_{\alpha}({\mathbf{r}},t)$ is the bosonic field for the
spin-$\alpha$ component (with ${\mathbf{r}}=(x,y)$ being the 2D coordinate)
and $M$ and $q$ are the mass and the fictitious charge of an atom. The gauge
field $\mathbf{A}_{\alpha}({\mathbf{r}})$ for spin-$\alpha$ bosons is given by

$\mathbf{A}_{\alpha}=\frac{B_{\alpha}}{2}{\mathbf{e}}_{z}\times{\mathbf{r}}=\epsilon_{\alpha}\frac{B}{2}(-y,x),$
(2)

where we assume $qB>0$ and $\epsilon_{\uparrow}=\epsilon_{\downarrow}=1$
$(\epsilon_{\uparrow}=-\epsilon_{\downarrow}=1)$ for parallel (antiparallel)
fields. The number of magnetic flux quanta piercing each component is given by
$N_{\mathrm{v}}=qBA/(2\pi\hbar)=A/(2\pi\ell^{2})$, where $A$ is the area of
the system and $\ell:=\sqrt{\hbar/qB}$ is the magnetic length. In the
Lagrangian (1), the numbers of spin-$\uparrow$ and $\downarrow$ atoms,
$N_{\uparrow}$ and $N_{\downarrow}$, are separately conserved. We introduce
the total filling factor $\nu:=N/N_{\mathrm{v}}$, where
$N:=N_{\uparrow}+N_{\downarrow}$ is the total number of atoms.

We assume a contact interaction between atoms. For simplicity, we set
$g_{\uparrow\uparrow}=g_{\downarrow\downarrow}\equiv
g>|g_{\uparrow\downarrow}|$ and $N_{\uparrow}=N_{\downarrow}$ in the
following. With these conditions, the system in parallel fields is invariant
under the interchange of the two components, while the system in antiparallel
fields is invariant under time reversal. To apply the LLL approximation, we
further assume that the scale of the interaction energy per atom,
$|g_{\alpha\beta}|n$, is much smaller than the Landau-level spacing
$\hbar\omega_{\mathrm{c}}$. Here, $n:=N_{\uparrow}/A=N_{\downarrow}/A$ is the
average density of $\uparrow$ or $\downarrow$ atoms, and
$\omega_{\mathrm{c}}:=qB/M$ is the cyclotron frequency.

###  2.2 Effective field theory

To obtain a low-energy effective description, it is useful to rewrite the
field as
$\psi_{\alpha}={\mathrm{e}}^{-\mathrm{i}\theta_{\alpha}}\sqrt{n_{\alpha}}$,
where $n_{\alpha}({\mathbf{r}},t)$ and $\theta_{\alpha}({\mathbf{r}},t)$ are
the density and phase variables, respectively. The Lagrangian density (1) is
rewritten in terms of these variables as

$\mathcal{L}=\sum_{\alpha}\left[\hbar
n_{\alpha}\dot{\theta}_{\alpha}-\frac{n_{\alpha}}{2M}(\hbar\nabla\theta_{\alpha}+q\mathbf{A}_{\alpha})^{2}-\frac{\hbar^{2}(\nabla
n_{\alpha})^{2}}{8n_{\alpha}M}\right]-\sum_{\alpha,\beta}\frac{g_{\alpha\beta}}{2}n_{\alpha}n_{\beta}.$
(3)

In the presence of vortices, the phase variables
$\\{\theta_{\alpha}({\mathbf{r}},t)\\}$ involve singularities. This motivates
us to decompose $\theta_{\alpha}$ into regular and singular contributions as
$\theta_{\alpha}=\theta_{\mathrm{reg},\alpha}+\theta_{\mathrm{sing},\alpha}$.
We also introduce the displacement $\mathbf{u}_{\alpha}({\mathbf{r}},t)$ of a
vortex from the equilibrium position. The derivatives of the singular part
$\theta_{\mathrm{sing},\alpha}$ of the phase can be related to the
displacement field ${\mathbf{u}}_{\alpha}$ as [23, 64]

$\hbar\dot{\theta}_{\mathrm{sing},\alpha}=-\frac{qB_{\alpha}}{2}({\mathbf{u}}_{\alpha}\times\dot{{\mathbf{u}}}_{\alpha})_{z},~{}~{}\hbar\nabla\theta_{\mathrm{sing},\alpha}+q{\mathbf{A}}_{\alpha}=qB_{\alpha}{\mathbf{e}}_{z}\times{\mathbf{u}}_{\alpha}-\frac{qB_{\alpha}}{2}\epsilon_{ij}u^{i}_{\alpha}\nabla
u^{j}_{\alpha},$ (4)

where $\epsilon_{ij}$ is an antisymmetric tensor with
$\epsilon_{xy}=-\epsilon_{yx}=+1$. The displacement
${\mathbf{u}}_{\alpha}({\mathbf{r}},t)$ also results in a change in the
elastic energy
$\int\mathrm{d}^{2}{\mathbf{r}}~{}\mathcal{E}_{\mathrm{el}}({\mathbf{u}}_{\alpha},\partial_{i}{\mathbf{u}}_{\alpha})$,
whose explicit form will be given in Sec. 2.3. The Lagrangian density can then
be expressed in terms of
$\\{\theta_{\mathrm{reg},\alpha},{\mathbf{u}}_{\alpha},n_{\alpha}\\}$ as

$\begin{split}\mathcal{L}=&\sum_{\alpha}\Biggl{[}\hbar
n_{\alpha}\dot{\theta}_{\mathrm{reg},\alpha}-\frac{qB_{\alpha}n_{\alpha}}{2}({\mathbf{u}}_{\alpha}\times\dot{{\mathbf{u}}}_{\alpha})_{z}-\frac{n_{\alpha}}{2M}\left(\hbar\nabla\theta_{\mathrm{reg},\alpha}+qB_{\alpha}\mathbf{e}_{z}\times{\mathbf{u}}_{\alpha}-\frac{qB_{\alpha}}{2}\epsilon_{ij}u_{\alpha}^{i}\nabla
u_{\alpha}^{j}\right)^{2}\\\ &~{}~{}~{}~{}~{}-\frac{\hbar^{2}(\nabla
n_{\alpha})^{2}}{8n_{\alpha}M}\Biggr{]}-\sum_{\alpha,\beta}\frac{\bar{g}_{\alpha\beta}}{2}n_{\alpha}n_{\beta}-\mathcal{E}_{\mathrm{el}}.\end{split}$
(5)

Henceforth, we ignore the term
$-\frac{qB_{\alpha}}{2}\epsilon_{ij}u_{\alpha}^{i}\nabla u_{\alpha}^{j}$ in
the round brackets as it only gives more than quadratic contributions to
$\mathcal{L}$ in terms of $\nabla\theta_{\mathrm{reg},\alpha}$ and
${\mathbf{u}}_{\alpha}$. We also omit the subscript “reg” in
$\theta_{\mathrm{reg},\alpha}$.

In the effective Lagrangian density (5), we have introduced renormalized
coupling constants $\bar{g}_{\alpha\beta}$ with
$\bar{g}_{\uparrow\uparrow}=\bar{g}_{\downarrow\downarrow}\equiv\bar{g}$ and
$\bar{g}_{\uparrow\downarrow}=\bar{g}_{\downarrow\uparrow}$. The necessity of
this renormalization has been overlooked in previous studies [45, 46] and can
be explained as follows. The introduction of the vortex displacement fields
$\\{u_{\alpha}({\mathbf{r}},t)\\}$ necessarily involves the coarse graining of
the theory. Namely, we smooth out details within the scale of the lattice
constants and instead focus on the physics at larger scales. Therefore, the
density variable $n_{\alpha}({\mathbf{r}},t)$ should now be understood as the
density averaged over the unit cell containing ${\mathbf{r}}$. The
renormalized coupling constants $\bar{g}_{\alpha\beta}$ can then be introduced
through the relation

$\int_{\mathrm{u.c.}}\mathrm{d}^{2}{\mathbf{r}}^{\prime}~{}g_{\alpha\beta}|\psi_{\alpha}({\mathbf{r}}^{\prime},t)|^{2}|\psi_{\beta}({\mathbf{r}}^{\prime},t)|^{2}=2\pi\ell^{2}\bar{g}_{\alpha\beta}n_{\alpha}({\mathbf{r}},t)n_{\beta}({\mathbf{r}},t),$
(6)

where $|\psi_{\alpha}|^{2}$ and $n_{\alpha}$ are the original and coarse-
grained densities of the spin-$\alpha$ atoms, respectively, and the
integration on the left-hand side is taken over the unit cell (with the area
$2\pi\ell^{2}$) that contains ${\mathbf{r}}$. As explained later in this
section and in Sec. 4.1, the renormalized constants $\bar{g}_{\alpha\beta}$
for describing the low-energy physics can be determined by calculating the
contribution of each interaction term to the mean-field ground-state energy.
As seen in Eq. (6), a positive (negative) correlation between the density
fluctuations $|\psi_{\alpha}|^{2}-n_{\alpha}$ and
$|\psi_{\beta}|^{2}-n_{\beta}$ leads to an enhanced (reduced) renormalized
coupling $|\bar{g}_{\alpha\beta}|>|g_{\alpha\beta}|$
($|\bar{g}_{\alpha\beta}|<|g_{\alpha\beta}|$). In particular, the
intracomponent coupling $g$ is always enhanced by the renormalization while
the intercomponent repulsion $g_{\uparrow\downarrow}>0$ (attraction
$g_{\uparrow\downarrow}<0$) is reduced (enhanced) by the renomalization owing
to the displacement (overlap) of vortices between the components; this can be
confirmed in Fig. 2 below. For a rotating scalar BEC with a repulsive coupling
$g>0$, a similar renormalization with $\bar{g}/g=1.1596$ has been discussed in
Refs. [7, 18, 66]. At high filling factors $\nu\gg 1$ and for not too large a
number of vortices $N_{\mathrm{v}}$, we can assume that the condensates are
only weakly depleted; we therefore have $|n_{\alpha}({\mathbf{r}},t)-n|\ll n$
for the coarse-grained density $n_{\alpha}$.

Because the displacement fields $\\{{\mathbf{u}}_{\alpha}\\}$ involve the mass
term $-{\mathbf{u}}_{\alpha}^{2}$ in Eq. (5), one may safely integrate them
out in the discussion of low-energy dynamics. Instead of performing the
integration directly, it is useful to derive the Euler-Lagrange equations for
$\\{{\mathbf{u}}_{\alpha}\\}$:

${\mathbf{u}}_{\alpha}-\epsilon_{\alpha}\ell^{2}{\mathbf{e}}_{z}\times\nabla\theta_{\alpha}-\frac{\epsilon_{\alpha}}{\omega_{\mathrm{c}}}{\mathbf{e}}_{z}\times\dot{{\mathbf{u}}}_{\alpha}+\frac{\ell^{2}}{n\hbar\omega_{\mathrm{c}}}\left[\frac{\partial\mathcal{E}_{\mathrm{el}}}{\partial{\mathbf{u}}_{\alpha}}-\partial_{j}\left(\frac{\partial\mathcal{E}_{\mathrm{el}}}{\partial\left(\partial_{j}{\mathbf{u}}_{\alpha}\right)}\right)\right]=0,$
(7)

where we have made the approximation $n_{\alpha}\approx n$. We can ignore the
third and fourth terms on the left-hand side in the LLL approximation with
$\hbar\omega,\mathcal{E}_{\mathrm{el}}/n\ll\hbar\omega_{\mathrm{c}}$, where
$\omega$ is the frequency of our concern. Introducing
${\mathbf{u}}_{\pm}:={\mathbf{u}}_{\uparrow}\pm{\mathbf{u}}_{\downarrow}$ and
$\theta_{\pm}:=\theta_{\uparrow}\pm\theta_{\downarrow}$, we can rewrite Eq.
(7) as

${\mathbf{u}}_{\pm}=\left\\{\begin{aligned}
&\ell^{2}\mathbf{e}_{z}\times\nabla\theta_{\pm}\qquad\text{(parallel
fields)};\\\
&\ell^{2}\mathbf{e}_{z}\times\nabla\theta_{\mp}\qquad\text{(antiparallel
fields)}.\end{aligned}\right.\qquad\text{}$ (8)

These relations indicate that the (anti)symmetric movement of vortices is
coupled to the (anti)symmetric component of the phase variables in parallel
fields while they are coupled in a crossed manner in antiparallel fields. By
substituting Eq. (8) into Eq. (5), we obtain the effective Lagrangian in terms
of $\theta_{\pm}$ and their conjugate momenta $\hbar
n_{\pm}=\hbar(n_{\uparrow}\pm n_{\downarrow})/2$. The Hamiltonian density is
then obtained as

$\mathcal{H}=\sum_{\nu=\pm}\left[\bar{g}_{\nu}n_{\nu}^{2}+\frac{\hbar^{2}(\nabla
n_{\nu})^{2}}{4nM}\right]+\mathcal{E}_{\mathrm{el}},~{}~{}~{}\bar{g}_{\pm}:=\bar{g}\pm\bar{g}_{\uparrow\downarrow}.$
(9)

where $\mathcal{E}_{\mathrm{el}}$ is expressed in terms of the phase variables
$\theta_{\pm}$ by using Eq. (8). The theory can be quantized by requiring the
canonical commutation relations

$[\theta_{\nu}({\mathbf{r}}),n_{\nu^{\prime}}({\mathbf{r}}^{\prime})]=\mathrm{i}\delta_{\nu\nu^{\prime}}\delta({\mathbf{r}}-{\mathbf{r}}^{\prime})~{}~{}~{}(\nu,\nu^{\prime}=\pm).$
(10)

In the present coarse-grained description, the mean-field ground state
corresponds to the uniform state with $n_{+}({\mathbf{r}})=n$,
$n_{-}({\mathbf{r}})=0$, and $\nabla\theta_{\pm}({\mathbf{r}})=0$. Therefore,
using Eq. (9), we obtain the mean-field ground-state energy density as

$\frac{E_{\mathrm{GS}}^{\mathrm{MF}}}{A}=\bar{g}_{+}n^{2}=(\bar{g}+\bar{g}_{\uparrow\downarrow})n^{2}=\frac{1}{2}\sum_{\alpha,\beta}\bar{g}_{\alpha\beta}n^{2}.$
(11)

In contrast, the same energy density obtained by the microscopic calculation
[the first term of Eq. (67) shown later] has the form
$E_{\mathrm{GS}}^{\mathrm{MF}}/A=(\beta
g+\beta_{\uparrow\downarrow}g_{\uparrow\downarrow})n^{2}$, where $\beta$ and
$\beta_{\uparrow\downarrow}$ are dimensionless constants that depend on the
lattice structure. The renormalized coupling constants are thus determined as
$\bar{g}=\beta g$ and
$\bar{g}_{\uparrow\downarrow}=\beta_{\uparrow\downarrow}g_{\uparrow\downarrow}$.

###  2.3 Elastic energy

The expression of the elastic energy density $\mathcal{E}_{\mathrm{el}}$ has
been determined in Refs. [45, 46], and we summarize it in the following. We
first note that the elastic energy must be invariant under a constant change
in ${\mathbf{u}}_{\alpha}({\mathbf{r}},t)$, i.e., translation of the lattices.
Therefore, to the leading order in the derivative expansion,
$\mathcal{E}_{\mathrm{el}}$ should be a function of
$\partial_{i}{\mathbf{u}}_{+}~{}(i=x,y)$ and ${\mathbf{u}}_{-}$, resulting in
the decomposition

$\mathcal{E}_{\mathrm{el}}=\mathcal{E}_{\mathrm{el}}^{(+)}(\partial_{i}{\mathbf{u}}_{+})+\mathcal{E}_{\mathrm{el}}^{(-)}({\mathbf{u}}_{-})+\mathcal{E}_{\mathrm{el}}^{(+-)}(\partial_{i}{\mathbf{u}}_{+},{\mathbf{u}}_{-}).$
(12)

To express $\mathcal{E}_{\mathrm{el}}^{(+)}$, it is useful to introduce

$w_{1}:=\partial_{x}u_{+}^{x}-\partial_{y}u_{+}^{y},~{}~{}w_{2}:=\partial_{y}u_{+}^{x}+\partial_{x}u_{+}^{y}.$
(13)

On the basis of a symmetry consideration, each term in Eq. (12) can be
expressed as

$\displaystyle\mathcal{E}_{\mathrm{el}}^{(+)}(\partial_{i}{\mathbf{u}}_{+})=\frac{gn^{2}}{2}\left(C_{1}w_{1}{}^{2}+C_{2}w_{2}{}^{2}+C_{3}w_{1}w_{2}\right),$
(14a)
$\displaystyle\mathcal{E}_{\mathrm{el}}^{(-)}({\mathbf{u}}_{-})=\frac{gn^{2}}{2\ell^{2}}\left[D_{1}\left(u_{-}^{x}\right)^{2}+D_{2}\left(u_{-}^{y}\right)^{2}+D_{3}u_{-}^{x}u_{-}^{y}\right],$
(14b)
$\displaystyle\mathcal{E}_{\mathrm{el}}^{(+-)}(\partial_{i}{\mathbf{u}}_{+},{\mathbf{u}}_{-})=\frac{gn^{2}}{2\ell}F_{1}\left(w_{1}u_{-}^{y}+w_{2}u_{-}^{x}\right).$
(14c)

For each of the vortex-lattice structures in Fig. 1(a)-(e), the dimensionless
elastic constants $\\{C_{1},C_{2},C_{3},D_{1},D_{2},D_{3},F_{1}\\}$ satisfy

$\begin{split}&\text{(a)}~{}C_{1}=C_{2}\equiv C>0,~{}D_{1}=D_{2}\equiv
D>0,~{}C_{3}=D_{3}=F_{1}=0;\\\ &\text{(b)}~{}C_{1}=C_{2}\equiv
C>0,D_{1}=D_{2}\equiv D>0,~{}C_{3}=D_{3}=0,~{}F_{1}\neq 0;\\\
&\text{(c)}~{}C_{1},C_{2},D_{1},D_{2}>0,~{}C_{3},D_{3}\neq 0,~{}F_{1}=0;\\\
&\text{(d)}~{}C_{1},C_{2}>0,~{}D_{1}=D_{2}\equiv
D>0,~{}C_{3}=D_{3}=F_{1}=0;\\\
&\text{(e)}~{}C_{1},C_{2}>0,~{}D_{1},D_{2}>0,~{}C_{3}=D_{3}=F_{1}=0.\end{split}$
(15)

Keçeli and Oktel [45] have determined the constants
$\\{C_{1},C_{2},C_{3},D_{1},D_{2},D_{3}\\}$ by calculating a change in the
mean-field energy under deformation of vortex lattices. In our previous work
[46], we have pointed out the presence of the $F_{1}$ term for interlaced
triangular vortex lattices (b), which was missed in Ref. [45]. For $\nu\gg 1$,
the elastic constants should be the same for the two types of fields because
of the exact correspondence of the GP energy functionals [42]. In Sec. 4.2, we
will determine all the elastic constants above as a function of
$g_{\uparrow\downarrow}/g$ by using the data of the excitation spectra.

###  2.4 Diagonalization of the effective Hamiltonian

We are now in a position to calculate the energy spectrum of the Hamiltonian
$H=\int d^{2}{\mathbf{r}}~{}{\cal H}$, where the Hamiltonian density ${\cal
H}$ is given by Eq. (9). We perform Fourier expansions

$\theta_{\nu}({\mathbf{r}})=\frac{1}{\sqrt{A}}\sum_{\mathbf{k}}\theta_{{\mathbf{k}},\nu}{\mathrm{e}}^{\mathrm{i}{\mathbf{k}}\cdot{\mathbf{r}}},~{}~{}n_{\nu}({\mathbf{r}})=\frac{1}{\sqrt{A}}\sum_{\mathbf{k}}n_{{\mathbf{k}},\nu}{\mathrm{e}}^{\mathrm{i}{\mathbf{k}}\cdot{\mathbf{r}}}~{}~{}(\nu=\pm),$
(16)

where the Fourier components satisfy

$[\theta_{{\mathbf{k}},\nu},n_{-{\mathbf{k}}^{\prime},\nu^{\prime}}]=\mathrm{i}\delta_{\nu\nu^{\prime}}\delta_{{\mathbf{k}}{\mathbf{k}}^{\prime}},~{}~{}\theta^{\dagger}_{{\mathbf{k}},\nu}=\theta_{-{\mathbf{k}},\nu},~{}~{}n^{\dagger}_{{\mathbf{k}},\nu}=n_{-{\mathbf{k}},\nu}~{}~{}(\nu,\nu^{\prime}=\pm).$
(17)

We note that the ${\mathbf{k}}={\mathbf{0}}$ component $n_{{\mathbf{0}},\pm}$
of the densities is related to the atom numbers as
$n_{{\mathbf{0}},\pm}=(N_{\uparrow}\pm N_{\downarrow})/(2\sqrt{A})$. The
Hamiltonian $H$ is then expressed as

$\begin{split}H=&\frac{n}{2}\sum_{{\mathbf{k}}}\begin{pmatrix}\theta_{-{\mathbf{k}},+}&\theta_{-{\mathbf{k}},-}\end{pmatrix}\begin{pmatrix}\Gamma_{\pm}({\mathbf{k}})&\pm\mathrm{i}\Gamma({\mathbf{k}})\\\
\mp\mathrm{i}\Gamma({\mathbf{k}})&\Gamma_{\mp}({\mathbf{k}})\end{pmatrix}\begin{pmatrix}\theta_{{\mathbf{k}},+}\\\
\theta_{{\mathbf{k}},-}\end{pmatrix}+\frac{1}{2n}\sum_{{\mathbf{k}}}\sum_{\nu=\pm}e_{{\mathbf{k}},\nu}n_{-{\mathbf{k}},\nu}n_{{\mathbf{k}},\nu},\end{split}$
(18)

where222We slightly change the definitions of $\Gamma_{\pm}({\mathbf{k}})$ and
$\Gamma({\mathbf{k}})$ from Ref. [46] by dividing them by $n$ so that they
have the dimension of energy.

$\displaystyle e_{{\mathbf{k}},\nu}$
$\displaystyle:=2\bar{g}_{\nu}n+\frac{\hbar^{2}{\mathbf{k}}^{2}}{2M},$ (19a)
$\displaystyle\Gamma_{+}({\mathbf{k}})$
$\displaystyle:=gn\ell^{4}\left[C_{1}\left(2k_{x}k_{y}\right)^{2}+C_{2}\left(k_{x}^{2}-k_{y}^{2}\right)^{2}-C_{3}(2k_{x}k_{y})\left(k_{x}^{2}-k_{y}^{2}\right)\right],$
(19b) $\displaystyle\Gamma_{-}({\mathbf{k}})$
$\displaystyle:=gn\ell^{2}\left(D_{1}k_{y}^{2}+D_{2}k_{x}^{2}-D_{3}k_{x}k_{y}\right),$
(19c) $\displaystyle\Gamma({\mathbf{k}})$
$\displaystyle:=\frac{1}{2}gn\ell^{3}F_{1}\left[(2k_{x}k_{y})k_{x}+(k_{x}^{2}-k_{y}^{2})k_{y}\right].$
(19d)

Here, the upper and lower signs correspond to the cases of parallel and
antiparallel fields, respectively.333The same sign rule applies to Eqs. (23),
(24), (27), (32), (34), (35), (90), (91), (43), (44), (51), and (56) below.

It is useful to decompose the Hamiltonian (18) into the zero-mode
(${\mathbf{k}}={\mathbf{0}}$) and oscillator-mode
(${\mathbf{k}}\neq{\mathbf{0}}$) parts. First, the zero-mode part is given by

$H^{\mathrm{zero}}=\sum_{\nu=\pm}\bar{g}_{\nu}n_{{\mathbf{0}},\nu}^{2}=\frac{\bar{g}_{+}}{4A}N^{2}+\frac{\bar{g}_{-}}{4A}(N_{\uparrow}-N_{\downarrow})^{2}.$
(20)

Thus, the zero-mode energy is specified by the atom numbers $N_{\uparrow}$ and
$N_{\downarrow}$. In our setting of balanced population
$N_{\uparrow}=N_{\downarrow}$, the zero-mode state is given by the product
state
$\left|N_{\uparrow}=N/2\right\rangle\left|N_{\downarrow}=N/2\right\rangle$,
which has no intercomponent entanglement.

Next, we discuss the oscillator-mode part $H^{\mathrm{osc}}$ of the
Hamiltonian (18). To treat this part, we perform canonical transformations in
two steps. The first transformation reads

$\tilde{\theta}_{{\mathbf{k}},+}=r_{\mathbf{k}}^{-1}\theta_{{\mathbf{k}},+},~{}\tilde{\theta}_{{\mathbf{k}},-}=r_{\mathbf{k}}\theta_{{\mathbf{k}},-},~{}\tilde{n}_{{\mathbf{k}},+}=r_{\mathbf{k}}n_{{\mathbf{k}},+},~{}\tilde{n}_{{\mathbf{k}},-}=r_{\mathbf{k}}^{-1}n_{{\mathbf{k}},-},$
(21)

where $r_{\mathbf{k}}:=(e_{{\mathbf{k}},+}/e_{{\mathbf{k}},-})^{1/4}$. Then,
$H^{\mathrm{osc}}$ is rewritten as

$H^{\mathrm{osc}}=\frac{n}{2}\sum_{{\mathbf{k}}\neq{\mathbf{0}}}\begin{pmatrix}\tilde{\theta}_{-{\mathbf{k}},+}&\tilde{\theta}_{-{\mathbf{k}},-}\end{pmatrix}M({\mathbf{k}})\begin{pmatrix}\tilde{\theta}_{{\mathbf{k}},+}\\\
\tilde{\theta}_{{\mathbf{k}},-}\end{pmatrix}+\frac{1}{2n}\sum_{{\mathbf{k}}\neq{\mathbf{0}}}\sum_{\nu=\pm}e_{\mathbf{k}}\tilde{n}_{-{\mathbf{k}},\nu}\tilde{n}_{{\mathbf{k}},\nu},$
(22)

where $e_{\mathbf{k}}:=\sqrt{e_{{\mathbf{k}},+}e_{{\mathbf{k}},-}}$. Here, the
$2\times 2$ matrix $M({\mathbf{k}})$ is given by

$M({\mathbf{k}}):=\begin{pmatrix}r_{\mathbf{k}}^{2}\Gamma_{\pm}({\mathbf{k}})&\pm\mathrm{i}\Gamma({\mathbf{k}})\\\
\mp\mathrm{i}\Gamma({\mathbf{k}})&r_{\mathbf{k}}^{-2}\Gamma_{\mp}({\mathbf{k}})\end{pmatrix}=\Gamma_{0}({\mathbf{k}})I\mp\Gamma({\mathbf{k}})\sigma_{y}+\Gamma_{z}({\mathbf{k}})\sigma_{z},$
(23)

where $I$ is the identity matrix, $(\sigma_{x},\sigma_{y},\sigma_{z})$ are the
Pauli matrices, and

$\Gamma_{0}({\mathbf{k}}):=\frac{1}{2}\left[r_{\mathbf{k}}^{2}\Gamma_{\pm}({\mathbf{k}})+r_{\mathbf{k}}^{-2}\Gamma_{\mp}({\mathbf{k}})\right],~{}~{}\Gamma_{z}({\mathbf{k}}):=\frac{1}{2}\left[r_{\mathbf{k}}^{2}\Gamma_{\pm}({\mathbf{k}})-r_{\mathbf{k}}^{-2}\Gamma_{\mp}({\mathbf{k}})\right].$
(24)

We then perform the second canonical transformation using the unitary matrix
$U({\mathbf{k}})$ as

$\begin{pmatrix}\tilde{\theta}_{{\mathbf{k}},+}\\\
\tilde{\theta}_{{\mathbf{k}},-}\end{pmatrix}=U({\mathbf{k}})\begin{pmatrix}\bar{\theta}_{{\mathbf{k}},1}\\\
\bar{\theta}_{{\mathbf{k}},2}\end{pmatrix},~{}~{}\begin{pmatrix}\tilde{n}_{{\mathbf{k}},+}\\\
\tilde{n}_{{\mathbf{k}},-}\end{pmatrix}=U({\mathbf{k}})\begin{pmatrix}\bar{n}_{{\mathbf{k}},1}\\\
\bar{n}_{{\mathbf{k}},2}\end{pmatrix}.$ (25)

We note that the second term of Eq. (22) is invariant under this
transformation if $U(-{\mathbf{k}})^{t}U({\mathbf{k}})=I$ for all
${\mathbf{k}}\neq{\mathbf{0}}$. It is therefore useful to choose
$U({\mathbf{k}})$ in such a way as to diagonalize the Hermitian matrix
$M({\mathbf{k}})$ as

$U^{-1}({\mathbf{k}})M({\mathbf{k}})U({\mathbf{k}})=\begin{pmatrix}m_{{\mathbf{k}},1}&0\\\
0&m_{{\mathbf{k}},2}\end{pmatrix},~{}~{}U({\mathbf{k}})={\mathrm{e}}^{\mathrm{i}\chi_{\mathbf{k}}\sigma_{x}/2}=\begin{pmatrix}\cos(\chi_{\mathbf{k}}/2)&\mathrm{i}\sin(\chi_{\mathbf{k}}/2)\\\
\mathrm{i}\sin(\chi_{\mathbf{k}}/2)&\cos(\chi_{\mathbf{k}}/2)\end{pmatrix},$
(26)

where444 Since $\Gamma_{z}({\mathbf{k}})=\Gamma_{z}(-{\mathbf{k}})$ and
$\Gamma(-{\mathbf{k}})=-\Gamma({\mathbf{k}})$, we find
$\chi_{-{\mathbf{k}}}=-\chi_{{\mathbf{k}}}$ and thus
$U(-{\mathbf{k}})^{t}=U({\mathbf{k}})^{\dagger}=U({\mathbf{k}})^{-1}$;
therefore, the aforementioned condition
$U(-{\mathbf{k}})^{t}U({\mathbf{k}})=I$ is met.

$\displaystyle
m_{{\mathbf{k}},j}:=\Gamma_{0}({\mathbf{k}})+(-1)^{j-1}\Lambda({\mathbf{k}})~{}~{}(j=1,2),$
(27a)
$\displaystyle\Lambda({\mathbf{k}})(\cos\chi_{\mathbf{k}},\sin\chi_{\mathbf{k}}):=\left(\Gamma_{z}({\mathbf{k}}),\Gamma({\mathbf{k}})\right).$
(27b)

In terms of the new set of canonical variables, the Hamiltonian (22) is
further rewritten as

$H^{\mathrm{osc}}=\frac{1}{2}\sum_{{\mathbf{k}}\neq{\mathbf{0}}}\sum_{j=1,2}\left(nm_{{\mathbf{k}},j}\bar{\theta}_{-{\mathbf{k}},j}\bar{\theta}_{{\mathbf{k}},j}+\frac{e_{\mathbf{k}}}{n}\bar{n}_{-{\mathbf{k}},j}\bar{n}_{{\mathbf{k}},j}\right).$
(28)

Here, $\bar{\theta}_{{\mathbf{k}},j}$ and $\bar{n}_{{\mathbf{k}},j}$ satisfy
$[\bar{\theta}_{{\mathbf{k}},j},\bar{n}_{-{\mathbf{k}}^{\prime},j^{\prime}}]=\mathrm{i}\delta_{jj^{\prime}}\delta_{{\mathbf{k}}{\mathbf{k}}^{\prime}}~{}({\mathbf{k}},{\mathbf{k}}^{\prime}\neq{\mathbf{0}};~{}j,j^{\prime}=1,2)$
as the transformations in Eqs. (21) and (25) leave the commutation relations
unchanged. Finally, introducing the (bogolon) annihilation and creation
operators

$\gamma_{{\mathbf{k}},j}:=\frac{1}{\sqrt{2}}\left(\sqrt{n}\zeta_{{\mathbf{k}},j}\bar{\theta}_{{\mathbf{k}},j}+\frac{\mathrm{i}\bar{n}_{{\mathbf{k}},j}}{\sqrt{n}\zeta_{{\mathbf{k}},j}}\right),~{}\gamma_{{\mathbf{k}},j}^{\dagger}:=\frac{1}{\sqrt{2}}\left(\sqrt{n}\zeta_{{\mathbf{k}},j}\bar{\theta}_{-{\mathbf{k}},j}-\frac{\mathrm{i}\bar{n}_{-{\mathbf{k}},j}}{\sqrt{n}\zeta_{{\mathbf{k}},j}}\right)~{}({\mathbf{k}}\neq{\mathbf{0}};~{}j=1,2),$
(29)

with $\zeta_{{\mathbf{k}},j}:=(m_{{\mathbf{k}},j}/e_{\mathbf{k}})^{1/4}$, we
can diagonalize the Hamiltonian (28) as

$H^{\mathrm{osc}}=\sum_{{\mathbf{k}}\neq{\mathbf{0}}}\sum_{j=1,2}E_{j}({\mathbf{k}})\left(\gamma_{{\mathbf{k}},j}^{\dagger}\gamma_{{\mathbf{k}},j}+\frac{1}{2}\right),~{}~{}E_{j}({\mathbf{k}}):=\sqrt{m_{{\mathbf{k}},j}e_{\mathbf{k}}}.$
(30)

The ground state $|0^{\mathrm{osc}}\rangle$ of this Hamiltonian is specified
by the condition that $\gamma_{{\mathbf{k}},\pm}|0^{\mathrm{osc}}\rangle=0$
for all ${\mathbf{k}}\neq{\mathbf{0}}$.

We now discuss the single-particle spectrum $E_{j}({\mathbf{k}})~{}(j=1,2)$ in
the long-wavelength limit $k\ell\ll 1$. In this limit, we can make the
approximation $r_{\mathbf{k}}\approx
r_{\mathbf{0}}=\left(\bar{g}_{+}/\bar{g}_{-}\right)^{1/4}$ and
$e_{\mathbf{k}}\approx e_{\mathbf{0}}=2(\bar{g}_{+}\bar{g}_{-})^{1/2}n$. Since
$\Gamma_{+}({\mathbf{k}})=\mathcal{O}(k^{4})$,
$\Gamma_{-}({\mathbf{k}})=\mathcal{O}(k^{2})$, and
$\Gamma({\mathbf{k}})=\mathcal{O}(k^{3})$ as seen in Eq. (19), we can
approximate $m_{{\mathbf{k}},j}$ in Eq. (27) as

$\displaystyle m_{{\mathbf{k}},1}\approx r_{\mathbf{0}}^{\mp
2}\Gamma_{-}({\mathbf{k}}),~{}~{}m_{{\mathbf{k}},2}\approx r_{\mathbf{0}}^{\pm
2}\left[\Gamma_{+}({\mathbf{k}})-\frac{\Gamma({\mathbf{k}})^{2}}{\Gamma_{-}({\mathbf{k}})}\right].$
(31)

By parametrizing the wave vector as
$(k_{x},k_{y})=k(\cos\varphi,\sin\varphi)$, $m_{{\mathbf{k}},j}$ can also be
expressed as

$m_{{\mathbf{k}},1}\approx r_{\mathbf{0}}^{\mp
2}gnD(\varphi)(k\ell)^{2},~{}~{}m_{{\mathbf{k}},2}\approx r_{\mathbf{0}}^{\pm
2}gnC(\varphi)(k\ell)^{4},$ (32)

where

$\displaystyle C(\varphi)$
$\displaystyle:=C_{1}\sin^{2}(2\varphi)+C_{2}\cos^{2}(2\varphi)-C_{3}\sin(2\varphi)\cos(2\varphi)-C_{4}\sin^{2}(3\varphi),$
(33a) $\displaystyle D(\varphi)$
$\displaystyle:=D_{1}\sin^{2}(\varphi)+D_{2}\cos^{2}(\varphi)-D_{3}\sin(\varphi)\cos(\varphi),$
(33b)

with $C_{4}:=F_{1}^{2}/4D_{1}$. We then obtain the low-energy spectrum as

$E_{j}({\mathbf{k}})\approx\sqrt{m_{{\mathbf{k}},j}e_{\mathbf{0}}}\approx\sqrt{2}gn(k\ell)^{j}f_{j}(\varphi)~{}~{}~{}(j=1,2),$
(34)

where the dependence on the angle $\varphi$ is expressed by the dimensionless
functions

$f_{1}(\varphi)=\sqrt{\bar{g}_{\mp}D(\varphi)/g},~{}~{}~{}f_{2}(\varphi)=\sqrt{\bar{g}_{\pm}C(\varphi)/g}.$
(35)

We thus find that low-energy modes with linear and quadratic dispersion
relations emerge with anisotropy that depends on the lattice structure.
Furthermore, the low-energy dispersion relations for parallel (P) and
antiparallel (AP) fields are related by proper rescaling as follows:

$E_{1}^{\mathrm{P}}({\mathbf{k}})/\sqrt{\bar{g}_{-}}=E_{1}^{\mathrm{AP}}({\mathbf{k}})/\sqrt{\bar{g}_{+}},~{}~{}~{}E_{2}^{\mathrm{P}}({\mathbf{k}})/\sqrt{\bar{g}_{+}}=E_{2}^{\mathrm{AP}}({\mathbf{k}})/\sqrt{\bar{g}_{-}}.$
(36)

Using the dimensionless functions $f_{j}(\varphi)~{}(j=1,2)$ for the two types
of fields, these relations can also be written as

$\begin{split}f_{1}^{\mathrm{P}}(\varphi)\sqrt{\frac{g}{\bar{g}_{-}}}=f_{1}^{\mathrm{AP}}(\varphi)\sqrt{\frac{g}{\bar{g}_{+}}}=\sqrt{D(\varphi)},~{}~{}~{}f_{2}^{\mathrm{P}}(\varphi)\sqrt{\frac{g}{\bar{g}_{+}}}=f_{2}^{\mathrm{AP}}(\varphi)\sqrt{\frac{g}{\bar{g}_{-}}}=\sqrt{C(\varphi)}.\end{split}$
(37)

In Sec. 4.2, we will confirm these relations through the numerical
calculations by the Bogoliubov theory for all the five vortex-lattice
structures in Fig. 1(a)-(e). While similar rescaling relations are also
discussed in Ref. [46], the importance of using the renormalized coupling
constants $\bar{g}_{\pm}$ is overlooked there. Without the renormalization,
the rescaling relations in Eqs. (36) and (37) are satisfied only for
overlapping vortex lattices where $\beta=\beta_{\uparrow\downarrow}$, as
confirmed numerically in Ref. [46].

###  2.5 Intercomponent entanglement

We now calculate the reduced density matrix (RDM) $\rho_{\uparrow}$ for the
spin-$\uparrow$ component, which is defined by starting from the ground state
$|0^{\mathrm{zero}}\rangle\otimes|0^{\mathrm{osc}}\rangle$ of the total system
and trancing out the degrees of freedom in the spin-$\downarrow$ component. We
then discuss the properties of the intercomponent entanglement. Because of the
decoupling of the zero and oscillator modes, the RDM takes the form of
$\rho_{\uparrow}=\rho_{\uparrow}^{\mathrm{zero}}\otimes\rho_{\uparrow}^{\mathrm{osc}}$.
As the zero-mode ground state
$\left|N_{\uparrow}=N/2\right\rangle\left|N_{\downarrow}=N/2\right\rangle$ is
a product state, there is no intercomponent entanglement in the zero-mode
part; the RDM in this part is given simply by
$\rho_{\uparrow}^{\mathrm{zero}}=|N_{\uparrow}=N/2\rangle\langle
N_{\uparrow}=N/2|$. Below we consider the oscillator-mode part
$\rho_{\uparrow}^{\mathrm{osc}}$.

For $\rho_{\uparrow}^{\mathrm{osc}}$, we introduce the following Gaussian
ansatz [61, 62, 63, 67, 68, 69]:

$\rho^{{\mathrm{osc}}}_{\uparrow}=\frac{1}{Z^{{\mathrm{osc}}}_{\mathrm{e}}}\mathrm{e}^{-H_{\mathrm{e}}^{\mathrm{osc}}},~{}Z^{\mathrm{osc}}_{\mathrm{e}}=\Tr\mathrm{e}^{-H^{{\mathrm{osc}}}_{\mathrm{e}}},~{}H^{\mathrm{osc}}_{\mathrm{e}}=\frac{1}{2}\sum_{{\mathbf{k}}\neq{\mathbf{0}}}\left(nF_{\mathbf{k}}\theta_{-{\mathbf{k}},\uparrow}\theta_{{\mathbf{k}},\uparrow}+\frac{G_{\mathbf{k}}}{n}n_{-{\mathbf{k}},\uparrow}n_{{\mathbf{k}},\uparrow}\right),$
(38)

where $F_{\mathbf{k}}$ and $G_{\mathbf{k}}$ are positive dimensionless
coefficients to be determined later and we assume
$F_{\mathbf{k}}=F_{-{\mathbf{k}}}$ and $G_{\mathbf{k}}=G_{-{\mathbf{k}}}$ for
convenience. By introducing annihilation and creation operators as

$\displaystyle\eta_{\mathbf{k}}=$
$\displaystyle\frac{1}{\sqrt{2}}\left[\sqrt{n}\left(\frac{F_{\mathbf{k}}}{G_{\mathbf{k}}}\right)^{1/4}\theta_{{\mathbf{k}},\uparrow}+\frac{\mathrm{i}}{\sqrt{n}}\left(\frac{G_{\mathbf{k}}}{F_{\mathbf{k}}}\right)^{1/4}n_{{\mathbf{k}},\uparrow}\right],$
(39a) $\displaystyle\eta_{\mathbf{k}}^{\dagger}=$
$\displaystyle\frac{1}{\sqrt{2}}\left[\sqrt{n}\left(\frac{F_{\mathbf{k}}}{G_{\mathbf{k}}}\right)^{1/4}\theta_{-{\mathbf{k}},\uparrow}-\frac{\mathrm{i}}{\sqrt{n}}\left(\frac{G_{\mathbf{k}}}{F_{\mathbf{k}}}\right)^{1/4}n_{-{\mathbf{k}},\uparrow}\right]~{}~{}~{}({\mathbf{k}}\neq{\mathbf{0}}),$
(39b)

the entanglement Hamiltonian $H_{\mathrm{e}}^{\mathrm{osc}}$ in Eq. (38) is
diagonalized as

$H^{{\mathrm{osc}}}_{\mathrm{e}}=\sum_{{\mathbf{k}}\neq{\mathbf{0}}}\xi_{\mathbf{k}}\left(\eta_{\mathbf{k}}^{\dagger}\eta_{\mathbf{k}}+\frac{1}{2}\right),$
(40)

where $\xi_{\mathbf{k}}:=\sqrt{F_{\mathbf{k}}G_{\mathbf{k}}}$ is the single-
particle ES.

The single-particle ES $\xi_{\mathbf{k}}$ and the coefficients
$F_{\mathbf{k}}$ and $G_{\mathbf{k}}$ can be determined in the following way
[61]. Using the relations in Eq. (39) and the Bose distribution function

$\Tr\left(\eta_{\mathbf{k}}^{\dagger}\eta_{\mathbf{k}}\rho^{\mathrm{osc}}_{\uparrow}\right)=\frac{1}{{\mathrm{e}}^{\xi_{\mathbf{k}}}-1}\equiv
f_{\mathrm{B}}(\xi_{\mathbf{k}}),$ (41)

we obtain the phase and density correlators as

$\displaystyle\Tr\left(\theta_{-{\mathbf{k}},\uparrow}\theta_{{\mathbf{k}},\uparrow}\rho^{\mathrm{osc}}_{\uparrow}\right)=\frac{1}{n}\left(\frac{G_{\mathbf{k}}}{F_{\mathbf{k}}}\right)^{1/2}\\!\left[f_{\mathrm{B}}(\xi_{\mathbf{k}})+\frac{1}{2}\right],$
(42a)
$\displaystyle\Tr\left(n_{-{\mathbf{k}},\uparrow}n_{{\mathbf{k}},\uparrow}\rho^{\mathrm{osc}}_{\uparrow}\right)=n\left(\frac{F_{\mathbf{k}}}{G_{\mathbf{k}}}\right)^{1/2}\\!\left[f_{\mathrm{B}}(\xi_{\mathbf{k}})+\frac{1}{2}\right].$
(42b)

We can then determine $f_{\mathrm{B}}(\xi_{\mathbf{k}})$ and
$F_{\mathbf{k}}/G_{\mathbf{k}}$ by requiring these correlators to be equal to
the same correlators calculated for the oscillator ground state
$|0^{\mathrm{osc}}\rangle$ of the total system. Details of this calculation
are described in A. In the long-wavelength limit $k\ell\ll 1$, we obtain

$\xi_{\mathbf{k}}\approx c(\varphi)\sqrt{k\ell},~{}~{}F_{\mathbf{k}}\approx
F(\varphi)(k\ell)^{2},~{}~{}G_{\mathbf{k}}\approx\frac{G(\varphi)}{k\ell},$
(43)

where the dependences on the angle $\varphi$ are expressed by the
dimensionless functions

$c(\varphi)=4\left[\frac{\bar{g}_{\mp}C(\varphi)}{\bar{g}_{\pm}D(\varphi)}\right]^{1/4},~{}F(\varphi)=4\sqrt{\frac{2gC(\varphi)}{\bar{g}_{\pm}}},~{}G(\varphi)=2\sqrt{\frac{2\bar{g}_{\mp}}{gD(\varphi)}}.$
(44)

We find that the ES shows a gapless square-root dispersion relation with
anisotropy that depends on the lattice structure. Furthermore, similarly to
the case of excitation spectra [see Eq. (36)], the single-particle ES
$\xi_{\mathbf{k}}$ for parallel (P) and antiparallel (AP) fields are related
by suitable rescaling as

$\left(\frac{\bar{g}_{+}}{\bar{g}_{-}}\right)^{1/4}\xi_{\mathbf{k}}^{\mathrm{P}}=\left(\frac{\bar{g}_{-}}{\bar{g}_{+}}\right)^{1/4}\xi_{\mathbf{k}}^{\mathrm{AP}}.$
(45)

Using the dimensionless functions $c(\varphi)$ for the two types of fields,
this relation can also be written as

$\left(\frac{\bar{g}_{+}}{\bar{g}_{-}}\right)^{1/4}c^{\mathrm{P}}(\varphi)=\left(\frac{\bar{g}_{-}}{\bar{g}_{+}}\right)^{1/4}c^{\mathrm{AP}}(\varphi)=4\left[\frac{C(\varphi)}{D(\varphi)}\right]^{1/4}.$
(46)

The entanglement Hamiltonian is given in the long-wavelength limit by
$H_{\mathrm{e}}^{\mathrm{osc}}$ in Eq. (38) with $F_{\mathbf{k}}$ and
$G_{\mathbf{k}}$ in Eq. (43). Using the fields
$\theta_{\uparrow}({\mathbf{r}})$ and $n_{\uparrow}({\mathbf{r}})$ in real
space, it can be expressed as

$\begin{split}H_{\mathrm{e}}=\int\mathrm{d}^{2}{\mathbf{r}}\int\mathrm{d}^{2}{\mathbf{r}}^{\prime}\Biggl{[}\frac{n\ell^{2}}{2}U_{F}({\mathbf{r}}-{\mathbf{r}}^{\prime})\nabla\theta_{\uparrow}({\mathbf{r}})\cdot\nabla\theta_{\uparrow}({\mathbf{r}}^{\prime})+\frac{1}{2n}U_{G}({\mathbf{r}}-{\mathbf{r}}^{\prime})n_{\uparrow}({\mathbf{r}})n_{\uparrow}({\mathbf{r}}^{\prime})\Biggr{]},\end{split}$
(47)

where we introduce the interaction potentials

$\begin{split}U_{F}({\mathbf{r}}-{\mathbf{r}}^{\prime})=\frac{1}{A}\sum_{{\mathbf{k}}}F(\varphi){\mathrm{e}}^{\mathrm{i}{\mathbf{k}}\cdot({\mathbf{r}}-{\mathbf{r}}^{\prime})},~{}~{}U_{G}({\mathbf{r}}-{\mathbf{r}}^{\prime})=\lim_{\alpha\to
0^{+}}\frac{1}{A}\sum_{{\mathbf{k}}\neq{\mathbf{0}}}\frac{G(\varphi)}{k\ell}{\mathrm{e}}^{-\alpha
k+\mathrm{i}{\mathbf{k}}\cdot({\mathbf{r}}-{\mathbf{r}}^{\prime})}.\end{split}$
(48)

Here, we use the convergence factor ${\mathrm{e}}^{-\alpha k}$ to regularize
the infinite sum for $U_{G}({\mathbf{r}}-{\mathbf{r}}^{\prime})$. For
simplicity, we consider the case of overlapping triangular lattices, in which
$C(\varphi)$ and $D(\varphi)$ [and thus $F(\varphi)$ and $G(\varphi)$ as well]
are constant; see (a) in Eq. (15). In this case, the potentials in Eq. (48)
are calculated as

$U_{F}({\mathbf{r}}-{\mathbf{r}}^{\prime})=F\delta({\mathbf{r}}-{\mathbf{r}}^{\prime}),~{}~{}U_{G}({\mathbf{r}}-{\mathbf{r}}^{\prime})=\frac{G}{2\pi\ell|{\mathbf{r}}-{\mathbf{r}}^{\prime}|}.$
(49)

For the calculation of $U_{G}({\mathbf{r}}-{\mathbf{r}}^{\prime})$, we refer
the reader to Appendix A of Ref. [61]. We note that
$\theta_{\uparrow}({\mathbf{r}})$ is the regular part of the superfluid phase
of the spin-$\uparrow$ component, and that its gradient is related to the
regular part of the superfluid velocity,
${\mathbf{v}}_{\mathrm{s},\uparrow}({\mathbf{r}})=-\frac{\hbar}{M}\nabla\theta_{\uparrow}({\mathbf{r}})$.
Therefore, the entanglement Hamiltonian (47) has a short-range interaction in
terms of the superfluid velocity ${\mathbf{v}}_{s,\uparrow}({\mathbf{r}})$ and
a long-range one in terms of the density $n_{\uparrow}({\mathbf{r}})$. If the
density interaction were short-ranged, the ES would show a phonon mode with a
linear dispersion relation. Therefore, the anomalous square-root dispersion
relation in Eq. (43) is closely related with the presence of a long-range
interaction in $H_{\mathrm{e}}$.

Using the single-particle ES $\xi_{\mathbf{k}}$ in Eq. (43), we can calculate
the intercomponent EE $S_{\mathrm{e}}$. For simplicity, we assume that
$\xi_{\mathbf{k}}$ is isotropic, i.e., $c(\varphi)$ is constant, as in the
case of overlapping triangular lattices; however, we expect that the result
holds qualitatively for all the lattice structures. With this simplification,
the EE is calculated as (see Appendix B of Ref. [61] for the derivation)

$S_{\mathrm{e}}=\frac{\sigma
A}{c^{4}\ell^{2}}-\frac{1}{2}\ln\frac{\sqrt{A}}{2\pi
c^{2}\ell}+O(1)=\frac{2\pi\sigma
N_{\mathrm{v}}}{c^{4}}-\frac{1}{4}\ln\frac{N_{\mathrm{v}}}{2\pi c^{4}}+O(1).$
(50)

Here, the leading contribution is given by the first term, which is
proportional to the area $A$ with a non-universal coefficient $\sigma$ that
depends, e.g., on the choice of the high-momentum cutoff. Besides, there is a
subleading logarithmic term with the universal coefficient (equal to $-1/2$
when written as a function of the linear system size $\sqrt{A}$), which is
identified through a careful examination of small-$k$ contributions and
therefore originates from the Nambu-Goldstone modes. The intercomponent EE per
flux in the thermodynamic limit is calculated as

$\lim_{N_{\mathrm{v}}\to\infty}\frac{S_{\mathrm{e}}}{N_{\mathrm{v}}}=\frac{2\pi\sigma}{c^{4}}=\frac{\pi\sigma\bar{g}_{\pm}D}{128\bar{g}_{\mp}C}.$
(51)

Because of the factor $\bar{g}_{\pm}/\bar{g}_{\mp}$ in Eq. (51), when the
intercomponent interaction $g_{\uparrow\downarrow}$ is repulsive (attractive),
the intercomponent EE is expected to be larger for the case of parallel
(antiparallel) fields. We note that the above calculation of the
intercomponent EE $S_{\mathrm{e}}$ is simple yet approximate as it is based on
the single-particle ES $\xi_{\mathbf{k}}$ for long wavelengths. In Sec. 4.3,
we will present numerical results on $S_{\mathrm{e}}$ based on the Bogoliubov
theory, in which $\xi_{\mathbf{k}}$ over the full Brillouin zone is taken into
account, and confirm the consistency with the field-theoretical predictions.

###  2.6 Intracomponent correlation functions

Here we calculate some intracomponent correlation functions, and discuss their
connections with the (long-wavelength) entanglement Hamiltonian
$H_{\mathrm{e}}$ obtained in the preceding section. Let $\langle{\cal
O}\rangle$ denote the expectation value of an operator ${\cal O}$ with respect
to the ground state $|0^{\mathrm{zero}}\rangle\otimes|0^{\mathrm{osc}}\rangle$
of the total system. If ${\cal O}$ acts only on the spin-$\uparrow$ component,
$\langle{\cal O}\rangle$ should be equal to $\Tr\left({\cal
O}{\mathrm{e}}^{-H_{\mathrm{e}}}\right)/\Tr{\mathrm{e}}^{-H_{\mathrm{e}}}$ as
far as long-distance properties are concerned. Our purpose here is to
investigate how the unusual long-range interactions in $H_{\mathrm{e}}$
manifest themselves in the correlation properties of the system.

Owing to the gapless ES $\xi_{\mathbf{k}}$, we can approximate the Bose
distribution function (41) as
$f_{\mathrm{B}}(\xi_{\mathbf{k}})\approx\xi_{\mathbf{k}}^{-1}=(F_{\mathbf{k}}G_{\mathbf{k}})^{-1/2}$
for sufficiently small $k$. Then, in the long-wavelength limit, Eq. (42) gives

$\displaystyle\langle\theta_{-{\mathbf{k}},\uparrow}\theta_{{\mathbf{k}},\uparrow}\rangle\approx\frac{1}{nF_{\mathbf{k}}}\approx\frac{1}{nF(\varphi)k^{2}\ell^{2}},~{}~{}~{}\langle
n_{-{\mathbf{k}},\uparrow}n_{{\mathbf{k}},\uparrow}\rangle\approx\frac{n}{G_{\mathbf{k}}}\approx\frac{nk\ell}{G(\varphi)}~{}~{}({\mathbf{k}}\neq{\mathbf{0}}),$
(52)

where we use Eq. (43). Therefore, the phase and density fluctuations are
directly related to the coefficients $F_{\mathbf{k}}$ and $G_{\mathbf{k}}$,
respectively, in the entanglement Hamiltonian (38). Equation (52) indicates
that in the long-wavelength limit $k\to 0$, the phase fluctuation diverges and
the density fluctuation is suppressed. From the viewpoint of the entanglement
Hamiltonian, suppression of the density fluctuation is a consequence of the
long-range interaction in terms of the density. The enhanced phase fluctuation
in the long-wavelength limit leads to a quasi-long-range order in the one-
particle density matrix, as we explain in the following.

The one-particle density matrix plays a key role in the characterization of
the Bose-Einstein condensation [70, 71, 72]. To analyze its behavior in our
course-grained description, we introduce the modified bosonic field
$\tilde{\psi}_{\uparrow}=e^{-i\theta_{\uparrow}}\sqrt{n_{\uparrow}}$. As
$\theta_{\uparrow}({\mathbf{r}})$ is the regular part of the superfluid phase
and $n_{\uparrow}({\mathbf{r}})$ is the course-grained density,
$\tilde{\psi}({\mathbf{r}})$ is expected to vary slowly over space. We
consider the modified one-particle density matrix

$\langle\tilde{\psi}_{\uparrow}({\mathbf{r}})^{\dagger}\tilde{\psi}_{\uparrow}({\mathbf{0}})\rangle=\left\langle\sqrt{n_{\uparrow}({\mathbf{r}})}{\mathrm{e}}^{i(\theta_{\uparrow}({\mathbf{r}})-\theta_{\uparrow}({\mathbf{0}}))}\sqrt{n_{\uparrow}({\mathbf{0}})}\right\rangle,$
(53)

which describes the slowly varying component of the ordinary one-particle
density matrix. Its long-distance behavior is determined dominantly by the
phase fluctuation as seen in the small-$k$ behavior of Eq. (52). Here we again
consider the case of overlapping triangular lattices where $F(\varphi)$ is
constant. By using Eq. (52), the phase correlation function in real space is
obtained as (see B for the derivation)

$\left\langle\left[\theta_{\uparrow}({\mathbf{r}})-\theta_{\uparrow}({\mathbf{0}})\right]^{2}\right\rangle=\frac{2}{A}\sum_{{\mathbf{k}}\neq{\mathbf{0}}}{\mathrm{e}}^{-\alpha
k}\left[1-\cos\left({{\mathbf{k}}\cdot{\mathbf{r}}}\right)\right]\langle\theta_{-{\mathbf{k}},\uparrow}\theta_{{\mathbf{k}},\uparrow}\rangle\approx\frac{1}{\pi
nF\ell^{2}}\ln\frac{r}{2\alpha}~{}~{}(r\gg\alpha),$ (54)

where we again introduce the convergence factor ${\mathrm{e}}^{-\alpha k}$ to
regularize the infinite sum. The modified one-particle density matrix is then
obtained as

$\langle\tilde{\psi}_{\uparrow}({\mathbf{r}})^{\dagger}\tilde{\psi}_{\uparrow}({\mathbf{0}})\rangle\approx
n\exp\left\\{-\frac{1}{2}\left\langle\left[\theta_{\uparrow}({\mathbf{r}})-\theta_{\uparrow}({\mathbf{0}})\right]^{2}\right\rangle\right\\}\approx
n\left(\frac{r}{2\alpha}\right)^{-\frac{1}{2\pi nF\ell^{2}}}.$ (55)

We thus have a quasi-long-range order in the one-particle density matrix.

For a finite system of area $A$, the particle density $n_{0}$ of the
condensate can be evaluated from the one-particle density matrix (55) at the
large separation $r=\sqrt{A}/2$. The density $n^{\prime}$ of the depletion is
then given by $n^{\prime}=n-n_{0}$. Assuming $n^{\prime}\ll n$ as required for
the Bogoliubov theory and the present effective theory, the fraction of
depletion is estimated as

$\frac{n^{\prime}}{n}\approx\frac{1}{2}\left\langle\left[\theta_{\uparrow}({\mathbf{r}})-\theta_{\uparrow}({\mathbf{0}})\right]^{2}\right\rangle\bigg{|}_{r=\sqrt{A}/2}\approx\frac{1}{2\pi
nF\ell^{2}}\ln\frac{\sqrt{A}}{4\alpha}=\frac{1}{2\nu}\left(\frac{\bar{g}_{\pm}}{2gC}\right)^{1/2}\ln\frac{\sqrt{A}}{4\alpha}.$
(56)

where we use Eq. (44) and $\nu=2nA/N_{\mathrm{v}}=4\pi n\ell^{2}$. For fixed
$\nu$, the fraction of quantum depletion of the condensate increases
logarithmically as a function of the vortex number
$N_{\mathrm{v}}=A/(2\pi\ell^{2})$. Furthermore, for an intercomponent
repulsion (attraction), it is larger for the case of parallel (antiparallel)
fields, indicating larger quantum fluctuations. This behavior is in accord
with the larger intercomponent EE given in Eq. (51). We note that a quasi-
long-range order in the one-particle density matrix and a logarithmic increase
of the fraction of depletion as in Eqs. (55) and (56) have also been discussed
for a vortex lattice in a scalar BEC [18, 20, 21, 22].

##  3 Bogoliubov theory

The Bogoliubov theory with the LLL approximation has been formulated for a
scalar BEC in Refs. [18, 21, 22]. In Ref. [46], we have applied this theory to
the present problem of binary BECs in parallel and antiparallel magnetic
fields. In Secs. 3.1 and 3.2, we summarize the formulation of Ref. [46]. In
particular, we give an expression of the ground-state energy [Eq. (67)] in
which a quantum correction due to zero-point fluctuations is included. In Sec.
3.3, we use this formulation to derive expressions of the intercomponent ES
and EE. Numerical results based on this formulation will be presented in Sec.
4.

###  3.1 LLL magnetic Bloch states and Hamiltonian

We employ the LLL magnetic Bloch states
$\\{\Psi_{{\mathbf{k}}\alpha}({\mathbf{r}})\\}$ [22, 73, 74, 75] as a
convenient single-particle basis for describing vortex lattices. The
expressions of these states are shown in Ref. [46], and we summarize their
main features in the following. Let $\mathbf{a}_{1}$ and $\mathbf{a}_{2}$ be
the primitive vectors of a vortex lattice, and let ${\mathbf{K}}_{\alpha}$ be
the pseudomomentum operator for a spin-$\alpha$ atom in a synthetic magnetic
field $B_{\alpha}~{}(\alpha=\uparrow,\downarrow)$. As expected for a “Bloch
state”, $\Psi_{{\mathbf{k}}\alpha}({\mathbf{r}})$ is an eigenstate of the
magnetic translation
${\mathrm{e}}^{-\mathrm{i}{\mathbf{K}}_{\alpha}\cdot{\mathbf{a}}_{j}/\hbar}$
with an eigenvalue
${\mathrm{e}}^{-i{\mathbf{k}}\cdot{\mathbf{a}}_{j}}~{}(j=1,2)$. By taking
$N_{\mathrm{v}}$ discrete wave vectors ${\mathbf{k}}$ consistent with the
boundary conditions of the system,
$\\{\Psi_{{\mathbf{k}}\alpha}({\mathbf{r}})\\}$ form a complete orthogonal
basis of the LLL manifold.555 In numerical calculations presented later, we
set
${\mathbf{k}}=\frac{n_{1}}{N_{\mathrm{v1}}}{\mathbf{b}}_{1}+\frac{n_{2}}{N_{\mathrm{v2}}}{\mathbf{b}}_{2}$
with $n_{j}\in\\{0,1,\dots,N_{\mathrm{v}j}-1\\}$ and
$N_{\mathrm{v1}}N_{\mathrm{v2}}=N_{\mathrm{v}}$. Here, ${\mathbf{b}}_{1}$ and
${\mathbf{b}}_{2}$ are the reciprocal primitive vectors as shown in Fig. 1.
Notably, $\Psi_{{\mathbf{k}}\alpha}({\mathbf{r}})$ has a periodic pattern of
zeros at [74]

${\mathbf{r}}=n_{1}\mathbf{a}_{1}+n_{2}\mathbf{a}_{2}+\frac{1}{2}({\mathbf{a}}_{1}+{\mathbf{a}}_{2})-\epsilon_{\alpha}\ell^{2}{\mathbf{e}}_{z}\times{\mathbf{k}},~{}~{}n_{1},n_{2}\in\mathbb{Z}.$
(57)

Therefore, $\Psi_{{\mathbf{k}}\alpha}({\mathbf{r}})$ represents a vortex
lattice with primitive vectors ${\mathbf{a}}_{1}$ and ${\mathbf{a}}_{2}$ for
any ${\mathbf{k}}$, and the locations of vortices (zeros) can be shifted by
varying ${\mathbf{k}}$. Vortex lattices of binary BECs in Fig. 1 are obtained
when spin-$\alpha$ bosons condense into
$\Psi_{{\mathbf{q}}_{\alpha},\alpha}({\mathbf{r}})$, where the wave vectors
${\mathbf{q}}_{\uparrow}$ and ${\mathbf{q}}_{\downarrow}$ are chosen in a way
consistent with the displacement $u_{1}{\mathbf{a}}_{1}+u_{2}{\mathbf{a}}_{2}$
between the components.

In the LLL approximation, the kinetic energy of each particle stays constant,
and therefore we can focus on the interaction Hamiltonian $H_{\mathrm{int}}$.
Using the LLL magnetic Bloch states
$\\{\Psi_{{\mathbf{k}}\alpha}({\mathbf{r}})\\}$ as the basis, it is
represented as

$H_{\mathrm{int}}=\frac{1}{2}\sum_{\alpha,\beta}\sum_{\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3},\mathbf{k}_{4}}V_{\alpha\beta}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3},\mathbf{k}_{4})b^{{\dagger}}_{\mathbf{k}_{1}\alpha}b^{{\dagger}}_{\mathbf{k}_{2}\beta}b_{\mathbf{k}_{3}\beta}b_{\mathbf{k}_{4}\alpha},$
(58)

where $b_{{\mathbf{k}}\alpha}$ is a bosonic annihilation operator for the
state $\Psi_{{\mathbf{k}}\alpha}({\mathbf{r}})$ and

$V_{\alpha\beta}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3},\mathbf{k}_{4})=g_{\alpha\beta}\int\mathrm{d}^{2}\mathbf{r}~{}\Psi^{\ast}_{\mathbf{k}_{1}\alpha}(\mathbf{r})\Psi^{\ast}_{\mathbf{k}_{2}\beta}(\mathbf{r})\Psi_{\mathbf{k}_{3}\beta}(\mathbf{r})\Psi_{\mathbf{k}_{4}\alpha}(\mathbf{r}).$
(59)

The expression of the interaction matrix element
$V_{\alpha\beta}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3},\mathbf{k}_{4})$
that is convenient for numerical calculations is given in Ref. [46].

###  3.2 Bogoliubov approximation

We now apply the Bogoliubov approximation [18, 21, 22, 70], assuming that the
condensation occurs at the wave vector ${\mathbf{q}}_{\alpha}$ in the
spin-$\alpha$ component. To this end, it is useful to introduce

$\tilde{b}_{{\mathbf{k}}\alpha}:=b_{{\mathbf{q}}_{\alpha}+{\mathbf{k}},\alpha},~{}~{}\tilde{V}_{\alpha\beta}({\mathbf{k}}_{1},{\mathbf{k}}_{2},{\mathbf{k}}_{3},{\mathbf{k}}_{4}):=V_{\alpha\beta}({\mathbf{q}}_{\alpha}+{\mathbf{k}}_{1},{\mathbf{q}}_{\beta}+{\mathbf{k}}_{2},{\mathbf{q}}_{\beta}+{\mathbf{k}}_{3},{\mathbf{q}}_{\alpha}+{\mathbf{k}}_{4}).$
(60)

By substituting

$\tilde{b}_{{\mathbf{0}}\alpha}\simeq\tilde{b}_{{\mathbf{0}}\alpha}^{\dagger}\simeq\sqrt{N_{\alpha}-\sum_{{\mathbf{k}}\neq{\mathbf{0}}}\tilde{b}_{{\mathbf{k}}\alpha}^{\dagger}\tilde{b}_{{\mathbf{k}}\alpha}}$
(61)

in Eq. (58) and retaining terms up to the second order in
$\tilde{b}_{{\mathbf{k}}\alpha}$ and
$\tilde{b}_{{\mathbf{k}}\alpha}^{\dagger}$ (${\mathbf{k}}\neq 0$), we obtain
the Bogoliubov Hamiltonian

$\begin{split}H_{\mathrm{int}}=&\frac{1}{2}\sum_{\alpha,\beta}N_{\alpha}N_{\beta}\tilde{V}_{\alpha\beta}({\mathbf{0}},{\mathbf{0}},{\mathbf{0}},{\mathbf{0}})-\frac{1}{2}\sum_{{\mathbf{k}}\neq{\mathbf{0}}}J({\mathbf{k}})\\\
&+\frac{1}{2}\sum_{{\mathbf{k}}\neq{\mathbf{0}}}\left(\tilde{b}_{{\mathbf{k}}\uparrow}^{\dagger},\tilde{b}_{{\mathbf{k}}\downarrow}^{\dagger},\tilde{b}_{-{\mathbf{k}},\uparrow},\tilde{b}_{-{\mathbf{k}},\downarrow}\right)\mathcal{M}({\mathbf{k}})\begin{pmatrix}\tilde{b}_{{\mathbf{k}}\uparrow}\\\
\tilde{b}_{{\mathbf{k}}\downarrow}\\\
\tilde{b}_{-{\mathbf{k}},\uparrow}^{\dagger}\\\
\tilde{b}_{-{\mathbf{k}},\downarrow}^{\dagger}\end{pmatrix},\end{split}$ (62)

where

$\begin{split}J({\mathbf{k}}):=\sum_{\alpha,\beta}N_{\beta}\left[\tilde{V}_{\alpha\beta}({\mathbf{k}},{\mathbf{0}},{\mathbf{0}},{\mathbf{k}})-\tilde{V}_{\alpha\beta}({\mathbf{0}},{\mathbf{0}},{\mathbf{0}},{\mathbf{0}})\right]+\sum_{\alpha}N_{\alpha}\tilde{V}_{\alpha\alpha}({\mathbf{k}},{\mathbf{0}},{\mathbf{k}},{\mathbf{0}}),\end{split}$
(63)

and the expression of the $4\times 4$ matrix $\mathcal{M}({\mathbf{k}})$ is
shown in Ref. [46].

To diagonalize Eq. (62), we perform the Bogoliubov transformation

$\begin{pmatrix}\tilde{b}_{{\mathbf{k}}\uparrow}\\\
\tilde{b}_{{\mathbf{k}}\downarrow}\\\
\tilde{b}_{-{\mathbf{k}},\uparrow}^{\dagger}\\\
\tilde{b}_{-{\mathbf{k}},\downarrow}^{\dagger}\end{pmatrix}=W({\mathbf{k}})\begin{pmatrix}\gamma_{{\mathbf{k}},1}\\\
\gamma_{{\mathbf{k}},2}\\\ \gamma_{-{\mathbf{k}},1}^{\dagger}\\\
\gamma_{-{\mathbf{k}},2}^{\dagger}\end{pmatrix},~{}W({\mathbf{k}})=\begin{pmatrix}{\cal
U}({\mathbf{k}})&{\cal V}^{*}(-{\mathbf{k}})\\\ {\cal V}({\mathbf{k}})&{\cal
U}^{*}(-{\mathbf{k}})\end{pmatrix}.$ (64)

Here, the paraunitary matrix $W({\mathbf{k}})$ is chosen to satisfy

$\tau_{3}\mathcal{M}({\mathbf{k}})W({\mathbf{k}})=W({\mathbf{k}})\mathrm{diag}(E_{1}({\mathbf{k}}),E_{2}({\mathbf{k}}),-E_{1}(-{\mathbf{k}}),-E_{2}(-{\mathbf{k}})),$
(65)

where $\tau_{3}:=\mathrm{diag}(1,1,-1,-1)$. Namely, $W({\mathbf{k}})$ and
$E_{j}({\mathbf{k}})~{}(j=1,2)$ are obtained by solving the right eigenvalue
problem of $\tau_{3}\mathcal{M}({\mathbf{k}})$. Equation (62) is then
diagonalized as

$\begin{split}H_{\mathrm{int}}=\frac{1}{2}\sum_{\alpha,\beta}N_{\alpha}N_{\beta}\tilde{V}_{\alpha\beta}({\mathbf{0}},{\mathbf{0}},{\mathbf{0}},{\mathbf{0}})-\frac{1}{2}\sum_{{\mathbf{k}}\neq{\mathbf{0}}}J({\mathbf{k}})+\sum_{{\mathbf{k}}\neq{\mathbf{0}}}\sum_{j=1,2}E_{j}({\mathbf{k}})\left(\gamma_{{\mathbf{k}}j}^{\dagger}\gamma_{{\mathbf{k}}j}+\frac{1}{2}\right).\end{split}$
(66)

We thus find that the Bogoliubov excitations (bogolons) are created by
$\gamma_{{\mathbf{k}}j}^{\dagger}~{}(j=1,2)$ and have the dispersion relations
$\\{E_{j}({\mathbf{k}})\\}$.

The ground state of Eq. (66) is given by the bogolon vacuum $|0\rangle$, which
is specified by the condition that $\gamma_{{\mathbf{k}}j}|0\rangle=0$ for all
${\mathbf{k}}\neq{\mathbf{0}}$ and $j=1,2$. The ground-state energy
$E_{\mathrm{GS}}$ (scaled by the interaction energy scale $gn^{2}A$) is
therefore given by

$\begin{split}\frac{E_{\mathrm{GS}}}{gn^{2}A}=&\frac{1}{2}\sum_{\alpha,\beta}\frac{A}{g}\tilde{V}_{\alpha\beta}({\mathbf{0}},{\mathbf{0}},{\mathbf{0}},{\mathbf{0}})+\frac{1}{gn\nu
N_{\mathrm{v}}}\sum_{{\mathbf{k}}\neq{\mathbf{0}}}\left[\sum_{j}E_{j}({\mathbf{k}})-J({\mathbf{k}})\right]\\\
=&\frac{1}{2}\sum_{\alpha,\beta}\frac{A}{g}\tilde{V}_{\alpha\beta}({\mathbf{0}},{\mathbf{0}},{\mathbf{0}},{\mathbf{0}})+\frac{1}{2\pi
gn\nu}\int_{\mathrm{BZ}}d^{2}{\mathbf{k}}\ell^{2}\left[\sum_{j}E_{j}({\mathbf{k}})-J({\mathbf{k}})\right].\end{split}$
(67)

Here, in the final expression, we take the thermodynamic limit
$N_{\mathrm{v}}\to\infty$ so that the sum is replaced by the integral over the
Brillouin zone as
$\frac{1}{N_{\mathrm{v}}}\sum_{{\mathbf{k}}\neq{\mathbf{0}}}\to\frac{1}{|{\mathbf{b}}_{1}\times{\mathbf{b}}_{2}|}\int\mathrm{d}^{2}{\mathbf{k}}=\frac{1}{2\pi}\int_{\mathrm{BZ}}\mathrm{d}^{2}{\mathbf{k}}\ell^{2}$;
this integral is convergent as the integrand is finite over the entire
Brillouin zone. The first term of Eq. (67) corresponds to the mean-field
ground-state energy, which has been analyzed by Mueller and Ho [36]. The other
term gives a quantum correction and is inversely proportional to the filling
factor $\nu$. In Sec. 4.5, we numerically calculate Eq. (67), and discuss how
the quantum correction affects the ground-state phase diagrams.

Using Eq. (64), we can calculate the following correlators in the ground
state:

$\displaystyle\langle
0|\tilde{b}_{{\mathbf{k}},\alpha}^{\dagger}\tilde{b}_{{\mathbf{k}},\alpha}|0\rangle=\sum_{j}|{\cal
V}_{\alpha,j}(-{\mathbf{k}})|^{2},~{}~{}\langle
0|\tilde{b}_{-{\mathbf{k}},\alpha}\tilde{b}_{-{\mathbf{k}},\alpha}^{\dagger}|0\rangle=\sum_{j}|{\cal
U}_{\alpha,j}(-{\mathbf{k}})|^{2},$ (68a) $\displaystyle\langle
0|\tilde{b}_{-{\mathbf{k}},\alpha}\tilde{b}_{{\mathbf{k}},\alpha}|0\rangle=\sum_{j}{\cal
U}_{\alpha,j}(-{\mathbf{k}}){\cal
V}_{\alpha,j}^{\ast}(-{\mathbf{k}}),~{}~{}\langle
0|\tilde{b}_{{\mathbf{k}},\alpha}^{\dagger}\tilde{b}_{-{\mathbf{k}},\alpha}^{\dagger}|0\rangle=\sum_{j}{\cal
U}_{\alpha,j}^{\ast}(-{\mathbf{k}}){\cal V}_{\alpha,j}(-{\mathbf{k}}).$ (68b)

Here, we have nonzero “anomalous” correlators in Eq. (68b) as the particle
numbers $N_{\uparrow}$ and $N_{\downarrow}$ are not conserved in the
Bogoliubov Hamiltonian (62). Using Eq. (68a), we can further calculate the
fraction of depletion $n^{\prime}/n$, which is equal for the two components,
as

$\frac{n^{\prime}}{n}=\frac{1}{N_{\alpha}}\sum_{{\mathbf{k}}\neq{\mathbf{0}}}\langle
0|\tilde{b}_{{\mathbf{k}},\alpha}^{\dagger}\tilde{b}_{{\mathbf{k}},\alpha}|0\rangle=\frac{2}{\nu
N_{\mathrm{v}}}\sum_{{\mathbf{k}}\neq{\mathbf{0}}}\sum_{j=1,2}|{\cal
V}_{\alpha,j}(-{\mathbf{k}})|^{2}~{}~{}(\alpha=\uparrow,\downarrow).$ (69)

As discussed in Sec. 2.6 [see Eq. (56)], this quantity is expected to diverge
logarithmically as a function of $N_{\mathrm{v}}$. We will confirm this
behavior numerically in Sec. 4.4. This diverging behavior comes from the
divergence of $|{\cal V}_{\alpha,2}(-{\mathbf{k}})|^{2}$ for
${\mathbf{k}}\to{\mathbf{0}}$ in Eq. (69). We note that the Bogoliubov theory
should be applied in the condition of weak depletion $n^{\prime}/n\ll 1$; this
condition is satisfied in typical experiments of ultracold atomic gases, where
$N_{\mathrm{v}}$ is at most of the order of 100 [12].

###  3.3 Intercomponent entanglement

For the RDM $\rho_{\uparrow}$ for the spin-$\uparrow$ component, we introduce
the following Gaussian ansatz [62, 63, 67, 68]:

$\rho_{\uparrow}=\frac{1}{Z_{\mathrm{e}}}{\mathrm{e}}^{-H_{\mathrm{e}}},~{}~{}H_{\mathrm{e}}=\frac{1}{2}\sum_{{\mathbf{k}}\neq{\mathbf{0}}}\left(\tilde{b}_{{\mathbf{k}},\uparrow}^{\dagger},\tilde{b}_{-{\mathbf{k}},\uparrow}\right)M_{\mathrm{e}}({\mathbf{k}})\begin{pmatrix}\tilde{b}_{{\mathbf{k}},\uparrow}\\\
\tilde{b}_{-{\mathbf{k}},\uparrow}^{\dagger}\end{pmatrix},~{}~{}Z_{\mathrm{e}}=\Tr~{}{\mathrm{e}}^{-H_{\mathrm{e}}},$
(70)

with

$M_{\mathrm{e}}({\mathbf{k}})=\begin{pmatrix}h_{\mathbf{k}}&-\lambda_{\mathbf{k}}\\\
-\lambda_{-{\mathbf{k}}}^{\ast}&h_{-{\mathbf{k}}}\end{pmatrix},~{}\lambda_{\mathbf{k}}=\lambda_{-{\mathbf{k}}}.$
(71)

By performing a Bogoliubov transformation

$\begin{pmatrix}\tilde{b}_{{\mathbf{k}}\uparrow}\\\
\tilde{b}_{-{\mathbf{k}},\uparrow}^{\dagger}\end{pmatrix}=W_{\mathrm{e}}({\mathbf{k}})\begin{pmatrix}\eta_{\mathbf{k}}\\\
\eta_{-{\mathbf{k}}}^{\dagger}\end{pmatrix},~{}~{}W_{\mathrm{e}}({\mathbf{k}})=\begin{pmatrix}\cosh\theta_{\mathbf{k}}&{\mathrm{e}}^{-\mathrm{i}\phi_{\mathbf{k}}}\sinh\theta_{\mathbf{k}}\\\
{\mathrm{e}}^{\mathrm{i}\phi_{\mathbf{k}}}\sinh\theta_{\mathbf{k}}&\cosh\theta_{\mathbf{k}}\end{pmatrix}$
(72)

with

$\cosh
2\theta_{\mathbf{k}}=\frac{h_{\mathbf{k}}+h_{-{\mathbf{k}}}}{2\tilde{h}_{\mathbf{k}}},~{}~{}{\mathrm{e}}^{-i\phi_{\mathbf{k}}}\sinh
2\theta_{\mathbf{k}}=\frac{\lambda_{\mathbf{k}}}{\tilde{h}_{\mathbf{k}}},~{}~{}\tilde{h}_{\mathbf{k}}=\sqrt{\frac{1}{4}\left(h_{\mathbf{k}}+h_{-{\mathbf{k}}}\right)^{2}-|\lambda_{\mathbf{k}}|^{2}},$
(73)

the entanglement Hamiltonian in $H_{\mathrm{e}}$ in Eq. (70) is diagonalized
as

$H_{\mathrm{e}}=\frac{1}{2}\sum_{{\mathbf{k}}\neq{\mathbf{0}}}\left(\xi_{\mathbf{k}}\eta_{\mathbf{k}}^{\dagger}\eta_{\mathbf{k}}+\xi_{-{\mathbf{k}}}\eta_{-{\mathbf{k}}}\eta_{-{\mathbf{k}}}^{\dagger}\right)=\sum_{{\mathbf{k}}\neq{\mathbf{0}}}\xi_{\mathbf{k}}\left(\eta_{\mathbf{k}}^{\dagger}\eta_{\mathbf{k}}+\frac{1}{2}\right),$
(74)

where

$\xi_{\mathbf{k}}:=\tilde{h}_{\mathbf{k}}+\frac{h_{\mathbf{k}}-h_{-{\mathbf{k}}}}{2}$
(75)

is the single-particle ES.

Using the relation (72) and the Bose distribution functions

$\Tr\left(\eta_{\mathbf{k}}^{\dagger}\eta_{\mathbf{k}}\rho_{\uparrow}\right)=\frac{1}{{\mathrm{e}}^{\xi_{\mathbf{k}}}-1}=f_{\mathrm{B}}(\xi_{\mathbf{k}}),~{}~{}\Tr\left(\eta_{-{\mathbf{k}}}\eta_{-{\mathbf{k}}}^{\dagger}\rho_{\uparrow}\right)=1+f_{\mathrm{B}}(\xi_{-{\mathbf{k}}})=-f_{\mathrm{B}}(-\xi_{-{\mathbf{k}}}),$
(76)

we obtain

$\displaystyle\Tr\left(\tilde{b}_{{\mathbf{k}},\uparrow}^{\dagger}\tilde{b}_{{\mathbf{k}},\uparrow}\rho_{\uparrow}\right)$
$\displaystyle=f_{\mathrm{B}}(\xi_{\mathbf{k}})\cosh^{2}\theta_{\mathbf{k}}-f_{\mathrm{B}}(-\xi_{-{\mathbf{k}}})\sinh^{2}\theta_{\mathbf{k}},$
(77a)
$\displaystyle\Tr\left(\tilde{b}_{-{\mathbf{k}},\uparrow}\tilde{b}_{-{\mathbf{k}},\uparrow}^{\dagger}\rho_{\uparrow}\right)$
$\displaystyle=-f_{\mathrm{B}}(-\xi_{-{\mathbf{k}}})\cosh^{2}\theta_{\mathbf{k}}+f_{\mathrm{B}}(\xi_{\mathbf{k}})\sinh^{2}\theta_{\mathbf{k}},$
(77b) $\displaystyle
2\Tr\left(\tilde{b}_{-{\mathbf{k}},\uparrow}\tilde{b}_{{\mathbf{k}},\uparrow}\rho_{\uparrow}\right)$
$\displaystyle=\left[f_{\mathrm{B}}(\xi_{\mathbf{k}})-f_{\mathrm{B}}(-\xi_{-{\mathbf{k}}})\right]{\mathrm{e}}^{-i\phi_{\mathbf{k}}}\sinh\left(2\theta_{\mathbf{k}}\right).$
(77c)

We require these to be equal to the correlators (68) with respect to the
Bogoliubov ground state. We can then express
$f_{\mathrm{B}}(\xi_{\mathbf{k}})$ in terms of the correlators (68) in the
following way. First, by taking the sum and the difference of Eqs. (77a) and
(77b), we have

$\displaystyle\langle\tilde{b}_{{\mathbf{k}}\uparrow}^{\dagger}\tilde{b}_{{\mathbf{k}}\uparrow}\rangle+\langle\tilde{b}_{-{\mathbf{k}},\uparrow}\tilde{b}_{-{\mathbf{k}},\uparrow}^{\dagger}\rangle$
$\displaystyle=\left[f_{\mathrm{B}}(\xi_{\mathbf{k}})-f_{\mathrm{B}}(-\xi_{-{\mathbf{k}}})\right]\cosh\left(2\theta_{\mathbf{k}}\right),$
(78a)
$\displaystyle\langle\tilde{b}_{{\mathbf{k}}\uparrow}^{\dagger}\tilde{b}_{{\mathbf{k}}\uparrow}\rangle-\langle\tilde{b}_{-{\mathbf{k}},\uparrow}\tilde{b}_{-{\mathbf{k}},\uparrow}^{\dagger}\rangle$
$\displaystyle=f_{\mathrm{B}}(\xi_{\mathbf{k}})+f_{\mathrm{B}}(-\xi_{-{\mathbf{k}}}),$
(78b)

where we take the shorthand notation $\langle\cdot\rangle:=\langle
0|\cdot|0\rangle$. Next, using Eqs. (77c) and (78a), we have

$f_{\mathrm{B}}(\xi_{\mathbf{k}})-f_{\mathrm{B}}(-\xi_{-{\mathbf{k}}})=\sqrt{\left(\langle\tilde{b}_{{\mathbf{k}}\uparrow}^{\dagger}\tilde{b}_{{\mathbf{k}}\uparrow}\rangle+\langle\tilde{b}_{-{\mathbf{k}},\uparrow}\tilde{b}_{-{\mathbf{k}},\uparrow}^{\dagger}\rangle\right)^{2}-4\big{|}\langle\tilde{b}_{-{\mathbf{k}},\uparrow}\tilde{b}_{{\mathbf{k}}\uparrow}\rangle\big{|}^{2}}.$
(79)

Lastly, using Eqs. (78b) and (79), we obtain

$f_{\mathrm{B}}(\xi_{\mathbf{k}})=\frac{1}{2}\left(\langle\tilde{b}_{{\mathbf{k}}\uparrow}^{\dagger}\tilde{b}_{{\mathbf{k}}\uparrow}\rangle-\langle\tilde{b}_{-{\mathbf{k}},\uparrow}\tilde{b}_{-{\mathbf{k}},\uparrow}^{\dagger}\rangle\right)+\sqrt{\frac{1}{4}\left(\langle\tilde{b}_{{\mathbf{k}}\uparrow}^{\dagger}\tilde{b}_{{\mathbf{k}}\uparrow}\rangle+\langle\tilde{b}_{-{\mathbf{k}},\uparrow}\tilde{b}_{-{\mathbf{k}},\uparrow}^{\dagger}\rangle\right)^{2}-\big{|}\langle\tilde{b}_{-{\mathbf{k}},\uparrow}\tilde{b}_{{\mathbf{k}}\uparrow}\rangle\big{|}^{2}},$
(80)

from which we can calculate the single-particle ES
$\xi_{\mathbf{k}}=\ln\left[1+f_{B}(\xi_{\mathbf{k}})^{-1}\right]$.

We can further use Eq. (80) to calculate the intercomponent EE as

$S_{\mathrm{e}}=\sum_{{\mathbf{k}}\neq{\mathbf{0}}}\left\\{-f_{\mathrm{B}}(\xi_{\mathbf{k}})\ln
f_{\mathrm{B}}(\xi_{\mathbf{k}})+\left[1+f_{\mathrm{B}}(\xi_{\mathbf{k}})\right]\ln\left[1+f_{\mathrm{B}}(\xi_{\mathbf{k}})\right]\right\\}.$
(81)

As discussed in Sec. 2.5 [see Eq. (50)], $S_{\mathrm{e}}$ is expected to show
a volume-law behavior with a subleading logarithimic correction. The EE per
flux quantum in the thermodynamic limit is then expressed in the integral form

$\lim_{N_{\mathrm{v}}\to\infty}\frac{S_{\mathrm{e}}}{N_{\mathrm{v}}}=\frac{1}{2\pi}\int_{\mathrm{BZ}}\mathrm{d}^{2}{\mathbf{k}}\ell^{2}\left\\{-f_{\mathrm{B}}(\xi_{\mathbf{k}})\ln
f_{\mathrm{B}}(\xi_{\mathbf{k}})+\left[1+f_{\mathrm{B}}(\xi_{\mathbf{k}})\right]\ln\left[1+f_{\mathrm{B}}(\xi_{\mathbf{k}})\right]\right\\}.$
(82)

We note that the correlators (68) are independent of $\nu$ once the lattice
structure is fixed. Therefore, the EE per flux quantum in Eq. (82) is also
independent of $\nu$ in a similar manner. In Sec. 4.3, we will calculate the
EE per flux quantum assuming the structure in the mean-field ground state.

##  4 Numerical results

In this section, we present numerical results that are obtained using the
Bogoliubov theory formulation of Sec. 3. In this formulation, one starts from
the lattice structure with the primitive vectors ${\mathbf{a}}_{1}$ and
${\mathbf{a}}_{2}$ and the displacement parameters $u_{1}$ and $u_{2}$ [see
Fig. 1(f)]. In the GP mean-field theory, these parameters are determined so as
to minimize the mean-field ground-state energy [the first term on the right-
hand side of Eq. (67)] for a fixed magnetic length $\ell$. As demonstrated by
Mueller and Ho [36], this mean-field analysis gives a rich phase diagram that
consists of five different vortex-lattice phases as shown in Fig. 1. In Secs.
4.1, 4.2, 4.3, and 4.4, we analyze renormalized coupling constants, excitation
spectra, intercomponent entanglement, and the fraction of depletion,
respectively, using the Bogoliubov theory based on the mean-field vortex
lattice structures. Therefore, the results in these sections correspond to the
case of $\nu=\infty$. As we lower the filling factor $\nu$, quantum
fluctuations are expected to affect the vortex lattice structures and the
ground-state phase diagrams. In Sec. 4.5, we investigate how quantum
fluctuations affect the ground-state phase diagrams for parallel and
antiparallel fields by calculating a quantum correction to the ground-state
energy [i.e., the term proportional to $\nu^{-1}$ in Eq. (67)].

###  4.1 Renormalized coupling constants

Figure 2:  Renormalization factors $\beta:=\bar{g}/g$ (black) and
$\beta_{\uparrow\downarrow}:=\bar{g}_{\uparrow\downarrow}/g_{\uparrow\downarrow}$
(green) for the intracomponent and intercomponent coupling constants. These
factors are calculated using Eq. (83) for the mean-field vortex lattice
structure for each $g_{\uparrow\downarrow}/g$, and take the same values for
parallel and antiparallel fields. We set $N_{\mathrm{v}}=97^{2}$, with which a
sufficient convergence to the thermodynamic limit is achieved. For overlapping
triangular lattices ($-1<g_{\uparrow\downarrow}/g<0$), we have
$\beta=\beta_{\uparrow\downarrow}$ and thus the two curves overlap; the value
$\beta=\beta_{\uparrow\downarrow}=1.1596$ in this region coincides with the
renormalization factor obtained for a scalar BEC [7, 18, 66]. Vertical dashed
lines indicate the transition points in the mean-field phase diagram [36].

We first determine the renormalized coupling constants $\bar{g}_{\alpha\beta}$
or equivalently, the renormalization factors
$\beta_{\alpha\beta}:=\bar{g}_{\alpha\beta}/g_{\alpha\beta}$ that are
introduced in Sec. 2.2. Comparing the mean-field ground-state energy [the
first term of Eq. (67)] with the corresponding field-theoretical expression
(11), we find

$\bar{g}_{\alpha\beta}=\beta_{\alpha\beta}g_{\alpha\beta}=AV_{\alpha\beta}({\mathbf{q}}_{\alpha},{\mathbf{q}}_{\beta},{\mathbf{q}}_{\beta},{\mathbf{q}}_{\alpha})=Ag_{\alpha\beta}\int\mathrm{d}^{2}\mathbf{r}~{}|\Psi_{{\mathbf{q}}_{\alpha},\alpha}(\mathbf{r})|^{2}|\Psi_{{\mathbf{q}}_{\beta},\beta}(\mathbf{r})|^{2}.$
(83)

This expression can also be obtained by substituting the condensate wave
function $\sqrt{N_{\alpha}}\Psi_{{\mathbf{q}}_{\alpha},\alpha}(\mathbf{r})$
into $\psi_{\alpha}({\mathbf{r}})$ in Eq. (6). As seen in this expression,
$\bar{g}_{\alpha\beta}$ is determined from the contribution of each
interaction term to the mean-field ground-state energy.

Figure 2 shows the renormalization factors
$\beta:=\beta_{\uparrow\uparrow}=\beta_{\downarrow\downarrow}$ and
$\beta_{\uparrow\downarrow}=\beta_{\downarrow\uparrow}$ calculated for the
mean-field vortex lattice structures. For overlapping triangular, interlaced
triangular, and square lattices, $\beta$ and $\beta_{\uparrow\downarrow}$ do
not depend on $g_{\uparrow\downarrow}/g$ as the lattice structures remain
unchanged in the concerned regions. For rhombic and rectangular lattices, in
contrast, $\beta$ and $\beta_{\uparrow\downarrow}$ do depend on
$g_{\uparrow\downarrow}/g$ as the inner angle $\theta$ and the aspect ratio
$b/a$ continuously vary for the former and the latter, respectively.

We have argued in Sec. 2.2 that the intracomponent coupling is always enhanced
by the renormalization. We indeed find $\beta>1$ in all the regions in Fig. 2.
We have also argued that the intercomponent repulsion
$g_{\uparrow\downarrow}>0$ (attraction $g_{\uparrow\downarrow}<0$) is reduced
(enhanced) by the renormalization owing to the displacement (overlap) of
vortices between the components. In Fig. 2, we indeed find
$\beta_{\uparrow\downarrow}<1$ ($\beta_{\uparrow\downarrow}>1$) for
$g_{\uparrow\downarrow}>0$ ($g_{\uparrow\downarrow}<0$). Furthermore, $\beta$
($\beta_{\uparrow\downarrow}$) monotonically increases (decreases) as a
function of $g_{\uparrow\downarrow}/g$ for $g_{\uparrow\downarrow}>0$. This
reflects the fact that with increasing $g_{\uparrow\downarrow}/g$, vortices in
different components tend to repel more strongly with each other at the cost
of increasing the intracomponent interaction energy.

###  4.2 Excitation spectrum and elastic constants

Figure 3:  Bogoliubov excitation spectra $\\{E_{j}({\mathbf{k}})\\}$ (in units
of $gn$) for different values of $g_{\uparrow\downarrow}/g$ corresponding to
the lattice structures (a)-(e) shown in Fig. 1. In each panel, both results
for parallel (blue) and antiparallel (red) magnetic fields are shown.
Calculations are done along dotted arrows in the lower panels of Fig. 1.
Figure 4:  Rescaled Bogoliubov excitation spectra for the five cases (a)-(e)
shown in Fig. 3. Blue curves show
$E_{1}^{\mathrm{P}}({\mathbf{k}})/(\sqrt{g\bar{g}_{-}}n)$ and
$E_{2}^{\mathrm{P}}({\mathbf{k}})/(\sqrt{g\bar{g}_{+}}n)$ for parallel (P)
fields while red curves show
$E_{1}^{\mathrm{AP}}({\mathbf{k}})/(\sqrt{g\bar{g}_{+}}n)$ and
$E_{2}^{\mathrm{AP}}({\mathbf{k}})/(\sqrt{g\bar{g}_{-}}n)$ for antiparallel
(AP) fields. Here, $E_{j}^{\mathrm{P/AP}}({\mathbf{k}})$ with $j=1$ and $2$
correspond to the upper and lower excitation bands, respectively. We can
confirm that the blue and red curves overlap at sufficiently low energies
around the $\Gamma$ point, which indicates the rescaling relations in Eq.
(36).  Figure 5:  Dimensionless functions
$f_{2}^{\mathrm{P}/\mathrm{AP}}(\varphi)$ (left) and
$f_{1}^{\mathrm{P}/\mathrm{AP}}(\varphi)$ (right) for parallel ($\mathrm{P}$;
blue) and antiparallel ($\mathrm{AP}$; red) fields for the same cases (a)-(e)
as in Figs. 3 and 4. These functions express the anisotropy of the low-energy
spectra as in Eq. (84), and they are calculated from the
$\\{E_{j}({\mathbf{k}})\\}$ along a circular path
${\mathbf{k}}=k(\cos\varphi,\sin\varphi)$ with $k=0.001a/\ell^{2}$ and
$\varphi\in[0,2\pi)$. With proper rescaling as in Eq. (37), the curves for
parallel and antiparallel fields are found to agree perfectly up to numerical
precision. Specifically, in the left panels,
$f_{2}^{\mathrm{P}}(\varphi)\sqrt{g/\bar{g}_{+}}$ and
$f_{2}^{\mathrm{AP}}(\varphi)\sqrt{g/\bar{g}_{-}}$ give the common function
$\sqrt{C(\varphi)}$ (black). In the right panels,
$f_{1}^{\mathrm{P}}(\varphi)\sqrt{g/\bar{g}_{-}}$ and
$f_{1}^{\mathrm{AP}}(\varphi)\sqrt{g/\bar{g}_{+}}$ give the common function
$\sqrt{D(\varphi)}$ (black).

As explained in Sec. 3.2 [see Eq. (65)], the excitation spectrum
$E_{j}({\mathbf{k}})~{}(j=1,2)$ can be obtained by numerically calculating the
right eigenvalues of the $4\times 4$ matrix
$\tau_{3}\mathcal{M}({\mathbf{k}})$. Figure 3 presents spectra obtained in
this way for all the lattice structures (a)-(e) shown in Fig. 1 and for both
parallel and antiparallel fields. In Ref. [46],666 There are errors in the
scales of some figures in Ref. [46]. Specifically, the numerical data for the
vertical axes in Figs. 2, 4, and 5 should be multiplied by $1/4$,
$1/\left(2\sqrt{2}\right)$, and $1/8$, respectively. These errors are
unrelated to the issue of renormalization discussed in the present paper. As
our understanding of the rescaling relations is now updated from Ref. [46],
Figs. 3, 5, and 6 of the present paper could be seen as improved versions of
these figures.  we discuss various unique features of these spectra such as
linear and quadratic dispersion relations at low energies and the emergence of
line and point nodes at high energies that are related to a fractional
translation symmetry. Here, we aim to demonstrate the rescaling relations in
Eqs. (36) and (37) which are predicted by the low-energy effective field
theory. In Ref. [46], the unrenormalized coupling constants $g_{\alpha\beta}$
are used for the rescaling, which leads to an incorrect conclusion that the
rescaling relations hold only for overlapping triangular lattices. Using the
renormalized coupling constants obtained in Sec. 4.1, we can demonstrate the
rescaling relations for all the five structures (a)-(e) shown in Fig. 1.

Figure 4 displays rescaled excitation spectra for the five cases (a)-(e) in
Fig. 3. Here, $\bar{g}_{\pm}:=\bar{g}\pm\bar{g}_{\uparrow\downarrow}=\beta
g\pm\beta_{\uparrow\downarrow}g_{\uparrow\downarrow}$ are used for the
rescaling, where the renormalization factors $\beta$ and
$\beta_{\uparrow\downarrow}$ are shown in Fig. 2. We can confirm that the
rescaling relations in Eq. (36) hold at sufficiently low energies around the
$\Gamma$ point. Interestingly, in Figs. 4(a) and (b), the rescaling relations
hold approximately up to high energies, which is beyond the scope of effective
field theory. At low energies, the spectra in Fig. 3 can be fit well by linear
and quadratic dispersion relations

$E_{j}({\mathbf{k}})=\sqrt{2}gn(k\ell)^{j}f_{j}(\varphi)~{}~{}(j=1,2),$ (84)

where the wave vector is parametrized as
${\mathbf{k}}=k(\cos\varphi,\sin\varphi)$ and $\\{f_{j}(\varphi)\\}$ are
dimensionless functions that characterize the anisotropy of the spectrum.
Figure 5 shows the functions $\\{f_{j}(\varphi)\\}$ obtained numerically for
the same cases as in Figs. 3 and 4. We find that with proper rescaling as in
Eq. (37), the curves for parallel and antiparallel fieds coincide perfectly up
to numerical precision, giving the funcions $\sqrt{C(\varphi)}$ and
$\sqrt{D(\varphi)}$ that are related to the elastic constants.

Figure 6:  Dimensionless elastic constants $C_{i}~{}(i=1,2,3,4)$ (left) and
$D_{i}~{}(i=1,2,3)$ (right) for the lattice structures (a)-(e) shown in Fig.
1. These constants are obtained by fitting the numerically obtained functions
$\sqrt{C(\varphi)}$ and $\sqrt{D(\varphi)}$ (common for parallel and
antiparallel fields as shown in Fig. 5) using Eq. (33). See Eq. (15) for the
symmetry constraints on the constants. Vertical dashed lines indicate the
transition points in the mean-field ground state.

By comparing the obtained $\sqrt{C(\varphi)}$ and $\sqrt{D(\varphi)}$ with the
analytical expressions in Eq. (33), we can determine the dimensionless elastic
constants $\\{C_{i}\\}$ and $\\{D_{i}\\}$. Figure 6 shows the determined
elastic constants as functions of $g_{\uparrow\downarrow}/g$, which are common
for parallel and antiparallel fields. In our previous work [46], we obtained
different elastic constants for the two types of fields as we missed the
necessity of using the renormalized coupling constants in relating the spectra
to the elastic constants. Figure 6 is essentially consistent with the elastic
constants determined by a different method by Keçeli and Oktel [45] except
that the constant $C_{4}$ in interlaced triangular lattices (b) was missed in
Ref. [45].

###  4.3 Intercomponent entanglement spectrum and entropy

Figure 7:  Left panels: single-particle ES
$\xi_{\mathbf{k}}^{\mathrm{P}/\mathrm{AP}}$ for parallel (P; blue) and
antiparallel (AP; red) fields for the same cases (a)-(e) as in Figs. 3, 4, and
5. Calculations are based on Eq. (80) and done along the paths indicated by
dotted arrows shown in the lower panels of Fig. 1. The spectra
$\xi_{\mathbf{k}}$ show divergence at some high-symmetry points and along
lines in the first Brillouin zone. In (d) and (e), in particular, divergence
occurs along the edges of the Brillouin zone for antiparallel fields and thus
the value of $\xi_{\mathbf{k}}$ is not shown there. Right panels:
dimensionless functions $c^{\mathrm{P}/\mathrm{AP}}(\varphi)$ that express the
anisotropy of the ES around the $\Gamma$ point as in Eq. (85). These are
calculated from $\xi_{\mathbf{k}}^{\mathrm{P}/\mathrm{AP}}$ along a circular
path ${\mathbf{k}}=k(\cos\varphi,\sin\varphi)$ with $k=0.001a/\ell^{2}$ and
$\varphi\in[0,2\pi)$. With proper rescaling, the curves for parallel and
antiparallel fields coincide perfectly up to numerical precision, confirming
the rescaling relation in Eq. (46). Namely,
$(\bar{g}_{+}/\bar{g}_{-})^{1/4}c^{\mathrm{P}}(\varphi)$ and
$(\bar{g}_{-}/\bar{g}_{+})^{1/4}c^{\mathrm{AP}}(\varphi)$ share the same
curves shown in black.

Here we present numerical results on the intercomponent ES and EE in the
Bogoliubov ground state. Calculations are based on the formulation in Sec.
3.3. The obtained results are compared with the field-theoretical results in
Sec. 2.5.

The left panels of Fig. 7 display the single-particle ES $\xi_{\mathbf{k}}$
for the same cases as in Figs. 3, 4, and 5. Around the $\Gamma$ point,
$\xi_{\mathbf{k}}$ can be fit well by the square-root dispersion relation

$\xi_{\mathbf{k}}=c(\varphi)\sqrt{k\ell}$ (85)

for ${\mathbf{k}}=k(\cos\varphi,\sin\varphi)$, where $c(\varphi)$ is a
dimensionless function that expresses the anisotropy. The right panels of Fig.
7 show the function $c(\varphi)$ determined from the data of
$\xi_{\mathbf{k}}$ along a circular path around the $\Gamma$ point. With
proper rescaling, the curves for parallel and antiparallel fields are found to
coincide perfectly up to numerical precision; furthermore, the rescaled curves
are found to agree accurately with $4[C(\varphi)/D(\varphi)]^{1/2}$ (not
shown), where $\sqrt{C(\varphi)}$ and $\sqrt{D(\varphi)}$ are shown in Fig. 5.
Thus, we can confirm the rescaling relation for the entanglement spectra in
Eq. (46).

In the left panels of Fig. 7, we also find that $\xi_{\mathbf{k}}$ diverges at
some high-symmetry points and along lines in the Brillouin zone. For square
(d) and rectangular (e) lattices in antiparallel fields, in particular,
divergence occurs along the edges of the Brillouin zone (thus
$\xi_{\mathbf{k}}$ is not shown along the paths $R\to M_{1}$ and $R\to
M_{2}$). In our previous work [46] (see Appendix D therein), it has been found
that at the $M_{1}$ and $M_{2}$ points for rhombic, square and rectangular
lattices and for both parallel and antiparallel fields, the Bogoliubov
Hamiltonian matrix $\mathcal{M}({\mathbf{k}})$ has the structure in which the
spin-$\uparrow$ and $\downarrow$ components are decoupled. Therefore,
$\xi_{\mathbf{k}}$ naturally diverges at these points. However, we have not
been able to find such a simple structure of $\mathcal{M}({\mathbf{k}})$ along
the edges of the Brillouin zone for square (d) and rectangular (e) lattices in
antiparallel fields.

Figure 8:  (a) Intercomponent EE per flux quantum,
$S_{\mathrm{e}}/N_{\mathrm{v}}$, as a function of $g_{\uparrow\downarrow}/g$
for parallel (blue) and antiparallel (red) fields for $N_{\mathrm{v}}=97^{2}$.
This is calculated from numerically obtained $\xi_{\mathbf{k}}$ using Eq.
(81). Vertical dashed lines indicate the transition points. (b) Intercomponent
EE $S_{\mathrm{e}}$ versus $N_{\mathrm{v}}$ for
$g_{\uparrow\downarrow}/g=0.75$. A fit to the form
$S_{\mathrm{e}}=\alpha_{1}N_{\mathrm{v}}-\alpha_{2}\ln
N_{\mathrm{v}}+\alpha_{3}$ gives
$(\alpha_{1},\alpha_{2},\alpha_{3})=(0.560,0.236,-0.70)$ for parallel fields
and $(\alpha_{1},\alpha_{2},\alpha_{3})=(0.130,0.243,-0.06)$ for antiparallel
fields.

Figure 8(a) shows the intercomponent EE per flux quantum as a function of the
ratio $g_{\uparrow\downarrow}/g$ for parallel (blue) and antiparallel (red)
fields. We find that for repulsive (attractive) $g_{\uparrow\downarrow}$, the
EE tends to be larger for parallel (antiparallel) fields, in consistency with
the field-theoretical result in Sec. 2.5 [see Eq. (51)]. This behavior is in
qualitative agreement with the exact diagonalization results in a quantum
(spin) Hall regime with $\nu=\mathcal{O}(1)$ [42, 52], in which product states
of nearly uncorrelated quantum Hall states are found to be robust for an
intercomponent attraction (repulsion) in the case of parallel (antiparallel)
fields.

In Figs. 8(b), we examine the scaling of the intercomponent EE
$S_{\mathrm{e}}$ as a function of $N_{\mathrm{v}}$. As seen in this figure,
the dominant part of the scaling is given by a volume-law behavior
$\alpha_{1}N_{\mathrm{v}}$; such a volume-law contribution is standard for an
extensive cut as discussed here, and has also been found elsewhere [63, 62,
76, 77, 78]. In agreement with the field-theoretical result in Eq. (50), we
find that the data are well fitted by the form
$S_{\mathrm{e}}=\alpha_{1}N_{\mathrm{v}}-\alpha_{2}\ln
N_{\mathrm{v}}+\alpha_{3}$ and that the coefficient $\alpha_{2}$ obtained by
the fitting is close to $1/4$.

###  4.4 Fraction of depletion

Figure 9:  (a) Fraction of depletion $n^{\prime}/n$ (scaled by $\nu^{-1}$) as
a function of $g_{\uparrow\downarrow}/g$ for parallel (blue) and antiparallel
(red) fields for $N_{\mathrm{v}}=5^{2}$ (dotted), $11^{2}$ (dashed), and
$19^{2}$ (solid). The calculations are based on Eq. (69). Vertical dashed
lines indicate the transition points. (b) $\nu n^{\prime}/n$ versus
$N_{\mathrm{v}}$ for $g_{\uparrow\downarrow}/g=0.75$. A fit to the form $\nu
n^{\prime}/n=\gamma_{1}\ln N_{\mathrm{v}}+\gamma_{2}$ gives
$(\gamma_{1},\gamma_{2})=(1.91,-1.53)$ for parallel fields and
$(\gamma_{1},\gamma_{2})=(1.09,-0.40)$ for antiparallel fields.

Figure 9 shows numerical results on the fraction of depletion $n^{\prime}/n$
(scaled by $\nu^{-1}$). As seen in Fig. 9(a), for an intercomponent repulsion
(attraction), this quantity tends to be larger for parallel (antiparallel)
fields, indicating stronger quantum fluctuations. This is in agreement with
the field-theoretical result in Eq. (56). At the transition point between
interlaced triangular and rhombic lattices, the fraction of depletion changes
discontinuously owing to a discontinuous change in the lattice structure.
Meanwhile, the fraction of depletion diverges at both the transition points
between rhombic, square, and rectangular lattices; this seems to be related to
rapid changes in the inner angle and the aspect ratio in the mean-field ground
state as shown in Fig. 11. In Fig. 9(b), we examine the scaling of $\nu
n^{\prime}/n$ as a function of $N_{\mathrm{v}}$. The data are well fitted by
the logarithmic form $\nu n^{\prime}/n=\gamma_{1}\ln
N_{\mathrm{v}}+\gamma_{2}$, in agreement with the field-theoretical result in
Eq. (56).

###  4.5 Ground-state phase diagrams

Figure 10:  (a) Ground-state energy [Eq. (67)] as a function of
$g_{\uparrow\downarrow}/g$. The mean-field energy (black and grey) is
appreciably changed by a quantum correction, as shown for both parallel (blue
and purple) and antiparallel (red and brown) fields for $\nu=20$. The mean-
field energy for the square lattices,
$E^{\mathrm{squ}}_{\mathrm{MF}}/(gn^{2}A)=1.18034+0.834627g_{\uparrow\downarrow}/g$,
is subtracted to emphasize the changes due to quantum fluctuations. Vertical
dashed lines indicate the transition points, and alternating colors correspond
to different phases. In particular, the transition points between rhombic-,
square- and rectangular lattices shift appreciably owing to quantum
corrections. In contrast, the transition point $g_{\uparrow\downarrow}/g=0$
between overlapping and interlaced triangular lattices remain unchanged by
quantum corrections. (b) Enlarged view of (a) around the transition point
between interlaced triangular and rhombic lattices, which also shows small
shifts due to quantum corrections.

Here we analyze how quantum fluctuations affect the ground-state phase
diagrams for parallel and antiparallel fields. The GP mean-field analyses [36,
37, 38, 39] have led to the five types of lattice structures that depend on
$g_{\uparrow\downarrow}/g$ as shown in Fig. 1. We assume that the same types
of structures appear in the presence of quantum fluctuations as well, and
examine a quantum correction to the ground-state energy.

Figure 11:  (a) The inner angle $\theta$ of rhombic lattices and (b) the
aspect ratio $b/a$ of rectangular lattices, plotted against
$g_{\uparrow\downarrow}/g$. The mean-field results (black) [36] are changed
appreciably by quantum corrections, as shown both for parallel (blue) and
antiparallel (red) fields. Insets show the definitions of $\theta$ and $b/a$.

Figure 10 shows the mean-field ground-state energy as well as those with
quantum corrections for parallel and antiparallel fields, where the filling
factor is $\nu=20$. Here, the energies for rhombic and rectangular lattices
are minimized with respect to the inner angle $\theta$ and the aspect ratio
$b/a$, respectively. As seen in Fig. 10 (a), the transition points between
rhombic, square and rectangular lattices shift appreciably owing to quantum
corrections; the shift occurs more significantly for parallel fields. In
contrast, the transition point $g_{\uparrow\downarrow}/g=0$ between
overlapping and interlaced triangular lattices remains unchanged by quantum
corrections. While the shift of the transition point between interlaced
triangular and rhombic lattices is not clearly seen in Fig 10(a), the shift
indeed occurs as seen in the enlarged plot in Fig. 10(b).

Figure 11(a) and (b) show the inner angle $\theta$ of rhombic lattices and the
aspect ratio $b/a$ of rectangular lattices, respectively, which are obtained
through the one-parameter minimization of the ground-state energy (67). In the
mean-field results, both the inner angle $\theta$ and the aspect ratio $b/a$
show rapid changes near the transition points to square lattices [46]. In
these regimes, the systems are expected to be highly susceptible to quantum
fluctuations. Indeed, changes in $\theta$ and $b/a$ are enhanced in these
regimes, which explains the large shifts of the transition points as
demonstrated in Fig. 10(a).

##  5 Summary and outlook

In this paper, we have presented a detailed comparative study of vortex
lattices of binary BECs in parallel and antiparallel synthetic magnetic
fields. Within the GP mean-field theory valid for high filling factors $\nu\gg
1$, the two types of fields are known to lead to the same phase diagram that
consists of a variety of vortex lattices [36, 37, 38, 39, 42]. We have
formulated an improved effective field theory for such vortex lattices by
introducing renormalized coupling constants for coarse-grained densities, and
studied properties of collective modes and ground-state intercomponent
entanglement. We have also performed numerical calculations based on the
Bogoliubov theory with the LLL approximation to confirm the field-theoretical
predictions. We have shown that the low-energy excitation spectra for the two
types of fields are related to each other by suitable rescaling using the
renormalized coupling constants [Eqs. (36) and (37) and Figs. 4 and 5]. By
calculating the intercomponent EE in the ground state, we have found that for
an intercomponent repulsion (attraction), the two components are more strongly
entangled in the case of parallel (antiparallel) fields [Eq. (51) and Fig.
8(a)], in qualitative agreement with recent numerical results for a quantum
(spin) Hall regime [42, 52]. As a by-product, we have also found that the ES
exhibits an anomalous square-root dispersion relation [Eq. (43) and Fig. 7],
and that the EE exhibits a volume-law scaling followed by a subleading
logarithmic term [Eq. (50) and Fig. 8(b)]. Finally, we have investigated the
effects of quantum fluctuations on the phase diagrams by calculating the
correction to the ground-state energy due to zero-point fluctuations in the
Bogoliubov theory [Eq. (67) and Fig. 10]. We have found that the boundaries
between rhombic-, square-, and rectangular-lattice phases shift appreciably
with a decrease in $\nu$.

We have seen a similarity between the regimes of high ($\nu\gg 1$) and low
($\nu=O(1)$) filling factors in the behavior of intercomponent entanglement.
It will be interesting to investigate how the two regimes are connected by
applying sophisticated numerical methods such as a variational wave function
[22] and the infinite density matrix renormalization group [51]. Furthermore,
the similarity between the two regimes suggests that the behavior of
intercomponent entanglement does not depend on the details of the systems and
can be universal for a wide range of Hamiltonians. In fact, it has been found
in lattice models that two coupled bosonic Laughlin states with opposite
chiralities (i.e., fractional quantum spin Hall states [44]) are more robust
against an intercomponent repulsion than the ones with the same chirality
[79]. The stability of fractional quantum spin Hall states against an
intercomponent repulsion has also been discussed in fermionic models [80, 81,
82, 83]. Comparative investigation of multicomponent systems in gauge fields
as conducted in the present work will be a useful approach for exploring
universal features of interacting topological states of matter.

## Acknowledgments

TY is grateful to Kazuya Fujimoto, Yuki Kawaguchi, and Terumichi Ohashi for
stimulating discussions during his stay at Nagoya University. This work was
supported by KAKENHI Grant No. JP18H01145 and No. JP18K03446, a Grant-in-Aid
for Scientific Research on Innovative Areas “Topological Materials Science”
(KAKENHI Grant No. JP15H05855) from the Japan Society for the Promotion of
Science (JSPS), and Keio Gijuku Academic Development Funds. TY was supported
by JSPS through the Program for Leading Graduate School (ALPS).

##  Appendix A Entanglement spectrum and entanglement Hamiltonian

Here we describe details of the calculation of the ES $\xi_{\mathbf{k}}$ and
the coefficients $F_{\mathbf{k}}$ and $G_{\mathbf{k}}$ in the entanglement
Hamiltonian in Sec. 2.5.

For preparation, we calculate the phase and density correlators in the
oscillator ground state $|0^{\mathrm{osc}}\rangle$ of the total system. By
using Eqs. (21), (25), and (29), the phase and the density of the
spin-$\uparrow$ component are expressed in terms of the bogolon operators as

$\displaystyle\theta_{{\mathbf{k}},\uparrow}$
$\displaystyle=\frac{1}{2}\left(r_{\mathbf{k}},r_{\mathbf{k}}^{-1}\right)U({\mathbf{k}})\begin{pmatrix}\bar{\theta}_{{\mathbf{k}},1}\\\
\bar{\theta}_{{\mathbf{k}},2}\end{pmatrix}=\frac{1}{2}\sum_{j=1,2}R_{{\mathbf{k}},j}\bar{\theta}_{{\mathbf{k}},j}=\frac{1}{2\sqrt{2n}}\sum_{j=1,2}R_{{\mathbf{k}},j}\zeta_{{\mathbf{k}},j}^{-1}\left(\gamma_{{\mathbf{k}},j}+\gamma_{-{\mathbf{k}},j}^{\dagger}\right),$
(86a) $\displaystyle n_{{\mathbf{k}},\uparrow}$
$\displaystyle=\left(r_{\mathbf{k}}^{-1},r_{\mathbf{k}}\right)U({\mathbf{k}})\begin{pmatrix}\bar{n}_{{\mathbf{k}},1}\\\
\bar{n}_{{\mathbf{k}},2}\end{pmatrix}=\sum_{j=1,2}R_{{\mathbf{k}},3-j}\bar{n}_{{\mathbf{k}},j}=\frac{1}{\mathrm{i}}\sqrt{\frac{n}{2}}\sum_{j=1,2}R_{{\mathbf{k}},3-j}\zeta_{{\mathbf{k}},j}\left(\gamma_{{\mathbf{k}},j}-\gamma_{-{\mathbf{k}},j}^{\dagger}\right).$
(86b)

Here we introduce

$R_{{\mathbf{k}},1}:=r_{\mathbf{k}}\cos\frac{\chi_{\mathbf{k}}}{2}+\mathrm{i}r_{\mathbf{k}}^{-1}\sin\frac{\chi_{\mathbf{k}}}{2},~{}~{}~{}R_{{\mathbf{k}},2}:=r_{\mathbf{k}}^{-1}\cos\frac{\chi_{\mathbf{k}}}{2}+\mathrm{i}r_{\mathbf{k}}\sin\frac{\chi_{\mathbf{k}}}{2},$
(87)

which satisfy $R_{-{\mathbf{k}},j}=R_{{\mathbf{k}},j}^{\ast}~{}(j=1,2)$. As
$|0^{\mathrm{osc}}\rangle$ is the vacuum for bogolons, i.e.,
$\gamma_{{\mathbf{k}},\pm}|0^{\mathrm{osc}}\rangle=0$ for all
${\mathbf{k}}\neq{\mathbf{0}}$, the correlators of the operators in Eq. (86)
are calculated as

$\displaystyle\langle
0^{\mathrm{osc}}|\theta_{-{\mathbf{k}},\uparrow}\theta_{{\mathbf{k}},\uparrow}|0^{\mathrm{osc}}\rangle=\frac{1}{8n}\sum_{j=1,2}|R_{{\mathbf{k}},j}|^{2}\zeta_{{\mathbf{k}},j}^{-2}=\frac{1}{8n\zeta_{{\mathbf{k}},1}\zeta_{{\mathbf{k}},2}}\left(|R_{{\mathbf{k}},1}|^{2}\frac{\zeta_{{\mathbf{k}},2}}{\zeta_{{\mathbf{k}},1}}+|R_{{\mathbf{k}},2}|^{2}\frac{\zeta_{{\mathbf{k}},1}}{\zeta_{{\mathbf{k}},2}}\right),$
(88a) $\displaystyle\langle
0^{\mathrm{osc}}|n_{-{\mathbf{k}},\uparrow}n_{{\mathbf{k}},\uparrow}|0^{\mathrm{osc}}\rangle=\frac{n}{2}\sum_{j=1,2}|R_{{\mathbf{k}},3-j}|^{2}\zeta_{{\mathbf{k}},j}^{2}=\frac{n\zeta_{{\mathbf{k}},1}\zeta_{{\mathbf{k}},2}}{2}\left(|R_{{\mathbf{k}},1}|^{2}\frac{\zeta_{{\mathbf{k}},2}}{\zeta_{{\mathbf{k}},1}}+|R_{{\mathbf{k}},2}|^{2}\frac{\zeta_{{\mathbf{k}},1}}{\zeta_{{\mathbf{k}},2}}\right).$
(88b)

We now require that the correlators (42) obtained from the Gaussian ansatz
(38) be equal to the ones (88) obtained for the oscillator ground state
$|0^{\mathrm{osc}}\rangle$. We then find

$\displaystyle f_{\mathrm{B}}(\xi_{\mathbf{k}})+\frac{1}{2}$
$\displaystyle=\sqrt{\langle
0^{\mathrm{osc}}|\theta_{-{\mathbf{k}},\uparrow}\theta_{{\mathbf{k}},\uparrow}|0^{\mathrm{osc}}\rangle\langle
0^{\mathrm{osc}}|n_{-{\mathbf{k}},\uparrow}n_{{\mathbf{k}},\uparrow}|0^{\mathrm{osc}}\rangle}=\frac{1}{4}\left(|R_{{\mathbf{k}},1}|^{2}\frac{\zeta_{{\mathbf{k}},2}}{\zeta_{{\mathbf{k}},1}}+|R_{{\mathbf{k}},2}|^{2}\frac{\zeta_{{\mathbf{k}},1}}{\zeta_{{\mathbf{k}},2}}\right),$
(89a) $\displaystyle\sqrt{\frac{F_{\mathbf{k}}}{G_{\mathbf{k}}}}$
$\displaystyle=\frac{1}{n}\sqrt{\frac{\langle
0^{\mathrm{osc}}|n_{-{\mathbf{k}},\uparrow}n_{{\mathbf{k}},\uparrow}|0^{\mathrm{osc}}\rangle}{\langle
0^{\mathrm{osc}}|\theta_{-{\mathbf{k}},\uparrow}\theta_{{\mathbf{k}},\uparrow}|0^{\mathrm{osc}}\rangle}}=2\zeta_{{\mathbf{k}},1}\zeta_{{\mathbf{k}},2}.$
(89b)

In the long-wavelength limit $k\ell\ll 1$, we have

$\begin{split}&\zeta_{{\mathbf{k}},1}\approx\left[\frac{gD(\varphi)}{2\bar{g}_{\pm}}\right]^{1/4}\sqrt{k\ell},~{}~{}\zeta_{{\mathbf{k}},2}\approx\left[\frac{gC(\varphi)}{2\bar{g}_{\mp}}\right]^{1/4}k\ell,~{}~{}\cos\chi_{\mathbf{k}}\approx\mp
1,\\\
&|R_{{\mathbf{k}},1}|^{2}=r_{\mathbf{k}}^{2}\frac{1+\cos\chi_{\mathbf{k}}}{2}+r_{\mathbf{k}}^{-2}\frac{1-\cos\chi_{\mathbf{k}}}{2}\approx
r_{\mathbf{0}}^{\mp 2}=\sqrt{\frac{\bar{g}_{\mp}}{\bar{g}_{\pm}}},\\\
&|R_{{\mathbf{k}},2}|^{2}=r_{\mathbf{k}}^{-2}\frac{1+\cos\chi_{\mathbf{k}}}{2}+r_{\mathbf{k}}^{2}\frac{1-\cos\chi_{\mathbf{k}}}{2}\approx
r_{\mathbf{0}}^{\pm 2}=\sqrt{\frac{\bar{g}_{\pm}}{\bar{g}_{\mp}}},\end{split}$
(90)

where Eq. (32) is used. Then, $f_{\mathrm{B}}(\xi_{\mathbf{k}})$ and
$\sqrt{F_{\mathbf{k}}/G_{\mathbf{k}}}$ in Eq. (89) are calculated as

$f_{\mathrm{B}}(\xi_{\mathbf{k}})\approx\frac{1}{4}\left[\frac{\bar{g}_{\pm}D(\varphi)}{\bar{g}_{\mp}C(\varphi)}\right]^{1/4}(k\ell)^{-1/2},~{}~{}~{}\sqrt{\frac{F_{\mathbf{k}}}{G_{\mathbf{k}}}}\approx\left[\frac{4g^{2}C(\varphi)D(\varphi)}{\bar{g}_{+}\bar{g}_{-}}\right]^{1/4}(k\ell)^{3/2}.$
(91)

from which we obtain Eq. (43).

##  Appendix B Phase correlation function

Here we describe the derivation of the phase correlation function (54). Using
Eq. (52), this correlation function is expressed as

$\langle\left[\theta_{\uparrow}({\mathbf{r}})-\theta_{\uparrow}({\mathbf{0}})\right]^{2}\rangle\approx\frac{2}{nF\ell^{2}}\left[G({\mathbf{0}};\alpha)-G({\mathbf{r}};\alpha)\right].$
(92)

Here, $G({\mathbf{r}};\alpha)$ is Green’s function for the 2D Poisson equation

$G({\mathbf{r}};\alpha)=\frac{1}{A}\sum_{{\mathbf{k}}\neq{\mathbf{0}}}\frac{1}{k^{2}}{\mathrm{e}}^{-\alpha
k+\mathrm{i}{\mathbf{k}}\cdot{\mathbf{r}}}=\frac{1}{A}\sum_{{\mathbf{k}}\neq{\mathbf{0}}}\frac{1}{k^{2}}{\mathrm{e}}^{-\alpha
k+\mathrm{i}kr\cos\theta},$ (93)

where we express the wave vector ${\mathbf{k}}$ using the polar coordinate
$(k,\theta)$ with $\theta$ being the angle relative to ${\mathbf{r}}$.
Although the logarithmic behavior of Green’s function for the 2D Poisson
equation is known, we derive it within the present regularization scheme using
the convergence factor ${\mathrm{e}}^{-\alpha k}$.

By differentiating Eq. (93) with respect to $r$, we have

$\begin{split}-\frac{\partial}{\partial
r}G({\mathbf{r}};\alpha)&=\frac{1}{A}\sum_{{\mathbf{k}}\neq{\mathbf{0}}}\frac{-\mathrm{i}\cos\theta}{k}{\mathrm{e}}^{-\alpha
k+\mathrm{i}kr\cos\theta}=\frac{1}{(2\pi)^{2}}\int_{0}^{2\pi}\\!\\!\\!\\!\mathrm{d}\theta\int_{0}^{\infty}\\!\\!\\!\\!\mathrm{d}k~{}(-\mathrm{i}\cos\theta){\mathrm{e}}^{-k(\alpha-\mathrm{i}r\cos\theta)}\\\
&=\frac{1}{(2\pi)^{2}}\int_{0}^{2\pi}\mathrm{d}\theta\frac{-\mathrm{i}\cos\theta}{\alpha-\mathrm{i}r\cos\theta}.\end{split}$
(94)

With a change of the integration variable as
$z={\mathrm{e}}^{\mathrm{i}\theta}$, the last integral can be written as a
contour integral along the unit circle:

$-\frac{\partial}{\partial
r}G({\mathbf{r}};\alpha)=\frac{1}{(2\pi)^{2}}\oint\frac{\mathrm{d}z}{\mathrm{i}z}\frac{z^{2}+1}{r(z^{2}+1)+2\mathrm{i}\alpha
z}=\frac{1}{(2\pi)^{2}\mathrm{i}r}\oint\mathrm{d}z\frac{z^{2}+1}{z(z-z_{+})(z-z_{-})},$
(95)

where $z_{\pm}=\mathrm{i}(-\alpha\pm\sqrt{r^{2}+\alpha^{2}})/r$ are the
locations of poles. Since $|z_{+}|<1<|z_{-}|$, the integral picks up the
residues at $z=0$ and $z_{+}$, leading to

$-\frac{\partial}{\partial r}G({\mathbf{r}};\alpha)=\frac{1}{2\pi
r}\left[\frac{1}{z_{+}z_{-}}+\frac{z_{+}^{2}+1}{z_{+}(z_{+}-z_{-})}\right]=\frac{1}{2\pi
r}\left(1-\frac{\alpha}{\sqrt{r^{2}+\alpha^{2}}}\right).$ (96)

Therefore, $G({\mathbf{0}};\alpha)-G({\mathbf{r}};\alpha)$ in Eq. (92) can be
calculated as

$\begin{split}G({\mathbf{0}};\alpha)-G({\mathbf{r}};\alpha)&=\frac{1}{2\pi}\lim_{a_{0}\to
0}\int_{a_{0}}^{r}\mathrm{d}r^{\prime}\left(\frac{1}{r^{\prime}}-\frac{\alpha}{r^{\prime}\sqrt{{r^{\prime}}^{2}+\alpha^{2}}}\right)\\\
&=\frac{1}{2\pi}\lim_{a_{0}\to
0}\left(\ln\frac{r}{a_{0}}+\mathrm{arsinh}~{}\frac{\alpha}{r}-\mathrm{arsinh}~{}\frac{\alpha}{a_{0}}\right)=\frac{1}{2\pi}\left(\ln\frac{r}{2\alpha}+\mathrm{arsinh}~{}\frac{\alpha}{r}\right).\end{split}$
(97)

By substituting this into Eq. (92) and taking the limit $r\gg\alpha$, we
obtain Eq. (54).

## References

## References

  * [1] Dalibard J, Gerbier F, Juzeli$\bar{\mathrm{u}}$nas G and Öhberg P 2011 Rev. Mod. Phys. 83(4) 1523
  * [2] Goldman N, Juzeli$\bar{\mathrm{u}}$nas G, Öhberg P and Spielman I B 2014 Reports on Progress in Physics 77 126401
  * [3] Zhang S L and Zhou Q 2017 Journal of Physics B: Atomic, Molecular and Optical Physics 50 222001
  * [4] Aidelsburger M, Nascimbene S and Goldman N 2018 Comptes Rendus Physique 19 394 – 432
  * [5] Galitski V, Juzeliūnas G and Spielman I B 2019 Physics Today 72 38–44
  * [6] Cooper N R 2008 Advances in Physics 57 539
  * [7] Fetter A L 2009 Rev. Mod. Phys. 81(2) 647
  * [8] Lin Y J, Compton R L, Jimenez-Garcia K, Porto J V and Spielman I B 2009 Nature 462 628
  * [9] Abrikosov A A 1957 Sov. Phys. JETP 5 1174 [Zh. Eksp. Teor. Fiz.32,1442(1957)]
  * [10] Abo-Shaeer J R, Raman C, Vogels J M and Ketterle W 2001 Science 292 476
  * [11] Engels P, Coddington I, Haljan P C and Cornell E A 2002 Phys. Rev. Lett. 89(10) 100403
  * [12] Schweikhard V, Coddington I, Engels P, Mogendorff V P and Cornell E A 2004 Phys. Rev. Lett. 92(4) 040404
  * [13] Tkachenko V K 1966 Soviet Journal of Experimental and Theoretical Physics 22 1282 [Zh. Eksp. Teor. Fiz.49,1875(1966)]
  * [14] Tkachenko V K 1966 Soviet Journal of Experimental and Theoretical Physics 23 1049 [Zh. Eksp. Teor. Fiz.50,1573(1966)]
  * [15] Tkachenko V K 1969 Soviet Journal of Experimental and Theoretical Physics 29 945 [Zh. Eksp. Teor. Fiz.56,1763(1969)]
  * [16] Sonin E B 1987 Rev. Mod. Phys. 59(1) 87
  * [17] Coddington I, Engels P, Schweikhard V and Cornell E A 2003 Phys. Rev. Lett. 91(10) 100402
  * [18] Sinova J, Hanna C B and MacDonald A H 2002 Phys. Rev. Lett. 89(3) 030403
  * [19] Baym G 2003 Phys. Rev. Lett. 91(11) 110402
  * [20] Baym G 2004 Phys. Rev. A 69(4) 043618
  * [21] Matveenko S I and Shlyapnikov G V 2011 Phys. Rev. A 83(3) 033604
  * [22] Kwasigroch M P and Cooper N R 2012 Phys. Rev. A 86(6) 063618
  * [23] Watanabe H and Murayama H 2013 Phys. Rev. Lett. 110(18) 181601
  * [24] Ho T L 2001 Phys. Rev. Lett. 87(6) 060403
  * [25] Sonin E B 2005 Phys. Rev. A 72(2) 021606
  * [26] Wilkin N K, Gunn J M F and Smith R A 1998 Phys. Rev. Lett. 80(11) 2265
  * [27] Cooper N R, Wilkin N K and Gunn J M F 2001 Phys. Rev. Lett. 87(12) 120405
  * [28] Cooper N R in Fractional Quantum Hall Effects, edited by B. I. Halperin and J. K. Jain (World Scientific, 2020)
  * [29] Cooper N R and Rezayi E H 2007 Phys. Rev. A 75(1) 013627
  * [30] Liu Z, Guo H L, Vedral V and Fan H 2011 Phys. Rev. A 83(1) 013620
  * [31] Schweikhard V, Coddington I, Engels P, Tung S and Cornell E A 2004 Phys. Rev. Lett. 93(21) 210403
  * [32] Lin Y J, Jiménez-García K and Spielman I B 2011 Nature 471 83
  * [33] Zhai H 2012 International Journal of Modern Physics B 26 1230001
  * [34] Wu Z, Zhang L, Sun W, Xu X T, Wang B Z, Ji S C, Deng Y, Chen S, Liu X J and Pan J W 2016 Science 354 83–88
  * [35] Beeler M C, Williams R A, Jimenez-Garcia K, LeBlanc L J, Perry A R and Spielman I B 2013 Nature 498 201 letter
  * [36] Mueller E J and Ho T L 2002 Phys. Rev. Lett. 88(18) 180403
  * [37] Kasamatsu K, Tsubota M and Ueda M 2003 Phys. Rev. Lett. 91(15) 150406
  * [38] Kasamatsu K, Tsubota M and Ueda M 2005 International Journal of Modern Physics B 19 1835
  * [39] Mingarelli L, Keaveny E E and Barnett R 2018 Phys. Rev. A 97(4) 043622
  * [40] Liu X J, Liu X, Kwek L C and Oh C H 2007 Phys. Rev. Lett. 98(2) 026602
  * [41] Fialko O, Brand J and Zülicke U 2014 New Journal of Physics 16 025006
  * [42] Furukawa S and Ueda M 2014 Phys. Rev. A 90(3) 033602
  * [43] Girvin S M and MacDonald A H in Perspectives in Quantum Hall Effects, edited by S. Das Sarma and A. Pinczuk (Wiley, New York, 1997)
  * [44] Bernevig B A and Zhang S C 2006 Phys. Rev. Lett. 96(10) 106802
  * [45] Keçeli M and Oktel M O 2006 Phys. Rev. A 73(2) 023611
  * [46] Yoshino T, Furukawa S, Higashikawa S and Ueda M 2019 New Journal of Physics 21 015001
  * [47] Furukawa S and Ueda M 2013 Phys. Rev. Lett. 111(9) 090401
  * [48] Wu Y H and Jain J K 2013 Phys. Rev. B 87(24) 245123
  * [49] Regnault N and Senthil T 2013 Phys. Rev. B 88(16) 161106
  * [50] Wu Y H and Jain J K 2015 Phys. Rev. A 91(6) 063623
  * [51] Geraedts S D, Repellin C, Wang C, Mong R S K, Senthil T and Regnault N 2017 Phys. Rev. B 96(7) 075148
  * [52] Furukawa S and Ueda M 2017 Phys. Rev. A 96(5) 053626
  * [53] Woo S J, Choi S, Baksmaty L O and Bigelow N P 2007 Phys. Rev. A 75(3) 031604
  * [54] Barnett R, Refael G, Porter M A and Büchler H P 2008 New Journal of Physics 10 043030
  * [55] Mason P and Aftalion A 2011 Phys. Rev. A 84(3) 033611
  * [56] Aftalion A, Mason P and Wei J 2012 Phys. Rev. A 85(3) 033614
  * [57] Kuopanportti P, Huhtamäki J A M and Möttönen M 2012 Phys. Rev. A 85(4) 043613
  * [58] Kumar R K, Tomio L, Malomed B A and Gammal A 2017 Phys. Rev. A 96(6) 063624
  * [59] Mingarelli L and Barnett R 2019 Phys. Rev. Lett. 122(4) 045301
  * [60] Lee T D, Huang K and Yang C N 1957 Phys. Rev. 106(6) 1135–1145
  * [61] Yoshino T, Furukawa S and Ueda M 2021 Phys. Rev. A 103(4) 043321
  * [62] Chen X and Fradkin E 2013 J. Stat. Mech. Theor. Exp. 2013 P08013
  * [63] Lundgren R, Fuji Y, Furukawa S and Oshikawa M 2013 Phys. Rev. B 88(24) 245137
  * [64] Moroz S, Hoyos C, Benzoni C and Son D T 2018 SciPost Phys. 5(4) 39
  * [65] Moroz S and Son D T 2019 Phys. Rev. Lett. 122(23) 235301
  * [66] Aftalion A, Blanc X and Dalibard J 2005 Phys. Rev. A 71(2) 023611
  * [67] Peschel I 2003 J. Physics A: Mathematical and General 36 L205–L208
  * [68] Peschel I and Eisler V 2009 Journal of Physics A: Mathematical and Theoretical 42 504003
  * [69] Metlitski M A and Grover T 2011 Entanglement entropy of systems with spontaneously broken continuous symmetry (Preprint arXiv:1112.5166)
  * [70] Pethick C J and Smith H 2008 Bose–Einstein Condensation in Dilute Gases 2nd ed (Cambridge University Press)
  * [71] Penrose O and Onsager L 1956 Phys. Rev. 104(3) 576–584
  * [72] Yang C N 1962 Rev. Mod. Phys. 34(4) 694–704
  * [73] Rashba E I, Zhukov L E and Efros A L 1997 Phys. Rev. B 55(8) 5306
  * [74] Burkov A A 2010 Phys. Rev. B 81(12) 125111
  * [75] Panfilov I, Patri A, Yang K and Burkov A A 2016 Phys. Rev. B 93(12) 125126
  * [76] Furukawa S and Kim Y B 2011 Phys. Rev. B 83(8) 085112
  * [77] Xu C 2011 Phys. Rev. B 84(12) 125119
  * [78] Mollabashi A, Shiba N and Takayanagi T 2014 Journal of High Energy Physics 2014 185
  * [79] Repellin C, Bernevig B A and Regnault N 2014 Phys. Rev. B 90(24) 245401
  * [80] Neupert T, Santos L, Ryu S, Chamon C and Mudry C 2011 Phys. Rev. B 84(16) 165107
  * [81] Li W, Sheng D N, Ting C S and Chen Y 2014 Phys. Rev. B 90(8) 081102
  * [82] Chen H and Yang K 2012 Phys. Rev. B 85(19) 195113
  * [83] Ghaemi P, Cayssol J, Sheng D N and Vishwanath A 2012 Phys. Rev. Lett. 108(26) 266801