#  _CHANDRA_ X-RAY OBSERVATIONS OF ABELL 119:
COLD FRONTS AND A SHOCK IN AN EVOLVED OFF-AXIS MERGER

Courtney B. Watson Institute for Astrophysical Research and Department of
Astronomy, Boston University, Boston, MA 02215, USA<EMAIL_ADDRESS>Elizabeth
L. Blanton Institute for Astrophysical Research and Department of Astronomy,
Boston University, Boston, MA 02215, USA  Scott W. Randall Harvard Smithsonian
Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA  Craig L.
Sarazin Department of Astronomy, University of Virginia, P.O. Box 400325,
Charlottesville, VA 22904-4325, USA  Arnab Sarkar Kavli Institute for
Astrophysics and Space Research, Massachusetts Institute of Technology, 77
Massachusetts Ave, Cambridge, MA 02139, USA  John A. ZuHone Harvard
Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138,
USA  E. M. Douglass Farmingdale State College—SUNY, 2350 Broadhollow Rd.,
Farmingdale, NY 11735, USA

(Accepted 08/07/2023)

###### Abstract

We present _Chandra_ X-ray observations of the dynamically complex galaxy
cluster Abell 119 ($z=0.044$). A119 is host to two NAT radio sources (0053-015
& 0053-016) whose tails are oriented parallel to each other despite
orthogonally oriented jet axes. Imaging and spectral analysis reveal X-ray
emission elongated along the NE-SW axis along with the presence of complex
structures, including surface brightness discontinuities, which suggest
possible merger activity along this axis. From radial profiles of the X-ray
surface brightness, temperature, pressure, and density, we identify two
surface brightness edges which are found to be cold fronts, possibly
associated with large-scale sloshing of ICM gas. We also identify a brightness
edge to the south which is found to be a shock front with Mach number
$M=1.21\pm 0.11$, consistent with a merger shock. In addition, previous
optical studies show alignment of optical substructures along the north-south
direction. The elongated X-ray emission, orientations of the NAT tails, and
alignment of optical substructure all suggest recent or on-going merger
activity in the NE-SW direction.

††journal: The Astrophysical Journal

##  1 Introduction

Galaxy clusters are the largest gravitationally bound objects in the universe,
containing hundreds to thousands of galaxies. In the current model of large-
scale structure formation, galaxy clusters grow through mergers of smaller
systems in a process called hierarchical structure formation (West et al.,
1995). Galaxy cluster mergers are some of the most energetic events in the
universe. During a merging event, this energy can be dispersed throughout the
gas of the intracluster medium (ICM), leaving observable impacts on the ICM
structure.

In merging systems, bulk motions of ICM gas can give rise to density
discontinuities observed as edges in the surface brightness in X-ray
observations of the hot ICM. Surface brightness edges typically seen in
merging systems are indicative of either shocks or cold fronts.

The formation of shock fronts resulting from cluster mergers should be fairly
common and has been seen in numerous simulations (Schindler & Muller, 1993;
Roettiger et al., 1993; Burns, 1998). However, while evidence of shocks has
been seen, i.e., regions of ICM showing shock-heating (Markevitch & Vikhlinin,
2001; Markevitch et al., 2003; Kempner & David, 2004; Govoni et al., 2004),
observation of shock fronts, with both a sharp edge in gas density and a
distinct jump in temperature, have only been seen in a handful of clusters.
The lack of direct observations of shock fronts is mostly because we have to
observe the merger at the right time (i.e. before the shock has propagated out
to the low-brightness outskirts) and from the right orientation (i.e the
tangent plane of the shock front edge must be close to perpendicular to the
plane of the sky). Despite the challenges of observing shock fronts, when we
are lucky enough to see/study them, they can provide important insight into
the merger history and dynamics of the cluster.

Cold fronts (or contact discontinuities), which are more easily observed than
shock fronts (Markevitch, 2010), have proven to be invaluable tools in
understanding the physical properties of the ICM (Markevitch & Vikhlinin,
2007) and can be used as a gauge of cluster merger activity (Owers et al.,
2009). Simulations of cluster mergers have shown that cold fronts can arise in
several different ways, depending on the merger conditions (Poole et al.,
2006; Ascasibar & Markevitch, 2006; ZuHone et al., 2010, 2016). However, cold
fronts resulting from cluster mergers can generally be classified into two
categories (Owers et al., 2009; Markevitch et al., 2000): the remnant core
(e.g. the Bullet cluster; Markevitch et al. (2002)), and sloshing spirals
(e.g. Abell 2052; Blanton et al. (2011) or Abell 98; Sarkar et al. (2023)).
The remnant core scenario results in a contact discontinuity between the
interface of the infalling subcluster’s cool core and the hotter surrounding
ICM. Cold fronts resulting from gas sloshing are typically seen in relaxed
clusters or following a minor merger event and arise from the cool gas at the
center being displaced from the central potential well, resulting in “sloshing
spirals” (Tittley & Henriksen, 2005; Markevitch et al., 2000; Roediger et al.,
2011; Paterno-Mahler et al., 2013). The observed morphology of sloshing
depends on the sloshing direction. If the sloshing direction is closer to the
plane of the sky, a spiral structure may be more evident in direct imaging.
However, a morphology of concentric arcs around the cluster center is more
likely to be observed when the sloshing direction is not in the plane of the
sky (Gastaldello et al., 2013).

The study of clusters hosting sloshing spirals is important not only for
informing the merger history of these clusters but because they also can
provide insight into the origin of non-cool core (NCC) clusters. It is
expected that cool core (CCs) systems should be host to sloshing cold fronts
Markevitch & Vikhlinin (2007), however two notable exceptions to this exist.
Abell 2142 (Rossetti et al., 2013) and Abell 1763 (Douglass et al., 2018) both
lack CCs but show evidence of off-axis merger activity through the presence of
sloshing cold fronts. The recent simulations of Valdarnini & Sarazin (2021)
also support this link between merger history and CC vs NCC systems. This
indicates off-axis mergers may play a more meaningful role than previously
assumed in the formation of NCC systems.

In addition to surface brightness edges in the hot ICM, cluster mergers can
sometimes distort radio sources into a bent-tail morphology. The radio
emission from these sources comes from outflows driven by the central active
galactic nuclei (AGN) of the host galaxy. These outflows result in radio-
emitting jets and lobes that are oriented perpendicular to the accretion disk.
The tails are bent via ram pressure resulting from relative motion between the
host galaxy and the ICM. One type of bent-tail radio source is the narrow-
angle tail (NAT), so-called because the radio lobes are typically bent back at
large angles (relative to the outflowing jet axis). NATs are generally thought
to be formed via ram-pressure resulting from the high-speed infall or orbital
motions of the host galaxy relative to the ICM. However, Bliton et al. (1998)
suggest that, in many cases, the orbital velocity of the host galaxy is simply
not high enough to produce the high degree of bending seen in NATs. They
present an alternative scenario of NAT formation due to the host cluster
experiencing a merger event. In this scenario, the cluster merger activity
induces bulk motion of the ICM, which has been shown to produce the ram-
pressure needed to bend the radio jets to the high angles of NATs. This has
since been supported by simulations which show that high-velocity bulk flow of
the ICM can produce bending of radio tails (Roettiger et al., 1993; Mendygral
et al., 2012).

In this paper, we present _Chandra_ X-ray observations of the galaxy cluster
Abell 119 (z=0.044; A119 hereafter). A119 lacks a cool core (Laganá et al.,
2019) and previous studies have reported global X-ray temperatures ranging
from 5.1 to 5.9 keV (Edge et al., 1990; David et al., 1993; Ebeling et al.,
1996; Markevitch et al., 1998; Elkholy et al., 2015). The brightest cluster
galaxy, UGC 579, is classified as a cD galaxy (Postman & Lauer, 1995; Saglia
et al., 1997). A119 hosts two NAT sources located E (0053-015), which we will
refer to as the “eastern” NAT, and SW (0053-016), which we will refer to as
the “western” NAT, of the central cD galaxy (Feretti et al., 1999). There is a
third radio source, 3C 29, located $\sim 22$′ from the cluster center, but
lies outside of our Chandra coverage and is therefore excluded from the
analysis in this paper. Previous studies have suggested that A119 is a
dynamically complex cluster, showing evidence of multiple substructures
(Fabricant et al., 1993; Ramella et al., 2007; Tian et al., 2012; Lee et al.,
2016). This paper aims to use the presented _Chandra_ X-ray observations, in
conjunction with results from previous optical and radio studies, to form a
picture of the merger history of A119.

Throughout this paper, we adopt a standard $\Lambda$CDM cosmology with
$H_{0}=70$ km s-1 Mpc-1, $\Omega_{M}=0.3$, and $\Omega_{\Lambda}=0.7$. At
$z=0.044$, the luminosity distance is $D_{L}=194.8$ Mpc, the angular size
distance is $D_{A}=178.7$ Mpc, and $1\arcsec=0.866$ kpc. Uncertainties
reported here are 1$\sigma$ confidence intervals unless otherwise noted.

##  2 _Chandra_ Observations & Data Reduction

We used two _Chandra_ pointings taken 2003 September 4 (ObsID 4180; 11.93 ks)
and 2007 August 20 (ObsID 7918; 45.04 ks). Both observations were taken using
the ACIS-I configuration in very faint (VFAINT) mode. The cluster emission
fills the majority of the four ACIS-I CCDs, so both observations also include
data from the S2 CCD for use in background flare removal.

The data were processed using CIAO (v4.12) and CALDB (v4.9.0) provided by the
_Chandra_ X-ray Center (CXC). Level-2 events files were produced using
`chandra_repro` with `check_vf_pha=yes`. Background flares were removed by
extracting light curves using `dmextract`, and detecting and removing flares
using `deflare`. After filtering for flares, the total exposure time of the
two observations was 46.9 ks. Blank-sky background files were selected from
CALDB, using `acis_background_lookup`, and reprojected to match the
observations. Background images were exposure corrected and normalized such
that the hard band count rate in the 10-12 keV band matched that of the
observations. Aspect histograms, instrument maps, and exposure maps for each
observation were created following the step-by-step guide outlined in the CIAO
thread _Multiple Chip ACIS Exposure Map_
111https://cxc.harvard.edu/ciao/threads/expmap_acis_multi/. Exposure maps were
then filtered to exclude regions with $<10\%$ the maximum exposure. Both
ObsIDs were stacked and reprojected to match the coordinates of ObsID 7918
using `reproject_events` to create a final combined image with a cleaned
exposure time of 46.9 ks.

Figure 1: (_Left_) Gaussian smoothed (using a 7.38″ radius Gaussian) and
(_Right_) adaptively smoothed background and exposure corrected composite
X-ray images in the 0.7-7 keV energy range with 1.4 GHz VLA FIRST (green) and
150 MHz TGSS (cyan) radio contours.

$\begin{array}[]{cc}\includegraphics[width=130.08731pt]{NAT15_zoom.pdf}&\includegraphics[width=151.76964pt]{NAT16_zoom.pdf}\end{array}$

Figure 2: Zoom in of the eastern (_Left_) and western (_Right_) NAT sources on
the adaptively smoothed X-ray image. Overlaid are the 1.4 GHz VLA FIRST
(green) and 150 MHz TGSS (cyan) radio contours.

##  3 Imaging Analysis

The left panel of Fig. 1 shows the background-subtracted and exposure-
corrected composite _Chandra_ image in the 0.7-7 keV band. The image was
smoothed using a 7.38″ radius Gaussian. In addition to Gaussian smoothing, we
use the CIAO routine `csmooth`, with a minimum signal-to-noise ratio of 4, to
adaptively smooth the X-ray image. The right panel of Fig. 1 shows the
adaptively smoothed, background and exposure corrected composite X-ray image
of A119 in the 0.7-7 keV band. The overall X-ray emission is fairly asymmetric
with an elongation to the NE resulting in a “teardrop” shape. The adaptively
smoothed X-ray image shows the possible presence of clumpy substructure within
the ICM.

