# Two types of superconducting pairs in stripe-ordered La2-xBaxCuO4 ($x=1/8$):
evidence from the resistivity measurements

Tianhao Ren<EMAIL_ADDRESS>Condensed Matter Physics and Materials Science
Division, Brookhaven National Laboratory, Upton, New York 11973, USA  Pedro M.
Lozano  Condensed Matter Physics and Materials Science Division, Brookhaven
National Laboratory, Upton, New York 11973, USA  Department of Physics and
Astronomy, Stony Brook University, Stony Brook, New York 11794-3800, USA
Qiang Li  Condensed Matter Physics and Materials Science Division, Brookhaven
National Laboratory, Upton, New York 11973, USA  Department of Physics and
Astronomy, Stony Brook University, Stony Brook, New York 11794-3800, USA
Genda Gu  Condensed Matter Physics and Materials Science Division, Brookhaven
National Laboratory, Upton, New York 11973, USA  Alexei M. Tsvelik  Condensed
Matter Physics and Materials Science Division, Brookhaven National Laboratory,
Upton, New York 11973, USA

###### Abstract

Recent angle-resolved $c$-axis resistivity measurements of the stripe-ordered
La2-xBaxCuO4 (LBCO) with $x=1/8$ revealed an unexpected dependence on the
direction of the in-plane magnetic field. We argue that these and other
available data for the $c$-axis transport point to the existence of
superconducting pairs of two different types in the $x=1/8$ LBCO below the
stripe ordering temperature. The pairs of one type carry finite momentum and
are confined to the Cu-O planes; the pairs of other type (probably the
conventional $d$-wave with zero momentum) propagate along narrow conducting
channels traversing the sample in the $c$-axis direction. The evidence for
this comes from the observed exponential temperature dependence of the
$c$-axis resistivity $\rho_{c}(T)$ which we attribute to the thermally excited
slips of the superconducting phase and flux flows. We present a simple theory
to fit the observed $\pi/2$-periodic dependence of $\rho_{c}$ on the direction
of the in-plane magnetic field and the other data.

##  I Introduction

“It is a riddle, wrapped in a mystery, inside an enigma; but perhaps there is
a key.” These words of Winston Churchill he uttered in an entirely different
context can be used to characterize the situation with the cuprates where
enigma indeed contains in itself many layers. These structurally simple
materials display a breathtaking complexity of properties, ranging from two-
dimensional superconductivity Li et al. (2007a); Tranquada et al. (2008); Berg
et al. (2007) to magnetically driven anomalous resistive states Li et al.
(2019); Shi et al. (2020a, b). In this paper we discuss the stripe-ordered
La2-xBaxCuO4 (LBCO) with $x=1/8$, with additional insights gained from recent
angle-resolved $c$-axis resistivity measurements in the presence of an in-
plane magnetic field Lozano et al. (2022).

The most remarkable property of the stripe-ordered LBCO, discovered back in
2007, is the two-dimensional superconductivity Li et al. (2007a); Tranquada et
al. (2008). The experimental measurements were performed on the charge ordered
low temperature tetragonal (LTT) phase in 1/8 doped LBCO. This phase has a
crystal structure with two Cu-O layers per unit cell at $T\leq T_{\text{co}}$
= 54 K, where $T_{\text{co}}$ is the charge ordering temperature. In
alternative layers, the charge stripes run along the $x$ or $y$ axis
orthogonally. In the next-nearest neighboring layer, the parallel stripes are
shifted by a half period to minimize the Coulomb interaction, further doubling
the number of layers per unit cell. At temperatures below the spin ordering
temperature $T_{\text{so}}$ = 42 K, the spins between each charge stripe have
antiferromagnetic order along the stripe direction. The charge stripes act as
anti-phase boundaries for neighboring spin stripes, further doubling the unit
cell size in the direction perpendicular to the stripes. Below
$T_{\text{onset}}<T_{\text{so}}$ the in-plane resistivity experiences a sharp
drop and the system experiences a crossover into the regime of strong 2D
superconducting fluctuations leading eventually to Berezinskii-Kosterlitz-
Thouless (BKT) transition at $T_{c}^{2D}\sim 16$ K. Existence of such a sharp
BKT transition is consistent with the theoretical expectations for layered
superconductors without interlayer Josephson coupling Blatter et al. (1994).
In mean field theories $T_{\text{onset}}$ would be a transition temperature,
but in reality it marks the entrance to the regime of phase incoherence
dominated by thermally excited vortices. Also, $T_{\text{so}}$ does not depend
on magnetic field, but $T_{\text{onset}}$ does, which strongly suggests the
existence of preformed pairs presumably located on the stripes, protected by
the spin gap Berg et al. (2007); Soto et al. (2015); Šopík et al. (2015). The
$c$-axis resistivity vanishes below $T_{3D}\sim 10$ K, while the 3D
superconductivity, characterized by the Meissner state, is achieved only below
$T_{c}^{3D}\sim 4$ K. The relation between different temperature scales are
shown in Fig. 1.

Figure 1: The relation between different temperature scales in the $x=1/8$
doped La2-xBaxCuO4 (LBCO), together with the transport regimes. Figure is
adapted from Ref. (3).

The theoretical explanation suggested soon after the original discovery Berg
et al. (2007) portrayed the superconducting state as a quasi-condensate of
pairs with a finite momentum, termed as pair density wave (PDW). Then its two
dimensional character originates from the fact that momenta of the pairs in
adjacent layers have different orientations which frustrates the interlayer
Josephson coupling. This conclusion is supported by the optical measurements
which established the decoupling of the superconducting orders in adjacent
copper oxide layers Agterberg and Tsunetsugu (2008). The theory describes the
stripe-ordered state of the doped antiferromagnet as an array of self-
assembled rivers of charges separated by insulating regions. Each charge river
is represented as a Luther-Emery liquid where the formation of the spin gap
leads to the superconducting pairing; the pairing fields from different rivers
are weakly coupled by the Josephson pair tunneling. To explain the finite
momentum of the pairs, the theory in Berg et al. (2007) posits the negative
sign of the Josephson matrix element - an assumption for which so far little
theoretical justification has been provided.

It was later found in Li et al. (2019); Shi et al. (2020a, b) that when the
magnetic field destroys superconductivity, a very peculiar resistive metallic
state emerges which persists down to lowest temperatures and inherits such
prominent feature of the superconductor as zero Hall response. The critical
magnetic field is relatively small; the electrical resistance gradually
increases with the magnetic field, and at field strength around $25-30$ T, the
sheet resistance reaches a plateau at $R_{\square}\approx 2\pi\hbar/2e^{2}$.
This is dubbed “ultra-quantum metal”, which lies between $T_{\text{onset}}$
and $T_{\text{so}}$, and the Josephson coupling between the stripes is
suppressed and the superconducting fluctuations are one dimensional. There is
a dichotomy in the theoretical explanation of this state between the
quasiparticle picture Tsvelik (2016, 2019); Lee (2014) and the Bose metal
picture Das and Doniach (1999); Phillips (2018); Ren and Tsvelik (2020), where
the transport is carried by quasiparticles in the former and by the incoherent
pairs in the latter.

The previous experimental and theoretical efforts have concentrated mostly on
understanding of physics of the $ab$ planes, while the $c$-axis transport and
the possible mechanism for the Meissner state have not attracted much
attention. This situation is changed by recent angle-resolved $c$-axis
resistivity measurements Lozano et al. (2022), which attract attention to the
puzzles related to the dramatic difference between the $c$-axis transport and
the $ab$-plane transport. The experimental findings strongly suggest the
existence of two subsystems in the stripe-ordered LBCO, one being the 2D
superconductivity in the form of PDW within the $ab$ planes, the other being
the 1D superconductivity traversing along the $c$ axis. This is what we will
describe in detail in the following sections.

Before proceeding, we want to make a clarification of the relation between the
current paper and Ref. (7) by the same group of authors. On the experimental
side, the data shown later in Fig. 2, 3, and 4 are taken from Ref. (7), and
analyzed further in the current paper. The data shown later in Fig. 5, 8, 9,
and 10 are exclusively reported in the current paper, although similar data as
that of Fig. 5 and 8 for different temperatures and magnetic fields are
already reported in Ref. (7). On the theoretical side, the existence of the
two types of superconducting pairs based on the observed angle-dependent
magneto-resistivity is briefly mentioned in Ref. (7) with a promise to discuss
it in detail in subsequent publications. The current paper serves as a follow-
up to Ref. (7) fulfilling this promise, providing theoretical details and
discussing further important implications.

##  II Experimental Observations

The experimental measurements are performed on the stripe-ordered LBCO, where
the in-plane resistivity $\rho_{ab}$ and the $c$-axis resistivity $\rho_{c}$
are measured in presence of an in-plane magnetic field. The measurements are
performed at different temperatures and different strengths of the magnetic
field, with the direction of the magnetic field varied within the $ab$ plane.
Before moving on to the behavior of the $c$-axis resistivity, we briefly
address the behavior of the in-plane resistivity $\rho_{ab}$. Due to the layer
decoupling, we can focus on a single layer for the leading-order behavior of
$\rho_{ab}$. Below $T_{\text{onset}}$, the preformed pairs are protected by
the spin gap. Above $T^{2D}_{c}$, there is a proliferation of vortices, which
are thermally activated to produce a resistivity as Halperin and Nelson
(1979); Minnhagen (1987)

$\rho_{ab}\propto\rho^{n}_{ab}\exp\left(-b\sqrt{\frac{T_{\text{onset}}-T}{T-T^{2D}_{c}}}\right),$
(2.1)

where $b$ is a dimensionless fitting parameter, and the normal resistivity
$\rho^{n}_{ab}$ receives a contribution from the freed quasiparticles due to
loss of pair coherence within the vortices. The temperature dependence of Eq.
(2.1) agrees qualitatively well with the experimental data Berg et al. (2007).
The magnetic field dependence of Eq. (2.1) hides itself in the normal
resistivity $\rho^{n}_{ab}$:

$\rho^{n}_{ab}=\frac{1+\omega^{2}_{c}\tau^{2}}{\sigma_{D}},$ (2.2)

where $\sigma_{D}$ is the Drude conductivity of the quasiparticles, $\tau$ is
the elastic scattering relaxation time, and $\omega_{c}\propto B_{\perp}$ is
the cyclotron frequency, with $B_{\perp}=B\cos\phi$ being the component of the
in-plane magnetic field perpendicular to the current. The longitudinal
component $B_{\parallel}=B\sin\phi$ exerts no force on the current, as can be
seen from the Lorentz formula $\bm{F}_{L}=\bm{j}\times\bm{B}$. This renders
$\rho_{ab}$ dependent on the angle between the in-plane magnetic field and the
current (but not on the stripe directions), with a period of $\pi$. This is
consistent with the experimental data, shown in Fig. 2.

Figure 2: Characteristic measurements of the angle dependent in-plane
resistivity $\rho_{ab}$ at $B=14$T, $T=30$K. $\rho_{ab}$ depends only on the
angle between the in-plane magnetic field and the current direction, but not
on the stripe directions. The angle dependence has a period of $\pi$. These
data are originally reported in Ref. (7).

Now we proceed to the behavior of the $c$-axis resistivity. As mentioned
above, in a zero magnetic field, the $c$-axis transport is resistive above
$T_{3D}\sim 10$ K. As shown in Fig. 3, $\rho_{c}(T)$ displays a characteristic
maximum around the stripe ordering temperature and a thermally activated form
$\rho_{c}\sim\exp(-\Delta/T)$ below $T_{\text{onset}}$. This brings to mind
the Ambegaokar-Halperin formula Ambegaokar and Halperin (1969); Anderson
(1966) describing the resistivity in thin superconducting wires, where the
dissipation is generated by thermally activated phase slips:

$\rho^{\text{slip}}_{c}\propto\rho^{n}_{c}\frac{\Omega}{T}\exp\left(-\Delta^{\text{slip}}/T\right),$
(2.3)

where $\rho^{n}_{c}$ is the normal resistivity due to the interlayer
quasiparticle hopping, $\Omega$ is a characteristic frequency, and the
activation energy $\Delta^{\text{slip}}$ can be estimated as $\Delta\sim
H_{c1}d^{2}\xi_{c}$ according to the Langer-Ambegaokar-McCumber-Halperin
theory Langer and Ambegaokar (1967); McCumber and Halperin (1970), provided
the thickness $d$ of the wire is much smaller than the $c$-axis coherence
length $\xi_{c}$, and $H_{c1}$ is the bulk lower critical magnetic field. The
temperature dependence of Eq. (2.3) agrees qualitatively well with our
experimental data shown in Fig. 3 and 4, where the same data of Fig. 3 is
plotted on semi-logarithmic scale in Fig. 4 for a better visualization of the
exponential dependence in Eq. (2.3).

Figure 3: The temperature dependence of the $c$-axis resistivity $\rho_{c}$ at
$B=0$T and $B=1$T, where the in-plane magnetic field is along the Cu-Cu bond.
$\rho_{c}(T)$ displays a maximum around the stripe ordering temperature, and a
thermally activated form below $T_{\text{onset}}$. The temperature scales are
also shown, where $T^{3D}_{c}$ (below 5K) and $T_{\text{co}}$ (above 45K) are
not shown. These data are originally reported in Ref. (7). Figure 4: The
temperature dependence of the $c$-axis resistivity at $B=0$T and $B=1$T (the
in-plane magnetic field is along the Cu-Cu bond) shown on semi-logarithmic
scale as $\ln\rho_{c}(T)$ vs $1/T$. The error bars are hidden inside the data
points and can be seen if one enlarges the figure. These data are originally
reported in Ref. (7).

The recent experiments measured $\rho_{c}$ in the presence of an in-plane
magnetic field Lozano et al. (2022). The measurements were performed both
above and below $T_{\text{so}}$ as a function of the angle $\theta$ between an
in-plane magnetic field and the stripe direction. Above $T_{\text{so}}$,
$\rho_{c}$ shows no angular dependence. Below $T_{\text{so}}$, $\rho_{c}$
exhibits an angular dependence with a period of $\pi/2$. The angular
dependence does not experience any drastic change at $T_{\text{onset}}$. The
$\rho_{c}$ minima appear at $0,\pi/2,\pi,3\pi/2$, where the magnetic field is
along the directions of the stripes. The $\rho_{c}$ maxima appear at
$\pi/4,3\pi/4,5\pi/4,7\pi/4$, where the magnetic field bisects the stripe
directions in adjacent layers, or along the nodal directions of the $d$-wave
order parameters. Characteristic experimental results are shown in Fig. 5.

Figure 5: Characteristic measurements of the angle-resolved $c$-axis
resistivity at $B=1$T, $T=7.5$K. $\rho_{c}$ depends on the angle between the
in-plane magnetic field and the stripe direction. $a$ and $b$ are along the
Cu-O-Cu bond, so the minima of $\rho_{c}$ appear where the magnetic field is
along the directions of the stripes. The angle dependence has a period of
$\pi/2$. Materials and methods for these original data are discussed in
Appendix C.

In the presence of magnetic field, there is an extra contribution to
$\rho_{c}$ due to the thermally activated flux flows Bardeen and Stephen
(1965); Blatter et al. (1994):

$\rho^{\text{flow}}_{c}\propto\rho^{n}_{c}\frac{B}{H_{c2}}\exp\left(-\Delta^{\text{flow}}/T\right),$
(2.4)

where $B$ is the strength of the magnetic field, $H_{c2}$ is the bulk upper
critical magnetic field, and the activation energy $\Delta^{\text{flow}}$ is
generally smaller than $\Delta^{\text{slip}}$ (this can be seen from Fig. 4).
Consequently, in the presence of magnetic field $\rho_{c}$ is dominated by
$\rho^{\text{flow}}_{c}$. The temperature dependence of Eq. (2.4) agrees
qualitatively well with the experimental data shown in Fig. 3 and 4, while the
angular dependence must be hidden inside the normal resistivity $\rho^{n}_{c}$
111Both Eq. (2.3) and (2.4) have the prefactor $\rho^{n}_{c}$, which encodes
the dependence on the direction of the magnetic field. We assume that the bulk
critical fields $H_{c1}$, $H_{c2}$ and the activation energies
$\Delta^{\text{slip}}$, $\Delta^{\text{flow}}$ all have negligible dependence
on the direction of the in-plane magnetic field.. If the quasiparticles
possess anisotropy compatible with the PDW configuration, the rotation of the
PDW direction along the $c$ axis will come with a rotation of the axis of
anisotropy, such that the quasiparticle hopping along the $c$ axis will
produce the observed angular dependence of $\rho^{n}_{c}$, which will be shown
later. This naturally explains the different angular dependences observed
between $\rho_{ab}$ and $\rho_{c}$.

##  III The Phenomenological Model

The transport data described above indicate that above the 3D superconducting
transition in the stripe-ordered LBCO there are two types of carriers
responsible for the in-plane and $c$-axis transport respectively. These
different pairs inhabit subspaces of different dimensions. One subspace, let
us call it A, consists of two-dimensional arrays of superconducting stripes,
where pairs with finite momentum condense. The superconducting pairs formed in
the stripes of a given plane remain confined to this plane which determines
the BKT character of the transition. The other subspace, call it B, includes
areas where the pairing has a different symmetry (it is most likely a
conventional $d$-wave one, which is common for cuprates Norman and Pépin
(2003) ) and hence the corresponding order parameter field does not couple to
the one of the stripes. We assume that the subspace B consists of narrow
channels oriented along the $c$ axis, probably formed by columnar defects.
These pairs are unable to achieve phase coherence due to the phase slips and
flux flows. Since due to the momentum mismatch there is no pair tunneling
between A and B subspaces (otherwise it will break momentum conservation), the
transport along the $c$-axis can remain resistive even below the BKT
transition as seen experimentally. Since according to Eq. (2.3) and (2.4) the
observed $\rho_{c}$ is proportional to the resistivity of the normal state, to
explain the angular dependence of the magnetoresistivity, we have to study the
normal state resistivity $\rho^{n}_{c}$.

Figure 6: The two subsystems of the stripe-ordered LBCO, where within the $ab$
plane, we have the peculiar PDW configuration, and along the $c$ axis we have
the superconducting wire. The superconducting wire is enlarged on the
righthand side for illustration of the fluxes produced by the in-plane
magnetic field.

To estimate $\rho^{n}_{c}$ we adopt a simplified model of quasiparticle
hopping to illustrate that the desired angular dependence of $\rho^{n}_{c}$
can be attributed to quasiparticles of the Dirac fermion type with anisotropy
compatible with the peculiar PDW configuration shown in Fig. 6. At present we
do not have any microscopic model to support this claim. However, we are
forced to introduce the Dirac quasiparticles by the weight of the evidence
coming from the experimental data. As demonstrated in Appendix B, the scenario
where quasiparticles has quadratic dispersion should be excluded. The
discussion below provides further details.

