© The Authors 2024

1 Introduction

The interaction of high-mass stars with their environment determines the rate of star formation and star formation efficiency of molecular clouds, and regulates the evolution of galaxies. Mechanical and radiative energy input from high-mass stars can stir up, ionize, and heat the gas, leading to the disruption of the natal molecular cloud within a few free-fall times (see Schneider et al. 2020, and references therein.) However, high-mass stars can also provide positive feedback, by compressing the gas, thus allowing gravity to overcome the turbulent and magnetic support that had prevented collapse and star formation.

The aim of the SOFIA FEEDBACK legacy program (Schneider et al. 2020; Tiwari et al. 2022; Schneider et al. 2023) was to quantify the effects of stellar feedback in a range of environments ranging from single OB stars, to small groups, to rich young clusters, and to mini-starburst complexes, by determining the physical conditions using observations of the far-infrared (FIR) [C II] and [O I] fine structure lines, combined with H₂ line observations from Spitzer, and ground-based millimeter wave observations.

Photon-dominated regions (PDRs), the surface layers of molecular clouds exposed to strong fluxes of UV photons, dominate the IR spectra of star-forming clouds and galaxies, and provide a unique tool to study in detail the physical and chemical processes in diverse environments ranging from diffuse clouds, molecular cloud surfaces, globules, protoplanetary disks, planetary nebulae, and starburst galaxies. The JWST Early Release Science program PDRs4All (Berné et al. 2022; Peeters et al. 2024; Berné et al. 2024) provides template datasets designed to identify key PDR characteristics in JWST spectra of star-forming regions in our Galaxy and beyond. As discussed in Berné et al. (2022), PDRs are strong emitters in the near and mid-IR, where their key diagnostic signatures include (a) the Aromatic Infrared Bands (AIBs) attributed to polycyclic aromatic hydrocarbons (PAHs), (b) continuum emission from very small carbonaceous grains (VSGs), (c) ro-vibrational and pure rotational lines of H₂, and (d) emission from atomic ions (e.g., S⁺, Si⁺, Fe⁺).

Hydrodynamical numerical simulations (e.g., Mackey et al. 2010) show that shadowing of ionizing radiation by an inhomogeneous density field can form so-called elephant trunks (ETs), or pillars of dense gas without the assistance of self-gravity or of ionization front and cooling instabilities. In some cases, the head of the ET is shown to recede from the radiation source more rapidly than the shadowed trunk which has not been exposed to radiation. For dynamically forming ETs, gas streams into the shadowed trunk past the head, and is thus moving faster than, or at a comparable speed to, the pillar’s head. This could result in gas compression, triggering collapse and formation of next generation of stars. A combined Chandra/optical study aimed at quantifying triggered star formation in IC 1396A (Getman et al. 2012) also argued for a radiation-driven implosion process persisting over several million years. These authors argued that the contribution of the triggered star formation for the entire HII region exceeds 14–25% today, and may have been higher over the lifetime of the HII region. Such triggering on the periphery of HII regions may thus be a significant mode of star formation in the Galaxy.

Magnetic fields have been suggested to have a significant influence on the development of pillars around H II regions (Williams 2007; Mackey et al. 2011). In the case of M 16, submillimeter dust polarization observations have shown that the magnetic field lines are aligned with the pillars, and perpendicular to and decoupled from the field of the surrounding photo-ionized cloud (Pattle et al. 2018).

The best-known example of such a region, M 16, has recently been imaged using NIRSpec and MIRI, revealing an incredible level of detail in these “Pillars of Creation” (Hester et al. 1996). Karim et al. (2023) presented SOFIA and groundbased spectroscopic observations of this region. They infered that the ionized, atomic, and molecular phases are in pressure equilibrium if the atomic gas is magnetically supported. The dense clumps at the tops of the pillars are currently supported by the magnetic field, but are likely to collapse within the photoveaporation time scales. The present study focuses on IC 1848, a similar region with a lower-UV-illumination, for which an extensive set of supporting FIR continuum and line observations from Spitzer, SOFIA, and Herschel already exists.

IC 1848 is a large HII region in Cassiopea, forming part of the W5 HII region-molecular cloud complex in the Perseus Arm, at a distance 1.9 kpc (Ishida et al. 1970), consistent within uncertainties with the recent distance estimate of 2075+44−42 pc to the Cas OB 6 association, formed by W3/W4/W5, derived from the Gaia EDR3 release (Maíz Apellániz et al. 2022). Past observations revealed three bright rimmed clouds on the southern edge of the HII region (Thompson et al. 2004) at a projected distance of 13 pc from HD17505, the primary O6V ionizing star. These filaments, referred to as G1 – G3 (Figure 1), are associated with recent or ongoing star formation, as evidenced by the presence of submillimeter cores and multiple YSO’s identified in 2MASS images. This strongly suggests that UV illumination may have induced the collapse of the dense molecular cores at the head of the three clouds and Thompson et al. (2004) argued that the overall morphology of the IC 1848 clouds is reasonably consistent with radiative-driven implosion models of cometary globules.

A recent study of IC 1848 (Chauhan et al. 2011) suggested that the filaments may have been caused by hydrodynamical instability of the ionization/shock front of the expanding H II region, as similar structures often show up in recent numerical simulations of the evolution of H II regions. These authors further hypothesized that this hydrodynamic instability mechanism makes a third mode of triggered star formation associated with H II regions, in addition to the two previously recognized mechanisms: collect-and-collapse of the shell accumulated around an expanding H II region, and radiation-driven implosion of bright rimmed clouds that were pre-existing cloud clumps. Accurate characterization of the physical conditions in the gas, through observations combined with state-of-the-art Meudon PDR and Cloudy modeling, would provide good quantitative constraints for such models.

To characterize the physical conditions and the velocity field in the IC 1848 filaments, observations of [C II] 158 µm, and [O I] 63 and 145 µm fine structure lines were carried out using the upGREAT instrument on SOFIA. Full velocity-resolved [C II] maps of the three filaments were obtained, along with [O I] observations at selected positions.

