#  Supplemental Material for “Lindbladian Many-Body Localization”

Ryusuke Hamazaki  Nonequilibrium Quantum Statistical Mechanics RIKEN Hakubi
Research Team, RIKEN Cluster for Pioneering Research (CPR), RIKEN iTHEMS,
Wako, Saitama 351-0198, Japan  Masaya Nakagawa  Department of Physics,
University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan  Taiki Haga
Department of Physics and Electronics, Osaka Metropolitan University, Osaka,
599-8531, Japan  Masahito Ueda  Department of Physics, University of Tokyo,
7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan  RIKEN Center for Emergent
Matter Science (CEMS), Wako, Saitama 351-0198, Japan  Institute for Physics of
Intelligence, University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033,
Japan

###### Contents

  1. I Relaxation dynamics
    1. I.1 Transverse and longitudinal relaxations
    2. I.2 Rigid decay rate for the dynamics of coherence
      1. I.2.1 Rigid decay rate deep in the Lindbladian MBL
      2. I.2.2 Comparison with the decay rate for a single-spin system
      3. I.2.3 Stability against additional interactions
  2. II Structures of $|\mathrm{Tr}[\hat{O}\hat{R}_{\alpha}]|,|\mathrm{Tr}[\hat{L}_{\alpha}^{\dagger}\hat{\rho}]|$, and $|\mathrm{Tr}[\hat{R}_{\alpha}^{\dagger}\hat{L}_{\alpha}]|^{-1}$
    1. II.1 System-size dependences of $\mathcal{C},\mathcal{B}_{\hat{\rho}},$ and $\mathcal{A}_{\hat{O}}$
    2. II.2 Non-Hermitian random-matrix prediction of eigenstates
    3. II.3 Effect of locality on $\mathrm{Tr}[\hat{O}\hat{R}_{\alpha}]$
    4. II.4 Similar hypothesis for classical dynamics
  3. III Phenomenology of the Lindbladian MBL phase
    1. III.1 Structure of the Lindblad generator
    2. III.2 Case of dephasing
    3. III.3 Case of damping
  4. IV Effects of dissipation on the Lindbladian MBL transition point
    1. IV.1 Phase diagram
    2. IV.2 Relation to the many-body localization without dissipation

##  I Relaxation dynamics

###  I.1 Transverse and longitudinal relaxations

Figure S-1:  Transverse and longitudinal dynamics from an initial state given
in Eq. (S-1) for the disordered Ising model with (left) dephasing and (right)
damping types of dissipation. (a) Transverse relaxation of spins. We consider
the time evolutions of $\braket{\hat{\sigma}_{L}^{x}(t)}$ (solid) and
$\braket{\hat{\sigma}_{L-1}^{x}(t)}$ (dashed) for $h=1.2$ (blue) and $h=10$
(red). For $\braket{\hat{\sigma}_{L}^{x}(t)}$, we find rapid many-body
decoherence for weak disorder and slower oscillatory decay (with time scale
$\propto\gamma^{-1}$) for strong disorder. On the other hand,
$\braket{\hat{\sigma}_{L-1}^{x}(t)}$ does not show significant oscillations
for all cases since the initial coherence at site $L-1$ is absent. (b)
Longitudinal relaxation of spins. We consider time evolutions of
$\braket{\hat{\sigma}_{L}^{z}(t)}$ (solid) and
$\braket{\hat{\sigma}_{L-1}^{z}(t)}$ (dashed) for $h=1.2$ (blue) and $h=10$
(red). For dephasing, $\braket{\hat{\sigma}_{L-1}^{z}(t)}$ exhibits a slow
decay. Similarly, $\braket{\hat{\sigma}_{L}^{z}(t)}$ decays slowly, although
small oscillation appears. For the damping case,
$\braket{\hat{\sigma}_{L-1}^{z}(t)}$ is almost constant around $-1$, and
$\braket{\hat{\sigma}_{L}^{z}(t)}$ exhibits a decay to $-1$ with a time scale
$\sim\gamma^{-1}$. We use a single disorder realization and take
$L=6,J=1,g=-0.9$, and $\gamma=0.1$.

Here, we numerically demonstrate the difference between longitudinal and
transverse (coherence) relaxation dynamics. We consider the time evolution
from the initial state with $q=1$ in the main text, i.e.,

$\displaystyle\ket{\psi_{q=1}}=(\ket{\downarrow}_{1}\cdots\ket{\downarrow}_{L-1})\left(\frac{\ket{\uparrow}_{L}+\ket{\downarrow}_{L}}{\sqrt{2}}\right).$
(S-1)

Figure S-1(a) shows the dynamics of $\braket{\hat{\sigma}_{L}^{x}(t)}$ and
$\braket{\hat{\sigma}_{L-1}^{x}(t)}$, which describe a transverse relaxation
of spins. For $\braket{\hat{\sigma}_{L}^{x}(t)}$, we find rapid many-body
decoherence for weak disorder and slower oscillatory decay (with time scale
$\propto\gamma^{-1}$) for strong disorder. As discussed in the main text, the
latter is attributed to the localization of off-diagonal degrees of freedom in
the Lindbladian spectrum. Note that $\braket{\hat{\sigma}_{L-1}^{x}(t)}$ does
not show significant oscillations for all cases since the coherence at site
$L-1$ is absent from the beginning of the dynamics.

Figure S-1(b) shows the dynamics of $\braket{\hat{\sigma}_{L}^{z}(t)}$ and
$\braket{\hat{\sigma}_{L-1}^{z}(t)}$, which describe the longitudinal
relaxation of spins. For small $h$, we find a rapid relaxation to a stationary
value for both dephasing and damping cases. On the other hand, for strong
disorder, the dynamics differ for these two different types of dissipation.
For dephasing, $\braket{\hat{\sigma}_{L-1}^{z}(t)}$ exhibits slow decay, which
is consistent with previous studies [1, 2, 3]. This is attributed to the
delocalization of diagonal degrees of freedom in the Lindbladian spectrum
(especially around the zero eigenvalue). Similarly,
$\braket{\hat{\sigma}_{L}^{z}(t)}$ decays slowly, although small oscillations
appear due to the small overlap between $\hat{\sigma}_{L}^{z}$ and
$\hat{\tau}_{L}^{x}$ (see Appendix IIIA for the definition of
$\hat{\tau}_{i}^{x}$). For the damping case, the stationary state is
approximately described by a state in which spins are all down. Then,
$\braket{\hat{\sigma}_{L-1}^{z}(t)}$ is almost constant around $-1$, and
$\braket{\hat{\sigma}_{L}^{z}(t)}$ exhibits a decay to $-1$ with a time scale
$\sim(2\gamma)^{-1}$. Note that this time scale for the longitudinal
relaxation $\sim(2\gamma)^{-1}$ in the strong-disorder regime is twice as long
as the time scale for the transverse relaxation $\sim(\gamma)^{-1}$ under
damping-type dissipation. While the factor two is known in the two-level
dynamics described by a quantum-optical Master equation with damping [4], it
is surprising that the same factor appears even for many-body systems showing
the Lindbladian MBL.

###  I.2 Rigid decay rate for the dynamics of coherence

####  I.2.1 Rigid decay rate deep in the Lindbladian MBL

Here, we numerically demonstrate that the decay rate for the dynamics of
coherence approaches $\sim 2\gamma$ (dephasing) or $\sim\gamma$ (damping) for
sufficiently strong disorder. For this purpose, let us consider the dynamics
of $\braket{\hat{\sigma}_{L}^{x}(t)}$ discussed in the previous subsection and
fit it with the following function:

$\displaystyle f(t)=e^{-at}(\cos bt+(1-c)\sin bt)$ (S-2)

with fitting parameters $a,b$, and $c$. The fitting function is motivated as
follows: (i) $\braket{\hat{\sigma}_{L}^{x}(0)}=1$, and (ii) an exponential
decay of $\hat{\tau}_{L}^{x}$ is accompanied by an oscillation at an almost
fixed frequency (corresponding to $2\tilde{h}_{l}$ in Eqs. (S-36) or (S-41))
in the deep Lindbladian MBL phase. Note that the fitting may not work well for
the delocalized phase. Even in the Lindbladian MBL phase, the fitting may not
work (e.g., an oscillation with multiple frequencies appears) due to the
difference between $\hat{\tau}_{l}$ and $\hat{\sigma}_{l}$. Nevertheless,
given that the fitting works well in the deep Lindbladian MBL phase with
sufficiently strong disorder, its decay rate $a$ will be close to $\sim
2\gamma$ (dephasing) or $\sim\gamma$ (damping).

Table S-1 shows the number of samples whose coefficient of determination $R$
for the fitting satisfies $R\geq 0.9$ among 100 samples. For sufficiently
large $h$, the fitting with Eq. (S-2) becomes better for many samples, as
expected. Note that the fit seems to work well for small $h$ as well, but this
is not due to localization, as the fitting parameter $a$ is far from $\gamma$.
Indeed, as shown in Fig. S-2, the fitting parameter $a$ averaged over the
samples with $R\geq 0.9$ decreases for $h\geq 2.0$ and approaches $\sim
2\gamma$ (dephasing) or $\sim\gamma$ (damping) as $h$ is increased.

|  $h=0.6$ |  1.2 |  2.0 |  3.0 |  4.0 |  5.0 |  7.5 |  10.0 |  20.0
---|---|---|---|---|---|---|---|---|---
dephasing ($\gamma=0.05$) | 4 | 16 | 11 | 12 | 29 | 43 | 55 | 74 | 80
dephasing ($\gamma=0.1$) | 8 | 23 | 16 | 16 | 28 | 42 | 58 | 69 | 79
damping ($\gamma=0.05$) | 18 | 31 | 12 | 8 | 13 | 27 | 46 | 65 | 77
damping ($\gamma=0.1$) | 59 | 45 | 15 | 9 | 16 | 26 | 45 | 62 | 76
Table S-1: Number of samples whose coefficient of determination $R$ for the
fitting satisfies $R\geq 0.9$ among 100 samples. For sufficiently large $h$,
the fitting in Eq. (S-2) becomes better for many samples. The fitting is made
with the use of the SciPy curve_fit function written in Python, and samples
that fail to converge within the number of calls maxfev=800 are neglected.
Figure S-2:  Fitting parameter $a$ averaged over the samples with $R\geq 0.9$,
where we use $\gamma=0.05$ (dashed) and $\gamma=1.0$ (solid). It decreases for
$h\geq 2.0$ and approaches $\sim 2\gamma$ for dephasing (blue) or $\sim\gamma$
for damping (red) as $h$ is increased.

####  I.2.2 Comparison with the decay rate for a single-spin system

Here, we compare the result in the previous subsection with the decay rate for
a single-spin system. For a single-spin system, the Lindblad equation reads

$\displaystyle\frac{\mathrm{d}\hat{\rho}}{\mathrm{d}t}=-i[g\hat{\sigma}_{1}^{x}+h_{1}\hat{\sigma}_{1}^{z},\hat{\rho}]+\frac{1}{2}\left[2\hat{\Gamma}_{1}\hat{\rho}\hat{\Gamma}_{1}^{\dagger}-{\left\\{\hat{\Gamma}_{1}^{\dagger}\hat{\Gamma}_{1},\hat{\rho}\right\\}}\right]$
(S-3)

with $\Gamma_{1}=\sqrt{\gamma}\hat{\sigma}_{1}^{z}$ for dephasing and
$\Gamma_{1}=\sqrt{\gamma/2}\hat{\sigma}_{1}^{-}$ for damping.

For simplicity, let us assume that $\gamma$ is small enough. In that case, the
four eigenvalues of the Liouvillian super-operator are given by

$\displaystyle
0,\quad-\frac{2g^{2}}{g^{2}+h_{1}^{2}}\gamma+\mathrm{O}(\gamma^{2}),\quad-\frac{g^{2}+2h_{1}^{2}}{g^{2}+h_{1}^{2}}\gamma\pm
2\sqrt{g^{2}+h_{1}^{2}}i+\mathrm{O}(\gamma^{2})$ (S-4)

for dephasing, and

$\displaystyle
0,\quad-\frac{g^{2}+2h_{1}^{2}}{g^{2}+h_{1}^{2}}\gamma+\mathrm{O}(\gamma^{2}),\quad-\frac{3g^{2}+2h_{1}^{2}}{2(g^{2}+h_{1}^{2})}\gamma\pm
2\sqrt{g^{2}+h_{1}^{2}}i+\mathrm{O}(\gamma^{2})$ (S-5)

for damping.

