Similar Scaling Relations for the Gas Content of Galaxies across Environments to z ~ 3.5

We study the effects of the local environment on the molecular gas content of a large sample of log($M_{*}$/$M_{\odot}$) $\gtrsim$ 10 star-forming and starburst galaxies with specific star-formation rates (sSFRs) on and above the main-sequence (MS) to $z$ $\sim$ 3.5. ALMA observations of the dust continuum in the COSMOS field are used to estimate molecular gas masses at $z$ $\approx$ 0.5-3.5. We also use a local universe sample from the ALFALFA HI survey after converting it to molecular masses. The molecular mass ($M_{ISM}$) scaling relation shows a dependence on $z$, $M_{*}$, and sSFR relative to the MS, but no dependence on environmental overdensity $\Delta$ ($M_{ISM}$ $\propto$ $\Delta^{0.03}$). Similarly, gas mass fraction (f$_{gas}$) and depletion timescale ($\tau$) show no environmental dependence to $z$ $\sim$ 3.5. At $\langle z\rangle$ $\sim$ 1.8, the average $\langle M_{ISM}\rangle$,$\langle$f$_{gas}\rangle$, and $\langle \tau \rangle$ in densest regions is (1.6$\pm$0.2)$\times$10$^{11}$ $M_{\odot}$, 55$\pm$2%, and 0.8$\pm$0.1 Gyr, respectively, similar to those in the lowest density bin. Independent of the environment, f$_{gas}$ decreases and $\tau$ increases with increasing cosmic time. Cosmic molecular mass density ($\rho$) in the lowest density bins peaks at $z$ $\sim$ 1-2 and this peak happens at $z$ $<$ 1 in densest bins. This differential evolution of $\rho$ across environments is likely due to the growth of the large-scale structure with cosmic time. Our results suggest that the molecular gas content and the subsequent star-formation activity of log($M_{*}$/$M_{\odot}$) $\gtrsim$ 10 star-forming and starburst galaxies is primarily driven by internal processes, not their local environment since $z$ $\sim$ 3.5


