© The Authors 2026

1 Introduction

Understanding how the Earth obtained its water and whether water-rich earthlike planets are common in the Universe is a fundamental cosmic origins question. Water forms efficiently in the cold interstellar medium via low-temperature ion-neutral chemistry or grain surface chemistry (van Dishoeck et al. 2013). Gaseous water and water-ice covered dust grains, either pristine or partially reprocessed, make their way through the subsequent stages of star formation, including dense cores and protostellar disks, before being incorporated into planetesimals in forming planetary systems (van Dishoeck et al. 2021). Understanding whether water and organics available for delivery to nascent planets in habitable zones around stars are supplied from their initial interstellar reservoir through a process that operates universally in forming planetary systems is a subject of active research. The possible links between the composition of interstellar and Solar System materials are of great interest for understanding the Solar System beginnings (Bockelée-Morvan et al. 2000; Drozdovskaya et al. 2019).

An investigation of the “water trail” was one of the main science themes of the Herschel Space Observatory (van Dishoeck et al. 2021). Herschel observations demonstrated that gas-phase water abundance is universally low, often orders of magnitude below the canonical value of 4 × 10⁻⁴ with respect to the molecular hydrogen, expected if all volatile oxygen is locked in water. Water vapor emission in molecular clouds was also found to be weak and difficult to detect away from embedded protostars. Melnick et al. (2011, 2020) concluded that water vapor is confined primarily within a few magnitudes of dense cloud surfaces. An exception are photon dominated regions (PDRs), surface layers of molecular clouds exposed to enhanced UV radiation.

Putaud et al. (2019) analyzed velocity-resolved observations of multiple water lines in the Orion Bar, one of the best studied PDRs, where warm chemistry with water ice desorption dominates. Using the Meudon PDR code (Le Petit et al. 2006), they derived the water abundance and ortho-para ratio in this region. They concluded that gas-phase water arises from a region deep into the cloud, corresponding to a visual extinction of A_V ~ 9. The ${\mathrm{H}}_2^{16}\mathrm{O}$ fractional abundance in this region is ~2 × 10⁻⁷, and the total water column density is (1.4 ± 0.8) × 10¹⁵ cm⁻². Contrary to earlier suggestions, a line-of-sight averaged ortho-para ratio was found to be consistent with a nuclear spin isomer repartition at the temperature of the water-emitting gas, 36 ± 2 K.

Pilleri et al. (2012) analyzed observations of water vapor in Mon R2, where the associated strongly irradiated PDR could be spatially resolved with Herschel. They derived a low mean abundance of o-H₂O of ~10⁻⁸ relative to H₂, and a higher abundance of ~10⁻⁷ in the high-velocity wings detected toward the H II region.

The Orion Bar is a dense PDR characterized by a very high flux of UV photons with G_o = (3–5) × 10⁴ in Habing units (Peeters et al. 2024) and a high thermal pressure of ~(1–3) × 10⁸ K cm⁻³ (Joblin et al. 2018; Putaud et al. 2019). Mon R2 is a relatively distant (830 pc) high-mass star forming region with a complex morphology and G_o = 5 × 10⁵ in its PDR. In the present work we analyze observations of water and other molecular and atomic tracers toward the Horsehead nebula, a region characterized by a simple geometry and much lower UV illumination with G₀ ~ 100 (Santa-Maria et al. 2023), which is more representative of the average UV-illuminated neutral gas in the Milky Way (e.g., Cubick et al. 2008). The archival data used in the analysis are described in Section 2, the results are presented in Section 3, and the PDR models used to interpret the observations are discussed in Section 4. Section 6 gives a summary and conclusions.

Fig. 1 (Left panel) SPIRE 350 μm image of the Horsehead PDR. The white circle shows the FWHM SPIRE beam (25.2″). Black circles show the locations at which 557 GHz o-H2O spectra were obtained, with the circle size corresponding to the FWHM HIFI beam (38.1″). (Right panel) Spectra of the 557 GHz o-H2O and CO (2–1) lines across the PDR (black and blue curves, respectively), labeled by the right ascension offsets with respect to the reference position. The CO 2–1 spectra have been scaled down by a factor of 200. The bottom-right panel shows spectra of the ground state o– and p–H2O lines at the central position.

2 Archival data

The 557 GHz 1₁₀–1₀₁ ground-state ortho-water line was observed in the frequency-switched (FSW) mode at six positions along a right ascension strip across the Horsehead PDR using the HIFI instrument (de Graauw et al. 2010) aboard the Herschel Space Observatory (Pilbratt et al. 2010) as a part of the open time program OT_dteyssie_2 (Herschel OBSID 1342203151). Spectra were taken at right ascension offsets from −40 to +60″, with 20" spacing, with respect to the center position α_J2000 = 05^h40^m53.7^s, δ_J2000 = −2°28′00″. In addition, FSW spectra of the same ortho-water line and the 1113 GHz 1₁₁–0₀₀ ground-state para-water line were observed at the center position as a part of guaranteed time program KPGT_vossenko_1 (OBSIDs 1342203191 and 1342217721). The on-source integration time is ~30 min per point for the 557 GHz line and 63 min for the 1113 Ghz line.

In our analysis, we used spectra obtained with the HIFI wide-band spectrometer (WBS) as a backend, which provided 1.1 MHz spectral resolution for H and V instrumental polarizations. The two polarizations were averaged together with equal weighting. The data downloaded from the ESA Herschel Science Archive¹ were reduced with version 14.1.0 of the HIFI pipeline. The latest values of the HIFI beam efficiencies of 0.635 at 557 GHz and 0.59 at 1113 GHz were applied to convert the spectra to the main beam brightness temperature scale.

The total observing time for the H₂O 1₁₀–1₀₁ strip varied between 20.2 and 33.1 min per position, resulting in an rms noise level between 3.0 and 4.6 mK in a single spectrometer channel. The integration times for the additional spectra toward the center positions were 27.8 and 63.1 min, for ortho- and para-water lines, respectively. However, these spectra are affected by baseline ripples and have relatively high rms noise levels of 9.3 and 18.8 mK, respectively.

In addition to the HIFI water spectra described above, we used in the analysis ground-based observations of the Horsehead nebula from Philipp et al. (2006), which include ~15″ resolution single-dish images of the [C I] 609 μm and CO (4–3) lines carried out using the CHAMP array at the Caltech Submillimeter Observatory, as well as observations of the (2–1) transitions of the CO isotopologs carried out with the IRAM 30 m telescope. We also used archival observations of [C II] downloaded from the SOFIA Science Archive at IRSA/IPAC².

3 Results

Figure 1 shows the overall morphology of the Horsehead PDR as traced by the 350 μm dust continuum emission observed with the SPIRE instrument on Herschel. The center position for the HIFI water observations is marked with a black asterisk, and the black circles show locations of the HIFI beams at the six positions along the strip. The white circle shows the full width at half maximum (FWHM) SPIRE beam at 350 μm.

The velocity-resolved HIFI water spectra of the 557 GHz along the east-west strip across the PDR are shown as black histograms in the top six panels on the right. The panels are labeled with the corresponding right ascension offset with respect to the reference position. CO (2–1) spectra from Philipp et al. (2006) at the same position are shown in blue. The water spectra are generally broader than CO spectra, indicating high line opacity. In particular, the water spectrum at the center (Δα = 0) position shows a self-absorption line profile. The bottom-right panel shows the additional spectra of the ortho- and para-water lines at the center position.

Figure 2 shows integrated line intensities of the various atomic and molecular tracers as a function of right ascension offset with respect to the reference position at a zero declination offset. The upper panel shows a clear stratification between the ionized and neutral carbon, with water emission peaking behind both tracers. The lower panels shows the distribution of water emission as compared to the three CO isotopologs. The water emission peaks at a similar location C¹⁸O, but extends further out into the ionized region.

To determine whether the observed stratification may be due to the differences in the angular resolution of the data, we convolved all images to the lowest spatial resolution of the HIFI o-H₂O data. The resulting strips are shown in Figure 3. The same stratification can still be seen, with water emission extending further out into the ionized region compared to other tracers. Integrated water line intensities between 7 and 14 km s⁻¹ at different positions are listed in Table 1 along with the corresponding uncertainties. Maximum line intensities of the various atomic and molecular tracers considered here in the HIFI beam are listed in Table 2.