For the strong disorder $h_{1}\gg g$, the eigenvalues are approximated by
$\sim 0,0,2\gamma\pm 2h_{1}i$ for dephasing and $\sim 0,-2\gamma,-\gamma\pm
2h_{1}i$ for damping. Thus, the longest lived modes have decay rates $2\gamma$
and $\gamma$ for dephasing and damping, respectively, as argued in the
previous subsection and the main text. We also note that the real parts for
all the eigenvalues are larger than $-2\gamma$ for any $h_{1}$ and $g$. This
fact suggests that the large value of $a$ for small $h$ in the fitting to
$f(t)$ (see Eq. (S-2) and Fig. S-2) is not attributed to the decay rate for a
single-spin system.

####  I.2.3 Stability against additional interactions

To strengthen our argument, we demonstrate that the decay rate for large $h$
is robust against additional interactions distinct from the Ising one. Here,
we consider the case where the Hamiltonian is perturbed by an XY-type
interaction, i.e.,

$\displaystyle\hat{H}^{\prime}=\hat{H}+J_{XY}\sum_{l=1}^{L-1}(\hat{\sigma}_{l}^{x}\hat{\sigma}_{l+1}^{x}+\hat{\sigma}_{l}^{y}\hat{\sigma}_{l+1}^{y}).$
(S-6)

We can perform the same fitting procedure based on Eq. (S-2) for the dynamics
of coherence under the Hamiltonian $\hat{H}^{\prime}$. First, Table S-2 shows
the number of samples whose coefficient of determination $R$ for the fitting
satisfies $R\geq 0.9$ among 100 samples. Second, in Fig. S-3, we show the
fitting parameter $a$ averaged over the samples with $R\geq 0.9$. These two
results are similar to the case without the XY interaction. In particular, $a$
decreases for $h\geq 2.0$ and approaches $\sim 2\gamma$ (dephasing) or
$\sim\gamma$ (damping) as $h$ is increased, which indicates the stability of
the decay rate against perturbation.

|  $h=0.6$ |  1.2 |  2.0 |  3.0 |  4.0 |  5.0 |  7.5 |  10.0 |  20.0
---|---|---|---|---|---|---|---|---|---
dephasing ($\gamma=0.05$) | 9 | 11 | 4 | 12 | 27 | 38 | 53 | 68 | 78
dephasing ($\gamma=0.1$) | 14 | 22 | 14 | 20 | 26 | 39 | 54 | 67 | 77
damping ($\gamma=0.05$) | 18 | 34 | 8 | 5 | 13 | 23 | 42 | 57 | 75
damping ($\gamma=0.1$) | 58 | 57 | 21 | 14 | 15 | 23 | 41 | 57 | 74
Table S-2: Number of samples whose coefficient of determination $R$ for the
fitting satisfies $R\geq 0.9$ among 100 samples for the case where the
Hamiltonian is perturbed by the XY interaction ($J_{XY}=0.2$). The fitting is
made with the use of the SciPy curve_fit function written in Python, and
samples that fail to converge within the number of calls maxfev=800 are
neglected. Figure S-3:  Fitting parameter $a$ averaged over the samples with
$R\geq 0.9$ for the case where the Hamiltonian is perturbed by the XY
interaction ($J_{XY}=0.2$). Here, we use $\gamma=0.05$ (dashed) and
$\gamma=1.0$ (solid). The parameter decreases for $h\geq 2.0$ and approaches
$\sim 2\gamma$ for dephasing (blue) or $\sim\gamma$ for damping (red) as $h$
is increased.

##  II Structures of
$|\mathrm{Tr}[\hat{O}\hat{R}_{\alpha}]|,|\mathrm{Tr}[\hat{L}_{\alpha}^{\dagger}\hat{\rho}]|$,
and $|\mathrm{Tr}[\hat{R}_{\alpha}^{\dagger}\hat{L}_{\alpha}]|^{-1}$

###  II.1 System-size dependences of $\mathcal{C},\mathcal{B}_{\hat{\rho}},$
and $\mathcal{A}_{\hat{O}}$

Figure S-4: Structures of
$|\mathrm{Tr}[\hat{O}\hat{R}_{\alpha}]|,|\mathrm{Tr}[\hat{L}_{\alpha}^{\dagger}\hat{\rho}]|$,
and $|\mathrm{Tr}[\hat{L}_{\alpha}^{\dagger}\hat{R}_{\alpha}]|^{-1}$ for one
sample, as a function of real and imaginary parts of $\lambda_{\alpha}$. While
$|\mathrm{Tr}[\hat{O}\hat{R}_{\alpha}]|$ and
$|\mathrm{Tr}[\hat{L}_{\alpha}^{\dagger}\hat{\rho}]|$ tend to decrease with
$L$, $|\mathrm{Tr}[\hat{L}_{\alpha}^{\dagger}\hat{R}_{\alpha}]|^{-1}$
increases with increasing $L$. We consider $L=7$ (blue) and $L=5$ (red) for
each plot. We use the dephasing-type dissipation and set $q=1$,
$J=1,g=-0.9,h=1.2$, and $\gamma=0.5$.  Figure S-5: Average values of
$|\mathrm{Tr}[\hat{L}_{\alpha}^{\dagger}\hat{R}_{\alpha}]|^{-1},|\mathrm{Tr}[\hat{L}_{\alpha}^{\dagger}\hat{\rho}]|,$
and $|\mathrm{Tr}[\hat{O}\hat{R}_{\alpha}]|$ for different values of $h$
($=0.6,1.2,1.8,2.4,4.0,5.0,10.0$ from red to orange) and for dephaing (top)
and damping (bottom) types of dissipation. We see that
$|\mathrm{Tr}[\hat{L}_{\alpha}^{\dagger}\hat{R}_{\alpha}]|^{-1}\propto
D=2^{L}$ and $\mathrm{Tr}[\hat{L}_{\alpha}^{\dagger}\hat{\rho}]\propto D^{-1}$
for small $h$. We use the Lindbladian with $L=7,J=1,g=-0.9$, and $\gamma=0.5$.
For the calculation of the average, we consider eigenvalues within the range
$\mathrm{Re}[\lambda_{\alpha}]\in[0.6\lambda^{r}_{\mathrm{min}},0.4\lambda^{r}_{\mathrm{min}}]$
and
$\mathrm{Im}[\lambda_{\alpha}]\in[0.6\lambda^{i}_{\mathrm{min}},0.4\lambda^{i}_{\mathrm{min}}]$,
where $\lambda^{r}_{\mathrm{min}}=\min_{\alpha}\mathrm{Re}[\lambda_{\alpha}]$
and $\lambda^{i}_{\mathrm{min}}=\min_{\alpha}\mathrm{Im}[\lambda_{\alpha}]$.
The number of samples used for the simulation is $10000\>(L=4,5),800\>(L=6),$
and $50\>(L=7)$.

We here discuss the system-size dependences of
$\mathcal{C},\mathcal{B}_{\hat{\rho}},$ and $\mathcal{A}_{\hat{O}}$ (see Eqs.
(3)-(5) in the main text) in the middle of the spectrum. Figure S-4 shows
$|\mathrm{Tr}[\hat{O}\hat{R}_{\alpha}]|,|\mathrm{Tr}[\hat{L}_{\alpha}^{\dagger}\hat{\rho}]|$,
and $|\mathrm{Tr}[\hat{R}_{\alpha}^{\dagger}\hat{L}_{\alpha}]|^{-1}$ for one
sample for system sizes $L=5$ and $L=7$. While
$|\mathrm{Tr}[\hat{O}\hat{R}_{\alpha}]|$ and
$|\mathrm{Tr}[\hat{L}_{\alpha}^{\dagger}\hat{\rho}]|$ tend to decrease with
$L$, $|\mathrm{Tr}[\hat{L}_{\alpha}^{\dagger}\hat{R}_{\alpha}]|^{-1}$
increases for increasing $L$. Figure S-5 shows finite-size scaling of these
factors averaged over many samples, where the disorder strength $h$ is varied.
We find that $|\mathrm{Tr}[\hat{L}_{\alpha}^{\dagger}\hat{R}_{\alpha}]|^{-1}$
is almost proportional to $D=2^{L}$ for small $h$, where $D$ is the dimension
of the Hilbert space. In addition, we also find
$|\mathrm{Tr}[\hat{L}_{\alpha}^{\dagger}\hat{\rho}]|\propto D^{-1}$ for small
$h$. These scalings agree with the prediction of the non-Hermitian RMT (see
Appendix IIB). On the other hand, for strong $h$, the scaling behavior differs
from the RMT for both cases. For $|\mathrm{Tr}[\hat{O}\hat{R}_{\alpha}]|$, the
scaling behavior is not straightforward because of the locality of the
operator and that of the Lindblad generator (see Appendix IIC).

###  II.2 Non-Hermitian random-matrix prediction of eigenstates

The expressions in Eqs. (3), (4), and (5) in the main text are motivated by
the non-Hermitian random-matrix theory: Let
$\ket{R_{\alpha}^{G}}/\ket{L_{\alpha}^{G}}$ be right/left eigenstates of an
$N\times N$ non-Hermitian real Ginibre random matrix. Then we obtain

$\displaystyle\braket{\mathcal{O}}{R_{\alpha}^{G}},\>\braket{L_{\alpha}^{G}}{\mathcal{O}}\sim\frac{\sqrt{\braket{\mathcal{O}}{\mathcal{O}}}}{\sqrt{N}}r^{G}_{\alpha}$
(S-7)

and

$\displaystyle\braket{L_{\alpha}^{G}}{R_{\alpha}^{G}}^{-1}\sim\sqrt{N}r^{G^{\prime}}_{\alpha},$
(S-8)

where $\ket{\mathcal{O}}=(\bra{\mathcal{O}})^{\dagger}$ is an arbitrary vector
with dimension $N$, and $r^{G}_{\alpha}$ and $r^{G^{\prime}}_{\alpha}$ are
random variables with zero mean and unit variance. The non-orthogonality of
eigenstates in Eq. $\eqref{northo}$ was studied previously in Refs. [5, 6, 7,
8, 9].

Let us assume that $\lambda_{\alpha}$ has a nonzero imaginary part. In this
case, we find that the phases of $r_{\alpha}^{G}$ and
$r^{G^{\prime}}_{\alpha}$ are randomly distributed over $[0,2\pi)$. The
squared absolute value of $r_{\alpha}^{G}$, $|r_{\alpha}^{G}|^{2}$, is
essentially given by the Porter-Thomas distribution. From this we can obtain

$\displaystyle
P_{G}(x):=P(|r_{\alpha}^{G}|=x)=\frac{\pi}{2}xe^{-\frac{\pi}{2}x^{2}}$ (S-9)

after a normalization to set the average of $|r_{\alpha}^{G}|$ unity.
Moreover, $P(|r^{G^{\prime}}_{\alpha}|=x)$ becomes a universal distribution
for non-Hermitian random matrices, and it is expected to take the form of

$\displaystyle
P_{G}^{\prime}(x):=P(|r_{\alpha}^{G^{\prime}}|=x)=\frac{32}{\pi^{2}x^{5}}e^{-\frac{4}{\pi
x^{2}}}.$ (S-10)

Note that this function is originally obtained from the analytical results of
a complex Ginibre ensemble [6]. While the analytical formula for a real
Ginibre ensemble is not known, we have numerically confirmed this
distribution. If $\lambda_{\alpha}$ is real, then $r_{\alpha}^{G}$ and
$r^{G^{\prime}}_{\alpha}$ can be taken to be real, and their distributions may
be different from the case of complex eigenvalues.

Now, we note that any matrix can have a vector representation through the
vectorization mapping (see, e.g., Ref. [10])

$\displaystyle\hat{A}=\sum_{ij}(\hat{A})_{ij}\ket{i}\bra{j}\rightarrow|A)=\sum_{ij}(\hat{A})_{ij}\ket{i}\otimes\ket{j}.$
(S-11)

Using this representation, we have a form similar to the left-hand sides of
Eqs. (S-7) and (S-8):

$\displaystyle\mathrm{Tr}[\hat{O}\hat{R}_{\alpha}]=({O}|{R}_{\alpha}),$ (S-12)
$\displaystyle\mathrm{Tr}[\hat{L}_{\alpha}^{\dagger}\hat{\rho}]=(L_{\alpha}|\rho),$
(S-13)

and

$\displaystyle\mathrm{Tr}[\hat{L}_{\alpha}^{\dagger}\hat{R}_{\alpha}]=({L}_{\alpha}|{R}_{\alpha}),$
(S-14)

