Resolved Molecular Gas and Star Formation Properties of the Strongly Lensed z = 2.26 Galaxy SDSS J0901+1814^∗

J0901 is lensed by a group of galaxies, which needs to be accounted for explicitly in order to reconstruct the galaxy's source-plane structure. Our lens model therefore comprises one component representing the group halo and others representing the group members. The former is described by an elliptical power-law density distribution, whose (spherical) convergence profile is given by

$$\begin{eqnarray}&&\kappa ({\boldsymbol{x}})=\displaystyle \frac{{b}^{2-\alpha }}{2| {\boldsymbol{x}}{| }^{2-\alpha }},\end{eqnarray} \tag{ 1 }$$

where b is the Einstein radius. The group members within two Einstein radii are represented by singular isothermal ellipsoids (SIEs) given by Equation (1) with α = 1. In this case, b not only represents the Einstein radius but also is related to the velocity dispersion σ_v by b ∝ ${\sigma }_{v}^{2}$.¹⁵

While the position and ellipticity of the group halo are allowed to vary, the group members' positions and ellipticities are fixed to the observed values. Additionally, a lognormal prior about the nominal Faber–Jackson relation (Faber & Jackson 1976) is placed on their velocity dispersions. For any two galaxies G₁ and G₂, Equation (1) and the Faber–Jackson relation give ${b}_{2}/{b}_{1}\propto {\sigma }_{v,2}^{2}/{\sigma }_{v,1}^{2}\propto \sqrt{{L}_{2}}/\sqrt{{L}_{1}}$, where L_i is the observed luminosity of G_i. Using mass as a proxy for luminosity, we set priors, noting that Gallazzi et al. (2006) find that the scatter in the logarithmic mass–velocity dispersion relation is ≈0.07 for early-type galaxies selected from the SDSS (Abazajian et al. 2004). We also note that the presence of a galaxy at the location of the southern image represents a unique challenge given its close proximity. Due to its small halo mass, fits with an SIE model are challenging since deflections due to that potential never reach zero. Since this is a smaller galaxy in a dense environment, its mass profile may be tidally truncated, and we therefore adopt a truncated, elliptical pseudo-Jaffe profile (Keeton 2001) to represent this component. The spherical convergence profile for this model is given by

$$\begin{eqnarray}&&\kappa ({\boldsymbol{x}})=\displaystyle \frac{b^{\prime} }{2}\left[{\left(| {\boldsymbol{x}}{| }^{2}+{s}^{2}\right)}^{-\tfrac{1}{2}}-{\left(| {\boldsymbol{x}}{| }^{2}+{a}^{2}\right)}^{-\tfrac{1}{2}}\right],\end{eqnarray} \tag{ 2 }$$

where s and a are the core and truncation radii, respectively. The truncated pseudo-Jaffe assumption allows us to explore truncated mass models but preserves more extended profile options in the limit that the truncation radius (a) approaches infinity. The best-fit lens model parameters for all components are listed in Table 4.

The data used to constrain the model consist of the HST F606W imaging and the integrated CO (3–2) intensity map (Figure 7). The pair of merging images composing the northern arc lie across a critical curve in the image plane and are more highly magnified than the southern and western images (Figure 8). A larger magnification can allow for a more detailed analysis, but only over the fraction of the source that has crossed the caustic. There is also a larger uncertainty associated with the source-plane reconstruction using the northern arc, as the magnification varies rapidly near the critical curve (Figure 8). For these reasons, we do not include the northern arc when constraining the lens model parameters or performing the source-plane reconstructions presented throughout.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In addition to optimizing the lens model parameters, we include a registration offset between these data sets (referenced to the CO (3–2) data). For each set of lens model parameters and registration offsets, a goodness-of-fit statistic is computed by multiplicatively combining the Bayesian evidence from the optical and radio. We use the framework described in Suyu et al. (2006), Vegetti & Koopmans (2009), and Tagore & Keeton (2014) to reconstruct the pixelated source distribution of J0901 in the source plane, as seen in each band. An irregular, adaptive source grid is used with priors on the sources' surface brightness in the form of curvature regularization; the Bayesian evidence is maximized at each step. After optimization, slight discrepancies between the optical data and the model remain. We add smoothly varying, nonparametric perturbations to the potential to compensate for limitations of the macro-model (Figure 8). These lens potential perturbations are at the 1%–2% level, which correspond to changes in the deflection angle of 100 mas or less. For the optical HST data, such changes are significant; however, because the beam size is ∼1'' in the radio bands, the effect on the CO data is negligible.

Lens modeling of interferometric maps is complicated by the imaging process, which does not conserve surface brightness, can be strongly affected by choices in mapping parameters (e.g., visibility weights), and yields noise that is correlated in the resulting image. All of these effects can potentially cause the lens model and source-plane reconstruction to diverge from reality. While a number of routines have been developed in recent years to constrain lens models using visibility data directly (e.g., Bussmann et al. 2012, 2013; Hezaveh et al. 2013, 2016; Rybak et al. 2015; Spilker et al. 2016; Dye et al. 2018), many rely on parametric source models, which are overly simplistic compared to the resolved observations we have for J0901. Recognizing that lens models derived from visibility data and from deconvolved maps are both fundamentally limited by incomplete sampling in the uv-plane, we prefer to exploit the well-resolved structure in our maps of J0901 to derive our lens model. We defer comparisons with source-plane reconstructions inferred from nonparametric visibility-based models to future work.

In order to account for the image-plane correlated noise in our lens modeling, we follow the noise scaling technique of Riechers et al. (2008). For an individual data set, we scale the noise (for input into the lens modeling code) by some factor greater than unity that could be determined and verified by comparing the statistical properties of noise residuals in areas where lensed features are present and absent. However, because we are comparing source reconstructions across various data sets with different noise properties, we fix the noise scaling. A large noise scaling factor allows the code to underfit the data in the formal reduced-χ² sense, since the code assumes that there is more noise in the data, which leads to a higher regularization strength. Qualitatively, this approach smooths the source over a larger physical scale, and the resulting source-plane beam is larger.

Our source-plane reconstructions yield a spatially varying synthesized beam/PSF. In Figure 9 we show a grid of the beam HWHMs overlaid on a contour plot of the source-plane reconstruction for the matched CO (3–2) integrated line map as an example of the variation in beam/PSF shape that results from delensing. Although the beam shape varies by a factor of a few over the entire reconstruction, the beam is smaller and more consistent in the direction of the emission for J0901. We therefore adopt surface-brightness-weighted average beams/PSFs when analyzing the spatial information for J0901; these have FWHMs of 02–03 (corresponding to physical scales of 1.7–2.6 kpc).

