# JWST Observations of Cosmic-Ray-excited H_2 in Barnard 68: Spatial Variations and Constraints on Cosmic-Ray Attenuation

David A Neufeld, Kedron Silsbee, Alexei V Ivlev, Shmuel Bialy, Brandt A.L Gaches, Marco Padovani, Sirio Belli, Thomas G Bisbas, Amit Chemke,

Benjamin Godard, et al.

<figure><img src="image_1.png" /></figure>

► To cite this version:

David A Neufeld, Kedron Silsbee, Alexei V Ivlev, Shmuel Bialy, Brandt A.L Gaches, et al. JWST Observations of Cosmic-Ray-excited H_2 in Barnard 68: Spatial Variations and Constraints on Cosmic-Ray Attenuation. The Astrophysical Journal, 2026, 998 (1), pp.71-10.3847/1538-4357/ae22cd hal-05495445

HAL Id: hal-05495445

https://hal.science/hal-05495445v1

Submitted on 20 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
OPEN ACCESS

<figure><img src="image_3.png" /></figure>

CrossMark

# JWST Observations of Cosmic-Ray-excited H₂ in Barnard 68: Spatial Variations and Constraints on Cosmic-Ray Attenuation

David A. Neufeld¹, Kedron Silsbee², Alexei V. Ivlev³, Shmuel Bialy⁴, Brandt A. L. Gaches⁵, Marco Padovani⁶, Sirio Belli⁷, Thomas G. Bisbas⁸, Amit Chemke⁴, Benjamin Godard⁹,¹⁰, James Muzerolle Page¹¹, and Christian Rab³,¹²

¹William H. Miller Department of Physics & Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA

²Physics Department, University of Texas at El Paso, El Paso 79968, USA

³Max-Planck-Institut für Extraterrestrische Physik, Giessenbachstrasse 1, 85748 Garching, Germany

⁴Technion—Israel Institute of Technology, Haifa, 3200003, Israel

⁵Faculty of Physics, University of Duisburg-Essen, Lotharstraße 1, 47057 Duisburg, Germany

⁶INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy

⁷Dipartimento di Fisica e Astronomia, Università di Bologna, 40125 Bologna, Italy

⁸Research Center for Astronomical Computing, Zhejiang Laboratory, Hangzhou, 311000, People's Republic of China

⁹Observatoire de Paris, Université PSL, Sorbonne Université, LERMA, 75014 Paris, France

¹⁰Laboratoire de Physique de l'Ecole Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université de Paris, 75005 Paris, France

¹¹Space Telescope Science Institute, Baltimore, MD 21218, USA

¹²University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität München, Scheinerstr. 1, 81679 Munich, Germany

Received 2025 October 10; revised 2025 November 17; accepted 2025 November 19; published 2026 February 3

## Abstract

We present James Webb Space Telescope (JWST) NIRSpec observations of the starless dark cloud Barnard 68 (B68) that reveal the spatially resolved signature of cosmic-ray-excited molecular hydrogen (CRXH₂) emissions for the first time. Following up on our initial detection of CRXH₂ emissions from B68, we now exploit JWST's sensitivity and spatial multiplexing to map CRXH₂ rovibrational lines across 16 sight lines through the cloud. By disentangling the CRXH₂ and UV-pumped H₂ components, we isolate the para-H₂-dominated spectrum attributable to cosmic-ray (CR) excitation. We find that there are significant spatial variations in the ratio of the CRXH₂ line intensity to the line-of-sight H₂ column density; these cannot be accounted for by dust extinction alone and demonstrate a clear attenuation of the CR flux with increasing shielding column. Modeling B68 as a Bonnor–Ebert sphere, we constrain both the unshielded CR ionization rate, ζH₂, and how it decreases with shielding column. At a reference depth of N(H₂) = 3 × 10²¹ cm⁻², we infer ζH₂ ≈ 1.4 × 10⁻¹⁶ s⁻¹, a factor of ≈3 higher than the average value derived from H⁺ absorption studies. These results provide the most direct probe to date of CR penetration into cold, dense gas, offering new constraints on both the microphysics of CR–H₂ interactions and the attenuation of low-energy CRs in molecular clouds. Our findings establish CRXH₂ emission as a powerful new diagnostic of the CR environment in interstellar space.

Unified Astronomy Thesaurus concepts: Cosmic rays (329); Molecular clouds (1072); Interstellar medium (847); Dark interstellar clouds (352); Infrared spectroscopy (2285); Interstellar line emission (844)

Materials only available in the online version of record: figure set

## 1. Introduction

Recent JWST observations of the starless dark cloud Barnard 68 (B68) have led to the first detection of cosmic-ray-excited H₂ (CRXH₂) emissions from interstellar space (S. Bialy et al. 2025; hereafter B25), thereby opening a new window on Galactic cosmic rays (CRs). Because CRXH₂ emissions originate within cold (~15 K) gas in the cloud interior, they exhibit a spectrum that is readily distinguished from UV-excited hydrogen (UVXH₂) emissions from warmer (~100 K) gas at the cloud surface: CRXH₂ exhibits a spectrum that is strongly dominated by para-H₂ (even- J) lines (S. Bialy 2020, hereafter B20; M. Padovani et al. 2022, hereafter P22) whereas UVXH₂ emits lines of both ortho- and para-H₂ with an ortho-to-para ratio (OPR) ≈ 1.6 (e.g., A. Sternberg 1988, hereafter S88). The spectrum of B68 measured with the NIRSpec multi-object spectrograph (MOS) on JWST showed (B25) a clear excess of CRXH₂ lines, implying an average cosmic-ray ionization rate (CRIR), ζH₂, of

1.7 × 10⁻¹⁶ s⁻¹. This value is the total rate of ionization of H₂, including secondary ionizations. B25 also observed a nearby “background” region, where the column density is much lower than at the B68 position. Because the CR signal scales with the column density (B20), we expect the background region to be dominated by UVXH₂ emissions, and indeed, the relative H₂ line intensities measured at the background position were fully consistent with the UV model predictions.

In the present study, we have extended our analysis to include spatial variations in the CRXH₂ emissions. The exquisite sensitivity of JWST/NIRSpec is sufficient to permit high signal-to-noise ratio (SNR) detections of H₂ lines at multiple positions across B68, allowing variations in ζH₂ to be inferred. These, in turn, provide important constraints on the propagation of CR from the diffuse outer regions of the B68 cloud to the dense interior.

In Section 2, we describe the data acquisition, reduction, and analysis methods that were adopted. In Section 3, we present the line intensities measured as a function of position on the sky. In Section 4, we compare the observed spatial variations in the line intensities with a parameterized model for the attenuation of CR that also includes the effects of dust extinction. A discussion follows in Section 5, and the conclusions of our study are summarized in Section 6.

<figure><img src="image_4.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Arrangement of NIRSpec MOS shutters across B68. White lines show the locus of the 382 open shutters used to obtain spectra across the cloud. The 16 averaged slitlets (~11" each) are marked with white dots. The background ("OFF") position, observed 30' north of B68, used the same configuration. The extinction map is adapted from J. Alves et al. (2025, private communication).</figcaption></figure>

## 2. Observations and Data Reduction

As discussed in B25, all the data were acquired in JWST program 5064, using the NIRSpec MOS with the G235H/F170LP grating/filter combination. The latter provides almost complete coverage of the 1.66–3.05 µm spectral region at a spectral resolving power, $\lambda/\Delta\lambda$ of 2700, sufficient to resolve individual rotational lines within the H$_2$ $v = 1-0$ Q-branch around 2.4 µm. The MOS allows spectra to be obtained in separate shutters that may be opened at commandable positions within a 3"/4 × 3"/6 region on the sky. For the present study, we arranged for 382 open shutters to lie on a linear locus across the B68 cloud. Because the four quadrants of the detector do not perfectly abut one another, there is a central gap in the spatial coverage of length 37" and a small (~0.1 µm) gap in the wavelength coverage around 2.25 µm. In Figure 1, we show the locus of the shutters (white lines) superposed on an extinction map of the B68 cloud (J. Alves et al. 2025, private communication). In a separate observation, we used the same arrangement of shutters to observe a background region (hereafter the "OFF" position) located 30' due north of B68.

