#  The Kirkwood-Riseman Polymer Model of Polymer Dynamics is Qualitatively
Correct

George D. J. Phillies<EMAIL_ADDRESS>508-754-1859 Department of Physics,
Worcester Polytechnic Institute,Worcester, MA 01609

###### Abstract

We use Brownian dynamics to show: For an isolated polymer coil, the Kirkwood-
Riseman model for chain motion is qualitatively correct. The Rouse model for
chain motion is qualitatively incorrect. The models are qualitatively
different. Kirkwood and Riseman say polymer coils in a shear field perform
whole-body rotation; in the Rouse model rotation does not occur. Our
simulations demonstrate that _in shear flow_ : Polymer coils rotate. Rouse
modes are cross-correlated. The amplitudes and relaxation rates of Rouse modes
depend on the shear rate. Rouse’s calculation only refers to a polymer coil in
a quiescent fluid, its application to a polymer coil undergoing shear being
invalid.

###### pacs:

83.80.Rs,66.10.cg,66.30.hk,83.86.Hf

##  I Introduction

Seven decades ago, Kirkwood and Risemankirkwood1948a , Rouserouse1953a , and
Zimmzimm1956a  advanced simple, seemingly transparent models for the dynamics
of dilute polymers in solution. A particular focus of their models was a
calculation of the polymeric contribution to a solution’s viscosity. The
models were similar in that each approximated a polymer coil as a series of
hydrodynamically active points (”beads”) held together by hydrodynamically
inert connectors (”springs”).

Here the similarity between the models ends. As is not uniformly recognized,
the Kirkwood-Riseman and Rouse-Zimm models give entirely contradictory
descriptions of polymer dynamics. In the Rouse and Zimm models, the connectors
(”springs”) between the beads create Hooke’s-law restoring forces that pull
the beads together. Thermal fluctuations in the solvent create random forces
on the polymer beads, tending on the average to drive the beads apart. The
competition between the spring and thermal forces determines the size of a
polymer coil. Rouse-Zimm polymer coils can translate; all other bead motions
are described by internal modes in which the relative positions of the beads
change. These internal modes are, in the Rouse-Zimm picture, responsible for
the polymer coil’s contribution to the viscosity.

In contrast, in the Kirkwood-Riseman model, the dynamics of a polymer coil are
approximated as comprising uniform translation and whole-body rotation. The
distances between polymer beads are fixed at their average values; within the
Kirkwood-Riseman model the distances between pairs of beads do not change with
time. According to Kirkwood and Riseman, the dominant contribution to a
polymer’s intrinsic viscosity arises from whole-body rotation. While Kirkwood
and Riseman recognized that polymer beads can move with respect to each other,
within their model these internal motions were neglected.

Many experimental aspects of polymer dynamics are described by the Rouse and
Zimm models. However, there does not appear to have been an equivalent effort
to test experiment against the Kirkwood-Riseman model, so there is no basis
for proposing that experiment favors the Rouse-Zimm model over the Kirkwood-
Riseman model.

The objective of this paper is to present a simulational test of the
predictions of the Rouse and Kirkwood-Riseman models. We actually advance
beyond the approximations of Kirkwood and Riseman’s model, by including the
fluctuations in bead positions that they neglected. We employ a wider range of
computational diagnostics than has sometimes been used in the past to
interpret polymer dynamics. In particular, we ask whether or not a chain in a
shear field is rotating. We ask if fluctuations in Rouse mode amplitudes are
cross-correlated, and if Rouse mode amplitudes or relaxation rates depend on
the rate of the applied shear.

The following Section of this paper outlines salient features of the Rouse and
Kirkwood-Riseman models. A further Section describes our simulation
procedures, including the physical quantities that we calculated. We then
outline our major results, revealing the relative validities of the Rouse and
Kirkwood-Riseman models as ttreatments of polymer dynamics during a
rheological experiment. Conclusions are presented. To anticipate our results,
we show that Kirkwood and Riseman were correct, and Rouse was incorrect.

##  II Rouse and Kirkwood-Riseman Models

This section present aspects of the Rouserouse1953a  and Kirkwood-
Risemankirkwood1948a  models. We begin with the more familiar Rouse model, and
then consider the Kirkwood-Riseman model.

Rouse’s original treatment was quite involved. As is often the case with novel
theoretical results, as time advances the key aspects of the calculation are
abstracted from the original structure. The presentation of Doi and
Edwardsdoi1986a  and the more extended development by Paddingpadding2005a  are
followed here. The Rouse model describes an isolated polymer in a solvent. The
polymer is approximated as a line of $N$ beads, each linked to the next by a
springlike connector. The bead positions are denoted
$(\bm{R}_{1},\bm{R}_{2},\ldots\bm{R}_{N})$. The beads are points having no
excluded volume; they are all free to move with respect to each other. Each
bead has a hydrodynamic drag coefficient $f$. The connectors do not interact
with the solvent.

The strength of the connectors is determined by the Gaussian statistics that
describe the shape of a random-walk polymer coil. In Rouse’s model, each bead
represents some substantial number of monomers. The distance along the polymer
chain from each bead to the next is sufficiently large that the bead-bead
distances $r_{i,i+1}=|\bm{R}_{i+1}-\bm{R}_{i}|$ have Gaussian distributions
$P(r_{i,i+1})\sim\exp(-\alpha r_{i,i+1}^{2})$.

Rouse implicitly explains that for each statistico-mechanical distribution
function $P(r_{ij})$, there is a corresponding potential of average force
$W(r_{ij})$, namely

$W(r_{ij})=-k_{B}T\ln(P(r_{ij})).$ (1)

Here $k_{B}$ is Boltzmann’s constant and T is the absolute temperature. The
potential of average force gives the average force between two adjoining beads
that are a distance $r_{ij}$ apart.

The calculations here use Rouse’s original potential of average force

$W(r_{ij})=\frac{1}{2}kr_{ij}^{2},$ (2)

for two beads $i$ and $j$ that adjoin along the polymer chain. In Rouse’s
model the force constant $k$ is determined by the mean-square bead separation
$b$, namely

$k=\frac{3k_{B}T}{b^{2}}.$ (3)

It is possible to use considerably more sophisticated forms for the potential
energy of the polymer chain. Note, for example, work of Tsalikis, et
al.tsalikis2017a , Perez-Aparicio, et al.perez2011a , and Kalathi, et
al.kalathi2014a . The objective here is to test the Rouse and Kirkwood-Riseman
models. To test the Rouse model, we must use Rouse’s potential. The Kirkwood-
Riseman model specifies only average interbead vectors, and does not invoke a
particular form for the interbead potential energy.

In the original Rouse model, a bead $i$ was also subject to a thermal force
$\mathbf{\cal{F}}_{i}(t)$ due to fluctuations in the solvent. Hydrodynamic
interactions between beads, and correlations between the thermal forces on
different beads, were neglected in the Rouse model but included in the
elsewise-similar Zimm model.

We can now write the equations of motion – Newton’s second law – for each bead
of the Rouse model. The drag force on each bead is large. On the time scales
of interest bead motions are massively overdamped. Bead inertia is therefore
neglected. If bead inertia vanishes, the total force on each bead must also
vanish. The direct forces (the spring forces) on each bead must therefore
cancel the hydrodynamic forces. The equations of motion for beads other than
the two end beads (beads 1 and $N$) are then

$f\frac{d\bm{R}_{i}}{dt}=-k(2\bm{R}_{i}-\bm{R}_{i-1}-\bm{R}_{i+1})+\mathbf{\cal{F}}_{i}(t),$
(4)

while for the first and last beads in the chain one has

$f\frac{d\bm{R}_{1}}{dt}=-k(\bm{R}_{1}-\bm{R}_{2})+\mathbf{\cal{F}}_{1}(t),$
(5)

and

$f\frac{d\bm{R}_{N}}{dt}=-k(\bm{R}_{N}-\bm{R}_{N-1})+\mathbf{\cal{F}}_{N}(t).$
(6)

It is generally the case that the bead positions $\bm{R}_{i}$ are not all
equal to each other. As a result, the spring forces on individual beads are
not zero, so the beads must be moving with respect to the solvent to create
the countervailing hydrodynamic forces.

The above are vector equations. They correspond to a total of $3N$ scalar
equations describing bead motions. The direction cosine for the $x$-component
of the force between beads $i$ and $i+1$ is
$(x_{i+1}-x_{i})/|\bm{R}_{i+1}-\bm{R}_{i}|$; corresponding forms give the $y$
and $z$ direction cosines for each force vector. The $N$ vector equations
given above therefore correspond to $3N$ scalar equations such as

$f\frac{dx_{i}}{dt}=-k(2x_{i}-x_{i-1}-x_{i+1})+{\cal F}_{xi},$ (7)

$x_{i}$ being the $x$-component of the bead coordinate of bead $i$ and ${\cal
F}_{xi}$ being the $x$ component of the thermal force on bead $i$. As
explained by Rouse in his original paper, the equations for the $x$
coordinates, for the $y$ coordinates, and for the $z$ coordinates are, except
for the coordinate label, the same as each other. Changing the $x$-component
of a particle’s position has no effect on the $y$ and $z$ components of the
forces on any particle, and correspondingly for displacements of a bead in the
$y$ or $z$ directions. The equations of motion for the $x$, $y$, and $z$
coordinates are thus completely uncoupled. The equations of motion therefore
partition into three sets of $N$ coupled equations, one set for each of the
three coordinate axes. Because each set of equations is the same as the
others, except for the label on the coordinates, only one set of $N$ equations
needs to be solved. The solutions for the other two sets of equations can be
obtained by a change of the coordinate label. The model is three-dimensional,
not one-dimensional; beads move in all three coordinate directions.

Equation 7 and the matching equations for beads $1$ and $N$ are a set of $N$
coupled linear differential equations whose coefficients are constants. The
solutions are a set of $N$ eigenmodes $Q_{n}$ describing motions parallel to
one of the three coordinate axes, each mode having a corresponding eigenvalue
$q_{n}$. One mode has eigenvalue $q_{0}=0$; that mode corresponds to uniform
translation of all beads in the $x$-direction. The other $N-1$ modes decay
exponentially ($\exp(-\Gamma_{n}t)$) in time; their relaxation rates
$\Gamma_{n}$ are

$\Gamma_{n}=\frac{8k\sin^{2}(n\pi/2N)}{f}$ (8)

with $n\in(1,N-1)$ being the mode label.

The normal mode amplitudes $C_{xn}(t)$ for the $x$-coordinate modes may be
calculated from the bead coordinates $x_{i}(t)$ via

$C_{xn}(t)=\frac{1}{N}\sum_{i=1}^{N}x_{i}(t)\cos\left(\frac{n\pi(i-1/2)}{N}\right).$
(9)

Entirely similar equations give the amplitudes $C_{yn}$ and $C_{zn}$ of the
$y$\- and $z$-coordinate modes. The inverse equations give the $x_{i}$ in
terms of the normal mode amplitudes as

$x_{i}(t)=C_{x0}(t)+2\sum_{n=1}^{N-1}C_{xn}(t)\cos\left(\frac{n\pi(i-1/2)}{N}\right).$
(10)

Totally similar equations describe the $y$ and $z$ modes. Standard
mathematical techniques show how the random forces ${\cal F}_{xi}(t)$ serve as
source terms, driving the fluctuations in the $C_{xn}(t)$.

There are three coordinate axes, so the relaxation rates $\Gamma_{n}$ are
three-fold degenerate. For each $n$, the same relaxation rate applies to all
three coordinate axes. The modes having degenerate eigenvalues are orthogonal;
their amplitudes fluctuate independently. The Rouse model thus has three
translational modes, each with eigenvalue zero, and $3N-3$ internal modes
(’internal’ in the sense that in each internal mode the beads move with
respect to each other) having non-zero eigenvalues.

On setting all but one of the $C_{xi}$ to zero, eq. 10 gives the
representation in bead position space of the eigenvector corresponding to
$C_{xi}$. The Rouse eigenvectors thus provide a set of $3N$ normal coordinates
that can replace the bead coordinates $\\{x_{i},y_{i},z_{i}\\}$ as a
specification of the polymer’s configuration. Eqs. 9 and 10 may be interpreted
as a pair of discrete Fourier transforms, in which $i-1/2$ plays the role of
the position coordinate, $n\pi/N$ is the wave vector, and $x_{i}$ and $C_{xn}$
are the amplitudes of the function and its transform at $i$ and $n$,
respectively.

Rouse uses the Rouse modes to describe the behavior of a polymer coil in a
shear flow. In Rouse’s calculation, a shear flow exerts forces on the polymer.
The polymer’s responses are described by the Rouse modes The Rouse solutions
therefore were taken by Rouse to be valid descriptions of polymer motion when
a shear flow is applied.

