In the main part of the paper we focused on the spin-spin connected
correlation function
$S^{z}_{m,n}(t)\equiv\braket{\sigma^{z}_{m}(t)\sigma^{z}_{n}(t)}-\braket{\sigma^{z}_{m}(t)}\braket{\sigma^{z}_{n}(t)}$,
with $m-n\sim t^{\alpha}$ and $\alpha\in[0,1]$. We studied three partitioning
protocols that were engineered in such a way to isolate the independent
contributions to $S^{z}_{m,n}$ and we argued that the contributions should be
summed independently for a generic partitioning protocol in which all the
dynamical fields have non-zero value over at least one of the homogeneous
states involved in the partitioning protocol. We saw that the leading order
for large time is obtained by considering the root density only, and
neglecting all the other fields, not only for $\alpha=0$, when the support of
the observable is fixed and finite (i.e. the GHD’s regime of validity), but
also for $\alpha<1/2$. We also gave an argument according to which the
crossover scaling $\alpha=1/2$ does not depend on the special observable we
are looking at. Indeed, when considering observables with larger (but finite)
support in the fermionic representation, the higher-order fields needed to
describe it give a contribution that decays faster to zero than the one coming
from the root density.

We point out that, for $\alpha=0$, the leading order from the fields different
from the root density is in general of the same order ($t^{-1}$) as the late-
time corrections from higher-order GHD, implying that, unless the homogeneous
states forming the partitioning protocol are stationary, a higher-order GHD
that relies on the root density only can not in general give the right
predictions. Note the important role of non-Gaussian correlations in the
initial state on the relaxation process. For instance, in global homogeneous
quenches from a Gaussian state, the GGE value of the spin-spin correlation
function in the Ising model is attained with corrections decaying as $t^{-3}$.
Non-Gaussian correlations slow this process down to $t^{-1}$.

The results above also show how Gaussianification takes place in partitioning
protocols. We argued that the expectation value of local observables at the
Euler scale can be described using a Gaussian state that depends on the
position of the observable. We also answered to the question how local is
“local enough” for the spin-spin connected correlation to Gaussianize. The
answer is that the distance between the spins should not grow faster than
$t^{1/2}$, since otherwise the non-Gaussian fields become dominant.

As for the generality of our arguments, the generalization to any quadratic
model is quite straightforward. As a matter of fact one needs to change the
dispersion relation, the Bogoliubov angle and, at most, introduce other angles
to diagonalize the Hamiltonian (see e.g. [41]), but the definition of the
dynamical fields would stay the same. Although we did not carry out the
computation in the general case, we argued that the qualitative behaviors of
the independent contributions to the spin-spin connected correlation would not
change, except at most, the one linked to the auxiliary field, which may decay
as $t^{-1}$ instead of $t^{-3+2\alpha}$, as it happens in the XY quantum spin
chain when the external magnetic field is such that the dispersion relation
has more than two stationary points. Note that this implies that a Gaussian
partitioning protocol in the Ising chain attains its locally quasi-stationary
state faster than in those other spin chains.

The generalization to other observables with 4-sites support in fermions is
also straightforward: the precise expression of the observable would change,
but, since our arguments rely only on the structure of time evolution and the
singularities that define the partitioning protocol, the conclusions would be
the same. We also discussed how observables that require a finite number of
higher order connected correlations can be included in the picture: $n+1$ new
independent dynamical fields should be introduced to account for the
$2n$-point fermionic connected correlation. We also gave a qualitative
argument on how to obtain their scaling.

It remains the question of what happens to local observables the are not local
in fermions, such as $\Braket{\sigma^{x}_{\ell}}$. We can not apply our
arguments to those observables because infinitely-many (in the thermodynamic
limit) dynamical fields are involved. The situation for those observables is
generally more complicated and it even turns out that there are situations in
which the usual GGE/GHD is not enough to describe the late-time behavior of
those observable, despite the model being free [72, 73]. We leave the
inclusion of those observables in our description as a prospect for a future
work.

## Acknowledgements

I am grateful to Maurizio Fagotti for giving me the idea behind this work. I
thank also Lenart Zadnik and Vanja Marić for useful discussions. This work was
supported by the European Research Council under the Starting Grant No. 805252
LoCoMacro.

##  Appendix A Complement to the root-density partitioning protocol

In the main text – Eqs. (65) and (66) – we identified two contributions to the
correlation matrix coming from the field $\rho_{x,t}(p)$ and we showed how to
compute one of them. In this section we use the same ideas to compute the
second contribution. Fig. 7 shows how one should deform the contours in this
case for the same rays used in Fig. 2. As can be inferred from the pictures,
the contour for the residue contribution can always be chosen to be outside
the unit circle. For both rays inside the lightcone, it comes in the form of
two curves that are symmetric with reflection with respect to the origin.
Calling $\gamma$ one of the two pieces, we have that the residue contribution
is

$\mathcal{R}^{\rho}_{2}[f(w_{1},w_{2})](z,t):=2\mathrm{i}\int_{\gamma}\frac{\mathrm{d}w}{2\pi}\mathrm{e}^{-2z\ln
w}f(w,w).$ (107)

For $\alpha=1$, the real part of $-2z\ln w$ is negative outside the unit
circle (recall $z>0$), so that such contribution is always exponentially
suppressed. As for the previous case, in this regime we can just deform the
two initial contours as if the pole was not there. In the end, the leading
contribution is given by a standard application of the saddle-point method to
the remaining double integral over $\mathcal{C}^{\prime}_{1}$ and
$\mathcal{C}^{\prime}_{2}$. When both rays are inside the lightcone, there are
four saddle points to be considered for each of those integrals, yielding a
leading order of $\mathcal{I}^{\rho}_{2}[f(w_{1},w_{2})](z,t)\sim t^{-1}$.
Whenever at least one of the two rays is outside the lightcone, at least one
integral of the double integral has no saddle points and the domain can be
deformed in such a way that the function under integration is exponentially
suppressed.

$x/t=0.3,z/t=3.0$$0$$0$$-1$$-1$$+1$$+1$

$x/t=-0.4,z/t=1.8$$0$$0$$-1$$-1$$+1$$+1$

$x/t=-1.3,z/t=0.2$$0$$0$$-1$$-1$$+1$$+1$

(a)

$x/t=0.8,z/t=0.6$$0$$0$$-1$$-1$$+1$$+1$

$x/t=1.3,z/t=0.2$$0$$0$$-1$$-1$$+1$$+1$

(b)

$x/t=-0.6,z/t=0.7$$0$$0$$-1$$-1$$+1$$+1$

$x/t=0.4,z/t=1.0$$0$$0$$-1$$-1$$+1$$+1$

$x/t=0.6,z/t=0.7$$0$$0$$-1$$-1$$+1$$+1$

(c)

Figure 7: Regions in which
$\operatorname{Re}(S_{(x-\frac{z}{2})/t}(a+\mathrm{i}b))>0$ (blue region) and
$\operatorname{Re}(S_{(x+\frac{z}{2})/t}(a+\mathrm{i}b))<0$ (orange region) in
function of $a$ (horizontal axis) and $b$ (vertical axis). The rest is the
same of Fig. 2.

In the case of $0<\alpha<1/2$, as we have already discussed, the residue
contribution is not negligible anymore, but gives instead what turns out to be
the leading contribution to the initial double-integral. As we have done for
the previous contribution, we linearize it around its extrema. Calling
$\mathrm{e}^{\mathrm{i}\bar{p}_{j}}+\delta
w_{j}\mathrm{e}^{\mathrm{i}\bar{p}_{j}}$, with $\delta w\sim t^{\alpha-1}$ and
$j\in\\{1,2\\}$, the extrema of the first piece of $\gamma$:

$\mathcal{R}^{\rho}_{2}[f(w_{1},w_{2})](z,t)\sim
2\mathrm{i}\sum_{j=1}^{2}(-1)^{j-1}\int_{0}^{\Delta}\frac{\mathrm{d}s}{2\pi}\mathrm{e}^{\mathrm{i}\bar{p}_{j}}\mathrm{e}^{-2z\ln(\mathrm{e}^{\mathrm{i}\bar{p}_{j}}+\delta
w_{j}\mathrm{e}^{\mathrm{i}\bar{p}_{j}}+s\mathrm{e}^{\mathrm{i}\bar{p}_{j}})}f(\mathrm{e}^{\mathrm{i}\bar{p}_{j}},\mathrm{e}^{\mathrm{i}\bar{p}_{j}}),$
(108)

where $\Delta\sim t^{-\omega}$, $\omega\in(\alpha/2,\alpha)$. Computing the
integral gives

$\mathcal{R}^{\rho}_{2}[f(w_{1},w_{2})](z,t)\sim-\mathrm{i}\sum_{j=1}^{2}(-1)^{j}\frac{1}{2\pi
z}\mathrm{e}^{\mathrm{i}\bar{p}_{j}}\mathrm{e}^{-\mathrm{i}2z\bar{p}_{j}}f(\mathrm{e}^{\mathrm{i}\bar{p}_{j}},\mathrm{e}^{\mathrm{i}\bar{p}_{j}}),$
(109)

where we have used $(1+\delta w+\Delta)^{-2z}\sim\mathrm{e}^{-2z\delta w}$, to
conclude that the term containing $\Delta$ goes exponentially to zero, and
$(1+\delta w)^{-2z}\sim 1$. This is the leading contribution to the integral.
Note that there is an implicit dependence on $x$, since $\bar{p}_{j}$ depends
on $x/t$.

For $1/2<\alpha<1$ instead, the residue contribution is exponentially
suppressed for the same reasons that we discussed in the main text and the
leading contribution comes from the remaining double integral. Using the same
notations introduced for the previous contribution, we get

$\mathcal{I}^{\rho}_{2}[f(w_{1},w_{2})](x,z,t)\sim\sum_{j=1}^{4}\int_{-\Delta}^{\Delta}\frac{\mathrm{d}s^{+}\mathrm{d}s^{-}}{(2\pi)^{2}}\mathrm{e}^{\mathrm{i}\phi^{+}_{j}}\mathrm{e}^{\mathrm{i}\phi^{-}_{j}}f(\mathrm{e}^{\mathrm{i}\bar{p}_{j}},\mathrm{e}^{\mathrm{i}\bar{p}_{j}})\times\\\
\times\frac{\mathrm{e}^{tS_{(x-\frac{z}{2})/t}(\mathrm{e}^{\mathrm{i}\bar{p}_{j}-\mathrm{i}\delta
p_{j}})-tS_{(x+\frac{z}{2})/t}(\mathrm{e}^{\mathrm{i}\bar{p}_{j}+\mathrm{i}\delta
p_{j}})-\frac{t}{2}|S^{\prime\prime}_{x/t}(\mathrm{e}^{\mathrm{i}\bar{p}_{j}+\mathrm{i}\delta
p_{j}})|(s^{-})^{2}-\frac{t}{2}|S^{\prime\prime}_{x/t}(\mathrm{e}^{\mathrm{i}\bar{p}_{j}-\mathrm{i}\delta
p_{j}})|(s^{+})^{2}}}{-2\mathrm{i}\mathrm{e}^{\mathrm{i}\bar{p}_{j}}\sin(\delta
p_{j})+s^{+}\mathrm{e}^{\mathrm{i}\phi^{+}_{j}}-s^{-}\mathrm{e}^{\mathrm{i}\phi^{-}_{j}}},$
(110)

where we took $\Delta\sim
t^{\omega},\omega\in(-1/2,\operatorname{min}\\{\alpha-1,-1/3\\})$ and the two
phases $\phi^{\pm}_{j}$ are approximately the same of the case discussed in
the main text. Performing the Gaussian integrals we get

$\mathcal{I}^{\rho}_{2}[f(w_{1},w_{2})](x,z,t)\sim\sum_{j=1}^{4}\frac{\mathrm{e}^{-tS_{(x+\frac{z}{2})/t}(\mathrm{e}^{\mathrm{i}\bar{p}_{j}+\mathrm{i}\delta
p_{j}})+tS_{(x-\frac{z}{2})/t}(\mathrm{e}^{\mathrm{i}\bar{p}_{j}-\mathrm{i}\delta
p_{j}})}}{-4\mathrm{i}\pi
t|S^{\prime\prime}_{x/t}(\mathrm{e}^{\mathrm{i}\bar{p}_{j}})|\mathrm{e}^{\mathrm{i}\bar{p}_{j}}\delta
p_{j}}\mathrm{e}^{\mathrm{i}\phi^{+}_{j}}\mathrm{e}^{\mathrm{i}\phi^{-}_{j}}f(\mathrm{e}^{\mathrm{i}\bar{p}_{j}},\mathrm{e}^{\mathrm{i}\bar{p}_{j}}).$
(111)

This is combined in the text we the other contribution
$\mathcal{I}_{1}^{\rho}$ to get the behavior of the correlation matrix
elements.

##  Appendix B Complement to the auxiliary-field partitioning protocol

In this section we compute the asymptotics of the integrals (96) in the
scaling limit $t\rightarrow+\infty$, with $x=\zeta
t,z=ct^{\alpha}>0,\alpha\in[0,1],\zeta\neq 0$. We can use the same arguments
that we used in the root-density case. The saddle-point landscape for the
first term of each integral is summarized in Fig. 8 and the situation is
similar for the second term (see Section B.4). As can be inferred from the
figure, it is always possible to deform $\mathcal{C}_{1}$ and
$\mathcal{C}_{2}$ to $\mathcal{C}^{\prime}_{1}$ and $\mathcal{C}^{\prime}_{2}$
such that the real part of the phase is negative almost everywhere. If
$x/t<0$, we do not get any residue contribution after the deformation. If
instead $x/t>0$, the two contours have to be fully exchanged and we get a
residue contribution in the form of an integral over a closed curve around the
origin. If at least one ray is outside the lightcone, meaning that
$|\zeta+\frac{z}{2t}|>v_{max}$ or $|\zeta-\frac{z}{2t}|>v_{max}$, then the
double integral that one gets after the deformation is exponentially
suppressed (the real part of the phase is negative everywhere) and we are left
only with the residue contribution, if present. When both the rays are inside
the lightcone, the function under the double-integral after the deformation is
exponentially suppressed except in the saddle points, where the real part of
the phase is zero. A standard application of the saddle point method shows
that the contribution of the double integral in this case is of order $\sim
t^{-1}$, independently from $\alpha$.

$x/t=-1.3,z/t=0.2$$0$$0$$-1$$-1$$+1$$+1$

$x/t=-0.4,z/t=1.8$$0$$0$$-1$$-1$$+1$$+1$

(a)

$x/t=1.3,z/t=0.2$$0$$0$$-1$$-1$$+1$$+1$

$x/t=1.3,z/t=2.2$$0$$0$$-1$$-1$$+1$$+1$

$x/t=0.8,z/t=0.8$$0$$0$$-1$$-1$$+1$$+1$

(b)

$x/t=-0.6,z/t=0.7$$0$$0$$-1$$-1$$+1$$+1$

$x/t=-0.4,z/t=1.0$$0$$0$$-1$$-1$$+1$$+1$

$x/t=0.4,z/t=1.0$$0$$0$$-1$$-1$$+1$$+1$

(c)

