Vortex lattices in binary Bose–Einstein condensates: collective modes, quantum fluctuations, and intercomponent entanglement

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

$$\begin{equation}\begin{aligned}\hfill \mathcal{L}& =\sum\limits _{\alpha ={\uparrow},{\downarrow}}\left[\frac{\text{i}\hslash }{2}({\psi }_{\alpha }^{\ast }{\dot{\psi }}_{\alpha }-{\dot{\psi }}_{\alpha }^{\ast }{\psi }_{\alpha })-\frac{1}{2M}\vert (-\text{i}\hslash \nabla -q{\mathbf{A}}_{\alpha }){\psi }_{\alpha }{\vert }^{2}\right]\hfill \\ \hfill & \quad -\sum\limits _{\alpha ,\beta ={\uparrow},{\downarrow}}\frac{{g}_{\alpha \beta }}{2}\vert {\psi }_{\alpha }{\vert }^{2}\vert {\psi }_{\beta }{\vert }^{2},\hfill \end{aligned}\end{equation} \tag{ 1 }$$

where ψ_α (r, t) is the bosonic field for the spin-α component (with r = (x, y) being the 2D coordinate), and M and q are the mass and the fictitious charge of an atom. This fictitious charge specifies the coupling strength between atoms and the synthetic gauge fields. The gauge field A_α (r) for spin-α bosons is given by

$$\begin{equation}{\mathbf{A}}_{\alpha }=\frac{{B}_{\alpha }}{2}{\mathbf{e}}_{z}\times \mathbf{r}={{\epsilon}}_{\alpha }\frac{B}{2}(-y,x),\end{equation} \tag{ 2 }$$

where we assume qB > 0 and _↑ = _↓ = 1  (_↑ = −_↓ = 1) for parallel (antiparallel) fields. The number of magnetic flux quanta piercing each component is given by N_v = qBA/(2πℏ) = A/(2πℓ²), where A is the area of the system and $\ell {:=}\sqrt{\hslash /qB}$ is the magnetic length. In the Lagrangian (1), the numbers of spin-↑ and ↓ atoms, N_↑ and N_↓, are separately conserved. We introduce the total filling factor ν := N/N_v, where N := N_↑ + N_↓ is the total number of atoms.

We assume a contact interaction between atoms. For simplicity, we set g_↑↑ = g_↓↓ ≡ g > |g_↑↓| and N_↑ = N_↓ 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_αβ |n, is much smaller than the Landau-level spacing ℏω_c. Here, n := N_↑/A = N_↓/A is the average density of ↑ or ↓ atoms, and ω_c := qB/M is the cyclotron frequency.

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

$$\begin{equation}\begin{aligned}\hfill \mathcal{L}& =\sum\limits _{\alpha }\left[\hslash {n}_{\alpha }{\dot{\theta }}_{\alpha }-\frac{{n}_{\alpha }}{2M}{(\hslash \nabla {\theta }_{\alpha }+q{\mathbf{A}}_{\alpha })}^{2}-\frac{{\hslash }^{2}{(\nabla {n}_{\alpha })}^{2}}{8{n}_{\alpha }M}\right]\hfill \\ \hfill & \quad -\sum\limits _{\alpha ,\beta }\frac{{g}_{\alpha \beta }}{2}{n}_{\alpha }{n}_{\beta }.\hfill \end{aligned}\end{equation} \tag{ 3 }$$

In the presence of vortices, the phase variables {θ_α (r, t)} involve singularities. This motivates us to decompose θ_α into regular and singular contributions as θ_α = θ_reg,α + θ_sing,α . We also introduce the displacement u_α (r, t) of a vortex from the equilibrium position. The derivatives of the singular part θ_sing,α of the phase can be related to the displacement field u_α as [23, 64]

$$\begin{equation}\hslash {\dot{\theta }}_{\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{g},\alpha }=-\frac{q{B}_{\alpha }}{2}{({\mathbf{u}}_{\alpha }\times {\dot{\mathbf{u}}}_{\alpha })}_{z},\end{equation} \tag{ 4a }$$

$$\begin{equation}\hslash \nabla {\theta }_{\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{g},\alpha }+q{\mathbf{A}}_{\alpha }=q{B}_{\alpha }{\mathbf{e}}_{z}\times {\mathbf{u}}_{\alpha }-\frac{q{B}_{\alpha }}{2}{{\epsilon}}_{ij}{u}_{\alpha }^{i}\nabla {u}_{\alpha }^{j},\end{equation} \tag{ 4b }$$

where _ij is an antisymmetric tensor with _xy = −_yx = +1. The displacement u_α (r, t) also results in a change in the elastic energy $\int {\mathrm{d}}^{2}\mathbf{r}\,{\mathcal{E}}_{\mathrm{e}\mathrm{l}}({\mathbf{u}}_{\alpha },{\partial }_{i}{\mathbf{u}}_{\alpha })$, whose explicit form will be given in section 2.3. The Lagrangian density can then be expressed in terms of {θ_reg,α , u_α , n_α } as

$$\begin{equation}\begin{aligned}\hfill \mathcal{L}=\;& \sum\limits _{\alpha }\left[\hslash {n}_{\alpha }{\dot{\theta }}_{\mathrm{r}\mathrm{e}\mathrm{g},\alpha }-\frac{q{B}_{\alpha }{n}_{\alpha }}{2}{({\mathbf{u}}_{\alpha }\times {\dot{\mathbf{u}}}_{\alpha })}_{z}\right.\hfill \\ \hfill & -\frac{{n}_{\alpha }}{2M}{\left(\hslash \nabla {\theta }_{\text{reg},\alpha }+q{B}_{\alpha }{\mathbf{e}}_{z}\times {\mathbf{u}}_{\alpha }-\frac{q{B}_{\alpha }}{2}{{\epsilon}}_{ij}{u}_{\alpha }^{i}\nabla {u}_{\alpha }^{j}\right)}^{2}\hfill \\ \hfill & \left.-\frac{{\hslash }^{2}{(\nabla {n}_{\alpha })}^{2}}{8{n}_{\alpha }M}\right]-\sum\limits _{\alpha ,\beta }\frac{{\bar{g}}_{\alpha \beta }}{2}{n}_{\alpha }{n}_{\beta }-{\mathcal{E}}_{\mathrm{e}\mathrm{l}}.\hfill \end{aligned}\end{equation} \tag{ 5 }$$