Polymer coils whose motions are described by Rouse’s model have one ill-
recognized property: They do not rotate. This property follows by comparison
with a standard problem in classical mechanics, namely the vibrational modes
of an isolated molecule. In general, an $N$-atom molecule has 3 translational
modes with eigenvalue zero, 3 rotational modes with eigenvalue zero, and
$3N-6$ internal vibrational modes. The internal modes are the modes that
change the distances between pairs of atoms. In translation and rotation the
distances between the atoms remain fixed. The Rouse problem only differs from
the molecular vibration problem in that the Rouse equations of motion are
overdamped, so the Rouse amplitudes relax exponentially at some rate
$\Gamma_{n}$ rather than oscillating at some frequency $\omega_{n}$. A polymer
coil is therefore like a vibrating isolated molecule in having a total of $3N$
modes. However, the $3N$ modes of the Rouse model include 3 translational
modes and $3N-3$ internal modes, for a total of $3N$ modes, leaving no modes
available for rotational motion.

The statement that Rouse chains do not rotate is not new. Rouse specifies in
his paper that a polymer coil under shear does not rotate, namely (his paper,
p. 1274, column 2) ”…since the velocity of the liquid has a nonvanishing
component only in the $x$ direction, the components $(\dot{y}_{j})_{\alpha}$
and $(\dot{z}_{j})_{\alpha}$ are zero.” $(\dot{y}_{j})_{\alpha}$ and
$(\dot{z}_{j})_{\alpha}$ are the velocities of bead $j$ in the $y$ and $z$
directions due to the shear. One or the other of them must be non-zero if the
chain is to be rotating. Rouse also argues his paper (p. 1274, column 2, top)
that ’…an atom at the junction between two submolecules…’ (springs) moves
’…with a velocity equal to that of the surrounding liquid…’ except for
Brownian motion, because, according to Rouse, otherwise there would be motion
of the solvent relative to the polymer chain, leading to energy dissipation.
If the beads only move with the liquid, then they can only be moving parallel
to the $x$-axis.

We now consider the Kirkwood-Riseman model. While both models refer to a line
of beads, the Kirkwood-Riseman model is radically different from the Rouse
model. The Kirkwood-Riseman model is based on three fundamental assumptions.
First, all distances between pairs of beads are treated as being their
statistico-mechanical average values; fluctuations and changes in these
distances are explicitly not included in the model. Second, the distribution
of beads around the chain center-of-mess is spherically symmetric. Third, the
system is massively overdamped, so that the inertia of the polymer coil is
negligible. These three assumptions completely define the chain dynamics, the
description of how a Kirkwood-Riseman polymer chain moves in solution.
Kirkwood and Riseman recognized that a polymer coil has internal modes
(”fluctuations”) so that polymer beads actually do move with respect to each
other, but these bead motions were specified as being not included in their
model.

The system is heavily overdamped, so its inertia is negligible. The total
force on the chain must therefore be zero. The moments of inertia of the chain
are negligible. The total torque on the chain must therefore also be zero. The
chain satisfies these two zero conditions by adjusting its linear velocity
$\bm{V}$ and its angular velocity $\bm{\Omega}$ until the total force and the
total torque on the chain both vanish.

Kirkwood and Riseman consider how a polymer coil moves in a shear field in
which the fluid velocity is

$\bm{u}_{i}=\bm{u}_{i}^{(0)}+Gy_{i}\bm{\hat{i}},$ (11)

Here $\bm{u}_{i}^{(0)}$ is a possible uniform motion of the fluid, $G$ is a
constant linear shear gradient, $y_{i}$ is the $y$-component of the vector
location of bead $i$, and $\bm{\hat{i}}$ is the unit vector parallel to the
$x$-axis. We do not follow Kirkwood’s notation closely. Kirkwood and Riseman
included in their model bead-bead hydrodynamic interactions as described by
the Oseen tensor. We return to these internal interactions below.

The Kirkwood-Riseman model describes a chain of $N$ beads whose coordinates
are $(\bm{R}_{1},\bm{R}_{2},\ldots\bm{R}_{N})$. Sequential bead positions form
a highly restricted random walk. Each bead is subject to a hydrodynamic force
$\bm{F}_{iH}$ exerted by the fluid. $\bm{F}_{iH}$ is determined by the bead
drag coefficient $f$, the velocity $\bm{v}_{i}$ of the bead, and the velocity
$\bm{u}_{i}$ that the fluid would have had, at the location of bead $i$, if
bead $i$ were not there, via

$\bm{F}_{iH}=f(\bm{u}_{i}-\bm{v}_{i}).$ (12)

$\bm{F}_{iH}$ is the hydrodynamic force on the bead, not the total force. The
total force on each bead, including the forces due to links to adjoining
beads, vanishes, so $\bm{F}_{iH}$ in general is non-zero. Correspondingly, the
bead and solvent velocities are in general not equal to each other.

In the Kirkwood-Riseman model, the velocity of bead $i$ is

$\bm{v}_{i}=\bm{V}+\bm{\Omega}\times\bm{s}_{i}.$ (13)

Here $\bm{s}_{i}$ is the vector from the chain center-of-mass to bead $i$, $V$
is a linear velocity, the same for each bead, and $\Omega$ is an angular
rotation rate, the same for each bead. Internal modes neglected in the model
would add to the right hand side of this equation an additional term
$\dot{\xi}_{i}$, the contribution of the internal modes to the bead velocity;
that term does not appear in the model.

Kirkwood and Riseman use the zero-total-force and zero-total-torque conditions
to determine $\bm{V}$ and $\bm{\Omega}$ in terms of $\bm{u}_{i}^{(0)}$ and
$G$, finding

$\bm{V}=\bm{u}_{i}^{(0)}+GY_{0}\bm{\hat{i}}$ (14)

and

$\bm{\Omega}=-\frac{G}{2}\bm{\hat{k}}.$ (15)

$Y_{0}$ is the y-coordinate of the polymer chain’s center of mass, and
$\bm{\hat{k}}$ is the unit vector in the $z$-direction. The model predicts
viscous dissipation because the bead velocity $\bm{v}_{i}$ and the solvent
velocity $\bm{u}_{i}$ cannot be equal at every point. For example, for most
beads $\bm{v}_{i}$ but not $\bm{u}_{i}$ will have a non-zero $y$-component

The Kirkwood-Riseman dynamic model is completely specified by eqs. 13-15. In
calculating the viscous dissipation, Kirkwood and Riseman include hydrodynamic
interactions between the polymer beads. Hydrodynamic interactions do affect
the polymer internal modes, which are not a part of the Kirkwood-Riseman
model, and do change the viscous dissipation created by the polymer coil.
However. these hydrodynamic interactions have no effect on the dynamic model
specified by eqs. 13-15. The dynamic model is not affected by intrachain
hydrodynamic interactions because bead-bead hydrodynamic interactions are
_internal_ forces, forces between different beads on the same chain. The total
force and the total torque exerted on a polymer chain by internal forces must
both vanish, an outcome guaranteed by Newton’s Third Law of Motion. Adding
hydrodynamic interactions has no effect on the motions described by eqs.
13-15.

In applying the Oseen tensor to describe bead-bead hydrodynamic interactions,
Kirkwood and Riseman took the distance between each pair of beads to be the
equilibrium average distance between those two beads. Fluctuations in those
interbead distances, and the time dependences of those fluctuations, were not
taken into account; i.e., Kirkwood and Riseman suppressed the time dependences
of the internal coordinates.

Considerably more could be said about the Kirkwood-Riseman and Rouse models
for polymer dynamics, but the above captures their descriptions of polymer
motion. The lack of complete consonance between the two models is apparent. In
the Rouse model, polymer beads partake of the polymer’s center-of-mass motion
and of internal modes that change the relative positions of the beads.
Contrariwise, in the Kirkwood-Riseman model, beads partake of center-of-mass
displacement and whole-body rotation, neither of which change any of the
distances between pairs of beads.

##  III Computational

A Brownian dynamics simulation for a polymer in a shear field was implemented.
The equations of motion resemble the Rouse equations of motion, but a shear
field has been added. The simulations reveal how a Rouse-like polymer coil
moves in a rheological experiment, as envisioned by Rouse and also by Kirkwood
and Riseman.

For beads $i\in(2,\ldots,N-1)$, we write these as

$\frac{d\bm{R}_{i}}{dt}=f^{-1}(-k(2\bm{R}_{i}-\bm{R}_{i-1}-\bm{R}_{i+1}))+Gy_{i}\mathbf{\hat{i}}+\mathbf{\cal{F}}_{i}(t),$
(16)

$\bm{\hat{i}}$ being the unit vector in the $x$-direction. For beads $1$ and
$N$ the equations of motion are

$\frac{d\bm{R}_{1}}{dt}=f^{-1}(-k(\bm{R}_{1}-\bm{R}_{2}))+Gy_{1}\bm{\hat{i}}+\mathbf{\cal{F}}_{1}(t),$
(17)

and

$\frac{d\bm{R}_{N}}{dt}=f^{-1}(-k(\bm{N}_{i}-\bm{R}_{N-1}))+Gy_{N}\bm{\hat{i}}+\mathbf{\cal{F}}_{N}(t).$
(18)

These equations describe chain motion relative to the chain center. The shear
force is directed in the $\bm{\hat{i}}$ direction, and changes linearly with
the distance in the $y$-direction from the chain center-of-mass. The thermal
forces $\mathbf{\cal{F}}_{i}(t)$ were generated using standard methods as
independent random variables having Gaussian distributions. The bead
displacements during a single time step $\Delta t$ are $\Delta
t\frac{d\bm{R}_{i}}{dt}$. These are re-evaluated after each time step to
compute the trajectory of each bead.

Multiple characteristic functions of chain behavior were determined. Most of
these functions were used as diagnostics to validate the core software. The
radius of gyration and mean-square radius of gyration were calculated. The
mean-square center-of-mass displacement was found to be linear in time, as
expected. The second $\langle x_{\alpha}x_{\beta}\rangle$ and fourth $\langle
x_{\alpha}^{2}x_{\beta}^{2}\rangle$ moments of the bead distribution around
the center of mass were calculated. Here $\alpha,\beta\in(1,3)$ designate
individual cartesian components of the vectors from the chain center to each
bead, the average being over all beads and all times.

Distribution functions for the nearest-neighbor distance, the second-nearest-
neighbor distance, the magnitude of the end-to-end vector, the distance from
the polymer center-of-mass to each bead, and the distances between all pairs
of beads were measured. The time autocorrelation functions
$\langle\bm{R}_{e}(0)\cdot\bm{R}_{e}(t)\rangle$ and
$\langle\bm{\hat{R}}_{e}(0)\cdot\bm{\hat{R}}_{e}(t)\rangle$ of the chain end-
to-end vector $\bm{R}_{e}=\bm{R}_{N}-\bm{R}_{1}$ and its unit vector
$\bm{\hat{R}}_{e}=(\bm{R}_{N}-\bm{R}_{1})/\mid\bm{R}_{N}-\bm{R}_{1})\mid$ were
obtained.

Using eq. 9, we calculated the time-dependent Rouse amplitudes
$C_{\alpha,n}(t)$ of the bead positions. For an $N$-bead system there are
$3(N-1)$ such components, plus the three $C_{\alpha,0}$ describing the polymer
center of mass position. We also calculated several time-dependent spatial
Fourier components

$a_{k,\alpha}(t)=\sum_{i=1}^{N}\cos(kr_{i\alpha}(t))$ (19)

of the bead locations. Here $k$ is the wavenumber for the transformation and
$r_{i\alpha}(t)$ is the $\alpha^{\rm th}$ Cartesian component of the location
of bead $i$, relative to the chain center of mass, at time $t$.

Finally, we calculated time-dependent Haar-likehaar1910a
waveletdaubechies1992a  components $c(n,\alpha,j)(t)$ and $d(n,\alpha,j)(t)$
of the particle positions. In this calculation, $n$ is the wavelet
decomposition level, $\alpha$ is again the Cartesian coordinate, and $j$
labels the wavelet location along the polymer chain. The maximum value of $j$
depends on $N$. For a $2^{m}$ bead polymer the upper limit on $j$ is $2^{m-n}$
with $m-n\geq 0$. The decomposition proceeds naturally if for some integer $m$
there are $2^{m}$ beads in the chain. The wavelet components are defined

$\displaystyle c(1,\alpha,j)$ $\displaystyle=$
$\displaystyle(r_{\alpha,2*j}+r_{\alpha,2*j-1})/2$ $\displaystyle
d(1,\alpha,j)$ $\displaystyle=$
$\displaystyle(r_{\alpha,2*j}-r_{\alpha,2*j-1})/2$ (20)

for $n=1$ and

