© Y.-C. Cheng et al. 2022

1 Introduction

Comets are providing interesting clues regarding the history of the Solar System. Through their mineralogical, chemical, and isotopic composition, comet nuclei document environmental conditions and processes occurring from the protostellar collapse phase to the protoplanetary disk phase. The nature and isotopic properties of their ices suggest inheritance from the presolar phase (Altwegg et al. 2019; Drozdovskaya et al. 2019), whereas the presence of high-temperature minerals shows that comets also incorporated materials that formed close to the Sun (Engrand et al. 2016, and references therein).

Abundance ratios of nuclear-spin isomers for cometary molecules having identical protons have been measured in a number of comets as possible cosmogonic indicators of environmental formation conditions (Mumma et al. 1987; Kawakita & Kobayashi 2009). About two dozen measurements have been obtained for water (Faggi et al. 2018, and references therein) and ammonia (Shinnaka et al. 2011, 2016). These two molecules present two nuclear-spin isomers, so-called ortho and para species. Most of the ortho-to-para ratio (OPR) determinations for water have been obtained from near-IR observations of ro-vibrational lines, using high spectral resolution. In the case of NH₃, the OPRs are estimated indirectly from observations of NH₂ in the optical domain. The nuclear-spin temperatures determined from OPRs of both H₂O and NH₃ are ~30 K on average (e.g. Faggi et al. 2018; Shinnaka et al. 2016). The meaning of these low spin temperatures remains debated. Laboratory studies of the nuclear-spin conversion of water during thermal desorption show fast re-equilibration in the gas phase to an OPR of 3 (which is equal to the statistical weight ratio and corresponds to a nuclear-spin temperature higher than ~60 K) (Hama et al. 2011, 2016, 2018; Sliter et al. 2011). However, the equilibrated water ortho-to-para ratio of 3, corresponding to a spin temperature higher than 60 K, is also commonly observed in comets.

In this paper, we present the measurement of the ortho-to-para ratio of water in the coma of comet 67P/Churyumov-Gerasimenko (67P) from observations performed with the high spectral resolution channel (H) of the Visible and Infrared Thermal Imaging Spectrometer (VIRTIS) on board Rosetta (Coradini et al. 2007). Because of the low activity of the comet (water production rate of typically 10²⁸ s⁻¹ at perihelion; e.g. Biver et al. 2019), this measurement has not been possible from ground-based telescopes. During most of the mission, VIRTIS-H acquired 2.5–2.9 μm spectra in which ro-vibrational lines of both ortho and para H₂O variants are clearly present. However, the OPR retrievals are affected by optical depth effects. Hence, a careful evaluation of these effects has been performed using lines from weak vibrational bands in combination with radiative transfer calculations (Debout et al. 2016). In Sect. 2, the observational dataset is presented. The analysis of the spectra and the results are described in Sect. 3. Radiative transfer calculations are presented in Sect. 4, followed by a summary and a discussion in Sect. 5. Results for individual VIRTIS-H data cubes are given in the appendices.

2 VIRTIS-H observations and selected datasets

VIRTIS is composed of two channels: VIRTIS-M, a spectro-imager with a visible (VIS) (0.25–1 μm) and an infrared (IR) (1–5 μm) channel operating at moderate spectral resolution (λ/Δλ = 70–380), and VIRTIS-H, a cross-dispersing spectrometer providing spectra with higher spectral resolution capabilities (λ/Δλ = 1300–3500) in eight orders of diffraction covering the range 1.9–5.0 μm (Drossart et al. 2000; Coradini et al. 2007).

To determine the OPR of water, we used VIRTIS-H spectra in order numbered 4 (see Bockelée-Morvan et al. 2019, for the list of orders), covering the wavelength range 2.432–3.077 μm in 432 spectral elements with a spectral sampling Δλ_samp = 0.0015 μm. VIRTIS-H order 4 has been found to be best suited for the OPR retrieval, as neighbouring orders either suffer from stray-light effects or have a less adequate wavelength coverage. Moreover, the cross-order absolute calibration has additional uncertainties that increase instrumental errors when several orders are in use. The nominal spectral resolving power in order 4 λ/Δλ is ~2200 at 2.67 μm, where Δλ is the grating resolution. However, the effective spectral resolution is lower due to undersampling with respect to the grating resolution (Δλ_samp = 1.25 Δλ). In addition, we used calibrated data cubes where detector odd-even column readout effects (Raponi et al. 2020) were removed by performing spectral smoothing with a boxcar average over three spectral channels, thereby reducing the spectral resolution by another factor. From these considerations, a reduction of the spectral resolving power by a factor ~2.7 with respect to the nominal value is expected (Debout 2015). For accurate OPR determinations through spectral fitting, the good knowledge of the instrument spectral response is a key requirement. By fitting high signal-to-noise ratio (S/N) VIRTIS-H H₂O spectra of the comet with synthetic spectra (Sects. 3.2–3.4) and minimising the χ-square statistics, the effective spectral resolution was determined to be equal to 779 at 2.7 μm (Δλ_eff = 0.00347 μm). The wavelength calibration (found to be off by up to 0.5 × Δλ_samp) was also improved by adjusting the wavelengths of the spectral channels with strong water lines to expected values.

VIRTIS-H acquired data cubes of typically 3–4h duration in various pointing modes. For the present study we considered limb observations during which the instrument stared at a fixed limb distance, with a line of sight (LOS) at an azimuth close to the comet-Sun line for most of them. The instantaneous field of view (FOV) of the VIRTIS-H instrument is 0.58 × 1.74 mrad² (58 m × 174 m for a Rosetta S/C distance to the comet of 100 km). The version of the calibration pipeline for the data is CALIBROS-1.2. The calibration process is detailed in Jacquinod et al. (2015) and summarized in Bockelée-Morvan et al. (2016).

The considered data are (1) dataset S1: a set of 102 data cubes (Table B.1) acquired from 3 June 2015 (heliocentric distance of r_h = 1.506 au) to 13 January 2016 (r_h = 2.120 au) (the perihelion of comet 67P was on 13 August 2015 at r_h = 1.242 au); the mean distances of the FOV to the comet centre (ρ) for these data cubes are between 1.65 and 7.13 km; (2) dataset S2: the average of 5 data cubes obtained at far limb distances (ρ = 171 to 444km) and mean r_h = 1.31 au; (3) dataset S3: 28 data cubes obtained at r_h = 1.561–2.730 au pre-perihelion (mean r_h = 2.292 au) that were averaged; (4) dataset S4: average of four data cubes obtained at r_h = 1.543 au post-perihelion and limb distances ρ = 12.0 km; and (5) datasets H1 and H2: two spectra with high S/N obtained by averaging numerous data with high-intensity H₂O lines acquired near perihelion.

The resulting spectra obtained from datasets S2, S3 and S4 provide OPR measurements for a low H₂O column density along the LOS (<5 × 10¹⁹ m⁻²). Spectra H1 and H2 allow us to measure the OPR ratio from the faint v₁ and hot bands near 2.750 μm and longwards of 2.774 μm. Geometric informations and exposure times for datasets S2–S4 and H1–H2 are provided in Table 1.

