$\mathbb{Z}_{2}$ phases and Majorana spectroscopy in paired Bose-Hubbard chains

While $p$-wave fermionic pairing is a subject of intense ongoing investigations, analogous bosonic phenomena are less well explored. Here we pose and address a variety of questions regarding pairing in spinless 1D bosonic systems: What new phases emerge due to pairing? What are the effects of interactions? What parallels exist to 1D fermionic counterparts? The fermionic Kitaev chain hosts gapped topological phases and mid-gap Majorana states; do bosonic systems exhibit related features?

The Bose-Hubbard model offers fertile ground to explore these questions; its rich phase diagram has been characterized for decades Fisher et al. (1989); Sachdev (2011); Jaksch et al. (1998); Kuklov et al. (2004); Sheshadri et al. (1993), sustaining a productive interplay between theory and experiment since the advent of direct realizations using ultracold atoms in optical lattices Jaksch et al. (1998); Greiner et al. (2002). The phase diagram hosts a gapless superfluid phase exhibiting number fluctuations and condensation as well as multiple gapped Mott phases. A commonly employed mapping between site occupations and spins serves a fruitful source of physical intuition and a connection to models of magnetism. The description reduces to a spin-$1/2$ system in the limit of strong interactions where the occupation number on each site can only fluctuate between particular numbers $n_{0}$ and $n_{0}+1$ Sachdev (2011); Sheshadri et al. (1993); Barankov et al. (2007); Matsubara and Matsuda (1956); Bruder et al. (1993); Altman and Auerbach (2002); Rokhsar and Kotliar (1991). Here, we focus on a homogeneous 1D chain of bosons hopping between neighboring lattice sites with tunneling strength $w$ and an on-site interaction $U$. The chemical potential $\mu$ is determined by the total number of bosons. The additional crucial ingredient is a pairing contribution that creates and annihilates pairs of bosons on neighboring sites with strength $\Delta$. This chain is described by the Hamiltonian

Here $\langle ij\rangle$ denotes a sum over nearest neighbor sites, $b_{i}^{\dagger}$ ($b_{i}$) is the creation (annihilation) operator for a boson on site $i$, and $\hat{n_{i}}=b_{i}^{\dagger}b_{i}$ is the number operator. The effect of the pairing term, $\Delta$, on 1D Bose-Hubbard physics is the main subject of this work.

Combining established insights stemming from boson-to-spin and spin-to-fermion maps in one-dimension, we demonstrate that pairing gives rise to remarkable new features in strongly interacting Bose-Hubbard chains. Most prominently, as depicted in Fig. 1, we find that pairing reduces the $U(1)$ symmetry of the unpaired system to an Ising symmetry, leading to a unique $\mathbb{Z}_{2}$ symmetry-broken phase. As with the usual condensate, which corresponds to the $U(1)$ symmetry-broken phase, this phase exhibits number fluctuations on each site. However, it is gapped and lacks the Goldstone modes associated with an ordinary superfluid. Moreover, its ground state is doubly degenerate, reflecting two very different superpositions of bosonic occupation on every site. Effectively, this subspace forms a gap-protected macroscopic qubit. We note a close relation to non-Abelian Majorana bound state physics: while this $\mathbb{Z}_{2}$ phase is not a topological phase that hosts such states, the two phases can be connected via a specific non-local mapping. Furthermore, since the energy spectrum of the Hamiltonian ${\cal H}$ is identical to that of the fermionic Kitaev chain, a Bose-Hubbard realization would enable experimental measurement of the spectrum for the entire phase diagram. In what follows, we present our analyses of these features, beginning for concreteness by briefly describing a possible realization.