Overlaid on the two X-ray images are radio contours of the two NAT sources.
The green contours are taken from the Very Large Array (VLA) Faint Images of
the Radio Sky at Twenty-cm (FIRST) survey at 1.4 GHz with a beam size of
$6.4\arcsec\times 5.4\arcsec$ (Becker et al., 1995). The cyan contours are
taken from the 150 MHz all-sky survey observed with Giant Metrewave Radio
Telescope (GMRT) as part of the TIFR GMRT Sky Survey (TGSS) project with a
$26.7\arcsec\times 24.9\arcsec$ beam size (Intema et al., 2017).

Fig. 2 shows zoomed in views of the adaptively smoothed X-ray image in the
region of the two NAT sources with the 150 MHz and 1.5 GHz radio contours
overlaid, showing the jet/lobe structure in more detail. The tails of the two
NAT sources are also orientated along the NE-SW direction, following the X-ray
emission. This is notable because while the tails are oriented parallel, the
jets actually leave their host galaxies in very different directions. In the
eastern NAT, the jets are leaving the host galaxy along the NE-SW plane, while
in the western NAT the jets are leaving in more of the SE-NW direction. The
150 MHz contours of the eastern NAT show a break in the radio emission that
appears to bisect an X-ray brightness edge to the south. In addition, the
eastern NAT could be displacing the surrounding X-ray emitting gas. There is a
prominent cavity of decreased surface brightness around the NAT that can be
seen in the adaptively smoothed X-ray image in Fig. 1 and the zoomed in view
in Fig. 2. The region of decreased X-ray surface brightness is localized to
the region of higher frequency radio emission, shown by the green contours in
Fig. 2. To asses the significance of the cavity, we plot the surface
brightness in a semi-circular region over the NAT source (see Fig. 3). The
bubble of decreased emission is found to be significant at the 3.3$\sigma$
level, relative to the average background of the two adjacent bins. The
interior edge of the X-ray deficient cavity is also observed as a surface
brightness edge at a radius of $\sim 135$″ in the radial profiles across the
NAT source, which are discussed later in §6.2 and shown in Fig. 11.

Fig. 4 shows an optical image from the Digitized Sky Survey (DSS) with
contours of X-ray emission (black), taken from the adaptively smoothed X-ray
image (right panel of Fig.1), and radio (green and blue) emission of the two
NAT sources are overlaid. The magenta ‘x’ in the optical image marks the peak
of the X-ray emission and shows that the central cD galaxy is slightly offset
to the east of the peak. Point sources are identified using the CIAO routine
`wavdetect`, using wavelet scales of 4, 6, 8, and 16 pixels (1 pix = 0.492″).
The detected point sources are indicated in Fig. 4 with red ellipses. After
removing the detected point sources, we use the CIAO routine `dmfilth` to
replace the removed point sources pixels with interpolated pixel values from
the surrounding area, creating an image of the diffuse emission which is used
in some of the following analysis.

$\begin{array}[]{cc}\includegraphics[width=182.1196pt]{gausssmooth_azprof.pdf}&\includegraphics[width=195.12767pt]{azprof.pdf}\\\
\end{array}$

Figure 3: (_Left_) X-ray image showing the region used to plot the azimuthal
surface brightness profile (_right_) across the eastern NAT source to
highlight the decreased X-ray emission surrounding the NAT discussed in §3.
Figure 4: Optical image from the Digitized Sky Survey (DSS) with the 1.4 GHz
(green) and 150 MHz (blue) radio contours as well as black contours showing
the X-ray emission from the adaptively smoothed image. Red ellipses mark point
sources detected by wavdetect. The magenta ‘x’ marks the peak of the X-ray
emission.

###  3.1 Highlighting X-ray Substructure

To enhance features seen in the X-ray image, we employed two techniques: beta-
model subtraction and Gaussian Gradient Magnitude (GGM) filtering. For both
methods, we use the unsmoothed, merged, background and exposure corrected
0.7-7 keV image of diffuse emission. These techniques help to highlight faint
structure within the ICM and are briefly described below.

####  3.1.1 Beta-model Residual Image

We fit an elliptical 2D $\beta$-model using the `beta2d` model in `Sherpa` to
a Gaussian smoothed ($\sigma=2\arcsec$) background and exposure corrected
image in the 0.7-7 keV range. The peak of the X-ray emission, as determined by
the 2D $\beta$-model fit, is found to be located roughly 41 kpc NNW of the
central cD galaxy (see Fig. 4).

The fitted 2D $\beta$-model is subtracted from the unsmoothed X-ray image and
the resulting residual image, smoothed using a $\sigma=15\arcsec$ Gaussian, is
shown in Fig. 5. The smoothed residual image shows the presence of brightness
edges seen in the Gaussian smoothed and adaptively smoothed X-ray images (Fig.
1). The model subtracted residual image shows a similar tear-drop like shape
seen in the X-ray image, with excess emission extending over 580″to the NNE
from the cluster center. Additionally, excess emission towards the south
suggest the presence of an X-ray brightness edge along the southern edge of
the cluster. These features are studied in more detail in §5.

Figure 5: Excess emission in the 0.7-7 keV energy range after subtraction of a
2D $\beta$-model. The residual image was smoothed using a $\sigma=15\arcsec$
Gaussian.

####  3.1.2 Gaussian Gradient Magnitude (GGM) Filter

A Gaussian Gradient Magnitude (GGM) filter was applied using the
`gaussian_gradient_magnitude` within `SCIPY`222http://scipy.org/. GGM
filtering calculates gradients in the 2D data, assuming Gaussian derivatives,
thereby enhancing surface brightness edges (Sanders et al., 2016). The scaling
is determined through the choice of $\sigma$ (the width of the Gaussian
assumed), with smaller values typically chosen to highlight small scale
structure while larger values of $\sigma$ can be used to highlight large scale
structure. Figure 6 shows the GGM filtered image with $\sigma=40$ pixels (1
pixel = 0.492″). This scaling was chosen to highlight large-scale structures,
particularly in the cluster outer regions where brightness edges are seen in
the X-ray image.

While we perform a more detailed analysis of the significance of the
brightness edges in §5, here we highlight some of the features seen in the GGM
filtered image. Several notable features are indicated in Fig. 6. In the far
north, there is an edge along the western part of the tear-drop shape (green
arrows). North of the cluster center there is an edge that has a concave shape
(cyan arrows), similar to the Kelvin-Helmholtz instabilities that have been
seen in earlier observations (e.g., Walker et al. (2017)) and in simulations
of sloshing spirals (Roediger et al., 2013). To the south, there is an edge
(yellow arrows) which could be bisecting the southern end of the eastern NAT’s
tail. A more detailed analysis of the surface brightness edges identified in
A119 is presented in §5.

Figure 6: Gaussian gradient magnitude filtered (0.7-7 keV) image with
$\sigma=40$ pixels (1 pixel $=$ 0.492″). The brightest regions correspond to
the steepest changes in X-ray surface brightness. Overlaid are the 1.4 GHz
(green) and 150 MHz (cyan) radio contours. Some notable features are
highlighted.

##  4 Spectral Analysis

In order to determine the thermodynamic properties of the ICM in A119, we
extract the X-ray spectra from our _Chandra_ observations. Unless otherwise
noted, all spectra were extracted separately for the two observations using
`specextract` with point sources excluded. Background spectra were taken from
the hard band normalized, reprojected blank-sky background observations.
Spectra were binned such that each energy bin contained a minimum of 40
background subtracted counts. Extracted spectra were fitted simultaneously in
the 0.7-7 keV range in `XSPEC` (Arnaud, 1996) using a thermal `APEC` model.
Galactic absorption was fixed at $N_{H}=3.2\times 10^{20}$ cm-2 (Dickey,
1990). We adopted the Asplund et al. (2009) solar abundance tables. The
temperature, metallicity, and normalization were allowed to vary.

###  4.1 Total Spectrum

To determine global ICM properties we extracted the spectrum from a 400″
circular region, centered on the cluster center. This was the largest region
size that could be used before reaching a chip edge. From a total of $\sim$
87,000 background subtracted counts, the global temperature was found to be
$6.0\pm 0.1$ keV with an abundance of $0.41\pm 0.04$ Z⊙. The fit was good with
$\chi^{2}=648.6$ with 597 degrees of freedom, corresponding to a reduced
$\chi^{2}=1.1$. From the global temperature value, the average sound speed of
the cluster was calculated to be $c_{s}=\sqrt{\gamma\frac{kT}{\mu
m_{p}}}=1265\pm 11$ km/s, where $\gamma=\frac{5}{3}$, $\mu=0.6$, and $m_{p}$
is the mass of the proton.

In addition to the global fit, we created profiles of the surface brightness,
projected temperature, pressure and density for the overall cluster emission.
Concentric circular annuli were defined extending out to a radius of 400″ to
match the region defined above. Surface brightness profiles were extracted
from the merged, background subtracted, and exposure corrected X-ray image
with point sources excluded. Projected temperature and abundance profiles were
created by extracting and fitting the spectrum within each circular annulus,
using the method described in §4.

The pressure and density profiles were calculated following the methods of
Kriss et al. (1983), Wong et al. (2008) and Blanton et al. (2009). The X-ray
surface brightness profiles were deprojected to determine the X-ray emissivity
and gas density, with both assumed to be constant within spherical shells. To
determine the radial variations in pressure, we used the deprojected density
with the fitted projected temperatures, assuming the fitted temperature from
the projected spectra is the temperature at that spherical radius.

Fig. 7 shows the radial profiles of (a) surface brightness, (b) density, (c)
pressure, (d) projected temperature, and (e) projected abundance. The
temperature shows some radial variation with an overall slight decline from
the center outward. The abundance profile is roughly constant with radius,
with some variations, including a drop at approximately 200″. The average
values of the temperature and abundance from the radial profiles is consistent
with the values found above for the total spectrum fit. While the profiles of
brightness, density, and pressure show some indication of a change in slope at
$\sim$ 180″, which would correspond to the slight increase in temperature seen
in panel (d), there is no notable presence of the edges discussed in §3.
Because there is no apparent indication of edge features in the global radial
profiles, we define a number of subregions in which we examine the edge
features seen in the X-ray images in more detail (see §5).

Figure 7: Radial profiles of (a) surface brightness, (b) deprojected density,
(c) deprojected pressure, (d) projected temperature, and (e) projected
abundance for the global cluster emission.

###  4.2 Temperature Map

We created a spectral map of A119, following the methods of Randall et al.
(2008, 2009), to examine the temperature distribution of the cluster. For each
spectral map pixel, a circular region was defined and allowed to increase in
radius until a minimum of 2600 net source counts is contained. Within this
circular radius, the spectrum was extracted and fitted following the method
described above (see start of §4).

Figure 8: Temperature map of A119, overlaid with 1.4 GHz (green) and 150 MHz
(blue) radio contours. The white contours are from the smoothed X-ray image
(see Fig.1).

Figure 8 shows the projected temperature map of A119 with 1.4 GHz FIRST and
150 MHz TGSS radio contours of the two NAT sources overlaid. The white
contours are taken from the smoothed X-ray image (see Fig.1). The temperature
map is binned such that the pixel scale is 3.94”/pixel. Errors for the
temperature map range from about 5% at the cluster center to about 21% in the
outskirts.

Cool ($\sim 6$ keV) gas surrounds the central region of the cluster and is
elongated along the NE-SW direction. The cool gas extends up to the north,
connecting to a large ($\sim 168$ kpc in diameter) cool ($\sim 4$ keV) spot in
the NE corner. The central cool gas also appears to extend out to a cold
($\sim 4$ keV) spot in the SE corner. To the north, there is a strip of cool
gas corresponding to the western edge of the tear-drop shape seen in the X-ray
image. Northwest of this strip is a clump of hot ($\sim 7$ keV) gas that is
137 kpc across. Clumps of hot gas surround the cluster center at radii of
$\sim 230$ kpc. The clumpy hot gas to the north of the core has temperatures
of $\sim 6-7$ keV while the arm of clumpy hot gas to the south of the core has
temperatures of $\sim 7-8$ keV. In order to determine the significance of the
features seen in the spectral maps, follow-up detailed spectral fitting is
performed on a number of regions of interest and is discussed further in §5.