Fig. 1 Continuum-subtracted VIRTIS-H spectrum of dataset H1 (blue curve) superimposed on the synthetic model that best fits the HB (2.774–2.910 μm) spectral region (orange curve). Dataset H1 is the weighted average of 22 data cubes obtained near perihelion (see Sect. 3.5.3 and Table 1 for geometric parameters and fitted parameters). The spectral range used to fit the HB region is indicated with a dotted green line. The bottom panel shows a synthetic optically thin spectrum at Trot = 116 K, with different colours for the v1 (green), v3 (grey), and hot bands (red).

3 Fitting method, data analysis and retrievals

3.1 H₂O spectrum in the 2.5–3μm range

A typical H₂O spectrum acquired by VIRTIS-H near perihelion is shown in Fig. 1. The intense lines in the 2.55–2.76 μm range correspond to the v₃ (001–000) fundamental vibrational band. The 2.77–2.91 μm region is dominated by the weak v₁ (100– 000) band, and v₁ + v₃ − v₁ (101–100) and v₂ + v₃ − v₂ (011–010) hot bands (Bockelée-Morvan & Crovisier 1989; Villanueva et al. 2012). These weak bands contribute slightly to the 2.55–2.76 μm water emission. In optically thin conditions, the total H₂O fluorescence emission rate (g-factor) in the 2.5–3.0 μm range at r_h = 1 au is 4.24 × 10⁻⁴ s⁻¹ for a rotational temperature in the ground vibrational state T_rot = 100 K (fluorescence model of Villanueva et al. 2012).

Weak emission lines from OH v(1–0) prompt emission have been observed in this spectral domain in cometary spectra observed by the Infrared Space Observatory (Crovisier et al. 1997, 1999), and ground-based telescopes (e.g. Bonev et al. 2006). The strongest OH line is expected at 2.8 μm, blended with water lines (Bockelée-Morvan & Crovisier 1989). Weak excess emission with respect to expected H₂O emission is observed at this wavelength in VIRTIS-H spectra with a high S/N (e.g. Fig. 1), consistent with the detection of OH prompt emission from comet 67P. Therefore, the 2.7991–2.8018 μm region of the VIRTIS-H spectra was excluded from the fit.

3.2 Synthetic H₂O spectra

Synthetic H₂O spectra were computed using the model developed by Crovisier (2009). This model computes the full fluorescence cascade of the water molecule excited by solar radiation, describing the population of the rotational levels in the ground vibrational state by a Boltzmann distribution at T_rot (the rotational temperature T_rot is expected to be representative of the gas kinetic temperature in the inner coma of 67P). This model assumes optically thin conditions. We verified that excitation by solar radiation scattered by the nucleus and by the thermal radiation of the nucleus can be neglected. The model, which includes fundamental bands and hot bands, uses the comprehensive H₂O ab initio database of Schwenke & Partridge (2000) and describes solar radiation as a blackbody. The synthetic spectra closely resemble those obtained by the model of Villanueva et al. (2012), which uses the BT2 ab initio database of Barber et al. (2006) and includes a more exact description of the solar radiation field including solar lines. However, in the 3 μm region, the excitation of the water bands is not affected by solar Fraunhofer lines (Vil-lanueva, personal communication). When we compare spectra from the models of Crovisier (2009) and Villanueva et al. (2012), the total water emission rates in the 2.5–3 μm region differ by about 1% (T_rot = 50 K) to 4.5% (T_rot = 150 K). The most significant differences occur in the hot bands, whose fluorescence excitation and emission rates are underestimated in the model of Crovisier (2009). The total emission rates in the 2.774–2.910 μm region (referred to as the hot-band domain in Sect. 3.4) are underestimated by 4.4%, 8.4%, and 9.1% for T_rot = 50, 100, and 150 K, respectively. We checked that OPR and T_rot retrievals are similar when the model of Villanueva et al. (2012) is used. First, by fitting synthetic fluorescence spectra from Villanueva et al. (2012) with our model as described in Sect. 3.4, the retrieved OPR and T_rot values are consistent within 1% with the correct values. A second test consisted of modifying the fitting procedure to fit VIRTIS-H spectra with synthetic spectra from Villanueva et al. (2012) instead of spectra from Crovisier (2009): the retrievals are consistent within the uncertainties for all the output parameters of the fitting procedure (see Sect. 3.4).

Figure 2 shows a synthetic H₂O spectrum at the effective spectral resolution of VIRTIS-H. The ortho and para lines are plotted separately.

3.3 Baseline subtraction and data combination

The water lines are observed atop an underlying continuum due to dust scattering and thermal emission (Bockelée-Morvan et al. 2019). This continuum presents significant short-term variations in the time interval covered by the data cubes. In order to remove the baseline of the spectra in the best possible way, we split the observing data into cells of 32 acquisitions (the acquisition exposure time is typically 3 s). The acquisitions constituting each cell were then averaged, and the baseline was fitted by a fifth-order polynomial function considering wavelength ranges with limited H₂O fluorescence emission. All the baseline-subtracted cells were then averaged using a sigma clipping algorithm (with threshold of 3) to dismiss noisy spectral pixels. Figure 3 shows the individual acquisitions for a selected data cube and the final background-subtracted spectrum (in blue). The average spectrum without background subtraction is also shown (orange line).

The standard deviation for each spectral channel was evaluated using five data cubes obtained at limb distances >900 km, which present no detectable H₂O emission. The statistics was performed on continuum-removed average spectra.

Fig. 2 Synthetic optically thin fluorescence H2O spectrum at the effective spectral resolution (and channel spacing) of VIRTIS-H. The blue and orange lines show the ortho and para lines, respectively. The total spectrum is plotted in grey. The ortho and para lines at infinite spectral resolution are shown at the bottom. The ortho-to-para ratio is equal to 3 and Trot = 100 K.

Fig. 3 Illustration of data combination and background subtraction. The grey spectra show the 1728 acquisitions of cube T1-00400765976 acquired on 13 September 2015 from 11h:46m:19.3s to 15h:24m:49.8s UT. The orange spectrum is the averaged spectrum, and the blue spectrum is the average spectrum obtained after subtraction of the continuum background, following the method explained in the text. This cube, which shows extreme variations in dust continuum because it was acquired during an outburst (Bockelée-Morvan et al. 2017), was not considered for this paper.

3.4 Fitting procedure

The fitting procedure consists of fitting the background-subtracted VIRTIS-H spectra with synthetic spectra parametrised by T_rot, OPR, and a normalisation factor, and convolved to the effective spectral resolution of VIRTIS-H (Sect. 2). The 2.590–2.760 and 2.774–2.910 μm ranges, which we refer to as the main-band (MB) and hot-band (HB) domains, respectively, were independently fitted. The fitting procedure also corrects for small baseline residuals. In a first step, the whole fitting procedure was applied to the HB domain (plus the 2.523–2.590 μm region with negligible water emission), using a third-order polynomial function for the baseline. Next, the MB region, without this baseline, was fitted. The derived OPR values are found to be very sensitive to baseline removal, and we optimised the removal in the best possible way.