Fig. 1 Spitzer/IRAC 8 µm image of the three filaments in IC 1848 (G1–G3, left-to-right, respectively). Bright IRAC compact sources are clearly visible at the head of the filaments (white saturated spots in the image). Cyan contours show the [C II] 158 µm integrated line intensity observed with SOFIA/upGREAT. Contour levels are 8, 12, 18, and 27 K km s−1, and typical 1σ noise levels are 3.3, 1.8, and 2.4 K km s−1 in G1–G3, respectively. Blue stars mark the location of the three brightest 2MASS/WISE sources at the heads of the filaments, and the white arrow shows the direction toward the O6V ionizing star HD 17505.

2 Observations

During six flights between December 2019 and April 2022, all departing from Palmdale, CA, the three trunks G1 – G3 (Figure 1) were observed with the array spectrometer upGREAT¹ (Risacher et al. 2018) aboard SOFIA (Young et al. 2012). All observations were carried out on the GREAT Consortium time (total flight time 7.5 hours, including all overheads). First, the filaments were mapped in fast total power on-the-fly mode with both polarizations of the 7-pixel low- frequency arrays (LFA) tuned to the carbon fine-structure line [C II] ³P_3/2 − ³P_1/2 158 µm (ƒ₀ = 1900.53690 GHz; Cooksy et al. 1986). The three mapped fields were over-sampled at 3″ spacing, with 0.2 sec integration time per dump. The high- frequency array (HFA), tuned to the [O I] ³P₁ − ³P₂ 63 µm (ƒ₀ = 4744.77749 GHz; Zink et al. 1991) transition, was operated in parallel, but post-flight analysis revealed a too weak line intensity for quantitative analysis. Therefore, we decided to revisit the three strongest [C II] condensations, in each filament, with targeted, deep(er) integrations in rather slow chopped mode (rate 0.626 Hz) for better stability. For these observations the LFA was operated in split mode, with the LFA-V tuned to [O I] ³P₀ − ³P₁ 145 µm (ƒ₀ = 2060.06886 GHz), while the LFA-H was still recording [C II]. The HFA was running in parallel, but due to the frequency dependent array geometry, only the central pixels of the two arrays were co-aligned on the sky.

Half-power main beam sizes (θ_mb: 14″.1 for [C II], 13″ for [O I] 145 µm and 6″.3 for [O I] 63 µm) and efficiencies (<η_mb> = 0.64 for LFAV, 0.66 for LFAH, and 0.68 for HFA, respectively), derived from planet observations, were provided by the GREAT consortium. The instrument parameters are described in Risacher et al. (2018). Both spectrometer arrays used NbN-based hot electron bolometer mixers, driven by either solid-state (LFA) or quantum cascade laser (HFA) local oscillator sources. The LFA operated 14 pixels, seven in each polarization, in a hexagonal arrangement with 31″.8 radial spacing around a central pixel. The HFA had the same symmetry, using the V polarization only, with 13″.8 radial spacing. Co-alignment between the central pixels was generally better than 1″ (Risacher et al. 2018). For this experiment we tracked on the central HFA pixel. Pointing was established by the telescope operators on nearby optical reference stars, to an accuracy of 1–2″. Prior to each flight series, the optical axis of the GREAT instrument had been aligned to these imagers by observations of planets. Fast Fourier Transform Spectrometers (FFTS4G, updated from Klein et al. 2012) offer 32k channels across the 0–4 GHz intermediate-frequency bands.

Raw data were amplitude calibrated with the KOSMA kali- brate software package (version 2022.04), following (Guan et al. 2012). The data were corrected for atmospheric extinction and calibrated in T_mb. The atmospheric transmission was smooth near the targets’ velocities, except for the pointed observations of [O I] 145 µm towards G3 (in April 2022) that were made impossible because of strong blending with a prominent telluric ozone feature.

Further processing of the calibrated data was performed with the GILDAS software². An initial inspection of the data cube was done in order to identify and eliminate unusable spectra with artefacts. Spectra with baseline ripples were identified by comparing the baseline noise with the radiometric noise deduced from the calibration and discarded. Once in main beam temperature, the spectra were smoothed to a resolution between 0.1 and 0.15 km s⁻¹, and baselines of orders one or two were removed.

3 Results and discussion

3.1 Morphology and kinematics

Figure 1 shows excellent correlation between the IRAC 8 µm continuum emission in IC 1848 (color image), likely dominated by the PAH emission, and the [C II] line emission (cyan contours). Bright IRAC point sources are located at the head of each pillar, which coincide with 2MASS point sources classified as YSO or YSO Candidate. Their 2MASS and WISE magnitudes are listed in Table 1. From their WISE colors (Koenig & Leisawitz 2014; Kang et al. 2017) the sources at the head of the G1 and G2 filaments can be identified as Class I, while that at the head of the G3 filament is at the boundary between Class I and Class II.

To estimate bolometric luminosities of the embedded stars from their 2MASS magnitudes, we used Eq. (5) of González Hernández & Bonifacio (2009). We used the K_S magnitudes, which are least affected by extinction, and adopted the Vega K magnitude and absolute integrated flux from Table 3 of González Hernández & Bonifacio (2009). The resulting bolometric luminosities are 0.3–0.9 L_⊙ (Table 1). These values are not corrected for extinction.

Figure 2, top panel, shows the morphology of the [C II] line emission toward the brightest G1 filament. A bright, narrow ridge of emission is seen on the west side of the filament, with weaker emission extending toward the east. Example [C II] spectra toward selected positions marked A–F are shown on the top-right. With the exception of the position F at the base of the filament, which shows a broad line with two velocity components, the remaining spectra show consistent velocities, with the head of the filament (position A) displaying slightly red- shifted emission compared to the other positions, possibly owing to the radiation pressure exerted by the impinging UV photons or the rocket effect of the evaporating gas (Bertoldi & McKee 1990), which causes the cloud to be accelerated away from the ionizing star. Bertoldi & McKee (1990) estimates the characteristic rocket velocity as 10–5 km s⁻¹. This is higher than the ≃ 1 km s⁻¹ velocity difference between the surface and the center of the G1 filament (Figure 2, lower-left), possibly due to projection effects. The [O I] velocities at the head of the G1 filament are consistent with [C II] (Figure A.6, positions G, H).