INTRODUCTION
It is now well established that the environment of galaxies influences their properties such as morphology, color, star-formation rate (SFR), and likely metallicity, stellar mass, dust, and gas content (e.g. ;Dressler 1980;Kauffmann et al. 2004;Peng et al. 2010;Scoville et al. 2013;Sobral et al. 2016;Darvish et al. 2016).
These studies along with observations of e.g. jelly-fish galaxies in clusters (Owers et al. 2012;Poggianti et al. 2016) suggest that this is likely due to the ram pressure stripping the interstellar medium of galaxies.
The situation is different for the molecular gas as it is less extended and much denser than HI and hence more bound to the host galaxy.
A number of studies support this picture, finding no environmental effects on the molecular gas in galaxies (Stark et al. 1986;Kenney & Young 1989;Lavezzi & Dickey 1998;Koyama et al. 2017).
However, others find a deficit (Fumagalli et al. 2009;Corbelli et al. 2012;Jablonka et al. 2013;Scott et al. 2013; or even an enhancement (Mok et al. 2016) of molecular gas in denser environments than the field.
Thanks to modern radio interferometers, the molecular gas studies based on CO observations are now performed to z ∼ 3-4 (e.g.; Daddi et al. 2010;Saintonge et al. 2011a;Tacconi et al. 2018), suggesting a higher gas mass fraction at higher-z and other potential dependence of the molecular gas on stellar mass (lower gas fraction in more massive systems) and SFR relative to the mainsequence (higher gas fraction for systems above the mainsequence; e.g.; Tacconi et al. 2013;Genzel et al. 2015). However, these studies often target field galaxies and the need for statistically large samples of high-z galaxies in dense environments is essential.
Nonetheless, there are a limited number of molecular gas studies in dense environments at high z, often with limited CO detections (Aravena et al. 2012;Wagg et al. 2012;Rudnick et al. 2017;Hayashi et al. 2017;Lee et al. 2017;Noble et al. 2017;Hayashi et al. 2018) and often with contradictory results likely because of small sample size, different selections, and a small dynamical range in stellar mass (M * ) and the sSFR with respect to the mainsequence (sSFR/sSFR MS ). For example, Noble et al. (2017) find a higher gas fraction and a similar depletion timescale in z ∼ 1.6 clusters than the field, Hayashi et al. (2018) find a higher gas fraction and depletion timescale for a z=1.46 cluster galaxies than the field, and Lee et al. (2017) find a similar gas fraction but likely shorter depletion timescale in denser regions at z ∼ 2.5.
One major problem with molecular mass estimates -Top: Distribution of stellar mass, sSFR relative to the main-sequence, and redshift in four environmental density quantiles of δ 4 (densest), δ 3 , δ 2 , and δ 1 for our high-z sample (z ≈ 0.5-3.5). The parent high-z sample from Scoville et al. (2017) is also shown with blue dashed lines. Note the similarity of distributions in different density quantiles. In all cases, K-S test p-value is > 0.1, ensuring a good control sample selection (Section 4). High-z samples are located on and above the main-sequence and have M * 10 10 M ⊙ . Bottom: Similar to the top panels but for our local universe sample at z ∼ 0.04. Note the similarity of M * , sSFR/sSFR M S , and z distributions in different environmental density bins, showing a good control sample selection (Section 4). Similar to our high-z sample, our local universe galaxies are located on and above the main-sequence with M * 10 10 M ⊙ .
based on CO observations is that they are timeconsuming, often leading to a small sample size. Moreover, converting the excited state CO fluxes into reliable estimates of the gas content in galaxies remains an issue (Carilli & Walter 2013). A much faster (minute-long vs. hour-long) alternative arises from the long-wavelength Rayleigh-Jeans (RJ) dust continuum method (e.g.; Santini et al. 2014;Scoville et al. 2014). This technique has been calibrated over a broad range of CO luminosities (Scoville et al. 2014), successfully leading to ISM mass estimates for 708 galaxies at z ∼ 0.3-4.5 with ALMA (Scoville et al. 2017). This large blind sample increases the chance of finding galaxies that happen to be in dense regions, even though the original selection is not designed to target dense environments a priori. Scoville et al. (2017) studied the molecular gas content of galaxies as a function of z, M * , and sSFR/sSFR MS , providing analytical expressions for gas mass, gas fraction, depletion timescale, etc. The missing parameter is the environment which will be investigated here.
In this paper, by using large samples of galaxies covering a broad range of environments, z, M * , and sSFR/sSFR MS , we address the following question: Given fixed z, M * , and sSFR/sSFR MS , does the local environment affect the molecular gas content of galaxies in the local universe and at higher redshift? For the highz sample, we rely on Scoville et al. (2017) large sample by dividing it into environmental bins. For the local universe, we rely on the ALFALFA HI sample (Haynes et al. 2011) after converting it to molecular gas and placing them into different density bins.
In Section 2 we present the data. The methods used to quantify the local environment of galaxies are developed in Section 3. Control sample selection is given in Section 4. Results are presented in Section 5, discussed in Section 6, and summarized in Section 7. Throughout this work, we assume a flat ΛCDM cosmology with H 0 =70 kms −1 Mpc −1 , Ω m =0.3, and Ω Λ =0.7, a Chabrier initial mass function (IMF; Chabrier 2003), and a CO-to-H 2 conversion factor of α CO =6.5 M ⊙ (K kms −1 pc 2 ) −1 (Scoville et al. 2014). The relevant values from the literature are converted to our adopted IMF and α CO .

High Redshift
The parent sample is from Scoville et al. (2017), based on ALMA observations of the long wavelength dust continuum for a sample of 708 galaxies at z ≈ 0.3-4.5 in the COSMOS field (Scoville et al. 2007). The sample galaxies have stellar masses of log(M * /M ⊙ ) 10 and sSFRs from the main-sequence to ∼ 50 times the MS ( Figure  1). Stellar masses are based on SED template fitting using optical to infrared bands (Laigle et al. 2016). SFRs are estimated using the UV and far-infrared luminosity (Lee et al. 2013(Lee et al. , 2015. For the M * dependence of the   Scoville et al. 2014Scoville et al. , 2017. Note that the parent sample selection is not designed a priori to target galaxies in dense environments. Therefore, we estimate the local environment of the parent sample (Section 3.1) and then place the galaxies into four local density quantiles: δ 1 , δ 2 , δ 3 , and δ 4 (highest density). This results in 102 galaxies in each density bin after controlling for their M * , z, and sSFR/sSFR MS (Section 4). They cover the redshift range z ≈ 0.5-3.5, log(sSFR/sSFR MS ) -0.2, and log(M * /M ⊙ ) 10 ( Figure 1).

Local Universe
The local universe sample is from the ALFALFA HI survey (Haynes et al. 2011) matched with the SDSS DR12 (Alam et al. 2015). Our SDSS galaxies have magnitudes r 17.7, are in the spectroscopic redshift range of 0.02 z 0.12, and they are located in the contiguous northern galactic cap (130.0 RA (deg) 240.0 and 0.0 Dec (deg) 60.0) (Darvish et al. 2018). The mass completeness limit of the SDSS sample is log(M limit * /M ⊙ ) ∼ 10.0 to z=0.12. To make our local universe selection as similar as possible to that of our high-z sample, we further select galaxies that are on and above the MS (log(sSFR/sSFR MS ) > -0.2) and with log(M * /M ⊙ ) 10 ( Figure 1). The local universe MS is based on Salim et al. (2007) (log(SFR)=0.65 log(M * /M ⊙ )−6.33). Stellar masses and SFRs are from the MPA-JHU DR8 catalog (Kauffmann et al. 2003). We convert HI to H 2 mass assuming M H2 /M HI =0.26 (Obreschkow & Rawlings 2009) and further multiply it by 1.36 to account for heavy element contribution 5 . We assign our ALFALFA/SDSS matched galaxies to different local density quantiles (Section 3.2), leading to 117 galaxies in each density quantile after controlling for their M * , z, and sSFR/sSFR MS (Section 4).

High Redshift
We use the slightly modified density field estimation of Darvish et al. (2017) in the COSMOS field. The local environment measurement relies on the adaptive kernel smoothing method (Scoville et al. 2013;Darvish et al. 2015b) estimated over a series of overlapping redshift slices (Darvish et al. 2015b). A mass-complete sample (log(M * /M ⊙ ) 10.4 similar to a volume-limited sample) is used for density estimation. The lower limit of log(M * /M ⊙ ) 10.4 roughly corresponds to the stellar mass completeness limit at the highest redshift of this study (Laigle et al. 2016). This selection is used to minimize the unrealistic underestimation of the density values at higher z due to the Malmquist bias. The density values for our lower mass sample galaxies ( log(M * /M ⊙ )=10-10.4) are estimated by interpolation of the density field to their positions. We define four density bins using the projected surface density quantiles from the low-to highdensity regions: δ 1 , δ 2 , δ 3 , and δ 4 . The dynamical range of environments for our sources is ∼ 2.1 dex, sampling those of the overall galaxy population. The range of den- -Top: Distribution of molecular mass, gas mass fraction, and depletion timescale in different environmental density quantiles for our high-z sample (z ≈ 0.5-3.5), after controlling for M * , sSFR/sSFR M S , and z in different density quantiles ( Figure 1). The parent high-z sample from Scoville et al. (2017) is also shown with blue dashed lines. Note the similarities between different density quantiles. We find no clear environmental dependence of M ISM , fgas, and τ for our high-z sample. In all cases, the K-S test p-value is > 0.1, showing that it is unlikely that the distributions of M ISM , fgas, and τ in different environmental density quantiles are drawn from different parent distributions. Bottom: Similar to the top panels but for our local universe sample at z ∼ 0.04. We do not find a significant environmental dependence of M ISM , fgas, and τ for our local universe sample, in line with our high-z results. sity (Σ) and overdensity (∆) values for the defined density quantiles at different redshifts are given in Table 1. Overdensity (∆) is defined as the density divided by the median density at each redshift (see e.g.; Darvish et al. 2017).
To compare the overdensity values with those of real galaxy groups/clusters, we use a sample of X-ray groups at z < 1.2 in the COSMOS field with measured positions, z, R 200 , and M 200 from Finoguenov et al. (2007). For each reliable group (flag=1), we determine the median overdensity of galaxies (used for density estimation) within R 200 of each group center and in redshift slices around the z of each group. The width of each redshift slice is similar to those used here for density estimation (also see Darvish et al. 2017). For a median M 200 =4.0 × 10 13 M ⊙ and 2.1 × 10 14 M ⊙ group, we obtain a median overdensity of 3.4 and 10.2, respectively. Unfortunately, no such catalog exists at z > 1.2 for comparison with our estimated density values.

Local Universe
We use the density measurements of Darvish et al. (2018) based on the SDSS DR12 with redshift, magnitude, and angular position cuts already defined in Section 2.2. Local density estimations are based on the projected comoving distance to the 10th nearest neighbor to each galaxy, considering only galaxies that are within the recessional velocity range of ∆v=c∆z=±1000 kms −1 to the galaxy of interest, and corrected for incompleteness due to the fiber collision and the Malmquist bias. Similar to the high-z sample, we place galaxies into four density bins. The dynamical range of environments here is ∼ 2.9 dex. Table 1 lists the range of density and overdensity values for the density quantiles.
Using different density estimators at high-and low-z (adaptive kernel smoothing versus 10th nearest neighbor) might lead to a potential bias in presenting the lowand high-z results. However, Darvish et al. (2017) compared the density values with the 10th nearest neighbor and adaptive kernel smoothing methods for a similar sample at high-z and found a good agreement (an offset of < 0.1 dex and a median absolute deviation of < 0.2 dex). Moreover, Darvish et al. (2015b) find an overall good agreement between the estimated density fields using different methods (including the adaptive kernel smoothing and 10th nearest neighbor) over ∼ 2 dex in overdensity values through simulations and also observational data. Hence, the selection of different estimators does not likely have a significant effect in the presented results.

CONTROL SAMPLE SELECTION
Since the estimated M ISM of galaxies likely depends on z, M * , and sSFR/sSFR MS (e.g.; Genzel et al. 2015;Scoville et al. 2017), we control for these parameters when we compare our high-density sample (δ 4 ) with the rest. To do this, we use our high-density sample as the reference 6 . Then, for each galaxy in this sample, we first search for any galaxy in the lower density samples (e.g., δ 3 ) that is within ± 0.2 dex of M * , ± 0.2 dex of sSFR/sSFR MS , and ± 0.2 (±0.005 for the local sample) of redshift of that galaxy. If more than one match is found, we then select the galaxy that has the closest cartesian distance to the galaxy of our interest in the 3D M * -sSFR/sSFR MS -z space. If no match is found, we remove that galaxy from the high-density sample. This results in 102 (117) galaxies in each density bin for our high-z (local universe) sample. Figure 1 shows the normalized distributions of M * , sSFR/sSFR MS , and z for the original high-density sample and all the other lower-density control samples for both the high-z and local universe galaxies. Kolmogorov-Smirnov (K-S) tests are performed to make sure that distributions for different density bins are unlikely to be drawn from different parent distributions in M * , sSFR/sSFR MS , and z. Large K-S p-values (p > 0.1) in all cases guarantee a good control sample selection.
We note that this approach of controlling samples in different environments might lead to some bias and misinterpretation of the results if there is some degree of association between denser environments and stellar mass and/or the sSFR/sSFR MS of star-forming and starburst systems. Therefore, careful analysis is required prior to investigating potential differences between control samples in different environmental density bins. In Appendix A, we discuss some relevant potential biases.

Scaling Relations with Environment
Prior to placing our sample galaxies into different density quantiles, we first fit a power-law dependence of M ISM of all galaxies as a function of z, sSFR/sSFR MS , M * (similar to Scoville et al. 2017), and overdensity (∆). For this analysis, we only rely on our ALMA dustcontinuum high-z sample but similar results are expected for the local universe sample. To find the best fit parameters, we use both a multi-dimensional non-linear Levenberg-Marquardt (LM) least square algorithm and a Monte Carlo Markov chain (MCMC) Bayesian approach (Metropolis-Hastings sampling). The observational uncertainties of M ISM , sSFR/sSFR MS , and M * are taken into account. The results of the LM and MCMC fitting are: Both LM and MCMC methods yield consistent results. As expected, the parametric dependence of M ISM on z, sSFR/sSFR MS , and M * is similar to what was found in Scoville et al. (2017). However, here, we explicitly investigate the potential dependence of M ISM on the environment of galaxies as well (with overdensity of galaxies as a measure of environment). From Equations 1 and 2, we find only a weak dependence of M ISM on the galaxy overdensity, with a power-law of the form M ISM ∝ ∆ 0.03 +0.01 −0.01 . We also investigate the scaling relations for the SFR (similar to Scoville et al. 2017) including the overdensity dependence. As noted by Scoville et al. (2017), we impose a linear dependence of the SFR on M ISM , with M ISM obtained directly from Equation 2 rather than going back to the observed M ISM values. This is essential in order to isolate the star-formation efficiency (SFR/M ISM ) variation with z, sSFR/sSFR MS , M * , and ∆ from the variation of the M ISM with the same parameters. The results of the LM and MCMC fits are: The LM and MCMC results are consistent with no dependence of SFR (and subsequently star-formation efficiency and depletion timescale τ =M ISM /SFR) on overdensity values.
However, we note that the fitting procedure is only applicable to the range of parameters probed here. A galaxy sample covering a much larger range of overdensity values (for example, sampled in dense core of clusters) and other physical parameters are needed to be able to more robustly investigate the dependence of M ISM on galaxy properties including their environment. It is also worth noting that the overdensity might not be fully independent of the other variables such as the stellar mass used in the fitting procedure. This dependence might be particularly important in the extreme regions of environments and high stellar masses. However, the potential relation between environment and stellar mass is still debated (see e.g.; Baldry et al. 2006and Darvish et al. 2016vs. von der Linden et al. 2010).

General Trends with Environment for Control
Samples We also compare galaxy properties in different density quantiles for the controlled samples. Table 2 summarizes the results in this paper. It contains the average and median gas properties of galaxies in different density quantiles and redshifts. After controlling for z, M * , and sSFR/sSFR MS , we compare the total molecu- a δ 1 , δ 2 , δ 3 , and δ 4 refer to density quantiles from the low-to high-density regions, respectively. b Sample size c Average redshift d Average molecular mass, assuming α CO =6.5 M ⊙ (K kms −1 pc 2 ) −1 e Median molecular mass f Average gas mass fraction g Median gas mass fraction h Average depletion timescale i Median depletion timescale lar mass (M ISM ), gas mass fraction (f gas ), and depletion timescale (τ ) in different environmental bins, as shown in Figure 2. We see no clear environmental dependence of M ISM , f gas , and τ in the local universe and out to z ∼ 3.5. We perform K-S test and in all cases, the K-S tow-tailed p-value is > 0.1 (corresponding to a significance of 1.6σ), indicating that it is unlikely that the distributions of M ISM , f gas , and τ in different density bins are drawn from different parent distributions.
These values are reasonably similar to those in the lowest density quantile within the uncertainties ((3.9±0.2)×10 9 M ⊙ , 14.5±0.6%, and 1.3±0.1 Gyr). The uncertainties here are the standard error of the mean. For the high-z sample ( z ∼ 1.8) and in densest regions, we obtain the average values -Redshift evolution of gas mass fraction fgas in different environmental density quantiles (empty circles), along with cluster and field quantities in the literature. Note that the literature studies have their own sample selections that might be different than ours and they are just shown for reference. fgas values in different environmental density bins show that independent of the environment, the gas mass fraction decreases with decreasing redshift and that our data points nearly follow the overall decline in fgas with cosmic time.  Table 2).

Environment and Redshift Dependence of Gas Mass Fraction
We further investigate the environmental dependence of M ISM , f gas , and τ in redshift bins. Note that we control for z, M * , and sSFR/sSFR MS distributions in different density bins for the redshift-binned data as well.
Even after dividing the sample into redshift bins, we still find no significant environmental dependence of M ISM , f gas , and τ within the errors. Figure 3 shows the redshift evolution of f gas in different density bins, along with cluster (Cybulski et al. 2016;Lee et al. 2017;Noble et al. 2017) and field (Tacconi et al. 2010;Saintonge et al. 2011b;Geach et al. 2011;Bauermeister et al. 2013;Kirkpatrick et al. 2014;Cybulski et al. 2016;Schinnerer et al. 2016;Scoville et al. 2017) quantities in the literature. Note that the literature studies have their own sample selections that might be different than ours and they are just shown for reference. Nonetheless, we find that independent of the environment (within the uncertainties), the gas fraction decreases with decreasing redshift. Moreover, our data points nearly follow the overall decline in f gas with cosmic time.  Tables 2 and 3) and  TABLE 3 Fraction of Volume (f V ) Occupied by and the Global Molecular Mass Density (ρ) in Different Environmental Density Quantiles a Fraction of the volume occupied by galaxies located in δ 1 quantile. b Global molecular mass density in δ 1 quantile.   Table 3, where dust-data refers to gas mass estimates based on a metallicity-dependent dust-togas ratio and dust mass estimates using the SFR, and the estimated dust temperature assuming a modified blackbody spectrum and an emissivity index of β=1.5; see section 2.2 of Genzel et al. 2015) are shown. We use the median sSFR/sSFR MS and M * of our sample in all density bins in Scoville et al. (2017) and Genzel et al. (2015) equations. The uncertainties include the model fitting errors and the dispersion in the sSFR/sSFR MS and M * observables (dominated by the latter). Within the uncertainties, our estimated τ values are consistent with these global trends, regardless of their local environment.

