Probing the interstellar medium of the quasar BRI 0952 − 0115 An analysis of [C ii ], [C ii ], CO, OH, and H 2 O

Context. The extent of the e ﬀ ect of active galactic nuclei (AGN) on their host galaxies at high-redshift is not apparent. The processes governing the co-eval evolution of the stellar mass and the mass of the central supermassive black hole, along with the e ﬀ ects of the supermassive black hole on the host galaxy, remain unclear. Studying this e ﬀ ect in the distant universe is a di ﬃ cult process as the mechanisms of tracing AGN activity can often be inaccurately associated with intense star formation and vice versa. Aims. Our aim is to better understand the processes governing the interstellar medium (ISM) of the quasar BRI0952 − 0952 at z = 4 . 432, speciﬁcally with regard to the individual heating processes at work and to place the quasar in an evolutionary context. Methods. We analyzed ALMA archival bands 3, 4, and 6 data and combined the results with high-resolution band-7 ALMA observations of the quasar. We detected [C i ](2–1), [C ii ]( 2 P 3 / 2 − 2 P 1 / 2 ), CO(5–4), CO(7–6), CO(12–11), OH 2 Π 1 / 2 (3 / 2 − 1 / 2), and H 2 O(2 11 − 2 02 ), and we report a tentative detection of OH + . We updated the lensing model and we used the radiative transfer code MOLPOP-CEP to produce line emission models which we compared with our observations. Results. We used the [C i ] line emission to estimate the total molecular gas mass in the quasar. We present results from the radiative transfer code MOLPOP-CEP constraining the properties of the CO emission and suggest di ﬀ erent possible scenarios for heating mechanisms within the quasar. We extended our results from MOLPOP-CEP to the additional line species detected in the quasar to place stronger constraints on the ISM properties. Conclusions. Modeling from the CO SLED suggests that there are extreme heating mechanisms operating within the quasar in the form of star formation or AGN activity; however, with the current data, it remains unclear which of the two is the preferred mechanism as both models reasonably reproduce the observed CO line ﬂuxes. The updated lensing model suggests a velocity gradient across the [C ii ] line, suggestive of ongoing kinematical processes within the quasar. We ﬁnd that the H 2 O emission in BRI0952 is likely correlated with star-forming regions of the ISM. We used the molecular gas mass from [C i ] to calculate a depletion time for the quasar. We conclude that BRI0952 − 0952 is a quasar with a signiﬁcant AGN contribution while also showing signs of extreme starburst activity, indicating that the quasar could be in a transitional phase between a starburst-dominated stage and an AGN-dominated stage.


Introduction
The formation and evolution of massive galaxies in the highredshift universe remains to be understood in detail.The correlation between the mass of the central supermassive black hole (SMBHs) and the velocity dispersion of the galaxy suggests that there is a co-evolution of these two properties (e.g., Kormendy & Ho 2013).This correlation has been confirmed by several studies (e.g., Magorrian et al. 1998;Ferrarese & Merritt 2000;Gebhardt et al. 2000;Häring & Rix 2004;Gültekin et al. 2009;Beifiori et al. 2012;Kormendy & Ho 2013;Bennert et al. 2015;Reines & Volonteri 2015).A common method of studying the interplay between the central active galactic nuclei and the ISM is by studying the host galaxies of active galactic nuclei (AGN) across cosmic time and, in particular, studying quasars and their host galaxies.High-redshift quasar host galaxies are often found to be dust-rich and star forming, possibly undergoing a starburst phase (e.g., Venemans et al. 2018Venemans et al. , 2020)).Thus, detailed analyses of the interstellar medium (ISM) properties of AGN host galaxies and other dust-rich starburst galaxies in the high-redshift uni-verse have become increasingly common to understand the early phases of the evolution of massive galaxies.Several studies have been conducted in this area, such as those by Chen et al. (2017), Circosta et al. (2019), Venemans et al. (2019), Li et al. (2020), Jarugula et al. (2021), andHagimoto et al. (2023).
The ISM consists of multiple gas phases, including molecular, neutral, and ionized gas phases.The heating and cooling mechanisms of these phases, as well as the gas and dust content, imply that studies considering a range of properties are required in order to establish a complete understanding of the physical conditions in a galaxy.At high redshift, the CO molecule is often employed as a tool for investigating the ISM due to the variation in density and temperature necessary to excite different population levels of CO (e.g., Carilli & Walter 2013).The low-J (J < 3) CO lines trace regions of the ISM with lower densities and temperatures, and are commonly used to trace the bulk molecular gas mass (e.g., Ivison et al. 2011).Mid-J (4 ≤ J ≤ 8) CO lines trace warmer and denser regions of the ISM.These lines have also been shown to linearly correlate with infrared luminosity Notes. (a) ALMA band for each line observation. (b) Date of ALMA observation. (c) Number of antennas used in observation. (d) Central frequency of the spectral window containing the specific line. (e) Channel width of the calibrated and imaged data for the specific line. (f) Synthesized beam of the spectral window. (g) Per-channel RMS of the spectral window.
resolution vary depending on the specific observation set.Additionally, all archival observations were performed in continuum mode, including those containing lines, which affects the resolution and noise properties between channels, but is negligible for the purposes of this paper.
Similar to the band 7 data, all calibration and imaging processes for the archival band 3, 4, and 6 data were performed in the Common Astronomy Software Application package (CASA McMullin et al. 2007).The data in each band were calibrated using the CASA pipeline calibration scripts in version CASA 4.5.1 for the band 4 and 6 data and CASA version 5.1.1 for the band 6 data.The pipeline includes calibration of the phase, bandpass, flux, and gain.The following phase calibrators were used; band 3: J1058+0133, J0948+0022, J1000+0005; band 4: J1058+0133, Ganymede, J0948+0022; band 6: J0854+2006, J1037-2934, J0948+0022.The pipeline calibrated data were checked to ensure quality.Imaging was performed in CASA using the TCLEAN task, cleaning in the region of expected emission using a circular mask with a radius of ∼ 20 ′′ from the band 7 data including both BRI 0952 and the companions detected in Kade et al. (2023).The images were cleaned down to a 2σ level for each band.Pipeline-provided images of the calibrated data from the ALMA data reduction pipeline were visually inspected to ensure data quality and identify line-free channels.Continuum subtraction was performed using CASA's task UVCONTSUB using line-free channels and a polynomial fit order of 0 for all bands.Continuum images were created using the line-free channels with a Briggs weighting scheme with a robust parameter of 0.5.Spectral cubes were created employing the same weighting scheme using the continuum subtracted measurement set, with the exception of the band 3 data which was imaged using natural weighting.As mentioned above, the band 3 data were of significantly higher angular resolution than the other bands.In order to improve data continuity, a taper was performed to the uv-data which effectively lowered the resolution of produced continuum images and line-emission cubes to approximately match that of the band 7 data.The band 3 data were also of significantly higher spectral resolution than the other bands (∼ 5.5 km s −1 ); we binned the data by a factor of 8× in order to improve the signal-to-noise ratio in each channel.To account for the quasar not being centered in the observational field, primary beam correction was applied to the band 3 data in addition to the taper.However, primary beam correction was not necessary for the band 4 and 6 data since the observations were centered on the source.A conservative estimate of 10% is used to represent the uncertainty in the absolute flux calibration (note that the fiducial value for bands 3 and 4 is ∼ 5% and for bands 6 and 7 is ∼ 10%) 1 .
Table 1 provides details on the different bands and spectral windows in which lines were detected including synthesized beam sizes, observing time, spectral resolution, and final image rms'.

