text:
\section{Introduction}
Superconductivity is perhaps the most spectacular example of collective behavior of electrons.
The microscopic BCS \cite{bardeen1957theory,bardeen1957microscopic} theory of the phenomena explains superconductivity
in terms of bound (Cooper) electron pairs due to attractive interactions, mediated in the best understood case by
phonons. Condensation of Cooper pairs gives rise to a non-Fermi-liquid state with spontaneously broken $U(1)$ symmetry
characterized by a homogenous complex order parameter $\Delta$.  One experimental manifestation of superconductivity is
a non-Ohmic junction resistance with a normal metal.
At subgap values of the bias voltage on the interface, the only way for a single electron to penetrate inside a
superconductor is to form a Cooper pair. Therefore, by charge conservation, transmission through the interface should be accompanied by reflection of a hole, a process known as Andreev \cite{andreev1965thermal} reflection.  Studies of
Andreev reflection can be used to probe \cite{bradleySuperfluidAndreev, beenakkerColloquim} a state of interest and continue to gather both theoretical
\cite{de1995andreev,diez2012andreev,beenakker2006specular,nilsson2008splitting,stanescu2011majorana,kundu2013transport,mazin1999degreeofspinpol,mazin2001probing} and experimental \cite{sasaki2011topological,kastalsky1991observation,pothier1994flux,park2005andreev} attention.

In recent years experimentalists have uncovered examples of both electron-electron (Cooper) and
electron-hole (excitonic) pairing in atomically thin two-dimensional materials\cite{xi2016ising, li2018nontrivial,hsu2020inversion,eisensteinQHE, simonet2017anomalous,tutucEnhancement}.
In general Cooper pairing is driven by
attractive effective interactions, and electron-hole pair formation by repulsive interactions between electrons
in different bands.  Exciton condensates \cite{keldysh1968collective, hanamura1977condensation, littlewood2004models}
are coherent states of electrons and holes bound into pairs by the Coulomb interaction, and are
described by a mean-field theory that is identical
to BCS theory apart from a particle-hole transformation.
In bulk materials the concept of exciton condensation is partly ambiguous since the ordered states are difficult to distinguish from density-wave or nematic states \cite{fradkin2010nematic}.
Although exciton condensates were predicted theoretically more than 50 years ago,
their experimental identification in bulk materials has suffered as a consequence,
although progress has been achieved recently\cite{kogar2017signatures,mazuz2019dynamical,fogler2014high, high2012condensation}.
Experimental advances with two-dimensional materials have solved the ambiguity problem
by making it possible to prepare devices with strong interactions between subsystems
located in different layers that separately have nearly perfectly conserved particle number
Interlayer exciton condensation breaks the independent gauge invariance of the individual layers,
{\it i.e.} when the exciton order parameter $x$ is nonzero, there is spontaneous phase coherence
between layers and we no longer can perform gauge transformation independently.   There are still only a few
examples of convincing experimental evidence for exciton condensation in systems without
a strong external magnetic field, but this situation seems likely to change.

Recently, superconductivity has been observed in magic angle twisted bilayer
graphene (tBG) \cite{bistritzer2011moire, cao2018unconventional}, a two-dimensional system with
a moir\'e superlattice that yields extremely flat conduction and valence bands.
The strongest superconductivity appears when the
valence band is doped away from half filling.  Although several models \cite{wu2018theory, lian2019twisted, kozii2019nematic, chichinadze2020nematic, khalaf2020charged} for its superconductivity have been proposed, the mechanism remains unknown.
Twisted bilayers are the simplest examples of a rich variety of graphene multilayer moir\'e superlattice systems that
have received experimental attention recently \cite{cao2018correlated,lu2019superconductors,yankowitz2019tuning,park2021tunable}. The motivation for this paper is to
consider the possibility of realizing, in this flexible family of strongly-correlated electron systems, an exotic state in which interlayer exciton condensation and superconductivity
occur simultaneously.  Our target system is two twisted bilayers separated by
a hexagonal boron nitride tunnel barrier.  The superconductivity of the individual
twisted bilayers is already established.  Here we ask the following questions: i) can interlayer
exciton condensation \cite{rickhaus2020density} occur in principle between two-dimensional electron systems that are
superconducting ii) if so, how would such a state be most unambiguously detected?  So far
research on the coexistence of two phases has been scarce \cite{sager2020potential} partially because of the lack of a physical system with a potential to have both.
In section II  we describe the model we study and relationships between its order parameters.
In section III we describe the phases we have identified and phase diagrams as a function of the strength of the
interlayer repuslive interaction strength $g_x$ and the intralayer attractive interaction strength
$g_s$.  In the section IV we propose an Andreev drag measurement which can be used to
identify states with both types of order, taking advantage \cite{su2008make} of the possibility of
separate contacting to individual layers.
This method of probing bilayer exciton condensates has been
very successfully exploited in quantum Hall excitonic superfluids \cite{liu2017quantum,eisenstein2014exciton,yoon2010interlayer,li2017excitonic,li2019pairing}
Finally, in section V we discuss prospects for observing these
states and the importance of deviations from the highly symmetric point that we considered in our explicit calculations.
\iffalse
There are two main scenarios for spatially indirect exciton condensation.
In the first, a system is symmetric under the combination of layer interchange and
an electron-hole transformation. In this case, an electron-hole transformation in one
of the layers maps exciton condensates to superconductors that pair electrons in different layers and
occurs for arbitrarily weak interactions.
Consequently, charge transmission at low voltage occurs via a process
similar to Andreev reflection: to transfer an electron through an interface in one of the layers,
an electron should be reflected from the bottom layer.

In a scenario with two identical layers with the same Fermi level, exciton condensation is similar to itinerant electron
ferromagnetism with the spin degree-of-freedom replaced by a layer degree of freedom.  The phase will appear if the interaction $\lambda_x$ between the electrons and holes in two different layers is strong enough $g_x\equiv \lambda_x \nu(\epsilon_F)\geq 1$, where $\nu(\epsilon_F)$ is the density of states on the Fermi level.
 In this case, the eigenstates of the mean-field Hamiltonian have layer spinors that are coherent combinations of the
 individual layers, for example symmetric and antisymmetric.
 When viewed in this basis, transmission through the interface will be similar to free electron reflection from a
 potential.
  \fi

