Quantum Gutzwiller approach for the two-component Bose-Hubbard model

V. E. Colussi1*, F. Caleffi2, C. Menotti1 A. Recati1,3

1 INO-CNR BEC Center and Dipartimento di Fisica, Università di Trento, Via
Sommarive 14, 38123 Povo, Trento, Italy

2 International School for Advanced Studies (SISSA), Via Bonomea 265, I-34136
Trieste, Italy

3 Trento Institute for Fundamental Physics and Applications, INFN, Via
Sommarive 14, 38123 Povo, Trento, Italy

*<EMAIL_ADDRESS>
## Abstract

We study the effects of quantum fluctuations in the two-component Bose-Hubbard
model generalizing to mixtures the quantum Gutzwiller approach introduced
recently in [Phys. Rev. Research 2, 033276 (2020)]. As a basis for our study,
we analyze the mean-field ground-state phase diagram and spectrum of
elementary excitations, with particular emphasis on the quantum phase
transitions of the model. Within the quantum critical regimes, we address both
the superfluid transport properties and the linear response dynamics to
density and spin probes of direct experimental relevance. Crucially, we find
that quantum fluctuations have a dramatic effect on the drag between the
superfluid species of the system, particularly in the vicinity of the paired
and antipaired phases absent in the usual one-component Bose-Hubbard model.
Additionally, we analyse the contributions of quantum corrections to the one-
body coherence and density/spin fluctuations from the perspective of the
collective modes of the system, providing results for the few-body
correlations in all the regimes of the phase diagram.

###### Contents

  1. 1 Introduction
  2. 2 Model and Theory
    1. 2.1 Lagrangian formulation within the Gutzwiller ansatz for mixtures
    2. 2.2 Ground state and excitations
      1. 2.2.1 Repulsive interaction ($U_{12}>0$)
      2. 2.2.2 Attractive interaction ($U_{12}<0$)
    3. 2.3 Quantum Gutzwiller theory for mixtures
  3. 3 Linear Response
    1. 3.1 Density and spin response
    2. 3.2 Current response and superfluid components
      1. 3.2.1 Superfluid regime
      2. 3.2.2 Phase transitions for $U_{12}>0$
      3. 3.2.3 Phase transitions for $U_{12}<0$
  4. 4 Correlation Functions
    1. 4.1 Coherence function
    2. 4.2 Density and spin fluctuations
  5. 5 Conclusions and Outlook
  6. A Quantum corrections within the QGW approach
    1. A.1 Evaluation protocol for the observables
    2. A.2 First-order fluctuations
    3. A.3 Second-order fluctuations
  7. B Linear response formalism within the QGW theory
    1. B.1 Linear response formalism
    2. B.2 Density response function
    3. B.3 Current response functions
      1. B.3.1 Current operator
      2. B.3.2 Intraspecies response
      3. B.3.3 Interspecies response
    4. B.4 Average kinetic energy

##  1 Introduction

The precision control available in ultracold atomic gas experiments has opened
up a platform where models of condensed matter physics can be simulated in a
relatively defect-free environment. In particular, ultracold atoms trapped in
optical lattices are described by Hubbard models [1, 2], which are of central
importance to the study of solid-state materials [3]. The reduction of three-
body losses in lattices makes possible the study of strongly-correlated
regimes absent in the continuum, such as the paradigmatic Mott insulator (MI)
to superfluid (SF) quantum phase transition for single-component bosons [4].
In the two-component Bose-Hubbard (BH) model a richer phase diagram emerges,
including the additional possibility of pair (PSF) and counterflow (CFSF)
superfluids, supersolidity, charge-density quasiorder, and peculiar magnetic
states [5, 6, 7, 8, 9, 10]. Such coupled superfluids can undergo also mutual
dissipationless transport with an induced entrainment or counterflow of one
component due to a non-zero superfluid velocity of the other. This phenomenon,
better known as superfluid drag, was first discussed by Andreev and Bashkin in
the context of three-fluid hydrodynamics [11], but is of universal relevance
to systems ranging from neutron-star matter [12, 13, 14, 15, 16] to
multicomponent superconductors [17, 18, 19] and ultracold atomic mixtures [20,
21, 22, 23, 24, 25, 26, 27, 28, 9]. Direct measurement of this effect has
however remained elusive, due in part to the low miscibility of superfluid 3He
and 4He, and recombination heating in strongly-interacting ultracold atomic
mixtures. Recently, the PSF and CFSF phase transitions of the two-component BH
model have emerged as promising candidates where the drag can saturate at its
maximum value [25, 26]. Still, a deeper understanding of the fundamental role
played by quantum fluctuations is needed to gain insight into the physics of
such strongly-correlated quantum critical regimes of Hubbard models at zero
temperature.

The Gutzwiller approach provides a simple solution of the BH model able to
interpolate between a coherent-state description of a dilute Bose-Einstein
condensate (BEC) and a MI with integer occupation of each lattice site [29,
30, 31]. Its generalization to the two-component BH model is based on the wave
function ansatz [32]

$|\Psi_{G}\rangle=\bigotimes_{\bf
r}\sum_{n_{1},n_{2}=0}^{\infty}c_{n_{1},n_{2}}({\bf r})\,|n_{1},n_{2},{\bf
r}\rangle\,,$ (1)

which is a site-factorized product of local Fock states with occupation of
$n_{i}$ atoms in the $i^{\text{th}}$ component at site ${\bf r}$ weighted by
normalized complex amplitudes $c_{n_{1},n_{2}}({\bf r})$, which can capture at
the mean-field level paired and antipaired ground states and superfluidity
[32]. The Gutzwiller description is however inadequate in strongly-interacting
and critical regimes, where quantum fluctuations are enhanced. Of course, also
the basic assumptions of the celebrated Bogoliubov theory of quantum
fluctuations about a coherent-state BEC break down in the vicinity of the
phase transitions [33, 34]. In this spirit, a canonical quantization procedure
$c\to\hat{c}$ was proposed in Ref. [35], which enables an order-by-order
layering of quantum fluctuations on top of the mean-field Gutzwiller ground
state, remaining valid throughout the phase diagram. Predictions of the
resultant quantum Gutzwiller theory (QGW) featured a successful comparison
across the MI-SF phase transition for single-component bosons to the available
quantum Monte Carlo (QMC) data, while requiring only modest computational
costs and providing a valuable description of the role of quantum fluctuations
through a semi-analytical formalism. Moreover, the QGW theory was used
recently to study the non-Markovian spin decoherence of a two-level impurity
embedded in a two-dimensional BH model [36]. It is important to note that, at
the level of Gaussian fluctuations, the QGW approach is analogous to slave
boson methods (c.f Refs. [37, 38, 39]). However, in the two-component BH model
comparatively little is known about time-dependent phenomena and the role of
quantum fluctuations in quantum critical regimes, motivating the establishment
of the QGW theory for bosonic mixtures.

In this work, we study a homogeneous system of two-component bosons on a
square lattice considering short-range intra and interspecies interactions,
which can be realized in optical lattices loaded with atoms of two different
species or internal states [40, 41, 42]. Although the derivations presented in
this work are completely general, our numerical findings are specific to two
dimensions where existing QMC results [25] make quantitative comparisons with
our predictions for the drag possible, and where there exists a strong
motivating analogy between the fermionic version of the problem and high-
temperature superconductivity [43]. We first study the rich phase diagram of
the model for both repulsive and attractive interspecies interactions, finding
counterflow and paired superfluid phases in addition to the MI and SF phases
which carry over from single-component bosonic systems. We show how also in
the case of mixtures, the QGW approach provides a straightforward way to
calculate linear response and correlation functions to a desired order in the
quantum fluctuations. This permits a systematic study of the role of quantum
corrections, which we first investigate by considering the linear response
dynamics of the system to density and spin perturbations. In this respect, we
highlight the experimentally relevant signatures of the onset of the PSF and
CFSF phases in the dynamical structure factor. We then focus on superfluid
transport in quantum critical regimes, finding in particular a large
interspecies drag comparable in magnitude with the superfluid density in the
vicinity of the CFSF and PSF phases. Furthermore, we address the one and two-
body correlation functions, focusing on the strongly-interacting regime, where
quantum fluctuations play a crucial role. Specifically, we find that the PSF
and CFSF transitions behave as channel-selective MI transitions with respect
to the spin and density degrees of freedom, respectively.

We have structured the present work in a relatively self-contained manner,
which we hope may prove useful also as a starting point and reference for
future studies. In Sec. 2, we first outline the time-dependent Gutzwiller
method for the two-component BH model, and then study the ground state and
elementary excitations of the system for both repulsive and attractive
interactions, expanding on Ref. [32] where possible. Subsequently, to include
the effects of quantum fluctuations, we introduce the QGW approach, which
generalizes the formalism of Ref. [35] to mixtures. The remainder of the paper
is focused on studying the role of quantum fluctuations in linear response
dynamics and superfluidity, addressed in Sec. 3, and few-body correlations,
discussed in Sec. 4. In these sections, general analytic formulas are derived
for the susceptibility, compressibility, density and spin sound speeds,
superfluid components (including the drag), one-body coherence function, and
spin and density pair correlation functions. In particular, the order-by-order
contributions of quantum fluctuations are made always explicitly clear. We
conclude in Sec. 5 and provide also an outlook for future studies. The
appendices contain additional details on the QGW approach and derivations of
the linear response formalism for mixtures.

##  2 Model and Theory

In Sec. 2.1, we analyze the two-component BH model using the
$\mathbb{C}$-number Gutzwiller ansatz. Subsequently, in Sec. 2.2.1 (2.2.2), we
discuss the ground state and elementary excitations for repulsive (attractive)
interactions between the two components of the system. In Sec. 2.3, we go
beyond the mean-field Gutzwiller approximation by generalizing the QGW method
introduced in Ref. [35] to Bose mixtures.

###  2.1 Lagrangian formulation within the Gutzwiller ansatz for mixtures

We start from the two-component Bose-Hubbard (BH) model [1]

$\hat{H}=\sum_{i=1}^{2}\sum_{\bf
r}\left[-J_{i}\sum_{j=1}^{d}\left(\hat{a}^{\dagger}_{i,{\bf
r}}\,\hat{a}_{i,{\bf r}+{\bf
e}_{j}}+\mathrm{h.c.}\right)+\frac{U_{i}}{2}\,\hat{n}_{i,{\bf
r}}\left(\hat{n}_{i,{\bf r}}-1\right)-\mu_{i}\,\hat{n}_{i,{\bf
r}}\right]+U_{12}\sum_{\bf r}\hat{n}_{1,{\bf r}}\,\hat{n}_{2,{\bf r}},$ (2)

where ${\bf r}$ is the site index in $d$ dimensions with unit vector ${\bf
e}_{j}$, while $J_{i}$, $\mu_{i}$, $U_{i}$ and $U_{12}$ are the species-
dependent tunneling coefficient between neighboring sites, chemical potential
and on-site intra/interspecies interaction strengths, respectively. In this
work, we examine only the $\mathbb{Z}_{2}$-symmetric case where
$J_{1}=J_{2}=J$, $\mu_{i}=\mu_{2}=\mu$ and $U_{1}=U_{2}=U$, keeping always
$\left|U_{12}/U\right|<1$ to avoid phase separation, as well as to prevent the
system from collapsing [44]. The bosonic creation and annihilation operators
$\hat{a}^{\dagger}_{i,{\bf r}}$ and $\hat{a}_{i,{\bf r}}$ create and destroy,
respectively, a particle of species $i$ at the lattice site ${\bf r}$, and are
related to the corresponding local density operator via $\hat{n}_{i,{\bf
r}}=\hat{a}^{\dagger}_{i,{\bf r}}\,\hat{a}_{i,{\bf r}}$. We define also the
total and spin densities as given by $\hat{n}_{d,{\bf
r}}=\sum_{i}\hat{n}_{i,{\bf r}}$ and $\hat{n}_{s,{\bf r}}=\hat{n}_{1,{\bf
r}}-\hat{n}_{2,{\bf r}}$, respectively. In the following, we consider a
uniform cubic lattice of volume $V$ composed by $I$ sites with lattice spacing
$a$111In this work, the QGW results are obtained on a lattice of size
$I=128^{2}$ unless otherwise specified. This lattice is sufficiently large
such that our calculations are free from finite-size effects., such that the
“bare” effective mass is $m\equiv\hbar^{2}/\left(2Ja^{2}\right)$.

The two-component Gutzwiller mean-field ansatz for the many-body wave function
introduced in Eq. (1) is a site-factorized product of local Fock states
weighted by normalized complex amplitudes $c_{n_{1},n_{2}}({\bf r})$ such that
$\sum_{n_{1},n_{2}}\left|c_{n_{1},n_{2}}({\bf r})\right|^{2}=1$. Within the
Gutzwiller ansatz, the local density is given by

$n_{i}({\bf r})=\langle\hat{n}_{i,{\bf
r}}\rangle=\sum_{n_{1},n_{2}}^{\infty}\left(n_{1}\,\delta_{i,1}+n_{2}\,\delta_{i,2}\right)\left|c_{n_{1},n_{2}}({\bf
r})\right|^{2}\,.$ (3)

In addition to the Mott insulator (MI) and superfluid (SF) phases in common
with the one-component case (cf. Ref. [29]), the two-component Bose mixture
exhibits also the possibility of counterflow superfluid (CFSF) and pair
superfluid (PSF) phases, see Refs. [7, 6, 10, 45]. To distinguish between
these phases, besides the one-body order parameters

$\psi_{1}({\bf r})=\langle\hat{a}_{1,{\bf
r}}\rangle=\sum_{n_{1},n_{2}}^{\infty}\sqrt{n_{1}}\,c^{*}_{n_{1}-1,n_{2}}({\bf
r})\,c_{n_{1},n_{2}}({\bf r})\,,$ (4a) $\psi_{2}({\bf
r})=\langle\hat{a}_{2,{\bf
r}}\rangle=\sum_{n_{1},n_{2}}^{\infty}\sqrt{n_{2}}\,c^{*}_{n_{1},n_{2}-1}({\bf
r})\,c_{n_{1},n_{2}}({\bf r})\,,$ (4b)

which are non-zero only in the SF phase, we introduce the pair and antipair
order parameters

$\psi_{\mathrm{P}}({\bf r})=\langle\hat{a}_{1,{\bf r}}\,\hat{a}_{2,{\bf
r}}\rangle-\langle\hat{a}_{1,{\bf r}}\rangle\,\langle\hat{a}_{2,{\bf
r}}\rangle=\sum_{n_{1},n_{2}}\sqrt{n_{1}\,n_{2}}\,c^{*}_{n_{1}-1,n_{2}-1}({\bf
r})\,c_{n_{1},n_{2}}({\bf r})-\psi_{1}({\bf r})\,\psi_{2}({\bf r})\,,$ (5a)
$\psi_{\mathrm{C}}({\bf r})=\langle\hat{a}_{1,{\bf
r}}\,\hat{a}^{\dagger}_{2,{\bf r}}\rangle-\langle\hat{a}_{1,{\bf
r}}\rangle\,\langle\hat{a}^{\dagger}_{2,{\bf
r}}\rangle=\sum_{n_{1},n_{2}}\sqrt{n_{1}(n_{2}+1)}\,c^{*}_{n_{1}-1,n_{2}+1}({\bf
r})\,c_{n_{1},n_{2}}({\bf r})-\psi_{1}({\bf r})\,\psi^{*}_{2}({\bf r})\,,$
(5b)

which identify univocally the PSF and CFSF phases, respectively. Notice that,
in constructing the pair/antipair order parameters, possible contributions of
the one-body order parameters have been explicitly removed to ensure that
$\psi_{\mathrm{P/C}}\neq 0$ reflects intrinsically particle + hole (CFSF) or
particle + particle (hole + hole) (PSF) correlations between the two species.

From the Hamiltonian expectation value
$\langle\Psi_{G}|\,\hat{H}\,|\Psi_{G}\rangle$, we can readily build a
Lagrangian for the Gutzwiller ansatz,

$\displaystyle\mathcal{L}[c,c^{*}]=$ $\displaystyle\,\sum_{\bf
r}\sum_{n_{1},n_{2}}\left\\{\frac{i\,\hbar}{2}\left[c^{*}_{n_{1},n_{2}}({\bf
r})\,\dot{c}_{n_{1},n_{2}}({\bf
r})-\mathrm{c.c.}\right]-H_{n_{1},n_{2}}\left|c_{n_{1},n_{2}}({\bf
r})\right|^{2}\right\\}$ $\displaystyle+J\sum_{i=1}^{2}\sum_{\bf
r}\sum_{j=1}^{d}\left[\psi^{*}_{i}({\bf r})\,\psi_{i}({\bf r}+{\bf
e}_{j})+\mathrm{c.c.}\right]\,,$ (6)

where the dots indicate temporal derivatives and

$H_{n_{1},n_{2}}=\sum_{i=1}^{2}\left[\frac{U}{2}\,n_{i}(n_{i}-1)-\mu\,n_{i}\right]+U_{12}\,n_{1}\,n_{2}\,.$
(7)

The classical Euler-Lagrange equations of motion for the Gutzwiller
amplitudes, with the complex conjugate parameters $c^{*}_{n_{1},n_{2}}({\bf
r})=\partial\mathcal{L}/\partial\dot{c}_{n_{1},n_{2}}({\bf r})$ playing the
role of canonical momenta, are given by the two-component time-dependent
Gutzwiller equations (2GE)

$\displaystyle i\,\hbar\,\dot{c}_{n_{1},n_{2}}({\bf
r})=H_{n_{1},n_{2}}\,c_{n_{1},n_{2}}-J\sum_{i=1}^{d}\bigg{\\{}\sqrt{n_{1}+1}\,c_{n_{1}+1,n_{2}}\left[\psi^{*}_{1}({\bf
r}+{\bf e}_{i})+\psi^{*}_{1}({\bf r}-{\bf e}_{i})\right]$
$\displaystyle+\sqrt{n_{1}}\,c_{n_{1}-1,n_{2}}({\bf r})\left[\psi_{1}({\bf
r}+{\bf e}_{i})+\psi_{1}({\bf r}-{\bf
e}_{i})\right]+\sqrt{n_{2}+1}\,c_{n_{1},n_{2}+1}({\bf
r})\left[\psi^{*}_{2}({\bf r}+{\bf e}_{i})+\psi^{*}_{2}({\bf r}-{\bf
e}_{i})\right]$ $\displaystyle+\sqrt{n_{2}}\,c_{n_{1},n_{2}-1}({\bf
r})\left[\psi_{2}({\bf r}+{\bf e}_{i})+\psi_{2}({\bf r}-{\bf
e}_{i})\right]\bigg{\\}}\,,$ (8)