The outputs of the fitting process are the band intensities and rotational temperatures deduced for the two domains. The OPR is a free parameter when the MB domain is fitted. It is a fixed parameter (equal to 3) when the HB domain is fitted, except for the analysis of the H1–H2 datasets with a high-S/N (Sect. 3.5.3). The column density $N_{{\text{H}}_2\text{O}}^{\text{MB}}\left(N_{{\text{H}}_2\text{O}}^{\text{HB}}\right)$ is inferred from the wavelength-integrated band intensities I^MB (Wm⁻² sr⁻¹) (I^HB) using $$I^{\text{MB}}=\frac{hv{g}_{\text{thin}}^{\text{MB}}}{4\pi }{N}_{{\text{H}}_2\text{O}}^{\text{MB}},$$(1) $$I^{\text{HB}}=\frac{hv{g}_{\text{thin}}^{\text{HB}}}{4\pi }{N}_{{\text{H}}_2\text{O}}^{\text{HB}},$$(2)

where $g_{\text{thin}}^{\text{MB}}$ and $g_{\text{thin}}^{\text{HB}}$ are the sum of the fluorescence emission rates (g-factors) of the individual lines in the MB and HB wavelengths domains, respectively. $g_{\text{thin}}^{\text{MB}}$ and $g_{\text{thin}}^{\text{HB}}$ were calculated from the optically thin model of the water fluorescence (Sect. 3.2) using the rotational temperatures $T_{\text{rot}}^{\text{MB}}$ and $T_{\text{rot}}^{\text{HB}}$ deduced from the fitting process, and scaled according to $r_{\text{h}}^{-2}$.

The presence of optical depth effects can be evaluated through the attenuation factor f_atten (<1 for an optically thick coma), $$f_{\text{atten}}=\frac{N_{{\text{H}}_2\text{O}}^{\text{MB}}}{N_{{\text{H}}_2\text{O}}^{\text{HB}}}.$$(3)

Indeed, in optically thick conditions, the effective g-factor g_eff, which relates the measured band intensity I (W m⁻² sr⁻¹) to the column density $N_{{\text{H}}_2\text{O}}$ along the LOS, is smaller than the fluorescence g-factor g_thin (Debout et al. 2016). Hence, the water column density derived from the MB using Eq. (1) underestimates the actual column density by the factor g_thin/g_eff. On the other hand, the intensity of the v₁ and hot bands contributing to the HB domain, which are typically one order of magnitude fainter than the main v₃ band and are therefore much less affected by optical depth effects (Sect. 4.2), provide a good estimate of the actual column density through Eq. (2). To first approximation, f_atten is therefore equal to the ratio g_eff/g_thin of the whole v₃ band.

Spectral fits to the HB and MB wavelength ranges are shown in Figs. 1, 4–6 for datasets H1, H2, S2, S3 and S4. The inferred parameters are given in Table 1. We note that the poor fit of the MB part of datasets H1 and H2 regarding relative line intensities (Fig. 4) is characteristics of optical depth effects. A good fit is obtained when optical effects are much less significant: HB part of H1 and H2 datasets (Figs. 1 and 5) and MB part of dataset S2, S3, and S4 (Fig. 6). The attenuation of the v₃ band due to opacity is clear on Figs. 1 and 5, especially for the H1 dataset. Spectral fits for selected¹ S1 data cubes are shown in Appendix A.

Fig. 4 Continuum-subtracted VIRTIS-H spectra (blue) superimposed on the synthetic model that best fits the MB (2.590–2.760 μm) spectral region (orange). Top: H1 dataset. Bottom: H2 dataset. The standard deviation for each spectral channel is shown at the bottom of the plots (and vertically offset with respect to the measured radiances for better clarity). The MB spectral region is shown by a horizontal green line.

3.5 Results of the fitting process

3.5.1 Dataset S1: Optically thick spectra

Dataset S1 consists of 102 individual data cubes covering the period from 3 June 2015 (r_h =1.506 au) to 13 January 2016 (r_h = 2.120 au). Most spectra (92 out of 102) have been obtained at azimuth angles with respect to the comet-Sun line PA_LOS of less than 45° (Table B.1). This dataset has been selected for its spectra with detected signal in the HB spectral region with an S/N > 10. During this time interval, the water production rate of comet 67P is seen to vary by a factor of ~30, from typically 3 × 10²⁶ to 8 × 10²⁷ s⁻¹ according to observations carried out with the Microwave Instrument for the Rosetta Orbiter (MIRO) (Biver et al. 2019). Table B.1 provides geometric information for all the individual data cubes of dataset S1, and it lists the results obtained from spectral fitting.

With limb distances spanning a few hundred metres to ~5 km above the nucleus surface, this dataset is a representative sample of the spectra obtained under highly variable LOS water column density (1.8 × 10²⁰–1.9 × 10²¹ m⁻²) and rotational temperature conditions (60–150K) (Fig. 7A). High (small) rotational temperatures are obtained for high (low) column densities, and correspond to small (large) limb distances (Fig. 8A). Rotational temperatures differ only slightly from the (number-density weighted) average of the gas temperature along the LOS (Debout 2015). This decrease in gas temperature with increasing altitude results from adiabatic cooling, and it was also measured by the MIRO instrument (Biver et al. 2019). A temperature dependence on perihelion distance is clearly observed (Fig. 8A).

Figure 7A shows that many water spectra are somewhat optically thick, with minimum f_atten values of ~0.4 (equivalent optical depth τ = 0.9, with τ = –ln(f_atten)). As expected for optical depth effects, there is a clear anti-correlation between f_atten and $N_{{\text{H}}_2\text{O}}^{\text{HB}}$. Because bands with lines in the HB spectral region are one order of magnitude fainter than the v₃ band, they are probably optically thin to first approximation, so that $N_{{\text{H}}_2\text{O}}^{\text{HB}}$ provides the effective water column density along the LOS. A temperature dependence of f_atten is also marginally observed (see the data points for column densities >10²¹ m⁻² in Fig. 7A). This trend is consistent with radiative transfer calculations (Sect. 4). We note that $T_{\text{rot}}^{\text{HB}}$ should be representative of the LOS average of H₂O rotational temperature, whereas $T_{\text{rot}}^{\text{MB}}$, although much more accurately derived, is expected to be higher than $T_{\text{rot}}^{\text{HB}}$ when optical depth effects are significant. Indeed, faint ro-vibrational lines excited from weakly populated high-energy rotational levels are relatively less affected by optical depths than lines formed through low-lying levels. The relative intensities of the lines in the spectra are then modified and mimic a warmer environment. This trend is observed, as is clear in Fig. 8B: the ratio $T_{\text{rot}}^{\text{MB}}/T_{\text{rot}}^{\text{HB}}$ increases with decreasing f_atten, reaching ~1.35 in the most optically thick cases corresponding to low f_atten values. For the highest values of f_atten, the rotational temperature measured from the lines of the v₃ band matches the value measured from the much thinner v1 and hot bands within 10%.