To investigate further the kinematics of the gas across the filament, we computed the first and second moment maps of the [C II] line emission. The bottom-left panel shows the mean [C II] velocity computed as $\overline{v}={{\displaystyle \sum }}^{\text{}}{v}_i{T}_i/{{\displaystyle \sum }}^{\text{}}{T}_i$, and the bottomright panel shows the velocity dispersion, ${\sigma }_v=\sqrt{\overline{v^2}-{\overline{v}}^2}$. The mean line velocity increases by up to ≃1 km s⁻¹ at the head and near the edges compared to the center of the filament, while the velocity dispersion remains uniform within the observational uncertainties. A hint of this velocity shift can also be seen in the position A spectrum in Figure 2. The base of the G1 filament is characterized by a significantly lower mean velocity and higher velocity dispersion, owing to the presence of a distinct velocity component at −40 km s⁻¹ (Figure 2, position F), which is likely a part of the extended molecular shell that can be seen in the IRAC image (Figure 1). Outside of this region, the velocity dispersion is quite uniform along the filament, in contrast to Pillar 2 in M 16, where the dispersion increases toward the tail (see Fig. 8 of Karim et al. 2023). This is likely due to the difference in the UV illumination between the two regios. The three M 16 pillars are illuminated by a handful of bright O5– 7 stars located at a ~3 pc projected distance (Hillenbrand et al. 1993). Karim et al. (2023) estimated the far-ultraviolet (FUV) field strength toward the M 16 pillars is in the range 1500–2500 Habing units. Consequently, the UV photons have smaller influence on the kinematics of the gas and the resulting degree of turbulence in IC 1848 compared to M 16.

Moment maps of the G2 and G3 filaments are shown in Figures A.1–A.2. Filament G2 shows an elliptical morphology, with the west side (Figure A.1, positions A and C) significantly redshifted compared to the east side (positions B, D, and E). There is some evidence for lower velocity dispersions on the west side, but the uncertainties are large. Filament G3 splits into two distinct parts, the head and the base, which display large differences in the mean velocity and velocity dispersion (Figure A.2). The [C II] spectrum at the base is double peaked, indicating the presence of a separate velocity component at about −43 km s⁻¹.

The kinematic structure of the three filaments is further displayed in the velocity channel maps (Figures A.3–A.5, filaments G1–G3, respectively). Based on the morphology of the Hα emission, Thompson et al. (2004) argued that although the three clouds lie close together on the sky, this arrangement may be a chance superposition, with G1 located behind the illuminating star and G2/G3 in front. Such a geometry could explain the ~5 km s⁻¹ velocity differences between G1 and the other two clouds. However, the observed velocity gradients observed along all three filaments, with the [C II] emission at the head being redshifted with respect to the tail, would seem to indicate that the three filaments are collocated on the front side of the illuminating star (see, e.g., Bertoldi & McKee 1990; Pound 1998). As discussed in Sect. 3.4, the observed [C II]/[O I] 63 µm line intensity ratio does not seem consistent with back-side illuminated PDR models of G2 and G3. The exact geometry of the three filaments with respect to the illuminating star remains, therefore, ill-defined.

Fig. 2 [C II] emission toward the G1 filament. (Top) Color image of the [C II] integrated line intensity in the −40 to − 33.5 km s−1 velocity range. The right panel shows spectra at selected positions labeled A–F in the image. Blue circles correspond to the FWHM SOFIA beam size at the [C II] frequency. (Bottom) Maps of the [C II] line center velocity and velocity dispersion (color images, left and right panels, respectively). Only pixels with [C II] intensities above 3.5σ (11.6 K km s−1) are shown. Black contours show the integrated line intensity, with the same contour levels as in Figure 1. Black star marks the location of the bright IRAC compact source in each panel.

3.2 Hydrogen column densities

In addition to the SOFIA and Spitzer data, far-infrared continuum images of IC 1848, obtained using PACS and SPIRE instruments, are available in the Herschel Science Archive, which allow characterization of the continuum spectral energy distribution (SED) over a broad wavelength range to constrain the dust temperature and the H₂ column density required for the PDR modeling (Sect. 3.4). Figure 3 shows the continuum SEDs at positions A–F in the G1 filament. All continuum images have been convolved to the 36″.7 resolution (full width at half maximum, FWHM) of the longest-wavelength 500 µm SPIRE channel. Blue curves show single-component modified blackbody fits, I_V = B_V(T_d) × [1 − exp(−τ₃₅₀ × (350µm/λ)^β)], to the PACS and SPIRE data only. Best fit parameters, the dust temperature (T_d), the 350 µm optical depth (τ₃₅₀), and the grain emissivity exponent (β) are listed in Table 2.

The SED toward position A is well described by a scaled blackbody (β = 0). This is consistent with the flux, even at the longest SPIRE wavelengths, being dominated by the stellar source seen in the IRAC images rather than dust emission. The remaining positions have steeper dust emission spectra characterized by average dust temperatures of ~28 K, low 350 µm optical depths, and shallow values of β ~ 0.8.

