JWST Reveals Widespread CO Ice and Gas Absorption in the Galactic Center Cloud G0.253+0.016

Images were processed with a slightly modified version of the JWST pipeline based on version 1.11.1 (Bushouse et al. 2023). The tweakreg command was run on long-wavelength NIRCam data using the VISTA variables in the Vía Láctea (VVV) DR2 catalog (Saito et al. 2012) as an astrometric reference instead of the Gaia catalog (there were too few stars detected in common by both Gaia and JWST, and most were saturated). We then created a reference catalog based on the F405N catalog, which had fewer saturated bright stars than F410M and therefore more good associations with the ground-based near-infrared data. In the make_reftable.py script, we cut the F405N catalog based on crowdsource quality flags (qf > 0.95, spread < 0.25, frac > 0.9). This reference catalog was then used as the input to tweakreg for the other filters. We found that the tweakreg pipeline did not adequately correct the image registration to the absolute coordinates we provided, so we manually cross-matched the catalogs and computed shifts using the realign_to_catalog and merge_a_to_b functions in align_to_catalogs.py.

The narrowband filters, particularly F187N, F212N, and F466N, exhibited significant "streaking" noise that is strongly evident in the low-signal regions of the image, i.e., the majority of the molecular cloud. This streaking was caused by 1/f noise in the detectors (STScI helpdesk ticket INC0181624). As a first preprocessing step, we performed "destreaking" on each detector following a method suggested by Massimo Robberto (private communication), in which:

• 1.  
  Each detector was split into four horizontal quadrants with width 512 pixels and height 2048 pixels.
• 2.  
  The median across the horizontal axis was calculated, resulting in a 2048 pixel array.
• 3.  
  In a slight departure from Massimo's method, we then smoothed the median array using a 1D median filter with length that varied depending on the filter. 'F410M': 15, 'F405N': 256, 'F466N': 55, 'F182M': 55, 'F187N': 256, 'F212N': 512. We found that the original method, which was to obtain a single constant at this step, turned the destreaking process into a high-pass filter and therefore removed significant extended emission.
• 4.  
  We subtracted the median array from each quadrant, then added back the smoothed median.

We evaluated the effectiveness of this process by eye. While the original destreaker completely removed the 1/f horizontal features, it also removed all of the extended background. The modified version removed most of the horizontal features while preserving the large-scale extended structure. It is likely that some intermediate-scale features (i.e., physical features comparable to 512 pixels across) are not recovered by this procedure, which will need to be accounted for in analysis of the extended emission.

For photometry of unsaturated stars, we used the crowdsource python package (Schlafly 2021). We used a point-spread function (PSF) model from webbpsf (Perrin et al. 2015). Because we were using mosaicked images, the webbpsf PSF is not a perfect representation of the data; each individual frame had to be shifted and drizzled to form our final images. We therefore used the webbpsf model smoothed by a small amount, 0.05 pixels, as our PSF. We experimented with other PSFs and found that this approach was good enough, but it is not quantitatively optimal. Furthermore, for the short-wavelength bands, we adopted the PSF for a single detector (NRCA1 or NRCB1 as appropriate) for the full frame, since webbpsf does not provide a tool to produce a PSF grid across the whole module.

We identified saturated stars so that we could ignore them, and stars too close to them (which are likely to be affected by the extended PSFs of saturated stars), for further analysis. To identify these stars, we measured the centroids of all regions in which either the data (FITS extension SCI) or the variance (FITS extension VAR_POISSON) was zero. We then fitted the PSF of these stars excluding the saturated pixels and a surrounding set of pixels identified through binary dilation. The detailed values of the dilation size are given in saturated_star_finding.py. While we measured both photometry and astrometry of these saturated stars, we use only the astrometry in subsequent sections, and only as a means to automatically exclude saturated stars and their nearest neighbors.

We assembled a catalog consisting of all sources found in any of our six filters. To assemble the coordinate list, we started with all coordinates in the F405N catalog, then for each other filter, we added all sources that did not have a match in the existing catalog within d < 015. The cross-match shows a large peak for matches within d < 01 with large tails at greater separation; we have not investigated the origin of these large-offset sources, but suspect low signal-to-noise sources, PSF artifacts, and features in the extended background may contribute. We then excluded all sources with a cross-match distance to the reference filter's coordinates d > 01 (the reference filter is F405N by default, but for sources with nondetections in F405N, a different reference filter was adopted).

