# Molecular Gas in the Nuclear Region of NGC 6240

Adalyn Fyhrie  CASA University of Colorado Boulder
389 UCB
Boulder, CO 80309, USA  Jason Glenn  NASA Goddard Space Flight Center
Code 665
8800 Greenbelt Road
Greenbelt, MD 20771, USA  CASA University of Colorado Boulder
389 UCB
Boulder, CO 80309, USA  Naseem Rangwala  NASA Ames Research Center
Space Science & Astrobiology Division
MS 245-6, Bldg 245, Rm 107F
Moffet Field, CA 94035, USA  Jordan Wheeler  CASA University of Colorado
Boulder
389 UCB
Boulder, CO 80309, USA  Sara Beck  Tel Aviv University
School of Physics and Astronomy
Ramat Aviv 69978, Israel  John Bally  CASA University of Colorado Boulder
389 UCB
Boulder, CO 80309, USA

###### Abstract

NGC 6240 is a luminous infrared galaxy in the local universe in the midst of a
major merger. We analyze high-resolution interferometric observations of warm
molecular gas using CO J = $3-2$ and $6-5$ in the central few kpc of NGC 6240
taken by the Atacama Large Millimeter Array. Using these CO line observations,
we model the density distribution and kinematics of the molecular gas between
the nuclei of the galaxies. Our models suggest that a disk model represents
the data poorly. Instead, we argue that the observations are consistent with a
tidal bridge between the two nuclei. We also observe high velocity redshifted
gas that is not captured by the model. These findings shed light on small-
scale processes that can affect galaxy evolution and the corresponding star
formation.

NGC 6240 — LIRG — LIME — Radiative Transfer — ALMA — CO — tidal bridge — AGN —
outflow

††journal: ApJ††software: CASA (McMullin et al. (2007)), LIME (Brinch &
Hogerheijde (2010))

##  1 Introduction

NGC 6240 (Wright et al. (1984); Thronson et al. (1990)) is a unique galaxy in
the local universe (z = 0.02448) in the midst of a major merger event (Fried &
Shulz, 1983) that is triggering high star formation rates (Genzel et al.
(1998); Tecza et al. (2000)) and active galactic nuclei (AGN) activity in its
two progenitor nuclei (Vignatti et al., 1999), separated with a projected
distance of around 1 kpc. This activity results in a far-infrared (FIR)
luminosity L${}_{FIR}\approx$1011.8 L⊙ (Sanders et al. (1988); Thronson et al.
(1990); Sanders & Mirabel (1996)) that classifies it as a luminous infrared
galaxy (LIRG), just below the threshold that would classify it as an ultra-
luminous infrared galaxy (ULIRG). Its luminosity is expected to cross this
threshold when a second starburst is triggered during final coalescence (Engel
et al., 2010). As such, NGC 6240 presents an excellent opportunity to study in
fine detail the processes that power ULIRG activity.

The nuclear region ($\textless$1 kpc) of LIRGs and ULIRGs is a critical
location to study as it is often the location of the highest star formation
rates and the origin of many stellar feedback processes, such as stellar winds
and AGN outflows. Correspondingly, in the most luminous LIRGs
(L${}_{IR}\textgreater$ 6$\times$1011 L⊙) the majority of the mid-IR emission
originates from this central kpc region (Alonso-Herrero, 2013). The dynamics
in this central region are complicated by the galaxy merger interactions that
trigger the processes that lead to the high IR luminosity of most LIRGs
(Lonsdale, Farrah & Smith, 2006). The nuclear molecular gas dynamics of NGC
6240 are no exception, with a concentration of molecular gas that peaks in
emission between the two nuclei while the stars and dust are found to be
concentrated around the two nuclei (Tacconi et al. (1999); Tecza et al.
(2000); Engel et al. (2010)).

The molecular gas dynamics of NGC 6240’s nuclear region have been studied and
modeled for decades without consensus. Tacconi et al. (1999) observed a
velocity gradient in CO J = $2-1$ in the nuclear region and modeled it as a
rotating disk. A similar disk model was used to describe the motions of HCN by
Scoville et al. (2014). This inter-nuclear disk model has both been used to
support other authors’ observations of molecular gas (e.g. Iono et al. (2007)
in CO J = $3-2$ and HCO+ ($4-3$)) and been claimed to be unphysical in the
context of other observations (e.g. Gerssen et al. (2004) in H$\alpha$+[N
II]). Alternate geometries have been proposed to explain this central
molecular gas, including a tidal bridge connecting the two nuclei (Engel et
al., 2010) and the origin site for a warm molecular outflow (Cicone et al.,
2018). Even more recently, Treister et al. (2020) present observations of the
nuclear region in CO J = $2-1$ in unprecedented detail with angular
resolutions of 0.03′′. They find the central region to be a clumpy
concentration of gas dominated by a high velocity outflow and a gas bridge
connecting the two nuclei and interacting with the stellar disk kinematics.

With the lack of consensus on the geometry of the inter-nuclear molecular gas,
updated detailed modeling of the central nuclear region using multiple new
line transition observations in high resolution is needed. New telescopes have
allowed the nuclear region of NGC 6240 to be observed in unprecedented detail
and have enabled this modeling. In this paper we use high resolution CO J =
$6-5$ and CO J = $3-2$ observations from the Atacama Large Millimeter Array
(ALMA) to model the velocity profile and density distribution of the central
molecular gas using the Line Modeling Engine (LIME) (Brinch & Hogerheijde,
2010). These higher-J transitions trace shock excited warm molecular gas that
was detected but unresolved in observations by Herschel (see, e.g. Kamenetzky
et al. (2014)). LIME produces emission line profiles based on three-
dimensional velocity and density distributions provided to the modeling code.
We alter these three-dimensional distributions to match the modeled to
observed emission. We find multiple fiducial models, with the model whose
emission most closely matches observations suggesting the central molecular
gas is flat and pancake-like with a high velocity dispersion compared to its
rotational velocity, suggesting that it is not a self-gravitating, rotating
disk but instead a transient tidal bridge connecting the two nuclei.

The residuals from the fiducial model definitively show high velocity gas
($\textgreater$ 300 km/s) that is not associated with the nuclear pancake-like
concentration of gas. This high-velocity gas exists in others groups’
observations and is potentially associated with outflows postulated therein
(Cicone et al. (2018); Müller-Sánchez et al. (2018);Treister et al. (2020)).
The study of the high velocity gas is outside the scope of this work, which
aims to describe the properties of the majority of the nuclear molecular gas.
The presented observations also reveal extended emission features, some of
which correspond to previously studied filamentary structures, and others of
which are new to these observations.

Section 2 presents the ALMA observations of CO J = $3-2$ and J = $6-5$ and
continuum observations at 345 and 678 GHz. We analyze these observations in
Section 3 and calculate the mass of the dust from the continuum emission,
highlight observed extended molecular gas emission features, and explore the
velocity structure of the gas. Section 4 presents a non-local thermodynamic
equilibrium model created using LIME that fits the nuclear molecular gas’s
velocity and density distributions. We find that the molecular gas between the
two nuclei is unlikely to be a self-gravitating disk, but could instead be a
tidal bridge in light of the model findings.

##  2 Observations

Observations of CO J = $3-2$ and J = $6-5$ were completed using ALMA with
baselines of 1.6 km and 460 m, respectively. These are the first observations
of the nuclear region of NGC 6240 in CO J = $6-5$ and the highest resolution
observations to date in CO J = $3-2$. The CO J = $3-2$ observations were
completed during Cycle 2 for project number 2013.1.00813.S. The CO J = $6-5$
observation was completed for project number 2015.1.00658.S during Cycle 3.
The longest baselines, full-width half-maxima (FWHM) of the beams, channel
widths, reported channel root mean squared (RMS), central frequency of
continuum observations, bandwidth (BW) of continuum observations, and RMS of
continuum observations are reported in Table 1. For the integrated moment 0
maps, $\sigma_{integrated}$ is calculated using line channel RMS values
$\sigma_{line~{}channel}$ from Table 1. The RMS in the moment map is
calculated using $\sigma_{integrated}$ =
$\sqrt{n_{chan}}\sigma_{line~{}channel}$ where $n_{chan}$ is the number of
channels included in the integration.

To check the reliability of the provided reduced data products, we re-imaged
the data using the Common Astronomy Software Applications (CASA) software
package and the National Radio Astronomy Observatory (NRAO) provided imaging
script. Re-imaging was completed with natural, uniform, and Briggs weighting
schemes and user-created masks. For CO J = $3-2$ no new structure emerged in
the molecular gas when re-imaging the data, nor could we improve upon the
noise. Therefore, we deemed the NRAO provided data products sufficient for
analyses of the CO J = $3-2$ molecular gas.

Re-imaging was required for the CO J = $6-5$ and both continuum observations
due to contamination in the continuum maps from extremely high-redshift with
respect to systematic ($\sim$ 400 $-$ 740 km/s) molecular gas. This CO line
contamination created a false continuum source between the two nuclei with the
same peak flux density as the southern nucleus. Upon re-imaging, this source
disappeared. We use Briggs weighting with robust parameter 0.5 for the re-
imaged maps to match the weighting scheme used for the NRAO-provided molecular
gas data products.

Below, we introduce all observations with further details of all figures
discussed in the following sections. At the distance of NGC 6240, 1′′
corresponds to 500 pc of projected distance. AGN locations from Hagiwara et
al. (2011) are included as crosses in the figures.

The two continuum maps are shown in Figure 1, revealing continuum emission in
the vicinity of the nuclei as observed by previous authors (e.g. Scoville et
al. (2014)). The concentration around the southern nucleus is much brighter
than the northern concentration, with flux densities reported in Table 2. In
detail, the morphologies appear to be different because of differing
synthesized beams and signal-to-noise ratios in the two maps; however, there
also differences that arise from real structure. In the south, the 678 GHz
emission has a peak on the 345 GHz nucleus and in addition a stronger emission
clump, $\sim$70$\%$ higher intensity, $\sim$0.4′′ north of the southern
nucleus and collinear with the two nuclei. The bright clump was seen also by
Scoville et al. (2014). This fidelity of continuum observations lends
confidence to the associated line observations presented in this paper. The
differences between the 678 and 345 GHz observations likely derives from a
gradient of dust temperature and optical depth along the axis between the two
nuclei. The 678 GHz source may arise in a hot spot on the axis between the two
nuclei or in a local fluctuation in the optical depth; it should be observed
further. The potential temperature and optical depth gradient, with
temperature increasing to the north, is supported by the CO J = $6-5$ moment 0
centroid corresponding better than the CO J = $3-2$ moment zero centroid to
the 678 GHz dust emission.

Moments 0, 1, and 2 (integrated flux density, average velocity, and velocity
dispersion) for CO J = $3-2$ are plotted in Figure 2 with the 345 GHz
continuum contours included. Similarly, Figure 3 shows the moments of CO J =
$6-5$ with the 678 GHz continuum contours. Observations in both J = $3-2$ and
J = $6-5$ show a concentration of molecular gas between the two nuclei, as
observed previously in other observations of the molecular gas (e.g. Tacconi
et al. (1999), Scoville et al. (2014), Treister et al. (2020)). Also similar
to previous observations, the CO moment 1 maps show a velocity gradient from
highly redshifted with respect to systematic ($\textless v\textgreater\sim$
350 to 400 km/s) to blueshifted with respect to systematic ($\textless
v\textgreater\sim$ -100 to -150 km/s) along a position angle of approximately
34∘. Extended structure including two dim features of molecular gas to the SE
and SW of the southern nucleus are observed in the CO J = $3-2$ maps that have
also been observed in other wavelength bands, such as H2 (Max et al., 2005),
Fe XXV, and H$\alpha$ (Wang et al., 2014) compared in Figure 7.Although these
lines trace different components of the ISM (molecular, shocked, ionized
respectively) they are all associated with star formation activity. We also
observe new extended emission to the north in the CO J = $3-2$ observation.

Channel maps for the CO J = $3-2$ and CO J = $6-5$ observations are shown in
Figure 4 and Figure 5, respectively. Radio velocities of the observations are
calculated by the CASA software relative to the central observed frequencies,
345.796 GHz for CO J = $3-2$ and 692.2 GHz for CO J = $6-5$. The systemic
velocity of NGC 6240 is taken as the velocity corresponding to the peak flux
density of the galaxy-integrated line profile from the observations of CO J =
$3-2$ and CO J = $6-5$: 7180 $\pm$ 20 km/s for CO J = $3-2$ and 7150 $\pm$ 20
km/s for CO J = $6-5$.

