Charge excitations across a superconductor-insulator transition

The theory of conventional superconductors describes the emergence of pairing among repulsive charges by means of a Fermi surface instability in a metal (or Fermi liquid) [1]. It leads to an effective attractive interaction whose hallmark is the formation of a gap for one-particle excitations, while zero-momentum pair excitations are gapless [2]. Yet, various cases are known to evade this paradigm, where the parent state does not have a well-defined Fermi surface, such as in the case of insulators. This superconductor-insulator transition (SIT) has been experimentally investigated in various low-dimensional systems, either via disorder tuning [3, 4, 5, 6], introducing magnetic fields [7, 8, 9], or introducing changes in the carrier density [10, 11, 12]. The development of analog quantum simulators, comprising ultracold atoms trapped in optical lattices [13, 14, 15], has opened a door to controllably study such SIT transitions in exquisitely clean strongly correlated systems.

Motivated by understanding the nature of SITs in clean strongly correlated systems, as well as their possible experimental exploration in experiments with ultracold fermionic atoms, in a previous work [16] two of us (R.M. and M.R., in collaboration with P. Nikolić) studied the SIT in the attractive Hubbard model in the presence of a staggered potential in the square lattice geometry. We showed that in the insulating phase in the strongly interacting regime, the lowest-energy charge excitations can be bosonic or fermionic, depending on the parameters chosen. We also showed that, for all the Hamiltonian parameters that could be studied using quantum Monte Carlo simulations, the SIT belongs to 3$d$-XY universality class. We concluded from those results that the lowest-energy charge excitations are bosonic in both the superconducting and insulating phases across the SIT, as argued using field-theory arguments in the weak-coupling limit [17, 18]. Such a bosonic-like insulator is akin to a pseudogap phase, where pre-formation of pairs occurs but such pairs are not phase correlated to exhibit superconductivity [16]. Lowest-energy bosonic charge excitations across an SIT were later argued to also occur in a model of two coupled triangular lattices [19, 20].

Here we study the SIT in the ground state of the attractive Hubbard model in the presence of a staggered potential in the honeycomb lattice geometry. Our main goal is to directly probe the nature of the lowest-energy charge excitations across the SIT, which we find to be bosonic, and study their crossover to being fermionic upon increasing the strength of the staggered potential in the insulating phase. For this, we compute the one-particle (fermionic) and two-particle (bosonic) gaps, which are more precise probes of the nature of the lowest energy excitations than the probes used in Ref. [16]. We also study the universality class of the SIT, which we find to be 3$d$-XY as in the square lattice case [16]. Our second goal, which together with potential experimental realizations motivated us to study the honeycomb lattice geometry, is to show that our conclusions are robust independently of the lattice geometry (Ref. [16] focused on the square-lattice geometry); and likely to be universal for SITs in clean systems. The fermionic model considered in this paper, except for the next-nearest-neighbor terms, was studied experimentally with ultracold atoms [21]. To properly account for the effects of quantum fluctuations and strong correlations in our model, we use two unbiased computational techniques, exact diagonalization and quantum Monte Carlo simulations.

The presentation is organized as follows. In Sec. 2, we introduce the attractive Hubbard model in the presence of a staggered potential in the honeycomb lattice geometry, present its phase diagram, and discuss limiting regimes that can be solved analytically. We also briefly introduce the unbiased computational techniques used to study this model. In Sec. 3, we discuss how the phase diagram of the model is obtained using exact diagonalization and quantum Monte Carlo simulations. This is where we show that the SIT belongs to 3$d$-XY universality class. Section 4 is devoted to studying the one-particle (fermionic) and two-particle (bosonic) charge excitations across the SIT and how the nature of the lowest-energy one changes from bosonic to fermionic upon increasing the strength of the staggered potential in the insulating phase. Our results are summarized in Sec. 5.

Our model of interest is the SU(2) honeycomb Hubbard model in the presence of a staggered potential [21],