As discussed in B25, all the observations were conducted using the NRSIRS2 readout mode, with 4 integrations of 20 groups each at five dither positions on B68, providing a total on-source time of 8.2 hr. The dither pattern shifted the pointing by five shutters along both spatial axes to reduce contamination from nearby stars and to average over detector defects. Measurements of the background region employed three dither positions, each consisting of a single 20-group integration, for a total exposure time of 1.2 hr.

All data were reprocessed beginning with the raw "uncal" files using the JWST calibration pipeline (v1.16.1) together with CRDS context 1303. Standard settings were applied throughout, except during outlier rejection, where weighting by "exptime" substantially reduced the number of residual outliers relative to the default configuration. The multi shutter array (MSA) metadata were adjusted so that the pipeline could align and extract spectra from each open shutter of the long slit across all dither positions in both target and background fields.

Average spectra were extracted separately for each of 16 slitlets of length ~11", each formed from ~24 adjacent shutters on B68; the central positions for these 16 slitlets are marked by white dots. We also obtained an average spectrum for the entire ensemble of slitlets, both for slitlets on B68 and for those on the OFF position. The difference spectrum, "B68 – OFF," was presented in B25 (their Figure 1(b)). A total of 21 spectral lines were securely detected in the ensemble average spectrum for B68. These comprise (i) 11 lines in the $v = 1-0$ band of H$_2$—one of which, $v = 1-0$ S(0), lies in the detector gap near 2.25 µm and can be observed only in a subset of the slitlets and was therefore excluded from the analysis; (ii) five lines in the $v = 2-1$ band of H$_2$; (iii) four lines in the $v = 3-2$ band of H$_2$; (iv) the Paschen alpha ($n = 4-3$) line of atomic hydrogen.

## 3. Results

### 3.1. Average H$_2$ Line Intensities

We obtained Gaussian fits to each of these lines, using the Levenberg-Marquardt algorithm to fit a flat baseline plus a Gaussian to each line. Here, there were four free parameters:
the line width (FWHM in kilometers per second), the integrated line intensity (erg cm⁻² s⁻¹ sr⁻¹), the line centroid (kilometers per second in the barycentric frame), and the continuum level (MJy sr⁻¹). Each fit was optimized over the barycentric velocity interval [−600, 600] km s⁻¹. In computing the fits to the v = 1–0 O(2) line, we corrected for contamination by the partially blended H Br-β line, which has a rest wavelength that is only 1.015 × 10⁻³ µm smaller than that of v = 1–0 O(2), corresponding to a blueshift of 116 km s⁻¹. This correction was made by assuming a Br-β line width and a barycentric centroid velocity identical to that observed for H Pa-α and a Br-β/Pa-α flux ratio of 0.134, the value expected for Case B recombination (P. J. Storey & D. G. Hummer 1995).

Figure 2 shows the line profiles (blue histogram) and fitted Gaussians (red) for 9 example lines, and Table 1 lists the fit parameters with 1σ statistical errors. Figure 3 shows the line FWHM (kilometers per second) and centroid velocities (kilometers per second in the barycentric frame) for the Gaussian fits shown in Figure 2 as a function of line wavelength. Here, para-H₂ lines appear in red and ortho-H₂ lines appear in blue. The dashed black line in the upper panel is a quadratic fit to the line FWHM. It is in good agreement with what is expected for an unresolved spectral line (P. Jakobsen et al. 2022, their Figure 5). The horizontal lines in the lower panel show the mean velocity centroids (−15.9, −10.2, −14.4 km s⁻¹, respectively), for the ortho-H₂ lines (blue), para-H₂ lines (red), and for all H₂ lines (black). There is no significant difference between the average centroid velocities for the ortho- and para-lines, and the line-to-line variation of velocity centroids is consistent with the expected systematic calibration errors<fnref n="13" /> of 15 km s⁻¹. Within the expected errors, the velocity centroids also agree with the systemic velocity of B68 determined by J. F. Alves et al. (2001; hereafter A01): +3.4 km s⁻¹ relative to the local standard of rest, corresponding to a barycentric velocity of −8.2 km s⁻¹.

Because the results shown in Figure 3 support our strong expectations that (1) the lines are spectrally unresolved with NIRSpec, and (2) the velocity centroids are the same for all lines, we decided to adopt these expectations as priors in the fitting process and subsequently refit all lines with the velocity centroid fixed at −14.4 km s⁻¹ and the line width constrained to the (wavelength-dependent) values shown by the dashed line in the upper panel of Figure 3. Only the continuum level and line intensity are now allowed to vary in the optimization of the fit. This approach reduces the uncertainties in the line intensities and is particularly advantageous for weak lines. For strong lines, no significant effects on the inferred line parameters could be discerned. The resultant line intensities are tabulated in Table 2 for the B68, OFF, and B68 − OFF spectra.

## 3.2. Separation of the UVXH₂ and CRXH₂ Emissions

Figure 4 (upper panel) shows the integrated line intensities for the OFF position, I(OFF) = 10⁻⁷I_7(OFF) erg cm⁻² s⁻¹ sr⁻¹, each divided by fUV, defined here as the fraction of the rovibrational UVXH₂ line emission that is predicted to emerge in each line. Blue and red points are for transitions of ortho- and para-H₂, respectively. Here, we used the values for fUV computed by S88 for a gas temperature of 100 K, and the blue and red

horizontal lines show the mean values of I_7(OFF)/fUV for ortho- and para-H₂. Figure 4 indicates that the OFF position line strengths are in good agreement with the predictions for pure UV fluorescence in gas with a temperature 100 K (in which the OPR is ≈1.6). Our estimate for the summed intensity for all H₂ rotational lines (both observed and unobserved) is (30 ± 2) × 10⁻⁷ erg cm⁻² s⁻¹ sr⁻¹.

Because CRXH₂ emissions originate primarily in the cold cloud interior, where the OPR is ≪1, they make a negligible contribution to the ortho-H₂ line strengths. For the difference spectrum, B68 − OFF, the ortho-H₂ lines are weaker than in either the B68 or OFF spectra, but the cancellation is incomplete; this presumably reflects differences in the incident UV field and/or the geometry (e.g., the degree of limb brightening). The overall ortho-H₂ line intensities in the difference spectrum are smaller than those in the OFF spectrum: on average, the ratio, Q, of the ortho-H₂ line intensities in the difference spectrum to those in the OFF spectrum is Q = 0.445. We may therefore estimate the CRXH₂ line emissions from B68 by subtracting 0.445 times the OFF line intensity from the intensity observed in the difference spectrum. Here, we assume that the UVXH₂ line ratios are the same at the two observed positions, B68 and OFF. The resultant intensities, attributable to CR excitation alone, are shown in the bottom panel of Figure 4. Unequivocal detections (with an SNR > 5) are obtained for CRXH₂ emission in the v = 1–0 O(2), O(4), Q(2), and S(2) lines. CRXH₂ emissions are also securely detected in the v = 1–0 Q(4) line and are tentatively detected (4 < SNR < 5) in the v = 2–1 O(2) and Q(2) lines. The strongest line, v = 1–0 O(2), has a CR contribution that is detected in the difference spectrum at the 28σ level.

## 3.3. Spatial Variation