For subsequent analysis, we then rejected all sources with magnitude errors σ_m > 0.1, "quality factor" qf ≤ 0.6, spread > 0.25, or fracflux < 0.8, all of which are values calculated by crowdsource. These choices select for round, pointlike, unblended stars. We also limited our analysis to sources with detections in all six bands.

The resulting cross-matched catalogs had rms positional offsets <002 (Table 1). We do not concern ourselves further with astrometry in this manuscript, but caution that, with these sizeable cross-match errors, our catalogs likely are not yet of sufficient quality to support measurements of proper motion.

We found 377,236 stars with a good measurement in at least one filter, and 56,146 with good measurements in all six filters. The number of good measurements found in each filter is given in Table 1 (these include sources with offsets from the reference filter d > 01).

For comparison of star locations to extinction features, we preferred to work with an image with stars removed. Note that, because of significant uncertainty in this process, we have used the star-subtracted images only for qualitative, not quantitative, analysis. A starless image is the natural residual of an image that has been processed through a PSF photometry routine that appropriately accounts for the non-point-source background. However, when we produced such images, they had substantial residual features, which were caused by a combination of an imperfect PSF model and oversubtraction of sources on extended backgrounds. To create a cleaner starless image, we took the difference between the narrowband and medium-band images after appropriately scaling the narrowband image. The F405N image was convolved with a 0.3 pixel Gaussian to better match the PSF of the F410M filter.

We produced a line-free F410M image, labeled 410m405, by using the following equation:

$$\begin{eqnarray}&&{S}_{410{\rm{m}}405}=\displaystyle \frac{{S}_{{\rm{F}}410{\rm{M}}}-{S}_{{\rm{F}}405{\rm{N}}}\times w}{1-w}\end{eqnarray} \tag{ 1 }$$

where w is the fractional bandwidth of F410M covered by F405N and S_filtername is the surface brightness in a given filter. This process effectively creates a continuum-only "notch" filter image. We produced a star-free F405N image, labeled 405m410, by subtracting the (theoretically continuum-only) 410m405 image from F405N. Note that the images are in units of surface brightness, MJy sr⁻¹, such that line emission in broadband filters is diluted (will have a lower surface brightness), while spectrally flat continuum sources (to a coarse approximation, stars) will have the same brightness in broad and narrow filters. We then produced a somewhat star-free F466N image, which we label 466m410, by subtracting the 410m405 image scaled by R = (4.66/4.10)⁻² = 0.77 (assuming a blackbody on the Rayleigh–Jeans tail, S_λ ∝ λ⁻²).

$$\begin{eqnarray}&&{S}_{466{\rm{m}}410}={S}_{{\rm{F}}466{\rm{N}}}-{S}_{410{\rm{m}}405}\times R.\end{eqnarray} \tag{ 2 }$$

This S_466m410 image has much greater residuals, since the wavelengths do not overlap and differences in dust extinction (and ice absorption; see below) render the subtraction somewhat poor. Nevertheless, the stars are largely removed, and in particular, their extended PSFs are mitigated. The subtraction process is recorded in the notebooks

• BrA_separation_nrca.ipynb,
• BrA_separation_nrcb.ipynb,
• F466N_separation_nrca.ipynb,
• and F466N_separation_nrca.ipynb.

This process still left significant residuals throughout both the 405m410 and 466m410 images. To further remove stars—at this stage, purely for aesthetic purposes—we identified the locations of significant residuals and masked them out, then interpolated across them. We performed this process iteratively, using larger masks for stars with more extended PSF features and smaller masks for more compact, fainter stars. The details of the process were largely decided "by hand," i.e., testing a small variation in a parameter (e.g., the mask size) and revising if it did not look good. We also created custom masks to remove residuals from extended PSFs. The masking process is recorded in the StarDestroyer_nrca.ipynb and StarDestroyer_nrcb.ipynb notebooks. After each module was fully star-subtracted, the images were merged in the Stich_A_to_B.ipynb notebook.

