The MASSIVE survey - XI. What drives the molecular gas properties of early-type galaxies

TL;DR

This study investigates the molecular gas properties of massive early-type galaxies, revealing that gas presence is independent of many galaxy parameters and is mainly supplied by gas-rich mergers, with some influence from environment and galaxy dynamics.

Contribution

First comprehensive analysis of molecular gas content in a representative sample of massive early-type galaxies, linking gas presence to galaxy dynamics and environment.

Findings

Detection rate of molecular gas is about 22-25%.

Slow rotators can host molecular gas, but less frequently than fast rotators.

Galaxies in dense environments have lower molecular gas fractions.

Abstract

In this paper we study the molecular gas content of a representative sample of 67 of the most massive early-type galaxies in the local universe, drawn uniformly from the MASSIVE survey. We present new IRAM-30m telescope observations of 30 of these galaxies, allowing us to probe the molecular gas content of the entire sample to a fixed molecular-to-stellar mass fraction of 0.1%. The total detection rate in this representative sample is 25$^{+5.9}_{-4.4}$%, and by combining the MASSIVE and ATLAS$^{\rm 3D}$ molecular gas surveys we find a joint detection rate of 22.4$^{+2.4}_{-2.1}$%. This detection rate seems to be independent of galaxy mass, size, position on the fundamental plane, and local environment. We show here for the first time that true slow rotators can host molecular gas reservoirs, but the rate at which they do so is significantly lower than for fast-rotators. Objects with a higher velocity dispersion at fixed mass (a higher kinematic bulge fraction) are less likely to have detectable molecular gas, and where gas does exist, have lower molecular gas fractions. In addition, satellite galaxies in dense environments have $\approx$0.6 dex lower molecular gas-to-stellar mass ratios than isolated objects. In order to interpret these results we created a toy model, which we use to constrain the origin of the gas in these systems. We are able to derive an independent estimate of the gas-rich merger rate in the low-redshift universe. These gas rich mergers appear to dominate the supply of gas to ETGs, but stellar mass loss, hot halo cooling and transformation of spiral galaxies also play a secondary role.

Figures

Tables

Table 1. Table 1 : Properties of the MASSIVE representative sample objects.

Galaxy	Distance	$M_{K_S}$	$L_{K_S}$	${\sigma }_e$	Re	log10(M${}_{H_2}{}^{})$	Ref
	(Mpc)	(mag)	(L⊙)	(km/s)	(kpc)	(M⊙)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
IC0310	77.5	-25.35	11.45	205.	5.7	= 9.02	1
NGC0057	76.3	-25.75	11.61	251.	10.0	8.60	2
NGC0227	75.9	-25.32	11.44	262.	5.4	8.52	2
NGC0315	70.3	-26.30	11.83	341.	8.5	8.81	6
NGC0383	71.3	-25.81	11.64	257.	7.1	= 9.23	1
NGC0410	71.3	-25.90	11.67	247.	10.9	8.66	2
NGC0467	75.8	-25.40	11.47	247.	9.0	= 8.77	2
NGC0499	69.8	-25.50	11.51	266.	5.3	8.62	2
NGC0507	69.8	-25.93	11.68	257.	13.0	8.68	2
NGC0533	77.9	-26.05	11.73	258.	15.3	8.67	2
NGC0547	74.0	-25.83	11.64	232.	7.1	8.76	2
NGC0665	74.6	-25.51	11.52	164.	4.9	= 9.18	1
NGC0708	69.0	-25.65	11.57	219.	16.5	= 8.83	3
NGC0741	73.9	-26.06	11.74	289.	9.6	8.66	2
NGC0890	55.6	-25.50	11.51	194.	8.2	8.79	2
NGC0910	79.8	-25.33	11.44	219.	9.9	8.12	2
NGC0997	90.4	-25.40	11.47	215.	10.3	= 9.26	1
NGC1060	67.4	-26.00	11.71	271.	12.0	7.89	4
NGC1129	73.9	-26.14	11.77	259.	10.8	8.81	6
NGC1132	97.6	-25.70	11.59	218.	14.6	8.61	1
NGC1167	70.2	-25.64	11.57	172.	10.1	= 8.52	4
NGC1453	56.4	-25.67	11.58	272.	7.9	8.91	2
NGC1497	87.8	-25.31	11.44	190.	7.9	= 9.10	1
NGC1573	65.0	-25.55	11.53	264.	7.9	8.66	2
NGC1600	63.8	-25.99	11.71	293.	12.7	8.41	7
NGC1684	63.5	-25.34	11.45	262.	9.0	= 9.20	1
NGC1700	54.4	-25.60	11.55	223.	6.0	8.31	7
NGC2256	79.4	-25.87	11.66	259.	16.8	8.54	2
NGC2258	59.0	-25.66	11.58	254.	10.1	8.81	1
NGC2274	73.8	-25.69	11.59	259.	10.1	8.58	2
NGC2320	89.4	-25.93	11.68	298.	8.4	= 8.41	5
NGC2418	74.1	-25.42	11.48	247.	7.1	8.68	7
NGC2513	70.8	-25.52	11.52	253.	8.2	8.49	2
NGC2672	61.5	-25.60	11.55	262.	4.3	8.65	2
NGC2693	74.4	-25.76	11.62	296.	5.5	8.40	2
NGC2783	101.4	-25.72	11.60	264.	18.7	8.85	8
NGC2832	105.2	-26.42	11.88	291.	10.8	8.59	2
NGC2892	101.1	-25.70	11.59	234.	11.4	8.84	3
NGC3158	103.4	-26.28	11.82	289.	8.1	8.81	2
NGC3805	99.4	-25.69	11.59	225.	7.9	8.30	9
NGC3816	99.4	-25.40	11.47	207.	8.8	8.67	2
NGC3842	99.4	-25.91	11.68	231.	11.6	8.15	9
NGC3862	99.4	-25.50	11.51	232.	19.2	= 8.49	10
NGC4055	107.2	-25.40	11.47	270.	7.1	= 8.72	6
NGC4073	91.5	-26.33	11.84	292.	10.2	8.90	2
NGC4472	16.7	-25.72	11.60	258.	14.3	7.25	11
NGC4486	16.7	-25.31	11.44	336.	5.7	= 6.70	1
NGC4649	16.5	-25.36	11.46	340.	5.4	7.83	1
NGC4839	102.0	-25.85	11.65	275.	14.4	8.88	2
NGC4874	102.0	-26.18	11.78	258.	15.8	8.86	2
NGC4889	102.0	-26.64	11.97	337.	16.3	8.93	2
NGC4914	74.5	-25.72	11.60	225.	11.3	8.79	2
NGC5208	105.0	-25.61	11.56	235.	9.3	= 9.49	1
NGC5252	103.8	-25.32	11.44	196.	7.9	8.86	1
NGC5322	34.2	-25.51	11.52	239.	3.3	7.76	11
NGC5353	41.1	-25.45	11.49	290.	4.8	$=$ 8.28	12
NGC5490	78.6	-25.57	11.54	282.	7.4	8.73	3
NGC5557	51.0	-25.46	11.50	223.	3.6	7.92	11
NGC6482	61.4	-25.60	11.55	291.	4.9	8.62	1
NGC7052	69.3	-25.67	11.58	266.	9.2	= 9.63	3
NGC7265	82.8	-25.93	11.68	206.	12.7	8.65	2

Table 2. Table 2 : continued

Galaxy	Distance	$M_{K_S}$	$L_{K_S}$	${\sigma }_e$	Re	log10(M${}_{H_2}{}^{})$	Ref
	(Mpc)	(mag)	(L⊙)	(km/s)	(kpc)	(M⊙)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
NGC7274	82.8	-25.39	11.47	244.	9.4	8.55	2
NGC7550	72.7	-25.43	11.48	224.	9.9	= 8.85	2
NGC7556	103.0	-25.83	11.64	243.	13.2	8.56	1
NGC7618	76.3	-25.44	11.49	265.	6.2	8.60	2
NGC7619	54.0	-25.65	11.57	277.	9.0	7.52	1
NGC7626	54.0	-25.65	11.57	250.	7.0	7.90	10

Table 3. Table 3 : Observational parameters and derived molecular gas masses/limits for the sample ETGs

Name	Dbeam1-0	Peak1-0	RMS1-0	${\delta }_{1-0}^{\mathrm{chan}}$	$\int S_{\nu }\delta V$1-0	Peak2-1	RMS2-1	${\delta }_{2-1}^{\mathrm{chan}}$	$\int S_{\nu }\delta V$2-1	$W_{20}$	log10(M${}_{H_2}$)
	(kpc)	(mJy)	(mJy)	(km s-1)	(Jy km s-1)	(mJy)	(mJy)	(km s-1)	(Jy km s-1)	(km s-1)	(${10}^8$ M⊙)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)