\section{Model}
We consider the simplest possible model that can have both interlayer exciton condensation and s-wave superconductivity.
The model assumes that  both the attractive intralayer and repulsive interlayer
interactions are momentum-independent, includes a layer degree of freedom ($l=t,b$),
and assumes either valley or spin singlet superconductivity
that is described by making a particle-hole transformation for one spin or valley in each layer,
The  Hamiltonian
\begin{multline}\label{model:full-hamiltonian}
	H
	=
	\sum_{pl\sigma}\xi_{pl}c^{\dagger}_{l\sigma}(p)c_{l\sigma}(p)
	\\+
	\frac{\lambda_s}{S}\sum_{pp'q}
	c^{\dagger}_{l\uparrow}(p+q)
	c^{\dagger}_{l\downarrow}(p'-q)
	c_{l\downarrow}(p')
	c_{l\uparrow}(p)
	\\+
	\frac{\lambda_x}{S}
	\sum_{pp'q}
	c^{\dagger}_{t\sigma}(p+q)
	c^{\dagger}_{b\sigma'}(p'-q)
	c_{b\sigma'}(p')
	c_{t\sigma}(p),
\end{multline}
where $c^\dagger_{l\sigma}$ $(c_{l\sigma})$ creates (annihilates) a
fermion in layer $l=t,b$ with pseudospin $\sigma=\uparrow,\downarrow$, and $S$ is the area of the sample.
In this work, we will focus on the case of layer-independent bands,
$\xi_{pl}\equiv \xi_p$, and reserve discussion of deviations from this symmetric limit to the end of the paper.
(As we explain in Section V coexistence does not occur for $\xi_{pt}\equiv -\xi_{pb}$.)
The pseudospin label refers to real spin in the case of spin-singlet superconductors and to
valley in the case of valley-singlet superconductors.  Below we refer to this degree of freedom as spin for simplicity.

Below we characterize the ordered states by their mean-field Hamiltonians,
and allow an exciton order parameter for each spin species,
\begin{equation}
	x_\sigma=\frac{\lambda_x}{S}\sum_p\<c_{t\sigma}(p)c_{b\sigma}^\dagger(p)\>,
\end{equation}
and superconducting order parameters both within and between layers,
\begin{equation}
	\Delta^{ll'}_{\sigma\sigma'}=-\frac{\lambda_{ll'}}{S}\sum_p\<c_{l\sigma}(p)c_{l\sigma'}(-p)\>,
\end{equation}
where $\lambda_{tt}=\lambda_{bb}=\lambda_s$ and $\lambda_{tb}=\lambda_{bt}=\lambda_x$.
The ordered state is therefore
characterized by six complex self-energies that vanish by symmetry if no symmetries are broken.

In a single-band superconductors, the pair potential can be chosen to be real because of global gauge invariance.
It follows that one real number, the gap, fully characterizes the superconducting phase.
In the present two-band system, we start with four complex pair potentials.
If we also allow excitonic particle-hole pairing for both spins, the ordered state is characterized by
12 real numbers.  Since only three phases can be chosen at will by exploiting the conservation of the number of particles of particles of
each spin in each layer, the energy can depend not only on absolute values of the order parameters,
but also on their phases, or, more specifically, on gauge invariant combinations of such phases.

 In the absence of interlayer superconducting coupling, there is only one gauge invariant combination of the phases: $\psi_t-\psi_b-\psi_{\uparrow}-\psi_{\downarrow}$, where $\psi_{\sigma=\uparrow,\downarrow}$ is the phase of the corresponding exciton pair potential and $\psi_{l=t,b}$ is the phase of the intralayer superconducting pairing potential. This combination changes sign  under  flip of layer. If we assume that the ground state preserves spin and layer invariance,
 we conclude that $\psi_t-\psi_b-\psi_{\uparrow}-\psi_{\downarrow}=\pi n$, where $n$ is an integer.
 In the gauge $\psi_b=0$, $\psi_{\uparrow}=0$, $\psi_{\downarrow}=0$, an even $n$ will correspond to a state with two supercondutving gaps having the same sign, while an odd $n$ will correspond to a state with superconducting gap that changes sign under the layer flip.

To extend this analysis to the case with interlayer pairing,
we identify four gauge invariant quantities that do not transform to themselves under layer or spin inversion:
$x_\uparrow \Delta^{bt}_{\uparrow\downarrow} {\Delta^{tt}_{\uparrow\downarrow}}^*$, $x_\downarrow \Delta^{bt}_{\downarrow\uparrow}{\Delta^{tt}_{\downarrow\uparrow}}^*$, $x^*_\uparrow \Delta^{tb}_{\uparrow\downarrow}{\Delta^{bb}_{\uparrow\downarrow}}^*$, and $x^*_\downarrow \Delta^{tb}_{\downarrow\uparrow}{\Delta^{bb}_{\downarrow\uparrow}}^*$
and require that their values do not change under these transformations. \footnote{One way to represent these quantities is to imagine space with 4 points corresponding to each kind of fermion. Then any gauge invariant quantity will correspond to a closed trajectory between these points.}
This condition restricts the absolute values of all self-energies
$|\Delta^t|=|\Delta^b|\equiv\Delta_d$, $|x_{\uparrow}|= |x_{\downarrow}|\equiv x$, $|\Delta^{tb}_{\uparrow\downarrow}|=|\Delta^{bt}_{\uparrow\downarrow}|\equiv \Delta_i$ as well as two gauge invariant combinations of phases
\begin{align}\label{phase-cond:1}
\psi_1-\psi_2=\psi_\downarrow-\psi_\uparrow-2\pi m,\\\label{phase-cond:2} \psi_t-\psi_b=\psi_\uparrow+\psi_\downarrow-2\pi k,
\end{align}
where $\psi_1$ is the phase of $\Delta^{bt}_{\uparrow\downarrow}$ and $\psi_2$ is the phase of  $\Delta^{tb}_{\uparrow\downarrow}$.
($d$, and $i$ in $\Delta_d$ and $\Delta_i$ are
intended to suggest {\it direct} and {\it indirect} pairing.) The system has one more independent gauge invariant combination of phases $\psi_t+\psi_b-\psi_1-\psi_2$. Unlike those we considered previously, this combination transforms to itself under
layer or spin flip, and therefore should be determined by energy minimization.

Using global gauge invariance related to conserved
particle numbers for each spin and valley $\psi_1$, $\psi_\uparrow$ and  $\psi_\downarrow$ can be chosen to be 0.  Then because of \eqref{phase-cond:1} $\psi_2=2\pi m$ and $\psi_-\equiv (\psi_t-\psi_b)/2=\pi k$. In what follows, we use this gauge.
This leaves us with
mean-field Hamiltonians that depend only on four parameters, three absolute values and a single phase,
whose values can be determined by solving self-consistent field equations.

\iffalse
leaves us with 6 real parameters, 3 of which can be gauged away because of number conservation. We conclude then it suffices to have 3 real parameters to describe the mean field state.

In single-band superconductors, the pairing order parameter can be chosen to be real because of global gauge invariance. Then there will be only one real number that characterizes the superconducting phase. In this simple case, one way to see this is to recall that the $U(2)$ Bogolyubov matrix has 4 independent real parameters. Among those, there is a single angle that characterizes $O(2)$ rotation as well as three complex phases. Because the Hamiltonian is gauge invariant both in new and the old basis we can always choose the gauge where all 3 phases are zero.

In the present system of two layers the situation is different and there are phases that cannot be gauged away.  In other words, the energy of the uniform ground state of the system will depend on combinations of order parameter phases.

The $U(4)$ Bogolyubov matrix has 16 independent real parameters. $6$ of them parametrize $O(4)$ rotation matrix. Among the remaining $10$ phases, $4+4$ can be made zeros and two remain. The energy can therefore depend on them. Now, if we assume that mirror (layer equivalence)  and spin-invariance is maintained, it follows that
$|\Delta^{t}_{\uparrow\downarrow}|=|\Delta^{b}_{\uparrow\downarrow}|$ and that $|x_{\uparrow}|=|x_{\downarrow}|$ as well as $|\Delta^{tb}_{\uparrow\downarrow}|=|\Delta^{bt}_{\uparrow\downarrow}|$. These 3 conditions will restrict reduce the number of independent rotations to $3$. We conclude then that 5 parameters (3 angles and 2 phases) fully characterize the mixed state.

The generation of phase masses is reminiscent to the case of the Josephson
junction, but in our case it has deeper consequences on the phase diagram itself.
Let's see in detail how global gauge transformations work in our system. In this work, we use a gauge where all but intralayer superconductivity order parameters are real and positive. In this gauge, $x$ gives mass to the difference of phases $2\psi_-=\psi_b-\psi_t$ of intralayer superconducting order parameters, and $\Delta_i$ gives mass to $2\psi_+=\psi_t+\psi_b$.
\fi

\section{Phase Diagram}
We saw previously that layer and spin symmetries restrict the value of $\psi_-$, but not
$\psi_+$. In this section, we allow both phases to take any value demonstrate explicitly that energies are minimized when spin/layer symmetries are respected: the four allowed
values of $\psi_-$ follow from the self-consistency equations.
The mean-field version of \eqref{model:full-hamiltonian} can be written as
\begin{equation}\label{phases:ham-matrix}
	H_{mf}=\sum_p \Psi^\dagger(p)
\mathcal{H}(p)
	\Psi(p),
\end{equation}
where $\Psi^\dagger(p)=\(c^\dagger_{t\uparrow},c_{t\downarrow},c^\dagger_{b\uparrow},c_{b\downarrow}\)$ is a vector in an extended Nambu space and the Hamiltonian matrix is
\begin{widetext}
\begin{equation}
  \mathcal{H}(p)=	\begin{pmatrix}
		\xi_p & \Delta_d e^{i(\psi_++\psi_-)} & x & \Delta_i \\
		\Delta_d e^{-i(\psi_++\psi_-)}&-\xi_p & \Delta_i & -x\\
		x & \Delta_i & \xi_p &\Delta_d e^{i(\psi_+-\psi_-)}\\
		\Delta_i & -x & \Delta_d e^{-i(\psi_+-\psi_-)} &-\xi_p
	\end{pmatrix}
\end{equation}.
\end{widetext}
The negative eigenvalues of the matrix have the form $\epsilon_{\pm}=-\(\eta^2\pm 2\alpha^2\)^{1/2}$, where $\eta^2=\xi^2+x^2+\Delta_i^2+\Delta^2_d$ and $\alpha^4=x^2 \Delta^2_d \sin^2(\psi_-)+2\Delta_d \Delta_i x\xi \cos(\psi_+)\cos(\psi_-)+\Delta_i^2\Delta_d^2\cos^2(\psi_+)+(x\xi)^2$.
Stable phases of the system are determined by minimization of the energy density.
At energy extrema,
\begin{equation} \label{phases:total-energy}
	E=\frac{2\Delta^2_d}{|\lambda_s|}-\frac{2\Delta_i^2}{\lambda_x}+\frac{2 x^2}{\lambda_x}+\sum_b\int d\xi \nu(\xi )\epsilon_{b}n_b(\xi),
\end{equation}
where the sum is over quasiparticle bands $b$, $\nu(\xi)$ is the quasiparticle density-of-states per spin and per layer,
and $n_b(\xi)$ is the Fermi occupation factor.
At zero temperature, only quasiparticle states with negative energies are filled.
\iffalse
 \begin{figure}[htb!]
	\includegraphics[width=\columnwidth]{Figures/energy_contour_vs_phi_fixed_orders_cust_cmap.png}
	\caption{Energy density \eqref{phases:total-energy} as a function of $\psi_-$ (x-axis) and $\psi_+$ (y-axis). Absolute values of all three order parameters fixed at some realistic values (on the main figure $\mu=0.1\Lambda$ $x=\mu$, $\Delta_d=0.01\Lambda$, $\Delta_i=\Delta_d/3$, and the density of states is the same as on Fig.\ref{fig:parallel-phase:underdoped} B), on the inset $\Delta_i=0$). The highest energy corresponds to red, the lowest energy corresponds to blue. We see that $\Delta_i$ lifts the energy degeneracy in $\psi_+$. Note that all minima are connected by $\psi_l\to \psi_l+2\pi n, \psi_l'\to\psi_l'+2\pi m$ for some integer $m$ and $n$ and therefore correspond to the same state. }
	\label{fig:parallel-phase:energy-vs-phi}
\end{figure}
\fi
Extrema of the energy are also extrema of Eq.~\eqref{phases:total-energy}, which we vary first
with respect to the phases $\psi_-$ and $\psi_+$.
We conclude that energy minima occur at extrema of $\alpha^2$, and that these
occur (independent of $\xi_p$)
when
\begin{equation}\label{phase-diagram:parallel-condition}
	\sin\(\psi_+\)=0 \text{ and } 	\sin\(\psi_-\)=0
	\end{equation}
	or
	\begin{equation}\label{phase-diagram:antiparallel-condition}
	\cos\(\psi_+\)=0 \text{ and } 	\cos\(\psi_-\)=0.
\end{equation}
The former condition is consistent with the conditions we derived for a mirror- and spin-symmetric state with $\Delta_i\neq0$.
We will call it a parallel phase.  The latter condition defines the antiparallel phase \cite{phase-condition-comment}.
We will see that the latter is consistent only with $\Delta_i=0$.
Note that each condition in fact corresponds to two different states. Namely, \eqref{phase-diagram:parallel-condition} can be satisfied with $(\psi_t,\psi_b)=(\pi,\pi)$, and $(0,0)$, that typically have different energies
when interlayer superconducting order is present.  Now we consider the antiparallel and parallel phases separately.

\subsection{Antiparallel phase}
We consider first extrema in which the phases $\psi_t$ and $\psi_b$ differ by $\pi$.
.
The main conclusion of this subsection is that for any sufficiently smooth density of states,
the energy of the antiparallel phase is always higher than that of a pure superconducting phase.
In other words this solution of the mean-field equations corresponds to a local energy maximum not
an energy minimum, and can therefore be discarded.

  \begin{figure*}[t]
	\includegraphics[width=18cm,height=5.5 cm]{Figures/order_parameters_lxstep0005_lsstep0005_random_den8nmin40.png}
	\caption{Order parameter dependence on dimensionless interaction constants $g_s=\nu \lambda_s$,
	and $g_x=\nu \lambda_x$ in the parallel phase,
	where $\nu$ is the parabolic band density-of-states.
	The exciton self-energy $x$
	Fermi energy is plotted in panel A.
	At small $|g_s|$ there is a phase transition at $g_x=1$, which is first order in our model, to a large
	$g_x$ state with $x \ne 0$.  The intralayer superconducting self-energy $\Delta_d$, plotted in
	panel B is always non zero.
	The interlayer superconducting self-energy $\Delta_i$ is plotted
	in panel C and acquires a finite value only when
	$x \ne 0$.  After the transition the purely intralayer superconducting phase
	to the mixed phase, $\Delta_d$ (B) tend to decrease. For large values of $|g_s|$ the excitonic transition is pushed to larger values of $g_x$. The inset in panel A shows the three self-energies as a function of $g_x$ at a fixed $g_s=-0.4$. These results were obtained for cutoff $\Lambda=10 \ meV$, $\epsilon_F=3 \ meV$, and a parabolic band model with a constant density of states. To treat both $f=0$ and $f=1$ minima on the equal footing, we allowed the interlayer superconducting gap $\Delta_i$ to take negative values. Self-energies are in units of meV. We minimized energy as a function of the Bogolyubov angles using the conjugate graduate method.}
	\label{fig:parallel-phase:underdoped}
\end{figure*}
In the antiparallel phase, variation of the energy density in Eq.~\ref{phases:total-energy} with respect to ${\Delta_i}$ combines with Eq.~\eqref{phase-diagram:antiparallel-condition} to yield the self-consistency equation
\begin{equation}
	\Delta_i
	=
	-\frac{\lambda_x \Delta_i}{4}\int\nu(\xi)\(\frac{1}{\epsilon_{+}}+\frac{1}{\epsilon_{-}}\).
\end{equation}
This equation clearly cannot be satisfied at any $\lambda_x>0$.  We conclude that
$\Delta_i=0$ in the antiparallel phase.
The quasiparticle band energies are therefore
\begin{equation}\label{antiparallel:band-energy}
	\pm \epsilon_{\pm}=\pm |\sqrt{\xi^2+\Delta_d^2}\pm x|\\
\end{equation}
The form of the band energies allows us to identify
the antiparallel phase as an exciton condensate (spontaneous interlayer phase coherent state)
formed on top of superconductors within each of the layers.
With ${\Delta_i}$ eliminated, the system of gap equations reduces to:
\begin{equation}\label{phase-diagram:antiparallel:x}
	x=
	\lambda_x
	\mathcal{N}
	(\sqrt{x^2-\Delta_d^2} \,),
\end{equation}
\begin{multline}\label{phase-diagram:antiparallel:delta}
	\Delta_d=-\frac{\lambda_s \Delta_d}{2}
	\int  \frac{\nu(\xi)}{\sqrt{\xi^2+\Delta_d^2}}+\\
	\frac{\lambda_s \Delta_d}{2}\int^{\sqrt{x^2-\Delta_d^2}}_{-\sqrt{x^2-\Delta_d^2}}\frac{\nu(\xi)}{\sqrt{\xi^2+\Delta_d^2}}.
\end{multline}
That is to say that the exchange self-energy is contributed by wavevectors that have different
occupation numbers for the two orthogonal quasiparticle layer spinors.
In Eq.~\eqref{phase-diagram:antiparallel:x}  $\mathcal{N}(\epsilon) \equiv 1/2\int_{-\epsilon}^{\epsilon} d\xi \ \nu(\xi)$.

\iffalse
A transition between two superconductors and the mixed phase becomes possible at the smallest value of the interlayer interaction strength $\lambda_x$ where the gap equations can be satisfied.  From equation \eqref{phase-diagram:antiparallel:x} it follows that at the transition point the value $\frac{x}{\mathcal{N}(\sqrt{x^2-\Delta_d^2})}$,  has the minimum possible value.  For any given $x$ and any continuous density of state it implies that $\Delta_d=0$, so that at a transition point the superconducting gap in each layer should close. It is clear then from \eqref{phase-diagram:antiparallel:delta} that at the same point $x$ should be nonzero to ensure the gap closing. We conclude then that the possible transition between two phases is of the first order.
\fi
A transition between pure superconductivity and an antiparallel
phase would occur at zero temperature if the latter phase were lower in energy.
We find that
\begin{multline}
\label{eq:energydiff}
	E_{m}-E_{sc}=
	\frac{2}{|\lambda_s|}\({\Delta_d}_m^2-{\Delta_d}_0^2\)
	\\-2\int_{-\mu}^\Lambda\nu(\xi)
	\(\sqrt{\xi^2+{\Delta_d}_m^2}-\sqrt{\xi^2+{\Delta_d}_0^2}\)
	\\+
	2\int^{\sqrt{x^2-{\Delta_d}^2_m}}_{-\sqrt{x^2-{\Delta_d}^2_m}}
	\nu(\xi)
	\sqrt{\xi^2+{\Delta_d}_m^2}
	-
	2x \mathcal{N}(\sqrt{x^2-{\Delta_d}^2_m})
	\leq
	0,
\end{multline}
 where we distinguish the superconducting gap ${\Delta_d}_m$ in the mixed phase from the superconducting gap ${\Delta_d}_0$ in the pure superconducting phase, and have used a self-consistency equation to reexpress the last term.
 Since the energy of the pure superconducting phase is minimized by ${\Delta_d}_0$,
 we can conclude that the sum of the first two lines in Eq.~\ref{eq:energydiff} is non-negative.
 Provided that the density of states does not significantly vary  on the scale of $\sqrt{x^2-{\Delta_d}^2_m}$,
 the third line can be approximated by
 \begin{multline}
 	2\nu_0{\Delta_d}_m^2\sinh^{-1}\(\frac{\sqrt{x^2-{\Delta_d}_m^2}}{{\Delta_d}_m}\),
 \end{multline}
where $\nu_0$ is the constant density-of-states, and the full energy difference is positive.
Hence, at least for a sufficiently smooth density of states,
the antiparallel phase will never be thermodynamically stable. We now analyze the parallel phase.

 \iffalse
 \begin{figure*}[t]
	\includegraphics[width=18cm]{Figures/order_parameters_lxstep01_lsstep002_mu5.png}
	\caption{Order parameter dependence on dimensionless interaction constants $g_s=\nu_{max}\lambda_s$, $g_x=\nu_{max}\lambda_x$ in the parallel phase for $\nu(\epsilon)=\nu_{max}\gamma^2/\(\gamma^2+(\xi-\xi_s)^2\)$with $\xi_s=0$ and $\gamma=1$.	Normalized exciton self-energy $x/\epsilon_F$ is plotted on A). For small absolute values of $g_s$, the exciton condensate appear at the point $g_x=1$ in accordance to the Stoner criterion. Interlayer superconductivity ${\Delta_i}$ is vanishingly since the Fermi energy is tuned to the peak in the density of states. In this calculation, $\Lambda=10 \ meV$, $\epsilon_F=5 \ meV$, $\gamma=1 \ meV$, $\alpha=1$, $\xi_s=0\ meV$. We performed a minimization of energy as a function of the Bogolyubov angles using conjugate graduate method.}
	\label{fig:parallel-phase:van-hove}
\end{figure*}
\fi
\subsection{Parallel phase}
The parallel phase is defined by the conditions $\sin(\psi_+)=0$, $\sin(\psi_-)=0$, so that $\psi_t$ and $\psi_b$ are equal to within a multiple of
$2\pi$.  This condition, however, leaves an ambiguity since the $\cos(\psi_+)\cos(\psi_-)$ factor
present in the dispersion relation can be either $+1$ or $-1$.
We distinguish these possibilities by introducing $f=0,1$ such that $\cos(\psi_+)\cos(\psi_-)=(-1)^f$.
The quasiparticle energies can be expressed as:
\begin{equation}\label{phase-diagram:parallel:dispersion}
	\epsilon_\pm=\sqrt{(\xi \pm x)^2+({\Delta_d} \pm (-1)^f {\Delta_i})^2}.
\end{equation}
We can interpret this phase as two superconductors: one with Cooper pairs formed out of symmetric combinations of layers, another -- from antisymmetric combinations.  Indeed, the mean-field Hamiltonian matrix \eqref{phases:ham-matrix} is block-diagonal in this basis. Note here that even though the transformation to the symmetric and antisymmetric states block-diagonalizes $\mathcal{H}(p)$, energy density $E$ \eqref{phases:total-energy} will still have terms that couple the two unless $|\lambda_s|=\lambda_x$.
Solutions with $f=0/1$ then correspond to phases with in which either symmetric or antisymmetric superconducting gaps are larger.
Variation of the full energy density \eqref{phases:total-energy} with respect to ${\Delta_i}$ yields
\begin{multline}\label{phase-diagram:parallel:energy-derivative:sigma}
	\frac{\delta E}{\delta {\Delta_i}}
	=
	-\frac{4{\Delta_i}}{\lambda_x}
	-
	{\Delta_i} \int
	\nu(\xi)
	\(
	\frac{1}{\epsilon_-}
	+
	\frac{1}{\epsilon_+}
	\)
	\\+
	(-1)^f
	\Delta_d \int
	\nu(\xi)
	\(
	\frac{1}{\epsilon_-}
	-
	\frac{1}{\epsilon_+}
	\)
\end{multline}
One conclusion we can immediately draw is that the ${\Delta_i}=0$ point is not an energy extremum
whenever both exciton condensates and superconductivity are present.
Note that the superconducting gap of the symmetric quasiparticles is $\Delta_d+ (-)^f{\Delta_i}$,
while for the antisymmetric quasiparticles is $\Delta_d- (-)^f {\Delta_i}$.
Because the exciton condensate breaks the symmetry between quasiparticle bands,
one should expect that $\Delta_d+{\Delta_i}\neq\Delta_d-{\Delta_i}$.
A complimentary explanation is that even though the interlayer interaction is repulsive,
the BCS instability present within each layer together with interlayer coherence
induces a response in the interlayer Cooper channel.

If we write the self-consistency equations in terms of
$\Delta_-\equiv\Delta_d-(-1)^f{\Delta_i}$, $\Delta_+\equiv \Delta_d}+(-1)^f{\Delta_i}$ and $x$ we obtain,
\begin{multline}\label{phase-diagram:parallel:delta+}
	\frac{1}{\lambda_-^2-\lambda_+^2}(\lambda_+\Delta_+-\lambda_-\Delta_-)
	\\=
	\frac{\Delta_+}{2}
	\int \frac{d\xi \nu(\xi)}{\sqrt{(\xi+x)^2+\Delta_+^2}},
	\end{multline}
\begin{multline}
	\label{phase-diagram:parallel:delta-}
\frac{1}{\lambda_-^2-\lambda_+^2}(\lambda_+\Delta_--\lambda_-\Delta_+)
	\\=
	\frac{\Delta_-}{2}
	\int \frac{d\xi \nu(\xi)}{\sqrt{(\xi-x)^2+\Delta_-^2}},
\end{multline}
\begin{multline}
	\label{phase-diagram:parallel:x}
		\frac{4}{\lambda_+-\lambda_-}x=\int \frac{d\xi\nu(\xi)(x-\xi)}{\sqrt{(\xi-x)^2+\Delta_-^2}}
	\\+
	\int \frac{d\xi\nu(\xi)(x+\xi)}{\sqrt{(\xi+x)^2+\Delta_+^2}}.
\end{multline}
Here $\lambda_+\equiv(\lambda_s+\lambda_x)/2$ and $\lambda_-\equiv(\lambda_s-\lambda_x)/2$. Note that even though the mean-field
Hamiltonian matrix is block diagonal, the energy density has a term of the form $\Delta_+\Delta_-$ because $\lambda_s\neq\lambda_x$.
We solve this system numerically for the constant density-of-states case.
The mixed state solution summarized in Fig.\ref{fig:parallel-phase:underdoped} exists for $\nu(0)g_x\geq 1$.
  \begin{figure}[htb!]
	\includegraphics[width=\columnwidth]{Figures/atb_curve_combined_potentials.png}
	\caption{Probabilities for an electron coming from the top lead at energy $\omega$ to be reflected as a hole/electron in the same layer ($A_{tt}$/$B_{tt}$), and as a hole/electron in the opposite layer ($A_{tb}$/$B_{tb}$).
	The quantity $A_{tb}-B_{tb}$, proportional to the differential conductance in the bottom layer,
	changes sign as function of $\omega$. At low energies, the conductance is dominated by electron-electron reflection.
	At energies around superconducting gaps, the conductance is dominated by Andreev reflection processes.
	The inset shows a schematic experimental device with voltage bias $V$ across the top layer.
	When mesoscopic effects are neglected, the top layer conductance of the device is $\sigma_{t}(V/2)/2$,
	where $\sigma_{t}(V)$ is the differential conductance between the top left lead and the top layer,
	and the transconductance driven by excitonic order is $\sigma_{b}(V/2)/2$, where
	$\sigma_{b}(V)$ is the conductance between the top left lead and the bottom layer.}
	\label{fig:exp:setup}
\end{figure}
  \begin{figure}[htb!]
	\includegraphics[width=\columnwidth]{Figures/sigmatb_dis_curve_plotDp0.1Dm0.15x0.8V4V.png}
	\caption{Influence of interface disorder on the difference between hole (Andreev process) and
	electron interlayer reflection coefficients $A_{tb}-B_{tb}$.
	Disorder is modeled by a delta-function potential at the interface $v_l\delta(z)$.
	On each subplot, the ratio of $\<v\>=(v_t+v_b)/2$ to the Fermi momentum is kept fixed. $\Delta v/p_F=(v_t-v_b)/(2p_F)$ is swept from -1.0 (red) to +0.8 (black) with step 0.1.}
	\label{fig:exp:disorder}
\end{figure}
In the numerical calculation, to treat both $f=0$ and $f=1$ minima on the equal footing, we allowed the interlayer superconducting gap $\Delta_i$ to take negative values. Then negative value of $\Delta_i$ corresponds to the $f=1$ minimum while positive values will correspond to the $f=0$ minimum. Our calculation shows that $f=1$ minimum tend to be more energetically favorable for larger values of intralayer attraction, which has a simple physical meaning: for $|\lambda_s|\approx \lambda_x$ the symmetric and the antisymmetric superconductors are effectively decoupled and thus larger gap correspond to a larger effective Fermi energy.
Along the first order phase boundary in $g_s$/$g_x$ space separating mixed and pure superconducting phases,
the two states have identical energies.  It follows that the differential energy changes
$dE_{sc}$ in the superconducting phase and $dE_{mix}$ in the mixed phase should be the same.
It follows that the slope of the phase boundary lines
\begin{equation}\label{phase-diagram:parallel:coexistence-curve}
	\frac{dg_x}{dg_s}
	=
	\frac{\frac{\d E_{sc}}{\d g_s}-\frac{\d E_{mix}}{\d g_s}}{\frac{\d E_{mix}}{\d g_x}}.
\end{equation}
Partial derivatives in the formula will be proportional to order parameters in each phase: $\d E_{sc,mix}/\d g_s=-2\Delta_{sc,mix}^2/g_s^2$, $\d E_{sc}/\d g_x=2(x^2-{\Delta_i}^2)/g_x^2\approx 2\mu^2/g_x^2$. Formation of the mixed state accompanied by the decrease in the superconducting gap at least for large enough $|g_s|$. Therefore, $\frac{dg_x}{dg_s}<0$ --  for larger attractive interaction $g_s$ transition will happen at a larger critical $g_x$ in agreement with Fig. \ref{fig:parallel-phase:underdoped}. \iffalse Assuming that $\Delta_{sc}^2-\Delta_{mix}^2=g_s^2 c$ with $c$ being weakly dependent on $g_s$ and $g_x$ at the interface, by integrating \label{phase-diagram:parallel:coexistence-curve} from point $(g_x,g_s)=(1,0)$ we get $g_x=(1+cg_s/\mu^2)^{-1}$ along the interface.\fi

\section{Non-local Andreev Reflection}

Having established that an exotic mixed state with both excitonic and Cooper pair order can exist,
we now discuss Andreev reflection experiments in separately contacted bilayers that
can be used to detect their presence. Consider metallic leads separately connected to the top and the bottom
layers of our system. An electron incident on the interface between the top lead and the top layer
has equal weight in symmetric and antisymmetric quasiparticle
channels $\ket{t}=(\ket{+}+\ket{-})2^{-1/2}$. Since the parallel state does not couple $\ket{+}$ with $\ket{-}$, we can
consider their transmission contributions separately, at least if disorder is neglected.
Incoming quasiparticles are reflected with probabilities $A_{tt}$ and $A_{tb}$ as holes,
and with probabilities $B_{tt}$ and $B_{tb}$ as electrons, in the top and bottom layers respectively.
The calculation of these reflection probabilities is discussed below.

If we apply $V$ voltage across the top part of the interface, the
top and bottom layer differential conductances will be connected to reflection coefficients in a normal fashion.
The total currents in the top and bottom layers are
\begin{equation}
I_{t/b} = e \nu_0 \int_{0}^{eV} d \omega \; \sigma_{t/b}(V)
\end{equation}
where
\begin{eqnarray}
 \sigma_t(V)&\equiv& dj_t/dV=e^2 v_f\nu_0 (1-B_{tt}(V)+A_{tt}(V)), \nonumber \\
 \sigma_b(V)&\equiv& dj_b/dV=e^2 v_f\nu_0 (A_{tb}(V)-B_{bt}(V)),
 \end{eqnarray}
$e$ is an electron charge, and $v_f\nu_0$ is the product of the Fermi velocity and the
density of states.  Note that the combination
 $A_{tb}-B_{tb}$ of cross-layer reflection probabilities can be measured experimentally.
 We derive expressions for the scatterring amplitudes in the Appendix using the Bogoliubov-de-Gennes equations. The amplitude to reflect an electron as a hole in either symmetric or antisymmetric states, at least in the limit $x\ll \mu$, around 1 at $\omega=0$ similarly to a single N-S junction. However, the amplitude to reflect an incoming electron as an electron is no longer negligible even in the $\Delta_d\ll \mu$ case. The electron reflection amplitude at $\omega=0$ is given by $b_\pm\approx\mp x/(2\mu\pm x)$.
 Because $b_+$ and $b_-$ have opposite signs, the probability to reflect an electron to the bottom layer $B_{tb}=|b_+-b_-|^2/4$ is substantilly higher than the probability of the corresponding Andreev process, which is almost zero.
For the the same reason electron reflection in the top layer is unimportant for small $\omega$ and instead the Andreev process dominates (Fig.\ref{fig:exp:setup}).
When energy is equal to the $+$ or $-$ gap and the gap $\Delta_\pm/\mu\ll 1$,
the corresponding momenta in the leads and device match, giving rise to peaks in Andreev reflection and minima in electron reflection.
As a result the Andreev process is dominant over electron reflection in the bottom layer.
Altogether, we expect $A_{tb}-B_{tb}$ to have a $\mu$-like shape, with negative values
around $\omega=0$ and positive peaks at the superconducting gaps (Fig.\ref{fig:exp:setup}).

\iffalse $A_{tb}-B_{tb}$ is a quantity one can measure experimentally.
 If we apply $V$ voltage across the top part of the interface, top and bottom differential conductance will be connected to reflection coefficients in a normal fashion: $\sigma_t(V)\equiv dj_t/dV=e^2v_f\nu_0 (1-B_{tt}(V)+A_{tt}(V))$, $\sigma_b(V)\equiv dj_b/dV=e^2v_f \nu_0 (A_{tb}(V)-B_{bt}(V))$, where $e$ is an electron charge, and $v_f\nu_0$ is a characteristic of a lead. As a result, we expect the differential conductance in the bottom layer to change sign several times as function of the applied voltage, with robust positive peaks at values corresponding to the superconducting gaps.
If we model disorder by a delta-function $v_l\delta(z)$ at the interface in each layer, $\<v\>=(v_t+v_b)/2$ will play a role identical to the disorder in a regular superconductor: by damping Andreev reflection at subgap energies and increasing electron-electron reflection in symmetric and antissymetric superconductors separately. A dimensionless quantity controlling its importance is $\<v\>/v_F$, where $v_F$ is the Fermi velocity in the lead. A difference $\Delta v=(v_t-v_b)/2$ will couple different superconductors at the interface. Its importance however, also controlled by $\Delta_d v/v_F$ and for values $\Delta_d v/v_F<0.5$ does not qualitatively change $\mu$-like shape of the differential conductance  (Fig.\ref{fig:exp:disorder}). Larger values of $\Delta_d v$ raise low-energy part of the differential conductance to positive energies. \fi

If we follow the standard practice \cite{blonder1982transition} of modelling
interface disorder by a delta-function potentials $v_l\delta(z)$, $\<v\>=(v_t+v_b)/2$ will play a role identical to that of  the disorder in a regular superconductor
by damping Andreev reflection at subgap energies and increasing electron-electron reflection in both symmetric and antisymetric channels. The dimensionless quantity controlling its importance is $\<v\>/v_F$, where $v_F$ is the Fermi velocity. Any difference $\Delta v=(v_t-v_b)/2$ between layers in the disorder parameter values
will couple the symmetric and antisymmetric quasiparticle channels.
The importance of these corrections
is however also controlled by $\Delta_d v/v_F$, and for values $\Delta_d v/v_F<0.5$  there is
no qualitative change in the $\mu$-like shape of the differential conductance  (Fig.\ref{fig:exp:disorder}). Sufficiently large values of $\Delta_d v/v_F$ raise the low-energy part of the differential conductance to positive values.

\section{Discussion}

In recent years experimenters have established moir\'e heterojunctions as attractive
platforms \cite{kennes2021moire} for new types of two-dimensional electron ground states.
In this Article, we have explored the phase diagram of bilayers with attractive effective
interactions ($\lambda_s$) within each layer and repulsive interactions ($\lambda_x$) between layers.
Two-dimensional bilayers with repulsive interlayer interactions can \cite{eisensteinQHE} have
exictonic insulator ground states that are counterflow superfluids and have spontaneous
interlayer phase coherence.
Our work is motivated by the discovery  of superconductivity, and hence attractive
effective interactions, in graphene bilayers \cite{cao2018unconventional,yankowitz2019tuning}
and trilayers \cite{park2021tunable}.
We therefore address the possibility of states
that have both superconductivity and counterflow superfluidity, and
study how the occurrence of such exotic states would be manifested by
generalized Andreev effects in separately contacted bilayers.

With this goal we have performed mean-field calculations for a model bilayer Hamiltonian with
two identical layers, and attractive intralayer and repulsive interlayer interactions.
We have neglected tunneling $\tau$ between the layers, assuming that the
bilayers are separated by dielectric layers.
Our calculations show that a mixed phase with both an interlayer exciton condensate and
superconductivity appears over a wide range of model parameters.
Interlayer coherence couples the two superconducting order parameters with
energy extrema occurring when the pair amplitudes are in phase (parallel) and out of phase (antiparallel).
We find that the energy of the antiparallel phase is higher than that of the pure
superconducting phase for any sufficiently smooth density of states.
In the parallel phase superconducting state,
interlayer coherence appears as a strong-coupling instability of states
with superconductivity in each layer, and
is established when $g_x \equiv \nu(\epsilon_F)\lambda_x>1$.
In twisted bilayer graphene we estimate that $\lambda_x\approx e^2/(4\pi k_F) \approx 100 \ meV nm^2$,
and that the density of states is around $\nu(\xi_s)\approx 10^{-2} \ nm^{-2} \ meV^{-1}$ \cite{xie2020weak}, implying that $g_x=1$ is within reach. Since the superconducting transition temperature
in tBG is $\sim 1 K$, the value we have chosen for $|g_s|\lessapprox 1$, is not unrealistic.
After interlayer coherence is established, further $\lambda_x$ increases
cause the superconducting order parameters to decay.
Note that self-consistent solution with both exciton order parameter
$x$ and direct intralayer pairing $\Delta_d$
nonzero also have non-zero inter-layer pairing $\Delta_i$.

As a convenience, we have limited our considerations to the case of two completely identical layers.
This condition is however not important for the stability of the phase. Consider for example the
case in which the Fermi energies in the two layers are slightly different: $\mu_{t,b}=\mu\pm\Delta\mu$.
This change induces a perturbation to the mean field Hamiltonian matrix of the
form $-\tau_z\sigma^z\Delta \mu$, where $\tau^z$ and $\sigma_z$ are the z-Pauli matrices
acting on layer and Nambu indices correspondingly.
After transformation to $+/-$ basis the perturbation acquire form $i\tau^y\sigma_z\mu$.
As a result, there will be no first-order contribution to the quasparticle energies.
At second order, only matrix elements between positive and negative energy states
will contribute. The leading order correction to the quasiparticle energies will be small,
$\sim (\Delta\mu^2/\sqrt{\Delta^2+x^2})$.  We conclude that small deviations from the perfect layer
symmetric case do not destroy the state, although they will change the mean-field-equation
solutions quantitatively.

\section{Acknowledgments}
We acknowledge support by  DOE grant DE-FG02-02ER45958. I.V.B. thanks Nemin Wei,  Ajesh Kumar and Rafi Bistritzer for discussions.

\section{Introduction}
Superconductivity is perhaps the most spectacular example of collective behavior of electrons.
The microscopic BCS \cite{bardeen1957theory,bardeen1957microscopic} theory of the phenomena explains superconductivity
in terms of bound (Cooper) electron pairs due to attractive interactions, mediated in the best understood case by
phonons. Condensation of Cooper pairs gives rise to a non-Fermi-liquid state with spontaneously broken $U(1)$ symmetry
characterized by a homogenous complex order parameter $\Delta$.  One experimental manifestation of superconductivity is
a non-Ohmic junction resistance with a normal metal.
At subgap values of the bias voltage on the interface, the only way for a single electron to penetrate inside a
superconductor is to form a Cooper pair. Therefore, by charge conservation, transmission through the interface should be accompanied by reflection of a hole, a process known as Andreev \cite{andreev1965thermal} reflection.  Studies of
Andreev reflection can be used to probe \cite{bradleySuperfluidAndreev, beenakkerColloquim} a state of interest and continue to gather both theoretical
\cite{de1995andreev,diez2012andreev,beenakker2006specular,nilsson2008splitting,stanescu2011majorana,kundu2013transport,mazin1999degreeofspinpol,mazin2001probing} and experimental \cite{sasaki2011topological,kastalsky1991observation,pothier1994flux,park2005andreev} attention.

In recent years experimentalists have uncovered examples of both electron-electron (Cooper) and
electron-hole (excitonic) pairing in atomically thin two-dimensional materials\cite{xi2016ising, li2018nontrivial,hsu2020inversion,eisensteinQHE, simonet2017anomalous,tutucEnhancement}.
In general Cooper pairing is driven by
attractive effective interactions, and electron-hole pair formation by repulsive interactions between electrons
in different bands.  Exciton condensates \cite{keldysh1968collective, hanamura1977condensation, littlewood2004models}
are coherent states of electrons and holes bound into pairs by the Coulomb interaction, and are
described by a mean-field theory that is identical
to BCS theory apart from a particle-hole transformation.
In bulk materials the concept of exciton condensation is partly ambiguous since the ordered states are difficult to distinguish from density-wave or nematic states \cite{fradkin2010nematic}.
Although exciton condensates were predicted theoretically more than 50 years ago,
their experimental identification in bulk materials has suffered as a consequence,
although progress has been achieved recently\cite{kogar2017signatures,mazuz2019dynamical,fogler2014high, high2012condensation}.
Experimental advances with two-dimensional materials have solved the ambiguity problem
by making it possible to prepare devices with strong interactions between subsystems
located in different layers that separately have nearly perfectly conserved particle number
Interlayer exciton condensation breaks the independent gauge invariance of the individual layers,
{\it i.e.} when the exciton order parameter $x$ is nonzero, there is spontaneous phase coherence
between layers and we no longer can perform gauge transformation independently.   There are still only a few
examples of convincing experimental evidence for exciton condensation in systems without
a strong external magnetic field, but this situation seems likely to change.

Recently, superconductivity has been observed in magic angle twisted bilayer
graphene (tBG) \cite{bistritzer2011moire, cao2018unconventional}, a two-dimensional system with
a moir\'e superlattice that yields extremely flat conduction and valence bands.
The strongest superconductivity appears when the
valence band is doped away from half filling.  Although several models \cite{wu2018theory, lian2019twisted, kozii2019nematic, chichinadze2020nematic, khalaf2020charged} for its superconductivity have been proposed, the mechanism remains unknown.
Twisted bilayers are the simplest examples of a rich variety of graphene multilayer moir\'e superlattice systems that
have received experimental attention recently \cite{cao2018correlated,lu2019superconductors,yankowitz2019tuning,park2021tunable}. The motivation for this paper is to
consider the possibility of realizing, in this flexible family of strongly-correlated electron systems, an exotic state in which interlayer exciton condensation and superconductivity
occur simultaneously.  Our target system is two twisted bilayers separated by
a hexagonal boron nitride tunnel barrier.  The superconductivity of the individual
twisted bilayers is already established.  Here we ask the following questions: i) can interlayer
exciton condensation \cite{rickhaus2020density} occur in principle between two-dimensional electron systems that are
superconducting ii) if so, how would such a state be most unambiguously detected?  So far
research on the coexistence of two phases has been scarce \cite{sager2020potential} partially because of the lack of a physical system with a potential to have both.
In section II  we describe the model we study and relationships between its order parameters.
In section III we describe the phases we have identified and phase diagrams as a function of the strength of the
interlayer repuslive interaction strength $g_x$ and the intralayer attractive interaction strength
$g_s$.  In the section IV we propose an Andreev drag measurement which can be used to
identify states with both types of order, taking advantage \cite{su2008make} of the possibility of
separate contacting to individual layers.
This method of probing bilayer exciton condensates has been
very successfully exploited in quantum Hall excitonic superfluids \cite{liu2017quantum,eisenstein2014exciton,yoon2010interlayer,li2017excitonic,li2019pairing}
Finally, in section V we discuss prospects for observing these
states and the importance of deviations from the highly symmetric point that we considered in our explicit calculations.
\iffalse
There are two main scenarios for spatially indirect exciton condensation.
In the first, a system is symmetric under the combination of layer interchange and
an electron-hole transformation. In this case, an electron-hole transformation in one
of the layers maps exciton condensates to superconductors that pair electrons in different layers and
occurs for arbitrarily weak interactions.
Consequently, charge transmission at low voltage occurs via a process
similar to Andreev reflection: to transfer an electron through an interface in one of the layers,
an electron should be reflected from the bottom layer.

