# Uncovering the Ghostly Remains of an Extremely Diffuse Satellite in the
Remote Halo of NGC 253 111This research is based on data collected at the
Subaru Telescope, which is operated by the National Astronomical Observatory
of Japan.

Sakurako Okamoto Subaru Telescope, National Astronomical Observatory of Japan,
650 North A’ohoku Place, Hilo, HI 96720, U.S.A.  National Astronomical
Observatory of Japan, Osawa 2-21-1, Mitaka, Tokyo, 181-8588, JAPAN  The
Graduate University for Advanced Studies, SOKENDAI, Osawa 2-21-1, Mitaka,
Tokyo 181-8588, JAPAN  Annette M.N. Ferguson Institute for Astronomy,
University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh, EH9 3HJ
U.K.  Nobuo Arimoto  National Astronomical Observatory of Japan, Osawa 2-21-1,
Mitaka, Tokyo, 181-8588, JAPAN  The Graduate University for Advanced Studies,
SOKENDAI, Osawa 2-21-1, Mitaka, Tokyo 181-8588, JAPAN  Itsuki Ogami The
Graduate University for Advanced Studies, SOKENDAI, Osawa 2-21-1, Mitaka,
Tokyo 181-8588, JAPAN  Rokas Žemaitis Institute for Astronomy, University of
Edinburgh, Royal Observatory, Blackford Hill, Edinburgh, EH9 3HJ U.K.  Masashi
Chiba Astronomical Institute, Tohoku University, 6-3 Aoba, Aramaki, Aoba-ku,
Sendai, Miyagi 980-8578, JAPAN  Mike J. Irwin Institute of Astronomy,
University of Cambridge, Madingley Road, Cambridge CB3 0HA, U.K.  In Sung Jang
Department of Astronomy & Astrophysics, University of Chicago, 5640 South
Ellis Avenue, Chicago, IL 60637, U.S.A.  Jin Koda Department of Physics and
Astronomy, Stony Brook University, Stony Brook, NY 11794-3800, U.S.A.  Yutaka
Komiyama Dept. of Advanced Sciences, Faculty of Science and Engineering, Hosei
University, 3-7-2 Kajino-cho, Koganei-shi, Tokyo 184-8584, JAPAN  Myung Gyoon
Lee Astronomy Program, Department of Physics and Astronomy, SNUARC, Seoul
National University, 1 Gwanak-ro, Gwanak-gu, Seoul 08826, Republic of Korea
Jeong Hwan Lee The Center for High Energy Physics, Kyungpook National
University, 80 Daehakro, Bukgu Daegu 41566, Republic of Korea  Michael Rich
Department of Physics and Astronomy, The University of California, Los
Angeles, CA 90095, U.S.A.  Masayuki Tanaka National Astronomical Observatory
of Japan, Osawa 2-21-1, Mitaka, Tokyo, 181-8588, JAPAN  The Graduate
University for Advanced Studies, SOKENDAI, Osawa 2-21-1, Mitaka, Tokyo
181-8588, JAPAN  Mikito Tanaka Dept. of Advanced Sciences, Faculty of Science
and Engineering, Hosei University, 3-7-2 Kajino-cho, Koganei-shi, Tokyo
184-8584, JAPAN

###### Abstract

We present the discovery of NGC253-SNFC-dw1, a new satellite galaxy in the
remote stellar halo of the Sculptor Group spiral, NGC 253. The system was
revealed using deep resolved star photometry obtained as part of the Subaru
Near-Field Cosmology Survey that uses the Hyper Suprime-Cam on the Subaru
Telescope. Although rather luminous ($\rm{M_{V}}=-11.7\pm 0.2$) and massive
($M_{*}\sim 1.25\times 10^{7}~{}\rm{M}_{\odot}$), the system is one of the
most diffuse satellites yet known, with a half-light radius of
$\rm{R_{h}}=3.37\pm 0.36$ kpc and an average surface brightness of $\sim 30.1$
mag arcmin-2 within the $\rm{R_{h}}$. The colour-magnitude diagram shows a
dominant old ($\sim 10$ Gyr) and metal-poor ($\rm{[M/H]}=-1.5\pm 0.1$ dex)
stellar population, as well as several candidate thermally-pulsing asymptotic
giant branch stars. The distribution of red giant branch stars is asymmetrical
and displays two elongated tidal extensions pointing towards NGC 253,
suggestive of a highly disrupted system being observed at apocenter.
NGC253-SNFC-dw1 has a size comparable to that of the puzzling Local Group
dwarfs Andromeda XIX and Antlia 2 but is two magnitudes brighter. While
unambiguous evidence of tidal disruption in these systems has not yet been
demonstrated, the morphology of NGC253-SNFC-dw1 clearly shows that this is a
natural path to produce such diffuse and extended galaxies. The surprising
discovery of this system in a previously well-searched region of the sky
emphasizes the importance of surface brightness limiting depth in satellite
searches.

Dwarf galaxy (416) — Galaxies (573) — Stellar Populations (1622) — Survey
(1672) — Photometry (1234)

††journal: ApJL††facilities: Subaru(HSC)

##  1 Introduction

The $\Lambda$CDM cosmological model predicts that galaxies form
hierarchically; large galaxies like the Milky Way originate from small over-
densities in the primordial matter distribution and grow via the agglomeration
of numerous smaller building blocks, some of which survive later merging and
represent the present-day dwarf satellites (e.g. Bullock & Johnston, 2005;
Pillepich et al., 2014). The numbers and properties of satellite galaxies
around massive hosts has emerged as a powerful way to test the $\Lambda$CDM
paradigm on small scales as well as the physics of galaxy formation (e.g.
Bullock & Boylan-Kolchin, 2017; Sales et al., 2022). As a result, there has
been much interest in characterising the satellite system of the Milky Way and
its nearest neighbours (e.g. Homma et al., 2019; Carlsten et al., 2022;
Doliva-Dolinsky et al., 2023). Although the sample size beyond the Local Group
(LG) is still limited, such observations are critical to benchmark results
from the Milky Way and M31. Recent wide-field imaging and/or survey campaigns
are having remarkable success in uncovering new satellites around Local Volume
galaxies (e.g. Chiboucas et al., 2013; Crnojević et al., 2016; Okamoto et al.,
2019; Mutlu-Pakdil et al., 2022) but these studies rapidly lose sensitivity at
surface brightnesses $\mu_{V}\gtrsim 28-29$ mag arcsec-2, raising questions
about what they might be missing.

NGC 253, lying at D $\sim 3.5$ Mpc (Radburn-Smith et al., 2011), is a barred
spiral galaxy located in the Sculptor group. It has a similar luminosity and
virial radius to the Milky Way and M31 (Karachentsev et al., 2013; Mutlu-
Pakdil et al., 2021), making it a prime object for comparative studies with
the LG spirals. An extended and structured stellar halo has long been known in
this galaxy, detected in both integrated light (e.g. Beck et al., 1982; Malin
& Hadley, 1997) and resolved star counts (Davidge, 2010; Bailin et al., 2011;
Greggio et al., 2014; Harmsen et al., 2017). It also hosts a rich population
of halo globular clusters (Cantiello et al., 2018).

In recent years, attention has focused on the faint satellite system of NGC
253. The PISCeS survey used Magellan/Megacam resolved star photometry to
discover five new dwarf satellites of NGC 253 (Sand et al., 2014; Toloba et
al., 2016; Mutlu-Pakdil et al., 2022), with two of them, Scl-MM-dw2 (hereafter
MM-dw2) and Scl-MM-dw3, independently discovered as NGC253-dw2 by Romanowsky
et al. (2016) and Donatiello II by Martínez-Delgado et al. (2021),
respectively. Martínez-Delgado et al. (2021) proposed two further dwarf galaxy
candidates from Dark Energy Survey data and Carlsten et al. (2022) suggested
two others from DECaLS data, both from integrated light analyses. Mutlu-Pakdil
et al. (2024) confirm three of these latter systems with Hubble Space
Telescope photometry and undertake a comprehensive study of the NGC 253
satellite system as a whole, finding it to be deficient in luminous satellites
compared systems of similar stellar mass.

In this Letter, we report the discovery of a highly-disrupted satellite galaxy
in the outer halo of NGC 253. Although luminous, the system is extremely
diffuse and thus has been missed by all previous searches in this part of the
sky. The discovery represents the first result from the “Subaru Near Field
Cosmology” survey (hereafter SNFC), a systematic imaging survey of a sample of
late-type galaxies in the Local Volume in order to study their resolved
stellar halos. In Section 2, we summarize the observations and data reduction
methodology. We derive the physical properties of the dwarf, which we name NGC
253-SNFC-dw1 (hereafter SNFC-dw1), in Section 3 and discuss and conclude in
Section 4.

##  2 Observation and data

Figure 1: Left: The HSC pointings around NGC 253 overlaid on the Digitized Sky
Survey image using the Aladin Sky Atlas(Bonnarel et al., 2000). The orange
tile indicates the pointing where SNFC-dw1 is discovered. The location of
SNFC-dw1 and the previously-known satellites within the footprint are
indicated. Right: The spatial distribution of RGB stars selected using the
polygon in Figure 2(a). The color of each point represents the $(g-i)_{0}$
color of the star with transparency. The dotted ellipses indicate the
$2\times\rm{R_{25}}$ and $4\times\rm{R_{25}}$ radii of NGC 253 (de Vaucouleurs
et al., 1991), assuming b/a=0.7. The solid and dashed circles show the
circular radius of 9 arcmin centered on SNFC-dw1 and a nearby reference field.
The $2.5\times\rm{R_{h}}$ of the two known dwarf satellites within the FOV are
also shown as solid lines and labeled. Lower middle: The smoothed density map
of metal-poor RGB stars centered on SNFC-dw1 and selected within the left part
of the solid polygon of Figure 2(b). The magenta cross shows the estimated
centroid of the structure, of which the coordinate is shown in Table1. The
dashed line shows the fitted $\rm{R_{h}}$ of the Sérsic profile. The black
points represent the TP-AGB candidates. Gray-shaded areas are the masked
regions affected by saturated foreground stars. The white arrow indicates the
direction toward the NGC 253 center.

We observed three pointings around NGC 253 using the Hyper Suprime-Cam (HSC)
on the Subaru 8.2m telescope (left panel of Figure 1). The HSC comprises 104
CCD detectors and provides a field-of-view (FOV) of 1.76 $\rm{deg^{2}}$ with a
pixel scale of 0.17″(Miyazaki et al., 2018). The observations were obtained in
queue-mode as part of the open-use SNFC intensive program during the period
2019-2022. For NGC 253, $g$\- and $i$-band images were acquired under seeing
conditions ranging from 0.7-1.0″. For the two northern fields, we obtained a
total exposure time of 12500 sec and 2750 sec in the $g$\- and $i$-bands,
respectively. In the southern field, we obtained the same exposure time in the
$i$-band but the $g$-band exposures amounted to only 6000 sec; furthermore,
the seeing was slightly worse in this field than in the northern fields. Our
exposure times are designed to provide sensitivity to a broad range of stellar
metallicities, including the most metal-rich halo stars.

The data were processed using the reduction pipeline hscPipe 8.4 (Bosch et
al., 2018), which is based on a software suite being developed for the Vera C.
Rubin Observatory data (Axelrod et al., 2010; Jurić et al., 2017; Ivezić et
al., 2019). Each frame was corrected for bias, dark, and flat-fielding, and
then the sky was modeled internally and subtracted with an aggressive 32-pixel
mesh. The frames in each filter were calibrated using Pan-STARRS1 catalog
(Magnier et al., 2013; Chambers et al., 2016), then mosaicked and coadded. The
final photometry is in the HSC filter system and in AB magnitudes. Forced
photometry based on the point spread function (PSF), CModel, and Kron model
was performed on the coadded images, with the $g$-band image providing the
initial detections. Extended sources were excluded by considering the PSF to
CModel flux ratio, following Pucha et al. (2019). In particular, we classified
an object as a point-source if the flux ratio is within $1\sigma$ of unity in
the $i$-band. The point-source detection limits of the coadded images were
$g\sim 27.9$ and $i\sim 26.3$ with SNR$=5$ for the northern two fields and
$g\sim 27.5$ and $i\sim 25.8$ with SNR$=5$ for the southern field. These
depths correspond to roughly 1.5 to 2 mag below the tip of red giant branch
(TRGB) at the distance of NGC 253 (Radburn-Smith et al., 2011).

The Galactic extinction corrections were taken from the Schlegel et al. (1998)
reddening map using the coefficients from Schlafly & Finkbeiner (2011). We
applied an extinction correction to each source individually according to its
location, assuming a Fitzpatrick (1999) reddening law with $\rm{R_{V}}=3.1$.
The mean reddening across the field is rather small, amounting to $E(B-V)\sim
0.0161$.