$\displaystyle c(n,\alpha,j)$ $\displaystyle=$
$\displaystyle(c(n-1,\alpha,2*j)+c(n-1,\alpha,2*j-1)/2$ (21) $\displaystyle
d(n,\alpha,j)$ $\displaystyle=$
$\displaystyle(c(n-1,\alpha,2*j)-c(n-1,\alpha,2*j-1))/2$ (22)

for $n>1$. The $d(n,\alpha,j)$ differ from the spatial Fourier components and
the Rouse components in that they are localized; they refer to the behavior of
specific parts of the polymer coil. In contrast, the $a_{k,\alpha}(t)$ and the
$C_{\alpha,n}(t)$ are both global variables, each depending on the relative
positions of all the beads in the chain.

For the $a_{\alpha,k}(t)$, $C_{\alpha,n}(t)$, and $d(n,\alpha,j)(t)$ the
temporal self correlation functions were evaluated. For the $a_{\alpha,k}(t)$
and $C_{\alpha,n}(t)$, we also calculated the temporal cross-correlation
functions, e.g., $\langle C_{\alpha,n}(t)C_{\beta,m}(t)\rangle$ for
$\alpha\neq\beta$ and/or $m\neq n$. There are $3(N-1)$ Rouse internal modes
and therefore $9(N-1)^{2}$ Rouse-Rouse self and cross-correlation functions.
In Rouse’s original model, if either $\alpha\neq\beta$ or $m\neq n$ or both,
the temporal crosscorrelation function vanishes.

How does one show that an object is performing whole-body rotation? For a
fluid velocity in the $x$ direction, with a non-zero velocity shear gradient
$dv_{x}/dy$, the induced angular velocity should on the average be parallel to
the $z$-axis. A simple test is advanced. If the beads are each taken to be
performing circular motion, the instantaneous angular rotation can be written

$\sum_{i=1}^{N}\bm{R}_{i}\times\bm{v}_{i}=\sum_{i=1}^{N}\bm{R}_{i}\times(\bm{\omega}\times\bm{R}_{i})$
(23)

The $z$ component of $\bm{L}$ is $\bm{\hat{k}}\cdot\bm{L}$. Applying standard
identities one obtains for the rotational velocity
$\bm{\omega}=\omega_{z}\bm{\hat{k}}$ around the $z$-axis

$\omega_{z}\sum_{i=1}^{N}(\langle(x_{i})^{2}\rangle+\langle(y_{i})^{2}\rangle)=\left\langle\sum_{i=1}^{N}x_{i}\frac{dy_{i}}{dt}-\sum_{i=1}^{N}y_{i}\frac{dx_{i}}{dt}\right\rangle.$
(24)

Corresponding forms describe rotation around the $x$ and $y$ axes. The
velocities are related to the bead displacements during a single time step
$\Delta t$, namely bead $i$’s displacements are $\Delta x_{i}=\Delta t\
dx_{i}/dt$, $\Delta y_{i}=\Delta t\ dy_{i}/dt$, and $\Delta z_{i}=\Delta t\
dz_{i}/dt$, so in evaluating the right-hand-side of eq. 24 we replace the
velocities with the single-step displacements.

For whole-body rotation the two terms on the right-hand-side of eq. 24 are
equal by symmetry. In the absence of rotation, the first sum on the right hand
side of the equation will average to zero. A polymer chain is not a solid
object that performs rigid-body motion, so $\bm{\omega}$ should not be
overinterpreted. Sablic, et al.sablic2017a  discuss rotation in terms of
Eckart frames, and note alternative definitions of rotation rates and their
physical implications.

Simulations were made for 8 and 16 bead chains at shear rates $G\in(0,0.15)$;
we treat here outcomes from 16 bead chains. In the simulations, we chose
$k=1$, $f=1$, nominal temperature $k_{B}T=1$, a basic time step $\Delta
t=0.001$, a unit diffusion step to be $\Delta r=(2k_{B}T\Delta T/f)^{1/2}$,
with a unit force $kr_{i,i+1}$ giving a displacement $\Delta t/f$. The
characteristic functions were computed every ten time steps. A simulation with
$\Delta t=0.0003$ gave very nearly the same results as a simulation using the
longer time step.

Calculations were performed on an 8-core 3.4 Ghz CPU and an Nvidia Tesla K-40
GPU using locally written software run under Simply Fortran 2 and PGI Fortran.
Our code and necessary files are found on ResearchGate.comresearchgate . In
production runs, polymer positions were advanced through $1\cdot 10^{8}$ time
steps. Prior to each production run, a $1\cdot 10^{7}$ or longer timestep
thermalization run was performed.

##  IV Results

We first consider the effect of shear on the polymer coil’s shape. As shown by
Figure 1, our results include both a small-shear region, in which the polymer
coil is not distorted significantly, and a large-shear region, in which the
polymer coil on the average is distorted by the shear. Figure 1 plots the
second moments $\langle x_{i}^{2}\rangle$, $\langle y_{i}^{2}\rangle$, and
$\langle x_{i}y_{i}\rangle$ against shear rate. At zero shear, $\langle
x_{i}^{2}\rangle=\langle y_{i}^{2}\rangle$. With increasing shear, the polymer
is stretched in the $x$ direction, but not in the $y$ or (not shown) $z$
directions, so that $\langle x_{i}^{2}\rangle>\langle y_{i}^{2}\rangle$. The
shear field creates a non-zero $\langle x_{i}y_{i}\rangle$ correlation that
increases nearly linearly with shear rate $G$. However, a shear $dv_{x}/dy$
has no effect on bead displacement in the $z$ direction, so $\langle
y_{i}z_{i}\rangle$ and $\langle x_{i}z_{i}\rangle$ remain equal to zero
regardless of the shear rate.

Figure 1: Effect of the shear rate $G$ on the polymer coil shape, from the
equal-time correlation functions $\langle(x_{i})^{2}\rangle$ (open circles),
$\langle(y_{i})^{2}\rangle$ (crosses), and $\langle x_{i}y_{i}\rangle$ (filled
circles).

We now examine the most fundamental question. Do simulated chains rotate when
placed in a shear field? Eq. 24 supplies the test. When the shear rate is
greater than zero, the right-hand-side of eq. 24 is non-zero. The polymer
chair rotates around the $z$-axis. We also evaluated the analogs of eq. 24 for
rotation around the $x$ and $y$ axes. Our shear field creates no rotation
around the $x$ or $y$ axes, to within the accuracy of the simulation. For
$G>0$, rotation in the $x-y$ plane should be clockwise, i.e., $\omega<0$, as
is found.

Is the motion actually rotational? For circular motion, the two terms on the
right-hand-side of eq. 24 should average to the same value. Figure 2 shows
that they do. We find

$\left\langle\sum_{i=1}^{N}x_{i}\frac{dy_{i}}{dt}\right\rangle=-\left\langle\sum_{i=1}^{N}y_{i}\frac{dx_{i}}{dt}\right\rangle.$
(25)

Rouse predicts that the left hand side of eq. 25 vanishes; it does not vanish.

Figure 2: The $\sum_{i=1}^{N}x_{i}\ dy_{i}/dt$ ($\circ$) and
$\sum_{i=1}^{N}y_{i}\ dx_{i}/dt$ ($\times$) contributions to $L$ as determined
at various shear rates $G$. As shown in eq. 25, for rotational motion the two
terms should be equal, as seen here to good approximation.

The nominal angular motion $L$

$L=\left\langle\sum_{i=1}^{N}x_{i}\frac{dy_{i}}{dt}\right\rangle-\left\langle\sum_{i=1}^{N}y_{i}\frac{dx_{i}}{dt}\right\rangle.$
(26)

Figure 3 gives $L$ and the rotation rate $\omega_{z}$ (eq. 24), as functions
of the applied shear $G$. $L$ is precisely linear in $G$ up to the largest $G$
that we examined, as predicted by Kirkwood and Riseman. Because the polymer
coil distorts when a shear is applied, $L$ and $\omega_{z}$ are not simply
linearly proportional to each other.

Figure 3: Nominal angular motion $|L|$ (eq. 26) and angular rotation rate
$|\omega_{z}|$ (eq. 24) for the 16-bead chain, as induced by the applied
shear.

We now go beyond Kirkwood and Riseman, and beyond Rouse. Kirkwood and Riseman
partitioned chain motions in shear into uniform translation, whole-body
rotation, and residual contributions of internal modes. They ignored internal
modes. Rouse implicitly assumed that mode relaxations are not perturbed by an
applied shear. Here we ask whether an applied shear actually affects the
internal modes, as represented by the Rouse amplitudes $C_{n\alpha}(t)$ and
their relaxation rates $\Gamma_{n\alpha}$. Rouse’s solutions indicate

$\langle
C_{n\alpha}(0)C_{m\beta}(t)\rangle=\delta_{mn}\delta_{\alpha\beta}\langle(C_{n\alpha}(0))^{2}\rangle\exp(-\Gamma_{n\alpha}t)$
(27)

According to Rouse’s analysis: The fluctuating amplitudes $C_{n\alpha}(t)$ are
uncorrelated. Modes with different $n$ fluctuate independently of each other.
Modes with the same $n$, but corresponding to different directions (different
$\alpha$), also fluctuate independently. The temporal correlation function for
each mode decays exponentially in time.

We first consider the autocorrelation functions $\langle
C_{n\alpha}(0)C_{n\alpha}(t)\rangle$. We obtained the $\Gamma_{n\alpha}$ and
$\langle(C_{nx}(0))^{2}\rangle$ as functions of the shear rate by fitting an
early-time segment of each $\langle C_{nx}(0)C_{n\alpha}(t)\rangle$ to a
single exponential. Figure 4a shows the decay rates $\Gamma_{nx}$ as functions
of the shear rate. Open circles mark the $n=1$ mode. Figure 4b shows the
corresponding mean-square average amplitudes $\langle(C_{nx}(0))^{2}\rangle$.

For some modes, the decay rates and initial amplitudes are significantly
shear-sensitive. For $n=3$ and 2, and much more markedly for $n=1$, the mode
relaxation rates $\Gamma_{nx}$ decrease with increasing shear rate, while the
corresponding mode amplitudes $\langle(C_{nx}(0))^{2}\rangle$ increase with
increasing shear rate. For $n>3$, $\Gamma_{nx}$ and
$\langle(C_{nx}(0))^{2}\rangle$ are very nearly independent of shear rate. The
relaxation rates and amplitudes for the $y$ and $z$ components of the Rouse
modes are independent of the shear rate. We did not explore the dependence of
this result on chain length. These non-trivial dependences of the mode
amplitudes and relaxation rates on shear rate are contrary to Rouse’s picture,
in which the $\Gamma_{nx}$ and $\langle(C_{nx}(0))^{2}\rangle$ are not
affected by solvent shear.

Figure 4: Effect of shear rate $G$ on the Rouse-Rouse time correlation
functions $\langle C_{nx}(0)C_{nx}(t)\rangle$, showing the dependences of (a)
$\Gamma_{nx}$ and (b)$\langle(C_{nx}(0))^{2}\rangle$ on $G$. The open circles
denote the $n=1$ mode, the $n=2$ to $n=7$ modes moving seriatim away from the
$n=1$ mode’s behavior. $\Gamma_{1x}$ depends on $G$ down to the smallest non-
zero $G$ that we studied.

When shear is applied, some Rouse modes become cross-correlated. Figure 5
shows the $xy$ cross-correlations $\langle C_{nx}(0)C_{ny}(t)\rangle$. These
cross-correlation functions vanish in the Rouse model. They are not zero in
our simulations. The corresponding $yz$ and $zx$ crosscorrelation functions
(not shown) do vanish no matter whether or not shear is applied, as do all
crosscorrelation functions $\langle C_{n\alpha}(0)C_{m\beta}(t)\rangle$ with
$n\neq m$. The cross-correlation functions are not exponentials; they first
increase and then fall off rapidly.

Figure 5: The Rouse-Rouse cross correlation functions $\langle
C_{nx}(0)C_{ny}(t)\rangle$ for a polymer coil in shear with (a) $G=0.005$ and
(b) $G=0.150$. Open circles mark $n=1$; crosses are (a) $n=6$ and (b) $n=5$.
Note the large change in the vertical scale between these two figures.

The time dependences of the $\langle C_{nx}(0)C_{ny}(t)\rangle$ are
qualitatively only little affected by the shear rate, but the initial
amplitudes $\langle C_{nx}(0)C_{ny}(0)\rangle$ depend strongly on $G$. Figure
6 shows this dependence. To reasonable approximation the initial amplitudes of
the cross-correlation functions are linear in the shear rate $G$.

Figure 6: The initial amplitudes $\langle C_{nx}(0)C_{ny}(0)\rangle$ of the
Rouse-Rouse cross correlation functions as functions of the shear $G$. Open
circles mark the $n=1$ modes. $\langle C_{nx}(0)C_{ny}(0)\rangle$ decreases
monotonically with increasing $n$. Straight lines are linear fits to $\langle
C_{nx}(0)C_{ny}(0)\rangle=aG+b$, $a$ and $b$ being fitting constants.

The x-y correlations are clearly driven by rotation. Rotational motion around
the $z$-axis will pump amplitude directly from $C_{nx}$ into the corresponding
$C_{ny}$, as may be seen by considering the rotation of a perfectly rigid
body. Whatever the amplitude $C_{nx}$ was at time 0, at the moment the body
has rotated through 90 degrees the component $C_{ny}$ is exactly equal in
magnitude to the initial component $C_{nx}$. If the $C_{nx}$ and $C_{ny}$ were
initially uncorrelated, rotation will cause the cross-correlation functions
$\langle C_{nx}(0)C_{ny}(t)\rangle$ to increase with increasing time. Indeed,
a close examination of the cross-correlation functions in Figure 5 suggests
the presence of such an increase at longer times.

Figure 7: Temporal autocorrelation functions for (a) several spatial fourier
components of the bead positions, and (b) a Haar-like wavelet decomposition of
the particle positions. These representative results refer to a 16-bead chain
and $G=0.03$. The wavelet components are $d(1,x,1)$ $(\square)$, $d(1,x,2)$
$(\circ)$, $d(1,x,3)$ $(\vartriangle)$, $d(1,x,4)$ $(\triangledown)$,
$d(2,x,1)$ $(+)$, $d(2,x,2)$ $(\times)$, and $d(3,x,1)$ $(\bullet)$.

Figure 7 shows representative measurements of two sets of collective
coordinates that could be used as alternatives to Rouse coordinates. It is not
claimed that either of these sets is necessarily the best possible choice for
a set of collective coordinates, but only that there are alternatives to
Rouse’s coordinates that may be worth examining. The polymer coil had 16
beads; the shear rate was $0.03$. The spatial Fourier components do not decay
as simple exponentials, in that they decay too slowly at longer times, but
there is no sign in them of multimodal behavior. The wavelet decompositions
provide measurements of true localized motions, showing that even in this
extremely simple model for polymer dynamics there is room for local
differentiation of structural relaxation. We show here only the components
corresponding to one-half of the full chain; the corresponding components for
the other half of the chain show exactly the same set of behaviors. In
considering the series $d(1,x,j)$ (open symbols) for $j\in(1,4)$, $d(1,x,1)$,
which relaxes the most rapidly, corresponds to the motions of the outer pair
of beads. The $d(1,x,j)$ for $j>2$ have clearly bimodal relaxations, speaking
to more complex chain dynamics nearer to the center of the polymer. The
$d(2,x,j)$ emphasize the differences between inner and outer beads of the
polymer. $d(2,x,1)$,which corresponds to the outer four beads of the polymer,
has a non-exponential but unimodal relaxation; $d(2,x,2)$, describing the four
beads nearest the chain center, shows a visibly bimodal relaxation.

##  V Discussion

This paper describes a simulational study of the motions of a polymer in a
shear field. Comparison was made with the Kirkwood-Riseman and Rouse
treatments of the dynamics of an isolated polymer chain. We show that the
Kirkwood-Riseman model of polymer dynamics, in which a polymer coil translates
and rotates when subject to the influence of a shear, is qualitatively
correct. The Rouse model, in which polymer coils do not rotate during
viscometric studies, is incorrect as applied to the viscosity increment of a
polymer in solution. We note several alternatives to Rouse coordinates that
could in principle serve as descriptions of polymer internal motions.

It is certainly legitimate to ask how the issues raised here were not already
noticed. It is not suggested here that there were past errors. Several
contributory factors are readily identified. First, while there are multiple
excellent presentations of the Rouse-Zimm model, e.g., refs. 10 and 11,
equivalent presentations of the Kirkwood-Riseman model more recent than their
original paper are far less common, so there is little familiarity with the
Kirkwood-Riseman model. Second, in the absence of shear, the two models
converge; computer simulations of polymer coils in unsheared liquids cannot
readily detect the disagreement between the models. Third, in order to
identify our issues, one would have needed to analyse a chain trajectory with
the correct diagnostic, e.g., eq. 25, but in the context of the Rouse-Zimm
model there is no rational reason to develop such a diagnostic. As a result,
in aopst studies many fine questions have been asked about the nature of
polymer dynamics, but not the questions answered here.

Larson and co-workersjain2008a ; dalal2013a  provide considerable evidence
that potential energies more precise than Rouse’s potential can cause a
chain’s dynamics to deviate from simple Rouse behavior. Jain and
Larsonjain2008a  made Brownian dynamics simulations of a string of polymer
beads to which stiff springs, bond-angle, and bond-torsion-angle forces were
added seriatim. They calculated the time autocorrelation functions for the
polymer end-to-end vector and the connector unit-vector autocorrelation
functions, the latter being averaged over all springs in the chain. Dalal and
Larsondalal2013a  extended these results, showing that adding side groups,
chain excluded-volume effects, and explicit treatment of solvent molecules
jointly lead to the experimentally-observed single-exponential relaxation for
short chains. They also note what they viewed as an interesting coincidence,
namely that the relaxation times for the orientation of the chain end-to-end
vector and the single-spring orientation vectors are very nearly the same.

The difficulty with the Rouse model is apparent on comparing eqs. 4-6 with
eqs. 16-18. The second set of equations include a solvent shear force
$Gy_{i}\bm{\hat{i}}$ on each bead. That force is absent from the first set of
equations, the equations solved by Rouse. Rouse calculated how a bead-spring
polymer coil would evolve in time in a quiescent fluid. In a quiescent fluid,
the polymer coil by symmetry has no tendency to rotate. When a fluid shear
field $Gy_{i}\bm{\hat{i}}$ is included in the calculation, the forces on the
beads are changed. The motions of the beads therefore also change, and are no
longer the motions described by Rouse. The Rouse model thus does not describe
polymer dynamics during a rheological experiment.However, in his original
paper, Rouse uses his quiescent-fluid solutions to calculate dissipation and
hence viscosity increment for a polymer in a shear flow, as though his
solutions were applicable under these conditions.

The above has focused on a polymer coil in an imposed macroscopic shear field,
as encountered in viscoelastic measurements. However, the fluctuation-
dissipation theorem gives us two other circumstances in which polymer coils
find themselves in shear fields:

First, consider a polymer coil performing Brownian motion. The fluctuation-
dissipation theorem indicates that if the chain diffuses through some distance
in a given time, the chain motions will have correlated fluid motions, exactly
as if the chain’s motion were being created by an imposed external force. That
fluid motion, the wake created in the solvent by the polymer, acts on other
polymer coils, causing them to move in turn. Because the fluid flow is not the
same everywhere, those other polymer coils are subject to a fluid shear field
which causes them to translate, rotate, and create fresh flow fields in the
surrounding solvent. This image of flow fields being scattered and re-
scattered by diffusing macromolecules forms the core of modern theoretical
treatments of the diffusion of interacting spherical colloidal
particlesphillies2016b , these theoretical treatments giving reasonably
accurate quantitative predictions for colloidal behavior. It should therefore
not be surprising that the same general approach is valid for interacting
polymer coils in solution. That is, via the fluctuation-dissipation theorem we
can extend the Kirkwood-Riseman model from treating a single isolated polymer
molecule to treat the hydrodynamic interactions between polymer molecules.
Indeed, there is a substantial development of polymer dynamics in non-dilute
solutions based on computing the hydrodynamic interactions between polymer
coilsphillies2016a ; phillies1988a ; phillies1993a ; phillies1997a ;
phillies1998a ; phillies2002c ; phillies2002b ; phillies2002a ; phillies2004a
.

Second, consider a polymer coil in a quiescent fluid. On the average, there is
no tendency for the molecule to rotate in any direction. However, the
fluctuating thermal forces on the polymer beads create evanescent fluctuating
torques on the molecule as a whole, causing the polymer end-to-end vector to
perform rotational diffusion, so that its later positions gradually become
decorrelated from its earlier positions. The end-to-end vector is a sum of the
individual bead-to-bead vectors, so there is a component of each bead-to-bead
vector that is correlated with the chain end-to-end vector. The bead-to-bead
vectors can only become completely uncorrelated on the time scale on which the
chain end-to-end vector relaxes. The result of whole-chain rotational
diffusion is that the spring unit-vector correlation functions will in part
relax on the time scale on which the chain end-to-end vector relaxes,
precisely as found by Dalal and Larsondalal2013a .

Rouse modes and the Rouse model are used in an extremely large number of
different contexts. I have not here generated a full list of contexts the
Rouse model is inappropriate, though clearly any theoretical problem in which
a polymer chain is placed in a shear field must be on that list. Nor have I
considered here any extensions to the Kirkwood-Riseman model.

## References

  * (1) Kirkwood, J. G.; Riseman, J. The Intrinsic Viscosities and Diffusion Coefficients of Flexible Molecules in Solution. _J. Chem. Phys._ 1948, 16, 565-573.
  * (2) Rouse, P. E. A Theory of the Linear Viscoelastic Properties of Dilute Solutions of Coiling Polymers. _J. Chem. Phys._ 1953, 21, 1272-1280.
  * (3) Zimm, B. H. Dynamics of Polymer Molecules in Dilute Solution: Viscoelasticity, Flow Birefringence, and Dielectric Loss. _J. Chem. Phys._ 1956, 24, 269-278.
  * (4) Tsalikis, D. G.; Koukoulas, T.; Mavrantzas, V. G.; Pasquino, R.; Vlassapoulos, D.; Pyckhout-Hintzen, W.; Wischnewski, A.; Monkenbusch, M.; Richter, D. Microscopic Structure, Conformation, and Dynamics of Ring and Linear Poly(ethylene oxide) from Detailed Atomistic Molecular Dynamics Simulations: Dependence on Chain Length and Direct Comparison with Experimental Data. _Macromolecules_ 2017, 50, 2565-2584.
  * (5) Perez-Aparicio, R.; Alvarez, F.; Arbe, A.; Willner, L.; Richter, D.; Falus, P.; Colmenero, J. Chain Dynamics of Unentangled Poly(ethylene-alt-propylene) Melts by Means of Neutron Scattering and Fully Atomistic Molecular Dynamics Simulations. _Macromolecules_ 2011, 44, 3129-3139.
  * (6) Kalathi, J. T.; Kumar, S. K.; Rubinstein, M.; Grest, G. S. Rouse Mode Analysis of Chain Relaxation in Homopolymer Melts. _Macromolecules_ 2014, 47, 6925-6931.
  * (7) Doi, M.; Edwards, S. F. _The Theory of Polymer Dynamics._ Oxford University Press: Oxford (1986).
  * (8) Padding, J. T. Theory of Polymer Dynamics. Advanced Courses in Macroscopic Physical Chemistry. (Han-sur-lesse winterschool 2005) https://www.scribd.com/document/181367522 /THEORY-OF-POLYMER-DYNAMICS-Paddings.
  * (9) Haar, A. Zur Theorie der orthogonalen Funktionensysteme. _Mathematische Annalen_ 1910, 69, 331–371.
  * (10) Daubechies, I. _Ten Lectures on Wavelets._ CBMS-NSF Regional Conf. Ser. Appl. Math., Soc. Indust. Appl. Math., Philadelphia (1992).
  * (11) Sablic, J.; Delgado-Buscalioi, R.; Praprotnick, M. Application of the Eckart Frame to Soft Matter: Rotation of Star Polymers under Shear Flow. Cornell University Library arxiv.org ArXiv:1707.09170v1[cond-mat.soft].
  * (12) www.researchgate.net/publication/320988644
  * (13) Jain, S.; Larson, R. G. Effects of Bending and Torsional Potentials on High-Frequency Viscoelasticity of Dilute Polymer Solutions. _Macromolecules_ 2008, 41, 3692-3700.
  * (14) Dalal, I. S.; Larson, R. G. Explaining the Absence of High-Frequency Viscoelastic Relaxation modes of Polymers in Dilute Solution. _Macromolecules_ 2013, 46, 1981-1992.
  * (15) Phillies, G. D. J. Diffusion in Crowded Solutions. _Advances in Chemical Physics_ 2016, 48, 277-358.
  * (16) Phillies, G. D. J. The Hydrodynamic Scaling Model for the Dynamics of Non-Dilute Polymer Solutions: A Comprehensive Review. Cornell University Library http://arxiv.org _arxiv.org/abs/1606.09302_ 2016\.
  * (17) Phillies, G. D. J. Quantitative Prediction of $\alpha$ in the Scaling Law for Self-Diffusion. _Macromolecules_ 1988, 21, 3101-3106.
  * (18) Phillies, G. D. J.; Kirkitelos, P. C. Higher-Order Hydrodynamic Interactions in the Calculation of Polymer Transport Properties. _J. Polymer Sci. B: Polymer Physics_ 1993, 31, 1785-1797.
  * (19) Phillies, G. D. J.; Lacroix, M.; Yambert, J. Probe Diffusion in Sodium Polystyrene Sulfonate - Water: Experimental Determination of Sphere-Chain Binary Hydrodynamic Interactions. _J. Phys. Chem. B_ 1997, 101, 5124-5130.
  * (20) Phillies, G. D. J. Derivation of the Universal Scaling Equation of the Hydrodynamic Scaling Model via Renormalization Group Analysis. _Macromolecules_ 1998, 31, 2317-2327.
  * (21) Phillies, G. D. J. Low-Shear Viscosity of Non-Dilute Polymer Solutions from a Generalized Kirkwood-Riseman Model. _J. Chem. Phys._ 2002, 116, 5857-5866.
  * (22) Phillies, G. D. J. Viscosity of Hard Sphere Suspensions. _J. Colloid Interface Sci._ 2002, 248, 528-529.
  * (23) Phillies, G. D. J. Self-Consistency of Hydrodynamic Models for the Low-Shear Viscosity and the Self-Diffusion Coefficient. _Macromolecules_ 2002, 35, 7414-7418.
  * (24) Merriam, S. C.; Phillies, G. D. J. Fourth-Order Hydrodynamic Contribution to the Polymer Self-Diffusion Coefficient. _J. Polymer Sci. B Polymer Physics_ 2004, 42, 1663-1670.


# Electronic correlations and superconducting instability in La3Ni2O7 under
high pressure

Frank Lechermann  Institut für Theoretische Physik III, Ruhr-Universität
Bochum, D-44780 Bochum, Germany  Jannik Gondolf  Institut für Theoretische
Physik III, Ruhr-Universität Bochum, D-44780 Bochum, Germany  Steffen Bötzel
Institut für Theoretische Physik III, Ruhr-Universität Bochum, D-44780 Bochum,
Germany  Ilya M. Eremin  Institut für Theoretische Physik III, Ruhr-
Universität Bochum, D-44780 Bochum, Germany

###### Abstract

Motivated by the report of superconductivity in bilayer La3Ni2O7 at high
pressure, we examine the interacting electrons in this system. First-
principles many-body theory is utilized to study the normal-state electronic
properties. Below 100 K, a multi-orbital non-Fermi liquid state resulting from
loss of Ni-ligand coherence within a flat-band dominated low-energy landscape
is uncovered. The incoherent low-temperature Fermi surface displays strong
mixing between Ni-$d_{z^{2}}$ and Ni-$d_{x^{2}-y^{2}}$ orbital character. In a
model-Hamiltonian picture, spin fluctuations originating mostly from the
Ni-$d_{z^{2}}$ orbital give rise to strong tendencies towards a
superconducting instability with $d_{x^{2}-y^{2}}$ order parameter. The
dramatic enhancement of $T_{\rm c}$ in pressurized La3Ni2O7 is due to stronger
Ni-$d_{z^{2}}$ correlations compared to those in the infinite-layer
nickelates.

Introduction.— In a recent finding, Sun et al. [1] reported superconductivity
near a temperature $T=80$ K in bulk single crystalline La3Ni2O7 at pressures
$p>14$ GPa. This adds a whole new chapter to the young research field of
superconducting (SC) nickelates, hosting high-$T_{\rm c}$ cuprate-akin NiO2
square-lattice planes. The field has been inaugurated by the discovery of
electron pairing in thin films of Sr-doped infinite-layer NdNiO2 with a
$T_{\rm c}\sim 15$ K in 2019 [2]. Early follow-up studies [2, 3, 4, 5, 6, 7,
8] detected similar SC phases in thin films of Pr1-xSrxNiO2, La1-xSrxNiO2 and
also in multilayer Nd6Ni5O12 thin films. While these reduced SC nickelates
share a common motif by the lack of apical oxygens (resulting from a
topotactic reaction) and Ni$(3d^{9-\delta})$ oxidation states, the
characteristics of bilayer La3Ni2O7 differ. It still holds the apical oxygens
and Ni formally has $3d^{7.5}$ configuration. Furthermore, a comparison to
high-$T_{\rm c}$ cuprates with their key Cu-$d_{x^{2}-y^{2}}$ single-orbital
character becomes quite stretched. Whereas there is an ongoing debate
concerning a dominant single Ni-$d_{x^{2}-y^{2}}$ [9, 10, 11, 12, 13, 14, 15,
16, 17, 18, 19, 20] versus dominant multi-orbital Ni-$e_{g}$ [21, 22, 23, 24,
25, 26, 27, 28] low-energy physics in reduced SC nickelates, the nominal hole
doping away from Ni$(3d^{9})$ is that large in La3Ni2O7 as to render Ni multi-
orbital physics inevitable. A Ni-$e_{g}$ multi-orbital picture for infinite-
layer nickelates results in a competition between SC instabilitites of varying
flavor [29].

On a wider scope, two further issues appear relevant. First, even if the
formal oxidation state reads accordingly, a Ni$(3d)$ occupation well below
$n_{d}=8$ is hardly occuring in known nickel oxides. Instead in most cases, a
$3d^{8}L$ state incorporating holes on ligand oxygen is realized [30, 31, 32,
33, 34, 35, 26], also accompanied by a smaller charge-transfer energy
$\Delta=\varepsilon_{d}-\varepsilon_{p}$ between Ni$(3d)$ and O$(2p)$.
Secondly, bilayer oxides from the $p$-layered Ruddlesden-Popper series
Ap+1TMpO3p+1 (A: rare-earth, alkaline-earth; TM: transition; metal) often
display a much more delicate normal-state low-energy physics than
corresponding single-laver systems. This is e.g. exemplified for ruthenates
[36, 37, 38] and iridates [39]. Previous theoretical accounts of bilayer
lanthanum nickelate focussed on the paramagnetic metal [40, 41, 42] at ambient
pressure. From density-functional theory (DFT), a charge-density wave state
was predicted [43], while DFT+Hubbard $U$ considerations [44] remark the
possible relevance of magnetically-ordered states.

In this work, we provide a theoretical description of the correlated
electronic structure of paramagnetic La3Ni2O7 under high pressure, by
employing the combination of DFT, self-interaction correction (SIC) and
dynamical mean-field theory (DMFT), i.e. the so-called DFT+sicDMFT approach
[45]. Moreover, a model Hamiltonian perspective via the random-phase
approximation (RPA) onto the possible superconducting instabilitites from spin
fluctuations is presented. We reveal an intriguing low-energy physics of
pressurized La3Ni2O7 in the normal state. The highly-correlated interplay
between Ni-$d_{z^{2}}$, Ni-$d_{x^{2}-y^{2}}$ and O$(2p)$, with the former
displaying partial flat-band character, gives rise to a distinct non-Fermi-
liquid (NFL) regime below $T<100$ K. The model-RPA investigation points to a
$d_{x^{2}-y^{2}}$ SC order parameter that would emerge from this multi-orbital
scenario. Most importantly, we argue that the Ni-$d_{z^{2}}$ orbital is much
more correlated in the bilayer case than in infinite-layer nickelates.
Consequently, it interacts in a much more concerted fashion with
Ni-$d_{x^{2}-y^{2}}$, which could explain the much higher $T_{c}$ in the
bilayer case. Note that Ni-$d_{x^{2}-y^{2}}$ is always close to half filling
in superconducting nickelates and alone cannot explain the difference in
$T_{\rm c}$ between these two classes of materials.

Theoretical Approach.— The charge self-consistent [46] DFT+sicDMFT framework
[45], where the Ni sites act as quantum impurities and Coulomb interactions on
oxygen enter by SIC on the pseudopotential level [47], is put into practise.
The DFT part consists of a mixed-basis pseudopotential code [48, 49, 50] and
SIC is applied to the O$(2s,2p)$ orbitals via weight factors $w_{p}$. While
the $2s$ orbital is fully corrected with $w_{p}=1.0$, the choice [47, 45, 22,
23, 24, 25, 26] $w_{p}=0.8$ is used for $2p$ orbitals. Continuous-time quantum
Monte Carlo in hybridization expansion [51] as implemented in the TRIQS code
[52, 53] solves the DMFT problem. A five-orbital general Slater-Hamiltonian,
parameterized by Hubbard $U=10$ eV and Hund exchange $J_{\rm H}=1$ eV [22, 23,
24, 25, 26], governs the correlated subspace defined by Ni projected-local
orbitals [54]. Crystallographic data are taken from experiment [1]. Further
calculational details are given in the supplementary material (with Refs [55,
56, 57, 58, 59]).

Results.—

Figure 1: k-integrated electronic spectrum from DFT+sicDMFT at $T=80$ K. (a)
$Fmmm$ crystal structure of bilayer La3Ni2O7 with $c$-axis vertical: La
(green), Ni(grey) and O(red). (b) Total spectral function of low-pressure
$Amam$ and high-pressure $Fmmm$ phase (inset: low-energy blow up). (c) Site-
and orbital-projected spectral function for the $Fmmm$ phase (inset: low-
energy blow up).

In experiment, there is a structural transition in La3Ni2O7 from a low-
pressure $Amam$ phase with finite NiO6 octahedral tilting to a high-pressure
$Fmmm$ phase without tilting [1]. The $Amam$ (space group 63) crystal
structure has four equivalent Ni sites in the primitive cell, in contrast to
two Ni sites for the $Fmmm$ (space group 69) one. Key feature of the present
latter structure (see Fig. 1a) is a rather small distance of $\sim 1.76$ Å
between Ni and apical oxygen within the bilayer. The calculations show (cf.
Fig. 1b) that the spectrum of the $Amam$ phase is (nearly) gapped at low-
energy, in line with the measured different transport properties [1]. For the
rest of the paper, we will restrict the discussion to the properties of the
high-pressure $Fmmm$ phase.

The electronic spectrum in Fig. 1c exhibits a metallic state with
Ni-$e_{g}\,\\{d_{z^{2}},d_{x^{2}-y^{2}}\\}$ and O$(2p)$ character at the Fermi
level $\varepsilon_{\rm F}$, but lacks strong quasiparticle (QP) signature.
The Ni-$t_{2g}\,\\{d_{xz},d_{yz},d_{xy}\\}$ manifold is mostly filled with a
peak at $\sim-1.5$ eV. The main O$(2p)$ weight peaks at $\sim-3.5$ eV. Sizable
O$(2p)$ weight in the unoccupied part of the spectrum points to ligand-hole
states. And indeed, the integrated projected spectral parts yield occupations
$n_{d}=7.98$ and $n_{p}=5.60$, resulting in a substantial content of 0.4 holes
per oxygen and a near Ni$(3d^{8})$ configuration. Thus maybe nonsurprisingly
for this high-pressure system, the degree of covalency is significant and
about 1.8 formula-unit-cell valence electrons have to reside in delocalized
states of La$(5d6s)$ and/or Ni$(4s)$ character. The DFT+sic computed charge
transfer energy $\Delta=3.6$ eV lies inbetween the infinite-layer values of
5.0 eV for NdNiO2 and of 1.3 eV for SrCuO2 [22]. Note in this respect that
while the Ni-$e_{g}$ character dominates at $\varepsilon_{\rm F}$, the
corresponding O$(2p)$ content is still larger than in NdNiO2, hinting to a
comparatively enhanced role of oxygen degrees of freedom at low-energy.

Figure 2: FL-like vs. NFL behavior with $T$ in the $Fmmm$ phase from
DFT+sicDMFT. (a) Top panel: total spectral function (inset: wide energy
scale). Middle panels: local Ni-$e_{g}$ spectral function (inset: wide energy
scale including local Ni-$t_{2g}$ spectra for $T=80$ K). Bottom panel:
Ni-$e_{g}$ hybridization function (inset: same protocol as for middle panels).
(b) Imaginary part of the Ni-$e_{g}$ self-energies $\Sigma(i\omega_{n})$ on
the Matsubara axis (colors according to middle panel of (a)).  Figure 3:
k-dependent spectral features of the $Fmmm$ phase from DFT (a,b) and
DFT+sicDMFT at $T=80$ K (c-f). (a) DFT band structure in Ni$(3d)$-fatspec
representation and Brillouin zone with high-symmetry lines. (b) DFT $k_{z}=0$
Fermi surface in Ni-fatspec representation. (c) Interacting spectral function
in large energy window, and (d) in small window around $\varepsilon_{\rm F}$.
(e) Same as (d) but in Ni-fatspec representation. (f) Interacting $k_{z}=0$
Fermi surface from left to right: regular intensity, 10$\times$-increased
intensity, Ni-fatspec representation.

As anticipated for a bilayer oxide, and already documented by the subtle
$Amam$ vs. $Fmmm$ low-energy difference, the electronic states near the Fermi
level are delicate. To reach a better understanding of the relevant coherence
scales, we therefore performed additional calculations at higher $T$. Figure 2
shows that the low-energy regime reacts sensitively to temperature. The total
spectral function (top panel of Fig. 2a) evolves from a QP-like structure at
room temperature to a flattened weight around the Fermi level at $T=80$ K. On
the local level (middle panel of Fig. 2a), it is first noted that
Ni-$d_{z^{2}}$ and Ni-$d_{x^{2}-y^{2}}$ are both half filled. Second, upon
lowering $T$, the near-$\varepsilon_{\rm F}$ $A_{\rm loc}(\omega)$ transforms
from QP-like, to pseudogap and finally to low-amplitude peak. A link between
total and local spectrum may be established via the hybridization function
$-{\rm Im}\Delta_{\rm h}$, shown in bottom panel of Fig. 2a. It displays a
(pseudo)gap at $T=80$ K, altogether rendering it obvious that a NFL state is
reached. The Ni-$e_{g}$ self-energies shown in Fig. 2b underline this picture
with a low-frequency upturn at low $T$. A Fermi-liquid (FL) fit to the
room-$T$ data yields effective masses $m^{*}_{z^{2}}=6.4$ and
$m^{*}_{x^{2}-y^{2}}=5.6$. But note that even the ambient system is not a good
FL. Though a linear-in-frequency regime holds for smallest Matsubara
$\omega_{n}$, the scattering rate $\sim-\lim_{\omega_{n}\rightarrow
0}\Sigma(i\omega_{n})$ remains quite large for the given $T$. And already the
$T=150$ K data displays further NFL tendencies. Because of half-filled
Ni-$e_{g}$ (one electron in each of the two orbitals) as well as the very low-
energy scale for non-QP formation, a sole $U$\- or $J_{\rm H}$-driven NFL
behavior seems not likely.

Let us thus turn to the k-resolved spectra to gain insight into the origin of
the NFL behavior. Figure 3(a,b) displays the DFT band structure and Fermi
surface to set the stage. The Ni fatspec representation [25] marks the
dominant Ni$(3d)$ character, showing that there are majorly four Ni-$e_{g}$
dispersions, associated with the two equivalent Ni sites in the unit cell,
governing the low-energy region. The inner two bands form a ($\alpha$)
electron pocket around $\Gamma$ and a ($\beta$) hole sheet opening towards
$X$. Note that those two Fermi-surface sheets are strongly
Ni-$d_{z^{2}}/d_{x^{2}-y^{2}}$ mixed. The upper (anti-bonding)
Ni-$d_{z^{2}}$-dominated band, also having sizable apical O$(2p)$ character,
is not crossing $\varepsilon_{\rm F}$. Instead, a self-doping mainly La-based
band mingles into the Ni-$e_{g}$ fourfold dispersion and gives rise to a large
electron pocket around $Z$. Finally, the lower (bonding)
Ni-$d_{z^{2}}$-dominated band forms flattened ($\gamma$) hole pockets around
$M$. The Ni-$t_{2g}$ character very weakly mixes into part of the Fermi-
surface sheets, but otherwise has major weight below $\varepsilon_{\rm F}$ and
does not play a key role for the low-energy physics.

With correlations at $T=80$ K, i.e. deep into the NFL regime, the low-energy
picture changes quite dramatically (see Fig. 3c-f). First, the
near-$\varepsilon_{\rm F}$ dispersions become generally very weak in
intensity. Figure 3d shows that while the dispersions away from the Fermi
level still keep some renormalized coherence, right at $\varepsilon_{\rm F}$
they dissolve. This is orthogonal to the understanding of a FL state and marks
the strong NFL nature of the pressurized bilayer system. Accordingly, the
(weakly $k_{z}$-dependent) interacting Fermi surface displayed in Fig. 3e
becomes very weak and blurry. Only when raising the representation intensity
(middle part of Fig. 3e), a hole-like ”sheet” structure around $M$ emerges. It
is reminiscent of the original DFT Ni-$d_{z^{2}}$ flat-band-based $\gamma$
sheet, but with stronger mixed Ni-$e_{g}$ character (cf. right part of Fig.
3e)). The enhanced intermixing presumably comes from a correlation-induced
meet up of with the Ni-$d_{x^{2}-y^{2}}$ branded $\beta$ sheet in $\Gamma$-$M$
direction. This is also supported from the disappearance of the self-doping
$Z$ pocket from the Fermi surface. Such strong correlation-induced shifts of
self-doping bands has already been observed in other nickelates [25, 60, 26].
It becomes intuitively obvious that all these very-low-energy ramifications in
the interacting regime have to strongly build up on the present flat-band
scenario. There, the introduced renormalizations create a large phase space
for intriguing quantum fluctuations, rendering robust QP formation impossible
[61].

