# Burton’s Curse: The Impact of Bulk Flows on the Galactic Longitude–Velocity
Diagram and the Illusion of a Continuous Perseus Arm

J. E. G. Peek Space Telescope Science Institute, 3700 San Martin Drive,
Baltimore, MD 21218, USA  Department of Physics & Astronomy, Johns Hopkins
University, Baltimore, MD 21218, USA  Kirill Tchernyshyov Department of
Astronomy, University of Washington, Seattle, WA, USA  Marc-Antoine Miville-
Deschenes AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot,
Sorbonne Paris Cité, F-91191 Gif-sur-Yvette, France

###### Abstract

In this work we demonstrate that the Perseus Arm is not a continuous structure
of molecular gas in the second quadrant. We first show that the observed,
distanced-resolved velocity structure of the Galaxy in the outer disk is
capable of creating illusory spiral arms, as was first proposed by Burton
(1971). Second, we measure the distances to a collection of CO clouds at
velocities consistent with the Perseus arm with $135^{\circ}<l<160^{\circ}$.
We find these distances using 3D dust maps from Green et al. (2019). We
determine that these molecular cloud do not preferentially lie at the distance
of a purported Perseus arm, but rather extend over 3 kpc in distance, with
some evidence for a closer, high pitch angle structure between 1 and 1.5 kpc
away. Finally, we demonstrate that velocity perturbations of the amplitude
found near the Perseus arm can wreak havoc on our interpretation of the
longitude–velocity diagram for more than half of the Milky Way disk.

Milky Way dynamics — Milky Way disk — Milky Way rotation — Spiral arms —
Interstellar medium — Interstellar dust extinction

##  1 Introduction

Much of the star formation in disk galaxies happens in their spiral arms, and
thus understanding spiral structure is key to understanding the process of
star formation on large scales. On the other hand, we know star formation to
be mediated by very small-scale stellar feedback and turbulent processes,
which can be studied in most detail within the Milky Way. Thus, understanding
the true nature of Milky Way spiral structure is critical to connecting the
environment of star formation to its eventual result.

While our need to accurately map Milky Way spiral structure has long been
recognized, our location within the Milky Way disk limits our ability to
discern the distances to the molecular clouds and star forming regions that
make up spiral arms. Emission tracers of the ionized, atomic, and molecular
gas that make up the arms contain no intrinsic distance information, and
stellar structures embedded within the arms are often obscured by intervening
material. Even _Gaia_ ’s state of the art, large scale precise parallactic
distance information to individual stars is largely limited to nearby
($\lesssim$ 2 kpc), relatively unobscured, and unconfused regions Gaia
Collaboration et al. (2018). Very Long Baseline Array (VLBI) parallax
measurements toward high-mass star formation regions (HMSFRs) provides a much
more limited set of stellar distances, but one that is unaffected by dust, and
traces the star forming regions of interest (Reid et al., 2014, Reid14).

Historically, astronomers have estimated distances to gas seen in emission by
taking advantage of the analytic relationship between distance and radial
velocity that is created by the differential rotation of the Milky Way disk
(van de Hulst et al., 1954). In particular, in the outer Milky Way there is a
one-to-one, monotonic relationship between distance and velocity that, under
the assumption of a known and unperturbed rotation curve, allows us to exactly
map out the structure of the disk beyond the solar circle (Levine et al.
(2006), Kalberla & Kerp (2009), Miville-Deschênes et al. (2017a)). This method
represents one of the only ways to connect individual regions of star
formation (as traced by e.g. HMSFRs) into larger spiral structures across the
Galactic disk. This linking is very important, as the pitch angle and
continuity of spiral structure is a key discriminant between different
theories of spiral arm formation and persistence (Reid14).

Burton (1971, B71) examined this velocity-based map-making method and pointed
out a fundamental flaw. By modeling a _completely smooth_ ISM with the
velocity field implied by the Lin et al. (1969) spiral structure density wave
theory, B71 was able to reproduce many of the observed features in the
Galactic Hi longitude–velocity diagram. The implication of this work is that
realistic spiral structure _velocity fields_ generate perturbations to the
longitude–velocity diagram that mimic in shape and amplitude those generated
by spiral arm density fluctuations themselves. Burton & Bania (1974) followed
up on this analysis, comparing distances to star clusters to the implied
distances from velocities in Hi . While these original theories of spiral
structure are by no means the only ones under consideration for the Milky Way,
recent work has shown that the Milky Way does indeed host line-of-sight
velocity fluctuations of more than 10 km/s stretching for many kpc, which may
be capable of generating these kinds of distortions (Tchernyshyov & Peek,
2017; Tchernyshyov et al., 2018, TPZ18).

The Perseus arm in the outer Milky Way has long been used as a test bed for
theories of spiral structure, as it is (a) a purported nearby spiral arm, (b)
in the outer Galaxy, (c) easily visible in the north, and (d) seemingly
continuous across a large swath of the Galactic disk. Early work argued that
the velocity and density field in this region was consistent with the two-arm
spiral shock model of spiral structure (Roberts, 1972). Subsequent
observations of the positions and motions of gas and stars in the Perseus arm
have also been found to agree with the Roberts (1972) model (Burton & Bania,
1974; Humphreys, 1976; Sparke & Dodd, 1978; Wouterloot & Habing, 1985; Sakai
et al., 2019), though often with caveats. Humphreys (1976) find that the line-
of-sight velocities of supergiant stars have the right signs, but twice the
expected amplitude. Sparke & Dodd (1978) note the Perseus arm gap discussed in
this paper, but conclude that it could in fact be the result of a non-linear
effect hypothesized to occur in the spiral shock model. Sakai et al. (2019,
Sakai19) find that the 3D velocities of HMSFRs in the Perseus arm are
consistent with the Roberts (1972) model but that HMSFRs in the rest of the
purported spiral arms are not consistent with spiral shock-like models. The
development of other theories of spiral structure continued, in particular
variations on swing amplification (Toomre, 1981; Sellwood & Carlberg, 1984).
Unlike Roberts (1972), these models did not make the sort of explicit
predictions for the Perseus arm and so were typically not compared with
observations.

TPZ18 showed that the velocity fields around Perseus star-forming regions are
inconsistent with these kinds of long-lived persistent spiral models, and that
material arm models, in which spiral arms are created and destroyed on less
than a winding time, are more consistent with the observed velocity fields.
These material arm models are incapable of making long, shallow pitch angle
spiral arms, which stands in contradiction to the observed length and tight
winding in the Perseus arm as identified in the longitude–velocity diagram.
Here, shallow pitch angle refers to the small angle made by the arm and a
constant radius circle (see e.g. Sellwood & Masters, 2021). New observations
of hot stars, Cepheids, and open clusters mined from Gaia EDR3 (Gaia
Collaboration et al., 2021) also seem to show a more loosely wound Perseus arm
(Poggio et al., 2021). Longer, more tightly wound spiral arms, as argued for
in Reid14, rely on sparse HMSFRs linked together by features in the molecular
(CO) and atomic (Hi) longitude–velocity diagrams. While the distances and
velocities of HMSFRs are not in doubt, can we really trust the links between
them, or are they an illusion, as described in B71?

In this work, we put the Perseus arm to the test, examining whether a key part
of this arm is indeed contiguous, or rather an illusion generated by Galactic
velocity fields as described in B71. The part of Perseus we examine here does
not have any known masers, but is important in connecting up star forming
regions in the second and third quadrants, and thus establishing the shallow
pitch angle of the purported arm. In Section 2 we describe the 3D dust maps,
CO, Hi, and HMSFR maser data we use. In Section 3 we show how previously
measured radial velocity fields provide for the possibility of spiral arm
illusions in the Perseus arm, and describe our method for measuring the
distance to individual molecular clouds in this region. In Section 4 we
discuss the implications of our measured distance distribution of molecular
clouds in the region under consideration and we conclude in Section 5.

##  2 Data

###  2.1 3D dust maps

We use the 3D dust data products described in Green et al. (2019, G19). In
brief, G19 uses the optical and near IR photometry of 799 million stars, along
with their _Gaia_ data release 2 Gaia Collaboration et al. (2018) parallax
measurements, to infer the combinations of stellar type, distance, and
foreground reddening that are possible for each star. By combining these
probability distribution functions for all stars near each line of sight, G19
produced a 3D dust map, with typical angular resolution 3.4′, and distance
bins of 0.125 distance moduli. We will refer to the total reddening profile
along the line of sight from this map as ${\rm
E\left(B-V\right)}\left(d\right)$ and to the radial differential reddening
profile as $\Delta{\rm E\left(B-V\right)}\left(d\right)$. Adjacent bins are
somewhat correlated, as a smooth spatial prior is applied in an iterative
scheme to the map. Five realizations of the map are provided and we access
them using the _dustmaps_ code (M. Green, 2018). Additionally, each pixel in
the map is labeled with its maximum and minimum reliable distance; lines of
sight with few stars or opaque clouds cannot reliably return dust information
beyond a certain point.

We note that G19 finds that the Reid14 HMSFRs are associated with regions of
dense dust, but that there is no clear evidence for the kind of shallow pitch
angle spiral structure found in Reid14 in the G19 dust maps themselves. We
believe this result is consistent with what we find in this work.

###  2.2 Kinetic Tomography

In order to show that the true velocity field of the Milky Way is capable of
producing the significant velocity crowding effects as described in B17 we use
the kinetic tomography map published in TPZ18. TPZ18 used the 1.53 $\mu$m
diffuse interstellar band absorption feature in stars from APOGEE (Majewski et
al., 2017; Zasowski et al., 2015) to measure the velocity field across the
Milky Way disk in the first three quadrants within about 4 kpc. These maps
confirmed the overall results of the previous maps (Tchernyshyov & Peek, 2017,
TP17), which showed that there are bulk flows of material in the Galaxy that
can stretch many kpc at tens of km s-1 beyond what is expected from Galactic
rotation. Additionally, we use the results of TP17 as an illustration of our
results in Section 4.

###  2.3 ISM gas tracer observations

We use the compilation of 21-cm Hydrogen (Hi ) and microwave carbon monoxide
(CO) observations described in TP17 to provide a template for contaminating
reddening, as described in §3.2. Practically speaking, this means we used CO
observations from Dame et al. (2001, DHT01) and Hi observations from Kalberla
et al. (2005).

###  2.4 Cloud identification

To select individual molecular clouds, we used the decomposition of the DHT01
map provided by (Miville-Deschênes et al., 2017a, MML17). MML17 used a
hierarchical cluster identification method to collect together Gaussian
velocity components of the DHT01 CO survey into a list of 8246 molecular
clouds across the Galactic disk. The data contain both a list of the 8246
clouds and their bulk properties, but also a list of the individual Gaussian
components associated with each cloud. This augmented and slightly modified
MML17 data set is described in Appendix A and can be found
athttps://doi.org/10.7910/DVN/QR9CFW (catalog 10.7910/DVN/QR9CFW).

###  2.5 High Mass Star Formation Regions

Finally, we also use the HMSFRs compiled in Reid14 and Sakai19. This
represents the largest single collection of HMSFRs with parallactic distances
ever assembled. These targets are very helpful for studying the global
structure of star formation in the Milky Way because they are found in areas
of active star formation and contain masers, whose parallaxes (and hence
distances) and radial velocities can be measured accurately using VLBI. This
combination is very powerful for tracing out the spatial and velocity
structure of the Milky Way’s star forming disk.

##  3 Methods & Analysis

In this section we develop methods and conduct analyses in pursuit of answers
to two questions. First, does the previously measured velocity field of the
Milky Way in the outer Galaxy allow for illusory spiral structure in the outer
Milky Way? Second, if so, is the dense material in the velocity feature
actually distributed over a much larger range of distances?

###  3.1 Illusory Spiral Structure is Possible

