#  Out-of-equilibrium Eigenstate Thermalization Hypothesis

Laura Foini<EMAIL_ADDRESS>IPhT, CNRS, CEA, Université Paris Saclay,
91191 Gif-sur-Yvette, France  Anatoly Dymarsky<EMAIL_ADDRESS>Department
of Physics and Astronomy, University of Kentucky, Lexington, KY 40503 USA
Silvia Pappalardi<EMAIL_ADDRESS>Institut für Theoretische
Physik, Universität zu Köln, Zülpicher Straße 77, 50937 Köln, Germany

(August 28, 2024)

###### Abstract

Understanding how out-of-equilibrium states thermalize under quantum unitary
dynamics is an important problem in many-body physics. In this work, we
propose a statistical ansatz for the matrix elements of non-equilibrium
initial states in the energy eigenbasis. The approach is inspired by the
Eigenstate Thermalisation Hypothesis (ETH) but the proposed ansatz exhibits
different scaling. Importantly, exponentially small cross-correlations between
the observable and the initial state matrix elements determine relaxation
dynamics toward equilibrium. We numerically verify scaling and cross-
correlation, point out the emergent universality of the high-frequency
behavior, and outline possible generalizations.

Introduction \- Over the past decades, the unitary evolution of nonequilibrium
states, including post-quench dynamics, has been a prominent subject in the
field of quantum dynamics. The mechanism for thermalization is now well
understood via the Eigenstate Thermalization Hypothesis (ETH) [1; 2; 3; 4].
The latter is a statistical ansatz for the matrix elements of physical
observables $\hat{A}$ is the energy eigenbasis
$\hat{H}|E_{i}\rangle=E_{i}|E_{i}\rangle$:

$A_{ij}=\mathcal{A}(E^{+})\delta_{ij}+e^{-S(E^{+})/2}f_{A}(E^{+},\omega_{ij})R_{ij}\
,$ (1)

with $E^{+}=(E_{i}+E_{j})/2$, $\omega_{ij}=E_{i}-E_{j}$ being the average
energy and frequency, $S(E)$ is thermodynamic entropy, and $R_{ij}$ is a
pseudorandom variable, such that $\overline{R_{ij}}=0$ and
$\overline{R_{ij}R_{ji}}=1$. Finally, $\mathcal{A}(E)$ and $f_{A}(E,\omega)$
are smooth functions of their arguments. This ansatz has proved to be
extremely successful in describing the equilibrium dynamics [5; 4] of physical
local Hamiltonians, as was shown by extensive numerical calculations [6; 7; 8;
9; 10; 11; 12; 13; 14]. Recently, the study of correlations between matrix
elements [15] has led to novel developments beyond the standard framework,
connecting ETH with Free Probability theory [15; 16; 17; 18], random matrix
universality [19; 20; 21; 22; 19; 23; 24; 25], conformal field theories [26]
and motivating the study of energy eigenvectors statistics [27; 28; 29; 30;
31; 32; 33; 34].

One of the central questions is how to extend the ETH framework to describe
non-equilibrium dynamics [2; 35; 36; 37; 38; 39; 40; 41]. In this work, we
propose a statistical ansatz for the matrix elements of the projector on the
initial out-of-equilibrium state $\Psi=\ket{\psi}\bra{\psi}$ written in the
eigenbasis of the Hamiltonian. Notably, the non-equilibrium dynamics are
encoded in the _correlations between the initial state and the observable’s
off-diagonal matrix elements_ , which we describe in our framework. After
introducing the ansatz and verifying its consistency, we discuss its
implications for the relaxation dynamics towards equilibrium and numerically
verify it in a non-integrable one-dimensional spin chain.

Set up \- The dynamics of a local observable can be written in the eigenbasis
of the Hamiltonian as

$\langle\psi|\hat{A}(t)|\psi\rangle=\sum_{ij}c_{i}c_{j}^{*}\,A_{ij}e^{i(E_{i}-E_{j})t}\
.$ (2)

with $c_{i}=\langle\psi|E_{i}\rangle$. The original ETH (1) is designed to
describe the stationary equilibrium point. In the absence of degeneracies, the
expectation value of $A$ eventually attains a stationary value

$\sum_{i}|c_{i}|^{2}\,A_{ii}=\langle\hat{A}\rangle_{\rm diag}\ ,$ (3)

which can be described by standard statistical mechanics. Namely, one
introduces the diagonal ensemble $\hat{\rho}_{\rm
diag}=\sum_{i}|c_{i}|^{2}|E_{i}\rangle\langle E_{i}|$ such that
$\langle\hat{A}\rangle_{\rm diag}=\text{Tr}\left(\hat{A}\,\hat{\rho}_{\rm
diag}\right)$ [3; 42]. In this work, we consider pure initial states with
extensive mean energy and sub-extensive energy fluctuations in the number of
degrees of freedom $N$

$\begin{array}[]{l}\displaystyle\langle\psi|\hat{H}|\psi\rangle=E_{0}\simeq
e_{0}N,\\\ \vspace{-0.2cm}\\\
\sqrt{\langle\psi|(\hat{H}-E_{0})^{2}|\psi\rangle}\simeq\delta_{e_{0}}N^{a},\qquad
a<1.\displaystyle\end{array}\ .$ (4)

For such initial states, the stationary value of $\langle\hat{A}(t)\rangle$ is
given by the microcanonical expectation, that, combined with ETH implies
thermalization, i.e. $\langle\hat{A}\rangle_{\rm
diag}\simeq\mathcal{A}(E_{0})$ [2; 3]. As a main example, we consider the case
$a=1/2$, satisfied if one performs a global quench, which also characterizes
equilibrium ensembles. Nonetheless, we will discuss the validity of our ansatz
also for other initial states (see the Discussions).

The fundamental object that we want to characterize is the projector on the
initial state written in the basis of the Hamiltonian

$\Psi_{ij}=c_{i}c_{j}^{\ast}=\langle
E_{j}|\psi\rangle\langle\psi|E_{i}\rangle\ .$ (5)

We will treat it as a pseudorandom object, analogously to $A_{ij}$ in ETH. A
crucial difference, in comparison with $\hat{A}$, is that this operator is of
rank one and that each off-diagonal matrix element is the product of two
pseudo-random numbers. This will radically change the scaling in the proposed
ansatz. Crucially, to capture the out-of-equilibrium dynamics, we will assume
that correlations exist between $\Psi$ and $\hat{A}$, when expressed in the
energy eigenbasis.

Ansatz \- We introduce an ansatz for the matrix $\Psi$,

$\Psi_{ij}\simeq\frac{e^{-\Phi(E_{i})}}{Z}\delta_{ij}+\frac{e^{-\frac{1}{2}(\Phi(E_{i})+\Phi(E_{j}))}}{Z}\,\tilde{R}_{ij}\
,$ (6)

where $\Phi(E)$ is a smooth function of energy and
$Z=\sum_{i}e^{-\Phi(E_{i})}$ the normalization, which defines the diagonal
ensemble:

$\overline{\Psi_{ii}}=|c_{i}|^{2}=\frac{e^{-\Phi(E_{i})}}{Z}\ .$ (7)

In Eq.(6), $\tilde{R}_{ij}$ are pseudorandom variables with zero average and
unit variance, i.e.

$\overline{\tilde{R}_{ij}}=0\
,\quad\overline{\tilde{R}_{ij}^{2}}=1,\qquad\text{for $i\neq j$}.$

Diagonal pseudo-random $\tilde{R}_{ii}$ also have zero average, but the
particular value of their variance may depend on the symmetry class of
$\hat{H}$ (e.g. GOE or GUE), as it is the case for $R_{ii}$ in standard ETH
[43; 19]. Given that $\Psi_{ii}$ is positive, $(1+R_{ii})\geq 0$. Since
$\Psi_{ij}$ is a product of two quasi-random numbers, this implies various
constraints on the joint properties of $\tilde{R}_{ij}$. Furthermore, these
variables are exponentially weakly correlated with $R_{ij}$ of the original
ETH ansatz (1),

$\overline{R_{ij}\tilde{R}_{ji}}=g_{A,\Psi}(e^{+}_{ij,},\omega_{ij})\,e^{-S(e^{+}_{ij})/2}\
,$ (8)

where $g_{A,\Psi}(e^{+}_{ij,},\omega_{ij})$ is an order one smooth function of
its variables, which describes the correlations crucial for non-equilibrium
dynamics. The existence of such correlations between the off-diagonal products
of $\Psi_{ij}A_{ji}$, can be shown by averaging over the (random) phases of
the eigenvectors of $\hat{H}$, see [44].

In particular, for the initial states conforming to Eq.(4), we will assume
that a large deviation scaling of the form $\Phi(E)=N\phi(e=E/N)$ applies,
such that the following is a _convex function_ 111This implies
$S^{\prime\prime}-\Phi^{\prime\prime}<0$, i.e.
$\phi^{\prime\prime}(e)>s^{\prime\prime}(e)$.,

$S(E)-\Phi(E)=N[s(e)-\phi(e)]\ .$ (9)

Summarising, in the _out-of-equilibrium ETH_ , observables and the initial
state look like pseudorandom matrices with smooth statistical properties
describing correlations or variance of the off-diagonal matrix elements,

$\displaystyle\overline{|A_{ij}|^{2}}$
$\displaystyle=e^{-Ns(e^{+})}\,|f_{A}(e^{+},\omega_{ij})|^{2},$ (10a)
$\displaystyle\overline{|\Psi_{ij}|^{2}}$
$\displaystyle=\frac{e^{-2N\phi(e^{+})}}{Z^{2}}\,e^{-\frac{\phi^{\prime\prime}(e^{+})}{4N}\omega_{ij}^{2}},$
(10b) $\displaystyle\overline{\Psi_{ij}A_{ij}}$
$\displaystyle=e^{-Ns(e^{+})}\,\frac{e^{-N\phi(e^{+})}}{Z}f_{A}(e^{+},\omega_{ij})\,g_{A,\Psi}(e^{+},\omega_{ij}),$
(10c)

where we made explicit the dependence on the system size $N$. Eq.(10b) applies
to the states conforming to (4) and it is obtained by expanding
$E_{i,j}=E^{+}\pm\omega_{ij}/2$ around $E^{+}$ using the assumption of large
deviation. A similar term,
$e^{-\frac{\phi^{\prime\prime}(e^{+})}{8N}\omega_{ij}^{2}}$, should appear
also in Eq.(10c), however, in the limit of large $N$ this can be neglected
because $f_{A}(e,\omega)$ is expected to decay at large frequencies.
Finally, we note that the product $\Psi_{ij}A_{ij}$ has large fluctuations
compared to the average (10c), see [44].

Consistency checks –  Let us first see how the ansatz (6) and in particular
the diagonal ensemble derived from that, satisfy the assumptions (4). The
large deviation scaling leads to an ensemble strongly peaked around the
characteristic (extensive) energy which maximizes (9) and with sub-extensive
fluctuations. In fact, in the large $N$ limit, the energy uncertainty reads:

$\displaystyle\Delta^{2}_{E_{0}}\equiv\langle\psi|(\hat{H}-E_{0})^{2}|\psi\rangle=\frac{1}{\Phi^{\prime\prime}(E_{0})-S^{\prime\prime}(E_{0})}\
.$ (11)

Owing to the extensivity in Eq.(9), this implies that $\Delta_{E}$ is sub-
extensive, in particular, for a post-quench state,

$\Delta_{E_{0}}=\delta_{e_{0}}\sqrt{N}\ ,$

where $\delta_{e_{0}}$ is an order-one constant, determined by the shape of
the large deviation 222Second derivative of an extensive quantity with respect
to the energy is
$\Phi^{\prime\prime}(E_{0})=\frac{\partial^{2}\Phi(E)}{\partial
E^{2}}\Big{|}_{E=E_{0}}=\frac{1}{N}\frac{\partial^{2}\phi(e)}{\partial
e^{2}}\Big{|}_{e=e_{0}}=\frac{1}{N}\phi^{\prime\prime}(e_{0}).$ Therefore
$\delta_{e_{0}}=1/\sqrt{\phi^{\prime\prime}(e_{0})-s^{\prime\prime}(e_{0})}$.
We shall now proceed to discuss a set of consistency checks to validate our
proposed approach.

_Normalization –_ The state normalisation $\text{Tr}\,\Psi=1$ is ensured by
the definition of $Z$. In the large $N$ limit, the “partition function” $Z$ in
Eq.(6) reads:

$\displaystyle Z$
$\displaystyle={\sqrt{2\pi}}{\Delta_{E_{0}}}e^{S(E_{0})-\Phi(E_{0})}\ ,$ (12)

with $\Delta_{E_{0}}$ given by Eq.(11). However one can show a stronger
property, namely that $\text{Tr}\,\Psi^{2}=1$. See Eq.(16) below at time zero.

_Projector –_ We now discuss an even tighter constraint: the projector
identity $\Psi^{2}=\Psi$ at the level of individual matrix elements. For our
ansatz, this turns out to be true in a statistical way in the thermodynamic
limit. In particular, thinking of the matrix elements $\Psi_{ij}$ as products
of two random variables we assume the following properties:

$\begin{array}[]{ll}\displaystyle\tilde{R}_{ik}\tilde{R}_{kj}=\tilde{R}_{ij}(1+\tilde{R}_{kk}),&\text{for
$i\neq j\neq k$},\\\ \vspace{-0.2cm}\\\
\displaystyle\tilde{R}_{ij}\tilde{R}_{ji}=(1+\tilde{R}_{ii}+\tilde{R}_{jj}+\tilde{R}_{ii}\tilde{R}_{jj}),&\text{for
$i\neq j$}.\end{array}$ (13)

At the leading order in $N$, this implies, see [44],