Artificial star tests were performed on the coadded images to estimate the
accuracy and completeness of the photometric catalog. Artificial point sources
were created using the PSF models generated by the pipeline and injected into
the coadded images as described in Ogami et al. (2024). The resulting images
were processed in the same way as the original ones. About 30,000 artificial
stars were added to the images in intervals of 1.0 mag from 20 mag to 24 mag,
and 0.5 mag from 24 mag to 28 mag, in both the $g$\- and $i$-bands. We
considered the artificial stars to be recovered if the PSF to CModel flux
ratio of the detected source was within $1\sigma$ of unity at the input
magnitude. The completeness was examined in elliptical annuli centered on NGC
253, assuming an axis ratio $\rm{b}/\rm{a}=0.7$ which is judged by eye to be a
good match to the shape of the outer halo. At the radial distances of
relevance for this paper (30-110 arcmin), the recovery fraction is almost
constant for a given input magnitude. We interpolated the recovery fractions
to produce a 2D histogram representing the completeness fraction as function
of magnitude and color, and used this to correct for incompleteness where
required.

Figure 2: De-reddened Hess diagrams of stellar objects in NGC 253, the SNFC-
dw1 region, a nearby reference field, and the field-subtracted SNFC-dw1
region. The bin sizes are 0.15 mag in the x-axis and 0.1 mag in the y-axis.
(a) Stars within an elliptical annulus of $2\times\rm{R_{25}}<$ r
$<4\times\rm{R_{25}}$ centered on NGC 253, assuming b/a$=0.7$. The solid
polygon indicates the area used to select RGB stars for the right panel of
Figure 1. PARSEC isochrones of age 10 Gyr and
$\rm{[M/H]}=-2.2,-1.8,-1.4,-1.0,-0.7,-0.4$ are shifted to NGC 253’s distance
and overlaid as solid lines (Radburn-Smith et al., 2011; Bressan et al.,
2012). (b) Stars within a circular radius of 9 arcmin of SNFC-dw1. A 10 Gyr
isochrone of $\rm{[M/H]}=-1.8$ is shifted to the estimated SNFC-dw1 distance
and overlaid. The dashed line in the solid polygon is used to separate the
metal-poor RGB stars. (c) Stellar objects in a reference field, selected as an
equal area region that lies at the same elliptical radius from NGC 253 as
SNFC-dw1. (d) The field-subtracted Hess diagram of SNFC-dw1, created by
subtracting panel (c) from (b). The gray color indicates negative values.

##  3 Properties of NGC253-SNFC-dw1

###  3.1 Star Counts and Stellar Populations

The right panel of Figure 1 shows the spatial distribution of RGB stars across
three HSC pointings around NGC 253. RGB stars are selected from stellar
objects within the black solid polygon on the color-magnitude diagram (CMD)
shown in Figure 2(a), the boundaries of which have been chosen to optimize the
selection of RGB stars at the distance of NGC 253, while limiting the number
of foreground and background contaminants. The color of each point in the map
denotes the $(g-i)_{0}$ color of the star, which can be interpreted as a rough
proxy for metallicity. We linearly convert the color and magnitude of each
star to the index (0 to 1) and assign the color (dark blue to orange). The
bluest (index$=0$) and reddest (index$=1$) colors correspond to
$\rm{[M/H]}=-2.2$ and $-0.4$, respectively, assuming 10 Gyr old age (Bressan
et al., 2012). The map reveals a very extended and asymmetric metal-rich
stellar halo that completely fills the central HSC pointing
($\textrm{R}_{proj}\sim 50$ kpc) and which will be the subject of a future
paper.

The map also reveals two very prominent metal-poor substructures at the edges
of the current SNFC survey coverage. While the northern feature is the
satellite galaxy, MM-Dw2, discovered independently by Toloba et al. (2016) and
Romanowsky et al. (2016), the southern over-density, which we call SNFC-dw1,
is new. Another known satellite, Scl-MM-dw1, is located near SNFC-dw1, but it
is hardly visible on the map due to its compactness ($\rm{R_{h}}=16.8\arcsec$)
(Sand et al., 2014). To investigate the properties of SNFC-dw1, we first set
the center by eye and isolate a surrounding circular region of 9 arcmin
radius, as shown by the black solid circle. We also select a reference field
of identical area at the same elliptical radius from NGC 253’s center as SNFC-
dw1, shown by the black dashed circle.

The lower-middle panel of Figure 1 presents the zoom-in view around SNFC-dw1.
It shows the smoothed density map of metal-poor RGB stars within the leftward
of the dashed line in the polygonal region in Figure 2(b). The kernel density
is estimated with the bandwidth of 1.′0, which is $\sim 1$ kpc at the distance
of NGC 253. There are two areas affected by narrow diffraction spikes from
saturated foreground stars where the pipeline failed to return reliable
photometry. To reduce any spurious effects, we masked these regions, shown
with gray shading, and re-populated them with stars drawn from the adjacent
area before estimating SNFC-dw1’s properties. The density map shows the
flattened core of SNFC-dw1 with a surrounding asymmetrical extension.
Specifically, the stellar density rapidly decreases on the southern side,
while on the northern side two elongated tidal extensions can be seen that
point towards NGC 253, indicated by the white arrow in the panel. This
morphology indicates that SNFC-dw1 is highly-disrupted and that we are likely
viewing it at apocenter, with both the leading and trailing tails projecting
on the same side of its main body.

Figure 2 shows the de-reddened Hess diagrams of stellar sources around NGC
253, SNFC-dw1, the reference field, and the field-subtracted SNFC-dw1. The
error bars in panels (a) to (c) show the median photometric error of stars in
the interval $-1.5<(g-i)_{0}<3.0$. The 50 percent completeness threshold is
shown as the dotted line. Panel (a) contains stars within an elliptical
annulus of $2\times\rm{R_{25}}<$ r $<4\times\rm{R_{25}}$ of NGC 253, adopting
the $\rm{R_{25}}=13.7\arcmin$ (de Vaucouleurs et al., 1991). PARSEC isochrones
of age 10 Gyr and $\rm{[M/H]}=-2.2,-1.8,-1.4,-1.0,-0.7,-0.4$ (Bressan et al.,
2012) are shifted to NGC 253’s distance and overlaid as solid lines. A well-
populated broad RGB is seen ($i_{0}>24$, $(g-i)_{0}>1$) that spans the full
range of isochrones shown. Another obvious feature is the blue plume at
$-0.5<(g-i)_{0}<0.7$, which we identify with unresolved high-z galaxies as
seen in other wide-field surveys of nearby galaxies (e.g. Barker et al., 2012;
Okamoto et al., 2019) as well as in the reference field (Figure 2(c)).
Comparing Figure 2(b) that shows stars within a circular radius of 9 arcmin
centered on SNFC-dw1 with Figure 2(c) that shows the equal area reference
field, the signature of a metal-poor RGB sequence is striking. This feature is
well-matched to a theoretical isochrone of $\rm{[M/H]}=-1.8$ and 10 Gyr old. A
tantalising feature in panel (b) is the plume of stars that sits above the RGB
selection box that we tentatively associate with a thermally-pulsing
asymptotic giant branch (TP-AGB) population. We select these stars using the
orange solid parallelogram in panel (b), and plot their spatial distribution
as the black points in the lower-middle panel of Figure 1. These sources are
indeed visually clustered on the main body of SNFC-dw1 as would be expected
for an instrinsic population. The field-subtracted Hess diagram in Figure 2(d)
confirms the clear RGB overdensity as well as the existence of TP-AGB stars in
the SNFC-dw1 field.

###  3.2 Distance

We derive the distance to SNFC-dw1 using the TRGB method (e.g. Lee et al.,
1993), following the procedure described in Okamoto et al. (2019). Stars of
$1.2<(g-i)_{0}<2.25$ that lie within 9 arcmin radius are used to derive the
$i$-band luminosity function (LF). Using the LF of point sources in the
reference field, we correct for contamination from the NGC 253 halo population
and from foreground/background objects. Then, we apply a Sobel filter to
detect a sharp transition at $i_{0}=24.18\pm 0.08$. We estimate the TRGB color
as $(g-i)_{0}=1.63\pm 0.08$ and apply equation (1) of Okamoto et al. (2019) to
calculate $\rm{M}_{\textrm{i,TRGB}}$ as $-3.62\pm 0.10$. Therefore, the
distance modulus to SNFC-dw1 is $(m-M)_{0}=27.79\pm 0.12$, corresponding to D
$=3.62\pm 0.2$ Mpc. To be able to directly compare, we also estimate the TRGB
distance of NGC 253 using our dataset. In this case, we use stars lying in the
elliptical annulus covering $3-4\times\rm{R_{25}}$ and correct for
contamination of foreground/background objects using the same reference field
but scaled up to cover the same area within the elliptical annulus. We find
$(m-M)_{0}=27.72\pm 0.14$, or D $=3.50\pm 0.22$ Mpc, which is in excellent
agreement with Radburn-Smith et al. (2011). This places SNFC-dw1 at a 3D
radius of $\sim 146$ kpc on the far side of NGC 253.

###  3.3 Luminosity, Structure and Metallicity

The total luminosity and color of SNFC-dw1 are estimated using the flux of
metal-poor RGB stars within the solid circle of the right panel in Figure 1.
The photometric completeness of a given star is inferred from comparing its
magnitude and color with the interpolated completeness curve. The flux is then
scaled to account for stars fainter than the RGB selection box using the best-
fit theoretical isochrone. The contaminant contribution is calculated by
applying the same approach to the reference field, and then subtracted. We
find $\rm{M_{g}}=-11.15\pm 0.2$ and $(g-i)_{0}=1.35\pm 0.3$, which corresponds
to $\rm{M_{V}}=-11.7\pm 0.2$ using the transformations in Komiyama et al.
(2018). Adopting equation (8) of Taylor et al. (2011), we calculate the
stellar mass of SNFC-dw1 to be $\sim 1.25\times 10^{7}~{}\rm{M}_{\odot}$.

Figure 3: Metallicity distribution of RGB stars derived using the $(g-i)_{0}$
color. The open histogram shows the completeness-corrected MDF, and the shaded
histogram illustrates the contaminant-subtracted MDF using the reference
field. The solid line shows a Gaussian fit to the contaminant-subtracted MDF.

Figure 3 shows the metallicity distribution function (MDF) of SNFC-dw1 derived
from the completeness-corrected color distribution of bright RGB stars. We use
RGB stars of $i_{0}<25.0$ to reduce the uncertainties and to keep the wide
spacing between isochrones in the $(g-i)_{0}$ color. RGB stars are selected
within the solid circle of the right panel of Figure 1 and within the solid
polygon of Figure 2. The metallicity is derived from the linear interpolation
of isochrones of $\rm{[M/H]}=-2.2$ to $-0.7$ with a fixed age of 10 Gyr and
extrapolation with the spline function up to $\rm{[M/H]}=-2.5$ in the same
manner as done in Okamoto et al. (2023). The completeness is corrected as
before according to the color and magnitude of each star. The same estimation
is done for the reference field in order to derive the contaminant
contribution to the MDF arising from the NGC 253 halo and
background/foreground objects. The contaminant-corrected MDF is shown as the
shaded histogram. The uncertainty is estimated by performing a Monte Carlo
simulation with N $=2000$ for each RGB star. A star is randomly re-sampled
from a Gaussian distribution with a width equal to the photometric error, and
the metallicity is redetermined. The standard deviation of the resulting
metallicity distribution is adopted as the uncertainty for each RGB star. The
estimated uncertainties increase from 0.03 to 0.14 dex with decreasing
metallicity from $\rm{[M/H]}=-0.4$ to $-2.4$ due to the narrow range in RGB
color between metal-poor isochrones. The weighted mean metallicity corrected
for this effect is shown as the error-weighted mean value
$\langle\rm{[M/H]}_{w}\rangle=-1.51\pm 0.1$ dex. The overall appearance of the
MDF is well-described by a single Gaussian shown as the solid line in Figure
3. The sudden drop on the metal-rich side ($\rm{[M/H]}>-1.3$) of the MDF
likely reflects an over-subtraction due to small number statistics in this
metallicity range, in which the NGC 253 halo dominates the contaminant signal.

Table 1: The properties of NGC253-SNFC-dw1

Parameter | Value
---|---
R.A.(J2000) | $00^{h}48^{m}39^{s}.68$
Dec. (J2000) | $-26\arcdeg 33\arcmin 48\arcsec.7$
Distance |  $3.62\pm 0.2$ Mpc
$(m-M)_{0}$ | $27.79\pm 0.12$
P.A.(a) | 104.8$\arcdeg$
$\epsilon^{(b)}$ | 0.06
$\rm{M_{g}}^{(c)}$ | $-11.15\pm 0.2$
$(g-i)_{0}^{(c)}$ | $1.35\pm 0.3$
$\rm{M_{V}}^{(c)}$ | $-11.7\pm 0.2$
$n^{(d)}$ | $0.45\pm 0.04$
$\rm{R_{h}}^{(e)}$ | $3.20\arcmin\pm 0.16$
$\rm{R_{h}}^{(e)}$ |  $3.37\pm 0.36$ kpc
$\langle\rm{[M/H]}_{w}\rangle^{(c)}$ |  $-1.51\pm 0.10$ dex
$\sigma_{\rm{[M/H]}}^{(c)}$ |  $0.38$ dex
$\langle\mu_{g}\rangle(r<\rm{R_{h}})^{(f)}$ |  $\sim 30.1$ mag arcsec-2
$\textrm{M}_{*}^{(g)}$ | $\sim 1.25\times 10^{7}\rm{M}_{\odot}$