which were derived previously in Ref. [32]. The 2GE are straightforward
generalizations of their one-component counterparts (cf. Ref. [29]) with the
additional contribution of the diagonal interspecies coupling $U_{12}$. In
order to explore solutions of the 2GE, we search first for the stationary
solutions $c_{n_{1},n_{2}}({\bf
r})=c^{0}_{n_{1},n_{2}}\,e^{-i\,\omega_{0}\,t}$ independent of the site index
${\bf r}$. The ground state energy is then given by

$\hbar\,\omega_{0}=-4\,d\,J\sum_{i=1}^{2}\left|\psi_{0,i}\right|^{2}+\sum_{n_{1},n_{2}}\left\\{\sum_{i=1}^{2}\left[\frac{U}{2}\,n_{i}(n_{i}-1)-\mu\,n_{i}\right]+U_{12}\,n_{1}\,n_{2}\right\\}\left|c^{0}_{n_{1},n_{2}}\right|^{2}\,,$
(9)

where the “$0$” sub/superscript indicates quantities evaluated with respect to
the ground state Gutzwiller amplitudes $c^{0}_{n_{1},n_{2}}$ obtained by
diagonalizing Eq. (2.1). For instance, the expression of the mean-field total
density is simply given by

$n_{0,d}=\sum_{n_{1},n_{2}}^{\infty}\left(n_{1}+n_{2}\right)\left|c^{0}_{n_{1},n_{2}}({\bf
r})\right|^{2}\,.$ (10)

In the remainder, it is assumed always that the ground state parameters
$c^{0}_{n_{1},n_{2}}$ are real numbers.

In order to study the linear-response dynamics of the system around the ground
state, we consider small perturbations around the stationary solution of the
form

$c_{n_{1},n_{2}}({\bf r})=\left[c^{0}_{n_{1},n_{2}}+c^{1}_{n_{1},n_{2}}({\bf
r},t)\right]e^{-i\,\omega_{0}\,t},$ (11)

which can be expanded in terms of plane waves as

$c^{1}_{n_{1},n_{2}}({\bf r},t)=\sum_{\bf k}\left[u_{{\bf
k},n_{1},n_{2}}\,e^{i({\bf k\cdot r}-\omega_{\bf k}\,t)}+v^{*}_{{\bf
k},n_{1},n_{2}}\,e^{-i({\bf k\cdot r}-\omega_{\bf k}\,t)}\right]\,.$ (12)

Linearizing the 2GE with respect to the amplitudes $u_{\bf k,n_{1},n_{2}}$ and
$v_{\bf k,n_{1},n_{2}}$, one obtains an eigenvalue problem

$\hbar\,\omega_{\bf k}\begin{pmatrix}\underline{u}_{\bf k}\\\
\underline{v}_{\bf k}\end{pmatrix}=\mathcal{L}_{\bf
k}\begin{pmatrix}\underline{u}_{\bf k}\\\ \underline{v}_{\bf
k}\end{pmatrix}\quad,\quad\hat{\mathcal{L}}_{\bf k}=\begin{pmatrix}H_{\bf
k}&K_{\bf k}\\\ -K_{\bf k}&-H_{\bf k}\end{pmatrix},$ (13)

where the vectors $\underline{u}_{\alpha,\mathbf{k}}$
$\left(\underline{v}_{\alpha,\mathbf{k}}\right)$ contain the particle (hole)
components $u_{\alpha,\mathbf{k},n_{1},n_{2}}$
$\left(v_{\alpha,\mathbf{k},n_{1},n_{2}}\right)$, respectively. The positive
eigenvalues $\omega_{\bf k}$ of the pseudo-Hermitian matrix
$\hat{\mathcal{L}}_{\bf k}$ describe the multi-branch excitation spectrum of
the fluctuations and are identified with the collective modes of the system.
The elements of the matrix blocks $H_{\bf k}$ and $K_{\bf k}$ have expressions

$\displaystyle H^{n_{1},n_{2},n^{\prime}_{1},n^{\prime}_{2}}_{\bf k}=$
$\displaystyle\left\\{\sum_{i=1}^{2}\left[\frac{U}{2}\,n_{i}(n_{i}-1)-\mu\,n_{i}\right]+U_{12}\,n_{1}\,n_{2}-\hbar\,\omega_{0}\right]\delta_{n_{1},n_{1}^{\prime}}\delta_{n_{2},n_{2}^{\prime}}$
(14a) $\displaystyle-J({\bf
0})\,\psi_{0,1}\left(\sqrt{n^{\prime}_{1}}\,\delta_{n^{\prime}_{1},n_{1}+1}+\sqrt{n_{1}}\,\delta_{n_{1},n^{\prime}_{1}+1}\right)\delta_{n_{2},n^{\prime}_{2}}$
$\displaystyle-J({\bf
0})\,\psi_{0,2}\left(\sqrt{n^{\prime}_{2}}\,\delta_{n^{\prime}_{2},n_{2}+1}+\sqrt{n_{2}}\,\delta_{n_{2},n^{\prime}_{2}+1}\right)\delta_{n_{1},n^{\prime}_{1}}$
$\displaystyle-J({\bf
k})\left[\sqrt{n_{1}+1}\,\sqrt{n^{\prime}_{1}+1}\,c^{0}_{n_{1}+1,n_{2}}\,c^{0}_{n^{\prime}_{1}+1,n^{\prime}_{2}}+\sqrt{n_{1}}\,\sqrt{n^{\prime}_{1}}\,c^{0}_{n_{1}-1,n_{2}}\,c^{0}_{n^{\prime}_{1}-1,n^{\prime}_{2}}\right]$
$\displaystyle-J({\bf
k})\left[\sqrt{n_{2}+1}\,\sqrt{n^{\prime}_{2}+1}\,c^{0}_{n_{1},n_{2}+1}\,c^{0}_{n^{\prime}_{1},n^{\prime}_{2}+1}+\sqrt{n_{2}}\,\sqrt{n^{\prime}_{2}}\,c^{0}_{n_{1},n_{2}-1}\,c^{0}_{n^{\prime}_{1},n^{\prime}_{2}-1}\right]\,,$
$\displaystyle K^{n_{1},n_{2},n^{\prime}_{1},n^{\prime}_{2}}_{\bf k}=$
$\displaystyle-J({\bf
k})\left[\sqrt{n_{1}+1}\,\sqrt{n^{\prime}_{1}}\,c^{0}_{n_{1}+1,n_{2}}\,c^{0}_{n^{\prime}_{1}-1,n^{\prime}_{2}}+\sqrt{n_{1}}\,\sqrt{n^{\prime}_{1}+1}\,c^{0}_{n_{1}-1,n_{2}}\,c^{0}_{n^{\prime}_{1}+1,n^{\prime}_{2}}\right]$
(14b) $\displaystyle-J({\bf
k})\left[\sqrt{n_{2}+1}\,\sqrt{n^{\prime}_{2}}\,c^{0}_{n_{1},n_{2}+1}\,c^{0}_{n^{\prime}_{1},n^{\prime}_{2}-1}+\sqrt{n_{2}}\,\sqrt{n^{\prime}_{2}+1}\,c^{0}_{n_{1},n_{2}-1}\,c^{0}_{n^{\prime}_{1},n^{\prime}_{2}+1}\right]\,,$

with $J({\bf k})=2\,d\,J-\epsilon({\bf k})$ and where

$\epsilon({\bf
k})=4\,J\sum_{j=1}^{d}\sin^{2}{\left(\frac{k_{j}\,a}{2}\right)}$ (15)

is the free-particle dispersion law on a $d$-dimensional lattice. It follows
that the dependence on $\mathbf{k}$ of the excitations is determined by the
variable

$x=\left[\frac{1}{d}\sum_{j=1}^{d}\sin^{2}{\left(\frac{k_{j}\,a}{2}\right)}\right]^{1/2}\,,$
(16)

which varies from 0 to 1 and scales as
$x\approx\left|\mathbf{k}\right|/\left(2\,d\right)$ at small momenta.

Consequently, the linear response dynamics of the one-body order parameter
$\psi_{i}({\bf r},t)=\psi_{0,i}+\psi^{1}_{i}({\bf r},t)$ can be expanded in
plane waves as

$\psi^{1}_{i}({\bf r},t)=\sum_{\bf k}\left[U_{i,{\bf k}}\,e^{i({\bf
k}\cdot{\bf r}-\omega_{\bf k}\,t)}+V^{*}_{i,{\bf k}}\,e^{-i({\bf k}\cdot{\bf
r}-\omega_{\bf k}\,t)}\right]\,,$ (17)

where we have introduced the particle-hole amplitudes

$U_{1,{\bf
k}}=\sum_{n_{1},n_{2}}\sqrt{n_{1}}\left(c^{0}_{n_{1}-1,n_{2}}\,u_{{\bf
k},n_{1},n_{2}}+c^{0}_{n_{1},n_{2}}\,v_{{\bf k},n_{1}-1,n_{2}}\right)\,,$
(18a) $V_{1,{\bf
k}}=\sum_{n_{1},n_{2}}\sqrt{n_{1}}\left(c^{0}_{n_{1},n_{2}}\,u_{{\bf
k},n_{1}-1,n_{2}}+c^{0}_{n_{1}-1,n_{2}}\,v_{{\bf k},n_{1},n_{2}}\right)\,,$
(18b) $U_{2,{\bf
k}}=\sum_{n_{1},n_{2}}\sqrt{n_{2}}\left(c^{0}_{n_{1},n_{2}-1}\,u_{{\bf
k},n_{1},n_{2}}+c^{0}_{n_{1},n_{2}}\,v_{{\bf k},n_{1},n_{2}-1}\right)\,,$
(18c) $V_{2,{\bf
k}}=\sum_{n_{1},n_{2}}\sqrt{n_{2}}\left(c^{0}_{n_{1},n_{2}}\,u_{{\bf
k},n_{1},n_{2}-1}+c^{0}_{n_{1},n_{2}-1}\,v_{{\bf k},n_{1},n_{2}}\right)\,.$
(18d)

It is possible also to study the linear response dynamics of the pair and
antipair order parameters by defining analogous pair/antipair amplitudes
$U_{\mathrm{P/C},{\bf k}}$ and $V_{\mathrm{P/C},{\bf k}}$, which are discussed
further in App. A.2.

Finally, the linear response of the local density for each species is given by

$n_{i}({\bf r},t)=n_{0,i}+\sum_{\bf k}\left[N_{i,{\bf k}}\,e^{i({\bf
k}\cdot{\bf r}-\omega_{\bf k}\,t)}+\text{c.c.}\right]\,,$ (19)

where we have introduced the perturbation amplitude

$N_{i,{\bf k}}=\sum_{n_{1},n_{2}}n_{i}\,c^{0}_{n_{1},n_{2}}\left(u_{{\bf
k},n_{1},n_{2}}+v_{{\bf k},n_{1},n_{2}}\right)\,.$ (20)

In the following two subsections, we proceed to analyze in detail the
properties of the ground state and excitations within the SF, MI, CFSF and PSF
phases for both a repulsive (Sec. 2.2.1) and an attractive (Sec. 2.2.2)
coupling between the two Bose components.

###  2.2 Ground state and excitations

We explore the phase diagram of the system in the presence of repulsive and
attractive interspecies interactions in Figs. 1 and 4. We preface this
discussion by noting that in contrast to the single-species BH model, one
finds also paired CFSF and PSF phases signalled by a non-zero pairing order
parameter $\psi_{0,\mathrm{C}}$ and $\psi_{0,\mathrm{P}}$, respectively.
Additionally, these order parameters may be nonzero in the vicinity of the
various phase transitions alongside a finite one-body condensate $\psi_{0,i}$,
which signals the entrance into the SF region. Ultimately, the pair/antipair
order parameters vanish in the limit $2\,d\,J/U\gg 1$ as expected.

####  2.2.1 Repulsive interaction ($U_{12}>0$)

Figure 1: Ground state phase diagram for repulsive interspecies interactions
(a) $U_{12}/U=0.5$ and (b) $U_{12}/U=0.9$. First and second-order transition
lines are indicated by dashed and solid lines respectively, separating the MI
(CFSF) lobes, identified by light blue (yellow) areas, from the SF region. The
O$(2)$ tip transitions are indicated by red dots, while the dot-shaped
illustrations depict the particle-hole antipairs coupling the two components
in the CFSF phase. (c) Behavior of the order parameters at the crossing of the
first-order MI-to-SF transition for fixed $n_{0,d}=2$ and $U_{12}/U=0.5$. (d)
Behavior of the order parameters at the crossing of the second-order CFSF-to-
SF transition for fixed $n_{0,d}=1$ and $U_{12}/U=0.9$.

(Mott Insulator) – In Figs. 1(a)-(b), we show the ground state phase diagram
for two values of $U_{12}>0$. In general, for strong enough $U$ we find MI
regions [light blue areas] at even total filling whose ground state is
$|n_{0,d}/2,n_{0,d}/2\rangle$ with energy

$\hbar\,\omega_{0}=U\frac{n_{0,d}}{2}\left(\frac{n_{0,d}}{2}-1\right)-\mu\,n_{0,d}+U_{12}\frac{n^{2}_{0,d}}{4}\,.$
(21)

Within the MI lobes, the excitation spectrum can be calculated analytically
from Eq. (13) as

$\hbar\,\omega_{\pm,{\bf k}}=\frac{1}{2}\left\\{\sqrt{J^{2}({\bf k})-2\,J({\bf
k})\,U(n_{0,d}+1)+U^{2}}\pm\left[J({\bf
k})-U(n_{0,d}-1)-2\,U_{12}+2\,\mu\right]\right\\}\,,$ (22)

which is a modification of the one-component result to include the mean-field
interaction energy between different species. This spectrum describes four
dispersive branches in total, a pair of degenerate particle (‘$+$’) branches
and a pair of degenerate hole (‘$-$’) branches.222On top of the particle-hole
excitations, Eq. (13) presents an infinity of non-zero, uncoupled diagonal
elements which describe non-dispersive bands with energy
$\hbar\,\omega_{n_{1},n_{2}}=\sum_{i=1}^{2}\left[\frac{U}{2}n_{i}(n_{i}-1)-\mu\,n_{i}\right]+U_{12}\,n_{1}\,n_{2}-\omega_{0}\,,$
(23) where the occupation indices $n_{1}$ and $n_{2}$ must be chosen not to
fall into the $4\times 4$ block that yields the dispersive bands (22).
Therefore, Eq. (23) applies to non-negative integers $n_{1}$ and $n_{2}$ such
that neither $(n_{1},n_{2})=(j,j\pm 1)$ nor $(j\pm 1,j)$ are satisfied, with
$j\in\mathbb{N}$. We note also that, as a consequence of the $\mathbb{Z}_{2}$
symmetry, these bands are doubly degenerate
$\omega_{n_{1},n_{2}}=\omega_{n_{2},n_{1}}$ for $n_{1}\neq n_{2}$ and singly
degenerate for $n_{1}=n_{2}$.

The second-order phase transition boundary between the MI and SF phases is
determined by the onset of a finite one-body order parameter $\psi_{0,i}$ and
the disappearance of the gap in the excitation spectrum ($\omega_{\bf
k=0}=0$). From Eq. (22), we find that this occurs for

$2\,d\left(\frac{J}{U}\right)_{c}=\frac{(n_{0,d}/2-\mu/U+U_{12}/U)(\mu/U-U_{12}/U-n_{0,\,d}/2+1)}{1+\mu/U-U_{12}/U}\,,$
(24)

which can be linked to the MI boundary of the one-component case via the
mapping $\mu\to\mu-U_{12}$. The maximal value of $2\,d\,(J/U)_{c}$ determines
the locations of the O$(2)$ tip transitions at

$2\,d\left(\frac{J}{U}\right)^{\mathrm{max}}_{c}=\left(\sqrt{\frac{n_{0,d}}{2}+1}-\sqrt{\frac{n_{0,d}}{2}}\right)^{2}\,.$
(25)

Again, we note that the chemical potential of the tip critical points is
similarly shifted from the one-component value and is given by

$\left(\frac{\mu}{U}\right)_{c}=\sqrt{\frac{n_{0,d}}{2}\left(\frac{n_{0,d}}{2}+1\right)}+\frac{U_{12}}{U}-1\,.$
(26)

Figure 2:  Excitation spectra across the MI-to-SF second-order (a) edge
transition at $\mu/U=1$ and (b) O$(2)$ transition at $(\mu/U)_{c}\approx
0.914$, as well as (c) in the superfluid phase for $\mu/U=(\mu/U)_{c}$,
choosing $U_{12}/U=0.5$ at fixed $n_{0,d}=2$ in $d=2$. For the edge transition
in panel (a), the chosen hopping energies are $2\,d\,J/U=0.12$ (red lines) and
$2\,d\,(J/U)_{c}\approx 0.167$ (blue lines); for the O$(2)$ transition in
panel (b), the hopping energies $2\,d\,J/U=0.12$ (red lines) and
$2\,d\,(J/U)^{\mathrm{max}}_{c}\approx 0.172$ (blue lines) are considered. The
hopping energy in the SF phase in panel (c) is $2\,d\,J/U=0.18$. The
parentheses in panels (a)-(b) refer to the degeneracy of the hybridized
particle and hole bands, while the spin and density Goldstone modes are
indicated explicitly.

In Fig. 2(a), we show how the band structure changes as the second-order MI-
to-SF transition is traversed through the edge of the lobe (namely, away from
the tip), specifically for $n_{0,d}=2$ at $2\,d\,(J/U)_{c}\approx 0.167$ for
$U_{12}/U=0.5$ and $\mu/U=1$. In general, in the MI phase, if the chemical
potential is set above (below) its value at the tip, the first two bands are a
pair of degenerate gapped particle (hole) bands of mixed spin or density
character, whereas the next two bands correspond to degenerate gapped hole
(particle) excitations. The doubly degenerate non-dispersive band
$\omega_{0,2}$ is also visible as a horizontal line. At the transition point,
the gap of the lowest pair of degenerate particle (hole) bands vanish, while
the degenerate hole (particle) bands remain gapped. Additionally, the gapless
modes are purely quadratic at low momenta, signalling the characteristic
vanishing of the speed of sound in common with the one-component case [29]. In
Fig. 2(b), we perform a similar analysis of the spectral features of the
second-order MI-to-SF transition when traversed through the tip critical point
at $2\,d\,(J/U)^{\mathrm{max}}_{c}\approx 0.172$ for $U_{12}/U=0.5$ and
$(\mu/U)_{c}\approx 0.914$. Most importantly, we observe that at the tip
transition point the gaps of the doubly degenerate particle and hole bands
vanish simultaneously. Also, their dispersion becomes degenerate and perfectly
linear at low momenta, which is characteristic of a finite sound speed at the
O$(2)$ transition [29].