Albeit the NFL character may be relevant for superconductivity, let us get a
first handle on SC instabilities from a weak-coupling perspective for coherent
Fermi-surface sheets and leave the discussion of the role of the NFL behavior
and its relevance for superconductivity to future studies. To this, we
constructed a $4\times 4$ maximally-localized Wannier [62] Hamiltonian for the
Ni-$e_{g}$ based DFT bands. It carries hopping integrals
$t_{ij}^{\ell\ell^{\prime}}$ for
$\ell,\ell^{\prime}=1d_{z^{2}},1d_{x^{2}-y^{2}},2d_{z^{2}},2d_{x^{2}-y^{2}}$
and lattice sites $i,j$. Here $1$ and $2$ refer to the two Ni sublattices in
the $Fmmm$ structure.

Figure 4: RPA results at $T=80$ K. (a) Spin susceptibility
$\chi_{\mathrm{S}}(\bf{q},\omega=0)$ (in 1/eV) for $U=1.5$ eV and $J_{\rm
H}=U/7$. (b) Leading solution of the linearized gap equation for same
interaction strength at $k_{z}=0$. (c) Phase diagram of the leading eigenvalue
solution. Level lines mark the eigenvalue magnitude, and the green dot case
(a,b). (d) Sketch for hole doping (red arrows) near the SC regime, relating
the flat-band scenarios of pressurized bilayer and reduced nickelates (see
text).