In a scenario with two identical layers with the same Fermi level, exciton condensation is similar to itinerant electron
ferromagnetism with the spin degree-of-freedom replaced by a layer degree of freedom.  The phase will appear if the interaction $\lambda_x$ between the electrons and holes in two different layers is strong enough $g_x\equiv \lambda_x \nu(\epsilon_F)\geq 1$, where $\nu(\epsilon_F)$ is the density of states on the Fermi level.
 In this case, the eigenstates of the mean-field Hamiltonian have layer spinors that are coherent combinations of the
 individual layers, for example symmetric and antisymmetric.
 When viewed in this basis, transmission through the interface will be similar to free electron reflection from a
 potential.
  \fi

\section{Model}
We consider the simplest possible model that can have both interlayer exciton condensation and s-wave superconductivity.
The model assumes that  both the attractive intralayer and repulsive interlayer
interactions are momentum-independent, includes a layer degree of freedom ($l=t,b$),
and assumes either valley or spin singlet superconductivity
that is described by making a particle-hole transformation for one spin or valley in each layer,
The  Hamiltonian
\begin{multline}\label{model:full-hamiltonian}
	H
	=
	\sum_{pl\sigma}\xi_{pl}c^{\dagger}_{l\sigma}(p)c_{l\sigma}(p)
	\\+
	\frac{\lambda_s}{S}\sum_{pp'q}
	c^{\dagger}_{l\uparrow}(p+q)
	c^{\dagger}_{l\downarrow}(p'-q)
	c_{l\downarrow}(p')
	c_{l\uparrow}(p)
	\\+
	\frac{\lambda_x}{S}
	\sum_{pp'q}
	c^{\dagger}_{t\sigma}(p+q)
	c^{\dagger}_{b\sigma'}(p'-q)
	c_{b\sigma'}(p')
	c_{t\sigma}(p),
\end{multline}
where $c^\dagger_{l\sigma}$ $(c_{l\sigma})$ creates (annihilates) a
fermion in layer $l=t,b$ with pseudospin $\sigma=\uparrow,\downarrow$, and $S$ is the area of the sample.
In this work, we will focus on the case of layer-independent bands,
$\xi_{pl}\equiv \xi_p$, and reserve discussion of deviations from this symmetric limit to the end of the paper.
(As we explain in Section V coexistence does not occur for $\xi_{pt}\equiv -\xi_{pb}$.)
The pseudospin label refers to real spin in the case of spin-singlet superconductors and to
valley in the case of valley-singlet superconductors.  Below we refer to this degree of freedom as spin for simplicity.

Below we characterize the ordered states by their mean-field Hamiltonians,
and allow an exciton order parameter for each spin species,
\begin{equation}
	x_\sigma=\frac{\lambda_x}{S}\sum_p\<c_{t\sigma}(p)c_{b\sigma}^\dagger(p)\>,
\end{equation}
and superconducting order parameters both within and between layers,
\begin{equation}
	\Delta^{ll'}_{\sigma\sigma'}=-\frac{\lambda_{ll'}}{S}\sum_p\<c_{l\sigma}(p)c_{l\sigma'}(-p)\>,
\end{equation}
where $\lambda_{tt}=\lambda_{bb}=\lambda_s$ and $\lambda_{tb}=\lambda_{bt}=\lambda_x$.
The ordered state is therefore
characterized by six complex self-energies that vanish by symmetry if no symmetries are broken.

In a single-band superconductors, the pair potential can be chosen to be real because of global gauge invariance.
It follows that one real number, the gap, fully characterizes the superconducting phase.
In the present two-band system, we start with four complex pair potentials.
If we also allow excitonic particle-hole pairing for both spins, the ordered state is characterized by
12 real numbers.  Since only three phases can be chosen at will by exploiting the conservation of the number of particles of particles of
each spin in each layer, the energy can depend not only on absolute values of the order parameters,
but also on their phases, or, more specifically, on gauge invariant combinations of such phases.

 In the absence of interlayer superconducting coupling, there is only one gauge invariant combination of the phases: $\psi_t-\psi_b-\psi_{\uparrow}-\psi_{\downarrow}$, where $\psi_{\sigma=\uparrow,\downarrow}$ is the phase of the corresponding exciton pair potential and $\psi_{l=t,b}$ is the phase of the intralayer superconducting pairing potential. This combination changes sign  under  flip of layer. If we assume that the ground state preserves spin and layer invariance,
 we conclude that $\psi_t-\psi_b-\psi_{\uparrow}-\psi_{\downarrow}=\pi n$, where $n$ is an integer.
 In the gauge $\psi_b=0$, $\psi_{\uparrow}=0$, $\psi_{\downarrow}=0$, an even $n$ will correspond to a state with two supercondutving gaps having the same sign, while an odd $n$ will correspond to a state with superconducting gap that changes sign under the layer flip.

To extend this analysis to the case with interlayer pairing,
we identify four gauge invariant quantities that do not transform to themselves under layer or spin inversion:
$x_\uparrow \Delta^{bt}_{\uparrow\downarrow} {\Delta^{tt}_{\uparrow\downarrow}}^*$, $x_\downarrow \Delta^{bt}_{\downarrow\uparrow}{\Delta^{tt}_{\downarrow\uparrow}}^*$, $x^*_\uparrow \Delta^{tb}_{\uparrow\downarrow}{\Delta^{bb}_{\uparrow\downarrow}}^*$, and $x^*_\downarrow \Delta^{tb}_{\downarrow\uparrow}{\Delta^{bb}_{\downarrow\uparrow}}^*$
and require that their values do not change under these transformations. \footnote{One way to represent these quantities is to imagine space with 4 points corresponding to each kind of fermion. Then any gauge invariant quantity will correspond to a closed trajectory between these points.}
This condition restricts the absolute values of all self-energies
$|\Delta^t|=|\Delta^b|\equiv\Delta_d$, $|x_{\uparrow}|= |x_{\downarrow}|\equiv x$, $|\Delta^{tb}_{\uparrow\downarrow}|=|\Delta^{bt}_{\uparrow\downarrow}|\equiv \Delta_i$ as well as two gauge invariant combinations of phases
\begin{align}\label{phase-cond:1}
\psi_1-\psi_2=\psi_\downarrow-\psi_\uparrow-2\pi m,\\\label{phase-cond:2} \psi_t-\psi_b=\psi_\uparrow+\psi_\downarrow-2\pi k,
\end{align}
where $\psi_1$ is the phase of $\Delta^{bt}_{\uparrow\downarrow}$ and $\psi_2$ is the phase of  $\Delta^{tb}_{\uparrow\downarrow}$.
($d$, and $i$ in $\Delta_d$ and $\Delta_i$ are
intended to suggest {\it direct} and {\it indirect} pairing.) The system has one more independent gauge invariant combination of phases $\psi_t+\psi_b-\psi_1-\psi_2$. Unlike those we considered previously, this combination transforms to itself under
layer or spin flip, and therefore should be determined by energy minimization.

Using global gauge invariance related to conserved
particle numbers for each spin and valley $\psi_1$, $\psi_\uparrow$ and  $\psi_\downarrow$ can be chosen to be 0.  Then because of \eqref{phase-cond:1} $\psi_2=2\pi m$ and $\psi_-\equiv (\psi_t-\psi_b)/2=\pi k$. In what follows, we use this gauge.
This leaves us with
mean-field Hamiltonians that depend only on four parameters, three absolute values and a single phase,
whose values can be determined by solving self-consistent field equations.

\iffalse
leaves us with 6 real parameters, 3 of which can be gauged away because of number conservation. We conclude then it suffices to have 3 real parameters to describe the mean field state.

In single-band superconductors, the pairing order parameter can be chosen to be real because of global gauge invariance. Then there will be only one real number that characterizes the superconducting phase. In this simple case, one way to see this is to recall that the $U(2)$ Bogolyubov matrix has 4 independent real parameters. Among those, there is a single angle that characterizes $O(2)$ rotation as well as three complex phases. Because the Hamiltonian is gauge invariant both in new and the old basis we can always choose the gauge where all 3 phases are zero.

In the present system of two layers the situation is different and there are phases that cannot be gauged away.  In other words, the energy of the uniform ground state of the system will depend on combinations of order parameter phases.

The $U(4)$ Bogolyubov matrix has 16 independent real parameters. $6$ of them parametrize $O(4)$ rotation matrix. Among the remaining $10$ phases, $4+4$ can be made zeros and two remain. The energy can therefore depend on them. Now, if we assume that mirror (layer equivalence)  and spin-invariance is maintained, it follows that
$|\Delta^{t}_{\uparrow\downarrow}|=|\Delta^{b}_{\uparrow\downarrow}|$ and that $|x_{\uparrow}|=|x_{\downarrow}|$ as well as $|\Delta^{tb}_{\uparrow\downarrow}|=|\Delta^{bt}_{\uparrow\downarrow}|$. These 3 conditions will restrict reduce the number of independent rotations to $3$. We conclude then that 5 parameters (3 angles and 2 phases) fully characterize the mixed state.

The generation of phase masses is reminiscent to the case of the Josephson
junction, but in our case it has deeper consequences on the phase diagram itself.
Let's see in detail how global gauge transformations work in our system. In this work, we use a gauge where all but intralayer superconductivity order parameters are real and positive. In this gauge, $x$ gives mass to the difference of phases $2\psi_-=\psi_b-\psi_t$ of intralayer superconducting order parameters, and $\Delta_i$ gives mass to $2\psi_+=\psi_t+\psi_b$.
\fi

\section{Phase Diagram}
We saw previously that layer and spin symmetries restrict the value of $\psi_-$, but not
$\psi_+$. In this section, we allow both phases to take any value demonstrate explicitly that energies are minimized when spin/layer symmetries are respected: the four allowed
values of $\psi_-$ follow from the self-consistency equations.
The mean-field version of \eqref{model:full-hamiltonian} can be written as
\begin{equation}\label{phases:ham-matrix}
	H_{mf}=\sum_p \Psi^\dagger(p)
\mathcal{H}(p)
	\Psi(p),
\end{equation}
where $\Psi^\dagger(p)=\(c^\dagger_{t\uparrow},c_{t\downarrow},c^\dagger_{b\uparrow},c_{b\downarrow}\)$ is a vector in an extended Nambu space and the Hamiltonian matrix is
\begin{widetext}
\begin{equation}
  \mathcal{H}(p)=	\begin{pmatrix}
		\xi_p & \Delta_d e^{i(\psi_++\psi_-)} & x & \Delta_i \\
		\Delta_d e^{-i(\psi_++\psi_-)}&-\xi_p & \Delta_i & -x\\
		x & \Delta_i & \xi_p &\Delta_d e^{i(\psi_+-\psi_-)}\\
		\Delta_i & -x & \Delta_d e^{-i(\psi_+-\psi_-)} &-\xi_p
	\end{pmatrix}
\end{equation}.
\end{widetext}
The negative eigenvalues of the matrix have the form $\epsilon_{\pm}=-\(\eta^2\pm 2\alpha^2\)^{1/2}$, where $\eta^2=\xi^2+x^2+\Delta_i^2+\Delta^2_d$ and $\alpha^4=x^2 \Delta^2_d \sin^2(\psi_-)+2\Delta_d \Delta_i x\xi \cos(\psi_+)\cos(\psi_-)+\Delta_i^2\Delta_d^2\cos^2(\psi_+)+(x\xi)^2$.
Stable phases of the system are determined by minimization of the energy density.
At energy extrema,
\begin{equation} \label{phases:total-energy}
	E=\frac{2\Delta^2_d}{|\lambda_s|}-\frac{2\Delta_i^2}{\lambda_x}+\frac{2 x^2}{\lambda_x}+\sum_b\int d\xi \nu(\xi )\epsilon_{b}n_b(\xi),
\end{equation}
where the sum is over quasiparticle bands $b$, $\nu(\xi)$ is the quasiparticle density-of-states per spin and per layer,
and $n_b(\xi)$ is the Fermi occupation factor.
At zero temperature, only quasiparticle states with negative energies are filled.
\iffalse
 \begin{figure}[htb!]
	\includegraphics[width=\columnwidth]{Figures/energy_contour_vs_phi_fixed_orders_cust_cmap.png}
	\caption{Energy density \eqref{phases:total-energy} as a function of $\psi_-$ (x-axis) and $\psi_+$ (y-axis). Absolute values of all three order parameters fixed at some realistic values (on the main figure $\mu=0.1\Lambda$ $x=\mu$, $\Delta_d=0.01\Lambda$, $\Delta_i=\Delta_d/3$, and the density of states is the same as on Fig.\ref{fig:parallel-phase:underdoped} B), on the inset $\Delta_i=0$). The highest energy corresponds to red, the lowest energy corresponds to blue. We see that $\Delta_i$ lifts the energy degeneracy in $\psi_+$. Note that all minima are connected by $\psi_l\to \psi_l+2\pi n, \psi_l'\to\psi_l'+2\pi m$ for some integer $m$ and $n$ and therefore correspond to the same state. }
	\label{fig:parallel-phase:energy-vs-phi}
\end{figure}
\fi
Extrema of the energy are also extrema of Eq.~\eqref{phases:total-energy}, which we vary first
with respect to the phases $\psi_-$ and $\psi_+$.
We conclude that energy minima occur at extrema of $\alpha^2$, and that these
occur (independent of $\xi_p$)
when
\begin{equation}\label{phase-diagram:parallel-condition}
	\sin\(\psi_+\)=0 \text{ and } 	\sin\(\psi_-\)=0
	\end{equation}
	or
	\begin{equation}\label{phase-diagram:antiparallel-condition}
	\cos\(\psi_+\)=0 \text{ and } 	\cos\(\psi_-\)=0.
\end{equation}
The former condition is consistent with the conditions we derived for a mirror- and spin-symmetric state with $\Delta_i\neq0$.
We will call it a parallel phase.  The latter condition defines the antiparallel phase \cite{phase-condition-comment}.
We will see that the latter is consistent only with $\Delta_i=0$.
Note that each condition in fact corresponds to two different states. Namely, \eqref{phase-diagram:parallel-condition} can be satisfied with $(\psi_t,\psi_b)=(\pi,\pi)$, and $(0,0)$, that typically have different energies
when interlayer superconducting order is present.  Now we consider the antiparallel and parallel phases separately.

\subsection{Antiparallel phase}
We consider first extrema in which the phases $\psi_t$ and $\psi_b$ differ by $\pi$.
.
The main conclusion of this subsection is that for any sufficiently smooth density of states,
the energy of the antiparallel phase is always higher than that of a pure superconducting phase.
In other words this solution of the mean-field equations corresponds to a local energy maximum not
an energy minimum, and can therefore be discarded.

  \begin{figure*}[t]
	\includegraphics[width=18cm,height=5.5 cm]
	{Figures/order_parameters_lxstep0005_lsstep0005_random_den8nmin40.png}
	\caption{Order parameter dependence on dimensionless interaction constants $g_s=\nu \lambda_s$,
	and $g_x=\nu \lambda_x$ in the parallel phase,
	where $\nu$ is the parabolic band density-of-states.
	The exciton self-energy $x$
	Fermi energy is plotted in panel A.
	At small $|g_s|$ there is a phase transition at $g_x=1$, which is first order in our model, to a large
	$g_x$ state with $x \ne 0$.  The intralayer superconducting self-energy $\Delta_d$, plotted in
	panel B is always non zero.
	The interlayer superconducting self-energy $\Delta_i$ is plotted
	in panel C and acquires a finite value only when
	$x \ne 0$.  After the transition the purely intralayer superconducting phase
	to the mixed phase, $\Delta_d$ (B) tend to decrease. For large values of $|g_s|$ the excitonic transition is pushed to larger values of $g_x$. The inset in panel A shows the three self-energies as a function of $g_x$ at a fixed $g_s=-0.4$. These results were obtained for cutoff $\Lambda=10 \ meV$, $\epsilon_F=3 \ meV$, and a parabolic band model with a constant density of states. To treat both $f=0$ and $f=1$ minima on the equal footing, we allowed the interlayer superconducting gap $\Delta_i$ to take negative values. Self-energies are in units of meV. We minimized energy as a function of the Bogolyubov angles using the conjugate graduate method.}
	\label{fig:parallel-phase:underdoped}
\end{figure*}
In the antiparallel phase, variation of the energy density in Eq.~\ref{phases:total-energy} with respect to ${\Delta_i}$ combines with Eq.~\eqref{phase-diagram:antiparallel-condition} to yield the self-consistency equation
\begin{equation}
	\Delta_i
	=
	-\frac{\lambda_x \Delta_i}{4}\int\nu(\xi)\(\frac{1}{\epsilon_{+}}+\frac{1}{\epsilon_{-}}\).
\end{equation}
This equation clearly cannot be satisfied at any $\lambda_x>0$.  We conclude that
$\Delta_i=0$ in the antiparallel phase.
The quasiparticle band energies are therefore
\begin{equation}\label{antiparallel:band-energy}
	\pm \epsilon_{\pm}=\pm |\sqrt{\xi^2+\Delta_d^2}\pm x|\\
\end{equation}
The form of the band energies allows us to identify
the antiparallel phase as an exciton condensate (spontaneous interlayer phase coherent state)
formed on top of superconductors within each of the layers.
With ${\Delta_i}$ eliminated, the system of gap equations reduces to:
\begin{equation}\label{phase-diagram:antiparallel:x}
	x=
	\lambda_x
	\mathcal{N}
	(\sqrt{x^2-\Delta_d^2} \,),
\end{equation}
\begin{multline}\label{phase-diagram:antiparallel:delta}
	\Delta_d=-\frac{\lambda_s \Delta_d}{2}
	\int  \frac{\nu(\xi)}{\sqrt{\xi^2+\Delta_d^2}}+\\
	\frac{\lambda_s \Delta_d}{2}\int^{\sqrt{x^2-\Delta_d^2}}_{-\sqrt{x^2-\Delta_d^2}}\frac{\nu(\xi)}{\sqrt{\xi^2+\Delta_d^2}}.
\end{multline}
That is to say that the exchange self-energy is contributed by wavevectors that have different
occupation numbers for the two orthogonal quasiparticle layer spinors.
In Eq.~\eqref{phase-diagram:antiparallel:x}  $\mathcal{N}(\epsilon) \equiv 1/2\int_{-\epsilon}^{\epsilon} d\xi \ \nu(\xi)$.

\iffalse
A transition between two superconductors and the mixed phase becomes possible at the smallest value of the interlayer interaction strength $\lambda_x$ where the gap equations can be satisfied.  From equation \eqref{phase-diagram:antiparallel:x} it follows that at the transition point the value $\frac{x}{\mathcal{N}(\sqrt{x^2-\Delta_d^2})}$,  has the minimum possible value.  For any given $x$ and any continuous density of state it implies that $\Delta_d=0$, so that at a transition point the superconducting gap in each layer should close. It is clear then from \eqref{phase-diagram:antiparallel:delta} that at the same point $x$ should be nonzero to ensure the gap closing. We conclude then that the possible transition between two phases is of the first order.
\fi
A transition between pure superconductivity and an antiparallel
phase would occur at zero temperature if the latter phase were lower in energy.
We find that
\begin{multline}
\label{eq:energydiff}
	E_{m}-E_{sc}=
	\frac{2}{|\lambda_s|}\({\Delta_d}_m^2-{\Delta_d}_0^2\)
	\\-2\int_{-\mu}^\Lambda\nu(\xi)
	\(\sqrt{\xi^2+{\Delta_d}_m^2}-\sqrt{\xi^2+{\Delta_d}_0^2}\)
	\\+
	2\int^{\sqrt{x^2-{\Delta_d}^2_m}}_{-\sqrt{x^2-{\Delta_d}^2_m}}
	\nu(\xi)
	\sqrt{\xi^2+{\Delta_d}_m^2}
	-
	2x \mathcal{N}(\sqrt{x^2-{\Delta_d}^2_m})
	\leq
	0,
\end{multline}
 where we distinguish the superconducting gap ${\Delta_d}_m$ in the mixed phase from the superconducting gap ${\Delta_d}_0$ in the pure superconducting phase, and have used a self-consistency equation to reexpress the last term.
 Since the energy of the pure superconducting phase is minimized by ${\Delta_d}_0$,
 we can conclude that the sum of the first two lines in Eq.~\ref{eq:energydiff} is non-negative.
 Provided that the density of states does not significantly vary  on the scale of $\sqrt{x^2-{\Delta_d}^2_m}$,
 the third line can be approximated by
 \begin{multline}
 	2\nu_0{\Delta_d}_m^2\sinh^{-1}\(\frac{\sqrt{x^2-{\Delta_d}_m^2}}{{\Delta_d}_m}\),
 \end{multline}
where $\nu_0$ is the constant density-of-states, and the full energy difference is positive.
Hence, at least for a sufficiently smooth density of states,
the antiparallel phase will never be thermodynamically stable. We now analyze the parallel phase.

 \iffalse
 \begin{figure*}[t]
	\includegraphics[width=18cm]{Figures/order_parameters_lxstep01_lsstep002_mu5.png}
	\caption{Order parameter dependence on dimensionless interaction constants $g_s=\nu_{max}\lambda_s$, $g_x=\nu_{max}\lambda_x$ in the parallel phase for $\nu(\epsilon)=\nu_{max}\gamma^2/\(\gamma^2+(\xi-\xi_s)^2\)$with $\xi_s=0$ and $\gamma=1$.	Normalized exciton self-energy $x/\epsilon_F$ is plotted on A). For small absolute values of $g_s$, the exciton condensate appear at the point $g_x=1$ in accordance to the Stoner criterion. Interlayer superconductivity ${\Delta_i}$ is vanishingly since the Fermi energy is tuned to the peak in the density of states. In this calculation, $\Lambda=10 \ meV$, $\epsilon_F=5 \ meV$, $\gamma=1 \ meV$, $\alpha=1$, $\xi_s=0\ meV$. We performed a minimization of energy as a function of the Bogolyubov angles using conjugate graduate method.}
	\label{fig:parallel-phase:van-hove}
\end{figure*}
\fi
\subsection{Parallel phase}
The parallel phase is defined by the conditions $\sin(\psi_+)=0$, $\sin(\psi_-)=0$, so that $\psi_t$ and $\psi_b$ are equal to within a multiple of
$2\pi$.  This condition, however, leaves an ambiguity since the $\cos(\psi_+)\cos(\psi_-)$ factor
present in the dispersion relation can be either $+1$ or $-1$.
We distinguish these possibilities by introducing $f=0,1$ such that $\cos(\psi_+)\cos(\psi_-)=(-1)^f$.
The quasiparticle energies can be expressed as:
\begin{equation}\label{phase-diagram:parallel:dispersion}
	\epsilon_\pm=\sqrt{(\xi \pm x)^2+({\Delta_d} \pm (-1)^f {\Delta_i})^2}.
\end{equation}
We can interpret this phase as two superconductors: one with Cooper pairs formed out of symmetric combinations of layers, another -- from antisymmetric combinations.  Indeed, the mean-field Hamiltonian matrix \eqref{phases:ham-matrix} is block-diagonal in this basis. Note here that even though the transformation to the symmetric and antisymmetric states block-diagonalizes $\mathcal{H}(p)$, energy density $E$ \eqref{phases:total-energy} will still have terms that couple the two unless $|\lambda_s|=\lambda_x$.
Solutions with $f=0/1$ then correspond to phases with in which either symmetric or antisymmetric superconducting gaps are larger.
Variation of the full energy density \eqref{phases:total-energy} with respect to ${\Delta_i}$ yields
\begin{multline}\label{phase-diagram:parallel:energy-derivative:sigma}
	\frac{\delta E}{\delta {\Delta_i}}
	=
	-\frac{4{\Delta_i}}{\lambda_x}
	-
	{\Delta_i} \int
	\nu(\xi)
	\(
	\frac{1}{\epsilon_-}
	+
	\frac{1}{\epsilon_+}
	\)
	\\+
	(-1)^f
	\Delta_d \int
	\nu(\xi)
	\(
	\frac{1}{\epsilon_-}
	-
	\frac{1}{\epsilon_+}
	\)
\end{multline}
One conclusion we can immediately draw is that the ${\Delta_i}=0$ point is not an energy extremum
whenever both exciton condensates and superconductivity are present.
Note that the superconducting gap of the symmetric quasiparticles is $\Delta_d+ (-)^f{\Delta_i}$,
while for the antisymmetric quasiparticles is $\Delta_d- (-)^f {\Delta_i}$.
Because the exciton condensate breaks the symmetry between quasiparticle bands,
one should expect that $\Delta_d+{\Delta_i}\neq\Delta_d-{\Delta_i}$.
A complimentary explanation is that even though the interlayer interaction is repulsive,
the BCS instability present within each layer together with interlayer coherence
induces a response in the interlayer Cooper channel.

If we write the self-consistency equations in terms of
$\Delta_-\equiv\Delta_d-(-1)^f{\Delta_i}$, $\Delta_+\equiv \Delta_d+(-1)^f{\Delta_i}$ and $x$ we obtain,
\begin{multline}\label{phase-diagram:parallel:delta+}
	\frac{1}{\lambda_-^2-\lambda_+^2}(\lambda_+\Delta_+-\lambda_-\Delta_-)
	\\=
	\frac{\Delta_+}{2}
	\int \frac{d\xi \nu(\xi)}{\sqrt{(\xi+x)^2+\Delta_+^2}},
	\end{multline}
\begin{multline}
	\label{phase-diagram:parallel:delta-}
\frac{1}{\lambda_-^2-\lambda_+^2}(\lambda_+\Delta_--\lambda_-\Delta_+)
	\\=
	\frac{\Delta_-}{2}
	\int \frac{d\xi \nu(\xi)}{\sqrt{(\xi-x)^2+\Delta_-^2}},
\end{multline}
\begin{multline}
	\label{phase-diagram:parallel:x}
		\frac{4}{\lambda_+-\lambda_-}x=\int \frac{d\xi\nu(\xi)(x-\xi)}{\sqrt{(\xi-x)^2+\Delta_-^2}}
	\\+
	\int \frac{d\xi\nu(\xi)(x+\xi)}{\sqrt{(\xi+x)^2+\Delta_+^2}}.
\end{multline}
Here $\lambda_+\equiv(\lambda_s+\lambda_x)/2$ and $\lambda_-\equiv(\lambda_s-\lambda_x)/2$. Note that even though the mean-field
Hamiltonian matrix is block diagonal, the energy density has a term of the form $\Delta_+\Delta_-$ because $\lambda_s\neq\lambda_x$.
We solve this system numerically for the constant density-of-states case.
The mixed state solution summarized in Fig.\ref{fig:parallel-phase:underdoped} exists for $\nu(0)g_x\geq 1$.
  \begin{figure}[htb!]
	\includegraphics[width=\columnwidth]{Figures/atb_curve_combined_potentials.png}
	\caption{Probabilities for an electron coming from the top lead at energy $\omega$ to be reflected as a hole/electron in the same layer ($A_{tt}$/$B_{tt}$), and as a hole/electron in the opposite layer ($A_{tb}$/$B_{tb}$).
	The quantity $A_{tb}-B_{tb}$, proportional to the differential conductance in the bottom layer,
	changes sign as function of $\omega$. At low energies, the conductance is dominated by electron-electron reflection.
	At energies around superconducting gaps, the conductance is dominated by Andreev reflection processes.
	The inset shows a schematic experimental device with voltage bias $V$ across the top layer.
	When mesoscopic effects are neglected, the top layer conductance of the device is $\sigma_{t}(V/2)/2$,
	where $\sigma_{t}(V)$ is the differential conductance between the top left lead and the top layer,
	and the transconductance driven by excitonic order is $\sigma_{b}(V/2)/2$, where
	$\sigma_{b}(V)$ is the conductance between the top left lead and the bottom layer.}
	\label{fig:exp:setup}
\end{figure}
  \begin{figure}[htb!]
	\includegraphics[width=\columnwidth]{Figures/sigmatb_dis_curve_plotDp0.1Dm0.15x0.8V4V.png}
	\caption{Influence of interface disorder on the difference between hole (Andreev process) and
	electron interlayer reflection coefficients $A_{tb}-B_{tb}$.
	Disorder is modeled by a delta-function potential at the interface $v_l\delta(z)$.
	On each subplot, the ratio of $\<v\>=(v_t+v_b)/2$ to the Fermi momentum is kept fixed. $\Delta v/p_F=(v_t-v_b)/(2p_F)$ is swept from -1.0 (red) to +0.8 (black) with step 0.1.}
	\label{fig:exp:disorder}
\end{figure}
In the numerical calculation, to treat both $f=0$ and $f=1$ minima on the equal footing, we allowed the interlayer superconducting gap $\Delta_i$ to take negative values. Then negative value of $\Delta_i$ corresponds to the $f=1$ minimum while positive values will correspond to the $f=0$ minimum. Our calculation shows that $f=1$ minimum tend to be more energetically favorable for larger values of intralayer attraction, which has a simple physical meaning: for $|\lambda_s|\approx \lambda_x$ the symmetric and the antisymmetric superconductors are effectively decoupled and thus larger gap correspond to a larger effective Fermi energy.
Along the first order phase boundary in $g_s$/$g_x$ space separating mixed and pure superconducting phases,
the two states have identical energies.  It follows that the differential energy changes
$dE_{sc}$ in the superconducting phase and $dE_{mix}$ in the mixed phase should be the same.
It follows that the slope of the phase boundary lines
\begin{equation}\label{phase-diagram:parallel:coexistence-curve}
	\frac{dg_x}{dg_s}
	=
	\frac{\frac{\d E_{sc}}{\d g_s}-\frac{\d E_{mix}}{\d g_s}}{\frac{\d E_{mix}}{\d g_x}}.
\end{equation}
Partial derivatives in the formula will be proportional to order parameters in each phase: $\d E_{sc,mix}/\d g_s=-2\Delta_{sc,mix}^2/g_s^2$, $\d E_{sc}/\d g_x=2(x^2-{\Delta_i}^2)/g_x^2\approx 2\mu^2/g_x^2$. Formation of the mixed state accompanied by the decrease in the superconducting gap at least for large enough $|g_s|$. Therefore, $\frac{dg_x}{dg_s}<0$ --  for larger attractive interaction $g_s$ transition will happen at a larger critical $g_x$ in agreement with Fig. \ref{fig:parallel-phase:underdoped}. \iffalse Assuming that $\Delta_{sc}^2-\Delta_{mix}^2=g_s^2 c$ with $c$ being weakly dependent on $g_s$ and $g_x$ at the interface, by integrating \label{phase-diagram:parallel:coexistence-curve} from point $(g_x,g_s)=(1,0)$ we get $g_x=(1+cg_s/\mu^2)^{-1}$ along the interface.\fi

\section{Non-local Andreev Reflection}

Having established that an exotic mixed state with both excitonic and Cooper pair order can exist,
we now discuss Andreev reflection experiments in separately contacted bilayers that
can be used to detect their presence. Consider metallic leads separately connected to the top and the bottom
layers of our system. An electron incident on the interface between the top lead and the top layer
has equal weight in symmetric and antisymmetric quasiparticle
channels $\ket{t}=(\ket{+}+\ket{-})2^{-1/2}$. Since the parallel state does not couple $\ket{+}$ with $\ket{-}$, we can
consider their transmission contributions separately, at least if disorder is neglected.
Incoming quasiparticles are reflected with probabilities $A_{tt}$ and $A_{tb}$ as holes,
and with probabilities $B_{tt}$ and $B_{tb}$ as electrons, in the top and bottom layers respectively.
The calculation of these reflection probabilities is discussed below.

If we apply $V$ voltage across the top part of the interface, the
top and bottom layer differential conductances will be connected to reflection coefficients in a normal fashion.
The total currents in the top and bottom layers are
\begin{equation}
I_{t/b} = e \nu_0 \int_{0}^{eV} d \omega \; \sigma_{t/b}(V)
\end{equation}
where
\begin{eqnarray}
 \sigma_t(V)&\equiv& dj_t/dV=e^2 v_f\nu_0 (1-B_{tt}(V)+A_{tt}(V)), \nonumber \\
 \sigma_b(V)&\equiv& dj_b/dV=e^2 v_f\nu_0 (A_{tb}(V)-B_{bt}(V)),
 \end{eqnarray}
$e$ is an electron charge, and $v_f\nu_0$ is the product of the Fermi velocity and the
density of states.  Note that the combination
 $A_{tb}-B_{tb}$ of cross-layer reflection probabilities can be measured experimentally.
 We derive expressions for the scatterring amplitudes in the Appendix using the Bogoliubov-de-Gennes equations. The amplitude to reflect an electron as a hole in either symmetric or antisymmetric states, at least in the limit $x\ll \mu$, around 1 at $\omega=0$ similarly to a single N-S junction. However, the amplitude to reflect an incoming electron as an electron is no longer negligible even in the $\Delta_d\ll \mu$ case. The electron reflection amplitude at $\omega=0$ is given by $b_\pm\approx\mp x/(2\mu\pm x)$.
 Because $b_+$ and $b_-$ have opposite signs, the probability to reflect an electron to the bottom layer $B_{tb}=|b_+-b_-|^2/4$ is substantilly higher than the probability of the corresponding Andreev process, which is almost zero.
For the the same reason electron reflection in the top layer is unimportant for small $\omega$ and instead the Andreev process dominates (Fig.\ref{fig:exp:setup}).
When energy is equal to the $+$ or $-$ gap and the gap $\Delta_\pm/\mu\ll 1$,
the corresponding momenta in the leads and device match, giving rise to peaks in Andreev reflection and minima in electron reflection.
As a result the Andreev process is dominant over electron reflection in the bottom layer.
Altogether, we expect $A_{tb}-B_{tb}$ to have a $\mu$-like shape, with negative values
around $\omega=0$ and positive peaks at the superconducting gaps (Fig.\ref{fig:exp:setup}).

\iffalse $A_{tb}-B_{tb}$ is a quantity one can measure experimentally.
 If we apply $V$ voltage across the top part of the interface, top and bottom differential conductance will be connected to reflection coefficients in a normal fashion: $\sigma_t(V)\equiv dj_t/dV=e^2v_f\nu_0 (1-B_{tt}(V)+A_{tt}(V))$, $\sigma_b(V)\equiv dj_b/dV=e^2v_f \nu_0 (A_{tb}(V)-B_{bt}(V))$, where $e$ is an electron charge, and $v_f\nu_0$ is a characteristic of a lead. As a result, we expect the differential conductance in the bottom layer to change sign several times as function of the applied voltage, with robust positive peaks at values corresponding to the superconducting gaps.
If we model disorder by a delta-function $v_l\delta(z)$ at the interface in each layer, $\<v\>=(v_t+v_b)/2$ will play a role identical to the disorder in a regular superconductor: by damping Andreev reflection at subgap energies and increasing electron-electron reflection in symmetric and antissymetric superconductors separately. A dimensionless quantity controlling its importance is $\<v\>/v_F$, where $v_F$ is the Fermi velocity in the lead. A difference $\Delta v=(v_t-v_b)/2$ will couple different superconductors at the interface. Its importance however, also controlled by $\Delta_d v/v_F$ and for values $\Delta_d v/v_F<0.5$ does not qualitatively change $\mu$-like shape of the differential conductance  (Fig.\ref{fig:exp:disorder}). Larger values of $\Delta_d v$ raise low-energy part of the differential conductance to positive energies. \fi

If we follow the standard practice \cite{blonder1982transition} of modelling
interface disorder by a delta-function potentials $v_l\delta(z)$, $\<v\>=(v_t+v_b)/2$ will play a role identical to that of  the disorder in a regular superconductor
by damping Andreev reflection at subgap energies and increasing electron-electron reflection in both symmetric and antisymetric channels. The dimensionless quantity controlling its importance is $\<v\>/v_F$, where $v_F$ is the Fermi velocity. Any difference $\Delta v=(v_t-v_b)/2$ between layers in the disorder parameter values
will couple the symmetric and antisymmetric quasiparticle channels.
The importance of these corrections
is however also controlled by $\Delta_d v/v_F$, and for values $\Delta_d v/v_F<0.5$  there is
no qualitative change in the $\mu$-like shape of the differential conductance  (Fig.\ref{fig:exp:disorder}). Sufficiently large values of $\Delta_d v/v_F$ raise the low-energy part of the differential conductance to positive values.

\section{Discussion}

In recent years experimenters have established moir\'e heterojunctions as attractive
platforms \cite{kennes2021moire} for new types of two-dimensional electron ground states.
In this Article, we have explored the phase diagram of bilayers with attractive effective
interactions ($\lambda_s$) within each layer and repulsive interactions ($\lambda_x$) between layers.
Two-dimensional bilayers with repulsive interlayer interactions can \cite{eisensteinQHE} have
exictonic insulator ground states that are counterflow superfluids and have spontaneous
interlayer phase coherence.
Our work is motivated by the discovery  of superconductivity, and hence attractive
effective interactions, in graphene bilayers \cite{cao2018unconventional,yankowitz2019tuning}
and trilayers \cite{park2021tunable}.
We therefore address the possibility of states
that have both superconductivity and counterflow superfluidity, and
study how the occurrence of such exotic states would be manifested by
generalized Andreev effects in separately contacted bilayers.

With this goal we have performed mean-field calculations for a model bilayer Hamiltonian with
two identical layers, and attractive intralayer and repulsive interlayer interactions.
We have neglected tunneling $\tau$ between the layers, assuming that the
bilayers are separated by dielectric layers.
Our calculations show that a mixed phase with both an interlayer exciton condensate and
superconductivity appears over a wide range of model parameters.
Interlayer coherence couples the two superconducting order parameters with
energy extrema occurring when the pair amplitudes are in phase (parallel) and out of phase (antiparallel).
We find that the energy of the antiparallel phase is higher than that of the pure
superconducting phase for any sufficiently smooth density of states.
In the parallel phase superconducting state,
interlayer coherence appears as a strong-coupling instability of states
with superconductivity in each layer, and
is established when $g_x \equiv \nu(\epsilon_F)\lambda_x>1$.
In twisted bilayer graphene we estimate that $\lambda_x\approx e^2/(4\pi k_F) \approx 100 \ meV nm^2$,
and that the density of states is around $\nu(\xi_s)\approx 10^{-2} \ nm^{-2} \ meV^{-1}$ \cite{xie2020weak}, implying that $g_x=1$ is within reach. Since the superconducting transition temperature
in tBG is $\sim 1 K$, the value we have chosen for $|g_s|\lessapprox 1$, is not unrealistic.
After interlayer coherence is established, further $\lambda_x$ increases
cause the superconducting order parameters to decay.
Note that self-consistent solution with both exciton order parameter
$x$ and direct intralayer pairing $\Delta_d$
nonzero also have non-zero inter-layer pairing $\Delta_i$.

As a convenience, we have limited our considerations to the case of two completely identical layers.
This condition is however not important for the stability of the phase. Consider for example the
case in which the Fermi energies in the two layers are slightly different: $\mu_{t,b}=\mu\pm\Delta\mu$.
This change induces a perturbation to the mean field Hamiltonian matrix of the
form $-\tau_z\sigma^z\Delta \mu$, where $\tau^z$ and $\sigma_z$ are the z-Pauli matrices
acting on layer and Nambu indices correspondingly.
After transformation to $+/-$ basis the perturbation acquire form $i\tau^y\sigma_z\mu$.
As a result, there will be no first-order contribution to the quasparticle energies.
At second order, only matrix elements between positive and negative energy states
will contribute. The leading order correction to the quasiparticle energies will be small,
$\sim (\Delta\mu^2/\sqrt{\Delta^2+x^2})$.  We conclude that small deviations from the perfect layer
symmetric case do not destroy the state, although they will change the mean-field-equation
solutions quantitatively.

\section{Acknowledgments}
We acknowledge support by  DOE grant DE-FG02-02ER45958. I.V.B. thanks Nemin Wei,  Ajesh Kumar and Rafi Bistritzer for discussions.

redpajama_set_name: RedPajamaArXiv