## Appendix A: Robustness of the joint-likelihood method using different line widths

One of the key assumptions is the width of the emission lines, which we fixed to 500km/s. Previous studies revealed a correlation of the total IR or line luminosity and the far-IR to mm line width that might originate from the regulation of gaseous disk rotation by gravity or (and) feedback from star formation or AGN (Bothwell et al. 2013; Goto & Toft 2015). The assumed line width generally matches the average of the DSFGs with ULIRG-HyLIRG luminosities in infrared, which is similar to the derived IR luminosity of our sample. However, observations also showed significant scatter in IR-luminous DSFGs. Our sample might be a typical example of the variety in the line width of luminous DSFGs, whose line FWHMs range from ~250km/s (HLS-2-1) to ~750km/s (HLS-3). In addition, the assumption of aa Gaussian line profile generally holds for most of our sources, but HLS-3, as described in Sect. 5.1, has a significant double-peaked feature in the detected emission line in band2.

The impact on the joint-analysis result of the mismatch between the real and assumed line width and profile is an uncertain prior. Thus, we made the following tests to verify if and how the results of this joint analysis might change with the assumption of different line widths. In addition to the default setting of a 500 km/s FWHM Gaussian line profile, we further performed the joint analysis with line profile FWHMs of 300km/s and 800 km/s, using the redshift and infrared luminosity derived from the fit with Béthermin et al. (2015a). We identified the redshift solutions of the four HLS sources with at least one detected line with these two different assumptions using the same method and criteria described as in Sect. 4.2. The results are listed in Table A.1.

From the results in Table A.1, we conclude that the redshift solutions from the joint analysis method are generally not sensitive to the assumptions on the emission line widths. For all sources but HLS-22, we find little variation in the redshift solutions using different line width assumptions. The differences of $\Delta z \sim 0.003$, as shown in Fig. A.1, mostly originate from the changes in the peak intensity of the emission lines, which could

lead to slight variations of $\chi^2(z)$. However, this small difference is still within the width of the emission lines, and will neither cause a false identification of the emission line, nor affect the analysis of the line fluxes and kinematics in Sect. 5.1 with the corresponding central frequency as an initial guess.

<table><caption>**Table A.1.** Redshift of HLS sources from the joint analysis with models of different line width</caption><thead><tr><th>Name</th><th>Z<sub>500km/s</sub></th><th>Z<sub>300km/s</sub></th><th>Z<sub>800km/s</sub></th></tr></thead><tbody><tr><td>HLS-2-1</td><td>5.241</td><td>5.242</td><td>5.243</td></tr><tr><td>HLS-2-2</td><td>5.129</td><td>5.129</td><td>5.132</td></tr><tr><td>HLS-3</td><td>3.123</td><td>3.122</td><td>3.125</td></tr><tr><td>HLS-22</td><td>3.036</td><td>2.436</td><td>3.034</td></tr></tbody></table>

The only case with a significant inconsistency in redshift from the test is HLS-22, where the procedure using the 300km/s line width strongly favors a redshift solution at $z = 2.436$. Given the frequency of the two emission lines detected with high S/N, we are confident about the redshift solution at $z = 3.036$ from the analysis with the model spectra with the 500km/s line width. Therefore, we checked the 300km/s model and the data at $z=2.436$ and found that the misidentification is caused by a strong noise spike at 100.64 GHz, as shown in Fig. A.2. The false identification suggests an increased sensitivity to narrow spikes in the spectrum with decreasing model line width. In our calculation of $\chi^2_{spec}(z)$ and correspondingly, the joint log-likelihood, their variation with redshift is dominated by the goodness of match between the model and data within the range of model line profiles. With a narrower line width in the model, the number of data points that dominate the variation of $\chi^2_{spec}(z)$ is smaller than in the cases with a wider line width. This makes the analysis with a narrow line width more sensitive to single spurious data points, such as the noise spike in the HLS-22 spectra, and leads to the misidentification in Fig. A.2. A less aggressive spectral binning along frequency and a preprocessing with sigma clipping could probably reduce this false identification in practice.