Note. — (a) Position angle from north to east. (b) Ellipticity
$\epsilon=1-b/a$ where $b/a$ is the axis ratio. (c) Measured using stars
within a circular radius of 9 arcmin at the SNFC-dw1. (d) Sérsic index. (e)
Half-light radius of the Sérsic profile. (f) Average surface brightness within
the half-light radius in g-band. (g) The stellar mass calculated using the
inferred total magnitude.

Although the extended regions of SNFC-dw1 are very irregular, the core
structure seen in the lower-middle panel of Figure 1 is well-defined and so we
proceed to estimate its structural properties. We compute the values using the
density-weighted zeroth, first and second moments of the metal-poor RGB
spatial distribution and list them in Table 1. The estimated shape shown as a
dashed line in Figure 1 seems to be more circular than the visual impression
of the density map, probably due to the existence of extended components
toward the north. We also fit the standard Sérsic profile via least-squares
minimization to the completeness-corrected metal-poor RGB star counts. The
radial profile is constructed by calculating the average number density of RGB
stars in a series of elliptical annuli, defined using the derived structural
parameters. We take the Poisson uncertainties in the RGB counts as the
uncertainty. Initially, the contaminant contribution is visually estimated as
the average number density at large radius ($>10\arcmin$) and we fit the
profile and estimate the half-light radius ($\rm{R_{h}}$). We then iterate on
the estimation of the contaminant level, this time using the average number
density between 4 and $6\times\rm{R_{h}}$. Finally, the profile is refit with
this contamination level. The resulting $\rm{R_{h}}=3.37\pm 0.36$ kpc is very
large for a dwarf galaxy and, in combination with the relatively low
luminosity, yields an extremely faint average surface brightness of
$\mu_{g}\sim 30.1$ mag arcsec-2 within $\rm{R_{h}}$.

Figure 4: Upper: The luminosity-size relation. The purple symbols are known
NGC 253 satellites (Martínez-Delgado et al., 2021; Mutlu-Pakdil et al., 2022).
The yellow squares are UDGs in the Coma cluster (van Dokkum et al., 2015). The
filled circles show LG dSphs from McConnachie (2012, updated in January 2021).
The open squares with names are particularly diffuse dSphs (Toloba et al.,
2016; Torrealba et al., 2019; Crnojević et al., 2019; Collins et al., 2020; Ji
et al., 2021; Savino et al., 2022; Žemaitis et al., 2023). Lower: The
luminosity-metallicity relation of nearby dSphs. The filled and open black
squares represent SNFC-dw1 with $\rm{[\alpha/Fe]}=0.0$ and 0.32 dex,
respectively. The solid and dotted lines show the Kirby et al. (2013) relation
with the $1\sigma$ deviation.

##  4 Discussion and Summary

Using data from the SNFC survey being conducted with the Subaru Telescope, we
have uncovered a new dwarf satellite of the nearby spiral NGC 253. The system
lies at a projected separation of 83.4 kpc from the center of NGC 253. We have
used the TRGB method to derive consistent distances to SNFC-dw1 and NGC 253,
finding them to lie at $3.62\pm 0.2$ Mpc and $3.5\pm 0.2$ Mpc, respectively.
This places SNFC-dw1 at a 3D radius of 146 kpc on the far side of the NGC 253
halo and well within its estimated virial radius of 330 kpc (Mutlu-Pakdil et
al., 2021), confirming it is a bound satellite.

We find SNFC-dw1 to be dominated by an old metal-poor population;
additionally, there is tentative evidence for an associated TP-AGB population.
As such stars have ages of typically $0.6-2$ Gyr (Marigo et al., 2017), their
presence suggests that the system was star-forming until relatively recently.
Although SNFC-dw1 has a comparable luminosity, stellar mass and metallicity to
bright classical Galactic dSphs, it is extremely diffuse ($\mu_{g}\sim 30.1$
mag arcsec-2) – a direct consequence of its abnormally large size with
$\rm{R_{h}}=3.37\pm 0.36$ kpc. Indeed, star counts can be traced over more
than 10 kpc and reveal a very asymmetric structure in the peripheral regions.
Notably, SNFC-dw1 has a size comparable to some of the most diffuse galaxies
currently known, including the puzzling Local Group systems Crater 2
($\rm{R_{h}}=1.1$ kpc), Antlia 2 ($\rm{R_{h}}=2.5$ kpc) and Andromeda XIX
($\rm{R_{h}}=3.4$ kpc) (Torrealba et al., 2016; Ji et al., 2021; Savino et
al., 2022).

Figure 4 compares SNFC-dw1 to other dwarf galaxies in terms of scaling
relationships. These plots include LG dSphs (McConnachie, 2012, updated in
January 2021), Ultra-diffuse galaxies (UDGs) in the Coma Cluster (van Dokkum
et al., 2015), some of the known NGC 253 satellites (Toloba et al., 2016;
Martínez-Delgado et al., 2021; Mutlu-Pakdil et al., 2022) and the particularly
diffuse dSphs (open squares) Antlia 2 (Torrealba et al., 2019; Ji et al.,
2021; Vivas et al., 2022), Andromeda XIX (Collins et al., 2020), the M81
satellite F8D1 (Žemaitis et al., 2023), the Centauras A satellite CenA-MM-dw1
(Crnojević et al., 2019), and the recently-discovered NGC 55-dw1 (McNanna et
al., 2024). The upper panel of Figure 4 shows the systems on the size-
luminosity plane. Only six systems have a similarly large size
($\rm{R_{h}}\gtrsim 1.5$ kpc) to SNFC-dw1 – Andromeda XIX and Antlia 2 and NGC
55-dw1, which are $\gtrsim 2$ magnitudes fainter, and the Sagittarius (Sgr)
dSph, CenA-MM-dw3 and F8D1, which are $\gtrsim 2$ magnitudes brighter. While
both Sgr, CenA-MM-dw3 and F8D1 show direct evidence for tidal disruption
(Crnojević et al., 2016; Majewski et al., 2003; Žemaitis et al., 2023), the
evidence in Antlia 2 and Andromeda XIX has thus far only been circumstantial
(e.g. Belokurov et al., 2006; Torrealba et al., 2019; Collins et al., 2020; Ji
et al., 2021). The situation is also not clear for NGC 55-dw1 which, in
addition to lacking clear tidal features, does not yet have a directly
measured distance (McNanna et al., 2024). Amongst the NGC 253 satellites, MM-
dw2 is most similar to SNFC-dw1 and it exhibits a highly elongated shape (see
Figure 1). However, the properties of MM-dw2 estimated from two independent
studies are rather different – $\rm{M_{V}}=-12.1$ and $\rm{R_{h}}=2.94$ kpc
using Magellan/Megacam (Toloba et al., 2016) and $\rm{M_{V}}=-10.7$ and
$\rm{R_{h}}=1.1$ kpc using Subaru/Suprime-Cam (Romanowsky et al., 2016). We
plot both values connected with a purple solid line in the upper panel of
Figure 4. The prominent tidal extensions emanating from SNFC-dw1 provide
unambiguous evidence of tidal disruption, demonstrating that this is a natural
pathway to produce such extremely diffuse and extended galaxies. This
conclusion is further supported by integrated light analyses of two more
distant dwarfs with similarly large sizes – NGC4449B associated with the dwarf
starburst galaxy NGC 4449 (Rich et al., 2012; Martínez-Delgado et al., 2012)
and HCC-087 in the Hydra I galaxy cluster, (Koch et al., 2012), both of which
are also known to exhibit S-shaped tidal features.

The lower panel of Figure 4 shows the luminosity-metallicity scaling relation
(LZR) defined by local dwarfs. To place SNFC-dw1 on this plot, we converted
the global metallicity $\rm{[M/H]}$ to $\rm{[Fe/H]}$ using equation (3) of
Salaris et al. (1993) assuming two cases, $\rm{[\alpha/Fe]}=0.0$ and 0.32 dex,
the latter inferred from recent spectroscopic studies of UDGs (Ferré-Mateu et
al., 2023). Broadly speaking, SNFC-dw1 obeys the same LZR as LG dSphs (Kirby
et al., 2013). If we assume an $\alpha$-enhanced abundance, it is somewhat
more metal-poor than expected for a galaxy of the same luminosity but it still
lies within $\sim 1\sigma$ of the relation, suggesting that it has not yet
lost a significant fraction of its original stellar mass. Whether or not this
finding is at odds with the highly disrupted appearance of SNFC-dw1 is
unclear. A measurement of the velocity dispersion would be extremely valuable
for constraining the degree of tidal stripping (e.g. Borukhovetskaya et al.,
2022) but such observations would be very challenging with current facilities
given the faintness of SNFC-dw1’s stars.

It is worth noting that SNFC-dw1 lies in an area of the sky that has been the
target of multiple dwarf galaxy searches in recent years, some of which have
used integrated light (Martínez-Delgado et al., 2021; Mutlu-Pakdil et al.,
2022) while others have used resolved stars (Carlsten et al., 2022). Although
it is the second most luminous of the six satellites projected within 150 kpc
of NGC 253 (Mutlu-Pakdil et al., 2024), and comparable to the brightest
classical dSph of the Milky Way and Andromeda, SNFC-dw1 has remained unearthed
until now due to its extremely diffuse nature. This emphasizes the importance
of surface brightness limiting depth in satellite searches, suggesting that
even bright-end incompleteness can be an issue at the typical depths of most
present-day studies. Fortunately, with the recent launch of ESA’s Euclid
satellite, the start of the Legacy Survey of Space and Time with the Vera C.
Rubin Observatory in 2025 and the Nancy Grace Roman Space Telescope launch in
2027, the era of very sensitive ($\mu>30$ mag arcsec-2) satellite searches
around large numbers of Local Volume galaxies will soon be upon us.

We are grateful to all the staff at Subaru Telescope and the HSC team. This
research is based on data collected at the Subaru Telescope, operated by the
National Astronomical Observatory of Japan. We are honored and grateful for
the opportunity to observe the Universe from Maunakea, which has cultural,
historical, and natural significance in Hawaii. This paper makes use of
software developed for Vera C. Rubin Observatory. We thank the Rubin
Observatory for making their code available as free software at
http://pipelines.lsst.io/. The Pan-STARRS1 Surveys (PS1) and the PS1 public
science archive have been made possible through contributions by the Institute
for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the
Max-Planck Society and its participating institutes, the Max Planck Institute
for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial
Physics, Garching, The Johns Hopkins University, Durham University, the
University of Edinburgh, the Queen’s University Belfast, the Harvard-
Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global
Telescope Network Incorporated, the National Central University of Taiwan, the
Space Telescope Science Institute, the National Aeronautics and Space
Administration under Grant No. NNX08AR22G issued through the Planetary Science
Division of the NASA Science Mission Directorate, the National Science
Foundation Grant No. AST-1238877, the University of Maryland, Eotvos Lorand
University (ELTE), the Los Alamos National Laboratory, and the Gordon and
Betty Moore Foundation. S.O. acknowledges support in part from JSPS Grant-in-
Aid for Scientific Research (18H05875, 20K04031, 20H05855). M.G.L. was
supported by the National Research Foundation grant funded by the Korean
Government (NRF-2019R1A2C2084019).