Figure 8: Regions in which
$\operatorname{Re}(S_{(x+\frac{z}{2})/t}(a+\mathrm{i}b))>0$ (blue region) and
$\operatorname{Re}(S_{(-x+\frac{z}{2})/t}(a+\mathrm{i}b))>0$ (orange region)
in function of $a$ (horizontal axis) and $b$ (vertical axis). We consider a
TFIC with $J=1/2$, so that $v_{max}=1$, and $h=2$. The black and red lines
represent respectively the curves $\mathcal{C}^{\prime}_{1}$ and
$\mathcal{C}^{\prime}_{2}$ obtained by deforming $\mathcal{C}_{1}$ and
$\mathcal{C}_{2}$ as detailed in the text. (a) Case in which the ray
$(x-\frac{z}{2})/t,z>0,$ is on the left of the lightcone. (b) Case in which
the ray $(x+\frac{z}{2})/t,z>0,$ is on the right of the lightcone and
$(x-\frac{z}{2})/t$ is either inside or to the right of it. (c) Case in which
the two rays $(x\pm\frac{z}{2})/t$ are both inside the lightcone.

The residue contribution, whenever it is present, comes in the form

$\mathcal{R}_{1}^{\psi}(z,t)\equiv-\frac{\mathrm{i}}{2}\operatorname{Re}\int_{\gamma}\frac{\mathrm{d}w}{2\pi
w}\psi^{R}(-2\mathrm{i}\ln(w))\left(\mathrm{e}^{2tS_{\frac{z}{2t}}(w)}-\mathrm{e}^{2tS_{-\frac{z}{2t}}(w)}\right)$
(112)

for $\left(\Gamma^{z;\psi}_{x}\right)_{11}$, and

$\mathcal{R}_{2}^{\psi}(z,t)\equiv\frac{\mathrm{i}}{2}\operatorname{Im}\int_{\gamma}\frac{\mathrm{d}w}{2\pi
w}\psi^{R}(-2\mathrm{i}\ln
w)\left(\mathrm{e}^{2tS_{\frac{z}{2t}}(w)}\mathrm{e}^{\mathrm{i}\theta(-2\mathrm{i}\ln
w)}-\mathrm{e}^{2tS_{-\frac{z}{2t}}(w)}\mathrm{e}^{-\mathrm{i}\theta(-2\mathrm{i}\ln
w)}\right)$ (113)

for $\left(\Gamma^{z;\psi}_{x}\right)_{21}$. Here $\gamma$ is a closed curve
in the annulus of analiticity around the unit circle. Note that for
$|z/t|>2v_{max}$ we can deform $\gamma$ (independently for each of the two
terms in a given integral) in such a way that the function under integration
is exponentially small; we conclude that in such a case the residue
contribution is exponentially small, consistently with the Lieb-Robinson
bound, since sites with $|z/t|>2v_{max}$ are causally disconnected. In all the
cases $|z/t|<2v_{max}$, instead, we have four saddle points on the unit
circle, in which the real part of the exponent is zero. We discuss this case
distinguishing between the different scaling regimes.

###  B.1 Linear scaling

Calling $\mathrm{e}^{\mathrm{i}\bar{p}_{j}^{\pm}},\forall j\in\\{1,2,3,4\\}$
the saddle points for the $S_{\pm\frac{z}{2t}}(w)$, the residue contributions,
that are nonzero only for $\zeta>0$, can be written as

$\mathcal{R}_{1}^{\psi}(z,t)\sim-\frac{\mathrm{i}}{2}\operatorname{Re}\sum_{j=1}^{4}\int_{-\infty}^{\infty}\frac{\mathrm{d}y}{2\pi}\\\
\left(\mathrm{e}^{2\mathrm{i}ts_{\frac{z}{2t}}(\bar{p}_{j}^{+})-t|s^{\prime\prime}_{\frac{z}{2t}}(\bar{p}_{j}^{+})|y^{2}}\mathrm{e}^{\mathrm{i}\phi^{+}_{j}}\mathrm{e}^{-\mathrm{i}\bar{p}_{j}^{+}}\psi^{R}(2\bar{p}_{j}^{+})-\mathrm{e}^{2\mathrm{i}ts_{-\frac{z}{2t}}(\bar{p}_{j}^{-})-t|s^{\prime\prime}_{-\frac{z}{2t}}(\bar{p}_{j}^{-})|y^{2}}\mathrm{e}^{\mathrm{i}\phi^{-}_{j}}\mathrm{e}^{-\mathrm{i}\bar{p}_{j}^{-}}\psi^{R}(2\bar{p}_{j}^{-})\right),$
(114)

$\mathcal{R}_{2}^{\psi}(z,t)\sim\frac{\mathrm{i}}{2}\operatorname{Im}\sum_{j=1}^{4}\int_{-\infty}^{\infty}\frac{\mathrm{d}y}{2\pi}\biggl{(}\mathrm{e}^{2\mathrm{i}ts_{\frac{z}{2t}}(\bar{p}_{j}^{+})-t|s^{\prime\prime}_{\frac{z}{2t}}(\bar{p}_{j}^{+})|y^{2}}\mathrm{e}^{\mathrm{i}\phi^{+}_{j}}\mathrm{e}^{\mathrm{i}\theta(2\bar{p}^{+}_{j})}\mathrm{e}^{-\mathrm{i}\bar{p}_{j}^{+}}\psi^{R}(2\bar{p}_{j}^{+})+\\\
-\mathrm{e}^{2\mathrm{i}ts_{-\frac{z}{2t}}(\bar{p}_{j}^{-})-t|s^{\prime\prime}_{-\frac{z}{2t}}(\bar{p}_{j}^{-})|y^{2}}\mathrm{e}^{\mathrm{i}\phi^{-}_{j}}\mathrm{e}^{-\mathrm{i}\theta(2\bar{p}^{-}_{j})}\mathrm{e}^{-\mathrm{i}\bar{p}_{j}^{-}}\psi^{R}(2\bar{p}_{j}^{-})\biggr{)},$
(115)

where

$\mathrm{e}^{2\mathrm{i}\phi^{\pm}_{j}-2\mathrm{i}\bar{p}_{j}^{\pm}-\mathrm{i}\frac{\pi}{2}+\mathrm{i}\pi\Theta(-s_{\pm\frac{z}{2t}}^{\prime\prime}(2\bar{p}^{\pm}_{j}))}=-1\qquad\Rightarrow\qquad\phi^{\pm}_{j}=\bar{p}_{j}^{\pm}-\frac{\pi}{4}+\frac{\pi}{2}\Theta(\epsilon^{\prime\prime}(2\bar{p}^{\pm}_{j}))+\pi\kappa^{\pm}_{j},$
(116)

with $\kappa^{\pm}_{j}\in\mathbb{Z}$ determined by the sign of circulation of
the integration contour. Performing the Gaussian integration we get

$\mathcal{R}_{1}^{\psi}(z,t)\sim-\frac{\mathrm{i}}{4\sqrt{\pi
t}}\operatorname{Re}\sum_{j=1}^{4}\left(\frac{\mathrm{e}^{2\mathrm{i}ts_{\frac{z}{2t}}(\bar{p}_{j}^{+})}\mathrm{e}^{\mathrm{i}\phi^{+}_{j}}\mathrm{e}^{-\mathrm{i}\bar{p}_{j}^{+}}\psi^{R}(2\bar{p}_{j}^{+})}{\sqrt{|s^{\prime\prime}_{\frac{z}{2t}}(\bar{p}_{j}^{+})|}}-\frac{\mathrm{e}^{2\mathrm{i}ts_{\frac{z}{2t}}(\bar{p}_{j}^{-})}\mathrm{e}^{\mathrm{i}\phi^{-}_{j}}\mathrm{e}^{-\mathrm{i}\bar{p}_{j}^{-}}\psi^{R}(2\bar{p}_{j}^{-})}{\sqrt{|s^{\prime\prime}_{\frac{z}{2t}}(\bar{p}_{j}^{-})|}}\right),$
(117)

$\mathcal{R}_{2}^{\psi}(z,t)\sim\frac{\mathrm{i}}{4\sqrt{\pi
t}}\operatorname{Im}\sum_{j=1}^{4}\\\
\left(\frac{\mathrm{e}^{2\mathrm{i}ts_{\frac{z}{2t}}(\bar{p}_{j}^{+})}\mathrm{e}^{\mathrm{i}\phi^{+}_{j}}\mathrm{e}^{\mathrm{i}\theta(2\bar{p}^{+}_{j})}\mathrm{e}^{-\mathrm{i}\bar{p}_{j}^{+}}\psi^{R}(2\bar{p}_{j}^{+})}{\sqrt{|s^{\prime\prime}_{\frac{z}{2t}}(\bar{p}_{j}^{+})|}}-\frac{\mathrm{e}^{2\mathrm{i}ts_{-\frac{z}{2t}}(\bar{p}_{j}^{-})}\mathrm{e}^{\mathrm{i}\phi^{-}_{j}}\mathrm{e}^{-\mathrm{i}\theta(2\bar{p}^{-}_{j})}\mathrm{e}^{-\mathrm{i}\bar{p}_{j}^{-}}\psi^{R}(2\bar{p}_{j}^{-})}{\sqrt{|s^{\prime\prime}_{\frac{z}{2t}}(\bar{p}_{j}^{-})|}}\right).$
(118)

At this point let us assign the numbers $\\{1,2,3,4\\}$ to the saddle points
that, for $z=0$, are in $\\{\pi/2,\pi,-\pi/2,0\\}$ respectively (since they
change continuously with $z/t$, the definition is not ambiguous). Then we have
that $\bar{p}^{\pm}_{1}=\bar{p}^{\pm}_{3}+\pi$,
$\bar{p}^{\pm}_{2}=\bar{p}^{\pm}_{4}+\pi$,
$\epsilon^{\prime\prime}(2\bar{p}^{\pm}_{1})<0$,
$\epsilon^{\prime\prime}(2\bar{p}^{\pm}_{2})>0$,
$\kappa^{\pm}_{1}=1=\kappa^{\pm}_{3}$, $\kappa^{\pm}_{2}=0=\kappa^{\pm}_{4}$.
Moreover, we have $\bar{p}_{1}^{+}=\pi-\bar{p}_{1}^{-}$ and
$\bar{p}_{2}^{+}=-\bar{p}_{2}^{-}$, and
$\sqrt{|s^{\prime\prime}_{\zeta}(p)|}=2\sqrt{|\epsilon^{\prime\prime}(2p)|}$.
Everything considered, the residue contribution gives the leading order for
the correlation matrix elements for $\alpha=1$ and $\zeta>0$:

$\displaystyle\left(\Gamma^{z;\psi}_{x}\right)_{11}=-\left(\Gamma^{z;\psi}_{x}\right)_{22}$
$\displaystyle\sim-\frac{\mathrm{i}}{2\sqrt{\pi
t}}\operatorname{Re}\sum_{j=1}^{2}\mathrm{e}^{\mathrm{i}\frac{\pi}{4}(5-2j)}\frac{\mathrm{e}^{2\mathrm{i}ts_{\frac{z}{2t}}(\bar{p}_{j}^{+})}\psi^{R}(2\bar{p}_{j}^{+})}{\sqrt{|\epsilon^{\prime\prime}(2\bar{p}_{j}^{+})|}},$
(119)
$\displaystyle\left(\Gamma^{z;\psi}_{x}\right)_{21}=-\left(\Gamma^{-z;\psi}_{x}\right)_{12}$
$\displaystyle\sim\frac{\mathrm{i}}{2\sqrt{\pi
t}}\operatorname{Im}\sum_{j=1}^{2}\mathrm{e}^{\mathrm{i}\frac{\pi}{4}(5-2j)}\frac{\mathrm{e}^{2\mathrm{i}ts_{\frac{z}{2t}}(\bar{p}_{j}^{+})}\mathrm{e}^{\mathrm{i}\theta(2\bar{p}^{+}_{j})}\psi^{R}(2\bar{p}_{j}^{+})}{\sqrt{|\epsilon^{\prime\prime}(2\bar{p}_{j}^{+})|}}.$

Note that the existence of a non-exponentially-suppressed residue contribution
depends only on the two conditions $x>0$, $|z/t|<2v_{max}$, which means that
we have a residue contribution also outside (to the right) of the lightcone.
For $-v_{max}<\zeta<0$, instead, the residue contribution is zero and, as
already pointed out, the leading order comes from the double integral over
$\mathcal{C}_{1}^{\prime}$ and $\mathcal{C}_{2}^{\prime}$.

###  B.2 Sub-linear scaling

Let us now look at the case in which $z\sim t^{\alpha}$, $\alpha\in(0,1)$,
starting with the case $\zeta>0$. The arguments above for $\alpha=1$ still
hold and the double integral that one gets after the deformation gives the
same $t^{-1}$ contribution, but there is a subtlety in the application of the
saddle-point-method to the residue contribution. The peculiarity is that the
four saddle point $\bar{p}_{j}^{+}$ belong to $\\{\mathrm{e}^{\mathrm{i}\delta
p_{4}},\mathrm{e}^{\mathrm{i}\frac{\pi}{2}+\mathrm{i}\delta
p_{1}},\mathrm{e}^{\mathrm{i}\pi+\mathrm{i}\delta
p_{2}},\mathrm{e}^{-\mathrm{i}\frac{\pi}{2}+\mathrm{i}\delta p_{3}}\\}$, with
$\delta p_{j}\sim t^{\alpha-1}$, and the function $\psi^{R}(p)$ is zero in $0$
and $\pi$. Importantly, this does not depend on the particular choice of the
state that we made: $\psi^{R}(p)$ essentially represents the field
$\psi_{x,t}(p)$ for a homogeneous case, for which it can be shown
$\psi_{x,t}(p)\propto\braket{b(p)b(-p)}$, that equals zero for $p=0$ and
$p=\pi$ for any state. Therefore, when it is computed in the saddle points, it
can be expanded and we get $\psi^{R}(2\bar{p}_{j}^{+})\sim 2\delta
p_{j}\psi^{\prime}(0)$ for $j\in\\{2,4\\}$ and $\psi^{R}(2\bar{p}_{j}^{+})\sim
2\delta p_{j}\psi^{\prime}(\pi)$ for $j\in\\{1,3\\}$; combined with the
$t^{-1/2}$ that one gets from the saddle point we can conclude that the
residue contribution in this case is of order $t^{-3/2+\alpha}$. Note that for
$\alpha<1/2$, this becomes smaller than the double-integral contribution
inside the lightcone, so that, for both rays inside the lightcone, it is
actually the double-integral that gives the leading contribution.

We should point out that in this case the TFIC is special with respect to
generic free chains. As a matter of fact, since in the generic case the
dispersion relation have additional stationary points that are different from
0 and $\pi$, the symmetries of the auxiliary field do not imply that it is
zero in such points, so that the generic decaying of the contribution under
consideration is $t^{-1/2}$ for any observable whose support scales as
$t^{\alpha}$, $\forall\alpha\in(0,1)$. That is the case e.g. of the XY quantum
spin chain with transverse magnetic field when the dispersion relation has
more than two stationary points.

###  B.3 Constant distance

Let us finally discuss what happens for $\alpha=0$. As the case
$\alpha\in(0,1)$, the only thing that is affected is the residue contribution
and only for spin chains whose dispersion relation has stationary points in
$\pi m$, with $m\in\mathbb{Z}$, like the TFIC. We have

$\mathcal{R}_{1}^{\psi}(z,t)=\mathrm{i}\operatorname{Re}\int_{-\pi}^{\pi}\frac{\mathrm{d}p}{2\pi}\psi^{R}(p)\mathrm{e}^{-2\mathrm{i}t\epsilon(p)}\sin(zp)\,,$
(120)
$\mathcal{R}_{2}^{\psi}(z,t)=-\mathrm{i}\operatorname{Im}\int_{-\pi}^{\pi}\frac{\mathrm{d}p}{2\pi}\psi^{R}(p)\mathrm{e}^{-2\mathrm{i}t\epsilon(p)}\sin\left(zp+\theta(p)\right).$
(121)

