\documentclass{article}
\usepackage[affil-it]{authblk}
\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
\usepackage{url}
\usepackage{hyperref}
\hypersetup{colorlinks=false,pdfborder={0 0 0}}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
\errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}
\usepackage[utf8]{inputenc}
\usepackage[english]{babel}
\iflatexml
\documentclass{emulateapj}
\fi
\usepackage{amsmath}
\newcommand{\HII}{\ensuremath{\mathrm{H}\,\mathrm{II}}}
\newcommand{\HI}{\ensuremath{\mathrm{H}\, \mathrm{I}}}
\newcommand{\Hmol}{\ensuremath{\mathrm{H}_2}}
\newcommand{\glon}{\ensuremath{\ell}}
\newcommand{\glat}{\ensuremath{\mathrm{b}}}
\newcommand{\RGC}{\ensuremath{R_{\mathrm{GC}}}}
\newcommand{\vlos}{\ensuremath{v_{\mathrm{los}}}}
\newcommand{\deg}{\ensuremath{^\circ}}
\begin{document}
\title{Kinetic Tomography I: (deprecated)}
\author{Kirill Tchernyshyov}
\affil{Affiliation not available}
\author{Josh Peek}
\affil{Space Telescope Science Institute}
\date{\today}
\maketitle
\selectlanguage{english}
\begin{abstract}
text%
\end{abstract}%
\section{Introduction}
Many open problems in star formation, molecular cloud evolution, and galaxy-scale gas dynamics remain open because it has not been possible to measure the most useful quantities for resolving them -- the 3D gas velocity vector and 3D gas density over an extended area of sky. A measurement of these fields would allow us to solve the continuity equation \cite{euler1757principes}, and measure the rate at which density is changing in any region of the Galaxy, and over what physical scale.
Flows of converging gas are a central part of theories of the formation of giant molecular clouds (GMCs) \cite{Vazquez_Semadeni_2007,Audit_2005}. Suggested mechanisms for forming GMCs out of diffuse gas include different sorts of thermodynamic, hydrodynamic, and magnetohydrodynamic instabilities, particularly in colliding flows \citep{Clark:2012bq,2014ApJ...790...37C,Heitsch06}; gas compression driven by large-scale structures such as spiral arms and expanding supergiant shells \citep{Roberts:1972bp,Bonnell:2006hn,Fujimoto:2014kh}; and collapse due to self-gravity \citep{Kim:2002da,2012MNRAS.425.2157D,VazquezSemadeni:2007cj}.
There is another set of mechanisms in which small, dense cloudlets form by one of the possibilities listed above and then grow, by cloudlet-cloudlet agglomeration \citep{Roberts:1987eb,Dobbs:2008ez,Tasker:2009gc} or diffuse interstellar medium (ISM) accretion \citep{Goldbaum:2011kj,Heitsch:2013jp}, into GMCs. It is very hard to distinguish observationally between these mechanisms using the standard observable of molecular clouds, emission maps from molecular tracers like carbon monoxide (CO) because the position of distance to gas along the line of sight can not be measured independently of its velocity. A method that could independently measure both the distance and the velocity of these flows could make significant progress in discriminating between these many theories of converging flows.
The importance of gas flows is not restricted to the areas immediately surrounding GMCs. Gas flows on much larger scales, driven by spiral arms, may dramatically effect the pattern of galaxy-wide star formation \cite{Roberts_1972,Bonnell_2006}. Spiral shocks have been seen in strongly tidally interacting two-arm spiral galaxies in the nearby universe \cite{Visser:1980vc,Visser:1980ud,Shetty_2007} using CO and neutral hydrogen (HI) observations, and have been implicated as a critical aspect of the the star formation process. Unfortunately, the resolution in HI required to test these theories, beyond very nearby galaxies with extreme two-armed spiral structure, is not yet observationally feasible \cite{Visser:1980ud}. Thus on larger scales we also need new methods to answer critical questions about how galaxies form molecular clouds and stars.
How can we construct a method to measure this critical quantity? Since transvere velocities have been historically derived from proper motions, which cannot be reliably measured for any phase of the ISM, true 3D velocity fields are likely to remain inaccessible.
For many of these open problems we have outlined, a measurement of the line-of-sight (radial) velocity as a function of 3D position, the 1-velocity field, would represent a huge step forward. With these measurements the inflow and outflow from structures can be measured statistically, and in the context of established models, if not fully empirically.
No single observable can simultaneously measure the distance, density, and radial velocity to each parcel of gas in the interstellar medium, a construct we call the position-position-distance-velocity (PPDV) 4-cube. Instead, we have a number of observables of the ISM, each of which can be considered a (possibly lossy) projection of this PPDV cube. We dub the process of reconstructing this PPDV 4-cube from some set of observations (projections) ``Kinetic Tomography''. There are a number of existing kinds of data we can use. The first are classical radio observations of the ISM in both diffuse tracers (HI) and denser tracers (CO). These are position-position-velocity (PPV) cubes, a specific projection of the PPDV 4-cube, and have been classically used to infer the density of the ISM in 3-space assuming some kind of Galactic rotation curve \citep[e.g.][]{Levine_2006}. Another observable is the position-position-distance (PPD) reddening cubes generated by examining the photometry of large numbers of stars and performing inference on the intervening dusty ISM. There has been dramatic progress in this field \cite{Marshall_2006, Lallement_2014, Green_2015}, which is key to the present investigation. Interstellar absorption lines toward stars also represent a projection of the PPDV 4-cube, and can simultaneously contain distance information, column density information, and velocity information about the intervening gas \cite{Welsh10,Zasowski_2014,2015MmSAI..86..521Z}.
There have been some attempts to construct if not PPDV 4-cubes then maps or point estimates of $\vlos$ as a function of distance.
Most such attempts have focused on individual spiral arms and used models of ISM flows around the spiral arms to try to directly invert PPV data (e.g \citealt{1972A&A....16..118S,Foster_2006}).
\citet{Reid_2016} have developed a very different approach, which combines probability distributions coming from the standard kinematic distance, various geometric hints, and possible associations of emitting gas with structures that have parallax-based distance measurements into a single sort of syncretic probability distribution for the distance.
Neither of these approaches uses the information available in reddening-based PPD cubes.
Conversely, the method we use does not incorporate the sort of information these techniques are based on.
In the method we describe in this work, we use large-area CO and HI PPV cubes and the \citet{Green_2015} PPD cube to reconstruct the ISM PPDV 4-cube. We perform a restricted version of the full tomographic reconstruction in which we assume each parcel of gas in the PPD cube is assigned to a single central radial velocity with some radial velocity width. This can be considered an inversion of the usual kinematic distance method, in which a radial velocity is converted to a distance using a Galactic rotation curve. Here we convert each distance along a line of sight to a velocity, but instead of forcing a given rotation curve, we allow deviations from a rotation curve to best fit the PPV data cube.
In this work, we describe our method, the map it produces, and our evaluation of this map's accuracy and precision.
In Section \ref{sec:data}, we describe the datasets we use to make and evaluate our PPPV map.
In Section \ref{sec:methods}, we give a detailed explanation of the PPPV mapping technique and a brief overview of how to fit a rotation curve to a PPPV map.
In Section \ref{sec:results}, we present the PPPV map and its corresponding rotation curve.
In Section \ref{sec:validation}, we compare the PPPV map to earlier work and discuss the probable effects of known systematics.
In Section \ref{sec:conclusion}, we conclude.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/gas-compare3/gas-compare3}
\caption{{The interstellar medium in the Galactic plane, shown in three tracers: reddening (top), HI (middle), and CO (bottom). Each shows the total column measured by each tracer, with a log intensity stretch. The regions are 360 degrees in Galactic longitude wide, centered on Galactic Center, and 60 degrees in Galactic latitude tall, centered on the Galactic plane.%
}}
\end{center}
\end{figure}
\section{Data}
\label{sec:data}
\subsection{HI and CO data}
Radio emission lines in HI and CO trace the two dominant constituents of the Galactic ISM, atomic and molecular gas. Ionized phases of the ISM do not contribute significantly to the column density and will therefore not trace the extinction measured in \cite{Green_2015}. 21-cm line emission from the hyperfine transition of HI is usually optically thin and its integral is an excellent tracer of HI column:
\begin{equation}\label{XHI}
N_{HI} = 1.8 \times 10^{18} \rm cm^{-2} \frac{ \int T_B dv}{\rm K~km/s}
\end{equation}
The 21-cm line can become optically thick, meaning equation \ref{XHI} underestimates the HI column, but this mostly happens in H$_2$ dominated regions \cite{Goldsmith_2007}.
We trace molecular gas here using the 115 GHz 1-0 rotational transition of CO. The integral of this emission line can be converted to total hydrogen column CITE using
\begin{equation}\label{XCO}
X_{CO} = .
\end{equation}
This conversion factor has number of known weaknesses, stemming from much more complex excitation and opacity effects, as well as real variation in relative population of CO and $H_2$ molecules. We will address the impacts of these weaknesses in \S \ref{sec:validation}.
For our CO data, we simply use the interpolated whole Galaxy cube provided by \cite{Dame_2001}. We post-process these data with a plus-shaped median smoothing kernel, to eliminate single-pixel artifacts in the data. This filtering procedure generates less than X\% effect in the total CO emission contained in the map. In this text we investigate the Galactic plane, and so we use the sky area provided in this data set, $-30^\circ < b < 30^\circ$, and over the full range of $l$ as our grid. These data have a radial velocity resolution of 1.3 km s$^{-1}$, over -320 km s$^{-1} < V_{LSR} <$ 320 km s$^{-1}$. We find that the native resolution and PPV extent of these data are appropriate for our investigation, and thus retain the exact pixelization of these data for our analysis.
For our HI data we use a combination of three large-area Galactic HI surveys. South of declination 0$^\circ$ we use data from the 16$^\prime$ resolution GASS survey \cite{Kalberla_2010}, from declination 0$^\circ$ to 38$^\circ$ we use unpublished data from the 4$^\prime$ resolution GALFA-HI survey \cite{Peek_2011} Data Release 2, North of 38$^\circ$ we default to the 36$^\prime$ resolution LAB Survey \cite{Kalberla_2005}. We regrid these data onto the $7.5^\prime X 1.3$km s$^(-1}$ pixels of the \cite{Dame_2001} CO map.
These two maps are then combined using the above equations \ref{XHI} and \ref{XHI} to make a single N_H data cube.
\subsection{Dust data}
Our extinction cube is derived entirely from the \cite{Green_2015} reddening data. \citet{Green_2015} use PanSTARRS photometry of 800 million stars to infer the cumulative reddening along the line of sight in 6.8$^\prime$ (NSIDE=512 HEALPix) pixels. The distance information provided by the \cite{Green_2015} map is in steps of half a distance modulus, from 63 pc to 63 kpc. We regrid these data onto the \cite{Dame_2001} $l-b$ grid, and difference them in distance to find the reddening between each distance modulus. This is then converted to a $N_H$ using the factor measured in \cite{Peek_2013},
\begin{equation}
N_H = E\left(B-V\right) 7 \times 10^{21}.
\end{equation}
\subsection{High-mass star forming region data}
\label{sec:data-HMSFR}
To check the accuracy of our method, we need measurements of the line-of-sight velocities and distances of clouds of gas.
CITET REID have measured the line-of-sight velocities, proper motions, and trigonometric parallaxes of water and methanol masers associated with 103 high-mass star forming regions (HMSFRs). Of these 103, 4 fall outside of region of the Milky Way we consider in our solution. We adopt the line-of-sight velocities, velocity uncertainties, and parallaxes of the remaining 99 HMSFRs as stated in Table 1 of CITET REID. Following the discussion (but not usage) in CITET BOVY, we assume that the parallax uncertainties correspond to distance uncertainties of approximately 10\%.
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/KT-method-diagram-paper3/KT-method-diagram-paper3}
\caption{\label{fig:diagram}} A schematic diagram of the inverse kinematic distance method.%
}
\end{center}
\end{figure}
\section{Kinetic tomography}
\label{sec:KT}
We have developed a procedure for deriving the distribution of interstellar matter in position-position-distance-velocity (PPDV) space from measurements of its distribution in position-position-distance (PPD) and position-position-velocity (PPV) space.
The technical term for deriving a multi-dimensional distribution from lower-dimensional measurements is tomographic reconstruction; in this context, lower-dimensional measurements are called projections.
As should be intuitively clear, tomographic reconstruction from two projections is, in general, not possible -- there are many more independent variables than observational constraints.
Our procedure for solving our specific case of the tomographic reconstruction problem is based around exploiting some simplifying assumptions about the structure of the ISM in PPDV space to reduce the number of independent variables.
For motivation, we first examine the structural assumptions behind a common technique for finding solutions to the even more underdetermined problem of reconstructing a PPDV-space ISM distribution from PPV-space measurements alone.
The technique we are referring to here is the widely known kinematic distance method (e.g. \citealt{Levine_2006}), which assigns distances to observed line-of-sight velocities based on an assumed Galactic rotation curve and geometry.
While the kinematic distance technique is usually presented as a way of converting PPV-space distributions to PPD-space, rather than PPDV-space, distributions, its underlying assumptions make the two equivalent.
These assumptions can be combined into a single statement: (1) a location in PPD space can be assigned a single line-of-sight velocity (2) according to an assumed rotation curve and Galactic geometry.
From these assumptions it follows that knowing the PPD-space distribution of the ISM is equivalent to knowing its PPDV-space distribution.
We have PPD-space measurements in addition to PPV-space measurements, allowing us to make a more relaxed version of this assumption.
Our version is that (1) a parcel of interstellar matter in PPD space can be assigned a Gaussian distribution of line-of-sight velocities (2) whose center is within a fixed range of the line-of-sight velocity predicted by an assumed rotation curve and Galactic geometry.
Here, a "parcel" of interstellar matter refers to the contents of a single PPD voxel.
For each voxel in observed PPD cube, we aim to assign a central line-of-sight velocity $\vlos$ and a line-of-sight Gaussian velocity width $\sigma_v$.
With these assumptions, a description of the ISM in PPDV space consists of a description of its PPD-space distribution (the observed PPD cube), a line-of-sight central velocity cube in PPD space ($\vlos(\glon, \glat, d)$), and a line-of-sight velocity width cube in PPD space ($\sigma_v (\glon, \glat, d)$).
Thus we have reduced our original problem of finding a PPDV cube which is consistent with our observed PPD and PPV cubes to finding a $\vlos(\glon, \glat, d)$ and $\sigma_v(\glon, \glat, d)$ pair consistent with our PPV observations.
%We have developed a procedure for deriving the distribution of interstellar matter in position-position-distance-velocity (PPDV) space from measurements of its distribution in position-position-distance (PPD) and position-position-velocity (PPV) space. The technical term for deriving a multi-dimensional distribution from lower-dimensional measurements is tomographic reconstruction; in this context, lower-dimensional measurements are called projections. As should be intuitively clear, tomographic reconstruction from two projections is, in general, not possible -- there are many more independent variables than observational constraints. The procedure we have devised to solve our specific case of the tomographic reconstruction problem is based around exploiting some simplifying assumptions about the structure of the ISM in PPDV space to bring the number of independent variables closer to the number of observational constraints.
%Our strategy for solving this problem is based on two major assumptions: (1) a parcel of interstellar matter can be assigned a single Gaussian distribution of line-of-sight velocities and (2) that the center of this distribution resides in a bounded subregion of the PPDV space. This first assumption is a relaxed version of a core assumption of the kinematic distance method --- that a distance can be converted into a velocity. In the kinematic distance method, a direct dependence of distance on observed line-of-sight velocity is derived from a assumed rotation curve \citep[e.g.][]{Levine_2006}. In our method, we instead assume that a parcel of gas has a single \emph{central} velocity, but allow for a dispersion about this center. This assumption inherently limits our method to resolving ISM flows larger than the PPD voxels. Our second assumption can also be considered a relaxation of the kinematic distance method's assumed fixed rotation curve. We assume that the line-of-sight velocity of a parcel must lie within some interval about the value expected from an assumed rotation curve.
%In our method, a parcel of interstellar material is associated to a single PPD voxel. For each of these voxels we aim to assign both a central line-of-sight velocity $\vlos$ and a line-of-sight Gaussian velocity width $\sigma_v$. Thus, we have reduced the description of the ISM in PPDV space to the description of a density field (the original PPD cube), a line-of-sight central velocity cube in PPD space($\vlos\left(\glon, \glat, d \right)$), and a line-of-sight velocity width cube in PPD space ($\sigma_v \left(\glon, \glat, d \right)$). Thus we have reduced our original problem of finding a PPDV cube which is consistent with our observed PPD and PPV cubes to finding a $\vlos\left(\glon, \glat, d \right)$ and $\sigma_v \left(\glon, \glat, d \right)$ pair consistent with our PPV observations.
%This method can be understood as a variation on CfR's extension of kinematic distances.
%The kinematic distance to a point at some on-sky position ($\ell$, b) and line-of-sight velocity $v$ is the distance at which a model of gas motion in that ($\ell$, b) direction is equal to $v$.
%The gas motion model is usually taken to be a galactic rotation curve.
%Kinematic distances have several well-known drawbacks.
%Firstly, they are not unique in directions where different distances are assigned the same line-of-sight velocity.
%When a flat rotation curve is assumed, sightlines inside the solar circle and towards the galactic anti-center do not have unique distances.
%Next, kinematic distances ignore the possibility of deviations from the gas motion model.
%Measurements of the distances and line-of-sight velocities of point-like objects such as masers show that deviations from azimuthally symmetric rotation, or streaming motions, are the norm rather than the exception.
%In some cases, streaming motions can mean that there is no distance corresponding to a given velocity.
%This is known to happen towards the galactic anti-center CITE and near the galactic center CITE.
%Finally, computing kinematic distances for an entire PPV dataset yields a PPP dataset, trading one map of three-dimensional structure for another.
%Despite these drawbacks, kinematic distances are an adequate estimate of actual distances over much of the galactic plane.
%CfR's method improves on kinematic distances by treating them as one of a number of context clues to the distance to an object.
%Other context clues include the physical size and height above the galactic plane that a distance would imply given the object's angular size and height above the galactic plane, the presence or absence of H I 21cm self-absorption in the direction of the object, and the distances to any other objects that are known to be close to the object of interest.
%If there are enough context clues, this approach can help resolve the degeneracy between distances that correspond to the same line-of-sight velocity and can be used to detect and measure streaming motions.
%If we think of kinematic distances as a somewhat unstable zeroth-order approximation to the distance corresponding to a point in a PPV dataset, CfR's method is a second-order correction with improved stability properties.
%From this perspective, building on CfR's method is unlikely to lead to significant improvements to their solution.
%Instead, we apply a rudimentary form of their approach to the complementary problem of assigning a velocity to a point in a PPP dataset.
%With the advent of large optical and near-infrared photometric surveys, it has become possible to build maps of the three-dimensional distribution of dust out several kiloparsecs from the Sun CITE.
%If we make some assumptions about the mixing and relative abundances of dust and gas in the ISM, we can consider these maps to be of the three-dimensional stucture of the ISM.
%To get a zeroth-order velocity for a point in one of these PPP maps, we can take inspiration from kinematic distances and evaluate the galactic rotation curve at its position.
%Even at this stage, the PPP to PPPV problem has an advantage over the PPV to PPPV problem -- the position to velocity assignment is unique, if not necessarily accurate.
%There are different possible approaches to improving this approximation.
%Our approach is to use a PPV dataset of the ISM over the same area as a guide.
%Reproducing this PPV dataset requires modestly adjusting the velocities assigned to the points in the PPP dataset.
%These adjusted velocity assignments are our PPPV map.
%We have checked this map against the maser measurements in CITET REID (R14).
%Reid et al. have measured the distance to and line-of-sight velocity of a number of masers.
%Our PPPV map evaluated at the on-sky positions and distances of the masers agrees with the masers' measured line-of-sight velocities.
%This is particularly heartening considering that many of the masers are known to deviate significantly from galactic rotation.
\subsection{Formalism}
\label{sec:KT-method}
In the introduction, we described our method for mapping the ISM in PPPV space as a kind of inversion of the usual idea of a kinematic distance. Instead of setting the distance of a voxel in a PPV cube based on the voxel's line-of-sight velocity and a rotation curve, set the velocity of a voxel in a PPP cube based on the voxel's distance and a rotation curve. Since a voxel is not a point, we need to assign it a distribution of velocites rather than a single velocity; we use a Gaussian. This operation is shown for all of the voxels along a single line of sight in Figure (first proj). If we apply this operation a large number of sightlines at a single galactic latitude, we get the $\ell$-$v$ map shown in figure (first lv).
The expected velocity-based $\ell$-$v$ map does not look like the observed $\ell$-$v$ map. Some of the differences between the two $\ell$-$v$ maps are caused by the low distance resolution of the PPP cube beyond a few kiloparsecs, while others come from us assigning the same fiducial width to every PPP voxel's velocity Gaussian. Many of the differences, however, are clearly caused by the expected velocities being incorrect. These velocity offsets affect not just single pointings but entire extended structures. For example, the line-of-sight velocity of the cloud at $\ell \approx 155^\circ$ is 0 km/sec in the flat rotation-derived map and 5 km/sec in the observed map. In this and many other cases, the necessary adjustments to the expected velocities can be picked out by eye.
Since deviations from the expected velocities are clearly necessary, we need to allow them. Practically, "allowing deviations" means varying the center, $\mu_k$, and width, $\sigma_k$, of each distance voxel $k$ to minimize the difference between the model and observed $\ell$-$v$ maps. Along a single line of sight, the difference is given by the expression
\begin{equation}
\sum_u ((\sum_k f(v_u; \mu_k, \sigma_k)) - y_k)^2,
\end{equation}
where $v_u$ and $y_u$ are the $u$th velocity and column density values along the velocity axis and $f$ is the pixel-convolved Gaussian function.
The objective function over all lines of sight is
\begin{equation}
\label{eqn:objective-nosprings}
\sum_i \sum_j \sum_u ((\sum_k f(v_u; \mu_{i,j,k}, \sigma_{i,j,k})) - y_{i,j,k})^2,
\end{equation}
where $i$ and $j$ are indices along the $\ell$ and $b$ directions.
Based on visual inspection of the flat rotation-derived and observed $\ell$-$v$ maps and physical arguments about plausible magnitudes for streaming motions (CITE), we place a bound of 40 km/sec on the magnitude of deviations from flat rotation. Explicitly, $\mu_k$ must be within 40 km/sec of the line-of-sight velocity corresponding to flat rotation at the position of distance voxel $i$. We have found that maximum deviations between 35 and 50 km/sec give similar, though not identical, results.
Figure (varproj) shows a typical [REVISE, FAKE NOT TYPICAL] solution for a single pointing, and Figure (vary lv) shows an $\ell$-$v$ map derived from an actual solution. Apart from a few small mismatches (e.g. at $\ell \approx 80^\circ$ and $v\approx$ $-40$ km/sec), this model gives a faithful reproduction of the observations.
For a given starting PPP 3-cube, there will generally be more than one solution that correctly reproduces the observed PPV 3-cube. Since we are interested in the per-PPP-voxel velocity assignments, we need a way of selecting the solution that is, in some sense, more likely to be the correct one. We do this by adding the constraint that velocity centers of PPP voxels that are adjacent in $\ell$ and $v$ should be similar. We implement this by adding the following term to the objective function in equation \ref{eqn:objective-nosprings}:
\begin{equation}
\lambda \times \sum_k \left( \sum_i (\mu_{i,j,k} - \mu_{i+1,j,k})^2 +
\sum_j (\mu_{i,j,k} - \mu_{i,j+1,k})^2 \right).
\end{equation}
Using this term is an example of Tikhonov regularization. As in equation \ref{eqn:objective-nosprings}, $i$, $j$, and $k$ are indices for the $\ell$, $b$, and distance axes. The strength of the regularization, i.e. the degree to which we require neighboring voxels to have similar centers, is set by $\lambda$. The higher the value of $\lambda$, the more regularized the solution.
In practice, regularizing the solution mainly affects the central velocities of relatively faint, low-mass PPP voxels. Along a typical sightline, there will be few high-mass voxels in both distance and velocity. The velocity of these high-mass voxels will be well-constrained and will not change from reprojection to reprojection. Moving one of these voxels to a different central velocity would greatly increase the difference bettween the model and observation. Conversely, low-mass voxels tend to be more interchangeable from reprojection to reprojection and more locally movable within a reprojection. As a result, including a regularization term favors solutions in which low-mass voxels' central velocities are similar to any near-by high-mass voxels' central velocities.
Figure (spring proj) shows the effect of including this constraint on a single line of sight solution, and figure (spring lv) shows the corresponding lv map. In both cases, the \emph{fit}, i.e. similarity of model and observation, is not as good as the fit of the unconstrained solution. We find the trade-off between slightly biased central velocities throughout the diffuse ISM on a local scale and a more unique and correct solution on a global scale to be acceptable. We also find that the regularized solution is in substantially better agreement with the reference dataset, which has complete 6-dimensional phase space measurements for a small number of discrete points throughout the galaxy (see section \ref{sec:KT-validation}).
The actual PPPV reconstruction, then, comes from evaluating the velocity Gaussians from the previous step at their positions in the PPP cube and not projecting down to PPV space. A more easily visualizable summary of the solution is the 3-cube of velocity centers, which gives the mean motion of the ISM as a function of on-sky position and distance.
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/cl-gd-smooth-maser-comp/maser-comparison}
\caption{\label{fig:maser_comparison}} All three panels show the residuals from flat rotation of the \protect\citet{2009ApJ...700..137R} HMSFR line-of-sight velocity observations on the y-axis and residuals from flat rotation of the values of the \protect\citet{1985ApJ...295..422C} rotation curve (left panel), the unregularized kinetic tomography solution (middle panel), and the regularized kinetic tomography solution (right panel) at the positions of the HMSFRs on the x-axis.%
}
\end{center}
\end{figure}
\subsection{Checking the KT solution}
\label{sec:KT-validation}
The CITET REID HMSFRs (see Section DATA.3) are embedded in and, presumably, moving with dense molecular gas, and can therefore serve as point-probes of the ISM's velocity field.
By comparing the HMSFRs' observed line-of-sight velocities to those a given technique assigns to their locations in PPD space, we can determine the accuracy of that technique.
By comparing the accuracies of the regularized and unregularized KT solutions, we can get a sense of the effectiveness of our regularization scheme.
By comparing the accuracy of both KT solutions to that of the CITE CLEMESN non-flat rotation curve, we can get a sense of how much accuracy moving beyond axisymmetry can add to our description of the ISM.
While the $\glon$ and $\glat$ of an HMSFR from CITET REID are known to effectively arbitrary precision, its distance is only known to within about 10\%.
To take this distance uncertainty into account when assigning a line-of-sight velocity to an HMSFR based on a KT solution, we take an average, weighted by the HMSFRs distance probability density function, of the KT solution's $\vlos(d)$ profile towards the HMSFR's $\glon$ and $\glat$.
A similar procedure gives us an estimate of the uncertainty, expressed as a standard deviation, of this $\vlos$ value.
When computing the line-of-sight velocity of an HMSFR according to a rotation curve, we ignore this distance uncertainty and simply use the best-fit distance.
We show a comparison of the HMSFRs' measured line-of-sight velocities and line-of-sight velocities assigned based on the CITET CLEMENS rotation curve, unregularized KT, and regularized KT in Figure (maser comp).
There is a very clear improvement from a rotation-curve to either version of KT and a slight but noticeable improvement from unregularized to regularized KT.
To get a crude quantitative estimate of the effect of regularization, we can compute the reduced $\chi$-squared values of the two sets of velocity estimates.
If we assume the regularization parameter counts against the number of degrees of freedom, the reduced $\chi$-squared values of the unregularized and regularized KT solutions are 5 and 3, respectively.
Much of the total discrepancy between the regularized KT solution and the HMSFR measurements is driven by 5 catastrpphic outliers.
If we remove these 5 outliers, leaving 94 HMSFRs, the reduced $\chi$-squared values of the unregularized and regularized KT solutions drop to 4 and 1.3, respectively.
We consider the advantage of regularzied over unregularized KT to be sufficient to adopt the regularized KT solution as \emph{the} KT solution, and will refer to it as such below.
At the positions of 94 of 99 HMSFRs, at least, KT performs remarkably well.
Next, we discuss the 5 catastrophic outliers.
%We have used the unregularized and regularized versions of the kinetic tomography technique to derive PPPV 4-cubes for the region of the galactic plane with $0^\circ < \ell < 220^\circ$ and $-10^\circ < b < 10^\circ$. %We will refer to these two solutions by the names of their methods, i.e. as the UKT and RKT solutions.
%These solutions can be thought of as either one 4-cube of the ISM density in PPPV space or as a set of three 3-cubes of ISM density, central velocity ($v_{los}$), and velocity dispersion in PPP space. One benefit of the ``three 3-cubes'' representation is that the accuracy of each of the three 3-cubes can be checked separately. The accuracy of the 3-cube of ISM density map is discussed in Schlafly+. In this section, we check the accuracy of the unregularized and regularized $v_{los}$ 3-cubes against the R14 dataset, which is described in section \ref{sec:data-HMSFR}. We find that both $v_{los}$ 3-cubes represent clear improvements over flat and \cite{Clemens_1985} rotation curves, and that the regularized solution is statistically consistent with the R14 measurements, with the exception of a handful of outliers.
%The quantities we are comparing are the line-of-sight velocity of high-mass star forming regions (HMSFR) from R14 with the $v_{los}$ 3-cube evaluated at the $\ell$, $b$, and distance of the HMSFRs. To compute the central velocity from the $v_{los}$ 3-cube, we take the $v_{los}$ profile of the $\ell$, $b$ cell that contains a HMSFR, linearly interpolate $v_{los}$ between the centers of the distance bins, and evaluate this piecewise linear $v_{los}$ profile at the distance of the HMSFR given by R14. We estimate the uncertainty of this velocity by propagating the distance uncertainty of the HMSFR to a $v_{los}$ uncertainty. This distance-derived uncertainty is added in quadrature with the radial velocity uncertainty quoted for the HMSFR in R14. The typical magnitude of this total uncertainty is 3-10 km/sec; Figure (HMSFR on profile) shows an HMSFR for which this uncertainty is particularly large. We also compute the $v_{los}$ predicted for each HMSFR by the \citet{Clemens_1985} rotation curve using the same procedure; we refer to these $v_{los}$ values as the \Clemens values.
%Figure (maser comp) compares the R14 measured HMSFR $v_{los}$ to the \Clemens, unregularized, and regularized $v_{los}$ values. There is a clear improvement in agreement from the \Clemnens to the unregularized, and from the unregularized to the regularized $v_{los}$. The unregularized $v_{los}$ captures the full amplitude of the discrepancy between the R14 measurements and a flat rotation curve while \Clemens rotation curve does not. As the regularized solution is clearly superior to the unregularized solution, we devote the rest of the text to discussing the regularized solution, and consider it our best measurement of the $v_{los}$ 3-cube.
%For a more quantitative assessment, we can analyze the distribution of \emph{standardized residuals} -- residuals divided by uncertainties. This distribution is shown in Panel (a) of Figure (5). For comparison, the figure also shows a normal distribution with mean zero and standard deviation one, a \emph{standard normal} distribution. If the uncertainties are accurate and errors are normally distributed, the standardized residuals should be a drawn from the standard normal distribution. Of the 99 HMSFRs that fall within our $\ell$ and $b$ range, 5 have residuals that are more than three standard deviations away from zero; we discuss these six in section \ref{sec:discussion-catastrophic}. The standardized residuals of the remaining 94 objects are consistent with being drawn from the standard normal distribution. We discuss what this means for the more global accuracy of the regularized solution in section \ref{sec:discussion-systematics}.
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/maser-pie/maser-pie-ann}
\caption{\label{fig:maser_pie}}
Colors show the difference between an observed or estimated $\vlos$ and the value predicted from a flat rotation curve. The background is an average of the PPD velocity field over $-2.5\deg \leq \glat \leq +2.5\deg$. The color of the inner part of each circle is the $\vlos$ of an HMSFR. The color of the outer ring of each circle is the value of the PPD velocity field at the position (including $\glat$) of the HMSFR.%
}
\end{center}
\end{figure}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/HMSFR-outliers/HMSFR-outliers}
\caption{\label{fig:outliers}}
\textbf{(a)} The distribution of residuals per uncertainty (\textit{standardized} residuals) between the observed and KT-derived line-of-sight velocities of HMSFRs (filled histogram), with an appropriately scaled standard normal distribution for reference (empty histogram). We consider HMSFRs with $\vert \Delta v / \sigma \vert > 3$ to be outliers. In all three panels, outliers are denoted in orange and non-outliers are denoted in blue.
\textbf{(b)} The distribution of the HMSFRs' Galactic longitudes. HMSFRs (outliers and non-outliers) with $\glon < 35^\circ$ are outlined with a dashed black line.
\textbf{(c)} The distribution of the HMSFRs' Galactocentric radii. HMSFRs (outliers and non-outliers) with $\glon < 35^\circ$ are outlined with a dashed black line.
\textbf{(d)} The locations of the outliers and non-outliers in heliocentric Cartesian coordinates, where $\mathrm{X_{HC}}$ increases towards $\glon=0^\circ$ and $\mathrm{Y_{HC}}$ increases towards $\glon=90^\circ$.
\textbf{(e)} The probability of finding m or more of the n=5 total outliers among M=31 HMSFRs drawn at random and without replacement from the full population of N=99. The $\glon < 35^\circ$ subsample contains 31 HMSFRs, total, and 5 outliers. The probability of this happening if an HMSFR's chance of being an outlier is independent of Galactic longitude is well below the $2\sigma$-equivalent threshold of $5\%$.
\textbf{(f)} The same probabilities as in (e), but for N=99, n=5, M=17 (solid lines) and N=31, n=5, M=17 (dashed lines). In both cases, the number M of HMSFRs in the subsample is the number of HMSFRs with $3 \text{ kpc} < \mathrm{R_{GC}} < 5\text{ kpc}$. The N=99 and N=31 cases correspond to the subsample being drawn from the full 99 HMSFR sample and from the $\glon < 35^\circ$ subsample. The N=99 case is clearly inconsistent with outliers being independent of $\mathrm{R_{GC}}$. The probability of the N=31 case is just below the 2$\sigma$ equivalent level, suggesting that the outliers' $\mathrm{R_{GC}}$ values stand out from those of all HMSFRs with $\glon < 35^\circ$.%
}
\end{center}
\end{figure}
\subsection{Catastrophic outliers}
\label{sec:discussion-catastrophic}
In Section \ref{sec:KT-validation}, we showed that the measured and KT-derived $\vlos$ values of 94 of the 99 \Reid HMSFRs are in agreement.
In this section, we discuss the five HMSFRs where they do not agree.
Panel (a) of Figure \ref{fig:outliers} shows the distribution of standardized residuals of the full 99 HMSFR sample along with a scaled standard normal distribution.
There are five HMSFRs, which have been indicated on the figure in orange, that are clearly separated from the rest.
These are the HMSFRs for which we consider the measured and KT-derived $\vlos$ values to disagree.
We will refer to these five as the outliers.
We have not been able to find an obvious reason for why these five HMSFRs, specifically, are outliers.
None of their properties are extreme -- they are not all clumped together; they are not particularly close to or far from the Sun; the total dust and gas columns along their sightlines are typical for HMSFRs in their vicinity; they are reasonably close to their nearest neighbors, none of which are outliers; the differences between their measured $\vlos$ values and those of their close neighbors are reasonable, as are the differences between the magnitudes of their velocity three-vectors; they do not all sit in the same spiral arm (\Reid).
There is no clear distinction between an outlier and its non-outlier neighbors.
While there does not appear to be a way to predict whether a specific HMSFR will definitely be an outlier, we have found that HMSFRs with $0^\circ < \glon < 35^\circ$ and Galactocentric radii ($\RGC$) between 3 and 5 kpc are more likely to be outliers.
The distributions of outlier and non-outlier $\glon$ and $\RGC$ are shown in panels (b) and (c) of Figure \ref{fig:outliers}, and the locations of the outlier HMSFRs and some nearby non-outliers are shown in panel (d).
By visual inspection, the outliers' $\glon$ and $\RGC$ distributions clearly deviate from those of the full sample.
We can quantify this deviation by assuming that any HMSFR is equally likely to be one of the five outliers and computing the probability of finding all five in a randomly chosen subsample of the same sizes as one of our two selections.
Since a single HMSFR can not be present in a sample twice, the appropriate distribution to use for this calculation is the hypergeometric distribution, which assumes that the subsample is chosen without replacement.
The probabilities of finding all five outliers in a random subsample the size of the $\glon$ and $\RGC$ selections are both well below 1\%.
Since every HMSFR in the $\RGC$ selection is also in the $\glon$ selection, this calculation does not tell us which selection is responsible for the enhance outlier probability.
To try to differentiate between the two, we can repeat the above exercise with the $\RGC$ selection as the subsample and the $\glon$ selection, rather than the full set of 99 HMSFRs, as the population.
The probability of finding all five outliers in a random 17-element subsample of the 31-element $\glon$ selection is about 4\%.
This probability is low, certainly low enough to suggest that there may be an enhancement in the outlier rate in the $\RGC$ selection over the $\glon$ selection's rate, but not enough for that to be a strong conclusion.
This conclusion is weak particularly considering that we have neglected the distance and outlier/non-outlier classification uncertainties.
We have attempted to distinguish between the two selections because they suggest different causal explanations for the increased likelihood of being an outlier.
If being in the $\glon$ selection is responsible for most of the effect, the likely cause is increased line-of-sight confusion.
The closer to $\glon$ of $0^\circ$ a sightline passes, the more matter it typically intersects near the Galactic center and, assuming the Galactic disk is roughly circular, the more matter it passes through on the opposite side of the Galaxy.
If instead being in the $\RGC$ selection is what matters, the cause is likely to be the complexity of ISM dynamics near the Galactic center.
The distance resolution of the PPD cube declines with distance, so it would not necessarily be surprising if we were unable to spatially resolve flows with small spatial scales.
The end of the Galactic bar is generally considered to fall right in the middle of the five outliers \cite{2016arXiv160207702B}.
If we treat sightlines with HMFSRs as typical, the implied catastrophic failure rate $\RGC <$ 5 kpc is XX\%.
Conversely, there are no outliers with $\RGC >$ 5 kpc or $\glon >$ 35$^\circ$.
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/figure-/KT-and-Clemens-rcurves}
\caption{\label{fig:rotation_curves}}
A piecewise-flat rotation curve derived from our kinetic tomography solution (black) and a rotation curve from \protect\citet{Clemens:1985dp} (red).%
}
\end{center}
\end{figure}
\section{The rotation curve of the Milky Way and spatially coherent streaming motions}
\label{sec:rotation_curve}
In this section, we analyze the KT-derived $\vlos$ field in the Galactic plane.
The primary purpose of this analysis is to establish that on a coarse scale, we basically agree with the Milky Way kinematics literature.
This serves as a complement to the previous section, where we established that we agree with fine-grained measurements of dense HMSFRs.
We first decompose the in-plane $\vlos$ field into a rotation curve and residual velocity field under different assumptions about the form of the rotation curve and compare the results to results from the Milky Way kinematics literature.
We find that the magnitude and spatial scale of spatially coherent streaming motions is such that we cannot reliably fit a rotation curve without a model for these streaming motions.
When we fit a non-flat rotation curve to the in-plane $\vlos$ field, these streaming motions are incorporated into the rotation curve.
This non-flat rotation curve is consistent with the \citet{Clemens} rotation curve, which also describes large-scale streaming motions.
When we instead attempt to fit a flat rotation curve while assuming residuals from flat rotation are spatially uncorrelated, we get results that change dramatically depending on which subset of the $\vlos$ field we use and which are inconsistent with estimates of the Galactic rotation rate from the literature.
Next, we discuss the residual velocity field.
We compare it to the two-dimensional residual velocity field of \citet{BrandBlitz} and find that the two agree.
We associate some of the clear large-scale streaming motions with features discussed in \citet{Clemens}, \citet{McClureGriffithsEarly}, and \citet{McClureGriffithsLate}.
Finally, we give a rough estimate of the spatial scale of the strongest fluctuations and compare this to power spectra of velocity fluctuations derived from measurements of gas \citep{Clemens,McClureGriffithsLate} and stellar kinematics \citep{Bovy2015}.
\subsection{The rotation curve}
\label{sec:rotation_fit}
We assume that the global motion of the Milky Way ISM can be entirely described by a rotation curve, $V(\RGC)$.
In particular, we assume that there are no global radial flows and that azimuthal and radial departures from the rotation curve average to zero over the Galaxy.
We allow for local, spatially independent departures from the rotation curve and model these departures as a Gaussian perturbation superimposed on the global motion.
We assume that the azimuthal and radial components of the velocity perturbation are uncorrelated and have standard deviations $\sigma_\phi$ and $\sigma_R$.
In addition to the motion of the ISM, we need to describe the location and motion of our observing site, the Sun.
We use $\RGC_{, \odot}$ for the Sun's Galactocentric radius and $V_{\phi, \odot}$ and $V_{R, \odot}$ for its azimuthal and radial motion relative to the Galactic center.
Because we will only be using the $\glat \approx 0 \deg$ portion of the $\vlos(\glon, \glat, d)$ cube, we can ignore the Sun's motion perpendicular to the Galactic plane.
To compute the likelihood of a parameter set $V_\phi (\RGC)$, $\sigma_\phi$, $\sigma_R$, $R_{\odot}$, $V_{\phi, \odot}$, and $V_{R, \odot}$, we first need to convert these parameters to a model $\hat{\vlos}(\glon, d)$.
The Galactocentric radius of a given $\glon$ and $d$ at $\glat = 0 \deg$ is
\begin{equation}
\RGC = \sqrt{\RGC^2_{,\odot} + d^2 - 2 \RGC_{, \odot} d \cos \glon}.
\end{equation}
The mean line-of-sight velocity at a $\glon$ and $d$ is
\begin{equation}
\label{eqn:vlos_mean}
\hat{\vlos}(\glon, d) = \frac{V_\phi (\RGC) \RGC_{,\odot} \sin\glon}{\RGC} - V_{\phi, \odot} \sin \glon + V_{R, \odot} \cos \glon.
\end{equation}
The standard deviation of the distribution of line-of-sight velocities about this mean is
\begin{equation}
\label{eqn:vlos_var}
\sigma^2_\mathrm{los} = \sqrt{\left(\sigma_\phi \frac{\RGC_{,\odot} \sin \glon}{\RGC} \right)^2 +
\left(\sigma_R \frac{\RGC_{,\odot} \cos \glon + d (1 - 2 \cos^2 \glon)}{\RGC} \right)^2}.
\end{equation}
To these model $\hat{\vlos}(\glon, d)$ values we compare the mass-weighted $\glat$-average of the KT-derived $\vlos(\glon, \glat, d)$ cube for $-2.5 < \glat < + 2.5$; we will call this $\glat-$averaged quantity the $\vlos(\glon, 0 \deg, d)$ \emph{map}.
Our likelihood function for a single $(\glon, 0\deg, d)$ voxel is a Gaussian with mean $\hat{\vlos}(\glon, d)$ and standard deviation $\sigma_\mathrm{los}(\glon, d)$.
The full likelihood function is the product of the single-voxel likelihood functions over all $(\glon, 0\deg, d)$ voxels.
We fit two sets of models, one in which we assume a flat rotation curve and let the Solar parameters vary and one in which we assume a piecewise-linear rotation curve with breaks every 0.25 kpc and fix the solar parameters to $\RGC=8.1$ kpc CITEP E.G. BOVY, $v_{\phi, \odot}=TK$, and $v_{R, \odot}=TK$ (CITE THAT STANDARDS PAPER THAT'S CITED IN IBID).
We use different portions of the $\vlos(\glon, 0\deg, d)$ map for the different models.
There are three competing considerations for choosing what portion of the map to use --- the quality of the KT solution drops beyond some maximum heliocentric distance, model parameters become less degenerate as the $\glon$ and $\RGC$ ranges of the data widen, and parts of the Galaxy are known to deviate from flat rotation (e.g. CITEALT WRONG ROTATION CURVE).
Motivated by the third consideration, our fiducial ranges for the flat rotation curve case are TK and TK.
In the piecewise-linear rotation curve case, the first and second considerations are more important, so we use DRANGE TK and RGCRANGE TK, which translates to a non-separable $\glon$-$d$ constraint.
For the fiducial $\glon$- and $d$-ranges, our results for the flat rotation curve model are $V_\phi=TK \pm TK$, $\sigma_\phi = TK \pm TK$, $\sigma_R = TK \pm TK$, $R_\odot = TK \pm TK$, $v_{\phi, \odot} = TK \pm TK$, and $v_{R, \odot} = TK \pm TK$.
The $V_\phi$ and $R_{\odot}$ values are in tension with the IAU-recommended values, as well as the results of CITET BOVY and CITET REID, if we take our quoted uncertanties to be accurate.
If we vary the outer limit of the $d$-range used in the fit from TK to TK, $V_\phi$ varies from TK to TK and $R_\odot$ varies from TK to TK.
This variation is smooth as a function of the outer limit of the $d$-range, indicating that the variation is not driven by a small problem region.
We believe that the cause of this variation is global model misspecification --- we have assumed that deviations from the rotation curve are spatially uncorrelated when they are in fact correlated, and on spatial scales that are of order the spatial extent of the $\vlos(\glon, 0\deg, d)$ map.
We will discuss this issue further in Section \ref{sec:rotation_discussion}.
Our best-fit piecewise-linear rotation curve is shown, along with the CITET CLEMENS rotation curve for comparison, in Figure \ref{fig:rotation_curves}.
Both curves can be qualitatively described as 3 ``bumps'' or ``wiggles'' about a mean circular velocity.
Pointwise, the two curves agree to within about 10 km/sec.
In terms of the qualitative description, the curves are in even better agreement --- the mean circular velocities and the locations of the bumps in the two curves are essentially the same, while the amplitudes of the bumps are greater by the same factor of $\approx 10$ km/sec for the KT-derived rotation curve.
The differences that are present can most likely be explained by differences in the sets of $\vlos$ measurements that go in to curves' derivations.
%where $V_\phi$ and $V_R$ are azimuthal (increasing clockwise) and radial (increasing outward from the Galactic center) velocities, $R_\odot$, $V_{\phi, \odot}$, and $V_{R, \odot}$ are the Galactocentric radius, azimuthal velocity, and radial velocity of the Sun, respectively, and $\RGC$ is the Galatocentric distance to the point specified by $\glon$ and $d$:
%Here, we have assumed that $V_\phi$ and $V_R$ depend only on $\RGC$, meaning that they are constant as a function of Galactic azimuth.
%We will assume that there are no azimuthally constant radial flows, meaning that $V_R(\RGC) \equiv 0\text{ km/sec}$.
%Following, e.g., \citet{Reid:2009jb} and \citet{Bovy_2009}, we allow for non-zero azimuthal and radial velocity dispersions, parametrized as independent normal distributions.
%In this parametrization, the line-of-sight velocity dispersion at a single location is a normal distribution with variance
%\begin{equation}
%\label{eqn:vlos_var}
%\sigma^2_\mathrm{los} = \left(\sigma_\phi \frac{R_\odot \sin \glon}{\RGC} \right)^2 +
%\left(\sigma_R \frac{R_\odot \cos \glon + d (1 - 2 \cos^2 \glon)}{\RGC} \right)^2,
%\end{equation}
%where $\sigma_\phi$ and $\sigma_R$ are the standard deviations of the azimuthal and radial velocity dispersions.
%We use two models for the rotation curve and other Galactic parameters.
%In the first case, we assume a flat rotation curve, expressed as a single radially constant $V_\phi(\RGC) \equiv V_\phi$, and solve for the best-fit $V_\phi$, $V_{\phi, \odot}$, $V_{R, \odot}$, $\RGC$, $\sigma_\phi$, and $\sigma_R$.
%In the second, we assume a continuous, piecewise linear rotation curve and fix the Galactic parameters to a fiducial set of values -- $\RGC = 8.1\text{ kpc}$, $V_{\phi, \odot} - V_{\phi} = 15.3 \text{ km/sec}$, and $V_{R, \odot} = - 10.3 \text{ km/sec}$.
%In the first case, we only use $\vlos$ measurements within the nearest 4 kpc.
%This avoids the inner 4 kpc of the Galaxy, where the rotation curve is known to be non-flat.
%In the second case, we use $\vlos$ measurements that are between 2 and 15 kpc from the Galactic center on the near side of the Galaxy.
%The breakpoints of the piecewise linear curve are placed every 0.25 kpc.
%Our results for the parameters of the first model are $\RGC=8.90 \pm 0.04 \text{ kpc}$, $V_\phi = 279 \pm 1.5 \text{ km/sec}$, $V_{\phi, \odot} - V_\phi = 28.5 \pm 0.1 \text{ km/sec}$, and $V_{R, \odot} = 6.2 \pm 0.1 \text{ km/sec}$.
%While the formal uncertainties on these values are small, we can tell that the systematic uncertainties must be large, particularly on $\RGC$ and $V_\phi$, because the results depend strongly on which portion of the measurements we include in the fit.
%The parameter values we give above use the measurements from the nearest 4 kpc of the Galactic disk.
%If we instead use measurements from the nearest 3 kpc, we find that $\RGC$ changes to $8.05 \pm 0.06 \text{ kpc}$ and $V_\phi$ remains unchanged, while if we use measurements from the nearest 6 kpc and exclude the inner Galaxy by only using sightlines with $\glon > 45^\circ$, we get $\RGC = 7.43 \pm 0.03 \text{ kpc}$ and $V_\phi = 223 \pm 1 \text{ km/sec}$.
%We believe the cause of these large changes in parameter values is that as we vary the region being analyzed, we include and exclude different large-scale streaming motions.
%Since these large-scale streaming motions are not included in the model specification, they end up dragging parameters in different directions.
%We discuss this further in Section \ref{sec:rotation_discussion}.
%The piecewise-linear rotation curve is shown in Figure \ref{fig:rotation_curves}, with the \Clemens rotation curve for comparison.
%These two non-flat rotation curves mostly agree, in the sense that deviations from flat rotation have the same sign, if not always the same magnitude.
\subsection{Spatially coherent streaming motions}
\label{sec:rotation_discussion}
The Milky Way velocity field has structure on four main spatial scales --- the size of the Galaxy i.e. the rotation curve; $\approx 2$ kpc (e.g. CITEALT CLEMENS for atomic and molecular gas, )
We are looking for motions with spatial scales larger than of order 100 parsecs.
[studies such as Clemens85, B&B93, Bovy+15, McG+16 have found high-amplitude motions on scales of about 200 pc, 2 kpc, and the Galaxy. there's also cloud-cloud motions, seen in high-density tracers (Clemens85). vertical averaging may have washed some stuff out; also most voxels have a long side length > 200pc. so, we don't expect to cleanly detect the 200 pc motions or cloud-cloud motions.]
In Figure \ref{fig:six_pies}, we show residuals between the $\vlos$ measurements and (panel b) the flat rotation curve derived from the nearest 4 kpc and (panel c) the piecewise-linear rotation curve.
For comparison, we also show residuals between the $\vlos$ measurements and (panel a) our fiducial flat rotation curve and Galactic parameters, (panel d) the \Clemens non-flat rotation curve, (panel e) the \cite{2012ApJ...759..131B} flat rotation curve, (panel f) and the \Reid A5 model flat rotation curve.
In every panel, there are obvious, multi-kpc regions with consistent non-zero velocity residuals.
Most of these regions persist across all six panels, though there is some panel-to-panel variation in the exact bounds and residual magnitude of each region.
For example, there is always a blue-shifted feature just outside the solar circle at $90^\circ \lesssim \glon \lesssim 180^\circ$, a smaller red-shifted feature $\approx 12$ kpc from the Galactic center at the high-$\glon$ boundary of the map, and another red-shifted feature just inside the solar circle.
Persistence across different sets of rotation curves and Galactic parameters provides evidence that the spatially coherent residuals are not caused by a misspecification of the global Galactic motion and geometry.
We can argue for the accuracy of most of the residual regions by taking the HMSFR measurements as truth and appealing to the regions' strong spatial coherence.
The argument goes as follows -- if the $\vlos$ map accurately predicts the velocities of HMSFRs that are spatially coincident with a specific residual region, it is reasonable to assume that the $\vlos$ values assigned to the rest of that region are also likely to be accurate.
Based on the discussion in Section \ref{sec:KT-validation}, we can claim that the HMSFR velocities are accurately predicted everywhere outside of the inner 4-5 kpc of the Galaxy or so.
A visual inspection of Figure \ref{fig:maser_pie} shows that there are HMSFRs in every residual region within about 6 kpc of the Sun.
We interpret the combination of these two facts as evidence that all of the residual regions shown in Figures \ref{fig:maser_pie} and \ref{fig:six_pies} outside of the inner 4-5 kpc of the galaxy are, in fact, actually present in the velocity field of the Milky Way ISM.
Once we have established that the spatially-structured residuals are real, we can ask what they are.
In particular, we can ask what fraction of the residuals can be ascribed to radial rotation curve variations and what fraction should instead be interpreted as streaming motions.
If we compare the residuals from the two radially-varying rotation curves in Figure \ref{fig:six_pies} (panels c and d) to the residuals from the four flat rotation curves (panels a, b, e, and f), the typical magnitudes of the former are clearly lower than typical magnitudes of the latter.
While we could interpret this difference in residual magnitudes as evidence that the residuals are mostly attributable to radial rotation curve variations, that interpretation would only be true for a very empirical definition of a rotation curve.
We can consider two definitions of a rotation curve: the azimuthal average of the azimuthal component of the Galactic velocity field (empirical) and the circular velocity in the Galactic gravitational potential (theoretical), both of which may be functions of Galactocentric radius.
If azimuthal departures from the theoretical rotation curve, i.e. azimuthal streaming motions, at some Galactocentric radius do not average to zero, then the empirical and theoretical rotation curves at that radius will not be equal.
We are particularly prone to cross-talk between spatially coherent streaming motions and radial variations in the rotation curve because of the limited azimuthal coverage of our $\vlos$ measurements.
At some radii, a streaming motion that is coherent over a few kiloparsecs is enough to fill most of the azimuth range in our map at those radii.
In the inner Galaxy, particularly near the tip of the Galactic bar, this danger has been pointed out by e.g. \citet{Chemin_2015}.
Those authors argue that most of the radial variations in the inner Galaxy in \Clemens, in particular, are actually streaming motions.
While our azimuth coverage is wider than that of \Clemens, the similarity in the rotation curve shapes (Figure rotation curves) suggests that the point still applies.
While \citet{Chemin_2015} assume that deriving the rotation curve in the outer Galaxy should be relatively free of hassle, extragalactic studies such as \citet{Meidt_2013} and (some of those Halpha ones) have shown that strongly spiral galaxies can have quite strong streaming motions far away from the bar.
So, the advice seems to be to just steer clear of assuming fluctuations in azimuthal velocities have much to do with the theoretical rotation curve -- it's probably all streaming motions.
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/six-pies2/six-pies}
\caption{\label{fig:six_pies}}
Maps of the difference between an average of the $\vlos$ field over $-2.5 \leq \glat +2.5$ and six sets of Galactic parameters. The dashed arcs correspond to distances, from right to left, of 4, 8, and 12 kpc from the Galactic center. The parameter sets are: (a) a fiducial flat rotation curve and Galactic parameter set; (b) the flat rotation curve and Galactic parameter set, defined in Section \ref{sec:rotation_fit}, that minimizes the square sum of the $\vlos$ residual field within 4 kpc of the Sun; (c) the fiducial values of the Sun's Galactocentric radius and motion relative to the local standard of rest with the non-flat rotation curve derived in Section \ref{sec:rotation_curve}; (d) the non-flat rotation curve derived in \protect\citet{Clemens:1985dp} assuming $\RGC=8.5 \text{ kpc}$ and $V_\phi(8.5 \text{ kpc})= 220\text{ km/sec}$; (e) the flat rotation curve and Galactic parameters derived in \protect\citet{Bovy_2012}; (f) the "fit A5" flat rotation curve and Galactic parameters derived in \protect\citet{Reid:2014km}.%
}
\end{center}
\end{figure}
\section{Discussion}
\label{sec:discussion}
\subsection{Unmet assumptions and their potential consequences}
\label{sec:discussion-systematics}
The procedure described in Section Method is based on four main assumptions about the shape of the ISM in PPDV space and the relation between the two input datasets.
The first assumption is that the distribution of the ISM, in absolute units (e.g. number of protons), in PPD and PPV space is a linear scaling of the input reddening and gas line emission cubes.
The second assumption is that PPD and PPV cubes are projections of the same region of PPDV space, or equivalently that all the ISM contained in one is also contained in the other.
The third assumption is that the PPV cube can be reproduced by reprojecting the PPD cube -- ``lifting'' to PPDV space and projecting the result back down -- to PPV space.
This lift is assumed to have a simple description, a one-to-one function convolved with a Gaussian kernel in velocity space.
The fourth assumption is that voxels in PPD space that share a $\glon$ or $\glat$ boundary have similar central line-of-sight velocities.
These assumptions fail to hold, to different degrees, in the actual case we are applying the procedure to.
The reddening-to-matter and $\atomHI$- and CO-to-matter conversion is not a single linear relationship due to both variations in physical quantities such as the dust-to-gas ratio and non-linearities in the transfer function from amount of matter to amount of reddening or flux (e.g. $\atomH$ self-absortion).
The PPV cube includes matter on the far side of the Galaxy while the PPD cube stops at or before the level of the Galactic center.
The PPV cube can not be a reprojection of the PPD cube if the two are not tracing the same matter.
Furthermore, there will be voxels in the PPD cube where the lift to PPDV space is not the simple one we have assumed.
Finally, since the velocity field of the ISM does not have to be smooth or even continuous, the first due to the influence of departures from axisymmetry such as the Galactic bar and the second due to the fact that the ISM can contain shocks.
Despite the many ways in which these assumptions can and do break, we have been able to empirically confirm that the solution is, in many places, correct (see Section \ref{sec:KT-validation}).
Evidently, the KT procedure is robust to some bending of the assumptions.
Judging by the clustering of outlier HMSFRs in the inner Galaxy, there are some parts of the Galaxy where the assumptions are bent enough to start to cause problems.
In particular, the inner Galaxy is particularly likely to have mismatches between the PPD and PPV and complicated, rapidly spatially variable velocity fields.
Even there, KT does not completely and utterly fail -- the error rate is $\sim 30 \%$, not $\sim 100 \%$.
The inner Galaxy falls into a gray region in the parameter space of assumption bending where there is just enough of it to noticeably, but not absolutely, degrade the solution quality.
The 30\% mentioned above is the error rate specifically among PPD voxels that contain HMSFRs, which are not necessarily typical voxels.
The error rate of an arbitrary inner Galaxy voxel could potentially be quite different.
While we may expect the error rate of a voxel that contains an HMSFR and others near it to be similar, this is merely a heuristic argument.
This argument is also not useful for parts of the space that are not near HMSFRs.
The HMSFRs are all in dense molecular gas within a few tens of parsecs of the Galactic plane, meaning that there is very little we can say about the error rate in diffuse gas, especially at high Galactic latitudes.
To get a quantitative and more generally applicable estimate of the error rate, we would need either a less density-biased set of independent $\vlos(\glon, \glat, d)$ measurements or a set of artificial injection tests.
By an "artificial injection test", we mean a numerical experiment in which we artificially observe a model galaxy's ISM, reconstruct the $\vlos$ field from these artificial observations using KT, and compare the reconstructed and input model $\vlos$ fields.
To the best of our knowledge, there are no currently available catalogs of these sorts of less density-biased measurements, ruling out option one.
The steps involved in artificial injection tests, particularly simulating galaxies at sufficiently high resolution for our purposes and producing artificial observations in a way that includes the non-trivial systematics in the actual observations, are complicated enough to put option two beyond the scope of this work.
Both of these options are plausible directions for future work.
The 1.527 $\mu$m diffuse interstellar band (DIB), for instance, has been mapped over much of the northern sky by the APOGEE survey \citep{2015ApJ...798...35Z}.
Observations of this DIB towards APOGEE stars with known distances could potentially be used to build catalogs or even maps of $\vlos(\glon, \glat, d)$ using an independent dataset.
Artificial injection tests are conceptually straightforward, though they do require a substantial investment of time and computational resources.
While neither option is easy, we would argue that some combination of the two will be necessary before the KT-derived $\vlos$ map can be trusted away from the HMSFRs of \Reid.
\subsection{On the astrophysical plausibility and implications of our gas streaming motions}
\label{sec:discussion-plausibility}
While the exact shapes and magnitudes of the streaming motions in our velocity map depend on the assumed galactic parameters (Galactocentric distance to the sun, rotation curve, solar motion relative to the rotation curve), there is no choice of Galactocentric distance, rotation curve, and solar motion that can completely remove them (see Section \ref{sec:rotation_discussion} and Figure \ref{fig:six_pies}).
If we limit ourselves to rotation curves that are flat outside of the inner galaxy (panels a, b, e, f of Figure \ref{fig:six_pies}), we see streaming motions on $\sim$ 1 kpc scales with magnitudes in excess of 15 km/sec.
These gas motions would be perfectly reasonable if the Milky Way were an obvious grand design spiral such as M51 \citep{Meidt_2013}.
Since it is not, we should first ask if they are astrophysically plausible in a more humble sort of spiral galaxy such as our own.
We will argue that these spatially extended, high magnitude motions are indeed plausible.
Within 2 kpc of the Sun, the RAVE survey has detected stellar streaming of comparable spatial extent and magnitude (S2012). They interpret these motions as evidence of the gravitational influence of the Perseus spiral arm.
While the dynamics of stars and gas in s spiral potential will not necessarily be the same, this spiral-induced stellar streaming offers a an explanation for gas streaming of an essentially appropriate magnitude.
Streaming motions 3-5 kpc away in the inner galaxy can be associated with the dynamical influence of the Galactic bar.
This interpretation is bolstered by the similarity of our non-flat rotation curve to that of Clemens (198?), whose shape and deviation magnitude 3-5 kpc from the Galactic center can be explained by the proximity of the bar to the first quadrant tangent points (Incorrect Rotation Curve).
Our rotation curve is in general also quite similar to many that have been measured by the GHASP survey (CITE).
The GHASP survey used $\Halpha$ measurements to derive rotation curves for NUMBER disk galaxies with NUMBER kpc of the Milky Way with masses ranging from NUMBER to NUMBER times the mass of the Milky Way.
These rotation curves often show the same sort of $\sim 15$ km/sec, $\sim 1$ kpc velocity ``corrugations'' that we see in our rotation curve.
The GHASP collaboration derived separate rotation curves for the advancing and receding sides of their galaxies.
Often, the central radii of the velocity corrugations are slightly different on different sides of a galaxy.
These shifts imply that these deviations from flat rotation are caused by non-axisymmetric perturbations, such as spiral arms and/or bars.
Combining these pieces of circumstantial intra- and extragalactic evidence, we can conclude that $\sim 1$ kpc-scale $\sim 15$ km/sec gas streaming motions in a non-grand design spiral galaxy can not immediately be dismissed as unphysical.
We can also loosely associate these motions with spiral arms and the bar at large (small) and small (large) Galactocentric (heliocentric) radii.
If these motions are real, what conclusions can we draw from them without fixing a galactic parameter set? One immediate implication of these streaming motions is that improving the accuracy of kinematic distances will require accounting for large-scale motions beyond the rotation curve. Junichi+ (20??), e.g., estimated the distance distortion induced by assuming a flat rotation curve in a model galaxy non-circular HMSFRs velocities comparable to those of the Reid+09 HMSFRs.
THERE ARE IMPLICATIONS
The agreement between the gas and expanded HMFSR sample's line-of-sight velocities (see Section US WINNING) suggests that the non-circularity that the non-circularity of the HMSFR's velocities is real.
An argument against that possibility appeared shortly after the publication of Reid+ (2009), an early analysis of an a subset of the full \Reid sample, in McMillan and Binney (2010); this argument has been repeated with some frequency in discussions of discrepancies between Galactic parameters derived from the HMSFRs and other tracers (the usual suspect and hopefully someone else too).
If the HMSFRs' motions are real and due to orbital ellipticity, the argument goes, then the implied ellipticity of their orbits is higher than that of stars that are young but already out of their birth clouds.
Since HMSFRs are temporaly closer to collisional, and hence dynamically colder, dense gas than the collisionless, and hence dynamically hotter, young stars, this dynamical state of things is not physically plausible.
McMillan and Binney (2010) suggest some alternatives, including a bias towards bluer line-of-sight velocities, to reduce the required orbital ellipticity and bring the velocity dispersions of the two ostensibly similar populations closer together.
The agreement between the expanded HMSFR sample's line-of-sight velocities and our (completely independently) derived map suggests that the HMSFRs deviations from flat rotation are real and not due to measurement biases.
The spatial scale and similarity to spiral- and bar-induced streaming of non-circular motions in our map also suggests a resolution to the velocity dispersion inconsistency. The HMSFRs are moving coherently, with the gas they are embedded in, rather than randomly. Once the gas and the embedded HMSFRs leave the perturbations they are currently passing through, their apparent velocity dispersion as a population will presumably fall to something closer to that of the current generation of young stars.
Conversely, the McMillan and Binney argument and the agreement between the HMSFR velocities and our map combine into yet another piece of circumstantial evidence in favor of our streaming motions being associated with spiral- and bar-induced perturbations.
If interpreting the HMSFR motions as stable elliptical orbits is out of the question, then interpreting the (if anything dynamically even colder) gas motions as stable elliptical motions must be even less acceptable.
\section{Conclusion}
\label{sec:conclusion}
we have an awesome map
a list of reasons this map is awesome
\selectlanguage{english}
\FloatBarrier
\bibliographystyle{plain}
\bibliography{bibliography/converted_to_latex.bib%
}
\end{document}