<figure><img src="image_11.png" /><figcaption>**Fig. A.1.** The output and comparison between the results of joint redshift analysis under different assumptions of line width. Upper row: Joint log-likelihood of the four HLS sources with model spectra with line widths of 300km/s, 500km/s, and 800km/s. Lower row: Comparison of the models of these line widths at the corresponding redshift solution and the observed source spectra.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. A.2.** Comparison of the observed spectra of HLS-22 and the model spectra with an FHWM of 300 km/s at the best redshift solution $z=2.436$. The feature identified as an emission line might be a missing glitch within the data.</figcaption></figure>

## Appendix B: Robustness of the joint-likelihood method with narrower frequency coverage

Millimeter spectral scans made by interferometers are widely used to blindly search for the emission lines from candidate high-redshift DSFGs, determine their spectroscopic redshift, and study the conditions of their cold ISM (Strandet et al. 2016; Fudamoto et al. 2017; Jin et al. 2019; Neri et al. 2020; Reuter et al. 2020). These spectral scan observations are designed to cover a continuous frequency range with several spectral setups. For ALMA and NOEMA, the default setup of the blind spectral scans at their current lowest-frequency band covers ~ 31 GHz. The earliest observations in 2018 of HLS sources blindly and continuously cover the spectra of HLS sources between 71 GHz and 102 GHz in NOEMA band1 with two setups, which follows this basic strategy of a blind redshift search. To test the joint-analysis method under more realistic conditions in large DSFG redshift survey projects, we applied

<table><caption>**Table B.1.** Redshift of the HLS sources from the joint analysis for the 31 GHz continuous spectra observed in 2018.</caption><thead><tr><td>Name</td><td>$Z_{true}$</td><td>$Z_{B15,2018}$</td><td>$Z_{MMPZ,2018}$</td></tr></thead><tbody><tr><td>HLS-2-1</td><td>5.241</td><td>5.241</td><td>5.242</td></tr><tr><td>HLS-2-2</td><td>5.129</td><td>5.131</td><td>5.129</td></tr><tr><td>HLS-3</td><td>3.123</td><td>N/A</td><td>N/A</td></tr><tr><td>HLS-22</td><td>3.036</td><td>4.380</td><td>1.690</td></tr></tbody></table>

the method to analyze these band1 spectra and compare their resulting redshifts with those in Sect. 4.2. The results from this analysis are presented in Table B.1 and Fig. B.1.

From the comparison of the redshift analysis using the early and full datasets, we find that the best redshift solutions of HLS-2-1 and HLS-2-2 remain stable, while the results of HLS-22 and HLS-3 are affected by the narrowed spectral coverage. This difference in robustness for different spectral coverage can be explained as follows. The joint-analysis method works equivalently to the automatic alignment and stacking of two or more lines in the spectra. If it is at the correct redshift and multiple lines are covered, this method can numerically boost the stacked S/N of emission lines, even if none of the single lines is detected with high significance. At a fixed coverage in frequency, we could expect that more CO lines are covered for sources at higher redshift. Taking our sample as an example, although the lines of HLS-2-1 and HLS-2-2 at $z \sim 5.2$ are only tentatively detected, their relatively high redshifts ensure that at least two CO and/or [CI] lines are covered by the spectral scan. In contrast, the narrowed spectral coverage leaves only one CO line in the spectral coverage of early observations of HLS-3 and HLS-22, leading to ambiguous redshift solutions regardless of whether the line is detected at high significance. With the comparison of the redshift robustness of these two groups of sources for different spectral coverage, we also emphasize that wide spectral coverage covering at least two strong molecular and/or atomic lines could be even more crucial in the redshift identification of DSFGs compared to reaching high sensitivity.