The H₂ column densities in a 36″.7 beam were computed using the formula N(H₂) = 2aρR_𝑔d /3m_H × τ₃₅₀/ Q₃₅₀ (Lis & Goldsmith 1990), where a = 0.1 µm is the mean grain radius, ρ = 3 g cm⁻³ is the mean density, R_𝑔d = 100 is the gas-to- dust mass ratio, m_H is the hydrogen mass, and Q₃₅₀ is the 350 µm grain emissivity coefficient. Extrapolating the 125 µm grain emissivity coefficient of Hildebrand (1983) (7.5 × 10⁻⁴) with a v² frequency dependence gives Q₃₅₀ = 1 × 10⁻⁴, corresponding to the grain mass opacity coefficient κ₃₅₀ = 3Q/4/a/ρ/R_d𝑔 = 0.025 cm² g⁻¹, including the gas-to-dust ratio (see Lis et al. 1998). The resulting H₂ column densities are listed in Table 2 together with the corresponding visual extinction values, computed assuming the conversion factor of N_H = 2.21 × 10²¹ cm⁻²¹ A_V(mag) (Güver & Feryal 2009). The typical visual extinction values characterizing the [C II] emitting gas are moderate, ranging from A_V ~ 10 near the head of the filament to ~4 at the base. The characteristic width of the G1 filament is ~90″, corresponding to 0.8 pc, or 2.6 × 10¹⁸ cm at a distance of 1.9 kpc. Assuming the same depth along the line of sight, the H₂ column densities toward the head of the filament correspond to a mean volume density of ~4500 cm⁻³ . This should be considered a lower limit, as narrower filamentary structures can be seen in the [C II] and dust continuum images. These values were used as input for the PDR models described below.

Modified blackbody fits toward positions A–E in the G2 and G3 filaments are shown in Figure A.7 in the Appendix, and best fit parameters are listed in Table 2. The dust temperatures, H₂ column densities, and the corresponding visual extinction values are comparable to those in the G1 filament. The dust temperature, A_V, and average density estimates derived here are in a reasonable agreement with those derived by Thompson et al. (2004) based on SCUBA and IRAS HIRES fluxes (see Table 5 of Thompson et al. 2004).

In addition to the H₂ column densities derived from dust SED fits, Table 2 lists the [C II] line intensities and the corresponding C⁺ column densities at the same positions derived using equations (26) and (43) of Goldsmith et al. (2012). PDR models (Sect. 3.4.2) suggest that the [C II] emission originates in the outer A_V ~ 1.5 molecular layer characterized by gas temperatures between ~50 and 160 K. Consequently, we used an average temperature of 100 K and an average H₂ density of 4500 cm⁻³ derived above in the C⁺ column density calculations. Typical C⁺ column densities in the G2 and G3 filaments are of order 1 − 2 × 10¹⁷ cm⁻², while those in the G1 filament are a factor of ~2 higher (Table 2). Assuming a C⁺ abundance of ~2.5 × 10⁻⁴ with respect to H₂ in the outer molecular layer, as suggested by PDR models, the [C II] emission traces on average ~17% of the total H₂ column density derived from the dust SED fits. The C⁺ column densities in IC 1848 are toward the low end of the values derived by Karim et al. (2023) in the M 16 pillars, ~2 × 10¹⁷ − 4.7 × 10¹⁸ cm⁻². This is once again likely due to the difference in the FUV illumination between the two regions.

Fig. 3 Spectral energy distributions toward positions A–F in the G1 filament based on Spitzer and Herschel observations. Blue curves show the modified blackbody fits to the Herschel PACS and SPIRE fluxes from which the dust temperatures and 350 µm optical depths used for the column density determination are derived. Typical flux calibration uncertainties are ~5%.

3.3 Atomic oxygen observations

In addition to the fully-sampled [C II] maps, the 63 µm and 145 µm [O I] fine structure lines were observed at selected positions. [O I] spectra toward the G1 filament with detections above the 4σ level for the integrated line intensity are shown in Figure A.6 (upper panel). The 63 µm line is detected at 6 positions, labeled G–L, while the 145 µm line is detected at 3 positions, labeled M–O (see Figure A.6). The [O I] integrated line intensities are listed in Table 3, together with [C II] line intensities at the same positions extracted from the fully-sampled map.

The average [O I] 63 µm/[C II] intensity ratio (in K km s⁻¹) for positions G–L is ~0.10, while the average [O I] 145 µm/[C II] intensity ratio for positions M–O is ~0.20. This suggests that the [O I] 145 µm line is typically stronger than the 63 µm line.The [O I] 63 µm line is detected at four positions in the G2 filament (see Table 1 and Figure. A.6, lower panel), and is not detected in the G3 filament. The [O I] 63 µm/[C II] intensity ratios in the G2 filament are in the same range as those observed in G1.

Given the sparse and limited spatial coverage, limited signal-to-noise ratio of the [O I] 145 µm spectra, and the difference in the beam size at the frequencies of the two oxygen lines, the [O I] 63 µm/145 µm intensity ratio cannot be determined with high accuracy. However, positions N and I are located nearby, 14″ apart along the ridge of the [C II] emission and the [O I] 63 µm/145 µm intensity ratio at these positions in ~0.7. An independent estimate of the [O I] 63 µm/145 µm is obtained by averaging [O I] spectra near the head of the G1 filament, within the white rectangle in Figure A.6. The [O I] 63 and 145 µm lines are both detected with a signal-to-noise ratio of 14 and 8, respectively. The resulting [O I] line intensities are reported in Table 3 (lines labeled “G1 Head”) together with the average [C II] integrated line intensity within the white rectangle.

The [O I] 63 µm/145 µm ratio is 0.67, consistent with the values derived above, further supporting the conclusion that the [O I] 63 µm line is weaker than the 145 µm line.

3.4 PDR modeling

To model the fine structure line emission in IC 1848 we used the Meudon PDR code, a publicly available, state-of-the-art astro-chemistry code (Le Petit et al. 2006). Our objective was to verify whether the observed line intensities are consistent with the geometry and the UV illumination of IC 1848. The Meudon code is well suited for our analysis, since it is one of very few astrochemistry codes that computes self-consistently the radiative transfer, the chemistry both in the gas phase and on grain surfaces, as well as the thermal processes and thermal balance at each point within the cloud. The code simulates a stationary plane-parallel slab of gas and dust illuminated on both sides by the radiation field from the far UV (FUV) to the radio domain. This radiation field is discretized on a wavelength grid (typically 500 000 points) and the radiative transfer equation is solved at each wavelength, and at each position in the cloud, considering (a) absorption in the continuum by the distribution of grains and by the ionization of atoms such as carbon, (b) absorption in the atomic lines of hydrogen, the Lyman and Werner transitions of H₂ as well as the pre-dissociating lines of CO isotopologues (in total more than 30 000 UV absorption lines are considered), (c) backscattering by dust, and (d) re-emission of photons as the fluorescent UV photons emitted by H and H₂ (see Goicoechea & Le Bourlot 2007).

