# Transition density matrices of Richardson-Gaudin states

Paul A. Johnson<EMAIL_ADDRESS>Hubert Fortin  Samuel Cloutier
Charles-Émile Fecteau  Département de chimie, Université Laval, Québec,
Québec, G1V 0A6, Canada

###### Abstract

Recently, ground state eigenvectors of the reduced Bardeen-Cooper-Schrieffer
Hamiltonian, Richardson-Gaudin (RG) states, have been employed as a
wavefunction ansatz for strong correlation. This wavefunction physically
represents a mean-field of pairs of electrons (geminals) with a constant
pairing strength. To move beyond the mean-field, one must develop the
wavefunction in the basis of all the RG states. This requires both practical
expressions for transition density matrices and an idea of which states are
most important in the expansion. In this contribution, we present expressions
for the transition density matrix elements and calculate them numerically for
half-filled picket fence models. There are no Slater-Condon rules for RG
states, though an analogue of the aufbau principle proves to be useful in
choosing which states are important.

##  I Introduction

Strongly-correlated electrons remain a problem in quantum chemistry. Active
space methods are effective and expedient when the number of important Slater
determinants is small, but quickly become infeasible as the number of
important Slater determinants gets large. In such cases, it is more productive
to start from a different qualitative picture, such as weakly-interacting
pairs of electrons (geminals). Geminals are an old idea for quantum
chemistry,[1, 2, 3] though recent work has been quite promising.[4, 5, 6, 7,
8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]

The group of Scuseria are developing a many-body theory for pairs based on the
anti-symmetrized geminal power (AGP) as a mean-field. [24, 25, 26, 27, 28, 29,
30] Their approach has the advantage that the wavefunction is simple and
represents a clear physical picture: condensation of pairs. Their principal
application thus far has been the reduced Bardeen-Cooper-Schrieffer (BCS)[31,
32] pairing model.

Our approach is to employ the eigenvectors of the reduced BCS Hamiltonian, the
Richardson-Gaudin (RG) states,[33, 34, 35, 36] directly as a wavefunction
ansatz. Variational mean-field calculations have been performed,[37] optimal
energy expressions have been derived,[38] and the connection to single
reference coupled-cluster methods has been explored.[39] The principal
advantage of our approach is clear: a complete set of eigenvectors with which
to construct the physical wavefunction. Thus it will be straightforward,
though technically difficult, to develop perturbation theories, selective
configuration interactions, or coupled-cluster like approaches analogous to
those built upon Hartree-Fock (HF). The next clear step in this direction is
to obtain the necessary transition density matrix (TDM) elements between RG
states. In the HF case, there are Slater-Condon rules which means that only
singly- and doubly-excited Slater determinants have any non-zero coupling
through the molecular Hamiltonian. This is not the case for RG states and the
purpose of this contribution is to both develop expressions for the TDM
elements and to evaluate them numerically to demonstrate which RG states have
the largest coupling with the ground state through the Hamiltonian. The
methods to develop the equations are quite similar, though trickier, than for
the reduced density matrix(RDM) elements[38] so the presentation will be
brief. RG states have been used previously as a wavefunction model in
condensed-matter theory[40] and for configuration interaction in nuclear
structure theory.[41]

A brief summary of the most relevant properties of RG states is presented in
section II. Expressions for TDM elements in both the physical basis (PB) and
the Gaudin basis (GB) are presented in section III. TDM elements and
transition probabilities were computed numerically and the results are shown
in section IV.

##  II RG Eigenvectors

The reduced BCS Hamiltonian is built from the set of pair objects in the PB:

$\displaystyle
S^{+}_{i}=a^{\dagger}_{i\uparrow}a^{\dagger}_{i\downarrow},\quad
S^{-}_{i}=a_{i\downarrow}a_{i\uparrow},\quad
S^{z}_{i}=\frac{1}{2}\left(a^{\dagger}_{i\uparrow}a_{i\uparrow}+a^{\dagger}_{i\downarrow}a_{i\downarrow}-1\right)$
(1)

where $a^{\dagger}_{i\sigma}$ creates an electron in spatial orbital $i$ with
spin $\sigma$ etc. $S^{+}_{i}$ creates a pair in spatial orbital $i$,
$S^{-}_{i}$ removes a pair from spatial orbital $i$, and $S^{z}_{i}$ scales a
pair-occupied spatial orbital by $\frac{1}{2}$ and an empty spatial orbital by
$-\frac{1}{2}$. These objects have the structure of su(2):

$\displaystyle\left[S^{z}_{i},S^{\pm}_{j}\right]=\pm\delta_{ij}S^{\pm}_{i},\quad\left[S^{+}_{i},S^{-}_{j}\right]=2\delta_{ij}S^{z}_{i}.$
(2)

Rather than $S^{z}_{i}$, it is useful to employ the number operator

$\displaystyle\hat{n}_{i}=2S^{z}_{i}+1$ (3)

which counts the number of electrons in the spatial orbital $i$. The reduced
BCS Hamiltonian is written in terms of these objects as

$\displaystyle\hat{H}_{BCS}=\frac{1}{2}\sum_{i}\varepsilon_{i}\hat{n}_{i}-\frac{g}{2}\sum_{ij}S^{+}_{i}S^{-}_{j}$
(4)

and expresses a competition between filling the lowest energy spatial orbitals
with isotropic scattering of pairs between the spatial orbitals. The
particular form of the interaction, $-\frac{g}{2}$, is a consequence of the
choice made for the definition of $S^{z}(u)$ in (5).[42]

The eigenvectors of the reduced BCS Hamiltonian are built from a set of pair
objects in the GB:

$\displaystyle S^{+}(u)=\sum_{i}\frac{S^{+}_{i}}{u-\varepsilon_{i}},\quad
S^{-}(u)=\sum_{i}\frac{S^{-}_{i}}{u-\varepsilon_{i}},\quad
S^{z}(u)=\frac{1}{g}-\sum_{i}\frac{S^{z}_{i}}{u-\varepsilon_{i}}$ (5)

where $u$ is an arbitrary complex number. For distinct $u_{1},u_{2}$ these
objects have commutation relations

$\displaystyle[S^{z}(u_{1}),S^{\pm}(u_{2})]=\pm\frac{S^{\pm}(u_{1})-S^{\pm}(u_{2})}{u_{1}-u_{2}},\quad[S^{+}(u_{1}),S^{-}(u_{2})]=2\frac{S^{z}(u_{1})-S^{z}(u_{2})}{u_{1}-u_{2}},$
(6)

and for $u_{1}=u_{2}=u$ these may be analytically continued to give[43]:

$\displaystyle[S^{z}(u),S^{\pm}(u)]=\pm\frac{\partial S^{\pm}(u)}{\partial
u},\quad[S^{+}(u),S^{-}(u)]=2\frac{\partial S^{z}(u)}{\partial u}.$ (7)

The vacuum is denoted $\ket{\theta}$ and is usually the state with no
particles. It may also contain unpaired electrons that do not participate in
the pairing dynamics of the reduced BCS model. The vacuum is destroyed by the
lowering operators:

$\displaystyle S^{-}_{i}\ket{\theta}$ $\displaystyle=0$ (8a) $\displaystyle
S^{-}(u)\ket{\theta}$ $\displaystyle=0$ (8b)

and is an eigenvector of the spin measurement operators

$\displaystyle S^{z}_{i}\ket{\theta}$ $\displaystyle=-\frac{1}{2}\ket{\theta}$
(9a) $\displaystyle S^{z}(u)\ket{\theta}$
$\displaystyle=\alpha(u)\ket{\theta}.$ (9b)

The eigenvalue $\alpha(u)$ is important, so we highlight its value

$\displaystyle\alpha(u)=\frac{1}{g}+\frac{1}{2}\sum_{i}\frac{1}{u-\varepsilon_{i}}.$
(10)

Linearly independent vectors are produced by acting with the raising operators
on the vacuum. States of the type

$\displaystyle\ket{\Phi}=S^{+}_{1}S^{+}_{2}\dots S^{+}_{M}\ket{\theta}$ (11)

are Slater determinants, specifically closed shell (restricted Hartree-Fock)
though the pairing scheme may be engineered to allow for more general cases.
Acting twice with one $S^{+}_{i}$ results in zero.

Eigenvectors of the reduced BCS Hamiltonian have the form

$\displaystyle\ket{\\{v\\}}=S^{+}(v_{1})S^{+}(v_{2})\dots
S^{+}(v_{M})\ket{\theta}$ (12)

where the set of _rapidities_ $\\{v\\}$ are a solution of Richardson’s
equations

$\displaystyle\frac{2}{g}+\sum_{i}\frac{1}{v_{a}-\varepsilon_{i}}+\sum_{b\neq
a}\frac{2}{v_{b}-v_{a}}=0.$ (13)

Here we fix our notation: $\\{u\\}$ will always denote an arbitrary set of
complex numbers, $\\{v\\}$ denotes _a_ solution of Richardson’s equations and
$\\{w\\}$ denotes the _ground state solution_ of Richardson’s equations. The
eigenvalues of (12) are

$\displaystyle\hat{H}_{BCS}\ket{\\{v\\}}=\sum_{a}v_{a}\ket{\\{v\\}}$ (14)

and thus just the sum of the rapidities defining the state.

Our long term goal is develop the electronic wavefunction in the basis of RG
eigenvectors in the hope of having a short expansion when the correct mean-
field picture is of weakly-interacting pairs of electrons. The RG mean-field
has already been reported,[37] along with optimal expressions for the density
matrix elements required to compute it.[38] To calculate pertubative
corrections, we will need to evaluate coupling of distinct RG eigenvectors
through the Coulomb Hamiltonian

$\displaystyle\hat{H}_{C}=\sum_{ij}h_{ij}\sum_{\sigma}a^{\dagger}_{i\sigma}a_{j\sigma}+\frac{1}{2}\sum_{ijkl}V_{ijkl}\sum_{\sigma\tau}a^{\dagger}_{i\sigma}a^{\dagger}_{j\tau}a_{l\tau}a_{k\sigma}.$
(15)

where $\sigma$ and $\tau$ are spin labels. The molecular integrals have been
calculated in an orthonormal spatial orbital basis $\\{\phi\\}$, and the two-
body integrals are in _physicists’_ notation:

$\displaystyle h_{ij}$ $\displaystyle=\int
d\mathbf{r}\phi^{*}_{i}(\mathbf{r})\left(-\frac{1}{2}\nabla^{2}-\sum_{I}\frac{Z_{I}}{|\mathbf{r}-\mathbf{R}_{I}|}\right)\phi_{j}(\mathbf{r})$
(16) $\displaystyle V_{ijkl}$ $\displaystyle=\int
d\mathbf{r}_{1}d\mathbf{r}_{2}\frac{\phi^{*}_{i}(\mathbf{r}_{1})\phi^{*}_{j}(\mathbf{r}_{2})\phi_{k}(\mathbf{r}_{1})\phi_{l}(\mathbf{r}_{2})}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}.$
(17)

RG states have a definite number of unpaired electrons, or _seniority_. Most
often, they have no unpaired electrons, though non-interacting unpaired
electrons may be added. We will henceforth assume that all RG states have
seniority zero. No generality is lost as adding unpaired electrons affects all
the RG states in precisely the same manner. Thus, the only non-zero
contributions arise from the seniority-conserving or seniority-zero piece of
the Coulomb Hamiltonian:

$\displaystyle\hat{H}_{SZ}=\sum_{i}h_{ii}\hat{n}_{i}+\frac{1}{4}\sum_{i\neq
j}\left(2V_{ijij}-V_{ijji}\right)\hat{n}_{i}\hat{n}_{j}+\sum_{ij}V_{iijj}S^{+}_{i}S^{-}_{j}.$
(18)

In the physical basis, we will therefore require the (unnormalized) TDM
elements:

$\displaystyle D^{vw}_{i}$
$\displaystyle=\frac{1}{2}\braket{\\{v\\}}{\hat{n}_{i}}{\\{w\\}}$ (19)
$\displaystyle D^{vw}_{ij}$
$\displaystyle=\frac{1}{4}\braket{\\{v\\}}{\hat{n}_{i}\hat{n}_{j}}{\\{w\\}}$
(20) $\displaystyle P^{vw}_{ij}$
$\displaystyle=\braket{\\{v\\}}{S^{+}_{i}S^{-}_{j}}{\\{w\\}}$ (21)

We will also evaluate TDM elements in the Gaudin basis. The Gaudin pair
creators are a rectangular linear transformation of the local pair operators:

$\displaystyle\begin{pmatrix}S^{+}(u_{1})\\\ S^{+}(u_{2})\\\ \vdots\\\
S^{+}(u_{M})\end{pmatrix}=\begin{pmatrix}\frac{1}{u_{1}-\varepsilon_{1}}&\frac{1}{u_{1}-\varepsilon_{2}}&\dots&\frac{1}{u_{1}-\varepsilon_{N}}\\\
\frac{1}{u_{2}-\varepsilon_{1}}&\frac{1}{u_{2}-\varepsilon_{2}}&\dots&\frac{1}{u_{2}-\varepsilon_{N}}\\\
\vdots&\vdots&\ddots&\vdots\\\
\frac{1}{u_{M}-\varepsilon_{1}}&\frac{1}{u_{M}-\varepsilon_{2}}&\dots&\frac{1}{u_{M}-\varepsilon_{N}}\\\
\end{pmatrix}\begin{pmatrix}S^{+}_{1}\\\ S^{+}_{2}\\\ \vdots\\\
S^{+}_{N}\end{pmatrix}$ (22)

or more succinctly

$\displaystyle S^{+}(\mathbf{u})=C\mathbf{S}^{+}.$ (23)

The Cauchy matrix $C$ has a known[44] explicit right-inverse $C^{R}$

$\displaystyle[C^{R}]_{ia}=(u_{a}-\varepsilon_{i})\prod_{k\neq
i}\left(\frac{u_{a}-\varepsilon_{k}}{\varepsilon_{i}-\varepsilon_{k}}\right)\prod_{b\neq
a}\left(\frac{\varepsilon_{i}-u_{b}}{u_{a}-u_{b}}\right)$ (24)

such that $CC^{R}=I_{M}$. The system of equations (23) has more unknowns than
equations, thus has infinitely-many solutions (provided it has any). A
particular solution of (23) is

