© The Authors 2025

1 Introduction

The planetary-mass companion GJ 504b is one of the few imaged planets to bridge the gap between the population of directly imaged young warm exoplanets ∼1000 K and our Solar System’s Jupiter at ∼130 K. Together with recent imaging of the planet Eps Ind b at ∼275 K (Matthews et al. 2024), GJ 504 b is one of the coldest planetary-mass companions (∼500 K) imaged to date. It orbits a solar-type star (spectral type G0V, Anderson et al. 2010) at a separation of 43 au, showcasing an orbit that is slightly beyond that of the nearby Neptune (i.e., 30 au). GJ 504 b exhibits bluer colors in the near-infrared (NIR) than any previously directly imaged exoplanet (J − H = −0.23 mag, Kuzuhara et al. 2013), but is redder than any observed brown dwarfs with a similar temperature or brightness, probing an unexplored parameter space of the color-magnitude diagram (Bonnefoy et al. 2018). This is shown in Fig. 1 with the color-magnitude diagram at NIR wavelengths. Its properties are rather typical of late T- type brown dwarfs, suggesting a largely cloud-free atmosphere. Indeed, Janson et al. (2013) confirmed the first detection of CH₄ in its atmosphere, as expected for a T-type object, which is hints at the presence of disequilibrium chemistry in its atmosphere.

Furthermore, the existence of this object represents a challenge for theories of planetary formation and evolution. To date, neither of the two most well-established planet formation scenarios can explain GJ 504 b’s current orbital separation and super-solar metallicity, raising fundamental questions about its origin (Skemer et al. 2016; Bonnefoy et al. 2018). In the core accretion model, planets form close to their stars, with a predicted typical outer boundary of ∼30 au. Forming planets slowly accrete dust, grains, and pebbles over a few millions of years, after which they accrete their gaseous envelope (Bodenheimer 2000). This scenario fails to explain the large orbital separation of GJ 504 b. The disk instability model mirrors the process of star and brown dwarf formation: planets are believed to form through the collapse of a massive protoplanetary disk due to gravitational instabilities, leading to fragmentation and subsequent planet formation (Marley et al. 2007). According to this model, planets are presumed to retain the chemical composition of the disk material at the location of their formation, meaning that the planet is expected to have the same metallicity as the host-star (Öberg et al. 2011; Mollière et al. 2022). The supersolar metallicity observed in GJ 504 b, along with its derived low mass of 4 M_Jup (Kuzuhara et al. 2013) cannot be explained with the disk instability model, unless it undergoes significant migration either inward or outward (Shibata et al. 2020; Turrini et al. 2021). The metallicity of the companion is more enriched ([M/H] ≈ +0.6) in metal than its parent star GJ 504 A ([M/H] ≈ +0.1 to 0.3, Skemer et al. 2016). Given that its atmosphere and the history of its formation and evolution are not well understood, previous studies have given rise to controversy over its planetary nature. Indeed, its derived mass relies on evolutionary models, based on the assumed age for this system.

The stellar gyrochronological and chromospheric activities of GJ 504 A have led to the determination of its characteristics as a rather young star (${160}_{-60}^{+350}$ Myr, Kuzuhara et al. 2013). Later studies indicated an older stellar age, between 1.5 and 4 Gyr: Fuhrmann & Chini (2015) and D’Orazi et al. (2017) suggested that the high levels of rotation and chromospheric activity, normally characteristic of a young stellar age, are due to the recent engulfment of a short-period hot Jupiter, further arguing in favor of an older system. Isochronal studies also show two solutions corresponding to the star: one being at a young age and the other at an older age (Bonnefoy et al. 2018). Finally, more recent studies favor the hypothesis of a young system, based on high- resolution spectra of the star and indicators of stellar activity (Šubjak et al. 2023; Di Mauro et al. 2022). However, they did not rule out the hypothesis of an older system. The degeneracies in the system’s age have led to highly disparate mass values for GJ 504 b; namely: $M={1.3}_{-0.3}^{+0.6}$ M_Jup or $M={23.3}_{-9}^{+10}$ M_Jup for the ages of young and old isochronal systems of 21±2 Myr and 4.0±1.8 Gyr, respectively (Bonnefoy et al. 2018). These values allow for such an object to be either one of the lowest-mass young planetary imaged companions known to date or, simply, an old brown dwarf.