NGC0057	8.1	–	7.66	52.	–	–	19.73	52.	–	–	< 3.94
NGC0227	8.1	–	6.50	52.	–	–	14.39	52.	–	–	< 3.31
NGC0410	7.6	–	10.27	52.	–	–	33.61	52.	–	–	< 4.62
NGC0467	8.1	17.94	6.09	78.	8.67 $\pm$ 1.57	30.31	12.29	72.	13.89 $\pm$ 2.95	370.	5.87 $\pm$ 1.06
NGC0499	7.4	–	8.21	52.	–	–	14.18	52.	–	–	< 4.17
NGC0507	7.4	–	11.11	52.	–	–	27.40	52.	–	–	< 4.79
NGC0533	8.3	–	10.94	52.	–	–	25.30	52.	–	–	< 4.71
NGC0547	7.9	–	10.70	52.	–	–	23.02	52.	–	–	< 5.74
NGC0741	7.9	–	9.44	52.	–	–	24.25	52.	–	–	< 4.57
NGC0890	5.9	–	12.67	52.	–	–	33.87	52.	–	–	< 6.12
NGC0910	8.5	–	4.84	52.	–	–	13.62	52.	–	–	< 1.32
NGC1453	6.0	–	14.33	52.	–	–	5.83	50.	–	–	< 8.07
NGC1573	6.9	–	8.14	52.	–	–	17.03	52.	–	–	< 4.58
NGC2256	8.5	–	9.37	52.	–	–	20.35	52.	–	–	< 3.50
NGC2274	6.0	–	6.83	52.	–	–	15.78	52.	–	–	< 3.81
NGC2513	6.9	–	6.41	52.	–	–	16.89	52.	–	–	< 3.09
NGC2672	8.5	–	10.06	52.	–	–	30.66	52.	–	–	< 4.46
NGC2693	7.9	–	7.60	52.	–	–	20.11	52.	–	–	< 2.54
NGC2832	11.2	–	7.86	52.	–	–	19.90	52.	–	–	< 3.85
NGC3158	11.0	–	6.56	52.	–	–	16.90	52.	–	–	< 6.42
NGC3816	10.6	–	4.96	52.	–	–	9.19	52.	–	–	< 4.69
NGC4073	9.7	–	9.06	52.	–	–	14.84	52.	–	–	< 7.92
NGC4839	10.9	–	10.14	52.	–	–	38.80	52.	–	–	< 7.50
NGC4874	10.9	–	7.94	52.	–	–	19.91	52.	–	–	< 7.31
NGC4889	10.9	–	9.31	52.	–	–	15.90	52.	–	–	< 8.56
NGC4914	7.9	–	6.77	52.	–	–	12.77	52.	–	–	< 6.22
NGC7265	8.8	–	9.00	52.	–	–	24.06	52.	–	–	< 4.42
NGC7274	8.8	–	5.90	52.	–	–	7.07	52.	–	–	< 3.57
NGC7550	7.7	64.85	12.89	26.	11.28 $\pm$ 1.20	153.54	27.41	26.	26.58 $\pm$ 2.56	180.	7.03 $\pm$ 0.75
NGC7618	8.1	–	6.62	52.	–	–	16.79	52.	–	–	< 4.02

Equations

$$M_{H_2}=2m_H\frac{{\lambda }^2}{2k_b}{X}_{\mathrm{CO}}{D}_L^2\int S_v\delta V,$$

$$\Delta \Psi =60^{\circ }\frac{|\cos (\Psi )|}{\cos 60^{\circ }}\frac{\Delta T}{t_{\mathrm{relax}}},$$

$$\Delta M_{H_2}={ϵ}_{sf}\frac{{\mathrm{M}}_{{\mathrm{H}}_2}}{t_{\mathrm{dyn}}}\Delta T,$$

$$\Delta M_{H_2}=f_{\mathrm{ML},\mathrm{iso}}\frac{{\mathrm{M}}_{\ast }}{10^{11}}\Delta T,$$

Full text

††pagerange: The MASSIVE survey – XI. What drives the molecular gas properties of early-type galaxies–References††pubyear: 2015

1 Introduction

Low redshift lenticular and elliptical galaxies have been shown in recent years to not be simple ‘red and dead’ objects, as was often assumed in the past. Various studies have shown that these early-type galaxies (ETGs) often contain cold gas (e.g. Wiklind & Rydbeck, 1986; Phillips et al., 1987; Knapp & Rupen, 1996) and associated star formation (Trager et al., 2000; Yi et al., 2005). These new stars typically form dynamically cold (and often chemically distinct; e.g. Kuntschner et al. 2010; McDermid et al. 2015) structures at the hearts of these systems. These objects provide an ideal laboratory in which to study the processes capable of regenerating galaxies.

This molecular material is important because it provides the fuel for both future star formation, and potentially for cold-mode active galactic nucleus (AGN) accretion (e.g. Hardcastle et al., 2007; Ruffa et al., 2019). Molecular gas has been shown to be present in at least 22% of local ETGs (Welch & Sage, 2003; Combes et al., 2007; Sage et al., 2007; Welch et al., 2010; Young et al., 2011). These volume limited surveys have shown little difference in molecular gas properties between lenticular and elliptical galaxies, but that the specific stellar angular momentum (a kinematic morphology indicator) seems to be important. For instance Young et al. (2011) show that the majority of the molecular gas in ETGs is present in fast-rotators. Such studies were, however, not able to probe significant numbers of slow-rotating ETGs.

Exactly where the molecular gas that regenerates these objects comes from is debated. Internally, the stellar mass-loss rates of massive ellipticals can approach several solar masses a year (e.g. Jungwiert et al., 2001), which, if it were able to cool, could easily maintain reservoirs like those observed. However, the detection rate of molecular gas seems to be independent of galaxy properties such as stellar mass, suggesting either that the majority of the gas in these systems has instead been accreted from external sources, or that some mechanism (such as AGN) may suppress the cooling of this material in a manner that leaves no dependence on galaxy mass.

Various observational studies (e.g. using the kinematic misalignment of gas reservoirs; Sarzi et al. 2006; Davis et al. 2011b; Davis & Bureau 2016, motions of large scale atomic gas distributions; Oosterloo et al. 2010; Serra et al. 2012, gas metallicities and dust-to-gas ratios; Smith et al. 2012; Davis et al. 2015 etc.) have suggested that minor mergers, and accretion of gas from close interactions/flybys dominate the supply of gas to fast-rotating ETGs. However, various open questions remain. Some theoretical studies have suggested that cooling from the hot halo could still be important in supplying the gas to ETGs (Lagos et al., 2014; Negri et al., 2014; Lagos et al., 2015). In addition, the lack of a change in detection rates of cold gas in cluster ETGs, where mergers should be rare, has yet to be explained.

The most massive early type galaxies ($M_{*}>10^{11.5}$M*⊙) provide the most stringent tests of the origin of this regenerating gas. The majority of the stars in ETGs are old (Greene et al., 2015), and thus the mass loss is dominated by low mass AGB stars. These most massive ETGs should thus have $>3$ M⊙* a year of material returned to their interstellar medium (ISM). The star-formation rates of such objects are typically \,\stackrel{{\scriptstyle\raisebox{-0.60275pt}{\textstyle<}}}{{\sim}}\,0.1 M*⊙* yr*-1* (e.g. Davis et al., 2016), and thus their gas reservoirs should be building up. These most massive galaxies also host the largest hot halos (e.g. Goulding et al., 2016), and are typically slow-rotators (e.g. Veale et al., 2017b).

Some of the brightest ETGs are the central galaxy in cluster/group environments, and thus have preferential access to cooling gas from their surrounding inter-group/cluster medium. Such objects are routinely found to have large molecular gas reservoirs (e.g. Edge, 2001; Salomé & Combes, 2003; Russell et al., 2016; Vantyghem et al., 2017; O’Sullivan et al., 2018; Tremblay et al., 2018; Russell et al., 2019), likely as a result of cooling from the hot reservoir (e.g. Gaspari et al., 2018; Babyk et al., 2018). However, not all massive ETGs live in such privileged positions within large halos. For instance, in a complete, volume limited sample of M_{*}\raisebox{-2.58334pt}{,\stackrel{{\scriptstyle\raisebox{-0.60275pt}{$\textstyle>$}}}{{\sim}},}10^{11.5}M_{\odot} galaxies nearly half are isolated, or satellites embedded in a larger halo (Veale et al., 2017b). It is thus crucial to study the molecular gas content of representative samples of massive ETGs in a systematic way.

We aim to address this problem in this work, as part of the MASSIVE project. MASSIVE is a survey of the $\approx$100 most massive galaxies within 108 Mpc, using a combination of wide-field and adaptive-optics-assisted integral field spectroscopy (IFS). This volume-limited survey targets the distinct stellar mass range (M_{*}\raisebox{-2.58334pt}{,\stackrel{{\scriptstyle\raisebox{-0.60275pt}{$\textstyle>$}}}{{\sim}},}10^{11.5}M_{\odot}; $M_{K}{\,<\,}{-}25.3$ mag) that has not been systematically studied to date. Full details of this survey are presented in Ma et al. (2014).

In Davis et al. (2016), paper III of the MASSIVE series, we presented a pilot study attempting to detect molecular gas in these high mass galaxies. We observed 15 MASSIVE galaxies with the IRAM-30m telescope and studied their gas masses, star-formation efficiencies and circular velocities. Within this sample, it seemed that the very most massive objects may have a lower probability of hosting cold gas. However, the sample of objects presented in Paper III were selected via various proxies of the presence of gas, and thus are not representative of the MASSIVE sample as a whole. In this paper we aim to address this, extending our survey to a representative subset of the MASSIVE galaxies. In this way we are able to draw more robust conclusions about the sources of gas in these objects, and how these may change with galaxy properties.

In Section 2 of this paper we present details of our sample selection and the properties of the target objects. Section 3 details the observation parameters and reduction. In Section 4 we present our observational results, and present a toy-model to help explain them in Section 5. Finally we discuss and conclude in Section 6.

2 Sample selection

MASSIVE is a volume-limited, multi-wavelength, spectroscopic and photometric survey of the most massive galaxies in the local universe. The full sample includes 116 galaxies in the northern sky with distance $D<108$ Mpc and absolute $K$-band magnitude $M_{K}{\,<\,}{-}25.3$ mag, corresponding to stellar masses M{}_{*}\raisebox{-2.58334pt}{,\stackrel{{\scriptstyle\raisebox{-0.60275pt}{$\textstyle>$}}}{{\sim}},}10^{11.5}M_{\odot}. The MASSIVE survey volume is more than an order of magnitude larger than that probed by ATLAS3D and there are only 6 overlapping galaxies in the two surveys. Wide-field IFS data, taken using the Mitchell Spectrograph (formerly called VIRUS-P; Hill et al. 2008) at McDonald Observatory, are being obtained for all objects in the MASSIVE survey, enabling us to derive stellar and gas parameters (and study the galaxies stellar populations; Greene et al. 2015) out to beyond $\sim$2 effective radii.

