The ESO SupJup Survey V: exploring atmospheric variability and orbit of the super-Jupiter AB Pictoris b with CRIRES+

ABSTRACT

A growing number of directly-imaged companions have been recently characterized, with robust constraints on carbon-to-oxygen ratios and even isotopic ratios. Many companions and isolated targets have also shown spectral variability. In this work, we observed the super-Jupiter AB Pictoris b across four consecutive nights using VLT/CRIRES+ as part of the ESO SupJup survey, exploring how the constraints on chemical composition and temperature profile change over time using spectral line shape variations between nights. We performed atmospheric retrievals of the high-resolution observations and found broadly consistent results across all four nights, but there were differences for some parameters. We clearly detect H|$_2$|O, |$^{12}$|CO, and |$^{13}$|CO in each night, but abundances varied by |${\sim} 2\sigma$|⁠, which was correlated to the deep atmosphere temperature profiles. We also found differences in the |$^{12}$|C|$/^{13}$|C ratios in each night by up to |${\sim} 3\sigma$|⁠, which seemed to be correlated with the cloud deck pressure. Our combined retrieval simultaneously analysing all nights together constrained broadly the average of each night individually, with the C/O|$=0.59\pm 0.01$|⁠, consistent with solar composition, and |$^{12}$|C|$/^{13}$|C |$= 102\pm 8$|⁠, slightly higher than the Interstellar Medium (ISM) and Solar System values. We also find a low projected rotational velocity, suggesting that AB Pictoris b is either intrinsically a slow rotator due to its young age or that the spin axis is observed pole-on with a |${\sim} 90^\circ$| misalignment with its orbit inclination. Future observations will be able to further explore the variability and orbit of AB Pictoris b as well as for other companions.

1 INTRODUCTION

Directly-imaged companions at wide orbital separations represent unique targets to characterize formation at the extremes of orbital separation. They are a new frontier in exoplanet astronomy and have been characterized in unprecedented detail thanks to developments in instrumentation and adaptive optics (see e.g. Currie et al. 2023). These companions have an internal source of heat from their formation and potentially deuterium burning, rather than strong incident stellar radiation as is the case for hot Jupiters. Therefore their temperature profiles are more steeply increasing with pressure, resulting in strong absorption features from spectroscopically active molecular and atomic species.

A range of chemical species have been detected in isolated objects as well as companions through retrievals, which constrain the atmospheric properties by performing Bayesian analyses and exploring a wide range of parameter space. The chemical constraints can be used to determine the carbon-to-oxygen (C/O) ratio and metallicity of these objects (e.g. Line et al. 2015; Burningham et al. 2017; Lavie et al. 2017; Line et al. 2017), key parameters in planet formation (Öberg et al. 2011; Madhusudhan 2012; Mordasini et al. 2016). Atmospheric retrievals have also been able to constrain cloud properties (e.g. Burningham et al. 2021; Vos et al. 2023) as well as chemical disequilibrium processes (Mollière et al. 2020; Regt et al. 2024) thanks to the often very high signal-to-noise observations. Such retrievals are now also opening up atmospheric surveys across the population of companions as well as comparative science with the compositions of close-in hot Jupiters (e.g. Welbanks et al. 2019; Xuan et al. 2024a).

A new paradigm shift in recent years has been to explore isotopic ratios (Mollière & Snellen 2019; Morley et al. 2019). The atmospheric |$^{12}$|C/|$^{13}$|C has been measured for close-in transiting planets, directly-imaged companions as well as field brown dwarfs (e.g. Line et al. 2021; Zhang et al. 2021a, b; Costes et al. 2024; Regt et al. 2024; Xuan et al. 2024c). These measurements have pointed to a range in values, with some objects nearer |$^{12}$|C/|$^{13}$|C|$\sim$|30, but others with |$^{12}$|C/|$^{13}$|C values |$\gtrsim$|100. Thanks to the arrival of JWST, we have been able to obtain high-precision constraints on the carbon isotope ratios (Gandhi et al. 2023b; Hood et al. 2024; Lew et al. 2024), and been able to explore spectroscopically weaker oxygen and nitrogen isotopes in planetary atmospheres for the first time (Barrado et al. 2023; Gandhi et al. 2023b; Xuan et al. 2024c). Recent work has also constrained the D/H ratio of a Y dwarf (Rowland et al. 2024), proving the exciting potential of isotopic ratios. The isotopic ratios are influenced by several fractionation processes in protoplanetary discs, such as isotope selective photodissociation, gas/ice partitioning, and isotopic exchange reactions, and are therefore a potential avenue to disentangle how companions form and migrate (e.g. Miotello et al. 2014; Zhang et al. 2021a; see Fig. 1). Recently, Landman et al. (2024) observed a directly-imaged companion, |$\beta$| Pictoris b, with the upgraded CRIRES+ (Cryogenic High-Resolution Infrared Echelle Spectrograph) on the VLT. Such facilities are ideal given their high spectral resolution allowing us to resolve many thousands of molecular lines in the spectrum, resulting in robust detections of a range of chemical species and isotopologues in the atmosphere.

Directly-imaged companions and isolated sub-stellar objects have shown significant variability over time (e.g. Artigau et al. 2009; Radigan et al. 2012, 2014; Biller et al. 2015, 2024; Metchev et al. 2015; Eriksson et al. 2019; Zhou et al. 2022; Liu et al. 2024). These sub-stellar mass objects are often fast rotators (e.g. Snellen et al. 2014; Bryan et al. 2020; Tannock et al. 2021; Hsu et al. 2024; Landman et al. 2024; Morris et al. 2024) and inhomogeneities which vary across the observed face as they rotate result in spectroscopic variability. Theoretical works modelling the L/T transition have proposed that clouds, thermochemical instabilities, thermal variations, and aurorae can all contribute to variability across the surface (Marley et al. 2012; Robinson & Marley 2014; Hallinan et al. 2015; Tremblin et al. 2016). Three-dimensional modelling has also shown that feedback between clouds and chemistry can drive inhomogeneities in the photosphere (e.g. Tan & Showman 2021; Lee et al. 2024). Younger, lower gravity objects have shown evidence for thick clouds (e.g. Barman et al. 2011; Currie et al. 2011), and even more variability compared with field brown dwarfs (Vos et al. 2019, 2022). Therefore, sub-stellar mass objects can often be dynamic and changing, but this gives us an excellent opportunity to explore how the physical and chemical properties vary over the planetary surface (e.g. Crossfield et al. 2014).

In this work, we characterize the atmosphere of the directly-imaged companion AB Pictoris b (hereafter AB Pic b) through observations made with the CRIRES+ spectrograph on the VLT as part of the SupJup survey (Program ID: 1110.C-4264, PI: Snellen). AB Pic b appears redder on a colour–magnitude diagram than other targets from the SupJup survey of similar spectral type (Regt et al. 2024), typical for young and lower gravity objects. First discovered by Chauvin et al. (2005), this L0-1 object has a temperature of |$\sim$|1700 K (Bonnefoy et al. 2010, 2014; Patience et al. 2012; Palma-Bifani et al. 2023) and orbits its 13.3 Myr K1V type host star (Booth et al. 2021) at a projected orbital separation of 270 au (5.378 arcsec on-sky separation), as shown in Fig. 1. Its mass is between the range for planets and brown dwarfs and thus presents an ideal target to characterize, particularly given that we have no such analogue in the Solar system at this mass or at such a wide orbital separation.