Now, since the dispersion relation is stationary in $0$ and $\pi$, the
expansion around those points gives

$\mathcal{R}_{1}^{\psi}(z,t)\sim\mathrm{i}\operatorname{Re}\biggl{(}\int_{-\infty}^{+\infty}\frac{\mathrm{d}p}{2\pi}{{\psi}^{R}}^{\prime}(0)\mathrm{e}^{-2\mathrm{i}t\epsilon(0)}\mathrm{e}^{-\mathrm{i}t\epsilon^{\prime\prime}(0)p^{2}}zp^{2}+\\\
+(-1)^{z}\int_{-\infty}^{+\infty}\frac{\mathrm{d}p}{2\pi}{\psi^{R}}^{\prime}(\pi)\mathrm{e}^{-2\mathrm{i}t\epsilon(\pi)}\mathrm{e}^{-\mathrm{i}t\epsilon^{\prime\prime}(\pi)p^{2}}zp^{2}\biggr{)}\\\
\sim\mathrm{i}\operatorname{Re}\biggl{(}\frac{1}{4\sqrt{\pi}}t^{-3/2}z\mathrm{e}^{-\mathrm{i}\frac{\pi}{4}}\left({{\psi}^{R}}^{\prime}(0)\mathrm{e}^{-2\mathrm{i}t\epsilon(0)}(\epsilon^{\prime\prime}(0))^{-3/2}+(-1)^{z}{\psi^{R}}^{\prime}(\pi)\mathrm{e}^{-2\mathrm{i}t\epsilon(\pi)}(\epsilon^{\prime\prime}(\pi))^{-3/2}\right)\biggr{)}\,,$
(122)

and similarly for $\mathcal{R}_{2}^{\psi}$. So the residue contribution in
this case decays as $t^{-3/2}$ and it gives the leading order only to the
right of the lightcone, since inside the lightcone the double integral over
$\mathcal{C}_{1}^{\prime}$ and $\mathcal{C}_{2}^{\prime}$ scales as $t^{-1}$.

###  B.4 Second contribution

In this section we show the deformations envisaged in the main text for the
second term of the integrals (96) representing the contribution of the
auxiliary field to the correlation matrix. They are reported in Fig. 9 and the
same considerations done for the first contribution hold: it is always
possible to deform $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ to
$\mathcal{C}^{\prime}_{1}$ and $\mathcal{C}^{\prime}_{2}$ such that the real
part of the phase is negative almost everywhere; we get a residue contribution
from the deformation only for $x/t>0$, when the two contours have to be fully
exchanged.

$x/t=-1.3,z/t=0.2$$0$$0$$-1$$-1$$+1$$+1$

$x/t=-0.4,z/t=1.8$$0$$0$$-1$$-1$$+1$$+1$

(a)

$x/t=1.3,z/t=0.2$$0$$0$$-1$$-1$$+1$$+1$

$x/t=1.3,z/t=2.2$$0$$0$$-1$$-1$$+1$$+1$

$x/t=0.8,z/t=0.8$$0$$0$$-1$$-1$$+1$$+1$

(b)

$x/t=-0.6,z/t=0.7$$0$$0$$-1$$-1$$+1$$+1$

$x/t=-0.4,z/t=1.0$$0$$0$$-1$$-1$$+1$$+1$

$x/t=0.4,z/t=1.0$$0$$0$$-1$$-1$$+1$$+1$

(c)

Figure 9: Regions in which
$\operatorname{Re}(S_{(x-\frac{z}{2})/t}(a+\mathrm{i}b))>0$ (blue region) and
$\operatorname{Re}(S_{(-x-\frac{z}{2})/t}(a+\mathrm{i}b))>0$ (orange region)
in function of $a$ (horizontal axis) and $b$ (vertical axis). The rest is the
same of Fig. 8.

##  Appendix C Non-Gaussian partitioning protocol

In this appendix, we consider the partitioning protocol (53) between a RITS
and an infinite-temperature state, described by the density matrix

$\hat{\rho}=Z^{-1}\mathrm{e}^{\beta\sum_{l=0}^{\infty}\sigma^{z}_{l}\sigma^{z}_{l+1}},$
(123)

with $\beta>0$. By writing explicitly the correlation functions of Bogoliubov
operators, we give an integral representation of the contribution from the
field $\xi_{x_{1},x_{2},t}(p_{1},p_{2})$ to the spin-spin connected
correlation function (47).

###  C.1 General expression for fermionic 4-point correlations

We can use the Jordan-Wigner transformation (2) to express the 2- and 4-point
correlation functions of Majorana fermions in terms of spin expectation
values. We assume that the state under consideration has the same symmetries
of our example, i.e. that it is invariant under rotations of $\pi$ around the
$x$ and the $y$ axis, and it is invariant under
$\sigma^{\alpha}_{\ell}\rightarrow-\sigma^{\alpha}_{\ell},\forall\ell\in\mathbb{Z}\land\forall\alpha\in\\{x,y\\}$.

Observe first that, whenever we consider the expectation value of a string of
Pauli matrices, it vanishes if any component of the string is a $\sigma^{x}$
or a $\sigma^{y}$, or if the number of $\sigma^{z}$’s is odd. For example,
assuming $m\leq n$, the only way in which

$\braket{a_{2m}a_{2n}}=\Braket{\sigma_{m}^{y}\left(\prod_{j=m}^{n-1}\sigma^{z}_{j}\right)\sigma^{y}_{n}}$
(124)

can satisfy all of the symmetry requirements above, is for $m=n$, while there
is no way that $\braket{a_{2m-1}a_{2n}}$ can satisfy them, leading to

$\braket{a_{i}a_{j}}_{0}=\delta_{ij}.$ (125)

This proves that the correlation matrix for this state is zero. The
computation of the 4-point correlation in Majorana operators is slightly more
complicated, but proceeds on the same lines. It can be shown that it is
nonzero only if two of the indices are odd and the other two are even:

$C_{2m_{1}+j_{1},2m_{2}+j_{2},2m_{3}+j_{3},2m_{4}+j_{4}}=\\\
=\left[\delta_{m_{1}m_{2}}\delta_{m_{3},m_{4}}\braket{\sigma_{m_{1}}^{z}\sigma_{m_{3}}^{z}}_{0}-\frac{1}{2}\left(\delta_{m_{1}m_{3}}\delta_{m_{2}m_{4}}+\delta_{m_{1}m_{4}}\delta_{m_{2}m_{3}}\right)\braket{\sigma_{m_{1}}^{z}\sigma_{m_{2}}^{z}}_{0}\right]\times\\\
\times(\sigma^{y}\otimes\sigma^{y})_{2(j_{1}-1)+j_{3},2(j_{2}-1)+j_{4}}+\\\
+\frac{1}{2}\left(\delta_{m_{1}m_{3}}\delta_{m_{2}m_{4}}-\delta_{m_{1}m_{4}}\delta_{m_{2}m_{3}}\right)\braket{\sigma_{m_{1}}^{z}\sigma_{m_{2}}^{z}}_{0}\times\\\
\times(\mathbb{I}\otimes\mathbb{I}-\sigma^{z}\otimes\sigma^{z}-\sigma^{x}\otimes\sigma^{x})_{2(j_{1}-1)+j_{3},2(j_{2}-1)+j_{4}},$
(126)

for any $m_{1},m_{2},m_{3},m_{4}\in\mathbb{Z}$ and
$j_{1},j_{2},j_{3},j_{4}\in\\{1,2\\}$.

At this point we can use the Bogoliubov transformation to write the Bogoliubov
4-point connected correlation function in terms of Majorana fermions using Eq.
(13):

$\displaystyle\braket{b^{\dagger}(p_{1})b(p_{2})b^{\dagger}(p_{3})b(p_{4})}_{c}=$
$\displaystyle\frac{1}{4}\sum_{\ell_{1},\ell_{2},\ell_{3},\ell_{4}\in\mathbb{Z}}w_{\ell_{1}}(p_{1})w^{*}_{\ell_{2}}(p_{2})w_{\ell_{3}}(p_{3})w^{*}_{\ell_{4}}(p_{4})\braket{a_{\ell_{1}}a_{\ell_{2}}a_{\ell_{3}}a_{\ell_{4}}}_{c},$
(127) $\displaystyle\braket{b(p_{1})b(p_{2})b(p_{3})b(p_{4})}_{c}=$
$\displaystyle\frac{1}{4}\sum_{\ell_{1},\ell_{2},\ell_{3},\ell_{4}\in\mathbb{Z}}w^{*}_{\ell_{1}}(p_{1})w^{*}_{\ell_{2}}(p_{2})w^{*}_{\ell_{3}}(p_{3})w^{*}_{\ell_{4}}(p_{4})\Braket{a_{\ell_{1}}a_{\ell_{2}}a_{\ell_{3}}a_{\ell_{4}}}_{c},$
$\displaystyle\braket{b^{\dagger}(p_{1})b(p_{2})b(p_{3})b(p_{4})}_{c}=$
$\displaystyle\frac{1}{4}\sum_{\ell_{1},\ell_{2},\ell_{3},\ell_{4}\in\mathbb{Z}}w_{\ell_{1}}(p_{1})w^{*}_{\ell_{2}}(p_{2})w^{*}_{\ell_{3}}(p_{3})w^{*}_{\ell_{4}}(p_{4})\Braket{a_{\ell_{1}}a_{\ell_{2}}a_{\ell_{3}}a_{\ell_{4}}}_{c}.$

Now the idea is to plug in the expression of
$C_{\ell_{1},\ell_{2},\ell_{3},\ell_{4}}=\Braket{a_{\ell_{1}}a_{\ell_{2}}a_{\ell_{3}}a_{\ell_{4}}}_{c}$
that we obtained above. Before doing so, it is convenient to split the sum as

$\sum_{\ell_{1},\ell_{2},\ell_{3},\ell_{4}\in\mathbb{Z}}f_{\ell_{1},\ell_{2},\ell_{3},\ell_{4}}C_{\ell_{1},\ell_{2},\ell_{3},\ell_{4}}=\\\
=\sum_{\ell_{1},\ell_{2},\ell_{3},\ell_{4}\in\mathbb{Z}}\biggl{(}f_{2\ell_{1}+1,2\ell_{2}+1,2\ell_{3}+2,2\ell_{4}+2}C_{2\ell_{1}+1,2\ell_{2}+1,2\ell_{3}+2,2\ell_{4}+2}+\\\
+f_{2\ell_{1}+1,2\ell_{2}+2,2\ell_{3}+1,2\ell_{4}+2}C_{2\ell_{1}+1,2\ell_{2}+2,2\ell_{3}+1,2\ell_{4}+2}\\\
+f_{2\ell_{1}+1,2\ell_{2}+2,2\ell_{3}+2,2\ell_{4}+1}C_{2\ell_{1}+1,2\ell_{2}+2,2\ell_{3}+2,2\ell_{4}+1}+\\\
+f_{2\ell_{1}+2,2\ell_{2}+1,2\ell_{3}+1,2\ell_{4}+2}C_{2\ell_{1}+2,2\ell_{2}+1,2\ell_{3}+1,2\ell_{4}+2}+\\\
+f_{2\ell_{1}+2,2\ell_{2}+1,2\ell_{3}+2,2\ell_{4}+1}C_{2\ell_{1}+2,2\ell_{2}+1,2\ell_{3}+2,2\ell_{4}+1}+\\\
+f_{2\ell_{1}+2,2\ell_{2}+2,2\ell_{3}+1,2\ell_{4}+1}C_{2\ell_{1}+2,2\ell_{2}+2,2\ell_{3}+1,2\ell_{4}+1}\biggr{)}=\\\
=\sum_{\ell_{1},\ell_{2}\in\mathbb{Z}}\biggl{(}f_{2\ell_{1}+1,2\ell_{2}+1,2\ell_{1}+2,2\ell_{2}+2}-f_{2\ell_{1}+1,2\ell_{1}+2,2\ell_{2}+1,2\ell_{2}+2}+f_{2\ell_{1}+1,2\ell_{1}+2,2\ell_{2}+2,2\ell_{2}+1}+\\\
+f_{2\ell_{1}+2,2\ell_{1}+1,2\ell_{2}+1,2\ell_{2}+2}-f_{2\ell_{1}+2,2\ell_{1}+1,2\ell_{2}+2,2\ell_{2}+1}+f_{2\ell_{1}+2,2\ell_{2}+2,2\ell_{1}+1,2\ell_{2}+1}+\\\
-f_{2\ell_{1}+1,2\ell_{2}+1,2\ell_{2}+2,2\ell_{1}+2}+f_{2\ell_{1}+1,2\ell_{2}+2,2\ell_{2}+1,2\ell_{1}+2}-f_{2\ell_{1}+1,2\ell_{2}+2,2\ell_{1}+2,2\ell_{2}+1}+\\\
-f_{2\ell_{2}+2,2\ell_{1}+1,2\ell_{2}+1,2\ell_{1}+2}+f_{2\ell_{2}+2,2\ell_{1}+1,2\ell_{1}+2,2\ell_{2}+1}-f_{2\ell_{2}+2,2\ell_{1}+2,2\ell_{1}+1,2\ell_{2}+1}\biggr{)}\braket{\sigma_{\ell_{1}}^{z}\sigma_{\ell_{2}}^{z}}_{0},$
(128)

where we have used that, in order to have a nonzero result, $C$ needs to have
two even and two odd indices, then we have used its explicit expression. Using
the explicit expression of the eigenvectors $w(p)$ – see Eq. (11) – and
introducing the (discrete) Fourier transform of the spin-spin correlation
function

$T(p,q):=\sum_{m,n\in\mathbb{Z}}\Braket{\sigma_{m}^{z}\sigma_{n}^{z}}_{0}\mathrm{e}^{-\mathrm{i}pm}\mathrm{e}^{-\mathrm{i}qn},$
(129)

we finally have

$\braket{b^{\dagger}(p_{1})b(p_{2})b^{\dagger}(p_{3})b(p_{4})}_{c}=\frac{1}{4}T(-p_{1}-p_{3},p_{2}+p_{4})\sin\left(\tfrac{\theta(p_{1})-\theta(p_{3})}{2}\right)\sin\left(\tfrac{\theta(p_{2})-\theta(p_{4})}{2}\right)\\\
+\frac{1}{4}T(-p_{1}+p_{2},-p_{3}+p_{4})\cos\left(\tfrac{\theta(p_{1})+\theta(p_{2})}{2}\right)\cos\left(\tfrac{\theta(p_{3})+\theta(p_{4})}{2}\right)\\\
-\frac{1}{4}T(-p_{1}+p_{4},p_{2}-p_{3})\cos\left(\tfrac{\theta(p_{1})+\theta(p_{4})}{2}\right)\cos\left(\tfrac{\theta(p_{2})+\theta(p_{3})}{2}\right),$
(130)