The objects observed in this paper were selected to form a representative sample, which matches the full distribution of stellar masses and velocity dispersions found in the MASSIVE parent catalog. In order to select an unbiased sample we divided our parent sample (which is itself complete, mass-selected and volume-limited) into bins of $K_{\rm s}$-band magnitude and velocity dispersion (0.25 magnitudes wide, 40 km/s high). We rank the objects in each cell by their distance. We then choose the closest 50% of objects in each cell, creating a fully representative sample of 57 objects, which match the full distribution of stellar masses and velocity dispersions found in the parent catalog by construction (see Figure 1). Twenty seven objects from this unbiased sample already have suitable observations (15 of these presented in Davis et al. 2016), and in this paper we present new observations of 30 objects.

In what follows we also include data on an additional 10 objects, such that we utilise all data available from the literature for the MASSIVE sample. Including this data slightly over-represents high mass/velocity dispersion systems (see Figure 1), but we found that their inclusion did not change any of our results. We thus here study a total sample of 67 objects, with observations that allow us to detect gas (or set an upper limit to the gas content) down to a fixed molecular-to-stellar mass fraction of 0.1%.

2.1 Lower mass ETGs

In order to study how the molecular gas properties of ETGs change as a function of mass we also include information from Young et al. (2011), the molecular gas survey of the ATLAS3D sample (Cappellari et al., 2011) of ETGs. ATLAS3D is a volume limited survey, covering all the morphologically classified ETGs with distances $<42$ Mpc, $|\delta-29^{\circ}|<35^{\circ}$ and $|b|>15^{\circ}$, that are brighter than a $K_{s}$-band magnitude of $-21.5$ mag (stellar mass M_{*}\raisebox{-2.58334pt}{,\stackrel{{\scriptstyle\raisebox{-0.60275pt}{$\textstyle>$}}}{{\sim}},}6\times 10^{9} M*⊙*).

Young et al. (2011) observed all of these galaxies with the IRAM-30m telescope, to a fixed sensitivity of 3 mK Ta∗, and thus provides us with a complete, volume limited sample (which is thus dominated by lower mass, fast-rotating systems) without incompleteness biases. See Ma et al. (2014) for a more in depth comparison of these two samples.

3 Data and Derived Quantities

The IRAM 30-m telescope at Pico Veleta, Spain, was used between both the 28th July – 9th September 2015 and the 1st June – 16th July 2018 to observe CO emission in our sample galaxies (proposals 067-15 and 064-18, PI Davis). We aimed to simultaneously detect CO(1-0) and CO(2-1), at rest-frequencies of 115.27 and 230.54 GHz respectively, in the 3 mm and 1 mm atmospheric windows. The beam full width at half-maximum (FWHM) of the IRAM-30m at the frequency of these lines is 21$.\!\!^{\prime\prime}$3 and 10$.\!\!^{\prime\prime}$7, corresponding to physical scales of between 5.7 and 10.8 kpc at the frequency of CO(1-0), and between 2.9 and 5.4 kpc at the frequency of CO(2-1), given the varying distance to these sources.

The Eight MIxer Receiver (EMIR) was used for observations in the wobbler switching mode, with reference positions offset by $\pm$100*′′* in azimuth. The Fourier Transform Spectrograph (FTS) back-end gave an effective total bandwidth of $\approx$4 GHz per window, and a raw spectral resolution of 200 kHz ($\approx$0.6 km s*-1* at 3 mm, $\approx$0.3 km s*-1*  at 1 mm). The Wideband Line Multiple Autocorrelator (WILMA) back-end was used simultaneously with the FTS.

The system temperatures ranged between 154 and 493 K at 3 mm and between 200 and 2160 K at 1 mm. The time on source ranged from 24 to 120 min, being weather-dependent, and was interactively adjusted by the observers to reach our required sensitivity (where we could detect objects with a molecular-to-stellar mass fraction of 0.1%).

The individual $\approx$6 minute scans were inspected, and the baseline removed, using a zeroth-, first- or second-order polynomial, depending on the scan. Scans with poor quality baselines or other problems were discarded. The good scans were averaged together, weighted by the inverse square of the system temperature. We consider emission lines where the integrated intensity has greater than a 3$\sigma$ significance (including the baseline uncertainty; Crocker et al. 2012) to be detected.

We convert the spectra from the observed antenna temperature (T${}_{a}^{*}$) into units of Janskys, utilising the point-source conversion as tabulated on the IRAM website111http://www.iram.es/IRAMES/mainWiki/EmirforAstronomers (5.9 Jy/K at 3mm, 8.0 Jy/K at 1mm).

3.1 Line detections

We detected line emission in two out of the 30 galaxies observed in this work (see Figure 2). This low detection rate reflects the fact that the selection criteria used in the previous biased survey did a good job in selecting the objects which were most likely to be detected. Table 3 lists the properties of the observed spectra, and the detected lines. The CO(1-0) detection in NGC0467 is somewhat marginal ($\approx$4.7$\sigma$), while NGC7550 is well detected ($\approx$10.4$\sigma$).

The ratio of the CO(1–0) and CO(2–1) lines we observe is affected by both the gas excitation temperature and the spatial distribution of the gas. If the observed CO emission were to fill the telescope beam at both frequencies (and the CO is not sub-thermally excited) we would expect a line ratio of one (measured in main beam temperature units). However, if the CO emission is compact compared to the beam then the measured intensity in the CO(2–1) line should be larger by up to a factor of 4 (as the beam at such frequencies covers a 4 times smaller area). In our two detected sources we find a line ratio (in beam temperature units) of 1.2 for NGC0467, and 1.75 in NGC7550. Thus, all else being equal, we expect these gas reservoirs to be somewhat extended, but not to completely fill the IRAM-30m beam. This is consistent with the molecular gas being co-spatial with the dust-discs present in these two sources, which are seen in absorption in optical imaging (see Fig 2).

3.1.1 H2 masses

We estimate H2 gas masses and limits from our CO observations in the standard manner, using the following equation

[TABLE]

where $m_{H}$ is the mass of a hydrogen atom, $\lambda$ is the wavelength, $k_{b}$ is Boltzmann’s constant, $D_{L}$ is the luminosity distance, $\int{S_{v}\delta V}$ is the integrated CO flux density and XCO is your CO-to-H2 conversion factor of choice in units of K km s*-1*.

As in Davis et al. (2016) we here use a Galactic XCO factor of 3$\times$1020 cm*-2* (K km s*-1*)-1 (Dickman et al., 1986). As these ETGs are massive, and mass correlates positively with metallicity, such a value seems reasonable. It is possible, however, that the gas in these galaxies has been accreted from a low-metallicity source. In such a case we would be underestimating the total gas mass.

In order to set upper-limits on the H2 mass of our undetected sources we use the 3$\sigma$ RMS on the CO(1-0) spectrum with 52 km s*-1* channels, assuming a total velocity width of 250 km s*-1*, to allow direct comparison to Young et al. (2011). The detected MASSIVE objects all have larger velocity widths than this, and if we assumed a velocity width of 500 km s*-1* and 100 km s*-1* channels our upper-limits would increase by a factor of $\sqrt{2}$. The gas masses and limits derived are presented in Table 3.

4 Results

Combining this work, and the objects from our pilot sample (and the literature) we now have access to the gas properties of a representative sample, drawn uniformly from the MASSIVE parent volume. Sixty seven of the MASSIVE ETGs have been observed in molecular gas to date (from this work, Davis et al. 2016 and other literature sources; see Table 3). Crucially these objects were not selected on any proxy for the presence of cold gas, and we use this fact in what follows to allow us to probe the average gas properties of this galaxy population.

4.1 Detection rate

From the 67 observed MASSIVE galaxies we detect line emission in 17; a detection rate of 25${}^{+5.9}_{-4.4}$%. This rate is very similar to the detection rate of of 22$\pm 3$% (56/259) found by Young et al. (2011) from a survey of all the lower mass ATLAS3D ETGs, and the 26% reported by in the smaller volume limited studies of Knapp & Rupen (1996) and Welch et al. (2010). Combining the MASSIVE and ATLAS3D surveys yields a joint detection rate of 22.4${}^{+2.4}_{-2.1}$% for ETGs of all masses.

Below we combine the MASSIVE and ATLAS3D surveys to explore if (and how) the molecular gas detection rate varies with the physical properties of ETGs. The detection rate is an interesting quantity, as it tells us the fraction of objects with molecular masses greater than $\approx$0.1% of the stellar mass (the detection threshold in this work). We move on in a later section to consider how the molecular masses themselves correlate with these quantities.

In what follows we utilise the $K_{s}$-band luminosity as a measure of the stellar mass of our objects. Dynamical masses are available for the ATLAS3D objects, but not (to date) for the MASSIVE galaxies. None of the results reported in this paper change if dynamical masses are used instead of $K_{s}$-band luminosities for the ATLAS3D objects.

We also note that the molecular gas observations of MASSIVE and ATLAS3D followed a different rationale. In this work, as described above, we observed each ETG such that we were sensitive to a molecular-to-stellar mass fraction of $\approx$0.1% (at 3$\sigma$). In ATLAS3D, however, each object was observed to the same RMS noise level, meaning that objects nearer to us have deeper observations than those further away. Here we combine these two surveys directly, in order to maximise our sample size, but we note that essentially none of the results we present would change if we instead removed the eleven objects from ATLAS3D with a molecular-to-stellar mass fraction of $<$0.1%. In the few cases where results depend on the choice to include these objects we highlight this explicitly in the text.

