© The Authors 2025

1 Introduction

Ground-based and New Horizons observations have shown Pluto’s tenuous, time-variable, N₂-dominated, atmosphere to be alluringly Titan- and Mars-like. With a surface pressure (currently ~12 μbar, Sicardy et al. 2021) comparable to Titan’s at 400 km altitude, Pluto possesses a chemically-rich (N₂, CH₄ and other hydrocarbons, CO and HCN) atmosphere with an extensive haze representing the end-product of atmospheric photochemistry (see reviews, e.g., in Stern et al. 2018; Summers et al. 2021; Mandt et al. 2021). Similar to Mars, Pluto possesses a complex climatic system driven by the redistribution cycles of volatile ices (N₂, CH₄, CO, see review in Young et al. 2021) that ultimately shapes Pluto’s landscapes on seasonal and geological timescales.

Ground-based observations of Pluto’s atmospheric composition have been conducted in the near-infrared (NIR), yielding the detection and monitoring of CH₄ (Young et al. 1997; Lellouch et al. 2009, 2015), and in the millimeter and submillimeter (mm/submm) range with ALMA, allowing for the detection of CO, HCN, and HNC (Lellouch et al. 2017, 2022) and offering some spatial resolution. Unsurprisingly, the most detailed picture of Pluto’s atmospheric composition was achieved by New Horizons/Alice in July 2015 (Gladstone et al. 2016; Young et al. 2018; Kammer et al. 2020; Steffl et al. 2020). These observations provided vertically-resolved profiles of N₂ and the major hydrocarbons (CH₄, C₂H₂, C₂H₄, and C₂H₆) and haze in broad (albeit species-dependent) portions of the atmosphere from ultraviolet (UV) solar and stellar occultations, as well as the detection of several airglow emissions. These results set the foundation for constraining photochemical models (see review, e.g., in Mandt et al. 2021, and discussions hereafter). However, except for a claim for CH₃C₂H (Steffl et al. 2020), higher order hydrocarbons or more complex nitriles have not been found yet, although they are prevalent in Titan’s atmosphere (see review, e.g., in Bézard et al. 2014).

While Pluto’s haze was exquisitely characterized by New Horizons in terms of its vertical distribution and optical and scattering properties, much less is known of its composition. Meanwhile, laboratory experiments on simulated Pluto environments have yielded indications on its bulk (atom-wise) composition (Jovanović et al. 2020). Given the low temperatures in the upper atmosphere (~70 K), Pluto’s haze, initially thought to be of similar composition as Titan’s due to the similarity of chemical pathways, has lately been proposed to contain a major organic ice component from the direct condensation of the primary photochemical products, and to be dominated (by mass) by HCN ice above 350 km and C₄H₂ ice below that altitude (Lavvas et al. 2021). Knowledge of the haze composition is especially important because it affects its radiative properties. The latter have been suspected to largely control Pluto’s atmospheric heat budget and, in particular, to be responsible for the colder than anticipated upper atmosphere (Zhang et al. 2017; Wan et al. 2021). This hypothesis has recently been confirmed through rotationally resolved measurements of Pluto’s thermal flux with JWST/MIRI imaging over 15–25 μm, allowing us to distinguish the thermal emission of haze from that of the surface (Bertrand et al. 2025).

The frigid temperatures of Pluto’s atmosphere (~37.5 K at the surface, 110 K at the stratopause, and ≲70 K above 300 km) have made its thermal emission undetectable to instruments so far, except at mm/submm wavelengths, where the Planck function effect is less detrimental. This has left the entire mid- to far-IR range (5–500 μm) of Pluto’s atmospheric emission as a virgin territory. Yet, ever since the Voyager epoch (1979–1989), it has been known that the 5–30 μm range of the Giant Planets and Titan is rich in signatures of gases (and condensates) associated with the fundamental vibration modes of many species, particularly hydrocarbons and nitriles. Despite expected fluxes that are typically ~10⁴ times smaller on Pluto than on Titan, the unsurpassed sensitivity of JWST/MIRI now allows us to open this spectral range for Pluto. We report here on our measurements obtained within the JWST GO-1 program #1658, along with their implications for our knowledge of Pluto’s atmosphere.

2 Observations and data reduction

JWST MIRI observations of the Pluto-Charon system were obtained within the JWST GO-1 program #1658 (Pluto’s climate system with JWST, PI Lellouch). Data were taken on May 4, 2023 (start time: 07:39 UT; end time: 16:41). The observations were carried out using the Medium-Resolution Spectrometer (MRS; Wells et al. 2015; Argyriou et al. 2023), which consists of four integral field units (Channels 1–4) and three grating settings (A, B, and C) that together provide continuous coverage of the 5–28 μm wavelength range.

Observations were acquired in the form of a four-point dither. We used the SLOWR1 readout mode, with Ngroups = 20 per integration. For each dither and grating position, five integrations were acquired, leading to 9938 s (2.76 h) exposure time per grating, namely, 8.28 h total exposure time. During the entire ~9 h time interval covered by the observations, Pluto rotated by 21° (sub-observer longitudes 353.3 E–332.2 E at start/end). Pluto’s sub-observer and heliocentric distances were 34.5058 au and 34.7602 au respectively at mid-observation point, and the sub-observer latitude was 59.6° N.

The uncalibrated data were downloaded from the Mikulski Archive for Space Telescopes and passed through the first two stages of Version 1.14.0 of the official JWST calibration pipeline (Bushouse et al. 2024) to produce flat-fielded, wavelength- and flux-calibrated data cubes. The necessary reference files were automatically drawn from context jwst_1229.pmap of the JWST Calibration Reference Data System. The optional residual_fringe step was manually activated in order to reduce the level of fringing in the extracted spectra. All other parameter settings in the pipeline were left at their default values. Pairs of dithered exposures were subtracted from one another to remove most of the background flux; the thermal background becomes significant in Channel 4 (i.e., wavelengths longer than ~18 μm). The negative source that arose from dither subtraction was masked prior to spectral extraction.

We used a circular spectral extraction aperture with a diameter of 5 pixels for both Pluto and Charon. There was non-negligible residual flux at wide separations, indicating imperfect removal of the background. We therefore established a background annulus between 12 and 30 pixels from the centroid and removed the average flux level within this region from the extracted spectra. Next, the spectra were passed through the additional defringing routine ifu_rfcorr contained within the extract_1d step of the JWST pipeline, which reduced the relative amplitude of the fringes to within a few percent.

At the time of the observations, Charon was located roughly 0.7″ north-east of Pluto. The wings of the Pluto’s and Charon’s point-spread functions (PSFs) overlap significantly, requiring special care in the spectral extraction process. Due to the partial blend of the sources and the coarse spatial sampling (0.13″–0.20″), a blind centroiding was unreliable. We therefore used the calculated positions from the JPL Horizons ephemerides service to center the extraction apertures for Pluto and Charon in each data cube.

In a manner analogous to the methodology employed by Protopapa et al. (2024) in their analysis of Charon’s NIR JWST observations, we estimated the contamination of the binary companions’ PSFs in each other’s extraction apertures by defining contamination apertures that are situated directly opposite the companion within the field of view. The flux contained within the contamination aperture was then subtracted from the background-subtracted flux of the object of interest to produce the decontaminated spectrum.

