# Pairing susceptibility in the multi-layer Hubbard model

Rayan Farid<EMAIL_ADDRESS>Department of Physics and Physical Oceanography,
Memorial University of Newfoundland and Labrador, St. John’s, Newfoundland &
Labrador, Canada A1B 3X7  J. P. F. LeBlanc<EMAIL_ADDRESS>Department of
Physics and Physical Oceanography, Memorial University of Newfoundland and
Labrador, St. John’s, Newfoundland & Labrador, Canada A1B 3X7

###### Abstract

We present a systematic study of the interaction, doping, and layer dependence
of the $d_{x^{2}-y^{2}}$-wave pairing susceptibility of the Hubbard model for
a stacked 2D square lattice. We perform a multi-index perturbative expansion
up to fourth-order to obtain coefficients in powers of the Hubbard $U$, the
inter-layer $V$, and the pair-hopping $J$ interactions. We evaluate the vertex
diagrams that contribute to the pairing susceptibility for $\ell=2,3,4$
layered models in the _U-V-J_ interaction space. This provides unprecedented
access to the pairing amplitudes, allowing us to identify the processes that
enhance or reduce pairing. We distinguish pairing within the diagonal channel,
$P^{\parallel}_{d}$, and off-diagonal channel, $P^{\perp}_{d}$, and find that,
in the absence of $J$, the qualitative behavior of the layered system is
equivalent to the single-layer model. In the presence of $J$, we show that
pairing is enhanced sublinearly with increasing $\ell$ and is primarily
mediated by the $P^{\perp}_{d}$ component and find which coefficients and
diagram sets are responsible. Finally, we construct a generalized
$\ell$-dependent equation for the d-wave pair susceptibility to speculate
pairing beyond $\ell=4$.

##  I Introduction

Understanding the mechanisms that give rise to superconductivity for models or
materials has been a lofty goal within condensed matter theory and the focus
of numerous studies. In the case of strongly correlated models, this pursuit
has driven advancement in numerical methods,[1, 2, 3, 4, 5, 6, 7] with the
Hubbard model playing a special role. This is due to its relevant early
application to understanding cuprates[8, 9] (LSCO) and the similarities
between the Hubbard phase diagram and the cuprate phase diagram[10, 11] both
for finite systems[12, 13, 14] and in the ground state,[15, 16]. Despite
successes, there remains a large gap between the physics of the 2D Hubbard
model and the complexity of layered cuprate materials. Further, the local
assumption of the Hubbard interaction seems to be insufficient to describe
superconducting ground states, and it has been shown that non-local attraction
gives rise to pairing and delicate interplay between charge-density wave and
superconducting ground states, stripes, and other exotic $T=0$ phases.[17, 18,
19] At finite temperatures, the importance of non-local interactions is also
seen, and together, these ideas motivate us to expand our considerations
beyond local interactions.[20, 21, 22, 23] For Hubbard-like models, adding to
the interaction space increases the computational complexity exponentially,
akin to a multi-band problem. Interestingly, there is no fundamental
difference between a general multi-band problem and one arising by coupling
two or more 2D planes.

In this work, we study $d$-wave pairing processes on infinite 2D planes
coupled via kinetic hopping, onsite and interplane repulsion, as well as both
local and non-local pair hopping processes by computing the uniform correlated
pairing susceptibility $P_{d}(\mathbf{q}=0)$ that is nothing more than the
vertex contributions of the zero-frequency particle-particle susceptibility
projected into an appropriate symmetry channel. We do this via a perturbative
approach using the state-of-the-art symbolic tool Algorithmic Matsubara
Integration (AMI).[24, 25] Initially, we produce the symbolic representation
of all the particle-particle topologies.[26] Then, we iteratively assign the
band index to each diagram symbolically, generating a complete set of
particle-particle diagram sets and the Matsubara sums are then resolved by
repeated application of the residue theorem. The spatial dependency is
stochastically integrated out, ensuring that our results are in the
thermodynamic limit and do not suffer from finite-size effects.

We categorize the computed diagram sets based on the type and number of
interaction lines present, resulting coefficients for powers of $U,V,$ and
$J$. These coefficients provide full access to the interaction space in the
weak coupling regime without additional computational expense, an advantage
unique to our symbolic approach. Leveraging this, we extensively probe the
interaction space and perform a comparative analysis of the dominant
coefficients to identify diagrammatic processes that enhance or suppress
pairing. We find that the presence of a local interaction has no qualitative
impact on coupled planes compared to a single 2D plane. With the inclusion of
local and non-local pair hopping processes, layer-dependent features begin to
appear, where the $P_{d}$ changes sublinearly as a function of layer number.
In an instructive manner, we extensively explore this up to a four-layered
system as a function of fixed chemical potential $\mu$ and identify the
diagram contributions to each pairing channel. In the supplemental materials
we examine further the role of inter-plane tunneling strength and
unconstrained $\mu$ phase space in pairing as well as study $q=(\pi,\pi)$ spin
susceptibility using the same approach.

##  II Models and Methods

###  II.1 Non-interacting Hamiltonian

We briefly summarize the tight-binding models for $\ell=2,3,$ and $4$. This is
necessary due to the large parameter space involved, which we discuss in
detail below.

####  II.1.1 Bilayer Hubbard Model

We consider a Bilayer ($\ell=2$) Hubbard Model with intra-planer nearest, $t$,
and next-nearest, $t^{\prime}$, single particle hopping as well as inter-
planer tunneling between the two planes $t_{z}$ [27]. Setting the lattice
constant $a=1$, the non-interacting dispersion of the intra-layer component is
given by

$\epsilon^{\parallel}(\textbf{k})=-2t[\cos(k_{x})+\cos(k_{y})]-4t^{\prime}[\cos(k_{x})\cos(k_{y})].$
(1)

To account for the anisotropic nature of cuprate splitting, we incorporate the
momentum-dependent inter-layer hybridization term $t_{\perp}(k)$ alongside an
isotropic component $t_{bs}$ given by

$t_{z}(k)=t_{bs}+t_{\perp}[\cos(k_{x})-\cos(k_{y})]^{2}.$ (2)

The $t_{z}(k)$ term is valid for any general $(\ell)$ layered Cuprates and is
derived from tight-binding parameters fitted to LDA calculations and confirmed
by experimental observation [28, 29, 30, 31]. Note that the c-axis lattice
constant is set to $c=1$ for simplification. This formulation leads to the
construction of the complete non-interacting Hamiltonian $H_{0}$:

$\displaystyle H_{0}^{\ell=2}$
$\displaystyle=\sum_{k\sigma}\vec{c}^{\dagger}_{k\sigma}\hat{\xi}_{k}\vec{c}_{k\sigma}$
$\displaystyle H_{0}^{\ell=2}$
$\displaystyle=\sum_{k}\begin{bmatrix}c^{1\dagger}_{k\sigma}&c^{2\dagger}_{k\sigma}\end{bmatrix}\begin{bmatrix}\epsilon_{\parallel}(k)&t_{z}(k)\\\
t_{z}(k)&\epsilon_{\parallel}(k)\end{bmatrix}\begin{bmatrix}c^{1}_{k\sigma}\\\
c^{2}_{k\sigma}\end{bmatrix},$ (3)

where $c^{1\dagger}_{k\sigma}(c^{1}_{k\sigma})$ and
$c^{2\dagger}_{k\sigma}(c^{2}_{k\sigma})$ represent the creation
(annihilation) operators of spin $\sigma$ electron with momenta $k$ in plane
index $i$ indicated by the superscript $1,2$. Diagonalizing $\hat{\xi}_{k}$
yields:

$H_{0}^{\ell=2}=\sum_{k\sigma}\begin{bmatrix}a^{1\dagger}_{k\sigma}&a^{2\dagger}_{k\sigma}\end{bmatrix}\begin{bmatrix}\tilde{\epsilon}^{1}(k)&0\\\
0&\tilde{\epsilon}^{2}(k)\end{bmatrix}\begin{bmatrix}a^{1}_{k\sigma}\\\
a^{2}_{k\sigma}\end{bmatrix}$

where $a^{1\dagger}_{k\sigma}(a^{1}_{k\sigma})$ and
$a^{2\dagger}_{k\sigma}(a^{2}_{k\sigma})$ represents the
creation(annihilation) operator of spin $\sigma$ the bonding ($k_{z}=0$) and
anti-bonding ($k_{z}=\pi$) momenta plane respectively. The $a^{1,2}_{k\sigma}$
eigenvector can be expressed as linear combination of uncorrelated basis
$a^{1,2}_{k\sigma}=\frac{1}{\sqrt{2}}(c^{1}_{k\sigma}\pm c^{2}_{k\sigma})$.
The eigenvalue of $a^{1,2}_{k\sigma}$ represents the full dispersion with
inter-planer hybridization given by
$\tilde{\epsilon}^{1,2}(\textbf{k})=\epsilon^{\parallel}(k)\pm
t_{z}(k)-\mu_{1,2}$ where $\mu_{1,2}$ is the chemical potential that controls
the filling in the respective bands. This reduces the bilayer Hamiltonian into
an effective two-orbital problem where the $k_{z}$ index can treated as a
separate band index. An isotropic bilayer band splitting would be obtained for
$t_{\perp}=0$ and $t_{bs}\neq 0$ case. Conversely, $t_{\perp}\neq 0$ and
$t_{bs}=0$ would result in highly anisotropic splitting with a maximum gap of
$4t_{\perp}$ at the anti-nodal point $k=[\pi,0]$ and dispersion-less in the
c-axis along diagonal lines $|k_{x}|=|k_{y}|$.

####  II.1.2 Trilayer and Quadlayer Hubbard model

Considering inter-planar single particle hopping between only the adjacent
planes and utilizing the same intra-layer dispersion defined in Eq. (1), the
trilayer $(\ell=3)$ Hamiltonian can be defined as

$\displaystyle H^{\ell=3}_{0}$
$\displaystyle=\sum_{k\sigma}\begin{bmatrix}c^{1\dagger}_{k\sigma}&c^{2\dagger}_{k\sigma}&c^{3\dagger}_{k\sigma}\end{bmatrix}\begin{bmatrix}\epsilon^{\parallel}(k)&t_{z}(k)&0\\\
t_{z}(k)&\epsilon^{\parallel}(k)&t_{z}(k)\\\
0&t_{z}(k)&\epsilon^{\parallel}(k)\\\
\end{bmatrix}\begin{bmatrix}c^{1}_{k\sigma}\\\ c^{2}_{k\sigma}\\\
c^{3}_{k\sigma}\end{bmatrix}$ (4)

Similar to the Bilayer case, the Hamiltonian above can be diagonalized on the
basis of the eigenvectors. The
$a^{1,3}_{k\sigma}=c^{1}_{k\sigma}\pm\sqrt{2}c^{2}_{k\sigma}+c^{3}_{k\sigma}$
eigenvectors yield a bonding $(k_{z}=0)$ and an anti-bonding $(k_{z}=\pi)$
band corresponding to the two outer plane in reciprocal space, given by
$\tilde{\epsilon}^{1,3}_{k\sigma}=\epsilon_{\parallel}(k)\pm\sqrt{2}t_{z}(k)-\mu_{1,3}$.
Meanwhile, the eigenvector
$a^{2}_{k\sigma}=(c^{1}_{k\sigma}-c^{3}_{k\sigma})/\sqrt{2}$ results in an
inner planar dispersion $(k_{z}=\frac{\pi}{2})$ given by
$\tilde{\epsilon}^{2}_{k\sigma}=\epsilon{{}^{\parallel}}(k)-\mu_{2}$ that is
independent of $t_{z}(k)$, thereby identical to the single-layer model. This
parameterization allows us to treat $k_{z}=[0,\frac{\pi}{2},\pi]$ momenta
labels as band indices where
$\tilde{\epsilon}^{1}_{k\sigma}\leq\tilde{\epsilon}^{2}_{k\sigma}\leq\tilde{\epsilon}^{3}_{k\sigma}$,
effectively recasting the trilayer model to a three-orbital problem. With the
inclusion of an additional layer, the quad-layer $(\ell=4)$ Hamiltonian is
constructed as