4.1.1 Fundamental plane

As highlighted in detail in Cappellari et al. (2013a) many of the properties of ETGs vary over the mass-size, and mass-velocity dispersion planes. These surfaces are projections of the ‘Fundamental Plane’ of early-type galaxies (Djorgovski & Davis, 1987; Dressler et al., 1987). We here consider the detection rate (and mass fraction) of molecular gas as a function of the various fundamental plane parameters.

On the top and right of Figure 3 are histograms, whose bins denote the fraction of ETGs detected by either survey within that bin when the distribution is projected down onto this axis. The error bars shown on each bin are the binomial 1$\sigma$ confidence limits. The black dashed line represents the total detection rate of 22%. This allows us to show that the detection rate is the same for galaxies in all parts of the parameter space.

In the top panel of Figure 3 we show the $K_{s}$-band luminosity of the combined MASSIVE and ATLAS3D CO survey galaxies, plotted against a combination of the effective radius (the radius containing half of the stellar light) and velocity dispersion (measured within that effective radius) which together form an edge-on view of the fundamental plane of ETGs. The observational quantities are all reproduced from

Cappellari et al. (2011, 2013b); Ma et al. (2014); Veale et al. (2017a), and the constants a, b and c (which are used to produce the dashed fit line) are those which produce an edge on view of the plane (Cappellari et al., 2013b). The approximate dividing line between the MASSIVE and ATLAS3D surveys is shown on the figure. Objects which were not detected in molecular gas by either survey are shown as open circles, while detections are shown with red filled symbols, whose size scales with the total H2 mass (between $10^{7}$ and $10^{9.5}$ M*⊙*).

We note that the luminosities and radii of massive ETGs are hard to measure in shallow infrared imaging, and may be somewhat underestimated in the 2MASS survey data we use here (e.g. Lauer et al., 2007). This issue will be addressed in detail in Quenneville et al., in prep, who will present details of a dedicated deep infrared imaging campaign. For consistency with ATLAS3D we use the 2MASS photometry here, but note that none of our conclusions would change if we instead used these new data.

The top panel of Figure 3 clearly shows that molecular gas is detected in ETGs in all parts of the fundamental plane. The detection rate of objects appears to be constant both as a function of galaxy mass and fundamental plane parameters, apart from at very low values of the fundamental plane axis. This signature is likely an artefact, as discussed in Section 4.1.3 below. No molecular gas is detected in the very most massive objects, but the low number statistics mean that we cannot reject the null hypothesis of no correlation between the detection rate and stellar mass. The amount of gas present (as indicated by the size of the symbols) also appears to not correlate with either axis (see Section 4.2).

4.1.2 Mass – Size plane

In the middle panel of Figure 3 we show the mass – size projection of the fundamental plane, plotting the $K_{s}$-band luminosity of the combined MASSIVE and ATLAS3D CO survey galaxies against their effective radius (the radius containing half of the stellar light).

Molecular gas is detected in ETGs in all parts of the mass – size plane. The detection rate of objects appears to be constant both as a function of galaxy mass and size. Once again, the amount of gas present (as indicated by the size of the symbols) also appears to not correlate with either axis (see Section 9).

4.1.3 Mass – Velocity Dispersion plane

In the bottom panel of Figure 3 we show the mass – velocity dispersion projection of the fundamental plane, plotting the $K_{s}$-band luminosity of the combined MASSIVE and ATLAS3D CO survey galaxies against their velocity dispersion measured within one effective radius.

This panel shows clearly that the excess at low fundamental plane parameters is caused by an excess of molecular gas detections at very low velocity dispersions, below $\approx$80 km s*-1*, where the detection rate of molecular gas increases to $\approx$50%. This signature comes entirely from the ATLAS3D survey, and was already reported in Young et al. (2011). It is likely due to misclassification of a small number of low-mass spiral/dwarf galaxies for which no good imaging exists, and thus contaminate the sample.

Despite the constant detection fraction as a function of velocity dispersion shown across most of the combined MASSIVE+ATLAS3D sample, it is clear that the molecular gas detections do not fill the mass – velocity dispersion plane uniformly (with a lack of detections in the top left of the figure). At fixed mass the detection rate of cold molecular gas is dependent on velocity dispersion, at least in the lower mass ATLAS3D sample. The dynamic range covered by MASSIVE is lower, so it is harder to tell if the same is true in the more massive objects. The molecular gas fraction at fixed mass also seems to vary across this plane (as discussed previously by Cappellari et al. 2013a).

In order to show this trend more clearly in Figure 4 we plot the $K_{s}$-band luminosity of the combined sample objects versus ${\eta_{\rm kin}\equiv\frac{\sqrt[3]{M}}{\sigma_{e}}}$, a combination of variables that allows us to rotate the mass – velocity dispersion plane to remove its mass-dependent slope. $\eta_{\rm kin}$ correlates with various measures of galaxy morphology (such as Sersic index, concentration index, and stellar mass surface density), and can be interpreted as a kinematic bulge fraction at fixed stellar mass (e.g Cappellari et al., 2013a). This Figure clearly shows that the gas detection rate varies significantly across this plane. At low $\eta_{\rm kin}$ the detection rate of cold gas is $<$10%, while it increases to $\approx$60% at high $\eta_{\rm kin}$ (low kinematic bulge fraction at fixed mass). The low $\sigma_{e}$ outliers discussed above do not drive this trend, as it persists when all objects with $\sigma_{e}<$80 km s*-1*are removed.

It is less clear if the detection rate in the MASSIVE sample alone correlates with $\eta_{\rm kin}$ in the same way as for the lower mass ATLAS3D sample. A Kolmogorov-Smirnov test provides some evidence that the distribution of $\eta_{\rm kin}$ values for MASSIVE detections differs significantly from that for ATLAS3D detections (at $>5\sigma$). This difference seems driven by the slow-rotators; when these objects are removed the $\eta_{\rm kin}$ distributions become statistically indistinguishable. We will discuss this correlation further in Section 6.

4.1.4 Colour

One would expect to find a correlation between the presence of gas in an ETG and its colour, especially in the Ultraviolet bands, which are sensitive to the presence of young stars. In Figure 5 we plot such a figure, showing the $K_{s}$-band luminosity of the combined MASSIVE and ATLAS3D CO survey galaxies against their Near Ultraviolet (NUV) minus $K_{s}$ colours (from GALEX and 2MASS, respectively). NUV magnitudes for the galaxies used in this work were obtained from the GALEX catalogue server, release GR7. Where multiple observations of the same target exist, we always used the deepest observation. Not all of our objects lie within the GALEX survey area, but all that do are detected in the NUV band. In what follows we have assumed that ignoring those ETGs that fall outside the GALEX survey footprint does not bias our results.

UV emission in ETGs can typically arise from two sources. The first is OB stars indicating ongoing star formation. The second is emission from evolved stars (the UV upturn effect; e.g. Yi et al. 2005). The UV-upturn is, however, important mostly in the FUV, and if even a small amount of young stars are present they will dominate the UV emission. In Figure 5 clearly the tail of low-mass galaxies with blue $NUV$-$K$ colours are much more likely to be molecular gas rich, as already shown in Young et al. (2014). We note that excluding galaxies with molecular-to-stellar mass fraction below 0.1% in ATLAS3D (where our combined survey is not complete) strengthens the colour trend, as these systems are preferentially found on the $NUV$-$Ks$ red sequence. The addition of the MASSIVE sample, however, shows how insensitive even ultraviolet colours are to the presence of gas (and star formation) in the most massive galaxies.

4.1.5 Specific angular momentum

A variety of authors have shown that classifying ETGs by their specific angular momentum allows one to predict their properties (and constrain their evolutionary histories) in a physically motivated way.

Young et al. (2011) were able to study the molecular gas detection rate of the ATLAS3D sample as a function of $\lambda_{Re}$, a simplified and dimensionless version of the luminosity-weighted specific angular momentum. They found that there was no dependence of the detection rate with $\lambda_{Re}$ amongst the fast-rotators ($\lambda_{Re}$/$\sqrt{\epsilon}$>0.31), but that the detection rate was very much lower amongst the slow-rotators. Indeed, none of the CO detected objects classified as slow-rotators were true bulge/dispersion dominated systems. NGC1222, is an ongoing merger (Beck et al. 2007), NGC4550 is an object with two counter-rotating stellar discs (denoted 2CR objects; Crocker et al. 2009; Krajnović et al. 2011), and NGC4476 is a complex case with a large kinematically decoupled core (Young, 2002; Davis et al., 2011b). This lead to speculation that the mechanisms that supply gas to ETGs were suppressed in true slow-rotators. This study was hampered, however, by a lack of true bulge/dispersion dominated systems within the survey volume. With MASSIVE we are in a prime position to investigate this further.

In the top panel of Figure 6 we show the $\lambda_{Re}$ vs $\epsilon$ diagram (Emsellem et al., 2007, 2011) for the MASSIVE galaxies in red, and ATLAS3D in grey. As before, objects which were not detected in molecular gas are shown as open circles, while detections are shown with filled circles, whose size scales with the total H2 mass (between $10^{7}$ and $10^{9.5}$ M*⊙*). $\lambda_{Re}$ and $\epsilon$ are taken from Emsellem et al. (2011) for ATLAS3D and Ene et al. (2018) for MASSIVE. The empirical divide between slow- and fast-rotators from Emsellem et al. (2011) is shown as a dashed line. As in Figure 3 we show histograms at the top and right of the plot, in order to show if the detection rate is the same for galaxies in all parts of the parameter space. For clarity we also present $\lambda_{Re}$/$\sqrt{\epsilon}$ vs $K_{s}$-band luminosity for the sample objects in the bottom panel of Figure 6, which allows a clearer separation of fast- and slow-rotators.

