Dense-gas tracers and carbon isotopes in five 2.5

The origin of the high SFR observed in high-z dusty star-forming galaxies is still unknown. Large fractions of dense molecular gas might provide part of the explanation, but there are few observational constraints on the amount of dense gas in high-z systems dominated by star formation. We present the results of our ALMA program targeting dense-gas tracers (HCN(5-4), HCO+(5-4), and HNC(5-4)) in 5 strongly lensed galaxies from the SPT SMG sample. We detected two of these lines (SNR>5) in SPT-125-47 at z=2.51 and tentatively detected all three (SNR~3) in SPT0551-50 at z=3.16. Since a significant fraction of our target lines is not detected, we developed a statistical method to derive unbiased mean properties taking into account both detections and non-detections. On average, the HCN(5-4) and HCO+(5-4) luminosities of our sources are a factor of ~1.7 fainter than expected, based on the local L'HCN(5-4)-LIR relation, but this offset corresponds to only ~2 sigma. We find that both the HCO+/HCN and HNC/HCN flux ratios are compatible with unity. The first ratio is expected for PDRs while the second is consistent with PDRs or XDRs and/or mid-IR pumping of HNC. Our sources are at the high end of the local relation between the star formation efficiency, determined using the LIR/[CI] and LIR/CO ratios, and the dense gas fraction, estimated using the HCN/[CI] and HCN/CO ratios. In SPT0125-47, we found that the velocity profiles of the lines tracing dense (HCN, HCO+) and lower-density (CO, [CI]) gas are similar. In addition to these lines, we obtained one robust and one tentative detection of 13CO(4-3) and found an average I12CO(4-3)/I13CO(4-3) flux ratio of 26.1$_{-3.5}^{+4.5}$, indicating a young but not pristine interstellar medium. We argue that the combination of large and slightly enriched gas reservoirs and high dense-gas fractions could explain the prodigious star formation in these systems.