In the simplified model, we consider two layers in the unit cell instead of
four 222The parallel stripes in the next-nearest neighboring layer are shifted
by a half period, doubling the number of layers per unit cell, but this does
not have an effect on the rotation of the axis of anisotropy, and it is the
stripe direction that matters for the oscillation of the magneto-resistivity,
as shown by the experimental data.. The simplest arrangement corresponds to
Dirac fermions with anisotropy compatible with the peculiar PDW configuration
as shown in Fig. 7, where the main axes of the cross sections are along the
directions of the stripes and rotated by 90 degrees in the nearest neighboring
layer. Although the $d$-wave superconducting regions do not have the stripe
order, they do still have the LTT structure, where there is an anisotropy of
the Cu-O hopping. In one direction, the Cu-O-Cu bonds are in line, while at 90
degrees there is a bend of $2\phi$ about 7 degrees (for an estimation of the
tilt angle, refer to Fig. 15(a) of [(26)]). The bend direction alternates from
one O site to the next, and rotates 90 degrees between adjacent layers. This
could be a potential mechanism responsible for the anisotropy investigated
here.

Figure 7: The quasiparticle of the Dirac fermion type with anisotropy
compatible with the peculiar PDW configuration, where black and red ellipses
correspond to the cross sections of the Dirac cones in nearest neighboring
layers.

Since the spin degree of freedom is irrelevant in the following discussion, we
focus on a model of spinless Dirac fermions. In terms of the momentum label
$k_{z}$ in the $c$ direction, the Hamiltonian can be written as

$\begin{split}&H_{\text{quasi}}=\sum_{\bm{k},k_{z}}\Psi^{\dagger}(\bm{k},k_{z})\mathcal{H}(\bm{k},k_{z})\Psi(\bm{k},k_{z}),\quad\Psi=\begin{pmatrix}\psi_{+}\\\
\psi_{-}\end{pmatrix},\\\
&\mathcal{H}(\bm{k},k_{z})=\begin{pmatrix}v_{x}k_{x}\sigma_{x}+v_{y}k_{y}\sigma_{y}&-t(1+e^{-\mbox{i}k_{z}a_{0}})\mathbb{I}\\\
-t(1+e^{+\mbox{i}k_{z}a_{0}})\mathbb{I}&v_{y}k_{x}\sigma_{x}+v_{x}k_{y}\sigma_{y}\end{pmatrix},\end{split}$
(3.1)

where $\bm{k}=(k_{x},k_{y})$ is the in-plane momentum, $\pm$ labels the two
layers in a unit cell and $a_{0}$ is the lattice constant in the $c$
direction. The anisotropy shown in Fig. 7 is characterized by
$\gamma^{2}\equiv v_{x}/v_{y}\neq 1$. Now we apply an in-plane magnetic field
in the $ab$ plane:

$\bm{B}=(B\cos\theta,B\sin\theta,0)~{}~{}\Rightarrow~{}~{}\bm{A}=Bz(\sin\theta,-\cos\theta,0),$
(3.2)

where $\theta$ is the angle between the in-plane magnetic field and the stripe
direction. It does not matter with respect to which layer in the unit cell
$\theta$ is defined, since the system is invariant under a rotation by 90
degrees. The in-plane magnetic field appears in the Hamiltonian through

$(k_{x},k_{y})~{}~{}\to~{}~{}(k_{x}-q_{e}zB\sin\theta,k_{y}+q_{e}zB\cos\theta),$
(3.3)

where $q_{e}$ is the electron charge and $z=-\mbox{i}\mbox{d}_{k_{z}}$ acts a
differential operator. The $c$-axis resistivity can be calculated via the Kubo
formula:

$\begin{split}&\frac{1}{\rho_{c}^{n}}=-\frac{q_{e}^{2}t^{2}}{2\pi}\int\frac{\mbox{d}^{2}k~{}\mbox{d}k_{z}\mbox{d}k^{\prime}_{z}}{(2\pi)^{4}}\int\frac{\mbox{d}y}{2T}~{}\text{sech}^{2}\left(\frac{y}{2T}\right)\\\
&\times\operatorname{Tr}\Big{[}\left(e^{-\mbox{i}k_{z}a_{0}+\mbox{i}k^{\prime}_{z}a_{0}}\right)\\\
&\quad\quad\quad\times\operatorname{Im}G^{R}_{-}(y;\bm{k};k_{z},k^{\prime}_{z})\operatorname{Im}G^{R}_{+}(y;\bm{k};k^{\prime}_{z},k_{z})\\\
&\quad\quad+\left(e^{\mbox{i}k_{z}a_{0}-\mbox{i}k^{\prime}_{z}a_{0}}\right)\\\
&\quad\quad\quad\times\operatorname{Im}G^{R}_{+}(y;\bm{k};k_{z},k^{\prime}_{z})\operatorname{Im}G^{R}_{-}(y;\bm{k};k^{\prime}_{z},k_{z})\Big{]},\end{split}$
(3.4)

where the Green’s functions are defined within each layer in a unit cell. For
a meaningful calculation, we introduce point disorders within each layer,
which can be of two types for the Dirac fermions:

$\mathcal{H}_{\text{dis}}(\bm{r})=\begin{pmatrix}V_{1}(\bm{r})+V_{2}(\bm{r})\sigma_{z}&0\\\
0&V_{1}(\bm{r})+V_{2}(\bm{r})\sigma_{z}\end{pmatrix},$ (3.5)

where the point disorders are short-range correlated:

$\mkern 1.5mu\overline{\mkern-1.5muV_{i}(\bm{r})\mkern-1.5mu}\mkern
1.5mu=0,~{}\mkern
1.5mu\overline{\mkern-1.5muV_{i}(\bm{r})V_{j}(\bm{r}^{\prime})\mkern-1.5mu}\mkern
1.5mu=\Delta_{i}\delta_{ij}\delta^{(2)}(\bm{r}-\bm{r}^{\prime}),~{}i,j=1,2.$
(3.6)

These point disorders can be caused by point defects and vortex cores within
the Cu-O plane. The effect of the point disorders can be collectively
characterized by a parameter
$\beta\equiv(\Delta_{1}+\Delta_{2})/(4v_{x}v_{y})$, and we are assuming weak
disorder $\beta\ll 1$. To the leading order in $t^{2}$, and under the
condition of large anisotropy ($\gamma^{2}\ll 1$ or $\gamma^{2}\gg 1$), the
$c$-axis resistivity $\rho^{n}_{c}(\theta)$ can be approximated by the
following expression:

$\begin{split}&\rho^{n}_{c}(\theta)=\frac{\rho_{c,0}}{g(\theta)+g(\pi/2-\theta)},\\\
&g(\theta)=\sqrt{\frac{\frac{1}{\gamma^{2}}\cos^{2}\theta+\gamma^{2}\sin^{2}\theta}{\frac{1}{\gamma^{2}}\sin^{2}\theta+\gamma^{2}\cos^{2}\theta}}\\\
&\quad\quad\quad\times
f\left(\xi\sqrt{\frac{1}{\gamma^{2}}\cos^{2}\theta+\gamma^{2}\sin^{2}\theta}\right),\\\
&f(x)=\int\mbox{d}y~{}\text{sech}^{2}y~{}\frac{y^{2}}{y^{2}+x^{2}},\end{split}$
(3.7)

where $\rho_{c,0}\propto\beta t^{-2}$ does not depend on $\theta$, and the
parameters $\gamma,\xi$ are defined as

$\gamma^{2}=v_{x}/v_{y},\quad\xi=q_{e}Ba_{0}\sqrt{v_{x}v_{y}}/(2T).$ (3.8)

The result respects the symmetry under $\gamma\to 1/\gamma$, which is required
since we have the freedom to name the directions. The fitting of Fig. 5
according to Eq. (3.7) is shown in Fig. 8.

Figure 8: The fitting of the $c$-axis resistivity according to Eq. (3.7), in
correspondence to Fig. 5. The blue solid line represents the theoretical
prediction, while the orange solid dots represent the experimental data.
Materials and methods for these original data are discussed in Appendix C.

It is clear that the $c$-axis resistivity varies periodically with the
direction of the in-plane magnetic field. The period of the resistivity
oscillation is $\pi/2$, and the minima of the resistivity lie at $k\pi/2$
while the maxima lie at $(2k+1)\pi/4$, where $k=0,1,2\cdots$. These features
agree well with the experimental data. In plotting Fig. 8, we assumed an
average velocity $v=\sqrt{v_{x}v_{y}}$ equal to 0.1% the speed of light,
roughly the velocity of the Dirac fermions in graphene. Accordingly, to
produce the observed amplitude of oscillation shown in Fig. 5, we obtain a
relatively large anisotropy $\gamma^{2}\approx 0.18$, or equivalently
$\gamma^{2}\approx 5.7$. In contrast, as discussed in Appendix B, the
conventional fermions with quadratic yet anisotropic dispersion cannot produce
the observed amplitude of oscillation shown in Fig. 5 with a realistic fermion
mass. In this sense, the experimental data favor the possibility of Dirac
fermions. By combining Eq. (3.7) together with Eq. (2.4), we can obtain the
dependence of the $c$-axis resistivity on the magnitude of the in-plane
magnetic field for a fixed temperature. We show in Fig. 9 and 10 the fitting
of our result to the experimental data, which utilizes the anisotropy
parameter $\gamma=0.42$ obtained previously. The consistency is reasonably
good given our simple model of anisotropy and approximate functional form for
large anisotropy.

Figure 9: The dependence of the $c$-axis resistivity on the magnitude of the
in-plane magnetic field, shown for a fixed direction of the magnetic field (in
fact, the direction of resistivity maximum). The blue solid line is the
theoretical prediction from Eq. (3.7) and (2.4), and the orange solid dots are
from the experimental data. Materials and methods for these original data are
discussed in Appendix C. Figure 10: The dependence of amplitude of the
$c$-axis resistivity oscillation on the magnitude of the in-plane magnetic
field. The blue solid line is the theoretical prediction from Eq. (3.7) and
(2.4), and the orange solid dots are from the experimental data. Materials and
methods for these original data are discussed in Appendix C.