It is clear from Figure 6 that with MASSIVE we significantly increase the number of slow-rotators with limits on their molecular gas content. We detect cold gas both in MASSIVE fast-rotators, and in some of the slow-rotators (NGC0708, NGC1684, NGC3862, NGC7052, NGC7550). All of these systems appear to be true bulge/dispersion dominated systems based on the data available, apart from one which has a kinematic twist (NGC1684; Ene et al., 2018). The detection rate of cold gas is consistent with being constant as a function of both $\lambda_{Re}$ and $\lambda_{Re}$/$\sqrt{\epsilon}$ in the fast-rotating population, but drops to $8.6_{-2.6}^{+3.0}$% in the slow-rotators. If one were to remove the merging/2CR systems from ATLAS3D this drops to $5.2_{-1.7}^{+2.9}$%. Overall it seems that slow-rotators are able to host cold molecular gas reservoirs, but the frequency with which they do so is significantly reduced. We will discuss the mechanisms potentially driving this in Section 6.

Slow rotators dominate at the highest masses, so given their lower detection rate one might expect the molecular gas detection rate to drop in the most massive galaxies. The histograms in the bottom panel of Figure 6 show that there is a drop in the detection rate at the highest masses, that could well be due to the dominance of slow rotators, but this is not significant given our sample size.

4.1.6 Environment

Environment is another potential driver of the molecular gas content in galaxies. As mentioned above, many of the mechanisms thought to be fuelling the cold gas in ETGs should be shut down in cluster environments (at least for satellites, which cannot get gas from cooling flows). In ATLAS3D Young et al. (2011) did not find any dependance of the molecular gas detection rate on environment in the single cluster with their volume. In MASSIVE we probe systems with dark matter halo masses from $\approx 6\times 10^{11}$ to $2\times 10^{15}$ M*⊙* (Veale et al., 2017b), giving us the ability to test this once more. MASSIVE includes the largest ETGs that lie in several different clusters, and a larger number of central brightest group galaxies (BGGs) which may have other mechanisms fuelling their gas reservoirs. Our sample does not, however, include many brightest cluster galaxies, and we refer the reader to other works on these distinct objects (e.g. Edge 2001; Salomé & Combes 2003).

In Figure 7 we show the global overdensity indicator $\delta_{g}$ (the luminosity-weighted galaxy density contrast calculated over a 5.7 Mpc Gaussian kernel), and the local overdensity indicator $\nu_{10}/<\nu_{10}>$ (the normalised luminosity density within a region enclosing the 10 nearest neighbours of each galaxy), both reproduced here from Veale et al. (2017b), plotted against the $K_{s}$-band luminosity of the combined MASSIVE and ATLAS3D CO survey galaxies.

The plots are laid out as above, with round coloured symbols denoting molecular gas detections, and the size of those symbols encoding the molecular gas mass.

In both panels there does seem to be a slight drop in the molecular gas detection rate at intermediate densities, but this is not significant. Overall we thus conclude that there is no strong dependance of the molecular gas detection rate on environment density.

The position a galaxy occupies within its halo can be equally as important as the density of galaxies. For example an object that is the brightest in its vicinity is likely to lie close to the centre of the halo, where gas cooling from the hot gas halo can accumulate. Satellite galaxies, however, do not have this advantage. To test if the molecular gas detection rate depends on the position of an ETG in its halo, we plot in Figure 8 the molecular gas detection rate for objects classified as brightest group galaxies, satellites and isolated by Veale et al. (2017b). Once again we show the global detection rate as a dashed line, in order to show that the detection rate is the same for galaxies in all environments.

As Figure 8 shows overall we do not find any significant trend in molecular gas detection rate with the position of each object within its halo. We do see a marginal enhancement of the detection rate in isolated galaxies, and a corresponding drop in the detection rate of satellite objects, but this is not hugely significant. These differences are detected at less than the $<2\sigma$ level, so we would need a significant enhancement in the sample size to determine if these effects are real.

4.1.7 Ionised gas emission

Ionised gas emission is present in the majority of early-type galaxies (e.g. in $\approx$75% of ATLAS3D objects; Sarzi et al., 2006; Davis et al., 2011b). Within MASSIVE at least 38% of objects show clear [OII] emission with an equivalent width $>$2 Å (Pandya et al., 2017).

Insufficient overlap between the MASSIVE and ATLAS3D samples (in terms of emission line selection/sensitivity) exists to compare the properties of the molecular emission with those of the ionised medium systematically. However, all of the molecular gas detected objects presented here that are also studied in Pandya et al. (2017) show detectable [OII] emission. The inverse is not true, however, as some objects without detectable molecular gas do have detectable ionised gas emission. The same trend was found in the lower mass ATLAS3D sample (Sarzi et al., 2006; Davis et al., 2011b). This is consistent with the expectation that an ionised gas layer would form in the outer regions of molecular clouds (where self-shielding is not yet efficient) and in Hii regions within the clouds, in addition to the presence of a dynamically hotter (Davis et al., 2011b; Levy et al., 2018), more diffuse component which can be maintained by ionisation from the older stellar population (e.g. Sarzi et al., 2010).

4.1.8 X-ray emission

Given the possible importance of the hot halo in both quenching and refuelling ETGs, we searched for a correlation between the properties of X-ray emission (that arises from cooling hot-halos) and the detection rate of cold gas. Suitable observations unfortunately are not available for all sample members. We combined the X-ray observations presented in Sarzi et al. (2013) and Goulding et al. (2016) for the available sample objects. We do not find any correlation between hot halo properties (X-ray luminosity, X-ray temperature, difference between the X-ray temperature and the Virial expectation) and the presence/absence of molecular gas, or the amount of gas present in this limited sub-sample of galaxies (see Figure 16 in Appendix B). This lack of correlation will be discussed further in Section 6.

4.2 Gas mass and gas fraction

In the previous section we attempted to determine the physical parameters of galaxies that might be important in controlling the acquisition of cold gas. Now we move on to understand if any of these processes control the amount of gas present, either in an absolute sense or as a function of the mass of the stars within each system.

4.2.1 Stellar mass

As argued above, most of the internal processes that can return material to the ISM of ETGs are related in some way to the mass of stars present. In addition many other gas supply/depletion processes in galaxy evolution also correlate with the stellar (or halo) mass of the system (e.g. the formation of virial shocks, the masses of supermassive black holes, etc). It is thus crucial to know if the mass of ISM material present in ETGs correlates directly with the stellar mass (and thus gas fractions are fairly constant across the ETG population). Previous studies have found no such dependance of the gas mass on stellar mass (e.g. Young et al., 2011), suggesting the gas is dominantly supplied by external processes. We adopt this as our null hypothesis, and will here determine if there is any significant evidence of variation from this within our new observations.

The combination of two different volume limited samples, which have observations that reach different sensitivity levels complicates our analysis. In the top panel of Figure 9 we show the H2 masses of detected ATLAS3Dand MASSIVE galaxies, plotted against $K_{s}$-band luminosity, while in the bottom panel we show the molecular gas fractions (here defined as $M_{\rm H_{2}}/L_{\rm Ks}$). The two surveys cover disjoint regions of the parameter space, and it seems at first glance as if the more massive objects have much higher H2 masses (and thus gas fractions that do not follow the induced trend with stellar mass). This effect is, however, entirely caused by the different sensitivity limits of the two surveys, as indicated by a dashed lines in the top panel of Figure 9. We note that in each case a small fraction of objects had deeper observations available, which we include for completeness. Removing these galaxies would not change any of the results of this paper.

In order to determine if there is any true deviation away from our null hypothesis we create a Monte-Carlo simulation. This simulation incorporates these detection limits in order to produce a sample of model data points which correspond to our null hypothesis (no dependance of the gas mass on stellar mass). In order to do so we start by drawing galaxy masses randomly from the observed mass distributions of the MASSIVE and ATLAS3D surveys. We then assign a random amount of gas to be present in each object, by drawing from the H2 mass function of Young et al. (2011), that was fitted with a Gaussian at the high mass end, and was assumed constant below the completeness limit down to an H2 mass of $10^{6}$ M*⊙*, as in Davis & Bureau (2016). Removing the low-mass cutoff does not affect our results. In the subsample representing ATLAS3D we assign each object a random distance between 10 and 42 Mpc, and consider the object as detected if the CO flux from that gas reservoir would be detectable at that distance. We note that the true distribution of distances is not uniform (e.g. $\approx$25% of ATLAS3D galaxies are in the Virgo cluster), but this is a second order effect. For the subset of galaxies representing MASSIVE we impose a detection limit of 0.1% of the stellar mass (with a small scatter of 0.15 dex to match the observations, some of which are slightly deeper due to weather variations, etc).

In this way we are able to generate a model of the distribution of galaxies within our parameter space (shown with grey shading representing galaxy density on Figure 9), which we can compare to the observed distribution to determine if we can reject the null hypothesis of no correlation between the variables. Both Kolmogorov-Smirnov and Mann–Whitney U tests are unable to reject the possibility that the observed and model distributions are drawn from the same parent distribution, with no evidence for variation greater than $\approx 1\sigma$. This suggests, once again, that the amount of H2 present in ETGs is not correlated with the stellar mass of the host galaxy, and thus that stellar mass loss is not the dominant source of this material.