Experimental access to the paired Bose-Hubbard model is most naturally attained using controlled systems of photons or bosonic atoms. Such systems have attracted theoretical interest due to the close relationship to the Kitaev chain: proposals and investigations include employing photonic modes as analogues to Majorana states Bardyn and İmamoǧlu (2012); McDonald et al. (2018); Menssen et al. (2020), and ultracold fermions Goldman et al. (2016); Nascimbène (2013); Jiang et al. (2011); Kraus et al. (2013) and bosons Sträter et al. (2016) for the realization of anyonic modes. In the context of cold atoms, experimental access to the bosonic Kitaev chain analogue is more straightforward than in the original fermionic case: an optical lattice configured as a biased zigzag chain offers a possible natural experimental realization of the Hamiltonian ${\cal H}$ in Eq. 1. Fig. 2 diagrams a section of such a potential, which can be created in a four-beam monochromatic optical lattice with a controllable relative phase between the beams. An appropriate choice of phases yields an array of quasi-1D zig-zag ladders having a tunable chemical potential offset $\mu_{0}$ between the legs. We define the higher-energy leg as our system of interest and the lower-energy leg as the “reservoir chain.” In a tight-binding picture, the hopping term $w$ is achieved through tunneling along the upper chain. The pairing term $\Delta$ is realized by setting $\mu_{0}=U_{r}$ where $U_{r}$ is the interaction energy of two bosons on the reservoir chain. At this offset, single-particle tunneling between the legs of the ladder is non-resonant, and resonant pair tunneling dominates. Similar resonant or near-resonant chemical potential offsets have been shown to enable control of superexchange dynamics and realization of spin Hamiltonians in optical lattices Sachdev et al. (2002); Simon et al. (2011); Meinert et al. (2013); Dimitrova et al. (2020). To maintain constant pairing $\Delta$ and an associated phase which we arbitrarily set to 0, the tunneling rate within the reservoir chain should be set at a value sufficient to stabilize superfluidity. Asymmetrical interchain hopping rates can be tuned by adjusting the lattice phase.

Now we show that approaches commonly employed to analyze the standard Bose-Hubbard chain can provide insight into the effects of the additional pairing term. We focus on the regime of $U$ and $\mu$ such that the average density on each site lies between $n_{0}$ and $n_{0}+1$. For larger enough interaction strength, $w/U\ll 1$, the Hilbert space can be restricted to the number-basis states $|n_{0}\rangle$ and $|n_{0}+1\rangle$ at each site. The most energetically relevant excluded states, $|n_{0}-1\rangle$ and $|n_{0}+2\rangle$, would contribute corrections of order of $w^{2}/U$. Following truncation, the bosonic chain can be mapped to an $XY$ (pseudo-)spin-$1/2$ chain in a transverse field Barankov et al. (2007); Sachdev et al. (2002); Rokhsar and Kotliar (1991); Simon et al. (2011); Sachdev (2011). The truncated Hilbert space may be represented by the spin-$1/2$ states, $|n_{0}+1\rangle=|\uparrow\rangle$ and $|n_{0}\rangle=|\downarrow\rangle$, the eigenstates of the operator $s^{z}$ having eigenvalues $\pm 1/2$. The tunneling term in the Bose-Hubbard Hamiltonian can be identified with raising and lowering spin-$1/2$ operators, $s^{+}$ and $s^{-}$, such that $b^{\dagger}_{i}b_{j}\to(n_{0}+1)s_{i}^{+}s_{j}^{-}$. The interaction and the potential energy terms are diagonal in the number basis at each site, and the boson number operator, $\hat{n}_{i}$, can be expressed in terms of the spin-1/2 matrix $s^{z}$: $\hat{n}\to n_{0}+1/2+s^{z}$. In the truncated spin-$1/2$ Hilbert space, the one-dimensional version of the paired Bose-Hubbard Hamiltonian of Eq. 1 takes the form

where we have the identifications $J_{x}\leftrightarrow 2(n_{0}+1)(w+\Delta)$, $J_{y}\leftrightarrow 2(n_{0}+1)(w-\Delta)$, and $h\leftrightarrow\mu-Un_{0}$. As shown in Fig. 1a and 1b, in the standard case without pairing, Mott states having integer boson filling $n_{0}$ and $n_{0}+1$ correspond to gapped phases in which spins are polarized along the $\pm\hat{z}$ directions, while the gapless superfluid phase at intermediate fillings corresponds to a ferromagnetic state which spontaneously breaks $U(1)$ symmetry.