$\displaystyle\mathbf{S}^{+}=C^{R}S^{+}(\mathbf{u})$ (25)

as is easily verified by substitution. By the exact same argument, for the
Gaudin pair removal operator we arrive at

$\displaystyle\mathbf{S}^{-}=C^{R}S^{-}(\mathbf{u}).$ (26)

It is helpful to rewrite $S^{z}(u)$ as

$\displaystyle
S^{z}(u)=\alpha(u)-\frac{1}{2}\sum_{i}\frac{\hat{n}_{i}}{u-\varepsilon_{i}},$
(27)

with $\alpha(u)$ the vacuum eigenvalue (10) to get

$\displaystyle\hat{\mathbf{n}}=-2C^{R}\left(S^{z}(\mathbf{u})-\alpha(\mathbf{u})\right).$
(28)

The seniority-zero Hamiltonian can then be written in the GB as:

$\displaystyle\hat{H}_{SZ}(\\{u\\})=E_{0}+\sum_{a}\tilde{h}_{a}S^{z}(u_{a})+\sum_{ab}\tilde{V}^{z}_{ab}S^{z}(u_{a})S^{z}(u_{b})+\sum_{ab}\tilde{V}^{p}_{ab}S^{+}(u_{a})S^{-}(u_{b})$
(29)

with the transformed integrals

$\displaystyle E_{0}$
$\displaystyle=2\sum_{ia}h_{ii}C^{R}_{ia}\alpha(u_{a})+\sum_{i\neq
j}\sum_{ab}(2V_{ijij}-V_{ijji})C^{R}_{ia}C^{R}_{jb}\alpha(u_{a})\alpha(u_{b})$
(30) $\displaystyle\tilde{h}_{a}$
$\displaystyle=-2\sum_{i}h_{ii}C^{R}_{ia}-2\sum_{i\neq
j}\sum_{b}(2V_{ijij}-V_{ijji})C^{R}_{ia}C^{R}_{jb}\alpha(u_{b})$ (31)
$\displaystyle\tilde{V}^{z}_{ab}$ $\displaystyle=\sum_{i\neq
j}(2V_{ijij}-V_{ijji})C^{R}_{ia}C^{R}_{jb}$ (32)
$\displaystyle\tilde{V}^{p}_{ab}$
$\displaystyle=\sum_{ij}V_{iijj}C^{R}_{ia}C^{R}_{jb}.$ (33)

The Hamiltonian may be transformed to any set of rapidities $\\{u\\}$, which
may satisfy Richardson’s equations or not. The transformation depends only
upon the Cauchy structure. We will choose the ground state rapidities
$\\{w\\}$. Thus, in the Gaudin basis, we will evaluate the unnormalized TDM
elements:

$\displaystyle Z^{vw}(w_{a})$
$\displaystyle=\braket{\\{v\\}}{S^{z}(w_{a})}{\\{w\\}}$ (34) $\displaystyle
Z^{vw}(w_{a},w_{b})$
$\displaystyle=\braket{\\{v\\}}{S^{z}(w_{a})S^{z}(w_{b})}{\\{w\\}}$ (35)
$\displaystyle P^{vw}(w_{a},w_{b})$
$\displaystyle=\braket{\\{v\\}}{S^{+}(w_{a})S^{-}(w_{b})}{\\{w\\}}$ (36)

The notation employed here is different from that we employed for the RDM
elements in the GB. As we must choose a set of rapidities in which to work,
the present notation is meant to emphasize this choice.

###  II.1 Nature of RG excited states

Slater determinants are uniquely defined by which orbitals are occupied, with
the lowest in energy generally being that with the lowest energy orbitals
filled. This is not the case for RG states, but there is a useful analogue of
the aufbau principle, and one of the aims of the present contribution is to
reinforce it. Numerical algorithms to solve Richardson’s equations[45, 46, 47,
48, 49, 50, 51] begin with a coupling strength $g$ of zero, for which the
reduced BCS Hamiltonian is diagonal in a basis of Slater determinants, before
turning on the interaction adiabatically to evolve to the interacting RG
solutions. Thus, we can identify an RG state with its Slater determinant
initial guess as the evolution is smooth. The ground state is associated with
the lowest energy Slater determinant, defined by filling the $M$ lowest sites.
It is natural to refer to this RG state as a string of $M$ 1’s and $N-M$ 0’s.
The next lowest energy RG states evolve from Slater determinants which are
single pair excitations on top of the lowest energy Slater determinant. We
will refer to these RG states as _single excitations_ or _singles_ , and label
them individually with strings of 1’s and 0’s based on the occupations of the
pairs at zero coupling. This nomenclature easily extends to doubles, triples,
etc. Strings of 1’s and 0’s uniquely label each RG state: Richardson
established that the $\binom{N}{M}$ distinct solutions of Richardson’s
equations correspond to the $\binom{N}{M}$ distinct eigenvectors of the
reduced BCS Hamiltonian, so provided that each site can only hold one pair
this labelling scheme is a one-to-one correspondence. In the case of
degeneracies, the situation is a little more complicated, but has been solved
as well.[47]

In an analogous Bethe ansatz for electrons,[52] Richardson’s equations reduce
to $N$ independent conditions. In that case, each state is uniquely labelled
in the usual occupied/virtual sense familiar to HF. However, _each RG state is
defined by a different set of rapidities_ , which are coupled by Richardson’s
equations. As a result it is not possible to write RG excited states as the
action of a simple second-quantized operator acting on the RG ground state.

##  III Transition density matrices

###  III.1 Scalar products

Slavnov’s theorem[53, 54] expresses the scalar product between two RG states
with at least one of the sets (in this case $\\{v\\}$) a solution of
Richardson’s equations, while the set $\\{u\\}$ is arbitrary:

$\displaystyle\braket{\\{v\\}}{\\{u\\}}=K(\\{v\\},\\{u\\})\det
J(\\{v\\},\\{u\\})$ (37)

where $K(\\{v\\},\\{u\\})$ is 1 over the determinant of a Cauchy matrix

$\displaystyle
K(\\{v\\},\\{u\\})=\frac{\prod_{ab}(v_{a}-u_{b})}{\prod_{a<b}(u_{a}-u_{b})(v_{b}-v_{a})}$
(38)

and the elements of the matrix $J$ are

$\displaystyle
J_{ab}=\frac{1}{(v_{a}-u_{b})^{2}}\left(\frac{2}{g}+\sum_{i}\frac{1}{(u_{b}-\varepsilon_{i})}-\sum_{c\neq
a}\frac{2}{(u_{b}-v_{c})}\right).$ (39)

$K$ and $J$ are understood to be functions of two sets of rapidities. As this
dependence will always be clear, it will be suppressed. Slavnov’s result
original result concerned the more general case of the six-vertex model, while
the present expression was first written by Zhou et al.[55]

There are two particular limits of (37) of immediate relevance. The first
occurs when $\\{u\\}\rightarrow\\{v\\}$ in which case we obtain the square of
the norm. The result is the determinant of the Gaudin matrix $G^{v}$, which is
the Jacobian of Richardson’s equations:

$\displaystyle\braket{\\{v\\}}{\\{v\\}}$ $\displaystyle=\det G^{v}$ (40)
$\displaystyle G^{v}_{ab}$
$\displaystyle=\begin{cases}\sum_{i}\frac{1}{(v_{a}-\varepsilon_{i})^{2}}-\sum_{c\neq
a}\frac{2}{(v_{a}-v_{c})^{2}},\quad&a=b\\\
\frac{2}{(v_{a}-v_{b})^{2}},\quad&a\neq b.\end{cases}$ (41)

The norm of the state $\\{w\\}$ is similarly obtained, substituting $\\{w\\}$
for $\\{v\\}$.

The second case occurs when $\\{u\\}\rightarrow\\{w\\}$. In this case, we
obtain the scalar product between distinct (non-degenerate) eigenvectors of a
Hermitian operator, which must therefore be zero. Using Richardson’s
equations, one can see that the columns of $J$ become linearly dependent:

$\displaystyle c_{1}J_{1}+c_{2}J_{2}+\dots+c_{N}J_{N}=0$ (42)

with coefficients

$\displaystyle c_{k}=\prod_{a}(w_{k}-v_{a})^{2}$ (43)

so the determinant of $J$, and hence the overlap $\braket{\\{v\\}}{\\{w\\}}$,
is zero.

###  III.2 Physical basis

We follow the same approach as for the RDM case,[56, 38] so we pass directly
to the results. By moving the PB operators to the right past the Gaudin pair
creators until they act upon the vacuum, the TDM elements are:

$\displaystyle D^{vw}_{i}$
$\displaystyle=\sum_{a}\frac{F^{i}_{a}}{w_{a}-\varepsilon_{i}}$ (44a)
$\displaystyle D^{vw}_{ij}$ $\displaystyle=\sum_{a\neq
b}\frac{F^{ij}_{ab}}{(w_{a}-\varepsilon_{i})(w_{b}-\varepsilon_{j})}$ (44b)
$\displaystyle P^{vw}_{ij}$
$\displaystyle=\sum_{a}\frac{F^{i}_{a}}{w_{a}-\varepsilon_{i}}-\sum_{a\neq
b}\frac{F^{ij}_{ab}}{(w_{a}-\varepsilon_{j})(w_{b}-\varepsilon_{j})}$ (44c)

where the common elements are the form factors:

$\displaystyle F^{i}_{a}$
$\displaystyle=\braket{\\{v\\}}{S^{+}_{i}}{\\{w\\}_{a}}$ (45) $\displaystyle
F^{ij}_{ab}$
$\displaystyle=\braket{\\{v\\}}{S^{+}_{i}S^{+}_{j}}{\\{w\\}_{a,b}}.$ (46)

The notation $\\{w\\}_{a}$ denotes the set $\\{w\\}$ without the element
$w_{a}$ etc. As in the RDM case, the form factors are evaluated as limits of
Slavnov’s theorem. First, notice that the PB operators are the residues of the
Gaudin pair creators

$\displaystyle
S^{+}_{i}=\lim_{u_{a}\rightarrow\varepsilon_{i}}(u_{a}-\varepsilon_{i})S^{+}(u).$
(47)

Now, keep the set $\\{u\\}$ as arbitrary, and evaluate the corresponding
residues of the scalar products

$\displaystyle\braket{\\{v\\}}{S^{+}_{i}}{\\{u\\}_{a}}=\lim_{u_{a}\rightarrow\varepsilon_{i}}(u_{a}-\varepsilon_{i})\braket{\\{v\\}}{\\{u\\}}$
(48)

and take the limit $\\{u\\}\rightarrow\\{w\\}$ to get:

$\displaystyle
F^{i}_{a}=\frac{\prod_{c}(v_{c}-\varepsilon_{i})}{\prod_{c}(v_{c}-w_{a})}\frac{\prod_{b\neq
a}(w_{a}-w_{b})}{\prod_{b\neq a}(\varepsilon_{i}-w_{b})}K\det J^{i}_{a}$ (49)

Here, $J^{i}_{a}$ is the matrix $J$ with the $a$th column replaced by the
column:

$\displaystyle\mathbf{q}_{i}=\begin{pmatrix}\frac{1}{(v_{1}-\varepsilon_{i})^{2}}\\\
\frac{1}{(v_{2}-\varepsilon_{i})^{2}}\\\ \vdots\\\
\frac{1}{(v_{M}-\varepsilon_{i})^{2}}\end{pmatrix}.$ (50)

Similarly,

$\displaystyle\braket{\\{v\\}}{S^{+}_{i}S^{+}_{j}}{\\{u\\}_{a,b}}=\lim_{u_{a}\rightarrow\varepsilon_{i}}\lim_{u_{b}\rightarrow\varepsilon_{j}}(u_{a}-\varepsilon_{i})(u_{b}-\varepsilon_{j})\braket{\\{v\\}}{\\{u\\}}$
(51)

and taking the limit $\\{u\\}\rightarrow\\{w\\}$ gives

$\displaystyle
F^{ij}_{ab}=\frac{(w_{a}-w_{b})}{(\varepsilon_{i}-\varepsilon_{j})}\frac{\prod_{c}(v_{d}-\varepsilon_{i})(v_{d}-\varepsilon_{j})}{\prod_{c}(v_{d}-w_{a})(v_{d}-w_{b})}\frac{\prod_{d\neq
a,b}(w_{a}-w_{d})(w_{b}-w_{d})}{\prod_{d\neq
a,b}(\varepsilon_{i}-w_{d})(\varepsilon_{j}-w_{d})}K\det J^{ij}_{ab},$ (52)

where $J^{ij}_{ab}$ is the matrix $J$ with the $a$th column replaced by the
$i$th version of (50) and the $b$th column replaced by the $j$th version of
(50).

Thus, to compute the TDM elements in the physical basis, we compute the form
factors (49) and (52) then compute the sums (44). In the RDM case, substantial
further simplification was possible due to Cramer’s rule and a theorem of
Jacobi[57] for scaled minors: dividing by the norm of the state yielded ratios
of determinants differing by one or two columns. In this case, the form
factors differ substantially from the Gaudin matrices. For $\\{u\\}$
arbitrary, there is a relevant Jacobi identity

$\displaystyle\frac{\det J^{ij}_{ab}}{\det J}=\frac{\det J^{i}_{a}}{\det
J}\frac{\det J^{j}_{b}}{\det J}-\frac{\det J^{i}_{b}}{\det J}\frac{\det
J^{j}_{a}}{\det J}$ (53)

but when $\\{u\\}\rightarrow\\{w\\}$, $J$ becomes singular and this identity
no longer holds. Therefore, while it may be possible to simplify the TDM
elements, for this contribution we content ourselves with (44) to evaluate
numerically.

The result for $D^{vw}_{ij}$ is valid provided $i\neq j$. The diagonal element
corresponds to $D^{vw}_{i}$ (and to $P^{vw}_{ii}$). It is not difficult to see
that $D^{vw}_{i}$ and $D^{vw}_{ij}$ sum to zero.

$\displaystyle\sum_{i}D^{vw}_{i}$ $\displaystyle=0$ (54)
$\displaystyle\sum_{ij}D^{vw}_{ij}$ $\displaystyle=0$ (55)

###  III.3 Gaudin basis