$H^{\ell=4}_{0}=\sum_{k\sigma}\begin{bmatrix}c^{1\dagger}_{k\sigma}&c^{2\dagger}_{k\sigma}&c^{3\dagger}_{k\sigma}&c^{4\dagger}_{k\sigma}\end{bmatrix}\\\
\begin{bmatrix}\epsilon^{\parallel}(k)&t_{z}(k)&0&0\\\
t_{z}(k)&\epsilon^{\parallel}(k)&t_{z}(k)&0\\\
0&t_{z}(k)&\epsilon^{\parallel}(k)&t_{z}(k)\\\
0&0&t_{z}(k)&\epsilon^{\parallel}(k)\\\
\end{bmatrix}\begin{bmatrix}c^{1}_{k\sigma}\\\ c^{2}_{k\sigma}\\\
c^{3}_{k\sigma}\\\ c^{4}_{k\sigma}\end{bmatrix}.$ (5)

The $\ell=4$ Hamiltonian is readily diagonalized to obtain four distinct
dispersions given by
$\tilde{\epsilon}^{1,2}_{k\sigma}=\epsilon^{\parallel}(k)-t_{z}(k)(\sqrt{5}\pm
1)/2-\mu_{1,2}$ and
$\tilde{\epsilon}^{3,4}_{k\sigma}=\epsilon^{\parallel}(k)+t_{z}(k)(\sqrt{5}\pm
1)/2-\mu_{3,4}$ parameterized by band indices
$k_{z}=[0,\frac{\pi}{3},\frac{2\pi}{3},\pi]$ such that
$\tilde{\epsilon}^{1}_{k\sigma}\leq\tilde{\epsilon}^{2}_{k\sigma}\leq\tilde{\epsilon}^{3}_{k\sigma}\leq\tilde{\epsilon}^{4}_{k\sigma}$.

####  II.1.3 Tight binding fitting parameters

The challenge in conducting a comparative study among $\ell=1,2,3,4$ models
lies in the vast parameter space one needs to probe. This is exacerbated by
the fact that LDA fits for tight-binding parameters in Cuprates vary depending
on which family of homologous series it belongs to. This necessitates adopting
the same tight-binding parameters across the $\ell$ layered models to
facilitate a direct and clear comparison. First, energies are normalized in
the unit of $t$ by setting $t=1$. Following existing literature on theoretical
studies of bilayer and trilayer systems, we fix $t^{\prime}/t=-0.3$ and
anisotropic inter-layer tunneling strengths within the range $0.05\leq
t_{\perp}/t\leq 0.150$ [32, 33, 34, 35, 36, 37]. Here, we set the isotropic
tunneling term $t_{bs}/t=0$ and primarily study the $t_{\perp}/t=0.125$
anisotropic case for all $\ell$-layered systems. The results for non-zero
isotropic tunneling $t_{bs}$ and the dependence of $t_{\perp}$ on pairing are
provided in the supplemental materials. We have fixed $\mu$ in each band to
reduce the parameter space further. The result for unconstrained $\mu$ in each
band is studied only for bilayer and trilayer cases in the supplemental
materials.

###  II.2 Interacting Hamiltonian

The interacting component of the Hamiltonian includes both intra-layer and
inter-layer components, represented as follows:

$H_{U}=\frac{U}{2}\sum_{i,\ell,\sigma\neq\sigma^{\prime}}n_{i\ell\sigma}n_{i\ell\sigma^{\prime}}+\frac{V}{2}\sum_{i,\ell\neq\ell^{\prime},\sigma}n_{i\ell\sigma}n_{i\ell^{\prime}\sigma^{\prime}}\\\
+J\sum_{i,\ell\neq\ell^{\prime},\sigma\neq\sigma^{\prime}}c^{\dagger}_{i\ell\sigma}c^{\dagger}_{i\ell\sigma^{\prime}}c_{i\ell^{\prime}\sigma^{\prime}}c_{i\ell^{\prime}\sigma}\\\
+J^{\prime}\sum_{i,\ell\neq\ell^{\prime},\sigma\neq\sigma^{\prime}}c^{\dagger}_{i\ell\sigma}c^{\dagger}_{i+\delta,\ell\sigma^{\prime}}c_{i+\delta,\ell^{\prime}\sigma^{\prime}}c_{i\ell^{\prime}\sigma}.$
(6)

Here, $n_{i\ell\sigma}$ denotes the number operator for electron with $\sigma$
spin in the $i$th site and $\ell$th layer and
$\delta\in\\{\hat{-x},\hat{-y},\hat{x},\hat{y}\\}$ refers to the nearest
neighbor site. $U$ represents the onsite local interaction between electrons
with opposite spins, and $V$ denotes the inter-layer interaction between two
electrons in adjacent layers, separated by a distance of the c-axis lattice
constant. We also incorporate inter-layer pair hopping of two configurations:
local _onsite_ pair hopping $J$ and a nonlocal _offsite_ pair hopping
$J^{\prime}$. In the offsite configuration, a pair of electrons with opposite
spins sits on the nearest neighbor and scatters to the next layer without any
spin flip. Both $J$ and $J^{\prime}$ could be considered as microscopic
mechanisms behind Josephson coupling, but rather than being a second-order
kinetic process, they are interpreted here as Coulombic processes [38, 39,
32]. This formulation allows for the scattering of a Cooper pair between the
next adjacent and next-next adjacent momenta layers via two or more inter-
layer pair hopping interactions in the $\ell=3$ and $\ell=4$ Hubbard model.
The long-range Coulombic process $J^{\prime}$, when Fourier transformed into
momentum space has a form factor
$J^{\prime}(q)=2(\cos(q_{x})+\cos(q_{y}))$.[32] We incorporate the effect of
$J^{\prime}$ inside every $J$ interaction via parameterization of $J$ as
$J(q)\rightarrow J(1+2J^{\prime}/J[\cos(q_{x})+\cos(q_{y})])$ which requires
fixing the ratio $J^{\prime}/J$ to a specific value.

###  II.3 Pair Correlation in Multi-layered Hubbard model

To investigate the pairing potential in the normal state of the layered
Hubbard model, we calculate the uniform ($\mathcal{\mathbf{Q}}=0$) and static
($\mathbf{\Omega}=0$) pairing susceptibility from linear response theory [40,
41]

$\scalebox{1.1}{{\hbox{\chi}}}_{g}=\frac{1}{\ell}\int^{\beta}_{0}d\tau\langle\Delta^{\dagger}_{g}(\tau)\Delta_{g}(0)\rangle
e^{i\Omega\tau},$ (7)

where $\beta$ is the inverse temperature $1/T$, $g$ denotes the symmetry
factor of the superconducting order parameter
$\Delta^{\dagger}_{g}=\sum_{\mathbf{k}}g(k)c^{\dagger}_{\mathbf{k}\downarrow}c^{\dagger}_{\mathbf{-k}\uparrow}$,
and $\tau$ represents the imaginary time. The summation over $k_{x}$, $k_{y}$,
and $k_{z}$ inside the expectation value is implicitly implied. The pair
correlation function is normalized by the number of layers ($\ell$). One can
perform explicit summation over $k_{z}$ momenta dependency within the
expectation value of the Eq. 7, resulting in

$\scalebox{1.1}{{\hbox{\chi}}}_{g}(\tau)=\frac{1}{\ell}\sum_{k_{z},k^{\prime}_{z}}\underbrace{\big{\langle}\Delta_{g}^{\dagger
k_{z}}(\tau)\Delta_{g}^{k^{\prime}_{z}}(0)\big{\rangle}}_{\scalebox{1.1}{{\hbox{\chi}}}^{k_{z},k^{\prime}_{z}}_{g}}.$
(8)

From this, we can distinguish the different channels of pairing
$\scalebox{1.1}{$\chi$}^{k_{z},k^{\prime}_{z}}_{g}$ that accounts for the
correlation between two time-ordered cooper pairs residing in $k_{z}$ and
$k^{\prime}_{z}$ momenta plane.

In the layered system three types of symmetry factor are of particular
interest: the $d_{x^{2}-y^{2}}$-wave $(\cos k_{x}-\cos k_{y})$ representing
intra-layer pairing, the $d^{(1)}_{z}$-wave $(\cos k_{z})$ for inter-layer
pairing between adjacent layers, and the $d^{(2)}_{z}$-wave $(\cos 2k_{z})$
for inter-layer pairing between the two outer planes relevant only in the
context of the trilayer and quadlayer model[42, 43, 37, 44, 45]. We primarily
focus on $d_{x^{2}-y^{2}}$-wave pairing; the results on $d^{(1)}_{z}$\- and
$d^{(2)}_{z}$-wave are not shown but discussed briefly in a later section.
Interpretation of $\scalebox{1.1}{$\chi$}_{g}$ is that a positive response
would signal an enhancement of the anomalous green’s function corresponding to
the symmetry factor $g(k)$. One can re-write the Eq. (7) in the form of the
Bethe-Salepeter equation for the particle-particle channel as

$\displaystyle\scalebox{1.1}{$\chi$}_{g}$
$\displaystyle=\frac{\bar{\scalebox{1.1}{$\chi$}}_{g}}{1-\Gamma\bar{\scalebox{1.1}{$\chi$}}_{g}},$
(9) $\displaystyle\bar{\scalebox{1.1}{$\chi$}}_{g}$
$\displaystyle=\frac{1}{\ell}\int^{\beta}_{0}d\tau\langle
g(k)^{2}c_{k\downarrow}(\tau)c^{\dagger}_{k\downarrow}(0)\rangle\langle
c_{-k\uparrow}(\tau)c^{\dagger}_{-k\uparrow}(0)\rangle,$ (10)

where uncorrelated susceptibility $\bar{\scalebox{1.1}{$\chi$}}_{g}$ is fully
dressed particle-particle bubble and $\Gamma_{g}$ is the vertex insertion
[46]. When the eigenvalue of the $\Gamma\bar{\scalebox{1.1}{$\chi$}}_{g}$
reaches unity $\scalebox{1.1}{$\chi$}_{g}$ diverges, and second-order phase
transition to a superconducting state for the symmetry factor $g$ is attained.
Since divergence to $\scalebox{1.1}{$\chi$}_{g}$ is an attribute of the vertex
$\Gamma_{g}$, the strength of the vertex component of the pair correlation,
$P_{g}$, can be utilized as an indicator of superconducting instability [40,
47], given by

$\displaystyle
P_{g}=\scalebox{1.1}{$\chi$}_{g}-\bar{\scalebox{1.1}{$\chi$}}_{g}$ (11)

We shall evaluate the expectation value of $P_{g}$ diagrammatically order-by-
order in powers of $U$,$V$, and $J$. To proceed with the expansion, we first
make a simplifying assumption that the non-interacting propagator can be
expressed in the diagonal eigenbasis of the Hamiltonian,

$\displaystyle G_{ab}=\frac{\delta_{ab}}{i\omega-\epsilon_{ab}},$ (12)

where $\delta_{ab}$ represents the Kronecker delta. In the next section, we
will make use of this property to construct our multiband perturbative
expansion symbolically.

Figure 1: (a) Schematic illustrating the assignment of initial band indices
$k_{z}$ and $k^{\prime}_{z}$ to the particle-particle diagrams corresponding
to pairing element $P^{k_{z},k^{\prime}_{z}}_{g}$. The blue bar denotes the
incoming and outgoing fermionic lines, where the symmetry factor $g(\vec{k})$
applies to the momenta. (b,c) Shows the particle-particle diagrams associated
with first and second-order coefficients we generate in the $P^{0,0}_{g}$ and
$P^{0,\pi}_{g}$ channel in $\ell=2$ system. $P^{\pi,\pi}_{d}$ channel has
identical coefficient and diagram sets as $P^{0,0}_{d}$, which can obtained by
substituting band index $0$ with $\pi$ and vice-versa.

###  II.4 Diagrammatic Expansion in powers of $U-V-J$