We observed AB Pic b over four consecutive nights, determining the spectral variability of the data over each of the nights and performing atmospheric retrievals to explore differences in the temperature profile, atmospheric composition, and cloud properties on the time-scale of |$\sim$|days. We use the variations in the spectral line shapes and depths over each of the nights to explore this variability, as high-resolution spectroscopy is most reliable when exploring the individual lines and their ratios rather than continuum fluxes (e.g. Pelletier et al. 2023; Maguire et al. 2024). Our observations also allow us to constrain the carbon isotopes given the strong CO absorption features in the K band. In addition to the chemical and physical properties of the atmosphere, the R|$\sim$|100 000 observations allow us to determine the rotational velocity and radial velocity of AB Pic b, providing an excellent opportunity to probe the orbital solution of the companion.

The next section discusses the methodology, describing the data analysis and atmospheric modelling, followed by the results and discussion for each night of observation individually as well as all four nights combined. Finally, we present our conclusions in Section 4.

2 METHODS

We observed AB Pic b with CRIRES+ (Kaeufl et al. 2004; Dorn et al. 2014, 2023) on the VLT’s UT3 as part of the ESO SupJup survey. We performed atmospheric retrievals using the HyDRA retrieval framework (Gandhi et al. 2019, 2023b). Below, we discuss the observations of AB Pic b, the reduction and telluric correction of the data and the modelling setup for the atmospheric retrievals.

2.1 Observations

AB Pic b was observed across four consecutive nights between 01-11-2022 to 04-11-2022, with each of the observations performed towards the end of the night when AB Pic system was at its lowest airmass (See Fig. 2 and Table 1). For all of the nights the seeing was relatively good (⁠|${\lesssim} 0.8$| arcsec) and the coherence times were high, particularly for nights two and three. This resulted in a clear signal of the companion, uncontaminated by starlight from AB Pic A given the |${\sim} 5.4$| arcsec separation, and with a good AO correction (see below). We observed AB Pic b into morning twilight for each night as well, but given the K-band observations this did not significantly impact the data quality. Our longest set of observations was for the first night, with 11 AB pairs, and our shortest was for night 2 with 6 AB pairs. We adopted a 300 s exposure (DIT), with NDIT set to 1 for every set of observations.

Fig. 3 shows the slit viewer position and orientation with respect to the system. We used the 0.4 arcsec slit for the observations, with a 6 arcsec nod throw, centred on the companion AB Pic b. This allowed us to maximize the companion’s light entering the slit. We adopted the SKY position for the derotator to ensure the position of the companion and the host star remained the same throughout the observations. We observed AB Pic b alone due to the large on-sky separation between its host star, which made placing both on the slit difficult. This configuration also allowed us to adopt longer exposure times which reduced the relative readout noise. We used the host AB Pic A, a 13.3 Myr old K1V star with H-band magnitude 7.09, as both the natural guide star for the adaptive optics and the slit viewer guide star given its brightness. The high coherence times and good seeing resulted in the actual spectral resolution derived for each night (shown in Table 1) being significantly higher than the quoted R|$\sim$|48 000 for K2166 observations with CRIRES+ . We calculated the spectral resolution of our observations by using the width of the trace on the detector, and assumed that the PSF was the same in all directions, as previously done by Holmberg & Madhusudhan (2022) and Nortmann et al. (2024). Each of the nights also had low precipitable water vapour, with relative humidity |${\lesssim} 5~{{\ \rm per\ cent}}$|⁠, resulting in low telluric H|$_2$|O contamination.

2.1.1 Data reduction

We reduced the data using the excalibuhr¹ pipeline (Zhang et al. 2024), a reduction pipeline for CRIRES+ observations of transiting and directly-imaged companions (see also Holmberg & Madhusudhan 2022). The first step of the pipeline is to create the master dark frames, flat, readout noise, and bad pixel map from the calibration frames. The traces of the spectral orders are then identified and the slit curvature is corrected using the Fabry–Pérot Etalon frames. The flux of the flat frame is then extracted and used to obtain the blaze function. We then subtract each AB nod pair to remove the sky background and dark current, and the normalised flat field is applied to each of the science frames. Each A and B nod position is then mean combined, and the tiled orders and slits in detector images are straightened. The 1D spectrum is extracted using an optimal extraction algorithm (Horne 1986) using an extraction aperture of 15 pixels. Note that the B position trace for the shortest wavelength observations is very close to the edge of the detector, but this was not significant for our observations as we only used the reddest 5 orders (see below).

2.1.2 Telluric standard

To correct for the telluric absorption we used standard stars, which were |$\mu$| Pictoris the first night, 53 Psc for the second, |$\lambda$| tau for the third, and |$\zeta$| Phe for the final night. These were used to refine the wavelength solution given that the actual observations of AB Pic b did not have sufficient signal-to-noise per pixel for high accuracy wavelength calibrations. We use a SkyCalc model (Noll et al. 2012; Jones et al. 2013) and applied a third order polynomial on the initial solution from the extraction of the standard star, and used this refined wavelength grid for our AB Pic b observations. Fig. 4 shows the fitted telluric spectrum with the SkyCalc model and the observations of |$\mu$| Pic taken on the first night. The spectral features show excellent agreement with the SkyCalc model, with the only difference occurring for the bluest order where there is heavy telluric contamination and the 2161–2171 nm range of the Brackett |$\gamma$| line, both of which were masked out for our analysis.

2.1.3 Extracted spectrum

To calibrate the flux we adopted the measured 2MASS and VLT/NACO brightness of K = 14.14 (Cutri et al. 2003; Chauvin et al. 2005) and GAIA distance of 50.14 pc (Gaia Collaboration 2020). We masked out the regions with heavy telluric contamination to avoid interference in the atmospheric retrievals. We also masked out the Brackett |$\gamma$| line in order 3. This masking was carried out independently for each night, as the seeing, precipitable water vapour and average airmass all varied from night to night.

2.2 Atmospheric modelling

We perform atmospheric retrieval analysis on the spectrum of AB Pic b to determine the chemical composition and temperature profile of the atmosphere for each of the four nights. This uses the the HyDRA retrieval framework (e.g. Gandhi et al. 2019, 2022). This has previously been adapted to retrieve the atmospheres of non-irradiated atmospheres (Gandhi et al. 2023b).

2.2.1 P–T profile

We use the method of Regt et al. (2024) to determine the P–T profile of the atmosphere, using a grid of P–T points and weighting the 3rd order derivative to reduce ringing/oscillatory behaviour. This is a likelihood penalty imposed for temperature profiles with high 3rd order derivatives (see Li & Cao 2022). This is given by