$\displaystyle[\Psi^{2}]_{ij}$
$\displaystyle\simeq\frac{e^{-\frac{1}{2}(\Phi(E_{i})+\Phi(E_{j}))}}{Z}\tilde{R}_{ij}\simeq[\Psi]_{ij}\quad
i\neq j,$ (14a) $\displaystyle[\Psi^{2}]_{ii}$
$\displaystyle\simeq\frac{e^{-\Phi(E_{i})}}{Z}\left(1+\tilde{R}_{ii}\right)\simeq[\Psi]_{ii}.$
(14b)

Therefore the ansatz preserves its structure upon multiplication.

Implications for the dynamics –  We now discuss the main motivation behind our
ansatz, designed to describe equilibration dynamics of physical observables.

_Fidelity decay –_ First of all, we show that our ansatz is consistent with
the expected behavior of the fidelity decay (survival probability) [47; 48;
49; 50], defined as

$\displaystyle|\langle\psi|\psi(t)\rangle|^{2}$
$\displaystyle=\left|\sum_{i}|c_{i}|^{2}e^{-iE_{i}t}\right|^{2}\ .$ (15)

By substituting sums with integrals, neglecting spectral correlations, and
using the out-of-equilibrium ETH ansatz in Eq.(10b), at the leading order in
$N$, one finds

$\displaystyle\begin{split}|\langle\psi|\psi(t)\rangle|^{2}\simeq\frac{1}{2\Delta_{E_{0}}\sqrt{\pi}}\int
d\omega e^{-\frac{1}{2}\frac{1}{2\Delta^{2}_{E_{0}}}\omega^{2}}e^{i\omega
t}=e^{-t^{2}\,\Delta^{2}_{E_{0}}}\ ,\end{split}$ (16)

where we have used the definition of the energy variance in Eq.(11), i.e.
$\Delta_{E_{0}}=1/\sqrt{\Phi^{\prime\prime}(E_{0})-S^{\prime\prime}(E_{0})}=\delta_{e_{0}}\sqrt{N}$.
Thus the large deviation ansatz in Eq.(6) is consistent with the Gaussian
decay, controlled by the energy variance of the initial state, in agreement
with the literature on global quenches, see e.g. Ref. [48]. Dynamical
behavior, different from Eq.(16), e.g. including an exponential decay, is
known to arise from the initial states which are different from Eq.(4) [47;
51; 52; 48; 49; 53], as discussed below.

_Relaxation dynamics –_ The primary purpose of this work is the study
relaxation dynamics in Eq.(2), namely:

$\delta
A_{\Psi}(t)=\langle\psi|\hat{A}(t)|\psi\rangle-\langle\hat{A}\rangle_{\rm
diag}=\sum_{i\neq j}c_{i}c_{j}^{*}\,A_{ij}e^{i(E_{i}-E_{j})t}\ .$ (17)

Plugging Eq.(10c) into Eq.(17), the standard ETH manipulations lead to

$\displaystyle\delta A_{\Psi}(t)$ $\displaystyle\simeq\int
d\omega\,f_{A}(e_{0},\omega)g_{A,\Psi}(e_{0},\omega)e^{\frac{\omega^{2}}{4N\delta^{2}_{e_{0}}}}e^{-i\omega
t}\ .$ (18)

Hence, the correlation between the initial state and the operator encodes the
Fourier transform of the relaxation:

$\tilde{\delta A_{\Psi}}(\omega)=f_{A}(e_{0},\omega)g_{A,\Psi}(e_{0},\omega)\
,$ (19)

where we have neglected the Gaussian frequency term for $N\gg 1$. Thus, the
out-of-equilibrium behavior is encoded in this function and will depend, in
general, on the correlations between the state and the observable.

The relaxation dynamics (17) share some properties with the (two-time)
dynamical correlations at thermal equilibrium for the same observable. One has
[5; 4]:

$C(t)=\frac{1}{2}\langle\\{A(t),A(0)\\}\rangle_{c}=\int\,{\rm
d}\omega\,e^{i\omega
t}\cosh\left(\frac{\beta\omega}{2}\right)f_{A}^{2}(e_{\beta},\omega)\ ,$ (20)

where $\langle\cdot\rangle=\mathrm{Tr}(e^{-\beta
H}\cdot)/\mathrm{Tr}(e^{-\beta H})$ and $e_{\beta}=\langle H\rangle/N$.
Therefore the ETH function $f_{A}(e,\omega)$ enters both Eqs. (18) and (20)
and its properties in the $\omega\to 0$ limit control the long-time behavior.
This fact is usually invoked in the literature, see e.g. [2; 5], and our
ansatz in Eq. (18) makes it explicit.
Similarly to $f_{A}^{2}(\omega)$, which has to decay exponentially at large
$\omega$ in $D\geq 2$ (superexponentially in 1D), high-frequency tail of
$g_{\Psi,A}(\omega)$ has to be exponentially suppressed for states $\Psi$
associated with local perturbations, see [44].

Let us now comment on some differences between Eq. (18) and (20). The
integrand in (20) is positive-definite. As a result, $C(t)$ necessarily decays
at early times. On the contrary, the integrand in Eq. (18) is not sign-
definite, hence $\delta A_{\Psi}(t)$ can both increase or decrease throughout
relaxation dynamics.

Numerical results - We test the predictions above in the case of the one-
dimensional Ising model with a tilted field

$H=\sum_{i=1}^{L}w\sigma_{i}^{x}+\sum_{i=1}^{L}h\sigma_{i}^{z}+\sum_{i=1}^{L-1}J\sigma_{i}^{z}\sigma_{i+1}^{z}$
(21)

with $w=\sqrt{5}/2$, $h=(\sqrt{5}+5)/8$ and $J=1$ and consider different local
single or two sites observables,

$\hat{A}=\sigma^{x}_{1},\,\hat{A}=\sigma^{z}_{1}\,\quad\text{or}\quad\hat{A}=\sigma^{z}_{1}\sigma^{z}_{2}\
.$ (22)

We consider simple out-of-equilibrium initial states, fully polarized states

$|\psi\rangle=\ket{\downarrow_{\alpha}\downarrow_{\alpha}\dots\downarrow_{\alpha}}$
(23)

in the $\alpha=z$ or $\alpha=y$ directions. We impose periodic boundary
conditions on the Hamiltonian in Eq.(21) and restrict the analysis to
translationally-invariant sector $k=0$ with positive parity reflection
symmetry. As a technical tool, we use the smoothed average of our energy-
resolved data as

$\overline{[f(x)]}_{\tau}=\frac{\sum_{n}f(x_{n})\delta_{\tau}(x-x_{n})}{\sum_{n}\delta_{\tau}(x-x_{n})}\
,$ (24)

where, $\delta_{\tau}(x)$ is a smoothed delta functions such that
${\lim_{\tau\to\infty}\delta_{\tau}(x)=\delta(x)}$. In the simulations, we
chose a Gaussian smoothing
$\delta_{\tau}(x)=e^{-\frac{\tau^{2}}{2}x^{2}}/\sqrt{2\pi/\tau^{2}}$.

Figure 1: Out-of-Equilibrium ETH of the fully polarized initial state
${|\psi\rangle=|\downarrow_{z}\dots\downarrow_{z}\rangle}$. (a) The rescaled
diagonal ensemble as a function of the energy density for different system
sizes $L=12,14,16,18$. In the inset, the initial energy
$E_{0}=\langle\psi|\hat{H}|\psi\rangle$ and the energy fluctuations $(\Delta
E_{0})^{2}=\langle\psi|\hat{H}^{2}|\psi\rangle-E_{0}^{2}$ are plotted as a
function of the inverse of the system size $1/L$. In both panels, the blue
dots correspond to individual overlaps for $L=16$. The smoothing parameter is
$\tau=4$. In the inset, the numerical fidelity decay up to $L=22$ is compared
with the prediction in Eq.(16) (dashed) without fitting parameter.

First, we establish that the initial states (23) are consistent with the
ansatz in Eq.(6). In Fig. 1a we plot the diagonal ensemble for different
length sizes $L=12,14,16,18$, showing that it obeys the large deviation
prediction $\overline{\Psi_{ii}}=\frac{e^{-L\phi(E_{i}/L)}}{Z}$. This is
confirmed by the inset, where we plot the scaling of the initial energy
$E_{0}=e_{0}L$ and variance
$\Delta_{E_{0}}/E_{0}=\delta_{e_{0}}/e_{0}\sqrt{L}$, c.f. Eq.(11). From a fit
of the data, we extract the dimensionless values $e_{0}=0.10$,
$\delta_{e_{0}}=1.12$. In panel (b), we study the fluctuations of the out-of-
equilibrium ETH functions in the frequency domain [cf. Eq.(10b)]. To address
the dependence on $L$, we re-scale $\overline{|\Psi_{ij}|^{2}}$ by the
diagonal matrix elements at energy $e_{0}$. With this choice, the ansatz in
Eq.(10b) is equivalent to

$\frac{\overline{\Psi_{ij}\Psi_{ji}}}{\overline{{|\Psi_{e^{+}e^{+}}|}}^{2}}\simeq
e^{-\frac{\phi^{\prime\prime}(e^{+})}{4L}\omega_{ij}^{2}},\quad
E^{+}=Ne^{+}={E_{i}+E_{j}\over 2}\ .$ (25)

In Fig.1b, we fix the energy density to be $e_{0}$ by rescricting the energy
indices $i,j$ of $\overline{\Psi_{ij}\Psi_{ji}}$ to
$|(E_{i}+E_{j})/2-E_{0}|\leq\sqrt{L}\delta_{e_{0}}$. The figure shows the
smoothed average (25) as a function of the energy difference
$\omega_{ij}=E_{i}-E_{j}$ rescaled by $\sqrt{L}$, for different system sizes.
For $L=16$ we also show individual values without smoothing (blue dots). The
plot confirms that this initial state has fluctuations that decay as a
Gaussian with a variance $1/\sqrt{L}$, consistent with Eq.(25). In the inset,
we also confirm the Gaussian decay of the fidelity upon increasing system size
[cf. Eq.(16)]: we plot $e^{-\delta_{e_{0}}^{2}t^{2}}$ with
$\delta_{e_{0}}=1.12$ without fitting parameter.

Figure 2: ETH correlations between the initial state and the observable
$A=\hat{\sigma}^{x}_{1}$. (a) Scaling with the system size of the ETH
predictions in Eqs.(10). The red and black dashed lines indicate
$(\dim\mathcal{H})^{-1}$ and $(\dim\mathcal{H})^{-2}$ respectively (b)
Smoothed averages $\overline{A_{ij}\Psi_{ij}}$ describing the Fourier
transform of the relaxation dynamics increasing system size. Smoothed
quantities at energy density $e_{0}=0.10$ with $\tau=4$.

We then proceed to establish the validity of the ansatz for the correlations
between the initial state and observable $\hat{A}$ in the energy eigenbasis.
In Fig.2, we focus on
${|\psi\rangle=|\downarrow_{z}\dots\downarrow_{z}\rangle}$ and
$\hat{A}=\hat{\sigma}_{1}^{x}$. In panel (a), we test the system size
dependency of Eqs.(10) at energy density $e_{0}$ for finite frequency
$\omega_{ij}=5.1$. As predicted by out-of-equilibrium ETH Eq.(10a), the
observable off-diagonal matrix elements $\overline{A_{ij}A_{ji}}$ decay as
$O(e^{-Ls(e_{0})})$, while both $\overline{|\Psi_{ij}|^{2}}$ and
$\overline{A_{ij}\Psi_{ji}}$ decay as $O(e^{-2L})$, cf. Eqs(10b)-(10c). The
red and black dashed lines indicate $(\dim\mathcal{H})^{-1}$ and
$(\dim\mathcal{H})^{-2}$ respectively. We checked that the same results hold
at zero or for other finite $\omega_{ij}$.
In Fig.2b, we consider

$\frac{\overline{A_{ij}\Psi_{ij}}}{\overline{\Psi_{e_{0}e_{0}}}^{2}}\simeq
f_{A}(e_{0},\omega)g_{A,\Psi}(e_{0},\omega_{ij})\ ,$ (26)

where $\Psi_{e_{0}e_{0}}$ is the same as in Eq.(25) and the right-end side
follows from Eq.(10c) and (7) 333Using that the diagonal ensemble (7) for
$e^{+}=e_{0}$ reads ${|c_{e_{0}}|^{2}=\frac{e^{-\Phi(E_{0})}}{Z}\simeq
e^{-S(E_{0})}}$ where one uses the normalization in Eq.(12). This quantity is
of order one, i.e. it remains finite in the thermodynamic limit. Its Fourier
transform yields equilibration dynamics, see Eq.(18). Note that we have
plotted the absolute value of Eq.(26), since the sign of $g(e,\omega)$
oscillates. This sign change is a characteristic feature of the out-of-
equilibrium dynamics, as was emphasized above.

Figure 3: Absolute value of the correlations between initial state and
observable $|g_{A,\Psi}(e_{0},\omega)|$ extracted using Eq.(27) for different
observables $A=\sigma_{1}^{x}$, $\sigma_{1}^{z}$ and
$\sigma_{1}^{z}\sigma_{2}^{z}$ as a function of frequency. (a) Results from
the initial state $\ket{\downarrow_{z}\dots\downarrow_{z}}$. (b) Results from
the initial state $\ket{\downarrow_{y}\dots\downarrow_{y}}$. Here $L=18$ and
$\tau=4$.

To better understand the behaviour of the function $g(e_{0},\omega)$, in Fig.3
we consider:

$\frac{\overline{A_{ij}\Psi_{ij}}}{\sqrt{\overline{A_{ij}^{2}}}\,\sqrt{\overline{|\Psi_{ij}|^{2}}}\sqrt{\overline{\Psi_{e_{0}e_{0}}}}}\simeq
g_{A,\Psi}(e_{0},\omega_{ij})\ ,$ (27)