## References

  * Axelrod et al. (2010) Axelrod, T., Kantor, J., Lupton, R. H., & Pierfederici, F. 2010, in Proc. SPIE, Vol. 7740, Software and Cyberinfrastructure for Astronomy, 774015, doi: 10.1117/12.857297
  * Bailin et al. (2011) Bailin, J., Bell, E. F., Chappell, S. N., Radburn-Smith, D. J., & de Jong, R. S. 2011, ApJ, 736, 24, doi: 10.1088/0004-637X/736/1/24
  * Barker et al. (2012) Barker, M. K., Ferguson, A. M. N., Irwin, M. J., Arimoto, N., & Jablonka, P. 2012, MNRAS, 419, 1489, doi: 10.1111/j.1365-2966.2011.19814.x
  * Beck et al. (1982) Beck, R., Hutschenreiter, G., & Wielebinski, R. 1982, A&A, 106, 112
  * Belokurov et al. (2006) Belokurov, V., Zucker, D. B., Evans, N. W., et al. 2006, ApJ, 642, L137, doi: 10.1086/504797
  * Bonnarel et al. (2000) Bonnarel, F., Fernique, P., Bienaymé, O., et al. 2000, A&AS, 143, 33, doi: 10.1051/aas:2000331
  * Borukhovetskaya et al. (2022) Borukhovetskaya, A., Navarro, J. F., Errani, R., & Fattahi, A. 2022, MNRAS, 512, 5247, doi: 10.1093/mnras/stac653
  * Bosch et al. (2018) Bosch, J., Armstrong, R., Bickerton, S., et al. 2018, PASJ, 70, S5, doi: 10.1093/pasj/psx080
  * Bressan et al. (2012) Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127, doi: 10.1111/j.1365-2966.2012.21948.x
  * Bullock & Boylan-Kolchin (2017) Bullock, J. S., & Boylan-Kolchin, M. 2017, ARA&A, 55, 343, doi: 10.1146/annurev-astro-091916-055313
  * Bullock & Johnston (2005) Bullock, J. S., & Johnston, K. V. 2005, ApJ, 635, 931, doi: 10.1086/497422
  * Cantiello et al. (2018) Cantiello, M., Grado, A., Rejkuba, M., et al. 2018, A&A, 611, A21, doi: 10.1051/0004-6361/201731325
  * Carlsten et al. (2022) Carlsten, S. G., Greene, J. E., Beaton, R. L., Danieli, S., & Greco, J. P. 2022, ApJ, 933, 47, doi: 10.3847/1538-4357/ac6fd7
  * Chambers et al. (2016) Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv e-prints, arXiv:1612.05560, doi: 10.48550/arXiv.1612.05560
  * Chiboucas et al. (2013) Chiboucas, K., Jacobs, B. A., Tully, R. B., & Karachentsev, I. D. 2013, AJ, 146, 126, doi: 10.1088/0004-6256/146/5/126
  * Collins et al. (2020) Collins, M. L. M., Tollerud, E. J., Rich, R. M., et al. 2020, MNRAS, 491, 3496, doi: 10.1093/mnras/stz3252
  * Crnojević et al. (2016) Crnojević, D., Sand, D. J., Spekkens, K., et al. 2016, ApJ, 823, 19, doi: 10.3847/0004-637X/823/1/19
  * Crnojević et al. (2019) Crnojević, D., Sand, D. J., Bennet, P., et al. 2019, ApJ, 872, 80, doi: 10.3847/1538-4357/aafbe7
  * Davidge (2010) Davidge, T. J. 2010, ApJ, 725, 1342, doi: 10.1088/0004-637X/725/1/1342
  * de Vaucouleurs et al. (1991) de Vaucouleurs, G., de Vaucouleurs, A., Corwin, Herold G., J., et al. 1991, Third Reference Catalogue of Bright Galaxies
  * Doliva-Dolinsky et al. (2023) Doliva-Dolinsky, A., Martin, N. F., Yuan, Z., et al. 2023, ApJ, 952, 72, doi: 10.3847/1538-4357/acdcf6
  * Ferré-Mateu et al. (2023) Ferré-Mateu, A., Gannon, J. S., Forbes, D. A., et al. 2023, MNRAS, 526, 4735, doi: 10.1093/mnras/stad3102
  * Fitzpatrick (1999) Fitzpatrick, E. L. 1999, PASP, 111, 63, doi: 10.1086/316293
  * Greggio et al. (2014) Greggio, L., Rejkuba, M., Gonzalez, O. A., et al. 2014, A&A, 562, A73, doi: 10.1051/0004-6361/201322759
  * Harmsen et al. (2017) Harmsen, B., Monachesi, A., Bell, E. F., et al. 2017, MNRAS, 466, 1491, doi: 10.1093/mnras/stw2992
  * Homma et al. (2019) Homma, D., Chiba, M., Komiyama, Y., et al. 2019, PASJ, 71, 94, doi: 10.1093/pasj/psz076
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111, doi: 10.3847/1538-4357/ab042c
  * Ji et al. (2021) Ji, A. P., Koposov, S. E., Li, T. S., et al. 2021, ApJ, 921, 32, doi: 10.3847/1538-4357/ac1869
  * Jurić et al. (2017) Jurić, M., Kantor, J., Lim, K. T., et al. 2017, in Astronomical Society of the Pacific Conference Series, Vol. 512, Astronomical Data Analysis Software and Systems XXV, ed. N. P. F. Lorente, K. Shortridge, & R. Wayth, 279, doi: 10.48550/arXiv.1512.07914
  * Karachentsev et al. (2013) Karachentsev, I. D., Makarov, D. I., & Kaisina, E. I. 2013, AJ, 145, 101, doi: 10.1088/0004-6256/145/4/101
  * Kirby et al. (2013) Kirby, E. N., Cohen, J. G., Guhathakurta, P., et al. 2013, ApJ, 779, 102, doi: 10.1088/0004-637X/779/2/102
  * Koch et al. (2012) Koch, A., Burkert, A., Rich, R. M., et al. 2012, ApJ, 755, L13, doi: 10.1088/2041-8205/755/1/L13
  * Komiyama et al. (2018) Komiyama, Y., Chiba, M., Tanaka, M., et al. 2018, ApJ, 853, 29, doi: 10.3847/1538-4357/aaa129
  * Lee et al. (1993) Lee, M. G., Freedman, W. L., & Madore, B. F. 1993, ApJ, 417, 553, doi: 10.1086/173334
  * Magnier et al. (2013) Magnier, E. A., Schlafly, E., Finkbeiner, D., et al. 2013, ApJS, 205, 20, doi: 10.1088/0067-0049/205/2/20
  * Majewski et al. (2003) Majewski, S. R., Skrutskie, M. F., Weinberg, M. D., & Ostheimer, J. C. 2003, ApJ, 599, 1082, doi: 10.1086/379504
  * Malin & Hadley (1997) Malin, D., & Hadley, B. 1997, PASA, 14, 52, doi: 10.1071/AS97052
  * Marigo et al. (2017) Marigo, P., Girardi, L., Bressan, A., et al. 2017, ApJ, 835, 77, doi: 10.3847/1538-4357/835/1/77
  * Martínez-Delgado et al. (2012) Martínez-Delgado, D., Romanowsky, A. J., Gabany, R. J., et al. 2012, ApJ, 748, L24, doi: 10.1088/2041-8205/748/2/L24
  * Martínez-Delgado et al. (2021) Martínez-Delgado, D., Makarov, D., Javanmardi, B., et al. 2021, A&A, 652, A48, doi: 10.1051/0004-6361/202141242
  * McConnachie (2012) McConnachie, A. W. 2012, AJ, 144, 4, doi: 10.1088/0004-6256/144/1/4
  * McNanna et al. (2024) McNanna, M., Bechtol, K., Mau, S., et al. 2024, ApJ, 961, 126, doi: 10.3847/1538-4357/ad07d0
  * Miyazaki et al. (2018) Miyazaki, S., Komiyama, Y., Kawanomoto, S., et al. 2018, PASJ, 70, S1, doi: 10.1093/pasj/psx063
  * Mutlu-Pakdil et al. (2021) Mutlu-Pakdil, B., Sand, D. J., Crnojević, D., et al. 2021, ApJ, 918, 88, doi: 10.3847/1538-4357/ac0db8
  * Mutlu-Pakdil et al. (2022) —. 2022, ApJ, 926, 77, doi: 10.3847/1538-4357/ac4418
  * Mutlu-Pakdil et al. (2024) —. 2024, arXiv e-prints, arXiv:2401.14457, doi: 10.48550/arXiv.2401.14457
  * Ogami et al. (2024) Ogami, I., Tanaka, M., Komiyama, Y., et al. 2024, arXiv e-prints, arXiv:2401.00668, doi: 10.48550/arXiv.2401.00668
  * Okamoto et al. (2019) Okamoto, S., Arimoto, N., Ferguson, A. M. N., et al. 2019, ApJ, 884, 128, doi: 10.3847/1538-4357/ab44a7
  * Okamoto et al. (2023) Okamoto, S., Arimoto, N., Ferguson, A. M. N., Irwin, M. J., & Žemaitis, R. 2023, ApJ, 952, 77, doi: 10.3847/1538-4357/acdad1
  * Pillepich et al. (2014) Pillepich, A., Vogelsberger, M., Deason, A., et al. 2014, MNRAS, 444, 237, doi: 10.1093/mnras/stu1408
  * Pucha et al. (2019) Pucha, R., Carlin, J. L., Willman, B., et al. 2019, ApJ, 880, 104, doi: 10.3847/1538-4357/ab29fb
  * Radburn-Smith et al. (2011) Radburn-Smith, D. J., de Jong, R. S., Seth, A. C., et al. 2011, ApJS, 195, 18, doi: 10.1088/0067-0049/195/2/18
  * Rich et al. (2012) Rich, R. M., Collins, M. L. M., Black, C. M., et al. 2012, Nature, 482, 192, doi: 10.1038/nature10837
  * Romanowsky et al. (2016) Romanowsky, A. J., Martínez-Delgado, D., Martin, N. F., et al. 2016, MNRAS, 457, L103, doi: 10.1093/mnrasl/slv207
  * Salaris et al. (1993) Salaris, M., Chieffi, A., & Straniero, O. 1993, ApJ, 414, 580, doi: 10.1086/173105
  * Sales et al. (2022) Sales, L. V., Wetzel, A., & Fattahi, A. 2022, Nature Astronomy, 6, 897, doi: 10.1038/s41550-022-01689-w
  * Sand et al. (2014) Sand, D. J., Crnojević, D., Strader, J., et al. 2014, ApJ, 793, L7, doi: 10.1088/2041-8205/793/1/L7
  * Savino et al. (2022) Savino, A., Weisz, D. R., Skillman, E. D., et al. 2022, ApJ, 938, 101, doi: 10.3847/1538-4357/ac91cb
  * Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103, doi: 10.1088/0004-637X/737/2/103
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525, doi: 10.1086/305772
  * Taylor et al. (2011) Taylor, E. N., Hopkins, A. M., Baldry, I. K., et al. 2011, MNRAS, 418, 1587, doi: 10.1111/j.1365-2966.2011.19536.x
  * Toloba et al. (2016) Toloba, E., Sand, D. J., Spekkens, K., et al. 2016, ApJ, 816, L5, doi: 10.3847/2041-8205/816/1/L5
  * Torrealba et al. (2016) Torrealba, G., Koposov, S. E., Belokurov, V., & Irwin, M. 2016, MNRAS, 459, 2370, doi: 10.1093/mnras/stw733
  * Torrealba et al. (2019) Torrealba, G., Belokurov, V., Koposov, S. E., et al. 2019, MNRAS, 488, 2743, doi: 10.1093/mnras/stz1624
  * van Dokkum et al. (2015) van Dokkum, P. G., Abraham, R., Merritt, A., et al. 2015, ApJ, 798, L45, doi: 10.1088/2041-8205/798/2/L45
  * Vivas et al. (2022) Vivas, A. K., Martínez-Vázquez, C. E., Walker, A. R., et al. 2022, ApJ, 926, 78, doi: 10.3847/1538-4357/ac43bd
  * Žemaitis et al. (2023) Žemaitis, R., Ferguson, A. M. N., Okamoto, S., et al. 2023, MNRAS, 518, 2497, doi: 10.1093/mnras/stac3133


# Simple Mitigation of Global Depolarizing Errors in Quantum Simulations

Joseph Vovrosh  Blackett Laboratory, Imperial College London, London SW7 2AZ,
United Kingdom  Kiran E. Khosla  Blackett Laboratory, Imperial College London,
London SW7 2AZ, United Kingdom  Sean Greenaway  Blackett Laboratory, Imperial
College London, London SW7 2AZ, United Kingdom  Christopher Self  Blackett
Laboratory, Imperial College London, London SW7 2AZ, United Kingdom  M. S. Kim
Blackett Laboratory, Imperial College London, London SW7 2AZ, United Kingdom
Johannes Knolle  Department of Physics TQM, Technische Universität München,
James-Franck-Straße 1, D-85748 Garching, Germany  Munich Center for Quantum
Science and Technology (MCQST), 80799 Munich, Germany  Blackett Laboratory,
Imperial College London, London SW7 2AZ, United Kingdom

###### Abstract

To get the best possible results from current quantum devices error mitigation
is essential. In this work we present a simple but effective error mitigation
technique based on the assumption that noise in a deep quantum circuit is well
described by global depolarizing error channels. By measuring the errors
directly on the device, we use an error model ansatz to infer error-free
results from noisy data. We highlight the effectiveness of our mitigation via
two examples of recent interest in quantum many-body physics: entanglement
measurements and real time dynamics of confinement in quantum spin chains. Our
technique enables us to get quantitative results from the IBM quantum
computers showing signatures of confinement, i.e. we are able to extract the
meson masses of the confined excitations which were previously out of reach.
Additionally, we show the applicability of this mitigation protocol in a wider
setting with numerical simulations of more general tasks using a realistic
error model. Our protocol is device-independent, simply implementable and
leads to large improvements in results if the global errors are well described
by depolarization.

##  I Introduction

Quantum computers are becoming large enough ($\sim 50-100$ qubits Preskill
(2018)) to, in principle, allow demonstrations of their ‘Quantum Advantage’
Boixo  _et al._ (2018); Arute  _et al._ (2019). However, the actual amount of
entanglement that can be generated in current devices is constrained by noise
and errors, limiting their ability to solve complex problems such as quantum
simulation. To address this, various _error mitigation_ strategies have
recently been developed to counteract noise and boost the fidelity of
experimental results.

Error mitigation differs from fault tolerance. Fault tolerant quantum
computers will eventually be able to suppress errors by encoding quantum
information in a redundantly large number of qubits as error correcting codes
Steane (1996); Shor (1995). In this way they will be able to execute
arbitrarily deep circuits by the repeated application of active error
corrections. Unfortunately, these encodings cannot be used in current devices
as they require smaller hardware errors and larger numbers of qubits than are
currently available. In contrast, error mitigation strategies are applied to
unencoded physical qubits. Rather than actively correcting errors, they aim to
estimate what the effect of the error was and infer the error-free result.
This current stage in the development of quantum computers has been dubbed the
noisy intermediate-scale quantum (NISQ) era Preskill (2018) and is expected to
last for the foreseeable future.