<figure><img src="image_12.png" /><figcaption>**Fig. B.1.** Result of the joint analysis of the four HLS sources with spectroscopic redshifts derived in Sect. 4.2, using only the 31 GHz NOEMA spectral scans observed in 2018. The first row shows the likelihood from the SED fittings and the joint log-likelihood of photometric and spectroscopic data, using the SED-fitting outputs with the Béthermin et al. (2015a) templates and MMPZ. The comparison of the observed spectra and the model spectra predicted by the IR luminosities from the two SED-fitting results are shown in the second and third row.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. A.2.** Comparison of the observed spectra of HLS-22 and the model spectra with an FHWM of 300 km/s at the best redshift solution $z=2.436$. The feature identified as an emission line might be a missing glitch within the data.</figcaption></figure>

## Appendix B: Robustness of the joint-likelihood method with narrower frequency coverage

Millimeter spectral scans made by interferometers are widely used to blindly search for the emission lines from candidate high-redshift DSFGs, determine their spectroscopic redshift, and study the conditions of their cold ISM (Strandet et al. 2016; Fudamoto et al. 2017; Jin et al. 2019; Neri et al. 2020; Reuter et al. 2020). These spectral scan observations are designed to cover a continuous frequency range with several spectral setups. For ALMA and NOEMA, the default setup of the blind spectral scans at their current lowest-frequency band covers ~ 31 GHz. The earliest observations in 2018 of HLS sources blindly and continuously cover the spectra of HLS sources between 71 GHz and 102 GHz in NOEMA band1 with two setups, which follows this basic strategy of a blind redshift search. To test the joint-analysis method under more realistic conditions in large DSFG redshift survey projects, we applied

<table><caption>**Table B.1.** Redshift of the HLS sources from the joint analysis for the 31 GHz continuous spectra observed in 2018.</caption><thead><tr><td>Name</td><td>$z_{true}$</td><td>$z_{B15,2018}$</td><td>$z_{MMPZ,2018}$</td></tr></thead><tbody><tr><td>HLS-2-1</td><td>5.241</td><td>5.241</td><td>5.242</td></tr><tr><td>HLS-2-2</td><td>5.129</td><td>5.131</td><td>5.129</td></tr><tr><td>HLS-3</td><td>3.123</td><td>N/A</td><td>N/A</td></tr><tr><td>HLS-22</td><td>3.036</td><td>4.380</td><td>1.690</td></tr></tbody></table>

the method to analyze these band1 spectra and compare their resulting redshifts with those in Sect. 4.2. The results from this analysis are presented in Table B.1 and Fig. B.1.

From the comparison of the redshift analysis using the early and full datasets, we find that the best redshift solutions of HLS-2-1 and HLS-2-2 remain stable, while the results of HLS-22 and HLS-3 are affected by the narrowed spectral coverage. This difference in robustness for different spectral coverage can be explained as follows. The joint-analysis method works equivalently to the automatic alignment and stacking of two or more lines in the spectra. If it is at the correct redshift and multiple lines are covered, this method can numerically boost the stacked S/N of emission lines, even if none of the single lines is detected with high significance. At a fixed coverage in frequency, we could expect that more CO lines are covered for sources at higher redshift. Taking our sample as an example, although the lines of HLS-2-1 and HLS-2-2 at $z \sim 5.2$ are only tentatively detected, their relatively high redshifts ensure that at least two CO and/or [CI] lines are covered by the spectral scan. In contrast, the narrowed spectral coverage leaves only one CO line in the spectral coverage of early observations of HLS-3 and HLS-22, leading to ambiguous redshift solutions regardless of whether the line is detected at high significance. With the comparison of the redshift robustness of these two groups of sources for different spectral coverage, we also emphasize that wide spectral coverage covering at least two strong molecular and/or atomic lines could be even more crucial in the redshift identification of DSFGs compared to reaching high sensitivity.