Figure 1 shows the merged full-frame image. The color choices were made only for aesthetic reasons. Figure 2(a) shows the same image, and Figure 2(b) shows the version with stars before removal. Figure 2 also shows a subset of the cataloged stars as green and blue X's, as will be described in Section 4.1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The first intriguing result from these data is that, in colors including F466N, which is our longest-wavelength filter, the stars that are most extincted appear too blue. In general, dust extinction causes reddening, i.e., the shorter-wavelength photons are more attenuated than the long-wavelength photons; in this case, we see the inverse effect happening. This feature is evident in color–magnitude and color–color diagrams (CMDs and CCDs; Figure 4). Figure 4(c) shows a CCD with an extinction vector from Chiar & Tielens (2006, hereafter CT06) overlaid, ¹³ demonstrating that colors including the F466N filter go in a direction not accounted for by normal dust-extinction-driven reddening. We see excess absorption in the F466N filter of roughly 0.035 magnitudes per A_V with substantial scatter.

The bluest stars in [F410M] – [F466N] ¹⁴ are seen at the edge of the cloud. Figure 2 shows the location of these bluest ([F410M] – [F466N] < −0.45 mag) stars. The most extincted stars, which are the reddest in most colors but are bluest in colors involving F466N, are primarily seen along the outskirts of the cloud. Figure 4(a) shows stars color-coded by [F187N] – [F405N] color. We use [F187N] – [F405N] color as a consistency check: this color is very well correlated with both [F182M] – [F212N] and [F212N] – [F410M] color in panel (b), indicating that the trend seen in panel (c) is not caused by F410M. The interior of the cloud appears relatively blue in Figure 4(a) because only low-extinction stars are detected in the shorter-wavelength filters, and we plot only stars with detections in all six filters in this figure. Figure 4(b) highlights that colors not involving the F466N filter are consistent with extinction.

The diffuse emission circumscribing The Brick, seen in F405N and F466N in Figure 1, is composed primarily of Brα (F405N) and Pfβ plus Hu (F466N) emission. Figure 3 shows where these lines reside with respect to the transmission profiles of the filters. There are no other expected sources of emission in this band, as there are no known features of polycyclic aromatic hydrocarbons in the 4–5 μm range and free–free emission is expected to be weaker by ∼100× in the narrow bands.

The ratio of hydrogen recombination lines is governed by simple rules under the assumption of Case B recombination, which is expected at moderate densities. The expected ratio of Pfβ/Brα under Case B recombination at electron temperature T_e ∼ 5000–10⁴ K is R_Pfβ/Brα = 0.202 (Storey & Hummer 1995). The ratio of Pfβ+Hu, the sum of the two lines in F466N, to Brα is R_{(Pfβ+Hu)/Brα} = 0.234. The average foreground extinction toward the CMZ (Launhardt et al. 2002; Nogueras-Lara et al. 2021) is A_V = 30. Using a CT06 extinction curve, in the absence of narrow spectral features, the ratio above rises to R_{(Pfβ+Hu)/Brα} (A_V = 30) = 0.269. At greater extinction, this ratio (which is equivalent to [F405N] – [F466N] color) is expected to rise (become redder). However, contrary to this expectation, we see the [F405N] – [F466N] color becoming more negative (bluer) along the edge of The Brick (Figures 1 and 2). We are therefore seeing that, in regions of greater extinction, the ratio is the inverse of what is expected from dust extinction alone.

4.2.1. CO Absorption of the F466N Recombination Lines

As shown in Figure 3, the ¹²CO lines in the F466N band can produce significant absorption against stellar continuum light. In this section, we evaluate whether CO gas can produce the observed stellar colors. We already saw in Section 4.2 that CO gas is unlikely to produce selective extinction of Pfβ+Hu. We find here that CO gas contributes to, but does not dominate, the total absorption in F466N.

We model this absorption as a function of temperature, column density, and line width assuming LTE conditions. Details of this modeling are given in an associated Jupyter notebook, COFundamentalModeling.ipynb, that can be found in the associated GitHub repository. The model implementation is in the pyspeckit-models package, which implements models compatible with pyspeckit (Ginsburg et al. 2022). We used transition and level tables from the exomol database (Tennyson et al. 2016) derived from Li et al. (2015) using Yurchenko et al. (2018) as an implementation reference.