TDM elements in the GB are obtained in a similar manner to that for the RDM
elements: first use the Gaudin algebra to write the TDM elements as a sum of
form factors, then evaluate the form factors. As the Gaudin algebra actions
are exactly the same as in the RDM case, they are listed in appendix A. The
form factors required are more tricky than for the RDM case, but the same
procedure is to be followed. Expressions for the form factors are listed in
appendix B.

The final expression for the $S^{z}(w_{a})$ TDM elements in the GB is:

$\displaystyle Z^{vw}(w_{a})$ $\displaystyle=K\sum_{c}\lambda^{a}_{c}\det
J^{a}_{c}$ (56)

with

$\displaystyle\lambda^{a}_{c}=\frac{\prod_{k}(v_{k}-w_{a})}{\prod_{k}(v_{k}-w_{c})}\frac{\prod_{l\neq
c}(w_{c}-w_{l})}{\prod_{l\neq a}(w_{a}-w_{l})},$ (57)

and the matrix $J^{a}_{c}$ is the matrix $J$ with the $c$th column replaced
with the derivative of the $a$th column with respect to $w_{a}$:
$\frac{\partial J_{a}}{\partial w_{a}}$. The derivatives are:

$\displaystyle\left(\frac{\partial J_{a}}{\partial w_{a}}\right)_{b}$
$\displaystyle=\frac{2}{(v_{b}-w_{a})^{3}}\left(\frac{2}{g}+\sum_{i}\frac{1}{(w_{a}-\varepsilon_{i})}-\sum_{c\neq
b}\frac{2}{(w_{a}-v_{c})}\right)$
$\displaystyle-\frac{1}{(v_{b}-w_{a})^{2}}\left(\sum_{i}\frac{1}{(w_{a}-\varepsilon_{i})^{2}}-\sum_{c\neq
b}\frac{2}{(w_{a}-v_{c})^{2}}\right)$ (58)

Clearly, $\lambda^{a}_{a}=1$. Further, with

$\displaystyle k(w_{a})=\sum_{c}\frac{1}{w_{a}-v_{c}}-\sum_{c\neq
a}\frac{1}{w_{a}-w_{c}}$ (59)

and the elements of the second derivatives of $J$

$\displaystyle\frac{1}{2}\left(\frac{\partial^{2}J_{a}}{\partial
w_{a}^{2}}\right)_{b}$
$\displaystyle=\frac{3}{(v_{b}-w_{a})^{4}}\left(\frac{2}{g}+\sum_{i}\frac{1}{(w_{a}-\varepsilon_{i})}-\sum_{c\neq
b}\frac{2}{(w_{a}-v_{c})}\right)$
$\displaystyle-\frac{2}{(v_{b}-w_{a})^{3}}\left(\sum_{i}\frac{1}{(w_{a}-\varepsilon_{i})^{2}}-\sum_{c\neq
b}\frac{2}{(w_{a}-v_{c})^{2}}\right)$
$\displaystyle+\frac{1}{(v_{b}-w_{a})^{2}}\left(\sum_{i}\frac{1}{(w_{a}-\varepsilon_{i})^{3}}-\sum_{c\neq
b}\frac{2}{(w_{a}-v_{c})^{3}}\right)$ (60)

the expressions for the other transition density matrix elements are:

$\displaystyle\frac{1}{K}Z^{vw}(w_{a},w_{a})$ $\displaystyle=\det
J^{\bar{a}}_{a}+k(w_{a})\det J^{a}_{a}$ $\displaystyle+\sum_{c\neq
a}\lambda^{a}_{c}\left(3\det
J^{\bar{a}}_{c}+\left(3k(w_{a})+\frac{4}{(w_{a}-w_{c})}\right)\det
J^{a}_{c}\right)$ $\displaystyle+\sum_{c,d\neq
a}\lambda^{a}_{c}\lambda^{a}_{d}\frac{(w_{a}-w_{c})(w_{a}-w_{d})}{(w_{d}-w_{c})}\det
J^{\bar{a}a}_{cd}$ (61)

$\displaystyle\frac{1}{K}Z^{vw}(w_{a},w_{b})$ $\displaystyle=\det
J^{ab}_{ab}+k(w_{a})\det J^{b}_{b}+k(w_{b})\det
J^{a}_{a}+\frac{\left(\lambda^{a}_{b}\det J^{a}_{b}-\lambda^{b}_{a}\det
J^{b}_{a}\right)}{(w_{a}-w_{b})}$ $\displaystyle+\sum_{c\neq
a,b}\lambda^{b}_{c}\left(\det
J^{ab}_{ac}+\left(k(w_{a})+\frac{3(w_{c}-w_{b})}{(w_{a}-w_{c})(w_{a}-w_{b})}\right)\det
J^{b}_{c}\right)$ $\displaystyle+\sum_{c\neq a,b}\lambda^{a}_{c}\left(\det
J^{ba}_{bc}+\left(k(w_{b})+\frac{3(w_{c}-w_{a})}{(w_{b}-w_{c})(w_{b}-w_{a})}\right)\det
J^{a}_{c}\right)$ $\displaystyle+\frac{1}{2}\sum_{c,d\neq
a,b}\lambda^{a}_{c}\lambda^{b}_{d}\frac{(w_{a}-w_{d})(w_{b}-w_{c})+(w_{a}-w_{c})(w_{b}-w_{d})}{(w_{a}-w_{b})(w_{d}-w_{c})}\det
J^{ab}_{cd}$ (62)

$\displaystyle\frac{1}{K}P^{vw}(w_{a},w_{a})$ $\displaystyle=-2\sum_{c\neq
a}\lambda^{a}_{c}\left(\det
J^{\bar{a}}_{c}+\left(k(w_{a})-\frac{2}{(w_{c}-w_{a})}\right)\det
J^{a}_{c}\right)$ $\displaystyle-\sum_{c,d\neq
a}\lambda^{a}_{c}\lambda^{a}_{d}\frac{(w_{a}-w_{c})(w_{a}-w_{d})}{(w_{d}-w_{c})}\det
J^{\bar{a}a}_{cd}$ (63)

$\displaystyle\frac{1}{K}P^{vw}(w_{a},w_{b})$
$\displaystyle=\left(\sum_{i}\frac{1}{(w_{b}-\varepsilon_{i})^{2}}-\sum_{c\neq
b}\frac{2}{(w_{c}-w_{b})^{2}}\right)\lambda^{a}_{b}(w_{a}-w_{b})\det
J^{a}_{b}+\frac{2}{w_{a}-w_{b}}\det J^{b}_{b}$ $\displaystyle-2\sum_{c\neq
a,b}\frac{(w_{c}-w_{a})}{(w_{c}-w_{b})}\lambda^{a}_{c}\left(\det
J^{ba}_{bc}+\left(k(w_{b})-\frac{(w_{c}-w_{a})+(w_{b}-w_{a})}{(w_{c}-w_{b})(w_{b}-w_{a})}\right)\det
J^{a}_{c}\right)$ $\displaystyle-2\sum_{c\neq a,b}\frac{\lambda^{b}_{c}\det
J^{b}_{c}}{w_{b}-w_{a}}+\sum_{c,d\neq
b}\lambda^{a}_{c}\lambda^{b}_{d}\frac{(w_{a}-w_{c})(w_{a}-w_{d})}{(w_{a}-w_{b})(w_{c}-w_{d})}\det
J^{ab}_{cd}$ (64)

$J^{ab}_{cd}$ is the matrix $J$ with the $c$th column replaced with
$\frac{\partial J_{a}}{\partial w_{a}}$ and the $d$th column replaced with
$\frac{\partial J_{b}}{\partial w_{b}}$. Similarly, $J^{\bar{a}}_{c}$ is $J$
with the $c$th column replaced by
$\frac{1}{2}\frac{\partial^{2}J_{a}}{\partial w_{a}^{2}}$, etc. Again, the
expressions for the TDM elements may not be optimal as further simplification
from Cramer’s rule and Jacobi’s identity are not possible as in the RDM case.

###  III.4 Consistency checks

The TDM elements in the GB have been verified numerically against direct
transformations of the corresponding elements in the PB. Specifically,

$\displaystyle
P^{vw}(w_{a},w_{b})=\sum_{ij}\frac{P^{vw}_{ij}}{(w_{a}-\varepsilon_{i})(w_{b}-\varepsilon_{j})}$
(65)

$\displaystyle
Z^{vw}(w_{a})=-\sum_{i}\frac{D^{vw}_{i}}{w_{a}-\varepsilon_{i}}$ (66)

and

$\displaystyle
Z^{vw}(w_{a},w_{b})=\sum_{ij}\frac{D^{vw}_{ij}-\frac{1}{2}D^{vw}_{i}-\frac{1}{2}D^{vw}_{j}}{(w_{a}-\varepsilon_{i})(w_{b}-\varepsilon_{j})}-\frac{1}{g}\sum_{i}\frac{D^{vw}_{i}}{w_{a}-\varepsilon_{i}}-\frac{1}{g}\sum_{j}\frac{D^{vw}_{j}}{w_{b}-\varepsilon_{j}}.$
(67)

The TDM elements in the GB satisfy sum rules, which also serve as a
consistency check. With the eigenvalue based variables

$\displaystyle W_{i}=\sum_{a}\frac{1}{\varepsilon_{i}-w_{a}},$ (68)

the sum of the $Z^{vw}(w_{a})$ elements is

$\displaystyle\sum_{a}Z^{vw}(w_{a})=\sum_{i}W_{i}D^{vw}_{i}$ (69)

the sum of the $P^{vw}(w_{a},w_{b})$ elements gives,

$\displaystyle\sum_{ab}P^{vw}(w_{a},w_{b})=\sum_{ij}W_{i}W_{j}P^{vw}_{ij}$
(70)

and the sum of the $Z^{vw}(w_{a},w_{b})$ elements is

$\displaystyle\sum_{ab}Z^{vw}(w_{a},w_{b})=\frac{2M}{g}\sum_{i}W_{i}D^{vw}_{i}+\sum_{ij}W_{i}W_{j}\left(D^{vw}_{ij}-\frac{1}{2}D^{vw}_{i}-\frac{1}{2}D^{vw}_{j}\right).$
(71)

Defining

$\displaystyle\Delta^{vw}(w_{a})=\braket{\\{v\\}}{\frac{\partial
S^{z}(w_{a})}{\partial w_{a}}}{\\{w\\}}=K\sum_{c}\lambda^{a}_{c}\left(\det
J^{\bar{a}}_{c}+k(w_{a})\det J^{a}_{c}\right),$ (72)

another consistency check can be performed since $\ket{\\{v\\}}$ and
$\ket{\\{w\\}}$ are both eigenvectors of $S^{2}(w_{a})$, so that

$\displaystyle\braket{\\{v\\}}{S^{2}(w_{a})}{\\{w\\}}=0.$ (73)

This is easily verified, as

$\displaystyle S^{2}(w_{a})$
$\displaystyle=S^{z}(w_{a})S^{z}(w_{a})+S^{+}(w_{a})S^{-}(w_{a})-\frac{\partial
S^{z}(w_{a})}{\partial w_{a}}$ (74)

and therefore,

$\displaystyle\braket{\\{v\\}}{S^{2}(w_{a})}{\\{w\\}}=Z^{vw}(w_{a},w_{a})+P^{vw}(w_{a},w_{a})-\Delta^{vw}(w_{a})=0$
(75)

by a direct computation.

##  IV Numerical Results

To understand which RG excited states would be most important in a
perturbation theory, we have numerically computed the _normalized_ TDM
elements in both the PB and the GB for a picket-fence model, i.e. the reduced
BCS model with constant energy spacing between the levels. The dynamics of
this model are a function of the ratio of the pairing strength $g$ to the
spacing $\Delta\varepsilon$. In our chosen convention (4), the pairing
interaction is $-\frac{g}{2}$ and thus a negative pairing interaction
corresponds to a repulsive interaction.

Rather than plot individual TDM elements, we present sums of the corresponding
transition probabilities, i.e. the sum of the absolute squares of the TDM
elements. Transition probabilities are real and positive definite, and are
thus useful to measure relative importance of RG excited states. Further, in a
Rayleigh-Schrödinger pertubation theory (RSPT) expansion, contributions from
excited states are weighted with the difference in zeroth-order energy
differences between the two states. Using the reduced BCS model as a starting
point, the energy denominators are the sums of differences of rapidities.
Therefore, we also plot transition probabilities with energy denominators to
judge which RG excited states would be important for RSPT.

All curves are plotted for small systems, half-filled 4- and 6-site models,
while excitations are grouped together for the half-filled 10-site model. The
size of the system is not a problem, but the number of states grows
combinatorially. Additional results for half-filled 8-site model are included
in the supplementary material.

###  IV.1 Physical Basis

We compute the TDM elements in the PB from the results of section III before
taking their absolute squares to plot transition probabilities. In particular,
we plot $\sum_{i}|D^{vw}_{i}|^{2}$, $\sum_{ij}|D^{vw}_{ij}|^{2}$ and
$\sum_{ij}|P^{vw}_{ij}|^{2}$, along with their weighting by energy
denominators: $\frac{\sum_{i}|D^{vw}_{i}|^{2}}{\sum_{a}w_{a}-v_{a}}$,
$\frac{\sum_{ij}|D^{vw}_{ij}|^{2}}{\sum_{a}w_{a}-v_{a}}$ and
$\frac{\sum_{ij}|P^{vw}_{ij}|^{2}}{\sum_{a}w_{a}-v_{a}}$.

Figure 1: a-c) Sums of transition probabilities for half-filled 4-site picket
fence model, with energy denominators included d-f).

Transition probabilities for the half-filled 4-site picket fence model are
shown in figure 1. In each plot the curves for the states 1001 and 0110 are
indistinguishable. The TDM elements and transition probabilities for these
states are _not_ identical, but their sums are. For the $D^{vw}_{i}$ elements,
the important contributions come from the single excitations, while the double
excitation gives almost no contribution. In the attractive regime the most
important elements come from the singles close to the Fermi level. We will
refer to these excitations as Fermi states. In the repulsive regime the most
important contribution is from 0101, which puts the lowest energy pair in the
highest level at zero coupling. We will refer to these types of states as
Rydberg states. For $D^{vw}_{ij}$ at weak positive coupling the lowest Fermi
state is the most important, but contributions from all excited states become
the same at large positive coupling. In the repulsive regime, the highest
Rydberg single is the most important, followed by the double excitation. For
$P^{vw}_{ij}$, the singles all pass through 1 at zero coupling, and are more
important than the double at all couplings.