where $(A|=\\{|A)\\}^{\dagger}$. Then, in analogy with the eigenstate
thermalization hypothesis, which states that eigenstates of generic
Hamiltonians obey the statistics of a random matrix, we conjecture that Eqs.
(3), (4) and (5) in the main text hold for generic super-operators with
$r_{\alpha}$ and $r_{\alpha}^{\prime}$ obeying the distributions $P_{G}$ and
$P_{G}^{\prime}$, respectively. Furthermore, owing to the $N$-scaling in Eq.
(S-8), we can also conjecture that
$\mathcal{C}(\lambda_{\alpha})\propto\sqrt{D^{2}}=D$ (note that $|R_{\alpha})$
and the other vectorized states have dimension $D^{2}=2^{2L}$). Numerical
simulations indicate that the similar scaling also occurs for
$\mathrm{Tr}[\hat{L}_{\alpha}^{\dagger}\hat{\rho}]\sim D^{-1}$ (especially for
pure states, where $(\rho|\rho)=1$) but not for
$\mathrm{Tr}[\hat{O}\hat{R}_{\alpha}]$ owing to the locality of our model (see
Appendix IIC).

###  II.3 Effect of locality on $\mathrm{Tr}[\hat{O}\hat{R}_{\alpha}]$

We here argue that $|\mathrm{Tr}[\hat{O}\hat{R}_{\alpha}]|$ is exponentially
suppressed as a function of $|\lambda_{\alpha}|$ for a local Lindbladian and a
local observable $\hat{O}$, unless many of $|R_{\alpha})$ are almost parallel.
To explain the idea behind this argument, we consider

$\displaystyle\sum_{|\lambda_{\alpha}|>\lambda}|\mathrm{Tr}[\hat{O}\hat{R}_{\alpha}]|^{2}$
(S-15)

for some positive $\lambda$. Using
$\mathcal{L}[\hat{R}_{\alpha}]=\lambda_{\alpha}\hat{R}_{\alpha}$, we find

$\displaystyle|\mathrm{Tr}[\hat{O}\hat{R}_{\alpha}]|^{2}=\frac{|\mathrm{Tr}[\hat{O}\mathcal{L}^{k}[\hat{R}_{\alpha}]]|^{2}}{|\lambda_{\alpha}|^{2k}}=\frac{|\mathrm{Tr}[\mathcal{L}^{{\dagger}k}[\hat{O}]\hat{R}_{\alpha}]|^{2}}{|\lambda_{\alpha}|^{2k}}$
(S-16)

for an arbitrary positive integer $k$. Here, by taking
$\mathcal{L}^{{\dagger}k}[\hat{O}]=:\hat{A}=\hat{A}^{\dagger}$, we have

$\displaystyle\sum_{\alpha}|\mathrm{Tr}[\mathcal{L}^{{\dagger}k}[\hat{O}]\hat{R}_{\alpha}]|^{2}$
$\displaystyle=\sum_{\alpha}\sum_{\mu\nu\xi\eta}A_{\mu\nu}^{*}R^{\alpha}_{\mu\nu}R_{\xi\eta}^{\alpha*}A_{\xi\eta}$
$\displaystyle=\sum_{\mu\nu\xi\eta}A_{\mu\nu}^{*}X_{\mu\nu;\xi\eta}A_{\xi\eta}$
$\displaystyle\leq\sum_{\mu\nu}|A_{\mu\nu}|^{2}||X||_{\mathrm{OP}}=\mathrm{Tr}[\hat{A}^{2}]||X||_{\mathrm{OP}}\leq
D||\hat{A}||_{\mathrm{\infty}}^{2}||X||_{\mathrm{OP}},$ (S-17)

where

$\displaystyle
X_{\mu\nu;\xi\eta}={\sum_{\alpha}R^{\alpha}_{\mu\nu}R_{\xi\eta}^{\alpha*}}$
(S-18)

and $||X||_{\mathrm{OP}}$ is defined as the operator norm for this
$D^{2}\times D^{2}$ matrix ($||\hat{A}||_{\mathrm{\infty}}$ is the operator
norm for the $D\times D$ matrix). We therefore have

$\displaystyle\sum_{|\lambda_{\alpha}|>\lambda}|\mathrm{Tr}[\hat{O}\hat{R}_{\alpha}]|^{2}$
$\displaystyle\leq\sum_{|\lambda_{\alpha}|>\lambda}\frac{|\mathrm{Tr}[\mathcal{L}^{{\dagger}k}[\hat{O}]\hat{R}_{\alpha}]|^{2}}{|\lambda|^{2k}}$
$\displaystyle\leq\sum_{\alpha}\frac{|\mathrm{Tr}[\mathcal{L}^{{\dagger}k}[\hat{O}]\hat{R}_{\alpha}]|^{2}}{|\lambda|^{2k}}$
$\displaystyle\leq
D||X||_{\mathrm{OP}}\frac{||\mathcal{L}^{{\dagger}k}[\hat{O}]||^{2}_{\infty}}{|\lambda|^{2k}}.$
(S-19)

Next, if $\mathcal{L}$ is a local super-operator and $\hat{O}$ is a local
operator, we can show that $||\mathcal{L}^{{\dagger}k}[\hat{O}]||_{\infty}$
can be bounded from above, using a technique similar to the case of the
Hamiltonian formalism [11]. For simplicity, we assume that $\mathcal{L}$ and
$\hat{O}$ are at most $s$-local. Then we have

$\displaystyle||\mathcal{L}^{{\dagger}k}[\hat{O}]||_{\infty}\leq
o_{l}(2\ell_{l}ks)^{k},$ (S-20)

where $o_{l}$ and $\ell_{l}$ represent the norms of the local operators, which
are independent of the system size. Then, by choosing $k=\lambda/2\ell_{l}se$,
we have

$\displaystyle\sum_{|\lambda_{\alpha}|>\lambda}|\mathrm{Tr}[\hat{O}\hat{R}_{\alpha}]|^{2}$
$\displaystyle\leq
Do_{l}^{2}e^{-\frac{\lambda}{\ell_{l}se}}||X||_{\mathrm{OP}}.$ (S-21)

Note that, for the Hermitian case,
$\hat{R}_{\alpha}=\ket{E_{a}}\bra{E_{b}}\>(\alpha=(a,b))$ and thus
$X_{\mu\nu;\xi\eta}=\delta_{\mu\xi}\delta_{\nu\eta}$, and hence
$||X||_{\mathrm{OP}}=1$. However, owing to non-Hermiticity,
$||X||_{\mathrm{OP}}$ can be larger; if we assume that all the $|R_{\alpha})$
are almost parallel, we have $||X||_{\mathrm{OP}}\sim D^{2}$. In this case,
the right-hand side of Eq. (S-21) becomes large, and the bound may be trivial.
On the other hand, we numerically find that $||X||_{\mathrm{OP}}$ does not
increase faster than $\propto{D}$ for our models. Thus, assuming this scaling
and noting that the number of elements of the order of $D^{2}$ contributes to
the sum in the left-hand side, our rough estimation becomes

$\displaystyle|\mathrm{Tr}[\hat{O}\hat{R}_{\alpha}]|\lesssim
o_{l}e^{-\frac{|\lambda_{\alpha}|}{2\ell_{l}se}}.$ (S-22)

Due to the presence of the factor
$e^{-\frac{|\lambda_{\alpha}|}{2\ell_{l}se}}$ on the right-hand side,
$|\mathrm{Tr}[\hat{O}\hat{R}_{\alpha}]|$ is suppressed in the middle of the
spectrum for our models, where $\mathrm{Re}[\lambda_{\alpha}]$ is proportional
to the volume of the system. Thus, eigenmodes in the middle of the spectrum do
not contribute to the dynamics of local observables. Note that this is similar
to the results in Refs. [12, 13], which focus on the local random Lindbladian
without Hamiltonian terms.

###  II.4 Similar hypothesis for classical dynamics

Here, we briefly discuss that eigenstates of classical stochastic systems
(i.e., eigenstates of a Markovian generator) can be described by random
matrices as in the Lindblad case. For a general classical Markovian dynamics,
we consider a $D\times D$ non-Hermitian matrix $M$ as a generator of the
dynamics. We then obtain $\vec{p}(t)=e^{Mt}\vec{p}(0)$, where
$\vec{p}=(p_{1},\cdots,p_{D})^{\sf{T}}$ is the probability vector in the state
space. By the spectral decomposition

$\displaystyle
M=\sum_{a}\lambda_{a}\frac{\vec{R}_{a}\vec{L}_{a}^{\dagger}}{(\vec{R}_{a},\vec{L}_{a})},$
(S-23)

where $(\vec{a},\vec{b})$ denotes the inner product and
$\mathrm{Re}[\lambda_{a}]\leq 0$, the time evolution for an observable
$\vec{\mathcal{O}}$ can be written as

$\displaystyle\braket{\mathcal{O}(t)}=(\vec{\mathcal{O}},\vec{p}(t))=\sum_{a}e^{\lambda_{a}t}\frac{(\vec{\mathcal{O}},\vec{R}_{a})(\vec{L}_{a},\vec{p})}{(\vec{R}_{a},\vec{L}_{a})}.$
(S-24)

For a sufficiently complex dynamics, we conjecture that assumptions similar to
the Lindblad case hold true, namely

$\displaystyle(\vec{\mathcal{O}},\vec{R}_{a})\sim\mathcal{A}_{\vec{\mathcal{O}}}(\lambda_{a})r_{\alpha},$
(S-25)
$\displaystyle(\vec{\mathcal{O}},\vec{R}_{a})\sim\mathcal{B}_{\vec{p}}(\lambda_{a})r_{\alpha}$
(S-26)

and

$\displaystyle(\vec{L}_{a},\vec{R}_{a})^{-1}\sim\mathcal{C}(\lambda_{a})r^{\prime}_{\alpha}.$
(S-27)

We note that $r_{\alpha}$ and $r_{\alpha}^{\prime}$ obey non-Hermitian random
matrix theory but $\mathcal{C}(\lambda_{a})$ will be proportional to
$\sqrt{D}$ rather than $D$. In addition, employing a derivation similar to
that for Eq. (S-22) in the Lindblad case, we conjecture that for local $M$ and
$\vec{\mathcal{O}}$,

$\displaystyle\mathcal{A}_{\vec{\mathcal{O}}}(\lambda_{a})\propto
e^{-\mathrm{O}(|\lambda_{a}|)}$ (S-28)

for a locally interacting classical Markovian generator $M$.

##  III Phenomenology of the Lindbladian MBL phase

###  III.1 Structure of the Lindblad generator

To analyze the effects of dissipation in the MBL phase, we consider the matrix
representation of the Lindbladian as in Appendix II. In accordance with the
vector representation of the matrix in Eq. (S-11), the Lindbladian super-
operator is mapped to an operator

$\displaystyle\tilde{\mathcal{L}}=-i(\hat{H}\otimes\mathbb{\hat{I}}-\hat{\mathbb{I}}\otimes\hat{H}^{T})+\sum_{l}\left(\hat{\Gamma}_{l}\otimes\hat{\Gamma}_{l}^{*}-\frac{1}{2}\hat{\Gamma}_{l}^{\dagger}\hat{\Gamma}_{l}\otimes\hat{\mathbb{I}}-\frac{1}{2}\hat{\mathbb{I}}\otimes\hat{\Gamma}_{l}^{T}\hat{\Gamma}_{l}^{*}\right).$
(S-29)

In the fully MBL phase, the Hamiltonian part can be diagonalized by a quasi-
local unitary operator $U$ as [14]

$\displaystyle\hat{H}=\sum_{l}\tilde{h}_{l}\hat{\tau}_{l}^{z}+\sum_{lm}J_{lm}\hat{\tau}_{l}^{z}\hat{\tau}_{m}^{z}+\sum_{lmn}J_{lmn}\hat{\tau}_{l}^{z}\hat{\tau}_{m}^{z}\hat{\tau}_{n}^{z}+\cdots,$
(S-30)

where the coupling constants $J_{lm},J_{lmn},\cdots$ decay exponentially as
the distance between any two sites increases. The local integral of motion
$\hat{\tau}_{l}^{z}=U\hat{\sigma}_{l}^{z}U^{\dagger}$ has a large overlap with
$\hat{\sigma}_{l}^{z}$:

$\displaystyle\hat{\tau}_{l}^{z}=\sum_{\alpha=x,y,z}B_{l}^{\alpha}\hat{\sigma}_{l}^{\alpha}+\sum_{jk}\sum_{\alpha,\beta=x,y,z}D_{l,jk}^{\alpha\beta}\hat{\sigma}_{j}^{\alpha}\hat{\sigma}_{k}^{\beta}+\text{
(higher order terms)},$ (S-31)

where $D_{l,jk}^{\alpha\beta}$ decays exponentially with respect to the
distances $|l-j|$ and $|l-k|$. For the case of the transverse- and
longitudinal-field Ising model, the first-order Schrieffer-Wolff
transformation [15, 16] suggests that $B_{l}^{z}=1-\mathrm{O}((g/h)^{2})$,
$B_{l}^{x,y}=\mathrm{O}(g/h)$ and
$D_{l,jk}^{\alpha\beta}=\mathrm{O}((g/h)^{2})$, if the rare resonant sites are
neglected. Conversely, we can expand
$\hat{\sigma}_{l}^{z}=U^{\dagger}\hat{\tau}_{l}^{z}U$ with the product of
$\hat{\tau}_{l}^{\alpha}$, where
$\hat{\tau}_{l}^{x}=U\hat{\sigma}_{l}^{x}U^{\dagger}$ and
$\hat{\tau}_{l}^{y}=U\hat{\sigma}_{l}^{y}U^{\dagger}$, as