This detailed solution of the radiative transfer equation gives the spectral energy density at all wavelengths and positions within the cloud, which is then used to compute locally precise photo-destruction rates by direct integration of the photo-absorption cross-sections (Heays et al. 2017) with the monochromatic energy density. This means that photo-destruction rates include the UV shielding of molecules by H Lyman transitions and H₂ Lyman and Werner bands. So, whereas most of other astrochemical codes rely on basic approximations (simple exponential fits) to treat absorption of UV photons, here, thanks to the detailed treatment of the photo-destruction processes, the destruction rates of molecules are determined directly, and so their densities as a function of the position are computed more precisely.

For the chemistry, the code considers approximately 200 chemical species linked by a chemical network of ~2000 chemical reactions. In the gas phase, in addition to the standard exothermic reactions that constitute the bulk of chemical pathways, state-to-state reactions are considered for key endothermic reactions. For the grain surface chemistry, a power law grain size distribution is considered (Mathis et al. 1977). The temperature for grains in each size bin is computed balancing absorption and emission of photons. The processes that are considered are the adsorption from the gas, reactions in the bulk and on the surface considering thermal and tunneling effects, desorption processes (thermal and non-thermal: photo-desorption, cosmic ray desorption, chemical desorption, i.e., release of molecules into the gas phase due to the exothermicity of some reactions). The Meudon PDR code has another strong advantage over most other astro-chemistry codes: it computes very precisely the formation of H₂ on dust grains. Indeed, it can either consider Langmuir-Hinshelwood and Eley-Rideal processes with the rate equation formalism (Le Bourlot et al. 2012), or even consider stochastic processes on dust (Le Petit et al. 2009; Bron et al. 2014, 2016).

3.4.1 ISMDB models

A particularly useful online tool is the interstellar medium dataBase (SMDB)³, which allows quick and easy access to pre-calculated theoretical models for various ISM conditions and objects, including models produced by the Meudon PDR code. Presently, ISMDB gives access to several thousand PDR models covering a wide range of astrophysical conditions. To prepare and interpret observations, ISMDB provides an inverse search service to allow queries on observables to find the best models that match observations, and gives access to line intensity maps in a parameter space (e.g., the total hydrogen nuclei volume density n_H = n(H) + 2 × n(H₂), the ISRF scaling factor χ in Mathis units⁴, Mathis et al. 1983). Once the best fit model is found, model results can be downloaded without the need for running the code on a local machine.

To analyze the observed fine structure line intensities in IC 1848, we used the ISMDB inverse search tool based on a grid of isochoric (constant density) models computed using the latest public version of the Meudon PDR code (1.5.4). We assumed front side illumination (enhanced illumination of side of cloud facing the observer), and the total thickness of the cloud corresponding to A_V = 10 mag, as derived above from the SED analysis. We note that the exact value of A_V is not important for the computations, as the fine structure line emission originates from the surface of the PDR (see Sect. 3.4.2 below). We verified that models with A_V = 5 and 10 mag predict essentially the same line intensities. The two key parameters varied in the fit are the gas density and the radiation field scaling intensity factor, χ, with the observed fine structure line intensities being input observables. We computed averages of the individual [C II] and [O I] 63 µm line intensities at positions G–L, and of the [C II] and [O I] 145 µm line intensities at positions M–O, and use these values as input for the ISMDB fit, assuming a 30% uncertainty as an estimate of both statistical and calibration uncertainties.

The χ² for the ISMDB fit to the combined data set (Figure 4; color map and thin black contours) shows a broad, shallow distribution with no unique solution. Since we do not have observations of the two [O I] fine structure lines at the same positions, we then fit separately the combinations of [C II]/[O I] 63 µm at positions G–L and [C II]/[O I] 145 µm at positions M–O. A fit to the observed [C II] and [O I] 145 µm line intensities at positions M–O gives a unique solution with n_H ≃ 3 × 10⁴ cm⁻³ and χ ≃ 200–250 (blue dot), while the observed [C II] and [O I] 63 µm line intensities at positions G–L are consistent with a range of parameters within the elongated red area. The two regions do not overlap; there is no unique solution that is consistent with both data sets. The [O I] 145 µm observations suggest higher gas densities compared to the [O I] 63 µm observations. For χ = 200–250, the best fit for positions G–L has a density n_H ≃ 1–3 × 10³ cm⁻³, an order of magnitude lower than the value derived from the [O I] 145 µm fit. These results show that the observed [O I] 63 µm/145 µm ratio cannot be explained easily with the PDR models. To bring the two solutions into agreement, we must lower the [O I] 145 µm line intensity by about a factor of 3. The best fit to the combined data set with such scaling is shown as the narrow magenta elongated area. We note that the mean gas density n_H ≃ 9 × 10³ cm⁻³ derived from the dust SED analysis is consistent with the southern tip of the best fit area; the corresponding χ value is ≃100.

The observed PACS fluxes toward the IC 1848 filaments allow an independent estimate of the UV radiation field strength (Kramer et al. 2008; Roccatagliata et al. 2013). We considered as an example position B near the head of the G1 filament, but away from the embedded IRAC point source. The PACS 70 and 160 µm fluxes at this position are 212 and 598 MJy sr⁻¹, respectively (Fig. 3). The corresponding effective filter widths read from the filter curves in the PACS online documentation are ~1.3 × 10¹² and 1.1 × 10¹² Hz, respectively. The UV field strength can then be estimated using Eq. (2) of Roccatagliata et al. (2013) by adding the 70 and 160 µm PACS intensities (Schneider et al. 2016). The resulting value is χ = 74 Habing units, corresponding to 96 Mathis units, consistent with the southern tip of the elongated best fit area in Figure 4.