where $\hat{c}^{\dagger}_{\boldsymbol{i}\sigma}$ ($\hat{c}_{\boldsymbol{i}\sigma}$) are the fermionic creation (annihilation) operators at site $\boldsymbol{i}$, with (pseudo-)spin $\sigma=\uparrow,\downarrow$, and $\hat{n}_{\boldsymbol{i}\sigma}=\hat{c}^{\dagger}_{\boldsymbol{i}\sigma}\hat{c}_{\boldsymbol{i}\sigma}$ is the corresponding (pseudo-)spin site-occupation operator. The nearest-neighbor, $\langle{\boldsymbol{i\,j}}\rangle$, and next-nearest neighbor (NNN), $\langle\langle{\boldsymbol{i\,j}}\rangle\rangle$, hopping parameters are denoted by $t$ and $t^{\prime}$, respectively; the strength of the on-site attractive interaction by $U<0$, and the strength of the staggered potential by $\Delta$ [$s_{\boldsymbol{i}}$ is either 0 or 1 depending on the sublattice to which ${\boldsymbol{i}}$ belongs to, see Fig. 1(a)]. In what follows, we focus on an average filling of one electron per site $N\equiv N_{\uparrow}+N_{\downarrow}=V$, with $N_{\uparrow}\equiv\sum_{\boldsymbol{i}}\langle\hat{n}_{\boldsymbol{i}\uparrow}\rangle=N_{\downarrow}\equiv\sum_{\boldsymbol{i}}\langle\hat{n}_{\boldsymbol{i}\downarrow}\rangle$, and $V=2L^{2}$ being the number of lattice sites [$L$ is the number of unit cells in each direction, see Fig. 1(a)]. We set $t=1$ as the energy scale, but write $t$ whenever it is helpful in equations.

Compared to the model studied experimentally in Ref. [21], Eq. (2) contains only an extra set of terms, namely, the NNN hoppings, which describe hoppings between sites that belong to the same sublattice. As will become apparent below, those terms introduce nontrivial changes to the solely nearest-neighbors model, and they can be introduced in experiments via a modulated shaking of the optical lattice [22].

To understand the interplay of the different terms in Eq. (2), it is useful to analyze several of its limiting regimes. First, let us discuss the noninteracting ($U=0$) limit, in which the Hamiltonian is diagonalizable in $\boldsymbol{k}$ space resulting in the following two bands:

with $f({\boldsymbol{k}})=2\cos(\sqrt{3}k_{y})+4\cos\left(\frac{3}{2}k_{x}\right)\cos\left(\frac{\sqrt{3}}{2}k_{y}\right)$. The first question one can ask is how large $\Delta$ needs to be for the ground state at half filling to transition between the semi-metal ground state for $\Delta=0$ to the band insulator for $\Delta\gg t$. The answer to that question depends on the NNN hopping amplitude $t^{\prime}$. For $t^{\prime}=0$, the Dirac cones of the upper and lower bands, touching at the high symmetry points K and K^′ of the first Brillouin zone, split, generating a direct gap ($=2\Delta$) for any nonzero value of $\Delta$. After including a nonzero $t^{\prime}$, the direct gap is still obtained so long as $t^{\prime}\leq t^{\prime}_{c}=t/3$ [see Fig. 1(b)]. As $t^{\prime}$ increases crossing $t^{\prime}_{c}$, the minimum of the upper band moves from the K,K^′ points to the $\Gamma$ point. For $t^{\prime}>t^{\prime}_{c}$, the gap becomes indirect (K,K${}^{\prime}\leftrightarrow\Gamma$) and the amplitude of the staggered potential $\Delta_{c}$ needed to turn the system into a band-insulator becomes non-vanishing [see Fig. 1(c)]. Specifically, $\Delta_{c}=\frac{9t^{\prime 2}-t^{2}}{2t^{\prime}}$, as depicted by the dashed-dotted line in the $U=0$ plane in Fig. 2.