$\displaystyle\hat{\sigma}_{l}^{z}=\sum_{\alpha=x,y,z}Z_{l}^{\alpha}\hat{\tau}_{l}^{\alpha}+\sum_{jk}\sum_{\alpha,\beta=x,y,z}G_{l,jk}^{\alpha\beta}\hat{\tau}_{j}^{\alpha}\hat{\tau}_{k}^{\beta}+\text{(higher
order terms)},$ (S-32)

where $G_{l,jk}^{\alpha\beta}$ decays exponentially with respect to the
distances $|l-j|$ and $|l-k|$, $Z_{l}^{z}=1-\mathrm{O}((g/h)^{2})$,
$Z_{l}^{x,y}=\mathrm{O}(g/h)$, and
$G_{l,jk}^{\alpha\beta}=\mathrm{O}((g/h)^{2})$. A similar expression holds for
$\hat{\sigma}_{l}^{-}$.

###  III.2 Case of dephasing

Let us first consider the case of dephasing. We assume that the dissipation
strength $\gamma$ is not so large. We can decompose the matrix representation
of the Lindbladian as

$\displaystyle\tilde{\mathcal{L}}=-i(\hat{H}\otimes\mathbb{\hat{I}}-\hat{\mathbb{I}}\otimes\hat{H}^{T})+\sum_{l}({Z_{l}^{z}})^{2}\gamma(\hat{\tau}_{l}^{z}\otimes\hat{\tau}_{l}^{z}-\hat{\mathbb{I}}\otimes\hat{\mathbb{I}})+\tilde{\mathcal{L}}_{P},$
(S-33)

where $\tilde{\mathcal{L}}_{P}$ is a perturbation term of the order of $g/h$,
which arises from the second- or higher-order terms in Eq. (S-32). Without
$\tilde{\mathcal{L}}_{P}$, the eigenstates are written as

$\displaystyle|{R}_{\alpha})=|{L}_{\alpha})=\ket{\tau_{1}\cdots\tau_{L}}\otimes\ket{\tau^{\prime}_{1}\cdots\tau^{\prime}_{L}}=|\phi^{\alpha_{1}}_{1})|\phi^{\alpha_{2}}_{2})\cdots|\phi^{\alpha_{L}}_{L}).$
(S-34)

Here,
$|\phi^{\alpha_{l}}_{l})=|\tau_{l},\tau_{l}^{\prime})=\ket{\tau_{l}}\otimes\ket{\tau_{l}^{\prime}}$
is the vector representation of
$\hat{\phi}_{l}^{\alpha_{l}}=\ket{\tau_{l}}\bra{\tau_{l}^{\prime}}$, where
$\ket{\tau_{l}}$ is the eigenstate of $\hat{\tau}_{l}^{z}$ with eigenvalue
$\tau_{l}\>(=\pm 1)$. Specifically, we take $|\phi^{1})=|+1,-1)$,
$|\phi^{2})=|-1,+1)$, $|\phi^{3})=|+1,+1)$, and $|\phi^{4})=|-1,-1)$. The
corresponding eigenvalue of $|R_{\alpha})$ is given by

$\displaystyle\lambda_{\alpha}=$
$\displaystyle-i\left(\sum_{l}\tilde{h}_{l}{\tau}_{l}+\sum_{lm}J_{lm}{\tau}_{l}{\tau}_{m}+\cdots\right)+i\left(\sum_{l}\tilde{h}_{l}{\tau}_{l}^{\prime}+\sum_{lm}J_{lm}{\tau}_{l}^{\prime}{\tau}_{m}^{\prime}+\cdots\right)$
$\displaystyle+\sum_{l}({Z_{l}^{z}})^{2}\gamma({\tau}_{l}{\tau}_{l}^{\prime}-1),$
(S-35)

which can also be written as

$\displaystyle\lambda_{\alpha}=-i\sum_{l}2\tilde{h}_{l}(\delta_{\alpha_{l},1}-\delta_{\alpha_{l},2})-i(\text{small
terms including
}J_{lm},\cdots)-\sum_{l}2{Z_{l}^{z}}^{2}\gamma(\delta_{\alpha_{l},1}+\delta_{\alpha_{l},2}).$
(S-36)

Let us consider the effect of $\tilde{\mathcal{L}}_{P}$. This perturbation
includes the bit flip terms, i.e., terms that can flip $\ket{+1}\>(\ket{-1})$
(or its product) to $\ket{-1}\>(\ket{+1})$ (or its product). To simplify the
notation, we write
$|\phi^{\alpha_{1}}_{1})|\phi^{\alpha_{2}}_{2})\cdots|\phi^{\alpha_{L}}_{L})$
as $|\alpha_{1}\alpha_{2}\cdots\alpha_{L})$. First, states such as
$|d_{1}\cdots d_{L})\>(d_{l}=3,4)$ [$d$ stands for diagonal degree of freedom]
are unstable under perturbation and mixed with the other states in the form
$|d^{\prime}_{1}\cdots d^{\prime}_{L})\>(d_{l}^{\prime}=3,4)$, since the
energy difference between these states is zero. Note that these states
correspond to the long-time longitudinal relaxation of a local integral of
motion, which can be treated by an effective classical rate equation [1, 3].
Next, consider a state that contains a large cluster of $d_{l}$ with size
$c\gg 1$, e.g., $[d_{i}\cdots d_{k}\cdots d_{i+c-1}]\>(d_{l}=3,4)$ with
$|i-k|,|i+c-1-k|\gg 1$. In this case, $d_{k}$ and nearby bits can be flipped
to another $d_{k}^{\prime}$ by perturbation because the eigenvalue difference
between these states is $\sim\max\\{J_{i,k},J_{k,i+c-1}\\}$, which is
exponentially small. Note that $d_{l}\>(d_{l}=3,4)$ corresponds the diagonal
elements $\ket{+1}\bra{+1}$ or $\ket{-1}\bra{-1}$ in the original density
matrix before vectorization. In this sense, some eigenstates exist for which
the diagonal degrees of freedom form large delocalized clusters.

On the other hand, the states $|\phi^{1}_{l})$ and $|\phi^{2}_{l})$, which
correspond to the off-diagonal degrees of freedom in the density matrix, are
typically stable under (first-order) perturbation. Here, we explain some of
the prototypical examples for which the transition processes are suppressed by
strong disorder. For the single flip from $o_{l}\>(o_{l}=1,2)$ [$o$ stands for
off-diagonal degree of freedom] to $d_{l}\>(d_{l}=3,4)$, while the transition
matrix element is $\sim\gamma(Z_{l}^{x/y})\sim g\gamma/\tilde{h}_{l}$, the
eigenvalue difference becomes $\sim|\pm
2i\tilde{h}_{l}+2(Z_{l}^{z})^{2}\gamma|$, so the ratio becomes $\sim
g\gamma/h^{2}\ll 1$. According to perturbation theory, this indicates that the
transition process is suppressed due to the large energy mismatch. For the
single flip from $o_{l}\>(o_{l}=1,2)$ to
$o^{\prime}_{l}\>(o^{\prime}_{l}=2,1)$ with the transition matrix element
$\sim\gamma{(Z_{l}^{x/y})}^{2}\sim g^{2}\gamma/\tilde{h}_{l}^{2}$, the
eigenvalue difference becomes $\sim|4i\tilde{h}_{l}|$, and therefore the ratio
again becomes small, of the order of $g^{2}\gamma/h^{3}\ll 1$. For another
example of the double flip from $|\cdots o_{l}\cdots d_{m}\cdots)$ to $|\cdots
d_{l}\cdots o_{m}\cdots)$ with the transition matrix element
$\mathrm{O}(g^{2}\gamma/{h}^{2})$ and exponentially small with respect to
$|l-m|$, the eigenvalue difference becomes $\sim|\pm
2i(\tilde{h}_{l}-\tilde{h}_{m})+2\\{(Z_{l}^{z})^{2}\gamma-(Z_{m}^{z})^{2}\gamma\\}|$,
and therefore the ratio is again small unless $\tilde{h}_{l}-\tilde{h}_{m}$ is
accidentally small.

In conclusion, for sufficiently strong disorder, the off-diagonal degrees of
freedom can be localized. The decay rate is then stabilized at

$\displaystyle\mathrm{Im}[\lambda_{\alpha}]=\sum_{l}2\gamma(\delta_{\alpha_{l},1}+\delta_{\alpha_{l},2})+\mathrm{O}\left(\left(\frac{g}{h}\right)^{2}\right).$
(S-37)

If there are $q$ localized off-diagonal degrees of freedom, the decay rate is
$\sim 2\gamma q$, as discussed in the main text.

The localization/delocalization of off-diagonal/diagonal degrees of freedom is
intuitively understood from the matrix representation of the Lindbladian. As
shown in Fig. S-6, the Lindbladian includes a term
$-ih_{l}(\hat{\tau}_{l}^{z}\otimes\hat{\mathbb{{I}}}-\hat{\mathbb{{I}}}\otimes\hat{\tau}_{l}^{z})$.
For off-diagonal degrees of freedom $|\phi^{1})$ and $|\phi^{2})$, they
experience net disorder $\mp 2ih_{l}$ and thus localize. On the other hand,
diagonal degrees of freedom $|\phi^{3})$ and $|\phi^{4})$ do not undergo
disorder and thus delocalize.

Figure S-6:  Schematic illustration of the Lindbladian in the matrix
representation. Black arrows show quasi-local bits in the $\hat{\tau}^{z}$
basis. For diagonal degrees of freedom (surrounded by dotted lines), the net
disorder vanishes and delocalization occurs. On the other hand, off-diagonal
degrees of freedom (surrounded by solid lines) undergo finite net disorder
$\mp 2ih_{l}$, and localization occurs.

To discuss the operator-space entanglement entropy (OSEE), we notice that
eigenstates in the middle of the spectrum have $\sim(1-\alpha)L\>\>(\alpha<1)$
localized off-diagonal degrees of freedom. Thus, the number of the remaining
diagonal degrees of freedom in those states is $\sim\alpha L$. Since the
localized degrees of freedom do not contribute to the OSEE (except for the
exponential tail coming from the difference between $\hat{\tau}_{l}$ and the
physical bits $\hat{\sigma}_{l}$), the half-chain OSEE for strong disorder is
typically bounded from above by $\sim\alpha L\log 2$. This is generally
smaller than the half-chain OSEE for the weak-disorder case, where the RMT
prediction leads to $\sim L\log 2$. We also note that, if the clusters of the
diagonal degrees of freedom are small and kept far apart from one another by
many localized off-diagonal clusters, they cannot create much entanglement. In
that case, the OSEE can be further reduced from $\sim\alpha L\log 2$ and
satisfy the area law $\sim\mathrm{O}(L^{0})$. Thus, depending on the number
and the position of the localized off-diagonal degrees of freedom, the OSEE
can satisfy either the volume law $\sim\alpha L\log 2$ or the area law
$\sim\mathrm{O}(L^{0})$. In the middle of the spectrum, we expect many
localized off-diagonal clusters and the area-law behavior of the OSEE. This is
consistent with the numerical simulation in the main text, but we leave it for
a future problem to investigate the precise evaluation for larger system
sizes. In any case, at the transition between these two phases, there appears
a large fluctuation of the OSEE for different eigenstates, leading to the
observed peak in the main text.

We also note that, even when a large delocalized cluster of diagonal degrees
of freedom exists, it may not have the power to delocalize the localized off-
diagonal degrees of freedom, which is distinct from the situation where the
large bath can delocalize the localized bits in the Hermitian disordered
systems [17]. To see this, let us consider a simple situation where a single
localized off-diagonal bit, $|\phi_{1}^{1})$ (or $|\phi_{1}^{2})$), is coupled
to a delocalized cluster of diagonal degrees of freedom $|D)$. For the
resonant transition from $|\phi_{1}^{1})\otimes|D)$ to
$|\phi_{1}^{\alpha})\otimes|D^{\prime})\>(\alpha=2,3,4)$ to occur under some
small perturbation, where $|D^{\prime})$ is another state for the delocalized
cluster, the difference between the imaginary parts of the energies of $D$ and
$D^{\prime}$ should be of the order of $h$. However, since even large
delocalized clusters have almost zero imaginary parts of energies, which
correspond to the eigenvalues in the original language, this resonance
condition is typically not satisfied.

###  III.3 Case of damping