Once the SF phase takes over, both at the edge and at the tip transition, the
excitation bands which become gapless hybridize into spin and density modes
(see Fig. 2(c)). These excitations correspond to the Goldstone modes that
result from the breaking of the two $U(1)$ symmetries of the model, one for
the density and the other for the spin channel. Their dispersion relation
approaches the Bogoliubov bands [29]

$\hbar\omega_{d/s,{\bf k}}=\sqrt{\epsilon({\bf k})\left[\epsilon({\bf
k})+n_{0,d}\,(U\pm U_{12})\right]}\,$ (27)

in the weakly-interacting limit $2\,d\,J/U\gg 1$. Close to the transition, the
first two gapped branches, which also display individual density and spin
character, are referred to as the Higgs modes of the system, in analogy to the
single-component BH model [46, 47, 48, 49]. The $\hbar\,\omega_{0,2}$ band
becomes dispersive and hybridizes into spin and density excitations as well.

The transition between MI and SF phase can also be of first-order, as
discussed in Refs. [6, 50, 32, 45]. In this case, the one-body order
parameters $\psi_{0,i}$ exhibit a discontinuity across the critical boundary,
as shown in Fig. 1(d). The behavior of the discontinuity at the first-order
transition was studied in the mean-field Gutzwiller analysis of Ref. [32],
where it was found that: (i) the jump of the order parameter increases with
$U_{12}$ and then rapidly goes to zero as the phase separation boundary
$U_{12}/U\sim 1$ is approached; (ii) the hopping window corresponding to the
first-order transition widens with increasing $U_{12}$ starting from the tip
and reaching $J=0$ when approaching the phase separation condition. Across the
first-order critical point, the structure of the excitation spectrum changes
discontinuously between the different spectral features explored before, such
that the modes in the MI and SF phases cannot be smoothly connected.

(Counterflow Superfluid) – In Figs. 1(a)-(b), the phase diagram displays also
CFSF phases [yellow areas] at odd total filling, characterized by a finite
antipair order parameter $\psi_{0,C}$ and whose size increases with larger
$U_{12}$. Within the CFSF lobes, one finds that the Fock states
$|\,(n_{0,d}+1)/2,(n_{0,d}-1)/2\,\rangle$ and
$|\,(n_{0,d}-1)/2,(n_{0,d}+1)/2\,\rangle$ are doubly degenerate with ground
state energy

$\hbar\,\omega_{0}=U\,(n_{0,d}-1)^{2}+U_{12}\,n_{0,d}\,(n_{0,d}-1)+\mu\,(1-2\,n_{0,d})\,.$
(28)

Comparing this energy with Eq. (21), we find the boundary between neighboring
MI and CFSF lobes at $J=0$, which in the case of the $n_{0,d}=2$ MI lobe and
$n_{0,d}=1$ CFSF lobe occurs at $\mu=U_{12}$. In order to obtain the correct
ground state, we symmetrize the antipair state as
$[\,|\,(n_{0,d}+1)/2,(n_{0,d}-1)/2\,\rangle+|\,(n_{0,d}-1)/2,(n_{0,d}+1)/2\,\rangle\,]/\sqrt{2}$
as predicted in Refs. [7, 46, 6, 51, 10, 9].

In the CFSF phase, single particle and hole excitations involve the subset of
states $\\{|(n_{0,d}+1)/2\pm
1,(n_{0,d}-1)/2\rangle;|(n_{0,d}+1)/2,(n_{0,d}-1)/2\pm 1\rangle\\}$, which
amount to three particle-like excitations and one (three) hole-excitations for
the CFSF phase at $n_{0,d}=1$ ($n_{0,d}\geq 3$). From those excitations, one
can construct three (four) modes belonging to the density channel, plus one
(two) belonging to the spin channel. The general behaviour is such that a pair
of particle and hole excitations in the density channel lower their energy
while moving towards the boundary of the CFSF lobe, with the particle (hole)
excitation closing the gap if the transition in crossed above (below) the tip
chemical potential. Exactly at the tip, the lowest particle and hole
excitations in the density channel close the gap simultaneously. Remarkably,
the distinction between the edge and tip critical points in the density
channel closely resembles the properties of the MI-to-SF transition.

For the CFSF phase at $n_{0,d}=1$, the excitation spectrum can be calculated
analytically from Eq. (13). In that case, the spin channel hosts only one
particle branch, decoupled from the rest and given by

$\hbar\,\omega_{\bf k}=U-\mu-J({\bf k})\,.$ (29)

This corresponds to a free-particle dispersion, shifted by the mean-field
interaction energy $U-\mu$. In the density sector, we extract two particle
branches and one hole branch, whose excitation energies are the solutions of
the equation

$J({\bf k})\left[U\,U_{12}-(\hbar\,\omega_{\bf
k}\pm\mu)^{2}\right]\mp(\hbar\,\omega_{\bf k}\pm\mu)(U\mp\hbar\,\omega_{\bf
k}-\mu)(U_{12}\mp\hbar\,\omega_{\bf k}-\mu)=0\,.$ (30)

Similarly to the spectrum of the MI phase, all the other excitations consist
in an infinite sequence of non-dispersive bands, the first two of which have
energies $\hbar\,\omega_{0,1}=\hbar\,\omega_{1,0}=-\mu-\hbar\,\omega_{0}=0$.
In particular, these branches correspond to cost-free excitations and mirror
the degeneracy of antipair states characterizing the mean-field ground state
of the CFSF phase333We note that the flat bands $\omega_{0,1}=\omega_{1,0}$
describe (trivial) excitations within the the antipair sector
$\\{u_{0,1},u_{1,0},v_{0,1},v_{1,0}\\}$ and cannot describe the excitation and
tunnelling of antipairs out of the ground state. Describing pairing collective
modes within our model leads to violated completeness relations (see App.
A.3). In Ref. [32], these bands were found to acquire a sound-like profile
when higher-order hopping processes are included perturbatively in CFSF phase.
These processes are however neglected in the 2GE (2.1), such that we find
higher gapped excitation modes, absent in that work, whose low-energy behavior
is strongly tied to the appearance of the SF phase and determines the one-body
correlations in the CSFS phase (see Sec. 4)..

The second-order phase transition boundary between the CFSF and SF phases is
determined from the closure of the smallest gap in the CFSF excitation
spectrum. From Eq. (30), we find that for $n_{0,d}=1$ the gap closing occurs
at

$2\,d\left(\frac{J}{U}\right)_{c}=\frac{\mu}{U}\frac{\left(1-\mu/U\right)\left(\mu/U-U_{12}/U\right)}{(\mu/U)^{2}-U_{12}/U}\,,$
(31)

whose maximal value gives the location of the tip of the CFSF lobe. As a
reference, for $U_{12}/U=0.9$, the location of the tip of the $n_{0,d}=1$ CFSF
lobe shown in Fig. 1(b) is located at $\\{(\mu/U)_{c}\approx
0.391,2\,d\,(J/U)^{\mathrm{max}}_{c}\approx 0.164\\}$.

Figure 3:  Excitation spectra across the CFSF-to-SF second-order for
$U_{12}/U=0.9$ at fixed $n_{0,d}=1$ in $d=2$ for: (a) the edge transition at
$\mu/U=0.5$; (b) the O$(2)$ transition at $(\mu/U)_{c}\approx 0.391$; (c) the
superfluid phase for $\mu/U=(\mu/U)_{c}$. For the edge transition in panel
(a), the chosen hopping energies are $2\,d\,J/U=0.1$ (red lines) and
$2\,d\,(J/U)_{c}\approx 0.154$ (blue lines); for the O$(2)$ transition in
panel (b), the hopping energies are $2\,d\,J/U=0.14$ (red lines) and
$2\,d\,(J/U)^{\mathrm{max}}_{c}\approx 0.162$ (blue lines). The hopping energy
in the SF phase in panel (c) is $2\,d\,J/U=0.172$. The particle (hole) and
spin (density) characters of gapped excitations, as well as the physical
nature of the Goldstone modes, are indicated explicitly by the labels “p”
(“h”) and “s” (“d”), respectively.

In Fig. 3(a), we show how the excitation spectrum appears as the second-order
CFSF-to-SF transition for $n_{0,d}=1$ and $U_{12}/U=0.9$ is crossed through
the edge point at $2\,d\,(J/U)_{c}\approx 0.154$ at $\mu/U=0.5$. Above (below)
the tip chemical potential, the gap is closed by the lowest of particle (hole)
band in the density channel, while the other bands remain gapped. More
specifically, the particle band in the spin channel corresponding to Eq. (29),
which in this case at $x=0$ is the third in ascending order, will never
participate in the gap closure. In Fig. 3(b), we consider the evolution of the
band structure across the tip transition at
$2\,d\,(J/U)^{\mathrm{max}}_{c}\approx 0.162$ for $U_{12}/U=0.9$ and
$(\mu/U)_{c}\approx 0.391$, while holding $n_{0,d}=1$ fixed. Indeed, we
observe that at the transition point the gaps of the two lowest-energy bands
in the density channel vanish, while the remaining modes retain a finite gap.
In the SF region, illustrated in Fig. 3(c), the two lowest bands in the
density channel participate in the creation of the density Goldstone and the
density Higgs excitation respectively. The spin Goldstone mode emerges from
the non-dispersive antipair band as described before.

####  2.2.2 Attractive interaction ($U_{12}<0$)

Figure 4: Ground state phase diagram for attractive interspecies interactions
(a) $U_{12}/U=-0.4$ and (b) $U_{12}/U=-0.7$, with the order of the transition
lines and O$(2)$ critical points within the MI lobes indicated as in Fig. 1.
The PSF phase is identified by the purple horizontal lines, with total filling
$n_{0,d}$ evaluated in the $\eta\to 0^{+}$ limit. The dot-shaped illustrations
depict the particle-particle pairs that develop between the two components of
the gas in the PSF phase. (c) Behavior of the order parameters along the PSF-
to-SF transition line at the border between the vacuum (green-shaded area) and
the $n_{0,d}=2$ MI lobes for $U_{12}/U=-0.7$ and fixed $n_{0,d}\approx 1.47$.
(d)-(e) Excitation spectrum and compressibility near the first-order vacuum to
superfluid transition for $\mu/U\approx-1.00956$, $2\,d\,J/U=1$, and
$U_{12}/U=-0.6$ with critical filling $n_{0,d}\approx 0.239$. The log-log
scale of the vertical axis of panel (e) reveals the quadratic power law of the
density Goldstone mode at all momenta.

(Mott Insulator) – In Figs. 4(a)-(b), we show the ground state phase diagram
for two different values of $U_{12}<0$. For attractive interactions, the MI
lobes present the same ground state and spectral properties as their repulsive
counterparts, with the spin/density character of the excitation modes being
reversed. Furthermore, the MI-to-SF criticality is found again to be of either
first or second order. However, at odds with the repulsive case, the first-
order boundaries appear initially at small $J/U$ rather than close to the tip,
as shown in Fig. 4(b) for $U_{12}/U=-0.4$, eventually spanning the entire lobe
boundary, as shown in Fig. 4(b) for $U_{12}/U=-0.7$. In addition, the boundary
between the vacuum lobe [green area] and the superfluid phase becomes also a
first-order transition line. Along this boundary, the density Goldstone mode
acquires a purely quadratic dispersion $\omega_{d,\bf k}\propto{\bf k}^{2}$ at
small momenta, as shown in Fig. 4(d), which indicates the vanishing of the
sound velocity of density excitations ($c_{d}\to 0$), while the spin sound
velocity remains finite ($c_{s}>0$). Accordingly, the (mean-field)
compressibility $\partial n_{0,d}/\partial\mu$ diverges due to the
discontinuity in the filling, as shown in Fig. 4(e). These behaviors indicate
that the system collapses along the first-order vacuum-to-SF transition
boundary. Along this line, the critical filling decreases for increasing
$2\,d\,J/U$, vanishing eventually in the deep SF regime (not shown) where the
transition becomes again of second order. As $U_{12}/U$ becomes more
attractive, the collapse transition line extends towards larger values of
$2\,d\,J/U$.

Figure 5:  Variation of (a) $n_{0,d}$ and $\psi_{0,\mathrm{P}}$ and (b)
$\alpha$ and $\beta$ as functions of $\mu/U$ in the vicinity of the PSF line
at $(\mu/U)_{c}=-0.35$ for fixed $2\,d\,J/U=0.12$ and $U_{12}/U=-0.7$. Solid
(dashed) lines correspond to the symmetry-breaking offset
$\eta=10^{-3}(2\times 10^{-3})$.

(Pair Superfluid) – The MI lobes shown in Figs. 4(a)-(b) are separated by
sharp transition lines [purple horizontal lines], extending towards
increasingly large values of $2\,d\,J/U$ as $U_{12}$ becomes more and more
attractive. On these lines at the boundary of the $n_{0,d}=2\,j$ and
$n_{0,d}=2\left(j+1\right)$ MI lobes, one finds that the states with the
lowest eigenenergies $|j+1,j+1\rangle$ and $|j,j\rangle$ are degenerate, thus
corresponding to a PSF phase with a finite pair coherence $\psi_{0,P}\neq 0$.
When imposed in the 2GE, such degeneracy condition pinpoints the location of
the PSF lines on a discrete set of critical chemical potentials $(\mu/U)_{c}$
through the equation

$-U_{12}-2\,j\,(U+U_{12})+2\,\mu_{c}=0\,,$ (32)

valid for any non-negative integer $j$. Thus, in general the ground state is
given by the superposition $\alpha\,|j+1,j+1\rangle+\beta\,|j,j\rangle$, where
the coefficients $\alpha$ and $\beta$ are restricted to lie on the unit circle
$|\alpha|^{2}+|\beta|^{2}=1$. However, these coefficients are undetermined
within the usual Gutzwiller approximation, as the individual states are
already $\mathbb{Z}_{2}$-symmetric.

In the QMC results of Ref. [25], $\alpha$ and $\beta$ were found to vary with
$\mu$ at fixed $U_{12}$ and $J$. The determination of $\alpha$ and $\beta$ in
the mean-field Gutzwiller theory requires the introduction of an ad hoc
perturbative order parameter $\psi_{0,i}=\eta>0$ along the PSF line mimicking
the tunneling of residual background fluctuations, such that ground states
with broken U$(1)$ symmetry can be accessed. As we show in Fig. 5(a), this
turns out to uniquely fix $\alpha$ [orange lines] and $\beta$ [light blue
lines] along the PSF lines described by Eq. (32) regardless of the value of
$\eta$, provided that $\eta$ remains sufficiently small. In particular, for
the PSF line separating the vacuum region from the $n_{0,d}=2$ MI lobe, we
obtain $\alpha\approx 0.858$ and $\beta\approx 0.513$, from which we derive
$n_{0,d}\approx 1.47$. Similarly, for the PSF line separating the $n_{0,d}=2$
and $n_{0,d}=4$ MI lobes, we find $\alpha\approx 0.794$ and $\beta\approx
0.608$, which gives $n_{0,d}\approx 3.260$ in the PSF phase. In this way, we
reliably obtain a PSF ground state with a well-defined value of the pair order
parameter $\psi_{0,\mathrm{P}}\neq 0$.

The symmetry-breaking offset $\eta\neq 0$ makes the PSF phase line to develop
a finite width, as shown in Fig. 5(b), where the variation of $n_{0,d}$ and
$\psi_{0,\mathrm{P}}$ is studied for fixed $\mu/U$. We find that $n_{0,d}$
[light blue lines] changes continuously between the filling of the two
neighboring MI lobes, which is qualitatively consistent with the QMC results
of Ref. [25], with the major difference being the size of the transition width
in the chemical potential $\mu/U$. Secondly, we also see that
$\psi_{0,\mathrm{P}}$ [orange lines] behaves smoothly, vanishing identically
as one enters the MI lobes and reaching a maximum in between. This reveals
that the PSF-to-MI transitions are of second-order within the present
approach. In general, the maximum of $\psi_{0,\mathrm{P}}$ is found to occur
for $\mu/U$ below the exact PSF line, becoming increasingly shifted for larger
$\eta$. Also, the maximum is located where $\alpha=\beta=1/\sqrt{2}$ [black
solid line], which can be seen easily by matching the results of the two
panels of Fig. 5. However, we observe that not only $\alpha$ and $\beta$, but
also $n_{0,d}$ and $\psi_{0,\mathrm{P}}$ are invariant with respect to the
choice of the value of $\eta$ if taken precisely along the PSF line
($\mu-\mu_{c}=0$). Moreover, as a function of $2\,d\,J/U$, the width of the
PSF region shrinks with decreasing $2\,d\,J/U$, collapsing onto the PSF line
in the strongly-interacting limit. Therefore, an infinitesimal symmetry-
breaking perturbation $\eta\to 0^{+}$ can be safely applied for practical
purposes (cf. Ref. [52]). In the remainder of this work, this limit is assumed
whenever the PSF phase is discussed, in virtue of its essential insensitivity
to the choice of $\eta$.

Figure 6:  Excitation spectra (a) along the PSF-to-SF second-order transition
with $n_{0,d}\approx 1.47$ and $(\mu/U)_{c}=-0.35$ and (b) in the SF phase for
$\mu/U\approx-0.35$ and $U_{12}/U=-0.7$ in $d=2$. In panel (a), the hopping
energies are $2\,d\,J/U=0.08$ (red lines) and $2\,d\,(J/U)_{c}\approx 0.131$
(blue lines). In the superfluid phase, we consider $2\,d\,J/U=0.14$. The
parentheses refer to the degeneracy of the hybridized quasiparticle bands.

Along the PSF line, Eq. (13) can be treated analytically to obtain the
excitation spectrum. For the $n_{0,d}\approx 1.47$ PSF phase separating the
vacuum from the $n_{0,d}=2$ MI lobe, one obtains the implicit equation