Apparent OPRs derived from MB fitting are strongly data dependent, reaching values as low as 1.2–1.5 for LOS sampling high column densities (Fig. 7B). These low apparent OPRs are due to opacity effects, because lines from ortho-H₂O species are more efficiently weakened. OPR and f_atten follow similar trends with the H₂O rotational temperature and column density (Figs. 7A, B). The decrease in optical depth effects with increasing rotational temperature is clearly visible in Fig. 7B: the trend for higher OPRs at higher rotational temperatures is clear over the full range of the water column density. The strong correlation between the derived OPR and the attenuation factor is shown in Fig. 7C.

At low opacity levels (f_atten > 0.65), that is, at low column densities, there is significant scatter in the derived OPR. The values lie between 2.0 and 2.6. This is presumably related to spectra with low S/N.

Fig. 5 Continuum-subtracted VIRTIS-H spectra (blue) superimposed on the synthetic model that best fits the HB (2.774–2.910 μm) spectral region (orange). Top: H1 dataset. Bottom: H2 dataset. The HB spectral region is shown by a horizontal green line. See the caption of Fig. 4 for other details.

3.5.2 OPRs derived in low column density cases

In the spectra of dataset S1, the H₂O MB is optically thick (f_atten < 0.65), affecting the determination of the OPR. For these spectra, the column density is higher than ~1.8 × 10²⁰ m⁻². Consequently, we attempted to derive the OPR from datasets S2, S3, and S4, for which the column density is lower by a factor 3 to 30, using the MB (the signal in the HB region is too weak to be analysed). The results are given in Table 1. The derived OPR varies between 2.3 (S4) and 2.8 (S2). The low OPR value for S4 (2.26 ± 0.07) suggests that this spectrum is more severely affected by opacity effects. The low-temperature conditions for this set $\left(T_{\text{rot}}^{\text{MB}}=48\text{\hspace{0.17em}K}\right)$ compared to those for set S3 $\left(T_{\text{rot}}^{\text{MB}}=68\text{\hspace{0.17em}K}\right)$, but similar column densities, is consistent with a higher optical depth for S4. The highest measured OPR value (2.79 ± 0.13) is for dataset S2, which samples a very low column density. The attenuation factor for these datasets is estimated in Sect. 4.2 using radiative transfer modelling.

Fig. 6 Continuum-subtracted VIRTIS-H spectra (blue) superimposed on the synthetic model that best fits the MB (2.590–2.760 μm) spectral region (orange). Results for the S2, S3 and S4 spectra. See the caption of Fig. 4 for other details.

Fig. 7 Water column density $\left(N_{{\text{H}}_2\text{O}}^{\text{HB}}\right)$, rotational temperature $\left(T_{\text{rot}}^{\text{HB}}\right)$, apparent OPR (from MB fitting), and attenuation factor (fatten) derived from the fitting of VIRTIS-H spectra by optically thin fluorescence spectra. The 102 spectra of dataset S1 are considered. Plots A–C show how these parameters vary with each other. The colour-coding is indicated at the top of each plot.

Fig. 8 Rotational temperatures derived from fitting the MB (2.59–2.76 μm) and HB (2.774–2.91 μm) spectral regions. (A) $T_{\text{rot}}^{\text{HB}}$ as a function of distance to nucleus centre ρ, with a colour-coding according to heliocentric distance; only LOSs with a position angle with respect to the comet-Sun line PALOS < 45° were considered; results from S3-S4 and H1-H2 spectra are included; the fitted curves are $T_{\text{rot}}^{\text{HB}}=63+182/\rho \left(\text{K}\right)$ (dashed red line for rh < 1.4 au), and $T_{\text{rot}}^{\text{HB}}=35+225/\rho \left(\text{K}\right)$ (dashed green line for 1.4 au < rh < 1.8 au). (B) Ratio $T_{\text{rot}}^{\text{MB}}/T_{\text{rot}}^{\text{HB}}$ as a function of the attenuation factor fatten with symbol colours according to $T_{\text{rot}}^{\text{HB}}$.

3.5.3 OPRs derived from the H₂O hot-band spectral region

Near the perihelion of 67P, lines in the 2.774–2.910 μm HB spectral region are detected with a high S/N (Fig. 1). These intrinsically weak lines from the v₁ and hot bands are thinner than ro-vibrational transitions from the v₃ band, by typically one order of magnitude to first approximation if we scale the opacity according to the relative pumping rates (Crovisier 2009; Villanueva et al. 2012). Hence, OPR values retrieved from the fit of HB spectral region are an alternative way to determine the effective H₂O ortho-to-para ratio for comet 67P.