We briefly discuss the damping case. To do this, we decompose the Hamiltonian
as $\hat{H}=\hat{H}_{NI}+\hat{H}_{I}$, where
$\hat{H}_{NI}=\sum_{l}\tilde{h}_{l}\hat{\tau}_{l}^{z}$ is the non-interacting
term of $\hat{\tau}_{l}^{z}$ and $\hat{H}_{I}$ is the interaction term. Then,
we have

$\displaystyle\tilde{\mathcal{L}}=-i(\hat{H}_{NI}\otimes\mathbb{\hat{I}}-\hat{\mathbb{I}}\otimes\hat{H}_{NI}^{T})+\sum_{l}\frac{{Z_{l}^{z}}^{2}\gamma}{4}(2\hat{\tau}_{l}^{-}\otimes\hat{\tau}_{l}^{-}-\hat{\tau}_{l}^{+}\hat{\tau}_{l}^{-}\otimes\hat{\mathbb{I}}-\hat{\mathbb{I}}\otimes\hat{\tau}_{l}^{+}\hat{\tau}_{l}^{-})+\tilde{\mathcal{L}}_{I}+\tilde{\mathcal{L}}_{P},$
(S-38)

where
$\tilde{\mathcal{L}}_{I}=-i(\hat{H}_{I}\otimes\mathbb{\hat{I}}-\hat{\mathbb{I}}\otimes\hat{H}_{I}^{T})$.
We first consider $\tilde{\mathcal{L}}_{I}+\tilde{\mathcal{L}}_{P}$ as a
perturbation.

The eigenstates of $\mathcal{L}-\mathcal{L}_{I}-\mathcal{L}_{P}$ can be
written as

$\displaystyle|{R}_{\alpha})=|\eta^{\alpha_{1}}_{1})|\eta^{\alpha_{2}}_{2})\cdots|\eta^{\alpha_{L}}_{L})$
(S-39)

and

$\displaystyle|{L}_{\alpha})=|\chi^{\alpha_{1}}_{1})|\chi^{\alpha_{2}}_{2})\cdots|\chi^{\alpha_{L}}_{L}),$
(S-40)

where $|\eta^{1})=|\chi^{1})=|+1,-1)$, $|\eta^{2})=|\chi^{2})=|-1,+1)$,
$|\eta^{3})=\frac{|+1,+1)-|-1,-1)}{\sqrt{2}}$, $|\chi^{3})=|+1,+1)$,
$|\eta^{4})=|-1,-1)$, and $|\chi^{4})=\frac{|+1,+1)+|-1,-1)}{\sqrt{2}}$. The
corresponding eigenvalue is given by

$\displaystyle-i\sum_{l}2\tilde{h}_{l}(\delta_{\alpha_{l},1}-\delta_{\alpha_{l},2})-i(\text{terms
including
}J_{lm},\cdots)-\sum_{l}{Z_{l}^{z}}^{2}\gamma(\delta_{\alpha_{l},1}+\delta_{\alpha_{l},2}+2\delta_{\alpha_{l},3}).$
(S-41)

Even if we add $\tilde{\mathcal{L}}_{I}+\tilde{\mathcal{L}}_{P}$, the states
$|\eta^{1}_{l})$ and $|\eta^{2}_{l})$ are robust under the first-order
perturbation. Indeed, just as the dephasing case, changing these states to
other states is typically energetically costly due to strong disorder. This
leads to the localization of these off-diagonal degrees of freedom and the
breakdown of non-Hermitian random-matrix-type universality.

We also note that, the longest-lived modes involving $q$ off-diagonal degrees
of freedom have stabilized decay rates $\sim\gamma
q(1+\mathrm{O}((g/h)^{2}))$. To see this, we next consider only
$\tilde{\mathcal{L}}_{P}$ as a perturbation as in the dephasing case. Without
perturbation, we notice that states such as

$\displaystyle|\eta^{4}_{1})|\eta^{4}_{2})\cdots|\eta^{4}_{L-q})|\eta^{1\text{
or }2}_{L-q+1})\cdots|\eta^{1\text{ or }2}_{L})$ (S-42)

are eigenstates of the unperturbed Lindbladian in the matrix representation
with the eigenvalue $\gamma q(1+\mathrm{O}((g/h)^{2}))$. Note that the
eigenstates involving $|\eta^{3})$ have a faster decay rate, so the states
represented in Eq. (S-42) are those having the longest lifetime with $q$ off-
diagonal degrees of freedom. When we introduce the perturbation
$\tilde{\mathcal{L}}_{P}$, the leading order of the change of the eigenvalue
is $\mathrm{O}(\gamma(g/h)^{2}))$, and the eigenvalue is still expressed as
$\gamma q(1+\mathrm{O}((g/h)^{2}))$.

##  IV Effects of dissipation on the Lindbladian MBL transition point

Here, we discuss how the Lindbladian MBL transition depends on the dissipation
strength $\gamma$.

Figure S-7:  Complex-spacing ratio $r$ (a) and the variance $\sigma$ of the
OSEE (b) as a function of $h$ for different values of $\gamma$. We show the
results for a dephasing type of dissipation. We find that $r$ changes from the
universal value of non-Hermitian RMT (class AI†) to the Poisson value as $h$
increases, unless $\gamma$ is too small or too large. On the other hand, the
critical value of $h$ depends on $\gamma$. For small $\gamma\>(\lesssim 0.5)$,
an increase in $\gamma$ increases the critical value of the transition, which
is read out from the peak of $\sigma$ (dotted circles in panel (b)). On the
other hand, a further increase in $\gamma$ decreases the critical value. When
$\gamma$ is large enough (e.g., $\gamma=1.75$), the peak no longer appears,
which indicates that the system belongs to the localized phase even for weak
$h$. We use the Lindbladian with $J=1,g=-0.9,L=6$ and eigenvalues/eigenstates
in the middile of the supectrum [18]. The number of samples used for the
simulation is $800$.  Figure S-8:  Complex-spacing ratio $r$ (a) and the
variance $\sigma$ of the OSEE (b) as a function of $h$ for different values of
$\gamma$. We show the results for a damping type of dissipation. We find that
$r$ changes from the universal value of non-Hermitian RMT (class A) to the
Poisson value as $h$ increases, unless $\gamma$ is too small or too large. On
the other hand, the critical value of $h$ depends on $\gamma$. For small
$\gamma\>(\lesssim 1)$, increasing $\gamma$ increases the critical value of
the transition, which is read out from the peak of $\sigma$ (dotted circles in
panel (b)). On the other hand, a further increase in $\gamma$ decreases the
critical value. When $\gamma$ is strong enough (e.g., $\gamma=5.0$), the peak
no longer appears, which indicates that the system belongs to the localized
phase even for weak $h$. We use the Lindbladian with $J=1,g=-0.9,L=6$ and
eigenvalues/eigenstates in the middile of the supectrum [18]. The number of
samples used for the simulation is $800$. Figure S-9:  (a) Schematic
illustration of the phase diagram of dissipative disordered systems. There are
two phases: a phase characterized by the non-Hermitian RMT and that
characterized by the Lindbladian MBL. (b) Phase diagram for the dephasing
case. Each transition point is determined from the maximum peak of the
variance of the OSEE for $L=6$, as shown in Fig. S-7(b). (c) Phase diagram for
the damping case. Each transition point is determined from the maximum peak of
the variance of the OSEE for $L=6$, as shown in Fig. S-8(b). While the phase
boundaries are quantitatively different for (b) and (c), both of them are
qualitatively described by the non-monotonic phase boundary given in (a).

###  IV.1 Phase diagram

We first discuss the possible phase diagram of our system. Figures S-7 and S-8
show the $h$-dependence of the complex-spacing ratio $r$ and that of the
variance $\sigma$ of the OSEE for different values of $\gamma$ for the cases
of dephasing and damping, respectively. We find that $r$ changes from the
universal value of non-Hermitian RMT to the Poisson value as $h$ increases, as
long as $\gamma$ is not too small. On the other hand, the critical value
depends on $\gamma$. For small $\gamma\>(\lesssim 0.5)$, increasing $\gamma$
increases the critical value of the transition, which is read out from the
peak of $\sigma$. On the other hand, a further increase in $\gamma$ decreases
the critical value. If $\gamma$ is strong enough, the peak no longer appears,
which indicates that the system belongs to the localized phase even for weak
$h$.

From these results, we can identify the phase diagram of our system with each
type of dissipation, as shown in Fig. S-9. Two remarks are in order here.
First, the transition points $h_{c}(\gamma)$ exhibit a non-monotonic behavior.
Second, a clear quantitative difference exists for the transition points of
dephasing-type and damping-type cases. These two points indicate a nontrivial
interplay between dissipation and disorder in this system.

Note that the phase diagrams in Fig. S-9 are obtained for relatively small
system sizes ($L=6$). When we increase the system size, the exact values of
the transition points may change. However, we conjecture that the qualitative
feature of the phase diagram remains the same for larger system sizes.

###  IV.2 Relation to the many-body localization without dissipation

Here, we discuss the relation between the Lindbladian MBL and the conventional
MBL that occurs in our Ising model $\hat{H}$ without dissipation (i.e.,
$\gamma=0$). While this short-range interacting model is believed to exhibit
MBL for strong disorder [15, 16], not much is investigated about the
transition behavior between delocalized and localized phases (see, e.g., Ref.
[19] for an exception).

Figure S-10:  Many-body localization transition in the Hermitian quantum Ising
model. (a) The spacing ratio $r_{\mathrm{H}}$ of (real) energy eigenvalues of
$\hat{H}$. It changes from $r_{\mathrm{H,AI}}\simeq 0.5307$ (the RMT
prediction corresponding to Hermitian class AI) to $r_{\mathrm{H,Po}}\simeq
0.3863$ (the Poisson value), as $h$ is increased. We use $L=11$ (green) and
$12$ (blue). (b) The variance of the entanglement entropy of energy
eigenstates. The value of the peak becomes larger as we increase the system
size from $L=6$ (red) to $12$ (blue). From the peak for $L=12$, we can
estimate the MBL transition point as $h_{\mathrm{H},c}=2.4\pm 0.2$. On the
other hand, the peak position for $L=6$ is slightly shifted to $2.2\pm 0.2$.
For both panels, we use the Hamiltonian with $J=1$ and $g=-0.9$. We use the
eigenstates in the middle of the spectrum (6.25%). The number of samples used
for the simulation is $10000\>(L=6,7,8,9,10)$ and $800\>(L=11,12)$.

We numerically investigate the MBL transition for $\hat{H}$ with the same
parameters in the Hamiltonian as those used in the main text. For this
purpose, we first calculate the spacing ratio of (real) energy eigenvalues of
$\hat{H}$ [20], defined by
$r_{\mathrm{H}}=\mathbb{E}_{a}\left[\frac{\min{(E_{a+1}-E_{a},E_{a}-E_{a-1})}}{\max{({E_{a+1}-E_{a},E_{a}-E_{a-1}})}}\right]$.
As shown in Fig. S-10(a), it changes from $r_{\mathrm{H,AI}}\simeq 0.5307$
(the RMT prediction corresponding to Hermitian class AI) to
$r_{\mathrm{H,Po}}\simeq 0.3863$ (the Poisson value), as $h$ is increased.

We next calculate the eigenstate variance $\sigma_{\mathrm{H}}$ of the
entanglement entropy
$S_{\mathrm{H},a}(\ket{E_{a}})=-\mathrm{Tr}[\hat{\rho}_{a}^{l}\ln\hat{\rho}_{a}^{l}]$
of the energy eigenstates $\ket{E_{a}}$, where $\hat{\rho}_{a}^{l}$ is the
reduced density matrix of $\ket{E_{a}}$ to a region $[1,l]\>(l=\lfloor
L/2\rfloor)$. As shown in Fig. S-10(b), $\sigma_{\mathrm{H}}$ has a peak when
$h$ is varied, from which we can read out the MBL transition point
$h_{\mathrm{H},c}$ [21]. We can estimate that the transition point as
$h_{\mathrm{H},c}=2.4\pm 0.2$ from the results up to $L=12$.

A natural question is whether the Lindbladian MBL transition point
$h_{c}(\gamma)$ coincides with the Hermitian MBL transition point
$h_{\mathrm{H},c}$ in the dissipationless limit $\gamma\rightarrow 0$. From
Fig. S-9, we can see that $h_{c}(\gamma\rightarrow 0)\leq 2.0\pm 0.2$. This is
different from the above estimate of $h_{\mathrm{H},c}=2.4\pm 0.2$ obtained
from $L=12$. On the other hand, since $h_{c}(\gamma)$ is numerically estimated
from $L=6$, it is more natural to compare the result with $h_{\mathrm{H},c}$
obtained from the data for $L=6$ 111Note that the true transition point should
be determined from the limit of $L\rightarrow\infty$. We here compare the
estimation of the transition points for the same system size to take account
of the degree of the finite-size effects appropriately.. From Fig. S-10(b),
this corresponds to $h_{\mathrm{H},c}=2.2\pm 0.2$, which is consistent with
the value of $h_{c}(\gamma)$ with small $\gamma$.

