# A novel approach for calculating galaxy rotation curves using spaxel cross-
correlation and iterative smoothing

Satadru Bag1, Arman Shafieloo1,2, Rory Smith3, Haeun Chung4, Eric V.
Linder5,6,7, Changbom Park8, Y. Sultan Abylkairov7, Khalykbek Yelshibekov7,9

1 Korea Astronomy and Space Science Institute, Daejeon 34055, Korea
2 University of Science and Technology, Daejeon 34113, Korea
3 Departamento de Física, Universidad Técnica Federico Santa María, Avenida
Vicuña Mackenna 3939, San Joaquín, Santiago, Chile
4 University of Arizona, Steward Observatory, 933 N Cherry Ave, Tucson, AZ
85721, USA
5 Berkeley Center for Cosmological Physics, University of California,
Berkeley, CA 94720, USA
6 Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA
7 Energetic Cosmos Laboratory, Nazarbayev University, Nur-Sultan 010000,
Kazakhstan
8 School of Physics, Korea Institute for Advanced Study, Seoul 02455, Korea
9 Department of Physics, University of California San Diego, La Jolla,
California 92093, USA
satadru<EMAIL_ADDRESS>
###### Abstract

Precise measurements of the internal dynamics of galaxies have proven of great
importance for understanding the internal dark matter distribution of
galaxies. We present a novel method for measuring the line-of-sight (LOS)
velocities across the face of galaxies by cross-correlation of spectral pixels
(spaxels) and an iterative method of smoothing. On simulated data the method
can accurately recover the input LOS velocities for different types of spectra
(absorption line dominated, emission line dominated, and differing shapes of
the continuum), and can handle stellar population radial gradients. Most
important of all, it continues to provide reliable measurements of LOS
velocities with reasonable uncertainties even when the spectra are very low
signal-to-noise (approaching $\sim 1$), which is a challenge for traditional
template-fitting approaches. We apply our method to data from a real MaNGA
galaxy as a demonstration and find promising results with good precision. This
novel approach can be complementary to existing methods primarily based on
template fitting.

###### keywords:

Galaxy: kinematics and dynamics – methods: statistical – methods: data
analysis

††pubyear: 2021††pagerange: A novel approach for calculating galaxy rotation
curves using spaxel cross-correlation and iterative smoothing–LABEL:lastpage

##  1 Introduction

The importance of measuring the line-of-sight velocities in astronomy cannot
be overstated since they shed light on several key questions in cosmology and
astrophysics pertaining to structure formation, galaxy dynamics in clusters
(that first revealed the existence of dark matter halos (Zwicky, 1933)),
distance measures, as well as internal galaxy (and globular cluster) dynamics.
In particular, the internal dynamics of galaxies has drawn major amount of
interests for decades as this is a fundamental quantity, a part of the
empirical scaling properties that almost all galaxies follow, such as the
Tully-Fisher (Luminosity vs rotational velocity) (Tully & Fisher, 1977),
Faber-Jackson (Luminosity vs velocity dispersion) (Faber & Jackson, 1976),
Fundamental plane (Luminosity vs $R_{\rm eff}$ vs velocity dispersion)
(Gudehus, 1973; Cole et al., 1994). Furthermore, if baryons are in virial
equilibrium with the galaxy’s total potential, the dynamics can be used to
measure the total mass of the system (the dynamical mass) and the mass-
distribution (Rubin et al., 1978; Bosma, 1981; van Albada et al., 1985; Thomas
et al., 2011; Loubser et al., 2020; Kretschmer et al., 2021); see Sofue &
Rubin (2001a) for a review. This is what historically first revealed that dark
matter dominates in the outskirts of galaxies (Rubin & Ford, 1970; Roberts &
Rots, 1973; Rubin et al., 1980). One can also measure the dynamics of baryons
that are not in equilibrium, e.g., gas outflows (Riera et al., 1988; Storchi-
Bergmann, 2010; Cazzoli et al., 2016; Couto et al., 2021). Therefore, the
precise measurement of the internal dynamics of galaxies remains as one of the
most important aspects of studying galaxy evolution.