In B71 it was shown that, in the inner Galaxy, significant spiral-like
structure could be generated using the velocity field expected from the
density-wave theory derived in Lin et al. (1969) with _no density variation at
all_. While the idea of a perfectly smooth ISM is inconsistent with
essentially every ISM observation, the analysis proved the point that “the
profiles are more sensitive to small variations in the streaming motions than
to small variations in the hydrogen density”. This revelation was key to the
movement toward studying CO as a tracer of spiral arm structure rather than Hi
(R. Benjamin, priv. comm).

To verify this result, and thus demonstrate that velocity crowding in the
Milky Way disk is possible, we apply the disk velocity field results of TPZ18
to a perfectly smooth disk. We lay down a smooth grid of points across the
second quadrant, out to a distance of 4 kpc where we believe the TPZ18 map to
be largely complete and accurate. Then we perform two experiments. First, we
use a standard flat, 220 km s-1 rotation curve to place each of these points
at a velocity, shown in panel 1 of Figure 1. Then, we use the radial velocity
field as measured by TPZ18 to place these same smoothly distributed points on
a second longitude–velocity diagram, shown in panel 2 of Figure 1.

It is quite clear from these synthetic longitude–velocity diagrams that the
Perseus arm in the second quadrant would be easily visible in the
longitude–velocity diagram given the TPZ18 map of the velocity field and a
perfectly smooth Galactic disk. While this “velocity crowding” situation
provides the opportunity for illusory features in the longitude–velocity
diagram, it does not mean that what we see in the maser-free sections of the
Perseus arm is illusory. It is certainly possible that dense material also
physically exists at these distances, and that the velocity crowding is only
piling on lower density, non-star-forming material at these velocities.
Indeed, the HMSFRs from Reid14 and various maps of young stars (e.g. Xu et
al., 2018) show that at least in some regions there do exist star forming
regions at the distance of the Perseus arm.

Figure 1: The Galactic longitude–velocity diagram in the Galactic midplane
($b=0$). For each panel Hi column is shown in grayscale and CO column is shown
in “magma” color map. The bounds of the regions we focus on in §3.2 are marked
in dotted gray lines. The velocity track of the purported Perseus arm is shown
in dashed green from Choi et al. (2014). Masers from Reid14 and Sakai19 are
shown as green dots. Overplotted in brown is a representation of a perfectly
smooth ISM. The brown points are smoothly distributed over the Galactic disk
within 4 kpc, projected using a flat rotation curve (top) and the velocity
field derived in Tchernyshyov et al. (2018) (bottom). This figure shows that
track of the Perseus arm can be generated by the velocity field derived in
Tchernyshyov et al. (2018) without density fluctuations.

###  3.2 Illusory Spiral Structure Impacts the Perseus Arm Identification

In §3.1 we have shown that velocity fields observed in the outer Milky Way
disk are sufficiently strong to create spiral arms in the longitude–velocity
diagram in the absence of any density enhancements in the disk at all. To
determine whether dense gas features are being erroneously connected into
longer spiral structures in practice, we seek to measure the true distances to
molecular clouds thought to reside in kinematic spiral arms. We focus our
attention on the region between $l$=135 and $1$=160 at velocities consistent
with the Perseus arm.

We focus on this region because it contains no HMSFRs with parallactic
distances and no evidence of young stars clustered at the expected distance
Reid14, but has a very clear spiral feature in both Hi and CO gas and is near
enough ($\sim$2 kpc) that the G19 reddening map is generally trustworthy. If
this feature exists only in velocity space, and not in distance space, we
believe there is currently no evidence of a continuous Perseus spiral arm
across the second and into the third quadrant of the outer Galaxy.

First we select a sample of clouds from MML17, which are part of the visually
identifiable Perseus arm. First we constrain ourselves in longitude:
$135^{\circ}<l<160^{\circ}$. Then we make overall constraints on velocity:
$-65{\rm~{}km~{}s}^{-1}<V_{\rm LSR}<-20{\rm~{}km~{}s}^{-1}$. Finally we fine-
tune the velocity boundaries to follow the arm in longitude–velocity space:
$V_{\rm LSR}/{\rm~{}km~{}s}^{-1}<(l-169^{\circ})/^{\circ}$ and
$0.6\left(V_{\rm
LSR}+35{\rm~{}km~{}s}^{-1}\right)/{\rm~{}km~{}s}^{-1}>(l-160^{\circ})/^{\circ}$
These constraints result in 114 CO clouds, shown in Figure 2. For each one of
these clouds we measure a distance using the G19 map. If the ISM consisted
only of these clouds, the procedure for determining their distance would be
relatively simple – measure the average $\Delta{\rm
E\left(B-V\right)}\left(d\right)$ along the line of sight and find the peak.
Works like Schlafly et al. (2014) and Zucker et al. (2020) have made very
detailed studies of distances to clouds and sub-clouds using some of the
techniques in G19, but these have typically been at closer distances or
further off the plane. Unfortunately, there are many other structures along
the line of sight, including the CO clouds in very strong Local Arm between 0
and -20 km s-1 (see figure 1) and the more diffuse phases of the ISM along the
line of sight not detected in DMT01.

To mitigate the effect of foreground ISM structures along the line of sight we
also choose a control region, which we call an “off”, consisting of the area
of sky more than 7 pixels away from any region of the cloud and less than 16
pixels in the DMT01 CO map (see left panels of Figure 3). Each pixel in the
DMT01 map is 7.5′, such that the off region consists of the area roughly 1∘ to
2∘ from the cloud in question. For each of the cloud and off region we find
all the corresponding dust pixels in G19 that have complete data out to at
least 4 kpc.

For both the cloud and the off region we make a weighting map consisting of
total column of material along the line of sight excepting the cloud itself.
To do this we use the fiducial column density conversion factors for 21-cm Hi
and microwave CO 1-0 transitions Bolatto et al. (2013)

:

$N\left(\rm H~{}\sc{I}\right)=1.8\times 10^{18}\frac{{\rm T_{B}}d\nu}{\rm
K~{}km~{}s^{-1}}$ (1)

and

$X\left(\rm CO\right)=2.0\times 10^{20}\frac{{\rm T_{B}}d\nu}{\rm
K~{}km~{}s^{-1}}$ (2)

apply them to the CO and Hi maps, sum them, and subtract the contribution from
the cloud according to MML17. These maps allow us to weight our average
measurement of the dust for the clouds and the offs by the inverse of the
unrelated ISM column. This weighting allows us to focus on sightlines that are
less contaminated by intervening and unrelated material, though in practice we
find using an unweighted average does not change the conclusions of this work.

We look to estimate the distance to each cloud using the dust reddening
measurements of G19. To compute this distance we construct the difference
between our on-cloud weighted $\Delta{\rm E\left(B-V\right)}\left(d\right)$
and the off-cloud weighted $\Delta{\rm E\left(B-V\right)}\left(d\right)$. In
the very near field, $d<0.7$ kpc, we find that it is impossible to fully
subtract off the effects of the Local Arm material, which is very strong and
quite variable. Thus, we focus our attention on the range 0.7 kpc $<d<$ 5 kpc
in our analysis. For each cloud we find the interval over this distance range
that captures 68.3% of the nonnegative cloud-off $\Delta{\rm
E\left(B-V\right)}\left(d\right)$ that covers the least distance, and use the
center of this range as our distance estimator (see right panels of Figure 3).
We find this method is much more resilient to the complex, non-Gaussian
distributions that are the result of the complex intervening ISM than e.g.
first moment measurements.

To measure errors we turn to block bootstrapping. Standard bootstrap analysis
(Efron, 1982) is extremely useful for situations where errors on samples are
uncorrelated but otherwise quite complex. Unfortunately, while the errors in
distance for this measurement are indeed complex (being totally dominated by
intervening clouds rather than measurement noise) they are by no means
independent: two adjacent pixels are very likely to be contaminated by the
same intervening cloud. Thus, standard bootstrapping is both technically
inapplicable and practically inaccurate. For this reason we turn to block
bootstrapping, a procedure of resampling with replacement, but for subregions
(blocks), rather than individual pixels Hall (1985).

To make subregions for a given region (cloud or off) we first find the average
position of the region in Galactic coordinates. Then we compute a clock-angle
on the surface of the sky for each of the pixels in the region around this
center position. Finally, we divided the region up into blocks based on the
pixels clock-angle rank, with the largest block having no more than one extra
pixel than the smallest block. This creates largely contiguous blocks of
pixels of equal size. We use 12 blocks for each region. As we perform this
block bootstrap we also bootstrap on the five realizations of the dust map. We
find that the map realization bootstrapping adds very little to the overall
errors, consistent with our expectation that intervening clouds, and not
measurement noise, are our largest source of uncertainty.

While we have built a system to find cloud distances that is robust as
possible in the presence of intervening material, some clouds are too small or
too obscured to find with this method. First, we reject 26 additional clouds
for having a $\pm 1\sigma$ error range that is more than a factor of 1.5 in
distance. Second, we reject 8 clouds for having very heavy-tailed dust
profiles, defined by having a 68.3% ($\pm 1\sigma$) interval 3.5 times larger
than their 38.3% ($\pm 0.5\sigma$) interval or having a 86.6% ($\pm
1.5\sigma$) interval 3.5 times larger than their 68.3% interval. This leaves
us a sample of 81 of 114 initial clouds which we report in Table 1 and map in
Figure 4. We note that this procedure makes it clear that we are biased away
from finding the most distant clouds, as they will tend to be fainter and
smaller overall. We expect that this sample is thus more complete closer to
the Sun. We see in Figure 4 that the bulk of CO clouds at the Perseus arm
velocity in this longitude range are not consistent with being at the distance
of the Perseus arm. The longitude–velocity diagram of the Perseus arm is shown
with distances assigned in Figure 5, clearly demonstrating a very wide range
of distances for clouds very near each other in velocity space.