Fig. 2 Integrated line intensities of atomic and molecular tracers computed over a 9.2–11.8 km s−1 velocity range across the PDR as a function of right ascension offset at Δδ = 0. (Top) [C II], [C I], and H2O (black, blue, and magenta, respectively). (Bottom) CO (2–1), C18O (2–1), CO (4–3), and H2O (black, blue, red, and magenta, respectively). Line intensities have been normalized to their maxima toward the PDR, with the corresponding scaling factors as labeled. The data are shown at the native resolution of the images, with the corresponding FWHM beam sizes marked for each tracer. The size of the squares corresponds approximately to ±1σ H2O observational uncertainties.

4 Pseudo 2-D PDR modeling of the Horsehead nebula

4.1 Previous approaches to model the Horsehead geometry with 1D plane-parallel codes

The Horsehead nebula has long served as a benchmark for the study of PDRs, owing to its well-defined, nearly edge-on geometry and wealth of multiwavelength observational data (e.g., Pety et al. 2005). Despite its apparent spatial simplicity, modeling the Horsehead poses specific challenges. Most PDR models, including those using the Meudon PDR code, adopt a 1D plane-parallel slab geometry in which the radiation field and physicochemical gradients vary the depth into the cloud, as the UV photons are gradually attenuated by dust. However, in the case of the Horsehead, the PDR is viewed nearly edge-on, meaning that the observer’s line of sight is perpendicular to the primary gradient direction in the model. As a result, direct comparisons between model predictions (computed along the depth axis) and observational quantities (integrated along the line of sight) require nontrivial adaptations, typically involving assumptions about the effective length of the PDR along the line of sight or geometrical corrections.

Habart et al. (2005), using SOFI observations of the H₂ 1–0 S(1) line at the NTT, modeled the Horsehead nebula as an edge-on, semi-infinite plane-parallel slab illuminated by an external radiation field and exhibiting a steep density gradient. This structure corresponded to an isobaric PDR with a thermal pressure of 4 × 10⁶ K cm⁻³. The computed emissivity was then scaled by the PDR length along the line of sight, assumed to be at most equal to the filament’s extent in the plane of the sky, i.e., approximately 0.1 pc.

Goicoechea et al. (2006) first used the output of Meudon PDR models as input for non-local thermodynamic equilibrium (non-LTE) radiative transfer models of CS and C¹⁸O emission lines, adapted to the edge-on geometry of the PDR. Similarly, Pabst et al. (2017) used the predictions of a face-on PDR model to calculate the expected [C II] 158 μm line emission. Both approaches compute the line-of-sight radiative transfer assuming 1D plane-parallel slabs.

Guzmán et al. (2011) also employed Meudon PDR models with a steep density gradient to interpret H₂CO transition lines observed with the IRAM-30m telescope, which were converted into chemical abundances. These derived abundances were compared to model predictions with and without surface chemistry, at two positions in the cloud referred to as the “core” and the “PDR”. A similar methodology (direct comparison between observed abundances and Meudon PDR model outputs) was applied in the interpretation of IRAM-30m observations of CF⁺, l-C₃H⁺, CH₃OH, and HCO in Pety et al. (2012); Guzmán et al. (2012, 2013).

Using sub-arcsecond angular resolution ALMA observations of the CO J = 3–2 line, Hernández-Vera et al. (2023) constrained the thermal pressure in the Horsehead nebula to the range (3.7–9.2) × 10⁶ K cm⁻³. This constraint was derived by comparing the observed spatial offsets between the C/CO and H/H₂ transitions with respect to the H/H⁺ transition, to those predicted by 1D Meudon PDR models. Only models within the specified thermal pressure range produced transition separations of the same order of magnitude as observed.

Fig. 3 Same as Figure 2 except that all tracers are convolved to the 38.1″ resolution of the HIFI H2O spectra.

4.2 The PDR wrapper approach

Here, we also use Version 7 of the Meudon PDR code³, but to compare the outputs of the code with spatially resolved observations of edge-on structures such as the Horsehead nebula, we developed a dedicated post-processing tool: the PDR wrapper. This tool maps the 1D Meudon PDR model onto a 2D plane defined by the lines of sight and the direction of radiation propagation, under the assumption of a simplified, shell-like geometry characterized by a fixed radius of curvature.

Figure 4 depicts a typical 1D PDR model geometry as used by the Meudon PDR code, with a two-component radiation field composed of a perpendicular beamed stellar radiation field from a nearby star and an isotropic mean interstellar radiation field (ISRF). We also see the chemical stratification, with the ionization front denoting the beginning of the PDR, followed by the dissociation front (H/H₂ transition) separating the atomic and molecular regions, and the C⁺/C/CO transition taking place deeper in the cloud, followed by the appearance of H₂O molecules.

An illustration of the wrapper’s conceptual geometry is provided in Figure 5. The left panel shows a real image of the Horsehead nebula, while the middle and right panels offer a rotated side-view artistic rendition of the cloud. In this representation, the stellar radiation field is incident from the top, while the lines of sight toward the observer extend horizontally to the right. The total set of lines of sight corresponds to a vertical column of pixels in the original image, effectively forming a cut across the nebula.

The shell-like structure superimposed in the right panel models the curved PDR surface, along which the output of a single 1D Meudon PDR model is mapped. This 1D structure is effectively “swept” along the arc, assuming a constant radius of curvature. Examples of such 1D models are drawn as red vertical lines on the figure.

A run of the Meudon PDR code produces a wide range of outputs. Among these, we make use of the density profiles of the chemical species of interest, as well as the profiles of their populations in individual quantum levels. Note that by using perpendicular illumination in the PDR code to simulate irradiation by the illuminating star, each point on the surface of the resulting sphere is illuminated perpendicularly. This approach provides a more realistic representation of a stellar-illuminated surface than would do a 1D spherical PDR model, in which the illumination is assumed to be isotropic.

Once the level populations were known throughout this 2D structure, we computed the line intensities as seen by an observer along many perpendicular lines of sight. This allowed us to reconstruct the spatial variation of the line intensities from the ionization front to the core of the cloud. The line intensities were determined by solving the radiative transfer equation along each line of sight, using the level populations interpolated from the Meudon PDR model, integrating from the far side of the cloud to the edge of the PDR on the observer’s side (see Appendix A for more details).

Lines of sight shown as solid arrows in Figure 5 intersect only physically meaningful regions of the cloud where the PDR model predicts densities and level populations. In contrast, dashed lines cross into central regions not covered by the model and are excluded from further analysis.

The wrapper reconstructs the spatial emission profiles as they would appear when observing a locally curved PDR edge-on, accounting for both curvature effects and radiative transfer along the lines of sight. This method enables the generation of synthetic spatial emission profiles while avoiding the computational cost of full 2D radiative transfer, which can be prohibitive for large parameter-space explorations while maintaining a very high level of process modeling.

Examples illustrating the wrapper geometry are shown in Figure 6, where a single Meudon PDR model (P_th = 4 × 10⁶ K cm⁻³, G₀ = 100 Habing units, Habing 1968) is projected onto shells with three different radii of curvature. The RGB composite maps simultaneously represent the abundances of the three main carbon-bearing species: red for ionized carbon (C⁺), green for atomic carbon (C), and blue for carbon monoxide (CO). This color-coding highlights the chemical stratification from the ionized outer layers to the molecular interior, with mixed regions appearing as intermediate colors such as yellow, cyan, or white. The observer is located at the bottom of each map, and the stellar radiation field enters from the right. These examples visualize the transition from ionized to atomic to molecular carbon, as well as the local spherical curvature along the line of sight. In this representation, radiative transfer is computed from top to bottom.

Fig. 4 Illustration of the geometry of 1D PDR models made with the Meudon PDR code. We show the two-component radiation field composed of a beamed part representing the stellar radiation field and an isotropic part representing the ISRF. This configuration yields the typical chemical stratification, with the ionization front denoting the beginning of the PDR, followed by the dissociation front, or H/H2 transition, and the C+/C/CO transition deeper in the cloud. H2O appears even deeper.