To each of the matrix elements
$\scalebox{1.1}{$\chi$}^{k_{z},k^{\prime}_{z}}_{g}$ we apply the ‘m’-th order
correction of the interacting Hamiltonian

$\scalebox{1.1}{{\hbox{\chi}}}^{k_{z},k^{\prime}_{z}}_{g}=\frac{(-1)^{m}}{m!\ell}\big{\langle}\mathcal{T}\Delta_{g}^{k_{z}}(\tau)\bigg{[}\int_{0}^{\beta}d\tau_{n}\prod^{m}_{n}H_{U}(\tau_{n})\bigg{]}\Delta_{g}^{k^{\prime}_{z}}(0)\big{\rangle}.$
(13)

For any given order $m$, the expression yields the product of $2m+2$ creation
and annihilation operators whose ensemble average must be evaluated. In a
generalized multi-band problem, the standard procedure is to Wick contract the
expression into a sum of all possible pairs, translating to a set of particle-
particle Feynman diagrams from which all the disconnected and topologically
indistinguishable diagrams are removed. The diagram sets are Fourier
transformed and depend on a set of momentum and frequency-conserving labels.
The Matsubara summation can be performed by repeated application of the
residue theorem. The rest of the momentum dependency is integrated over
stochastically via uniform distribution Monte Carlo sampling. While Matsubara
sum can be done readily by hand for selected low-order topologies, the task
becomes prohibitive at higher order due to the exponentially growing diagram
set and complexity of the analytic expression. Therefore, we use the recently
developed symbolic tool Algorithmic Matsubara integration (AMI), which
automates the Matsubara frequency summation using residue theorem[48, 24, 49].
The symbolic AMI tool has found success in evaluating single-particle and two-
particle correlation functions in single-band problems [50, 51, 52, 53, 21,
26].

The challenge in applying AMI to a generalized multi-band problem lies in
generating a symbolic representation of Feynman diagrams with correct band
indices and interactions assigned to them. This has been accomplished for a
generalized multi-band model using a determinantal method that is suited to
AMI.[54] However, in the determinantal methods, it is impossible to separate
the coefficients for series expansions based on $U$, $V$, and $J$ as the
interactions become intertwined during the determinant construction. This
necessitates the exhaustive computation of the entire diagram sets for each
interaction parameter set. To circumvent this, we directly generate in advance
all the distinct connected $\scalebox{1.1}{$\chi$}^{k_{z},k^{\prime}_{z}}_{g}$
particle-particle Feynman topologies, avoiding the determinant construction
for Wick’s decomposition altogether. We then exclude diagrams that are
topologically indistinguishable, do not possess vertex components, or contain
any Hartree insertions to filter out the required topologies to compute
$P^{k_{z},k^{\prime}_{z}}_{g}$ at fixed $\mu$. For every generated topology,
if one were to naively assign every possible combination of band indexes to
the propagators, it would produce numerous diagrams containing at least one or
more interactions corresponding to a zero matrix element of the $H_{U}$
interaction. These diagram sets have no contribution, so sampling and
computing them would result in a large, unnecessary computational cost. In
order to avoid this, we introduce a new iterative scheme where we assign the
band indexes based on only _non-zero_ matrix elements of interaction, ensuring
that we directly generate and compute only the contributing diagram sets for
$P^{k_{z},k^{\prime}_{z}}_{g}$.

To do so, we first diagonalize the Hamiltonian to obtain the eigenbasis that
forms our band indices. This allows us to construct the diagonal basis of our
non-interacting propagator given by Eq. 12, ensuring only a single band
indices can be assigned to each propagator. Similarly, the non-zero matrix
elements of $H_{U}$ are parameterized on the basis of these band indices. For
the matrix element of the pair correlation $P^{k_{z},k_{z}^{\prime}}_{g}$, we
already know that two time-ordered order parameters belong to band indices
$k_{z}$ and $k^{\prime}_{z}$. Diagrammatically, this pertains to assigning
$k_{z}$ and $k^{\prime}_{z}$ to the two incoming fermionic lines and two
outgoing lines, respectively, as shown by the schematics in Fig. 1(a). For the
remaining $2m-2$ unassigned fermionic lines, we iteratively assign all
possible combinations of band indices, with the selection criterion being that
the interaction sets must correspond to non-zero matrix elements of the
interaction. If an interaction not corresponding to $H_{U}$ emerges at any
stage, the sampling process is immediately halted, and that particular
combination of band indices is dismissed. We then proceed to the next possible
combination. This process is repeated until the complete diagram set with
correct band indices and non-zero interactions is generated for a given
topology.

This scheme can be automated for all topologies, quickly constructing the
entire symbolic representation of the multiband diagrammatic expansion of
$P^{k_{z},k^{\prime}_{z}}_{g}$ that is amenable to AMI. Furthermore, during
the sampling process, one can easily book-keep the number of $U$,$V$, and $J$
interactions present in the diagram sets. Consequently, diagrams can be sorted
and summed based on the number of each interaction present. This allows one to
extract the coefficient of $P^{k_{z},k^{\prime}_{z}}_{g}$ for a multi-power
series expansion in the form

$\displaystyle P^{k_{z},k^{\prime}_{z}}_{g}$
$\displaystyle=\sum_{i,j,k}a_{[i,j,k]}U^{i}V^{j}J^{k}$
$\displaystyle=a_{[1,0,0]}U+a_{[0,1,0]}V+a_{[0,0,1]}J+a_{[2,0,0]}U^{2}$
$\displaystyle\quad+a_{[0,2,0]}V^{2}+a_{[0,0,2]}J^{2}+a_{[1,1,0]}UV+a_{[1,0,1]}UJ$
$\displaystyle\quad+a_{[0,1,1]}VJ+a_{[1,1,1]}UVJ+a_{[3,0,0]}U^{3}+a_{[0,3,0]}V^{3}$
$\displaystyle\quad+a_{[0,0,3]}J^{3}\ldots$ (14)

Here $a_{[i,j,k]}$ is a function of $k_{z}$ and $k^{\prime}_{z}$ band indices
that represent the coefficient for a particular combination of the powers of
$U$, $V$, and $J$, as indicated by the indices $i$, $j$, and $k$ for the $m$th
order, where $i+j+k=m$. From a single computation of non-zero coefficients,
$P^{k_{z},k^{\prime}_{z}}_{g}$ can be evaluated for any value of $U,V,$ and
$J$ interaction strength. Moreover, these coefficients provide an
understanding of the underlying microscopic mechanisms in which pairing occurs
in a previously unexplored manner.

The total correlated pairing susceptibility is computed as the sum of all
vertex $P^{k_{z},k^{\prime}_{z}}_{g}$ channels [23] normalized by $\ell$

$P^{tot}_{g}=\frac{1}{\ell}\sum_{k_{z},k^{\prime}_{z}}P^{k_{z},k^{\prime}_{z}}_{g}.$
(15)

From the $P^{tot}_{g}$, we distinguish two different components of pairing: an
intra-plane component ($P^{\parallel}_{g}$) given by summing the diagonal
channels $P^{\parallel}_{g}=\frac{1}{\ell}\sum_{k_{z}}P^{k_{z},k_{z}}_{g}$ and
inter-plane component ($P^{\perp}_{g}$) that consists of a summation of all
the off-diagonal channels $P^{\perp}_{g}=\frac{1}{\ell}\sum_{k_{z}\neq
k^{\prime}_{z}}P^{k_{z},k^{\prime}_{z}}_{g}$ with the total being
$P^{tot}_{g}=P^{\parallel}_{g}+P^{\perp}_{g}.$

Given the form of the interaction term, in order for $P^{\perp}_{g}$ to have a
non-zero contribution in $P^{tot}_{g}$, a finite $J$ interaction is required.
For a bilayer system ($\ell=2$) consisting of band indices $k_{z}=[0,\pi]$,
one obtains

$\displaystyle P^{\parallel}_{g}(\ell=2)=(P^{0,0}_{g}+P^{\pi,\pi}_{g})/2$ (16)
$\displaystyle P^{\perp}_{g}(\ell=2)=(P^{0,\pi}_{g}+P^{\pi,0}_{g})/2$ (17)

and their summation

$P^{tot}_{g}(\ell=2)=(P^{0,0}_{g}+P^{\pi,\pi}_{g}+2P^{0,\pi}_{d})/2.$ (18)

Here we reduce the computation space of $P^{\perp}_{g}$ by taking advantage of
the symmetry $P^{k_{z},k^{\prime}_{z}}_{g}=P^{k^{\prime}_{z},k_{z}}_{g}$.
Extending this to trilayer $(\ell=3)$ case with band indices
$k_{z}=[0,\frac{\pi}{2},\pi]$ , one obtains

$\displaystyle
P^{\parallel}_{g}(\ell=3)=(P^{0,0}_{g}+P^{\frac{\pi}{2},\frac{\pi}{2}}_{g}+P^{\pi,\pi}_{g})/3$
(19) $\displaystyle
P^{\perp}_{g}(\ell=3)=2(P^{0,\frac{\pi}{2}}_{g}+P^{\frac{\pi}{2},\pi}_{g}+P^{0,\pi}_{g})/3$
(20)

and their summation

$P^{tot}_{d}(\ell=3)=(P^{0,0}_{d}+P^{\frac{\pi}{2},\frac{\pi}{2}}_{d}+P^{\pi,\pi}_{d})/3\\\
+2(P^{0,\frac{\pi}{2}}_{d}+P^{\frac{\pi}{2},\pi}_{d}+P^{0,\pi}_{d})/3.$ (21)

Eqs. (16)$\to$(21) lead us to see that in any $\ell$-layered system, there are
a total $(\ell)^{2}$ pairing channels of which there are $\ell$ diagonal
$P^{k_{z},k_{z}}_{g}$ channels that contribute to $P^{\parallel}_{g}$
component and $\ell(\ell-1)$ off-diagonal $P^{k_{z},k^{\prime}_{z}}_{g}$
channels that contribute to $P^{\perp}_{g}$ component. Upon summation and
normalization by $1/\ell$, each of $P^{k_{z},k_{z}}_{g}$ channels inside
$P^{\parallel}_{g}$ are averaged out while the $P^{k_{z},k^{\prime}_{z}}_{g}$
channels collectively contribute a net additive effect in $P^{\perp}_{d}$ with
increasing layers. Thus, $P^{\parallel}_{g}$ can be regarded as intensive and
$P^{\perp}_{g}$ as extensive property with respect number of layers.

We summarize the number of Feynman diagrams belonging to each of the
inequivalent pairing channels and $P^{tot}_{d}$ for $m$ th order in Tab. 1 and
Tab. 2 for the $\ell=2$ and $\ell=3$ system. Note that the integral
dimensionality in two-dimensional $(d=2)$ stacked lattice scales as $dm+d$.
Here, we would like to highlight the enormity of the computational expense
necessary to compute $P^{tot}_{d}$. For $\ell=2$, in the third order
$P^{tot}_{d}$, one needs to perform eight-dimensional momentum integrals over
the 276 particle-particle vertex diagrams for each data point. This complexity
balloons to 3558 diagrams with an additional two integrals when progressing to
the fourth order. This is further exacerbated by the fact that the $\ell=3$
system possesses more diagrams than the $\ell=2$ system due to increased spin-
orbital indices and interaction matrix size.

$\ell=2$
order($m$) $P^{0,0}_{g}$ $P^{0,\pi}_{g}$ $P^{\pi,\pi}_{g}$ $P^{tot}_{g}$ $1$ 1
1 1 4 $2$ 5 6 5 22 $3$ 68 70 68 276 $4$ 869 910 869 3558

Table 1: Table summarizing the total number of non-zero Feynman diagrams
involved in each inequivalent pairing channel $P^{k_{z},k^{\prime}_{z}}_{g}$
and the sum of all channels $P^{tot}_{g}$ in $\ell=2$ system (including
equivalent channels).