Environment and Redshift Dependence of Molecular Mass Density
We also investigate the redshift evolution and environmental dependence of the global molecular mass density (ρ). For this, we rely on the overall molecular mass density at each redshift (see Scoville et al. 2017) and scale it with the fraction of galaxies in each density quantile (25%) and the inverse of the fraction of the volume assigned to each density quantile: where ρ(δ, z) is the molecular mass density for the density quantile of δ and redshift of z, f V (δ, z) is the fraction of the volume assigned to the density quantile δ at z, and Φ(M ISM , z) is the molecular mass function at z. f V (δ, z) can be measured via: In the above equation, V (δ, z) is the volume assigned to the density quantile δ at z and V total (z) is the total volume assigned to all quantiles at that redshift. V (δ, z) is estimated using all the galaxies N (used in estimating the density field) that are located in density quantile δ at z. σ i is the projected surface density of the galaxy i (in Mpc −2 ), and ∆L i is the comoving length (in Mpc) that corresponds to the redshift slice (∆z) within which the projected surface density of the galaxy i is estimated  . We also show the estimation of ρ by directly using the stellar mass functions in densest (HD) and the least dense (LD) regions from the VIPERS (Davidzon et al. 2016) and ORELSE (Tomczak et al. 2017) surveys (only the integral term in Equation 7 is used). The cluster study of ) and some observations and simulations for the overall field galaxies are also shown (Keres et al. 2003;Berta et al. 2013;Sargent et al. 2013;Keating et al. 2016;Decarli et al. 2016;Scoville et al. 2017). At each redshift, denser environments have a higher molecular mass density than the less-dense regions, primarily because of a volume effect as they occupy a smaller volume of the total volume of a survey. Our ρ measurements in the lowest density bins peak around z ∼ 1-2 and their redshift evolution has a similar shape to that of Scoville et al. (2017) and some field studies. The peak in the ρ measurements in the highest density quantiles occurs at z < 1. This differential redshift evolution in low-and high-density environments is likely the result of the growth of the large-scale structure with cosmic time. The SFRD (scaled by a factor of 2000) from equation (15) of Madau & Dickinson (2014) is also shown, depicting a resemblance to the overall evolution of ρ.