Fig. 5 Illustration of the geometry of the PDR wrapper. (Left panel) Image of the Horsehead nebula as seen by the Euclid space telescope (Credits: ESA/Euclid/Euclid Consortium/NASA, image processing by J.-C. Cuillandre (CEA Paris-Saclay), G. Anselmi). Notice the reference frame, with the x-axis corresponding to the direction of the lines of sight. We rotate the nebula so that the x–axis points to the right. (Middle panel) Artistic interpretation of the rotation of the Horsehead nebula, generated using OpenAI’s DALL·E model based on the first panel image. (Right panel) Geometry of the wrapping 2D model, with the star radiation field coming from the top, and examples of Meudon PDR model output overplotted as vertical red lines. Each of these lines correspond to a PDR model as represented on Fig. 4. Notice that we use the same model for each position along the cut of the cloud. The multiple lines of sight correspond to different pixels along a vertical 1D cut in the Horsehead nebula. We emphasize that only the lines of sight shown as solid lines are correctly accounted for by the model, as the central section of the nebula outside of the crescent-shaped area is not included in the computations. The geometry is assumed to be circular with a radius of curvature that remains to be determined.

Fig. 6 Illustration of the wrapper geometry, showing RGB composite maps of carbon-bearing species for different radii of curvature. The color encodes the local abundances as follows: red for ionized carbon (C+), green for neutral atomic carbon (C), and blue for carbon monoxide (CO). Areas where these species coexist appear as color blends (e.g., cyan for C and CO, yellow for C+ and C, white for equal contributions). The left panel shows a curvature radius of RC = 0.05 pc, the middle to RC = 0.2 pc, and the right to RC = 0.5 pc. These maps highlight the local spherical curvature in the plane of the sky, with the observer located at the bottom of each panel and the stellar radiation field entering from the right. Radiative transfer is computed along vertical lines from top to bottom in this configuration. The abundances originate from a single Meudon PDR model (Pth = 4 × 106 K cm−3, G0 = 100). The usual edge-on representation of abundances of this model can be found in Figure 10.

4.3 Model grid and physical ingredients

As discussed in Section 4.5, reproducing the H₂O emission line proved challenging. For this reason, we separate the analysis into two parts: the reproduction of the C⁺/C/CO transitions, based on both line intensities and spatial profiles, is presented here, while the modeling of the H₂O line and its spatial distribution is treated independently.

For the first part of the analysis, we constructed a grid of isobaric Meudon PDR models coupled with the wrapper tool. The explored parameter space includes:

• thermal pressures P_th = 0.5, 1, 2, 3, 4, 5, 6, 8, and 10 × 10⁶ K cm⁻³; and

• curvature radii R_C = 0.02, 0.05, 0.1, 0.2, 0.5, and 1.0 pc.

A standard value for the interstellar radiation field of G₀ = 100 Habing units is used for the models, as the line intensities of [C II] 158 μm, [C I] 609 μm, and low-J CO transitions depend only weakly on G₀ in low-excitation PDRs such as the Horsehead nebula (see, e.g., Kaufman et al. 1999). This choice of G₀ is motivated by multiple studies of the Horsehead nebula (e.g., Abergel et al. 2003; Habart et al. 2005; Goicoechea et al. 2009; Santa-Maria et al. 2023; Zannese et al. 2025) and the spectral type of the illuminating star, σ Ori, O 9.5 V at a projected distance of 3.5 pc. We further discuss the weak G₀ dependence in Appendix C.

The models incorporate an exact radiative transfer treatment for H₂ lines. This means that H₂ self-shielding is calculated precisely, properly accounting for the overlap of UV pumping lines from both H and H₂. This formalism also ensures that the model accurately considers the shielding of the CO UV lines by the H₂ Lyman and Werner transitions.

The chemical network considers 220 species and more than 3300 reactions including a surface chemistry network for the chemistry of carbon and oxygen and the formation of mantles of H₂O on the grains. The radiation field is modeled using a synthetic spectra from the POLLUX database (Palacios et al. 2010) to simulate the σ-Ori stellar spectrum. The selected parameters are T_eff = 33 000 K, log (g) = 4.26 and R = 5.26 R_⊙. Turbulent broadening is included with a velocity dispersion σ_turb = 0.38 km s⁻¹, following Hernández-Vera et al. (2023), derived from ALMA observations of HCO⁺ in the Horsehead nebula. The impact of this choice on the predicted line intensities is discussed in Sect. 4.5.

4.4 Model fitting procedure and main results

The synthetic emission profiles produced by the wrapper were compared to observational data using a χ² minimization method implemented in the pyismtools Python suite. The fitting simultaneously considered both the maximum intensities and, for the first time, the spatial separation between the emission peaks along the PDR. For the minimization, we did not use the raw output of the wrapping procedure, but rather the profiles convolved to match the native angular resolution of each tracer image. The minimum χ² value is obtained for P_th = 4.1 × 10⁶ K cm⁻³ and R_C = 0.057 pc, being consistent with previously published, with R_C being around half the typical 0.1 pc width usually assumed in plane-parallel models of the Horsehead PDR (Figure 7).

A comparison between the best model and the observations is shown in Fig. 8. The bottom panel displays the maximum intensities of the emission lines, demonstrating excellent agreement. Model predictions are plotted as empty circles with 40% error bars, accounting for geometrical uncertainties and the complexity of the parameter space. In the top panels, we compare the full spatial emission profiles from our model to the observed profiles, even though the minimization relied solely on maximum intensities and separations between various tracers. The left panel shows the raw model profiles produced by the Wrapper, while the right panel presents the same profiles convolved to the native angular resolution of the respective tracers, as used in the minimization. As stated before, we emphasize that our modeling focuses solely on reproducing the initial rise in intensity for each emission line, since the geometrical model is physically consistent only within the computed width of the PDR.

Notably, even though only maximum positions were used for fitting, the complete spatial profiles are well reproduced. An exception is observed for the CO lines and their isotopologs, for which the observed rise in intensity near the illuminated edge is smoother than predicted. This discrepancy remains unexplained and may be due to convolution of multiple CO fronts (eg., Hernández-Vera et al. 2023).

Fig. 7 A map of χ2 values obtained from the minimization procedure using maximum intensities and spatial separations. The best-fit model, corresponding to Pth = 4.1 × 106 K cm−3, and RC = 0.057 pc, is denoted by the orange circle on the map. The constraint on RC remains weak, with the χ2 surface relatively flat in a wide valley.

4.5 H₂O emission

Figures 8 and 9 show the spatial emission profiles (after the wrapping procedure, and convolution to the native angular resolution of each tracer) and the maximum intensities including the H₂O 1₁₀–1₀₁ line, for the model: P_th = 4 × 10⁶ K cm⁻³, G₀ = 100, and R_C = 0.05 pc. For water, the model predicts an intensity about 20 times higher than the observed value, and a spatial distribution of the emission much more recessed into the cloud. Models with lower G₀ values provide a somewhat better fit to the observations. However, the model water line intensity is stil a factor of seven higher than the observations (see Fig. C.2 in Appendix C).

There are several possible explanations for these discrepancies. First, H₂O excitation is very sensitive to radiative transfer effects (Poelman & Spaans 2006; Gonzalez Garcia et al. 2008), such as nonlocal pumping by warm dust emission. Also, the line is strongly optically thick, with optical depths of the order of 10³, so the emergent intensity is determined by the local excitation temperature at the τ ~ 1 surface. The location of this surface is highly sensitive to both the H₂O abundance profile and the line width, which is controlled by turbulent broadening. A smaller turbulent velocity leads to a narrower line, increasing the opacity, and thus shifting the τ ~ 1 surface outward, closer to the cloud edge where the excitation temperature is lower, thereby decreasing the emergent intensity.

