# Mass Ejection Driven by Sudden Energy Deposition in Stellar Envelopes

Nicholas J. Corso Department of Astronomy, Center for Astrophysics and
Planetary Science, Cornell University, Ithaca, NY 14853, USA  Dong Lai
Department of Astronomy, Center for Astrophysics and Planetary Science,
Cornell University, Ithaca, NY 14853, USA  Tsung-Dao Lee Institute, Shanghai
Jiao-Tong University, Shanghai, China

###### Abstract

A number of stellar astrophysical phenomena, such as tidal novae and planetary
engulfment, involve sudden injection of sub-binding energy in a thin layer
within the star, leading to mass ejection of the stellar envelope. We use a 1D
hydrodynamical model to survey the stellar response and mass loss for various
amounts ($E_{\mathrm{dep}}$) and locations of the energy deposition. We find
that the total mass ejection has a nontrivial dependence on $E_{\mathrm{dep}}$
due to the varying strengths of mass ejection events, which are associated
with density/pressure waves breaking out from the stellar surface. The rapid
occurrence of multiple breakouts may present a unique observational signature
for sudden envelope heating events in stars.

##  1 Introduction

Mass ejection from a star resulting from sudden energy infusion is a common
occurrence during stellar evolution. A well-known example is the thermal
instability associated with shell Helium burning (“He-shell flashes”) in AGB
stars, leading to thermal relaxation oscillations of the stellar envelope and
the associated mass loss (Schwarzschild & Härm, 1965). Another example is
classical novae, luminous eruptions associated with white dwarfs accreting
matter from nondegenerate binary companions (Gallagher & Starrfield, 1978;
Chomiuk et al., 2021). Under appropriate conditions (e.g. a certain range of
accretion rates), the accreted layer undergoes unstable (runaway) nuclear
burning upon reaching a critical mass, and the resulting energy release causes
the accreted envelope to expand, ultimately leading to its ejection.

Enhanced mass loss due to energy deposition in stellar envelopes has also been
discussed in the context of massive stars prior to supernova explosion. Such
elevated mass loss may help explain the luminosity evolution and variability
of Type IIn SNe as a result of the interaction between SN ejecta and
circumstellar material. Quataert & Shiode (2012) (also Shiode & Quataert 2014)
suggest that vigorous core convection in pre-SN stars excite internal gravity
waves that leak out as acoustic waves in the stellar envelope, and the
dissipation of the waves then leads to mass ejection of the outer layer and
pre-SN outbursts. A number of follow-up works have examined various aspects of
this wave-heating scenario (Fuller, 2017; Fuller & Ro, 2018; Leung et al.,
2021a, b; Wu & Fuller, 2021). Energy deposition by a companion during common
envelope evolution may be an alternative mechanism of driving the required
mass loss (Chevalier, 2012). The response of the stellar envelope to a steady
energy injection has been studied by Quataert et al. (2016), Kuriyama &
Shigeyama (2020), and Ko et al. (2022).

Sudden energy release in stellar envelopes can also occur in other contexts.
One example is “tidal novae” in coalescing binary white dwarfs prior to merger
(Fuller & Lai, 2012). In these systems, the time-varying tidal force from the
companion excites internal gravity waves, which then propagate outward and
dissipate in the outer envelope of the white dwarf (e.g. Su et al. 2020). The
heating rate increases rapidly as the binary orbit decays, and at some point
the heating becomes significant enough to ignite hydrogen burning, creating a
thermonuclear runaway that releases even more energy. Another very different
example involves an evolved star engulfing an orbiting giant planet or brown
dwarf (O’Connor et al., 2023). As the planet spirals inward due to gas drag,
it plows through the stellar envelope and releases energy. As the inspiral
speeds up, the largest amount of energy release occurs somewhat suddenly, just
before the planet “dissolves” (e.g., by the tidal force from the stellar
core).

Figure 1: Ejected mass as a function of the dimensionless energy deposition,
$\epsilon=E_{\mathrm{dep}}R_{\ast}/(GM_{\ast}^{2})$ for the deposition radius
$\tilde{r}_{\mathrm{dep}}=0.7$. The solid line with circular points represents
the total ejected mass, while the dashed line with triangular points depicts
only the mass ejected from the first breakout (measured at
$t=2.5t_{\mathrm{dyn}}$). The separation between the two lines indicates the
presence of multiple mass breakout events. Our fiducial resolution has
$N_{\mathrm{cells}}=3000$ (maroon). The dots of different color depicts
results at different resolutions. Also included is a dotted line of spectral
slope $3.2$, representing the scaling of the mass loss in the first breakout
with $\epsilon$.

All of the above examples share the common feature in that a certain amount of
energy is released quickly for a short period of time in a narrow region of
the stellar envelope. The amount of energy release and the location of heat
deposition depend on the specifics of the problem. Some previous works have
explored the response of the star (particularly mass loss) as a function of
the amount and radial location of energy deposition, including Dessart et al.
(2010) and Quataert et al. (2016). However, the former generally concentrates
on the energy deposition close to the core with an amount comparable to the
binding energy of the star. The latter, meanwhile, focuses on very shallow
energy depositions occurring over a relatively long timescale, such that a
quasi-steady state is achieved, leading to a mostly hydrostatic stellar
expansion with a strong wind.

In this study, we study the idealized problem where we assume that a sub-
binding energy deposition occurs instantaneously in a narrow region of the
star. This generally requires that the energy injection timescale be shorter
than the local dynamical time. We determine the total mass loss as a function
of the amount and radial location of the energy deposition. Using 1D
hydrodynamics while ignoring the effects of radiation transport, we hope to
obtain the “base result” on how mass loss scales with the energy deposition,
and lay the groundwork for future explorations where additional physical
effects can be added.

##  2 Problem Setup and Method

Our initial star is a polytrope (satisfying the equation of state
$P\propto\rho^{\gamma}$) with mass $M_{\ast}$ and radius $R_{\ast}$. We
discretize its radial profile using cells of equal steps in radius. These
correspond to non-uniform Lagrangian mass cells that are used in the ensuing
hydrodynamical simulations. Unless stated otherwise, we use 3000 cells. We
choose $\gamma=5/3$ and express mass and radius in units of $M_{\ast}$ and
$R_{\ast}$.

Figure 2: Evolution of the radii of the outermost mass shells as a function of
$t$ for energy injection $\epsilon=0.005,0.01,0.032,0.05$ (left to right, top
to bottom), all with $\tilde{r}_{\mathrm{dep}}=0.7$. The evolution is shown up
to $t=20t_{\mathrm{dyn}}$ to concentrate on the behavior of the second
breakout, even though all simulations were run until at least
$t=45t_{\mathrm{dyn}}$. All cells with external mass $\Delta M\leq
10^{-2}M_{\ast}$ are shown, but the cells at half logarithmic steps in $\Delta
M$ are represented with thicker lines for ease of reference. In each plot, the
thick blue dash-dotted line represents the transition from bound to unbound
material, such that all mass external to the cell is considered ejected.