##  5 Radial Profiles

The edges seen in the X-ray and GGM images could be indicators of cold fronts
or shocks. As mentioned in the introduction, cold fronts arise at the
boundaries of cool gas moving through hotter, less dense surrounding gas. Cold
fronts present as drops in the X-ray surface brightness (and therefore
density), with an accompanying jump in temperature. The changes in the
temperature and density tend to have similar amplitudes, resulting in
continuous pressure across the front. Shock fronts, on the other hand, follow
the Rankine-Hugoniot jump conditions, so across the front, we would expect to
see a drop in brightness, density, temperature, and pressure.

Figure 9: Gaussian smoothed (_top left_ ; see Fig. 1 left panel) and
adaptively smoothed (_top right_ ; see Fig. 1 right panel) X-ray images and
projected temperature map (_bottom left_) showing the regions defined for
extracting the radial profiles (labeled CF1-2 and SF1; see §5). The solid blue
(cold fronts) and red (shock) lines mark the location of the brightness edges
identified from the fits to the radial profiles. In order to aid
visualization, included is a cartoon schematic (_bottom right_) showing the
merger shock and a possible scenario of the CF regions connecting to form a
sloshing spiral structure. Overlaid on all panels are the 1.4 GHz (green) and
150 MHz (blue) radio contours of the two NAT sources.

To determine whether any of the X-ray structures seen in A119 are shocks or
cold fronts, we defined three wedge-shaped regions for further study (shown in
Fig. 9). Within each wedge region, we extract surface brightness profiles in
the 0.7-7 keV band from the merged, background subtracted, and exposure
corrected X-ray image with point sources excluded. Projected temperature
profiles were created by extracting the spectra from the same wedge regions
used to get surface brightness profiles, but using fewer annuli in an attempt
to minimize errors. From each annulus, the spectrum is extracted and fitted
using the method described in §4. Pressure and density profiles were created
following the deprojection method used to create the global profiles described
in §4.1.

We fit each surface brightness profile with a broken power-law model assuming
a 3D electron density profile defined as (Sarkar et al., 2022):