Table 1: Clouds with measured dust distances.  Cloud | MML17 Cloud | $l$ | $b$ | VLSR | Area | Av. Flux | Dust Distance
---|---|---|---|---|---|---|---
Number | Number | [∘] | [∘] | [km s-1] | [deg2] | [K km s-1] | [kpc]
1 | 2168 | 137.65 | 1.41 | -38.8 | 1.31 | 2.8 | $1.92^{+0.09}_{-0.06}$
2 | 2178 | 136.60 | 1.17 | -36.6 | 0.78 | 4.4 | $1.86^{+0.09}_{-0.03}$
3 | 2182 | 148.13 | 0.25 | -33.8 | 1.92 | 2.6 | $1.39^{+0.03}_{-0.08}$
4 | 2190 | 144.39 | -1.20 | -31.4 | 2.27 | 3.7 | $1.84^{+0.02}_{-0.16}$
5 | 2201 | 151.76 | -1.19 | -29.3 | 2.62 | 4.2 | $2.12^{+0.09}_{-0.34}$
6 | 2702 | 136.75 | 1.10 | -41.8 | 0.23 | 3.8 | $2.20^{+0.31}_{-0.06}$
7 | 2712 | 137.00 | 1.37 | -39.1 | 1.03 | 3.4 | $1.76^{+0.09}_{-0.01}$
8 | 2731 | 145.50 | 0.86 | -35.6 | 2.06 | 4.2 | $1.26^{+0.03}_{-0.09}$
9 | 2736 | 150.52 | -0.94 | -32.7 | 3.20 | 2.4 | $1.46^{+0.08}_{-0.06}$
10 | 2749 | 151.30 | -0.53 | -25.7 | 1.16 | 5.5 | $2.31^{+0.12}_{-0.59}$
11 | 2757 | 150.58 | -0.88 | -23.2 | 0.41 | 3.6 | $1.85^{+0.09}_{-0.35}$
12 | 3397 | 140.85 | -1.17 | -40.4 | 0.66 | 4.4 | $1.15^{+0.26}_{-0.09}$
13 | 3398 | 141.56 | -1.22 | -41.0 | 1.86 | 1.7 | $1.30^{+0.08}_{-0.05}$
14 | 3417 | 140.68 | 0.81 | -38.2 | 1.30 | 2.0 | $1.94^{+0.09}_{-0.06}$
15 | 3422 | 156.78 | -1.85 | -35.6 | 1.44 | 2.5 | $1.51^{+0.03}_{-0.00}$
16 | 3429 | 144.80 | 1.01 | -35.9 | 0.41 | 4.6 | $1.23^{+0.08}_{-0.04}$
17 | 3430 | 146.03 | 1.01 | -36.4 | 0.17 | 4.6 | $1.20^{+0.16}_{-0.07}$
18 | 3436 | 156.95 | -2.94 | -32.3 | 1.52 | 2.9 | $1.48^{+0.07}_{-0.04}$
19 | 3440 | 157.55 | 0.85 | -33.0 | 1.91 | 2.4 | $1.45^{+0.03}_{-0.06}$
20 | 3444 | 145.24 | -1.08 | -30.5 | 1.69 | 2.5 | $1.32^{+0.24}_{-0.14}$
21 | 3454 | 146.24 | 1.25 | -29.6 | 1.09 | 3.0 | $1.28^{+0.03}_{-0.09}$
22 | 3477 | 154.37 | -1.15 | -20.7 | 0.91 | 2.6 | $1.69^{+0.09}_{-0.45}$
23 | 4324 | 135.64 | -0.13 | -40.9 | 1.66 | 1.5 | $2.28^{+0.05}_{-0.07}$
24 | 4338 | 142.80 | -0.96 | -41.7 | 1.80 | 1.4 | $1.06^{+0.05}_{-0.02}$
25 | 4356 | 144.87 | 1.56 | -41.9 | 0.55 | 1.9 | $1.39^{+0.03}_{-0.06}$
26 | 4360 | 138.90 | -0.41 | -37.8 | 2.92 | 1.2 | $2.66^{+0.04}_{-0.06}$
27 | 4363 | 139.87 | 0.29 | -40.3 | 0.67 | 1.5 | $1.46^{+0.22}_{-0.31}$
28 | 4386 | 150.43 | 0.42 | -35.2 | 1.81 | 2.1 | $1.53^{+0.03}_{-0.04}$
29 | 4388 | 151.23 | 0.73 | -35.9 | 1.08 | 1.4 | $3.39^{+0.06}_{-0.34}$
30 | 4394 | 152.13 | 1.57 | -34.8 | 2.05 | 1.3 | $3.43^{+0.07}_{-0.08}$
31 | 4396 | 155.59 | 2.54 | -36.0 | 0.97 | 1.5 | $3.58^{+0.03}_{-0.19}$
32 | 4402 | 147.78 | -0.53 | -33.9 | 1.25 | 1.9 | $1.32^{+0.10}_{-0.07}$
33 | 4403 | 156.45 | -0.37 | -33.0 | 3.52 | 1.4 | $1.38^{+0.05}_{-0.07}$
34 | 4404 | 151.16 | -0.32 | -33.6 | 0.11 | 2.4 | $2.44^{+0.16}_{-0.19}$
35 | 4407 | 149.11 | -0.16 | -33.8 | 0.41 | 1.8 | $1.53^{+0.03}_{-0.15}$
36 | 4409 | 149.31 | 0.75 | -36.4 | 1.98 | 1.6 | $1.56^{+0.05}_{-0.04}$
37 | 4410 | 146.58 | 0.47 | -33.4 | 1.61 | 2.1 | $1.39^{+0.03}_{-0.12}$
38 | 4414 | 155.42 | 1.35 | -33.3 | 2.11 | 1.4 | $1.55^{+0.36}_{-0.10}$
39 | 4417 | 157.68 | 1.80 | -33.8 | 0.78 | 1.2 | $1.58^{+0.04}_{-0.06}$
40 | 4427 | 152.97 | -0.85 | -30.1 | 1.66 | 1.9 | $2.60^{+0.03}_{-0.24}$
41 | 4461 | 158.44 | -1.98 | -26.7 | 0.91 | 1.9 | $1.81^{+0.04}_{-0.03}$
42 | 4463 | 148.40 | -1.58 | -25.0 | 1.05 | 1.5 | $1.89^{+0.13}_{-0.09}$
43 | 4476 | 153.15 | -1.55 | -24.3 | 2.08 | 2.5 | $1.76^{+0.06}_{-0.04}$
44 | 4494 | 157.80 | -2.24 | -20.3 | 1.44 | 2.1 | $1.79^{+0.07}_{-0.03}$
45 | 4495 | 159.02 | -2.30 | -20.2 | 1.88 | 2.3 | $2.04^{+0.01}_{-0.08}$
46 | 4496 | 156.89 | -2.35 | -21.5 | 0.55 | 2.0 | $1.51^{+0.04}_{-0.03}$
47 | 5877 | 136.62 | 0.29 | -54.7 | 1.09 | 1.5 | $2.15^{+0.07}_{-0.06}$
48 | 5892 | 145.19 | 2.83 | -58.7 | 0.14 | 1.0 | $4.12^{+0.12}_{-0.98}$
49 | 5909 | 135.63 | 1.81 | -55.9 | 0.08 | 0.9 | $3.78^{+0.17}_{-0.27}$
50 | 5935 | 138.82 | -1.83 | -48.0 | 1.19 | 1.8 | $1.12^{+0.08}_{-0.01}$
51 | 5954 | 137.17 | 3.12 | -51.4 | 0.20 | 1.2 | $3.98^{+0.09}_{-0.38}$
52 | 5960 | 152.34 | -0.70 | -47.7 | 0.16 | 0.9 | $2.30^{+0.06}_{-0.22}$
53 | 5981 | 138.30 | -0.75 | -45.4 | 1.02 | 1.2 | $2.74^{+0.06}_{-0.04}$
54 | 5988 | 140.66 | 0.32 | -45.2 | 0.28 | 1.1 | $1.71^{+0.04}_{-0.10}$
55 | 6000 | 136.18 | 2.12 | -45.4 | 0.48 | 1.1 | $4.32^{+0.04}_{-0.16}$
56 | 6047 | 141.21 | 1.58 | -40.6 | 0.09 | 1.0 | $1.81^{+0.08}_{-0.13}$
57 | 6057 | 147.50 | -2.28 | -39.0 | 0.42 | 1.0 | $1.52^{+0.09}_{-0.27}$
58 | 6070 | 145.96 | 2.62 | -37.5 | 0.47 | 1.1 | $3.94^{+0.18}_{-0.20}$
59 | 6074 | 142.15 | 3.31 | -38.2 | 0.27 | 0.8 | $2.01^{+0.05}_{-0.32}$
60 | 6078 | 143.78 | -3.38 | -35.2 | 0.17 | 4.2 | $1.07^{+0.07}_{-0.10}$
61 | 6097 | 144.88 | 2.34 | -39.5 | 0.12 | 1.1 | $4.17^{+0.08}_{-0.23}$
62 | 6103 | 148.12 | -2.35 | -32.4 | 0.12 | 1.7 | $1.74^{+0.07}_{-0.08}$
63 | 6125 | 148.82 | 2.11 | -32.0 | 0.12 | 1.5 | $1.65^{+0.06}_{-0.22}$
64 | 6131 | 157.30 | 2.88 | -34.7 | 0.08 | 1.2 | $3.35^{+0.33}_{-0.17}$
65 | 6135 | 154.24 | -0.91 | -31.8 | 0.23 | 1.0 | $1.88^{+0.01}_{-0.59}$
66 | 6142 | 158.75 | 0.01 | -30.4 | 0.39 | 0.9 | $0.99^{+0.17}_{-0.09}$
67 | 6143 | 152.87 | 0.22 | -32.8 | 0.14 | 0.8 | $2.43^{+0.16}_{-0.09}$
68 | 6155 | 156.55 | 2.45 | -27.7 | 0.94 | 1.0 | $3.50^{+0.07}_{-0.06}$
69 | 6163 | 155.74 | -2.73 | -27.6 | 0.08 | 1.2 | $1.84^{+0.09}_{-0.08}$
70 | 6164 | 147.01 | -2.44 | -29.8 | 0.09 | 1.2 | $1.10^{+0.10}_{-0.03}$
71 | 6165 | 154.91 | -2.03 | -27.0 | 0.22 | 1.2 | $1.88^{+-0.01}_{-0.25}$
72 | 6170 | 155.76 | -1.22 | -27.1 | 0.30 | 0.8 | $1.52^{+0.03}_{-0.09}$
73 | 6171 | 159.26 | -0.85 | -28.8 | 0.53 | 0.8 | $1.33^{+0.33}_{-0.12}$
74 | 6183 | 150.21 | -2.51 | -25.5 | 0.38 | 1.3 | $1.39^{+0.29}_{-0.10}$
75 | 6185 | 156.00 | -2.41 | -24.5 | 0.09 | 1.1 | $1.46^{+0.14}_{-0.04}$
76 | 6222 | 159.36 | 0.88 | -24.7 | 0.09 | 0.9 | $0.93^{+0.12}_{-0.07}$
77 | 6238 | 147.11 | -2.34 | -22.8 | 0.08 | 1.1 | $1.22^{+0.23}_{-0.10}$
78 | 6240 | 155.49 | -1.88 | -22.9 | 0.19 | 0.9 | $1.55^{+0.09}_{-0.01}$
79 | 9099 | 138.41 | 3.42 | -54.5 | 0.08 | 0.8 | $3.95^{+0.10}_{-0.30}$
80 | 9139 | 153.68 | 1.54 | -32.0 | 0.08 | 0.8 | $1.42^{+0.32}_{-0.17}$
Figure 2: Circles represent CO clouds from MML17. CO clouds shown in orange
represent the ones we investigate, in purple those not under investigation.
Green dots are HMSFRs from Reid14 and Sakai19. Figure 3: Three clouds
demonstrating three main results from the analysis. The left column of panels
shows the on region representing the cloud from MML17 (blue) and the “off”
region in orange. The CO spectra representing these regions are shown in the
second column, with the cloud velocity in black dotted line and gray region
representing Perseus velocities. The third column shows the average reddening
at each distance for the regions in blue and orange, and their difference in
red, with the expected distance to the Perseus arm at the dashed black line.
The top cloud is an example where the method finds a result that is consistent
with expected distance of the arm, the middle cloud is a successful distance
at an inconsistent distance, and the bottom cloud is an example of where the
method fails.  Figure 4: Distribution of 81 MML17 clouds determined by
reddening analysis, shown in orange. Gray regions are not examined by the
method. Dashed green line represents the expected position of the Perseus arm
from Choi et al. (2014). HMSFRs from Reid14 and Sakai19 are shown as green
dots. In blue we show the stars from Xu et al. (2018) that have parallaxes
detected at $>$5 $\sigma$ within this region. A number of features match
between the stars and clouds. The clouds, like the stars, do not seem to bear
any relation to the arm model, except very near the HMSFRs at $l\sim
140^{\circ}$ as expected.  Figure 5: Circles are CO clouds from MML17. Color
bar represents the distance found through reddening the analysis. Green dots
are Reid14 and Sakai19 HMSFRs. The green dashed line is the Perseus arm model
from Choi et al. (2014). Clouds that are very near each other in the
longitude–velocity diagram often have very different distances

##  4 Discussion