To test the sensitivity of the results to turbulent velocity dispersion, we also considered the smaller value of σ_turb = 0.225 km s⁻¹ derived from the analysis of Gerin et al. (2009), who derived this value for the Horsehead from HCO observations combining Plateau de Bure Interferometer (PdBI) and IRAM-30m data. In this case the predicted intensity decreases by about 24%, from 7.58 × 10⁻⁷ erg cm⁻² s⁻¹ sr⁻¹ to 5.75 × 10⁻⁷ erg cm⁻² s⁻¹ sr⁻¹. Both predictions remain far above the observed value of 3.67 × 10⁻⁸ erg cm⁻² s⁻¹ sr⁻¹. For the other tracers, we find that decreasing the turbulent velocity has almost no effect on [C II] and [C I], but it reduces the CO line intensities by ~20%, thereby slightly improving the agreement with the observations.

In addition, the intensity of optically thick lines varies non linearly with the abundance profile. A small change in the predictions of the density of H₂O can result in significant variations in line intensities. We checked that our model is perfectly in line with the classical PDR model of the chemistry of H₂O (Hollenbach et al. 2009), except that surface chemical data (e.g., binding energies) have been updated and the modeling of some processes has been modernized (Cuppen et al. 2017). Also, we know that the model reproduces H₂O line intensities in more intense PDRs such as the Orion Bar (Putaud et al. 2019). This may suggest that particular processes may be at work in the Horsehead, possibly involving grains with unusual properties, whether in composition or size distribution.

The detailed analysis of dust and H₂ emission recently observed by JWST in the Horsehead nebula could help clarify this issue (Abergel et al. 2024). These findings also highlight the need to refine the modeling of chemical processes on interstellar grains, particularly by acknowledging that the physical conditions of these grains likely differ from those under which binding energies and process efficiencies, such as photo-desorption yields, are measured in laboratory experiments (Hacquard et al. 2024). For instance, given the diversity of interstellar grain types, it would be more realistic to adopt distributions of binding energies for adsorbed species instead of single fixed values. However, this would require the development of new formalisms to accurately simulate surface chemistry. Last but not least, another factor to consider is the geometry of the emission and nonlocal radiative transfer effects: the discrepancy between models and observations may well arise from scattering of the radiation emitted by water molecules, as discussed in the following section.

Fig. 8 (Top) Spatial profiles of normalized emission line intensities observed in the Horsehead nebula, compared with predictions from a Meudon PDR model with Pth = 4 × 106 K cm−3 and G0 = 100, wrapped in a locally spherical cloud with a curvature radius RC = 0.05 pc. These conditions correspond to commonly adopted literature values. The left panels show the raw profiles from the wrapping procedure, while the right panels present the same profiles convolved to the native angular resolution of each tracer. In those figures, the radiation field comes from the right. (Bottom) Observed maximum intensities of various tracers (colored crosses) compared to model predictions (circles with error bars). All tracers are well reproduced, except for H2O, where the model overestimates the observed maximum intensity by an order of magnitude.

Fig. 9 Spatial profiles of normalized emission line intensities observed in the Horsehead nebula, compared to predictions from Meudon PDR model with Pth = 4 × 106 K cm−3 and G0 = 100, wrapped in a locally spherical cloud with curvature radius RC = 0.05 pc. The profiles are convolved to the native angular resolution of the tracer maps. On this figure, the radiation field comes from the right. We observed a significant shift between the observed H2O profile and the model, the latter being much deeper in the cloud.

5 Scattered o-H₂O 1₁₀–1₀₁ 557 GHz line emission by low-density warm envelopes

In this section, we describe a plausible scenario explaining why the o-H₂O 1₁₀–1₀₁ 557 GHz line emission is more extended toward the rim of the PDR than the ¹²CO (2–1) emission, and why it nearly coincides with the [CII] 158 μm emission when convolved to the same low angular resolution (~38″), which traces lower-density gas at the rim of the PDR (Pabst et al. 2017; Bally et al. 2018). Firstly, we note that a steep gas density gradient exists in the Horsehead, with densities varying from the PDR/H II interface (n_H ≈ a few 10³ cm⁻³, T ≳ 100 K) to the more shielded, cold molecular cores (n_H ≈ a few 10⁵ cm⁻³, T ≈ 20 K) (e.g., Habart et al. 2005; Pety et al. 2007; Gerin et al. 2009; Goicoechea et al. 2009; Hernández-Vera et al. 2023). The nearly edge-on orientation of the PDR allows these different environments to be partially spatially resolved.

Our PDR models (Sect. 4.5) predict that the peak of the water vapor abundance, reaching 2.6 × 10⁻⁷, occurs at A_V ~ 3 mag, well beyond the CO/C transition layer of the PDR, as can be seen on Figure 10 showing the edge-on Meudon PDR model output (thus without the wrapping). As a reference, models from Putaud et al. (2019) yielded a comparable abundance of H₂O of 2 × 10⁻⁷, that occurred at A_V ~ 9 mag in the Orion Bar, due to the substantially stronger radiation field. Our water abundance peak roughly corresponds to the point where the rate of gas-phase O atoms sticking to grains equals the photo-desorption rate of water molecules from icy grain surfaces. At this depth, the gas temperature is approximately 40 K, and the H₂O peak lies deeper than the ¹²CO abundance peak. In this region of the PDR, the H₂O abundance profile closely resembles that of C¹⁸O. Hence, contrary to observations, one would expect the o-H₂O 557 GHz emission to be as extended as the C¹⁸O (2–1) emission. In contrast to C¹⁸O, however, H₂O starts to form close to the PDR surface, ahead of the CO/C transition, where the C⁺ abundance is much greater than that of CO, temperatures are high, and densities are lower than at the H₂O abundance peak. With this in mind, one can approximate the H₂O line radiative transfer problem in the Horsehead using a two-component model: a denser, moderately cold PDR component with a large H₂O column density (in the line of sight), surrounded by a lower-density, warm envelope with a lower H₂O column density. This envelope represents the bright [CII]158 μm but faint ¹²CO (2–1) emission edge of the PDR (e.g., Pabst et al. 2017; Hernández-Vera et al. 2023).

Owing to the very high critical density for collisional excitation, n_cr (1₁₀–1₀₁) ≃ 10⁸ cm⁻³, and the elevated 557 GHz line opacity for moderate H₂O column densities, the excitation of the 1₁₀–1₀₁ transition is very subthermal (n_H ≪ n_cr) (see, e.g., Flagey et al. 2013), which implies T_ex ≪ T, and is prone to nonlocal radiative effects such as line resonant scattering by low-density envelopes containing water vapor. In such envelopes, the local o-H₂O 557 GHz line emission is negligible because T_ex → T_cmb. However, when o-H₂O 557 GHz line photons emitted by the deeper PDR layers diffuse through the envelope, they are self-absorbed and re-emitted by water molecules, which increases T_ex in the low-density gas and scatters the emission to greater distances. This produces an impression of a more extended spatial distribution of H₂O emission. While this scattering process has been invoked to explain the extended distribution and anomalous J = 1–0 line intensity ratios of very polar molecules such as HCO⁺ or HCN in dark clouds (i.e., rotational transitions with a high n_cr; Langer et al. 1978; Walmsley et al. 1982; Cernicharo et al. 1984; Gonzalez-Alfonso & Cernicharo 1993; Goicoechea et al. 2022), all these studies required cold envelopes. The interesting property of H₂O is that the critical density of the ground-state 1₁₀–1₀₁ transition is at least 3 orders of magnitude higher than that of HCN and HCO⁺ J = 1–0. As we demonstrate below, this implies that the excitation temperature of the 1₁₀–1₀₁ transition remains very low even at high gas temperatures, provided that the gas density is significantly lower than the effective critical density of this transition. In this case, a warm/hot envelope containing moderately water vapor abundances will produce self-absorption and resonant scattering due to very low excitation temperature of the 1₁₀–1₀₁ in the envelope. This nonlocal radiative effect cannot be reproduced by standard (local) escape probability methods.