A second important limit that can be promptly understood is the $t^{\prime}=0$ case for $U\neq 0$. In this regime, see Ref. [16] for a similar discussion in the context of the square lattice, a particle-hole transformation on only one of the components of the (pseudo-)spin, say, the $\downarrow$-component, $\hat{c}_{{\boldsymbol{i}},\downarrow}\leftrightarrow(-1)^{s_{\boldsymbol{i}}}\hat{c}_{{\boldsymbol{i}},\downarrow}^{\dagger}$ ($s_{\boldsymbol{i}}$ is either 0 or 1, depending on the sublattice to which ${\boldsymbol{i}}$ belongs), maps the original attractive Hubbard model onto the repulsive Hubbard model in the presence of a staggered Zeeman field, i.e., $\Delta\sum_{i}(-1)^{s_{\boldsymbol{i}}}(\hat{n}_{\boldsymbol{i}\uparrow}+\hat{n}_{\boldsymbol{i}\downarrow})\leftrightarrow\Delta\sum_{i}(-1)^{s_{\boldsymbol{i}}}(\hat{n}_{\boldsymbol{i}\uparrow}-\hat{n}_{\boldsymbol{i}\downarrow})$. Since a nonzero staggered field ($\Delta\neq 0$) in the repulsive Hubbard model in the honeycomb lattice explicitly breaks SU(2) symmetry, the antiferromagnetic Mott insulator for $U_{c}\gtrsim 3.8$ [23, 24, 25, 26, 27, 28] becomes an $S_{z}$ antiferromagnet. Recalling the mapping between magnetic and charge/pair degrees of freedom in the repulsive and attractive Hubbard models under the particle-hole transformation [29],

one finds that any nonvanishing value of the staggered potential $\Delta$ in the attractive Hubbard model breaks the supersolid state at $|U|>|U_{c}|$, leading to an insulating ground state with a different $\langle\hat{n}_{\boldsymbol{i}}\rangle$ in the two sublattices that make the bipartite honeycomb lattice. An important point to keep in mind is that this difference in the site occupations in the sublattices results from having the staggered potential. It is not the result of a spontaneous symmetry-breaking process, i.e., the insulating ground state is a Mott insulator that does not break symmetries of the model. For $|U|<|U_{c}|$, as for $U=0$ and $|U|>|U_{c}|$, we expect that any non-vanishing value of the staggered potential $\Delta$ leads to an insulating ground state.

For $t^{\prime}>t^{\prime}_{c}$ and $\Delta=0$, we mentioned before that the ground state at $U=0$ is a metal. This means that adding attractive on-site interactions $U<0$ results in a superconducting state. Hence, in the regime with $t^{\prime}>t^{\prime}_{c}$ and $U<0$, a nonzero value $\Delta_{c}$ of the staggered potential is needed to drive the SIT. In the absence of unbiased analytical techniques to tackle this regime for arbitrary values of $U$, $\Delta$, and $t^{\prime}$, we carry out numerical calculations to obtain the phase diagram of Hamiltonian (2) as characterized by the surface $\Delta_{c}(U,t^{\prime})$. A compilation of the results obtained using Krylov-based exact diagonalization [30, 31] in a lattice with $V=2L^{2}=18$ sites is shown in Fig. 2(a). This lattice, which has $L=3$ and is depicted in Fig. 1(a), contains the Brillouin zone corner as a valid momentum point, i.e., it captures the effects of low-energy excitations about the high-symmetry K points. This makes it optimal to determine the phase diagram within the lattice sizes that we can study using exact diagonalization ¹¹1See, e.g., Ref. [65] for a discussion of finite-size effects in topological transitions in commensurate vs. incommensurate lattices.. Qualitatively similar results were obtained on a 16-sites lattice; see Appendix A.