$\ell=3$
order($m$) $P^{0,0}_{g}$ $P^{\frac{\pi}{2},\frac{\pi}{2}}_{g}$
$P^{0,\frac{\pi}{2}}_{g}$ $P^{0,\pi}_{g}$ $P^{tot}_{g}$ $1$ 1 1 1 0 6 $2$ 5 8
6 1 39 $3$ 68 123 79 11 529 $4$ 904 1820 1100 177 7478

Table 2: Table summarizing the total number of non-zero Feynman diagrams
involved in each inequivalent pairing channel $P^{k_{z},k^{\prime}_{z}}_{d}$
and the sum of all channels $P^{tot}_{d}$ in $\ell=3$ system (including
equivalent channels) . $P^{\pi,\pi}_{d}$ and $P^{\frac{\pi}{2},\pi}_{d}$
channels have same number of diagram as the the $P^{0,0}_{d}$ and
$P^{0,\frac{\pi}{2}}_{d}$ channels respectively.

##  III Result and Discussion

In this section, we study the total correlated susceptibility $P^{tot}_{d}$
along with its intra-plane $P^{\parallel}_{d}$ and inter-plane $P^{\perp}_{d}$
components in the $d_{x^{2}-y^{2}}$ ($d$) symmetry channel for $\ell=2,3$ and
$4$ systems. We conduct a detailed analysis of the $a_{[i,j,k]}$ coefficients
for the different channels, enabling us to identify the diagram sets that
enhance or suppress two pairing components over the _U-V-J_ interaction space.
We also construct a generalized equation that speculates how
$P^{tot}_{d}(\ell)$ changes as a function of layer beyond $\ell=4$. We
primarily truncate our expansion at the third order and operate at high
temperature and weak coupling limit by setting $\beta=5$ and $U\leq 3.0$ with
$V,J<U$. While temperature plays a key role, we restrict this parameter choice
to ensure our perturbative expansion remains controlled, as demonstrated in
our prior single-band 2D Hubbard model study [26]. We examine the effect of
fourth-order corrections on selected data sets due to prohibitive
computational expense. Guided by previous studies, we fix the ratio of
$J^{\prime}/J=-0.5$, such that Cooper pairs remain in phase between the layers
[32, 37]. Here, we only study the anisotropic tunneling case with
$t_{\perp}=0.125$ and $t_{bs}=0$. We study the effect of both isotropic and
anisotropic tunneling in the supplemental materials. We find that inter-plane
tunneling strength adversely affects $P^{tot}_{d}$, and there is no note-
worthy difference between $t_{bs}$ and $t_{\perp}$. Numerous theoretical works
have extensively studied single-particle properties and two particles with
both isotropic and anisotropic tunneling, using various approaches [55, 56,
57, 58, 59, 60, 61, 32, 33, 34, 35, 36, 37, 62].

###  III.1 Bilayer result

We present the $a_{[i,j,k]}$ coefficients of the
$P^{0,0}_{d}$,$P^{\pi,\pi}_{d}$ and $P^{0,\pi}_{d}$ channels for $\ell=2$
system as a function of $\mu$ for $U,V,J=1$ up to third order in Fig. 2. The
coefficients are arranged in columns based on their relative magnitude, from
left (largest) to right (smallest). The schematic in Fig. 1(b,c) depicts all
the first and second-order diagram sets corresponding to the non-zero
coefficients for the $P^{0,0}_{d}$ and $P^{0,\pi}_{d}$ channels. It’s worth
noting that the $P^{\pi,\pi}_{d}$ channel shares identical sets of
coefficients and corresponding diagrams as the $P^{0,0}_{d}$ channel. We
remind the readers that since the offsite pair hopping term $J^{\prime}$ is
incorporated inside the onsite $J$ terms with the parameterization
$J(q)\rightarrow J[1+2J^{\prime}/J(\cos(q_{x})+\cos(q_{y}))]$, the value of
$J^{\prime}$ is fixed to the ratio $J^{\prime}=-0.5J$. The $\mu$ in the
$k_{z}=[0,\pi]$ bands is fixed such that doping in the two bands is equal.
Multi-power series expansion is performed to $a_{[i,j,k]}$ coefficients using
the Eq.II.4 to evaluate each of the $P^{k_{z},k^{\prime}_{z}}_{d}$ channel in
the _U-V-J_ interaction space. These coefficients will be utilized later to
identify diagram sets that enhance or reduce pairing in the two components.

Figure 2:  The $a_{[i,j,k]}$ coefficients of $P^{k_{z},k^{\prime}_{z}}_{d}$ as
a function of $\mu$ in the $\ell=2$ model determined by grouping particle-
particle Feynman diagrams based on the number of $U$,$V$, and $J$ interactions
present. The first row (a-c) shows intra-layer pairing in the bonding plane
$P^{0,0}_{d}$, the second row (d-f) in the anti-bonding plane, and the third
row (g-h) the off-diagonal $P^{0,\pi}_{d}$ components, representing inter-
plane scattering. Here, we set $U,V,J=1$ with $J^{\prime}/J=-0.5$.

Using the Eq. 18, Eq. 19 and Eq. 20 on the three channels of pairing, one can
obtain two pairing components $P^{\parallel}_{d}$ and $P^{\perp}_{d}$ and
their sum, $P^{tot}_{d}$. In order to understand the role interaction in
$P^{tot}_{d}$ pairing, we first plot truncated third order $P^{\parallel}_{d}$
and $P^{\perp}_{d}$ as a function of $\mu$ in Fig. 3. We systematically probe
the interaction space in the three panels by discretely varying one of $U$,
$V$, and $J$ while keeping the others fixed. Firstly, it is observed that the
peak structure resides in the negative $\mu$ region. In Fig. 3(a,d), when we
keep $V/U=0.25$ and $J/U=0.10$ fixed, we observe that as $U$ increases, the
pairing amplitude in $P^{\parallel}_{d}$ increases, while in $P^{\perp}_{d}$
it remains relatively unchanged. However, $P^{\perp}_{d}$ shows a shift in
peak structure and location, with the formation of multiple valleys
encompassing the $\mu=0$ to $\mu=-2$ region, especially noticeable in the
cases of $U=3.0$ and $U=3.5$. It’s important to note that the $U$ dependent
features are masked by the fixed ratio $J/U$, as increasing $U$ also increases
$J$, which affects $P^{\perp}_{d}$. When the absolute value of $J$ is fixed,
increasing $U$ weakens the $P^{\perp}_{d}$ component. In Fig. 3(b,d), we
examine the $V/U$ dependence of the pairing component with fixed $U=3$ and
$J/U=0.10$. The strength and peak structures of $P^{\parallel}_{d}$ and
$P^{\perp}_{d}$ are largely unaffected for the four choices of $V$. Finally,
in Fig. 3(b,d), we increase $J/U$ while keeping $U=3$ and $V/U=0.25$ fixed and
notice an overall enhancement in both $P^{\parallel}_{d}$ and $P^{\perp}_{d}$
component. While $P^{\perp}_{d}$ has no contribution at $J=0$, a significantly
larger enhancement is observed with increasing $J$ compared to the
$P^{\parallel}_{d}$ case.

Figure 3:  $P^{\parallel}_{d}$ and $P^{\perp}_{d}$ component obtained for
$\ell=2$ system by performing multi-power-series expansion on the coefficients
in Fig. 2 is plotted as functions of the chemical potential $\mu$. Panels
(a,d) vary $U/t$ with fixed $V/U=0.25$ and $J/U=0.10$, panels (b,e) vary $V/U$
with fixed $U/t=3.0$ and $J/U=0.10$ and panels (c,f) vary $J/U$ with fixed
$U/t=3.0$ and $V/U=0.25$.

The cumulative effect of $P^{\parallel}_{d}$ and $P^{\perp}_{d}$ in
interaction space is studied by plotting $P^{tot}_{d}$ as a function of $\mu$
in Fig. 4(a-c), while maintaining the same parameter selection as in Fig. 3.It
is observed that $P^{tot}_{d}$ shows a peak structure centered around
$\mu=-0.6$ in all cases. The strength of the peak is positively correlated to
the values of $U$ and $J$, where the former mainly originates from the
$P^{\parallel}_{d}$ component and the latter is influenced by both
$P^{\parallel}_{d}$ and $P^{\perp}_{d}$ components.

Given the enormity of the diagram space presented in Tab. 1, we resort to
computing the fourth order for only a single data point. There are no
conceptual hurdles evaluating beyond the fourth order but a computational
difficulty associated with evaluating many Feynman diagrams with growing
integral dimensionality that scales as $2m+2$. We apply the fourth-order
correction at $\mu=-0.6$, where the third-order $P^{tot}_{d}$ is at its peak,
and examine it as a continuous function of $U,V$, and $J$ in Fig. 4(d-f). In
Fig. 4(d), we notice that $P^{tot}_{d}$ increases with increasing $U$, with
the onset of saturation only appearing for $U>3$ in the fourth-order case.
This follows the pattern observed in the single-layer case, where the fourth
and fifth-order terms have negative contributions that suppress pairing [26].
On the other hand, increasing $V$ has a contrasting effect, where the
truncated third order slightly enhances $P^{tot}_{d}$, while the fourth-order
correction suppresses it at the large $V$ limit. In Fig. 4(e), we find a
significant increase in both $P^{\parallel}_{d}$ and $P^{\perp}_{d}$
components as $J/U$ increases. This growth is quadratic with a prominent
linear component for both the truncated third and fourth-order expansions.
This pattern can be attributed to the quadratic scaling of $P^{\parallel}_{d}$
and the linear scaling of $P^{\perp}_{d}$ with respect to $J$ (not shown),
indicating a substantial contribution primarily from first- and then second-
order diagrams. This explains why the fourth-order correction has a negligible
impact even when a large $J/U=0.4$ ratio is used.

Figure 4:  $P^{tot}_{d}$ in $\ell=2$ system obtained by summing the
$P^{\perp}_{d}$ and $P^{\parallel}_{d}$ components in Fig. 3 plotted as
function of $\mu$ (a-c). The blacked dotted line at $\mu=-0.6$ shows where
truncated 3rd order $(m=3)$ $P^{tot}_{d}$ is peaked for $U/t=3$ case. Fourth-
order correction $(m=4)$ is applied to $P^{tot}_{d}$ and plotted as a function
of (d) $U/t$, (e)$V/U$ and (f) $J/U$ at $\mu=-0.6$.

###  III.2 Analysis of coefficients in $\ell=2$

To dissect the microscopic mechanism driving bilayer pairing, we now take a
closer look at $a_{[i,j,k]}$ of all the inequivalent
$P^{k_{z},k^{\prime}_{z}}_{d}$ channels that are present in Fig. 2.
Identifying the dominant coefficient makes it possible to pinpoint the key
particle-particle diagrams and the combinations of $U,V$, and $J$ that
positively contribute to pairing. Since we utilize a $U$ interaction that is
significantly higher than $V$ and $J$, we first focus on the dominant
coefficients consisting of only $U$ interactions given by $a_{[i,0,0]}$
coefficients with $V,J=0$. For $P^{0,0}_{d}$ and $P^{\pi,\pi}_{d}$ case, the
$a_{[1,0,0]}$ coefficient consisting of a first-order ladder diagram has no
contribution. This is because the factorizability of ladder diagrams makes the
symmetry projection of $d_{x^{2}-y^{2}}$ independent of each other, resulting
in a zero contribution for $\mathcal{\mathbf{Q}}=0$ pairing. The largest
contribution originates from the second-order $a_{[2,0,0]}$ coefficient
consisting of a ladder and a crossed-interaction line topology shown in Fig.
1(b) and followed by $a_{[3,0,0]}$ coefficient consisting of the same 3rd-
order diagrams set that appears in single layer Hubbard model case [26].
Summing the $a_{[i,0,0]}$ coefficients for respective $P^{k_{z},k_{z}}_{d}$
diagonal channels would represent the single band case, and the results are
consistent with our previous single-layer study. While $a_{[2,0,0]}$ and
$a_{[3,0,0]}$ coefficients are relatively smaller, they constitute by far the
biggest contribution when expanded in powers of $U>2$.