Due to the changing size of the PSF as a function of wavelength, the fraction of the source flux that lies outside of the fixed circular aperture varies throughout the spectral range. To correct for this wavelength-dependent flux loss, we applied the technique described in Rivkin et al. (2023) and computed empirical correction vectors for each channel and grating setting. In short, we extracted the spectrum from the MIRI MRS observation of the A-type standard star del UMi (Program #1536; PI: K. Gordon) with the same five-pixel circular aperture. After dividing the measured spectrum by the corresponding CALSPEC model (Bohlin et al. 2014), we fit cubic polynomial trends to the resultant ratio arrays to obtain the correction vectors. These vectors were then applied to the extracted flux of Pluto and Charon to arrive at the final spectra. There were small, but measurable offsets between the adjacent segments of the spectra obtained with different grating settings, indicating a low level of rotational variation across the duration of the observations. These were accounted for by sequentially matching the median flux levels within the overlapping wavelength ranges.

The resulting Pluto spectrum is shown in Fig. 1. The MIRI MRS extracted fluxes over 15–21 μm are nicely consistent with, though slightly below (~10%), the Pluto fluxes determined from MIRI imaging photometry at similar Pluto longitude (342° E). Longwards of 21 μm and particularly in band 4C (24.40–27.90 μm), fluxes start to systematically deviate below the levels indicated by MIRI F2550W photometry. To avoid distorting the spectral shape over 18–21 μm, we chose to not attempt recalibrating the MRS spectrum against the MIRI photometry, and here we focus on the 4.9–20 μm part of the spectrum, where all atmospheric and solid state signatures occur. The Charon data will be presented elsewhere (Souza Feliciano et al., in prep.).

Figure 1 reveals a very rich spectrum, showing numerous surface and atmospheric features. Surface features include a broad absorption over 5–8 μm, dominated by CH₄ ice, and two narrower bands at 8.65 μm (due to CH₃D ice) and 10.5 μm (attributed to C₂H₄ ice). The reflectance features caused by surficial ices will be discussed elsewhere (Souza Feliciano et al., in prep.). Atmospheric features include prominent fluorescent emissions due to CH₄ and CH₃D, as well as strong “Titan-like” bands of C₂H₆, C₂H₂, CH₃C₂H, and C₄H₂. A broad, two-peaked emission over 14.3–16.7 μm is also clearly detected, which we attribute to haze emission.

The noise level indicated by the pipeline is formally 0.002–0.003 mJy over 5–18 μm, increasing to ~0.015 mJy at 23 μm. However, these noise levels certainly underestimate the true error bars. Based on the standard deviation of the spectrum after subtraction of a continuum model (red curve in Fig. 1) outside of atmospheric bands, we estimated a noise level of ~0.010 mJy over 5–7 μm, ~0.007 mJy over 9.2–11.5 μm, and ~0.010 mJy over 16.5–18 μm. This corresponds to continuum S/N of ~80 at 5 μm, ~25 at 11 μm and ~15 at 17 μm. The peak S/N on the gas emissions reaches ~50 for CH₃C₂H and C₄H₂, ~80 for C₂H₆ and up to ~ 180 for C₂H₂. In an astronomical perspective, the JWST/MIRI Pluto spectrum is comparable in quality to the ISO spectrum of Titan obtained in 1997 (Coustenis et al. 2003), even though the Pluto flux levels are of the order of 1 mJy, compared to tens of Jy to ~200 Jy for Titan.

Fig. 1 Overview of the MIRI MRS spectrum of Pluto, with atmospheric gas features identified (black). The light blue points show the Pluto fluxes from MIRI imaging photometry at a Pluto longitude of 342° E (Bertrand et al. 2025). Red curves indicate a model without atmospheric features. The model includes a solar component (blue) and a thermal component (green). Solid (resp. dashed) lines correspond to the diurnal-only (resp. CH4 seasonal) case for surface temperature. The I/F reflectivity for the solar reflected component is shown in the inset. It indicates absorptions attributed to CH4, CH3D, and C2H4 ices.

3 Analysis

The data were analyzed using our radiative transfer code, initially developed for analyzing Pluto ALMA observations (Lellouch et al. 2022) and fully accounting for the spherical geometry associated with the large extent of Pluto’s atmosphere compared to its radius. The model covers 0–1300 km altitude range with an altitude-varying layer thickness. Originally designed for the thermal range and for a gas-only atmosphere, the model was expanded to include thermal emission from haze and the contribution of a solar reflected component. For the solar spectrum, we used a continuum model taken from R.L. Kurucz¹. Weak lines are present in the Sun spectrum over 5–6 μm. Although they are not obviously seen in the Pluto spectrum, they were included in the model by using the line-to-continuum ratio spectrum provided by the ACE-FTS atlas of Hase et al. (2010). Pluto’s surface I/F reflectivity was adjusted to fit the overall spectral shape, including the ice features in the solar-reflected dominated spectral range (see Fig. 1, inset). The surface I/F reflectivity is constrained only up to ~13 μm (the onset of the C₂H₂ gas emission); longwards of 13 μm, it is assumed to retain its value at 13 μm (0.5).

Thermal emission from the surface was calculated based on the surface thermophysical and climate models from Bertrand et al. (2025). Those physically based models replicate Pluto’s thermal light curves observed by MIRI imaging in F1800W, F2100W, and F2550W filters (data also obtained in this program). Bertrand et al. considered both diurnal-only and CH₄ seasonal cycle temperature scenarios. Both scenarios reproduce the thermal light curve contrasts (i.e., the differential fluxes to the mean), but lead to absolute thermal fluxes over 18–25.5 μm that differ by a factor of ~1.5 between the two models.

We have considered both scenarios. The CH₄ seasonal cycle case is conceptually preferable as it includes more physics (i.e., seasonal temperature cycles and the cooling associated with methane sublimation). However, the diurnal-only case leads to the highest surface temperatures and as such corresponds to the most conservative (minimum) case of haze emission. Both models lead to a spectrally variable surface brightness temperature, rising from 46.4 K (resp. 46.0 K) at 28 μm to 50.2 (resp. 47.6 K) at 12.5 μm for the diurnal-only (resp. CH₄ seasonal) models. The cross-over point of the solar reflected and thermal components occurs near 17.3–18.2 μm (see Fig. 1).

For determining the abundance of detected gases and for deriving upper limits, we generally used photochemically-based profiles as templates. As reviewed in Mandt et al. (2021), numerous photochemical models of Pluto’s atmosphere have been developed or revived in the wake of the New Horizons exploration (Wong et al. 2017; Luspay-Kuti et al. 2017; Mandt et al. 2017; Krasnopolsky 2020a,b; Lavvas et al. 2021). In addition to differences in the adopted chemical pathways and reaction rates, the various models may differ in several aspects: (i) some of them include coupled ion-neutral photochemistry while others are restricted to neutral chemistry; (ii) models altogether invoke a wide range of eddy diffusion coefficients (10³ to 10⁶ cm² s⁻¹); and (iii) models may disagree on the fate of the secondary products (mostly hydrocarbons and nitriles), which can be atmospheric condensation, sticking onto aerosols, or condensation at the surface. Here, we adopted for reference the mixing ratio profiles from Lavvas et al. (2021). This model assumes a CH₄ mixing ratio at the surface of 0.4% and a constant with altitude eddy mixing profile of K_ZZ = 10³ cm² s⁻¹. It includes a self-consistent approach to Pluto’s gas-haze system coupling photochemistry and haze microphysics and is able to reproduce the line-of-sight (LOS) densities of CH₄, C₂H₂, and C₂H₆, as well as the haze opacities observed by New Horizons. Upon comparison to observations, we then allowed for the scaling of the modeled gas profiles by constant factors (except for C₂H₂, as explained in Sect. 3.3.2). For the detected gases, we also calculated mixing ratios assuming a vertically uniform distribution.

We also expanded the atmospheric radiative transfer code to include thermal emission from a haze with a specified vertical distribution and adjustable optical depth (see Sect. 3.2). For the thermal profile (which was not measured independently in our 2023 data), we used the New Horizons profile (Lavvas et al. 2021; Young et al. 2018, their Fig. 9). Although thermal profiles obtained from ALMA (Lellouch et al. 2017, 2022) are slightly different from the New Horizons profile, the use of the latter is justified by its reference character and for consistency with the Lavvas et al. (2021) gas profiles. Local thermodynamical equilibrium (LTE) was nominally assumed throughout the atmosphere when performing the radiative transfer calculations, but as detailed in Sect. 3.1, this approach was modified upon modeling of the fluorescent emissions. Spectroscopic parameters for the different species were taken from the GEISA (Delahaye et al. 2021) or HITRAN (Gordon et al. 2022) databases. Synthetic spectra were calculated at a 0.0005 cm⁻¹ step and convolved to the wavelength-dependent resolving power of MIRI MRS (Jones et al. 2023).

3.1 Methane and CH₃D fluorescence

The strong emission from CH₄ and CH₃D gases, in the 7.25–8.25 and 8.25–9.0 μm ranges respectively, had not been anticipated upon proposal writing as standard radiative modeling (i.e., assuming LTE) predicted vanishingly small emission levels (0.01 mJy in the strongest lines). The emission in this range must therefore reflect non-LTE fluorescence. Non-LTE effects in planetary atmospheres refer to situations in which excited molecular states are not exclusively populated by collisions, but driven by other external factors (e.g., solar pumping). This results in an over- or under-population of these states, which can be characterized by vibrational temperatures (T_vib) in excess of, or below, the atmosphere local gas (kinetic) temperatures. Methane non-LTE emission at 3.3 μm (ν₃ band complex) has been observed on Jupiter and Saturn (Drossart et al. 1999; Encrenaz et al. 1999; Sánchez-López et al. 2022, and references therein), and Titan (Kim et al. 2000; Geballe et al. 2003; García-Comas et al. 2011). For these objects, observational evidence for non-LTE effects in the dyad region near 6.2–7.7 μm (ν₂/ν₄) is less direct, though significant departures from LTE are expected to occur at 400–500 km (10–1 μbar) on Titan (Vinatier et al. 2007; García-Comas et al. 2011) and (1–0.1) μbar on Jupiter (Sánchez-López et al. 2022; Rodríguez-Ovalle et al. 2024), leading to a reduction of the ν₄ radiances for these objects below the 150–170 K LTE radiation. The above pressure estimates indicate that non-LTE effects may be important throughout Pluto’s atmosphere, boosting in this case the ν₄ radiances above the 100–110 K stratospheric temperature. To our knowledge, CH₃D non-LTE emission at 8.6 μm (ν₆ band) has not been reported before in solar system objects.

The adaptation of, for instance, the Titan CH₄ non-LTE models, to Pluto, as well as the development of non-LTE models for CH₃D, have been left to future works. Here, we used semiempirical models aimed at identifying and fitting the spectral structure observed in the Pluto spectrum. Specifically, we created a “non-LTE” version of our code, in which the Planck function at altitude z in the radiative transfer equation was calculated not at the local kinetic temperature, T_kin(z), but at some T^⋆(z) temperature, and we adjusted the T^⋆(z) profile for each band to fit the overall spectral structure and the radiance level of the emissions. In contrast, the gas kinetic temperature was retained upon calculating the line intensities from the values at T₀ = 296 K given in HITRAN; namely, LTE was assumed for the rotational distribution. Then, T^⋆(z) is related to the vibrational temperature profile T_vib(z) of the upper state of the corresponding band². In the models, we adopted the CH₄ vertical profile from New Horizons (Young et al. 2018; Lavvas et al. 2021), with 0.4% CH₄ at the surface, increasing with altitude due to diffusive separation. For ¹³CH₄, we adopted a telluric ¹²C/¹³C ratio of 89, and assumed the same T_vib(z) as for ¹²CH₄.

The detailed structure of the 7.25–8.25 μm spectrum (see Figs. 2, A.1, A.2) indicates emission from the ν₄ (with separate contributions of ¹²CH₄ and ¹³CH₄), 2ν₄–ν₄, and 3ν₄–2ν₄ bands of CH₄, while all features between 8.25–9.0 μm can be ascribed to the ν₆ band of CH₃D. To account for the observed flux densities, we found it necessary to invoke T_vib temperatures that were significantly enhanced over T_kin down to low altitude levels (Fig. 3). Contribution functions, defined as $$CF(z)=B_{\nu }(T_{vib}(z))\frac{dexp(-{\tau }_{\nu })}{dln(p)}$$(1)

are shown globally for the ν₄ emission in Fig. B.1 and at some illustrative wavelengths in Fig. B.5. They indicate that the non-LTE CH₄ emission probes the atmosphere up to ~450 km altitude and possibly down to the first atmospheric scale heights, depending on the T_vib values there.

Searching for the simplest possible solutions, we first explored isothermal T_vib profiles above 25 km altitude (corresponding to the stratopause in the kinetic temperature profile) at some value T_vib,25. At levels below 25 km, we assumed a T_vib − T_kin difference equal to its value at 25 km. This somewhat arbitrary description aims at representing the fact that near-surface layers remain closer to LTE (i.e., having lower T_vib) than the middle atmosphere does. In doing so, we fine-tuned T_vib,25 separately for the ν₄, 2ν₄ and 3ν₄ levels, until the radiances in the different bands were fit. This approach yielded respective values of T_vib,25 = 125.5 K, 175 K, and 239 K for the ν₄, 2ν₄ and 3ν₄ levels (green dashed, red and dark blue lines in Fig. 3), with a ≲ 1 K uncertainty on these numbers. Therefore, T_vib (3ν₄) ≫ T_vib(2ν₄) ≫ T_vib(ν₄), which is similar to the Titan case (see Fig. 3 in García-Comas et al. 2011). Furthermore, we found that a constant T_vib(ν₄) profile above 25 km did not yield a satisfactory fit of the ν₄ lines, overestimating the emissions in the R branch over 7.5–7.65 μm and underestimating them in the P branch over 7.85–8.05 μm (see the top panel of Fig. A.2).

Exploring vertically variable T_vib(ν₄) profiles and searching again for solutions with minimum complexity, a considerably better agreement was obtained with a T_vib(ν₄) profile decreasing from 129.5 K over 25–200 km to 116 K at 400 km and above (solid green line in Fig. 3). With this profile and including all components (as well as a very weak CH₃D contribution, see below), the overall fit of the 7.2–8.25 μm region is shown in Fig. 2, top panel. Individual band contributions are displayed in Fig. A.1 and their reality further evidenced in Fig. A.2. We note that while overall satisfactory, the model does not completely match the observed flux levels within the noise level. We finally stress that although the inferred T_vib are sensitive to the assumed CH₄ profile, the effect is dwarfed by the extreme sensitivity of fluorescent emission to T_vib. For example, even a factor of two change upward (resp. downward) in the CH₄ mixing ratio can be compensated by T_vib changes of order −1.0 K (+1.5 K) only. Ultimately, when more physically based models of CH₄ fluorescence become available, it will (in principle) be possible to use our observations to constrain Pluto’s CH₄ profile.

The CH₃D mixing ratio³ profile (i.e., the Pluto gas-phase D/H) is currently unknown. To model the CH₃D non-LTE ν₆ emission at 8.6 μm, we assumed that the T_vib profile of the ν₆ level of CH₃D is identical to that of the ν₄ band of CH₄. Under this assumption, which is un-demonstrated but not unreasonable given that these two fundamental bands occur in nearby wavelength ranges, an excellent fit of the 8.25–9.0 μm emissions (Fig. 2, bottom panel) can be obtained with a D/H ratio equal to two times the SMOW telluric value (i.e., with D/H = 3.11 × 10⁻⁴). However, the D/H ratio and the T_vib(CH₃D, ν₆) profile are considerably degenerate, given the extreme temperature sensitivity of the radiance and the fact that the CH₃D emission is mostly optically thin. For example, a virtually identical model can be achieved by assuming a one-time telluric D/H ratio and a T_vib(CH₃D, ν₆) profile only slightly (by 5.5 K) warmer than T_vib(ν₄). Therefore, we are not able to report a definite D/H ratio in Pluto’s atmospheric methane.

In addition to the ν₆, CH₃D exhibits a weak ν₃ band, centered at 1306.8 cm⁻¹ (7.652 μm), in the blue shoulder of the Q-branch of the CH₄ ν₄ band. As shown in Fig. A.2, last panel, this band seems detected in our observations, although with a contrast that is two times smaller than expected based on that of the CH₃D ν₆, assuming the same T_vib profile for the ν₃ and ν₆. Model calculations of this band have therefore been divided by 2 upon showing the individual band contributions in Fig. A.1 and A.2. A modeling of the excitation of this band from first principles will also be warranted.

Finally, we note that JWST NIRSpec spectra of Pluto (Pinilla-Alonso et al. 2023; Pinilla-Alonso et al., in prep.) show fluorescence of CH₄ at 3.3 μm (ν₃ and ν₃+ν₄−ν₄ bands) and of CH₃D at 4.54 μm (ν₂ band). A global excitation model explaining simultaneously all the fluorescent bands of CH₄ and CH₃D will need to be developed to pin-down the D/H ratio in Pluto’s gas methane.

Fig. 2 CH4 and CH3D fluorescence. Observations (in black) are compared to best model fits (in blue). The model includes contributions from the ν4, 2ν4–ν4, and 3ν4–2ν4 bands of CH4, and from the ν6 and ν3 bands of CH3D (see text and Appendix A for details), as well as a solar reflected component with the I/F reflectivity shown in Fig. 1. CH4 and 13CH4 (resp. CH3D) entirely dominate the emission in the upper (lower) panels.

Fig. 3 Pluto’s kinetic temperature profile from New Horizons (Tkin, solid black line) and inferred vibrational temperature profiles for CH4, CH3D, and C2H2, color-coded by band. For CH4, Tvib for the fundamental (ν4) and first two overtone (2ν4 and 3ν4) levels, are shown. The dashed green curve is an alternative Tvib profile for the ν4 level which does not satisfy observations (see Fig. A.2, top panel). The solid green curve (nominal Tvib) for the ν4 CH4 band was also used for calculating the CH3D ν6 emission. The dashed-dotted light blue curve is the retrieved Tvib for the ν5 band of C2H2, assuming the C2H2 Lavvas et al. (2021) vertical profile. The black dashed-dotted and dotted lines are the equatorial and north polar kinetic temperature profiles retrieved by Lellouch et al. (2022) from ALMA 2017 observations.

3.2 Haze

Longwards of 11 μm, continuum emission is progressively affected by thermal emission. Based on the 15–25 μm light curves observed by JWST/MIRI imaging, Bertrand et al. (2025) showed that haze is a significant contributor to Pluto’s emission in this spectral range, superimposing on surface (thermal and reflected) and gas emission⁴. Bertrand et al. (2025) estimated vertically-integrated, filter-averaged haze opacities ranging from a few times 10⁻⁵ at 25 μm to a few times 10⁻⁴ at 15 μm. The present MIRI MRS spectrum now provides spectrally-resolved information on the haze emission.

To recover the haze spectrum, we included haze emission in our LTE radiative transfer code. For that, haze was assumed to extend down to Pluto’s surface with a density (or opacity) scale height H_haze = 50 km (Gladstone et al. 2016; Cheng et al. 2017), and a temperature as a function of altitude identical to that of the gas, an assumption shown to be valid up to 700 km by Zhang et al. (2017). The only free parameter for haze was therefore the vertically-integrated haze optical depth as a function of wavelength, τ_haze(λ). These assumptions enabled us to calculate the haze opacity in each layer, which was added to that of the gas, before radiative transfer was computed. Calculations were done for a series of wavelength-independent values of τ_haze between 0 and 1 × 10⁻³, and the comparison to the observed spectrum in regions unaffected by gas emission yielded the spectral haze optical depth τ_haze from 11 to 22 μm. The τ_haze spectrum was smoothed with a 0.1-μm boxcar and its error bars were estimated from the dispersion of values within the boxcar. Such retrievals were performed with the two surface temperature models. As expected from Fig. 1, results for the two models are virtually identical shortward of ~15 μm and progressively diverge at longer wavelengths. The inferred optical depths are shown in Fig. 4, and compared to the imaginary indices of relevant ices. As anticipated, the order of magnitude τ_haze agrees with that inferred from MIRI imaging. As a consequence of haze optical thinness, results are also very insensitive to the precise haze vertical distribution (i.e., scale height). The spectral structure revealed by MIRI MRS shows at least two distinct peaks in τ_haze, at 686 cm⁻¹ (14.58 μm) and 647 cm⁻¹ (15.45 μm), with respective widths of ~30 and ~50 cm⁻¹. A third peak near 735 cm⁻¹ (13.60 μm) with ~10 cm⁻¹ width is also suggested, while features at wavenumbers >770 cm⁻¹ (λ < 13 μm) are unreliable given the dominance of the solar reflected component there. The three features we here identify are much broader than is characteristic for gas signatures (e.g., ~0.3 cm⁻¹ for individual C₂H₆ lines – defined by the instrument’s resolution, ~1 cm⁻¹ for the Q-branch of C₄H₂ and ~2–3 cm⁻¹ for the strong Q-branch of C₂H₂, see hereafter). They can therefore be safely ascribed to solid-state features.

Fig. 4 Haze optical depth (multiplied by 103), compared to imaginary refractive indices of several plausible icy compounds at relevant temperatures. The solid thick line with 2-σ error bars corresponds to the haze optical depth inferred by using the CH4 seasonal surface temperature model. The thin solid line, which diverges from the previous one longward of 15 μm, makes use of the diurnal-only surface temperature model. References are Moore et al. (2010) for HCN, Ugelow & Anderson (2022) for HC3N, Hudson & Yarnall (2022) for C6H6 and Khanna et al. (1988) for C4H2. Also shown are the imaginary refractive indices of Titan’s haze (Vinatier et al. 2012) and of Titan’s tholin analogues (Imanaka et al. 2012). The retrieved haze optical depth is unreliable shortwards of ~13 μm (thin dashed line).

Fig. 5 Ethane emission over 11.7–12.8 μm. Black: observations. Red: model using a C2H6 profile equal to 1.05 times the Lavvas et al. (2021) profile and a telluric C2H6/13C12CH6 ratio. The inset shows a zoom on the 13C12CH6 lines, and includes models in which the C2H6/13C12CH6 ratio is multiplied (blue curve) or divided (green curve) by 2.

Fig. 6 Methylacetylene (CH3C2H) and diacetylene (C4H2) emission over 15.4–16.2 μm. Black: observations. Red: model using CH3C2H and C4H2 profiles following the Lavvas et al. (2021) profiles but divided by factors of 9.5 and 4.5, respectively.

3.3 Gas composition

With the surface and haze model established, we turn to modeling of the minor gas emission longwards of 12 μm. We nominally assumed that LTE prevails for these emissions, although, as discussed later, the case for non-LTE emission from C₂H₂ was also envisaged and actually favored.

3.3.1 Ethane, methylacetylene and diacetylene

The MIRI MRS spectrum shows prominent emission from the ν₉ band of ethane (C₂H₆), now routinely observed in the spectra of the four Giant Planets and Titan. The comparison to radiative transfer calculations indicates this emission can be fit with a vertically uniform ethane mixing ratio of 1.60 × 10⁻⁵. Using the C₂H₆ profile from the photochemical model of Lavvas et al. (2021), we found that the latter profile needs to be rescaled by a mere factor of 1.05±0.1 to match the observation (Fig. 5), where the error bar is derived from fit fine-tuning in different parts of the spectrum. With such a mixing ratio profile, the C₂H₆ band is essentially optically thin (maximum optical depth = 0.40 at infinite spectral resolution for a ray intersecting the disk at 0.5 R_Pluto). Contribution functions, shown in Fig. B.2, indicate that the C₂H₆ band mostly probes the stratopause region – as a result of the highest temperatures there – with contribution functions at all wavelengths being centered in a layer located at ${30}_{-12}^{+33}$ km (at half-maximum). Therefore, the spectrum does not yield vertically-resolved information on the C₂H₆ vertical profile, and we best determine a (1.35±0.15) × 10⁻⁵ ethane mixing ratio at 30 km altitude. The model shown in Fig. 5 includes the ¹³C¹²CH₆ isotope in natural abundance. Pluto’s ethane spectrum indicates a likely detection of ¹³C¹²CH₆ lines (set inset), with an abundance consistent with the terrestrial ¹²C/¹³C ratio, but this determination is probably accurate to within a factor of ~2 only.

Figure 6 shows the detection of methylacetylene (CH₃C₂H) and diacetylene (C₄H₂) emission over 15.4–16.2 μm, where the 15.92 μm peak is due to C₄H₂ while CH₃C₂H accounts for all other features. Like for ethane, these features are well known in the spectra of outer planets, but are here observed at Pluto for the first time. Assuming uniform vertical distributions, they can be fit for CH₃C₂H and C₄H₂ mixing ratios of 5.5 × 10⁻⁸ and 5.2 × 10⁻⁹, respectively. In addition, we find that the observed emissions are much weaker than anticipated based on the Lavvas et al. (2021) model, and indicate abundances (9.5±1.0) and (4.5±0.5) times smaller, for CH₃C₂H and C₄H₂ respectively, than in the model. This corresponds to column densities⁵ of (3.7±0.4) × 10¹⁴ and (2.1±0.2) × 10¹³ cm⁻², for CH₃C₂H and C₄H₂ respectively. Again, these emissions are optically thin (maximum optical depth = 0.038 and 0.018 for CH₃C₂H and C₄H₂ respectively), preventing vertical information from being retrieved, as confirmed by examining the contribution functions (Fig. B.2). Considering again the half-maximum of the contribution functions as indicative of the probed region, we determine q(C₄H₂) = (1.1±0.1) × 10⁻⁸ at ${30}_{-7}^{+14}$ km and q(CH₃C₂H) = (4.6±0.5) × 10⁻⁸ at ${55}_{-30}^{+80}$ km.

All the above numbers and error bars are nominally based on the New Horizons thermal profile and the associated Lavvas et al. (2021) mixing ratio profiles and do not include the effect of temperature uncertainties. For illustrative purposes, we refitted the MIRI MRS spectrum with the disk-averaged temperature profiles for 2015 and 2017 presented in Lellouch et al. (2017) and Lellouch et al. (2022), respectively, which generally bracket the New Horizons profile in the broad stratopause region of interest. Considering these profiles leads to abundance changes of C₂H₆, CH₃C₂H, and C₄H₂ by ±25%, ±30%, and ±10%, respectively, from their nominal values.

Fig. 7 Acetylene (C2H2) and its isotope emission over 12.8–14.8 μm. Black: observations. Colored curves: models. All models assume a terrestrial 12C/13C ratio. Top: LTE models (i.e,. using the Tkin profile) for various C2H2 profiles. Light blue: reference profile from Lavvas et al. (2021); red: vertically uniform C2H2 with 4 × 10−6 mixing ratio; green: slope profile from Fig. 8. Residuals for this latter model are shown. Bottom: Non-LTE models (adjusted Tvib) using the C2H2 profile from Lavvas et al. (2021). Green curve: nominal solution, using the dashed-dotted light blue Tvib(z) profile in Fig. 3. Observation minus model residuals are shown for this model, as well as for models using Tvib(z) + 1.5 K (dark blue) and Tvib(z) – 1.5 K (red). Insets: Zoom on the region of the ν5 band of C2HD at 14.75 μm. Models use either the C2H2 slope profile (top) or the Tvib solution profile (bottom), and D/H ratios in C2H2 equal to one, three, and five times the telluric value.

3.3.2 Acetylene and mono-deuterated acetylene

The 12.8–14.8 μm range of Pluto’s spectrum is dominated by emission from acetylene C₂H₂ and its isotopic variant ¹²C¹³CH₂. To zeroth order, the spectral structure is reminiscent of Titan’s spectrum (Coustenis et al. 2003), but (in addition to the lack of C₃H₈ and HCN signatures here, see below), a significant difference is that the ¹²C¹³CH₂ Q-branch at 13.65 μm, and its P and R branches embedded with those of ¹²C¹²CH₂ over 13.1–14.1 μm are much more evident in Pluto than in Titan.

Unlike the case of ethane, the nominal acetylene profile from Lavvas et al. (2021) leads to a marked underestimate (by essentially a factor of two) of the observed 12.9–14.7 μm radiances (blue curve in Fig. 7, top). In fact we were unable to achieve a fully satisfactory fit of the entire C₂H₂ emission. The mean (assumed vertically constant) C₂H₂ indicated by the data in a least-square sense is 4 × 10⁻⁶, but the fit is very poor, with the Q-branch strongly overestimated (red curve in Fig. 7, top). The strongest C₂H₂ lines are optically thick (maximum optical depth = 75), indicating that vertically-resolved information can be retrieved. Exploring altitude-varying profiles (with constant d(log q)/dz gradient over 20–400 km, and C₂H₂ constant above 400 km), we obtained a best fit for d(log q)/dz = 0.0064±0.0004 km⁻¹ and q_{100 km} = (5.8±0.4) × 10⁻⁶ (profile shown as the thin solid blue lines with data points in Fig. 8). The profile is less steep than the Lavvas et al. (2021) profile over 50–200 km. Although seemingly paradoxical, the fact that it produces a larger flux than the latter is related to the negative mesospheric temperature gradient, which leads to stronger absorption cores within the emission lines for the Lavvas et al. (2021) profile, diminishing the overall emission contrast. Contribution functions (see Fig. B.4) indicate that the 20–200 km altitude range is probed. However the model underestimates the high J lines, especially near 12.9–13.3 μm (green curve and residuals in Fig. 7, top).

Difficulties in fitting the C₂H₂ ν₅ band down to observational errors have been encountered previously in Titan (Coustenis et al. 2003, their Fig. 9), and in Jupiter’s auroral regions, both in high-resolution ground-based spectra (Sinclair et al. 2018, their Fig. 16) and in JWST MIRI data (Rodríguez-Ovalle et al. 2024, their Fig. 14), where models underestimate the emission cores. In Jupiter, the problem is tentatively attributed to non-LTE effects in the low pressure (<10⁻² mbar) regions probed by the C₂H₂ strong lines. This explanation is not obviously viable at Jupiter, because, as pointed out by Rodríguez-Ovalle et al. (2024), non-LTE effects there are likely to lead to vibrational temperatures lower than the kinetic temperature; namely, this acts to decrease the expected radiances below the LTE case. In Pluto, the explanation may be more plausible, because solar excitation may enhance the upper atmosphere T_vib above T_kin, as readily demonstrated by the strong CH₄ 7.7 μm emission (Sect. 3.1). For C₂H₂, assuming that solar excitation can occur through the ν₃ band near 3.0 μm, we calculated a fluorescence equilibrium temperature⁶ T_eq.fluo = 74 K. This result is indeed a few degrees above Pluto’s upper atmosphere temperature, which is in the range of 65–70 K (see Fig. 3).

In an alternative approach, we therefore attempted to fit the acetylene emission by using the nominal C₂H₂ profile from Lavvas et al. (2021), this time adjusting the T_vib profile as a function of altitude, as we did previously for methane. The best fit T_vib profile in a least-square sense is shown in Fig. 3, dashed-dotted blue curve. These T_vib are conservatively determined to within ±1.5 K, although this approach does not improve the overall fit, as problems in fitting simultaneously the Q-branch and the distant P, R lines remain (Fig. 7, bottom). We are therefore left with two imperfect “end-member” models, one calling for a shallower variation of C₂H₂ over 30–200 km compared to Lavvas et al. (2021), and one calling for non-LTE fluorescent excitation of the C₂H₂ ν₅ band, leading to T_vib temperatures in excess of T_kin. Because the Lavvas et al. (2021) model yields a good fit of the C₂H₆ and C₂H₂ LOS column densities measured by New Horizons/Alice (Young et al. 2018; Lavvas et al. 2021, their Fig. 5), while the slope model underestimates the C₂H₂ mixing ratio gradient over 100–200 km compared to the observational profiles of Young et al. (2018, in particular, their Fig. 21), we at this point tentatively favor the non-LTE model.

The T_vib temperatures for C₂H₂ are consistently warmer than the adopted T_kin profile, by 4–6 K below 200 km, and up to ~10 K at 400 km. They are also warmer than the disk-averaged 2015 and 2017 profiles from ALMA (Lellouch et al. 2017, 2022), and warmer than the temperatures at all latitudes in the Lellouch et al. (2022) thermal field (their Fig. 8; see also Fig. 3 above). Under LTE, these various temperature profiles would underpredict the ν₅ Q-branch of C₂H₂ by at least 40%. Whether the enhanced T_vib for C₂H₂ relative to the New Horizons T_kin profile might actually represent the atmospheric temperature profile in 2023 is doubtful and will require further study. Current general circulation models (Toigo et al. 2015; Forget et al. 2017) predict that variability of Pluto’s temperature profile over a decade timescale is restricted to a few K or less, and that latitudinal variations are similarly small (<1 K), except at the winter pole. However, these predictions hold for a gas-only radiative control of the atmosphere, which we now know is not the case (Bertrand et al. 2025). Furthermore, the 59.6° N sub-observer latitude in the MIRI observations preferentially probes high spring/summer latitudes, tentatively found to be warmer than the rest of the atmosphere in 2017 (Lellouch et al. 2022), though below 50 km altitude only. Progress on this complex issue will require (i) developing physically based NLTE models for C₂H₂; (ii) extending the GCM models to include haze radiative effects; and (iii) remeasuring Pluto’s thermal field with ALMA to search for posssible variations since 2017.

The long-wavelength part of the C₂H₂ ν₅ region shows a distinct, spectrally resolved, feature at 14.75 μm (see insets in Fig. 7), which is due to the ν₅ band of C₂HD, previously observed only in Titan (Coustenis et al. 2008). From comparison with synthetic spectra in which the D/H ratio in C₂H₂ is varied (see insets in Fig. 7), we infer D/H (in C₂H₂) = (3.5±0.5) × terrestrial for the C₂H₂ slope models and (2.5±0.5) × terrestrial for the enhanced T_vib models (assuming T_vib to be the same for C₂H₂ and C₂HD). We conservatively conclude that deuterium is enriched in Pluto’s C₂H₂, by a factor of (3±1) over the telluric value (SMOW, 1.56 × 10⁻⁴).

Fig. 8 Volume mixing ratio profiles from Lavvas et al. (2021), compared to inferences from this work. Arrows indicate upper limits.

3.4 Other compounds

We also searched (unsuccessfully, however) for signatures of other atmospheric compounds previously detected (C₂H₄, HCN; Young et al. 2018; Lellouch et al. 2017) in Pluto’s atmosphere or expected to be present in some amounts from photochemical models (C₃H₈, CO₂, HC₃N, C₆H₆; e.g., Lavvas et al. 2021). The corresponding spectral regions are shown in Fig. C.1, and the upper limits on the abundances at the relevant altitudes are reported in Table 1 and Fig. 8. HCN is searched for through its ν₂ band Q-branch at 14.04 μm, embedded within the P-branch multiplets of C₂H₂ and ¹²C¹³CH₂. A feature is actually observed there and can be fit with a HCN profile equal to twice the modeled one. However, given the complications associated with the C₂H₂ fit (see Sect. 3.3.2), with significant fit residuals in this spectral region, and the fact that the HCN feature superimposes to one of the ¹²C¹³CH₂ multiplets, we view this detection of HCN as tentative and prefer to regard it as an upper limit. Upper limits on HC₃N (from the lack of a 15.07 μm band) and C₂H₄ (at 10.53 and 10.85 μm) are two and four times the Lavvas et al. (2021) profile. The more constraining upper limit we obtained is on propane (C₃H₈), which, based on the absence of the ν₂₁ band at 13.37 μm, is found to be at least five times less abundant than in the model. Finally, upper limits on C₆H₆ and CO₂, based on the lack of signatures at 14.84 and 14.98 μm respectively, are rather unconstraining (50 and 500 times, respectively, the model profile) and not reported in Fig. 8. Abundance results or upper limits on the various species, in terms of mixing ratio, multiplicative factor from the Lavvas et al. (2021) profiles, and associated column densities are gathered in Table 1. For HC₃N, the upper limit (3.0×10¹³ cm⁻²) is coincidentally equal to that previously obtained from ALMA (Lellouch et al. 2022).

4 Implications for Pluto’s chemistry and haze

4.1 Gas photochemistry

Due to the strong temperature dependence of the Planck function in the MIR, and the fact that, except for C₂H₂, the emissions are optically thin, the observed spectrum primarily samples the warmest atmospheric region (>95 K) at 15–100 km altitude. As summarized in Sect. 3.3 and Table 1, our results for C₂ hydrocarbons indicate a generally good agreement with the New Horizons/Alice measurements and with the Lavvas et al. (2021) model. For C₂H₆, the Lavvas et al. (2021) model already yielded a very good fit of the LOS column densities observed by New Horizons/Alice down to 40 km altitude (see their Fig. 5), demonstrating the role of heterogeneous chemistry and desorption of C₂H₆ gas on the surface of aerosols. Here, we determine a measurement of the C₂H₆ mixing ratio near 30 km which confirms these findings. As explained before, the fit of the C₂H₂ emission is imperfect and ambiguous, so at this point our results on C₂H₂ cannot be used to further constrain models. For C₂H₄, our upper limit is comfortably consistent (<4 times) with the Lavvas et al. (2021) profile, but the model profile shows discrepancies (at the factor of ~2 level) with the UV-measured LOS column densities. The comparison with Fig. 21 of Young et al. (2018) indicates that our new constraint (C₂H₄ < 3.4 × 10⁻⁷ near 30 km) is consistent with the New Horizons/Alice C₂H₄ profile, and that a detection of C₂H₄ by JWST MIRI was probably not far from being achieved.

Our most novel results are for the C₃ and C₄ hydrocarbons, with a clear detection of methylacetylene (CH₃C₂H) and diacetylene (C₄H₂) and a stringent upper limit on propane (C₃H₈), with abundances that are 9.5 times, 4.5 times, and at least 5 times lower than in the Lavvas et al. (2021) model. For CH₃C₂H, the derived column density (3.7 × 10¹⁴ cm⁻²) is consistent with the upper limit (<8.5 × 10¹⁴ cm⁻²) obtained in the mm range from ALMA (Lellouch et al. 2022). In contrast, it is sharply at odds (a factor of seven less) with the 2.5 × 10¹⁵ cm⁻² column density inferred by Steffl et al. (2020) from New Horizons/Alice⁷. Our clear detection of CH₃C₂H confirms the premise in Lellouch et al. (2022) that the claimed identification of Steffl et al. (2020) was spurious; in fact the 1535–1550 Å feature in the Alice spectrum attributed to CH₃C₂H by these authors may be dominated by CO absorption (Tommi Koskinen, priv. comm, 2023).

A possible explanation for the lower abundances of CH₃C₂H, C₄H₂ and C₃H₈ over 15–100 km compared to the Lavvas et al. (2021) model, already put forward for CH₃C₂H by Lellouch et al. (2022), is that the model overestimates the re-sublimation of the corresponding ices in the lower atmosphere. The initial ALMA observations (Lellouch et al. 2017), near concurrent with the New Horizons flyby, revealed the presence of HCN in Pluto’s atmosphere. The satellite (hyperfine structure) emissions of the main J = 4–3 HCN line were sensitive to the 50–150 km altitude range; to reproduce them a release of HCN from the condensed phase of the organic ice haze particles was necessary. This re-sublimation process was treated as a free parameter in Lavvas et al. (2021). Subsequent ALMA observations obtained in 2017 (Lellouch et al. 2022) revealed reduced HCN satellite emissions, thus requiring a smaller re-sublimation contribution. As the re-sublimation affects all condensed ices in the haze particles, reducing its contribution leads to smaller gas phase abundances. For CH₃C₂H, C₄H₂, and C₃H₈, turning off re-sublimation in the model reduces their column abundances by factors of 24, 1.4, and 4. This modification brings the simulated profiles closer to the retrieved abundance constraints (Table 1), but a discrepancy remains for C₄H₂, which continues to be overpredicted in the model by a factor of ~3. We note that the sublimation law for C₄H₂ in the relevant temperature range is very uncertain, as measurements are restricted to 127–152 K and reliable only within a factor of 2.5–3 (see Fray & Schmitt 2009).

A face value comparison with the models of Krasnopolsky (2020a,b) can also be performed. For CH₃C₂H in the relevant altitude range (20–135 km), Krasnopolsky (2020a) predicts a ~10⁻⁷ mixing ratio, 50 times less than observed, but the updated model in Krasnopolsky (2020b) has a mean CH₃C₂H of 3 × 10⁻⁶ mixing ratio, closer to observations. For C₄H₂ at 25–45 km, the Krasnopolsky (2020a) model value is ~2 × 10⁻⁸, within a factor of two to the observed; for C₃H₈ at 15–130 km, the model value is ~5 × 10⁻⁷, consistent with the 1.2 × 10⁻⁶ upper limit. This comparison is rough because the vertical profiles in Krasnopolsky (2020a,b) are very different (much steeper over 0–200 km) from those of Lavvas et al. (2021) that we used as templates, and a more definite comparison would involve testing their actual profiles (their Fig. 8 and Fig. 1, respectively).

4.2 Haze composition

The initial post New Horizons view that Pluto’s haze might be primarily composed of Titan-like organic aerosols (e.g., Wong et al. 2017; Gao et al. 2017; Zhang et al. 2017) was challenged by Lavvas et al. (2021). Based on their photochemical-microphysical modeling, these authors argued that Pluto’s haze might include a large hydrocarbon and nitrile ice component, with C₄H₂ ice as the prime contributor, followed by HCN and C₆H₆ ices. Figure 4 shows the optical depth spectrum of Pluto’s haze, compared to imaginary refractive indices of several relevant ices (C₄H₂, C₆H₆, HCN, HC₃N) and to those inferred for Titan from Cassini/CIRS (Vinatier et al. 2012) and for laboratory tholin analogues (Imanaka et al. 2012).

The FIR spectrum of Titan’s haze shows features at 630, 700, 750, and 800–1000 cm⁻¹, of variable widths, the most prominent of which being attributed to C–H bending vibrations in methyl and methylene groups (see a more detailed discussion in Vinatier et al. 2012). Pluto’s haze spectrum does not show such features, being instead characterized by narrower peaks at 686 cm⁻¹ (14.58 μm), 647 cm⁻¹ (15.45 μm), and possibly 735 cm⁻¹ (13.60 μm). Interestingly, the first one falls close to the 676–696 cm⁻¹ emission in Titan’s southern stratospheric polar cloud observed by Cassini/CIRS in nadir and limb mode during southern autumn (Vinatier et al. 2018). This emission, with an asymmetric maximum at 685–690 cm⁻¹, depending on the altitude and season, was at least partly attributed to benzene ice, whose ν₄ C-H asymmetric bending mode falls at 682 cm⁻¹, with small particle size (<1.5 μm radius). As shown in Fig. 4, C₆H₆ ice could also contribute to Pluto’s haze 686 cm⁻¹ emission, but only to some extent as the match in wavelength is not perfect⁸, and as C₆H₆ ice also exhibits a weaker band near 706 cm⁻¹ (ν₈ ring deformation, see also Mouzay et al. 2021a), unobserved in the JWST/MIRI spectrum. Similarly the broader 647 cm⁻¹ feature in Pluto’s haze spectrum is close to – but not quite at the position of – the ν₄ C≡C-H band of C₄H₂ ice (Khanna et al. 1988), and the latter actually shows a double peak structure (653/663 cm⁻¹) not present in the data.

As a further check of the plausibility of such assignments, we made simple estimates of the associated mass loading. Fan et al. (2022) inferred that Pluto’s haze particle size distribution is bimodal, consisting of the combination of ~ 1-μm radius two- dimensional aggregates (with 20-nm size monomers) and a smaller-sized population made of ~80 nm spheres, with similar relative contributions to the total mass. We performed Mie calculations of the extinction coefficient (Q_ext) for these two populations, considering the 15.3 μm peak of the C₄H₂ ice feature (where the real and imaginary indices are n = 1.5 and k = 0.8, see Khanna et al. 1988). For 1-μm and 80 nm radius particles, Q_ext ~0.6 and ~0.045, and the observed optical depth (~4.5 × 10⁻⁴), along with a scale height of 50 km, implies a nearsurface C₄H₂ ice number density of ~0.004 cm⁻³ and ~8 cm⁻³ for the two particle sizes respectively. These numbers compare reassuringly well with (and bracket) the estimated overall haze number density (0.8–0.3 cm⁻³ for the aggregates and 10 cm⁻³ for the spheres; Gladstone et al. 2016; Fan et al. 2022). We also note that the imaginary refractive index at the C₄H₂ peak is 2.5 times smaller than at the C₆H₆ peak. This typically compensates for the ~4 times larger contribution of C₄H₂ ice to the haze mass flux compared to C₆H₆ ice in the Lavvas et al. (2021) model, indicating that the C₄H₂ and C₆H₆ ice features should be of comparable contrast or area, which is indeed the case in the observed spectrum.

Given that the assignment of solid state and ice features in the thermal IR has typically proven difficult⁹, and that other plausible species such as HCN and HC₃N ice do not match the observations, we provisionally suggest that Pluto’s mid-IR haze features have more resemblance to ice spectra (including C₄H₂ and C₆H₆ ice) than to Titan or Titan-like tholins. While this generally gives reasonable support to the Lavvas et al. (2021) model, we recommend laboratory measurements of the mid- to far-IR spectrum of binary (or more) mixtures ices (see e.g., Mouzay et al. 2021b, for C₆H₆-HCN mixtures) to strengthen the identification.

Finally, we note that based on New Horizons/LORRI phase curves in four visible/near-IR filters, Hillier et al. (2021) attempted to recover the photometric properties of Pluto’s haze and surface, and proposed that, in terms of single-scattering albedo (see their Fig. 6), Pluto’s haze is more similar to Titan’s organic haze than to Triton’s icy haze. As these authors recognize, however, for an optically thin haze, there is a full degeneracy between optical depth and single scattering albedo; furthermore their analysis did not take into account the aggregate nature of the particles. Therefore, we feel that the spectroscopic evidence from JWST is much more diagnostic of the composition of Pluto’s aerosols. In any case, our study does not directly constrain the optical properties of the haze in the visible/near-IR range.

4.3 D/H ratio

C₂H₂ is a product of the photochemistry of CH₄. The latter is the main reservoir of atmospheric deuterium, and its gas phase abundance is defined by seasonal sublimation/condensation cycles and escape. Therefore, the (D/H)_C2H2 ratio results from the combination of: (i) the unknown D/H ratio in Pluto’s CH₄ ice, itself diagnostic of the origin and evolution of methane (Glein et al. 2024); (ii) possible fractionation effects on atmospheric methane during sublimation and condensation (at low temperatures, the vapor pressure of CH₃D is lower than that of CH₄; Calado et al. 1997) and as a result of escape (CH₄ and H escape more readily than CH₃D and D due to their smaller masses); and (iii) kinetic and photolytic-induced fractionation effects, resulting from the higher energy of the C-D bond compared with the C-H bond. The latter effect tends to increase the atmospheric (D/H)_CH4 over time (see Pinto et al. 1986; Lunine et al. 1999, for the Titan case) and to enhance (D/H)_C2H2 with respect to (D/H)_CH4 in the atmosphere.

Using JWST/NIRSpec, Grundy et al. (2024) obtained the first measurement of (D/H) in the CH₄ ice of Eris and Makemake (1.60±0.32 and 1.86±0.38 times SMOW, respectively); this mild enrichment, as opposed to much larger (D/H) values in methane or refractory organics in comet 67P/Churyumov-Gerasimenko, indicates an endogenic (e.g., from hydrothermal activity, metamorphic reactions involving organic matter, etc.), rather than a primordial, origin of methane in these two TNOs (Glein et al. 2024). A similar situation may hold on Pluto, whereby interior models suggest possible high temperatures and the potential for a subsurface ocean (Bierson et al. 2020; Nimmo & McKinnon 2021). Hence, a (D/H) ratio in Pluto’s CH₄ ice of ~2 × SMOW seems plausible.

At Titan, the atmospheric (D/H)_C2H2 was recently shown to be fully consistent with (D/H)_CH4 (Bézard et al. 2024), indicating no significant fractionation of deuterium in acetylene compared with methane, with implications for the kinetic isotope effect (KIE) in photolytical and chemical routes to C₂H₂. Photochemical models of Pluto need to be extended to account for isotopes, but we can provisionally assume that a similar conclusion holds for Pluto’s atmosphere, namely, (D/H)_C2H2 ~ (D/H)_CH4.

We are not yet able to report a definite D/H ratio in Pluto’s atmospheric CH₄. In contrast, we believe that our D/H value in C₂H₂, (D/H)_C2H2 = (3±1) × SMOW, is robust, despite the imperfect fit of the entire C₂H₂ ν₅ band and the possible non-LTE effects on that band. Given the large error bar, this result is consistent with a scenario where the D/H in CH₄ ice, CH₄ gas and C₂H₂ gas are all equal; however, it alternatively might favor the case for deuterium fractionation occurring between the surface and the atmosphere during sublimation-condensation cycles, enriching the atmosphere in D/H in both (and perhaps equally) CH₄ and C₂H₂¹⁰. Such surface fractionation would combine with atmospheric escape to set an enhanced D/H in the gas versus the ice phase.

5 Summary and prospects

Deep observations of the Pluto-Charon system, recorded in May 2023 with JWST’s MIRI MRS instrument over 4.9–27 μm, have yielded a high S/N spectrum of Pluto, opening up a new spectral range for the study of Pluto’s atmospheric gas and haze composition. Based on radiative transfer modeling of these data, the main results of this study are as follows.

• Pluto’s mid-infrared (MIR) spectrum consists of the combination of solar light reflected off Pluto’s surface, gas thermal and non-thermal emission, and haze emission.

• The reflected solar component shows signatures of CH₄, CH₃D, and C₂H₄ ices at 5–8.5, 8.6, and 10.4 μm, respectively.

• Pluto’s haze far-infrared (FIR) emission, first evidenced by findings in photometric MIRI imaging data (Bertrand et al. 2025), is spectrally characterized over 13–20 μm, presenting emission peaks at 15.45, 14.58, and possibly 13.60 μm. Although identifications are qualitative at this point, this points to a significant contribution of pure or mixed ices (e.g., C₄H₂, C₆H₆, and HCN) to the haze composition, supporting a recently proposed scenario (Lavvas et al. 2021) and in contrast to the composition of Titan’s haze.

• Unexpectedly strong emissions due to CH₄ and CH₃D gases are detected at 7.7 and 8.6 μm, which appear to be caused by non-LTE fluorescent emission in the ν₄, 2ν₄–ν₄, and 3ν₄–2ν₄ bands of CH₄ as well as the ν₆ (and possibly ν₃) bands of CH₃D. The vibrational temperatures for the CH₄ states are typically 115–130 K for ν₄ (altitude-dependent), 175 K for 2ν₄, and 240 K for 3ν₄, significantly above the 65–110 K kinetic temperatures in Pluto’s stratosphere/mesosphere.

• Emissions due to C₂H₆, C₂H₂, C₂HD, CH₃C₂H, and C₄H₂ gases are prominently detected over 12–16 μm, yielding the first unambiguous detection of the latter three species in Pluto’s atmosphere. They broadly sound the stratopause region (altogether 15–100 km). Assuming LTE emission conditions, the inferred C₂H₆ abundance is in excellent agreement with predictions based on the photochemicalmicrophysical model of Lavvas et al. (2021), but the CH₃C₂H and C₄H₂ abundances are a factor of about five and ten times smaller, respectively, than in those models. This suggests that the resublimation of the corresponding ices in the lower atmosphere is overestimated by the model. A similar conclusion may be reached based on the non-detection of C₃H₈.

• The C₂H₂ emission is about two times stronger than anticipated. This may result from either (i) a shallower slope of the C₂H₂ vertical profile over 30–200 km compared to models or, preferably, (ii) non-LTE effects on C₂H₂, leading to vibrational temperatures in excess of kinetic temperatures by 4–10 K. In both scenarios, the C₂H₂ emission is imperfectly reproduced by our radiative transfer models, but the unambiguous C₂HD emission at 14.75 μm indicates a (D/H)_C2H2 ratio of (3±1) x SMOW, which likely results from the combination of an enhanced D/H ratio in Pluto’s CH₄ ice, as well as additional fractionation effects between the surface and atmosphere.

As a take-away conclusion, these “Voyager-like” JWST/MIRI-MRS data yield a new perspective on Pluto’s surface and atmospheric gas and haze composition, and on the importance of atmospheric non-LTE effects. Further progress in analysis, particularly in relation to the D/H ratio, will include (i) a modeling of the CH₃D and CH₄ ice features in the MIRI (Souza Feliciano et al., in prep.) and NIRSpec (Pinilla-Alonso et al., in prep.) spectra; and (ii) the development of a full-blown, non-LTE model, for CH₄, CH₃D, and C₂H₂ fluorescence in the MIRI and NIRSpec ranges. This should ultimately yield the D/H ratio in three different reservoirs, which features a unique case for a planetary body beyond Earth. From the theoretical point of view, similarly to recent modeling developments for Mars (Vals et al. 2022), the deuterium cycle should be included in Pluto GCM and climate models, as, for instance, described by Forget et al. (2017); Bertrand et al. (2020). In another avenue, the direct measurement of the gas and haze thermal emission from Pluto afforded by the MIRI MRS spectrum will provide strong constraints on the heat budget of Pluto’s atmosphere, warranting a revisiting of earlier models (Yelle & Lunine 1989; Strobel et al. 1996; Strobel & Zhu 2017; Zhang et al. 2017; Wan et al. 2021).

Acknowledgements

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 #1658. NPA acknowledges the Ministry of Science, Innovation, 118 and Universities (MICIU) in Spain and the State Agency for Research (AEI) for 119 funding through the ATRAE program, project ATR2023-145683. We thank Manuel López-Puertas for useful discussions.

Appendix A Details on fluorescence modeling

Fig. A.1 Components of the fluorescence in the 7.2–8.2 μm range: ν4 band, shown separately for 12CH4 and 13CH4, 2ν4–ν4, and 3ν4–2ν4 bands. The CH3D contribution (ν3 band), calculated with the same model as for the ν6 band, has been divided by 2 (see text).

Fig. A.2 Model sensitivity for the CH4 fluorescence in some selected spectral intervals. Top panel: sensitivity to Tvib (ν4) profile (non-isothermal above 25 km vs isothermal at 125.5 K). Other panels: sensitivity to model omitting individual fluorescence contributions one-by-one: 13CH4, 2ν4–ν4, 3ν4–2ν4, CH3D ν3. In each panel, the nominal model including all contributions is shown in blue, and the alternative model in red.

Appendix B Contribution functions

Fig. B.1 Contribution functions for the ν4 non-LTE emission of CH4. They are calculated for the Tvib(z) solution profile (solid green line in Fig. 3), a LOS intersecting the Pluto disk at r = 0.5 RPluto, and convolved to the instrumental resolution.

Fig. B.2 Contribution functions for the ν9 band of C2H6. They are calculated for the solution C2H6 abundance and a LOS intersecting the Pluto disk at r = 0.5 RPluto, and convolved to the instrumental resolution.

Fig. B.3 Contribution functions for CH3C2H and C4H2. They are calculated for the solution CH3C2H and C4H2 abundances, a LOS intersecting the Pluto disk at r = 0.5 RPluto, and convolved to the instrumental resolution.

Fig. B.4 Contribution functions for the ν5 band of C2H2, in the “non-LTE” solution case. They are calculated for the Tvib(z) solution profile (dashed-dotted light blue line in Fig. 3), the nominal C2H2 profile from Lavvas et al. (2021), and a LOS intersecting the Pluto disk at r = 0.5 RPluto. and convolved to the instrumental resolution

Fig. B.5 Contribution functions at selected wavenumbers for the different atmospheric species, calculated as described in Fig. B.1–B.4.

Appendix C Upper limits

Fig. C.1 Tentative detection of HCN and upper limits on C2H4, C3H8, C6H6, CO2, and HC3N. For each species, the number in parentheses indicates the factor by which the nominal profile from Lavvas et al. (2021) has been multiplied. In all panels, model excluding (resp. including) the considered species is shown in blue (resp. in red). In the first panel (HCN), observation minus model residuals are also shown.

References

All Tables

All Figures

Fig. 1 Overview of the MIRI MRS spectrum of Pluto, with atmospheric gas features identified (black). The light blue points show the Pluto fluxes from MIRI imaging photometry at a Pluto longitude of 342° E (Bertrand et al. 2025). Red curves indicate a model without atmospheric features. The model includes a solar component (blue) and a thermal component (green). Solid (resp. dashed) lines correspond to the diurnal-only (resp. CH4 seasonal) case for surface temperature. The I/F reflectivity for the solar reflected component is shown in the inset. It indicates absorptions attributed to CH4, CH3D, and C2H4 ices.

In the text

Fig. 2 CH4 and CH3D fluorescence. Observations (in black) are compared to best model fits (in blue). The model includes contributions from the ν4, 2ν4–ν4, and 3ν4–2ν4 bands of CH4, and from the ν6 and ν3 bands of CH3D (see text and Appendix A for details), as well as a solar reflected component with the I/F reflectivity shown in Fig. 1. CH4 and 13CH4 (resp. CH3D) entirely dominate the emission in the upper (lower) panels.

In the text

Fig. 3 Pluto’s kinetic temperature profile from New Horizons (Tkin, solid black line) and inferred vibrational temperature profiles for CH4, CH3D, and C2H2, color-coded by band. For CH4, Tvib for the fundamental (ν4) and first two overtone (2ν4 and 3ν4) levels, are shown. The dashed green curve is an alternative Tvib profile for the ν4 level which does not satisfy observations (see Fig. A.2, top panel). The solid green curve (nominal Tvib) for the ν4 CH4 band was also used for calculating the CH3D ν6 emission. The dashed-dotted light blue curve is the retrieved Tvib for the ν5 band of C2H2, assuming the C2H2 Lavvas et al. (2021) vertical profile. The black dashed-dotted and dotted lines are the equatorial and north polar kinetic temperature profiles retrieved by Lellouch et al. (2022) from ALMA 2017 observations.

In the text

Fig. 4 Haze optical depth (multiplied by 103), compared to imaginary refractive indices of several plausible icy compounds at relevant temperatures. The solid thick line with 2-σ error bars corresponds to the haze optical depth inferred by using the CH4 seasonal surface temperature model. The thin solid line, which diverges from the previous one longward of 15 μm, makes use of the diurnal-only surface temperature model. References are Moore et al. (2010) for HCN, Ugelow & Anderson (2022) for HC3N, Hudson & Yarnall (2022) for C6H6 and Khanna et al. (1988) for C4H2. Also shown are the imaginary refractive indices of Titan’s haze (Vinatier et al. 2012) and of Titan’s tholin analogues (Imanaka et al. 2012). The retrieved haze optical depth is unreliable shortwards of ~13 μm (thin dashed line).

In the text

	Fig. 5 Ethane emission over 11.7–12.8 μm. Black: observations. Red: model using a C2H6 profile equal to 1.05 times the Lavvas et al. (2021) profile and a telluric C2H6/13C12CH6 ratio. The inset shows a zoom on the 13C12CH6 lines, and includes models in which the C2H6/13C12CH6 ratio is multiplied (blue curve) or divided (green curve) by 2.
In the text

	Fig. 6 Methylacetylene (CH3C2H) and diacetylene (C4H2) emission over 15.4–16.2 μm. Black: observations. Red: model using CH3C2H and C4H2 profiles following the Lavvas et al. (2021) profiles but divided by factors of 9.5 and 4.5, respectively.
In the text

Fig. 7 Acetylene (C2H2) and its isotope emission over 12.8–14.8 μm. Black: observations. Colored curves: models. All models assume a terrestrial 12C/13C ratio. Top: LTE models (i.e,. using the Tkin profile) for various C2H2 profiles. Light blue: reference profile from Lavvas et al. (2021); red: vertically uniform C2H2 with 4 × 10−6 mixing ratio; green: slope profile from Fig. 8. Residuals for this latter model are shown. Bottom: Non-LTE models (adjusted Tvib) using the C2H2 profile from Lavvas et al. (2021). Green curve: nominal solution, using the dashed-dotted light blue Tvib(z) profile in Fig. 3. Observation minus model residuals are shown for this model, as well as for models using Tvib(z) + 1.5 K (dark blue) and Tvib(z) – 1.5 K (red). Insets: Zoom on the region of the ν5 band of C2HD at 14.75 μm. Models use either the C2H2 slope profile (top) or the Tvib solution profile (bottom), and D/H ratios in C2H2 equal to one, three, and five times the telluric value.

In the text

	Fig. 8 Volume mixing ratio profiles from Lavvas et al. (2021), compared to inferences from this work. Arrows indicate upper limits.
In the text

	Fig. A.1 Components of the fluorescence in the 7.2–8.2 μm range: ν4 band, shown separately for 12CH4 and 13CH4, 2ν4–ν4, and 3ν4–2ν4 bands. The CH3D contribution (ν3 band), calculated with the same model as for the ν6 band, has been divided by 2 (see text).
In the text

Fig. A.2 Model sensitivity for the CH4 fluorescence in some selected spectral intervals. Top panel: sensitivity to Tvib (ν4) profile (non-isothermal above 25 km vs isothermal at 125.5 K). Other panels: sensitivity to model omitting individual fluorescence contributions one-by-one: 13CH4, 2ν4–ν4, 3ν4–2ν4, CH3D ν3. In each panel, the nominal model including all contributions is shown in blue, and the alternative model in red.

In the text

	Fig. B.1 Contribution functions for the ν4 non-LTE emission of CH4. They are calculated for the Tvib(z) solution profile (solid green line in Fig. 3), a LOS intersecting the Pluto disk at r = 0.5 RPluto, and convolved to the instrumental resolution.
In the text

	Fig. B.2 Contribution functions for the ν9 band of C2H6. They are calculated for the solution C2H6 abundance and a LOS intersecting the Pluto disk at r = 0.5 RPluto, and convolved to the instrumental resolution.
In the text

	Fig. B.3 Contribution functions for CH3C2H and C4H2. They are calculated for the solution CH3C2H and C4H2 abundances, a LOS intersecting the Pluto disk at r = 0.5 RPluto, and convolved to the instrumental resolution.
In the text

	Fig. B.4 Contribution functions for the ν5 band of C2H2, in the “non-LTE” solution case. They are calculated for the Tvib(z) solution profile (dashed-dotted light blue line in Fig. 3), the nominal C2H2 profile from Lavvas et al. (2021), and a LOS intersecting the Pluto disk at r = 0.5 RPluto. and convolved to the instrumental resolution
In the text

	Fig. B.5 Contribution functions at selected wavenumbers for the different atmospheric species, calculated as described in Fig. B.1–B.4.
In the text

Fig. C.1 Tentative detection of HCN and upper limits on C2H4, C3H8, C6H6, CO2, and HC3N. For each species, the number in parentheses indicates the factor by which the nominal profile from Lavvas et al. (2021) has been multiplied. In all panels, model excluding (resp. including) the considered species is shown in blue (resp. in red). In the first panel (HCN), observation minus model residuals are also shown.

In the text