We simulate the “instantaneous” heat deposition in the star by increasing the
specific internal energy of the gas in a thin layer. We specify this
deposition using two dimensionless parameters. Namely, we deposit a total
energy of $E_{\mathrm{dep}}=\epsilon\left(GM_{\ast}^{2}/R_{\ast}\right)$ in a
layer with the radius range
$r\in[\tilde{r}_{\mathrm{dep}}-\Delta\tilde{r}_{\mathrm{dep}}/2,\tilde{r}_{\mathrm{dep}}+\Delta\tilde{r}_{\mathrm{dep}}/2]R_{\ast}$.
Throughout this study, we choose the width of the deposition layer to be
$\Delta\tilde{r}_{\mathrm{dep}}=0.05$. We consider the dimensionless energy in
the range $\epsilon\in[5\times 10^{-4},10^{-1}]$ and the deposition radius in
the range $\tilde{r}_{\mathrm{dep}}\in[0.6,1)$. Note that by raising the
energy in the deposition layer instantaneously, we require that the energy is
injected on a timescale shorter than the dynamical time,
$t_{\mathrm{dyn}}=\left(R_{\ast}^{3}/GM_{\ast}\right)^{1/2}$.

We run all simulations using a 1D first-order Lagrangian hydrodynamics code
constructed using the scheme described by Bodenheimer et al. (2007) with the
inclusion of self-gravity and artificial viscosity. All simulations last for
at least $45t_{\mathrm{dyn}}$. We choose this duration to provide sufficient
time for the density/pressure waves generated by the initial energy injection
and reflected by the stellar core and surface to affect the energetics of the
stellar envelope. In some cases, it is necessary to run for a longer duration
(up to $\sim 150t_{\mathrm{dyn}}$) in order for the flow to reach a quasi-
steady state and the mass ejection to be reliably determined.

To determine the mass ejection following the energy injection, we consider the
Bernoulli parameter for each mass shell,

$B=\frac{1}{2}v^{2}+\Phi_{\mathrm{g}}+\left(\frac{\gamma}{\gamma-1}\right)\frac{P}{\rho},$
(1)

where $\Phi_{\mathrm{g}}$ is the gravitational potential and $v$ is the radial
velocity. Note that $B$ is constant for barotropic flows in steady state; this
is generically not satisfied for our system after the energy deposition.
However, if the simulation is run for a sufficiently long time and the
outflowing mass has reached a steady state, the ejected material approximately
settles into a barotropic structure with no shock dissipation. As such, we
measure $B$ across the profile for every simulation timestep, and after the
system has reached its quasi-steady state, we consider a surface layer with
$B>0$ to be unbound and therefore ejected.

Figure 3: Same as Figure 2 for $\epsilon=3.2\times 10^{-2}$ and
$\tilde{r}_{\mathrm{dep}}=0.8$. In this case, waves close to the unbound layer
prevent a quasi-steady state from forming at $t\lesssim 45t_{\mathrm{dyn}}$,
and an accurate determination of the ejected mass requires long-term
integration.

##  3 Results

###  3.1 Fiducial Energy Deposition Radius

Our fiducial runs assume the energy deposition radius
$\tilde{r}_{\mathrm{dep}}=0.7$. Figure 1 shows the mass ejection,
$M_{\mathrm{ej}}$, as a function of
$\epsilon=E_{\mathrm{dep}}/\left(GM_{\ast}^{2}/R_{\ast}\right)$ for our
fiducial runs (solid line). While $M_{\mathrm{ej}}$ generally increases with
$\epsilon$, a nontrivial “wavy” behavior emerges in the considered domain,
revealing critical $\epsilon$ regions where the trend falls off, only to
quickly revive. Included in this figure is another curve labeled “First
Breakout,” which shows $M_{\mathrm{ej}}$ measured at $t=2.3t_{\mathrm{dyn}}$.
In our simulations, distinct mass ejection events (hereby termed “breakouts”)
occur, and $t=2.3t_{\mathrm{dyn}}$ generally divides the first and second
breakouts, allowing us to observe the former separately. Despite the
irregularity in the total $M_{\mathrm{ej}}$ vs. $\epsilon$ relation, the mass
loss from the first breakout, $M_{\mathrm{ej}}^{\mathrm{(first)}}$, smoothly
follows the power law,
$M_{\mathrm{ej}}^{\mathrm{(first)}}\propto\epsilon^{3.2}$. Figure 1 also shows
the results from the simulations with different resolutions (characterized by
the number of radial cells, $N_{\mathrm{cells}}$), indicating the level of
convergence of our results.

We contextualize the nontrivial trend of the $M_{\mathrm{ej}}$ vs. $\epsilon$
relation using Figure 2, which depicts the time evolution of the radii of the
outermost mass shells until $t=20t_{\mathrm{dyn}}$ for a variety of test
values of $\epsilon$. By tracing the cells of unbound material, we see that in
each case the initial wave breakout ejects some mass. In most cases, a second
wave breakout ejects more mass. However, in the other, often more energetic,
ones, waves continually disturb the region surrounding the transition from
bound to unbound material, preventing it from settling into a quasi-steady
state. We ran these simulations, such as the one depicted in Figure 3, for
much longer durations, upwards of $t=120t_{\mathrm{dyn}}$, until we achieved a
quasi-steady state. In this example, we identify more than two breakout
events. We observe this behavior for $\epsilon$ close to the critical regions,
where the energy deposition is distributed more evenly among the rebounded
material.

Figure 4: Profiles of radius, radial velocity, density, and specific internal
energy (top to bottom) as a function of the external mass at time steps
$t=0.2,2.2,3.6,8.0,45.0t_{\mathrm{dyn}}$ (yellow, orange, red, purple, black)
for our fiducial run with $\epsilon=10^{-2}$ and
$\tilde{r}_{\mathrm{dep}}=0.7$. The specific time steps are chosen to
highlight the key stages in the time evolution of the system. In order, they
are: (i) the initial generation of the ingoing and outgoing density waves;
(ii) the rebound of the initially ingoing wave after the initially outgoing
wave has ejected material; (iii) the generation of the shock of the second
breakout as material that failed to eject during the first breakout crashes
back down onto the stellar surface; (iv) the escape of the second breakout and
beginning of pulsations in the remnant star; (v) the weakening of remnant
oscillations in the underdense region. Bound material is represented by solid
lines, while unbound material is represented by dotted lines.  Figure 5:
Ejected mass as a function of $\epsilon$ for
$\tilde{r}_{\mathrm{dep}}=0.6,0.7,0.8$ (bright to dark in increasing
$\tilde{r}_{\mathrm{dep}}$). As in Figure 1, the total mass ejected (solid
line, circular points) is distinguished from the mass ejected by the first
breakout alone (dashed line, triangular points), and the first breakout is
described by a power law of spectral slope 3.2.