<figure><img src="image_13.png" /><figcaption>**Fig. B.1.** Result of the joint analysis of the four HLS sources with spectroscopic redshifts derived in Sect. 4.2, using only the 31 GHz NOEMA spectral scans observed in 2018. The first row shows the likelihood from the SED fittings and the joint log-likelihood of photometric and spectroscopic data, using the SED-fitting outputs with the Béthermin et al. (2015a) templates and MMPZ. The comparison of the observed spectra and the model spectra predicted by the IR luminosities from the two SED-fitting results are shown in the second and third row.</figcaption></figure>
## Appendix C: Cross-validation and difference between different SED modeling

The application of the joint-analysis framework to HLS sources largely relies on the current knowledge of the far-IR SED of high-z galaxies. However, although Herschel provides estimates of the mean far-IR SEDs and the redshift evolution of the main population of star-forming galaxies, these results are also limited by significant source confusion, especially in SPIRE data at longer wavelengths. Moreover, current studies revealed some DSFGs with an apparently low dust temperature (Jin et al. 2019), as well as a significant warm-dust contribution in some star-burst galaxies (Eisenhardt et al. 2012; Wu et al. 2012; Fan et al. 2016). These results suggest that a large variation in far-IR SEDs could exist in high-z DSFG populations.

Our choices of template are not free from these issues, and we therefore adjusted our joint-analysis framework to the results of two different far-IR SED templates and modeling framework, and we cross-validated the results of the two. The analysis with the Béthermin et al. (2015a) templates and MMPZ mostly shows consistent redshift solutions. This suggests the relative stability of the joint-analysis method with input information from different SED-fitting results.

However, some discrepancies on HLS-22 when using typical blind spectral scan conditions in Appendix B are also found, which caused us to make an additional check on its origin. From the comparison of the derived IR luminosity in Fig. 3, and the comparison of the model and data in Fig. 4 and Fig. B.1, we note that the estimated infrared luminosity and line fluxes from MMPZ are systematically lower than those from Béthermin et al. (2015a). At the correct redshift, the predicted line fluxes of Béthermin et al. (2015a) match the observed line fluxes better than MMPZ. In contrast, MMPZ generally returns more accurate photometric redshifts, especially on HLS-3, where the photometric redshift from Béthermin et al. (2015a) significantly deviates from the spectroscopic redshift. However, as indicated by the low dust temperatures of the HLS sample, it is possible that the properties of far-IR emission of these galaxies are not representative for high-redshift star-forming galaxies. Thus, we decided not to prefer a dust template and far-IR SED fitting in our framework, and we recommend a cross-validation of the redshift solutions from various methods in application.

In addition, the faint emission of HLS sources in the SPIRE bands introduces large uncertainties on the constraints of the source SEDs around the peak of the far-IR emission. As a result, this might contribute to the difference in IR luminosities derived from methods with different prior constraints (Casey 2020). These issues also further suggest the importance of matching observations at ALMA band 8-10 frequencies to properly reconstruct the far-IR SED, as well as to estimate the IR luminosity and star formation rate of high-z DSFGs selected by mm surveys.

## Appendix D: Impact of the scatter of the $L_{FIR}-L_{CO}$ scaling relation

The redshift from the joint analysis was derived based on the goodness of match between the emission line model and observed spectra. As mentioned in Sect. 4.2, in this approach, we predict the expected fluxes of spectral lines based on the $L_{FIR}$ from the SED template fitting using the best-fit scaling relation between $L_{FIR}-L_{line}$ from the literature (Greve et al. 2014; Liu et al. 2015; Valentino et al. 2018). However, these scaling relations are subject to substantial scatter of up to a factor of a few in observations. To test the possible impact of this scatter