to obtain an order one quantity, which encodes the correlations in Eq.(8). The
results are shown in Fig.3 for the fully polarized states in Eq.(23) along the
directions $\alpha=z,y$ in panels (a) and (b) respectively, or the three
different operators in Eq.(22). The plot shows that the
$g_{A,\Psi}(e_{0},\omega)$ may still decay, albeit slowly, as a function of
frequency. The most notable fact is that, for different observables, the
smooth functions $g$ have approximately the same behavior at large
frequencies, which does not depend on the observable. This seems to indicate
that the large frequency behavior and the oscillations in the
$g_{\psi}(e_{0},\omega)$ reflect physics _of the initial state_.

Discussion and conclusions \- In this paper, we have introduced a new ansatz
for out-of-equilibrium dynamics, which predicts correlations between the
initial state and observables when written in the energy eigenbasis.

Let us remark that our results describe _a wide class of initial states_ ,
including for example products $|a\rangle_{A}|b\rangle_{B}$ of energy
eigenstates $|a\rangle_{A}$ and $|b\rangle_{B}$ of subsystems $A$ and $B$,
that have recently motivated studies of the eigenstate correlations [29; 30;
31; 32; 33; 34]. For local Hamiltonians, these states have only intensive
energy fluctuations $\Delta_{E}^{2}=O(1)$. While obeying Eq.(4), they do not
have the form of a large deviation (9) and their survival probability (15) is
known to decay exponentially in time [48]. Nevertheless, this class of states
is naturally included in our ansatz on the state-observable correlations in
Eqs.(10), with a difference in Eq.(10b) which generically reads

$\overline{|\Psi_{ij}|^{2}}=e^{-\Phi(E^{+}+\omega/2)+\Phi(E^{+}-\omega/2)}/Z^{2}.$
(28)

These state are discussed in [44], where we verify numerically the general
scaling with the system size of Eqs.(10), and comment on the relation with the
literature.

Our work opens a series of perspectives. At long times, hydrodynamic modes are
expected to play a dominant role in equilibration dynamics [55; 56; 57], and
it would be valuable to investigate how hydrodynamic description can be
incorporated into the non-equilibrium ETH. Additionally, one could explore how
the current ETH framework applies to integrable systems that equilibrate to a
generalized Gibbs ensemble [11; 42; 58; 59] or in the presence of many-body
quantum scars [60; 61; 62; 63]. Further questions concern the natural
generalization of our ansatz to cross-correlations between different states
$c_{i}^{a}=\langle E_{i}|\psi_{a}\rangle$ for multiple choices of $a$ and $i$
and in the higher-order correlations between states and observables [16; 33;
34].

###### Acknowledgements.

We thank J. Kurchan for inspiring discussions and collaboration on related
topics. We thank A. Rosch and D. Abanin for useful discussions. AD is
supported by the NSF grant PHY 2310426. S.P. acknowledges support by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy - Cluster of Excellence Matter and Light for
Quantum Computing (ML4Q) EXC 2004/1 -390534769.

## References

  * Deutsch [1991] J. M. Deutsch, Physical Review A 43, 2046 (1991).
  * Srednicki [1999] M. Srednicki, Journal of Physics A: Mathematical and General 32, 1163 (1999).
  * Rigol  _et al._ [2008a] M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854–858 (2008a).
  * D’Alessio  _et al._ [2016] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, Advances in Physics 65, 239 (2016).
  * Khatami  _et al._ [2013] E. Khatami, G. Pupillo, M. Srednicki, and M. Rigol, arXiv preprint arXiv:1304.7279 (2013).
  * Prosen [1999] T. Prosen, Phys. Rev. E 60, 3949 (1999).
  * Rigol  _et al._ [2008b] M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854 (2008b).
  * Biroli  _et al._ [2010] G. Biroli, C. Kollath, and A. M. Läuchli, Physical review letters 105, 250401 (2010).
  * Ikeda  _et al._ [2013] T. N. Ikeda, Y. Watanabe, and M. Ueda, Physical Review E 87, 012125 (2013).
  * Steinigeweg  _et al._ [2013] R. Steinigeweg, J. Herbrych, and P. Prelovšek, Physical Review E 87, 012118 (2013).
  * Alba [2015] V. Alba, Physical Review B 91, 155123 (2015).
  * Beugeling  _et al._ [2015] W. Beugeling, R. Moessner, and M. Haque, Physical Review E 91, 012144 (2015).
  * Luitz [2016] D. J. Luitz, Physical Review B 93, 134201 (2016).
  * Fritzsch and Prosen [2021] F. Fritzsch and T. Prosen, Phys. Rev. E 103, 062133 (2021).
  * Foini and Kurchan [2019a] L. Foini and J. Kurchan, Physical Review E 99, 042139 (2019a).
  * Pappalardi  _et al._ [2022] S. Pappalardi, L. Foini, and J. Kurchan, Phys. Rev. Lett. 129, 170603 (2022), arXiv:2204.11679 [cond-mat.stat-mech] .
  * Pappalardi  _et al._ [2023] S. Pappalardi, F. Fritzsch, and T. Prosen, arXiv preprint arXiv:2303.00713 (2023).
  * Fava  _et al._ [2023] M. Fava, J. Kurchan, and S. Pappalardi, arXiv preprint arXiv:2308.06200 (2023).
  * Foini and Kurchan [2019b] L. Foini and J. Kurchan, Phys. Rev. Lett. 123, 260601 (2019b).
  * Richter  _et al._ [2020] J. Richter, A. Dymarsky, R. Steinigeweg, and J. Gemmer, Physical Review E 102 (2020).
  * Wang  _et al._ [2021] J. Wang, M. H. Lamann, J. Richter, R. Steinigeweg, A. Dymarsky, and J. Gemmer, arXiv preprint arXiv:2110.04085 (2021).
  * Brenes  _et al._ [2021] M. Brenes, S. Pappalardi, M. T. Mitchison, J. Goold, and A. Silva, Physical Review E 104 (2021).
  * Jafferis  _et al._ [2023] D. L. Jafferis, D. K. Kolchmeyer, B. Mukhametzhanov, and J. Sonner, Physical Review X 13, 031033 (2023).
  * Jafferis  _et al._ [2022] D. L. Jafferis, D. K. Kolchmeyer, B. Mukhametzhanov, and J. Sonner, arXiv preprint arXiv:2209.02131 (2022).
  * Wang  _et al._ [2023] J. Wang, J. Richter, M. H. Lamann, R. Steinigeweg, J. Gemmer, and A. Dymarsky, arXiv preprint arXiv:2310.20264 (2023).
  * de Boer  _et al._ [2024a] J. de Boer, D. Liska, and B. Post, arXiv preprint arXiv:2405.13111 (2024a).
  * Dymarsky  _et al._ [2018] A. Dymarsky, N. Lashkari, and H. Liu, Phys. Rev. E 97, 012140 (2018), arXiv:1611.08764 [cond-mat.stat-mech] .
  * Chan  _et al._ [2019] A. Chan, A. De Luca, and J. T. Chalker, Phys. Rev. Lett. 122, 220601 (2019).
  * Huang [2019] Y. Huang, Nuclear Physics B 938, 594 (2019).
  * Lu and Grover [2019] T.-C. Lu and T. Grover, Phys. Rev. E 99, 032111 (2019).
  * Murthy and Srednicki [2019] C. Murthy and M. Srednicki, Physical Review E 100, 022131 (2019).
  * Shi  _et al._ [2023] Z. D. Shi, S. Vardhan, and H. Liu, Phys. Rev. B 108, 224305 (2023).
  * Hahn  _et al._ [2023] D. Hahn, D. J. Luitz, and J. Chalker, arXiv preprint arXiv:2309.12982 (2023).
  * Jindal and Hosur [2024] S. Jindal and P. Hosur, arXiv preprint arXiv:2401.13829 (2024).
  * Richter  _et al._ [2019] J. Richter, J. Gemmer, and R. Steinigeweg, Physical Review E 99, 050104 (2019).
  * Dymarsky [2019] A. Dymarsky, Physical Review B 99, 224302 (2019).
  * Knipschild and Gemmer [2020] L. Knipschild and J. Gemmer, Physical Review E 101, 062205 (2020).
  * Reimann and Gemmer [2020] P. Reimann and J. Gemmer, Physica A: Statistical Mechanics and its Applications 552, 121840 (2020), tributes of Non-equilibrium Statistical Physics.
  * Dymarsky [2022] A. Dymarsky, Physical Review Letters 128, 190601 (2022).
  * Capizzi  _et al._ [2024] L. Capizzi, J. Wang, X. Xu, L. Mazza, and D. Poletti, (2024), arXiv:2405.16975 [quant-ph] .
  * de Boer  _et al._ [2024b] J. de Boer, D. Liska, B. Post, and M. Sasieta, Journal of High Energy Physics 2024, 1 (2024b).
  * Essler and Fagotti [2016] F. H. Essler and M. Fagotti, Journal of Statistical Mechanics: Theory and Experiment 2016, 064002 (2016).
  * Mondaini and Rigol [2017] R. Mondaini and M. Rigol, Physical Review E 96, 012157 (2017).
  * [44] See the Supplemental Material for additional details on the derivations.
  * Note [1] This implies $S^{\prime\prime}-\Phi^{\prime\prime}<0$, i.e. $\phi^{\prime\prime}(e)>s^{\prime\prime}(e)$.
  * Note [2] Second derivative of an extensive quantity with respect to the energy is  $\Phi^{\prime\prime}(E_{0})=\frac{\partial^{2}\Phi(E)}{\partial E^{2}}\leavevmode@ifvmode{\hbox{$\left|\vbox to1.5pt{}\right.$}}_{E=E_{0}}=\frac{1}{N}\frac{\partial^{2}\phi(e)}{\partial e^{2}}\leavevmode@ifvmode{\hbox{$\left|\vbox to1.5pt{}\right.$}}_{e=e_{0}}=\frac{1}{N}\phi^{\prime\prime}(e_{0}).$ Therefore $\delta_{e_{0}}=1/\sqrt{\phi^{\prime\prime}(e_{0})-s^{\prime\prime}(e_{0})}$.
  * Gorin  _et al._ [2006] T. Gorin, T. Prosen, T. H. Seligman, and M. Žnidarič, Physics Reports 435, 33 (2006).
  * Torres-Herrera and Santos [2014] E. Torres-Herrera and L. F. Santos, Physical Review E 89, 062110 (2014).
  * Schiulaz  _et al._ [2019] M. Schiulaz, E. J. Torres-Herrera, and L. F. Santos, Phys. Rev. B 99, 174313 (2019).
  * Rai  _et al._ [2023] K. S. Rai, J. I. Cirac, and Á. M. Alhambra, arXiv preprint arXiv:2307.05200 (2023).
  * Santos  _et al._ [2012a] L. F. Santos, F. Borgonovi, and F. Izrailev, Physical review letters 108, 094102 (2012a).
  * Santos  _et al._ [2012b] L. F. Santos, F. Borgonovi, and F. M. Izrailev, Phys. Rev. E 85, 036209 (2012b).
  * Micklitz  _et al._ [2022] T. Micklitz, A. Morningstar, A. Altland, and D. A. Huse, Phys. Rev. Lett. 129, 140402 (2022).
  * Note [3] Using that the diagonal ensemble (7) for $e^{+}=e_{0}$ reads ${|c_{e_{0}}|^{2}=\frac{e^{-\Phi(E_{0})}}{Z}\simeq e^{-S(E_{0})}}$ where one uses the normalization in Eq.(12).
  * Mukerjee  _et al._ [2006] S. Mukerjee, V. Oganesyan, and D. Huse, Physical Review B 73, 035113 (2006).
  * Lux  _et al._ [2014] J. Lux, J. Müller, A. Mitra, and A. Rosch, Physical Review A 89, 053608 (2014).
  * Michailidis  _et al._ [2023] A. A. Michailidis, D. A. Abanin, and L. V. Delacrétaz, arXiv preprint arXiv:2310.10564 (2023), https://doi.org/10.48550/arXiv.2310.10564.
  * Vidmar and Rigol [2016] L. Vidmar and M. Rigol, Journal of Statistical Mechanics: Theory and Experiment 2016, 064007 (2016).
  * Essler and de Klerk [2023] F. Essler and A. de Klerk, arXiv preprint arXiv:2307.12410 (2023).
  * Serbyn  _et al._ [2021] M. Serbyn, D. A. Abanin, and Z. Papić, Nature Physics 17, 675 (2021).
  * Moudgalya  _et al._ [2022] S. Moudgalya, B. A. Bernevig, and N. Regnault, Reports on Progress in Physics 85, 086501 (2022).
  * Chandran  _et al._ [2023] A. Chandran, T. Iadecola, V. Khemani, and R. Moessner, Annual Review of Condensed Matter Physics 14, 443 (2023).
  * Gotta  _et al._ [2023] L. Gotta, S. Moudgalya, and L. Mazza, Phys. Rev. Lett. 131, 190401 (2023).
  * Maillard  _et al._ [2019] A. Maillard, L. Foini, A. L. Castellanos, F. Krzakala, M. Mézard, and L. Zdeborová, Journal of Statistical Mechanics: Theory and Experiment 2019, 113301 (2019).
  * Avdoshkin and Dymarsky [2020] A. Avdoshkin and A. Dymarsky, Physical Review Research 2, 043234 (2020).
  * Deutsch [2010] J. M. Deutsch, New Journal of Physics 12, 075021 (2010).

——————————

Supplemental Material:

In the Supplementary Material, we provide additional analysis and background
calculations to support the results in the main text: RMT-based argument for
scaling, consistency checks and a bound on high frequency behavior of the out-
of-equilibrium ETH ansatz; discussion of bipartite states, etc.

##  1 A simple scaling drawn from RMT