Adding the energy denominator to $D^{vw}_{i}$, the state 1010 is the most
important at weak couplings but becomes asymptotically the same as 1001 and
0110 in both attractive and repulsive limits. The Rydberg state 0101 becomes
the most important in the repulsive regime. The double excitation gives almost
no contribution. All transition probabilities are finite at all couplings.
Including the energy denominator for $D^{vw}_{ij}$ makes 1010 the only
important state in the repulsive regime: the numerator remains constant but
the denominator goes to zero linearly. This is because the rapidities for the
ground state and for the state 1010 asymptotically approach each other at
large negative couplings. In the attractive regime everything remains finite
as the rapidities for each state remain distinct (and in fact become large).
The same is true for $P^{vw}_{ij}$.

Figure 2: Sums of transition probabilities of half-filled 6-site picket fence
model for a-c) single, d-f) double and g-i) triple excitations.

Transition probabilities for the half-filled 6-site model are shown in figure
2: the first row shows the single excitations, the second row the doubles, and
the third row the lone triple excitation. As for the half-filled 4 site model,
there are overlapping curves. In particular, each state is degenerate with the
state obtained by reading the 1’s and 0’s backwards then exchanging 1’s for
0’s. For example, 110010 gives the same contribution as 101100, 110001 gives
the same contribution as 011100 etc. Concretely, the green curves are hidden
beneath the cyan curves, the red curves are hidden behing the maroon curves,
and the yellow curves are hidden behind the orange curves. For $D^{vw}_{i}$,
only the singles matter, with the Fermi states being the most important in the
attractive limit and the Rydberg states the most important in the repulsive
regime. The same states are the most important for $D^{vw}_{ij}$ but now the
highest Rydberg double excitation is the second most important contribution in
the repulsive limit. The Fermi doubles are important in the attractive limit.
For $P^{vw}_{ij}$ all the singles are important, along with all the doubles in
the repulsive limit and the Fermi doubles in the attractive limit. The triple
excitation gives no important contributions.

Figure 3: Sums of transition probabilities with energy denominators of half-
filled 6-site picket fence model for a-c) single, d-f) double and g-i) triple
excitations.

With energy denominators, shown in figure 3, $D^{vw}_{i}$ is dominated by the
singles and is finite at all couplings. $D^{vw}_{ij}$ and $P^{vw}_{ij}$ are
dominated by the singles, with a few going linearly to minus infinity.

Figure 4: a-c) Sums of transition probabilities for half-filled 10-site picket
fence model, with energy denominators included for d-f).

Transition probabilities grouped by excitation order for the half-filled
10-site picket fence model are shown in figure 4. For $D^{vw}_{i}$ only the
singles matter at all, and with energy denominators included everything
remains finite at all couplings. For $D^{vw}_{ij}$ and $P^{vw}_{ij}$ the
singles are the most important, the doubles are less important, while triples
and above give almost no contribution. With energy denominators, the singles
diverge linearly at strong negative couplings.

Before passing to the GB, we summarize the principal results for the PB.
$D^{vw}_{i}$ is entirely dominated by the singles, and with an energy
denominator remains finite at all couplings. $D^{vw}_{ij}$ and $P^{vw}_{ij}$
are dominated by the singles, though the doubles are relevant. With energy
denominators, the singles for both $D^{vw}_{ij}$ and $P^{vw}_{ij}$ go to minus
infinity linearly in the repulsive limit as the states become near degenerate
while the numerators remain constant. Both are well behaved in the attractive
limit.

###  IV.2 Gaudin Basis

Directly plotting TDM elements or transition probabilities in the GB leads to
divergences at critical points. At these points, here characterized by a
specific $g$, the nature of the solutions of Richardson’s equations changes
from two distinct real numbers to a complex conjugate pair. _At_ a critical
point, two rapidities are equal to one of the single-particle energies
$\varepsilon$ so the TDM elements themselves diverge. However, in the
transformation to the GB, the elements of the inverse of the Cauchy matrix
manage these divergences so that the result is always finite. Thus, we may
absorb some of the transformation using the factor

$\displaystyle q(w_{a})=\frac{\prod_{i}(w_{a}-\varepsilon_{i})}{\prod_{b\neq
a}(w_{a}-w_{b})}$ (76)

to regularize the TDM elements:

$\displaystyle\bar{Z}^{vw}(w_{a})$ $\displaystyle=q(w_{a})Z(w_{a})$ (77)
$\displaystyle\bar{Z}^{vw}(w_{a},w_{b})$
$\displaystyle=q(w_{a})q(w_{b})Z(w_{a},w_{b})$ (78)
$\displaystyle\bar{P}^{vw}(w_{a},w_{b})$
$\displaystyle=q(w_{a})q(w_{b})P(w_{a},w_{b}).$ (79)

Specifically, we plot normalized $\sum_{a}|\bar{Z}^{vw}(w_{a})|^{2}$,
$\sum_{ab}|\bar{Z}^{vw}(w_{a},w_{b})|^{2}$ and
$\sum_{ab}|\bar{P}^{vw}(w_{a},w_{b})|^{2}$.

Figure 5: Sums of regularized Gaudin basis transition probabilities for half-
filled 4-site picket fence model.

Regularized transition probabilities for the half-filled 4-site picket fence
model are presented in figure 5. At large positive $g$, the attractive pairing
regime, the rapidities become large, and thus the transition probabilities
grow like large even polynomials. We therefore plot the results
logarithmically. The cusps in the plots occur at the critical points where the
nature of the rapidities change. Incorporating the energy denominator does not
change the appearance of the plots. For all three types of TDM elements, the
singles 1010, 1001 and 0110 are the most important in the attractive case, and
the double 0011 becomes more important than the single 0101 in the strongly
attractive limit. Each state is important in the strongly repulsive limit, but
to different TDM elements.

Figure 6: Sums of regularized Gaudin basis transition probabilities for half-
filled 6-site picket fence model for a-c) single, d-f) double and g-i) triple
excitations.

Regularized transition probabilities for the half-filled 6-site picket fence
model are presented in figure 6. In the repulsive regime, the three most
important states are singles with a pair in the highest level. In the
attractive regime, the Fermi singles are dominant, while the remaining singles
are about even with the most important doubles. The triple 000111 gives no
meaningful contribution.

Figure 7: Sums of regularized Gaudin basis transition probabilities for half-
filled 10-site picket fence model.

Regularized transition probabilities grouped by excitation level for the half-
filled 10-site model are shown in figure 7. At all couplings the singles give
the largest contribution, followed by the doubles etc.

##  V Conclusion

The purpose of this contribution is to demonstrate that while there are not
Slater-Condon rules for the TDM elements of RG states, the situation is still
manageable without drastic approximations. In both the physical and Gaudin
bases, contributions from the single pair excitations are by far the most
important, double pair excitations should be included, but triples and beyond
are much less important, especially in the physical basis. The dominant
excitations in the attractive regime are the Fermi states while for the
repulsive pairing regime the Rydberg states are most important. In the
physical basis, $|D^{vw}_{ij}|^{2}$ and $|P^{vw}_{ij}|^{2}$ singles with
energy denominators go linearly to minus infinity in the purely repulsive
limit. In the Gaudin basis, all transition probabilities behave like
polynomials in the attractive regime. Practically, these results suggest that
for a pertubation theory it would be better to use the physical basis for
attractive pairing strengths, and the Gaudin basis for repulsive pairing
strengths. Generally, there is a clear decrease in importance with respect to
excitation level, though the physical basis may be preferable as it avoids the
rectangular transformation to the Gaudin basis. In the following
contributions, we will employ these results to develop a perturbation theory
starting from an RG mean-field.

##  VI Acknowledgements

We thank the Natural Sciences and Engineering Research Council of Canada for
their financial support. This research was enabled in part by support from
Calcul Québec and Compute Canada. We also thank Alexandre Faribault for useful
discussions concerning critical points of Richardson’s equations.

##  Appendix A Gaudin Algebra actions

The actions of Gaudin algebra elements on an RG eigenvector were reported
previously,[38] so we present only the final results here. For a set of
$\\{w\\}$ that satisfy Richardson’s equations, the relevant results are:

$\displaystyle S^{z}(w_{a})\ket{\\{w\\}}=\frac{\partial S^{+}(w_{a})}{\partial
w_{a}}\ket{\\{w\\}_{a}}+\sum_{b\neq
a}\frac{S^{+}(w_{a})}{w_{a}-w_{b}}\ket{\\{w\\}_{b}}$ (80)

$\displaystyle S^{z}(w_{a})S^{z}(w_{a})\ket{\\{w\\}}$
$\displaystyle=\frac{1}{2}\frac{\partial^{2}S^{+}(w_{a})}{\partial
w_{a}^{2}}\ket{\\{w\\}_{a}}+\sum_{c,d\neq
a}\frac{S^{+}(w_{a})S^{+}(w_{a})}{(w_{a}-w_{c})(w_{a}-w_{d})}\ket{\\{w\\}_{c,d}}$
$\displaystyle+\sum_{c\neq a}\frac{3}{w_{a}-w_{c}}\frac{\partial
S^{+}(w_{a})}{\partial w_{a}}\ket{\\{w\\}_{c}}+\sum_{c\neq
a}\frac{S^{+}(w_{a})}{(w_{a}-w_{c})^{2}}\ket{\\{w\\}_{c}}$ (81)

$\displaystyle S^{z}(w_{a})S^{z}(w_{b})\ket{\\{w\\}}$
$\displaystyle=\frac{\partial S^{+}(w_{a})}{\partial w_{a}}\frac{\partial
S^{+}(w_{b})}{\partial
w_{b}}\ket{\\{w\\}_{a,b}}-\frac{3}{(w_{a}-w_{b})^{2}}\ket{\\{w\\}}$
$\displaystyle+\sum_{c\neq a,b}\left(\frac{1}{w_{b}-w_{c}}\frac{\partial
S^{+}(w_{a})}{\partial
w_{a}}S^{+}(w_{b})\ket{\\{w\\}_{a,c}}+\frac{1}{w_{a}-w_{c}}\frac{\partial
S^{+}(w_{b})}{\partial w_{b}}S^{+}(w_{a})\ket{\\{w\\}_{b,c}}\right)$
$\displaystyle+\sum_{c\neq
a,b}\left(\frac{2}{(w_{a}-w_{b})(w_{b}-w_{c})}S^{+}(w_{a})\ket{\\{w\\}_{c}}+\frac{2}{(w_{b}-w_{a})(w_{a}-w_{c})}S^{+}(w_{b})\ket{\\{w\\}_{c}}\right)$
$\displaystyle+\frac{1}{2}\sum_{c,d\neq
a,b}\left(\frac{1}{(w_{a}-w_{c})(w_{d}-w_{b})}+\frac{1}{(w_{a}-w_{d})(w_{b}-w_{c})}\right)S^{+}(w_{a})S^{+}(w_{b})\ket{\\{w\\}_{c,d}}$
$\displaystyle+\frac{1}{(w_{a}-w_{b})^{2}}S^{+}(w_{a})\ket{\\{w\\}_{b}}+\frac{1}{(w_{a}-w_{b})^{2}}S^{+}(w_{b})\ket{\\{w\\}_{a}}.$
(82)

$\displaystyle S^{+}(w_{a})S^{-}(w_{a})\ket{\\{w\\}}$
$\displaystyle=-2\left(\frac{\partial\alpha(w_{a})}{\partial
w_{a}}+\sum_{c\neq a}\frac{1}{(w_{c}-w_{a})^{2}}\right)\ket{\\{w\\}}$
$\displaystyle+2\sum_{c\neq a}\frac{1}{w_{c}-w_{a}}\frac{\partial
S^{+}(w_{a})}{\partial w_{a}}\ket{\\{w\\}_{c}}-2\sum_{c\neq
a}\frac{S^{+}(w_{a})}{(w_{c}-w_{a})^{2}}\ket{\\{w\\}_{c}}$
$\displaystyle-\sum_{c,d\neq
a}\frac{S^{+}(w_{a})S^{+}(w_{a})}{(w_{a}-w_{c})(w_{a}-w_{d})}\ket{\\{w\\}_{c,d}}.$
(83)

$\displaystyle S^{+}(w_{a})S^{-}(w_{b})\ket{\\{w\\}}$
$\displaystyle=-2\left(\frac{\partial\alpha(w_{b})}{\partial
w_{b}}+\sum_{c\neq
b}\frac{1}{(w_{c}-w_{b})^{2}}\right)S^{+}(w_{a})\ket{\\{w\\}_{b}}$
$\displaystyle+2\sum_{c\neq b}\frac{S^{+}(w_{a})}{w_{c}-w_{b}}\frac{\partial
S^{+}(w_{b})}{\partial w_{b}}\ket{\\{w\\}_{b,c}}-2\sum_{c\neq
b}\frac{S^{+}(w_{a})}{(w_{c}-w_{b})^{2}}\ket{\\{w\\}_{c}}$
$\displaystyle-\sum_{c,d\neq
b}\frac{S^{+}(w_{a})S^{+}(w_{b})}{(w_{b}-w_{c})(w_{b}-w_{d})}\ket{\\{w\\}_{c,d}}$
(84)

##  Appendix B Gaudin Basis Form factors

Form factors are calculated in exactly the same manner as for the RDM case:
begin from Slavnov’s theorem, set the free variables $u_{a}=w_{a}+h$ for some
small $h$, use the geometric series to expand in powers of $h$ and take the
limit $h\rightarrow 0$. For form factors of derivatives, first take the
derivative of Slavnov’s theorem with respect to one of the free $\\{u\\}$
before performing the same limiting procedure. In all cases, terms
proportional to negative powers of $h$ will vanish identically, as they are
determinants with repeated columns, while those with positive powers of $h$
vanish in the limit. The remaining desired terms are the result. Two
expressions are presented for each form factor: the first as a limit of
Slavnov’s theorem, while the second is a consistency check based on the
physical basis form factors, $F^{i}_{a}$ and $F^{ij}_{ab}$, which we also list
here for convenience. Each form factor has been verified numerically with the
two expressions.