Section 3.1 and Figure 1 show that the velocity field toward the Perseus Arm
in the second quadrant of the Milky Way is able to produce the structure-
emulating velocity crowding effect discussed by Burton (1971). From §3.2 and
Figure 4 we can see that this velocity crowding is indeed solely responsibly
for the bridge of molecular clouds that connect star forming regions at
$l\simeq 140^{\circ}$ to those in the anti-center in the longitude–velocity
diagram. From this analysis we find no evidence for any kind of overdensity of
molecular clouds at the distance of the Perseus arm over this range in
longitude. The analyses conducted in Section 3.1 and Section 3.2 produce
consistent results, in that clouds in the Perseus component of the l-v diagram
are actually distributed over a very large range of distances. Indeed, the
clouds generally line up with the sample of O and B stars described in Xu et
al. (2018). This result is consistent with the lack of evidence for a stellar
overdensity in this region in Poggio et al. (2021). Further, the wedge of
molecular clouds we find closer than the purported Perseus arm is somewhat
consistent with the “Cepheus Spur” of young stars highlighted most recently in
Pantaleoni González et al. (2021), although we note the overlap is not
complete. In Figure 5 the impact of velocity crowding is made clear; clouds
with very different distances can be put at very similar places in the
longitude–velocity diagram. The velocity bridge we have interrogated here is,
to our knowledge, the only connection between the star forming region that
ends at $l\simeq 135^{\circ}$ and the cluster near Galactic anti-center,
$l\simeq 180^{\circ}$. Additionally, we show the cloud distance results in
distance–velocity space broken up over five regions of Galactic longitude in
Figure 6. In this figure we include the velocities derived in TP17 to
illustrate both how the final results of this work consistent with TP17 and
how many distances correspond to the same narrow range of velocities over many
sightlines. This figure very clearly shows that it is very difficult to
directly relate a velocity to a distance over this area of sky in the velocity
range we are examining. At velocities above about -30 km s-1 the TP17 curves
show a distance consistently closer than the purported track of the Perseus
arm, with velocities closer to zero typically closer to the Sun, but at lower
velocities there is essentially no consistently measurable relationship
between distance and velocity.

Figure 6: Distance-velocity diagrams for five degree sections of the Galactic
plane. The standard 220 km s-1 flat rotation curve is shown in blue. The model
of the Perseus arm presented in Choi et al. (2014) is shown as a green dot for
the central longitude of the region. The orange lines represent the results of
TP17. These distance–velocity lines are a grid of samples across the longitude
selection and within -2.5${}^{\circ}<b<2.5^{\circ}$, spaced at 0.5∘ intervals.
The purple dots and error bars are the results of our analysis presented in
Table 1.

To understand the scope of problem this effect presents to our disk-mapping
efforts across the whole of the disk we chart the steepness of the
distance–velocity relationship for an assumed flat, 220 km s-1 disk in Figure
7. The red line shows the boundary of 13 km s-1 / kpc, about the steepest
expected gradient in the region where we see a velocity crowding effect
generate the illusion in the longitude–velocity diagram. This map shows that
nearly the entire outer Galaxy and about half of the inner Galaxy are
susceptible to these illusory features given bulk flows of the scale seen
toward Perseus. This provides a serious warning to investigators using the
longitude–velocity diagram to link star-forming features in large swathes of
the Milky Way disk.

Figure 7: A top-view of the Galactic disk, with the sun to the left and
Galactic center at (0, 0). The arc of the Perseus arm defined in Choi et al.
(2014) is shown in green. Regions with a shallower distance–velocity slope
expected from a flat rotation curve than seen toward Perseus shown in color
scale, regions with steeper slope shown in gray scale, separated by red line.
Masers from Reid14 are shown in black. Almost all of the outer disk could be
impacted at a level similar to the region we consider, as could inner-disk
regions near the tangency points and on the far side of the Galaxy. This
diagram implies that most connections we draw between masers using l-v diagram
to form spiral arms are subject to “Burton’s Curse”.

##  5 Conclusion

In this work we have shown that, as originally proposed in B71, real velocity
deviations from a flat Galactic rotation curve can a have significant impact
on our ability to discern spiral structure in the Galactic longitude–velocity
diagram (Figure 1). In particular, we demonstrated that the longitude–velocity
bridge of molecular gas connecting the star forming regions at $l=135^{\circ}$
and $l=180^{\circ}$ is entirely an artifact of these velocity deviations
(Figure 4). The amplitude of these velocity variations can create confusion in
the longitude velocity diagram for the bulk of the Galactic disk (Figure 7).

We believe these results have significant implications for the study of
Galactic structure. On the smaller scale, the Perseus arm star is often used
as a proving ground for theories of star formation. This work shows that one
cannot simply lump gas together at Perseus arm velocities in the second
quadrant and expect to be probing only a spiral arm; gas at the Perseus arm
velocity can extend over at least 3 kpc in radial distance. On larger scales,
we believe this work implies that Galactic spiral arm identification,
especially when many discrete star-forming regions are connected using the
longitude–velocity diagram, requires true 3D evidence. Luckily, recent dust
tomography work (e.g. G19) has been able to probe gas out to many kpc from the
Sun. To extend this work, and reveal the true spiral structure of the Galaxy,
we will need even deeper observations and more sophisticated stellar and dust
inferences schema. We hope this motivates the pursuit of deep and precise
observations by upcoming large-angle imaging surveys, such as those conducted
with the Vera Rubin Observatory and the Nancy Grace Roman Space Telescope.

ACKNOWLEDGEMENTS

This work was supported in part by the US National Science Foundation under
grant 1616177. The authors profusely thank Bob Benjamin for long discussions,
guidance, and wisdom without which this work would not have been possible. The
authors acknowledge Paris-Saclay University’s Institut Pascal programs “The
Self-Organized Star Formation Process” and “The Grand Cascade” and the
Interstellar Institute for hosting discussions that nourished the development
of the ideas behind this work. The authors thank an anonymous referee whose
report has significantly enhanced this work. JEGP thanks Eddie Schlafly for
pointing out the Perseus arm discrepancy in the dust maps to begin with.

##  Appendix A Corrected and augmented version of the MML17 molecular clouds
catalog

The analysis in this paper relies on a corrected version of the molecular
cloud catalog of Miville-Deschênes et al. (2017b).

###  A.1 Astrometry error correction

An error in the coordinates of the cloud catalog published in MML17 has been
identified after publication: GLAT is offset by one beam (7.5 arcmin), and the
central velocity is offset by one channel (1.3 km/s). Correcting this
astrometry error resulted in a slight increase in the number of clouds for
which a kinematic distance could be estimated, from 8107 to 8246 clouds. The
current study uses this corrected version of the catalog.

###  A.2 Extended products

The original catalog provided in MML17 contains only clouds for which a
kinematic distance could be estimated. These clouds correspond to 89% of the
total CO emission, while the whole cloud data set contains 98% of the
emission. To extend the original product we are providing a catalog that
contains the full set of 9710 clouds that encapsulate 98% of the 12CO (J=1-0)
emission. Clouds for which a kinematic distance could not be estimated do not
have mass, density, or physical size estimates.

In addition to this extended cloud catalog, we used the full set of Gaussian
components used to describe the 12CO emission of the Milky Way disk to compute
a series of extended products. These include

  * •

The whole Gaussian data set used to described the CO emission.

  * •

All 9710 reconstructed CO emission position-position-velocity cubes, one per
cloud

  * •

2D maps of the number of clouds and the number of Gaussian on each line of
sight.

  * •

A cube providing the list of Cloud ID at each sky position

  * •

A reconstruction of the whole 12CO data cube, only based on the Gaussian
components.

All these products are publicly available on Dataverse at
https://doi.org/10.7910/DVN/QR9CFW (catalog 10.7910/DVN/QR9CFW).

## References

  * Bolatto et al. (2013) Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, Annual Review of Astronomy and Astrophysics, 51, 207
  * Burton (1971) Burton, W. B. 1971, Astronomy and Astrophysics, 10, 76
  * Burton & Bania (1974) Burton, W. B., & Bania, T. M. 1974, A&A, 33, 425
  * Choi et al. (2014) Choi, Y. K., Hachisuka, K., Reid, M. J., et al. 2014, ApJ, 790, 99
  * Dame et al. (2001) Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, ApJ, 547, 792
  * Efron (1982) Efron, B. 1982, The jackknife, the bootstrap and other resampling plans, Vol. 38 (Society for industrial and applied mathematics Philadelphia)
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1
  * Gaia Collaboration et al. (2021) —. 2021, A&A, 649, A1
  * Green et al. (2019) Green, G. M., Schlafly, E., Zucker, C., Speagle, J. S., & Finkbeiner, D. 2019, ApJ, 887, 93
  * Hall (1985) Hall, P. 1985, Stochastic Processes and their Applications, 20, 231
  * Humphreys (1976) Humphreys, R. M. 1976, ApJ, 206, 114
  * Kalberla et al. (2005) Kalberla, P. M. W., Burton, W. B., Hartmann, D., et al. 2005, A&A, 440, 775
  * Kalberla & Kerp (2009) Kalberla, P. M. W., & Kerp, J. 2009, ARA&A, 47, 27
  * Levine et al. (2006) Levine, E. S., Blitz, L., & Heiles, C. 2006, ApJ, 643, 881
  * Lin et al. (1969) Lin, C. C., Yuan, C., & Shu, F. H. 1969, ApJ, 155, 721
  * M. Green (2018) M. Green, G. 2018, The Journal of Open Source Software, 3, 695
  * Majewski et al. (2017) Majewski, S. R., Schiavon, R. P., Frinchaboy, P. M., et al. 2017, AJ, 154, 94
  * Miville-Deschênes et al. (2017a) Miville-Deschênes, M.-A., Murray, N., & Lee, E. J. 2017a, ApJ, 834, 57
  * Miville-Deschênes et al. (2017b) Miville-Deschênes, M.-A., Salomé, Q., Martin, P. G., et al. 2017b, A&A, 599, A109
  * Pantaleoni González et al. (2021) Pantaleoni González, M., Maíz Apellániz, J., Barbá, R. H., & Reed, B. C. 2021, MNRAS, 504, 2968
  * Poggio et al. (2021) Poggio, E., Drimmel, R., Cantat-Gaudin, T., et al. 2021, arXiv, arXiv:2103.01970
  * Reid et al. (2014) Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2014, ApJ, 783, 130
  * Roberts (1972) Roberts, W. W. J. 1972, ApJ, 173, 259
  * Sakai et al. (2019) Sakai, N., Reid, M. J., Menten, K. M., Brunthaler, A., & Dame, T. M. 2019, ApJ, 876, 30
  * Schlafly et al. (2014) Schlafly, E. F., Green, G., Finkbeiner, D. P., et al. 2014, ApJ, 789, 15
  * Sellwood & Carlberg (1984) Sellwood, J. A., & Carlberg, R. G. 1984, ApJ, 282, 61
  * Sellwood & Masters (2021) Sellwood, J. A., & Masters, K. L. 2021, arXiv e-prints, arXiv:2110.05615
  * Sparke & Dodd (1978) Sparke, L. S., & Dodd, R. J. 1978, MNRAS, 182, 1
  * Tchernyshyov & Peek (2017) Tchernyshyov, K., & Peek, J. E. G. 2017, AJ, 153, 8
  * Tchernyshyov et al. (2018) Tchernyshyov, K., Peek, J. E. G., & Zasowski, G. 2018, AJ, 156, 248
  * Toomre (1981) Toomre, A. 1981, Structure and Evolution of Normal Galaxies, 111
  * van de Hulst et al. (1954) van de Hulst, H. C., Muller, C. A., & Oort, J. H. 1954, Bulletin of the Astronomical Institutes of the Netherlands, 12, 117
  * Wouterloot & Habing (1985) Wouterloot, J. G. A., & Habing, H. J. 1985, A&A, 151, 297
  * Xu et al. (2018) Xu, Y., Bian, S. B., Reid, M. J., et al. 2018, A&A, L15
  * Zasowski et al. (2015) Zasowski, G., Menard, B., Bizyaev, D., et al. 2015, ApJ, 798, 35
  * Zucker et al. (2020) Zucker, C., Speagle, J. S., Schlafly, E. F., et al. 2020, A&A, 633, A51