To investigate the impact of line scattering on the emergent H₂O line intensities, we use a nonlocal and non-LTE Monte Carlo radiative transfer code (see Goicoechea et al. 2006, 2022) and develop a “simple” two-component model of the water vapor line emission adapted to the physical conditions known to exist in the Horsehead PDR. We used the H₂O–H₂ inelastic collision rate coefficients of Daniel et al. (2011, 2012). Following a large list of observational studies of this source, we adopt n_H=6×10⁴ cm⁻³ and T =40 K for the dense PDR component, and n_H=6×10³ cm⁻³, T =60 K, and for the extended envelope near the PDR/HII region interface. We found a good match to the observed line intensities by adopting a (beam-averaged) line-of-sight o-H₂O column densities, N(o-H₂O) ≃ 4.3×10¹³ cm⁻² and N(o-H₂O) ≃ 1.7×10¹³ cm⁻² in these two components, respectively. These models include thermal, microturbulent, and line opacity broadening. To reproduce the observed line-widths, we adopt σ_turb ≃ 0.4 km s⁻¹ in the dense PDR and σ_turb ≃ 1 km s⁻¹ in the extended envelope. An increased velocity dispersion in the envelope diminishes the depth of the 557 GHz self-absorption line dip. Figure 11 shows the resulting line and radial T_ex (1₁₀–1₀₁) profiles. The blue curves represent a model of the dense PDR component-only, while the red curves represent a model of the low-density envelope-only. The black curves represent the two-component model.

Due to the low gas density, the envelope-only model shows highly subthermal (T_ex ≃ 3.2 K) line emission, despite the moderate gas temperatures. Thus, the predicted o-H₂O 557 GHz line intensities are very faint (a few millikelvins) and would remain undetectable by Herschel/HIFI. The dense-PDR-only alone, with and order of magnitude higher gas density and a lower temperature (closer to the energy of the upper level 1₁₀), leads to more efficient collisional excitation and thus increased excitation temperatures (T_ex ~ 4.5–5.5 K). In addition, the o-H₂O 557 GHz line from the dense component is optically thick, which produces line opacity broadening and line-trapping (i.e., an increase of T_ex at low radii despite the gas density is constant). Compared to the envelope-only model, the two-component model (black curves) has different excitation conditions in the envelope region (R > R_densePDR). The absorption and reemission of o-H₂O 557 GHz line photons has two consequences: (i) It reduces the apparent o-H₂O 557 GHz line intensity toward the dense component, and (ii) it increases T_ex (1₁₀–1₀₁) in the envelope to levels that produce detectable scattered o-H₂O 557 GHz line emission at R > R_densePDR. This is a purely radiative effect and leads to o-H₂O 557 GHz line emission beyond the dense PDR component. This scenario would imply that our detection of the o-H₂O 557 GHz emission at the rim of the PDR is scattered emission from deeper layers of the PDR.

Figure 12 shows the predicted radial o-H₂O 1₁₀–1₀₁ integrated line intensities, with the black curves representing the two-component model. This model results in detectable o-H₂O 557 GHz emission that is more extended than the dense PDR region. It also results in an decrease in the o-H₂O 557 GHz line intensity for line of sights toward the core (e.g., for R/R_densePDR = 0). This simple two-component model reasonably fits the observed o-H₂O 557 GHz line intensities across the PDR (magenta-filled circles in Fig. 12). We recall that a necessary condition for this radiative effect to occur is the presence of an extended, low-density envelope (lower excitation) containing moderate amounts of water vapor, as predicted by our PDR photochemical model. The size of the H₂O scattered emission region and the degree of line self-absorption ultimately depends in the exact physical conditions and water vapor abundance in the different components. We note that resonant scattering does not affect the excitation or emission of low-critical-density transitions (e.g., low-J CO, [C II], or [C I] lines). Indeed, most of the C¹⁸O 2–1 emission originates in the dense PDR interior, while any C¹⁸O present in the envelope will be excited collisionally (T_ex ≃ T_k) and not radiatively. As the line emission from the dense PDR is optically thin, this emission is not self-absorbed. For C⁺, the [C II] 158 μm line (also collisionally excited) arises mainly from the foreground envelope. These two cases, both with n_cr of a few 10³ cm⁻³, are presented in Appendix B.

In summary, our simple model captures the average physical conditions and H₂O column densities within the large Herschel/HIFI beam. Higher angular resolution observations of water vapor will be necessary to properly resolve the size of each component and to account for the gradients in physical conditions and water abundance.

Fig. 10 Abundances of key species with respect to hydrogen nuclei (top panel) and gas temperature (bottom panel) profiles of the Horsehead nebula model with Pth = 4 × 106 K cm−3, G0 = 100, computed with the Meudon PDR code (no wrapping or convolution procedure is applied to these profiles). The H/H2 transition can be seen to take place at around AV ~ 0.06 mag, the C+/C/CO transition is at AV ~ 1.7 mag, and the peak of H2O abundance is reached at an AV ~ 3 mag.

Fig. 11 Two-component, nonlocal radiative transfer models adapted to the physical conditions in the Horsehead. The upper panel shows the o-H2O 557 GHz synthetic line profiles (continuous curves) and observed line profiles (histograms) at the emission peak. The lower panel shows the predicted excitation temperatures of the 110–101 transition. The blue curves represent the model of the dense PDR component alone. The red curves represent the model of the low-density envelope alone, while the black curves represent the two-component model.

Fig. 12 Same as the lower panel of Fig. 11 but showing the predicted radial o-H2O 110–101 integrated line intensities. The magenta-filled circles represent the observed line intensities across the Horsehead PDR, with the peak position at Δα = +20″.

6 Summary and conclusions

We have analyzed archival Herschel observations of water emission toward the Horsehead PDR, along with supporting ground-based and airborne observations of CO isotopologs and fine structure lines of ionized and atomic carbon. The main results of this study can be summarized as follows:

• Water emission in the Horsehead nebula is very weak and, surprisingly, extends outward beyond other PDR tracers, such as CO isotopologs or [C I] 609 μm, and as far out as [C II] 158 μm;

• We developed a PDR wrapper to model line emission from 3D objects illuminated by a flux of UV photons. PDR modeling shows that the [C II], [C I], and CO isotopolog intensities and spatial profiles provide strong constraints on the thermal pressure (P_th = 4 × 10⁶ K cm⁻³), but not on G₀. A model with published source parameters provides a reasonable fit to the data, while underestimating the intensities of the [C II] and CO J = 4–3 lines by a factor of ~2.5. Spatial profiles are also well reproduced, except for CO isotopologs, where the increase on the illuminated side of the PDR is steeper than observed;

• Water vapor abundance in the PDR model reaches 3.6 × 10⁻⁷ at A_V ~ 3 mag. However, the ground state o-H₂O 557 GHz line is systematically overestimated by the models by at least a factor of 7 for any values of the model parameters. This line has a very high optical depth, ~10, and the emergent line intensity is sensitive to radiative transfer effects and the assumed turbulent line width. Nevertheless, the PDR modeling suggests that some ingredients describing the processes that lead to the formation and destruction of water ice on grain surfaces are missing in the model;

• We demonstrate that scattering of photons by water molecules in a low-density envelope surrounding the dense PDR provides a plausible explanation for the observed spatial extent of the water emission. While an isobaric Meudon PDR model naturally includes a low-density, warm envelope, the local escape and radiative transfer formalism in the code does not treat correctly the resonant scattering of water photons by the outer, lower-density gas layers.

The PDR wrapper developed here will be a powerful tool for analyzing future high angular resolution observations of molecular tracers in the Horsehead nebula obtained using ground-based millimeter-wave interferometers. Such observations will much better delineate the spatial stratification among various PDR tracers than the single-dish observaions presented here, providing more stringent constraints on the physical parameters.

Acknowledgements

This research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). D.C.L. acknowledge financial support from the National Aeronautics and Space Administration (NASA) Astrophysics Data Analysis Program (ADAP). HIFI was designed and built by a consortium of institutes and university departments from across Europe, Canada, and the United States (NASA) under the leadership of SRON, Netherlands Institute for Space Research, Groningen, The Netherlands, and with major contributions from Germany, France, and the US. This research has made use of the NASA/IPAC Infrared Science Archive, which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology. JRG thanks the Spanish MCINN for funding support under grant PID2023-146667NB-I00. This work was supported by the Thematic Action “Physique et Chimie du Milieu Interstellaire” (PCMI) of INSU Programme National “Astro”, with contributions from CNRS Physique & CNRS Chimie, CEA, and CNES. This research was achieved using the POLLUX database (pollux.oreme.org) operated at LUPM (Université de Montpellier – CNRS, France) with the support of the PNPS and INSU. V.M. and J.R.G. thank the Spanish MCINN for funding support under grant PID2023-146667NB-I00. We thank an anonymous referee for helpful comments.