The velocity structure of the inter-nuclear gas is most apparent in the higher
sensitivity CO J = $3-2$ channel map, showing emission at $\sim$-500 km/s with
respect to systematic near the southern nucleus and emission at $\sim$ 650
km/s with respect to systematic located between the two nuclei. The CO J =
$6-5$ channel map shows a similar velocity structure to that of the CO J =
$3-2$. This velocity structure is best resolved in the nuclear region between
the two AGN, where the gas emission is brightest.

Table 1: Parameters of observations received from ALMA.  Parameter | CO J = $3-2$ observation | CO J = $6-5$ observation
---|---|---
Configuration (Baseline) | Extended (1.6 km) | Extended (460 m)
Peak Molecular Flux Density Location | (16:52:58.898, +02.24.03.672) | (16:52:58.898, +02.24.03.756)
Peak Flux Density | 100 Jy/beam km/s | 308 Jy/beam km/s
Beam FWHM | 0.34 x 0.15′′ | 0.32 x 0.23′′
Beam PA | -67o | -83o
Line channel width | 20 km/s | 20 km/s
Line channel RMS $\sigma_{line~{}channel}$ | 2.52 mJy/beam | 14.5 mJy/beam
Continuum band center $\nu$ | 345 GHz | 678 GHz
Bandwidth | 7.5 GHz | 6 GHz
Continuum RMS $\sigma_{cont}$ | 0.21 mJy/beam | 2.1 mJy/beam
Figure 1: Left: Colors: The 345 GHz continuum observation. Solid White
Contours: 4, 6, 10, and 14 $\sigma_{cont}$ contours. Black Crosses: The
locations of the two known AGN (Hagiwara et al., 2011). The beam FWHM contour
is plotted in the lower right of the image. Right: Colors: The continuum
observation at 678 GHz. White Contours: 5, 10, and 15 $\sigma_{cont}$
contours. Black Crosses: The locations of the two known AGN. The beam FWHM
contour is plotted in the lower right of the image, 1′′ corresponds to 500 pc,
east is left, and north is up.” Figure 2: Upper Left: Full map of moment 0 of
CO J = $3-2$ with dashed white contours corresponding to 10, 25, 50, 100, 200,
300, 400, 500, and 600 $\sigma_{integrated}$ of the integrated line emission,
solid black contours show the 345 GHz continuum at 4, 6, 10, and 14
$\sigma_{cont}$, and crosses show the location of the two AGN (Hagiwara et
al., 2011). Upper Right: The same, but zoomed to show only the nuclear region.
Lower Left: Moment 1 in the nuclear region with channels below 5
$\sigma_{line~{}channel}$ masked out prior to moment calculation. Lower Right:
Moment 2 in the nuclear region with channels below 5 $\sigma_{line~{}channel}$
masked out prior to moment calculation. The beam FWHM contour is plotted in
the lower right of each image, 1′′ corresponds to 500 pc, east is left, and
north is up. Velocities are calculated relative to the average radio velocity
of CO J = $3-2$ for the entire observation, 7180 km/s, calculated relative to
its rest frequency 345.796 GHz.  Figure 3: Upper Left: Moment 0 of CO J =
$6-5$, with dashed contours corresponding to 10, 25, 50, 100, 200
$\sigma_{integrated}$ of the integrated line emission and solid contours show
the 678 GHz continuum at 5, 10, and 15 $\sigma_{cont}$. Channels below
5$\sigma$ are masked out of the integration. Upper right: Moment 1 of CO J =
$6-5$, with pixels below 5 $\sigma_{line~{}channel}$ masked out prior to
moment calculation. Lower panel: Moment 2 of CO J = $6-5$, with pixels below 5
$\sigma_{line~{}channel}$ masked out prior to moment calculation. The beam
FWHM contour is plotted in the lower right of each image. Channels below
5$\sigma_{line~{}channel}$ are masked out of the integration for all maps, 1′′
corresponds to 500 pc, east is left, and north is up. Velocities are
calculated relative to the average velocity of CO J = $6-5$ for the entire
observation, 7150 km/s calculated relative to the CO J= $6-5$ rest frequency
691.473 GHz.  Figure 4:  The channel map of primary beam corrected, continuum
subtracted CO J = $3-2$ observations show every third line channel (each panel
is separated by 60 km/s). White contours: 10, 25, 50, 75, and 100 times
$\sigma_{line~{}channel}$ of the CO J = $3-2$ observations. White crosses: the
locations of the AGN. The beam FWHM contour is plotted in the bottom right of
each figure, 1′′ corresponds to 500 pc, east is left, and north is up.
Velocities are calculated relative to the average velocity of CO J = $3-2$ for
the entire observation, 7180 km/s calculated relative to its rest frequency
345.796 GHz.  Figure 5:  The channel map of primary-beam-corrected, continuum-
subtracted CO J = $6-5$ showing every third line channel. White contours: 5,
10, 25, 50, and 75 times $\sigma_{line~{}channel}$ of the CO J = $6-5$
observations. White crosses: The locations of the AGN. The beam FWHM contour
is plotted in the bottom right of each figure, 1′′ corresponds to 500 pc, east
is left, and north is up. Velocities are calculated relative to the average
velocity of CO J = $6-5$ for the entire observation, 7150 km/s calculated
relative to its rest frequency 691.473 GHz.

##  3 Continuum and Extended Molecular Gas Emission Analysis

###  3.1 Mass from Continuum

To calculate the mass of the dust $M_{d}$ from the continuum observations we
use

$M_{d}=S_{\nu}D_{L}^{2}/\kappa_{\nu}B_{\nu}(T)$ (1)

(Casey, 2012), where $S_{\nu}$ is the flux density in the continuum frequency
band, $D_{L}$ is the luminosity distance of 108 Mpc, $\kappa_{\nu}$ is the
dust mass opacity coefficient, and $B_{\nu}(T)$ is the blackbody emission in
the continuum frequency band for a dust temperature T. For this calculation,
we use the 345 GHz as the measure of $S_{\nu}$ because it should have a lower
optical depth than the 678 GHz continuum observation. The total $S_{345GHz}$
for this observation is 27 mJy, 18$\%$ of the galaxy-integrated flux density
measured at 850 $\mu$m by SCUBA of 150 mJy (Klaas et al., 2001). This flux
density recovery is comparable to Scoville et al. (2014) who measured 18-24
mJy at 340 GHz for a comparable beam size and sensitivity with ALMA in Cycle
0. The SCUBA beam is as large as our primary beam of the ALMA observations
with a diameter of 15′′ (Klaas et al., 2001), and as such we can expect that
their measurement includes extended emission that is lost in our high-
resolution observations. Therefore we expect our value of $S_{\nu}$ to be
lower than that measured with SCUBA.

For $B_{\nu}(T)$ we choose a dust temperature of 56 K, the dust temperature
fit in Kamenetzky et al. (2014) using a greybody fit to Herschel-SPIRE, IRAS,
Planck, SCUBA, and ISO photometry for NGC 6240. This dust temperature is a
galaxy-averaged property since Kamenetzky et al. (2014) used observations with
beam FWHM between $\sim$ 17′′ and 45′′, at least 100 times larger than the
ALMA observations. It is likely that the dust temperature in NGC 6240’s
nuclear region is higher than the galaxy-averaged temperature due to the
influences of concentrated star formation and AGN luminosity. However, we do
not have flux density measurements across sufficient wavelengths in this
central region to independently calculate the dust temperature. From James et
al. (2002) we use $\kappa_{850}$ = 0.07 m2kg-1 and
$\kappa_{\nu}\propto\nu^{2}$ to find $\kappa_{870}$ = 0.067 m2kg-1.

The calculated masses for the regions outlined in Figure 6 are tabulated in
Table 2. The total mass of the dust is calculated to be 1.2$\times$107 M⊙. As
a check on this total dust mass we compare to the galaxy-integrated dust mass
of 5$\times$107 M⊙ found in Kamenetzky et al. (2014) from their dust SED model
described earlier in this section. Our calculated dust mass is 24$\%$ of this
value, consistent with filtered flux from the ALMA observations.

Using a gas-to-dust mass ratio of 100 we can convert the dust masses from
Table 2 to a total gas mass, 1.2$\times$109 M⊙. The mass values derived from
the continuum emission are similar to the values derived from the sub-mm
continuum (235 GHz) in Treister et al. (2020) of 2.8$\times$109 M⊙. The gas
mass is converted to a column density $N_{H_{2}}$ by dividing the gas mass by
the mass of molecular hydrogen, the projected size of the region on the sky,
and multiplying by the fraction of the total mass assumed to be molecular
Hydrogen (0.73, assuming solar metallicity Hollenbach & Thronson (1987)). The
average column density for the entire nuclear region is 3$\times$1022 cm-2,
while the concentrations around the two nuclei both have higher column
densities of 1.5$\times$1023 cm-2. This value is consistent with the findings
of Tacconi et al. (1999) of $N$(H2) $\sim$ 1–2$\times$1023 cm-2.

Figure 6: Regions for which we calculate the dust mass. Color map: Extended
configuration 345 GHz continuum map. Dashed contours: CO J = $3-2$ integrated
line emission. Masses for these regions are tabulated in Table 2. The beam
FWHM contour for the 345 GHz continuum is plotted in the bottom right.

Table 2: Integrated continuum flux densities $S_{\nu}$, derived dust masses from the 345 GHz continuum calculated using Equation 1, gas masses assuming a gas-to-dust mass ratio of 100, and associated column densities of the regions outlined in Figure 6.  Region | Description | 678 GHz $S_{\nu}$ [Jy] | 345 GHz $S_{\nu}$ [Jy] |  $M_{d}$ [M⊙] |  $M_{gas}$[M⊙] |  $N_{H_{2}}$ [cm-2]
---|---|---|---|---|---|---
1 | Entire Continuum Observation | 0.27 | 0.027 | 1.2$\times$107 | 1.2$\times$109 | 4.0$\times$1022
2 | Northern Concentration | 0.03 | 0.002 | 9.0$\times$105 | 9.0$\times$107 | 2.1$\times$1023
3 | Southern Concentration | 0.116 | 0.013 | 5.8$\times$106 | 5.8$\times$108 | 2.4$\times$1023

###  3.2 Extended Molecular Gas Emission

A dim trail of gas extending directly north of the central concentration,
which we call the “Northern Finger”, is observed in the higher-sensitivity CO
J = $3-2$ observation but not in CO J = $6-5$ or either continuum observation.
It is labeled in the CO J = $3-2$ moment 0 map in Figure 2. It is not
spatially coincident with observed H2 V = $1-0$ S(1) and S(5) emission (Max et
al. (2005), their Figure 12a). It dimly appears in their Figure 2a and 2c,
false color images including Keck $K^{\prime}$ band, $H$ band, $J$ band, F814W
filter from WFPC2 on $Hubble$ $Space$ $Telescope$ (HST), and the F450W filter
from WFPC2 on $HST$. Its southern half is spatially coincident with [O III]
$\lambda$ = 5007 and H$\alpha$ emission (Müller-Sánchez et al. (2018), their
Figure 1).

The Finger’s average velocity is 44 km/s with a FWHM of $\sim$ 300 km/s,
measured with an integrated spectrum of the emission in this region. The low
average velocity suggests this gas is unlikely to be an outflow. A nuclear
inflow would cause $\sigma_{v}$ to be elevated near the nuclei, for which
there is no conclusive evidence. The low velocity is consistent with a tidal
tail or bridge but the high FWHM does not conform to this idea. The finger
extends from the northern nucleus towards the northern dust lane of the
galaxy, which is interpreted as a tidal tail by other authors (Gerssen et al.
(2004); Yun et al. (2001)). The top edge of the finger is spatially coincident
with the eastern edge of the dust lane, suggesting a possible correspondence.
The dust lane is much wider and extends much farther beyond the extent of the
finger – approximately 7′′ wide (roughly East/West) and over 20′′ long
(roughly North/South), beyond the extent of our observations.