In this work, we present the analysis of new observations from the James Webb Space Telescope (JWST) MidInfraRed instrument (MIRI) obtained within the framework of the ExoMIRI Guaranteed Time Observations (GTO) program 1277 (PI: P.-O. Lagage). They are also part of MIRIco, a EU/US coordinated observing effort with the MIRI coronagraphs between programs 1194, 1277, and 1241. In Sect. 2, we present the observational parameters, data reduction, and use of a reference star library. In Sect. 3, we describe the atmospheric characterization of the object. We discuss our results in Sect. 4 and our conclusions in Sect. 5.

Fig. 1 Color-magnitude diagram of the field brown dwarfs (M, L, and T dwarfs are shown in color), objects with confirmed youth or low gravity, and well-known directly imaged exoplanets (photometry from the UltracoolSheet, Best et al. 2024). GJ 504 b is at a unique position between the sequence of young low-gravity objects and old T-type brown dwarfs.

2 Observations and data reduction

2.1 Program observations

The system was observed with all three of MIRI’s 4-Quadrant Phase-Masks coronagraphs (4QPM, Rouan et al. 2000), along with the paired filters F1065C, F1140C, and F1550C. These filters are centered at 10.575, 11.30, and 15.50 µm, with a ∼5% bandwidth, proving the first mid-infrared (MIR) images of the system. Background observations are included for each filter to mitigate the “glowstick” effect identified during commissioning (Boccaletti et al. 2022). The background is observed using two dithers, which are then averaged to optimize its subtraction. No dedicated reference stars were observed during this sequence of observations. In fact, the purpose of this GTO sequence was precisely to test to which level of contrast the diffraction pattern can be subtracted out with one or several other reference stars from our program or from other programs. With a favorable separation of ∼2.5″, we were expecting the companion to be detected even if starlight subtraction would be degraded. As a result, we anticipate that stellar subtraction may not be optimal; however, this does not hinder the ability to detect the companion. The observation parameters are summarized in Table 1. The system is observed in the F1065C filter with twice more integrations (720 s instead of 360 s for the F1140 and F1550C filters), as GJ 504 b is expected to appear fainter at this specific wavelength, due to ammonia absorption.

2.2 Data reduction

The uncalibrated data were retrieved (_uncal files) from the Mikulski Archive for Space Telescopes MAST (Marston et al. 2018). The data reduction in this work was carried out in a similar way to that of previous datasets, as described in Boccaletti et al. (2024) and Mâlin et al. (2024). We ran stage 1 of JWST pipeline¹ (Bushouse et al. 2022), which applies essential detector-level corrections to all exposure types to obtain a corrected count-rate image (_rates files). Then, we ran stage 2 to subtract the background contribution and apply photometric calibration (_cal files). The flat-field correction is skipped to avoid increasing noise and the glowstick effect, as attenuation is more important at the edge of each quadrant. We checked that the impact on the photometry was lower than 2%, so much smaller than the other sources of noise (Boccaletti et al. 2024). Finally, we applied a σ-clipping function to correct the remaining bad pixels and NaN values. Any pixel with a value greater than 3σ compared to the median of its closest neighbors was replaced with this median value. The images were rotated to align north with the top of the image using the position angle (angle of ∼121.7° for the F1065C filter and 118.5° for the F1140C and F1550C filters). The MIRI coronagraphic images are fully dominated by diffraction, which we aim to subtract, as described in the following sections (see Sect. 2.4).

Fig. 2 Coronagraphic images of GJ 504 in each filter (F1065, F1140, and F1550) after the subtraction of an optimized-reference star. The coronagraphic center is illustrated with the small star in green and GJ 504 b is the bright source in the upper right corner, as indicated by the arrow in the F1065C image (first panel). An asinh color scale is used to show more details. Images in linear scale are available in the appendix.

2.3 Building a library of references

We built a library of reference star observations to estimate and remove the stellar diffraction pattern. We used all the available observations, at the time of this work (all reference star observations are publicly available in MAST), obtained with MIRI coronagraphic mode during the commissioning program 1037, ERS program 1386, and the GTO and GO programs (i.e., GTO programs 1277, 1194, 1413, 1411, and 1241, along with GO Cycle 1 programs 1668, 2153, 2243, and 2538). We also added reference observations from the GO cycle 2 available until June 1, 2024 (programs from 3254 and 3662). This information is summarized in the appendix. The result is a total of 10 stars observed in the F1065C filter, 18 stars in the F1140C, and 10 stars in the F1550C. We note that filter F1140C is the most requested one in MIRI observation programs. All reference star observations were reduced in the exact same way as described in Sect. 2.2. Each reference star has a spectral type that is similar to that of the target of each of the above observing programs, and it does not necessarily correspond to the spectral type of GJ 504 A. All the reference stars observations were carried out using the small-grid-dither (SGD) strategy, which includes a small offset (by steps of 10 mas) between each of the dithered observations to account for the fact that the coronagraph center is not perfectly measured (Lajoie et al. 2016). Most reference stars are observed with nine dithers, but some observational data sets use only five dithers. Taking these different dither positions into account, there are a total of 70, 146, and 66 observations in the F1065C, F1140C, and F1550C filters, respectively.