$\braket{b(p_{1})b(p_{2})b(p_{3})b(p_{4})}_{c}=\frac{1}{4}T(p_{1}+p_{3},p_{2}+p_{4})\sin\left(\tfrac{\theta(p_{1})-\theta(p_{3})}{2}\right)\sin\left(\tfrac{\theta(p_{2})-\theta(p_{4})}{2}\right)\\\
-\frac{1}{4}T(p_{1}+p_{2},p_{3}+p_{4})\sin\left(\tfrac{\theta(p_{1})-\theta(p_{2})}{2}\right)\sin\left(\tfrac{\theta(p_{3})-\theta(p_{4})}{2}\right)\\\
-\frac{1}{4}T(p_{1}+p_{4},p_{2}+p_{3})\sin\left(\tfrac{\theta(p_{1})-\theta(p_{4})}{2}\right)\sin\left(\tfrac{\theta(p_{2})-\theta(p_{3})}{2}\right),$
(131)

$\braket{b^{\dagger}(p_{1})b(p_{2})b(p_{3})b(p_{4})}_{c}=\frac{\mathrm{i}}{4}T(-p_{1}+p_{2},p_{3}+p_{4})\cos\left(\tfrac{\theta(p_{1})+\theta(p_{2})}{2}\right)\sin\left(\tfrac{\theta(p_{3})-\theta(p_{4})}{2}\right)\\\
-\frac{\mathrm{i}}{4}T(-p_{1}+p_{3},p_{2}+p_{4})\cos\left(\tfrac{\theta(p_{1})+\theta(p_{3})}{2}\right)\sin\left(\tfrac{\theta(p_{2})-\theta(p_{4})}{2}\right)\\\
+\frac{\mathrm{i}}{4}T(-p_{1}+p_{4},p_{2}+p_{3})\cos\left(\tfrac{\theta(p_{1})+\theta(p_{4})}{2}\right)\sin\left(\tfrac{\theta(p_{2})-\theta(p_{3})}{2}\right)\,.$
(132)

###  C.2 Integral representation of the spin-spin correlation function

Let us now specialize to our case of interest. Using

$\braket{\sigma_{m}^{z}\sigma_{n}^{z}}_{0}=(\tanh\beta)^{|m-n|}\Theta_{m}\Theta_{n}\,,$
(133)

we can compute explicitly the function $T(p,q)$ of the previous section:

$T(p,q)=\sum_{m,n\in\mathbb{Z}}\mathrm{e}^{-\mathrm{i}pm}\mathrm{e}^{-\mathrm{i}qn}(\tanh\beta)^{|m-n|}\Theta_{m}\Theta_{n}=\int_{-\pi}^{\pi}\frac{\mathrm{d}k}{2\pi}\frac{T_{\beta}(k)}{(1-\mathrm{e}^{-\mathrm{i}(k+p-\mathrm{i}0)})(1-\mathrm{e}^{-\mathrm{i}(q-k-\mathrm{i}0)})},$
(134)

where

$T_{\beta}(k):=\sum_{j\in\mathbb{Z}}\left(\tanh\beta\right)^{|j|}\mathrm{e}^{-\mathrm{i}kj}=\frac{1}{\cosh^{2}(\beta)}\frac{1}{\tanh^{2}(\beta)-2\tanh(\beta)\cos(k)+1}$
(135)

is the (discrete) Fourier transform of the homogeneous spin-spin correlation
function; note that a straightforward generalization of our problem to a
different initial state is obtained simply by changing $T_{\beta}(k)$. Now we
can compute the initial value of all the dynamical fields; for example:

$\xi^{phys}_{x_{1},x_{2},0}(p_{1},p_{2})=\int_{-\pi}^{+\pi}\frac{\mathrm{d}^{2}q\mathrm{d}k}{(2\pi)^{3}}\mathrm{e}^{2\mathrm{i}x_{1}q_{1}}\mathrm{e}^{2\mathrm{i}x_{2}q_{2}}\frac{T_{\beta}(k)}{4}\times\\\
\times\biggl{(}\frac{\sin\left(\frac{\theta(p_{1}-q_{1})-\theta(p_{2}-q_{2})}{2}\right)\sin\left(\frac{\theta(p_{1}+q_{1})-\theta(p_{2}+q_{2})}{2}\right)}{(1-\mathrm{e}^{-\mathrm{i}(k-p_{1}+q_{1}-p_{2}+q_{2}-\mathrm{i}0)})(1-\mathrm{e}^{-\mathrm{i}(p_{1}+q_{1}+p_{2}+q_{2}-k-\mathrm{i}0)})}\\\
+\frac{\cos\left(\frac{\theta(p_{1}-q_{1})+\theta(p_{1}+q_{1})}{2}\right)\cos\left(\frac{\theta(p_{2}-q_{2})+\theta(p_{2}+q_{2})}{2}\right)}{(1-\mathrm{e}^{-\mathrm{i}(k+2q_{1}-\mathrm{i}0)})(1-\mathrm{e}^{-\mathrm{i}(2q_{2}-k-\mathrm{i}0)})}+\\\
-\frac{\cos\left(\frac{\theta(p_{1}-q_{1})+\theta(p_{2}+q_{2})}{2}\right)\cos\left(\frac{\theta(p_{1}+q_{1})+\theta(p_{2}-q_{2})}{2}\right)}{(1-\mathrm{e}^{-\mathrm{i}(k-p_{1}+q_{1}+p_{2}+q_{2}-\mathrm{i}0)})(1-\mathrm{e}^{-\mathrm{i}(p_{1}+q_{1}-p_{2}+q_{2}-k-\mathrm{i}0)})}\biggr{)}\\\
=\int_{-\pi}^{+\pi}\frac{\mathrm{d}^{2}q\mathrm{d}k}{(2\pi)^{3}}\frac{\mathrm{e}^{2\mathrm{i}x_{1}q_{1}}\mathrm{e}^{2\mathrm{i}x_{2}q_{2}}}{4(1-\mathrm{e}^{-\mathrm{i}(2q_{1}+k-\mathrm{i}0)})(1-\mathrm{e}^{-\mathrm{i}(2q_{2}-k-\mathrm{i}0)})}F_{\beta}(p_{1},p_{2},q_{1},q_{2},k),$
(136)

with

$F_{\beta}(p_{1},p_{2},q_{1},q_{2},k):=\\\
T_{\beta}(k+q_{1}+p_{1}+p_{2}-q_{2})\sin\left(\tfrac{\theta(p_{1}-q_{1})-\theta(p_{2}-q_{2})}{2}\right)\sin\left(\tfrac{\theta(p_{1}+q_{1})-\theta(p_{2}+q_{2})}{2}\right)+\\\
+T_{\beta}(k)\cos\left(\tfrac{\theta(p_{1}-q_{1})+\theta(p_{1}+q_{1})}{2}\right)\cos\left(\tfrac{\theta(p_{2}-q_{2})+\theta(p_{2}+q_{2})}{2}\right)+\\\
-T_{\beta}(k+p_{1}+q_{1}-p_{2}-q_{2})\cos\left(\tfrac{\theta(p_{1}-q_{1})+\theta(p_{2}+q_{2})}{2}\right)\cos\left(\tfrac{\theta(p_{1}+q_{1})+\theta(p_{2}-q_{2})}{2}\right).$
(137)

We can now use Eq. (46) to get the time-dependent expression of the field

$\xi_{x_{1},x_{2},t}(p_{1},p_{2})=\\\
\int_{-\pi}^{+\pi}\frac{\mathrm{d}^{2}q\mathrm{d}k}{(2\pi)^{3}}\frac{\mathrm{e}^{2\mathrm{i}x_{1}q_{1}}\mathrm{e}^{2\mathrm{i}x_{2}q_{2}}\mathrm{e}^{-\mathrm{i}t(\epsilon(p_{1}+q_{1})-\epsilon(p_{1}-q_{1})+\epsilon(p_{2}+q_{2})-\epsilon(p_{2}-q_{2}))}}{4(1-\mathrm{e}^{-\mathrm{i}(2q_{1}+k-\mathrm{i}0)})(1-\mathrm{e}^{-\mathrm{i}(2q_{2}-k-\mathrm{i}0)})}F_{\beta}(p_{1},p_{2},q_{1},q_{2},k),$
(138)

and finally compute the contribution of the field $\xi$ to the spin-spin
correlation using its expression (48), the parametrization of the symbol (39)
and the inversion relation (38):

$\displaystyle
S^{z,\xi}_{mn}(t)=-4\int_{-\pi}^{\pi}\frac{\mathrm{d}^{2}p\mathrm{d}^{2}k}{(2\pi)^{4}}\sum_{y_{1},y_{2}\in\mathbb{Z}/2}\mathrm{e}^{2\mathrm{i}k_{1}(m-1-y_{1})}\mathrm{e}^{2\mathrm{i}k_{2}(n-1-y_{2})}\xi_{y_{1},y_{2},t;ee}(p_{1},p_{2})\times$
(139)
$\displaystyle\times\left(\mathrm{e}^{-\mathrm{i}\frac{\theta(p_{1}+k_{1})}{2}\sigma^{z}}\otimes\mathrm{e}^{-\mathrm{i}\frac{\theta(p_{2}+k_{2})}{2}\sigma^{z}}\right)\left(\sigma^{y}\otimes\sigma^{y}\right)\left(\mathrm{e}^{\mathrm{i}\frac{\theta(p_{1}-k_{1})}{2}\sigma^{z}}\otimes\mathrm{e}^{\mathrm{i}\frac{\theta(p_{2}-k_{2})}{2}\sigma^{z}}\right)|_{1,4}$
$\displaystyle=4\int_{-\pi}^{\pi}\frac{\mathrm{d}^{2}p\mathrm{d}^{2}k}{(2\pi)^{2}4}\sum_{y_{1},y_{2}\in\mathbb{Z}/2}\mathrm{e}^{2\mathrm{i}k_{1}(m-1-y_{1})}\mathrm{e}^{2\mathrm{i}k_{2}(n-1-y_{2})}\xi_{y_{1},y_{2},t;ee}(p_{1},p_{2})\times$
$\displaystyle\times\mathrm{e}^{-\mathrm{i}\frac{\theta(p_{1}-k_{1})+\theta(p_{1}+k_{1})+\theta(p_{2}-k_{2})+\theta(p_{2}+k_{2})}{2}}$
$\displaystyle=4\int_{-\pi}^{\pi}\frac{\mathrm{d}^{2}p\mathrm{d}^{2}k}{(2\pi)^{4}}\sum_{y_{1},y_{2}\in\mathbb{Z}/2}\mathrm{e}^{2\mathrm{i}k_{1}(m-1-y_{1})}\mathrm{e}^{2\mathrm{i}k_{2}(n-1-y_{2})}\xi_{y_{1},y_{2},t}(p_{1},p_{2})$
$\displaystyle\cos\left(\frac{\theta(p_{1}-k_{1})+\theta(p_{1}+k_{1})}{2}\right)\cos\left(\frac{\theta(p_{2}-k_{2})+\theta(p_{2}+k_{2})}{2}\right)$
$\displaystyle=\int_{-\pi}^{\pi}\frac{\mathrm{d}^{2}p\mathrm{d}^{2}q\mathrm{d}k}{(2\pi)^{5}}\frac{\mathrm{e}^{2\mathrm{i}mq_{1}}\mathrm{e}^{2\mathrm{i}nq_{2}}\mathrm{e}^{-\mathrm{i}t(\epsilon(p_{1}+q_{1})-\epsilon(p_{1}-q_{1})+\epsilon(p_{2}+q_{2})-\epsilon(p_{2}-q_{2}))}}{(1-\mathrm{e}^{-\mathrm{i}(2q_{1}+k-\mathrm{i}0)})(1-\mathrm{e}^{-\mathrm{i}(2q_{2}-k-\mathrm{i}0)})}\mathrm{e}^{-2\mathrm{i}q_{1}}\mathrm{e}^{-2\mathrm{i}q_{2}}$
$\displaystyle\cos\left(\tfrac{\theta(p_{1}-q_{1})+\theta(p_{1}+q_{1})}{2}\right)\cos\left(\tfrac{\theta(p_{2}-q_{2})+\theta(p_{2}+q_{2})}{2}\right)F_{\beta}(p_{1},p_{2},q_{1},q_{2},k),$

where we used the Moyal product’s integral representation (31) adapted to the
symbol $\tilde{C}$. We then change variables as

$p_{1}\rightarrow p_{1}+k/2\,,\qquad p_{2}\rightarrow p_{2}+q_{2}\,,\qquad
q_{1}\rightarrow q_{1}+k/2\,,\qquad q_{2}\rightarrow q_{2}-q_{1}\,,\qquad
k\rightarrow-2q_{1}\,,$ (140)

obtaining

$\displaystyle S^{z,\xi}_{mn}(t)=$ (141)
$\displaystyle\int_{-\pi}^{\pi}\frac{\mathrm{d}^{2}p\mathrm{d}^{2}q\mathrm{d}k}{(2\pi)^{5}}\frac{\mathrm{e}^{2\mathrm{i}q_{1}(m-n)}\mathrm{e}^{\mathrm{i}mk}\mathrm{e}^{2\mathrm{i}nq_{2}}\mathrm{e}^{-\mathrm{i}t(\epsilon(p_{1}+q_{1}+k)-\epsilon(p_{1}-q_{1})+\epsilon(p_{2}-q_{1}+2q_{2})-\epsilon(p_{2}+q_{1}))}}{(1-\mathrm{e}^{-\mathrm{i}(k-\mathrm{i}0)})(1-\mathrm{e}^{-\mathrm{i}(2q_{2}-\mathrm{i}0)})}\mathrm{e}^{-\mathrm{i}k}$
$\displaystyle\mathrm{e}^{-2\mathrm{i}q_{2}}\cos\left(\tfrac{\theta(p_{1}-q_{1})+\theta(p_{1}+q_{1}+k)}{2}\right)\cos\left(\tfrac{\theta(p_{2}+q_{1})+\theta(p_{2}-q_{1}+2q_{2})}{2}\right)G_{\beta}(p_{1},p_{2},q_{1},2q_{2},k),$

where

$G_{\beta}(p_{1},p_{2},q_{1},q_{2},k):=\\\
T_{\beta}(k+p_{1}+p_{2})\sin\left(\tfrac{\theta(p_{1}-q_{1})-\theta(p_{2}+q_{1})}{2}\right)\sin\left(\tfrac{\theta(p_{1}+q_{1}+k)-\theta(p_{2}-q_{1}+q_{2})}{2}\right)+\\\
+T_{\beta}(2q_{1})\cos\left(\tfrac{\theta(p_{1}-q_{1})+\theta(p_{1}+q_{1}+k)}{2}\right)\cos\left(\tfrac{\theta(p_{2}+q_{1})+\theta(p_{2}-q_{1}+q_{2})}{2}\right)+\\\
-T_{\beta}(k+p_{1}-p_{2}-q_{2})\cos\left(\tfrac{\theta(p_{1}-q_{1})+\theta(p_{2}-q_{1}+q_{2})}{2}\right)\cos\left(\tfrac{\theta(p_{1}+q_{1}+k)+\theta(p_{2}+q_{1})}{2}\right).$
(142)

Note that $G_{\beta}$ is periodic in all its entries with period $2\pi$. In
the same way, we could derive an integral representation of the contribution
of the other non-Gaussian fields. However, as discussed in the main text,
those fields give a sub-leading contribution.

##  Appendix D Numerical simulations

Numerical simulations for Figures (3), (4) and (6) have been performed
considering finite-size chains of length $2L$, in which $\mathcal{H}$ is a
$2L\times 2L$ matrix. The part of state entering the expression of the spin-
spin connected correlation (48) is described by the objects $\Gamma$ and $C$
defined in the main text, which have dimensions $2L\times 2L$ and $2L\times
2L\times 2L\times 2L$ respectively. Time-evolution has been done exactly using