Figure 7:  Figure 10 from Max et al. (2005), “Keck AO narrowband image of the
nuclear region in the H2 1-0 $S$(1) emission line […] This image is 7.′′5 wide
by 5.′′5 tall and has a log color map. The eight contour levels are on a log
scale, with the first contour at 0.2 times the maximum flux (at the South1
nucleus) and the last contour at 0.02 times the maximum flux.”, included as a
comparison to the CO J = $3-2$ moment 0 map in Figure2.

There is also extended structure surrounding the two nuclei that align with H2
V = $1-0$ S(1) and S(5) concentrations presented in Max et al. (2005) (their
Figures 10 and 11, reproduced in our Figure 7): a molecular concentration to
the NE of the northern nucleus, a concentration to the SW of the southern
nucleus, and dim, diffuse emission extending to the SE of the southern
nucleus. The faint arms of H2 extending SE and SW observed in Max et al.
(2005) align with these faint arms extending SE and SW from the central
concentration and are also bright in Fe XXV and H$\alpha$ (Max et al., 2005;
Wang et al., 2014). Wang et al. (2014) postulates that these faint molecular
arms are molecular gas entrapped and shocked by the superwind caused by a
vigorous starburst in the southern nucleus. Similarly, Max et al. (2005) argue
that these arms are a thin layer of gas at the edges of soft X-ray bubbles
observed in Komossa et al. (2003), where a starburst driven wind is “driving
shocks or ionization fronts into the interstellar medium and surrounding
molecular clouds”. We therefore conclude that the CO arms to the SE and SW of
the southern nucleus are likely associated with the same starburst-driven
superwind shocked gas.

##  4 Molecular Gas Between the Nuclei: a Test Model Motivated by History and
Observations

In Tacconi et al. (1999) and in works since (e.g. Iono et al. (2007); Scoville
et al. (2014)), the velocity gradient observed in the central molecular
concentration was interpreted and modeled as a disk of gas between the two
nuclei. The presence or absence of a molecular disk in the nuclear region of
NGC 6240 is the basis for other authors’ arguments regarding important stellar
feedback processes, such as the outflow studied in Cicone et al. (2018). The
disk model has been claimed to be unphysical in other studies of molecular gas
in NGC 6240 (e.g. Gerssen et al. (2004)), but no further modeling has
definitively proven or disproven the internuclear disk model. More generally,
the nuclear regions of (U)LIRGs are critical to understand as they are often
the location of the highest star formation rates and the origin of many
feedback processes that affect star formation, such as stellar winds and AGN
outflows (Alonso-Herrero, 2013).

For these reasons, we model this nuclear region as a concentration of gas with
exponentially decreasing velocity and density profiles from the central peak
using the observations in CO J = $3-2$ and J = $6-5$ to explore in detail the
validity of this disk interpretation. The modeled gas concentration is not
constrained to a self-gravitationally bound disk geometry, rather, it is a
general concentration whose resultant velocity and density profiles can be
compared to the disk interpretation. Our fiducial model finds the molecular
gas concentration between the two nuclei is likely flat and pancake-like, with
a high velocity dispersion that could cause it to dissipate as quickly as
$\sim$ 0.3 Myr. Its geometry suggests this transient structure could be a
tidal bridge connecting the two nuclei.

Figure 8:  Contours of high and moderate velocity redshifted and blueshifted
gas for CO J = $3-2$ and J = $6-5$, plotted over corresponding integrated line
emission (moment 0). Contours are 0.2, 0.4, 0.6, and 0.8 of the maximum in
each map. Beam FWHM contours are plotted in the bottom right of panels b and
d. Crosses denote locations of the AGN. a: CO J = $3-2$ high velocity gas; red
contours: [250, 700] km/s, dashed blue contours: [-500,-200] km/s b: CO J =
$6-5$ high velocity gas; red contours: [250, 740] km/s, dashed blue contours:
[-420,-200] km/s c: CO J = $3-2$ moderate velocity gas; red contours: [0, 250]
km/s, dashed blue contours: [-200,0] km/s d: CO J = $6-5$ moderate velocity
gas; red contours: [0, 250] km/s, dashed blue contours: [-200,0].

###  4.1 Code for Generating the Test Model: the Line Modeling Engine (LIME)

To generate the molecular gas model that tests the validity of the disk
interpretation, we use the Line Modeling Engine (LIME) (Brinch & Hogerheijde,
2010). LIME uses non-local thermodynamic equilibrium radiative transfer and
molecule rotational energy level populations calculations to predict the line
and continuum emission from molecular clouds. The user defines a 3-D model
describing the density distribution of hydrogen molecules of the source, then
assigns 3-D temperature and velocity distributions to those molecules. The
gas-to-dust mass ratio, the abundance of CO relative to hydrogen, and the 3-D
distribution of the dust temperature are also set by the user. Finally, the
user defines the distance to the source to obtain the proper distance scale
per pixel and observed flux density.

The code does not require that the source is in local thermodynamic
equilibrium, and instead solves for population levels iteratively until the
model populations have converged at all grid points. After convergence, LIME
ray-traces photons to obtain an image of the modeled source at a user-defined
observing angle. This simulation methodology allows great flexibility in both
geometry and kinematics of the simulated source and minimizes the assumptions
made about the source to generate the observed line profiles.

To compare the observed data to the simulated model images, we use the CASA
package to smooth and continuum subtract the simulated images. The resulting
images have the same beam size as the observations. We must also account for
the redshift of the galaxy, which is not included in the LIME simulations. To
do this, we measure the central velocity of the galaxy-integrated line profile
from the observations of CO J = $3-2$ and CO J = $6-5$: 7180 $\pm$ 20 km/s for
CO J = $3-2$ and 7150 $\pm$ 20 km/s for CO J = $6-5$. These velocities are
subtracted from the observed moment 1 maps in Figures 2 and 3, which we can
compare to the simulations that are simulated with zero redshift.

###  4.2 Parameters of the Test Model

Using LIME, we model the 3-D central gas concentration with density and
velocity distributions with exponentially decaying density and velocity
profiles. Our 3-D model for the disk density follows the function

$n_{H2}=n_{H_{2},0}\exp(-r/r_{sh})\exp(-|z|/z_{sh}),$ (2)

where $r=\sqrt{x^{2}+y^{2}}$, $r_{sh}$ is the radial scale height of the
density, $z_{sh}$ is the vertical scale height, and $n_{H_{2},0}$ is the
central H2 number density of the gas concentration. The total gas mass within
the extent of our ALMA observations is set to a simulation parameter
$g_{mass}$, from which $n_{H_{2},0}$ is calculated based on the simulation
geometry. We also allow an overall density asymmetry by multiplying the
density on one half of the gas concentration (x$\textgreater$0) by a model
parameter $a_{n}$. For the CO J = $3-2$ model we also allow the extent of the
gas to differ for x$\textgreater$0 and x$\textless$0 by varying the scale
height on either side of the modeled concentration.

We model the velocity with

$v(r)=v_{circ}(1-\exp(-r/r_{circ})),$ (3)

where $v_{circ}$ is the circular velocity and $r_{circ}$ is the radial scale
length of the circular velocity. We use the exponential profile because sside
from the high velocity gas component, we find no sufficient justification for
deviation from a rotation curve that rapidly becomes flat with radius, (which
would be appropriate in the case of a self-gravitating disk).

In Equations 2 and 3, parameters except $r$ and $z$ are fit parameters that
can be tuned to match the resultant modeled line profiles to the observed
data. In addition to those parameters, the LIME model requires a gas
inclination, position angle, dust temperature $T_{dust}$, turbulent velocity
$v_{turb}$, gas temperature $T$, and gas-to-dust mass ratio (set to 100). All
available parameters, fiducial fit results, and acceptable ranges of fiducial
fit parameters can be found in Table 3 and are discussed in Section 4.3.

The position angle of the gas concentration is modeled as 34∘, the angle
between the maximum of the highly blueshifted gas and that of the highly
redshifted gas from Figure 8. This is close to the position angle used by
Tacconi et al. (1999) of 40∘ for their disk models. Thirty-four degrees was
chosen following methodologies laid out in Cicone et al. (2018) and Tacconi et
al. (1999). Following the methodology presented in Cicone et al. (2018), we
separate the observations into “quiescent” [-200, 250] km/s gas and high
velocity gas. Figure 8 shows the contours of the observed highly redshifted
(CO J = $3-2$: [250, 700] km/s, CO J = $6-5$: [250, 740] km/s) and highly
blueshifted (CO J = $3-2$: [-500, -200] km/s, CO J = $6-5$: [-420, -200] km/s)
gas in panels a and b. The centroids between these two velocity extremes (34∘)
is at a steeper angle than between the two nuclei ($\sim$19∘). The contours of
moderate velocity, or “quiescent” gases of velocity [-200, 250] km/s are shown
in in panels c and d. The velocity gradient observed in moment 1 for CO J =
$3-2$ and J = $6-5$ remains clearly present in the high velocity gas, with
distinct centroids of emission at a PA of $\approx$ 34∘. The lower velocity
gas has much more overlap between the redshifted and blueshifted components,
without as clear of a separation between centroids of emission. This lower
velocity gas aligns approximately with the semi-major axis of the central gas
concentration at a PA of $\approx$ 0∘, along the north/south axis. Due to the
distinct gradient present in the high velocity gas that more closely resembles
a traditional disk, we choose the position angle based on the highest velocity
gas. However, we find during modeling that the position angle does not
influence the fit parameters outside of acceptable ranges already incorporated
into the model.

The center of the modeled gas concentration is placed at the maximum of the
observed integrated line emission, following the methodology of Tacconi et al.
(1999). If the gas were a simple disk, the maximum of the line emission would
be expected to align closely with the turnover point of the velocity gradient.
As noted in Tacconi et al. (1999) and confirmed in this paper’s observations,
this is not the case. This offset could be caused by a number of phenomena,
including that the gas is not a simple disk. The offset could also be
explained by a gas concentration with a highly asymmetric mass distribution
leading to an offset in the maximum line emission. Optical depth effects could
also displace the velocity gradient turnover point from the maximum of the
line emission.

###  4.3 Test Model Fitting Strategy and Fiducial Fits

In NGC 6240 the CO J = $6-5$ emission has been found to be dominated by
relatively hot gas while the CO J = $3-2$ emission is dominated by by
cool/warm gas (Kamenetzky et al., 2014). The upper levels of the CO J = $6-5$
and J = $3-2$ CO lines are at very different energies, with the J = $6$ level
at $\sim$116 K and the J = $3$ level at $\sim$33 K; it is clear that J = $6-5$
line will arise in warmer gas than does the J = $3-2$. The true excitation
picture is of course more complicated than a simple model of ‘warm’ and ‘cool’
components. The J = $3-2$ line is to some extent sensitive to warm and hot
gas: it is observed in quite hot regions (e.g. Consiglio et al. (2017)) and in
starburst galaxies is usually in excess of what the true ‘cold’ gas component
predicts (Kamenetzky et al., 2014). But although there is probably some
overlap between the J = $6-5$ and J = $3-2$ emitting gas, in bulk they
preferentially probe two different temperature regimes.

We therefore generate separate models for a ‘warm’ component matching the J =
$6-5$ profile and a ‘cool’ component matching the J = $3-2$ line. Models of
the two temperature components have the same velocity structure, gas
inclination, dust temperature and position angle. The gas mass, density and
temperature distribution are different for the two models. In effect we
require the cool and warm gas to have the same morphology, but permit the gas’
extent and density to vary within that morphology.

Due to the simpler nature of the CO J = $6-5$ line profiles compared to those
of CO J = $3-2$, we first use CO J = $6-5$ observations to constrain the
velocity structure, inclination and position angle, and dust temperature. We
then applied the fiducial fit parameters from the CO J = $6-5$ model to find
the gas density distribution and gas temperature of CO J = $3-2$.

The model parameters are found by eye due to the complexity of the nuclear gas
structure and dynamics, the large parameter space, and the number of fit
parameters available to the model. The parameters are chosen to minimize the
apparent differences between the modeled line profiles and the observed line
profiles extracted at points separated by a beam FWHM along the major and
minor axes of the modeled gas. Once a fiducial fit is found, ranges on
acceptable fit parameters are found by varying one parameter at a time until
the model’s line profiles are no longer acceptably close to the observed line
profiles. “Acceptable” fits are determined by a combination of factors: first,
the modeled line profiles on average cannot appear to differ from the observed
line profiles by more than $\sim$ 25$\%$ in height, width, or central
velocity. The fits focus on the brightest region along the gas’ major axis.
Second, the shape of the modeled line profiles should approximately match that
of the observed line profiles.