2.4 Stellar diffraction subtraction

Using the reference star library (described in Sect. 2.3), we used a principal component analysis (PCA) to reconstruct an optimized reference image. We choose the number of PCA components to remove to optimize the signal-to-noise ratio (S/N) for GJ 504 b. We noticed that the stellar residuals can be quadrantdependent, as illustrated in appendix, whereby one quadrant has a larger amount of flux compared to the other three. Therefore, we independently apply the PCA method to each of the four quadrants (method referred to as 4Q-PCA) to achieve the most effective subtraction. This final result is shown in Fig. 2 using the complete library of reference. The contrast limits curves for each filter are shown in Fig. 3 (plain line). They represent the achievable contrast sensitivity at 5σ as a function of the separation for each filter. We also show the contrast limits obtained with a single reference image for comparison (Fig. 3, dashed line). Using a larger reference library provides a gain of at least one order of magnitude in terms of contrast limits at all separations with the F1065C and F1140C. The gain in contrast is even larger at shorter separations.

Even if some stellar diffraction residuals persist at a shorter distance from the coronagraph center, the contrast sensitivity achieved at 5σ is limited by the background noise farther than 1″ , consistent with the MIRI commissioning results (Boccaletti et al. 2022). The profiles in Fig. 3 still show perceptible “bumps,” due to stellar diffraction residuals. The F1065C filter provides poorer contrast than the F1140C at separation closer to 1.5″, which may be due to the fact that the library contains fewer references in the F1065C filter. The gain at F1550C in using a reference library is lower than for the two other filters because the background noise is higher, so the diffraction residuals are not a hard limit even with fewer frames than for F1140C. Moreover, the library features heterogeneous exposures times, hence the shortest exposures drive the overall background noise. In addition, we also tested different subsets of the reference library and various algorithms to suppress the diffraction pattern. These include the linear optimization of the reference observations, in the same spirit as the algorithm Locally Optimized Combination of Images (LOCI, Lafreniere et al. 2007), but with a single optimization zone: a ring centered at the coronagraph center from 0.3″ to 8″. The companion was detected, regardless of the reference observations used to subtract stellar diffraction, as shown in the appendices. The S/N values for planet detection are indicated in the appendix, with the noise measured as the standard deviation of the flux in a ring of 4 λ/D (~8 pixels) at the separation of the planet. The classical PCA provided the highest S/N, but this method has the strongest residuals in the center of the images. The corresponding contrast curves measured with the different subtraction methods are presented in the appendix. The method described above using the entire library and the 4Q-PCA method provides the largest contrast.

Fig. 3 Contrast curve at 5σ computed for each filter image. All curves are computed using the 4Q-PCA stellar subtraction method. The plain lines correspond to the subtraction with the entire library and the dashed lines correspond to the case of using only the reference star from the GTO 1413 program. The shaded regions correspond to the inner working angle of the coronagraphs.

3 Atmospheric characterization

3.1 Extraction of the photometry

