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

We consider a system of a 2D pseudospin-$\tfrac{1}{2}$ Bose gas having two hyperfine spin states (labeled by $\alpha =\uparrow ,\,\downarrow$). The spin-α component is subject to a synthetic magnetic field B_α in the z direction. In the case of a gas rotating with an angular frequency Ω, parallel fields ${B}_{\uparrow }={B}_{\downarrow }=2M{\rm{\Omega }}/q$ are induced in the two components in the rotating frame of reference, where M and q are the mass and the fictitious charge of a neutral atom. An optical dressing technique of [49], in contrast, can be used to produce antiparallel fields ${B}_{\uparrow }=-{B}_{\downarrow }$. We focus on a central region of the system where the atomic density is sufficiently uniform and the effect of the harmonic potential can be ignored. In the second-quantized form, the Hamiltonian of the system is given by

$$\begin{eqnarray}\begin{array}{rcl}H & = & {H}_{\mathrm{kin}}+{H}_{\mathrm{int}}\\ & = & \displaystyle \sum _{\alpha =\uparrow ,\downarrow }\displaystyle \int {{\rm{d}}}^{2}{\bf{r}}\,{\hat{\psi }}_{\alpha }^{\dagger }({\bf{r}})\displaystyle \frac{{({\bf{p}}-q{{\bf{A}}}_{\alpha })}^{2}}{2M}{\hat{\psi }}_{\alpha }({\bf{r}})+\displaystyle \sum _{\alpha ,\beta }\displaystyle \frac{{g}_{\alpha \beta }}{2}\displaystyle \int {{\rm{d}}}^{2}{\bf{r}}\,{\hat{\psi }}_{\alpha }^{\dagger }({\bf{r}}){\hat{\psi }}_{\beta }^{\dagger }({\bf{r}}){\hat{\psi }}_{\beta }({\bf{r}}){\hat{\psi }}_{\alpha }({\bf{r}}),\end{array}\end{eqnarray} \tag{ 1 }$$

where ${\bf{r}}=(x,y)$ is the coordinate on the 2D plane, ${\bf{p}}=-{\rm{i}}{\hslash }({\partial }_{x},{\partial }_{y})$ is the momentum, and ${\hat{\psi }}_{\alpha }({\bf{r}})$ is the bosonic field operator for the spin-α component satisfying the commutation relations $[{\hat{\psi }}_{\alpha }({\bf{r}}),{\hat{\psi }}_{\beta }^{\dagger }({\bf{r}}^{\prime} )]={\delta }_{\alpha \beta }{\delta }^{(2)}({\bf{r}}-{\bf{r}}^{\prime} )$ and $[{\hat{\psi }}_{\alpha }({\bf{r}}),{\hat{\psi }}_{\beta }({\bf{r}}^{\prime} )]=[{\hat{\psi }}_{\alpha }^{\dagger }({\bf{r}}),{\hat{\psi }}_{\beta }^{\dagger }({\bf{r}}^{\prime} )]=0$. The gauge field for the spin-α component is given by

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

where we assume B > 0 and ${\epsilon }_{\uparrow }={\epsilon }_{\downarrow }=1$ (${\epsilon }_{\uparrow }=-{\epsilon }_{\downarrow }=1$) for parallel (antiparallel) fields. For a 2D system of area A, the number of magnetic flux quanta piercing each component (or the number of vortices) is given by ${N}_{{\rm{v}}}=A/(2\pi {{\ell }}^{2})$, where ${\ell }=\sqrt{{\hslash }/{qB}}$ is the magnetic length. The total number of atoms is given by $N={N}_{\uparrow }+{N}_{\downarrow }$, where N_α is the number of spin-α bosons.

In the Hamiltonian (1), we assume a contact interaction between atoms. For a gas tightly confined in a harmonic potential with frequency ${\omega }_{z}$ in the z direction, the effective coupling constants in the 2D plane are given by ${g}_{\alpha \alpha }={a}_{\alpha }\sqrt{8\pi {{\hslash }}^{3}{\omega }_{z}/M}$ and ${g}_{\uparrow \downarrow }\,=\,{g}_{\downarrow \uparrow }\,=\,{a}_{\uparrow \downarrow }\sqrt{8\pi {{\hslash }}^{3}{\omega }_{z}/M}$,³ where a_α and a_↑↓ are the s-wave scattering lengths between like and unlike bosons, respectively, in the 3D space. For simplicity, we set ${g}_{\uparrow \uparrow }\,=\,{g}_{\downarrow \downarrow }\,\equiv \,g\gt 0$ and ${N}_{\uparrow }={N}_{\downarrow }$ in the following discussions. We further assume that the synthetic magnetic fields B_α are sufficiently high or the interactions are sufficiently weak so that the energy scales of the interaction per atom, $| {g}_{\alpha \beta }| n$, are much smaller than the Landau-level spacing ${\hslash }{\omega }_{{\rm{c}}}:= {\hslash }{qB}/M$, where $n:= {N}_{\uparrow }/A={N}_{\downarrow }/A$ is the density of atoms of each component. In this situation, it is legitimate to employ the LLL approximation in which the Hilbert space is restricted to the lowest Landau level [10, 16, 17].

When the filling factor $\nu \equiv N/{N}_{v}$ is sufficiently high $(\nu \gg 1)$, the system is well described by the GP mean-field theory. In this theory, the GP energy functional $E[{\psi }_{\uparrow },{\psi }_{\downarrow }]$ is introduced by replacing the field operator ${\hat{\psi }}_{\alpha }({\bf{r}})$ by the condensate wave function ${\psi }_{\alpha }({\bf{r}})$ in the Hamiltonian (1); then, the functional is minimized under the conditions $\int {{\rm{d}}}^{2}{\bf{r}}| {\psi }_{\alpha }{| }^{2}={N}_{\alpha }$ ($\alpha =\uparrow ,\,\downarrow$) to determine the ground-state wave functions $\{{\psi }_{\alpha }({\bf{r}})\}$. Using the LLL wave functions which have periodic patterns of zeros and are equivalent to the LLL magnetic Bloch states described in section 2.2, Mueller and Ho [41] have obtained a rich ground-state phase diagram for the parallel-field case, which consists of five different vortex-lattice structures as shown in the upper panels of figure 1. Notably, the GP energy functionals for the parallel- and antiparallel-field cases are related to each other as ${E}_{\mathrm{antiparallel}}[{\psi }_{\uparrow },{\psi }_{\downarrow }]={E}_{\mathrm{parallel}}[{\psi }_{\uparrow },{\psi }_{\downarrow }^{* }]$ [45]. This implies that within the GP theory, the ground-state wave function of one case can be obtained from that of the other through the complex conjugation of the spin-$\downarrow$ component⁴ . Therefore, BECs in antiparallel fields also exhibit a rich variety of vortex-lattice structures as shown in figure 1 in the same way as BECs in parallel fields.

To describe the excitation properties of a vortex lattice, it is important to choose the basis states consistent with the periodicity of the lattice. Following [38, 57, 58], we utilize the LLL magnetic Bloch states for this purpose. Let a₁ and a₂ be the primitive vectors of a vortex lattice as shown in figure 1(f). These vectors satisfy

$$\begin{eqnarray}&&{({{\bf{a}}}_{1}\times {{\bf{a}}}_{2})}_{z}=2\pi {{\ell }}^{2}=A/{N}_{{\rm{v}}},\end{eqnarray} \tag{ 3 }$$

which implies the presence of one vortex in each component per unit cell. The reciprocal primitive vectors are then given by

$$\begin{eqnarray}&&{{\bf{b}}}_{1}=-{{\bf{e}}}_{z}\times {{\bf{a}}}_{2}/{{\ell }}^{2},\,{{\bf{b}}}_{2}={{\bf{e}}}_{z}\times {{\bf{a}}}_{1}/{{\ell }}^{2},\end{eqnarray} \tag{ 4 }$$

which satisfy ${{\bf{a}}}_{i}\cdot {{\bf{b}}}_{j}=2\pi {\delta }_{{ij}}\,(i,j=1,2)$. Using the pseudomomentum for a spin-α particle

$$\begin{eqnarray}&&{{\bf{K}}}_{\alpha }={\bf{p}}-q{{\bf{A}}}_{\alpha }+q{{\bf{B}}}_{\alpha }\times {\bf{r}}={\bf{p}}+{\epsilon }_{\alpha }\displaystyle \frac{{qB}}{2}{{\bf{e}}}_{z}\times {\bf{r}},\end{eqnarray} \tag{ 5 }$$

we introduce the magnetic translation operator as ${T}_{\alpha }({\bf{s}})={{\rm{e}}}^{-{\rm{i}}{{\bf{K}}}_{\alpha }\cdot {\bf{s}}/{\hslash }}$ [59]. We note that the pseudomomentum ${{\bf{K}}}_{\alpha }=({K}_{\alpha ,x},{K}_{\alpha ,y})$ satisfies the commutation relation [K_α,x, K_α,y] = −i_α ℏ²/ℓ². Starting from the most localized symmetric LLL wave function ${c}_{0}({\bf{r}})={{\rm{e}}}^{-{{\bf{r}}}^{2}/4{{\ell }}^{2}}/\sqrt{2\pi {{\ell }}^{2}}$, we construct a set of LLL wave functions by multiplying two translation operators as

$$\begin{eqnarray*}&&{c}_{{\bf{m}}\alpha }({\bf{r}})={T}_{\alpha }({m}_{1}{{\bf{a}}}_{1}){T}_{\alpha }({m}_{2}{{\bf{a}}}_{2}){c}_{0}({\bf{r}})=\displaystyle \frac{{(-1)}^{{m}_{1}{m}_{2}}}{\sqrt{2\pi {{\ell }}^{2}}}\exp \left[-\displaystyle \frac{1}{4{{\ell }}^{2}}{({\bf{r}}-{{\bf{r}}}_{m})}^{2}-\displaystyle \frac{{\rm{i}}{\epsilon }_{\alpha }}{2{{\ell }}^{2}}{({\bf{r}}\times {{\bf{r}}}_{m})}_{z}\right],\end{eqnarray*}$$