$n_{e}(r)\propto{\left\\{\begin{matrix}\left(\frac{r}{r_{\text{edge}}}\right)^{-\alpha_{1}}\quad\text{
if }r<r_{\text{edge}}\\\
\frac{1}{\text{jump}}\left(\frac{r}{r_{\text{edge}}}\right)^{-\alpha_{2}}\quad\text{
if }r\geq r_{\text{edge}}\end{matrix}\right.}$ (1)

where $n_{e}(r)$ is the 3D electron density at a given radius $r$,
$r_{\text{edge}}$ is the distance of the putative edge from the cluster
center, $jump$ is the density jump factor, and $\alpha_{1}$ and $\alpha_{2}$
are the power-law slopes inside and outside of the edge. We use this density
model to estimate the emission measure profile and project it onto the sky
plane, and fit the observed surface brightness profile by using the least
square fitting technique. We let the parameters, $\alpha_{1}$, $\alpha_{2}$,
$r_{\text{edge}}$, and the $jump$, free to vary during fitting.

Figure 10: Radial profiles of X-ray surface brightness, projected temperature,
and deprojected density and pressure for the regions discussed in §5 and shown
in Fig. 9. The surface brightness profiles include the best-fit model (red
line) and the 3D density profile (inset) obtained from the fit. The dashed
lines in all panels denote the best-fit edge radius and the shaded regions are
the 1$\sigma$ errors on the fitted radius.

Figure 10 shows the radial profiles of surface brightness, projected
temperature, deprojected pressure and density, the best-fit surface brightness
model, and the 3D density profiles (inset) for each region shown in Fig. 9.

##### Region CF1

The radial profiles of region CF1 show a drop in X-ray brightness and
deprojected density at 470″. Across the edge, the temperature increases from
$5.3^{+0.8}_{-0.7}$ to $8.5^{+1.7}_{-1.4}$ while the pressure remains
constant. The behavior of the radial profiles suggests that the NNW edge of
the “tear-drop” shape is a cold front. The best-fit power-law indices across
the edge are $\alpha_{1}=0.32\pm 0.04$ and $\alpha_{2}=0.94\pm 0.24$. The
best-fit density jump factor is $1.4\pm 0.2$. The edge radius, obtained from
the fit, is $r_{\text{edge}}=482.3\pm 18.2$″, and is shown by the dashed line
in Fig. 10 (upper left panel). Due to the asymmetric nature of the cluster
emission, the center of origin for the region labeled CF1 was chosen to match
the curvature of the feature of interest, rather than centered on the cluster
center. Thus, the profile for region CF1 has radial coordinates that are not
relative to the cluster center.

##### Region CF2

There is a drop in X-ray brightness seen near $\sim$ 225″. At this radius, the
temperature increases from $5.0^{+0.6}_{-0.5}$ to $6.4^{+0.8}_{-0.7}$, while
the pressure remains roughly constant across the edge, consistent with the
presence of a cold front. The best-fit power-law indices across the edge are
$\alpha_{1}=0.38\pm 0.13$ and $\alpha_{2}=0.79\pm 0.03$. The best-fit density
jump factor is $1.21\pm 0.10$. The edge radius obtained from the fit is
$r_{\text{edge}}=228.2\pm 3.4$″, and is indicated by the dashed line in Fig.
10 (upper right panel). This is consistent with the brightness edge north of
the cluster center that is noted in §3.1.2 and Fig. 6 (cyan arrows). Note that
while the bright spot seen just inside the edge (blue arc) of CF2 in Fig. 9
was not detected as a point source, we define a circular region containing
that clump and exclude it from the surface brightness and spectral
extractions. This was done in order to ensure proper detection of the cold
front in that region.

##### Region SF1

The radius of curvature for region SF1 was chosen to better fit the observed
edge rather than the cluster center. Therefore, while the surface brightness
profile of SF1 was fit using the radii from the off-centered coordinates, the
profiles for SF1 plotted in Fig. 10 (bottom panel) show the radius in terms of
distance from the cluster center. The radial profiles for region SF1 show a
drop in brightness, temperature, density, and pressure at a radius of 265″
from the cluster center. The temperature drops from $7.6^{+1.2}_{-0.9}$ to
$5.5^{+0.6}_{-0.5}$ while the pressure decreases by a factor of
$1.8^{+0.5}_{-0.4}$, consistent with a shock front. The best-fit power-law
indices are $\alpha_{1}=0.04\pm 0.6$ and $\alpha_{2}=1.14\pm 0.23$. The
density jump factor, obtained from the fit, is $1.31\pm 0.16$. The best-fit
edge radius is at a distance of $266\pm 10$″ from the cluster center, and is
indicated by the dashed line in Fig. 10 (bottom panel). This is consistent
with the brightness edge south of the cluster center seen in the GGM filtered
image (see §3.1.2 and yellow arrows in Fig. 6).

##  6 Discussion

###  6.1 Cold Fronts and a Shock

Figure 9 shows the cold fronts and shock front that were identified from the
radial profiles discussed above. We have two regions, CF1 and CF2, whose
radial profiles exhibit behaviors consistent with cold fronts. The exact
geometry of the connection between the two cold fronts in A119 is not
perfectly clear. Unfortunately, when examining regions just south of the
cluster core, we were unable to confidently identify a cold front in this
area. This was mostly due to the inability to disentangle a potential cold
front from both the X-ray deficient cavity near the eastern NAT and what could
possibly be regions of shock heated gas. However, we propose that the two cold
fronts may be connected through a sloshing spiral structure, beginning at
large radii in the north with CF1, and spiraling inwards in a clockwise
direction to the south, continuing towards the east, and finally NNW to region
CF2. This is shown schematically in Fig. 9 where we indicate with dashed lines
how the two cold fronts might be portions of a large-scale sloshing spiral
structure. The presence of a sloshing cold front could explain a few of the
observed peculiarities of A119. The possible sloshing spiral structure could
correspond to the elongated X-ray emission (i.e. the “teardrop”) seen in the
NE. The regions of displaced cool gas seen in the temperature map could have
been the result of sloshing induced motions of the ICM. A sloshing spiral
could also account for the sudden drop in abundance seen in Fig. 7 as the
sloshing motion tends to displace the metal-rich gas from the center, bringing
it in contact with the metal-poor gas of the surrounding ICM. Our proposed
scenario for how a sloshing spiral could have formed in A119 as a result of an
off-axis merger is explored further in §6.4.

Region SF1 shows compelling evidence of the presence of a shock front located
$\sim$250″ south of the cluster center. From the Rankine Hugoniot jump
conditions, we calculate the Mach number, $M$, of the shock front in terms of
the density jump (Landau & Lifshitz, 1959):

$M=\sqrt{\frac{2r}{\gamma+1-r(\gamma-1)}}$ (2)

where $r=\frac{\rho_{1}}{\rho_{2}}$, the subscripts 1 and 2 refer to the
density upstream (pre-shock) and downstream (post-shock), respectively, and
$\gamma=\frac{5}{3}$ is the adiabatic index assuming ideal gas. Using the
density jump obtained from the fit to the radial profile of region SF1,
$r=\frac{\rho_{1}}{\rho_{2}}=1.31\pm 0.16$, we calculate a Mach number of
$\mathcal{M}=1.21\pm 0.11$. Using the sound speed of the cluster (see §4.1),
we estimate the velocity of the shock front to be
$v_{\text{shock}}=c_{s}\mathcal{M}\simeq 1530\pm 140$ km/s. A temperature drop
from $7.6^{+1.2}_{-0.9}$ to $5.5^{+0.6}_{-0.5}$ with a pressure decreases of a
factor of $1.8^{+0.5}_{-0.4}$ is also measured across the shock (see §5).

A shock along the southern outskirts of the cluster could have formed as a
result of merger activity. This merger shock, when included in the overall
picture with our suggested sloshing spiral structure described above, suggests
A119 could have experienced recent or ongoing off-axis merger activity along
the N-S axis. This scenario is further explored through comparison with
simulations of cluster mergers in §6.4.

###  6.2 Velocity Estimates of Tailed Radio Galaxies

As discussed in the introduction, the tails of NAT radio sources are bent via
ram-pressure, which causes the radio jets to be bent backward. This can arise
from high orbital velocities of the host galaxy, bulk motion of the ICM, or
some combination of the two. A119 hosts two NAT sources which, notably, have
their tails oriented in roughly the same direction, along the NE-SW axis. The
similar tail orientations, along with both tails being aligned with the
overall NE-SW elongation of the X-ray emission, suggest a possible merger
origin for the two NATs. To determine whether this scenario is feasible, we
calculate the velocity of the host galaxies relative to the ICM that would be
required in order to produce the observed radio jet/lobe bending.

To estimate the velocities of the two NAT sources, we follow the methods of
Douglass et al. (2008, 2011) and Paterno-Mahler et al. (2013). If we assume
that the radio lobes are bent only due to the ram pressure from the
surrounding ICM, and assuming pressure balance between the lobes and the ICM,
we can apply Euler’s equation

$\frac{\rho_{r}v^{2}_{r}}{r_{c}}=\frac{\rho_{\text{\tiny
ICM}}v^{2}_{g}}{r_{r}}$ (3)

where $\rho_{r}$ is the mass density of the lobe, $r_{r}$ is the radius of the
lobe, $v_{r}$ is the velocity of the plasma in the lobe, $r_{c}$ is the radius
of curvature, $\rho_{\text{\tiny ICM}}$ is the mass density of the ICM, and
$v_{g}$ is the velocity of the galaxy relative to the ICM. For the mass
density of the ICM, we take $\rho_{\text{\tiny ICM}}=1.92n_{e}\mu m_{p}$,
where $\mu$ is the mean atomic mass of the cluster, taken to be 0.6, and
$n_{e}$ is the gas density found from the density profiles. If we assume in
situ particle acceleration in the lobes (O’Donoghue et al., 1993), then the
condition used to determine the internal density required to produce the
observed luminosity of the lobes is

$L_{rad}=\frac{1}{2}\epsilon\pi r^{2}_{r}\rho_{r}v^{3}_{r}$ (4)

where $\epsilon$ is the radiative efficiency and has typical values between
$0.001$ and $0.1$ (Birzan et al., 2004). The radiative efficiency measures the
conversion of kinetic energy into the observed radiation. It is defined as the
ratio of total radio luminosity, $L_{rad}$, to kinetic luminosity,
$L_{kin}=4PV/t$ (for relativistic jets; Churazov et al. (2001)). To calculate
the kinetic luminosity, we can estimate the ICM pressure from our pressure
profiles. The volume of the NAT region was calculated assuming an elliptical
geometry. We assume a typical AGN repetition rate of $t=5\times 10^{7}$ yr
(Birzan et al., 2004). The total radio luminosity is calculated by integrating
the flux between $\nu_{1}=10^{7}$ Hz and $\nu_{2}=10^{11}$ Hz

$L_{rad}=4\pi
D^{2}_{L}S_{\nu_{0}}\int^{\nu_{2}}_{\nu_{1}}\bigg{(}\frac{\nu}{\nu_{0}}\bigg{)}^{-\alpha}d\nu$
(5)

For both NAT sources, we use a spectral index of $\alpha=0.7$, calculated for
these sources in Feretti et al. (1999) from 1.4 GHz, 4.9 GHz, and 8.4 GHz VLA
observations. Reference fluxes and frequencies are taken from the VLA FIRST
catalogs.

Combining the luminosity condition of Eq. 4 with Euler’s equation from Eq. 3,
the velocity of the host galaxy relative to the ICM that would be required
bend the lobes is

$v_{g}=\bigg{(}\frac{2L_{rad}}{\epsilon\pi\rho_{\text{\tiny
ICM}}v_{r}r_{c}r_{r}}\bigg{)}^{1/2}$ (6)

For the velocity of the plasma in the radio lobes, $v_{r}$, we define upper
and lower limits of 0.08c to 0.2c (O’Neill et al., 2019; Terni De Gregory et
al., 2017). To calculate $\epsilon$ and $\rho_{\text{\tiny ICM}}$, we need to
determine estimates of the pressure and density of the ICM surrounding the NAT
sources. For the eastern NAT, we estimate the surrounding ICM pressure and
density from the radial profile shown in Fig. 11, between radii of $\sim
162-188\arcsec$. The pressure and density of the ICM surrounding the western
NAT is estimated from the radial profile shown in Fig. 12 between radii of
$\sim 321-365\arcsec$. Note that we assume a plane of sky orientation for the
two NAT sources and therefore, due to projection effects, the actual distances
from the NAT sources to the cluster center may be different than the assumed
distances.

For each of the NAT sources, we use Eq. 5 to calculate the total radio
luminosity. From the total radio and kinetic luminosities, we calculate the
efficiency ($\epsilon$). These are then used in Eq. 6 to determine the
velocity of the host galaxy relative to the ICM that would be required to
produce the observed bending of the NAT. The results for the individual NAT
sources are summarized below.

$\begin{array}[]{cc}\includegraphics[width=260.17464pt]{Eastern.pdf}&\includegraphics[width=173.44534pt]{csmooth_nat15_zoom_velbins.pdf}\end{array}$

Figure 11: Radial profile across the eastern NAT source used to calculate the
velocities in §6.2. Dashed lines denote the radial extent of the bin,
containing the NAT host galaxy, from which the pressure and density are
estimated.

$\begin{array}[]{cc}\includegraphics[width=260.17464pt]{Western.pdf}&\includegraphics[width=173.44534pt]{csmooth_nat16_zoom_velbins.pdf}\end{array}$

Figure 12: Radial profile across the western NAT source used to calculate the
velocities in §6.2. Dashed lines denote the radial extent of the bins,
containing the NAT host galaxy, from which the pressure and density are
estimated.

##### Eastern NAT (0053-015)

The total radio luminosity of 0053-015 was calculated to be $3.5\times
10^{40}$ erg s-1 using a reference flux of $47.5$ mJy at a frequency of 1435
MHz. To calculate the kinetic luminosity, we need to determine the pressure of
the overlying ICM around the NAT, as well as an estimate of the volume of the
NAT. The pressure of the surrounding ICM at the radius of the eastern NAT is
found to be $P=2.2\times 10^{-11}$ dyne cm-2. The volume of the NAT region was
calculated assuming an ellipsoidal geometry with $a=16.5$ kpc and $b=13.6$
kpc. We calculate a kinetic luminosity of $2.1\times 10^{43}$ erg s-1. From
the radio and kinetic luminosities, we find an efficiency of
$\epsilon=0.0016$. The mass density of the ICM was calculated to be
$\rho_{\text{\tiny ICM}}=2.1\times 10^{-27}$ g cm-3, where the gas density of
the surrounding ICM at the radius of the NAT is $n_{e}=1.1\times 10^{-3}$
cm-3. The radius of the lobe, measured from the 1.4 GHz maps provided in
Feretti et al. (1999), was found to be $r_{r}=6.6$ kpc, measured from the
eastern lobe which displays a greater degree of bending. The radius of
curvature was found to be $r_{c}=19.8$ kpc, again taken from the curvature of
the eastern lobe and measured from the 1.4 GHz map. It should be noted that
the bending in the eastern-most jet occurs near the edge of the visible
stellar halo of the galaxy. This could suggest that the galaxy has retained
some of its interstellar medium (ISM) which could be protecting the jet until
it reaches the ISM/ICM interface, at which point the jet is swept back. Thus
our measurement of the radius of curvature for this source should be
considered an upper limit.

From the values discussed above, we find that the relative velocity, in the
plane of the sky, between the host galaxy and the ICM (Eq. 6) is between 290
km s-1 and 450 km s-1, consistent with values previously seen in the sloshing-
induced bending of radio lobes (Paterno-Mahler et al., 2013). Additionally,
Feretti et al. (1999) gives the peculiar radial velocity, along the line of
sight, relative to the average cluster velocity as $|\Delta v|=1697$ km s-1
indicating a significant velocity component from the galaxy motion itself.

##### Western NAT (0053-016)

The total radio luminosity of 0053-016 was calculated as $2.3\times 10^{41}$
erg s-1 using a reference flux of 314.2 mJy at a frequency of 1435 MHz. The
kinetic luminosity was found to be $1.1\times 10^{44}$ erg s-1, where the
pressure of the surrounding ICM at the radius of the western NAT is
$P=2.3\times 10^{-11}$ dyne cm-2, and the volume was calculated assuming an
ellipsoidal geometry with $a=32.2$ kpc and $b=22.3$ kpc. From the total radio
and kinetic luminosities, we find an efficiency of $\epsilon=0.0020$.

From the density profile, we estimate the gas density of the surrounding ICM
at the radius of the NAT to be $n_{e}=1.1\times 10^{-3}$ cm-3, giving an ICM
mass density of $\rho_{\text{\tiny ICM}}=2.1\times 10^{-27}$ g cm-3. From the
1.4 GHz maps of Feretti et al. (1999), we measure $r_{r}=1.9$ kpc and
$r_{c}=12.3$ kpc.

Using the values above in Eq. 6, we calculate the velocity, in the plane of
the sky, of the host galaxy relative to the ICM to be between 1572 km s-1 and
2486 km s-1. These high relative velocities suggest that the bending of the
western NAT’s lobes is likely not solely induced from ICM sloshing motions.
Additionally, Feretti et al. (1999) gives the peculiar radial velocity, along
the line of sight, relative to the average cluster velocity as $|\Delta
v|=415$ km s-1, indicating the predominant galaxy motion is in the plane of
the sky. The western NAT appears to lie at a much further radius from the
cluster core and we would expect higher galaxy orbital velocities. Therefore,
there is likely some combination of both ICM bulk flow and galaxy orbital
motion that contribute to the shaping of the western NAT.

###  6.3 Optical Substructure in A119

Figure 13: Adaptively smoothed X-ray image showing locations of three optical
subclusters reported in previous studies. (Fabricant et al., 1993; Kriessler &
Beers, 1997; Ramella et al., 2007; Lee et al., 2016).

Figure 13 shows the locations of three subclusters identified in previous
optical studies (Fabricant et al., 1993; Kriessler & Beers, 1997; Ramella et
al., 2007; Lee et al., 2016). The substructures are oriented in a nearly
straight line along the N-S axis, providing further evidence of a merger
occurring along this direction. In addition to the three subclusters shown,
Tian et al. (2012) identify possible substructure located nearly directly
foreground or background to the core. They suggest this is indicative of
merger activity in the cluster core that is occurring along the line of sight.

Along with the alignment of optical substructures, the overall galaxy
distribution of A119 has also been shown to follow the NE-SW trend. Lee et al.
(2016) found that within a 0.64 Mpc radius of the cluster center the galaxy
distribution of cluster members is elongated to the NE. They further show that
the large-scale distribution of galaxies with radial velocities similar to
that of cluster galaxies form a filamentary structure that also extends to the
N-NE. This filament appears to extend nearly 8 Mpc north, connecting A119 to
another galaxy cluster, A116 ($z=0.066$) (see Fig. 11 of Lee et al. (2016)).

It is interesting to note that, if A119 is connected to A116 through a large-
scale filament, the orientation of the two NAT tails are also aligned along
this same axis. This has been seen previously in the alignment of jets from
wide-angle tail (WAT) radio sources with the orientation of the supercluser in
which the cluster belong. Novikov et al. (1999) showed 98% confidence in the
alignment of WAT jets and supercluster axes and they suggest this could be the
result of remnant drainage of matter along the large-scale filamentary
structures connecting the neighboring clusters. While the Novikov et al.
(1999) study focused on WAT rather than NAT sources, intriguingly, we observe
a similar occurrence in A119 of the two NAT tails having general alignment
with the large scale filamentary structure seen in Lee et al. (2016).

###  6.4 Comparison With Simulations

Cold fronts and shock fronts are produced by merger activity. Sloshing cold
fronts, in particular, are produced by an encounter of a massive cluster with
a much smaller subcluster that passes by with a non-negligible impact
parameter (Ascasibar & Markevitch, 2006). To further refine the merger
scenario for A119, we examined the idealized binary merger simulations from
the Galaxy Cluster Merger Catalog333https://gcmc.hub.yt (ZuHone et al. (2017);
hereafter GCMG) which are of this type. Our observation appears to be
qualitatively reproduced by a simulation of a 1:10 mass-ratio merger with an
initial impact parameter of $\sim$1 Mpc, a simulation originally presented in
ZuHone (2011). While the simulation results provide useful comparisons to our
observations, it should be noted that this simulation was not originally
designed to reproduce the features in A119.

In this simulation, the subcluster initially passes by the main cluster with a
closest approach of a $\sim$few hundred kpc. This first passage initiates the
sloshing motions. Because the subcluster is on a bound orbit with the main
cluster, it returns to the main cluster core $\sim 3.7$ Gyr later, narrowly
missing the core this time, and drives a shock front. This shock front
coincides roughly with the semicircular southern shock (the red line shown in
the cartoon schematic of Fig. 9) in A119.

The particular epoch of this simulation that best matched our observations is
$\sim 0.2$ Gyr after this second core passage, which has a shock front
surrounding the sloshing cold
fronts444http://gcmc.hub.yt/fiducial/1to10_b1/0255.html. This is shown in
Figure 14, which shows a gradient-filtered image of the X-ray surface
brightness and projected temperature from the simulation. It has been rotated
so that the positions of the cold and shock fronts approximately match those
in A119, and that the merger trajectory occurs in the NE-SW direction, as in
our preferred merger scenario. The relatively close match of this simulation
to the observed features of A119 indicates that a small subcluster approached
from the SW, passed by the main cluster from the south, traveled to the NNW
before turning around and passing the main cluster a second time from the N.
Currently, the whereabouts of the subcluster is unknown, but depending on the
exact line-of-sight angle, it may be co-aligned with the main cluster on the
sky or may be located to the SW. A simulation more tailored to match A119
would be necessary to resolve this question.