Zoom In Zoom Out Reset image size

The lens model uncertainties are explored via Markov chain Monte Carlo modeling for the CO (3–2) data only to save computational time. As the CO (3–2) moment map was the primary input used to constrain the lens model, this method accurately captures the uncertainties in lens model parameters. Magnification factors are then derived for the individual maps by delensing the emission for the distribution of model parameters. The magnification factor uncertainties thus take into account uncertainties in both the surface brightness of the source and the lens model parameters.

With the lens model optimized, we perform source reconstructions of the integrated line maps, the individual velocity channel maps, and the 2.2 μm continuum map. We present the natural-resolution source-plane reconstructions of the CO, Hα, and [N ii] lines for J0901 in Figures 10 and 11. In Table 5, we present the 50th percentile magnification factors and 68% confidence intervals derived using the "natural" resolution data and the magnification factors derived from the "matched" resolution data, for each image separately and in aggregate.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 5.  Magnification Factors

Transition	Map	North	South	West	Total
CO (1–0)	Natural	${10.2}_{-0.9}^{+1.2}$	${7.4}_{-0.5}^{+0.6}$	${5.3}_{-0.4}^{+0.4}$	${22.7}_{-1.5}^{+2.1}$
	Matched	${20.9}_{-2.9}^{+4.3}$	${15.1}_{-1.9}^{+2.7}$	${11.1}_{-1.5}^{+2.0}$	${47.2}_{-5.7}^{+8.8}$
CO (3–2)	Natural	${14.2}_{-1.6}^{+1.8}$	${10.4}_{-1.3}^{+1.4}$	${5.5}_{-0.7}^{+0.8}$	${30.1}_{-3.2}^{+3.7}$
	Matched	${14.1}_{-1.5}^{+1.8}$	${10.7}_{-1.1}^{+1.2}$	${5.7}_{-0.6}^{+0.7}$	${30.6}_{-2.9}^{+3.3}$
${\rm{H}}\alpha$	Natural	${11.8}_{-1.9}^{+2.2}$	${11.4}_{-1.3}^{+1.7}$	${6.3}_{-0.7}^{+0.9}$	${29.6}_{-3.1}^{+4.0}$
	Matched	${12.2}_{-1.7}^{+2.0}$	${11.9}_{-1.7}^{+1.9}$	${5.7}_{-0.7}^{+0.9}$	${29.9}_{-3.3}^{+4.3}$
[N ii]	Natural	${11.2}_{-2.1}^{+2.7}$	${11.5}_{-2.0}^{+3.3}$	${4.5}_{-0.8}^{+1.1}$	${27.2}_{-4.1}^{+5.8}$
	Matched	${9.7}_{-1.1}^{+1.2}$	${8.8}_{-1.0}^{+1.3}$	${3.5}_{-0.4}^{+0.5}$	${21.9}_{-2.1}^{+2.5}$
$2.2\,\mu {\rm{m}}$	Natural	${11.1}_{-3.7}^{+5.4}$	${16.7}_{-4.2}^{+9.5}$	${8.6}_{-1.9}^{+3.8}$	${37.1}_{-8.1}^{+16.1}$
	Matched	${8.9}_{-3.3}^{+4.6}$	${16.9}_{-5.1}^{+12.1}$	${7.9}_{-2.3}^{+5.3}$	${33.7}_{-8.9}^{+19.2}$

Download table as:  ASCIITypeset image

While the CO, Hα, and [N ii] lines all show two emission peaks in the southern arc in the image plane, those peaks do not correspond to one another across all lines. In the source-plane reconstructions, the CO peaks remain distinct, but the Hα and [N ii] peaks do not. The two peaks seen in the VLT maps are nearly aligned with the positions of the average dynamical center determined from the channelized source reconstructions (see Section 4.2) and are potentially multiple images of the same region within J0901 caused by a foreground member of the lensing group. However, the two peaks may also have disappeared on reconstruction because of the degree of regularization (i.e., the smoothness prior may have "won" over fitting the data because of noise or flaws in the lens model), and/or because the CO and HST data used to constrain the model may not have much power over the relatively small region encompassed by the two VLT peaks.

We also reconstruct J0901 in the source plane for the individual channel maps (Figure 12). The prominent velocity gradient observed in the image plane is also apparent in the reconstructed channel maps. The well-resolved and smooth velocity gradient seen in all lines suggests that J0901 is likely a disk galaxy, despite the two bright peaks seen in the integrated line maps. We extract the spectra from the reconstructed channel maps and compare the line profiles to the observed profiles from the image plane (Figure 13). Since the per-channel magnification factors were not computed to include the northern image, we use the sum of the southern and western observed spectra scaled by the mean per-channel magnification factor in order to understand what effects differential lensing might have on the line profile.¹⁶ We extracted the source-plane spectra in apertures defined by the S/N > 2 regions in the corresponding integrated line maps. We note that this method is not a perfect match to the procedure used to extract the image-plane spectra; a more perfect match would require delensing the image-plane aperture for each channel. Since the area occupied by a channel's emission varies with velocity (as expected, particularly when considering the variation in magnification factor), the aperture defined by the integrated line map reconstruction may miss some emission in individual channels. However, this method is adequate for revealing any dramatic or velocity-correlated differential lensing effects.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Differential lensing does not appear to strongly affect the shape of the line profile of J0901 in the southern and western images. Since it is the bright northern image that only captures a portion of J0901's source-plane structure (and thus only a portion of the velocity structure), one might suspect that any distortions of the line profile are most likely to appear in analyses that include the northern image. However, it is the northern image's CO spectral profile that shows the double-peaked structure typical of rotating disks (Figure 2), which is perhaps only hinted at in the combined spectrum of the southern and western images and their reconstruction (Figure 13). While the spatial structure of the least distorted western image best matches the source-plane reconstruction, as expected, the delensed Hα and [N ii] spectral lines appear to peak at redder wavelengths than seen in the observed spectrum of the western image. In addition, some of the internal structure of J0901 is multiply imaged within the southern arc due to a foreground lensing group member. We are therefore unable to firmly constrain the intrinsic profiles of the Hα and [N ii] spectral lines for J0901.

