Complex distribution and velocity field of molecular gas in NGC 1316 as revealed by Morita Array of ALMA

TL;DR

This study uses ALMA to map molecular gas in NGC 1316, revealing complex structures and disturbed kinematics indicative of recent external gas accretion and interaction with nuclear jets.

Contribution

First detailed kpc-scale molecular gas distribution and velocity field in NGC 1316 using ALMA, highlighting recent external gas accretion and jet interaction effects.

Findings

Molecular gas forms shell, blob, and clumps in the galaxy center.

Total molecular gas mass is approximately 5.62 x 10^8 solar masses.

Disturbed velocity field suggests recent external gas accretion (<1 Gyr).

Abstract

We present the results of $^{12}$CO($J$=1-0) mosaicing observations of the cD galaxy NGC 1316 at kpc-resolution performed with the Morita Array of the Atacama Large Millimeter/submillimeter Array (ALMA). We reveal the detailed distribution of the molecular gas in the central region for the first time: a shell structure in the northwest, a barely resolved blob in the southeast of the center and some clumps between them. The total molecular gas mass obtained with a standard Milky-Way CO-to-H$_2$ conversion factor is $(5.62 \pm 0.53)\times10^8$ M$_\odot$, which is consistent with previous studies. The disturbed velocity field of the molecular gas suggests that the molecular gas is injected very recently ($<1$ Gyr) if it has an external origin and is in the process of settling into a rotating disk. Assuming that a low-mass gas-rich galaxy has accreted, the gas-to-dust ratio and H$_2$-to-HI ratio are unusually low ($\sim 28$) and high ($\sim 5.6$), respectively. To explain these ratios, additional processes should be taken into accounts such as an effective dust formation and conversion from atomic to molecular gas during the interaction. We also discuss the interaction between the nuclear jet and the molecular gas.

Full text

1 Introduction

Galaxy major/minor merger and Active Galactic Nucleus (AGN) play an important role in galaxy evolution, especially for massive early-type galaxies (ETGs, e.g., Croton et al., 2006; Guo et al., 2011; Oogi & Habe, 2013; Oogi et al., 2016; Rodriguez-Gomez et al., 2016). Mergers are considered to drive an effective gas inflow into the central region of galaxies and may activate star formation and/or AGN there (e.g., Hopkins et al., 2008). Once the AGN is ignited, its glaring radiation and powerful jet are believed to suppress star formation by heating and/or blowing up the cold interstellar medium (ISM, e.g., Croton et al., 2006; Fabian, 2012). Therefore, it is important to investigate cold ISM properties of merging and/or post-merger galaxies with AGN.

NGC 1316 is one of the best targets to study the galaxy merger and the interaction between the AGN jet and ISM. This is a cD galaxy111 NGC 1316 is classified as a lenticular galaxy in RC3 (de Vaucouleurs et al., 1991), but it is confirmed that this galaxy is a typical D-type galaxy, and probably a cD galaxy based on deep and wide optical observation (Schweizer, 1980; Iodice et al., 2017).

located on the outskirt of the Fornax cluster in the south-west direction ($\sim 4$ degree from the center or $\sim 2$ times virial radius, Drinkwater et al., 2001). This galaxy, also known as Fornax A, is the third-nearest radio-bright galaxy after NGC 5128 (Centaurus A) and M 87. It has prominent radio lobes spanning 33 arcmin at the position angle (P.A.) of 110 degree (Ekers et al., 1983; McKinley et al., 2015). The unusually low X-ray luminosity of the AGN of NGC 1316 suggests a declined activity in the last 0.1 Gyr (Iyomoto et al., 1998; Kim & Fabbiano, 2003; Lanz et al., 2010). In the central region, there is an S-shaped nuclear radio jet, which is considered to be bent by interaction with ISM (Geldzahler & Fomalont, 1984). Observations in X-ray and vibration-rotation H2 line suggest an interaction between AGN jet and ISM of NGC 1316 (Kim et al., 1998; Roussel et al., 2007). A peculiar PAH spectrum of this galaxy (an unusually high intensity ratio of the 11.3 $\mu$m to the 7.7 $\mu$m emission lines and an absence of the 3.3 $\mu$m emission) is observed and considered to be due to an effective destruction of smaller PAHs through sputtering in the hot plasma (Iyomoto et al., 1998; Smith et al., 2007; Kaneda et al., 2007).

There is much observational evidence of a rich history of merging events for NGC 1316. Contrary to the central cD galaxy of the Fornax cluster, NGC 1399, NGC 1316 is surrounded by many late-type galaxies and considered to be in an early phase of cD-galaxy formation (Iodice et al., 2017). Deep optical observations revealed complex structures such as ripples and loops in the outer regions, as well as prominent dust extinction in the central region of the galaxy (see Figure 3a): a shell structure in the northwest, a blob in the south-east, and several dust patches between them (Schweizer, 1980, 1981; Grillmair et al., 1999; Carlqvist, 2010; Iodice et al., 2017). The regions with high dust extinction are bright in far-infrared (FIR) wavelengths (Lanz et al., 2010; Galametz et al., 2012, 2014; Duah Asabere et al., 2016). Schweizer (1980) suggested the presence of the “polar” rapidly rotating ionized disk based on optical spectroscopic observations, although the author presented only a picture of 2D spectrum but not a rotation curve. The author also claimed that the central dust patches are associated with this ionized gas disk.

Emission lines associated with rotational transition of molecules are powerful tools to investigate the velocity field of the dust patches which are traced with continuum emission, since dust and molecular gas are generally well-mixed. In NGC 1316, CO emission is detected in the vicinity of the dust shell in the northwest and the dust blob in the southeast with single-dish observations (Sage & Galletta, 1993; Horellou et al., 2001), whereas atomic gas has not been detected down to $10^{8}$ M⊙ (Horellou et al., 2001). Total molecular gas mass is reported to be $(4.0-5.4)\times 10^{8}$ M⊙ (corrected for the choice of distance and CO-to-H2 conversion factor to match that of this study, Wiklind & Henkel, 1989; Sage & Galletta, 1993; Horellou et al., 2001). The single-dish CO velocity varies by $\sim 460$ km s*-1* when moving from 72” northwest to 35” southeast of the center. This variation is smaller than that measured for the ionized gas by Schweizer (1980), which is $\sim 700$ km s*-1* when moving from 54” northwest to 35” southeast. However, the complex dust structure has not been sufficiently resolved with these single-dish observations.