# Separating the impact of nuclear skin and nuclear deformation on elliptic
flow and its fluctuations in high-energy isobar collisions

Jiangyong Jia<EMAIL_ADDRESS>Department of Chemistry, Stony
Brook University, Stony Brook, NY 11794, USA  Physics Department, Brookhaven
National Laboratory, Upton, NY 11976, USA  Giuliano Giacalone
<EMAIL_ADDRESS>Institut für Theoretische Physik,
Universität Heidelberg, Philosophenweg 16, 69120 Heidelberg, Germany  Chunjian
Zhang  Department of Chemistry, Stony Brook University, Stony Brook, NY 11794,
USA

###### Abstract

Bulk nuclear structure properties, such as radii and deformations, leave
distinct signatures in the final states of relativistic heavy-ion collisions.
Collisions of isobars, in particular, offer an easy route to establish clear
correspondences between the structure of colliding nuclei and the final state
observables. Here we investigate the impact of nuclear skin and nuclear
deformations on elliptic flow ($v_{2}$) and its fluctuations in high-energy
96Ru+96Ru and 96Zr+96Zr collisions, for which experimental data is available.
We show that the difference in skin thickness between these isobars impacts
the intrinsic ellipticity of the collision systems, or reaction-plane flow,
$v_{2}^{\mathrm{rp}}$. In contrast, differences in nuclear deformations impact
only the fluctuations of $v_{2}$ around $v_{2}^{\mathrm{rp}}$. Through isobar
collisions, one can separate the influence of nuclear skin and nuclear
deformations in the $v_{2}$ data. This is a significant step towards assessing
the consistency of nuclear phenomena across energy scales.

###### pacs:

25.75.Gz, 25.75.Ld, 25.75.-1

The structure of most of atomic nuclei is characterized by prominent bulk
properties, reflecting collective correlations in many-body systems held
together by the strong force. Mapping out such properties and how they evolve
across the Segré chart is one of the main goals of nuclear physics Nakatsukasa
_et al._ (2016); Nazarewicz (2016). Information about the collective features
of nuclei is typically inferred via spectroscopic and scattering experiments
conducted at low energies. In the past five years or so, an expanded list of
nuclei utilized in collision experiments at ultra-relativistic energies has
lead to the identification of fingerprints of collective nuclear properties
also in such processes Shou  _et al._ (2015); Goldschmidt  _et al._ (2015);
Giacalone  _et al._ (2018); Giacalone (2020a, b); Giacalone  _et al._ (2021a,
b); Jia  _et al._ (2022a); Jia (2022a); Bally  _et al._ (2022); Jia (2022b).
The angular (azimuthal) particle distributions emitted in high-energy nuclear
collisions carry, in particular, direct information about the structure of the
colliding ions, to the extent that, for certain collision configurations, the
shape of the particle distributions in momentum space can be related directly
to the shape of the colliding ions at the time of scattering.

This is made possible by the nearly-ideal fluid nature of the quark-gluon
plasma (QGP) created in high-energy collisions. In a hydrodynamic picture, the
emergence of final-state anisotropies in the azimuthal particle spectra stems
from the presence of spatial anisotropies in the initial conditions of the
fluid expansion Gale  _et al._ (2013); Heinz and Snellings (2013); Romatschke
and Romatschke (2019). These anisotropies are in turn sourced by the positions
of the nucleons that populate the colliding ions at the time of interaction.
Contrary to low-energy scattering experiments, where one can only access the
average nuclear charge distributions, high-energy heavy-ion collisions grant
access to the spatial positions of all nucleons on an event-by-event basis,
and consequently probe multi-nucleon correlations in the nuclear ground
states. These correlations can then be measured experimentally using so-called
multi-particle correlation techniques Jia (2022a, b). The important question
is, then, to which extent established information from low-energy nuclear
physics can provide a consistent picture of the phenomena observed at high-
energy colliders. In this paper we make significant progress in this
direction.

In heavy-ion collisions, the deformation of the colliding ions becomes
manifest mainly in the limit of fully-overlapping configurations (central
collisions), where the shape of the colliding bodies is fully resolved.
However, nuclear structure effects can be cleanly isolated as well over the
full centrality range by comparing two isobaric collision systems Giacalone
_et al._ (2021a); Jia (2022a, b). Isobar nuclei have the same mass number,
therefore, any visible difference in the observables measured in two isobaric
systems must originate from differences in their structure, which impacts the
initial condition and evolution of the QGP. Collisions of isobars, 96Ru+96Ru
and 96Zr+96Zr collisions, have been performed at the BNL Relativistic Heavy
Ion Collider (RHIC). They demonstrate this argument. Ratios of observables
taken between 96Ru+96Ru and 96Zr+96Zr collisions show significant and
centrality-dependent departures from unity Abdallah  _et al._ (2022). Most
models of heavy ion collisions describe the data by parameterizing the nucleon
density within the colliding ions according to a Woods-Saxon (WS) profile,

$\displaystyle\rho(r,\theta,\phi)$
$\displaystyle\propto\frac{1}{1+e^{[r-R_{0}\left(1+\beta_{2}Y_{2}^{0}(\theta,\phi)+\beta_{3}Y_{3}^{0}(\theta,\phi)\right)]/a_{0}}},$
(1)

which contains four structure parameters, nuclear skin $a_{0}$, half-width
radius $R_{0}$, quadrupole deformation $\beta_{2}$, and octupole deformation
$\beta_{3}$. Model studies have established that isobar ratios are indeed
controlled by parameter differences,
$\Delta\beta_{2}^{2}=\beta_{\mathrm{2Ru}}^{2}-\beta_{\mathrm{2Zr}}^{2}$ ,
$\Delta\beta_{3}^{2}=\beta_{\mathrm{3Ru}}^{2}-\beta_{\mathrm{3Zr}}^{2}$,
$\Delta a_{0}=a_{\mathrm{0Ru}}-a_{\mathrm{0Zr}}$ and $\Delta
R_{0}=R_{\mathrm{0Ru}}-R_{\mathrm{0Zr}}$ Zhang and Jia (2022).

Many observables have been shown to present a sensitivity to the nuclear
profile parameters, such as the mean transverse momentum $p_{\mathrm{T}}$ Xu
_et al._ (2021), its fluctuations Jia (2022b), the spectator neutron number
Liu  _et al._ (2022), Fourier plane correlations Zhao  _et al._ (2022); Jia
_et al._ (2022b), and shape-size correlations Giacalone (2020a); Bally  _et
al._ (2022); Jia  _et al._ (2022a). In the present study, the key observable
is the simple elliptical asymmetry of the azimuthal particle distributions,
known as elliptic flow $V_{2}=v_{2}e^{2i\Psi}$, which characterizes the
anisotropic flow of particles along the direction $\Psi$ with an amplitude
$v_{2}$ and periodicity of $\pi$. $V_{2}$ emerges as a hydrodynamic response
to the elliptical shape of the region of overlap between two colliding ions.
In isobar collisions, the ratio of $v_{2}$ between 96Ru+96Ru and 96Zr+96Zr
collisions shows a complex non-monotonic centrality dependence, which can be
explained as a combined effect from differences in the WS parameters between
96Ru and 96Zr Jia and Zhang (2021). In this paper, we show that the impact of
the deformations parameters ($\beta_{2}$ and $\beta_{3}$) can be fully
disentangled from that of the radial profile parameters ($a_{0}$ and $R_{0}$),
and point out the consequences of such a result.

We start from Fig. 1. We parametrize the plane transverse to the collision
axis (transverse plane) with Cartesian coordinates where the $x$ direction is
along the impact parameter direction. For events at a given centrality, the
distribution of $V_{2}\equiv(v_{2x},v_{2y})$ is approximately a two-
dimensional Gaussian Voloshin  _et al._ (2008)

$\displaystyle
p(v_{2\mathrm{x}},v_{2\mathrm{y}})=\frac{1}{\pi\delta^{2}}\exp\biggl{[}-\frac{(v_{2\mathrm{x}}-v_{2}^{\mathrm{rp}})^{2}+v_{2\mathrm{y}}^{2}}{\delta^{2}}\biggr{]}.$
(2)

The displacement along $x$, $v_{2}^{\mathrm{rp}}$, is the so-called reaction
plane flow, associated with the average elliptic geometry, whereas the
fluctuation around the intrinsic geometry, $\delta$, is the variance of
elliptic flow due to, e.g., fluctuations in the positions of the colliding
nucleons. Our point is that, to leading order, a change in the radial profile
of the nucleus, determined by either $a_{0}$ or $R_{0}$ [see Fig. 1(a)],
modifies the intrinsic ellipticity, $v_{2}^{\mathrm{rp}}$, with little impact
on the flow fluctuation. On the other hand, in the presence of nuclear
deformations [see Fig. 1(b)], the random orientation of the colliding nuclei
implies an increase in $\delta$, with little impact on $v_{2}^{\mathrm{rp}}$
111Note that a change in nuclear structure also impacts the distribution
$p(N_{\mathrm{part}})$, such that the events with the same $N_{\mathrm{part}}$
correspond to slightly different centrality and vice versa Jia  _et al._
(2022c). This secondary effect leads to a small correlation between
$\varepsilon_{2}^{\mathrm{rp}}$ and the amount of deformation, but this is a
subleading effect compared to that discussed in Fig. 1.. Now, if $p(V_{2})$ is
a Gaussian, the root-mean-squared elliptic flow, experimentally accessible via
a two-particle correlation, is
$v_{2}\\{2\\}=\sqrt{(v_{2}^{\mathrm{rp}})^{2}+\delta^{2}}$, while higher-order
cumulants of $v_{2}$, measured via multi-particle correlations, are all
identical and equal to $v_{2}^{\mathrm{rp}}$,

$v_{2}\\{4\\}=v_{2}\\{6\\}=\ldots=v_{2}\\{\infty\\}=v_{2}^{{\mathrm{rp}}}.$
(3)

In this limit, the fluctuation of $v_{2}$ can be measured as

$\delta^{2}=v_{2}\\{2\\}^{2}-v_{2}\\{4\\}^{2}.$ (4)

In the following, we demonstrate our argument about the sensitivity of
$v_{2}^{\mathrm{rp}}$ and $\delta$ to the nuclear structure parameters by
means of transport model calculations.

Figure 1:  (a) Sketch of a collision of spherical nuclei with different skin
thickness, $a_{0}$. The influence of the nuclear skin is primarily on the
average ellipticity along the $x$ direction, i.e., an increase in the reaction
plane flow, $v_{2}^{\mathrm{rp}}$. (b) Collisions of deformed nuclei with
$\beta_{2}=0.25$ and random orientations (four orientations for each nucleus
as shown by solid lines). The random orientations mostly result in an enhanced
ellipticity fluctuation, $\delta$.