References

Appendix A Presentation of the PDR wrapper

A.1 Aims of the code

The Meudon PDR code (Le Petit et al. 2006) models interstellar clouds as 1D, plane-parallel slabs, infinite and homogeneous along their surface. It computes the steady-state gas and chemistry structure along the normal to the surface, as a function of optical depth, and predicts line intensities emerging from the illuminated side. Level excitation is determined at each position in the cloud by solving the non-LTE problem, taking into account radiative and collisional transitions, as well as chemical formation and destruction in specific levels. The algorithm that computes the stationary level populations is intimately coupled to the solution of the radiative transfer equation, as described in Gonzalez Garcia et al. (2008). The algorithm considers non local pumping in the continuum (dust emission plus CMB). For example, the radiation field emitted by warm dust at the edge of the PDR can pump transitions deeper in the cloud where dust is cold. Concerning line emission, the algorithm uses a on the spot approximation meaning that the line can be re-absorbed, but only locally. For H₂O, the code considers inelastic collisions with H (Daniel et al. 2015), He (scaling of the collision rates with H), ortho and para H₂ (Daniel et al. 2011), e⁻ (Faure et al. 2004).

Edge-on clouds illuminated from the side reveal PDR stratification through spatially offset line emissions. Although their complex geometry requires simplifications, the Meudon code allows for a detailed treatment of microphysics, chemistry, and radiative transfer. Extending this approach to a full 2D model with the same level of refinement in the treatment of physical processes would, however, be computationally prohibitive.

When the PDR width is small compared to the cloud radius, radiative transfer can be simplified by mapping the 1D model onto a 2D spherical geometry. This is the basis of the PDR Wrapper, which reconstructs line emission by solving radiative transfer along each line of sight through a spherical cloud cross-section.

A.2 Geometry

We assumed a locally spherical curvature for the cloud surface, with center of curvature O and radius R_C. The PDR is modeled as the outermost shell of thickness d_PDR ≪ R_C.

The radiation field is assumed constant and perpendicular to the line of sight. The cloud has a curvature radius R_C and a constant horizontal thickness d_PDR. It is modeled as the shell between two arcs of radius R_C, with centers shifted horizontally by d_PDR (Fig. A.1). For a point on the outer arc making an angle α with the horizontal, the depth d from the illuminated edge is $$d=R_{\mathrm{C}}(\cos \alpha -1)+b,$$(A.1)

with sin α = s/R_C. The expression for s_max remains unchanged, while: $$s_{\mathrm{m}\mathrm{i}\mathrm{n}}(b)=\{\begin{array}{ll}0& \text{\hspace{0.17em}if\hspace{0.17em}}b<d_{\mathrm{P}\mathrm{D}\mathrm{R}},\\ \sqrt{R_{\mathrm{C}}^2-{(R_{\mathrm{C}}-b+d_{\mathrm{P}\mathrm{D}\mathrm{R}})}^2}& \text{\hspace{0.17em}otherwise\hspace{0.17em}}.\end{array}$$(A.2)

A.3 Radiative transfer

Each image pixel integrates emission and absorption along a single line of sight. Our model solves the radiative transfer equation along multiple, parallel, and independent lines of sight across the cloud.

Fig. A.1 Horizontal geometry used for the wrapping of Meudon PDR models. The beamed ISRF remains perpendicular to the line of sight.

The goal is to compute the intensity of a user-selected line at frequency ν_ul, where u and l are the upper and lower energy levels, respectively. Accounting for spontaneous emission (A_ul), stimulated emission (B_ul), and absorption (B_lu), the radiative transfer equation reads $$\frac{\mathrm{d}{I}_{\nu }}{\mathrm{d}s}=\frac{\varphi (\nu ,{\nu }_{ul},s)}{4\pi }h\nu [A_{ul}{n}_u(s)+(B_{ul}{n}_u(s)-B_{lu}{n}_l(s))I_{\nu }],$$(A.3)

where n_u(s) and n_l(s) are the level populations along the line of sight.

The normalized line profile ϕ is assumed Doppler-broadened due to thermal motions and turbulence, and is given by: $$\varphi (\nu ,{\nu }_{ul},s)=\frac{1}{{\sigma }_{{\nu }_{ul}}(s)\sqrt{2\pi }}\exp [-\frac{{(\nu -{\nu }_{ul})}^2}{2{\sigma }_{{\nu }_{ul}}^2(s)}],$$(A.4) $${\sigma }_{{\nu }_{ul}}=\frac{\sqrt{2}}{2}\frac{{\nu }_{ul}}{c}\sqrt{\frac{2k_{\mathrm{B}}T}{m}+{\sigma }_{\text{turb}}^2},$$(A.5)

where the normalization relation is $${\int }_0^{+\mathrm{\infty }}\varphi (\nu ,{\nu }_{ul},s)\mathrm{d}\nu =1.$$(A.6)

The Einstein relations (valid when working with I_ν) are: $$B_{ul}=\frac{c^2}{2h{\nu }_{ul}^3}{A}_{ul},B_{lu}=\frac{g_u}{g_l}{B}_{ul}.$$(A.7)

Rewriting the transfer equation as: $$\frac{\mathrm{d}{I}_{\nu }}{\mathrm{d}s}={ϵ}_{\nu }-{\alpha }_{\nu }{I}_{\nu },$$(A.8)

where we define the emissivity and absorption coefficients: $${ϵ}_{\nu }=\frac{\varphi (\nu ,{\nu }_{ul},s)}{4\pi }h\nu A_{ul}{n}_u(s),$$(A.9) $${\alpha }_{\nu }=\frac{\varphi (\nu ,{\nu }_{ul},s)}{4\pi }h\nu [B_{lu}{n}_l(s)-B_{ul}{n}_u(s)].$$(A.10)

The formal solution is $$I_{\nu }(s)={\int }_{s_{min}}^s{ϵ}_{\nu }\left(s^{\mathrm{\prime }}\right)e^{-{\tau }_{\nu }\left(s^{\mathrm{\prime }}\right)}\mathrm{d}{s}^{\mathrm{\prime }},{\tau }_{\nu }(s)={\int }_{s_{min}}^s{\alpha }_{\nu }\left(s^{\mathrm{\prime }}\right)\mathrm{d}{s}^{\mathrm{\prime }}.$$(A.11)

Finally, the total emergent intensity along the line of sight (impact parameter b) is obtained by frequency integration: $$I(b)={\int }_0^{+\mathrm{\infty }}{I}_{\nu }(b)\mathrm{d}\nu .$$(A.12)

A.4 Convolution with instrumental PSF

The PSF (point spread function) characterizes the response of an instrument to a point-like source. It defines the instrument angular resolution. To model the instrumental response to the spatial profiles, we convolve I(b) or N(b) with a Gaussian kernel: $$g(x)=\frac{1}{\sigma \sqrt{2\pi }}\exp (-\frac{x^2}{2{\sigma }^2}),$$(A.13)

where x is the distance to the center of the PSF, and σ the standard deviation of the PSF, linked to the full width half maximum (FWHM) of the instrument: $$\sigma =\frac{\mathrm{F}\mathrm{W}\mathrm{H}\mathrm{M}}{2\sqrt{2\ln 2}}.$$(A.14)

The convolution product is thus: $$(I\ast g)(b)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}I(b^{\mathrm{\prime }})g(b-b^{\mathrm{\prime }})\mathrm{d}{b}^{\mathrm{\prime }}.$$(A.15)

In practice, the integral is limited to the range of b plus or minus 3 times the standard deviation σ of the gaussian kernel. Outside of the range of b, the intensity is taken as a linear extrapolation of the profile.

Appendix B Two-component model of the C¹⁸O 2–1 and [C II] 158 μm line emission in the Horsehead