where |$\mathbf {D}_3\mathbf {C}$| is the 3rd order general difference matrix, which computes the P–T points on the knots and weights them by their separation in log pressure for our unevenly spaced P–T points (see Gandhi et al. 2023b). We chose the 3rd order derivative to apply the penalty as we expect some curvature of the temperature profile with pressure, as otherwise this would be suppressed/rejected in our Bayesian analysis. The P–T knots we chose are not evenly spaced, and sample the expected photosphere around 1 bar more closely. Similarly to Line et al. (2015), we include the factor |$\gamma$| as a free parameter in the retrievals. This parameter divides the penalty applied to the 3rd order derivative, and so allows the likelihood penalty to be suppressed at large |$\gamma$| if the data warrants a temperature profile with more oscillatory features. Hence we have 9 free parameters for the temperature profile, with the temperature retrieved at various pressures throughout the atmosphere (see Table 2), and an additional parameter for the P–T penalty factor |$\gamma$|⁠.

2.2.2 Atmospheric opacity

The atmospheric model includes opacity arising from a range of sources. First, we include the line absorption from the molecular species in the atmosphere. This uses a line list to calculate the cross-section for each species. We use the ExoMol line list data base (Tennyson et al. 2016) for H|$_2$|O (Polyansky et al. 2018), HCN (Harris et al. 2006; Barber et al. 2014), and NH|$_3$| (Coles et al. 2019), and the HITEMP data base (Rothman et al. 2010) for CH|$_4$| (Hargreaves et al. 2020) and CO and its isotopologues (Li et al. 2015). We also use the line list from Huang et al. (2013, 2017) for CO|$_2$|⁠. For each pressure and temperature grid point every spectral line is broadened into a Voigt profile (e.g. Gandhi & Madhusudhan 2017) using the H|$_2$| and He pressure broadening coefficients (Gandhi et al. 2020). This is then summed for every line to give the overall cross-section from each species. We retrieve the volume mixing ratio of each chemical species as a free parameter, with vertically constant abundance profiles. This has been shown to be the most optimal setup from previous work (Regt et al. 2024). In addition to molecular line absorption, we also include opacity arising from collisionally induced absorption due to H|$_2$|–H|$_2$| and H|$_2$|–He interactions (Richard et al. 2012).

In addition to the opacity due to gaseous chemical species, we also include a cloud opacity model based on Mollière et al. (2020). The cloud opacity |$\kappa _\mathrm{cl}$| at a pressure P and wavelength |$\lambda$| is given by

with |$P_\mathrm{cl}$| and |$\alpha _\mathrm{cl}$| representing the cloud deck pressure and power law respectively, both of which are set to free parameters for the retrieval. We determine |$\kappa _{\mathrm{cl},0}(\lambda )$| through

where |$\kappa _{\mathrm{cl, 2\, \mu m}}$| and |$\kappa _{\mathrm{cl, 2.5\, \mu m}}$| are free parameters in the retrieval. The opacity at intermediate wavelengths is then interpolated between 2 and 2.5 |$\mu$|m.

2.2.3 Likelihood evaluation

We adopt a method similar to Ruffio et al. (2019) and Gibson et al. (2020) to determine the likelihood (see Regt et al. 2024, for further details). For each order and detector, we construct the covariance matrix |$\mathbf {\Sigma }$| as

with a and l representing the gaussian process amplitude and lengthscale, which are both free parameters in the retrieval. The diagonal elements consist of the overall error |$\sigma _i$| for each data point, and the off diagonal elements give the correlated component. Here,

We also define the arithmetic mean of the variances of elements i and j,

The overall log likelihood for a covariance matrix |$\mathbf {\Sigma }$| is then given by

where N is the number of data points. The residual |$\mathbf {r}$| is defined

where |$\mathbf {d}$| refers to the data, |$\mathbf {m}$| refers to the model, and the scale factor |$\phi$| that is applied to the model for each detector and order except for the first (see Regt et al. 2024). This scale factor is calculated through

The |$s^2$| term represents a scaling of the overall covariance for each order and detector. This is given by

2.2.4 Retrieval setup

We run our retrieval for each night separately, and perform a combined retrieval which analyses the four nights simultaneously. In addition to the chemistry, temperature profile, clouds and GP evaluation, we also retrieve the radius of the companion and the surface gravity, |$\log (\mathrm{g})$|⁠, as free parameters (see Table 2). This is because for our directly-imaged companion the mass and radius uncertainties are quite large as they are based on evolutionary models (Chauvin et al. 2005). We generate our spectral models at a spectral resolution of |$R=300\,000$|⁠, and convolve them down to the instrumental resolution for each night, which ranges between |$R=93\,200$| and |$R=112\,800$| (see Table 1). Once the instrumental broadening has been applied to the model, any additional broadening of the spectrum will be due to the companion’s rotation. We therefore leave the rotational velocity of AB Pic b and the |$\epsilon$| (the limb-darkening coefficient) as free parameters. It is therefore important in this case that the actual spectral resolution of the observations be calculated effectively, as otherwise the rotational broadening will be biased. In addition, we retrieve the radial velocity shift of the spectrum due to Doppler motion of Ab Pic b. Given slight differences in the barycentric velocity as well as wavelength calibration between nights, we retrieve separate velocity shifts for each night in the combined retrieval of the four nights. We apply the barycentric correction for each night on the spectrum to obtain the |$\mathrm{RV}$|⁠, and hence this value represents the radial velocity of the companion in the heliocentric frame.

Altogether, we have 33 free parameters for the retrieval on each individual night, and 51 free parameters for the combined retrieval of the four nights together. To obtain the detection significances of each species we also perform separate retrievals with the species removed. For our analysis we only use the reddest five orders, with wavelengths between |${\sim} 2064$| nm and |${\sim} 2475$| nm, given the heavy telluric contamination in the bluest two orders. We perform our retrievals using the MultiNest Nested Sampling algorithm (Feroz & Hobson 2008; Feroz et al. 2009, 2013; Buchner et al. 2014), with 500 live points.

2.2.5 High-pass filtering

During our retrievals we high-pass filter the data and the model to remove any long range wavelength dependence over the order/detector. This ensures that the retrieved constraints are driven entirely by the individual line features of spectrally active species. We adopted this approach because high-resolution spectroscopy is inherently more reliable in obtaining line ratios and features than the continuum (e.g. Gandhi et al. 2023a; Maguire et al. 2024), particularly so given that the observing conditions varied over our 4 nights of observations (see Table 1). Hence, it is more difficult to determine whether any deviations in the non high-pass filtered spectra between nights are definitively due to variability and not due to changing observing conditions. We therefore apply a boxcar rolling average which is 100 pixels width to remove the long range structure within each order and detector.

The use of the high-pass filter does affect the constraints on a number of atmospheric parameters from our retrieval. First, the radius of the planet is less certain because we cannot use the continuum flux level to determine its value. The radius is therefore constrained indirectly through line shape variations and the log(g) value, which is also coupled to the chemical abundances of each species. As the continuum level is lost, the absolute photospheric temperatures also become more uncertain. The temperature gradients are intrinsically coupled to chemical abundances given that line features can be fit with higher abundances with shallower temperature gradients or lower abundances with more steep gradients. Therefore, the absolute values of the radius, log(g), temperature and abundances become less precise but remain accurate and reliable. However, abundance ratios are very robust given that the line ratios of the many molecular lines in the observations are not altered, and any changes to the temperature, log(g) or radius affect the lines from all of the species in a similar way. Therefore, compositional ratios such as C/O and |$^{12}$|C/|$^{13}$|C will still be tightly constrained even though absolute abundances may not be, as is commonly seen in both ground-based high-resolution spectroscopy and JWST observations with NIRSpec (e.g. Pelletier et al. 2023; Gandhi et al.2023b).

