11institutetext: Racah Institute of Physics. The Hebrew University of
Jerusalem. Jerusalem 91904, Israel

# Type IIP supernova SN 2016X in radio frequencies

R. Ruiz-Carmona , 11<EMAIL_ADDRESS>I. Sfaradi 11 A. Horesh
11

(Received September 15, 1996; accepted March 16, 1997)

###### Abstract

Context. The study of radio emission from core-collapse supernovae (SNe)
probes the interaction of the ejecta with the circumstellar medium (CSM) and
reveals details of the mass-loss history of the progenitor.

Aims. We report observations of the type IIP supernova SN 2016X during the
plateau phase, at ages between 21 and 75 days, obtained with the Karl G.
Jansky Very Large Array (VLA) radio observatory.

Methods. We modelled the radio spectra as self-absorbed synchrotron emission,
and we characterised the shockwave and the mass-loss rate of the progenitor.
We also combined our results with previously reported X-ray observations to
verify the energy equipartition assumption.

Results. The properties of the shockwave are comparable to other type IIP
supernovae. The shockwave expands according to a self-similar law $R\propto
t^{m}$ with $m=0.76\pm 0.08$, which is notably different from a constant
expansion. The corresponding shock velocities are approximately 10700 – 8000
km s-1 during the time of our observations. The constant mass-loss rate of the
progenitor is $\dot{M}=$ (7.8 $\pm$ 0.9)
$\times\,10^{-7}\,\alpha^{-8/19}\,(\epsilon_{B}/0.1)^{-1}M_{\odot}$ yr-1, for
an assumed wind velocity of 10 km s-1. We observe spectral steepening in the
optically thin regime at the earlier epochs, and we demonstrate that it is
caused by electron cooling via the inverse Compton effect. We show that the
shockwave is characterised by a moderate deviation from energy equipartition
by a factor of $\epsilon_{e}/\epsilon_{B}\approx 28$, being the second type
IIP supernova to show such a feature.

###### Key Words.:

stars: mass-loss – circumstellar matter – supernovae: general – supernovae:
individual (ASASSN-16at, SN 2016X) – galaxies: individual (UGC 08041) – radio
continuum: general

##  1 Introduction

Massive stars experience mass loss via stellar winds or due to eruptive events
that can produce enhanced mass-loss rates towards the end of their lives
(Foley et al., 2007; Pastorello et al., 2007; Smith et al., 2010; Margutti et
al., 2014; Moriya et al., 2014; Ofek et al., 2014; Bruch et al., 2020;
Strotjohann et al., 2021). Mass loss enriches the circumstellar medium (CSM),
affects the subsequent evolution of the star, and influences the phenomenology
of some supernovae depending on the amount of mass lost at the moment of
explosion (e.g. Smith et al. 2009; van Loon 2010).

These massive stars die in core-collapse supernova (SN) explosions, once their
cores cannot be supported by thermonuclear reactions (see Smartt 2015 for a
review). SN explosions eject large amounts of mass at high velocities, that
collide with the pre-existing CSM and produce shocks, which in turn accelerate
particles to relativistic velocities and enhance the local magnetic field. As
a result, synchrotron radiation is generated, which can be prominent in radio
and X-ray wavelengths (Chevalier, 1982). Conversely, this emission can be used
to constrain the mass-loss history of the progenitor, and this provides key
information on the latest stages of stellar evolution.

The radio emission from SNe has been theoretically modelled (Chevalier, 1982,
1998; Weiler et al., 2002; Chevalier et al., 2006), and observationally
studied in several cases (e.g. Weiler et al. 1990; Chevalier 1998; Stockdale
et al. 2003; Ryder et al. 2004; Soderberg et al. 2005; Horesh et al. 2013a;
Margutti et al. 2014; Corsi et al. 2014; Bietenholz et al. 2021). Radio
observations at early phases after the SN explosion, in combination with X-ray
observations, have revealed properties of the CSM shock, provided estimates of
the mass-loss rate of the progenitor, and constrained the fraction of post-
shock energy converted to particle acceleration and magnetic field
amplification (e.g. Soderberg et al. 2010; Horesh et al. 2013a, b; Margutti et
al. 2014; Yadav et al. 2014; Kamble et al. 2016; Horesh et al. 2020).

Type IIP SNe (SNe IIP) are the most abundant SN explosions; over half of all
core-collapse SNe in several volume-limited samples are SNe IIP (Smartt 2009
and references therein). SNe IIP exhibit a long plateau phase in the optical
light curve, associated with the hydrogen envelope of the progenitor star. The
progenitors of SNe IIP are red supergiant stars (RSGs; Chevalier et al. 2006).
RSGs are characterised by masses in the range 8–25 M⊙ and radii in the range
500–1500 R⊙ (Heger et al. 2003; Hirschi et al. 2004, Chevalier 1976; Arnett
1980), although pre-explosion imaging has only identified progenitors towards
the lower mass so far (Van Dyk et al., 2003; Li et al., 2005; Maund & Smartt,
2005). RSGs are believed to sustain slow winds that drive moderate mass loss
at rates $\dot{M}=10^{-6}$–$10^{-4}\,\rm{M}_{\odot}$ with a wind velocity of
approximately $w\approx 10$ km s-1 (Smith, 2014). However, the mass loss from
these stars has been measured empirically in a small number of cases.

In this paper we present radio observations of the type IIP supernova SN
2016X. We model the radio spectra with the standard SN–CSM interaction model
(Chevalier, 1981) and derive the properties of the shock and the mass-loss
rate of the progenitor.

The paper is structured as follows. In Section 2 we summarise the previously
published results on SN 2016X. In Section 3 we present our radio observations
in detail. In Section 4 we then describe the SN–CSM interaction model and
analyse the radio spectra to characterise the shock and calculate the mass-
loss rate of the progenitor. Electron cooling and its effects are discussed in
Section 5, while in Section 6 we discuss the possible deviation from the
energy equipartition assumption. Finally, we discuss our results and summarise
our conclusions in Section 7.

##  2 SN 2016X

SN 2016X was discovered on 20.6 January 2016 by the All Sky Automated Survey
for SuperNovae (ASASSN-16at; Bock et al. 2016), and classified as a SN IIP
using optical spectroscopy acquired within the first three days after
discovery (Hosseinzadeh et al., 2016; Zheng & Zhang, 2016). The adopted time
of the explosion is 18.9 January 2016 (MJD = 57405.92 $\pm$ 0.57; Huang et al.
2018), based on a non-detection ($V>18.0$ mag) on 18.4 January and a first
detection ($V\sim 16.6$ mag) on 19.5 January. The host galaxy, UGC 08041, is
at a distance of 15.2 Mpc (Tully-Fisher distance, Sorce et al. 2014), and we
adopt this as the distance to SN 2016X.

Huang et al. (2018) and Bose et al. (2019) present high-cadence optical and
ultraviolet light curves and spectra from the early phases after explosion
until the nebular phase. The radius of the progenitor star, as revealed by the
rise time, is in the range 860 – 990 R⊙ (Gall et al., 2015). Using mass-radius
relations for RSGs, the mass of the progenitor is estimated to be 18.5–19.7
M⊙. The quasi-bolometric luminosity is comparable to other SNe IIP, but it
declines quickly during the plateau phase. The study of absorption features in
the spectral lines H$\alpha$, H$\beta$, and Fe ii 5169Å, 5018Å, and 4924Å
shows that photospheric velocities are high compared to other SNe IIP such as
SN 1999em and SN 2005cs, but are comparable to SN 2013ej and SN 2014cx
(Leonard et al., 2002; Pastorello et al., 2009; Huang et al., 2015, 2016),
evolving rapidly from $\sim$9000 km s-1 at $t$ = 20 days to $\sim$6000 km s-1
at $t$ = 80 days for H$\alpha$.

SN 2016X was also observed in X-rays with the Neil Gehrels Swift Observatory
and the Chandra X-ray Observatory in the 0.3–10 keV energy band. The observed
X-ray luminosity was $L_{X}\sim$ 2.5 $\times$ $10^{39}$ erg s-1 at $t\sim$2–5
days, and declined by a factor of 5 by $t\sim$ 19 days (Bose et al., 2019).
The count rate is too low to perform any spectral analysis on the X-ray
observations.

Finally, it is noteworthy that SN 2016X showed double-peaked Balmer lines with
strong asymmetries during the nebular phase. This phenomenology has been
associated with bipolar ejecta or aspherical distribution of the CSM (Bose et
al., 2019; Utrobin & Chugai, 2019), and possibly with the formation of dust in
the ejecta (Inserra et al., 2011; Bose et al., 2019).

##  3 Observations

We observed SN 2016X ($\alpha$ = 12h 55m 15${}^{\rm{s}}.61$, $\delta$ = 12∘
05’ 59”.35, J2000) with the National Radio Astronomy Observatory (NRAO) Karl
G. Jansky Very Large Array (VLA) under the programme 15B-395 (PI A. Horesh).
Observations were acquired at four epochs: 9.50 February, 15.47 February, 2.56
March, and 3.84 April 2016, which correspond to 21.58, 27.55, 43.54, and 75.42
days after the explosion, adopting the ephemeris of the explosion by Bose et
al. (2019). Observations were carried out in bands S (2–4 GHz), C (4–8 GHz), X
(8–12 GHz), Ku (12–18 GHz), and K (18–26 GHz) in the three later epochs, while
only C- and K-band observations were obtained in the first epoch. The antenna
array was in configuration C.

We reduced the observational data with the Common Astronomy Software
Applications package (casa; McMullin et al. 2007). The reduction was done
manually, as additional flagging was needed in most bands. The radio sources J
1246–0730 and 3C 286 were used as phase and flux calibrators, respectively. We
imaged the field of SN 2016X in several sub-bands depending on the signal-to-
noise ratio of the detection and the predominance of radio frequency
interference (RFI) using the task clean interactively. We extracted the flux
of the source by fitting a two-dimensional Gaussian with the task imfit, and
we added in quadrature the contribution of three factors to estimate the error
bars on the flux: (a) the error provided by the fitting routine, (b) the RMS
of the image estimated in an empty region around the source, and (c) 10% of
the peak flux to account for calibration errors. We adopted three times the
RMS of the image as a 3$\sigma$ upper limit on the flux for non-detections. In
the S-band images, SN 2016X appears embedded in the side lobe of a bright
source in the field. We do not use the S-band data in this work since the
extracted flux in that band is artificially high. A summary of the
observations is presented in Table 1, and images at 4.8 and 22.0 GHz are
represented in Figure 1.

Figure 1: Images of SN 2016X at the four epochs of our observations, obtained
at 4.8 GHz (top row) and 22.0 GHz (bottom row). The flux on the images is
indicated by the colour bars on the right, and the beam size is represented at
the bottom left of every image. Table 1: Summary of radio observations of SN
2016X.

Epoch | Central Frequency | Bandwidth | Flux
---|---|---|---
(days) | (GHz) | (GHz) | ($\mu$Jy)
+21.6 | 4.8 | 1.0 | 220 $\pm$ 41
+21.6 | 7.4 | 1.0 | 659 $\pm$ 72
+21.6 | 18.8 | 1.7 | 659 $\pm$ 77
+21.6 | 22.0 | 1.5 | 454 $\pm$ 67
+21.6 | 23.6 | 1.6 | 519 $\pm$ 79
+21.6 | 25.2 | 1.7 | 495 $\pm$ 61
+27.6 | 4.8 | 1.0 | 404 $\pm$ 47
+27.6 | 7.4 | 1.0 | 739 $\pm$ 73
+27.6 | 8.4 | 0.8 | 735 $\pm$ 81
+27.6 | 9.1 | 0.8 | 631 $\pm$ 78
+27.6 | 9.9 | 0.9 | 628 $\pm$ 66
+27.6 | 10.9 | 1.2 | 575 $\pm$ 68
+27.6 | 13.3 | 1.0 | 457 $\pm$ 53
+27.6 | 14.3 | 1.0 | 455 $\pm$ 56
+27.6 | 15.3 | 1.0 | 402 $\pm$ 52
+27.6 | 16.3 | 1.0 | 394 $\pm$ 54
+27.6 | 17.4 | 1.3 | 341 $\pm$ 48
+27.6 | 22.0 | 8.0 | 233 $\pm$ 46
+43.5 | 4.8 | 1.0 | 537 $\pm$ 65
+43.5 | 7.4 | 1.0 | 475 $\pm$ 57
+43.5 | 8.4 | 0.8 | 443 $\pm$ 58
+43.5 | 9.1 | 0.8 | 437 $\pm$ 58
+43.5 | 9.9 | 0.9 | 345 $\pm$ 42
+43.5 | 10.9 | 1.2 | 298 $\pm$ 44
+43.5 | 13.6 | 1.7 | 220 $\pm$ 32
+43.5 | 15.4 | 1.7 | 211 $\pm$ 33
+43.5 | 17.1 | 1.8 | 174 $\pm$ 32
+43.5 | 22.0 | 8.0 | 145 $\pm$ 26
+75.4 | 4.8 | 1.0 | 284 $\pm$ 45
+75.4 | 7.4 | 1.0 | 193 $\pm$ 33
+75.4 | 8.4 | 0.8 | 194 $\pm$ 41
+75.4 | 9.1 | 0.8 | 163 $\pm$ 27
+75.4 | 9.9 | 0.9 | 165 $\pm$ 31
+75.4 | 10.9 | 1.2 | 162 $\pm$ 29
+75.4 | 15.4 | 5.3 | 93 $\pm$ 16
+75.4 | 22.0 | 8.0 |  $<90$ $\pm$ 30