$a_{i}(t)=\sum_{j=1}^{2L}(\mathrm{e}^{-\mathrm{i}t\mathcal{H}})_{ij}a_{j}(0),$
(143)

which implies

$\Gamma_{m,n}(t)=\sum_{m^{\prime},n^{\prime}=1}^{2L}(\mathrm{e}^{-\mathrm{i}t\mathcal{H}})_{m,m^{\prime}}(\mathrm{e}^{-\mathrm{i}t\mathcal{H}})_{n,n^{\prime}}\Gamma_{m^{\prime},n^{\prime}}(0)$
(144)

and

$C_{m,n,r,s}(t)=\sum_{m^{\prime},n^{\prime},r^{\prime},s^{\prime}=1}^{2L}(\mathrm{e}^{-\mathrm{i}t\mathcal{H}})_{m,m^{\prime}}(\mathrm{e}^{-\mathrm{i}t\mathcal{H}})_{n,n^{\prime}}(\mathrm{e}^{-\mathrm{i}t\mathcal{H}})_{r,r^{\prime}}(\mathrm{e}^{-\mathrm{i}t\mathcal{H}})_{s,s^{\prime}}C_{m^{\prime},n^{\prime},r^{\prime},s^{\prime}}(0).$
(145)

We always stop the simulation at a time such that the finite size effects have
only exponentially-small consequences on the expectation values that we are
looking at, which is always possible thanks to the Lieb-Robinson bounds. In
the case of Gaussian protocols, the initial state has been actually considered
sharp, and not smooth as described in Sections 5.2.1 and 5.2.2, since it is
easier to construct. Therefore the state that we are simulating differs from
the one described in the main text around the junction. However, the
homogeneous states associated with the two partitioning protocols are the same
and the difference arguably induces only sub-leading effects. That is why such
a simulation can still be used to check the leading-order predictions.

## References

  * [1] B. Bertini, M. Collura, J. De Nardis and M. Fagotti,  _Transport in out-of-equilibrium $xxz$ chains: Exact profiles of charges and currents_,  Phys. Rev. Lett. 117, 207201 (2016),  10.1103/PhysRevLett.117.207201.
  * [2] O. A. Castro-Alvaredo, B. Doyon and T. Yoshimura,  _Emergent hydrodynamics in integrable quantum systems out of equilibrium_ ,  Phys. Rev. X 6, 041065 (2016),  10.1103/PhysRevX.6.041065.
  * [3] B. Doyon and T. Yoshimura,  _A note on generalized hydrodynamics: inhomogeneous fields and other concepts_ ,  SciPost Phys. 2, 014 (2017),  10.21468/SciPostPhys.2.2.014.
  * [4] L. Piroli, J. De Nardis, M. Collura, B. Bertini and M. Fagotti,  _Transport in out-of-equilibrium xxz chains: Nonballistic behavior and correlation functions_ ,  Phys. Rev. B 96, 115124 (2017),  10.1103/PhysRevB.96.115124.
  * [5] E. Ilievski and J. De Nardis,  _Microscopic origin of ideal conductivity in integrable quantum models_ ,  Phys. Rev. Lett. 119, 020602 (2017),  10.1103/PhysRevLett.119.020602.
  * [6] E. Ilievski and J. De Nardis,  _Ballistic transport in the one-dimensional hubbard model: The hydrodynamic approach_ ,  Phys. Rev. B 96, 081118 (2017),  10.1103/PhysRevB.96.081118.
  * [7] G. Misguich, K. Mallick and P. L. Krapivsky,  _Dynamics of the spin- $\frac{1}{2}$ heisenberg chain initialized in a domain-wall state_,  Phys. Rev. B 96, 195151 (2017),  10.1103/PhysRevB.96.195151.
  * [8] J.-M. Stéphan,  _Return probability after a quench from a domain wall initial state in the spin-1/2 xxz chain_ ,  Journal of Statistical Mechanics: Theory and Experiment 2017(10), 103108 (2017),  10.1088/1742-5468/aa8c19.
  * [9] M. Ljubotina, M. Žnidarič and T. Prosen,  _Spin diffusion from an inhomogeneous quench in an integrable system_ ,  Nature Communications 8(1), 16117 (2017),  10.1038/ncomms16117.
  * [10] V. B. Bulchandani, R. Vasseur, C. Karrasch and J. E. Moore,  _Bethe-boltzmann hydrodynamics and spin transport in the xxz chain_ ,  Phys. Rev. B 97, 045407 (2018),  10.1103/PhysRevB.97.045407.
  * [11] J. De Nardis, D. Bernard and B. Doyon,  _Hydrodynamic diffusion in integrable systems_ ,  Phys. Rev. Lett. 121, 160603 (2018),  10.1103/PhysRevLett.121.160603.
  * [12] M. Collura, A. De Luca and J. Viti,  _Analytic solution of the domain-wall nonequilibrium stationary state_ ,  Phys. Rev. B 97, 081111 (2018),  10.1103/PhysRevB.97.081111.
  * [13] B. Doyon, T. Yoshimura and J.-S. Caux,  _Soliton gases and generalized hydrodynamics_ ,  Phys. Rev. Lett. 120, 045301 (2018),  10.1103/PhysRevLett.120.045301.
  * [14] B. Doyon,  _Exact large-scale correlations in integrable systems out of equilibrium_ ,  SciPost Phys. 5, 054 (2018),  10.21468/SciPostPhys.5.5.054.
  * [15] A. Bastianello, B. Doyon, G. Watts and T. Yoshimura,  _Generalized hydrodynamics of classical integrable field theory: the sinh-Gordon model_ ,  SciPost Phys. 4, 045 (2018),  10.21468/SciPostPhys.4.6.045.
  * [16] B. Bertini and L. Piroli,  _Low-temperature transport in out-of-equilibrium xxz chains_ ,  Journal of Statistical Mechanics: Theory and Experiment 2018(3), 033104 (2018),  10.1088/1742-5468/aab04b.
  * [17] B. Bertini, L. Piroli and P. Calabrese,  _Universal broadening of the light cone in low-temperature transport_ ,  Phys. Rev. Lett. 120, 176801 (2018),  10.1103/PhysRevLett.120.176801.
  * [18] L. Mazza, J. Viti, M. Carrega, D. Rossini and A. De Luca,  _Energy transport in an integrable parafermionic chain via generalized hydrodynamics_ ,  Phys. Rev. B 98, 075421 (2018),  10.1103/PhysRevB.98.075421.
  * [19] V. Alba,  _Towards a generalized hydrodynamics description of rényi entropies in integrable systems_ ,  Phys. Rev. B 99, 045150 (2019),  10.1103/PhysRevB.99.045150.
  * [20] J. D. Nardis, D. Bernard and B. Doyon,  _Diffusion in generalized hydrodynamics and quasiparticle scattering_ ,  SciPost Phys. 6, 049 (2019),  10.21468/SciPostPhys.6.4.049.
  * [21] B. Doyon,  _Generalized hydrodynamics of the classical toda system_ ,  Journal of Mathematical Physics 60(7), 073302 (2019),  10.1063/1.5096892,  https://doi.org/10.1063/1.5096892.
  * [22] V. B. Bulchandani, X. Cao and J. E. Moore,  _Kinetic theory of quantum and classical toda lattices_ ,  Journal of Physics A: Mathematical and Theoretical 52(33), 33LT01 (2019),  10.1088/1751-8121/ab2cf0.
  * [23] A. Bastianello, V. Alba and J.-S. Caux,  _Generalized hydrodynamics with space-time inhomogeneous interactions_ ,  Phys. Rev. Lett. 123, 130602 (2019),  10.1103/PhysRevLett.123.130602.
  * [24] U. Agrawal, S. Gopalakrishnan and R. Vasseur,  _Generalized hydrodynamics, quasiparticle diffusion, and anomalous local relaxation in random integrable spin chains_ ,  Phys. Rev. B 99, 174203 (2019),  10.1103/PhysRevB.99.174203.
  * [25] V. Alba, B. Bertini and M. Fagotti,  _Entanglement evolution and generalised hydrodynamics: interacting integrable systems_ ,  SciPost Phys. 7, 005 (2019),  10.21468/SciPostPhys.7.1.005.
  * [26] M. Mestyán, B. Bertini, L. Piroli and P. Calabrese,  _Spin-charge separation effects in the low-temperature transport of one-dimensional fermi gases_ ,  Phys. Rev. B 99, 014305 (2019),  10.1103/PhysRevB.99.014305.
  * [27] B. Bertini, L. Piroli and M. Kormos,  _Transport in the sine-gordon field theory: From generalized hydrodynamics to semiclassics_ ,  Phys. Rev. B 100, 035108 (2019),  10.1103/PhysRevB.100.035108.
  * [28] M. Mestyán and V. Alba,  _Molecular dynamics simulation of entanglement spreading in generalized hydrodynamics_ ,  SciPost Phys. 8, 055 (2020),  10.21468/SciPostPhys.8.4.055.
  * [29] M. Collura, A. De Luca, P. Calabrese and J. Dubail,  _Domain wall melting in the spin- $\frac{1}{2}$ xxz spin chain: Emergent luttinger liquid with a fractal quasiparticle charge_,  Phys. Rev. B 102, 180409 (2020),  10.1103/PhysRevB.102.180409.
  * [30] B. Bertini, F. Heidrich-Meisner, C. Karrasch, T. Prosen, R. Steinigeweg and M. Žnidarič,  _Finite-temperature transport in one-dimensional quantum lattice models_ ,  Rev. Mod. Phys. 93, 025003 (2021),  10.1103/RevModPhys.93.025003.
  * [31] A. Bastianello, B. Bertini, B. Doyon and R. Vasseur,  _Introduction to the special issue on emergent hydrodynamics in integrable many-body systems_ ,  Journal of Statistical Mechanics: Theory and Experiment 2022(1), 014001 (2022),  10.1088/1742-5468/ac3e6a.
  * [32] C. N. Yang and C. P. Yang,  _Thermodynamics of a one‐dimensional system of bosons with repulsive delta‐function interaction_ ,  Journal of Mathematical Physics 10(7), 1115 (1969),  10.1063/1.1664947,  https://doi.org/10.1063/1.1664947.
  * [33] M. Takahashi,  _Thermodynamics of One-Dimensional Solvable Models_ , chap. Finite temperature integral equations for un-nested systems,  Cambridge University Press (1999).
  * [34] J.-S. Caux,  _The quench action_ ,  Journal of Statistical Mechanics: Theory and Experiment 2016(6), 064006 (2016),  10.1088/1742-5468/2016/06/064006.
  * [35] V. Alba, B. Bertini, M. Fagotti, L. Piroli and P. Ruggiero,  _Generalized-hydrodynamic approach to inhomogeneous quenches: correlations, entanglement and quantum effects_ ,  Journal of Statistical Mechanics: Theory and Experiment 2021(11), 114004 (2021),  10.1088/1742-5468/ac257d.
  * [36] M. Schemmer, I. Bouchoule, B. Doyon and J. Dubail,  _Generalized hydrodynamics on an atom chip_ ,  Phys. Rev. Lett. 122, 090601 (2019),  10.1103/PhysRevLett.122.090601.
  * [37] N. Malvania, Y. Zhang, Y. Le, J. Dubail, M. Rigol and D. S. Weiss,  _Generalized hydrodynamics in strongly interacting 1d bose gases_ ,  Science 373(6559), 1129 (2021),  10.1126/science.abf0147,  https://www.science.org/doi/pdf/10.1126/science.abf0147.
  * [38] F. Møller, C. Li, I. Mazets, H.-P. Stimming, T. Zhou, Z. Zhu, X. Chen and J. Schmiedmayer,  _Extension of the generalized hydrodynamics to the dimensional crossover regime_ ,  Phys. Rev. Lett. 126, 090602 (2021),  10.1103/PhysRevLett.126.090602.
  * [39] B. Bertini and M. Fagotti,  _Determination of the nonequilibrium steady state emerging from a defect_ ,  Phys. Rev. Lett. 117, 130402 (2016),  10.1103/PhysRevLett.117.130402.
  * [40] M. Fagotti,  _Higher-order generalized hydrodynamics in one dimension: The noninteracting test_ ,  Phys. Rev. B 96, 220302 (2017),  10.1103/PhysRevB.96.220302.
  * [41] M. Fagotti,  _Locally quasi-stationary states in noninteracting spin chains_ ,  SciPost Phys. 8, 48 (2020),  10.21468/SciPostPhys.8.3.048.
  * [42] V. Eisler and Z. Rácz,  _Full counting statistics in a propagating quantum front and random matrix spectra_ ,  Phys. Rev. Lett. 110, 060602 (2013),  10.1103/PhysRevLett.110.060602.
  * [43] J. Viti, J.-M. Stéphan, J. Dubail and M. Haque,  _Inhomogeneous quenches in a free fermionic chain: Exact results_ ,  Europhysics Letters 115(4), 40011 (2016),  10.1209/0295-5075/115/40011.
  * [44] M. Kormos,  _Inhomogeneous quenches in the transverse field Ising chain: scaling and front dynamics_ ,  SciPost Phys. 3, 020 (2017),  10.21468/SciPostPhys.3.3.020.
  * [45] V. Eisler,  _Universality in the full counting statistics of trapped fermions_ ,  Phys. Rev. Lett. 111, 080402 (2013),  10.1103/PhysRevLett.111.080402.
  * [46] D. S. Dean, P. Le Doussal, S. N. Majumdar and G. Schehr,  _Finite-temperature free fermions and the kardar-parisi-zhang equation at finite time_ ,  Phys. Rev. Lett. 114, 110402 (2015),  10.1103/PhysRevLett.114.110402.
  * [47] D. S. Dean, P. L. Doussal, S. N. Majumdar and G. Schehr,  _Universal ground-state properties of free fermions in a d-dimensional trap_ ,  Europhysics Letters 112(6), 60001 (2015),  10.1209/0295-5075/112/60001.
  * [48] E. Wigner,  _On the quantum correction for thermodynamic equilibrium_ ,  Phys. Rev. 40, 749 (1932),  10.1103/PhysRev.40.749.
  * [49] J. E. Moyal and M. S. Bartlett,  _Quantum mechanics as a statistical theory_ ,  Proceedings of the Cambridge Philosophical Society 45(1), 99 (1949),  10.1017/S0305004100000487.
  * [50] E. Bettelheim, A. G. Abanov and P. Wiegmann,  _Nonlinear quantum shock waves in fractional quantum hall edge states_ ,  Phys. Rev. Lett. 97, 246401 (2006),  10.1103/PhysRevLett.97.246401.
  * [51] E. Bettelheim and L. Glazman,  _Quantum ripples over a semiclassical shock_ ,  Phys. Rev. Lett. 109, 260602 (2012),  10.1103/PhysRevLett.109.260602.
  * [52] I. V. Protopopov, D. B. Gutman, P. Schmitteckert and A. D. Mirlin,  _Dynamics of waves in one-dimensional electron systems: Density oscillations driven by population inversion_ ,  Phys. Rev. B 87, 045112 (2013),  10.1103/PhysRevB.87.045112.
  * [53] M. Gluza, J. Eisert and T. Farrelly,  _Equilibration towards generalized Gibbs ensembles in non-interacting theories_ ,  SciPost Phys. 7, 38 (2019),  10.21468/SciPostPhys.7.3.038.
  * [54] C. Murthy and M. Srednicki,  _Relaxation to gaussian and generalized gibbs states in systems of particles with quadratic hamiltonians_ ,  Phys. Rev. E 100, 012146 (2019),  10.1103/PhysRevE.100.012146.
  * [55] E. H. Lieb and D. W. Robinson,  _The finite group velocity of quantum spin systems_ ,  Commun. Math. Phys. 28, 251 (1972),  10.1007/BF01645779.
  * [56] P. Ruggiero, P. Calabrese, B. Doyon and J. Dubail,  _Quantum generalized hydrodynamics_ ,  Phys. Rev. Lett. 124, 140603 (2020),  10.1103/PhysRevLett.124.140603.
  * [57] P. Ruggiero, P. Calabrese, B. Doyon and J. Dubail,  _Quantum generalized hydrodynamics of the tonks–girardeau gas: density fluctuations and entanglement entropy_ ,  Journal of Physics A: Mathematical and Theoretical 55(2), 024003 (2021),  10.1088/1751-8121/ac3d68.
  * [58] L. Landau,  _Theory of the superfluidity of helium ii_ ,  Phys. Rev. 60, 356 (1941),  10.1103/PhysRev.60.356.
  * [59] J. D. Nardis, B. Doyon, M. Medenjak and M. Panfil,  _Correlation functions and transport coefficients in generalised hydrodynamics_ ,  Journal of Statistical Mechanics: Theory and Experiment 2022(1), 014002 (2022),  10.1088/1742-5468/ac3658.
  * [60] B. Doyon, G. Perfetto, T. Sasamoto and T. Yoshimura,  _Emergence of hydrodynamic spatial long-range correlations in nonequilibrium many-body systems_ ,  10.48550/ARXIV.2210.10009 (2022).
  * [61] M. Fagotti,  _Charges and currents in quantum spin chains: late-time dynamics and spontaneous currents_ ,  Journal of Physics A: Mathematical and Theoretical 50(3), 034005 (2016),  10.1088/1751-8121/50/3/034005.
  * [62] M. Fagotti,  _On conservation laws, relaxation and pre-relaxation after a quantum quench_ ,  Journal of Statistical Mechanics: Theory and Experiment 2014(3), P03016 (2014),  10.1088/1742-5468/2014/03/P03016.
  * [63] M. Gaudin,  _Une démonstration simplifiée du théorème de wick en mécanique statistique_ ,  Nuclear Physics 15, 89 (1960),  10.1016/0029-5582(60)90285-6.
  * [64] D. Karevski,  _Scaling behaviour of the relaxation in quantum chains_ ,  The European Physical Journal B - Condensed Matter and Complex Systems 27(1), 147 (2002),  10.1140/epjb/e20020139.
  * [65] T. Platini and D. Karevski,  _Scaling and front dynamics in ising quantum chains_ ,  The European Physical Journal B - Condensed Matter and Complex Systems 48(2), 225 (2005),  10.1140/epjb/e2005-00402-2.
  * [66] A. De Luca, J. Viti, D. Bernard and B. Doyon,  _Nonequilibrium thermal transport in the quantum ising chain_ ,  Phys. Rev. B 88, 134301 (2013),  10.1103/PhysRevB.88.134301.
  * [67] V. Eisler, F. Maislinger and H. G. Evertz,  _Universal front propagation in the quantum Ising chain with domain-wall initial states_ ,  SciPost Phys. 1, 014 (2016),  10.21468/SciPostPhys.1.2.014.
  * [68] B. Bertini, M. Fagotti, L. Piroli and P. Calabrese,  _Entanglement evolution and generalised hydrodynamics: noninteracting systems_ ,  Journal of Physics A: Mathematical and Theoretical 51(39), 39LT01 (2018),  10.1088/1751-8121/aad82e.
  * [69] V. Alba and M. Fagotti,  _Prethermalization at low temperature: The scent of long-range order_ ,  Phys. Rev. Lett. 119, 010601 (2017),  10.1103/PhysRevLett.119.010601.
  * [70] S. Bocini and J.-M. Stéphan,  _Non-probabilistic fermionic limit shapes_ ,  Journal of Statistical Mechanics: Theory and Experiment 2021(1), 013204 (2021),  10.1088/1742-5468/abcd34.
  * [71] P. Calabrese, F. H. L. Essler and M. Fagotti,  _Quantum quenches in the transverse field ising chain: II. stationary state properties_ ,  Journal of Statistical Mechanics: Theory and Experiment 2012(07), P07022 (2012),  10.1088/1742-5468/2012/07/p07022.
  * [72] V. Zauner, M. Ganahl, H. G. Evertz and T. Nishino,  _Time evolution within a comoving window: scaling of signal fronts and magnetization plateaus after a local quench in quantum spin chains_ ,  Journal of Physics: Condensed Matter 27(42), 425602 (2015),  10.1088/0953-8984/27/42/425602.
  * [73] V. Eisler and F. Maislinger,  _Front dynamics in the XY chain after local excitations_ ,  SciPost Phys. 8, 37 (2020),  10.21468/SciPostPhys.8.3.037.