In this study, we present new results of 12CO($J$=1-0) mapping observations of NGC 1316 with the Morita Array (Compact array) of the Atacama Large Millimeter/submillimeter Array (ALMA) as part of ALMA 12CO($J$=1-0) survey toward 64 Fornax galaxies (project code: 2017.1.00129.S, PI: Kana Morokuma-Matsui). We also conducted optical spectroscopic observations with the Low-Resolution Imaging Spectrometer (LRIS) at the Keck Observatory to compare the kinematics of molecular gas and ionized gas, since the previous optical study presented only a picture of 2-D spectrum. The content of this paper is as follows: ALMA and Keck observations and data reductions are described in section 2, the spatial distribution of molecular gas and velocity structures of both molecular and ionized gas are presented in section 3. We discuss the kinematic nature of the molecular gas and the possible cause for the nuclear jet bending are discussed in section 4 and section 5, respectively, and summarize this study in section 6. We adopt the distance to NGC 1316 of $20.8\pm 0.5$ Mpc (Cantiello et al., 2013) throughout the paper and other basic parameters of NGC 1316 are summarized in Table 1. At this distance, 1*′′* corresponds to $\sim 100$ pc.

2 Observations and data analysis

We describe our ALMA observations and data analysis in this section, focusing on NGC 1316. Detailed information of the whole 64 galaxy samples will be presented in the forthcoming overview paper for the survey (Morokuma-Matsui et al. in prep.). Optical spectroscopic observations with Keck/LRIS are also described in this section. Both the radio and optical observations are summarized in Table 2.

2.1 ALMA 7-m and total-power array observations

The 7-m antennae (7M) and 12-m total power (TP) array observations in Band 3 were carried out during cycle 5 as part of the ALMA survey toward 64 Fornax galaxies.