Continuum
We have detected continuum emission in bands 3, 4, 6, and 7.For details on how we extracted the band 7 continuum data, refer to Kade et al. (2023).For the archival data, we extracted continuum fluxes using CASA's IMFIT routine following Kade et al. (2023).Continuum fluxes are provided in Table 2 and images are shown in Fig. 1.As the spatial resolution of the band 4 and 6 data is lower than that of the band 3 and 7 data, the two images of the quasar are not resolved in the continuum images presented in Fig. 1.While extended emission is visible in the band 4 continuum data, it is improbable that it arises from the south-western companion as the continuum sensitivity is poorer in the band 4 data than in the band 7 data, and the south-western companion was not detected in the latter (see Kade et al. 2023).Without higher resolution continuum data it is not possible to draw robust conclusions about this emission feature.  a) Continuum in band 3 is only detected in the northern image of the quasar, therefore we report a 3σ upper limit for the southern image of the quasar.
Although the angular resolution of both the band 3 and band 7 data are sufficient to distinguish between the two images (Img-S and Img-N; shown in Fig. 1) of BRI 0952, we do not detect Img-S in the band 3 data.This non-detection can be attributed to two factors.Firstly, the relatively faint continuum of Img-S is compounded with the high angular resolution of the obser-vations, making it difficult to detect.Secondly, the frequencies covered by the band 3 data, which probe the long-wavelength end of the Rayleigh-Jeans tail, are intrinsically fainter than those covered by the band 7 data, which probe closer to the peak of the Rayleigh-Jeans regime.This difference in frequency coverage, combined with the native angular resolution of the band 3 observations, may have resulted in the fainter flux from Img-S to go undetected.We note that Img-S is also not visible in the non-tapered continuum image, shown in
We do not detect the companion sources, Comp-N or Comp-SW, reported in [C ii] emission in Kade et al. (2023) in any of the additional line species reported in this paper.This is not surprising as the companion sources are very faint in comparison to the quasar (< 10% the flux of the quasar in [C ii] ), and the band 7 [C ii] data are the data of highest quality among the observations used in this paper.Although the band 3 data has higher native spectral and angular resolution than the band 7 data, it has a shorter integration time and lower sensitivity.This, combined with the significantly lower observed brightness of CO(5-4) in relation to [C ii] limits the possibility of detecting emission from Comp-N and Comp-SW.Higher angular resolution data with longer integration times than that of the band 4 and 6 data would be required to detect these sources.In general, the observations reveal compact unresolved emission.
We perform a regional spectral extraction for each species detected, meaning we extract the spectrum for a given species from a specific region, shown by the black ellipses in Fig. 2, at all σ levels.We opt to extract the spectra from these regions rather than performing a σ cut as this method has the benefit of preserving emission at σ levels lower than the chosen cut level; this is particularly important towards the edges of emission regions where the lower-brightness surface emission is distributed2 .This is crucial for fainter the lines, specifically [C i], H 2 O, and OH + lines.We choose the region of extraction to cover the extent of the emission and, where possible, exclude the location of the companions detected in Kade et al. (2023).In some cases it was not possible to entirely exclude the spatial regions in which the companions manifest in [C ii], however, given the relative faintness of the companions in [C ii] compared to BRI 0952, if there was any contribution from the companions to the flux of the line, it would likely be negligible.Although it would be beneficial to use the same region of extraction for all species, specifically the region used in Kade et al. (2023), the differing resolutions results We calculate the rms of the spectra individually for each species as follows.We obtained the per-channel rms of the spectral cube by sampling twelve distinct emission-free regions chosen to be the same size as those utilized for extracting the spectra of the specific species.Fig. 3, and B.1 present the spectra of each species detected in BRI 0952, centered on the redshifted systemic velocity of the line with the exception of the [C i] line which is offset from the CO(7-6) line by its relative velocity difference.The moment-0 and moment-1 maps were created from the same regions as used to extract the spectra, however, here we remove flux below a 2σ level in the moment-1 maps for clarity, shown in Fig. 2 along with the velocity range used for extraction.Table 3 provides details on the detected emission lines.
We calculate integrated flux densities and subsequently line luminosities using the following equations from Solomon et al. (1997): (2)

CO emission
The CO(5-4) observations were carried out in an extended configuration with a maximum recoverable scale of ∼ 23 kpc.We detect CO(5-4) emission at 1.50±0.32Jy km s −1 , which is within 1.5σ of the previously reported value of 0.91 ± 0.11 Jy km s −1 from Guilloteau et al. (1999).This suggests that any filtered out diffuse emission is negligible3 .Therefore, it is unlikely that we are missing significant flux.
The spectrum of the CO(7-6) emission exhibits a 'bump' in the blue part of its spectrum at velocities < −100 km s −1 .We fit this line with two Gaussian profiles to account for this.Since the current data is limited by the angular resolution, a spatial correlation between the bump and Img-N or Img-S is not feasible.We must therefore explore other methods of investigation.One possible explanation for the bump is that it is an effect of the gravitational lensing.One simple method of investigation is to extract the spectra from only the very central region of Img-N of the quasar where the flux of the CO line is highest.Hence, we extract the spectra from Img-N at levels >10σ; however, in this analysis the bump remains.A second investigation would be to use the magnification factor corrected spectra; the process for this is described in Section 3.3; similarly, in this investigation, the bump remains.It is possible that this feature originates from the gravitational lensing given that the lensing corrected spectra is still biased by the extrapolation of the [C ii] magnification factor per channel to the CO(7-6) spectra (Section 3.3).Higher angular resolution data would be necessary to further investigate this possibility.An additional explanation is that this emission originates from one of the companion galaxies detected in Kade et al. (2023), however we discard this possibility as the [C ii] emission from Comp-N and Comp-SW is located at very similar systemic velocities as that of BRI 0952.We suggest that this bump provides an indication of the gas kinematics operating in BRI 0952.The slight offset or asymmetry in the line profile of the CO(7-6) line is reminiscent of the line profile of CO(6-5) detected in the strongly lensed starburst galaxy G09v1.97 (Yang et al. 2019) in which the authors argued that the complex line profile was an indication of two interacting galaxies.It remains to be seen if this is the precise case for BRI 0952 but it is suggestive of an on-going kinematical process, be it rotation or interaction.
There appears to be an elongation of the CO(12-11) emission with respect to the semi-major axis of the emission compared to other emission lines, however the large and extended beam with a similar position angle to the extension of the  data suggests that this is merely an affect of the beam.
We bin the H 2 O spectra by a factor two to improve the per channel signal-to-noise ratio.We note that the H 2 O line follows the trend observed by Omont et al. (2013); Yang et al. (2017); Stanley et al. (2021); Li et al. (2020) of having a similar line profile to that of the CO(7-6) line with a full width half max (FWHM) of 257 ± 74 km s −1 to the 259 ± 20 km s −1 of the CO(7-6) line.The CO(5-4) line is significantly broader than the CO(7-6) line, indicating the possibility that the CO(5-4) line originates from a different region than that of the H 2 O and CO(7-6) emission.We show the profile of the normalized H 2 O , CO(5-4), CO(7-6), and [C i] spectra in Fig. 5.We include the [C i] emission as it is likely that these lines trace similar regions of the ISM, see Sections 3.5 and 5.1.We note that the line profile of the H 2 O emission does not exhibit the so-called 'bumps' that we find in the CO spectra.One possible explanation for this is the relative faintness of the H 2 O line or the possibility that the 'bumps' arise from, for example, an interaction with a companion source.Alternatively, if the H 2 O originates from a different region, specifically a smaller and more compact region, of the galaxy than the CO(5-4) and CO(7-6) lines perhaps this feature would simply not be present in the spectra regardless of data quality.Indeed, Quinatoa et al. (in prep) find different H 2 O and CO line profiles and postulate that the difference is due to a discrepancy in spatial distribution of the emitting regions.

OH
We detect OH emission towards BRI 0952.As this is a doublet we fit the spectra in a slightly different manner than described for the other lines.We fit an individual Gaussian for both lines of the doublet, but we fix the peak velocity to be the systemic velocity of each line in the doublet.We do not fix the FWHM of the Gaussian to be the same value, but both have the same boundary conditions.We do this in order to determine if there are outflow signatures in the spectra in the form of red-or blueshifted emission from the systemic velocity.We do not find that the Gaussian fits are offset from the observed emission, see the bottom left image in Fig. B.1, providing no indication of in-or out-flowing gas.We find that the individual lines composing the doublet are nearly symmetrical with very similar FWHMs and line fluxes, see Table 3, in agreement with the detected CO lines.

OH +
We report a tentative detection of OH + emission towards BRI 0952.This emission arises solely from Img-N in the quasar, see Fig. 2. We extract the spectra from a region smaller than that used for the [C ii] and OH emission to limit the effect of the noise and we bin the spectra by a factor of two to improve the per channel signal-to-noise ratio.The spectrum of the OH + is shown in
Here we expand upon the lensing model described in Kade et al. (2023).While other papers, such as Apostolovski et al. (2019), have detected multiple atomic or molecular species and modelled their emission individually using gravitational lensing modelling codes like Visilens (Spilker et al. 2016), our paper focuses solely on modelling the [C ii] emission.This approach is necessary due to the varying data quality obtained from a combination of archival data and targeted high-resolution observations.The [C ii] line emission provides the highest quality observation of line emission in BRI 0952 in terms of signal-to-noise, sensitivity, and angular resolution.Moreover, this emission is thought to originate from both colder neutral regions of the ISM and photo-dissociation regions (PDRs), thus likely tracing a large extent of the galaxy (e.g., Carilli & Walter 2013;Pavesi et al. 2018).As a result, we use the [C ii] emission to determine the magnification factor across the [C ii] line, which is then extrapolated to the other atomic and molecular species detected in BRI 0952 based on their systemic velocities.This approach is preferred because it enables us to leverage the high-quality [C ii] observations to derive magnification factors for other species without having to model each species individually using lowerquality data.
It is important to note that our approach of extrapolating magnification factors from the [C ii] emission to other species does not take into account that different emission lines likely originate from different locations within the galaxy.For instance, high-excitation CO lines may not necessarily originate in the same locations where [C ii] emission is present.Nevertheless, as long as the other line species detected follow the bulk motion and systematic kinematics of the entire galaxy, as traced by [C ii], this method of extrapolating the magnification factor is valid.Biases in the magnification factor between species would arise in cases where individual species exhibit peculiar motions such as outflows or shocks which break from the bulk motion of the gas.While important to consider, with the current data available, specifically with regard to resolution limitations, the errors associated with performing a lensing analysis on lines with insufficient resolution to resolve the two images of the quasar would lead to an unnecessary increase in the systematic errors associated with those methods.

