Supplementary Information for:

Exact numerical solution of the
classical and quantum Heisenberg spin glass

Nikita Kavokine, Markus Müller, Antoine Georges and Olivier Parcollet
\-

###### Contents

  1. 1 Solution of the mean-field Heisenberg spin glass
    1. 1.1 Model definition
    2. 1.2 Replica trick
    3. 1.3 Replica symmetry breaking
    4. 1.4 Internal energy
    5. 1.5 Doped quantum spin glass
  2. 2 Numerical procedures
    1. 2.1 Solution of the Parisi equations
    2. 2.2 Solution of the quantum impurity problem
    3. 2.3 General solution procedure
  3. 3 Supplementary figures

##  1 Solution of the mean-field Heisenberg spin glass

In this section, we describe the analytical procedure that reduces the mean-
field (infinite dimensional) Heisenberg spin glass problem to a self-
consistent single spin problem, amenable to an exact numerical solution. We
make use of the replica method, as introduced by Sherrington and Kirkaptrick
in the classical case [1] and by Bray and Moore in the quantum case [2]. We
then follow Parisi’s seminal replica symmetry breaking construction [3, 4],
which was generalized to the quantum case of the transverse field Sherrington-
Kirkpatrick model [5, 6], but not to the quantum Heisenberg spin glass.

###  1.1 Model definition

We consider an all-to-all interacting system of $N$ spins $\mathbf{S}_{i}$,
described by the Hamiltonian

$H=-\sum_{i<j}J_{ij}\mathbf{S}_{i}\cdot\mathbf{S}_{j}.$ (1)

The $J_{ij}$ are independent random variables with distribution
$\mathbb{G}(J_{ij}|J^{2}/N)$, where $\mathbb{G}(x|v)$ is the normalized
Gaussian with zero mean and variance $v$. We will consider both a classical
and a quantum version of the model. In the classical case, the
$\mathbf{S}_{i}$ are vectors constrained to reside on the $\ell$-dimensional
sphere of radius $S$. In the quantum case, the $\mathbf{S}_{i}$ are operators,
forming a representation of spin $S$ of the group SU(2); we denote $\ell=3$
the number of spin components.

###  1.2 Replica trick

Our aim is to compute $f$, the thermodynamic limit of the free energy per site
at inverse temperature $\beta$. To this end, we introduce the quantity
$f_{n}(N)\equiv(-1/\beta N)\log(\overline{Z^{n}})$, where the bar denotes
averaging over the disorder. $Z^{n}$ is the partition function of $n$ replicas
of the $N$-spin system, with the same realization of the disorder. Allowing
$n$ to take non-integer values through analytical continuation, we find

$\frac{f_{n}(N)}{n}=-\frac{1}{\beta nN}\log\left(\overline{e^{n\log
Z}}\right)\underset{n\to 0}{=}-\frac{1}{\beta N}\overline{\log Z}.$ (2)

Since the free energy is self-averaging in the thermodynamic limit,

$f=\lim_{N\to\infty}\lim_{n\to 0}\frac{f_{n}(N)}{n}.$ (3)

We will assume that the order of the limits can be reversed, and first carry
out the computation of $f_{n}$ in the thermodynamic limit.

####  1.2.1 Classical case

We average over the random couplings the partition function of $n$ replicas of
the spin system, labeled by the index $a$:

$\displaystyle\overline{Z^{n}}$
$\displaystyle=\sum_{\\{\mathbf{S}_{i}^{a}\\}}\int_{-\infty}^{+\infty}\prod_{i<j}\mathrm{d}J_{ij}\,\mathbb{G}(J_{ij}|J^{2}/N)\exp\left[\beta\sum_{a}\sum_{i<j}J_{ij}\mathbf{S}_{i}^{a}\cdot\mathbf{S}_{j}^{a}\right]$
(4)
$\displaystyle=\sum_{\\{\mathbf{S}_{i}^{a}\\}}\exp\left[\frac{\beta^{2}J^{2}}{2N}\sum_{i<j}\sum_{a,b}(\mathbf{S}_{i}^{a}\cdot\mathbf{S}_{j}^{a})(\mathbf{S}_{i}^{b}\cdot\mathbf{S}_{j}^{b})\right]$
(5)
$\displaystyle=\sum_{\\{\mathbf{S}_{i}^{a}\\}}\exp\left[\frac{\beta^{2}J^{2}}{4N}\sum_{\alpha,\beta}\sum_{a,b}\left(\sum_{i}S_{i,\alpha}^{a}\cdot
S_{i,\beta}^{b}\right)^{2}-n^{2}\frac{\beta^{2}J^{2}NS^{4}}{4}\right].$ (6)

Here, the indices $\alpha,\beta$ denote the spin components. We may drop in
the following the term proportional to $n^{2}$, since it vanishes in the $n\to
0$ limit. We now perform a Hubbard-Stratonovitch transformation to decouple
the sum over spins:

$\overline{Z^{n}}=\int\prod_{a,b,\alpha,\beta}\frac{\mathrm{d}Q^{\alpha\beta}_{ab}}{\sqrt{4\pi/(N\beta^{2}J^{2})}}e^{-N\mathcal{F}[\\{Q^{\alpha\beta}_{ab}\\}]},$
(7)

with

$\mathcal{F}[\\{Q^{\alpha\beta}_{ab}\\}]=\frac{\beta^{2}J^{2}}{4}\sum_{a,b}\sum_{\alpha,\beta}(Q_{ab}^{\alpha\beta})^{2}-\log\left[\sum_{\\{\mathbf{S}^{a}\\}}\exp\left(\frac{\beta^{2}J^{2}}{2}\sum_{a,b}\sum_{\alpha,\beta}Q^{\alpha\beta}_{ab}\,S^{a}_{\alpha}\cdot
S^{b}_{\beta}\right)\right].$ (8)

In the thermodynamic limit, the integral in Eq. (7) is dominated by the saddle
point $\\{Q^{\alpha\beta\star}_{ab}\\}$:

$f_{n}(N)=(-1/\beta
N)\log(\overline{Z^{n}})\underset{N\to\infty}{\longrightarrow}\frac{1}{\beta}\mathcal{F}[\\{Q_{ab}^{\alpha\beta\star}\\}],$
(9)

with $\partial\mathcal{F}/\partial
Q_{ab}^{\alpha\beta}|_{Q_{ab}^{\alpha\beta\star}}=0,\forall(a,b,\alpha,\beta)$.
The latter amounts to