The H1 dataset averages 22 data cubes that satisfy the following geometric conditions: r_h < 1.4 au, 2.0 < ρ < 4.0 km, and position angle PA_LOS < 10° with respect to the comet-Sun line. The H2 dataset is a combination of 9 data cubes acquired from 30 May 2015 to 1 June 2015 (r_h = 1.52 au, 8.3 < ρ < 9.4 km, and PA_LOS = 3°. The HB domain of the H₂O spectrum is well detected for both sets (cf. Fig. 1 for H1). Fitted parameters are given in Table 1. The OPR determination with the best accuracy gives 2.62 ± 0.16 (H1 dataset), in agreement with retrievals obtained from S2 and S3. The attenuation factor of the HBs is estimated in Sect. 4.2 from radiative transfer modelling.

4 Radiative transfer calculations and 67P water OPR

In the previous section, the VIRTIS-H data were analysed using optically thin synthetic H₂O fluorescence spectra. Variations in retrieved parameters OPR, T_rot, and f_atten with column density and correlations between these parameters have been qualitatively explained by optical depth effects. In this section, the results of radiative transfer calculations are compared with the data and are used to derive the water OPR in the coma of 67P. Our approach consists of computing synthetic spectra with a model considering optical depth effects both for the excitation of H₂O and for the received radiation. These optically thick synthetic spectra are then analysed with a fitting procedure similar to that used for VIRTIS-H data (Sect. 3.4), that is, by adjusting optically thin synthetic spectra to derive the apparent OPR and rotational temperature and the attenuation factor. The approach consisting of directly fitting optically thick spectra to ViRTIS-H measured spectra was not possible, in particular, because of the CPU time requested per model run, the variety of the coma conditions, and the large number of data cubes.

4.1 Radiative transfer model

We used the dedicated radiative transfer model of Debout et al. (2016), which was developed for the purpose of analysing H₂O, CO₂ and CO vibrational spectra of comet 67P from VIRTIS. Only the v₃ band is considered for the H₂O model. The calculations are made for an isotropic coma, with physical parameters (number density, gas kinetic temperature and expansion velocity) derived from fluid models. As shown by Debout et al. (2016), optical depth effects in the solar excitation rate vary not only with the cometocentric radial distance, but also with azimuthal angle due to the unidirectionality of the exciting source (the Sun). They are maximised along the comet-Sun direction and minimised in perpendicular directions. The received radiation is also altered by opacity effects, which are maximum (minimum) when the phase angle (S/C-comet Sun angle) is 0° (90°). These results can be explained by the structure of the coma gas flow.

The main input parameters of the model are the total production rate and the nucleus surface temperature, which determine the radial profiles of the gas kinetic temperature and velocity (respectively decreasing and increasing with the distance to nucleus), and the gas number density (see Debout et al. 2016). We performed calculations for water production rates from 2 × 10²⁶ to 4 × 10²⁸ s⁻¹, heliocentric distances of 2 and 1.3 au and surface temperatures of 284 K (r_h = 2 au) and 353 K (r_h = 1.3 au). The observational geometric conditions we considered are terminator orbits and an azimuthal angle PALOS = 0^°, which correspond to most of the observing circumstances for VIRTIS-H data. Calculations were made for limb distances from 125 m to 20 km above the surface. Hence the set of computed spectra covers a wide range of column density and gas temperature conditions. We considered ortho-to-para ratios of 2.5 and 3. The spectral resolution of the synthetic spectra was 750, that is, similar to VIRTIS-H spectra.

These H₂O synthetic spectra (referred to as RT spectra) were analysed with a fitting procedure similar to that used for the VIRTIS-H data. However, because the RT spectra only include the v₃ band, the synthetic optically thin fluorescence spectra we used for the fitting only included this band. The fitting was performed in the MB spectral region, with $T_{\text{rot}}^{\text{RT}}$ and OPR^RT as output parameters. The attenuation factor $f_{\text{atten}}^{\text{RT}}$ was deduced from the ratio of the integrated intensity of the whole v₃ band to that of the expected intensity under optically thin conditions. We defined $T_{\text{LOS}}^{\text{mean}}$ as the rotational temperature obtained by not considering opacity effects in the radiative transfer model, because it is close to the LOS weighted-average of the gas temperature. $T_{\text{LOS}}^{\text{mean}}$ corresponds to $T_{\text{rot}}^{\text{HB}}$, measured on the VIRTIS data.

Figure 9 shows an optically thick RT spectrum and the fitted optically thin spectrum. The strong similarity with the fit obtained for the MB part of H1 spectrum (Fig. 4, top) is remarkable.

Fig. 9 Optically thick RT spectrum (H2O v3 band) superimposed on the fitted optically thin spectrum. The input parameters for the calculation of the RT spectrum are OPR = 3, $N_{{\text{H}}_2\text{O}}=1.3\times {10}^{21}{\text{m}}^{-2}$, an LOS mean temperature (so-called $T_{\text{LOS}}^{\text{mean}}$) equal to 130 K, PALOS = 0°, and a terminator orbit. The results of the fit are apparent OPR (i.e., OPRRT) of 1.95, attenuation factor $f_{\text{atten}}^{\text{RT}}=0.61$, and rotational temperature $T_{\text{rot}}^{\text{RT}}=157\text{K}$.

4.2 Comparison of RT model and data outputs

The apparent OPR (OPR^RT) is shown in Fig. 10 as a function of water column density and mean LOS temperature $T_{\text{LOS}}^{\text{mean}}$ for the two assumptions of the effective OPR. This figure can be compared directly to the data (Fig. 7), taking into account that the colour-coding is different. The same trends are observed: the values of the apparent OPRs converge to the value of the effective OPR for water column densities significantly below 10²⁰ m⁻². Hence, this confirms that the determination of the effective 67P’s OPR requires using results from the low-column density datasets (Table 1). The extent to which these data are affected by residual opacity effects must be examined.

4.2.1 Direct comparison of optical depth effects in modelled and observed spectra

In the radiative transfer model, the assumption is made that the water coma is isotropic, which is not realistic. The water distribution observed around the nucleus of 67P looks like a broad fan approximately oriented towards the Sun (Biver et al. 2019). Moreover, the assumed velocity field, though physically realistic for an isotropic coma, should deviate significantly from the one in the complex coma of 67P. For these reasons, one can expect optical depth effects to be quantitatively different in the synthetic and real H₂O 67P’s coma. Figure 11 shows that this is indeed the case. We plot in this figure (dashed black line) a fit of $N_{{\text{H}}_2\text{O}}^{\text{MB}}$ versus $N_{{\text{H}}_2\text{O}}^{\text{HB}}$ from dataset S1. This curve is a pair of two linear regressions, one obtained by fitting data with $N_{{\text{H}}_2\text{O}}^{\text{MB}}\le 3.2\times {10}^{20}{\text{m}}^{-2}$ (Eq. (4)) and the other with $N_{{\text{H}}_2\text{O}}^{\text{MB}}>3.2\times {10}^{20}{\text{m}}^{-2}$ (Eq. (5)),

$$N_{{\text{H}}_2\text{O}}^{\text{HB}}=1.761\times N_{{\text{H}}_2\text{O}}^{\text{MB}}-1.666\times {10}^{19},$$(4) $$N_{{\text{H}}_2\text{O}}^{\text{HB}}=2.482\times N_{{\text{H}}_2\text{O}}^{\text{MB}}-2.710\times {10}^{20}.$$(5)

The trend of $N_{{\text{H}}_2\text{O}}^{\text{MB}}$ versus $N_{{\text{H}}_2\text{O}}^{\text{HB}}$ is nicely reproduced by this two-line fit.

In Fig. 11, the results from the radiative transfer model are plotted by blue dots. Opacity effects are underestimated by the model, as is best seen for $N_{{\text{H}}_2\text{O}}>5.0\times {10}^{20}{\text{m}}^{-2}$. The attenuation factor f_atten as measured from the data is about 18% lower than for the model results. Since it is beyond the scope of this paper to change the description of the water coma in the radiative transfer model, we empirically assigned the model outputs to a column density 2.3 times lower. When this factor f = 2.3 is applied, the attenuation factor given by the model (red dots in Fig. 11) becomes consistent with the data (dashed line) for $N_{{\text{H}}_2\text{O}}>5.0\times {10}^{20}{\text{m}}^{-2}$. A slightly different factor, f = 3.0, is found to reproduce the data for $N_{{\text{H}}_2\text{O}}<5.0\times {10}^{20}{\text{m}}^{\text{-2}}$ best. Interestingly the knee observed in the $N_{{\text{H}}_2\text{O}}^{\text{MB}}$ versus $N_{{\text{H}}_2\text{O}}^{\text{HB}}$ correlation is also apparent in the model results (Fig. 11). This behaviour is expected when the band opacity τ, defined as f_atten = e^−τ, varies proportionally to $N_{{\text{H}}_2\text{O}}$.

Fig. 10 Results from radiative transfer calculations. Apparent OPR as a function of LOS mean temperature $\left(T_{\text{LOS}}^{\text{mean}}\right)$ and column density (A, B), and as a function of the attenuation factor (C). Plots A and B refer to calculations with an effective OPR of 2.5 and 3, respectively, with a colour-coding according to the derived apparent OPR. For C plot, the colour-coding is according to the column density (in base-10 log scale).

Fig. 11 Direct comparison of optical depth effects derived from model and observations. The column density (blue dots) obtained by fitting the synthetic RT spectra $\left(f_{\text{atten}}^{\text{RT}}\times N_{{\text{H}}_2\text{O}}\right)$ is plotted as a function of the effective column density $N_{{\text{H}}_2\text{O}}$. The values plotted as red dots are obtained by considering that the radiative transfer model underestimates the opacity effects by a factor f = 2.3 (see text). The dashed line shows the two-line fit of $\left(N_{{\text{H}}_2\text{O}}^{\text{MB}},N_{{\text{H}}_2\text{O}}^{\text{HB}}\right)$ values from dataset S1 (see text). The plain line shows optically thin conditions.

Fig. 12 Opacity $\tau =-\ln \left(f_{\text{atten}}^{\text{RT}}\right)$ as a function of column density (dots). Here, we assume that the radiative transfer model underestimates the opacity effects, and a factor of f = 2.3 is applied (see text), so that the values correspond to the red dots plotted in Fig. 11. Colour-coding is according to the LOS mean temperature $T_{\text{LOS}}^{\text{mean}}$, and the colour scale is given at the right. A polynomial fitting of degree 2 is shown for various temperatures (35, 55, 70, 80, 100, 120, and 140 K). The dashed line corresponds to the two-line fit of $\left(N_{{\text{H}}_2\text{O}}^{\text{MB}},N_{{\text{H}}_2\text{O}}^{\text{HB}}\right)$ values from dataset S1 shown in Fig. 11.

4.2.2 Opacity as a function of column density and LOS temperature

Figure 12 displays the opacity as a function of water column density for the model outputs (dots), applying f = 2.3 as described above, and for the data (dashed line). A colour-coding is used for the model outputs to show how the opacity varies with mean coma temperature along the LOS. The opacity increases with decreasing LOS mean temperature. There is a non-linear correlation between opacity and water column density at fixed temperature; instead the trend with column density is correctly described applying a polynomial function of degree two (Fig. 12). Comparing in this figure the VIRTIS results from dataset S1 (dashed line corresponding to the two-line fit) to the model, the agreement is satisfactory overall. The opacities derived from the data vary weakly (τ in the range 0.5–0.8) with the column density, especially for low column densities, because the rotational temperature in the coma of 67P is lower on average for low column densities (Fig. 7).

4.3 H₂O OPR for 67P: Final steps

4.3.1 Attenuation factors of the S2–S4 main-band region

We used the opacity variation with column density and rotational temperature shown in Fig. 12 (i.e. the polynomial functions shown by the solid lines) to estimate the attenuation factor f_atten for datasets S2–S4. Following previous discussions, calculations were also made for a correction factor f= 3, as the data suggest that at low column densities the opacity is still underestimated with a model using f = 2.3. For example, for column densities ~3.0 × 10²⁰ m⁻², the opacity derived from the data is 0.50 for a rotational temperature of 85 K, whereas the models with f = 2.3 and f = 3.0 give τ = 0.43 and 0.51, respectively. However, the derived f_atten values are not very different for the two models (see the small error bars of f_atten for S2–S4 in Fig. 13). Mean f_atten values are 0.92, 0.85, and 0.79 for S2, S3, and S4, respectively.

Fig. 13 Apparent OPR as a function of the attenuation factor fatten. Results from datasets S1 and S2–S4 are shown as grey and black dots, respectively. Results from datasets H1–H2 are shown by black squares. Only S1 data with an S/N on IHB > 20 are considered. For S2–S4 and H1–H2, the error bars for fatten represent the range of values obtained for f= 2.3 and f= 3. The arrows drawn for H1–H2 indicate a possible underestimation of fatten for these datasets. The curves show the polynomial regressions that fit the OPRRT versus $f_{\text{atten}}^{\text{RT}}$ values shown in Fig. 10C: red and green represent OPR = 2.5 and 3, respectively.

4.3.2 Attenuation factors of the H1–H2 hot-band region

To estimate the attenuation factor for the sum of the emission lines in the HB region (about 50% are from v₁ and 47% are from v₁ + v₃ − v₁ and v₂ + v₃ − v₂), one has to consider that optical thickness has two sources: (i) attenuation of the solar pump (ASP), and ii) absorption of emitted photons (AEP) along the ray path towards the observer. These two effects scale proportionally to g_lu excitation and g_ul emission rates (g-factors), respectively, where u and l are the upper and lower vibrational states of the bands, respectively. Concerning ASP, the weighted average of the g_lu excitation rates of v₁ and hot bands is ~15 times lower than the v₃ excitation rate (Villanueva et al. 2012). As for AEP, the g_ul emission rate for v₁ is ~10 times lower than for v₃ (Villanueva et al. 2012). For the hot-bands, AEP is very weak since the lower vibrational states (e.g., v₁ for the v₁ + v₃ − v₁ hot-band) are weakly populated (Bockelée-Morvan et al. 2004). The attenuation factor for H1 and H2 was calculated from the polynomial functions applying a multiplicative factor of 1/15 to the column density. The derived attenuation factors (~0.88 for H1, ~0.95 for H2) might be slightly underestimated since AEP is weak for the hot bands.