Adding local Coulomb interactions, the effective Hamiltonian reads

$H=\sum_{i\neq
j,\ell\ell^{\prime},\sigma}t_{ij}^{\ell\ell^{\prime}}\,c^{\dagger}_{i\ell\sigma}c^{\hfill}_{j\ell^{\prime}\sigma}+\sum_{i}\left(H^{(i)}_{\rm
int}+H^{(i)}_{\rm orb}\right),$ (1)

The onsite interaction $H^{(i)}_{\rm int}$ has Slater-Kanamori form, i.e.
includes density-density terms as well as pair-hopping and spin-flip terms,
parameterized again by $U$ and $J_{\rm H}$. Note that within the downfolded
model, the $U$ value is smaller than in the comprehensive DFT+sicDMFT
treatment. The remaining non-interacting onsite $H^{(i)}_{\rm orb}$ carries
crystal-field terms via onsite levels $\varepsilon_{\ell}$. In order to
investigate the pairing symmetry, we employ a diagrammatic, i.e. RPA treatment
of interactions starting from the effective band structure. This procedure
provides a reliable first insight into the leading pairing instabilities in
layered unconventional superconductors, including nickelates [56, 29]. Further
details are presented in the supplementary material.

We calculate the RPA spin (charge) susceptibility $\chi_{S}$ $(\chi_{C})$
based on a two-dimensional cut through first and second Brillouin zone (BZ) to
automatically include important scatterings from the $k_{z}=0$ and $k_{z}=1/2$
sheets of the BZ. Calculating $\chi_{S}(\mathbf{q},\omega=0)$ at $T=80$ K and
noting the flatness of the $M$-pocket (see Fig. 3(a)), the real part resembles
band nesting and is not restricted only to the Fermi surface. The resulting
peak structure is shown in Fig. 4a, with similarity to previous findings [57].
By solving the linearized gap equation for a range of $U$\- and $J_{\rm
H}$-values we plot the leading SC instability in Fig. 4c. The largest and most
likely part of the parameter space is occupied by a solution with a cuprate-
like $B_{1g}$ instability of $d_{x^{2}-y^{2}}$ symmetry on $\beta$ and
$\gamma$ bands. This solution is presented in Fig. 4b. At the same time, the
origin of the $d_{x^{2}-y^{2}}$ gap is not directly related to the
Ni-$d_{x^{2}-y^{2}}$ orbital. Observe that the gap magnitude is about 3 times
higher in the $\gamma$ band compared to $\alpha$ and $\beta$ band,
highlighting the importance of the flat Ni-$d_{z^{2}}$ band around the $M$
point. For lower interactions, a solution with $d_{xy}$ symmetry is dominant,
whereas a triplet solution does not become leading within the inspected
parameter space. The less likely solutions are discussed in the supplementary
material. Note that the $s_{\pm}$-wave solution is also a sub-leading one. The
total pairing strength given by the leading eigenvalue is strongly driven by
$U$, and generally, superconductivity would sensitively react to the level of
of Ni-$d_{x^{2}-y^{2}}$/$d_{z^{2}}$ incoherence. This issue should be a
subject of further theoretical and experimental study.

