Physical constraints on the extended interstellar medium of the z=6.42 quasar J1148+5251: [CII] 158um, [NII] 205um and [OI] 146um observations

We report new Northern Extended Millimeter Array (NOEMA) observations of the [CII], [NII] and [OI] atomic fine structure lines and dust continuum emission of J1148+5251, a z=6.42 quasar, that probe the physical properties of its interstellar medium (ISM). The radially-averaged [CII] and dust continuum emission have similar extensions (up to $\theta = 2.51^{+0.46}_{-0.25}\ \rm{arcsec}$, corresponding to $r= 9.8^{+3.3}_{-2.1}\ \rm{kpc}$ accounting for beam-convolution), confirming that J1148+5251 is the quasar with the largest [CII]-emitting has reservoir known at these epochs.Moreover, if the [CII] emission is examined only along its NE-SW axis, a significant excess ($>5.8\sigma$) of [CII] emission (with respect to the dust) is detected. The new wide--bandwidth observations enable us to accurately constrain the continuum emission, and do not statistically require the presence of broad [CII] line wings that were reported in previous studies. We also report the first detection of the [OI] and (tentatively) [NII] emission lines in J1148+5251. Using Fine Structure Lines (FSL) ratios of the [CII], [NII], [OI] and previously measured [CI] emission lines, we show that J1148+5251 has similar ISM conditions compared to lower--redshift (ultra)-luminous infrared galaxies. CLOUDY modelling of the FSL ratios exclude X--ray dominated regions (XDR) and favours photodissociation regions (PDR) as the origin of the FSL emission. We find that a high radiation field ($10^{3.5-4.5}\,G_0$), high gas density ($n \simeq 10^{3.5-4.5}\, \rm{cm}^{-3}$) and HI column density of $10^{23} \,\rm{cm^{-2}}$ reproduce the observed FSL ratios well.