To probe the spatial variation of the CRXH₂ emissions, we obtained line intensities for each of the 16 individual slitlets, corresponding to 16 positions along the slit (Figure 1). Here, we adopted the method used above to compute the intensities plotted in the lower panel of Figure 4, but we now implemented it separately for each slitlet (with account taken for spatial variations in the value of Q). In Figure 5, results are shown for the three strongest CRXH₂ lines. The intensities attributable to CRXH₂, after subtraction of the UVXH₂ contribution, are shown as a function of the position along the slit.

In Figure 6, we show the ratio of these intensities to the column density of H₂ along the line of sight,

$$N(\text{H}_2) = 10^{22} N_{22} \text{ cm}^{-2}.$$

The latter was obtained from the extinction map, assuming that the gas is fully molecular with N(H₂) = 1/2 N<sub>H</sub> and adopting an N<sub>22</sub>/A<sub>V</sub> ratio of 0.094.<fnref n="14" /> The column densities determined from

<footnote n="14">¹⁴ This value corresponds to N<sub>H</sub>/A<sub>V</sub> = 1.87 × 10<sup>21</sup> cm<sup>-2</sup>, the canonical N<sub>H</sub>/A<sub>V</sub> ratio for the diffuse ISM. This ratio is based on the N<sub>H</sub>/E(B-V) ratio obtained from UV absorption measurements of H and H₂ (R. C. Bohlin et al. 1978) and an assumed A<sub>V</sub>/E(B-V) ratio of 3.1. X-ray absorption observations (e.g., H. Zhu et al. 2017 and references therein) yield similar values for N<sub>H</sub>/A<sub>V</sub> in the diffuse ISM. For the dense ISM, some studies have favored a smaller value of N<sub>H</sub>/A<sub>V</sub>, guided by dust models from B. T. Draine (2003), N. J. Evans et al. (2009) assumed N<sub>H</sub>/A<sub>V</sub> = 1.37 × 10<sup>21</sup> cm<sup>-2</sup>, about 27% smaller than the value we adopt here. On the other hand, direct measurements of H₂ in the Taurus molecular cloud, obtained by observing its weak near-IR quadrupole transitions in absorption (J, H. Lacy et al. 2017), yielded a mean N(H₂)/E(J-K) ratio of 5.4 × 10<sup>21</sup> cm<sup>-2</sup>, from which a N(H₂)/A<sub>V</sub> ratio of 1.0 × 10<sup>21</sup> cm<sup>-2</sup> was inferred. This, then, is slightly (6%) larger than the ratio we assume.</footnote>

<footnote n="13">¹³ https://jwst-docs.stsci.edu/jwst-calibration-status/nirspec-calibration-status, downloaded 2025 May 2.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Example H₂ rovibrational line profiles from B68 (black histograms) with Gaussian fits (red curves).</figcaption></figure>
<table><caption>**Table 1**
Gaussian Fit Results for Detected H₂ and H I Lines in B68</caption><thead><tr><th>Transition</th><th>λ<sub>rest</sub><br>(μm)</th><th>Intensity<br>(10<sup>-7</sup> erg cm<sup>-2</sup> s<sup>-1</sup> sr<sup>-1</sup>)</th><th>v<sub>cent</sub><br>(km s<sup>-1</sup>)</th><th>FWHM<br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Pa-α</td><td>1.87561</td><td>2.90 ± 0.076</td><td>-20 ± 2</td><td>147 ± 4</td></tr><tr><td>1-0 S(3)</td><td>1.95756</td><td>0.59 ± 0.040</td><td>-28 ± 5</td><td>168 ± 12</td></tr><tr><td>1-0 S(2)</td><td>2.03376</td><td>0.59 ± 0.056</td><td>-21 ± 5</td><td>134 ± 13</td></tr><tr><td>1-0 S(1)</td><td>2.12183</td><td>0.97 ± 0.029</td><td>-18 ± 2</td><td>139 ± 4</td></tr><tr><td>1-0 Q(1)</td><td>2.40659</td><td>0.95 ± 0.036</td><td>-12 ± 2</td><td>111 ± 4</td></tr><tr><td>1-0 Q(2)</td><td>2.41344</td><td>0.92 ± 0.082</td><td>-21 ± 4</td><td>112 ± 11</td></tr><tr><td>1-0 Q(3)</td><td>2.42373</td><td>0.66 ± 0.030</td><td>-11 ± 3</td><td>133 ± 6</td></tr><tr><td>1-0 Q(4)</td><td>2.43749</td><td>0.39 ± 0.029</td><td>-6 ± 4</td><td>125 ± 10</td></tr><tr><td>3-2 S(0)</td><td>2.50144</td><td>0.33 ± 0.022</td><td>9 ± 4</td><td>143 ± 10</td></tr><tr><td>2-1 Q(1)</td><td>2.55099</td><td>0.49 ± 0.020</td><td>-25 ± 2</td><td>133 ± 6</td></tr><tr><td>2-1 Q(2)</td><td>2.55851</td><td>0.43 ± 0.021</td><td>-12 ± 3</td><td>127 ± 6</td></tr><tr><td>2-1 Q(3)</td><td>2.56983</td><td>0.20 ± 0.020</td><td>-10 ± 5</td><td>99 ± 11</td></tr><tr><td>1-0 O(2)</td><td>2.62688</td><td>1.55 ± 0.095</td><td>-20 ± 4</td><td>135 ± 9</td></tr><tr><td>3-2 Q(1)</td><td>2.71025</td><td>0.19 ± 0.038</td><td>-18 ± 9</td><td>104 ± 23</td></tr><tr><td>3-2 Q(2)</td><td>2.71862</td><td>0.15 ± 0.015</td><td>-11 ± 4</td><td>90 ± 10</td></tr><tr><td>2-1 O(2)</td><td>2.78616</td><td>0.29 ± 0.032</td><td>-5 ± 4</td><td>87 ± 10</td></tr><tr><td>1-0 O(3)</td><td>2.80252</td><td>0.73 ± 0.045</td><td>-12 ± 3</td><td>102 ± 7</td></tr><tr><td>3-2 O(2)</td><td>2.96206</td><td>0.14 ± 0.018</td><td>-8 ± 5</td><td>81 ± 11</td></tr><tr><td>2-1 O(3)</td><td>2.97406</td><td>0.28 ± 0.023</td><td>-11 ± 4</td><td>99 ± 9</td></tr><tr><td>1-0 O(4)</td><td>3.00387</td><td>0.53 ± 0.021</td><td>-10 ± 2</td><td>87 ± 4</td></tr></tbody></table>

Note. Listed are the transition, rest wavelength, integrated line intensity, velocity centroid, and FWHM.

the extinction map are in acceptable agreement with independent determinations based on submillimeter dust emission (M. Nielb-ock et al. 2012; A. Roy et al. 2014), although the latter are also dependent on the assumed dust properties (i.e., the submillimeter opacity law).

The decline of $I_{-7}/N_{22}$ with $N_{22}$ has two possible causes that must be disentangled: (1) a decline in $\zeta_{H_2}$ with increasing depth into the cloud and (2) dust extinction of the emergent $H_2$ emissions. The next section describes a method for separating these two effects.

## 4. Cloud Model and Inferred Attenuation of CR

### 4.1. Cloud Model

To analyze the spatial variation of the observed CRXH₂ emission, we used a simple model in which the B68 cloud was assumed to be a fully molecular Bonnor–Ebert sphere. This assumption is motivated by the pioneering study of A01, who analyzed an extinction map they obtained for B68. They found that a Bonnor–Ebert sphere with a dimensionless outer radius $\xi_{out} = 6.9$ provided an excellent fit to the azimuthally averaged extinction as a function of projected angular distance. The dimensionless radius, $\xi$, is related to the physical radius, $r$, the isothermal sound speed, $c_s$, and the central $H_2$ density $n_0$, by the expression