By comparing different values of $\epsilon$, the trend introduced in Figure 1
comes into focus. Most importantly, we see that at the critical values of
$\epsilon$ with less ejected material, rebounding waves within the star do not
contain enough energy to push through the unbound material and create a second
breakout. As $\epsilon$ approaches these critical regions, the delay between
the first and second breakouts increases. By contrast, at the $\epsilon$
values where the ejected mass locally peaks, the delay between the first and
second breakouts is minimized. In other words, the timing of this delay
appears to be tightly related to how much mass is ejected in the second
breakout. We also note that above the critical region located at
$\epsilon\simeq 0.03$, more than one wave rebound is required to produce the
second breakout.

Figure 4 depicts the radius, radial velocity, density, and specific internal
energy profiles as a function of $\Delta M$ (the amount of mass external to a
shell) for the simulation initiated using our fiducial energy and deposition
radius. We choose five key moments in the simulation to paint a more detailed
picture of the stellar response to the heat deposition. We see that the energy
deposition produces ingoing and outgoing density/pressure waves. The outgoing
wave directly ejects the outer layers, causing the first breakout, while the
ingoing wave reflects off the stellar core before proceeding to the surface to
eject more material, causing the second breakout. Although much of the energy
is carried off by unbound material from these two breakouts, matter that fails
to unbind falls back to produce more ingoing waves, and these waves will
approach the surface again after reflecting off the core. However, by this
point, they have dissipated enough energy that they are no longer able to
unbind more material from the star. Instead, they continuously reflect,
causing the star to pulsate on a dynamical timescale. Remaining at the site of
the original heat deposition is a remnant hot bubble, which had rapidly
expanded at the start of the simulation. Additionally, a hot extended
atmosphere lingers, composed of material just beneath the ejected layer.

Figure 6: Same as Figure 5, but for $\tilde{r}_{\mathrm{dep}}=0.9,0.93,0.95$.
The shaded boxes correspond to the region of $\Delta M$ at the initial heat
deposition, with the $\epsilon$ boundaries constrained to highlight its
correspondence with the ‘high-energy’ regime. Upon reaching this $\Delta M$
region, the slope of the curve flattens, with the black dotted line of
spectral slope $1$ representing the new scaling.

###  3.2 $\tilde{r}_{\mathrm{dep}}$-Dependence

In Figure 5 we repeat the analysis performed in Figure 1, now for the energy
deposition radius $\tilde{r}_{\mathrm{dep}}=0.6,0.7,0.8$. The three curves
follow a similar trend, exhibiting the nontrivial “wavy” behavior shown in
Figure 1. However, for $\tilde{r}_{\mathrm{dep}}=0.8$, we note that the wavy
behavior is weaker compared to the other cases, and the second breakout always
occurs. It should also be noted that in some cases, often with higher
$\epsilon$ values, more than two breakouts may occur. As in Figure 1, the mass
ejected from the first breakout exhibits a rather consistent power law
dependence on $\epsilon$, i.e. $M_{\mathrm{ej}}\propto\epsilon^{3.2}$. The
steepness of this slope arises both from the direct energy increase as well as
the increasing capability of the injected energy to reach directly to the
surface to unbind material rather than distribute to the rest of the star.

Figure 6 again repeats this analysis, now for larger
$\tilde{r}_{\mathrm{dep}}$ (energy deposition closer to the stellar surface).
In doing so, we reveal a new, ‘high-energy’ regime. Specifically, if the
energy deposited is sufficiently large, the first breakout may eject most of
the mass that had contained the original heat deposition layer. As a result,
the remnant bubble at this layer may expand supersonically, forming a shock.
If the energy injection is larger still, this shock, or even part of the
heated remnant itself, may efficiently escape, causing the amount of ejected
mass to noticeably increase. Afterward, the internal waves struggle to eject
more material since the envelope has expanded so much, suppressing the second
breakout almost entirely. The strong power law dependence of the first
breakout mass ejection on $\epsilon$ transforms to a modest
$M_{\mathrm{ej}}^{\mathrm{(first)}}\propto\epsilon$, and this slope decreases
still at even higher values of $\epsilon$. Thus, in this ‘high-energy’ regime,
the first breakout ‘saturates,’ such that increasing the deposited energy
mostly energizes material that is already able to unbind at lower energies.
Figure 6 highlights this by showing the mass external to the inner radius of
the heat deposition layer. In the ‘high-energy’ regime, the curves asymptote
to this mass, indicating that the shell fully escapes without any interior
material.

##  4 Discussion and Conclusions

We have used 1D hydrodynamical simulations to model the evolution of a star
following a sudden energy deposition in the stellar interior, to determine the
stellar mass loss as a function of the amount ($E_{\mathrm{dep}}$) and
location ($\tilde{r}_{\mathrm{dep}}$) of the energy deposition. We find that
the ejected mass generally has a nontrivial dependence on $E_{\mathrm{dep}}$
(see Figs. 5 and 6). The mass ejection is associated with the breakouts of
density/pressure waves that are generated by the energy deposition and the
ensuing wave reflections from the stellar core and surface. While the mass
ejection from the first wave breakout has a simple power-law dependence on the
energy deposition ($M_{\mathrm{ej}}^{(\mathrm{first})}\propto
E_{\mathrm{dep}}^{3.2}$), the total mass loss generally exhibits a wavy
behavior due to the breakout of a second wave reflected off the stellar core.
The strength of the second breakout depends on the time delay between the two
breakouts, and in some cases there may be no second breakout. For sufficiently
large energy deposition close to the stellar surface, the envelope exterior to
the deposition radius is efficiently ejected, and the total mass loss has a
much weaker dependence on $E_{\mathrm{dep}}$.

The nontrivial dependence of $M_{\mathrm{ej}}$ on $E_{\mathrm{dep}}$ implies
that in astrophysical situations involving sudden energy deposition (such as
tidal nova and planet engulfment), one should exercise caution when trying to
infer the total energetics of the event from observations. The multiple mass
breakouts may be a unique feature of the event. Of course, since our study
does not include any treatment of radiation, further works are needed to
clarify the observational signatures of such energy deposition events.

Beyond the inclusion of radiation physics, there are a number of extensions of
our work. First, we have adopted a simplified stellar model (i.e. polytrope)
in this paper, and it will be useful to consider more realistic stellar models
in future works. Second, we have assumed sudden energy injection in this
paper. It will be important to consider spreading the energy injection over a
finite duration of time, which may result in a more gradual mass ejection (see
Quataert et al. 2016). Extending the duration may blur the distinction between
the multiple mass breakout events. Finally, we have focused on 1D flows in
this paper. Expanding the simulations to 2D or 3D would help to address
anisotropies that may emerge in the system, particularly in the situations
where the energy deposition is not spherically symmetric.

This research was partially supported by the NSF grants AST-2107796 and
DGE-2139899.