$\displaystyle\left\\{J^{2}({\bf
k})\left[\left(8\,\beta^{4}-8\,\beta^{2}+1\right)\mu^{2}+\left(1-2\,\beta^{2}\right)^{2}U^{2}+2\left(4\,\beta^{4}-2\,\beta^{2}-1\right)\mu\,U-\hbar^{2}\omega^{2}_{\bf
k}\right]\right.$ $\displaystyle\left.+2\,J({\bf
k})\left[\mu\,U^{2}+2\,\beta^{2}\,\mu^{2}\,U-2\left(\beta^{2}-1\right)U\,\omega^{2}_{\bf
k}+\left(2\,\beta^{2}-1\right)\mu\left(\mu^{2}-\hbar^{2}\omega^{2}_{\bf
k}\right)\right]\right.$
$\displaystyle\left.+\left(\mu^{2}-\hbar^{2}\omega^{2}_{\bf
k}\right)\left[\left(\mu+U\right)^{2}-\hbar^{2}\omega^{2}_{\bf
k}\right]\right\\}^{2}=0\,,$ (33)

which describes two distinct gapped bands, both doubly degenerate and with a
mixed quasiparticle characteristic. In analogy with the CFSF phase, all the
remaining excitations have a flat band dispersion; in particular, the zero-
energy modes $\hbar\,\omega_{0,0}=\hbar\,\omega_{1,1}=0$ describe cost-free
density fluctuations due to the formation of local particle-hole pairs in the
PSF ground state predicted by Gutzwiller mean-field theory.

Once again, the second-order phase transition from the PSF and to the SF
phases is identified by the closure of the gap in the excitation spectrum.
From Eq. (2.2.2), we find that along the $n_{0,d}\approx 1.47$ PSF line, the
critical hopping strength is given by

$2\,d\left(\frac{J}{U}\right)_{c}=\frac{\mu}{U}\frac{1+\mu/U}{\mu/U\left(2-2\,\beta^{2}\right)-\left(1+\mu/U\right)\left(2\,\alpha\,\beta+1\right)}\,.$
(34)

In Fig. 6(a), we show how the excitation spectrum of the PSF phase of density
$n_{0,d}\approx 1.47$ evolves as the second-order transition point
$2\,d\,(J/U)_{c}\approx 0.131$ is reached for $U_{12}/U=-0.7$ along the
critical line at $(\mu/U)_{c}=-0.35$. At the PSF-to-SF transition point, only
the energy gap of the lowest quasiparticle branches vanishes, such that the
band structure appears reminiscent of the MI-to-SF and CFSF-to-SF second-order
tip transitions in Fig. 2(b) and Fig. 3(b) respectively but with the two
gapless modes being degenerate. More distinctly, the O$(2)$ nature of the PSF-
to-SF transition is suggested in Fig. 6(b), where we observe that the gapless
quasiparticle band splits into the Goldstone and Higgs modes active in the
spin channel, while the density Goldstone mode emerges from the non-dispersive
band $\omega_{0,0}=\omega_{1,1}$.

It is important to note that the flatness of the lowest spin bands
$\omega_{0,1}$ and $\omega_{1,0}$ for CFSF and density bands $\omega_{0,0}$
and $\omega_{1,1}$ for PSF is an artifact resulting from the mean-field
approximation. A more careful treatment of non-local fluctuations would alter
these modes to produce linear Goldstone dispersions as a result of the broken
$U(1)$ spin or density symmetries, respectively.

###  2.3 Quantum Gutzwiller theory for mixtures

In this subsection, we briefly review the quantum Gutzwiller (QGW) theory
developed in Ref. [35] for the one-component BH model and extend it to the
two-component case. The quantization of the Gutzwiller theory provides us with
a simple tool for going beyond the mean-field approximation, which proves to
be of key importance in determining the structure of non-trivial correlation
functions and the superfluid components of the system.

Following the derivation outlined in Ref. [35], the essential idea behind the
QGW method is to include quantum fluctuations on top of the Gutzwiller state
via a canonical quantization [53, 54] of the coordinates of the Lagrangian
(2.1). Specifically, we promote the Gutzwiller variables and their conjugate
momenta to operators as
$c_{n_{1},n_{2}}{\left(\mathbf{r}\right)}\to\hat{c}_{n_{1},n_{2}}{\left(\mathbf{r}\right)}$
and impose equal-time commutation relations between them, namely

$\left[\hat{c}_{n_{1},n_{2}}{\left(\mathbf{r}_{1}\right)},\hat{c}^{\dagger}_{m_{1},m_{2}}{\left(\mathbf{r}_{2}\right)}\right]=\delta_{\mathbf{r}_{1},\mathbf{r}_{2}}\,\delta_{n_{1},n_{2}}\,\delta_{m_{1},m_{2}}\,.$
(35)

We stress that, in the same way as the Gutzwiller ansatz (1) assigns a weight
to each local configuration, the corresponding Gutzwiller operators cannot be
decomposed into single-species operators without overlooking a relevant
fraction of the interspecies correlations: this sharply contrasts with
Bogoliubov’s theory [55, 33], where the quantum fields of different species
are always treated separately. It follows that the two-component QGW approach
can take into accurate account local pair correlations, while quantization
allows for an approximate view on the non-local quantum correlations missed by
mean-field theory.

In analogy with the number-conserving approaches in dilute ultracold atomic
gases [56, 57], the next step consists in expanding the operators
$\hat{c}_{n_{1},n_{2}}$ around their ground state $\mathbb{C}$-values
$c^{0}_{n_{1},n_{2}}$ as

$\hat{c}_{n_{1},n_{2}}{\left(\mathbf{r}\right)}=\hat{A}{\left(\mathbf{r}\right)}\,c^{0}_{n_{1},n_{2}}+\delta\hat{c}_{n_{1},n_{2}}{\left(\mathbf{r}\right)}\,.$
(36)

The normalization operator $\hat{A}{\left(\mathbf{r}\right)}$ is a functional
of the fluctuation fields $\delta\hat{c}_{n_{1},n_{2}}\left(\mathbf{r}\right)$
and is defined in order to ensure that the physical constraint
$\sum_{n_{1},n_{2}}\hat{c}^{\dagger}_{n_{1},n_{2}}{\left(\mathbf{r}\right)}\,\hat{c}_{n_{1},n_{2}}{\left(\mathbf{r}\right)}=\hat{\mathds{1}}$
is fulfilled exactly. By restricting ourselves to local fluctuations
orthogonal to the ground state
$\sum_{n_{1},n_{2}}\delta\hat{c}^{\dagger}_{n_{1},n_{2}}{\left(\mathbf{r}\right)}\,c^{0}_{n_{1},n_{2}}=0$
as usual, one has

$\hat{A}{\left(\mathbf{r}\right)}=\left[\hat{\mathds{1}}-\sum_{n_{1},n_{2}}\delta\hat{c}^{\dagger}_{n_{1},n_{2}}{\left(\mathbf{r}\right)}\,\delta\hat{c}_{n_{1},n_{2}}{\left(\mathbf{r}\right)}\right]^{1/2}\,.$
(37)

Physically, $\hat{A}{\left(\mathbf{r}\right)}$ in Eqs. (36)-(37) serves to
deplete the weight of the ground state parameters $c^{0}_{n}$ in response to
the onset of quantum fluctuations. This role is especially important in
determining quantum correlations in the strongly-interacting regime.

The Fourier transform of the Gutzwiller operators is given by

$\delta\hat{c}_{n_{1},n_{2}}{\left(\mathbf{r}\right)}\equiv\frac{1}{\sqrt{I}}\sum_{\mathbf{k}}e^{i\mathbf{k}\cdot\mathbf{r}}\,\delta\hat{C}_{n_{1},n_{2}}{\left(\mathbf{k}\right)}\,.$
(38)

Inserting Eq. (38) in the quantized Gutzwiller Hamiltonian

$\hat{H}_{QGW}=\sum_{\bf
r}\left\\{-J\sum_{i=1}^{2}\sum_{j=1}^{d}\left[\hat{\psi}^{\dagger}_{i}{\left({\bf
r}\right)}\,\hat{\psi}_{i}{\left(\bf
r+e_{j}\right)}+\mathrm{h.c.}\right]+\sum_{n_{1},n_{2}}H_{n_{1},n_{2}}\,\hat{c}^{\dagger}_{n_{1},n_{2}}{\left({\bf
r}\right)}\,\hat{c}_{n_{1},n_{2}}{\left({\bf r}\right)}\right\\}$ (39)

and keeping only terms up to the quadratic order in the fluctuations, we
obtain

$\hat{H}^{\left(2\right)}_{QGW}=E_{0}+\frac{1}{2}\sum_{\mathbf{k}}[\delta\underline{\hat{C}}^{\dagger}{\left(\mathbf{k}\right)},-\delta\underline{\hat{C}}{\left(-\mathbf{k}\right)}]\,\hat{\mathcal{L}}_{\mathbf{k}}\begin{bmatrix}\delta\underline{\hat{C}}{\left(\mathbf{k}\right)}\\\
\delta\underline{\hat{C}}^{\dagger}{\left(-\mathbf{k}\right)}\end{bmatrix}\,,$
(40)

where $E_{0}=\hbar\omega_{0}+z\,J\sum_{i=1}^{2}\left|\psi_{0,i}\right|^{2}$ is
the mean-field ground state energy, the vector
$\delta\underline{\hat{C}}(\mathbf{k})$ collects the Fock components
$\delta{\hat{C}}_{n_{1},n_{2}}(\mathbf{k})$ and
$\hat{\mathcal{L}}_{\mathbf{k}}$ coincides with the pseudo-Hermitian matrix
introduced in Eq. (13). The quadratic form (40) can be easily diagonalized by
a suitable Bogoliubov rotation of the Gutzwiller operators in terms of the
many-body excitation modes of the system,

$\delta\hat{C}_{n_{1},n_{2}}{\left(\mathbf{k}\right)}=\sum_{\alpha}u_{\alpha,\mathbf{k},n_{1},n_{2}}\,\hat{b}_{\alpha,\mathbf{k}}+\sum_{\alpha}v^{*}_{\alpha,-\mathbf{k},n_{1},n_{2}}\,\hat{b}^{\dagger}_{\alpha,-\mathbf{k}}\,.$
(41)

This allows a recasting of the Hamiltonian as

$\hat{H}^{\left(2\right)}_{QGW}=\sum_{\alpha}\sum_{\mathbf{k}}\hbar\omega_{\alpha,\mathbf{k}}\,\hat{b}^{\dagger}_{\alpha,\mathbf{k}}\hat{b}_{\alpha,\mathbf{k}}\,,$
(42)

where each mode $\hat{b}_{\alpha,\mathbf{k}}$ corresponds to a different
collective excitation with frequency $\omega_{\alpha,\mathbf{k}}$, labeled by
its branch index $\alpha$ and momentum $\mathbf{k}$. Bosonic commutation
relations between the annihilation and creation operators
$\hat{b}_{\alpha,\mathbf{k}}$ and $\hat{b}^{\dagger}_{\alpha,\mathbf{k}}$,

$\left[\hat{b}_{\alpha,\mathbf{k}},\hat{b}^{\dagger}_{\beta,\mathbf{p}}\right]=\delta_{\mathbf{k},\mathbf{p}}\,\delta_{\alpha,\beta}\,,$
(43)

are enforced by choosing the usual Bogoliubov normalization condition

$\underline{u}^{*}_{\alpha,\mathbf{k}}\cdot\underline{u}_{\beta,\mathbf{k}}-\underline{v}^{*}_{\alpha,-{\mathbf{k}}}\cdot\underline{v}_{\beta,-\mathbf{k}}=\delta_{\alpha,\beta}\,.$
(44)

The effective, quadratic description of the two-component BH model in terms of
its quasiparticle excitations provides a simple and versatile tool for the
computation of any expectation value. Drawing on the quantization procedure
described above and the Bogoliubov rotation (41), the evaluation of a generic
observable $\left\langle\hat{O}{\left[\hat{a}_{i,{\bf
r}},\hat{a}^{\dagger}_{i,{\bf r}}\right]}\right\rangle$ amounts to the
application of a straightforward calculation protocol first outlined in Ref.
[35] and included for completeness in App. A.1. Briefly, in a similar fashion
as the QGW quantization of the BH Hamiltonian (39), each observable is
projected along its Gutzwiller representation
$\hat{\mathcal{O}}{\left[\hat{c},\hat{c}^{\dagger}\right]}$ and expanded in
terms of the quantum fluctuations. Therefore, expectation values are
calculated via the simple application of Wick’s theorem to the quantized
excitations $\hat{b}_{\alpha,\mathbf{k}}$ with respect to the ground state of
the system.

In the following sections, the QGW quantization protocol is applied to
calculate a range of experimentally relevant one-body and two-body observables
in the two-component BH model. These observables are calculated always
assuming a zero-temperature vacuum of the quasiparticle modes of the theory
$|\Omega\rangle$, formally defined by
$\hat{b}_{\alpha,\mathbf{k}}|\Omega\rangle$ for each
$\left(\alpha,\mathbf{k}\right)$. We clarify that the quantization procedure
does not lead to a substantial change in the ground state and spectral
properties discussed in Secs. 2.2.1-2.2.2, but rather yields a transparent
formalism for accessing quantum fluctuations in a systematic manner. In order
to make our discussion of the QGW result consistent, from now on we calculate
the lattice filling $n_{i}$ by always including the second-order quantum
corrections as detailed in App. A. Further details of the application of the
quantization protocol to a specific observable can be found in App. A.

##  3 Linear Response

In this section, we investigate the linear response of the two-component BH
system to density/spin (Sec. 3.1) and current (Sec. 3.2) probes in two
dimensions. In the latter subsection, we then calculate the superfluid
components of the system within the linear response formalism outlined in Ref.
[27]. In particular, in Sec. 3.2.1, our QGW results for the collisionless drag
in the superfluid regime are compared against existing results available in
the literature. In Secs. 3.2.2 and 3.2.3, we study the current response
functions and superfluid components across the various phase transitions
discussed in Secs. 2.2.1 and 2.2.2.

Throughout the subsections below, we provide directly the semi-analytical
expressions of the response functions calculated within the QGW approach,
whose derivation is explicitly carried out in App. B. In general terms, we
consider the effect of an external field applied at a fixed frequency
$\omega$. The local operator associated with the perturbation is denoted by
$\hat{G}_{\bf r}$, and the local operator whose linear response dynamics we
are interested in studying is denoted by $\hat{F}_{\bf r}$.

###  3.1 Density and spin response

Figure 7: Normalized dynamical structure factors $\tilde{S}_{\hat{F}}({\bf
k},\omega)=S_{\hat{F}}({\bf k},\omega)/S_{\hat{F}}({\bf k})$ for (a)-(a′) the
one-component density channel $\hat{F}_{\bf r}=\hat{n}_{1,{\bf r}}$, (b)-(b′)
the total density channel $\hat{F}_{\bf r}=\hat{n}_{1,{\bf r}}+\hat{n}_{2,{\bf
r}}$ and (c)-(c′) the spin channel $\hat{F}_{\bf r}=\hat{n}_{1,{\bf
r}}-\hat{n}_{2,{\bf r}}$ in the SF phase in the vicinity of CFSF (a)-(c) and
PSF (a′)-(c′) transitions for $d=2$. For panels (a)-(c), the parameters are
$U_{12}/U=0.9$, $\mu/U=0.391$, $2\,d\,J/U=0.172$ and $n_{d}=1$, see Fig. 3(c).
For (a′)-(c′) panels, the parameters are $U_{12}/U=-0.7$, $\mu/U\approx-0.35$,
$2\,d\,J/U=0.14$ and $n_{d}\approx 1.47$, see Fig. 6(c).

##### Dynamic structure factors

In a two-component system, the most interesting response functions are the
density and the spin or magnetisation response functions, corresponding to
$\hat{G}_{\bf r}=\hat{F}_{\bf r}=\hat{n}_{1,{\bf r}}+\hat{n}_{2,{\bf r}}$ and
$\hat{G}_{\bf r}=\hat{F}_{\bf r}=\hat{n}_{1,{\bf r}}-\hat{n}_{2,{\bf r}}$,
respectively. Such response functions are related to the density and spin
structure factors of the system. Let us start from the species-resolved
density response function, i.e. $\hat{G}_{\bf r}=\hat{F}_{\bf
r}=\hat{n}_{i,{\bf r}}$, which, up to second-order in the quantum
fluctuations, reads (see App. B.2 for the derivation details)

$\chi_{\hat{n}_{i}}({\bf
q},\omega)=\frac{2}{\hbar}\,\sum_{\alpha}\frac{N^{2}_{i,\alpha,{\bf
q}}\,\omega_{\alpha,{\bf
q}}}{\left(\omega+i\,0^{+}\right)^{2}-\omega^{2}_{\alpha,{\bf q}}},$ (45)

at zero temperature. The response functions for the total density and spin
channels are obtained by simple extensions of Eq. (45), namely

$\chi_{\hat{n}_{d}}({\bf
q},\omega)=\frac{2}{\hbar}\,\sum_{\alpha}\frac{\left(N_{1,\alpha,{\bf
q}}+N_{2,\alpha,{\bf q}}\right)^{2}\omega_{\alpha,{\bf
q}}}{\left(\omega+i\,0^{+}\right)^{2}-\omega^{2}_{\alpha,{\bf q}}},$ (46)

and

$\chi_{\hat{n}_{s}}({\bf
q},\omega)=\frac{2}{\hbar}\,\sum_{\alpha}\frac{\left(N_{1,\alpha,{\bf
q}}-N_{2,\alpha,{\bf q}}\right)^{2}\omega_{\alpha,{\bf
q}}}{\left(\omega+i\,0^{+}\right)^{2}-\omega^{2}_{\alpha,{\bf q}}}\,,$ (47)

respectively. At zero temperature, the imaginary part of the response
functions is proportional to the corresponding dynamic structure factors via
the relation $S_{\hat{F}}({\bf q},\omega)=-\Im[\chi_{\hat{F}}({\bf
q},\omega)/\pi]$. Useful information on and from the dynamical structure
factors is provided by their energy momenta, a.k.a. sum rules

$m^{p}_{\hat{F}}(\mathbf{q})=\int_{0}^{+\infty}d\omega\,\omega^{p}\,S_{\hat{F}}({\mathbf{q}},\omega)\,.$
(48)

In particular, $m^{0}_{\hat{F}}(\mathbf{q})=S_{\hat{F}}(\mathbf{q})$ is the
so-called static structure factor. Within our approximation in Eqs. (45)-(47),
the dynamic structure factors are simply given by a sum of Dirac delta
functions at the values of the quasiparticle energies
$\omega_{\alpha,\mathbf{q}}$,