Turning our attention to the dominant coefficients at $U,V,J=1$, off-diagonal
$P^{0,\pi}_{d}$ channel includes leading coefficients $a_{[0,0,1]}$
represented by a single first-order ladder diagram followed by $a_{[0,1,1]}$,
whose diagram sets consist of a single pair hopping interaction $J$ as shown
in Fig. 1(c). The physical interpretation of these coefficients is the pairing
contribution that originates from a Cooper pair hopping from one plane to the
adjacent plane via a single $J$ interaction. In addition, the diagonal
$P^{0,0}_{d}$ and $P^{\pi,\pi}_{d}$ channels have $a_{[0,0,2]}$ as the leading
coefficient, which consists of a single second-order ladder diagram, and
$a_{[0,1,2]}$ as subdominant coefficients with diagrams sets also involving
two $J$ interactions. These coefficients suggest pairing contribution from a
Cooper pair hopping to the nearest plane and then returning to the starting
plane via two $J$ interactions.

When comparing all the coefficients presented above, $a_{[0,0,1]}$ coefficient
in the off-diagonal $P^{0,\pi}_{d}$ channel stands out as the largest
coefficient that is roughly three times greater than the next leading
coefficient $a_{[0,0,2]}$ in $P^{\pi,\pi}_{d}$ diagonal channel. Unlike
$a_{[i,0,0]}$ coefficients expanded in $U/t=3$, we employ a $J$ value less
than unity, indicating that higher-order coefficients involving only $J$ will
make increasing lesser contribution. As a result, lower-order coefficients
make a greater contribution to $P^{tot}_{d}$ when expanded in powers of $J$.
The fact suggests that the rapid growth of $P^{tot}_{d}$ with increasing $J$
primarily originates from the $P^{\perp}_{d}$ component that is linear with
$J$ and a secondary contribution from the $P^{\parallel}_{d}$ component that
scales quadratically $(J^{2})$.

Also, it’s noteworthy that the other two subdominant $a_{[0,1,2]}$ and
$a_{[0,1,1]}$ coefficients contain a single $V$ interaction. Implying that the
presence of finite $V$ should facilitate the extent to which $P^{tot}_{d}$
grows as a function of $J$. However, when examining the remaining coefficients
involving $V$, we observe numerous cancellations among coefficients with
alternating signs, which cumulatively tunes down the enhancement in
$P^{tot}_{d}$ and introduces a large stochastic uncertainty.

Now, a contrasting picture emerges where ladder diagrams involving $U$ in
$a_{[i,0,0]}$ coefficients having zero contribution while those with
$a_{[0,0,j]}$ have the largest contribution. This is a direct consequence of
the incorporation of $J^{\prime}$ with $J$ where $2(\cos(q_{x})+\cos(q_{y}))$
form factor is applied to momentum transfer, $\mathbf{q}$, associated with $J$
vertex correction in ladder diagrams. As a result, the ladder diagrams are no
longer factorizable and yield a significant contribution to $P^{tot}_{d}$,
particularly the first-order ladder diagram. The contribution of these ladder
diagrams is pinned to the ratio and sign of $J^{\prime}/J$ employed. In the
$J^{\prime}/J=0$ case, the ladder diagrams mimic the factorizable nature of
$a_{[i,0,0]}$ ladder diagrams, leading to zero contribution. The non-ladder
diagrams results in marginal enhancement to both $P^{\parallel}_{d}$ and
$P^{\perp}_{d}$ components. With a positive $J^{\prime}/J=0.5$, the signs of
ladders diagrams inside $a_{[0,0,j]}$ for every odd $j$ is flipped, in
particular the leading $a_{[0,0,1]}$ coefficient in $P^{0,\pi}_{d}$, driving
$P^{\perp}_{d}$ repulsive with increasing $J$. Moreover, the magnitude of
$J^{\prime}/J$ utilized alters the extent of enhancement or suppression
experienced. Therefore, we argue that $J^{\prime}/J<0$ is a key requirement to
increase pairing in multilayer systems over that of a single-layer system.
Since these ladder diagrams are not part of spin and charge diagram sets, $J$
dependency of $P^{tot}_{d}$ will not be reflected on spin and charge
susceptibility. Therefore, $q=(\pi,\pi)$ antiferromagnetic spin fluctuation
that mediates the formation of $d$-wave anomalous greens self-energy is not
the mechanism behind the enhancement of $P^{\perp}_{d}$ with $J$ [63]. To
demonstrate this, we have calculated the total intra-layer staggered spin
susceptibility on the bilayer system in the supplemental materials.

In the previous section, we have discussed that increasing the value of $U$
has a detrimental effect on $P^{\perp}_{d}$ when the values of $J$ and $V$ are
fixed. But the pairing amplitude remains relatively unchanged when the ratio
of $J/U$ and $V/U$ is fixed, as shown in Fig. 3(d). This can be explained by
inspecting the $a_{[1,1,1]}$ and $a_{[2,0,1]}$ coefficients in Fig. 2(h) with
negative amplitude in the $P^{0,\pi}_{d}$ channel. When expanded in powers of
$U$, $V$, and $J$ interactions, they form a comparable negative contribution
that competes with the positive influence of $a_{[0,0,1]}$ and $a_{[0,1,1]}$
coefficients. These opposing effects cancel each other out, resulting in a
comparable $P^{\perp}_{d}$ amplitude.

###  III.3 Trilayer Result

In the bilayer case, our analysis has revealed that $a_{[0,0,2]}$ and
$a_{[0,1,2]}$ are the two leading dominant coefficients in the diagonal
channel. These coefficients are responsible for enhancing $P^{\parallel}_{d}$
with increasing $J$. While coefficients $a_{[2,0,0]}$ and $a_{[3,0,0]}$ are
small they become the leading contributor to $P^{\parallel}_{d}$ when expanded
with $U\geq 2$. Furthermore, in the off-diagonal channel for $P^{\perp}_{d}$,
the essential features are captured by first-order $a_{[0,0,1]}$ followed by
second-order $a_{[0,1,1]}$ coefficients when restricted to small values of
$V,J<1$.

In order to understand the pairing process in the trilayer model, we will
first compare the dominant coefficients $a_{[i,j,k]}$ in $\ell=3$ with the
coefficients in $\ell=2$ as a function of $\mu$. The raw coefficient for
$\ell=3$ with $U,V,J=1$, and $J^{\prime}/J=-0.5$ across all the inequivalent
pairing channels $P^{k_{z},k^{\prime}_{z}}_{g}$ is provided in the
supplemental materials. By simply comparing the dominant coefficients with the
$\ell=2$ case, we can qualitatively predict how the pairing in $\ell=3$ might
behave.

Figure 5: The contribution of the coefficient in $\ell=2,3$ to
$P^{\parallel}_{d}$ for $U=3,V=0.75,J=0.3$ as function of $\mu$ for (a)
$a_{[0,0,2]}$ with dashed line and $a_{[2,0,0]}$ with solid line (b)
$a_{[0,1,2]}$ a with dashed line and $a_{3,0,0}$ with solid line. (c)
Comparison of $P^{\parallel}_{d}$ between $\ell=2$ and $\ell=3$ as function of
$\mu$ with same choice of interaction strength.

Upon examining the $\ell=3$ coefficients in the $P^{0,0}_{d}$,
$P^{\frac{\pi}{2},\frac{\pi}{2}}_{d}$, and $P^{\pi,\pi}_{d}$ diagonal
channels, it is found that the non-zero coefficient sets are shared, including
the dominant coefficient with those in the $P^{0,0}_{d}$ and $P^{\pi,\pi}_{d}$
channels for $\ell=2$ as shown in Fig. 2(a-f). Although there are slight
variations in magnitude and peak locations due to the relative position of the
van Hove singularity to the Fermi surface (see supplemental) when these
channels are summed and averaged over $\ell$, the resulting
$P^{\parallel}_{d}$ shows no significant differences. To illustrate this, we
plot the overall contribution of leading coefficients $a_{[2,0,0]}$,
$a_{[3,0,0]}$, $a_{[0,0,2]}$, and $a_{[0,1,2]}$ to $P^{\parallel}_{d}$ at
$U=3.0$, $V=0.75$ and $J=0.3$ for the $\ell=2,3$ system in Fig. 5(a,b). It is
noted that the positive peak of $a_{[2,0,0]}$ and $a_{[3,0,0]}$ in $\ell=3$ is
slightly smaller than in the $\ell=2$ case at $U=3$, suggesting that the
$P^{tot}_{d}$ in $\ell=3$ is expected to be attenuated compared to $\ell=2$
when $U=3$ and $V,J=0$. Looking at $a_{[0,0,2]}$ and $a_{[0,1,2]}$, although
being the dominant coefficients at $U,V,J=1$, they are now reduced to
secondary contributions. A contrasting scenario is observed where
$a_{[0,0,2]}$ and $a_{[0,1,2]}$ are slightly larger in $\ell=3$ compared to
the $\ell=2$ case. Along with other coefficients providing non-negligible
contributions in the presence of finite $J$, the weakening of $a_{[2,0,0]}$
and $a_{[3,0,0]}$ coefficients in the $\ell=3$ system is counteracted,
resulting in an equivalent $P^{\parallel}_{d}$ between $\ell=2$ and $\ell=3$
at $J=0.3$, as shown in Fig. 5(c).

The adjacent-plane, off-diagonal, channels, $P^{0,\frac{\pi}{2}}_{d}$ and
$P^{\frac{\pi}{2},\pi}_{d}$, also retain the same leading coefficients
consisting of identical diagram topologies as in the $P^{0,\pi}_{d}$ in
$\ell=2$ case. To asses their impact on $P^{\perp}_{d}$, we now plot the
overall contribution of $a_{[0,0,1]}$ and $a_{[0,1,1]}$ coefficient for
$U=3.0$, $V=0.75$ and $J=0.3$ in Fig. 6(a). It is evident that these
coefficients in $\ell=3$ contribute significantly more in $P^{\perp}_{d}$ than
in $\ell=2$. This increased contribution is attributed to the presence of two
pairs of adjacent off-diagonal channels ($P^{0,\pi/2}_{d}$,$P^{\pi/2,0}_{d}$ &
$P^{\pi/2,\pi}_{d}$,$P^{\pi,\pi/2}_{d}$) containing these coefficients in
$\ell=3$ as opposed to a single pair ($P^{0,\pi}_{d}$,$P^{\pi,0}_{d}$) in
$\ell=2$. The coefficient sets in these channels have comparable pairing
amplitude for both $\ell=2$ and $\ell=3$. Consequently, when the increased
number of channels is summed and normalized by $1/\ell$, a net additive effect
is yielded.

The only difference between $\ell=2$ and $\ell=3$ lies in the presence of the
next-adjacent off-diagonal $P^{0,\pi}_{d}$ channels, whose non-zero
coefficients $a_{[0,0,2]}$ and $a_{[0,1,2]}$ have positive contributions,
albeit weaker than dominant $a_{[0,0,1]}$ and $a_{[0,1,1]}$ coefficients. Note
that $a_{[0,0,2]}$ consists of a single second-order ladder diagram while
$a_{[0,1,2]}$ has collections of multiple topologies. We plot the contribution
of these coefficients from the next-adjacent off-diagonal channels to
$P^{\perp}_{d}$ in Fig. 6(b) and show that they further help enhance
$P^{\perp}_{d}$ in $\ell=3$. Therefore, we arrive at a conclusion that in the
presence of $J$, $P^{\perp}_{d}(\ell=3)>P^{\perp}_{d}(\ell=2)$ as demonstrated
in Fig. 6(c). Nevertheless, our comparative analysis of the coefficients
indicates diagrammatic processes that enhance the pairing in the bilayer to be
similar to those in the trilayer case.