<table><caption>**Table D.1.** Redshift of HLS sources from the joint analysis when using different numbers of offsets for all $L_{FIR}-L_{line}$ scaling relations.</caption><thead><tr><th>Name</th><th>z<sub>best,med</sub></th><th>z<sub>best,-1.0σ</sub></th><th>z<sub>best,-0.5σ</sub></th><th>z<sub>+0.5σ</sub></th><th>z<sub>+1.0σ</sub></th></tr></thead><tbody><tr><td>HLS-2-1</td><td>5.241</td><td>5.241</td><td>5.241</td><td>5.241</td><td>5.241</td></tr><tr><td>HLS-2-2</td><td>5.129</td><td>5.128</td><td>5.128</td><td>5.128</td><td>6.305</td></tr><tr><td>HLS-3</td><td>3.123</td><td>3.123</td><td>3.123</td><td>3.122</td><td>0.649</td></tr><tr><td>HLS-22</td><td>3.036</td><td>3.036</td><td>3.036</td><td>3.036</td><td>5.724</td></tr></tbody></table>

on our analysis and the robustness of the joint analysis method against them, we first checked the output best-redshift solution after adding a systematic offset to all $L_{FIR}-L_{line}$ scaling relations when generating the predicted spectrum models. In this test, we shifted the predicted CO/[CI] line fluxes by four different systematic offsets corresponding to ±0.5 and ±1.0 times of the $1\sigma$ scatter of the scaling relations. The exact values of the $1\sigma$ scatter for the considered lines are given in Table. 7.

The best-redshift solutions for HLS sources (except for HLS-4) after applying these four different offsets in $L_{FIR}-L_{line}$ conversion are listed in Table D.1. We find that all offsets in line flux, except for the $+1.0\sigma$, result in best-redshift solutions that are similar to the analysis using the median (zero offset) scaling relations. This suggests a good robustness of our joint-analysis method against the existing scatter of $L_{FIR}-L_{line}$ scaling relations in observations. As for the test with $+1.0\sigma$ offset, we further determined the reason for the discrepancy in the best-redshift solution. For HLS-22 and HLS-2-2, the application of the $+1.0\sigma$ offset leads to mismatches between the model and the data due to glitches or noise spikes in the spectra (see Fig. A.2 as an example, where the noise spike matches CO(5-4) at the best-redshift solution of HLS-22 here). For HLS-3, we find this unlikely low redshift solution after applying the $+1.0\sigma$ offset as the code assigns the only strong emission line in the spectral scan at 139.746 GHz to be CO(2-1). This suggests the stronger demand that the spectral scan be wide enough to cover more than one strong spectral line in the redshift confirmation of DSFGs with moderate redshift (i.e., $z \sim 2-3$).

To test the self-consistency of the joint-analysis method, we placed the measured CO line fluxes and far-IR luminosities of HLS sources at their best-redshift solutions (the $z_{best,med}$ in Table D.1) in the corresponding $L_{FIR}-L_{CO}$ diagrams to determine whether they follow the scaling relations used for line flux predictions. The results are shown in Fig. D.1. The plotted CO and far-IR luminosities were derived using the Gaussian-fit line fluxes and modified blackbody fitting (see Table. 9 and Table. 10). In addition to the average $L_{FIR}-L_{CO}$ correlations, we also show the samples from Cañamas et al. (2018) with multiple line transitions from the same sources as a comparison to our sources.

From Fig. D.1, we find that most of our sources fall within $+2\sigma$ of the scaling relation used in our analysis, which is also consistent with the regions occupied by bright submm galaxies in Cañamas et al. (2018). The only exception is HLS-3, which is also highlighted in Fig. D.1. In either the $L_{FIR}-L'_{CO(5-4)}$ or the $L_{FIR}-L'_{CO(4-3)}$ diagram, HLS-3 falls well below the scaling relation even when we consider the scatter of these scaling relations. However, we also note that its SPIRE photometry is one of the poorest in all of the four HLS sources. As the SPIRE bands close to the peak wavelength of SED predominantly constrain the IR luminosity, it is likely that the IR luminosity of HLS-3 is much less constrained than the other HLS sources, especially compared to HLS-2 and HLS-4.
## Appendix C: Cross-validation and difference between different SED modeling