3 RESULTS AND DISCUSSION

In this section, we discuss the constraints on the atmospheric chemistry, temperature profile, cloud deck and other parameters from the retrievals of AB Pic b. We first discuss the variation in the spectra over the nights, then report on the retrievals for each night of observations individually. We then discuss the retrieval over all of the nights combined. Finally, we explore the spin axis and obliquity of AB Pic b from our retrieved projected rotation velocity and the revisited orbital solution.

3.1 Spectral line variations in the data between nights

Fig. 5 shows the spectra of each night with the median spectrum combining all 4 nights subtracted. This shows that there are some differences between nights, particularly for the reddest detectors for night 2. This night had the lowest number of AB pairs and hence the lowest signal-to-noise, which may make the blaze correction more likely to be biased. This night also had the lowest spectral resolution, so it is not surprising that this is the night with the largest deviation from the median. The fact that the most strongly affected regions are in the reddest detector for all of the orders indicates that this effect is instrumental, and does not represent night-to-night atmospheric variations. Given that absolute fluxes and continuum values are challenging with ground-based observations at high resolution, it is more difficult to construct a variability in absolute flux over time such as that for other targets with JWST (e.g. Biller et al. 2024; McCarthy et al. 2024). Therefore, we adopt an approach to look for line depth and shape variations across the nights, which are significantly more stable (see Section 2.2.5). This may limit the degree of variability we are sensitive to, but it should provide the most accurate constraints even though their uncertainty may be higher.

The variations in the line features themselves are small and therefore more difficult to detect. We looked for differences in the spectral lines between each night by subtracting the median from the high-pass filtered data, and performing a cross-correlation against our best-fitting spectral model, shown in Fig. 6. This ensures that we are sensitive to the variations of the individual spectral lines, and the use of this cross-correlation function (CCF) offers a powerful tool to simultaneously explore many thousands of spectral lines at once. For night 3, we see a positive correlation, indicating that the spectral line features are deeper for this night, and vice versa for nights 1 and 4. Therefore, the data do point to variations in the spectra between nights, albeit with a relatively modest |$|\mathrm{CCF}|\lesssim 4$| for each night. Note that if some lines are stronger but others weaker then the overall cross-correlation will be average over these, which may underestimate the variability. On the other hand, using the median-subtracted high-pass filtered data is also prone to errors if there are strong differences in the spectral resolution and seeing between the nights, which will overestimate the variations in the spectra. To obtain a more complete picture of the atmosphere of AB Pic b, we therefore retrieve the various atmospheric parameters for each night.

3.2 Atmospheric retrieval of the individual nights

Fig. 7 shows the best-fitting atmospheric retrieval model against the data for night 1, with the other nights shown in the appendix, Figs. A1, A2 and A3. These show an excellent overall fit to all of the 5 orders we used for our analysis, and therefore that our model is flexible enough to capture the spectral features well. The residuals also indicate that there are no significant lines which are not well fit and hence that the sources of opacity which we have included in our model are sufficient. Note that we have masked regions with heavy telluric contamination, in particular the second order with the broad Brackett |$\gamma$| line from the telluric standard, as discussed in Section 2.1.2.

Fig. 8 and Table 3 show the constraints on the volume mixing ratios for H|$_2$|O and |$^{12}$|CO, and the derived values for the C/O, and the |$^{12}$|C/|$^{13}$|C ratio from the |$^{13}$|CO abundance constraints. These show good consistency across the nights, but there is some variation on the order of |$\sim$|1–2|$\sigma$| for the chemical abundances, and |${\sim} 3\sigma$| for the C/O and isotope ratios. Night 1 shows the lowest abundances for both H|$_2$|O and |$^{12}$|CO, and night 4 shows the highest. All appear to have broadly similar uncertainties, but night 1 has the highest precision given its longest observing time with 11 AB pairs (see Table 1). While the absolute values of the abundances have a wide uncertainty, the C/O and isotope ratios are relatively tightly constrained given that we are more sensitive to the relative line depth between species rather than the absolute value given the high-pass filter applied to the data. The C/O ratio is largely determined by the CO and H|$_2$|O, as these are the dominant species of the atmosphere. This is consistent with chemical models at such temperatures (e.g. Lodders & Fegley 2002; Madhusudhan 2012), which show that CO is the dominant carrier of carbon and CO and H|$_2$|O are the dominant carriers of oxygen at temperatures |$\gtrsim$|1200 K. We included HCN, CH|$_4$|⁠, NH|$_3$|⁠, and CO|$_2$| in our models but these did not show any significant constraints for any of the nights.

3.2.1 ¹²C/¹³C ratio

The |$^{13}$|CO isotopes are strongly detected at |$\gt 4.5\sigma$| for all nights, as shown in Table 3. We verified our detection of |$^{13}$|CO through a CCF (see e.g. Zhang et al. 2021b; Regt et al. 2024), as shown in Fig. 9. This confirms that the detection of |$^{13}$|CO for each night is robust, and follows a similar trend as the detections from the retrievals, with night 3 having the strongest detection and night 2 having the weakest. This is despite the first night, night 1, having a much longer observing time. Such differences in the detection significance and isotope ratio could be driven by patchy clouds in the atmosphere, which may be varying on the time-scale of |$\sim$|1 d (e.g. Tan & Showman 2021; Lee et al. 2024) and affecting the detectability of |$^{13}$|CO. This also backs up the cross-correlation of the median subtracted data in Fig. 6, which showed that night 3 has slightly stronger spectral lines, potentially due to a partially clearer atmosphere with a deeper cloud deck, thereby increasing the line strengths of trace species such as |$^{13}$|CO.

We do see differences in the constrained isotopologue constraints over our four nights of observations (see Fig. 8). Night 3 shows the lowest |$^{12}$|CO/|$^{13}$|CO ratio of all of the nights, |$^{12}$|CO/|$^{13}$|CO |$= 70^{+16}_{-9}$| and night 4 shows the highest, with |$^{12}$|CO/|$^{13}$|CO |$= 120^{+16}_{-15}$|⁠. Fig. 10 also seems to indicate a correlation between the |$^{12}$|CO/|$^{13}$|CO abundance and the cloud deck pressure across the four nights. Night 3 constrains the cloud deck at a higher pressure than the other nights, while simultaneously constraining the lowest |$^{12}$|CO/|$^{13}$|CO abundance. If on night 3 the atmosphere was less cloudy and we were more sensitive to higher pressures where the |$^{13}$|CO lines are generated, the constrained |$^{13}$|CO abundance may be higher and hence reduce the |$^{12}$|C/|$^{13}$|C ratio. At the same time, an increased sensitivity to the higher pressures would also increase the detection significance as we see for this night, and explain the line depths being strongest for this night in Fig. 6. The |$^{13}$|CO would be more susceptible to the effect of clouds given these lines are generated at deeper pressures than the |$^{12}$|CO. The |$\tau =1$| surface histograms in Fig. 11 confirm that |$^{13}$|CO is most sensitive to pressures only |${\lesssim} 1$| dex above the cloud deck, whereas the |$^{12}$|CO is most sensitive to |$\sim$|mbar pressures. However, we should note that different spectral noise realisations and/or the effectiveness of the telluric fit for each of the nights may also contribute to the observed difference to the isotopologue ratios between the nights. Therefore, more observations and modelling are required to definitively assess if this trend is significant and if patchy/broken clouds are present in the atmosphere.