We simulate the dynamics of the QGP using the multi-phase transport model
(AMPT) Lin  _et al._ (2005). We use AMPT v2.26t5 in the string-melting mode at
$\mbox{$\sqrt{s_{\mathrm{NN}}}$}=200$ GeV with a partonic cross section of 3.0
$m$b Ma and Bzdak (2014); Bzdak and Ma (2014). This model successfully
describes the isobar ratios of $v_{2}$, $v_{3}$ and $N_{\mathrm{ch}}$ measured
by the STAR collaboration Giacalone  _et al._ (2021b); Zhang and Jia (2022).
We simulate generic isobar 96X+96X collisions with five choices of
$\beta_{2}$, $\beta_{3}$, $R_{0}$ and $a_{0}$, as listed in Table 1. This
allows us to define ratios that isolate the influence of the nuclear structure
parameters step-by-step, e.g., Case1/Case2 includes the effect of $\beta_{2}$,
Case1/Case3 includes the effect of $\beta_{2}$ and $\beta_{3}$, and so on. The
cumulants of elliptic flow are calculated within the multi-particle cumulant
framework Bilandzic  _et al._ (2011, 2014), for hadrons with
$0.2<p_{\mathrm{T}}<2$ GeV. The two-particle cumulant $v_{2}\\{2\\}$ is
obtained by correlating particles in $0<\eta<2$ with those in $-2<\eta<0$ to
suppress short-range correlations that do not emerge from the collective
expansion of the system Jia  _et al._ (2017). $v_{2}\\{4\\}$ is instead free
of such contributions by construction, and is thus obtained from all particles
with $|\eta|<2$. Additionally, we calculate the true $v_{2}^{\mathrm{rp}}$
from the azimuthal correlation of particles relative to the impact parameter,
and the true flow fluctuation as
$\delta_{\mathrm{rp}}^{2}=v_{2}\\{2\\}^{2}-(v_{2}^{\mathrm{rp}})^{2}\;$. The
simulated events are binned in classes defined by the number of nucleons that
participate in the interaction $N_{\mathrm{part}}$, which is a good enough
proxy for the experimentally-defined collision centrality.

|  $R_{0}$ (fm) |  $a_{0}$ (fm) | $\beta_{2}$ | $\beta_{3}$
---|---|---|---|---
Case1 96Ru | 5.09 | 0.46 | 0.162 | 0
Case2 | 5.09 | 0.46 | 0.06 | 0
Case3 | 5.09 | 0.46 | 0.06 | 0.20
Case4 | 5.09 | 0.52 | 0.06 | 0.20
Case5 96Zr | 5.02 | 0.52 | 0.06 | 0.20
Ratios |  $\frac{\mathrm{\textstyle\small Case1}}{\mathrm{\textstyle\small Case2}}$ , $\frac{\mathrm{\textstyle\small Case1}}{\mathrm{\textstyle\small Case3}}$ , $\frac{\mathrm{\textstyle\small Case1}}{\mathrm{\textstyle\small Case4}}$ , $\frac{\mathrm{\textstyle\small Case1}}{\mathrm{\textstyle\small Case5}}$
Table 1:  Nuclear structure parameters used in the simulations of 96Ru+96Ru
and 96Zr+96Zr collisions. Case1 and Case5 represent, respectively, our full
parameterizations of 96Ru and 96Zr.

In Fig. 2, we show our results for $v_{2}\\{2\\}$, $v_{2}\\{4\\}$, and
$\delta$, averaged between 96Ru+96Ru and 96Zr+96Zr collisions, and how they
compare with STAR data. The comparison with data is generally good. In central
collisions, AMPT correctly predicts a sign change for $v_{2}\\{4\\}$,
qualitatively similar to STAR data but with a smaller magnitude. This sign
change can be attributed to volume fluctuations associated with the event-
activity variable, such as $N_{\mathrm{part}}$ or $N_{\mathrm{ch}}$, which
could have a strong impact on multi-particle cumulants Zhou and Jia (2018);
Jia  _et al._ (2020, 2022a). Indeed, we observe a stronger sign change pattern
when the events are binned in $N_{\mathrm{ch}}$ 222The $N_{\mathrm{ch}}$ is
defined as number of charged particles with $p_{\mathrm{T}}>0.1$ GeV and
$|\eta|<0.5$ similar to STAR analysis Abdallah  _et al._ (2022)., which can be
naturally attributed to the fact that $N_{\mathrm{part}}$ has a better
centrality resolution than $N_{\mathrm{ch}}$. In peripheral collisions, the
model underpredicts significantly the value of $v_{2}\\{4\\}$, while correctly
capturing the measured $\delta$. This implies that AMPT has a too small
$v_{2}^{\mathrm{rp}}$. Probably, this is due to the fact that particle
production in AMPT scales like $N_{\mathrm{part}}$, which is known to lead to
significantly smaller $v_{2}^{\mathrm{rp}}$ compared to other models in closer
agreement to the energy deposition scaling expected in high-energy quantum
chromodynamics Giacalone  _et al._ (2017a). A recent calculation of
$v_{2}\\{2\\}$ in isobar collisions by Nijs and van der Schee Nijs and van der
Schee (2021), for instance, starting from the
T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo model Moreland  _et al._ (2015), does
not have this problem.

Figure 2:  Values of $v_{2}\\{2\\}$ (a) $v_{2}\\{4\\}$ and
$v_{2}^{{\mathrm{rp}}}$ (b), $\delta$ and $\delta_{\mathrm{rp}}$ (c) as a
function of $N_{\mathrm{part}}$, averaged between Ru+Ru and Zr+Zr collisions.
The AMPT results are compared to the STAR data taken from Fig. 23 in Ref.
Abdallah  _et al._ (2022). For the STAR data, we also approximate
$v_{2}^{{\mathrm{rp}}}$ by $v_{2}\\{\mathrm{zdc}\\}$, and
$\delta_{\mathrm{rp}}$ by $\delta_{\mathrm{zdc}}$, as discussed in the text.

Isobar ratios represent nevertheless more robust predictions than single-
system estimates. We move on and calculate the ratios of these quantities to
assess the impact of the nuclear structure. For an observable, $\mathcal{O}$,
the ratio is calculated at a given $N_{\mathrm{part}}$ as

$\displaystyle
R_{\mathcal{O}}(N_{\mathrm{part}})=\frac{\mathcal{O}_{\mathrm{Ru}}(N_{\mathrm{part}})}{\mathcal{O}_{\mathrm{Zr}}(N_{\mathrm{part}})}\;.$
(5)

The results are shown in the left column of Fig. 3. Figure 3(a) shows the
complex centrality dependence of $R_{v_{2}\\{2\\}}$, arising from both
deformation and radial profile parameters. In contrast, $R_{v_{2}\\{4\\}}$
[Fig. 3(b)] is sensitive mostly to $a_{0}$, whereas $R_{\delta}$ [Fig. 3(c)]
is sensitive mostly to $\beta_{n}$. Therefore, the behavior of
$R_{v_{2}\\{2\\}}$ can be fully decomposed into a part sensitive to the
nuclear skin and a part sensitive to the nuclear deformations. One can
establish the following relation,

$\displaystyle
R_{v_{2}\\{2\\}}^{2}=R_{\delta}^{2}+(R_{v_{2}\\{4\\}}^{2}-R_{\delta}^{2})r\;,\;r=v_{2}\\{4\\}^{2}/v_{2}\\{2\\}^{2}$
$\displaystyle R_{v_{2}\\{2\\}}\approx
R_{\delta}+(R_{v_{2}\\{4\\}}-R_{\delta})r\;,$ (6)

where the second line is obtained by assuming all ratios are close to unity.
In central collisions where $r\approx 0$, the behavior of $R_{v_{2}\\{2\\}}$
is dominated by $R_{\delta}$. In mid-central collisions, the non-monotonic
behavior of $R_{v_{2}\\{2\\}}$ in the top-left panel of Fig. 3 results from
the interplay of $R_{v_{2}\\{4\\}}$ and $R_{\delta}$. This is our main result.
The right column of Fig. 3 shows, in addition, the ratios of the true
intrinsic ellipticities, $R_{v_{2}^{\mathrm{rp}}}$ [Fig. 3(d)], and the true
flow fluctuation, $R_{\delta_{\mathrm{rp}}}$ [Fig. 3(e)], obtained from
$\delta_{\mathrm{rp}}^{2}=v_{2}\\{2\\}^{2}-(v_{2}^{\mathrm{rp}})^{2}$. We see
that, for a difference in skin thickness of
$a_{\mathrm{0Ru}}-a_{\mathrm{0Zr}}=0.06$ fm, the value of
$v_{2}^{{\mathrm{rp}}}$ is enhanced by about 10% in 96Ru+96Ru collisions. The
impact of $\beta_{n}$ on $R_{v_{2}^{\mathrm{rp}}}$ is subleading, as expected.
In addition, we see that the value of $v_{2}^{\mathrm{rp}}$ varies more
strongly when changing nuclear structure parameters than that of
$v_{2}\\{4\\}$. As a result, the ratio $R_{\delta_{\mathrm{rp}}}$ shows
somewhat larger dependence on nuclear structure than $R_{\delta}$.

One further comment is in order. Elliptic flow emerges, event-by-event, as a
response to the initial ellipticity of the system. This quantity, denoted by
$\mathcal{E}_{2}$, is usually quantified as the second Fourier harmonic of the
initial density distribution Teaney and Yan (2011). With the excellent
approximation of a linear scaling, $V_{2}\propto\mathcal{E}_{2}$, on an event-
by-event basis, the ratios of observables analyzed in Fig. 3 can be estimated
starting solely from the knowledge of $\mathcal{E}_{2}$ and its fluctuations.
In the supplemental material we show that the behaviors observed in Fig. 3 do
indeed originate largely from the initial state.

Figure 3:  Isobar ratios $R_{v_{2}\\{2\\}}$ (a), $R_{v_{2}\\{4\\}}$ (b),
$R_{\delta}$ (c), $R_{v_{2}^{\mathrm{rp}}}$ (d), and
$R_{\delta_{\mathrm{rp}}}$ (e) plotted as a function of $N_{\mathrm{part}}$.
They are compared with STAR data from Fig. 23 of Ref. Abdallah  _et al._
(2022). These ratios are re-calculated at matching $N_{\mathrm{part}}$, which
are slightly different from the published ratios at matching centrality as
described in detail in Ref. Jia  _et al._ (2022c).

Before concluding, let us comment on the fact that the STAR collaboration has
also measured an approximation of $v_{2}^{\mathrm{rp}}$ by correlating
particles with the spectator neutrons in the zero-degree calorimeters (ZDC).
In the previous figure Fig. 2(b), the resulting elliptic flow,
$v_{2}\\{\mathrm{zdc}\\}$, is compared to the measured $v_{2}\\{4\\}$. The
$v_{2}\\{\mathrm{zdc}\\}$ is smaller than $v_{2}\\{4\\}$ in peripheral
collisions, but is above it toward central collisions. The overall centrality
dependent trends is similar to AMPT $v_{2}^{\mathrm{rp}}$, showing that
$v_{2}\\{{\mathrm{zdc}}\\}$ is likely a good proxy for $v_{2}^{\mathrm{rp}}$
at least in peripheral collisions. In Fig. 2(c), we also define the
corresponding
$\delta_{\mathrm{zdc}}=\sqrt{v_{2}\\{2\\}^{2}-(v_{2}\\{\mathrm{zdc}\\})^{2}}$,
and compare the results to the measured $\delta$. They are close to each other
in peripheral collisions. In Figure 3, then, we see that
$R_{v_{2}\\{\mathrm{zdc}\\}}$ agrees with $R_{v_{2}^{\mathrm{rp}}}$ in
peripheral collisions [Fig. 3(d)], but gradually deviates and looses
sensitivity to WS parameters towards central collisions. This could be due to
a strong decorrelation between the spectator plane and the reaction plane
ALICE Collaboration (2022) occurring when the number of spectator neutrons is
small. In contrast, $R_{\delta_{\mathrm{zdc}}}$ shows a good agreement with
$R_{\delta_{\mathrm{rp}}}$ [Fig. 3(d)].