The presence of the pairing term renders the $XY$ couplings anisotropic, indicative of the aforementioned reduction of the gapless $U(1)$ phase to a gapped phase of the Ising $\mathbb{Z}_{2}$ universality class. In order to study this effect of the pairing term, it is useful to recall features of the standard unpaired case, where $\Delta=0$ and $J_{x}=J_{y}\equiv J$. The $U(1)$ symmetry associated with particle number can be seen in Eq. 1 by noting its invariance under the transformation $b^{\dagger}\rightarrow b^{\dagger}e^{i\phi},b\rightarrow be^{-i\phi}$. The superfluid phase associated with the breaking of this symmetry exists even beyond the spin-1/2 limit where $w/U\ll 1$; for large enough tunneling, it includes occupation numbers beyond the $n_{0}$ and $n_{0}+1$ Mott lobes. In the spin-1/2 limit described by Eq. 2, a large enough transverse field strength $h$, or equivalently, deviation of the chemical potential $\mu$ from the value $Un_{0}$, tips the system into the $n_{0}$ or $n_{0}+1$ Mott states, depending on the sign of the field. For field strengths lower than this critical value, the system is in the superfluid state. At a mean-field level, which is more appropriate for higher dimensions but provides good insights into the phases for any dimension, this behavior can be captured by replacing the pseudospin on each site with an expectation value parametrized by $\langle{\bf S}_{i}\rangle=\frac{1}{2}(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$. The Hamiltonian of Eq. 2 then takes the form of a spin on each site in an effective magnetic field, $H_{\mathrm{mf}}=-{\bf S}_{i}\cdot{\bf B}_{\mathrm{mf}}$, where for a chain, ${\bf B}_{\mathrm{mf}}=(J\sin\theta\cos\phi,J\sin\theta\sin\phi,h)$. The equilibrium configuration of the pseudospin corresponds to minimizing the associated energy, yielding $\cos\theta=h/J$ and no constraints on the value of $\phi$. In the Mott phase, the pseudospins are completely polarized along the $z$ direction, i.e. $\langle s^{z}_{i}\rangle=\pm 1/2$, allowing the identification of $\mu_{\pm}=Un_{0}\pm 2w(n_{0}+1)$, the values of the chemical potential at the boundaries of the Mott states having $n_{0}$ and $n_{0}+1$ bosons per site. In the superfluid phase, the mean-field ground state comprises a superposition on each site:

Spontaneous $U(1)$ symmetry breaking entails making a choice of the continuous parameter $\varphi$, giving rise to a gapless Goldstone mode in the excitation spectrum.

Compared with this unpaired Bose-Hubbard chain, in which number fluctuations on a given site are due to nearest neighbor hopping, the pairing term in Eq. 1 allows for the insertion and depletion of bosons in pairs at neighboring sites. The system no longer respects the invariance under $b^{\dagger}\rightarrow b^{\dagger}e^{i\phi},b\rightarrow be^{-i\phi}$ associated with number conservation. The only exception is for the choices $\phi=0,\pi$, reflecting pariwise processes as opposed to those of single bosons. Thus, in the presence of the pairing term, $U(1)$ symmetry is reduced to $\mathbb{Z}_{2}$ Ising symmetry. In the spin-chain limit of Eq. 2, the reduction to the Ising symmetry is reflected in the identification $\Delta\leftrightarrow J_{x}-J_{y}$. The inequality of $J_{x}$ and $J_{y}$ leads to the reduction to Ising symmetry; for $J_{x}>J_{y}$, the ferromagnetic coupling favors ordering along the eigenstates of the Pauli spin $\sigma^{x}$ on each site and for $J_{x}<J_{y}$ along $\sigma^{y}$. In terms of spin symmetry properties, one can see that Eq. 2 is invariant under the transformation $s_{i}^{x}\rightarrow-s_{i}^{x};s_{i}^{y}\rightarrow-s^{i}_{y};s_{i}^{z}\rightarrow s_{i}^{z}$.

Repeating the mean-field argument above in the presence of pairing establishes salient features of the $\mathbb{Z}_{2}$ phase. Energy minimization is no longer independent of $\phi$ due to the presence of the non-vanishing $\Delta$ term. For $J_{x}>J_{y}$, the phase is pinned to values $\phi=0,\pi$, while for $J_{x}<J_{y}$, it is pinned to $\phi=\pi/2,3\pi/2$. In these two cases, the average density, captured by the expectation value $\langle s^{z}_{i}\rangle$, now varies as $\cos\theta=h/J_{x}$ and $\cos\theta=h/J_{y}$, respectively. Hence, once more, the system enters into the regular $n_{0}\>(\theta=\pi/2)$ or $n_{0}+1\>(\theta=0)$ Mott phases for large enough effective field. However, for couplings that allow a solution in the range $0<\theta<\pi/2$ , the symmetry-broken $\mathbb{Z}_{2}$ phase is different from the usual condensate phase. On the one hand, it too exhibits number fluctuations such as those depicted in Eq. 3. On the other, the pairing term pins the phase such that no continuous symmetry is broken. As seen from the mean-field wavefunction, depending on the relative magnitudes of the anisotropic couplings, the phase $\phi$ is pinned to one of two possible values that are consistent with eigenstates of $s^{x}$ or $s^{y}$, respectively, translating to two macroscopically different ground state wavefunctions for the bosonic chain. Furthermore, the absence of $U(1)$ symmetry breaking implies the absence of a low-lying Goldstone mode; the system is gapped.

Thus, a remarkable feature of this bosonic $\mathbb{Z}_{2}$ phase is that its correlated many-body ground state is doubly degenerate and protected by a gap. As a simple example, in the limit of equal hopping and pairing, $\Delta=w$, the system is exactly in the transverse Ising limit $J_{y}=0$. If the system is now set to half filling, $h=\mu-Un_{0}=0$, then the mean-field states described above become the exact ground states with all spins either pointing along $+x$ or $-x$. In terms of bosons, every site has equal superposition of $n_{0}$ and $n_{0}+1$ bosons and the degeneracy corresponds to all sites having a symmetric or antisymmetric superposition.

Having discussed the possible $\mathbb{Z}_{2}$ phases of the paired Bose-Hubbard chain and their ground state properties, we proceed to describe the entire energy spectrum in the regime of interest. To this end, the transverse $XY$ spin chain of Eq. 2 can be diagonalized via the standard Jordan-Wigner transformation Lieb et al. (1961); Sachdev (2011); DeGottardi et al. (2013). Here, the spin operators on a given site $j$ along the chain are expressed in terms of fermionic creation and annihilation operators ($f^{\dagger},f$) through the transformation $s^{x}_{j}=f_{j}^{\dagger}f_{j}-1/2$, $s^{+}_{j}=f^{\dagger}_{j}e^{-i\pi\sum_{l<j}N_{l}}$, where $N_{l}$ is the fermion occupation number on site $l$. The spin chain Hamiltonian, in terms of the fermions, transforms to a one-dimensional tight-binding representation of a p-wave superconductor, known as the Kitaev chain Kitaev (2001); Zhang et al. (2019); Aasen et al. (2016); Sau et al. (2020), whose chief feature is that it supports topologically robust isolated Majorana fermionic bound states. The Kitaev chain consists of spinless fermions experiencing nearest neighbor hopping $\tilde{w}$ and pairing $\tilde{\Delta}$, and is described by the Hamiltonian

In terms of the original Bose-Hubbard parameters, we can identify $\tilde{w}=w(n_{0}+1)$ and $\tilde{\Delta}=-\Delta(n_{0}+1)$. The effective magnetic field experienced by the transverse spin chain translates to $\tilde{\mu}=\mu-Un_{0}$, namely, an effective chemical potential for fermions.

As can be seen by comparing the fermionic Hamiltonian above with the paired bosonic Hamiltonian of Eq.1, nearest neighbor hopping and pairing terms in one system directly translate to those in the other. While the two systems are very similar in form, the nature of the phases and states they exhibit are radically different. The $\mathbb{Z}_{2}$ aspect, while common to both, reflects local order in the bosonic system and topological order in the fermionic system. The bosonic $\mathbb{Z}_{2}$ phase discussed above corresponds to a topological phase in the fermionic chain that supports Majorana bound states at its ends. The ground state degeneracy is common to both systems, and in the fermionic case, corresponds to the pair of Majorana end bound states forming a Dirac fermionic state that can either be occupied or unoccupied. The $\mathbb{Z}_{2}$ degree of freedom is thus associated with electron parity. Spontaneous symmetry breaking in the topological phase thus corresponds to picking odd or even parity, or a specific superposition. The Mott-insulating phase in the bosonic system translates in the fermionic system to a topologically trivial non-degenerate phase having no isolated Majorana bound states.

Not only is the ground state double degeneracy common to the symmetry-broken $\mathbb{Z}_{2}$ phases in the bosonic and fermionic systems described by Eq. 1 and Eq. 4, respectively, as expected of basis-independent properties, the full energy spectrum is identical across the entire range of parameters. Specifically, diagonalizing the quadratic fermionic Hamiltonian by transforming into momentum space gives the energy dispersion

The dispersion is gapped across the range of parameters, except for the lines i) $\tilde{\mu}=2\tilde{w}$, ii) $\tilde{\mu}=-2\tilde{w}$ and iii) $\tilde{\Delta}=0,|\tilde{\mu}|<2\tilde{w}$. Along these lines, the dispersion vanishes at specific momenta; the lines demarcate phase boundaries between different gapped phases. In terms of paired boson physics, the line along $\tilde{\Delta}=0$ corresponds to the standard superfluid-Mott insulator phase diagram in the limit of large interaction for a given hopping strength $w$. As depicted in Fig.1, Phases I and II correspond to the doubly-degenerate $\mathbb{Z}_{2}$ phases. Phases III and IV correspond to the bosonic Mott phases and their topologically trivial fermionic counterparts. For any given point in the phase diagram, Eq. 5 provides the energy dispersion, which can be easily modified to a discrete spectrum for finite-sized systems. The momentum modes, while corresponding to free fermions, are highly non-trivial in terms of bosons due to the non-local mapping between the two bases. The energy spectrum, however, is common to both. Furthermore, in the $\mathbb{Z}_{2}$ phases, the two degenerate states populate the center of the gap. Finite size effects are expected to split this degeneracy.