Figure 5 shows example optical depth spectra overlaid on the transmission profile of the F466N filter. The v = 1←0 and J = 4←3, 3←2, 2←1, 1←0 R-branch transitions and the v = 1←0 J = 0←1 P-branch transition all lie within the range where F466N has >50% of peak transmission (we use ← in the transition names to indicate that these are absorption lines). The maximum absorption in this filter occurs for temperatures between 10 and 20 K, while we expect gas temperatures near 50–100 K (Immer et al. 2012; Ginsburg et al. 2016; Krieger et al. 2017). At greater temperatures, a large fraction of CO molecules are in states that only produce transitions outside of the F466N band, reducing the absorption for a fixed assumed column density. It is possible that some of the CO gas is at moderate densities (n(H₂) ≲ 10⁴ cm⁻³) and therefore is subthermally excited, which would concentrate the CO molecules into the lower-J levels, thereby reducing the effect of high gas temperature. Nevertheless, the LTE models shown in Figure 5 capture the range of expected behavior.

Zoom In Zoom Out Reset image size

We model the CO gas absorption for the expected range of line width and column density in the Galactic Center. For narrow line widths, such as those caused by thermal broadening at T < 100 K, the absorption is negligible. In the Galactic Center, there is significant Doppler broadening that is generally attributed to turbulence. The total line width in the cloud may range from σ ∼ 5 to 20 km s⁻¹ (Henshaw et al. 2019). CO column densities will span the full range from effectively zero (since CO is destroyed by UV at A_V ≲ 2) to a few ×10¹⁹ cm⁻² (Rathborne et al. 2015) assuming X_CO = 10⁻⁴. In our observations, we detect stars at wavelengths short of 2 μm only at intermediate column densities, most likely below A_V < 80 mag (N(H₂) < 10²³ cm⁻²), since dust extinction hides the stellar continuum at greater column density at the current level of sensitivity.

Figure 6 summarizes the modeling results. Given the plausible range of column density and line width, the total CO gas absorption in F466N may range from ∼1% to at most ≲20%.

Zoom In Zoom Out Reset image size

At the column densities where the absorption is easily detectable (fractional absorption ≳0.1 results in a change in magnitude Δm ≳ 0.1; green line in Figure 6), change in line width dominates over the change in column density or temperature. The foreground Galactic disk clouds, which have narrow lines, produce relatively little absorption; we therefore argue that intervening material between us and the Galactic Center is not primarily responsible for the blue [F410M] – [F466N] color of the stars. These models also imply that, even at very extreme column densities (N_CO > 10¹⁹ cm⁻²), "normal" Galactic disk clouds with σ < 10 km s⁻¹ will produce minimal CO absorption in the F466N band, while typical clouds in the Galaxy Center and Galactic bar with σ ∼ 10–50 km s⁻¹ will produce readily detectable absorption. However, even for very broad lines (σ = 50 km s⁻¹) at high column (N_CO ∼ 10^19.5 cm⁻²), CO gas produces ≲1 magnitude (<50%) of absorption in the F466N band.

While we show above that CO gas can produce a substantial amount of absorption, the observed absorption depths reach levels difficult to explain with gas alone. Figure 6 shows that CO column densities N(CO) > 10^19.5 cm⁻², implying N(H₂) > 3 × 10²³ cm⁻², are required to explain the ≳2 magnitudes of F466N absorption shown in Figure 4. Such high column densities are rare in The Brick (Rathborne et al. 2014a), occurring primarily in the inner regions (see their Figure 2) and in dense cores (Walker et al. 2021), while the high-extinction and high-CO-absorption stars we detect are primarily in the outskirts (Figure 2). Stars behind these high column densities would be too extincted to detect in the shorter-wavelength band; the highest extinction we report in Section 4.4.2 is A_V ∼ 80 mag, or N(H₂) ∼ 10²³ cm⁻². Additionally, in Section 4.2, we showed that CO gas is unlikely to absorb recombination lines. We therefore examine the possibility that ice absorption is responsible for the observed F466N deficits.