$S_{\hat{n}_{d}}({\bf
q},\omega)=\frac{1}{\hbar}\,\sum_{\alpha}\left(N_{1,\alpha,{\bf
q}}+N_{2,\alpha,{\bf
q}}\right)^{2}\left[\delta{\left(\omega-\omega_{\alpha,{\bf
q}}\right)}-\delta{\left(\omega+\omega_{\alpha,{\bf q}}\right)}\right]\,,$
(49) $S_{\hat{n}_{s}}({\bf
q},\omega)=\frac{1}{\hbar}\,\sum_{\alpha}\left(N_{1,\alpha,{\bf
q}}-N_{2,\alpha,{\bf
q}}\right)^{2}\left[\delta{\left(\omega-\omega_{\alpha,{\bf
q}}\right)}-\delta{\left(\omega+\omega_{\alpha,{\bf q}}\right)}\right]\,,$
(50)

and therefore the sum rules are easily determined.

In Ref. [32], the normalized dynamical structure factor
$\tilde{S}_{\hat{F}}({\bf q},\omega)=S_{\hat{F}}({\bf
q},\omega)/S_{\hat{F}}({\bf q})$ (cf. Ref. [33]) has been analyzed in the SF
regime for both repulsive and attractive interspecies interactions, for which
it was found that: (i) the low-momentum part of the gapped modes does not
respond significantly to any of the density-type probes, in agreement with the
single-component case [46, 29]; (ii) the density and spin Goldstone modes are
strongly excited by the total density and spin probes, respectively; (iii) the
single-species density fluctuations are strongest for the lowest gapless mode.
In [32] the density response in the PSF and CFSF regimes were also considered,
although for an effective Hamiltonian including ad hoc perturbative hopping
processes in the strongly-interacting limit of the model.

In Fig. 7, we show our results for the normalized dynamical structure factors
in the immediate vicinity of the CFSF (a)-(c) and PSF (a′)-(c′) transitions,
confirming the qualitative picture of Ref. [32] and displaying detailed
signatures of the pairing phase transitions. These calculations have been
performed for a total filling $n_{d}=1$ which includes the second-order
quantum corrections to the mean-field filling $n_{0,d}$ as detailed in App.
A.3. We carry this notation throughout the remainder of the work. Starting
from panels (a)-(a′), we illustrate that single-species perturbations are
sufficient for testing the proximity of the antipaired and paired phases, as
$\tilde{S}_{\hat{F}}({\bf q},\omega)$ receives a dominant contribution from
the lowest-lying Goldstone mode, while both the density and spin Goldstone
excitations are found to have approximately the same weight in the deep SF
region. More precisely, close to the CFSF (PSF) transition, the spin (density)
mode (which softens at the transition) dominates the system response and
enhances the amplitude of the structure factor at low energies for all
momenta. Moving to the total density channel, in panel (b) we see that,
despite the CFSF phase being dominated by spin fluctuations, the structure
factor receives an increasing contribution by the density Goldstone mode at
low momenta as the critical point is approached; as expected, the same mode
controls entirely the total density response of the system close to the PSF
transition, depicted in panel (b′). The situation is reversed in the case of
the spin channel, considered in panels (c)-(c′): here, the projection of the
structure factor over the spin Goldstone mode acts as a marker of both the
CFSF and the PSF transitions, albeit on different momenta and energy ranges.

Figure 8: Upper panels: dependence of the density (blue circles) and spin (red
circles) sound velocities in the vicinity of the O$(2)$ (a) CFSF and (b) PSF
transitions. Lower panels: dependence of the compressibility (blue diamonds)
and susceptibility (red diamonds) close to the same phase transitions. The
data in (a) and (a′) correspond to fixed $n_{d}=1$ and $U_{12}/U=0.9$, while
the data in (b) and (b′) are derived for fixed $n_{d}\approx 1.47$ and
$U_{12}/U=-0.7$; the calculations are always performed for $d=2$.

##### Static response and sum rules

The energy moments of the dynamic structure factors allow to obtain a number
of important quantities and identities that characterize the system. First of
all, the uniform limit of the inverse-energy-weighted sum rule
$m^{-1}_{\hat{F}}(\mathbf{q})$ gives the static response of the system to the
selected perturbation. Using the QGW expressions of Eq. (46), we obtain the
relations

$m^{-1}_{d}(\mathbf{q})=\frac{1}{\hbar}\sum_{\alpha}\frac{\left(N_{1,\alpha,{\bf
q}}+N_{2,\alpha,{\bf q}}\right)^{2}}{\omega_{\alpha,{\bf q}}}\underset{{\bf
q\to 0}}{=}\frac{\kappa}{2}\,,$ (51) $\quad
m^{-1}_{s}(\mathbf{q})=\frac{1}{\hbar}\sum_{\alpha}\frac{\left(N_{1,\alpha,{\bf
q}}-N_{2,\alpha,{\bf q}}\right)^{2}}{\omega_{\alpha,{\bf q}}}\underset{{\bf
q\to 0}}{=}\frac{\chi}{2}$ (52)

for the compressibility $\kappa$ and the susceptibility $\chi$ of the two-
component BH system, where the subscripts $d$ and $s$ are just shorthand
notations for $\hat{F}=\hat{G}=\hat{n}_{d/s}$. In particular, we note that the
former compressibility relation generalizes the result of Ref. [29] to
mixtures.

Our results for the compressibility and the spin susceptibility are shown in
Fig. 8 in the vicinity of the CFSF transition [panel (a)] and the PSF
transition [panel (b)]. Notably, we find that the spin susceptibility
(compressibility) diverges near the CFSF (PSF) transition. Close to the CFSF
regime, this finding parallels the decreasing energy cost to produce spin
excitations (see Fig. 3), which corresponds to an increase in the response of
the system towards magnetic perturbations. On the other hand, in the PSF
regime, the divergence of the compressibility corresponds to the decreasing
energy of density excitations (see Fig. 6), hence the increasing sensitivity
of the system to density fluctuations. Additionally, we find that the
compressibility (susceptibility) tends to vanish near the CFSF (PSF)
transition due to the opening of a spectral gap in the density (spin) channel.

The divergence of the static response functions – which would suggest an
instability of the system towards phase separation or collapse – is due the
lack of a proper inclusion of pairing quantum correlations when describing the
CFSF and the PSF phases within the Gutzwiller approximation: specifically, it
is simply related to the presence of a zero-energy flat dispersion relation
for the density and the spin modes, respectively. As shown in Fig. 8(a’)-(b’),
the appearance of such modes reflects into the vanishing behaviour of the spin
(density) sound velocities at the CSF (PSF) critical points.

An explicit relationship between sound velocities and static response
functions can be directly uncovered by making use again of sum rules. Indeed,
upon approaching the CFSF (PSF) phase, the low-momentum response function is
exhausted by the spin (density) Goldstone mode, as shown in Fig. 7. Thus, the
sum rules satisfy the general relation

$m^{p}_{d(s)}(\mathbf{q}\rightarrow
0)\simeq\omega^{p-k}_{d(s),\mathbf{q}}\,m^{k}_{d(s)}\simeq\left(c_{d(s)}|\mathbf{q}|\right)^{p-k}m^{k}_{d(s)}(\mathbf{q})\,,$
(53)

from which a number of identities can be derived. For instance, by considering
$p=0$ and $k=-1$, one can write the identities

$m^{0}_{d}({\bf q},\omega)=\left(N_{1,d,{\bf q}}+N_{2,d,{\bf
q}}\right)^{2}\underset{{\bf q\to 0}}{=}\frac{\kappa}{2}c_{d}\left|{\bf
q}\right|$ (54)

and

$m^{0}_{s}({\bf q},\omega)=\left(N_{1,s,{\bf q}}-N_{2,s,{\bf
q}}\right)^{2}\underset{{\bf q\to 0}}{=}\frac{\chi}{2}c_{s}\left|{\bf
q}\right|\,,$ (55)

which we verified numerically and generalize the single-component identity
given in Ref. [29] relating the Goldstone mode structure factor to the
compressibility and the sound speed. The latter equations are moreover not
independent. In fact, within the present Bogoliubov-like approach, the $p=1$
sum rule (a.k.a. $f$-sum rule) is given its exact value [33]:
$m^{1}_{d(s)}(\mathbf{q})=\langle[\delta\hat{n}_{d(s)}(\mathbf{q}),[\hat{H},\delta\hat{n}_{d(s)}(\mathbf{q})]]\rangle\propto{\bf
q}^{2}$ for $|{\bf q}|\rightarrow 0$. From this, using the single-mode sum
rule relation $m^{1}_{d(s)}(\mathbf{q})=\omega^{2}_{d(s),\bf
q}\,m^{-1}_{d(s)}(\mathbf{q})$, we immediately infer that $c^{2}_{d(s)}\propto
1/\kappa\,(1/\chi)$. Therefore, if the sound dispersion becomes flat, the
corresponding static response must diverge.

###  3.2 Current response and superfluid components

In this subsection, we consider the transverse (i) intraspecies current
response with $\hat{G}_{\bf r}=\hat{F}_{\bf r}=\hat{j}_{i}$ and (ii)
interspecies current response with $\hat{F}_{\bf r}=\hat{j}_{1}$ and
$\hat{G}_{\bf r}=\hat{j}_{2}$, evaluated at $\omega=0$ along $x$-directed
links of the square lattice. Here, $\hat{j}_{i}$ is the current operator
referred to $i^{\text{th}}$ species taken in the uniform limit ${\bf q\to 0}$,
that is

$\hat{j}_{i,{\bf q\to 0}}|_{x}=\frac{\hbar}{m\,a}\sum_{\bf
k}\sin(k_{x}\,a)\,\hat{a}^{\dagger}_{i,{\bf k}}\,\hat{a}_{i,{\bf k}}\,.$ (56)

In App. B.3, the intra/interspecies current response functions are derived
within the QGW formalism, giving the results

$\displaystyle\chi^{T}_{\hat{j}_{i},\hat{j}_{i}}({\bf q\to
0},\omega=0)=-\frac{\hbar}{\left(m\,a\right)^{2}}\sum_{\alpha,\beta}\sum_{\bf
k}\frac{\left(U_{i,\alpha,{\bf k}}\,V_{i,\beta,{\bf k}}-U_{i,\beta,{\bf
k}}\,V_{i,\alpha,{\bf k}}\right)^{2}}{\omega_{\alpha,{\bf
k}}+\omega_{\beta,{\bf k}}}\sin^{2}(k_{x}\,a)\,,$ (57)
$\displaystyle\chi^{T}_{\hat{j}_{1},\hat{j}_{2}}({\bf q\to
0},\omega=0)=-\frac{\hbar}{\left(m\,a\right)^{2}}\sum_{\alpha,\beta}\sum_{\bf
k}\frac{\prod_{i=1}^{2}\left(U_{i,\alpha,{\bf k}}\,V_{i,\beta,{\bf
k}}-U_{i,\beta,{\bf k}}\,V_{i,\alpha,{\bf k}}\right)}{\omega_{\alpha,{\bf
k}}+\omega_{\beta,{\bf k}}}\sin^{2}(k_{x}\,a)\,,$ (58)

respectively. The above equations naturally generalize the findings of Ref.
[27] provided by the Bogoliubov approximation. In that case, one has only spin
and density Goldstone modes given by Eq. (27) and consequently the current
response functions satisfy the relation
$\chi^{T}_{\hat{j}_{1},\hat{j}_{1}}({\bf q\to
0},\omega=0)=-\chi^{T}_{\hat{j}_{1},\hat{j}_{2}}({\bf q\to 0},\omega=0)$
exactly. Due to the presence of additional excitation bands with a substantial
spectral weight for strong enough interactions, within the QGW approach, the
very same equality is approximately fulfilled only in the deep SF phase.
Physically, the violation of aforementioned identity reflects the breaking of
Galilean invariance (c.f. Ref. [27]), which allows for a non-zero normal
component of the gas contributing to superfluidity even at zero temperature,
as we discuss in the following.

We proceed to introduce the relevant superfluid quantities in the two-species
BH model. For a homogeneous two-component superfluid system, the relations
between the mass current densities $m_{i}\,{\bf j}_{i}$ of each species and
the velocities of the gas components are governed by the two-fluid model
originally introduced in Ref. [11] and read

$m_{1}\,{\bf j}_{1}=\rho_{n,1}\,{\bf v}_{n}+\rho_{s,1}\,{\bf
v}_{1}+\rho_{12}\,{\bf v}_{2}\,,\\\ $ (59a) $m_{2}\,{\bf
j}_{2}=\rho_{n,2}\,{\bf v}_{n}+\rho_{s,2}\,{\bf v}_{2}+\rho_{12}\,{\bf
v}_{1}\,,$ (59b)

where $m_{i}$ is the bare mass of the $i^{\text{th}}$ species of the system,
$\rho_{s,i}$ and ${\bf v}_{i}$ are respectively the superfluid (mass)
densities and velocities for each component, $\rho_{12}$ is the so-called drag
interaction and $\rho_{n,i}$ are the normal (mass) densities of the system,
which are assumed to flow both at the same velocity ${\bf v}_{n}$. The
physical role of the superfluid drag, going under the name of Andreev-Bashkin
effect, has a self-evident explanation: a superflow in one component can be
induced by the collisionless drag from the superflow in the other component of
the system and vice versa. For a continuous system, by Galilean invariance
Eqs. (59) are supplemented by a close relationship
$m_{i}\,N_{i}/V=\rho_{n,i}+\rho_{s,i}+\rho_{12}$ between the normal part of
the system and the superfluid densities, such that at zero temperature, where
$\rho_{n,i}=0$, the whole volume density of the system $(N_{1}+N_{2})/V$
participates in the superfluid flow [58]. On a lattice, the breaking of
Galilean invariance requires the mass current densities to satisfy a different
transformation rule in presence of a vector potential acting as a probe, such
that the density $N_{i}/V$ is replaced by $-K_{i}/\left(2\,J\,a^{d}\right)$,
where $K_{i}=\big{\langle}\hat{K}_{i}\big{\rangle}$ is the local kinetic
energy referring to the $i^{\text{th}}$ species along the direction of the
phase twist induced by the vector potential [27]. In particular, the local
kinetic operator acting along $x$-directed links of a square lattice is given
by

$\hat{K}_{i}({\bf r})|_{x}=-J\left(\hat{a}^{\dagger}_{i,{\bf r}+{\bf
e}_{x}}\,\hat{a}_{i,{\bf r}}+\mathrm{h.c.}\right)\,.$ (60)

Most importantly, we also remark that the normal component $n_{n,i}$ may not
vanish at zero temperature on a lattice. Therefore, recovering a result
derived in Ref. [21], for the two-component BH model we obtain

$-\frac{K_{i}|_{x}}{2J}=n_{n,i}+n_{s,i}+n_{12}\,,$ (61)

where we have rescaled the superfluid (normal) number densities
$n_{s(n),i}\equiv a^{d}\rho_{s(n),i}/m$ and the drag $n_{12}\equiv
a^{d}\,\rho_{12}/m$ to facilitate a direct comparison to the filling factors
on a lattice. The average kinetic energy $K_{i}|_{x}$ within the QGW theory is
calculated explicitly in App. B.4 and reads

$K_{i}|_{x}\approx-2\,J\left|\psi_{0,i}\right|^{2}-\frac{2\,J}{I}\sum_{\alpha}\sum_{\bf
k}\left|V_{i,\alpha,{\bf k}}\right|^{2}\,\cos(k_{x}\,a)\,.$ (62)

To quantify how the superfluid drag impacts on the effective mass of the gas
components, we define the dimensionless parameter $\xi^{*}_{i}$ through the
equations

$\displaystyle n_{n,i}+n_{s,i}=-\xi^{*}_{i}\frac{\,K_{i}|_{x}}{2J}\,,$ (63)
$\displaystyle n_{12}=-\frac{K_{i}|_{x}}{2J}\left(1-\xi^{*}_{i}\right)\,,$
(64)

where $\xi^{*}=m/m^{*}$ reflects the extent to which the bare effective mass
$m=\hbar^{2}/\left(2\,J\,a^{2}\right)$ on a lattice is renormalized by the
interaction between the superfluid flows to produce the renormalized effective
mass $m^{*}$. Explicitly, when $\xi^{*}_{i}>1$ ($\xi^{*}_{i}<1$) the
renormalized effective mass is smaller (larger) than the bare effective mass.
This extends the concept of the renormalized effective mass discussed in Refs.
[11, 22] to the lattice, as Eq. (64) is exactly analogous to Eq. (4) in the
second work. We note that in those works, translational invariance ensures
that the renormalized effective mass remains always larger than the bare
effective mass due to the guaranteed positivity of $n_{12}$. On a lattice, the
inclusion of $n_{n,i}$ in Eq. (63), which is trivially absent for a continuous
superfluid at zero temperature, adds additional complexity while ensuring that
the variation of $\xi^{*}_{i}$ is due solely to the collisionless drag.

Having introduced the relevant superfluid components, we now resort to the
formal results of Ref. [27] to relate them to the current-current response
functions of the two-component BH model. Explicitly, one has

$\displaystyle n_{s,i}=\frac{m}{I}\chi^{T}_{\hat{j}_{i},\hat{j}_{i}}({\bf q\to
0},\omega=0)-\frac{K_{i}|_{x}}{2J}\,,$ (65) $\displaystyle
n_{12}=\frac{m}{I}\chi^{T}_{\hat{j}_{1},\hat{j}_{2}}({\bf q\to
0},\omega=0)\,.$ (66)

A third relation provides a sum rule for the total normal fraction of the
system $n_{n}=n_{n,1}+n_{n,2}$ as the total transverse current response
function, reading

$n_{n}=-\frac{m}{I}\left[\chi^{T}_{\hat{j}_{1},\hat{j}_{1}}({\bf q\to
0},\omega=0)+\chi^{T}_{\hat{j}_{2},\hat{j}_{2}}({\bf q\to
0},\omega=0)+2\,\chi^{T}_{\hat{j}_{1},\hat{j}_{2}}({\bf q\to
0},\omega=0)\right]\,.$ (67)

From Eq. (67), we notice how, upon exchanging $\hat{j}_{1}$ with
$\hat{j}_{2}$, opposite values of the current response functions
$\chi^{T}_{\hat{j}_{1},\hat{j}_{1}}({\bf q\to
0},\omega=0)=-\chi^{T}_{\hat{j}_{1},\hat{j}_{2}}({\bf q\to 0},\omega=0)$
results in $n_{n}=0$, as would be predicted by the Bogoliubov approximation,
in clear contrast with the correct physics of a strongly-interacting lattice
system at zero temperature.

It is worth mentioning that measuring the superfluid drag is a challenging
task experimentally, however recently it has been proposed to use fast
response to directly access the entrainment in continuous systems [59, 27]. We
note that the same proposals can be applied in the presence of a lattice.