## References

  * Bodenheimer et al. (2007) Bodenheimer, P., Laughlin, G. P., Różyczka, M., & Yorke, H. W. 2007, Numerical Methods in Astrophysics: An Introduction (Boca Raton, Fl: CRC Press), doi: 10.1201/9781420011869
  * Chevalier (2012) Chevalier, R. A. 2012, The Astrophysical Journal Letters, 752, L2, doi: 10.1088/2041-8205/752/1/L2
  * Chomiuk et al. (2021) Chomiuk, L., Metzger, B. D., & Shen, K. J. 2021, ARA&A, 59, 391, doi: 10.1146/annurev-astro-112420-114502
  * Dessart et al. (2010) Dessart, L., Livne, E., & Waldman, R. 2010, MNRAS, 405, 2113, doi: 10.1111/j.1365-2966.2010.16626.x
  * Fuller (2017) Fuller, J. 2017, Monthly Notices of the Royal Astronomical Society, 470, 1642, doi: 10.1093/mnras/stx1314
  * Fuller & Lai (2012) Fuller, J., & Lai, D. 2012, ApJ, 756, L17, doi: 10.1088/2041-8205/756/1/L17
  * Fuller & Ro (2018) Fuller, J., & Ro, S. 2018, Monthly Notices of the Royal Astronomical Society, 476, 1853, doi: 10.1093/mnras/sty369
  * Gallagher & Starrfield (1978) Gallagher, J. S., & Starrfield, S. 1978, Annual Review of Astronomy and Astrophysics, 16, 171, doi: 10.1146/annurev.aa.16.090178.001131
  * Ko et al. (2022) Ko, T., Tsuna, D., Takei, Y., & Shigeyama, T. 2022, The Astrophysical Journal, 930, 168, doi: 10.3847/1538-4357/ac67e1
  * Kuriyama & Shigeyama (2020) Kuriyama, N., & Shigeyama, T. 2020, A&A, 635, A127, doi: 10.1051/0004-6361/201937226
  * Leung et al. (2021a) Leung, S.-C., Fuller, J., & Nomoto, K. 2021a, ApJ, 915, 80, doi: 10.3847/1538-4357/abfcbe
  * Leung et al. (2021b) Leung, S.-C., Wu, S., & Fuller, J. 2021b, ApJ, 923, 41, doi: 10.3847/1538-4357/ac2c63
  * O’Connor et al. (2023) O’Connor, C. E., Bildsten, L., Cantiello, M., & Lai, D. 2023, ApJ, 950, 128, doi: 10.3847/1538-4357/acd2d4
  * Quataert et al. (2016) Quataert, E., Fernández, R., Kasen, D., Klion, H., & Paxton, B. 2016, MNRAS, 458, 1214, doi: 10.1093/mnras/stw365
  * Quataert & Shiode (2012) Quataert, E., & Shiode, J. 2012, MNRAS, 423, L92, doi: 10.1111/j.1745-3933.2012.01264.x
  * Schwarzschild & Härm (1965) Schwarzschild, M., & Härm, R. 1965, ApJ, 142, 855, doi: 10.1086/148358
  * Shiode & Quataert (2014) Shiode, J. H., & Quataert, E. 2014, The Astrophysical Journal, 780, 96, doi: 10.1088/0004-637X/780/1/96
  * Su et al. (2020) Su, Y., Lecoanet, D., & Lai, D. 2020, MNRAS, 495, 1239, doi: 10.1093/mnras/staa1306
  * Wu & Fuller (2021) Wu, S., & Fuller, J. 2021, The Astrophysical Journal, 906, 3, doi: 10.3847/1538-4357/abc87c


# Benchmarking Information Scrambling

Joseph Harris  Theoretical Division, Los Alamos National Laboratory, Los
Alamos, New Mexico 87545, USA  Department of Applied Mathematics and
Theoretical Physics, University of Cambridge, Cambridge, CB3 0WA, UK  School
of Physics and Astronomy, University of Edinburgh, Edinburgh, EH9 3FD, UK  Bin
Yan  Theoretical Division, Los Alamos National Laboratory, Los Alamos, New
Mexico 87545, USA  Center for Nonlinear Studies, Los Alamos National
Laboratory, Los Alamos, New Mexico 87545, USA  Nikolai A. Sinitsyn
Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico
87545, USA

###### Abstract

Information scrambling refers to the rapid spreading of initially localized
information over an entire system, via the generation of global entanglement.
This effect is usually detected by measuring a temporal decay of the out-of-
time order correlators. However, in experiments, decays of these correlators
suffer from fake positive signals from various sources, e.g., decoherence due
to inevitable couplings to the environment, or errors that cause mismatches
between the purported forward and backward evolutions. In this work, we
provide a simple and robust approach to single out the effect of genuine
scrambling. This allows us to benchmark the scrambling process by quantifying
the degree of the scrambling from the noisy backgrounds.

In complex many-body systems, initially localized information quickly spreads
throughout the entire system — a process known as information scrambling.
Though information is ultimately conserved, it gets encoded into global
entanglement among many degrees of freedom, and hence becomes inaccessible by
local measurement. Information scrambling was originally studied in the
context of black hole physics Hayden and Preskill (2007); Kitaev (2015);
Maldacena  _et al._ (2016), and has since emerged as a field with a wide-
ranging impact across different areas in physics, e.g., quantum chaos in many-
body systems Rozenbaum  _et al._ (2017); Roberts and Yoshida (2017); Cotler
_et al._ (2017); Lin and Motrunich (2018); García-Mata  _et al._ (2018); Yan
_et al._ (2020); Yan and Sinitsyn (2020), phase transition Sahu  _et al._
(2019); Choi  _et al._ (2020), and quantum machine learning Holmes  _et al._
(2021). Considerable effort has also been made in probing this effect in
various experimental systems Landsman  _et al._ (2019); Li  _et al._ (2017);
Gärttner  _et al._ (2017).

Information scrambling is usually measured by the temporal decay of the out-
of-time order correlators (OTOCs) Kitaev (2015); Larkin and Ovchinnikov
(1969); Swingle (2018), defined as

$\langle W^{\dagger}(t)V^{\dagger}W(t)V\rangle,$ (1)

where the average is taken over a given quantum state. $W$ and $V$ are local
operators, usually considered to act on distinct subsystems. $W(t)$ is the
Heisenberg evolution of $W$, which becomes a global operator as scrambling
proceeds, causing decay of the OTOCs.