Let us discuss a simple example that exhibits the scaling in Eq.(10). Consider
an observable
$A=\sum_{\alpha}\lambda_{\alpha}|\lambda_{\alpha}\rangle\langle\lambda_{\alpha}|$
and as an initial state we will take an eigenvector of such observable
$|\psi\rangle=|\lambda_{\psi}\rangle$. Let us suppose that afterward, the
state evolves under a ${\mathcal{D}}\times{\mathcal{D}}$ Hamiltonian that is
drawn from a rotationally invariant ensemble, i.e. $P(H)=P(U^{-1}HU)$ where
$U$ is arbitrary orthogonal (or unitary) matrix, for instance a GOE or GUE
ensemble. With this choice, the Hamiltonian eigenvectors $|E_{j}\rangle$, in
the basis of the observable, i.e. $\langle E_{j}|\lambda_{\alpha}\rangle$, are
represented by random orthogonal or unitary matrices. The properties of the
matrix elements of a given observable $A$ in such random basis have been
discussed in [64]. In the large ${\mathcal{D}}$ limit assuming to initialise
the dynamics in $|\psi\rangle=|\lambda_{\psi}\rangle$, some eigenvector of the
observable $A$, the properties of this toy “out-of-equilibrium” ETH can be
easily derived

$\displaystyle\overline{|A_{ij}|^{2}}$
$\displaystyle=\frac{1}{{\mathcal{D}}}\,(\langle A^{2}\rangle-\langle
A\rangle^{2})$ (S1a) $\displaystyle\overline{|\Psi_{ij}|^{2}}$
$\displaystyle=\frac{1}{{\mathcal{D}}^{2}}\qquad\text{for $i\neq j$}$ (S1b)
$\displaystyle\overline{\Psi_{ij}A_{ij}}$
$\displaystyle=\frac{1}{{\mathcal{D}}^{2}}\left[\lambda_{\psi}-\langle
A\rangle\right]$ (S1c)

where $\langle\bullet\rangle=\frac{1}{{\mathcal{D}}}\mathrm{Tr}(\bullet)$.
This is a particularly simple example of the ansatz (10) discussed in the main
text. In term of normalized fluctuations this means
$\overline{R_{ij}\tilde{R}_{ij}}\simeq{\mathcal{D}}^{-{1/2}}$, as in (8). This
examples illustrates the difference in scaling between
$\overline{|A_{ij}|^{2}}$ and $\overline{|\Psi_{ij}|^{2}}$. The first quantity
has rank ${\mathcal{D}}$, leading to
$\sum_{ij}\overline{|A_{ij}|^{2}}=O({\mathcal{D}})$, while the second has rank
one, $\sum_{ij}\overline{|\Psi_{ij}|^{2}}=O(1)$. Similarly
$\sum_{ij}\overline{\Psi_{ij}A_{ij}}=O(1)$, which is consistent with the
scaling above.

##  2 A constraint on the ansatz

Let us justify Eqs. (13). As we stressed several times, contrary to the
standard ETH ansatz for observables, the matrix that we have are chracterising
has rank 1. In particular, calling
$z_{i}=\frac{1}{\sqrt{Z}}e^{-\frac{1}{2}\Phi(E_{i})}$ and following the
notation of the main text we have:

$\begin{array}[]{c}\displaystyle
1+\tilde{R}_{ii}=\frac{|c_{i}|^{2}}{z_{i}^{2}}\\\ \vspace{-0.2cm}\hfil\\\
\displaystyle\tilde{R}_{ij}=\frac{c_{i}c_{j}^{\ast}}{z_{i}z_{j}}\qquad\text{for
$i\neq j$}\end{array}$ (S2)

Taking products:

$\tilde{R}_{ik}\tilde{R}_{kj}=\frac{c_{i}|c_{k}|^{2}c_{j}^{\ast}}{z_{i}|z_{k}|^{2}z_{j}}=\tilde{R}_{ij}(1+\tilde{R}_{kk})\qquad\text{for
$i\neq j\neq k$}$ (S3)

and similarly

$\tilde{R}_{ij}\tilde{R}_{ji}=\frac{|c_{i}|^{2}|c_{j}|^{2}}{|z_{i}|^{2}|z_{j}|^{2}}=(1+\tilde{R}_{ii})(1+\tilde{R}_{jj})\qquad\text{for
$i\neq j$}$ (S4)

Let us now see how these constraints imply that $\Psi$ is a projector by
proving Eqs. (14).

We start by evaluating the off-diagonal with $i\neq j$:

$\displaystyle\begin{split}[\Psi^{2}]_{ij}&=\sum_{k}\Psi_{ik}\Psi_{kj}=\Psi_{ii}\Psi_{ij}+\Psi_{ij}\Psi_{jj}+\sum_{k:k\neq
i\neq j}\Psi_{ik}\Psi_{kj}\\\
&=\frac{1}{Z^{2}}e^{-\frac{1}{2}(\Phi(E_{i})+\Phi(E_{j}))}\left[e^{-\Phi(E_{i})}\tilde{R}_{ii}\tilde{R}_{ij}+e^{-\Phi(E_{j})}\tilde{R}_{jj}\tilde{R}_{ij}\right]+\frac{1}{Z}e^{-\frac{1}{2}(\Phi(E_{i})+\Phi(E_{j}))}\sum_{k:k\neq
i\neq j}\frac{e^{-\Phi(E_{k})}}{Z}\tilde{R}_{ik}\tilde{R}_{kj}\end{split}$
(S5)

where from the first to the second line we have inserted the ansatz (10) in
the individual matrix elements. The first term is subleading $O(e^{-2N})$,
while in the second term, we can substitute the ansatz of Eq.(S3) and obtain

$\displaystyle\begin{split}[\Psi^{2}]_{ij}&\simeq\frac{1}{Z}e^{-\frac{1}{2}(\Phi(E_{i})+\Phi(E_{j}))}\sum_{k:k\neq
i\neq
j}\frac{e^{-\Phi(E_{k})}}{Z}\tilde{R}_{ij}(1+\tilde{R}_{kk})=\frac{1}{Z}e^{-\frac{1}{2}(\Phi(E_{i})+\Phi(E_{j}))}\tilde{R}_{ij}\sum_{k:k\neq
i\neq j}\Psi_{kk}\\\
&=\frac{1}{Z}e^{-\frac{1}{2}(\Phi(E_{i})+\Phi(E_{j}))}\tilde{R}_{ij}=\Psi_{ij}\
,\end{split}$ (S6)

where we used $\sum_{k:k\neq i\neq
j}\Psi_{kk}\simeq\sum_{k}\Psi_{kk}=\langle\psi|\psi\rangle=1$ which shows the
first equation in (14). Similar manipulations can be done on the diagonal
elements, leading to

$\displaystyle\begin{split}[\Psi^{2}]_{ii}&=\sum_{k}\Psi_{ik}\Psi_{ki}=\Psi_{ii}\Psi_{ii}+\sum_{k:k\neq
i}\Psi_{ik}\Psi_{ki}\\\
&=\frac{e^{-2\Phi(E_{i})}}{Z^{2}}(1+\tilde{R}_{ii})^{2}+\frac{1}{Z}e^{-\Phi(E_{i})}\sum_{k:k\neq
i}\frac{e^{-\Phi(E_{k})}}{Z}\tilde{R}_{ik}\tilde{R}_{ki}\\\
&\simeq\frac{1}{Z}e^{-\Phi(E_{i})}(1+\tilde{R}_{ii})\sum_{k:k\neq
i}\frac{e^{-\Phi(E_{k})}}{Z}(1+\tilde{R}_{kk})=\Psi_{ii}\sum_{k\neq
i}\Psi_{kk}=\Psi_{ii}\ ,\end{split}$ (S7)

where, from the second to the third line we have used the fact that the first
term is subleading and the ansatz in Eq.(S4).

##  3 Fluctuations of state-observable correlations

At the level of single matrix elements, the product of the initial state and
the observable has large fluctuations in the system size. In fact

$\Psi_{ij}A_{ij}\simeq\overline{\Psi_{ij}A_{ji}}+\sqrt{\overline{|\Psi_{ij}|^{2}}}\,\sqrt{\overline{|A_{ij}|^{2}}}\,\xi_{ij}$
(S8)

with $\xi_{ij}$ some random variable with average zero and fluctuations order
one. Here, the amplitude of the fluctuations is larger than the average:

$\sqrt{\overline{|\Psi_{ij}|^{2}}\,\overline{|A_{ij}|^{2}}}\gg\overline{\Psi_{ij}A_{ji}}\
,$

since $\sqrt{\overline{|\Psi_{ij}|^{2}}\,\overline{|A_{ij}|^{2}}}\simeq
e^{-3S/2}$ and $\overline{\Psi_{ij}A_{ji}}\sim e^{-2S}$. However, when
computing physical observables, one has to sum over many indices, and, due to
the presence of randomness, the fluctuations become negligible. This is
analogous to what happens to high-order products of matrix elements in
standard ETH, which also possess large fluctuations [15]. These fluctuations
contribute, at least for finite systems sizes, to $\langle A(t)\rangle\langle
A(-t)\rangle$, A detailed understanding of their influence on the dynamics is
left to future work.

##  4 A bound on high frequency tail of $g(\omega)$

As a starting point we introduce

$Z(\beta)=\langle\psi|e^{-\beta H}|\psi\rangle,$ (S9)

generalizing $Z(0)$ defined in (6).

To constrain $g_{\psi,A}$ we use the approach similar to one used [65], which
bounds on high-frequency tail of $f_{A}$,

$|f_{A}(\tilde{e},\omega)|\leq
O\left(e^{-(\tilde{\beta}/4+\beta^{*})\omega}\right),\qquad\omega\rightarrow\infty.$
(S10)

where $\beta^{*}$ is an $O(1)$ constant defined by local model parameters.
Here $O\left(\dots\right)$ means that possible pre-exponential
$\omega$-dependent factors are ignored. Finally, temperature $\tilde{\beta}$
is associated with energy density $\tilde{e}$,
$S^{\prime}(N\tilde{e})=\tilde{\beta}$. We now consider

$C_{\Psi}^{\beta}(t)\equiv{\langle\psi|e^{-\beta H}A(t)|\psi\rangle\over
Z(\beta)}.$ (S11)

We now use the following inequality

$|\langle\psi|A|\psi^{\prime}\rangle|\leq|A||\psi||\psi^{\prime}|,$ (S12)

where $|A|$ is an infinity norm of the operator $A$, meaning the largest (by
absolute value) eigenvalue when $A$ is hermitian, or largest singular value
when $A$ is not hermitian. Taking $|\psi\rangle=\ket{\psi}$ and
$\bra{\psi^{\prime}}=\bra{\psi}e^{-\beta H}$ we arrive at

$\left|\int
d\omega\,f_{A}(\tilde{e},\omega)g_{\psi,A}(\tilde{e},\omega)e^{\omega(it-\beta/2)}\right|\leq|A(t)|{Z^{1/2}(2\beta)Z^{1/2}(0)\over
Z(\beta)}.$ (S13)

Here $\tilde{e}$ is the energy density where the main contribution to the
integral in (S9) comes from, $\tilde{e}=-\partial/\partial\beta\ln
Z(\beta)/N$. We can now redefine $t\rightarrow t-i\beta/2$,

$\left|\int
d\omega\,f_{A}(\tilde{e},\omega)g_{\psi,A}(\tilde{e},\omega)e^{i\omega
t}\right|\leq|A(t-i\beta/2)|{Z^{1/2}(2\beta)Z^{1/2}(0)\over Z(\beta)}.$ (S14)

The LHS is an even function of $t$, while $|A(t)|$ is analytic within the
strip $|\,\mathrm{Im}(t)|\leq\beta^{*}$ [65]. Thus the RHS of (S14) is
analytic inside the strip
$-\beta^{*}+\beta/2\leq\,\mathrm{Im}(t)\leq\beta^{*}+\beta/2$. Because the LHS
is even, it has to be analytic inside a wider strip
$-\beta^{*}-\beta/2\leq\,\mathrm{Im}(t)\leq\beta^{*}+\beta/2$. For the
integral over $\omega$ to converge, taking into account the bound (S10) we
find

$|g_{\psi,A}(\tilde{e},\omega)|\leq
O\left(e^{-(\beta/2-\tilde{\beta}/4)\omega}\right){Z^{1/2}(2\beta)Z^{1/2}(0)\over
Z(\beta)},\qquad\omega\rightarrow\infty.$ (S15)

Here $\beta$ is a free parameter, it determines the “saddle point” (mean
energy density) $\tilde{e}(\beta)$, where the integral in (S9) is saturated,
which in turn defines $\tilde{\beta}$. Parameter $\beta^{*}$, which
characterizes the model does not appear in (S15). When $\beta=0$, mean energy
density $\tilde{e}=e_{0}$, and $\tilde{\beta}$ is the effective temperature of
state $\psi$.

The logic behind free parameter $\beta$ appearing in bound (S15) is exactly as
in [65], this is a parameter to optimize over, to find the best possible
bound. For the large deviation states (9), the factor ${Z^{1/2}(2\beta)\over
Z(\beta)}$ grows extensively, $\ln({Z^{1/2}(2\beta)\over Z(\beta)})\sim
N\,O(\beta^{2})$, not leading to a meaningful bound in the thermodynamic
limit.

For the states with energy of order one, e.g. discussed in the next section,
effective energy density $\tilde{e}=e_{0}$ is $\beta$-independent,
$\tilde{\beta}=\beta_{0}$, at least so for the parameter $\beta$ smaller than
certain value $\beta<\lambda$, see (S20). In this regime
${Z^{1/2}(2\beta)Z^{1/2}(0)/Z(\beta)}$ is of order one, and
$g_{\psi,A}(\tilde{e},\omega)$ for large $\omega$ decays exponentially,
bounded by $e^{-(\lambda/2-\beta_{0}/4)\omega}$, provided $\beta_{0}\leq
2\lambda$.

##  5 Bipartite energy eigenstates