To discuss behavior in the thermodynamic limit, let us assume that the OSEE of
the Lindbladian eigenstates is continuous as a function of $\gamma$. If that
is the case, the OSEE can be obtained as

$\displaystyle\lim_{\gamma\rightarrow
0}S_{\alpha}(\gamma)=S_{\alpha}(\gamma=0)=S_{\mathrm{H},a}+S_{\mathrm{H},b},$
(S-43)

where we have assumed $\hat{R}_{\alpha}(\gamma=0)=\ket{E_{a}}\bra{E_{b}}$ with
$\lambda_{\alpha}=-i(E_{a}-E_{b})$. For simplicity, we focus on the
statistical behavior of $S_{\alpha}(\gamma=0)$ for eigenstates satisfying
$-i\lambda_{\alpha}\sim\omega\in\mathbb{R}$. The total number of such
eigenstates becomes $\int dEe^{\mathcal{S}(E)+\mathcal{S}(E+\omega)}$, where
$\mathcal{S}(E)$ is the thermodynamic entropy of the system. While this number
contains multiple energy scales, for large system sizes, the contribution
coming from $E=\mathcal{E}$ that satisfies

$\displaystyle\frac{\mathrm{d}\mathcal{S}}{\mathrm{d}E}(\mathcal{E})+\frac{\mathrm{d}\mathcal{S}}{\mathrm{d}E}(\mathcal{E}+\omega)=0$
(S-44)

dominates within the saddle-point approximation. Then, $S_{\alpha}(\gamma=0)$
for almost all $\alpha$ behaves as $S_{\mathrm{H},a}+S_{\mathrm{H},b}$ with
$E_{a}\simeq\mathcal{E}$ and $E_{b}\simeq\mathcal{E}+\omega$. For example, the
average of $S_{\alpha}(\gamma=0)$ over eigenstates satisfies
$S(\lambda=i\omega)\simeq
S_{\mathrm{H}}(\mathcal{E})+S_{\mathrm{H}}(\mathcal{E}+\omega)$, and
similarly, the variance satisfies
$\sigma(\lambda=i\omega)\simeq\sigma_{\mathrm{H}}(\mathcal{E})+\sigma_{\mathrm{H}}(\mathcal{E}+\omega)$,
provided that there is no correlation between fluctuations of
$S_{\mathrm{H},a}$ and $S_{\mathrm{H},b}$.

From the above discussion, we can relate $h_{c}(\gamma\rightarrow 0)$ to
$h_{\mathrm{H},c}$ if we further make either of the following two assumptions:
(i) $\omega=\mathrm{o}(L)$, or (ii) there is no mobility edge, i.e., the MBL
transition occurs for all energy scales simultaneously. If we assume (i), we
have $S(\lambda=i\omega)\simeq 2S_{\mathrm{H}}(\mathcal{E})$ and
$\sigma(\lambda=i\omega)\simeq 2\sigma_{\mathrm{H}}(\mathcal{E})$. We also
have $\frac{\mathrm{d}\mathcal{S}}{\mathrm{d}E}(\mathcal{E})=0$, indicating
that $\mathcal{E}$ corresponds to the energy whose entropy is maximal. Then,
the MBL transition (i.e., the transition that $S_{\mathrm{H}}$ changes to the
area law and that $\sigma_{\mathrm{H}}$ has a peak) at $\mathcal{E}$ leads to
the Lindbladian MBL. In other words, $h_{\mathrm{H},c}$ at infinite
temperature and $h_{c}(\gamma\rightarrow 0)$ with small $\omega=\mathrm{o}(L)$
coincide in the thermodynamic limit. Instead, if we assume (ii), whose
possibility is discussed in Ref. [23], the MBL transition point
$h_{\mathrm{H},c}$ does not depend on $E$. Thus, in crossing
$h_{\mathrm{H},c}$, $S_{\mathrm{H}}(\mathcal{E})$ and
$S_{\mathrm{H}}(\mathcal{E}+\omega)$ change to the area law, and
$\sigma_{\mathrm{H}}(\mathcal{E})$ and
$\sigma_{\mathrm{H}}(\mathcal{E}+\omega)$ have a peak. This directly leads to
the change from the volume to the area laws for $S(\lambda)$ and the emergence
of a peak for $\sigma(\lambda)$. In other words, $h_{\mathrm{H},c}$ and
$h_{c}(\gamma\rightarrow 0)$ coincide in the thermodynamic limit. It is a
future problem to numerically investigate if either (or both) of the above two
scenarios holds for large system sizes.

## References

  * Fischer  _et al._ [2016] M. H. Fischer, M. Maksymenko, and E. Altman, Phys. Rev. Lett. 116, 160401 (2016).
  * Levi  _et al._ [2016] E. Levi, M. Heyl, I. Lesanovsky, and J. P. Garrahan, Phys. Rev. Lett. 116, 237203 (2016).
  * Medvedyeva  _et al._ [2016] M. V. Medvedyeva, T. Prosen, and M. Znidaric, Phys. Rev. B 93, 094205 (2016).
  * Breuer  _et al._ [2002] H.-P. Breuer, F. Petruccione, _et al._ ,  _The theory of open quantum systems_ (Oxford University Press on Demand, 2002).
  * Chalker and Mehlig [1998] J. T. Chalker and B. Mehlig, Phys. Rev. Lett. 81, 3367 (1998).
  * Bourgade and Dubach [2019] P. Bourgade and G. Dubach, Probability Theory and Related Fields , 1 (2019).
  * Mehlig and Chalker [2000] B. Mehlig and J. T. Chalker, Journal of Mathematical Physics 41, 3233 (2000).
  * Janik  _et al._ [1999] R. A. Janik, W. Nörenberg, M. A. Nowak, G. Papp, and I. Zahed, Phys. Rev. E 60, 2699 (1999).
  * Fyodorov [2018] Y. V. Fyodorov, Communications in Mathematical Physics 363, 579 (2018).
  * Yoshioka and Hamazaki [2019] N. Yoshioka and R. Hamazaki, Phys. Rev. B 99, 214306 (2019).
  * Mori  _et al._ [2016] T. Mori, T. Kuwahara, and K. Saito, Phys. Rev. Lett. 116, 120401 (2016).
  * Wang  _et al._ [2020] K. Wang, F. Piazza, and D. J. Luitz, Phys. Rev. Lett. 124, 100604 (2020).
  * Sommer  _et al._ [2021] O. E. Sommer, F. Piazza, and D. J. Luitz, Phys. Rev. Research 3, 023190 (2021).
  * Huse  _et al._ [2014] D. A. Huse, R. Nandkishore, and V. Oganesyan, Phys. Rev. B 90, 174202 (2014).
  * Imbrie [2016a] J. Z. Imbrie, Phys. Rev. Lett. 117, 027201 (2016a).
  * Imbrie [2016b] J. Z. Imbrie, Journal of Statistical Physics 163, 998 (2016b).
  * De Roeck and Huveneers [2017] W. De Roeck and F. m. c. Huveneers, Phys. Rev. B 95, 155129 (2017).
  * [18] Specifically, let us consider the real part $\lambda_{\alpha}^{r}$ and the absolute value $|\lambda_{\alpha}^{i}|$ of the imaginary part of $\lambda_{\alpha}$, and sort each of them independently as $a_{1}\leq\cdots\leq a_{D^{2}}=0$ and $0\leq b_{1}\leq\cdots\leq b_{D^{2}}$. We use only those modes whose eigenvalues satisfy $a_{3D^{2}/8}\leq\lambda_{\alpha}^{r}\leq a_{5D^{2}/8}$ and $b_{3D^{2}/8}\leq|\lambda_{\alpha}^{i}|\leq b_{5D^{2}/8}$.
  * Yao  _et al._ [2020] Z. Yao, C. Liu, P. Zhang, and H. Zhai, Phys. Rev. B 102, 104302 (2020).
  * Atas  _et al._ [2013] Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux, Phys. Rev. Lett. 110, 084101 (2013).
  * Kjäll  _et al._ [2014] J. A. Kjäll, J. H. Bardarson, and F. Pollmann, Phys. Rev. Lett. 113, 107204 (2014).
  * Note [1] Note that the true transition point should be determined from the limit of $L\rightarrow\infty$. We here compare the estimation of the transition points for the same system size to take account of the degree of the finite-size effects appropriately.
  * De Roeck  _et al._ [2016] W. De Roeck, F. Huveneers, M. Müller, and M. Schiulaz, Phys. Rev. B 93, 014203 (2016).


# Absence of jet quenching
in peripheral nucleus–nucleus collisions

Constantin Loizides
LBNL, Berkeley, CA, 94720, USA
Andreas Morsch
CERN, 1211 Geneva 23, Switzerland

###### Abstract

Medium effects on the production of high-$p_{\rm T}$ particles in
nucleus–nucleus (AA) collisions are generally quantified by the nuclear
modification factor ($R_{\rm AA}$), defined to be unity in absence of nuclear
effects. Modeling particle production including a nucleon–nucleon impact
parameter dependence, we demonstrate that $R_{\rm AA}$ at midrapidity in
peripheral AA collisions can be significantly affected by event selection and
geometry biases. Even without jet quenching and shadowing, these biases cause
an apparent suppression for $R_{\rm AA}$ in peripheral collisions, and are
relevant for all types of hard probes and all collision energies. Our studies
indicate that calculations of jet quenching in peripheral AA collisions should
account for the biases, or else they will overestimate the relevance of parton
energy loss. Similarly, expectations of parton energy loss in light–heavy
collision systems based on comparison with apparent suppression seen in
peripheral $R_{\rm AA}$ should be revised. Our interpretation of the
peripheral $R_{\rm AA}$ data would unify observations for lighter collision
systems or lower energies where significant values of elliptic flow are
observed despite the absence of strong jet quenching.

Medium effects on the production of high-$p_{\rm T}$ particles are in general
quantified by the nuclear modification factor

$R_{\rm AA}=\frac{Y_{\rm AA}}{N_{\rm coll}Y_{\rm pp}}=\frac{Y_{\rm AA}}{T_{\rm
AA}\ \sigma_{\rm pp}}$ (1)

defined as the ratio of the per-event yield $Y_{\rm AA}$ measured in
nucleus–nucleus (AA) collisions to the yield of an equivalent incoherent
superposition of $N_{\rm coll}$ binary pp collisions. The number of binary
collisions depends on the overlap between the two colliding nuclei quantified
by the nuclear overlap $T_{\rm AA}$. It is expected that in the absence of
nuclear effects $R_{\rm AA}$ is unity. However, strictly speaking this holds
only for centrality integrated measurements. In this case $N_{\rm coll}$ is
given by $N_{\rm coll}=A^{2}\sigma_{\rm pp}/\sigma_{\rm AA}$, where
$\sigma_{\rm pp}$ and $\sigma_{\rm AA}$ are, respectively, the pp and AA
inelastic cross-sections. As will be outlined in the following in more detail,
centrality classification can lead to the selection of AA event samples for
which the properties of the binary NN collisions deviate from unbiased pp
collisions. In this case $R_{\rm AA}$ can deviate from unity even in the
absence of nuclear effects. There are two main origins for selection biases.
Firstly, the spatial distribution of nucleons bound in nuclei in the plane
transverse to the beam direction differs from those of protons in a beam
leading to a bias on the NN impact parameter. Secondly, centrality selection
is based on measurements related to bulk, soft particle production. The
selection can bias the mean multiplicity of individual NN collisions and in
case of a correlations between soft and hard particle production the yield of
hard processes in AA collisions.

In the optical Glauber model [1], the nuclear overlap is obtained from the
nuclear density distributions and the impact parameter $b$ between two nuclei,
which is the only parameter characterising a collision. Instead, Monte Carlo
(MC) Glauber models [1] take into account collision-by-collision fluctuations
at fixed impact parameter by allowing the positions of the nucleons in the
nuclei to vary. The number of binary collisions is obtained by assuming that
the nucleons move on straight trajectories and a collision is counted if the
nucleon–nucleon (NN) impact parameter $b_{\rm NN}$ is below a certain
threshold (usually given by the inelastic NN cross section). For each
simulated AA collision the MC Glauber determines $N_{\rm coll}$, and for each
of the $N_{\rm coll}$ nucleon–nucleon collisions the impact parameter $b_{\rm
NN}^{i}$ and the respective collision position $(x^{i},y^{i})$ in the
transverse plane. In this way, such calculations provide important information
about the energy density distribution including its event-by-event
fluctuations in the initial state of AA collisions, which can be used as input
for hydrodynamic calculations. However, for the evaluation of the nuclear
modification factor the information about the individual NN collisions is
usually ignored. An impact parameter dependent NN profile can also be enabled
in the GLISSANDO model [2], but is not (yet) available in the widely-used
standard Glauber MC [3, 4].