The application of the joint-analysis framework to HLS sources largely relies on the current knowledge of the far-IR SED of high-z galaxies. However, although Herschel provides estimates of the mean far-IR SEDs and the redshift evolution of the main population of star-forming galaxies, these results are also limited by significant source confusion, especially in SPIRE data at longer wavelengths. Moreover, current studies revealed some DSFGs with an apparently low dust temperature (Jin et al. 2019), as well as a significant warm-dust contribution in some star-burst galaxies (Eisenhardt et al. 2012; Wu et al. 2012; Fan et al. 2016). These results suggest that a large variation in far-IR SEDs could exist in high-z DSFG populations.

Our choices of template are not free from these issues, and we therefore adjusted our joint-analysis framework to the results of two different far-IR SED templates and modeling framework, and we cross-validated the results of the two. The analysis with the Béthermin et al. (2015a) templates and MMPZ mostly shows consistent redshift solutions. This suggests the relative stability of the joint-analysis method with input information from different SED-fitting results.

However, some discrepancies on HLS-22 when using typical blind spectral scan conditions in Appendix B are also found, which caused us to make an additional check on its origin. From the comparison of the derived IR luminosity in Fig. 3, and the comparison of the model and data in Fig. 4 and Fig. B.1, we note that the estimated infrared luminosity and line fluxes from MMPZ are systematically lower than those from Béthermin et al. (2015a). At the correct redshift, the predicted line fluxes of Béthermin et al. (2015a) match the observed line fluxes better than MMPZ. In contrast, MMPZ generally returns more accurate photometric redshifts, especially on HLS-3, where the photometric redshift from Béthermin et al. (2015a) significantly deviates from the spectroscopic redshift. However, as indicated by the low dust temperatures of the HLS sample, it is possible that the properties of far-IR emission of these galaxies are not representative for high-redshift star-forming galaxies. Thus, we decided not to prefer a dust template and far-IR SED fitting in our framework, and we recommend a cross-validation of the redshift solutions from various methods in application.

In addition, the faint emission of HLS sources in the SPIRE bands introduces large uncertainties on the constraints of the source SEDs around the peak of the far-IR emission. As a result, this might contribute to the difference in IR luminosities derived from methods with different prior constraints (Casey 2020). These issues also further suggest the importance of matching observations at ALMA band 8-10 frequencies to properly reconstruct the far-IR SED, as well as to estimate the IR luminosity and star formation rate of high-z DSFGs selected by mm surveys.

## Appendix D: Impact of the scatter of the $L_{FIR}-L_{CO}$ scaling relation

The redshift from the joint analysis was derived based on the goodness of match between the emission line model and observed spectra. As mentioned in Sect. 4.2, in this approach, we predict the expected fluxes of spectral lines based on the $L_{FIR}$ from the SED template fitting using the best-fit scaling relation between $L_{FIR}-L_{line}$ from the literature (Greve et al. 2014; Liu et al. 2015; Valentino et al. 2018). However, these scaling relations are subject to substantial scatter of up to a factor of a few in observations. To test the possible impact of this scatter

<table><caption>**Table D.1.** Redshift of HLS sources from the joint analysis when using different numbers of offsets for all $L_{FIR}-L_{line}$ scaling relations.</caption><thead><tr><th>Name</th><th>z<sub>best,med</sub></th><th>z<sub>best,-1.0σ</sub></th><th>z<sub>best,-0.5σ</sub></th><th>z<sub>+0.5σ</sub></th><th>z<sub>+1.0σ</sub></th></tr></thead><tbody><tr><td>HLS-2-1</td><td>5.241</td><td>5.241</td><td>5.241</td><td>5.241</td><td>5.241</td></tr><tr><td>HLS-2-2</td><td>5.129</td><td>5.128</td><td>5.128</td><td>5.128</td><td>6.305</td></tr><tr><td>HLS-3</td><td>3.123</td><td>3.123</td><td>3.123</td><td>3.122</td><td>0.649</td></tr><tr><td>HLS-22</td><td>3.036</td><td>3.036</td><td>3.036</td><td>3.036</td><td>5.724</td></tr></tbody></table>