The epoch of the observations is estimated from the explosion time (MJD =
57405.92 $\pm$ 0.57; Huang et al. 2018). The flux is calculated fitting a 2D
Gaussian profile to the images of the sub-bands specified by their central
frequency and bandwidth. The error bars are calculated adding in quadrature
the error of the fitting routine, the RMS of the image, and 10% of the peak
flux in order to account for calibration errors.

##  4 Analysis of radio observations

We use the SN–CSM interaction model by Chevalier (1981) to study the radio
emission from SN 2016X. We make the following assumptions in our analysis: (i)
we model the radio flux as self-absorbed synchrotron radiation, as described
in Chevalier (1998) and Chevalier et al. (2006), although we note that other
absorption mechanisms such as free-free absorption by ionised gas may also
play a role in the case of SNe IIP (see Section 4.3); (ii) we assume that the
density profile of the CSM decreases with radius as $\rho_{\rm{CSM}}\propto
R^{-2}$; (iii) we assume that only half of the spherical volume inside the
radius $R$ is contributing to the radio emission at all times (i.e. the
filling factor is $f=0.5$); (iv) we do not include the effects of electron
cooling in our modelling; and (v) we assume that there is energy equipartition
between particles and fields, that is, the fraction of the post-shock energy
converted into particle acceleration, $\epsilon_{e}$, and the fraction of the
post-shock energy converted into magnetic field amplification, $\epsilon_{B}$,
are equal and approximately 0.1. Therefore, the equipartition parameter
$\alpha$ becomes unity ($\alpha\equiv\epsilon_{e}/\epsilon_{B}$ = 1). We note
that deviations from the assumptions above will affect the shockwave and mass-
loss estimates, as we investigate in Sections 5 and 6.

The energy distribution of the relativistic electrons follows a power law,
$N(E)=N_{0}E^{-\gamma}$, where $E$ is the electron energy and $N_{0}$ is a
constant. Under the energy equipartition assumption we can relate the constant
$N_{0}$ and the magnetic field strength, $B$,

$N_{0}=\frac{1}{8\pi}\alpha B^{2}(\gamma-2)E_{l}^{\gamma-2},$ (1)

where $E_{l}=0.51$ MeV is the rest mass energy of the electron.

The complete expression that describes the synchrotron emission from SNe is
given in Equation 1 of Chevalier (1998). The observable flux in the optically
thick regime is

$F_{\nu}=\frac{\pi
c_{5}}{c_{6}}\,\frac{R^{2}}{D^{2}}\,B^{-1/2}\,\left(\frac{\nu}{2c_{1}}\right)^{5/2},$
(2)

and in the optically thin regime yields

$F_{\nu}=\frac{1}{6}\,\alpha
fc_{5}\,(\gamma-2)\,E_{l}^{\gamma-2}\,\frac{R^{2}}{D^{2}}\,B^{(\gamma+5)/2}\,\left(\frac{\nu}{2c_{1}}\right)^{-(\gamma-1)/2}.$
(3)

Here $c_{1}$ is a constant (Chevalier, 1998); $c_{5}$ and $c_{6}$ are
functions of the electron energy distribution power-law index $\gamma$
(Pacholczyk, 1970); $D$ is the distance to the supernova; $R$ is the radius of
the forward shock; $B$ is the magnetic field strength; and $f$ is the filling
factor. The spectral slope in the optically thin part of the spectrum,
$\beta$, is a function of the electron energy distribution power-law index:
$\beta=-(\gamma-1)/2$.

It can be seen from Equations 2 and 3 that the transition from optically thick
to optically thin emission creates a peak in the spectrum. By characterising
the observed spectral peak by its frequency, $\nu_{\rm{peak}}$, and its flux,
$F_{\rm{peak}}$, we can obtain the following expressions for the radius of the
shock and the magnetic field strength (in cgs units, see Chevalier 1998):

$R=\frac{2c_{1}}{\nu_{\rm{peak}}}\,\left(\frac{6\pi\,c_{6}^{\gamma+5}\,(0.5\,{F}_{\rm{peak}}\,D^{2})^{\gamma+6}\,}{c_{5}^{\gamma+6}\alpha\,f\,(\gamma-2)\,E_{l}^{\gamma-2}}\right)^{1/(2\gamma+13)},$
(4)

$B=\frac{\nu_{\rm{peak}}}{2c_{1}}\,\left(\frac{706\,c_{5}}{\left(\alpha\,f\,(\gamma-2)\,E_{l}^{\gamma-2}\right)^{2}\,c_{6}^{3}\,F_{\rm{peak}}\,D^{2}}\right)^{2/(2\gamma+13)}.$
(5)

The mass-loss rate of the progenitor, $\dot{M}$, divided by the wind velocity,
$w$, can be calculated using the fact that the magnetic energy density is a
fraction $\epsilon_{B}$ of the post-shock energy,

$\frac{B^{2}}{8\pi}=\frac{9}{8}\,\epsilon_{B}\,\rho_{\rm{CSM}}\,V_{\rm{sh}}^{2},$
(6)

where $\rho_{\rm{CSM}}=\dot{M}/4\pi\,w\,R^{2}$ is the density of the CSM, and
$V_{\rm{sh}}$ is the expansion velocity of the shock. Assuming constant
expansion velocity ($V_{\rm{sh}}=R/t$), and using Equation 5 and 6, we obtain

$\frac{\dot{M}}{w}=\frac{1}{9\,\epsilon_{B}\,c_{1}^{2}}\left(\frac{706\,c_{5}}{c_{6}^{3}\left(\alpha
f(\gamma-2)E_{l}^{\gamma-2}\right)^{2}\,F_{\rm{peak}}\,D^{2}}\right)^{4/2\gamma+13}\nu_{\rm{peak}}^{2}\,t^{2}.$
(7)

###  4.1 Single-epoch spectral modelling

In this section we analyse the radio spectrum for individual epochs. The free
parameters of the SN–CSM interaction model (Chevalier, 1998) are the power-law
index of the electron energy distribution, $\gamma$; the radius of the forward
shock, $R$; and the magnetic field strength, $B$. Theoretical considerations
and observations of other SNe IIP indicate that $\gamma\approx 3.0$ (e.g.
Chevalier 1998 and references within). This implies that the spectral index in
the optically thin regime equals $\beta=-1$, although steepening of the
spectral slope is possible due to electron cooling. Since the simple SN–CSM
interaction model described in Chevalier (1998) does not account for electron
cooling, any modelling of a spectral index steeper than $\beta=-1$ when
cooling is in effect may incorrectly point to large values of the electron
energy distribution power-law index, $\gamma$. Since our analysis suggests
that electron cooling is important for SN 2016X, as we demonstrate in Section
5, we fix the electron energy distribution power-law index to $\gamma=3.0$,
and fit for the optically thin spectral index, $\beta$, instead. Transforming
Equation 1 of Chevalier (1998), our model reads

$F_{\nu}=\frac{\pi
c_{5}}{c_{6}}\,\frac{R^{2}}{D^{2}}\,B^{-1/2}\,\left(\frac{\nu}{2c_{1}}\right)^{5/2}\left[1-\exp\left(-\frac{\nu}{\nu_{1}}\right)^{\beta-5/2}\right],$
(8)

where $\nu_{1}$ is the frequency at which the optical depth equals unity for
self-absorbed synchrotron:

$\nu_{1}=2\,c_{1}\left(\frac{c_{6}}{6\pi}\,\alpha\,f\,E_{l}\,R\,B^{9/2}\right)^{2/7}.$
(9)

We implement a Markov chain Monte Carlo (MCMC) algorithm sampling over the
optically thin spectral index, $\beta$, the radius of the shock, $R$, and the
magnetic field strength, $B$. We determine the best values and the error bars
of the free parameters from the posterior probability distributions. We also
estimate the confidence intervals of the frequency and flux of the spectral
peak from the Markov chains. The velocity of the shock is calculated assuming
constant expansion, $V_{\rm{sh}}=R/t$, and the ratio of the mass-loss rate to
the wind velocity is calculated with Equation 7. The results of the analysis
of individual epochs are gathered in Table 2, and the best-fit models are
represented in Figure 2.

We note that the optically thick part of the spectrum and the spectral peaks
in the later epochs, especially the fourth one, are not well constrained by
the modelling, which results in larger error bars in the estimates of the
physical quantities. We observe a steepening of the spectral slope $\beta$ in
the first three epochs, while the spectral index approaches the expected value
$\beta\approx-1$ in the last epoch, which we interpret as a consequence of
electron cooling (see Section 5). The radius of the shock expands from
$R\approx 2.6\,\times\,10^{15}$ cm at $t=21.6$ days to $R\approx
5.0\,\times\,10^{15}$ cm at $t=75.4$ days. The resulting shock velocities,
assumed constant in this section ($V_{\rm{sh}}=R/t$), are in the range 13900 –
7800 km s-1, and are higher than for other SNe IIP observed at ages between 8
and 81 days, which typically show expansion velocities below 10000 km s-1
(e.g. Beswick et al. 2005; Chevalier et al. 2006; Yadav et al. 2014). The high
expansion velocities we obtained can indicate that the free expansion phase,
where $R\propto t$, does not extend to the time of our observations. Our
estimates of the shock radius are compatible with a power law $R\propto t^{m}$
with index m = 0.8 $\pm$ 0.1, excluding the last epoch. Finally, we obtain a
weighted average mass-loss rate of $\dot{M}$ = (4.4 $\pm$ 0.3)
$\times\,10^{-7}\,M_{\odot}$ yr-1, for an assumed wind velocity of $w=10$ km
s-1.

Table 2: Results of the single-epoch spectral analysis.

Epoch | $\chi^{2}$ | dof | $\beta$ | $\nu_{\rm{peak}}$ | $F_{\rm{peak}}$ | Magnetic | Shock | Shock | $\dot{M}/w$
---|---|---|---|---|---|---|---|---|---
field | radius | velocity
(days) |  |  |  | (GHz) | (mJy) | (G) | ($10^{15}$ cm) | ($10^{3}$ km s-1) | $\left(\rm{10^{-8}\frac{M_{\odot}\,yr^{-1}}{km\,s^{-1}}}\right)$
+21.6 | 2.2 | 3 | $-1.26^{+0.40}_{-0.50}$ | $11.0^{+0.4}_{-0.5}$ | $0.93^{+0.22}_{-0.16}$ | $1.3^{+0.1}_{-0.1}$ | $2.6^{+0.1}_{-0.1}$ | $13.9^{-0.6}_{-0.6}$ | $4.9^{+0.4}_{-0.4}$
+27.6 | 1.2 | 9 | $-1.35^{+0.16}_{-0.16}$ | $7.6^{+0.3}_{-0.3}$ | $0.72^{+0.04}_{-0.04}$ | $1.0^{+0.1}_{-0.1}$ | $3.2^{+0.2}_{-0.2}$ | $13.5^{+0.7}_{-0.7}$ | $4.1^{+0.3}_{-0.3}$
+43.5 | 1.9 | 7 | $-1.25^{+0.19}_{-0.20}$ | $5.3^{+0.5}_{-0.8}$ | $0.57^{+0.06}_{-0.05}$ | $0.7^{+0.1}_{-0.1}$ | $4.2^{+1.1}_{-0.6}$ | $11.2^{+2.9}_{-1.4}$ | $5.1^{+1.1}_{-1.5}$
+75.4 | 1.9 | 4 | $-1.05^{+0.21}_{-0.25}$ | $3.5^{+1.4}_{-1.0}$ | $0.32^{+0.10}_{-0.06}$ | $0.5^{+0.2}_{-0.1}$ | $5.0^{+3.0}_{-1.8}$ | $7.8^{+4.7}_{-2.9}$ | $7.7^{+7.6}_{-3.9}$

We use a model with $\gamma=3$ and with free parameters for the optically thin
spectral index, $\beta$; the radius of the shock, $R$; and the magnetic field,
$B$. The columns in the table (from left to right) are the epoch of the
observations; the $\chi^{2}$; the number of degrees of freedom of the fitting;
the spectral slope in the optically thin regime, $\beta$; the frequency and
flux of the observed spectral peak, $\nu_{\rm{peak}}$ and $F_{\rm{peak}}$; the
magnetic field strength, $B$; the radius of the shock $R$; the shock expansion
velocity assuming constant velocity $V_{\rm{sh}}$; and the ratio of the mass-
loss rate to the velocity of the wind, $\dot{M}/w$ (Equation 7).

Figure 2: Single-epoch spectral modelling of SN2016X. The shaded regions
represent the 1$\sigma$ confidence intervals of the best models.

###  4.2 Multiple-epoch spectral modelling

In this section we analyse the radio spectra at all epochs including
expressions for the temporal evolution of the radius of the shock and the
magnetic field. We assume that the electron energy distribution power-law
index ($\gamma=3$) and the filling factor ($f=0.5$) remain constant for the
entire evolution of the supernova, as observations have shown in other SNe
(e.g. Weiler et al. 1991).

We introduce a self-similar solution for the evolution of the radius of shock,
$R=R_{0}\,(t/t_{0})^{m}$, where $m$ is a number between 2/3 and 1 (Chevalier,
1982). We also introduce a self-similar solution for the evolution of the
magnetic field: $B=B_{0}\,(t/t_{0})^{-b}$. The value of the exponent $b$ is
indicative of the amplification mechanism of the magnetic field: $b\approx 1$
points towards turbulent motions in the shocked regions, while $b\approx m$
implies compression of the circumstellar magnetic fields (Chevalier, 1982).