# Statistical Analysis on Random Quantum Sampling by Sycamore and Zuchongzhi
Quantum Processors

Sangchul Oh<EMAIL_ADDRESS>Department of Chemistry, Department of
Physics and Astronomy, and Purdue Quantum Science and Engineering Institute,
Purdue University, West Lafayette, IN, USA  Sabre Kais  Corresponding author:
<EMAIL_ADDRESS>Department of Chemistry, Department of Physics and Astronomy,
and Purdue Quantum Science and Engineering Institute, Purdue University, West
Lafayette, IN, USA

###### Abstract

Random quantum sampling, a task to sample bit-strings from a random quantum
circuit, is considered one of suitable benchmark tasks to demonstrate the
outperformance of quantum computers even with noisy qubits. Recently, random
quantum sampling was performed on the Sycamore quantum processor with 53
qubits [Nature 574, 505 (2019)] and on the Zuchongzhi quantum processor with
56 qubits [Phys. Rev. Lett. 127, 180501 (2021)]. Here, we analyze and compare
statistical properties of the outputs of random quantum sampling by Sycamore
and Zuchongzhi. Using the Marchenko-Pastur law and the Wasssertein distances,
we find that quantum random sampling of Zuchongzhi is more closer to classical
uniform random sampling than those of Sycamore. Some Zuchongzhi’s bit-strings
pass the random number tests while both Sycamore and Zuchongzhi show similar
patterns in heatmaps of bit-strings. It is shown that statistical properties
of both random quantum samples change little as the depth of random quantum
circuits increases. Our findings raise a question about computational
reliability of noisy quantum processors that could produce statistically
different outputs for the same random quantum sampling task.

#### Introduction.—

A quantum computer is believed to simulate quantum systems much better [1] and
to solve some computational tasks exponentially faster [2, 3] than a classical
computer. A demonstration of the outperformance of a quantum computer, called
quantum supremacy [4] or quantum advantage, is considered one of the important
milestones in developing practical quantum computers. Random quantum sampling
[5] is regarded as a good candidate for demonstrating quantum advantage with
noisy intermediate-scale quantum (NISQ) computers available these days.
Recently, quantum advantage has been claimed for random quantum sampling on
the Sycamore quantum processor with 53 superconducting qubits [6] and the
Zuchongzhi quantum processor with 56 superconducting qubits [7], and for Boson
sampling with optical qubits [5, 8, 9].

Quantum advantages of the Sycamore and Zuchongzhi quantum processors over
classical computers for random quantum sampling were verified using the linear
cross-entropy benchmarking (XEB) fidelity, whose values were estimated
slightly larger than zero [6, 7]. Both quantum processors are made of two-
dimensional arrays of superconducting transmon qubits, have similar error
rates, executed the same random quantum circuit, and scored similar XEB
values. So, one may speculate that the output bit-strings generated by the two
noisy quantum processors would be statistically close. However, it is unknown
whether two noisy quantum processors with similar values of the XEB fidelity
sample statistically similar bit-strings or not. In this paper, we analyze the
statistical closeness of the outputs of the two noisy quantum processors using
the NIST random number tests [10], the Marchenko-Pastur distribution of
eigenvalues of random matrices of bit-strings [11], and the Wasserstein
distance between samples [12]. We show that Sycamore’s random bit-strings are
farther away from classical uniform random bit-strings than Zuchongzhi’s
outputs, that is, the two outputs of bit-strings are noticeably different.

#### Random quantum sampling.—

Let us first summarize random quantum sampling implemented on the Sycamore and
Zuchongzhi quantum processors [6, 7]. The task of random quantum sampling is
to sample bit-strings $x=a_{1}\cdots a_{n}\in\\{0,1\\}^{n}$ by applying a
random quantum circuit $U$ on the initial state $\ket{0}$ of $n$ qubits
followed by the measurement. Both Sycamore and Zuchongzhi quantum processors
executed the same random quantum circuit $U$ composed of $m$ cycles as
follows. A quantum circuit $U_{k}$ at the $k$-th cycle consists of single-
qubit gates $R$ selected randomly from the set
$\\{\sqrt{X},\sqrt{Y},\sqrt{W}\\}$ on all qubits and deterministic two-qubit
gates on the pair of qubits selected in the sequence of the coupler activation
patterns of 2-dimensional superconducting qubits. After $m$ cycles, a final
single-qubit gate $R$ is applied before the measurement, so the whole random
quantum circuit with cycle $m$ is given by $U=RU_{m}U_{m-1}\cdots U_{1}$. The
output of the random quantum circuit is a bit-string $x$ that is sampled from
the probability $p(x)=|\bra{x}U\ket{0}|^{2}$. By implementing the same random
quantum circuit $M$ times, a collection of $M$ random bit-strings ${\cal
D}=\\{x_{1},\dots,x_{M}\\}$, i.e., an $M\times n$ binary bit array, is
obtained.

Figure 1:  Heatmaps of quantum random sampling (a) by the Sycamore quantum
processor with $n=53$ qubits, $m=20$ cycles, and sample numbers $M=10^{6}$,
(b) by the Zuchongzhi quantum processor with $n=56$ qubits, $m=20$ cycles, and
$M=10^{6}$ samples, and (c) of classical uniform random sampling with $n=56$
and $M=10^{6}$ are plotted. The average of finding bit 1 is denoted by
$p_{1}$. Random quantum samples, (a) and (b), show the stripe patterns at
specific qubit indices while (b) and (c) pass the NIST random number tests.

To verify that a noisy quantum processor is performing well random quantum
sampling [6, 7], the linear cross-entropy benchmarking (XEB) fidelity $F_{\rm
XEB}$ is introduced

$F_{\rm XEB}=2^{n}\cdot\frac{1}{M}\sum_{x\in\cal D}p(x)-1\,,$ (1)

where the ideal probability $p(x)=|\bra{x}U\ket{0}|^{2}$ of finding a bit-
string $x$ is computed on a classical computer using Schrödinger or Feynman
simulators and the bit-strings ${\cal D}=\\{x_{1},\dots,x_{M}\\}$ are
generated by a quantum processor. It is known that $F_{\rm XEB}=1$ if a
quantum processor implements a random quantum circuit without errors, and
$F_{\rm XEB}=0$ if bit-strings are sampled from a classical uniform
distribution. The estimated XEB fidelity of Sycamore is $F_{\rm XEB}=(2.24\pm
0.021)\times 10^{-3}$ for 53 qubits, 20 cycles, and $M=30\times 10^{6}$
samples over 10 circuit instances [6]. For the Zuchongzhi quantum processor
[7], the estimated XEB fidelity is $F_{\rm XEB}=(6.62\pm 0.72)\times 10^{-4}$
for 56 qubits, 20 cycles, and $M=1.9\times 10^{7}$ samples.

While the XEB fidelity could serve as a benchmark for quantum advantage of
random quantum sampling, it has some limitations. The XEB fidelity stems from
the Kullback-Leibler divergence or the cross entropy of an empirical
probability distribution $\tilde{p}(x)$ of the data ${\cal D}$ from the ideal
probability distribution $p(x)$ [13]. As the number of qubits increases, it is
very difficult to calculate both probability distributions. The number of
samples required to construct the empirical probability $\tilde{p}(x)$
increases exponentially because the range of $x\in[0,2^{n}-1]$ does so. While
the XEB fidelity needs only an ideal probability $p(x)$ but not
$\tilde{p}(x)$, it is still not scalable because in quantum advantage regime a
supercomputer cannot calculate the ideal probability $p(x)$ in Eq. (1) and no
other NISQ processors could do. Recently, the limitation of the XEB fidelity
as a benchmark for quantum advantage have been pointed out by Gao et. al.
[14]. More importantly, the XEB fidelity could not give any clue about
statistical properties of bit-strings of random quantum sampling implemented
on NISQ processors.

#### Comparison of random quantum sampling of Sycamore and Zuchongzhi.—

A simple way of comparing the performance of two quantum processors is to
compare directly their outputs of a task, here random quantum sampling. Both
Sycamore and Zuchongzhi have similar noise levels. The averages of single-
qubit gate errors, two-qubit gates errors, and the readout error of Sycamore
is about 0.15%, 0.36%, and 3.1%, respectively [6]. The average errors of
Zuchongzhi are about 0.14% for single-qubit gates, 0.59% for two-qubit gates,
and 4.52% for readout, respectively [7]. Both processors implemented the same
random quantum circuit as described before and obtained similar values of XEB
fidelity. So, one may expect that their outputs would be statistically close
each other. However, that is not the case.

Figure 2: (a) For $n=51$ qubits, the outliers of the Marchenko-Pastur
distributions of eigenvalues of classical random bit-strings (green),
Sycamore’s quantum random bit-strings (red), and Zuchongzhi’s quantum random
bit-strings (blue) are plotted. The inset of (a) depicts the bulk distribution
of eigenvalues of the Marchenko-Pastur law given by Eq. (3). (b) The distances
of the outlier peaks of Sycamore (red solid circle ), Zuchongzhi (blue solid
square), and classical random sampling (green triangle) from $n/4$ are plotted
as a function of the qubit number $n$. (c) The distances of the peaks of
Sycamore with $n=53$ (solid red circle) and Zuchongzhi with $n=56$ (blue solid
square) from $n/4$ is plotted as a function of cycle $m$. The dot-dash lines
in (b) and (c) represent the average distance of classical random sampling
from $n/4$.

We first examine whether the output bit-strings of random quantum sampling are
random. Random numbers have many practical applications in Monte-Carlo
simulation, statistics, cryptography, etc.. Given a random quantum circuit
$U$, the probability $p(x)=|\bra{x}U\ket{0}|^{2}$ of finding bit-string $x$ is
not uniform due to the entanglement and interference, so some bit-strings are
sampled more likely than others. However, there is no physical ground that the
output of bit-strings ${\cal D}=\\{x_{1},x_{2},\dots,x_{M}\\}$ would contain
more bit 1 than 0 or vice versa. Bits $0$ and $1$ are equally probable if
there is no error. To see this, the output of bit-strings ${\cal D}$ is sliced
into the collection of $n\times n$ binary arrays, with which the heatmaps are
plotted as shown in Fig. 1. Fig. 1 (a) depicts the heat map of Sycamore’s bit-
strings for $n=53$ qubits, $m=20$ cycles, and $M=3\times 10^{6}$ samples [15].
Fig. 1 (b) displays the heatmap of Zuchongzhi’s bit-strings for $n=56$ qubits,
$m=18$ cycles, and $M=3\times 10^{6}$ samples. Fig. 1 (c) shows the heatmap of
uniform random bits sampled from a classical computer for $n=56$ and
$M=3\times 10^{6}$. As depicted in Fig. 1, the heatmaps of random quantum
sampling on Sycamore and Zuchongzhi quantum processors show bright and darks
stripes at some qubit indices while the classical uniform random sampling does
not. One may suspect that stripe patterns could be caused by readout errors.
However, for both Sycamore and Zuchongzhi data sets, the locations of bright
or dark stripes in the heatmaps do not coincide with the indices of qubits
with high readout errors.