Fig. 4 χ2 values for a grid of isochoric ISMDB models compared to the observed fine structure line intensities in the G1 filament. Blue circle and the elongated red area show best fits, respectively, to the [C II]/[O I] 145 µm data at positions M–O and [C II]/[O I] 63 µm data at positions G–L. The elongated magenta area shows the best fit to all data with [O I] 145 µm line intensities scaled down by a factor of 3, as described in the text. All models assume normal incidence. Increasing the incidence angle to 60° provides comparable fit, but shifts the best-fit solution to higher densities and lower χ values (magenta area).

Fig. 5 Temperature, density, and abundance profiles in the outer PDR layers of Model 1027.

3.4.2 PDR models illuminated by a star

The ISMDB models discussed above assume isotropic radiation and a normal incidence angle. To investigate the model results in more detail, we ran a series of offline PDR models with different geometries. The results are summarized in Table 4. The input parameters included the spectral type of the illuminating star (O5V is the closest spectral type considered in the Meudon code), distance D (negative values imply front-side illumination and positive values back-side illumination), inclination angle θ, A_V, and gas density n_H. The UV radiation field of an OV5 star at a distance of 13 pc corresponds to χ ~ 35 Mathis units (F. Le Petit, private communication), lower than the values derived above. We note, however, that using a stellar spectral type as input parameter for the Meudon PDR code is different from using an ISRF scaling factor χ in the ISMDB analysis. χ is a scaling factor for an isotropic ISRF, whereas the radiation field from a star is plane-parallel. A plane-parallel (stellar) UV radiation field penetrates deeper into the cloud than an isotropic radiation field, resulting in higher intensities of FIR fine structure lines. The difference in the predicted line intensities can sometimes be significant.

Table 4 shows predicted fine structure line intensities and line ratios for several PDR models with different parameters computed using an off-line version of the Meudon PDR code (1.5.4). Models 1027–1028 correspond to normal incidence at the closest possible (projected) distance from the illuminating star (highest possible illumination). Models 1029–1030 correspond to more physical geometries, with inclination angles θ = 30 and 60°, and the distance increased accordingly as the projected distance divided by sin θ. Models 1027 and 1029 are consistent with the observed line intensities of the [C II] and [O I] 63 µm lines, assuming front side illumination, although the [O I] 63 µm/[C II] line ratio is slightly lower than the observed values. However, the [O I] 145 µm/[C II] line ratio for all models is an order of magnitude lower than the observed values, as suggested earlier by the ISMDB models. The temperature, density, and abundance profiles in the outer PDR layers of Model 1027 are shown in Figure 5.

While models with back-side illumination produce [O I] 63/145 µm line ratios below unity, which is explained by foreground absorption of the optically thick [O I] 63 µm line (Goldsmith et al. 2021), the predicted intensities of both lines are much too low compared to [C II]. Therefore, we do not consider back-side illumination and optically thick [O I] 63 µm emission as a viable solution for the IC 1848 pillars.

Since under most conditions encountered here, the [O I] 145 µm line has a higher critical density than do the [C II] and [O I] 63 µm lines, we considered two additional models (1033–1034) with the same geometry as models 1029-1030, but with density increased by a factor of 5. While the intensity of the [O I] 145 µm line increases by a factor of 1.7, it is still well outside of the observed range. The [O I] lines in IC 1848 are weak and the S/N in the observed spectra is limited. Moreover, the [O I] 145 µm line is close to a telluric ozone line, which affects the calibration of the data. However, we conclude that the apparent discrepancy between the intensities of the two [O I] lines cannot be explained by instrumental effects.

Based on past experience (Joblin et al. 2018; Wu et al. 2018), the Meudon PDR code may underestimate the [C II] line intensity in astronomical sources. This is explained by [C II] contributions from gas in front or behind the PDR - the model only reproduces the emission associated with the PDR. For oxygen lines, the plane-parallel geometry of the code can sometimes lead to problems, because photons from optically thick lines are artificially trapped in the model, whereas in reality they can escape in directions perpendicular to the line of sight. This line trapping effect is more important for the optically thick 63 µm line, with the optically thin 145 µm line being little affected, and can artificially decrease the [O I] 63/145 µm line ratio. When comparing PDR models to observations, a scaling factor, κ, is sometime used to take into account that (a) several PDR layers along the line of sight can be present within the beam (κ > 1), (b) the PDR is inclined (κ > 1, an effect that we attempted to include directly in the analysis presented here), and (c) the PDR emission does not fill the beam (κ < 1). While the [C II] emission in IC 1848 is spatially extended, beam filling factors for different tracers may vary. Since we have only sparse single-point observations of the [O I] lines with a limited S/N, we cannot ascertain the extent of the [O I] emission, which is typically less extended than [C II] given its higher critical density.

Röllig et al. (2007) presented a comparison of numerical codes developed for modeling the physics and chemistry of UV illuminated regions. Their Figure 15 shows the surface brightness of the main fine-structure cooling lines computed using eight such models. With the exception of the Meijerink code, all models give the [O I] 63/145 µm intensity ratio greater than unity (in main beam brightness temperature units). The [O I] 63/145 µm ratio predicted by the Meudon PDR code is higher compared to other codes (e.g., a factor of ~1.8 higher that the KOSMA-τ model prediction). The particular model used for this comparison has a high UV illumination (χ = 10⁵) and the results may not be representative for all physical conditions. However, the differences among predictions of the various PDR models do not appear to be large enough to account for the discrepancy in the [O I] 63/145 µm intensity ratio between the model predictions and our IC 1848 observations. Given the modeling and observational uncertainties (limited S/N, calibration uncertainty), a factor 2-3 agreement with the observations is usually considered satisfactory. However, our observations suggest that the [O I] 63/145 µm intensity ratio is a sensitive probe of the physical conditions (density or pressure, χ) and chemistry in moderate A_V and UV illumination regions such as IC 1848.