Figure 14: A cluster merger simulation which qualitatively matches our
preferred scenario for A119. _Left_ : Gradient-filtered X-ray surface
brightness (0.5-7.0 keV band). _Right_ : Projected spectroscopic-like
(Mazzotta et al., 2004) gas temperature. We include an inset of the cartoon
schematic from Fig.9 to aid in comparison of the features observed in A119
with the simulation results.

##  7 Summary & Conclusions

We have presented new _Chandra_ observations of the dynamically complex
cluster Abell 119. The X-ray emission from the ICM is asymmetric with an
elongation to the NE, resulting in a tear-drop like shape clearly seen in the
_Chandra_ image. Residual images show evidence of underlying structure in the
ICM that is explored in depth using radial profiles of the X-ray brightness,
temperature, density and pressure.

From the extracted radial profiles we identify two cold fronts which we
propose could be connected to form a sloshing spiral structure (see the
cartoon schematic in Fig. 9). Spectral maps show regions of cool gas along the
possible spiral, consistent with the cooler gas of the cluster core being
pulled out to larger radii via sloshing. If this is a sloshing cold front,
A119 could also provide useful insight into the role of off-axis mergers in
the formation of NCC clusters.

A surprising discovery was the presence of a shock front located $\sim 250$″
from the core and just outside of the potential sloshing cold front. From the
temperature profile we measure a temperature drop from $7.6^{+1.2}_{-0.9}$ to
$5.5^{+0.6}_{-0.5}$ across the shock. We calculate a Mach number for the shock
of $\mathcal{M}=1.21\pm 0.11$, corresponding to a velocity of
$v_{\text{shock}}=1530\pm 140$ km/s, when using the sound speed of the cluster
calculated in §4.1.

Another remarkable feature of this cluster is the orientation of the two NAT
sources with each other as well as the asymmetric X-ray emission. Both NAT
tails are oriented nearly parallel to each other despite the fact that the
jets are actually leaving their hosts in very different directions. In
addition to this, both tails are aligned with the NE-SW elongation of the
X-ray emission. We examine this further by using our results from the radial
profiles to calculate the velocities of the NAT host galaxies relative to the
ICM. Velocities of the eastern NAT are found to be consistent with values
expected from sloshing-induced bending. However, Feretti et al. (1999) reports
a significant line of sight velocity suggesting there is some contribution to
the bending of the tail from the galaxy motion itself. The western NAT was
found to have fairly high relative velocities. This, along with its distance
from the cluster center, suggest there is some combination of ICM bulk flow
with the host galaxy infall motions that is working to produce the observed
bending.

Following the NE-SW alignment trend of the X-ray emission and NAT tails, Lee
et al. (2016) showed that the optical distribution of cluster members is
likewise elongated to the north and the presence of optical substructures
oriented along the N-S direction. Beyond this, there is evidence of cluster
galaxies forming a filamentary structure which extends nearly 8 Mpc to the
N-NE which appears to connects A119 to A116. The alignment of the two NAT
sources with this filament could also be suggestive of cluster winds along the
supercluster axis (Novikov et al., 1999).

Comparison of our observations with the simulations of ZuHone et al. (2011)
helped to form a clearer picture of a possible merger scenario which we
propose could explain both the cold fronts and the shock in A119. We find that
the sloshing cold front produced as a result of the first passage of a much
smaller perturbing subcluster approaching from the SW. The simulations also
suggest that the shock front in the south could then be the result of a
secondary passage several Gyr later, with the subcluster now coming in from
the north. While these simulations were not originally designed to reproduce
the exact dynamics of A119, they provide a very useful visual comparison to
our observations.

A119 has proven to be a fairly complex system, playing host to two NAT
sources, a potential sloshing spiral, a merger shock, and possible connection
with a neighboring cluster through large-scale filamentary structures.
Untangling these features and investigating how they are connected to or
interacting with the cluster is important to further our understanding of
cluster formation and evolution. Our results show alignment of the elongated
X-ray emission, the optical substructures, and the flow directions of the
radio jets/lobes of both NATs. This, with the comparison to simulations, all
supports our picture of recent or on-going merger activity occurring in the
NE-SW direction.

We thank the anonymous referee for their helpful comments that greatly helped
improve this manuscript. CBW was supported by the National Aeronautics and
Space Administration through SAO Chandra Cycle 21 Award GO0-21123X, XMM-Newton
Cycle 17 Award 80NSSC19K0743, and STScI HST Cycle 27 Award GO-15994. CBW was
also supported by Massachusetts Space Grant Consortium Awards 580379, 633061,
699842, 742790, 793224, and 858617. Additional support for this work was
partially provided by the Chandra X-ray Center through NASA contract
NAS8-03060, the Smithsonian Institution, and by the Chandra X-ray Observatory
grant GO9-20112X.

## References

  * Arnaud (1996) Arnaud, K. A. 1996, Astron. Data Anal. Softw. Syst. V, 101, 17
  * Ascasibar & Markevitch (2006) Ascasibar, Y., & Markevitch, M. 2006, Astrophys. J., 650, 102, doi: 10.1086/506508
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, Annu. Rev. Astron. Astrophys., 47, 481, doi: 10.1146/annurev.astro.46.060407.145222
  * Becker et al. (1995) Becker, R. H., White, R. L., & Helfand, D. J. 1995, ApJ, 450, 559
  * Birzan et al. (2004) Birzan, L., Rafferty, D. A., McNamara, B. R., Wise, M. W., & Nulsen, P. E. J. 2004, Astrophys. J., 607, 800, doi: 10.1086/383519
  * Blanton et al. (2011) Blanton, E. L., Randall, S. W., Clarke, T. E., et al. 2011, ApJ, 737, doi: 10.1088/0004-637X/737/2/99
  * Blanton et al. (2009) Blanton, E. L., Randall, S. W., Douglass, E. M., et al. 2009, ApJ, 697, 95, doi: 10.1088/0004-637X/697/2/L95
  * Bliton et al. (1998) Bliton, M., Rizza, E., Burns, J. O., Owen, F. N., & Ledlow, M. J. 1998, MNRAS, 301, 609, doi: 10.1111/j.1365-8711.1998.01973.x
  * Burns (1998) Burns, J. O. 1998, Science (80-. )., 280, 400, doi: 10.1126/science.280.5362.40
  * Churazov et al. (2001) Churazov, E., Bruggen, M., Kaiser, C. R., Bohringer, H., & Forman, W. 2001, Astrophys. J., 554, 261, doi: 10.1086/321357
  * David et al. (1993) David, L. P., Slyz, A., Jones, C., Forman, W., & Vrtilek, S. D. 1993, ApJ, 412, 479, doi: 10.1086/172936
  * Dickey (1990) Dickey, J. M. 1990, ARA&A, 28, 215, doi: 10.1146/annurev.aa.28.090190.001243
  * Douglass et al. (2011) Douglass, E. M., Blanton, E. L., Clarke, T. E., Randall, S. W., & Wing, J. D. 2011, ApJ, 743, doi: 10.1088/0004-637X/743/2/199
  * Douglass et al. (2008) Douglass, E. M., Blanton, E. L., Clarke, T. E., Sarazin, C. L., & Wise, M. 2008, ApJ, 673, 763, doi: 10.1086/523886
  * Douglass et al. (2018) Douglass, E. M., Blanton, E. L., Randall, S. W., et al. 2018, ApJ, 868, 121, doi: 10.3847/1538-4357/aae9e7
  * Ebeling et al. (1996) Ebeling, H., Voges, W., Böhringer, H., et al. 1996, Mon. Not. R. Astron. Soc., 281, 799, doi: 10.1093/mnras/281.3.799
  * Edge et al. (1990) Edge, A. C., Stewart, G. C., Fabian, A. C., & Arnaud, K. A. 1990, MNRAS, 245, 559
  * Elkholy et al. (2015) Elkholy, T. Y., Bautz, M. W., & Canizares, C. R. 2015, Astrophys. J., 805, 1, doi: 10.1088/0004-637X/805/1/3
  * Fabricant et al. (1993) Fabricant, D., Kurtz, M., Geller, M., & Zabludoff, A. 1993, Astron. J., 105, 788, doi: 10.1017/CBO9781107415324.004
  * Feretti et al. (1999) Feretti, L., Dallacasa, D., Govoni, F., et al. 1999, Astron. Astrophys., 344, 472\.  https://arxiv.org/abs/9902019
  * Gastaldello et al. (2013) Gastaldello, F., Di Gesu, L., Ghizzardi, S., et al. 2013, Astrophys. J., 770, doi: 10.1088/0004-637X/770/1/56
  * Govoni et al. (2004) Govoni, F., Markevitch, M., Vikhlinin, A., et al. 2004, Astrophys. J., 605, 695, doi: 10.1086/382674
  * Intema et al. (2017) Intema, H. T., Jagannathan, P., Mooley, K. P., & Frail, D. A. 2017, Astron. Astrophys., 598, doi: 10.1051/0004-6361/201628536
  * Kempner & David (2004) Kempner, J. C., & David, L. P. 2004, Mon. Not. R. Astron. Soc., 349, 385, doi: 10.1111/j.1365-2966.2004.07534.x
  * Kriessler & Beers (1997) Kriessler, J. R., & Beers, T. C. 1997, Astron. J., 113, 80, doi: 10.1017/CBO9781107415324.004
  * Kriss et al. (1983) Kriss, G. A., Cioffi, D. F., & Canizares, C. R. 1983, ApJ, 272, 439
  * Laganá et al. (2019) Laganá, T. F., Durret, F., & Lopes, P. A. 2019, Mon. Not. R. Astron. Soc., 484, 2807, doi: 10.1093/mnras/stz148
  * Landau & Lifshitz (1959) Landau, L. D., & Lifshitz, E. M. 1959, Fluid Mechanics
  * Lee et al. (2016) Lee, Y., Rey, S.-C., Hilker, M., Sheen, Y.-K., & Yi, S. K. 2016, Astrophys. J., 822, 92, doi: 10.3847/0004-637x/822/2/92
  * Markevitch (2010) Markevitch, M. 2010, doi: 10.1142/9789814374552_0018
  * Markevitch et al. (1998) Markevitch, M., Forman, W. R., Sarazin, C. L., & Vikhlinin, A. 1998, Astrophys. J., 503, 77, doi: 10.1086/305976
  * Markevitch et al. (2002) Markevitch, M., Gonzalez, A. H., David, L., et al. 2002, Astrophys. J., 567, L27, doi: 10.1086/339619
  * Markevitch & Vikhlinin (2001) Markevitch, M., & Vikhlinin, A. 2001, Astrophys. J., 10, 95
  * Markevitch & Vikhlinin (2007) —. 2007, Phys. Rep., 443, 1, doi: 10.1016/j.physrep.2007.01.001
  * Markevitch et al. (2000) Markevitch, M., Ponman, T. J., Nulsen, P. E. J., et al. 2000, ApJ, 541, 542, doi: 10.1086/309470
  * Markevitch et al. (2003) Markevitch, M., Mazzotta, P., Vikhlinin, A., et al. 2003, Astrophys. J., 586, L19, doi: 10.1086/374656
  * Mazzotta et al. (2004) Mazzotta, P., Rasia, E., Moscardini, L., & Tormen, G. 2004, Mon. Not. R. Astron. Soc., 354, 10, doi: 10.1111/j.1365-2966.2004.08167.x
  * Mendygral et al. (2012) Mendygral, P. J., Jones, T. W., & Dolag, K. 2012, ApJ, 750, doi: 10.1088/0004-637X/750/2/166
  * Novikov et al. (1999) Novikov, D. I., Melott, A. L., Wilhite, B. C., et al. 1999, Mon. Not. R. Astron. Soc., 304, 5, doi: 10.1046/j.1365-8711.1999.02428.x
  * O’Donoghue et al. (1993) O’Donoghue, A. A., Eilek, J. A., & Owen, F. N. 1993, ApJ, 408, 428, doi: 10.1086/172600
  * O’Neill et al. (2019) O’Neill, B. J., Jones, T. W., Nolting, C., & Mendygral, P. J. 2019, Astrophys. J., 884, 12, doi: 10.3847/1538-4357/ab40b1
  * Owers et al. (2009) Owers, M. S., Nulsen, P. E., Couch, W. J., & Markevitch, M. 2009, Astrophys. J., 704, 1349, doi: 10.1088/0004-637X/704/2/1349
  * Paterno-Mahler et al. (2013) Paterno-Mahler, R., Blanton, E. L., Randall, S. W., & Clarke, T. E. 2013, ApJ, 773, doi: 10.1088/0004-637X/773/2/114
  * Poole et al. (2006) Poole, G. B., Fardal, M. A., Babul, A., et al. 2006, Mon. Not. R. Astron. Soc., 373, 881, doi: 10.1111/j.1365-2966.2006.10916.x
  * Postman & Lauer (1995) Postman, M., & Lauer, T. R. 1995, ApJ, 440, 28
  * Ramella et al. (2007) Ramella, M., Biviano, A., Pisani, A., et al. 2007, A&A, 51, 39
  * Randall et al. (2008) Randall, S., Nulsen, P., Forman, W. R., et al. 2008, Astrophys. J., 688, 208, doi: 10.1086/592324
  * Randall et al. (2009) Randall, S. W., Jones, C., Kraft, R., Forman, W. R., & O’Sullivan, E. 2009, Astrophys. J., 696, 1431, doi: 10.1088/0004-637X/696/2/1431
  * Roediger et al. (2011) Roediger, E., Brüggen, M., Simionescu, A., et al. 2011, Mon. Not. R. Astron. Soc., 413, 2057, doi: 10.1111/j.1365-2966.2011.18279.x
  * Roediger et al. (2013) Roediger, E., Kraft, R. P., Forman, W. R., Nulsen, P. E., & Churazov, E. 2013, Astrophys. J., 764, doi: 10.1088/0004-637X/764/1/60
  * Roettiger et al. (1993) Roettiger, K., Burns, J., & Loken, C. 1993, ApJ, 407, L53, doi: 10.1086/186804
  * Rossetti et al. (2013) Rossetti, M., Eckert, D., De Grandi, S., et al. 2013, Astron. Astrophys., 556, 1, doi: 10.1051/0004-6361/201321319
  * Saglia et al. (1997) Saglia, R. P., Burstein, D., Baggley, G., et al. 1997, Mon. Not. R. Astron. Soc., 292, 499, doi: 10.1093/mnras/292.3.499
  * Sanders et al. (2016) Sanders, J. S., Fabian, A. C., Russell, H. R., Walker, S. A., & Blundell, K. M. 2016, Mon. Not. R. Astron. Soc., 460, 1898, doi: 10.1093/mnras/stw1119
  * Sarkar et al. (2022) Sarkar, A., Randall, S., Su, Y., et al. 2022, Astrophys. J. Lett., 935, L23, doi: 10.3847/2041-8213/ac86d4
  * Sarkar et al. (2023) —. 2023, Astrophys. J., 944, 132, doi: 10.3847/1538-4357/acae9f
  * Schindler & Muller (1993) Schindler, S., & Muller, E. 1993, AAP, 272, 137.  https://ui.adsabs.harvard.edu/abs/1993A{&}A...272..137S
  * Terni De Gregory et al. (2017) Terni De Gregory, B., Feretti, L., Giovannini, G., et al. 2017, Astron. Astrophys., 608, 1, doi: 10.1051/0004-6361/201730878
  * Tian et al. (2012) Tian, J. T., Yuan, Q. R., Zhou, X., et al. 2012, Res. Astron. Astrophys., 12, 1381, doi: 10.1088/1674-4527/12/10/002
  * Tittley & Henriksen (2005) Tittley, E. R., & Henriksen, M. 2005, Astrophys. J., 618, 227, doi: 10.1086/425952
  * Valdarnini & Sarazin (2021) Valdarnini, R., & Sarazin, C. L. 2021, Mon. Not. R. Astron. Soc., 504, 5409, doi: 10.1093/mnras/stab1126
  * Walker et al. (2017) Walker, S. A., Hlavacek-Larrondo, J., Gendron-Marsolais, M., et al. 2017, Mon. Not. R. Astron. Soc., 468, 2506, doi: 10.1093/mnras/stx640
  * West et al. (1995) West, M. J., Jones, C., & Forman, W. 1995, ApJL, 451, L5, doi: 10.1086/309673
  * Wong et al. (2008) Wong, K., Sarazin, C. L., Blanton, E. L., & Reiprich, T. H. 2008, Astrophys. J., 682, 155, doi: 10.1086/588272
  * ZuHone (2011) ZuHone, J. A. 2011, Astrophys. J., 728, doi: 10.1088/0004-637X/728/1/54
  * ZuHone et al. (2017) ZuHone, J. A., Kowalik, K., Öhman, E., Lau, E., & Nagai, D. 2017, Astrophys. J. Suppl. Ser., 234, 4, doi: 10.3847/1538-4365/aa99db
  * ZuHone et al. (2010) ZuHone, J. A., Markevitch, M., & Johnson, R. E. 2010, Astrophys. J., 717, 908, doi: 10.1088/0004-637X/717/2/908
  * ZuHone et al. (2011) ZuHone, J. A., Markevitch, M., & Lee, D. 2011, Astrophys. J., 743, doi: 10.1088/0004-637X/743/1/16
  * ZuHone et al. (2016) ZuHone, J. A., Miller, E. D., Simionescu, A., & Bautz, M. W. 2016, Astrophys. J., 821, 6, doi: 10.3847/0004-637x/821/1/6