Consider a tensor product Hilbert space
$\mathcal{H}=\mathcal{H}_{A}\otimes\mathcal{H}_{B}$ and a Hamiltonian
$\hat{H}=\hat{H}_{A}+\hat{H}_{B}+\hat{H}_{AB}$, describing an interaction of a
(sub)system $A$ with a “bath” $B$. A particular example to keep in mind is a
spin-chain split into two parts, interacting through a local term $H_{AB}$. We
start with a pair $\ket{a}$, $\ket{b}$ of the eigenstates of $\hat{H}_{A}$,
$\hat{H}_{B}$ respectively and consider a state $|\psi\rangle=|ab\rangle$.
Decomposition of this state in the eigenbasis $\ket{E_{i}}$ of $\hat{H}$
defines the projector (6),

$\Psi^{(ab)}_{ij}=\langle E_{i}|ab\rangle\langle
ab|E_{i}\rangle=c^{(ab)}_{i}c^{(ab)*}_{j}\ .$ (S16)

These initial states have been the focus of a great attention lately [66; 27;
29; 30; 31; 32; 33; 34], , especially in relation to bipartite entanglement.
The statistical properties of (S16) can be formulated in terms of the so-
called _Ergodic Bipartition_ (EB) ansatz [34], which postulates that on
average

$\overline{\Psi_{ii}^{(ab)}}=\overline{|c^{(ab)}_{i}|^{2}}=e^{-S(E_{a}+E_{b})}F(E_{i}-E_{a}-E_{b})\
,$ (S17)

where $F(E_{i}-E_{a}-E_{b})$ is a narrowly-peaked function around $E_{i}\simeq
E_{a}+E_{b}$. More accurately instead of $E_{a}+E_{b}$ in the expression above
one should use mean energy $E_{0}$ of $|ab\rangle$, as we do below. More
detailed properties of $F(x)$ for 1D systems with local interactions, the
Lorentzian shape at small $x$ and exponential suppression at large $x$,