$$ \xi = \frac{r}{c_s} \sqrt{4\pi G\mu n_0}, $$

where G is the gravitational constant and $\mu = 4.58 \times 10^{-24}$ g is the mean mass per $H_2$ molecule (including helium). The best-fit value for $\xi_{out}$ obtained by A01 was close to—but somewhat larger than—the critical value of 6.45 above which a (nonmagnetized) cloud becomes unstable. The volume density, $n_{H_2}$, was obtained as a function of the distance from

the center of the cloud, $r$, by solving the Lane–Emden equation for an isothermal sphere in hydrostatic equilibrium.

In the following analysis, we adopt a Cartesian coordinate system with the x-axis in the plane of the sky along the slit, the z-axis along the line of sight with z values increasing with distance from the observer, and the y-axis perpendicular to x and z. The origin is at the center of the cloud, and the slit is assumed to lie at $y = 0$. The radial coordinate, $r$, is then $(x^2 + y^2 + z^2)^{1/2}$. The total column density along the line of sight, $N_{H_2}$, was obtained as a function of projected distance, $x$, by integrating the volume density

$$ N_{H_2}(x) = \int n_{H_2}(r) dz = \int n_{H_2}([x^2 + z^2]^{1/2}) dz. $$

In Figure 7, the results for $n_{H_2}(r)$ (blue curve) and $N_{H_2}(x)$ (red curve) are shown as a function of $r/r_{out}$ or $x/r_{out}$, where $r_{out}$ is the outer radius of the cloud.

The quantities plotted in Figure 7 are dimensionless. To completely specify the properties of the B68 cloud, two of the following four dimensional quantities must be specified: the outer radius, $r_{out}$, the maximum $H_2$ column density, $N(0)$, the central density, $n(0)$, or the sound speed, $c_s$. The first two quantities are closest to the observations, and we adopt $r_{out} = 1.25 \times 10^4$ au $= 1.87 \times 10^{17}$ cm and $N(0) = 2.82 \times 10^{22}$ cm$^{-2}$, following A01. These values imply a central density of $1.92 \times 10^5$ cm$^{-3}$ and an isothermal sound speed, $c_s$, of 0.233 km s$^{-1}$; the latter corresponds to a gas temperature of 18 K if thermal pressure is dominant.<fnref n="15" />

The $H_2$ emission from such a sphere was computed as a function of $x$, with the inclusion of CR attenuation and dust extinction. We treated the effects of CR attenuation by adopting a simple parameterization for the reduction in the CRIR due to CR propagating through a shielding $H_2$ column density of $N_s$:

$$ \zeta_{H_2} = \frac{\zeta_0}{(1 + N_s/N_0)^\alpha} = \zeta_0 F_s(N_s), $$

where $\zeta_0$, $N_0$, and $\alpha$ are parameters that will be adjusted to optimize the fit to the data. We note here that the shielding column, $N_s$, is different from the column density along a given sight line; the former varies with position along the sight line and the CR propagation direction. The power-law index $\alpha$ describes the rate at which $\zeta_{H_2}$ decreases with shielding column in the limit of large $N_s$, while the parameter $N_0$ is the shielding column below which $\zeta_{H_2}$ is constant. If the intrinsic energy spectrum of the CR at the site of their production is a power law extending to low energies, $N_0$ is naturally interpreted as the column density through which the CR have passed in propagating from their production site to the surface of the B68 molecular cloud.

For the transition from state u to l, the emergent intensity of CRXH₂ at projected distance $x < r_{out}$ from the cloud center is

$$ I_{ul}(x) = \frac{b_u \alpha_{ul} h c}{4\pi \lambda_{ul}} \int_{-z_0}^{+z_0} \zeta_{H_2}(x, z') n_{H_2}([x^2 + z'^2]^{1/2}) e^{-\tau_d} dz' $$

<footnote n="15">¹⁵ This value significantly exceeds the temperature estimate of 10 K obtained by S. Hotzel et al. (2002) from $NH_3$ line ratios. This discrepancy might be resolved if the source distance is considerably smaller than that assumed by A01 or if nonthermal support contributes significantly: while the $NH_3$ line widths measured by Hotzel et al. rule out turbulence as a significant source of nonthermal support, magnetic fields may be important here (R. Kandori et al. 2020).</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Measured line widths (upper panel) and centroid velocities (lower panel) as a function of wavelength for all detected H₂ lines. Para-H₂ lines are shown in red, and ortho-H₂ lines in blue. The dashed curve in the upper panel is a quadratic fit to the unresolved instrumental profile. In the lower panel, horizontal lines show mean velocity centroids for ortho-, para-, and all H₂ lines; the scatter is consistent with calibration uncertainties.</figcaption></figure>

where z is the distance along the line of sight from the cloud midplane; $z_0 = (r_{out}^2 - x^2)^{1/2}$ is the half thickness of the cloud at projected distance, x;

$$ \tau_d(x, z) = \sigma_d \int_{-z_0}^{z} n_{H_2} dz' $$

is the extinction optical depth from the cloud surface to position z along the line of sight; and $\sigma_d$ is the dust extinction cross-section per H₂ molecule at the wavelength, $\lambda_{ul}$, of the H₂

line. To relate $A_\lambda$ to the visual extinction, $A_V$, we adopted the KP5 extinction curve tabulated in K. M. Pontoppidan et al. (2024): this extinction curve is considered most appropriate to dark clouds such as B68. For simplicity, our treatment adopts a “foreground screen” geometry and thus does not include the diffuse scattered CRXH₂ component. The quantities $b_u$ and $\alpha_{ul}$ are specific to each transition (B25; see their Table S1): they are respectively the excitation rate to state u divided by $\zeta_{H_2}$ and the probability that radiative decay of state u leads to state l.
<table><caption>**Table 2**
Average Line Intensities (10⁻⁷ erg cm⁻² s⁻¹ sr⁻¹) with Line Kinematics Fixed</caption><thead><tr><th>Transition</th><th>λ<br>(μm)</th><th>B68</th><th>OFF</th><th>B68 - OFF</th></tr></thead><tbody><tr><td>Pa-α</td><td>1.87561</td><td>2.976 ± 0.083</td><td>1.905 ± 0.138</td><td>1.071 ± 0.103</td></tr><tr><td>1-0 S(3)</td><td>1.95756</td><td>0.534 ± 0.035</td><td>0.401 ± 0.037</td><td>0.133 ± 0.031</td></tr><tr><td>1-0 S(2)</td><td>2.03376</td><td>0.615 ± 0.044</td><td>0.226 ± 0.026</td><td>0.389 ± 0.047</td></tr><tr><td>1-0 S(1)</td><td>2.12183</td><td>0.963 ± 0.024</td><td>0.623 ± 0.038</td><td>0.340 ± 0.038</td></tr><tr><td>1-0 Q(1)</td><td>2.40659</td><td>0.993 ± 0.031</td><td>0.600 ± 0.040</td><td>0.393 ± 0.034</td></tr><tr><td>1-0 Q(2)</td><td>2.41344</td><td>0.943 ± 0.066</td><td>0.293 ± 0.058</td><td>0.650 ± 0.028</td></tr><tr><td>1-0 Q(3)</td><td>2.42373</td><td>0.613 ± 0.023</td><td>0.448 ± 0.049</td><td>0.166 ± 0.043</td></tr><tr><td>1-0 Q(4)</td><td>2.43749</td><td>0.375 ± 0.023</td><td>0.151 ± 0.023</td><td>0.224 ± 0.027</td></tr><tr><td>3-2 S(0)</td><td>2.50144</td><td>0.269 ± 0.025</td><td>0.143 ± 0.029</td><td>0.127 ± 0.025</td></tr><tr><td>2-1 Q(1)</td><td>2.55099</td><td>0.430 ± 0.022</td><td>0.300 ± 0.035</td><td>0.130 ± 0.031</td></tr><tr><td>2-1 Q(2)</td><td>2.55851</td><td>0.392 ± 0.016</td><td>0.168 ± 0.031</td><td>0.224 ± 0.031</td></tr><tr><td>2-1 Q(3)</td><td>2.56983</td><td>0.214 ± 0.017</td><td>0.132 ± 0.028</td><td>0.083 ± 0.028</td></tr><tr><td>1-0 O(2)</td><td>2.62688</td><td>1.280 ± 0.036</td><td>0.290 ± 0.034</td><td>0.990 ± 0.027</td></tr><tr><td>3-2 Q(1)</td><td>2.71025</td><td>0.184 ± 0.028</td><td>0.169 ± 0.034</td><td>0.015 ± 0.027</td></tr><tr><td>3-2 Q(2)</td><td>2.71862</td><td>0.161 ± 0.012</td><td>0.100 ± 0.023</td><td>0.061 ± 0.024</td></tr><tr><td>2-1 O(2)</td><td>2.78616</td><td>0.302 ± 0.028</td><td>0.107 ± 0.033</td><td>0.195 ± 0.028</td></tr><tr><td>1-0 O(3)</td><td>2.80252</td><td>0.701 ± 0.034</td><td>0.545 ± 0.033</td><td>0.156 ± 0.031</td></tr><tr><td>3-2 O(2)</td><td>2.96206</td><td>0.146 ± 0.014</td><td>0.155 ± 0.020</td><td>-0.010 ± 0.019</td></tr><tr><td>2-1 O(3)</td><td>2.97406</td><td>0.254 ± 0.017</td><td>0.215 ± 0.037</td><td>0.039 ± 0.035</td></tr><tr><td>1-0 O(4)</td><td>3.00387</td><td>0.512 ± 0.019</td><td>0.208 ± 0.032</td><td>0.304 ± 0.036</td></tr></tbody></table>