The angular dependence of the $c$-axis resistivity in the simplified model
described above can be understood as follows. In order to conduct along $c$
axis, the Dirac fermions have to hop between layers. From one layer to
another, the probability of such hopping is proportional to the amount of
overlapping between the two cross sections of the Dirac cones determined by
the temperature, as shown in Fig. 7. The in-plane magnetic field has the
effect to shift the cross sections in different layers by different amount
along the direction perpendicular to the in-plane magnetic field, thus
changing their overlap, as shown in Fig. 11. If the Dirac cone in each layer
is isotropic, such shifting will also be isotropic, then the overlapping
between the cross sections and consequently the hopping is also isotropic,
resulting in an angle-independent $c$-axis resistivity. Otherwise, if there is
anisotropy in the Dirac cone, the amount of overlapping between the cross
sections changes with the shifting direction, resulting in an angle-dependent
hopping and $c$-axis resistivity. The above argument is illustrated
schematically in Fig. 11. The periodicity in the angular dependence is fixed
by the rotational symmetry - in each layer, the stripe configuration gives us
a $C_{2}$ symmetry, while in each unit cell, the orthogonal stripe directions
of the neighboring layers enhance the symmetry to $C_{4}$. As a result, we
will have a four-fold periodicity, namely the $\pi/2$ period in the $c$-axis
resistivity. As a bonus point, the anisotropy in the Dirac cone as well as the
hopping along $c$ axis relaxes the quasi-2D nature of the system inherited
from the layer structure, to the extent that the weak localization effect is
negligible, leaving us a finite $c$-axis resistivity.

Figure 11: The schematic illustration of the mechanism of the oscillation
observed in the $c$-axis resistivity in presence of an in-plane magnetic
field. The tips of the Dirac cones in nearest neighboring layers are shifted
by the in-plane magnetic field. The magnitude of the shift is fixed by the
strength of the magnetic field and the lattice constant, while the direction
of the shift varies with the direction of the magnetic field. If there is
anisotropy, such a shift will produce a varying overlap between the cross
sections, and thus a varying $c$-axis resistivity.

##  IV Discussion

In the current paper, we have argued that the recent angle-resolved $c$-axis
resistivity measurements of the stripe-ordered LBCO, together with temperature
dependence, imply the existence of two weakly coupled subsystems containing
superconducting pairs of different symmetry. Pairs of the first type are
located on charge stripes and are confined to $ab$ planes with the peculiar
PDW configuration first proposed in Berg et al. (2007), such that the
interlayer Josephson coupling is nullified. The pairs of the other type must
have a different symmetry, they are most likely conventional $d$-wave pairs,
whose presence is due to doping inhomogeneities. The one-dimensional character
of the $c$-axis transport indicates that these pairs must traverse the sample
along narrow channels along the $c$ direction without interruption, likely
formed by columnar defects. The layer decoupling can be destroyed by defects
in the charge-stripe order. Then, with decreasing temperature, the two
subsystems merge into one and first go through a glassy transition to a
dissipationless state, and then evolve into a Meissner state.

The phenomenological model of the two-subsystem picture consists of preformed
pairs moving in the charge stripes, and pairing of quasiparticles confined in
the columnar defects. It can provide a good, qualitative explanation of the
temperature dependence of the experimental data. In the mean time, the
observed dependence of the $c$-axis resistivity on the direction of the in-
plane magnetic field can be explained by introducing anisotropy of the normal
state quasiparticles compatible with the peculiar in-plane stripe order, and
the experimental data favor the possibility of the quasiparticles being Dirac
fermions (which we will discuss further at the end of this section). As a
consequence, when we apply a magnetic field strong enough to destroy the
coherence between the preformed pairs, and free the Dirac fermions to
propagate in three dimensions, we will obtain a resistive state in presence of
weak disorder. Such a resistive state has a zero Hall response due to the
particle-hole symmetry of the Dirac fermions, which is consistent with the
experimental findings Li et al. (2019). This gives support to the
quasiparticle picture over the Bose metal picture about the nature of such a
resistive state.

The Dirac fermion picture favored by the experimental data as proposed in this
paper is unusual, so we emphasize the logic behind our reasoning here. To
explain the magneto-resistivity, we have two options at hand, one is a
quasiparticle mechanism, and the other is a Bose metal. We can rule out the
later option since the incoherent pairs would not react to the direction of
the magnetic field. As a consequence there would be no oscillation of the
magneto-resistivity as seen in Fig. 5. As far as the quasiparticle mechanism
is concerned, the first candidate is the conventional quasiparticles with
quadratic dispersion. However, in Appendix B we have shown that in order to
explain the experimental data, the quasiparticle mass must be vanishingly
small. This brings us to the next candidate on the list, which is the Dirac
fermions. We have shown in the main text that the Dirac fermions produce a
surprising agreement with the experimental data (the oscillation pattern with
respect to the direction of the in-plane magnetic field, and the oscillation
amplitude with respect to the strength of the in-plane magnetic field), given
the simplicity of our phenomenology. Therefore, our result provides a minimal
phenomenological model that is consistent with the experimental data. Although
we are forced to introduce the Dirac quasiparticles by the weight of the
experimental data, we admit that at present we cannot suggest any microscopic
justification for this hypothesis. To date, Dirac fermions are found in
topological insulators and semimetals like graphene. For $d$-wave
superconducting cuprates, spinons with linear dispersion are proposed for the
undoped system Affleck and Marston (1988), and Bogoliubov quasiparticles
(effectively of spin $1/2$ and charge $0$) can present a Dirac fermion type
dispersion along the nodal direction of the Brillouin zone Vafek and
Vishwanath (2014), but the normal state quasiparticles (of spin $1/2$ and
charge $e$) with Dirac fermion type dispersion still evade direct or indirect
experimental observations. To our knowledge, the only exception is the
observation of a closely related type-II Dirac fermion in cuprates with a
different chemical substitution, La1.77Sr0.23CuO4, where the same Cu and O
bands are involved and the observed tilted linear dispersion is attributed to
the tight-binding band structure Horio et al. (2018). Therefore, we hope that
our work may encourage theoretical efforts in the direction of Dirac fermions
in the stripe-ordered LBCO.

The phenomenological theory discussed in this paper has relations with other
theories previously proposed. The theoretical description in Tsvelik (2019)
provides a possible microscopic mechanism for the preformed pairs moving on
the charge stripes and quasiparticles hopping between layers, using a strongly
correlated Kondo lattice model, but the quasiparticles are conventional
fermions with quadratic dispersion. On the other hand, Dirac fermions have
been proposed Mulligan and Raghu (2016); Mulligan (2017) to explain the
magnetically induced resistive state in disordered superconducting films, but
the Dirac fermions are of the composite type. Thus, a consistent microscopic
mechanism for the phenomenological theory discussed in this paper requires
further investigations. On the experimental side, a definitive experimental
evidence should come from the angle-resolved photoemission spectroscopy
(ARPES) identified with the normal state band structure, as has been done for
La1.77Sr0.23CuO4. The position of the Dirac point should be tunable by
changing the chemical substitution and the doping level, which is also subject
to experimental observations. In addition to the Dirac fermions, the formation
mechanism, the geometric property, and the distribution of the narrow
superconducting channels that host the Dirac fermions also call for a further
investigation. Other experimental observations such as thermal transport
measurements are need to obtain a full understanding of these detailed
properties.

##  V Acknowledgement

We are grateful to John Tranquada for his suggestion of the possible
anisotropy mechanism and Dragana Popovic for extensive discussions and
valuable suggestions. This work was supported by Office of Basic Energy
Sciences, Material Sciences and Engineering Division, U.S. Department of
Energy (DOE) under Contracts No. DE-SC0012704. This paper was written in part
while AMT was a visitor at KITP, which is supported by by the National Science
Foundation under Grant No. NSF PHY-1748958.

##  Appendix A Quasiparticles with Linear Dispersion

Here we show details of the calculation of the $c$-axis normal conductivity
$\sigma_{c}$ due to hopping of quasiparticles with linear dispersion (namely
the Dirac fermions), using the simplified model introduced in the main text.
The inverse of $\sigma_{c}$ will be the $c$-axis normal resistivity
$\rho^{n}_{c}$ discussed in the main text.

###  A.1 Green’s Function for Linear Dispersion

The Hamiltonian for the Dirac fermions is

$\begin{split}&H_{\text{quasi}}=\sum_{\bm{k},k_{z}}\Psi^{\dagger}(\bm{k},k_{z})\mathcal{H}(\bm{k},k_{z})\Psi(\bm{k},k_{z}),~{}\Psi=\begin{pmatrix}\psi_{+}\\\
\psi_{-}\end{pmatrix},\\\
&\mathcal{H}(\bm{k},k_{z})=\begin{pmatrix}v_{x}k_{x}\sigma_{x}+v_{y}k_{y}\sigma_{y}&-t(1+e^{-\mbox{i}k_{z}a_{0}})\mathbb{I}\\\
-t(1+e^{+\mbox{i}k_{z}a_{0}})\mathbb{I}&v_{y}k_{x}\sigma_{x}+v_{x}k_{y}\sigma_{y}\end{pmatrix},\end{split}$
(A.1)

where $\bm{k}=(k_{x},k_{y})$, $\pm$ labels the two layers in a unit cell, and
$a_{0}$ is the lattice constant in the $c$ direction. In presence of the in-
plane magnetic field, we need to perform the following substitution:

$\begin{split}&k_{x}\to\tilde{k}_{x}=k_{x}-q_{e}zB\sin\theta\\\
&k_{y}\to\tilde{k}_{y}=k_{y}+q_{e}zB\cos\theta,\end{split}$ (A.2)

where $z=-\mbox{i}\mbox{d}_{k_{z}}$ in the momentum representation. The
Matsubara Green’s function $\hat{G}\equiv-\left<\Psi\Psi^{\dagger}\right>$ is
determined by the following equation:

$\begin{split}\left[\mbox{i}\omega_{n}-\mathcal{H}(\bm{k};k_{z})\right]\hat{G}(\mbox{i}\omega_{n};\bm{k};k_{z},k^{\prime}_{z})=\delta(k_{z}-k^{\prime}_{z})\begin{pmatrix}\mathbb{I}&0\\\
0&\mathbb{I}\end{pmatrix}.\end{split}$ (A.3)

The matrix Green’s function $\hat{G}$ can be partitioned as

$\begin{split}&\hat{G}=\begin{pmatrix}G_{+}&F\\\
\bar{F}&G_{-}\end{pmatrix},\\\
&\begin{cases}G_{+}=-\left<\psi_{+}\psi^{\dagger}_{+}\right>,&G_{-}=-\left<\psi_{-}\psi^{\dagger}_{-}\right>\\\
F=-\left<\psi_{+}\psi^{\dagger}_{-}\right>,&\bar{F}=-\left<\psi_{-}\psi^{\dagger}_{+}\right>\end{cases}.\end{split}$
(A.4)

To evaluate the conductivity to leading order in $t$, we can just keep the
zeroth order Green’s functions:

$\begin{split}&F^{(0)}(\mbox{i}\omega_{n};\bm{k};k_{z},k^{\prime}_{z})=0,\quad\bar{F}^{(0)}(\mbox{i}\omega_{n};\bm{k};k_{z},k^{\prime}_{z})=0,\\\
&G^{(0)}_{+}(\mbox{i}\omega_{n};\bm{k};k_{z},k^{\prime}_{z})\\\
=&-\frac{1}{2\pi}\int\mbox{d}z\frac{\mbox{i}\omega_{n}+v_{x}\tilde{k}_{x}\sigma_{x}+v_{y}\tilde{k}_{y}\sigma_{y}}{\omega^{2}_{n}+m^{2}+v^{2}_{x}\tilde{k}^{2}_{x}+v^{2}_{y}\tilde{k}^{2}_{y}}e^{\mbox{i}(k_{z}-k^{\prime}_{z})z},\\\
&G^{(0)}_{-}(\mbox{i}\omega_{n};\bm{k};k_{z},k^{\prime}_{z})\\\
=&-\frac{1}{2\pi}\int\mbox{d}z\frac{\mbox{i}\omega_{n}+v_{y}\tilde{k}_{x}\sigma_{x}+v_{x}\tilde{k}_{y}\sigma_{y}}{\omega^{2}_{n}+m^{2}+v^{2}_{y}\tilde{k}^{2}_{x}+v^{2}_{x}\tilde{k}^{2}_{y}}e^{\mbox{i}(k_{z}-k^{\prime}_{z})z},\end{split}$
(A.5)

and the physical Green’s function $G^{R}_{\pm}$ can be obtained by the
analytic continuation $\mbox{i}\omega_{n}\to\omega+\mbox{i}0^{+}$. Now we
introduce the point disorders within each layer mentioned in the main text:

$\begin{split}&\mathcal{H}_{\mathrm{dis}}(\boldsymbol{r})=\left(\begin{array}[]{cc}V_{1}(\boldsymbol{r})+V_{2}(\boldsymbol{r})\sigma_{z}&0\\\
0&V_{1}(\boldsymbol{r})+V_{2}(\boldsymbol{r})\sigma_{z}\end{array}\right),\\\
&\overline{V_{i}(\boldsymbol{r})}=0,\quad\overline{V_{i}(\boldsymbol{r})V_{j}\left(\boldsymbol{r}^{\prime}\right)}=\Delta_{i}\delta_{ij}\delta^{(2)}\left(\boldsymbol{r}-\boldsymbol{r}^{\prime}\right),\end{split}$
(A.6)

where $i,j$ takes the values $1,2$. The point disorders introduce a self-
energy to the Dirac fermions as

$\begin{split}\Sigma^{R}(\omega;\bm{k};k_{z},k^{\prime}_{z})=&\Delta_{1}\int\frac{\mbox{d}^{2}k}{(2\pi)^{2}}G^{R}_{\pm}(\omega;\bm{k};k_{z},k^{\prime}_{z})\\\
+&\Delta_{2}\int\frac{\mbox{d}^{2}k}{(2\pi)^{2}}\sigma_{z}G^{R}_{\pm}(\omega;\bm{k};k_{z},k^{\prime}_{z})\sigma_{z}.\end{split}$
(A.7)

The calculation of $\Sigma^{R}$ involves the following integral, which can be
approximated as Shon and Ando (1998)

$\begin{split}&\int\frac{\mbox{d}^{2}k}{(2\pi)^{2}}\frac{\omega+\mbox{i}0^{+}+v_{x}\tilde{k}_{x}\sigma_{x}+v_{y}\tilde{k}_{y}\sigma_{y}}{(\omega+\mbox{i}0^{+})^{2}-\left(v^{2}_{x}\tilde{k}^{2}_{x}+v^{2}_{y}\tilde{k}^{2}_{y}\right)}\\\
\approx&\frac{A}{2(2\pi)^{2}}\int_{-\infty}^{\infty}\mbox{d}x\left<\frac{\omega+\mbox{i}0^{+}+v_{x}\tilde{k}_{x}\sigma_{x}+v_{y}\tilde{k}_{y}\sigma_{y}}{(\omega+\mbox{i}0^{+})^{2}-x}\right>_{\omega},\\\
=&\frac{A}{2(2\pi)^{2}}\int_{-\infty}^{\infty}\mbox{d}x\frac{\omega+\mbox{i}0^{+}}{(\omega+\mbox{i}0^{+})^{2}-x}=-i\frac{|\omega|}{4v_{x}v_{y}},\end{split}$
(A.8)

where the average is over the cross section
$v^{2}_{x}\tilde{k}^{2}_{x}+v^{2}_{y}\tilde{k}^{2}_{y}=\omega^{2}$, and
$A=\nu(|\omega|)/|\omega|=2\pi/(v_{x}v_{y})$ with $\nu(|\omega|)$ being the
density of states. In the end, the disorder averaged Green’s function can be
represented as

$\begin{split}&\mkern
1.5mu\overline{\mkern-1.5muG^{R}_{+}(\omega;\bm{k};k_{z},k^{\prime}_{z})\mkern-1.5mu}\mkern
1.5mu\\\
=&-\frac{1}{2\pi}\int\mbox{d}z\frac{\omega+\mbox{i}\Gamma+v_{x}\tilde{k}_{x}\sigma_{x}+v_{y}\tilde{k}_{y}\sigma_{y}}{-(\omega+\mbox{i}\Gamma)^{2}+v^{2}_{x}\tilde{k}^{2}_{x}+v^{2}_{y}\tilde{k}^{2}_{y}}e^{\mbox{i}(k_{z}-k^{\prime}_{z})z},\\\
&\mkern
1.5mu\overline{\mkern-1.5muG^{R}_{-}(\omega;\bm{k};k_{z},k^{\prime}_{z})\mkern-1.5mu}\mkern
1.5mu\\\
=&-\frac{1}{2\pi}\int\mbox{d}z\frac{\omega+\mbox{i}\Gamma+v_{y}\tilde{k}_{x}\sigma_{x}+v_{x}\tilde{k}_{y}\sigma_{y}}{-(\omega+\mbox{i}\Gamma)^{2}+v^{2}_{y}\tilde{k}^{2}_{x}+v^{2}_{x}\tilde{k}^{2}_{y}}e^{\mbox{i}(k_{z}-k^{\prime}_{z})z},\end{split}$
(A.9)

where $\Gamma(\omega)=|\omega|(\Delta_{1}+\Delta_{2})/(4v_{x}v_{y})$. In the
following calculations, we assume that $\Gamma(\omega)\ll|\omega|$, which
means a weak disorder $\beta\equiv(\Delta_{1}+\Delta_{2})/(4v_{x}v_{y})\ll 1$.

###  A.2 Kubo Formula for $c$-axis Conductivity

We focus on the DC conductivity obtained as setting $\bm{k}\to\bm{0}$ and
taking the limit $\omega\to 0$ in the end. The current operator with
$\bm{k}\to\bm{0}$ is

$\begin{split}&J_{z}(\bm{0},\omega)=q_{e}t\sum_{\bm{k}}\left[\Psi^{\dagger}(\omega;\bm{k},k_{z})\Sigma_{z}(k_{z})\Psi(\omega;\bm{k},k_{z})\right],\\\
&\Sigma_{z}(k_{z})=\begin{pmatrix}0&\mbox{i}e^{-\mbox{i}k_{z}a_{0}}\mathbb{I}\\\
-\mbox{i}e^{\mbox{i}k_{z}a_{0}}\mathbb{I}&0\end{pmatrix}.\end{split}$ (A.10)

Then the $c$-axis conductivity can be derived from the Kubo formula

$\operatorname{Re}\sigma_{c}(\omega)=\frac{1-e^{-\beta\omega}}{2\omega
V}\int_{-\infty}^{\infty}\mbox{d}t~{}e^{\mbox{i}\omega
t}\left<J^{\dagger}_{z}(\bm{0},t)J_{z}(\bm{0},0)\right>,$ (A.11)

where $V=V_{2D}V_{z}$ is the size of the system. In the leading order
($\sigma_{c}\propto t^{2}$), we obtain

$\sigma_{c}(\omega)=\frac{\mbox{i}}{\omega}q_{e}^{2}t^{2}\Pi_{zz}(\omega),$
(A.12)

where $\Pi_{zz}(\omega)$ is obtained from its Matsubara representation
$\Pi_{zz}(\mbox{i}\Omega_{n})$:

$\begin{split}&\Pi_{zz}(\mbox{i}\Omega_{n})\\\
=&\frac{1}{\pi^{2}V}\sum_{\bm{k}}\sum_{k_{z},k^{\prime}_{z}}\iint\mbox{d}\xi\mbox{d}\xi^{\prime}\frac{\tanh(\xi/2T)-\tanh(\xi^{\prime}/2T)}{i\Omega_{n}+\xi-\xi^{\prime}}\\\
&\times\operatorname{Tr}\Big{[}\left(e^{-\mbox{i}k_{z}a_{0}+\mbox{i}k^{\prime}_{z}a_{0}}\right)\\\
&\quad\quad\quad\times\operatorname{Im}G^{R}_{-}(\xi;\bm{k};k_{z},k^{\prime}_{z})\operatorname{Im}G^{R}_{+}(\xi^{\prime};\bm{k};k^{\prime}_{z},k_{z})\\\
&\quad\quad+\left(e^{\mbox{i}k_{z}a_{0}-\mbox{i}k^{\prime}_{z}a_{0}}\right)\\\
&\quad\quad\quad\times\operatorname{Im}G^{R}_{+}(\xi;\bm{k};k_{z},k^{\prime}_{z})\operatorname{Im}G^{R}_{-}(\xi^{\prime};\bm{k};k^{\prime}_{z},k_{z})\Big{]}.\end{split}$
(A.13)

Then the $c$-axis conductivity in the DC limit $\omega\to 0$ is

$\begin{split}\sigma_{c}=&-\frac{q_{e}^{2}t^{2}}{\pi
V}\sum_{\bm{k}}\sum_{k_{z},k^{\prime}_{z}}\int\frac{\mbox{d}y}{2T}~{}\text{sech}^{2}\left(\frac{y}{2T}\right)\\\
&\times\operatorname{Tr}\Big{[}\left(e^{-\mbox{i}k_{z}a_{0}+\mbox{i}k^{\prime}_{z}a_{0}}\right)\\\
&\quad\quad\quad\times\operatorname{Im}G^{R}_{-}(y;\bm{k};k_{z},k^{\prime}_{z})\operatorname{Im}G^{R}_{+}(y;\bm{k};k^{\prime}_{z},k_{z})\\\
&\quad\quad+\left(e^{\mbox{i}k_{z}a_{0}-\mbox{i}k^{\prime}_{z}a_{0}}\right)\\\
&\quad\quad\quad\times\operatorname{Im}G^{R}_{+}(y;\bm{k};k_{z},k^{\prime}_{z})\operatorname{Im}G^{R}_{-}(y;\bm{k};k^{\prime}_{z},k_{z})\Big{]},\end{split}$
(A.14)

and we need to perform the disorder average for the expression in Eq. (A.14).
Under the weak disorder condition $\beta\ll 1$, we can ignore the vertex
correction and replace the individual Green’s function $G^{R}_{\pm}$ with its
disorder average $\mkern
1.5mu\overline{\mkern-1.5muG^{R}_{\pm}\mkern-1.5mu}\mkern 1.5mu$ determined in
Eq. (A.9).

###  A.3 Calculation of the $c$-axis Conductivity

We consider the evaluation of Eq. (A.14) with $G^{R}_{\pm}$ replaced with
$\mkern 1.5mu\overline{\mkern-1.5muG^{R}_{\pm}\mkern-1.5mu}\mkern 1.5mu$. The
integration over $k_{z},k_{z}^{\prime}$ can be carried out straightforwardly,
while the integration over $\bm{k}$ can be carried out using the strategy
shown in Eq. (A.8). Finally, we obtain

$\sigma_{c}=\frac{q_{e}^{2}t^{2}}{2\pi
v_{x}v_{y}}\int\frac{\mbox{d}y}{2T}~{}\text{sech}^{2}\left(\frac{y}{2T}\right)\left[I(y,\theta)+I^{*}(y,\theta)\right]$
(A.15)

where the function $I(y)$ has the following expression:

$\begin{split}I(y>0,\theta)=&\left(\left<\frac{y^{2}+\Gamma^{2}+v_{x}v_{y}\left(k_{x}k^{\prime}_{x}+k_{y}k^{\prime}_{y}\right)}{4y\Gamma-
iF_{+}}\right>_{+}\right.\\\
+&\left.\left<\frac{y^{2}+\Gamma^{2}+v_{x}v_{y}\left(k_{x}k^{\prime}_{x}+k_{y}k^{\prime}_{y}\right)}{4y\Gamma+iF_{-}}\right>_{-}\right),\\\
I(y<0,\theta)=&\left(\left<\frac{y^{2}+\Gamma^{2}+v_{x}v_{y}\left(k_{x}k^{\prime}_{x}+k_{y}k^{\prime}_{y}\right)}{4|y|\Gamma+iF_{+}}\right>_{+}\right.\\\
+&\left.\left<\frac{y^{2}+\Gamma^{2}+v_{x}v_{y}\left(k_{x}k^{\prime}_{x}+k_{y}k^{\prime}_{y}\right)}{4|y|\Gamma-
iF_{-}}\right>_{-}\right).\end{split}$ (A.16)

The subscript $\pm$ denotes the two layers within a unit cell, and the average
is parameterized by $\alpha\in[0,2\pi)$ according to the following rules:

$\begin{split}+:&\quad k_{x}=\frac{|y|\cos\alpha}{v_{x}},\quad
k_{y}=\frac{|y|\sin\alpha}{v_{y}};\quad
k^{\prime}_{x}=-q_{e}Ba_{0}\sin\theta+\frac{|y|\cos\alpha}{v_{x}},\quad
k^{\prime}_{y}=q_{e}Ba_{0}\cos\theta+\frac{|y|\sin\alpha}{v_{y}},\\\ &\quad
F_{+}=v^{2}_{y}\left(-q_{e}Ba_{0}\sin\theta+|y|\cos\alpha/v_{x}\right)^{2}+v^{2}_{x}\left(q_{e}Ba_{0}\cos\theta+|y|\sin\alpha/v_{y}\right)^{2}-y^{2},\\\
-:&\quad k_{x}=q_{e}Ba_{0}\sin\theta+\frac{|y|\cos\alpha}{v_{y}},\quad
k_{y}=-q_{e}Ba_{0}\cos\theta+\frac{|y|\sin\alpha}{v_{x}};\quad
k^{\prime}_{x}=\frac{|y|\cos\alpha}{v_{y}},\quad
k^{\prime}_{y}=\frac{|y|\sin\alpha}{v_{x}},\\\ &\quad
F_{-}=v^{2}_{x}\left(q_{e}Ba_{0}\sin\theta+|y|\cos\alpha/v_{y}\right)^{2}+v^{2}_{y}\left(-q_{e}Ba_{0}\cos\theta+|y|\sin\alpha/v_{x}\right)^{2}-y^{2}.\end{split}$
(A.17)

Introducing the following parameters:

$\begin{split}&\epsilon_{0}\equiv vq_{e}Ba_{0},\quad v_{x}=v\gamma,\quad
v_{y}=v/\gamma,\\\ &\beta=(\Delta_{1}+\Delta_{2})/(4v^{2}),\end{split}$ (A.18)

we can express $I(y,\theta)+I^{*}(y,\theta)=J(y,\theta)+J(y,\pi/2-\theta)$ as

$\begin{split}&J(y,\theta)=\left<\frac{8\beta
y^{2}\left[(1+\beta^{2})y^{2}+P\right]}{(4\beta
y^{2})^{2}+Q^{2}}\right>_{\alpha\in[0,2\pi)},\\\
&P(\alpha)=\left(\frac{y^{2}\cos^{2}\alpha}{\gamma^{2}}+\gamma^{2}y^{2}\sin^{2}\alpha+\gamma\epsilon_{0}|y|\cos\theta\sin\alpha\right.\\\
&\quad\quad\quad\quad\left.-\frac{\epsilon_{0}|y|\sin\theta\cos\alpha}{\gamma}\right),\\\
&Q(\alpha)=\frac{1}{\gamma^{2}}\left(\frac{|y|\cos\alpha}{\gamma}-\epsilon_{0}\sin\theta\right)^{2}\\\
&\quad\quad\quad+\gamma^{2}\left(\gamma|y|\sin\alpha+\epsilon_{0}\cos\theta\right)^{2}-y^{2}.\\\
\end{split}$ (A.19)

An important feature of $\sigma_{c}$ expressed in the above formulas is that
it is symmetric under $\gamma\to 1/\gamma$, since the system is invariant
under rotation by 90 degrees, or in another word, we have the freedom to
define $x,y$-axis in neighboring layers. The analytical and numerical
evaluation of the average over $\alpha$ presents obstacles due to the singular
behavior of the integrand. As a practical approximation, we work in the limit
$\gamma\ll 1$ and $\beta\ll 1$, and take into account only the dominating
contributions near the singular points of the integrand. The singular points
are determined by

$Q\approx
0\quad\Rightarrow\quad\cos\alpha_{1,2}\approx\frac{\gamma\epsilon_{0}\sin\theta}{|y|},\quad\alpha_{1}+\alpha_{2}=\pi,$
(A.20)

and the function $J(y,\theta)$ is approximated as

$\begin{split}J(y,\theta)\approx&\frac{1}{2\pi}\int\mbox{d}\alpha\frac{8\beta
y^{2}\left[y^{2}+P(\alpha_{1})\right]}{(4\beta
y^{2})^{2}+[Q^{\prime}(\alpha_{1})]^{2}(\alpha-\alpha_{1})^{2}}\\\
+&\frac{1}{2\pi}\int\mbox{d}\alpha\frac{8\beta
y^{2}\left[y^{2}+P(\alpha_{2})\right]}{(4\beta
y^{2})^{2}+[Q^{\prime}(\alpha_{2})]^{2}(\alpha-\alpha_{2})^{2}}\\\
\approx&\frac{1}{\gamma^{6}}\left|\frac{\cos\theta}{\sin\theta}\right|\frac{y^{2}}{y^{2}+\epsilon^{2}_{0}\cos^{2}\theta/\gamma^{2}}.\end{split}$
(A.21)