The photometry was extracted as described in Boccaletti et al. (2024) and Mâlin et al. (2024). We use both MJy/sr calibrated images _cal and _rates images, to which we applied a photometric calibration based on the contrast with the host star. The website whereistheplanet (Wang et al. 2021) was used to provide an estimate of the position at the observation date, based on past observations. We use WebbPSF (Perrin et al. 2014) to simulate the planet PSF models, taking into account the appropriate filter and mask configurations for MIRI coronagraphs. The planet’s position relative to the 4QPM axis was taken into account by specifying its position in the detector coordinates. The residuals between the PSF model and the data are minimized using the Nelder-Mead algorithm (Nelder & Mead 1965). The position of the PSF model is also optimized with two free parameters, which shift the model to x- and y-positions. In the case of PCA-based stellar subtraction methods, the planet’s photometry may be biased. Consequently, we modeled the planet’s PSF before applying the PCA algorithm and we then minimized the residuals. This reduces the uncertainties in planet photometry measured on images obtained with the various stellar subtraction methods. The best-model PSF is presented in the appendix. The flux is extracted on the best-fit PSF model image for each filter. The attenuation due to the coronagraph mask is evaluated with WebbPSF simulations: we measured the ratio between two simulated PSFs, the first at the planet’s position (within the detector frame) and the second at a position unaffected by coronagraph attenuation. The transmission at GJ 504 b’s position is 0.84 in F1065C, 0.93 in F1140C, and 0.77 in F1550C. Due to different position angles of the telescope, the companion is slightly closer to the transition of the 4QPM axis for the observation at F1065C, justifying the fact that the attenuation is more important than at F1140C. Each flux value is divided by this attenuation factor in order to recover the emitted flux of the object. There are two main sources of uncertainties on the flux extracted for the companion. First, the stellar subtraction can have an impact; thus, we repeated the same procedure for each method of stellar subtraction (with the entire library, as shown in the appendix) We obtain consistent values with a variation of 0.8, 1.3, and 6.9% in the F1065C, F1140C, and F1550C filters, respectively. This uncertainty is labeled σ_stellarsub in the following. Secondly, the PSF normalization from DN to (W/m²/µm) is another source of uncertainty, which has not yet been documented. Therefore, we compared these photometric values with those derived with a contrast approach. First, we used the method from Boccaletti et al. (2024), based on a contrast measurement and estimation of the stellar flux using target acquisition observations. We added the method based on simulation with WebbPSF to estimate the stellar flux as in Mâlin et al. (2024). The different results obtained for the flux measurement are presented in the appendix. All methods provided consistent values, but the PSF normalization to obtain physical flux units remains the main source of uncertainty, with a variation of 4.1, 12.6, and 10.7% for each of the three filters, respectively. We refer to this uncertainty value as σ_{PS Fnorm}. Finally, we measured the standard deviation of the residual stellar flux in a ring at the planet separation, to check whether the remaining stellar diffraction may be an additional source of uncertainty in the flux measurement. These standard deviations are on the order of 0.1% of the flux of the planet in each filter and can consequently be neglected. The final photometric flux values are indicated in Table 2. They correspond to the averaged flux on the three normalization methods, measured using images from Fig. 2, which provide the highest S/N values for the detection of GJ 504 b. We took the final uncertainty as the quadratic sum of the two main sources of uncertainties (Eq. (1)). $${\sigma }_{flux}=\sqrt{{\sigma }_{stellarsub}^2+{\sigma }_{PSFnorm}^2}.$$(1)

The dispersion between the flux values obtained by the different methods is larger in the F1140C filter than in the F1065C, even though the detection occurs with a higher S/N at this wavelength.

Based on the best fit PSF model, we determined the astrometry relative to the center of the coronagraphic mask to be (∆RA, ∆DEC) = (−1.75 ± 0.02″, 1.77 ± 0.01″), averaged over the three filters. The uncertainties correspond to the dispersion of the measurement over the three filters. This is consistent with previous predictions within the uncertainties and corresponds to a separation of 2.48± 0.02″ (i.e., 44.4 ± 0.7 au at 17.56 pc).

3.2 Atmospheric characterization of the planet

We used the self-consistent 1D atmospheric model Exo-REM, developed to simulate the atmosphere of young giant planets (Baudino et al. 2015) and to understand the L–T transition (Charnay et al. 2018). The grid parameter ranges are described in Table 3. Exo-REM includes a cloud model that takes into account the microphysics (iron and silicate clouds with supersaturation parameter S=0.003). Disequilibrium chemistry is also included in the model. The sources of opacity include collision-induced absorption of H₂–H₂ , H₂–He, H₂ O–H₂O, and H₂O–N₂ , the ro- vibrational bands of molecules (H₂O, CH₄, CO, CO₂, NH₃, PH₃, TiO, VO, H₂S, HCN, and FeH). Line lists are given in Blain et al. (2021). The NIR photometric values were taken from the literature, coming from ground-based instruments: Subaru/CIAO (Kuzuhara et al. 2013; Janson et al. 2013), LBTI/LMIRcam (Skemer et al. 2016), and VLT/SPHERE (Bonnefoy et al. 2018). The points in Fig. 4. represent the photometry values colored- coded by instrument. No additional scaling factors are used in between instruments.

We used a forward-modeling process, using the python package species² (Stolker 2023) to measure the atmospheric parameters and their posterior distributions. We ran the Bayesian analysis with 5000 live points, without prior constraints on the atmospheric parameters. We ran the analysis twice, once with only the NIR photometry and then again, adding the MIRI photometric points. The best-fit parameters obtained are summarized in Table 4. The best-fit model using NIR and MIR photometry is shown in Fig. 4. The posterior distributions for each parameter are presented in the appendix.