Many of the parameters can be well constrained by the observed line profiles,
for example $r_{sh}$ by the rate of decay of the line profile heights with
radius and $v_{turb}$ by their widths. The singly peaked profiles and lack of
gradient in their central velocities constrains the inclination to nearly
face-on models. Inclinations between 45∘ and 135∘ are strongly disfavored
because they showed doubly peaked profiles. Other parameters are more
difficult to constrain. $v_{circ}$ and $r_{circ}$ are the least constrained
because they are anti-correlated (see Equation 3) leading to large fit ranges
as they trade off against each other. However, $v_{circ}$ is not wholly
unconstrained. It cannot be so high as to shift the line centers, but must be
high enough to create the singly peaked profiles. The temperature $T$ is
constrained by the brightness of the lines, with higher temperatures (and
therefore higher CO rotational excitation) corresponding to brighter lines.

All available parameters, fit values, and acceptable ranges of model
parameters are reported in Table 3. The fiducial fit model whose emission
profiles most closely match the observed profiles is a nearly face-on pancake-
like distribution of gas, with a scale height $z_{sh}$ of 20-60 pc and a
radial scale length of 390 pc at an inclination of 22.5∘ from face-on. The
line profiles extracted from these models along the semi-major axis of the gas
concentration for CO J = $6-5$ and J = $3-2$ are compared to those from the
data in Figure 9.

Table 3: Fiducial LIME model values for CO J = $3-2$ and J = $6-5$ (column 2), as well as the ranges in parameter space that resulted in acceptable models (column 3). Velocity structure, gas inclination, and dust temperature is shared between the two models while gas mass, density distribution, and gas temperature vary between the two. Gas masses are calculated by integrating the modeled density over the observed extent of the gas (the high signal-to-noise ratio areas in Figures 2 and 3).  Variable Name | Fit Value | Range
---|---|---
Shared Fiducial Model Values
Inclination | 22.5∘ |  $0^{\circ}-45^{\circ}$, $135^{\circ}-180^{\circ}$
Position Angle | 33.7∘ | $30^{\circ}-40^{\circ}$
Tdust | 30 K |  $20$ K $-~{}100$ K
vturb | 140 km/s |  $100$ K $-$ $160$ km/s
vcirc | 100 km/s | 50 km/s $-$ 300 km/s
rcirc | 200 pc | 0 pc $-$ 500 pc
CO J = $6-5$ Fiducial Model Values
Gas Mass | 6.0$\times$108 M⊙ |  $4.0\times 10^{8}-7.8\times 10^{8}$ M⊙
rsh | 390 pc | 240 pc $-$ 500 pc
zsh | 20 pc | 18 pc $-$ 60 pc
$a_{n}$ | 1.0 | 0.6 $-$ 1.0
T | 2,000 K | 1,500 K $-$ 3,000 K
CO J = $3-2$ Fiducial Model Values
Gas Mass | 7.6$\times$108 M⊙ |  $2.9\times 10^{8}$ M${}_{\odot}-9.2\times 10^{8}$ M⊙
rsh | 350 pc | 250 pc $-$ 400 pc
rsh,x<0 | 250 pc | 200 $-$ 400 pc
zsh | 13 pc | 10 pc $-$ 16 pc
$a_{n}$ | 1.0 | $0.8-1.8$
T | 600 K | 400 K $-$ 1500 K
Figure 9: The line profiles for the fiducial fit CO J = $6-5$ (top) and J =
$3-2$ (bottom) LIME models compared to the observed CO J = $6-5$ and J = $3-2$
line profiles extracted along 34∘, the angle between the maximum of the highly
blueshifted and redshifted gas. Dashed Lines: model line profiles at the
location color-coded in the inset image. Solid Lines: Observed line profiles
at the location color-coded in the inset image. Inset Image: The color map is
the moment 0 map for the corresponding observation, and contours show the
model’s integrated line emission.

A wide range of parameter space resulted in fits of similar quality to our
fiducial model. As an example, we include a thicker-disk model with the
parameters in Table 4 with the corresponding modeled line profiles in Figure
10. This model is thicker than the models presented in Table 3, with a
vertical scale height $z_{sh}$ of 60 pc instead of 20 pc, a gas concentration
that is more horizontally concentrated ($r_{sh}$ = 270 pc instead of 390 pc),
and with a lower ratio of turbulent to circular velocity ($v_{turb}/v_{circ}$
= 0.44 instead of 1.4). While this thicker model’s parameters are all within
the thin disk’s acceptable ranges and still capture the majority of the
observed emission (line profile width, intensity, and in some cases line
profile shape), the line profiles are more doubly peaked and less
representative of the observed line profiles.

Table 4: Fiducial LIME model values for the thicker modeled CO J = $6-5$ gas concentration compared to the thinner model. Parameters not listed are the same as in Table 3. The thick model’s line profiles are presented in Figure 10. The thin model’s values are taken from Table 3 and the line profiles are presented in Figure 9.  Variable Name | Thick Disk Fit Value | Thin Disk Fiducial Fit Value
---|---|---
zsh | 60 pc | 20 pc
rsh | 270 pc | 390 pc
vturb | 110 km/s | 140 km/s
vcirc | 250 km/s | 100 km/s
rcirc | 10 pc | 200 pc
Gas Mass | 7.0$\times$108 M⊙ | 6.0$\times$108 M⊙
T | 2,500 K | 2,000 K
Figure 10: The line profiles for the thicker fiducial fit CO J = $6-5$ LIME
model compared to the observed CO J = $6-5$ line profiles extracted along 34∘,
the angle between the maximum of the highly blueshifted and redshifted gas.
The model parameters of the gas concentration are reported in Table 4. Dashed
Lines: model line profiles at the location color-coded in the inset image.
Solid Lines: Observed line profiles at the location color-coded in the inset
image. Inset Image: The color map is the moment 0 map for the corresponding
observation, and contours show the model’s integrated line emission.

###  4.4 Fiducial Model Parameters in Context

We place the model parameters in a physical context by comparing them to the
findings of Kamenetzky et al. (2014) who use Herschel Space Observatory and
other supplementary observations (see their Table 16) of CO J = $1-0$ to
$13-12$ to model gas temperatures, H2 column densities, gas masses, and H2
number densities for NGC 6240. We find our CO J = $3-2$ model has a
temperature and mass consistent with their calculations. The CO J = $6-5$
model is on the high end of both temperature and mass, but the total modeled
mass of the gas is consistent with other findings for NGC 6240.

Our models indicate the CO J = $6-5$ emission is dominated by relatively hot
(T = 2,000 K) gas, and the CO J = $3-2$ emission by warm (T = 600 K) gas.
Kamenetzky et al. (2014) modeled the CO emission as a combination of cold (16
K with 1$\sigma$ range [5, 50] K) and hot (1,260 K with 1$\sigma$ range [790,
2,000] K) gas. They found the CO J = $3-2$ emission is dominated by a mixture
of both cold and hot gas, while the CO J = $6-5$ emission is produced
primarily by the hot gas component. Our CO J = $3-2$ fiducial model
temperature of 600 K lies between the temperatures of Kamenetzky et al.
(2014)’s cool and hot components, supporting the argument that the CO J =
$3-2$ emission is created by a combination of gas temperatures. Also in line
with their model, our model finds the CO J = $6-5$ emission to be dominated by
hot (2,000 K) gas, at the upper bound of their modeled 1$\sigma$ range on gas
temperature. Our model temperatures may be on the higher end due to the much
smaller ALMA beam sizes that probe the nuclear region while Kamenetzky et al.
(2014)’s models incorporate Herschel data that combines emission from both
nuclear and extended regions. The nuclear region is likely to be much warmer
than the extended gas due to star formation and AGN activity concentrated in
the central kpcs. The high gas temperatures are supported by the presence of
prominent H${}_{2}1-0S(1)$ emission that requires temperatures $\textgreater$
1,000 K to be excited (Max et al. (2005); Meijerink et al. (2013)).

The model fit dust temperature range of [20, 100] contains Kamenetzky et al.
(2014)’s best fit dust temperature of 56 K that we used in our continuum mass
calculations (Section 3.1).

Kamenetzky et al. (2014) also models the masses of these cool and warm
components, finding a cool gas mass of 2$\times$10${}^{9}\pm$ 5$\times$108 M⊙
and a warm gas mass of 4$\times$10${}^{8}\pm$ 108 M⊙. Our total modeled mass
for CO J = $3-2$ is 7.6$\times$108 M⊙, higher than Kamenetzky et al. (2014)’s
modeled warm gas mass but lower than their cool gas mass. This is again
consistent with the theory that CO J = $3-2$ is a mixture of cool and warm gas
components. The ranges of CO J=$6-5$ gas mass modeled by LIME spans a factor
of two (Table 3) that contain Kamenetzky et al. (2014)’s masses. However, our
modeled hot gas mass (captured by simulations of CO J = $6-5$) is
6$-$7$\times$108 M⊙, higher than their best-fit warm gas mass by around
50$\%$. Despite the high mass of the hot gas in the model, the total modeled
mass of $\sim$1.4$\times$109 M⊙ (CO J = $3-2$ plus CO J = $6-5$) is close to
the nuclear region’s mass calculated from the continuum in Section 3.1
(1.2$\times$109 M⊙, region 1, Figure 6). This presents a possible conundrum:
the dust continuum observations in Section 3.1 suggest that much of the dust
emission is not captured by these observations, while the modeled gas masses
do appear to capture the majority of the gas mass. It is possible that this
discrepancy comes down to the high $L_{CO}$/$L_{FIR}$ ratio in NGC 6240,
around 10 times higher than other nearby galaxies (Kamenetzky et al., 2014).
The CO emission could be dominated by the shocked gas observed with ALMA,
while there is more dust and cold (less luminous) gas outside of the ALMA
observations that would included in Herschel observations. In other words, due
to the concentrated nuclear excitation of CO in NGC 6240, the ALMA
observations efficiently captured CO emission, but not the extended dust
emission.

###  4.5 Transience and Stability of Fiducial Models

In this section we examine the modeled velocity dispersion and circular
velocities in the context of the model geometry, finding that the model
suggests a transient structure that is stable against collapse.

The ratio of our average rotational velocity to our modeled velocity
dispersion is $\textless v\textgreater$/ $\sigma\sim$ 0.7 for the thin
fiducial models presented in Figure 9. According to Tacconi et al. (1999) this
is in the range that indicates a disk that must be geometrically thick.
However, the fiducial models are quite thin, with $z_{sh}$ of only 20 pc and
$r_{sh}$ of 250$-$390 pc. This aspect ratio is unlikely for a self-gravitating
structure with such high velocity dispersion since the dispersion would have
the effect of “puffing up” the disk. It is also possible that with such a thin
concentration of gas, the modeled high velocity dispersion could partially be
accounted for by shear. That is, the gas concentration could be undergoing
tidal disruption in the plane of the sky.

The thicker fiducial model presented in Table 4 and Figure 10 has a higher
ratio of $\textless v\textgreater$/ $\sigma\sim$ 2.3 and a thicker geometry
($z_{sh}$ = 60 pc, $r_{sh}$ = 270 pc). In both senses this model is less
extreme than the thinner model, though it visibly does not fit the observed
line profiles as well. That the thicker model would correspond to a lower
velocity dispersion may seem counterintuitive, expecting that higher velocity
dispersions would “puff up” the gas. However, if one invokes shear in the
plane of the sky this would artificially inflate the velocity dispersion for
thinner and more optically thin gas concentrations.

We can also conduct a back-of-the-envelope calculation to determine the
lifetime of the gas concentration given its modeled velocity dispersion and
vertical extent. For the thin model, the velocity dispersion is 140 km/s
(assuming no contribution from shear) and the vertical extent of the
concentration is $\sim$ 40 pc (two vertical scale heights). Gas moving at 140
km/s would travel this distance after 300,000 years. This is an extremely
short timescale in the context of galaxy mergers that occur over timescales of
Gyr, suggesting the model represents a highly transient structure. For the
thicker model, gas traveling at a speed equal to the velocity dispersion of
110 km/s would travel two scale heights’ distance (120 pc) in $\sim$ 1 Myr,
still a short timescale in the context of galaxy mergers. Such short lifetimes
are unlikely for a self-gravitating structure.