The physical basis form factors are:

$\displaystyle
F^{i}_{a}=\frac{\prod_{c}(v_{c}-\varepsilon_{i})}{\prod_{c}(v_{c}-w_{a})}\frac{\prod_{b\neq
a}(w_{a}-w_{b})}{\prod_{b\neq a}(\varepsilon_{i}-w_{b})}K\det J^{i}_{a}$ (85)

$\displaystyle
F^{ij}_{ab}=\frac{(w_{a}-w_{b})}{(\varepsilon_{i}-\varepsilon_{j})}\frac{\prod_{c}(v_{d}-\varepsilon_{i})(v_{d}-\varepsilon_{j})}{\prod_{c}(v_{d}-w_{a})(v_{d}-w_{b})}\frac{\prod_{d\neq
a,b}(w_{a}-w_{d})(w_{b}-w_{d})}{\prod_{d\neq
a,b}(\varepsilon_{i}-w_{d})(\varepsilon_{j}-w_{d})}K\det J^{ij}_{ab}$ (86)

In the Gaudin basis, the form factors required to compute transition density
matrices are:

$\displaystyle\braket{\\{v\\}}{S^{+}(w_{a})}{\\{w\\}_{c}}$
$\displaystyle=\frac{\prod_{k}(v_{k}-w_{a})}{\prod_{k}(v_{k}-w_{c})}\frac{\prod_{l\neq
c}(w_{c}-w_{l})}{\prod_{l\neq a,c}(w_{a}-w_{l})}K\det J^{a}_{c}$ (87)
$\displaystyle=\sum_{i}\frac{F^{i}_{c}}{w_{a}-\varepsilon_{i}}$ (88)

$\displaystyle\braket{\\{v\\}}{\frac{\partial S^{+}(w_{a})}{\partial
w_{a}}}{\\{w\\}_{a}}$ $\displaystyle=K\det J^{a}_{a}$ (89)
$\displaystyle=-\sum_{i}\frac{F^{i}_{a}}{(w_{a}-\varepsilon_{i})^{2}}$ (90)

$\displaystyle\braket{\\{v\\}}{S^{+}(w_{a})S^{+}(w_{a})}{\\{w\\}_{c,d}}$
$\displaystyle=\frac{\prod_{k}(v_{k}-w_{a})^{2}\prod_{l\neq
c,d}(w_{c}-w_{l})(w_{d}-w_{l})}{\prod_{k}(v_{k}-w_{c})(v_{k}-w_{d})\prod_{l\neq
a,c,d}(w_{a}-w_{l})^{2}}(w_{c}-w_{d})\frac{}{}K\det J^{\bar{a}a}_{cd}$ (91)
$\displaystyle=\sum_{ij}\frac{F^{ij}_{cd}}{(w_{a}-\varepsilon_{i})(w_{a}-\varepsilon_{j})}$
(92)

$\displaystyle\braket{\\{v\\}}{S^{+}(w_{a})S^{+}(w_{b})}{\\{w\\}_{c,d}}$
$\displaystyle=\frac{\prod_{k}(v_{k}-w_{a})(v_{k}-w_{b})\prod_{l\neq
c,d}(w_{c}-w_{l})(w_{d}-w_{l})}{\prod_{k}(v_{k}-w_{c})(v_{k}-w_{d})\prod_{l\neq
a,c,d}(w_{a}-w_{l})\prod_{m\neq b,c,d}(w_{b}-w_{m})}$
$\displaystyle\frac{(w_{c}-w_{d})}{(w_{a}-w_{b})}K\det J^{ab}_{cd}$ (93)
$\displaystyle=\sum_{ij}\frac{F^{ij}_{cd}}{(w_{a}-\varepsilon_{i})(w_{b}-\varepsilon_{j})}$
(94)

$\displaystyle\braket{\\{v\\}}{\frac{\partial S^{+}(w_{a})}{\partial
w_{a}}}{\\{w\\}_{c}}$ $\displaystyle=\frac{\prod_{k}(v_{k}-w_{a})\prod_{l\neq
c}(w_{c}-w_{l})}{\prod_{k}(v_{k}-w_{c})\prod_{l\neq
a,c}(w_{a}-w_{l})}K\left(\det
J^{\bar{a}}_{c}+\left(k(w_{a})+\frac{1}{(w_{a}-w_{c})}\right)\det
J^{a}_{c}\right)$ (95)
$\displaystyle=-\sum_{i}\frac{F^{i}_{c}}{(w_{a}-\varepsilon_{i})^{2}}$ (96)

$\displaystyle\braket{\\{v\\}}{\frac{\partial S^{+}(w_{a})}{\partial
w_{a}}S^{+}(w_{b})}{\\{w\\}_{a,c}}$
$\displaystyle=\frac{\prod_{k}(v_{k}-w_{b})\prod_{l\neq
c}(w_{c}-w_{l})}{\prod_{k}(v_{k}-w_{c})\prod_{l\neq b,c}(w_{b}-w_{l})}K$
$\displaystyle\left(\det
J^{ab}_{ac}+\left(k(w_{a})+\frac{(w_{c}-w_{b})}{(w_{a}-w_{c})(w_{a}-w_{b})}\right)\det
J^{b}_{c}\right)$ (97)
$\displaystyle=-\sum_{ij}\frac{F^{ij}_{ac}}{(w_{a}-\varepsilon_{i})^{2}(w_{b}-\varepsilon_{j})}$
(98)

$\displaystyle\braket{\\{v\\}}{\frac{\partial S^{+}(w_{a})}{\partial
w_{a}}\frac{\partial S^{+}(w_{b})}{\partial w_{b}}}{\\{w\\}_{a,b}}$
$\displaystyle=K\det J^{ab}_{ab}+Kk(w_{a})\det J^{b}_{b}+Kk(w_{b})\det
J^{a}_{a}$ (99)
$\displaystyle=\sum_{ij}\frac{F^{ij}_{ab}}{(w_{a}-\varepsilon_{i})^{2}(w_{b}-\varepsilon_{j})^{2}}$
(100)

$\displaystyle\frac{1}{2}\braket{\\{v\\}}{\frac{\partial^{2}S^{+}(w_{a})}{\partial
w^{2}_{a}}}{\\{w\\}_{a}}$ $\displaystyle=Kk(w_{a})\det J^{a}_{a}+K\det
J^{\bar{a}}_{a}$ (101)
$\displaystyle=\sum_{i}\frac{F^{i}_{a}}{(w_{a}-\varepsilon_{i})^{3}}$ (102)

## References

  * [1] A. C. Hurley, J. Lennard-Jones, and J. A. Pople.  Proceedings of the Royal Society, A220:446, 1953.
  * [2] D. M. Silver.  The Journal of Chemical Physics, 50:5108, 1969.
  * [3] D. M. Silver, E. L. Mehler, and K. Ruedenberg.  The Journal of Chemical Physics, 52:1174, 1970.
  * [4] A. J. Coleman.  International Journal of Quantum Chemistry, 63:23, 1997.
  * [5] P. R. Surján.  An Introduction to the Theory of Geminals.  Springer, Berlin, 1999.
  * [6] M. Kobayashi, A. Szabados, H. Nakai, and P. Surjan.  Journal of Chemical Theory and Computation, 6:2024, 2010.
  * [7] P. R. Surján, Á Szabados, P. Jeszenski, and T. Zoboki.  Journal of Mathematical Chemistry, 50:534, 2012.
  * [8] E. Neuscamman.  Physical Review Letters, 109:203001, 2012.
  * [9] P. A. Limacher, P. W. Ayers, P. A. Johnson, S. De Baerdemacker, D. Van Neck, and P. Bultinck.  Journal of Chemical Theory and Computation, 9:1394, 2013.
  * [10] T. Stein, T. M. Henderson, and G. E. Scuseria.  The Journal of Chemical Physics, 140:214113, 2014.
  * [11] K. Boguslawski, P. Tecmer, P. W. Ayers, P. Bultinck, S. De Baerdemacker, and D. Van Neck.  Physical Review B, 98:201106(R), 2014.
  * [12] K. Boguslawski, P. Tecmer, P. Bultinck, S. De Baerdemacker, D. Van Neck, and P. W. Ayers.  Journal of Chemical Theory and Computation, 10:4873, 2014.
  * [13] K. Boguslawski, P. Tecmer, P. A. Limacher, P. A. Johnson, P. W. Ayers, P. Bultinck, S. De Baerdemacker, and D. Van Neck.  Journal of Chemical Theory and Computation, 140:214114, 2014.
  * [14] P. Tecmer, K. Boguslawski, P. A. Johnson, M. Chan, T. Verstraelen, and P. W. Ayers.  Journal of Physical Chemistry, A118:9058, 2014.
  * [15] T. M. Henderson, G. E. Scuseria, J. Dukelsky, A. Signoracci, and T. Duguet.  Physical Review C, 89:054305, 2014.
  * [16] T. M. Henderson, I. W. Bulik, T. Stein, and G. E. Scuseria.  The Journal of Chemical Physics, 141:244104, 2014.
  * [17] J. J. Shepherd, T. M. Henderson, and G. E. Scuseria.  The Journal of Chemical Physics, 144:094112, 2014.
  * [18] I. W. Bulik, T. M. Henderson, and G. E. Scuseria.  Journal of Chemical Theory and Computation, 11:3171, 2015.
  * [19] K. Boguslawski, P. Tecmer, and Ö. Legeza.  Physical Review B, 94:155126, 2016.
  * [20] K. Boguslawski and P. Tecmer.  Journal of Chemical Theory and Computation, 13:5966, 2017.
  * [21] Zs. É. Mihálka, P. Surján, and Á. Szabados.  Journal of Chemical Theory and Computation, 16:892, 2019.
  * [22] N. Vu and A. E. DePrince III.  The Journal of Chemical Physics, 152:244103, 2020.
  * [23] A. Nowak, Ö. Legeza, and K. Boguslawski.  arxiv:2010.01934.
  * [24] T. M. Henderson and G. E. Scuseria.  The Journal of Chemical Physics, 151:051101, 2019.
  * [25] A. Khamoshi, T. M. Henderson, and G. E. Scuseria.  The Journal of Chemical Physics, 151:184103, 2019.
  * [26] T. M. Henderson and G. E. Scuseria.  The Journal of Chemical Physics, 153:084111, 2020.
  * [27] R. Dutta, T. M. Henderson, and G. E. Scuseria.  Journal of Chemical Theory and Computation, 16:6358, 2020.
  * [28] G. Harsha, T. M. Henderson, and G. E. Scuseria.  The Journal of Chemical Physics, 153:124115, 2020.
  * [29] A. Khamoshi, F. A. Evangelista, and G. E. Scuseria.  Quantum Science and Technology, 6:014004, 2020.
  * [30] A. Khamoshi, G. P. Chen, T. M. Henderson, and G. E. Scuseria.  arXiv:2012.02377.
  * [31] J. Bardeen, Cooper L. N., and J. R. Schrieffer.  Physical Review, 106:162, 1957.
  * [32] J. Bardeen, L. N. Cooper, and J. R. Schrieffer.  Physical Review, 108:1175, 1957.
  * [33] R. W. Richardson.  Physics Letters, 3:277, 1963.
  * [34] R. W. Richardson and N. Sherman.  Nuclear Physics, 52:221, 1964.
  * [35] R. W. Richardson.  Journal of Mathematical Physics, 6:1034, 1965.
  * [36] M. Gaudin.  Journal de Physique, 37:1087, 1976.
  * [37] P. A. Johnson, C.-É. Fecteau, F. Berthiaume, S. Cloutier, L. Carrier, M. Gratton, P. Bultinck, S. De Baerdemacker, D. Van Neck, P. Limacher, and P. W. Ayers.  The Journal of Chemical Physics, 153:104110, 2020.
  * [38] C.-É. Fecteau, H. Fortin, S. Cloutier, and P. A. Johnson.  The Journal of Chemical Physics, 153:164117, 2020.
  * [39] C.-É. Fecteau, F. Berthiaume, M. Khalfoun, and P. A. Johnson.  Journal of Mathematical Chemistry, 2020.
  * [40] P. W. Claeys, J.-S. Caux, D. Van Neck, and S. De Baerdemacker.  Physical Review B, 96:155149, 2017.
  * [41] S. De Baerdemacker, P. W. Claeys, J.-S. Caux, D. Van Neck, and P. W. Ayers.  arXiv:1712.01673.
  * [42] E. K. Sklyanin.  Journal of Soviet Mathematics, 47:2473, 1989.
  * [43] G. Ortiz, R. Somma, J. Dukelsky, and S. Rombouts.  Nuclear Physics B, 707:421, 2005.
  * [44] S. Schechter.  Mathematical Tables and Other Aids to Computation, 13:73, 1959.
  * [45] S. Rombouts, D. Van Neck, and J. Dukelsky.  Physical Review C, 69:061303(R), 2004.
  * [46] A. Faribault, O. El Araby, C. Sträter, and V. Gritsev.  Physical Review B, 83:235124, 2011.
  * [47] O. El Araby, V. Gritsev, and A. Faribault.  Physical Review B, 85:115130, 2012.
  * [48] X. Guan, K. D. Launey, M. Xie, L. Bao, F. Pan, and J. P. Draayer.  Physical Review C, 86:024313, 2012.
  * [49] W. V. Pogosov.  Journal of Physics: Condensed Matter, 24:075701, 2012.
  * [50] S. De Baerdemacker.  Physical Review C, 86:044332, 2012.
  * [51] P. W. Claeys, S. De Baerdemacker, M. Van Raemdonck, and D. Van Neck.  Physical Review B, 91:155102, 2015.
  * [52] L. Carrier, C.-É. Fecteau, and P. A. Johnson.  International Journal of Quantum Chemistry, 120:e26255, 2020.
  * [53] N. A. Slavnov.  Theoretical and Mathematical Physics, 79:502, 1989.
  * [54] S. Belliard and N. A. Slavnov.  Journal of High Energy Physics, 2019:103, 2019.
  * [55] H.-Q. Zhou, J. Links, R. H. McKenzie, and M. D. Gould.  Physical Review B, 65:060502, 2002.
  * [56] G. Gorohovsky and E. Bettelheim.  Physical Review B, 84:224503, 2011.
  * [57] R. Vein and P. Dale.  Determinants and Their Applications in Mathematical Physics.  Springer-Verlag, New York, 1999.