Introduction
Traditionally, the molecular gas in high-redshift galaxies has been probed by observations of the rotational lines of CO (e.g., . Systematic CO line surveys have shown an increasing molecular gas fraction with redshift (e.g., Tacconi et al. 2010Tacconi et al. , 2013Saintonge et al. 2013;Carilli & Walter 2013;Dessauges-Zavadsky et al. 2015;Aravena et al. 2016;Keating et al. 2016). Other methods based on the galaxy dust content found similar results (Magdis et al. 2012;Scoville et al. 2014;Béthermin et al. 2015;Scoville et al. 2016;Schinnerer et al. 2016). These analyses suggest that high star forma-tion activity in massive galaxies observed at high redshift is related to their larger gas content. Furthermore, an increase in the star formation efficiency, that is, the star formation rate relative to the total gas mass (as traced by CO or dust), seems necessary to account for the prodigious star formation in the most extreme systems (e.g., Engel et al. 2010;Genzel et al. 2010;Daddi et al. 2010;Magdis et al. 2011;Tan et al. 2014;Hodge et al. 2015;Silverman et al. 2015).
In contrast, we have much less information about the amount of dense, actively star-forming gas in high-redshift galaxies. The dense gas is usually traced using the rotational lines of highdipole molecules such as HCN, which have H 2 critical densities 10 4 cm −3 (Shirley 2015). While ? argue that in some Galactic clouds the HCN(1-0) line traces gas densities similar to those traced by high-J CO lines (∼ 10 3 cm −3 ), the higher J transitions of HCN are still considered to be our best indicators of dense gas in galaxies.
As gas in the interstellar medium collapses to form stars, it passes though the density regime where the hydrogen is predominantly molecular and is dense enough to excite rotational transitions of HCN and HCO + into emission (Larson 1994). For stars to form, the density must be sufficiently high in the starforming cloud cores that self-gravity dominates over tidal shear. For much of the volume of a typical galaxy, the threshold density for collapse in spite of tidal shear and the threshold density for excitation of high dipole moment molecules into emission are approximately the same (Stark et al. 1989). The luminosity of the HCN line is then a good measure of the star formation rate. In the central kiloparsec of large galaxies, however, the density of the interstellar medium may be high enough to excite HCN but not high enough to resist the very much higher tidal shear in those regions of high differential rotation; the HCN line may then have a large beam-filling factor and be bright in emission from an extended, non-cloudy molecular gas, even though no star formation can take place.
Observations of Milky Way clumps in HCN(1-0), as well as in a variety of other dense gas tracers, by Wu et al. (2005Wu et al. ( , 2010 found a roughly linear relation between the SFR and the densegas mass (M dense ) consistent with the galaxy-integrated measurements. They arrive at the same interpretation as Gao & Solomon (2004), namely that the SFR of a galaxy scales linearly with the number of dense-gas "units" in a galaxy, with the dense-gas star formation efficiency being constant. Similarly, observations toward high visual extinction lines of sight by Lada et al. (2010Lada et al. ( , 2012 and Heiderman et al. (2010) support this interpretation. In contrast, Bigiel et al. (2016) showed that the increase in densegas fraction towards the center of M 51 leads to a decrease in the dense-gas star formation efficiency (∝ L IR /L' HCN ), since there will be fewer overdense regions able to gravitationally collapse if the average density is higher.
This interpretation, however, as well as the claimed linearity of the SFR-dense gas relation is debated in the literature. For example, Bussmann et al. (2008a) observe sub-linear IR-HCN(3-2) luminosity relations for nearby galaxies, while García-Burillo et al. (2012) found a slightly super-linear IR-HCN(1-0) luminosity relation. Both of these sets of observations can be reconciled with simulations (e.g., Narayanan et al. 2008bNarayanan et al. , 2011. Similarly, star formation models suggest that the SFR-dense gas relationship is primarily set by the gas density probability distribution function (PDF) in galaxies. As a result, a range of SFR-dense gas mass slopes may be expected, depending on the effective density of the tracer and the exact form of the gas density PDF (e.g., Krumholz & Thompson 2007;Narayanan et al. 2008b;Hopkins et al. 2013;?;Onus et al. 2018). In a high-pressure turbulent interstellar medium (ISM), the average gas density is expected to be higher at all physical scales with a larger fraction of the total gas mass residing at higher densities (> 10 4 cm −3 ). This type of high-pressure turbulent ISM is more likely to occur in vigorously star-forming regions than in more quiescent ISM conditions (Papadopoulos 2010). This is supported by many numerical simulations, which tend to predict a dramatic increase in the dense-gas fraction during major mergers causing a short (∼100 Myr) boost of the star formation efficiency (e.g., Renaud et al. 2014). However, a recent study by Fensch et al. (2017) suggests that this phenomenon could not be efficient in mergers of two gas-rich galaxies at high redshift. Observational constraints are thus essential to test this result. HCN, HCO + , and CS studies of local (ultra-)luminous infrared galaxies (ULIRGS) have found high gas fractions (e.g., Gao & Solomon 2004;García-Burillo et al. 2012;Zhang et al. 2014) -but so far this has not been firmly established at high redshift.
In this paper, we present deep ALMA observations of HCN(5-4) in a sample of five lensed dusty star-forming galaxies (DSFGs) between z = 2.5 and z = 4 from the South Pole Telescope (SPT) sample (Vieira et al. 2010(Vieira et al. , 2013. We chose to target HCN(5-4), since the HNC(4-3) line is not observable beyond z = 3.2 and we want to observe the same transition in all the sources of our sample. The transition has a typical effective excitation density of ∼10 6 cm −3 (Shirley 2015) and is therefore a bone fide tracer of the densest molecular gas in galaxies. In addition to HCN(5-4), and since they can be observed at the same time, our observations also targeted the HCO + (5-4) and HNC(5-4) lines, which are dense-gas tracers in their own right. However, these two lines are harder to interpret because of their more complex chemistry. HCO + , being an ion, is dependent on the ionization of the dense molecular clouds and is thus a less direct tracer of these high densities (Papadopoulos 2007). The line ratio between HCN and HNC is dependent on the physical conditions and varies from ∼1 in dense dark clouds (Hirota et al. 1998) to ∼10 −2 in hot and dense regions of star formation like Orion (Schilke et al. 1992). In addition to the above lines, we were able to simultaneously target the J = 4 − 3 transition of the in massive short-lived stars. Thus, detecting this line is a clue of previous star formation episodes in these galaxies (Hughes et al. 2008;Henkel et al. 2010).
In Sect. 2 we present our observations and our data reduction and line extraction methods. We then discuss the properties of dense-gas tracers in Sect. 3. In Sect. 4, we compare the properties of dense and lower density molecular gas tracers (CO, [CI]). Finally, we present and interpret our measurements of the CO isotopes in Sect. 5. In this paper, we assume a Planck Collaboration et al. (2016) cosmology and a Chabrier (2003) initial mass function (IMF).

Observations
This paper is based on ALMA cycle-4 observations (2016.1.00065.S, PI: Béthermin) of five DSFGs from the SPT sample (Vieira et al. 2013;Weiß et al. 2013;Strandet et al. 2016), which were selected for having bright apparent infrared luminosities (L IR ≥ 5 × 10 13 L ) and high-quality ancillary data. ALMA bands 3 and, when needed, 4 were tuned to the redshifted frequencies of the HCN(5-4) transition. When possible, spare spectral windows of the correlator were placed such as to cover 13 CO(4-3), HCO + (5-4), HNC(5-4), and [CI](1-0) (this latest line is only observable in SPT0125-47). Each spectral window covers 1.875 GHz and contains 240 channels (coarsest frequency domain mode resolution with an online spectral averaging by a factor of 16). The previously mentioned lines are covered by two contiguous spectral windows in one of the side bands.
The sensitivities of the observations were determined to reach 5 σ (10 σ for SPT0125-47) using the L' HCN -L IR of Zhang et al. (2014). Since our sources are gravitationally lensed and have an extension of ∼1 arcsec (Spilker et al. 2016), we requested the most compact configuration of the array to avoid spreading the flux of our sources over several synthesized beams and thus maximizing our chances to detect the integrated emission of our objects. The characteristics of our survey are summarized in Table 1. The observations were performed on December 29th, 2016 for SPT0125-47 and between January 7th, 2017 and January 9th, 2017 for the other sources. During our observations, the precipitable water vapor level (PWV) was between 2.8 and 5.5 mm.

Data reduction and extraction of the spectra
We analyzed our data with the CASA software (McMullin et al. 2007). They were initially calibrated by the standard ALMA pipeline at the observatory. In addition, we manually flagged some antennae and spectral windows with poor bandpass, phase, or amplitude calibration. In particular, the antenna DA48 had y and z position offsets larger than 2 cm in all the January 2017 observations and its data were flagged systematically. The overall quality of the rest of the data is very good.
The data were imaged using the CLEAN algorithm, and a natural weighting scheme was chosen to maximize the pointsource sensitivity. We used cleaning thresholds corresponding to 3 σ in a given channel using the map standard deviation around our source to estimate the noise. Since our lines are very broad and the signal-to-noise ratio (S/N) is low, the channels were again rebinned at the imaging stage. Because SPT0125-47 has several 5 σ detections, we only needed a rebinning by a factor of four. Similarly, since SPT0300-46 and SPT0551-50 have at least one line above 3 σ, we used a factor of six. For the remaining sources, we set the rebinning factor to eight.
We decided not to subtract the continuum directly in the uv plane. Instead, the continuum was subtracted later at the line extraction stage (see Sect. 2.3). This choice was motivated by the presence of numerous broad lines in the two contiguous spectral windows of interest. The only good area of the spectrum without line contamination is between the HCO + and the HNC line, but is too narrow to accurately constrain the slope of the continuum (see discussion in Appendix B). The imaging of the other side band, which is free of detected lines except for SPT0125-47, showed that, at the depth of our observations, the continuum varies significantly over 2×1.875 GHz. It is thus necessary to use a first-order subtraction of the continuum.
Before extracting the spectra, we checked that our sources are not significantly extended. We produced high S/N images by combining all the channels of the four spectral windows (mode multi-frequency synthesis). We fitted a two-dimensional Gaussian model and found the width of the Gaussian to be consistent with that of the synthesized beam. Since the continuum emission of our sources is compatible with a point source at our ALMA resolution, we thus extracted their spectrum from the ALMA cube at the centroid of this Gaussian model. Using this method, we implicitly assume that the size of the regions emitting the dense-gas lines is not much more extended than the continuum. The extracted spectrum corresponds to the two contiguous spectral windows in the side band of HCN(5-4) and is presented in Fig. 1.

Line extraction
We extracted the lines by fitting the following model to the whole spectra. We assumed that the line profiles are Gaussian. We set a positivity prior on the line fluxes and allowed the full width at half maximum (FWHM) to vary between 200 and 1000 km/s, which is typical for these types of sources (e.g., Aravena et al. 2016). We allowed a velocity offset up to 200 km/s compared with the expected center of the line based on the systemic redshift estimated from the band-3 spectral scan of CO. Finally, we assumed a first-order baseline for the continuum emission of the sources. All the lines ( 13 CO(4-3), HCN(5-4), HCO + (5-4), HNC(5-4), and CN(4-3)) and the continuum baseline are fitted simultaneously. This allows us to estimate the degeneracies between the baseline parameters and the line properties. In Appendix B, we explain why this method rather than the classical extraction in moment-zero maps was preferred.
We checked a posteriori from the residuals that these are fair assumptions. We found no significant feature in the residual spectra. We computed the Pearson correlation coefficient between two neighboring channels and found that the correlation is lower than 0.2 for all our five residual spectra. Consequently, we can use these residuals to estimate the uncertainties on the line properties derived with our fit using a bootstrapping method. We thus took our best-fit model and added the residuals after randomly reshuffling the channels. Then, we refitted the result and saved the obtained best-fit parameters. We repeated this 10 000 times. The error bars on the parameters are derived by computing the standard deviation of all these realizations. By construction, this method takes into account the degeneracies between parameters, and in particular between the baseline determination and the line fluxes. We estimated the signal-to-noise ratio (S/N) by dividing the best-fit peak flux density by its uncertainty. Because of the combined uncertainties on the peak flux Table 1. Summary of the characteristics of our observations. Each source was observed in a separate ALMA scheduling block (SB). Two similar SBs were used for SPT0551-50, which required a longer integration. The total observing time and the time on source are t obs and t on , respectively. PWV is the average precipitable water vapor level during our observations. The beam size in the  Table 2. Summary of the properties of our sources. We report the properties of lines of our detections and tentative detections (S/N≥3). Upper limits correspond to 3 σ. The blend between HNC(5-4) and CN(4-3) in SPT0551-50 can be considered as a tentative detection as justified in Fig. 1  Notes. (a) w e assumed the same FWHM and ∆v for all lines to fit these low S/N sources. (b) T he magnifications come from Spilker et al. (2016). We have no good lens model for SPT0551-50 and we thus assume the median magnification of the SPT sample of 6.3 (Spilker et al. 2016). (c) T he average sample luminosity is determined combining detections and non-detections using the method described in Sects. 2.5 and 2.6. (d) : The S/N of the blend of HCN(5-4) and CN(4-3) is computed using the total flux of the two lines. (e) T he 13 CO(4-3) line of SPT0125-47 is extracted assuming that the three lines have the same width and velocity offset as discussed in Sect. 2.3.

Fig. 1.
Continuum-subtracted spectra of our five sources (blue filled histograms). Two spectral windows were combined to produce them. The best fit of the baseline is subtracted from the data. The red solid line is our best-fit model and the black dashed line is the best-fit decomposition of the blend between HNC(5-4) and CN(4-3). Only lines with S/N>3 are shown. The S/N of the blend between HNC(5-4) and CN(4-3) in SPT0551-50 only has an S/N of 2.6 using our deblending method, which could be underestimated. However, its S/N determined in Appendix B using the moment-zero map is 3.3. We thus included the line in the plot. The three bottom panels show sources with a low S/N. Their spectra were rebinned by a factor of two in the figure to allow a better visualization of the faint lines. density and the width, the relative uncertainty on the line flux is slightly higher than on the peak.
Only SPT0125-47 and SPT0551-50 have sufficiently bright lines to fit them directly with reasonable constraints on the FWHM and the velocity offset for each line. For the other sources, even if there is systematically a positive signal at the position of the lines, the S/N is often lower than 3 and the width and velocity of the lines cannot be constrained. We thus assume in our fit that these two properties are similar for all the lines. The low-frequency wing of the SPT0125-47 13 CO(4-3) line is not perfectly fitted, when the width and the velocity offset of each line are independent parameters. The FWHM is smaller than for other lines: 252±45 versus 475±128 and 496±119 for HCN(5-4) and HCO + (5-4), respectively. However, this narrower profile is unlikely to be real (see Sect. 5.2) and might be caused by the noise. To determine the best fit of 13 CO(4-3), we thus performed another fit assuming that the three lines have the same width and velocity offset. The flux is then slightly higher but consistent at 1 σ with the previous value: 1.15±0.22 versus 0.92±0.15 Jy km/s.
The results of our line extraction are summarized in Table 2 and the best fit is shown in Fig. 1 (red solid line). When a line is not detected at ≥3 σ, we derived a 3-σ upper limit by summing the best-fit flux measured in the spectrum and three times the standard deviation of the flux in our multiple bootstrap realizations. Adding the signal present in the spectrum is crucial to obtain a reliable upper limit and is often forgotten in the literature. Otherwise, a source detected with an S/N of 2.9999 would have an ∼50 % probability of having a flux above the incorrectly computed 3-σ upper limit.
When HNC(5-4) and CN(4-3) are both present in the spectra, a special method is used to extract them, since they are severely blended. The deblending of these two lines is discussed in Sect. 2.6.
Previously, HCN(5-4), HCO + (5-4), and HNC(5-4) have been detected at high redshift only in quasars. Danielson et al. (2011) attempted to detect them in the eyelash star-forming galaxy, but obtained only upper limits. The other transitions of these three molecules have been found in star-forming galaxies only below z = 2.5. Our two detections in SPT0125-47 and our tentative detection in SPT0551-50 push the high-z observations of these dense-gas tracers in galaxies dominated by star formation to higher redshift. When we were finishing this paper, we became aware that Chentao Yang et al. were also working on HCN detections at z∼3 in NCv1.143 and G09v1.97 and they kindly provided us with their measurements (private communication 1 ).
The 13 CO molecule has rarely been detected in star-forming galaxies at high redshift. Danielson et al. (2013) found it in SMM J2135-0102 at z = 2.3. In addition, Weiß et al. (2013) reported two possible detections in the initial SPT SMG redshift surveys: SPT0529-54 at z = 3.36, and SPT0532-50 at z = 3.39. Our study doubles the number of 13 CO detections at high redshift.

Unbiased average luminosities and line ratios
Given the relatively low line detection rate of our sample, we seek sample-averaged line luminosities and ratios without being biased towards the brightest sources. To this end, we used the bootstrap analysis results presented in Sect. 2.3. In Fig. 2, we illustrate our method in the case of the average HCN(5-4) and HCO + (5-4) luminosities and the flux ratio between these two lines. The probability density distribution (PDF) of the HCN(5-4) (upper left panels) and HCO + (5-4) (upper right panel) luminosities of each source are represented as colored histograms. Concerning the >3 σ lines of SPT0125-47 and SPT0551-50, their PDFs exclude clearly the null hypothesis (L' = 0) as is expected for a (tentative) detection.
However, even for non-detections, the mode of the distribution is systematically above zero. Since these PDFs are quasi-Gaussian, the mode is very close to the best-fit value. This is the case for the eight 13 CO(4-3), HCN(5-4), or HCO + (5-4) lines present in our spectral windows. The case of HNC(5-4) and CN(4-3) is not discussed here but in Sect. 2.6, since they are blended and the interpretation is more complicated. If the luminosity of these lines were strictly zero, the probability to have such a result would be (1/2) 2 = 0.3 %, since there would be a 50 % chance for each individual PDF to have a negative mode. However, their flux is unlikely to be zero and could correspond to 1-3 σ, since we planned our observation to attempt detections and the lines should not be too far from the detection threshold. Negative values of the mode of the distribution would thus be rather unlikely and would request a 1-3 σnegative fluctuation of the noise in our observations. It is thus not so surprising that the modes of the PDFs of our non-detections are systematically positive. Indeed, the PDFs of such non-detections thus contain weak but potentially useful information, if we combine several objects. However, of course, these lines cannot be qualified individually as detections, since their PDFs do not exclude negative values and our measurements are thus compatible with a zero luminosity.
We determined the average luminosity of the sample using both detections and non-detections by combining their PDFs. Since our observations of each source are independent, we can assume that the measured luminosities are independent. The PDF of the sum of the luminosities of our sources can then be computed by convolving the PDFs of their individual luminosities, where p i (L i ) is the PDF of the luminosity of the i-th source. The average is then computed by dividing this sum by N, the number of sources. In practice, we do not really need to compute this convolution analytically. We can just randomly draw luminosities from our 10 000 bootstrap realizations for each of our sources and then sum them. We performed this operation 10 000 times to obtain the PDF of the average luminosity of the sample. In Fig. 2, the results are shown as gray filled histograms. As expected, the width of the PDF of the mean luminosity is significantly narrower than the PDFs of individual objects. We can also remark that the mean luminosity is clearly detected, since its PDF completely excludes zero. We do not apply a positivity prior to perform our fits of the spectra. Indeed, in the hypothetical case of a line with a null flux for all sources, the bootstrap realizations would have a zero flux or a positive flux depending on the realization of the noise. The average flux would thus be positive even if the flux is zero for all sources. Practically, the absence of positivity prior would have a minor impact, since the probability of a negative flux is small (<10%), except for HCN(5-4) in SPT0300-46 and 13 CO(4-3) in SPT0125-50. Indeed, our unbiased mean luminosity estimates Table 3. Summary of the CO and [CI] data used in this paper. All these data, except the new detection of [CI](1-0) in SPT0125-47, are ancillary and were extracted from the redshift-search program (Weiß et al. 2013) by Bothwell et al. (2016 (Bothwell et al. 2016) no data Notes. (a) Converted from the CO(3-2) flux using the mean flux ratio measured by Spilker et al. (2014) by stacking. Table 4. Line flux ratios measured in our sources. The mean ratio is computed combining all the sources for which these two lines were observed (see the description of the method in Sects. 2.5 and 2.6).
Line flux ratio SPT0103-45 SPT0125-47 SPT0125-50 SPT0300-46 SPT0551-50 Mean ratio Notes. For the ratios between two dense-gas tracers, we provided values only when the two lines are detected at more than 2 σ. Otherwise, the ratio is often compatible at 2 σ with both zero and infinity. Concerning the 12 CO/ 13 CO, only SPT0125-50 is observed and not detected and we consequently derived a lower limit on this ratio, since 12 CO is well detected. For the ratio between HCN(5-4) and CO(4-3) or [CI](1-0), we derived upper limits, since the denominator (CO or [CI]) is often detected. and the one derived with the positivity prior agree at better than 10 %, that is, 0.5 σ. This method is close to a stacking method, except that we first measure noisy fluxes and then average them instead of the contrary. However, since the lines of the various sources have different widths and the baselines are difficult to subtract reliably, we preferred this approach to a standard stacking. Indeed, the stacked lines would have had a very peculiar shape, since these are the sum of lines with various widths. It would have been impossible to use our procedure to fit the lines and the baseline simultaneously. Finally, as explained in Appendix B, flux measurements based on stacked moment maps are also unreliable, since the baseline subtraction is potentially affected by biases due to the crowded spectral windows. A stacking procedure should not be used in this case, since the systematic biases could be similar for all objects and this bias would stay roughly constant with the number of sources, while the noise would decrease in 1/ √ N. The measurement would thus be dominated by systematic effects and not the instrumental noise.
We used a very similar method to derive mean line flux ratios. We first computed the PDFs of the mean of the line fluxes from our bootstrap realizations using the same method as previously described. However, to avoid giving more importance to low-z sources in the computation of the mean, we weighted the sources by the square of the luminosity distance. We then computed the ratio between these mean line fluxes. The PDF of the ratio between two random variables is where p A (rX) and p B (X) are the PDF of the flux of the line A and B respectively, and r is the line ratio in luminosity. In practice, the uncertainties on the ratio are determined by drawing realizations from the PDF of the mean luminosity of each line. The results are presented in Fig. 2. The distribution is clearly asymmetric and we thus used the 16th and 84th percentile to produce the error bars in Table 2. With this method, the contribution of the sources to the mean ratio is weighted by their luminosity. We tried to compute first the PDFs of the line flux ratios of individual sources and then average them, but this failed to get meaningful results. For low-S/N sources, the distribution is very broad with very large negative and positive outliers. Indeed, the PDF of the line flux in the denominator is compatible with zero and there are thus realizations for which the ratios are tending to +∞ or −∞. In Appendix C, we present the simulation used to validate this method.
2.6. Deblending of HNC(5-4) and CN(4-3) and determination of unbiased mean ratios HNC(5-4) is blended with CN(4-3). To our knowledge, the only previous individual detection of this blend at high redshift was performed by Guélin et al. (2007) in the lensed quasar APM 08279+5255 (see also Riechers et al. 2007b about the detection of CN(3-2) in the Cloverleaf). It is also detected in the stacked spectrum of all the SPT SMG sources observed by ALMA in cycle 0 (Spilker et al. 2014). Guélin et al. (2007) found that HNC(5-4) is 1.74 times brighter than CN(4-3) by fitting simultaneously the profile of the two lines in their spectrum. This last source has a different nature than ours, but it shows that, even if HNC(5-4) might dominate the flux, the CN(4-3) contamination cannot be neglected. It is thus important to deblend the two lines in order to put a constraint on HNC(5-4).
To do so, we refitted the data with a slightly different method than the one described in Sect. 2.3. We fitted simultaneously all the lines in the spectrum including HNC(5-4) and CN(4-3) and forced all the lines to have the same width and velocity offset (including the sources with good S/N). The width and the velocity of the blended lines are thus strongly constrained by the other lines without being completely fixed. Since CN(4-3) and HNC(5-4) are not at the same exact frequency, these additional constraints are sufficient to extract information from the data. We also have to apply a positivity prior on the flux of HNC(5-4) and CN(4-3), since the degeneracies between the two fluxes tend to produce negative values to overfit noise patterns. Even under these assumptions, the uncertainties on the ratio for an individual source remain very high, with relative uncertainties higher than 50%. We thus derived the mean luminosities of the two lines and the mean line flux ratio between HNC and CN using the method presented in Sect. 2.5. We found a mean HNC(5-4)/CN(4-3) ratio of 1.60 +1.74 −0.82 . This value is close to the value found in the APM 08279+5255 quasar by Guélin et al. (2007).

CO and [CI] lines extracted from ancillary data
We want to compare our dense-gas tracers with lines from cold gas at lower density (CO, [CI]) found by our previous redshift search programs (Weiß et al. 2013;Strandet et al. 2016). The line fluxes of the [CI](1-0) line and the CO transitions covered by the redshift-search spectral scans were extracted in Bothwell et al. (2016). For most of our sources, the CO(4-3) line falls in a frequency window covered by the redshift-search data. We thus chose to use this transition available for most of our sources as the reference one for CO in our analysis. For SPT0125-47, which is at lower redshift, only CO(3-2) is available. We derived the expected CO(4-3) flux using the line ratio measured in the stacked spectrum of the SPT SMG sources derived by Spilker et al. (2014, I CO(4−3) /I CO(3−2) = 0.7). Finally, we detected [CI](1-0) in SPT0125-47 using our new ALMA data (see Appendix A). The ancillary data used in this paper are summarized in Table 3.

Scaling relation between the HCN and HCO + flux and the infrared luminosity
The relation between the dense-gas content, traced by the HCN(1-0) line, and the SFR, traced by L IR , was found to be linear in the local Universe by Gao & Solomon (2004). A linear correlation between L IR and L' HCN(1−0) is consistent with the simple physical picture in which the giant molecular clouds (GMCs), traced by HCN, convert a fixed fraction of their mass into stars before being disrupted (e.g., Faucher-Giguère et al. 2013; Grudić et al. 2018). However, further studies found a sublinear slope for the J=3-2 transition (Bussmann et al. 2008b;Juneau et al. 2009), which could have an impact on the interpretation of the physics of the star formation in infrared-luminous objects and active galactic nuclei (AGN, e.g., Narayanan et al. 2008a). More recently, after considering careful aperture corrections, Zhang et al. (2014) found a linear slope in the local Universe for the J=4-3 transition. The L IR -L' CO relations inferred for nearby galaxies have been found to be linear for J = 6 − 5 (Greve et al. 2014;Liu et al. 2015;Kamenetzky et al. 2016;Yang et al. 2017) and, possibly remain linear up to transitions as high as J = 12−11 (Liu et al. 2015). Extreme galaxies, however, such as the local ULIRG population and high-z starbursts, show sublinear L IR -L' CO relations for J = 7 − 6 and higher, due to large amounts of energy being injected, likely via mechanical heating, into a warm, dense, and non-star-forming ISM component (Greve et al. 2014;Kamenetzky et al. 2016). At high redshift, the L IR -L' HCN relation is poorly constrained due to the small number of detections, all of which probe the high luminosity regime. However, Riechers et al. (2007a) concluded, based on the couple of detections and upper limits obtained at that time, that the L IR /L' HCN ratio must be higher in high-z starbursts and quasars. This trend agrees with the theoretical model of, for example, Krumholz & Thompson (2007), which predicts a superlinear trend for HCN and HCO + at high infrared luminosity. Similarly, Narayanan et al. (2011) predicted a similar superlinear trend for CO at very high L IR . According to them, the median density of galaxies approaches the effective density of the molecular tracers, and the L IR -L mol relation will approach the L IR -M gas relation, which in both models has a 1.5 exponent and is thus superlinear. On average, in high L IR systems such as DSFGs, the median gas density would be much higher than in the lower luminosity systems investigated by Gao & Solomon (2004) for which they found a linear trend.
To test if our galaxies follow the relation of Zhang et al. (2014), we converted the L' HCN(5−4) of our objects into L'  assuming the line ratio measured in the local Universe by Mills et al. (2013, i.e., L' HCN(5−4) / L' HCN(4−3) = 0.73), which is also compatible with the measurements of Knudsen et al. (2007) in NGC253. We also corrected our sources for the magnification based on Spilker et al. (2016). The results are presented in Fig. 3  (upper panel). In addition to our detections and upper limits (black squares), we also show the average L' HCN(4−3) and L IR of our sample (red stars). We also compared the mean properties of our sample with the mean L' HCN(4−3) measured by Spilker et al. (2014) from a stacked spectrum of six SPT sources in a slightly lower redshift range than our targets. To place their data in the diagram, we divided the mean L'  and L IR of their six sources by their mean magnification, where a samplemedian magnification of µ=6.3 was assumed for sources with unknown magnification factor. Our new measurements are compatible with the previous results obtained by stacking.
All our objects have a small deficit of HCN luminosity compared to the Zhang et al. (2014) relation (black line), perhaps similar to the deficit in HCN(2-1) relative to infrared (IR) observed in high-z quasars by Riechers et al. (2007a) (after accounting for the AGN contribution to the IR luminosity). This explains a posteriori why we did not reach our initial goal of 5 σ detections in our four z>3 sources (see Sect. 2.1). From the Zhang et al. (2014) relation, we expect a mean L' HCN(5−4) of 11.1×10 8 K km/s pc 2 (after converting to the 5-4 transition) based on the infrared luminosity of our sources. We measured an average flux of (7.1±1.6)×10 8 K km/s pc 2 . This corresponds to a deficit by a factor of 1.6 (0.20 dex) corresponding to a 2.5 σ difference. However, Zhang et al. (2014) used the Sanders et al. (2003) method to derive L IR , while another method was used for the SPT SMG sample (Blain et al. 2003, model 1). We estimated the median factor between the two methods by fitting the photometric points used by Zhang et al. (2014) with the SPT SMG method. We found that their method derives luminosities that are a factor of 1.27 higher. When we take this into account, the tension increases up to a factor of 2.0 (0.30 dex). In contrast, the HCN(5-4) detection of NCv1.143 and the upper limit towards G09v1.97 (provided by C. Yang) are both, when converted to HCN(4-3), consistent with the observed local relation.
We checked if the deficit found with our sample could be explained by sample variance. The scatter around the Zhang et al. (2014) relation is ∼0.28 dex and we thus expect a 1-σ sample variance for five objects of 0.28 / √ 5 = 0.125 dex. If we combine this with our measurement uncertainties, the significance of the HCN deficit has thus only a ∼2 σ significance. It is therefore not possible to draw any firm conclusions, especially considering the uncertainties involved in converting our HCN(5-4) fluxes to HCN(4-3). Apart from intrinsic scatter in the 5-4/4-3 ratio, systematic effects could also be introduced by applying a locally-determined ratio to high-z starburst galaxies. More data will be needed in the future.
We performed a similar analysis for the HCO + . The conversion factor from the 5-4 to the 4-3 transition based on Mills et al.  Zhang et al. (2014) relation and the mean infrared luminosity of our sample, we expect a mean L' HCO + (5−4) of 6.9×10 8 K km/s pc 2 (after applying the correction to homogenize the L IR , see above) and we found (6.7±1.3)×10 8 K km/s pc 2 , which is in excellent agreement. Within the scatter, the HCO + detections provided by C. Yang are also in agreement with the relation.

The HCO + to HCN flux ratio
We compared the HCO + /HCN J = 5 − 4 flux ratios of our SPT sources (see Table 4) with HCO + /HCN ratios of low-and high-z galaxies from the literature (Fig. 4). The main advantage of using flux ratios is that they cancel out the magnification factor µ if the differential magnification (Hezaveh et al. 2012;Serjeant 2012) is negligible (see Sect. 4.3). While some of the HCO + /HCN ratios from the literature were for transitions other than J = 5 − 4, we proceeded to compare them with our ratios under the assumption that the spectral line energy distributions (SLEDs) of the HCO + and HCN are similar. In addition to the average value (1.00 +0.23 −0.19 ) of our sample, we put in our diagram only the sources that have at least a 3 σ signal at the position of each line (SPT0125-47 and SPT0551-50). For the other sources with <3 σ signal, both the numerator (HCO + flux) and the denominator (HCN flux) are compatible with zero at 3 σ (see Sect. 2.5). The PDF of their ratio is thus compatible with both zero and infinity at 3 σ and deriving upper or lower limits does not make sense. Nevertheless, we checked the average line ratio derived for these three other sources and found that it is compatible at 1 σ with the average value derived for the fives sources and the two individual measurements, but with very large uncertainties (0.98 +0.73 −0.40 ). Similar to the local star-forming samples of Graciá-Carpio et al. (2008) and Zhang et al. (2014), the mean ratio of our sample is compatible with unity. This is consistent with the ratio derived for the J=4-3 transitions using the stacked spectra of Spilker et al. (2014), but with improved uncertainties. The mean ratio of our sample is 2 σ higher than in the two high-z starforming galaxies of Oteo et al. (2016). NCv1.143 (Yang et al.) has similar values to these two objects, but G09v1.97 has a much higher value (∼2). However, all these high-z sources are in the intrinsic scatter of the local relation. As in the low-z Universe, the ratio varies significantly within the population of lensed DS-FGs. Braine et al. (2017) found that low-metallicity regions of local galaxies (<0.5 Z ) have a high HCO + /HCN flux ratio (∼2) instead of a ratio close to unity. Finding a unity ratio in our sources could be consistent with an already mature ISM at early cosmic times, but a larger statistical sample will be necessary to confirm or not this possibility.
Recently, Imanishi et al. (2016) and Izumi et al. (2016) found a HCO + /HCN flux ratio (∼0.5) for both the J = 3 − 2 and J = 4 − 3 transitions that is lower in local AGNs than in starforming galaxies 3 and proposed that this quantity could be used to identify AGNs. Izumi et al. (2016) discussed two explanations for the low HCO + /HCN: an enhanced abundance of HCN compared with HCO + coming from the complex chemical and radiative mechanisms involving these molecules in the neighborhood of an AGN or a systematically higher gas density around AGNs. We found a ratio close to unity, which is consistent with our objects being star-formation dominated. However, we should interpret this simple diagnostic with caution, since some AGNs of Imanishi et al. (2016) and the APM08279+5255 quasar at z = 3.9 have flux ratios close to unity.

The HNC to HCN flux ratio
Other important diagnostics can be performed using the isomer ratio between HNC and HCN. HNC traces gas of similar density to HCN, but the observed HNC/HCN flux ratio is close to unity in dark clouds and up to 100 times smaller in hot environments (Schilke et al. 1992;Hirota et al. 1998). In addition, Aalto et al. (2007) showed that the HNC/HCN flux ratio can be above unity in local starbursts, which is not intuitive because HNC should be more easily destroyed than HCN by the strong radiation fields and high temperatures in these objects. They proposed two possible explanations: HNC is excited by mid-IR pumping of its rotational levels or X-ray dissociation regions (XDRs) have an impact on the abundance of HNC. Since then, the case of Arp220 has been extensively investigated and more signs of HCN or HNC pumping have been identified (e.g., Aalto et al. 2015;Galametz et al. 2016). Finally, a similar scenario was also discussed in Weiß et al. (2007) to explain HCN luminosity in the high-redshift quasar APM08279+5255. Line flux ratio between HNC and HCN as a function of the infrared luminosity. The red star is the mean ratio found in our sources and the gray diamond is the mean ratio found by Spilker et al. (2014) using a stacking of the cycle-0 SPT SMG spectra. The black square is the upper limit determined for SPT0551-50. For the other sources, the individual ratio cannot be constrained, since both lines are too weak to derive an upper or a lower limit. The blue filled circles represent the local sample of Costagliola et al. (2011, 1-0 transition). The green triangles are the 3-2 transitions in the Arp220, Mrk231, and NGC4418 starbursts (Aalto et al. 2007). The purple crosses are the two z∼1.5 lensed star-forming galaxies of Oteo et al. (2016) and the yellow plus is the APM08279+5255 quasar (Guélin et al. 2007, see also Riechers et al. 2010 for the 6-5 transition). The orange downwards triangles are two H-ATLAS sources provided in Chentao Yang's thesis (private communication).
Because of the blending of HNC(5-4) with CN(4-3), we derived only a mean HNC/HCN ratio using the method described in Sect. 2.6, and found a ratio compatible with unity (1.03 +0.59 −0.39 ). This ratio is compatible with the one obtained by stacking of the HNC(4-3) and HCN(4-3) lines by Spilker et al. (2014, gray diamond in Fig. 5). Since the HCN(4-3) line is not blended, it is thus reassuring to find similar values. Our sources are at the border between the regime dominated by photodissociation regions (PDRs) and the domain, where XDR and/or mid-IR pumping are necessary to explain the line ratios. In Fig. 5, we compare the mean ratio found for our sample (derived using the method described in Sect. 2.6) with local and distant samples. The mean HNC/HCN ratio of our sample is a factor of approximately three above the local IRAM/30m sample of Costagliola et al. (2011), but lower by a factor of approximately two than the starbursts of Aalto et al. (2007). Our average measurements are only 1 σ above the measurements of APM08279+5255 by Guélin et al. (2007), SDP.9 of Oteo et al. (2017), and NCv1.143 in Yang's thesis, but an order of magnitude above the upper limit for SDP.11 (Oteo et al. 2017).

Dense gas versus lower density tracers (CO, [CI])
The origin of the strong star formation in the most extreme highredshift starbursts is a source of intense debates (e.g., Engel et al. 2010;Daddi et al. 2010;Hayward et al. 2011Hayward et al. , 2013Carilli & Walter 2013;Casey et al. 2014;Narayanan et al. 2015). Large gas reservoirs are not sufficient to explain the SFR of the most extreme systems and a temporary increase of the star formation efficiency, measured relative to the total gas mass as traced by CO(1-0) or CO(2-1), is necessary (e.g., Daddi et al. 2010;Genzel et al. 2010;Sargent et al. 2014). Gao et al. (2007) and Daddi et al. (2010) suggested that the increase of the star formation efficiency in these objects is linked to an increase of the densegas fraction (DGF). The dense gas fraction is thus one of the keys to understanding the nature of this type of sources. We note, however, that an increase in the dense-gas star formation efficiency is not required in this picture.

Dense gas fraction versus infrared luminosity
In this section we compare the flux ratio between HCN(5-4), which probes the dense gas, and [CI](1-0), which is thought to be a tracer of the bulk gas reservoir. For the three sources in our sample with [CI](1-0) measurements, we therefore adopt this ratio as a proxy for the DGF.
To build a reference sample in the local Universe, we combined the HCN(4-3) sample of Zhang et al. (2014) with the [CI] fluxes from the local Herschel/spectral and photometric imaging receiver (SPIRE) spectroscopic sample of Rosenberg et al. (2015). We adopted the aperture-corrected integrated [CI] fluxes published by Rosenberg et al. (2015), and, similarly, the aperture-corrected integrated HCN(4-3) fluxes published by Zhang et al. (2014). Finally, the HCN(4-3) fluxes were converted into HCN(5-4) using the conversion described in Sect. 3.1.
In the left panel of Fig. 6, we show the ratio between HCN(5-4) and [CI](1-0) versus L IR for our three sources observed in [Ci](1-0). Also shown is their mean flux ratio, as well as the upper limit derived from the stacked spectra of SPT sources (Spilker et al. 2014). Our sources are seen to be consistent with the stacked SPT value, but at the limit of the upper envelope of the local reference sample.
Given the lack of low-J CO transitions for some of our sources, we are not able to gauge the dense gas fraction using HCN(mid-J)/CO(low-J) ratio. Instead we opted for the J = 4−3 CO transition, since it is available for most of our sources, and examine the HCN(4-3)/CO(4-3) ratio. HCN(4-3) and CO(4-3) have similar upper level energies (E J /k B ∼ 40 − 55 K) and both trace dense gas, albeit HCN(4-3) has a ∼ 400× higher critical density than CO(4-3) (∼ 8 × 10 6 cm −3 vs. ∼ 10 4 cm −3 ). For the local reference sample we again adopt Zhang et al. (2014), and use their directly measured HCN(4-3) fluxes and the CO(4-3) fluxes from Rosenberg et al. (2015) to form the HCN(4-3)/CO(4-3) ratios. In the right panel of Fig. 6, we compare the HCN(4-3)/CO(4-3) flux ratios of our sources, where we have converted HCN(5-4) to HCN(4-3). There are not a lot of HCN(4-3)/CO(4-3) measurements of high-z sources available in the literature, and we therefore extended our comparison to HCN(5-4)/CO(5-4) and HCN(3-2)/CO(3-2) in order to facilitate a comparison with other high-z starbursts. Our mean HCN(4-3)/CO(4-3) ratio is in the scatter of the local values and compatible with the upper limit on G09v1.97 (Yang's thesis). In contrast, our sources are on average a factor of 1.6 lower than the stacking measurement of Spilker et al. (2014) on the full cycle-0 SPT SMG sample (∼2 σ difference) and a factor of approximately three below the high-z measurements of Oteo et al. (2016) in SDP.9, Danielson et al. (2011) in SMM J2135-0102, and NCv1.143 in Yang's thesis. Gao et al. (2007) found a correlation between L' HCN(1−0) / L' CO(1−0) and L IR / L' CO(1−0) , which they interpreted as a cor-  (1-0), tracing the dense-gas fraction, as a function of the ratio between the infrared and the [CI](1-0) luminosity, tracing the star formation efficiency. The data comes from the same references as in Fig. 6. The red solid line is a linear relation between the dense-gas fraction and the star formation efficiency (DGF ∝ SFE), which normalization has been set to match the mean value of our sample (red star).

A link between dense-gas fraction and star formation efficiency
relation between the DGF and the SFE (with respect to the total gas mass). In Fig. 7, we performed a similar diagnostic using [CI](1-0) instead of CO(1-0) and HCN(5-4) instead of HCN(1-0). Our SPT sources are consistent with the trend of DGF versus SFE found in our local reference sample described in Sect. 4.1.
A similar result is found if we use CO(4-3) instead of [CI](1-0) (see Fig. C.1 in appendix). This seems to indicate that the correlation between the DGF and the SFE is still valid for the most star-forming systems at high redshift and for the densest gas probed by the high-J transitions of HCN. However, an artificial correlation can appear in diagrams representing A/C versus B/C. In Appendix D, we confirm that it is not the case for our analysis using two different total gas mass tracers in the x and y axis.
The relation linking our local reference sample and the SPT sources suggests that DGF ∝ SFE (solid red line in Fig. 7). This result agrees with the suggestion by Gao et al. (2007) and Daddi et al. (2010) that the high SFEs in starbursts are directly connected to their DGF. The link between these two quantities is not surprising, if, as suggested by the various physical models cited previously, the amount of dense gas, rather than the total gas reservoir, drives the SFR. Unfortunately, our observations do not allow us to form a conclusion about the mechanism that causes these high DGFs. The sources in our sample are at the high end of the SFE and DGF distribution of our local reference sample. Their impressive SFRs are thus caused by a combination of large gas reservoirs ) coupled with a high DGF.

Similarity of the line profiles in SPT0125-47 and differential lensing
In SPT0125-47, dense-gas tracers are detected at more than 5 σ and it is thus possible to compare the line profile of HCN and HCO + with [CI] and CO (see Fig. 8). Our high-S/N [CI] detection and the ancillary CO(3-2) line detection from the redshift search have a similar asymmetric profile with a much broader redshifted tail. A similar asymmetry is found for HCN(5-4) and HCO + (5-4). This suggests that molecular regions from the lowest to the highest density are distributed in the same way across this object. The similarity of the profiles also suggests that differential lensing (Hezaveh et al. 2012;Serjeant 2012) should not be too strong in this object. Concerning the other objects of our sample, we do not have a sufficient S/N to be able to check the similarity of the profiles between the 12 CO lines and the dense-gas lines. However, as discussed in Gullberg et al. (2015), most of the SPT Fig. 8.
Comparison of the velocity profile of the dense-gas lines (HCN(5-4) in blue solid line and HCO + (5-4) in red dashed line, and other gas tracers (CO(3-2) in green dotted line and [CI](1-0) in yellow dot-dash line. We plotted the channel RMS uncertainties of HCN(5-4). HCO + (5-4) has similar channel uncertainties (not plotted). 12 CO and [CI] have a much better S/N and their uncertainties can be neglected. We normalized all the lines to have S ν dv = 500 Jy km/s. With this normalization, the peak flux of a rectangular line with a typical 500 km/s width is unity.
sources have for instance similar CO and [CII] profiles. Except if the dense-gas lines have a fundamentally different spatial distribution than lower density tracers, we have no good reason to expect a priori any strong differential lensing effect. We thus neglected this effect in this paper. The only way to measure the impact of the differential lensing on integrated fluxes would be to perform high-resolution imaging of dense-gas lines, but it is out of reach of ALMA in a reasonable amount of time.

Properties of the CO isotopic lines
5.1. Flux ratios between 13 CO and 12 CO As explained in the Introduction, the [ 12 C]/[ 13 C] abundance ratio has been proposed as a diagnostic of the evolutionary stage of a galaxy, since different nuclear reactions produce 12 C and 13 C; the former is produced via triple alpha nuclear processes in young massive stars, while the latter is produced in the CNO cycle in evolved asymptotic giant branch (AGB) stars (Wilson & Rood 1994).
A high [ 12 C]/[ 13 C] abundance ratio, therefore, could indicate a chemically young and largely unprocessed ISM (Hughes et al. 2008;Henkel et al. 2010). However, other physical mechanisms could result in a high [ 12 C]/[ 13 C] ratio, namely star formation from a top-heavy IMF in young starbursts (e.g., Romano et al. 2017;Zhang et al. 2018) -with the latter being the result of extreme cosmic-ray-dominated, star formation regions rather than a metal-poor gas. In intense starburst environments, the cosmic ray heating may be so severe that dense, deeply embedded molecular cores are significantly heated, resulting in a raise of the Jeans mass floor (Papadopoulos et al. 2011).
The [ 12 C]/[ 13 C] ratio is typically constrained indirectly from observations of CO, HCN, and HCO + and their 13 C isotopologues. This approach, however, is complicated by isotope fractionation where chemical reactions involving 13 C are energetically favored over the same reactions involving 12 C (Watson Fig. 9. Flux ratio between the 13 CO and 12 CO as a function of the infrared luminosity. The same transition of both lines is used to compute these ratios. Our sources are represented by black filled squares and their average properties by a red star. The two serendipitous detections obtained during the first redshift-search campaign in SPT0529-54 and SPT0532-50 (Weiß et al. 2013) are shown using orange pentagons. The gray diamond shows the average ratio measured by Spilker et al. (2014) using a stacking of the full redshift-search SPT sample. The local sample from Costagliola et al. (2011) is plotted with a blue filled circle. The yellow triangle and the green cross are the high-z measurements of Danielson et al. (2013) in SMM J2135-0102 and Henkel et al. (2010) in the Cloverleaf quasar. Fig. 10.
Comparison of the velocity profiles of the 13 CO(4-3) (blue solid line) and 12 CO(3-2) lines (green dotted line). The uncertainties on the 12 CO spectrum are not plotted, since they are negligible compared with the 13 CO ones. We normalized all the lines to have S ν dv = 500 Jy km/s. With this normalization, the peak flux of a line with a typical 500 km/s width is unity. The red dashed line is the 13 CO(3-2) profile measured by Zhang et al. (2018Zhang et al. ( ). et al. 1976Langer et al. 1984). The expectation from isotope fractionation is that [ 12 CO]/[ 13 CO] provides a lower limit to [ 12 C]/[ 13 C] (Langer et al. 1984;Tunnard & Greve 2016).
Furthermore, in intense far-UV environments, selective photodissociation can lead to an increase in [ 12 CO]/[ 13 CO], since 13 CO is readily destroyed given its low abundance, while 12 CO will be able to self-shield (Bally & Langer 1982). However, this effect is expected to play a role only in diffuse gas (A V ∼1) and it is thus doubtful whether selective photodissociation plays a significant role in dusty galaxies, where FUV light undergoes heavy extinction (Casoli et al. 1992;Papadopoulos et al. 2014).
Finally, it is non-trivial to infer 12 CO-to-13 CO abundance ratios from their line intensity ratios. This is because the 12 CO lines are optically thick, while the rarer 13 CO lines tend to be optically thin (e.g., Casoli et al. 1992;Aalto et al. 1995). Optical depth effects can thus complicate the picture, since in the case where 12 CO is optically thick the 12 CO/ 13 CO line intensity ratio will be an upper limit to the [ 12 CO]/[ 13 CO] abundance ratio.
In Fig. 9 (left), we compare our measurements of the 12 CO versus 13 CO ratio with the literature. We found a mean ratio of 26.1 +4.5 −3.5 , which is two times higher than the measurements of Spilker et al. (2014) using the stacked spectra of all SPT SMGs observed by ALMA during the cycle 0. Our sample was selected because of their high apparent luminosity and could be biased towards objects with a stronger contribution of the recent star formation to the ISM enrichment. It is therefore not surprising to find a lower 13 CO abundance, and thus a higher 12 CO/ 13 CO ratio, in our subsample. Weiß et al. (2013) reported a detection of 13 CO(4-3) in SPT0529-54 and a tentative detection in SPT0532-50 in the initial SPT SMG redshift-search program. The fluxes of these two lines, not published in the original paper, are 1.36±0.24 Jy km/s and 1.07±0.32 Jy km/s , respectively, while the 12 CO(4-3) fluxes are 8.32±0.34 Jy km/s and 14.89±0.41 Jy km/s. These two sources have much lower ratios than the new objects presented in this paper. This is not surprising, since these serendipitous detections are biased towards sources with particularly bright 13 CO(4-3) lines or with a positive fluctuation of the noise at the position of this line, while the 12 CO lines are brighter and thus systematically detected. Consequently, the 12 CO/ 13 CO ratios of these serendipitous detections are biased toward lower values. In Appendix E, we present a simulation illustrating this effect.
Our mean ratio corresponds to the upper envelope of the local sample of Costagliola et al. (2011) in the LIRG regime, but is lower than previous high-z measurements on individual objects of Henkel et al. (2010, Cloverleaf) and Danielson et al. (2013. Despite our sources having a deficit of 13 CO compared to the average values in the local Universe, their isotopic flux ratio is compatible with some local LIRGs. However, it is hard to know if it is really due to similar 13 CO abundances in local LIRGs and DSFGs or if the 13 CO abundance evolves and is compensated by optical depth effects. The presence of 13 CO shows that our objects are clearly not experiencing a first giant starburst fed by pristine gas.

Comparison of 13 CO and 12 CO line profiles in SPT0125-47
If the various parts of a galaxy have very different ages and ISM enrichment, we might observe very different 12 CO and 13 CO profiles. In particular, we expect a stronger enrichment of 13 CO in the center of galaxies, which harbor older stellar populations responsible for producing this isotope. Furthermore, levels of accretion of pristine cosmological gas onto the galaxy are lower in the central regions, resulting in less dilution of the enriched gas already there. In SPT0125-47, however, we found very similar 12 CO and 13 CO line profiles, suggesting the abundance ratio of the two molecules remains roughly constant throughout. In Fig. 10, we compare the velocity profiles of 13 CO(4-3) and 12 CO(3-2). In addition, we included the 13 CO(3-2) measurement of Zhang et al. (2018). Overall, the 12 CO and 13 CO profiles are rather similar. The two channels at the peak of 13 CO(4-3) are 1 σ and 2 σ above the 12 CO(3-2), but this is not the case for 12 CO(3-2). This small tension thus seems to be caused by noise. The similarity of the three profiles suggest that the maturity of the ISM across the different regions of SPT0125-47 is relatively homogeneous. These similar profiles are also a reassuring indication that there should be no strong differential magnification affecting this source.

Conclusion
In this paper, we presented the main results of our ALMA program targeting dense-gas lines (HCN(5-4), HCO + (5-4), and HNC(5-4)) in five strongly-lensed DSFGs at 2.5< z <4. We obtained in total two detections (S/N>5) in SPT0125-47 and four tentative detections (S/N∼3) of dense-gas lines in other sources. In addition, our observations yielded two detections (one of which is tentative), as well as two upper limits, of 13 CO(4-3). We also detected [CI](1-0) with an S/N of 36 in SPT0125-47. Finally, we developed a method to derive unbiased average properties of our sample to compare them with the low-and highredshift literature. From this information, we draw the following conclusions: -The average HCN(5-4) luminosity of our sample is formally a factor of ∼ 1.7 lower than what is expected from the local L' HCN -L IR linear relation of Zhang et al. (2014). However, if we take into account the sample variance, this tension decreases to 2 σ. Furthermore, we cannot rule out the introduction of systematic uncertainties in the conversion from HCN(5-4) to HCN(4-3) when comparing to the Zhang et al. (2014) sample. A similar trend, however, was reported in high-z quasars for high-J transitions by Riechers et al. (2007a) and is also expected by such theoretical models as that of Krumholz & Thompson (2007) for high-SFR systems. -Our sample has a mean I HCO + (5−4) /I HCN(5−4) flux ratio of 1.00 +0.23 −0.19 . This ratio is close to unity and is similar to what is found in local star-forming systems. Based on low-redshift observations, we would have expected a higher ratio for lowmetallicity galaxies or a lower ratio in AGN-dominated systems.
-After deblending the HNC(5-4) and CN(4-3) lines, we found an average I HNC(5−4) /I HCN(5−4) flux ratio of 1.03 +0.59 −0.39 , placing our objects at the border between the PDR-dominated regime corresponding to a ratio below unity and the highratio regime caused by XDR and/or mid-IR pumping of HNC.
-We inferred the HCN(5-4)/[CI](1-0) ratio for our sources and found it to be higher than what is measured in LIRGs and ULIRGs in the local Universe. To the extent that this ratio is a proxy of the dense-gas fraction, this finding would suggest that our SPT sources have a larger dense-gas fraction than local (U)LIRGs. In the diagram showing the dense-gas fraction versus the star formation efficiency relative to the bulk gas mass, our sources are located on the same linear relation as the local galaxies. -In SPT0125-47, we compared the profiles of our 5.5-σ detections of HCN(5-4) and HCO + (5-4), our 36-σ detection of [CI](1-0), and the CO(4-3) line from the original redshiftsearch observation. All these lines have the same asymmetric profile. This suggests that the low-and high-density molecular gas is distributed in the same way across this galaxy. This similarity of the line profiles also suggests the absence of significant differential lensing.
-Concerning the CO isotopologues, our sample has a mean I12 CO(4−3) /I13 CO(4−3) flux ratio of 26.1 +4.5 −3.5 , which is 30% lower than the previously observed high-z galaxies but 50% higher than typical local LIRGs. Even if the interpretation of this single line ratio is ambiguous, it tends to indicate that the ISM of our objects has already been enriched by intermediate-mass stars.
Even if the interpretation of these lines is non-trivial, our study allows us to draw a first preliminary portrait of the dense ISM in the SPT SMG sources. Previous works had already established that they had large gas reservoirs ), but we can now confirm that they have also a high DGF, as traced by HCN(5-4)/ [CI](1-0). A high DGF naturally leads to a higher star formation efficiency measured relative to the bulk gas mass, that is, a high L IR /L'  or, as in our case, a high L IR / L' [CI](1−0) . The combination of large reservoirs and high SFEs might be sufficient to explain the prodigious SFRs observed for the SPT sources. The ratios between dense-gas tracers are compatible with PDR-dominated objects and no evidence of AGN activity has been found so far.
Combining our sample with the previously observed highredshift galaxies dominated by star formation, we currently have dense-gas data for less than ten systems. Our analysis revealed some small tensions with other studies. However, their statistical significance is often below 2-σ. We need larger samples and higher S/N ratios to identify whether these differences are real and if they are linked to the physical properties of the galaxies. Our study shows that dense-gas tracers can be detected in ∼1 h with ALMA in highly-magnified objects and we can realistically hope to build a sample containing 20-30 objects in the coming years. Unfortunately, HCN is a bit fainter than expected from the observations at low redshift, but this bad news is totally compensated for by the fact that HCO + and HNC are at least as bright and can be detected by the same observations. Notes. (a) The S/N of our fitting method is determined by computing S peak / σ S peak . derived using our own method are slightly higher than the ones obtained with the moment maps. This is expected, since it takes into account the impact of the uncertainties on the baseline level and on the line width, while the moment-zero method assumes implicitly that the continuum is subtracted perfectly and uses a fixed velocity window to integrate the line flux. While they are slightly higher, the uncertainties provided by our fitting method are thus probably more realistic. Concerning the S/N ratios, the estimates based on the moment-zero maps are close to the ones derived using the ratio between the peak line flux density determined with our method and the uncertainty on it (S peak / σ S peak ). In SPT0125-47, the S/N and the line flux of [CI](1-0) based on the method presented in Appendix A is between the momentzero and the fitting method and we will keep this intermediate value for analysis. Indeed, since the [CI] flux is used to compute ratios with much more uncertain quantities, these small differences do not have any significant impact on our results. Our Gaussian fit and the moment-zero method agree at the 3% level. This confirms that the Gaussian approximation is a reasonable assumption for lines detected at ∼5 σ. In this case, the systematic effect induced by the non-Gaussian profile of [CI] is only 0.15 σ.

Appendix C: Validation of our method deriving mean luminosities and ratios
In Sects. 2.5 and 2.6, we presented a statistical approach to derive mean luminosities and line ratios without being biased towards detections. In this appendix, we present the simulation used to validate our method. We generated 500 simulated spectra at z∼3.2 of two spectral windows covering 13 CO(4-3), HCN(5-4), HCO + (5-4), HNC(5-4), and CN(4-3) (104.6-108.35 GHz). This redshift was chosen because it allows us to cover all the line at the same time and it is close to the mean redshift of the real sample. With these 500 spectra in hand, we reduced the statistical uncertainties on the derived mean quantities by a factor of ten compared to our real sample. This allows us to check if small systematic biases are affecting our approach.
For each source, we drew the fluxes of the 13 CO(4-3), HCN(5-4), HCO + (5-4), and HNC(5-4) lines from an uniform distribution between 0 Jy km/s and 0.3 Jy km/s, which is the typical range in our real sample if we leave out the exceptionally bright SPT0125-47. We used a maximum flux of 0.187 Jy km/s for CN(4-3) to be consistent with the mean observed flux ratio with HNC(5-4) (1.6, see Sect. 2.6). For all the lines of a given source, we used the same line width drawn from a uniform dis- Table B.1. Results of our simulation validating our method to derive mean luminosities and ratios (see Appendix C). For line fluxes, the true sample mean value X true is the mean of the 500 injected fluxes of our simulated sample. As justified in Sect. 2.5, the mean line flux ratios are derived by dividing the mean flux of a given line by the mean flux of another line. The fluxes measured using our line extraction tool are used to derive X mes . Xmes Xtrue must be close to unity if the method is not biased.
σ full ( Xmes ) Xtrue is the uncertainty on this ratio caused by the instrumental noise. Finally, by drawing 1000 subsamples of five sources from our 500 objects without withdrawal, we estimated the mean uncertainty on the measured-versus-true ratio for a subsample of five objects σ 5 ( Xmes ) Xtrue as is the case for our real observations.

Observable
X true X mes X mes X true σ full ( X mes ) X true σ 5 ( X mes ) X true Line fluxes in Jy km/s I13  0 tribution centered on 370 km/s and with a half width of 130 km/s as measured by Aravena et al. (2016) for low-J CO lines in the SPT SMG sample. Finally, we added a random Gaussian noise of 0.2 mJy per channel of 31.25 GHz (resolution of SPT0125-47 after rebinning). We extracted the lines in these simulated spectra using the method presented in Sect. 2.3. For HNC(5-4) and CN(4-3), following Sect. 2.6, we performed a second run assuming a positivity prior and no velocity offset to obtain a better deblending. The median S/N of the 13 CO(4-3), HCN(5-4), and HCO + (5-4) lines is 2.9, 2.8, and 3.1, respectively. Roughly half of our simulated sample is thus not even tentatively detected. In Table B.1, we compare the true fluxes injected into our simulation with the measured ones. The maximal difference between these two quantities is 6.3 %. This demonstrates that our method is reasonably accurate to derive the mean flux of a sample even if a large fraction of the objects is not detected. This difference of 6.3 % obtained for 13 CO(4-3) is significant at 4 σ, when we consider the uncertainty corresponding to the instrumental noise in our simulated spectra. There are thus residual systematic effects, but they remain small compared to the typ-  Ratio between HCN(4-3) and CO(4-3), tracing the densegas fraction, as a function of the ratio between the infrared and the CO(4-3) luminosity, tracing the star formation efficiency. The red solid line is a linear relation between the dense-gas fraction and the star formation efficiency (DGF ∝ SFE), which normalization has been set to match the mean value of our sample.
ical uncertainties on our measurements based on four or five sources ( 20%). In addition, we found that the measured fluxes of > 3 σ (tentative) detections are overestimated on average by 28 %, 41 %, and 13 % for 13 CO(4-3), HCN(5-4), and HCO + (5-4), respectively. Our statistical approach is thus clearly better than computing naively the mean flux of the detected sources.
Finally, we compared the flux ratio derived using our statistical method and the true ratio (Table B.1). We found no significant deviation. The small systematic effects affecting mean line fluxes could compensate for each other. The ratios including blended lines are also well recovered. This is a good demonstration of the accuracy of our approach including for blended lines (Sect. 2.6).