As a further test of the models in the context of a self-gravitating
structure, we can compare the velocity dispersion and circular velocities to
the escape velocity from the modeled gas concentration. For this calculation,
we use the escape velocity of $\sqrt{2GM/r}$ for the entire gas concentration
with $M=$1.4$\times$109 M⊙. At the radial scale height of the thinner model
$r_{sh}$ = 390 pc, $v_{esc}$ is 170 km/s. This value is similar to the thin
models’ velocity dispersion of 140 km/s, and is less than the circular
velocity (100 km/s) added to the modeled dispersion. At the radial scale
height of the thicker model $r_{sh}$ = 270 pc, $v_{esc}$ is 210 km/s. This
exceeds the thicker model’s turbulent velocity of 110 km/s but is less than
its circular velocity alone ($v_{circ}$ = 250 km/s). This means both models’
gas velocity can exceed the escape velocity of the gas concentration,
indicating it is extremely unlikely this gas is a self-gravitating disk.

Given the fiducial LIME model, we can also check the stability of the gas
concentration by calculating the Toomre parameter $Q$ = $\sigma_{r}\kappa$ /
$\pi$ G $\Sigma_{gas}$. Here $\kappa$ = $\sqrt{3}$ $v_{max}$ / $R$ is the
epicyclic frequency, $\sigma_{r}$ is the line-of-sight velocity dispersion,
and $\Sigma_{gas}$ is the mass surface density of the gas (Toomre, 1964). $Q$
for the thinner LIME model is 2.3 for CO J = $6-5$ and 2.5 for the CO J =
$3-2$ model. A $Q$ $\textgreater$ 1 means the gas is stable against collapse
at this time, with the majority of models within the acceptable parameter
ranges fitting into this category. These high values of $Q$ are consistent
with this highly turbulent system. The majority of the luminosity being
generated in this central mass concentration is then unlikely to be due to
star formation, consistent with other papers that argue for heating from
shocks and superwinds originating in starbursts around the nuclei, not in the
central region between the two nuclei (e.g. Tecza et al. (2000); Max et al.
(2005); Engel et al. (2010)).

###  4.6 Model Residuals: What is Captured by the Model?

Figures 11 and 12 show the moment 0 and 1 normalized residuals (upper panels)
and moment 1 absolute residual (lower panel) between ALMA data and the thinner
fiducial LIME models for CO J = $6-5$ and J = $3-2$ (presented in Table 3 and
Figure 9). The normalized fractional residuals are calculated by subtracting
the model from the observation and dividing by the observation for each pixel.
We present the residuals for the thinner and not the thicker models because
they most closely represent the observed line profiles.

The residuals for each moment take a similar shape for both emission lines.
The low values of the moment 0 fractional residual in the brightest central
region indicate that the test model captures the majority of the line emission
in this central region. However, the model misses the northernmost portion of
the observed gas concentration, especially apparent in the CO J = $6-5$
residual. The test model does not accurately describe the
$\textless$v$\textgreater$ data, as evidenced by the large values and amount
of structure that remains in the fractional residuals of moment 1 for both CO
J = $6-5$ and J = $3-2$. The non-normalized residual of
$\textless$v$\textgreater$ (observed - model), shown in the bottom panel of
each figure, remains largely unchanged from the observed values. This
indicates that the model does not reproduce the velocity structure well at
all, and explains why the normalized residuals are all close to one.

The residuals shed light on what components of NGC 6240’s molecular gas are
captured by the test model, if any. There are aspects of the gas that are
clearly not captured, for example the modeled line profiles in Figure 9 miss a
shoulder of highly redshifted emission ($\textgreater$ 300 km/s) in the
northeast plotted therein as dark brown and black line profiles. This emission
corresponds to the highly redshifted gas visible in the average velocity maps
($\sim$ 300 km/s in Figure 2 and $\sim$ 400 km/s in Figure 3) and in the
channel maps (up to 665 km/s in Figure 4 and up to 724 km/s in Figure 5). It
also corresponds to the highly redshifted gas in Figure 8 with an average
velocity of $\sim$400 km/s when isolated from the quiescent gas, or $\sim$250
km/s when including the quiescent gas. This gas is dim compared to the
majority of the line emission, and as such does not augment the normalized
fractional residual moment 0 maps above $\sim$0 despite not being described by
the models.

Despite missing this high velocity gas, the low values for the integrated line
emission residuals in the central/southern portion of the central gas
concentration indicate that the models capture the majority of the CO emission
in this, the brightest and most massive part of the ALMA observations. We
calculate that the model captures 96$\%$ of the observed emission, a value
found by dividing the sum of the modeled CO J = $6-5$ emission divided by the
sum of the observed emission within the 5$\sigma$ contour area. Capturing the
majority of the emission indicates we accurately model the gas mass and gas
temperature, which are degenerate parameters that each increase total emission
when increased. Constraints on gas mass and temperature are provided by the
modeled line profile shapes with exceedingly high masses or low temperatures
resulting in doubly-peaked profiles. Accurate models of mass and temperature
do not directly imply accurate models of gas distribution. The other
parameters important for this are $r_{sh}$, well constrained by the visible
extent of the gas, and $z_{sh}$. $z_{sh}$ is somewhat difficult to constrain
because the gas is nearly face-on, but thick gas results in doubly peaked
profiles constraining $z_{sh}$ to small values (tens of parsecs). Density is
directly calculated from the gas mass, $r_{sh}$, and $z_{sh}$. Our constraints
on these values and the gas temperature indicate we have successfully modeled
the physical distribution and density of the gas. Our models found that this
central gas region is pancake-like: quite thin (tens of pc) with respect to
its extent (hundreds of pc).

While the fractional residuals of the line emission appear to accurately
capture the distribution of the gas, the velocity residuals tell a different
story. The normalized fractional velocity residuals (upper right panel of each
figure) are comprised almost entirely of values close to one, with absolute
residuals (lower panel of each figure) close to the average velocities of the
observed gas. These large residuals suggest the simple exponentially decaying
velocity profile (like that of a disk) does not accurately capture the
kinematics of the gas. When taken with the arguments in Section 4.4 that the
modeled velocity dispersion indicates a transient structure whose velocities
can exceed the escape velocity of the gas concentration, we suggest that the
nuclear molecular gas emission is dominated by a thin pancake-like gas
concentration that is not rotating like a disk. This conclusion is also
consistent with the findings of Section 4.4, finding our masses and
temperatures fit within current observations but the high velocity dispersion
relative to the circular velocity is unlikely for such a thin geometry if the
gas were to be a self-gravitating disk.

Figure 11: Normalized fractional residuals (color) between CO J = $6-5$ data
and simulation for moment 0 (upper left) and moment 1 (upper right). The lower
panel shows the non-normalized residual in moment 1. Dashed contours
correspond to the moment 0 data, solid contours correspond to the 678 GHz
continuum emission, and crosses show the locations of the two AGN. Beam FWHM
contours are shown in the bottom right of each panel.  Figure 12: Normalized
fractional residuals (color) between CO J = $3-2$ data and simulation for
moment 0 (upper left) and moment 1 (upper right). The lower panel shows the
non-normalized residual in moment 1. Dashed contours correspond to the moment
0 data, solid contours correspond to the 345 GHz continuum emission, and
crosses show the locations of the two AGN. Beam FWHM contours are shown in the
bottom right of each panel.

##  5 Is the Central Molecular Gas a Tidal Bridge?

Figure 13: The average velocity and dispersion of CO $6-5$ (red) and CO 3-2
(black) extracted along the semimajor axis of the central gas region (0∘
solid) and along the axis of the modeled gas (33.7∘ blue dashed).

The models indicate that the gas distribution has a vertical thickness of tens
of parsecs with a horizontal extent of hundreds of parsecs. At the location of
the brightest molecular gas emission the fractional moment 0 residuals between
the observations and the models are close to zero. Capturing the majority of
the emission indicates we accurately model the gas mass and gas temperature,
which are degenerate parameters that each increase total emission when
increased. However, as argued in Section 4.6, the molecular gas kinematics are
not well captured by the disk model as demonstrated by the high fractional and
absolute residual values for $\textless v\textgreater$ (moment 1).
Additionally, our calculations and discussion in Section 4.5 indicate the gas
is unlikely to be a self-gravitating disk, consistent with arguments from
Treister et al. (2020), Engel et al. (2010) and Cicone et al. (2018) mentioned
in the introduction. The contra-indication of a gravitationally stable disk is
true for both the thin and thicker fiducial models. Therefore, we must
consider other possibilities to explain the molecular gas kinematics. The
extended distribution of the model could suggest a tidal bridge between the
two nuclei, stretched by tidal forces along the semi-major axis of the
observed CO emission at PA = 0∘. Other scenarios are possible, such as gas
infall and outflow and multiple disks. However, it would not be surprising for
there to be tidally disruption of dissipative gas in the center of two merging
galaxies that have recently undergone a first pass. Furthermore, NGC 6240 does
not exhibit the clear evidence for multiple disks as in other double-nucleus
galaxies, such as Arp 220 (Wheeler et al., 2020). Given that there is already
evidence for a tidal bridge between the nuclei, we explore this
interpretation.

If the gas is indeed a tidal bridge, its low density means it is unlikely to
be in virial equilibrium and could result in sub-thermal excitation. The
X-factor (the relationship between $N(H_{2})$ and the CO luminosity) only
applies for virialized, denser clouds, and can be radically modified for non-
virialized, thin clouds (Tacconi et al., 2008). Physically this can manifest
as higher than expected CO line luminosity compared to the observed H2 and
dust luminosity as well as offsets between the peak of CO line emission and
the peak of H2 and dust emission (Engel et al., 2010; Tacconi et al., 2008).
Our observations show an offset between CO and dust emission peaks, and
comparison to H${}_{2}(1-0)~{}S(1)$ (Max et al., 2005) in Figure 7 also shows
an offset between H2 and CO emission. This sub-thermal and therefore optically
thinner emission could also help to explain the anomalously high
$L_{CO}/L_{FIR}$ observed in NGC 6240, around 10 times higher than other
nearby galaxies (Kamenetzky et al., 2014). To check for the possibility of
sub-thermal emission we compare the observed density to the critical densities
$n_{crit}$ for the observed lines. The modeled average $H_{2}$ density is
$\sim$ 1$\times$103 cm-3 for both the thicker and thinner fiducial models. The
critical density $n_{crit}$ for CO J = $3-2$ emission of 3.6$\times$104 cm-3
and 2.9$\times$105 cm-3 for CO J = $6-5$. $n_{crit}$ for both transitions is
1-2 orders of magnitude higher than the modeled density, indicating that sub-
thermal emission is likely in this central region.

The presence of a tidal bridge or ribbon between the two nuclei is also
supported by observations of H2 $1-0$ S(1) and S(5) presented in Max et al.
(2005) and recreated in our Figure 7. They observe a ribbon of H2 with a
reverse S shape (their Figure 10) extending between the northern and southern
nuclei and postulate it could be material flowing along a bridge connecting
the two progenitor galaxy nuclei. They compare this geometry to the tidal
bridges predicted by computer simulations (e.g. Barnes & Hernquist (1991); see
also Figs. 9 and 10 in Barnes & Hernquist (1996); Figs. 1 and 4 in Barnes
(2002)). One important note is that these simulations show tidal bridges on
scale of $\sim$ 5-40 kpc, while the projected separation between the nuclei in
NGC 6240 is only $\sim$ 1 kpc. Nonetheless, smaller scale simulations (e.g.
Hopkins et al. (2013)) also show significant mass flowing between galaxy
nuclei on scales $\lesssim$ 5 kpc. While the irregular H2 morphology does not
follow the morphology of CO, this could be explained by an altered CO-to-H2
conversion factor in this nuclear region (Tacconi et al., 2008) causing CO to
be brightened with respect to H2 and dust along the N/S axis as described
above. The differences in the two observations are unlikely to be caused by
differences in sensitivities or beam sizes, as the H2 observations show
offsets in peak brightness located at the southern nucleus, farther than 1.5
beam FWHM away from the peak brightness observed in the CO observations.