However, it is difficult to distinguish between scrambling and decoherence:
the latter causes leakage of the system information to the environment and, in
general, induces decay of the OTOC as well. Protocols that measure the OTOC
often involve forward and backward evolution of the system, which in practice
do not exactly match each other due to operational errors. This can also cause
decay of the OTOC signals. Understanding the behavior of scrambling in
presence of decoherence and operational errors is an active and ongoing
research in the field Swingle and Yunger Halpern (2018); Knap (2018); Syzranov
_et al._ (2018); Zhang  _et al._ (2019); González Alonso  _et al._ (2019);
Yoshida and Yao (2019); Joshi  _et al._ (2020); Touil and Deffner (2021);
Domínguez and Álvarez (2021); Zanardi and Anand (2021). In this line of
effect, the first experiment that provided a positive verification of
scrambling was performed with trapped ions Landsman  _et al._ (2019), using a
teleportation-based protocol Yoshida and Yao (2019). This approach uses two
copies of the system and an entangled input state (Bell state) between the
copies, requiring sophisticated engineering of the system and therefore
hindering its practical applications. More recently, the field has witnessed
an increasing number of studies of information scrambling in various
experimental platforms, e.g., superconductors Zhao  _et al._ (2021), trapped-
ions Joshi  _et al._ (2020), and cloud-based quantum computers Braumüller  _et
al._ (2021); Mi  _et al._ (2021); Geller  _et al._ (2021). However, it remains
a challenge to design a simple and robust protocol for benchmarking the true
signals of scrambling from a noisy background. Thus a solution to this problem
is in demand and is beneficial for the development of the field, especially
for experimental studies.

In this work, we propose a solution to this task. Our approach is based on a
novel quantum butterfly effect Yan and Sinitsyn (2020); Sinitsyn and Yan
(2020) that can unambiguously distinguish scrambling and non-scrambling
dynamics. Due to the global entanglement generated by scrambling, information
becomes more robust against local perturbations, and hence, when scrambled and
disturbed by local perturbations, can still be partially recovered through a
reversed unscrambling process. In contrast, decoherence and errors only temper
this effect and weaken the signal produced by true scrambling. In the
following, we develop the underlying theory, introduce the benchmarking
protocol, and apply it to a model problem of fast scrambling in the presence
of decoherence.

Figure 1: In the Lorenz picture of the butterfly effect, one compares two
trajectories evolved under the same Hamiltonian, but with slightly different
initial conditions. In the Bradbury picture, perturbation is applied in the
past. These two pictures are equivalent in the classical world but have subtle
differences in quantum dynamics.

In chaotic classical dynamics, small perturbations in the initial conditions
can trigger dramatic changes in the time evolution. This effect is well-known
as the Lorenz butterfly effect. A few years before Lorenz, Bradbury in his
famous story Bradbury (1952) introduced a variant of the butterfly effect. In
this picture, one starts with an initial state in the present, evolves
backward in time, applies a small perturbation, and then travels forward back
to the present.

These two pictures offer equivalent descriptions of the same classical
butterfly effect (Fig. 1). However, they exhibit subtle differences in quantum
dynamics Yan and Sinitsyn (2020); Cao and Scaffidi (2021): the overlap between
the two trajectories (wavefunctions) in the Lorenz picture remains a constant
during time evolution, due to the unitarity of (isolated) quantum dynamics. On
the other hand, in the Bradbury’s picture, the overlap between the two states
at time $t=0$ — the initial input state and the final output state after the
backward and forward evolution loop — does decay as a function of the
evolution time. Moreover, asymptotically, the output state always contains
partial information of the initial one, with the amount of information
determined by the type of perturbation. This is in sharp contrast with
classical chaotic dynamics, which on average smear all the initial information
over the entire accessible phase space. For these reasons, we propose to call
this phenomenon the quantum anti-butterfly effect.

More precisely, the perturbation in the ‘past’ can be described by a general
quantum channel $\Lambda$:

$\rho\rightarrow\Lambda(\rho)=\sum_{k}M^{\dagger}_{k}\rho
M_{k},\quad\sum_{k}M_{k}M_{k}^{\dagger}=\mathbb{I},$ (2)

where $\mathbb{I}$ is the identity operator, and $M_{k}$ are the Kraus
operators. The initial state $\rho$, after the Bradbury’s process
aforementioned, becomes

$\rho(t)=U_{t}\Lambda(U^{\dagger}_{t}\rho U_{t})U^{\dagger}_{t},$ (3)

where $U_{t}$ is the evolution operator for a time $t$. This process is also
recognized as the quantum twirling channel of the perturbation $\Lambda$
Dankert  _et al._ (2009).

The fidelity of the output state of the twirling channel, with respect to the
initial state, is

$F(t)\equiv{\rm Tr}\left[\rho(t)\cdot\rho\right]=\sum_{k}{\rm
Tr}\left[M^{\dagger}_{k}(t)\rho M_{k}(t)\rho\right].$ (4)

This quantity is the central object for the protocol developed in the
following section.

After a long time evolution with a chaotic Hamiltonian, when the evolution
operator becomes sufficiently random, the asymptotic state can on average be
described as

$\rho_{\rm as}=p\rho+(1-p)\frac{\mathbb{I}}{d},$ (5)

where $d$ the dimension of the total Hilbert space and the probability $p$ is
determined by the error channel, namely,

$p=\frac{\sum_{k}{\rm Tr}M_{k}{\rm Tr}M^{\dagger}_{k}-1}{d^{2}-1}$ (6)

The corresponding asymptotic value of the fidelity (4) has the simple
expression

$F_{\rm as}=p{\rm Tr}\left[\rho^{2}\right]+(1-p)/d.$ (7)

At this point, we would like to emphasize several remarks:

i) The asymptotic state (5) is obtained by averaging over an assemble of
random unitaries with respect to Haar measure 111Technically, the same result
can be obtained for an ensemble of unitaries which is as random as a unitary
2-design, since the average is performed for an expression in terms of the
second moment of the unitary.. However, the fluctuation from the mean is
exponentially small in the size of the total system. Consequently, for a
single unitary randomly drawn from such an ensemble, deviations from this
averaged behavior are exponentially suppressed. Therefore, we can treat (5) as
the asymptotic state for a typical complex evolution.

ii) The fidelity (4) is identified as a sum of special types of OTOCs between
the state $\rho$ and the Kraus operators. Hence, scrambling causes decay of
the fidelity in the same manner as to the OTOCs. This particular type of OTOC
has been used Gärttner  _et al._ (2017) and demonstrated to have various
benefits Lewis-Swan  _et al._ (2019). As will be shown in the following, this
quantity can also acquire a large universal asymptotic value (7), which will
response further to decoherence and errors. This unique feather lies in the
core of our protocol for singling out information scrambling.

Inspired by the above theory of the anti-butterfly effect, we propose the
following protocol to detect and benchmark information scrambling. In this
protocol, the total system is sent through a Bradbury process, and only the
fidelity of a subsystem is measured:

  1. 1.

Initialize the total system such that a small subsystem is prepared in state
$\rho_{S}$.

  2. 2.

Evolve the system forward for a time $t$; apply a perturbation, and then
evolve the system backward for the same time $t$.

  3. 3.

Measure the same subsystem, and evaluate its fidelity with respect to the
initial state $\rho_{S}$.