Substituting this into Eq. (A.15), we obtain

$\begin{split}&\sigma_{c}(\theta)\approx\frac{q^{2}_{e}t^{2}}{2\pi\beta
v_{x}v_{y}}\frac{1}{\gamma^{6}}\left[K(\theta)+K(\pi/2-\theta)\right]\\\
&K(\theta)=\left|\frac{\cos\theta}{\sin\theta}\right|\int\frac{\mbox{d}y}{2T}~{}\text{sech}^{2}\left(\frac{y}{2T}\right)~{}\frac{y^{2}}{y^{2}+\epsilon_{0}^{2}\cos^{2}\theta/\gamma^{2}}.\end{split}$
(A.22)

This approximate result misses the information from the opposite limit
$\gamma\gg 1$, which can be cured by performing a symmetrization according to
the required symmetry under $\gamma\to 1/\gamma$. Finally, we obtain

$\begin{split}&\sigma_{c}(\theta)\approx\frac{q^{2}_{e}t^{2}}{2\pi\beta
v_{x}v_{y}}\left(\gamma^{2}+\frac{1}{\gamma^{2}}\right)^{3}\left[g(\theta)+g(\pi/2-\theta)\right],\\\
&g(\theta)=\sqrt{\frac{(1/\gamma^{2})\cos^{2}\theta+\gamma^{2}\sin^{2}\theta}{(1/\gamma^{2})\sin^{2}\theta+\gamma^{2}\cos^{2}\theta}}\\\
&\quad\quad\quad\times
f\left(\xi\sqrt{(1/\gamma^{2})\cos^{2}\theta+\gamma^{2}\sin\theta}\right),\\\
&f(x)=\int\mbox{d}y~{}\text{sech}^{2}y~{}\frac{y^{2}}{y^{2}+x^{2}},\end{split}$
(A.23)

where we have introduced the parameter $\xi=\epsilon_{0}/(2T)$. The fitting of
Eq. (A.23) to the experimental data is shown in the main text, and we obtain
$\gamma^{2}=v_{x}/v_{y}\approx 0.18$ for the averaged velocity equal to 0.1%
the speed of light.

##  Appendix B Quasiparticles with Quadratic Dispersion

Here we discuss the $c$-axis normal conductivity $\sigma_{c}$ due to hopping
of quasiparticles with quadratic dispersion, where the inverse of $\sigma_{c}$
will be the $c$-axis normal resistivity $\rho^{n}_{c}$. We will show the
reason why such conventional quasiparticles with quadratic dispersion fail to
explain the experimental data.

###  B.1 Green’s Function for Quadratic Dispersion

The Hamiltonian for the conventional fermions is

$\begin{split}&H_{\text{quasi}}=\sum_{\bm{k},k_{z}}\Psi^{\dagger}(\bm{k},k_{z})\mathcal{H}(\bm{k},k_{z})\Psi(\bm{k},k_{z}),\quad\Psi=\begin{pmatrix}\psi_{+}\\\
\psi_{-}\end{pmatrix},\\\
&\mathcal{H}=\begin{pmatrix}\epsilon_{+}(\bm{k})&-t(1+e^{-\mbox{i}k_{z}a_{0}})\\\
-t(1+e^{\mbox{i}k_{z}a_{0}})&\epsilon_{-}(\bm{k})\end{pmatrix},\\\
&\epsilon_{+}(k_{x},k_{y})=\epsilon_{-}(k_{y},k_{x})=\frac{k^{2}_{x}}{2m_{x}}+\frac{k^{2}_{y}}{2m_{y}}-\mu,\end{split}$
(B.1)

where in presence of the in-plane magnetic field, we need to perform the
substitution shown in Eq. (A.2). The Matsubara Green’s function
$\hat{G}\equiv-\left<\Psi\Psi^{\dagger}\right>$ is determined by

$\left[\mbox{i}\omega_{n}-\mathcal{H}(\bm{k},k_{z})\right]\hat{G}(i\omega_{n};\bm{k};k_{z},k^{\prime}_{z})=\delta(k_{z}-k^{\prime}_{z})\mathbb{I},$
(B.2)

where the matrix Green’s function again can be partitioned as

$\hat{G}=\begin{pmatrix}G_{+}&F\\\ \bar{F}&G_{-}\end{pmatrix},$ (B.3)

and to evaluate the conductivity to leading order in $t$, we can just keep the
zeroth order Green’s functions:

$\begin{split}&F^{(0)}(\mbox{i}\omega_{n};\bm{k};k_{z},k^{\prime}_{z})=0,\quad\bar{F}^{(0)}(\mbox{i}\omega_{n};\bm{k};k_{z},k^{\prime}_{z})=0,\\\
&G^{(0)}_{+}(\mbox{i}\omega_{n};\bm{k};k_{z},k^{\prime}_{z})=\frac{1}{2\pi}\int\mbox{d}z\frac{1}{\mbox{i}\omega_{n}-\epsilon_{+}(\bm{k},z)}e^{\mbox{i}(k_{z}-k^{\prime}_{z})z},\\\
&G^{(0)}_{-}(\mbox{i}\omega_{n};\bm{k};k_{z},k^{\prime}_{z})=\frac{1}{2\pi}\int\mbox{d}z\frac{1}{\mbox{i}\omega_{n}-\epsilon_{-}(\bm{k},z)}e^{\mbox{i}(k_{z}-k^{\prime}_{z})z}.\end{split}$
(B.4)

where
$\epsilon_{\pm}(\bm{k},z)\equiv\epsilon_{\pm}(k_{x}-q_{e}zB\sin\theta,k_{y}+q_{e}zB\cos\theta)$.
We then take into account the effect of the disorder within each layer:

$\mkern 1.5mu\overline{\mkern-1.5muV\mkern-1.5mu}\mkern 1.5mu=0,\quad\mkern
1.5mu\overline{\mkern-1.5muV(\bm{r})V(\bm{r}^{\prime})\mkern-1.5mu}\mkern
1.5mu=\Delta\delta^{(2)}(\bm{r}-\bm{r}^{\prime}).$ (B.5)

Consequently, the disorder-averaged physical Green’s function can be presented
as

$\mkern
1.5mu\overline{\mkern-1.5muG^{R}_{\pm}(\omega;\bm{k};k_{z},k^{\prime}_{z})\mkern-1.5mu}\mkern
1.5mu=\frac{1}{2\pi}\int\mbox{d}z\frac{e^{\mbox{i}(k_{z}-k^{\prime}_{z})z}}{\omega+\mbox{i}/(2\tau)-\epsilon_{\pm}(\bm{k},z)},$
(B.6)

and $\tau$ is the elastic scattering relaxation time, which can be taken as a
constant $\tau=1/(2\pi\nu\Delta)$ with $\nu$ being the constant density of
states in 2D, if we consider low enough temperatures.

###  B.2 Calculation of the $c$-axis Conductivity

We can still use the Kubo formula shown in Eq. (A.14), where the trace
operation becomes trivial. In case of weak disorder such that the vertex
correction can be ignored, we just replace the Green’s function $G^{R}_{\pm}$
with the disorder-averaged ones $\mkern
1.5mu\overline{\mkern-1.5muG^{R}_{\pm}\mkern-1.5mu}\mkern 1.5mu$. Similar in
spirit to the strategy shown in Eq. (A.8), here we replace the integration
over $\bm{k}$ with an integration over energy, for example:

$\begin{split}&\int\frac{\mbox{d}^{2}k}{(2\pi)^{2}}\frac{f(\bm{k})}{\omega+i/(2\tau)-\epsilon_{+}(\bm{k},z)}\\\
\approx~{}&\nu\int\mbox{d}\zeta\left<\frac{f(\bm{k})}{\omega+i/(2\tau)-\zeta}\right>_{\text{FS}_{+}},\end{split}$
(B.7)

where the average is over the Fermi surface determined by
$\epsilon_{+}(\bm{k},z)=0$. Following this strategy, in the leading order in
$t^{2}$, the $c$-axis conductivity is obtained as

$\sigma_{c}(\theta)=\left<\frac{4\nu
q^{2}_{e}t^{2}\tau}{1+\tau^{2}\mu^{2}f^{2}(\theta)}+\frac{4\nu
q^{2}_{e}t^{2}\tau}{1+\tau^{2}\mu^{2}f^{2}(\pi/2-\theta)}\right>_{\text{FS}},$
(B.8)

where the function $f(\theta)$ takes the expression

$\begin{split}&f(\theta)=\left(\sqrt{\xi\gamma}\sin\theta+\gamma\cos\varphi\right)^{2}\\\
&\quad\quad\quad+\left(-\sqrt{\frac{\xi}{\gamma}}\cos\theta+\frac{1}{\gamma}\sin\varphi\right)^{2}-1,\\\
&\gamma=\sqrt{\frac{m_{y}}{m_{x}}},\quad\epsilon_{0}=\frac{\left(q_{e}Ba_{0}\right)^{2}}{\sqrt{m_{x}m_{y}}},\quad\xi=\frac{\epsilon_{0}}{2\mu},\end{split}$
(B.9)

and $\left<\cdots\right>_{\text{FS}}$ is simply the average over
$\varphi\in[0,2\pi]$. A schematic plot of Eq. (B.8) and (B.9) in terms of the
resistivity $\rho^{n}_{c}=1/\sigma_{c}$ is shown in Fig. 12. As can be seen
from Fig. 12, although it presents the desired period and locations of
maxima/minima, it has a qualitatively different shape near the maxima and
minima. To make it even worse, in order to produce the amplitude of
oscillation in $\rho^{n}_{c}$ observed in experiments, we have to adopt an
effective electron mass on the order of $10^{-5}$ the bare electron mass,
which is unrealistic. In this sense, the experimental data favor the Dirac
fermions over the conventional fermions.