For the 7M observations, from nine to eleven antennae were employed with baselines ranging from 8.9 m to 48.9 m (uv range from 2.5 k$\lambda$ to 12.0 k$\lambda$, angular scale222https://almascience.org/about-alma/alma-basics from $13^{\prime\prime}$ to $74^{\prime\prime}$) during the observing campaign from October 16 to December 21 in 2017. The number of pointings for the NGC 1316 mosaic is ten. The 64 galaxies were divided into two scheduling blocks (SBs) and NGC 1316 is included in the SB named as ${\tt NGC1316\_a\_03\_7M}$. Each SB consists of 33 execution blocks (EBs). The objects J0334-4008 and J0522-3627 were observed as a phase calibrator and a bandpass/flux calibrator, respectively. In order to cover the whole sample of Fornax galaxies with a single correlator setup, the systemic velocity for each source was fixed to 1500 km s*-1*. The 12CO($J$=1-0) emission line (rest frequency $\nu_{\rm rest}$ of 115.271202 GHz) was covered in one of the upper-side band (USB) spectral windows (SPWs) whose bandwidth and resolution are 1875 MHz ($4902.5$ km s*-1*) and 1.128 MHz ($2.952$ km s*-1*), respectively. The other USB SPW (centered on 113.3962 GHz) and two lower-side band SPWs (centered on 103.2712 GHz and 101.3962 GHz) were mainly used for continuum observations.

The TP on-the-fly observations were carried out from May 29 to 31 in 2018 with three or four antennae. The TP mapping area covers the whole area observed with the 7M array. The SB name of the TP observation of NGC 1316 is ${\tt NGC1316\_a\_03\_TP}$, which consists of five EBs. The correlator and SPW setups are the same as that of 7M observations.

2.2 ALMA data analysis with CASA

Data calibration and imaging were conducted with the standard ALMA data analysis package, the Common Astronomy Software Applications (CASA, McMullin et al., 2007; Petry & CASA Development Team, 2012). The absolute flux and gain fluctuations of the 7M data were calibrated with the ALMA Science Pipeline (version of r40896 of Pipeline-CASA51-P2-B) in the CASA 5.1.1 package. The flux accuracy of the ALMA 7M band-3 data is reported to be better than 5 % (ALMA proposer’s guide). The obtained fluxes of the phase calibrator at individual SPWs are consistent within the errors with the values estimated using the other ALMA measurements of the same source on the closest date to our observations (e.g., our measurement of 0.46 Jy on Oct. 16 2017 vs the other measurements of $0.45\pm 0.03$ Jy which is estimated based on the flux measurements at 91.46 GHz, 103.5 GHz, and 343.48 GHz on Oct. 11 2017)333We utilized getALMAFluxForMS of the “Analysis Utilities” (https://casaguides.nrao.edu/index.php/Analysis_Utilities).. The 7M 12CO($J$=1-0) mosaic data cube was generated with TCLEAN task in CASA version 5.4 with options of Briggs weighting with a robust parameter of 0.5, auto-multithresh mask with standard values for 7M data provided in the CASA Guides for auto-masking444https://casaguides.nrao.edu/index.php/Automasking_Guide (sidelobethreshold of 1.25, noisethreshold of 5.0, minbeamfrac of 0.1, lownoisethreshold of 2.0, and negativethreshold of 0.0), and niter of 10000. The 7M and TP data are combined basically on the image plane with FEATHER task in CASA to account for zero spacing information. The achieved synthesized beam was $15^{\prime\prime}.2\times 7^{\prime\prime}.7$ ($1.5~{}{\rm kpc}\times 0.8~{}{\rm kpc}$ at the distance of NGC 1316) with a P.A. of 89∘. The beam size is a few times smaller than the $43^{\prime\prime}$ beam of the 15 m Swedish-ESO Submillimeter Telescope used in the previous CO observations of NGC 1316 (Wiklind & Henkel, 1989; Sage & Galletta, 1993; Horellou et al., 2001). After the combination, we verified that the 7M+TP data matched the spectral profile of the TP data (Figure 1). The achieved rms noise of the box area centered on the galaxy center with a width of 80 arcsec and with a height of 105 arcsec with the final velocity resolution of 9.9 km s*-1* was 12 mJy beam*-1*, after the primary beam correction.

2.3 Moment maps: CO emission search with SoFiA

We searched for CO emission using the Source Finding Application (SoFiA; Serra et al., 2015). We set up SoFiA to:

• •

normalize the cube by the local noise level in the mosaic; the noise varies by a factor of $\sim 5$ as a function of position on the sky (but only by a factor $1.4$ in the region where we detect CO emission) and by a factor $1.6$ along the frequency axis;

• •

convolve the cube with a set of smoothing kernels, and build a detection mask which includes voxels outside the 4 $\sigma$ range in at least one of the convolved cubes; the smoothing kernels are circular Gaussians on the sky and box functions in velocity, and we use all possible combinations of Gaussian FWHM of 0, 3 pixels and box width of 0, 3, 7 channels;

• •

construct individual objects by merging detected voxels with a friends-of-friends algorithm using a merging radius of 2 pixels along RA and Dec axes and 3 channels along the frequency axis;

• •

remove from the detection mask all objects smaller than 3 pixels along RA or Dec axes or 3 channels along the frequency axis;

• •

remove from the detection mask all objects with an integrated S/N below 3;

• •

dilate the detection mask of the remaining objects by at most 2 pixels and 1 channel to include faint emission at the edge of the detected objects.

The moment maps shown in Figure 2 are obtained considering only voxels included in the resulting SoFiA detection mask. In moment-1 and -2 images, we further blank pixels with a CO surface brightness below $0.9$ Jy beam*-1* km s*-1*, which corresponds to 3 $\sigma$ when averaging 6 channels together ($3.0\times 12\times 10^{-3}$ Jy beam${}^{-1}\times 9.9$ km s${}^{-1}\times\sqrt{(}6)\sim 0.9$ Jy beam*-1* km s*-1*). The continuum map in Figure 2 was generated with the whole SPWs except for the 12CO($J$=1-0) and CN($N$=1-0) line-detected channels. Note that the rest frequencies of nine hyperfine components of $N$=1$\rightarrow$0 of CN molecule are covered in the 113 GHz SPW and a channel range for a tentative detection at 112.85 GHz which probably corresponds to CN($N$=1-0, $J$=3/2-1/2, $F$=5/2-3/2) whose rest frequency of $113.499639$ GHz was excluded when generating continuum image555$<1\sigma$ for 9.9 km s*-1* resolution data without smoothing and $\sim 2\sigma$ for 9.9 km s*-1* resolution data with a 7-ch boxcar smoothing..

2.4 Keck/LRIS observations and data reduction

We obtained optical spectra with LRIS (Oke et al., 1995; Rockosi et al., 2010) installed at the 10-m Keck-I telescope on September 16, 2018 (UT, PI: S. Perlmutter). We used the 560 dichroic mirror. The 600/4000 grism and the 400/8500 grating with a grating angle of 22.9 degree were used for the blue ($3100~{}\AA\lesssim\lambda\lesssim 5600~{}\AA$) and red channels ($5500~{}\AA\lesssim\lambda\lesssim 1~{}\mu$m), respectively. With this set-up, multiple major atomic lines can be observed, such as [O ii]${\lambda 3727}$, H$\beta(\lambda 4861)$, [O iii]${\lambda 5007}$, H$\alpha(\lambda 6563)$, and [N ii]$\lambda 6548,6583$. The width of the slit is 1.0 arcsec. These configurations provide spectral resolution of $R\equiv\lambda/\Delta\lambda\sim 600-1,000$, corresponding to $\sim 300-500$ km s*-1*. The pixel samplings are 2.18 Å pixel*-1* and 2.32 Å pixel*-1* in the blue and red channels, respectively. The slit is set to cross the galaxy center with P.A. of 158.56 deg. The exposure time is 480 sec. We adopt 14 arcsec for the gap width between the two CCDs for both the blue and red channels666https://www2.keck.hawaii.edu/realpublic/inst/lris/longslit\_geometry.html. This is a rough value but does not change our conclusion.

The LRIS data are reduced as follows. We first subtract bias images from the raw data and flat-field bias-subtracted data with domeflat spectra. Wavelength is calibrated with sky emission lines for the red channel data. In the blue channel, since the number of sky lines used for wavelength calibration is not large enough, lamp data is used. Note that the telescope altitude when the lamp data is obtained is different from the one during the observation. Recession velocity is calculated in the same definition as that in ALMA data. Then, we add velocity offsets ($\sim-180$ km s*-1*) to the blue channel data so that [N ii] (in the red channel) and [O ii] (in the blue channel) emission line structures are consistent with each other. Uncertainty in absolute wavelength calibration of the data is $\lesssim 2$ Å, roughly corresponding to an absolute value of velocities of $\sim 100$ km s*-1*. Therefore, we discuss only relative velocity differences between the different lines.

3 Total mass, spatial distribution and velocity field of molecular gas

3.1 Total molecular gas mass

The total molecular gas mass is estimated to be $(5.62\pm 0.53)\times 10^{8}$ M⊙ with the standard CO-to-H2 conversion factor of the Milky-Way of $X_{\rm CO}=2.0\times 10^{20}$ cm*-2* (K km s*-1*)-1 (Bolatto et al., 2013). The conversion from Jy beam*-1* to Kelvin is done based on the Rayleigh-Jeans approximation with an equation of $T=1.222\times 10^{6}\frac{I}{\nu^{2}\theta_{maj}\theta_{min}}$, where $T$ is the brightness temperature in Kelvin, $\nu$ is the observing frequency in GHz, $\theta_{maj}$ and $\theta_{min}$ are half-power beam widths along the major and minor axes in arcsec, respectively and $I$ is the brightness in Jy beam*-1*.

The molecular gas mass is calculated by summing the fluxes of the original moment-0 map (without any masks) for a frequency range of 114.497-114.73 GHz (or a channel range of 191-252 ch with velocity resolution of 9.9 km s*-1* or a velocity range of 1406.8-2013.55 km s*-1* with a rest frame frequency of $115.271202$ GHz) enclosed with a box centered on (R.A., Dec.)${}_{\rm ICRS}=(03^{h}22^{m}41^{s}.448,-37^{\circ}12^{\prime}16^{\prime\prime}.2)$ with a width of 120 arcsec and a height of 160 arcsec to cover the whole CO line-detected region. The error for the molecular mass is calculated with a noise map as $\sigma_{rms}\sqrt{npts}({\rm pixel\ scale})^{2}$, where $\sigma_{rms}$ and $npts$ are the rms and the number of pixels of the box used to calculate total flux, respectively. The noise map is generated with CO line-free channels at the lower and higher frequency-sides of the CO channels. Note that the number of channels used to generate the noise map is the same as the number of channels to calculate the moment 0. The helium contribution to mass is accounted for by multiplying with a factor of $1.36$. The obtained molecular gas mass is comparable to the previous single-dish observations ($(4.0-5.4)\times 10^{8}$ M⊙, Wiklind & Henkel, 1989; Sage & Galletta, 1993; Horellou et al., 2001).

3.2 Spatial distribution

For the first time, our ALMA data clearly show the shell structure (hereafter “Shell”, see Figure 2) in the northwest (NW) side and a large concentration (hereafter “Blob” in Figure 2) in the southeast (SE) side. A small clump was detected along the line connecting the NW Shell and the SE Blob (“Clump” in Figure 2). There are also a weak emission at the center and an extended structure in the east-west direction just above the SE Blob (“Extended” in Figure 2). The fluxes and molecular gas mass for each component are summarized in Table 3.2. In Figure 3, we can see that the spatial distribution of molecular gas traced with 12CO($J$=1-0) excellently matches the dust patches visible in the optical (Grillmair et al., 1999) and FIR images (Xilouris et al., 2004; Temi et al., 2005; Lanz et al., 2010; Duah Asabere et al., 2016), as suggested in the previous studies using lower spatial-resolution CO data (Horellou et al., 2001).

At the central position, a weak and narrow CO emission is detected with a peak intensity of $26.8$ mK at $\sim 1720$ km s*-1* in our $9.9$ km s*-1* velocity-resolution data ($4.2\sigma$, Figure 4), which is consistent within the error with the systemic kinematic-LSR velocity in radio definition of $1732\pm 10$ km s*-1*, which is measured with optical spectroscopy (Longhetti et al., 1998). Horellou et al. (2001) reported a very broad emission ($\sim 500$ km s*-1* at the base) at the central position with a $43^{\prime\prime}$ beam. As they stated, the CO line width is broad because their $43^{\prime\prime}$ beam encloses parts of the NW clump and the SE extended structure.

3.3 Velocity field: not only a simple rotating disk?

In this section, we introduce the velocity field of molecular gas in section 3.3.1 and compare with the velocity field of ionized gas obtained with our Keck/LRIS observations in section 3.3.2.

3.3.1 CO velocity field

The high spatial resolution and sensitivity of the ALMA data reveals a very complex velocity structure which cannot be explained only by a “polar” rotating disk with kinematical major axis along P.A. $\sim-20$ deg (the axis along which the CO velocity range is the largest). Figure 2 (b) and (c) show the moment-1 and moment-2 maps, respectively. Our ALMA data clearly show a monotonic trend in velocity along the NW Shell, $\sim 1600$ km s*-1* at the east edge and $\sim 2000$ km s*-1* at the west edge. This is not what we expect for a rotating disk, i.e, spider diagram, where the velocity increases (decreases) when moving towards the kinematical major axis and decreases (increases) when moving away from it for the receding (approaching) side.

The velocity dispersion of the molecular gas is roughly in the range of $20-50$ km s*-1* at kpc-resolution. There are some areas with as high velocity dispersion as $>80$ km s*-1* and a very high value of $\sim 150$ km s*-1* is found at the west edge of the SE Blob. The former area is an overlapping region of the Shell ($\sim 2000$ km s*-1*) and Clump ($\sim 1740$ km s*-1*) components. The latter is due to multiple velocity components around $\sim 1530$ km s*-1*, $\sim 1730$ km s*-1* and $\sim 1814$ km s*-1* (Figure 5, “Offset” of 5 arcsec in the PV cut “F”). Again, this trend in the moment-2 map is not what we expect for a rotating disk. We often see a high velocity dispersion along the kinematical major axis, since there is an abrupt change in the direction of the velocity vector of gas within the beam along the axis.

In Figure 5, the position-velocity (PV) maps are presented. The velocity difference between the NW Shell and the SE Blob ($\sim 400$ km s*-1*, see PV cut “B”, a line connecting the intensity peaks of NW Shell and SE Blob) is consistent with the value reported in previous CO studies (Horellou et al., 2001; Lanz et al., 2010). The followings are new findings with our ALMA data. The PV cut “C” (a line passing through the Shell, Clump, the galaxy center, Extended, and Blob) suggests that some components seem to follow regular rotation but others do not (e.g., “Clump”). This system may be a combination of a (nearly edge-on) rotating disk/ring and components with disturbed kinematics. The relative amount of the disturbed gas to the rotating gas is large, which is contrary to the other radio galaxies with dust or molecular gas (Kotanyi & Ekers, 1979; de Koff et al., 2000; de Ruiter et al., 2002; Ruffa et al., 2019).

In the NW Shell, there are two velocity components ($\sim 1800$ km s*-1* and $\sim 1875$ km s*-1*) near the brightest point (PV cuts along “B” and “C”). In addition, there are at least two velocity components at the east edge of the NW Shell. We can see along the PV cut “A” and “E” that the NW Shell spans both velocities below and above systemic velocity, which is also one of the features unlikely for a rotating disk.

The SE Blob was barely resolved with our beam and likely to have velocity gradient of $\sim 125$ km s*-1* along northeast to southwest direction (PV cut “F”). In the SE Blob, there are multiple velocity components, the main body around $1500$ km s*-1*, and very weak components around $1650$ km s*-1*, $1725$ km s*-1*, and $1825$ km s*-1* (“Offset” of around $-35$ arcsec in the PV cut plot along “C” in Figure 5). With the sensitivity of our data, it seems that the $1725$ km s*-1* and $1825$ km s*-1* components are discrete structures but the $1650$ km s*-1* component is physically connected with the main body of the SE Blob.

The “Extended” component above the SE Blob has a continuous velocity gradient, $\sim 1475$ km s*-1* at the east-side edge and $\sim 1775$ km s*-1* at the west-side edge (PV cut “D”). The “Clump” seems to have a velocity gradient in the north-south direction (PV cut “C” and Figure 2b). These CO components revealed with ALMA likely correspond to the dust patches recognized in ultraviolet, optical (Figure 3a or Figure 5 of Iodice et al., 2017) and FIR data (Figure 3e or Figures 2, 5 of Lanz et al. 2010; Figure 5 of Duah Asabere et al. 2016).

3.3.2 Comparison between CO and ionized gas

We found that [O ii], [O iii] and [N ii] are detected in emission and Ca ii HK are detected in absorption. Figure 6 shows H$\alpha$ and [N ii] spectra where H$\alpha$ emission seems to be detected. H$\alpha$ is detected but very weak maybe due to it is a combination of the emission of ionized gas and stellar absorption. Therefore, Figure 3c has been supposed to be an H$\alpha$+[N ii] image, but the contribution from [N ii] is likely to be more dominant than H$\alpha$.

The PV maps of [O ii], Ca ii H, [O iii], and [N ii] are presented in Figure 7. Note that the intensity of each line in this figure has not been calibrated. The Ca ii H absorption line is not resolved in our data with velocity resolution of $300-500$ km s*-1* and its velocity is almost constant around the galaxy center. On the other hand, the other lines show a velocity gradient, $\sim 1500$ km s*-1* at the offset of 25 arcsec from the galaxy center (SE) and $\sim 1900$ km s*-1* at $-50$ arcsec offset (NW). Although the velocity resolution is as coarse as $300-500$ km s*-1*, the accuracy of determining the central velocity at each pixel is estimated to be $\lesssim 100$ km s*-1*. Note that Schweizer (1980) reported a larger velocity gradient of $\sim 700$ km s*-1* but with a different slit angle (P.A. of 142∘). There also exist some irregularities which cannot be explained only by rotation, such as a velocity depression at $-15$ arcsec.

Our optical data showed that the kinematics of the ionized gas cannot be explained only by rotation, as we found in the molecular gas data, or the kinematic major axis is largely different from the slit angle of $158.56^{\circ}$. Though the velocity resolution of the ALMA and Keck data are very different ($\sim 10$ km s*-1* vs $\sim 100$ km s*-1*), the velocity gradient of the molecular gas at the same P.A. is roughly consistent with the one seen in ionized gas (Figure 7). The systemic velocity of $1732\pm 10$ km s*-1* recorded in the NED (Table 1) is derived from the central [O ii] velocity (Longhetti et al., 1998) and consistent with our measurement. However, there is a significant offset of a few tens km s*-1* in the velocity of the stellar absorption Ca ii and this is considered to be a real systemic velocity of NGC 1316. It is important to obtain integral field unit (IFU) data, compare the velocity fields of ionized gas and molecular gas in more detail, and also investigate the excitation mechanism of the ionized gas.

4 Implications for the kinematic nature of the molecular gas

In this section, we discuss the kinematic nature of the molecular gas observed in the central region of NGC 1316: inflow (section 4.2) or outflow (section 4.3).

4.1 Relative line-of-sight location of each component

Once the relative line-of-sight location of the NW Shell, the SE Blob, and the nucleus (“Nucleus”, marked with a cross in Figure 2) is determined (Figure 8), more detailed kinematic information of the molecular gas can be derived. Our ALMA observations revealed that larger recession velocity of the NW Shell ($\sim 1850$ km s*-1*) and smaller recession velocity of the SE Blob ($\sim 1500$ km s*-1*) than the systemic velocity of NGC 1316 (1732 km s*-1*, Longhetti et al., 1998). Assuming gas in this region is not rotating, the case of Shell-Nucleus-Blob order from the observer suggests inflow motion (case A) the Blob-Nucleus-Shell case suggests outflow motion (case B) and the cases where both the NW Shell and SE Blob are located on the same side with respect to Nucleus suggest a coexistence of inflow and outflow motions (case C).

Our ALMA data and the ancillary data seems to favor the case A. In Figure 3, both the CO and 160 $\mu$m dust emission are stronger at the peak of SE Blob than that of NW Shell (panel (e)), whereas the dust extinction is stronger at the NW shell than the SE Blob (panel (a)) and the distribution of stellar component is smooth (panel (b)). This suggests that there is a larger amount of ISM at the SE Blob than at the NW Shell whereas the extinction by the SE Blob is weaker than the NW Shell. However, $160\mu$m brightness does not only reflect the amount of dust but also the dust temperature and the strength of the back ground far-ultraviolet radiation field, i.e., heating source of the dust. Therefore we consider the case A most plausible, while we leave the possibility of the other cases in the following discussions. The inclination of the nuclear jet would provide us further hints on the relative line-of-sight location of the Shell, the Blob, and the Nucleus but it is difficult to estimate it with our current data set.

4.2 Inflow?

The previous studies on NGC 1316 consider that the dust and molecular gas of the galaxy is injected via wet minor mergers (e.g., Sage & Galletta, 1993; Horellou et al., 2001; Lanz et al., 2010), since ETGs are generally gas-poor. Generally, the gas-rich ETGs tend to show disturbed morphology in deep optical images (van Dokkum, 2005; Duc et al., 2015), suggesting recent merging events as observed in NGC 1316. Lanz et al. (2010) estimated the dust mass of NGC 1316 to be $2.0\times 10^{7}$ M⊙ (corrected for the different choice of the distance) from $24$ $\mu$m, $70$ $\mu$m, and $160$ $\mu$m data taken with the Multiband Imaging Photometer (MIPS) on the Spitzer space telescope. They showed that the dust-to-stellar mass ratio of NGC 1316 is $>100$ times larger than the value expected from the empirical relation between the dust-to-stellar mass ratio and stellar mass of elliptical galaxies and concluded that the dust currently present in NGC 1316 was injected by a merger galaxy. They also claimed that the spatial coincidence between dust and CO suggests a common origin for them.

In the case of A or C as described above (Figure 8), the molecular gas is likely (at least partially) inflowing. Especially, the multi-wavelength data of NGC 1316 seem to favor the case A. In this section, we discuss the infalling galaxy (section 4.2.1) and its merging timescale (section 4.2.2).

4.2.1 Infalling galaxy (galaxies)

For NGC 1316, Lanz et al. (2010) estimated the stellar mass of the merged galaxy to be $(1-6)\times 10^{10}$ M⊙ which injected the dust and gas to the central region based on the dust mass and the dust-to-stellar mass ratios of Sa-Sm galaxies (assuming that the merged galaxy is an Sa-Sm galaxy).

Stellar masses of infalling galaxies can be also estimated with a gas-to-dust ratio (GDR) and an H2/H i ratio combined with an empirical relations between stellar mass and GDR and H2/H i ratio. Davis et al. (2015) found a high GDR of $2.0-730$ (the median of 315) and a low H2/H i ratio of $0.1-2.0$ (the median of 0.4) in ETGs with prominent dust lanes (corrected for the use of different $\alpha_{\rm CO}$), suggesting a recent merger with a lower mass companion (predicted stellar mass ratio of $40:1$). Note that the stellar masses of Davis et al. (2015) sample are $10^{10.5-11.2}$ M⊙ (the median of $10^{10.8}$ M⊙), which is smaller than that of NGC 1316 of $10^{11.8}$ M⊙. One of the nearby ETGs with a prominent dust lane, NGC 5128 (Centaurus A) whose stellar mass of $1.1\times 10^{11}$, has a GDR of $103$ (Parkin et al., 2012) and an H2/H i ratio of $1.06$ (Charmandaris et al., 2000).

For NGC 1316, given the literature data of dust mass ($2.0\times 10^{7}$ M⊙, Draine et al., 2007) and the upper limit of H i mass ($<10^{8}$ M⊙, Horellou et al., 2001), the obtained GDR and H2/H i ratio of the dust patches are respectively $\sim 28$ and $\gtrsim 5.6$, which are roughly $\sim 10$ times lower and higher than the median values found in the ETGs with dust lane. The inferred stellar mass of accreted galaxy to NGC 1316 using the empirical relation of these ratios and stellar mass is $>10^{11.5}$ M⊙ (Rémy-Ruyer et al., 2014; Bothwell et al., 2014; Davis et al., 2015). Considering an ordered stellar distribution of the main body of NGC 1316 (Schweizer, 1980; Iodice et al., 2017), it is unlikely that NGC 1316 is now experiencing a major merger. In the following paragraphs, we discuss three possible explanations on the observed low GDR and high H2/H i ratio of NGC 1316.

First, it is possible that the atomic gas of the infalling galaxy had been selectively stripped by ram pressure from hot halo of NGC 1316. Cortese et al. (2016) found that H i-deficient cluster galaxies are poorer in atomic but richer in molecular hydrogen if normalized to their dust content. This trend can be considered as a consequence of the selective stripping of the components distributed at the outer radius of galaxies by ram pressure of the inter-galactic hot gas. Cortese et al. (2016) showed that the low GDR and high H2/H i ratio are reproduced with a numerical simulations of ram-pressure gas stripping with models of Bekki (2013, 2014b, 2014a). If we assume that only atomic gas of the infalling galaxy was selectively and completely stripped by ram pressure from hot halo of NGC 1316 and H2/H i ratio of the progenitor was 0.1 (0.01) as dwarf galaxies, the GDR and H2/H i ratio of the infalling galaxy are used to be $\sim 280$ ($\sim 2800$) and $\sim 0.56$ ($\sim 0.056$), respectively. These ratios correspond to galaxies with stellar masses of $\sim 10^{9}$ M⊙ ($<10^{7}$ M⊙) and $\sim 10^{10.5}$ M⊙ ($\sim 10^{9}$ M⊙), respectively. However, it should be noted that NGC 1316 is located at far from the center of the Fornax cluster and its own hot halo is not significantly prominent (Feigelson et al., 1995; Kaneda et al., 1995; Tashiro et al., 2001; Isobe et al., 2006; Tashiro et al., 2009; Seta et al., 2013). In addition, NGC 1316 is considered to be now infalling onto the Fornax cluster for the first time (Drinkwater et al., 2001). Therefore, the ram pressure stripping alone may not be able to explain the observed low GDR and high H2/H i ratio of NGC 1316.

Second, it is claimed that the dust formation and conversion from atomic gas to molecular gas are enhanced in galaxy interactions (Young & Knezek, 1989; Sanders & Mirabel, 1996; Nakanishi et al., 2006; Kaneko et al., 2017). Sage & Galletta (1993) discussed that the high H2/H i ratio of their sample ETGs including NGC 1316 may be due to interaction or infall events. Conversion from atomic to molecular gas occurs primarily on the surface of dust grains (Hollenbach & Salpeter, 1971; Cazaux & Tielens, 2004) down to metallicities as low as $\sim 10^{-5}$ of solar value (Omukai et al., 2010). Once the starburst is triggered due to the galaxy interaction, the dust is formed and blasted into interstellar space via core-collapse supernovae within the lifetime of massive stars (e.g., Cernuschi et al., 1967; Hoyle & Wickramasinghe, 1970; Gall et al., 2014). In this case, the timescale of the conversion, $t_{\rm chem}$ is calculated as $t_{\rm chem}=\frac{1}{n_{\rm H}R_{\rm H_{2}}}$, where $n_{\rm H}({\rm cm^{-3}})$ is the number density of H atoms, and $R_{\rm H_{2}}\ ({\rm cm^{3}\ s^{-1}})$ is the rate constant of H2 formation. $R_{\rm H_{2}}$ is reported to be roughly in the range of $10^{-17}-10^{-16}\ {\rm cm^{3}\ s^{-1}}$ (Jura, 1975; Gry et al., 2002; Browning et al., 2003; Habart et al., 2003). If we adopt a typical volume density of cold neutral medium of $\sim 10~{}{\rm cm^{-3}}$ as $n_{\rm H}$, $t_{\rm chem}$ becomes $\sim 3\times 10^{7}$ yr. Considering that the density should increase during the galaxy merger, this is an upper limit. This is smaller than the dynamical time of the system ($3\times 10^{8}$ yr, see Section 4.2.2) so this process is likely to be occurred.

Third possibility is the ionization of atomic gas by AGN-related radiation (Horellou et al., 2001) or a shock produced by jet. It is known that there is a currently low-activity AGN with low-ionization nuclear emission-line region (LINER) with a low-power nuclear radio jet in the central region of NGC 1316. Based on the comparison between NGC 5128 and NGC 1316, Horellou et al. (2001) discussed that the high H2/H i ratio may be due to stronger effects from the nuclear activity in NGC 1316, e.g., brighter in X-ray for NGC 1316 than NGC 5128 and a perpendicular jet to the dust lane in NGC 5128 (Espada et al., 2009) but not in NGC 1316. The ionization timescale is generally much shorter than the dynamical time of galaxies. The fractional ionization is determined by ionization rate and recombination rate. The recombination timescale ($t_{\rm r}$) is calculated as $t_{\rm r}=\frac{1}{n_{\rm e}\alpha_{\rm A}}$ (where $n_{\rm e}$ is the volume electron density and $\alpha_{\rm A}$ is the total recombination coefficient, Avrett & Loeser, 1988). A typical $t_{\rm r}$ is $\sim 10^{5}/n_{\rm e}$ yr. Considering a typical $n_{\rm e}$ value of $10^{8}{\rm cm}^{-3}$ in “broad-line region (BLR)”, $10^{3}-10^{6}{\rm cm}^{-3}$ in “narrow-line region (NLR)”, $t_{\rm r}$ is much smaller than the galactic timescales ($3\times 10^{8}$ yr, see Section 4.2.2). In case of the cocoon in Cygnus A, $n_{\rm e}\sim 10^{-3}-10^{-1}{\rm cm}^{-3}$ (Snios et al., 2018) and $t_{\rm r}$ becomes $10^{6-8}$ yr. However, it is unlikely that this is the only case, since it is unclear why the ionization of atomic gas is enhanced while the dissociation of molecular gas to atomic gas is not enhanced. It may be possible that the dust shielding works effectively at the area where the molecular gas exist (large $A_{V}$). In addition, if the most atomic gas is ionized, strong H$\alpha$ emission is expected to be observed, but no prominent H$\alpha$ emission was detected in our Keck/LRIS observations, as previous optical study suggested (Schweizer, 1980). It is necessary to observationally investigate the properties of the photo-dissociation region (PDR) of the central region of NGC 1316.

In summary, the second possibility (effective dust formation and conversion from atomic to molecular gas) seems to be a most plausible reason for the observed low GDR and high H2/H i ratio of NGC 1316, and the other two mechanisms (ram-pressure stripping and ionization of atomic gas) may partially contribute to these unusual ratios.

4.2.2 Timescales

Lauer et al. (1995) considered NGC 1316 is in the earliest stage of the “settling sequence” of dust in the ETGs (Tran et al., 2001). Goudfrooij et al. (2001b, a) estimated merger age to be $\approx 3$ Gyr by measuring the age distribution of globular clusters in NGC 1316 with Hubble Spacec Telescope (HST). Sesto et al. (2016, 2018) also measured the age of star clusters and found a younger cluster population ($\sim 1-2$ Gyr) in addition to the intermediate-age clusters ($\sim 5$ Gyr) which is the dominant population. Lanz et al. (2010) estimated a lower limit on the merger age to be 22 Myr by calculating the free-fall time of the central molecular gas to the galaxy center.

High-resolution ALMA data of NGC 1316 revealed a complex spatial distribution and velocity field of molecular gas, suggesting that the molecular gas has not been settled into a steady state yet, which has been already claimed in the previous dust extinction observations (Lauer et al., 1995). Theoretically predicted timescale for gas to settle into a disk, $t_{\rm relax}$ ranges from 108 yr to Hubble time depending on models, and a typical value is $10^{9}$ yr (e.g., Gunn, 1979; Tubbs, 1980; Tohline et al., 1982; Steiman-Cameron & Durisen, 1988; Habe & Ikeuchi, 1985; Christodoulou et al., 1992; Christodoulou & Tohline, 1993). $t_{\rm relax}$ in theoretical models depends mainly on four factors (Tohline et al., 1982; Christodoulou et al., 1992; West, 1994): (1) the degree of misalignment between the angular momentum vector of the gas disk and the symmetry axis of the galaxy; (2) the shape of the galactic potential; (3) the viscosity of the gas and the efficiency of the dissipative force, and (4) the distance from the galactic center. $t_{\rm relax}$ is expected to be longer for larger misalignment for (1), more spherically symmetric potential for (2), smaller viscosity for (3), and larger distance for (4).

Lake & Norman (1983) investigated the star (collisionless component) and gas (dissipative component) orbits in triaxial potential and provided a formula to estimate $t_{\rm relax}$ with a dynamical time, $t_{\rm dyn}$ and an eccentricity of the potential, $\epsilon$ as $t_{\rm relax}=t_{\rm dyn}/\epsilon$. They considered that the accreted and disrupted gas is organized into a sequence of tube orbits, and the tube orbits settled to the equatorial plane via a differential precession damped by a dissipative nature of gas. This process is expected to occur on the timescale of the differential precession of the orbits. For NGC 1316, $t_{\rm relax}$ is estimated to be $8\times 10^{8}$ yr with $t_{\rm dyn}$ of $3\times 10^{8}$ yr ($\sim 100$ km s*-1* at $5$ kpc, Bosma et al., 1985; D’Onofrio et al., 1995; Arnaboldi et al., 1998; Bedregal et al., 2006) and $\epsilon$ of 0.4 (Fig. 7 of Iodice et al., 2017). Here, we assume that the eccentricity of the potential is the same as that of stellar component. Note that the timescale is an upper limit since gas viscosity is not taken into account. However, the effect of the viscosity on $t_{\rm relax}$ is claimed to be as small as $t_{\rm relax}\propto$ (viscosity coefficient)-1/3 (Steiman-Cameron & Durisen, 1988). Varnas (1990) showed that $t_{\rm relax}$ is changed by a factor of $\sim 2$ when changing the viscosity parameter by a factor of 10. Therefore, our data suggest that the central molecular gas is injected within $\sim 1$ Gyr.

The $t_{\rm relax}$ of $\sim 1$ Gyr is also supported by a gas depletion time by star formation, $t_{\rm dep}$ of NGC 1316. $t_{\rm dep}$ is defined as a total cold gas mass divided by star formation rate (SFR). If the $t_{\rm dep}$ is much shorter than the $t_{\rm relax}$, it is unlikely that the molecular gas survive for $t_{\rm relax}$. SFR of NGC 1316 is estimated to be $0.30-0.77$ M⊙ yr*-1* with total infrared luminosity that an AGN contribution is not considered (Duah Asabere et al., 2016). Thus, $t_{\rm dep}$ of NGC 1316 is calculated to be $\sim 1$ Gyr or longer, which is comparable to local late-type galaxies (e.g., Kennicutt, 1998; Bigiel et al., 2011; Kennicutt & Evans, 2012). The star formation properties of NGC 1316 will be investigated in the forthcoming paper (Morokuma-Matsui et al. in prep.).

However, it is possible that the $t_{\rm relax}$ for NGC 1316 is longer than $\sim 1$ Gyr. Davis & Bureau (2016) claimed that a longer $t_{\rm relax}\sim 100~{}t_{\rm dyn}\approx 1-5$ Gyr is required to explain the observed histogram of the difference between the projected angular momenta of stellar and molecular gas of ETGs. Numerical simulations showed that $t_{\rm relax}$ becomes as long as Hubble time under the spherically symmetric potential (Christodoulou et al., 1992). We assumed here that the eccentricity of potential is the same as that of stellar component but it is non-trivial that they are the same. In addition, the true SFR (without AGN contribution) may be much smaller than the adopted value for NGC 1316. In order to constrain the timescale of the merger event for NGC 1316, it is needed to conduct numerical simulations with dust and molecular gas formation as well as AGN feedbacks.

4.3 Outflow?

There is a number of galaxies with AGN (jet)-driven outflows of molecular gas (e.g., Alatalo et al., 2011; Combes et al., 2013; Cicone et al., 2014; García-Burillo et al., 2014). In the cases of B or C (Figure 8), the motion of molecular gas of NGC 1316 can be (partially) explained by outflow. As previous studies with dust data claimed (Geldzahler & Fomalont, 1984; Lanz et al., 2010), the nuclear jet seems to bend at just south of the NW Shell and just north of the SE Blob (Figure 3f), implying a interaction between ISM and nuclear jet (see Section 5). In addition, we showed that the observed GDR and H2-to-H i ratio of NGC 1316 are comparable to those of galaxies with stellar mass of $>10^{11.5}$ M⊙, which is comparable to the stellar mass of NGC 1316. The molecular gas fraction, $f_{\rm mol}=\frac{M_{\rm mol}}{M_{\rm mol}+M_{\star}}$ (where $M_{\star}$ is stellar mass) of NGC 1316 is $\sim 0.6$ %, which is not so high compared to those of ETGs with similar stellar masses (Young et al., 2014). The smaller velocity gradient (outflow velocity) of molecular gas in NGC 1316 ($\sim|200|$ km s*-1*) compared to the typical galaxies with molecular outflow ($>|500|$ km s*-1*) cannot be a strong objection against the internal-origin scenario, because it is possible that the nuclear jet has a small inclination ($\lesssim 24$ degree). The small inclination of the nuclear jet is suggested from the symmetric structure and the intensity ratio between the north and in the south. Therefore, these properties may prefer outflow scenario for molecular gas in NGC 1316.

However, there are also some facts against the outflow scenario for NGC 1316. First, the simplest speculations on the relative line-of-sight locations of the NW Shell, the SE Blob and the Nucleus from the dust-extinction and the CO or dust emission data is the case A, suggesting inflow motion, as shown in Figure 8. Second, the dust-to-stellar mass ratio of NGC 1316 is unusually high (i.e., dust rich) compared to the other early-type galaxies with similar stellar masses (Lanz et al., 2010). Finally, the steep velocity gradient along the NW Shell (PV cut “A” in Figure 5) is difficult to reproduce in the jet-induced molecular outflow.

In reality, the situation is not so simple and all the factors discussed in sections 4.2 and 4.3 may play a role to make the ISM in the current physical status. The gas may be injected by other galaxies and blown out by AGN jet, or the internal origin gas is first blown out by AGN jet then flows back to the nucleus. In addition, it should be noted that the flux ratio of CO to Herschel/PACS $160\mu$m data ($\sim$ GDR) at the NW Shell is roughly two times higher than that of the SE Blob. This may suggest that the degree of jet-ISM interaction is different for the NW Shell and the SE Blob or that it is not a single galaxy which injected the molecular gas and dust. Again, with these observational boundary conditions, it is important to investigate the detailed history of galaxy merger of NGC 1316 with numerical simulations of galaxies with molecular gas and dust formation and AGN feedbacks.

5 Bending of the nuclear jet

Jet bending has been observed in some radio galaxies (e.g., Ekers et al., 1978; Smith & Norman, 1981; Wilson & Ulvestad, 1982). Figure 3(f) showed that the nuclear jet bends at the vicinity where the dust and molecular gas emission are strong, regardless of the origin of the molecular gas of NGC 1316. Three scenarios have been claimed to explain the bending jets: (1) precession of the jet nozzle (e.g., Ekers et al., 1978), (2) buoyant effect (pressure gradient of the ISM, e.g., Smith & Norman, 1981; Henriksen et al., 1981; Bridle et al., 1981) or (3) ram pressure effect from ISM (e.g., Wilson & Ulvestad, 1982). Geldzahler & Fomalont (1984) considers that (1) and (3) are unlikely: for scenario (1), the symmetry of the nuclear jet is not correct for the precession and the orientation of the nuclear jet appears to be associated with kpc-scaled phenomenon (e.g., the X-ray cavity Lanz et al., 2010); for scenario (3), the direction of the rotation of the gas disk implied from the ionized gas and the dust extinction data is opposite to the one expected in the case of ram pressure.

Our data clearly showed that the jet bends at the denser regions (higher pressure), toward sparser regions (lower pressure) of molecular gas, which supports the (2) buoyant effect scenario. Whether or not the (3) ram pressure scenario is ruled out depends on the kinematics of molecular gas. In all cases (A, B, and C), the direction of the jet bending is consistent with the flow direction of the molecular gas (from east to west for the NW Shell and vice versa for the SE Blob). Thus, our ALMA data support that both scenarios (2) and (3).

6 Summary

We conducted a 12CO($J$=1-0) mapping observation of NGC 1316 as part of the ALMA survey of 64 Fornax galaxies and also optical spectroscopic observation with Keck/LRIS. The obtained results and implications based on the data are as follows:

• •

The obtained total molecular gas mass ($5.6\times 10^{8}$ M⊙ with the standard Milky-Way $X_{\rm CO}$) and the velocity gradient between the NW Shell and SE Blob ($\sim 400$ km s*-1*) are consistent with previous single-dish observations within the error (section 3 and Figure 5).

• •

12CO($J$=1-0) distribution coincides with dust structures seen in the previous high-resolution optical and FIR observations (Figure 3). For the first time, our high spatial resolution and deep CO data reveal the NW shell structure, some clumps between the NW Shell and the SE Blob, and the SE extended structure just above the Blob (section 3.2).

• •

We found some features disfavoring a simple rotating disk scenario for the central cold ISM of NGC 1316, such as the crossing of the systemic velocity of NGC 1316 in the NW Shell, the multiple velocity components in the NW Shell (PV cut “A” in Figure 5), the monotonic velocity variation in the NW Shell and the SE Blob (PV cuts “A” and “F”), and the velocity deviations from the interpolation between the velocities of the NW Shell and the SE Blob (PV cut “C”, section 3.3).

• •

Keck/LRIS observation confirms that H$\alpha$ emission is very weak in NGC 1316, as the previous optical study indicates.

• •

Despite some disagreements in kinematics of molecular and ionized gas, the velocity of ionized gas estimated with [O ii] line is roughly consistent with that of molecular gas: $\sim 1500$ km s*-1* at 25 arcsec in SE and $\sim 1900$ km s*-1* at 50 arcsec in NW from the galaxy center. Although the accuracy of determining the central velocity of ionized gas is not so high ($\sim 100$ km s*-1*), the obtained PV map of [O ii] cannot be explained only by rotation.

• •

If the observed molecular gas has an external origin, the complex spatial distribution and velocity field suggest a recent merger of $\lesssim 1$ Gyr. The inferred stellar mass of the merged galaxy from the observed H2-to-H i ($\sim 5.6$) and dust-to-gas ratios (GDR$\sim 28$) is $>10^{11.5}$ M⊙, while it is unlikely that NGC 1316 is now experiencing a major-major considering the ordered distribution of its stellar component. To explain the observed H2-to-H i ratio and GDR, additional processes should be taken into account such as an effective dust formation and conversion from atomic to molecular gas during the interaction (section 4.3).

• •

The nuclear jet bends at the NW Shell and the SE Blob suggesting an interaction between the jet and ISM. Our data support the scenario that the nuclear jet is bent due to the buoyant effect and/or ram pressure from the background ISM (section 5).

{ack}

We thank the anonymous referee for his/her comments which improve our paper. We are grateful to Dr. Nao Suzuki and Dr. Saul Perlmutter for kindly providing us with Keck/LRIS data of NGC 1316. KMM acknowledges Dr. Asao Habe and Dr. Hiroshi Nagai for meaningful discussions on the evolution of gas in the elliptical galaxies and the interpretation of the observed kinematics of molecular gas in NGC 1316, respectively. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2017.1.00129.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. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. This project has received funding from the European Research Council under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 679627; project name Fornax). Parts of this research were conducted with the support of Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. T.M. was supported by JSPS KAKENHI Grant Number JP 16H02158. D.E. was supported by JSPS KAKENHI Grant Number JP 17K14254. F.E. was supported by JSPS KAKENHI Grant Number JP 17K14259.