Discussion.— We have shown that the peculiar correlated electronic structure
and concomitant SC instability of La3Ni2O7 originates from the interplay of
half-filled Ni-$e_{g}$ orbitals within a Ni-$d_{z^{2}}$-created flat-band
scenario. The role of O$(2p)$ is enhanced compared to reduced nickelates, yet
a decisive role cannot be deduced from this initial theory study. But note
that the Ni-O distance along $c$ within the bilayer turns out remarkably
small, thus inter-site Ni-Ni self-energies may not be negligible. Those could
e.g. be addressed in a two-site, two-orbital cluster-DMFT study, which however
is beyond the present scope. Comparing to the phenomenology of reduced
nickelates, a line can be drawn between these and the bilayer system as
sketched in Fig. 4d for a modelized single-Ni unit-cell system: In the reduced
systems, the hole doping relevant for superconductivity occurs mainly in the
Ni-$d_{z^{2}}$ upper-branch flat-band part around $k_{z}=1/2$ [23, 24, 25,
26], corresponding to a well-filled Ni-$d_{z^{2}}$ orbital. For pressurized
La3Ni2O7 however, the hole doping takes place in the lower-branch flat-band
part around $k_{z}=0$. There, Ni-$d_{z^{2}}$ is close to half-filling, much
more correlated and therefore more on par with Ni-$d_{x^{2}-y^{2}}$. This
should be the reason for the different $T_{\rm c}$ in the unlike nickelates.
One may speculate that this different flat-band doping regime can also be
realized in reduced multilayer nickelates via tailored doping protocols even
at ambient pressure.

##  I Acknowledgements

The work is supported by the German Research Foundation within the bilateral
NSFC-DFG Project ER 463/14-1. Computations were performed at the Ruhr-
University Bochum and the JUWELS Cluster of the Jülich Supercomputing Centre
(JSC) under project miqs.