Figure 5 shows the gas temperature, density, and abundance profiles near the cloud surface for Model 1027. PDR chemistry computations show that ionized carbon is only abundant in the outer A_V ~ 1.5 layer, transitioning to neutral carbon, CO, and other carbon-bearing species in deeper layers. Atomic oxygen extends somewhat deeper into the cloud up to A_V ~ 2.5, before being transformed into CO, water, and other oxygen-bearing species. The gas in the [C II] emitting region is molecular, with temperatures between 50 and 150 K.

Fig. 6 Equally weighted average [O I] spectra toward the head of the G1 filament (within the white rectangle in Figure A.6). Averages include [O I] 63 µm spectra at 14 individual positions and [O I] 145 µm spectra at 3 individual positions. The excess noise on the right hand side of the [O I] 145 µm spectrum is due to increased atmospheric opacity produced by a nearby telluric ozone line.

4 Summary and conclusions

Using the upGREAT instrument on SOFIA, we imaged the [C II] 158 µm fine structure line in the bright-rimmed clouds G1–G3 in IC 1848 and obtained pointed observations of the [O I] 63 and 145 µm fine structure lines at several positions to characterize the morphology, velocity field, and the physical conditions in the region. The fine structure line intensities were analyzed using the Meudon PDR code. The main results of this study can be summarized as follows:

1. Velocity-resolved [C II] spectra show evidence of a velocity shift at the head of the brightest G1 filament, possibly caused by radiation pressure from the impinging UV photons.

2. Archival Herschel PACS and SPIRE images imply H₂ column densities in the range 10²¹–10²² cm⁻², corresponding to maximum visual extinction A_V ≃ 10. Assuming that the line-of-sight depth of the filaments is the same as their characteristic width in the plane of the sky, the average H₂ volume density is ≃4500 cm⁻³.

3. The [C II] emission traces ~17% of the H₂ column density derived from the dust SED fits.

4. The [O I] 63 and 145 µm emission is generally weak in this low-χ PDR. The 63 µm line is detected above 4σ level at 6 positions in G1 and 4 positions in G2, while the 145 µm line is detected only at 3 positions in G1. Both lines are detected at a high S/N in average spectra toward the head of the G1 filament (Figure 6), implying a low average [O I] 63/145 µm line intensity (in K kms⁻¹) ratio of ~0.67.

5. PDR models are unable to explain the observed line intensities of the two [O I] fine structure lines in the IC 1848 region. To obtain a good solution, the intensity of the [O I] 145 µm line would have to be a factor of 3 lower compared to the observations. The [O I] lines in IC 1848 are weak and the S/N in the observed spectra is limited. However, the discrepancy between the intensities of the two [O I] lines is not likely caused by instrumental/calibration effects.

The apparent discrepancy between the model intensities of the two [O I] lines and IC 1848 observations may be due to the limited signal-to-noise ratio or trapping of the [O I] 63 µm photons in a plane-parallel geometry model. Nevertheless, our observations suggest that the [O I] 63/145 µm intensity ratio is a sensitive probe of the physical conditions (density or pressure, χ) and chemistry in moderate A_V and UV illumination regions such as IC 1848. Antarctic balloons are the only current FIR facilities that would allow extending this work to other sources. However, some of the proposed NASA FIR Astrophysics Probe concepts include high-resolution spectroscopic capabilities in this wavelength range.

Acknowledgements

Based on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy (SOFIA). SOFIA was jointly operated by the Universities Space Research Association, Inc. (USRA), under NASA contract NAS2-97001, and the Deutsches SOFIA Institut (DSI) under DLR contract 50 OK 0901 to the University of Stuttgart. GREAT was a development of the MPI für Radioastronomie and the KOSMA/Universität zu Köln, in cooperation with the DLR Institut für Optische Sensorsysteme, financed by the participating institutes, by the German Aerospace Center (DLR) under grants 50 OK 1102, 1103 and 1104, and within the Collaborative Research Centre 956, funded by the Deutsche Forschungsgemeinschaft (DFG). Part of 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., P.F.G., and Y.S. acknowledge financial support from the National Aeronautics and Space Administration (NASA) Astrophysics Data Analysis Program (ADAP). We thank Franck Le Petit and William D. Langer and an anonymous referee for helpful comments.

Appendix A G2 and G3 filaments

Figures A.1 and A.2 show moment maps and [C II] spectra at selected positions toward the G2 and G3 filaments, respectively. Figures A.3–A.5 show the channel maps of the [C II] emission in the filaments G1–G3, respectively. Figure A.6 shows [O I] spectra toward the G1 and G2 filaments. Figure A.7 shows the continuum SEDs at selected positions toward the G2 and G3 filaments.