4.3.3 H₂O OPR for 67P

The apparent OPRs as a function of f_atten for all datasets (selecting S1 data with the highest S/N in the HB domain) and for the RT models with OPR = 2.5 and 3 are shown in Fig. 13. For clarity, the results from the radiative transfer model are displayed using the polynomial regression that fits the OPR^RT versus values shown in Fig. 10C (instead of plotting the individual results of Fig. 10C). This figure shows that the data points for datasets S2–S4 and H1–H2 are consistent with the RT model results obtained for OPR =3.0. In order to estimate the OPR in the coma of 67P and its uncertainty, we searched for the polynomial regression that provides the best fit to these data points. Specifically, we assumed that the polynomial regression linking OPR^RT/OPR versus $f_{\text{atten}}^{\text{RT}}$ that we obtain for OPR = 3 (green curve in Fig. 13), is valid for OPR values slightly different from OPR = 3. The water OPR ratio in 67P is finally estimated to 2.94 ± 0.06. It is worth mentioning that the same result was obtained regardless of whether we considered the H1 and H2 data points because their overall weight is low.

5 Summary and discussion

Through spectral fitting, we have investigated the H₂O ortho-to-para ratio in the coma of comet 67P using 2.5–3 μm spectra obtained with the VIRTIS-H instrument on board Rosetta from 3 June 2015 to 13 January 2016. We took advantage of the detection of faint lines from the v₁, v₁ + v₃ − v₁ and v₂ + v₃ − v₂ bands beyond 2.774 μm to assess optical depth effects affecting the strong v₃ band. Most of the spectra are affected by optical thickness, and the derived apparent value for the OPR is found to decrease with increasing opacity. In the most severe conditions (H₂O column densities $N_{{\text{H}}_2\text{O}}>{10}^{21}{\text{m}}^{\text{-2}}$), the intensity of the v₃ band is reduced by a factor of ~2.5 and the apparent OPR is found to be as low as ~1.3. Optical thickness also affects the ground-state rotational temperatures derived from the v₃ ro-vibrational lines: the values are found to be up to 50 K higher than those obtained from spectral fitting of the faint v₁ and hot-bands. Observed trends with H₂O column density and mean gas coma temperature along the LOS are consistent with radiative transfer calculations, although quantitative differences persist, presumably related to the simplistic description of the coma of 67P that was used for these calculations.