.
In order to evaluate the integral in Equation 5, we rely on the stellar mass function (Φ(M * )) instead of the unknown Φ(M ISM ) and the stellar mass dependence of M ISM given in Equation 2 (also see Scoville et al. 2017). As seen in Equation 2, M ISM also depends on sSFR/sSFR MS and very weakly on the overdensity ∆. Since the majority of active galaxies lie on the main-7 For the local universe sample ∆z=±∆v/c=±0.0033 (Darvish et al. 2018) and at high-z, ∆z=±1.5σ ∆z/(1+z) , where σ ∆z/(1+z) is the median photo-z uncertainty at each redshift (Darvish et al. 2017) sequence (with starbursts making up only a few percent of the active population; e.g.; Rodighiero et al. 2011), we assume (sSFR/sSFR MS )=1 in Equation 2 when estimating ρ. Assuming a 5% contribution of starburst galaxies and an average (sSFR/sSFR MS )=10 for them changes the estimated ρ by only a factor of ∼ 1.03. We ignore the M ISM weak dependence on ∆ since changing ∆ by even a factor of 100 results in only a slight ρ increase of a factor of ∼ 1.15. We also do not include the contribution of passive galaxies in ρ. Given the above assumptions and knowing that Φ(M ISM , z)dM ISM =Φ(M * , z)dM * , Equa-tion 5 is simplified to: In Equation 7, we use the stellar mass function of Ilbert et al. (2013) for star-forming galaxies at z=0.2-4.0. The integral is performed at M * =10 10 to 10 12 M ⊙ because we only used M * 10 10 M ⊙ active galaxies in modelling M ISM (z, sSF R/sSF R MS , M * , ∆). Using M * =10 9 for the lower limit of integration changes the estimated ρ by a factor of ∼ 2.6. For the local universe stellar mass function, we also use that of Ilbert et al. (2013) at z=0.2-0.5.
The estimated f V and ρ values in different density quantiles and at different redshifts are given in Table  3. The uncertainties for f V are estimated by slightly modifying the width of the redshift slices within which projected surface densities are estimated (0.5 × and 2 × of the ∆z adopted in this work; see Section 3) and recalculating the new f V . The maximum difference between the new f V and the original f V is used as this volume-related uncertainty. The uncertainties of ρ have contributions from the uncertainties of f V and those of the parameters in Equation 2 and the parameters of the stellar mass function. Figure 5 shows the estimated molecular mass density as a function of redshift for different density quantiles. We also estimate it by directly using the star-forming stellar mass functions in densest (HD) and the least dense (LD) regions from the VIPERS (Davidzon et al. 2016) and ORELSE (Tomczak et al. 2017) surveys (only the integral term in Equation 7 is used). Note that the definition of the environment and how densest and lowest density regions are defined in these surveys are not necessarily the same as ours. We also show the cluster study of ) alongside some observational and computational studies for the overall field galaxies (Keres et al. 2003;Berta et al. 2013;Sargent et al. 2013;Keating et al. 2016;Decarli et al. 2016;Scoville et al. 2017).
At each redshift, denser environments have a higher molecular mass density than the less-dense regions. This is mainly due to a volume effect as the denser regions occupy a smaller volume of the total volume of a survey. Our ρ measurements in the lowest density bins (δ 1 and δ 2 ) peak around z ∼ 1-2 and the redshift evolution has a similar shape to that of Scoville et al. (2017) and some field studies as shown in Figure 5. However, the peak in the ρ measurements in the highest environmental density bins (δ 3 and δ 4 ) is shifted to lower redshift of z < 1. This is the result of the growth of the largescale structure with cosmic time as reflected in the z evolution of the measured fraction of volume occupied by densest regions in Table 3. Densest environments occupy a smaller regions of space at lower redshifts than similar environments at higher z due to mergers of structures extended over larger volumes at higher redshifts (see e.g.; Chiang et al. 2013) and late-time assembly of galaxy groups and clusters. Therefore, the decrease in the global molecular mass at z 2 is partly compensated by the decrease in the volume assigned to dense regions at z 2, causing the shift in the peak of the measured ρ in dense regions relative to that of the less dense environments. At z=0.65 and in densest regions (δ 4 ), ρ is maximized at 338.0 +109.0 −104.2 ×10 6 M ⊙ Mpc −3 and in the lowest density quantile (δ 1 ), it peaks at z=1.75 with a value of 28.2 +2.8 −2.9 ×10 6 M ⊙ Mpc −3 (Table 3). We also mention that the molecular mass densities between densest and the least dense regions based on Davidzon et al. (2016) and Tomczak et al. (2017) stellar mass functions coincide with the midpoint of our measurement (between δ1 and δ 4 ) and also that of Scoville et al. (2017). This further supports the accuracy of our ρ estimation in different environmental bins.
Interestingly, the redshift dependence of ρ for the lowest density quantiles (and the overall ρ) resembles the shape of the global star-formation rate density (SFRD) of the universe (Madau-Lilly plot; e.g.; Madau & Dickinson 2014 and the references therein). For example, our data points in δ 1 density quantile rise as (1+z) ∼3.4±1.3 and decline as (1+z) ∼−2.9 with increasing redshift. This is in agreement with rising (1+z) 2.7 and declining (1+z) −2.9 SFRD from Madau & Dickinson (2014) (labelled as 2000×SFRD in Figure 5). This tempting result indicates that the evolution of the cosmic SFRD is likely governed by the evolution in the cosmic gas content of galaxies, with minimal environmental effects.
Several studies have found that the shape of the cosmic SFRD in dense environments of clusters and protoclusters is similar to that of the field galaxies or likely has a peak at redshift of z 2, slightly higher than that in the general field (e.g.; see Clements et al. 2014;Kato et al. 2016 andsimulations of Chiang et al. 2017). The peak of the cosmic SFRD in dense environments seems to shift in an opposite sense compared with the shape of the cosmic molecular gas mass density in dense regions found in our study. This apparent difference might be due to a nontrivial and likely an environmental-dependent conversion factor between the global SFRD and that of the molecular mass density. Nonetheless, the details of such study are beyond the scope of this paper.