We implement a Markov chain Monte Carlo algorithm, sampling over the spectral
slope at individual epochs, $\beta$, the parameters describing the temporal
evolution of the radius of the shock, $R_{0}$ and $m$, and the parameters
describing the temporal evolution of the magnetic field, $B_{0}$ and $b$. We
evaluate $R_{0}$ and $B_{0}$ at the first epoch of our observations by fixing
$t_{0}$ to $t=21.6$ days. We estimate the values and the error bars of the
free parameters from the posterior probability distributions, for which we
obtain $R_{0}=(2.6\pm 0.1)\,\times\,10^{15}$ cm, $m=0.76^{+0.09}_{-0.07}$,
$B_{0}=(1.25\pm 0.06)$ G, and $b=1.00\pm 0.06$ (with a $\chi^{2}=7.8$ and
degrees of freedom = 27). We also derive the features of the spectral peaks at
individual epochs from the Markov chains. We estimate the magnetic field
strength, the shock radius and the expansion velocity from their temporal
evolution ($V_{\rm{sh}}=dR/dt=R_{0}\,m\,t^{m-1}$ for the shock velocity). We
calculate the ratio of the mass-loss rate to the wind velocity with a modified
version of Equation 7 that includes a factor $m^{-2}$ to account for the
evolution of the shock radius. The results of the multiple-epoch analysis are
summarised in Table 3, and the best-fit model is represented in Figure 3.

The results of the multiple-epoch analysis are similar to those in Section
4.1, except for the slower expansion velocities in the range 10700 – 8000 km
s-1. The temporal evolution of the magnetic field ($b\approx 1$) suggests that
turbulent motions are responsible for the amplification mechanism. We obtain a
weighted average mass-loss rate of $\dot{M}=$ (7.8 $\pm$ 0.9)
$\times\,10^{-7}\,M_{\odot}$ yr-1, taking into account the error bars and for
an assumed wind velocity of $w=10$ km s-1. We also observe spectral steepening
during the first three epochs (see Section 5). Since the spectral peaks are
more tightly constrained with the multiple-epoch analysis and the error bars
are generally lower, we use the results of this section in Sections 5 and 6.

Table 3: Results of the multiple-epoch spectral analysis.

Epoch | $\beta$ | $\nu_{\rm{peak}}$ | $F_{\rm{peak}}$ | Magnetic | Shock | Shock | $\dot{M}/w$
---|---|---|---|---|---|---|---
field | radius | velocity
(days) |  | (GHz) | (mJy) | (G) | ($10^{15}$ cm) | ($10^{3}$ km s-1) | $\left(\rm{10^{-8}\frac{M_{\odot}\,yr^{-1}}{km\,s^{-1}}}\right)$
+21.6 | $-1.08^{+0.15}_{-0.15}$ | $10.7^{+0.4}_{-0.3}$ | $0.86^{+0.05}_{-0.05}$ | $1.25^{+0.06}_{-0.06}$ | $2.6^{+0.1}_{-0.1}$ | $10.7^{+1.1}_{-0.9}$ | $8.1^{+1.6}_{-1.4}$
+27.6 | $-1.42^{+0.14}_{-0.15}$ | $7.8^{+0.2}_{-0.2}$ | $0.73^{+0.03}_{-0.03}$ | $0.98^{+0.04}_{-0.04}$ | $3.2^{+0.1}_{-0.1}$ | $10.1^{+1.2}_{-1.0}$ | $7.2^{+1.6}_{-1.4}$
+43.5 | $-1.15^{+0.10}_{-0.10}$ | $5.0^{+0.2}_{-0.2}$ | $0.55^{+0.03}_{-0.03}$ | $0.62^{+0.03}_{-0.03}$ | $4.5^{+0.3}_{-0.2}$ | $9.1^{+1.4}_{-1.1}$ | $7.9^{+2.1}_{-1.9}$
+75.4 | $-1.02^{+0.13}_{-0.13}$ | $2.9^{+0.2}_{-0.2}$ | $0.40^{+0.05}_{-0.04}$ | $0.36^{+0.03}_{-0.03}$ | $6.7^{+0.7}_{-0.5}$ | $8.0^{+1.7}_{-1.2}$ | $8.5^{+2.9}_{-2.4}$

We use a model with $\gamma=3$ and with free parameters for the optically thin
spectral indexes, $\beta$; the parameters of the temporal evolution of the
shock radius, $R_{0}$ and $m$ ($R=R_{0}\,(t/t_{0})^{m}$, $t_{0}$=21.6 days);
and temporal evolution of the magnetic field, $B_{0}$ and $b$
($B=B_{0}\,(t/t_{0})^{-b}$, $t_{0}$=21.6 days).The columns in the table (from
left to right) are the epoch of observation; the optically thin spectral
index, $\beta$; the frequency of the spectral peak, $\nu_{\rm{peak}}$; the
flux of the spectral peak, $F_{\rm{peak}}$; the magnetic field strength, $B$;
the radius of the shock $R$; the expansion velocity, $V_{\rm{sh}}$; and the
ratio of the mass-loss rate to the velocity of the wind, $\dot{M}/w$ (Equation
7).

Figure 3: Multiple-epoch spectral modelling of SN2016X. The shaded regions
represent the 1$\sigma$ confidence intervals of the best model.

###  4.3 Absorption mechanism

Radio emission from SNe is affected by absorption at early times. Our radio
model presented in Section 4 includes self-absorbed synchrotron radiation, but
free-free absorption by ionised gas in the CSM can also play an important
role. Here we investigate the significance of the free-free absorption for SN
2016X.

Free-free absorption has been studied in detail (e.g. Chevalier 1982; Weiler
et al. 2002; Chevalier & Fransson 2006), and examined for individual SNe on
several occasions (e.g. Horesh et al. 2012; Nayana et al. 2018). The optical
depth for free-free absorption is calculated as

$\tau_{\rm{ffa}}=0.083\,T_{\rm{e}}^{-1.35}\,\nu^{-2.1}\,\rm{EM},$ (10)

where $T_{\rm{e}}$ is the temperature of the electrons in the CSM, the
frequency is expressed in GHz, and $\rm{EM}$ is the integral of the square of
the electron density along the line of sight expressed in cm-6 pc, known as
emission measure. The EM from a radius, $R$, to infinity can be calculated as

$\rm{EM}=\int_{R}^{\infty}n_{0}^{2}\left(\frac{r}{R_{0}}\right)^{-4}\,dr=\frac{1}{3}\,n_{0}^{2}\,R_{0}\left(\frac{R}{R_{0}}\right)^{-3},$
(11)

where $n_{0}$ is the electron density at radius $R_{0}$.

Expressing $n_{0}$ as a function of $\dot{M}/w$ and substituting the radius by
its time evolution solution found in Section 4.2 ($R=R_{0}(t/t_{0})^{m}$), the
free-free optical depth becomes

$\tau_{\rm{ffa}}=8.1\times
10^{16}\,\left(\frac{\dot{M}}{w}\right)^{2}\,\nu^{-2.1}\,t^{-2.28},$ (12)

where we use $T_{e}=10^{5}$ K (Lundqvist & Fransson 1988, Chevalier et al.
2006). For our first epoch, we obtain
$\tau_{\rm{ffa}}[\nu=4.8,8.4,22.0\,\rm{GHz}]=[0.018,0.006,0.001]$, which shows
that the free-free optical depth is not relevant for our observations of SN
2016X.

##  5 Electron cooling

Relativistic electrons can undergo synchrotron or inverse Compton cooling, and
affect the spectra and light curves of SNe. In particular inverse Compton
cooling can be important for some SNe IIP during the plateau phase due to the
high density of photospheric photons (Chevalier et al., 2006). Our spectral
analysis in Section 4 suggests that electron cooling may be significant for SN
2016X; therefore, we study the characteristic timescales and frequencies of
the cooling mechanisms in this section.

First, we compare the timescales of synchrotron and the inverse Compton
cooling, $t_{\rm{synch}}$ and $t_{\rm{IC}}$, to the timescale of the expansion
of the radio-emitting shell, $t_{\rm{SN}}$. We use Equations 5 and 6 from
Chevalier et al. (2006),

$t_{\rm{synch}}/t_{\rm{SN}}=2.0\times
10^{-6}\,\left(\epsilon_{B}\,\frac{\dot{M}}{w}\right)^{-3/4}\,\left(\frac{t}{\nu}\right)^{1/2},$
(13)

$t_{\rm{IC}}/t_{\rm{SN}}=3.6\times
10^{35}\,\frac{V_{\rm{sh}}^{2}}{L_{\rm{bol}}}\left(\epsilon_{B}\,\frac{\dot{M}}{w}\right)^{1/4}\,\left(\frac{t}{\nu}\right)^{1/2},$
(14)

where $L_{\rm{bol}}$ is the bolometric luminosity of the supernova in erg s-1,
the velocity of the shock is expressed in km s-1, the ratio of the mass-loss
rate to the wind velocity is expressed in $M_{\odot}$ yr-1/ km s-1, times are
expressed in days, and frequencies are expressed in GHz. We approximate the
bolometric light curve presented in Huang et al. (2018) by a power law,
$L_{\rm{bol}}(t)\propto t^{\delta}$, with $\delta=-0.59$, between phases $t$ =
30 days ($L_{\rm{bol}}(30)=6.3\times 10^{41}$ erg s-1) and $t$ = 80 days
($L_{\rm{bol}}(80)=3.5\,\times\,10^{41}$ erg s-1). Assuming
$\epsilon_{B}=0.1$, we obtain $t_{\rm{synch}}/t_{\rm{SN}}\approx$ 2.4 – 9.5
and $t_{\rm{IC}}/t_{\rm{SN}}\approx$ 0.4 – 1.7 for the frequencies and the
times of our observations ($\nu$ = 4.8 – 22.0 GHz, $t$ = 21.6 –75.4 days). We
note that the shorter cooling timescales are obtained at earlier times and
higher frequencies. These timescales indicate that inverse Compton cooling can
be important for SN 2016X, which in turn may indicate low values of the mass-
loss rate, $\dot{M}$, and $\epsilon_{B}$ (Chevalier et al., 2006).

Second, we estimate the cooling frequencies for both mechanisms. The
synchrotron cooling frequency is

$\nu_{\rm{synch}}=1.8\times
10^{-8}\,\pi\,\frac{c\,m_{e}\,q_{e}}{\sigma_{T}^{2}}\,B^{-3}\,t^{-2}\quad\rm{GHz},$
(15)

where $m_{e}$ and $q_{e}$ are the mass and charge of the electron, $c$ is the
speed of light, $\sigma_{T}$ is the Thomson cross-section, and all variables
are expressed in cgs units. The inverse Compton cooling frequency is

$\nu_{\rm{IC}}=1.3\times
10^{71}\,\frac{V_{\rm{sh}}^{4}}{L_{\rm{bol}}^{2}}\left(\epsilon_{B}\,\frac{\dot{M}}{w}\right)^{1/2}\quad\rm{GHz},$
(16)

where the bolometric luminosity is expressed in erg s-1, the velocity of the
shock is expressed in km s-1, and the ratio of the mass-loss rate to the wind
velocity is expressed in $M_{\odot}$ yr-1/ km s-1.

We obtain characteristic synchrotron cooling frequencies above
$\nu_{\rm{synch}}\gtrsim$ 200 GHz for all our epochs, so our observations are
not affected by synchrotron cooling. We obtain inverse Compton cooling
frequencies in the range $\nu_{\rm{IC}}$ = 2.8 – 14.7 GHz, which can affect
our observations. We interpret the steepening of the optically thin spectral
index at the earlier epochs, especially at $t=27.6$ days for which
$\beta\approx-1.4$, as the signature of inverse Compton cooling in the radio
spectrum.

##  6 X-ray emission and energy equipartition

In this section we study the origin of the X-ray emission from SN 2016X, and
we combine the X-ray data with our radio analysis results to verify the energy
equipartition assumption. We use the X-ray luminosity registered by Chandra,
$L_{X}=4.7^{+1.4}_{-1.7}\,\times\,10^{38}$ erg s-1 in the 0.3–10 keV band at
phase $t=18.9$ days (Bose et al., 2019), and we assume that the X-ray
luminosity remains constant until the time of our first radio observations at
$t=21.6$ days.

First, we consider thermal emission from the circumstellar and reverse shocks
as the source of the X-ray emission. We estimate the temperatures behind the
circumstellar shock, $T_{\rm{cs}}$, and the reverse shock, $T_{\rm{rs}}$, as
(Fransson et al., 1996)

$T_{\rm{cs}}=22.7\,m^{2}\,V_{\rm{sh}}^{2}=1.5^{+0.7}_{-0.5}\times
10^{9}\,\rm{K}\qquad[m=0.76],$ (17)
$T_{\rm{rs}}=\left(1-\frac{1}{m}\right)^{2}\,T_{\rm{cs}}=1.5^{+3.0}_{-1.1}\times
10^{8}\,\rm{K}\quad\;[m=0.76],$ (18)

where $m$ is the exponent of the evolution of the shock radius ($R\propto
t^{m}$), and $V_{\rm{sh}}$ is expressed in km s-1.

With these temperatures, we calculate the X-ray emission per unit energy
(Equations 8 and 9 from Chevalier et al. 2006),

$\frac{dL_{\rm{cs}}}{dE}=4.0\times
10^{47}\,\frac{E_{X}^{-0.23}\,T_{\rm{cs}}^{-0.08}}{\exp{(0.0116\,E_{X}/T_{\rm{cs}})}}\,\left(\dfrac{\dot{M}}{w}\right)^{2}V_{\rm{sh}}^{-1}\,t^{-1},$
(19)

$\frac{dL_{\rm{rs}}}{dE}=5.9\times
10^{48}\,\frac{T_{\rm{rs}}^{-0.024}}{\exp{(0.0116/T_{\rm{rs}}})}\,\left(\dfrac{\dot{M}}{w}\right)^{2}V_{\rm{sh}}^{-1}\,t^{-1}$
(20)