In order to estimate a gas mass for J0901, we use the magnification-corrected natural CO (1–0) line luminosity derived from all three images, obtaining ${M}_{\mathrm{gas}}=({1.6}_{-0.2}^{+0.3})\,\times {10}^{11}({\alpha }_{\mathrm{CO}}/4.6)\,{M}_{\odot }$ (Solomon & Barrett 1991). We use the Milky Way CO-to-H₂ conversion factor due to J0901's disk-like ordered rotation (Figure 12), but it is also the value favored by the Narayanan et al. (2012) continuous metallicity and surface-brightness-dependent version of the CO-to-H₂ conversion factor. The metallicity-dependent form of the CO-to-H₂ conversion factor presented in Genzel et al. (2015) and Tacconi et al. (2018) yields a slightly lower value of α_CO = 3.8 M_⊙ K⁻¹ km⁻¹ s pc⁻². However, the inferred gas mass is consistent with our Milky Way α_CO-derived mass within the uncertainties. We note that there is also some uncertainty on the metallicity of J0901 (see Section 4.4). As a sanity check on the ∼30% difference between the CO (1–0) and CO (3–2) lines' magnification factors, we also calculate M_gas using the CO (3–2) map and its corresponding magnification (corrected for excitation using our measured global r_3,1 without magnification correction) and find M_gas = (1.2 ± 0.3) × 10¹¹(α_CO/4.6) M_⊙; this value is consistent with the CO (1–0)-derived gas mass and therefore gives additional credibility to the difference in the two lines' magnification factors (at least for the natural-resolution images).

Following Scoville et al. (2016), we also use the 877 μm (observed frame) SMA continuum detection as an alternative probe of the gas mass. This method relies on the adoption of a dust temperature; Scoville et al. (2014, 2016) recommend against using dust temperatures derived from multiband SED fits (T_dust = 36 K in the case of J0901; Saintonge et al. 2013), since they are luminosity weighted and thus biased toward the hotter components of the ISM that do not make up the bulk of the mass, and instead recommend the adoption of T_dust = 25 K. Both values result in ∼2–3× lower ISM masses than the CO-derived gas masses (M_mol = (4.8 ± 1.3) × 10¹⁰ M_⊙ for T_dust = 36 K and M_mol = (7.1 ± 2.0) × 10¹⁰ M_⊙ for T_dust = 25 K, when corrected by the CO (3–2) "natural" magnification factor). These continuum-derived ISM masses suggest that a lower value of α_CO ∼ 1.4–2 would be more appropriate for J0901 (closer to values derived for low-metallicity systems, or to the canonical value used for local U/LIRGs). However, since we do not independently derive a magnification factor for the 877 μm continuum data because of its low angular resolution and S/N, there is some additional uncertainty in the continuum-derived ISM mass and implied CO-to-H₂ conversion factor.

Given the uncertainty in α_CO for J0901, we adopt the magnification-corrected, natural-resolution CO (1–0)-derived value of ${M}_{\mathrm{gas}}=({1.6}_{-0.2}^{+0.3})\times {10}^{11}({\alpha }_{\mathrm{CO}}/4.6)\,{M}_{\odot }$, carrying the uncertainty in α_CO as a free parameter. Even with conversion factor uncertainties, we note that the gas mass of J0901 is comparable to those of other galaxies selected at submillimeter wavelengths but larger than those of other UV-selected high-redshift galaxies (e.g., Riechers et al. 2010).

Adopting the dust mass from Saintonge et al. (2013), corrected to our CO (3–2) magnification factor, we obtain a dust-to-gas mass ratio of ${({4.7}_{-1.2}^{+1.4})\times {10}^{-3}({\alpha }_{\mathrm{CO}}/4.6)}^{-1}$ for J0901. This ratio is within the normal range for disk galaxies in the local universe (e.g., Draine et al. 2007) but is a bit low for those with the same metallicity (as seen for the high-redshift galaxies in Saintonge et al. 2013). However, the dust-to-gas mass ratio strongly depends on the assumed CO-to-H₂ conversion factor, as well as the properties of dust adopted by the Draine & Li (2007) dust models. Lower CO-to-H₂ conversion factors would increase the dust-to-gas mass ratio by a factor of ∼5, bringing it more in line with the dust-to-gas ratios of systems where authors tend to adopt those lower values (i.e., SMGs and U/LIRGs; e.g., Santini et al. 2010).

Using our new magnification factors and Hα measurements, we determine improved SFRs for J0901. We use the SFR scaling factor from Hao et al. (2011)/Murphy et al. (2011) (as compiled in Kennicutt & Evans 2012) scaled to a Kroupa (2001) initial mass function (IMF). We find ${\mathrm{SFR}}_{{\rm{H}}\alpha }={14.5}_{-1.7}^{+2.1}\,{{\rm{M}}}_{\odot }\,{\mathrm{yr}}^{-1}$ using the total L_Hα and native magnification factor without correction for obscuration. Hainline et al. (2009) measured the Hα and Hβ lines for two regions within J0901, finding extreme obscuration corrections from Hα/Hβ that would increase the SFR by a factor of ≳20. However, that ratio could have been affected by the coincidence of a skyline with the Hβ emission. Using the total infrared (TIR) luminosity (L_TIR from 8 to 1000 μm) derived from the Draine et al. (2007) fits to J0901's dust SED in Saintonge et al. (2013) (${L}_{\mathrm{TIR}}={1.80}_{-0.41}^{+0.42}\times {10}^{12}\,{L}_{\odot }$ assuming our new magnification factor for the native-resolution CO (3–2) data) and our choice in IMF yields ${\mathrm{SFR}}_{\mathrm{TIR}}={268}_{-61}^{+63}\,{{\rm{M}}}_{\odot }\,{\mathrm{yr}}^{-1}$, comparable to the expected value based on the Hβ extinction correction to Hα. Kennicutt & Evans (2012)/Kennicutt et al. (2009) also give an alternative method for correcting Hα to account for obscured star formation using the observed L_TIR, but this method yields a much smaller value of ${\mathrm{SFR}}_{{\rm{H}}\alpha }+\mathrm{TIR}\,={103}_{-20}^{+21}\,{{\rm{M}}}_{\odot }\,{\mathrm{yr}}^{-1}$ (where we have corrected the luminosities for the different magnification factors for Hα and TIR as above). This hybrid method for calculating obscured SFRs involves a number of assumptions that may not apply to galaxies in the early universe and was calibrated using galaxies with infrared luminosities lower than that of J0901 (albeit with similar L_TIR/L_Hα ratios and attenuation levels). We therefore adopt ${\mathrm{SFR}}_{\mathrm{TIR}}={268}_{-61}^{+63}\,{{\rm{M}}}_{\odot }\,{\mathrm{yr}}^{-1}$ for our subsequent analysis, since it likely accounts for the bulk of the star formation in J0901 and is not likely contaminated by significant emission from the AGN (Fadely et al. 2010).