The atmospheric parameters measured by adding MIRI photometric points are consistent to the 1σ level with the atmospheric parameters previously estimated, based solely on near-infrared data. The mass value is derived from the surface gravity and radius values; and the luminosity is measured from the temperature and radius. Expanding the wavelength range reduces uncertainties in the radius measurement by at least a factor 3, consequently enhancing the precision of the logarithmic measure of luminosity by a factor ~3.5. These parameters are evaluated independently of the assumptions about the age of the system and the evolution models. The mass values seem to be more consistent with the young age hypothesis; however, it relies on the surface gravity measurement, which is not confidently reliable with only photometric points. Indeed, the surface gravity is embedded in the shape of the lines rather than on the continuum.

NH₃ absorption explains the lower flux at F1065C than at F1140C. We generated a second Exo-REM atmospheric grid based on the best-fit parameters and varying the volume mixing ratio (vmr) of NH₃ from 5⋅10⁻¹⁰ to 1⋅10⁻⁴. Following Danielski et al. (2018), we measured the detection level of NH₃ by comparing it to the model with the lowest NH₃ abundance, set at 5⋅10⁻¹⁰, which has a negligible effect on the spectrum. It is expressed as $$S/N_{\left(N{H}_3\right)}=\frac{F_{noN{H}_3,F1065C}-F_{obs,F1065C}}{{\sigma }_{tot}},$$(2)

where σ_tot is the quadratic sum of the relative uncertainties on the measured flux, $${\sigma }_{tot}=\sqrt{{\sigma }_{F1065C}^2+{\sigma }_{F1140C}^2}.$$(3)

Comparing the flux difference between F1065C and F1140C to the photometric errors (derived in Sect. 3.1) yields a result of $S/N_{\left(N{H}_3\right)}=3\sigma$. However, these uncertainties correspond to an absolute photometric precision; hence, the outcome for NH₃ is quite conservative and must be taken as a lower limit. Instead, given it is a relative measurement, a more reliable estimate of the NH₃ detection level should be derived from each photometric method independently. This means zeroing the term ${\sigma }_{PSFnorm}^2$ in Eq. (1). Taking this precaution provides a more reliable estimate at 12.5 σ. Finally, we calculate the χ² values for each model in this second atmospheric grid. Models with a vmr of NH₃ in the range –5.37 and –5.22 (in log₁₀) fall within the 1 σ confidence level of the χ² minimisation (see the figure in the appendix).

4 Discussion

4.1 Stellar subtraction

To remove the stellar diffraction, we reconstructed an optimized reference image subtracted from the data using traditional algorithms previously developed for observations with ground-based instruments (such as PCA and linear combination of references images). These algorithms are adapted for MIRI coronagraph data. Moreover, we used individual sets of references from a specific program and then the entire reference library. We found that using a single set of reference is not effective enough to subtract the stellar diffraction pattern from coronagraphic images, at least in the particular case of GJ 504 observations. We note that the bright residuals that could be assimilated as point sources are highly dependent on the reference dataset and the stellar diffraction subtraction methods used. Hence, they cannot be considered as real new point sources.

The reference star from the GTO 1413 program is of a spectral type and magnitude similar to GJ 504 A and it was observed a few days after the GJ 504 system. It represented an optimistic dataset to provide a good stellar subtraction. However, a small wavefront drift occurred between these two observations and the wavefront error increased by ~6 nm (measures available with WebbPSF), which could have had an impact on the quality of this data set. The quality of stellar subtraction is improved when we increase the number of datasets in the reference library, thereby achieving more favorable contrasts. Previous works conducted within this GTO 1277 program (systems HR 8799 and HD 95086; Boccaletti et al. 2024; Mâlin et al. 2024) indicated that the use of a reference library did not yield an improved stellar subtraction. For these observations, dedicated reference stars were observed together with the scientific target. Furthermore, the library used was built with fewer references available, coming only from the ERS, GTO and commissioning programs. We tested our updated library, containing all GTO and GO programs available (cycle 1 and part of cycle 2) for these two datasets, but it does not provide improved results. We conclude that obtaining reference observations captured in the same sequence as the target observations remains the optimal method for mitigating stellar diffraction in coronagraphic images, especially for closer-in planets (such as the planets from the systems HD 95086 and HR 8799). However, the accumulation of a larger reference library in the next few years could provide better results in the future, allowing for the capture and removal of diffraction residuals at separations shorter than 1″ .