Figure 1: Results from the IBM device Paris par  of the second order Rényi
entanglement entropy. Quench dynamics before and after error mitigation of the
second order Réyni entropy for the transverse field Ising model with varying
longitudinal field strengths. Here, $J=1$, $h_{x}=0.5$ and $n=6$. Note, $800$
random unitaries were used to collect this data giving an uncertainty in the
purity of $\delta\rho\sim O(10^{-2})$. Clearly the mitigation allows results
to go from qualitative agreement to quantitative agreement with the results
obtained numerically through exact diagonalization (ED).

In the past few years, error mitigation has been a thriving research direction
as more and more quantum devices become available. Error mitigation strategies
typically address measurement errors qis (2019); Jattana  _et al._ (2020);
Maciejewski  _et al._ (2020); Bravyi  _et al._ (2020); Chen  _et al._ (2019);
McClean  _et al._ (2020a); Bonet-Monroig  _et al._ (2018); Sagastizabal  _et
al._ (2019); Tran  _et al._ (2020) or the algorithms and gates employed for
digital quantum simulation Huggins  _et al._ (2019); McArdle  _et al._ (2019);
Dumitrescu  _et al._ (2018); McClean  _et al._ (2017); Colless  _et al._
(2018); Kandala  _et al._ (2019); Otten and Gray (2019a, b). While some of
these techniques have proven useful for certain quantum algorithms, their
general success is hampered by the fact that they either rely on a high level
of control of the quantum device itself, or are specifically designed for a
given quantum simulation problem, i.e. exploiting specific symmetries Smith
_et al._ (2019a).

One promising direction is to employ machine learning algorithms Strikis  _et
al._ (2020); Zlokapa and Gheorghiu (2020); Czarnik  _et al._ (2020) for error
mitigation. Generally speaking, these methods train classical computers to
predict the error seen in a quantum device and use the results to infer the
error-free quantum simulations. While successful, these methods require a
large increase in the classical computational overhead and are somewhat
uncontrolled. The latter drawback is also true for popular protocols based on
the idea of increasing errors in the device systematically and then
extrapolating back to the zero error case Kandala  _et al._ (2019); Li and
Benjamin (2017); Endo  _et al._ (2018); Temme  _et al._ (2017); Gustafson  _et
al._ (2019); Giurgica-Tiron  _et al._ (2020); Cai (2020); He  _et al._ (2020).
In general, how to tune the error rates varies from device to device and
reliable fitting requires expert knowledge of the specific hardware.

Here, we propose a new protocol for gate error mitigation combining a raft of
desirable features: it is easily implementable on any quantum device with
little increase in workload, it is well rooted in a mathematical description
of the errors, and it is suitable for any quantum algorithm of interest.
Furthermore, as we show via specific examples, it can lead to to large
improvements in the performance of quantum devices, for example see Fig. 1 for
Rényi entropy results.

The paper is organised as follows: First, we derive an ansatz for the density
matrix resulting from the action of a noisy quantum simulator with global
depolarizing errors. Following this, we explain how this ansatz can be used
for a general error mitigation protocol. We then demonstrate its effectiveness
by studying real time dynamics of the transverse field Ising model (TFIM) with
a longitudinal field. In that context, some of the authors recently showed
that signatures of confinement and entanglement spreading can be observed on
the IBM quantum computer Vovrosh and Knolle (2020). With our new error
mitigation protocol we are able to measure the meson masses of confinement
induced bound states directly on the IBM device and can quantitatively measure
entanglement spreading previously out of reach. Furthermore, we corroborate
these results with classical simulations of error models taken from an IBMQ
device in a more general setting to show the wide applicability of our
mitigation scheme. Finally, we close with a discussion and outline future
applications.

##  II Error Mitigation Protocol

Error channels can be conveniently modelled through Kraus operators Nielsen
and Chuang (2002) which define a completely positive map on the density matrix

$\rho\rightarrow\sum_{i}K_{i}\rho K^{\dagger}_{i}\;\;\text{ such that
}\;\;\sum_{i}K^{\dagger}_{i}K_{i}=1.$ (1)

For a single qubit, choosing four Kraus operators to be proportional to the
Pauli operators (identity included) defines a depolarizing channel. The
proportionality constants are related to the individual error probabilities
and must satisfy the identity constraint to preserve the trace of $\rho$.

Here, we concentrate on such depolarizing errors which are inevitably present
in any digital quantum simulator platform and can be treated without much
specific knowledge about the device performance (which also fluctuates over
time). Moreover, our focus on depolarizing errors has the advantage of being
treatable mathematically in a controlled way as detailed below. In fact, given
a deep enough quantum circuit (with sufficient qubits) for self-averaging of
incoherent errors to occur, this depolarizing error model is a good
approximation to the physical errors in a device (more details are given in
the supplementary material). Last but not least, the large improvement in the
quality of results on the IBM device motivate our choice of depolarizing
errors a posteriori.

An n-single qubit depolarizing error channel can be modelled via

$\mathcal{E}^{\otimes
n}(\rho)=(1-p)^{n}\rho+\sum_{\alpha\in[x,y,z]}\sum_{j=1}^{n}(1-p)^{n-1}\frac{p}{3}\sigma_{\alpha}^{j}\rho\sigma_{\alpha}^{j}+...$
(2)

where $\mathcal{E}$ is the error channel, $p$ is the probability of an error
occurring for each qubit (assumed to be equal for each qubit, and for each
Pauli error) and ‘…’ indicates higher-order terms corresponding to errors on
multiple qubits Nielsen and Chuang (2002). One important feature of this
mathematical formulation is that the second term, which describes the
depolarizing error, commutes with any unitary operator. Consequently, the
error on the $i^{th}$ qubit in a quantum circuit with purely depolarizing
errors is

$\mathcal{E}^{i}(\rho)=(1-p_{i})\rho+p\Tr_{i}[\rho]\otimes\frac{\mathcal{I}_{i}}{2}$
(3)

in which $\rho$ is the density matrix, $\mathcal{I}_{i}/2$ is the maximally
mixed (i.e. completely depolarized) state for the $i^{th}$ qubit, $\Tr_{i}$ is
the partial trace over the $i^{th}$ qubit and $p_{i}$ is the error on the
$i^{th}$ qubit.

Instead of dealing with all combinations of single qubit errors, we
approximate the total error channel of Eq. (2), under the assumption of
symmetric depolarization Eq. (3), as an effective depolarizing channel on the
_entire_ quantum state

$\rho=(1-p_{\mathrm{tot}})\rho_{\mathrm{exact}}+p_{\mathrm{tot}}\frac{\mathcal{I}^{\otimes
n}}{2^{n}}$ (4)

where the effective total error probability is $p_{\mathrm{tot}}$. In
principle $p_{tot}$ is well approximated by $\prod_{i}(1-p_{i})$, however, we
do not make that identification here. We will soon show that
$p_{\mathrm{tot}}$ can be measured directly on the device. The simplicity of
this ansatz allows it to be easily calculated. It has already been shown to be
useful when mitigating measurement error Vermersch  _et al._ (2018), however
we later demonstrate that it is powerful tool for mitigating global
depolarizing errors which do not themselves originate from local depolarizing
errors. Furthermore, we later numerically demonstrate that this ansatz does
not rely on the assumption of single- and two-qubit depolarizing errors.

The many partial traces over single qubits, which would have conserved some
coherence in the remaining qubits, have been replaced by the maximally mixed
state $\mathcal{I}^{\otimes n}/2^{n}$ over the global quantum state,
destroying all coherence. We stress that even though this may not be a good
approximation for a single layer of qubit errors, it becomes a reasonable
approximation for the error channel of a many layer, many qubit circuit. Eq.
(4) is our basic ansatz for an effective error model after a many-layered
unitary circuit and $\rho_{\mathrm{exact}}$ is the exact density operator
without noise (see Supplementary Material for more details). With this ansatz
one can analytically calculate the effect of errors on a measured observable,
$\hat{O}$, via

$\langle\hat{O}\rangle=(1-p_{\mathrm{tot}})\Tr[\hat{O}\rho_{\mathrm{exact}}]+\frac{p_{\mathrm{tot}}}{2^{n}}\Tr[\hat{O}\big{]}.$
(5)

Our approach estimates $p_{tot}$ in order to apply error mitigation. We
propose and test two approaches for finding $p_{tot}$, the first based on
estimating the purity of the final state and the second based on studying
specific observables. To estimate the purity, we employ the recent protocol
for obtaining the trace of the reduced density matrix squared
$\Tr[\rho_{A}^{2}]$ via randomised measurements van Enk and Beenakker (2012);
Elben  _et al._ (2018), where $A$ is a subspace of the full density matrix.
This randomised measurement scheme has been successfully implemented in
trapped ion quantum simulators Brydges  _et al._ (2019) and recently by some
of us on the IBM quantum computer Vovrosh and Knolle (2020). We stress the
present mitigation is a far simpler task compared to inverting the quantum
error channel to tomographically reconstruct the error free quantum state
$\rho$.

As current quantum devices initialise systems in pure states that are then
manipulated with unitary transformations, $\Tr[\rho^{2}]$ over the full
Hilbert space should lead to a result that is identically one. However, since
the noisy implementation of quantum circuits will in general deviate from
unitarity, after a given quantum circuit is run on a quantum processor this
will generally not be the case. Instead, with Eq.(4) we expect that

$\Tr[\rho^{2}]=(1-p_{\mathrm{tot}})^{2}+\frac{p_{\mathrm{tot}}(1-p_{\mathrm{tot}})}{2^{n-1}}+\frac{p^{2}_{tot}}{2^{n}}$
(6)

Using the fact that $\rho_{\mathrm{exact}}$ is pure. Now, given that the left
hand side, $\Tr[\rho^{2}]$, can be measured directly on the device Vovrosh and
Knolle (2020); Brydges  _et al._ (2019), this quadratic equation can be solved
to obtain the total error $p_{\mathrm{tot}}$ Vermersch  _et al._ (2018); Elben
_et al._ (2020). We stress that correlated errors in quantum circuits do not
increase entropy and thus $p_{\mathrm{tot}}$ obtained via this method should
really be understood as the global depolarizing error probability. While this
mitigation approach does not address coherent errors, it could additionally be
combined with other techniques such as twirling Bennett  _et al._ (1996a, b);
Knill  _et al._ (2008); Emerson  _et al._ (2007).

Therefore, with $p_{\mathrm{tot}}$ extracted and $\langle\hat{O}\rangle$
measured the only unknown quantity in Eq. (5) is the desired error-free
observable $\Tr[\hat{O}\rho_{\mathrm{exact}}]$. Note, we assume that
$\Tr[\hat{O}]$ can be calculated which for most practical cases should be the
case, e.g. see our application examples below.

Putting all steps together, we finally obtain our general error mitigation
protocol:

  1. 1.

Prepare the quantum state of interest by running a quantum circuit and measure
$\Tr(\rho^{2})$, e.g. via randomized measurements Elben  _et al._ (2018);
Brydges  _et al._ (2019); Vovrosh and Knolle (2020).

  2. 2.

Use the results to obtain values for $p_{\mathrm{tot}}$ via Eq. (6).

  3. 3.

Prepare the quantum state again and measure the desired observable
$\langle\hat{O}\rangle$. 111As errors on the quantum devices can vary over
time (see Ref. Smith  _et al._ (2019a)), the determination of
$p_{\mathrm{tot}}$ has to be recalibrated accordingly.

  4. 4.

Use Eq. (5) with the measured value of $p_{\mathrm{tot}}$ to obtain the
desired $\Tr[\hat{O}\rho_{\mathrm{exact}}]$.

An alternative approach to estimating $p_{tot}$ is to consider specific
observables whose expectation values are known. For example, consider the time
dynamics of a system in which our quantum circuit approximates the time
evolution operator $U(t)=\exp{-iHt}$ where we wish to measure $\langle
O(t)\rangle$ for a range of times. We can tune the circuit such that
${t\ast(E_{max})=\epsilon<<1}$ where $E_{max}$ is the largest eigenvalue of
the Hamiltonian (shifted so $E_{min}=0$), so that our quantum circuit now
approximates the identity operation. Assuming that $\langle O(t=0)\rangle$ is
known, we can use the measurements from the quantum device and Eq.(5) to
obtain $p_{tot}$. Note, for infinite dimensional Hamiltonian with unbounded
eigenvalues, the approximation is slightly more subtle, but does not apply to
qubits. We show that, while this method is more efficient, it does not
discriminate between coherent and incoherent errors in the purity.

Our protocol can be applied to essentially any quantum circuit and quantum
simulation device. In the following, we choose a representative example from
condensed matter physics as a first application. We showcase the effectiveness
of our technique by presenting previously unobtainable results for confinement
and entanglement dynamics in spin chains.

##  III Application to spin chain confinement