Fig. A.1 [C II] emission toward the G2 filament. (Top) Color image of the [C II] integrated line intensity in the −45 to −37 km s−1 velocity range. The right panel shows spectra at selected positions labeled A–E in the image. Blue circles correspond to the FWHM SOFIA beam size at the [C II frequency. (Bottom) Maps of the [C II] line center velocity and velocity dispersion (color images, left and right panels, respectively). Only pixels with [C II] intensities above 3.5σ (6.3 K km s−1) are shown. Black contours show the integrated line intensity, with the same contour levels as in Figure 1. Black stars mark the location of the bright IRAC compact source.

Fig. A.2 [C II] emission toward the G3 filament. (Top) Color image of the [C II] integrated line intensity in the −44 to −38 km s−1 velocity range. The right panel shows spectra at selected positions labeled A–E in the image. Blue circles correspond to the SOFIA beam at the [C II frequency. (Bottom) Maps of the [C II] line center velocity and velocity dispersion (color images, left and right panels, respectively). Only pixels with [C II] intensities above 3.5σ (8.4 K km s−1) are shown. Black contours show the integrated line intensity, with the same contour levels as in Figure 1. Black stars mark the location of the bright IRAC compact source.

Fig. A.3 Channel maps of the [C II] emission in the G1 filament. Center velocity of each channel is labeled in the upper-right corner.

Fig. A.4 Channel maps of the [C II] emission in the G2 filament. Center velocity of each channel is labeled in the lower-right corner.

Fig. A.5 Channel maps of the [C II] emission in the G3 filament. Center velocity of each channel is labeled in the upper-right corner.

Fig. A.6 [O I] spectra (black histograms) toward positions G–O in the G1 filament (top) and positions F–G in the G2 filament (bottom). Red histograms show the corresponding [C II] spectra divided by 10. Positions at which the spectra are taken are marked on the SPIRE 350 µm images shown in the left panels.

Fig. A.7 Spectral energy distributions toward positions A–E in the G2 and G3 filaments, top and bottom respectively. Blue curves show the modified blackbody fits to the Herschel PACS and SPIRE fluxes from which the dust temperatures and 350 µm optical depths used for the column density determination are derived.

References

All Tables

All Figures

Fig. 1 Spitzer/IRAC 8 µm image of the three filaments in IC 1848 (G1–G3, left-to-right, respectively). Bright IRAC compact sources are clearly visible at the head of the filaments (white saturated spots in the image). Cyan contours show the [C II] 158 µm integrated line intensity observed with SOFIA/upGREAT. Contour levels are 8, 12, 18, and 27 K km s−1, and typical 1σ noise levels are 3.3, 1.8, and 2.4 K km s−1 in G1–G3, respectively. Blue stars mark the location of the three brightest 2MASS/WISE sources at the heads of the filaments, and the white arrow shows the direction toward the O6V ionizing star HD 17505.

In the text

Fig. 2 [C II] emission toward the G1 filament. (Top) Color image of the [C II] integrated line intensity in the −40 to − 33.5 km s−1 velocity range. The right panel shows spectra at selected positions labeled A–F in the image. Blue circles correspond to the FWHM SOFIA beam size at the [C II] frequency. (Bottom) Maps of the [C II] line center velocity and velocity dispersion (color images, left and right panels, respectively). Only pixels with [C II] intensities above 3.5σ (11.6 K km s−1) are shown. Black contours show the integrated line intensity, with the same contour levels as in Figure 1. Black star marks the location of the bright IRAC compact source in each panel.

In the text

	Fig. 3 Spectral energy distributions toward positions A–F in the G1 filament based on Spitzer and Herschel observations. Blue curves show the modified blackbody fits to the Herschel PACS and SPIRE fluxes from which the dust temperatures and 350 µm optical depths used for the column density determination are derived. Typical flux calibration uncertainties are ~5%.
In the text

Fig. 4 χ2 values for a grid of isochoric ISMDB models compared to the observed fine structure line intensities in the G1 filament. Blue circle and the elongated red area show best fits, respectively, to the [C II]/[O I] 145 µm data at positions M–O and [C II]/[O I] 63 µm data at positions G–L. The elongated magenta area shows the best fit to all data with [O I] 145 µm line intensities scaled down by a factor of 3, as described in the text. All models assume normal incidence. Increasing the incidence angle to 60° provides comparable fit, but shifts the best-fit solution to higher densities and lower χ values (magenta area).

In the text

	Fig. 5 Temperature, density, and abundance profiles in the outer PDR layers of Model 1027.
In the text

	Fig. 6 Equally weighted average [O I] spectra toward the head of the G1 filament (within the white rectangle in Figure A.6). Averages include [O I] 63 µm spectra at 14 individual positions and [O I] 145 µm spectra at 3 individual positions. The excess noise on the right hand side of the [O I] 145 µm spectrum is due to increased atmospheric opacity produced by a nearby telluric ozone line.
In the text

Fig. A.1 [C II] emission toward the G2 filament. (Top) Color image of the [C II] integrated line intensity in the −45 to −37 km s−1 velocity range. The right panel shows spectra at selected positions labeled A–E in the image. Blue circles correspond to the FWHM SOFIA beam size at the [C II frequency. (Bottom) Maps of the [C II] line center velocity and velocity dispersion (color images, left and right panels, respectively). Only pixels with [C II] intensities above 3.5σ (6.3 K km s−1) are shown. Black contours show the integrated line intensity, with the same contour levels as in Figure 1. Black stars mark the location of the bright IRAC compact source.

In the text

Fig. A.2 [C II] emission toward the G3 filament. (Top) Color image of the [C II] integrated line intensity in the −44 to −38 km s−1 velocity range. The right panel shows spectra at selected positions labeled A–E in the image. Blue circles correspond to the SOFIA beam at the [C II frequency. (Bottom) Maps of the [C II] line center velocity and velocity dispersion (color images, left and right panels, respectively). Only pixels with [C II] intensities above 3.5σ (8.4 K km s−1) are shown. Black contours show the integrated line intensity, with the same contour levels as in Figure 1. Black stars mark the location of the bright IRAC compact source.

In the text

	Fig. A.3 Channel maps of the [C II] emission in the G1 filament. Center velocity of each channel is labeled in the upper-right corner.
In the text

	Fig. A.4 Channel maps of the [C II] emission in the G2 filament. Center velocity of each channel is labeled in the lower-right corner.
In the text

	Fig. A.5 Channel maps of the [C II] emission in the G3 filament. Center velocity of each channel is labeled in the upper-right corner.
In the text

	Fig. A.6 [O I] spectra (black histograms) toward positions G–O in the G1 filament (top) and positions F–G in the G2 filament (bottom). Red histograms show the corresponding [C II] spectra divided by 10. Positions at which the spectra are taken are marked on the SPIRE 350 µm images shown in the left panels.
In the text

	Fig. A.7 Spectral energy distributions toward positions A–E in the G2 and G3 filaments, top and bottom respectively. Blue curves show the modified blackbody fits to the Herschel PACS and SPIRE fluxes from which the dust temperatures and 350 µm optical depths used for the column density determination are derived.
In the text