We choose to measure the fidelity of a subsystem, since evaluating that of the
total system is generally infeasible for large systems. One can verify that
the asymptotic value of the fidelity (7) also applies to any subsystems.

For example, for a (large) spin-$1/2$ system with scrambling dynamics, we
prepare a single target spin in any pure state and choose the perturbation
channel as a projective measurement on any single spin. The fidelity of the
target spin, computed with (7) and (6), saturates to $\sim 0.75$. On the other
hand, in presence of decoherence and errors, the asymptotic state of the
target spin will become random, hence resulting in $\sim 0.5$ fidelity. If we
chooses to measure the fidelity of the total system, the corresponding values
of the asymptotic fidelity would be $\sim 0.5$ and zero for ideal scrambling
and decoherence respectively.

We also note the similarity between our approach and randomized benchmarking
(RB) protocol Emerson  _et al._ (2005); Magesan  _et al._ (2011) for quantum
computers. Indeed, our approach inherits the main benefits of RB, i.e., it is
scalable and is independent with errors in state initialization and readout.

In practice, subsystem state initialization and measurement may not be
directly accessible. For these systems, the corresponding steps in the
protocol can be adapted as

  1. 1’.

Measure a local observable with a discrete set of outcomes.

  2. 3’.

Measure the same observable, and evaluate the coincide rate, i.e., the
probability that the two measurements produce the same outcome.

It can be verified that the coincidence rate is identical to the fidelity of
the same subsystem on which the measurements are performed.

Figure 2: Measurements of the OTOC. (a) Quantum circuit for the interference
protocol (left). Structure of a single evolution step of the fast scrambling
model. (b) OTOC for the cases of scrambling and decoherence (Table 1), which
produces similar signals that can hardly be distinguished.

As an illustration, we study a fast scrambling model recently proposed in Ref.
Belyansky  _et al._ (2020). We perform simulations of this model under
decoherence and compare our protocol with a standard OTOC measurement.

The evolution unitary in the fast scrambling model consists of repetitive
layers of circuit evolution. Each layer (Fig. 2a, right) is composed of Haar
random single-qubit unitaries, immediately followed by a global entangling
gate, i.e., which for $n$ qubits is given by

$\exp{\left(-i\frac{g}{2\sqrt{n}}\sum_{i<j}Z_{i}Z_{j}\right)},$ (8)

where $g$ is a constant parameter that controls the scrambling strength. This
building block can be viewed as a trotterization of the evolution generated by
a spin Hamiltonian with strong random local fields (generating single-gate
random rotations) and all-to-all two-body couplings (creating the global
entangling gate).

To simulate the effect of a noisy environment, we introduce errors in each
layer of the single-qubit gates — after each single-qubit Haar random gate, a
Pauli $X$-gate and a Pauli $Z$-gate are applied independently with probability
$p$. Note that this error channel of the system qubits can be extended to a
unitary evolution in an enlarged Hilbert space including ancillary qubits.
Hence, the error model describes a decoherence process as well. We now have
two parameters, $g$ and $p$, that control the strength of scrambling and
decoherence, respectively. Table 1 lists the parameters we considered, namely,
two cases for ideal scrambling (S1 and S2), two cases for strong decoherence
(D1 and D2), and the cases (I) where scrambling and decoherence become
comparable and compete with each other.

Table 1: Cases compared for the fast scrambling model. The rate of scrambling increases with the value of $g$. The strength of the decoherence increase with the value of $p$. cases | S1 | S2 | D1 | D2 | I
---|---|---|---|---|---
$g$ | 1 | 2 | 0.5 | 0.5 | 1
$p$ | 0 | 0 | 0.025 | 0.1 | 0.001

The OTOC measurement is achieved using the interference protocol developed in
Ref. Swingle  _et al._ (2016), as shown in the quantum circuit diagram in Fig.
2a. The intermediate $X$-gate is placed on the $i$-th qubit $q_{i}$. The final
measurement of the average value $\langle\sigma^{0}_{x}\rangle$ (or
$\langle\sigma^{0}_{y}\rangle$) on the ancillary qubit $q_{0}$ then gives the
real (or imaginary) part of the OTOC
$\langle\sigma^{0}_{x}(t)\sigma^{i}_{x}\sigma^{0}_{x}(t)\sigma^{i}_{x}\rangle$,
where $t$, in analog to time, is the number of layers in the forward (and
hence the backward) evolution.

Figure 3: Measurement of the recovery signal (fidelity) of the benchmark
protocol. (a) Schematic of the quantum circuit. (b) Fidelity decay for the
cases of scrambling and decoherence (Table 1). (c, d) Fidelity as a function
of the number of layers ($t_{2}$) in the backward evolution, for fixed numbers
of layers $t_{1}$ in the forward evolution. (e, g) Density plot for the
fidelity.

Figure 2b compares the evolution of the OTOC for ideal scrambling (without
noise) and weak scrambling with strong decoherence, which exhibit similar
decay curves. Hence, these two situations are practically indistinguishable
from the OTOC measurements.

We now examine the performance of our benchmark protocol for the fast
scrambling model under the same conditions. The circuit diagram for this
protocol is shown in Fig. 3a. Here, all the qubits are prepared in the
computational $|0\rangle$. The recovery signal is obtained by measuring the
fidelity of the first qubit $q_{1}$ with respect to its initial state. We also
specify the intermediate perturbation between the forward and backward
evolution as a projective measurement on a single qubit (other than the
qubit-1) along the $Z$-axis. Note that since the evolution unitary contains
Haar-random single-qubit gates in each layer, the measurement in a fixed
direction is equivalent to random projective measurements. As discussed in the
foregoing section, for this particular perturbation channel, the expected
value of fidelity (4) is $0.75$ for ideal scrambling unitaries. On the other
hand, in the case of strong decoherence, when the qubits eventually lose their
coherence information, the fidelity would be trivially $0.5$. Hence, the cases
of scrambling with and without decoherence exhibit distinct asymptotic values
of fidelity. This is clearly demonstrated in Fig. 3b.

To vitalize the emergence of the recovery signal produced by scrambling, we
also performed simulations with different numbers of layers in the forward and
backward evolution. That is, the fidelity is measured after $t_{1}$ and
$t_{2}$ layers of forward and backward evolution respectively, and if
$t_{2}>t_{1}$, an additional $t_{2}-t_{1}$ layers in the backward evolution
are chosen as independently random. Figures 3(c, d) depict the fidelity as a
function of $t_{2}$ for fixed values of $t_{1}$. The density plots for the
fidelity scanned through various $t_{1}$ and $t_{2}$ are shown in Fig. 3(e,
f). The recovery signal emerges in a finite window around the peak
$t_{1}=t_{2}$. The width of the peak reflects the time scale for local
dissipation. For the current fast scrambling model, due to the random single-
qubit rotation in every single layer, the recovery signal disappears as soon
as $t_{2}$ is one layer away from $t_{1}$.