INTRODUCTION
Luminous quasar activity is a key process of galaxy evolution. Indeed, massive outflows driven by the radiation pressure generated by the accretion of gas onto the ular, when they probe the early phase of co-evolution between the first galaxies and their central black hole. The advent of large optical and infrared surveys has enabled the discovery of quasars up to z ∼ 7.5 (Bañados et al. 2018;Yang et al. 2020;Wang et al. 2021), with several hundreds at z > 6 (see Bosman 2020, for an upto-date list). These early quasars harbour SMBHs with M • 10 8 M and accrete gas at or near the Eddington limit for most of their life, challenging models of SMBH formation and growth (e.g., De Rosa et al. 2014;Bañados et al. 2018;Mazzucchelli et al. 2017;Wang et al. 2021). Because the bright quasar light outshines that of the host in the optical and near-infrared, the galaxies hosting early luminous quasars have remained relatively mysterious until the advent of modern (sub-)millimeter observatories.
Numerous observations of z > 6 quasars targeting the bright Far-Infrared (FIR) [C II] 158 µm emission line have revealed their host galaxies to be infrared luminous, dusty and actively forming stars with estimated rates of 10 2 − 10 3 M yr −1 (e.g., Walter et al. 2003Walter et al. , 2009aMaiolino et al. 2005Maiolino et al. , 2012Bañados et al. 2015;Decarli et al. 2018;Venemans et al. 2018Venemans et al. , 2020Novak et al. 2019Novak et al. , 2020Wang et al. 2013Wang et al. , 2019Yang et al. 2020). [C II] 158 µm kinematics also show that most quasar hosts are massive galaxies (M * ∼ 10 10 M ) displaying a variety of morphologies such as stable disks, bulge-dominated galaxies and mergers with nearby companions (Shao et al. 2017Wang et al. 2013Wang et al. , 2019Decarli et al. 2019a,b;Neeleman et al. 2019Neeleman et al. , 2021. The reports of broad [C II] 158 µm line wings in the z=6.42 quasar J1148+5251 (Maiolino et al. 2012;Cicone et al. 2015) spurred the use of the [C II] 158 µm emission line to identify quasar outflow signatures in the early Universe. Such features have however remained rare, and stacking analyses have led to contradictory results (Bischetti et al. 2019;Novak et al. 2020). Recently, Izumi et al. (2021a,b) reported broad [C II] 158 µm line wings in two low-luminosity quasars at z = 6.72 and z = 7.07.
Two decades of [C II] 158 µm and CO studies have shown that early luminous quasars provide an unparalleled observational window into the physics of the earliest (and most massive) galaxies in the Universe. However, multi-line studies using lines other than [C II] 158 µm and CO transitions have been much rarer until now. Since different fine structure lines (FSL) trace different gas densities and excitation levels, only in combination can they probe the ionized and neutral atomic gas phases, and the excitation source(s) of the gas (e.g., Carilli & Walter 2013 369 µm , all accessible at z ∼ 6 with (sub-)millimeter arrays such as the Atacama Large Millimeter Array (ALMA) or the Northern Extended Millimeter Array (NOEMA). Moreover, these lines have been observed with Herschel in large samples of local (ultra)-luminous infrared galaxies -(U)LIRGs, which can which can be readily compared to z ∼ 6 quasars hosts (e.g., Díaz-Santos et al. 2017;Herrera-Camus et al. 2018a).
Detections of FSL other than [C II] 158 µm in z > 6 quasars are still relatively recent Hashimoto et al. 2019;Novak et al. 2019;Li et al. 2020), and a complex picture is emerging from these first results. Emission lines probing the neutral phase ([O I] 145 µm and [C I] 369 µm ) show good agreement between the line ratios and line-to-far-infrared (FIR) ratios of distant quasars and local (U)LIRGs (Novak et al. 2019;Li et al. 2020). Whilst Novak et al. (2019) and Li et al. (2020) report potentially high [O III] 88 µm to [C II] 158 µm ratios (∼ 2) in J1342+0928 and J2310+1855, respectively, this is not the case for the quasar J2100-1715 and its companion galaxy  which have ratios similar to the average of the local population of LIRGS and AGNs (∼ 0.1 − 1; Herrera-Camus et al. 2018a). These first results are however difficult to interpret since the origin of FSL is not always clearly determined and can be linked to different phases (as is the case for the [C II] 158 µm emission line). Clearly, more FIR multi-line studies of z > 6 quasars are needed to understand the ISM of their host galaxies.
SDSS J1148+5251 is one of the earliest high-redshift quasars discovered in the SDSS survey (Fan et al. 2003, z=6.4189), and harbours a 3 × 10 9 M SMBH (Willott et al. 2003). Being the redshift record-holder for many years after its discovery, it was extensively observed with the Very Large Array and the IRAM Plateau de Bure Interferometer (PdBI) and was the first object detected in CO and [C II] 158 µm at z > 5 (Walter et al. , 2004Bertoldi et al. 2003b,a;Maiolino et al. 2005;Riechers et al. 2009;Walter et al. 2009a,b). These pioneering studies probed the host galaxy Star Formation Rate (SFR), dust and ISM properties. Additionally, Maiolino et al. (2012) and Cicone et al. (2015) reported the presence of a broad [C II] 158 µm emission (σ v = 900 km s −1 ) component in the PdBI data, suggesting the presence of an outflow as well as spatially extended [C II] 158 µm emission (up to r ∼ 30 kpc). In this paper, we return to J1148+5251 with a new set of NOEMA observations, capitalizing on larger bandwidths and more antennas, thus improving on the image fidelity as compared to earlier PdBI observations. The new observations targeted atomic fine structure emission lines ([O I] (Riechers et al. 2009) that trace the neutral/molecular gas. Thanks to the wide spectral coverage of the new NOEMA correlator PolyFix and the upgraded NOEMA array, these observations achieved a high fidelity that resulted in deep [C II] 158 µm observations and tight constraints on the underlying dust continuum.
The structure of this paper is as follows. We present in Section 2 and 3 the new observations of the [C II] 158 µm , [N II] 205 µm and [O I] 146 µm emission lines in J1148+5251 as well as the FIR continuum observed between 200 and 280 GHz. We focus on the [C II] 158 µm emission line to re-assess the evidence for a broad velocity component and investigate its spatial extension in Section 4. In Section 5, we derive ISM properties from the strength of the atomic fine-structure emission lines observed, before concluding our study in Section 6. Throughout this paper, we assume a concordance cosmology with H 0 = 70 km s −1 Mpc −1 , Ω M = 0.3, Ω Λ = 0.7. At the redshift of the target (z = 6.42), 1" corresponds to 5.62 proper kpc.
Imaging and cleaning was performed using the latest version of MAPPING/GILDAS (jan2021a). The dirty maps were obtained from the visibilities without tapering and using natural weighting. The data were not primary beam corrected 3 . Cleaning was performed down to 2σ (where σ is the rms noise in the dirty map) using a circular clean region of radius r = 5". The reason for choosing such a wide radius are the earlier reports of extended [C II] 158 µm emission (Maiolino et al. 2012;Cicone et al. 2015). An additional clean region with radius r = 2" was added on the NW source reported by Leipski et al. (2010) and Cicone et al. (2015) which is also detected in the NOEMA data 4 . The final products were created using the following procedure. Firstly, data cubes with 50 km s −1 channels were produced to search for significant emission lines. The [C II] 158 µm , [N II] 205 µm , [O I] 146 µm , H 2 O (5 23 − 4 32 and 3 22 − 3 13 , ν rest = 258.7 GHz), CO(14-13) and CO(13-12) emission lines were fitted with a single Gaussian profile in order to estimate their FWHM. Continuum maps were created using all channels at least 1.25×FWHM away from each emission line. The continuum, determined from the line-free channels using an order 1 interpolation (GILDAS UV BASELINE routine) over the ∼ 7.6 GHz sidebands, was subtracted in the uv plane to create continuum-subtracted cubes.
In order to determine the significance of emission lines, velocity-integrated emission line maps ("line maps") were created by integrating channels over 1.2 times the FWHM of the [C II] 158 µm line (i.e., 482 km s −1 ) at the redshifted frequency of the line. We use the [C II] 158 µm redshift (z = 6.4189) to determine redshifted frequency of all lines. Such maps, assuming the line is Gaussian, contain by definition 84% of the total flux 5 (see for a short derivation appendix A of Novak et al. 2020) . All total line fluxes measured from the line maps and reported in this paper are accounting for this effect. Additionally, all continuum and line fluxes in this paper were computed using the residual scaling method (e.g., Jorsater & van Moorsel 1995;Walter & Brinks 1999;Walter et al. 2008;Novak et al. 2019).
In order to determine the aperture needed to recover most of the flux of the emission lines and the dust continuum, a curve of growth approach was adopted. We show in Appendix A that all line and continuum fluxes reach a maximum or plateau at an aperture radius r = 3", which corresponds to 16.9 kpc at z = 6.42. A nominal aperture radius of 3" is thus adopted throughout the paper. The line subtraction procedure described above was repeated using the r = 3" aperture to obtain final datacubes. Additionally, we have investigated the use of Multiscale cleaning in appendix B. We conclude that the [C II] 158 µm emission and the 259 − 272 GHz continuum are better recovered using Multiscale cleaning, and we therefore use Multiscale cleaning with Gildas/MAPPING for these data throughout the paper.

Dust continuum emission
We first present the FIR continuum maps in Figure 1. The FIR continuum is clearly detected in all four sidebands. The measured continuum flux densities are tabulated in Table 1. Due to the upgraded bandwidth of the NOEMA PolyFix correlator, these continuum measurements have higher sensitivity than previous observations at ∼ 260 GHz (Walter et al. 2009a;Maiolino et al. 2012;Cicone et al. 2015 Figure 1, effective bandwidth 7.68 GHz (masking the edges and central baseband gap) and integrated in a circular aperture with radius r = 3".
Assuming optically thin dust emission at λ > 40µm (e.g., Beelen et al. 2006), we use a modified black body model for the dust emission and correct both for contrast and CMB heating as prescribed by Da Cunha et al. (2013). The dust mass is derived assuming an opacity κ νrest = κ ν0 (ν rest /ν 0 ) β with ν 0 = c/(125µm) and κ 0 = 2.64m 2 kg −1 following Dunne et al. (2003) with β being the dust spectral emissivity index. Our purpose is primarily to measure the FIR luminosity to constrain the SFR in J1148+5251. Therefore, we omit data points at ν rest > 1000 GHz (λ rest 125 µm) where contamination by the quasar non-thermal and torus emission becomes significant (e.g., Leipski et al. 2010Leipski et al. , 2014. The dust SED model uses three free parameters (total dust mass M d , dust emissivity index β, dust temperature T d ) and is fitted using MCMC with the emcee package (Foreman-Mackey et al. 2013) The resulting bestfit and observational constraints are shown in Figure 2 and the posterior probability distribution of the dust SED parameters is displayed in Appendix C. The median dust mass is 3.2 × 10 8 M , the median dust SED index is β = 1.77 +0.30 −0.26 , and the median dust temperature is moderately high (T d = 51.3 K), in agreement with earlier studies (e.g., Beelen et al. 2006;Leipski et al. 2010;Cicone et al. 2015), which is not surprising considering that most of the constraining power comes from observations at lower and higher frequencies than those reported in this work.  The dotted circle centred on J1148+5251 indicates the aperture radius of 3" adopted throughout this paper for flux density estimates. An additional source, reported previously by Leipski et al. (2010) and Cicone et al. (2015) is visible in the three higher-frequency bands at ∼ 10" to the north-west (black dotted circle). The contours are logarithmic (−4, −2, 2, 4, 8, 16, 32)σ (rms). The colour scaling is log-linear, the threshold being at 3σ.  We integrate the modified black body to derive the total infrared (IR, 8 − 1000 µm) and far-infrared (FIR, 42.5 − 122.5µm, e.g., Helou et al. 1985) luminosities. The total infrared luminosity is L IR = (20.9 ± 6.8) × 10 12 L and the far-infrared luminosity L FIR = (13.4 ± 2.4) × 10 12 L , in agreement with earlier studies of J1148+5251. Assuming that dust heating is dominated by young stars, the infrared luminosity can be converted to a SFR using the Kennicutt (1998) and Kennicutt & Evans (2012) conversions, giving SFR= (1830 ± 595) − (2090 ± 680)M yr −1 , respectively. This is in agreement with earlier studies which found that J1148+5251 is in an intense starburst phase (Maiolino et al. 2005;Walter et al. 2009a;Maiolino et al. 2012;Cicone et al. 2015). We detect [C II] 158 µm and [O I] 146 µm at 42σ and 5.3σ (where the SNR is calculated using the peak surface brightness in the line maps and the pixel rms level).

Fine structure line detections
[N II] 205 µm is only marginally detected (3.7σ) and its peak emission is potentially offset from the rest-frame UV, dust, [C II] 158 µm , and [O I] 146 µm emission. In the central pixel of the [C II] 158 µm emission, [N II] 205 µm is formally undetected. We defer further discussion of the [N II] 205 µm emission line to section 5.2. In our subsequent analysis of the ISM of J1148+5251, we focus on the extended aperture-integrated fluxes. All the fluxes and derived line luminosities are tabulated in Table 2. For comparison, we have also measured the flux density of the [C II] 158 µm emission by fitting the apertureintegrated spectrum with a Gaussian, and show a comparison between our work and earlier studies for various aperture sizes in Appendix D. In summary, we find a good agreement between the line map flux estimates and the spectrum best-fits values, as well as good agreement between previous studies and the data presented here,  Figure 5 shows the total and continuum-subtracted [C II] 158 µm spectra (extracted in a r = 3" aperture), as well as single and double Gaussian fits to the data. We found an emission feature at 262.2 GHz which was subsequently masked to determine and subtract the continuum. An emission line map of the feature does not reveal any convincing emission, and we are not aware of any strong molecular or ionised line which could correspond to this redshifted frequency. We checked that masking this feature or not does not impact any of the results that follow. We find that a single Gaussian describes the [C II] 158 µm emission well both in the total and continuum-subtracted spectra, and that the single Gaussian residuals are consistent with the rms noise (Fig. 5). For the total spectrum, the best-fit single Gaussian and continuum model has a monochromatic continuum flux density of S ν = 4.6 ± 0.2 mJy, [C II] 158 µm observed frequency ν [CII],obs = 256.158 ± 0.008 GHz,  Table 2. Atomic fine structure line flux measurements from the line maps (1.2×FWHMCII = 482 km s −1 ), flux-corrected by 1/0.84. Upper limits are given at the 3σ level. The first 3 lines give luminosity for the total integrated flux and the last 3 for the peak [C II]158 µm surface brightness. All integrated fluxes are derived applying residual-scaling correction, and luminosities are computed following the definitions of Solomon et al. (1997) [ = 2.9 ± 3.0 Jy km s −1 . The χ 2 difference between the single-and double-Gaussian models is only ∆χ 2 = 3.66 for three additional degrees of freedom, which corresponds to a p-value of p = 0.7 or evidence at the 0.52σ level for a broad component. Equivalently, the Bayesian Information Criterion (BIC) difference is BIC 1 − BIC 2 = −12.0, where a lower BIC implies a better model, and we adopt a threshold ∆BIC 12 > 10 for strong significance to prefer the more complex model (e.g. Kass & Raftery 1995). We also perform an Anderson-Darling test on the residuals to find any deviations from normally distributed residuals (with a variance equal to the measured rms squared). We find no significant deviation, for either fit, at the p > 0.15 level. We thus do not find evidence for a broad spectral [C II] 158 µm component in J1148+5251 based on the new NOEMA data alone.
For the continuum-subtracted spectrum, the best-fit single and double Gaussian models yield consistent results within the uncertainties (for more details, including the impact of different aperture sizes, see appendix D), and the residuals are nearly identical (Fig. 5 [OI] 145 m, r=1" 2). The 1σ noise level per channel for a r = 1" aperture is shown in red, and a single Gaussian fit is shown in orange (the best-fit parameters and uncertainties in the upper right corner). a double Gaussian is ∆χ 2 = 2.24, which for an increase of 3 free parameters gives a p-value of 0.46 of rejecting the simple model in favor of a double Gaussian, in agreement with the analysis of the total spectrum. The BIC difference is BIC 1 − BIC 2 = −14.0, e.g. the single-Gaussian model is strongly favoured, and an Anderson-Darling test on the residuals also concludes that they  . Total (top panels) and continuum-subtracted (bottom panels) spectrum of the [C II]158 µm line (black) extracted in a r = 3" aperture using residual scaling and 1-and 2-components Gaussian fits (orange line and dashed-dotted blue, respectively). Thin dashed dark blue lines show the two components of the double Gaussian model. Both models are fitted to the entire frequency range to the exception of the overlap between the two sidebands that is coincident with an atmospheric absorption feature (hatched in grey) and the 262 GHz feature (see text). The right panels show a zoomed-in version of the left plots and the residuals of the two models. The shaded red area in the bottom right panels show the 1σ noise level, computed using the rms noise in each channel in mJy beam −1 with sigma-clipping scaled accordingly for the number of beams in the r = 3" aperture. The median noise level at ν obs > 263 GHz is 2.36 ± 0.14 mJy, whereas a direct rms estimate from the continuum-subtracted spectrum at the same frequencies gives 2.00 ± 0.15 mJy. We find no evidence for a broad [C II]158 µm emission line in the new NOEMA data (see text for details).
are consistent with the measured noise at the p > 0.15 level. We note that, regardless of the aperture chosen or the use of residual scaling, the addition of a broad [C II] 158 µm component to the spectral fit is not preferred by the NOEMA data (see Appendix D). We have checked that the continuum fluxes in the lower and higher frequency sidebands are consistent within ∼ 2%. This ensures that the continuum of J1148+5251, determined from the full ∼ 7.6 GHz sideband, is consistent with that determined solely from the single ∼ 3.8 GHz baseband containing the [C II] 158 µm line. We also show in Appendix E that the choice of the line masking width and the flux calibration differences between the two ∼ 3.8 GHz basebands do not impact the recovered The presence of a broad [C II] 158 µm wing reported in previous works (Maiolino et al. 2012;Cicone et al. 2015) and its absence in the NOEMA data seems puzzling at first. Indeed Cicone et al. (2015) reach a sensitivity in the [C II] 158 µm line comparable to ours (0.46 mJy beam −1 vs. 0.45 mJy beam −1 per 100 km s −1 channel) and have similar angular resolution (1.3"×1.2" vs. 1.97" × 1.59"). However, their observations were performed with a smaller number of antennas (6, corresponding to 15 baselines, where the new observations were done with 8-9, corresponding to 28-36 baselines), using a correlator that provided a spectral coverage (∼ 3.6 GHz) twice as narrow than the new observations (∼ 7.7 GHz) around the [C II] 158 µm line. Hence the discrepancy cannot stem from sensitivity or resolution issues and we expect that due to our slightly larger beam, larger number of antennas and higher imaging fidelity we would be more sensitive, if anything, to extended [C II] 158 µm structures in J1148+5251 compared to earlier data. The wider spectral coverage of the new Polyfix also enables us to better constrain the dust continuum emission and the exact shape of a potential broad [C II] 158 µm emission line.
In Appendix F we discuss in detail the origin of these discrepancies which can be ascribed to a combination of (i) continuum subtraction and (ii) residual-scaling effects. We show that when taking these into account, the previous PdBI and new NOEMA data are consistent. We present the combined data in Appendix F, and find that, using the same approach as used here, broad [C II] 158 µm line wings are only tentatively recovered in the merged PdBI and NOEMA dataset. We finally note that the statistical methodology used here to assess the presence of an outflow is arguably quite conservative. Indeed, with the sensitivity and spectral coverage of the data at hand, only extremely powerful outflows would produce a broad line component whose significance is such that a ξ 2 or BIC criterion would prefer a double component fit over a single-component one. In low-redshift studies, it is common to use multiple diagnostics (PV diagrams, other outflow tracers) to assume the presence of an outflow first and then fit the spectral profile with two components (see, e.g., the discussion in Cicone et al. 2014). In any case, these new NOEMA data, thanks to the improved continuum subtraction and analysis methodology, indicate a much less prominent broad component than previously reported, and so rule out the presence of an extremely strong highvelocity outflow. We have however provided above the 2-Gaussian best-fit results in case future observations provide ancillary evidence in favour of an outflow.

Spatially extended [C II] 158 µm emission
The [C II] 158 µm emission is spatially extended (Fig.  3). In Figure 6 we show the dust and [C II] 158 µm radial profiles to investigate their spatial extension. The dust continuum emission (at ∼ 259 GHz) is as extended (within 1σ errors) as the [C II] 158 µm emission (up to θ = 2.51 +0.46 −0.25 arcsec, at the 3σ level, see Figure 6, leftmost panel), suggesting that both trace star-forming material with physical conditions similar to those that we normally ascribe to the ISM of galaxies rather than an extended [C II] 158 µm halo 7 or outflow. This corresponds to a radius of r = 9.8 +3.3 −2.1 kpc once accounting for beam-convolution. This is twice larger than any of the 27 quasars observed in [C II] 158 µm by Venemans et al. (2020), confirming earlier reports (Maiolino et al. 2012;Cicone et al. 2015) that J1148+5251 is an outlier in terms of [C II] 158 µm emission extension. This result holds as well if the dust extension is measured from other spectral setups (e.g., Fig. 1).
The extension of [C II] 158 µm is however asymmetric, with a prominent NE-SW axis (Fig. 3, although this is less pronounced with a larger channel width, see Fig.  8, first panel). In the second and third panel of Figure  6 we compare the dust continuum and the [C II] 158 µm radial surface brightness profile for the radially averaged case, the NE-SW axis and the NW-SE axis (see Fig. 3). Whilst we find no evidence for an extended [C II] 158 µm halo when averaging radially (or along the NW-SE axis), along the NE-SW axis the [C II] 158 µm emission is significantly more extended than the dust continuum (5.8σ). NE + SW Figure 6. [C II]158 µm (red) and dust continuum emission (blue) profiles for J1148+5251, which both extend up to θ = 2.51 +0.46 −0.25 arcsec at the 3σ level, corresponding to r = 9.8 +3.3 −2.1 (accounting for beam-convolution). The shaded areas show the corresponding 1 sigma noise level. A point source profile is shown in black and the dirty (before cleaning and residual-scaling) profile is shown in dashed-dotted red. The three panels show, from left to right the radial profile integrated, over the full 2π, only in the North-West and South-East quadrant, and finally only in the North-East and South-East quadrant (e.g. along the extended emission seen in Figure 3). In the latter case, the [C II]158 µm emission is more extended than the FIR continuum (5.8σ significance).
We finally explore the extension of the [C II] 158 µm emission and the dust continuum directly in the uv plane. To that end, we use the UV FIT and UV CIRCLE routines in GILDAS/MAPPING to fit the visibilities with a point source and 2D Gaussian emission model, and then bin the modelled and observed visibilities radially. We plot the real part of the visibilities against the uv radius in Figure 7. In such plots, a point source gives a constant flux density at all uv radii, while a Gaussian emission model yields a Gaussian profile centered at r = 0. We find good agreement between the observed visibilities and a composite emission model comprising a point source and a 2D Gaussian. We fit both a circular Gaussian and an elliptical Gaussian to the dust and [C II] 158 µm continuum visibilities. For the [C II] 158 µm emission, the best-fit model gives corrected velocity-integrated flux of 4.2 ± 0.5 Jy km s −1 in the point source component, and 5.0 ± 0.5 Jy km s −1 in the extended Gaussian, in good agreement with the flux derived from the line map (Table 2) and the fitted spectrum (see Appendix D). The elliptical Gaussian models yields velocity-integrated fluxes of 5.0 ± 0.3 Jy km s −1 and 5.5 ± 0.4 Jy km s −1 for the point source and resolved components, respectively. In either case, we are consistent with the large fraction of emission coming from the unresolved component reported in Cicone et al. (2015) but not with the overall flux, which is possibly due to continuum subtraction differences as discussed in Appendix F. Additionally, the unresolved flux is consistent with that measured at higher resolution data (Walter et al. 2009a), suggesting the difference in total flux is due to the extended component being resolved out in high-resolution configurations.
The FWHM of the [C II] 158 µm Gaussian component is FWHM = 1.6" ± 0.2" (9 ± 1 kpc), or minor/major axes a = 3.6"±0.3" (20±2 kpc), b = 1.8"±0.2" (10±1 kpc) for the elliptical model. This is in agreement with the scale of the 3σ extension directly measured on the cleaned image. The elliptical Gaussian has a major/minor axis difference and angle (PA=53 ± 4 deg) in agreement with the asymmetric [C II] 158 µm emission discussed above. However, a likelihood ratio test does not prefer the elliptical Gaussian model over the circular one (p-value = 0.83, e.g. 0.96σ significance). We do not find any difference between the elliptical Gaussian and circular Gaussian model for the dust emission (see Fig. 7). The best-fit model for the dust emission gives a flux density of 1.0 ± 0.7 mJy and 3.3 ± 0.7 mJy for the unresolved/resolved components, with the sum in agreement with the measurement aperture-integrated flux density (Table 1). The dust uv profile (Fig. 7, bottom panel) does not show clear evidence for a flattening at large uv distances that is characteristic of a point source. The Gaussian component has a FWHM of r = 0.7" ± 0.1" (3.9 ± 0.6 kpc).  (we show in Appendix F how we recover these structures in the previous PdBI data, although at a lower significance level). Such features are not recovered in the new NOEMA data when averaging over a similar velocity interval (Fig. 8, left, see also Appendix F). The [C II] 158 µm emission is concentrated mostly within 3" (∼ 11kpc) no significant emission is recovered beyond that radius 8 . We further show in Appendix G that no significant velocity structure is detected in J1148+5251 in the NOEMA data, neither by inspecting the channel maps (with channel width ∆v = 120 km s −1 ) nor by means of a kinematical analysis. Following Weiß et al. (2003Weiß et al. ( , 2005, in the optically thin limit, the total mass of ions X emitting photons with rest-frame frequency ν ij stemming from a transition between the i and j levels can be derived from the observed luminosity as where k is the Boltzmann constant, c the speed of light, h the Plank constant, m X is the mass of a single ion X, Q(T ex ) = i g i e −Ti/Tex is the partition function of the species, with g i the statistical weight of level i,  the energy of (above-ground) level i, T ex is the excitation temperature of the i → j transition, and L [Xij ] is the observed integrated source brightness temperature of the line in K km s −1 pc −2 . Note that Eq.  (2) where Q(T ex ) = 2 + 4e 91.2/Tex is the [C II] 158 µm partition function. The optically thin limit assumption, while widespread in the literature, is uncertain at higherredshift since optical depths measurements, although pointing to moderate values, are scarce (e.g., Neri et al. 2014;Gullberg et al. 2015). Lagache et al. (2018) and Vallini et al. (2015) note that in the optically thick limit, only emission from PDRs would reach the observer. As we will show in Section 5.4, most of the [C II] 158 µm emission in J1148+5251 comes indeed from PDRs, and hence we do no expect the optically thin limit assumption to affect our results. Assuming an excitation temperature T ex = 50K (from the CO modelling, Riechers et al. 2009), we find M C + = (5.8 ± 0.3) × 10 7 M . This is five times the neutral carbon mass (M C = 1.1 × 10 7 M ) measured by Riechers et al. (2009).

[N II] 205 µm emission
We report a [N II] 205 µm marginal detection in J1148+5251 at SNR= 3.7 and an integrated (r = 3") luminosity L [N II] = (0.4±0.2)×10 9 L . This is in slight tension with the 3σ limit (< 0.4 × 10 9 L ) reported by Walter et al. (2009b) using the IRAM 30m telescope. Nonetheless, we would expect better image fidelity with the new NOEMA facility. We consider the [N II] 205 µm detection marginal (3.7σ) and urge caution in interpretations that rely upon it.
The marginal (1.59"/8.9kpc ∼ 1 beam) offset between [N II] 205 µm and [C II] 158 µm in J1148+5251 could be of interest. On the one hand, spatial offsets between low and high-ionization lines have been reported in other high-redshift quasars, most often between [C II] 158 µm and [O III] 88 µm (e.g., Novak et al. 2019). Spatial offsets have also been predicted in theoretical simulations (Katz et al. 2017(Katz et al. , 2019 where they arise from different gas phases with different temperatures and densities. Although nitrogen and carbon have a similar ionization level, Katz et al. (2019) show that the [N II] 205 µm and [C II] 158 µm can arise from different gas phases, with [C II] 158 µm originating in lower temperature and higher density regions 9 , which could explain the offset of [N II] 205 µm towards the outskirts of J1148+5251. On the other hand, offsets often seen at low-resolution in highredshift galaxies can disappear in higher-resolution data once fainter emission components are detected (e.g. see HZ10 in Pavesi et al. 2016Pavesi et al. , 2019. Following Ferkinhoff et al. (2010Ferkinhoff et al. ( , 2011, we derive the minimum H + for the [N II] 205 µm luminosity observed. To do so, we assume high densities and a high temperature (as found around O and B stars), such that all nitrogen in the HII regions is ionized. We use Eq. 1 to derive the mass of N + in J1148+5251 from the observed luminosity, with A 10 = 2.1 × 10 −6 the Einstein coefficient of the 3 P 1 → 3 P 0 transition, g 1 = 3 the statistical weight of the 3 P 1 emitting level, ν 10 = 1461.1 GHz the restframe frequency, and g t 9 the partition function. We can then derive the minimum H + mass by assuming the upper limit on the ionized nitrogen to ionized hydrogen ratio (i.e. χ(N + ) = N + /H + ) to be the total nitrogen abundance ratio χ(N ), such that (3) We adopt the abundance value for HII regions ξ(N ) = 9.3 × 10 −5 from Savage & Sembach (1996). Hence we estimate an ionized hydrogen mass M (H + ) ≥ (2.1 ± 1.0) × 10 9 M (2σ level). Using the H 2 gas mass from the CO luminosity (Riechers et al. 2009), the ionizedto-molecular gas ratio is M (H + )/M (H 2 ) > 0.1. This is significantly higher than what is found by Using again Eq. 1 for the [O I] 146 µm 3 P 0 → 3 P 1 transition (T 0 = 329 K, g 0 = 1, Q(T ex ) = 5 + 3e −329/Tex + e −228/Tex ), we derive a neutral oxygen mass Assuming an excitation temperature T ex = 50 K (from the CO modelling, Riechers et al. 2009), we find M O = (9 ± 5) × 10 8 M . This estimate is extremely sensitive to the excitation temperature, e.g. ranging from 0.7 × 10 7 M at 200 K to 8 × 10 8 M at T ex = 50K. Note that the above expression only applies in the optically thin limit which is probably not the case of [O I] 146 µm , therefore underestimating the neutral oxygen mass.

CLOUDY modelling of the fine-structure line ratios
We now use the FSL ratios to determine the physical properties of the ISM of J1148+5251. In order to do so, we make use of CLOUDY (Ferland et al. 2017), a spectral synthesis code designed to simulate the spectra of astrophysical plasmas used to study both low and high redshift galaxies sub-mm lines. The grid of models used in this work were generated to study both highredshift quasar hosts and their companions in Pensabene et al. (2021), which we refer to for further details. To summarise briefly, the model grids includes both PDR and XDR predictions for total hydrogen column densities N H /[cm −2 ] = 10 23 − 10 24 , total hydrogen number density 1 < log n/[cm −3 ] < 6, and local FUV radiation field 2 < log G/[G 0 ] < 6 in units of Habing flux (for the PDRs) or X-ray flux −2 < log F X /[ergs −1 cm −2 < 2] < 2 (for the XDRs). For each combination of parameters, the flux of various lines of interest is predicted and the ISM conditions can be determined by comparing to the observed line ratios.
Whereas ratio is primarily a function of the N + /C + abundance ratio (e.g., Oberst et al. 2006). Photoionization models (e.g., Oberst et al. 2006;Pavesi et al. 2016;Croxall et al. 2017 The ratio of the luminosities measured in this work for J1148+5251 is 27 ± 13 for the aperture-integrated (r = 3") flux and > 100(2σ level) at the peak of the [C II] 158 µm emission. In either case, most of the [C II] 158 µm emission (88% − 97%) comes from the neutral phase. Given the resolution of the NOEMA data, we use the r = 3" aperture-integrated line luminosities to compute line luminosity ratios and compare to the grid of CLOUDY models. In order to include properly the 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 Log    measurements can be explained by the larger aperture, and is comparable to that seen for [C II] 158 µm (see Appendix A). We show the radiation field and density predictions of CLOUDY for the observed FSL ratios in Figs. 10 and 11. We find that the different line ratios are not perfectly reproduced for a single density and radiation field in the chosen grid of models. However, the [C II] 158 µm neutral / [C I] 369 µm luminosity ratio (47 ± 3) excludes XDR models which cannot reproduce such high ratios (the maximum being ∼ 15) (Fig. 10, e.g., Venemans et al. 2017a,b;Novak et al. 2019;Pensabene et al. 2021). Besides, our XDR model grid cannot reproduce all the FSL-to-FIR luminosity ratios with a single hydrogen number density and X-ray flux. We note that this result, based on the analysis of the fine structure lines, is in tension with Gallerani et al. (2014) who concluded that an XDR component is needed in J1148+5251 to explain their reported CO(17-16) line. However, as noted by the authors, this line is potentially contaminated by a nearby OH + emission line, and thus additional observations of high-J CO lines are needed to clarify the situation.
The PDR models in agreement with the observed ratios (except [O I] 146 µm / [C I] 369 µm ) have an HI column density of N HI = 10 23 cm −2 , high radiation fields (10 3.5−4.5 G 0 ) and moderate hydrogen number densities (n 10 3.5−4.5 cm −3 ), commensurable with other studies of high-redshift quasars (e.g., Novak et al. 2019;Pensabene et al. 2021 discrepancy could stem from the different quality and resolution of the data, as well as the low SNR of both lines. Another possibility is that a fraction of the [C II] 158 µm emission could come from collisional excitation in the cold neutral medium, rather than solely from the PDRs as is the case for [O I] 146 µm and [C I] 369 µm , leading to a small discrepancy between ratios involving the [C II] 158 µm line and others in our PDR-only models. This hypothesis is supported by the spatial extent of the [C II] 158 µm line emission (Section 4.2), which is beyond what is normally observed for the disk component of z ∼ 6 massive star-forming galaxies and other quasar hosts (e.g. Venemans et al. 2020 and the underlying dust continuum emission. The high-fidelity data, together with the large instantaneous bandwidths that NOEMA provides, enabled us to revisit the properties of the quasar's host galaxy and derive the physical conditions of its interstellar medium. The main conclusions of this paper are as follows: • A uv plane analysis confirms the presence of an extended [C II] 158 µm and dust emission component (FWHM = 1.6"±0.2" (9±1 kpc)) accounting for ∼ 50 − 60% of the total [C II] 158 µm and most of the dust emission, in agreement with earlier studies (Maiolino et al. 2012;Cicone et al. 2015). J1148+5251 thus remains an outlier with very extended [C II] 158 µm and dust emission compared to other z > 6 quasars (e.g., Venemans et al. 2020).
• The [C II] 158 µm emission shows a similar spatial extent as the FIR dust when averaged radially (up to θ = 2.51 +0.46 −0.25 arcsec, corresponding to r = 9.8 +3.3 −2.1 physical kpc, accounting for the beam convolution). However, if the [C II] 158 µm emission is examined only along its NE-SW axis, a significant (5.8σ) excess [C II] 158 µm emission (w.r.t. to the dust) is detected.
• The [C II] 158 µm line profile can be fitted with a single Gaussian with a FWHM of 408 ± 23 km s −1 .
The new NOEMA data has enabled us to rassess the significance of the broad [C II] 158 µm line wings reported in previous studies. We find that an additional broad [C II] 158 µm component is not statistically required to describe the NOEMA data (or the merged NOEMA and previous PdBI data).
• The FSL line ratios are consistent with the trends observed in local (U)LIRGs. We find that a large fraction (∼ 90%) of the [C II] 158 µm emission originates in the neutral phase (PDR) of the gas.
• We have compared CLOUDY models to the observed FSL ratios in J1148+5251. The [C II] 158 µm /[C I] 369 µm ratio and the multiple FSL-to-FIR ratios exclude XDR models in favour of PDRs. We find good agreement for models that have a high radiation field (10 3.5−4.5 G 0 ), a moderate hydrogen number densities (n 10 3.5−4.5 cm −3 ) and HI column density N HI = 10 23 cm −2 .
The results described here highlight the importance of large instantaneous bandwidths when observing high-redshift quasars (or galaxies) to search for weak extended emission of atomic or molecular lines. Our findings enabled a renewed view on the host galaxy of the J1148+5251 quasar shedding light on the feedback activity and providing new constraints on the excitation conditions of its interstellar medium that appear similar to what is found in local ULIRGs. Higher angular resolution and sensitivity data, in particular in the [C I] 369 µm emission line, would be required to put additional constraints on the gas density and the ionisation source of J1148+5251 and further explore its properties. In this appendix, we investigate what is the size of the aperture that is necessary to recover most of the continuum or [C II] 158 µm line fluxes. Figure 12 shows the line flux density of [C II] 158 µm , [O I] 146 µm and [N II] 205 µm as well as the continuum flux density as a function of increasing aperture radius. We find that all line/continuum fluxes reach a maximum or plateau at an aperture radius r = 3", which corresponds to 16.9 kpc at z = 6.42. The [C II] 158 µm flux presents tentative evidence for additional flux up to 4" (∼ 1 Jy km s −1 ), but this is within the 1σ errors. Note that at large radii where there is no more cleaned flux, residual-scaling can become numerically unstable. Throughout this paper, an aperture of r = 3" is therefore adopted for measurements unless specified otherwise. In this appendix, we briefly detail our experiments with different cleaning methods for our interferometric data. Högbom cleaning (Högbom 1974) is one of the standard method for cleaning interferometric data. It relies on iteratively finding peaks in the data and subtracting the dirty beam at that location until the residuals reach a desired level. It is particularly efficient for point sources, but struggles with large-scale emission which it tries to reconstruct using a multitude of point sources. In that case, so-called multi-scale algorithms which convolve the beam with various Gaussians to subtract larger scales are preferable (e.g., Wakker & Schwarz 1988;Cornwell 2008). Figure 13 shows the clean map, dirty map and residuals for Högbom and Multiscale clean on the [C II] 158 µm map (integrated over 482 km s −1 ). Clearly, the Högbom clean residuals show a flat excess of 2σ flux filling the 3" aperture which would not be expected if the source was a single (or a limited number of) point source(s). On the contrary, the residuals of the Multiscale algorithm are closer to zero on average. Therefore, we chose the Multiscale algorithm for all [C II] 158 µm -derived quantities and images in this paper. Figure 14 shows the clean, dirty and residual maps for the continuum maps and the other FSL maps for a Högbom clean. For the lower frequency continuum maps, [O I] 146 µm , and [N II] 205 µm maps expect the residuals are wellbehaved and do not require the use of multiscale cleaning. For the 259, 274 GHz continuum, some residuals are seen and the multiscale algorithm is adopted for those in the paper, albeit changing only the final continuum flux by 2%.   14.5 ± 0.9 9.41 ± 0.32 30 ± 3 20.1 ± 0.7 23 ± 2 23.5 ± 1.2 34 ± 4 24.7 ± 2.0 Iν [Jy km s −1 ] 5.3 ± 0.5 3.7 ± 0.2 11.0 ± 1.5 8.1 ± 0.4 14 ± 3 10.0 ± 0.8 13 ± 3 11.0 ± 1.0 Table 3. [C II]158 µm line fluxes measured from a Gaussian fit the aperture-integrated spectra in this work and previous studies. For aperture radius and study, we give the best-fit velocity width σv, the peak line flux density S peak and the integrated flux Iν . For the Cicone et al. (2015, C15)  In table 3, we compare the [C II] 158 µm fluxes measured from the extracted aperture-integrated [C II] 158 µm spectra for various aperture sizes between previous studies and this work. This work's spectra and best-fit Gaussians are displayed in Figure 16 for circular apertures with radii r = 1", 2", 4" and a central-pixel-only spectrum. The r = 3" case is already presented in the main text (Fig. 5). For Cicone et al. (2015), we report their best-fit parameters and fluxes for the "narrow" (−200, 200) km s −1 component defined in their study. We find good agreement for all apertures between this work and previous studies. This implies that most of the excess flux in larger apertures (r 3") described in Cicone et al. (2015) is due to the reported presence of blue/red-shifted components with large spatial offsets which are not recovered in this work, as discussed in Section 4.

E. CONTINUUM SUBTRACTION AND MASKING OF THE [C II] 158 µm LINE
In this appendix, we provide details on the continuum subtraction procedures, focussing on the [C II] 158 µm spectral setup. Figure 17 shows the aperture-integrated [C II] 158 µm spectrum for various half-width masking regions ranging  r=4" Figure 16. Continuum-subtracted aperture-integrated [C II]158 µm spectra (black) for circular apertures with radii r = 1", 2", 4". The best fit Gaussian is shown in orange and the best-fit parameters are displayed in the upper-right corner of each plot.
model in favor of a double Gaussian emission profile. Finally, Figure 18 shows the dust continuum measurement as a function of the masking half-width. The continuum is predictably higher for small masking regions, but has converged at the masking half-width adopted in this paper (1.25 × FWHM([CII])).

F. PREVIOUS PDBI DATA AND MERGE WITH THE NEW NOEMA DATA
In this appendix, we analyse the previous PdBI observations of J1148+5251 published in Cicone et al. (2015). We do not recalibrate their data, and follow our analysis method presented in Sections 2 and 3. Based on the previous PdBI data, we show the total spectrum extracted in a r = 3" aperture applying residual scaling corrections (e.g., Jorsater & van Moorsel 1995;Walter & Brinks 1999;Walter et al. 2008;Novak et al. 2019) in Figure 19 (where we also show the new NOEMA spectrum for comparison). Clearly the two datasets are compatible, and a single-Gaussian+continuum fit gives similar continuum levels and [C II] 158 µm peak fluxes within the ∼ 10 − 20% amplitude calibration errors. We note that the noise in the previous PdBI spectrum is not uniform, as different datasets with different frequency coverages and antenna configuration were stiched together (Cicone et al. 2015). This leads to increased noise in the ranges ν 254.7 GHz and ν 258.0 GHz. The best-fit continuum to the total spectrum in the r = 3" (4.6 ± 0.3 mJy) is perfectly consistent with the new measurement (4.5 ± 0.2 mJy) and the best-fit dust SED (Section 3.1), and in tension with the 3.3 mJy continuum value (at 256 GHz) published in Cicone et al. (2015) derived from a best-fit model in the UV plane to line-free channels. We therefore infer that this tension is probably due to an unfortunate choice of the continuum channels or an unstable UV-plane fitting routine in the old pipeline.
We then show in Figure 20 different spectra extracted from the previous PdBI dataset with and without residual scaling. It can be clearly seen that the continuum flux at the edges of the band drops significantly in larger apertures when residual scaling is not used. This subtle effect is due to the fact that the continuum is only detected at the 2σ level in most line-free channels (e.g. at the edges of the band). When cleaning down to a typical threshold 2σ, there is in fact no or little flux above that threshold and therefore no or few clean components are added to the final "clean" map (which is the sum of the clean Gaussian components and the residual map which does not have units of "clean beam" but "dirty beam") for channels that only contain continuum emission. As a consequence, these channels are dominated by the residual map which imprints the dirty beam pattern. With large apertures, e.g. r > 3" in this case, one starts to integrate over negative sidelobes of the synthesized beam, decreasing the integrated flux significantly. This effect does not play a dominant role where the [C II] 158 µm line is present, as in those channels the flux is dominated by the actual "clean" components. This can then lead, when combined with an under-subtracted continuum, to the illusion of a strong broad component in the final spectrum. However, with the residual scaling approach taken into account, this broad component is less significant (although still tentatively detected, see Fig. 21) and the results become compatible with the new NOEMA data.
We now comment on the [C II] 158 µm structure reported in Cicone et al. (2015). We show in Figure 22 the collapsed [C II] 158 µm channels in the range of (−1400, +1200) km s −1 based on the previous PdBI data only. We find a similar structure as previously reported, but at a different significance level (in our case not exceeding the 3σ   NOEMA, r=3", with residual scaling Figure 19. Total spectrum (black) of the [C II]158 µm line extracted in r = 3" apertures using residual scaling in the previous PdBI data (left) and the new NOEMA data (right). Single Gaussian (+continuum) fits to both dataset are shown in orange and the parameters are given in the upper left corner. The best-fit continuum (constant) and single-Gaussian model parameters are consistent between the two datasets. We note that the increase in noise (shown as red line) in the PdBI data at low (< 254.7 GHz) and high (> 258 GHz) frequencies is due to the way the data were frequency-stitched at the time.
level in the extended regions). We therefore speculate that the rms might have been underestimated by the earlier GILDAS/MAPPING software release, which would explain the numerous −6 and −3σ regions in the previously published [C II] 158 µm map. We conclude this appendix by presenting the results of our analysis of the merged PdBI and NOEMA [C II] 158 µm spectrum. We show in Figure   PdBI, r=4", no residual scaling Figure 20. Total spectra of the [C II]158 µm line with (first row) and without (second row) residual scaling extracted in r = 2" and r = 4" apertures (first and second column respectively) in the PdBI data. Single Gaussian (+continuum) fits to both dataset are shown in orange. In the absence of residual scaling, the continuum flux drops away from the [C II]158 µmline in the larger apertures (see text for details). Even though no continuum was subtracted in the UV-plane, the continuum flux approaches ∼ 0 mJy at the band edges in the r = 4" aperture (without residual scaling), creating the illusion of a broad emission line.

G. CHANNEL AND MOMENTS MAPS OF THE [C II] 158 µm EMISSION
In this appendix, we present additional visualizations of the [C II] 158 µm emission velocity structure in J1148+5251. Firstly, we present channel maps in Figure 24 with a channel width of 117 km s −1 . No significant emission (> 3σ) is detected at large radii or at velocity offsets > 400 km s −1 from the peak of the [C II] 158 µm line (cf. Cicone et al. 2015). Secondly, we present the integrated flux, mean velocity and velocity dispersion (so-called moment maps) in Figure 25. These maps are generated using Qubefit  using all [C II] 158 µm voxels detected at SN> 3 and standard parameters. We find no kinematic evidence for a bulge-dominated dispersion or a rotating disk model.  PdBI, r=3", with residual scaling Figure 21. Total spectrum of the [C II]158 µm line (black) presented as in Figure 5, using the previous PdBI data, following the methodology used in this paper. Both single-Gaussian and double-Gaussian fits are consistent with that performed on the NOEMA data only (see Fig. 5). With a ∆BIC12 = 3.0 < 10, a broad [C II]158 µm component is only tentatively detected.  PdBI+NOEMA, r=3", with residual scaling Figure 23. Total spectrum of the [C II]158 µm line (black) presented as in Figure 5, but using the merged PdBI and NOEMA dataset. Both single-Gaussian and double-Gaussian fits are consistent with that performed on the NOEMA data only (see Fig.  5  . The colour scaling is log-linear, the threshold being at 3σ (rms). The colour scaling is log-linear, the threshold being at 3σ (rms).