Line-of-sight velocities of baryons are typically measured by Doppler shift of
lines in the spectra. For atomic gas the primary line of interest was the 21cm
line emission (Ewen & Purcell, 1951; Catinella et al., 2006; de Blok et al.,
2008; Das et al., 2020). For stellar components, and ionised gas, one can
consider their spectral lines (H$\alpha$ and NII are especially bright lines
in the optical part of the spectrum). Traditionally this was done with long
slit spectroscopy, resulting in a 1D rotation curve (Amram et al., 1992; Vogt
et al., 1993; Sofue & Rubin, 2001b; Márquez et al., 2004; Torres-Flores et
al., 2011; Erroz-Ferrer et al., 2015). Spectra also provides valuable
information on other properties of stars (other stellar dynamics (e.g.,
dispersion), stellar populations, ages, metallicities, alpha abundances).
Recently, IFUs have become very important. Here spectra are measured in a 2D
grid of spaxels across the face of the galaxy. Many large area surveys, e.g.
CALIFA (https://califa.caha.es/) (Sánchez et al., 2012, 2016b), MANGA
(https://www.sdss.org/surveys/manga/) (Bundy et al., 2015; Drory et al., 2015;
Yan et al., 2016; Wake et al., 2017), SAMI (http://sami-survey.org/) (Scott et
al., 2018; Croom et al., 2021), and in the future Hector (Bryant et al.,
2016), have databases of thousands of galaxies.

Traditionally matching the spectra is done with template fitting (Cappellari &
Emsellem, 2004; Cid Fernandes et al., 2005; Ocvirk et al., 2006; Walcher et
al., 2006; Koleva et al., 2009; Sánchez et al., 2016a; Cappellari, 2017). We
get a lot of information from the fit - dynamics, stellar populations, etc
(Acquaviva et al., 2012; Cappellari, 2017; Boquien et al., 2019; Johnson et
al., 2021). But results are strongly influenced by assumptions built into the
templates. Also, when S/N is low, even when it is still possible to measure
stellar populations, it becomes difficult to measure line-of-sight velocities
with the template fitting approach. S/N can become low in low surface
brightness galaxies, or in the outskirts of any galaxy. To overcome this, IFU
studies attempt to sum together the spectra from multiple spaxels to try to
increase S/N using Voronoi tessellation fields (Aurenhammer, 1991; Cappellari
& Copin, 2003; García-Lorenzo et al., 2015; Ferreras et al., 2019; Roberts-
Borsani et al., 2020; Ge et al., 2021). However, there is still a limit to how
far out in radius they can go or how faint galaxies can be. For example, most
galaxies in SAMI only have detected stellar dynamics out to one or two
effective radii (Croom et al., 2021). If the aim is to study dark matter
content, this is limiting as it is only beyond $\sim 1R_{\rm eff}$ that the
dark matter begins to dominate the potential (Cappellari et al., 2013). This
in turn curtails the possibility for new discoveries from the existing data.

In this article, we seek to explore a template independent method, proposing a
novel approach to calculate the velocity differences between the pairs of
spaxels. Specifically, we estimate the Doppler shifts using the cross-
correlation of the respective spectra, with robust control over both
regularization (smoothing) and selection of the particular parts of the
spectra. Our method is not dependent on identifying individual features
(emission/absorption lines) in the spectra, and of course independent of any
template. Therefore, apart from providing a useful crosscheck on traditional
methods, it can handle difficult observational conditions (e.g. variation in
spectra types, low S/N, etc), as well as realistic complex non-ideal data
(incomplete spectra, missing spaxels etc). We demonstrate that this approach
allows success even for lower signal to noise spectra than usual, and more
robustness as it takes into account the whole spectra, not just predominantly
a few features.

The paper is organized as follows. In the next section, we explain the novel
technique based on iterative smoothing and computing the cross-correlation
between a pair of spectra. We validate the approach on a variety of
simulations in section 3 after demonstrating the method on a pair of spaxels.
In the following section, we push the method to limit by systematically
testing on data with different noise levels from high to very low S/N. In
section 5, we validate the method on a MaNGA galaxy and compare the results
with that from Marvin. We conclude and put useful discussion in section 6.

##  2 Methodology

Figure 1: The galactic plane is shown by the shaded surface. Suppose the line
of sight (LOS) has an angle $\alpha$ with the rotational axis. The velocity
components along the projection of LOS for a point on the galactic plane would
be $v(r)\sin\beta$. This can be further projected along the LOS by multiplying
with the $\sin\alpha$ term.

For a galaxy observed with a spectrograph, we obtain spaxel (spectral pixel)
for different parts of the galaxy, whether a 2D array in the case of an
integral field spectrograph or a line from slit spectroscopy. By comparing one
spaxel to another, one can measure wavelength shifts of spectral features in
common.

The wavelength shift from spaxel 1 to 2 is directly related to the line of
sight velocity difference through the Doppler shift, so

$\Delta V=c\ \frac{\lambda_{2}-\lambda_{1}}{\lambda_{1}}\;.$ (1)

In order to relate this to the galaxy rotation curve (GRC) we must take into
account the inclination of the galaxy with respect to the line of sight, with

$\displaystyle V_{1}$ $\displaystyle=V_{\rm
sys}+v(r_{1})\sin\alpha\,\sin\beta_{1}\;,$ (2) $\displaystyle V_{2}$
$\displaystyle=V_{\rm sys}+v(r_{2})\sin\alpha\,\sin\beta_{2}\;,$

where the angles are as illustrated in Figure 1, $V_{\rm sys}$ is the line of
sight velocity of the central spaxel, or in the simplest case of the whole
galaxy, and $v(r)$ is the rotational velocity of a spaxel at distance $r$ from
the galactic centre.

The velocity difference between spaxels as measured by the wavelength shift
therefore allows one to build up the galaxy rotation curve $v(r)$ through

$\Delta
V=\sin\alpha\,\left[(v(r_{2})\sin\beta_{2}-v(r_{1})\sin\beta_{1}\right]\;,$
(3)

given the values of the angles $\alpha$ and $\beta$. One repeats the process
for many spaxel pairs, and hence $\Delta V$’s, to construct the GRC.

In this paper we focus on introduction of our technique for measuring $\Delta
V$; the extraction of $v(r)$ from those is well established in the literature.
Our basic approach is to cross-correlate the spectra from spaxel pairs. Since
observed spectra have noise, one needs to avoid spurious features from the
noise overwhelming the true wavelength shift. In addition, spectra are
sometimes incomplete, missing wavelength regions due to, for example,
intervening atmospheric lines.

Therefore we first smooth one spectrum of the pair while leaving that of the
other intact; we will then also switch which of the pair is smoothed to
compare results and prevent a single high noise spectrum from distorting the
estimate. We cross-correlate the spectra pair as a function of offset $\delta
V_{\rm off}$ for a wide range of offsets, building up a cross-correlation
function $r(\delta V_{\rm off})$. The mode of this curve is expected to be the
actual $\Delta V$.

In the next subsections, we go into the technical details of how we carry out
the smoothing, estimate the cross-correlation function, and then use it to
obtain $\Delta V$.

###  2.1 Iterative smoothing of the spectra

We smooth an observed spectrum, $F_{\rm obs}(\lambda_{i})$, iteratively with a
Gaussian kernel following Shafieloo et al. (2006); Shafieloo (2007); Shafieloo
& Clarkson (2010); Aghamousa & Shafieloo (2015); these papers establish the
criteria for stable, accurate results. The smoothed spectrum in the $n$th step
is obtained from the previous step as

$F^{s}_{n}(\lambda)=F^{s}_{n-1}(\lambda)+\frac{1}{N(\lambda)}\sum_{i}\frac{\left(F_{\rm
obs}(\lambda_{i})-F^{s}_{n-1}(\lambda_{i})\right)}{{\sigma_{\rm
obs}}^{2}(\lambda_{i})}\times\exp{\left[-\frac{(\lambda-\lambda_{i})^{2}}{2\Delta^{2}}\right]}$
(4)

where the normalisation term $N(\lambda)$ is given by

$N(\lambda)=\sum_{i}\left(\frac{1}{{\sigma_{\rm
obs}}^{2}(\lambda_{i})}\right)\times\exp{\left[-\frac{(\lambda-\lambda_{i})^{2}}{2\Delta^{2}}\right]}\;.$
(5)

Here $F^{s}_{n-1}(\lambda)$ is the smoothed spectrum obtained in the previous
step, i.e. at the $n-1$ step. We start with an initial guess for the smoothed
spectrum, say $F^{s}_{0}(\lambda)$= constant in the first step, and continue
smoothing for $N_{\rm it}$ number of iterations. For sufficiently number of
iteration, the final smoothed spectra should become independent of the initial
guess. Note that we have two parameters in this smoothing algorithm: the
smoothing scale $\Delta$ and the number of iterations $N_{\rm it}$. These two
parameters should be fixed according to the problem at hand. For example,
$\Delta\gtrsim\delta\lambda$ where $\delta\lambda$ is the average wavelength
resolution in the observed spectrum and $N_{\rm it}\gtrsim 5$ work well for
the typical examples we consider in this work. Testing by simulations as
described in Section 3 establishes that the final results are properly
insensitive to the choice of $\Delta$ and $N_{\rm it}$ when they are within
reasonable ranges. Also note that in this method we can obtain the smoothed
spectrum $F^{s}$ at any desired values of $\lambda$ in between the observed
wavelengths $\\{\lambda_{i}\\}$. Thus this smoothing algorithm also serves the
purpose of interpolation needed for cross-correlations with any arbitrary
Doppler (wavelength) shift.

###  2.2 Estimating $\Delta V$ using cross-correlation

For different relative velocities ($\Delta V$) we cross-correlate the spectra
considering the corresponding Doppler shift in the wavelength (related through
equation (1)), with one of them being smoothed using (4). (We also tried
initial “supersmoothing”, using a very large $\Delta_{\rm ini}$ to remove the
continuum trend; this did not give any advantage over the results we present.)
We define the weighted cross-correlation between the spectra of two spaxels
(say A and B) as a function of $\Delta V$ as

$r_{A^{s}B}(\Delta V)\equiv F^{s}_{A}(\lambda+\Delta\lambda)\otimes
F_{B}(\lambda)=\frac{\sum_{i}w_{i}\left[F^{s}_{A}(\lambda_{i}+\Delta\lambda)-\left<F^{s}_{A}(\lambda_{i}+\Delta\lambda)\right>_{w}\right]\left[F_{B}(\lambda_{i})-\left<F_{B}(\lambda_{i})\right>_{w}\right]}{\sqrt{\sum_{i}w_{i}\left[F^{s}_{A}(\lambda_{i}+\Delta\lambda)-\left<F^{s}_{A}(\lambda_{i}+\Delta\lambda)\right>_{w}\right]^{2}}~{}\sqrt{\sum_{i}w_{i}\left[F_{B}(\lambda_{i})-\left<F_{B}(\lambda_{i}\right>_{w}\right]^{2}}}\;,$
(6)

where $\Delta\lambda=\lambda(\Delta V/c)$. The index $i$ runs over the
wavelengths sampled in the spectra. We employ inverse variance weighting
$w_{i}=1/{\sigma^{2}_{B}}_{i}$, where ${\sigma_{B}}_{i}$ is the uncertainty of
the second, unsmoothed spectrum. The correlation coefficient $r_{AB}(\Delta
V)$ should have a maximum when $\Delta V$ is the same as the true velocity
difference between the two spaxels.

###  2.3 The algorithm to estimate $\Delta V$

Given these ingredients, the algorithm for obtaining the velocity difference
between a pair of spaxels A and B is:

  * •

Smooth the spectra observed in one of the spaxels, say A, calling the smoothed
spectrum $F_{A}^{s}$.

  * •

Cross-correlate $F_{A}^{s}$ and $F_{B}$ for a large number of relative
velocities ($\Delta V$). For added robustness through crosschecks, instead of
cross-correlating the whole spectra to give a single answer, we divide the
(unsmoothed) spectrum of B into $N_{B}$ number of bins to utilize information
from different parts of the spectrum independently. We determine the cross-
correlation (6) for each bin of $F_{B}$ separately, i.e. we calculate separate
estimates $r_{A^{s}B_{j}}(\Delta V)$ for the $j$th bin where
$j\in\\{1,N_{B}\\}$.

  * •

Find the global maxima of the correlation coefficients $r_{A^{s}B_{j}}(\Delta
V)$ for each bin. Let the maximum for the $j$th bin be $r^{\rm
max}_{A^{s}B_{j}}$, occurring at the velocity difference $\Delta V_{AB_{j}}$.
We will have $N_{B}$ such velocity estimations corresponding to $N_{B}$ bins.

  * •

Swap the two spaxels and repeat the same procedure, i.e. evaluate
${r_{B^{s}A}}_{j}(\Delta V)$ and find $\Delta V_{BA_{j}}$ (always using the
variance of the unsmoothed spectrum as the weight factor), again separately
for each of the $N_{B}$ bins. We now have $2$$N_{B}$ number of velocity
estimations. We expect $\Delta V_{BA}\approx-\Delta V_{AB}$ for all the bins;
this is another crosscheck.

  * •

Estimate the velocity difference from these $2$$N_{B}$ estimations. As
different parts of the spectra may have different signal to noise ratios,
different bins can give rise to $\Delta V$ estimation with different degrees
of accuracy. To filter out noise we only accept estimates satisfying the
symmetry requirement: the velocity estimations $|\Delta V_{AB}|$ and $|\Delta
V_{BA}|$ are close to each other, with $|\Delta V_{AB_{j}}+\Delta
V_{BA_{j}}|\leq 0.05\,|\Delta V_{AB_{j}}-\Delta V_{BA_{j}}|$, or $\leq 10$
km/s111 Note that this criterion could be set according to the problem in
hand, considering several factors like the signal to noise ratio of the
spectra..

  * •

From the acceptable estimates determine the final velocity difference $\Delta
V_{AB}$ between the spaxels by averaging over these estimates. Estimate the
uncertainty $\sigma_{\Delta V_{AB}}$ using the standard deviations.

Cross-correlating one unsmoothed spectrum with another smoothed spectrum helps
to increase the effective signal to noise. Without smoothing either spectrum,
the correlation coefficient curve $r(\Delta V)$ exhibits wiggles because of
random matching of the two spectral noises. This may alter the locations of
the maxima of $r(\Delta V)$ for different bins leading to erroneous
estimations, especially for low signal to noise cases. Conversely, smoothing
both spectra simultaneously leads to neglect of some uncertainties in the
data, and hence affects the realistic estimation of the final $\Delta V$. By
smoothing each spectrum in sequence in addition we can have a first, basic
consistency test between $\Delta V_{AB}$ and $\Delta V_{BA}$ to make sure that
we are not fitting noise and both estimations are meaningfully consistent.

##  3 Validating the method using simulation

Our aim is to test the new technique’s ability to measure the shape of a model
galaxy’s rotation curve in diverse observational circumstances. Galaxies have
different spectra, e.g. red and dead galaxies have absorption lines, while
star-forming galaxies may be dominated by emission lines. And galaxies often
have differing spectra as a function of radius (e.g. some spiral galaxies have
red and dead centres and star formation in their disks). We need to ensure our
method is able to successfully recover the input rotation curve even when
different types of spectra are considered. Therefore, we test the method
through simulations of spectra of galaxies in different scenarios, e.g. with
absorption vs emission line dominated spectra, with different noise levels and
strengths of features, radially varying spectral properties etc.

In this first work, for simplicity we always simulate the model galaxy to be
edge-on and the spaxels are placed on the major axis perpendicular to the line
of sight (so $\alpha=\beta=\pi/2$). That is, we demonstrate the procedure in a
1D situation, with the full 2D galaxy data treated in a follow-up paper; this
only affects determination of $v(r)$, not $\Delta V$. The velocity difference
between two spaxels then simply becomes the difference in their rotational
velocities, i.e. $\Delta V_{AB}=v(r_{A})-v(r_{B})$. We assume an exponential
disk for the model galaxy’s stellar disk with an effective radius of $1.5$
kpc. The spaxels are distributed from $-5$ to 5 kpc, and there are 21 spaxels
in total such that spaxel 11 is centred on the disk centre. For the line-of-
sight dynamics of the stars in each spaxel, we assume a Polyex model
(Giovanelli & Haynes, 2002; Catinella et al., 2006):

$v_{\rm PE}(r)=v_{0}\left(1-e^{-r/r_{\rm PE}}\right)\left(1+\frac{\gamma
r}{r_{\rm PE}}\right)$ (7)

where $v_{0}$, $r_{\rm PE}$ and $\gamma$ are free parameters that can be
varied to alter the shape of the stellar dynamics at each radius. When we
carry out validation fits, we take two cases: one where the input dynamics is
known (so are these parameters) for testing, in the other it is blinded. We
will find that this does not affect the fit quality.

We now build the stellar spectrum of each spaxel. In order to test with
realistic spectra, we generate the spectra in the simulations using the
observed spectra of ‘benchmark’ spaxels from two rather different MaNGA
galaxies, 7991–12701 and 8952–9102222We obtain the IFU data of MaNGA galaxies
from Marvin (Cherinka et al., 2019) DR15 given at
https://dr15.sdss.org/marvin/.. The spectra in all spaxels for the former
galaxy are dominated by strong emission lines; spectra in two spaxels
including the central one are shown in figure 2. On the other hand, the
spectra of the MaNGA galaxy 8952–9102 show different properties in different
spaxels: the central ones are absorption line dominated while the edge ones
have both weak absorption and weak emission lines. Spectra from some of these
spaxels are shown in figure 3. We then apply a Doppler shift to the spectra,
according to the dynamics that the Polyex model provided. We also inject white
noise into the MANGA galaxy spectra so that the signal-to-noise in each spaxel
will reduce as the surface-brightness of the galaxy’s disk becomes fainter
exponentially in the disk outskirts, as occurs in real galaxies.

We carry out the validation procedure in steps, from initially assuming
spectral homogeneity and noise homogeneity, then relaxing these conditions one
by one. Our model should be considered a toy model for mocking up real IFU
observations of galaxies. It allows us to test the method’s ability to recover
the input dynamics in a controlled manner, with a realistic set of galaxy
spectra, a reasonable choice of rotation curve shapes, and with differing but
controlled amounts of noise injected.

We measure the signal to noise ratio (S/N hereafter) of a whole spectrum
following Stoehr et al. (2008). Figure 4 shows how S/N of the spectra vary
across the spaxels for the different tests that we pursue in this section. The
details of the tests involving a variety of different observational conditions
are explained below in the respective subsections.

Figure 2: The MaNGA galaxy 7991–12701 has emission line dominated spectra,
shown here from the two spaxels (37,37) and (37,12), given by their
coordinates in the 2D IFU grid. The flux observed in MaNGA is in the standard
unit of [$10^{-17}$ erg/cm2/s/$\AA$] which has been used as the unit of flux
in the relevant figures throughout the article. The two spectra exhibit a
similar set of features on different continua, and have a slight wavelength
shift with respect to each other that is the manifestation of their velocity
difference ($\Delta V\sim 168$ km/s, so the 0.06% shift is difficult to see by
eye). The spectra in the central spaxel (37,37) has been used in the
simulation for demonstrating the method in section 3.1.1 and for the test
presented in section 3.1. The test demonstrated in section 3.3 uses both these
spectra in the simulation.  Figure 3: The MaNGA galaxy 8952–9102 has
absorption line or weak emission line dominated spectra, shown here from the
three spaxels (32,32), (32,37) and (32,43). The spectra for the central spaxel
(32,32) has been used for the test demonstrated in section 3.2. The spectrum
from the spaxels (32,37) and (32,43) are used for the test in section 3.4.
Figure 4: Signal to noise ratio (S/N) for fluxes in the spaxels for different
tests illustrated in sections 3.1–3.4. The solid and dotted curves represent
two different realizations for each test. Note in the mixed spectra cases
(green and red curves) we challenge the method with step function transition
in spectra type; we further challenge it in the lowest S/N case (red curve) by
flipping the sign of the transition between the two realizations.

###  3.1 Simulation with emission line spectrum

We begin validation by simulating spectra in 21 spaxels located on the major
axis of a galaxy, but fixing the properties of all the spectra in different
spaxels to that of the spectrum of the central spaxel of MaNGA 7991-12701
(shown in orange in figure 2). Each is shifted according to the line of sight
velocities input in the simulation. The spectra have been observed typically
within the wavelength range $\lambda=[3621.6,10339.5]$ Å with an average
interval of $\delta\lambda\approx 1.47$ Å. However, we remove all the data
points above $\lambda>9850$ Å due to the telluric contamination arising from
the earth’s atmosphere. For this first stage we also keep the noise level low
and the same for all the spaxels. Note that even with a uniform noise level,
the spectra simulated in the outer/edge spaxels would have worse signal to
noise due to the outer disk of a galaxy being fainter.

####  3.1.1 Determining the velocity differences

Since we are working with a line of simulated spaxels along the major axis of
a galaxy, we can refer to them simply by their second coordinate; this also
serves to distinguish simulated spaxels (with a single number) from MaNGA’s
observed spaxels (denoted with the 2D IFU coordinates). Consider simulated
spectra in a pair of spaxels, $A=7$ and $B=15$, shown in figure 5. They are
dominated by strong emission lines and derived from the same underlying
spectra that is observed in the central spaxel of the MaNGA 7991-12701 galaxy.
The true velocity difference between these spaxels is $326.83$ km/s (the
corresponding Doppler shift is evident from the inset which zooms into the
spectral region of H$\alpha$ and NII emission lines). For each of $N_{B}=4$
wavelength bins, we calculate the cross-correlation between one spaxel and the
other (smoothed) one as a function of applied velocity shifts.

Figure 5: Spectra simulated for the initial tests in Section 3.1, in two
spaxels, 7 and 15, based on the spectrum observed in the central spaxel of
MaNGA 7991-12701 (orange curve in figure 2). The simulated spectra are shifted
according to the input line of sight velocity difference (as shown in the
inset which zooms into the spectral range of H$\alpha$ and NII emission lines)
and sampled in equispaced wavelengths. After discarding the $\lambda>9850$ Å
region due to telluric contamination, the rest of the data has been equally
divided into 4 bins (each bin encloses the same number of datapoints), shown
by the vertical dashed lines.

Figure 6 displays the correlation coefficient as a function of the velocity
difference, $\Delta V$, between the pair of spaxels. The left and right panels
show $r_{A^{s}B}$ and $r_{B^{s}A}$ respectively for $4$ bins. For this
illustration we use our baseline $\Delta=1.5$ Å and $N_{\rm it}=10$, later
showing robustness to variation of these parameters. In both panels we find
that $r$ exhibits global maxima at very similar values of $\Delta V$ (taking
proper account of the sign), not only for $A^{s}B$ and $B^{s}A$, but also for
different wavelength bins. Importantly, although the second and the fourth bin
do not contain any strong feature, as evident from figure 5, the correlation
$r$ shows a prominent peak at the correct $\Delta V$ values; this arises from
the presence of many small features in these bins, and holds for both panels
of figure 6. This illustrates the key benefit of this method in that it does
not rely on single or few strong identifiable features in the spectra.

The maxima in the left and right panels arise at near-identical absolute
values of $\Delta V$, i.e. $\Delta V_{BA}\approx-\Delta V_{AB}$ for all the
bins as expected for an accurate shift estimation. Since both the spectra are
derived from the same spectra (the central one of MaNGA 7991-12701) and have
little noise, all the maxima in this first example have very high values for
the correlation coefficient, $r\approx 1$. Following the algorithm described
above in section 2.3, we find all 4 bins satisfy the selection criteria. Thus
we have $4\times 2=8$ ‘good’ estimations of $\Delta V$ from both the panels.

(Note that the secondary maxima in the third bin are readily identified as due
to correlation between the $H\alpha$ line (6564.6 Å restframe) and the NII
line (6585.3 Å restframe), giving a $\pm 944$ km/s offset relative to the
truth. This is an example of using not just statistics, but astrophysical
knowledge, in assessing the data.)

From these eight maximum correlations we compute the mean and standard
deviation for the velocity difference. The final result is $\Delta
V_{AB}\approx 326.94\pm 0.50$ km/s (for our baseline $\Delta=1.5$ Å, $N_{\rm
it}=10$, and 4 bins). This is wholly consistent with the simulation input
$\Delta V_{7,15}=326.83$ km/s, demonstrating validation for this first test.

In table 1, we present the estimations of $\Delta V$ for the same pair of
spaxels as in figure 6 but for various choices of $\Delta$, $N_{\rm it}$, and
$N_{B}$ (number of bins). All our estimations corresponding to these different
choices are highly accurate with very small uncertainty and consistent with
each other demonstrating the robustness of this approach. We choose as our
baseline $\Delta=1.5$ Å, $N_{\rm it}=10$, and $N_{B}=4$ for the best
combination of accuracy, precision, and computational time efficiency.

(a) $7\bigotimes 15$ (b) $15\bigotimes 7$

Figure 6: The correlation coefficient $r(\Delta V)$ plotted as a function of velocity difference $\Delta V$ between a pair of spaxels, $A=7$ and $B=15$ from the simulation based on the central spaxel of MaNGA galaxy 7991-12701. The left and right panels show $r_{A^{s}B}$ and $r_{B^{s}A}$ respectively for $4$ bins, calculated using (6). All the bins have maxima at similar $\Delta V$ in either panel and at the opposite sign value of the other panel – important crosschecks.  $\Delta$ (in Å) | Estimated velocity difference ($\Delta V$) in km/s
---|---
$4$ bins |  $8$ bins
$N_{\rm it}=10$ | $N_{\rm it}=20$ | $N_{\rm it}=10$ | $N_{\rm it}=20$
$1.5$ | $326.9\pm 0.5$ | $326.8\pm 0.8$ | $326.8\pm 0.6$ | $326.7\pm 1.0$
$2.0$ | $326.8\pm 0.6$ | $326.6\pm 1.1$ | $326.7\pm 0.7$ | $326.5\pm 1.2$
$3.0$ | $327.0\pm 0.8$ | $326.9\pm 0.7$ | $326.8\pm 0.6$ | $326.7\pm 0.6$
$4.0$ | $327.4\pm 1.7$ | $327.4\pm 1.7$ | $327.0\pm 1.5$ | $327.0\pm 1.4$
Table 1: Velocity difference between the spaxels 7 and 15 (from the simulation test in section 3.1) for different values of smoothing scale $\Delta$, number of iterations $N_{\rm it}$, and number of bins $N_{B}$ considered. The true velocity difference between these two spaxels is $326.8$ km/s. The estimations of $\Delta V$ for all the above choices of $\\{\Delta,N_{\rm it},N_{B}\\}$ match extremely well with the truth. However, we choose $\Delta=1.5$ Å, $N_{\rm it}=10$, and $N_{B}=4$ as our base line for the best combination of accuracy, precision, and computational time efficiency.  Pair of spaxels (A,B) | Our estimation of $\Delta V_{AB}$ in km/s | True $\Delta V_{AB}$ in km/s | residual in km/s
---|---|---|---
$(1,2)$ | $5.90\pm 0.25$ | $6.21$ | $-0.31\pm 0.25$
$(1,21)$ | $453.15\pm 0.69$ | $452.59$ | $0.56\pm 0.69$
$(2,13)$ | $330.27\pm 0.78$ | $329.91$ | $0.36\pm 0.78$
$(4,21)$ | $432.08\pm 0.70$ | $431.20$ | $0.88\pm 0.70$
$(5,6)$ | $12.79\pm 0.70$ | $13.18$ | $-0.39\pm 0.70$
$(7,11)$ | $163.76\pm 0.44$ | $163.48$ | $0.28\pm 0.44$
$(7,15)$ | $326.94\pm 0.50$ | $326.83$ | $0.11\pm 0.50$
$(9,14)$ | $251.14\pm 0.78$ | $251.37$ | $-0.23\pm 0.78$
$(10,21)$ | $292.45\pm 0.34$ | $291.98$ | $0.46\pm 0.34$
$(11,15)$ | $163.01\pm 0.85$ | $163.34$ | $-0.33\pm 0.85$
$(11,17)$ | $194.32\pm 0.71$ | $193.96$ | $0.36\pm 0.71$
$(13,19)$ | $104.20\pm 0.89$ | $103.51$ | $0.69\pm 0.89$
$(16,19)$ | $31.15\pm 0.31$ | $31.55$ | $-0.40\pm 0.31$
Table 2: We compare our estimation of velocity difference $\Delta V_{AB}$
between various pair of spaxels with the corresponding truth values used in
the simulation. Note that spaxel 11 is the central spaxel in the simulation,
so (1,21) spans from one edge of the galaxy to the other.

To assess the technique further we estimate the velocity difference for
several other pairs of spaxels spanning the range of 21 spaxels (with spaxel
11 the central one). In table 2, our estimated $\Delta V$ values are compared
with the corresponding true velocity differences used in the simulation. Again
we find excellent matches between our estimations and the truths, with
uncertainty below $1$ km/s in all these cases. Note that such a small
uncertainty arises due to the homogeneity of the spectral properties and (low)
noise across all the spectra. The remaining subsections of this section will
validate the method on simulations under more challenging observational
circumstances, before proceeding to real observed data in section 5.

####  3.1.2 Combining information from different spaxel pairs

Before proceeding further we briefly explain here how we construct the galaxy
rotation curves. Without loss of generality, we can assume the central spaxel
has zero rotational velocity, i.e. we compute the velocity of each spaxel on
the diameter (more precisely along the major axis of the galaxy) relative to
the central velocity and construct the galaxy rotation curve 333 By setting
the velocity of the central spaxel to zero, without loss generality, we are
essentially measuring the LoS velocities of the spaxels relative to the
reference central spaxel. In reality, the central spaxel may not be
identified. In such cases one can assume the symmetry in the rotation curve to
determine the centre and estimate the regularised spaxel velocities. .

However, the spectrum of the central spaxel itself could be noisy or
distorted, so to avoid such an issue we use the estimated velocity differences
between all the pairs, i.e. $\Delta V_{i,j}$ where $i,j$ represents all the
spaxels (at least all the spaxels with good S/N)444When spaxels have poor S/N
we still include them as $i$ spaxels but not as $j$ spaxels, i.e. we do not
measure good spaxels relative to bad. This does not occur for tests in section
3 but see section 5.2 for further details.. We then fit for the velocities of
all the spaxels simultaneously using all the $\Delta V_{i,j}$ measurements,
i.e. all spaxel pairs. Defining the velocity of the central spaxel to be zero,
we thus have $N_{\rm spax}-1$ free parameters in the fit, and up to $N_{\rm
spax}(N_{\rm spax}-1)/2$ measurements, where $N_{\rm spax}$ is the total
number of spaxels. Since $N_{\rm spax}$ can be large in principle (e.g. around
a few thousand for the MaNGA IFU observations on a 2D grid), we employ
Hamiltonian Monte Carlo (HMC) sampling with the help of the pystan package
(Carpenter et al., 2017; Stan Development Team, 2017). The fit takes care of
the consistency between the velocities of the spaxels, i.e. the triangle
inequality $\vec{\Delta}V_{AB}+\vec{\Delta}V_{BC}=\vec{\Delta}V_{AC}$, or
alerts us to inconsistency by not converging or giving a poor $\chi^{2}$.

####  3.1.3 Constructing the galaxy rotation curves

Now we complete the construction of the galaxy rotation curve for the galaxy
in the case where the spectra in all 21 spaxels are generated from the
observed spectrum in the central spaxel of MaNGA 7991-12701. All pairs of
spaxels have good signal to noise ratio (S/N $\sim 60$ as evident from figure
4) so we include the $\Delta V$ estimated for all the pairs into the fit.

To ensure the method is fairly tested, we conduct each test in pairs
considering two realizations (as described above in equation 7), each with
independent input parameters and noise realizations. The first test is
conducted such that the input rotation curve is known, and then followed up by
a second test which is similar except the input rotation curve is not revealed
until the results have been collected (a blind test). In the top left panel of
figure 7 we compare our estimated galaxy rotation curve with the true one for
the blind case555The results for the ‘known cases’ are statistically same as
the blind cases and, therefore, are not shown separately in this article.. We
find an excellent match between our estimations and the true velocities. The
uncertainties in estimations are so small ($\sim 0.1$ km/s) that they are not
visible in the $v(r)$ plot. The subplot at the top of the panel shows the
residuals of our estimations, which are typically $\mathcal{O}(0.1)$ km/s .
Since all the $\Delta V_{ij}$ estimations, which are correlated, are
considered in the fit, “chi-by-eye” is not accurate. We compute
$\chi^{2}\equiv\sum_{ij}({\bf V_{\rm est}-V_{\rm true}})_{i}\cdot{\bf
Cov}^{-1}_{ij}\cdot({\bf V_{\rm est}-V_{\rm true}})^{T}_{j}$ where ${\bf Cov}$
is the covariance matrix of the fit velocities. One expects $\chi^{2}\sim 20$
since the velocity of the central spaxel out of $21$ is fixed to zero. The
$\chi^{2}$ value is $22.29$, indicating that our method is both accurate and
statistically robust for estimating the galaxy rotation curve.

In the following sections we increase the difficulty by changing the spectrum,
noise, and adding inhomogeneity.

(a) Emission line dominated (section 3.1) (b) Absorption line dominated
(section 3.2) (c) Inhomogeneous spectral properties:
emission line dominated (section 3.3) (d) Inhomogeneous spectral properties:
weak emission vs absorption line dominated (section 3.4)

Figure 7: The four panels present the results of the tests described in
sections 3.1 – 3.4 clockwise from the top-left. The top two panels compare the
our estimated galaxy rotation curves with the corresponding truths when the
spectra in different spaxels are generated from a single base spectrum – top
left panel (a): emission line dominated (the orange spectrum in figure 2), top
right panel (b): absorption line dominated (the blue spectrum in figure 3).
The bottom two panels corresponds to the tests when the spectral properties
change like a step function along the major axis as described in sections 3.3
and 3.4. The analyses are carried out blindly and when the truths are revealed
we find that the velocities are recovered accurately for all the cases. The
uncertainty in $v(r)$ estimation is so small that the errorbars are not
visible in the plots. The subplot at the top of each panel shows the residual
for each estimation along with the errorbars. The residuals are very small,
especially for the top two cases with homogeneous spectral properties. The
$\chi^{2}$ values for all four cases are close to $20$, as expected for $20$
degrees of freedom. Note that there is (modest) correlation between the fit
velocities; the $\chi^{2}$ estimation considers the whole covariance matrix,
not just the diagonal terms.

###  3.2 Absorption line dominated spectra

Next we test how our method performs when the galaxy spectra do not contain
any sharp emission line (in other words, if the galaxy is not star-forming),
rather they are dominated by absorption lines. We again simulate spectra in
$21$ spaxels along the major axis of a galaxy. For simplicity, the spectral
properties for all spaxels here are fixed to that of the central spaxel
$(32,32)$ of the MaNGA galaxy 8592–9102, shown in figure 3 by the blue colour.
The comparison between our estimation of rotational velocities of the spaxels
with the truth is presented in the top right panel of figure 7. Again we find
an excellent match with the truth and the small errorbars are not visible in
the $v(r)$ plot. The residuals, shown in the top subplot, are again quite
small, however the residuals and their uncertainties here are slightly larger
as compared to the case with the emission line dominated spectra because the
emission lines are typically stronger than the absorption lines. E.g., the
average residual in the absorption line dominated case is $\sim 0.22$ km/s
which is slightly larger compared to that of the emission line dominated case
($\sim 0.14$ km/s). The $\chi^{2}$ value for this case (taking the full
covariance matrix of the fit velocities into account) $22.32$ for 20 degrees
of freedom, is again reasonable.

###  3.3 Different spectral properties: emission line dominated

In the previous two subsections, we tested the method on the simple scenarios
when the spectral properties are kept the same across all the spaxels (only
shifted according to the velocity differences along line of sight). But in
reality, galaxies can go from having absorption to emission lines as one moves
across them, because of stellar population gradients. So here we test the
scenarios when the spectral properties are inhomogeneous. To challenge the
technique, we take a radical step function in radius, where “inner” spaxels
have a different spectrum than “outer” spaxels. This and the next subsection
take two different models for the base spectra.

Here, both the inner and outer spaxels are simulated from strong emission line
dominated spectra, but now only the inner spaxels come from the central spaxel
($37,37$) of the MaNGA galaxy 7991–12701, while the outer spaxels are
simulated with the spaxel ($37,12$) of 7991–12701. These two emission line
dominated model spectra have been shown in figure 2. While they have the same
emission lines, their continua are quite different; the former has red slope
continuum whereas the latter has a bluer continuum. This subsection tests our
method on this case of spectral property inhomogeneity.

We compare our estimation of rotational velocities (of the spaxels) against
the true velocities for this scenario in bottom left panel of figure 7. Since
there is a sharp transition of the spectrum, we get slightly stronger
correlations in the fit velocities among the inner spaxels (8-14, the central
one is spaxel 11) and the outer ones. This leads to the pattern apparent in
the residual plot shown at the top subplot. Nevertheless, even in this extreme
step function case the $\chi^{2}$ is $21.12$ for $20$ degrees of freedom.

###  3.4 Different spectral properties: absorption line vs weak emission line
dominated spectra

Next we test the method with simulations again involving a step function in
spectral properties, now with the inner spaxels based on weak absorption line
dominated spectra and the outer spaxels based on weak emission line dominated
spectra. Specifically, we use the (32,37) spaxel of the MaNGA galaxy
8952–9102, shown as the orange spectrum (with weak absorption lines) in figure
3, and the (32,43) spaxel of the same galaxy, the green spectrum (with weak
emission lines) in that figure. Note that these two base spectra have much
weaker features compared to the emission line dominated spectra shown in
figure 2 and analysed in section 3.3.

The comparison between our estimation of the rotational velocities with the
true ones appears in the bottom right panel of 7. Our estimation matches very
well with the truth even with the step function variation in the spectral
properties among spaxels where the features are weak and of different types
(absorption vs emission). Again, due to the sharp transition in the spectra we
get some structure in the fit as evident from the pattern in the residuals
shown in the top subplot. Nevertheless, the $\chi^{2}$ value is good: $21.9$.
Thus, even in this more inhomogeneous case, the method succeeds.

##  4 Validating on simulations with different noise levels

Figure 8: Emission line dominated spectra are shown for four sets (used in
section 4.1) with different noise levels – low, medium, high, and very high –
from the top to the bottom panel. Each panel shows spectra in three spaxels,
two edge spaxels and the central one. S/N for these four sets are as low as
$59.06,~{}18.56,~{}3.97,~{}0.78$, as shown in the respective panels.

(a) Emission line dominated spectra (b) Absorption line dominated spectra

Figure 9: The variation of signal to noise ratio (S/N) of the spectra across
the spaxels is shown for the four sets with different noise levels – low,
medium, high, and very high. The left panel is for simulations from an
emission line dominated spectrum (section 4.1); the right panel for an
absorption line dominated spectrum (section 4.2). Note that S/N of the three
spaxels in each panel of figure 8 are shown by the middle (central spaxel) and
the two boundary dots (edge spaxels) on the corresponding curves in the left
panel. Same applies for figure 11 and the right panel.

In the previous section we apply the method to different scenarios while
keeping the noise level low, giving high S/N in the spectra across the spaxels
(S/N $\gtrsim 10$ in the worst cases). Because of that our estimations are not
only in extremely good agreement with the truth in every case but have quite
small uncertainty. In this section we systematically test the method on the
spectra with different noise levels.

###  4.1 Emission line dominated

We again simulate spectra in the 21 spaxels generated from the same emission
line dominated spectrum, namely of the central spaxel ($37,37$) of the MaNGA
galaxy 7991–12701 (shown by orange in figure 2), now with four different noise
levels. Figure 8 shows the spectra in three spaxels (two edge spaxels and the
central one) for the four sets with four noise levels – low (similar to in the
previous section), medium, high, and very high – from top to bottom panel.
Notice that the spectra of the outer spaxels have much larger noise, because
of being faint, most evident in the high and very high noise sets shown in the
bottom two panels. The spectra in the edge spaxels still contain the sharp
emission lines, however, the continua appear extremely noisy. The S/N in the
spectra of the spaxels are shown in the left panel of figure 9; note the rapid
drop for the medium and especially for the high and very high noise sets going
radially outwards. The spectra in the edge spaxels have S/N $\sim 4$ for the
high noise set and S/N as low as $\sim 0.8$ for the very high noise set. Note
that for all the sets we keep the true velocities the same so that the results
can be compared easily.

For the higher noise sets we find that a slightly larger smoothing scale gives
better accuracy and robust results. Thus, we use $\Delta=2.0$ Å for the medium
and higher noise sets while keep $\Delta=1.5$ Å for the low noise set, same as
before. In the three panels of figure 10 we compare the estimated rotational
velocities for each of the four noise level sets with the true velocities
(shown by the red dots). We find that our estimated velocities for the low and
medium noise levels virtually coincide with the truth, again with very little
uncertainty. Remarkably, even for the higher noise level cases the velocities
estimated are good, $\chi^{2}=24.9$ and $\chi^{2}=20.4$ – mostly the outermost
spaxels are slightly off and the uncertainties increase (recall that S/N drops
down to $\sim 4$ or $\sim 1$ in the edge spaxels). Therefore, our method can
estimate the rotational velocities accurately for the emission line dominated
spectra even when the noise in the data is significant.

(a) low noise set (S/N $\gtrsim 59.06$) (b) medium noise set (S/N $\gtrsim
18.56$)

(c) high noise set (S/N $\gtrsim 3.97$) (d) very high noise set (S/N $\gtrsim
0.78$)

Figure 10: Velocity reconstruction for spaxels with emission line dominated
spectra, for different noise levels.

###  4.2 Absorption line dominated

Next we test the method on absorption line dominated spectral data with four
sets having different noise levels. We use the weak absorption line dominated
model spectrum from the $(32,37)$ spaxel of the MaNGA galaxy 8952–9102, shown
by orange colour in figure 3. The spectra in three spaxels (the two edge
spaxels and the central one) are shown in figure 11, for the low, medium,
high, and very high noise level cases. The right panel of figure 9 portrays
how the S/N of the spectra varies across the spaxels for the four sets; again,
S/N decays significantly at the outer-disk spaxels for the noisier sets,
reaching as low as $\sim 2$ and $\sim 0.7$ in the two highest noise cases.
Indeed, it is difficult to discern by eye any clear features in the bottom
panels of figure 11.

To provide more robust estimations for spectra with lower S/N, we use
$\Delta=4.0$ Å for analysing the medium and higher noise sets but keep
$\Delta=1.5$ Å for the low noise set. The rotational velocities estimated for
the four sets with different noise levels are compared with the truth (red
dots) in the panels of the figure 12. The method is again successful, even for
such low S/N, though the uncertainties increase noticeably in the higher noise
cases.

Figure 11: Absorption line dominated spectra are shown for four sets (used in
section 4.2) with different noise levels – low, medium, high, and very high –
from the top to the bottom panel. Each panel shows spectra in three spaxels,
two edge spaxels and the central one. S/N for these four sets are as low as
$37.39,~{}9.95,~{}2.05,~{}0.66$, as shown in the respective panels. Notice
that in the bottom two panels barely any feature is visible.

(a) low noise set (S/N $\gtrsim 37.39$) (b) medium noise set (S/N $\gtrsim
9.95$)

(c) high noise set (S/N $\gtrsim 2.05$) (d) very high noise set (S/N $\gtrsim
0.66$)

Figure 12: Velocity reconstructions for spaxels with absorption line dominated
spectra, with different noise levels.

Let us now focus on the results of the very high noise sets from the two
cases, with emission line dominated spectra (test demonstrated in section 4.1)
and with absorption line dominated spectra (test demonstrated in section 4.2).
In figure 13, we explicitly show the absolute values of the residuals of
velocities (for spaxels) as a function of S/N for these two sets. The red
vertical line represents S/N $=4$ below which the traditional template fitting
approach struggles to estimate rotational velocities. Note that the velocity
estimations, along with the uncertainties, have been shown in figures 10(d)
and 12(d) respectively for these two sets.

We find that our estimations continue to be reasonably accurate well beyond
S/N $=4$ for both sets. For the case with emission line dominated spectra, we
get residual $\lesssim 20$ km/s for this set when S/N $\sim 0.8$ only (due to
the strong emission lines barely emerging from the noisy continua in the
bottom panel of figure 8). Even for the set with absorption line dominated
spectra the results are accurate till S/N $\sim 1.3$, below which we may get
quite large deviations from the true velocities but still the results are
consistent with the truth within $2\sigma$ as evident from figure 12(d). As
demonstrated in appendix 5, our results for the absorption lines dominated
sets with four different noise levels are superior to the velocity estimations
using the traditional fitting approach based on Penalized Pixel-Fitting (pPXF)
(Cappellari & Emsellem, 2004; Cappellari, 2017) and MILES stellar spectral
library (Vazdekis et al., 2010).

Figure 13: The figure shows the absolute residuals as a function of S/N for
the very high noise sets from the tests with the emission line dominated
spectra (shown by blue stars) and absorption line dominated spectra (shown by
green dots). The corresponding velocity estimations are already shown in
figures 10(d) and 12(d) with the uncertainties. The red vertical line
represents S/N$=4$ below which the traditional template fitting approach
struggles to estimate rotational velocities.

##  5 Application to an observed MaNGA galaxy

Finally, we calculate the galaxy rotation curve for the MaNGA galaxy
7991–12701 itself. We estimate the velocity differences between spaxels along
a major axis on the IFU hexagon, from spaxel ($37,12$) to spaxel
($37,64$)666We find maximal variation in velocity difference along this major
axis. The emission line dominated spectra in some of the spaxels along this
axis have been shown in figures 14 and 2. for simplicity. The spectra have
been observed with good signal to noise ratio along this axis; S/N $\sim 47$
for the spectrum in the central spaxel (37,37) and it falls down to $\sim 14$
at the outskirts. We apply the same approach to the observed data as we had to
simulated data.

###  5.1 Velocity difference

An example of the actual spectra for a pair of spaxels, $A=(37,37)$ and
$B=(37,42)$, is plotted in figure 14, both are seen to be emission line
dominated with the former having a redder continuum.

Figure 14: Observed spectra from the spaxels (37,37) and (37,42) of the MaNGA
galaxy 7991-12701 are dominated by emission lines (the flux is measured in the
unit of [$10^{-17}$ erg/cm2/s/$\AA$]). We discard the data with $\lambda>9850$
Å, i.e. beyond the red vertical dashed line, for telluric contamination and
divide the rest of the data into 4 bins, shown by the vertical dashed lines.
Note that for these observed spectra the wavelength separation increases with
wavelength, leading to different bin widths.

In figure 15, we show the correlation coefficient as a function of the
velocity difference, $r(\Delta V$), between the pair of spaxels. The left and
right panels show $r_{A^{s}B}$ and $r_{B^{s}A}$ respectively for $4$ bins. For
this illustration we choose $\Delta=1.5$ Å and $N_{\rm it}=10$. In either
panel we find that $r$ for different bins exhibits global maxima at similar
values of $\Delta V$ (and $\Delta V_{BA}\approx-\Delta V_{AB}$). More
importantly, although bin 4 does not contain any strong feature, as evident
from figure 14, due to small spectral features $r$ shows a prominent peak at a
consistent $\Delta V$ in both panels of figure 15 for this bin. This
illustrates the key benefit of this method that it does not rely on strong
identifiable features; rather makes use of the whole spectra. Thus we have
$4\times 2=8$ ‘good’ estimations of $\Delta V$ from both the panels. We
estimate the velocity difference between these two spaxels by taking the mean
and standard deviation of these eight individual estimations, obtaining
$\Delta V_{AB}\approx 80.29\pm 9.82$ km/s.

(a) $37,37\bigotimes 37,42$ (b) $37,42\bigotimes 37,37$

Figure 15: The correlation coefficient $r(\Delta V)$ is plotted as a function
of velocity difference ($\Delta V$) between a pair of spaxels, $A=(37,37)$ and
$B=(37,42)$ from the MaNGA galaxy 7991-12701. The left and right panel show
$r_{A^{s}B}$ and $r_{B^{s}A}$ respectively for $4$ bins, calculated using (6).
All the bins have maxima in similar $\Delta V$ in either panel and the peak
values are the negative of each other, as expected.

In table 3, we present the estimations of $\Delta V$ for the same pair of
spaxels as in figure 15 but for a number of different $\Delta$, $N_{\rm it}$,
and number of bins. The $\Delta V$ estimations corresponding to different
values of these parameters are consistent with each other that demonstrates
the robustness of this approach. Since the spectra in the spaxels we consider
here have good S/N ($\gtrsim 14$), we adopt $\Delta=1.5$, $N_{\rm it}=10$ and
4 bins in the analysis of this MaNGA galaxy, same as what we used for the
simulations in section 3.

Note the value of $\Delta V$ obtained from the stellar velocity map by the
Marvin team is $85.89\pm 3.34$ km/s, consistent with our estimations. There is
nothing to say what the true answer is. However, the Marvin results come from
their full 2D analysis, giving extra stability to the results and enforcing
the full triangle equality (that the vector sum of velocity differences among
a triangle of spaxels is zero). Having demonstrated our correlation method
here, we plan to apply it to the 2D data in the follow-up paper.

$\Delta$ (in Å) | 4 bins | 8 bins
---|---|---
$N_{\rm it}=10$ | $N_{\rm it}=20$ | $N_{\rm it}=10$ | $N_{\rm it}=20$
$1.5$ | $80.3\pm 9.9$ | $78.8\pm 10.1$ | $83.5\pm 10.1$ | $83.2\pm 10.6$
$2.0$ | $84.0\pm 9.5$ | $82.5\pm 9.9$ | $82.0\pm 10.5$ | $84.0\pm 9.5$
$3.0$ | $85.1\pm 8.4$ | $85.7\pm 8.5$ | $84.9\pm 8.9$ | $85.2\pm 8.4$
$4.0$ | $82.8\pm 9.2$ | $83.7\pm 8.6$ | $84.7\pm 10.6$ | $85.4\pm 10.1$
Table 3: Velocity difference between MaNGA 7991-12701 spaxels (37,37) and
(37,42) for different values of smoothing scale $\Delta$, number of iterations
$N_{\rm it}$, and number of wavelength bins considered. The value quoted in
Marvin DR15 between these two spaxels is $85.89\pm 3.34$ km/s .

We list several more velocity differences from spaxel pairs of MaNGA galaxy
7991–12701 in table 4, along with the $\Delta V$ values obtained from the
stellar velocity map from Marvin.

Our estimations are somewhat consistent with the Marvin values but tend to be
slightly higher than that from the Marvin stellar map. However, note that two
estimations based on two different approaches have uncertainties of similar
order of magnitude. In the simulation studies our method estimated $\Delta V$
quite accurately, without a bias or underestimation of uncertainties, so it is
not clear what the true answer is. The full 2D data analysis in our follow-up
paper, now that we have demonstrated the method, will give a more parallel
comparison to Marvin.

Pair of spaxels (A,B) | Our estimation of $\Delta V_{AB}$ in km/s |  $\Delta V_{AB}$ in km/s from Marvin
---|---|---
($37$,$12$) , ($37$,$36$) | $137.45\pm 3.65$ | $130.45\pm 6.39$
($37$,$16$) , ($37$,$37$) | $153.84\pm 4.31$ | $152.60\pm 5.95$
($37$,$20$) , ($37$,$41$) | $223.99\pm 4.29$ | $216.36\pm 5.33$
($37$,$22$) , ($37$,$38$) | $167.38\pm 4.91$ | $171.15\pm 5.16$
($37,26$) , ($37,32$) | $32.43\pm 3.61$ | $39.54\pm 4.78$
($37,26$) , ($37,37$) | $135.69\pm 4.62$ | $128.94\pm 4.45$
($37$,$28$) , ($37$,$36$) | $100.48\pm 2.72$ | $98.94\pm 4.15$
($37,32$) , ($37,37$) | $98.62\pm 5.19$ | $89.40\pm 3.37$
($37,37$) , ($37,42$) | $80.29\pm 9.82$ | $85.89\pm 3.34$
($37,37$) , ($37,48$) | $143.79\pm 6.34$ | $131.70\pm 3.65$
($37$,$37$) , ($37$,$58$) | $171.08\pm 2.81$ | $166.68\pm 5.19$
($37$,$38$) , ($37$,$60$) | $155.67\pm 5.25$ | $151.91\pm 5.85$
($37$,$41$) , ($37$,$64$) | $90.19\pm 4.29$ | $86.77\pm 5.42$
Table 4: We compare our estimation of the velocity difference $\Delta V_{AB}$
between various pairs of spaxels $A$ and $B$ with that from Marvin DR15. Here
the smoothing scale $\Delta=1.5$ Å, $N_{\rm it}=10$, and we divided the
spectra into $4$ wavelength bins.

###  5.2 Constructing the galaxy rotation curve

Now we complete constructing the galaxy rotation curve for this MaNGA galaxy,
using the spectra in the spaxels along the major axis on the IFU hexagon, from
spaxel ($37,12$) to spaxel ($37,64$). Calculating $\Delta V$ for all the pairs
of spaxels along this major axis would be more computationally expensive than
warranted now for our purpose of simple demonstration. Thus we choose the 11
central spaxels as anchor spaxels (i.e. the $j$ spaxels of section 3.1.2), and
add 23 spaxels more over the range where the data quality is good to the $i$
spaxel set, giving a total of 34 spaxels in the HMC analysis.

In figure 16 we compare our estimated velocities (relative to the central
spaxel) with that from the Marvin for the selected spaxels. As noted in the
previous section, we find reasonable agreement between Marvin results and our
estimations, with some trend toward larger velocities. Again we find that the
uncertainties in our estimations are comparable to that from Marvin across the
spaxels; the errorbars are too small to be visible for most of the central
spaxels. Uncertainties in the outer disk region are larger, but we do not find
a decrease in velocity as Marvin does. Our forthcoming 2D analysis may shed
further light.

Figure 16: Our estimations (black dots with error bars) of velocity along the
line of sight of the spaxels on the chosen axis with respect to the central
spaxel $(37,37)$ of the MaNGA galaxy 7991–12701. The uncertainties are so
small that the errorbars are not visible for most of the spaxels. We compare
to the velocities given in the Marvin database (red dots with error bars).
Note that S/N $\approx 47.5$ at the centre and gradually falls down to $\sim
14$ at outskirts.

##  6 Conclusions and discussion

Galaxy rotation curves, or the internal dynamics in general, provide one of
the primary lines of evidence for the existence of dark matter. Besides merely
indicating the presence of dark matter, one can map the velocity-radius
relation $v(r)$ directly to the mass distribution with radius. Traditionally,
rotational velocities are obtained by first fitting a template to the spectra
of different regions, e.g. spaxels along a galaxy major axis, and then
calculating the Doppler shift between spaxels. The Doppler shifts are then
translated into the velocity differences between the pairs of spaxels. This
works, but has dependence on the template accuracy for that particular galaxy
type and spectral noise properties.

In this article, we present a novel and template-free method to calculate the
galaxy rotation curves based on cross-correlation between the spectra in the
IFU data, with each member of the comparison pair alternately smoothed to
avoid spurious features from noise. We demonstrate that one can achieve
accurate and precise velocity difference measurements with this method. We
then globally optimize the array of velocity difference through Hamiltonian
Monte Carlo to construct the full galaxy rotation curves for galaxies, with
spectra simulated in a variety of observational conditions.

The test conditions include emission line dominated spectra with red/blue
continuum, absorption line dominated spectra, switching between different
types of spectra, injection of various levels of noise, etc. The method
appears to be promising and the results are precise with small uncertainties
in most cases, in particular for spaxels with S/N $\gtrsim 4.0$. Remarkably,
even for very noisy spectra with S/N $\sim 1$ we are able to recover the true
rotational velocities with reasonable accuracy.

After these validation tests we apply the method to data from the observed
MaNGA galaxy 7991–12701. Our estimations of velocity differences (between
different pairs of spaxels) appear to be slightly higher than what was
obtained from the Marvin velocity field, but for most cases they are
consistent with each other. Comparing the 1D rotation curve (i.e. the
velocities of the spaxels with respect to the central one) between our method
and Marvin, we find that the results are mostly consistent with each other
(figure 16). However, for the outer disk our estimated rotation curve does not
show the slight reduction in rotational velocities seen in Marvin. The results
are sufficiently promising to pursue further.

One of the chief advantages of this method is that it does not rely on strong
identifiable features in the spectra. Rather, cross-correlation can utilize
many small features in the spectra, and in conjunction with smoothing, remain
fairly insensitive to noise. We have found good results even for quite noisy
spectra.

There are several aspects to pursue further, going beyond the present proof of
concept. We will enhance the pipeline to analyse the whole 2D IFU data of a
number of MaNGA galaxies, not just a 1D slice on a diameter/major axis. This
should also add further stability to the results through the triangle equality
within our HMC analysis. We adopted here a simplified algorithm to construct
the galaxy rotation curve for illustrating our results; the general conversion
from a well measured velocity field to constructing the galaxy rotation curve
has additional elements. We also plan to focus on low surface-brightness
galaxies where the traditional template fitting approach struggles, especially
at the outskirts. Further investigation will show the degree to which this new
method can improve upon, or nicely complement, traditional template fitting
approaches.

## Acknowledgement

The authors acknowledge that the high performance computing facility at the
Korea Institute of Science and Technology Information (KISTI), assignment no.
KSC-2020-CRE-0153, has been used in this project. SB and AS thank Adarsh
Ranjan for crucial helps and Alex G. Kim for useful discussions at different
phases of the project. S.B. also thanks Brian Cherinka and Maria Argudo-
Fernández for useful discussions and explanations regarding the MaNGA data.
YSA, EL, and KY were supported in part by the Energetic Cosmos Laboratory and
EL by the U.S. Department of Energy, Office of Science, Office of High Energy
Physics, under contract no. DE-AC02-05CH11231. KY thanks Yessenov Foundation
for funding his stay in South Korea. YSA and KY also thank KASI for
hospitality during the early work.

## Data availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

##  Appendix A Comparing with Penalized Pixel-Fitting (pPXF)

In this section we compare our results with that from the traditional fitting
approach based on Penalized Pixel-Fitting (pPXF) (Cappellari & Emsellem, 2004;
Cappellari, 2017), using the MILES templates from Vazdekis et al. (2010).
Since obtaining the fitting results for all the cases is beyond the scope of
this article, we restrict ourselves to the data sets with absorption line
dominated spectra. In particular, we consider the four data sets with
different noise levels as described in section 4.2. Our results for these
tests have been presented in figure 12 already.

(a) low noise set (S/N $\gtrsim 37.39$) (b) medium noise set (S/N $\gtrsim
9.95$)

(c) high noise set (S/N $\gtrsim 2.05$) (d) very high noise set (S/N $\gtrsim
0.66$)

Figure 17: Comparing the velocity reconstructions from our approach (blue) and
pPXF (black) for the absorption line dominated spectra data sets with
different noise levels (this test and our results are explained in section 4.2
in more details). Note that the errorbar from pPXF for the rightmost spaxel in
the right-bottom panel is unrealistically high ($\sim 9000$km/s, truncated by
hand in the plot) which indicates that the automated pPXF is unstable for such
low S/N cases. Also note that pPXF can overestimates the uncertainty on the
fit velocity as vividly evident in the low-noise case presented in the top-
left panel.

As evident from figure 17, the velocities estimated from pPXF are in good
agreement with the truths for all four sets, especially for lower noise sets
(higher S/N cases) presented in the top panels. However, for the higher noise
sets, pPXF results becomes less accurate which can be seen clearly in residual
plots presented in the top subplots of the bottom two panels. Nevertheless,
our estimations are superior than that from pPXF, as expected, in all four
cases. The uncertainties from pPXF are much larger than ours for all the sets,
although we find that pPXF overestimates the uncertainties for many spaxels
(especially for the low noise level cases in the top panels 777When we apply
pPXF on the base MaNGA spectrum that was used for the simulation and has
similar noise level as the low-noise case, we obtain velocity errors of the
same order ($\sim 3$ km/s) as in the top-left panel of figure 17.). Also, we
notice that the automated pPXF is not robust in the very high noise cases as
it suffers from occasional failings as well as produces unrealistically huge
errors sometimes (bottom-right panel).

We compare the accuracy of the velocities estimated from the two approaches in
terms of the following quantities

$\epsilon_{r}=\frac{1}{N_{\rm spax}}\sum_{i}\left|v_{\rm est}-v_{\rm
true}\right|~{},~{}~{}\epsilon_{p}=\frac{1}{N_{\rm
spax}}\sum_{i}\left|\frac{v_{\rm est}-v_{\rm true}}{v_{\rm
true}}\right|~{},~{}~{}b=\frac{1}{N_{\rm spax}}\sum_{i}\left(v_{\rm
est}-v_{\rm true}\right)\;,$ (8)

which are the average residual, the average percentage error and the bias
respectively. Table 5 compares the results of these two approaches for the
four sets in terms of $\epsilon_{r}$, $\epsilon_{p}$ and bias ($b$). While
both method shows insignificant bias, the precision (in terms of
$\epsilon_{r}$, $\epsilon_{p}$) of the spaxel cross-correlation approach is
found be always better than that of pPXF.

Set | Average residual ($\epsilon_{r}$) in km/s | Percentage error ($\epsilon_{p}$) | Bias ($b$) in km/s
---|---|---|---
This work | pPXF | This work | pPXF | This work | pPXF
Low noise | $0.16$ | $0.99$ | $0.1\%$ | $0.6\%$ | $1.37\cdot 10^{-3}$ | $-2.66\cdot 10^{-3}$
Medium noise | $0.72$ | $3.29$ | $0.42\%$ | $1.89\%$ | $-1.62\cdot 10^{-3}$ | $-2.21\cdot 10^{-3}$
High noise | $6.14$ | $9.20$ | $3.62\%$ | $5.33\%$ | $1.35\cdot 10^{-2}$ | $6.68\cdot 10^{-3}$
Very high noise | $25.24$ | $28.60$ | $14.52\%$ | $16.50\%$ | $2.87\cdot 10^{-2}$ | $3.42\cdot 10^{-2}$
Table 5: We compare the results of the two approaches – spaxel cross-
correlation vs pPXF – for the four sets with different noise levels in terms
of average residual, average percentage error and bias.

This exercise thus demonstrates that our approach based on spaxel cross-
correlation provides more precise estimation of velocities as it use the
information from all parts of the spectra. However, this comes at the expense
of moderate amount associated computation cost. pPXF takes roughly $0.5$
seconds on a 32-thread CPU to fit one MaNGA-like spectrum, whereas it takes
$\sim 30$ sec to obtain the velocity difference between two such spectra using
spaxel cross-correlation on the same computer. Nevertheless, the robustness
for low S/N spectra and higher precision in our approach possibly justify the
tread off in some particular use cases.

## References

  * Acquaviva et al. (2012) Acquaviva V., Gawiser E., Guaita L., 2012, in Tuffs R. J., Popescu C. C., eds, The Spectral Energy Distribution of Galaxies - SED 2011 Vol. 284, The Spectral Energy Distribution of Galaxies - SED 2011. pp 42–45 (arXiv:1111.4243), doi:10.1017/S1743921312008691
  * Aghamousa & Shafieloo (2015) Aghamousa A., Shafieloo A., 2015, Astrophys. J., 804, 39
  * Amram et al. (1992) Amram P., Le Coarer E., Marcelin M., Balkowski C., Sullivan W. T. I., Cayatte V., 1992, A&AS, 94, 175
  * Aurenhammer (1991) Aurenhammer F., 1991, ACM Comput. Surv., 23, 345–405
  * Boquien et al. (2019) Boquien M., Burgarella D., Roehlly Y., Buat V., Ciesla L., Corre D., Inoue A. K., Salas H., 2019, Astron. Astrophys., 622, A103
  * Bosma (1981) Bosma A., 1981, Astron. J., 86, 1825
  * Bryant et al. (2016) Bryant J. J., et al., 2016, in Evans C. J., Simard L., Takami H., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 9908, Ground-based and Airborne Instrumentation for Astronomy VI. p. 99081F (arXiv:1608.03921), doi:10.1117/12.2230740
  * Bundy et al. (2015) Bundy K., et al., 2015, Astrophys. J., 798, 7
  * Cappellari (2017) Cappellari M., 2017, Mon. Not. Roy. Astron. Soc., 466, 798
  * Cappellari & Copin (2003) Cappellari M., Copin Y., 2003, Mon. Not. Roy. Astron. Soc., 342, 345
  * Cappellari & Emsellem (2004) Cappellari M., Emsellem E., 2004, Publ. Astron. Soc. Pac, 116, 138
  * Cappellari et al. (2013) Cappellari M., et al., 2013, Mon. Not. Roy. Astron. Soc., 432, 1709
  * Carpenter et al. (2017) Carpenter B., Gelman A., Hoffman M. D., Lee D., Goodrich B., Betancourt M., et al.. 2017, Journal of Statistical Software, 76, 1
  * Catinella et al. (2006) Catinella B., Giovanelli R., Haynes M. P., 2006, Astrophys. J., 640, 751
  * Cazzoli et al. (2016) Cazzoli S., Arribas S., Maiolino R., Colina L., 2016, Astron. Astrophys., 590, A125
  * Cherinka et al. (2019) Cherinka B., et al., 2019, Astron. J., 158, 74
  * Cid Fernandes et al. (2005) Cid Fernandes R., Mateus A., Sodré L., Stasińska G., Gomes J. M., 2005, Mon. Not. Roy. Astron. Soc., 358, 363
  * Cole et al. (1994) Cole S., Aragon-Salamanca A., Frenk C. S., Navarro J. F., Zepf S. E., 1994, Mon. Not. Roy. Astron. Soc., 271, 781
  * Couto et al. (2021) Couto G. S., Hughes T. M., Boquien M., Ibar E., Viaene S., Leiton R., Xue Y., 2021, Astron. Astrophys., 654, A128
  * Croom et al. (2021) Croom S. M., et al., 2021, Mon. Not. Roy. Astron. Soc., 505, 991
  * Das et al. (2020) Das M., McGaugh S. S., Ianjamasimanana R., Schombert J., Dwarakanath K. S., 2020, Astrophys. J., 889, 10
  * Drory et al. (2015) Drory N., et al., 2015, Astron. J., 149, 77
  * Erroz-Ferrer et al. (2015) Erroz-Ferrer S., et al., 2015, Mon. Not. Roy. Astron. Soc., 451, 1004
  * Ewen & Purcell (1951) Ewen H. I., Purcell E. M., 1951, Nature, 168, 356
  * Faber & Jackson (1976) Faber S. M., Jackson R. E., 1976, Astrophys. J., 204, 668
  * Ferreras et al. (2019) Ferreras I., et al., 2019, Mon. Not. Roy. Astron. Soc., 489, 608
  * García-Lorenzo et al. (2015) García-Lorenzo B., et al., 2015, Astron. Astrophys., 573, A59
  * Ge et al. (2021) Ge J., Mao S., Lu Y., Cappellari M., Long R. J., Yan R., 2021, Mon. Not. Roy. Astron. Soc., 507, 2488
  * Giovanelli & Haynes (2002) Giovanelli R., Haynes M. P., 2002, Astrophys. J., 571, L107
  * Gudehus (1973) Gudehus D. H., 1973, Astron. J., 78, 583
  * Johnson et al. (2021) Johnson B. D., Leja J., Conroy C., Speagle J. S., 2021, Astrophys. J. Suppl., 254, 22
  * Koleva et al. (2009) Koleva M., Prugniel P., Bouchard A., Wu Y., 2009, Astron. Astrophys., 501, 1269
  * Kretschmer et al. (2021) Kretschmer M., Dekel A., Freundlich J., Lapiner S., Ceverino D., Primack J., 2021, Mon. Not. Roy. Astron. Soc., 503, 5238
  * Loubser et al. (2020) Loubser S. I., Babul A., Hoekstra H., Bahé Y. M., O’Sullivan E., Donahue M., 2020, Mon. Not. Roy. Astron. Soc., 496, 1857
  * Márquez et al. (2004) Márquez I., et al., 2004, Astron. Astrophys., 416, 475
  * Ocvirk et al. (2006) Ocvirk P., Pichon C., Lançon A., Thiébaut E., 2006, Mon. Not. Roy. Astron. Soc., 365, 46
  * Riera et al. (1988) Riera A., Mampaso A., Vilchez J. M., Manchado A., Phillips J. P., 1988, in Bianchi L., Gilmozzi R., eds, Mass Outflows from Stars and Galactic Nuclei Vol. 142, Mass Outflows from Stars and Galactic Nuclei. p. 285, doi:10.1007/978-94-009-2941-8_32
  * Roberts & Rots (1973) Roberts M. S., Rots A. H., 1973, Astron. Astrophys., 26, 483
  * Roberts-Borsani et al. (2020) Roberts-Borsani G. W., Saintonge A., Masters K. L., Stark D. V., 2020, Mon. Not. Roy. Astron. Soc., 493, 3081
  * Rubin & Ford (1970) Rubin V. C., Ford W. Kent J., 1970, Astrophys. J., 159, 379
  * Rubin et al. (1978) Rubin V. C., Ford W. K. J., Thonnard N., 1978, Astrophys. J., 225, L107
  * Rubin et al. (1980) Rubin V. C., Ford W. K. J., Thonnard N., 1980, Astrophys. J., 238, 471
  * Sánchez et al. (2012) Sánchez S. F., et al., 2012, Astron. Astrophys., 538, A8
  * Sánchez et al. (2016a) Sánchez S. F., et al., 2016a, Rev. Mex. Astron. Astrofis., 52, 171
  * Sánchez et al. (2016b) Sánchez S. F., et al., 2016b, Astron. Astrophys., 594, A36
  * Scott et al. (2018) Scott N., et al., 2018, Mon. Not. Roy. Astron. Soc., 481, 2299
  * Shafieloo (2007) Shafieloo A., 2007, Mon. Not. Roy. Astron. Soc., 380, 1573
  * Shafieloo & Clarkson (2010) Shafieloo A., Clarkson C., 2010, Phys. Rev. D, 81, 083537
  * Shafieloo et al. (2006) Shafieloo A., Alam U., Sahni V., Starobinsky A. A., 2006, Mon. Not. Roy. Astron. Soc., 366, 1081
  * Sofue & Rubin (2001a) Sofue Y., Rubin V., 2001a, Ann. Rev. Astron. Astrophys., 39, 137
  * Sofue & Rubin (2001b) Sofue Y., Rubin V., 2001b, Ann. Rev. Astron. Astrophys., 39, 137
  * Stan Development Team (2017) Stan Development Team 2017, PyStan: the Python interface to Stan, http://mc-stan.org
  * Stoehr et al. (2008) Stoehr F., et al., 2008, in Argyle R. W., Bunclark P. S., Lewis J. R., eds, Astronomical Society of the Pacific Conference Series Vol. 394, Astronomical Data Analysis Software and Systems XVII. p. 505
  * Storchi-Bergmann (2010) Storchi-Bergmann T., 2010, arXiv e-prints, p. arXiv:1001.2480
  * Thomas et al. (2011) Thomas J., et al., 2011, Mon. Not. Roy. Astron. Soc., 415, 545
  * Torres-Flores et al. (2011) Torres-Flores S., Epinat B., Amram P., Plana H., Mendes de Oliveira C., 2011, Mon. Not. Roy. Astron. Soc., 416, 1936
  * Tully & Fisher (1977) Tully R. B., Fisher J. R., 1977, Astron. Astrophys., 500, 105
  * Vazdekis et al. (2010) Vazdekis A., Sánchez-Blázquez P., Falcón-Barroso J., Cenarro A. J., Beasley M. A., Cardiel N., Gorgas J., Peletier R. F., 2010, Mon. Not. Roy. Astron. Soc., 404, 1639
  * Vogt et al. (1993) Vogt N. P., Herter T., Haynes M. P., Courteau S., 1993, Astrophys. J., 415, L95
  * Wake et al. (2017) Wake D. A., et al., 2017, Astron. J., 154, 86
  * Walcher et al. (2006) Walcher C. J., Boker T., Charlot S., Ho L. C., Rix H.-W., Rossa J., Shields J. C., van der Marel R. P., 2006, Astrophys. J., 649, 692
  * Yan et al. (2016) Yan R., et al., 2016, Astron. J., 152, 197
  * Zwicky (1933) Zwicky F., 1933, Helvetica Physica Acta, 6, 110
  * de Blok et al. (2008) de Blok W. J. G., Walter F., Brinks E., Trachternach C., Oh S. H., Kennicutt R. C. J., 2008, Astron. J., 136, 2648
  * van Albada et al. (1985) van Albada T. S., Bahcall J. N., Begeman K., Sancisi R., 1985, Astrophys. J., 295, 305


$d=4+\varepsilon$ | $\dfrac{g^{2}+(M^{2}+\mu^{2})^{2}}{16\pi^{2}}$ | 0
$d=6+\varepsilon$ | $-\dfrac{(M^{2}+\mu^{2})^{3}+g^{2}(3M^{2}+\mu^{2})}{192\pi^{3}}$ | $\dfrac{g^{2}}{96\pi^{3}}\mu^{2}$
$d=8+\varepsilon$ | $\dfrac{5g^{4}+5(M^{2}+\mu^{2})^{4}+g^{2}(30M^{4}+20M^{2}\mu^{2}+6\mu^{4})}{15360\pi^{4}}$ | $-\dfrac{g^{2}(M^{2}+\mu^{2})}{384\pi^{4}}\mu^{2}+\dfrac{g^{2}}{960\pi^{4}}\mu^{4}$
Table 1: Divergent terms in dimensional regularization in even dimensions
$d=2n+\varepsilon$ (with $n>1$), and their $\mu$ dependent part. See eqs.
(182) and (183) for the expression of $M^{2}$ and $g^{2}$ in terms of the
interaction potential $V_{\rm int}(\phi)$ and their $\mu$ dependence. Figure
3:  One loop Feynman diagrams up to four interaction points. Diagrams $a_{0})$
and $a_{1})$ are divergent in dimension $d>2$, but the corresponding
divergence is always independent of external momenta. Diagram $a_{n})$ is
marginally (logarithmically) divergent in $d=2n$ dimensions. Derivative
counterterms in the $\mu=0$ theory arise from diagrams $a_{n})$ with $n\geq 2$
when they are at least quadratically divergent. This corresponds to $d\geq 6$.