Henceforth, we ignore the term $-\frac{q{B}_{\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 ∇θ_reg,α and u_α . We also omit the subscript 'reg' in θ_reg,α .

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_α (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_α (r, t) should now be understood as the density averaged over the unit cell containing r. The renormalized coupling constants ${\bar{g}}_{\alpha \beta }$ can then be introduced through the relation

$$\begin{equation}{\int }_{\text{u.c.}}{\mathrm{d}}^{2}{\mathbf{r}}^{\prime }\,{g}_{\alpha \beta }\vert {\psi }_{\alpha }({\mathbf{r}}^{\prime },t){\vert }^{2}\vert {\psi }_{\beta }({\mathbf{r}}^{\prime },t){\vert }^{2}=2\pi {\ell }^{2}{\bar{g}}_{\alpha \beta }{n}_{\alpha }(\mathbf{r},t){n}_{\beta }(\mathbf{r},t),\end{equation} \tag{ 6 }$$

where |ψ_α |² and n_α are the original and coarse-grained densities of the spin-α atoms, respectively, and the integration on the left-hand side is taken over the unit cell (with the area 2πℓ²) that contains r. As explained later in this section and in section 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 equation (6), a positive (negative) correlation between the density fluctuations |ψ_α |² − n_α and |ψ_β |² − n_β leads to an enhanced (reduced) renormalized coupling $\vert {\bar{g}}_{\alpha \beta }\vert > \vert {g}_{\alpha \beta }\vert$ $(\vert {\bar{g}}_{\alpha \beta }\vert < \vert {g}_{\alpha \beta }\vert )$. In particular, the intracomponent coupling g is always enhanced by the renormalization while the intercomponent repulsion g_↑↓ > 0 (attraction g_↑↓ < 0) is reduced (enhanced) by the renormalization owing to the separation (overlap) of vortices between the two components; this can be confirmed in figure 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 references [7, 18, 66]. At high filling factors ν ≫ 1 and for not too large a number of vortices N_v, we can assume that the condensates are only weakly depleted; we therefore have |n_α (r, t) − n| ≪ n for the coarse-grained density n_α .

Zoom In Zoom Out Reset image size

Because the displacement fields {u_α } involve the mass term $-{\mathbf{u}}_{\alpha }^{2}$ in equation (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 {u_α }:

$$\begin{equation}\begin{aligned}\hfill & {\mathbf{u}}_{\alpha }-{{\epsilon}}_{\alpha }{\ell }^{2}{\mathbf{e}}_{z}\times \nabla {\theta }_{\alpha }-\frac{{{\epsilon}}_{\alpha }}{{\omega }_{\text{c}}}{\mathbf{e}}_{z}\times {\dot{\mathbf{u}}}_{\alpha }\hfill \\ \hfill & \hspace{25.0pt}+\frac{{\ell }^{2}}{n\hslash {\omega }_{\text{c}}}\left[\frac{\partial {\mathcal{E}}_{\mathrm{e}\mathrm{l}}}{\partial {\mathbf{u}}_{\alpha }}-{\partial }_{j}\left(\frac{\partial {\mathcal{E}}_{\mathrm{e}\mathrm{l}}}{\partial \left({\partial }_{j}{\mathbf{u}}_{\alpha }\right)}\right)\right]=0,\hfill \end{aligned}\end{equation} \tag{ 7 }$$

where we have made the approximation n_α ≈ n. We can ignore the third and fourth terms on the left-hand side in the LLL approximation with $\hslash \omega ,{\mathcal{E}}_{\mathrm{e}\mathrm{l}}/n\ll \hslash {\omega }_{\text{c}}$, where ω is the frequency of our concern. Introducing u_± := u_↑ ± u_↓ and θ_± := θ_↑ ± θ_↓, we can rewrite equation (7) as

$$\begin{equation}{\mathbf{u}}_{\pm }=\left\{\begin{aligned}\hfill & {\ell }^{2}{\mathbf{e}}_{z}\times \nabla {\theta }_{\pm }\hspace{25.0pt}(\text{parallel}\;\text{fields});\hfill \\ \hfill & {\ell }^{2}{\mathbf{e}}_{z}\times \nabla {\theta }_{\mp }\hspace{25.0pt}(\text{antiparallel}\;\text{fields}).\hfill \end{aligned}\right.\hspace{25.0pt}\end{equation} \tag{ 8 }$$

These relations indicate that the (anti)symmetric movement of vortices is coupled to the (anti)symmetric component of the phase variables for parallel fields, while they are coupled in a crossed manner for antiparallel fields. By substituting equation (8) into equation (5), we obtain the effective Lagrangian in terms of θ_± and their conjugate momenta ℏn_± = ℏ(n_↑ ± n_↓)/2. The Hamiltonian density is then obtained as

$$\begin{equation}\mathcal{H}=\sum\limits _{\nu =\pm }\left[{\bar{g}}_{\nu }{n}_{\nu }^{2}+\frac{{\hslash }^{2}{(\nabla {n}_{\nu })}^{2}}{4nM}\right]+{\mathcal{E}}_{\mathrm{e}\mathrm{l}},\quad {\bar{g}}_{\pm }{:=}\bar{g}\pm {\bar{g}}_{{\uparrow}{\downarrow}},\end{equation} \tag{ 9 }$$

where ${\mathcal{E}}_{\mathrm{e}\mathrm{l}}$ is expressed in terms of the phase variables θ_± by using equation (8). The theory can be quantized by requiring the canonical commutation relations

$$\begin{equation}[{\theta }_{\nu }(\mathbf{r}),{n}_{{\nu }^{\prime }}({\mathbf{r}}^{\prime })]=\mathrm{i}{\delta }_{\nu {\nu }^{\prime }}\delta (\mathbf{r}-{\mathbf{r}}^{\prime })\quad (\nu ,{\nu }^{\prime }=\pm ).\end{equation} \tag{ 10 }$$

In the present coarse-grained description, the mean-field ground state corresponds to the uniform state with n₊(r) = n, n₋(r) = 0, and ∇θ_±(r) = 0. Therefore, using equation (9), we obtain the mean-field ground-state energy density as

$$\begin{equation}\frac{{E}_{\text{GS}}^{\text{MF}}}{A}={\bar{g}}_{+}{n}^{2}=(\bar{g}+{\bar{g}}_{{\uparrow}{\downarrow}}){n}^{2}=\frac{1}{2}\sum\limits _{\alpha ,\beta }{\bar{g}}_{\alpha \beta }{n}^{2}.\end{equation} \tag{ 11 }$$

In contrast, the same energy density obtained by the microscopic calculation (the first term of equation (67) shown later) has the form ${E}_{\text{GS}}^{\text{MF}}/A=(\beta g+{\beta }_{{\uparrow}{\downarrow}}{g}_{{\uparrow}{\downarrow}}){n}^{2}$, where β and β_↑↓ 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}}$.

The expression of the elastic energy density ${\mathcal{E}}_{\mathrm{e}\mathrm{l}}$ has been determined in references [45, 46], and we summarize it in the following. We first note that the elastic energy must be invariant under a constant change in u_α (r, t), i.e., translation of the lattices. Therefore, to the leading order in the derivative expansion, ${\mathcal{E}}_{\mathrm{e}\mathrm{l}}$ should be a function of ∂_i u₊ (i = x, y) and u₋, resulting in the decomposition

$$\begin{equation}{\mathcal{E}}_{\mathrm{e}\mathrm{l}}={\mathcal{E}}_{\mathrm{e}\mathrm{l}}^{(+)}({\partial }_{i}{\mathbf{u}}_{+})+{\mathcal{E}}_{\mathrm{e}\mathrm{l}}^{(-)}({\mathbf{u}}_{-})+{\mathcal{E}}_{\mathrm{e}\mathrm{l}}^{(+-)}({\partial }_{i}{\mathbf{u}}_{+},{\mathbf{u}}_{-}).\end{equation} \tag{ 12 }$$

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

$$\begin{equation}{w}_{1}{:=}{\partial }_{x}{u}_{+}^{x}-{\partial }_{y}{u}_{+}^{y},\quad {w}_{2}{:=}{\partial }_{y}{u}_{+}^{x}+{\partial }_{x}{u}_{+}^{y}.\end{equation} \tag{ 13 }$$

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

$$\begin{equation}{\mathcal{E}}_{\mathrm{e}\mathrm{l}}^{(+)}({\partial }_{i}{\mathbf{u}}_{+})=\frac{g{n}^{2}}{2}\left({C}_{1}{{w}_{1}}^{2}+{C}_{2}{{w}_{2}}^{2}+{C}_{3}{w}_{1}{w}_{2}\right),\end{equation} \tag{ 14a }$$

$$\begin{equation}{\mathcal{E}}_{\mathrm{e}\mathrm{l}}^{(-)}({\mathbf{u}}_{-})=\frac{g{n}^{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],\end{equation} \tag{ 14b }$$

$$\begin{equation}{\mathcal{E}}_{\mathrm{e}\mathrm{l}}^{(+-)}({\partial }_{i}{\mathbf{u}}_{+},{\mathbf{u}}_{-})=\frac{g{n}^{2}}{2\ell }{F}_{1}\left({w}_{1}{u}_{-}^{y}+{w}_{2}{u}_{-}^{x}\right).\end{equation} \tag{ 14c }$$

For each of the vortex-lattice structures in figures 1(a)–(e), the dimensionless elastic constants {C₁, C₂, C₃, D₁, D₂, D₃, F₁} satisfy

$$\begin{equation}\begin{aligned}\hfill & (\text{a})\quad {C}_{1}={C}_{2}\equiv C > 0,\quad {D}_{1}={D}_{2}\equiv D > 0,\quad {C}_{3}={D}_{3}={F}_{1}=0;\hfill \\ \hfill & (\text{b})\quad {C}_{1}={C}_{2}\equiv C > 0,\quad {D}_{1}={D}_{2}\equiv D > 0,\quad {C}_{3}={D}_{3}=0,\quad {F}_{1}\ne 0;\hfill \\ \hfill & (\text{c})\quad {C}_{1},{C}_{2},{D}_{1},{D}_{2} > 0,\quad {C}_{3},{D}_{3}\ne 0,\quad {F}_{1}=0;\hfill \\ \hfill & (\mathrm{d})\quad {C}_{1},{C}_{2} > 0,\quad {D}_{1}={D}_{2}\equiv D > 0,\quad {C}_{3}={D}_{3}={F}_{1}=0;\hfill \\ \hfill & (\text{e})\quad {C}_{1},{C}_{2} > 0,\quad {D}_{1},{D}_{2} > 0,\quad {C}_{3}={D}_{3}={F}_{1}=0.\hfill \end{aligned}\end{equation} \tag{ 15 }$$

Keçeli and Oktel [45] have determined the constants {C₁, C₂, C₃, D₁, D₂, D₃} 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₁ term for interlaced triangular vortex lattices (b), which was missed in reference [45]. For ν ≫ 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 section 4.2, we will determine all the elastic constants above as a function of g_↑↓/g by using the data of the excitation spectra.

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

$$\begin{equation}{\theta }_{\nu }(\mathbf{r})=\frac{1}{\sqrt{A}}\sum\limits _{\mathbf{k}}{\theta }_{\mathbf{k},\nu }\,{\mathrm{e}}^{\mathrm{i}\mathbf{k}\cdot \mathbf{r}},\end{equation} \tag{ 16a }$$

$$\begin{equation}{n}_{\nu }(\mathbf{r})=\frac{1}{\sqrt{A}}\sum\limits _{\mathbf{k}}{n}_{\mathbf{k},\nu }\,{\mathrm{e}}^{\mathrm{i}\mathbf{k}\cdot \mathbf{r}}\quad (\nu =\pm ),\end{equation} \tag{ 16b }$$

where the Fourier components satisfy

$$\begin{equation}[{\theta }_{\mathbf{k},\nu },{n}_{-{\mathbf{k}}^{\prime },{\nu }^{\prime }}]=\mathrm{i}{\delta }_{\nu {\nu }^{\prime }}{\delta }_{\mathbf{k}{\mathbf{k}}^{\prime }},\end{equation} \tag{ 17a }$$

$$\begin{equation}{\theta }_{\mathbf{k},\nu }^{{\dagger}}={\theta }_{-\mathbf{k},\nu },\quad {n}_{\mathbf{k},\nu }^{{\dagger}}={n}_{-\mathbf{k},\nu }\quad (\nu ,{\nu }^{\prime }=\pm ).\end{equation} \tag{ 17b }$$

We note that the k = 0 component n_0,± 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{equation}\begin{aligned}\hfill H=\;& \frac{n}{2}\sum\limits _{\mathbf{k}}\left(\begin{matrix}\hfill {\theta }_{-\mathbf{k},+}\hfill & \hfill {\theta }_{-\mathbf{k},-}\hfill \end{matrix}\right)\left(\begin{matrix}\hfill {{\Gamma}}_{\pm }(\mathbf{k})\hfill & \hfill \pm \mathrm{i}{\Gamma}(\mathbf{k})\hfill \\ \hfill \mp \mathrm{i}{\Gamma}(\mathbf{k})\hfill & \hfill {{\Gamma}}_{\mp }(\mathbf{k})\hfill \end{matrix}\right)\left(\begin{matrix}\hfill {\theta }_{\mathbf{k},+}\hfill \\ \hfill {\theta }_{\mathbf{k},-}\hfill \end{matrix}\right)\hfill \\ \hfill & +\frac{1}{2n}\sum\limits _{\mathbf{k}}\sum\limits _{\nu =\pm }{e}_{\mathbf{k},\nu }{n}_{-\mathbf{k},\nu }{n}_{\mathbf{k},\nu },\hfill \end{aligned}\end{equation} \tag{ 18 }$$

where ⁶

$$\begin{equation}{e}_{\mathbf{k},\nu }{:=}2{\bar{g}}_{\nu }n+\frac{{\hslash }^{2}{\mathbf{k}}^{2}}{2M},\end{equation} \tag{ 19a }$$

$$\begin{align}\hfill {{\Gamma}}_{+}(\mathbf{k})& {:=}gn{\ell }^{4}\left[{C}_{1}{\left(2{k}_{x}{k}_{y}\right)}^{2}+{C}_{2}{\left({k}_{x}^{2}-{k}_{y}^{2}\right)}^{2}\right.\hfill \\ \hfill & \quad \left.-{C}_{3}(2{k}_{x}{k}_{y})\left({k}_{x}^{2}-{k}_{y}^{2}\right)\right],\hfill \end{align} \tag{ 19b }$$

$$\begin{equation}{{\Gamma}}_{-}(\mathbf{k}){:=}gn{\ell }^{2}\left({D}_{1}{k}_{y}^{2}+{D}_{2}{k}_{x}^{2}-{D}_{3}{k}_{x}{k}_{y}\right),\end{equation} \tag{ 19c }$$

$$\begin{equation}{\Gamma}(\mathbf{k}){:=}\frac{1}{2}gn{\ell }^{3}{F}_{1}\left[(2{k}_{x}{k}_{y}){k}_{x}+({k}_{x}^{2}-{k}_{y}^{2}){k}_{y}\right].\end{equation} \tag{ 19d }$$

In equation (18), the upper and lower signs correspond to the cases of parallel and antiparallel fields, respectively ⁷ .

It is useful to decompose the Hamiltonian (18) into the zero-mode (k = 0) and oscillator-mode (k ≠ 0) parts. First, the zero-mode part is given by

$$\begin{equation}{H}^{\mathrm{z}\mathrm{e}\mathrm{r}\mathrm{o}}=\sum\limits _{\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}.\end{equation} \tag{ 20 }$$

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

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

$$\begin{equation}{\tilde{\theta }}_{\mathbf{k},+}={r}_{\mathbf{k}}^{-1}{\theta }_{\mathbf{k},+},\quad {\tilde{\theta }}_{\mathbf{k},-}={r}_{\mathbf{k}}{\theta }_{\mathbf{k},-},\end{equation} \tag{ 21a }$$

$$\begin{equation}{\tilde{n}}_{\mathbf{k},+}={r}_{\mathbf{k}}{n}_{\mathbf{k},+},\quad {\tilde{n}}_{\mathbf{k},-}={r}_{\mathbf{k}}^{-1}{n}_{\mathbf{k},-},\end{equation} \tag{ 21b }$$

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

$$\begin{equation}\begin{aligned}\hfill {H}^{\mathrm{o}\mathrm{s}\mathrm{c}}& =\frac{n}{2}\sum\limits _{\mathbf{k}\ne \mathbf{0}}\left(\begin{matrix}\hfill {\tilde{\theta }}_{-\mathbf{k},+}\hfill & \hfill {\tilde{\theta }}_{-\mathbf{k},-}\hfill \end{matrix}\right)M(\mathbf{k})\left(\begin{matrix}\hfill {\tilde{\theta }}_{\mathbf{k},+}\hfill \\ \hfill {\tilde{\theta }}_{\mathbf{k},-}\hfill \end{matrix}\right)\hfill \\ \hfill & \quad +\frac{1}{2n}\sum\limits _{\mathbf{k}\ne \mathbf{0}}\sum\limits _{\nu =\pm }{e}_{\mathbf{k}}{\tilde{n}}_{-\mathbf{k},\nu }{\tilde{n}}_{\mathbf{k},\nu },\hfill \end{aligned}\end{equation} \tag{ 22 }$$

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

$$\begin{equation}\begin{aligned}\hfill M(\mathbf{k})& {:=}\left(\begin{matrix}\hfill {r}_{\mathbf{k}}^{2}{{\Gamma}}_{\pm }(\mathbf{k})\hfill & \hfill \pm \mathrm{i}{\Gamma}(\mathbf{k})\hfill \\ \hfill \mp \mathrm{i}{\Gamma}(\mathbf{k})\hfill & \hfill {r}_{\mathbf{k}}^{-2}{{\Gamma}}_{\mp }(\mathbf{k})\hfill \end{matrix}\right)\hfill \\ \hfill & ={{\Gamma}}_{0}(\mathbf{k})I\mp {\Gamma}(\mathbf{k}){\sigma }_{y}+{{\Gamma}}_{z}(\mathbf{k}){\sigma }_{z},\hfill \end{aligned}\end{equation} \tag{ 23 }$$

where I is the identity matrix, (σ_x , σ_y , σ_z ) are the Pauli matrices, and

$$\begin{equation}{{\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],\end{equation} \tag{ 24a }$$

$$\begin{equation}{{\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].\end{equation} \tag{ 24b }$$

We then perform the second canonical transformation using the unitary matrix U(k) as

$$\begin{equation}\left(\begin{matrix}\hfill {\tilde{\theta }}_{\mathbf{k},+}\hfill \\ \hfill {\tilde{\theta }}_{\mathbf{k},-}\hfill \end{matrix}\right)=U(\mathbf{k})\left(\begin{matrix}\hfill {\bar{\theta }}_{\mathbf{k},1}\hfill \\ \hfill {\bar{\theta }}_{\mathbf{k},2}\hfill \end{matrix}\right),\quad \left(\begin{matrix}\hfill {\tilde{n}}_{\mathbf{k},+}\hfill \\ \hfill {\tilde{n}}_{\mathbf{k},-}\hfill \end{matrix}\right)=U(\mathbf{k})\left(\begin{matrix}\hfill {\bar{n}}_{\mathbf{k},1}\hfill \\ \hfill {\bar{n}}_{\mathbf{k},2}\hfill \end{matrix}\right).\end{equation} \tag{ 25 }$$

We note that the second term of equation (22) is invariant under this transformation if U(−k)^t U(k) = I for all k ≠ 0. It is therefore useful to choose U(k) in such a way as to diagonalize the Hermitian matrix M(k) as

$$\begin{equation}{U}^{-1}(\mathbf{k})M(\mathbf{k})U(\mathbf{k})=\left(\begin{matrix}\hfill {m}_{\mathbf{k},1}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {m}_{\mathbf{k},2}\hfill \end{matrix}\right),\end{equation} \tag{ 26a }$$

$$\begin{equation}U(\mathbf{k})={\mathrm{e}}^{\mathrm{i}{\chi }_{\mathbf{k}}{\sigma }_{x}/2}=\left(\begin{matrix}\hfill \mathrm{cos}({\chi }_{\mathbf{k}}/2)\hfill & \hfill \mathrm{i}\,\mathrm{sin}({\chi }_{\mathbf{k}}/2)\hfill \\ \hfill \mathrm{i}\,\mathrm{sin}({\chi }_{\mathbf{k}}/2)\hfill & \hfill \mathrm{cos}({\chi }_{\mathbf{k}}/2)\hfill \end{matrix}\right),\end{equation} \tag{ 26b }$$

where ⁸

$$\begin{equation}{m}_{\mathbf{k},j}{:=}{{\Gamma}}_{0}(\mathbf{k})+{(-1)}^{j-1}{\Lambda}(\mathbf{k})\quad (j=1,2),\end{equation} \tag{ 27a }$$

$$\begin{equation}{\Lambda}(\mathbf{k})(\mathrm{cos}\,{\chi }_{\mathbf{k}},\mathrm{sin}\,{\chi }_{\mathbf{k}}){:=}\left({{\Gamma}}_{z}(\mathbf{k}),{\Gamma}(\mathbf{k})\right).\end{equation} \tag{ 27b }$$

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

$$\begin{equation}{H}^{\mathrm{o}\mathrm{s}\mathrm{c}}=\frac{1}{2}\sum\limits _{\mathbf{k}\ne \mathbf{0}}\sum\limits _{j=1,2}\left(n{m}_{\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).\end{equation} \tag{ 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 }_{j{j}^{\prime }}{\delta }_{\mathbf{k}{\mathbf{k}}^{\prime }}\enspace (\mathbf{k},{\mathbf{k}}^{\prime }\ne \mathbf{0};\enspace j,{j}^{\prime }=1,2)$ as the transformations in equations (21) and (25) leave the commutation relations unchanged. Finally, introducing the (bogolon) annihilation and creation operators

$$\begin{equation}{\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),\end{equation} \tag{ 29a }$$

$$\begin{equation}{\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)\quad (\mathbf{k}\ne \mathbf{0};\enspace j=1,2),\end{equation} \tag{ 29b }$$

with ${\zeta }_{\mathbf{k},j}{:=}{({m}_{\mathbf{k},j}/{e}_{\mathbf{k}})}^{1/4}$, we can diagonalize the Hamiltonian (28) as

$$\begin{equation}{H}^{\mathrm{o}\mathrm{s}\mathrm{c}}=\sum\limits _{\mathbf{k}\ne \mathbf{0}}\sum\limits _{j=1,2}{E}_{j}(\mathbf{k})\left({\gamma }_{\mathbf{k},j}^{{\dagger}}{\gamma }_{\mathbf{k},j}+\frac{1}{2}\right),\end{equation} \tag{ 30a }$$

$$\begin{equation}{E}_{j}(\mathbf{k}){:=}\sqrt{{m}_{\mathbf{k},j}{e}_{\mathbf{k}}}.\end{equation} \tag{ 30b }$$

The ground state $\left\vert {0}^{\mathrm{o}\mathrm{s}\mathrm{c}}\right\rangle$ of this Hamiltonian is specified by the condition that ${\gamma }_{\mathbf{k},\pm }\left\vert {0}^{\mathrm{o}\mathrm{s}\mathrm{c}}\right\rangle =0$ for all k ≠ 0.

We now discuss the single-particle spectrum E_j (k) (j = 1, 2) in the long-wavelength limit kℓ ≪ 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 equation (19), we can approximate m_k,j in equation (27) as

$$\begin{equation}{m}_{\mathbf{k},1}\approx {r}_{\mathbf{0}}^{\mp 2}{{\Gamma}}_{-}(\mathbf{k}),\quad {m}_{\mathbf{k},2}\approx {r}_{\mathbf{0}}^{\pm 2}\left[{{\Gamma}}_{+}(\mathbf{k})-\frac{{\Gamma}{(\mathbf{k})}^{2}}{{{\Gamma}}_{-}(\mathbf{k})}\right].\end{equation} \tag{ 31 }$$

By parametrizing the wave vector as (k_x , k_y ) = k(cos φ, sin φ), m_k,j can also be expressed as

$$\begin{equation}{m}_{\mathbf{k},1}\approx {r}_{\mathbf{0}}^{\mp 2}gnD(\varphi ){(k\ell )}^{2},\quad {m}_{\mathbf{k},2}\approx {r}_{\mathbf{0}}^{\pm 2}gnC(\varphi ){(k\ell )}^{4},\end{equation} \tag{ 32 }$$

where

$$\begin{align}\hfill C(\varphi )& {:=}{C}_{1}\,{\mathrm{sin}}^{2}(2\varphi )+{C}_{2}\,{\mathrm{cos}}^{2}(2\varphi )\hfill \\ \hfill & \enspace \enspace \enspace \enspace \enspace -{C}_{3}\,\mathrm{sin}(2\varphi )\mathrm{cos}(2\varphi )-{C}_{4}\,{\mathrm{sin}}^{2}(3\varphi ),\hfill \end{align} \tag{ 33a }$$

$$\begin{align}\hfill D(\varphi )& {:=}{D}_{1}\,{\mathrm{sin}}^{2}(\varphi )+{D}_{2}\,{\mathrm{cos}}^{2}(\varphi )-{D}_{3}\,\mathrm{sin}(\varphi )\mathrm{cos}(\varphi ),\hfill \end{align} \tag{ 33b }$$

with ${C}_{4}{:=}{F}_{1}^{2}/4{D}_{1}$. We then obtain the low-energy spectrum as

$$\begin{equation}{E}_{j}(\mathbf{k})\approx \sqrt{{m}_{\mathbf{k},j}{e}_{\mathbf{0}}}\approx \sqrt{2}gn{(k\ell )}^{j}{f}_{j}(\varphi )\quad (j=1,2),\end{equation} \tag{ 34 }$$

where the dependence on the angle φ is expressed by the dimensionless functions

$$\begin{equation}{f}_{1}(\varphi )=\sqrt{{\bar{g}}_{\mp }D(\varphi )/g},\quad {f}_{2}(\varphi )=\sqrt{{\bar{g}}_{\pm }C(\varphi )/g}.\end{equation} \tag{ 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:

$$\begin{equation}{E}_{1}^{\text{P}}(\mathbf{k})/\sqrt{{\bar{g}}_{-}}={E}_{1}^{\text{AP}}(\mathbf{k})/\sqrt{{\bar{g}}_{+}},\end{equation} \tag{ 36a }$$

$$\begin{equation}{E}_{2}^{\text{P}}(\mathbf{k})/\sqrt{{\bar{g}}_{+}}={E}_{2}^{\text{AP}}(\mathbf{k})/\sqrt{{\bar{g}}_{-}}.\end{equation} \tag{ 36b }$$

Using the dimensionless functions f_j (φ) (j = 1, 2) for the two types of fields, these relations can also be written as

$$\begin{equation}{f}_{1}^{\text{P}}(\varphi )\sqrt{\frac{g}{{\bar{g}}_{-}}}={f}_{1}^{\text{AP}}(\varphi )\sqrt{\frac{g}{{\bar{g}}_{+}}}=\sqrt{D(\varphi )},\end{equation} \tag{ 37a }$$

$$\begin{equation}{f}_{2}^{\text{P}}(\varphi )\sqrt{\frac{g}{{\bar{g}}_{+}}}={f}_{2}^{\text{AP}}(\varphi )\sqrt{\frac{g}{{\bar{g}}_{-}}}=\sqrt{C(\varphi )}.\end{equation} \tag{ 37b }$$

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

We now calculate the reduced density matrix (RDM) ρ_↑ for the spin-↑ component, which is defined by starting from the ground state $\left\vert {0}^{\mathrm{z}\mathrm{e}\mathrm{r}\mathrm{o}}\right\rangle \otimes \left\vert {0}^{\mathrm{o}\mathrm{s}\mathrm{c}}\right\rangle$ of the total system and trancing out the degrees of freedom in the spin-↓ 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}}^{\text{zero}}\otimes {\rho }_{{\uparrow}}^{\text{osc}}$. As the zero-mode ground state $\left\vert {N}_{{\uparrow}}=N/2\right\rangle \left\vert {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{z}\mathrm{e}\mathrm{r}\mathrm{o}}=\left\vert {N}_{{\uparrow}}=N/2\right\rangle \left\langle {N}_{{\uparrow}}=N/2\right\vert$. Below we consider the oscillator-mode part ${\rho }_{{\uparrow}}^{\text{osc}}$.

For ${\rho }_{{\uparrow}}^{\mathrm{o}\mathrm{s}\mathrm{c}}$, we introduce the following Gaussian ansatz [61–63, 67–69]:

$$\begin{equation}{\rho }_{{\uparrow}}^{\mathrm{o}\mathrm{s}\mathrm{c}}=\frac{1}{{Z}_{\text{e}}^{\mathrm{o}\mathrm{s}\mathrm{c}}}\,{\mathrm{e}}^{-{H}_{\text{e}}^{\mathrm{o}\mathrm{s}\mathrm{c}}},\quad {Z}_{\text{e}}^{\mathrm{o}\mathrm{s}\mathrm{c}}=\mathrm{Tr}\,{\mathrm{e}}^{-{H}_{\text{e}}^{\mathrm{o}\mathrm{s}\mathrm{c}}},\end{equation} \tag{ 38a }$$

$$\begin{equation}{H}_{\text{e}}^{\text{osc}}=\frac{1}{2}\sum\limits _{\mathbf{k}\ne \mathbf{0}}\left(n{F}_{\mathbf{k}}{\theta }_{-\mathbf{k},{\uparrow}}{\theta }_{\mathbf{k},{\uparrow}}+\frac{{G}_{\mathbf{k}}}{n}{n}_{-\mathbf{k},{\uparrow}}{n}_{\mathbf{k},{\uparrow}}\right),\end{equation} \tag{ 38b }$$

where F_k and G_k are positive dimensionless coefficients to be determined later and we assume F_k = F_−k and G_k = G_−k for convenience. By introducing annihilation and creation operators as

$$\begin{equation}{\eta }_{\mathbf{k}}=\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],\end{equation} \tag{ 39a }$$

$$\begin{align}\hfill {\eta }_{\mathbf{k}}^{{\dagger}}=\;& \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]\hfill \\ \hfill & (\mathbf{k}\ne \mathbf{0}),\hfill \end{align} \tag{ 39b }$$

the entanglement Hamiltonian ${H}_{\mathrm{e}}^{\mathrm{o}\mathrm{s}\mathrm{c}}$ in equation (38) is diagonalized as

$$\begin{equation}{H}_{\mathrm{e}}^{\mathrm{o}\mathrm{s}\mathrm{c}}=\sum\limits _{\mathbf{k}\ne \mathbf{0}}{\xi }_{\mathbf{k}}\left({\eta }_{\mathbf{k}}^{{\dagger}}{\eta }_{\mathbf{k}}+\frac{1}{2}\right),\end{equation} \tag{ 40 }$$

where ${\xi }_{\mathbf{k}}{:=}\sqrt{{F}_{\mathbf{k}}{G}_{\mathbf{k}}}$ is the single-particle ES.

The single-particle ES ξ_k and the coefficients F_k and G_k can be determined in the following way [61]. Using the relations in equation (39) and the Bose distribution function

$$\begin{equation}\mathrm{Tr}\left({\eta }_{\mathbf{k}}^{{\dagger}}{\eta }_{\mathbf{k}}{\rho }_{{\uparrow}}^{\mathrm{o}\mathrm{s}\mathrm{c}}\right)=\frac{1}{{\mathrm{e}}^{{\xi }_{\mathbf{k}}}-1}\equiv {f}_{\mathrm{B}}({\xi }_{\mathbf{k}}),\end{equation} \tag{ 41 }$$

we obtain the phase and density correlators as

$$\begin{equation}\mathrm{Tr}\left({\theta }_{-\mathbf{k},{\uparrow}}{\theta }_{\mathbf{k},{\uparrow}}{\rho }_{{\uparrow}}^{\mathrm{o}\mathrm{s}\mathrm{c}}\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],\end{equation} \tag{ 42a }$$

$$\begin{equation}\mathrm{Tr}\left({n}_{-\mathbf{k},{\uparrow}}{n}_{\mathbf{k},{\uparrow}}{\rho }_{{\uparrow}}^{\mathrm{o}\mathrm{s}\mathrm{c}}\right)=n{\left(\frac{{F}_{\mathbf{k}}}{{G}_{\mathbf{k}}}\right)}^{1/2}\!\left[{f}_{\mathrm{B}}({\xi }_{\mathbf{k}})+\frac{1}{2}\right].\end{equation} \tag{ 42b }$$

We can then determine f_B(ξ_k ) and F_k /G_k by requiring these correlators to be equal to the same correlators calculated for the oscillator ground state $\left\vert {0}^{\text{osc}}\right\rangle$ of the total system. Details of this calculation are described in appendix A. In the long-wavelength limit kℓ ≪ 1, we obtain

$$\begin{equation}{\xi }_{\mathbf{k}}\approx c(\varphi )\sqrt{k\ell },\quad {F}_{\mathbf{k}}\approx F(\varphi ){(k\ell )}^{2},\quad {G}_{\mathbf{k}}\approx \frac{G(\varphi )}{k\ell },\end{equation} \tag{ 43 }$$

where the dependences on the angle φ are expressed by the dimensionless functions

$$\begin{equation}c(\varphi )=4{\left[\frac{{\bar{g}}_{\mp }C(\varphi )}{{\bar{g}}_{\pm }D(\varphi )}\right]}^{1/4},\end{equation} \tag{ 44a }$$

$$\begin{equation}F(\varphi )=4\sqrt{\frac{2gC(\varphi )}{{\bar{g}}_{\pm }}},\quad G(\varphi )=2\sqrt{\frac{2{\bar{g}}_{\mp }}{gD(\varphi )}}.\end{equation} \tag{ 44b }$$

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 equation (36)), the single-particle ES ξ_k ^P for parallel (P) fields and that ξ_k ^AP for antiparallel (AP) fields are related by suitable rescaling as

$$\begin{equation}{\left(\frac{{\bar{g}}_{+}}{{\bar{g}}_{-}}\right)}^{1/4}{\xi }_{\mathbf{k}}^{\text{P}}={\left(\frac{{\bar{g}}_{-}}{{\bar{g}}_{+}}\right)}^{1/4}{\xi }_{\mathbf{k}}^{\text{AP}}.\end{equation} \tag{ 45 }$$

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

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

The entanglement Hamiltonian is given in the long-wavelength limit by ${H}_{\mathrm{e}}^{\mathrm{o}\mathrm{s}\mathrm{c}}$ in equation (38) with F_k and G_k in equation (43). Using the fields θ_↑(r) and n_↑(r) in real space, it can be expressed as

$$\begin{equation}\begin{aligned}\hfill {H}_{\mathrm{e}}& =\int {\mathrm{d}}^{2}\mathbf{r}\int {\mathrm{d}}^{2}{\mathbf{r}}^{\prime }\left[\frac{n{\ell }^{2}}{2}{U}_{F}(\mathbf{r}-{\mathbf{r}}^{\prime })\nabla {\theta }_{{\uparrow}}(\mathbf{r})\cdot \nabla {\theta }_{{\uparrow}}({\mathbf{r}}^{\prime })\right.\hfill \\ \hfill & \quad \left.+\frac{1}{2n}{U}_{G}(\mathbf{r}-{\mathbf{r}}^{\prime }){n}_{{\uparrow}}(\mathbf{r}){n}_{{\uparrow}}({\mathbf{r}}^{\prime })\right],\hfill \end{aligned}\end{equation} \tag{ 47 }$$

where we introduce the interaction potentials

$$\begin{equation}{U}_{F}(\mathbf{r}-{\mathbf{r}}^{\prime })=\frac{1}{A}\sum\limits _{\mathbf{k}}F(\varphi ){\mathrm{e}}^{\mathrm{i}\mathbf{k}\cdot (\mathbf{r}-{\mathbf{r}}^{\prime })},\end{equation} \tag{ 48a }$$

$$\begin{equation}{U}_{G}(\mathbf{r}-{\mathbf{r}}^{\prime })=\underset{\alpha \to {0}^{+}}{\mathrm{lim}}\,\frac{1}{A}\sum\limits _{\mathbf{k}\ne \mathbf{0}}\frac{G(\varphi )}{k\ell }\,{\mathrm{e}}^{-\alpha k+\mathrm{i}\mathbf{k}\cdot (\mathbf{r}-{\mathbf{r}}^{\prime })}.\end{equation} \tag{ 48b }$$

Here, we use the convergence factor e^−αk to regularize the infinite sum for U_G (r − r'). For simplicity, we consider the case of overlapping triangular lattices, in which C(φ) and D(φ) (and thus F(φ) and G(φ) as well) are constant; see (a) in equation (15). In this case, the potentials in equation (48) are calculated as

$$\begin{equation}{U}_{F}(\mathbf{r}-{\mathbf{r}}^{\prime })=F\delta (\mathbf{r}-{\mathbf{r}}^{\prime }),\end{equation} \tag{ 49a }$$

$$\begin{equation}{U}_{G}(\mathbf{r}-{\mathbf{r}}^{\prime })=\frac{G}{2\pi \ell \vert \mathbf{r}-{\mathbf{r}}^{\prime }\vert }.\end{equation} \tag{ 49b }$$

For the calculation of U_G (r − r'), we refer the reader to appendix A of reference [61]. We note that θ_↑(r) is the regular part of the superfluid phase of the spin-↑ component, and that its gradient is related to the regular part of the superfluid velocity, ${\mathbf{v}}_{\mathrm{s},{\uparrow}}(\mathbf{r})=-\frac{\hslash }{M}\nabla {\theta }_{{\uparrow}}(\mathbf{r})$. Therefore, the entanglement Hamiltonian (47) has a short-range interaction in terms of the superfluid velocity v_s,↑(r) and a long-range one in terms of the density n_↑(r). If the density-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 equation (43) is closely related with the presence of a long-range interaction in H_e.

Using the single-particle ES ξ_k in equation (43), we can calculate the intercomponent EE S_e. For simplicity, we assume that ξ_k is isotropic, i.e., c(φ) 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 reference [61] for the derivation)

$$\begin{equation}\begin{aligned}\hfill {S}_{\mathrm{e}}& =\frac{\sigma A}{{c}^{4}{\ell }^{2}}-\frac{1}{2}\,\mathrm{ln}\,\frac{\sqrt{A}}{2\pi {c}^{2}\ell }+O(1)\hfill \\ \hfill & =\frac{2\pi \sigma {N}_{\mathrm{v}}}{{c}^{4}}-\frac{1}{4}\,\mathrm{ln}\,\frac{{N}_{\mathrm{v}}}{2\pi {c}^{4}}+O(1).\hfill \end{aligned}\end{equation} \tag{ 50 }$$

Here, the leading contribution is given by the first term, which is proportional to the area A with a non-universal coefficient σ 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 determined 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

$$\begin{equation}\underset{{N}_{\mathrm{v}}\to \infty }{\mathrm{lim}}\,\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}.\end{equation} \tag{ 51 }$$

Because of the factor ${\bar{g}}_{\pm }/{\bar{g}}_{\mp }$ in equation (51), when the intercomponent interaction g_↑↓ 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_e is simple yet approximate as it is based on the single-particle ES ξ_k for long wavelengths. In section. 4.3, we will present numerical results on S_e based on the Bogoliubov theory, in which ξ_k over the full Brillouin zone is taken into account, and confirm the consistency with the field-theoretical predictions.

Here we calculate some intracomponent correlation functions, and discuss their connections with the (long-wavelength) entanglement Hamiltonian H_e obtained in the preceding section. Let $\langle \mathcal{O}\rangle$ denote the expectation value of an operator $\mathcal{O}$ with respect to the ground state $\left\vert {0}^{\mathrm{z}\mathrm{e}\mathrm{r}\mathrm{o}}\right\rangle \otimes \left\vert {0}^{\mathrm{o}\mathrm{s}\mathrm{c}}\right\rangle$ of the total system. If $\mathcal{O}$ acts only on the spin-↑ component, $\langle \mathcal{O}\rangle$ should be equal to $\mathrm{Tr}\left(\mathcal{O}{\mathrm{e}}^{-{H}_{\mathrm{e}}}\right)/\mathrm{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_e manifest themselves in the correlation properties of the system.

Owing to the gapless ES ξ_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, equation (42) gives

$$\begin{equation}\langle {\theta }_{-\mathbf{k},{\uparrow}}{\theta }_{\mathbf{k},{\uparrow}}\rangle \approx \frac{1}{n{F}_{\mathbf{k}}}\approx \frac{1}{nF(\varphi ){k}^{2}{\ell }^{2}},\end{equation} \tag{ 52a }$$

$$\begin{equation}\langle {n}_{-\mathbf{k},{\uparrow}}{n}_{\mathbf{k},{\uparrow}}\rangle \approx \frac{n}{{G}_{\mathbf{k}}}\approx \frac{nk\ell }{G(\varphi )}\quad (\mathbf{k}\ne \mathbf{0}),\end{equation} \tag{ 52b }$$

where we use equation (43). Therefore, the phase and density fluctuations are directly related to the coefficients F_k and G_k , respectively, in the entanglement Hamiltonian (38). Equation (52) indicates that in the long-wavelength limit k → 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–72]. To analyze its behavior in our course-grained description, we introduce the modified bosonic field ${\tilde{\psi }}_{{\uparrow}}={\text{e}}^{-\text{i}{\theta }_{{\uparrow}}}\sqrt{{n}_{{\uparrow}}}$. As θ_↑(r) is the regular part of the superfluid phase and n_↑(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

$$\begin{equation}\langle {\tilde{\psi }}_{{\uparrow}}{(\mathbf{r})}^{{\dagger}}{\tilde{\psi }}_{{\uparrow}}(\mathbf{0})\rangle =\left\langle \sqrt{{n}_{{\uparrow}}(\mathbf{r})}{\mathrm{e}}^{\text{i}({\theta }_{{\uparrow}}(\mathbf{r})-{\theta }_{{\uparrow}}(\mathbf{0}))}\sqrt{{n}_{{\uparrow}}(\mathbf{0})}\right\rangle ,\end{equation} \tag{ 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 equation (52). Here we again consider the case of overlapping triangular lattices where F(φ) is constant. By using equation (52), the phase correlation function in real space is obtained as (see appendix B for the derivation)

$$\begin{equation}\begin{aligned}\hfill \left\langle {\left[{\theta }_{{\uparrow}}(\mathbf{r})-{\theta }_{{\uparrow}}(\mathbf{0})\right]}^{2}\right\rangle & =\frac{2}{A}\sum\limits _{\mathbf{k}\ne \mathbf{0}}{\mathrm{e}}^{-\alpha k}\left[1-\mathrm{cos}\left(\mathbf{k}\cdot \mathbf{r}\right)\right]\langle {\theta }_{-\mathbf{k},{\uparrow}}{\theta }_{\mathbf{k},{\uparrow}}\rangle \hfill \\ \hfill & \approx \frac{1}{\pi nF{\ell }^{2}}\,\mathrm{ln}\,\frac{r}{2\alpha }\quad (r\gg \alpha ),\hfill \end{aligned}\end{equation} \tag{ 54 }$$

where we introduce the same convergence factor e^−αk as before to regularize the infinite sum. The modified one-particle density matrix is then obtained as

$$\begin{equation}\begin{aligned}\hfill \langle {\tilde{\psi }}_{{\uparrow}}{(\mathbf{r})}^{{\dagger}}{\tilde{\psi }}_{{\uparrow}}(\mathbf{0})\rangle & \approx n\,\mathrm{exp}\left\{-\frac{1}{2}\left\langle {\left[{\theta }_{{\uparrow}}(\mathbf{r})-{\theta }_{{\uparrow}}(\mathbf{0})\right]}^{2}\right\rangle \right\}\hfill \\ \hfill & \approx n{\left(\frac{r}{2\alpha }\right)}^{-\frac{1}{2\pi nF{\ell }^{2}}}.\hfill \end{aligned}\end{equation} \tag{ 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₀ of the condensate can be evaluated from the one-particle density matrix (55) at the large separation $r=\sqrt{A}/2$. The density n' of the depletion is then given by n' = n − n₀. Assuming n' ≪ n as required for the Bogoliubov theory and the present effective theory, the fraction of depletion is estimated as

$$\begin{equation}\begin{aligned}\hfill \frac{{n}^{\prime }}{n}& \approx {\left.\frac{1}{2}\left\langle {\left[{\theta }_{{\uparrow}}(\mathbf{r})-{\theta }_{{\uparrow}}(\mathbf{0})\right]}^{2}\right\rangle \right\vert }_{r=\sqrt{A}/2}\hfill \\ \hfill & \approx \frac{1}{2\pi nF{\ell }^{2}}\mathrm{ln}\frac{\sqrt{A}}{4\alpha }=\frac{1}{2\nu }{\left(\frac{{\bar{g}}_{\pm }}{2gC}\right)}^{1/2}\,\mathrm{ln}\,\frac{\sqrt{A}}{4\alpha }.\hfill \end{aligned}\end{equation} \tag{ 56 }$$

where we use equation (44) and ν = 2nA/N_v = 4πnℓ². For fixed ν, the fraction of quantum depletion of the condensate increases logarithmically as a function of the number of vortices N_v = A/(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 equation (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 equations (55) and (56) have also been discussed for a vortex lattice in a scalar BEC [18, 20–22].

We employ the LLL magnetic Bloch states {Ψ_{k α} (r)} [22, 73–75] as a convenient single-particle basis for describing vortex lattices. The expressions of these states are shown in reference [46], and we summarize their main features in the following. Let a₁ and a₂ be the primitive vectors of a vortex lattice, and let K_α be the pseudomomentum operator for a spin-α atom in a synthetic magnetic field B_α  (α =  ↑, ↓). As expected for a 'Bloch state', Ψ_{k α} (r) is an eigenstate of the magnetic translation ${\mathrm{e}}^{-\mathrm{i}{\mathbf{K}}_{\alpha }\cdot {\mathbf{a}}_{j}/\hslash }$ with an eigenvalue ${\mathrm{e}}^{-\text{i}\mathbf{k}\cdot {\mathbf{a}}_{j}}\enspace (j=1,2)$. By taking N_v discrete wave vectors k consistent with the boundary conditions of the system, {Ψ_{k α} (r)} form a complete orthogonal basis of the LLL manifold ⁹ . Notably, Ψ_{k α} (r) has a periodic pattern of zeros at [74]

$$\begin{equation}\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},\quad {n}_{1},{n}_{2}\in \mathbb{Z}.\end{equation} \tag{ 57 }$$

Therefore, Ψ_{k α} (r) represents a vortex lattice with primitive vectors a₁ and a₂ for any k, and the locations of vortices (zeros) can be shifted by varying k. Vortex lattices of binary BECs in figure 1 are obtained when spin-α bosons condense into ${{\Psi}}_{{\mathbf{q}}_{\alpha },\alpha }(\mathbf{r})$, where the wave vectors q_↑ and q_↓ are chosen in a way consistent with the displacement u₁ a₁ + u₂ a₂ 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_int. Using the LLL magnetic Bloch states {Ψ_{k α} (r)} as the basis, it is represented as

$$\begin{equation}{H}_{\text{int}}=\frac{1}{2}\sum\limits _{\alpha ,\beta }\;\sum\limits _{{\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}_{{\mathbf{k}}_{1}\alpha }^{{\dagger}}{b}_{{\mathbf{k}}_{2}\beta }^{{\dagger}}{b}_{{\mathbf{k}}_{3}\beta }{b}_{{\mathbf{k}}_{4}\alpha },\end{equation} \tag{ 58 }$$

where b_{k α} is a bosonic annihilation operator for the state Ψ_{k α} (r) and

$$\begin{equation}\begin{aligned}\hfill & {V}_{\alpha \beta }({\mathbf{k}}_{1},{\mathbf{k}}_{2},{\mathbf{k}}_{3},{\mathbf{k}}_{4})\hfill \\ \hfill & \hspace{25.0pt}={g}_{\alpha \beta }\int {\mathrm{d}}^{2}\mathbf{r}\,{{\Psi}}_{{\mathbf{k}}_{1}\alpha }^{\ast }(\mathbf{r}){{\Psi}}_{{\mathbf{k}}_{2}\beta }^{\ast }(\mathbf{r}){{\Psi}}_{{\mathbf{k}}_{3}\beta }(\mathbf{r}){{\Psi}}_{{\mathbf{k}}_{4}\alpha }(\mathbf{r}).\hfill \end{aligned}\end{equation} \tag{ 59 }$$

The expression of the interaction matrix element V_αβ (k₁, k₂, k₃, k₄) that is convenient for numerical calculations is given in reference [46].

We now apply the Bogoliubov approximation [18, 21, 22, 70], assuming that the condensation occurs at the wave vector q_α in the spin-α component. To this end, it is useful to introduce

$$\begin{equation}{\tilde{b}}_{\mathbf{k}\alpha }{:=}{b}_{{\mathbf{q}}_{\alpha }+\mathbf{k},\alpha },\end{equation} \tag{ 60a }$$

$$\begin{equation}{\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}).\end{equation} \tag{ 60b }$$

By substituting

$$\begin{equation}{\tilde{b}}_{\mathbf{0}\alpha }\simeq {\tilde{b}}_{\mathbf{0}\alpha }^{{\dagger}}\simeq \sqrt{{N}_{\alpha }-\sum\limits _{\mathbf{k}\ne \mathbf{0}}{\tilde{b}}_{\mathbf{k}\alpha }^{{\dagger}}{\tilde{b}}_{\mathbf{k}\alpha }}\end{equation} \tag{ 61 }$$

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

$$\begin{equation}\begin{aligned}\hfill {H}_{\text{int}}=\;& \frac{1}{2}\sum\limits _{\alpha ,\beta }{N}_{\alpha }{N}_{\beta }{\tilde{V}}_{\alpha \beta }(\mathbf{0},\mathbf{0},\mathbf{0},\mathbf{0})-\frac{1}{2}\sum\limits _{\mathbf{k}\ne \mathbf{0}}J(\mathbf{k})\hfill \\ \hfill & +\frac{1}{2}\sum\limits _{\mathbf{k}\ne \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})\left(\begin{matrix}\hfill {\tilde{b}}_{\mathbf{k}{\uparrow}}\hfill \\ \hfill {\tilde{b}}_{\mathbf{k}{\downarrow}}\hfill \\ \hfill {\tilde{b}}_{-\mathbf{k},{\uparrow}}^{{\dagger}}\hfill \\ \hfill {\tilde{b}}_{-\mathbf{k},{\downarrow}}^{{\dagger}}\hfill \end{matrix}\right),\hfill \end{aligned}\end{equation} \tag{ 62 }$$

where

$$\begin{equation}\begin{aligned}\hfill J(\mathbf{k}){:=}& \sum\limits _{\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]\hfill \\ \hfill & +\sum\limits _{\alpha }{N}_{\alpha }{\tilde{V}}_{\alpha \alpha }(\mathbf{k},\mathbf{0},\mathbf{k},\mathbf{0}),\hfill \end{aligned}\end{equation} \tag{ 63 }$$

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

To diagonalize equation (62), we perform the Bogoliubov transformation

$$\begin{equation}\left(\begin{matrix}\hfill {\tilde{b}}_{\mathbf{k}{\uparrow}}\hfill \\ \hfill {\tilde{b}}_{\mathbf{k}{\downarrow}}\hfill \\ \hfill {\tilde{b}}_{-\mathbf{k},{\uparrow}}^{{\dagger}}\hfill \\ \hfill {\tilde{b}}_{-\mathbf{k},{\downarrow}}^{{\dagger}}\hfill \end{matrix}\right)=W(\mathbf{k})\left(\begin{matrix}\hfill {\gamma }_{\mathbf{k},1}\hfill \\ \hfill {\gamma }_{\mathbf{k},2}\hfill \\ \hfill {\gamma }_{-\mathbf{k},1}^{{\dagger}}\hfill \\ \hfill {\gamma }_{-\mathbf{k},2}^{{\dagger}}\hfill \end{matrix}\right),\end{equation} \tag{ 64a }$$

$$\begin{equation}W(\mathbf{k})=\left(\begin{matrix}\hfill \mathcal{U}(\mathbf{k})\hfill & \hfill {\mathcal{V}}^{\ast }(-\mathbf{k})\hfill \\ \hfill \mathcal{V}(\mathbf{k})\hfill & \hfill {\mathcal{U}}^{\ast }(-\mathbf{k})\hfill \end{matrix}\right).\end{equation} \tag{ 64b }$$

Here, the paraunitary matrix W(k) is chosen to satisfy

$$\begin{equation}\begin{aligned}\hfill & {\tau }_{3}\mathcal{M}(\mathbf{k})W(\mathbf{k})\hfill \\ \hfill & \hspace{25.0pt}=W(\mathbf{k})\mathrm{diag}({E}_{1}(\mathbf{k}),{E}_{2}(\mathbf{k}),-{E}_{1}(-\mathbf{k}),-{E}_{2}(-\mathbf{k})),\hfill \end{aligned}\end{equation} \tag{ 65 }$$

where τ₃ := diag(1, 1, −1, −1). Namely, W(k) and E_j (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{equation}\begin{aligned}\hfill {H}_{\text{int}}& =\frac{1}{2}\sum\limits _{\alpha ,\beta }{N}_{\alpha }{N}_{\beta }{\tilde{V}}_{\alpha \beta }(\mathbf{0},\mathbf{0},\mathbf{0},\mathbf{0})-\frac{1}{2}\sum\limits _{\mathbf{k}\ne \mathbf{0}}J(\mathbf{k})\hfill \\ \hfill & \quad +\sum\limits _{\mathbf{k}\ne \mathbf{0}}\sum\limits _{j=1,2}{E}_{j}(\mathbf{k})\left({\gamma }_{\mathbf{k}j}^{{\dagger}}{\gamma }_{\mathbf{k}j}+\frac{1}{2}\right).\hfill \end{aligned}\end{equation} \tag{ 66 }$$

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

The ground state of equation (66) is given by the bogolon vacuum $\left\vert 0\right\rangle$, which is specified by the condition that ${\gamma }_{\mathbf{k}j}\left\vert 0\right\rangle =0$ for all k ≠ 0 and j = 1, 2. The ground-state energy E_GS (scaled by the interaction energy scale gn² A) is therefore given by

$$\begin{equation}\begin{aligned}\hfill \frac{{E}_{\mathrm{G}\mathrm{S}}}{g{n}^{2}A}& =\frac{1}{2}\sum\limits _{\alpha ,\beta }\frac{A}{g}{\tilde{V}}_{\alpha \beta }(\mathbf{0},\mathbf{0},\mathbf{0},\mathbf{0})\hfill \\ \hfill & +\frac{1}{gn\nu {N}_{\mathrm{v}}}\sum\limits _{\mathbf{k}\ne \mathbf{0}}\left[\sum\limits _{j}{E}_{j}(\mathbf{k})-J(\mathbf{k})\right]\hfill \\ \hfill & =\frac{1}{2}\sum\limits _{\alpha ,\beta }\frac{A}{g}{\tilde{V}}_{\alpha \beta }(\mathbf{0},\mathbf{0},\mathbf{0},\mathbf{0})\hfill \\ \hfill & +\frac{1}{2\pi gn\nu }{\int }_{\text{BZ}}{d}^{2}\mathbf{k}{\ell }^{2}\left[\sum\limits _{j}{E}_{j}(\mathbf{k})-J(\mathbf{k})\right].\hfill \end{aligned}\end{equation} \tag{ 67 }$$

Here, in the final expression, we take the thermodynamic limit N_v → ∞ so that the sum is replaced by the integral over the Brillouin zone as $\frac{1}{{N}_{\mathrm{v}}}{\sum }_{\mathbf{k}\ne \mathbf{0}}\to \frac{1}{\vert {\mathbf{b}}_{1}\times {\mathbf{b}}_{2}\vert }\int {\mathrm{d}}^{2}\mathbf{k}=\frac{1}{2\pi }{\int }_{\text{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 on the right-hand side of equation (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 ν. In section 4.5, we numerically calculate equation (67), and discuss how the quantum correction affects the ground-state phase diagrams.

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

$$\begin{equation}\left\langle 0\right\vert {\tilde{b}}_{\mathbf{k},\alpha }^{{\dagger}}{\tilde{b}}_{\mathbf{k},\alpha }\left\vert 0\right\rangle =\sum\limits _{j}\vert {\mathcal{V}}_{\alpha ,j}(-\mathbf{k}){\vert }^{2},\end{equation} \tag{ 68a }$$

$$\begin{equation}\left\langle 0\right\vert {\tilde{b}}_{-\mathbf{k},\alpha }{\tilde{b}}_{-\mathbf{k},\alpha }^{{\dagger}}\left\vert 0\right\rangle =\sum\limits _{j}\vert {\mathcal{U}}_{\alpha ,j}(-\mathbf{k}){\vert }^{2},\end{equation} \tag{ 68b }$$

$$\begin{equation}\left\langle 0\right\vert {\tilde{b}}_{-\mathbf{k},\alpha }{\tilde{b}}_{\mathbf{k},\alpha }\left\vert 0\right\rangle =\sum\limits _{j}{\mathcal{U}}_{\alpha ,j}(-\mathbf{k}){\mathcal{V}}_{\alpha ,j}^{\ast }(-\mathbf{k}),\end{equation} \tag{ 68c }$$

$$\begin{equation}\left\langle 0\right\vert {\tilde{b}}_{\mathbf{k},\alpha }^{{\dagger}}{\tilde{b}}_{-\mathbf{k},\alpha }^{{\dagger}}\left\vert 0\right\rangle =\sum\limits _{j}{\mathcal{U}}_{\alpha ,j}^{\ast }(-\mathbf{k}){\mathcal{V}}_{\alpha ,j}(-\mathbf{k}).\end{equation} \tag{ 68d }$$

Here, we have nonzero 'anomalous' correlators in equations (68c) and (68d) as the particle numbers N_↑ and N_↓ are not conserved in the Bogoliubov Hamiltonian (62). Using equations (68a) and (68b), we can further calculate the fraction of depletion n'/n, which is equal for the two components, as

$$\begin{equation}\begin{aligned}\hfill \frac{{n}^{\prime }}{n}& =\frac{1}{{N}_{\alpha }}\sum\limits _{\mathbf{k}\ne \mathbf{0}}\left\langle 0\right\vert {\tilde{b}}_{\mathbf{k},\alpha }^{{\dagger}}{\tilde{b}}_{\mathbf{k},\alpha }\left\vert 0\right\rangle \hfill \\ \hfill & =\frac{2}{\nu {N}_{\mathrm{v}}}\sum\limits _{\mathbf{k}\ne \mathbf{0}}\sum\limits _{j=1,2}\vert {\mathcal{V}}_{\alpha ,j}(-\mathbf{k}){\vert }^{2}\quad (\alpha ={\uparrow},{\downarrow}).\hfill \end{aligned}\end{equation} \tag{ 69 }$$

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

For the RDM ρ_↑ for the spin-↑ component, we introduce the following Gaussian ansatz [62, 63, 67, 68]:

$$\begin{equation}{\rho }_{{\uparrow}}=\frac{1}{{Z}_{\mathrm{e}}}\,{\mathrm{e}}^{-{H}_{\mathrm{e}}},\quad {Z}_{\mathrm{e}}=\mathrm{Tr}\,{\mathrm{e}}^{-{H}_{\mathrm{e}}},\end{equation} \tag{ 70a }$$

$$\begin{equation}{H}_{\mathrm{e}}=\frac{1}{2}\sum\limits _{\mathbf{k}\ne \mathbf{0}}\left({\tilde{b}}_{\mathbf{k},{\uparrow}}^{{\dagger}},{\tilde{b}}_{-\mathbf{k},{\uparrow}}\right){M}_{\mathrm{e}}(\mathbf{k})\left(\begin{matrix}\hfill {\tilde{b}}_{\mathbf{k},{\uparrow}}\hfill \\ \hfill {\tilde{b}}_{-\mathbf{k},{\uparrow}}^{{\dagger}}\hfill \end{matrix}\right),\end{equation} \tag{ 70b }$$

with

$$\begin{equation}{M}_{\mathrm{e}}(\mathbf{k})=\left(\begin{matrix}\hfill {h}_{\mathbf{k}}\hfill & \hfill -{\lambda }_{\mathbf{k}}\hfill \\ \hfill -{\lambda }_{-\mathbf{k}}^{\ast }\hfill & \hfill {h}_{-\mathbf{k}}\hfill \end{matrix}\right),\quad {\lambda }_{\mathbf{k}}={\lambda }_{-\mathbf{k}}.\end{equation} \tag{ 71 }$$

By performing a Bogoliubov transformation

$$\begin{equation}\left(\begin{matrix}\hfill {\tilde{b}}_{\mathbf{k}{\uparrow}}\hfill \\ \hfill {\tilde{b}}_{-\mathbf{k},{\uparrow}}^{{\dagger}}\hfill \end{matrix}\right)={W}_{\mathrm{e}}(\mathbf{k})\left(\begin{matrix}\hfill {\eta }_{\mathbf{k}}\hfill \\ \hfill {\eta }_{-\mathbf{k}}^{{\dagger}}\hfill \end{matrix}\right),\end{equation} \tag{ 72a }$$

$$\begin{equation}{W}_{\mathrm{e}}(\mathbf{k})=\left(\begin{matrix}\hfill \mathrm{cosh}\,{\theta }_{\mathbf{k}}\hfill & \hfill {\mathrm{e}}^{-\mathrm{i}{\phi }_{\mathbf{k}}}\,\mathrm{sinh}\,{\theta }_{\mathbf{k}}\hfill \\ \hfill {\mathrm{e}}^{\mathrm{i}{\phi }_{\mathbf{k}}}\,\mathrm{sinh}\,{\theta }_{\mathbf{k}}\hfill & \hfill \mathrm{cosh}\,{\theta }_{\mathbf{k}}\hfill \end{matrix}\right)\end{equation} \tag{ 72b }$$

with

$$\begin{equation}\mathrm{cosh}\,2{\theta }_{\mathbf{k}}=\frac{{h}_{\mathbf{k}}+{h}_{-\mathbf{k}}}{2{\tilde{h}}_{\mathbf{k}}},\quad {\mathrm{e}}^{-\text{i}{\phi }_{\mathbf{k}}}\,\mathrm{sinh}\,2{\theta }_{\mathbf{k}}=\frac{{\lambda }_{\mathbf{k}}}{{\tilde{h}}_{\mathbf{k}}},\end{equation} \tag{ 73a }$$

$$\begin{equation}{\tilde{h}}_{\mathbf{k}}=\sqrt{\frac{1}{4}{\left({h}_{\mathbf{k}}+{h}_{-\mathbf{k}}\right)}^{2}-\vert {\lambda }_{\mathbf{k}}{\vert }^{2}},\end{equation} \tag{ 73b }$$

the entanglement Hamiltonian in H_e in equation (70) is diagonalized as

$$\begin{equation}\begin{aligned}\hfill {H}_{\mathrm{e}}& =\frac{1}{2}\sum\limits _{\mathbf{k}\ne \mathbf{0}}\left({\xi }_{\mathbf{k}}{\eta }_{\mathbf{k}}^{{\dagger}}{\eta }_{\mathbf{k}}+{\xi }_{-\mathbf{k}}{\eta }_{-\mathbf{k}}{\eta }_{-\mathbf{k}}^{{\dagger}}\right)\hfill \\ \hfill & =\sum\limits _{\mathbf{k}\ne \mathbf{0}}{\xi }_{\mathbf{k}}\left({\eta }_{\mathbf{k}}^{{\dagger}}{\eta }_{\mathbf{k}}+\frac{1}{2}\right),\hfill \end{aligned}\end{equation} \tag{ 74 }$$

where

$$\begin{equation}{\xi }_{\mathbf{k}}{:=}{\tilde{h}}_{\mathbf{k}}+\frac{{h}_{\mathbf{k}}-{h}_{-\mathbf{k}}}{2}\end{equation} \tag{ 75 }$$

is the single-particle ES.

Using the relation (72) and the Bose distribution functions

$$\begin{equation}\mathrm{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}}),\end{equation} \tag{ 76a }$$

$$\begin{equation}\mathrm{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}}),\end{equation} \tag{ 76b }$$

we obtain

$$\begin{equation}\mathrm{Tr}\left({\tilde{b}}_{\mathbf{k},{\uparrow}}^{{\dagger}}{\tilde{b}}_{\mathbf{k},{\uparrow}}{\rho }_{{\uparrow}}\right)={f}_{\mathrm{B}}({\xi }_{\mathbf{k}}){\mathrm{cosh}}^{2}\,{\theta }_{\mathbf{k}}-{f}_{\mathrm{B}}(-{\xi }_{-\mathbf{k}}){\mathrm{sinh}}^{2}\,{\theta }_{\mathbf{k}},\end{equation} \tag{ 77a }$$

$$\begin{equation}\mathrm{Tr}\left({\tilde{b}}_{-\mathbf{k},{\uparrow}}{\tilde{b}}_{-\mathbf{k},{\uparrow}}^{{\dagger}}{\rho }_{{\uparrow}}\right)=-{f}_{\mathrm{B}}(-{\xi }_{-\mathbf{k}}){\mathrm{cosh}}^{2}\,{\theta }_{\mathbf{k}}+{f}_{\mathrm{B}}({\xi }_{\mathbf{k}}){\mathrm{sinh}}^{2}\,{\theta }_{\mathbf{k}},\end{equation} \tag{ 77b }$$

$$\begin{equation}2\,\mathrm{Tr}\left({\tilde{b}}_{-\mathbf{k},{\uparrow}}{\tilde{b}}_{\mathbf{k},{\uparrow}}{\rho }_{{\uparrow}}\right)=\left[{f}_{\mathrm{B}}({\xi }_{\mathbf{k}})-{f}_{\mathrm{B}}(-{\xi }_{-\mathbf{k}})\right]{\mathrm{e}}^{-\text{i}{\phi }_{\mathbf{k}}}\,\mathrm{sinh}\left(2{\theta }_{\mathbf{k}}\right).\end{equation} \tag{ 77c }$$

We require these to be equal to the correlators (68) with respect to the Bogoliubov ground state. We can then express f_B(ξ_k ) in terms of the correlators (68) in the following way. First, by taking the sum and the difference of equations (77a) and (77b), we have

$$\begin{equation}\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 =\left[{f}_{\mathrm{B}}({\xi }_{\mathbf{k}})-{f}_{\mathrm{B}}(-{\xi }_{-\mathbf{k}})\right]\mathrm{cosh}\left(2{\theta }_{\mathbf{k}}\right),\end{equation} \tag{ 78a }$$

$$\begin{equation}\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 ={f}_{\mathrm{B}}({\xi }_{\mathbf{k}})+{f}_{\mathrm{B}}(-{\xi }_{-\mathbf{k}}),\end{equation} \tag{ 78b }$$

where we take the shorthand notation $\langle \cdot \rangle {:=}\left\langle 0\right\vert \cdot \left\vert 0\right\rangle$. Next, using equations (77c) and (78a), we have

$$\begin{equation}\begin{aligned}\hfill & {f}_{\mathrm{B}}({\xi }_{\mathbf{k}})-{f}_{\mathrm{B}}(-{\xi }_{-\mathbf{k}})\hfill \\ \hfill & \hspace{25.0pt}=\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{\left\vert \langle {\tilde{b}}_{-\mathbf{k},{\uparrow}}{\tilde{b}}_{\mathbf{k}{\uparrow}}\rangle \right\vert }^{2}}.\hfill \end{aligned}\end{equation} \tag{ 79 }$$

Lastly, using equations (78b) and (79), we obtain

$$\begin{align}\hfill {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)\hfill \\ \hfill & \quad +\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}-{\left\vert \langle {\tilde{b}}_{-\mathbf{k},{\uparrow}}{\tilde{b}}_{\mathbf{k}{\uparrow}}\rangle \right\vert }^{2}},\hfill \end{align} \tag{ 80 }$$

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

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

$$\begin{equation}{S}_{\mathrm{e}}=\sum\limits _{\mathbf{k}\ne \mathbf{0}}\left\{-{f}_{\mathrm{B}}({\xi }_{\mathbf{k}})\mathrm{ln}\,{f}_{\mathrm{B}}({\xi }_{\mathbf{k}})+\left[1+{f}_{\mathrm{B}}({\xi }_{\mathbf{k}})\right]\mathrm{ln}\left[1+{f}_{\mathrm{B}}({\xi }_{\mathbf{k}})\right]\right\}.\end{equation} \tag{ 81 }$$

As discussed in section 2.5 (see equation (50)), S_e is expected to show a volume-law behavior with a subleading logarithmic correction. The EE per flux quantum in the thermodynamic limit is then expressed in the integral form

$$\begin{equation}\begin{aligned}\hfill \underset{{N}_{\mathrm{v}}\to \infty }{\mathrm{lim}}\,\frac{{S}_{\mathrm{e}}}{{N}_{\mathrm{v}}}& =\frac{1}{2\pi }{\int }_{\text{BZ}}{\mathrm{d}}^{2}\mathbf{k}\,{\ell }^{2}\left\{-{f}_{\mathrm{B}}({\xi }_{\mathbf{k}})\mathrm{ln}\,{f}_{\mathrm{B}}({\xi }_{\mathbf{k}})\right.\hfill \\ \hfill & \quad \left.+\left[1+{f}_{\mathrm{B}}({\xi }_{\mathbf{k}})\right]\mathrm{ln}\left[1+{f}_{\mathrm{B}}({\xi }_{\mathbf{k}})\right]\right\}.\hfill \end{aligned}\end{equation} \tag{ 82 }$$

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

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 section 2.2. Comparing the mean-field ground-state energy (the first term of equation (67)) with the corresponding field-theoretical expression (11), we find

$$\begin{equation}\begin{aligned}\hfill {\bar{g}}_{\alpha \beta }& ={\beta }_{\alpha \beta }{g}_{\alpha \beta }=A{V}_{\alpha \beta }({\mathbf{q}}_{\alpha },{\mathbf{q}}_{\beta },{\mathbf{q}}_{\beta },{\mathbf{q}}_{\alpha })\hfill \\ \hfill & =A{g}_{\alpha \beta }\int {\mathrm{d}}^{2}\mathbf{r}\vert {{\Psi}}_{{\mathbf{q}}_{\alpha },\alpha }(\mathbf{r}){\vert }^{2}\vert {{\Psi}}_{{\mathbf{q}}_{\beta },\beta }(\mathbf{r}){\vert }^{2}.\hfill \end{aligned}\end{equation} \tag{ 83 }$$

This expression can also be obtained by substituting the condensate wave function $\sqrt{{N}_{\alpha }}{{\Psi}}_{{\mathbf{q}}_{\alpha },\alpha }(\mathbf{r})$ into ψ_α (r) in equation (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 β := β_↑↑ = β_↓↓ and β_↑↓ = β_↓↑ calculated for the mean-field vortex lattice structures. For overlapping triangular, interlaced triangular, and square lattices, β and β_↑↓ do not depend on g_↑↓/g as the lattice structures remain unchanged in the concerned regions. For rhombic and rectangular lattices, in contrast, β and β_↑↓ do depend on g_↑↓/g as the inner angle θ and the aspect ratio b/a continuously vary for the former and latter lattices, respectively.

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

As explained in section 3.2 (see equation (65)), the excitation spectrum E_j (k) (j = 1, 2) can be obtained by numerically calculating the right eigenvalues of the 4 × 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 figure 1 and for both parallel and antiparallel fields. In reference [46], ¹⁰ 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 equations (36) and (37) which are predicted by the low-energy effective field theory. In reference [46], the unrenormalized coupling constants g_αβ 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 section 4.1, we can demonstrate the rescaling relations for all the five structures (a)–(e) shown in figure 1.

Zoom In Zoom Out Reset image size

Figure 4 displays rescaled excitation spectra for the five cases (a)–(e) in figure 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 β and β_↑↓ are shown in figure 2. We can confirm that the rescaling relations in equation (36) hold at sufficiently low energies around the Γ point. Interestingly, in figures 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 figure 3 can be fit well by linear and quadratic dispersion relations

$$\begin{equation}{E}_{j}(\mathbf{k})=\sqrt{2}gn{(k\ell )}^{j}{f}_{j}(\varphi )\quad (j=1,2),\end{equation} \tag{ 84 }$$

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

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

By comparing the obtained $\sqrt{C(\varphi )}$ and $\sqrt{D(\varphi )}$ with the analytical expressions in equation (33), we can determine the dimensionless elastic constants {C_i } and {D_i }. Figure 6 shows the determined elastic constants as functions of g_↑↓/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₄ for interlaced triangular lattices (b) was overlooked in reference [45].

Zoom In Zoom Out Reset image size

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

The left panels of figure 7 display the single-particle ES ξ_k for the same cases as in figures 3–5. Around the Γ point, ξ_k can be well fitted by the square-root dispersion relation

$$\begin{equation}{\xi }_{\mathbf{k}}=c(\varphi )\sqrt{k\ell }\end{equation} \tag{ 85 }$$

for k = k(cos φ, sin φ), where c(φ) is a dimensionless function that expresses the anisotropy. The right panels of figure 7 show the function c(φ) determined from the data of ξ_k along a circular path around the Γ 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(φ)/D(φ)]^1/2 (not shown), where $\sqrt{C(\varphi )}$ and $\sqrt{D(\varphi )}$ are shown in figure 5. Thus, we can confirm the rescaling relation for the entanglement spectra in equation (46).

Zoom In Zoom Out Reset image size

In the left panels of figure 7, we also find that ξ_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. In our previous work [46] (see appendix D therein), it has been found that at the M₁ and M₂ 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-↑ and ↓ components are decoupled. Therefore, ξ_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) shows the intercomponent EE per flux quantum as a function of the ratio g_↑↓/g for parallel (blue) and antiparallel (red) fields. We find that for repulsive (attractive) g_↑↓, the EE tends to be larger for parallel (antiparallel) fields, in consistency with the field-theoretical result in section 2.5 (see equation (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.

Zoom In Zoom Out Reset image size

In figure 8(b), we examine the scaling of the intercomponent EE S_e as a function of N_v. As seen in this figure, the dominant part of the scaling is given by a volume-law behavior α₁ N_v; such a volume-law contribution is standard for an extensive cut as discussed here, and has also been found elsewhere [62, 63, 76–78]. In agreement with the field-theoretical result in equation (50), we find that the data are well fitted by the form S_e = α₁ N_v − α₂ ln N_v + α₃ and that the coefficient α₂ obtained by the fitting is close to 1/4.

Figure 9 shows numerical results on the fraction of depletion n'/n (scaled by ν⁻¹). As seen in figure 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 equation (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 figure 11. In figure 9(b), we examine the scaling of νn'/n as a function of N_v. The data are well fitted by the logarithmic form νn'/n = γ₁ ln N_v + γ₂, in agreement with the field-theoretical result in equation (56).

Zoom In Zoom Out Reset image size

Here we analyze how quantum fluctuations affect the ground-state phase diagrams for parallel and antiparallel fields. The GP mean-field analyses [36–39] have led to the five types of lattice structures that depend on g_↑↓/g as shown in figure 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 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 ν = 20. Here, the energies for rhombic and rectangular lattices are minimized with respect to the inner angle θ and the aspect ratio b/a, respectively. As seen in figure 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_↑↓/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 figure 10(a), the shift indeed occurs as seen in the enlarged plot in figure 10(b).

Zoom In Zoom Out Reset image size

Figures 11(a) and (b) show the inner angle θ 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 θ 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 θ and b/a are enhanced in these regimes, which explains the large shifts of the transition points as demonstrated in figure 10(a).

Zoom In Zoom Out Reset image size