Furthermore, the PCA analysis provides better subtraction (such as for HIP 65426 system, Carter et al. 2023) for GJ 504, in contrast to the HR 8799 and HD 95086 systems for which the optimized linear combination of reference yields a better subtraction. This is likely due to the fact that the systems GJ 504 and HIP 65426 do not contain a warm inner disk (or multiple bright planets and prominent background objects), which prevents a straightforward PCA-based stellar subtraction. In conclusion, the best algorithm for subtracting the stellar diffraction in MIRI data strongly depends on the architecture of the system itself.

The asymmetries in the stellar residuals visible in the coronagraphic image due to the 4QPM led us to apply the PCA in each quadrant independently. Future algorithms adapted specifically to the 4QPM coronagraph could improve the stellar subtraction. In the case of the GJ 504 system, using the full library provides all the diversity available. We note that we did not get any improvement in the performance when we selected samples from the library. For more challenging observations (fainter targets or closer separations), we could achieve a better performance by evaluating the choice of references in the library, avoiding those containing an inner disk, and optimizing the regions of the field of view where the residuals should be minimized.

Fig. 4 Exo-REM best-fit model (black curve) and the lighter black models come from the posterior distribution from the best fit. The points represent the photometric values from VLT/SPHERE (green), Suburu/HiCIAO (purple), LBTI/LMIRcam (orange), and Subaru/ICRS (blue). The MIRI photometry is shown in red. Photometric points are plotted at the filter’s central wavelength, with the width indicated. The bottom plot represent the residuals between the photometric points and the best fits models: all point agree within a 3 σ precision. Right: zoom on the MIRI data and the NH3 absorption, fluxes values are displayed with a linear scale.

4.2 Spectral characterization

We present the detection of NH₃ in the atmosphere of a directly imaged planetary-mass companion at 12.5 σ using the relative measurement of the flux between the two filters F1065C centered on the absorption of ammonia and F1140C in the continuum of the spectrum. Using a conservative measurement of the flux and uncertainties that takes several methods into account, we estimated a lower limit of detection of 3σ. Even though it has been detected in many T-type isolated brown dwarfs (for example, Suárez & Metchev 2022), this molecule has been previously inferred in the atmosphere of only one directly imaged planet (detection at 2.7 σ for 51 Erib using retrieval analysis with its NIR data by Whiteford et al. 2023). Using pre-computed ExoREM grids and varying the volume mixing ratio of NH₃ , we measured a vmr of 10^−5.3±0.07 for GJ 504 b. For comparison, this is at least an order of magnitude lower than in Jupiter’s atmosphere (Taylor & Atreya 2004).

The NH₃ is sensitive to gravity (higher NH₃ abundance with higher surface gravity, as shown in Fig. 5) and unlike the CH₄/CO ratio, its ratio NH₃/N₂ is insensitive to mixing, making this molecule an interesting proxy for gravity (Zahnle & Marley 2014). In the case of a planetary mass object for which T_eff and the radius can be measured more accurately, the measurement of the abundance of NH₃ could therefore provide an additional estimate of the mass, since mass is directly linked to radius and surface gravity. The measured value of NH₃ abundance for GJ 504 b is more consistent with an object of a surface gravity log𝑔 < 4, according to Exo-REM models (Fig. 5). As a result, GJ 504 b photometry appears more consistent with a planetary-mass object rather than a brown dwarf. The presence of clouds in the atmosphere may affect this measurement, albeit to a lesser extent than metallicity. The detection of NH₃ is also promising to allow for the measurement of the isotopic ratio ¹⁴N/¹⁵N, as done for isolated Y-type brown dwarfs with MIRI/MRS spectra (Barrado et al. 2023; Kühnle et al. 2024). Young stars and consequently their planets should be more strongly enriched in the ¹⁵N isotope (Adande & Ziurys 2012). Determining the ¹⁴N/¹⁵N isotopic ratio, the C/N ratio (Turrini et al. 2021; Pacetti et al. 2022), and C/O ratio (e.g., Öberg et al. 2011; Madhusudhan 2019) can provide important constraints on the formation location and pathway of GJ 504 b.