Figure 12: The schematic plot of the $c$-axis resistivity according to Eq.
(B.8) and Eq. (B.9), in correspondence to Fig. 4 in the main text.

##  Appendix C Materials and Methods

Here we discuss the materials and methods for obtaining the experimental data.
Single crystals of LBCO with $x=1/8$ studied here were grown in an infrared
image furnace by the floating-zone technique. They are pieces from the same
cylindrical crystal used previously to characterize two-dimensional
fluctuating superconductivity Li et al. (2007b). Single crystal samples were
cut and aligned into slabs, then fixed on a 0.5 mm thick sapphire substrate.
The imperfection in the sample alignment, estimated from X-ray diffraction, is
less than 0.5 degrees. For transport measurements, current contacts were made
at the ends of the longest dimension of crystals to ensure uniform current
flow, while the voltage contacts were made on both the top and side of the
crystals. We used a low-temperature contact annealing procedure Li et al.
(2007b) leading to low contact resistance $(<0.2~{}\Omega)$ that allows us to
measure the resistivity over seven orders of magnitude. The angular-resolved
magnetoresistance (ARMR) was measured using the 4-point probe in-line method
in a Quantum Design Physical Property Measurement System (PPMS) equipped with
a 14 T superconducting magnet. The resistivity measurements have been
performed with the current applied along either the in-plane direction or the
c-direction using DC and AC transport options with a current range of 50
$\mu\mathrm{A}-1\mathrm{~{}mA}$. Both DC and AC methods produced the same
results. The data shown are from the AC transport measurements
$(17\mathrm{~{}Hz})$. For crystal alignment with the magnetic field,
horizontal and vertical sample rotators with an angular resolution of $\sim
0.1^{\circ}$ were used. Temperature-dependent ARMR data were taken from $1.8$
K to $300\mathrm{~{}K}$, at various fields up to $14\mathrm{~{}T}$. ARMR data
at fixed temperatures and magnetic fields were taken in-situ with a vertical
sample rotator as a function of the in-plane magnetic field angles $(\theta)$
in a range of $-15^{\circ}$ to $355^{\circ}$.

## References

  * Li et al. (2007a) Q. Li, M. Hücker, G. D. Gu, A. M. Tsvelik, and J. M. Tranquada, Phys. Rev. Lett. 99, 067001 (2007a), URL https://link.aps.org/doi/10.1103/PhysRevLett.99.067001.
  * Tranquada et al. (2008) J. M. Tranquada, G. D. Gu, M. Hücker, Q. Jie, H.-J. Kang, R. Klingeler, Q. Li, N. Tristan, J. S. Wen, G. Y. Xu, et al., Phys. Rev. B 78, 174529 (2008), URL https://link.aps.org/doi/10.1103/PhysRevB.78.174529.
  * Berg et al. (2007) E. Berg, E. Fradkin, E.-A. Kim, S. A. Kivelson, V. Oganesyan, J. M. Tranquada, and S. C. Zhang, Phys. Rev. Lett. 99, 127003 (2007), URL https://link.aps.org/doi/10.1103/PhysRevLett.99.127003.
  * Li et al. (2019) Y. Li, J. Terzic, P. Baity, D. Popović, G. Gu, Q. Li, A. Tsvelik, and J. Tranquada, Science Advances 5, eaav7686 (2019), URL https://www.science.org/doi/abs/10.1126/sciadv.aav7686.
  * Shi et al. (2020a) Z. Shi, P. G. Baity, J. Terzic, T. Sasagawa, and D. Popović, Nature Communications 11, 3323 (2020a), URL https://doi.org/10.1038/s41467-020-17138-z.
  * Shi et al. (2020b) Z. Shi, P. G. Baity, T. Sasagawa, and D. Popović, Science Advances 6 (2020b), URL https://advances.sciencemag.org/content/6/7/eaay8946.
  * Lozano et al. (2022) P. M. Lozano, T. Ren, G. D. Gu, A. M. Tsvelik, J. M. Tranquada, and Q. Li, Phys. Rev. B 106, 174510 (2022), URL https://link.aps.org/doi/10.1103/PhysRevB.106.174510.
  * Blatter et al. (1994) G. Blatter, M. V. Feigel’man, V. B. Geshkenbein, A. I. Larkin, and V. M. Vinokur, Rev. Mod. Phys. 66, 1125 (1994), URL https://link.aps.org/doi/10.1103/RevModPhys.66.1125.
  * Soto et al. (2015) F. Soto, C. Carballeira, J. M. Doval, J. Mosqueira, M. V. Ramallo, A. Ramos-Álvarez, D. Sóñora, J. C. Verde, and F. Vidal, Superconductor Science and Technology 28, 065015 (2015), URL https://dx.doi.org/10.1088/0953-2048/28/6/065015.
  * Šopík et al. (2015) B. Šopík, J. Chaloupka, A. Dubroka, C. Bernhard, and D. Munzar, New Journal of Physics 17, 053022 (2015), URL https://dx.doi.org/10.1088/1367-2630/17/5/053022.
  * Agterberg and Tsunetsugu (2008) D. F. Agterberg and H. Tsunetsugu, Nature Physics 4, 639 (2008), URL https://doi.org/10.1038/nphys999.
  * Tsvelik (2016) A. M. Tsvelik, Phys. Rev. B 94, 165114 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.94.165114.
  * Tsvelik (2019) A. M. Tsvelik, Proceedings of the National Academy of Sciences 116, 12729 (2019), ISSN 0027-8424, URL https://www.pnas.org/content/116/26/12729.
  * Lee (2014) P. A. Lee, Phys. Rev. X 4, 031017 (2014), URL https://link.aps.org/doi/10.1103/PhysRevX.4.031017.
  * Das and Doniach (1999) D. Das and S. Doniach, Phys. Rev. B 60, 1261 (1999), URL https://link.aps.org/doi/10.1103/PhysRevB.60.1261.
  * Phillips (2018) P. W. Phillips, Topological Phase Transitions and New Developments (2018), arXiv: 1802.00451, URL http://arxiv.org/abs/1802.00451.
  * Ren and Tsvelik (2020) T. Ren and A. M. Tsvelik, New Journal of Physics 22, 103021 (2020), URL https://doi.org/10.1088/1367-2630/abbc2b.
  * Halperin and Nelson (1979) B. I. Halperin and D. R. Nelson, Journal of Low Temperature Physics 36, 599 (1979), URL https://doi.org/10.1007/BF00116988.
  * Minnhagen (1987) P. Minnhagen, Rev. Mod. Phys. 59, 1001 (1987), URL https://link.aps.org/doi/10.1103/RevModPhys.59.1001.
  * Ambegaokar and Halperin (1969) V. Ambegaokar and B. I. Halperin, Phys. Rev. Lett. 22, 1364 (1969), URL https://link.aps.org/doi/10.1103/PhysRevLett.22.1364.
  * Anderson (1966) P. W. Anderson, Rev. Mod. Phys. 38, 298 (1966), URL https://link.aps.org/doi/10.1103/RevModPhys.38.298.
  * Langer and Ambegaokar (1967) J. S. Langer and V. Ambegaokar, Phys. Rev. 164, 498 (1967), URL https://link.aps.org/doi/10.1103/PhysRev.164.498.
  * McCumber and Halperin (1970) D. E. McCumber and B. I. Halperin, Phys. Rev. B 1, 1054 (1970), URL https://link.aps.org/doi/10.1103/PhysRevB.1.1054.
  * Bardeen and Stephen (1965) J. Bardeen and M. J. Stephen, Phys. Rev. 140, A1197 (1965), URL https://link.aps.org/doi/10.1103/PhysRev.140.A1197.
  * Norman and Pépin (2003) M. R. Norman and C. Pépin, Reports on Progress in Physics 66, 1547 (2003), URL https://dx.doi.org/10.1088/0034-4885/66/10/R01.
  * Hücker et al. (2011) M. Hücker, M. v. Zimmermann, G. D. Gu, Z. J. Xu, J. S. Wen, G. Xu, H. J. Kang, A. Zheludev, and J. M. Tranquada, Phys. Rev. B 83, 104506 (2011), URL https://link.aps.org/doi/10.1103/PhysRevB.83.104506.
  * Affleck and Marston (1988) I. Affleck and J. B. Marston, Phys. Rev. B 37, 3774 (1988), URL https://link.aps.org/doi/10.1103/PhysRevB.37.3774.
  * Vafek and Vishwanath (2014) O. Vafek and A. Vishwanath, Annual Review of Condensed Matter Physics 5, 83 (2014), eprint https://doi.org/10.1146/annurev-conmatphys-031113-133841, URL https://doi.org/10.1146/annurev-conmatphys-031113-133841.
  * Horio et al. (2018) M. Horio, C. E. Matt, K. Kramer, D. Sutter, A. M. Cook, Y. Sassa, K. Hauser, M. Månsson, N. C. Plumb, M. Shi, et al., Nature Communications 9, 3252 (2018), URL https://doi.org/10.1038/s41467-018-05715-2.
  * Mulligan and Raghu (2016) M. Mulligan and S. Raghu, Phys. Rev. B 93, 205116 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.93.205116.
  * Mulligan (2017) M. Mulligan, Phys. Rev. B 95, 045118 (2017), URL https://link.aps.org/doi/10.1103/PhysRevB.95.045118.
  * Shon and Ando (1998) N. Shon and T. Ando, Journal of the Physical Society of Japan 67, 2421 (1998), URL https://doi.org/10.1143/JPSJ.67.2421.
  * Li et al. (2007b) Q. Li, M. Hücker, G. D. Gu, A. M. Tsvelik, and J. M. Tranquada, Phys. Rev. Lett. 99, 067001 (2007b), URL https://link.aps.org/doi/10.1103/PhysRevLett.99.067001.