Tidal bridges and filaments have been shown to have relatively small line-
widths ($\sim$ 50 – 100 km/s) in tidal dwarf galaxies (Braine et al., 2001).
In larger galaxies the linewidths and velocities are larger, for example the
CO FWHM is $\sim$ 200 km/s with $\textless v\textgreater\sim$ 130 km/s in the
bridge between the Taffy galaxies (Braine et al., 2003). In Arp 194 (total
dynamical mass $\textgreater$ 1011 M⊙), the dispersion is $\sim$ 25-125 km/s
with $|\textless v\textgreater|\textless$ 150 km/s in the bridge connecting
the two galaxies (Zasov et al. (2016), their Figure 3). In Figure 13 we plot
the velocity and dispersion along the major axis of the observed gas (PA = 0∘)
and the major axis of the modeled gas (PA = 33.7∘) to compare to these
observed line widths and velocities. For both extractions, the velocities all
lie below approximately $|$100$|$ km/s, with a structure that moves smoothly
from redshifted to blueshifted, as expected for a tidal bridge. The total
velocity dispersion shows values of 150-180 km/s, slightly higher than the
bridge in Arp 194 but within range of the Taffy galaxy bridge. Within the
context of these observations a tidal bridge remains a plausible explanation
of the internuclear molecular gas.

##  6 Possible Fates of the Molecular Gas

The geometry of the merger that formed NGC 6240 is discussed in detail in
Engel et al. (2010). They propose that the merger has a geometry that tends
towards coplanar/prograde because of the extended tidal tails observed in NGC
6240. Lotz et al. (2008) presents simulations of equal-mass gas-rich mergers
with a variety of geometries, including coplanar prograde encounters like NGC
6240. All coplanar encounters of gas-rich Sbc-type galaxies simulated in Lotz
et al. (2008) have two peaks of star formation, one around 1 Gyr and another
stronger starburst around 2 Gyr after the merger begins. Consistently, Tecza
et al. (2000) argued that we are observing NGC 6240 shortly after the first
encounter triggered an initial starburst, as did Engel et al. (2010) who
argued NGC 6240 “has currently elevated levels of star formation compared to a
quiescent galaxy; and will experience another, likely stronger, peak in star
formation rate in the near future when the galaxies coalesce”. Engel et al.
(2010) also make an important note that NGC 6240 is currently barely below the
ULIRG classification of L${}_{IR}\gtrsim$ 1012 L⊙, and will likely breach that
threshold once the second starburst is triggered, supported by the stronger
second starburst predicted by the simulated merger models in Lotz et al.
(2008).

It is possible that the tidal bridge will fall onto the nuclei prior to the
second pass and final coalescence, possibly triggering the second starburst or
feeding further AGN activity. We calculate the free-fall time of the gas of
mass $m$ onto a nucleus of mass $M$ using

$t_{ff}=\frac{\pi}{2}\frac{R^{3/2}}{\sqrt{2G(M+m)}},$ (4)

where $R$ is the distance between the gas and the nucleus. As this is an
order-of-magnitude calculation, we split the molecular gas mass in half and
assume half will fall to the northern nucleus and half to the southern. We use
the nuclei masses of 1.3$\times$1010 M⊙ for the southern nucleus and
2.5$\times$109 M⊙ for the northern (Engel et al., 2010). The projected radius
from the maximum of the integrated line emission to the northern nucleus is
590 pc and 145 pc to the southern. This means $t_{ff,northern}$ =
4.3$\times$106 years and $t_{ff,southern}$ = 2.5$\times$106 years. That is,
the free-fall time of the gas onto the nuclei is a few Myr while the time
between the first pass and the galaxies’ maximum separation prior to the
second pass is likely to be $\sim$ 400 Myr (Lotz et al., 2008). Therefore, it
is likely that this gas will fall onto the nuclei prior to the next pass,
possibly adding to the current nuclear starburst. One issue with this
interpretation is the high modeled velocity dispersion and small vertical
extent that indicate a transient structure that will dissipate after $\sim$0.3
Myr (see argument at the end of Section 4.4). Therefore, it is possible that
the nuclear concentration of gas will dissipate prior to streaming onto the
nuclei.

##  7 Summary

NGC 6240 presents an interesting test case of a galaxy merger between first
pass and final coalescence, an intermediate and turbulent stage of galaxy
evolution. It provides a detailed example to study what comes just before
galaxies evolve into ULIRGs, as it is just below the classification threshold
of a ULIRG but is likely to evolve into a ULIRG when the second, stronger,
merger-induced starburst is triggered. We presented high-resolution ALMA
observations of CO J = $6-5$ and J = $3-2$, the first observations of the
nuclear region of NGC 6240 in CO J = $6-5$ and the highest resolution
observations to date in CO J = $3-2$. We observe similar morphology to
previous CO observations, notably, a concentration of gas between the two
nuclei that is distinct from the continuum that is itself centered around the
nuclei (Tacconi et al., 1999; Scoville et al., 2014, etc.). We model the
molecular gas density and velocity distributions using LIME and find a thin,
pancake-like distribution of gas whose velocities and velocity dispersions
indicate a transient concentration that is unlikely to be a self-gravitating
disk. The model captures the majority of the gas’ emission but fails to
capture the gas’ kinematics, as demonstrated by the residual maps. We instead
argue that the nuclear region observation is consistent with superposed
emission from a tidal bridge and highly redshifted gas. This work demonstrates
the importance of high-resolution multi-line observations when trying to
disentangle the effects of energetic gas acceleration mechanisms, star
formation, and tidal forces in the central regions of major mergers.

We argue that the majority of the central molecular gas concentration is a
tidal bridge connecting the two nuclei of the progenitor galaxies. Our
fiducial models show that this central gas region is likely quite thin
(vertical scale height of 20 to 60 pc) with respect to its extent (horizontal
scale height of 240 to 500 pc). That this central gas is a bridge connecting
the nuclei is supported by Engel et al. (2010) who argue for an altered CO-
to-H2 conversion factor in the central region that would be exacerbated by a
drawn-out, thin, and diffuse gas bridge. The H2 observations presented in Max
et al. (2005) also support the idea of a tidal bridge, which are in turn
motivated by simulations from Barnes & Hernquist (1991) and Barnes & Hernquist
(1996) that show mergers can result in material flowing along bridges
connecting progenitor galaxies. Kinematic arguments from Treister et al.
(2020) suggest the gas bridge connects the two nuclei with clear streaming
kinematics. Gas that is subject to gravitational torques, such as that in a
tidal bridge in the nuclear region, will likely fall into the nuclear regions
by the time of final coalescence (Souchay et al., 2012). Therefore, the
molecular gas in the central region will likely fall into the nuclear regions
of NGC 6240 and contribute to the second starburst that will turn NGC 6240
into a bona fide ULIRG. However, the high velocity dispersion with respect to
the vertical extent of the gas means it is also possible that the structure
will dissipate prior to this streaming. These observations and models shed
light on one mechanism for star formation in major gas-rich mergers, that is,
small-scale tidal bridges forming between progenitor galaxy nuclei that may
ultimately feed into the nuclear regions.