# Fast simulation of quantum algorithms using circuit optimization

Gian Giacomo Guerreschi<EMAIL_ADDRESS>Intel Labs

###### Abstract

Classical simulators play a major role in the development and benchmark of
quantum algorithms and practically any software framework for quantum
computation provides the option of running the algorithms on simulators.
However, the development of quantum simulators was substantially separated
from the rest of the software frameworks which, instead, focus on usability
and compilation. Here, we demonstrate the advantage of co-developing and
integrating simulators and compilers by proposing a specialized compiler pass
to reduce the simulation time for arbitrary circuits. While the concept is
broadly applicable, we present a concrete implementation based on the Intel
Quantum Simulator, a high-performance distributed simulator. As part of this
work, we extend its implementation with additional functionalities related to
the representation of quantum states. The communication overhead is reduced by
changing the order in which state amplitudes are stored in the distributed
memory. We then implement a compiler pass to exploit the novel functionalities
by introducing special instructions governing data movement as part of the
quantum circuit. Those instructions target unique capabilities of simulators
and have no analogue in actual quantum devices. To quantify the advantage, we
compare the time required to simulate random circuits with and without our
optimization. The simulation time is typically halved.

Introduction —  The field of quantum computation has recently graduated from
scientific research to technology development. One of the most visible changes
has been the creation of increasingly sophisticated software frameworks to
enable the execution of algorithms on actual quantum devices Green2013 ;
Wecker2014 ; Steiger2018 ; Qiskit ; Cirq ; smith2016practical ; Killoran2019 ;
Orquestra ; Braket . A major component of the frameworks is the compiler,
sometimes also called transpiler, mapper or scheduler depending on its
specific role Javadiabhari2014 ; Cowtan2019 ; Almudever2020 . The compiler is
required to translate a relatively abstract algorithm into instructions
executable by the electronic controller of the quantum device. Every advanced
framework also includes the option of running the algorithm with a simulator,
i.e. classical software specialized in simulating quantum circuits.

Simulators have at least two important applications. On one hand they allow to
benchmark quantum algorithms without the limitations of existing hardware, in
particular finite coherence time Smelyanskiy2016 ; Pednault2017 ;
DeRaedt2019_48qub ; Haner2017 ; Khammassi2017 ; Jones2019 ; HuaweiHiQ ;
Villalonga2020 . On the other hand, when equipped with realistic noise models,
simulators help characterizing hardware noise and decoherence and, possibly,
inform on its origin OBrien2017 ; Villalonga2019 . In spite of their
importance, no specialized optimization is integrated in the popular compiler
frameworks that mostly deal with restrictions on the qubit connectivity and on
the type of available gates. Noticeably, ProjectQ compiler Steiger2018
accepts multi-qubit operations in the algorithm’s description without forcing
their decomposition in 1- and 2-qubit gates, a feature suitable for simulation
since multi-qubit gates are not native in most quantum hardware.

In this work, we present a new compiler pass dedicated to reducing the
simulation time for arbitrary quantum circuits. It takes advantage of the
specific implementation of the simulator and is, therefore, a backend-aware
compilation pass in contrast to hardware-agnostic passes earlier in the
compiler chain Haener2018 ; Cowtan2019 . We adopt the Intel Quantum Simulator
(IQS) Smelyanskiy2016 ; Guerreschi2020 , a high-performance simulator
exhibiting both multi-threading and distributed parallelization, and extend it
by introducing a flexible way to represent the quantum state as a distributed
vector. The functionality is available to all users of IQS via its open-source
distribution at https://github.com/iqusoft/intel-qs. Then we create a compiler
pass minimizing the communication overhead incurred when multiple computing
nodes are involved in the simulation. The optimization is based on the
distinction between gates that require inter-node communication and those that
do not, minimizing the number of the former ones by means of a greedy search.

Previous works on distributed simulators have proposed similar methods to
reduce the communication overhead DeRaedt2007 ; Haner2017 ; HuaweiHiQ , while
ad-hoc circuit manipulation have been used to simulate, with a hybrid
Schrödinger-Feynman simulator developed by NASA Villalonga2019 , the kind of
random circuits proposed by Google for quantum supremacy Arute2019 .

To evaluate the efficacy of the optimization, we compare the time to simulate
random circuits with and without the proposed compiler pass. We vary the ratio
of 1- and 2-qubit gates in the circuits and report a reduction by one order of
magnitude of the gates requiring inter-node communication. This translates in
the overall simulation time being typically halved.

Extending Intel Quantum Simulator —  IQS is a state-of-the-art Schrödinger
simulator which represents the full state of $n$ qubits by storing all its
$2^{n}$ complex amplitudes. All operations required by quantum circuits,
namely state preparation and several 1- and 2-qubit gates, are implemented as
linear manipulation of the amplitudes. Since the state vector grows
exponentially in size with the number of qubits, it soon exceeds the memory
available in a single computing node. Hence the necessity of distributing the
state across multiple nodes, each with its own local memory.

The implementation is presented in references Smelyanskiy2016 ; Guerreschi2020
, together with the communication pattern of the Message Passing Interface
(MPI) protocol. Here, we provide a high-level summary of the properties
exploited by the compiler pass. State preparation and measurement are done
with little or no communication and efficiently parallelized over computing
nodes and processor’s cores (via OpenMP). Gate operations are limited to those
gates that update amplitudes in pairs: arbitrary 1-qubit gates and 2-qubit
gates in the form of controlled 1-qubit gates. Notice that the 2-qubit gate
known as CNOT is part of this set and, thus, IQS can simulate all circuits
(arbitrary 1-qubit gates and CNOT allow universal quantum computation).

We divide the qubits into two distinct groups of local, respectively global,
qubits Haner2017 . The number of local vs global qubits depends on the number
$k$ of MPI processes, specifically there are $\lfloor\log_{2}k\rfloor$ global
qubits and $m=n-\lfloor\log_{2}k\rfloor$ local qubits. Inter-process
communication is required by all gates acting on any of the global qubits
(apart from the special case of a controlled 2-qubit gate with global control
and local target), and the communication overhead is about an order of
magnitude. Apart from the communication overhead, the cost of 1- and 2-qubit
gates is not significantly different. Figure 1 shows the time required to
simulate 1- and 2-qubit gates as a function of the involved qubits. In the
numerical experiments, the first $m=28$ qubits are local and is evident that
the communication overhead starts from gates on qubits with index 28 or
greater. The simulations were run on Frontera supercomputer at the Texas
Advanced Computing Center in which each computing node is equipped with 2
sockets of Intel® Xeon® Processor 8280 CPU (24 cores per socket).