The fraction of the total SFR that can be accounted for by our Hα measurements is consistent with the SFR_UV/SFR_IR derived in Saintonge et al. (2013): ${\mathrm{SFR}}_{{\rm{H}}\alpha }/{\mathrm{SFR}}_{\mathrm{TIR}}={0.054}_{-0.014}^{+0.015}$ versus SFR_UV/SFR_IR = 0.040 ± 0.007. Since J0901 is known to have an AGN (Hainline et al. 2009) on the basis of its high [N ii]/Hα line ratio and large Hα FWHM, it is possible that the Hα-determined SFR is contaminated by emission from the AGN; IFU observations of the Hα emission from the nuclear region of J0901 obtained using adaptive optics show signs of a broad low-level outflow once disk rotation is corrected for (Genzel et al. 2014). However, for the emission from both the disk and nucleus analyzed here, the Hα FWHM is no wider than one would expect based on single-Gaussian fits to the double-peaked CO line profiles (at least for the line profile derived from the sum of the three images). It seems likely that most of the Hα emission is due to star formation and that some emission from the AGN, near the systemic redshift, masks J0901's double-peaked profile (particularly given the slightly poorer ∼40 km s⁻¹ velocity resolution of the VLT data and ∼150 km s⁻¹ CO peak separations) but contributes only a small amount to the total Hα luminosity. Higher S/N would be necessary to do a pixel-by-pixel decomposition of the broad- and narrow-line emission components to correct for the Hα emission from the AGN, as done for the nucleus in Genzel et al. (2014).

If we rescale the stellar mass from Saintonge et al. (2013) to use the same Kroupa IMF that we assume for our SFR and apply our Hα-determined magnification factor, we find that J0901 has ${M}_{\star }=({9.5}_{-2.8}^{+3.8})\times {10}^{10}\,{{\rm{M}}}_{\odot }$. Combined with the TIR-derived SFR, J0901 has a specific SFR of $\mathrm{sSFR}={2.8}_{-1.1}^{+1.3}\,{\mathrm{Gyr}}^{-1}$. Since we have simply corrected the Saintonge et al. (2013) derived values by our new magnification factors (the CO and Hα magnification factors are very similar), choice of IMF, and TIR/SFR conversion factor, J0901 still falls along the star-forming MS (e.g., Noeske et al. 2007; Speagle et al. 2014), with an upward offset of just ${0.27}_{-0.16}^{+0.20}\,\mathrm{dex}$. We also compare J0901's sSFR to the bimodal MS and starburst (SB) populations parameterized in Sargent et al. (2012)/Rodighiero et al. (2011), who find an MS scatter of 0.188 dex and a second Gaussian peak for starbursts offset by $\mathrm{log}(\langle {\mathrm{sSFR}}_{\mathrm{SB}}\rangle /\langle {\mathrm{sSFR}}_{\mathrm{MS}}\rangle )=0.59$ with a 0.243 dex scatter. In this scheme, J0901 falls between the distributions for MS and starbursts at ${0.22}_{-0.16}^{+0.20}\,\mathrm{dex}$, but with considerable uncertainty. Based on these parameterizations of the MS, J0901 appears to be a massive but otherwise "normal" MS galaxy that falls a little to the high side of the sSFR distribution.

Using our delensed images, we can measure the physical size of J0901 and its dynamical mass. Despite the complications potentially introduced by the spatially varying resolution that results from the delensing, the size of J0901 is quite robust. Gaussian fits to the delensed integrated CO emission maps (without accounting for beam/resolution effects) are consistent for the two lines, with major- and minor-axis FWHMs of 11 ± 01 and 085 ± 005, respectively (position angle of 82° ± 7°). The VLT observations have slightly smaller and more elliptical delensed angular sizes, (10 ± 01) × (060 ± 002) for Hα and (068 ± 006) × (022 ± 002) for [N ii], at position angles similar to those of the CO lines. At these angular scales, the adopted beam/PSF values do not significantly affect the source sizes, and both the convolved and deconvolved (reported) source sizes are consistent within their uncertainties.

In order to estimate the dynamical mass, we analyze our delensed three-dimensional data using the Bayesian Markov chain Monte Carlo tool GalPaK^3D (Bouché et al. 2015), which constrains parametric fits to galaxy morphologies and dynamics while accounting for instrumentation-induced correlations in both the spatial and spectral directions. For the parametric model, we assume either a Gaussian or exponential intensity distribution originating from an inclined thick disk with a rotation profile of v(r) = v_circtan⁻¹(r/r_v) and intrinsic velocity dispersion σ_v. In Table 6, we list the best-fit parameters for both models, and the resulting dynamical mass estimates using M_dyn = 233.5 (2r_1/2)v_circ² (from the standard M_dyn = rv²/G with units of the dynamical mass, half-light radius, and circular velocity set to solar masses, parsecs, and kilometers per second, respectively). For the radius, we use twice the half-light radius since that is a reasonable approximation for the radius that encompasses 90% of the emission for both assumptions of Gaussian and exponential flux profiles. We fit dynamical models to the CO (1–0), CO (3–2), and Hα data for both the Gaussian and exponential flux distributions in order to estimate systematic uncertainties caused by model assumptions that may not accurately describe the underlying emission. Attempts to fit the native-resolution reconstruction of the [N ii] maps did not converge. We suspect that this failure is due to a combination of factors, including models that poorly describe the observed emission (which might be expected if the [N ii] emission is mostly associated with the central AGN), reconstructed velocity channels that are limited in number and do not fully trace the broad emission wings, and the lower S/N of these data. The fit to the reconstructed Hα map for the assumption of a Gaussian intensity distribution converges to circular velocities significantly larger than that of the other emission lines and flux profiles, likely for the same reasons that the [N ii] does not converge at all. The subunity reduced χ² values for both Hα fits are due to the small number of reconstructed velocity channels (11) and the large number of model parameters being fit (10).

Table 6.  Kinematic Fit Parameters

Model	Parameter	Transition
		CO (1–0)	CO (3–2)	Hα
Exponential disk	R.A.	09h01m223518	09h01m223533	09h01m223523
	Decl.	+18°14'314922	+18°14'314944	+18°14'314387
	r1/2	4.83 kpc	4.76 kpc	3.65 kpc
	i	36°	34°	23°
	P.A.	51°	47°	43°
	rv	0.09 kpc	0.67 kpc	0.56 kpc
	vcirc	188 km s−1	230 km s−1	345 km s−1
	σv	38 km s−1	39 km s−1	60 km s−1
	${\chi }_{\mathrm{red}}^{2}$	1.2	1.6	0.97
	Mdyn	0.8 × 1011 M⊙	1.2 × 1011 M⊙	2.0 × 1011 M⊙