The atmospheric parameters measured (temperature, radius, and luminosity) are compared to ATMO evolutionary models that include disequilibrium chemistry (Phillips et al. 2020) in Fig. 6. A smaller radius is measured when the MIR photometric points are added, together with lower uncertainties (red points, in comparison to the black points, Fig. 6). The comparison of luminosity and measured radius corresponds to isochrones between 400 Myr and 1 Gyr (at 1 σ). Compared with the Sonora evolution models (Marley et al. 2021), we find that the radius and effective temperature measured from the atmospheric fit correspond to isochrones slightly older, from ~500 Myr to 1.5 Gyr. The same trend is observed when the measured values are compared with the effective temperature expected from evolution models, rather than with the radius value. The MIR photometry seems to place GJ 504 b within the older age range. However, this is insufficient to confirm the nature of GJ 504 b, as this range of stellar ages translates into a wide range of masses. Indeed, for the isochrones between 400 Myr to 1 Gyr, the masses from ~ 1 to 17 M_Jup are consistent with the measured radius.

The mass measured from the surface gravity of the best- fit atmospheric model corresponds to the low-mass hypothesis. One can argue that the log𝑔 value is not reliable when fitting only a few photometric points, as this parameter is measured with higher uncertainties. At MIR wavelengths, surface gravity influences the shape of the spectral lines rather than the spectral continuum. Finally, the radius value might still be inconsistent with evolutionary models. Even with the addition of MIR information, atmospheric measurements still provide a slightly smaller radius than that predicted by the evolution models. Therefore, confirming the nature of GJ 504 b based on radius measurement may not be reliable. A deeper atmospheric analysis including different atmospheric models is outside the scope of this study, but should be carried out in the future with the coming JWST observation of GJ 504 b providing high S /N spectra at both NIR and MIR wavelengths.

Fig. 5 Abundance of the NH3 according to Exo-REM models as a function of the Teƒƒ for several values of log𝑔. The dashed lines represent the expectation for an atmosphere with a super-solar metallicity, and the plain line with a solar metallicity.

Fig. 6 Isochrones from 20 Myr to 4 Gyr: evolution of the luminosity of a planet (log scale and relative to the Sun’s luminosity) as a function of the radius (top panel) or the temperature (bottom panel) of a planet according to ATMO evolutionary models.

5 Conclusion and perspective

In this paper, we present our analysis of the MIRI coronagraphic images of one of the coldest directly imaged planetary-mass companion to date. Our findings are as follows:

• The stellar diffraction was subtracted from the coronagraph image using a reference library built on previous JWST observations. For the first time, we used a large library with all the observations available up to June 2024. This provides better performance than the previous library (used in Boccaletti et al. 2024; Mâlin et al. 2024), which only contained ERS and early GTO observations Even though this method does not perfectly remove the stellar residuals at shorter separation, it holds promise for achieving a better contrast when no reference observations are available;

• We confidently detected the presence of NH₃ in the atmosphere of a late T-type planetary mass companion;

• The atmospheric parameters were constrained with lower uncertainties, owing to photometric values measured at MIR wavelengths. The luminosity and radius measurement are more accurate by at least a factor of 3;

• The addition of MIRI photometry enabled the measurement of a smaller radius of R = 1.08 R_Jup. This value is in agreement with the isochrones from 400 Myr to 1 Gyr, according to evolutionary models;

• Although the MIR photometry is still not sufficient to narrow down the mass and confirm the planetary nature of GJ 504 b; these MIRI observations show the first MIR data and are valuable for forthcoming spectroscopic observations of this object.

The GJ 504 system has recently been observed with both integral-field spectrographs of JWST (NIRSpec/IFU, GTO 2778, PI: M. Perrin and MIRI/MRS, GO 3647, PI: P. Patapis). Furthermore, recent observations from VLT/CRIRES (PI: F. Kiefer) will provide high-resolution data at NIR wavelengths. The MIRI photometry presented here offers new constraints on this T-dwarf planetary mass companion, which will prove useful in analyses of the forthcoming datasets. The wealth of spectroscopic data coming will allow for an improved characterization of its atmosphere, while confirming its nature and to allowing us to better understand its formation and evolution pathway.

Data availability