3.2.2 Temperature profile

The retrieved P–T profile for each night is shown in Fig. 11. In the photosphere at pressures |$\sim$|1 bar, the agreement between each model is remarkable, with very tight constraints of |${\lesssim} 10$| K in the temperature at 1|$\sigma$| confidence thanks to the very high signal-to-noise and spectral resolution of the CRIRES+ observations. In the upper atmosphere, we are less sensitive to the temperature as there are fewer spectral lines which are capable of probing these lower pressures, and hence the uncertainty in the temperature increases as we go up the atmosphere. Similarly, going to deeper pressures below the photosphere also increases the uncertainty in the temperature profile. Fig. 11 shows that the uncertainty increases rapidly as we go to higher pressures than the cloud deck, as we are no longer sensitive to these deep layers of the atmosphere in the spectrum. There is however a degeneracy between the photospheric temperature gradient just above the cloud deck and the absolute abundances of H|$_2$|O, |$^{12}$|CO, and |$^{13}$|CO. A steeper temperature profile such as that for night 1, where the temperature increases more quickly going to higher pressures, results in lower abundances of the chemical species. This is because the data can be fit by having a steeper thermal gradient and lowering the abundances correspondingly to match the line features in the spectrum. On the other hand, the shallower temperature gradient for the other nights require higher abundances. However, we do note that this effect is small and remains within |${\sim} 1\sigma$| for each of the nights.

We also compared our P–T profiles to the best-fitting solution from Palma-Bifani et al. (2023) from Exo-REM self-consistent models (Charnay et al. 2018). For pressures |$\lesssim$|1 bar, the self-consistent model agrees remarkably well with the retrieved P–T profile for each night, |$\lesssim$|200 K difference from the median best fit for each night. For higher pressures however, the P–T profiles do show differences, as the self-consistent model becomes adiabatic. This may be driven by several factors. Firstly, we are not sensitive to the very deep pressures (P|$\gtrsim 3$| bar) as they are below the photosphere and below the cloud deck, and hence the uncertainty in the P–T profile increases. However, given our use of the P–T penalty factor for large deviations in the second derivative, the range in P–T profiles in the retrieval may be more limited. Additionally, as we are observing an average over the companion’s atmosphere, patchy clouds/inhomogeneities may alter our result compared to the self-consistent model. Nevertheless, the agreement in the regions with the photosphere is encouraging, and the photospheric temperature gradients, which our retrieval is most sensitive to, shows good consistency.

3.2.3 Clouds

The cloud deck pressure shows some slight variation between the nights (see Fig. 10). Night 3 indicates the highest pressure for the cloud layer, whereas night 4 indicates the cloud is |$\sim$|0.2 dex lower pressure. The cloud pressure shows some correlation with the |$^{12}$|C/|$^{13}$|C ratio, which is lowest for night 3 and highest for night 4, as previously discussed in Section 3.2.1. This also correlates well with the slightly stronger spectral features in night 3 as seen in Fig. 6. However, the differences in cloud deck constraints are relatively small and within |${\sim} 1\sigma$|⁠, and overall the atmosphere of AB Pic b shows good consistency across the 4 nights of observations.

We show the retrieved cloud opacity at 2 and 2.5 |$\mu$|m in Fig. 12. This shows clearly that there is a preference for a cloud deck whose opacity is decreasing with wavelength in the 2–2.5 |$\mu$|m range. For all 4 nights the |$\kappa _{\mathrm{cl, 2\, \mu m}}$| is consistently higher than |$\kappa _{\mathrm{cl, 2.5\, \mu m}}$|⁠. As the data is high-pass filtered, this difference must be driven by the line feature strengths being altered by the presence of the cloud opacity. Our retrievals with this wavelength-dependent cloud opacity were statistically preferred and gave significantly more reliable results over a model with a single wavelength-independent cloud opacity, which constrained very high CO|$_2$| abundances near |$\log (\mathrm{CO_2}) \sim -4$| in order to match the continuum. This highlights that with significant wavelength coverage using a grey cloud model may not always be applicable. Generally, cloud decks originating from Fe, Ti, and Si condensates are expected to have a drop in opacity in this wavelength range as long as the particle size remains |${\lesssim} 1$| |$\mu$|m (e.g. Wakeford & Sing 2015; Pinhas & Madhusudhan 2017), consistent with our result. However, note that other chemical species not included in our retrieval with opacity near 2 |$\mu$|m may also be able to produce this effect.

3.2.4 Other retrieved parameters

Fig. 13 shows the constraints on other retrieval parameters. We find good constraints on the radius of AB Pic b. Each of the nights shows a relatively well constrained radius of |$\sim$|1.85 R|$_\mathrm{J}$|⁠. However, we note that the radius is degenerate with the photospheric temperature and the data is high-pass filtered to remove the influence of the continuum, which may result in biases occurring on its retrieved value as noted in Section 2.2.5). The log(g) shows only a weak constraint near |$\log (\mathrm{g/cms^{-2}}) \sim 4.4$|⁠. This results in a mass constraint from AB Pic b of |$25^{+16}_{-9}$| M|$_\mathrm{J}$| for night 1, but a higher value with a wider constraint of |${\sim} 30^{+20}_{-18}$|M|$_\mathrm{J}$| for the other nights because the log(g) value is less well constrained for nights 2, 3, and 4. We are only weakly sensitive to the log(g) and it is degenerate with the chemical abundances, and this combined with the uncertainties on the radius makes placing significant constraints on the mass challenging. The mass constraints from each night indicate that AB Pic b is more likely a brown dwarf, but with such wide uncertainty it still falls within 1|$\sigma$| uncertainty of planetary-mass objects. Our constrained mass does appear to be slightly higher than previous studies from evolutionary models which found |$\sim$|10 M|$_\mathrm{J}$| (e.g. Bonnefoy et al. 2014), and hence further analyses, potentially without using high-pass filtering on the data, may be required to robustly determine the mass of AB Pic b and other companions from atmospheric retrievals.