Gaussian	R.A.	09h01m223515	09h01m223527	09h01m223519
	Dec.	+18°14'314887	+18°14'314925	+18°14'314506
	r1/2	4.07 kpc	3.95 kpc	3.14 kpc
	i	30°	27°	17°
	P.A.	51°	45°	42°
	rv	0.05 kpc	0.95 kpc	0.92 kpc
	vcirc	218 km s−1	306 km s−1	500 km s−1
	σv	39 km s−1	35 km s−1	57 km s−1
	${\chi }_{\mathrm{red}}^{2}$	1.1	1.5	0.86
	Mdyn	0.9 × 1011 M⊙	1.7 × 1011 M⊙	3.7 × 1010 M⊙

Note. Since GalPaK^3D does not produce meaningful uncertainties and the assumed models may not accurately reflect the underlying emission and dynamics of J0901, the best-fit values should be treated as approximate. As discussed in the text, we do not consider the fit of the Hα kinematics for a Gaussian flux profile to be credible.

Download table as:  ASCIITypeset image

Using the five consistent best-fit models for the three successfully fit lines, we calculate a mean dynamical center for J0901 of R.A. 09^h01^m223523 and decl. +18°14'314813. We then use the lens model to project the position of the dynamical center to the image plane; these positions are shown as black crosses in Figures 1, 3, 5, and 6. As the coordinates of the mean dynamical center are outside the (primary) lensing caustic, that position only appears in the southern and western images. For the southern image, the foreground member of the lensing group creates two subimages of the mean dynamical center position. As the two peaks of emission in the VLT Hα, [N ii], and continuum maps are nearly aligned with the image-plane positions of the average dynamical center, these peaks may correspond to multiple images of nuclear emission associated with the central AGN (higher angular resolution observations are necessary to confirm whether these peaks are multiple images or unrelated internal structures).

From these fits, J0901 appears to be consistent with a relatively face-on disk with a half-light radius of ∼4.25 kpc (consistent with sizes from the Gaussian fits we previously derived from the delensed integrated line maps). This size is consistent with what has been found for other star-forming galaxies with similar masses and redshifts (e.g., van der Wel et al. 2014). The circular velocity is somewhat degenerate with the source size and inclination angle, so the best-fit models either find higher circular velocities with lower inclination angles or lower circular velocities with large inclination angles. On average (neglecting the more questionable fit to the Hα data), we find v_circ ≈ 260 km s⁻¹ and i ≈ 30°. The Rhoads et al. (2014) measurement of v_circ = (120 ± 7)/sin(i)km s⁻¹ is consistent with our average best-fit circular velocity and inclination angle. Based on the models' best-fit circular velocities and velocity dispersions, the molecular gas kinematics appear to be consistent with other high-z disks (e.g., Tacconi et al. 2013), with v_circ/σ_v ∼ 6.

These models yield an average dynamical mass estimate of ∼1.3 × 10¹¹ M_⊙ (again, neglecting the likely unphysical fit to the Hα data for an assumed Gaussian intensity distribution). All of the five best-fit models' dynamical mass estimates are lower than the total baryonic mass of ${2.6}_{-0.3}^{+0.5}\times {10}^{11}\,{M}_{\odot }$ that we infer from our adopted gas and stellar masses. However, adopting a lower value of the CO-to-H₂ conversion factor significantly alleviates this tension, dropping the total baryonic mass to ${1.2}_{-0.3}^{+0.4}\times {10}^{11}\,{M}_{\odot }$ for α_CO = 0.8. Intermediate values of the CO-to-H₂ conversion factor (favored by metallicity-dependent models, for example) could also be possible if new constraints on the lensing of the stellar mass tracers yield larger magnification factors, or if the dynamical mass is evaluated out to a larger radius (than our assumed value of 2r_1/2) that captures more of the CO emission. Better models of the lensing potential, morphology, and dynamics of J0901 (from data with higher resolution and/or S/N, and/or models that more closely match the true flux distribution and kinematics) may also alleviate some of the tension with the baryonic mass estimates. Models of low-inclination systems are particularly sensitive to assumptions of azimuthal symmetry that may not be valid for J0901 or many rotating systems in the early universe; lower inclination angles (which would imply higher circular velocities) may also alleviate tensions between the baryonic and dynamical masses.

We examine the spatially resolved Schmidt–Kennicutt relation (Schmidt 1959; Kennicutt 1998) for J0901 using the Hα and CO maps smoothed to the same spatial resolution. We use the L_Hα–SFR conversion factor given in Kennicutt & Evans (2012), which assumes a Kroupa (2001) IMF. The Hα brightness has not been corrected for extinction. Properly accounting for spatially varying extinction can significantly affect the slope of the Schmidt–Kennicutt relation (Genzel et al. 2013), but our current dust continuum data lack the spatial resolution for us to effectively perform such a correction. A global correction for the extinction (as in Sharon et al. 2013) would simply offset the relation to higher SFR surface densities (discussed further below).

In order to fit the Schmidt–Kennicutt relation in J0901 to a power law, we follow the methodology of Blanc et al. (2009) and Leroy et al. (2013) and perform a Bayesian analysis, since standard orthogonal least-squares regression fits are biased by clipping of the molecular gas and star formation surface densities at a chosen significance level. While the full methodology is presented in Blanc et al. (2009) and Leroy et al. (2013), in short, we iteratively calculate the SFR surface density for a random sample (with replacement) of the observed pixelized molecular gas surface densities in J0901 for a grid of potential normalization factors (A), indices (n), and intrinsic scatter values (σ) in the equation

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{SFR}}}{1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{kpc}}^{-2}}=A{\left(\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{gas}}}{1000{M}_{\odot }{\mathrm{pc}}^{-2}}\right)}^{n}\times {10}^{{ \mathcal N }(0,\sigma )},\end{eqnarray} \tag{ 4 }$$