The above simulations demonstrate that our protocol can unambiguously
distinguish between scrambling and decoherence. In general, both of these two
competing factors contribute to the decay of the fidelity, but they make the
fidelity saturate to different values. This gives rise to a two-stage decay:
In the early scrambling stage, fidelity decay is influenced by both scrambling
and decoherence, until the information is fully scrambled and the fidelity
reaches the saturation value (7). In the latter decoherence stage, fidelity
further decays to a smaller value, at a rate determined purely by decoherence.
The appearance of the two-stage decay indicates the presence of scrambling.

Figure 4: Fidelity decay for case I (Table 1) with both scrambling and
decoherence. Squares are numerical data. The red solid curve is the best fit
to the ansatz (10). To visualize the two stages of the fidelity decay, we also
plotted (black dashed) the long time pure exponential part of the ansatz (10),
$a_{2}\exp{(-\lambda_{d}t)}+F^{d}_{\rm as}$, which clearly departs from the
early decay. Insets: $F(t)-F^{d}_{\rm as}$ in semi-log scale, which show a
single exponential decay at long times, and a sum of two exponential decays at
early times. At the very beginning, fidelity decay takes a quadratic form Yan
and Zurek (2021).

Suppose $\lambda_{s}$ and $\lambda_{d}$ are the exponential decay rates
corresponding to scrambling and decoherence, respectively. In the extreme
cases of strong scrambling and strong decoherence, the decay of the fidelity
can be described as

$F(t)=\begin{cases}\left(1-F^{s}_{\rm as}\right)e^{-\lambda_{s}t}+F^{s}_{\rm
as},&\lambda_{s}\gg\lambda_{d}\\\ \left(1-F^{d}_{\rm
as}\right)e^{-\lambda_{d}t}+F^{d}_{\rm
as},&\lambda_{s}\ll\lambda_{d}\end{cases}$ (9)

where $F^{s}_{\rm as}$ and $F^{d}_{\rm as}$ are the asymptotic values of the
fidelity induced by scrambling and decoherence, respectively. In the
intermediate regime, such that the scrambling and decoherence are comparable,
we propose an ansatz of the fidelity decay:

$F(t)=\left(a_{1}e^{-\lambda_{s}t}+a_{2}\right)e^{-\lambda_{d}t}+F^{d}_{\rm
as}.$ (10)

With this, we can fit the observed fidelity at long times (when the system
gets sufficiently scrambled and the first exponential term vanishes) to a pure
exponential function and extract the decoherence rate $\lambda_{d}$. Then the
scrambling rate $\lambda_{s}$ can be extracted by fitting the early time
fidelity to the ansatz (10). We applied this procedure to case I (Table 1),
where both scrambling and decoherence contribute substantially to the fidelity
decay. Figure 4 shows simulations of the fidelity decay, which is described
accurately by ansatz (10). This allows us to separate the scrambling rate
$\lambda_{s}=0.216$ and the decoherence rate $\lambda_{d}=0.040$, and hence
positively verify the presence of scrambling. The scrambling rate is
determined purely by the underlying scrambling dynamics, and should not be
altered by the intermediate perturbation. We verified this with different
types of perturbations (2) and recovery probabilities (6). The extracted
scrambling rates are the same with a precision to the second decimal place.

To summarize, we have developed a protocol to benchmark information
scrambling. This approach distinguishes between information scrambling and
fake positive signals, produced by decoherence and operational errors in
experiments unambiguously. It can be also used to quantify the degree of
scrambling from the noisy background. Our method requires only a single loop
of forward and backward evolution, and hence can be applied to any system with
access to time-reversing the dynamics.

###### Acknowledgements.

This work was carried out under support of the U.S. Department of Energy
(DOE), Office of Science, Basic Energy Sciences, Materials Sciences and
Engineering Division, Condensed Matter Theory Program. BY also acknowledges
support from the Center for Nonlinear Studies and the U.S DOE under the LDRD
program in Los Alamos. JH was supported by the U.S. DOE through a quantum
computing program sponsored by the Los Alamos National Laboratory (LANL)
Information Science & Technology Institute. This research used quantum
computing resources provided by the LANL Institutional Computing Program,
which is supported by the U.S. DOE National Nuclear Security Administration
under Contract No. 89233218CNA000001.