DISCUSSION
For M * 10 10 M ⊙ star-forming and starburst galaxies, we find no significant evidence for the environmental dependence of their molecular gas mass, gas mass fraction, and depletion timescale since z ∼ 3.5. This seems to be at odds with some studies showing that some environmentally driven processes such as ram-pressure is responsible for stripping the gas and dust content of galaxies in dense environments, resulting in a lower gas mass and gas fraction than their field counterparts (see references in Boselli & Gavazzi 2014).
Ram-pressure stripping is particularly effective in lowmass 10 9 M ⊙ galaxies (e.g.; Fillingham et al. 2016) as the gravitational bounding force of more massive systems might be strong enough to keep its gas content. However, there is evidence for ram-pressure stripping occurring in more massive galaxies as well (e.g.; Poggianti et al. 2017). Although observations of the HI deficiencies in cluster galaxies suggest that ram-pressure stripping is significant (e.g.; Cayatte et al. 1990;Gavazzi et al. 2005), it seems that the molecular gas is less vulnerable to stripping as it is much denser and the observations of the radial distribution of HI and H 2 gas in nearby galaxies show a much extended HI than H 2 , making molecu-lar gas more bound to the host galaxy and less prone to stripping. This is supported by the lack of environmental effects on the molecular gas content of galaxies as seen in some studies (e.g.; Kenney & Young 1989;Koyama et al. 2017). However, others find lower (e.g.; Fumagalli et al. 2009) molecular gas in denser environments, suggesting that ram pressure is still strong enough on stripping the H 2 gas. The study by Mok et al. (2016) find an excess of molecular gas in denser environments although the difference between H 2 mass in dense and less dense regions is not significant (within ∼ 2σ in best cases) according to their tables 2 and 3.
Since the estimation of molecular gas for our localuniverse sample relies on the HI gas, the lack of environmental dependence on the molecular gas content of our low-z sample is also indicative of the environmental independence of the atomic gas at z ∼ 0. This might be because of the relatively high mass of our sample galaxies to be affected by ram pressure and/or due to our selection as we are only investigating star-forming and starburst systems.
Selection biases can also play an important role. Gas content of galaxies might depend on redshift, stellar mass, and the sSFR of the galaxy relative to the main-sequence (e.g.; Scoville et al. 2017). If z, M * , or sSFR/sSFR MS distribution of galaxies in dense environments is (intrinsically or due to selections) different than the field samples, this automatically causes a bias which would likely lead to misinterpretation of the results. For example, there might be a correlation between massive systems and denser environments (e.g.; Bolzonella et al. 2010;Darvish et al. 2016). Many studies also found a lower sSFR in denser environments (especially at lower redshifts) due to a higher fraction of quiescent galaxies and/or a lower SFR of galaxies in denser environments (e.g.; Baldry et al. 2006;Kauffmann et al. 2004;Darvish et al. 2016). These automatically bias the denser environment samples to higher M * and lower sSFR/sSFR MS distributions, both leading to lower gas content for galaxies in denser environments. The role of selection biases is recently highlighted by Koyama et al. (2017) who found no environmental effects on molecular gas and star-formation efficiency at z ∼ 0 after controlling for potential biases.
As suggested by some studies, the environment seems to only control the fraction of quiescent/starforming galaxies (Peng et al. 2010;Sobral et al. 2011;Darvish et al. 2014Darvish et al. , 2016; i.e., dense environments increase the likelihood of a galaxy to become quenched, particularly at lower redshifts. However, as long as a galaxy is forming stars, its star-formation activity is not much affected by its host environment. This picture is consistent with our results as we deal with star-forming and starburst galaxies (on the main-sequence and above) whose molecular gas content (and the subsequent starformation activity) is primarily driven by processes (regardless of its physics) other than their local environment. One direct consequence of this scenario is that environmental quenching of galaxies (if any) in terms of gas removal or consumption is a fast process. Otherwise, it would leave its imprint as lower gas content for star-forming and active systems in dense environments.
We also note that if galaxies in dense environments happen to preferentially populate a region of the 3D M * -sSFR/sSFR MS -z space that is, for some reason, sparse in galaxies (in less dense bins), then the control sample selection procedure causes a bias, particularly if a large number of them are jettisoned. This is because the procedure (Section 4) would automatically discard these galaxies as they would not have a representative counterpart in the M * -sSFR/sSFR MS -z space in lower density bins.
We also found no evidence for a different depletion timescale in denser environments than the field to z ∼ 3.5 consistent with e.g. Koyama et al. (2017) at z ∼ 0 and Noble et al. (2017) at z ∼ 1.6. However, some studies have suggested that denser environments accelerate (decelerate) the consumption of molecular gas, leading to a shorter (longer) τ in denser regions (e.g.; Lee et al. 2017vs. Mok et al. 2016. However, as noted by Scoville et al. (2017), if the star-forming gas is in self-gravitating giant molecular clouds (GMCs), then the internal structure of the GMCs determines the physics of the star formation, and the star-forming gas does not know if it is in a less or a more massive galaxy or in our case, if its host galaxy is located in a less or a more dense environment. In other words, the local environment does not seem to influence the amount of M ISM that would go into forming those GMCs.
We note that although our samples cover a large dynamical range of environments, to fully understand the role of extreme environments on the molecular gas content of galaxies, a dedicated survey that targets a large sample of galaxies (e.g. dust-continuum observations with ALMA) in extremely dense core of confirmed structures at high-z is crucial.

SUMMARY
We investigate the role of environmental density on the molecular gas content of a large sample of M * 10 10 M ⊙ star-forming and starburst galaxies to z ∼ 3.5. Similar to Scoville et al. (2017), we derive the scaling relations for the molecular mass M ISM and SFR efficiency as a function of redshift (z), sSFR relative to the mainsequence (sSFR/sSFR MS ), stellar mass (M * ), and also galaxy overdensity (∆). We also investigate the redshift evolution of the global molecular mass density (ρ) in different environmental density quantiles. The key results from this paper are: 1. We find no dependence of the M ISM (subsequently, gas mass fraction f gas ) and star-formation efficiency (subsequently, depletion time τ ) scaling relations on galaxy overdensity. M ISM approximately increases as a power of 0.03 for overdensity ∆. The power term for the star-formation rate efficiency as a function of ∆ is 0.004.
3. At each redshift, denser environments have a higher molecular mass density than the less-dense regions. This is mainly due to a volume effect as the denser regions occupy a smaller volume of the total volume of a survey. Our ρ measurements in the lowest density bins (δ 1 and δ 2 ) peak around z ∼ 1-2, resembling the evolution of the global star-formation rate density. However, the peak in the global ρ is shifted to z < 1 for densest environmental bins (δ 3 and δ 4 ). The differential evolution in the peak of the global molecular mass density across different environments is likely the result of the large-scale structure growth with cosmic time. At z=0.65 and in densest regions (δ 4 ), ρ is maximized at 338.0 +109.0 −104.2 ×10 6 M ⊙ Mpc −3 and in the lowest density quantile (δ 1 ), it peaks at z=1.75 with a value of 28.2 +2.8 −2.9 ×10 6 M ⊙ Mpc −3 (see Table 3).

A. THE EFFECTS OF CONTROL SAMPLE SELECTION
The control sample selection presented in Section 4 might lead to some biases in interpreting the environmental dependence of the gas content of star-forming and active galaxies. For example, if galaxies in dense environments happen to preferentially populate a region of the 3D M * -sSFR/sSFR MS -z space that is, for some reason, sparse in galaxies (in less dense bins), then the control sample selection procedure would lead to a bias, particularly if a large number of them are thrown away in the control sample selection. This is because the procedure would automatically discard these galaxies as they would not have a representative counterpart in the M * -sSFR/sSFR MS -z space in less dense environments.
It is also likely that star-forming and starburst galaxies intrinsically present different distributions of stellar mass and sSFR/sSFR MS in different density bins, further biasing the study towards a lack of environmental dependence for the molecular gas. Here, we only control for redshift and present and compare the distribution of stellar mass and sSFR in different density quantiles before controlling them, as shown in Figures 6, 7, 8, and 9.
According to Figures 7 and 9 and given the K-S pvalues, the sSFR/sSFR MS distributions show no significant difference ( 1.1σ in all cases) in different environments for both the local and high-z samples. This is to some extend expected because our sample comprises only active star-forming and starburst galaxies. For these galaxies, a large number of studies including observations and simulations found no significant environmental dependence of their average SFR, sSFR, and the main-sequence over a broad redshift range (see e.g.; Peng et al. 2010 δ 2 δ 1 δ 3 δ 1 δ 4 δ 1 δ 3 δ 2 δ 4 δ 2 δ 4 δ 3 Cen 2014; Ricciardelli et al. 2014;Vogelsberger et al. 2014;Sobral et al. 2015;Duivenvoorden et al. 2016;Hung et al. 2016). Even for studies that found an environmental dependence of SFR for star-forming galaxies (e.g.; von der Linden et al. 2010;Vulcani et al. 2010;Patel et al. 2011;Haines et al. 2013;Tran et al. 2015;Erfanianfar et al. 2016;Darvish et al. 2017;Wang et al. 2018), the reduction in the mean SFR for star-forming and active systems in dense environments is often small (∼ 0.1-0.3 dex) and mostly applies to satellite galaxies. However, stellar mass distributions show statistically significant differences between the least dense and densest environmental bins, particularly for our high-z sample. The largest difference is at ∼ 2.1σ (local universe sample) and ∼ 4.3σ (high-z sample) levels. It is not clear to us whether this is due to an intrinsic difference between stellar masses in different environments or the result of sample selection (for instance how the ALMA sources were originally selected for observation). The potential intrinsic correlation between massive galaxies and dense environments are found in several studies (e.g.; Darvish et al. 2016;Kawinwanichakij et al. 2017) although others found no such relations (e.g.; von der Linden et al. 2010). Nonetheless, the presented results using stellar mass controlled samples should be considered with caution. We also note that even without controlling for M * , sSFR/sSFR MS , and z, our multivariable fits presented in Section 5.1 show no overdensity dependence for our sample of M * 10 10 M ⊙ starforming and starburst galaxies.