## References

  * Sun  _et al._ [2023] H. Sun, M. Huo, X. Hu, J. Li, Y. Han, L. Tang, Z. Mao, P. Yang, B. Wang, J. Cheng, D.-X. Yao, G.-M. Zhang, and M. Wang, (2023), arXiv:2305.09586 .
  * Li  _et al._ [2019] D. Li, K. Lee, B. Y. Wang, M. Osada, S. Crossley, H. R. Lee, Y. Cui, Y. Hikita, and H. Y. Hwang, Nature 572, 624 (2019).
  * Li  _et al._ [2020] D. Li, B. Y. Wang, K. Lee, S. P. Harvey, M. Osada, B. H. Goodge, L. F. Kourkoutis, and H. Y. Hwang, Phys. Rev. Lett. 125, 027001 (2020).
  * Zeng  _et al._ [2020] S. Zeng, C. S. Tang, X. Yin, C. Li, M. Li, Huiong, W. Liu, G. J. Omar, H. Jani, Z. S. Lim, K. Han, D. Wan, P. Yang, S. J. Pennycook, A. T. S. Wee, and A. Ariando, Phys. Rev. Lett. 125, 147003 (2020).
  * Osada  _et al._ [2020] M. Osada, B. Y. Wang, B. H. Goodge, K. Lee, H. Yoon, K. Sakuma, D. Li, M. Miura, L. F. Kourkoutis, and H. Y. Hwang, Nano Letters 20, 5735 (2020).
  * Osada  _et al._ [2021] M. Osada, B. Y. Wang, B. H. Goodge, S. P. Harvey, K. Lee, D. Li, L. F. Kourkoutis, and H. Y. Hwang, Advanced Materials 33, 2104083 (2021).
  * Pan  _et al._ [2021] G. A. Pan, D. F. Segedin, H. LaBollita, Q. Song, E. M. Nica, B. H. Goodge, A. T. Pierce, S. Doyle, S. Novakov, D. C. Carrizales, A. T. N’Diaye, P. Shafer, H. Paik, J. T. Heron, J. A. Mason, A. Yacoby, L. F. Kourkoutis, O. Erten, C. M. Brooks, A. S. Botana, and J. A. Mundy, Nature Materials 21, 160 (2021).
  * Zeng  _et al._ [2022] S. Zeng, C. Li, L. E. Chow, Y. Cao, Z. Zhang, C. S. Tang, X. Yin, Z. S. Lim, J. Hu, P. Yang, and A. Ariando, Science Advances 8, eabl9927 (2022).
  * Wu  _et al._ [2020] X. Wu, D. Di Sante, T. Schwemmer, W. Hanke, H. Y. Hwang, S. Raghu, and R. Thomale, Phys. Rev. B 101, 060504 (2020).
  * Zhang  _et al._ [2020] G.-M. Zhang, Y.-f. Yang, and F.-C. Zhang, Phys. Rev. B 101, 020501 (2020).
  * Karp  _et al._ [2020] J. Karp, A. S. Botana, M. R. Norman, H. Park, M. Zingl, and A. Millis, Phys. Rev. X 10, 021061 (2020).
  * Leonov  _et al._ [2020] I. Leonov, S. L. Skornyakov, and S. Y. Savrasov, Phys. Rev. B 101, 241108 (2020).
  * Adhikary  _et al._ [2020] P. Adhikary, S. Bandyopadhyay, T. Das, I. Dasgupta, and T. Saha-Dasgupta, Phys. Rev. B 102, 100501 (2020).
  * Kitatani  _et al._ [2020] M. Kitatani, L. Si, O. Janson, R. Arita, Z. Zhong, and K. Held, npj Quantum Materials 5, 10.1038/s41535-020-00260-y (2020).
  * Been  _et al._ [2021] E. Been, W.-S. Lee, H. Y. Hwang, Y. Cui, J. Zaanen, T. Devereaux, B. Moritz, and C. Jia, Phys. Rev. X 11, 011050 (2021).
  * Geisler and Pentcheva [2021] B. Geisler and R. Pentcheva, Phys. Rev. Res. 3, 013261 (2021).
  * Gu  _et al._ [2020] Y. Gu, S. Zhu, X. Wang, J. Hu, and H. Chen, Communications Physics 3, 10.1038/s42005-020-0347-x (2020).
  * Plienbumrung  _et al._ [2022] T. Plienbumrung, M. Daghofer, M. Schmid, and A. M. Oleś, Phys. Rev. B 106, 134504 (2022).
  * Jiang  _et al._ [2020] M. Jiang, M. Berciu, and G. A. Sawatzky, Phys. Rev. Lett. 124, 207004 (2020).
  * Held  _et al._ [2022] K. Held, L. Si, P. Worm, O. Janson, R. Arita, Z. Zhong, J. M. Tomczak, and M. Kitatani, Frontiers in Physics 9, 10.3389/fphy.2021.810394 (2022).
  * Werner and Hoshino [2020] P. Werner and S. Hoshino, Phys. Rev. B 101, 041104 (2020).
  * Lechermann [2020a] F. Lechermann, Phys. Rev. B 101, 081110 (2020a).
  * Lechermann [2020b] F. Lechermann, Phys. Rev. X 10, 041002 (2020b).
  * Lechermann [2021] F. Lechermann, Phys. Rev. Mater. 5, 044803 (2021).
  * Lechermann [2022a] F. Lechermann, Phys. Rev. B 105, 155109 (2022a).
  * Lechermann [2022b] F. Lechermann, Electronic Structure 4, 015005 (2022b).
  * Petocchi  _et al._ [2020] F. Petocchi, V. Christiansson, F. Nilsson, F. Aryasetiawan, and P. Werner, Phys. Rev. X 10, 041047 (2020).
  * Kang and Kotliar [2021] C.-J. Kang and G. Kotliar, Phys. Rev. Lett. 126, 127401 (2021).
  * Kreisel  _et al._ [2022] A. Kreisel, B. M. Andersen, A. T. Rømer, I. M. Eremin, and F. Lechermann, Phys. Rev. Lett. 129, 077002 (2022).
  * Demourgues  _et al._ [1993] A. Demourgues, F. Weill, B. Darriet, A. Wattiaux, J. Grenier, P. Gravereau, and M. Pouchard, Journal of Solid State Chemistry 106, 330 (1993).
  * Mizokawa  _et al._ [2000] T. Mizokawa, D. I. Khomskii, and G. A. Sawatzky, Phys. Rev. B 61, 11263 (2000).
  * Park  _et al._ [2012] H. Park, A. J. Millis, and C. A. Marianetti, Phys. Rev. Lett. 109, 156402 (2012).
  * Johnston  _et al._ [2014] S. Johnston, A. Mukherjee, I. Elfimov, M. Berciu, and G. A. Sawatzky, Phys. Rev. Lett. 112, 106404 (2014).
  * Subedi  _et al._ [2015] A. Subedi, O. E. Peil, and A. Georges, Phys. Rev. B 91, 075128 (2015).
  * Bisogni  _et al._ [2016] V. Bisogni, S. Catalano, R. J. Green, M. Gibert, R. Scherwitzl, Y. Huang, V. N. Strocov, P. Zubko, S. Balandeh, J.-M. Triscone, G. Sawatzky, and T. Schmitt, Nature Communications 7, 10.1038/ncomms13017 (2016).
  * Perry  _et al._ [2001] R. S. Perry, L. M. Galvin, S. A. Grigera, L. Capogna, A. J. Schofield, A. P. Mackenzie, M. Chiao, S. R. Julian, S. I. Ikeda, S. Nakatsuji, Y. Maeno, and C. Pfleiderer, Phys. Rev. Lett. 86, 2661 (2001).
  * Borzi  _et al._ [2007] R. A. Borzi, S. A. Grigera, J. Farrell, R. S. Perry, S. J. S. Lister, S. L. Lee, D. A. Tennant, Y. Maeno, and A. P. Mackenzie, Science 315, 214 (2007).
  * Lechermann  _et al._ [2020] F. Lechermann, Q. Han, and A. J. Millis, Phys. Rev. Res. 2, 033490 (2020).
  * King  _et al._ [2013] P. D. C. King, T. Takayama, A. Tamai, E. Rozbicki, S. M. Walker, M. Shi, L. Patthey, R. G. Moore, D. Lu, K. M. Shen, H. Takagi, and F. Baumberger, Phys. Rev. B 87, 241106 (2013).
  * Zhang  _et al._ [1994] Z. Zhang, M. Greenblatt, and J. Goodenough, Journal of Solid State Chemistry 108, 402 (1994).
  * Taniguchi  _et al._ [1995] S. Taniguchi, T. Nishikawa, Y. Yasui, Y. Kobayashi, J. Takeda, S. ichi Shamoto, and M. Sato, Journal of the Physical Society of Japan 64, 1644 (1995).
  * Wu  _et al._ [2001] G. Wu, J. J. Neumeier, and M. F. Hundley, Phys. Rev. B 63, 245120 (2001).
  * Seo  _et al._ [1996] D.-K. Seo, W. Liang, M.-H. Whangbo, Z. Zhang, and M. Greenblatt, Inorganic Chemistry 35, 6396 (1996).
  * Pardo and Pickett [2011] V. Pardo and W. E. Pickett, Phys. Rev. B 83, 245128 (2011).
  * Lechermann  _et al._ [2019] F. Lechermann, W. Körner, D. F. Urban, and C. Elsässer, Phys. Rev. B 100, 115125 (2019).
  * Grieger  _et al._ [2012] D. Grieger, C. Piefke, O. E. Peil, and F. Lechermann, Phys. Rev. B 86, 155121 (2012).
  * Körner and Elsässer [2010] W. Körner and C. Elsässer, Phys. Rev. B 81, 085324 (2010).
  * Elsässer  _et al._ [1990] C. Elsässer, N. Takeuchi, K. M. Ho, C. T. Chan, P. Braun, and M. Fahnle, Journal of Physics: Condensed Matter 2, 4371 (1990).
  * Lechermann  _et al._ [2002] F. Lechermann, F. Welsch, C. Elsässer, C. Ederer, M. Fähnle, J. M. Sanchez, and B. Meyer, Phys. Rev. B 65, 132104 (2002).
  * Meyer  _et al._ [1998] B. Meyer, C. Elsässer, F. Lechermann, and M. Fähnle,  _FORTRAN 90 Program for Mixed-Basis-Pseudopotential Calculations for Crystals_ , Max-Planck-Institut für Metallforschung, Stuttgart (1998).
  * Werner  _et al._ [2006] P. Werner, A. Comanac, L. de’ Medici, M. Troyer, and A. J. Millis, Phys. Rev. Lett. 97, 076405 (2006).
  * Parcollet  _et al._ [2015] O. Parcollet, M. Ferrero, T. Ayral, H. Hafermann, I. Krivenko, L. Messio, and P. Seth, Computer Physics Communications 196, 398 (2015).
  * Seth  _et al._ [2016] P. Seth, I. Krivenko, M. Ferrero, and O. Parcollet, Comput. Phys. Commun. 200, 274 (2016).
  * Amadon  _et al._ [2008] B. Amadon, F. Lechermann, A. Georges, F. Jollet, T. O. Wehling, and A. I. Lichtenstein, Phys. Rev. B 77, 205112 (2008).
  * Anisimov  _et al._ [1993] V. I. Anisimov, I. V. Solovyev, M. A. Korotin, M. T. Czyżyk, and G. A. Sawatzky, Phys. Rev. B 48, 16929 (1993).
  * Graser  _et al._ [2009] S. Graser, T. A. Maier, P. J. Hirschfeld, and D. J. Scalapino, New Journal of Physics 11, 025016 (2009).
  * Luo  _et al._ [2023] Z. Luo, X. Hu, M. Wang, W. Wú, and D.-X. Yao, (2023), arXiv:2305.15564 .
  * Jarrell and Gubernatis [1996] M. Jarrell and J. E. Gubernatis, Physics Reports 269, 133 (1996).
  * Vidberg and Serene [1977] H. J. Vidberg and J. W. Serene, J Low Temp Phys 29, 179 (1977).
  * Worm  _et al._ [2022] P. Worm, L. Si, M. Kitatani, R. Arita, J. M. Tomczak, and K. Held, Phys. Rev. Mater. 6, L091801 (2022).
  * Sayyad  _et al._ [2020] S. Sayyad, E. W. Huang, M. Kitatani, M.-S. Vaezi, Z. Nussinov, A. Vaezi, and H. Aoki, Phys. Rev. B 101, 014501 (2020).
  * Marzari  _et al._ [2012] N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Rev. Mod. Phys. 84, 1419 (2012).

Supplemental Material

##  II Details for the DFT+sicDMFT calculations

For the DFT part, a mixed-basis pseudopotential framework in the local-density
approximation (LDA) is put into practise. A $9\times 9\times 9$ for $Fmmm$ and
a $5\times 5\times 5$ $k$-point mesh for $Amam$ is utilized and the plane-wave
cutoff energy is generally set to $E_{\rm cut}=16$ Ry. Local basis orbitals
are introduced for La$(5d)$, Ni$(3d)$ and O$(2s,2p)$. The role of possible
spin-orbit effects is neglected in the crystal calculations. The DMFT
correlated subspace on each Ni site is governed by a rotational-invariant
five-orbital Slater Hamiltonian. It is applied to the Ni$(3d)$ projected-local
orbitals. The projection is performed on the Kohn-Sham bands above the 7(14)
bands of dominant O$(2s)$ character for $Fmmm$($Amam$). The projection window
spans 34(68) bands, including the KS states of dominant Ni$(3d)$ and O$(2p)$,
as well as one additional KS band for each La site in the primitive cell. A
Hubbard $U=10$ eV and a Hund exchange $J_{\rm H}=1.0$ eV are chosen to
parametrize the local interacting Hamiltonian. No explicit Coulomb
interactions beyond DFT are introduced on La. The fully-localized-limit
double-counting scheme [1] is applied. Continuous-time quantum Monte Carlo in
hybridization expansion as implemented in the TRIQS code is used to solve the
multisite DMFT problem. Up to $1.5\cdot 10^{9}$ Monte-Carlo sweeps are
performed to reach convergence. A Matsubara mesh of 1025(2049) frequencies is
used to account for the higher(lower)-temperature regime. In detail, this
means that 1025 frequencies are used for $T>100$ K and 2049 frequencies are
utilized for $T<100$ K. For the analytical continuation from Matsubara space
onto the real-frequency axis, the Maximum-entropy method [2] is used for the
k-integrated spectra (by continuation of the Bloch Green’s function) and the
Padé method [3] is employed for the k-resolved spectra (by continuation of the
local self-energies). Note that in the present work, the DMFT hybridization
function, reading $\Delta_{\rm hyb}(\omega)=\omega+\mu-G_{\rm
loc}^{-1}(\omega)-\Sigma_{\rm imp}(\omega)$, is constructed from using the
Padé-method continued local self-energies.

The maximally-localized Wannier construction was performed by initial
projection onto the Ni-$e_{g}$ orbitals in the $Fmmm$ primitive cell. The
comparison of the downfolded low-energy Wannier bands with full LDA bands is
depicted in Fig. 5b along with the original LDA Fermi surface in Fig. 5a. Note
that the La-based self-doping band forms electron pockets around the $Z$
point. Those pockets are excluded from the Wannier construction, which is
justified by their disappearance in the interacting regime described by the
comprehensive DFT+sicDMFT approach.

Figure 5: Additional DFT data for the $Fmmm$ phase. (a) Fermi surface with
$\alpha,\beta$ sheets (red), $\gamma$ sheet (magenta) and the La-based self-
doping sheet (blue). (b) Complete bandstructure (lightblue) and downfolded
four-band Ni-$e_{g}$ based maximally-localized Wannier dispersion (black) for
the low-energy region.

##  III Calculation of RPA susceptibilities and solution linearized gap
equation

We begin by calculating the multiorbital charge and spin susceptibilities in
the random phase approximation, which enter the pairing vertex for the
linearized gap equation for the paramagnetic case. This procedure has first
been used in ref. [4]. In Matsubara formalism the bare susceptibility can be
written in terms of non-interacting Green’s functions

$G_{l_{1},l_{2}}(\mathbf{k},i\omega_{n})=\sum_{j}\frac{a_{l_{1}}(j,\mathbf{k})a_{l_{2}}^{*}(j,\mathbf{k})}{i\omega_{n}-\epsilon_{j\mathbf{k}}},$
(2)

where $a_{l}(j,\mathbf{k})$ denotes the matrix elements of the unitary
transformation from orbital to band space of the non-interacting part of the
Hamiltonian given in eq. 1 of the main text. It connects the the j-th band
$\epsilon_{j}$ with the l-th orbital. After performing the sum over the
Matsubara frequencies, the analytical continuation of the bare susceptibility
follows as