Visilens Modelling
We utilize the python-based code Visilens (Spilker et al. 2016) to model each channel containing [C ii] emission.Visilens is designed specifically for modeling radio observations of gravitationally lensed systems.The lens parameters used in our study are from Kade et al. (2023), where a full description of these parameters is provided.Although Visilens does not generate a source plane reconstruction in the form of a spectral cube, we can perform a spectral reconstruction by modeling each channel in which [C ii] emission appears and adjusting for the bestfit magnification factor for that specific channel.Fig. 3 and B.1 displays the magnification factor corrected spectra compared to non-corrected spectra.
We investigate four different methods to model the [C ii] emission, wherein we describe the source as having a Sérsic profile parameterized by its position relative to the lens (x S , y S ), flux density, Sérsic index, half-light radius, axis-ratio, and position angle.Rather than optimizing each channel per model separately, we fit every channel using the same set of parameters.Although it is possible that individual channel optimization would marginally improve the accuracy of the model, we conclude that this would unnecessarily introduce user error without sufficient physical motivation.Below we describe the different models: 1. Model 1 (model 1 in Fig. 6 and C.1): We let the parameters vary within a reasonable range around best-fit values from Kade et al. (2023).2. Model 2(model 2 in Fig. 6 and C.1): The position of the source and Sérsic index were fixed.We fixed the Sérsic index to n = 2.5 (the average value of n in model 1), and the source position to the best-fit source position in each channel from the first model.The remaining parameters are allowed to vary in the same range as in the first model.3. Model 3 (model 3 in Fig. 6 and C.1): The position of the source and the Sérsic index were fixed.The Sérsic index was fixed to n = 2.5 as in model 2. In order to determine the source position in each channel we fit a linear regression of the best-fit source position from the first model and adopted those values as the positions in this model.The remaining parameters were allowed to vary in the same range as in the first model.4. Model 4 (model 4 in Fig. 6 and C.1): The parameters were set to be the same as in model 1, however we bin the input data by a factor of 2 in order to increase data quality.This is particularly important in the high-and low-velocity wings of the [C ii] line not only due to the relatively lower signal to noise in these channels but it is also of special import due to the wing-like structures noted in Kade et al. (2023) which may be indicative of differential magnification and could affect the lensing correction.We note that we attempted to fit the [C ii] emission using a Gaussian source profile.This model is parameterized by its position relative to the lens (x S , y S ), flux density, and Gaussian width.Although this model is more simplistic compared to the Sérsic source parameterization, the physical motivation for one model over the other is not sufficiently clear so as to exclude a Gaussian source being a good description of the light profile of BRI 0952.However, the fit from VISILENS for this model was sufficiently poor that we disregarded it.
We show the magnification factor as a function of velocity for the different methods described above in Fig. 6.We show the change in individual parameters (source position, flux, position angle, etc.) in Fig. C.1.We note that there is no large discrepancy for in the output best-fit parameters between the different models.Most notably, the magnification factor across the velocity range does not change greatly depending on the model used.This is significant for models 1 and 4, specifically in the wings of the [C ii] line where the emission has lower signal-to-noise ratios.Despite the low emission levels, these models remain consistent with each other, indicating that the variation in magnification factor with velocity is not solely influenced by the lower flux in the wings.Due to the small discrepancies between different methodologies, we choose the most simple and reproducible to use throughout the remainder of this paper; namely, model 1 and model 4. We use the magnification factors obtained from the first model to reconstruct the spectra.An important conclusion from Fig. 6 is that the magnification factor is not symmetric or stable across the [C ii] line; therefore, conclusions drawn from the shape of line profiles should only be done using a spectra which has been corrected for magnification.The methodology used to correct the line profiles is described in Section 3.3.2.
We investigate the change in the position and half-light radius of the best fit Sérsic source model across the [C ii] line.We plot the best-fit position and half-light radius as a function of velocity in Fig. 7 using model 4 to improve data quality through binning and ensure clarity.We find clear evidence of a velocity gradient across the line profile, suggesting that the galaxy may be rotating.This result is reminiscent of the results from Apostolovski et al. ( 2019) and Yang et al. (2019).

Spectral Reconstruction
We extrapolated the lensing factor based on the velocity relative to the center of the line.To correct the spectrum of the other line species, we used a closest value finding algorithm to search for the velocity bin in the [C ii] data which is nearest the velocity bin of the particular line and that specific channel was then corrected for the magnification factor in the closest [C ii] bin.We overplot the magnification corrected spectra in Fig. 3, where the corrected spectra are shown in yellow and the non-corrected spectra are shown in gray.
We limit this correction to the lines in which the line profile provides relevant information and high signal-to-noise ratios; that is, the [C ii], CO, and OH lines.We note a slight offset of the peak of both OH lines in the magnification corrected spectra.This could be a possible indication of gas motion, but given the errors associated with this approach and the limitations of the data, we do not investigate further.The reconstructed CO(7-6) line retains its 'bump' in the blue part of the spectrum and it remains unclear what the source of this asymmetry could be.The CO(12-11) line seems to exhibit a double-horned shape, which is could be a feature indicative of rotation or absorption, however, given the low angular-resolution and the errors associated with this type of lensing magnification correction, we do not investigate this feature further.Additionally, the broad [C ii] wings reported in Kade et al. (2023) remain, confirming that these are not due to differential lensing as hypothesized in Kade et al. (2023).Although the true origin of the wings remains unclear, this confirms that they do not arise from the lensing.We present the line luminosity values using µ = 3.92 (Kade et al. 2023) and the per-channel µ correction in Table 3.It is worth noting that the magnification-corrected luminosities generally align within the margin of error comparing between the two methods.For the subsequent sections of this paper we discuss the fluxes and luminosities of the [C ii] and CO lines employing the magnificationcorrected value obtained through the per-channel µ methodology.
This approach suffers from the limitations of the method, as described in 3.3.1.Given the constraints imposed by the differing qualities of the ALMA archival data we continue with this method for the remainder of the analysis.For a full analysis of the magnification factor of each individual line to be performed, a consistent survey would need to be done of the CO and [C ii] lines with the same sensitivity, angular, and spectral resolutions.

Gas Mass
A common procedure to estimate a galaxy's gas mass is to use CO(1-0) and assume an α CO to convert from CO to H 2 whereby a total gas mass can be estimated.In cases where CO(1-0) is not observed, a conversion factor is used to translate to CO(1-0) luminosity which is subsequently converted to a total gas mass using the above procedure.However, there are a number of important caveats to these approaches.First, α CO has been shown to vary to a significant amount due to its dependencies on gas density and temperature as well as with galaxy type (e.g., Bolatto et al. 2013, and references therein).In addition, conversions between any observed above-ground state CO line down to the ground state are also dependant on the conversion factor (e.g., Boogaard et al. 2020).A substitute tracer for gas mass is the ground state [C i] (1-0) line at 492 GHz.Studies have shown that this line is linearly correlated with CO(1-0) (e.g., Carilli & Walter 2013, and references therein).Hence, we can estimate the total gas mass using the equation from Weiß et al. (2003); where Q(T ex ) is the partition function defining the level difference between populations in the upper and ground state of . This is given by the following: where T 1 = 23.6K and T 2 = 62.5 K corresponding to the [C i] (1-0) and [C i] (2-1) transitions respectively.For the purposes of this paper we adopt a excitation temperature of T ex = 30 K following Walter et al. (2011) and Jarugula et al. (2021).Using this method we find a [C i] gas mass of (2.1±1.7)×10 6M ⊙ .
Converting from the above calculation to an H 2 gas mass can be done using a [C i] /H 2 ratio corresponding to an assumed abundance of [C i] .The chosen abundance ratio value can affect the total molecular gas mass by a factor ∼ 10 and may vary depending on the galaxy type or region from which the [C i] emission arises.For this reason, we choose to use [C i] /H 2 = (8.4± 3.5) × 10 −5 from Walter et al. (2011) which is calibrated from a sample of quasars and submillimeter galaxies at z > 2. However, it should be noted that this value may change in high extinction regions (Walter et al. 2011, and references therein).
From this result we obtain an H 2 gas mass of (6.2 ± 4.9) × 10 9 M ⊙ .This result is similar to the H 2 gas mass reported by Guilloteau et al. (1999) where they find M H 2 ∼ (2 − 3) × 10 9 M ⊙ .However, it is significantly lower than the H 2 gas mass reported by Hughes et al. (1997) where they find a value of 1.35×10 11 M ⊙ .However, the methodology varies widely with the different values.Guilloteau et al. (1999) use a conversion from the CO(5-4) to gas mass where Hughes et al. (1997) calculate the dust mass through continuum measurements (which are likely contaminated by the AGN) and subsequently convert the obtained dust mass to an H 2 mass.We note that the gas mass obtained through the [C i] (2-1) emission is similar to the dynamical mass from the [C ii] emission reported in Kade et al. (2023).Given the previously mentioned caveats, this value should be interpreted as a proxy and not a robust measurement.