The three-dimensional dependence of $\zeta_{H_2}$ depends on the direction in which the CR travels. We considered four cases. In cases X, Y, and Z, CR are assumed to propagate along the x-, y-, and z-axes, respectively, with equal fluxes in the ± directions. In Case Q (quasi-isotropic), we assumed equal contributions for all six directions, ±x, ±y, and ±z. Expressions for $\zeta_{H_2}$ are given in the Appendix.

## 4.2. CRIR and Attenuation of CR

For each of the four cases described above—X, Y, Z, and Q—we used the spatial variation of the observed intensity in the 1–0 O(2) line to constrain the adjustable parameters $\zeta_0$, $N_0$, and $\alpha$. Since our observations primarily probe the shielding column density range from 3 to $10 \times 10^{21} \text{ cm}^{-2}$, we adopt a reference column density of $N_{\text{ref}} = 3 \times 10^{21} \text{ cm}^{-2}$ and introduce the parameter $\zeta_{\text{ref}}$, defined as the value of $\zeta_{H_2}$ for a shielding column density of $N_{\text{ref}}$:

$$ \zeta_{\text{ref}} = \frac{\zeta_0}{(1 + N_{\text{ref}}/N_0)^\alpha} $$

Using $\chi^2$ as a measure of goodness of fit, and adopting flat priors on $\log_{10}(\zeta_{\text{ref}})$, $\log_{10}(N_0)$, and $\alpha$ in the intervals [-17, -13], [19, 22], and [0, 2] respectively, we obtained posterior probability distributions for the adjustable parameters. These are represented in Figure Set 8, which also shows (panel (a)) the best fit to $I_{-7}/N_{22}$ for the strongest CRXH$_2$ line ($v = 1-0$ O(2)). The data points here apply to the CRXH$_2$ component alone, the UVXH$_2$ having been subtracted as described in Section 3.

The posterior probability distributions for $\log_{10} N_0$, $\alpha$, and $\log_{10}(\zeta_{\text{ref}})$ are plotted in panels (c), (d), and (e), respectively. The joint probability distribution for $\log_{10} N_0$ and $\log_{10}(\zeta_{\text{ref}})$ is shown in panel (b), where red, green, and blue shaded regions indicate the $1\sigma$, $2\sigma$, and $3\sigma$ error “ellipses.” Clearly, there is a considerable degeneracy in $\log_{10} N_0$ and $\alpha$, as expected, since increasing $\alpha$ or decreasing $\log_{10} N_0$ both increase the degree of attenuation. In Figure Set 8, Cases X, Y, Z, and Q appear, respectively, as Figures 8.1 (shown in the printed version), 8.2,

8.3, and 8.4. The best-fit parameters are shown by a black dot in panel (b) and are tabulated in Table 3.

Given the priors discussed above and the assumption that $\zeta_{H_2}(N_s)$ is proportional to $(1 + N_s/N_0)^{-\alpha}$, we may compute the posterior probability distribution for $\zeta_{H_2}(N_s)$ as a function of $N_s$. Results are shown in Figure 9, where the black curve shows the median value of $\zeta_{H_2}$ (i.e., the value, V, for which the posterior probability of $\zeta_{H_2} < V$ equals 0.5). The dashed and dotted curves show, respectively, the $1\sigma$ and $2\sigma$ bounds on $\zeta_{H_2}(N_s)$ (i.e., the values for which the posterior probability of $\zeta_{H_2} < V$ is 0.023, 0.16, 0.84, and 0.972 from bottom to top). As expected, the uncertainties in $\zeta_{H_2}(N_s)$ become larger when $N_s$ lies outside the range that is primarily probed in this study.

# 5. Discussion

## 5.1. Attenuation of CR

Figure Set 8 shows clearly that CR are attenuated as they enter the interior of B68. The spatial variation of the line intensities rules out a model ($\alpha = 0$) where $\zeta_{H_2}$ is constant, and the decline of $I_{-7}/N_{22}$ with $N_{22}$ is the result of dust extinction alone. Because of the degeneracy between $\alpha$ and $N_0$, the probability distribution for $\alpha$ is quite broad, but values smaller than 0.31 are robustly excluded ($P \le 1 \times 10^{-2}$) in all four cases (X, Y, Z, and Q).

Although the best-fit parameters we obtain are quite similar in all four geometries, we can identify case X (propagation direction along the slit) as the most accurate representation of our observations of B68, given observational estimates of the magnetic field direction and the assumption that CR propagate along field lines. These estimates, obtained by R. Kandori et al. (2020) from near-IR polarimetry, indicated that the magnetic field direction, projected onto the plane of the sky, was at a position angle of $47^\circ \pm 5^\circ$ east of north and therefore close to our slit position angle of $57^\circ$. R. Kandori et al. (2020) were also able to derive an estimated inclination angle of $70^\circ \pm 10^\circ$ to the line of sight (i.e., concluded the field is mostly in the plane of the sky). This latter estimate is more model-dependent
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Separation of UVXH₂ and CRXH₂ emission. Top: integrated line intensities at the OFF position, divided by predicted UV fluorescence fractions, from the 100 K models of S88. Ortho-H₂ and para-H₂ transitions are shown in blue and red. Bottom panel: residual B68 line strengths attributable to CR excitation.</figcaption></figure>

than the position angle estimate, because it assumes a specific magnetic field configuration for the source.

The attenuation of the CRIR with the column density that we infer provides important insights into the possible mechanisms of CR transport in dense molecular clouds, and may also allow us to put constraints on the energy spectrum of

interstellar CRs. Here, we summarize and briefly discuss major signatures expected for different transport regimes.