Wrapping up, we have found that the nuclear radial profile parameters, i.e.,
nuclear skin thickness, $a_{0}$, and half-density radius, $R_{0}$, mostly
influence the magnitude of $v_{2}$ along the impact parameter direction, i.e.,
the reaction-plane elliptic flow, $v_{2}^{\mathrm{rp}}$. Conversely, the
nuclear deformations, $\beta_{n}$, mostly influence the elliptic flow
fluctuation, $\delta$. Based on our simulations, we conclude that the measured
isobar ratio of $v_{2}\\{4\\}$ is determined by
$a_{\mathrm{0Ru}}-a_{\mathrm{0Zr}}$, while the measured isobar ratio of
$\delta$ arises from the interplay of
$\beta_{\mathrm{2Ru}}^{2}-\beta_{\mathrm{2Zr}}^{2}$ and
$\beta_{\mathrm{3Ru}}^{2}-\beta_{\mathrm{3Zr}}^{2}$. This new result has now
to be combined with the information from previous studies Zhang and Jia
(2022); Nijs and van der Schee (2021), showing that the isobar ratio of
triangular flow coefficients, $v_{3}\\{2\\}$, is dominated by
$\beta_{\mathrm{3Ru}}^{2}-\beta_{\mathrm{3Zr}}^{2}$. Therefore, We have three
ratios, $R_{v_{3}\\{2\\}}$, $R_{v_{2}\\{4\\}}$, and $R_{\delta}$, that provide
separate constraints on three features of the colliding nuclei, namely $\Delta
a_{0}$, $\Delta\beta_{2}^{2}$, and $\Delta\beta_{3}^{2}$. A pressing question
in nuclear phenomenology is whether the structure-induced phenomena observed
in high-energy collisions are consistent or not with the expectations from
low-energy nuclear experiments and theories. These results represent, then, a
stepping stone to the achievement of this goal. Finally, let us briefly
comment on the fact that our conclusions do not make any specific use of the
Gaussian Ansatz for $v_{2}$ fluctuations that we have used to motivate our
analysis. The fluctuations of $v_{2}$ are, in fact, non-Gaussian, especially
in peripheral collisions, where $v_{2}^{\mathrm{rp}}$ is large and one becomes
sensitive to the bound $v_{2}<1$ Giacalone  _et al._ (2017b); Bhalerao  _et
al._ (2019). It would be interesting to generalize this study to higher-order
cumulants, $v_{2}\\{4,6,8\\}$, and see how nuclear structure impacts these
quantities in isobar collisions (see Ref. Giacalone (2019) for a study in
collisions of uranium-238 nuclei). In the supplemental material, we show
results for $R_{v_{2}\\{4\\}}$, $R_{v_{2}\\{6\\}}$ and $R_{v_{2}\\{8\\}}$, and
also study the fine splitting of these cumulants at the level of the
eccentricity fluctuations. Our preliminary conclusion is that, for such
observables, there is no obvious separation of nuclear structure effects, as
in the case of $v_{2}\\{2\\}$, although this is work in progress.

Acknowledgements: This research of J.J and C.Z is supported by DOE DE-
FG02-87ER40331. The research of G.G. is funded by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germanys
Excellence Strategy EXC2181/1-390900948 (the Heidelberg STRUCTURES Excellence
Cluster), within the Collaborative Research Center SFB1225 (ISOQUANT, Project-
ID 273811115). We acknowledge Somadutta Bhatta and Jean-Yves Ollitrault for
useful discussions.

## Supplemental material

A similar analysis is performed for the eccentricity, $\varepsilon_{2}$, and
its fluctuations, where we decompose
$\varepsilon_{2}\\{2\\}\equiv\sqrt{\left\langle\varepsilon_{2}^{2}\right\rangle}$
into a reaction plane component and a fluctuation component. The results are
shown in Fig. 4 with identical layout as in Fig. 3. Specifically, we calculate
the four-particle cumulant $\varepsilon_{2}\\{4\\}$ and the standard reaction
plane eccentricity Miller  _et al._ (2007), $\varepsilon_{2}^{\mathrm{rp}}$,
as well as the associated fluctuations defined by
$\delta_{\varepsilon_{2}}\equiv\sqrt{\varepsilon_{2}\\{2\\}^{2}-\varepsilon_{2}\\{4\\}^{2}}$
and
$\delta_{\varepsilon_{2}^{\mathrm{rp}}}\equiv\sqrt{\varepsilon_{2}\\{2\\}^{2}-(\varepsilon_{2}^{\mathrm{rp}})^{2}}$
to match the corresponding final state quantities in Fig. 3. We see that the
isobar ratios of these initial state estimators already qualitatively or even
quantitatively reproduce the AMPT result in Fig. 3, with only a few
exceptions. Compared to Fig. 3, we see in particular that the values of
$R_{\varepsilon_{2}\\{4\\}}$ are systematically smaller than
$R_{v_{2}\\{4\\}}$. However, the values of $R_{\varepsilon_{2}^{\mathrm{rp}}}$
agree quantitatively with $R_{v_{2}^{\mathrm{rp}}}$ in all four cases. The
stronger dependence of $R_{\delta_{\varepsilon}}$ on WS parameters, compared
to $R_{\delta}$, may then be attributed to the role of non-Gaussianities,
which are known to be larger in the distribution of $\varepsilon_{2}$ than in
that of $v_{2}$, as $v_{2}$ fluctuations are largely smeared by the
hydrodynamic expansion. We also see that
$R_{\delta_{\varepsilon_{2}^{\mathrm{rp}}}}$ agrees rather nicely with
$R_{\delta_{\mathrm{rp}}}$ in Fig. 3, which is expected as both
$\varepsilon_{2}^{\mathrm{rp}}$ and $v_{2}^{\mathrm{rp}}$ are defined relative
to the impact parameter, and therefore is less affected by non-Gaussianities.
However, both $R_{\delta_{\varepsilon_{2}}}$ and
$R_{\delta_{\varepsilon^{\mathrm{rp}}}}$ show some dependence on $a_{0}$ in
peripheral collisions, which are absent in $R_{\delta}$ and
$R_{\delta_{\mathrm{rp}}}$ in Fig. 3. This is once more likely due to a
smearing effect from the hydrodynamic expansion, washing out the primordial
non-Gaussianities.

Figure 4:  The isobar ratios $R_{\varepsilon_{2}\\{2\\}}$ (a),
$R_{\varepsilon_{2}\\{4\\}}$ (b), $R_{\delta_{\varepsilon}}$ (c),
$R_{\varepsilon_{2}^{\mathrm{rp}}}$ (e), and
$R_{\delta_{\varepsilon^{\mathrm{rp}}}}$ (f) plotted as a function of
$N_{\mathrm{part}}$. They are compared with the same STAR data shown in Fig.
3.

Moving to higher-order fluctuations, we also calculate the ratios of higher-
order cumulants, $R_{v_{2}\\{6\\}}$ and $R_{v_{2}\\{8\\}}$, and compare them
to $R_{v_{2}\\{4\\}}$. In the Gaussian limit, these ratios should all be
identical. However, a characteristic fine splitting is observed
experimentally, $v_{2}\\{4\\}\gtrsim v_{2}\\{6\\}\gtrsim v_{2}\\{8\\}$ Aad
_et al._ (2013); Sirunyan  _et al._ (2019); Acharya  _et al._ (2018); Aaboud
_et al._ (2020), reflecting the non-Gaussianity of the distribution of
$\varepsilon_{2}$ Giacalone  _et al._ (2017b); Bhalerao  _et al._ (2019). It
is interesting to study whether this fine splitting is affected by the nuclear
structure differences. The results are shown in Fig. 5. Unfortunately, the
statistical precision of our AMPT results does not allow for a definitive
answer to this question. On the other hand, the ratios of higher-order
cumulants for $\varepsilon_{2}$ can be calculated with high precision to
provide a useful guidance. The results, shown in Fig. 6, bring in a sense both
good and bad news. The good news is that the cumulant splittings are impacted
by the nuclear structure parameters. The bad news is that both $\beta_{n}$ and
$a_{0}$ seem to affect these ratios at similar level, by reducing the cumulant
ratios $\varepsilon_{2}\\{6\\}/\varepsilon_{2}\\{4\\}$ and
$\varepsilon_{2}\\{8\\}/\varepsilon_{2}\\{6\\}$. We note that the reduction is
larger for increasing $N_{\mathrm{part}}$, and may even change the overall
trends of these ratios. While more work may be needed on the conceptual side,
it would still be interesting to study whether these effects also survive to
the final state and leave similar imprints on $v_{2}\\{4,6,8\\}$, both in
simulations and in experiment.

Figure 5:  The isobar ratios $R_{v_{2}\\{4\\}}$, $R_{v_{2}\\{6\\}}$ and
$R_{v_{2}\\{8\\}}$ (left column) calculated at matching $N_{\mathrm{part}}$
and plotted as a function of $N_{\mathrm{part}}$. The right column show the
ratios calculated based on eccentricity $\varepsilon_{2}$,
$R_{\varepsilon_{2}\\{4\\}}$, $R_{\varepsilon_{2}\\{6\\}}$ and
$R_{\varepsilon_{2}\\{8\\}}$. Figure 6:  The ratio of higher-order
eccentricity cumulants $\varepsilon_{2}\\{6\\}/\varepsilon_{2}\\{4\\}$ (a) and
$\varepsilon_{2}\\{8\\}/\varepsilon_{2}\\{6\\}$ (b), and corresponding double
ratios, $R_{\varepsilon_{2}\\{6\\}/\varepsilon_{2}\\{4\\}}$ (c) and
$R_{\varepsilon_{2}\\{8\\}/\varepsilon_{2}\\{6\\}}$ (d) to quantify the
influence of nuclear structure parameters.