In order to understand this, consider the $\mu=0$ theory. We take $V_{\rm
int}(\phi)$ to be an arbitrary ${\rm U}(1)$ invariant polynomial potential,
built from even powers of $\phi$. Since we are considering (possibly) non-
renormalizable interactions, the counterterms needed to renormalize the theory
will in general include derivative operators, even if the interaction terms we
started from are non-derivative. Introducing a chemical potential is
equivalent to promoting ordinary derivatives to $\mu$ dependent covariant
derivatives, so that derivative counterterms will generate $\mu$ dependent
counterterms as well.212121We are grateful to Riccardo Rattazzi for a
discussion on this point. By a simple power counting argument it is easy to
see that no derivative counterterms are needed at one loop in $d=4$, for
arbitrary interactions. In order to see this notice that the (amputated)
diagrams which depend on external momenta, $a_{n})$ with $n\geq 2$ in Fig. 3,
are convergent or at most logarithmically divergent in $d=4$. As a
consequence, the corresponding divergencies are independent of external
momenta and can be renormalized by non-derivative polynomial counterterms.
This is a generalization of the familiar statement that at one loop there is
no wave-function renormalization in $d=4$. In the case of even ($\rm U(1)$
symmetric) interactions the wave-function renormalization is still absent at
one-loop for $d=2n$, but starting from $d=6$ derivative interaction
counterterms are needed. Notice that interaction terms are always irrelevant
in $d\geq 6$, so that the theories we consider are non-renormalizable and the
derivative operators we are discussing are always generated.