Much larger lattices can be studied by means of projective quantum Monte Carlo (PQMC) simulations. Within this approach, the ground state $|\Psi_{0}\rangle$ is obtained projecting a trial wave function $|\Psi_{T}\rangle$, via $|\Psi_{0}\rangle=\lim_{\Theta\to\infty}e^{-\Theta\hat{H}}|\Psi_{T}\rangle$, where $\Theta$ is a projector parameter. This approach works provided that the overlap between the trial wave function and the ground-state of Hamiltonian (2) is nonzero, $\langle\Psi_{0}|\Psi_{T}\rangle\neq 0$, and that $|\Psi_{0}\rangle$ is nondegenerate [33]. The expectation values of operators $\hat{O}$ in the ground state can then be written as

For our simulations, we considered two trial wave functions, the half-filled Fermi sea of the noninteracting part of the Hamiltonian and a Hartree-Fock solution. We found the latter to converge more quickly with increasing $\Theta$ so all results reported are obtained with that trial wave function for $\Theta t=40$, which is sufficiently large for the convergence of the expectation values of all observables studied here. Since our PQMC calculations are carried out within the canonical ensemble for $N_{\uparrow}=N_{\downarrow}$, they do not suffer from the sign problem [34, 33, 35, 36]. The phase diagram $\Delta_{c}(U,t^{\prime})$ extracted from the PQMC simulations of lattices with up to $L=15$ is reported in Fig. 2(b). It is similar to the one obtained using exact diagonalization in small lattices.

In the next section, we discuss in detail how the phase diagrams reported in Fig. 2 are obtained using exact diagonalization and PQMC simulations.

Before explaining in detail how we locate the phase transition using exact diagonalization and PQMC simulations and how we probe its universality class, let us briefly comment on the overall structure of the phase diagrams in Fig. 2. Increasing attractive interactions results in smaller Cooper pairs: In the limit $|U|\to\infty$ they fit in a site, becoming hardcore bosons [37]. Consequently, the magnitude of the staggered potential needed to transition between the superfluid and insulating ground states decreases as the increasingly local fermionic pairs at stronger interactions become pinned at the $-\Delta$ sites in the lattice. On the other hand, the NNN hopping terms act to counteract the pinning promoted by the staggered potential, so increasing the magnitude of $t^{\prime}$ results in the need for a stronger staggered potential to induce the superfluid-insulator transition. Notably, the phase diagram also indicates that the transition at nonvanishing interactions is smoothly connected to the noninteracting metal-band insulator one, marked by the dash-dotted lines at $U=0$ in both panels of Fig. 2. We also note that, for values of $t^{\prime}>t/3$, the noninteracting regime at $\Delta=0$ no longer corresponds to a semi-metal but rather to a metal with a finite density of states at the Fermi level. Hence, the inclusion of attractive interactions results in pairing and superconductivity.

The 3$d$-XY universality class of the transition advances that the lowest energy charge excitations across the transition are bosonic [16]. With increasing $\Delta$ in the insulating phase, the lowest-energy charge excitations must cross over from bosonic to fermionic (at $\Delta=\infty$ they are fermionic for any finite $U$).

In this section, we study the nature of the lowest-energy charge excitations across the SIT and their behavior with increasing $\Delta$ in the insulating phase. For this, we compute the one-particle ($m=1$, fermionic) and two-particle ($m=2$, bosonic) gaps [60, 61],

where $E_{0}(x)$ is the ground-state energy with $x$-fermions. For $m=2$, we always add/remove one fermion with (pseudo-)spin $\uparrow$ and one with (pseudo-)spin $\downarrow$, namely, a pair. Hence, in what follows, we call the $m=2$ gap the pair gap.