In variance to the standard MC Glauber approach, the HIJING model [5] takes
into account the possibility of multiple hard scatterings (multiple parton
interactions) in the same NN collision. As in MC models for pp collisions like
for example PYTHIA [6], the mean number of hard scatterings per collison
depends on the NN impact parameter. While the NN collisions are still modeled
as incoherent, the production rate of hard processes is not proportional to
$N_{\rm coll}$ but to

$N_{\rm hard}=N_{\rm coll}\cdot\left.N^{\rm hard}_{\rm NN}\right|_{\rm
C}/\left<N^{\rm hard}_{\rm NN}\right>\,,$ (2)

where $\left.N^{\rm hard}_{\rm NN}\right|_{\rm C}$ is the average number of
hard scatterings in a NN collision for a given centrality selection and
$\left<N^{\rm hard}_{\rm NN}\right>$ is its unbiased average value. Similarly,
LABEL:Jia:2009mq describes an extension of the optical Glauber model, in which
the nuclear overlap function is obtained from a convolution between the
product of the thickness functions of the two nuclei and the nucleon–nucleon
overlap function.

These extensions have important consequences for the AA impact parameter
dependence of hard processes. With respect to standard Glauber $N_{\rm coll}$
scaling, the number of hard processes is suppressed in peripheral collisions
due to a simple geometrical bias. With increasing AA impact parameter the
phase space for collisions increases $\propto b$ whereas the nuclear density
decreases leading to an increased probability for more peripheral than average
NN collisions.

A further consequence arises if the yield of hard and soft processes are
correlated via the common $b_{\rm NN}$ and centrality selection is based on
soft particle production (multiplicity or summed energy) [8, 9]. In this case
for a given centrality class, the NN collisions can be biased towards higher
or lower than average impact parameters. The event selection bias is in
particular important when fluctuations of the centrality estimator due to
$b_{\rm NN}$ are of similar size as the dynamic range of $N_{\rm coll}$, as in
pA collisions.

In contrast, centrality measurements based on zero-degree energy should not
introduce any selection bias, while the geometric bias could still play a
role. In the so called hybrid method, described in LABEL:Adam:2014qja, the pPb
centrality selection is based on zero-degree neutral energy in the Pb-going
directions (slow neutrons) and $N_{\rm coll}$ is determined from the measured
charged particle multiplicity $M$ according to $N_{\rm coll}=\langle N_{\rm
coll}\rangle\cdot M/\langle M\rangle$, where $\langle N_{\rm coll}\rangle$ and
$\langle M\rangle$ are, respectively, the centrality averaged number of
collisions and multiplicty. In case soft and hard particle yields are affected
in the same way, the selection bias would cancel in the nuclear modification
factor.

In peripheral AA collisions, one can expect the selection bias to be relevant,
in addition to the geometric bias, as demonstrated in figure 8 of [9], which
quantifies the ratio between the average multiplicity of the centrality
estimator and the average multiplicity per average ancestor of the Glauber
fit. To illustrate its potential effect on peripheral $R_{\rm AA}$ we use
PHENIX data in 80–92% central AuAu collisions at $\sqrt{s_{\rm NN}}=0.2$ TeV
[10, 11] and CMS data in 70–90% central PbPb collisions at $\sqrt{s_{\rm
NN}}=5.02$ TeV [12]. Above $5.25$ GeV/$c$ the PHENIX data from 2008 and 2012
were averaged using the quadratic sum of statistical and systematic
uncertainties of the original measurements as weights. The PHENIX data are
shown in Fig. 1 up to 10 GeV/$c$ and the CMS data in Fig. 2 up to 30 GeV/$c$.
The error bars represent statistical, while the shaded boxes the systematic
uncertainties. The vertical box around 1 at 0.5 GeV/$c$ denotes the global
normalization uncertainty, which is dominated by the uncertainties on
determining the centrality and $N_{\rm coll}$ (or $T_{\rm AA}$) from Glauber.
As indicated in the figures, constant functions were fit to the PHENIX and CMS
data between $3$–$17$ GeV/$c$ and between $10$–$100$ GeV/$c$, respectively,
yielding a value of $0.80\pm 0.03$ and $0.74\pm 0.02$ with a reduced
$\chi^{2}<1$ using statistical and systematic uncertainties (ignoring the
normalization uncertainty) added in quadrature. Using a linear fit instead of
a constant would in both cases result in a slope consistent with 0. For PbPb
at $5.02$ TeV, this is distinctively different for the $50$–$70$% centrality
class, where $R_{\rm AA}$ exhibits a significant slope of about $0.003$
GeV2/$c^{2}$, indicating that $R_{\rm AA}\sim 1$ is reached at around $125$
GeV/$c$, although parton energy loss should play a stronger role than in the
more peripheral class.

The data are compared to HIJING (v1.383) calculations without jet quenching
and shadowing and a toy model called HG-PYTHIA, which is based on the HIJING
Glauber model for the initial state and PYTHIA [6] as explained below. Besides
the jet quenching and shadowing settings, all other settings in HIJING were
used as set by default, except the minimum $p_{\rm T}$ of hard or semi-hard
scatterings, which was set to $2.3$ (instead of $2.0$) GeV for PbPb collisions
at $5.02$ TeV.

Figure 1: $R_{\rm AA}$ versus $p_{\rm T}$ in 80–92% central AuAu collisions at
$\sqrt{s_{\rm NN}}=0.2$ TeV. The PHENIX data from [10, 11], which were
averaged as explained in the text, are compared to HG-PYTHIA and HIJING
calculations. For details, see text. Figure 2: $R_{\rm AA}$ versus $p_{\rm T}$
in 70–90% central PbPb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV. The CMS
data from [12] are compared to HG-PYTHIA and HIJING calculations. For details,
see text.

HIJING accounts for fluctuations of $N^{\rm hard}_{\rm NN}$ via an NN overlap
function $T_{\rm NN}$ that depends on $b_{\rm NN}$. The probability for an
inelastic NN collision is given by

${\rm d}\sigma_{\rm inel}=2\pi\,b_{\rm NN}\,{\rm d}b_{\rm
NN}\,\left[1-e^{-\left(\sigma_{\rm soft}+\sigma_{\rm hard}\right)\,T_{\rm
NN}(b_{\rm NN})}\right]\,,$ (3)

where $\sigma_{\rm soft}$ is the geometrical soft cross-section of $57$ mb
related to the nucleon size and $\sigma_{\rm hard}$ the energy-dependent pQCD
cross-section for $2\rightarrow 2$ parton scatterings. At $\sqrt{s_{\rm
NN}}=0.2$ and $5.02$ TeV, $\sigma_{\rm hard}=11.7$ and $166.2$ mb,
respectively, while below $0.04$ TeV it is only about $1$ mb. The NN overlap
is approximated by the Eikonal function as $T_{\rm
NN}\propto\xi^{3}\,K_{3}(\xi)$ with $\xi\propto b_{\rm NN}/\sqrt{\sigma_{\rm
soft}}$. If there is at least one hard NN collision, the number of hard
scatterings or multiple parton interactions (MPI) per NN interactions is
distributed as

$P(N^{\rm hard}_{\rm NN})=\frac{{\langle N^{\rm hard}_{\rm NN}\rangle}^{N^{\rm
hard}_{\rm NN}}}{N^{\rm hard}_{\rm NN}!}e^{-\langle N^{\rm hard}_{\rm
NN}\rangle}$ (4)

with the average number of hard scatterings determined by $b_{\rm NN}$ as

$\langle N^{\rm hard}_{\rm NN}\rangle=\sigma_{\rm hard}\,T_{\rm NN}(b_{\rm
NN})\,.$ (5)

For a soft NN collision $N^{\rm hard}_{\rm NN}=0$. The average number of hard
collisions per NN collisions decreases from $\left<N^{\rm hard}_{\rm
NN}\right>=1.77$ at $\sqrt{s_{\rm NN}}=5.02$ TeV to $0.28$ at $0.2$ TeV and
becomes negligible ($<0.05$) below $0.04$ TeV. The total number of hard
scatterings for a AA collision is then obtained by looping over all NN
collisions in the MC Glauber, i.e. $N_{\rm hard}=\sum_{i=1}^{N_{\rm
coll}}\left(N^{\rm hard}_{\rm NN}\right)_{i}$.

In HIJING, $R_{\rm AA}$ at high $p_{\rm T}$ for central collisions is not
unity even in absence of jet quenching and shadowing, possibly due to the
implementation of adhoc energy–momentum conservation for large $N^{\rm
hard}_{\rm NN}$ and dueto the treatment of multiple scattering to model the
“Cronin” effect [13]. This is illustrated in Fig. 4, which shows $R_{\rm AA}$
versus $p_{\rm T}$ in $0$–$5$% central AuAu and PbPb collisions at
$\sqrt{s_{\rm NN}}=0.2$ and $5.02$ TeV, respectively, calculated using HIJING
without jet quenching and without shadowing. At high $p_{\rm T}$, a
suppression by about a factor $2$ is obtained. Such a strong suppression
despite nuclear effects in the model calculation is not consistent with direct
photon measurements [14, 15], which find $R_{\rm AA}\approx 1$ with an
uncertainty of about $\pm 25$% in central collisions.

To overcome these features and to study a baseline corresponding to an
incoherent and unconstrained superposition of NN collisions, we developed HG-
PYTHIA. For a given MC Glauber event calculated in HIJING, we run PYTHIA
(v6.28 with Perugia 2011 tune) to generate pp events with exactly
$\left(N^{\rm hard}_{\rm NN}\right)_{i}$ multiple parton interactions for each
NN collision $i$, where $i$ runs from $0$ to $N_{\rm coll}$. The generated
particles from all PYTHIA events are then combined and treated like a single
AA event in the further analysis.

The model calculations were performed at $\sqrt{s_{\rm NN}}=0.2$ TeV for AuAu
and $5.02$ TeV for PbPb collisions. The respective pp references were obtained
with the same MPI model from HIJING as for AA collisions. The centrality
selection was done similarly to the data, namely by ordering events based on
the charged particle multiplicity in $3<|\eta|<4$ for the low and
$2.5<|\eta|<5$ for the high energy data, respectively. The multiplicity in the
models approximately scales with $N_{\rm hard}$ and hence can not be expected
to be in agreement with the data. However, by ordering the events according to
multiplicity for the centrality selection, this deficiency is not relevant for
the present study. The average number of binary collisions needed for the
calculation of $R_{\rm AA}$ was obtained directly from the model for the
selected range in multiplicity and was found to be about $N_{\rm coll}=4.1$
and $8.6$, respectively, similar to the values reported bythe experiments.

The results of the model calculations are compared to the data in Figs. 1 and
2. As expected due to the selection bias, both calculations are below 1 even
in absence of any initial or final state effect in the model. The models
describe the data at both collision energies, despite a factor of 25
difference in their center-of-mass energies. HG-PYTHIA has the tendency to
overshoot the data, while HIJING to undershoot the data, and appears to
decrease with increasing $p_{\rm T}$. Presumably this originates from the
adhoc energy–momentum conservation and strong multiple scattering implemented
in HIJING, which is purposely neglected in HG-PYTHIA, and hence the difference
between the two models may be taken as a extreme range of model predictions.
The experimental results have large global normalization uncertainties of
about $15$ to $20$%, which however are largely correlated with those of the
calculations, since in all cases the same types of MC Glauber simulation are
used. However, the PHENIX data may be affected by finite trigger inefficiency
in the peripheral bin, whose correction as in the case of dA collisions [16]
would work in the opposite direction than the multiplicity bias [17].

Figure 3: $R_{\rm AA}$ versus centrality for AuAu and PbPb collisions at
$\sqrt{s_{\rm NN}}=0.2$ and $\sqrt{s_{\rm NN}}=5.02$ TeV, respectively,
calculated with HG-PYTHIA, where $R_{\rm AA}\equiv N_{\rm hard}/(N_{\rm
coll}\,\left<N^{\rm hard}_{\rm NN}\right>)$ by construction. For details, see
text.