Atomic Carbon vs. CO(7-6)
The ratio of mid-J CO lines to [C i] can also be conceptualized as a means of distinguishing the amount of bulk molecular mass (traced by [C i] in lower density environments) versus that of molecular mass in dense environments (traced by mid-J CO lines).This was suggested in Papadopoulos & Geach (2012) for CO(4-3) and [C i] (1-0), however studies have extended this method to CO(7-6) and [C i] (2-1) (e.g., Andreani et al. 2018).
The CO(7-6) line is commonly taken to be a reliable tracer of star-formation due to its specific excitation conditions (e.g., Greve et al. 2014;Lu et al. 2015;Yang et al. 2017); with excitation potential T ex = 154 K and n = 4.5 × 10 5 cm −3 (Carilli & Walter 2013) this line originates from significantly denser and warmer regions of the ISM than the [C i] (2-1) line which has an excitation potential of T ex = 63 K and n = 1.2×10 3 cm −3 (Carilli & Walter 2013).
An additional powerful tracer of colder and more diffuse regions of the ISM are the [C i] lines.Studies have shown that [C i] line emission excitation temperature and carbon abundance with respect to hydrogen in low-and high-redshift sources does seem to exhibit a strong evolution with redshift, nor does there appear to be a strong discrepancy in [C i] emission between galaxy types (i.e., SMG or AGN; Walter et al. 2011).[C i] emission has the additional benefit of being less affected by radiation from the CMB and cosmic rays (e.g., Zhang et al. 2016;Bisbas et al. 2015), therefore the [C i] lines have become a common proxy for calculations of the bulk molecular gas mass in galaxies at both low and high redshift (e.g., ) and the ratio of [C i] to mid-J CO lines has been used as a tracer for the amount of dense versus diffuse gas (e.g., Andreani et al. 2018).
The excitation conditions of [C i] (2-1) and CO(7-6) are sufficiently different so as to make the ratio between the two a meaningful comparison regarding their excitation properties.Indeed, Andreani et al. (2018) find this ratio to be approximately unity for disc-like galaxies such as the Milky Way whereas it may be up to a factor 10 higher in for example merger-driven galaxies.We find a ratio of CO(7-6)/[C i] (2-1) ∼ 4.8.This is significantly higher than the values found in Andreani et al.Notes. (a) Corrected for lensing using µ = 3.92 from Kade et al. (2023). (b) Corrected for lensing using the per-channel magnification factor from the [C ii] emission as described in Section 3.3. (c) As the OH line is a doublet, we fit it with two single Gaussian components with a fixed distance between peaks, see Section 3.2.3 for further discussion.
(2018), suggesting that a merger driven scenario could be occurring.Scholtz et al. (2023) reported detections of [C i] (2-1) and CO(7-6) in a sample of z ∼ 2.3 extremely red quasars, from their reported values the average of their is ∼ 2. The value we find for BRI 0952 is significantly higher than the ratio found for the dusty star-forming galaxy (DSFG) SPT0311-58 (z = 6.9;Jarugula et al. 2021) and the quasar J2348-3054 (z = 6.9;Venemans et al. 2017).

Radiative Transfer Modeling
We Given these parameters, MOLPOP-CEP solves the level population problem between each zone in the slab and outputs an emergent flux which can be directly compared to observations.Asensio Ramos & Elitzur (2018) suggested that even when the parameters are held constant across the entire slab, dividing the slabs into a number of zones helped achieve more accurate results.Hence, we divide the slab into 10 zones and employ a uniform setup in each zone throughout the slab following Asensio Ramos & Elitzur (2018) and Li et al. (2020).This approach is particularly important for optically thick lines since the strength of the line emission depends on the distance into the cloud and therefore the the level populations of these lines is dependent on their position within the slab geometry.
We generate a grid of models by varying the first three parameters, the zone width (∆L), the zone gas volume density n(H 2 ), and the zone kinetic temperature (T kin ); these parameters provide the most direct means of identifying different regions of the ISM.The zone width was allowed to vary from 16 cm to 19 cm in steps of 0.25 cm.The zone gas volume density was allowed to vary between 3 cm −3 to 6 cm −3 in steps of 0.25 cm −3 .The zone kinetic temperature was allowed to vary between 50 K to 200 K in steps of 25 K.The molecular abundance (X species = [species]/[H 2 ]) was set to 10 −4 for CO and 10 −7 for H 2 O, which are typical in molecular gas region (e.g., Weiß et al. 2003;Liu et al. 2017).The final grid consists of a total of 13×13×9 = 1521 models.
It is important to note the interdependence of the zone width (∆L) and zone volume density n(H 2 ) and the column density within the entire slab (N mol ).The total column density, defined as the sum of the column densities within each slab, is proportional to ∆L through the following: where the subscript denotes the specific atom or molecule under consideration and the factor of 10 represents the number of zones.As N species is in terms of both the zone column density and the zone width, there is a degeneracy between these two parameters which is important to consider when comparing different models and therefore throughout the remainder of the paper we use N species when presenting our results.
We also take into account the effect of the CMB at the source redshift (z = 4.432).At this redshift the CMB temperature is T CMB = 14.742K using the relation Hence, the CMB is sufficiently hot to act as a heating source for the CO emission at the redshift of the quasar (e.g., Asensio Ramos & Elitzur 2018; Li et al. 2020).This is particularly important for line transitions with lower excitation temperature as those lines are more heavily affected by this phenomenon (e.g., Obreschkow et al. 2009;da Cunha et al. 2013a).
Apart from allowing for variation of the molecular gas physical parameters within the slab, MOLPOP-CEP also offers the flexibility to incorporate various radiation fields.These radiation fields include the CMB (characterized by the CMB temperature at the specific redshift), blackbody radiation, dust heating, (characterized by the dust temperature and optical depth), and the inclusion of any other radiation field described by its SED.
In regard to the CO SLED we consider heating from PDRs and XDRs while for the H 2 O emission we consider heating via dust.

Heating via PDRs
We investigate the contribution of the stellar heating via PDRs to the CO excitation in the quasar.Should the CO-emitting cloud be situated in a star-forming region, then the radiation from nearby young massive stars (mainly type O and B stars) should contribute to the excitation, namely to the mid-J lines.To account for the contribution from a stellar radiation, we used the stellar population synthesis code Starburst99 (Leitherer et al. 1999) to generate an SED using the same assumptions as Vallini et al. (2019) and assuming an SFR of ∼ 3000 M ⊙ yr −1 .The stellar radiation field strength is generally given in terms of the interstellar FUV flux, G 0 (e.g., Vallini et al. 2019): where G 0,MW = 1.6 × 10 −3 erg cm −2 s −1 (Habing 1968).We can express the strength of the radiation field from the SED in terms of G 0 and scale this value accordingly to increase the radiation field strength.For example, assuming an SFR of 3000 M ⊙ yr −1 , and setting SFR MW = 1 M ⊙ yr −1 (following Vallini et al. 2019), the scale of the bolometric stellar radiation at the molecular source is set to be 4.8 × 10 −3 W m −2 .In our analysis we consider the stellar radiation SED scaled as log(G 0 ) = 0.0, 1.0, 2.0, 3.0, and 4.0.

Heating via X-rays from the AGN
An alternative heating scenario for the gas in BRI 0952 is heating via the AGN.We model the AGN radiation in the form of XDRs using an SED for an extreme super-Eddington Type-1 quasar as described in Jin et al. (2017) from CLOUDY covering the frequency range 2.99 × 10 16 -5.98× 10 19 Hz4 .MOLPOP-CEP requires that the flux level of the SED (corresponding to the XDR flux) at the emitting region be specified in terms of both the overall luminosity and the distance to the emitting region.Here we assume a radiation source of luminosity 10 13 L ⊙ at distances of 500 pc, 1.0 kpc, and 1.5 kpc to the emitting cloud, normalized to bolometric energy density 1.45 × 10 −1 W m −25 .

One-Component Models of the ISM
We model the CO SLED of BRI 0952 for each of these heating mechanisms, namely PDRs and XDRs, individually.To evaluate the goodness of the fit for each of the models we use a χ 2 test, where the χ 2 value is a weighted sum of squared deviations.In order to compare the accuracy of the different heating mechanisms we consider only those models with χ 2 ≤ 1σ where 1σ = 2.30 for the degrees of freedom in our model.We show the best fit and 1σ models in which G 0 is scaled to represent different stellar radiation field strengths and for AGN heating at 1.0 kpc and 1.5 kpc in D.1.We discard models with AGN radiation modeled at a distance of 0.5 kpc as we do not find a reasonable fit to the data in this case.Both heating scenarios fail to reproduce the observed CO fluxes; generally XDR-heating models are able to reproduce the CO(7-6) and CO(12-11) but fail to reproduce the CO(5-4) emission while PDR models reproduce only the CO(12-11) emission and struggle to reproduce the lower-J CO line emission.In the stellar heating regime, the best-fit model corresponds to log(G 0 ) = 3.0 and T kin = 75.0K and log(N CO /cm −2 ) = 18.0.We note that there is no difference in the χ 2 value between log(G 0 ) = 0.0, 1.0, and 2.0.In the XDR regime, we find a preference for the XDR to be situated at a distance of 1.5 kpc from the emitting region where the best-fit model has the same physical conditions as the log(G 0 ) = 3.0 model.Given the inability of these models to reproduce the observed CO SLED, we investigate additional heating scenarios below.