There is evidence that The Brick contains some ice, but also that CO is not entirely frozen out. Pure CO ice forms at low temperatures, T < 20 K (Hudgins et al. 1993). The average dust temperature in The Brick is close to 20 K (Tang et al. 2021), so it is probable that some of the volume of The Brick is cold enough to freeze CO. The Brick exhibits signs of substantial freezeout in its center based on gas observations (Rathborne et al. 2014b), but still has substantial gas-phase CO detected (Ginsburg et al. 2016; Rigby et al. 2016; Eden et al. 2020). It is likely that much of the observed gas-phase CO is on the cloud surface, while further into the interior, CO is more completely frozen out.

4.4.1. CO Ice Modeling

To model CO ice absorption, we convolved a given filter transmission curve with a synthetic stellar model spectrum. We began with a 4000 K PHOENIX stellar atmosphere (Husser et al. 2013) as the base model, then examined the fractional flux lost in the F466N band as a function of CO column density. We retrieved optical constants for pure CO ice and CO mixed with OCS and CH₄ in a 20:1 ratio from Hudgins et al. (1993) via the JPL Optical Constants Database. ¹⁵ We do not know which ice mixture is most appropriate for our data, so we chose to show all available laboratory mixtures in which CO was the primary constituent. We also considered the possibility that CO is embedded in other ices (e.g., H₂O and CH₃OH, Pontoppidan et al. 2003; Boogert et al. 2008), but found little practical difference from pure CO ice when using the optical constants of Hudgins et al. (1993) and Rocha et al. (2016).

Figure 7 shows the effects of CO ice absorption ¹⁶ : at N(CO) ≈ 10¹⁹ cm⁻², we expect ≈0.5–1 mag of absorption from the ice band. The greatest observed column density within The Brick, based on dust emission observations from the Atacama Large Millimeter/submillimeter Array (ALMA) with ∼3'' resolution, is N(H₂) ∼ 5 × 10²³ cm⁻² (Rathborne et al. 2014a), which implies an upper limit on the ice column density N(CO) < 5 × 10¹⁹ cm⁻² if we assume the CO/H₂ ratio is 10⁻⁴ and the gas-to-dust mass ratio is 100. If all of the CO is frozen out into pure CO ice, in the highest column-density line of sight, the absorption in F466N could just about reach 1.2 mag (∼65%).

Zoom In Zoom Out Reset image size

Figure 8 shows similar plots for the F405N and F410M filters for CO₂ ice to demonstrate that CO₂ ice can have some effect, but less so than CO, on our observed colors.

Zoom In Zoom Out Reset image size

4.4.2. CO Ice as a Function of Extinction

We compared the inferred extinction to the CO ice column to obtain a coarse estimate of how CO ice varies with A_V . We use [F182M] – [F212N] color to estimate A_V by using the CT06 extinction curve, which we justify by comparing to ground-based GALACTICNUCLEUS colors in Appendix A.

To measure the CO absorption, we used the [F410M] – [F466N] color. Figure 4(c) shows that the [F410M] – [F466N] color becomes bluer at greater extinction. We calculated the absorption from pure CO ice in the F466N band to obtain a mapping from N(CO) to [F410M] – [F466N] color. We dereddened our measured [F410M] – [F466N] color using the A_V computed above with the CT06 extinction curve. Note that the assumed extinction curve introduces a strong systematic uncertainty: using a Fritz et al. (2011) extinction curve would reduce the estimated A_V by a factor of 1.8. Figure 9 shows the resulting N(CO) as a function of A_V . Note that this curve assumes that 100% of the deficit in the [F466N] band is produced by ice, which is an upper limit and likely an overestimate because CO gas (Section 4.3) also contributes, especially at low column density.

Zoom In Zoom Out Reset image size

The red curve in Figure 9 shows the maximum possible CO at each A_V adopting standard values of N(H₂) = 2.21 × 10²¹ A_V (Güver & Özel 2009) and CO abundance relative to hydrogen X_CO = 10⁻⁴. Some of the data points reside above this curve, suggesting that one or both of these assumptions may be incorrect, which we will evaluate further in Section 5.1.

While there is a good overall correlation between N(CO) and A_V , the dispersion at any given A_V is large, greater than one order of magnitude. Such a large scatter indicates a wide range of conditions, with many lines of sight containing little CO at the lowest observed A_V .