The one-particle and pair gaps in Eq. (10) can be straightforwardly computed using exact diagonalization. Similarly, since for $m=2$ we have that $N_{\uparrow}=N_{\downarrow}$ so there is no sign problem, the pair gap in Eq. (10) can also be computed using PQMC simulations in much larger system sizes. PQMC runs into the sign problem for $m=1$ in Eq. (10), because $N_{\uparrow}\neq N_{\downarrow}$. Hence, we follow a different approach to compute the one-particle gap. We probe the decay of the appropriate time-displaced correlation function [24]. Specifically, the one-particle gap is extracted using the imaginary-time $\tau$ displaced Green’s functions $G_{+}({\bf k},\tau)=\langle\hat{c}_{{\bf k},\sigma}(\tau)\hat{c}^{\dagger}_{{\bf k},\sigma}(0)\rangle$ and $G_{-}({\bf k},\tau)=\langle\hat{c}^{\dagger}_{{\bf k},\sigma}(\tau)\hat{c}_{{\bf k},\sigma}(0)\rangle$, with $\hat{c}^{\dagger}_{{\bf k},\sigma}(\tau)\equiv e^{\tau\hat{H}}c^{\dagger}_{{\bf k},\sigma}(0)e^{-\tau\hat{H}}$, which describe the particle and hole excitations with respect to the Fermi energy, respectively. At large $\tau$’s, they decay as $G_{\pm}({\bf k},\tau)\propto e^{-\tau\delta_{\pm}^{(1)}({\bf k})}$. By comparing the two branches over different momenta ${\bf k}$, the one-particle gap is obtained as $\delta^{(1)}=\min_{\bf k}[\delta_{+}^{(1)}({\bf k})]+\min_{\bf k}[\delta_{-}^{(1)}({\bf k})]$ (see Appendix C for further details about this analysis).

Focusing on $t^{\prime}=0.6$, Figs. 6(a)–6(c) display the one-particle and pair gaps ($m=1$ and $2$) obtained using PQMC in lattices with $L=6$, 9, and 12 (main panels) and using exact diagonalization in a lattice with $L=3$ (insets). The PQMC and exact diagonalization results are qualitatively similar. The pair gap vanishes in the superconducting phase. It then becomes nonzero, with a magnitude that increases with increasing $\Delta$, once the SIT (marked by the dashed line) is crossed. On the other hand, the one-particle gap is nonvanishing in both the superconducting and insulating phases and exhibits a minimum at the SIT. Hence, as advanced, the pair gap is smaller than the one-particle gap in both the superconducting and insulating phases across the SIT. With increasing $\Delta$, those gaps cross in the insulating phase at a value of $\Delta$ that, deep in the strongly interacting regime, increases with increasing $|U|$. The PQMC results for different system sizes show that finite-size effects are small in the gaps computed for those system sizes, and even more so as $|U|$ and $\Delta$ increase. This suggests that our PQMC estimation of the crossing points does not suffer from significant finite-size effects.

Compiling results like those depicted in Fig. 6 allows us to obtain exact diagonalization (inset in Fig. 7) and PQMC (main panel in Fig. 7) phase diagrams in the ($U,\Delta$) plane for $t^{\prime}=0.6$. They highlight that, with increasing $|U|$, there is an increase in the extension of the region in which the lowest-energy charge excitation in the insulating phase is bosonic. The observed weak non-monotonicity of the crossover curve as the interaction strength decreases reflects the bosonic character of the lowest-energy charge excitations close to the SIT.

We used unbiased numerical calculations to study the SIT in the ground state of the attractive honeycomb Hubbard model in the presence of a staggered potential. We directly showed that the lowest-energy charge excitations are bosonic across the transition, and cross over to fermionic in the insulating phase. The former is consistent with the finding that the SIT in this model belongs to the 3$d$-XY universality class. The results obtained in the honeycomb lattice are qualitatively similar to those in the square lattice [16], suggesting that our findings are universal for SITs in clean systems. For example, we expect similar results for the attractive Hubbard model in the triangular lattice if one adds a positive (negative) local potential to one site in each triangle when the filling is $n=4/3$ ($n=2/3$).

Given that, in the absence of nearest-neighbor hoppings, the Hamiltonian considered here has already been simulated in optical lattice experiments [21], we expect that our findings can readily be tested in such experiments. An interesting open question for both theory and experiments is what happens in three dimensions, in which the critical temperature that triggers the onset of superconductivity is nonzero at half-filling [62, 63]. Exploring the interplay between the temperature and the parameters considered here in two dimensions may help improve our understanding of the role of the preformation of pairs in the finite-temperature realm.