Figure 6: The contribution of the coefficient in $\ell=2,3$ to $P^{\perp}_{d}$
for $U=3,V=0.75,J=0.3$ as function of $\mu$ for (a) $a_{[0,0,1]}$ with dashed
line and $a_{[0,1,1]}$ with solid line (b) $a_{[0,1,2]}$ a with dashed line
and $a_{3,0,0}$ with solid line. Note that The coefficients in (b) are absent
in $\ell=2$ (c) Comparison of $P^{\perp}_{d}$ between $\ell=2$ and $\ell=3$ as
a function of $\mu$ with the same choice of interaction strength.

For completeness, we plot $P^{tot}_{d}$ as a function $\mu$ in Fig. 7(a-c)
with same parameter choice as the $\ell=2$ in Fig. 4(a-c). Apart from
exhibiting greater sensitivity to $J$, there are no particular notable
distinctions. We also assess the effect of fourth-order correction on a single
data point owing to the massive computational expense that is more severe in
the trilayer case as indicated by Tab. 2. We apply fourth-order correction to
$\mu=-0.6$ in Fig. 7(d-f), where the peak is located to maintain consistency
with the $\ell=2$ case. Our results indicate only one notable distinction
compared to the $\ell=2$ system. With increasing $V$, the fourth correction
exhibits a significant weakening of $P^{tot}_{d}$ where it is fully suppressed
for the maximum $V/U=0.5$ we utilize, in contrast to $P^{tot}_{d}(\ell=3)$ in
contrast to weakly dependent $V$ features in $P^{tot}_{d}(\ell=2)$. This
suggests that overlooking the $V$ interaction as it can done for bilayer cases
is not ideal, especially for large $V$. Of interesting note, for large
$J/U=0.4$, $\ell=3$ system shows an approximate seven-fold increase in
$P^{tot}_{d}$ compared to a fourfold increase observed in the $\ell=2$ system.

Using the truncated fourth order $P^{tot}_{d}$ results on $\ell=2,3$ systems,
we can establish a region of $U/t\leq 3$ and $V/U\leq 0.3$ on interaction
space in which our third-order perturbative expansion remains controlled, and
fourth-order corrections are minimum. While $J/U$ has a wider region of
validity, we restrict ourselves to $J\leq 1/6$ so that it remains realistic
and comparable to Hund’s coupling. Confining ourselves to this parameter space
allows us to study layer dependency with third-order expansion in the
proceeding section.

We have also studied the effect of $P^{tot}_{d}$ in the _U-V-J_ space when the
$\mu$ in each band is unconstrained for $\ell=2$ and $\ell=3$ in the
supplemental section. We have found that pairing is peaked at $\mu$ in
proximity to the location of van Hove’s singularity in their respective bands.
However, enhancement received by fine-tuning $\vec{\mu}$ space at $J=0$ is
superseded by enhancement due to the presence of finite $J$ for the fixed
$\mu=-0.6$ case we study here. Nevertheless, they can work in tandem to
enhance $P^{tot}_{d}$ in the layered systems.

Finally, we conclude our discussion with remarks on the inter-layer pairing
$d^{(1)}_{z}$ and $d^{(2)}_{z}$ in $\ell=2$ and $\ell=3$ system. In the
context of trilayer cuprates, inter-layer pairing alongside Cooper pair
hopping is considered a possible mechanism behind enhancement in $T_{c}$
observed experimentally [45, 64]. In our calculation (not shown), however, we
have found inter-layer pairing ($d_{z}$) to be entirely repulsive in the
truncated third order in the weak coupling _U-V-J_ space at $t_{\perp}=0.125$
and $\beta=5$. However, upon the fourth-order correction, a small region of
the attractive region appears around $U>3.25$ that is very sensitive to $V$
positively and $J$ negatively for both $P_{d^{(1)}_{z}}$ and
$P_{d^{(2)}_{z}}$. Given that intra-layer pairing is a feature of the
intermediate coupling regime, higher-order corrections above the fourth order
are necessary to attain convergence, which is beyond the scope of
computational capacity. Therefore, we leave it for future research to
determine the role of intra-layer pairing in multi-layer systems.

Figure 7:  $P^{tot}_{d}$ for $\ell=3$ obtained by summing the $P^{\perp}_{d}$
and $P^{\parallel}_{d}$ components is plotted as a function of $\mu$ (a-c).
The blacked dotted line at $\mu=-0.6$ shows where truncated 3rd order $(m=3)$
$P^{tot}_{d}$ is peaked for $U/t=3$ case. Fourth-order correction $(m=4)$ is
applied to $P^{tot}_{d}$ and plotted as a function of (d) $U/t$, (e)$V/U$ and
(f) $J/U$ at $\mu=-0.6$.

###  III.4  Quadlayer ($\ell=4$) result and beyond

Having analyzed the $\ell=2$ and $\ell=3$ pairing, extending beyond allows us
to generate a more comprehensive picture of the layer-dependent features in
$P^{tot}_{d}$. To do so, we first extend our study by stacking one additional
layer to the $\ell=3$ system, resulting in the so-called quad-layer model
($\ell=4$) system. The $P^{\parallel}_{d}$ component in $\ell=4$ system
consists of four diagonal channels arising from four momenta planes
$P^{0,0}_{d}$, $P^{\frac{\pi}{3},\frac{\pi}{3}}_{d}$,
$P^{\frac{2\pi}{3},\frac{2\pi}{3}}_{d}$ and $P^{\pi,\pi}_{d}$. The
$P^{\perp}_{d}$ comprises of three equivalent adjacent pair scattering off-
diagonal channels
$P^{0,\frac{\pi}{3}}_{d}$,$P^{\frac{\pi}{3},\frac{2\pi}{3}}_{d}$ and
$P^{\frac{2\pi}{3},\pi}_{d}$ and two inquivalent next adjacent
$P^{0,\frac{2\pi}{3}}_{d}$ and $P^{\frac{\pi}{3},\pi}_{d}$. There also exists
$P^{0,\pi}_{d}$ next-next adjacent channel where the only possible coefficient
is $a_{[0,0,3]}$ within the truncated third order. We found $a_{[0,0,3]}$ to
have a negligible contribution, especially when expanded in powers of $J^{3}$
with $J<1$. So it can safely be ignored. Summing all the channels results in

$\begin{split}P^{tot}_{d}(l=4)&=(P^{0,0}_{d}+P^{\frac{\pi}{3},\frac{\pi}{3}}_{d}+P^{2\frac{\pi}{3},2\frac{\pi}{3}}_{d}+P^{\pi,\pi}_{d})/4\\\
&\quad+(P^{0,\frac{\pi}{3}}_{d}+P^{\frac{\pi}{3},\frac{2\pi}{3}}_{d}+P^{2\frac{\pi}{3},\pi}_{d})/2\\\
&\quad+(P^{0,\frac{2\pi}{3}}_{d}+P^{\frac{\pi}{3},\pi}_{d})/2.\end{split}$
(22)

here the first term represent $P^{\parallel}_{d}$ and the last two terms
represents $P^{\perp}_{d}$.

To quantify the layer dependence, we plot the $P^{\parallel}_{d}$,
$P^{\perp}_{d}$ and $P^{tot}_{d}$ up to $\ell=4$ system as a function of $\mu$
with a fixed $U=3.0$ and $V=0.75$ for three choices of $J=0.0$, $J=0.3$ and
$J=0.5$ depicted in the three panels of Fig. 8. In the absence of $J=0$
interaction, we see that $P^{tot}_{d}$ consisting only of $P^{\parallel}_{d}$
component is only slightly attenuated with the addition of layers Fig. 8(a-c).
This suggests that without $J$ hopping interaction, the pairing process in the
layered model is qualitatively similar to those found in the single-layer
Hubbard model. What sets the layered model apart is the presence of $J$
interaction, in particular attractive nonlocal pair hopping $J^{\prime}$, as
evident in Fig. 8(f,i). When a small $J=0.3$ is utilized, enhancement received
by $P^{\parallel}_{d}$ offsets the weak attenuation observed in $J=0$ and
magnitude of $P^{\parallel}_{d}$ becomes comparable among the four $\ell$ we
study. For a larger value of $J=0.5$, we observe that the $P^{\parallel}_{d}$
for $\ell>1$ is larger and separates itself from the $\ell=1$ case. However,
among $\ell=2,3,4$, $P^{\parallel}_{d}$ remains roughly the same as in the
case of $J=0.3$. This can be attributed to diagonal channels being averaged
out by $\ell$ in calculating $P^{\parallel}_{d}$, meaning that enhancement
with $J$ is experienced equally across all $\ell>1$ systems. Thus, we conclude
that $P^{\parallel}_{d}$ plays no role in the layer dependency. On the
contrary, $P^{\perp}_{d}$ sees a substantial growth with increasing layers
that is dependent on the strength of $J$ used as portrayed in Fig. 8(e,h).
However, this growth is sublinear, indicated by the saturating tendency with
$\ell$.

Figure 8:  $P^{\parallel}_{d}$, $P^{\perp}_{d}$ and $P^{tot}_{d}$ is plotted
for up to $\ell=4$ as a function of $\mu$. With fixed $(U,V)=(3.0,0.75)$,
three choices of $J$ interaction is utilized $J=0.0$ (a-c),$J=0.30$ (d-f) and
$J=0.50$ (g-i).

Figure 9:  (a,b) averaged adjacent inter-plane scattering term
$\bar{P}^{\perp}_{\delta\ell=1}$ and averaged next-adjacent inter-plane
scattering term $\bar{P}^{\perp}_{\delta\ell=2}$ plotted as a function of
$\mu$ for up to $\ell=4$. The vertical dotted line represents the point
$\mu=-0.6$ where these two terms are roughly constant among the $\ell$. (c)
$\bar{P}^{tot}_{d}$ as a function of $\ell$ originating from the Eq. 23. The
horizontal purple dashed line shows the $\ell\to\infty$ limit.

To speculate the layer dependency of $P^{tot}_{d}$ beyond $l>4$ for finite
$J$, we first make an assumption that there is no layer dependency
$P^{\parallel}_{d}(\ell)\approx\textit{const}$ and that layer dependent
features $P^{tot}_{d}$ are determined by only the
$P^{\perp}_{d}(\ell)\approx(P_{d}^{tot}(\ell)-P_{d}^{\parallel})$ component.
This implies that Cooper pair hopping between the planes via $J$ and
$J^{\prime}$ interactions are the only physical processes that give rise to
layer resolved features. This enables us to approximate
$P^{\perp}_{d}\approx\Delta P_{d}(\ell)$, where $\Delta P_{d}(\ell)$ is a
change in pairing as a function of $\ell$. Our prior analysis of the
coefficient in $\ell=3$ showed that different $P^{k_{z},k^{\prime}_{z}}_{d}$
scattering processes decay in amplitude drastically based on their $k_{z}$
separation in momenta planes. Scattering to the adjacent plane is dominant;
the next adjacent plane is subdominant, and beyond, the contributions are
negligible for small $J\leq 0.5$. However, not all the same $k_{z}$ separation
in $P^{k_{z},k^{\prime}_{z}}_{d}$ channels have the equivalent magnitudes and
peak locations. Should these channels be constant on average and knowing that
there are $2(\ell-1)$ and $2(\ell-2)$ adjacent and next adjacent planes, a
generalized equation that estimates $\Delta P_{d}(\ell)$ in the large $\ell$
limit can be formulated. To do this, we introduce new variables
$\bar{P}^{\perp}_{\delta\ell=1}$ and $\bar{P}^{\perp}_{\delta\ell=2}$ in Fig.
9(a,b) defined to be the average adjacent $(\delta\ell=1)$ and next adjacent
($\delta\ell=2$) pair scattering channel respectively. For instance, in
$\ell=4$, we take the average of three adjacent off-diagonal channels present
$\bar{P}^{\perp}_{\delta\ell=1}(\ell=4)=(P^{0,\pi/3}_{d}+P^{\pi/3,2\pi/3}_{d}+P^{2\pi/3,\pi}_{d})/3$
and average of two adjacent off-diagonal channel
$\bar{P}^{\perp}_{\delta\ell=2}(\ell=4)=(P^{0,2\pi/3}_{d}+P^{\pi/3,\pi}_{d})/2$.
In Fig. 9(a,b), we plot $\bar{P}^{\perp}_{\delta\ell=1}$ and
$\bar{P}^{\perp}_{\delta\ell=2}$ for $\ell=2,3,4$ systems and find existence
of wide range of $\mu$ where the curves overlap within the errorbar. This
suggests that two variables can be approximated as a constant value. This
allows us to formulate an equation as a combinatorics problem that estimates
the change in pairing $\Delta P_{d}(\ell)$ as a function of layer, given by