in erg s-1 keV-1 units. Here $E_{X}=4.9$ keV is the average energy of the
X-ray photon, the temperatures are expressed in $10^{9}$ K units, the ratio of
the mass-loss rate to the wind velocity is expressed in $M_{\odot}$ yr-1/ km
s-1, the velocity of the shock is expressed in km s-1, and the age of the SN
is expressed in days. Adding these two contributions, we obtain a thermal
X-ray emission of $\sim 6.8\times 10^{35}$ erg s-1 keV-1, for which the
uncertainty is dominated by the estimates of the temperatures. This can
roughly account for 1% of the X-ray emission observed with Chandra.

X-ray thermal emission can be absorbed by a dense shell of cool low-ionised
plasma in the contact discontinuity region if the reverse shock is radiative
(Fransson et al., 1996). Radiative reverse shocks are plausible when the
supernova presents large density gradients (i.e. $n\gtrsim 9$ with
$\rho_{\rm{SN}}\propto R^{-n}$), but we find $n=(2m-3)/(m-1)\approx 6.2$ for
SN 2016X. Moreover, the presence of a layer of cool absorbing plasma also
requires short cooling times for free-free emission behind the reverse shock;
however, for SN 2016X we find

$t_{\rm{cool}}=5\times
10^{-19}\,\left(\dfrac{\dot{M}}{w}\right)^{-1}V_{\rm{sh}}^{3}\,t^{2}=3650\,\rm{days}.$
(21)

We therefore conclude that absorption of X-ray emission from the reverse shock
is not relevant for SN 2016X.

Since the thermal emission cannot explain the observed X-ray flux alone, we
consider emission via inverse Compton scattering of photospheric photons. The
X-ray emission per unit energy due to inverse Compton scattering is (Equation
11 from Chevalier et al. 2006)

$\frac{dL_{\rm{IC}}}{dE}=8.8\,\Gamma_{\rm{min}}\,\epsilon_{e}\,E_{X}^{-1}\,L_{\rm{bol}}\,\frac{\dot{M}}{w}\,V_{\rm{sh}}\,t^{-1}\quad\rm{erg}\,\rm{s}^{-1}\,\rm{keV}^{-1},$
(22)

where $\Gamma_{\rm{min}}$ is the minimum Lorentz factor of the electron
population, and $L_{\rm{bol}}=10^{42}$ erg s-1 is the bolometric luminosity at
$t=21.6$ days (Huang et al., 2018). We assume that the minimum Lorentz factor
is $\Gamma_{\rm{min}}=1$ as the population of non-relativistic electrons can
be significant for all subtypes of SNe II (Chevalier & Fransson, 2006). We
obtain an X-ray luminosity via inverse Compton scattering of $6.2\times
10^{36}$ erg s-1 keV-1, which accounts for approximately 14% of the observed
emission with Chandra. We note that the X-ray luminosity could have decreased
at the epoch of our first radio observations. We also note that these
calculations are sensitive to the energy equipartition assumption as a
deviation from equipartition leads to an underestimation of the mass-loss rate
divided by the wind velocity.

Hence, we investigate potential deviations from energy equipartition combining
the X-ray and radio observations. From Equations 5 and 6 it can be seen that

$S_{*}\equiv\frac{\dot{M}}{w}\,\epsilon_{B}\,\alpha^{8/19}=2.4\times
10^{18}\,\left(\frac{\nu_{\rm{peak}}\,t}{m}\right)^{2}\,\left(F_{\rm{peak}}\,D^{2}\right)^{-4/19}.$
(23)

Using the results of Section 4.2, we obtain $S_{*}=5.2\times 10^{12}$ g cm-1.
Assuming that inverse Compton cooling is responsible for the observed X-ray
flux, we can rewrite Equation 22 substituting the mass-loss rate term by
$S_{*}$,

$E_{X}\,\frac{dL_{\rm{L_{X}}}}{dE}=8.8\,\Gamma_{\rm{min}}\,S_{*}\,\alpha^{11/19}\,L_{\rm{bol}}\,V_{\rm{sh}}\,t^{-1},$
(24)

from which we can obtain the ratio $\alpha=\epsilon_{e}/\epsilon_{B}$. For
$\Gamma_{\rm{min}}=1$,

$\alpha=\left(0.114\,E_{X}\,\frac{dL_{\rm{X}}}{dE}\,\frac{t}{S_{*}\,V_{\rm{sh}}\,L_{\rm{bol}}}\right)^{19/11},$
(25)

we obtain $\alpha=28.4$, and assuming $\epsilon_{e}=0.1$, we get
$\epsilon_{B}=0.0035$. This result represents a moderate departure from the
energy equipartition assumption compared to other SNe where $\alpha\gtrsim
100,$ such as SN 2011dh (Soderberg et al., 2012; Horesh et al., 2013b), SN
2012aw (Yadav et al., 2014), or SN 2020oi (Horesh et al., 2020). This
deviation from the energy equipartition assumption implies that our radio
analysis overestimated the shock radii by $\sim$19%, and the magnetic field
strengths by a factor of $\sim$2, and underestimated the ratio of the mass-
loss rate to the wind velocity by a factor of $\sim$6.9.

##  7 Discussion

###  7.1 Radio analysis

We presented radio observations of the type IIP supernova SN 2016X obtained
with the VLA. We performed a spectral analysis of these observations using the
SN–CSM interaction model and modelling the emission as self-absorbed
synchrotron radiation. We modelled both individual epochs and multiple epochs
introducing self-similar solutions for the evolution of the radius of the
shock and the magnetic field. The results of the two analyses are compatible
with each other, although the multiple-epoch analysis better constrains the
spectral peaks, which is expected given the limited number of data points for
individual epochs.

The expansion properties of SN 2016X are comparable to other SNe IIP. We
calculated from the multiple-epoch analysis that the radius of the shock
expands from (2.6 $\pm$ 0.1) $\times 10^{15}$ cm to (6.7 $\pm$ 0.7) $\times
10^{15}$ cm during the epochs of our observations, evolving according to a
power law, $R\propto t^{m}$, with index $m=0.76^{+0.09}_{-0.07}$. The
resulting expansion velocities are between 10700 $\pm$ 1100 km s-1 and 8000
$\pm$ 1700 km s-1. The expansion velocities are slightly higher than for other
SNe IIP observed at the same ages (e.g. Chevalier et al. 2006), although this
result depends on the equipartition assumption. The observed evolution of the
magnetic field can be explained by the turbulent motion of particles in the
shocked regions.

Assuming a wind velocity of $w=10$ km s-1 and a self-similar temporal
evolution of the shock, we obtained a mass-loss rate of $\dot{M}=$ (7.8 $\pm$
0.9) $\times\,10^{-7}\,\alpha^{-8/19}\,(\epsilon_{B}/0.1)^{-1}M_{\odot}$ yr-1.
This mass-loss rate is compatible with the expectations of RSGs and comparable
to other SNe IIP, such as SN 1999gi (Schlegel, 2001), SN 2004dj (Nayana et
al., 2018), or SN 2011ja (Chakraborti et al., 2013), if equipartition is
assumed.

We assessed the significance of free-free absorption by estimating the free-
free optical depth at our first epoch of observations. For our lowest
frequency, the free-free optical depth is $\tau_{\rm{ffa}}(\nu=4.8\,\rm{GHz})$
= 0.018. Hence, we conclude that self-absorbed synchrotron is the dominant
absorption mechanism.

Bose et al. (2019) interpreted double-peaked Balmer lines in the SN 2016X
spectrum as a result of bipolar ejecta or asymmetric wind distribution. Our
radio analysis incorporates possible deviations from spherical symmetries
through the volume filling factor, $f$. However, the flat functional
dependence of the physical parameters on the filling factor, $R\propto
f^{-1/19}$ and $B\propto f^{-4/19}$, does not allow for an accurate
determination. For instance, reducing the filling factor by half only changes
the result of the shock radius by less than 4%. Since the CSM density for SN
2016X is proportional to $R^{-2.1}$, in agreement with the density profile of
spherically symmetric wind ($\propto R^{-2}$), we see no evidence of
asymmetric plasma distribution in our radio data. It should also be borne in
mind that the only evidence that supports the claim of asphericity (i.e. the
double-peaked optical lines) only appears when the nebular phase is well
established ($t>200$ days), which may not be relevant for our radio
observations.

###  7.2 Electron cooling and equipartition

The characteristic cooling timescales and frequencies suggest that inverse
Compton electron cooling is important for SN 2016X. The signature of inverse
Compton cooling is evident in the radio spectrum, where we observe a
steepening of the optically thin spectral slope at early times, when electron
cooling is most important. This behaviour has been observed in other SNe of
different types (e.g. SN 2002ap, Björnsson & Fransson 2004; SN 2012aw, Yadav
et al. 2014; SN 2020oi, Horesh et al. 2020).

Combining our results from the radio analysis with X-ray observations reported
by Bose et al. (2019), we find a moderate deviation from energy equipartition
by a factor $\alpha=28.4$. Assuming that the fraction of the energy invested
in particle acceleration is $\epsilon_{e}=0.1$, the fraction of the post-shock
energy deposited in magnetic field intensity is $\epsilon_{B}=0.0035$. This
deviation is modest compared to other SNe, where $\alpha\gtrsim 100$ (e.g. SN
2011dh Soderberg et al. 2012; Horesh et al. 2013b; SN 2012aw Yadav et al.
2014; SN 2020oi Horesh et al. 2020). This deviation from equipartition implies
that our results overestimate the shock radii and magnetic field and
underestimate the mass-loss rate. The corrected mass-loss rate is $\dot{M}=$
(5.4 $\pm$ 0.6) $\times\,10^{-6}\,M_{\odot}$ yr-1, which is similar to that of
SN 1999em (Pooley et al., 2002), SN 2002hh (Nymark et al., 2006), or SN 2004et
(Misra et al., 2007).

Since the free-free optical depth depends on the mass-loss rate, we review our
calculations here. In Figure 4 we present the optical depth as a function of
frequency for our epochs of observations. The maximum value of the free-free
optical depth is $\tau_{\rm{ffa}}(\nu=4.8\,\rm{GHz},t=21.6\,\rm{d})=0.85$,
which is below unity. Additionally, we can check if free-free absorption is
important by verifying whether the shell velocities obtained from the radio
model are lower than the velocities derived from optical spectral lines. For
SN 2016X, we obtained shock velocities of 10700 – 8000 km s-1, while Huang et
al. (2018) inferred velocities in the range $\sim$9000 – 6000 km s-1 between
$t$=20 and $t$=80 days from H$\alpha$ lines. With these results, we argue that
the self-absorbed synchrotron is the dominant absorption mechanism.

Figure 4: Free-free optical depth as a function of frequency for the epochs of
observation of SN 2016X.

###  7.3 SN 2016X compared to other SNe IIP

SN 2016X is one of the few SNe IIP that are detected in both radio and X-rays,
and it is arguably one of the best monitored in radio bands since we collected
36 data points in several frequencies and four different epochs. The other SNe
IIP with similar or better radio coverage are SN 2004dj (Nayana et al., 2018),
SN 2004et (Misra et al., 2007), and SN 2012aw (Yadav et al., 2014).

In this section we compare the properties of SN 2016X to other SNe IIP
observed in radio and X-ray bands. We present an overview of these SNe IIP in
our Figure 5, updating the sample in Nayana et al. (2018)’s Figure 5. SN 2016X
is the most luminous SN IIP in the radio, with $L_{4.8\,\rm{GHz}}\approx
1.5\times 10^{26}$ erg s-1 Hz-1 ($t=43.5$ d). The X-ray luminosity of SN 2016X
is average, larger than for SN 1999gi, SN 1999em, SN 2004dj, or SN 2004et,
similar to SN 2011ja and SN 2002hh, and smaller than for SN 2012aw, SN
2017eaw, and SN 2013ja. While the range of radio luminosities is restricted to
roughly one order of magnitude, the values of X-ray luminosity span four
orders of magnitude.

The mass-loss rate of SN IIP is well constrained within $\sim$1–10
$\times\,10^{-6}\,M_{\odot}$ yr-1 for an assumed wind velocity of 10 km s-1.
SN 2016X exhibits a mass-loss rate towards the lower limit if equipartition is
assumed, but comparable to the highest mass-loss rates taking into account
deviations from equipartition.

The radio spectrum of SN 2016X is also affected by electron cooling via
inverse Compton scattering, such as SN 2004dj, SN 2004et, and SN 2012aw, which
may suggest that it is a common process for SNe IIP. SN 2016X also has in
common with SN 2012aw the fact that both exhibit deviations from energy
equipartition, although the deviation is more severe for SN 2012aw
($\alpha\sim 112$). The mass-loss rate of SN 2016X is higher than that of SN
2012aw by a factor of $\sim$3\. The case of SN 2016X illustrates the
importance of obtaining early radio and X-ray observations of young SNe.

Figure 5: Radio and X-ray properties of SNe IIP. The radio luminosity at
$\sim$5 GHz♮ is measured at the peak of the light curve for SN 1999em, SN
2004dj, SN 2004et, SN 2012aw, and SN 2016X; and limits on the flux peak are
used for the rest of the SNe§. The X-ray luminosity is taken as close to the
peak as possible in the energy ranges $\sim$0.4–8 keV for SN 1999em, SN
2002hh, SN 2004dj, SN 2004et, SN 2013ej, and $\sim$0.3–10 keV for the rest of
the SNe. The mass-loss rate is estimated from radio observations, except for
SN 1999gi, SN 2013ej, and SN 2017eaw. For SN 2016X we plot the mass-loss rate
assuming equipartition (empty symbol) and departure from equipartition (filled
symbol). The vertical lines are not error bars, but rather ranges of possible
values reported by different authors (e.g. Nayana et al. 2018). References:
[1] Pooley et al. 2002, [2] Chevalier et al. 2006, [3] Schlegel 2001, [4]
Pooley & Lewin 2002, [5] Beswick et al. 2005, [6] Beswick et al. 2005, [7]
Chakraborti et al. 2012, [8] Beswick et al. 2004, [9] Misra et al. 2007, [10]
Chakraborti et al. 2013, [11] Chakraborti et al. 2016, [12] Bietenholz et al.
2021, [13] Immler & Brown 2012, [14] Yadav et al. 2014, [15] Argo et al. 2017,
[16] Grefensetette et al. 2017, [17] Szalai et al. 2019, [18] This work.
Notes:  ♮ SN 2013ej was observed in radio at 9 GHz. § SN 1999gi was not
detected in radio, so it is plotted at the arbitrary radio luminosity
$L_{\rm{radio}}=10^{25}$ erg s-1 Hz-1.