We also retrieve the projected rotational velocity (⁠|$v\sin i$|⁠) of AB Pic b. This is achieved by fitting for a rotational broadening kernel on top of the instrumental broadening for each night. This shows that |$v\sin i$| is relatively low at |$\sim$|3.7 km/s. Night 2 does show a higher speed of rotation, and night 1 shows a lower value, but all of them indicate that the spectral lines are not significantly broadened due to the companion’s rotation. This gives AB Pic b one of the lowest projected rotation velocities measured for companions. This is in contrast to other targets (e.g. Snellen et al. 2014; Hsu et al. 2024; Landman et al. 2024; Regt et al. 2024) which showed strong rotational broadening signatures. Such a low |$v\sin i$| may be an indicator that AB Pic b is either almost pole on as we observe it and/or that it has a low rotation speed due to its young age, as discussed in Section 3.4.2. While we do see some variability in our retrieved atmospheric parameters (⁠|$|\mathrm{CCF}|\lesssim 4$|⁠), the lack of substantial variability across the 4 nights of observations may indicate that we are observing the companion along its axis of rotation, as the variability is likely to be lower for such a geometry given the same face is always being presented (Vos et al. 2022). However, we should note that the measurement of the |$v\sin i$| value is highly dependent on the spectral resolution assumed, which in turn is strongly dependent on the seeing given the 0.4 arcsec slit used for our observations. Whilst we have accounted for this to determine the resolution for each night, if the actual spectral resolution is higher than that assumed, our |$v\sin i$| may also be higher. Previous work by Holmberg & Madhusudhan (2022) report a similar spectral resolution to our nights using the 0.4 arcsec slit with AO in the K band for the LTT-9779 system, indicating that our values are consistent with expectations from the performance of CRIRES+ . Other works using the 0.2 arcsec slit have also found higher spectral resolutions due to the adaptive optics reducing the point spread function below the slit width (e.g. Nortmann et al. 2024; Pelletier et al. 2024). In addition, the thermal and pressure broadening of spectral lines can also influence the measured value, particularly when the |$v\sin i$| is so low.

The Doppler velocity measured for the companion, |$\mathrm{RV}$|⁠, varies between night 1 and the other nights, with almost 1 km/s difference. This indicates that our wavelength solution may not be well calibrated to the absolute level. The radial velocity of the star, AB Pic A, is |${\sim} 22.65$| km/s (Soubiran et al. 2018). The four nights do therefore indicate that AB Pic b has a line of sight velocity towards us in its orbit, by |${\sim} 1{\!-\!}2$| km/s. This is consistent with a bound orbit around AB Pic A, but given that we did not observe AB Pic A simultaneously with our observations of AB Pic b, it is difficult to further quantify this. Future work, constraining the velocity shift of both the star and the companion, may offer better insight into the radial velocity of AB Pic b. Measuring the |$\mathrm{RV}$| more accurately also opens up the opportunity to uncover close binary systems, such as the recent discovery of Gliese 229 Ba and Gliese 229 Bb (Xuan et al. 2024b).

3.3 Combined nights

Fig. 14 shows the posterior distributions for the chemical species, the derived C/O and |$^{12}$|C/|$^{13}$|C ratios, and the constrained temperature profile for the retrieval combining all 4 nights of observations. The combined retrieval allows for us to determine the effect of analysing all of the observations assuming the same thermal, chemical and planetary structure. However, we still retrieve different GP amplitudes and Doppler velocities for each night. The results are consistent with all parameters near the average of the 4 nights individually. We also obtain tighter constraints on the parameters than each night on its own. the H|$_2$|O abundance is lower than the |$^{12}$|CO by |$\sim$|0.15 dex, indicating that the companion’s C/O ratio is |$0.59^{+0.01}_{-0.01}$|⁠, consistent with the value measured for the Solar system (Asplund et al. 2021) and the individual nights. The overall abundances of H|$_2$|O and |$^{12}$|CO indicate that the atmosphere is slightly above solar metallicity, with [O/H] = 0.54 |$\pm$| 0.08 and [C/H] = 0.60 |$\pm$| 0.08. Both the overall metallicity and C/O ratio are in excellent agreement with previous works (e.g. Bonnefoy et al. 2010; Palma-Bifani et al. 2023) using grids of self-consistent models.

The P–T profile for our combined retrieval is shown in the right panel of Fig. 14. The profile is in good agreement within |$2\sigma$| of the individual nights at all pressures, but there is a slight downward trend of the temperature at the lowest pressures. This may be caused by slight differences in some of the strongest lines in the spectrum, as these probe the lowest pressures and are therefore the most sensitive in the observations. The 0.1 bar photosphere temperature of |$1714^{+5}_{-5}$| K is however in excellent agreement with each night and previous works (Bonnefoy et al. 2010, 2014; Patience et al. 2012; Palma-Bifani et al. 2023). For almost all pressures our P–T profile is above the CO/CH|$_4$| equilibrium boundary, which is consistent with our strong CO constraints but only upper limits on CH|$_4$| owing to its lower abundance at these high temperatures.

Our cloud deck pressure is also shown in Fig. 14 along with the condensation curves for various refractory species at solar abundance (Visscher et al. 2010; Wakeford et al. 2017). We find that our temperature profile constraint most closely matches the CaTiO|$_3$| condensation curve at the location of the cloud deck. However, we note that given the semigrey cloud model we are using, Fe and Mg|$_2$|SiO|$_4$| clouds are also potential sources of our cloud opacity, especially if present at super-solar abundances where they are expected to condense at higher temperatures.

The carbon isotope ratio shows a constraint of |$^{12}$|C/|$^{13}$|C |$= 102^{+8}_{-8}$|⁠, in line with the average values obtained for each of the nights individually. This value is slightly higher than that for the Solar system (see Fig. 14), but much higher than the value for the ISM. Fig. 15 shows the C/O and isotope ratios for a range of close-in transiting exoplanets, directly-imaged companions and field brown dwarfs. This indicates that AB Pic b has a |$^{12}$|C/|$^{13}$|C ratio in between that observed for other targets. There appears to be no strong trend over the sample as yet, but future studies exploring isotopic ratios in planetary atmospheres from ground and space may yield insights as to the compositional diversity. We should note however that the individual nights of AB Pic b did show some significant differences in the |$^{12}$|C/|$^{13}$|C ratio, ranging from values close to the ISM for night 3 to that significantly greater than the Solar system value for night 4. Therefore, the combined retrieval’s value should be treated with caution as it is likely affected by atmospheric variability, as discussed in Section 3.2.1.

3.3.1 Other chemical species

We find only upper limits for the other chemical species we retrieve. In addition to |$^{13}$|CO, we retrieved the abundance of C|$^{18}$|O, which has recently been detected in the atmospheres of directly-imaged companions (Gandhi et al. 2023b; Xuan et al. 2024c). However, for our retrievals we find |$\log (\mathrm{C^{18}O}) \lesssim -6.7$| at 2|$\sigma$| confidence. As the |$^{16}$|O/|$^{18}$|O abundance ratio is expected to be |${\sim} 500$| (Wilson 1999), it is more difficult to detect given its much weaker features compared with |$^{13}$|CO and hence more observations are needed to constrain this in AB Pic b. We also find no clear constraints for CH|$_4$|⁠, NH|$_3$| or CO|$_2$| (main isotopologues only), with only an upper limits near |$\log (\mathrm{Volume \, Mixing \, Ratio}) \sim -6$| at 2 |$\sigma$| confidence. There is a small hint of HCN with a peak in the posterior, which dominated by the data from night 4. However, given that the other nights showed no indication of HCN, we are unable to conclusively detect the species. The CH|$_4$| abundance for L-type objects at such high temperatures is expected to be low, and our upper limit is consistent with predictions from chemical modelling (Lodders & Fegley 2002; Visscher & Moses 2011). The CO|$_2$| is also relatively unconstrained, as the 2 bluest orders covering |${\sim} 1.9{\!-\!}2.05$| |$\mu$|m covering strong CO|$_2$| absorption were not used in our analysis given the very strong telluric CO|$_2$| at those wavelengths. Hence constraining CO|$_2$| in exoplanetary atmospheres from ground-based observations may prove more challenging, as noted for H|$_2$|O (e.g. Pelletier et al. 2021). We calculate the C/O ratio of the atmosphere using the total carbon content of the atmosphere divided by the total oxygen content, and so the uncertainty on all of the unconstrained trace species is taken into account when determining the C/O ratio. As the H|$_2$|O and CO are significantly more abundant than the upper limits for these chemical species, the uncertainty on the unconstrained other species does not significantly alter the C/O ratio constraints for the combined retrieval or for each night individually.