$\Delta
P_{d}(\ell)=2\left[\bar{P}^{\perp}_{\delta\ell=1}\big{(}1-\frac{1}{\ell}\big{)}_{\ell>1}+\bar{P}^{\perp}_{\delta\ell=2}\big{(}1-\frac{2}{\ell}\big{)}_{\ell>2}\right].$
(23)

In Fig. 9(c), we plot Eq.23 as a function of $\ell$ at $U=3$, $V=0.75$ and
$J=0.3$. Here we utilize $\bar{P}^{\perp}_{\delta\ell=1}$ and
$\bar{P}^{\perp}_{\delta\ell=2}$ at $\mu=-0.6$ where the two variables are
roughly equal as indicated by the dotted line in Fig. 9(a,b). As a reference
to readers, we also plot results for $P^{\perp}_{d}$ for $\ell=1,2,3,4$. We
emphasize here that this equation should not be treated as a fitting to the
data. But rather to capture how paring is expected to change as a function of
$\ell>4$. Our model indicates that $\Delta P_{d}$ grows as a function of the
layer sublinearly before saturating at the 3D limit where
$\lim_{\ell\to\infty}\Delta
P_{d}(\ell)=2\left[\bar{P}^{\perp}_{\delta\ell=1}+\bar{P}^{\perp}_{\delta\ell=2}\right]$.
Therefore, one cannot attain a continuous phase transition at a fixed
temperature just by stacking layers to the 3D limit despite the additive
effect of off-diagonal pairing channels for the parameter space we explored.
Nevertheless, the Eq.23 is reminiscent of the $\Delta
T_{c}(\ell)=\textit{const}(1-1/\ell)$ relation, derived by Legget from
interplane Coulombic interaction in identically doped layered Cuprates [65].
Nishiguchi _et al._ also reported a sublinear increase in $T_{c}$ using
fluctuation exchange approximation (FLEX) up to $\ell=3$ system by considering
pair hopping of various configurations. Chakravarty also derived a similar
relation via the inter-layer tunneling mechanism [66]. This monotonic growth
was in contradiction to experimental studies, where $T_{c}$ increases with the
number of the $CuO$ layers in the unit cell up to $\ell=3$, then gradually
decreases beyond $\ell>3$ and eventually saturating[67, 68]. We have observed
that, similar to the way in which attractive nonlocal interactions between
nearest neighbors enhance pairing in a single-layer model [20, 22, 69, 70],
our model extends this phenomenon to multi-layer systems. Our results indicate
an attractive nonlocal pair hopping interaction plays a critical role in
enhancing pairing fluctuations through layer stacking within the parameter
space we have explored.

###  III.5 Summary and Outlook

In this work, we have introduced a new symbolic technique to iteratively
assign the band indexes to topologies to generate complete symbolic
representations of diagram sets responsible for
$P^{k_{z},{k^{\prime}_{z}}}_{d}$. We resolved Matsubara frequency summations
analytically using algorithmic Matsubara integration and resolve momentum
integrations stochastically. For each $P^{k_{z},{k^{\prime}_{z}}}_{d}$, we
categorized the diagram sets based on interaction to obtain a set of
coefficients in the multi-power-series expansion giving full access to weak-
coupling $U,V,J$ space. We have summed the diagonal and off-diagonal channels
separately to obtain two independent components of pairing,
$P^{\parallel}_{d}$ and $P^{\perp}_{d}$ with total physical pairing
susceptibility being $P^{tot}_{d}=P^{\parallel}_{d}+P^{tot}_{d}$. We have
probed $P^{\parallel}_{d}$, $P^{\perp}_{d}$ and $P^{tot}_{d}$ as a function of
$\mu$ for up to an $\ell=4$ system in the presence of local $U,V,J$, and non-
local $J^{\prime}/J=-0.5$ interaction space.

From the analysis of coefficients in $\ell=2$ and $\ell=3$ diagonal channels,
we have determined that coefficients $a_{[2,0,0]}$ and $a_{[3,0,0]}$
consisting of only Hubbard interaction are the dominant contributor to
$P^{\parallel}_{d}$ at $U=3$. The coefficients $a_{[0,0,2]}$ and $a_{[0,1,2]}$
facilitate the enhancement in $P^{\parallel}_{d}$ component that scaling
quadratically with $J$. When $J=0$, we showed that $P^{\parallel}_{d}$ is
weakly attenuated with $\ell$ but qualitative features remain similar to
$\ell=1$ system. For a finite $J$, the enhancement received across the layer
is the same and exhibits weak to no layer-resolved features beyond $\ell=1$
for $P^{\parallel}_{d}$. This is because the number of diagonal channels is
always equal to $\ell$, so all the channels are averaged out when the
normalization factor of $1/\ell$ is applied.

The contribution of $P^{\perp}_{d}$ in small $J/U\leq 1/6$ (comparable to
Hund’s coupling) predominately comes from dominant $a_{[0,0,1]}$ and
$a_{[0,0,2]}$ coefficients in adjacent and next adjacent off-diagonal channels
consisting of the first order and second order ladder diagram respectively.
Scattering at and beyond $J^{3}$ is negligible, and we see that an attractive
non-local interaction $J^{\prime}/J$ is a crucial requirement for the
coefficients to be dominant. The channels possessing these coefficients grow
super-linearly with $\ell$, yielding a net additive effect even when
normalized by $1/\ell$, and are thus responsible for layer-dependent features
in pairing. However, $P^{\perp}_{d}$ grows sublinearly as a function for up to
$\ell=4$. To speculate beyond $\ell=4$, we use the observation that the
averaged adjacent $\bar{P}^{\perp}_{\delta\ell=1}$ and next adjacent
$\bar{P}^{\perp}_{\delta\ell=2}$ can be approximated as constants at
$\mu=-0.6$. Knowing the number of $\bar{P}^{\perp}_{\delta\ell=1}$ and
$\bar{P}^{\perp}_{\delta\ell=2}$ elements present for each $\ell$ system, one
can formulate a generalized equation for $\Delta P_{d}(\ell)$ that predicts
pairing trends beyond $\ell=4$. This model indicates a well-defined 3D limit
where pairing saturates.

Although this behavior does not reflect $T_{c}(\ell)$ dependence of Cuprates,
our results unequivocally suggest that both local and non-local pair hopping,
even when small, plays a significant role in promoting ($J^{\prime}\leq 0$) or
suppressing pairing ($J^{\prime}>0$) in coupled 2D systems in weak coupling
limit.

In the supplemental section, we have also calculated $q=(\pi,\pi)$ spin
susceptibility for $\ell=2$ and found little to no dependency in $V,J$, and
$J^{\prime}$ interaction space. From this, we concluded that spin fluctuations
do not mediate the layer-resolved features in $P^{\perp}_{d}$.

Finally, we emphasize that our approach has potential applications in
understanding pairing processes in any material that has a low-energy multi-
band effective model with any arbitrary long-range interaction. In particular,
our methodologies can be readily extended to study single and two-particle
properties in more realistic systems like the single layer three-band Emery
model for cuprates or Sr2RuO4, enabling us to identify key diagrammatic
processes and coefficients within weak coupling regimes [71, 72, 73, 74].

##  IV Acknowledgement

We acknowledge the support of the Natural Sciences and Engineering Research
Council of Canada (NSERC) RGPIN-2022-03882. Our codes make use of the
ALPSCore[75, 76] and AMI libraries[48, 49].