Given the lack of correlation between gas mass and stellar mass, it is clear that we would expect to find a weak correlation between the gas-to-stellar mass fraction and stellar mass, as seen in the bottom panel of Figure 9. This correlation arises simply because the typical amount of gas in each object stays the same, but the mass of the galaxy changes (see e.g. Young et al., 2011; Agius et al., 2015). We repeated the test described above for the gas-to-stellar mass fractions, and again show the results as grey contours. Again Kolmogorov-Smirnov and Mann–Whitney U tests are unable to reject the possibility that the observed and model distributions are drawn from the same parent distribution, suggesting no additional correlation between these variables is detectable, beyond that induced by the correlated axes.

4.2.2 Specific angular momentum

We showed in Section 4.1.5 above that the specific angular momentum of galaxies correlates with the detection rate of cold gas. Slow rotating ETGs seem to be detected far less often than their fast-rotating counterparts. These detected slow rotating systems have mean gas masses of (6.5$\pm 1.6)\times 10^{8}$ M*⊙, while detected fast-rotating galaxies have an average molecular gas mass of (4.5$\pm 1.2)\times 10^{8}$ M⊙*. Slow rotators are systematically more massive than fast rotators, so despite this higher average gas content, the median gas fraction of detected slow rotators is actually lower than fast rotators (0.23$\pm$0.08% vs 0.3$\pm$0.1%). However, given the significant scatter seen, our low number statistics, and the different selection functions of our two surveys these differences are not significant at more than $\approx 1\sigma$. We thus conclude that there is no strong evidence of the specific angular momentum itself correlating with the amount of cold gas present in ETGs.

4.2.3 Environment

As discussed above, environment is another potential driver of the molecular gas content of galaxies. In Figure 10 we show the global overdensity indicator $1+\delta_{g}$ (from Veale et al. 2017b) for our objects plotted against their H2 masses (top panel), and H2 mass fractions (bottom panel). The MASSIVE detections are shown as red points, ATLAS3D detections as blue symbols, with the non-detections shown as upper-limits with an open grey triangle.

Both of these panels appear to show strong trends, with the most isolated objects (with low values of $1+\delta_{g}$) having larger H2 masses, and thus higher molecular-to-stellar mass fractions. We stress, however, that one does need to exercise caution in drawing strong conclusions from this. The galaxies in dense regions of the ATLAS3D volume are typically closer than those in the void regions, due to the presence of the Virgo Cluster. This means that the ATLAS3D observations were more sensitive in this region. If one were to remove all ATLAS3D detections with H2 masses less than $10^{8}$ M*⊙*(indicated with blue crosses in Figure 10), to ensure that the sample is complete, the trend in the top panel appears less significant. We tested this using a generalised Kendall’s $\tau$ correlation test (implemented using the bhkmethod task in iraf), which properly accounts for the significant number of upper-limits in our study. This test was unable to reject the null-hypothesis, and thus suggests that the correlation between environment and the $H_{2}$ mass is not statistically significant ($\tau$=-0.062, probability=0.16).

In the bottom plot, showing molecular-to-stellar mass fractions, however, the trend that exists is more robust. A generalised Kendall’s $\tau$ correlation test that properly takes into account upper-limits (as discussed above) finds a correlation strength of $\tau$=-0.11, with a probability of 0.01, suggesting the result is significant at $\approx$3$\sigma$. This still does not mean that there is a true correlation between environment and molecular-to-stellar mass fraction, however, as the most massive ETGs ($>10^{11}$ M*⊙) are preferentially found in dense environments (see e.g. Fig 7). We are able to test if this affects our conclusions by only including ETGs with $<10^{11}$ M⊙*, as in our sample there is no strong correlation between mass and environment at these lower masses. When restricting ourselves to this lower mass range we find no sign of a weakening of the correlation (in fact it becomes stronger; $\tau$=-0.1422, P=0.008), suggesting there is some true correlation present between the molecular-to-stellar mass fractions and environment. This correlation is not as apparent when comparing with the local over-density indicator $\nu_{10}/<\nu_{10}>$, suggesting whatever is driving this correlation acts over large scales.

A clue to the likely processes driving this correlation can be found by splitting our galaxies by their positions within their dark matter halos. Figure 11 shows the Kaplan–Meier estimators for the gas properties of the combined sample galaxies split by environment, with isolated galaxies shown in green, satellite galaxies in orange and brightest group galaxies in blue (as classified by Veale et al. 2017b). The estimated distribution functions ($S_{K}$) of the cumulative molecular gas masses are shown in the top panel, and for the molecular gas mass to stellar luminosity ratio in the bottom panel. These estimators include the effects of upper-limits (censored data).

In the top panel BGG and isolated systems are found to have similar molecular gas distribution functions, and within the errors cannot be shown to be significantly different. Satellite galaxies, however, have lower average molecular gas masses than BGG/isolated galaxies. The bottom panel of Figure 11 shows that isolated galaxies have the highest molecular gas fractions on average (as they are typically lower stellar mass systems), while BGG and satellites have fractionally less gas. The molecular gas to stellar mass fractions of isolated objects are $\approx$0.6 dex higher on average than satellite and BGG’s. This suggests that group/cluster environments do affect the amount of gas present in ETGs, if not the total detection rates.

4.2.4 Kinematic bulge fraction

Given the strong dependance of the gas detection rate on the kinematic bulge fraction at fixed stellar mass (which we parameterised with $\eta_{\rm kin}\equiv\frac{\sqrt[3]{M}}{\sigma_{e}}$), one might expect to find a correlation between the molecular gas mass (or molecular gas-to-stellar mass fraction) and this parameter. Indeed Cappellari et al. (2013a) report a correlation between molecular-to-stellar gas mass fraction and the position of ETGs on the mass-velocity dispersion plane, after adaptive smoothing of the gas mass fractions of the lower mass ATLAS3D ETGs. Using our combined sample a generalised Kendall’s $\tau$ correlation test (which includes upper-limits) confirms that both molecular gas mass and molecular-to-stellar gas mass fraction correlate with $\eta_{\rm kin}$ at high significance (correlation strength $\tau$=0.21, P$<0.001$).

It seems that the ability to host a cold gas reservoir correlates with $\eta_{\rm kin}$, and as this kinematic bulge fraction changes so does the mean amount of gas present. Despite this, the scatter in H2 masses (and mass fractions) at fixed $\eta_{\rm kin}$ is very large, so this parameter does not give us the ability to predict with confidence the amount of gas present in any individual object. This is discussed further in Section 6.

5 Analytic toy model

In this section we attempt to synthesise the knowledge gained from the correlations discussed above via a simple analytic toy model. This model is based on that presented in Davis & Bureau (2016), but extended to include the effect of environment, transformation of spiral galaxies, and stellar mass loss. The aim of this exercise is to determine in a simplistic way the relative importance of the various processes we suspect of being involved in supplying molecular gas to ETGs, and is not meant to reproduce every property of these systems. For instance, we completely ignore atomic gas in this analysis, where clearly it may be important, especially in field environments (e.g. Serra et al. 2012). Placing such models within the correct cosmological context is better done within semi-analytic or full-hydrodynamic simulations, and is beyond the scope of this work.

5.1 Model details

Our aim is to create a model that can match the overall molecular gas detection rate in ETGs, the molecular gas detection rates in different environments, as well as the mass and mass-fraction functions taken from this work. In addition, we aim to reproduce the gas kinematic misalignment statistics as a function of environment from Davis et al. (2011b). Our model includes a range of physical processes, and each is parameterised using simple analytic expressions. We summarise below the salient details of the model, and highlight how our implementation differs from Davis & Bureau (2016).

The model begins at a lookback time of 7.7 Gyr ($z$=1), and follows the evolution of an arbitrary number of galaxies. These systems have masses drawn uniformly from the observed mass distribution of the combined ATLAS3D and MASSIVE sample. From the assigned stellar masses we calculate the maximum circular velocity at which the gas can rotate by assuming the galaxies follow the broken $K$-band CO Tully-Fisher relation of Davis et al. (2016), and determine the radial extent of the gas, and hence the dynamical time at its outer edge, by drawing uniformly from a log-normal fit to the observed distribution of molecular disc radii in Davis et al. (2013). Molecular gas reaches beyond the turnover of the rotation curve in $\approx 70\%$ of ATLAS3D ETGs (Davis et al., 2011a, 2013), and all of the interferometrically mapped MASSIVE ETGs (Davis et al., in prep) so this should provide a good estimate of the rotation velocity of the gas at the edges of the gas discs for the bulk of the population. We note that using the un-broken Tully-Fisher relation from Davis et al. (2011a) would not change our conclusions.

Initially all of these systems are classed as isolated, and 27% of them are designated as ETGs (as found for ETGs of all masses at $z$=1 in Tamburri et al. 2014). Although this is somewhat unrealistic, as the red-sequence in clusters is already in place by these redshifts (e.g. Faber et al., 2007) the majority of the satellite galaxies in this study reside in groups, which are still assembling. Our results do not strongly depend on the redshift range we choose to probe, and extending the model to start at $z=2$–3 (to better match the mean light weighted ages of our MASSIVE ETGs; Greene et al. 2015) would not change our results. Initially all the ETGs are gas free.

The evolution of these systems is governed by a linked set of equations (described in more detail below), which allow us to roughly approximate the action of gas rich mergers, star formation, stellar mass loss, environmental quenching and strangulation. These allow us to track the gas content, misalignment of the gas reservoir from the stellar kinematic axis and morphology of each system as a function of time.

In each time-step we randomly select galaxies that will accrete gas, such that the average number of gas-rich mergers per galaxy per unit time is equal to $R_{\rm acc}$ (which we here treat as a free parameter). In this simple model, $R_{\rm acc}$ is independent of redshift, galaxy mass and morphology, which may not be true in reality.