Three spectra (S2, S3, and S4), corresponding to the average of data cubes fulfilling low-opacity conditions $N_{{\text{H}}_2\text{O}}<5\times {10}^{19}{\text{m}}^{\text{-2}}$, were used to infer the effective water OPR. In addition, we combined numerous spectra with high S/N (dataset H1), to measure the OPR from lines of the v₁ and hot-bands. According to our radiative transfer calculations, the range of OPR values derived from these datasets (2.26 to 2.79) corresponds to small to moderate optical depth effects. The H₂O OPR in comet 67P was finally estimated to be 2.94 ± 0.06, which is consistent with the statistical value of 3.

The H₂O OPR has been measured in numerous comets (see Faggi et al. 2018; Faure et al. 2019, for the most recent compilations). For several comets, values are near 2.5, corresponding to a spin temperature of ~30K. However, in other comets, estimates are consistent with the statistical value. Faure et al. (2019) found that whereas the mean value in a sample of 19 comets is 2.62 ± 0.03, the median value is 2.86, which is not far from the statistical value. Unequilibrated and equilibrated values are found for both Oort cloud and Jupiter-family comets. Most OPR determinations are from ground-based observations of ro-vibrational lines from H₂O hot bands near 2.9 μm. These determinations should not be significantly affected by optical depth effects, but a quantitative investigation remains to be performed. Interestingly, spin temperatures measured for NH₃ are also clustered near 30 K, and a good correlation between spin temperatures of H₂O and NH₃ is observed in the range of spin temperatures from 24 K to >40 K (Shinnaka et al. 2016).

Whereas in the past, OPRs in comets were believed to be of cosmogonic significance, that is, related to the formation history of cometary ices (e.g. Bockelée-Morvan et al. 2004), their real meaning is now unclear. Sliter et al. (2011) measured the OPR of H₂O monomers thermally desorbed from ice. They isolated H₂O in an Ar matrix at 4 K and left the condensed sample for about a day to obtain almost purely para isomers of H₂O. Subsequently, the Ar matrix was sublimated by rapid heating, and the IR spectrum of the water vapour sublimated from the ice at 250 K was recorded. It was found that the relative intensities of the rovibrational lines were equivalent to those of a H₂O vapour at room temperature. Thus, it was concluded that the conversion from para to ortho in ice is considerably fast. Hama et al. (2016, 2018) prepared H₂O ice at 10 K in two different ways: vapour deposition of H₂O and hydrogenation of O₂. In both cases, the OPR of H₂O photo or thermally desorbed from ice at 10 K and 150 K, respectively, was found to be the high-temperature-limit value of 3. Theoretical studies suggest that nuclear-spin conversion results from intermolecular proton-proton magnetic dipolar interactions and occurs on very short time scales (~10⁻⁵−10⁻⁴ s; Buntkowsky et al. 2008) in ice. In other words, the OPR of H₂O desorbed from cometary nuclei should have the statistical value of 3, regardless of the past formation process of the ice.

As discussed by Shinnaka et al. (2016), some processes in the coma might change the OPR after the molecules sublimated from the cometary surface. Nuclear-spin conversions by radiative transitions are strictly forbidden, but might occur by interactions with protonated ions such as H₃O⁺, water clusters, ice grains, and paramagnetic materials (e.g. O₂ molecules and some minerals in dust particles) in the near nucleus colli-sional region (Irvine et al. 2000; Hama et al. 2013). Bonev et al. (2007, 2008) measured the water OPR as a function of come-tocentric distance in comets 73P-B/Schwassmann-Wachmann 3 and C/2004 Q2 (Machholz) and did not observe any change in the ranges 5–30 km (73P) and 50–800 km (C/2004 Q2). For these two comets, the OPR is close to the statistical value. Most of the VIRTIS-H limb data have been acquired within ρ =10 km from comet centre. OPR variations that could be due to nucleus-spin conversion, if any, cannot be brought out from these close-nucleus data, which are affected by optical depth effects. However, the sensitive measurement obtained at ρ = 300 km near perihelion (dataset S2, at r_h = 1.31 au) is consistent with the other low-opacity data acquired much closer to the nucleus (Table 1), and does not support significant ortho-to-para conversion in the inner coma of comet 67P. Some comets present unequilibrated OPR values at close cometocentric distances. For example, a value of 2.59 ± 0.13 was measured for comet 103P/Hartley 2 from near-IR data sampling projected distances < 66 km from the nucleus (Bonev et al. 2013). Hence, nuclear-spin conversion might occur in the very inner coma of some comets (and possibly on the surface by interaction with paramagnetic impurities) by a process that remains to be identified. This is not observed for comet 67P, however.

Acknowledgements

The authors would like to thank the following institutions and agencies, which supported this work: Italian Space Agency (ASI - Italy), Centre National d’Etudes Spatiales (CNES – France), Deutsches Zentrum für Luft- und Raumfahrt (DLR – Germany), National Aeronautic and Space Administration (NASA – USA). VIRTIS was built by a consortium from Italy, France and Germany, under the scientific responsibility of the Istituto di Astrofisica e Planetologia Spaziali of INAF, Rome (IT), which lead also the scientific operations. The VIRTIS instrument development for ESA has been funded and managed by ASI, with contributions from Observatoire de Meudon financed by CNES and from DLR. The instrument industrial prime contractor was former Officine Galileo, now Leonardo company in Campi Bisenzio, Florence, IT. The authors wish to thank the Rosetta Science Ground Segment and the Rosetta Mission Operations Centre for their fantastic support throughout the early phases of the mission. The VIRTIS calibrated data are available through the ESA’s Planetary Science Archive (PSA) Web site and NASA Planetary Data System (PDS). Y.-C. Cheng aknowledges funding from IRIS OCAV (Université PSL). The work of M. Roos was supported by the program DIM-ACAV of Région Ile de France. With fond memories of Angioletta Coradini, conceiver of the VIRTIS instrument, our leader and friend, and to Michel Combes, one of the first architects of the VIRTIS-H channel.

Appendix A Example of spectral fits to S1 cubes

Fig. A.1 Spectral fits to the MB (left) and HB (right) spectral region for individual data cubes from dataset S1. See the caption to Figs 4 and 5. The observation identification number is given above the plots. Observing and model output parameters are listed in Table B.1.

Fig. A.2 Spectral fits to the MB (left) and HB (right) spectral regions for individual data cubes from dataset S1. See the captions to Figs 4 and 5. The observation identification number (ObsID) is given above the plots. Observing and model output parameters are listed in Table B.1.

Appendix B Results from spectral fits of dataset S1

References

All Tables

All Figures

Fig. 1 Continuum-subtracted VIRTIS-H spectrum of dataset H1 (blue curve) superimposed on the synthetic model that best fits the HB (2.774–2.910 μm) spectral region (orange curve). Dataset H1 is the weighted average of 22 data cubes obtained near perihelion (see Sect. 3.5.3 and Table 1 for geometric parameters and fitted parameters). The spectral range used to fit the HB region is indicated with a dotted green line. The bottom panel shows a synthetic optically thin spectrum at Trot = 116 K, with different colours for the v1 (green), v3 (grey), and hot bands (red).