If the transport of CRs penetrating the core were in the regime of free streaming (M. Padovani et al. 2009, 2018, hereafter P18), their attenuation would depend critically on the slope of the low-energy (nonrelativistic) spectrum of interstellar CR protons.
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Spatial variation of CRXH₂ line intensities across B68. Shown are three of the strongest lines, plotted as a function of position along the slit (see Figure 1). Intensities are corrected for UV contamination using the OFF position, and are given in units of 10⁻⁷ erg cm⁻² s⁻¹ sr⁻¹.</figcaption></figure>

Assuming a CR proton energy spectrum of the form $j(E) \propto E^a/(E_0 + E)^b$ (A. V. Ivlev et al. 2015), with $E_0 = 650$ MeV and $b = a + 2.7$ (P18), the attenuation is predicted to be practically negligible for sufficiently hard spectra with $a > \sim -0.2$. For smaller $a$, the CR attenuation index $\alpha$ decreases approximately linearly with the spectral index, $a$ (see Equation (33) in K. Silsbee & A. V. Ivlev 2019). However, such a spectrum would yield a substantially higher reference value of the CRIR than the derived $\zeta_{\text{ref}}$, as depicted in the bottom panel of Figure 10. At the same time, the model spectrum $\mathcal{L}$ (representing the Voyager data) would substantially underpredict the reference value and lead to practically no attenuation of the CRIR. The black curves plotted in Figure 9 assume the free-streaming regime, apply specifically for an assumed value of 650 MeV for

$E_0$, and include only the effects of ionization by CR protons.<fnref n="16" /> They assume that the total ionization rate is enhanced by a factor of 1.51 due to heavy elements and a factor of 1.67 due to secondary ionizations.