Each isolated galaxy has a fixed probability of entering a dense environment during any given time-step, which is constrained such that the correct fraction of satellite and BGGs (based on Veale et al. 2017b) is reached at $z$ $=$ [math]. Once an object enters a dense environment mergers are prevented entirely.

For those systems that do experience mergers, the amount of gas accreted is chosen randomly from the H2 mass function presented in Davis & Bureau (2016). As discussed in that work we assume that after coalescence each one of our mergers results in the creation of a disc of H2 in the centre of the remnant, with a mass that is drawn from the observed distribution, and that this final state is obtained instantaneously. The angle at which gas is accreted onto galaxies is also chosen uniformly, sampling the full range [math]-$180^{\circ}$. If gas already exists in a galaxy, the interaction between the two reservoirs is roughly included by mass-weighting the resulting position angle of the gas (thus lying between the newly accreted and old material).

The relaxation of the gas is tracked at each time-step. The gas relaxes at a rate that has a cosine dependance on the current misalignment angle (i.e. the rate increases as the disk approaches the galaxy plane, where the torque is higher; e.g. Tohline et al. 1982).

[TABLE]

where $\Delta\Psi$ is the change in angle in any given time-step, $\Psi$ is the current misalignment angle, $\Delta T$ is the duration of each time-step, and $t_{\rm relax}$ is a free parameter. The total time taken to relax is calculated such that the gas disc relaxes into the plane from a misalignment of $60\degr$ in a time $t_{\rm relax}$, which we here fix as 25$t_{\rm dyn}$, following the results of van de Voort et al. (2015) and Davis & Bureau (2016).

Star formation in our model is parameterised as follows:

[TABLE]

where M${}_{\rm H_{2}}$ is the current gas mass, $\Delta{M}_{H_{2}}$ the change in the gas mass at any given time-step, $t_{\rm dyn}$ is the dynamical time, and $\Delta T$ is the duration of the time-step. $\epsilon_{sf}$ defines the efficiency of star formation per dynamical time, which we set to 10%, equivalent to the standard star formation law derived from local star-forming galaxies (Kennicutt, 1998). We note that various authors have suggested that ETGs actually have low star formation efficiencies (e.g. Saintonge et al., 2011; Martig et al., 2013; Davis et al., 2014). This was explored in detail in Davis & Bureau (2016), but for simplicity we do not explore this effect here.

We include AGN feedback in our model in a very simplistic way. Each system has a fixed probability of experiencing an AGN outburst, set such that the average galaxy will experience such an episode once every $t_{\rm AGN}$ years. When this happens, we simply remove 99% of the gas in the galaxy instantaneously. This treatment cannot hope to capture the true complexity of AGN feedback processes, but at least roughly may help us estimate the duty cycle of AGN driven quenching events.

In this work we have added two additional pathways by which ETGs can increase their gas reservoir. The first of these is stellar mass loss. As discussed above, the old low-mass AGB stars, which dominate the evolved stellar populations in ETGs, return $\approx 10^{-11}$ M*⊙* yr*-1* (per solar mass of stars in the galaxy) of gas back into the ISM. This material likely shocks and enters the hot halo of the galaxy initially, but being metal loaded and in the dense region of the halo close to disk, it would be expected to cool back onto the galaxy. We here include this in our model by allowing a fraction of this material (fML,iso - which we here treat as a free parameter) to cool back during each time-step.

[TABLE]

where $\mathrm{M}_{*}$ is the stellar mass of the galaxy at that time, and the other parameters are as defined above.

In satellite galaxies, however, this cooling is likely to be suppressed as the hot halo of the system is ram-pressure stripped. Once an object becomes a satellite we reduce the fraction of the mass loss which cools to zero over a timescale ($t_{\rm strip}$), which here we fix to 1 Gyr following the results of Steinhauser et al. (2016). The results we derive are not strongly dependant on the exact value we assume for this parameter. In contrast, for BGGs we allow cooling to continue with a higher efficiency (fML,BGG, which we again leave as a free parameter) due to the higher likelihood of cooling-flows in massive groups/clusters. This material is always assumed to enter the galaxy in a direction which is aligned with the stellar rotation, and if any misaligned gas exists in the system it interacts and torques it towards co-rotation (again via mass-weighting the two resulting position angles).

Due to the shutdown of mergers, and the reduction of cooling mass loss in satellite galaxies due to halo stripping, a second mechanism of producing gas rich satellite ETGs is clearly needed. We know from the morphology-density relation (Dressler, 1980) that environmental effects can morphologically transform and quench galaxies. We include this here by giving each spiral galaxy in our model a chance to transition into a ETG when it is in a cluster environment. Note that we do this in a very simple way, that is agnostic to the physical mechanism that is actually causing the quenching. This probability is set such that we reproduce the relative fraction of spiral and ETGs in dense environments at $z=0$ from Cappellari et al. (2011). The gas fraction of these spirals as they transition is treated as a free parameter, and parameterised as M${}_{\rm H_{2}}$/Mspiral. This gas is assumed always to be co-rotating within its host galaxy.

At the end of the simulation, we are thus left with the molecular gas mass and kinematic position angle of each galaxy. We convolve the kinematic misalignment angles with a Gaussian kernel to match the observational errors, and apply the observational limits of the ATLAS3D and MASSIVE molecular gas surveys to define objects that are detectable.

We can then compare the output to the observational constraints at $z=0$. The constraints we use are the overall molecular gas detection rate, the molecular gas detection rates in different environments, the mass and mass-fraction functions (all taken from this work), and the gas kinematic misalignment statistics as a function of environment from Davis et al. (2011b). We combine these constraints and generate a log-likelihood for each model. We were then able to run a Markov-Chain Monte Carlo (MCMC) fitting process (using the Gibb’s sampler with adaptive stepping from Davis et al. 2017) in order to find the best fit values for the free parameters (the gas rich merger rate, the fraction of stellar mass loss that cools in BGG and isolated galaxies, and the initial gas fraction of morphologically transformed spiral galaxies), and their uncertainties.

5.2 Model results

The results of the modelling process describe above are shown in Figure 12. While such a simple model is never going to be able to capture all the complexity of galaxy evolution, we do obtain a reasonable match to the observational quantities, with the best model having a reduced $\chi^{2}$ of 1.03.

One of the key parameters of our model is the rate of gas rich mergers/accretion events. As Figure 12 shows our constraint on this parameter is degenerate with the fraction of stellar mass loss that cools in isolated galaxies (and somewhat with the typical AGN duty cycle) but only accretion rates in the range $\approx$0.3-1.2 Gyr*-1* are allowed by the model. The 1$\sigma$ estimate of the gas rich merger rate is 0.69${}_{-0.18}^{+0.14}$ Gyr*-1*, which is in reasonable agreement with other studies. For instance, Lotz et al. (2011) presented a range of estimates from different techniques, and estimated a total merger rate of $0.2$ to $0.5$ Gyr*-1* with $75\%$ of these mergers being minor (mass ratios of 1:4 to 1:10). Theoretical studies such as that of Hopkins et al. (2010) find a major+minor merger rate of $\approx 0.2$ Gyr*-1* for objects of $10^{11}$ M*⊙, but predict that this merger rate scales positively with galaxy mass. The more massive galaxies in our sample (with stellar masses $\approx 10^{12}$ M⊙) are predicted to have merger rates of up to $0.9$ Gyr-1*. However, not all of these mergers are likely to be gas rich, especially at low $z$.

As mentioned above, the fraction of expected stellar mass loss material that cools is somewhat degenerate in our model with the merger rate. Our best estimate is 6.2$\pm$2% of the total old-star stellar mass cools to the molecular phase, equivalent to a return rate per solar mass of stars of $\approx(6\pm 2)\times 10^{-13}$ yr*-1*. If gas rich mergers/accretion events are rare (or AGN outbursts happen frequently) then this rate could be as high as $\approx 2\times 10^{-12}$ yr*-1*, which is still an order of magnitude below the known rate of mass loss from AGB stars. This suggests the majority of this material is unable to cool (perhaps due to AGN maintenance mode feedback which is not modelled here; e.g. Bower et al. 2006; Croton et al. 2006), and the hot halos of these systems are building up over time. We note that the material that cools here doesn’t have to be stellar mass loss, as this formalism would apply equally to any cooling from a kinematically aligned hot halo. If the hot halo were kinematically misaligned from the galaxy in some cases, as suggested in Lagos et al. (2015), then this fraction could increase. In brightest cluster/group galaxies, the model suggests a return of $\approx$14$\pm$2% of the total old-star stellar mass material, or $\approx 1.4\times 10^{-12}$ yr*-1*. This suggests even where cooling flows may be present, the majority of the stellar mass loss material must stay in the halo.

In our model the spiral galaxies that transition into ETGs in clusters (due to any quenching mechanism) have a mean gas fraction of M${}_{\rm H_{2}}$/M${}_{\rm spiral}\approx$10% at the moment they enter the ETG population. This is consistent with the gas fraction of typical sprial galaxies (e.g. Young & Scoville 1991), suggesting only a minimal reduction in the gas fraction of these galaxies is required prior to morphological transformation.

Another parameter in the model is the typical duty cycle for galaxy quenching via an AGN. This is fairly poorly constrained by the model itself, and Figure 12 shows that clear covariances exist between this timescale and the other parameters. If AGN driven quenching events are common (occurring $<$2 Gyrs apart on average), then the gas rich merger rate can be higher, a larger fraction of the stellar mass loss of isolated and BGG galaxies can cool, and spiral galaxies can transition with sizeable gas reservoirs, all without violating our observational constraints. If however the timescale is longer (similar to or greater than our best estimate of $\approx 4.7^{+1.1}_{-1.35}$ Gyr) then the constraints on our other parameters become tighter. We note that the timescale we derive is very much longer than the typical AGN duty cycle ($\approx 10^{6}$-$10^{8}$ yr: e.g. Soker et al. 2001;

Mazzotta et al. 2002; Hopkins et al. 2006; Shabala et al. 2008). This suggests that only a small fraction of AGN outbursts go on to significantly affect the gas reservoir of their host galaxies.