Returning to experimental realizations, we note that all major ingredients of the proposed scheme have already been demonstrated in cold atom experiments: tunable hopping, interaction-induced resonant tunneling, and of course tight-binding Hamiltonians Greiner et al. (2002); Simon et al. (2011). The standard Bose-Hubbard model has been a particularly fruitful focus for cold atom experimental work: since the initial realization of the superfluid-Mott insulator phase transition with ultracold rubidium in an optical lattice Greiner et al. (2002), a wide variety of experiments have probed aspects ranging from compressibility Gemelke et al. (2009) to string ordering Endres et al. (2011). A cold-atom manifestation of the paired Bose-Hubbard model would build upon these and related developments of advanced characterization techniques such as spectroscopic probes Campbell et al. (2006); Stöferle et al. (2004); Amato-Grill et al. (2019). For revealing the Kitaev chain spectrum, modulation spectroscopy of atoms in optical lattices would be a natural approach. Bragg spectroscopy can reveal energy-quasimomentum relations Clément et al. (2009); Ernst et al. (2010), complemented with new techniques based on position-space Bloch oscillations Geiger et al. (2018); Fujiwara et al. (2019) and exotic drives Rajagopal et al. (2019). In principle, wavefunction features of the $\mathbb{Z}_{2}$ phase highlighted in this work could be observed in interferometric signatures in time-of-flight images of a multiple-chain sample.

In summary, we have investigated the physics of Bose-Hubbard chains with nearest-neighbor pairing, motivated by a possible experimental realization using an interaction anti-blockade in zig-zag optical lattices. Physical insights into the resulting features are revealed by mapping to an anisotropic XY spin chain and from there to a fermionic chain. The model exhibits unusual gapped $\mathbb{Z}_{2}$ phases and a spectrum identical to that of the Kitaev chain, including doubly-degenerate zero-energy Majorana bound states. Future work suggested by these results could include considerations of larger spin, higher dimensions, and inhomogeneities, as well as the exploration of possible connections to $XY$ lattice models of fractionalization in gapless $U(1)$ versus gapped $\mathbb{Z}_{N}$ phases Gorantla et al. (2020).