##  8 Summary

In this paper we presented radio observations at frequencies 4.8 – 25.0 GHz
obtained at four different epochs at ages in the range 21.6–75.4 days, which
makes it one of the best monitored SN IIP in radio frequencies. The
observations can be well modelled with a self-absorbed synchrotron model,
introducing self-similar expressions for the time evolution of the shock
radius and magnetic field. The results are consistent with the interaction of
ejecta with density profile $\rho_{\rm{SN}}\propto R^{-6.2}$ with a CSM with
density profile $\rho_{\rm{w}}\propto R^{-2.1}$. The expansion of the shock
radius shows moderate deceleration, $R\propto t^{-0.76}$.

We observed steeping of the spectral index in the optically thin regime at
earlier epochs, which we interpreted as the signature of electron cooling via
inverse Compton scattering. We combined our radio data with X-ray
observations, and we found a mild deviation from energy equipartition, with
$\alpha=\epsilon_{e}/\epsilon_{B}\approx 28.4$. This is the second SN IIP to
exhibit this phenomenology, after SN 2012aw.

###### Acknowledgements.

We thank the referee for a thorough assessment of our work and for helping
improve this paper. RRC acknowledges support from the Lady Davis Trust
Fellowship and the Emily Erskine Endowment Funds at the Hebrew University of
Jerusalem. A.H. is grateful for the support by grants from the Israel Science
Foundation, the US-Israel Binational Science Foundation (BSF), and the I-CORE
Program of the Planning and Budgeting Committee and the Israel Science
Foundation. The VLA is a component of the National Radio Astronomy
Observatory, a facility of the National Science Foundation operated under
cooperative agreement by Associated Universities, Inc. We have made extensive
use of the data reduction software casa (McMullin et al., 2007), and the
python packages: astropy (Astropy Collaboration et al., 2013), numpy
(Oliphant, 2006), matplotlib (Hunter, 2007), emcee (Foreman-Mackey et al.,
2013), and corner (Foreman-Mackey, 2016).

## References

  * Argo et al. (2017) Argo, M., Torres, M. P., Beswick, R., & Wrigley, N. 2017, The Astronomer’s Telegram, 10472, 1
  * Arnett (1980) Arnett, W. D. 1980, ApJ, 237, 541
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Beswick et al. (2005) Beswick, R. J., Fenech, D., Thrall, H., et al. 2005, IAU Circ., 8572, 1
  * Beswick et al. (2004) Beswick, R. J., Muxlow, T. W. B., Argo, M. K., Pedlar, A., & Marcaide, J. M. 2004, IAU Circ., 8435, 3
  * Beswick et al. (2005) Beswick, R. J., Muxlow, T. W. B., Argo, M. K., et al. 2005, The Astrophysical Journal, 623, L21
  * Beswick et al. (2005) Beswick, R. J., Muxlow, T. W. B., Argo, M. K., et al. 2005, ApJ, 623, L21
  * Bietenholz et al. (2021) Bietenholz, M. F., Bartel, N., Argo, M., et al. 2021, The Astrophysical Journal, 908, 75
  * Bietenholz et al. (2021) Bietenholz, M. F., Bartel, N., Argo, M., et al. 2021, ApJ, 908, 75
  * Björnsson & Fransson (2004) Björnsson, C.-I. & Fransson, C. 2004, ApJ, 605, 823
  * Bock et al. (2016) Bock, G., Shappee, B. J., Stanek, K. Z., et al. 2016, The Astronomer’s Telegram, 8566, 1
  * Bose et al. (2019) Bose, S., Dong, S., Elias-Rosa, N., et al. 2019, The Astrophysical Journal, 873, L3
  * Bruch et al. (2020) Bruch, R. J., Gal-Yam, A., Schulze, S., et al. 2020, arXiv [2008.09986]
  * Chakraborti et al. (2016) Chakraborti, S., Ray, A., Smith, R., et al. 2016, ApJ, 817, 22
  * Chakraborti et al. (2013) Chakraborti, S., Ray, A., Smith, R., et al. 2013, ApJ, 774, 30
  * Chakraborti et al. (2013) Chakraborti, S., Ray, A., Smith, R., et al. 2013, The Astrophysical Journal, 774, 30
  * Chakraborti et al. (2012) Chakraborti, S., Yadav, N., Ray, A., et al. 2012, ApJ, 761, 100
  * Chevalier (1976) Chevalier, R. A. 1976, ApJ, 207, 872
  * Chevalier (1981) Chevalier, R. A. 1981, The Astrophysical Journal, 251, 259
  * Chevalier (1982) Chevalier, R. A. 1982, The Astrophysical Journal, 259, 302
  * Chevalier (1998) Chevalier, R. A. 1998, The Astrophysical Journal, 499, 810
  * Chevalier & Fransson (2006) Chevalier, R. A. & Fransson, C. 2006, The Astrophysical Journal, 651, 381
  * Chevalier et al. (2006) Chevalier, R. A., Fransson, C., & Nymark, T. K. 2006, The Astrophysical Journal, 641, 1029
  * Corsi et al. (2014) Corsi, A., Ofek, E. O., Gal-Yam, A., et al. 2014, ApJ, 782, 42
  * Foley et al. (2007) Foley, R. J., Smith, N., Ganeshalingam, M., et al. 2007, The Astrophysical Journal, 657, L105
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Fransson et al. (1996) Fransson, C., Lundqvist, P., & Chevalier, R. A. 1996, ApJ, 461, 993
  * Gall et al. (2015) Gall, E. E. E., Polshaw, J., Kotak, R., et al. 2015, Astronomy & Astrophysics, 582, A3
  * Grefensetette et al. (2017) Grefensetette, B., Harrison, F., & Brightman, M. 2017, The Astronomer’s Telegram, 10427, 1
  * Heger et al. (2003) Heger, A., Fryer, C. L., Woosley, S. E., Langer, N., & Hartmann, D. H. 2003, The Astrophysical Journal, 591, 288
  * Hirschi et al. (2004) Hirschi, R., Meynet, G., & Maeder, A. 2004, Astronomy & Astrophysics, 425, 649
  * Horesh et al. (2013a) Horesh, A., Kulkarni, S. R., Corsi, A., et al. 2013a, The Astrophysical Journal, 778, 63
  * Horesh et al. (2012) Horesh, A., Kulkarni, S. R., Fox, D. B., et al. 2012, The Astrophysical Journal, 746, 21
  * Horesh et al. (2020) Horesh, A., Sfaradi, I., Ergon, M., et al. 2020, The Astrophysical Journal, 903, 132
  * Horesh et al. (2013b) Horesh, A., Stockdale, C., Fox, D. B., et al. 2013b, Monthly Notices of the Royal Astronomical Society, 436, 1258
  * Hosseinzadeh et al. (2016) Hosseinzadeh, G., Arcavi, I., McCully, C., et al. 2016, The Astronomer’s Telegram, 8567, 1
  * Huang et al. (2016) Huang, F., Wang, X., Zampieri, L., et al. 2016, The Astrophysical Journal, 832, 139
  * Huang et al. (2015) Huang, F., Wang, X., Zhang, J., et al. 2015, The Astrophysical Journal, 807, 59
  * Huang et al. (2018) Huang, F., Wang, X.-F., Hosseinzadeh, G., et al. 2018, Monthly Notices of the Royal Astronomical Society, 475, 3959
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
  * Immler & Brown (2012) Immler, S. & Brown, P. J. 2012, The Astronomer’s Telegram, 3995, 1
  * Inserra et al. (2011) Inserra, C., Turatto, M., Pastorello, A., et al. 2011, Monthly Notices of the Royal Astronomical Society, 417, 261
  * Kamble et al. (2016) Kamble, A., Margutti, R., Soderberg, A. M., et al. 2016, The Astrophysical Journal, 818, 111
  * Leonard et al. (2002) Leonard, D., Filippenko, A., Gates, E., et al. 2002, Publications of the Astronomical Society of the Pacific, 114, 35
  * Li et al. (2005) Li, W., Van Dyk, S. D., Filippenko, A. V., & Cuillandre, J.-C. 2005, PASP, 117, 121
  * Lundqvist & Fransson (1988) Lundqvist, P. & Fransson, C. 1988, A&A, 192, 221
  * Margutti et al. (2014) Margutti, R., Milisavljevic, D., Soderberg, A. M., et al. 2014, ApJ, 780, 21
  * Margutti et al. (2014) Margutti, R., Milisavljevic, D., Soderberg, A. M., et al. 2014, The Astrophysical Journal, 797, 107
  * Maund & Smartt (2005) Maund, J. R. & Smartt, S. J. 2005, MNRAS, 360, 288
  * McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127
  * Misra et al. (2007) Misra, K., Pooley, D., Chandra, P., et al. 2007, MNRAS, 381, 280
  * Moriya et al. (2014) Moriya, T. J., Maeda, K., Taddia, F., et al. 2014, Monthly Notices of the Royal Astronomical Society, 439, 2917
  * Nayana et al. (2018) Nayana, A. J., Chandra, P., & Ray, A. K. 2018, ApJ, 863, 163
  * Nymark et al. (2006) Nymark, T. K., Fransson, C., & Kozma, C. 2006, A&A, 449, 171
  * Ofek et al. (2014) Ofek, E. O., Sullivan, M., Shaviv, N. J., et al. 2014, The Astrophysical Journal, 789, 104
  * Oliphant (2006) Oliphant, T. E. 2006, A Guide to NumPy, USA
  * Pacholczyk (1970) Pacholczyk, A. G. 1970, Radio astrophysics. Nonthermal processes in galactic and extragalactic sources (San Francisco, California, USA: W. H. Freeman)
  * Pastorello et al. (2007) Pastorello, A., Smartt, S. J., Mattila, S., et al. 2007, Nature, 447, 829
  * Pastorello et al. (2009) Pastorello, A., Valenti, S., Zampieri, L., et al. 2009, Monthly Notices of the Royal Astronomical Society, 394, 2266
  * Pooley & Lewin (2002) Pooley, D. & Lewin, W. H. G. 2002, IAU Circ., 8024, 2
  * Pooley et al. (2002) Pooley, D., Lewin, W. H. G., Fox, D. W., et al. 2002, ApJ, 572, 932
  * Ryder et al. (2004) Ryder, S. D., Sadler, E. M., Subrahmanyan, R., et al. 2004, MNRAS, 349, 1093
  * Schlegel (2001) Schlegel, E. M. 2001, ApJ, 556, L25
  * Schlegel (2001) Schlegel, E. M. 2001, The Astrophysical Journal, 556, L25
  * Smartt (2009) Smartt, S. J. 2009, Annual Review of Astronomy and Astrophysics, 47, 63
  * Smartt (2015) Smartt, S. J. 2015, Publications of the Astronomical Society of Australia, 32, e016
  * Smith (2014) Smith, N. 2014, Annual Review of Astronomy and Astrophysics, 52, 1
  * Smith et al. (2009) Smith, N., Hinkle, K. H., & Ryde, N. 2009, The Astronomical Journal, 137, 3558
  * Smith et al. (2010) Smith, N., Miller, A., Li, W., et al. 2010, The Astronomical Journal, 139, 1451
  * Soderberg et al. (2010) Soderberg, A. M., Brunthaler, A., Nakar, E., Chevalier, R. A., & Bietenholz, M. F. 2010, The Astrophysical Journal, 725, 922
  * Soderberg et al. (2005) Soderberg, A. M., Kulkarni, S. R., Berger, E., et al. 2005, ApJ, 621, 908
  * Soderberg et al. (2012) Soderberg, A. M., Margutti, R., Zauderer, B. A., et al. 2012, ApJ, 752, 78
  * Sorce et al. (2014) Sorce, J. G., Tully, R. B., Courtois, H. M., et al. 2014, Monthly Notices of the Royal Astronomical Society, 444, 527
  * Stockdale et al. (2003) Stockdale, C. J., Weiler, K. W., Van Dyk, S. D., et al. 2003, ApJ, 592, 900
  * Strotjohann et al. (2021) Strotjohann, N. L., Ofek, E. O., Gal-Yam, A., et al. 2021, The Astrophysical Journal, 907, 99
  * Szalai et al. (2019) Szalai, T., Vinkó, J., Könyves-Tóth, R., et al. 2019, ApJ, 876, 19
  * Utrobin & Chugai (2019) Utrobin, V. P. & Chugai, N. N. 2019, Monthly Notices of the Royal Astronomical Society, 490, 2042
  * Van Dyk et al. (2003) Van Dyk, S. D., Li, W., & Filippenko, A. V. 2003, PASP, 115, 1
  * van Loon (2010) van Loon, J. T. 2010, in Astronomical Society of the Pacific Conference Series, Vol. 425, Hot and Cool: Bridging Gaps in Massive Star Evolution, ed. C. Leitherer, P. D. Bennett, P. W. Morris, & J. T. Van Loon, 279
  * Weiler et al. (1991) Weiler, K. W., Dyk, S. D. v., Discenna, J. L., Panagia, N., & Sramek, R. A. 1991, The Astrophysical Journal, 380, 161
  * Weiler et al. (2002) Weiler, K. W., Panagia, N., Montes, M. J., & Sramek, R. A. 2002, Annual Review of Astronomy and Astrophysics, 40, 387
  * Weiler et al. (1990) Weiler, K. W., Panagia, N., & Sramek, R. A. 1990, ApJ, 364, 611
  * Yadav et al. (2014) Yadav, N., Ray, A., Chakraborti, S., et al. 2014, The Astrophysical Journal, 782, 30
  * Zheng & Zhang (2016) Zheng, X.-M. & Zhang, J.-J. 2016, The Astronomer’s Telegram, 8584, 1