In the text

	Fig. 2 Synthetic optically thin fluorescence H2O spectrum at the effective spectral resolution (and channel spacing) of VIRTIS-H. The blue and orange lines show the ortho and para lines, respectively. The total spectrum is plotted in grey. The ortho and para lines at infinite spectral resolution are shown at the bottom. The ortho-to-para ratio is equal to 3 and Trot = 100 K.
In the text

Fig. 3 Illustration of data combination and background subtraction. The grey spectra show the 1728 acquisitions of cube T1-00400765976 acquired on 13 September 2015 from 11h:46m:19.3s to 15h:24m:49.8s UT. The orange spectrum is the averaged spectrum, and the blue spectrum is the average spectrum obtained after subtraction of the continuum background, following the method explained in the text. This cube, which shows extreme variations in dust continuum because it was acquired during an outburst (Bockelée-Morvan et al. 2017), was not considered for this paper.

In the text

Fig. 4 Continuum-subtracted VIRTIS-H spectra (blue) superimposed on the synthetic model that best fits the MB (2.590–2.760 μm) spectral region (orange). Top: H1 dataset. Bottom: H2 dataset. The standard deviation for each spectral channel is shown at the bottom of the plots (and vertically offset with respect to the measured radiances for better clarity). The MB spectral region is shown by a horizontal green line.

In the text

	Fig. 5 Continuum-subtracted VIRTIS-H spectra (blue) superimposed on the synthetic model that best fits the HB (2.774–2.910 μm) spectral region (orange). Top: H1 dataset. Bottom: H2 dataset. The HB spectral region is shown by a horizontal green line. See the caption of Fig. 4 for other details.
In the text

	Fig. 6 Continuum-subtracted VIRTIS-H spectra (blue) superimposed on the synthetic model that best fits the MB (2.590–2.760 μm) spectral region (orange). Results for the S2, S3 and S4 spectra. See the caption of Fig. 4 for other details.
In the text

Fig. 7 Water column density $\left(N_{{\text{H}}_2\text{O}}^{\text{HB}}\right)$, rotational temperature $\left(T_{\text{rot}}^{\text{HB}}\right)$, apparent OPR (from MB fitting), and attenuation factor (fatten) derived from the fitting of VIRTIS-H spectra by optically thin fluorescence spectra. The 102 spectra of dataset S1 are considered. Plots A–C show how these parameters vary with each other. The colour-coding is indicated at the top of each plot.

In the text

Fig. 8 Rotational temperatures derived from fitting the MB (2.59–2.76 μm) and HB (2.774–2.91 μm) spectral regions. (A) $T_{\text{rot}}^{\text{HB}}$ as a function of distance to nucleus centre ρ, with a colour-coding according to heliocentric distance; only LOSs with a position angle with respect to the comet-Sun line PALOS < 45° were considered; results from S3-S4 and H1-H2 spectra are included; the fitted curves are $T_{\text{rot}}^{\text{HB}}=63+182/\rho \left(\text{K}\right)$ (dashed red line for rh < 1.4 au), and $T_{\text{rot}}^{\text{HB}}=35+225/\rho \left(\text{K}\right)$ (dashed green line for 1.4 au < rh < 1.8 au). (B) Ratio $T_{\text{rot}}^{\text{MB}}/T_{\text{rot}}^{\text{HB}}$ as a function of the attenuation factor fatten with symbol colours according to $T_{\text{rot}}^{\text{HB}}$.

In the text

Fig. 9 Optically thick RT spectrum (H2O v3 band) superimposed on the fitted optically thin spectrum. The input parameters for the calculation of the RT spectrum are OPR = 3, $N_{{\text{H}}_2\text{O}}=1.3\times {10}^{21}{\text{m}}^{-2}$, an LOS mean temperature (so-called $T_{\text{LOS}}^{\text{mean}}$) equal to 130 K, PALOS = 0°, and a terminator orbit. The results of the fit are apparent OPR (i.e., OPRRT) of 1.95, attenuation factor $f_{\text{atten}}^{\text{RT}}=0.61$, and rotational temperature $T_{\text{rot}}^{\text{RT}}=157\text{K}$.

In the text

Fig. 10 Results from radiative transfer calculations. Apparent OPR as a function of LOS mean temperature $\left(T_{\text{LOS}}^{\text{mean}}\right)$ and column density (A, B), and as a function of the attenuation factor (C). Plots A and B refer to calculations with an effective OPR of 2.5 and 3, respectively, with a colour-coding according to the derived apparent OPR. For C plot, the colour-coding is according to the column density (in base-10 log scale).

In the text

Fig. 11 Direct comparison of optical depth effects derived from model and observations. The column density (blue dots) obtained by fitting the synthetic RT spectra $\left(f_{\text{atten}}^{\text{RT}}\times N_{{\text{H}}_2\text{O}}\right)$ is plotted as a function of the effective column density $N_{{\text{H}}_2\text{O}}$. The values plotted as red dots are obtained by considering that the radiative transfer model underestimates the opacity effects by a factor f = 2.3 (see text). The dashed line shows the two-line fit of $\left(N_{{\text{H}}_2\text{O}}^{\text{MB}},N_{{\text{H}}_2\text{O}}^{\text{HB}}\right)$ values from dataset S1 (see text). The plain line shows optically thin conditions.

In the text

Fig. 12 Opacity $\tau =-\ln \left(f_{\text{atten}}^{\text{RT}}\right)$ as a function of column density (dots). Here, we assume that the radiative transfer model underestimates the opacity effects, and a factor of f = 2.3 is applied (see text), so that the values correspond to the red dots plotted in Fig. 11. Colour-coding is according to the LOS mean temperature $T_{\text{LOS}}^{\text{mean}}$, and the colour scale is given at the right. A polynomial fitting of degree 2 is shown for various temperatures (35, 55, 70, 80, 100, 120, and 140 K). The dashed line corresponds to the two-line fit of $\left(N_{{\text{H}}_2\text{O}}^{\text{MB}},N_{{\text{H}}_2\text{O}}^{\text{HB}}\right)$ values from dataset S1 shown in Fig. 11.

In the text

Fig. 13 Apparent OPR as a function of the attenuation factor fatten. Results from datasets S1 and S2–S4 are shown as grey and black dots, respectively. Results from datasets H1–H2 are shown by black squares. Only S1 data with an S/N on IHB > 20 are considered. For S2–S4 and H1–H2, the error bars for fatten represent the range of values obtained for f= 2.3 and f= 3. The arrows drawn for H1–H2 indicate a possible underestimation of fatten for these datasets. The curves show the polynomial regressions that fit the OPRRT versus $f_{\text{atten}}^{\text{RT}}$ values shown in Fig. 10C: red and green represent OPR = 2.5 and 3, respectively.

In the text

	Fig. A.1 Spectral fits to the MB (left) and HB (right) spectral region for individual data cubes from dataset S1. See the caption to Figs 4 and 5. The observation identification number is given above the plots. Observing and model output parameters are listed in Table B.1.
In the text

	Fig. A.2 Spectral fits to the MB (left) and HB (right) spectral regions for individual data cubes from dataset S1. See the captions to Figs 4 and 5. The observation identification number (ObsID) is given above the plots. Observing and model output parameters are listed in Table B.1.
In the text