## References

  * Hayden and Preskill (2007) P. Hayden and J. Preskill, “Black holes as mirrors: quantum information in random subsystems,” J. High Energy Phys. 2007, 120 (2007).
  * Kitaev (2015) A. Y. Kitaev, “A simple model of quantum holography, in proceedings of the kitp program: Entanglement in strongly-correlated quantum matter,” KITP, Santa Barbara (2015).
  * Maldacena  _et al._ (2016) J. Maldacena, S. H. Shenker, and D. Stanford, “A bound on chaos,” J. High Energy Phys. 2016, 106 (2016).
  * Rozenbaum  _et al._ (2017) E. B. Rozenbaum, S. Ganeshan, and V. Galitski, “Lyapunov exponent and Out-of-Time-Ordered correlator’s growth rate in a chaotic system,” Phys. Rev. Lett. 118, 086801 (2017).
  * Roberts and Yoshida (2017) D. A. Roberts and B. Yoshida, “Chaos and complexity by design,” J. High Energy Phys. 2017, 121 (2017).
  * Cotler  _et al._ (2017) J. Cotler, N. Hunter-Jones, J. Liu, and B. Yoshida, “Chaos, complexity, and random matrices,” J. High Energy Phys. 2017, 48 (2017).
  * Lin and Motrunich (2018) C-J Lin and O. I. Motrunich, “Out-of-time-ordered correlators in a quantum ising chain,” Phys. Rev. B 97, 144304 (2018).
  * García-Mata  _et al._ (2018) I. García-Mata, M. Saraceno, R. A. Jalabert, A. J. Roncaglia, and D. A. Wisniacki, “Chaos signatures in the short and long time behavior of the Out-of-Time ordered correlator,” Phys. Rev. Lett. 121, 210601 (2018).
  * Yan  _et al._ (2020) B. Yan, L. Cincio, and W. H. Zurek, “Information scrambling and loschmidt echo,” Phys. Rev. Lett. 124, 160603 (2020).
  * Yan and Sinitsyn (2020) B. Yan and N. A. Sinitsyn, “Recovery of damaged information and the Out-of-Time-Ordered correlators,” Phys. Rev. Lett. 125, 040605 (2020).
  * Sahu  _et al._ (2019) S. Sahu, S. Xu, and B. Swingle, “Scrambling dynamics across a Thermalization-Localization quantum phase transition,” Phys. Rev. Lett. 123, 165902 (2019).
  * Choi  _et al._ (2020) S. Choi, Y. Bao, X-L Qi, and E. Altman, “Quantum error correction in scrambling dynamics and Measurement-Induced phase transition,” Phys. Rev. Lett. 125, 030505 (2020).
  * Holmes  _et al._ (2021) Z. Holmes, A. Arrasmith, B. Yan, P. J. Coles, A. Albrecht, and A. T. Sornborger, “Barren plateaus preclude learning scramblers,” Phys. Rev. Lett. 126, 190501 (2021).
  * Landsman  _et al._ (2019) K. A. Landsman, C. Figgatt, T. Schuster, N. M. Linke, B. Yoshida, N. Y. Yao, and C. Monroe, “Verified quantum information scrambling,” Nature 567, 61–65 (2019).
  * Li  _et al._ (2017) J. Li, R. Fan, H. Wang, B. Ye, B. Zeng, H. Zhai, X. Peng, and J. Du, “Measuring Out-of-Time-Order correlators on a nuclear magnetic resonance quantum simulator,” Phys. Rev. X 7, 031011 (2017).
  * Gärttner  _et al._ (2017) M. Gärttner, J. G. Bohnet, A. Safavi-Naini, M. L. Wall, J. J. Bollinger, and A. M. Rey, “Measuring out-of-time-order correlations and multiple quantum spectra in a trapped-ion quantum magnet,” Nat. Phys. 13, 781 (2017).
  * Larkin and Ovchinnikov (1969) A. Larkin and Y. N. Ovchinnikov, “Quasiclassical method in the theory of superconductivity,” Sov. Phys., JETP 28, 1200 (1969).
  * Swingle (2018) B. Swingle, “Unscrambling the physics of out-of-time-order correlators,” Nat. Phys. 14, 988–990 (2018).
  * Swingle and Yunger Halpern (2018) B. Swingle and N. Yunger Halpern, “Resilience of scrambling measurements,” Phys. Rev. A 97, 062113 (2018).
  * Knap (2018) M. Knap, “Entanglement production and information scrambling in a noisy spin system,” Phys. Rev. B 98, 184416 (2018).
  * Syzranov  _et al._ (2018) S. V. Syzranov, A. V. Gorshkov, and V. Galitski, “Out-of-time-order correlators in finite open systems,” Phys. Rev. B 97 (2018), 10.1103/PhysRevB.97.161114.
  * Zhang  _et al._ (2019) Y-L Zhang, Y. Huang, and X. Chen, “Information scrambling in chaotic systems with dissipation,” Phys. Rev. B 99, 014303 (2019).
  * González Alonso  _et al._ (2019) J. R. González Alonso, N. Yunger Halpern, and J. Dressel, “Out-of-Time-Ordered-Correlator quasiprobabilities robustly witness scrambling,” Phys. Rev. Lett. 122, 040404 (2019).
  * Yoshida and Yao (2019) B. Yoshida and N. Y. Yao, “Disentangling scrambling and decoherence via quantum teleportation,” Phys. Rev. X 9, 011006 (2019).
  * Joshi  _et al._ (2020) M. K. Joshi, A. Elben, B. Vermersch, T. Brydges, C. Maier, P. Zoller, R. Blatt, and C. F. Roos, “Quantum information scrambling in a Trapped-Ion quantum simulator with tunable range interactions,” Phys. Rev. Lett. 124, 240505 (2020).
  * Touil and Deffner (2021) A. Touil and S. Deffner, “Information scrambling versus decoherence—two competing sinks for entropy,” PRX Quantum 2, 010306 (2021).
  * Domínguez and Álvarez (2021) F. D. Domínguez and G. A. Álvarez, “Dynamics of quantum information scrambling under decoherence effects,”  (2021), arXiv:2107.03870 [quant-ph] .
  * Zanardi and Anand (2021) P. Zanardi and N. Anand, “Information scrambling and chaos in open quantum systems,” Phys. Rev. A 103, 062214 (2021).
  * Zhao  _et al._ (2021) S. K. Zhao, Z-Y Ge, and Z. et. al. Xiang, “Probing operator spreading via floquet engineering in a superconducting circuit,”  (2021), arXiv:2108.01276 [quant-ph] .
  * Braumüller  _et al._ (2021) J. Braumüller, A. H. Karamlou, Y. Yanay, and et. al., “Probing quantum information propagation with out-of-time-ordered correlators,”  (2021), arXiv:2102.11751 [quant-ph] .
  * Mi  _et al._ (2021) X. Mi, P. Roushan, C. Quintana, and et. al., “Information scrambling in computationally complex quantum circuits,”  (2021), arXiv:2101.08870 [quant-ph] .
  * Geller  _et al._ (2021) M. R. Geller, A. Arrasmith, Z. Holmes, B. Yan, P. J. Coles, and A. Sornborger, “Quantum simulation of operator spreading in the chaotic ising model,”  (2021), arXiv:2106.16170 [quant-ph] .
  * Sinitsyn and Yan (2020) N. A. Sinitsyn and B. Yan, “The quantum butterfly non-effect,” Scientific American (2020).
  * Bradbury (1952) R. Bradbury, “A sound of thunder,” Collier’s Weekly June 28 (1952).
  * Cao and Scaffidi (2021) X. Cao and T. Scaffidi, “Origin and limit of the recovery of damaged information by time reversal,” Phys. Rev. A 103, L020401 (2021).
  * Dankert  _et al._ (2009) Christoph Dankert, Richard Cleve, Joseph Emerson, and Etera Livine, “Exact and approximate unitary 2-designs and their application to fidelity estimation,” Phys. Rev. A 80, 012304 (2009).
  * Note (1) Technically, the same result can be obtained for an ensemble of unitaries which is as random as a unitary 2-design, since the average is performed for an expression in terms of the second moment of the unitary.
  * Lewis-Swan  _et al._ (2019) R J Lewis-Swan, A Safavi-Naini, J J Bollinger, and A M Rey, “Unifying scrambling, thermalization and entanglement through measurement of fidelity out-of-time-order correlators in the dicke model,” Nat. Commun. 10, 1581 (2019).
  * Emerson  _et al._ (2005) J. Emerson, R. Alicki, and K. Życzkowski, “Scalable noise estimation with random unitary operators,” J. Opt. B 7, S347 (2005).
  * Magesan  _et al._ (2011) E. Magesan, J. M. Gambetta, and J. Emerson, “Scalable and robust randomized benchmarking of quantum processes,” Phys. Rev. Lett. 106, 180504 (2011).
  * Belyansky  _et al._ (2020) R. Belyansky, P. Bienias, Y. A. Kharkov, A. V. Gorshkov, and B. Swingle, “Minimal model for fast scrambling,” Phys. Rev. Lett. 125, 130601 (2020).
  * Swingle  _et al._ (2016) B. Swingle, G. Bentsen, M. Schleier-Smith, and P. Hayden, “Measuring the scrambling of quantum information,” Phys. Rev. A 94, 040302 (2016).
  * Yan and Zurek (2021) B. Yan and W. H. Zurek, “Decoherence factor as a convolution: an interplay between a gaussian and an exponential coherence loss,”  (2021), arXiv:2110.09463 [quant-ph] .