This paper makes use of the following ALMA data: ADS/JAO.ALMA 2013.1.00813.S.,
2015.1.00658.S. ALMA is a partnership of ESO (representing its member states),
NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA
(Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of
Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The
National Radio Astronomy Observatory is a facility of the National Science
Foundation operated under cooperative agreement by Associated Universities,
Inc.

## References

  * Alonso-Herrero (2013) Alonso-Herrero, A. 2013 arxiv e-prints , arXiv:1302.2033
  * Barnes (2002) Barnes, J. E. 2002, MNRAS, 333, 481
  * Barnes & Hernquist (1991) Barnes, J. E., & Hernquist, L. E. 1991, ApJ, 370, L65
  * Barnes & Hernquist (1996) Barnes, J. E., & Hernquist, L. E. 1996, ApJ, 471, 115
  * Braine et al. (2001) Braine, J., Duc, P.-A., Lisenfeld, U., et al. 2001, A&A 378, 51–69
  * Braine et al. (2003) Braine, J., Davoust, E., Zhu, M., Lisenfeld, U., Motch, C., & Seaquist, E. R. 2003, A&A 408, L13-L16
  * Brinch & Hogerheijde (2010) Brinch, C. & Hogerheijde, M. R. 2010, A&A, 523, A25.
  * Casey (2012) Casey, C. 2012 MNRAS, 425, 3094–3103
  * Cicone et al. (2018) Cicone, C., Severgnini, P., Papadopoulos, P. et al. 2018 ApJ863, 143
  * Consiglio et al. (2017) Consiglio, S.M., et al. 2017 ApJ850, 54
  * Engel et al. (2010) Engel, H., Davies, R. I., Genzel, R., et al. 2010, A&A, 524,A56
  * Fried & Shulz (1983) Fried, J. W., & Schulz, H. 1983, A&A, 118, 166
  * Genzel et al. (1998) Genzel, R., et al. 1998, ApJ, 498, 579
  * Gerssen et al. (2004) Gerssen, J., van der Marel, R. P., Axon, D., et al. 2004, AJ, 127, 75
  * Hagiwara et al. (2011) Hagiwara, Y., Baan, W. A., & Klockner, H. 2011, AJ, 142, 17
  * Hollenbach & Thronson (1987) Hollenbach, D.J., & Thronson, Harley A., eds. of Interstellar Processes, D. Reidel Publishing Company, Dordrecht, 1987
  * Hopkins et al. (2013) Hopkins, P. F., Cox, T. J., Hernquist, L., Narayanan, D., Hayward, C. C., & Murray, N. 2013, MNRAS 430, 1901–1927
  * Iono et al. (2007) Iono, D., Wilson, C. D., Takakuwa, S., et al. 2007, ApJ, 659, 283
  * James et al. (2002) A. James, A., Dunne, L., Eales, S. & Edmunds, M. G. 2002 MNRAS, 335, 753-761
  * Kamenetzky et al. (2014) Kamenetzky, J., Rangwala, N., Glenn, J., Maloney, P. R., & Conley, A. 2014 ApJ795, 174
  * Klaas et al. (2001) Klaas, U., Haas, M., Miller, S. A. H., et al. 2001, A$\&$A, 379, 823
  * Komossa et al. (2003) Komossa, S., Burwitz, V., Hasinger, G., Predehl, P., Kaastra, J., & Ikebe, Y. 2003, ApJ, 582, L15
  * Lonsdale, Farrah & Smith (2006) Lonsdale, C. J., Farrah, D. & Smith, H. E. “Ultraluminous Infrared Galaxies” Astrophysics Update 2, Springer Praxis Books. ISBN 978-3-540-30312-1. Praxis Publishing Ltd, Chichester, UK, 2006, p. 285
  * Lotz et al. (2008) Lotz, J. M., Jonsson, P., Cox, T. J. & Primack, J. R. 2008, MNRAS, 391, 3, 1137–1162
  * Max et al. (2005) Max, C. E., Canalizo, G., Macintosh, B. A., et al. 2005, ApJ, 621, 738
  * McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., Young, W., Golap, K. 2007, Astronomical Data Analysis Software and Systems XVI (ASP Conf. Ser. 376), ed. R. A. Shaw, F. Hill, and D. J. Bell (San Francisco, CA: ASP), 127
  * Meijerink et al. (2013) Meijerink, R., Kristensen, L. E., Weiß, A., et al. 2013, ApJL, 762, L16
  * Müller-Sánchez et al. (2018) Müller-Sánchez, F., Nevin, R., Comerford, J. M., Davies, R. I., Privon, G. C., & Treister, E. 2018 Nature, 556, 345–348
  * Sanders & Mirabel (1996) Sanders, D. B., & Mirabel, I. F. 1996, ARA&A, 34, 749
  * Sanders et al. (1988) Sanders, D. B., Soifer, B. T., Elias, J. H., Madore, B. F., Matthews, K., Neugebauer, G., & Scoville, N. Z. 1988, ApJ, 325, 74
  * Scoville et al. (2014) Scoville, N., Sheth, K., Walter, F., et al. 2015, AJ 800, 70
  * Souchay et al. (2012) Souchay, J., Mathis, S., & Tokieda, T. 2012 “Tides in Astronomy and Astrophysics”, Lecture Notes in Physics 861, DOI: 10.1007/978-3-642-32961-6
  * Tacconi et al. (2008) Tacconi, L. J., Genzel, R., Smail, I., et al. 2008, ApJ, 680, 246
  * Tacconi et al. (1999) Tacconi, L. J., Genzel, R., Tecza, M., Gallimore, J. F., Downes, D., & Scoville, N. Z. 1999, ApJ, 524, 732
  * Tecza et al. (2000) Tecza, M., Genzel, R., Tacconi, L. J., Anders, S., Tacconi-Garman, L. E., & Thatte, N. 2000, ApJ, 537, 178
  * Thronson et al. (1990) Thronson, H. A., Jr., Majewski, S., Descartes, L., & Hereld, M. 1990, ApJ, 364, 456
  * Toomre (1964) Toomre, A. 1964, ApJ, 139, 1217
  * Treister et al. (2020) Treister, E., Messias, H., Privon, G. C., et al. 2020, arXiv:2001.00601v1
  * Vignatti et al. (1999) Vignati, P., et al. 1999, A&A, 349, L57
  * Wang et al. (2014) Wang, J., Nardini, E., Fabbiano, G., et al. 2014, ApJ, 781, 55
  * Wheeler et al. (2020) Wheeler, J., Glenn, J., Rangwala, N., Fyhrie, A. 2020 arxiv e-prints , arXiv:2006.09262
  * Wright et al. (1984) Wright, G. S., Joseph, R. D., & Meikle, W. P. S. 1984, Nature, 309, 430
  * Yun et al. (2001) Yun, M. S., & Hibbard, J. E. 2001, in Gas and Galaxy Evolution, eds. J. E. Hibbard, M. Rupen, & J. H. van Gorkom ASP Conf. Proc., 240, 866
  * Zasov et al. (2016) Zasov, A. V., Saburova, A. S., Egorov, O. V., Afanasiev, V. L. 2016, MNRAS, 462, 4, 3419–3429


More Zero-Free Regions for Fractional Hypergeometric Zeta Functions

Demessie Ergabus Birmechu and Hunduma Legesse Geleta

Addis Ababa University, Department of mathematics

Email<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
Abstract. Some zero free regions were known on the right-half of the complex
plane in the form of vertical strips for fractional hypergeometric zeta
functions. In this paper we describe and demonstrate zero free regions on the
left-half of the complex plane for fractional hypergeometric zeta functions.
The fractional hypergeometric zeta function of order $a$ has no zeros to the
left-half of the complex plane except the trivial zeros on the real axis.
Key Words and Phrases:
Riemann zeta function, Zero free regions, fractional hypergeometric zeta
function

##  1 Introduction

In [2], it was described and demonstrated that zero free region on the right
half-plane $\mathcal{H}=\\{\sigma+it\in\mathbb{C}:\sigma>1\\}$ in the form of
vertical strips $V_{a}=\\{\sigma+it\in\mathbb{C}:1\leq\sigma<2-a\\}$ for
fractional hypergeometric zeta functions $\zeta_{a}(s)$ and
$V_{N}=\\{1\leq\sigma<2\\}$ for hypergeometric zeta functions $\zeta_{N}(s),$
where $``a"$ is a positive real number between $0$ and $1$ and $``N"$ is a
natural number. In [4], it was shown that the fractional hypergeometric zeta
functions $\zeta_{a}(s)$ is zero free for infinitely many positive real number
$a$ in the neighborhood of 1. Specifically there is a positive number $\delta$
so that $\zeta_{a}(s)$ is zero free where the classical Riemann zeta function
$\zeta(s)$ is zero free for $a\in(1-\delta,1+\delta).$

In [3, 5], it was shown that $\zeta_{a}(s)$ has infinitely many simple poles
at $1,~{}0,~{}-1,~{}-2,\ldots,$ and infinitely many zeros at
$1-a,-a,-(1+a),-(2+a),\ldots$ which they called trivial zeros for fractional
hypergeometric zeta function. Observe that if $a>1,$ is fixed, then all the
trivial zeros for fractional hypergeometric zeta function are on the negative
real axis and if $0<a<1,$ then $1-a$ is the only positive trivial zero for
fractional hypergeometric zeta function. The question whether or not these
zeros are the only zeros of $\zeta_{a}(s)$ is unsettled. The authors in [2] in
their conclusion part suggested that it is possible to extend zero-free
regions for these families of zeta functions to both the right and left half
of the complex plane with some evidence but no proof. In this paper motivated
by their suggestions therefore, we show that the fractional hypergeometric
zeta functions have no zeros for $a>1$ on $\\{\sigma+it:\sigma<1\\}$ and for
$0<a<1$ on $\\{\sigma+it:\sigma<0\\}$ except the trivial zeros aforementioned.
For $0<a<1$ at present time we are not sure whether or not $1-a$ is the only
zero on the critical strip of the Riemann zeta function
$\\{\sigma+it:0<\sigma<1\\}$ for the fractional hypergeometric zeta functions.
To prove these results we use the method of analytic continuation for
$\zeta_{a}(s)$ strip-by-strip as shown in [5] to the left-half of the complex
plane, positivity results on oscillatory integrals and monotonicity of real-
valued functions see [1].

Throughout this paper we use the following notations for vertical strips
between the poles of the fractional hypergeometric zeta function. For each
natural number $n$ let us define the vertical strip $V_{n}$ as follows:

$V_{n}=\\{s\in\mathbb{C}:~{}s=\sigma+it,t\in\mathbb{R}~{}\text{and}~{}1-n<\sigma<2-n\\}$
$V_{n}^{+}=\\{s\in\mathbb{C}:~{}s=\sigma+it,t>0~{}\text{and}~{}1-n<\sigma<2-n\\}$
$V_{n}^{-}=\\{s\in\mathbb{C}:~{}s=\sigma+it,t<0~{}\text{and}~{}1-n<\sigma<2-n\\}$
$V_{n}=V_{n}^{+}\cup V_{n}^{-}\cup\mathbb{R}.$

For each natural number $n$ let $F_{n,a}(s)$ be the analytic continuation of
$\displaystyle{\frac{\Gamma(s+a-1)}{\Gamma(a+1)}\zeta_{a}(s)=\int_{0}^{\infty}\frac{x^{s+a-2}}{a\gamma(a,x)e^{x}}dx}$
on $V_{n}.$ Then we observe that

$\Im\left(\frac{\Gamma(s+a-1)}{\Gamma(a+1)}\zeta_{a}(s)\right)=\Im\left(F_{n,a}(s)\right)$

on $V_{n}.$ So to prove that $\zeta_{a}(s)$ has no zeros on $V_{n},$ except
for the trivial zeros mentioned it is enough to show
$\Im\left(F_{n,a}(s)\right)$ has no zeros on $V_{n}^{+}$ and $V_{n}^{-}.$ We
also consider two cases for $a,$ when $0<a<1$ and $a>1.$

The main results of this paper are the following:

######  Theorem 1.1.

For fixed positive number $a\in(0,1)$, $\zeta_{a}(s)$ has no zeros on $V_{n}$
except for infinitely many “trivial” zeros on the left side of $\sigma=0$, one
in each of the intervals $I_{n}=[-n,1-n],$ where $n\in\mathbb{N}.$

######  Theorem 1.2.

For fixed positive real number $a>1,$ $\zeta_{a}(s)$ has no zeros on $V_{n}$
except for infinitely many “trivial” zeros on the left side of $\sigma=0$, one
in each of the intervals $I_{n}=[1-n,2-n],$ where $n\in\mathbb{N}.$

In the first result theorem 1.1 the trivial zero $``1-a"$ is not included, as
we are not sure wheher or not it is the only positive root in the critical
strip of the classical Riemann zeta function. We present this issue on
conclusion part as conjecture based on some evidences.

The structure of the present work is as follows. In section 2 we review some
of the main results obtained so far regarding fractional hypergeometric zeta
functions which we think important for the coming sections. In section 3 we
reveal and prove our main results and demonstrate that $\zeta_{a}(s)$ are
zero-free on left-half of the complex plane except the aforementioned trivial
zeros. In section 4 we give some concluding remarks and discussion.

##  2 Preliminaries

In this section, we review basic terms and results which we will use to prove
the main results of this paper. As our work is the continuation of [2] and [4]
first we review some basic results concerning about fractional hypergeometric
zeta functions, and second analytic continuation, to state and prove our
results in perspective.

######  Definition 2.1.

The Riemann zeta function $\zeta(s)$ is a function of a complex variable
$s=\sigma+it$ defined as an infinite sum given below, convergent for all
$\sigma>1,$

$\zeta(s)=\sum_{n=1}^{\infty}\frac{1}{n^{s}}.$

This is its Dirichlet series representation.

For $\sigma>1,$ the classical zeta function $\zeta(s)$ is also defined by the
integral given below,

$\zeta(s)=\frac{1}{\Gamma(s)}\int_{0}^{\infty}\frac{x^{s-1}}{e^{x}-1}dx,$

where $\Gamma(s)=\int_{0}^{\infty}x^{s-1}e^{-x}dx$ is the gamma function. For
$\sigma>1,$ $\zeta(s)$ is also defined by

$\zeta(s)=\prod_{p\in\mathbb{P}}\left(1-p^{-s}\right)^{-1},$

where the product runs over all prime numbers $\mathbb{P}.$ This product
formula is called Euler’s product formula for the Riemann’s zeta function.
Since a convergent infinite product of non-zero factors is not zero, the zeta
function does not vanish in the right half of complex plane. Therefore, this
product formula is one of the important tools to show that the Riemann zeta
function is zero free in the right half of the complex plane (actually for
$\sigma>1$.)

The functional equation is also another important representation of the zeta
function to locate the zeros of zeta. It is given by

$\zeta(s)=2(2\pi)^{s-1}\zeta(1-s)\Gamma(1-s)\sin\left(\frac{\pi s}{2}\right)$

and this is the celebrated functional equations for the Riemann zeta function.
The reflection principle

$\zeta(\overline{s})=\overline{\zeta(s)}~{}~{}\text{for}~{}~{}s\in\mathbb{C}$

provides a further functional equation for the Riemann zeta function. The
functional equation, together with the reflection principle, evokes a strong
symmetry of the Riemann zeta function with respect to the so called critical
line $\sigma=\frac{1}{2}.$

######  Definition 2.2.

The points $s=-2,-4,-6,\ldots$ are called the “trivial” zeros of the zeta
function $\zeta(s)$ and the strip $\\{s\in\mathbb{C}:0\leq\sigma\leq 1\\}$ is
called the critical strip.

Regarding the zeros inside the critical strip, it is conjectured that these
nontrivial zeros all must be located on the critical line at
$\Re(s)=\frac{1}{2}.$ This conjecture is known as Riemann’s Hypothesis. As a
generalization of the Riemann zeta function $\zeta(s)$ via integral
representation we have the following definition:

######  Definition 2.3.

[5, Definition 2.1]. The Fractional Hypergeometric zeta function
$\zeta_{a}(s)$ is defined and analytic for all positive real number $a$ and
$\sigma>1$ as

$\zeta_{a}(s)=\frac{\Gamma(a+1)}{\Gamma(s+a-1)}\int_{0}^{\infty}\frac{x^{s+a-2}e^{-x}}{a\gamma(a,x)}dx,$

where $\Gamma(s)$ is the Gamma function defined by

$\Gamma(s)=\int_{0}^{\infty}x^{s-1}e^{-x}dx$

and $\gamma(a,x)$ is the lower incomplete Gamma function of the form

$\gamma(a,x)=\int_{0}^{x}e^{-t}t^{a-1}dt.$

Observe that when $a=N,$ natural number we get the Hypergeometric zeta
functions $\zeta_{N}(s)$. In [5] the Fractional Hypergeometric zeta functions
$\zeta_{a}(s)$ has zeros at $s=1-a,-a,-(1+a),-(2+a),\ldots,$ and poles at
$s=1,0,-1,-2,-3,\ldots.$ These zeros are called the trivial zeros of
$\zeta_{a}(s).$

Concerning zero free regions for fractional hypergeometric zeta function
$\zeta_{a}(s),$ the following results were known.

######  Theorem 2.4.

[2, Theorem 6]. Let $0<a<1$ be fixed. Then $\zeta_{a}(s)\neq 0$ in the
vertical strip $V_{a}=\\{s=\sigma+it\in\mathbb{C}:1\leq\sigma<2-a\\}.$

######  Theorem 2.5.

[4, Theorem 4.2] There is a positive number $\delta$ such that the fractional
hypergeometric zeta functions of order $``a",$ $\zeta_{a}(s)$ is zero free for
$a\in(1-\delta,1+\delta).$

## Positivity Properties of Integrals

######  Proposition 2.6.

[1, Proposition 3.1] Let the function $h(r)\geq 0$ on
$(0,\infty),h\in\mathit{\text{L}^{1}_{loc}(0,\infty)},h$ decreasing and
strictly decreasing on some open sub-interval of the form
$\left(\frac{k\pi}{t},\frac{(k+1)\pi}{t}\right)$ for some non-negative integer
$k,$ and $t$ positive real number. Then

$\int_{0}^{\infty}h(r)\sin(tr)dr>0.$

Moreover,

$\int_{0}^{T}h(r)\sin(tr)dr>0$

for any $T>0$ provided, $h$ satisfies the above assumptions with $(0,\infty)$
replaced by $(0,T).$

######  Corollary 2.7.

[1, Corollary 3.3] Let $t>0$ and $\tilde{x}_{t,k}\geq 1$ be such that
$t\ln\tilde{x}_{t,k}=2\pi k,$ for some positive integer $k$ and let $g\geq 0$
on
$[\tilde{x}_{t,k},\infty),~{}g\in\mathit{L^{1}_{loc}(\tilde{x}_{t,k},\infty)}$
such that $x\mapsto xg(x)$ is decreasing on $[\tilde{x}_{t,k},\infty),$
strictly decreasing on
$\left(\ln\frac{j\pi}{t},~{}\ln\frac{(j+1)\pi}{t}\right)$ for some positive
integer $j.$ Then

$\int_{\tilde{x}_{t,k}}^{\infty}g(x)\sin(t\ln x)dx>0.$

Moreover,

$\int_{\tilde{x}_{t,k}}^{T}g(x)\sin(t\ln x)dx>0$

for any $T>\tilde{x}_{t,k},$ whenever $xg(x)$ is decreasing in
$[\tilde{x}_{t,k},~{}T],$ strictly decreasing on
$\left(\ln\frac{j\pi}{t},~{}\ln\frac{(j+1)\pi}{t}\right)$ for some positive
integer $j$ such that $\ln\frac{j\pi}{t}\geq\tilde{x}_{t,k}$ and
$\ln\frac{(j+1)\pi}{t}\leq T.$

## Analytic Continuation

Next we review the analytic continuation. Analytic continuation of an analytic
function is a process of extending the domain of the function to a large
domain.

######  Definition 2.8 (Analytic Continuation).

Let $f$ be analytic in domain $D_{1}$ and $g$ be analytic in domain $D_{2}.$
If $D_{1}\cap D_{2}\neq\emptyset$ and $f(z)=g(z)$ for all $z$ in $D_{1}\cap
D_{2},$ then we call $g$ a direct analytic continuation of $f$ to $D_{2}.$

In [5], the analytic continuation of $\zeta_{a}(s)$ has been shown strip-by-
strip in stages. Concerning the analytic continuation of $\zeta_{a}(s),$ the
following results were also known.

######  Theorem 2.9.

[5, Theorem 3.1] For $0<\Re(s)=\sigma<1,$

$\zeta_{a}(s)=\frac{\Gamma(a+1)}{\Gamma(s+a-1)}\left[\Gamma(s-1)+\int_{0}^{\infty}\left(\frac{1}{a\gamma(a,x)}-\frac{1}{x^{a}}\right)x^{s+a-2}e^{-x}dx\right].$

We denoted in the introduction part the analytic continuation of
$\displaystyle{\frac{\Gamma(s+a-1)}{\Gamma(a+1)}\zeta_{a}(s)}$ by $F_{n,a}(s)$
for each natural number $n$, on $V_{n}.$
For examples, on $V_{1}$ we have,

$F_{1,a}(s)=\Gamma(s-1)+\int_{1}^{\infty}\left[\frac{1}{a\gamma(a,x)}-\frac{1}{x^{a}}\right]x^{s+a-2}e^{-x}dx.$

On $V_{2}$ we have,

$F_{2,a}(s)=\Gamma(s-1)+\frac{a}{1+a}\Gamma(s)+\int_{1}^{\infty}\left[\frac{1}{a\gamma(a,x)}-\frac{1}{x^{a}}-\frac{a}{(1+a)x^{a-1}}\right]x^{s+a-2}e^{-x}dx.$

In general on each $V_{n}$ we have $F_{n,a}(s)$ as an analytic continuation of
$\displaystyle{\frac{\Gamma(s+a-1)}{\Gamma(a+1)}\zeta_{a}(s)}.$

##  3 Extending Zero Free Region to the left side of the complex plane for
$\zeta_{a}(s)$

In this section, we prove our main results and demonstrate that the fractional
hypergeometric zeta function $\zeta_{a}(s)$ has no zeros on the left-half of
the complex plane except the trivial zeros.

######  Theorem 3.1.

For fixed positive number $a\in(0,1)$, $\zeta_{a}(s)$ has no zeros on $V_{n}$
except for infinitely many “trivial” zeros on the left side of $\sigma=0$, one
in each of the intervals $I_{n}=[-n,1-n],$ where $n\in\mathbb{N}.$

###### Proof.

By analytic continuation on each $V_{n}.$ we have,

$\displaystyle\frac{\Gamma(s+a-1)}{\Gamma(a+1)}\zeta_{a}(s)=~{}$
$\displaystyle\int_{0}^{\infty}\frac{x^{s+a-2}e^{-x}}{a\gamma(a,x)}dx$
$\displaystyle=~{}$ $\displaystyle F_{n,a}(s).$

So we have,

$\displaystyle\Im\left(\frac{\Gamma(s+a-1)}{\Gamma(a+1)}\zeta_{a}(s)\right)=~{}$
$\displaystyle\Im\left(\int_{0}^{\infty}\frac{x^{s+a-2}e^{-x}}{a\gamma(a,x)}dx\right)$
$\displaystyle=~{}$ $\displaystyle\Im\left(F_{n,a}(s)\right).$

To apply positivity properties of the oscillatory integral, for
$x\in(0,\infty),$ put

$h(x)=\frac{x^{\sigma+a-2}}{a\gamma(a,x)e^{x}}>0.$

Then since $\sigma<0$ and $a\in(0,1)$ we have $\sigma+a-2<0$ on each $V_{n}.$
Thus, the function $h(x)$ is non-negative,
$\displaystyle{h(x)\in\mathit{L^{1}_{loc}(0,\infty)}}$ and strictly decreasing
on any sub-interval of $(0,\infty).$ Hence on each $V_{n}^{+},$

$\int_{0}^{\infty}\frac{x^{\sigma+a-2}e^{-x}}{a\gamma(a,x)}\sin(t\ln x)dx>0.$

This implies that

$\Im\left(F_{n,a}(s)\right)>0.$

Which in turn implies that,

$\Im\left(\frac{\Gamma(s+a-1)}{\Gamma(a+1)}\zeta_{a}(s)\right)>0.$

Thus, $\zeta_{a}(s)\neq 0,$ on $V_{n}^{+}.$ But, then by reflection principle
$\displaystyle{\zeta_{a}(s)=\overline{\zeta_{a}(\overline{s})}},$ so that
$\zeta_{a}(s)\neq 0,$ on $V_{n}^{-}$ as well. Therefore, $\zeta_{a}(s)\neq 0$
on the left-half of the complex plane except the aforementioned trivial zeros
on the real axis. ∎

######  Theorem 3.2.

For fixed positive real number $a>1,$ $\zeta_{a}(s)$ has no zeros on $V_{n}$
except for infinitely many “trivial” zeros on the left side of $\sigma=0$, one
in each of the intervals $I_{n}=[1-n,2-n],$ where $n\in\mathbb{N}.$

###### Proof.

By analytic continuation on each $V_{n}.$ we have,

$\displaystyle\frac{\Gamma(s+a-1)}{\Gamma(a+1)}\zeta_{a}(s)=~{}$
$\displaystyle\int_{0}^{\infty}\frac{x^{s+a-2}e^{-x}}{a\gamma(a,x)}dx$
$\displaystyle=~{}$ $\displaystyle F_{n,a}(s).$

So we have,

$\displaystyle\Im\left(\frac{\Gamma(s+a-1)}{\Gamma(a+1)}\zeta_{a}(s)\right)=~{}$
$\displaystyle\Im\left(\int_{0}^{\infty}\frac{x^{s+a-2}e^{-x}}{a\gamma(a,x)}dx\right)$
$\displaystyle=~{}$ $\displaystyle\Im\left(F_{n,a}(s)\right).$

To apply positivity properties of the oscillatory integral, for
$x\in(0,\infty),$ put

$h(x)=\frac{x^{\sigma+a-2}}{a\gamma(a,x)e^{x}}>0.$

Then since $\sigma<0$ and $a>1$ we have to consider intervals on which
$\sigma+a-2<0$ on each $V_{n}.$ On such intervals the function $h(x)$ is non-
negative, $\displaystyle{h(x)\in\mathit{L^{1}_{loc}(0,\infty)}}$ and strictly
decreasing on any sub-interval of $(0,\infty).$ Hence on each $V_{n}^{+},$

$\int_{0}^{\infty}\frac{x^{\sigma+a-2}e^{-x}}{a\gamma(a,x)}\sin(t\ln x)dx>0.$

This implies that

$\Im\left(F_{n,a}(s)\right)>0.$

Which in turn implies that,

$\Im\left(\frac{\Gamma(s+a-1)}{\Gamma(a+1)}\zeta_{a}(s)\right)>0.$

Thus, $\zeta_{a}(s)\neq 0,$ on $V_{n}^{+}.$ But, then by reflection principle
$\displaystyle{\zeta_{a}(s)=\overline{\zeta_{a}(\overline{s})}},$ so that
$\zeta_{a}(s)\neq 0,$ on $V_{n}^{-}$ as well. Therefore, $\zeta_{a}(s)\neq 0$
on the left-half of the complex plane except the aforementioned trivial zeros
on the real axis.

∎

##  4 Conclusion

In this paper we described zero free regions on the left-half of the complex
plane. Moreover for $0<a<1$ we showed that there are infinitely many trivial
zeros one in each interval $(1-n,2-n)$ for $n$ natural number. If $a>1,$ then
on each interval $(-n,-(n-1))$ we have exactly one trivial zero. Generally, if
$n<a<(n+1),$ then on each interval $[-n,-(n-1)]$ we have exactly one trivial
zero of fractional hypergeometric zeta function of order $``a".$ But for
$0<a<1,$ we are not sure whether or not $``1-a"$ is the only zero of
$\zeta_{a}(s)$ on $I_{1}=(0,1).$ However, by analyzing the paper by [4] we
conjecture that there are more zeros of $\zeta_{a}(s)$ on $(0,1)$ besides the
trivial zero $1-a$ and we call these zeros the nontrivial zeros of the
fractional hypergeometric zeta functions. We expect and hope that a proof of
such conjecture may shed light either to prove or disprove the existence of
zeros of classical Riemann zeta function in the critical strip.

## References

  * [1] Albeverio, Sergio, and Christof Cebulla. “Müntz formula and zero free regions for the Riemann zeta function.” Bulletin des sciences mathematiques 131(2007), No. 1, 12-38.
  * [2] D. E. Birmechu, & H. L. Geleta, (2022), New Zero-Free Regions for Hypergeometric Zeta and Fractional Hypergeometric Zeta Fun ctions. Journal of Mathematics, 2022.
  * [3] H. L. Geleta, 2014, Fractional Hypergeometric Zeta Function, Ph.D. Thesis, Addis Ababa University, http://etd.aau.edu.et/handle/123456789/8790.
  * [4] H. L. Geleta, (2022), Exploring zero free regions for fractional hypergeometric zeta functions. Indian Journal of Pure and Applied Mathematics, 1-8.
  * [5] H. L. Geleta and A. Hassen, Fractional Hypergeometric zeta function, Ramanujan J. 41(2016), 421-436.