where ${ \mathcal N }(0,\sigma )$ is a normal distribution with mean zero and standard deviation σ. For each possible combination of Schmidt–Kennicutt relation parameters, we grid the resulting model values of Σ_SFR and Σ_gas and compare them to a grid of the measured values to calculate a χ² value. As in Leroy et al. (2013), we apply a 2σ cut in gas mass surface density before comparing the grids of the observed and model data points in order to define a clear y-axis; since this cut is applied after the data are simulated, it does not bias the selection of the best-fit model in the same way as more conventional linear fitting algorithms. For each of the three Schmidt–Kennicutt parameters, we fit polynomials to the shape of their χ² values (taking the minimum χ² along the complementary parameters' axes, collapsing the model grid to a distribution of χ² values for each parameter separately) and use the polynomials' minima as the best-fit values of the parameters. We then perform this comparison multiple times, each time removing a pixel at random, perturbing the grid on which we compare the source and model, and perturbing the emission for both tracers by both the additive statistical uncertainty (on a per-pixel basis) and the multiplicative flux calibration uncertainty (applied to all pixels). Our best-fit values for the Schmidt–Kennicutt relation and their uncertainties (both statistical and systematic) are given by the mean and standard deviation of the resulting distribution of fitted parameters.

In addition to the different denominator of Equation (4) (which must be accounted for in comparisons to other Schmidt–Kennicutt studies and is chosen to reduce the fitting covariance), our implementation of the algorithm differs from that of Blanc et al. (2009) and Leroy et al. (2013) in the following ways: (1) We randomly draw 10⁴ values of Σ_gas for calculating model values of Σ_SFR and allow repeats, but we only perform the iterative fitting routine for 100 perturbations of the model/source (due to computational/time limits). (2) Since our Σ_gas and Σ_SFR uncertainties are both dominated by measurement uncertainties, we additively perturb both the model gas and SFR surface densities. (3) We sample the Schmidt–Kennicutt parameters in Δlog(A) = 0.03 from −1.25 ≥ log(A) ≥ −0.5, Δn = 0.05 from 0.8 ≥ n ≥ 2.3, and Δσ = 0.015 from 0 ≥ σ ≥ 0.3. (4) We assume a 10% flux calibration uncertainty for the Hα and CO data. (5) Our minimum χ² curves are fit by a third-order polynomial rather than a second-order polynomial as in Leroy et al. (2013) in order to better fit our skewed χ² curves.

In Figures 19 and 20 we show the Schmidt–Kennicutt relation using the native-resolution integrated CO (1–0) and CO (3–2) maps. In order to determine whether differential lensing affects the observed Schmidt–Kennicutt relation in J0901, we analyze each image of J0901 separately and all three images combined. We also compare these results to a Schmidt–Kennicutt analysis using the matched maps (not shown) and the source-plane reconstructions of the matched maps for both CO lines (Figure 21). Table 7 lists the best-fit parameters of the Schmidt–Kennicutt relation in Equation (4). The power-law fits are roughly consistent with superlinear indices of n ≈ 1.5 for both CO transitions, with a mean value of $\bar{n}=1.54\pm 0.13$, although individual fits for the image-plane analyses range from n = 1.38 to 1.73.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We note that in Figures 19–22 we show much smaller Σ_SFR and Σ_gas bins in our Schmidt–Kennicutt plots (0.025 dex) than we use in the fitting analysis (0.05–0.2 dex, randomly assigned for each iteration of the Schmidt–Kennicutt fitting routine). These smaller bins show some structure in the pixel densities due to differences in the brightness distributions of our Hα and CO maps that are amplified by the beams/PSFs (causing neighboring pixels to be correlated and blending real source structure together, in some cases causing misaligned peaks of emission). The pixel density structures are particularly apparent in the analysis of the source-plane reconstructed images since correlations in the observed images can be warped and exaggerated by the delensing process. The larger bins we use in the fitting analysis (plus the Monte Carlo perturbations in the model realizations) do not show these structures, and therefore these structures do not bias our power-law fits. The numbers of correlated pixels and corresponding numbers of independent beams/PSFs used in the fits are listed in Table 7 as N_pix and N_ind, respectively.

Zoom In Zoom Out Reset image size

The smoothing and inner-radius uv clipping used to create the "matched" data set lower the best-fit Schmidt–Kennicutt index for the CO (1–0) line only (except for the southern image, where the index is unchanged). We suspect that the best-fit index for the CO (3–2) data is unchanged between the "native"- and "matched"-resolution maps because considerably less smoothing (and no uv clipping) is required to create its matched map. Therefore, while the CO (1–0) line might better represent the true distribution of the total molecular gas mass, one must compare the matched maps in order to make fair comparisons between how the choice of molecular gas tracer affects the Schmidt–Kennicutt relation.

We find significant differences between the Schmidt–Kennicutt indices determined from the image-plane data and from the source-plane reconstructions; the best-fit indices for the source-plane reconstructions are much lower with $\bar{n}=1.24\pm 0.02$. It is difficult to explain this difference since lensing conserves surface brightness. We do expect the distributions of pixels in the parameter space of SFR versus gas mass to change owing both to the existence of multiple images of the same region within the source and to the larger number of samples per region within the source given the uniform image-plane pixel size. These effects should be particularly strong for the northern and southern images since they cross critical curves where the magnification is highest. The western image is the least distorted of the three images, and although it suffers from the lowest S/N, it shows a similar distribution of pixels in the parameters space of SFR versus gas mass to the source-plane reconstruction: the highest concentration of pixels is at the highest SFR and gas mass surface densities. The lack of difference in the Schmidt–Kennicutt index between the three images despite these lensing effects is reassuring since it indicates that the Schmidt–Kennicutt relation does not change much between different regions within J0901 (at least at the spatial resolutions probed here). The similar pixel distribution for the western image and source-plane reconstruction is similarly reassuring since it indicates that our lens modeling is working correctly. However, neither of these points explains why the best-fit index differs between the image-plane and source-plane analyses.

We note that there is one significant methodological difference between the image-plane and source-plane Schmidt–Kennicutt analyses that does influence the best-fit index. The source-plane reconstructions for the Hα data are more compact than for the CO data, and thus the gas and SFR surface densities are largely uncorrelated beyond the edges of the Hα emission. The uncorrelated pixels, if left in the Schmidt–Kennicutt analysis, result in indices of n ∼ 2, which is even steeper than found in the image plane. Clipping the data at a fixed value of the gas mass or SFR surface density to remove uncorrelated data would both bias the fits and leave too few pixels to support a robust fit. Instead, we only include pixels that have S/N ≥ 2 for both the gas and SFR maps. Due to the spatially varying noise in the source-plane reconstructions, this significance cut does not result in a constant surface brightness cut. We rederived the fits with S/N cuts between 2σ and 4σ and found no difference in the best-fit values of the Schmidt–Kennicutt parameters, so we believe that this method is reliable for the source-plane reconstructions. Therefore, we do not think that this methodological difference causes the difference in the Schmidt–Kennicutt index between the image and source planes.

One potential way to resolve the difference in the Schmidt–Kennicutt index between the source and image planes would be to observe J0901 at higher angular resolution. While the fit to the western image is determined by a large number of pixels (∼800–1000), those pixels only correspond to a small number of independent resolution elements (≲10). If the western image is expected to better represent the image-plane structure and have a different distribution of pixels in the parameter space of SFR versus gas mass relative to the other images, then more high-S/N independent data points would be helpful. While our calculated uncertainties for the best-fit correlations fold in the effects of excluding individual pixels, they do not fold in the effects of excluding entire resolution elements. Therefore, the uncertainties quoted in Table 7 are likely underestimates of the true uncertainties, particularly for the western image.

We also see a significant difference in the observed scatter of the Schmidt–Kennicutt relation, σ, between the source-plane and image-plane data. While the uncertainty in the scatter is likely an underestimate as described above, we suspect that the difference in the scatter is an artifact of the lens modeling noise scaling and regularization. As discussed in Section 4.1, we scale up the assumed image noise during the lens modeling in order to allow the code to underfit the flux distributions and thus minimize the effects of the correlated noise in the interferometric data. Since this process increases the regularization strength, it effectively smooths the source-plane emission and thus increases the source-plane S/N while decreasing the scatter of the SFR and gas mass surface densities.

Lastly, we find no significant systematic differences in the best-fit indices between the matched data sets for the two CO lines using either the image-plane data or the source-plane reconstructions, except in the case of the southern image. However, given that the native-resolution CO (1–0) and CO (3–2) data produce significantly different indices (n ∼ 1.7 using the CO (1–0) data and n ∼ 1.4 for the CO (3–2) data) and that the source-plane reconstruction produces yet a different index (n ∼ 1.2), the potential differences between indices for the two CO lines are inconclusive. Since we apply a global excitation correction in determining Σ_gas based on our measured r_3,1 values, we find consistent offsets (A) in the linear fits between the two CO transitions, regardless of map (native, matched, or reconstructed) or lens-plane image (north, south, west, or total).

In Figure 22 we compare the results for J0901 to other high-redshift galaxies in which resolved pixel-by-pixel Schmidt–Kennicutt analyses have been performed. Only eight high-redshift galaxies besides J0901 have been analyzed on a pixel-by-pixel basis on the Schmidt–Kennicutt relationship: seven SMGs and one normal disk galaxy (Genzel et al. 2013; Sharon et al. 2013; Rawle et al. 2014; Hodge et al. 2015; Cañameras et al. 2017; Gómez et al. 2018; Tadaki et al. 2018; but see also Freundlich et al. 2013; Sharda et al. 2018). We also compare J0901 to the sample of local disk galaxies from Leroy et al. (2013). For these comparisons, we convert all measurements to the same (Kroupa 2001) IMF. The position of J0901 (or any galaxy) on the star formation relation strongly depends on the assumed value of the CO-to-H₂ conversion factor. While α_CO is still uncertain for high-redshift galaxies, different authors' choices in CO-to-H₂ conversion factors are often justified based on metallicity or dynamical arguments. Therefore, we do not correct gas measurements to the same α_CO, but instead show a horizontal bar to indicate how gas mass surface density measurements would scale for the range of possible conversion values (0.7 ≤ α_CO ≤ 4.6). For GN20, EGS 13011166, PLCK G244.8+54.9, and HATLAS J084933, where the CO measurements were not made for the J = 1–0 transition, we include additional factors in this rough systematic molecular gas uncertainty to account for the unknown gas excitation; we allow the ratios to the ground state to be as low as those found for the Milky Way (r_2,1 = 0.50, r_3,1 = 0.26, r_4,1 = 0.15, or r_7,1 = 0.015; Fixsen et al. 1999) and as high as thermalized (r_2,1 = r_3,1 = r_4,1 = r_7,1 = 1.0).

One significant source of uncertainty in the SFE for J0901 is that we do not have comparable spatial resolution tracers of obscured star formation. Our dust map from the SMA is not sufficiently resolved to map local variations of the dust column, and we do not have Hβ maps or enough resolved multiband optical/UV data to perform a spatially resolved SED analysis as in Genzel et al. (2013). In Figure 22, the contours are for star formation traced by Hα but scaled to account for the total SFR as inferred from L_TIR; this rescaling is analogous to applying a uniform global extinction correction (the same technique used to correct the SFR surface density for J14011 in Sharon et al. 2013). The obscured star formation as probed by the total infrared luminosity is significantly higher than the Hα-traced star formation and moves J0901 to higher SFEs, making J0901 significantly offset from the Schmidt–Kennicutt relation found for local normal disk galaxies. For comparison, we also show how the locus of points would change for different extinction corrections to the Hα-derived SFR surface density in Figure 22, including no extinction corrections, the total infrared-corrected Hα emission to measure the SFR as in Kennicutt & Evans (2012), and the global extinction value calculated from the Hα/Hβ ratio in Hainline et al. (2009) (using the standard assumption of case B recombination). Using the total infrared-corrected Hα emission to measure the SFR moves J0901 to higher SFEs, but not as high as using L_TIR alone. The global extinction value calculated from the Hα/Hβ ratio is highly uncertain since the Hβ line was coincident with a skyline; however, using this line ratio to correct for extinction (and thus obscured star formation) produces a large SFR, comparable to that calculated using the infrared.

Global extinction corrections preserve the index of the Schmidt–Kennicutt relationship, but different extinction laws and patchy/localized extinction could significantly change the correlation's slope. Genzel et al. (2013) find that the index of the Schmidt–Kennicutt relation for EGS 13011166 varied between 0.8 ≤ n ≤ 1.7 for the different extinction corrections they explore. Such extinction corrections are particularly challenging for starburst SMGs where nearly all of the star formation is expected to be obscured. For GN20, HLS0918, AzTEC-1, and J084933, their SFR surface densities are inferred from maps of the continuum emission near the peak of the dust SED (∼170 μm rest frame) scaled to their L_TIR-determined SFRs (Rawle et al. 2014; Hodge et al. 2015); a similar L_TIR-determined SFR scaling is used to infer Σ_SFR for G244, but they scale continuum emission from farther down the Rayleigh–Jeans tail of the dust SED (∼750 μm rest frame; Cañameras et al. 2017), which may better trace dust mass than the SFR. Accounting for additional obscured star formation in J0901 (or unobscured star formation in the case of SMGs) may change the index of the Schmidt–Kennicutt relation, but J0901 would remain at elevated SFE relative to the local relation (regardless of the choice of CO-to-H₂ conversion factor). Resolved dust maps or extinction corrections are necessary to more firmly pin down the index of the Schmidt–Kennicutt relationship for UV-bright high-redshift galaxies and for J0901 specifically.

A second source of uncertainty in the position of J0901 relative to the star formation law is that we do not correct for contaminating AGN emission. While Fadely et al. (2010) determine that the AGN in J0901 is not a significant contributor to its FIR luminosity, the AGN may contribute to the Hα luminosity used in our pixelized analysis. If the AGN is producing Hα emission in excess of that expected from star formation (Genzel et al. 2014 suggest that the AGN is responsible for ∼10% of the Hα emission), then some regions of J0901 would have an unexpectedly large SFE, which could either globally bias the Σ_SFR–Σ_gas relation to higher SFEs (e.g., if the AGN emission were uncorrelated with the molecular gas) or bias the fit to steeper slopes (if the AGN were fueled by molecular gas, the excess Hα emission could correspond to high molecular gas surface brightness). In order to determine whether there are regions in J0901 that might be affected by the AGN, we examined SFE as a function of r_3,1 and the [N ii]/Hα ratio, since gas fueling the AGN may be at higher density, in a higher excitation state, and/or shocked. While the distribution of SFE has a tail toward higher values, we found no significant correlation between SFE and r_3,1 or SFE and metallicity. Absent some additional tracer of AGN-affected Hα emission in J0901, we err on the side of using all pixels in the Schmidt–Kennicutt analysis. An alternative method for identifying and excluding Hα emission from the AGN would be to perform a velocity decomposition for each SINFONI pixel; the broader Hα line profile could be associated with the AGN rather than star formation, although a narrow-line AGN component might still masquerade as star formation. However, our current data lack the S/N (and likely the spatial resolution) to do such a decomposition.

Adopting α_CO = 4.6 M_⊙ (K km s⁻¹ pc²)⁻¹ and scaling to the IR-determined SFR, we find that J0901 appears slightly offset to higher SFEs relative to the "sequence of disks" (e.g., Daddi et al. 2010b; Genzel et al. 2010), while lower values of α_CO move J0901 to the "sequence of starbursts," in line with the scenario that the distinction between such sequences is at least partly a product of the assumed α_CO factors. If we assume that the CO-to-H₂ conversion factors are correct for all of the resolved high-redshift sources, J0901 appears to fall along or slightly below a track of high-redshift starbursts with a net index that is potentially superunity (Figure 22). However, given the variety of assumptions involved (extinction corrections, excitation corrections, and α_CO), J0901 and the other eight individual galaxies studied to date may lie within the normal scatter of SFEs.

Individual Schmidt–Kennicutt indices range from 1.0 ≲ n ≲ 2.0 for these high-redshift sources. Given that many of the sources explored here use low-J CO lines (CO (1–0) or CO (2–1)) and that we find that there is no clear excitation difference for J0901 up to CO (3–2), we do not think that the range of indices reflects the critical densities of different observed gas tracers. For the local disk galaxies in Leroy et al. (2013), there is also a wide range of Schmidt–Kennicutt indices (all mapped in the CO (2–1) line), and it is only the distribution of indices that peaks at n ∼ 1. Wei et al. (2010) also find a range of Schmidt–Kennicutt indices (n ∼ 1.6–1.9, mapped in the CO (1–0) line) for nearby low-mass E/S0 galaxies with a median index of n ∼ 1.2. While some of the apparent variation in local galaxies is due to other factors (like spatially varying CO-to-H₂ conversion factors), some of the scatter is real, and we suspect that this is also the case for high-redshift galaxies. Therefore, larger samples of high-redshift galaxies need to be analyzed in a uniform way if we are to say conclusively whether they have a different Schmidt–Kennicutt index or higher SFE, or if there are differences between galaxy populations (like starbursts versus normal galaxies).

While different choices of molecular gas tracers can complicate comparisons between studies of the Schmidt–Kennicutt relation, these tracers' density sensitivities are valuable tools for probing the underlying volumetric "Schmidt law" (SFR ∝ ${\rho }_{\mathrm{gas}}^{n}$; Schmidt 1959). In a series of hydrodynamic galaxy simulations with 3D non-LTE radiative transfer modeling, Narayanan et al. (2008) and Narayanan et al. (2011) demonstrate that the change in the Schmidt–Kennicutt index with CO rotational line (for the surface density or integrated versions of the relation) differs depending on the index of the underlying volumetric star formation law. They argue that the cold gas less directly involved in star formation will be underluminous in higher-excitation emission lines; therefore, while the intrinsic star formation relation using a cold gas tracer might have an index of n = 1.5, higher-excitation emission lines would trace less mass per unit of star formation, resulting in observed indices closer to n = 1.

Variation in the power-law index with gas tracer has been seen in the integrated form of the Schmidt–Kennicutt relation (e.g., Sanders et al. 1991; Yao et al. 2003; Gao & Solomon 2004; Narayanan et al. 2005; Bussmann et al. 2008; Graciá-Carpio et al. 2008; Bayet et al. 2009; Iono et al. 2009; Juneau et al. 2009; see Tacconi et al. 2013; Sharon et al. 2016), but these studies do not observe all tracers for the same galaxies, nor do they examine spatially resolved star formation properties. Greve et al. (2014) analyze the integrated CO and FIR properties for a comprehensive sample of local U/LIRGs and high-redshift SMGs (although not every CO line is detected in every galaxy) and find strong trends in the integrated form of the Schmidt–Kennicutt index with critical density. They find n ∼ 1 for CO rotational transitions ≲J_upper = 6 and decreasing indices for higher-excitation lines, a pattern that does not match the predictions of Narayanan et al. (2008). However, this discrepancy is not entirely robust: Kamenetzky et al. (2016) perform a similar analysis using largely the same sample and find near-unity indices for all CO lines up to CO (13–12), unless AGN host galaxies were included in the analysis. For high-redshift galaxies with surface density measurements, neither J0901 nor the (current) high-redshift sources exhibit the change in index with CO rotational transition predicted by Narayanan et al. (2011) for any of the underlying potential Schmidt laws. Given the variation in indices seen for local disk galaxies (Leroy et al. 2013), it seems likely that intrinsic variations between galaxies will mask the population average when only small numbers of observations are available at high redshift. Larger samples of galaxies will therefore need to be mapped in multiple gas tracers (in addition to efforts addressing the extinction corrections mentioned previously) in order to properly test the Narayanan et al. (2011) surface density predictions and determine the underlying volumetric star formation law.