An ideal testing ground for NISQ devices is that of quench dynamics in
spin-$\frac{1}{2}$ systems. A global quantum quench is a sudden change to the
systems Hamiltonian, which induces non-equilibrium dynamics. Already one
dimensional spin chains can show a wide variety of physical phenomena of
interest, for example confinement of domain wall excitations Kormos  _et al._
(2017); Liu  _et al._ (2018); Vovrosh and Knolle (2020), quantum many-body
scars Turner  _et al._ (2018); Moudgalya  _et al._ (2018); Ho  _et al._
(2019); Mark  _et al._ (2020); Shibata  _et al._ (2020), or novel fracton
excitations Nandkishore and Hermele (2019); Pai and Pretko (2020); Sala  _et
al._ (2020). All of these show up in the time evolution which is challenging
to simulate on classical computers as the Hilbert space grows exponentially
$2^{n}$ with the number of spins $n$. As spin-$\frac{1}{2}$ systems directly
map onto physical qubits, quantum computers are ideally suited for studying
the rich physics of spin chains. Recently, first digital quantum simulation
results have been reported Cervera-Lierta (2018); Zhukov  _et al._ (2018);
Francis  _et al._ (2019); Smith  _et al._ (2019a, b); Vovrosh and Knolle
(2020) but in order to obtain results which are out of reach by classical
simulations and to probe non-trivial quantum many body physics better error
mitigation techniques are needed.

We concentrate on the one-dimensional TFIM with an additional longitudinal
field given by the following Hamiltonian

$H=-J\bigg{[}\sum_{i}\sigma_{i}^{z}\sigma_{i+1}^{z}+h_{x}\sum_{i}\sigma_{i}^{x}+h_{z}\sum_{i}\sigma^{z}\bigg{]},$
(7)

where $J$ is the Ising exchange of nearest neighbour spins $\sigma_{i}$ and
$h_{x/z}$ are the relative strengths of the transverse and longitudinal fields
respectively. For $h_{z}=0$ the TFIM can be solved exactly via Jordan-Wigner
transformation and its fermionic excitations are related to free domain wall
motion. When turning on the longitudinal field, $h_{z}\neq 0$, a confining
potential between these fermions is introduced. The attraction between
fermions grows linearly with their separation, reminiscent of quark
confinement in QCD. The result is the formation of ‘mesonic’ bound states of
domain wall excitations.

Figure 2: Results from the IBM machine before and after mitigation. Quench
dynamics of local magnetization from Trotterized time evolution of the TFIM
with longitudinal field before and after error mitigation. Here, $J=1$,
$h_{x}=0.5$, $h_{x}=0.75$ and $n=7$. Data is shown for, $N_{t}$, the number of
trotter steps $N_{t}=5,6$. Results go from qualitative agreement to
quantitative agreement. After fitting a cosine function (dashed green) to the
mitigated data (green dots) the dominant frequency is clearly captured by the
IBM device, for details see Supplementary Material. Note, more details of the
circuit composition of the evolution operator can be found in Vovrosh and
Knolle (2020)

The confining potential between fermions induces non-ergodic behaviour James
_et al._ (2019), which manifests itself through persistent oscillations in the
local magnetization, $\langle\sigma_{i}^{\alpha}\rangle$
($\alpha\in\\{x,y,z\\}$) and a slowing down of the entanglement spreading.

###  III.1 Example 1: Measuring Meson Masses

The frequencies of the oscillating magnetization can be mapped directly to the
energy of the domain wall bound states Liu  _et al._ (2018), which have a
large overlap with chosen initial states. These so called meson masses are
defined as the energy difference between the lowest excited states and the
ground state.

Recent work using a trapped ion quantum simulator to simulate the long-range
TFIM model showed how, by choosing a variety of initial states, the meson
masses can be measured through the persistent oscillations of local
magnetization Tan  _et al._ (2021). Previous attempts to perform a similar
measurement on a digital quantum computer have failed for the short ranged
TFIM, Eq. (7), because the results are too noisy to resolve the smaller
amplitude of oscillations Vovrosh and Knolle (2020). As our main result, we
show that our new error mitigation enables us to obtain the meson masses from
the IBM device.

We are mainly interested in the time dependence of the local magnetization,
which further simplifies with $\Tr[\sigma^{\alpha}]=0$ in Eq. (5) to

$\langle\sigma_{i}^{\alpha}\rangle=(1-p_{\mathrm{tot}})\langle\sigma_{i}^{\alpha}\rangle_{\mathrm{exact}}.$
(8)

The time dependence can be calculated by applying a quantum circuit from a
Trotterisation of the time evolution operator (for details of the
implementation and quantum circuits see Refs. Smith  _et al._ (2019a); Vovrosh
and Knolle (2020)). Increasing the number of trotter steps, $N_{t}$, leads to
a deeper circuit and the ensuing increase in errors results in a peculiar
dampening of the magnetization dynamics which can be removed via our error
mitigation. In Fig. 2 we display the results of the local magnetization
dynamics. Our mitigation technique not only allows us to measure the main
oscillation frequency on the IBM device but in fact to do so with quantitative
agreement with the analytical results (the latter is explained in Ref. Vovrosh
and Knolle (2020)).

Figure 3: Results from the IBM device Toronto tor  of meson masses. (a,b) The
quench dynamics of the $z$-axis local magnetization is shown for different
initial states. Here, $J=1$, $h_{x}=0.75$, $h_{x}=0.75$ and $n=7$. (a) Results
calculated via exact diagonalizion. (b) Mitigated results from the IBM device.
Here, clear dominant oscillations are extracted that quantitatively agree with
the analytically derived values. (c) An illustration of how the masses are
extracted from the analytically derived energy levels Vovrosh and Knolle
(2020). (d) A comparison of the masses obtained from the IBM device and the
analytically derived values for varying $h_{z}$ showing this quantitative
agreement. Note, more details of the circuit composition of the evolution
operator can be found in Vovrosh and Knolle (2020)

We note that an additional simplification can be used for the local
magnetization for added efficiency — no additional measurement of
$\Tr[\rho^{2}]$ was required because the result of $\langle\sigma^{z}\rangle$
is known at $tJ=0$ and thus $p_{\mathrm{tot}}$ can be inferred from the
corresponding measurement on the IBM device after running the time evolution
circuit for $tJ\approx 0$. This big simplification avoids the costly
randomized measurement scheme and should be generally applicable for
observables whose value is known at $tJ=0$. However, as this does not
distinguish correlated noise from the uncorrelated noise assumed in our
ansatz, it is possible for over-estimations in the mitigation to occur,
resulting in non-physical results. An example of this is seen in Fig. 2 at
time $tJ=0.1$ in which the local magnetization is measured to be greater than
unity.

We have obtained a range of results with different number of trotter steps.
Under the assumption that $(1-p_{\mathrm{tot}})$ scales as
$(1-p_{T})^{N_{T}}$, where $p_{T}$ is the error in one trotter step, we can
extrapolate the results back to the error-free case, see green data points in
Fig. 2. We are then in a position to suppress the noise to a level which
enables us to extract different meson masses on the IBM device by a basic fit
of the main oscillation frequency. In Fig. 3 we show the data for the first
three masses obtainable by starting from different initial states, see insets.
Here, 7 spins are mapped onto 5 qubits resulting in a circuit with $5N_{T}+5$
single qubit gates and $8N_{T}$ CNOT gates with $N_{T}=5,6$. Remarkably, we
find quantitative agreement between the mitigated results and the theoretical
predictions for the scaling of the meson masses with the transverse and
longitudinal fields Vovrosh and Knolle (2020).

###  III.2 Example 2: Entanglement Spreading

As a second example of we study the suppression of half chain entanglement
entropy spreading after the spin chain quench. For $h_{z}=0$ the entanglement
entropy is expected to increase linearly from the ballistic spreading of free
fermionic excitations Fagotti and Calabrese (2008). However, with a non-zero
longitudinal field this growth is suppressed in a characteristic fashion due
to confinement Kormos  _et al._ (2017). To observe this we have implemented
the randomized measurement protocol of Ref. Brydges  _et al._ (2019) on the
IBM device for measuring the second order Réyni entanglement entropy. In Ref.
Vovrosh and Knolle (2020) two of us obtained qualitative agreement for
entanglement dynamics of six spins compared to an exact diagonalization (ED)
calculation, but for a quantitative agreement a large shift of the results was
needed.

By using the ansatz in Eq. (4) we can see that the effect of global
depolarizing errors on measurements of $\Tr[\rho_{A}^{2}]$, where $A$ is a
subsystem in consideration, is

$\displaystyle\Tr[\rho_{A}^{2}]=$
$\displaystyle(1-p_{\mathrm{tot}})^{2}\Tr[\rho_{A,exact}^{2}]$ (9)
$\displaystyle+\frac{p_{\mathrm{tot}}(1-p_{\mathrm{tot}})}{2^{n_{A}-1}}+\frac{p_{tot^{2}}}{2^{n_{A}}}.$

If $p_{\mathrm{tot}}$ is known, $\Tr[\rho_{A,exact}^{2}]$ can be extracted and
the second order Rényi entropy measurement is calculated via

$S^{(2)}(\rho_{\mathrm{exact}})=-\log_{2}(\Tr[\rho_{A,exact}^{2}]).$ (10)

Fig. 1 shows how this mitigation protocol eliminates the error in the second
order Rényi entropy results. With our error mitigation protocol and Eq. (9) we
now obtain quantitative agreement with ED results for six spins and can
account for the large shift of the results.

##  IV Applications to more general examples

To demonstrate that our method works more generally, we simulate example tasks
of measuring the expectation value of different operators. We choose a
layered, brickwork circuit consisting of a layer of single qubit rotations
followed by a layer of entangling CNOTs — the same circuit structure as for
Trotterized evolution (see Supplementary material). However, instead of time
parameterized Trotter evolution, we simply choose random single qubit angles,
corresponding to different states, and therefore different expectation values.

In order to obtain a larger amount of data to support our mitigation
technique, we numerically simulate the layered circuits with noisy gates.
Qiskit’s circuit simulator allows one to simulate arbitrary gate-based noise
models by specifying independent Kraus operators for each single- and two-
qubit gates. Using this noisy simulator, we can directly test how gate-level
non-depolarizing errors can result in an effective global circuit-level
depolarizing error, and how well our mitigation technique works in the
presence of local non-depolarizing errors. To simulate a realistic noise
model, we take the Kraus operators directly from the ibmq_santiago backend san
. This error model goes beyond the single gate depolarizing assumption, by
including (asymmetric) thermal relaxation.

In order to test local, non-local, single, and many Pauli-string operators, we
have chosen the following operators: single $Z$ operator (local, single term),
TFIM Hamiltonian (local, many terms), Random Pauli string (non-local, single
term), and the molecular Hamiltonian of a H4 (non-local, many terms). Here
many ’terms refers’ to the number of non-commuting Pauli strings that must be
measured to construct the expectation value, and locality refers to only
containing single and two-qubit Pauli strings. Note the four hydrogen
Hamiltonian (see Supplementary for details) is chosen as a simple molecular
test case with a large Hilbert space, and where we can ignore complications
from freezing out orbitals.

Figure 4 shows unmitigated and mitigated expectation values for each operator,
and clearly demonstrates the added value of purity measurements (see
Supplementary for simulation details). For expectation value mitigation (i.e.
using a single expectation value to calibrate $p_{tot}$), a single extra
circuit is needed. For trace mitigation (using Eq. (9) to calibrate $p_{tot}$)
uses 500 extra circuits to find $p_{tot}$; these are split into five groups to
reduce the bias. Once $p_{tot}$ is calibrated for a given circuit depth, this
single value is used for error mitigation of the same circuit structure, but
for different single qubit parameter angles.

As the main results we find that our global depolarising assumption is an
effective error channel for the whole circuit, and, crucially, it does not
require each single- and two-qubit errors to be themselves purely
depolarising.

Figure 4: Error mitigation simulations with non-depolarizing errors.
Unmitigated (black diamonds), trace mitigated (green circles) and single
operator mitigated (red circles) expectation values for different operators,
with the solid blue line showing perfect mitigation. The numerical simulation
uses the general error model taken from ibmq_santiago. Expectation values are
evaluated with respect to parameterized brick-like circuits with (CNOT) depth
18 (see supp for details), and thirty random parameter points per plot.
Operators are (a) Hamiltonian of 4 linear-equidistant hydrogen atoms (6
qubits), (b) Single $Z$ Pauli term (8 qubits), (c) TFIM Hamiltonian (6 qubits)
and (d) Single $XYZIXX$ string (6 qubits). These operators are chosen to span
local/global Pauli strings, and contain single/many non-commuting
decompositions.

##  V Discussion

In this work we have proposed an error mitigation technique which is simple to
implement but which retains the mathematical rigour of more complicated
techniques. Our protocol is directly applicable to any quantum simulation
whose measurements are basic expectation values, as exemplified by our results
for the time evolution of the local magnetization in spin chain dynamics.
Though our error mitigation method may not replace full state tomography, we
expect it to be useful in measuring more complicated quantities beyond simple
observables, as corroborated by our results for the entanglement entropy. An
interesting avenue of future research will be an application to variational
quantum eigensolver (VQE) problems Kandala  _et al._ (2017) and other quantum
circuits.

The computational overhead of implementing our protocol is dominated by the
cost of evaluating $\Tr(\rho^{2})$ on the quantum device. The number of
corresponding measurements is $N_{u}N_{m}$ where $N_{u}$ is the number of
randomised unitaries and $N_{m}$ is the number of random measurements. Within
the randomised scheme of Ref. Elben  _et al._ (2018) it grows exponentially
with the system size Brydges  _et al._ (2019), which poses a potential problem
for large quantum computers but is easily feasible for currently available
NISQ devices. Alas, in some cases the costly randomized measurement scheme can
be avoided entirely, i.e. in our benchmark example of spin chain dynamics a
single measurement of a local observable whose exact value is known was
sufficient. This should be true generally for quantum circuits which can be
tuned to be close to the identity. However, this scalable simplification for
obtaining $p_{\mathrm{tot}}$ potentially faces the problem that it includes
correlated errors which makes the randomised measurements preferable as long
as it is feasible.