####  3.2.1 Superfluid regime

Figure 9:  Collisionless drag versus the interspecies coupling strength
$U_{12}/U$ for fixed $n_{d}\approx 0.5$ and $2\,d\,J/U=0.4$ in $d=2$. Here, we
compare the QGW prediction (orange solid line) with the Bogoliubov result from
Refs. [21, 27] (light-blue dashed line) and the QMC data (yellow dotted line)
from Ref. [25] evaluated on a lattice of size $I=10^{2}$. The hollow square
and circle for the QGW and Bogoliubov results, respectively, indicate the
point of the collapse.

In the deep SF regime, the current response functions and the average kinetic
energy match their expressions given by the mean-field Gutzwiller
approximation, namely $\chi^{T}_{\hat{j}_{i},\hat{j}_{j}}({\bf q\to
0},\omega=0)=0$ and $K_{0,i}|_{x}=-2\,J\left|\psi_{0,i}\right|^{2}$
respectively. This leads to equal superfluid and condensate densities, namely
$n_{s,i}=\left|\psi_{0,i}\right|^{2}$, as well as a vanishing drag $n_{12}=0$
and a trivial renormalization of the effective mass $\xi^{*}\approx 1$.

At intermediate $2\,d\,J/U$, the superfluid drag was calculated using quantum
Monte Carlo (QMC) simulations in Ref. [25]. The comparison of the QGW results
with the Bogoliubov predictions (see Refs. [21, 27]) and the QMC data is shown
in Fig. 9, for a hopping energy equal to $2\,d\,J/U=0.4$ and a fixed total
filling $n_{d}\approx 0.5$. We remark here that the accuracy of our drag
estimation hinges on the correct calculation of the total filling, which must
include second-order quantum corrections accounted for by the QGW theory (see
App. A.3). For repulsive interactions, these corrections are always less than
$\sim 10\%$ and therefore can be essentially neglected. By contrast, for
attractive interactions quantum corrections increase for larger
$\left|U_{12}\right|$ and can be as large as $\sim 25\%$ due to the diverging
compressibility near the collapse transition (see Sec. 2.2.2). We find that
including the quantum corrections has the effect of shifting the collapse
transition towards more attractive $U_{12}$ than one would obtain holding
instead the mean-field filling constant. This fact restricts the calculation
of the drag to well before the Bogoliubov prediction for the collapse point
($U_{12}/U\leq-1$) [27, 26].

The QGW results of Fig. 9 [orange solid line] are obtained on the same
$I=10^{2}$ lattice considered in the QMC calculations of Ref. [25] [yellow
dotted line]. We find that in general the drag increases for smaller lattice
sizes, in contrast to the results of that work. Furthermore, we also observe
that the collapse point is shifted towards less attractive $U_{12}$ for
smaller lattices. Additionally, we note that, regardless of the sign of
$U_{12}$, the collisionless drag remains always positive in the SF regime, and
consequently the renormalized effective mass is larger than the bare one.
Physically, this indicates that the dressing effect of individual particle
motions of one species due to the superfluid flow of the other component has a
particle (rather than hole) character.

On the whole, the QGW predictions underestimate the QMC results, but
qualitatively reproduce the asymmetry of the superfluid drag with respect to
$U_{12}$ and are generally comprised within the estimated statistical
uncertainty of the QMC data. On the other hand, the Bogoliubov results [light-
blue dashed line] are symmetric with respect to $U_{12}$, do not capture the
first-order transition and lie well below the uncertainty window of the QMC
points. In this regard, the QGW approach can be viewed as a major improvement
over the standard Bogoliubov treatment of quantum fluctuations in presence of
strong correlations [34, 21, 27]. Indeed, it is well known that Bogoliubov’s
theory underestimates the current response functions, as it takes into account
only the excitation vertex of the Goldstone modes in Eqs. (57)-(58) and
neglects the contribution of all the other excitation modes, which acquire a
sizeable spectral weight away from the deep SF limit.

####  3.2.2 Phase transitions for $U_{12}>0$

In this subsection, we study the superfluid components across the various
phase transitions characterized in Sec. 2.2.1 for repulsive interspecies
interactions. To that end, we begin by discussing the first and second-order
MI-to-SF transitions, passing to the second-order CFSF-to-SF transition in the
second place.

Figure 10:  Transport quantities in two dimensions across the MI-to-SF second
(left panel) and first (right panel) order transitions at fixed $\mu/U=1.4$
and $\mu/U=1$ for $U_{12}/U=0.5$ and $U_{12}/U=0.9$, corresponding to the
$n_{d}=2$ MI lobes in Fig. 1(a) and Fig. 1(b), respectively. The one-component
densities $n_{i}$ include the corrections of the second-order quantum
fluctuations. The tan-shaded area indicates the MI region. Thick black and
grey lines indicate the QGW predictions for the superfluid density $n_{s,1}$
and the (vanishingly small) superfluid drag $n_{12}$ respectively, while the
black dashed line refers to the normal component $n_{n,1}$. Pink dotted and
blue dashed lines are the contributions to $n_{s,1}$ from the average kinetic
energy $K_{i}|_{x}$ and the intraspecies current response function
$\chi^{T}_{\hat{j}_{i},\hat{j}_{i}}$ respectively. (Insets) Renormalization of
the effective mass across the MI-to-SF transitions due to the collisionless
drag.

(Mott Insulator to Superfluid) – In Fig. 10(a) and Fig. 10(b), we show the QGW
results for the relevant transport quantities, which include the superfluid
components, the intraspecies current response, and the average kinetic energy
across the second and first-order MI-to-SF transitions of the $n_{0,d}=2$
lobes for $U_{12}/U=0.5$ and $U_{12}/U=0.9$, respectively. Qualitatively, the
behavior of the superfluid fraction exhibits a number of features in common
with the QGW results of Ref. [35] for the single-component BH model.
Specifically, in the SF regime, the superfluid density $n_{s,i}$ [black solid
line] remains always larger than the condensate fraction and approaches the
total density of the corresponding species in the deep SF limit. Furthermore,
$n_{s,i}$ vanishes discontinuously (continuously) at the first-order (second-
order) critical point and is exactly zero in the MI phase as in the single
component case. As found in the single-component BH model [35], this latter
feature is ensured by the exact cancellation between the current response
function $\chi^{T}_{\hat{j}_{i},\hat{j}_{i}}$ [pink dashed line], given by Eq.
(57), and the contribution of zero-point fluctuations to the average kinetic
energy $K_{i}|_{x}$ [cyan dashed line] provided by Eq. (62). As expected, both
the quantities tend towards zero in the strongly-interacting limit
$2\,d\,J/U\to 0$. The collisionless drag [gray solid line] remains always on
the order of a few percent of $n_{s,i}$, reaching the maximal value close to
the MI-to-SF critical point and vanishing entirely within the MI phase.
Because the superfluid drag is small, the renormalization of the effective
mass is also negligible and the normal components $n_{n,i}$ [black dashed
line] are dominated by the contribution of the intraspecies current response
function $\chi^{T}_{\hat{j}_{i},\hat{j}_{i}}$, as one can see by comparing
Eqs. (66)-(67).

It is now worth commenting on the relative weight of the various excitation
vertices that contribute the most to the results displayed in Fig. 10. In the
SF regime, the spin-density Goldstone vertex makes the largest contribution to
the expression of $\chi^{T}_{\hat{j}_{i},\hat{j}_{i}}$ in Eq. (57), followed
by smaller terms coming from the vertex between the Goldstone mode and the
Higgs mode with density character. This is in strong contrast to the single
component case treated in Ref. [35], which showed that
$\chi^{T}_{\hat{j}_{i},\hat{j}_{i}}$ was nearly saturated by the Goldstone-
Higgs vertex, since in this case only one Goldstone excitation is present.
Similarly, the interspecies current response
$\chi^{T}_{\hat{j}_{1},\hat{j}_{2}}$, and therefore the superfluid drag, is
nearly saturated by the vertex involving spin and density Goldstone modes. On
the other hand, the zero-point fluctuations in the average kinetic energy
$K_{i}|_{x}$ [cyan dashed line] are dominated by the contribution of the spin
Goldstone mode, in addition to the mean-field effect of the order parameter
$\psi_{0,i}$.

Inside the MI phase, the spectral summation of the intraspecies response
$\chi^{T}_{\hat{j}_{i},\hat{j}_{i}}$ is saturated by the vertices of particle-
hole excitations, while the average kinetic energy $K_{i}|_{x}$ is dominated
by the contribution of the lowest particle or hole bands, depending on the
chemical potential. As the MI-to-SF transition is crossed, these excitations
turn into the Goldstone modes, such that the spectral content of $K_{i}|_{x}$
changes smoothly across the criticality. An analogous reasoning applies to
$\chi^{T}_{\hat{j}_{i},\hat{j}_{i}}$. On the other hand, we note that the
interspecies response $\chi^{T}_{\hat{j}_{1},\hat{j}_{2}}$ is not fully
saturated by the vertices involving only the first few particle-hole bands,
requiring also the contribution of higher-energy excitations in order to give
a perfect cancellation of the superfluid drag in the MI lobe.

Figure 11: Transport quantities in two dimensions across the CFSF-to-SF
second-order edge transition, traversed at fixed chemical potential
$\mu/U=0.5$ for $U_{12}/U=0.9$, corresponding to the $n_{d}=1$ CFSF lobe in
Fig. 1(b). The one-component filling $n_{i}$ includes the corrections of the
second-order quantum fluctuations. The tan-shaded area indicates the MI
region. Thick black and grey lines indicate the QGW predictions for the
superfluid density $n_{s,1}$ and the superfluid drag $n_{12}$ respectively,
while the black dashed line refers to the normal component $n_{n,1}$. Pink
dotted and blue dashed lines are the contributions to $n_{s,1}$ from the
average kinetic energy $K_{i}|_{x}$ and the intraspecies current response
function $\chi^{T}_{\hat{j}_{i},\hat{j}_{i}}$ respectively. (Inset)
Renormalization of the effective mass across the CFSF transition due to the
collisionless drag.

(Counterflow Superfluid to Superfluid) – The QGW results for the superfluid
components, the intraspecies current response, and the average kinetic energy
across the second-order CFSF-to-SF transition of the $n_{0,d}=1$ lobe for
$U_{12}/U=0.9$ are shown in Fig. 11.

In the SF region close to the CFSF phase, the results are equivalent to the
ones found in the SF phase close to the MI lobe previously discussed. In the
CFSF regime, we find that the superfluid drag reaches the negative saturation
threshold $n_{12}/\sqrt{n_{s,1}\,n_{s,2}}=-100\%$, a result that strongly
recalls the QMC calculations in Ref. [25], performed in $d=2$ as well.
Notably, we also observe that in the CFSF phase, even though the drag is
saturated, the net superfluidity $n_{s,1}+n_{12}=0$ vanishes. This mirrors the
fact that counterflow superfluidity occurs when particles and holes of
different species flow along counter-directed paths, such that equal and
opposite current densities for each component are established [6, 25].
Importantly, within the QGW formalism this perfect balance is due solely to
quantum fluctuations, as the mean-field Gutzwiller theory trivially predicts
$n_{s,1}=n_{12}=0$. Qualitatively speaking, such finding agrees with the
statement in Ref. [32] that superfluidity in the CFSF phase arises through
second-order hopping processes not seized by mean-field theory; however, the
ad hoc introduction of such processes performed in that work is not explicitly
comprised in the O${\left(J^{2}\right)}$ contribution of quantum fluctuations
to Eqs. (57) and (58), which is thus a genuine result of the QGW quantum
theory. The collisionless drag remains large and negative throughout the CFSF-
to-SF transition, so that the renormalized effective mass becomes
significantly less than its bare value, particularly at the critical point. We
remark here that the onset of a negative superfluid drag indicates that a
traveling particle transports, in addition to it’s own bare mass, holes of the
other species resulting in a reduced effective mass. As the deep SF regime is
approached, the drag $n_{12}$ changes sign leading to an increased effective
mass as discussed in Sec. 3.2.1.

As before, we comment on the relative weight of the various excitations
vertices that contribute to the quantities shown in Fig. 11. In the CFSF
regime, the average kinetic energy $K_{i}|_{x}$ is almost saturated by the
gapped hole mode with the lowest energy. Correspondingly, both the
intraspecies $\chi^{T}_{\hat{j}_{i},\hat{j}_{i}}$ and interspecies
$\chi^{T}_{\hat{j}_{1},\hat{j}_{2}}$ current response functions are dominated
by the all-to-all vertices of the lowest four particle-hole bands. It is
curious to observe that the spectral weight of these particle-hole excitations
is sufficient to obtain the $-100\%$ saturation of the superfluid drag, even
when the antipair excitations are described by flat bands $\omega_{0,1}$,
which do not contribute to the transport physics in our theory. As the
transition point is crossed into the SF regime, $K_{i}|_{x}$ is saturated by
the density Goldstone mode, while the response functions
$\chi^{T}_{\hat{j}_{i},\hat{j}_{i}}$ and $\chi^{T}_{\hat{j}_{1},\hat{j}_{2}}$
get a large contribution by the vertices between the density Goldstone mode
and the first gapped modes of the SF phase.

####  3.2.3 Phase transitions for $U_{12}<0$

In this last subsection, we study the superfluid components across the the
second-order PSF-to-SF transition described in Sec. 2.2.2. The analysis of the
MI-to-SF critical behavior given in the previous subsection is found to remain
valid also in the case $U_{12}<0$, with the only difference that the physical
roles of the spin and density Goldstone modes are swapped. Therefore, in the
following we proceed to investigate in detail the behavior of the superfluid
components across the PSF-to-SF transition.

Figure 12: Transport quantities in two dimensions across the PSF-to-SF phase
transition along the $(\mu/U)_{c}=-0.35$ critical line, which separates the
vacuum region from the $n_{d}=2$ MI lobe for $U_{12}/U=-0.7$, see Fig. 4(b).
The one-component filling $n_{i}$ include the corrections of the second-order
quantum fluctuations. The tan-shaded area indicates the MI region. Thick black
and grey lines indicate the QGW predictions for the superfluid density
$n_{s,1}$ and the superfluid drag $n_{12}$ respectively, while the black
dashed line refers to the normal component $\rho_{n,1}$. Pink dotted and blue
dashed lines are the contributions to $n_{s,1}$ from the average kinetic
energy $K_{i}|_{x}$ and the intraspecies current response function
$\chi^{T}_{\hat{j}_{i},\hat{j}_{i}}$ respectively. (Inset) Renormalization of
the effective mass across the PSF transition due to the collisionless drag.

(Pair Superfluid to Superfluid) – The QGW predictions for the transport
quantities across the second-order PSF-to-SF for $U_{12}/U=-0.7$ with
$(\mu/U)_{c}=-0.35$ and $n_{0,d}\approx 1.47$ are presented in Fig. 12.

In the PSF regime, we find that the superfluid drag fulfills the positive
saturation condition $n_{12}/\sqrt{n_{s,1}\,n_{s,2}}=+100\%$, which is again
compatible with the QMC findings in Ref. [25]. Unlike the CFSF phase, however,
here the saturation of the drag results in net non-zero superfluidity
$n_{s,1}+n_{12}>0$. Physically speaking, this can be interpreted from the
point of view of pairs of particles of different species following co-directed
paths and resulting in pair superfluidity [25]. Similarly to the QGW
representation of the CFSF phase, this result can be directly ascribed to the
quantum fluctuations captured by our theory. Indeed, the collisionless drag
remains large and positive along the whole PSF line and through the transition
into the SF phase. Consequently, the effective mass is strongly renormalized,
so as to be significantly larger than the bare mass, in particular at the
transition point. We note that this increase of the effective mass is much
larger than in the SF and MI regimes due to the tendency of a traveling
particle to transport, in addition to its own mass, particles of the other
species that consequently follow co-directed paths and increasing the dressing
effect of the medium. We note that as the system approaches the deep SF
regime, the drag remains positive but becomes increasingly small such that the
effective mass always remains slightly larger than the bare mass as expected
from Sec. 3.2.1.

We conclude by inspecting the spectral composition of the various quantities
shown in Fig. 12. On the verge of the PSF side of the transition, the average
kinetic energy $\langle\hat{K}_{i}|_{x}\rangle$ is strongly dominated by the
lowest quasiparticle band (with hybrid spin and density character). The same
observation partially applies to both the current response functions
$\chi^{T}_{\hat{j}_{i},\hat{j}_{i}}$ and $\chi^{T}_{\hat{j}_{1},\hat{j}_{2}}$,
in which the coupling between different quasiparticle bands plays a major
role. As for the CFSF phase in Fig. 11, it is surprising to verify that these
gapped quasiparticle modes are sufficient to reproduce the superfluid drag
saturation of the PSF phase, such that they can be regarded as the only
excitations responsible for the transport of particle pairs. Once the critical
point is reached and the SF phase develops, $K_{i}|_{x}$ gets a major
contribution from the spin Goldstone mode, while the current response
functions are saturated by the coupling of the same mode with the Higgs-like
branches that populate the high-energy part of the spectrum, in strict analogy
with the case of repulsive interactions.

##  4 Correlation Functions

In the current section, we study how quantum fluctuations due to the
collective modes determine the one and two-body correlation functions of the
two-component BH model across its phase diagram in two dimensions ($d=2$).

###  4.1 Coherence function

In analyzing the single-particle coherence functions, we restrict ourselves to
the investigation of intraspecies correlations, as the approximate description
of the CFSF and PSF phases within the Gutzwiller framework [32] is expected to
miss non-local coherence effects involving strongly-correlated pairs, in
contrast with the results for the current response functions studied in Sec.
3.2.

The normalized single-particle coherence function for the $i^{\text{th}}$
species is defined as

$g^{(1)}_{i}({\bf r})\equiv\frac{\langle\hat{a}^{\dagger}_{i,{\bf
r}}\,\hat{a}_{i,{\bf 0}}\rangle}{\langle\hat{a}^{\dagger}_{i,{\bf
0}}\,\hat{a}_{i,{\bf 0}}\rangle}\,.$ (68)

The QGW quantization scheme maps the microscopic operator $\hat{a}_{i,{\bf
r}}$ into the effective Bose field $\hat{\psi}_{i}({\bf r})$, which carries
both a macroscopic contribution due to condensation and the effect of short-
range quantum correlations on the one-body coherence. In strict analogy with
Bogoliubov’s theory, the first-order expansion of $\hat{\psi}_{i}({\bf r})$ in
terms of quantum fluctuations reads