In this Appendix, we present the results of applying the simple two-component model to C¹⁸O and C⁺. These two species serve as examples of line emission arising predominantly from the dense PDR and from the extended PDR envelope, respectively. In contrast to H₂O, the C¹⁸O 2–1 and [C II] 158 μm lines are low-critical-density transitions, with n_cr of only a few ×10³ cm⁻³. This means that they are collisionally excited (T_ex ≃ T_k) in the envelope and therefore do not scatter line emission originating from the dense PDR, which we demonstrate here. These representative models are designed to reproduce the observed line peak intensities and their approximate spatial distribution qualitatively. Specifically, C¹⁸O 2–1 peaks toward the dense PDR, while [C II] 158 μm peaks further ahead. The physical conditions and line-widths are the same as those of the H₂O models (Sect. 5). Following the PDR model, the adopted C¹⁸O abundances relative to H nuclei are 1.5 × 10⁻⁷ in the dense PDR and 1 × 10⁻⁸ in the envelope. For C⁺, we adopt 2×10⁻⁴ in the envelope and 2×10⁻⁶ in the dense PDR. These choices successfully reproduce the observed intensity levels, indicating that the simple two-component model is globally realistic.

Figure B.1 shows the predicted radial profiles of the C¹⁸O 2–1 and [C II] 158 μm line intensities for the dense PDR-only model (blue squares), the envelope-only model (red squares), and the combined dense PDR plus envelope model (black squares). This figure is analogous to Fig. 12 for the o-H₂O 557 GHz line. As C¹⁸O 2–1 and [C II] 158 μm lines are collisionally excited, the envelope does not produce scattering of line photons emitted in the dense PDR (i.e., the red and black curves in the envelope are the same). Furthermore, the C¹⁸O 2–1 line is optically thin, so the emission from the dense PDR is not self-absorbed. For C⁺, the abundance in the dense PDR is much lower than toward the envelope, so any emission along the line of sight originates from the foreground envelope.

Fig. B.1 Two-component, nonlocal non-LTE radiative transfer model of the Horsehead applied to C18O and C+. The horizontal magenta line shows the observed maximum line intensities of the C18O 2–1 (upper panel) and [C II] 158 μm (lower panel) lines (Table 2).

Appendix C Weak dependence of the model results on G₀

We discuss here the dependence of the model line intensities on the interstellar radiation field scaling factor, G₀. As shown below, the best fit PDR wrapper model suggests a weak dependence of the model line intensities on G₀, with the best fit G₀ value that appears low compared to the nominal value of G₀ ~ 100 derived in previous studies of the Horsehead nebula (Abergel et al. 2003; Habart et al. 2005; Goicoechea et al. 2009; Santa-Maria et al. 2023), and implied by the spectral type of the ionizing star, σ Ori, O 9.5 V at a projected distance of 3.5 pc. While observations suggest that the ionizing star is located behind the nebula, the angle is quite oblique (see Fig. 13 of Abergel et al. 2024). Therefore, the actual distance is expected to be close to the projected distance. Zannese et al. (2025) also favor G₀ = 100 in their analysis of new JWST observations of H₂ line emission in the Horsehead nebula. Analysis of observations of tens of H₂ lines obtained with the IGRINS spectrometer on Gemini South (Piluso et al., in prep.) will provide additional strong constraints on the value of G₀ in this PDR.

In principle, the difference in the G₀ values may be due to the new modeling approach using the PDR wrapper or due to the particular selection of tracers used in our analysis. To investigate the G₀ dependence of our model results, we first ran the minimization using the standard 1D modeling approach, using the direct output intensities of the plane-parallel isobaric Meudon PDR model grid. We verified that the χ² distribution is very flat in the parameter space corresponding to the Horsehead physical conditions, and it favors low pressure and G₀ values. We further note that Philipp et al. (2006) also derived a low value of G₀ (~ 25 in Draine units, or ~ 43 in Habing units) from their isochoric models of the Horsehead nebula based on the same CO isotopolog and [C I] observations as used here, but not including the newer SOFIA [C II] observations.

To further investigate the dependence of model results on G₀, we constructed an extended grid of isobaric Meudon PDR models coupled with the wrapper tool, explicitly including G₀ as a free parameter. The explored parameter space covers:

• thermal pressures P_th = 0.5, 1, 2, 3, 4, 5, 6, 8, and 10 × 10⁶ K cm⁻³

• interstellar radiation fields G₀ = 30, 60, 70, 90, 100, and 110 in Habing units; and

• curvature radii R_C = 0.02, 0.05, 0.1, 0.2, 0.5, and 1.0 pc.

The resulting χ² maps (Fig. C.1), obtained from a minimization based on both line intensities and the spatial separation of emission peaks along the PDR (see Sections 4.3 and 4.4), provide strong constraints on the thermal pressure, which must lie between 2.5 × 10⁶ and ~ 6 × 10⁶ K cm⁻³. By contrast, the constraints on G₀ and R_C are much weaker: the χ² surface (bottom center panel) is relatively flat, with only a mild preference for lower values. The formal minimum is found for P_th = 3.2 × 10⁶ K cm⁻³, G₀ = 42, and R_C = 0.032 pc. However, as discussed above, such a low G₀ value is inconsistent with previous studies of the Horsehead, whereas G₀ = 100 is more consistent with the literature and remains perfectly compatible with our constraints.

To illustrate the low sensitivity of χ² to G₀, we show in Fig. C.2 an alternative model adopting literature values for the parameters, most notably G₀ = 100. This model provides a fit of comparable quality to the best-fit solution, while also being more consistent with previous work, with P_th = 4 × 10⁶ K cm⁻³, G₀ = 100, and R_C = 0.05 pc, i.e., about half the typical 0.1 pc width usually assumed in plane-parallel models of the Horsehead PDR. The best-fit model does improve the [C II] 158 μm intensity by ~ 30%, but this has no significant impact on the other tracers (excluding H₂O, which was not included in the minimization).

The weak sensitivity to G₀ is expected. The intensities of [C II] 158 μm, [C I] 609 μm, and low-J CO transitions vary only weakly with G₀ in the range 30 ≤ G₀ ≤ 100, at thermal pressures typical of the Horsehead (P_th ~ a few 10⁶ K cm⁻³). This behavior, already noted in earlier PDR studies (Tielens & Hollenbach 1985a,b; Kaufman et al. 1999), can be understood as follows: (i) in Horsehead conditions, the [C II] line is close to local thermodynamical equilibrium and therefore saturated, so its intensity is largely set by dust opacity and remains nearly constant with G₀; (ii) the [C I] column density is itself almost insensitive to the FUV field; and (iii) the low-J CO lines are optically thick, with intensities controlled by the temperature at the C⁺/C/CO transition, which varies little with G₀. As a result, the χ² maps display a nearly vertical structure in the (P_th, G₀) plane: the data constrain the thermal pressure robustly, but provide little leverage on the incident FUV field. Adopting G₀ = 100, a value supported by previous studies of the Horsehead, therefore does not affect our conclusions.

We further verified that models using the PDR wrapper with the distance information about the stratification of the peak emission of various tracers provide a tighter fit for the model parameters – in particular without the distance information, G₀ is much less constrained. This suggests that observations at higher angular resolution and including additional tracers would be required to better constrain G₀.

The objective of the present manuscript is to study the water vapor emission in the Horsehead nebula, not to determine the best value of G₀. While models with lower G₀ values provide a somewhat better fit to the water observations (Fig. C.2), the model water line intensity is stil a factor of 7 higher than the observations. Photon scattering in a low-density warm envelope surrounding the dense PDR is thus still required to explain our Herschel/HIFI observations of the 557 GHz water line, regardless of the exact value of G₀ value used.

Fig. C.1 Maps of χ2 values obtained from the minimization procedure using maximum intensities and spatial separations. The best-fit solution, corresponding to Pth = 3.2 × 106 K cm−3, G0 = 42, and RC = 0.032 pc, is marked by a red circle. Constraints on G0 and RC remain weak, with the χ2 surface appearing relatively flat (bottom center panel), except for very small RC values. For comparison, a model adopting literature values (Pth = 4 × 106 K cm−3, G0 = 100) is shown by a green circle in the G0–Pth panel.

Fig. C.2 Observed maximum intensities of various tracers (colored crosses) compared with model predictions (symbols with error bars). Squares show the best-fit model from the minimization described in Appendix C, while circles represent the model obtained in Section 4.4 using literature parameters with G0 = 100. All tracers are well reproduced, except for H2O, for which the models over-predict the observed maximum intensity by about an order of magnitude. The error bars of 40% on models represent geometrical uncertainties and the complexity of the parameter space.