## References

  * Nakatsukasa  _et al._ (2016) Takashi Nakatsukasa, Kenichi Matsuyanagi, Masayuki Matsuo, and Kazuhiro Yabana, “Time-dependent density-functional description of nuclear dynamics,” Rev. Mod. Phys. 88, 045004 (2016), arXiv:1606.04717 [nucl-th] .
  * Nazarewicz (2016) Witold Nazarewicz, “Challenges in Nuclear Structure Theory,” J. Phys. G 43, 044002 (2016), arXiv:1603.02490 [nucl-th] .
  * Shou  _et al._ (2015) Q. Y. Shou, Y. G. Ma, P. Sorensen, A. H. Tang, F. Videbæk, and H. Wang, “Parameterization of Deformed Nuclei for Glauber Modeling in Relativistic Heavy Ion Collisions,” Phys. Lett. B 749, 215–220 (2015), arXiv:1409.8375 [nucl-th] .
  * Goldschmidt  _et al._ (2015) Andy Goldschmidt, Zhi Qiu, Chun Shen, and Ulrich Heinz, “Collision geometry and flow in uranium \+ uranium collisions,” Phys. Rev. C 92, 044903 (2015), arXiv:1507.03910 [nucl-th] .
  * Giacalone  _et al._ (2018) Giuliano Giacalone, Jacquelyn Noronha-Hostler, Matthew Luzum, and Jean-Yves Ollitrault, “Hydrodynamic predictions for 5.44 TeV Xe+Xe collisions,” Phys. Rev. C 97, 034904 (2018), arXiv:1711.08499 [nucl-th] .
  * Giacalone (2020a) Giuliano Giacalone, “Observing the deformation of nuclei with relativistic nuclear collisions,” Phys. Rev. Lett. 124, 202301 (2020a), arXiv:1910.04673 [nucl-th] .
  * Giacalone (2020b) Giuliano Giacalone, “Constraining the quadrupole deformation of atomic nuclei with relativistic nuclear collisions,” Phys. Rev. C 102, 024901 (2020b), arXiv:2004.14463 [nucl-th] .
  * Giacalone  _et al._ (2021a) Giuliano Giacalone, Jiangyong Jia, and Vittorio Somà, “Accessing the shape of atomic nuclei with relativistic collisions of isobars,” Phys. Rev. C 104, L041903 (2021a), arXiv:2102.08158 [nucl-th] .
  * Giacalone  _et al._ (2021b) Giuliano Giacalone, Jiangyong Jia, and Chunjian Zhang, “Impact of Nuclear Deformation on Relativistic Heavy-Ion Collisions: Assessing Consistency in Nuclear Physics across Energy Scales,” Phys. Rev. Lett. 127, 242301 (2021b), arXiv:2105.01638 [nucl-th] .
  * Jia  _et al._ (2022a) Jiangyong Jia, Shengli Huang, and Chunjian Zhang, “Probing nuclear quadrupole deformation from correlation of elliptic flow and transverse momentum in heavy ion collisions,” Phys. Rev. C 105, 014906 (2022a), arXiv:2105.05713 [nucl-th] .
  * Jia (2022a) Jiangyong Jia, “Shape of atomic nuclei in heavy ion collisions,” Phys. Rev. C 105, 014905 (2022a), arXiv:2106.08768 [nucl-th] .
  * Bally  _et al._ (2022) Benjamin Bally, Michael Bender, Giuliano Giacalone, and Vittorio Somà, “Evidence of the triaxial structure of 129Xe at the Large Hadron Collider,” Phys. Rev. Lett. 128, 082301 (2022), arXiv:2108.09578 [nucl-th] .
  * Jia (2022b) Jiangyong Jia, “Probing triaxial deformation of atomic nuclei in high-energy heavy ion collisions,” Phys. Rev. C 105, 044905 (2022b), arXiv:2109.00604 [nucl-th] .
  * Gale  _et al._ (2013) Charles Gale, Sangyong Jeon, and Bjoern Schenke, “Hydrodynamic Modeling of Heavy-Ion Collisions,” Int. J. Mod. Phys. A28, 1340011 (2013), arXiv:1301.5893 [nucl-th] .
  * Heinz and Snellings (2013) Ulrich Heinz and Raimond Snellings, “Collective flow and viscosity in relativistic heavy-ion collisions,” Ann. Rev. Nucl. Part. Sci. 63, 123–151 (2013), arXiv:1301.2826 [nucl-th] .
  * Romatschke and Romatschke (2019) Paul Romatschke and Ulrike Romatschke,  _Relativistic Fluid Dynamics In and Out of Equilibrium_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2019) arXiv:1712.05815 [nucl-th] .
  * Abdallah  _et al._ (2022) Mohamed Abdallah _et al._ (STAR), “Search for the chiral magnetic effect with isobar collisions at $\sqrt{s_{NN}}$=200 GeV by the STAR Collaboration at the BNL Relativistic Heavy Ion Collider,” Phys. Rev. C 105, 014901 (2022), arXiv:2109.00131 [nucl-ex] .
  * Zhang and Jia (2022) Chunjian Zhang and Jiangyong Jia, “Evidence of Quadrupole and Octupole Deformations in Zr96+Zr96 and Ru96+Ru96 Collisions at Ultrarelativistic Energies,” Phys. Rev. Lett. 128, 022301 (2022), arXiv:2109.01631 [nucl-th] .
  * Xu  _et al._ (2021) Hao-jie Xu, Wenbin Zhao, Hanlin Li, Ying Zhou, Lie-Wen Chen, and Fuqiang Wang, “Probing nuclear structure with mean transverse momentum in relativistic isobar collisions,” (2021), arXiv:2111.14812 [nucl-th] .
  * Liu  _et al._ (2022) Lu-Meng Liu, Chun-Jian Zhang, Jia Zhou, Jun Xu, Jiangyong Jia, and Guang-Xiong Peng, “Probing neutron-skin thickness with free spectator neutrons in ultracentral high-energy isobaric collisions,” (2022), arXiv:2203.09924 [nucl-th] .
  * Zhao  _et al._ (2022) Shujun Zhao, Hao-jie Xu, Yu-Xin Liu, and Huichao Song, “Probing the nuclear deformation with three-particle asymmetric cumulant in RHIC isobar runs,” (2022), arXiv:2204.02387 [nucl-th] .
  * Jia  _et al._ (2022b) Jiangyong Jia, Giuliano Giacalone, and Chunjian Zhang, “Precision tests of the nonlinear mode coupling of anisotropic flow via high-energy collisions of isobars,” (2022b), arXiv:2206.07184 [nucl-th] .
  * Jia and Zhang (2021) Jiangyong Jia and Chun-Jian Zhang, “Scaling approach to nuclear structure in high-energy heavy-ion collisions,” (2021), arXiv:2111.15559 [nucl-th] .
  * Voloshin  _et al._ (2008) Sergei A. Voloshin, Arthur M. Poskanzer, Aihong Tang, and Gang Wang, “Elliptic flow in the Gaussian model of eccentricity fluctuations,” Phys. Lett. B 659, 537–541 (2008), arXiv:0708.0800 [nucl-th] .
  * Note (1) Note that a change in nuclear structure also impacts the distribution $p(N_{\mathrm{part}})$, such that the events with the same $N_{\mathrm{part}}$ correspond to slightly different centrality and vice versa Jia  _et al._ (2022c). This secondary effect leads to a small correlation between $\varepsilon_{2}^{\mathrm{rp}}$ and the amount of deformation, but this is a subleading effect compared to that discussed in Fig. 1.
  * Lin  _et al._ (2005) Zi-Wei Lin, Che Ming Ko, Bao-An Li, Bin Zhang, and Subrata Pal, “A Multi-phase transport model for relativistic heavy ion collisions,” Phys. Rev. C72, 064901 (2005), arXiv:nucl-th/0411110 [nucl-th] .
  * Ma and Bzdak (2014) Guo-Liang Ma and Adam Bzdak, “Long-range azimuthal correlations in proton–proton and proton–nucleus collisions from the incoherent scattering of partons,” Phys. Lett. B739, 209–213 (2014), arXiv:1404.4129 [hep-ph] .
  * Bzdak and Ma (2014) Adam Bzdak and Guo-Liang Ma, “Elliptic and triangular flow in $p$+Pb and peripheral Pb+Pb collisions from parton scatterings,” Phys. Rev. Lett. 113, 252301 (2014), arXiv:1406.2804 [hep-ph] .
  * Bilandzic  _et al._ (2011) Ante Bilandzic, Raimond Snellings, and Sergei Voloshin, “Flow analysis with cumulants: Direct calculations,” Phys. Rev. C 83, 044913 (2011), arXiv:1010.0233 [nucl-ex] .
  * Bilandzic  _et al._ (2014) Ante Bilandzic, Christian Holm Christensen, Kristjan Gulbrandsen, Alexander Hansen, and You Zhou, “Generic framework for anisotropic flow analyses with multiparticle azimuthal correlations,” Phys. Rev. C 89, 064904 (2014), arXiv:1312.3572 [nucl-ex] .
  * Jia  _et al._ (2017) Jiangyong Jia, Mingliang Zhou, and Adam Trzupek, “Revealing long-range multiparticle collectivity in small collision systems via subevent cumulants,” Phys. Rev. C 96, 034906 (2017), arXiv:1701.03830 [nucl-th] .
  * Zhou and Jia (2018) Mingliang Zhou and Jiangyong Jia, “Centrality fluctuations in heavy-ion collisions,” Phys. Rev. C 98, 044903 (2018), arXiv:1803.01812 [nucl-th] .
  * Jia  _et al._ (2020) Jiangyong Jia, Chunjian Zhang, and Jun Xu, “Centrality fluctuations and decorrelations in heavy-ion collisions in a Glauber model,” Phys. Rev. Res. 2, 023319 (2020), arXiv:2001.08602 [nucl-th] .
  * Note (2) The $N_{\mathrm{ch}}$ is defined as number of charged particles with $p_{\mathrm{T}}>0.1$ GeV and $|\eta|<0.5$ similar to STAR analysis Abdallah  _et al._ (2022).
  * Giacalone  _et al._ (2017a) Giuliano Giacalone, Jacquelyn Noronha-Hostler, and Jean-Yves Ollitrault, “Relative flow fluctuations as a probe of initial state fluctuations,” Phys. Rev. C 95, 054910 (2017a), arXiv:1702.01730 [nucl-th] .
  * Nijs and van der Schee (2021) Govert Nijs and Wilke van der Schee, “Inferring nuclear structure from heavy isobar collisions using Trajectum,” (2021), arXiv:2112.13771 [nucl-th] .
  * Moreland  _et al._ (2015) J. Scott Moreland, Jonah E. Bernhard, and Steffen A. Bass, “Alternative ansatz to wounded nucleon and binary collision scaling in high-energy nuclear collisions,” Phys. Rev. C 92, 011901 (2015), arXiv:1412.4708 [nucl-th] .
  * Teaney and Yan (2011) Derek Teaney and Li Yan, “Triangularity and Dipole Asymmetry in Heavy Ion Collisions,” Phys. Rev. C 83, 064904 (2011), arXiv:1010.1876 [nucl-th] .
  * Jia  _et al._ (2022c) Jiangyong Jia, Gang Wang, and Chunjian Zhang, “Impact of event activity variable on the ratio of observables in isobar collisions,” (2022c), arXiv:2203.12654 [nucl-th] .
  * ALICE Collaboration (2022) ALICE Collaboration, “Elliptic flow of charged particles at midrapidity relative to the spectator plane in Pb-Pb and Xe-Xe collisions,” (2022), arXiv:2204.10240 [nucl-ex] .
  * Giacalone  _et al._ (2017b) Giuliano Giacalone, Li Yan, Jacquelyn Noronha-Hostler, and Jean-Yves Ollitrault, “Skewness of elliptic flow fluctuations,” Phys. Rev. C 95, 014913 (2017b), arXiv:1608.01823 [nucl-th] .
  * Bhalerao  _et al._ (2019) Rajeev S. Bhalerao, Giuliano Giacalone, and Jean-Yves Ollitrault, “Kurtosis of elliptic flow fluctuations,” Phys. Rev. C 99, 014907 (2019), arXiv:1811.00837 [nucl-th] .
  * Giacalone (2019) Giuliano Giacalone, “Elliptic flow fluctuations in central collisions of spherical and deformed nuclei,” Phys. Rev. C 99, 024910 (2019), arXiv:1811.03959 [nucl-th] .
  * Miller  _et al._ (2007) Michael L. Miller, Klaus Reygers, Stephen J. Sanders, and Peter Steinberg, “Glauber modeling in high energy nuclear collisions,” Ann. Rev. Nucl. Part. Sci. 57, 205–243 (2007), arXiv:nucl-ex/0701025 .
  * Aad  _et al._ (2013) Georges Aad _et al._ (ATLAS), “Measurement of the distributions of event-by-event flow harmonics in lead-lead collisions at = 2.76 TeV with the ATLAS detector at the LHC,” JHEP 11, 183 (2013), arXiv:1305.2942 [hep-ex] .
  * Sirunyan  _et al._ (2019) Albert M Sirunyan _et al._ (CMS), “Non-Gaussian elliptic-flow fluctuations in PbPb collisions at $\sqrt{\smash[b]{s_{{}_{\text{NN}}}}}=5.02$ TeV,” Phys. Lett. B 789, 643–665 (2019), arXiv:1711.05594 [nucl-ex] .
  * Acharya  _et al._ (2018) S. Acharya _et al._ (ALICE), “Energy dependence and fluctuations of anisotropic flow in Pb-Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ and 2.76 TeV,” JHEP 07, 103 (2018), arXiv:1804.02944 [nucl-ex] .
  * Aaboud  _et al._ (2020) Morad Aaboud _et al._ (ATLAS), “Fluctuations of anisotropic flow in Pb+Pb collisions at $\sqrt{{\mathrm{s}}_{\mathrm{NN}}}$ = 5.02 TeV with the ATLAS detector,” JHEP 01, 051 (2020), arXiv:1904.04808 [nucl-ex] .