In order to check our diagrammatic argument, consider the example of the
$\lambda\phi^{4}$ theory in $d=6$. The $\mu$ dependent divergence in the one-
loop effective potential can be read from Tab. 1 and is given by

$\Xi^{(1)}(\phi;\mu)\Big{|}_{d,{\rm
div}}=\dfrac{1}{\varepsilon}\dfrac{\lambda^{2}}{24\pi^{3}}\mu^{2}\phi^{4}.$
(189)

On the other hand, in the $\mu=0$ theory, the only diagram that can have
divergencies dependent on external momenta is $a_{2})$ in Fig. 3. Computing
the divergent part of this diagram we find that the following derivative
counterterm is needed for $\mu=0$:

$\mathcal{L}_{\rm
ct}\supset-\dfrac{1}{\varepsilon}\dfrac{\lambda^{2}}{24\pi^{3}}(\partial_{\nu}\Phi)^{\dagger}(\partial^{\nu}\Phi)\Phi^{\dagger}\Phi,$
(190)

where we only wrote the divergent part. Promoting derivatives to $\mu$
dependent covariant derivatives $D_{\nu}=\partial_{\nu}-i\mu\delta_{\nu}^{0}$,
we find that a $\mu$ dependent counterterm is generated that exactly cancels
the divergence of eq. (189).