The fact that our basic assumption of global depolarizing errors leads to such
large improvements in results is in itself remarkable. The basic conclusion is
that the total error on the IBM device is close to a global depolarizing error
at least for our choice of problems. We note our error ansatz, Eq. (4) is an
approximation even for single qubit depolarizing errors, let alone the more
complex channels that are no-doubt present in physical devices. Nevertheless
the ansatz works remarkably well for correcting physical errors of large/deep
circuits and is a no-lose addition to quantum simulation protocols. We suggest
this is because the depolarizing channel is a good approximation _per gate_ in
a many gate quantum circuit, even if this approximation breaks down for single
layers of gates. Of course, our mitigation falls short of accounting for large
coherent or correlated errors and it will be a worthwhile endeavour to think
about a controlled extension of our basic ansatz Eq. (4) for the density
matrix of a NISQ device. In particular, whether we can extend it to
incorporate other aspects of error channels Cai (2020).

As an application of our error mitigation we have presented previously
unobtainable quantitative results for confinement and entanglement dynamics of
a quantum spin chain. We have been able to extract the first meson masses of
confinement induced bound states and observed the corresponding halting of
entanglement spreading. In that context, an ambitious next step would be to
use the error mitigation to extend times that can be simulated on a quantum
computer to probe the confinement induced slow-thermalisation Liu  _et al._
(2018) or to see meson scattering events that have recently been predicted
Surace and Lerose (2020); Karpov  _et al._ (2020); Milsted  _et al._ (2020).

In general, we expect that our simple error mitigation brings us closer to
exploiting the quantum advantage of available NISQ devices for real world
practical applications.

###### Acknowledgements.

Acknowledgements.— We are grateful for discussions with Hongzheng Zhao, Adam
Smith and Peter Haynes. We acknowledge the Samsung Advanced Institute of
Technology Global Research Partnership, travel support via the Imperial-TUM
flagship partnership and the use of IBM Quantum services for this work. The
views expressed are those of the authors, and do not reflect the official
policy or position of IBM or the IBM Quantum team. This work is also supported
by the UK Hub in Quantum Computing and Simulation, part of the UK National
Quantum Technologies Programme with funding from UKRI EPSRC grant
EP/T001062/1.

## References

  * Preskill (2018) J. Preskill, Quantum 2, 79 (2018).
  * Boixo  _et al._ (2018) S. Boixo, S. V. Isakov, V. N. Smelyanskiy, R. Babbush, N. Ding, Z. Jiang, M. J. Bremner, J. M. Martinis, and H. Neven, Nature Physics 14, 595 (2018).
  * Arute  _et al._ (2019) F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandao, D. A. Buell, _et al._ , Nature 574, 505 (2019).
  * Steane (1996) A. M. Steane, Physical Review Letters 77, 793 (1996).
  * Shor (1995) P. W. Shor, Physical review A 52, R2493 (1995).
  * (6) “ibmq_paris v1.6.23, ibm quantum team. retrieved from https://quantum-computing.ibm.com, (2020),” .
  * qis (2019) “Qiskit: An open-source framework for quantum computing,”  (2019).
  * Jattana  _et al._ (2020) M. S. Jattana, F. Jin, H. De Raedt, and K. Michielsen, Quantum Information Processing 19, 1 (2020).
  * Maciejewski  _et al._ (2020) F. B. Maciejewski, Z. Zimborás, and M. Oszmaniec, Quantum 4, 257 (2020).
  * Bravyi  _et al._ (2020) S. Bravyi, S. Sheldon, A. Kandala, D. C. Mckay, and J. M. Gambetta, arXiv preprint arXiv:2006.14044 (2020).
  * Chen  _et al._ (2019) Y. Chen, M. Farahzad, S. Yoo, and T.-C. Wei, Physical Review A 100, 052315 (2019).
  * McClean  _et al._ (2020a) J. R. McClean, Z. Jiang, N. C. Rubin, R. Babbush, and H. Neven, Nature Communications 11, 1 (2020a).
  * Bonet-Monroig  _et al._ (2018) X. Bonet-Monroig, R. Sagastizabal, M. Singh, and T. O’Brien, Physical Review A 98, 062339 (2018).
  * Sagastizabal  _et al._ (2019) R. Sagastizabal, X. Bonet-Monroig, M. Singh, M. A. Rol, C. Bultink, X. Fu, C. Price, V. Ostroukh, N. Muthusubramanian, A. Bruno, _et al._ , Physical Review A 100, 010302 (2019).
  * Tran  _et al._ (2020) M. C. Tran, Y. Su, D. Carney, and J. M. Taylor, arXiv preprint arXiv:2006.16248 (2020).
  * Huggins  _et al._ (2019) W. J. Huggins, J. McClean, N. Rubin, Z. Jiang, N. Wiebe, K. B. Whaley, and R. Babbush, arXiv preprint arXiv:1907.13117 (2019).
  * McArdle  _et al._ (2019) S. McArdle, X. Yuan, and S. Benjamin, Physical review letters 122, 180501 (2019).
  * Dumitrescu  _et al._ (2018) E. F. Dumitrescu, A. J. McCaskey, G. Hagen, G. R. Jansen, T. D. Morris, T. Papenbrock, R. C. Pooser, D. J. Dean, and P. Lougovski, Physical review letters 120, 210501 (2018).
  * McClean  _et al._ (2017) J. R. McClean, M. E. Kimchi-Schwartz, J. Carter, and W. A. De Jong, Physical Review A 95, 042308 (2017).
  * Colless  _et al._ (2018) J. I. Colless, V. V. Ramasesh, D. Dahlen, M. S. Blok, M. Kimchi-Schwartz, J. McClean, J. Carter, W. De Jong, and I. Siddiqi, Physical Review X 8, 011021 (2018).
  * 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).
  * Otten and Gray (2019a) M. Otten and S. K. Gray, npj Quantum Information 5, 1 (2019a).
  * Otten and Gray (2019b) M. Otten and S. K. Gray, Physical Review A 99, 012338 (2019b).
  * Smith  _et al._ (2019a) A. Smith, M. Kim, F. Pollmann, and J. Knolle, npj Quantum Information 5, 1 (2019a).
  * Strikis  _et al._ (2020) A. Strikis, D. Qin, Y. Chen, S. C. Benjamin, and Y. Li, arXiv preprint arXiv:2005.07601 (2020).
  * Zlokapa and Gheorghiu (2020) A. Zlokapa and A. Gheorghiu, arXiv preprint arXiv:2005.10811 (2020).
  * Czarnik  _et al._ (2020) P. Czarnik, A. Arrasmith, P. J. Coles, and L. Cincio, arXiv preprint arXiv:2005.10189 (2020).
  * Li and Benjamin (2017) Y. Li and S. C. Benjamin, Physical Review X 7, 021050 (2017).
  * Endo  _et al._ (2018) S. Endo, S. C. Benjamin, and Y. Li, Physical Review X 8, 031027 (2018).
  * Temme  _et al._ (2017) K. Temme, S. Bravyi, and J. M. Gambetta, Physical review letters 119, 180509 (2017).
  * Gustafson  _et al._ (2019) E. Gustafson, Y. Meurice, and J. Unmuth-Yockey, Physical Review D 99, 094503 (2019).
  * Giurgica-Tiron  _et al._ (2020) T. Giurgica-Tiron, Y. Hindy, R. LaRose, A. Mari, and W. J. Zeng, arXiv preprint arXiv:2005.10921 (2020).
  * Cai (2020) Z. Cai, arXiv preprint arXiv:2007.01265 (2020).
  * He  _et al._ (2020) A. He, B. Nachman, W. A. de Jong, and C. W. Bauer, arXiv preprint arXiv:2003.04941 (2020).
  * Vovrosh and Knolle (2020) J. Vovrosh and J. Knolle, arXiv:2001.03044 (2020).
  * Nielsen and Chuang (2002) M. A. Nielsen and I. Chuang, “Quantum computation and quantum information,” (2002).
  * Vermersch  _et al._ (2018) B. Vermersch, A. Elben, M. Dalmonte, J. I. Cirac, and P. Zoller, Physical Review A 97, 023604 (2018).
  * van Enk and Beenakker (2012) S. van Enk and C. Beenakker, Physical review letters 108, 110503 (2012).
  * Elben  _et al._ (2018) A. Elben, B. Vermersch, M. Dalmonte, J. I. Cirac, and P. Zoller, Physical review letters 120, 050406 (2018).
  * Brydges  _et al._ (2019) T. Brydges, A. Elben, P. Jurcevic, B. Vermersch, C. Maier, B. P. Lanyon, P. Zoller, R. Blatt, and C. F. Roos, Science 364, 260 (2019).
  * Elben  _et al._ (2020) A. Elben, B. Vermersch, R. van Bijnen, C. Kokail, T. Brydges, C. Maier, M. K. Joshi, R. Blatt, C. F. Roos, and P. Zoller, Physical Review Letters 124, 010504 (2020).
  * Bennett  _et al._ (1996a) C. H. Bennett, D. P. DiVincenzo, J. A. Smolin, and W. K. Wootters, Physical Review A 54, 3824 (1996a).
  * Bennett  _et al._ (1996b) C. H. Bennett, G. Brassard, S. Popescu, B. Schumacher, J. A. Smolin, and W. K. Wootters, Physical review letters 76, 722 (1996b).
  * Knill  _et al._ (2008) E. Knill, D. Leibfried, R. Reichle, J. Britton, R. B. Blakestad, J. D. Jost, C. Langer, R. Ozeri, S. Seidelin, and D. J. Wineland, Physical Review A 77, 012307 (2008).
  * Emerson  _et al._ (2007) J. Emerson, M. Silva, O. Moussa, C. Ryan, M. Laforest, J. Baugh, D. G. Cory, and R. Laflamme, Science 317, 1893 (2007).
  * Note (1) As errors on the quantum devices can vary over time (see Ref. Smith  _et al._ (2019a)), the determination of $p_{\mathrm{tot}}$ has to be recalibrated accordingly.
  * Kormos  _et al._ (2017) M. Kormos, M. Collura, G. Takács, and P. Calabrese, Nature Physics 13, 246 (2017).
  * Liu  _et al._ (2018) F. Liu, R. Lundgren, P. Titum, G. Pagano, J. Zhang, C. Monroe, and A. V. Gorshkov, arXiv preprint arXiv:1810.02365 (2018).
  * Turner  _et al._ (2018) C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić, Nature Physics 14, 745 (2018).
  * Moudgalya  _et al._ (2018) S. Moudgalya, S. Rachel, B. A. Bernevig, and N. Regnault, Physical Review B 98, 235155 (2018).
  * Ho  _et al._ (2019) W. W. Ho, S. Choi, H. Pichler, and M. D. Lukin, Physical review letters 122, 040603 (2019).
  * Mark  _et al._ (2020) D. K. Mark, C.-J. Lin, and O. I. Motrunich, Physical Review B 101, 195131 (2020).
  * Shibata  _et al._ (2020) N. Shibata, N. Yoshioka, and H. Katsura, Physical Review Letters 124, 180604 (2020).
  * Nandkishore and Hermele (2019) R. M. Nandkishore and M. Hermele, Annual Review of Condensed Matter Physics 10, 295 (2019).
  * Pai and Pretko (2020) S. Pai and M. Pretko, Physical Review Research 2, 013094 (2020).
  * Sala  _et al._ (2020) P. Sala, T. Rakovszky, R. Verresen, M. Knap, and F. Pollmann, Physical Review X 10, 011047 (2020).
  * Cervera-Lierta (2018) A. Cervera-Lierta, Quantum 2, 114 (2018).
  * Zhukov  _et al._ (2018) A. Zhukov, S. Remizov, W. Pogosov, and Y. E. Lozovik, Quantum Information Processing 17, 223 (2018).
  * Francis  _et al._ (2019) A. Francis, J. Freericks, and A. Kemper, arXiv preprint arXiv:1909.05701 (2019).
  * Smith  _et al._ (2019b) A. Smith, B. Jobst, A. G. Green, and F. Pollmann, arXiv preprint arXiv:1910.05351 (2019b).
  * James  _et al._ (2019) A. J. James, R. M. Konik, and N. J. Robinson, Physical review letters 122, 130603 (2019).
  * Tan  _et al._ (2021) W. Tan, P. Becker, F. Liu, G. Pagano, K. Collins, A. De, L. Feng, H. Kaplan, A. Kyprianidis, R. Lundgren, _et al._ , Nature Physics , 1 (2021).
  * (63) “ibmq_toronto v1.1.6, ibm quantum team. retrieved from https://quantum-computing.ibm.com, (2020),” .
  * Fagotti and Calabrese (2008) M. Fagotti and P. Calabrese, Physical Review A 78, 010306 (2008).
  * (65) “ibmq_santiago v1.3.19, ibm quantum team. retrieved from https://quantum-computing.ibm.com, (2021),” .
  * 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).
  * Surace and Lerose (2020) F. M. Surace and A. Lerose, arXiv:2011.10583 (2020).
  * Karpov  _et al._ (2020) P. I. Karpov, G.-Y. Zhu, M. P. Heller, and M. Heyl1, arXiv:2011.11624 (2020).
  * Milsted  _et al._ (2020) A. Milsted, J. Liu, J. Preskill, and G. Vidal, “Collisions of false-vacuum bubble walls in a quantum spin chain,” (2020), arXiv:2012.07243 [quant-ph] .
  * Note (2) Assuming many qubits, but not necessarily all qubits have gates applied.
  * McClean  _et al._ (2020b) J. R. McClean, N. C. Rubin, K. J. Sung, I. D. Kivlichan, X. Bonet-Monroig, Y. Cao, C. Dai, E. S. Fried, C. Gidney, B. Gimby, P. Gokhale, T. Häner, T. Hardikar, V. Havlíček, O. Higgott, C. Huang, J. Izaac, Z. Jiang, X. Liu, S. McArdle, M. Neeley, T. O’Brien, B. O’Gorman, I. Ozfidan, M. D. Radin, J. Romero, N. P. D. Sawaya, B. Senjean, K. Setia, S. Sim, D. S. Steiger, M. Steudtner, Q. Sun, W. Sun, D. Wang, F. Zhang, and R. Babbush, Quantum Sci. Technol. 5, 034014 (2020b).
  * Bravyi  _et al._ (2017) S. Bravyi, J. M. Gambetta, A. Mezzacapo, and K. Temme, arXiv (2017), 1701.08213 .