Two-Component Models of the ISM
It is perhaps not surprising that models that employ a single heating mechanism are insufficient to reproduce the observations of the low-, mid-, and high-J CO transitions due to the different excitation requirements of these transitions.Studies of both local and high-redshift galaxies with high-J detections similar to the CO(12-11) line in BRI 0952 have been unable to successfully reproduce their observations without the use of a two component model (e.g., Weiß et al. 2007;van der Werf et al. 2010;Gallerani et al. 2014;Yang et al. 2019;Li et al. 2020).More specifically, in the case of the local quasar Mrk231 (van der Werf et al. 2010), the high-redshift quasars APM08279+5255 (z = 3.9) (Weiß et al. 2007) and J1148+5251 at (z = 6.4) (Gallerani et al. 2014) an AGN component contributing XDRs was a required to successfully model the CO SLED.However, as noted in Weiß et al. (2007), the hot component need not be directly correlated with the AGN and could instead be a result of for example intense star-formation.
We explore two scenarios for the two-component models of the ISM; in these scenarios, the heating mechanisms are represented by either two PDRs with varying stellar radiation field intensities, or a combination of a PDR and XDR component.In the former case, one PDR component is modeled as log(G 0 ) = 2.0, 3.0, and 4.0 to represent a stronger starburst region while the other component is modeled as log(G 0 ) = 0.0.In the latter, the XDR is modeled as in Section 4.2 at distances of 1.0 kpc and 1.5 kpc from the emitting region, and the PDR is modeled as log(G 0 ) = 0.0.In both cases the PDR, log(G 0 ) = 0.0, component is expected to contribute more to the lower-J transitions and the scaled PDR or XDR to the higher-J transitions.The latter of these models is of particular interest as it is physically motivated by the high SFR found in Kade et al. (2023) which, combined with the clear AGN signatures in the SED (Kade et al. 2023), suggests that there is both AGN and star-formation activity influencing the ISM of the quasar.We require the PDR component (modeled as log(G 0 )= 0.0) to have n < 4.5 cm −3 , physically motivated in our models to represent regions of the ISM where the gas physical properties are not as extreme 6 .We expect this component to contribute primarily to the lower-J transitions.Additionally, we limit the scaled PDR and XDR component to have T kin ≤ 150 K.Although this temperature range goes above the typically expected range for giant molecular clouds (GMCs), we include the possibility of gas with an intrinsically higher temperatures which could be due to non-radiative processes such as Fig. 8. CO SLEDs using a two-component description of the ISM in which one component is represented by PDRs scaled to log(G 0 ) = 0.0 (dashed blue line) and the other either PDRs scaled to log(G 0 ) = 2.0, 3.0, and 4.0 or an XDR at distances from the line emitting region of 1.0 kpc and 1.5 kpc (dashed red line).The solid black line shows the best-fit composite model and the lighter grey lines show other models within 1σ.We note that although the best-fit models appear to have been normalized to the CO(12-11) line, no such normalization was performed.
shocks from intense star-formation, turbulence, or cosmic rays (e.g., Ao et al. 2013).Given the companion sources, signatures of on-going star formation, and outflow detection, it is likely that these non-radiative mechanisms are also at work in the quasar.We show the best-fit models for these scenarios in Fig. 8; given the high number of models when combining all scenarios for the two components we consider good models to be within 1σ.The results of our analysis show that two-component model of the ISM provides a better fit to the observed CO emission than a single-component model.
The limited number of CO lines observed in BRI 0952 does not allow for an interpretation of the exact heating mechanism in the warm component.Purely in terms of the best χ 2 fit, the best model is composed of a combination of PDR and XDR heating where the XDR is located at a distance of 1.0 kpc from the emitting region.In this case the physical conditions of the two components are T kin,PDR = 125 K, log(N CO,PDR /cm −2 ) = 17, T kin,XDR = 50 K, and log(N CO,XDR /cm −2 ) = 18.75.However, given the very small difference in χ 2 between this model and the best-fit pure PDR two-component model of χ 2 = 1.33 for the former and χ 2 = 1.36 for the latter, either PDR or AGN heating can reproduce the CO SLED.We find that the best fit model for all five different two-component models have very similar physical conditions, suggesting that these conditions may be indicative of the intrinsic properties of the gas in BRI 0952.However, the wide range for each parameter that still provides a fit within the χ 2 ≤ 1σ range (see Table E.1) suggests that the current data is not capable of providing tight constraints on the true intrinsic properties of the gas.Further, without additional constraints specifically regarding the turnover point of the CO SLED and the behavior of the CO SLED at J > 12 a preferred distinction cannot be made between these heating mechanisms.However, the partial wing detection of , not included in the radiative transfer modeling given the tentativeness of the detection, suggests that there is still significant emission out to J upper = 16 (see Section 5.2) and therefore the favored scenario would be that of a combination of XDR and PDR heating as those models show significant emission out to J upper = 15.The implications of this are further discussed in Section 5.3.
In order to better understand the similarity in best-fit physical parameters between different two-component models, we investigate the likelihood distribution for each parameter twocomponent regime.We follow the approach of Ward et al. (2003), Kamenetzky et al. (2011), andGonzález-Alfonso et al. (2021) for calculating the likelihood distributions for each model given a vector, a, containing the different model parameters.The Bayesian likelihood of a specific model given a set of measurements, M, and parameters a is given by, P(a|M, σ) = P(a)P(M|a, σ) where σ is the measurement error and P(a) is the prior probability density function.The probability density function, P(M|a) , is given by the following, Here we define the prior probability function P(a) as P(a) = 1 for T kin and ∆L in the log(G 0 ) = 0.0 component and for the scaled PDR or XDR we define P(a) = 1 for ∆L and n(H 2 ) and as the following for T kin , We show the relative likelihood for each parameter (n(H 2 ), ∆L, and T kin ) in both the scaled PDR or XDR component and the PDR, log(G 0 ) = 0.0, component in Fig. 9.It is clear that the preferred kinetic temperature and zone gas volume density for both components remains relatively constant across the different heating scenarios.The only parameter that shows variation across the different heating scenarios is the zone width, wherein the log(G 0 ) = 0.0 PDR component favors lower zone widths in the XDR regime but is relatively flat for the PDR+PDR case.The small discrepancies between different heating scenarios provide further evidence that there is little statistical difference between scenarios; however, we maintain that the most physically motivated description is that of a combination of XDR and PDR heating mechanisms.

Modeling of [C ii] and [C i]
We extend out modeling of the species detected to the [C i] and [C ii] line detections using the same MOLPOP-CEP models as used for the CO emission.We evaluate the optimal models derived from various two-component scenarios explored for the CO emission, taking into account models with χ 2 < 1σ, shown in Fig. 10.Generally, the best-fit models from the CO underestimate both the [C i] and [C ii] emission.We find that only the PDR+PDR model, composed of one component with log(G 0 ) = 0.0 and one with log(G 0 ) = 4.0, reproduces the [C i] and [C ii] emission within the margin of error.We note that in all twocomponent scenarios, the scaled PDR or XDR component contributes more to the emission than the log(G 0 ) = 0.0 PDR component, indicating that the majority of the [C ii] and [C i] emission comes from the more intense heating source.In the case of the [C ii] emission, this suggests, given that [C ii] emission traces colder and more diffuse gas including the neutral ISM, the bestfit CO models likely do not describe the intrinsic properties of the gas responsible for the [C ii] emission.Indeed, the low level of emission from the PDR log(G 0 ) = 0.0 component also suggests that applying the best-fit CO model directly to the [C ii] does not provide a reasonable model of the intrinsic conditions responsible for the [C ii] emission in BRI 0952.Regarding the [C i] emission, the larger contribution from the scaled PDR or XDR component is more expected as [C i] is generally thought to trace the warmer molecular gas (see Section 3.4).However, the best-fit CO models also do not reproduce the [C i] emission with the exception of the single aforementioned model.However, as we show in Fig. 10, the additional models within χ 2 < 1σ reproduce the [C i] and [C ii] emission.This is not surprising given the large range in physical conditions that fall within 1σ for the CO emission, as given in Table E 2021) suggests that very high densities and radiation fields would be necessary to reproduce the ratios found in this work.This supports our findings from MOLPOP-CEP suggesting that intense PDR radiation is the only model able to reproduce the [C ii] and [C i] emission detected in BRI 0952; however, we note that the radiation field strengths from Pensabene et al. (2021) necessary to reproduce the ratio found here are significantly higher than those we utilized in our MOLPOP-CEP modeling.Overall, these results suggest that the [C ii] and [C i] = 2.0, 3.0, and 4.0 or an XDR at distances from the line emitting region of 1.0 kpc and 1.5 kpc (in bold in the legend).These models correspond to those shown in Fig. 8. emission is associated with on-going star formation in the quasar host galaxy rather than being associated with AGN activity.