# Quantum simulation of dynamical phase transitions in noisy quantum devices

Younes Javanmard<EMAIL_ADDRESS>Ugne Liaubaite
<EMAIL_ADDRESS>Tobias J. Osborne<EMAIL_ADDRESS>hannover.de Luis Santos<EMAIL_ADDRESS>Institut für Theoretische
Physik, Leibniz Universität Hannover, Appelstraße 2, 30167 Hannover, Germany

###### Abstract

Zero-noise extrapolation provides an especially useful error mitigation method
for noisy intermediate-scale quantum devices. Our analysis, based on matrix
product density operators, of the transverse-field Ising model with
depolarizing noise, reveals both advantages and inherent problems associated
with zero-noise extrapolation when simulating non-equilibrium many-body
dynamics. On the one hand, interestingly, noise alters systematically the
behavior of the Loschmidt echo at the dynamical phase transition times,
doubling the number of non-analytic points, and hence inducing an error that,
inherently, cannot be mitigated. On the other, zero-noise extrapolation may be
employed to recover quantum revivals of the Loschmidt echo, which would be
completely missed in the absence of mitigation, and to retrieve faithfully
noise-free inter-site correlations. Our results, which are in good agreement
with those obtained using quantum simulators, reveal the potential of matrix
product density operators for the investigation of the performance of quantum
devices with a large number of qubits and deep noisy quantum circuits.

Quantum computers may outperform their classical counterparts by exploiting
the dimensionality of Hilbert space, which increases exponentially with the
system size Nielsen and Chuang (2010); Feynman (1982); Lidar. and Brun (2013);
Lloyd (1996). They enable a large number of possible applications, including
quantum simulation, optimization, machine learning, and more Biamonte  _et
al._ (2017); Kandala  _et al._ (2017); Nielsen and Chuang (2010). Despite
major technical challenges, recent years have witnessed a rapid experimental
progress on quantum computing using various platforms, including
superconducting qubits Wendin (2017), trapped ion technologies Bruzewicz  _et
al._ (2019); Cirac and Zoller (1995); Blatt and Roos (2012), neutral atoms
Jaksch  _et al._ (1998); Bloch  _et al._ (2012), semiconductor spin qubits
Burkard  _et al._ (2021), and photonics Slussarenko and Pryde (2019);
Weedbrook  _et al._ (2012).

Currently available quantum information processing devices are severely
limited by noise and decoherence. For this reason, they are referred to as
noisy intermediate-scale quantum (NISQ) devices Preskill (2018). These devices
are limited to shallow quantum circuits, but, even then, noise can lead to
faulty estimates. Although quantum error correction can, in principle, reduce
or eliminate the effect of decoherence it is not possible to implement such
schemes on present NISQ devices due to the complex entanglement resource
required.

In response to the challenge of quantum noise, quantum error mitigation
techniques have been developed for NISQ devices Temme  _et al._ (2017); Endo
_et al._ (2018); Li and Benjamin (2017) to ameliorate the deleterious effects
of decoherence in a near-term compatible way. The two major general-purpose
error mitigation methods are probabilistic error cancellation van den Berg
_et al._ (2022) and zero-noise extrapolation (ZNE) Li and Benjamin (2017);
Mari  _et al._ (2021); Lowe  _et al._ (2021). Although there have been several
studies on the general performance of mitigation techniques Takagi  _et al._
(2022a); Tsubouchi  _et al._ (2022); Takagi  _et al._ (2022b), the behavior of
ZNE for large-scale quantum simulations is not well explored.

In this paper, we investigate, by means of matrix product density operators
(MPDO) Verstraete  _et al._ (2004); Zwolak and Vidal (2004); Cheng  _et al._
(2021); Werner  _et al._ (2016); Schollwöck (2011), the use of ZNE for the
study of out-of-equilibrium quantum many-body systems. In particular, we
evaluate the simulation of the dynamics of the transverse-field Ising model by
means of a NISQ device with depolarizing noise, characterized by a noise
channel $\mathcal{E}(\rho)=(1-p)\rho+p\frac{I}{2^{N}}$ Nielsen and Chuang
(2010); Lidar. and Brun (2013), with $\rho$ the density matrix, $p$ the
probabilistic error rate that depends on the device and the circuit, and $N$
the number of qubits. Although there are other possible noise channels, the
depolarizing noise model often appropriately describes the average noise for
large circuits involving many qubits and gates Urbanek  _et al._ (2021). Our
results indicate that, below a given noise and evolution-time threshold, ZNE
may faithfully recover quantum revivals of the Loschmidt echo, which would be
otherwise basically lost, as well as the inter-site correlation dynamics.
However, we show that the noise-free dynamics during the time windows
associated to dynamical quantum phase transitions (DQPTs) is inherently
irretrievable using ZNE, due to the systematic noise-induced doubling of non-
analytic points in the return rate. Our results, which are in good agreement
with direct quantum simulations, show the potential of MPDO calculations for
the study of NISQ devices with many qubits and deep circuits.

### Model.–

We consider in the following the transverse-field Ising model (TFIM):

$\displaystyle\mathcal{H}=-\frac{1}{2}\left[\sum^{N-1}_{i=1}\left(J_{z}\sigma_{i}^{z}\sigma_{i+1}^{z}+J_{x}\sigma_{i}^{x}\sigma_{i+1}^{x}\right)+\sum^{N}_{i=1}h_{x}\sigma_{i}^{x}\right],$
(1)

where $\sigma_{i}^{x,y,z}$ are the spin-$\frac{1}{2}$ Pauli matrices at the
different lattice sites $i=1,\dots,N$, $h_{x}$ is the transverse field, and
$J_{z}$ ($J_{x}$) denotes the $z$ ($x$) spin-spin coupling. For convenience,
we assume open-boundary conditions. For $J_{x}=0$, the model reduces to the
integrable transverse-field Ising chain, which is exactly solvable by a
Jordan-Wigner mapping into non-interacting fermions Pfeuty (1970); Osborne and
Nielsen (2002). To make the model generic, we add a weak $J_{x}$ term, which
becomes equivalent to a two-particle interaction in the fermionic language,
rendering the model non-integrable. In this work, we consider $J_{z}=J$,
$J_{x}=0.1J$, and $h_{x}=0.1J$.

We are particularly interested in the non-equilibrium dynamics of Model (1)
after a quench. Starting with a state $|\psi_{0}\rangle$, the Loschmidt echo
is defined as the probability $\Lambda(t)=|\langle\psi_{0}|e^{-i{\cal
H}t}|\psi_{0}\rangle|^{2}$ to find the time-evolved system in the initial
state. The return-rate function,
$\lambda(t)=-\lim\limits_{N\to\infty}\frac{1}{N}\log\Lambda(t)$, is the real-
time analogous of the free-energy per site. Non-analiticities of $\lambda(t)$
at critical times indicate the presence of DQPTs Heyl  _et al._ (2013); Heyl
(2014, 2018, 2015); Karrasch and Schuricht (2013); Dborin  _et al._ (2022);
Halimeh and Zauner-Stauber (2017); Homrighausen  _et al._ (2017); Lang  _et
al._ (2018); Van Damme  _et al._ (2022a, b), which have been recently realized
experimentally Jurcevic  _et al._ (2017); Wu  _et al._ (2022).

### Master equation.–

In the presence of dissipation, the evolution of the density matrix $\rho$ of
the system is given by the master equation Lidar. and Brun (2013); Nielsen and
Chuang (2010),

$\displaystyle\partial_{t}\rho=$
$\displaystyle\mathcal{L}(t)[\rho]=-i[\mathcal{H},\rho(t)]$
$\displaystyle+\sum_{\begin{subarray}{c}i\\\
\mu=x,y,z\end{subarray}}{\gamma_{i}\left(\sigma^{\mu}_{i}(t)\rho\sigma^{\mu\dagger}_{i}(t)-\frac{1}{2}\\{\sigma^{\mu\dagger}_{i}(t)\sigma^{\mu}_{i}(t),\rho\\}\right)}$
(2)

where the first line provides the unitary part of the time evolution, whereas
the second accounts for dissipation. The Lindblad (jump) operators for
depolarizing noise are given by the Pauli operators $\sigma_{i}^{x,y,z}$, and
encode the system-reservoir couplings, $\gamma_{i}$ characterizes the
dissipation strength at site $i$, and $[]$ and $\\{\\}$ denote commutator and
anti-commutator, respectively. For a time step $\delta t$, the probabilistic
error is related to the noise strength through $p_{\delta
t}=1-e^{-4\gamma\delta t}$ (see the Supplemental Material). Assuming a
constant $\gamma$ for each jump operator, we may re-write Eq. (2) in the form:

$\displaystyle\partial_{t}\rho_{\gamma}(t)={\cal
L}[\rho_{\gamma}]=-i[\mathcal{H}(t),\rho_{\gamma}(t)]+\gamma\mathcal{D}(\rho_{\gamma}(t))$
(3)

### Zero-noise extrapolation.–

ZNE is a relatively simple error-mitigation technique, easy to implement in
post-processing, and which requires neither knowledge of the noise, nor a
gate-level control of the device. The method assumes the ability to scale the
noise parameter to particular values $\alpha\gamma$ with $\alpha>1$. Assuming
a constant noise rate, this may be achieved by a proper stretching of the time
evolution by the factor $\alpha$, while reducing ${\cal H}\to{\cal H}/\alpha$,
resulting in the same overall Hamiltonian evolution, whereas the noise
increases to $\alpha\gamma$. This may be realized in practice by changing the
duration and strength of the pulses that implement the desired operations.
Alternatively one may employ unitary folding, where operations are added that
do not affect the Hamiltonian evolution but increase the noise, e.g. for a
unitary gate $\mathcal{u}$ one may use
$\mathcal{U}\rightarrow\mathcal{U}\mathcal{U}\mathcal{U}^{\dagger}$ Temme  _et
al._ (2017); Kandala  _et al._ (2019); Endo  _et al._ (2021); Giurgica-Tiron
_et al._ (2020); LaRose  _et al._ (2022)(see the Supplemental Materia).

For extrapolation to zero noise value, $\gamma=0$, we employ the Richardson
extrapolation technique Burden and Faires (2010); Giurgica-Tiron  _et al._
(2020). Given a sequence of noise parameters $\alpha_{i}\gamma$ with
$1=\alpha_{0}<\alpha_{1}<...<\alpha_{n}$, the estimated noiseless expectation
value for an operator $M$ is $\langle
M_{E}\rangle(\gamma=0)=\sum_{k=0}^{n}\beta_{k}\langle
M\rangle(\alpha_{k}\gamma)$. The coefficients $\beta_{k}$ satisfy the
relations $\sum^{n}_{k=0}\beta_{k}=1$ and
$\sum^{n}_{l=0}\beta_{l}\alpha_{l}^{k}=0$, giving $\beta_{k}=\prod_{i\neq
k}\frac{\alpha_{i}}{\alpha_{i}-\alpha_{k}}$ Endo  _et al._ (2021); Temme  _et
al._ (2017); Kandala  _et al._ (2019).

### Matrix Product Density Operators.–