The CO ice models can qualitatively explain most of the observed F466N deficits. However, neither the gas nor ice models appear to quantitatively explain the most deeply absorbed sources. In both the gas and ice absorption cases, it is possible to achieve substantial absorption in the F466N band, enough to be easily detected, but less than the 1–2 mag seen in Figure 4 if typical CO abundance and gas-to-dust ratios are used.

One possible explanation lies in our assumptions about how to convert column densities between dust, molecular hydrogen, and CO: if the CO/H₂ ratio is greater than the assumed 10⁻⁴, or the gas-to-dust ratio is less than 100 (e.g., Giannetti et al. 2017), the total CO column could be greater. Increasing the CO abundance or decreasing the gas-to-dust ratio, both of which are plausible because gas in the CMZ has a higher metallicity than the solar neighborhood, would both have the effect of shifting the red dashed line upward in Figure 9. These changes would therefore increase the maximum allowed CO abundance and explain the large observed F466N absorption.

The presence of significant quantities of CO ice in The Brick highlights the fact that the dust is significantly colder than the gas. Gas temperatures in the CMZ generally (Ginsburg et al. 2016; Krieger et al. 2017), and The Brick specifically (Johnston et al. 2014), are observed to be high, T ≳ 50 K, and in many locations T > 100 K. Freezing of pure CO into ice is expected to occur at dust temperatures T ≲ 20 K, though CO can be integrated into H₂O and CH₃OH ices that freeze out at higher temperatures (T ≳ 80 K; Garrod & Herbst 2006; Boogert et al. 2015). The dust temperatures observed in The Brick have been in the range T ∼ 20–30 K, albeit at lower resolution (Marsh et al. 2016; Tang et al. 2021), which is somewhat too warm for pure CO freezeout but cold enough to freeze other molecules. However, dust temperature measurements are biased toward warmer dust, since it is brighter, so it is likely that colder dust is present deep inside The Brick.

Both an excess of CO in Galactic Center gas and freezeout during gravitational collapse may result in a change in the effective equation of state of the gas. In the dense molecular medium (n ≳ 10³ cm⁻³) that comprises The Brick, CO is the dominant gas-phase coolant (Ginsburg et al. 2016). If the CO abundance is greater than in the solar neighborhood (Section 5.1), we expect more efficient cooling in the lower-density outskirts of CMZ clouds. By contrast, as the clouds collapse to greater density, there may be a point at which the CO has frozen out to the point that it is no longer the dominant coolant, but where the densities are still too low for dust to be an efficient coolant. We suggest that systematic variations in the cooling function as a function of density or column density should be explored in future simulations of CMZ cloud thermodynamics like those in Clark et al. (2013).

The prevalence of CO ice in our own Galactic Center hints that ice is likely widespread in galactic centers generally. At least in the local universe, then, JWST observations using the long-wavelength narrowband filters should carefully treat CO absorption in addition to extinction effects. Our Figure 9 gives an empirical tool to link extinction and ice, though we caution that the large scatter demonstrated in that plot limits the usefulness of the linear relation given in its legend.

The easy detection of ices in the F466N narrowband filter also creates opportunities to better understand dust and ice in the interstellar medium and to understand cloud structures. While NIRSpec will be capable of studying hundreds of stars with ice absorption in detail, NIRCam observations can easily measure tens of thousands of sightlines simultaneously, enabling detailed correlation analyses like those shown in Figure 9. By adding comparable-resolution gas observations from ALMA, it should be possible to trace the freezeout of CO from gas to ice in detail. With a few other bands, such as F300M and F335M, it will be possible to track H₂O and CH₃OH ice and determine when and how much CO is incorporated into H₂O ice, which freezes at a substantially higher temperature.

The overlap of the CO ice feature with Pfβ+Hu also opens the possibility of making ∼01 resolution maps of CO ice absorption against the diffuse ionized emission. From the Paα/Brα ratio, we can determine the dust extinction on a per-pixel basis, which will enable specific measurement of CO ice absorption from the now known Brα/Pfβ ratio. Since the CO ice feature affects the Pfβ line, but CO gas does not, this approach will also allow us to distinguish whether ice or gas is the dominant absorber on most sightlines.