The H 2 O -IR luminosity relation
Sub-mm and far-infrared observations of water with telescopes such as Herschel, ALMA, and NOEMA have shown that water can be a powerful tool for probing the ISM properties of galax-ies (e.g., González-Alfonso et al. 2010;van der Werf et al. 2011;Falstad et al. 2015;Yang et al. 2016;König et al. 2017;Liu et al. 2017).Indeed, studies have shown that water lines can be, in some cases, of equal strength as mid-J CO lines (Yang et al. 2013;Omont et al. 2013).In addition, there has been mounting evidence that there is a correlation between H 2 O and infrared luminosity (e.g., Yang et al. 2013Yang et al. , 2016;;Omont et al. 2013;Jarugula et al. 2021;Pensabene et al. 2022).This correlation is due, in part, to the connection between the excitation of certain water transitions and infrared pumping (González-Alfonso et al.Liu et al. (2017).As these conditions are the same as those required to excite CO(7-6), it is not surprising that the line profile of the water line is quite similar to that of the CO(7-6) line as both likely originate from similar regions within the galaxy (Fig. 5).
We make a very brief analysis of the conditions under which this water line forms.We evaluate the temperature of the dust using two different methods under the assumption that the dust is the main heating mechanism responsible for exciting this molecule.First, we consider the dust temperature from the SED as the continuum measurements in this paper mainly probe the colder dust.Secondly, we compare our SED-estimated dust temperature to MOLPOP-CEP modeling of the H 2 O line where the H 2 O emission is modeled using different dust temperatures.We note that Kade et al. (2023) fit an SED; however, the dust temperature was not a parameter fit in the previous iteration of the SED.As an alternative to estimate the temperature of the colder dust, we use the python package (dust emissivity, Robitaille et al. ( 2021)) to fit a modified black body function where the modified blackbody function is of the form, where S ν is the flux density, h is the Planck constant, τ ν is the opacity at frequency ν, T is the temperature, and κ is the Boltzmann constant.The functional form of the equation does not including heating from the CMB, however including the background contribution from the CMB in the form B ν [T = T CMB(z) ] does not alter the output dust temperature; at z = 4.432 it is unlikely that the CMB would have a large affect on the shape of the SED (e.g., da Cunha et al. 2013b) and therefore we chose not to include it in our modified blackbody fit.This functional form allows for variation in the dust temperature and the β value.We explored the parameter space of both of these values using a Markov chain Monte Carlo (MCMC) approach in which the dust temperature was allowed to vary between 14.742 K and 120 K where the lower limit is set by the CMB temperature at z = 4.432 and β was allowed to vary between 1.3 and 2.5.Both are typical parameter ranges for high-redshift galaxies (e.g., Beelen et al. 2006;Leipski et al. 2013;Pensabene et al. 2022).The best-fit parameters returned from this method were T dust ∼ 45 K and β ∼ 2.0, with a large margin of error.We note that the dust temperature from the SED fitting cannot be well constrained with the available data, as can be shown in Fig. 11 by the wide range of temperature and beta values that produce a reasonable fit.
In the case of the MOLPOP-CEP models, we consider the only heating mechanism to be the dust temperature and assume an opacity of τ = 0.37 .We model the H 2 O 2 11 − 2 02 line using three different dust temperatures, T dust = 25, 50, 75 K and using the same grid of parameters for ∆L, T kin , and n(H 2 ) as for the CO (Section 4).We evaluated the models using the same methodology as for the CO.We show the best-fit model and other good models within 1σ in Fig. 12.We find that all of the explored dust temperatures are able to reproduce the observed H 2 O flux, however, correlating the kinetic temperature, volume density, and zone width with modeling from Liu et al. (2017), which suggests that H 2 O (2 11 − 2 02 ) comes from regions with n(H) ∼ 10 5 − 10 6 cm −3 and T dust ∼ 40 − 70 K, giving an indication that T dust ∼ 50 K may be the most physically motivated model.This is in close agreement with the dust temperature from the SED of T dust ∼ 45 K.However, this analysis and limited by a lack of continuum measurements and a single H 2 O detection.Thus the dust temperature should therefore only be seen as a first estimate.
We compare the ratio of the H 2 O/L IR to the L IR for BRI 0952 to other low-and high-redshift galaxies in Fig. 13.We distinguish between sources with a known AGN contribution and those without.The majority of the sample of high-z AGN are extremely red quasars from Scholtz et al. (2023;in prep), shown in unfilled blue triangles, and the L IR should be considered as an upper limit.We show BRI 0952's position relative to other high-and low-redshift galaxies using both the star-formation dominated infrared luminosity and the infrared luminosity which includes the contribution from the AGN.For sources with reported far-infrared luminosities we correct these to L IR as L FIR ∼ 0.75× L IR (e.g., Decarli et al. 2017).Due to the discrepancies between sources and samples in reporting IR vs. FIR luminosities, we include a conservative indicative error of 30% on the L IR .Kade et al. (2023) decomposed the L IR of BRI 0952 into AGN versus star-formation components, respectively and we consider both components when comparing to other low-and high-redshift sources.We find that BRI 0952 follows the L H 2 O vs. L IR trend from Yang et al. (2013) only when using the decomposed L IR SF .When using the L IR AGN , the quasar lies distinctly below the trend line.However, it should be noted that there is a relative lack of detections of the H 2 O (2 11 − 2 02 ) line in high-redshift AGN.BRI 0952 lies with other AGN-dominated sources from Scholtz et al. (2023;in prep).Considering IR-pumping with the current SED decomposition, the 101 µm photons needed to pump the H 2 O 2 11 − 2 02 line could originate from either AGN or SF activity.However, given the placement of BRI 0952 in relation to other high-redshift quasars and star-forming galaxies and models of the H 2 O 2 11 − 2 02 line, there are strong indications that the line emission is correlated with star formation rather than AGN activity. 1 by assuming that the peak line flux is that which is measured in the spectra (∼ 2.0 mJy) and a line width of 250 km s −1 , consistent with the observed lower-J CO lines.We note that the partial detection is not robust enough for further analysis; however, we include the line as a lower limit when comparing the CO SLEDs as it is indicative of highly excited emission in BRI 0952.We limit our comparison sample to those galaxies that have detections of CO(5-4) in order to remove as much bias from the normalization as possible; however, J1148+0875 is normalized to CO(6-5) and APM 08279+5255 is normalized to an inferred value for the CO(5-4) transition by averaging the adjacent CO transitions (namely the CO(4-3) and CO(6-5) transitions).The flux of the CO(12-11) line in BRI 0952 is elevated compared to the average of the class III galaxies studies in Rosenberg et al. (2015) in which extreme heating mechanisms (AGN or mechanical processes) were present, and significantly higher than class I and II sources.The CO SLED of BRI 0952 follows that of the median CO SLED of galaxies with mid-IR AGN contributions of ≥ 50% from Kirkpatrick et al. (2019); however, it should be noted that Kirkpatrick et al. (2019) point out that the statistical difference between the CO SLEDs of galaxies with mid-IR AGN contributions < 50% and those with ≥ 50% is minimal.Comparisons between the CO transitions detected in BRI 0952 and those reported in Kirkpatrick et al. (2019) are also limited as the median CO SLED from Kirkpatrick et al. (2019) does not extend beyond J upper = 7.
The flux of the CO(12-11) transition observed in BRI 0952 is significantly higher than other samples which extend out to the J upper = 12 transition.One exception is the quasar APM 08279+5255 at z = 3.9 (Weiß et al. 2007).The CO SLED of BRI 0952 broadly follows the shape of the SLED of APM 08279+5255; however with only three robustly observed CO transitions in BRI 0952 it is not possible to determine if the CO SLED turns over between J = 7 and J = 12 as is the case for APM 08279+5255, where Weiß et al. (2007) show that a two component model including an XDR of the ISM is a best-fit to the CO SLED.On the other hand, while the CO(5-4) emission detected in BRI 0952 roughly match those of the quasar J2310+1855 at z = 6.003 (Li et al. 2020), the SLED of J2310+1855 does not exhibit highly excited high-J transitions and the flux of the CO(12-11) emission is significantly lower than that observed in BRI 0952.Regarding the turnover point of BRI 0952, without the partial detection of CO(16-15) the CO SLED could match the shape of either APM 08279+5255 or J2310+1855 but the partial detection of such a high-J CO line raises uncertainty about the turnover point and suggests that the shape may follow that of the quasar J1148+5251 at z = 6.4where the model for that CO SLED also required an XDR region (Gallerani et al. 2014).Without additional CO detections between J upper = 7 and J upper = 12 as well as at J upper > 12, precise conclusions about the shape of the CO SLED of BRI 0952 cannot be drawn.We highlight again the uncertainties of correcting the CO line emission for magnification using the [C ii] emission, as discussed in Section 3.3.It is possible that the CO(12-11) emission has a higher magnification factor, depending on the location of the emission, than the [C ii] emission, causing the CO(12-11) flux to be lower and the turnover point of the CO SLED to shift to lower-J transitions.However, This change would have to be significantly different to affect the conclusions of this paper.

The ISM of BRI 0952
In a commonly accepted paradigm of massive galaxy evolution, massive galaxies merge together causing burst in star formation and a subsequent transition to a quasar, whereby the galaxy then turns into a elliptical galaxy (e.g., Hopkins et al. 2008).Throughout this process the star formation in the galaxy is 'quenched' by either feedback processes, be it from mechanical processes such as intense star formation or AGN feedback (e.g., Fabian 2012;Ishibashi & Fabian 2016).In high-redshift galaxies, it is difficult to quantify the exact contribution of the underlying heating mechanisms present in the ISM.Indeed, it is also challenging to place constraints on the precise evolutionary stage of individual galaxies.This analysis of both archival and high-resolution data allows for an in-depth analysis of a variety of probes of the ISM in the quasar BRI 0952.However, as we have shown throughout the paper and discuss further here, even with these probes many facets of the ISM remain unclear.
BRI 0952 is a radio-quiet quasar (Omont et al. 1996) which has been shown through SED decomposition to host a dominant AGN contribution (Kade et al. 2023).Star formation tracers including dust (IR luminosity) and [C ii] emission suggest that the quasar also hosts intense star formation of a few ×100 − 1000 M ⊙ yr −1 .The three CO line detections reported here give the CO SLED a distinctive shape.The CO(5-4) and CO(7-6) emission are comparable, but the CO(12-11) emission is brighter than either of the lower-J transitions.This, coupled with the partial detection of , is indicative of an extreme underlying heating mechanism within the quasar.MOLPOP-CEP modeling required a two-component model of the ISM in order to reproduce the observed CO emission where the preferred scenario included a PDR and XDR contribution, suggesting that both star formation and the AGN play an active role in heating the gas of the quasar, but without additional CO line detections we cannot come to strong conclusions about individual effect of these mechanisms.
As an indication of AGN activity we consider the tentative OH + emission.In the quasar Mrk231, strong OH + emission was an indication of AGN activity (van der Werf et al. 2010) and has been detected in the high-redshift quasar J1148+5251 (Gallerani et al. 2014).Radiative transfer modelling from González-Alfonso et al. ( 2013) of NGC4418 and Arp220 suggested that the primary production mechanism for H + , a precursor to OH + , was through XDRs or cosmic rays rather than FUV ionizing photons.This result lends strength to the AGN in BRI 0952 being responsible for the OH + emission.
The total molecular gas mass, inferred from the [C i] (2-1) luminosity, suggests that the quasar has a lower total molecular gas mass (of a factor ∼ 100 in all cases) compared to higher-redshift star-forming or star-bursting galaxies such as HFLS3 (z = 6.34 Riechers et al. 2013), SPT0311-58 (z = 6.9 Jarugula et al. 2021), AzTEC-3 (z = 5.3 Riechers et al. 2020), and Mambo-9 (z = 5.85 Casey et al. 2019).The total molecular gas mass is also lower than the value found for quasars and star-forming galaxies at z ∼ 3 − 4 (e.g., Lu et al. 2018;Gururajan et al. 2022).The relatively lower molecular gas mass, coupled with the high starformation rate, suggests that the quasar could be using up its star-forming fuel supply at a tremendous rate.
We investigate the possibility that BRI 0952 has begun to quench its star formation by calculating the depletion time of BRI 0952.The depletion time of a galaxy, defined as t dep = M gas /SFR, is the length of time for which a galaxy can maintain its observed rate of star formation.The depletion time of the local main-sequence galaxies has been extensively studied (e.g., Saintonge et al. 2011;Tacconi et al. 2018Tacconi et al. , 2020)), and there exists a main-sequence trend for local galaxies, although there is ongoing debate about the evolution of depletion time with redshift.We calculate the depletion time of BRI 0952 using both the L IR and [C ii] SFRs and find t dep L IR ∼ 8 Myr and t dep [CII] ∼ 25 Myr.This is significantly lower than DSFGs at z ∼ 3, those found for the z = 6.9 DSFGs studied in Jarugula et al. (2021), and the depletion timescales found for the sample in Stanley et al. (2021).This relatively low depletion time provides an indication that a quenching process may be operating within the quasar.
BRI 0952 is also known to have two companions detected in [C ii] emission within 2.5 arcseconds of the northern image of the quasar where faint velocity gradients between both companions and the quasar have been noted (Kade et al. 2023).In addition, the high ratio of [C i] /CO(7-6) of ∼ 4.8 suggests that the source may be merger driven as lower ratios are associated with disc-like rotating galaxies.The lensing-corrected spectral line profiles of the CO(7-6) and CO(12-11) emission appear complex and the velocity gradient across the [C ii] line profile from lensing (Fig. 7) provide indications of the kinematics operating within BRI 0952, be it regular rotation or merger activity.The [C ii] emission line profile exhibits a broad component indicative of an outflow rate of Ṁout = 98 ± 19 M ⊙ yr −1 ; although it should be noted that the outflow cannot be confirmed without observations of a robust outflow tracer such as P-Cygni absorption line profiles in for example the OH 119 µm or OH + (1 1 − 0 1 ) lines (e.g., Spilker et al. 2020;Riechers et al. 2021).Considering the [C ii] outflow as a source of mass loss, the depletion time would be shortened with respect to the value reported above.(Gallerani et al. 2014), and APM 08279+5255 (Weiß et al. 2007;Bradford et al. 2011).J1148+5251 has been normalized to the CO(6-5) transition.APM 08279+5255 has been normalized to an inferred value for the CO(5-4) transition by averaging the adjacent CO transitions (namely the CO(4-3) and CO(6-5) transitions).We show BRI 0952 including a lower limit on the CO(16-15) from the partial line wing detected in the band 7 data.
Given the high SFR, low gas mass and thereby low depletion time, highly excited CO SLED, outflow evidence, and companions we postulate that BRI 0952 is in a transitionary stage between a galaxy dominated by star-formation activity to one dominated by AGN activity.However, we note the limitations of the data used in this analysis.Although we have a number of observations of atomic and molecular transitions observed towards BRI 0952, these observations are not standardized and suffer from angular resolution limitations, particularly when extrapolating between the data sets.Additionally, our analysis of the CO SLED has shown that without further detections of CO between J upper = 7 and J upper = 12 as well as confirmation of the CO(16-15) emission, it is not possible to robustly determine the heating mechanisms at work in the quasar.Hence, a standardized analysis of CO in BRI 0952 would be necessary to leverage the full power of the CO SLED for the quasar and coming to strong conclusions about the contribution of star formation versus AGN activity to the excitation would likely require an extension beyond CO.