$\displaystyle(\chi_{0})_{l_{1}l_{4}}^{l_{2}l_{3}}(\mathbf{q},\omega+i\delta)$
$\displaystyle=\frac{1}{N}\sum_{\mathbf{k},ij}a_{l_{1}}(j,\mathbf{k}+\mathbf{q})a^{*}_{l_{2}}(j,\mathbf{k}+\mathbf{q})$
(3)
$\displaystyle\times{a}_{l_{3}}(i,\mathbf{k}){a}^{*}_{l_{4}}(i,\mathbf{k})\frac{f(\epsilon_{i\mathbf{k}})-f(\epsilon_{j\mathbf{k}+\mathbf{q}})}{(\epsilon_{j\mathbf{k}+\mathbf{q}}-\epsilon_{i\mathbf{k}})-\omega-i\delta}.$

The on-site interaction part of the Hamiltonian includes intraorbital
(interorbital) Hubbard repulsion $U$ ($U^{\prime}$), Hund’s coupling $J$ and
pair-hopping $J^{\prime}$:

$\displaystyle H^{(i)}_{int}$
$\displaystyle=\frac{U}{2}\sum_{il\sigma}c^{\dagger}_{il\sigma}c^{\dagger}_{il\bar{\sigma}}c_{il\bar{\sigma}}c_{il{\sigma}}$
(4)
$\displaystyle+\,U^{\prime}\sum_{il<l^{\prime}\sigma\sigma^{\prime}}c^{\dagger}_{il\sigma}c^{\dagger}_{il^{\prime}\sigma^{\prime}}c_{il^{\prime}{\sigma^{\prime}}}c_{il{\sigma}}$
$\displaystyle+\,J\sum_{il<l^{\prime}\sigma\sigma^{\prime}}c^{\dagger}_{il^{\prime}\sigma}c^{\dagger}_{il\sigma^{\prime}}c_{i^{\prime}\sigma^{\prime}}c_{il{\sigma}}$
$\displaystyle+\,\frac{J^{\prime}}{2}\sum_{il\neq
l^{\prime}\sigma}c^{\dagger}_{il\sigma}c^{\dagger}_{il\bar{\sigma}}c_{il^{\prime}\bar{\sigma}}c_{il^{\prime}{\sigma}},$

with the creation operator $c^{\dagger}_{l\sigma}$ creating an quasiparticle
with orbital $l$ and spin $\sigma$. Assuming spin-rotational invariance for
the interactions enforces two relations $U^{\prime}=U-2J$ and $J^{\prime}=J$.
Thus, the interaction can be parameterized in terms of $U$ and $J$. Using the
RPA formalism different orbital-matrix Dyson equations for charge and spin
susceptibilities, which are for a spin-diagonal non-interacting Hamiltonian
equivalent on the bare level, arise:

$\displaystyle\chi_{S/C}$
$\displaystyle=\left[1\mp\chi_{0}U_{S/C}\right]^{-1}\chi_{0}.$ (5)

The physical observables follow from an orbital contraction of
$l_{1}$($l_{2}$) with $l_{4}$($l_{3}$). The final contraction of the orbitals
can be done orbital resolved. The on-site interaction matrices in the Dyson
equations for orbitals residing on the same atom are given by

$\displaystyle(U_{C})_{l_{2}l_{3}}^{l_{1}l_{4}}=\begin{cases}U&l_{1}=l_{2}=l_{3}=l_{4}\\\
2J-U^{\prime}&l_{1}=l_{2}\neq l_{3}=l_{4}\\\ 2U^{\prime}-J&l_{1}=l_{4}\neq
l_{2}=l_{3}\\\ J^{\prime}&l_{1}=l_{3}\neq l_{2}=l_{4}\\\
0&otherwise.\end{cases}$ (6)
$\displaystyle(U_{S})_{l_{2}l_{3}}^{l_{1}l_{4}}=\begin{cases}U&l_{1}=l_{2}=l_{3}=l_{4}\\\
U^{\prime}&l_{1}=l_{2}\neq l_{3}=l_{4}\\\ J&l_{1}=l_{4}\neq l_{2}=l_{3}\\\
J^{\prime}&l_{1}=l_{3}\neq l_{2}=l_{4}\\\ 0&otherwise.\end{cases}$ (7)

In our calculation the summation in eq. 3 is performed at constant $k_{z}$.
However, we use the first and second Brillouin zones (BZ) for integration such
that the integration includes the $k_{z}=0$ sheet and the $k_{z}=\pi$ sheet of
the BZ which follows from the face-centered cubic structure of the lattice. We
justify this treatment below. For the solution of the linearized gap equation
presented in the main work the integration was performed on a $320\times 320$
$k$-mesh together with a small broadening $\delta=10$ meV for numerical
purposes.

To ensure that we stay in the paramagnetic regime, we first checked the
parameters for which the matrix $(1-\chi_{0}U_{S})$ has zero eigenvalues,
which corresponds to an antiferromagnetic instability. The location of the
instability is a curve in $U$-$J$-parameter space. We therefore show the spin
susceptibility as a function of $U$ with fixed $U/J$ and as function of $J/U$
for fixed $U$ in Fig. 6(a) and Fig. 6(b) at the strongest peak, located at
$\mathbf{q}/\pi\approx(0,0.45)$. The instability is driven by mixed
$d_{z^{2}}$ and $d_{x^{2}-y^{2}}$ and pure $d_{z^{2}}$ orbital contributions.
The importance of the mixed contributions is tuned with $J$ and can become
dominant in the high $J$ limit.

By showing Fig. 6(c) we justify that we use two-dimensional integration by
comparing the charge and spin susceptibility with those in Fig. 6(d) for which
the full three-dimensional integration is performed along the same path
(adding additionally $\Gamma$-$Z$ for completeness). Note that due to
symmetry, the $Z$ point and the $(\pi,\pi,0)$ point are identical and that the
chosen path is not inside the first BZ. All essential peaks agree, only the
peak at $(\pi,0,0)$ is enhanced but not becoming too important.

The orbital resolved contributions to the spin susceptibility presented in the
constant $k_{z}$-plane are shown in Fig. 7. For the chosen values of $U=1.5$
eV and $J=U/7$, the strongest contributions originate again from the
$d_{z^{2}}$ orbital. The dominance of the $d_{z^{2}}$ orbital over
$d_{x^{2}-y^{2}}$ is likely due to the flatness of the $d_{z^{2}}$ pocket
around the $M$ point. To really decipher the role of the orbitals, however,
interactions should be taken into account for the band structure. The shape of
the orbital resolved spin susceptibility is agreeing with the results obtained
in Ref. [5]. For the two dimensional cuts a $256\times 256$ $k$-mesh is used
for the integration.

Figure 6: RPA susceptibility $(\chi(\mathbf{q},\omega=0))$ calculated at
$T=80$ K. (a) Spin susceptibility and its orbital resolved contributions
calculated at $q/\pi\approx(0,0.45)$ for $U=1.5$ eV and $J=U/7$. (b) As in (a)
but as function of $J$ for constant $U=1.35$. (c) and (d) Comparison of spin
and charge susceptibilities $\chi_{S}$ and $\chi_{C}$ calculated from two-
dimensional and three-dimensional approach, respectively.

The fully momentum-symmetrized (-antisymmetrized) interaction vertex in the
singlet (triplet) channel reads

$\displaystyle(\Gamma_{\mathrm{s}})_{l_{1}l_{4}}^{l_{2}l_{3}}(\mathbf{k},\mathbf{k^{\prime}})$
$\displaystyle=\left[\frac{3}{2}U_{S}\frac{\chi_{S}(\mathbf{k}+\mathbf{k^{\prime}})+\chi_{S}(\mathbf{k}-\mathbf{k^{\prime}})}{2}U_{S}\right.$
(8)
$\displaystyle\left.-\frac{1}{2}\frac{U_{C}(\chi_{C}(\mathbf{k}+\mathbf{k^{\prime}})+\chi_{C}(\mathbf{k}-\mathbf{k^{\prime}})}{2}U_{C}\right.$
$\displaystyle\left.+\frac{1}{2}U_{S}+\frac{1}{2}U_{C}\right]_{l_{1}l_{4}}^{l_{2}l_{3}}$
$\displaystyle(\Gamma_{\mathrm{t}})_{l_{1}l_{4}}^{l_{2}l_{3}}(\mathbf{k},\mathbf{k^{\prime}})$
$\displaystyle=\left[-\frac{1}{2}U_{S}\frac{\chi_{S}(\mathbf{k}+\mathbf{k^{\prime}})-\chi_{S}(\mathbf{k}-\mathbf{k^{\prime}})}{2}U_{S}\right.$
(9)
$\displaystyle\left.-\frac{1}{2}U_{C}\frac{\chi_{C}(\mathbf{k}+\mathbf{k^{\prime}})-\chi_{C}(\mathbf{k}-\mathbf{k^{\prime}})}{2}U_{C}\right.$
$\displaystyle\left.-\frac{1}{2}U_{S}+\frac{1}{2}U_{C}\right]_{l_{1}l_{4}}^{l_{2}l_{3}}$

and projected back to band space via unitary transform

$\displaystyle(\Gamma_{s/t})_{\mu\nu}(\mathbf{k},\mathbf{k^{\prime}})$
$\displaystyle=\mathrm{Re}\sum_{l_{1}l_{2}l_{3}l_{4}}a_{l_{1}}(\nu,\mathbf{k})a^{*}_{l_{2}}(\mu,\mathbf{k^{\prime}})$
(10)
$\displaystyle\times{a}_{l_{3}}(\nu,-\mathbf{k}){a}^{*}_{l_{4}}(\mu,-\mathbf{k^{\prime}})(\Gamma_{\mathrm{s/t}})_{l_{1}l_{4}}^{l_{2}l_{3}}(\mathbf{k},\mathbf{k^{\prime}}).$

The linearized gap equation

$\displaystyle\lambda_{i}^{s/t}\Delta_{i}^{s/t}(\mu,\mathbf{k})$
$\displaystyle=\left[-\frac{1}{V_{\mathrm{BZ}}}\sum_{\nu}\right.$ (11)
$\displaystyle\left.\int_{\mathrm{FS_{\nu}}}\mathrm{d}S^{\prime}\frac{(\Gamma_{s/t})_{\mu\nu}(\mathbf{k,k^{\prime}})}{|\mathbf{v}_{\mathrm{F}}(\nu,\mathbf{k^{\prime}})|}\right]\Delta_{i}^{s/t}(\nu,\mathbf{k^{\prime}})$

can then be solved in the singlet and triplet channel separately. The
integration is carried out over $\mathbf{k^{\prime}}$ where the bands $\nu$
cross the Fermi level. $V_{\mathrm{BZ}}$ is the Brillouin zone volume and
$v_{\mathrm{F}}(\nu,\mathbf{k^{\prime}})$ the Fermi velocity at
$\mathbf{k^{\prime}}$ in band $\nu$. The largest eigenvalue $\lambda_{i}$
corresponds to the leading instability $i$ with superconducting gap
$\Delta_{i}(\mu,\mathbf{k})$.

We solve the linearized gap equation by taking into account the Fermi surface
cut through the first and second Brillouin zone to include scatterings between
$k_{z}=0$ and $k_{z}=1/2$ sheets. The above integral is thus adopted to two
dimensions. The Fermi surface is parameterized by 513 points.

The leading instabilities for $U=1.5$ eV and $J=U/7$ are shown in Fig. 8,
revealing the two dominating instabilities being $d$-wave with
$d_{x^{2}-y^{2}}$-symmetry for the leading and $d_{xy}$-symmetry for the
subleading solution. In both cases the gapping of the $\alpha$\- and $\beta$\-
band is lower than for the $\gamma$-band. This is true across all coupling
parameters $U$, $J$ that exhibit a $d$-wave solution. Note, how the next
subdominant instability is in the triplet channel and thus trivially
degenerate, i.e. Fig. 8c and Fig. 8d appear rotated by $\pi/2$ relative to
each other.

Figure 7: RPA spin susceptibility Re$(\chi_{S}(\mathbf{q},\omega=0))$
calculated at $T=80$ K for $U=1.5$ eV and $J=U/7$. (a) Total contribution as
in the main text figure 4(a). (b)-(d) Different orbital resolved contributions
to Re$(\chi_{S}(\mathbf{q},\omega=0))$. Colorbars are in units of 1/eV. Figure
8: Solutions to the linearized gap equation for $k_{z}=0$ cuts at $U=1.5$ eV
and $J=U/7$. Leading instability with $d_{x^{2}-y^{2}}$-symmetry as in Fig.4b
of the main text (a) and subleading instability with $d_{xy}$-symmetry (b).
The second subleading instability (c,d) is degenerate $p$-wave in the triplet
channel.

##  IV References

[1] V. I. Anisimov, I. V. Solovyev, M. A. Korotin, M. T. Czy$\dot{\rm z}$yk,
and G. A. Sawatzky, Phys. Rev. B 48, 16929 (1993).
[2] M. Jarrell and J. E. Gubernatis, Physics Reports 269, 133 (1996).
[3] H. J. Vidberg and J. W. Serene, J Low Temp Phys 29, 179 (1977).
[4] S. Graser, T. A. Maier, P. J. Hirschfeld, and D. J. Scalapino, New Journal
of Physics 11, 025016 (2009).
[5] Z. Luo, X. Hu, M. Wang, W. Wu, and D.-X. Yao, arXiv preprint
arXiv:2305.15564 (2023).