All Tables

All Figures

Fig. 1 (Left panel) SPIRE 350 μm image of the Horsehead PDR. The white circle shows the FWHM SPIRE beam (25.2″). Black circles show the locations at which 557 GHz o-H2O spectra were obtained, with the circle size corresponding to the FWHM HIFI beam (38.1″). (Right panel) Spectra of the 557 GHz o-H2O and CO (2–1) lines across the PDR (black and blue curves, respectively), labeled by the right ascension offsets with respect to the reference position. The CO 2–1 spectra have been scaled down by a factor of 200. The bottom-right panel shows spectra of the ground state o– and p–H2O lines at the central position.

In the text

Fig. 2 Integrated line intensities of atomic and molecular tracers computed over a 9.2–11.8 km s−1 velocity range across the PDR as a function of right ascension offset at Δδ = 0. (Top) [C II], [C I], and H2O (black, blue, and magenta, respectively). (Bottom) CO (2–1), C18O (2–1), CO (4–3), and H2O (black, blue, red, and magenta, respectively). Line intensities have been normalized to their maxima toward the PDR, with the corresponding scaling factors as labeled. The data are shown at the native resolution of the images, with the corresponding FWHM beam sizes marked for each tracer. The size of the squares corresponds approximately to ±1σ H2O observational uncertainties.

In the text

	Fig. 3 Same as Figure 2 except that all tracers are convolved to the 38.1″ resolution of the HIFI H2O spectra.
In the text

Fig. 4 Illustration of the geometry of 1D PDR models made with the Meudon PDR code. We show the two-component radiation field composed of a beamed part representing the stellar radiation field and an isotropic part representing the ISRF. This configuration yields the typical chemical stratification, with the ionization front denoting the beginning of the PDR, followed by the dissociation front, or H/H2 transition, and the C+/C/CO transition deeper in the cloud. H2O appears even deeper.

In the text

Fig. 5 Illustration of the geometry of the PDR wrapper. (Left panel) Image of the Horsehead nebula as seen by the Euclid space telescope (Credits: ESA/Euclid/Euclid Consortium/NASA, image processing by J.-C. Cuillandre (CEA Paris-Saclay), G. Anselmi). Notice the reference frame, with the x-axis corresponding to the direction of the lines of sight. We rotate the nebula so that the x–axis points to the right. (Middle panel) Artistic interpretation of the rotation of the Horsehead nebula, generated using OpenAI’s DALL·E model based on the first panel image. (Right panel) Geometry of the wrapping 2D model, with the star radiation field coming from the top, and examples of Meudon PDR model output overplotted as vertical red lines. Each of these lines correspond to a PDR model as represented on Fig. 4. Notice that we use the same model for each position along the cut of the cloud. The multiple lines of sight correspond to different pixels along a vertical 1D cut in the Horsehead nebula. We emphasize that only the lines of sight shown as solid lines are correctly accounted for by the model, as the central section of the nebula outside of the crescent-shaped area is not included in the computations. The geometry is assumed to be circular with a radius of curvature that remains to be determined.

In the text

Fig. 6 Illustration of the wrapper geometry, showing RGB composite maps of carbon-bearing species for different radii of curvature. The color encodes the local abundances as follows: red for ionized carbon (C+), green for neutral atomic carbon (C), and blue for carbon monoxide (CO). Areas where these species coexist appear as color blends (e.g., cyan for C and CO, yellow for C+ and C, white for equal contributions). The left panel shows a curvature radius of RC = 0.05 pc, the middle to RC = 0.2 pc, and the right to RC = 0.5 pc. These maps highlight the local spherical curvature in the plane of the sky, with the observer located at the bottom of each panel and the stellar radiation field entering from the right. Radiative transfer is computed along vertical lines from top to bottom in this configuration. The abundances originate from a single Meudon PDR model (Pth = 4 × 106 K cm−3, G0 = 100). The usual edge-on representation of abundances of this model can be found in Figure 10.

In the text

	Fig. 7 A map of χ2 values obtained from the minimization procedure using maximum intensities and spatial separations. The best-fit model, corresponding to Pth = 4.1 × 106 K cm−3, and RC = 0.057 pc, is denoted by the orange circle on the map. The constraint on RC remains weak, with the χ2 surface relatively flat in a wide valley.
In the text

Fig. 8 (Top) Spatial profiles of normalized emission line intensities observed in the Horsehead nebula, compared with predictions from a Meudon PDR model with Pth = 4 × 106 K cm−3 and G0 = 100, wrapped in a locally spherical cloud with a curvature radius RC = 0.05 pc. These conditions correspond to commonly adopted literature values. The left panels show the raw profiles from the wrapping procedure, while the right panels present the same profiles convolved to the native angular resolution of each tracer. In those figures, the radiation field comes from the right. (Bottom) Observed maximum intensities of various tracers (colored crosses) compared to model predictions (circles with error bars). All tracers are well reproduced, except for H2O, where the model overestimates the observed maximum intensity by an order of magnitude.

In the text

Fig. 9 Spatial profiles of normalized emission line intensities observed in the Horsehead nebula, compared to predictions from Meudon PDR model with Pth = 4 × 106 K cm−3 and G0 = 100, wrapped in a locally spherical cloud with curvature radius RC = 0.05 pc. The profiles are convolved to the native angular resolution of the tracer maps. On this figure, the radiation field comes from the right. We observed a significant shift between the observed H2O profile and the model, the latter being much deeper in the cloud.

In the text

Fig. 10 Abundances of key species with respect to hydrogen nuclei (top panel) and gas temperature (bottom panel) profiles of the Horsehead nebula model with Pth = 4 × 106 K cm−3, G0 = 100, computed with the Meudon PDR code (no wrapping or convolution procedure is applied to these profiles). The H/H2 transition can be seen to take place at around AV ~ 0.06 mag, the C+/C/CO transition is at AV ~ 1.7 mag, and the peak of H2O abundance is reached at an AV ~ 3 mag.

In the text

Fig. 11 Two-component, nonlocal radiative transfer models adapted to the physical conditions in the Horsehead. The upper panel shows the o-H2O 557 GHz synthetic line profiles (continuous curves) and observed line profiles (histograms) at the emission peak. The lower panel shows the predicted excitation temperatures of the 110–101 transition. The blue curves represent the model of the dense PDR component alone. The red curves represent the model of the low-density envelope alone, while the black curves represent the two-component model.

In the text

	Fig. 12 Same as the lower panel of Fig. 11 but showing the predicted radial o-H2O 110–101 integrated line intensities. The magenta-filled circles represent the observed line intensities across the Horsehead PDR, with the peak position at Δα = +20″.
In the text

	Fig. A.1 Horizontal geometry used for the wrapping of Meudon PDR models. The beamed ISRF remains perpendicular to the line of sight.
In the text

	Fig. B.1 Two-component, nonlocal non-LTE radiative transfer model of the Horsehead applied to C18O and C+. The horizontal magenta line shows the observed maximum line intensities of the C18O 2–1 (upper panel) and [C II] 158 μm (lower panel) lines (Table 2).
In the text

Fig. C.1 Maps of χ2 values obtained from the minimization procedure using maximum intensities and spatial separations. The best-fit solution, corresponding to Pth = 3.2 × 106 K cm−3, G0 = 42, and RC = 0.032 pc, is marked by a red circle. Constraints on G0 and RC remain weak, with the χ2 surface appearing relatively flat (bottom center panel), except for very small RC values. For comparison, a model adopting literature values (Pth = 4 × 106 K cm−3, G0 = 100) is shown by a green circle in the G0–Pth panel.

In the text

Fig. C.2 Observed maximum intensities of various tracers (colored crosses) compared with model predictions (symbols with error bars). Squares show the best-fit model from the minimization described in Appendix C, while circles represent the model obtained in Section 4.4 using literature parameters with G0 = 100. All tracers are well reproduced, except for H2O, for which the models over-predict the observed maximum intensity by about an order of magnitude. The error bars of 40% on models represent geometrical uncertainties and the complexity of the parameter space.

In the text