We use MPDOs to simulate noisy quantum circuits. MPDOs are an extension of
matrix product states for simulating thermal or dissipative systems Verstraete
_et al._ (2004); Zwolak and Vidal (2004); Schollwöck (2011). We employ an
approach based on the work of Zwolak and Vidal Zwolak and Vidal (2004). Based
on Choi’s isomorphism
$\left|\psi\right\rangle\left\langle\phi\right|=\left|\psi\right\rangle\left|\phi\right\rangle=|{\psi\otimes\phi}\rrangle$,
one can vectorize the density matrix
$\rho=\left|\psi\right\rangle\left\langle\psi\right|\rightarrow|{\rho}\rrangle=|{\psi\otimes\psi}\rrangle$.
Similar to the matrix-product states representation for a quantum state
Montangero (2018); Javanmard  _et al._ (2018); Schollwöck (2011); Cirac  _et
al._ (2021)
$\left|\psi\right\rangle=\sum_{\\{\sigma\\}}{T^{\sigma_{1}}_{1}\ldots
T^{\sigma_{N}}_{N}\left|\sigma_{1}\ldots\sigma_{N}\right\rangle}$, we have in
the MPDO representation
$|{\rho}\rrangle=\sum_{\\{\sigma\sigma^{\prime}\\}}{M^{\sigma_{1}\sigma^{\prime}_{1}}_{1}\ldots
M^{\sigma_{N}\sigma^{\prime}_{N}}_{N}\left|\sigma_{1}\ldots\sigma_{N}\right\rangle\left\langle\sigma^{\prime}_{1}\ldots\sigma^{\prime}_{N}\right|}$,
with $\sigma_{j}$ possible single-site states, and
$M^{\sigma_{i},\sigma^{\prime}_{i}}_{i}=T^{\sigma_{i}}_{i}\otimes(T^{\sigma^{\prime}_{i}}_{i})^{\star}$.
By reformulating the Liouvillian in a vectorized form $\mathcal{L}_{\\#}$, and
for a time-independent $\mathcal{L}_{\\#}$, the time evolution is given by
$|{\rho(T)}\rrangle=e^{T\mathcal{L}_{\\#}}|{\rho(0)}\rrangle$. Here we
consider the case where $\mathcal{L}_{\\#}$ is a sum of nearest-neighbor
interaction terms $\mathcal{L}_{\\#}[\rho]=\sum_{\langle
i,j\rangle}\mathcal{L}^{[i,j]}_{\\#}[\rho]$. This structure allows for using
time-evolving block decimation to simulate the time evolution (see the
Supplemental Material).

In this work, we consider a chain with $N=32$ spins initially prepared in the
state $\left|\psi_{0}\right\rangle=\left|+\right\rangle^{\bigotimes N}$, with
$|+\rangle=H|\downarrow\rangle=\left(|\downarrow\rangle+|\uparrow\rangle\right)/\sqrt{2}$,
with $\\{|\uparrow\rangle,|\downarrow\rangle\\}$ the eigenstates of
$\sigma^{z}$, and $H$ the Hadamard gate. We evolve this state using
Hamiltonian (1). In our MPDO simulation, we use second-order Suzuki-Trotter
decomposition Schollwöck (2011); Hatano and Suzuki (2005) with time step
$\delta t=0.01/J$. We discard states with Schmidt’s values less than $10^{-5}$
during the algorithm run with maximum bond dimension $\chi_{max}=200$.
Although one can reach longer times, we stopped at $t_{max}=14/J$.

### Quantum simulation.–

In addition to our MPDO calculations, we have evaluated the noisy dynamics
using an IBMQ quantum simulator, which employs either the depolarizing noise
model with the one-qubit probabilistic noise parameter $p=0.001$
(corresponding in our MPDO calculations to $\gamma=0.025J$), or a noise model
simulating device noise. The latter was simulated using the qiskit Aer noise
package for the 7-qubit IBMQ Nairobi device IBM-Quantum (2021). The time
evolution given by the unitary operator $\mathcal{U}(t)=e^{-i\mathcal{H}t}$ is
evaluated using Suzuki-Trotter decomposition into $M$ time-steps $\delta
t_{QS}=t/M$, $\mathcal{U}(t)=[\mathcal{U}(\delta t_{QS})]^{M}$, with
$\mathcal{U}(\delta t_{QS})=e^{-i\delta t_{QS}\mathcal{H}_{zz}}e^{-i\delta
t_{QS}\mathcal{H}_{xx}}e^{-i\delta t_{QS}\mathcal{H}_{x}}+\mathcal{O}(\delta
t^{2}_{QS})$ Smith  _et al._ (2019); Fauseweh and Zhu (2021), with ${\cal
H}_{xx}$ (${\cal H}_{zz}$) the $x$ ($z$) spin-spin interactions, and ${\cal
H}_{x}$ the part of the Hamiltonian given by the transversal field. The
Loschmidt amplitude $\langle\psi_{0}|\mathcal{U}(t)|\psi_{0}\rangle=\langle
0^{\otimes N}|H^{\otimes N}\mathcal{U}(t)H^{\otimes N}|0^{\otimes N}\rangle$
can be evaluated using a gate-based quantum computer with the circuit of Fig.
1. We use the mitiq library LaRose  _et al._ (2022) in order to scale the
quantum circuits for the ZNE task.

Figure 1: Quantum circuit employed for the quantum simulation. The gate
parameters are given by $\theta=-2\delta t_{QS}J_{z}$, $\phi=-2\delta
t_{QS}J_{x}$, $\gamma=-2\delta t_{QS}h_{x}$. The part implementing the
evolution operator $\mathcal{U}(\delta t_{QS})$ is repeated $M$ times to reach
the desired time $t$ (see text). Figure 2: Rate function $\lambda(t)$ as a
function of time for different depolarization noise levels, $\gamma$. We show
the results for different re-scaling factors $\alpha$, as well as the
extrapolated value to zero noise. Loschmidt echo $\Lambda(t)=e^{-N\lambda(t)}$
for depolarizing noise for different noise strength $\gamma$. We show the
results for the ideal value and the extrapolated value to zero noise. For
long-enough times and large-enough noise, the Loschmidt echo is bounded by a
fully mixed state, which results in $\lambda=\ln{2}$, indicated by the red
dashed line in each subplot.  Figure 3: Quantum simulation results of the
behavior of the rate function $\lambda(t)$ and the Loschmidt echo
$\Lambda(t)$. (a,d) and (b,e) show the results using the depolarizing noise
model with one-qubit error $p=0.001$, for $N=10$ and $N=6$ qubits,
respectively. (c,f) show the simulation results with the noise model of the
IBMQ Nairobi 7-qubit device. In all the figures we depict the exact
diagonalization (ED) results, the unmitigated evolution, as well as the
mitigated curve. In Figs. (b,e) and (c,f) we show as well the results for
different noise scaling factors $\alpha$. Figure 4: Two-point correlation
$C_{zz}(t)=\langle\sigma_{\frac{N}{2}}^{z}\sigma_{\frac{N}{2}+1}^{z}\rangle$.
(a, b) show MPDO simulations results for depolarizing noise with noise
strength $\gamma$ for $N=32$. (c, d) shows the quantum simulator results for
$N=6$ for $\alpha=1.0,1.5,2.0$. (c) shows the results for depolarizing noise
channel with $p=0.001$ and (d) shows the results using the device noise of the
IBMQ Nairobi. We show the results for the ideal, the unmitigated and the
extrapolated value to zero noise. Figure 5: Error rate
$\varepsilon=\left|\frac{C^{zz}_{\mathrm{ideal}}-C^{zz}_{\mathrm{mitigated}}}{C^{zz,max}_{\mathrm{ideal}}-C^{zz,min}_{\mathrm{ideal}}}\right|$
as a function of noise strength $\gamma$ and time $tJ$. The red curve shows
the threshold $\varepsilon<0.04$, below which the simulation faithfully
retrieves the dynamics of Model (1).

### Loschmidt echo.–

Figure 2 shows the return rate $\lambda(t)$ for different values of $\gamma$
versus the ideal (noise-free) case. For a low noise $\gamma\leq 0.01J$,
extrapolation recovers very well the ideal results, away from the DQPTs, all
the way up to the maximal time considered $t=14/J$. For larger noise, ZNE may
be significantly successful below a given noise and evolution-time threshold.
In the right panels of Fig. 2 we show the Loschmidt echo $\Lambda(t)$. Note
that ZNE recovers the first revival up to $\gamma=0.02J$, although for an
unmitigated noise the revival is basically lost for $\gamma\gtrsim 0.005J$.
For longer times and/or larger $\gamma$, ZNE fails eventually to recover, even
qualitatively, the noise-free results. When accumulating noise, the system
reaches a fully mixed state, $\rho_{s}(t)\equiv\frac{I}{2^{N}}$, such that
$\operatorname{Tr}(\rho_{0}\rho_{s})\equiv\frac{1}{2^{N}}$ and the return rate
is equivalent to $\lambda_{s}\rightarrow\ln{2}$.

The results in the vicinity of the DQPTs are considerably more affected by the
noise. The expected peaks in $\lambda(t)$ associated to the DQPTs are not
recovered even for a very low noise. Interestingly, for any noise, we observe
rather a minimum of $\lambda(t)$, and hence a spurious revival of the
Loschmidt echo, within the time window at which the DQPT occurs in the noise-
free system. This minimum results from the noise-induced doubling of the non-
analytic points of $\lambda(t)$, which seem to indicate noise-induced DQPTs.
As a result, ZNE cannot recover (even qualitatively) the noise-free DQPT. This
points to an inherent limitation of ZNE when applied to the quantum simulation
of many-body dynamics.

Figure 3 shows our results obtained from the quantum simulator. For the case
of depolarizing noise the results resemble those on the MPDO calculations (see
Figs. 3 (a,d) and (b,e)). Note in particular that, as already observed in our
MPDO calculations, using ZNE allows for the basically perfect recovery of the
first revival of the Loschmidt echo, which is almost unobservable for the
unmitigated case. Our quantum simulation results for up to $10$ spins,
although significantly noisy in the DQPT region, suggest as well the doubling
effect mentioned above (see Fig. 3(a)). The MPDO study of systems with
variable number of spins indicate that the doubling, and the associated noise-
induced spurious revival, become more prominent for a larger number of spins.
Hence, finite-size effects may explain the differences with the MPDO results
in the DQPT window.

In contrast, the device-noise simulations, see Figs. 3 (c,f), show that the
noise strength in the considered device is large, and the results approach the
fully-mixed state. As a result, ZNE can faithfully recover the exact results
only for very short times ($t<2/J$), while failing at later times.

### Two-site correlation.–

The analysis of the two-site correlation function in the middle of the spin
chain,
$C^{zz}(t)=\langle\sigma_{\frac{N}{2}}^{z}\sigma_{\frac{N}{2}+1}^{z}\rangle$,
provides further information about the ZNE performance. The zero-noise
evolution is characterized by oscillations of period $\pi/J$, which damp to a
finite value. Figure 4 shows $C^{zz}(t)$ for different values of $\gamma$. We
depict for each case the noise-free result $C_{\mathrm{ideal}}^{zz}(t)$, as
well as the mitigated value $C_{\mathrm{mitigated}}^{zz}(t)$. For short times
(shallow circuits), $Jt<2$, ZNE is successful for all $\gamma$ up to $0.1J$.
For small values of the noise strength, $\gamma\leq 0.02J$, ZNE can recover
the noise-free results up to the longest time considered. Note that,
remarkably, this is so despite of the strong noise-induced damping observed
already for $\gamma=0.005J$. For large $\gamma$ values noise results at long-
enough times in a vanishing correlation, rendering the ZNE eventually unable
to recover the noise-free result.

In Fig. 5, we depict the error rate, defined as
$\varepsilon=\left|\frac{C^{zz}_{\mathrm{ideal}}-C^{zz}_{\mathrm{mitigated}}}{C^{zz,max}_{\mathrm{ideal}}-C^{zz,min}_{\mathrm{ideal}}}\right|$,
as a function of both noise strength $\gamma$ and time $tJ$. Our results show
a threshold, given by the red curve in Fig. 5, such that the ZNE is reliable
below this threshold. In Fig. 5, for a given $\gamma$ we have determined the
threshold time as that in which $\varepsilon<0.04$. One can, however, employ
other thresholds that would provide a similar picture.

Our quantum simulation results for $N=6$, using depolarizing noise, confirm
the significant improvement of the $C_{zz}(t)$ calculation when using ZNE (see
Fig. 4 (c)). In contrast, ZNE produces a more modest improvement for the
device-noise calculation (see Fig. 4 (d)).

### Conclusions.–

We have investigated the performance of zero-noise extrapolation for the
digital quantum simulation of the dynamics of many-body quantum systems. By
means of matrix-product density operators, we have evaluated the particular
case of the simulation of a quenched transverse-field Ising model in the
presence of depolarizing noise. Our results show that zero-noise extrapolation
allows to recover noise-free results, for both the Loschmidt echo and the
inter-site correlation dynamics, below a given noise and evolution-time
threshold. However, the method seems to be inherently unable to recover the
ideal Loschmidt dynamics at dynamical phase transitions due to the appearance
of a noise-induced revival. Our matrix-product density analysis, backed up by
direct quantum simulations, provides a blueprint for the study of the
performance of ZNE on noisy quantum processors with large number of qubits and
deep circuits. Future work will address different noise models using the MPDO
simulations for different quantum computing tasks, such as quantum
optimization, and will benchmark these results with those on the actual
device, developing a device-specific noise analysis. Another possibility would
be to consider the role of cross-talks between qubits in the simulations and
show how to mitigate them. Furthermore, one could expand such simulations to
higher dimensions and a more realistic quantum chip geometry.

###### Acknowledgements.

### Acknowledgments.–

We acknowledge funding by the Volkswagen foundation and the Ministry of
Science and Culture of Lower Saxony through _Quantum Valley Lower Saxony Q1
(QVLS-Q1)_ , and by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) – under Germany’s Excellence Strategy – EXC-2123 Quantum-Frontiers
– 390837967.

## References

  * Nielsen and Chuang (2010) M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information: 10th Anniversary Edition_ (Cambridge University Press, 2010).
  * Feynman (1982) R. P. Feynman, Int J Theor Phys 21, 467 (1982).
  * Lidar. and Brun (2013) D. A. Lidar. and T. A. Brun,  _Quantum Error Correction_ (Cambridge University Press, 2013).
  * Lloyd (1996) S. Lloyd, Science 273, 1073 (1996).
  * Biamonte  _et al._ (2017) J. Biamonte, P. Wittek, N. Pancotti, P. Rebentrost, N. Wiebe, and S. Lloyd, Nature 549, 195 (2017).
  * Kandala  _et al._ (2017) A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Nature 549, 242 (2017).
  * Wendin (2017) G. Wendin, Rep. Prog. Phys. 80, 106001 (2017).
  * Bruzewicz  _et al._ (2019) C. D. Bruzewicz, J. Chiaverini, R. McConnell, and J. M. Sage, Applied Physics Reviews 6, 021314 (2019).
  * Cirac and Zoller (1995) J. I. Cirac and P. Zoller, Phys. Rev. Lett. 74, 4091 (1995).
  * Blatt and Roos (2012) R. Blatt and C. F. Roos, Nature Phys 8, 277 (2012).
  * Jaksch  _et al._ (1998) D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Phys. Rev. Lett. 81, 3108 (1998).
  * Bloch  _et al._ (2012) I. Bloch, J. Dalibard, and S. Nascimbène, Nature Phys 8, 267 (2012).
  * Burkard  _et al._ (2021) G. Burkard, T. D. Ladd, J. M. Nichol, A. Pan, and J. R. Petta,  (2021), 10.48550/arXiv.2112.08863, comment: Rev. Mod. Phys. - Corrections and comments welcome, arXiv:2112.08863 [cond-mat, physics:physics, physics:quant-ph] .
  * Slussarenko and Pryde (2019) S. Slussarenko and G. J. Pryde,  (2019), 10.1063/1.5115814, comment: 21 pages, 3 figures. A brief review on some topics in photonic quantum computing with lots of references to longer specialist reviews. Close to published version, arXiv:1907.06331 [quant-ph] .
  * Weedbrook  _et al._ (2012) C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, Rev. Mod. Phys. 84, 621 (2012).
  * Preskill (2018) J. Preskill, Quantum 2, 79 (2018).
  * Temme  _et al._ (2017) K. Temme, S. Bravyi, and J. M. Gambetta, Phys. Rev. Lett. 119, 180509 (2017).
  * Endo  _et al._ (2018) S. Endo, S. C. Benjamin, and Y. Li, Phys. Rev. X 8, 031027 (2018).
  * Li and Benjamin (2017) Y. Li and S. C. Benjamin, Phys. Rev. X 7, 021050 (2017).
  * van den Berg  _et al._ (2022) E. van den Berg, Z. K. Minev, A. Kandala, and K. Temme,  (2022), 10.48550/arXiv.2201.09866, arXiv:2201.09866 [quant-ph] .
  * Mari  _et al._ (2021) A. Mari, N. Shammah, and W. J. Zeng, Phys. Rev. A 104, 052607 (2021).
  * Lowe  _et al._ (2021) A. Lowe, M. H. Gordon, P. Czarnik, A. Arrasmith, P. J. Coles, and L. Cincio, Phys. Rev. Research 3, 033098 (2021).
  * Takagi  _et al._ (2022a) R. Takagi, S. Endo, S. Minagawa, and M. Gu,  (2022a), 10.48550/arXiv.2109.04457, arXiv:2109.04457 [quant-ph] .
  * Tsubouchi  _et al._ (2022) K. Tsubouchi, T. Sagawa, and N. Yoshioka,  (2022), 10.48550/arXiv.2208.09385, arXiv:2208.09385 [quant-ph] .
  * Takagi  _et al._ (2022b) R. Takagi, H. Tajima, and M. Gu,  (2022b), 10.48550/arXiv.2208.09178, arXiv:2208.09178 [quant-ph] .
  * Verstraete  _et al._ (2004) F. Verstraete, J. J. García-Ripoll, and J. I. Cirac, Phys. Rev. Lett. 93, 207204 (2004).
  * Zwolak and Vidal (2004) M. Zwolak and G. Vidal, Phys. Rev. Lett. 93, 207205 (2004).
  * Cheng  _et al._ (2021) S. Cheng, C. Cao, C. Zhang, Y. Liu, S.-Y. Hou, P. Xu, and B. Zeng, Phys. Rev. Research 3, 023005 (2021).
  * Werner  _et al._ (2016) A. H. Werner, D. Jaschke, P. Silvi, M. Kliesch, T. Calarco, J. Eisert, and S. Montangero, Phys. Rev. Lett. 116, 237201 (2016).
  * Schollwöck (2011) U. Schollwöck, Annals of Physics 326, 96 (2011), january 2011 Special Issue.
  * Urbanek  _et al._ (2021) M. Urbanek, B. Nachman, V. R. Pascuzzi, A. He, C. W. Bauer, and W. A. de Jong, Phys. Rev. Lett. 127, 270502 (2021).
  * Pfeuty (1970) P. Pfeuty, Annals of Physics 57, 79 (1970).
  * Osborne and Nielsen (2002) T. J. Osborne and M. A. Nielsen, Phys. Rev. A 66, 032110 (2002).
  * Heyl  _et al._ (2013) M. Heyl, A. Polkovnikov, and S. Kehrein, Phys. Rev. Lett. 110, 135704 (2013).
  * Heyl (2014) M. Heyl, Phys. Rev. Lett. 113, 205701 (2014).
  * Heyl (2018) M. Heyl, Rep. Prog. Phys. 81, 054001 (2018).
  * Heyl (2015) M. Heyl, Phys. Rev. Lett. 115, 140602 (2015).
  * Karrasch and Schuricht (2013) C. Karrasch and D. Schuricht, Phys. Rev. B 87, 195104 (2013).
  * Dborin  _et al._ (2022) J. Dborin, V. Wimalaweera, F. Barratt, E. Ostby, T. E. O’Brien, and A. G. Green, Nat Commun 13, 5977 (2022).
  * Halimeh and Zauner-Stauber (2017) J. C. Halimeh and V. Zauner-Stauber, Phys. Rev. B 96, 134427 (2017).
  * Homrighausen  _et al._ (2017) I. Homrighausen, N. O. Abeling, V. Zauner-Stauber, and J. C. Halimeh, Phys. Rev. B 96, 104436 (2017).
  * Lang  _et al._ (2018) J. Lang, B. Frank, and J. C. Halimeh, Phys. Rev. B 97, 174401 (2018).
  * Van Damme  _et al._ (2022a) M. Van Damme, J.-Y. Desaules, Z. Papić, and J. C. Halimeh,  (2022a), 10.48550/arXiv.2210.02453, arXiv:2210.02453 [cond-mat, physics:quant-ph] .
  * Van Damme  _et al._ (2022b) M. Van Damme, T. V. Zache, D. Banerjee, P. Hauke, and J. C. Halimeh,  (2022b), 10.48550/arXiv.2203.01337, arXiv:2203.01337 [cond-mat, physics:hep-lat, physics:quant-ph] .
  * Jurcevic  _et al._ (2017) P. Jurcevic, H. Shen, P. Hauke, C. Maier, T. Brydges, C. Hempel, B. P. Lanyon, M. Heyl, R. Blatt, and C. F. Roos, Phys. Rev. Lett. 119, 080501 (2017).
  * Wu  _et al._ (2022) L.-N. Wu, J. Nettersheim, J. Feß, A. Schnell, S. Burgardt, S. Hiebel, D. Adam, A. Eckardt, and A. Widera,  (2022), 10.48550/arXiv.2208.05164, arXiv:2208.05164 [cond-mat, physics:quant-ph] .
  * Kandala  _et al._ (2019) A. Kandala, K. Temme, A. D. Córcoles, A. Mezzacapo, J. M. Chow, and J. M. Gambetta, Nature 567, 491 (2019).
  * Endo  _et al._ (2021) S. Endo, Z. Cai, S. C. Benjamin, and X. Yuan, J. Phys. Soc. Jpn. 90, 032001 (2021).
  * Giurgica-Tiron  _et al._ (2020) T. Giurgica-Tiron, Y. Hindy, R. LaRose, A. Mari, and W. J. Zeng, in  _2020 IEEE International Conference on Quantum Computing and Engineering (QCE)_ (2020) pp. 306–316.
  * LaRose  _et al._ (2022) R. LaRose, A. Mari, S. Kaiser, P. J. Karalekas, A. A. Alves, P. Czarnik, M. E. Mandouh, M. H. Gordon, Y. Hindy, A. Robertson, P. Thakre, M. Wahl, D. Samuel, R. Mistri, M. Tremblay, N. Gardner, N. T. Stemen, N. Shammah, and W. J. Zeng, Quantum 6, 774 (2022).
  * Burden and Faires (2010) R. L. Burden and J. D. Faires,  _Numerical Analysis_ (Cengage Learning, 2010).
  * Montangero (2018) S. Montangero,  _Introduction to Tensor Network Methods: Numerical Simulations of Low-Dimensional Many-Body Quantum Systems_ (Springer International Publishing, Cham, 2018).
  * Javanmard  _et al._ (2018) Y. Javanmard, D. Trapin, S. Bera, J. H. Bardarson, and M. Heyl, New J. Phys. 20, 083032 (2018).
  * Cirac  _et al._ (2021) J. I. Cirac, D. Pérez-García, N. Schuch, and F. Verstraete, Rev. Mod. Phys. 93, 045003 (2021).
  * Hatano and Suzuki (2005) N. Hatano and M. Suzuki,  _Quantum Annealing and Other Optimization Methods_ ,  Lecture Notes in Physics, 37 (2005).
  * IBM-Quantum (2021) IBM-Quantum, “Qiskit: An open-source framework for quantum computing,”  (2021).
  * Smith  _et al._ (2019) A. Smith, M. S. Kim, F. Pollmann, and J. Knolle, npj Quantum Information 5, 106 (2019).
  * Fauseweh and Zhu (2021) B. Fauseweh and J.-X. Zhu, Quantum Inf Process 20, 138 (2021).

## Supplemental Material: Quantum simulation of dynamical phase transitions in
noisy quantum devices

We provide additional details and numerical results supplementing the
conclusions from the main text.

##  I Matrix product density operators

The vectorized Liouvillian is given by

$\displaystyle\mathcal{L}_{\\#}$ $\displaystyle\equiv-i(\mathcal{H}\otimes
I-I\otimes\mathcal{H}^{T})$ $\displaystyle+\sum_{\mu}{(L_{\mu}\otimes
L^{*}_{\mu}-\frac{1}{2}L^{\dagger}_{\mu}L_{\mu}\otimes I-\frac{1}{2}I\otimes
L^{*}_{\mu}L^{T}_{\mu})}$ (S1)

Fig. S1 shows different steps for MPDO time evolution

Figure S1: MPDO steps to evolve a system based on the Liouvilian master
equation. Figure S2: Discarded error of the MPDO simulation for the main
results in the paper.

##  II Depolarizing noise channel

We comment in the following on the relation between the probabilistic error
$p$ of the depolarizing channel, and the dissipation strength $\gamma$ in the
Lindblad master equation:

$\displaystyle\dot{\rho}=\gamma\sum_{\mu=x,y,z}\left(-\rho+\sigma^{\mu}\rho\sigma^{\mu}\right)=2\gamma(I-2\rho).$
(S2)

Here we have employed the identity
$\sum_{\mu}\sigma^{\mu}M\sigma^{\mu}=2\operatorname{Tr}(M)I-M$, valid for any
$2\times 2$ matrix $M$, and $\operatorname{Tr}(\rho)=1$. Using the
parametrization: $\rho(t)=\frac{1}{2}(I+\vec{r}(t)\cdot\vec{\sigma})$, we can
rewrite the master equation in the form:

$\displaystyle\dfrac{1}{2}\sum_{\mu}\dot{r}_{\mu}\sigma^{\mu}=-2\gamma\sum_{\mu}r_{\mu}\sigma^{\mu},$
(S3)

and hence $r_{\mu}(t)=e^{-4\gamma t}r_{\mu}(0)$. Then:

$\displaystyle\rho(\delta t)=\frac{1}{2}(I+e^{-4\gamma\delta
t}\vec{r}(0)\cdot\vec{\sigma})=\left(1-p(\delta
t)\right)\rho(0)+\frac{I}{2}p(\delta t),$ (S4)

which has the form of a depolarizing channel with $p(\delta
t)=1-e^{-4\gamma\delta t}$.

##  III Benchmarking the MPDO simulations vs exact diagonalization

In order to benchmark our MPDO results, we compare in the following the
results obtained using exact diagonalisation and those from the MPDO
calculation. We consider a small system site ($N=8$) and different noise
strengths. We consider the integrable TFIM:

$\displaystyle\mathcal{H}=\sum_{i}{(\sigma_{i}^{z}\sigma_{i+1}^{z}+h\sigma_{i}^{z})},$
(S5)

with $J=1$ and $h=0.1$, and we quench from the same initial state as in the
main text. Our results for the return rate depicted in Fig. S3 show an
excellent agreement between exact diagonalization and MPDO calculations. Note
also that compared to our MPDO calculations with $N=32$ sites, we do not
observe for $N=8$ the same behavior within the DQPT region, which rather
resembles more closely our quantum simulation results. This indicates that the
differences in the main text between the quantum simulation and MPDO results
in the vicinity of the DQPT are due to the different system sizes.

Figure S3: Dynamics of model (S5) for $N=8$ sites. Comparison of the return
rate $\lambda(t)$ evaluated with exact diagonalization (ED) and with MPDO for
different noise strengths. In our MPDO calculations we have employed a maximum
bond dimension $\chi_{max}=200$.

##  IV Noise scaling in the quantum simulation

Noise scaling in gate-based quantum computation is based on unitary folding,
which consists in replacing a unitary operation $U$ by $U\rightarrow
U(U^{\dagger}U)^{n}$ with $n$ being a positive integer. As $U^{\dagger}U=I$,
folding does not have any logical effect. However, the noise is increased with
the prolongation of the circuit. The mitiq library is used to scale the
quantum circuits for zero noise extrapolation. The scaling is done using
folding at random introduced in LaRose  _et al._ (2022); Giurgica-Tiron  _et
al._ (2020). Before scaling the circuits, one needs to transpile the quantum
circuits using a set of universal basis gates, so that the depth of the
circuit ran is precisely known; Here we have used basis gates
$=[u1,~{}u2,~{}u3,~{}cx]$.

The value of the scaling factor $\alpha$ is important for the performance of
ZNE. There are two options for unitary folding, circuit-folding and layer- or
gate-based folding. The methods differ in the values of scaling factor
$\alpha$ that can be reached. In circuit-folding the entire circuit is folded,
and the values of $\alpha$ are of form $\alpha=1+2n$ where $n$ is the number
of identities added. In order to get finer scaling, i.e. $\alpha$ being able
to reach all real numbers, one can scale part of a circuit (folding only a
selected number of layers from the entire circuit). The values of $\alpha$ are
then of form $\alpha=1+(2k)/d$ where $k$ is any natural number and $d$ is the
depth of the circuit. In this work we have employed layer-based folding. For
more detailed information see Ref. Giurgica-Tiron  _et al._ (2020).