As a last remark, we note that the classic argument on the non-renormalization
of conserved currents (see _e.g._ Collins:1984xc , p. 162) advocated in
Benson:1991nj  to justify the absence of $\mu$ dependent divergencies in $d=4$
is not applicable at finite $\mu$. In fact, at finite chemical potential, an
identically conserved current can be constructed with the aid of the object
$\mu\delta_{\nu}^{0}$, invalidating the non-renormalization theorem. The
current is simply given by $\tilde{J}^{\nu}\leavevmode\nobreak\
=\leavevmode\nobreak\
(\partial^{\nu}\partial_{0}-\delta^{\nu}_{0}\square)\phi^{2}$.

##  Appendix F More details on the derivation of the scaling relations of
section 5

In this appendix the variable $X$ will be always set at its background value
$X\rightarrow\bar{X}=\mu^{2}$. We start from the relations
$f^{2}_{\pi}(Q)=Q/\mu=2P^{\prime}(X)$ and $Q=2P^{\prime}(X)\sqrt{X}$, and the
sound speed

$c_{s}^{2}=\dfrac{P^{\prime}(X)}{P^{\prime}(X)+2P^{\prime\prime}(X)X}.$ (191)

The derivatives of $f_{\pi}^{2}$ are:

$\displaystyle\dfrac{{\rm d}f_{\pi}^{2}}{{\rm d}Q}=\dfrac{{\rm
d}(2P^{\prime}(X))}{{\rm d}X}\dfrac{{\rm d}X}{{\rm
d}Q}=2P^{\prime\prime}(X)\dfrac{1}{2P^{\prime\prime}(X)\sqrt{X}+P^{\prime}(X)/\sqrt{X}}=\dfrac{1-c_{s}^{2}}{\sqrt{X}},$
(192) $\displaystyle\dfrac{{\rm d^{2}}f_{\pi}^{2}}{{\rm d}Q^{2}}=-\dfrac{{\rm
d}c_{s}^{2}}{{\rm d}Q}\dfrac{1}{\sqrt{X}}+(1-c_{s}^{2})\dfrac{{\rm d}}{{\rm
d}X}\left(\dfrac{1}{\sqrt{X}}\right)\dfrac{{\rm d}X}{{\rm d}Q}=-\dfrac{{\rm
d}c_{s}^{2}}{{\rm
d}Q}\dfrac{1}{\sqrt{X}}-\dfrac{c_{s}^{2}}{Q}\dfrac{(1-c_{s}^{2})}{\sqrt{X}}.$
(193)