on our analysis and the robustness of the joint analysis method against them, we first checked the output best-redshift solution after adding a systematic offset to all $L_{FIR}-L_{line}$ scaling relations when generating the predicted spectrum models. In this test, we shifted the predicted CO/[CI] line fluxes by four different systematic offsets corresponding to ±0.5 and ±1.0 times of the $1\sigma$ scatter of the scaling relations. The exact values of the $1\sigma$ scatter for the considered lines are given in Table. 7.

The best-redshift solutions for HLS sources (except for HLS-4) after applying these four different offsets in $L_{FIR}-L_{line}$ conversion are listed in Table D.1. We find that all offsets in line flux, except for the $+1.0\sigma$, result in best-redshift solutions that are similar to the analysis using the median (zero offset) scaling relations. This suggests a good robustness of our joint-analysis method against the existing scatter of $L_{FIR}-L_{line}$ scaling relations in observations. As for the test with $+1.0\sigma$ offset, we further determined the reason for the discrepancy in the best-redshift solution. For HLS-22 and HLS-2-2, the application of the $+1.0\sigma$ offset leads to mismatches between the model and the data due to glitches or noise spikes in the spectra (see Fig. A.2 as an example, where the noise spike matches CO(5-4) at the best-redshift solution of HLS-22 here). For HLS-3, we find this unlikely low redshift solution after applying the $+1.0\sigma$ offset as the code assigns the only strong emission line in the spectral scan at 139.746 GHz to be CO(2-1). This suggests the stronger demand that the spectral scan be wide enough to cover more than one strong spectral line in the redshift confirmation of DSFGs with moderate redshift (i.e., $z \sim 2-3$).

To test the self-consistency of the joint-analysis method, we placed the measured CO line fluxes and far-IR luminosities of HLS sources at their best-redshift solutions (the $z_{best,med}$ in Table D.1) in the corresponding $L_{FIR}-L_{CO}$ diagrams to determine whether they follow the scaling relations used for line flux predictions. The results are shown in Fig. D.1. The plotted CO and far-IR luminosities were derived using the Gaussian-fit line fluxes and modified blackbody fitting (see Table. 9 and Table. 10). In addition to the average $L_{FIR}-L_{CO}$ correlations, we also show the samples from Cañamas et al. (2018) with multiple line transitions from the same sources as a comparison to our sources.

From Fig. D.1, we find that most of our sources fall within $+2\sigma$ of the scaling relation used in our analysis, which is also consistent with the regions occupied by bright submm galaxies in Cañamas et al. (2018). The only exception is HLS-3, which is also highlighted in Fig. D.1. In either the $L_{FIR}-L'_{CO(5-4)}$ or the $L_{FIR}-L'_{CO(4-3)}$ diagram, HLS-3 falls well below the scaling relation even when we consider the scatter of these scaling relations. However, we also note that its SPIRE photometry is one of the poorest in all of the four HLS sources. As the SPIRE bands close to the peak wavelength of SED predominantly constrain the IR luminosity, it is likely that the IR luminosity of HLS-3 is much less constrained than the other HLS sources, especially compared to HLS-2 and HLS-4.
<figure><img src="image_13.png" /><figcaption>**Fig. D.1.** Comparison of the $L_{FIR}-L_{CO}$ correlations of different transitions (Greve et al. 2014; Liu et al. 2015) and our sample based on measurements from our observations at the best-redshift solutions. The sources with upper limits on the line luminosities are presented as leftward triangles.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. D.1.** Comparison of the $L_{FIR}-L_{CO}$ correlations of different transitions (Greve et al. 2014; Liu et al. 2015) and our sample based on measurements from our observations at the best-redshift solutions. The sources with upper limits on the line luminosities are presented as leftward triangles.</figcaption></figure>