where ${{\bf{r}}}_{m}={m}_{1}{{\bf{a}}}_{1}+{m}_{2}{{\bf{a}}}_{2}$ with ${\bf{m}}=({m}_{1},{m}_{2})\in {{\mathbb{Z}}}^{2}$. Here, T_α(m₁a₁) and T_α(m₂a₂) commute with each other since every unit cell is pierced by one magnetic flux quantum as seen in equation (3); this property justifies the application of Bloch's theorem. By superposing ${c}_{{\bf{m}}\alpha }({\bf{r}})$ for ${N}_{{\rm{v}}}$ possible translations ${\bf{m}}$ on a torus, we can construct the LLL magnetic Bloch state as [57]

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{{\bf{k}}\alpha }({\bf{r}})=\displaystyle \frac{1}{\sqrt{{N}_{{\rm{v}}}\zeta ({\bf{k}})}}\displaystyle \sum _{{\bf{m}}}{c}_{{\bf{m}}\alpha }({\bf{r}}){{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {{\bf{r}}}_{{\rm{m}}}}\end{eqnarray} \tag{ 6 }$$

with the normalization factor

$$\begin{eqnarray}&&\zeta ({\bf{k}})=\displaystyle \sum _{{\bf{m}}}{(-1)}^{{m}_{1}{m}_{2}}{{\rm{e}}}^{-{{\bf{r}}}_{{\bf{m}}}^{2}/4{{\ell }}^{2}-{\rm{i}}{\bf{k}}\cdot {{\bf{r}}}_{{\bf{m}}}}.\end{eqnarray} \tag{ 7 }$$

This state is an eigenstate of ${T}_{\alpha }({{\bf{a}}}_{j})$ with an eigenvalue ${{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {{\bf{a}}}_{j}}$.

The LLL magnetic Bloch state ${{\rm{\Psi }}}_{{\bf{k}}\alpha }({\bf{r}})$ represents a vortex lattice with a periodic pattern of zeros for any value of the wave vector ${\bf{k}}$ ⁵ . Indeed, by rewriting equation (6) as

$$\begin{eqnarray*}\begin{array}{rcl}\sqrt{{N}_{{\rm{v}}}\zeta ({\bf{k}})}{{\rm{\Psi }}}_{{\bf{k}}\alpha }({\bf{r}}) & = & \displaystyle \sum _{{\bf{m}}}{c}_{{\bf{m}}\alpha }^{* }({\bf{r}})\exp \left[{\rm{i}}\left(-\displaystyle \frac{{\epsilon }_{\alpha }}{{{\ell }}^{2}}{{\bf{e}}}_{z}\times {\bf{r}}+{\bf{k}}\right)\cdot {{\bf{r}}}_{{\bf{m}}}\right]\\ & = & \displaystyle \sum _{{\bf{m}}}{c}_{{\bf{m}}\alpha }^{* }({\bf{r}})\exp \left\{-\displaystyle \frac{{\rm{i}}{\epsilon }_{\alpha }}{{{\ell }}^{2}}[{{\bf{e}}}_{z}\times ({\bf{r}}+{\epsilon }_{\alpha }{{\ell }}^{2}{{\bf{e}}}_{z}\times {\bf{k}})]\cdot {{\bf{r}}}_{{\bf{m}}}\right\}\end{array}\end{eqnarray*}$$

and comparing it with the complex conjugate of the Perelomov overcompleteness equation ${\sum }_{{\bf{m}}}{(-1)}^{{m}_{1}+{m}_{2}}{c}_{{\bf{m}}\alpha }({\bf{r}})=0$ [60], we find that ${{\rm{\Psi }}}_{{\bf{k}}\alpha }({\bf{r}})$ has zeros at [58]

$$\begin{eqnarray}&&{\bf{r}}={{\bf{r}}}_{{\bf{n}}}+\displaystyle \frac{1}{2}({{\bf{a}}}_{1}+{{\bf{a}}}_{2})-{\epsilon }_{\alpha }{{\ell }}^{2}{{\bf{e}}}_{z}\times {\bf{k}},\,{\bf{n}}=({n}_{1},{n}_{2})\in {{\mathbb{Z}}}^{2}.\end{eqnarray} \tag{ 8 }$$

When one describes a triangular vortex lattice of a scalar BEC using a LLL magnetic Bloch state, the choice of the wave vector ${\bf{k}}$ is arbitrary once the primitive vectors ${{\bf{a}}}_{1}$ and ${{\bf{a}}}_{2}$ are set appropriately. This is because a change in ${\bf{k}}$ only leads to a translation of zeros as seen in equation (8). The vortex lattices of two-component BECs in figure 1 can also be described by the LLL magnetic Bloch states ${{\rm{\Psi }}}_{{{\bf{q}}}_{\alpha },\alpha }({\bf{r}})\,(\alpha =\uparrow ,\downarrow );$ however, the wave vectors ${{\bf{q}}}_{\uparrow }$ and ${{\bf{q}}}_{\downarrow }$ have to be chosen in a way consistent with the displacement ${u}_{1}{{\bf{a}}}_{1}+{u}_{2}{{\bf{a}}}_{2}$ between the components (see figure 1(f)). One useful choice is

$$\begin{eqnarray}\begin{array}{rcl}{{\bf{q}}}_{\uparrow } & = & +\displaystyle \frac{{\epsilon }_{\uparrow }}{2{{\ell }}^{2}}{{\bf{e}}}_{z}\times ({u}_{1}{{\bf{a}}}_{1}+{u}_{2}{{\bf{a}}}_{2})=\displaystyle \frac{{\epsilon }_{\uparrow }}{2}(-{u}_{2}{{\bf{b}}}_{1}+{u}_{1}{{\bf{b}}}_{2}),\\ {{\bf{q}}}_{\downarrow } & = & -\displaystyle \frac{{\epsilon }_{\downarrow }}{2{{\ell }}^{2}}{{\bf{e}}}_{z}\times ({u}_{1}{{\bf{a}}}_{1}+{u}_{2}{{\bf{a}}}_{2})=\displaystyle \frac{{\epsilon }_{\downarrow }}{2}(+{u}_{2}{{\bf{b}}}_{1}-{u}_{1}{{\bf{b}}}_{2}).\end{array}\end{eqnarray} \tag{ 9 }$$

Here, we displace the spin-$\uparrow$ component by $\tfrac{1}{2}({u}_{1}{{\bf{a}}}_{1}+{u}_{2}{{\bf{a}}}_{2})$ and the spin-$\downarrow$ component by $-\tfrac{1}{2}({u}_{1}{{\bf{a}}}_{1}+{u}_{2}{{\bf{a}}}_{2})$ instead of displacing only one of the components. This is useful for avoiding zeros of the normalization factor $\zeta ({\bf{k}})$ at some high-symmetry points in the first Brillouin zone [57] ⁶ .

Using the magnetic Bloch states (6), we expand the field operator as ${\hat{\psi }}_{\alpha }({\bf{r}})={\sum }_{{\bf{k}}}{{\rm{\Psi }}}_{{\bf{k}}\alpha }({\bf{r}}){b}_{{\bf{k}}\alpha }$, where ${\bf{k}}$ runs over the first Brillouin zone, and ${b}_{{\bf{k}}\alpha }$ is a bosonic annihilation operator satisfying $[{b}_{{\bf{k}}\alpha },{b}_{{\bf{k}}^{\prime} \alpha ^{\prime} }^{\dagger }]={\delta }_{{\bf{kk}}^{\prime} }{\delta }_{\alpha \alpha ^{\prime} }$. Substituting this expansion into the Hamiltonian, we obtain

$$\begin{eqnarray}\begin{array}{rcl}H & = & {H}_{\mathrm{kin}}+{H}_{\mathrm{int}}\\ & = & \displaystyle \frac{{\hslash }{\omega }_{{\rm{c}}}}{2}({\hat{N}}_{\uparrow }+{\hat{N}}_{\downarrow })+\displaystyle \frac{1}{2}\displaystyle \sum _{\alpha ,\beta }\displaystyle \sum _{{{\bf{k}}}_{1},{{\bf{k}}}_{2},{{\bf{k}}}_{3},{{\bf{k}}}_{4}}{V}_{\alpha \beta }({{\bf{k}}}_{1},{{\bf{k}}}_{2},{{\bf{k}}}_{3},{{\bf{k}}}_{4}){b}_{{{\bf{k}}}_{1}\alpha }^{\dagger }{b}_{{{\bf{k}}}_{2}\beta }^{\dagger }{b}_{{{\bf{k}}}_{3}\beta }{b}_{{{\bf{k}}}_{4}\alpha },\end{array}\end{eqnarray} \tag{ 10 }$$

where ℏω_c/2 is the LLL single-particle zero-point energy and ${\hat{N}}_{\alpha }={\sum }_{{\bf{k}}}{b}_{{\bf{k}}\alpha }^{\dagger }{b}_{{\bf{k}}\alpha }$ is the number operator for the spin-α component. The interaction matrix element ${V}_{\alpha \beta }({{\bf{k}}}_{1},{{\bf{k}}}_{2},{{\bf{k}}}_{3},{{\bf{k}}}_{4})$ is given by

$$\begin{eqnarray}&&{V}_{\alpha \beta }({{\bf{k}}}_{1},{{\bf{k}}}_{2},{{\bf{k}}}_{3},{{\bf{k}}}_{4})={g}_{\alpha \beta }\int {{\rm{d}}}^{2}{\bf{r}}\,{{\rm{\Psi }}}_{{{\bf{k}}}_{1}\alpha }^{* }({\bf{r}}){{\rm{\Psi }}}_{{{\bf{k}}}_{2}\beta }^{* }({\bf{r}}){{\rm{\Psi }}}_{{{\bf{k}}}_{3}\beta }({\bf{r}}){{\rm{\Psi }}}_{{{\bf{k}}}_{4}\alpha }({\bf{r}}).\end{eqnarray} \tag{ 11 }$$

As described in appendix B, this matrix element is calculated to be

$$\begin{eqnarray}&&{V}_{\alpha \beta }({{\bf{k}}}_{1},{{\bf{k}}}_{2},{{\bf{k}}}_{3},{{\bf{k}}}_{4})={\delta }_{{{\bf{k}}}_{1}+{{\bf{k}}}_{2},{{\bf{k}}}_{3}+{{\bf{k}}}_{4}}^{{\rm{P}}}\displaystyle \frac{{g}_{\alpha \beta }}{2A}\displaystyle \frac{{S}_{\alpha \beta }({{\bf{k}}}_{1},{{\bf{k}}}_{2},{{\bf{k}}}_{3})}{\sqrt{\zeta ({{\bf{k}}}_{1})\zeta ({{\bf{k}}}_{2})\zeta ({{\bf{k}}}_{3})\zeta ({{\bf{k}}}_{4})}}.\end{eqnarray} \tag{ 12 }$$

Here, ${\delta }_{{\bf{kk}}^{\prime} }^{{\rm{P}}}:= {\sum }_{{\bf{G}}}{\delta }_{{\bf{k}},{\bf{k}}^{\prime} +{\bf{G}}}$ is the periodic Kronecker's delta with ${\bf{G}}$ running over the reciprocal lattice vectors. In the case of parallel fields, the function ${S}_{\alpha \beta }({{\bf{k}}}_{1},{{\bf{k}}}_{2},{{\bf{k}}}_{3})$ does not depend on α or β, and is given by

$$\begin{eqnarray}\begin{array}{rcl}S({{\bf{k}}}_{1},{{\bf{k}}}_{2},{{\bf{k}}}_{3}) & = & \displaystyle \sum _{{\bf{p}}\in \{0,1\}{}^{2}}{(-1)}^{{p}_{1}{p}_{2}}{{\rm{e}}}^{-{{\bf{r}}}_{{\bf{p}}}^{2}/4{{\ell }}^{2}+{\rm{i}}{{\bf{k}}}_{3}\cdot {{\bf{r}}}_{{\bf{p}}}}\tilde{\zeta }({{\bf{k}}}_{1}+{{\bf{k}}}_{2}-2{{\bf{k}}}_{3}+({{\bf{r}}}_{{\bf{p}}}\times {{\rm{e}}}_{z}-{\rm{i}}{{\bf{r}}}_{{\bf{p}}})/2{{\ell }}^{2})\\ & & \times \zeta ({{\bf{k}}}_{1}+({{\bf{r}}}_{{\bf{p}}}\times {{\rm{e}}}_{z}+{\rm{i}}{{\bf{r}}}_{{\bf{p}}})/4{{\ell }}^{2})\zeta ({{\bf{k}}}_{2}+({{\bf{r}}}_{{\bf{p}}}\times {{\rm{e}}}_{z}+{\rm{i}}{{\bf{r}}}_{{\bf{p}}})/4{{\ell }}^{2}),\end{array}\end{eqnarray} \tag{ 13 }$$

where

$$\begin{eqnarray}&&\tilde{\zeta }({\bf{k}}):= \displaystyle \sum _{{\bf{m}}}{{\rm{e}}}^{-{{\bf{r}}}_{{\bf{m}}}^{2}/2{{\ell }}^{2}-{\rm{i}}{\bf{k}}\cdot {{\bf{r}}}_{{\bf{m}}}}.\end{eqnarray} \tag{ 14 }$$

In the case of antiparallel fields, ${S}_{\alpha \beta }({{\bf{k}}}_{1},{{\bf{k}}}_{2},{{\bf{k}}}_{3})$ depends on α and β, and is given in terms of $S({{\bf{k}}}_{1},{{\bf{k}}}_{2},{{\bf{k}}}_{3})$ defined above by

$$\begin{eqnarray}\begin{array}{rcl}{S}_{\uparrow \uparrow }({{\bf{k}}}_{1},{{\bf{k}}}_{2},{{\bf{k}}}_{3}) & = & S({{\bf{k}}}_{1},{{\bf{k}}}_{2},{{\bf{k}}}_{3}),\ {S}_{\downarrow \downarrow }({{\bf{k}}}_{1},{{\bf{k}}}_{2},{{\bf{k}}}_{3})=S{(-{{\bf{k}}}_{1},-{{\bf{k}}}_{2},-{{\bf{k}}}_{3})}^{* },\\ {S}_{\uparrow \downarrow }({{\bf{k}}}_{1},{{\bf{k}}}_{2},{{\bf{k}}}_{3}) & = & S({{\bf{k}}}_{1},-{{\bf{k}}}_{3},-{{\bf{k}}}_{2}),\ {S}_{\downarrow \uparrow }({{\bf{k}}}_{1},{{\bf{k}}}_{2},{{\bf{k}}}_{3})=S{(-{{\bf{k}}}_{1},{{\bf{k}}}_{3},{{\bf{k}}}_{2})}^{* }.\end{array}\end{eqnarray} \tag{ 15 }$$

At high filling factors, the condensate is only weakly depleted and therefore we can apply the Bogoliubov approximation [36–38, 61]⁷ . Provided that the condensation occurs at the wave vector ${{\bf{q}}}_{\alpha }$ in the spin-α component, it is useful to introduce

$$\begin{eqnarray}&&{\tilde{b}}_{{\bf{k}}\alpha }:= {b}_{{{\bf{q}}}_{\alpha }+{\bf{k}},\alpha },\,{\tilde{V}}_{\alpha \beta }({{\bf{k}}}_{1},{{\bf{k}}}_{2},{{\bf{k}}}_{3},{{\bf{k}}}_{4}):= {V}_{\alpha \beta }({{\bf{q}}}_{\alpha }+{{\bf{k}}}_{1},{{\bf{q}}}_{\beta }+{{\bf{k}}}_{2},{{\bf{q}}}_{\beta }+{{\bf{k}}}_{3},{{\bf{q}}}_{\alpha }+{{\bf{k}}}_{4}).\end{eqnarray} \tag{ 16 }$$

By setting

$$\begin{eqnarray}&&{\tilde{b}}_{0\alpha }\simeq {\tilde{b}}_{0\alpha }^{\dagger }\simeq \sqrt{{N}_{\alpha }-\displaystyle \sum _{{\bf{k}}\ne 0}{\tilde{b}}_{{\bf{k}}\alpha }^{\dagger }{\tilde{b}}_{{\bf{k}}\alpha }}\end{eqnarray} \tag{ 17 }$$

and retaining terms up to the second order in ${\tilde{b}}_{{\bf{k}}\alpha }$ and ${\tilde{b}}_{{\bf{k}}\alpha }^{\dagger }$ (${\bf{k}}\ne 0$), we obtain the following Bogoliubov Hamiltonian:

$$\begin{eqnarray}\begin{array}{rcl}{H}_{\mathrm{int}} & = & \displaystyle \frac{1}{2}\displaystyle \sum _{\alpha ,\beta }{N}_{\alpha }{N}_{\beta }{\tilde{V}}_{\alpha \beta }(0,0,0,0)-\displaystyle \frac{1}{2}\displaystyle \sum _{{\bf{k}}\ne 0}\displaystyle \sum _{\alpha }[{h}_{\alpha }({\bf{k}})+{\omega }_{\alpha \alpha }({\bf{k}})]\\ & & +\displaystyle \frac{1}{2}\displaystyle \sum _{{\bf{k}}\ne 0}({\tilde{b}}_{{\bf{k}}\uparrow }^{\dagger },{\tilde{b}}_{{\bf{k}}\downarrow }^{\dagger },{\tilde{b}}_{-{\bf{k}},\uparrow },{\tilde{b}}_{-{\bf{k}},\downarrow }){ \mathcal M }({\bf{k}})\left(\begin{array}{c}{\tilde{b}}_{{\bf{k}}\uparrow }\\ {\tilde{b}}_{{\bf{k}}\downarrow }\\ {\tilde{b}}_{-{\bf{k}},\uparrow }^{\dagger }\\ {\tilde{b}}_{-{\bf{k}},\downarrow }^{\dagger }\end{array}\right).\end{array}\end{eqnarray} \tag{ 18 }$$

Here, the matrix ${ \mathcal M }({\bf{k}})$ is given by

$$\begin{eqnarray}{ \mathcal M }({\bf{k}})=\left(\begin{array}{cccc}{h}_{\uparrow }({\bf{k}})+{\omega }_{\uparrow \uparrow }({\bf{k}}) & {\omega }_{\uparrow \downarrow }({\bf{k}}) & {\lambda }_{\uparrow \uparrow }({\bf{k}}) & {\lambda }_{\uparrow \downarrow }({\bf{k}})\\ {\omega }_{\downarrow \uparrow }({\bf{k}}) & {h}_{\downarrow }({\bf{k}})+{\omega }_{\downarrow \downarrow }({\bf{k}}) & {\lambda }_{\downarrow \uparrow }({\bf{k}}) & {\lambda }_{\downarrow \downarrow }({\bf{k}})\\ {\lambda }_{\uparrow \uparrow }^{* }({\bf{k}}) & {\lambda }_{\downarrow \uparrow }^{* }({\bf{k}}) & {h}_{\uparrow }(-{\bf{k}})+{\omega }_{\uparrow \uparrow }(-{\bf{k}}) & {\omega }_{\downarrow \uparrow }(-{\bf{k}})\\ {\lambda }_{\uparrow \downarrow }^{* }({\bf{k}}) & {\lambda }_{\downarrow \downarrow }^{* }({\bf{k}}) & {\omega }_{\uparrow \downarrow }(-{\bf{k}}) & {h}_{\downarrow }(-{\bf{k}})+{\omega }_{\downarrow \downarrow }(-{\bf{k}})\end{array}\right),\end{eqnarray} \tag{ 19 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{h}_{\alpha }({\bf{k}}) & := & \displaystyle \sum _{\beta }{N}_{\beta }[{\tilde{V}}_{\alpha \beta }({\bf{k}},0,0,{\bf{k}})-{\tilde{V}}_{\alpha \beta }(0,0,0,0)],\\ {\omega }_{\alpha \beta }({\bf{k}}) & := & \sqrt{{N}_{\alpha }{N}_{\beta }}{\tilde{V}}_{\alpha \beta }({\bf{k}},0,{\bf{k}},0),\,{\lambda }_{\alpha \beta }({\bf{k}}):= \sqrt{{N}_{\alpha }{N}_{\beta }}{\tilde{V}}_{\alpha \beta }({\bf{k}},-{\bf{k}},0,0).\end{array}\end{eqnarray} \tag{ 20 }$$

To diagonalize the Bogoliubov Hamiltonian (18), we perform the Bogoliubov transformation

$$\begin{eqnarray}&&\left(\begin{array}{c}{\tilde{b}}_{{\bf{k}}\uparrow }\\ {\tilde{b}}_{{\bf{k}}\downarrow }\\ {\tilde{b}}_{-{\bf{k}},\uparrow }^{\dagger }\\ {\tilde{b}}_{-{\bf{k}},\downarrow }^{\dagger }\end{array}\right)=W({\bf{k}})\left(\begin{array}{c}{\gamma }_{{\bf{k}},1}\\ {\gamma }_{{\bf{k}},2}\\ {\gamma }_{-{\bf{k}},1}^{\dagger }\\ {\gamma }_{-{\bf{k}},2}^{\dagger }\end{array}\right),W({\bf{k}})=\left(\begin{array}{cc}{ \mathcal U }({\bf{k}}) & {{ \mathcal V }}^{* }(-{\bf{k}})\\ { \mathcal V }({\bf{k}}) & {{ \mathcal U }}^{* }(-{\bf{k}})\end{array}\right).\end{eqnarray} \tag{ 21 }$$

Here, $W({\bf{k}})$ is a paraunitary matrix satisfying

$$\begin{eqnarray}&&{W}^{\dagger }({\bf{k}}){\tau }_{3}W({\bf{k}})=W({\bf{k}}){\tau }_{3}{W}^{\dagger }({\bf{k}})={\tau }_{3}:= \mathrm{diag}(1,1,-1,-1),\end{eqnarray} \tag{ 22 }$$

which ensures the invariance of the bosonic commutation relation. If the matrix $W({\bf{k}})$ is chosen to satisfy

$$\begin{eqnarray}&&{W}^{\dagger }({\bf{k}}){ \mathcal M }({\bf{k}})W({\bf{k}})=\mathrm{diag}({E}_{1}({\bf{k}}),{E}_{2}({\bf{k}}),{E}_{1}(-{\bf{k}}),{E}_{2}(-{\bf{k}})),\end{eqnarray} \tag{ 23 }$$

the Bogoliubov Hamiltonian is diagonalized as

$$\begin{eqnarray}&&{H}_{\mathrm{int}}=\displaystyle \frac{1}{2}\displaystyle \sum _{\alpha ,\beta }{N}_{\alpha }{N}_{\beta }{\tilde{V}}_{\alpha \beta }(0,0,0,0)-\displaystyle \frac{1}{2}\displaystyle \sum _{{\bf{k}}\ne 0}\displaystyle \sum _{\alpha }[{h}_{\alpha }({\bf{k}})+{\omega }_{\alpha \alpha }({\bf{k}})]+\displaystyle \sum _{{\bf{k}}\ne 0}\displaystyle \sum _{i=1,2}{E}_{i}({\bf{k}})\left({\gamma }_{{\bf{k}}i}^{\dagger }{\gamma }_{{\bf{k}}i}+\displaystyle \frac{1}{2}\right).\end{eqnarray} \tag{ 24 }$$

By multiplying equation (23) from the left by $W({\bf{k}}){\tau }_{3}$ and using equation (22), one finds

$$\begin{eqnarray}&&{\tau }_{3}{ \mathcal M }({\bf{k}})W({\bf{k}})=W({\bf{k}})\mathrm{diag}({E}_{1}({\bf{k}}),{E}_{2}({\bf{k}}),-{E}_{1}(-{\bf{k}}),-{E}_{2}(-{\bf{k}})).\end{eqnarray} \tag{ 25 }$$

Therefore, the excitation energies ${E}_{i}({\bf{k}})\,(i=1,2)$ can be obtained as the right eigenvalues of ${\tau }_{3}{ \mathcal M }({\bf{k}})$.

With the Bogoliubov Hamiltonian (24), the field operator shows the following time evolution:

$$\begin{eqnarray}&&{\hat{\psi }}_{\alpha }({\bf{r}},t)\simeq \sqrt{{N}_{\alpha }}{{\rm{\Psi }}}_{{{\bf{q}}}_{\alpha },\alpha }({\bf{r}})+\displaystyle \sum _{{\bf{k}}\ne 0}{{\rm{\Psi }}}_{{{\bf{q}}}_{\alpha }+{\bf{k}},\alpha }({\bf{r}})\displaystyle \sum _{i=1,2}[{{ \mathcal U }}_{\alpha i}({\bf{k}}){{\rm{e}}}^{-{{\rm{i}}{E}}_{i}({\bf{k}})t/{\hslash }}{\gamma }_{{\bf{k}}i}+{{ \mathcal V }}_{\alpha i}^{* }(-{\bf{k}}){{\rm{e}}}^{{{\rm{i}}{E}}_{i}(-{\bf{k}})t/{\hslash }}{\gamma }_{-{\bf{k}},i}^{\dagger }].\end{eqnarray} \tag{ 26 }$$

If we replace ${\gamma }_{{\bf{k}}i}$ and ${\gamma }_{{\bf{k}}i}^{\dagger }$ by c-numbers, we may view this equation as the classical time evolution of a condensate wave function ${\psi }_{\alpha }({\bf{r}},t)$. In particular, by setting ${\gamma }_{{\bf{k}}i},{\gamma }_{{\bf{k}}i}^{\dagger }\to c\sqrt{{N}_{\alpha }}=c\sqrt{{nA}}\ne 0$ (with c being a real constant) for the specific mode $({\bf{k}},i)$, we obtain

$$\begin{eqnarray}&&\displaystyle \frac{{\psi }_{\alpha }({\bf{r}},t)}{\sqrt{n}}=\sqrt{A}{{\rm{\Psi }}}_{{{\bf{q}}}_{\alpha },\alpha }({\bf{r}})+c\sqrt{A}[{{\rm{\Psi }}}_{{{\bf{q}}}_{\alpha }+{\bf{k}},\alpha }({\bf{r}}){{ \mathcal U }}_{\alpha i}({\bf{k}}){{\rm{e}}}^{-{{\rm{i}}{E}}_{i}({\bf{k}})t/{\hslash }}+{{\rm{\Psi }}}_{{{\bf{q}}}_{\alpha }-{\bf{k}},\alpha }({\bf{r}}){{ \mathcal V }}_{\alpha i}^{* }({\bf{k}}){{\rm{e}}}^{{{\rm{i}}{E}}_{i}({\bf{k}})t/{\hslash }}].\end{eqnarray} \tag{ 27 }$$

This can be used to show how the density profiles $| {\psi }_{\alpha }({\bf{r}},t){| }^{2}/n\,(\alpha =\uparrow ,\downarrow )$ and the vortex positions change in time in the concerned mode $({\bf{k}},i)$. In doing so, it is useful to use the representation of $\sqrt{A}{{\rm{\Psi }}}_{{\bf{k}},\alpha }({\bf{r}})$ in terms of Jacobi's theta function (equation (A.2) in appendix A) as this function is supported in various computing systems⁸ .

We use the formulation described above to numerically calculate the Bogoliubov excitation spectrum $\{{E}_{i}({\bf{k}})\}$ in the following way. For a given wave vector ${\bf{k}}$, we calculate the matrix ${ \mathcal M }({\bf{k}})$ in equation (19) by using equations (12), (13), and (15). We note that each of the functions $\zeta (\cdot )$ and $\tilde{\zeta }(\cdot )$ used in equation (13) involves an infinite sum but only with respect to two integer variables (see equations (7) and (14)), which can numerically be taken with high accuracy. We then calculate the right eigenvalues of ${\tau }_{3}{ \mathcal M }({\bf{k}})$ to obtain $\{{E}_{i}({\bf{k}})\}$.

Figure 2 presents the obtained energy spectra for all the lattice structures in figure 1 and for both the parallel- and antiparallel-field cases. In all the cases, we find that there appear two modes with linear and quadratic dispersion relations at low energies around the Γ point. Furthermore, we find anisotropy of the coefficients of these dispersion relations. For example, such anisotropy can clearly be seen along the path M₁ → Γ → R for (c) rhombic, (d) square, and (e) rectangular lattices. We discuss such anisotropy in detail in later sections.

Zoom In Zoom Out Reset image size

To gain some physical insight into the low-energy excitation modes, we present in figure 3 the density profiles of the modes with quadratic (i = 2) and linear (i = 1) dispersion relations at ${\bf{k}}=(0.2a/{{\ell }}^{2},0)$ for (b) interlaced triangular lattices in parallel fields. As seen in this figure, vortices move perpendicularly to ${\bf{k}}$ relative to the ground state. Furthermore, spin-$\uparrow$ and $\downarrow$ vortices show in-phase (anti-phase) oscillations in the i = 2 (i = 1) mode. Specifically, around k_x = 0, both spin-$\uparrow$ and $\downarrow$ vortices move in the $-y$ direction in the i = 2 mode (upper panels of figure 3) while they move in opposite directions ($\mp y$) in the i = 1 mode (lower panels). Similar results are also obtained in the antiparallel-field case (not shown). These features are consistent with those obtained from the effective field theory described in section 3.

Zoom In Zoom Out Reset image size

Apart from the low-energy features, the spectra in figure 2 also exhibit unique structures of band touching at some high-symmetry points or along lines in the Brillouin zone. In particular, the spectra for (c) rhombic, (d) square, and (e) rectangular lattices in parallel fields exhibit line nodes, whose locations in the Brillouin zones are shown in figure 2(f). This can be understood as a consequence of a 'fractional' translation symmetry⁹ [62, 63]. Namely, in these cases, the system is invariant under the product ${{ \mathcal T }}^{({\rm{P}})}$ of the translation by ${{\bf{a}}}_{3}/2$ and the spin reversal $\uparrow \,\leftrightarrow \,\downarrow$, where ${{\bf{a}}}_{3}:= {{\bf{a}}}_{1}+{{\bf{a}}}_{2}$. Since the unitary operator ${{ \mathcal T }}^{({\rm{P}})}$ commutes with the Bogoliubov Hamiltonian and ${({{ \mathcal T }}^{({\rm{P}})})}^{2}$ gives the translation by ${{\bf{a}}}_{3}$, the Bloch states at ${\bf{k}}$ can be chosen to be the eigenstates of ${{ \mathcal T }}^{({\rm{P}})}$ with ${{ \mathcal T }}^{({\rm{P}})}| {w}_{{\bf{k}}}^{\pm }\rangle =\pm {{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {{\bf{a}}}_{3}/2}| {w}_{{\bf{k}}}^{\pm }\rangle$. For a smooth change ${\bf{k}}\to {\bf{k}}+{{\bf{b}}}_{i}\,(i=1,2)$, the two eigenstates must be swapped, indicating the occurrence of an odd number of degeneracies. In figure 2(f), we can indeed confirm that starting from any point other than the line nodes, the degeneracy occurs once or three times for the above changes of ${\bf{k}}$. The emergence of point nodes at the M₁ and M₂ points for the same lattices ((c)–(e)) in antiparallel fields can be understood by considering the symmetry under the product ${{ \mathcal T }}^{(\mathrm{AP})}$ of the time reversal and the translation by ${{\bf{a}}}_{3}/2$. Since ${({{ \mathcal T }}^{(\mathrm{AP})})}^{2}$ is equal to the translation by ${{\bf{a}}}_{3}$, we have ${({{ \mathcal T }}^{(\mathrm{AP})})}^{2}={{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {{\bf{a}}}_{3}}$ in the subspace with the wave vector ${\bf{k}}$. The Kramers degeneracy thus occurs at time-reversal-invariant momenta with ${{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {{\bf{a}}}_{3}}\ne 1$, which is the case for ${\bf{k}}={{\bf{b}}}_{1}/2$ and ${{\bf{b}}}_{2}/2$ (M₁ and M₂ points). In appendix D, we further discuss some other features of the spectra at high-symmetry points, such as the coincidence of the excitation energies between the two types of fields at the M₁ and M₂ points in figures 2(c)–(e) by using the numerical data of the Bogoliubov Hamiltonian matrix ${ \mathcal M }({\bf{k}})$ and the density profiles of the excitation modes.

The Lagrangian density of the two-component BECs corresponding to the Hamiltonian (1) is given by [61]

$$\begin{eqnarray}&&{ \mathcal L }=\displaystyle \sum _{\alpha }\left[\displaystyle \frac{{\rm{i}}{\hslash }}{2}({\psi }_{\alpha }^{\dagger }{\dot{\psi }}_{\alpha }-{\dot{\psi }}_{\alpha }^{\dagger }{\psi }_{\alpha })-\displaystyle \frac{1}{2M}| (-{\rm{i}}{\hslash }{\rm{\nabla }}-q{{\bf{A}}}_{\alpha }){\psi }_{\alpha }{| }^{2}\right]-\displaystyle \sum _{\alpha ,\beta }\displaystyle \frac{{g}_{\alpha \beta }}{2}| {\psi }_{\alpha }{| }^{2}| {\psi }_{\beta }{| }^{2},\end{eqnarray} \tag{ 28 }$$

where ${\psi }_{\alpha }({\bf{r}},t)$ is the bosonic field for the spin-α component. To describe the low-energy properties of the BECs, it is useful to decompose the field as ${\psi }_{\alpha }=\sqrt{{n}_{\alpha }}\exp (-{\rm{i}}{\theta }_{\alpha })$, where ${n}_{\alpha }({\bf{r}},t)$ and ${\theta }_{\alpha }({\bf{r}},t)$ are the density and phase variables, respectively. Substituting this into equation (28) and keeping only the leading terms in the derivative expansion, we obtain

$$\begin{eqnarray}&&{ \mathcal L }={\mu }_{\uparrow }{n}_{\uparrow }+{\mu }_{\downarrow }{n}_{\downarrow }-\displaystyle \frac{g}{2}({n}_{\uparrow }^{2}+{n}_{\downarrow }^{2})-{g}_{\uparrow \downarrow }{n}_{\uparrow }{n}_{\downarrow },\end{eqnarray} \tag{ 29 }$$

where

$$\begin{eqnarray}&&{\mu }_{\alpha }={\hslash }{\dot{\theta }}_{\alpha }-\displaystyle \frac{1}{2M}{({\hslash }{\rm{\nabla }}{\theta }_{\alpha }+q{{\bf{A}}}_{\alpha })}^{2}\end{eqnarray} \tag{ 30 }$$

is an effective chemical potential for the spin-α component. Introducing ${n}_{\pm }:= {n}_{\uparrow }\pm {n}_{\downarrow }$ and ${g}_{\pm }:= g\pm {g}_{\uparrow \downarrow }$, we can rewrite equation (29) as

$$\begin{eqnarray}&&{ \mathcal L }=-\displaystyle \frac{{g}_{+}}{4}{n}_{+}^{2}-\displaystyle \frac{{g}_{-}}{4}{n}_{-}^{2}+\displaystyle \frac{{\mu }_{\uparrow }+{\mu }_{\downarrow }}{2}{n}_{+}+\displaystyle \frac{{\mu }_{\uparrow }-{\mu }_{\downarrow }}{2}{n}_{-}.\end{eqnarray} \tag{ 31 }$$

By integrating out ${n}_{\pm }({\bf{r}},t)$, we obtain the effective Lagrangian for the phase variables $\{{\theta }_{\alpha }({\bf{r}},t)\}$ as

$$\begin{eqnarray}&&{ \mathcal L }=\displaystyle \frac{1}{4{g}_{+}}{({\mu }_{\uparrow }+{\mu }_{\downarrow })}^{2}+\displaystyle \frac{1}{4{g}_{-}}{({\mu }_{\uparrow }-{\mu }_{\downarrow })}^{2}.\end{eqnarray} \tag{ 32 }$$

In the presence of vortices, the phase variables $\{{\theta }_{\alpha }({\bf{r}},t)\}$ involve singularities. It is thus useful to decompose θ_α into regular and singular parts as θ_α = θ_reg,α + θ_sing,α. Since the singular part θ_sing,α varies rapidly in space, it is not a convenient variable for a coarse-grained description over long length scales. To describe the long-wavelength physics, it is useful to start from the vortex-lattice ground state (as in figure 1) and to consider small displacement of vortices from the equilibrium positions. Specifically, we introduce the displacement vector field u_α(r, t) = r − X_α(r, t), where r is the equilibrium position of the vortex and X_α is the position at time t. The derivatives of the singular part θ_sing,α of the phase are related to the displacement ${{\bf{u}}}_{\alpha }$ as [39]

$$\begin{eqnarray*}&&{\hslash }{\dot{\theta }}_{\mathrm{sing},\alpha }=-\displaystyle \frac{{{qB}}_{\alpha }}{2}{({{\bf{u}}}_{\alpha }\times {\dot{{\bf{u}}}}_{\alpha })}_{z},\,{\hslash }{\rm{\nabla }}{\theta }_{\mathrm{sing},\alpha }+q{{\bf{A}}}_{\alpha }={{qB}}_{\alpha }{{\bf{e}}}_{z}\times {{\bf{u}}}_{\alpha }-\displaystyle \frac{{{qB}}_{\alpha }}{2}\displaystyle \sum _{i,j}{\epsilon }_{{ij}}{u}_{\alpha }^{i}{\rm{\nabla }}{u}_{\alpha }^{j},\end{eqnarray*}$$

where _ij is an antisymmetric tensor with _xy = − _yx = +1. The effective chemical potential in equation (30) can then be expressed in terms of $\{{\theta }_{\mathrm{reg},\alpha },{{\bf{u}}}_{\alpha }\}$ as

$$\begin{eqnarray*}&&{\mu }_{\alpha }={\hslash }{\dot{\theta }}_{\mathrm{reg},\alpha }-\displaystyle \frac{{{qB}}_{\alpha }}{2}{({{\bf{u}}}_{\alpha }\times {\dot{{\bf{u}}}}_{\alpha })}_{z}-\displaystyle \frac{1}{2M}{\left({\hslash }{\rm{\nabla }}{\theta }_{\mathrm{reg},\alpha }+{{qB}}_{\alpha }{{\bf{e}}}_{z}\times {{\bf{u}}}_{\alpha }-\displaystyle \frac{{{qB}}_{\alpha }}{2}\displaystyle \sum _{{ij}}{\epsilon }_{{ij}}{u}_{\alpha }^{i}{\rm{\nabla }}{u}_{\alpha }^{j}\right)}^{2}.\end{eqnarray*}$$

One should also note that the displacement ${{\bf{u}}}_{\alpha }({\bf{r}},t)$ leads to a change in the elastic energy $\int {{\rm{d}}}^{2}{\bf{r}}\,{{ \mathcal E }}_{\mathrm{el}}({{\bf{u}}}_{\alpha },{\partial }_{i}{{\bf{u}}}_{\alpha })$. Here, the form of the elastic energy density ${{ \mathcal E }}_{\mathrm{el}}$ depends on the type of a lattice as discussed in the next section and appendix E. The effective Lagrangian in terms of $\{{\theta }_{\mathrm{reg},\alpha },{{\bf{u}}}_{\alpha }\}$ is then obtained as

$$\begin{eqnarray}&&{{ \mathcal L }}_{\mathrm{eff}}=\displaystyle \frac{1}{4{g}_{+}}{({\mu }_{\uparrow }+{\mu }_{\downarrow })}^{2}+\displaystyle \frac{1}{4{g}_{-}}{({\mu }_{\uparrow }-{\mu }_{\downarrow })}^{2}-{{ \mathcal E }}_{\mathrm{el}}.\end{eqnarray} \tag{ 33 }$$

Here, the difference from equation (32) occurs because the rapidly varying {θ_sing,α} have been replaced by the slowly varying $\{{{\bf{u}}}_{\alpha }\}$ via coarse graining.

The ground state of $H-{\mu }_{0}({N}_{\uparrow }+{N}_{\downarrow })$ is given by ${\theta }_{\mathrm{reg},{\alpha }}={\mu }_{0}t/{\hslash }$ and ${{\bf{u}}}_{\alpha }=0$. To discuss the low-energy properties, it is therefore useful to introduce ${\varphi }_{\alpha }={\mu }_{0}t/{\hslash }-{\theta }_{\mathrm{reg},\alpha }$ and expand the Lagrangian (33) in terms of $\{{\varphi }_{\alpha },{{\bf{u}}}_{\alpha }\}$. Keeping only the quadratic terms in these variables, we obtain

$$\begin{eqnarray}&&{{ \mathcal L }}_{\mathrm{eff}}=\displaystyle \frac{{{\hslash }}^{2}{\dot{\varphi }}_{+}^{2}}{4{g}_{+}}+\displaystyle \frac{{{\hslash }}^{2}{\dot{\varphi }}_{-}^{2}}{4{g}_{-}}-\displaystyle \frac{{\mu }_{0}}{{g}_{+}}\displaystyle \sum _{\alpha }\left[\displaystyle \frac{{{qB}}_{\alpha }}{2}{({{\bf{u}}}_{\alpha }\times {\dot{{\bf{u}}}}_{\alpha })}_{z}+\displaystyle \frac{1}{2M}{({\hslash }{{\bf{e}}}_{z}\times {\rm{\nabla }}{\varphi }_{\alpha }+{{qB}}_{\alpha }{{\bf{u}}}_{\alpha })}^{2}\right]-{{ \mathcal E }}_{\mathrm{el}},\end{eqnarray} \tag{ 34 }$$

where ${\varphi }_{\pm }:= {\varphi }_{1}\pm {\varphi }_{2}$. Because $\{{{\bf{u}}}_{\alpha }\}$ have the mass term $-{{\bf{u}}}_{\alpha }^{2}$, one can expect that they can safely be integrated out in the discussion of low-energy dynamics. To do so, it is useful to derive the Euler–Lagrange equations for $\{{{\bf{u}}}_{\alpha }\}$:

$$\begin{eqnarray}&&{{\bf{u}}}_{\alpha }+{\epsilon }_{\alpha }{{\ell }}^{2}{{\bf{e}}}_{z}\times {\rm{\nabla }}{\varphi }_{\alpha }-\displaystyle \frac{{\epsilon }_{\alpha }}{{\omega }_{{\rm{c}}}}{{\bf{e}}}_{z}\times \dot{{\bf{u}}}+\displaystyle \frac{{g}_{+}{{\ell }}^{2}}{{\mu }_{0}{\hslash }{\omega }_{{\rm{c}}}}\left[\displaystyle \frac{\partial {{ \mathcal E }}_{\mathrm{el}}}{\partial {{\bf{u}}}_{\alpha }}-\displaystyle \sum _{j}{\partial }_{j}\left(\displaystyle \frac{\partial {{ \mathcal E }}_{\mathrm{el}}}{\partial ({\partial }_{j}{{\bf{u}}}_{\alpha })}\right)\right]=0,\end{eqnarray} \tag{ 35 }$$

where we use the cyclotron frequency ω_c = qB/M and the magnetic length ${\ell }=\sqrt{{\hslash }/{qB}}$. The third and fourth terms on the left-hand side can be ignored in the LLL approximation (${\hslash }\omega ,| {g}_{\alpha \beta }| n\ll {\hslash }{\omega }_{{\rm{c}}}$, where ω is the frequency of our interest). Similar relations are also found in hydrodynamic theory [25–29, 31–33, 44]. Introducing ${{\bf{u}}}_{\pm }:= {{\bf{u}}}_{\uparrow }\pm {{\bf{u}}}_{\downarrow }$, equation (35) can be rewritten as

$$\begin{eqnarray}{{\bf{u}}}_{\pm }=\left\{\begin{array}{ll} & -{{\ell }}^{2}{{\bf{e}}}_{z}\times {\rm{\nabla }}{\varphi }_{\pm }\qquad (\mathrm{parallel}\ \mathrm{fields});\\ & -{{\ell }}^{2}{{\bf{e}}}_{z}\times {\rm{\nabla }}{\varphi }_{\mp }\qquad (\mathrm{antiparallel}\ \mathrm{fields}).\end{array}\right.\end{eqnarray} \tag{ 36 }$$

These relations indicate that the vortex displacements ${{\bf{u}}}_{\pm }$ and the phases ${\varphi }_{\pm }$ are coupled in an opposite manner between the parallel- and antiparallel-field cases. Namely, the symmetric u₊ (antisymmetric u₋) is coupled to the symmetric φ₊ (antisymmetric φ₋) in parallel fields, while they are coupled in a crossed manner in antiparallel fields. Equation (36) also indicates that the vortex displacement is perpendicular to the wave vector ${\bf{k}}$, which is consistent with the results shown in figure 3.

Substituting equation (36) into the Lagrangian density (34), we obtain the Lagrangian density in terms of {φ_±}, which can be used to determine the excitation spectrum. For this purpose, we need to determine the form of the elastic energy density ${{ \mathcal E }}_{\mathrm{el}}$, which is done next.

Since the elastic energy is invariant under a uniform change in ${{\bf{u}}}_{+}({\bf{r}},t)$ (i.e. translation of the lattices), ${{ \mathcal E }}_{\mathrm{el}}$ should be a function of ${\partial }_{i}{{\bf{u}}}_{+}(i=x,y)$ and ${{\bf{u}}}_{-}$ to the leading order in the derivative expansion. We therefore introduce the form

$$\begin{eqnarray}&&{{ \mathcal E }}_{\mathrm{el}}={{ \mathcal E }}_{\mathrm{el}}^{(+)}({\partial }_{i}{{\bf{u}}}_{+})+{{ \mathcal E }}_{\mathrm{el}}^{(-)}({{\bf{u}}}_{-})+{{ \mathcal E }}_{\mathrm{el}}^{(+-)}({\partial }_{i}{{\bf{u}}}_{+},{{\bf{u}}}_{-}).\end{eqnarray} \tag{ 37 }$$

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

$$\begin{eqnarray}&&{w}_{0}:= {\partial }_{x}{u}_{+}^{x}+{\partial }_{y}{u}_{+}^{y},\,{w}_{1}:= {\partial }_{x}{u}_{+}^{x}-{\partial }_{y}{u}_{+}^{y},\,{w}_{2}:= {\partial }_{y}{u}_{+}^{x}+{\partial }_{x}{u}_{+}^{y}.\end{eqnarray} \tag{ 38 }$$

In the LLL regime, the vortex density stays constant, and therefore w₀ = 0; this can also be confirmed by using equation (36). From a symmetry consideration (see appendix E), each term in equation (37) can be expressed as

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal E }}_{\mathrm{el}}^{(+)}({\partial }_{i}{{\bf{u}}}_{+}) & = & \displaystyle \frac{{{gn}}^{2}}{2}({C}_{1}{w}_{1}{}^{2}+{C}_{2}{w}_{2}{}^{2}+{C}_{3}{w}_{1}{w}_{2}),\\ {{ \mathcal E }}_{\mathrm{el}}^{(-)}({{\bf{u}}}_{-}) & = & \displaystyle \frac{{{gn}}^{2}}{2{{\ell }}^{2}}[{D}_{1}{({u}_{-}^{x})}^{2}+{D}_{2}{({u}_{-}^{y})}^{2}+{D}_{3}{u}_{-}^{x}{u}_{-}^{y}],\\ {{ \mathcal E }}_{\mathrm{el}}^{(+-)}({\partial }_{i}{{\bf{u}}}_{+},{{\bf{u}}}_{-}) & = & \displaystyle \frac{{{gn}}^{2}}{2{\ell }}{F}_{1}({w}_{1}{u}_{-}^{y}+{w}_{2}{u}_{-}^{x}),\end{array}\end{eqnarray} \tag{ 39 }$$

where $n:= {N}_{\uparrow }/A={N}_{\downarrow }/A$ is the average number density of each component. For each of the vortex lattices in figures 1(a)–(e), the dimensionless elastic constants {C₁, C₂, C₃, D₁, D₂, D₃, F₁} satisfy

$$\begin{eqnarray}&&\begin{array}{l}({\rm{a}})\,{C}_{1}={C}_{2}\equiv C\gt 0,{D}_{1}={D}_{2}\equiv D\gt 0,{C}_{3}={D}_{3}={F}_{1}=0;\\ ({\rm{b}})\,{C}_{1}={C}_{2}\equiv C\gt 0,{D}_{1}={D}_{2}\equiv D\gt 0,{C}_{3}={D}_{3}=0,{F}_{1}\ne 0;\\ ({\rm{c}})\,{C}_{1},{C}_{2},{D}_{1},{D}_{2}\gt 0,{C}_{3},{D}_{3}\ne 0,{F}_{1}=0;\\ ({\rm{d}})\,{C}_{1},{C}_{2}\gt 0,{D}_{1}={D}_{2}\equiv D\gt 0,{C}_{3}={D}_{3}={F}_{1}=0;\\ ({\rm{e}})\,{C}_{1},{C}_{2}\gt 0,{D}_{1},{D}_{2}\gt 0,{C}_{3}={D}_{3}={F}_{1}=0.\end{array}\end{eqnarray} \tag{ 40 }$$

Keçeli and Oktel [44] have considered an elastic energy consisting of ${{ \mathcal E }}_{\mathrm{el}}^{(\pm )}$ above, but have not included ${{ \mathcal E }}_{\mathrm{el}}^{(+-)}$. Therefore, in their work, the symmetric and antisymmetric displacements ${{\bf{u}}}_{\pm }$ are decoupled from each other in collective modes. In our analysis in appendix E, ${{ \mathcal E }}_{\mathrm{el}}^{(+-)}$ is found to be allowed by symmetry for interlaced triangular lattices. As shown below, this part crucially changes the low-energy spectrum, and explains the anisotropy of the spectrum for the concerned lattice structure.

We note that within the mean-field theory, the elastic energy density ${{ \mathcal E }}_{\mathrm{el}}$ should take the same form (equations (37) and (39)) for the parallel- and antiparallel-field cases because of the exact correspondence of the GP energy functionals between the two cases [45]. The dimensionless elastic constants are also expected to take the same values between the two cases. However, as we will see in section 4, the elastic constants estimated from the numerical results of the energy spectra are different between the two cases. We discuss this puzzling issue in section 4.2.

The Lagrangian density in terms of φ_± is obtained by substituting equation (36) into equation (34) and using the above ${{ \mathcal E }}_{\mathrm{el}}$. After performing the Fourier transformation ${\varphi }_{\pm }({\bf{r}},t)={\sum }_{{\bf{k}}}\int \tfrac{{\rm{d}}\omega }{2\pi \sqrt{A}}{{\rm{e}}}^{{\rm{i}}({\bf{k}}\cdot {\bf{r}}-\omega t)}{\varphi }_{\pm }({\bf{k}},\omega )$, we obtain the action

$$\begin{eqnarray}&&S=\displaystyle \sum _{{\bf{k}}}\int \displaystyle \frac{{\rm{d}}\omega }{2\pi }\displaystyle \frac{1}{2}({\varphi }_{+}(-{\bf{k}},-\omega ),{\varphi }_{-}(-{\bf{k}},-\omega )){\rm{i}}{G}{({\bf{k}},\omega )}^{-1}\left(\begin{array}{c}{\varphi }_{+}({\bf{k}},\omega )\\ {\varphi }_{-}({\bf{k}},\omega )\end{array}\right),\end{eqnarray} \tag{ 41 }$$

where

$$\begin{eqnarray}{\rm{i}}{G}{({\bf{k}},\omega )}^{-1}=\left(\begin{array}{cc}\displaystyle \frac{{{\hslash }}^{2}{\omega }^{2}}{2{g}_{+}}-{{\rm{\Gamma }}}_{\pm }({\bf{k}}) & \pm {\rm{i}}{\rm{\Gamma }}({\bf{k}})\\ \mp {\rm{i}}{\rm{\Gamma }}({\bf{k}}) & \displaystyle \frac{{{\hslash }}^{2}{\omega }^{2}}{2{g}_{-}}-{{\rm{\Gamma }}}_{\mp }({\bf{k}})\end{array}\right)\end{eqnarray} \tag{ 42 }$$

is the inverse of Green's function in Fourier space with

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Gamma }}}_{+}({\bf{k}}) & = & {{gn}}^{2}{{\ell }}^{4}[{C}_{1}{(2{k}_{x}{k}_{y})}^{2}+{C}_{2}{({k}_{x}^{2}-{k}_{y}^{2})}^{2}-{C}_{3}(2{k}_{x}{k}_{y})({k}_{x}^{2}-{k}_{y}^{2})],\\ {{\rm{\Gamma }}}_{-}({\bf{k}}) & = & {{gn}}^{2}{{\ell }}^{2}({D}_{1}{k}_{y}^{2}+{D}_{2}{k}_{x}^{2}-{D}_{3}{k}_{x}{k}_{y}),\\ {\rm{\Gamma }}({\bf{k}}) & = & \displaystyle \frac{1}{2}{{gn}}^{2}{{\ell }}^{3}{F}_{1}(3{k}_{x}^{2}{k}_{y}-{k}_{y}^{3}).\end{array}\end{eqnarray} \tag{ 43 }$$

In equation (42) (and equations (44), (45), and (47) below), the upper and lower of the double signs correspond to the parallel- and antiparallel-field cases, respectively.

The excitation spectrum corresponds to the poles of the Green's function, and can thus be obtained by solving the equation $\det [{\rm{i}}{G}{({\bf{k}},\omega )}^{-1}]=0$. Since ${{\rm{\Gamma }}}_{-}({\bf{k}})\gg {\rm{\Gamma }}({\bf{k}})\gg {{\rm{\Gamma }}}_{+}({\bf{k}})$ for $k{\ell }\ll 1$, we obtain the low-energy dispersion relations as

$$\begin{eqnarray}&&{E}_{2}({\bf{k}})=\sqrt{2{g}_{\pm }\left({{\rm{\Gamma }}}_{+}({\bf{k}})-\displaystyle \frac{{\rm{\Gamma }}{({\bf{k}})}^{2}}{{{\rm{\Gamma }}}_{-}({\bf{k}})}\right)},\,{E}_{1}({\bf{k}})=\sqrt{2{g}_{\mp }{{\rm{\Gamma }}}_{-}({\bf{k}})}.\end{eqnarray} \tag{ 44 }$$

Using equation (43) and the fact that ${{\rm{\Gamma }}}_{-}({\bf{k}})$ is isotropic when ${F}_{1}\ne 0$ (see equation (40)), we obtain the following explicit expressions

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{E}_{2}({\bf{k}})}{\sqrt{2}{gn}} & = & {\left(\displaystyle \frac{{g}_{\pm }}{g}\right)}^{\tfrac{1}{2}}{{\ell }}^{2}{\left[{C}_{1}{(2{k}_{x}{k}_{y})}^{2}+{C}_{2}{({k}_{x}^{2}-{k}_{y}^{2})}^{2}-{C}_{3}(2{k}_{x}{k}_{y})({k}_{x}^{2}-{k}_{y}^{2})-{C}_{4}\displaystyle \frac{{(3{k}_{x}^{2}{k}_{y}-{k}_{y}^{3})}^{2}}{{k}^{2}}\right]}^{\tfrac{1}{2}},\\ \displaystyle \frac{{E}_{1}({\bf{k}})}{\sqrt{2}{gn}} & = & {\left(\displaystyle \frac{{g}_{\mp }}{g}\right)}^{\tfrac{1}{2}}{\ell }{({D}_{1}{k}_{y}^{2}+{D}_{2}{k}_{x}^{2}-{D}_{3}{k}_{x}{k}_{y})}^{1/2}\end{array}\end{eqnarray} \tag{ 45 }$$

with ${C}_{4}:= {F}_{1}{}^{2}/4{D}_{1}$. We thus find the emergence of quadratic and linear dispersion relations whose anisotropy reflects the symmetry of each lattice structure. Furthermore, we find that the modes with the quadratic and linear dispersion relations originate mainly from the symmetric and antisymmetric parts ${{\bf{u}}}_{\pm }$ of the vortex displacement, respectively (we, however, note that these two parts are mixed slightly in the case of interlaced triangular lattices owing to ${F}_{1}\ne 0$). This explains the in-phase (anti-phase) oscillations of the i = 2 (i = 1) mode found in figure 3.

To discuss the anisotropy further, we parametrize the wave vector in terms of polar coordinates as ${\bf{k}}=k(\cos \theta ,\sin \theta )\,(k{\ell }\ll 1)$ and introduce the dimensionless functions {f_i(θ)} via

$$\begin{eqnarray}&&{E}_{i}({\bf{k}})=\sqrt{2}{gn}{(k{\ell })}^{i}{f}_{i}(\theta ),\,i=1,2.\end{eqnarray} \tag{ 46 }$$

Using the dispersion relations (45) obtained from the effective field theory, these functions are calculated as

$$\begin{eqnarray}\begin{array}{rcl}{f}_{2}(\theta ) & = & \sqrt{\displaystyle \frac{{g}_{\pm }}{g}}{[{C}_{1}{\sin }^{2}(2\theta )+{C}_{2}{\cos }^{2}(2\theta )-{C}_{3}\sin (2\theta )\cos (2\theta )-{C}_{4}{\sin }^{2}(3\theta )]}^{1/2},\\ {f}_{1}(\theta ) & = & \sqrt{\displaystyle \frac{{g}_{\mp }}{g}}{({D}_{1}{\sin }^{2}\theta +{D}_{2}{\cos }^{2}\theta -{D}_{3}\sin \theta \cos \theta )}^{1/2}.\end{array}\end{eqnarray} \tag{ 47 }$$

In this result (and also in equations (44) and (45)), the dependence on the type of synthetic fields occurs only in the coefficients $\sqrt{{g}_{\pm }/g}$. This observation leads to the following remarkable relations:

$$\begin{eqnarray}&&{f}_{2}^{{\rm{P}}}(\theta )\sqrt{\displaystyle \frac{g}{{g}_{+}}}={f}_{2}^{\mathrm{AP}}(\theta )\sqrt{\displaystyle \frac{g}{{g}_{-}}},\,{f}_{1}^{{\rm{P}}}(\theta )\sqrt{\displaystyle \frac{g}{{g}_{-}}}={f}_{1}^{\mathrm{AP}}(\theta )\sqrt{\displaystyle \frac{g}{{g}_{+}}},\end{eqnarray} \tag{ 48 }$$

where the superscripts P and AP refer to the parallel- and antiparallel-field cases, respectively. Namely, the functions $\{{f}_{i}^{{\rm{P}}/\mathrm{AP}}(\theta )\}$ for the two types of synthetic fields are related to each other by simple rescaling. While these rescaling relations are expected for all the lattice structures within the effective field theory, we show in the next section that the relations hold only for overlapping triangular lattices and break down for the other lattices.

For (a) overlapping triangular lattices, by using equations (40) and (47), the analytic expressions of $\{{f}_{i}^{{\rm{P}}/\mathrm{AP}}(\theta )\}$ for parallel (P) and antiparallel (AP) fields are obtained as

$$\begin{eqnarray}&&{f}_{2}^{{\rm{P}}}(\theta )=\sqrt{\displaystyle \frac{{g}_{+}}{g}C},\,{f}_{1}^{{\rm{P}}}(\theta )=\sqrt{\displaystyle \frac{{g}_{-}}{g}D},\,{f}_{2}^{\mathrm{AP}}(\theta )=\sqrt{\displaystyle \frac{{g}_{-}}{g}C},\,{f}_{1}^{\mathrm{AP}}(\theta )=\sqrt{\displaystyle \frac{{g}_{-}}{g}D}.\end{eqnarray} \tag{ 49 }$$

Notably, these functions show no dependence on θ in the effective field theory.

In numerical calculations, we obtain $\{{f}_{i}^{{\rm{P}}/\mathrm{AP}}(\theta )\}$ from the data of the Bogoliubov excitation spectra along a circular path ${\bf{k}}=k(\cos \theta ,\sin \theta )$ with sufficiently small k and arbitrary θ ∈ [0, 2π). Figure 4(a) presents numerical results for ${g}_{\uparrow \downarrow }/g=-0.2$. We find that the functions $\{{f}_{i}^{{\rm{P}}/\mathrm{AP}}(\theta )\}$ stay constant to a good accuracy consistent with the analytical expressions (49). The figure also shows the rescaled functions (defined by the left- and right-hand sides of equation (48)), clearly demonstrating the rescaling relations (48). The dimensionless elastic constants C and D thus take the same values for the two types of fields and are plotted as functions of ${g}_{\uparrow \downarrow }/g$ in figure 5(a). Both constants are linear functions of ${g}_{\uparrow \downarrow }/g$, which is consistent with the fact that the elastic energy is a linear function of ${g}_{\uparrow \downarrow }/g$ for a fixed vortex-lattice structure (see also figure 4 of [44]). Thus the numerical results are consistent with the effective field theory in the case of overlapping triangular lattices.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have performed similar analyses for interlaced lattices as shown in figures 4(b)–(e). The functions $\{{f}_{i}^{{\rm{P}}/\mathrm{AP}}(\theta )\}$ displayed in the figure show anisotropy except in the right panels for (b) interlaced triangular and (d) square lattices. These behaviors are consistent with the analytical results in equations (40) and (47). Indeed, we can fit the numerical data perfectly using equation (47) if we determine ${C}_{i}^{{\rm{P}}/\mathrm{AP}}\,(i=1,2,3,4)$ and ${D}_{i}^{{\rm{P}}/\mathrm{AP}}\,(i=1,2,3)$ separately for parallel or antiparallel fields. Figure 5(b)–(e) presents the determined constants {C_i} and {D_i}. We note that the constant C₄, which is newly introduced in this work and originates from the coupling between the symmetric and antisymmetric vortex displacements ${{\bf{u}}}_{\pm }$, is indeed nonvanishing for (b) interlaced triangular lattices.

However, the rescaling relations (48) derived from the effective field theory do not hold in figures 4(b)–(e). It can also be seen in different values of the constants $\{{C}_{i}^{{\rm{P}}/\mathrm{AP}}\}$ and $\{{D}_{i}^{{\rm{P}}/\mathrm{AP}}\}$ between the parallel- and antiparallel-field cases in figures 5(b)–(e). The difference between the two cases tends to increase with increasing the ratio ${g}_{\uparrow \downarrow }/g\gt 0$. Furthermore, the constants for (d) square lattices show nonlinear dependences on ${g}_{\uparrow \downarrow }/g$, which is inconsistent with the expected linear dependences for a fixed vortex-lattice structure (see figure 6 of [44]). These results cannot be explained within our effective field theory.

As discussed in the last paragraph of section 3.3, the elastic constants should take the same values between the parallel- and antiparallel-field cases because of the exact correspondence of the GP energy functionals between the two cases [45]. Therefore, a possible insufficiency of our effective field theory may be ascribed to the way the elastic constants are related to the coefficients in the dispersion relations. We infer that the derivative expansions and the coarse graining of the variables done in the derivation of the effective Lagrangian should be improved for interlaced vortex lattices which have a finite displacement between the components.

In the case of parallel fields (${\epsilon }_{\uparrow }={\epsilon }_{\downarrow }=1$), we can drop the subscript α in ${T}_{\alpha }({\bf{s}})$ and ${\tilde{T}}_{\alpha }$. To express the fractional translation, it is useful to modify the basis slightly from the magnetic Bloch states introduced in section 2.2. For the spin-$\downarrow$ component, we use the same magnetic Bloch states as discussed in section 2.2. For the spin-$\uparrow$ component, we define ${{\rm{\Psi }}}_{{\bf{k}}+{\bf{q}},\uparrow }({\bf{r}})$ by operating $\tilde{T}$ on ${{\rm{\Psi }}}_{{\bf{k}}-{\bf{q}},\downarrow }({\bf{r}})$ as

$$\begin{eqnarray}&&\tilde{T}{{\rm{\Psi }}}_{{\bf{k}}-{\bf{q}},\downarrow }({\bf{r}})={{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {{\bf{a}}}_{3}/2}{{\rm{\Psi }}}_{{\bf{k}}+{\bf{q}},\uparrow }({\bf{r}}).\end{eqnarray} \tag{ C.3 }$$

Using $T({{\bf{a}}}_{j})\tilde{T}={{\rm{e}}}^{-[{\bf{K}}\cdot {{\bf{a}}}_{j},{\bf{K}}\cdot {{\bf{a}}}_{3}/2]/{{\hslash }}^{2}}\tilde{T}T({{\bf{a}}}_{j})=-\tilde{T}T({{\bf{a}}}_{j})\,(j=1,2)$, we can confirm that ${{\rm{\Psi }}}_{{\bf{k}}+{\bf{q}},\uparrow }({\bf{r}})$ defined in this way has the expected momentum:

$$\begin{eqnarray*}&&T({{\bf{a}}}_{j}){{\rm{\Psi }}}_{{\bf{k}}+{\bf{q}},\uparrow }({\bf{r}})=-{{\rm{e}}}^{-{\rm{i}}({\bf{k}}-{\bf{q}})\cdot {{\bf{a}}}_{j}}{{\rm{\Psi }}}_{{\bf{k}}+{\bf{q}},\uparrow }({\bf{r}})={{\rm{e}}}^{-{\rm{i}}({\bf{k}}+{\bf{q}})\cdot {{\bf{a}}}_{j}}{{\rm{\Psi }}}_{{\bf{k}}+{\bf{q}},\uparrow }({\bf{r}}).\end{eqnarray*}$$

Furthermore, by operating $\tilde{T}$ on ${{\rm{\Psi }}}_{{\bf{k}}+{\bf{q}},\uparrow }({\bf{r}})$, we have

$$\begin{eqnarray}&&\tilde{T}{{\rm{\Psi }}}_{{\bf{k}}+{\bf{q}},\uparrow }({\bf{r}})={{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {{\bf{a}}}_{3}/2}{\tilde{T}}^{2}{{\rm{\Psi }}}_{{\bf{k}}-{\bf{q}},\downarrow }({\bf{r}})={{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {{\bf{a}}}_{3}/2}{{\rm{\Psi }}}_{{\bf{k}}-{\bf{q}},\downarrow }({\bf{r}}).\end{eqnarray} \tag{ C.4 }$$

Equations (C.3) and (C.4) indicate that the operator $\tilde{T}$ has the role of interchanging ${{\rm{\Psi }}}_{{\bf{k}}-{\bf{q}},\downarrow }({\bf{r}})$ and ${{\rm{\Psi }}}_{{\bf{k}}+{\bf{q}},\uparrow }({\bf{r}})$ with the multiplication of the same phase factor ${{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {{\bf{a}}}_{3}/2}$, which is a useful feature of the present basis. In this representation, one can show

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{V}}_{\alpha \beta }({{\bf{k}}}_{1},{{\bf{k}}}_{2},{{\bf{k}}}_{3},{{\bf{k}}}_{4}) & = & \displaystyle \int {\rm{d}}{\bf{r}}{\rm{d}}{\bf{r}}^{\prime} \,{{\rm{\Psi }}}_{{{\bf{k}}}_{1}+{{\bf{q}}}_{\alpha },\alpha }^{* }({\bf{r}}){{\rm{\Psi }}}_{{{\bf{k}}}_{2}+{{\bf{q}}}_{\beta },\beta }^{* }({\bf{r}}^{\prime} ){g}_{\alpha \beta }{\delta }^{(2)}({\bf{r}}-{\bf{r}}^{\prime} ){{\rm{\Psi }}}_{{{\bf{k}}}_{3}+{{\bf{q}}}_{\beta },\beta }({\bf{r}}^{\prime} ){{\rm{\Psi }}}_{{{\bf{k}}}_{4}+{{\bf{q}}}_{\alpha },\alpha }({\bf{r}})\\ & = & {{\rm{e}}}^{-{\rm{i}}({{\bf{k}}}_{1}+{{\bf{k}}}_{2}-{{\bf{k}}}_{3}-{{\bf{k}}}_{4})\cdot {{\bf{a}}}_{3}/2}\displaystyle \int {\rm{d}}{\bf{r}}{\rm{d}}{\bf{r}}^{\prime} {[\tilde{T}{{\rm{\Psi }}}_{{{\bf{k}}}_{1}+{{\bf{q}}}_{\bar{\alpha }},\bar{\alpha }}({\bf{r}})]}^{* }{[\tilde{T}{{\rm{\Psi }}}_{{{\bf{k}}}_{2}+{{\bf{q}}}_{\bar{\beta }},\bar{\beta }}({\bf{r}}^{\prime} )]}^{* }\\ & & \,\times {g}_{\alpha \beta }{\delta }^{(2)}({\bf{r}}-{\bf{r}}^{\prime} )[\tilde{T}{{\rm{\Psi }}}_{{{\bf{k}}}_{3}+{{\bf{q}}}_{\bar{\beta }},\bar{\beta }}({\bf{r}}^{\prime} )][\tilde{T}{{\rm{\Psi }}}_{{{\bf{k}}}_{4}+{{\bf{q}}}_{\bar{\alpha }},\bar{\alpha }}({\bf{r}})]\\ & = & {{\rm{e}}}^{-{\rm{i}}({{\bf{k}}}_{1}+{{\bf{k}}}_{2}-{{\bf{k}}}_{3}-{{\bf{k}}}_{4})\cdot {{\bf{a}}}_{3}/2}{\tilde{V}}_{\bar{\alpha }\bar{\beta }}({{\bf{k}}}_{1},{{\bf{k}}}_{2},{{\bf{k}}}_{3},{{\bf{k}}}_{4}),\end{array}\end{eqnarray} \tag{ C.5 }$$

where the bars on α and β indicate the spin reversal $\uparrow \,\leftrightarrow \,\downarrow$ and we use the invariance of the interaction ${g}_{\alpha \beta }{\delta }^{(2)}({\bf{r}}-{\bf{r}}^{\prime} )$ under the translation and the spin reversal (${g}_{\alpha \beta }={g}_{\bar{\alpha }\bar{\beta }}$).

For a single particle, we define the fractional translation as the wave function changes by $\tilde{T}$ in equations (C.3) and (C.4) followed by the spin reversal ${\sigma }_{x}$. For many particles, the fractional translation operator ${{ \mathcal T }}^{({\rm{P}})}$ can be expressed in the second-quantized form as

$$\begin{eqnarray}{{ \mathcal T }}^{({\rm{P}})}({\tilde{b}}_{{\bf{k}}\uparrow }^{\dagger },{\tilde{b}}_{{\bf{k}}\downarrow }^{\dagger },{\tilde{b}}_{-{\bf{k}},\uparrow },{\tilde{b}}_{-{\bf{k}},\downarrow }){{ \mathcal T }}^{({\rm{P}})\dagger }={{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {{\bf{a}}}_{3}/2}({\tilde{b}}_{{\bf{k}}\uparrow }^{\dagger },{\tilde{b}}_{{\bf{k}}\downarrow }^{\dagger },{\tilde{b}}_{-{\bf{k}},\uparrow },{\tilde{b}}_{-{\bf{k}},\downarrow })\left(\begin{array}{cc}{\sigma }_{x} & 0\\ 0 & {\sigma }_{x}\end{array}\right).\end{eqnarray} \tag{ C.6 }$$

Using equation (C.5), one can confirm that the Bogoliubov Hamiltonian (18) is invariant under ${{ \mathcal T }}^{({\rm{P}})}$. The ground state $| \mathrm{GS}\rangle$ is obtained as the vacuum annihilated by the Bogolon annihilation operators ${\gamma }_{{\bf{k}},j}\,(j=1,2)$ in equation (21). The single-particle excitations ${\gamma }_{{\bf{k}},j}^{\dagger }| \mathrm{GS}\rangle \,(j=1,2)$ can be used for the Bloch states $| {w}_{{\bf{k}}}^{\pm }\rangle$ in the argument of section 2.5.

In the case of antiparallel fields (${\epsilon }_{\uparrow }=-{\epsilon }_{\downarrow }=1$), we again modify the basis slightly from the magnetic Bloch states introduced in section 2.2. While we use the same magnetic Bloch states as in section 2.2 for the spin-$\downarrow$ component, we define ${{\rm{\Psi }}}_{{\bf{k}}+{\bf{q}},\uparrow }({\bf{r}})$ for the spin-$\uparrow$ component via

$$\begin{eqnarray}&&{\tilde{T}}_{\uparrow }{{\rm{\Psi }}}_{-{\bf{k}}+{\bf{q}},\downarrow }^{* }({\bf{r}})={{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {{\bf{a}}}_{3}/2}{{\rm{\Psi }}}_{{\bf{k}}+{\bf{q}},\uparrow }({\bf{r}}).\end{eqnarray} \tag{ C.7 }$$

Using ${T}_{\alpha }({{\bf{a}}}_{j}){\tilde{T}}_{\alpha }=-{\tilde{T}}_{\alpha }{T}_{\alpha }({{\bf{a}}}_{j})$ and ${T}_{\alpha }^{* }({{\bf{a}}}_{j})={T}_{\bar{\alpha }}({{\bf{a}}}_{j})\,(j=1,2;\alpha =\uparrow ,\downarrow )$, we can confirm that ${{\rm{\Psi }}}_{{\bf{k}}+{\bf{q}},\uparrow }({\bf{r}})$ defined in this way has the expected momentum:

$$\begin{eqnarray*}&&{T}_{\uparrow }({{\bf{a}}}_{j}){{\rm{\Psi }}}_{{\bf{k}}+{\bf{q}},\uparrow }({\bf{r}})=-{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {{\bf{a}}}_{3}/2}{\tilde{T}}_{\uparrow }{[{T}_{\downarrow }({{\bf{a}}}_{j}){{\rm{\Psi }}}_{-{\bf{k}}+{\bf{q}},\downarrow }({\bf{r}})]}^{* }={{\rm{e}}}^{-{\rm{i}}({\bf{k}}+{\bf{q}})\cdot {{\bf{a}}}_{j}}{{\rm{\Psi }}}_{{\bf{k}}+{\bf{q}},\uparrow }({\bf{r}}).\end{eqnarray*}$$

We also find

$$\begin{eqnarray}&&{\tilde{T}}_{\downarrow }{{\rm{\Psi }}}_{-{\bf{k}}+{\bf{q}},\uparrow }^{* }({\bf{r}})={\tilde{T}}_{\downarrow }{[{{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {{\bf{a}}}_{3}/2}{\tilde{T}}_{\uparrow }{{\rm{\Psi }}}_{{\bf{k}}+{\bf{q}},\downarrow }^{* }({\bf{r}})]}^{* }={{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {{\bf{a}}}_{3}/2}{{\rm{\Psi }}}_{{\bf{k}}+{\bf{q}},\downarrow }({\bf{r}}).\end{eqnarray} \tag{ C.8 }$$

In this representation, one can show

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{V}}_{\alpha \beta }({{\bf{k}}}_{1},{{\bf{k}}}_{2},{{\bf{k}}}_{3},{{\bf{k}}}_{4}) & = & {{\rm{e}}}^{-{\rm{i}}({{\bf{k}}}_{1}+{{\bf{k}}}_{2}-{{\bf{k}}}_{3}-{{\bf{k}}}_{4})\cdot {{\bf{a}}}_{3}/2}\displaystyle \int {\rm{d}}{\bf{r}}{\rm{d}}{\bf{r}}^{\prime} {[{\tilde{T}}_{\alpha }{{\rm{\Psi }}}_{-{{\bf{k}}}_{1}+{\bf{q}},\bar{\alpha }}^{* }({\bf{r}})]}^{* }{[{\tilde{T}}_{\beta }{{\rm{\Psi }}}_{-{{\bf{k}}}_{2}+{\bf{q}},\bar{\beta }}^{* }({\bf{r}}^{\prime} )]}^{* }\\ & & \,\times {g}_{\alpha \beta }{\delta }^{(2)}({\bf{r}}-{\bf{r}}^{\prime} )[{\tilde{T}}_{\beta }{{\rm{\Psi }}}_{-{{\bf{k}}}_{3}+{\bf{q}},\bar{\beta }}^{* }({\bf{r}}^{\prime} )][{\tilde{T}}_{\alpha }{{\rm{\Psi }}}_{-{{\bf{k}}}_{4}+{\bf{q}},\bar{\alpha }}^{* }({\bf{r}})]\\ & = & {{\rm{e}}}^{-{\rm{i}}({{\bf{k}}}_{1}+{{\bf{k}}}_{2}-{{\bf{k}}}_{3}-{{\bf{k}}}_{4})\cdot {{\bf{a}}}_{3}/2}{\tilde{V}}_{\bar{\alpha }\bar{\beta }}^{* }(-{{\bf{k}}}_{1},-{{\bf{k}}}_{2},-{{\bf{k}}}_{3},-{{\bf{k}}}_{4}).\end{array}\end{eqnarray} \tag{ C.9 }$$

We define the fractional translation as the time reversal followed by the translation by $\tilde{T}$. Here, the time reversal involves the complex conjugation, the wave vector reversal ${\bf{k}}\to -{\bf{k}}$ (about ${\bf{q}}$), and the spin reversal $\uparrow \,\leftrightarrow \,\downarrow$. In the second-quantized form, the fractional translation operator ${{ \mathcal T }}^{(\mathrm{AP})}$ for many particles is represented as

$$\begin{eqnarray}{{ \mathcal T }}^{(\mathrm{AP})}({\tilde{b}}_{-{\bf{k}},\uparrow }^{\dagger },{\tilde{b}}_{-{\bf{k}},\downarrow }^{\dagger },{\tilde{b}}_{{\bf{k}}\uparrow },{\tilde{b}}_{{\bf{k}}\downarrow }){{ \mathcal T }}^{(\mathrm{AP})\dagger }={{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {{\bf{a}}}_{3}/2}({\tilde{b}}_{{\bf{k}}\uparrow }^{\dagger },{\tilde{b}}_{{\bf{k}}\downarrow }^{\dagger },{\tilde{b}}_{-{\bf{k}},\uparrow },{\tilde{b}}_{-{\bf{k}},\downarrow })\left(\begin{array}{cc}{\sigma }_{x} & 0\\ 0 & {\sigma }_{x}\end{array}\right).\end{eqnarray} \tag{ C.10 }$$

Since ${{ \mathcal T }}^{(\mathrm{AP})}$ is antiunitary, we find

$$\begin{eqnarray*}&&{({{ \mathcal T }}^{(\mathrm{AP})})}^{2}({\tilde{b}}_{{\bf{k}}\uparrow }^{\dagger },{\tilde{b}}_{{\bf{k}}\downarrow }^{\dagger },{\tilde{b}}_{-{\bf{k}},\uparrow },{\tilde{b}}_{-{\bf{k}},\downarrow }){({{ \mathcal T }}^{(\mathrm{AP})\dagger })}^{2}={{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {{\bf{a}}}_{3}}({\tilde{b}}_{{\bf{k}}\uparrow }^{\dagger },{\tilde{b}}_{{\bf{k}}\downarrow }^{\dagger },{\tilde{b}}_{-{\bf{k}},\uparrow },{\tilde{b}}_{-{\bf{k}},\downarrow }),\end{eqnarray*}$$

by which we can confirm that ${({{ \mathcal T }}^{(\mathrm{AP})})}^{2}$ is indeed equal to the translation by ${{\bf{a}}}_{3}$. By using equation (C.9), we can also confirm that the Bogoliubov Hamiltonian (18) is invariant under ${{ \mathcal T }}^{(\mathrm{AP})}$.

Finally, we note that in the above argument, we have used ${\sigma }_{x}$ rather than the more standard one ${\rm{i}}{\sigma }_{y}$ for the spin part of the time reversal. If we define ${\tilde{{ \mathcal T }}}^{(\mathrm{AP})}$ by replacing ${\sigma }_{x}$ by ${\rm{i}}{\sigma }_{y}$ in equation (C.10), the original Hamiltonian (10) in the LLL basis is invariant under ${\tilde{{ \mathcal T }}}^{(\mathrm{AP})}$. However, the Bogoliubov Hamiltonian (18) obtained after the breaking of U(1) × U(1) symmetry as in equation (17) is not invariant under ${\tilde{{ \mathcal T }}}^{(\mathrm{AP})}$ because of the presence of the terms ${\tilde{b}}_{{\bf{k}}\alpha }^{\dagger }{\tilde{b}}_{-{\bf{k}},\alpha }^{\dagger }$ and ${\tilde{b}}_{-{\bf{k}},\alpha }{\tilde{b}}_{{\bf{k}}\alpha }$. Namely, the mixing of a particle and a hole in the Bogoliubov theory is in conflict with time-reversal symmetry in the standard form (see [64] for a different type of conflict between condensation and time-reversal symmetry).