<footnote n="16">¹⁶ While it was shown by P18 that primary electrons may substantially enhance the ionization rate at low column densities, this was predicated upon the assumption that the energy spectrum of electrons measured by Voyager may be extrapolated in energy space with a constant power-law slope to arbitrarily low energy. There are two reasons to believe that this overestimates the role of primary electrons. First, models of CR acceleration produce a power law in momentum space, not energy space. Second, energy losses due to interaction with the interstellar medium are likely to further suppress the flux of CR electrons at low energies. The possible effect of CR electrons on the variation of the ionization rate with the shielding column will be investigated elsewhere.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Ratio of CRXH₂ line intensity to H₂ column density, *I*<sub>-7</sub>/*N*<sub>22</sub>, as a function of *N*<sub>22</sub>. The decrease with increasing column density reflects the combined effects of CR attenuation and dust extinction.</figcaption></figure>

Diffusive propagation of CRs, associated with preexisting MHD turbulence that is able to resonantly scatter CR particles penetrating molecular clouds, generally predicts a steeper attenuation of the CRIR for a given slope of the interstellar spectrum (K. Silsbee & A. V. Ivlev 2019). While this feature makes the diffusive transport scenario more favorable—in the sense that the data plotted in Figure 9 can be fit better without any adjustment to *E*₀—the very possibility of efficient scattering of nonrelativistic particles by preexisting turbulence in dense cores remains an open question. The turbulent cascade must operate down to extremely small resonant scales, corresponding to the CR gyroradius (≈0.1 au for a 1 GeV proton), i.e., much beyond the ambipolar diffusion scale at which Alfvén waves are expected to be damped (see, e.g.,

R. Kulsrud & W. P. Pearce 1969; T. C. Mouschovias et al. 2011). Furthermore, the resonant MHD turbulence is carried by ions, whose density in dense cores can be dramatically reduced compared to diffuse envelopes (where the ionization fraction is set by photoionization); as a result, the diffusive regime of CR propagation may cease to operate at column densities over a few times 10<sup>21</sup> cm<sup>-2</sup> (see Figure 1 in K. Silsbee & A. V. Ivlev 2019).

While streaming into the clouds, interstellar CRs themselves can resonantly excite MHD waves in diffuse gas surrounding dense molecular clouds (J. Skilling & A. W. Strong 1976; A. V. Ivlev et al. 2018). This phenomenon can lead to self-modulation of penetrating CRs: particles are efficiently scattered at the self-generated waves and, as a result, the CR
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Structure of the Bonnor–Ebert sphere model used to describe B68. Blue curve: radial density profile, normalized to the central density. Red curve: normalized column density as a function of projected radius. The adopted model is consistent with the extinction profile of J. F. Alves et al. (2001).</figcaption></figure>

<table><caption>**Table 3**
Best-fit (Median) Parameters from Likelihood Analysis</caption><thead><tr><th>CR Propagation Direction</th><th>α</th><th>log<sub>10</sub>(N<sub>0</sub> cm<sup>-2</sup>)</th><th>ζ<sub>H2</sub> at N<sub>s</sub> = N<sub>ref</sub> (s<sup>-1</sup>)</th><th>Min. Reduced χ<sup>2</sup></th></tr></thead><tbody><tr><td>x (perp. to line of sight, || slit)</td><td>0.47<sup>+0.16</sup><sub>-0.08</sub></td><td>20.07<sup>+0.66</sup><sub>-0.70</sub></td><td>(1.41<sup>+0.09</sup><sub>-0.07</sub>) × 10<sup>-16</sup></td><td>1.53</td></tr><tr><td>y (perp. to line of sight, ⊥ slit)</td><td>0.49<sup>+0.18</sup><sub>-0.08</sub></td><td>20.20<sup>+0.90</sup><sub>-0.83</sub></td><td>(1.39<sup>+0.08</sup><sub>-0.07</sub>) × 10<sup>-16</sup></td><td>1.61</td></tr><tr><td>z (along line of sight)</td><td>0.47<sup>+0.15</sup><sub>-0.08</sub></td><td>20.12<sup>+0.90</sup><sub>-0.77</sub></td><td>(1.42<sup>+0.12</sup><sub>-0.08</sub>) × 10<sup>-16</sup></td><td>1.64</td></tr><tr><td>xyz (quasi-isotropic; all directions)</td><td>0.48<sup>+0.20</sup><sub>-0.08</sub></td><td>20.16<sup>+0.86</sup><sub>-0.79</sub></td><td>(1.41<sup>+0.12</sup><sub>-0.07</sub>) × 10<sup>-16</sup></td><td>1.59</td></tr></tbody></table>

**Note.** Parameters are derived from fits to the spatial trend of $I_{-7}/N_{22}$ in B68 using the attenuation law $\zeta_{H2} = \zeta_0 [1 + N_s(H_2)/N_0]^{-\alpha}$. Quoted $\zeta_{H2}$ values are evaluated at the reference shielding column $N_s(H_2) = N_{ref} = 3 \times 10^{21} \text{ cm}^{-2}$. The values given are the median values with 68.3% confidence limits.

spectrum in the cloud is reduced compared to the interstellar spectrum at lower energies. Recently, this propagation mechanism was tested with GeV gamma-ray observations of nearby giant molecular clouds, showing excellent quantitative agreement with the theory (D. O. Chernyshov et al. 2024). For nonrelativistic CRs, the self-modulation is predicted to dramatically reduce the flux of penetrating particles starting from column densities of a few times $10^{21} \text{ cm}^{-2}$, which may lead to a very steep attenuation of CRIR.

In a separate paper (E. Makarenko et al. 2025, in preparation), we will present a detailed comparative analysis of different propagation models for B68.

## 5.2. CRIR

As indicated in Table 3, the CRIR at the reference shielding column of $3 \times 10^{21} \text{ cm}^{-2}$ is $\zeta_{H2} \sim 1.4 \times 10^{-16} \text{ s}^{-1}$, including both primary and secondary ionizations. This value, which is only very weakly dependent on the assumed CR propagation direction, is somewhat higher than the values typically obtained from recent analyses of the $H_3^+$ abundances in diffuse and translucent molecular clouds. Taken together, the studies of M. Obolentseva et al. (2024) and N. Indriolo et al. (2025) have identified 16 sight lines that exhibit $H_3^+$ absorption and

yield $\zeta_{H2}$ values ranging from $2.7 \times 11.5 \times 10^{-17} \text{ s}^{-1}$. The average value in these 16 clouds was $5.3 \times 10^{-17} \text{ s}^{-1}$, a factor of 3 smaller than the value determined for B68 (although the largest value determined from $H_3^+$ absorption observations was smaller by only 18%). Properly accounting for CR attenuation in the $H_3^+$ method would probably enhance the discrepancy somewhat since the typical shielding column densities probed in M. Obolentseva et al. (2024) and N. Indriolo et al. (2025) are a few times lower than our reference shielding column. The two methods being compared here, the CRXH$_2$ method and the $H_3^+$ method, are entirely independent and undoubtedly have different systematics. In particular, the CRIR determination from $H_3^+$ is proportional to the assumed volume density in the absorbing cloud, the latter being derived from extinction maps. If density estimates are taken from measurements of $C_2$ absorption line ratios instead of from extinction maps, the CRIR inferred from $H_3^+$ increases systematically by a typical factor of 2–3 (N. Indriolo et al. 2025), bringing the range inferred from $H_3^+$ into comfortable agreement with our estimate for B68. In any case, it remains unclear to what extent any differences between the CRIR derived in B68 and those typically obtained from $H_3^+$ observations reflect a true variation in the CRIR and to what extent it is a systematic discrepancy between the different methods. Further
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Posterior distributions for CR attenuation parameters, derived from fits to the $v = 1-0$ O(2) line. (a) Observed ratio $I_{-7}/N_{22}$ compared to the best-fit model. (b) Joint posterior distribution of $\alpha$ and $\log_{10} N_0$; best-fit values shown by black dot. (c)–(e) Marginalized distributions for $\log_{10} N_0$, $\alpha$, and $\log_{10} \zeta_{\text{ref}}$. Results apply to case X (propagation direction parallel to the slit). Similar results for alternative CR propagation geometries are shown in Figures 8.2–8.4 (online version only).
(The complete figure set (4 images) is available in the online article.)</figcaption></figure>

observations using the CRXH₂ method are urgently needed to extend the analysis beyond the single source considered in the present study.

### 5.3. Para-H₂ Line Ratios

Our observations of multiple para-H₂ lines provide a valuable probe of the excitation mechanism. In particular, the $v = 2-1$ O(2) / $v = 1-0$ O(2) intensity ratio, which we will refer to as $R_{21}$, is a powerful discriminant between CR and UV excitation. For the difference spectrum, B68–OFF, the intensities given in Table 2 above yield $R_{21} = 0.197 \pm 0.029$. This value lies between the predictions for UVXH₂ and CRXH₂, which differ significantly because of the different physical mechanisms involved: UV radiation pumps electronic

states, leading to a radiative cascade, whereas CRXH₂ emissions are excited directly by CR secondary electrons that excite to $v = 1$ much more efficiently than to $v = 2$. UVXH₂ emissions yield $R_{21} = 0.58$ (S88), a value that we confirmed with the Meudon photon-dominated region (PDR) code (F. Le Petit et al. 2006) and found to vary by less than 3% over many orders of magnitude in the assumed gas pressure. Because $R_{21}$ involves lines of para-H₂ alone, the UVXH₂ prediction is independent of the OPR in the UV-excited gas. For CRXH₂, the expected ratio is $R_{21} = 0.038$, based on cross sections given by M. C. Zammit et al. (2017) and L. H. Scarlett et al. (2023). Thus, the value of $R_{21}$ we observed strongly suggests that both CRXH₂ and UVXH₂ contribute to the signal measured in the difference spectrum. The pure UVXH₂
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Estimated $\zeta_{H_2}(N_s)$ as a function of $N_s$ (see the text for details).</figcaption></figure>

prediction is ruled out at the $13\sigma$ level, providing compelling evidence for the presence of a $\text{CRXH}_2$ component.

5.4. *Approximations Made in the Analysis*

We complete the discussion section with a review of several approximations made in the analysis.

1. We have assumed the B68 cloud to be spherically symmetric. While the projected extinction map shows only modest asymmetries, the structure along the line of sight is unknown.

2. We have adopted specific assumptions about the properties of the dust. In particular, we used the KP5 extinction curve tabulated in K. M. Pontoppidan et al. (2024), which includes the effects of ice absorption features, and assumed an $N_{22}/A_V$ ratio of 0.094.

3. In computing the attenuation of the $\text{CRXH}_2$ emissions, we applied a simple extinction correction: this treatment neglects the contribution of scattered radiation to the observed flux.

4. We assumed the gas within B68 to be fully molecular.

5. We assumed that $\text{CRXH}_2$ emission makes a negligible contribution to the $\text{ortho-H}_2$ line intensities we observed, and that the relative line strengths of the $\text{H}_2$ emissions attributable to the $\text{UVXH}_2$ component in B68 are the same as those in the background region.

In a future study, we will use detailed PDR models and dust radiative transfer models with scattering to remove the need for approximations 3, 4, and 5.

6. Summary

Our JWST/NIRSpec observations of the starless dark cloud B68 have revealed new insights into the excitation of molecular hydrogen by CRs. The main results of this study are as follows:

1. Multiple rovibrational $\text{H}_2$ lines attributable to CR excitation were detected across B68. Their spectra, dominated by para-$\text{H}_2$ transitions, clearly distinguish them from UV-pumped fluorescence.

2. By comparing the B68 spectra to a nearby OFF position, where the relative line strengths are in excellent agreement with models for UV excitation, we obtained line intensities that are robustly attributable to CR excitation alone.

3. The $\text{CRXH}_2$ line intensities vary systematically with position across the cloud, declining relative to the column density with increasing depth. This trend cannot be explained by dust extinction alone and demonstrates the attenuation of CRs within B68.

4. A Bonnor–Ebert sphere model was used to reproduce the observed spatial dependence of the $\text{CRXH}_2$ lines. Fits to the data require attenuation of the CRIR with depth, with
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Comparison of observed CR attenuation with theoretical predictions for the free-streaming regime. Top panel: dependence of the attenuation index, $\alpha$, on the spectral slope, $a$, of the low-energy CR proton spectrum. The horizontal shaded band indicates the 68% confidence interval for $\alpha$ inferred from our fits to B68. Bottom panel: predicted CRIR at the reference column density, $\zeta_{\text{ref}}$, for representative model spectra (P18), compared to the value inferred for B68 (vertical shaded band).</figcaption></figure>

values of the attenuation parameter $\alpha < 0.31$ strongly excluded.

5. At a reference column density of $N(\text{H}_2) = 3 \times 10^{21} \text{ cm}^{-2}$, the best-fit models yield a CRIR of $\zeta_{\text{H}_2} \sim 1.4 \times 10^{-16} \text{ s}^{-1}$, a factor of ≈3 higher than the average values derived from $H_3^\alpha$ absorption in diffuse and translucent clouds.

6. These results establish $\text{CRXH}_2$ emission as a powerful new diagnostic of CR penetration in dense molecular gas, with important consequences for our understanding of ionization, chemistry, and the role of low-energy CRs in star-forming regions.

## Acknowledgments

This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy,

Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with program #5064. We are grateful to J. Alves for providing an unpublished extinction map of B68 in digital format. Support for D.A.N. and K.S. was provided by NASA through a grant from the Space Telescope Science Institute. M.P. acknowledges the INAF grant 2023 MERCATOR ("Multiwavelength signatuRes of Cosmic rays in sTar-forming Regions") and the INAF grant 2024 ENERGIA ("ExplorNG low-Energy cosmic Rays through theoretical Investigations at INAF"). S.B. acknowledges financial support from the Israeli Science Foundation grant 2071540, the German-Israeli Science Foundation Grant 1568, and the Technion. B.A.L.G. is supported by the German Research Foundation (DFG) in the form of an Emmy Noether Research Group—DFG project #542802847 (GA 3170/3-1). C.H.R. acknowledges the support of the Deutsche Forschungsgemeinschaft (DFG; German Research Foundation) Research Unit "Transition discs"—325594231. C.H.R. is grateful for support from the Max Planck Society.
Software: JWST Calibration Pipeline (v1.16.1) with Calibration Reference Data System context 1303 (H. Bushouse et al. 2024), astropy (Astropy Collaboration et al. 2013, 2018, 2022), SciPy, NASA Astrophysics Data System (ADS), OpenAI GPT-5.

## Data Availability

The data underlying this article are publicly available from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute. These observations were obtained as part of JWST program 5064 (PI: Bialy). The specific observations analyzed can be accessed via doi:10.17909/1bw2-ga13.

## Appendix
### Expressions for $\zeta_{H_2}$

The expressions for $\zeta_{H_2}$ in the four cases discussed in Section 4 are as follows:

1. Case X: CR propagation in the plane of the sky along the slit direction, with an equal contribution in the $-x$ and $+x$ directions. Here, $\zeta_{H_2}(x, z) = \zeta_0 [F_s(N_{+x}) + F_s(N_{-x})]/2$, with

$$N_{+x} = \int_{-x_0}^{x} n_{H_2} ([x'^2 + z^2]^{1/2}) dx'$$

$$N_{-x} = \int_{x}^{x_0} n_{H_2} ([x'^2 + z^2]^{1/2}) dx'$$

and $x_0 = (r_{out}^2 - z^2)^{1/2}$. $F_s$ is the adopted shielding function ($1 + N_s/N_0)^{-\alpha}$ (Section 4), and $n_{H_2}(r)$ is the $H_2$ particle density as a function of the radial coordinate.

2. Case Y: CR propagation in the plane of the sky perpendicular to the slit direction, with an equal contribution in the $-y$ and $+y$ directions. Here $\zeta_{H_2}(x, z) = \zeta_0 [F_s(N_{+y}) + F_s(N_{-y})]/2$, with

$$N_{+y} = N_{-y} = \int_{0}^{y_0} n_{H_2} ([x^2 + y'^2 + z^2]^{1/2}) dy'$$

and $y_0 = (r_{out}^2 - x^2 - z^2)^{1/2}$.

3. Case Z: CR propagation along the line of sight, with an equal contribution in the $-z$ and $+z$ directions. Here $\zeta_{H_2}(x, z) = \zeta_0 [F_s(N_{+z}) + F_s(N_{-z})]/2$, with

$$N_{+z} = \int_{-z_0}^{z} n_{H_2} ([x^2 + z'^2]^{1/2}) dz'$$

$$N_{-z} = \int_{z}^{z_0} n_{H_2} ([x^2 + z'^2]^{1/2}) dz'$$

and $z_0 = (r_{out}^2 - x^2)^{1/2}$.

4. Case Q: Quasi-isotropic CR propagation, with an equal contribution in the $\pm x$, $\pm y$, and $\pm z$ directions. Here, $\zeta_{H_2}(x, z) = \zeta_0 [F_s(N_{+x}) + F_s(N_{-x}) + F_s(N_{+y}) + F_s(N_{-y}) + F_s(N_{+z}) + F_s(N_{-z})]/6$.

## ORCID iDs

David A. Neufeld https://orcid.org/0000-0001-8341-1646

Kedron Silsbee https://orcid.org/0000-0003-1572-0505

Alexei V. Ivlev https://orcid.org/0000-0002-1590-1018

Shmuel Bialy https://orcid.org/0000-0002-0404-003X

Brandt A. L. Gaches https://orcid.org/0000-0003-4224-6829

Marco Padovani https://orcid.org/0000-0003-2303-0096

Sirio Belli https://orcid.org/0000-0002-5615-6018

Thomas G. Bisbas https://orcid.org/0000-0003-2733-4580

Benjamin Godard https://orcid.org/0000-0003-0060-8887

James Muzerolle Page https://orcid.org/0000-0002-5943-1222

Christian Rab https://orcid.org/0000-0003-1817-6576

## References

Alves, J. F., Lada, C. J., & Lada, E. A. 2001, Natur, 409, 159

Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Bialy, S. 2020, CmPhy, 3, 32

Bialy, S., Chemke, A., Neufeld, D. A., et al. 2025, arXiv:2508.20168

Bohlin, R. C., Savage, B. D., & Drake, J. F. 1978, ApJ, 224, 132

Bushouse, H., Eisenhamer, J., Dencheva, N., et al. 2024, JWST Calibration Pipeline, v1.16.1, Zenodo, doi:10.5281/zenodo.14153298

Chernyshov, D. O., Ivlev, A. V., & Kiselev, A. M. 2024, PhRvD, 110, 043012

Draine, B. T. 2003, ApJ, 598, 1017

Evans, N. J., Dunham, M. M., Jørgensen, J. K., et al. 2009, ApJS, 181, 321

Hotzel, S., Harju, J., Juvela, M., Mattila, K., & Haikala, L. K. 2002, A&A, 391, 275

Indriolo, N., Obolentseva, M., Ivlev, A. V., et al. 2025, arXiv:2511.13915

Ivlev, A. V., Dogiel, V. A., Chernyshov, D. O., et al. 2018, ApJ, 855, 23

Ivlev, A. V., Padovani, M., Galli, D., & Caselli, P. 2015, ApJ, 812, 135

Jakobsen, P., Ferruit, P., Alves de Oliveira, C., et al. 2022, A&A, 661, A80

Kandori, R., Tamura, M., Saito, M., et al. 2020, PASJ, 72, 8

Kulsrud, R., & Pearce, W. P. 1969, ApJ, 156, 445

Lacy, J. H., Sneden, C., Kim, H., et al. 2017, ApJ, 838, 66

Le Petit, F., Nehmé, C., Le Bourlot, J., et al. 2006, ApJS, 164, 506

Mouschovias, T. C., Ciolek, G. E., & Morton, S. A. 2011, MNRAS, 415, 1751

Nielbock, M., Launhardt, R., Steinacker, J., et al. 2012, A&A, 547, A11

Obolentseva, M., Ivlev, A. V., Silsbee, K., et al. 2024, ApJ, 973, 142

Padovani, M., Bialy, S., Galli, D., et al. 2022, A&A, 658, A189

Padovani, M., Galli, D., & Glassgold, A. E. 2009, A&A, 501, 619

Padovani, M., Ivlev, A. V., Galli, D., & Caselli, P. 2018, A&A, 614, A111

Pontoppidan, K. M., Evans, N., Bergner, J., et al. 2024, RNAAS, 8, 68

Roy, A., André, P., Palmeirim, P., et al. 2014, A&A, 562, A138

Scarlett, L. H., Rehill, U. S., Zammit, M. C., et al. 2023, PhRvA, 107, 062804

Skilling, J., & Strong, A. W. 1976, A&A, 53, 253

Silsbee, K., & Ivlev, A. V. 2019, ApJ, 879, 14

Storey, P. J., & Hummer, D. G. 1995, MNRAS, 272, 41

Sternberg, A. 1988, ApJ, 332, 400

Zammit, M. C., Savage, J. S., Fursa, D. V., et al. 2017, PhRvA, 95, 022708

Zhu, H., Tian, W., Li, A., et al. 2017, MNRAS, 471, 3494