SHG_Supplementray_230513.pdf

††thanks: Deceased

# Spontaneous breaking of mirror symmetry beyond critical doping in Pb-Bi2212

Saegyeol Jung  Center for Correlated Electron Systems, Institute for Basic
Science, Seoul, 08826, Korea  Department of Physics and Astronomy, Seoul
National University, Seoul, 08826, Korea  Byeongjun Seok  Center for
Correlated Electron Systems, Institute for Basic Science, Seoul, 08826, Korea
Department of Physics and Astronomy, Seoul National University, Seoul, 08826,
Korea  Chang jae Roh  Center for Correlated Electron Systems, Institute for
Basic Science, Seoul, 08826, Korea  Department of Physics and Astronomy, Seoul
National University, Seoul, 08826, Korea  Donghan Kim  Center for Correlated
Electron Systems, Institute for Basic Science, Seoul, 08826, Korea  Department
of Physics and Astronomy, Seoul National University, Seoul, 08826, Korea
Yeonjae Lee  Center for Correlated Electron Systems, Institute for Basic
Science, Seoul, 08826, Korea  Department of Physics and Astronomy, Seoul
National University, Seoul, 08826, Korea  San Kang  Center for Correlated
Electron Systems, Institute for Basic Science, Seoul, 08826, Korea  Department
of Physics and Astronomy, Seoul National University, Seoul, 08826, Korea
Shigeyuki Ishida  National Institute of Advanced Industrial Science and
Technology (AIST), Tsukuba, Ibaraki 305-8568, Japan  Shik Shin  Institute for
Solid State Physics (ISSP), The University of Tokyo, Kashiwa, Chiba 277-8581,
Japan  Hiroshi Eisaki  National Institute of Advanced Industrial Science and
Technology (AIST), Tsukuba, Ibaraki 305-8568, Japan  Tae Won Noh  Center for
Correlated Electron Systems, Institute for Basic Science, Seoul, 08826, Korea
Department of Physics and Astronomy, Seoul National University, Seoul, 08826,
Korea  Dongjoon Song  Corresponding autor<EMAIL_ADDRESS>Center for
Correlated Electron Systems, Institute for Basic Science, Seoul, 08826, Korea
Stewart Blusson Quantum Matter Institute, The University of British Columbia,
2355 East Mall $|$ Vancouver BC $|$ V6T 1Z4 Canada  Changyoung Kim
Corresponding autor<EMAIL_ADDRESS>Center for Correlated Electron
Systems, Institute for Basic Science, Seoul, 08826, Korea  Department of
Physics and Astronomy, Seoul National University, Seoul, 08826, Korea

###### Abstract

Identifying ordered phases and their underlying symmetries is the first and
most important step toward understanding the mechanism of high-temperature
superconductivity; critical behaviors of ordered phases are expected to be
correlated with superconductivity. Efforts to find such ordered phases have
been focused on symmetry breaking in the pseudogap region while the Fermi
liquid-like metal region beyond the so-called critical doping $p_{c}$ has been
regarded as a trivial disordered state. Here, we used rotational anisotropy
second harmonic generation and uncovered a broken mirror symmetry in the Fermi
liquid-like phase in (Bi,Pb)2Sr2CaCu2O8+δ with $p=0.205>p_{c}$. By tracking
the temperature evolution of the symmetry-breaking response, we verify an
order parameter-like behavior with the onset temperature $T_{up}$ at which the
strange metal to Fermi liquid-like-metal crossover takes place. Complementary
angle-resolved photoemission study showed that the quasiparticle coherence
between $\mathrm{CuO_{2}}$ bilayers is enhanced in proportion to the symmetry-
breaking response as a function of temperature, indicating that the change in
metallicity and symmetry breaking are linked. These observations contradict
the conventional quantum disordered scenario for over-critical-doped cuprates
and provide new insight into the nature of the quantum critical point in
cuprates.

##  I 1\. Introduction

A long-standing question in the study of cuprate superconductors is how the
normal-state properties are related to the high-transition temperature
(high-$T_{c}$) superconductivity [1, 2]. One of the most exotic normal-state
characteristics is a strange metal behavior with linear temperature dependence
of the resistivity $\rho\sim T^{n}(n=1)$ [3]. In the experimentally
established doping temperature phase diagram of cuprates (Fig. 1a), the
strange metal phase with n = 1 is located above the regions of underdoped
pseudogap and overdoped Fermi liquid (FL)-like metal phases with $T$-sublinear
($n<1$ at $T<T^{*}$) and $T$-superlinear ($n>1$ at $T<T_{up}$) resistivities,
respectively [4, 5] (Fig. 1b , Supplementary Fig. 1). The prevailing view is
that the $T$-linear resistivity persists down to the zero temperature near the
critical doping $p_{c}$ of the underlying pseudogap-FL phase boundary,
resulting in a V-shaped strange metallic region from the putative pseudogap
quantum critical point (QCP) [6, 7]. Therefore, it has been suggested that the
strange metal behavior arises from the competition between the quantum and
thermal fluctuations, i.e., the so-called “QCP scenario” [8].

In the generic QCP scenario, the under-critical-doped pseudogap and over-
critical-doped FL-like metal phases play the roles of the quantum ordered and
disordered states, respectively, indicating that the transition from the
strange metal to the pseudogap phase at $T^{*}$ is a phase transition with a
corresponding symmetry breaking, while that to the FL-like metal phase at
$T_{up}$ is a smooth crossover. Indeed, various symmetry breaking across the
pseudogap phase boundary has been observed, which is consistent with the QCP
scenario [9, 10, 11, 12]. On the other hand, the over-critical-doped FL-like
metal region has been investigated less systematically, and it is unclear
whether it definitively represents a quantum disordered state. In contrast to
the long-standing belief, ordered phenomena, such as charge order and
ferromagnetism, have been observed in the heavily overdoped region [13, 14,
15]. The discovery of ordered phases is surprising not only because it is
incompatible with the QCP scenario, but also because the collective bosonic
mode attributed to the order parameter fluctuations can mediate the fermion-
fermion interaction, as well as superconductivity. Yet, it still remains an
open question whether the FL-like metal phase of cuprates is genuinely related
to a symmetry broken state with a hidden order or not.

To detect subtle symmetry breaking, rotational anisotropy second harmonic
generation (RA-SHG) has been used to research various quantum materials [16,
17, 18]. Specifically, as it measures the high-rank (n$>$2) optical
susceptibility tensor, it is useful in finding both electronically and
magnetically driven symmetry breaking induced by multipolar orders. Indeed,
previous RA-SHG studies on undoped Sr2CuO2Cl2 and under-critical-doped
YBa2Cu3Oy revealed hidden symmetry breaking inside antiferromagnetic and
pseudogap regime, respectively [19, 20]. To explore the symmetry evolution in
the over-critical-doped region using RA-SHG, (Bi,Pb)2Sr2CaCu2O8+δ (Pb-Bi2212)
is suitable because its crystallographic structure has higher symmetry
compared to other cuprates with monoclinic structures, such as LSCO, YBCO, and
pristine Bi2212 [19, 21]. For example, Fig. 1e shows RA-SHG intensity
$I_{SS}^{2\omega}$ of Pb-($T_{c}=79$ K) and pristine Bi2212 ($T_{c}=92$ K) as
a function of the $\varphi$, angle of the incident plane rotated around the
$c$-axis, measured in the $\mathrm{S_{in}{\text{-}}S_{out}}$ (SS) geometry at
room temperature (Fig. 1d). While mirror symmetry along the crystal axis is
broken in pristine Bi2212 due to monoclinic distortion [22], the SHG data of
Pb-Bi2212 is quite isotropic, which is beneficial for identifying the
emergence of small anisotropy in the SHG response.

Here, we performed a polarization- and temperature-dependent RA-SHG study,
along with complementary angle-resolved photoemission spectroscopy (ARPES)
measurements of the over-critical-doped Pb-Bi2212 with $p\sim 0.205$ and
$T_{up}\sim 225$ K. Our SHG investigation revealed the unexpected discovery of
mirror symmetry breaking with onset temperature coincident with the strange
metal to FL-like metal crossover. Moreover, we found that the symmetry
breaking is accompanied by order parameter-like enhancement of not only the
SHG response, but also the photoemission quasiparticle coherence with the same
transition temperature $T_{up}$. Our observations suggest that the ground
state of the FL-like metal in the $p>p_{c}$ region is a quantum ordered state
with broken mirror symmetry, which contradicts the conventional QCP scenario.

Figure 1: Orthorhombic crystallographic structure of Pb-Bi2212 with $p>p_{c}$.
(a) Temperature versus doping phase diagram of pristine Bi2212 extracted from
the temperature deviation from the linear behavior of in-plane resistivity
$\rho_{ab}$(T) (green diamonds [23], red triangles [24], and blue squares
[25]). The doping measured in this study is indicated by the grey line. (b)
Temperature-dependent resistivity curve for Pb-Bi2212 with $p\sim 0.205$. The
grey circle marks $T_{up}$, the onset of upturn deviation from linear
resistivity. (c) Crystal structure of Pb-Bi2212. Pb substitution alleviates
the structural modulation that exists along the $a$-axis. (d) Schematic
illustration of the RA-SHG setup. Linear polarized light, which is parallel
(P) or perpendicular (S) to the incidence plane, with frequency $\omega$
focused obliquely onto the (001) surface of the sample. $\varphi$ is the angle
of incident plane rotated around the $c$-axis. (e) RA-SHG patterns of Pb-
Bi2212 and pristine Bi2212 acquired in $\mathrm{S_{in}{\text{-}}S_{out}}$ (SS)
polarization geometry at room temperature. Data were fitted to the bulk
electric quadrupole induced SHG from $mmm$ and $2/m$, respectively. Figure 2:
Spontaneous mirror symmetry breaking in Pb-Bi2212. RA-SHG data of Pb-Bi2212
collected under four polarization geometries: (a,e) PP, (b,f) SP, (c,g) PS,
and (d,h) SS at $T\sim 295$ K ($T>T_{up}$) (a–d) and $T\sim 70$ K ($T<T_{up}$)
(e–h). All data sets are presented with the same intensity scale, in which the
maximum intensity in the PP geometry at $T\sim 295$ K is set to 1. A data set
of room and low temperatures for each polarization was obtained from same
samples (see Methods for details). The data taken at room temperature are
fitted to bulk EQ-induced SHG from $mmm$ (red lines), whereas the low-
temperature data are fitted to the coherent superposition of bulk EQ-induced
SHG from $mmm$ and $mm^{\prime}m^{\prime}$ (blue lines), as described in the
main text.

##  II 2\. Result

Fig. 2 shows experimental and simulated RA-SHG results (open circles and solid
lines, respectively) obtained with four polarization geometries at room
temperature and near 70 K: $\mathrm{P_{in}{\text{-}}P_{out}}$ (PP),
$\mathrm{P_{in}{\text{-}}S_{out}}$ (PS), $\mathrm{S_{in}{\text{-}}P_{out}}$
(SP), and $\mathrm{S_{in}{\text{-}}S_{out}}$ (SS). As there is no consensus on
the point group of Bi-based cuprates, we first determined the crystallographic
point group of Pb-Bi2212 at room temperature [26]. For all of the RA-SHG polar
data at room temperature, two mirror symmetries along the crystal axes $a$ and
$b$ (two mirror planes of $m_{ac}$ and $m_{bc}$) are clearly identified. We
simulated RA patterns from bulk crystallographic point group candidates within
these two mirror symmetries and twofold rotational symmetry: centrosymmetric
$mmm$ and noncentrosymmetric $mm2$. We found that bulk electric quadrupole
(EQ, $Q_{ij}^{2\omega}\sim\chi_{ijkl}^{EQ}E_{k}E_{l})$) SHG derived from the
$mmm$ point group $I_{SS}^{2\omega}\propto
sin^{2}2\varphi|\mathrm{A}~{}sin^{2}\varphi+\mathrm{B}~{}cos^{2}\varphi|^{2}$
well fit eight lobes in the SS channel, where A and B are linear combinations
of $\chi_{ijkl}^{EQ}$ susceptibility tensor components. Moreover, the fitting
results in Figs. 2a–d indicate excellent agreement in the other polarization
channels with an EQ SHG response from the centrosymmetric $mmm$ orthorhombic
point group. In contrast, we excluded $mm2$ from the point group candidates
because the electric dipole (ED,
$P_{i}^{2\omega}\sim\chi_{ijk}^{ED}E_{j}E_{k})$ SHG, which is the dominant
contribution from noncentrosymmetric $mm2$, is forbidden in the SS
polarization channel data (Methods and Supplementary Section 2).

Next, we turned our attention to the temperature-dependent symmetry change.
Figs. 2e–h show the RA-SHG at low temperature $T\sim$ 70 K below $T_{up}$ as
well as $T_{c}$. Compared to the room temperature results, PP and SP channels
show larger SHG intensities while retaining the symmetries. In contrast to PP
and SP cases, PS and SS channels exhibit significant amplitude modulation of
the lobes across the crystal axes, which is not observed at room temperature.
This constitutes clear evidence that the mirror symmetry of the system is
broken at low temperature, while the onset temperature of this symmetry
breaking remains questionable. We note that the signature of the mirror
symmetry breaking is detected only in the polarization channels with S output,
at which room temperature data are mostly described by in-plane tensor
elements ($\chi_{ijkl}^{EQ}$ elements with $i,j,k,l$ = $x$ or $y$). This
implies that the low-temperature SHG response mainly comes from the in-plane
symmetry breaking.