$\hat{\psi}_{i}({\bf
r})\approx\psi_{0,i}+\frac{1}{\sqrt{I}}\sum_{\alpha}\sum_{\bf
k}\left[U_{i,\alpha,{\bf k}}\,e^{i\,{\bf k}\cdot{\bf r}}\,\hat{b}_{\alpha,{\bf
k}}+V^{*}_{i,\alpha,{\bf k}}\,e^{-i\,{\bf k}\cdot{\bf
r}}\,\hat{b}^{\dagger}_{\alpha,{\bf k}}\right]\,,$ (69)

where the quasiparticle (quasihole) amplitudes $U_{i,\alpha,{\bf k}}$
$\left(V_{i,\alpha,{\bf k}}\right)$ have been introduced in Eqs. (18) and
saturate the Bogoliubov normalization condition

$\sum_{\alpha}\left(\left|U_{i,\alpha,{\bf
k}}\right|^{2}-\left|V_{i,\alpha,{\bf k}}\right|^{2}\right)=1$ (70)

throughout the whole phase diagram. Most importantly, the completeness
relation (70) implies that the QGW Bose field (69) satisfies the usual Bose
commutation relation $\left[\hat{\psi}_{i}({\bf
r}_{1}),\hat{\psi}^{\dagger}_{i}({\bf r}_{2})\right]=\delta_{{\bf r}_{1},{\bf
r}_{2}}$ identically, thus justifying the physical interpretation of
$\hat{\psi}_{i}({\bf r})$ (see App. A.2 for a detailed derivation). Applying
our evaluation protocol to Eq. (68) and inserting the operator expansion (69),
the single-particle coherence function can be recast into the form

$g^{(1)}_{i}({\bf r})=\frac{\langle\hat{\psi}^{\dagger}_{i}({\bf
r})\,\hat{\psi}_{i}({\bf 0})\rangle}{\langle\hat{\psi}^{\dagger}_{i}({\bf
0})\,\hat{\psi}_{i}({\bf
0})\rangle}\approx\frac{\left|\psi_{0,i}\right|^{2}+I^{-1}\sum_{\alpha}\sum_{\bf
k}\left|V_{i,\alpha,{\bf k}}\right|^{2}\cos({\bf k}\cdot{\bf
r})}{\left|\psi_{0,i}\right|^{2}+I^{-1}\sum_{\alpha}\sum_{\bf
k}\left|V_{i,\alpha,{\bf k}}\right|^{2}}\,,$ (71)

which is a straightforward generalization of the one-component result in Ref.
[35]. In the numerator of the right-hand side of Eq. (71), the first term
reflects the long-range order of the one-body density matrix in the SF phase,
while the second term reproduces the destructive interference of quantum
fluctuations at finite distances, such that only the condensate fraction
$\left|\psi_{0,i}\right|^{2}$ survives in the ${\bf r}\to\infty$ limit (see
the relevant discussion in Sec. 2.2 of Ref. [52]).

Figure 13: One-body coherence function $g^{(1)}_{i}({\bf r})$ for $d=2$ and
$U_{12}/U=0.9$ evaluated at fixed $\mu/U=1.4$ across the (a) MI to (b) SF
first-order transition. Figure 14: One-body coherence function $g^{(1)}({\bf
r})$ for $d=2$ and $U_{12}/U=0.9$ evaluated at fixed $\mu/U=0.5$ (a-b) and
$\mu/U\approx 0.39$ (c-d) across the edge and O(2) CFSF-to-SF transitions,
respectively. In panel (c), the approaching of the PSF-to-SF transition from
the PSF side is also shown for $\left(\mu/U\right)_{c}=-0.35$ and
$U_{12}/U=-0.7$.

We do not show here explicit results for the second-order MI-to-SF transition,
because it presents the same features as in the one-component case discussed
in [35]. Specifically, the QGW approach drastically improves mean-field theory
by predicting the onset of off-site coherence in the strongly-interacting
regime: in the MI phase, the one-body correlations are generally suppressed
exponentially as $g^{(1)}({\bf r})\sim\exp(-r/\xi)$ with a finite coherence
length $\xi$, whereas in the SF phase $g^{(1)}({\bf r})$ always decays as a
power law. More generally, in the deep SF limit ($2\,d\,J/U\gg 1$) the
spectral sum in Eq. (71) is almost saturated by the density and spin Goldstone
modes only and the behavior of a weakly-interacting gas is recovered [34, 35].
As a strongly-correlated SF develops, the contribution of other excitation
modes to the quantum depletion becomes relevant and the Bogoliubov predictions
are naturally amended by the QGW scenario.

In the two-component system, we find that the very same behavior carries over
also to the first-order MI-to-SF transition, shown in Fig. 13, and to the
second-order CFSF-to-SF transition and the PSF-to-SF critical point,
illustrated in Fig. 14. In particular, in analogy with the physics of the MI-
to-SF transition in the one-component system [35], the QGW theory is able to
capture the different critical behaviors of the CFSF-to-SF transition
depending on whether this is approached at integer or non-integer filling,
namely across either the tip of the edge of the CFSF lobes. Upon reaching the
SF phase through the edge of the CFSF lobe [Fig. 14(a), from purple to green
lines], the correlation length $\xi$ grows monotonically but remains bounded.
As soon as one enters the SF phase, the long-range behavior of
$g^{(1)}_{i}({\bf r})$ changes abruptly into a power-law scaling [Fig. 14(b)].
On the contrary, when approaching the SF phase at the tip of the CFSF region,
the correlation length $\xi$ diverges [Fig. 14(c)] and a power-law dependence
for $g^{(1)}_{i}({\bf r})$ gradually sets in [Fig. 14(d)]. In panel (c), the
evolution of $g^{(1)}_{i}({\bf r})$ along the PSF critical line is also shown,
displaying an analogous behavior.

As in the single particle case, this difference in the behavior of
$g^{(1)}_{i}({\bf r})$ is due to the different spectral properties of the
collective excitations close to the edge/tip critical points of the CFSF lobe.
At the edge transition, either a particle or a hole excitation out of the four
dispersive branches becomes gapless. Since our description of the short-
distance coherence in the CFSF phase relies on virtual particle-hole
excitations via Eq. (71), the exponential decay of $g^{(1)}_{i}({\bf r})$ is
dominated by the gap of the particle (hole) excitation that remains finite at
the transition. On the contrary, at the tip critical point, both the lowest-
energy particle and hole modes become gapless (before turning into the density
and spin Goldstone modes on the SF side of the transition), which explains the
divergent coherence length [60].

###  4.2 Density and spin fluctuations

Figure 15:  Local density-density correlation $g^{(2)}_{d}({\bf r}={\bf 0})$
evaluated for $d=2$ as a function of $2\,d\,J/U$ across (a) the first-order
MI-to-SF edge transition (blue dotted lines) for $U_{12}/U=0.9$ and
$\mu/U=1.4$ at fixed $n_{d}=2$, (b) the second-order MI-to-SF edge transition
(red solid line) for $U_{12}/U=0.5$ and $\mu/U=1$ at fixed $n_{d}=2$, (c) the
second-order CFSF-to-SF edge transition (purple dashed line) for
$U_{12}/U=0.9$ and $\mu/U=0.5$ at fixed $n_{d}=1$, and (d) the second-order
PSF-to-SF transition (pink dashed-dotted line) for $U_{12}/U=-0.7$ and
$(\mu/U)_{c}=-0.35$ with $n_{d}=1.47$. The point located on each curve
indicates the location of the respective phase transition. Figure 16: Local
spin-spin correlation $g^{(2)}_{s}({\bf r}={\bf 0})$ evaluated for $d=2$ as a
function of $2\,d\,J/U$ across the same phase transitions as in Fig. 15.

In this subsection, we address the local structure of two-body correlations
for both the density and spin channels, focusing on their behavior across the
quantum phase transitions of the system as determined by second-order quantum
fluctuations. Crucially, the QGW method treats local and non-local observables
separately. Whereas non-local two-body correlations $G^{(2)}_{d/s}{\left({\bf
r\neq 0}\right)}=\langle\hat{n}_{d/s}{\left({\bf r\neq
0}\right)}\,\hat{n}_{d/s}{\left({\bf 0}\right)}\rangle$ are directly related
to the Fourier transform of the static structure factors $S_{d/s}{(\bf q)}$
(and share then the pathologies discussed in Sec. 3.1), on-site fluctuations
can be always computed as expectation values of individual local operators, as
shown explicitly in the following. Ultimately, the Gutzwiller ansatz (1) is a
local approximation of the ground state, and therefore we expect predictions
presented in this subsection for local two-body correlations to remain
accurate even in the vicinity of the pair and anti-paired phases.

With these caveat in mind, we first consider the two-body correlation function
for an individual species

$G^{(2)}_{i}{\left({\bf r=0}\right)}=\langle\hat{n}^{2}_{i}{\left({\bf
0}\right)}\rangle\longrightarrow\langle\hat{\mathcal{D}}_{i}{\left({\bf
0}\right)}\rangle\,,$ (72)

where we have introduced the QGW square density operator

$\hat{\mathcal{D}}_{i}({\bf
r})=\sum_{n_{1},n_{2}}\left(n^{2}_{1}\,\delta_{i,1}+n^{2}_{2}\,\delta_{i,2}\right)\hat{c}^{\dagger}_{n_{1},n_{2}}({\bf
r})\,\hat{c}_{n_{1},n_{2}}({\bf r})\,,$ (73)

which, importantly, is distinct from the square of the density operator
$\hat{\mathcal{N}}_{i}({\bf r})$ (see App. A.2). Expanding
$\hat{\mathcal{D}}_{i}{\left({\bf 0}\right)}$ up to second-order in the
fluctuations and calculating its quantum average, we obtain

$G^{(2)}_{i}{\left({\bf
0}\right)}=\left(1-F\right)d_{0,i}+\frac{1}{I}\sum_{\alpha}\sum_{\bf
k}\sum_{n_{1},n_{2}}\left(n^{2}_{1}\,\delta_{i,1}+n^{2}_{2}\,\delta_{i,2}\right)\left|v_{\alpha,{\bf
k},n_{1},n_{2}}\right|^{2}\,,$ (74)

at zero temperature, which is a generalization of the result found in Ref.
[35] for the one-component BH model. Here,
$d_{0,i}=\sum_{n_{1},n_{2}}\left(n^{2}_{1}\,\delta_{i,1}+n^{2}_{2}\,\delta_{i,2}\right)\left|c^{0}_{n_{1},n_{2}}\right|^{2}$
is the mean-field value of the square density, while

$F=\langle\hat{A}^{2}{\left({\bf
0}\right)}\rangle=\frac{1}{I}\sum_{\alpha}\sum_{\bf
k}\sum_{n_{1},n_{2}}\left|v_{\alpha,{\bf k},n_{1},n_{2}}\right|^{2}$ (75)

is the control parameter of the theory, see App. (A.3) for a detailed
derivation. The quantization protocol corrects the mean-field value of local
observables in a two-fold way. On the one hand, the second term of Eq. (74)
makes a positive contribution due to quantum fluctuations; on the other hand,
the quantity $F$, which measures the magnitude of quantum fluctuations,
renormalizes the mean-field value $d_{0,i}$.

Next, we calculate the two-body correlation function between different species
$G^{(2)}_{12}{\left({\bf r=0}\right)}=\langle\hat{n}_{1}{\left({\bf
r=0}\right)}\,\hat{n}_{2}{\left({\bf 0}\right)}\rangle$, which is given by

$G^{(2)}_{12}{\left({\bf 0}\right)}=\langle\hat{\mathcal{D}}_{12}{\left({\bf
0}\right)}\rangle\,,$ (76)

where we have defined the composite density operator

$\hat{\mathcal{D}}_{12}{\left({\bf
r}\right)}=\sum_{n_{1},n_{2}}n_{1}\,n_{2}\,\hat{c}^{\dagger}_{n_{1},n_{2}}({\bf
r})\,\hat{c}_{n_{1},n_{2}}({\bf r})\,.$ (77)

We find the result

$G^{(2)}_{12}{\left({\bf
0}\right)}=\left(1-F\right)d_{0,12}+\frac{1}{I}\sum_{\alpha}\sum_{\bf
k}\sum_{n_{1},n_{2}}n_{1}\,n_{2}\left|v_{\alpha,{\bf
k},n_{1},n_{2}}\right|^{2}\,,$ (78)

where
$d_{0,12}=\sum_{n_{1},n_{2}}n_{1}\,n_{2}\left|c^{0}_{n_{1},n_{2}}\right|^{2}$
is the mean-field prediction.

Having outlined the form of single-species and pair correlations, the on-site
two-body correlation functions for the total density and spin channels can be
obtained directly from

$G^{(2)}_{d/s}{\left({\bf 0}\right)}=\langle\left[\hat{n}_{1}({\bf
0})\pm\hat{n}_{2}({\bf 0})\right]^{2}\rangle=G^{(2)}_{1}{\left({\bf
0}\right)}+G^{(2)}_{2}{\left({\bf 0}\right)}\pm 2\,G^{(2)}_{12}{\left({\bf
0}\right)}.$ (79)

For convenience, we analyze the normalized density and spin variances
$g^{(2)}_{d/s}{\left({\bf 0}\right)}$. Estimating the variances amounts to
shifting $G^{(2)}_{d/s}{\left({\bf 0}\right)}$ by $n^{2}_{d/s}$; additionally,
because $n_{s}=0$ in our model, we normalize the correlation functions by the
squared mean density $n^{2}_{d}$ to produce

$g^{(2)}_{d/s}{\left({\bf 0}\right)}=\frac{G^{(2)}_{d/s}{\left({\bf
0}\right)}-n^{2}_{d/s}}{n^{2}_{d}}\,.$ (80)

Our calculations for the density correlation $g^{(2)}_{d}({\bf 0})$ across the
various phase transitions of the model are shown in Fig. 15. In the SF region,
we observe the typical antibunching $g^{(2)}_{d}({\bf 0})<1$ of local
fluctuations due to the onsite interaction $U$. Moving towards the MI phase
[panels (a)-(b)], the qualitative features of density correlations are
strongly reminiscent of the behavior of a single-component BH model [35],
except for the the first-order MI-to-SF transition, at which the antibunching
factor shows a discontinuity. In particular, in the MI region, where mean-
field theory would predict $g^{(2)}_{d}({\bf 0})=0$, the QGW approach is able
to account for the virtual excitation of doublon-hole pairs, which leads to
the scaling $g^{(2)}_{i}({\bf 0})\propto\left(J/U\right)^{2}$ at low $J$, in
excellent agreement with perturbative calculations in the strongly-interacting
limit [31, 32]. Remarkably, we observe that the CFSF phase [panel (c)] shares
the same properties of the MI state in the density channel. This can be
understood as a consequence of the similarity between the spectral structure
of the two phases. In the CFSF phase, density fluctuations build on the
lowest-lying, gapped particle-hole excitations, which have a strong density
character exactly as their counterparts in the MI phase. By contrast,
$g^{(2)}_{d}({\bf 0})$ exhibits a quite different behavior across the PSF-to-
SF transition. Instead of being suppressed, sizeable density fluctuations
survive in the whole PSF region and saturate to a finite value at low $J/U$.
This result clearly agrees with the physical picture of the PSF phase [panel
(d)], where the formation of local pairs is explicitly favoured, differently
from the CFSF state, in which particles belonging to different species repel
each other. However, we believe that the independence of $g^{(2)}_{d}({\bf
0})$ on $J/U$ is a by-product of the inability of our theory to capture
density excitations in the PSF phase.

A complementary view on two-body correlations is provided by the spin
fluctuations $g^{(2)}_{s}({\bf 0})$, which are reported in Fig. 16 along the
same phase transition paths of Fig. 15. Across the MI-to-SF transition [panels
(a)-(b)], spin correlations mimic the behavior of density fluctuations,
meaning that the Mott interactions freezes equally both density and spin
degrees of freedom. This is not the case for the paired phases, where the
density and spin channels decouple. On the one hand, despite its Mott-like
character in the density channel, the CFSF phase [panel (c)] is characterized
by a large $g^{(2)}_{s}({\bf 0})$, indicating that counterflow superfluidity
is linked to the creation of local magnetic moments with large spin
fluctuations, such as in a spinful Mott state [61, 62]. Notice that, however,
$g^{(2)}_{s}({\bf 0})<1$ for finite values of $J/U$, signalling that the local
moments possess a finite stiffness due to the particle-hole excitations of the
CFSF phase. On the other hand, we find that $g^{(2)}_{s}({\bf 0})$ is strongly
suppressed in the PSF phase [panel (d)], where the spin degrees of freedom
interact repulsively, in analogy again with the physics of the MI-to-SF
transition [5].

##  5 Conclusions and Outlook

In this work, we have studied the properties of quantum fluctuations in the
two-component BH model at zero temperature, for both repulsive and attractive
interspecies interactions. Expanding on the mean-field ground-state phase
diagram first reported in Ref. [32], we have analyzed the band structure of
the model throughout the entire phase diagram with particular attention paid
to its quantum phase transitions and the limitations inherent in the local
nature of the Gutzwiller ansatz for the many-body wave function. Through a
canonical quantization of the variational parameters of the ansatz, we have
generalized to bosonic mixtures the single-component QGW theory of Ref. [35],
which permits a study of the quantum corrections to obsevables in a
systematic, order-by-order fashion in the spirit of the Bogoliubov theory of
dilute weakly-interacting gases. Importantly, we have shown how the QGW method
provides a comprehensive description of both local and non-local correlations
across the whole phase diagram of the model and, in particular, across its
quantum critical regimes.

We have first illustrated that the formation of pair correlations can be
directly connected to the analytical behavior of the compressibility and spin
susceptibility, which reflect the strength of critical fluctuations upon
reaching the PSF and CFSF transitions respectively. Notably, by the
application of spectral sum rules [33] together with the fluctuation-
dissipation theorem, we are able to relate the physical picture of the
response functions to the sound velocities of the Goldstone modes and the
static structure factors of the SF state. These results indicate that
experimental probes of strongly-correlated paired phases can yield direct
information on the spectral properties of the system in quantum critical
regimes.

Within the QGW approach, we have also studied the properties of superfluid
transport, finding an interspecies collisionless drag whose origin is due
solely to quantum fluctuations. In this respect, we have compared
quantitatively to QMC predictions [25] within the SF phase over a range of
interspecies couplings, and study the matrix of superfluid components across
the various phase transitions. In particular, we have found that the drag is
saturated in the vicinity of the PSF and CFSF phases, where strong pair
correlations prevail over single-particle coherence. Moreover, we have offered
a clear interpretation of the superfluid components in terms of multi-mode
scattering processes involving the collective excitations of the system,
including not only the density and spin Goldstone modes, but also the Higgs
modes and gapped quasiparticle excitations appearing at strong interactions.