## References

  * LeBlanc  _et al._ [2015] J. P. F. LeBlanc, A. E. Antipov, F. Becca, I. W. Bulik, G. K.-L. Chan, C.-M. Chung, Y. Deng, M. Ferrero, T. M. Henderson, C. A. Jiménez-Hoyos, E. Kozik, X.-W. Liu, A. J. Millis, N. V. Prokof’ev, M. Qin, G. E. Scuseria, H. Shi, B. V. Svistunov, L. F. Tocchio, I. S. Tupitsyn, S. R. White, S. Zhang, B.-X. Zheng, Z. Zhu, and E. Gull (Simons Collaboration on the Many-Electron Problem), Phys. Rev. X 5, 041041 (2015).
  * Motta  _et al._ [2017] M. Motta, D. M. Ceperley, G. K.-L. Chan, J. A. Gomez, E. Gull, S. Guo, C. A. Jiménez-Hoyos, T. N. Lan, J. Li, F. Ma, A. J. Millis, N. V. Prokof’ev, U. Ray, G. E. Scuseria, S. Sorella, E. M. Stoudenmire, Q. Sun, I. S. Tupitsyn, S. R. White, D. Zgid, and S. Zhang (Simons Collaboration on the Many-Electron Problem), Phys. Rev. X 7, 031059 (2017).
  * Qin  _et al._ [2020] M. Qin, C.-M. Chung, H. Shi, E. Vitali, C. Hubig, U. Schollwöck, S. R. White, and S. Zhang (Simons Collaboration on the Many-Electron Problem), Phys. Rev. X 10, 031016 (2020).
  * Chung  _et al._ [2020] C.-M. Chung, M. Qin, S. Zhang, U. Schollwöck, and S. R. White (The Simons Collaboration on the Many-Electron Problem), Phys. Rev. B 102, 041106 (2020).
  * Williams  _et al._ [2020] K. T. Williams, Y. Yao, J. Li, L. Chen, H. Shi, M. Motta, C. Niu, U. Ray, S. Guo, R. J. Anderson, J. Li, L. N. Tran, C.-N. Yeh, B. Mussard, S. Sharma, F. Bruneval, M. van Schilfgaarde, G. H. Booth, G. K.-L. Chan, S. Zhang, E. Gull, D. Zgid, A. Millis, C. J. Umrigar, and L. K. Wagner (Simons Collaboration on the Many-Electron Problem), Phys. Rev. X 10, 011041 (2020).
  * Qin  _et al._ [2022] M. Qin, T. Schäfer, S. Andergassen, P. Corboz, and E. Gull, Annual Review of Condensed Matter Physics 13, 275 (2022).
  * Schäfer  _et al._ [2021] T. Schäfer, N. Wentzell, F. Šimkovic, Y.-Y. He, C. Hille, M. Klett, C. J. Eckhardt, B. Arzhang, V. Harkov, F. m. c.-M. Le Régent, A. Kirsch, Y. Wang, A. J. Kim, E. Kozik, E. A. Stepanov, A. Kauch, S. Andergassen, P. Hansmann, D. Rohe, Y. M. Vilk, J. P. F. LeBlanc, S. Zhang, A.-M. S. Tremblay, M. Ferrero, O. Parcollet, and A. Georges, Phys. Rev. X 11, 011058 (2021).
  * Coldea  _et al._ [2001] R. Coldea, S. M. Hayden, G. Aeppli, T. G. Perring, C. D. Frost, T. E. Mason, S.-W. Cheong, and Z. Fisk, Phys. Rev. Lett. 86, 5377 (2001).
  * LeBlanc  _et al._ [2019] J. P. F. LeBlanc, S. Li, X. Chen, R. Levy, A. E. Antipov, A. J. Millis, and E. Gull, Phys. Rev. B 100, 075123 (2019).
  * Rybicki  _et al._ [2016] D. Rybicki, M. Jurkutat, S. Reichardt, C. Kapusta, and J. Haase, Nature Communications 7, 11413 (2016).
  * Pelc  _et al._ [2020] D. Pelc, M. J. Veit, C. J. Dorow, Y. Ge, N. Barišić, and M. Greven, Phys. Rev. B 102, 075114 (2020).
  * Terletska  _et al._ [2021] H. Terletska, S. Iskakov, T. Maier, and E. Gull, Phys. Rev. B 104, 085129 (2021).
  * Gull  _et al._ [2013] E. Gull, O. Parcollet, and A. J. Millis, Phys. Rev. Lett. 110, 216405 (2013).
  * Sordi  _et al._ [2011] G. Sordi, K. Haule, and A.-M. S. Tremblay, Phys. Rev. B 84, 075161 (2011).
  * Zheng and Chan [2016] B.-X. Zheng and G. K.-L. Chan, Phys. Rev. B 93, 035126 (2016).
  * Deng  _et al._ [2015] Y. Deng, E. Kozik, N. V. Prokof’ev, and B. V. Svistunov, EPL (Europhysics Letters) 110, 57001 (2015).
  * Jiang  _et al._ [2024] Y.-F. Jiang, T. P. Devereaux, and H.-C. Jiang, Phys. Rev. B 109, 085121 (2024).
  * Jiang  _et al._ [2020] Y.-F. Jiang, J. Zaanen, T. P. Devereaux, and H.-C. Jiang, Phys. Rev. Res. 2, 033073 (2020).
  * Tang  _et al._ [2024] T. Tang, D. Jost, B. Moritz, and T. P. Devereaux, arXiv:2405.11445 (2024).
  * Jiang and Devereaux [2019] H.-C. Jiang and T. P. Devereaux, Science 365, 1424 (2019), https://www.science.org/doi/pdf/10.1126/science.aal5304 .
  * Gazizova and LeBlanc [2023] D. Gazizova and J. P. F. LeBlanc, Phys. Rev. B 108, 165149 (2023).
  * Peng  _et al._ [2023] C. Peng, Y. Wang, J. Wen, Y. S. Lee, T. P. Devereaux, and H.-C. Jiang, Phys. Rev. B 107, L201102 (2023).
  * Mai  _et al._ [2021] P. Mai, G. Balduzzi, S. Johnston, and T. A. Maier, Phys. Rev. B 103, 144514 (2021).
  * Taheridehkordi  _et al._ [2019] A. Taheridehkordi, S. H. Curnoe, and J. P. F. LeBlanc, Phys. Rev. B 99, 035120 (2019).
  * Taheridehkordi  _et al._ [2020a] A. Taheridehkordi, S. H. Curnoe, and J. P. F. LeBlanc, Phys. Rev. B 101, 125109 (2020a).
  * Farid  _et al._ [2023] R. Farid, M. Grandadam, and J. P. F. LeBlanc, Phys. Rev. B 107, 195138 (2023).
  * Hubbard [1963] J. Hubbard, Proc. R. Soc. London, Ser. A 276, 238 (1963).
  * Andersen  _et al._ [1995] O. Andersen, A. Liechtenstein, O. Jepsen, and F. Paulsen, Journal of Physics and Chemistry of Solids 56, 1573 (1995), proceedings of the Conference on Spectroscopies in Novel Superconductors.
  * Andersen  _et al._ [1994] O. K. Andersen, O. Jepsen, A. I. Liechtenstein, and I. I. Mazin, Phys. Rev. B 49, 4145 (1994).
  * Chakravarty  _et al._ [1993] S. Chakravarty, A. Sudbø, P. W. Anderson, and S. Strong, Science 261, 337 (1993).
  * Markiewicz  _et al._ [2005] R. S. Markiewicz, S. Sahrakorpi, M. Lindroos, H. Lin, and A. Bansil, Phys. Rev. B 72, 054519 (2005).
  * Nishiguchi  _et al._ [2013] K. Nishiguchi, K. Kuroki, R. Arita, T. Oka, and H. Aoki, Phys. Rev. B 88, 014509 (2013).
  * Gan  _et al._ [2008] J. Y. Gan, M. Mori, T. K. Lee, and S. Maekawa, Phys. Rev. B 78, 094504 (2008).
  * Nishiguchi  _et al._ [2017] K. Nishiguchi, S. Teranishi, and K. Kusakabe, Journal of the Physical Society of Japan 86, 084707 (2017), https://doi.org/10.7566/JPSJ.86.084707 .
  * Nishiguchi  _et al._ [2018] K. Nishiguchi, S. Teranishi, K. Kusakabe, and H. Aoki, Phys. Rev. B 98, 174508 (2018).
  * Iwano and Yamaji [2022] A. Iwano and Y. Yamaji, Journal of the Physical Society of Japan 91, 094702 (2022), https://doi.org/10.7566/JPSJ.91.094702 .
  * Zegrodnik and Spałek [2017] M. Zegrodnik and J. Spałek, Phys. Rev. B 95, 024507 (2017).
  * Kusakabe [2009] K. Kusakabe, Journal of the Physical Society of Japan 78, 114716 (2009), https://doi.org/10.1143/JPSJ.78.114716 .
  * Kusakabe [2012] K. Kusakabe, Journal of Physics and Chemistry of Solids 73, 1546 (2012), 16th International Symposium on Intercalation Compounds (ISIC 16).
  * Chen  _et al._ [2015] X. Chen, J. P. F. LeBlanc, and E. Gull, Phys. Rev. Lett. 115, 116402 (2015).
  * Rohringer  _et al._ [2012] G. Rohringer, A. Valli, and A. Toschi, Phys. Rev. B 86, 125114 (2012), arXiv:1202.2796 [cond-mat.str-el] .
  * Sigrist and Ueda [1991] M. Sigrist and K. Ueda, Rev. Mod. Phys. 63, 239 (1991).
  * Bulut  _et al._ [1992] N. Bulut, D. J. Scalapino, and R. T. Scalettar, Phys. Rev. B 45, 5577 (1992).
  * Tsuei and Kirtley [2000] C. C. Tsuei and J. R. Kirtley, Rev. Mod. Phys. 72, 969 (2000).
  * Luo  _et al._ [2023] X. Luo, H. Chen, Y. Li, Q. Gao, C. Yin, H. Yan, T. Miao, H. Luo, Y. Shu, Y. Chen, C. Lin, S. Zhang, Z. Wang, F. Zhang, F. Yang, Q. Peng, G. Liu, L. Zhao, Z. Xu, T. Xiang, and X. J. Zhou, Nature Physics 19, 1841 (2023).
  * Salpeter and Bethe [1951] E. E. Salpeter and H. A. Bethe, Phys. Rev. 84, 1232 (1951).
  * White  _et al._ [1989] S. R. White, D. J. Scalapino, R. L. Sugar, N. E. Bickers, and R. T. Scalettar, Phys. Rev. B 39, 839 (1989).
  * Elazab  _et al._ [2022] H. Elazab, B. McNiven, and J. LeBlanc, Computer Physics Communications 280, 108469 (2022).
  * Burke and LeBlanc [2023] M. D. Burke and J. P. F. LeBlanc, arXiv , 2311.17189 (2023).
  * McNiven  _et al._ [2021] B. D. E. McNiven, G. T. Andrews, and J. P. F. LeBlanc, Phys. Rev. B 104, 125114 (2021).
  * McNiven  _et al._ [2022] B. D. E. McNiven, H. Terletska, G. T. Andrews, and J. P. F. LeBlanc, Phys. Rev. B 106, 035145 (2022).
  * Taheridehkordi  _et al._ [2020b] A. Taheridehkordi, S. H. Curnoe, and J. P. F. LeBlanc, Phys. Rev. B 102, 045115 (2020b).
  * Burke  _et al._ [2023] M. D. Burke, M. Grandadam, and J. P. F. LeBlanc, Phys. Rev. B 107, 115151 (2023).
  * Assi and LeBlanc [2024] I. Assi and J. P. F. LeBlanc, Phys. Rev. B 109, 125143 (2024).
  * Maier and Scalapino [2011] T. A. Maier and D. J. Scalapino, Phys. Rev. B 84, 180513 (2011).
  * Karakuzu  _et al._ [2021] S. Karakuzu, S. Johnston, and T. A. Maier, Phys. Rev. B 104, 245109 (2021).
  * Kato and Kuroki [2020] D. Kato and K. Kuroki, Phys. Rev. Res. 2, 023156 (2020).
  * Caplan and Orgad [2023] Y. Caplan and D. Orgad, Phys. Rev. B 108, 165131 (2023).
  * Bouadim  _et al._ [2008] K. Bouadim, G. G. Batrouni, F. Hébert, and R. T. Scalettar, Phys. Rev. B 77, 144527 (2008).
  * Okamoto and Maier [2008] S. Okamoto and T. A. Maier, Phys. Rev. Lett. 101, 156401 (2008).
  * Goldberger  _et al._ [2024] D. Goldberger, Y. Fridman, E. Gull, E. Eidelstein, and G. Cohen, Phys. Rev. B 109, 085133 (2024).
  * Eremin  _et al._ [2007] I. Eremin, D. K. Morr, A. V. Chubukov, and K. Bennemann, Phys. Rev. B 75, 184534 (2007).
  * Dong  _et al._ [2022] X. Dong, E. Gull, and A. Millis, Nature Physics 18 (2022), 10.1038/s41567-022-01710-z.
  * Ideta  _et al._ [2021] S. Ideta, S. Johnston, T. Yoshida, K. Tanaka, M. Mori, H. Anzai, A. Ino, M. Arita, H. Namatame, M. Taniguchi, S. Ishida, K. Takashima, K. M. Kojima, T. P. Devereaux, S. Uchida, and A. Fujimori, Phys. Rev. Lett. 127, 217004 (2021).
  * Leggett [1999] A. J. Leggett, Phys. Rev. Lett. 83, 392 (1999).
  * Chakravarty [1998] S. Chakravarty, The European Physical Journal B - Condensed Matter and Complex Systems 5, 337 (1998).
  * Schilling  _et al._ [1993] A. Schilling, M. Cantoni, J. D. Guo, and H. R. Ott, Nature 363, 56 (1993).
  * Iyo  _et al._ [2007] A. Iyo, Y. Tanaka, H. Kito, Y. Kodama, P. M. Shirage, D. D. Shivagan, H. Matsuhata, K. Tokiwa, and T. Watanabe, Journal of the Physical Society of Japan 76, 094711 (2007), https://doi.org/10.1143/JPSJ.76.094711 .
  * Chen  _et al._ [2021] Z. Chen, Y. Wang, S. N. Rebec, T. Jia, M. Hashimoto, D. Lu, B. Moritz, R. G. Moore, T. P. Devereaux, and Z.-X. Shen, Science 373, 1235 (2021), https://www.science.org/doi/pdf/10.1126/science.abf5174 .
  * Jiang [2022] M. Jiang, Phys. Rev. B 105, 024510 (2022).
  * Rømer  _et al._ [2021] A. T. Rømer, P. J. Hirschfeld, and B. M. Andersen, Phys. Rev. B 104, 064507 (2021).
  * Yuan  _et al._ [2023] A. C. Yuan, E. Berg, and S. A. Kivelson, Phys. Rev. B 108, 014502 (2023).
  * Jiang [2023] H.-C. Jiang, Phys. Rev. B 107, 214504 (2023).
  * Kowalski  _et al._ [2021] N. Kowalski, S. S. Dash, P. Sémon, D. Sénéchal, and A.-M. Tremblay, Proceedings of the National Academy of Sciences 118, e2106476118 (2021), https://www.pnas.org/doi/pdf/10.1073/pnas.2106476118 .
  * Gaenko  _et al._ [2016] A. Gaenko, E. Gull, A. E. Antipov, L. Gamper, G. Carcassi, J. Paki, R. Levy, M. Dolfi, J. Greitemann, and J. P. LeBlanc, “Alpscore: Version 0.5.4,”  (2016).
  * Wallerberger  _et al._ [2018] M. Wallerberger, S. Iskakov, A. Gaenko, J. Kleinhenz, I. Krivenko, R. Levy, J. Li, H. Shinaoka, S. Todo, T. Chen, X. Chen, J. P. F. LeBlanc, J. E. Paki, H. Terletska, M. Troyer, and E. Gull, arXiv:1811.08331 (2018).