To find the onset temperature of symmetry breaking, we fixed $\varphi$ to the
maximum lobe intensity angle and tracked the temperature dependence of the SHG
intensity $I^{2\omega}(T)$ in fine temperature steps. We obtained the $\Delta
I^{2\omega}(T)$ by subtracting the linear-in temperature background from the
$I^{2\omega}(T)$, where the linear-in temperature background is attributed to
the linearly decreasing change in the crystallographic $\chi_{ijkl}^{EQ}$
elements [27, 28, 20] (Supplementary Fig. 2). Figs. 3a–d show $\Delta
I^{2\omega}(T)$ for all four geometries normalized to their room temperature
values before subtracting the background (Supplementary Fig. 3). With changing
temperature, the SHG intensity in all polarization geometries abruptly
increases below $T\sim 225$ K which is close to $T_{up}$ observed in the
resistivity data instead of $T_{c}$. This feature is a strong indicator of
spontaneous symmetry breaking across the border from a strange metal to FL-
like metal phase, in contrast to the long-standing conjecture that $T_{up}$ is
associated with a smooth crossover to the disordered state. Moreover, except
for the $\mathrm{S_{in}{\text{-}}P_{out}}$ data, $\Delta I^{2\omega}(T)$ shows
obvious order parameter-like behavior below $T_{up}$, indicating that there is
an order parameter lowering the symmetry of the system.

An important open question is whether the $T_{up}$ coincidently meets the
onset temperature of the symmetry breaking, or whether they are mutually
connected. As $T_{up}$ is an indication of a change in dynamics of charge
carriers, the possible emergence of the order parameter is expected to impact
the quasiparticle self-energy, which is reflected in the photoemission
spectral function. Therefore, we performed complementary ARPES measurements
with fine temperature steps to observe the effects of symmetry breaking on the
quasiparticles. Figs. 4a–d show ARPES spectra and the second energy
derivatives, respectively, above $T_{up}$ at 280 K and below $T_{up}$ at $140$
K, along the Brillouin zone boundary. Compared to the blurred spectral
distribution in the $280$ K data, the spectrum at 140 K showed clear bilayer
splitting with certain features of the bonding and anti-bonding band. We
traced the temperature evolution of the split peaks through the energy
distribution curve (EDC) at the anti-nodal point (Fig. 4e, inset). In Fig. 4e,
while EDC near room temperature shows a shoulder-like feature around -50 meV,
the split peaks become clearly discernible with decreasing temperature, which
implies that the doped holes on the two $\mathrm{CuO_{2}}$ bilayers in the
unit cell propagate more coherently at lower temperature.

To quantify the quasiparticle coherence and track its temperature evolution,
we focus on the temperature dependence of the spectral weight transfer to the
coherent split peaks, $W_{coh}(T)$ (Inset of Fig. 4f). In detail, the $\Delta
W_{coh}(T)$ is obtained by integrating the spectral weight difference between
EDC of 300 K and each temperature for the energy range of -75
meV$<E{\text{-}}E_{F}<$0 with subsequent subtraction of the linear-in
temperature background, which is considered to be the simple impurity
scattering factor (Supplementary Fig. 3). In Fig. 4f, the $\Delta W_{coh}$ is
plotted along with $\rho-\rho_{Linear}$ and $\Delta I_{SS}^{2\omega}$, where
$\rho_{Linear}$ is the linear component of resistivity obtained by fitting the
$\rho(T)$ above $T_{up}$ (dashed line in Fig. 1b); thus, $\rho-\rho_{Linear}$
is a measure of the superlinear component. Remarkably, the $\Delta W_{coh}$
and $\Delta I_{SS}^{2\omega}$ overlap with each other, showing almost the same
order parameter-like temperature dependence with the transition temperature
$T$ = $T_{up}$, where the superlinear component of resistivity emerges.

A previous ARPES study proposed an association between the sudden enhancement
of quasiparticle coherence and the strange metal to FL-like metal crossover,
while the FL-like metal was still assumed to be disordered [25]. In addition
to this long-standing view, our precise quantification of the spectral weight
transfer and comparative study with SHG show that the sudden enhancement of
coherence is intertwined with the order parameter of the broken symmetry.
Therefore, the results presented here suggest that the strange metal to FL-
like metal transition at $T_{up}$ in Pb-Bi2212 is a phase transition
accompanied by mirror symmetry breaking.

Figure 3: Temperature dependent RA-SHG pattern amplitude for (a) PP, (b) SP,
(c) PS, and (d) SS in a heating cycle. Data were obtained with $\varphi$ fixed
to the maximum SHG lobe at low temperature. Data were first normalized to the
values at $T\sim 295$ K. Then, a high-temperature linear background was
subtracted. The error bars represent the standard deviation of the intensity
over 12 independent measurements in the PP, SP, and PS channels. The error
bars in the SS channel denote the standard deviation of 60 measurements while
heating by 5 K with a ramping rate of 2 K/min (Methods). Blue lines on the
data are guides for the eye. The width of the shaded grey interval shows the
uncertainty in the transition temperature.  Figure 4: Temperature dependence
of electronic structure in Pb-Bi2212 across $\bm{T_{up}}$. (a–b) ARPES spectra
(left) and second energy derivatives (right) taken above $T_{up}$ at 280 K and
(c–d) below $T_{up}$ at 140 K, along the Brillouin zone boundary, as denoted
by the solid grey arrow in the inset of (e). All ARPES data were divided by
the energy resolution-convoluted Fermi-Dirac distribution function, and the
momentum-independent background was then subtracted (see Methods for details).
(e) Temperature dependence of EDC at the anti-node (grey dot in the inset of
(e)). Each EDC was normalized to the total area over the energy range -200
meV$<E{\text{-}}E_{F}<$0\. Black rectangles denote the bilayer splitting of
antibonding band (AB) and bonding band (BB). The inset shows a schematic of
the Fermi surface of Pb-Bi2212 composed of the bilayer splitting. (f)
$\rho-\rho_{Linear}$ as a function of temperature, where $\rho_{Linear}$ is
the linear fitting of $\rho(T)$ between 250 and 300 K (top). Temperature
dependence of the integrated spectral weight difference over the energy range
-75 meV$<E{\text{-}}E_{F}<$0 between 300 K and each temperature ($W_{coh}$,
orange area in the inset) with subtracting linear-in temperature background,
in comparison with the normalized SHG intensity in SS geometry from Fig. 3d
(bottom).

##  III 4\. Discussion

We attempted to narrow down the point group symmetry candidates of the FL-like
phase by surveying both crystallographic and magnetic subgroups of the $mmm$
point group, and categorized possible sources of the low-temperature RA-SHG
response to the ED contribution from noncentrosymmetric subgroups, and the MD
or EQ contribution from centrosymmetric subgroups. Table 1 shows the highest-
rank subgroups of the $mmm$ point group and the possible SHG enhancement in
each polarization geometry. As susceptibility tensors $\chi_{ijk}^{ED}$,
$\chi_{ijk}^{MD}$, and $\chi_{ijkl}^{EQ}$ are invariant under the
corresponding symmetry operation of each subgroup, the SHG response can occur
in certain polarization geometries (Supplementary Section 2). Therefore, all
noncentrosymmetric subgroup candidates can be excluded due to the forbidden ED
SHG enhancement in SS geometry, which contradicts our observation in Fig. 3d.
As the SHG enhancement at $T_{up}$ occurs in all polarization geometries, we
identified the subgroup candidates from Table 1: the EQ process from $2/m$
crystallographic subgroup and EQ process from $mm^{\prime}m^{\prime}$ magnetic
subgroup, where $m^{\prime}$ denotes the combination of mirror operation and
time reversal.

We further fitted the low-temperature polar data with the two subgroup
candidates of $2/m$ and $mm^{\prime}m^{\prime}$, and crosschecked whether the
subgroup candidates indeed reflect the symmetry of the system. For the
crystallographic subgroup $2/m$, all of the data with four different
polarizations are well fitted with
$(P_{i}^{2\omega}\sim\chi_{ijkl}^{EQ}\nabla_{j}E_{k}E_{l})$ where
$\chi_{ijkl}^{EQ}$ is the nonvanishing susceptibility tensor from $2/m$
(Supplementary Fig. 4). In the case of $mm^{\prime}m^{\prime}$, $c$-type
contribution $\chi_{ijkl}^{EQ(c)}$ is allowed because time reversal symmetry
is broken. Therefore, we fit the data with
$P_{i}^{2\omega}\sim(\chi_{ijkl}^{EQ(i)}+\chi_{ijkl}^{EQ(c)})\nabla_{j}E_{k}E_{l}$,
where $\chi_{ijkl}^{EQ(i)}$ is the time-invariant tensor from crystallographic
point group $mmm$ and $\chi_{ijkl}^{EQ(c)}$ is the time-noninvariant tensor
from the magnetic point group $mm^{\prime}m^{\prime}$. The solid lines in
Figs. 2e–h show that a coherent superposition of $\chi_{ijkl}^{EQ(i)}$ from
$mmm$ and $\chi_{ijkl}^{EQ(c)}$ from $mm^{\prime}m^{\prime}$ well reproduces
two characteristic features of low-temperature S-output data: the modulation
of lobe amplitude, which is a marker of mirror symmetry breaking, and
shrinkage along the $a$-axis (Supplementary Fig. 5). These results showed that
$2/m$ and $mm^{\prime}m^{\prime}$ are suitable candidates to explain both the
temperature dependence of SHG intensity and low-temperature polar data.

Based on these two sub point groups, we discuss the potential physical origins
for the order parameter. We first examined the monoclinic structural
transition, which enabled the $mmm$ to evolve into the $2/m$ symmetry [29,
30]. It is worth noting that, based on our ARPES results, the low-temperature
ordering that reduces the symmetry is correlated with enhancement of
quasiparticle coherence along the $c$-axis, at least to the distance between
the $\mathrm{CuO_{2}}$ bilayers. Therefore, if monoclinic distortion is the
origin of the symmetry breaking, it would impact the $c$-axis lattice
constant. However, our complementary X-ray diffraction (XRD) study, which
focused on the d spacing measurement along the $c$-axis, did not identify any
significant anomalies in the temperature evolution of the $c$-axis lattice
constant (Supplementary Fig. 4). Moreover, previous XRD and neutron scattering
experiments with various Bi-based cuprates found no evidence of such
structural transition near $T_{up}$ [31, 32]. Nevertheless, we still cannot
fully rule out unidentified structural transitions. Therefore, further in-
depth study of the structural transition is required.

Another possibility is that the enhancement of SHG at $T_{up}$ could be
attributed to a hitherto unknown bulk order parameter described by the
$mm^{\prime}m^{\prime}$. Indeed, $mm^{\prime}m^{\prime}$ point group
incorporates $A_{2g}$ inversion symmetric magnetic order, which is compatible
with ferroic order, for example. Intriguingly, several fingerprints of
ferromagnetic fluctuations have been reported in the heavily overdoped regime
of various cuprates [13, 33]. In addition, coupling of such magnetic
fluctuations with optical phonons is theoretically able to generate the
ferroic order with a quasistatic magnetoelectric quadrupole [34]. Detection of
such multipolar orders via, for example, magnetic neutron scattering and muon
spin measurements represents an important direction for further study [35].

Lastly, our comprehensive study with SHG, ARPES, and resistivity promotes a
better understanding of the strange metal phase. We highlight that the
$T$-superlinear component of resistivity $\rho-\rho_{Linear}$ becomes zero as
the coherence of the bilayer-split quasiparticle peaks and multipolar order
which can be described by the $z$-axis component of the magnetization are
significantly suppressed above $T_{up}$ (see Fig. 4). This implies that the
decoherence of doped holes between the $\mathrm{CuO_{2}}$ layers and the
resulting two-dimensional confinement is an important ingredient of the
strange metal behavior in the over-critical-doped regime [36]. These results
provide experimental support for the extensive theoretical descriptions of
non-Fermi liquid behaviors from QCP in two-dimensional metals including
cuprates [37, 38].

[0.25cm] Point group | Process | PP | SP | PS | SS |
---|---|---|---|---|---|---
$222$ | ED | O | O | O | – |
$mm2$ | ED | O | O | O | – |
$2/m$ | MD | O | O | – | O |
[0.25cm] $2/m$ |  EQ | O | O | O | O
$mmm^{\prime}$ | ED | O | O | O | – |
$mm^{\prime}m^{\prime}$ | MD | O | O | – | O |
[0.25cm] $mm^{\prime}m^{\prime}$ |  EQ | O | O | O | O
$m^{\prime}m^{\prime}m^{\prime}$ | ED | O | O | O | – |

Table 1: Summary of polarization geometry dependent SHG for the highest-rank
subgroups of the $mmm$ point group. SHG process in light grey background
indicate possible candidates which are consistent with both the temperature
dependence of SHG intensity and low temperature polar data. See Supplementary
section 2 for detailed equations.

##  IV Acknowledgements

We are grateful to Y. S. Kim, R. Noguchi, S. S. Huh, H. Y. Choi and A. Hallas
for their helpful discussions and useful comments. We appreiciate the
technical support on fitting process from B.T. Fichera. This work was
conducted under the ISSP-CCES Collaborative Programme and was supported by the
Institute for Basic Science in Republic of Korea (Grant Numbers IBS-R009-G2
and Grant No. IBSR009-D1 ). This work was supported by the National Research
Foundation of Korea(NRF) grant funded by the Korea government(MSIT) (No.
2022R1A3B1077234) and the JSPS KAKENHI (No. JP19H0582*).