The appendices are available on Zenodo: zenodo.14517224.

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 #1277. Part of this work was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under contrast with NASA (80NM0018D0004). This publication makes use of VOSA, developed under the Spanish Virtual Observatory (https://svo.cab.inta-csic.es) project funded by MCIN/AEI/10.13039/501100011033/ through grant PID2020- 112949GB-I00. VOSA has been partially updated by using funding from the European Union’s Horizon 2020 Research and Innovation Programme, under Grant Agreement nº 776403 (EXOPLANETS-A) This work has benefitted from The UltracoolSheet at http://bit.ly/UltracoolSheet, maintained by Will Best, Trent Dupuy, Michael Liu, Aniket Sanghi, Rob Siverd, and Zhoujian Zhang, and developed from compilations by Dupuy & Liu (2012), Dupuy & Kraus (2013), Deacon et al. (2014), Liu et al. (2016), Best et al. (2018), Best et al. (2021), Sanghi et al. (2023), and Schneider et al. (2023). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France (DOI : 10.26093/cds/vizier). The original description of the VizieR service was published in 2000, A&AS 143, 23 Software: python, astropy, numpy, scipy matplotlib. M.M., A.B., P.-O.L., C.C. acknowledges funding support by CNES. B.V. thanks the European Space Agency (ESA) and the Belgian Federal Science Policy Office (BELSPO) for their support in the framework of the PRODEX Programme. O.A. is a Senior Research Associate of the Fonds de la Recherche Scientifique – FNRS. OA thanks the European Space Agency (ESA) and the Belgian Federal Science Policy Office (BELSPO) for their support in the framework of the PRODEX Programme. D.B. is supported by Spanish MCIN/AEI/10.13039/501100011033 grant PID2019-107061GB-C61 and No. MDM-2017-0737. L.D. acknowledges funding from the KU Leuven Interdisciplinary Grant (IDN/19/028), the European Union H2020-MSCA-ITN- 2019 under Grant no. 860470 (CHAMELEON) and the FWO research grant G086217N. MPIA acknowledges support from the Federal Ministry of Economy (BMWi) through the German Space Agency (DLR). J.P. acknowledges financial support from the UK Science and Technology Facilities Council, and the UK Space Agency. G.O. acknowledge support from the Swedish National Space Board and the Knut and Alice Wallenberg Foundation. P.P. thanks the Swiss National Science Foundation (SNSF) for financial support under grant number 200020_200399. P.R. thanks the European Space Agency (ESA) and the Belgian Federal Science Policy Office (BELSPO) for their support in the framework of the PRODEX Programme. E.C. acknowledges funding by the European Union (ERC, ESCAPE, project No 101044152). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. C.D. acknowledges support from the INAF initiative “IAF Astronomy Fellowships in Italy”, grant name GExoLife. E.v.D. acknowledges support from A-ERC grant 101019751 MOLDISK. The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant No. 140. TRG is grateful for support from the Carlsberg Foundation via grant No. CF20-0534. T.P.R. acknowledges support from the ERC 743029 EASY. Support from SNSA is acknowledged.

References

All Tables

All Figures

	Fig. 1 Color-magnitude diagram of the field brown dwarfs (M, L, and T dwarfs are shown in color), objects with confirmed youth or low gravity, and well-known directly imaged exoplanets (photometry from the UltracoolSheet, Best et al. 2024). GJ 504 b is at a unique position between the sequence of young low-gravity objects and old T-type brown dwarfs.
In the text

Fig. 2 Coronagraphic images of GJ 504 in each filter (F1065, F1140, and F1550) after the subtraction of an optimized-reference star. The coronagraphic center is illustrated with the small star in green and GJ 504 b is the bright source in the upper right corner, as indicated by the arrow in the F1065C image (first panel). An asinh color scale is used to show more details. Images in linear scale are available in the appendix.

In the text

	Fig. 3 Contrast curve at 5σ computed for each filter image. All curves are computed using the 4Q-PCA stellar subtraction method. The plain lines correspond to the subtraction with the entire library and the dashed lines correspond to the case of using only the reference star from the GTO 1413 program. The shaded regions correspond to the inner working angle of the coronagraphs.
In the text

Fig. 4 Exo-REM best-fit model (black curve) and the lighter black models come from the posterior distribution from the best fit. The points represent the photometric values from VLT/SPHERE (green), Suburu/HiCIAO (purple), LBTI/LMIRcam (orange), and Subaru/ICRS (blue). The MIRI photometry is shown in red. Photometric points are plotted at the filter’s central wavelength, with the width indicated. The bottom plot represent the residuals between the photometric points and the best fits models: all point agree within a 3 σ precision. Right: zoom on the MIRI data and the NH3 absorption, fluxes values are displayed with a linear scale.

In the text

	Fig. 5 Abundance of the NH3 according to Exo-REM models as a function of the Teƒƒ for several values of log𝑔. The dashed lines represent the expectation for an atmosphere with a super-solar metallicity, and the plain line with a solar metallicity.
In the text

	Fig. 6 Isochrones from 20 Myr to 4 Gyr: evolution of the luminosity of a planet (log scale and relative to the Sun’s luminosity) as a function of the radius (top panel) or the temperature (bottom panel) of a planet according to ATMO evolutionary models.
In the text