We count the number of bit $1$ contained in the output bit-strings ${\cal
D}=\\{x_{1},\cdots,x_{M}\\}$. The average of finding bit 1 is denoted by
$p_{1}$. Most Zuchongzhi outputs show $p_{1}$ is greater than 1/2, as shown in
Supplementary Material [16], while all Sycamore outputs have $p_{1}$ less than
1/2 [17]. As depicted in Fig. 1, Zuchongzhi’s bit-strings with $n=56$ qubits,
$m=18$ cycles has $p_{1}=0.50094$ while Sycamore’s bit-strings for $n=53$
qubits and $m=20$ cycles has $p_{1}=0.48383$. Note that the readout errors of
Sycamore and Zuchongzhi are reported to be 3.1% and 4.52%, respectively. It is
unclear why Sycamore and Zuchongzhi are quite different in the probability
$p_{1}$ of finding bit $1$ in random quantum sampling.

To see the randomness of bit-strings, we perform the NIST random number tests
[10] for Zuchongzhi’s outputs. We find that some Zuchongzhi data sets pass the
NIST random number tests while all Sycamore data sets do not [16, 17]. Note
that the probability $p_{1}$ of finding bit 1 is related to the frequency test
of the NIST random number tests. It is interesting to see that some Zuchongzhi
data sets pass the NIST random number tests while all random quantum sampling
data of Sycamore and Zuchongzhi show the stripe patterns unlike the classical
uniform random samples, as shown in Fig. 1. This implies random quantum
sampling may be used to generate quantum random numbers [18, 19].

Both Sycamore and Zuchongzhi showed similar exponential decays of the XEB
fidelity $F_{\rm XEB}$ as a function of qubit number $n$ or cycle $m$ [6, 7].
One obtains $F_{\rm XEB}=0$ if bit-strings are sampled from the classical
uniform random distribution $p_{\rm cl}(x)=1/2^{n}$. Using the Marchenko-
Pastur distribution [11] and the Wasserstein distances [12, 20], we examine
how the random bit-string outputs of Sycamore and Zuchongzhi are far away from
classical uniform random samples as a function of qubit number $n$ and cycle
$m$.

The signal out of randomness can be captured by the outliers of the Marchenko-
Pastur distribution of random matrices. This method has been applied to the
covariance matrices in finance [21, 22] in predicting ligand affinity [23],
and in denoising MRI [24] and single-cell data [25]. We use the Marchenko-
Pastur distribution of eigenvalues of random bit-strings to measure the
distance among Sycamore’s outputs, Zuchongzhi’s output, and classical uniform
random samples. The data ${\cal D}=\\{x_{1},x_{2},\dots,x_{M}\\}$, an binary
$M\times n$ array, is sliced into the collection of $k\times n$ random
matrices $X$. All entries of $X$ are assumed to be independent and identically
binary random variable $\\{0,1\\}$ with the mean $\mu_{X}=1/2$ and the
variance $\sigma_{X}^{2}=1/4$. We calculate the empirical distribution of
eigenvalues of $n\times n$ matrix $\frac{1}{k}X^{t}\cdot X$, as shown in Fig.
2. The empirical eigenvalue distribution is composed of the two parts: the
bulk corresponding to the noise or random and the outliers representing the
signal. To understand this, let consider the transformed random matrix
$Y=2X-J$ with $J$ of all entries 1. The matrix $\frac{1}{k}X^{t}\cdot X$ is
written as

$\displaystyle\frac{1}{k}X^{t}\cdot X=\frac{1}{4k}\left(Y^{t}\cdot
Y+X^{t}\cdot J+X^{t}\cdot Z+J^{t}\cdot J\right).$ (2)

The matrix $Y$ has the mean $\mu_{Y}=0$ and the variance $\sigma^{2}_{Y}=1$.
The eigenvalue distribution of the first term in Eq. (2),
$\frac{1}{k}Y^{t}\cdot Y$, follows the Marchenko-Pastur distribution [11]

$\displaystyle\rho(\lambda)=\frac{1}{2\pi\sigma^{2}\gamma}\frac{\sqrt{(\lambda_{+}-\lambda)(\lambda-\lambda_{-})}}{\lambda}\,,$
(3)

where $\gamma=n/k$ is the rectangular ratio and
$\lambda_{\pm}=\sigma^{2}(1\pm\sqrt{\gamma})^{2}$ are the upper and lower
bounds. Here we take $k=2n$, i.e., $\gamma=1/2$. The upper limit is given by
$\lambda_{+}=1+\sqrt{1/2}$. By considering the scaling factor $1/4$, the upper
limit of $(1/4k)\cdot Y^{t}\cdot Y$ is given by $(1+\sqrt{1/2})/4\approx
0.7285$. The last term of Eq. (2), $(1/4k)J^{t}\cdot J$ has the two
eigenvalues, 0 and $n/4$. So the outliers are located around $n/4$.

Fig. 2 (a) plots the empirical distribution of eigenvalues of
$\frac{1}{k}X^{t}\cdot X$ made of Sycamore bit-strings with $n=51$, $m=14$,
and $M=10^{6}$ (red color), Zuchongzhi bit-strings with $n=51$, $m=10$ and
$M=10^{6}$ (blue color), and classical uniform random bits with $M=10^{6}$
(green color). The outlier peak of Sycamore is located at the left of $n/4$
while the outlier peaks of Zuchongzhi and classical random sampling are
located at right of $n/4$. Surprisingly, Zuchongzhi’s bit-strings are more
closer to the classical random bit strings than Sycamore’s ones while Sycamore
and Zuchongzhi have similar the values of the XEB fidelity. Fig. 2 (b) plots
the distances of the outlier peaks from $n/4$ as a function of qubit number
$n$. Sycamore becomes farther away from $n/4$ as $n$ increases. On the other
hands, Zuchongzhi’s distance from $n/4$ changes a little. The XEB fidelity of
Sycamore and Zuchongzhi decays exponentially as a function of qubit number
$n$. Fig. 2 (c) plots the peak distance from $n/4$ as a function of cycles for
Sycamore $n=53$ and for Zuchongzhi $n=56$. This result is also in contrast
with the decay behavior of the XEB fidelity as a function of cycle $m$.

The statistical distances among Sycamore, Zuchongzhi, and classical random
sampling, based on the Marchenko-Pastur distribution, can be further confirmed
by calculating the Wasserstein distances. The 1-th Wasserstein distance
between two probability distribution $p(x)$ and $q(x)$ [12] is defined by

$\displaystyle
W_{p}(p,q)=\inf_{\pi\in\Pi(p,q)}{\mathbb{E}}_{(x,y)\sim\pi}\left[||x-y||\right]\,,$
(4)

where $\Pi(p,q)$ is the set of all joint distribution $\pi(x,y)$ whose
marginal distributions are $p(x)$ and $p(y)$, respectively. Given two samples,
$\\{x_{1},x_{2},\dots,x_{M}\\}$ and $\\{y_{1},y_{2},\dots,y_{M}\\}$, $W(p,q)$
can be calculated directly without calculating the empirical distributions
$p(x)$ and $q(x)$ [20]. The Wasserstein distance is a true metric on a
probability space, so the relative Wasserstein distances among Sycamore,
Zuchongzhi, and classical random sampling, give rise to the triangle
inequality. Fig. 3 (a) plot the triangle relation among the three data sets
for $n=18,24,30,36,42,45,48,51$. All Zuchongzhi data sets are close to
classical random sampling than Sycamore. This is consistent with the outlier
peak distances from $n/4$ of the Marchenko-Pastur distribution shown in Fig. 2
(b). Fig. 3 (b) plots the Wasserstein distance of Sycamore $n=53$ and
Zuchongzhi $n=56$ from the classical random bit-strings as a function of cycle
$m$. This result is consistent with the behavior of the relative distances of
the outliers of the Marchenko-Paster distribution as a function of $m$, shown
in Fig. 2 (c). Note that the Wasserstein distances between the bit-strings
with different cycles of Sycamore (or Zuchongzhi) are less than the
Wasserstein distance of them from the classical random sampling.

Figure 3: (a) For qubits $n=18,24,30,36,42,45,48,51$, the relative positions
of Sycamore (red circle), Zuchongzhi (blue square), and classical random
sampling (black circle at origin) form triangles whose edges are the
Wasserstein distances. For clarity, triangles are rotated. (b) for Sycamore
bit-strings with $n=51$ and Zuchongzhi with $n=56$, the Wasserstein distances
from the classical random sampling are plotted as a function of cycle $m$.

#### Conclusion.—

In this paper, we analyzed the statistical properties of random bit-strings
sampled from Sycamore and Zuchongzhi quantum processors using the heatmaps,
the NIST random number tests, the Marchenko-Pastur distribution, and the
Wasserstein distances. Both Sycamore and Zuchongzhi exhibit stripe patterns in
the heatmaps of random bit-strings. Some of Zuchongzhi’s data pass the NIST
random number tests. This may open up a possibility of the use of the random
quantum circuit as a quantum random number generator. The Marchenko-Pastur
distribution and the Wasserstein distances of random bit-strings shows that
Zuchongzhi random quantum sampling is statistically more closer to classical
random sampling than Sycamore’s, while both have similar error rates and
similar values of the XEB fidelity. The distances of random bit-strings of
Sycamore from the classical uniform sampling as a function of the number of
qubit $n$ is quite different from those of Zuchongzhi. The distances of
Sycamore’s and Zuchongzhi’s bit-strings from the classical random bit-strings
remain almost constant as the cycle number $m$ increases. This is in contrast
with the exponential decay of the XEB fidelity as a function of $m$.

The results found here raise a question about computational reliability and
verification of noisy quantum processors [26]. Two classical computers
performing the same task will produce the same outputs (within statistical
errors for stochastic calculation). One may expect that two NISQ processors
with similar error rates would produce statistically close outputs. However,
the results here show an exceptional case. In order to ensure the reliability
of quantum computing with NISQ quantum processors, it is necessary to have
good benchmark tools which can verify the output of quantum calculation on a
classical computer, until fully-error corrected quantum computers are
available.

###### Acknowledgements.

We would like to thank the Google quantum team for making their random quantum
sampling by Sycamore available. We also would like to thank Y. Wu and J.-W.
Pan for providing us the data on Zuchongzhi [7]. This material is based upon
work supported by the U.S. Department of Energy, Office of Science, National
Quantum Information Science Research Centers. We also acknowledge the National
Science Foundation under award number 1955907.

## References

  * Feynman [1982] R. P. Feynman, International Journal of Theoretical Physics 21, 467 (1982).
  * Shor [1997] P. W. Shor, SIAM Journal on Computing 26, 1484 (1997).
  * Harrow  _et al._ [2009] A. W. Harrow, A. Hassidim, and S. Lloyd, Phys. Rev. Lett. 103, 150502 (2009).
  * Preskill [2012] J. Preskill, Quantum computing and the entanglement frontier (2012), arXiv:1203.5813 [quant-ph] .
  * Aaronson and Arkhipov [2013] S. Aaronson and A. Arkhipov, Theory of Computing 9, 143 (2013).
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao, D. A. Buell, B. Burkett, Y. Chen, Z. Chen, B. Chiaro, R. Collins, W. Courtney, A. Dunsworth, E. Farhi, B. Foxen, A. Fowler, C. Gidney, M. Giustina, R. Graff, K. Guerin, S. Habegger, M. P. Harrigan, M. J. Hartmann, A. Ho, M. Hoffmann, T. Huang, T. S. Humble, S. V. Isakov, E. Jeffrey, Z. Jiang, D. Kafri, K. Kechedzhi, J. Kelly, P. V. Klimov, S. Knysh, A. Korotkov, F. Kostritsa, D. Landhuis, M. Lindmark, E. Lucero, D. Lyakh, S. Mandrà, J. R. McClean, M. McEwen, A. Megrant, X. Mi, K. Michielsen, M. Mohseni, J. Mutus, O. Naaman, M. Neeley, C. Neill, M. Y. Niu, E. Ostby, A. Petukhov, J. C. Platt, C. Quintana, E. G. Rieffel, P. Roushan, N. C. Rubin, D. Sank, K. J. Satzinger, V. Smelyanskiy, K. J. Sung, M. D. Trevithick, A. Vainsencher, B. Villalonga, T. White, Z. J. Yao, P. Yeh, A. Zalcman, H. Neven, and J. M. Martinis, Nature 574, 505 (2019).
  * Wu  _et al._ [2021] Y. Wu, W.-S. Bao, S. Cao, F. Chen, M.-C. Chen, X. Chen, T.-H. Chung, H. Deng, Y. Du, D. Fan, M. Gong, C. Guo, C. Guo, S. Guo, L. Han, L. Hong, H.-L. Huang, Y.-H. Huo, L. Li, N. Li, S. Li, Y. Li, F. Liang, C. Lin, J. Lin, H. Qian, D. Qiao, H. Rong, H. Su, L. Sun, L. Wang, S. Wang, D. Wu, Y. Xu, K. Yan, W. Yang, Y. Yang, Y. Ye, J. Yin, C. Ying, J. Yu, C. Zha, C. Zhang, H. Zhang, K. Zhang, Y. Zhang, H. Zhao, Y. Zhao, L. Zhou, Q. Zhu, C.-Y. Lu, C.-Z. Peng, X. Zhu, and J.-W. Pan, Phys. Rev. Lett. 127, 180501 (2021).
  * Zhong  _et al._ [2020] H.-S. Zhong, H. Wang, Y.-H. Deng, M.-C. Chen, L.-C. Peng, Y.-H. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, P. Hu, X.-Y. Yang, W.-J. Zhang, H. Li, Y. Li, X. Jiang, L. Gan, G. Yang, L. You, Z. Wang, L. Li, N.-L. Liu, C.-Y. Lu, and J.-W. Pan, Science 370, 1460 (2020).
  * Zhong  _et al._ [2021] H.-S. Zhong, Y.-H. Deng, J. Qin, H. Wang, M.-C. Chen, L.-C. Peng, Y.-H. Luo, D. Wu, S.-Q. Gong, H. Su, Y. Hu, P. Hu, X.-Y. Yang, W.-J. Zhang, H. Li, Y. Li, X. Jiang, L. Gan, G. Yang, L. You, Z. Wang, L. Li, N.-L. Liu, J. J. Renema, C.-Y. Lu, and J.-W. Pan, Phys. Rev. Lett. 127, 180502 (2021).
  * Bassham  _et al._ [2010] L. Bassham, A. Rukhin, J. Soto, J. Nechvatal, M. Smid, S. Leigh, M. Levenson, M. Vangel, N. Heckert, and D. Banks, A statistical test suite for random and pseudorandom number generators for cryptographic applications (2010).
  * Marchenko and Pastur [1967] V. A. Marchenko and L. A. Pastur, Mat. Sb. (N.S.) 72, 507 (1967).
  * Villani [2008] C. Villani,  _Optimal Transport: Old and New_, Grundlehren der mathematischen Wissenschaften (Springer Berlin Heidelberg, 2008).
  * Boixo  _et al._ [2018] S. Boixo, S. V. Isakov, V. N. Smelyanskiy, R. Babbush, N. Ding, Z. Jiang, M. J. Bremner, J. M. Martinis, and H. Neven, Nature Physics 14, 595 (2018).
  * Gao  _et al._ [2021] X. Gao, M. Kalinowski, C.-N. Chou, M. D. Lukin, B. Barak, and S. Choi, Limitations of linear cross-entropy as a measure for quantum advantage (2021), arXiv:2112.01657 [quant-ph] .
  * Martinis et al. [2021] J. M. Martinis et al., Quantum supremacy using a programmable superconducting processor, dryad, dataset (2021).
  * [16] See Supplemental Material for the details of calculation.
  * Oh and Kais [2021] S. Oh and S. Kais, Non-randomness of google’s quantum supremacy benchmark (2021), arXiv:2110.06046 [quant-ph] .
  * Herrero-Collantes and Garcia-Escartin [2017] M. Herrero-Collantes and J. C. Garcia-Escartin, Rev. Mod. Phys. 89, 015004 (2017).
  * Tamura and Shikano [2021] K. Tamura and Y. Shikano, in  _International Symposium on Mathematics, Quantum Theory, and Cryptography_ , edited by T. Takagi, M. Wakayama, K. Tanaka, N. Kunihiro, K. Kimoto, and Y. Ikematsu (Springer Singapore, Singapore, 2021) pp. 17–37.
  * Flamary  _et al._ [2021] R. Flamary, N. Courty, A. Gramfort, M. Z. Alaya, A. Boisbunon, S. Chambon, L. Chapel, A. Corenflos, K. Fatras, N. Fournier, L. Gautheron, N. T. Gayraud, H. Janati, A. Rakotomamonjy, I. Redko, A. Rolet, A. Schutz, V. Seguy, D. J. Sutherland, R. Tavenard, A. Tong, and T. Vayer, Journal of Machine Learning Research 22, 1 (2021).
  * Laloux  _et al._ [1999] L. Laloux, P. Cizeau, J.-P. Bouchaud, and M. Potters, Phys. Rev. Lett. 83, 1467 (1999).
  * Laloux  _et al._ [2000] L. Laloux, P. Cizeau, M. Potters, and J.-P. Bouchaud, International Journal of Theoretical and Applied Finance 03, 391 (2000).
  * Lee  _et al._ [2016] A. A. Lee, M. P. Brenner, and L. J. Colwell, Proceedings of the National Academy of Sciences 113, 13564 (2016).
  * Veraart  _et al._ [2016] J. Veraart, E. Fieremans, and D. S. Novikov, Magnetic Resonance in Medicine 76, 1582 (2016), https://onlinelibrary.wiley.com/doi/pdf/10.1002/mrm.26059 .
  * Aparicio  _et al._ [2020] L. Aparicio, M. Bordyuh, A. J. Blumberg, and R. Rabadan, Patterns 1, 100035 (2020).
  * Eisert  _et al._ [2020] J. Eisert, D. Hangleiter, N. Walk, I. Roth, D. Markham, R. Parekh, U. Chabaud, and E. Kashefi, Nature Reviews Physics 2, 382 (2020).