## References

  * [1] Keimer, B., Kivelson, S. A., Norman, M. R., Uchida, S. & Zaanen, J.  From quantum matter to high-temperature superconductivity in copper oxides.  _Nature_ 518, 179–186 (2015).
  * [2] Fradkin, E., Kivelson, S. A. & Tranquada, J. M.  Colloquium: Theory of intertwined orders in high temperature superconductors.  _Reviews of Modern Physics_ 87, 457 (2015).
  * [3] Varma, C. M.  Colloquium: Linear in temperature resistivity and associated mysteries including high temperature superconductivity.  _Reviews of Modern Physics_ 92, 031001 (2020).
  * [4] Naqib, S., Cooper, J., Tallon, J. & Panagopoulos, C.  Temperature dependence of electrical resistivity of high-tc cuprates—-from pseudogap to overdoped regions.  _Physica C: Superconductivity_ 387, 365–372 (2003).
  * [5] Sterpetti, E., Biscaras, J., Erb, A. & Shukla, A.  Comprehensive phase diagram of two-dimensional space charge doped $\mathrm{Bi_{2}Sr_{2}CaCu_{2}O_{8+\delta}}$.  _Nature communications_ 8, 1–8 (2017).
  * [6] Hussey, N. E.  Phenomenology of the normal state in-plane transport properties of high-$t_{c}$ cuprates.  _Journal of Physics: Condensed Matter_ 20, 123201 (2008).
  * [7] Cooper, R. _et al._ Anomalous criticality in the electrical resistivity of $\mathrm{La_{2-x}Sr_{x}CuO_{4}}$.  _Science_ 323, 603–607 (2009).
  * [8] Proust, C. & Taillefer, L.  The remarkable underlying ground states of cuprate superconductors.  _Annual Review of Condensed Matter Physics_ 10, 409–429 (2019).
  * [9] Xia, J. _et al._ Polar kerr-effect measurements of the high-temperature $\mathrm{YBa_{2}Cu_{3}O_{6+x}}$superconductor: Evidence for broken symmetry near the pseudogap temperature.  _Physical Review Letters_ 100, 127002 (2008).
  * [10] Sato, Y. _et al._ Thermodynamic evidence for a nematic phase transition at the onset of the pseudogap in $\mathrm{YBa_{2}Cu_{3}O_{y}}$.  _Nature Physics_ 13, 1074–1078 (2017).
  * [11] Zhang, J. _et al._ Discovery of slow magnetic fluctuations and critical slowing down in the pseudogap phase of $\mathrm{YBa_{2}Cu_{3}O_{y}}$.  _Science advances_ 4, eaao5235 (2018).
  * [12] Ishida, K. _et al._ Divergent nematic susceptibility near the pseudogap critical point in a cuprate superconductor.  _Journal of the Physical Society of Japan_ 89, 064707 (2020).
  * [13] Kurashima, K. _et al._ Development of ferromagnetic fluctuations in heavily overdoped $\mathrm{(Bi,Pb)_{2}Sr_{2}CuO_{6+\delta}}$ copper oxides.  _Physical review letters_ 121, 057002 (2018).
  * [14] Peng, Y. _et al._ Re-entrant charge order in overdoped $\mathrm{(Bi,Pb)_{2.12}Sr_{1.88}CuO_{6+\delta}}$ outside the pseudogap regime.  _Nature materials_ 17, 697–702 (2018).
  * [15] Miao, H. _et al._ Charge density waves in cuprate superconductors beyond the critical doping.  _npj Quantum Materials_ 6, 1–6 (2021).
  * [16] Zhao, L. _et al._ Evidence of an odd-parity hidden order in a spin–orbit coupled correlated iridate.  _Nature Physics_ 12, 32–36 (2016).
  * [17] Harter, J., Zhao, Z., Yan, J.-Q., Mandrus, D. & Hsieh, D.  A parity-breaking electronic nematic phase transition in the spin-orbit coupled metal $\mathrm{Cd_{2}Re_{2}O_{7}}$.  _Science_ 356, 295–299 (2017).
  * [18] Jin, W. _et al._ Observation of a ferro-rotational order coupled with second-order nonlinear optical fields.  _Nature Physics_ 16, 42–46 (2020).
  * [19] Zhao, L. _et al._ A global inversion-symmetry-broken phase inside the pseudogap region of yba2cu3oy.  _Nature Physics_ 13, 250–254 (2017).
  * [20] Torre, A. d. l. _et al._ Mirror symmetry breaking in a model insulating cuprate.  _Nature Physics_ 17, 777–781 (2021).
  * [21] Frison, R. _et al._ Crystal symmetry of stripe-ordered la 1.88 sr 0.12 cuo 4.  _Physical Review B_ 105, 224113 (2022).
  * [22] Kan, X. & Moss, S.  Four-dimensional crystallographic analysis of the incommensurate modulation in a $\mathrm{Bi_{2}Sr_{2}CaCu_{2}O_{8}}$ single crystal.  _Acta Crystallographica Section B: Structural Science_ 48, 122–134 (1992).
  * [23] Usui, T. _et al._ Doping dependencies of onset temperatures for the pseudogap and superconductive fluctuation in $\mathrm{Bi_{2}Sr_{2}CaCu_{2}O_{8+\delta}}$ , studied from both in-plane and out-of-plane magnetoresistance measurements.  _journal of the physical society of japan_ 83, 064713 (2014).
  * [24] Watanabe, T., Fujii, T. & Matsuda, A.  Anisotropic resistivities of precisely oxygen controlled single-crystal $\mathrm{Bi_{2}Sr_{2}CaCu_{2}O_{8+\delta}}$ : Systematic study on“spin gap”effect.  _Physical review letters_ 79, 2113 (1997).
  * [25] Kaminski, A. _et al._ Crossover from coherent to incoherent electronic excitations in the normal state of $\mathrm{Bi_{2}Sr_{2}CaCu_{2}O_{8+\delta}}$.  _Physical review letters_ 90, 207003 (2003).
  * [26] Ivanov, A. A. _et al._ Local noncentrosymmetric structure of $\mathrm{Bi_{2}Sr_{2}CaCu_{2}O_{8+y}}$ by x-ray magnetic circular dichroism at cu k-edge xanes.  _Journal of Superconductivity and Novel Magnetism_ 31, 663–670 (2018).
  * [27] Pavlov, V., Pisarev, R., Kirilyuk, A. & Rasing, T.  Observation of a transversal nonlinear magneto-optical effect in thin magnetic garnet films.  _Physical review letters_ 78, 2004 (1997).
  * [28] Ron, A., Zoghlin, E., Balents, L., Wilson, S. & Hsieh, D.  Dimensional crossover in a layered ferromagnet detected by spin correlation driven distortions.  _Nature communications_ 10, 1–6 (2019).
  * [29] Axe, J. _et al._ Structural phase transformations and superconductivity in $\mathrm{La_{2-x}Ba_{x}CuO_{4}}$.  _Physical review letters_ 62, 2751 (1989).
  * [30] Büchner, B., Breuer, M., Freimuth, A. & Kampf, A. P.  Critical buckling for the disappearance of superconductivity in rare-earth-doped $\mathrm{La_{2-x}Sr_{x}CuO_{4}}$.  _Physical review letters_ 73, 1841 (1994).
  * [31] Yang, Z. _et al._ Thermal expansion of $\mathrm{Bi_{2.2}Sr_{1.8}CaCu_{2}O_{x}}$ superconductor single crystals.  _Journal of superconductivity_ 8, 233–239 (1995).
  * [32] Asahi, T., Suzuki, H., Nakamura, M., Takano, H. & Kobayashi, J.  Thermal expansion of superconducting bi 2 sr 2 cacu 2 o 8.  _Physical Review B_ 55, 9125 (1997).
  * [33] Sarkar, T. _et al._ Ferromagnetic order beyond the superconducting dome in a cuprate superconductor.  _Science_ 368, 532–534 (2020).
  * [34] Fechner, M., Fierz, M. J., Thöle, F., Staub, U. & Spaldin, N. A.  Quasistatic magnetoelectric multipoles as order parameter for pseudogap phase in cuprate superconductors.  _Physical Review B_ 93, 174419 (2016).
  * [35] Lovesey, S., Chatterji, T., Stunault, A., Khalyavin, D. & van der Laan, G.  Direct observation of anapoles by neutron diffraction.  _Physical Review Letters_ 122, 047203 (2019).
  * [36] Boebinger, G. S.  An abnormal normal state.  _Science_ 323, 590–591 (2009).
  * [37] Stockert, O., Löhneysen, H. v., Rosch, A., Pyka, N. & Loewenhaupt, M.  Two-dimensional fluctuations at the quantum-critical point of cecu 6- x au x.  _Physical review letters_ 80, 5627 (1998).
  * [38] Metlitski, M. A. & Sachdev, S.  Quantum phase transitions of metals in two spatial dimensions. ii. spin density wave order.  _Physical Review B_ 82, 075128 (2010).
  * [39] Kim, S. _et al._ A compact and stable incidence-plane-rotating second harmonics detector.  _Review of Scientific Instruments_ 92, 043905 (2021).
  * [40] Fichera, B.  ShgPy.  https://bfichera.github.io/shgpy/ (2020).

##  V Methods

##  VI Material growth

Single crystals of $\mathrm{Bi_{1.7}Pb_{0.5}Sr_{2}CaCu_{2}O_{8+\delta}}$ were
grown by the floating-zone (FZ) method. The $\mathrm{Pb}$ substitution
suppresses the monoclinic superstructure modulation which comes from from the
$\mathrm{BiO_{2}}$ layer. The $T_{c}$ was defined as the onset temperature of
the Meissner effect, which was measured using the DC magnetic susceptibility
module of a Physical Property Measurement System from Quantum Design. Hole
doping $p$ was determined via the empirical parabolic relation. $T_{c}$ =
$T_{c,max}[1-82.6(p-0.16)^{2}]$, where $T_{c,max}$ = 95 K for the batch of
crystals studied here. The crystals used in our experiments were aligned with
high accuracy by X-ray Laue diffraction. $a$-axis, which is the
$\mathrm{BiO_{2}}$ modulation direction, was assigned by measuring the
suppressed superstructure replica of the main electronic band with ARPES
(Supplementary Fig. 6).

##  VII RA-SHG measurements

The RA-SHG measurements were obtained in a rotating incident plane setup based
on Fresnel rhomb prism (FR600QM; Thorlabs) [39]. The light source was a
Ti:Sapphire laser operating with a repetition rate of 250 kHz, pulse width of
150 fs, and center wavelength of 840 nm (RegA 9000; Coherent). The incidence
angle was fixed by specification of Fresnel rhomb prism to $\theta\sim
23^{\circ}$. The incident light was focused to a spot size of $\sim$ 30$\mu$m,
with a fluence of $\sim$ 2 mJ/$cm^{2}$ maintained to avoid laser-induced
heating effects on the sample only except the polar data in Figs. 2 d and h
(Supplementary Fig. 7). SS geometry polar data in Figs. 2d and h were obtained
with a fluence of $\sim$ 3 mJ/$cm^{2}$. Each data set at 295 K and 70 K in
Fig. 2 was obtained from same samples, cleaved in air before measurement and
pumped down to the pressure below $2\times 10^{-6}$ Torr . Additionally, each
data set of different polarization geometry was obtained from different
samples. Reflected SHG light was collected using a photomultiplier tube with a
lock-in amplifier synchronized to the 2-kHz pulsed output, with a time
constant of 3 s. Further details of the experimental setup are described in
Ref [39]. We carefully aligned the setup to retain the same incident angle
while rotating the incident plane, and confirmed the parallelity of incident
polarization and output polarization with a gold mirror and GaAs (110),
respectively (Supplementary Fig. 8). When measuring the temperature dependence
of SHG with SS geometry in Fig. 3d, we continuously measured the SHG intensity
while increasing the temperature with a ramping rate of 2 K/min to overcome
the weakest intensity in SS geometry. For the other polarization geometries,
PP, SP, and PS, we measured SHG intensity after stabilizing each temperature.
The data in this article were not only reproduced by scanning along the sample
surface, but also across multiple single crystals grown in different batches
(Supplementary Fig. 9).

##  VIII Fitting procedure.

For fitting the RA-SHG data, we used the open source software ShgPy provided
by Brian Fichera [40]. The high- ($T>T_{up}$) and low- ($T<T_{up}$)
temperature SHG data were fitted to the expression,
$I^{2\omega}\sim|\chi_{ijkl}^{EQ(i)}\nabla_{j}E_{k}E_{l}|^{2}$,
$I^{2\omega}\sim|(\chi_{ijkl}^{EQ(i)}+\chi_{ijkl}^{EQ(c)})\nabla_{j}E_{k}E_{l}|^{2}$,
where $\chi_{ijkl}^{EQ(i)}$ is the time-invariant tensor from crystallographic
point group $mmm$, and $\chi_{ijkl}^{EQ(c)}$ is the time-noninvariant tensor
from magnetic point group $mm^{\prime}m^{\prime}$. RA-SHG data in all
polarization geometries (PP, SP, PS, SS) were fitted simultaneously with 15
nonzero independent tensor elements from $mmm$ ($\chi_{xxxx}$, $\chi_{xxyy}$,
$\chi_{xxzz}$, $\chi_{xyyx}$, $\chi_{xzzx}$, $\chi_{yxyx}$, $\chi_{yyxx}$,
$\chi_{yyyy}$, $\chi_{yyzz}$, $\chi_{yzzy}$, $\chi_{zxzx}$, $\chi_{zyzy}$,
$\chi_{zzxx}$, $\chi_{zzyy}$, and $\chi_{zzzz}$) and 13 nonzero independent
tensor elements from $mm^{\prime}m^{\prime}$ ($\chi_{xxyx}$, $\chi_{xyxx}$,
$\chi_{xyyy}$, $\chi_{xyzz}$, $\chi_{xzzy}$, $\chi_{yxxx}$, $\chi_{yxyy}$,
$\chi_{yxzz}$, $\chi_{yyyx}$, $\chi_{yzzx}$, $\chi_{zxzy}$, $\chi_{zyzx}$, and
$\chi_{zzyx}$). Here, we assumed that all tensor components were real to
reduce the complexity of the computation. The explicit fitting functions are
summarized in Supplementary Section 2. Note that the fitting parameter values
are not unique solutions of the fitting equations. Therefore, the fitting
results only indicate whether a specific point group can reproduce four RA-SHG
polar data set.

##  IX ARPES measurements.

ARPES measurements were performed with a lab-based system equipped with a
He-$I\alpha$ discharge lamp ($hv$=21.2 eV) and Scienta DA30 electron analyzer.
The energy resolution of the system was set to 12 meV . For the data used in
Fig. 4f, the sample was cleaved at 300 K with a chamber pressure around
$9\times 10^{-11}$ Torr and measured upon cooling. For the data presented in
Fig. 4, the ARPES spectra were divided by the energy resolution-convoluted
Fermi-Dirac distribution function for each temperature, and a momentum-
independent background was obtained by integrating the featureless momentum
range -0.66 $\mathrm{\AA^{-1}}$$<k_{y}<$-0.6 $\mathrm{\AA^{-1}}$ was
subtracted. The EDCs in Fig. 4e were obtained in different experiment from the
data in Fig. 4f. ARPES data in Fig. 4e were acquired with longer measurement
time and larger temperature steps for a high S/N ratio.

in 1,…,0 See pages \x, of SHG_Supplementray_230513.pdf