Figure 1: Simulations of $n=35$ qubits using 128 MPI processes, each running
on one node of Frontera supercomputer. Each process uses 24 OpenMP threads.
(Top) Time to execute a 1-qubit gate as a function of the qubit involved. When
the gate is executed on qubit with index $m=n-\lfloor\log_{2}k\rfloor$ (with
$k$ being the number of MPI processes), the communication between the MPI
tasks is happening between sockets of the same node. For qubits with index
larger than $m$, communication is between distinct nodes. (Bottom) Time to
execute a 2-qubit gate as a function of the involved qubits. CNOT is the
controlled Pauli X gate.

In the original IQS implementation it is not possible to choose which qubits
are local and which are global, by default the first $m$ qubits are the local
ones. Consider an arbitrary $n$-qubit state, it can be written in the
computational basis as:

$\displaystyle\ket{\psi}=\sum_{i\in\\{0,1\\}^{n}}\alpha_{i}\ket{i_{0}}\ket{i_{1}}\dots\ket{i_{n-1}}$

where $i=[i_{n-1}\dots i_{1}i_{0}]$ can be seen as the binary representation
of integer $i=\sum_{q=0}^{n-1}i_{q}2^{q}$. IQS stores the quantum state as a
distributed vector whose $i$-th component corresponds to amplitude
$\alpha_{i}$. To provide opportunities for fast simulation, we want to
increase the flexibility of this representation and, specifically, we allow to
change the order in which the amplitudes are stored. Consider an arbitrary
permutation $\sigma$ of the qubit indices: under such permutation, qubit $q$
will be associated with bit $\sigma(q)$ of the index of the vector’s
component. Therefore, amplitude $\alpha_{i}$ is stored as component
$j(i,\sigma)=\sum_{q=0}^{n-1}i_{q}2^{\sigma(q)}$ of the state vector.

We introduce the new function PermuteQubits$(\sigma^{\prime})$ to update the
representation of quantum states in IQS. Consider that the current qubit
permutation is $\sigma$ and that we want to update it to $\sigma^{\prime}$.
The qubits are reordered in three steps: 1) local qubits (according to
$\sigma$) are reordered among themselves, 2) global qubits (according to
$\sigma$) are reordered among themselves, and 3) pairs formed by one local and
one global qubits (according to either $\sigma$ or, with opposite roles,
$\sigma^{\prime}$) are exchanged. The cost of the three steps is analogous to
that of a single local gate, a single global gate, and one global gate per
exchanged pair respectively. We discuss the implementation of PermuteQubits in
Appendix A, together with the pseudo-code to divide every transformation from
qubit permutation $\sigma$ to $\sigma^{\prime}$ in the three steps. As part of
this contribution, we added the distributed implementation of SWAP gates to
IQS.

The property of local vs global qubits is still valid, but the distinction now
depends on permutation $\sigma$. The local qubits are those with index $q$
such that $\sigma(q)<m$. There still are $m=n-\lfloor\log_{2}k\rfloor$ local
qubits and $\lfloor\log_{2}k\rfloor$ global qubits.

Compiler pass to optimize circuit simulation —  From the benchmarks reported
in Figure 1, it is clear that gates acting on global qubits takes much longer
to simulate than gates on local qubits alone. We propose a compiler pass that
goes through an arbitrary quantum circuit and eliminates the communication
overhead from most gates by adding qubit-reordering instructions.
Communication is required to execute such instructions in IQS but, when they
are chosen carefully, the overhead is reduced. As we quantify in the next
section, even for unstructured circuits the simulation time is halved.

[t] Compiler pass specialized for IQS simulation

1:

2:

3:

4:

5:

6:$G\leftarrow$ DAG of circuit

7:$\sigma\leftarrow$identity permutation

8:$C\leftarrow\\{\\}$

9:while $G$.vertices$\neq\emptyset$ do

10: for $v$ in $G$.vertices do

11: if ($v$ has no incoming edges) $\land$ ($v$ acts on local qubits according
to $\sigma$) then

12: add ApplyGate($v$) to $C$

13: remove $v$ from $G$

14: end if

15: end for

16: $\tilde{\sigma}\leftarrow$identity permutation

17: $\tilde{a}\leftarrow 0$

18: for $l$ in local qubits according to $\sigma$ do

19: for $g$ in global qubits according to $\sigma$ do

20: $\sigma^{\prime}\leftarrow\sigma\circ(l\,g)$

21: $a^{\prime}\leftarrow$ count $G$.vertices on local qubits according to
$\sigma^{\prime}$

22: if $a^{\prime}>\tilde{a}$ then

23: $\tilde{\sigma}\leftarrow\sigma^{\prime}$

24: $\tilde{a}\leftarrow a^{\prime}$

25: end if

26: end for

27: end for

28: if $\tilde{a}>0$ then

29: $\sigma\leftarrow\tilde{\sigma}$

30: add PermuteQubits($\sigma$) to $C$

31: else

32: $v\leftarrow$ vertex of $G$ without incoming edges

33: add $v$ to $C$

34: remove $v$ from $G$

35: end if

36:end while

The compiler pass is described as Algorithm 1 below. It assumes that the
quantum circuits is provided in terms of a directed acyclic graph (DAG) in
which each vertex corresponds to a quantum gate and each edge to a logical
dependency: the target vertex (i.e. gate) of the edge must be performed after
the source vertex. The DAG may or may not be created considering that quantum
gates may commute, but in our case we include commutativity to eliminate
unnecessary edges. While common Guerreschi2018 ; Lao2019a ; Itoko2020 , this
is not the only way to describe circuits as DAGs Sivarajah2020a . In
describing the compiler pass, we will use the terms gate and vertex
interchangeably.

The algorithm starts by scheduling all gates that act on local qubits only and
do not have unresolved logical dependencies (lines 5:10). The latter condition
corresponds to vertices without incoming edges. When it is not possible to
schedule any further gate, the algorithm considers all possible qubit
permutations obtained by exchanging one local and one global qubit according
to the current qubit permutation $\sigma$ (lines 11:22). There are $m(n-m)\leq
n^{2}/4$ permutations to consider. Each permutation is evaluated depending on
the number of gates that it allows to simulate on local qubits, and the
permutation with the highest count is recorded as $\tilde{\sigma}$ in line 18.
The compiler add instruction PermuteQubit$(\tilde{\sigma})$ if it removes
communication from at least one gate, otherwise schedule the first gate
without incoming edges even if it requires communication (lines 23:30). The
algorithm iterates until all gates are scheduled.

We observe that the compiler pass uses a simplified interpretation of the
simulation times reported in Figure 1. In particular, it considers only
whether a qubit is local or global and neglects the fine-grained pattern of
the actual simulation time. In addition, it does not consider the actual value
of the communication overhead, but is only aware that gates on global qubits
takes much longer than gates on local qubits. For example, the plot in Figure
1(top) would be abstracted as a step function: low simulation time until qubit
$m-1$ and high simulation time from qubit $m$ onward. By updating the qubit
permutation $\sigma$, the compiler is able to remove certain communication
overhead and, therefore, minimize the overall simulation time. More advanced
compiler passes may use the full information from simulator’s benchmarks.

Numerical results —  To quantify the advantage provided by the compiler pass,
we estimate the time to simulate random circuits before and after their
optimization. We consider random circuits of the following form: fix the
number of gates, every gate has a probability $p$ to be a 2-qubit gate
otherwise with probability $1-p$ it is a 1-qubit gate. Due to stochasticity,
the number of two-qubit gates is not strictly determined. The specific type of
gate is selected uniformly at random in the sets $\\{H,Y\\}$ and $\\{CNOT\\}$
for 1- and 2-qubit gates respectively. Finally the qubits involved in each
gate are chosen uniformly at random. Recall that the simulation time is
largely independent on the specific type of gate and the choice is therefore
arbitrary and inconsequential.

The random circuits we consider are simply a sequence of gates chosen at
random, and do not aim at reproducing a Haar-random unitary on $n$ qubits.
They have been chosen to reflect the common situation in which programmers of
quantum algorithms are not familiar with the design of the simulator and,
thus, do not tailor their code (via gate selection or qubit indexing) for a
specific simulator. In addition, while random circuits are not associated with
a specific application, they are representative of algorithms like the Quantum
Approximate Optimization Algorithm to solve combinatorial problems on class of
random instances Farhi2014  or the dynamical evolution of spin systems with
random interactions Salathe2015 .

Figure 2: Effect of the compiler pass to optimize quantum circuits for
simulation with IQS. (Top) Simulation time of random circuits of 1050 gates on
$n=35$ qubits. The timings have been estimated for the same configuration used
in Figure 1: 128 computing nodes of Frontera supercomputer at TACC. There are
$m=28$ local qubits. Datapoints represent the average over 20 random circuits
and the vertical bars indicate one standard deviation. (Bottom) Reduction of
the simulation time between original and optimized circuit. The blue line is
based on direct benchmarks of IQS, the red line is obtained by assuming that
operations requiring communications have an overhead of $8\times$ with respect
to operations implementable locally. This overhead is in line with the
benchmarks in Figure 1.

In the first study, we fix the number of qubits to $n=35$ and vary the
probability $p$. Despite only a fraction $7/35=1/5$ of the qubits being
global, the simulation time is reduced by between 32.3% (for $p=1$) and 56.0%
(for $p=0$), according to the formula $1-t_{\text{opt}}/t_{\text{orig}}$ in
which $t_{\text{orig}}$ is the simulation time for the original circuit and
$t_{\text{opt}}$ that for the optimized circuit. The results are shown in
Figure 2 and have been obtained by averaging over 20 random circuits of
$(30\,n)$ gates on $n=35$ qubits. The diminished advantage when $p$ increases
is not related to the different probability that a 2-qubit gate requires
communication compared to a 1-qubit gate, in fact for CNOT gates communication
is required if and only if the target is global irrespective of the control
qubit. We attribute the diminished reduction to the fact that a 2-qubit gate
(requiring communication) excludes a larger part of the remaining circuit from
execution since it has, on average, more outgoing edges in the DAG description
of the circuit (see condition in line 6 of Algorithm 1).

The simulation time has been estimated by summing up the time required to
simulate every gate separately. We start by benchmarking the execution of each
type of gate acting on every qubit or pair of qubits. We also benchmark the
time to reorganize the classical data representing the quantum state following
a call to PemuteQubits. This can be achieved by benchmarking the time to
execute a single SWAP gate on every pair of qubits since each new permutation
is related to the previous one by the composition with a single 2-cycle or
transposition (see line 15 in Algorithm 1). Notice that, unlike from the
circuit optimization that is based on the dicotomy local vs global qubits, to
estimate the simulation time of both the original and optimized circuits we
consider the full information content provided by benchmarks like those in
Figure. 1.

Figure 3: Effect of the compiler pass to optimize quantum circuits for
simulation with IQS. Fraction of gates requiring communication for the
original (blue) and optimized (red) circuit. (Top) Fixing $n=50$, we vary the
percentage of global qubits from 10% to 90%. (Bottom) Fixing the ratio global
vs total qubits to 1:5, we vary $n$. In both panels, datapoints are averaged
over 20 random circuits and the vertical bars indicate one standard deviation.
The random circuits have $(30\,n)$ gates.

In the second study, we fix the probability of 2-qubit gates to $p=0.3$ and
vary the number of local and global qubits. This choice corresponds to about
two 1-qubit gates every 2-qubit operation. In Figure 3 we report the number of
gates causing communication overhead (here including also the qubit reordering
operations) as a fraction of the total number of original operations. In the
top panel, the horizontal axis represents the percentage of qubits that are
global when the total number of qubits is $n=50$. This simulation size is at
the threshold of simulability, even considering using all computing nodes of
current supercomputers. However, we show in the bottom panel that, apart from
small-size effects for $n\leq 30$, the behavior is independent on $n$. This
allows us to estimate the advantages of the optimization at extremely large
scale (even beyond the simulability threshold expected around 50 qubits).
Specifically, when $1/5$ of the qubits are global, the fraction of gates
requiring communications is reduced from $\sim 20\%$ to $\sim 3.5\%$
independently of the number of qubits. Assuming an overhead of about $8\times$
for the communication (see bottom panel of Figure 2), this would correspond to
a reduction of simulation time by
$1-t_{\text{opt}}/t_{\text{orig}}=1-(0.965+0.035\times 8)/(0.8+0.2\times
8)\simeq 48\%$. Notice that our simple optimization may be improved, but with
limited gain since even the impossible case of completely eliminating
communication would lead to a reduction of $1-t_{\text{no-
comm}}/t_{\text{orig}}=1-1/(0.8+0.2\times 8)\simeq 58\%$ in this scenario.