Finally, we discuss the relative importance of the geometry and selection
biases on $R_{\rm AA}$. In HG-PYTHIA, $R_{\rm AA}$ at high $p_{\rm T}$
approaches $\left.N^{\rm hard}_{\rm NN}\right|_{\rm C}=N_{\rm hard}/N_{\rm
coll}/\left<N^{\rm hard}_{\rm NN}\right>$ from Eq. 2 by construction. The
ratio of $N_{\rm hard}/N_{\rm coll}$ for AuAu and PbPb collisions divided by
the average number of $N_{\rm hard}$ in pp collisions at $\sqrt{s_{\rm
NN}}=0.2$ and $\sqrt{s_{\rm NN}}=5.02$ TeV, respectively, was calculated, and
plotted versus centrality class in $5$% intervals in Fig. 3. The events were
either ordered according to impact parameter probing only the geometry bias,
or according to charged particle multiplicity in the same rapidity ranges as
before, probing both the geometry and the selection biases. It turns out that
the geometry bias sets in at mid-central collisions, reaching up to $10$ to
$20$% for most peripheral collisions. The additional effect of the selection
bias becomes noticeable above the 60% percentile and is significant above the
80% percentile, reaching up to $75\%$ and $40$% for most peripheral
collisions. The onset of the selection bias is steeper for the lower collision
energy, as expected if the fluctuations of the multiplicity and the dynamic
range of $N_{\rm coll}$ are similar. For peripheral classes of $85$% and
beyond, another bias becomes increasingly important, namely the multiplicity
selection in AA collisions effectively cuts into the pp cross section,
suppressing the rate of hard processes relative to minimum bias pp collisions.
While at center-of-mass energies of below 40 GeV [18, 19], the effect from
multiple interactions can be neglected, the latter, so called “jet-veto” [9],
bias is strongly enhanced due to auto-correlations, since centrality
estimators at mid-rapidity were used. Hence, the ratio of yields in central to
peripheral collisions, $R_{\rm CP}$, which is intended to quantify nuclear
modification in absence of pp reference data, would be larger than unity in
absence of nuclear effects.

In conclusion, our model studies which employ particle production including a
nucleon–nucleon impact parameter dependence in AA collisions but no nuclear
effects, demonstrate that $R_{\rm AA}$ at midrapidity in peripheral AA
collisions can be significantly affected by event selection and geometry
biases, similar to those previously discussed in the context of pPb collisions
at the LHC [9]. The apparent suppression induced by the biases for peripheral
$R_{\rm AA}$ results are present at all high energy collisions, and for very
peripheral collisions expected to be stronger for lower energy collisions. In
principle, all hard probes should be similarly affected; examples are the most
peripheral $R_{\rm AA}$ of inclusive J/$\psi$ [20] or jets [21] at LHC
energies. Furthermore, the discussion of the onset of parton energy loss
studied with the RHIC beam energy scan [18, 19] should take the biases into
account when considering the ratio of yields in central to peripheral
collisions, which may lead to an artificial enhancement of $R_{\rm CP}$ in
absence of nuclear effects. Calculations that attempt to address parton energy
loss in peripheral collisions have to account for the effect, or else they
will overestimate the suppression caused by parton energy loss. In particular,
to assess the multiplicity bias, the correlation between soft and hard
particle production has to be modelled and centrality selection has to be
performed in the same way as in the experiment. For example, in the model
calculation [22], the extracted transport coefficient changes by an order of
magnitude from $0.1$ GeV${}^{2}/$fm, ie. similar to cold nuclear matter, for
$R_{\rm AA}\sim 1$ in $80$–$92$% to $1$ GeV2/fm, ie. hot nuclear matter, for
$R_{\rm AA}\sim 0.8$ in $60$–$70$% AuAu collisions at $\sqrt{s_{\rm NN}}=0.2$
TeV. Hence, our studies indicate that expectations of parton energy loss in
light–heavy collision systems [23, 24, 25], which were perhaps guided by
comparing with apparent suppression seen in peripheral $R_{\rm AA}$, should be
revisited.

Our interpretation of the peripheral $R_{\rm AA}$ data also implies that large
values of elliptic flow ($v_{2}$), as measured in peripheral AA collisions
[26], can arise in a system that does not exhibit jet quenching, and hence in
accordance with similar observations in low-energy heavy-ion collisions [27],
or in light–heavy or even light–light collisions (see references in [28]).

To further study the effect of the biases, it would be useful measure
peripheral $R_{\rm AA}$ (or integrated values above $\sim 10$ GeV/$c$) in
finer ($5$% wide) centrality intervals. The trend of the progressively more
peripheral $R_{\rm AA}$ values should clearly demonstrate the effect of the
biases, as one would expect the $R_{\rm AA}$ values to decrease, and not to
rise, as expected if parton energy loss decreases. Further experimental
studies, such as measurements of $R_{\rm AA}$ for reference processes like
photon and intermediate vector boson for which no medium modification is
expected, in peripheral AA collisions, would be needed to experimentally
constrain the biases, but they are not straightforward due to the reduced
parton luminosity, and due to significant contamination from ultra-peripheral
collisions and electro-magnetic processes. Finally, semi-inclusive
measurements as [29, 30] can be used to constrain the energy loss in
peripheral AA and light–heavy collision systems, as they do not rely on the
modeling of the Glauber parameters.

## Acknowledgments

We thank Jamie Nagle and Yen-Jie Lee for fruitful discussions. The work of C.
Loizides is supported by the U.S. Department of Energy, Office of Science,
Office of Nuclear Physics, under contract number DE-AC02-05CH11231.

Figure 4: $R_{\rm AA}$ versus $p_{\rm T}$ in $0$–$5$% central AuAu and PbPb
collisions at $\sqrt{s_{\rm NN}}=0.2$ and $5.02$ TeV, respectively calculated
using HIJING without jet quenching and without shadowing.

## References

  * [1] M. L. Miller, K. Reygers, S. J. Sanders, and P. Steinberg, “Glauber modeling in high energy nuclear collisions,” Ann. Rev. Nucl. Part. Sci. 57 (2007) 205–243, arXiv:nucl-ex/0701025 [nucl-ex].
  * [2] M. Rybczynski, G. Stefanek, W. Broniowski, and P. Bozek, “GLISSANDO 2 : GLauber Initial-State Simulation AND mOre…, ver. 2,” Comput. Phys. Commun. 185 (2014) 1759–1772, arXiv:1310.5475 [nucl-th].
  * [3] B. Alver, M. Baker, C. Loizides, and P. Steinberg, “The PHOBOS Glauber Monte Carlo,” arXiv:0805.4411 [nucl-ex].
  * [4] C. Loizides, J. Nagle, and P. Steinberg, “Improved version of the PHOBOS Glauber Monte Carlo,” SoftwareX 1-2 (2015) 13–18, arXiv:1408.2549 [nucl-ex].
  * [5] X.-N. Wang and M. Gyulassy, “HIJING: A Monte Carlo model for multiple jet production in pp, pA and AA collisions,” Phys. Rev. D44 (1991) 3501–3516.
  * [6] T. Sjostrand, S. Mrenna, and P. Z. Skands, “PYTHIA 6.4 Physics and Manual,” JHEP 05 (2006) 026, arXiv:hep-ph/0603175 [hep-ph].
  * [7] J. Jia, “Influence of the nucleon-nucleon collision geometry on the determination of the nuclear modification factor for nucleon-nucleus and nucleus-nucleus collisions,” Phys. Lett. B681 (2009) 320–325, arXiv:0907.4175 [nucl-th].
  * [8] ALICE Collaboration, B. Abelev et al., “Centrality determination of PbPb collisions at $\sqrt{s_{\rm NN}}=2.76$ TeV with ALICE,” Phys. Rev. C88 no. 4, (2013) 044909, arXiv:1301.4361 [nucl-ex].
  * [9] ALICE Collaboration, J. Adam et al., “Centrality dependence of particle production in pPb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV,” Phys. Rev. C91 no. 6, (2015) 064905, arXiv:1412.6828 [nucl-ex].
  * [10] PHENIX Collaboration, A. Adare et al., “Suppression pattern of neutral pions at high transverse momentum in AuAu collisions at $\sqrt{s_{\rm NN}}=200$ GeV and constraints on medium transport coefficients,” Phys. Rev. Lett. 101 (2008) 232301, arXiv:0801.4020 [nucl-ex].
  * [11] PHENIX Collaboration, A. Adare et al., “Neutral pion production with respect to centrality and reaction plane in AuAu collisions at $\sqrt{s_{\rm NN}}=200$ GeV,” Phys. Rev. C87 no. 3, (2013) 034911, arXiv:1208.2254 [nucl-ex].
  * [12] CMS Collaboration, V. Khachatryan et al., “Charged-particle nuclear modification factors in PbPb and pPb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV,” JHEP 04 (2017) 039, arXiv:1611.01664 [nucl-ex].
  * [13] D. Antreasyan, J. W. Cronin, H. J. Frisch, M. J. Shochet, L. Kluberg, P. A. Piroue, and R. L. Sumner, “Production of hadrons at large transverse momentum in 200, 300 and 400 GeV pp and pn collisions,” Phys. Rev. D19 (1979) 764–778.
  * [14] PHENIX Collaboration, S. Afanasiev et al., “Measurement of direct photons in AuAu collisions at $\sqrt{s_{\rm NN}}=200$ GeV,” Phys. Rev. Lett. 109 (2012) 152302, arXiv:1205.5759 [nucl-ex].
  * [15] CMS Collaboration, S. Chatrchyan et al., “Measurement of isolated photon production in pp and PbPb collisions at $\sqrt{s_{\rm NN}}=2.76$ TeV,” Phys. Lett. B710 (2012) 256–277, arXiv:1201.3093 [nucl-ex].
  * [16] PHENIX Collaboration, A. Adare et al., “Centrality categorization for $R_{\rm p(d)A}$ in high-energy collisions,” Phys. Rev. C90 no. 3, (2014) 034902, arXiv:1310.4793 [nucl-ex].
  * [17] J. Nagle Private communication, 2017.
  * [18] X. Luo, “Exploring the QCD phase structure with Beam Energy Scan in heavy-ion collisions,” Nucl. Phys. A956 (2016) 75–82, arXiv:1512.09215 [nucl-ex].
  * [19] STAR Collaboration, L. Adamczyk et al., “Beam energy dependence of jet-quenching effects in AuAucollisions at $\sqrt{s_{\rm NN}}$ = 7.7, 11.5, 14.5, 19.6, 27, 39, and 62.4 GeV,” arXiv:1707.01988 [nucl-ex].
  * [20] ALICE Collaboration, J. Adam et al., “J/$\psi$ suppression at forward rapidity in PbPb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV,” Phys. Lett. B766 (2017) 212–224, arXiv:1606.08197 [nucl-ex].
  * [21] CMS Collaboration, V. Khachatryan et al., “Measurement of inclusive jet cross sections in pp and PbPb collisions at $\sqrt{s_{\rm NN}}=2.76$ TeV,” Phys. Rev. C96 no. 1, (2017) 015202, arXiv:1609.05383 [nucl-ex].
  * [22] A. Dainese, C. Loizides, and G. Paic, “Leading-particle suppression in high energy nucleus-nucleus collisions,” Eur. Phys. J. C38 (2005) 461–474, arXiv:hep-ph/0406201 [hep-ph].
  * [23] X. Zhang and J. Liao, “Jet quenching and its azimuthal anisotropy in AA and possibly high multiplicity pA and dA collisions,” arXiv:1311.5463 [nucl-th].
  * [24] K. Tywoniuk, “Is there jet quenching in pPb?,” Nucl. Phys. A926 (2014) 85–91.
  * [25] C. Shen, C. Park, J.-F. Paquet, G. S. Denicol, S. Jeon, and C. Gale, “Direct photon production and jet energy-loss in small systems,” Nucl. Phys. A956 (2016) 741–744, arXiv:1601.03070 [hep-ph].
  * [26] CMS Collaboration, S. Chatrchyan et al., “Multiplicity and transverse momentum dependence of two- and four-particle correlations in pPb and PbPb collisions,” Phys. Lett. B724 (2013) 213–240, arXiv:1305.0609 [nucl-ex].
  * [27] STAR Collaboration, L. Adamczyk et al., “Measurement of elliptic flow of light nuclei at $\sqrt{s_{NN}}=$ 200, 62.4, 39, 27, 19.6, 11.5, and 7.7 GeV at the BNL Relativistic Heavy Ion Collider,” Phys. Rev. C94 no. 3, (2016) 034908, arXiv:1601.07052 [nucl-ex].
  * [28] C. Loizides, “Experimental overview on small collision systems at the LHC,” Nucl. Phys. A956 (2016) 200–207, arXiv:1602.09138 [nucl-ex].
  * [29] ALICE Collaboration, J. Adam et al., “Measurement of jet quenching with semi-inclusive hadron-jet distributions in central PbPb collisions at $\sqrt{s_{\rm NN}}=2.76$ TeV,” JHEP 09 (2015) 170, arXiv:1506.03984 [nucl-ex].
  * [30] STAR Collaboration, L. Adamczyk et al., “Measurements of jet quenching with semi-inclusive hadron+jet distributions in AuAu collisions at $\sqrt{s_{\rm NN}}$ = 200 GeV,” arXiv:1702.01108 [nucl-ex].