## Supplementary material:
Statistical Analysis on Random Quantum Sampling by Sycamore and Zuchongzhi
Quantum Processors

### Sangchul Oh and Sabre Kais

Supplemental material presents the heatmaps and the NIST random number tests
for random bit-strings of Zuchongzhi. The Wasserstein distances among
different cycles for Sycamore and Zuchongzhi are presented.

Figure 4: Heatmaps of random bit-strings of Zuchongzhi with $n=56$ and
$m=10,12,14,16,18,20$.

Figure 5: Heatmaps of random bit-strings of Zuchongzhi with
$n=15,18,21,...,51$ and $m=10$.

File name |  $M$xxx |  $p_{1}$ |  NIST random number tests
---|---|---|---
`15_10_F.csv` | 1000000 |  0.5042914666666667 |  Nonrandom
`18_10_F.csv` | 1000000 |  0.5039382222222222 |  Nonrandom
`21_10_F.csv` | 1000000 |  0.5016010476190477 |  Nonrandom
`24_10_F.csv` | 1000000 |  0.502242875 |  Random
`27_10_F.csv` | 1000000 |  0.5015492592592593 |  Random
`30_10_F.csv` | 1000000 |  0.5014525333333333 |  Random
`33_10_F.csv` | 1000000 |  0.5014795757575757 |  Nonrandom
`36_10_F.csv` | 1000000 |  0.5012370277777778 |  Nonrandom
`39_10_F.csv` | 1000000 |  0.5014506153846154 |  Random
`42_10_F.csv` | 1000000 |  0.5007865238095238 |  Nonrandom
`45_10_F.csv` | 1000000 |  0.5006884222222222 |  Random
`48_10_F.csv` | 1000000 |  0.5009762291666666 |  Random
`51_10_F.csv` | 1000000 |  0.5008632549019608 |  Nonrandom
`54_10_F.csv` | 1050000 |  0.5009445502645503 |  Nonrandom
`56_10_F.csv` | 1170000 |  0.5007047619047619 |  Random
Table 1: NIST random number tests for random bit-strings of Zuchongzhi with
$n=15,...,56$ and $m=10$. $M$ stands for the number of bit-strings and $p_{1}$
for the probability of finding bit 1. 5 data files in directories, a1, a2, a3,
a4, a5, are merged into a single file. Nonrandom means that bit-strings fail
some of the NIST random tests.

    Test Data File:/home/soh/Desktop/data/Zuchongzhi_nxx_m10/56_10_F_1M.csv

    Type of Test                                        P-Value                 Conclusion
    01. Frequency Test (Monobit)                        0.04725583993159622     Random
    02. Frequency Test within a Block                   0.7377345740502045      Random
    03. Run Test                                        0.9013653598338198      Random
    04. Longest Run of Ones in a Block                  0.5576699582990581      Random
    05. Binary Matrix Rank Test                         0.13573737861315063     Random
    06. Discrete Fourier Transform (Spectral) Test      0.4912971242158931      Random
    07. Non-Overlapping Template Matching Test          0.7861450901488796      Random
    08. Overlapping Template Matching Test              0.3562201622350592      Random
    09. Maurer’s Universal Statistical test             0.18901932905865293     Random
    10. Linear Complexity Test                          0.6012867491687948      Random
    11. Serial test:
                                                        0.9287873362845678      Random
                                                        0.6473140941950127      Random
    12. Approximate Entropy Test                        0.9658862333540649      Random
    13. Cummulative Sums (Forward) Test                 0.08552956685396966     Random
    14. Cummulative Sums (Reverse) Test                 0.09230410462727737     Random
    15. Random Excursions Test:
                    State       Chi Squared             P-Value                 Conclusion
                    -4          2.2829376648618633      0.8087695808555937      Random
                    -3          2.9604                  0.706091109222725       Random
                    -2          5.783950617283951       0.3278123367522591      Random
                    -1          6.75                    0.23990688227010742     Random
                    +1          2.916666666666667       0.7128316836671243      Random
                    +2          1.876543209876543       0.8659503247833513      Random
                    +3          2.6755999999999998      0.7498477047163237      Random
                    +4          4.022629459947245       0.5461624839081891      Random
    16. Random Excursions Variant Test:
                    State       COUNTS                  P-Value                 Conclusion
                    -9.0        6                       0.5286121252556878      Random
                    -8.0        9                       0.5761501220305789      Random
                    -7.0        17                      0.7793054481711047      Random
                    -6.0        21                      0.8961247779728454      Random
                    -5.0        12                      0.5637028616507731      Random
                    -4.0        11                      0.47819532084058625     Random
                    -3.0        14                      0.5186050164287256      Random
                    -2.0        17                      0.5596689271994115      Random
                    -1.0        19                      0.47048642205878966     Random
                    +1.0        29                      0.47048642205878966     Random
                    +2.0        33                      0.4532547047537364      Random
                    +3.0        49                      0.1065831695748876      Random
                    +4.0        56                      0.08085559837005228     Random
                    +5.0        60                      0.08326451666355045     Random
                    +6.0        64                      0.08172275229865938     Random
                    +7.0        63                      0.11846489387154538     Random
                    +8.0        58                      0.20511768103340522     Random
                    +9.0        39                      0.5995101785600032      Random

Table 2: Result of the NIST random number tests for Zuchongzhi’s data 56-10-F.csv. File name |  $M$xxx |  $p_{1}$ |  NIST random number tests
---|---|---|---
`56_12_F_s11.csv` | 1000000 |  0.5004180892857143 |  Nonrandom
`56_12_F_s12.csv` | 1000000 |  0.4997393035714286 |  Random
`56_12_F_s13.csv` | 1000000 |  0.5016109285714285 |  Nonrandom
`56_12_F_s15.csv` | 1000000 |  0.5005064107142857 |  Random
`56_12_F_s3.csv` | 1000000 |  0.5004619464285714 |  Random
`56_12_F_s4.csv` | 1000000 |  0.5015308571428572 |  Nonrandom
`56_12_F_s5.csv` | 1000000 |  0.5012426964285714 |  Random
`56_12_F_s6.csv` | 1000000 |  0.5019179107142857 |  Random
`56_12_F_s8.csv` |  980000 |  0.5011627733236151 |  Random
`56_12_F_s9.csv` | 1000000 |  0.5012499464285715 |  Nonrandom
`56_14_F_s11.csv` | 1370000 |  0.502104027632951 |  Nonrandom
`56_14_F_s12.csv` | 1440000 |  0.5003947048611112 |  Nonrandom
`56_14_F_s13.csv` | 1700000 |  0.5007865021008403 |  Nonrandom
`56_14_F_s15.csv` | 1120000 |  0.5003659757653062 |  Nonrandom
`56_14_F_s3.csv` | 1100000 |  0.5011199512987013 |  Nonrandom
`56_14_F_s4.csv` | 1290000 |  0.5018366832779624 |  Nonrandom
`56_14_F_s5.csv` | 1310000 |  0.5014550572519084 |  Nonrandom
`56_14_F_s6.csv` | 1400000 |  0.5016173852040816 |  Nonrandom
`56_14_F_s8.csv` | 1480000 |  0.5015874396718146 |  Nonrandom
`56_14_F_s9.csv` | 1390000 |  0.5016998843782117 |  Nonrandom
`56_16_F_s11.csv` | 3560000 |  0.501555216693419 |  Nonrandom
`56_16_F_s12.csv` | 3770000 |  0.5006029840848807 |  Random
`56_16_F_s13.csv` | 3680000 |  0.5012373495729814 |  Nonrandom
`56_16_F_s15.csv` | 2870000 |  0.5008621329019413 |  Nonrandom
`56_16_F_s3.csv` | 2890000 |  0.5017474913494809 |  Nonrandom
`56_16_F_s4.csv` | 3200000 |  0.5013628404017857 |  Nonrandom
`56_16_F_s5.csv` | 2920000 |  0.5025822773972602 |  Nonrandom
`56_16_F_s6.csv` | 3570000 |  0.5023631752701081 |  Nonrandom
`56_16_F_s8.csv` | 3710000 |  0.5017546592221794 |  Nonrandom
`56_16_F_s9.csv` | 3420000 |  0.5018155127401838 |  Nonrandom
`56_18_F_s11.csv` | 10080000 |  0.5012828125 |  Nonrandom
`56_18_F_s12.csv` | 9730000 |  0.5006819630010277 |  Random
`56_18_F_s13.csv` | 11530000 |  0.5014127137281625 |  Nonrandom
`56_18_F_s15.csv` | 9530000 |  0.5014926154249738 |  Random
`56_18_F_s3.csv` | 7660000 |  0.5021444307161507 |  Nonrandom
`56_18_F_s4.csv` | 7510000 |  0.5024223606619745 |  Nonrandom
`56_18_F_s5.csv` | 7960000 |  0.5020297963029433 |  Nonrandom
`56_18_F_s6.csv` | 9390000 |  0.5027202723261829 |  Nonrandom
`56_18_F_s8.csv` | 8620000 |  0.501850381173351 |  Nonrandom
`56_18_F_s9.csv` | 7470000 |  0.5023868091413272 |  Nonrandom
`56_20_F_s11.csv` | 27900000 |  0.5018168573988735 |  Nonrandom
`56_20_F_s12.csv` | 25870000 |  0.5013624744602132 |  Nonrandom
`56_20_F_s13.csv` | 26580000 |  0.5017120135977642 |  Nonrandom
`56_20_F_s15.csv` | 20100000 |  0.5013940840440654 |  Nonrandom
`56_20_F_s3.csv` | 19100000 |  0.502459691473448 |  Nonrandom
`56_20_F_s4.csv` | 21510000 |  0.5024990527661553 |  Nonrandom
`56_20_F_s5.csv` | 22570000 |  0.5026209222102664 |  Nonrandom
`56_20_F_s6.csv` | 24510000 |  0.5030497945445008 |  Nonrandom
`56_20_F_s8.csv` | 19840000 |  0.5022362093173963 |  Nonrandom
`56_20_F_s9.csv` | 23510000 |  0.5026191111684998 |  Nonrandom
Table 3: NIST random number tests for random bit-strings of Zuchongzhi with $n=56$ and $m=12,14,16,18,20$. $M$ stands for the number of bit-strings and $p_{1}$ for the probability of finding bit 1. Nonrandom means that bit-strings fail some of the NIST random tests. |  classical n53 |  n53-m12 |  n53-m14 |  n53-m16 |  n53-m18 |  n53-m20
---|---|---|---|---|---|---
classical n53 |  0.0 |  0.02495338 |  0.02459099 |  0.02560033 |  0.02514710 |  0.0249869
n53-m12 |  0.02495337 |  0.0 |  0.00083951 |  0.00040303 |  0.00050475 |  0.0005789
n53-m14 |  0.02459099 |  0.00083951 |  0.0 |  0.00099459 |  0.00128885 |  0.0013602
n53-m16 |  0.02560033 |  0.00040303 |  0.00099459 |  0.0 |  0.00066794 |  0.0007389
n53-m18 |  0.02514710 |  0.00050475 |  0.00128885 |  0.00066794 |  0.0 |  0.0002826
n53-m20 |  0.02498690 |  0.00057892 |  0.00136027 |  0.00073890 |  0.00028265 |  0.0
Table 4: For $n=53$, the Wasserstein distances among classical random bit-strings and Sycamore random bit-strings with cycle $m=12,14,16,18,20$. The Wasserstein distance between the classical random bit-strings and Sycamore’s random bit-strings is much larger than those between Sycamore’s random bit-strings with different cycle. |  classical n56 |  n56-m12 |  n56-m14 |  n56-m16 |  n56-m18 |  n56-m20
---|---|---|---|---|---|---
classical n56 |  0.0 |  0.00173514 |  0.00183101 |  0.00171989 |  0.00161961 |  0.00174550
n56-m12 |  0.00173514 |  0.0 |  0.00071221 |  0.00023922 |  0.00056779 |  0.00034071
n56-m14 |  0.00183101 |  0.00071221 |  0.0 |  0.00085164 |  0.00125725 |  0.00073764
n56-m16 |  0.00171989 |  0.00023922 |  0.00085164 |  0.0 |  0.00045474 |  0.00025549
n56-m18 |  0.00161961 |  0.00056779 |  0.00125725 |  0.00045474 |  0.0 |  0.00056125
n56-m20 |  0.00174550 |  0.00034071 |  0.00073764 |  0.00025549 |  0.00056125 |  0.0
Table 5: For $n=56$, the Wasserstein distances among classical random bit-
strings and Zuchongzhi random bit-strings with cycle $m=12,14,16,18,20$.