6 Discussion & Conclusions

In this paper we studied a representative sample of 67 of the most massive ETGs in the universe, drawn uniformly from the MASSIVE survey. We presented new observations of the molecular gas content of 30 of these galaxies, which combined with measurements from the literature allows us to probe the molecular gas content of the entire sample down to a fixed molecular-to-stellar mass fraction of 0.1%. Two of these objects observed in this work were detected, hosting between 5 and 7$\times$108 M*⊙* of molecular gas.

The total detection rate of cold molecular gas in the MASSIVE representative sample we study here is 25${}^{+5.9}_{-4.4}$%. This rate is very similar to the detection rate of of 22$\pm 3$% found by Young et al. (2011) for lower mass ATLAS3D ETGs. By combining the MASSIVE and ATLAS3D molecular gas surveys we find a joint detection rate of 22.4${}^{+2.4}_{-2.1}$%. This detection rate seems to be independent of mass, radius, position of the galaxy on the fundamental plane, and (both local and global) environment. The detection rate of cold gas does increase in galaxies with blue ultraviolet-to-infrared colours, reflecting the presence of underlying residual star-formation. In the MASSIVE galaxies alone, however, little effect is seen, showing how insensitive even ultraviolet colours are to the presence of gas in the most massive galaxies. As such, red sequence selections should be treated with caution, as they do not guarantee a purely passive galaxy sample.

The position of a galaxy in the mass-velocity plane does seem to affect the molecular gas detection rate. Objects with a higher velocity dispersion at fixed mass (a higher kinematic bulge fraction at fixed stellar mass) are less likely to be detected in molecular gas. We showed that the parameter $\eta_{\rm kin}\equiv\frac{\sqrt[3]{M}}{\sigma_{e}}$ correlates strongly with the molecular gas detection rate, at least in fast rotating ETGs. Fast-rotating objects with high $\eta_{\rm kin}$ have a detection rate of $\approx$60%, while those with low $\eta_{\rm kin}$ are detected $<10$% of the time. Slow rotators do not seem to follow the same trend, suggesting the kinematic bulge fraction is less important in these objects. Assuming it holds, this correlation between $\eta_{\rm kin}$ and molecular gas could act as a predictor for future surveys which hope to find molecular gas in fast-rotating ETGs (e.g. at higher redshift; Spilker et al. 2018).

This strong correlation between kinematic bulge fraction (at fixed mass) and the ability to host cold gas must encode information about the mechanisms (re-)fuelling massive galaxies. For instance, AGN feedback is invoked in galaxy models in order to stop cooling from the halo. Galaxies with higher kinematic bulge fraction at fixed stellar mass also likely have a higher SMBH mass (given the strong black hole mass – velocity dispersion correlation;

e.g. McConnell & Ma 2013; van den Bosch 2016). This past black hole growth will have injected significant energy/heat into the surroundings, which in turn could have acted to suppress halo cooling (see e.g. Martín-Navarro et al., 2016; Terrazas et al., 2017). The growth of a bulge through the action of bars or dry mergers, for instance, may also have helped evacuate or use up any gas remaining in these systems. Future works will be required to separate out, and determine the balance between these varying competing mechanisms.

The other parameter that appears to play a role in the ability of ETGs to host gas is the stellar specific angular momentum. Fast rotating ETGs appear be be detected at approximately the same rate across the $\lambda_{\rm R_{e}}/\sqrt{\epsilon}$ vs L${}_{\rm K_{s}}$ plane. While we show here (for the first time) that true slow rotators can host molecular gas reservoirs, the rate at which they do so is significantly lower than for fast rotators (\raisebox{-2.58334pt}{,\stackrel{{\scriptstyle\raisebox{-0.60275pt}{$\textstyle<$}}}{{\sim}},}8% in SR vs $\approx$22% in FR).

The amount of molecular gas in ETGs does not appear to correlate with galaxy mass, which induces a weak correlation between gas-to-stellar mass ratio and galaxy mass. Environment does, however, appear to make a larger difference, with galaxies in denser environments having $\approx$0.5 dex lower molecular gas-to-stellar mass ratios. Brightest group galaxies appear to be somewhat less affected, likely because they are able to access additional gas from cooling flows. The cause of this lower gas fraction in satellite galaxies, at least in our toy model, is strangulation (Larson et al., 1980), as the gas accretion rate is strongly reduced for satellite galaxies (e.g. van de Voort et al., 2017). Ram-pressure stripping (Gunn & Gott, 1972) may also play a role in the low gas fractions of the observed systems. For instance, it has been shown that the gas reservoirs in ETGs are less extended in dense environments (Davis et al., 2013). Despite this, it is expected that ram-pressure affects the molecular gas in massive galaxies to a lesser degree than strangulation, as the gas is tightly bound at the centre of a deep potential well.

Three possible scenarios are often discussed to explain the presence of molecular gas in ETGs:

The gas could be a remnant of the star forming reservoir of the galaxy before it was morphologically transformed. 2. 2.

The gas could have been regenerated from an internal source, such as stellar mass loss or cooling of the hot halo. 3. 3.

The gas could be acquired externally, from minor/major mergers and cold accretion.

The first of these possibilities is explored somewhat in our toy model, where in cluster environments we include the transformation of spirals into ETGs. In our model around a quarter of the low redshift ETG population in clusters are the remnants of transformed spiral galaxies, which have maintained a detectable gas reservoir since they were morphologically transformed. While these systems are present, they do not dominate the population of gas rich ETGs, and thus other mechanisms also seem to be important.

The constant detection rate of molecular gas as a function of stellar mass we find in this work once again suggests that the stellar mass loss is not a dominant source of ISM in ETGs. As mentioned above the gas return rate should correlate linearly with stellar mass at old ages, leading more massive galaxies to be detected more frequently, which is at odds with our data. In addition, if mass loss were the dominant process then it is much harder to explain the lack of gas in 75% of ETGs without another process (such as AGN feedback) conspiring to keep the detection rate of cold gas constant. Various other studies have come to the same conclusion, but we can now extend this result to even the most massive systems. Indeed, the toy model outlined above suggests $<15\%$ of AGB star mass loss can cool into the ISM of isolated/satellite ETGs, and $<20\%$ in BGGs, in order to not violate the observational constraints. In our models the gas rich merger rate is $\approx$0.6 Gyr*-1*, and the typical accreted gas mass is $\approx$3$\times 10^{7}$ M*⊙* Gyr*-1*, and thus accreted material dominates over in-situ production of gas by a factor of at least 10,000.

The only other internal reservoir of gas that could cool to form a molecular reservoir is the galaxy’s hot halo. This is known to be an important mechanism for the central galaxies in galaxy clusters and groups, and ALMA observations have begun to resolve this process in a variety of systems in recent years (e.g. Russell et al., 2016; Vantyghem et al., 2017; O’Sullivan et al., 2018; Tremblay et al., 2018). However, the importance of halo cooling in isolated systems is less well understood. The correlations we find between the molecular gas detection rate and both kinematic bulge fraction and stellar specific angular momentum raise the possibility that the hot halo is also important in these systems, as cooling rates are predicted to be enhanced in more disky (kinematically colder) galaxies (Negri et al. 2014; Juráňová et al. 2019). If the hot-halo is a dominant source of cold gas then one requires a large number of halos to be kinematically misaligned from their host galaxy in order to match the fraction of kinematically misaligned gas discs in ETGs (e.g. Lagos et al., 2015). A misaligned hot halo could potentially show increased cooling rates due to collisions between stellar mass loss material and the halo. Given the above, one might expect to find a correlation between these parameters and the X-ray emission that arises from cooling hot-halos, but no such correlation was found in this work. This could be due to a real lack of correlation, to the low number of probed systems, or to the difficulty of probing the colder gas $\approx 10^{5}$ K gas in galaxy halos. This issue will be probed further in a future work of this series (Goulding et al. in prep).

In summary, in this work we show that molecular gas reservoirs are common across the ETG population. The kinematic bulge fraction and specific angular momentum of galaxies affect their ability to host such reservoirs, while the amount of gas present is also affected by environment. Mergers appear to dominate the supply of gas, but stellar mass loss, hot halo cooling and transformation of spiral galaxies may also play a role, Further observational and theoretical work will be required to fully disentangle the importance of each of these effects. For instance theoretical works may be needed to predict the observational signatures of hot halo cooling in isolated ETGs. Resolved observations (with e.g. ALMA) will reveal the angular momentum of the molecular gas, and its spatial distribution. These observations can also reveal if this molecular material is important in fuelling AGN, which are very common in massive ETGs, but are often thought to be fuelled by hot-mode accretion. A high-resolution X-ray spectrometer (to replace Hitomi after its unfortunate destruction) would also reveal the rotation of hot halos for the first time, and thus their importance in refuelling galaxies across cosmic time.

Acknowledgments

TAD acknowledges support from a Science and Technology Facilities Council Ernest Rutherford Fellowship, and thanks the referee for helpful comments which improved this paper.

The MASSIVE survey is supported in part by NSF AST-1411945 and AST- 1411642. This paper is based on observations carried out with the IRAM Thirty Meter Telescope, and we thank Claudia Marka and Wonju Kim, the pool supervisors for observational assistance. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). This research also 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.

The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation Grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation.

Appendix B X-rays and Ionised gas