Conclusions
We report detections of CO(5-4), CO(7-6),  1.We further improve the lens model from Kade et al. (2023) using [C ii] observations.We perform lensing modeling using the code Visilens Spilker et al. (2016) and obtain a magnification factor per channel across the [C ii] line which we use to correct the spectral line profiles of the [C ii] , CO(5-4), CO(7-6), and CO(12-11) emission.The broad component in the lensing corrected [C ii] line profile remains clear, confirming the likely presence of an outflow and indicating that lensing effects are not the cause of this component.The lens modeling suggests a velocity gradient across the line profile of the [C ii] emission that could be indicative of ordered rotation or originate from either interactions with the companions found in Kade et al. (2023).2. We detect H 2 O (2 11 − 2 02 ) in the quasar and confirm that the quasar exhibits a similar trend to those previously found in the literature regarding a correlation between L IR and H 2 O line luminosity.We use the L IR from star formation and that from AGN and show that in the case of L IR SF the source lies in a similar region to local AGN and star-bursting galaxies (such as ULIRGs) whereas in the case of L IR AGN the source falls within the region of high-z AGN-dominated sources, but below the trend found for previous star-formation dominated sources.This suggests that the H 2 O emission likely originates from star-forming regions of the ISM. 3. We detect [C i] (2-1) emission in the quasar, allowing for a calculation of the total gas mass within the quasar.Through this approach we find an H 2 gas mass of (6.2±4.9)×1011 M ⊙ , similar to previously reported values for the quasar.We calculate the depletion time using the SFR from both L IR and [C ii] and find t dep L IR ∼ 8 Myr and t dep [CII] ∼ 25 Myr, respectively.4. We use the radiative transfer modeling code MOLPOP-CEP to model the excitation conditions of BRI 0952, specifically with regard to the CO SLED.We find that a twocomponent model of the ISM gives the best fit to the observed data, specifically a combination of PDR and XDR heating where the XDR is located at a distance of 1.0 kpc from the emitting region with physical conditions T kin,PDR = 125 K, log(N CO,PDR cm −2 ) = 17, T kin,XDR = 50 K, and log(N CO,XDR /cm −2 ) = 18.75.However, the precise heating mechanism cannot be determined with the current data as both stellar heating and XDRs can reproduce the current observations.5. We compare the CO SLED of BRI 0952 to other local and high-redshift sources and find that the CO SLED appears to increase up to , but the turnover point is unclear with only the three observed transitions.We note that the CO(12-11) transition is bright compared to other detections, but could follow the trend of other high-redshift quasars where  represents the peak of the SLED.
Combined, we suggest that BRI 0952 is in a transitionary evolution phase, moving from a starburst-dominated phase to an AGN-dominated one.The improved resolution and technical capabilities of ALMA will allow for increased studies focusing on a variety of atomic and molecular species, as well as the kinematics and environment of high-redshift galaxies and eventually allow for a full decomposition of the true influence of star formation versus that of AGN activity.CO SLEDs showing different heating mechanisms including CMB heating, PDRs (log(G 0 ) = 0.0, 1.0, 2.0, 3.0, 4.0), and AGN heating via X-rays at 1.0 kpc and 1.5 kpc.The solid black curve is the best-fit model and the lighter grey lines show the models within the lowest 1% of the χ 2 values for the different models.Note that the method we employ for isolating the best-fit models is a percentile cut, and it is for this reason that a different number of models are shown for different heating mechanisms.