3.4 Orbit and obliquity

Our derived value of |$v\sin i$| is much lower than that reported in Table 4 of Palma-Bifani et al. (2023) (⁠|$v\sin i \sim 73$| km/s), used in their section 5.1 to derive the spin axis of AB Pic b. The CRIRES+ spectral resolution for each night is much higher than the SINFONI resolution (R|$\sim 3000$|⁠), allowing us to resolve many more individual lines, and we carefully calibrated the wavelength solution and corrected for instrumental broadening (see Section 3.3). Taking this into account and considering that our |$v\sin i$| measurement is stable over the four nights, we decided to re-estimate the spin axis inclination (⁠|$i_\mathrm{ p}$|⁠) of AB Pic b. Additionally, a new SPHERE-IRDIS observation was taken in October 2023 (program 0112.C-2045(A); PI P. Nogueira), which we used to refine the orbital fitting to revisit the orbital inclination (⁠|$i_o$|⁠) of the system. We used these two revised angles to gain further insight into the system’s obliquity, following the method of Bryan et al. (2021).

3.4.1 Revised orbital solution

We present an updated orbital fit for AB Pic b using the OFTI (Orbits for the Impatient) module of orbitize!, a Bayesian rejection sampling method designed for efficiently fitting the orbits of long-period exoplanets (Blunt et al. 2017). For this analysis, we incorporated the data from VLT/NaCo observations in 2003, 2004 (Chauvin et al. 2005), and 2009 (Rameau et al. 2013), as well as the VLT/SPHERE-IRDIS observations from 2015 already used in Palma-Bifani et al. (2023) and a new epoch from 2023. We reporcessed for this work both SPHERE-IRDIS epochs using PACO (Chomez et al. 2023) to extract precise astrometry. We report all the astrometry values in Table 4. In orbitize!, we employed the default priors, as detailed in the documentation, and ran the sampler for 100 000 orbits. We draw the best 50 orbital solutions in Fig. 16. The colourmap in the left panel of this figure indicates the epoch at which AB Pic b is likely to be at the given location with the grey arrow indicating the direction of motion. We can observe that AB Pic b is moving towards its host star. To the right, we zoom in at the location of AB Pic b and see the measured separation and position angle as a function of time together with the best 50 orbits in grey behind.

Adding the new 2023 VLT/SPHERE-IRDIS data point does not generally impact most posterior distributions, meaning most parameters have the same probability density as reported in Palma-Bifani et al. (2023). The impact of this new epoch are a slightly looser constraint on the semimajor axis and a tighter one on the inclination of the orbit. The new semimajor axis is |$307^{+305}_{-115}$| au instead of |$190^{+200}_{-50}$| au, where the errors represent asymmetric 1|$\sigma$| confidence intervals assuming Gaussian distributions. The new orbital inclination is |$98^{+12}_{-5}$| degrees, which suggests a highly inclined orbit but not fully edge-on. We compared the orbital inclination posterior density obtained by Palma-Bifani et al. (2023) with the one obtained in this work, as shown in the first panel of Fig. 17 (in grey and light pink, respectively). In the following section, we use the new inclination distribution for the re-estimation of the obliquity.

3.4.2 AB Pic b inclination and obliquity

To compute the spin axis of the companion, we use the open-source python package exospin.²This tool calculates the probability density distribution of the spin axis (⁠|$i_\mathrm{ p}$|⁠) and obliquity (⁠|$\Psi$|⁠) of a companion. The method for computing the |$i_\mathrm{ p}$| is based on Bryan et al. (2020), following the same strategy already implemented for AB Pic b in Palma-Bifani et al. (2023).

We must extract the inclination angle component from the projected rotational velocity (⁠|$v\sin i$|⁠) measurement to determine the companion’s spin axis. For this, we compute the equatorial rotational velocity of AB Pic b by assuming solid body rotation, using our derived radius distribution from the combined retrieval and the rotational period reported in Zhou et al. (2019) of 2.1 h. Since v and |$v\sin i$| are not independent, we use equation (10) from Bryan et al. (2020) to derive the probability density distribution for |$\cos (i_\mathrm{ p})$|⁠, from which we obtain the distribution of |$i_\mathrm{ p}$|⁠. Our measured |$i_\mathrm{ p}$|⁠, reported in the second panel of Fig. 17, is |$1.8 \pm 1.3$| degrees, indicating that from our line-of-sight AB Pic b is pole-on. This angle is degenerate and symmetric around 90 degrees, so we observe a secondary peak centred at 178.2 degrees.

Finally, we measure the obliquity from the constraints on the orbital and the companion’s inclination. We avoided estimating the true obliquity, as we cannot measure the companion’s orientation in a three-dimensional space. In addition, assuming a uniform distribution for the companion’s longitude of ascending node (⁠|$\Omega$|⁠) to compute equation (11) from Bryan et al. (2020) would only increase the uncertainties of our retrieved obliquity. Therefore, we decided to report only the projected obliquity, which is the difference between the inclination vectors for both the companion and orbit, assuming they lie in the same plane. The measured projected obliquity is |$86.4 \pm 18.6$| degrees, observable in the right panel of Fig. 17 in cyan. Note that this projected obliquity refers to the alignment of the companion’s spin axis to the orbital inclination, not to the star’s spin axis, as is the case when measuring the Rossiter–McLaughlin effect (Triaud 2018).

Our low value for the |$v\sin i$| may be an indication that the orbit of AB Pic b is pole-on. For such a scenario, we would not expect significant modulation of the spectrum that we could attribute to a rotational period (<0.14 per cent for i|$_\mathrm{ p}$| <3 degrees). This may mean the rotation period is higher than that given in Zhou et al. (2019), who find a fast rotation of 2.1 h at significance level of 1.6|$\sigma$|⁠. Such a fast rotation would also mean the companion is very close to its break-up velocity and would be difficult to achieve. Hence, we repeated our analysis of i|$_p$| for two additional rotational periods of 8.4 and 16.8 h. We report the results in the second panel of Fig. 17. If AB Pic b was rotating with a period of 16.8 h, the updated value of i|$_\mathrm{ p}$| would be |$7.9 \pm 5.6$| degrees, and the updated projected obliquity would be |$74.3 \pm 18.5$| degrees. In Palma-Bifani et al. (2023), the authors reported a misalignment of the spin axis of AB Pic b with respect to the orbital normal axis. A significant misalignment could mean AB Pic b is rolling in its orbit and has a Uranus-like orbit and obliquity (see Fig. 1), such as that recently reported for VHS 1256 b (Poon et al. 2024).