The derivative of the sound speed is instead:

$\dfrac{{\rm d}c_{s}^{2}}{{\rm d}Q}=\dfrac{{\rm d}c_{s}^{2}}{{\rm
d}X}\dfrac{{\rm d}X}{{\rm d}Q}=\dfrac{{\rm d}c_{s}^{2}}{{\rm
d}X}\dfrac{1}{2P^{\prime\prime}(X)\sqrt{X}+P^{\prime}(X)/\sqrt{X}}.$ (194)

In the limit $Q\rightarrow 0$, using that $2P^{\prime}(X)=\mathcal{O}(Q)$, we
obtain:

$\dfrac{{\rm d}c_{s}^{2}}{{\rm
d}X}=\dfrac{P^{\prime\prime}(X)}{P^{\prime}(X)+2P^{\prime\prime}(X)X}+\mathcal{O}(Q)=\dfrac{1}{2X}+\mathcal{O}(Q)\implies\dfrac{{\rm
d}c_{s}^{2}}{{\rm d}Q}=\dfrac{1}{4X^{3/2}P^{\prime\prime}(X)}+\mathcal{O}(Q).$
(195)

Using that for $Q\rightarrow 0$, $\sqrt{X}=\mu_{\rm
crit}+\mathcal{O}(Q)=m_{\rm pole}+\mathcal{O}(Q)$ we arrive at the result:

$\dfrac{{\rm d}c_{s}^{2}}{{\rm d}Q}=\dfrac{1}{4m_{\rm
pole}^{3/2}P^{\prime\prime}(m_{\rm pole}^{2})}+\mathcal{O}(Q).$ (196)

Higher order terms can be straightforwardly computed by rewriting $c_{s}^{2}$
as

$c_{s}^{2}=\dfrac{Q}{Q+4P^{\prime\prime}(\mu^{2})\mu^{3}},$ (197)

and then expressing $\mu(Q)$ recursively through

$\mu(Q)=m_{\rm pole}\exp\left(\int_{0}^{Q}\frac{c_{s}^{2}}{Q^{\prime}}{\rm
d}Q^{\prime}\right).$ (198)

##  Appendix G Ground state and stability for the ${\rm O}(N)$ model with
quartic interactions

We want to study the minimum of the large $N$ effective potential in the ${\rm
O}(N)$ model with quartic interactions in $d=3$, given by eq. (125). The
conditions for minimization are given in eqs. (126). First, it is easy to see
that the minimum of the effective potential always corresponds to $S_{i}=0$.
Moreover it is convenient to define the $\rm SO(2)$ invariant variable
$t=(s_{1}^{2}+s_{2}^{2})/N$. The stationary points of $V_{\rm eff}(s_{i})$ as
a function of $s_{i}$ correspond through the chain rule to either $t=0$ or to
real stationary points of $V_{\rm eff}(t)$, satisfying $V^{\prime}_{\rm
eff}(t)=0$. As already remarked in Coleman:1974jh , it can be convenient to
solve the condition (126a) for $\chi$ and plug back in eq. (125) to study the
ordinary effective potential, without the auxiliary parameter $\chi$. From eq.
(126a) we have:

$f(t)\equiv\sqrt{\chi}=\dfrac{\lambda_{N}}{8\pi}\left(\sqrt{1+\dfrac{64\pi^{2}}{\lambda_{N}}\left(t+\dfrac{m^{2}}{\lambda_{N}}\right)}-1\right),$
(199)

where we selected the solution with $\sqrt{\chi}>0$. We obtain

$\dfrac{V_{\rm
eff}(t)}{N}=\dfrac{1}{2}\left(f^{2}-\mu^{2}\right)t-\dfrac{\left(f^{2}-m^{2}\right)^{2}}{4\lambda_{N}}-\dfrac{f^{3}}{12\pi},$
(200)

from which it follows

$\dfrac{V^{\prime}_{\rm
eff}(t)}{N}=\dfrac{1}{2}\left(f^{2}-\mu^{2}\right)+ff^{\prime}t-\dfrac{\left(f^{2}-m^{2}\right)}{\lambda_{N}}ff^{\prime}-\dfrac{f^{2}f^{\prime}}{4\pi}=\dfrac{1}{2}\left(f^{2}-\mu^{2}\right),$
(201)

where we suppressed the explicit $t$ dependence in $f(t)$ and used
$f^{2}-m^{2}=\lambda_{N}t-\lambda_{N}f/4\pi$. The stationary points of $V_{\rm
eff}(t)$ correspond to real (hence physical) values of $t$ only when

$\mu^{2}\geq\mu_{\rm crit}^{2}=f(0)^{2},$ (202)

from which eq. (128) follows. For values of $\mu^{2}\leq\mu_{\rm crit}^{2}$
the potential $V_{\rm eff}(s_{i})$ is minimized for $t=0$, corresponding to
the unbroken phase. Moreover we see that in $d=3$ and for $\mu^{2}>\mu_{\rm
crit}^{2}$ the potential has only one stationary point for $t>0$ (a minimum)
and is always bounded from below, thus stable. The situation is different in
$d=4$, where due to logarithmic terms the potential is unbounded from below,
signalling an instability Gross:1974jv ; Coleman:1974jh .

Notice also that by rescaling eq. (201) we can also derive the value of the
pole mass in the unbroken phase of the UV theory at $\mu=0$. This is given by

$m^{2}_{\rm pole}=\dfrac{2}{N}V^{\prime}_{\rm
eff}(t;\mu=0)=f(0)^{2}=\left(\sqrt{\dfrac{\lambda_{N}^{2}}{64\pi^{2}}+m^{2}}-\dfrac{\lambda_{N}}{8\pi}\right)^{2}.$
(203)

In particular, it follows that $\mu_{\rm crit}^{2}=m_{\rm pole}^{2}$. This
result is non-perturbative in $\lambda_{N}$, at leading order in the $1/N$
expansion.

## References

  * (1) A. Nicolis and F. Piazza, “Spontaneous Symmetry Probing,” JHEP 06 (2012) 025, arXiv:1112.5174 [hep-th].
  * (2) R. Shankar, “Renormalization group approach to interacting fermions,” Rev. Mod. Phys. 66 (1994) 129–192, arXiv:cond-mat/9307009.
  * (3) D. Vollhardt and P. Wolfle, The superfluid phases of helium 3.  Dover Publications, 2013.
  * (4) S. Minwalla, A. Mishra, and N. Prabhakar, “Fermi seas from Bose condensates in Chern-Simons matter theories and a bosonic exclusion principle,” JHEP 11 (2020) 171, arXiv:2008.00024 [hep-th].
  * (5) M. Geracie, M. Goykhman, and D. T. Son, “Dense Chern-Simons Matter with Fermions at Large N,” JHEP 04 (2016) 103, arXiv:1511.04772 [hep-th].
  * (6) R. Ciccone, L. Di Pietro, and M. Serone, “Inhomogeneous Phase of the Chiral Gross-Neveu Model,” Phys. Rev. Lett. 129 no. 7, (2022) 071603, arXiv:2203.07451 [hep-th].
  * (7) M. Thies, “From relativistic quantum fields to condensed matter and back again: Updating the Gross-Neveu phase diagram,” J. Phys. A 39 (2006) 12707–12734, arXiv:hep-th/0601049.
  * (8) S. Hellerman, D. Orlando, S. Reffert, and M. Watanabe, “On the CFT Operator Spectrum at Large Global Charge,” JHEP 12 (2015) 071, arXiv:1505.01537 [hep-th].
  * (9) F. Strocchi, Symmetry Breaking.  Springer Berlin, Heidelberg, 2008.
  * (10) A. Joyce, A. Nicolis, A. Podo, and L. Santoni, “Integrating out beyond tree level and relativistic superfluids,” JHEP 09 (2022) 066, arXiv:2204.03678 [hep-th].
  * (11) A. Nicolis, A. Podo, and L. Santoni. Work in progress.
  * (12) N. Dondi, S. Hellerman, I. Kalogerakis, R. Moser, D. Orlando, and S. Reffert, “Fermionic CFTs at large charge and large N,” arXiv:2211.15318 [hep-th].
  * (13) J. I. Kapusta, “Bose-Einstein Condensation, Spontaneous Symmetry Breaking, and Gauge Theories,” Phys. Rev. D 24 (1981) 426–439.
  * (14) J. Bernstein and S. Dodelson, “Relativistic Bose gas,” Phys. Rev. Lett. 66 (1991) 683–687.
  * (15) K. M. Benson, J. Bernstein, and S. Dodelson, “Phase structure and the effective potential at fixed charge,” Phys. Rev. D 44 (1991) 2480–2497.
  * (16) T. Brauner, “Spontaneous symmetry breaking in the linear sigma model at finite chemical potential: One-loop corrections,” Phys. Rev. D 74 (2006) 085010, arXiv:hep-ph/0607102.
  * (17) A. Sharma, G. Kartvelishvili, and J. Khoury, “Finite temperature description of an interacting Bose gas,” Phys. Rev. D 106 no. 4, (2022) 045025, arXiv:2204.02423 [hep-th].
  * (18) R. Brunetti, K. Fredenhagen, and N. Pinamonti, “Algebraic approach to Bose-Einstein Condensation in relativistic Quantum Field Theory. Spontaneous symmetry breaking and the Goldstone Theorem,” Annales Henri Poincare 22 no. 3, (2021) 951–1000, arXiv:1911.01829 [math-ph].
  * (19) S. R. Coleman, “There are no Goldstone bosons in two-dimensions,” Commun. Math. Phys. 31 (1973) 259–264.
  * (20) N. D. Mermin and H. Wagner, “Absence of ferromagnetism or antiferromagnetism in one-dimensional or two-dimensional isotropic Heisenberg models,” Phys. Rev. Lett. 17 (1966) 1133–1136.
  * (21) A. Nicolis and A. Podo. Work in progress.
  * (22) G. Badel, G. Cuomo, A. Monin, and R. Rattazzi, “Feynman diagrams and the large charge expansion in $3-\varepsilon$ dimensions,” Phys. Lett. B 802 (2020) 135202, arXiv:1911.08505 [hep-th].
  * (23) G. Jona-Lasinio, “Relativistic field theories with symmetry breaking solutions,” Nuovo Cim. 34 (1964) 1790–1795.
  * (24) S. R. Coleman and E. J. Weinberg, “Radiative Corrections as the Origin of Spontaneous Symmetry Breaking,” Phys. Rev. D 7 (1973) 1888–1910.
  * (25) R. Jackiw, “Functional evaluation of the effective potential,” Phys. Rev. D 9 (1974) 1686.
  * (26) S. Weinberg, The quantum theory of fields. Vol. 2: Modern applications.  Cambridge University Press, 8, 2013.
  * (27) D. T. Son, “Low-energy quantum effective action for relativistic superfluids,” arXiv:hep-ph/0204199.
  * (28) P. Creminelli, O. Janssen, and L. Senatore, “Positivity bounds on effective field theories with spontaneously broken Lorentz invariance,” JHEP 09 (2022) 201, arXiv:2207.14224 [hep-th].
  * (29) S. R. Coleman, R. Jackiw, and H. D. Politzer, “Spontaneous Symmetry Breaking in the O(N) Model for Large N*,” Phys. Rev. D 10 (1974) 2491.
  * (30) D. J. Gross and A. Neveu, “Dynamical Symmetry Breaking in Asymptotically Free Field Theories,” Phys. Rev. D 10 (1974) 3235.
  * (31) L. F. Abbott, J. S. Kang, and H. J. Schnitzer, “Bound States, Tachyons, and Restoration of Symmetry in the 1/N Expansion,” Phys. Rev. D 13 (1976) 2212.
  * (32) J. Zinn-Justin, Quantum field theory and critical phenomena, vol. 77 of International Series of Monographs on Physics.  Oxford University Press, 4, 2021.
  * (33) L. Alvarez-Gaume, D. Orlando, and S. Reffert, “Large charge at large N,” JHEP 12 (2019) 142, arXiv:1909.02571 [hep-th].
  * (34) D. Orlando, S. Reffert, and T. Schmidt, “Following the flow for large N and large charge,” Phys. Lett. B 825 (2022) 136881, arXiv:2110.07616 [hep-th].
  * (35) S. Giombi, E. Helfenberger, and H. Khanchandani, “Long range, large charge, large N,” JHEP 01 (2023) 166, arXiv:2205.00500 [hep-th].
  * (36) A. Monin, D. Pirtskhalava, R. Rattazzi, and F. K. Seibold, “Semiclassics, Goldstone Bosons and CFT data,” JHEP 06 (2017) 011, arXiv:1611.02912 [hep-th].
  * (37) G. Badel, G. Cuomo, A. Monin, and R. Rattazzi, “The Epsilon Expansion Meets Semiclassics,” JHEP 11 (2019) 110, arXiv:1909.01269 [hep-th].
  * (38) L. A. Gaumé, D. Orlando, and S. Reffert, “Selected topics in the large quantum number expansion,” Phys. Rept. 933 (2021) 1–66, arXiv:2008.03308 [hep-th].
  * (39) T. D. Lee and C. N. Yang, “Many-Body Problem in Quantum Mechanics and Quantum Statistical Mechanics,” Phys. Rev. 105 (1957) 1119–1120.
  * (40) T. D. Lee, K. Huang, and C. N. Yang, “Eigenvalues and Eigenfunctions of a Bose System of Hard Spheres and Its Low-Temperature Properties,” Phys. Rev. 106 (1957) 1135–1145.
  * (41) E. Braaten and A. Nieto, “Renormalization effects in a dilute Bose gas,” Phys. Rev. B 55 (1997) 8090–8093, arXiv:hep-th/9609047.
  * (42) E. Braaten and A. Nieto, “Quantum corrections to the ground state of a trapped Bose-Einstein condensate,” Phys. Rev. B 56 (1997) 14745–14765, arXiv:cond-mat/9707199.
  * (43) S. Fournais and J. P. Solovej, “The energy of dilute bose gases,” Annals of Mathematics 192 no. 3, (2020) 893–976.
  * (44) S. Fournais and J. P. Solovej, “The energy of dilute bose gases ii: The general case,” Inventiones mathematicae (2022) 1–132.
  * (45) H. Matsumoto, H. Umezawa, N. Yamamoto, and N. J. Papastamatiou, “Nonparticle Like Goldstone Modes in the Spontaneous Breakdown of Symmetry,” Phys. Rev. D 34 (1986) 3217–3222.
  * (46) L. Alberte and A. Nicolis, “Spontaneously broken boosts and the Goldstone continuum,” JHEP 07 (2020) 076, arXiv:2001.06024 [hep-th].
  * (47) Z. Komargodski, M. Mezei, S. Pal, and A. Raviv-Moshe, “Spontaneously broken boosts in CFTs,” JHEP 09 (2021) 064, arXiv:2102.12583 [hep-th].
  * (48) A. Nicolis and F. Piazza, “Implications of Relativity on Nonrelativistic Goldstone Theorems: Gapped Excitations at Finite Charge Density,” Phys. Rev. Lett. 110 no. 1, (2013) 011602, arXiv:1204.1570 [hep-th]. [Addendum: Phys.Rev.Lett. 110, 039901 (2013)].
  * (49) J. C. Collins, Renormalization: An Introduction to Renormalization, The Renormalization Group, and the Operator Product Expansion, vol. 26 of Cambridge Monographs on Mathematical Physics.  Cambridge University Press, Cambridge, 1986.