Fig. 1 .
Fig. 1.ALMA bands 3, 4, 6, and 7 continuum images created using line free channels.The contours are shown at −3, 3, 4, 5, 6, 7, 8, 9, and 10σ levels.The purple 'X's show the position of the companion sources, Comp-N and Comp-SW, detected in Kade et al. (2023) as well as the position of Img-N and Img-S, the two images of the quasar BRI 0952.The synthesized beam is shown in the bottom left of each image.

Fig. 2 .
Fig. 2. Moment-0 and moment-1 maps of atomic and molecular species detected in BRI 0952.Rows one and three show the moment-0 maps for the entire field.The purple 'X's show the position of Img-N and Img-S of the quasar as well as the location of the two companions, Comp-N and Comp-SW, detected in Kade et al. (2023).Contours are at −3, −2, 2, 3, 4, 5, 6, 7, 8, 9, and 10σ levels.The black ellipses show the region of spectral extraction and the velocity range used for the extraction is show at the top of the image.Rows two and four show the moment-1 maps within the region of spectral extraction where the flux is limited to only that above a 2σ level.The synthesized beam is shown in the bottom left of each image.

Fig. 3 .
Fig.3.Right: Lensing magnification corrected spectra, corrected using the process described in Section 3.3.The gray background profile is the spectra before lensing correction.Similar to the left panel, single Gaussian fits are shown in red and double Gaussian fits are shown in blue.The dashed gray line represents the rms.The bottom panel shows the residuals from the Gaussian fits where the red corresponds to single fits and the blue to fits using two Gaussian profiles and the gray shaded region represents the rms.The OH redshift-axis is calibrated to the 1834.74735GHz line.Note that the red residuals in the bottom panel are shown to be narrower than the blue only the increase clarity.

Fig. 4 .
Fig. 4. Left: Spectra of molecular and atomic species detected in BRI 0952, not corrected for gravitational lensing.The top panel shows the observed spectra.Single Gaussian fits are shown in red and double Gaussian fits are shown in blue.The dashed gray line represents the rms.The bottom panel shows the residuals from the Gaussian fits where the red corresponds to single fits and the blue to fits using two Gaussian profiles and the gray shaded region represents the rms.The H 2 O and OH + spectra have been binned by a factor two in frequency/velocity.
Fig. B.1.Higher sensitivity data would be necessary to confirm this detection.

Fig. 5 .
Fig. 5. Comparison of the normalized H 2 O, CO(5-4) and CO(7-6) spectra.The CO lines are shown in pink and purple and the H 2 O line is shown in blue.The [C i] line appears at ∼ 400 km s −1 in the CO(7-6) spectra, and similarly the CO(7-6) line appears in the [C i] line profile at ∼> 750 km s −1 .The line profile of the H 2 O line broadly follows that of both CO lines as well as the [C i] line.

Fig. 6 .
Fig. 6.Magnification factor as a function of velocity across the [C ii] line.The different models are shown in the legend which corresponds to those described in Section 3.3.1.It is clear that each model reproduces approximately the same results across the line where the magnification factor peaks at velocities redder than the systemic velocity of the [C ii] line.

Fig. 7 .
Fig. 7. Best-fit position and half-light radius of BRI 0952 across the [C ii] line in the image plane from model 4 (Section 3.3.1).The black dashed lines show the critical and caustic lines.
used the radiative transfer code MOLPOP-CEP (Asensio Ramos & Elitzur 2018) to model the line species detected in BRI 0952.MOLPOP-CEP allows for exact calculations of multilevel line emission for any atomic or molecular line in the Leiden Atomic and Molecular Database database (Elitzur & Asensio Ramos 2006; Asensio Ramos & Elitzur 2018), and has the added advantage that the code is able to easily incorporate radiation fields from AGN, such as XDRs.MOLPOP-CEP utilizes a slab geometry where the line emitting region is divided into a number of zones as a function of optical depth where the input physical parameters can be either uniform or vary throughout the slab.The physical parameters of the slab are as follows: (1) zone width ∆L; (2) gas volume density within the zone n(H 2 ); (3) kinetic temperature T kin : (4) molecular abundance X mol ; and (5) local line width (which corresponds to the line absorption/emission profile at each point in the geometry).
.1, leaving the intrinsic physical conditions of the [C i] and [C ii] emitting regions unconstrained.Previous works by Venemans et al. (2017), Novak et al. (2019), and Pensabene et al. (2021) have utilized the [C ii] /[C i] line luminosity ratio as a diagnostic tool to distinguish between PDR and XDR heating regimes.We find [C ii] /[C i] ∼ 100, well above the expected ratio in the case of XDR heating found by Venemans et al. (2017) and Pensabene et al. (2021) (see Fig. 6 and Fig. 8, respectively).In the PDR case, modeling from Pensabene et al. (

Fig. 9 .
Fig. 9. Relative likelihood of each input parameter to MOLPOP-CEP in different heating scenarios.The left column shows the relative likelihood of each parameter in the PDR component modeled by log(G 0 ) = 0.0 (in bold in the legend) and the right shows the same for PDRs scaled to log(G 0 ) = 2.0, 3.0, and 4.0 or an XDR at distances from the line emitting region of 1.0 kpc and 1.5 kpc (in bold in the legend).These models correspond to those shown in Fig. 8.

Fig. 10 .
Fig. 10.Best-fit MOLPOP-CEP CO model applied to the [C i] and [C ii] emission.The different two-component model is given in the x-axis of the plot, the squares show the absolute lowest χ 2 model from the CO emission (i.e., the black lines in Fig.8), the triangle shows the contribution from the scaled PDR or XDR model contributing to the high-J transitions, likewise the circle shows the contribution from the PDR log(G 0 ) = 0.0 component, and the lightly colored circles show models within 1σ of the lowest χ 2 .Note that the error on the magnification factor is not propagated through to the error on the intensity.

Fig. 11 .
Fig. 11.Top: Modified black-body SED of BRI 0952 used to obtain the dust temperature.The best-fit model returns a dust temperature of T dust ∼ 45 K and β ∼ 2.0.Bottom: Corner plot from the MCMC fitting showing the 16th, 50th, and 84th quartiles.

Fig. 12 .
Fig. 12. MOLPOP-CEP modeling of the H 2 O (2 11 − 2 02 ) line as a function of different modeled dust temperatures.The grey circle, square, and triangle show the absolute best fit model, respectively.The smaller colored circles show the 1σ distribution.Note that the error on the magnification factor is not propagated through to the error on the intensity.
, [C i] (2-1), OH 2 Π 1/2 (3/2 − 1/2), H 2 O (2 11 − 2 02 ), a tentative detection of OH + , and a partial wing detection of CO(16-15) in the highredshift quasar BRI 0952 at z = 4.432.We update the lensing model from Kade et al. (2023) using the [C ii] line emission to extend the lensing model to the other line detections reported in this paper.We have performed an analysis of the ISM of the quasar using the radiative transfer tool MOLPOP-CEP to investigate the heating mechanisms at work.Our conclusions are provided below.

Fig
Fig. D.1.CO SLEDs showing different heating mechanisms including CMB heating, PDRs (log(G 0 ) = 0.0, 1.0, 2.0, 3.0, 4.0), and AGN heating via X-rays at 1.0 kpc and 1.5 kpc.The solid black curve is the best-fit model and the lighter grey lines show the models within the lowest 1% of the χ 2 values for the different models.Note that the method we employ for isolating the best-fit models is a percentile cut, and it is for this reason that a different number of models are shown for different heating mechanisms.

Table 1 .
Details of Observations Band a Date of Obs.b N ant

Table 3 .
Line Properties from our observations.