On the other hand, another potential scenario is that AB Pic b is intrinsically a slow rotator. The inclination of the orbit is degenerate if the rotational period is not well known, and hence could exceed 16.8 h. Companions such as GQ Lupi b have been shown to be slow rotators potentially due to their young age (Schwarz et al. 2016). AB Pic A is a |$\sim$|13.3 Myr star (Booth et al. 2021), and the large radius that we measure for AB Pic b could be an indication that the planet has recently formed. Using our measured radius, if the companion’s spin axis and orbital inclination are aligned, the orbital period could be as high as |$\sim$|65 h. As the planet contracts over time and potentially accretes more material, it may spin up to higher velocities and its period would reduce. Note however that given our use of the high-pass filtering on the observations the constrained radius of AB Pic b may not be as reliable, and a smaller radius would result in a lower orbital period.

4 CONCLUSIONS

We have observed the directly-imaged companion AB Pic b across 4 consecutive nights with the CRIRES+ spectrograph on the VLT. Due to the sensitivity of ground-based high-resolution observations to line depths and shapes and the difficulty in obtaining accurate continuum fluxes, we focus on the variations in the spectral lines to determine variability. The observations showed variations in the line features over the 4 nights, which were relatively small but cross-correlation of the median subtracted data across all nights against a model template was able to detect the differences in the spectra, with |$|\mathrm{CCF}|\lesssim 4$|⁠. We performed atmospheric retrievals across each of the nights individually as well as combined retrievals which simultaneously analysed all four nights together. These retrievals confirmed that the atmosphere of AB Pic b showed some variability across the nights, with differences of |$2{\!-\!}3\sigma$| uncertainty for a number of parameters.

We constrained the abundance of H|$_2$|O, |$^{12}$|CO, and |$^{13}$|CO in the atmosphere of AB Pic b, and found that there was some variation in the abundances at the 1–2|$\sigma$| level, which was correlated with the temperature gradient at higher pressures. We found a steeper temperature gradient in the deeper atmosphere required a lower abundance of all of the chemical species in order to fit the strength of the observed spectral lines. These differences, while small, do show that constraining absolute abundances from high-resolution spectra is challenging, as seen in other works (e.g. Pelletier et al. 2023; Gandhi et al. 2023a; Maguire et al. 2024).

We also found some variation of the |$^{12}$|C/|$^{13}$|C ratio across the nights. Night 4 constrained a value of |$^{12}$|C/|$^{13}$|C |$= 120^{+16}_{-15}$|⁠, but night 3 indicated a lower value of |$^{12}$|C/|$^{13}$|C |$= 70^{+16}_{-9}$|⁠. This did not seem to be correlated with the detection significances (see Table 3), but was most strongly correlated with the location of the cloud deck. On night 3 the cloud deck was constrained at deeper pressures than for the other nights. This may indicate that the |$^{12}$|C/|$^{13}$|C ratio was lower for this night as the |$^{13}$|CO features in the spectrum were stronger and required a deeper cloud deck. The |$^{13}$|CO is less abundant than |$^{12}$|CO and its spectral features are therefore generated at lower altitudes, hence it is more susceptible to slight variations in the location of the cloud deck as well as the presence of patchy clouds. This shows promise for further studies to investigate and could potentially result in even larger differences in the retrieved |$^{12}$|C/|$^{13}$|C ratios between observations if a companion shows significant variability.

We performed a combined retrieval on all four nights simultaneously. This was a way to assess the effect of combining the different data sets which show some variability. Generally, we found that the values for the parameters were the average of the constraints from the individual nights. The H|$_2$|O and |$^{12}$|CO abundances were slightly above solar composition, in line with previous studies for similar targets (e.g. Xuan et al. 2024a). The C/O ratio was consistent with solar composition, but the |$^{12}$|C/|$^{13}$|C isotope ratio was slightly higher than that of the Solar system and higher than that seen for the ISM (Wilson 1999; Milam et al. 2005). The isotopic abundance may be used as a way to distinguish formation scenarios (Zhang et al. 2021a), and does appear to have significant variations for such companions (Xuan et al. 2024a). Therefore, it is key that we understand the role atmospheric variability plays in determining isotopic ratios. Further observations across a range of Super-Jupiters will also determine whether there are any population level trends across the sample of companions, as well to the sample of field brown dwarfs.

We retrieved the rotational broadening, |$v\sin i$|⁠, of AB Pic b, which we found to be just |$\sim$|3.7 km/s for the four nights. This allowed us to revisit the obliquity measurement from Palma-Bifani et al. (2023). We found that either the |$i_\mathrm{ p}$| is close to 0, indicating that AB Pic b is oriented pole-on towards us, and/or that the companion is a slow rotator. The lack of substantial variability in the spectral lines across our four nights of observations may suggest that the companion’s rotation axis is aligned with our line of sight, as the same face is always presented. Combining this with the orbital inclination of AB Pic b, this may mean that the rotation axis of AB Pic b has a |${\sim} 90^\circ$| misalignment to its orbit and hence a Uranus-like orientation. However, for such a geometry the determination of the period would be more challenging, and Zhou et al. (2019) found a period of |$\sim$|2.1 h but at |$\lt 2\sigma$| confidence. On the other hand, younger companions have been shown to possess intrinsically lower rotation velocities (and therefore higher periods) (Schwarz et al. 2016), and AB Pic A is |${\sim} 13.3$| Myr old so the companion could potentially still be forming. The retrieved radius of |$\sim$|1.9 R|$_\mathrm{J}$| may be an indication that this companion has only recently formed and hence the rotation speed could increase over time as it contracts, although it should be noted that our radius value may not be reliable as the data is high-pass filtered. Therefore, a robust determination of the period of AB Pic b would resolve these two scenarios and open an exciting avenue for further studies into its dynamic history.

Atmospheric variability has been shown to be a key feature of L and T dwarfs, particularly near the L/T transition (e.g. Radigan et al. 2014; Liu et al. 2024). Younger, lower mass objects have shown even more variation in their spectra over time (Vos et al. 2022). Hence, patchy clouds and dynamics are vital to understand and study if we wish to obtain reliable estimates on the chemical composition and other atmospheric properties. High-resolution spectroscopy offers a powerful tool for such studies, enabling Doppler imaging through time-series observations over much shorter time-scales (Crossfield et al. 2014). In addition, a range of general circulation models have shown how cloud radiative feedback can result in photospheric temperatures varying by |$\gtrsim$|100 K (Tan & Showman 2021). In future, higher signal-to-noise observations through instruments such as the ELT’s METIS (Brandl et al. 2021) may help determine the variability of AB Pic b and other companions, opening up the exciting opportunity to study the variation and evolution of planetary atmospheres over time.

ACKNOWLEDGEMENTS

The authors would like to acknowledge the University of Warwick Scientific Computing Research Technology Platform (SCRTP) for assistance in the research described in this paper, performed with the Avon HPC cluster. AZ acknowledges support from ANID – Millennium Science Initiative Program – Center Code NCN2021_080. We thank Alex Sánchez-López for helpful discussions and feedback. We thank the anonymous referee for a careful review of our manuscript.

DATA AVAILABILITY

The data underlying this article are publicly available from the ESO science archive.

Footnotes

REFERENCES

APPENDIX A: BEST-FITTING RESIDUALS FOR OTHER NIGHTS