Conclusion —  Quantum circuits are compiled for execution on actual quantum
hardware, but an important feature of quantum computing frameworks is the
possibility of simulating algorithms on classical computers. We propose a new
compiler pass specialized for circuit simulation. Our results demonstrate that
the co-development and integration between circuit compiler and simulator is
largely beneficial and reduces the simulation time in half. The integration
requires dedicated functions in both software programs: on one hand, we extend
Intel Quantum Simulator, a high-performance, distributed and open-source
simulator, with flexible data structures and qubit re-ordering
functionalities; on the other hand, we propose a novel compiler pass based on
the minimization of the number of operation requiring communication between
computing nodes. We believe that compiler passes specialized for circuit
simulation will become an essential component of every quantum computation
framework.

Acknowledgements —  The author acknowledges the Texas Advanced Computing
Center (TACC) at The University of Texas at Austin for providing HPC resources
that have contributed to the research results reported in this work.
http://www.tacc.utexas.edu The author thanks Fabio Baruffa, Justin Hogaboam,
and Nicolas P.D. Sawaya for helpful discussions.

## References

  * [1] Alexander S. Green, Peter LeFanu Lumsdaine, Neil J. Ross, Peter Selinger, and Benoit Valiron.  Quipper: A scalable quantum programming language.  Proceedings PLDI 13, 48:333–342, 2013.
  * [2] David Wecker and Krysta M. Svore.  LIQUi$|\rangle$: A software design architecture and domain-specific language for quantum computing.  arXiv:1402.4467, 2014.
  * [3] Damian S. Steiger, Thomas Häner, and Matthias Troyer.  ProjectQ: An Open Source Software Framework for Quantum Computing.  Quantum, 2:49, 2018.
  * [4] The Qiskit Contributors.  Qiskit: An open-source framework for quantum computing, 2019.
  * [5] The Cirq Contributors.  Cirq, a python framework for creating, editing, and invoking noisy intermediate scale quantum (NISQ) circuits.  https://github.com/quantumlib/Cirq.
  * [6] Robert S. Smith, Michael J. Curtis, and William J. Zeng.  A practical quantum instruction set architecture.  arXiv:1608.03355, 2016.
  * [7] Nathan Killoran, Josh Izaac, Nicolás Quesada, Ville Bergholm, Matthew Amy, and Christian Weedbrook.  Strawberry Fields: A Software Platform for Photonic Quantum Computing.  Quantum, 3:129, 2019.
  * [8] The Orquestra Contributors.  Orquestra.  https://www.zapatacomputing.com/orquestra/.
  * [9] The Braket Contributors.  Braket.  https://aws.amazon.com/braket/.
  * [10] Ali Javadiabhari, Shruti Patil, Daniel Kudrow, Jeff Heckey, Alexey Lvov, Frederic T. Chong, and Margaret Martonosi.  ScaffCC: A framework for compilation and analysis of quantum computing programs.  CF ’14 Proceedings of the 11th ACM Conference on Computing Frontiers, 2014.
  * [11] Alexander Cowtan, Silas Dilkes, Ross Duncan, Alexandre Krajenbrink, Will Simmons, and Seyon Sivarajah.  On the qubit routing problem.  Leibniz International Proceedings in Informatics, LIPIcs, 135(5):5:1—-5:32, 2019.
  * [12] Carmen G. Almudever, Lingling Lao, Robert Wille, and Gian G. Guerreschi.  Realizing Quantum Algorithms on Real Quantum Computing Devices.  Proceedings of the 2020 Design, Automation and Test in Europe Conference and Exhibition, DATE 2020, pages 864–872, 2020.
  * [13] Mikhail Smelyanskiy, Nicolas P. D. Sawaya, and Alán Aspuru-Guzik.  qHiPSTER: The quantum high performance software testing environment.  arXiv:1601.07195, 2016.
  * [14] Edwin Pednault, John A. Gunnels, Giacomo Nannicini, Lior Horesh, Thomas Magerlein, Edgar Solomonik, Erik W. Draeger, Eric T. Holland, and Robert Wisnieff.  Breaking the 49-qubit barrier in the simulation of quantum circuits.  arXiv:1710.05867, 2017.
  * [15] Hans De Raedt, Fengping Jin, Dennis Willsch, Madita Willsch, Naoki Yoshioka, Nobuyasu Ito, Shengjun Yuan, and Kristel Michielsen.  Massively parallel quantum computer simulator, eleven years later.  Computer Physics Communications, 237:47–61, apr 2019.
  * [16] Thomas Häner and Damian S. Steiger.  0.5 petabyte simulation of a 45-qubit quantum circuit.  In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis on - SC 17. ACM Press, 2017.
  * [17] N. Khammassi, I. Ashraf, X. Fu, C.G. Almudever, and K. Bertels.  QX: A high-performance quantum computer simulation platform.  In Design, Automation & Test in Europe Conference & Exhibition (DATE), 2017. IEEE, mar 2017.
  * [18] Tyson Jones, Anna Brown, Ian Bush, and Simon C. Benjamin.  QuEST and high performance simulation of quantum computers.  Scientific Reports, 9(1), jul 2019.
  * [19] The Huawei HiQ Team.  Huawei hiq: A high-performance quantum computing simulator and programming framework.  http://hiq.huaweicloud.com.
  * [20] Benjamin Villalonga, Dmitry Lyakh, Sergio Boixo, Hartmut Neven, Travis S. Humble, Rupak Biswas, Eleanor G. Rieffel, Alan Ho, and Salvatore Mandrà.  Establishing the quantum supremacy frontier with a 281 Pflop/s simulation.  Quantum Science and Technology, 5(3):034003, 2020.
  * [21] Thomas E. O’Brien, B. Tarasinski, and Leo DiCarlo.  Density-matrix simulation of small surface codes under current and projected experimental noise.  npj Quantum Information, 3:39, 2017.
  * [22] Benjamin Villalonga, Sergio Boixo, Bron Nelson, Christopher Henze, Eleanor Rieffel, Rupak Biswas, and Salvatore Mandrà.  A flexible high-performance simulator for verifying and benchmarking quantum circuits implemented on real hardware.  npj Quantum Information, 5(1):1–16, 2019.
  * [23] Thomas Häner, Damian S. Steiger, Krysta M. Svore, and Matthias Troyer.  A Software Methodology for Compiling Quantum Programs.  Quantum Science and Technology, 3:020501, 2018.
  * [24] Gian Giacomo Guerreschi, Justin Hogaboam, Fabio Baruffa, and Nicolas P. D. Sawaya.  Intel Quantum Simulator: A cloud-ready high-performance simulator of quantum circuits.  Quantum Science and Technology, 5:034007, 2020.
  * [25] K. De Raedt, K. Michielsen, H. De Raedt, B. Trieu, G. Arnold, M. Richter, Th. Lippert, H. Watanabe, and N. Ito.  Massively parallel quantum computer simulator.  Computer Physics Communications, 176(2):121–136, jan 2007.
  * [26] Frank et al. Arute.  Quantum supremacy using a programmable superconducting processor.  Nature, 574:505–510, 2019.
  * [27] Gian Giacomo Guerreschi and Jongsoo Park.  Two-step approach to scheduling quantum circuits.  Quantum Science and Technology, 3(4):045003, 2018.
  * [28] Lingling Lao, Daniel M. Manzano, Hans van Someren, Imran Ashraf, and Carmina G. Almudever.  Mapping of quantum circuits onto NISQ superconducting processors.  arXiv:1908.04226, 2019.
  * [29] Toshinari Itoko, Rudy Raymond, Takashi Imamichi, and Atsushi Matsuo.  Optimization of quantum circuit mapping using gate transformation and commutation.  Integration, 70:43–50, 2020.
  * [30] Seyon Sivarajah, Silas Dilkes, Alexander Cowtan, Will Simmons, Alec Edgington, and Ross Duncan.  tiket : A Retargetable Compiler for NISQ Devices.  arXiv:2003.10611, 2020.
  * [31] Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.  A quantum approximate optimization algorithm.  arXiv:1411.4028, 2014.
  * [32] Y. Salathé, M. Mondal, M. Oppliger, J. Heinsoo, P. Kurpiers, A. Potočnik, A. Mezzacapo, U. Las Heras, Lucas Lamata, Enrique Solano, S. Filipp, and Andreas Wallraff.  Digital quantum simulation of spin models with circuit quantum electrodynamics.  Physical Review X, 5:021027, 2015.

$\bullet$

##  Appendix A Appendix:
Transformation between two qubit permutations

Consider the current qubit permutation $\sigma$ such that amplitude
$\alpha_{i}$ of the computational basis state
$\ket{i_{n-1}}\dots\ket{i_{1}}\ket{i_{0}}$ corresponds to the component
$j=\sum_{q=0}^{n-1}i_{q}2^{\sigma(q)}$ of the state vector. We want to update
the qubit permutation to $\sigma^{\prime}$ taking into account that $m$ of the
$n$ qubits are local and the remaining $n-m$ are global. The IQS method
implementing the update proceeds in three steps: The first one reorders local
qubits only, the second one reorders global qubits only, and the third step
exchanges pairs formed by one local and one global qubits.

The three steps have separate implementations: reordering local qubits is
obtained by copying and rearranging data in each local memory independently,
reordering global qubits is obtained by swapping the local memory of two
distinct processes, and exchanging a local and a global qubit corresponds to
the execution of a SWAP gate from the point of view of data movement. A SWAP-
like gates is implemented via the same communicate-compute-communicate
approach of controlled 1-qubit gates described in [13], but with a different
pattern of inter-process communication.

Below, we provide pseudocode to divide the transformation
$\sigma\mapsto\sigma^{\prime}$ in three transformations that, composed
sequentially, correspond to the desired one. For clarity of terminology, we
will say that qubit $q$ is mapped to position $\sigma(q)$ and, viceversa, that
position $p$ is mapped to qubit $\sigma^{-1}(p)$. Inverting the permutation is
efficient and takes $\mathcal{O}(n)$ operations. Algorithm A constructs the
intermediate permutations such that
$\sigma\mapsto\sigma_{1}\mapsto\sigma_{2}\mapsto\sigma^{\prime}$ satisfies the
criteria expressed above.

[h] Decompose transformation $\sigma\mapsto\sigma^{\prime}$

1:

2:$\sigma\leftarrow$ current qubit permutation

3:$\sigma^{\prime}\leftarrow$ desired qubit permutation

4:$\sigma_{1}\leftarrow$ identity permutation

5:$\sigma_{2}\leftarrow$ identity permutation

6:for $p=0,1,\dots,m-1$ do

7: $p^{\prime}\leftarrow\sigma^{\prime}(\sigma^{-1}(p))$

8: while $p^{\prime}\geq m$ do

9: $p^{\prime}\leftarrow\sigma^{\prime}(\sigma^{-1}(p^{\prime}))$

10: end while

11: $\sigma_{1}^{-1}(p^{\prime})\leftarrow\sigma^{-1}(p)$

12:end for

13:for $p=m,m+1,\dots,n-1$ do

14: $p^{\prime}\leftarrow\sigma^{\prime}(\sigma_{1}^{-1}(p))$

15: while $p^{\prime}<m$ do

16: $p^{\prime}\leftarrow\sigma^{\prime}(\sigma_{1}^{-1}(p^{\prime}))$

17: end while

18: $\sigma_{2}^{-1}(p^{\prime})\leftarrow\sigma_{1}^{-1}(p)$

19:end for