$Q_{ab}^{\alpha\beta\star}=\langle S^{a}_{\alpha}\cdot
S^{b}_{\beta}\rangle_{H_{\rm loc}[Q^{\star}},~{}~{}~{}~{}~{}H_{\rm
loc}[Q^{\star}]=-\frac{\beta
J^{2}}{2}\sum_{a,b}\sum_{\alpha,\beta}Q_{ab}^{\alpha\beta\star}S^{a}_{\alpha}\cdot
S^{b}_{\beta}.$ (10)

####  1.2.2 Quantum case

In the quantum case, we use an imaginary time path integral formalism to
express the partition function:

$\overline{Z^{n}}=\sum_{\\{\mathbf{S}_{i}^{a}\\}}\int_{-\infty}^{+\infty}\prod_{i<j}\mathrm{d}J_{ij}\,\mathbb{G}(J_{ij}|J^{2}/N)\exp\left[\beta\sum_{a}\sum_{i<j}J_{ij}\int_{0}^{1}\mathrm{d}\tau\,\mathbf{S}_{i}^{a}(\tau)\cdot\mathbf{S}_{j}^{a}(\tau)\right],$
(11)

where the notation $\sum_{\mathbf{S}}[.]$ is now taken to mean the path
integral $\int[D\mathbf{S}(\tau)]e^{-\mathcal{S}_{0}[\mathbf{S}(\tau)]}[.]$,
where $\mathcal{S}_{0}[\mathbf{S}(\tau)]$ is the imaginary time action
describing an isolated spin, which we do not need to specify explicitly. We
may then proceed as in the classical case:

$\displaystyle\overline{Z^{n}}$
$\displaystyle=\sum_{\\{\mathbf{S}_{i}^{a}\\}}\exp\left[\frac{\beta^{2}J^{2}}{2N}\int_{0}^{1}\mathrm{d}\tau\mathrm{d}\tau^{\prime}\sum_{i<j}\sum_{a,b}(\mathbf{S}_{i}^{a}(\tau)\cdot\mathbf{S}_{j}^{a}(\tau^{\prime}))(\mathbf{S}_{i}^{b}(\tau)\cdot\mathbf{S}_{j}^{b}(\tau^{\prime}))\right]$
(12)
$\displaystyle=\sum_{\\{\mathbf{S}_{i}^{a}\\}}\exp\left[\frac{\beta^{2}J^{2}}{4N}\int_{0}^{1}\mathrm{d}\tau\mathrm{d}\tau^{\prime}\sum_{\alpha,\beta}\sum_{a,b}\left(\sum_{i}S_{i,\alpha}^{a}(\tau)\cdot
S_{i,\beta}^{b}(\tau^{\prime})\right)^{2}-n^{2}\frac{\beta^{2}J^{2}N}{4}S^{2}(S+1)^{2}\right]$
(13)

We may again drop the term proportional to $n^{2}$ as it vanishes in the $n\to
0$ limit. The Hubbard-Stratonovitch transformation yields

$\overline{Z^{n}}=\int\prod_{a,b,\alpha,\beta}\int_{0}^{1}\mathrm{d}\tau\mathrm{d}\tau^{\prime}\,\frac{\mathrm{d}Q_{ab}^{\alpha,\beta}(\tau,\tau^{\prime})}{\sqrt{4\pi/(N\beta^{2}J^{2})}}e^{-N\mathcal{F}[\\{Q^{\alpha\beta}_{ab}(\tau,\tau^{\prime})\\}]},$
(14)

with

$\begin{split}\mathcal{F}[\\{Q_{ab}\\}]=&\frac{\beta^{2}J^{2}}{4}\int_{0}^{1}\mathrm{d}\tau\mathrm{d}\tau^{\prime}\sum_{a,b}\sum_{\alpha,\beta}(Q_{ab}^{\alpha\beta}(\tau,\tau^{\prime}))^{2}\\\
&-\log\left[\sum_{\\{\mathbf{S}^{a}\\}}\exp\left(\frac{\beta^{2}J^{2}}{2}\int_{0}^{1}\mathrm{d}\tau\mathrm{d}\tau^{\prime}\sum_{a,b}\sum_{\alpha,\beta}Q^{\alpha\beta}_{ab}(\tau,\tau^{\prime})\,S^{a}_{\alpha}(\tau)\cdot
S^{b}_{\beta}(\tau^{\prime})\right)\right].\end{split}$ (15)

The saddle point approximation in the thermodynamic limit yields

$f_{n}(N)=(-1/\beta
N)\log(\overline{Z^{n}})\underset{N\to\infty}{\longrightarrow}\frac{1}{\beta}\mathcal{F}[\\{Q_{ab}^{\alpha\beta\star}(\tau,\tau^{\prime})\\}],$
(16)

with $\partial\mathcal{F}/\partial
Q_{ab}^{\alpha\beta}(\tau,\tau^{\prime})|_{Q_{ab}^{\alpha\beta\star}(\tau,\tau^{\prime})}=0,\forall(a,b,\alpha,\beta,\tau,\tau^{\prime})$.
The latter amounts to

$Q^{\alpha\beta\star}_{ab}(\tau,\tau^{\prime})=\langle\mathrm{T}\,S^{a}_{\alpha}(\tau)\cdot
S^{b}_{\beta}(\tau^{\prime})\rangle_{\mathcal{S}_{\rm loc}[Q^{\star}]},$ (17)
$\mathcal{S}_{\rm
loc}[Q^{\star}]=\sum_{a}\mathcal{S}_{0}[\mathbf{S}^{a}(\tau)]-\frac{\beta^{2}J^{2}}{2}\int_{0}^{1}\mathrm{d}\tau\mathrm{d}\tau^{\prime}\sum_{a,b}\sum_{\alpha,\beta}Q_{ab}^{\alpha\beta\star}(\tau,\tau^{\prime})S^{a}_{\alpha}(\tau)\cdot
S^{b}_{\beta}(\tau^{\prime}).$ (18)

####  1.2.3 Interpretation of the $Q_{ab}$

We have reduced the initial lattice problem to a self-consistent local problem
of $n$ coupled replicas in the limit $n\to 0$, which remains to be defined. We
now establish the link between the spin-spin correlation functions of the
local problem and those of the initial lattice problem. This will allow us to
simplify some of the coefficients $Q_{ab}$ that define the local problem. For
generality, we present here the reasoning in the quantum case.

Let $\mathcal{S}$ be the action corresponding to the initial lattice problem.
We supplement it with a source term:

$\mathcal{S}\mapsto\mathcal{S}+\beta^{2}h\int_{0}^{1}\mathrm{d}\tau\mathrm{d}\tau^{\prime}\sum_{i}S_{i\alpha}(\tau)S_{i\beta}(\tau^{\prime}).$
(19)

Then, in the thermodynamic limit,

$\begin{split}\int_{0}^{1}\mathrm{d}\tau\mathrm{d}\tau^{\prime}\,\overline{\langle\mathrm{T}S_{i\alpha}(\tau)S_{i\beta}(\tau^{\prime})\rangle|_{\mathcal{S}}}&=\partial
f/\partial h|_{h=0}\\\ &=\lim_{n\to 0}\partial(f_{n}(N\to\infty)/n)/\partial
h|_{h=0}\\\ &=-\lim_{n\to 0}\frac{1}{n}\left.\frac{\partial}{\partial
h}\right|_{h=0}\log\sum_{\\{\mathbf{S}^{a}\\}}e^{-\mathcal{S}_{\rm
loc}[Q^{\star}]-h\int_{0}^{1}\mathrm{d}\tau\mathrm{d}\tau^{\prime}\sum_{a}S^{a}_{\alpha}(\tau)S^{a}_{\beta}(\tau^{\prime})}\\\
&=\lim_{n\to
0}\frac{1}{n}\sum_{a}\int_{0}^{1}\mathrm{d}\tau\mathrm{d}\tau^{\prime}\langle\mathrm{T}S^{a}_{\alpha}(\tau)S^{a}_{\beta}(\tau^{\prime})\rangle|_{\mathcal{S}_{\rm
loc}[Q^{\star}]}\\\ &=\lim_{n\to
0}\int_{0}^{1}\mathrm{d}\tau\mathrm{d}\tau^{\prime}\langle\mathrm{T}S^{a}_{\alpha}(\tau)S^{a}_{\beta}(\tau^{\prime})\rangle|_{\mathcal{S}_{\rm
loc}[Q^{\star}]}.\end{split}$ (20)

Since the correlation function on the left hand side vanishes for two
different spin components, the $Q_{aa}^{\alpha\beta\star}$ vanish if
$\alpha\neq\beta$. We then set $Q(\tau-\tau^{\prime})\equiv
Q_{aa}^{\alpha\beta\star}(\tau,\tau^{\prime})$. Consider now two copies of the
initial system, governed by the action

$\mathcal{S}_{\rm
II}=\mathcal{S}[\\{\mathbf{S}_{i}^{1}\\}]+\mathcal{S}[\\{\mathbf{S}_{i}^{2}\\}]$
(21)

By applying the replica method to this action in the same way as above, we
find the correspondence

$\overline{\langle\mathrm{T}S^{1}_{i\alpha}(\tau)S^{2}_{i\beta}(\tau^{\prime})\rangle|_{\mathcal{S_{\rm
II}}}}=\lim_{n\to 0,a\neq
b}\langle\mathrm{T}S^{a}_{\alpha}(\tau)S^{b}_{\beta}(\tau^{\prime})\rangle|_{\mathcal{S}_{\rm
loc}[Q^{\star}]}$ (22)

As the copies 1 and 2 are independent, the right hand side of the above
equation, and therefore the $Q_{ab}^{\alpha\beta\star},a\neq b$, cannot depend
on $\tau,\tau^{\prime}$. Similarly, in the absence of quadrupolar order in the
initial system, the $Q_{ab}^{\alpha\beta\star}$ vanish when $a\neq b$ and
$\alpha\neq\beta$. The reasoning is analogous in the classical case. Thanks to
the normalization of the classical spins, Eq. (8) may be further simplified by
integrating over the $Q_{aa}^{\alpha\alpha}$ in Eq. (7).

Let us summarize our results so far.

In the classical case, the free energy per spin is

$f=-\frac{\beta J^{2}S^{4}}{4\ell}+\lim_{n\to 0}\frac{1}{\beta
n}\left[\frac{\ell\beta^{2}J^{2}}{2}\sum_{a<b}Q_{ab}^{2}-\log Z_{\rm loc}^{\rm
Cl}\right],$ (23)

with

$Z_{\rm loc}^{\rm Cl}=\sum_{\\{\mathbf{S}^{a}\\}}e^{-\beta H_{\rm
loc}},~{}~{}~{}~{}~{}Q_{ab}=\frac{1}{\ell}\langle\mathbf{S}^{a}\cdot\mathbf{S}^{b}\rangle_{H_{\rm
loc}},~{}~{}~{}~{}~{}H_{\rm loc}=-\beta
J^{2}\sum_{a<b}Q_{ab}\mathbf{S}^{a}\cdot\mathbf{S}^{b}.$ (24)

In the quantum case,

$f=\frac{\ell\beta J^{2}}{4}\int_{0}^{1}\mathrm{d}\tau Q(\tau)^{2}+\lim_{n\to
0}\frac{1}{\beta
n}\left\\{\frac{\ell\beta^{2}J^{2}}{2}\sum_{a<b}Q_{ab}^{2}-\log Z_{\rm
loc}^{\rm Q}\right\\},$ (25)

with

$Z_{\rm loc}^{\rm Q}=\int\prod_{a}[D\mathbf{S}^{a}(\tau)]e^{-\mathcal{S}_{\rm
loc}},~{}Q_{ab}=\frac{1}{\ell}\langle\mathrm{T}\,\mathbf{S}^{a}(\tau)\cdot\mathbf{S}^{b}(\tau^{\prime})\rangle_{\mathcal{S}_{\rm
loc}},~{}Q(\tau)=\frac{1}{\ell}\langle\mathrm{T}\,\mathbf{S}^{a}(\tau)\cdot\mathbf{S}^{a}(0)\rangle_{\mathcal{S}_{\rm
loc}},$ (26)

and

$\begin{split}\mathcal{S}_{\rm
loc}=\sum_{a}\mathcal{S}_{0}[\mathbf{S}^{a}(\tau)]-\frac{\beta^{2}J^{2}}{2}\int_{0}^{1}\mathrm{d}\tau\mathrm{d}\tau^{\prime}Q(\tau-\tau^{\prime})\sum_{a}\mathbf{S}^{a}(\tau)\cdot\mathbf{S}^{a}(\tau^{\prime})\\\
-\beta^{2}J^{2}\sum_{a<b}Q_{ab}\int_{0}^{1}\mathrm{d}\tau\mathbf{S}^{a}(\tau)\int_{0}^{1}\mathrm{d}\tau^{\prime}\mathbf{S}^{b}(\tau^{\prime}).\end{split}$
(27)

We now need to specify an ansatz for the matrix $Q_{ab}$ which allows for
analytic continuation to non-integer matrix dimensions.

###  1.3 Replica symmetry breaking

####  1.3.1 Recursive construction

We use Parisi’s replica symmetry breaking (RSB) ansatz for the matrix
$Q_{ab}$. The matrix $Q^{\mathrm{RSB}_{k}}$ at the $k^{\rm th}$ stage of RSB
is defined by the choice of two sequences: $(q_{p})_{p=0..k}$ and
$(m_{p})_{p=0..k}$, the latter being strictly decreasing, with $m_{0}=n$. It
is built "from the inside out" according to the following prescription:

  1. 1.

$Q^{\mathrm{RSB}_{k}}_{k}=(q_{k}-q_{k-1})U_{m_{k}}$;

  2. 2.

$\forall
p\in[0,k-1],Q^{\mathrm{RSB}_{k}}_{p}=\mathrm{diag}_{m_{p}/m_{p+1}}(Q^{\mathrm{RSB}_{k}}_{p+1})+(q_{p}-q_{p-1})U_{m_{p}}$;

  3. 3.

$Q^{\mathrm{RSB}_{k}}_{0}\equiv Q^{\mathrm{RSB}_{k}}$.

Here $U_{m}$ is the matrix of size $m$ filled with ones of size,
$\mathrm{diag}_{r}(A)$ is the matrix with $r$ times the block $A$ on the
diagonal, and we set $q_{-1}=0$. We note that this construction places $q_{k}$
instead of 0 on the diagonal, and the corresponding term needs to be
subtracted in the end. In the classical case, the local Hamiltonian at the
$k^{\rm th}$ stage of RSB is therefore

$H_{\rm loc}^{\mathrm{RSB}_{k}}=nq_{k}\frac{\beta J^{2}S^{2}}{2}-\frac{\beta
J^{2}}{2}\sum_{ab}(Q^{\mathrm{RSB}_{k}})^{ab}\,\mathbf{S}^{a}\cdot\mathbf{S}^{b},$
(28)

which we may write as

$H_{\rm loc}^{\mathrm{RSB}_{k}}=\sum_{a}H_{k}[\mathbf{S}_{a}]-\frac{\beta
J^{2}}{2}\sum_{ab}(Q^{\mathrm{RSB}_{k}})^{ab}\,\mathbf{S}^{a}\cdot\mathbf{S}^{b}.$
(29)

In the quantum case, we introduce the notation
$\overline{\mathbf{S}}\equiv\int_{0}^{1}\mathrm{d}\tau\,\mathbf{S}(\tau)$.
Then, the local action at the $k^{\rm th}$ stage of RSB is

$\mathcal{S}_{\rm
loc}^{\mathrm{RSB}_{k}}=\sum_{a}\mathcal{S}_{0}[\mathbf{S}^{a}(\tau)]-\int_{0}^{1}\mathrm{d}\tau\mathrm{d}\tau^{\prime}\frac{\beta^{2}J^{2}(Q(\tau^{\prime}-\tau)-q_{k})}{2}\sum_{a}\mathbf{S}^{a}(\tau)\cdot\mathbf{S}^{a}(\tau^{\prime})-\frac{\beta^{2}J^{2}}{2}\sum_{ab}(Q^{\mathrm{RSB}_{k}})^{ab}\,\overline{\mathbf{S}}^{a}\cdot\overline{\mathbf{S}}^{b}.$
(30)

For simplicity, we denote

$\mathcal{S}_{\rm
loc}^{\mathrm{RSB}_{k}}=\sum_{a}\mathcal{S}_{k}[\mathbf{S}^{a}(\tau)]-\frac{\beta^{2}J^{2}}{2}\sum_{ab}(Q^{\mathrm{RSB}_{k}})^{ab}\,\overline{\mathbf{S}}^{a}\cdot\overline{\mathbf{S}}^{b}.$
(31)

We now establish a recursion between the local partition functions $Z_{\rm
loc}$ at steps $p$ and $p+1$ of the recursive construction. We use the quantum
notation, but the relation will be identical in the classical case. At step
$p$ the local partition function under an external field $\mathbf{h}$ reads:

$Z_{p}^{\mathrm{RSB}_{k}}(\mathbf{h})=\int\prod_{a}[D\mathbf{S}^{a}]\exp\left[-\sum_{a}\mathcal{S}_{k}[\\{\mathbf{S}^{a}\\}]+\frac{\beta^{2}J^{2}}{2}\sum_{ab}(Q^{\mathrm{RSB}_{k}}_{p})^{ab}\,\overline{\mathbf{S}}^{a}\cdot\overline{\mathbf{S}}^{b}+\beta\mathbf{h}\sum_{a}\overline{\mathbf{S}}^{a}\right].$
(32)

We note that $\mathbf{h}$ is introduced only as an auxiliary term, and we will
set $\mathbf{h}=0$ in the end. In particular it does not represent an actual
external field applied on the initial system: the matrix $Q_{ab}$ would then
no longer be rotationally invariant. At the first step ($p=k$), we find

$\begin{split}Z_{k}^{\mathrm{RSB}_{k}}(\mathbf{h})&=\int\prod_{a}[D\mathbf{S}^{a}]\exp\left[-\sum_{a=1}^{m_{k}}\mathcal{S}_{k}[\\{\mathbf{S}^{a}\\}]+\frac{\beta^{2}J^{2}}{2}(q_{k}-q_{k-1})\left(\sum_{a=1}^{m_{k}}\overline{\mathbf{S}}^{a}\right)^{2}+\beta\mathbf{h}\sum_{a=1}^{m_{k}}\overline{\mathbf{S}}^{a}\right]\\\
&=\int\mathrm{d}\mathbf{h}_{k}\,\mathbb{G}(\mathbf{h}_{k}|2J^{2}(q_{k}-q_{k-1}))\left(\int[D\mathbf{S}]e^{-\mathcal{S}_{k}[\mathbf{S}]+\beta(\mathbf{h}+\mathbf{h}_{k})\overline{\mathbf{S}}}\right)^{m_{k}}\\\
&\equiv\int\mathrm{d}\mathbf{h}_{k}\,\mathbb{G}(\mathbf{h}_{k}|2J^{2}(q_{k}-q_{k-1}))z_{k}(\mathbf{h}+\mathbf{h}_{k})^{m_{k}}.\end{split}$
(33)

At step $p$, we may consider separately the effect of the two terms in
$Q^{\mathrm{RSB}_{k}}_{p}$:

$\begin{split}Z_{p}^{\mathrm{RSB}_{k}}(\mathbf{h})=\int\prod_{a}[D\mathbf{S}^{a}]\exp\left[-\sum_{a=1}^{m_{p}}\mathcal{S}_{k}[\\{\mathbf{S}^{a}\\}]+\frac{\beta^{2}J^{2}}{2}(q_{p}-q_{p-1})\left(\sum_{a=1}^{m_{p}}\overline{\mathbf{S}}^{a}\right)^{2}\right.\\\
\left.+\frac{\beta^{2}J^{2}}{2}\sum_{\alpha=1}^{m_{p}/m_{p+1}}\sum_{a,b\in\alpha}(Q^{\mathrm{RSB}_{k}}_{p+1})^{ab}\overline{\mathbf{S}}^{a}\cdot\overline{\mathbf{S}}^{b}+\beta\mathbf{h}\sum_{a=1}^{m_{k}}\overline{\mathbf{S}}^{a}\right].\end{split}$
(34)

Here, $\alpha$ labels the blocks of the matrix $Q^{\mathrm{RSB}_{k}}_{p}$.
After a Hubbard-Stratonovitch transformation on the squared sum, we find the
recurrence relation

$Z_{p}^{\mathrm{RSB}_{k}}(\mathbf{h})=\int\mathrm{d}\mathbf{h}_{p}\,\mathbb{G}(\mathbf{h}_{p}|2J^{2}(q_{p}-q_{p-1}))\left[Z^{\mathrm{RSB}_{k}}_{p+1}(\mathbf{h}+\mathbf{h}_{p})\right]^{m_{p}/m_{p+1}}.$
(35)

####  1.3.2 Limit $n\to 0,k\to\infty$

We now take the simultaneous limit $n\to 0$ and $k\to\infty$, keeping
$m_{k}=1$. Then, $m_{p}$ becomes a continuous variable $x$ between 0 and 1,
and $q_{p}\to q(x)$, such that $m_{p}\to x$. $q(x)$ is the Parisi order
parameter. The corresponding limit of the recurrence relation (35) may be
taken if one uses the formal identity

$\left[\exp\left(\frac{1}{2}v\bm{\nabla}^{2}\right)\cdot
f\right](\mathbf{h})=\int\mathrm{d}\mathbf{h}^{\prime}\mathbb{G}(\mathbf{h}^{\prime}|v)f(\mathbf{h}+\mathbf{h}^{\prime}),$
(36)

that holds for any sufficiently regular function $f$.
$Z_{p}^{\mathrm{RSB}_{k}}(\mathbf{h})$ then becomes a continuous function
$\zeta(x,\mathbf{h})$, satisfying

$\zeta(x-\mathrm{d}x,\mathbf{h})=\exp\left(\frac{1}{2}J^{2}\mathrm{d}q(x)\bm{\nabla}^{2}\right)\zeta(x,\mathbf{h})^{1-\mathrm{d}x/x}$
(37)

This may be recast into a partial differential equation:

$\frac{\partial\zeta}{\partial
x}=-\frac{1}{2}J^{2}\frac{\mathrm{d}q}{\mathrm{d}x}\bm{\nabla}^{2}\zeta+\frac{1}{x}\zeta\log\zeta,$
(38)

with boundary condition $\zeta(1,h)=z_{\infty}(\mathbf{h})$ as defined below.
At this point, a distinction needs to be made between the classical and the
quantum case. In the classical case, the "single-replica" part of the local
Hamiltonian $H_{k}$ does not depend on the spin. We may hence split it off in
the partition function, so that applying Eq. (35) at $p=0$ yields

$\lim_{n\to 0}\frac{1}{n}\log Z_{\rm loc}^{\rm
Cl}=-\frac{\beta^{2}J^{2}S^{2}}{2}q(1)+\int\mathrm{d}\mathbf{h}\,\mathbb{G}(\mathbf{h}|2J^{2}q(0))\lim_{x\to
0}\left[\frac{1}{x}\log\zeta(x,\mathbf{h})\right],$ (39)

and the boundary condition is explicitly

$z_{\infty}^{\rm Cl}(\mathbf{h})=\int\mathrm{d}\mathbf{S}e^{-\beta
H_{\infty}(\mathbf{h})}=2\pi\int_{0}^{\pi}\mathrm{d}\theta\sin\theta e^{-\beta
hS\cos\theta}=\frac{2\pi}{\beta hS}(e^{\beta hS}-e^{-\beta hS}).$ (40)

We may define $\phi(x,\mathbf{h})\equiv(1/x)\log\zeta(x,\mathbf{h})$, which
satisfies

$\frac{\partial\phi}{\partial
x}=-\frac{1}{2}J^{2}\frac{\mathrm{d}q}{\mathrm{d}x}\left(\bm{\nabla}^{2}\phi+x(\bm{\nabla}\phi)^{2}\right).$
(41)

We also anticipate that $q(0)=0$, so that the Gaussian in Eq. (39) becomes a
$\delta$ function.

Then, introducing our results into Eq. (23) the free energy of the classical
Heisenberg spin glass is given by

$f=-\frac{\beta J^{2}S^{4}}{4\ell}-\frac{\ell\beta
J^{2}}{4}\int_{0}^{1}\mathrm{d}x\,q(x)^{2}+\frac{\beta
J^{2}S^{2}}{2}q(1)-\frac{1}{\beta}\phi(0,0),$ (42)

where $\phi(x,h)$ satisfies Eq. (41) with $\phi(1,\mathbf{h})=\log z^{\rm
Cl}_{\infty}(\mathbf{h})$.

We have used that, with the RSB ansatz [7],

$\lim_{n\to
0}\sum_{a<b}Q_{ab}^{2}=-\frac{1}{2}\int_{0}^{1}\mathrm{d}x\,q(x)^{2}.$ (43)

This provides the solution of the classical spin glass problem given a
function $q(x)$.

In the quantum case, the single-replica part of the action cannot be
simplified and we obtain

$\lim_{n\to 0}\frac{1}{n}\log Z_{\rm loc}^{\rm
Q}=\int\mathrm{d}\mathbf{h}\,\mathbb{G}(\mathbf{h}|2J^{2}q(0))\lim_{x\to
0}\left[\frac{1}{x}\log\zeta(x,\mathbf{h})\right],$ (44)

and

$z_{\infty}^{\rm
Q}(\mathbf{h})=\int[D\mathbf{S}]e^{-\mathcal{S}_{\infty}[\mathbf{S}(\tau),\mathbf{h}]},$
(45)

with

$\mathcal{S}_{\infty}(\mathbf{h})=\mathcal{S}_{0}[\mathbf{S}(\tau)]-\int_{0}^{1}\mathrm{d}\tau\mathrm{d}\tau^{\prime}\frac{\beta^{2}J^{2}(Q(\tau-\tau^{\prime})-q(1))}{2}\mathbf{S}(\tau)\cdot\mathbf{S}(\tau^{\prime})+\beta\mathbf{h}\cdot\mathbf{S}$
(46)

Then, introducing our results into Eq. (25) the free energy of the quantum
Heisenberg spin glass is given by

$f=\frac{\ell\beta J^{2}}{4}\left[\int_{0}^{1}\mathrm{d}\tau
Q(\tau)^{2}-\int_{0}^{1}\mathrm{d}x\,q(x)^{2}\right]-\frac{1}{\beta}\phi(0,0),$
(47)

where $\phi(x,h)$ satisfies Eq. (41) with $\phi(1,\mathbf{h})=\log z^{\rm
Q}_{\infty}(\mathbf{h})$.

This provides the solution of the quantum spin glass problem given $q(x)$ and
$Q(\tau)$.

####  1.3.3 Local observables and two point functions

We now determine the functions $q(x)$ and $Q(\tau)$ by applying the RSB ansatz
to the self-consistency conditions in Eqs. (24) and (26). We go back to the
recursive construction in 1.3.1 and consider the observable
$R^{a}(\tau,\tau^{\prime})=\mathbf{S}^{a}(\tau)\mathbf{S}^{a}(\tau^{\prime})$,
which is local in replica space. For every $\tau,\tau^{\prime}$, we define the
quantity $R^{\mathrm{RSB}_{k}}_{p}(\mathbf{h})$, which corresponds to the
average value of $R$ computed with the matrix $Q^{\mathrm{RSB}_{k}}$ built up
to step $p$. $R^{\mathrm{RSB}_{k}}_{p}$ is only slightly different from
$Z^{\mathrm{RSB}_{k}}_{p}$: indeed, at a given step, only one of the blocks is
perturbed by the insertion of $R$, while the others are not. Therefore, one
has the recurrence relation

$R_{p}^{\mathrm{RSB}_{k}}(\mathbf{h})=\int\mathrm{d}\mathbf{h}_{p}\,\mathbb{G}(\mathbf{h}_{p}|2J^{2}(q_{p}-q_{p-1}))\frac{R^{\mathrm{RSB}_{k}}_{p+1}(\mathbf{h}+\mathbf{h}_{p})Z^{\mathrm{RSB}_{k}}_{p+1}(\mathbf{h}+\mathbf{h}_{p})^{m_{p}/m_{p+1}}}{Z^{\mathrm{RSB}_{k}}_{p}(\mathbf{h})}.$
(48)

The situation is slightly different for an observable that is non-local in
replica space, in particular the spin-spin correlator
$\langle\overline{\mathbf{S}}^{a}\overline{\mathbf{S}}^{b}\rangle^{\mathrm{RSB}_{k}}$.
When $k$ is sufficiently large, $a$ and $b$ are always in different blocks at
$k$-RSB. If at step $p$ they are still in different blocks, then
$\langle\overline{\mathbf{S}}^{a}\overline{\mathbf{S}}^{b}\rangle^{\mathrm{RSB}_{k}}_{p}=(\langle\overline{\mathbf{S}}\rangle^{\mathrm{RSB}_{k}}_{p})^{2}$.
Let $p_{0}$ by the first step where $\overline{\mathbf{S}}^{a}$ and
$\overline{\mathbf{S}}^{b}$ are in the same block. Then for $p<p_{0}$, the
spin-spin correlator satisfies the same recurrence relation as a local
observable:

$\langle\overline{\mathbf{S}}^{a}\overline{\mathbf{S}}^{b}\rangle^{\mathrm{RSB}_{k}}_{p}(\mathbf{h})=\int\mathrm{d}\mathbf{h}_{p}\,\mathbb{G}(\mathbf{h}_{p}|2J^{2}(q_{p}-q_{p-1}))\frac{\langle\overline{\mathbf{S}}^{a}\overline{\mathbf{S}}^{b}\rangle^{\mathrm{RSB}_{k}}_{p+1}(\mathbf{h}+\mathbf{h}_{p})Z^{\mathrm{RSB}_{k}}_{p+1}(\mathbf{h}+\mathbf{h}_{p})^{m_{p}/m_{p+1}}}{Z^{\mathrm{RSB}_{k}}_{p}(\mathbf{h})}.$
(49)

In the limit $n\to 0,k\to\infty$, $R_{p}^{\mathrm{RSB}_{k}}(\mathbf{h})$
becomes a function $\rho(x,\mathbf{h})$, satisfying

$\zeta(x-\mathrm{d}x,\mathbf{h})\rho(x-\mathrm{d}x,\mathbf{h})=\exp\left(\frac{1}{2}J^{2}\mathrm{d}q(x)\bm{\nabla}^{2}\right)\rho(x,\mathbf{h})\zeta(x,\mathbf{h})^{1-\mathrm{d}x/x}.$
(50)

The corresponding PDE is

$\frac{\partial\rho}{\partial
x}=-\frac{1}{2}J^{2}\frac{\mathrm{d}q}{\mathrm{d}x}\left(\bm{\nabla}^{2}\rho+2(\bm{\nabla}\log\zeta)\cdot(\bm{\nabla}\rho)\right),$
(51)

with the boundary condition
$\rho(1,\mathbf{h})_{\tau,\tau^{\prime}}=\langle\mathbf{S}(\tau)\mathbf{S}(\tau^{\prime})\rangle_{\infty,\mathbf{h}}$,
where $\langle\cdot\rangle_{\infty,\mathbf{h}}$ is the mean value with respect
to $H_{\infty}(\mathbf{h})$ or $\mathcal{S}_{\infty}(\mathbf{h})$. If we now
define $\beta\mathbf{s}(x,\mathbf{h})=(1/x)\bm{\nabla}\log\zeta$, we find that
it satisfies the same equation as $\rho$, with the boundary condition
$\mathbf{s}(1,\mathbf{h})=\langle\mathbf{S}\rangle_{\infty}$. Then $\rho$
satisfies

$\frac{\partial\rho}{\partial
x}=-\frac{1}{2}J^{2}\frac{\mathrm{d}q}{\mathrm{d}x}\left(\bm{\nabla}^{2}\rho+2\beta
x(\mathbf{s}\cdot\bm{\nabla})\rho\right)$ (52)

Using Eq. (48) with $p=0$,

$Q(\tau-\tau^{\prime})=\int\mathrm{d}\mathbf{h}\mathbb{G}(\mathbf{h}|2J^{2}q(0))\rho(0,\mathbf{h})_{\tau,\tau^{\prime}}=\rho(0,0)_{\tau,\tau^{\prime}}.$
(53)

The spin-spin correlation function
$\langle\mathbf{S}^{a}\mathbf{S}^{b}\rangle^{\mathrm{RSB}_{k}}_{p}$ becomes a
continuous function $\nu(y,x,\mathbf{h})$, where $y$ is the continuum analogue
of the index $p_{0}$ defined above. For $y<x$, $\nu$ satisfies Eq. (52), with
the boundary condition $\nu(y=x,x,\mathbf{h})=\mathbf{s}^{2}(x,\mathbf{h})$.
The self-consistency conditions (24) and (26) become

$q(x)=\frac{1}{\ell}\int\mathrm{d}\mathbf{h}\,\mathbb{G}(\mathbf{h}|2J^{2}q(0))\nu(0,x,\mathbf{h})=\frac{1}{\ell}\nu(0,x,0).$
(54)

This completes in principle our solution of the spin glass problem. One only
needs to compute the average magnetization and the time correlation function
in a single spin problem governed by $H_{\infty}$ or $\mathcal{S}_{\infty}$.
This yields, through integration of Eq. (52), the functions $q(x)$ and
$Q(\tau)$, which allow one in turn to compute the free energy, starting from
the partition function associated with $H_{\infty}$ or $\mathcal{S}_{\infty}$.
However, these results are formulated much more conveniently by introducing
the probability distribution of the local field.

####  1.3.4 Probability distribution of the local field

We introduce the probability distribution $\mathbb{P}(\mathbf{h})$ of the
local field acting on the effective impurity problem. Let us define the
function $P(x,\mathbf{h}|x^{\prime},\mathbf{h}^{\prime})$ as the solution of
the "backward" eq. (52) with boundary condition
$P(x^{\prime},\mathbf{h}|x^{\prime},\mathbf{h}^{\prime})=\delta(\mathbf{h}-\mathbf{h}^{\prime})$.
For any $x^{\prime}\geq x$, it has the Markovian property

$P(x,\mathbf{h}|1,\mathbf{h}^{\prime\prime})=\int\mathrm{d}\mathbf{h}^{\prime}P(x,\mathbf{h}|x^{\prime},\mathbf{h}^{\prime})P(x^{\prime},\mathbf{h}^{\prime}|1,\mathbf{h}^{\prime\prime}).$
(55)

Differentiating with respect to $x^{\prime}$ yields

$0=\int\mathrm{d}\mathbf{h}^{\prime}\left[\partial_{x^{\prime}}P(x,\mathbf{h}|x^{\prime},\mathbf{h}^{\prime})P(x^{\prime},\mathbf{h}^{\prime}|1,\mathbf{h}^{\prime\prime})+P(x,\mathbf{h}|x^{\prime},\mathbf{h}^{\prime})\partial_{x^{\prime}}P(x^{\prime},\mathbf{h}^{\prime}|1,\mathbf{h}^{\prime\prime})\right],$
(56)

and since $P$ satisfies (52), we obtain

$\begin{split}\int\mathrm{d}\mathbf{h}^{\prime}\,&P(x^{\prime},\mathbf{h}^{\prime}|1,\mathbf{h}^{\prime\prime})\partial_{x^{\prime}}P(x,\mathbf{h}|x^{\prime},\mathbf{h}^{\prime})=\\\
&-\int\mathrm{d}\mathbf{h}^{\prime}\,P(x,\mathbf{h}|x^{\prime},\mathbf{h}^{\prime})\frac{J^{2}}{2}\dot{q}(x^{\prime})\left(\bm{\nabla}_{\mathbf{h}^{\prime}}^{2}P(x^{\prime},\mathbf{h}^{\prime}|1,\mathbf{h}^{\prime\prime})+2\beta
x^{\prime}(\mathbf{s}\cdot\bm{\nabla}_{\mathbf{h}^{\prime}})P(x,\mathbf{h}^{\prime}|1,\mathbf{h}^{\prime\prime})\right).\end{split}$
(57)

Integrating by parts,

$\begin{split}\int\mathrm{d}\mathbf{h}^{\prime}\,&P(x^{\prime},\mathbf{h}^{\prime}|1,\mathbf{h}^{\prime\prime})\partial_{x^{\prime}}P(x,\mathbf{h}|x^{\prime},\mathbf{h}^{\prime})=\\\
&\int\mathrm{d}\mathbf{h}^{\prime}\,P(x^{\prime},\mathbf{h}^{\prime}|1,\mathbf{h}^{\prime\prime})\frac{J^{2}}{2}\dot{q}(x^{\prime})\left(\bm{\nabla}_{\mathbf{h}^{\prime}}^{2}P(x,\mathbf{h}|x^{\prime},\mathbf{h}^{\prime})-2\beta
x^{\prime}\bm{\nabla}_{\mathbf{h}^{\prime}}(\mathbf{s}(x^{\prime},\mathbf{h}^{\prime})\cdot
P(x,\mathbf{h}|x^{\prime},\mathbf{h}^{\prime})\right).\end{split}$ (58)

Since this is true for any $\mathbf{h}^{\prime\prime}$,
$P(x,\mathbf{h}|x^{\prime},\mathbf{h}^{\prime})$ is found to satisfy the
"forward" equation

$\frac{\partial P}{\partial
x^{\prime}}=\frac{J^{2}}{2}\dot{q}(x^{\prime})\left(\bm{\nabla}_{\mathbf{h}^{\prime}}^{2}P-2\beta
x^{\prime}\bm{\nabla}_{\mathbf{h}^{\prime}}(\mathbf{s}\cdot P\right)).$ (59)

We now rewrite eq. (53) as

$\langle
R\rangle_{\mathcal{S}}=\int\mathrm{d}\mathbf{h}^{\prime}\mathrm{d}\mathbf{h}\,\mathbb{G}(\mathbf{h}|2J^{2}q(0))P(0,\mathbf{h}|1,\mathbf{h}^{\prime})\langle
R\rangle_{\mathcal{S}_{1}(\mathbf{h}^{\prime})}\equiv\int\mathbb{P}(1,\mathbf{h}^{\prime})\langle
R\rangle_{\mathcal{S}_{1}(\mathbf{h}^{\prime})},$ (60)

having defined

$\mathbb{P}(x^{\prime},\mathbf{h}^{\prime})=\int\mathrm{d}\mathbf{h}\,\mathbb{G}(\mathbf{h}|2J^{2}q(0))P(0,\mathbf{h}|x^{\prime},\mathbf{h}^{\prime}).$
(61)

$\mathbb{P}$ then satisfies eq. (59) with boundary condition
$\mathbb{P}(0,\mathbf{h})=\mathbb{G}(\mathbf{h}|2J^{2}q(0))=\delta(\mathbf{h})$.
Similarly, for the Parisi order parameter we find

$q(x)=\frac{1}{\ell}\int\mathrm{d}\mathbf{h}\,\mathbb{P}(x,\mathbf{h})\mathbf{s}^{2}(x,\mathbf{h})$
(62)

####  1.3.5 Final formulation of the solution

We are now in position to formulate the solution in its final form.

We have reduced the solution of the fully connected Heisenberg spin glass
problem to the solution of a single spin problem. It is governed in the
classical case by the Hamiltonian
$H_{\infty}(\mathbf{h})=-\mathbf{h}\cdot\mathbf{S}$ and in the quantum case by
the action

$\mathcal{S}_{\infty}(\mathbf{h})=\mathcal{S}_{0}[\mathbf{S}(\tau)]-\int_{0}^{1}\mathrm{d}\tau\mathrm{d}\tau^{\prime}\frac{\beta^{2}J^{2}(Q(\tau-\tau^{\prime})-q(1))}{2}\mathbf{S}(\tau)\cdot\mathbf{S}(\tau^{\prime})-\beta\int_{0}^{1}\mathbf{h}\cdot\mathbf{S}(\tau).$
(63)

The Parisi order parameter and local susceptibility are obtained as

$q(x)=\frac{1}{\ell}\int\mathrm{d}\mathbf{h}\,\mathbb{P}(x,\mathbf{h})\mathbf{s}(x,\mathbf{h})^{2}~{}~{}~{}~{}\mathrm{and}~{}~{}~{}~{}~{}Q(\tau)=\int\mathrm{d}\mathbf{h}\,\mathbb{P}(1,\mathbf{h})\langle\mathrm{T}\mathbf{S}(\tau)\mathbf{S}(0)\rangle_{\infty,\mathbf{h}},$
(64)

where $\langle\cdot\rangle_{\infty,\mathbf{h}}$ denotes the mean value with
respect to $H_{\infty}(\mathbf{h})$ or $\mathcal{S}_{\infty}(\mathbf{h})$.
$\mathbb{P}(x,\mathbf{h})$ and $\mathbf{s}(x,\mathbf{h})$ are obtained by
solving

$\displaystyle\frac{\partial\mathbb{P}}{\partial
x}=\frac{J^{2}}{2}\frac{\mathrm{d}q}{\mathrm{d}x}\left(\bm{\nabla}^{2}\mathbb{P}-2\beta
x\bm{\nabla}(\mathbf{s}\cdot\mathbb{P})\right)~{}~{}\text{with}~{}~{}\mathbb{P}(0,\mathbf{h})=\delta(\mathbf{h});$
(65) $\displaystyle\frac{\partial\mathbf{s}}{\partial
x}=-\frac{J^{2}}{2}\frac{\mathrm{d}q}{\mathrm{d}x}\left(\bm{\nabla}^{2}\mathbf{s}+2\beta
x(\mathbf{s}\cdot\bm{\nabla})\mathbf{s}\right)~{}~{}\text{with}~{}~{}\mathbf{s}(1,\mathbf{h})=\langle\mathbf{S}\rangle_{\infty,\mathbf{h}}.$
(66)

Once $q(x)$ and $Q(\tau)$ are known, the free energy can be evaluated
according to Eqs. (42) and (47).

In practice, the local problem is solved trivially in the classical case:

$\langle S\rangle_{\infty,\mathbf{h}}=S\left(\frac{1}{\mathrm{tanh}(\beta
Sh)}-\frac{1}{\beta Sh}\right),$ (67)

and the self-consistency needs only to be established between Eqs. (64)
through (66). In the quantum case, the local problem is still a many-body
problem, which can only be solved numerically. Moreover, the local action
depends self-consistently on the observables $q(x)$ and $Q(\tau)$, so that
convergence must be reached for two interconnected self-consistency loops (see
Fig. 1 of the main text).

###  1.4 Internal energy

The internal energy per spin $U/N$ can in principle be obtained by numerically
differentiating the free energy with respect to temperature. However, it can
be obtained in a more straightforward way by carrying out the differentiation
analytically before replica symmetry breaking. Going back to Eqs. (9) and
(16), we compute the internal energy as

$\frac{U}{N}=-\frac{\partial(\beta f)}{\partial\beta}=-\lim_{n\to
0}\frac{1}{n}\frac{\partial\mathcal{F}[\\{Q_{ab}^{\star}\\}]}{\partial\beta}.$
(68)

Since the $Q_{ab}$ were introduced as integration variables, they may be
arbitrarily rescaled. If we define $\tilde{Q}_{ab}=\beta^{2}J^{2}Q_{ab}$, the
"complicated" terms in Eq. (8) and (15) no longer have an explicit $\beta$
dependence. We then obtain, in the classical case,

$\frac{\partial\mathcal{F}_{\rm Cl}}{\partial\beta}=-\frac{n\beta
J^{2}S^{4}}{2\ell}-\frac{\ell}{\beta^{3}J^{2}}\sum_{a<b}\tilde{Q}_{ab}^{2}=-\frac{n\beta
J^{2}S^{4}}{2\ell}-\ell\beta J^{2}\sum_{a<b}Q_{ab}^{2}.$ (69)

We then use the property

$\lim_{n\to
0}\sum_{a<b}Q_{ab}^{2}=-\frac{1}{2}\int_{0}^{1}\mathrm{d}x\,q(x)^{2},$ (70)

to obtain

$\frac{U_{\rm Cl}}{N}=-\frac{\beta
J^{2}S^{4}}{2\ell}\left(1-\frac{\ell^{2}}{S^{4}}\int_{0}^{1}\mathrm{d}x\,q(x)^{2}\right).$
(71)

Similarly, in the quantum case,

$\begin{split}\frac{\partial\mathcal{F}_{\rm
Q}}{\partial\beta}&=-\frac{n\ell}{2\beta^{3}J^{2}}\int_{0}^{1}\mathrm{d}\tau\mathrm{d}\tau^{\prime}\tilde{Q}_{aa}(\tau,\tau^{\prime})^{2}-\frac{\ell}{\beta^{3}J^{2}}\sum_{a<b}\tilde{Q}_{ab}^{2}\\\
&=-\frac{n\ell\beta
J^{2}}{2}\int_{0}^{1}\mathrm{d}\tau\mathrm{d}\tau^{\prime}Q_{aa}(\tau,\tau^{\prime})^{2}-\ell\beta
J^{2}\sum_{a<b}Q_{ab}^{2},\end{split}$ (72)

so that

$\frac{U_{\rm Q}}{N}=-\frac{\ell\beta
J^{2}}{2}\left(\int_{0}^{1}\mathrm{d}\tau\,Q(\tau)^{2}-\int_{0}^{1}\mathrm{d}x\,q(x)^{2}\right).$
(73)

###  1.5 Doped quantum spin glass

Our solution is readily generalized to a "doped" Heisenberg spin glass model.
Physically, the randomly-interacting spins are now being carried by electrons,
which can hop between the sites of a fully connected lattice and experience an
on-site repulsion $U$. The corresponding Hamiltonian is

$H=H_{0}+U\sum_{i}n_{i\uparrow}n_{i\downarrow}-\sum_{i<j}J_{ij}\overline{\mathbf{S}}_{i}\cdot\overline{\mathbf{S}}_{j},$
(74)

with
$H_{0}=-\sum_{ij,\sigma}(t_{ij}+\mu\delta_{ij})c^{\dagger}_{i\sigma}c_{j\sigma}$.
Here the $c_{i\sigma}$ and $c^{\dagger}_{i\sigma}$ are the electronic
annihilation and creation operators, $\sigma=\uparrow,\downarrow$ labels the
spin state, $n_{i\sigma}=c^{\dagger}_{i\sigma}c_{i\sigma}$ is the occupation
number on site $i$, $U$ is the on-site electron-electron interaction, $\mu$ is
the chemical potential and the $t_{ij}$ are hopping amplitudes, randomly
distributed with variance $t^{2}/N$. The spins are defined in terms of the
fermionic operators as
$S_{i}^{\alpha}=(1/2)\sum_{\sigma,\sigma^{\prime}}c_{i\sigma}^{\dagger}\bm{\sigma}^{\alpha}_{\sigma\sigma^{\prime}}c_{i\sigma^{\prime}}$,
with $\bm{\sigma}^{\alpha}$ the Pauli matrices.

The solution of this model proceeds along the same steps as above, and here we
only outline the main differences. We represent the partition function as a
path integral over the Grassman variables
$(c_{i\sigma}^{\dagger}(\tau),c_{i\sigma}(\tau))$:

$Z=\int\prod_{i}[Dc_{i\sigma}^{\dagger}(\tau)][Dc_{i\sigma}(\tau)]e^{-\int_{0}^{\beta}\mathrm{d}\tau\left[\sum_{i,\sigma}c^{\dagger}_{i\sigma}(\tau)\partial_{\tau}c_{i\sigma}(\tau)+H(\tau)\right]}.$
(75)

Note that we now let the imaginary time run between 0 and $\beta$. We then
compute $\overline{Z^{n}}$, with the disorder average now being carried out
over both the couplings $J_{ij}$ and the hoppings $t_{ij}$. The Hubbard-
Stratonovitch transformation now introduces integration variables $Q_{ab}$ and
$\Delta_{ab}$, that decouple respectively the spin-spin interaction term and
the hopping term. Upon performing the relevant simplifications (for instance,
the $\Delta_{a\neq b}$ vanish), and dropping indices for simplicity,

$\overline{Z^{n}}\propto\int[DQ][D\Delta]\exp\left[-N\left\\{\int_{0}^{\beta}\mathrm{d}\tau\mathrm{d}\tau^{\prime}\left(\frac{J^{2}(Q^{ab})^{2}}{4}+\frac{(\Delta^{a})^{2}}{t^{2}}\right)-\log\int[Dc]e^{-\mathcal{S}_{\rm
loc}[Q,\Delta]}\right\\}\right],$ (76)

with

$\begin{split}&\mathcal{S}_{\rm
loc}=\int_{0}^{\beta}\mathrm{d}\tau\sum_{a\sigma}[c_{\sigma}^{a\dagger}(\tau)(\partial_{\tau}-\mu)c_{\sigma}^{a}(\tau)+Un^{a}_{\uparrow}(\tau)n^{a}_{\downarrow}(\tau)]+\dots\\\
&+\int_{0}^{\beta}\mathrm{d}\tau\mathrm{d}\tau^{\prime}\sum_{a}\left[\Delta(\tau^{\prime}-\tau)\sum_{\sigma}c^{a\dagger}_{\sigma}(\tau^{\prime})c^{a}_{\sigma}(\tau)-\frac{J^{2}Q(\tau^{\prime}-\tau)}{2}\mathbf{S}^{a}(\tau)\cdot\mathbf{S}^{a}(\tau^{\prime})\right]-J^{2}\sum_{a<b}Q^{ab}\mathbf{S}^{a}(\tau)\cdot\mathbf{S}^{b}(\tau^{\prime}).\end{split}$
(77)

The saddle point equations read

$\displaystyle\Delta(\tau)=t^{2}G(\tau),~{}~{}~{}G(\tau)=-\langle\mathrm{T}c^{a}_{\sigma}(\tau)c^{\dagger
a}_{\sigma}(0)\rangle;$ (78) $\displaystyle
Q(\tau-\tau^{\prime})=\frac{1}{\ell}\langle\mathrm{T}\mathbf{S}^{a}(\tau)\cdot\mathbf{S}^{a}(\tau^{\prime})\rangle_{\mathcal{S}_{\rm
loc}};$ (79) $\displaystyle
Q_{ab}=\frac{1}{\ell}\langle\mathbf{S}^{a}(\tau)\cdot\mathbf{S}^{b}(\tau)\rangle_{\mathcal{S}_{\rm
loc}}.$ (80)

Upon replica symmetry breaking, the original lattice problem is reduced to a
self-consistent impurity problem, governed by the action

$\begin{split}&\mathcal{S}_{\infty}(\mathbf{h})=\int_{0}^{\beta}\mathrm{d}\tau\sum_{\sigma}[c_{\sigma}^{\dagger}(\tau)(\partial_{\tau}-\mu)c_{\sigma}(\tau)+Un_{\uparrow}(\tau)n_{\downarrow}(\tau)]-\int_{0}^{\beta}\mathrm{d}\tau\mathbf{h}\,\cdot\mathbf{S}(\tau)\\\
&+\int_{0}^{\beta}\mathrm{d}\tau\mathrm{d}\tau^{\prime}\left[\Delta(\tau^{\prime}-\tau)\sum_{\sigma}c^{\dagger}_{\sigma}(\tau^{\prime})c_{\sigma}(\tau)-\frac{J^{2}(Q(\tau^{\prime}-\tau)-q(1))}{2}\mathbf{S}(\tau)\cdot\mathbf{S}(\tau^{\prime})\right].\\\
\end{split}$ (81)

The self-consistency conditions read

$q(x)=\frac{1}{\ell}\int\mathrm{d}\mathbf{h}\,\mathbb{P}(x,\mathbf{h})\mathbf{s}(x,\mathbf{h})^{2},~{}~{}~{}~{}~{}Q(\tau)=\int\mathrm{d}\mathbf{h}\,\mathbb{P}(1,\mathbf{h})\langle\mathrm{T}\mathbf{S}(\tau)\mathbf{S}(0)\rangle_{\infty,\mathbf{h}},$
(82)

and

$\Delta(\tau)=t^{2}\int\mathrm{d}\mathbf{h}\,\mathbb{P}(x,\mathbf{h})\frac{1}{2}\sum_{\sigma}G_{\sigma}(\tau,\mathbf{h}),~{}~{}~{}\mathrm{with}~{}~{}~{}G_{\sigma}(\tau)=-\langle\mathrm{T}c_{\sigma}(\tau)c^{\dagger}_{\sigma}(0)\rangle_{\infty,\mathbf{h}}.$
(83)

$\mathbb{P}(x,\mathbf{h})$ and $\mathbf{s}(x,\mathbf{h})$ are obtained by
solving Eqs. (65) and (66).

We note that we have used in Eq. (83) the spin-symmetrized Green’s function
$\overline{G}(\tau)=(G_{\uparrow}(\tau)+G_{\downarrow}(\tau))/2$. In the
absence of an external magnetic field, it is equal to the Green’s function of
either spin. In the presence of a field $\mathbf{h}$, $G_{\uparrow}$ and
$G_{\downarrow}$ depend on both the magnitude and direction of the
$\mathbf{h}$, but $\overline{G}$ depends only on $||\mathbf{h}||$.

Indeed, let us apply a unitary transformation $U$ to the basis of spin states
$\\{|\sigma\rangle\\}$. This defines new basis vectors
$\\{|\sigma^{\prime}\rangle\\}$. We may define creation and annihilation
operators for spin up and spin down electrons in this new basis. They are
given by

$c^{\dagger}_{\sigma^{\prime}}=\sum_{\sigma=\uparrow,\downarrow}U_{\sigma^{\prime}\sigma}c^{\dagger}_{\sigma}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\text{and}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}c_{\sigma^{\prime}}=\sum_{\sigma=\uparrow,\downarrow}U^{*}_{\sigma^{\prime}\sigma}c_{\sigma}$
(84)

The spin up Green’s function in the new basis is

$\displaystyle G^{\prime}_{\uparrow}(\tau)$
$\displaystyle=-\langle(U^{*}_{\uparrow\uparrow}c_{\uparrow}(\tau)+U_{\uparrow\downarrow}^{*}c_{\downarrow}(\tau))(U_{\uparrow\uparrow}c^{\dagger}_{\uparrow}(0)+U_{\uparrow\downarrow}c^{\dagger}_{\downarrow}(0))\rangle$
(85)
$\displaystyle=|U_{\uparrow\uparrow}|^{2}G_{\uparrow}(\tau)+|U_{\uparrow\downarrow}|^{2}G_{\downarrow}(\tau)-U^{*}_{\uparrow\uparrow}U_{\uparrow\downarrow}\langle
c_{\uparrow}(\tau)c^{\dagger}_{\downarrow}(0)\rangle-
U_{\uparrow\uparrow}U_{\uparrow\downarrow}^{*}\langle
c_{\downarrow}(\tau)c^{\dagger}_{\uparrow}(0)\rangle$ (86)

Similarly,

$\displaystyle G^{\prime}_{\downarrow}(\tau)$
$\displaystyle=-\langle(U^{*}_{\downarrow\uparrow}c_{\uparrow}(\tau)+U_{\downarrow\downarrow}^{*}c_{\downarrow}(\tau))(U_{\downarrow\uparrow}c^{\dagger}_{\uparrow}(0)+U_{\downarrow\downarrow}c^{\dagger}_{\downarrow}(0))\rangle$
(87)
$\displaystyle=|U_{\downarrow\uparrow}|^{2}G_{\uparrow}(\tau)+|U_{\downarrow\downarrow}|^{2}G_{\downarrow}(\tau)-U^{*}_{\downarrow\uparrow}U_{\downarrow\downarrow}\langle
c_{\uparrow}(\tau)c^{\dagger}_{\downarrow}(0)\rangle-
U_{\downarrow\uparrow}U_{\downarrow\downarrow}^{*}\langle
c_{\downarrow}(\tau)c^{\dagger}_{\uparrow}(0)\rangle$ (88)

Now, the unitarity of $U$ imposes that
$|U_{\uparrow\uparrow}|^{2}+|U_{\uparrow\downarrow}|^{2}=1$ and
$U_{\uparrow\uparrow}U_{\uparrow\downarrow}^{*}+U_{\downarrow\uparrow}U_{\downarrow\downarrow}^{*}=0$,
so that

$G^{\prime}_{\uparrow}(\tau)+G^{\prime}_{\downarrow}(\tau)=G_{\uparrow}(\tau)+G_{\downarrow}(\tau).$
(89)

The spin-symmetrized Green’s function is thus rotationally-invariant even
under an external magnetic field.

##  2 Numerical procedures

###  2.1 Solution of the Parisi equations

####  2.1.1 Integral equation formulation

The Parisi equations Eqs. (65) and (66) are non-linear partial differential
equations that need to be solved numerically. We introduce the Green’s
function $G(x,\mathbf{h}|x^{\prime},\mathbf{h}^{\prime})$ of their linear
part, which satisfies

$\left\\{\begin{array}[]{lr}\dfrac{\partial G}{\partial
x}=\dfrac{J^{2}}{2}\dfrac{\mathrm{d}q}{\mathrm{d}x}\bm{\nabla}^{2}G&\text{for}~{}x\geq
x^{\prime}\\\ &\\\ \dfrac{\partial G}{\partial
x}=-\dfrac{J^{2}}{2}\dfrac{\mathrm{d}q}{\mathrm{d}x}\bm{\nabla}^{2}G&\text{for}~{}x\leq
x^{\prime}\end{array}\right..$ (90)

with
$G(x^{\prime},\mathbf{h}|x^{\prime},\mathbf{h}^{\prime})=\delta(\mathbf{h}-\mathbf{h}^{\prime})$.
$G$ is simply a Gaussian:

$G(x,\mathbf{h}|x^{\prime},\mathbf{h}^{\prime})=\frac{1}{[2\pi
J^{2}|q(x)-q(x^{\prime})|]^{3/2}}\exp\left({-\frac{(\mathbf{h}-\mathbf{h}^{\prime})^{2}}{2J^{2}|q(x)-q(x^{\prime})|}}\right).$
(91)

If we now treat the non-linear terms as source terms, Eqs. (65) and (66) can
be formally integrated according to

$\mathbb{P}(x,\mathbf{h})=\frac{1}{(2\pi
J^{2}q(0))^{3/2}}e^{-h^{2}/2J^{2}q(0)}-\beta
J^{2}\int_{0}^{x}\mathrm{d}x^{\prime}\,x^{\prime}\dot{q}(x^{\prime})\int\mathrm{d}\mathbf{h}^{\prime}\,G(x,\mathbf{h}|x^{\prime},\mathbf{h}^{\prime})\bm{\nabla}(\mathbf{s}\cdot\mathbb{P})|_{x^{\prime},\mathbf{h}^{\prime}}$
(92)

and

$\mathbf{s}(x,\mathbf{h})=\int\mathrm{d}\mathbf{h}^{\prime}G(x,\mathbf{h}|1,\mathbf{h}^{\prime})\langle\mathbf{S}\rangle_{\infty,\mathbf{h}^{\prime}}+\beta
J^{2}\int_{x}^{1}\mathrm{d}x^{\prime}\,x^{\prime}\dot{q}(x^{\prime})\int\mathrm{d}\mathbf{h}^{\prime}\,G(x,\mathbf{h}|x^{\prime},\mathbf{h}^{\prime})(\mathbf{s}\cdot\bm{\nabla})\mathbf{s}|_{x^{\prime},\mathbf{h}^{\prime}}$
(93)

We now use the fact that $\mathbf{s}(x,\mathbf{h})$ is always in the direction
of $\mathbf{h}$: $\mathbf{s}(x,\mathbf{h})=s(x,h)\hat{\mathbf{h}}$. Then,

$\bm{\nabla}(\mathbf{s}\cdot\mathbb{P})|_{x^{\prime},\mathbf{h}^{\prime}}=\frac{1}{h^{\prime
2}}\frac{\partial(h^{\prime 2}s\mathbb{P})}{\partial h^{\prime}}$ (94)

and

$(\mathbf{s}\cdot\bm{\nabla})\mathbf{s}|_{x^{\prime},\mathbf{h}^{\prime}}=s\frac{\partial
s}{\partial h^{\prime}}\hat{\mathbf{h}}^{\prime}.$ (95)

Hence, we may introduce angle-integrated Green’s functions:

$\begin{split}\overline{G}(x,h&|x^{\prime},h^{\prime})=(h^{\prime})^{2}\int\mathrm{d}\hat{\mathbf{h}}^{\prime}\,G(x,\mathbf{h}|x^{\prime},\mathbf{h}^{\prime})\\\
&=\frac{1}{[2\pi
J^{2}|q(x)-q(x^{\prime})|]^{1/2}}\frac{hh^{\prime}}{h^{2}}\left[\exp\left({-\frac{(h-h^{\prime})^{2}}{2J^{2}|q(x)-q(x^{\prime})|}}\right)-\exp\left({-\frac{(h+h^{\prime})^{2}}{2J^{2}|q(x)-q(x^{\prime})|}}\right)\right];\end{split}$
(96)
$\begin{split}\overline{G}_{s}(x,h&|x^{\prime},h^{\prime})=(h^{\prime})^{2}\int\mathrm{d}\hat{\mathbf{h}}^{\prime}\,G(x,\mathbf{h}|x^{\prime},\mathbf{h}^{\prime})(\hat{\mathbf{h}}^{\prime}\cdot\hat{\mathbf{h}})\\\
&=\frac{1}{[2\pi
J^{2}|q(x)-q(x^{\prime})|]^{1/2}}\frac{1}{h^{2}}\left[(J^{2}|q(x)-q(x^{\prime})|+hh^{\prime})\exp\left({-\frac{(h+h^{\prime})^{2}}{2J^{2}|q(x)-q(x^{\prime})|}}\right)\right.\\\
&\left.+(hh^{\prime}-J^{2}|q(x)-q(x^{\prime})|)\exp\left({-\frac{(h-h^{\prime})^{2}}{2J^{2}|q(x)-q(x^{\prime})|}}\right)\right].\end{split}$
(97)

Eqs. (92) and (93) then reduce to one-dimensional integral equations:

$\mathbb{P}(x,h)=\frac{1}{(2\pi J^{2}q(x))^{3/2}}e^{-h^{2}/2J^{2}q(x)}-\beta
J^{2}\int_{0}^{x}\mathrm{d}x^{\prime}\,x^{\prime}\dot{q}(x^{\prime})\int\mathrm{d}h^{\prime}\,\overline{G}(x,h|x^{\prime},h^{\prime})\frac{1}{h^{\prime
2}}\frac{\partial(h^{\prime 2}s\mathbb{P})}{\partial h^{\prime}}$ (98)

and

$s(x,h)=\int\mathrm{d}h^{\prime}\,\overline{G}_{s}(x,h|1,h^{\prime})\langle
S\rangle_{\infty,\mathbf{h}^{\prime}}+\beta
J^{2}\int_{x}^{1}\mathrm{d}x^{\prime}\,x^{\prime}\dot{q}(x^{\prime})\int\mathrm{d}h^{\prime}\,\overline{G}_{s}(x,h|x^{\prime},h^{\prime})\left.s\frac{\partial
s}{\partial h^{\prime}}\right|_{x^{\prime},h^{\prime}}.$ (99)

####  2.1.2 Numerical procedure

The variable $x$ is discretized, and Eqs. (98) and (99) are solved on
successive $x$ points. We typically use $128$ to $256$ $x$ points in the
region where $q(x)$ is non-constant. When $q(x)$ and $q(x^{\prime})$ are very
close, the Green’s functions $\overline{G}(x,h|x^{\prime},h^{\prime})$ are
strongly peaked around $h=h^{\prime}$. However, the result of the convolution
varies much more slowly as a function of $h$. Therefore, we sample
$h^{\prime}$ on a very fine grid (typically 20000 points) and $h$ on a much
coarser Chebyshev grid (typically 256 points), and then use barycentric
interpolation to obtain the result of the convolution on the fine grid for the
next integration step.

When the non-linear term is large (typically, when $q(x)$ is large), the
numerical integration procedure is unstable. Indeed, at step $n$, the non-
linear term contains a contribution from the $n^{\rm th}$ numerical derivative
of the initial condition, which develops an oscillatory instability with a
frequency on the order of the Chebyshev grid spacing. As long as the grid is
much finer than the structure in $\mathbb{P}(x,h)$ and $s(x,h)$, the spurious
oscillations may be filtered out. We use Savitzky-Golay filtering, which
amounts to splitting the $h$ interval into panels, and fitting
$\mathbb{P}(x,h)$ or $s(x,h)$ with a second degree polynomial in $h$ on each
of the panels. The panel size is determined adaptively, as a fraction of the
value of $h$ where the non-linear term has a maximum. With this procedure, and
given our level of discretization, the filtered-out part is $\lesssim 10^{-5}$
of the remainder at every $x$ step, which roughly sets the precision level of
our solver. We further check that the normalization $\int_{0}^{\infty}4\pi
h^{2}\mathbb{P}(x,h)\mathrm{d}h=1$ and the sum rule [8]

$1-\frac{\ell}{S^{2}}\int_{0}^{1}\mathrm{d}x\,q(x)=T$ (100)

in the classical case are satisfied to a precision better than $10^{-4}$ after
self-consistency has been reached.

###  2.2 Solution of the quantum impurity problem

We solve the quantum impurity problem defined by the action in Eq. (81), and
its restriction in Eq. (63), by using a continuous-time Monte Carlo procedure
based on a hybridization expansion in the segment picture (CTSEG algorithm)
[9, 10]. Briefly, the algorithm is based on the double expansion of the
hybridization ($\Delta$) and spin-spin interaction ($Q$) terms. A
configuration is defined by the order of the expansion and by the choice of
points in imaginary time. Several possible updates or "moves" are defined for
the configuration, which are accepted or rejected according to the Metropolis
prescription. Further details of our implementation will be described
elsewhere.

###  2.3 General solution procedure

####  2.3.1 Classical case

In the classical case, the magnetization
$\langle\mathbf{S}\rangle_{\infty,\mathbf{h}}$ is known analytically and does
not depend on the spin glass parameters. Therefore, we start from an initial
guess for $q(x)$ and iteratively solve Eqs. (66) and (65) until convergence.
As self-consistency is approached, $q(x)$ develops a singularity: $q(x)=q(1)$
for $x>x_{c}$. Yet, our numerical solution cannot produce a singularity in the
mathematical sense. Two procedures can be envisioned for finding the
breakpoint $x_{c}$.

  1. 1.

We may apply a threshold to $\mathrm{d}q/\mathrm{d}x$. Concretely, once a new
$q(x)$ has been determined from Eq. (64), we enforce that if
$q(1)-q(x)<\epsilon$, $q(x)=q(1)$, with typically $\epsilon=10^{-5}$, on the
order of our solver’s accuracy. This will result in the breakpoint moving as
the solution is iterated, and eventually converging.

  2. 2.

We may iterate the equations until convergence for a range of fixed values of
$x_{c}$, and then choose the $x_{c}$ for which the sum rule in Eq. (100) is
most accurately satisfied.

We found that procedure 1 converges in a reasonable number of iterations at
low temperatures $T/J\leq 0.04$, while procedure 2 is required for higher
temperatures. At zero temperature, since $x_{c}/T$ logarithmically diverges as
$T\to 0$, there is no breakpoint in the finite interval of $x/T$ in which we
solve the equations.

####  2.3.2 Quantum case

In the quantum case, the magnetization
$\langle\mathbf{S}\rangle_{\infty,\mathbf{h}}$ is obtained by the solving the
action in Eq. (81) or (63) with quantum Monte Carlo for a range of values of
$h$. Furthermore, the action depends on the spin glass parameters and must be
updated until self-consistency is reached. We start from an initial guess for
$q(x)$ and $Q(\tau)$ (and $\Delta(\tau)$ in the doped case), and determine
$\langle\mathbf{S}\rangle_{\infty,\mathbf{h}}$ (as well as the spin-spin
correlation function, and possibly the electron Green’s function) on the
impurity, for typically 30 values of $h$ sampled on a Chebyshev grid. The
observables at arbitrary values of $h$ are then inferred by barycentric
interpolation. Then, given $\langle\mathbf{S}\rangle_{\infty,\mathbf{h}}$,
Eqs. (66) and (65) are iterated $N$ times. An arbitrarily large value of $N$
is not necessarily best for rapid convergence of the whole self-consistent
procedure; we found empirically that $N=5$ is adequate. Procedure 1 (defined
above) was sufficient to obtain reasonable convergence of the breakpoint. Once
$q(x)$ and $\mathbb{P}(h)$ are obtained, we determine the $Q(\tau)$ and
$\Delta(\tau)$ that define a new impurity action, and repeat the procedure
until convergence of $q(x)$, $Q(\tau)$ and $\Delta(\tau)$.

The Parisi equations (66) and (65) can be solved directly at zero temperature,
but not the quantum impurity problem. However, the input for the Parisi
equations is only the magnetization, which we find to be almost converged with
temperature for $\beta J>100$ (Fig. S5). We use the magnetization obtained at
the lowest accessible temperature $\beta J=200$ as an input for the zero-
temperature Parisi equations.

##  3 Supplementary figures

Figure 1: Parisi order parameter $q(x)$ for the classical Heisenberg spin
glass, for a range of temperature values. Figure 2: Local field probability
distribution $\mathbb{P}(h)$ for the classical Heisenberg spin glass, for a
range of temperature values. Figure 3: Parisi order parameter $q(x)$ for the
quantum Heisenberg spin glass, for a range of temperature values. We did not
attempt to reach convergence for the breakpoint at high temperature ($\beta
J<50$), as its position has a negligible influence on the other properties of
the spin glass (particularly, the scaling of $Q(\tau)$).  Figure 4: Local
field probability distribution $\mathbb{P}(h)$ for the quantum Heisenberg spin
glass, for a range of temperature values. Figure 5: Magnetization $\langle
S\rangle_{\infty,\mathbf{h}}$ of the quantum Heisenberg spin glass, for a
range of temperature values. Figure 6: Magnetization $\langle
S\rangle_{\infty,\mathbf{h}}$ versus temperature, for three magnetic field
values. At the lowest accessible temperatures, the magnetization is almost
converged to its zero temperature shape. Figure 7: Rescaled spin
susceptibility in the quantum spin glass phase $\chi(\tau)/\chi(\beta/2)$, for
a range of temperature values. It is well described by the conformal scaling
form $\chi(\tau)/\chi(\beta/2)=1/\sin(\pi\tau/\beta)^{\theta}$. The scaling
functions are plotted for $\theta=1$ and $\theta=3/2$ (black lines). Figure 8:
Parisi order parameter $q(x)$ for the doped Heisenberg spin glass at $\beta
J=50$, for a range of values of the doping $p$. Figure 9: Imaginary part of
the electron Green’s function for the doped Heisenberg spin glass at $\beta
J=50$, for a range of values of the doping $p$. Figure 10: Rescaled spin
susceptibility in the metallic spin glass phase $\chi(\tau)/\chi(\beta/2)$, at
$\beta J=50$, for a range of values of the doping $p$. It is reasonably
described by the conformal scaling form
$\chi(\tau)/\chi(\beta/2)=1/\sin(\pi\tau/\beta)^{\theta}$. The scaling
functions are plotted for $\theta=1$ and $\theta=3/2$ (black lines). Figure
11: Scaling exponent $\theta$, as a function of temperature in the metallic
spin glass phase, for a range of values of the doping $p$. Figure 12:
Breakpoint $x_{c}(T)$ vs $T$. $x_{c}(T)/T$ is a characteristic inverse energy
scale governing the distribution of free energies of excited states. It scales
logarithmically with $T$ in both the classical and the quantum case. Dashed
lines are guides to the eye. Inset: definition of the breakpoint $x_{c}$
($q(x>x_{c})=q(1)$).

## References

  * [1] Sherrington, D. & Kirkpatrick, S.  Solvable model of a spin-glass.  _Physical Review Letters_ 35, 1792–1796 (1975).
  * [2] Bray, A. J. & Moore, M. A.  Replica theory of quantum spin glasses.  _Journal of Physics C: Solid State Physics_ 13, L655 (1980).
  * [3] Parisi, G.  A sequence of approximated solutions to the S–K model for spin glasses.  _Journal of Physics A: Mathematical and General_ 13, L115–L121 (1980).
  * [4] Duplantier, B.  Comment on Parisi’s equation for the SK model for spin glasses.  _Journal of Physics A: General Physics_ 14, 283–285 (1981).
  * [5] Andreanov, A. & Müller, M.  Long-range quantum ising spin glasses at T=0: Gapless collective excitations and universality.  _Physical Review Letters_ 109, 177201 (2012).
  * [6] Kiss, A., Zaránd, G. & Lovas, I.  Exact Solution for the Transverse Field Sherrington-Kirkpatrick Spin Glass Model with Continuous-Time Quantum Monte Carlo Method (2023).  ArXiv:2306.07337 [cond-mat].
  * [7] Gabay, M., Garel, T. & De Dominicis, C.  Symmetry breaking a la Parisi in the n-component SK model of a spin glass.  _Journal of Physics C: Solid State Physics_ 15, 7165–7180 (1982).
  * [8] Crisanti, A. & Rizzo, T.  Analysis of the infinite-replica symmetry breaking solution of the Sherrington-Kirkpatrick model.  _Physical Review E_ 65, 046137 (2002).
  * [9] Otsuki, J.  Spin-boson coupling in continuous-time quantum Monte Carlo.  _Physical Review B - Condensed Matter and Materials Physics_ 87, 125102 (2013).
  * [10] Ayral, T.  PhD thesis. Nonlocal Coulomb interactions and electronic correlations: novel many-body approaches. Chapter 11. (2015).


# Evidence of mini-jet emission in a large emission zone from a magnetically-
dominated gamma-ray burst jet

S.-X. Yi1 Email<EMAIL_ADDRESS>C.-W. Wang1,2 E-mail<EMAIL_ADDRESS>B.
Zhang3,4 E-mail<EMAIL_ADDRESS>S.-L. Xiong1 E-mail<EMAIL_ADDRESS>S.-N. Zhang1,2 W.-J. Tan1,2 J.-C. Liu1,2 W.-C. Xue1,2 Y.-Q. Zhang1,2 C.
Zheng1,2 R. Moradi1 Y. Wang5,6,7 P. Zhang1,8 Z.-H. An1 C. Cai9 P.-Y. Feng1,2
K. Gong1 D.-Y. Guo1 Y. Huang1 B. Li1 X.-B. Li1 X.-Q. Li1 X.-J. Liu1 Y.-Q. Liu1
X. Ma1 W.-X. Peng1 R. Qiao1 L.-M. Song1 J. Wang1 P. Wang1 Y. Wang1,2 X.-Y.
Wen1 S. Xiao10 Y.-B. Xu1 S. Yang1 Q.-B. Yi1,11 D.-L. Zhang1 F. Zhang1 H.-M.
Zhang1 J.-P. Zhang1,2 Z. Zhang1 X.-Y. Zhao1 Y. Zhao1,12 S.-J. Zheng1

###### Abstract

The origin of the prompt emission of gamma-ray bursts (GRBs) is still subject
to debate because of the not-well-constrained jet composition, location of the
emission region, and mechanism with which $\gamma$-rays are produced[1]. For
the bursts whose emission is dominated by non-thermal radiation, two leading
paradigms are internal shock model invoking collisions of matter-dominated
shells[2] and internal-collision-induced magnetic reconnection and turbulence
(ICMART) model invoking collisions of magnetically-dominated shells[3]. These
two models invoke different emission regions[2, 3] and have distinct
predictions on the origin of light curve variability[4, 5] and spectral
evolution[6, 7]. The second brightest GRB in history, GRB230307A[8, 9, 10,
11], provides an ideal laboratory to study the details of GRB prompt emission
thanks to its extraordinarily high photon statistics and the single broad
pulse shape characterized by an energy-dependent fast-rise-exponential-decay
(FRED) profile[12, 13, 14, 15, 16]. Here we demonstrate that its broad pulse
is composed of many rapidly variable short pulses, rather than being the
superposition of many short pulses on top of a slow component. Such a feature
is consistent with the ICMART picture, which envisages many mini-jets due to
local magnetic reconnection events in a large emission zone far from the GRB
central engine[3, 5, 17, 18, 19, 20], but raises a great challenge to the
internal shock model that attributes fast and slow variability components to
shocks at different radii with the emission being the superposition of various
components[21]. The results provide strong evidence for a Poynting-flux-
dominated jet composition of this bright GRB.

Key Laboratory of Particle Astrophysics, Institute of High Energy Physics,
Chinese Academy of Sciences, Beijing 100049, China.

University of Chinese Academy of Sciences, Chinese Academy of Sciences,
Beijing 100049, China.

Nevada Center for Astrophysics, University of Nevada Las Vegas, NV 89154, USA.

Department of Physics and Astronomy, University of Nevada Las Vegas, NV 89154,
USA.

ICRANet, Piazza della Repubblica 10, I-65122 Pescara, Italy

ICRA, Dipartamento di Fisica, Sapienza Università di Roma, Piazzale Aldo Moro
5, I-00185 Rome, Italy

INAF, Osservatorio Astronomico d’Abruzzo, Via M. Maggini snc, I-64100, Teramo,
Italy

College of Electronic and Information Engineering, Tongji University, Shanghai
201804, China

College of Physics and Hebei Key Laboratory of Photophysics Research and
Application, Hebei Normal University, Shijiazhuang, Hebei 050024, China

Guizhou Provincial Key Laboratory of Radio Astronomy and Data Processing,
Guizhou Normal University, Guiyang 550001, China

School of Physics and Optoelectronics, Xiangtan University, Yuhu District,
Xiangtan, Hunan, 411105, China

School of Computer and Information, Dezhou University, Dezhou 253023, China

The extremely bright GRB 230307A was firstly reported by the Gravitational
wave high-energy electromagnetic counterpart all-sky monitor (GECAM)[9] with
trigger time of 15:44:06.650 UT on 7 March 2023 (denoted as T0). The Lobster
Eye Imager for Astronomy (LEIA, the pathfinder of the Einstein Probe mission)
also observed its prompt emission in the soft X-ray band[11]. While the
burst’s duration (measured as $T_{90}$) of about 41 seconds aligns with the
category of long-duration GRBs, the association of a kilonova[10] and its
unique properties in the prompt emission strongly suggest its origin from a
compact binary merger event[11].

Despite of the extraordinary brightness, unsaturated data record by GECAM
enable us to accurately characterize the temporal properties of GRB 230307A
(see Methods). The GECAM time-binned net count rate as a function of time in
different energy bands (multi-band light curves hereafter) show a typical FRED
profile (left panel of Figure 1). Therefore, we use a conventional
parameterization of the FRED profile [22] (“Norris05" hereafter)

$L(t)\propto\frac{1}{\exp(\frac{\tau_{\rm{r}}}{t-t_{s}}+\frac{t-t_{s}}{\tau_{\rm{d}}})},$
(1)

to fit the light curves in each energy bands. In the above formulation,
$t_{\rm{s}}$ is the starting instance of the pulse, $\tau_{\rm{r}}$ and
$\tau_{\rm{d}}$ are the rising and decaying time scales, respectively. The
fitting results are listed in Table 1. The peak time $t_{\rm{p}}$ and the
width $w$ of the pulse are therefore defined as:
$t_{\rm{p}}=t_{\rm{s}}+\sqrt{\tau_{\rm{r}}\tau_{\rm{d}}}$ and
$w=\tau_{\rm{r}}+\tau_{\rm{d}}$. It can be seen in figure 2 that there is a
clear energy dependence in both $w$ and $\tilde{t}_{\rm{p}}\equiv
t_{\rm{p}}-t_{\rm{s}}$, as commonly observed in GRB broad pulses. In the
energy range from 6 keV to 300 keV, the $w-E$ and $\tilde{t}_{\rm{p}}-E$
relations are both in a power law function, with the power index
-0.36${}^{+0.05}_{-0.06}$ and -0.37${}^{+0.05}_{-0.06}$, respectively. Above
$\sim 300$ keV, the energy dependence on both $w$ and $\tilde{t}_{\rm{p}}$
saturate. This highly correlated $w-E$ and $\tilde{t}_{\rm{p}}-E$ relationship
implies self-similar profiles across different energy bands, a feature that
has also been reported for other bursts in the literature[22, 14, 23, 24]. For
the details of the broad feature fitting see Methods.

Figure 1: Temporal analysis of the GECAM multi-band light curves of GRB
230307A.
Left: The black histograms are binned net count rates (with background
subtracted); The blue solid curves are best-fit FRED model with the Norris05
formulation, while the blue dashed curves are the best-fit with our new FRED
formulation. Gray shadowed regions (the precursor and the dip) are ignored in
the fitting. All error bars represent 1$\sigma$ uncertainties of the net count
rates.
Upper right: The residuals of the net light curves after fitting with the new
FRED formulation.
Lower right: The cross-correlation between the multi-band residuals and that
in the 6-30 keV channel. The vertical dotted line indicates the zero lag time.
Figure 2: The fitted FRED formulation parameters as a function of energy. The
vertical error bars indicate the 1$\sigma$ uncertainties of the fitted
parameters, while the horizontal error bars indicate the ranges of the energy
bins.

In order to further demonstrate the self-similar feature of the slow varying
profile, we smooth the light curves in a non-parametrical way, and manifest
their shape identity after energy-dependent re-scaling in time (see Methods).
We stack all the re-scaled and smoothed light curves in the left panel of
extended data Figure 1. It can be obviously shown that, the shape of the
profiles in different energy bands can be "stretched" in time domain into an
identical shape with an energy dependent scaling factor. The scaling factor as
a function of energy is plotted in the right panel of extended data Figure 1.

The self-similarity inspired us to propose another formulation of FRED
profile, with one less parameter than that of Norris05:

$L(t)\propto\frac{t-t_{\rm{s}}}{\tau}\exp(-\frac{t-t_{\rm{s}}}{\tau}).$ (2)

We fit the multi-band light curves again with the new FRED formulation. The
results of fitting are listed in Methods. The sole time scale parameter $\tau$
shows a similar energy dependence on energy, $\tau\propto
E^{-0.35^{+0.05}_{-0.05}}$, from 6 to 300 keV, and shows a hint of shallower
energy dependence above $\sim 300$ keV (Figure 2). In the new formulation, one
can define peak time as $t_{\rm p}=t_{s}+\tau_{E}$, where we denote the $\tau$
in the $E$ channel with $\tau_{E}$. The width of the FRED profile also scales
as $\tau_{E}$. Therefore, the found $\tau-E$ dependence can naturally result
in the $w-E$ and $\tilde{t}_{\rm{p}}-E$ relations.

Beside the self-similar FRED profile of the broad feature of GRB0307A’s light
curve, the extraordinary brightness of the burst allows us to study the fast
varying temporal structures. From the residuals of the light curves (upper
right panel of Figure 1), one can see that there are many rapidly-varying
short pulses. One noticeable feature is that even the broad pulse has a clear
energy dependence, the short time spikes and dips appear to align at the same
instances across different energy bands. We further demonstrate this alignment
of the fast temporal features by cross-correlations among multi-band light
curve residuals. As shown in the lower right panel of Figure 1, the peaks are
perfectly aligned across the full energy band. This immediately excludes the
attempt to explain the self-similar stretchable feature of the multi-band
light curves as relativistic time-dilation, as a time-dilation effect would
stretch both slow and fast varying pulses together.

We show that some of these prominent short duration structures can be fitted
with individual short pulses (see Methods). The negative residuals, most
prominently the dip at around $\sim$18 s in all energy bands, essentially
excludes the possibility that short time structures are added onto a broad
pulse component. It is also difficult to explain this dip as absorption, since
the effective optical depth shows an energy dependence which no known
absorption mechanism can reproduce (see Methods). The most likely explanation
is that, the broad pulse of the light curves is composed of many fast pulses,
and the dips are gaps between successive short pulses (see Methods).

The aforementioned results serve as important tests on the GRB prompt emission
models. The internal shock models interpret short-time variabilities as a
result of collisions of pairs of shells. As a result, the typical emission
region is at[2, 4] $R_{\rm IS}\sim\Gamma^{2}c\delta t\simeq(3\times 10^{14}\
{\rm cm})(\Gamma/100)^{2}(\delta t/1\ {\rm s})$ from the central engine, where
$\Gamma\sim 100$ is the bulk Lorentz factor of the GRB, and $\delta t\sim 1$ s
is the typical variability time scale of the fast component. Within such a
scenario, the broad pulse should be defined as the history of the central
engine activity, which should be the same for all energy bands. The clear
energy-dependence of the broad pulse abandons this interpretation and has to
appeal to another internal shock at a much larger emission radius defined by
$R_{2}\sim\Gamma^{2}c\Delta t\sim(1.2\times 10^{16}\ {\rm
cm})(\Gamma/100)^{2}(\Delta t/40\ {\rm s})$ to interpret the slow component,
where $\Delta t\sim 40$ s is the timescale of the broad pulse, which is
essentially the GRB duration itself. Within such a picture, the observed
emission should be the superposition between the fast and slow components, but
the very deep, achromatic dip at around 18 s has essentially ruled out this
possibility. We therefore conclude that the internal shock model is clearly
disfavored.

The ICMART model, on the other hand, offers a natural interpretation to the
data. Within this model[3, 5, 17], only the slow component with a
characteristic timescale $\Delta t$ is related to the central engine activity
duration and hence, the typical thickness of the magnetic blobs. The collision
site between two magnetically-dominated shells is at $R_{\rm
ICMART}=R_{2}\sim(1.2\times 10^{16}\ {\rm cm})(\Gamma/100)^{2}(\Delta t/40\
{\rm s})$. The rapid-varying fast component with timescale $\delta t$, on the
other hand, is related to the emission timescale of local mini-jets within the
global emission zone. The broad pulse is simply the superposition of emission
of many mini-jets. No underlying slow component is needed. Significant dips
are allowed, especially during the decay phase when high-latitude emission
starts to play a role. The location of the dip is consistent with such a
phase[11]. Because the fast pulses originate from local mini-jet events, the
broad-band emission of the fast component is related to the dynamics of the
mini-jets and therefore perfectly aligned in time.

The general trend of energy-dependent peak time and duration is naturally
expected within the framework of the ICMART model. Unlike the internal shock
model that attributes different short pulses as emission from different
emitting fluids, the ICMART model suggests that the emission of the entire
broad pulse originates from the same fluid as it expands in space. The
magnetic field strength in the emission region decays with radius due to
expansion[25], leading to rolling down of the characteristic synchrotron
emission frequency with time. As result, higher energy emission peaks earlier
and lasts shorter than lower energy emission[7, 26]. The general trend of
model prediction matches the observations well[7], even though the self-
similarity revealed from the data may require special model parameters, which
may be related to the details of magnetic turbulence and reconnection in the
emission region, and thus provide new insights into the physics of gamma-ray
bursts.

##  1 Method

###  1.1 Observations

GECAM is a dedicated all-sky gamma-ray monitor constellation funded by the
Chinese Academy of Sciences, and now consists of three telescopes, i.e.
GECAM-A and GECAM-B[27] micro-satellties launched together on December 10,
2020, and GECAM-C (also called High Energy Burst Searcher, HEBS)[28] onboard
SATech-01 experimental satellite launched on July 27, 2022. There are two
kinds of detectors in each GECAM telescope: Gamma-Ray Detectors (GRDs) and
Charged Particle Detectors (CPDs). GRDs are the main detector of GECAM, each
of which is composed of a scintillator and an array of SiPMs. There are 25
GRDs onboard each of GECAM-A and GECAM-B, and 12 GRDs onboard GECAM-C.

GRB 230307A triggered GECAM-B in real-time at 15:44:06.650 UT on 7 March 2023
(denoted as T0) and also detected by GECAM-C and other gamma-ray monitors
(e.g. Fermi/GBM, Konus-Wind), while GECAM-A was offline at that time. The
real-time alert data was transmitted instantly with the Global Short Message
Communication of Beidou satellite navigation system[29, 30], and processed by
the automatic pipeline of GECAM[31], based on which the extreme brightness of
GRB 230307A was firstly reported by GECAM to the community[9], initiating many
multi-wavelength follow-up observations[32, 33].

Thanks to the dedicated design of instrument[28, 34], neither GECAM-B nor
GECAM-C suffered from data saturation during the whole burst of GRB 230307A
despite of its extreme brightness[11]. High quality of GECAM data allow us to
accurately measure the temporal and spectral properties of GRB 230307A. GRD04
of GECAM-B and GRD01 of GECAM-C are selected for the analysis of light curves
because of their smallest incident angle to the direction of GRB 230307A.
These two detectors both operate in two readout channels: high gain (HG) and
low gain (LG), which are independent in terms of data processing,
transmission, and dead-time.

At the time of this burst, GECAM-C GRDs have a lower energy detection
threshold of about 6 keV (owing to less radiation damage on SiPM) while
GECAM-B GRDs have a relatively higher energy detection threshold of about 30
keV. For GRD04 of GECAM-B, the energy range of HG channel data are used from
about 30 keV to 300 keV while the energy range of LG channel data are used
from about 300 keV to 1000 keV. For GRD01 of GECAM-C, only HG channel data are
used with the energy range from 6 keV to 30 keV. Though the response of GRD01
of GECAM-C for 6-15 keV and GRD04 of GECAM-B for 300-700 keV is affected by
the electronics, this does not have any effect on the analysis of light
curves.

The background of GECAM-B is estimated by fitting the data from $T_{0}$-50 s
to $T_{0}$-5 s and $T_{0}$+160 s to $T_{0}$+200 s with the first order
polynomials. The background of GECAM-C is estimated by fitting the data from
from $T_{0}$-20 s to $T_{0}$-1 s and $T_{0}$+170 s to $T_{0}$+600 s with a
combination of the first and second-order exponential polynomials[11].

###  1.2 Time-binned light curve and broad pulse fitting

We divide the data into eight energy bands (1st to 8th channels), namely 6-30
keV, where the data are from GECAM-C; 30-70 keV, 70-100 keV, 100-150 keV,
150-200 keV, 200-300 keV, 300-500 keV and 500-1000 keV, where the data are
from GECAM-B. The counts are binned every one second, and the background is
subtracted to calculate the net count rate. The multiband net count rates are
then fitted to both FRED formulation with Norris05 and our new FRED function
(equations 1,2). The fittings are done with the maximum likelihood method,
which can be expressed as ${\rm ln}\mathcal{L}=\sum_{i}\frac{({\rm
DATA}_{i}-{\rm MODEL}_{i})^{2}}{\sigma_{i}^{2}}$ by assuming a Gaussian
distribution of the observed net count rate, where the subscript i runs over
all time bins from -4 s to 100 s (excluding 0 s to 1 s and 16 s to 20 s). The
posteriors of the fitted parameters are found with a Monte Carlo Markov Chain
(MCMC) method. The fitted parameters are listed in extended data Table 1.

In the upper-right panel of Figure 1, we plot the the FRED parameters, namely
$w$ and $\tilde{t}_{\rm{p}}$ of Norris05 profile and $\tau$ of the new FRED
profile, as functions of energy. The vertical error bars are the uncertainties
of the fitted parameters, and the horizontal error bars indicate the width of
each band. We fit power law relations between the FRED parameters and their
corresponding energies. The power law fitting include the points from the 1st
to the 6th channels (6-300 keV). The power index are
-0.36${}^{+0.05}_{-0.06}$, -0.37${}^{+0.05}_{-0.06}$ and
-0.35${}^{+0.05}_{-0.05}$ for $w-E$, $\tilde{t}_{\rm{p}}-E$ and $\tau-E$
relations respectively. The uncertainties of the power index are found with
Monte-Carlo samplings of the FRED parameters and $E$ from their distributions.
The distributions of the FRED parameters are assumed to be Gaussian with their
fitted uncertainties, while the distribution of the energy is the
corresponding energy deposition spectrum in each band.

Table 1: Fitting results for the net light curves of GRB 230307A. All errors represent the 1$\sigma$ uncertainties. | Norris05 | New FRED
---|---|---
Energy range | $\tau_{\rm r}$ | $\tau_{\rm d}$ | norm |  ts | $\tau$ | norm |  ts
(keV) | (s) | (s) |  (counts$\cdot$s-1) | (s) | (s) |  (counts$\cdot$s-1) | (s)
6-30 | $11.11^{+0.26}_{-0.26}$ | $19.76^{+0.11}_{-0.10}$ | $2.49^{+0.05}_{-0.05}\times 10^{4}$ | $-3.42^{+0.09}_{-0.09}$ | $13.74^{+0.03}_{-0.03}$ | $1.47^{+0.004}_{-0.004}\times 10^{4}$ | $-9.78^{+0.28}_{-0.29}\times 10^{-1}$
30-70 | $6.84^{+0.13}_{-0.13}$ | $13.61^{+0.06}_{-0.07}$ | $2.08^{+0.04}_{-0.03}\times 10^{4}$ | $-1.79^{+0.04}_{-0.04}$ | $9.43^{+0.02}_{-0.02}$ | $1.33^{+0.004}_{-0.004}\times 10^{4}$ | $-5.16^{+0.08}_{-0.08}\times 10^{-1}$
70-100 | $5.46^{+0.13}_{-0.12}$ | $10.73^{+0.06}_{-0.06}$ | $1.28^{+0.03}_{-0.03}\times 10^{4}$ | $-1.39^{+0.04}_{-0.04}$ | $7.56^{+0.02}_{-0.02}$ | $8.05^{+0.03}_{-0.03}\times 10^{3}$ | $-5.06^{+0.06}_{-0.05}\times 10^{-1}$
100-150 | $4.55^{+0.10}_{-0.10}$ | $9.52^{+0.05}_{-0.05}$ | $1.52^{+0.03}_{-0.03}\times 10^{4}$ | $-1.23^{+0.03}_{-0.03}$ | $6.67^{+0.02}_{-0.02}$ | $9.92^{+0.03}_{-0.04}\times 10^{3}$ | $-5.08^{+0.04}_{-0.04}\times 10^{-1}$
150-200 | $4.51^{+0.13}_{-0.12}$ | $8.34^{+0.06}_{-0.06}$ | $1.10^{+0.03}_{-0.03}\times 10^{4}$ | $-1.25^{+0.04}_{-0.03}$ | $5.98^{+0.02}_{-0.02}$ | $6.60^{+0.03}_{-0.03}\times 10^{3}$ | $-5.15^{+0.04}_{-0.05}\times 10^{-1}$
200-300 | $3.99^{+0.13}_{-0.12}$ | $7.68^{+0.06}_{-0.06}$ | $9.96^{+0.28}_{-0.26}\times 10^{3}$ | $-1.07^{+0.04}_{-0.04}$ | $5.48^{+0.02}_{-0.02}$ | $6.12^{+0.03}_{-0.03}\times 10^{3}$ | $-4.98^{+0.07}_{-0.03}\times 10^{-1}$
300-500 | $4.37^{+0.14}_{-0.15}$ | $7.08^{+0.06}_{-0.06}$ | $1.17^{+0.04}_{-0.04}\times 10^{4}$ | $-1.23^{+0.04}_{-0.04}$ | $5.19^{+0.02}_{-0.02}$ | $6.35^{+0.03}_{-0.03}\times 10^{3}$ | $-5.12^{+0.05}_{-0.05}\times 10^{-1}$
500-1000 | $5.23^{+0.29}_{-0.28}$ | $6.08^{+0.09}_{-0.09}$ | $6.66^{+0.44}_{-0.40}\times 10^{3}$ | $-1.38^{+0.07}_{-0.07}$ | $4.73^{+0.03}_{-0.03}$ | $2.74^{+0.02}_{-0.02}\times 10^{3}$ | $-4.93^{+0.29}_{-0.15}\times 10^{-1}$

###  1.3 Profile self-similarity checking

In the left panel of extended data Figure 1, we stack the time re-scaled
smoothed light curves, to show that the broad features of the light curves in
different energy bands share the identical shapes. The smoothed light curves
are convolution between a Gaussian kernel and the original light curves. The
sigma of the Gaussian kernel is 1.5 s. We use the smoothed light curve in the
500-1000 keV band as the template, and we fit the template to light curves in
other bands, which are scaled by multiplying scaling factors $\alpha_{\rm{E}}$
to the their time argument $T-T_{0}$. It is intuitively shown that the self-
similarity of the broad profile, and this conclusion is FRED formulation
independent. The scaling factor as a function of energy is plotted in the
right panel of extended data Figure 1, which again shows a power law
dependence.

Figure 1: Left: The smoothed and time-rescaled GECAM multiband light curves of
GRB 230307A.
Right: The scaling factor of the Gaussian smoothed light curve as a function
of energy. The vertical error bars indicate the 1$\sigma$ uncertainties of the
fitted parameters, while the horizontal error bars indicate the ranges of the
energy bins.

###  1.4 Fast pulses

After the FRED profiles are removed from the multi-band light curves, the
residuals show fast pulses (Upper right panel of Figure 1). It is shown
clearly that the fast temporal structures are aligned in time across all
energy ranges. Here, we further test this temporal alignment with cross-
correlation between the residuals. In the lower right panel of Figure 1, we
plot the normalized correlation coefficient as a function of lag time between
the residuals in a certain band and that in band 6-30 keV. As we can see, all
the cross-correlation peak at lag time zero, indicating the alignment of the
fast pulses.

Here we further demonstrate that some of the prominent fast structures can be
identified as small pulses. For instance, we show that the spike at $\sim 3.5$
s can be fitted with a pulse profile. We employ our new FRED formulation to
fit the pulse from 3 s to 3.75 s, and find the typical width-$E$ relation
($\tau-E$ in the new FRED formulation) in GRB pulses. In order to show this
conclusion is independent on the pulse profile formulation, we use Gaussian
function as an alternative profile to fit the spike, and obtain the same
width-$E$ relation ($\sigma-E$ with the Gaussian pulse formulation. See the
lower right panel of extended data Figure 2). It can be seen in the left panel
of extended data Figure 2 that the FRED profile results in a better fit than
the Gaussian profile. The fitted parameters are list in extended data Table 2.
It is also interesting to observe that this small pulse does not show spectrum
lag of its peak time as the broad pulse. We demonstrate this by cross-
correlating between the small pulses in multi-bands and that in 30-70 keV
channel. The lag time together with their uncertainties as a function of
energy is shown in the upper right panel of extended data Figure 2, where all
lag times are almost zero, and we can clearly conclude that the lag time does
not depend on energy like the broad pulse.

Table 2: Light curves fitting results for the fast varying temporal structure of light curves at about 3.5 s . All errors represent the 1$\sigma$ uncertainties. | Gaussian | New FRED
---|---|---
Energy range | $\mu$ | $\sigma$ | norm | $\tau$ | norm |  ts
(keV) | (s) | (s) |  (counts$\cdot$s-1) | (s) |  (counts$\cdot$s-1) | (s)
30-70 | $3.36^{+0.01}_{-0.01}$ | $3.20^{+0.11}_{-0.09}\times 10^{-1}$ | $4.78^{+0.09}_{-0.09}\times 10^{3}$ | $2.88^{+0.09}_{-0.08}\times 10^{-1}$ | $1.37^{+0.03}_{-0.03}\times 10^{4}$ | $2.99^{+0.01}_{-0.01}$
70-100 | $3.38^{+0.01}_{-0.01}$ | $2.54^{+0.07}_{-0.07}\times 10^{-1}$ | $3.64^{+0.08}_{-0.08}\times 10^{3}$ | $2.54^{+0.07}_{-0.06}\times 10^{-1}$ | $1.02^{+0.02}_{-0.02}\times 10^{4}$ | $3.02^{+0.003}_{-0.003}$
100-150 | $3.39^{+0.004}_{-0.004}$ | $2.49^{+0.06}_{-0.05}\times 10^{-1}$ | $5.09^{+0.09}_{-0.09}\times 10^{3}$ | $2.28^{+0.07}_{-0.07}\times 10^{-1}$ | $1.47^{+0.03}_{-0.03}\times 10^{4}$ | $3.07^{+0.01}_{-0.01}$
150-200 | $3.40^{+0.004}_{-0.004}$ | $2.03^{+0.04}_{-0.03}\times 10^{-1}$ | $4.08^{+0.09}_{-0.09}\times 10^{3}$ | $1.93^{+0.05}_{-0.05}\times 10^{-1}$ | $1.15^{+0.03}_{-0.03}\times 10^{4}$ | $3.11^{+0.003}_{-0.004}$
200-300 | $3.38^{+0.004}_{-0.004}$ | $1.90^{+0.04}_{-0.03}\times 10^{-1}$ | $4.24^{+0.09}_{-0.10}\times 10^{3}$ | $1.76^{+0.04}_{-0.04}\times 10^{-1}$ | $1.23^{+0.03}_{-0.03}\times 10^{4}$ | $3.11^{+0.003}_{-0.004}$
300-500 | $3.39^{+0.003}_{-0.004}$ | $1.87^{+0.03}_{-0.03}\times 10^{-1}$ | $4.46^{+0.09}_{-0.09}\times 10^{3}$ | $1.73^{+0.04}_{-0.04}\times 10^{-1}$ | $1.30^{+0.03}_{-0.03}\times 10^{4}$ | $3.12^{+0.003}_{-0.003}$
500-1000 | $3.37^{+0.004}_{-0.004}$ | $1.46^{+0.03}_{-0.03}\times 10^{-1}$ | $2.42^{+0.08}_{-0.08}\times 10^{3}$ | $1.42^{+0.04}_{-0.04}\times 10^{-1}$ | $6.72^{+0.22}_{-0.21}\times 10^{3}$ | $3.13^{+0.003}_{-0.003}$

Figure 2: Left: Pulse profile fitting to the fast varying temporal structure
at $\sim 3.5$ s. Red and blue dashed lines correspond to new FRED formulation
and Gaussian profile respectively. The shadowed region indicates the range of
data involved in the fitting. Upper right: Cross-correlation between the light
curves of small pulse in multi-bands as functions of the energy. Lower right:
The fitted width ($\tau$ for new FRED and $\sigma$ for Gaussian profile) of
pulse at $\sim 3.5$ s as function of energy.

The dip at $\sim 18$ s is another noticeable feature of the multi-band light
curves. A natural attempt is to attribute such a dip to some sort of
absorption or geometrical blocking. We define the effective “optical depth"
$\tau_{\rm{op},E}$ in different energy bands as:

$\exp(-\tau_{\rm{op},E})=\frac{C_{\rm{dip},E}}{C_{\rm{slow},E}},$ (3)

where $C_{\rm dip}$ is the net count rate at the bottom of the dip, which is
found with a negative Gaussian fitting superposed on the FRED of broad pulse
fitting from 17 s to 19.5 s (upper left panel of Extended Data Figure 3); and
$C_{\rm slow}$ is the net count rate of the broad feature. The results of the
fitting with a negative Gaussian is tabulate in extended data Table 3. It is
intriguing to find that there is a power law energy dependence of the
effective optical depth, and the best fit power law index is $\sim 0.4$ (the
lower left panel of extended data Figure 3). This energy dependence of the
optical depth challenges the absorption picture, for there is no known
absorption mechanism whose cross section proportional to energy to the order
of $\sim 0.4$.

Figure 3: Upper left: fitting the gap with a FRED profile and a negative
Gaussian function. The shadowed region indicates the range of data involved in
the fitting. Lower left: the effective optical depth defined with equation 3,
as a function of energy. Upper right: fitting the gap with two successive
short pulses. The dashed red curves indicate the earlier pulse component, and
the dashed blue curves indicate the rising edge of the later pulse. Cyan
curves indicating the summation of these two components. Lower right: the time
scale factor $\tau$ for new FRED of the earlier pulse component in the 18-s
dip, as a function of energy.

Here we view the dip as gaps between two successive fast temporal structures.
Therefore, we fit a pulse profile (new FRED formulation) and a rising edge of
another pulse (Norris05) to the dip at each energy bands from $\sim$17 s to
19.5 s (the upper right panel of extended data Figure 3). The fitted
parameters are list in extended data Table 4. The time scale $\tau$ of the
earlier pulses as a function of energy is plotted in the lower right panel of
extended data Figure 3, which shows clearly a typical power law energy
dependence of pulse width.

Table 3: Light curves fitting results for the dip and optical depth. All errors represent the 1$\sigma$ uncertainties. Energy range | $\mu$ | $\sigma$ | $C_{{\rm slow},E}$ | Optical Depth
---|---|---|---|---
(keV) | (s) | (s) |  (counts$\cdot$s-1) |
6-30 | $18.27^{+0.01}_{-0.01}$ | $4.68^{+0.10}_{-0.09}\times 10^{-1}$ | $5.08^{+0.01}_{-0.01}\times 10^{3}$ | $0.94^{+0.03}_{-0.02}$
30-70 | $18.25^{+0.01}_{-0.01}$ | $5.51^{+0.10}_{-0.10}\times 10^{-1}$ | $3.62^{+0.01}_{-0.01}\times 10^{3}$ | $1.55^{+0.04}_{-0.04}$
70-100 | $18.21^{+0.01}_{-0.01}$ | $6.13^{+0.14}_{-0.13}\times 10^{-1}$ | $1.68^{+0.01}_{-0.01}\times 10^{3}$ | $2.07^{+0.09}_{-0.08}$
100-150 | $18.25^{+0.01}_{-0.01}$ | $6.99^{+0.15}_{-0.15}\times 10^{-1}$ | $1.68^{+0.01}_{-0.01}\times 10^{3}$ | $2.25^{+0.10}_{-0.09}$
150-200 | $18.19^{+0.01}_{-0.01}$ | $6.70^{+0.18}_{-0.18}\times 10^{-1}$ | $9.06^{+0.05}_{-0.05}\times 10^{2}$ | $2.55^{+0.18}_{-0.15}$
200-300 | $18.19^{+0.01}_{-0.02}$ | $6.77^{+0.21}_{-0.20}\times 10^{-1}$ | $6.90^{+0.05}_{-0.05}\times 10^{2}$ | $3.35^{+0.50}_{-0.33}$
300-500 | $18.18^{+0.02}_{-0.02}$ | $7.26^{+0.24}_{-0.23}\times 10^{-1}$ | $6.23^{+0.05}_{-0.05}\times 10^{2}$ | $2.67^{+0.29}_{-0.20}$
500-1000 | $18.12^{+0.04}_{-0.04}$ | $8.08^{+0.54}_{-0.47}\times 10^{-1}$ | $2.11^{+0.03}_{-0.03}\times 10^{2}$ | $3.65^{+1.10}_{-0.64}$
Table 4: Light curves fitting results for the gap with two successive fast temporal structures. All errors represent the 1$\sigma$ uncertainties. | pluse before gap | pluse after gap | additional constant
---|---|---|---
Energy range |  ts | $\tau$ | norm |  ts | $\tau_{\rm r}$ | norm | constant
(keV) | (s) | (s) |  (counts$\cdot$s-1) | (s) | (s) |  (counts$\cdot$s-1) |  (counts$\cdot$s-1)
6-30 | $16.92^{+0.06}_{-0.10}$ | $3.19^{+0.63}_{-0.37}\times 10^{-1}$ | $9.22^{+0.87}_{-0.49}\times 10^{3}$ | $18.13^{+0.07}_{-0.08}$ | $1.29^{+0.25}_{-0.18}$ | $1.13^{+0.19}_{-0.13}\times 10^{4}$ | $1.54^{+0.20}_{-0.35}\times 10^{3}$
30-70 | $17.00^{+0.03}_{-0.04}$ | $2.74^{+0.30}_{-0.22}\times 10^{-1}$ | $7.24^{+0.32}_{-0.28}\times 10^{3}$ | $18.23^{+0.04}_{-0.07}$ | $1.10^{+0.19}_{-0.08}$ | $7.81^{+1.02}_{-0.47}\times 10^{3}$ | $5.20^{+1.05}_{-1.48}\times 10^{2}$
70-100 | $17.03^{+0.02}_{-0.03}$ | $2.48^{+0.24}_{-0.21}\times 10^{-1}$ | $3.42^{+0.15}_{-0.14}\times 10^{3}$ | $18.21^{+0.05}_{-0.10}$ | $1.16^{+0.31}_{-0.12}$ | $3.76^{+0.74}_{-0.32}\times 10^{3}$ | $1.11^{+0.46}_{-0.55}\times 10^{2}$
100-150 | $17.15^{+0.02}_{-0.02}$ | $2.21^{+0.17}_{-0.14}\times 10^{-1}$ | $3.15^{+0.14}_{-0.13}\times 10^{3}$ | $18.23^{+0.05}_{-0.07}$ | $1.23^{+0.27}_{-0.16}$ | $3.44^{+0.67}_{-0.37}\times 10^{3}$ | $86.96^{+14.91}_{-27.78}$
150-200 | $17.07^{+0.01}_{-0.02}$ | $2.16^{+0.13}_{-0.13}\times 10^{-1}$ | $1.97^{+0.10}_{-0.10}\times 10^{3}$ | $17.89^{+0.18}_{-0.26}$ | $2.25^{+1.28}_{-0.71}$ | $2.98^{+2.08}_{-0.90}\times 10^{3}$ | $17.23^{+19.08}_{-12.19}$
200-300 | $17.04^{+0.02}_{-0.03}$ | $2.03^{+0.12}_{-0.10}\times 10^{-1}$ | $1.63^{+0.08}_{-0.09}\times 10^{3}$ | $17.42^{+0.24}_{-0.21}$ | $5.09^{+1.71}_{-1.55}$ | $6.52^{+5.00}_{-2.79}\times 10^{3}$ | $3.32^{+3.69}_{-2.44}$
300-500 | $17.05^{+0.02}_{-0.04}$ | $1.95^{+0.16}_{-0.13}\times 10^{-1}$ | $1.59^{+0.10}_{-0.09}\times 10^{3}$ | $17.30^{+0.25}_{-0.15}$ | $5.86^{+1.23}_{-1.67}$ | $6.76^{+3.92}_{-2.92}\times 10^{3}$ | $11.75^{+7.42}_{-7.67}$
500-1000 | $17.01^{+0.06}_{-0.18}$ | $1.99^{+0.34}_{-0.24}\times 10^{-1}$ | $4.43^{+0.76}_{-0.55}\times 10^{2}$ | $17.66^{+0.32}_{-0.38}$ | $3.69^{+2.80}_{-1.90}$ | $9.28^{+16.28}_{-5.48}\times 10^{2}$ | $1.91^{+2.61}_{-1.42}$

## Data Availability

The processed data are presented in the tables and figures of the paper, which
are available upon reasonable request.

## Code Availability

Upon reasonable requests, the code (mostly in Python) used to produce the
results and figures will be provided.

## References

  * [1] Zhang, B.  The Physics of Gamma-Ray Bursts  (2018).
  * [2] Rees, M. J. & Meszaros, P.  Unsteady Outflow Models for Cosmological Gamma-Ray Bursts.  ApJ 430, L93 (1994).  astro-ph/9404038.
  * [3] Zhang, B. & Yan, H.  The Internal-collision-induced Magnetic Reconnection and Turbulence (ICMART) Model of Gamma-ray Bursts.  ApJ 726, 90 (2011).  1011.1197.
  * [4] Kobayashi, S., Piran, T. & Sari, R.  Can Internal Shocks Produce the Variability in Gamma-Ray Bursts?  ApJ 490, 92 (1997).  astro-ph/9705013.
  * [5] Zhang, B. & Zhang, B.  Gamma-Ray Burst Prompt Emission Light Curves and Power Density Spectra in the ICMART Model.  ApJ 782, 92 (2014).  1312.7701.
  * [6] Daigne, F. & Mochkovitch, R.  Gamma-ray bursts from internal shocks in a relativistic wind: temporal and spectral properties.  MNRAS 296, 275–286 (1998).  astro-ph/9801245.
  * [7] Uhm, Z. L. & Zhang, B.  Toward an Understanding of GRB Prompt Emission Mechanism. I. The Origin of Spectral Lags.  ApJ 825, 97 (2016).  1511.08807.
  * [8] Fermi GBM Team.  GRB 230307A: Fermi GBM Final Real-time Localization.  GRB Coordinates Network  33405, 1 (2023).
  * [9] Xiong, S., Wang, C., Huang, Y. & Gecam Team.  GRB 230307A: GECAM detection of an extremely bright burst.  GRB Coordinates Network  33406, 1 (2023).
  * [10] Levan, A.  et al.  JWST detection of heavy neutron capture elements in a compact object merger.  arXiv e-prints  arXiv:2307.02098 (2023).  2307.02098.
  * [11] Sun, H.  et al.  Magnetar emergence in a peculiar gamma-ray burst from a compact star merger.  arXiv e-prints  arXiv:2307.05689 (2023).  2307.05689.
  * [12] Fenimore, E. E., in ’t Zand, J. J. M., Norris, J. P., Bonnell, J. T. & Nemiroff, R. J.  Gamma-Ray Burst Peak Duration as a Function of Energy.  ApJ 448, L101 (1995).  astro-ph/9504075.
  * [13] Shen, R.-F., Song, L.-M. & Li, Z.  Spectral lags and the energy dependence of pulse width in gamma-ray bursts: contributions from the relativistic curvature effect.  MNRAS 362, 59–65 (2005).  astro-ph/0505276.
  * [14] Liang, E.-W.  et al.  Temporal Profiles and Spectral Lags of XRF 060218.  ApJ 653, L81–L84 (2006).  astro-ph/0610956.
  * [15] Kazanas, D., Titarchuk, L. G. & Hua, X.-M.  The Duration–Photon Energy Relation of Gamma-Ray Bursts and Its Interpretations.  ApJ 493, 708–714 (1998).  astro-ph/9709180.
  * [16] Qin, Y. P., Dong, Y. M., Lu, R. J., Zhang, B. B. & Jia, L. W.  Relationship between the Gamma-Ray Burst Pulse Width and Energy Due to the Doppler Effect of Fireballs.  ApJ 632, 1008–1020 (2005).  astro-ph/0411365.
  * [17] Shao, X. & Gao, H.  Gamma-Ray Burst Prompt Emission Spectrum and E p Evolution Patterns in the ICMART Model.  ApJ 927, 173 (2022).  2201.03750.
  * [18] Lyutikov, M. & Blandford, R.  Gamma Ray Bursts as Electromagnetic Outflows.  arXiv e-prints  astro–ph/0312347 (2003).  astro-ph/0312347.
  * [19] Narayan, R. & Kumar, P.  A turbulent model of gamma-ray burst variability.  MNRAS 394, L117–L120 (2009).  0812.0018.
  * [20] Lazarian, A., Zhang, B. & Xu, S.  Gamma-Ray Bursts Induced by Turbulent Reconnection.  ApJ 882, 184 (2019).  1801.04061.
  * [21] Hascoët, R., Daigne, F. & Mochkovitch, R.  Accounting for the XRT early steep decay in models of the prompt gamma-ray burst emission.  A&A 542, L29 (2012).  1206.6770.
  * [22] Norris, J. P.  et al.  Long-Lag, Wide-Pulse Gamma-Ray Bursts.  ApJ 627, 324–345 (2005).  astro-ph/0503383.
  * [23] Hakkila, J., Giblin, T. W., Norris, J. P., Fragile, P. C. & Bonnell, J. T.  Correlations between Lag, Luminosity, and Duration in Gamma-Ray Burst Pulses.  ApJ 677, L81 (2008).  0803.1655.
  * [24] Peng, Z. Y., Zhao, X. H., Yin, Y., Bao, Y. Y. & Ma, L.  Energy-dependent Gamma-Ray Burst Pulse Width Due to the Curvature Effect and Intrinsic Band Spectrum.  ApJ 752, 132 (2012).  1204.3420.
  * [25] Uhm, Z. L. & Zhang, B.  Fast-cooling synchrotron radiation in a decaying magnetic field and $\gamma$-ray burst emission mechanism.  Nature Physics  10, 351–356 (2014).  1303.2704.
  * [26] Uhm, Z. L., Zhang, B. & Racusin, J.  Toward an Understanding of GRB Prompt Emission Mechanism. II. Patterns of Peak Energy Evolution and Their Connection to Spectral Lags.  ApJ 869, 100 (2018).  1801.09183.
  * [27] Li, X. Q.  et al.  The technology for detection of gamma-ray burst with GECAM satellite.  Radiation Detection Technology and Methods  6, 12–25 (2021).
  * [28] Zhang, D.  et al.  The performance of SiPM-based gamma-ray detector (GRD) of GECAM-C.  Nuclear Instruments and Methods in Physics Research A  1056, 168586 (2023).  2303.00537.
  * [29] Zhao, X.-Y.  et al.  The In-Flight Realtime Trigger and Localization Software of GECAM.  arXiv e-prints  arXiv:2112.05101 (2021).
  * [30] Guo, S.  et al.  Integrated Navigation and Communication Service for LEO Satellites Based on BDS-3 Global Short Message Communication.  IEEE Access  11, 6623–6631 (2023).
  * [31] Huang, Y.  et al.  The GECAM Real-Time Burst Alert System.  arXiv e-prints  arXiv:2307.04999 (2023).
  * [32] Bom, C. R.  et al.  GRB 230307A: SOAR/Goodman detection of the possible host galaxy.  GRB Coordinates Network  33459, 1 (2023).
  * [33] Levan, A. J.  et al.  GRB 230307A: JWST observations consistent with the presence of a kilonova.  GRB Coordinates Network  33569, 1 (2023).
  * [34] Liu, Y.  et al.  The sipm array data acquisition algorithm applied to the gecam satellite payload.  Radiation Detection Technology and Methods  6, 70–77 (2021).

S.-X.Y. thanks the insightful discussion with Prof. Kinwah Wu on the origin of
the self-similarity of the multi-wavelength lightcurves of the burst. S.-X.Y.
also discussed with Dr. X.L. Wang on the potential physical implication of
these interesting behaviour in data. This work is supported by the National
Key R&D Program of China (2021YFA0718500). S.-X.Y. acknowledges support from
the Chinese Academy of Sciences (grant Nos. E329A3M1 and E3545KU2). S.-L.X.
acknowledges the support by the National Natural Science Foundation of China
(Grant No. 12273042). The GECAM (Huairou-1) mission is supported by the
Strategic Priority Research Program on Space Science (Grant No. XDA15360000)
of Chinese Academy of Sciences. ChatGPT3.5 was used to check and correct a few
sentences regarding their English grammar and wording, without any changes on
their original meanings.

S.-L.X. and S.-X.Y. initiated the study. S.-X.Y. led the investigation and
wrote the manuscript. B.Z. and S.-X.Y. proposed the physical implications of
the data. S.-L.X. led the GECAM observation and data analysis. B.Z., C.-W.W.,
S.-L.X. contributed to the writing of the manuscript. C.-W.W. is the major
contributor of GECAM data analysis. W.-J.T. investigated the nature of the dip
$\sim 18$ s. S.-N.Z., R.M., S.-L.X. contributed to the theoretical
interpretation of self-similarity of the broad pulse. W.-C.X. and J.-C.L.
performed background analysis for GECAM-C. C.Z. and Y.-Q.Z. performed
calibration analysis for GECAM. J.-C.L. performed the data saturation
assessment for GECAM. P.Z. developed the data analysis tools for GECAM. Y.W.
contributed to preparing of the manuscript. All authors participated in the
discussion and proofreading of the manuscript, or contributed to the
development and operation of GECAM, which is important for this work.

The authors declare that they have no competing financial interests.

Correspondence and requests for materials should be addressed to S.-X.Y.
(sxyi@ihep.ac.cn), C.-W.W. (cwwang@ihep.ac.cn), B.Z.(bing.zhang@unlv.edu),
S.-L.X. (xiongsl@ihep.ac.cn).