Finally, we have shown that the QGW theory gives an accurate account of the
role of quantum fluctuations in few-body correlations throughout the phase
diagram of the system. In particular, we have demonstrated how the critical
behaviors of the one-body coherence function are analogous to those found for
MI-SF transition in the single-component BH model in Ref. [35]. Remarkably, at
the two-body level, we have also found that the CFSF/PSF phase transitions
strongly mirror the MI transition physics in the density/spin channels,
respectively. Throughout our analysis, we have highlighted how quantum
correlations closely link with the character of collective modes in distinct
interaction regimes.

The QGW approach developed in this work proves to be a flexible, semi-
analytical method to study the rich phase diagram of bosonic mixtures in an
optical lattice, including the effects of quantum fluctuations while remaining
computationally inexpensive. Although only the zero-temperature case is
considered in this work, we also note that finite temperature effects can be
incorporated in a straightforward manner. In particular, this method is
expected to be capable of describing also non-linear interactions between
quasiparticles, for instance in Landau/Beliaev-type decay processes [63, 64].

The generalization of the present theory to multi-component Bose-Fermi and
Bose-Bose mixtures and to different trapping geometries, where novel types of
superfluid drag are predicted [28, 65], poses intriguing problems for future
research. Additionally, an improvement in the description of non-local
fluctuations, crucial for instance to introduce hopping-induced correlations
into the paired phases, appears possible within a cluster extension of the
Gutzwiller theory (c.f. Refs. [66, 67, 68, 69, 70, 71]). Such an extension is
required, for example, to describe long-range interactions or magnetic
ordering in supersolid phases where translational symmetry is spontaneously
broken (c.f. Refs. [5, 72, 73]). Relaxing the $\mathbb{Z}_{2}$-symmetry
constraint also opens exciting questions on the imbalanced two-component BH
model, whose fermionic counterpart is currently the subject of intense
interest in the ultracold community (c.f. Refs. [74, 75]).

From a more formal perspective, the systematic way in which quantum
fluctuations can be incorporated in the QGW theory raises interesting
curiosities about its diagrammatic representation in quantum field theory.
Conversely, such a connection might enable the translation of the
comparatively sizeable literature of diagrammatic techniques into the language
of QGW theory to diagnose possible issues in the method or introduce concepts
already well-known in that context (c.f. [76, 77]).

## Acknowledgements

The authors thank D. Contessi, K. Sellin, and E. Babaev for fruitful
discussions.

##### Funding information

This project has received financial support from Provincia Autonoma di Trento
and the Italian MIUR through the PRIN2017 project CEnTraL (Protocol Number
20172H2SC4).

##  Appendix A Quantum corrections within the QGW approach

In this appendix, we provide additional technical details on the structure of
quantum fluctuations within the QGW theory and outline the explicit derivation
of the quantum correlations incorporated by the quantization procedure.

###  A.1 Evaluation protocol for the observables

Here, we provide a brief recap of the QGW quantization protocol of Ref. [35]
for evaluating an expectation value of a given observable
$\hat{O}{\left[\hat{a}_{i,{\bf r}},\hat{a}^{\dagger}_{i,{\bf r}}\right]}$:

  1. 1.

Determine the average value of the observable
$\mathcal{O}{\left[c,c^{*}\right]}=\big{\langle}\Psi_{G}\big{|}\,\hat{O}\,\big{|}\Psi_{G}\big{\rangle}$
in terms of the $\mathbb{C}$-valued Gutzwiller variables.

  2. 2.

Define the operator
$\hat{\mathcal{O}}{\left[\hat{c},\hat{c}^{\dagger}\right]}$ by promoting the
Gutzwiller parameters in $\mathcal{O}{\left[c,c^{*}\right]}$ to the
corresponding operators without modifying their ordering.

  3. 3.

Expand the operator $\hat{\mathcal{O}}$ order by order in the fluctuations
$\delta\hat{c}_{n_{1},n_{2}}$ and $\delta\hat{c}^{\dagger}_{n_{1},n_{2}}$,
including the contribution due to the normalization operator $\hat{A}$
properly. Using Eq. (41), this step provides a convenient decoding of the
original microscopic operator $\hat{O}$ in terms of the quantized modes of the
theory.

  4. 4.

Taking advantage of the quadratic character of the QGW Hamiltonian (42),
invoke Wick’s theorem to compute the expectation value of products of the
quasiparticle operators $\hat{b}_{\alpha,\mathbf{k}}$ on Gaussian states –
such as ground or thermal states obtained from $H^{\left(2\right)}_{QGW}$.

In the following, we apply the above quantization procedure to the single and
two-particle observables addressed in this work and clarify the order-by-order
derivation of the quantum corrections accessed by our quantum theory.

###  A.2 First-order fluctuations

The spectral features of the first-order fluctuations can be estimated either
by using the linearized 2GE (2.1) as done in Ref. [32] or alternatively
employing the QGW protocol, which has the formal advantage of associating the
creation/annihilation of a well-defined collective mode with a given spectral
amplitude, as outlined in this subsection.

(Density) – Within the QGW formalism, the species-resolved density
$\hat{n}_{i}({\bf r})$ is mapped into the operator

$\hat{\mathcal{N}}_{i}({\bf
r})=\sum_{n_{1},n_{2}}\left(n_{1}\,\delta_{i,1}+n_{2}\,\delta_{i,2}\right)\hat{c}^{\dagger}_{n_{1},n_{2}}({\bf
r})\,\hat{c}_{n_{1},n_{2}}({\bf r})\,.$ (A.1)

Expanding the $\hat{c}$’s to lowest order in the fluctuations, one finds

$\hat{\mathcal{N}}_{i}({\bf r})\approx
n_{0,i}+\delta_{1}\hat{\mathcal{N}}_{i}({\bf
r})=n_{0,i}+\sum_{n_{1},n_{2}}\left(n_{1}\,\delta_{i,1}+n_{2}\,\delta_{i,2}\right)c^{0}_{n_{1},n_{2}}\left[\delta\hat{c}_{n_{1},n_{2}}({\bf
r})+\delta\hat{c}^{\dagger}_{n_{1},n_{2}}({\bf r})\right]\,,$ (A.2)

where
$n_{0,i}=\sum_{n_{1},n_{2}}\left(n_{1}\,\delta_{i,1}+n_{2}\,\delta_{i,2}\right)\left|c^{0}_{n_{1},n_{2}}\right|^{2}$
is the mean-field density of the $i^{\text{th}}$ species. Using Eq. (41), the
first-order operator $\delta_{1}\hat{\mathcal{N}}_{i}({\bf r})$ can be
expressed straightforwardly in terms of the quasiparticle operators as

$\delta_{1}\hat{\mathcal{N}}_{i}({\bf
r})=\frac{1}{\sqrt{I}}\sum_{\alpha}\sum_{\bf k}N_{i,\alpha,{\bf
k}}\left(e^{i\,{\bf k}\cdot{\bf r}}\,\hat{b}_{\alpha,{\bf k}}+e^{-i\,{\bf
k}\cdot{\bf r}}\,\hat{b}^{\dagger}_{\alpha,{\bf k}}\right)$ (A.3)

where the spectral weight of single-mode density excitations $N_{i,\alpha,{\bf
k}}$ is given by Eq. (20).

(Square density) – According to the QGW mapping, the upgrade of the square
density operator reads $\hat{d}_{i}({\bf r})=\hat{n}^{2}_{i}({\bf
r})\to\hat{\mathcal{D}}_{i}({\bf r})$ as given by Eq. (73). Specifically, up
to the lowest order in the fluctuations, one obtains

$\hat{\mathcal{D}}_{i}({\bf
r})=\sum_{n_{1},n_{2}}\left(n^{2}_{1}\,\delta_{i,1}+n^{2}_{2}\,\delta_{i,2}\right)\hat{c}^{\dagger}_{n_{1},n_{2}}({\bf
r})\,\hat{c}_{n_{1},n_{2}}({\bf r})\approx
d_{0,i}+\delta_{1}\hat{\mathcal{D}}_{i}({\bf r})\,.$ (A.4)

where
$d_{0,i}=\sum_{n_{1},n_{2}}\left(n^{2}_{1}\,\delta_{i,1}+n^{2}_{2}\,\delta_{i,2}\right)\left|c^{0}_{n_{1},n_{2}}\right|^{2}$
is the mean-field square density of the $i^{\text{th}}$ species. The expansion
of the first-order operator $\delta_{1}\hat{\mathcal{D}}_{i}({\bf r})$ in
terms of the quasiparticle excitations is

$\displaystyle\delta_{1}\hat{\mathcal{D}}_{i}({\bf r})$
$\displaystyle=\sum_{n_{1},n_{2}}\left(n^{2}_{1}\,\delta_{i,1}+n^{2}_{2}\,\delta_{i,2}\right)c^{0}_{n_{1},n_{2}}\left[\delta\hat{c}_{n_{1},n_{2}}({\bf
r})+\delta\hat{c}^{\dagger}_{n_{1},n_{2}}({\bf r})\right],$ (A.5)
$\displaystyle=\frac{1}{\sqrt{I}}\sum_{\alpha}\sum_{\bf k}D_{i,\alpha,{\bf
k}}\left(e^{i\,{\bf k}\cdot{\bf r}}\,\hat{b}_{\alpha,{\bf k}}+e^{-i\,{\bf
k}\cdot{\bf r}}\,\hat{b}^{\dagger}_{\alpha,{\bf k}}\right),$

where we have defined the amplitude

$D_{i,\alpha,{\bf
k}}=\sum_{n_{1},n_{2}}n^{2}_{i}\,c^{0}_{n_{1},n_{2}}\left(u_{\alpha,{\bf
k},n_{1},n_{2}}+v_{\alpha,{\bf k},n_{1},n_{2}}\right)\,.$ (A.6)

(One-body boson field) – The quantization protocol of the QGW method upgrades
the order parameter of the one-body condensate to an effective boson field as
$\psi_{i}({\bf r})\to\hat{\psi}_{i}({\bf r})$. It is important to observe that
$\hat{\psi}_{i}({\bf r})$ is a macroscopic object (related to the coherence
properties of the system) and is distinct from the microscopic field operators
$\hat{a}_{i,{\bf r}}$ defined within second quantization at the beginning of
Sec. 2.1. To lowest order in the fluctuations, we find

$\hat{\psi}_{1}({\bf
r})=\sum_{n_{1},n_{2}}\sqrt{n_{1}}\,\hat{c}^{\dagger}_{n_{1}-1,n_{2}}({\bf
r})\,\hat{c}_{n_{1},n_{2}}({\bf
r})\approx\psi_{0,1}+\delta_{1}\hat{\psi}_{1}({\bf r})\,,$ (A.7)
$\hat{\psi}_{2}({\bf
r})=\sum_{n_{1},n_{2}}\sqrt{n_{2}}\,\hat{c}^{\dagger}_{n_{1},n_{2}-1}({\bf
r})\,\hat{c}_{n_{1},n_{2}}({\bf
r})\approx\psi_{0,2}+\delta_{1}\hat{\psi}_{2}({\bf r})\,,$ (A.8)

where $\psi_{0,i}$ is the one-body order parameter of the $i^{\text{th}}$
species at the mean-field level – see Eq. (4) – and

$\delta_{1}\hat{\psi}_{i}({\bf r})=\frac{1}{\sqrt{I}}\sum_{\alpha}\sum_{\bf
k}\left[U_{i,\alpha,{\bf k}}\,e^{i\,{\bf k}\cdot{\bf r}}\,\hat{b}_{\alpha,{\bf
k}}+V^{*}_{i,\alpha,{\bf k}}\,e^{-i\,{\bf k}\cdot{\bf
r}}\,\hat{b}^{\dagger}_{\alpha,{\bf k}}\right]\,.$ (A.9)

The previous result clearly shows that the spectral amplitudes
$U_{i,\alpha,{\bf k}}$ $\left(V_{i,\alpha,{\bf k}}\right)$, introduced in Eqs.
(18), quantify the particle (hole) character of the excitation
$\left(\alpha,{\bf k}\right)$ for the $i^{\text{th}}$ component of the system.
When dealing with the calculation of expectation values in momentum space, it
is often convenient to replace the microscopic bosonic fields by the
corresponding expression in terms of quasiparticle operators as

$\hat{a}_{i,{\bf k}}\longrightarrow\sqrt{I}\,\psi_{0,i}\,\delta_{\bf
k,0}+\sum_{\alpha}\left(U_{i,\alpha,{\bf k}}\,\hat{b}_{\alpha,{\bf
k}}+V^{*}_{i,\alpha,{\bf k}}\,\hat{b}^{\dagger}_{\alpha,{\bf k}}\right)\,.$
(A.10)

The structure of the particle (hole) amplitudes $U_{i,\alpha,{\bf k}}$
$\left(V_{i,\alpha,{\bf k}}\right)$ is strictly related to the bosonic
statistics of the collective excitations, since

$\left[\delta_{1}\hat{\psi}_{i}({\bf
r}),\delta_{1}\hat{\psi}^{\dagger}_{i}({\bf s})\right]=\frac{1}{I}\sum_{\bf
k}e^{i\,{\bf k}\cdot({\bf r}-{\bf
s})}\sum_{\alpha}\left(\left|U_{i,\alpha,{\bf
k}}\right|^{2}-\left|V_{i,\alpha,{\bf k}}\right|^{2}\right)=\delta_{\bf
r,s}\,,$ (A.11)

where last equality is verified provided that
$\sum_{\alpha}\left(\left|U_{i,\alpha,{\bf
k}}\right|^{2}-\left|V_{i,\alpha,{\bf k}}\right|^{2}\right)=1$. In practice,
we find that the previous identity is numerically satisfied if a sufficiently
large number of excitation branches are included into the $\alpha$-summation.

(Pair/antipair field) – The quantized pair and antipair fields read

$\hat{\psi}_{\mathrm{P}}({\bf
r})=\sum_{n_{1},n_{2}}\sqrt{n_{1}\,n_{2}}\,\hat{c}^{\dagger}_{n_{1}-1,n_{2}-1}({\bf
r})\,\hat{c}_{n_{1},n_{2}}({\bf r})-\hat{\psi}_{1}({\bf
r})\,\hat{\psi}_{2}({\bf r})\,,$ (A.12) $\hat{\psi}_{\mathrm{C}}({\bf
r})=\sum_{n_{1},n_{2}}\sqrt{n_{1}(n_{2}+1)}\,\hat{c}^{\dagger}_{n_{1}-1,n_{2}+1}({\bf
r})\,\hat{c}_{n_{1},n_{2}}({\bf r})-\hat{\psi}_{1}({\bf
r})\,\hat{\psi}^{\dagger}_{2}({\bf r})\,,$ (A.13)

Including first-order fluctuations only, the pairing operators can be expanded
as

$\hat{\psi}_{\mathrm{P}}({\bf
r})\approx\psi_{0,\mathrm{P}}+\delta_{1}\hat{\psi}_{\mathrm{P}}({\bf r})\,,$
(A.14) $\hat{\psi}_{\mathrm{C}}({\bf
r})\approx\psi_{0,\mathrm{C}}+\delta_{1}\hat{\psi}_{\mathrm{C}}({\bf r}),$
(A.15)

where the mean-field quantities $\psi_{0,\mathrm{P}}$ and
$\psi_{0,\mathrm{C}}$ are respectively the pair and antipair order parameters
given by Eqs. (5) and

$\displaystyle\delta_{1}\hat{\psi}_{\mathrm{P}}({\bf r})=$
$\displaystyle\sum_{n_{1},n_{2}}\sqrt{n_{1}\,n_{2}}\left[c^{0}_{n_{1},n_{2}}\,\delta\hat{c}^{\dagger}_{n_{1}-1,n_{2}-1}({\bf
r})+c^{0}_{n_{1}-1,n_{2}-1}\,\delta\hat{c}_{n_{1},n_{2}}({\bf r})\right]$
(A.16) $\displaystyle-\psi_{0,2}\,\delta_{1}\hat{\psi}_{1}({\bf
r})-\psi_{0,1}\,\delta_{1}\hat{\psi}_{2}({\bf r})\,,$
$\displaystyle\delta_{1}\hat{\psi}_{\mathrm{C}}({\bf r})=$
$\displaystyle\sum_{n_{1},n_{2}}\sqrt{n_{1}(n_{2}+1)}\left[c^{0}_{n_{1},n_{2}}\,\delta\hat{c}^{\dagger}_{n_{1}-1,n_{2}+1}({\bf
r})+c^{0}_{n_{1}-1,n_{2}+1}\,\delta\hat{c}_{n_{1},n_{2}}({\bf r})\right]$
(A.17) $\displaystyle-\psi_{0,2}\,\delta_{1}\hat{\psi}_{1}({\bf
r})-\psi_{0,1}\,\delta_{1}\hat{\psi}^{\dagger}_{2}({\bf r})\,.$

Recasting the $\delta\hat{c}$’s in the quasiparticle basis, we can rewrite the
first-order expansions (A.16)-(A.17) in a compact, suggestive form

$\displaystyle\delta_{1}\hat{\psi}_{\mathrm{P/C}}({\bf r})=$
$\displaystyle\frac{1}{\sqrt{I}}\sum_{\alpha}\sum_{\bf
k}\left[U_{\mathrm{P/C},\alpha,{\bf k}}\,e^{i\,{\bf k}\cdot{\bf
r}}\,\hat{b}_{\alpha,{\bf k}}+V^{*}_{\mathrm{P/C},\alpha,{\bf k}}\,e^{-i\,{\bf
k}\cdot{\bf r}}\,\hat{b}^{\dagger}_{\alpha,{\bf k}}\right]$ (A.18)
$\displaystyle-\psi_{0,2}\,\delta_{1}\hat{\psi}_{1}({\bf
r})-\psi_{0,1}\,\delta_{1}\hat{\psi}^{\dagger}_{2}({\bf r})\,,$

where the particle-hole amplitudes $U_{\mathrm{P/C},\alpha,{\bf k}}$ and
$V_{\mathrm{P/C},\alpha,{\bf k}}$ are explicitly given by

$U_{\mathrm{P},\alpha,{\bf
k}}=\sum_{n_{1},n_{2}}\sqrt{n_{1}\,n_{2}}\left(c^{0}_{n_{1}-1,n_{2}-1}\,u_{\alpha,{\bf
k},n_{1},n_{2}}+c^{0}_{n_{1},n_{2}}\,v_{\alpha,{\bf