$\displaystyle
F(x)\propto\left\\{\begin{array}[]{lr}(x^{2}+\Delta^{2})^{-1},&x\ll\sigma,\\\
e^{-|x|\lambda},&x\gg\sigma,\end{array}\right.$ (S20)

where $\sigma,\Delta,\lambda$ are model-dependent local (finite in the
thermodynamic limit) parameters, were outlined in [32].

From the definition (S16) it is clear the ergodic bipartition ansatz (S17) is
the diagonal part of the out-of-equilibrium ETH ansatz (6) applied to a
particular initial states. We now discuss how our approach encompasses the
properties of (S17). Starting from (7), using saddle point approximation we
find,

$\displaystyle\overline{\Psi_{ii}}=\frac{e^{-\Phi(E_{i})}}{Z}\approx
e^{-S(E_{0})}\frac{e^{-\Phi(E_{i})+\Phi(E_{0})}}{\sqrt{2\pi}\Delta_{E_{0}}},$
(S21)

where $E_{0}$ is the mean energy of state $|ab\rangle$,
$E_{0}=E_{a}+E_{b}+\Delta_{ab}$. Here $\Delta_{ab}=\langle
ab|\hat{H}_{AB}|ab\rangle$ and

$\Delta_{E_{0}}^{2}=\langle
ab|\hat{H}_{AB}^{2}|ab\rangle-(\Delta_{ab})^{2}=\mathcal{O}(1)$

are of order one, i.e. remain finite in the thermodynamic limit, while $E_{0}$
and $E_{a}+E_{b}$ are extensive. From here follows

$F(E_{i}-E_{0})\approx\frac{e^{-\Phi(E_{i})+\Phi(E_{0})}}{\sqrt{2\pi}\Delta_{E_{0}}}$
(S22)

which is sharply peaked around mean energy $E_{i}\approx E_{0}$ with the
variance of order one.

Integrating (S21) over $dE_{b}$ with the weight $e^{S_{B}(E_{b})}$ and the
constraint $E_{0}\approx E_{a}+E_{b}$ readily gives diagonal approximation to
reduced density matrix [27]

${\rm Tr}_{B}(|E_{i}\rangle\langle E_{i})\approx\int
dE_{a}\,e^{S_{A}(E_{a})}e^{-S(E_{i})+S_{B}(E_{i}-E_{a})}|a\rangle\langle a|,$
(S23)

from where von Neumann and Renyi entropy follow via saddle point
approximation. Justification of the diagonal approximation to evaluate entropy
was recently addressed in [34] by considering the statistical properties of

$c_{i}^{(ab)}c^{(a^{\prime}b)*}_{i}c_{i}^{(a^{\prime}b^{\prime})}c_{i}^{(ab^{\prime})*},$
(S24)

which can be understood extending the ansatz to four different states
$\ket{\psi_{1}}=\ket{ab}$, $\ket{\psi_{2}}=\ket{a^{\prime}b}$,
$\ket{\psi_{3}}=\ket{a^{\prime}b^{\prime}}$,
$\ket{\psi_{4}}=\ket{ab^{\prime}}$, with $i_{1}=i_{2}=i_{3}=i_{4}$.

To illustrate the behavior of $\Phi(E_{i})$ we consider, c.f. (S20),

$\ln\overline{\Psi_{ii}}\simeq-S(E_{0})-{1\over 2}\ln(2\pi\Delta
E^{2})-{(E_{i}-E_{0})^{2}\over 2\Delta^{2}_{E_{0}}}+\dots$ (S25)

and note that for $E-E_{0}$ of order one, higher-order corrections to (S25)
are also of order one and can not be neglected [32]. It corresponds to the
first two terms of the Taylor expansion of (S20) in $(E_{i}-E_{0})^{2}$. We
plot $\phi(E_{i})\equiv\ln\overline{\Psi_{ii}}+S(E_{0})+{1\over
2}\ln(2\pi\Delta E^{2})$ as a function of $E_{i}-E_{0}$ numerically in the
left panel of Fig. S1 for the tilted field Ising model (21) of size
$L=2L^{\prime}+1$ with the parameters $J=-1,w=1.05,h=0.4$. The eigenstates
$|a\rangle,|b\rangle$ are chosen to be the ground state and the most excited
states of the subsystems of size $L^{\prime}$ and $L^{\prime}+1$
correspondingly. With this choice of $\psi_{0}$ the total energy $E_{0}$ is
very close to the middle of the spectrum. To obtain $\overline{\Psi_{ii}}$ we
use smoothed average (24) with $\tau=2$. The entropy

$\displaystyle
S(E)=\ln\Omega(E_{0}),\qquad\Omega(E)={\kappa\,L!\over(L/2-\kappa
E)!(L/2+\kappa E)!}$ (S26)

is evaluated using binomial analytic approximation for the density of states
of the titled filed Ising model, see Appendix A of [27], with $\kappa={1\over
2}\sqrt{J^{2}+w^{2}+h^{2}-1/L}$. Mean energy $E_{0}$ and energy variance
$\Delta^{2}_{E0}$ are evaluated numerically for each $L$. The plot in S1 shows
good collapse of $\phi(E_{i}-E_{0})$ for different values of $L$, and is well
described by ${(E_{i}-E_{0})^{2}\over 2\Delta^{2}_{E_{0}}}$ for small
$|E_{i}-E_{0}|$.

The out-of-equilibrium ETH anzats predicts the off-diagonal matrix elements

$\varphi(\omega)=\ln\overline{\Psi_{ij}^{2}}-2\ln\overline{\Psi_{ee}}\simeq-{\omega^{2}\over
4\Delta^{2}_{E_{0}}},\quad|\omega|\lesssim\Delta_{E_{0}},$ (S27)

plotted in the right panel of Fig. S1. Here $\overline{\Psi_{ee}}$ denotes the
average $\overline{\Psi_{ii}}$ over a narrow shell around $E_{i}=E_{0}$ of
size $\Delta_{E_{0}}$. Similarly, $\Psi_{ij}^{2}$ in (S27) is averaged only
over pairs $E_{i},E_{j}$ satisfying
$|(E_{i}+E_{j})/2-E_{0}|\leq\Delta_{E_{0}}$.

Figure S1: Left panel: averaged diagonal matrix elements (S21), with the
entropy subtracted $\phi(E_{i})\equiv\ln\overline{\Psi_{ii}}+S(E_{0})+{1\over
2}\ln(2\pi\Delta E^{2})$, for different system sizes, superimposed with
${-(E_{i}-E_{0})^{2}/(2\Delta^{2}_{E_{0}})}$ with the value of
$\Delta^{2}_{E_{0}}$ for $L=15$ (green dashed line). Blue points show raw data
(un-averaged value of $\ln\Psi_{ii}$, with the entropy subtracted) for $L=15$.
Right panel: average of the off-diagonal matrix elements (S27) for different
system sizes, superimposed with ${-\omega^{2}/(4\Delta^{2}_{E_{0}})}$ with the
value of $\Delta^{2}_{E_{0}}$ for $L=13$ (green dashed line).

Next, we study cross-correlations of $\psi$ with different operators $A$,
sitting at the L(eft) edge, R(ight) edge, or the M(iddle) site of the chain,

$A=\sigma^{z,x}_{L}=\sigma^{z,x}_{site=1}\ ,\quad
A=\sigma^{z,x}_{M}=\sigma^{z,x}_{site=L^{\prime}+1}\ ,\quad
A=\sigma^{z,x}_{R}=\sigma^{z,x}_{site=L}\ .$ (S28)

We plot $f_{A}(\omega)g_{\psi,A}(\omega)$ for different $A$ in Fig. S2. In the
second row of Fig. S2, we display the scaling with the system size of
Eqs.(10), finding a good agreement with the predictions at zero and finite
frequency ($\omega=2.1$ shown in the Figure).

Figure S2: First row: from left to right, plots of
$f_{A}(\omega)g_{\psi,A}(\omega)$ for
$A=\sigma^{z}_{L},\sigma^{z}_{M},\sigma^{z}_{R}$ correspondingly. Second row:
scaling of Eqs.(10) contrasted with $(\dim\mathcal{H})^{-1}$ (dashed red line)
and $(\dim\mathcal{H})^{-2}$ (dashed black line).

Finally, in Fig. S3 we illustrate approximate independence of
$g_{\psi,A}(\omega)$ at high frequencies, on the choice of the observable $A$.

Figure S3: In the absence of translational invariance, the function $g$ is
independent of the direction of the observable, but not on the site.


# Bayesian inference of strange star equation of state using the GW170817 and
GW190425 data

Zhiqiang Miao  Department of Astronomy, Xiamen University, Xiamen, Fujian
361005, China<EMAIL_ADDRESS>Jin-Liang Jiang  Key Laboratory of Dark Matter
and Space Astronomy, Purple Mountain Observatory, Chinese Academy of Sciences,
Nanjing 210023, China  School of Astronomy and Space Science, University of
Science and Technology of China, Hefei, Anhui 230026, China  Ang Li
Department of Astronomy, Xiamen University, Xiamen, Fujian 361005, China  Lie-
Wen Chen  School of Physics and Astronomy, Shanghai Key Laboratory for
Particle Physics and Cosmology, and Key Laboratory for Particle Astrophysics
and Cosmology (MOE), Shanghai Jiao Tong University, Shanghai 200240, China

###### Abstract

The observations of compact star inspirals from LIGO/Virgo provide a valuable
tool to study the highly uncertain equation of state (EOS) of dense matter at
the densities in which the compact stars reside. It is not clear whether the
merging stars are neutron stars or quark stars containing self-bound quark
matter. In this work, we explore the allowed bag-model-like EOSs by assuming
the merging stars are strange quark stars (SQSs) from a Bayesian analysis
employing the tidal deformability observational data of the GW170817 and
GW190425 binary mergers. We consider two extreme states of strange quark
matter, either in nonsuperfluid or color-flavor locked (CFL) and find the
results in these two cases essentially reconcile. In particular, our results
indicate that the sound speed in the SQS matter is approximately a constant
close to the conformal limit of $c/\sqrt{3}$. The universal relations between
the mass, the tidal deformability and the compactness are provided for the
SQSs. The most probable values of the maximum mass are found to be $M_{\rm
TOV}=2.10_{-0.12}^{+0.12}~{}(2.15_{-0.14}^{+0.16})\,M_{\odot}$ for normal
(CFL) SQSs at a $90\%$ confidence level. The corresponding radius and tidal
deformability for a $1.4\,M_{\odot}$ star are $R_{\rm
1.4}=11.50_{-0.55}^{+0.52}~{}({11.42}_{-0.44}^{+0.52})~{}\rm km$ and
$\Lambda_{1.4}={650}_{-190}^{+230}~{}({630}_{-150}^{+220})$, respectively. We
also investigate the possibility of GW190814’s secondary component $m_{2}$ of
mass $2.59_{-0.09}^{+0.08}\,M_{\odot}$ being an SQS, and find that it could be
a CFL SQS with the pairing gap $\Delta$ larger than $244~{}\rm MeV$ and the
effective bag parameter $B_{\rm eff}^{1/4}$ in the range of $170$ to $192$
MeV, at a $90\%$ confidence level.

Compact objects (288); Gravitational waves (678); High energy
astrophysics(739)

††software: Bilby (Ashton et al., 2019, v0.5.5,
https://git.ligo.org/lscsoft/bilby/), PyMultiNest (Buchner, 2016, v2.6,
https://github.com/JohannesBuchner/PyMultiNest).

##  1 Introduction

It is generally believed that the degree of freedom of dense matter is
hadronic around the nuclear saturation density, $n_{0}\approx 0.16$ fm-3, and
the color-flavor locked (CFL) state is expected to be the ground state of
three-flavor quark matter at asymptotic densities (Alford et al., 2008). The
phase states of cold QCD matter at intermediate densities ($\sim
1-10~{}n_{0}$) are unfortunately unknown (see discussions in, e.g., Kurkela et
al., 2010, 2014; Gorda et al., 2018). One key point is still unclear, i.e.,
whether compact stars are gravity-bound neutron stars (NSs) or indeed self-
bound quark stars (QSs)? After decades of speculation (Bodmer, 1971; Witten,
1984), QSs still serve as viable alternative physical model for compact stars
(Baym et al., 1985; Glendenning, 1990; Weber, 2005; Li et al., 2016; Zhou et
al., 2018; Bombaci et al., 2021; Cao et al., 2020; Traversi et al., 2021;
Sedaghat et al., 2021).

The recently observed binary star merger events (Abbott et al., 2017, 2018,
2019, 2020a, 2020b) have greatly promoted the study of the equation of state
(EOS) of dense stellar matter (pressure $p$ as a function of the density $n$,
or energy density $e$), restricting its stiffness or the degree of freedom of
dense matter in the density regime achieved inside compact stars (possibly up
to $\approx 8-10\,n_{0}$). After the release of the gravitational-wave
observational data from LIGO/Virgo, in-depth studies have been performed (see
recent reviews by Baiotti, 2019; Li et al., 2020; Chatziioannou, 2020;
Dietrich et al., 2021). However, the available mass and radius constraints
derived from the two merger events, GW170817 (Abbott et al., 2017) and
GW190425 (Abbott et al., 2020a), commonly assume NS EOSs, with or without
quark deconfinement phase transition in their cores (e.g., Bauswein et al.,
2017; Abbott et al., 2018; Most et al., 2018; Fattoyev et al., 2018; Abbott et
al., 2019; Fasano et al., 2019; Miao et al., 2020; Li et al., 2021b). Note
that the quark matter in NS cores is not self-bound (Miao et al., 2020; Li et
al., 2021b; Xia et al., 2021) and it only appears through a phase transition
from hadronic matter at high densities. Since the results assuming NSs cannot
be used to study QS properties which have a sharp surface due to the self-
bound feature of their EOS and that the parameterized NS EOS cannot catch this
behavior (see also discussions in Zhang (2020)), an updated parameter space
for QS EOSs based on the LIGO/Virgo observations is desirable.

As an endeavor in this direction, in the present study, we perform a Bayesian
analysis of the GW170817 and GW190425 data by assuming the merging stars are
strange quark stars (SQSs). We restrict ourselves to only two extreme cases of
normal nonsuperfluid strange quark matter (SQM) and CFL quark matter since it
is not known how far the CFL phase extends towards lower density.

This paper is organized as follows. In Section 2, we introduce the EOS model
for describing the SQSs. Section 3 presents the observational constraints
employed and the Bayesian analysis method that we apply. Our results are
presented in Section 4 and summarized in Section 5.

##  2 Normal and CFL SQS Models

SQSs exist when one applies the Bodmer-Witten hypothesis, namely self-bound
SQM is the physical nature of all compact stars (Bodmer, 1971; Witten,
1984).111Note that there are also discussions in the literature suggesting
that non-SQM can have a lower bulk energy per baryon than normal nuclei and
SQM (e.g., Holdom et al., 2018; Cao et al., 2020; Zhang & Mann, 2021). Various
approaches have been attempted for a model description of self-bound quark
matter (Xia et al., 2021) since it is not attainable directly by solving QCD
(Kurkela et al., 2010, 2014; Gorda et al., 2018). Among them, the MIT bag
model is the most widely used (see some of the latest studies, e.g., Li et
al., 2020; Roupas et al., 2021; Lopes et al., 2021; Joshi et al., 2021).
Others include the Nambu-Jona-Lasinio (NJL) model (Buballa, 2005), the
density-dependent quark masses (Li et al., 2010, 2011), a confining quark
matter model (Dey et al., 1998; Chu & Chen, 2017), an interacting quark-matter
model (Zhang & Mann, 2021) and the vector interaction enhanced-bag model
(Klähn & Fischer, 2015), as an incomplete list.

The SQM is composed of up ($u$), down ($d$) and strange ($s$) quarks with the
charge neutrality maintained by the inclusion of electrons. In our
calculations below, we shall adopt the bag model with and without
superfluidity. For ease of discussion later, we first briefly recall some of
the basic formulae and definitions. The expressions for the grand canonical
potential per unit volume in the MIT bag model is written as (Alcock et al.,
1986; Haensel et al., 1986; Alford et al., 2005):

$\Omega=\sum_{{i=u,\ d,\ s,\
e^{-}}}\Omega_{i}^{0}+\frac{3(1-a_{4})}{4\pi^{2}}\mu^{4}+B_{\mathrm{eff}}+\frac{3m_{s}^{4}-48\Delta^{2}\mu^{2}}{16\pi^{2}}\
,$ (1)

where $\Omega_{i}^{0}$ is the grand canonical potential for particle type $i$
described as ideal Fermi gas. The baryon chemical potential is $\mu_{\rm
B}=\mu_{u}+\mu_{d}+\mu_{s}$ and $\mu=\mu_{\rm B}/3$ is the average quark
chemical potential. The total baryon number density can be expressed as
$n=(n_{u}+n_{d}+n_{s})/3$ with $n_{i}=-(\partial\Omega/\partial\mu_{i})_{\rm
V}$. The effective bag constant $B_{\rm eff}$, accounting for the QCD vacuum’s
contributions, is usually regarded as a phenomenological parameter. The
parameter $a_{4}$ characterizes the QCD corrections due to gluon-mediated
interactions between quarks (Fraga et al., 2001; Alford et al., 2005;
Bhattacharyya et al., 2016; Li et al., 2017), with $a_{4}=1$ corresponding to
no QCD corrections (Fermi gas approximation). $\Delta$ is the CFL pairing gap
of an order of tens to 100 MeV (Alford et al., 1998; Rapp et al., 1998).
$\Delta=0$ corresponds to the nonsuperfluid case. There may be other
superfluid phases in the intermediate densities (Alford & Sedrakian, 2017;
Zhang & Mann, 2021), e.g., the two-flavor color-superconducting phase, a
gapless CFL phase, which are not considered here and will be studied in
future. Finite quark mass was shown to have a trivial influence on the results
(Li et al., 2020; Zhou et al., 2018), and we neglect the masses of the up and
down quarks while fix the strange quark mass as $m_{s}=100\,{\rm MeV}$. The
EOS of SQM, which is characterized by the independent parameters $B_{\rm eff}$
and $a_{4}$ ($\Delta$ as well if quark superfluity is included), can be
calculated using basic thermodynamic relations, and the stars’ global
properties (mass $M$, radius $R$, tidal deformability $\Lambda$) can then be
obtained (see details in e.g., Zhou et al., 2018; Zhu et al., 2018; Miao et
al., 2020).

##  3 Observational Constraints and Bayesian Analysis

###  3.1 Bayesian Analysis

Considering the maximal spin observed in Galactic pulsars, Abbott et al.
(2019) inferred that the 90% credible intervals for the component masses of
the GW170817 event lie between $1.16$ and $1.60\,M_{\odot}$, and Abbott et al.
(2020a) reported that the corresponding component masses range from $1.46$ to
$1.87\,M_{\odot}$ for the GW190425 event. We treat both events as SQS-SQS
mergers. By exploiting their tidal deformability measurements, we use Bayesian
inference to determine posteriors on the EOS parameter spaces of normal and
CFL SQSs.

Assuming that the noise in the LIGO/Virgo detectors is Gaussian and
stationary, the likelihood of a gravitational-wave event used to perform
Bayesian inference is often expressed as

$L(d|\vec{\theta}_{\rm GW})\\!\propto\\!{\rm
Exp}\\!\left(-2\int\frac{|d(f)-h(\vec{\theta}_{\rm GW},f)|^{2}}{S_{\rm
n}(f)}df\right),$ (2)

where $S_{\rm n}(f)$, $d(f)$, and $h(\vec{\theta}_{\rm GW},f)$, respectively,
denote the power spectral density (PSD), the frequency domain data, and the
frequency domain waveform generated using parameter set $\vec{\theta}_{\rm
GW}$. The tidal deformability $\Lambda$ encoded in the gravitational-wave
strain data can be mapped from the mass through the EOS. Thus the EOS
parameters together with component masses can be incorporated to construct the
gravitational-wave parameters $\vec{\theta}_{\rm GW}$ (e.g., Jiang et al.,
2020; Tang et al., 2020). We take the publicly available strain
data222https://www.gw-openscience.org/eventapi and PSDs333The PSD of GW170817
can be found at https://doi.org/10.7935/KSX7-QQ51, while for GW190425, the PSD
can be found in the parameter estimation sample release at
https://dcc.ligo.org/LIGO-P2000026/public. of GW170817 and GW190425, together
with the waveform model IMRPhenomD_NRTidal (Dietrich et al., 2017) to do the
analysis. We follow all the data analysis details of GW170817 and GW190425
described in (Abbott et al., 2019) and (Abbott et al., 2020a), respectively,
except that we only consider the low-spin case and do not consider the
calibration error which causes a minor difference in our results. When
combining the two events with the analysis, we encounter a problem whereby
many sampling parameters make it computationally expensive to convergence in
the Nest sampling algorithm. So we take the interpolated likelihood tables of
Hernandez Vivanco et al. (2020) into the analysis, which marginalized over all
the other parameters except masses and tidal deformabilities. We find that
using the interpolated likelihood tables is well consistent with those of
using gravitational-wave data in analyzing a single event. Then by employing
the python-based Bilby (Ashton et al., 2019) and Pymultinest (Buchner, 2016)
packages, we simultaneously infer the gravitational-wave parameters and the
EOS parameters.

###  3.2 Priors and Constraints

To improve the converging rate of the nest sampling, we marginalize the
coalescence phase parameter in the likelihood and fix the source’s sky
location determined by the electromagnetic observations (Abbott et al., 2017;
Levan et al., 2017). As for the priors of the other parameters in
$\vec{\theta}_{\rm GW}$, we take a similar choice as presented in Tang et al.
(2020). For the parameters $\vec{\theta}_{\rm EOS}$ that construct the EOS of
SQSs, following our previous studies (Zhou et al., 2018), we assign wide
boundaries to them as $B_{\rm eff}^{1/4}\in[125,150]\,{\rm MeV}$,
$\Delta\in[0,100]\,{\rm MeV}$, and $a_{4}\in[0.4,1]$ as theoretically
estimated, with which both uniform and logarithmic uniform distributions are
investigated. For technical reasons, the lower bound of the logarithmic
uniform distribution cannot be zero; thus, we set a reasonable lower bound
$0.1$ for $\Delta$ in the logarithmic uniform case.

As mentioned earlier, we assume all compact stars are self-bound SQSs and that
they are composed of charge-neutral bulk SQM. Consequently, two stability
constraints should be adopted: first, the energy per baryon for non-SQM should
satisfy $(E/A)_{\rm ud}\geq 934\,{\rm MeV}$ to guarantee the observed
stability of atomic nuclei; second, $(E/A)_{\rm uds}\leq 930\,{\rm MeV}$ is
required, according to the hypothesis that SQM is absolutely stable (Bodmer,
1971; Witten, 1984). see also Li et al. (2010, 2011). In addition, the
causality condition for the SQM EOS is guaranteed for all the SQS
calculations. The mass measurement of massive pulsars establishes a lower
bound on the maximum mass of SQSs. Only the EOSs that support a $M_{\rm TOV}$
larger than this lower bound can pass this constraint, while others will be
rejected. We adopt the largest mass measured through Shapiro delay, $M=2.08\pm
0.07\,M_{\odot}$ (68% confidence level) of MSP J0740+6620 (Cromartie et al.,
2020; Fonseca et al., 2021), to place the $M_{\rm TOV}$ constraint. A
Gaussian-like likelihood is used to encapsulate the mass measurement of the
$2.08\,\,M_{\odot}$ pulsar. More details on the implementation of the $M_{\rm
TOV}$ limit can be found in our previous study (Li et al., 2021b).

##  4 Results and Discussion

Figure 1: Posterior distributions ($90\%$ confidence level) of the pressure
vs. chemical potential relation (left panel) and the EOS (right panel) for
normal and CFL QSs from the GW170817 and GW190425 data, together with the pQCD
result from Kurkela et al. (2010) in the left panel and the NS result based on
GW170817 from LIGO/Virgo (Abbott et al., 2018) in the right panel. The central
densities of normal and CFL maximum-mass stars are marked in the left panel
with blue and red symbols. An additional upper horizontal axis are indicated
in the left panel, showing the corresponding energy density within the most
favored EOS parameter sets in the QS analysis.

###  4.1 EOS and the Pressure vs. Chemical Potential Relation

Table 1: Most probable intervals of the EOS parameters ($90\%$ confidence
level) constrained by the joint GW170817+GW190425 analysis for the two priors:
Uniform ($U$) and logarithmic uniform (${\rm log}U$) distributions.

Parameters | Prior Type |  | Joint Analysis
---|---|---|---
$B_{\rm eff}^{1/4}/\rm MeV$ |  $U(125,150$) | Normal | $137.6_{-4.4}^{+4.8}$
|  | CFL | $144.3_{-7.9}^{+5.0}$
|  ${\rm log}U(125,150$) | Normal | $137.3_{-4.2}^{+4.7}$
|  | CFL | $138.5_{-5.3}^{+8.3}$
$a_{4}$ |  $U(0.4,1$) | Normal | $0.70_{-0.12}^{+0.16}$
|  | CFL | $0.74_{-0.20}^{+0.22}$
|  ${\rm log}U(0.4,1$) | Normal | $0.69_{-0.11}^{+0.15}$
|  | CFL | $0.70_{-0.13}^{+0.18}$
$\Delta/\rm MeV$ |  $U(0,100$) | CFL | ${67.6}_{-54.3}^{+28.1}$
|  ${\rm log}U(0.1,100$) | CFL | $6.6_{-6.4}^{+73.7}$

The most probable values of the EOS parameters and their $90\%$ confidence
boundaries, constrained jointly by the GW170817 and GW190425 data, are
reported in Table 1. The corresponding results of the EOS (including the sound
velocity) and the SQS properties are shown in Figs. 1-3.

From Table 1, we see that two EOS parameters ($B_{\rm eff}^{1/4}$ and $a_{4}$)
are relatively well constrained: $B_{\rm eff}^{1/4}\approx 137~{}\rm MeV$
($144~{}\rm MeV$) and $a_{4}\approx 0.70$ ($0.74$) for the normal (CFL)
matter, which depend weakly on the prior choice. It is worth mentioning that
the “standard value” of the bag parameter is $(144\,\rm MeV)^{4}$ ($56\,\rm
MeV/fm^{3}$) from reproducing the mass spectrum of light hadrons and heavy
mesons (Degrand et al., 1975; Haxton & Heller, 1980), and the estimated value
from lattice calculations at zero chemical potential is $\sim(212\,\rm
MeV)^{4}$ ($\approx 262\,\rm MeV/fm^{3}$) (Braun-Munzinger & Stachel, 1996).
The inferred $a_{4}$ values here are close to the value suggested in Fraga et
al. (2001). However, the color superconductivity gap $\Delta$ is poorly
constrained by the observed global star properties. Nevertheless, the $\Delta$
could be constrained through future tidal deformability measurement of massive
QSs (close to the maximum mass) as discussed in Li et al. (2020).

In Fig. 1, we report the posterior distributions of the SQS EOS as well as the
pressure vs. chemical potential ($p$-$\mu_{B}$) relation. Also shown in the
$p$-$\mu_{\rm B}$ plot is the results from perturbative QCD (Kurkela et al.,
2010), which is only applicable above $\mu_{\rm B}\sim 2.6~{}\rm GeV$. It is
seen that the observational data can constrain the EOS effectively in the
nonperturbative realm where the properties of compact stars are relevant. It
is common in the literature to analyze the gravitational-wave data by assuming
all compact stars are NSs; the resulting EOS from LIGO/Virgo (Abbott et al.,
2018) is included in the right panel of Fig. 1. It is seen that SQS EOSs
become stiff earlier than the NS ones (Abbott et al., 2018) at low densities,
indicating the sharp surface of SQSs is of a relatively low density (Li et
al., 2017), as low as $\sim 0.22~{}\rm fm^{-3}$ ($\approx 1.4n_{0}$), to the
$90\%$ confidence level. However, at high densities, SQS EOSs are softer than
the NS ones (Abbott et al., 2018; Li et al., 2021b).

Figure 2: Posterior distributions of the mass-radius relation (left), the
mass-compactness relation (middle) and $\Lambda\beta^{6}$ (right) as functions
of the mass for normal and CFL SQSs from the joint GW170817 and GW190425
analysis, to the 90% confidence level. Also shown in the left panel are the
mass and radius analysis of PSR J0030+0451 and MSP J0740+6620 based on NS EOSs
from the Neutron Star Interior Composition Explorer (Miller et al., 2019;
Riley et al., 2019, 2021; Miller et al., 2021).  Figure 3: Posterior
distributions of the squared sound speed $c^{2}_{\rm s}$ (in units of $c^{2}$)
as a function of the density, for normal and CFL SQSs from the joint GW170817
and GW190425 analysis, to the 90% confidence level. Also shown is the
conformal limit of $c/3$.

###  4.2 SQS Mass, Radius and Tidal deformability

Generally, in the bag model, the stability of quark matter is dominated by the
vacuum term ($B_{\rm eff}$) and the perturbative interaction term ($a_{4}$)
(Farhi & Jaffe, 1984), with the quark pairing $\Delta$ helping to lower the
energy (Alford et al., 2001; Rajagopal & Wilczek, 2001; Lugones & Horvath,
2002, 2003; Zhou et al., 2018). As a consequence, the introduction of pairing
(as an extra degree of freedom) allows wider parameter spaces for the model
parameters ($B_{\rm eff},a_{4}$), the EOS, and SQS properties, as one can
observe in Table 1 and Figs. 1-3. In particular, the CFL matter (described in
a three-parameter model) allows more parameter space for SQS EOSs than the
normal SQM matter (described in a two-parameter model) (Li et al., 2020).
However, their predictions for SQS mass, radius, and tidal deformability agree
considerably well with each other, as demonstrated in Fig. 2. As a result, the
SQS’s mass and tidal deformability can be described nearly universally as
functions of the mass/radius (or equivalently the compactness $\beta\equiv
M/R$): $M(\beta)=13.07\times\beta^{1.32}$ and $\Lambda\beta^{6}(M)=0.02\times
M^{2.25}/(\exp(0.53\times M^{1.88})-1)$, with their coefficient of
determination being $r^{2}=0.993$ and $r^{2}=0.989$; One can also evaluate the
stars’ important properties, for example, the maximum mass is $M_{\rm
TOV}=2.10_{-0.12}^{+0.12}~{}(2.15_{-0.14}^{+0.16})\,M_{\odot}$, the radius for
a $1.4\,M_{\odot}$ star is $R_{\rm
1.4}=11.50_{-0.55}^{+0.52}~{}({11.42}_{-0.44}^{+0.52})~{}\rm km$ and the
corresponding tidal deformability is
$\Lambda_{1.4}={650}_{-190}^{+230}~{}({630}_{-150}^{+220})$ for normal (CFL)
SQSs, at a $90\%$ confidence level. The maximum mass around $2.15\,M_{\odot}$
(with an upper bound of $2.31\,M_{\odot}$) is very close to our previous
theoretical calculations (Zhou et al., 2018) with only the GW170817 constraint
considered. These results are potentially useful for the discussion of the
possible transition from NSs to QSs in the two-family scenario (Drago et al.,
2007; Bombaci et al., 2009; Drago et al., 2014; De Pietri et al., 2019; Char
et al., 2019; Drago & Pagliara, 2020), where NSs coexist with QSs and may
transit to QSs in their evolution as the central engines of short gamma-ray
bursts (Cheng & Dai, 1996; Bombaci & Datta, 2000; Li et al., 2016, 2017).

###  4.3 Sound Speed in Dense Matter

The earlier EOS stiffening in SQSs, compared to in NSs, has profound
influences on the density behavior of the squared speed of sound $c^{2}_{s}$
(in units of $c^{2}$) in the two kinds of stellar matter, which is shown in
Fig. 3. There have been many previous discussions saying that, to fulfill the
two-solar-mass lower-mass limit for compact stars, the $c_{s}$ value in NS
matter should substantially exceed the conformal limit of $c/\sqrt{3}$
(Kurkela et al., 2014; Bedaque & Steiner, 2015; Alsing et al., 2018), even
close to $\sim 0.9c$ in some studies (Tews et al., 2018), at around $5n_{0}$.
Such requirements are confirmed by both theoretical calculations (Xia et al.,
2021) and statistical analysis (Li et al., 2021b; Landry et al., 2020).
Considerable efforts are undertaken in the literature to pursue the underlying
mechanism explaining the rapid growth in $c_{s}$ (McLerran & Reddy, 2019).
However, from our present study of SQSs with a bag-model-like EOS, it is found
that $c_{s}$ is essentially a constant close to $c/\sqrt{3}$ (see also
Traversi & Char, 2020; Traversi et al., 2021)); A heavy compact star does not
necessarily demand a superconformal $c_{s}$ as long as the EOS stiffening
happens early. In fact, we generally have two kinds of EOSs in our model: the
first kind has a maximum sound speed at a low density corresponding to the
zero-pressure point, while the second kind reaches its maximum $c/\sqrt{3}$ at
asymptotic density. This is related to the uncertain superfluid phases
mentioned above, controlled in the grand canonical potential by the sign of a
combined term, $m_{s}^{2}-4\Delta^{2}$ (named coefficient $a_{2}$ in Alford et
al. (2005)). A negative $a_{2}$ corresponds to the first kind while a positive
one corresponds to the second kind (see more discussions in Alford et al.,
2005; Zhang & Mann, 2021). This causes the different $c_{s}$ behavior in the
normal and CFL cases, although they both give an approximately constant sound
speed close to the conformal limit. We mention here that the QS study from
gravitational wave data is still at an early stage with simple EOS modeling,
and it will be interesting to see how our results change if more sophisticated
quark-matter models are applied. For instance, it will be fascinating to
explore SQSs in the future by adopting the vector interaction enhanced-bag
model (Klähn & Fischer, 2015), where the effects of dynamical chiral symmetry
breaking and vector repulsion are included, based on the NJL model with vector
interaction.

###  4.4 Studies on CFL SQSs with Enlarged Parameter Spaces for Effective Bag
and Superconducting Gap

In the analysis above, the parameters of the effective bag $B_{\rm eff}$ and
CFL pairing gap $\Delta$ are varied in their theoretically estimated regions,
following our previous study (Zhou et al., 2018). However, some larger
$\Delta$ values are also used in the literature (Alford & Reddy, 2003). In
particular, much higher values for $B_{\rm eff}$ and $\Delta$ are used in
recent strange matter studies when discussing the possibility of GW190814’s
secondary component $m_{2}$ of mass $2.50-2.67\,M_{\odot}$ (Abbott et al.,
2020b) as a static CFL SQS (Roupas et al., 2021). It was found that this is
possible with the combination of a large bag parameter ($B_{\rm
eff}^{1/4}>159~{}\rm MeV$) and pairing parameter ($\Delta>200~{}\rm MeV$).

By employing the mass measurement of GW190814’s secondary component,
$2.59_{-0.09}^{+0.08}\,M_{\odot}$ (at the $90\%$ credible level), as the lower
bound on the maximum mass (instead of the MSP J0740+6620 measurement used
above), we extend the ranges of both parameters and repeat the analysis. In
particular, the upper boundary of $B_{\rm eff}^{1/4}$ ($\Delta$) is increased
from $150~{}\rm MeV$ ($100~{}\rm MeV$) to $250~{}\rm MeV$ ($500~{}\rm MeV$).
Our results suggest $\Delta>244\,{\rm MeV}$ and $170\,{\rm MeV}<B_{\rm
eff}^{1/4}<192\,{\rm MeV}$, to the 90% confidence level, in the context of
GW190814’s secondary component being a CFL SQS. In addition, we find that, for
a $2.6\,M_{\odot}$ star like GW190814’s secondary component, the radius is
$11.35\,{\rm km}<R_{2.6}<13.67\,{\rm km}$ and the corresponding tidal
deformability is $4.6<\Lambda_{2.6}<26.3$. Upcoming measurements on them, when
available, should shed light on the nature of compact objects with such large
masses.

##  5 Summary

We have performed a Bayesian analysis on the tidal deformability observations
of the gravitational-wave events GW170817 and GW190425 by assuming the merging
stars are SQSs based on bag-model-like EOSs. The results indicate that the
sound speed in the SQS matter is approximately a constant close to the
conformal limit of $c/\sqrt{3}$, in stark contrast with its rapid growth
behavior in NS matter. The SQSs are found to have a maximum mass of at most
$M_{\rm TOV}=2.31\,M_{\odot}$ at a $90\%$ confidence level; it is
$2.10_{-0.12}^{+0.12}\,M_{\odot}$ and $2.15_{-0.14}^{+0.16}\,M_{\odot}$ for
the normal and CFL matter, respectively. Several universal relations between
the observed properties of SQSs are provided. We also find that GW190814’s
secondary component could be a CFL quark star if the effective bag constant
and CFL pairing are large enough.

We are thankful to Alessandro Drago, Zheng Cao, and Toru Kojo for helpful
discussions. The work is supported by National SKA Program of China (No.
2020SKA0120300), the National Natural Science Foundation of China (grants Nos.
11873040 and 11625521), the science research grants from the China Manned
Space Project (No. CMS-CSST-2021-B11) and the Youth Innovation Fund of Xiamen
(No. 3502Z20206061). This research has made use of data and software obtained
from the Gravitational Wave Open Science Center (https://www.gw-
openscience.org), a service of LIGO Laboratory, the LIGO Scientific
Collaboration, and the Virgo Collaboration. LIGO is funded by the U.S.
National Science Foundation. Virgo is funded by the French Centre National de
Recherche Scientifique (CNRS), the Italian Istituto Nazionale della Fisica
Nucleare (INFN) and the Dutch Nikhef, with contributions by Polish and
Hungarian institutes.

## References

  * Abbott et al. (2017) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, Phys. Rev. Lett., 119, 161101. doi:10.1103/PhysRevLett.119.161101
  * Abbott et al. (2018) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018, Phys. Rev. Lett., 121, 161101. doi:10.1103/PhysRevLett.121.161101
  * Abbott et al. (2019) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, PhRvX, 9, 011001. doi:10.1103/PhysRevX.9.011001
  * Abbott et al. (2020a) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2020, ApJ, 892, L3. doi:10.3847/2041-8213/ab75f5
  * Abbott et al. (2017) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, ApJ, 848, L12. doi:10.3847/2041-8213/aa91c9
  * Abbott et al. (2020b) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020, ApJ, 896, L44. doi:10.3847/2041-8213/ab960f
  * Alcock et al. (1986) Alcock, C., Farhi, E., & Olinto, A. 1986, ApJ, 310, 261. doi:10.1086/164679
  * Alford et al. (2001) Alford, M., Rajagopal, K., Reddy, S., et al. 2001, Phys. Rev. D, 64, 074017. doi:10.1103/PhysRevD.64.074017
  * Alford et al. (1998) Alford, M., Rajagopal, K., & Wilczek, F. 1998, PhLB, 422, 247. doi:10.1016/S0370-2693(98)00051-3
  * Alford & Sedrakian (2017) Alford, M. & Sedrakian, A. 2017, Phys. Rev. Lett., 119, 161104. doi:10.1103/PhysRevLett.119.161104
  * Alford et al. (2005) Alford, M., Braby, M., Paris, M., et al. 2005, ApJ, 629, 969. doi:10.1086/430902
  * Alford & Reddy (2003) Alford, M. & Reddy, S. 2003, Phys. Rev. D, 67, 074024. doi:10.1103/PhysRevD.67.074024
  * Alford et al. (2008) Alford, M. G., Schmitt, A., Rajagopal, K., et al. 2008, Reviews of Modern Physics, 80, 1455. doi:10.1103/RevModPhys.80.1455
  * Alsing et al. (2018) Alsing, J., Silva, H. O., & Berti, E. 2018, MNRAS, 478, 1377. doi:10.1093/mnras/sty1065
  * Ashton et al. (2019) Ashton, G., Hübner, M., Lasky, P. D., et al. 2019,Bilby: Bayesian inference library, Astrophysics Source Code Library, ascl:1901.011
  * Ashton et al. (2019) Ashton, G., Hübner, M., Lasky, P. D., et al. 2019, ApJS, 241, 27. doi:10.3847/1538-4365/ab06fc
  * Baiotti (2019) Baiotti, L. 2019, Progress in Particle and Nuclear Physics, 109, 103714. doi:10.1016/j.ppnp.2019.103714
  * Bauswein et al. (2017) Bauswein, A., Just, O., Janka, H.-T., et al. 2017, ApJ, 850, L34. doi:10.3847/2041-8213/aa9994
  * Baym et al. (1985) Baym, G., Kolb, E. W., McLerran, L., et al. 1985, PhLB, 160, 181. doi:10.1016/0370-2693(85)91489-3
  * Bedaque & Steiner (2015) Bedaque, P. & Steiner, A. W. 2015, Phys. Rev. Lett., 114, 031103. doi:10.1103/PhysRevLett.114.031103
  * Bhattacharyya et al. (2016) Bhattacharyya, S., Bombaci, I., Logoteta, D., et al. 2016, MNRAS, 457, 3101. doi:10.1093/mnras/stw206
  * Bodmer (1971) Bodmer, A. R. 1971, Phys. Rev. D, 4, 1601. doi:10.1103/PhysRevD.4.1601
  * Bombaci et al. (2021) Bombaci, I., Drago, A., Logoteta, D., et al. 2021, Phys. Rev. Lett., 126, 162702. doi:10.1103/PhysRevLett.126.162702
  * Bombaci et al. (2009) Bombaci, I., Logoteta, D., Panda, P. K., et al. 2009, PhLB, 680, 448. doi:10.1016/j.physletb.2009.09.039
  * Bombaci & Datta (2000) Bombaci, I. & Datta, B. 2000, ApJ, 530, L69. doi:10.1086/312497
  * Braun-Munzinger & Stachel (1996) Braun-Munzinger, P. & Stachel, J. 1996, Nucl. Phys. A, 606, 320. doi:10.1016/0375-9474(96)00198-4
  * Buballa (2005) Buballa, M. 2005, Phys. Rep., 407, 205. doi:10.1016/j.physrep.2004.11.004
  * Buchner (2016) Buchner, J. 2016, PyMultiNest: Python interface for MultiNest, Astrophysics Source Code Library, ascl:1606.005
  * Cao et al. (2020) Cao, Z., Chen, L.-W., Chu, P.-C., et al. 2020, arXiv:2009.00942
  * Char et al. (2019) Char, P., Drago, A., & Pagliara, G. 2019, in AIP Conf. Proc. 2127, Xiamen-CUSTIPEN Workshop on the Equation of State of Dense Neutron-rich Matter in the Era of Gravitational Wave Astronomy, ed. A. Li et al. (Melville, NY: AIP). doi:10.1063/1.5117816
  * Chatziioannou (2020) Chatziioannou, K. 2020, General Relativity and Gravitation, 52, 109. doi:10.1007/s10714-020-02754-3
  * Cheng & Dai (1996) Cheng, K. S. & Dai, Z. G. 1996, Phys. Rev. Lett., 77, 1210. doi:10.1103/PhysRevLett.77.1210
  * Chu & Chen (2017) Chu, P.-C. & Chen, L.-W. 2017, Phys. Rev. D, 96, 083019. doi:10.1103/PhysRevD.96.083019
  * Cromartie et al. (2020) Cromartie, H. T., Fonseca, E., Ransom, S. M., et al. 2020,NatAs, 4, 72. doi:10.1038/s41550-019-0880-2
  * De Pietri et al. (2019) De Pietri, R., Drago, A., Feo, A., et al. 2019, ApJ, 881, 122. doi:10.3847/1538-4357/ab2fd0
  * Degrand et al. (1975) Degrand, T., Jaffe, R. L., Johnson, K., et al. 1975, Phys. Rev. D, 12, 2060. doi:10.1103/PhysRevD.12.2060
  * Dey et al. (1998) Dey, M., Bombaci, I., Dey, J., et al. 1998, PhLB, 438, 123. doi:10.1016/S0370-2693(98)00935-6
  * Dietrich et al. (2017) Dietrich, T., Bernuzzi, S., & Tichy, W. 2017, Phys. Rev. D, 96, 121501. doi:10.1103/PhysRevD.96.121501
  * Dietrich et al. (2021) Dietrich, T., Hinderer, T., & Samajdar, A. 2021, General Relativity and Gravitation, 53, 27. doi:10.1007/s10714-020-02751-6
  * Drago et al. (2014) Drago, A., Lavagno, A., & Pagliara, G. 2014, Phys. Rev. D, 89, 043014. doi:10.1103/PhysRevD.89.043014
  * Drago & Pagliara (2020) Drago, A. & Pagliara, G. 2020, Phys. Rev. D, 102, 063003. doi:10.1103/PhysRevD.102.063003
  * Drago et al. (2007) Drago, A., Lavagno, A., & Parenti, I. 2007, ApJ, 659, 1519. doi:10.1086/512112
  * Farhi & Jaffe (1984) Farhi, E. & Jaffe, R. L. 1984, Phys. Rev. D, 30, 2379. doi:10.1103/PhysRevD.30.2379
  * Fasano et al. (2019) Fasano, M., Abdelsalhin, T., Maselli, A., et al. 2019, Phys. Rev. Lett., 123, 141101. doi:10.1103/PhysRevLett.123.141101
  * Fattoyev et al. (2018) Fattoyev, F. J., Piekarewicz, J., & Horowitz, C. J. 2018, Phys. Rev. Lett., 120, 172702. doi:10.1103/PhysRevLett.120.172702
  * Fonseca et al. (2021) Fonseca, E., Cromartie, H. T., Pennucci, T. T., et al. 2021, ApJ, 915, L12. doi:10.3847/2041-8213/ac03b8
  * Fraga et al. (2001) Fraga, E. S., Pisarski, R. D., & Schaffner-Bielich, J. 2001, Phys. Rev. D, 63, 121702. doi:10.1103/PhysRevD.63.121702
  * Glendenning (1990) Glendenning, N. K. 1990, MPLA, 5, 2197. doi:10.1142/S021773239000250X
  * Gorda et al. (2018) Gorda, T., Kurkela, A., Romatschke, P., et al. 2018, Phys. Rev. Lett., 121, 202701. doi:10.1103/PhysRevLett.121.202701
  * Haensel et al. (1986) Haensel, P., Zdunik, J. L., & Schaefer, R. 1986, A&A, 160, 121
  * Haxton & Heller (1980) Haxton, W. C. & Heller, L. 1980, Phys. Rev. D, 22, 1198. doi:10.1103/PhysRevD.22.1198
  * Hernandez Vivanco et al. (2020) Hernandez Vivanco, F., Smith, R., Thrane, E., et al. 2020, MNRAS, 499, 5972. doi:10.1093/mnras/staa3243
  * Holdom et al. (2018) Holdom, B., Ren, J., & Zhang, C. 2018, Phys. Rev. Lett., 120, 222001. doi:10.1103/PhysRevLett.120.222001
  * Jiang et al. (2020) Jiang, J.-L., Tang, S.-P., Wang, Y.-Z., et al. 2020, ApJ, 892, 55. doi:10.3847/1538-4357/ab77cf
  * Joshi et al. (2021) Joshi, S., Sau, S., & Sanyal, S. 2021, Journal of High Energy Astrophysics, 30, 16. doi:10.1016/j.jheap.2021.03.001
  * Klähn & Fischer (2015) Klähn, T. & Fischer, T. 2015, ApJ, 810, 134. doi:10.1088/0004-637X/810/2/134
  * Kurkela et al. (2014) Kurkela, A., Fraga, E. S., Schaffner-Bielich, J., et al. 2014, ApJ, 789, 127. doi:10.1088/0004-637X/789/2/127
  * Kurkela et al. (2010) Kurkela, A., Romatschke, P., & Vuorinen, A. 2010, Phys. Rev. D, 81, 105021. doi:10.1103/PhysRevD.81.105021
  * Landry et al. (2020) Landry, P., Essick, R., & Chatziioannou, K. 2020, Phys. Rev. D, 101, 123007. doi:10.1103/PhysRevD.101.123007
  * Levan et al. (2017) Levan, A. J., Lyman, J. D., Tanvir, N. R., et al. 2017, ApJ, 848, L28. doi:10.3847/2041-8213/aa905f
  * Li et al. (2020) Li, A., Miao, Z.-Q., Jiang, J.-L., et al. 2020, arXiv:2009.12571. doi:10.1093/mnras/stab2029
  * Li et al. (2021b) Li, A., Miao, Z.-Q., Han, S., et al. 2021, ApJ, 913, 27. doi:10.3847/1538-4357/abf355
  * Li et al. (2011) Li, A., Peng, G.-X., & Lu, J.-F. 2011, Research in Astronomy and Astrophysics, 11, 482. doi:10.1088/1674-4527/11/4/010
  * Li et al. (2010) Li, A., Xu, R.-X., & Lu, J.-F. 2010, MNRAS, 402, 2715. doi:10.1111/j.1365-2966.2009.16091.x
  * Li et al. (2016) Li, A., Zhang, B., Zhang, N.-B., et al. 2016, Phys. Rev. D, 94, 083010. doi:10.1103/PhysRevD.94.083010
  * Li et al. (2017) Li, A., Zhu, Z.-Y., & Zhou, X. 2017, ApJ, 844, 41. doi:10.3847/1538-4357/aa7a00
  * Li et al. (2020) Li, A., Zhu, Z.-Y., Zhou, E.-P., et al. 2020, Journal of High Energy Astrophysics, 28, 19. doi:10.1016/j.jheap.2020.07.001
  * Lopes et al. (2021) Lopes, L. L., Biesdorf, C., Marquez, K. D., et al. 2021, Phys. Scr, 96, 065302. doi:10.1088/1402-4896/abef35
  * Lugones & Horvath (2003) Lugones, G. & Horvath, J. E. 2003, A&A, 403, 173. doi:10.1051/0004-6361:20030374
  * Lugones & Horvath (2002) Lugones, G. & Horvath, J. E. 2002, Phys. Rev. D, 66, 074017. doi:10.1103/PhysRevD.66.074017
  * McLerran & Reddy (2019) McLerran, L. & Reddy, S. 2019, Phys. Rev. Lett., 122, 122701. doi:10.1103/PhysRevLett.122.122701
  * Miao et al. (2020) Miao, Z., Li, A., Zhu, Z., et al. 2020, ApJ, 904, 103. doi:10.3847/1538-4357/abbd41
  * Miller et al. (2021) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2021, arXiv:2105.06979
  * Miller et al. (2019) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, ApJ, 887, L24. doi:10.3847/2041-8213/ab50c5
  * Most et al. (2018) Most, E. R., Weih, L. R., Rezzolla, L., et al. 2018, Phys. Rev. Lett., 120, 261103. doi:10.1103/PhysRevLett.120.261103
  * Rajagopal & Wilczek (2001) Rajagopal, K. & Wilczek, F. 2001, Phys. Rev. Lett., 86, 3492. doi:10.1103/PhysRevLett.86.3492
  * Rapp et al. (1998) Rapp, R., Schäfer, T., Shuryak, E., et al. 1998, Phys. Rev. Lett., 81, 53. doi:10.1103/PhysRevLett.81.53
  * Riley et al. (2019) Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, ApJ, 887, L21. doi:10.3847/2041-8213/ab481c
  * Riley et al. (2021) Riley, T. E., Watts, A. L., Ray, P. S., et al. 2021, arXiv:2105.06980
  * Roupas et al. (2021) Roupas, Z., Panotopoulos, G., & Lopes, I. 2021, Phys. Rev. D, 103, 083015. doi:10.1103/PhysRevD.103.083015
  * Sedaghat et al. (2021) Sedaghat, J., Zebarjad, S. M., Bordbar, G. H., et al. 2021, arXiv:2104.00544
  * Tang et al. (2020) Tang, S.-P., Jiang, J.-L., Gao, W.-H., et al. 2020, ApJ, 888, 45. doi:10.3847/1538-4357/ab5959
  * Tews et al. (2018) Tews, I., Margueron, J., & Reddy, S. 2018, Phys. Rev. C, 98, 045804. doi:10.1103/PhysRevC.98.045804
  * Traversi et al. (2021) Traversi, S., Char, P., Pagliara, G., et al. 2021, arXiv:2102.02357
  * Traversi & Char (2020) Traversi, S. & Char, P. 2020, ApJ, 905, 9. doi:10.3847/1538-4357/abbfb4
  * Weber (2005) Weber, F. 2005, Progress in Particle and Nuclear Physics, 54, 193. doi:10.1016/j.ppnp.2004.07.001
  * Witten (1984) Witten, E. 1984, Phys. Rev. D, 30, 272. doi:10.1103/PhysRevD.30.272
  * Xia et al. (2021) Xia, C., Zhu, Z., Zhou, X., et al. 2021, ChPhyC, 45, 055104. doi:10.1088/1674-1137/abea0d
  * Zhang (2020) Zhang, C. 2020, Phys. Rev. D, 101, 043003. doi:10.1103/PhysRevD.101.043003
  * Zhang & Mann (2021) Zhang, C. & Mann, R. B. 2021, Phys. Rev. D, 103, 063018. doi:10.1103/PhysRevD.103.063018
  * Zhou et al. (2018) Zhou, E.-P., Zhou, X., & Li, A. 2018, Phys. Rev. D, 97, 083015. doi:10.1103/PhysRevD.97.083015
  * Zhu et al. (2018) Zhu, Z.-Y., Zhou, E.-P., & Li, A. 2018, ApJ, 862, 98. doi:10.3847/1538-4357/aacc28