##  VI Supplementary Information

###  VI.1 Justification of a global depolarizing Error Model

Here we justify why such a simple global depolarizing error model is
applicable. For simplicity we consider single qubit gate errors, modelled as
depolarizing errors, before moving onto multiqubit gates. The error channel
for an operator acting only on qubit $i$ is,
$\mathcal{E}_{i}[\rho]=(1-p)U_{i}\rho
U_{i}^{\dagger}+p\mbox{Tr}_{i}[\rho]\otimes\mathcal{I}_{i}/2$ where $U_{i}$
acts locally on qubit $i$. That is with probability $1-p$ the gate was
successfully applied, and with probability $p$ the qubit suffered a completely
depolarizing error, and ends up in the maximally mixed state
$\mathcal{I}/2=\mbox{diag}[1/2,1/2]$, thereby destroying quantum correlations
between qubit $i$ and the reset of the state. This error model does not
commute with entangling gates (e.g. consider the reduced state
$\mbox{Tr}_{i}[\rho]$ depends on the entanglement between qubit $i$ and the
rest of the state). While applying single gates to multiple different qubits
commute, the resulting expression quickly becomes difficult to track,
$\mathcal{E}_{i}\circ\mathcal{E}_{j}[\rho]=(1-p_{i})(1-p_{j})U_{i}U_{j}\rho
U_{j}^{\dagger}U_{i}^{\dagger}+p_{i}(1-p_{j})\mbox{Tr}_{i}[U_{j}\rho
U_{j}^{\dagger}]\otimes\mathcal{I}_{i}/2+p_{j}(1-p_{i})\mbox{Tr}_{j}[U_{i}\rho
U_{i}^{\dagger}]\otimes\mathcal{I}_{j}/2+p_{i}p_{j}\mbox{Tr}_{ij}[\rho]\otimes\mathcal{I}_{i}/2\otimes\mathcal{I}_{j}/2$.

Our first step is to approximate a layer of single qubit gates 222Assuming
many qubits, but not necessarily all qubits have gates applied as

$\displaystyle\mathcal{E}_{\mathrm{layer}}=(1-p^{\prime})U\rho
U^{\dagger}+p^{\prime}\mathcal{I}/2^{n},$ (11)

for $n$ qubits. We don’t set $(1-p^{\prime})\neq\prod_{i}(1-p_{i})$, but
rather introduce an effective $p^{\prime}$ to account for components in the
total sum that have some reasonable overlap with $U\rho U^{\dagger}$. E.g. for
each $j$, the state $\mbox{Tr}_{j}[U\rho U^{\dagger}]\otimes\mathcal{I}_{j}/2$
may have a large overlap with $U\rho U^{\dagger}$, (and to a lesser extent
$\mbox{Tr}_{ij}$ for each $i,j$), and we’d like to include such contributions
without worrying about combinatorial multiplicity, or the state-dependent
details of all the partial traces. We can follow the same steps to arrive at
Eq. (11) for entangling gates $U_{ij}$ once we note that the 2-qubit
depolarizing error channel be written as
$\mathcal{E}_{ij}[\rho]=(1-p)U_{ij}\rho
U_{ij}^{\dagger}+p\mbox{Tr}_{ij}[\rho]\otimes\mathcal{I}_{ij}/4$, where
$\mathcal{I}_{ij}/4=\mbox{diag}[1,1,1,1]/4$ is the maximally mixed two qubit
state.

We can construct the quantum channel for the full circuit. The non-unitary
part of Eq. (11) commutes with itself over different layers meaning the a full
circuit error channel is approximated as

$\displaystyle\mathcal{E}_{\mathrm{circ}}[\rho]=(1-p_{\mathrm{tot}})U_{\mathrm{circ}}\rho
U_{\mathrm{circ}}^{\dagger}+p_{\mathrm{tot}}\mathcal{I}_{2}^{\otimes n}/2^{n}$
(12)

where $(1-p_{\mathrm{tot}})=\prod_{l}(1-p_{l}^{\prime})$ (product over
effective layer $p^{\prime}$’s) holds due to commuting non-unitary terms.
While this is not equal to $\prod_{g}(1-p_{g})$ (product over gate error
probabilities), the latter is a good zeroth-order estimate of this term.

Finally we note the total map Eq. (12) unravels to an effective error model
for single (i.e. one or two qubit) gates, given by:
$\mathcal{E}_{i}[\rho]=(1-p_{i}^{\prime\prime})U_{i}\rho
U_{i}^{\dagger}+p_{i}^{\prime\prime}\mathcal{I}_{2}^{\otimes n}/2^{n}$. Note
the all-or-nothing description of the map: either the gate is implemented
perfectly, or the _entire_ quantum state is destroyed, not just part of the
state as is the case for $\mbox{Tr}_{i}[\rho]\otimes\mathcal{I}_{i}/2$ etc.
While this seems quite an unphysical model for a single gate, it is a good
model of the effective gate error per gate in a multi-gate circuit.

###  VI.2 Cosine Fitting

In order to extract the frequencies from the local magnetization results
presented in Fig. 2 and Fig. 3, a function that contained a singular cosine
was used, namely,

$\langle\sigma^{z}(t)\rangle\sim Ae^{-dt}\cos(\omega t)+c_{1}t+c_{2}.$ (13)

Here, $A$ is the amplitude of the oscillations, $D$ allows for decay of the
amplitude, $\omega$ is the oscillation frequency and $c_{1}$ and $c_{2}$ give
a time dependent shift. While this function can not capture long term
behaviour of $\langle\sigma^{z}\rangle$, it is able to extract out the
dominant oscillations for short time dynamics.

###  VI.3 Statistical Errors Introduced by the Error Mitigation

The error when calculating $p_{tot}$ from $\Tr(\rho^{2})$ is given by

$\delta
p_{tot}=\pm(1-p_{tot})\frac{\delta\Tr(\rho^{2})}{2(2^{n}-\Tr(\rho^{2}))},$
(14)

where $\delta p_{tot}$ is the error in $p_{tot}$ and $\delta\Tr(\rho^{2})$ is
the error in the measurement of $\Tr(\rho^{2})$. Note that
$\delta\Tr(\rho^{2})$ can be calculated by resampling e.g. Jackknife or
Bootstrap Brydges  _et al._ (2019); Vovrosh and Knolle (2020). Thus, when
using the proposed error mitigation on an observable $\hat{O}$ this results
additional uncertainty on top of the measurement errors. More precisely,

$\delta\langle\hat{O}\rangle=\pm\Bigg{(}\frac{\delta_{meas}\langle\hat{O}\rangle}{1-p_{tot}}+\frac{\delta
p_{tot}\big{(}\langle\hat{O}\rangle-2^{-n}\Tr(\hat{O})\big{)}}{(1-p_{tot})^{2}}\Bigg{)},$
(15)

where $\delta\langle\hat{O}\rangle$ is the error in the mitigated value of
$\langle\hat{O}\rangle$ and $\delta_{meas}\langle\hat{O}\rangle$ is the error
in the measurement of $\langle\hat{O}\rangle$. Note, as the local
magnetisation results presented in this paper bypassed the measurement of
$Tr(\rho^{2})$, the error contributing to $\delta p_{tot}$ is just shot noise.
Thus, $\delta\langle\sigma_{i}^{z}\rangle$ is of the on the order $\sim
O(10^{-3})$ and is not displayed in Fig. 2 or Fig. 3.

###  VI.4 Numerical simulations

All numerical simulations are performed using qiskit qis (2019) with the
santaiago error model (as of 12${}^{t}h$ March 2021). For the hydrogen
Hamiltonian, openfermion McClean  _et al._ (2020b) was used to compute the
orbital integrals (for a row of four H-atoms with 0.1 nm inter-atomic
distance), and then map the problem to qubits using the symmetry conserving
Bravyi-Kitaev transformation Bravyi  _et al._ (2017), finally the
$\mathds{Z}_{2}$ symmetries (corresponding to particle spin) are removed,
reducing the number of qubits by two. For the TFIM, the Hamiltonian was
reduced to two non-commuting Pauli strings from which the total Hamiltonian
can be evaluated. Finally the single $Z$, and $XYZIXX$ (chosen randomly)
operators were measured trivially. We do not directly account for State-
Preparation and Measurement (SPAM) errors, but we note that we can (and do)
include SPAM errors as contributing to $p_{tot}$ in our channel ansatz.

We numerically simulate circuits consisting of layers of parameterized _u3_
gates, generic single-qubit rotation gates with 3 Euler angles, followed by an
entangling layer. The circuit depth is defined by the number of entangling
(CNOT) layers. There is a final _u3_ layer before measurement. The
parameterized circuit now creates a state $\rho(\theta)$, and varying $\theta$
changes the state and corresponding operator expectation values. For the
numerical simulations we consider how our mitigation performs for different
depth circuits, and for the four tasks shown in the main text. For each depth
the calibration $p_{tot}$ must be done independently, as $p_{tot}$ is an
effective parameter for the whole quantum channel and cannot simple be
calculated a priori for different depths.

Calibration.— Two independent methods of calibrating $p_{tot}$ are considered.
Firstly, via single expectation value: direct simulation of a circuit with
fixed parameters, which is then compared to the error free expectation value
to find $p_{tot}$. Secondly, via trace: from directly measuring
$\mathrm{Tr}[\rho^{2}]$ and estimating $p_{tot}$ via Eq. (9).

For single expectation calibration, a single fixed parameter circuit was used
to find $\langle O\rangle$ (for each operator $O$ considered). This result was
compared to the error free simulation of the exact same circuit to estimate
$p_{tot}$. To evaluate these expectation values 40960 shots were used for each
non-commuting term in the decomposition of the given operator into Pauli
strings. For trace mitigation, an extra 500 circuits are used to estimate
$p_{tot}$ using 8196 shots per circuit. This was broken up into five lots of
100 circuits to estimate $p_{tot}$ at five different (random) circuit
parameterizations. For each parameterization, the method of Ref. Elben  _et
al._ (2018) is used to find five independent estimates of $p_{tot}$, which are
then averaged to reduce the bias in the estimate of $p_{tot}$ for the circuit.

Mitigation.— Once $p_{tot}$ is estimated either from the single value or trace
measurement, it is used to mitigate further simulated expectation values. New
expectation values are obtained by varying the parameters in the circuit,
while keeping the circuit structure fixed. For each depth and each task,
thirty different (parameterized) circuits are simulated, and mitigation
applied. Figure 5 shows the mitigation results for different circuit depths.
While the full trace-estimation is robust showing a significant improvement in
every case, the single expectation value mitigation is more variable. This is
due to the vanishing expectation value of Pauli strings for random circuits,
which in turn increases the effect of shot noise when inverting Eq. 8. We note
this shot noise problem is absent when calibrating a Pauli string (or single
operator) with a (noise free) unit expectation value, as in the main text.

Figure 5: Complete mitigation simulations with non-depolarizing errors.
Scaling of our mitigation schemes (circles) with circuit depth for four
different operator expectation values. Single operator (red circles) show some
improvement, while trace mitigated (green circles) show a robust improvement
over the unmitigated results. As in the main text, numerical simulation uses
the ibmq_santiago error model, and circuits are layers _u3_ and CNOT, with the
depth equal to the number of entangling layers. Note there are some
combinations (e.g. TFIM at depth 14 etc) where the single expectation value
mitigation found $p_{tot}=1$, and therefore did not improve results. Likewise
there are other combinations (e.g. TFIM depth 22) where the single expectation
value mitigation is omitted because it gave an unphysical estimate of
$p_{tot}<0$. Inset numbers show $p_{tot}$ as estimated for the trace
estimation.