# One year of AU Mic with HARPS - I. Measuring the masses of the two transiting planets

Norbert Zicher, Oscar Barragán, Baptiste Klein, Suzanne Aigrain, James E. Owen, Davide Gandolfi, Anne-Marie Lagrange, Luisa Maria Serrano, Laurel Kaye, Louise Dyregaard Nielsen, et al.

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

To cite this version:

Norbert Zicher, Oscar Barragán, Baptiste Klein, Suzanne Aigrain, James E. Owen, et al. One year of AU Mic with HARPS - I. Measuring the masses of the two transiting planets. Monthly Notices of the Royal Astronomical Society, 2022, 512, pp.3060-3078. 10.1093/mnras/stac614. insu-03705346

HAL Id: insu-03705346

https://insu.hal.science/insu-03705346

Submitted on 27 Jun 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# One year of AU Mic with HARPS – I. Measuring the masses of the two transiting planets

Norbert Zicher<sup>1,★</sup> Oscar Barragán,<sup>1</sup> Baptiste Klein<sup>1</sup>,<sup>1</sup> Suzanne Aigrain,<sup>1</sup> James E. Owen<sup>1,2</sup>

Davide Gandolfi<sup>1,3</sup> Anne-Marie Lagrange,<sup>4</sup> Luisa Maria Serrano<sup>1,3</sup> Laurel Kaye,<sup>1</sup>

Louise Dyregaard Nielsen<sup>1,5</sup> Vinesh Maguire Rajpaul<sup>1,6</sup> Antoine Grandjean,<sup>4</sup> Elisa Goffo<sup>3,7</sup>
and Belinda Nicholson<sup>1,8</sup>

Accepted 2022 March 2. Received 2022 February 24; in original form 2021 November 16

## ABSTRACT

The system of two transiting Neptune-sized planets around the bright, young M-dwarf AU Mic provides a unique opportunity to test models of planet formation, early evolution, and star-planet interaction. However, the intense magnetic activity of the host star makes measuring the masses of the planets via the radial velocity (RV) method very challenging. We report on a 1-yr, intensive monitoring campaign of the system using 91 observations with the HARPS spectrograph, allowing for detailed modelling of the ~600 m s⁻<fnref n="1" /> peak-to-peak activity-induced RV variations. We used a multidimensional Gaussian Process framework to model these and the planetary signals simultaneously. We detect the latter with semi-amplitudes of K<sub>b</sub> = 5.8 ± 2.5 m s⁻<fnref n="1" /> and K<sub>c</sub> = 8.5 ± 2.5 m s⁻<fnref n="1" />, respectively. The resulting mass estimates, M<sub>b</sub> = 11.7 ± 5.0 M<sub>⊕</sub> and M<sub>c</sub> = 22.2 ± 6.7 M<sub>⊕</sub>, suggest that planet b might be less dense, and planet c considerably denser than previously thought. These results are in tension with the current standard models of core-accretion. They suggest that both planets accreted a H/He envelope that is smaller than expected, and the trend between the two planets' envelope fractions is the opposite of what is predicted by theory.

**Key words:** techniques: radial velocities – techniques: spectroscopic – planets and satellites: fundamental parameters – stars: activity – stars: individual: AU Microscopii – starspots.

## 1 INTRODUCTION

Planets orbiting young stars offer a unique window into the formation and evolution of planets and planetary systems. The first few hundreds of Myr, when the planets evolve most rapidly, and their observable parameters are still affected by initial conditions, offer the most sensitive tests of theoretical models. Planets whose masses and radii can be measured directly are particularly valuable in that respect, but only a few are known to date. This is in large part due to the rapid rotation and enhanced magnetic activity of the host stars, which hinders both transit detection and radial velocity (RV) follow-up.

In recent years, the K2 (Howell et al. 2014) and TESS (Ricker et al. 2015) space missions have enabled the detection of a number of planets transiting young stars in nearby open clusters and associations (including David et al. 2016a,b, 2019a; Mann et al. 2016a,b, 2017, 2018, 2020; Gaidos et al. 2017; Pepper et al. 2017;

Vanderburg et al. 2018; Rizzuto et al. 2018, 2020; Newton et al. 2019, 2021; Kossakowski et al. 2021), but most of them currently lack mass determinations. Dedicated spectroscopic surveys have also uncovered a number of non-transiting planets around very young stars (including Quinn et al. 2012; Donati et al. 2016; Yu et al. 2017), but these lack radius estimates. The combination of measured radius and mass gives access to the planet's mean density, which in turn contains clues regarding its bulk and atmospheric composition. Specifically, young planets (≤100 Myr) with orbital periods ≤30 d and known mean densities are particularly valuable, as this is the period range where planets are most strongly affected by atmospheric mass-loss over their lifetimes (e.g. Lopez & Fortney 2013; Owen 2019). Therefore, young planets provide the opportunity to probe the amount of primordial gas these planets accreted from their parent protoplanetary discs before the majority of it is lost. In addition, for those planets that are either close enough to their star, or have low enough mass that mass-loss is important on time-scales comparable to their age then constraints can be placed on their initial entropy (Owen 2020). Knowledge of the amount of primordial gas these

<footnote n="1"><sup>1</sup>Sub-department of Astrophysics, Department of Physics, University of Oxford, Oxford OX1 3RH, UK</footnote>

<footnote n="2"><sup>2</sup>Astrophysics Group, Blackett Laboratory, Imperial College London, London SW7 2AZ, UK</footnote>

<footnote n="3"><sup>3</sup>Dipartimento di Fisica, Università degli Studi di Torino, Via Pietro Giuria 1, I-10125 Torino, Italy</footnote>

<footnote n="4"><sup>4</sup>Université Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="5"><sup>5</sup>Geneva Observatory, University of Geneva, Chemin des Mailettes 51, CH-1290 Versoix, Switzerland</footnote>

<footnote n="6"><sup>6</sup>Astrophysics Group, Cavendish Laboratory, University of Cambridge, J. J. Thomson Avenue, Cambridge CB3 0HE, UK</footnote>

<footnote n="7"><sup>7</sup>Thüringer Landessternwarte Tautenburg, Sternwarte 5, D-07778 Tautenburg, Germany</footnote>

<footnote n="8"><sup>8</sup>Centre for Astrophysics, University of Southern Queensland, QLD 4350 Toowoomba, Australia</footnote>

<footnote>* E-mail: norbert.zicher@physics.ox.ac.uk</footnote>
planets accreted, and in what thermodynamic state, is extremely valuable as the origin of close-in planets remains uncertain (e.g. Bean, Raymond & Owen 2021). To date, there are only a few planets with well-determined masses and radii orbiting stars with well-determined ages below 1 Gyr (Mann et al. 2017; Barragán et al. 2019b, 2022a; David et al. 2019b; Plavchan et al. 2020; Klein et al. 2021; Kossakowski et al. 2021; Suárez Mascareño et al. 2021).

Two of them, K2-100b (Mann et al. 2017; Barragán et al. 2019b) and TOI-1201b (Kossakowski et al. 2021), orbit stars with ages in the range 600–800 Myr. K2-100b is a $3.88 \pm 0.16 R_⊕$, $21.0 \pm 6.2 M_⊕$ planet on a $1.67 d$ orbit around an early G-dwarf in the ~750 Myr-old Praesepe open cluster, whose transits were discovered in data from the K2 space mission (Howell et al. 2014), and whose mass was measured using HARPS-North (Pepe et al. 2010). Photoevaporation models suggest that K2-100b is likely still losing its atmosphere (see Barragán et al. 2019b, and references therein), but any estimate of both its original and its final composition depends critically on assumptions regarding the UV flux of the host star and how it has evolved throughout its lifetime. TOI-1201b is a $2.42 \pm 0.09 R_⊕$, $6.3 \pm 0.9 M_⊕$ planet on a $2.49 d$ orbit around an early M-dwarf whose gyrochronological age lies in the range 600–800 Myr and which is a likely kinematic member of the Hyades supercluster. Its transits were discovered in data from the TESS space mission (Ricker et al. 2015), and its mass was measured using CARMENES (Quirrenbach et al. 2014). Both host stars are relatively faint in the optical, so additional follow-up to refine the planets' masses or detect their atmosphere is best pursued with infrared instruments.

By contrast, another of these planets orbits AU Mic, the brightest M-dwarf in the Southern sky ($V = 8.73$ Torres et al. 2006). Plavchan et al. (2020, hereafter P20) discovered the transits of AU Mic b in the first month of data from the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015), reporting an orbital period of $8.46 d$ and a planetary radius of $4.2 \pm 0.2 R_⊕$. Due to its young age ($22 \pm 3$ Myr; Mamajek & Bell 2014), the AU Mic system provides a much more direct test of planet formation and evolution models, but its intense magnetic activity also exacerbates the challenges of detecting both the transits and the RV signal of any planetary companions. In addition, AU Mic hosts an edge-on debris disc (e.g. Kalas, Liu & Matthews 2004; Boccaletti et al. 2015), which shows peculiar fast-moving features, still unexplained (Boccaletti et al. 2018). The TESS light curve exhibits significant variability (~0.1 mag) due to star-spots, from which P20 determined the star's rotation period, $P_{rot} = 4.863 d$, as well as frequent flares, which hinder the detection and modelling of the transits. P20 combined data from several optical and near-infrared spectrographs, obtained over several years (starting long before the discovery of the transits) in an effort to constrain the planet's mass, but obtained only an upper limit. A more intensive observing campaign using the SPIRou near-infrared spectropolarimeter (Donati et al. 2020), focusing on several consecutive rotational cycles, enabled Klein et al. (2021, hereafter K21) to measure a mass of $17.1_{-4.5}^{+4.7} M_⊕$ for AU Mic b, recently refined at $20.1_{-1.6}^{+1.7} M_⊕$ in Cale et al. (2021).

P20 also reported an additional transit-like event from a possible second transiting planet. TESS reobserved the system 2 yr later, enabling Martioli et al. (2021, hereafter M21) to confirm the detection of AU Mic c, with a radius of $3.24 \pm 0.16 R_⊕$ planet and a period of $18.859 d$, though no more than a $5\sigma$ upper limit of $20.13 M_⊕$ has been reported for its mass (Cale et al. 2021). Note that the existence of more than one transiting planet around AU Mic is not altogether surprising, given that AU Mic hosts an edge-on debris disc, and that observations of the transits of AU Mic b have shown its orbit to be

aligned (including Hirano et al. 2020; Martioli et al. 2020; Addison et al. 2021; Szabó et al. 2021). The parameters of AU Mic and its two transiting planets that are relevant in this work, as reported in the literature at the time of writing, are listed in Table 1.

Together with AU Mic's youth and proximity, the fact that it hosts not one but two transiting planets opens up the possibility of comparative exoplanetology. It is thus the target of a number of ongoing ground- and space-based observing campaigns aiming to measure the planets' masses more precisely, as well as to detect the planet's atmosphere and exosphere in transmission.<fnref n="1" /> Knowing the planet's mass is important not only to understand its bulk composition, but also to interpret any transmission spectrum observations, as it helps break degeneracies in the atmospheric retrieval process (Batalha et al. 2019).

In this work, we report on the results of a 10-month intensive monitoring campaign on AU Mic with HARPS, which aims to characterize the activity-induced RV signal of the host star sufficiently to detect the two transiting planets. We use HARPS because activity signals are well characterized in the optical, where they are significantly strong signals. As our approach relies on modelling the effects of stellar activity alongside the planetary signal, rather than filtering one to reveal the other, working in a regime where activity signals are strong can be an advantage rather than an impediment.

This manuscript is part of a series of papers under the project GPRV: Overcoming stellar activity in radial velocity planet searches funded by the European Research Council (ERC, P.I. S. Aigrain), and will be followed by two companion papers (namely Klein et al. 2022, and Zicher et al., in preparation). This paper is structured as follows. Section 2 gives details of the observations and data reduction, as well as the identification of epochs affected by stellar flares, which were discarded from the rest of the analysis. We describe the framework we use to model the RVs and activity indicators in Section 3.1, where we detail the simulations done, using the same framework, to optimize the observing strategy ahead of time. We discuss our results and confront them with theoretical models of early planet evolution in Section 4. Finally, we conclude and discuss future prospects for the characterization of planetary signals orders of magnitude below the stellar noise in Section 5.

## 2 RV OBSERVATIONS

AU Mic has been monitored by a number of optical and near-IR RV spectrographs over the past decade, both before and after the discovery of the transiting planets, including HARPS, iShell, HIRES, CARMENES, SPIRou, CHIRON, IRD, TRES, Minerva Australis (P20, K21; Cale et al. 2021, hereafter C21). However, none of the individual data sets so far have sufficient time-sampling to constrain both the activity signals and the orbits of the transiting planets on their own. Combining data from multiple instruments is possible, and can in principle be advantageous, if the chromaticity of the activity signal is explicitly accounted for when modelling the data (C21). However, combining data from multiple instruments with different wavelength ranges is also challenging because each instrument has different zero-point and noise characteristics, as well as qualitatively distinct activity indicators. Furthermore, the sparse time-sampling of many of these data sets compared to the rotation period of AU Mic, as well as the orbital period of both planets, makes them of limited use to constrain the planetary masses and orbits. The data sets with the best time-sampling to date were those presented in C21, but even

<footnote n="1">¹See e.g. HST Programme 15836, PI E. Newton.</footnote>
<table><caption>**Table 1.** Parameters of the AU Michost star and transiting planets gathered from the literature.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Reference/Notes</th></tr></thead><tbody><tr><td colspan="3">Au Mic stellar parameters</td></tr><tr><td>Distance</td><td>9.7248 ± 0.0046 pc</td><td>Gaia Collaboration (2021)</td></tr><tr><td>T<sub>eff</sub></td><td>3700 ± 100 K</td><td>Afram & Berdyugina (2019)</td></tr><tr><td>Radius (R<sub>*</sub>)</td><td>0.75 ± 0.03 R<sub>☉</sub></td><td>White et al. (2019)</td></tr><tr><td>Mass (M<sub>*</sub>)</td><td>0.50<sup>+0.03</sup><sub>-0.03</sub> M<sub>☉</sub></td><td>P20</td></tr><tr><td>log g</td><td>4.39 ± 0.03</td><td>Computed from R<sub>*</sub> & M<sub>*</sub></td></tr><tr><td>Luminosity</td><td>0.09 ± 0.02 L<sub>☉</sub></td><td>Plavchan et al. (2009)</td></tr><tr><td>Age</td><td>22 ± 3 Myr</td><td>Mamajek & Bell (2014)</td></tr><tr><td>P<sub>rot</sub></td><td>4.86 ± 0.01 d</td><td>P20</td></tr><tr><td>i<sub>orb</sub></td><td>89.5 ± 0.4°</td><td>P20</td></tr><tr><td>v sin i</td><td>7.8 ± 0.3 km s<sup>-1</sup></td><td>K21</td></tr><tr><td>Linear limb dark. coef</td><td>0.2348</td><td>Claret (2018)</td></tr><tr><td>Quadratic limb dark. coef</td><td>0.3750</td><td>Claret (2018)</td></tr><tr><td colspan="3">Au Mic b</td></tr><tr><td>Transit epoch T<sub>0</sub> (BJD<sub>TDB</sub> −2450000)</td><td>8330.39051 ± 0.00015 d</td><td>M21</td></tr><tr><td>Period</td><td>8.463000 ± 0.000002 d</td><td>M21</td></tr><tr><td>Semimajor axis</td><td>0.0645 ± 0.0013 au</td><td>Kepler's law</td></tr><tr><td>Impact parameter</td><td>0.18 ± 0.11 R<sub>*</sub></td><td>M21</td></tr><tr><td>Planet radius</td><td>4.07 ± 0.17 R<sub>⊕</sub></td><td>M21</td></tr><tr><td>Velocity semi-amplitude</td><td>8.5<sup>+2.3</sup><sub>-2.2</sub> m s<sup>-1</sup></td><td>K21</td></tr><tr><td>Planet Mass</td><td>17.1<sup>-4.5</sup><sub>+4.7</sub> M<sub>⊕</sub></td><td>K21</td></tr><tr><td>Equilibrium temperature</td><td>593 ± 21 K</td><td>M21</td></tr><tr><td colspan="3">Au Mic c</td></tr><tr><td>Transit epoch T<sub>0</sub> (BJD<sub>TDB</sub> −2450000)</td><td>8342.2223 ± 0.00015 d</td><td>M21</td></tr><tr><td>Period</td><td>18.859019 ± 0.000016 d</td><td>M21</td></tr><tr><td>Semimajor axis</td><td>0.1101 ± 0.0022 au</td><td>Kepler's law</td></tr><tr><td>Impact parameter</td><td>0.51 ± 0.21 R<sub>*</sub></td><td>M21</td></tr><tr><td>Planet radius</td><td>3.24 ± 0.16 R<sub>⊕</sub></td><td>M21</td></tr><tr><td>Velocity semi-amplitude</td><td>0.6 < K<sub>c</sub> < 9.5 m s<sup>-1</sup></td><td>K21</td></tr><tr><td>Planet Mass</td><td>2.2 < M<sub>c</sub> < 25 M<sub>⊕</sub></td><td>M21</td></tr><tr><td>Equilibrium temperature</td><td>454 ± 16 K</td><td>M21</td></tr></tbody></table>

in that study, the largest number of nights on which AU Mic was observed by a given instrument in a given season was about 40. We therefore opted to carry out a dedicated, intensive survey using a single instrument, namely HARPS, and focused exclusively on this new data set in this work. We defer any attempt to combine the new data with archival data sets to future work.

## 2.1 Previous HARPS observations

The system has been observed since late 2013 with HARPS (under programmes 192.C-0224, 098.C-0739, 099.C-0205, 0104.C-0418, PI Lagrange), as part of a wider, long-term monitoring survey aiming to detect Jovian planets at intermediate separations (Grandjean et al. 2020) around targets of the SPHERE direct imaging survey (Beuzit et al. 2019). In total there are 54 archival spectra taken on 29 individual nights between 2013 October 2 and 2019 November 2. As these observations were primarily intended to search for giant planets beyond the snow line, their time sampling is too sparse to constrain the orbits of the transiting planets effectively, but we extracted the RVs and activity indicators from the ESO archive when planning our own programme, in order to evaluate the magnitude of the activity-induced RV variations and to estimate the number of new observations required. We note that neither this HARPS programme, nor the SPHERE survey have reported any planet detections at wide separations around AU Mic so far (Lannier et al. 2017). In future work, we plan to combine both new and archival HARPS data sets to place limits on the presence of additional planets in the system (Zicher et al. in preparation).

## 2.2 New HARPS observations

The new HARPS observations were collected during ESO Periods 106 and 107 under programmes 0105.C-0288² and 0106.C-0852 (PIs Aigrain and Zicher). In total, we obtained 91 individual spectra: 22 observations on 20 individual nights in P106 (between 2020 November 15 and 2020 December 9), and 69 observations on 49 individual nights in P107 (between 2021 May 24 and 2021 September 22). All observations were taken in High-Accuracy Mode (HAM), with resolution R = 115 000 (Mayor et al. 2003). The exposure time was initially set to 900 s, but increased to 1200 s as the seeing became more variable from 2021 June (onset of Chilean winter), resulting in a typical signal-to-noise ratio (SNR) at 550 nm of 80 to 100. Both programmes were carried out under the auspices of the HARPS time-share programme organized by F. Bouchy and X. Dumusque.

## 2.3 Data reduction and time-series extraction

The data were reduced using version 3.8 of HARPS Data Reduction Software (DRS). DRS RVs were extracted from the Cross-Correlation Function (CCF) obtained using an M4 mask. We also extracted RVs using two alternative pipelines which have been shown to outperform the DRS for M-dwarfs: SERVAL (Zechmeister et al. 2018) and HARPS-TERRA (Anglada-Escudé & Butler 2012). SERVAL and HARPS-TERRA both compute a template spectrum from

<footnote n="2">²The P105 observations were delayed because of the COVID-19 pandemic and taken in P107.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Mean-subtracted RV time-series as extracted with SERVAL (red circles), TERRA (blue diamonds), and the DRS (yellow squares). Note the gap in the x-axis corresponding to a ~5-month break in the observations. The vertical lines show the observations affected by flares (purple dashed) and clouds (green dotted), which are also marked by 'x' symbols, and were discarded from the rest of the analysis (see Section 2.4 for details).</figcaption></figure>

<table><caption>**Table 2.** Radial velocity observations and activity indicator time-series. Observations which were affected by flares or clouds are flagged in the 'Remarks' column. The full version of this table is available in machine-readable format as part of the supplementary material.</caption><thead><tr><th>Time<br>BJDTDB - 2450000</th><th>RV<sub>SERVAL</sub><br>km s<sup>-1</sup></th><th>σ<sub>RV, SERVAL</sub><br>km s<sup>-1</sup></th><th>RV<sub>TERRA</sub><br>km s<sup>-1</sup></th><th>σ<sub>RV, TERRA</sub><br>km s<sup>-1</sup></th><th>DLW<br>1000 m<sup>2</sup> s<sup>-2</sup></th><th>σ<sub>DLW</sub></th><th>S<sub>HK</sub></th><th>σ<sub>SHK</sub></th><th>SNR<br>@550 nm</th><th>Remarks</th></tr></thead><tbody><tr><td>7223.643264</td><td>0.1360</td><td>0.0060</td><td>0.2258</td><td>0.0059</td><td>-12.2629</td><td>5.5325</td><td>7.8457</td><td>0.0237</td><td>71.2</td><td></td></tr><tr><td>7223.652882</td><td>0.1333</td><td>0.0058</td><td>0.2252</td><td>0.0058</td><td>-3.2790</td><td>4.8515</td><td>7.7974</td><td>0.0238</td><td>70.6</td><td></td></tr><tr><td>7333.530396</td><td>0.0065</td><td>0.0024</td><td>0.1057</td><td>0.0023</td><td>-40.2712</td><td>4.7079</td><td>7.8161</td><td>0.0160</td><td>129.1</td><td></td></tr><tr><td>7333.541066</td><td>0.0110</td><td>0.0024</td><td>0.1066</td><td>0.0025</td><td>-41.6257</td><td>4.8444</td><td>7.7178</td><td>0.0168</td><td>122.3</td><td></td></tr><tr><td>7493.886516</td><td>-0.0004</td><td>0.0040</td><td>0.1059</td><td>0.0048</td><td>43.9644</td><td>3.3771</td><td>8.8909</td><td>0.0134</td><td>146.9</td><td></td></tr><tr><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td></td></tr></tbody></table>

the observations themselves, rather than using a pre-existing digitized mask to compute the CCF as in the DRS. In the case of SERVAL, the template is constructed by co-adding the observed spectra in the barycentric rest frame. In the case of HARPS-TERRA, the template is initially computed by co-adding high-SNR observations, allowing preliminary estimates of the RVs to be made, then a new, higher SNR template is produced by co-adding all the spectra in the stellar rest frame, and the RVs are evaluated again. The median RV uncertainties for the DRS, SERVAL, and HARPS-TERRA were 5.2, 3.6, and 3.4 m s⁻¹, respectively, while the Root Mean Square (RMS) RV variations are 121, 130, and 132 m s⁻¹, respectively. Even though all three time-series appear qualitatively similar (see Fig. 1), we observe a much better agreement between SERVAL and HARPS-TERRA RVs (median difference of ~5 m s⁻¹), than between either and the DRS (median difference between HARPS-TERRA and DRS RVs ~20 m s⁻¹). As outlined above, this is not surprising, as HARPS-TERRA and SERVAL both construct a template from the observed spectra, whereas the DRS M4 mask is optimized for older, less active and somewhat cooler stars than AU Mic.

We also extracted a number of ancillary time-series which might serve as activity indicators. These fall into two categories: those that quantify changes in the mean line profile, and those that measure chromospheric emission in the cores of certain spectral lines. In the former category, we include the full width at half-maximum (FWHM) and Bisector Inverse Slope (BIS) of the CCF provided by the DRS, as well as the Differential Line Width (DLW) computed by the SERVAL pipeline. The DLW is obtained by correlating the fit residuals with the second derivative of the template (Zechmeister et al. 2018). The chromospheric indicators we consider include the

Mt Wilson $S_{HK}$ and $H_α$ indices computed by the HARPS-TERRA pipeline,<fnref n="3" /> which measure emission in the cores of the Ca II H and K and the $H_α$ lines, respectively. As presented in Sections 2.4 and 3.3, the chromospheric activity indicators were used to identify observations affected by flares, while the line-shape indicators were most useful to model the activity-induced signals in the RVs.

All the HARPS RVs and activity indicators after 2015 June (after the HARPS fibre upgrade (Lo Curto et al. 2015) are given in Table 2 – the full version of which is available in machine-readable format in the supplementary online material. The activity indicator time-series are shown in Fig. 2.

## 2.4 Removing flares and clouds

Along with displaying significant rotational modulation of star-spots, the ground-based and TESS light curves of AU Mic show frequent white-light flares (e.g. Hebb et al. 2007, P20, M21). These flares can significantly affect the measured radial velocities as they distort the line profiles. As our current models are not equipped to account for the RV effects of flares, we flagged the observations affected by flares and discarded them from the analysis. The procedure used to identify observations strongly affected by flares is described in detail in the companion paper Klein et al. (2022). In short, we computed a chromospheric emission metric for seven different chromospheric lines, namely Ca II H and K (resp. 3968.47 and

<footnote n="3">³Note that the chromospheric indicators provided by the SERVAL pipeline and the DRS are entirely consistent with those provided by HARPS-TERRA, with Pearson correlation coefficients systematically larger than 0.99.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Activity indicators versus time (left) and versus the RVs produced by the corresponding pipeline (right). From top to bottom: DRS FWHM and BIS, SERVAL DLW, HARPS-TERRA $S_{HK}$, and $H\alpha$ index. The colour indicates the pipeline used to extract each indicator: blue for HARPS-TERRA, red for SERVAL, and yellow for the DRS. Significant outliers in either $S_{HK}$ or $H\alpha$, which indicate that the corresponding observations were affected by flaring, are marked by purple dashed vertical lines in the left column. Along with one observation affected by clouds (green dotted vertical line), these observations (marked as ‘x’ symbols in all panels) were discarded from the rest of the analysis (see Section 2.4 for details).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Lomb–Scargle periodograms of the three RVs time-series (top), DLW and FWHM (second panel), BIS (third panel), S-index and $H_α$ (fourth panel), and the window function of the observations (bottom). The vertical crimson lines denote the stellar rotation period and its first two harmonics; the vertical green lines denote 1-d aliases of the stellar rotation period and the vertical black dotted lines denote the orbital periods of AU Mic b and c.</figcaption></figure>

3933.66 Å), $H_α$ (6562.808 Å), $H_β$ (4861.363 Å), Na I D1 and D2 (resp. 5895.92 and 5889.95 Å), and He I D3 (5875.62 Å) using the method of Zechmeister et al. (2018) and the integrations windows of Gomes da Silva et al. (2011). As the effects of stellar flares might change from one chromospheric line to the other, we defined a master index by taking the median-normalized average of all chromospheric emission indices. We finally applied a 3σ-clipping process to the resulting time-series, flagging a total of six observations. These were deemed likely to be affected by a stellar flare, and were removed from the subsequent analysis (see Fig. 2). Note that the RVs affected by flares do not stand out as outliers in the raw RVs, but would do so after subtracting the activity model described in Section 3.1.

After removing the observations affected by flares, one substantial outlier (4 − σ) remained when looking at the residual time-series (after subtraction of the best-fitting stellar activity plus planetary signals model, as described in Section 3). We inspected images from the Danish all-sky camera at La Silla, which showed that variable, high cirrus were present that entire night. We then inspected all three observations of AU Mic taken that night, and found the spectra to be strongly contaminated by the Moon. We thus removed them from further analysis. Although only one of the three observations taken that night was a significant RV outlier, the spectra were all affected by the Moon so we deemed it safer to remove them all. We checked that no other spectra taken on other nights were affected in the same

way by dividing each spectrum by the median spectrum across all observations. The three observations from that night were the only ones which showed strong evidence of Moon contamination.

All the observations affected by flares or clouds are indi-
cated in Table 2. The total number of epochs remaining and
included in our final analysis, shown in Fig. 1, was therefore
82.

2.5 Correlations and periodograms

In the right-hand column of Fig. 2, we show the correlation between each activity indicator and the RVs extracted using the corresponding pipeline. Each plot is annotated with the corresponding Pearson correlation coefficient. Although all the time-series are dominated by activity signals, this does not result in a clear correlation between the activity indicators and the RVs, with the exception of the BIS, which is anticorrelated with the RVs. This is fully expected, as the other activity indicators depend only on the projected area of the active regions, their contrasts, and (for FWHM and DLW) on the absolute value of the surface RV at the location of the active regions relative to the stellar rest frame, but the RVs and BIS also depend on the sign of this quantity.

Lomb–Scargle periodograms of the RVs and activity indicators
time-series are shown in Fig. 3, along with that of the observation
window function. Both RV and DLW time-series show strong peaks at the rotation period and its first harmonic. There is also a peak in both periodograms near, but not exactly at, the second harmonic ($P_{rot}/3$). Peak splitting around the rotation period and its harmonics can indicate differential rotation or active region evolution (see e.g. Aigrain, Pont & Zucker 2012, K21). The raw RVs do not display significant peaks at the orbital period of either planet. There are no prominent peaks in the window function, although there is a little excess power around, but not exactly at, the period of AU Mic c. These periodograms were calculated using the `astropy` package (Astropy Collaboration 2013, 2018).

## 3 MODELLING THE RVS

In this section, we describe the method used to analyse the data, which involved jointly modelling activity and planet signals in the RVs and the activity indicators simultaneously, using the multidimensional Gaussian Process (GP) model introduced by Rajpaul et al. (2015) as implemented in the `pyaneti` https://github.com/oscaribv/pyaneti. package (Barragán et al. 2022b, hereafter B22).

### 3.1 Activity model

We model the activity signals in the RVs, $V(t)$, and one or two activity indicators, $\alpha(t)$ and optionally $\beta(t)$ as

$$ 
\begin{aligned}
V(t) &= A_V G(t) + B_V \dot{G}(t) + C_V \\
\alpha(t) &= A_{\alpha} G(t) + C_{\alpha}, \\
\beta(t) &= A_{\beta} G(t) + B_{\beta} \dot{G}(t) + C_{\beta},
\end{aligned}
\quad (1) 
$$

where the function $G(t)$ is a latent variable, loosely representing the projected area of the visible hemisphere covered by active regions, and is modelled as a quasi-periodic GP with covariance function:

$$ \gamma(t_i, t_j) = \exp \left[ - \frac{\sin^2[\pi(t_i - t_j)/P_{GP}]}{2\lambda_p^2} - \frac{(t_i - t_j)^2}{2\lambda_c^2} \right], \quad (2) $$

where $P_{GP}$ is the characteristic period of the GP, and corresponds to the stellar rotation period, $\lambda_p$ the inverse harmonic complexity, and $\lambda_c$ is the evolution time-scale, and is related to the lifetime of the active regions. The parameters $A_V, B_V, A_\alpha, A_\beta$, and $B_\beta$ are free parameters, which control the relationship between the latent GP variable and the observables, while $C_V, C_\alpha$, and $C_\beta$ are constant offsets for each time-series, which are also free parameters. The interested reader is referred to Rajpaul et al. (2015) and B22 for a full description of the activity model.

Chromospheric activity indicators, such as $S_{HK}$ or $H_\alpha$, depend primarily on the fraction of the visible disc that is covered in active regions, and are thus expected to be approximately proportional to $G(t)$ and were thus modelled under the form $\alpha(t)$ in equation (1). The same is true of the FWHM and DLW, which measure the width of the spectral lines. On the other hand, the BIS, which measures the asymmetry of the spectral lines, is expected to depend on $G(t)$, and was thus modelled like $\beta(t)$ in equation (1).

### 3.2 Joint activity and planet model

To model the activity signal alongside the signals of the known transiting planets, we subtract from the RVs the sum of two Keplerian signals, and compute the likelihood of the GP applied to the residuals (see equation 7 in B22). This introduces a further 6 or 10 free parameters into the model (depending on whether the orbits are assumed to be circular or not).

The measurement uncertainties on the RVs and DLWs at each epoch were accounted for by adding a term to the diagonal of the GP covariance matrix. In addition, we add a separate 'jitter' term to the diagonal of the covariance for each time-series (in effect, a constant term added in quadrature to the formal uncertainties on each observation, see equation 13 in B22). This term absorbs any imperfections in our activity plus Keplerian models and ensures that the resulting uncertainties are propagated to the final parameter estimates.

### 3.3 Choice of RV version and activity indicator(s)

We performed a number of runs using different versions of the RV time-series, and different combinations of activity indicators. We focused on the DRS FWHM & BIS and the SERVAL DLW because their periodograms (shown in Fig. 3) display prominent peaks at the rotation period of AU Mic and its harmonics. The $S_{HK}$ of the chromospheric activity index has a more complex periodogram, and proved to be less useful in constraining the activity signal in the RVs. For the RV time-series, we tried using both the DRS and the SERVAL versions, but we note that the HARPS-TERRA RVs would give results that are essentially identical to those obtained with the SERVAL RVs.

We tried the following models:

1) $V(t) = RV(DRS)$, $\alpha(t) = FWHM$ and $\beta(t) = BIS$;

2) $V(t) = RV(SERVAL)$, $\alpha(t) = DLW$ and $\beta(t) = BIS$;

3) $V(t) = RV(SERVAL)$ and $\alpha(t) = DLW$;

All three models give mutually consistent results, although model 1 does not yield a $3\sigma$ detection of AU Mic c, whereas models 2 and 3 do. In other words, the SERVAL RVs and DLWs result in slightly better sensitivity to the planetary signals compared to the DRS RVs and FWHMs. The results of models 2 and 3 are essentially equivalent. In principle, including the BIS in the modelling should help constrain $G(t)$, but in this specific case the BIS does not seem to provide additional information that the RV and DLW time-series do not already contain. For the remainder of this paper, we therefore adopt model 3 as our fiducial model, since it is simple and has three fewer free parameters. For completeness, we show the results for the three runs in Appendix B.

### 3.4 Exploration of the parameter space

The full model has 20 parameters: five per Keplerian signal, eight for the activity model, and two jitter parameters. These are listed in Table 3, which gives the prior adopted for each parameter. The period and time of transit for each planet are tightly constrained by the TESS light curve and we adopted Gaussian priors based on the ephemerides reported by M21. We also adopt a uniform prior between 4.8 and 4.9 d for the GP period $P_{GP}$, based on the rotation period reported by P20. For the orbital eccentricities of the planets, we adopt a beta distribution prior, as advocated by Van Eylen et al. (2019) for multiplanet systems. For all the other parameters, we adopted minimally informative priors.

We then explored the parameter space with a Markov chain Monte Carlo (MCMC) sampler (see Barragán, Gandolfi & Antoniciello 2019a, for details) to evaluate the joint posterior distribution over all the parameters. We iterated 250 independent Markov chains in sets of 5000 steps, checking for convergence using the autocorrelation length of the chains after each set of 5000 steps. If the chains were not converged, we iterated for a further set of 5000 steps, repeating the process until convergence was reached. We then used the last
<table><caption>**Table 3.** Priors and posterior median and confidence intervals for the parameters of the joint RV and DLW model. The posterior distributions from which we infer these parameters are shown in Fig. A1.</caption><thead><tr><th>Parameter</th><th>Prior<sup>a</sup></th><th>Posterior value<sup>b</sup></th></tr></thead><tbody><tr><td colspan="3">AU Mic b's parameters</td></tr><tr><td>Orbital period P<sub>orb</sub> (d)</td><td>N[8.463000, 0.000002]</td><td>8.463000 ± 0.000002</td></tr><tr><td>Transit epoch T<sub>0</sub> (BJD<sub>TDB</sub> −2450000)</td><td>N[8330.39051, 0.00015]</td><td>8330.39051 ± 0.00015</td></tr><tr><td>Eccentricity e</td><td>B[1.52, 29]<sup>c</sup></td><td>0.04<sup>+0.045</sup><sub>−0.025</sub></td></tr><tr><td>Angle of periastron ω (deg)</td><td>U[0, 360]</td><td>179<sup>+128</sup><sub>−125</sub></td></tr><tr><td>Doppler semi-amplitude variation K (m s<sup>−1</sup>)</td><td>U[0, 50]</td><td>5.8 ± 2.5</td></tr><tr><td colspan="3">AU Mic c's parameters</td></tr><tr><td>Orbital period P<sub>orb</sub> (d)</td><td>N[18.859019, 0.000016]</td><td>18.859019 ± 0.0000016</td></tr><tr><td>Transit epoch T<sub>0</sub> (BJD<sub>TDB</sub> −2450000)</td><td>N[8342.2223, 0.0005]</td><td>8342.22231 ± 0.00050</td></tr><tr><td>Eccentricity e</td><td>B[1.52, 29]<sup>c</sup></td><td>0.041<sup>+0.047</sup><sub>−0.026</sub></td></tr><tr><td>Angle of periastron ω (deg)</td><td>U[0, 360]</td><td>153<sup>+124</sup><sub>−94</sub></td></tr><tr><td>Doppler semi-amplitude variation K (m s<sup>−1</sup>)</td><td>U[0, 50]</td><td>8.5 ± 2.5</td></tr><tr><td colspan="3">Multidimensional GP parameters</td></tr><tr><td>GP Period P<sub>GP</sub> (d)</td><td>U[4.5, 5.5]</td><td>4.8571<sup>+0.0037</sup><sub>−0.0027</sub></td></tr><tr><td>λ<sub>p</sub></td><td>U[0.1, 2]</td><td>0.449<sup>+0.049</sup><sub>−0.043</sub></td></tr><tr><td>λ<sub>e</sub> (d)</td><td>U[10, 500]</td><td>108 ± 15</td></tr><tr><td>A<sub>v</sub> (m s<sup>−1</sup>)</td><td>U[0, 100]</td><td>7.5<sup>+10.9</sup><sub>−05.5</sub></td></tr><tr><td>B<sub>v</sub> (m s<sup>−1</sup> d)</td><td>U[−1000, 1000]</td><td>91<sup>+27</sup><sub>−18</sub></td></tr><tr><td>A<sub>α</sub> (1000 m<sup>2</sup> s<sup>−2</sup>)</td><td>U[0, 1000]</td><td>38.7<sup>+12.3</sup><sub>−08.5</sub></td></tr><tr><td colspan="3">Other parameters</td></tr><tr><td>C<sub>v</sub> (Offset RV, km s<sup>−1</sup>)</td><td>U[−0.82, 0.71]</td><td>−0.0982<sup>+0.0047</sup><sub>−0.0043</sub></td></tr><tr><td>C<sub>α</sub> (Offset DLW, m<sup>2</sup> s<sup>−2</sup>)</td><td>U[−1, 1]</td><td>−0.003<sup>+0.018</sup><sub>−0.019</sub></td></tr><tr><td>Jitter term σ<sub>j,RV</sub> (m s<sup>−1</sup>)</td><td>J[0, 1000]</td><td>10.3<sup>+1.6</sup><sub>−1.4</sub></td></tr><tr><td>Jitter term σ<sub>j,DLW</sub> (1000 m<sup>2</sup> s<sup>−2</sup>)</td><td>J[0, 1000]</td><td>30.7<sup>+2.7</sup><sub>−2.4</sub></td></tr><tr><td colspan="3">Derived Parameters</td></tr><tr><td colspan="3">AU Mic b's derived Parameters</td></tr><tr><td>Planet mass M<sub>p</sub> (M<sub>⊕</sub>)</td><td>–</td><td>11.7 ± 5.0</td></tr><tr><td>Planet density ρ<sub>p</sub> (g cm<sup>−3</sup>)</td><td>–</td><td>0.97 ± 0.43</td></tr><tr><td colspan="3">AU Mic c's derived Parameters</td></tr><tr><td>Planet mass M<sub>p</sub> (M<sub>⊕</sub>)</td><td>–</td><td>22.2 ± 6.7</td></tr><tr><td>Planet density ρ<sub>p</sub> (g cm<sup>−3</sup>)</td><td>–</td><td>3.66 ± 1.28</td></tr></tbody></table>

Note. $^a$U[a, b] refers to uniform priors between a and b, $N[μ, σ]$ to Gaussian priors with mean μ and standard deviation σ, $B[a, b]$ to a beta distribution with shape parameters a and b, and $J[a, b]$ is the modified Jeffrey's prior as defined by Gregory (2005, equation 16). $^b$Inferred parameters and errors are defined as the median and 68.3 per cent credible interval of the posterior distribution. $^c$Beta distribution to inform eccentricity sampling using the beta distribution for multi planetary systems as defined by Van Eylen et al. (2019).

5000 steps, thinning the chains by a factor of 10 so that samples are uncorrelated, to create the final posterior distribution (corresponding to 125 000 independent samples for each parameter). The resulting 'corner' plot, showing 1D and 2D posterior densities for all the parameters, is shown in the Appendix (Fig. A1). Except for the longitude of periastron ω, which is unconstrained for both planets, the marginal posterior distributions for each of the sampled parameters are unimodal, indicating that there are no pathological degeneracies between the parameters.

## 3.5 Results

We report fitted values and upper and lower uncertainties for each parameter in the final column of Table 3, using the median and 16th and 84th percentiles of the marginal posterior distributions. The inferred model is shown alongside the data in Fig. 4. Fig. 5 shows the phase-folded individual RV signal for each planet. The signals of AU Mic b and AU Mic c are detected at the 2.3σ and 3.4σ level, respectively, with amplitudes of 5.8 ± 2.5 and 8.5 ± 2.5 m s⁻¹. While AU Mic c is detected at high confidence, the detection of AU Mic b is more tentative in this data set. The eccentricities are consistent with zero.

The GP hyper-parameters are well constrained, and consistent with prior information about the rotation and activity behaviour of AU Mic. The GP period is $P_{GP} = 4.8571^{+0.004}_{-0.003}$ d, consistent with the rotation period $P_{rot} = 4.863 \pm 0.010$ d reported by P20. The GP evolution time-scale $\lambda_e = 108 \pm 15$ d is 20 times the rotation period and is consistent with estimates from previous studies where a GP was used to describe the RVs only (P20, K21, C21). The inverse harmonic complexity $\lambda_p = 0.449^{+0.049}_{-0.043}$ is relatively low. This can arise for a number of reasons: if the distribution of active regions on the stellar surface was relatively complex at the time of the observations, or if the latent variable $G(t)$ displays beat patterns caused by either differential rotation or the evolution of individual active regions (we note K21 detected solar-like differential rotation using Doppler Imaging with SPIRou).

We now examine the behaviour of the parameters $A_v$, $B_v$, and $A_\omega$, which control the relationship between the observables and the underlying latent variable $G(t)$. As shown in Table 3 and Fig. A1, the posterior for $A_v$ is consistent with zero at the 1.5σ level, while those for $B_v$ and $A_\alpha$ are significantly non-zero. In other words, the RV time-series behaves primarily like $\dot{G}(t)$, while the DLW time-series behaves primarily like $G(t)$. This is consistent with the prediction
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Radial velocity and DLW time-series (red symbols with error bars) with the best-fitting model (solid lines) and 3σ confidence intervals (shaded grey areas). Note the break in the x-axis corresponding to a ~5-month break in the observations. The top panel shows the RV data with the full model in black. The activity and Keplerian components are also shown separately, in blue and cyan respectively, vertically offset for clarity. The second panel shows the data after subtracting the activity model, with the Keplerian component overplotted. The final RV residuals are shown in the third panel. The fourth panel shows the DLW time-series, with the activity model in black, and the DLW residuals are shown in the fifth panel. In all the panels, the nominal error bars are in solid colour, and the error bars taking into account the jitter are semitransparent.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Phase-folded RV signals for AU Mic b (left) and AU Mic c (right) after subtracting the systemic velocity, the activity signal, and the other planet. The grey circles with error bars show the HARPS measurements, and the error bars taking into account the jitter are shown with lighter grey. The black line show the best-fitting Keplerian model. The red circles show binned RVs to facilitate comparison between data and models.</figcaption></figure>

of Aigrain et al. (2012) and with what has been observed for other young, rapidly rotating, active stars including K2-100 (Barragán et al. 2019b) and HD 73583 (Barragán et al. 2022a).

The jitter terms, $\sigma_{RV}$ and $\sigma_{DLW}$, can be used as diagnostics of the effectiveness (or otherwise) of our model at explaining the full data set. In both cases, the posteriors medians are ~10 times smaller than the amplitude of the dominant activity term ($B_V$ and $A_\alpha$, respectively), but roughly twice the nominal measurement uncertainties. This illustrates the fact that our activity model explains most of the RV and DLW variations, but is not perfect.

A widely used test of the quality of a model fit is to compare the root-mean-square (RMS) before and after subtracting the best-fitting model. Caution must be used in interpreting such a comparison, as a low residual RMS can result from overfitting as readily as from a particularly good model. On the other hand, an anomalously high residual RMS would be a robust indication that the model cannot fully explain the observations. The initial RMS of the RV observations is $137 \text{ m s}^{-1}$, compared to $7.5 \text{ m s}^{-1}$ after subtracting our best fit, combined activity, and planetary signal model. The DLW time-series has an initial RMS of $45 \times 10^3 \text{ m}^2 \text{ s}^{-2}$ which shrinks to $30 \times 10^3 \text{ m}^2 \text{ s}^{-2}$ after subtracting the best-fitting model (which, in the case of DLW, includes activity only). Thus, the RMS of both time-series is reduced by 94 per cent for the RVs and 33 per cent for the DLW, indicating that our model explains most of the variability seen in the RVs and a significant fraction of that seen in the DLWs. The residual RMS is none the less significant compared to the median formal uncertainty of each time-series ($3.6 \text{ m s}^{-1}$ and $4.8 \text{ m s}^{-2}$, respectively). This highlights that our model, while useful, is incomplete. The remaining variability might be due to stellar signals which our model cannot account for, or to instrumental systematics. Either way, it is absorbed by the jitter terms which, as

we can see from Table 3, have values similar to the residual RMS (as one would expect).

We also test if the planetary signals are recovered individually in the data. We ran three different set-ups: one including only planet b, one only planet c, and the last one including both signals. The results of this test are summarized in Appendix C. We found that we are able to recover each planetary signal individually, even if the signal of the other planet is not included in the model. In addition, we found that the model producing the highest likelihood is the one including both planetary signals. However, we note that if we use standard model comparison parameters such as Akaike Information Criterion and Bayesian Information Criterion, the preferred model is the one including AU Mic c signal only. This is expected given that the detection of AU Mic b is not well constrained. None the less, we adopt the two-planet model given that this is physically motivated. We know a priori that both planets are present from the transit observations.

# 4 DISCUSSION

In this section, we discuss the implications of the results presented in Section 3.5 for the masses, orbits, composition, and evolution of AU Mic b and c.

## 4.1 Planet masses

From our estimates for $K_b = 5.8 \pm 2.5 \text{ m s}^{-1}$ and $K_c = 8.5 \pm 2.5 \text{ m s}^{-1}$, we directly obtain mass estimates for the two planets of $M_b = 11.7 \pm 5.0 M_\oplus$ and $M_c = 22.2 \pm 6.7 M_\oplus$. This represents a $3.4\sigma$ detection of planet c, but only a tentative, $2.3\sigma$ detection of planet b.
K21 reported a mass estimate of $17.1 \pm 4.6 M_{\oplus}$ for planet b, which is consistent at the $1\sigma$ level with our estimate. We note that K21 did not include AU Mic c in their analysis, since the second planet was not confirmed at the time.

More recently, C21 reported a mass of $20.1 \pm 1.6 M_{\oplus}$ for planet b and a $5\sigma$ upper limit for planet c of $M_c < 20.1 M_{\oplus}$, using a combination of data from multiple optical and near-IR spectrographs, and a model which includes only the RVs, but exploits the chromaticity of the activity signal. Our results are not consistent with these values. Although it is difficult to identify the exact source of the discrepancy, several factors are likely at play. The first is the difference in the wavelength coverage; including both optical and near-IR data should in principle give a better handle on the activity signals, though it also necessitates an increase in the number of free parameters in the model. Another difference is the time-sampling of the data sets used, which is much sparser in the case of C21 than in this work. Sparser time sampling severely limits the ability to model the activity signal. The third is the fact that we model an activity indicator simultaneously with the RVs, which again should help constrain the activity component of the model. We note that C21 performed two different analyses, with different assumptions about the chromaticity of the activity signal, and that these gave mutually inconsistent results, particularly for the mass of planet c.

Another difference is that C21 use a Gaussian prior with mean 0.19 and standard deviation 0.04 for the eccentricity of planet b, derived from *Spitzer* secondary eclipse observations. As the latter are not yet published, we did not incorporate this prior in the analysis presented in Section 3.5, but we did carry out an additional model run with the same eccentricity prior for planet b as used by C21. This did not appreciably affect our mass estimates for either planet.

Finally, exploring the parameter space for such complex, activity plus Keplerian models, is in general challenging. To check how robust our results are to the details of the model, we carried out simple tests including only planet b or only planet c, and found that the resulting RV semi-amplitudes were consistent with the values obtained from the model including both planets. We also carried out a number of additional tests using simulations, described in Section 4.4, to further establish confidence in our results. Overall, however, it is clear that continued RV monitoring of the system, as well as independent analyses of all the available data sets by different teams, would be desirable to resolve the discrepancies between the published mass estimates.

## 4.2 TTV analysis

As AU Mic b and c are close to a 2:1 mean-motion resonance, transit timing variations (TTVs) can be used to constrain their masses and eccentricities. Szabó et al. (2021) and Gilbert et al. (2021) report tentative TTVs for AU Mic b at the level of 3–4 min and 80 s, respectively, using *TESS* and *CHEOPS* observations. Although the measurement of transit times is somewhat difficult because of AU Mic’s frequent flaring, these results indicate that strong TTVs, with amplitudes in excess of ~10 min, are not present in the system. This can be used to provide a sanity check on the mass and eccentricity estimates derived in this work and in previous RV analyses, and to motivate future TTV follow-up.

Using the *TTVFast* package (Deck et al. 2014) we forward modelled the TTVs expected from a two-planet configuration as observed in the AU Mic system. TTV amplitudes are most sensitive to planetary mass and eccentricity, so we varied those on a grid over a plausible range for these parameters, noting the amplitude of the TTV signal in each case. Our predictions are summarized in Appendix D. Even for the most extreme mass ratio and eccentricities compatible

with the results of our RV analysis, the maximum amplitude of the TTVs never exceeds a few minutes, with a superperiod of 80 d. We therefore conclude that our mass estimates are consistent with existing transit observations.

## 4.3 Dynamical analysis

We used *mercury6* (Chambers 1999) to test the dynamical stability of the AU Mic system. We used the masses and orbital parameters reported in Table 3, assuming that both planets have coplanar orbits. Starting from this configuration, we evolved the system for 1 Gyr with steps of 0.5 d per integration. We found that the eccentricity and semimajor axis of both planets display periodic fluctuations with a period of ~1000 yr, but these fluctuations are contained and the system appears stable over the full duration of the simulation. The eccentricity of AU Mic b oscillates between 0.05 and 0.25, and its semimajor axis varies by $8 \times 10^{-5}$ au. The eccentricity of AU Mic c remains <0.10 and the maximum departure of its semimajor axis from the starting value is $1.3 \times 10^{-4}$ au. We conclude that the orbital parameters and masses derived for AU Mic planets are consistent with a dynamically stable system.

## 4.4 Tests to establish detection robustness

It has been shown that the combination of complex models and the window function of the observations can create spurious planet-like signals in RV time-series, specially for active stars (e.g. Rajpaul, Aigrain & Roberts 2016). In order to check the reliability of our RV detection, we performed numerical simulations similar to Barragán et al. (2019b). We used *citlalatonac* (B22) to simulate an RV and DLW time-series containing only the best-fitting activity model as obtained from the real data. We decided to use the median values of our derived parameters given that the posteriors are quasi-Gaussian unimodal distributions. We then added correlated noise using a squared exponential kernel with a length-scale of 1 d, and the same amplitude as the jitter term obtained from the real data, plus white noise for each observation according to its nominal measurement uncertainty. We did this 100 times to obtain 100 simulated activity-only time-series, with similar noise properties and the same time-sampling as the real data.

We modelled each activity-only synthetic data set using a two-planet and two-dimensional GP configuration as described in Section 3.4. For each simulation, we plot the posterior over the semi-amplitudes for both ‘planets’ in the left column of Fig. 6, compared to the posterior obtained from the real data set. We then count the fraction of the simulations where the recovered semi-amplitude for planets b and c equals or exceeds the median value obtained from the actual observations. For planet b, this occurred 11 per cent if the time in the 100 simulations. For planet c, it occurred in 2 per cent of the cases. Taking Poisson counting errors into account, this is fully consistent with the confidence intervals that we derived for the real observations.

We then repeat the same exercise, creating another 100 synthetic data sets, but this time, injecting two Keplerian signals with the median planet parameters obtained from the real data set, as reported in Table 3. Again, we modelled these synthetic data sets using *pyaneti* https://github.com/oscaribv/pyaneti with the same configuration as described in Section 3.4. The resulting posteriors are shown in the right column of Fig. 6. In this case, we are interested in two distinct questions. First, how often are AU Mic b and AU Mic c detected at a given confidence level from these simulations? For this purpose, we consider that a detection has occurred if the median of
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Results of the activity-only (left) and activity plus Keplerian (right) simulations (see Section 4.4 for details). In each panel, the red histogram shows the posterior derived from the original data, and the thin grey lines show the posteriors obtained from the simulations.</figcaption></figure>

the posterior is larger than $2\sigma$, where $\sigma$ is half the interval between the 16.5th and 83.5th percentiles. We find that AU Mic b is detected 42 per cent of the time, and AU Mic c 88 per cent of the time. This confirms that, if the planet masses are similar to the maximum a posteriori values we derived from the real data, more data are needed for a robust detection of planet b. Second, we asked what fraction of the time the recovered semi-amplitude is within $2\sigma$ of the injected value? We find that this happens 91 per cent and 89 per cent of the time, for AU Mic b and AU Mic c, respectively. For Gaussian posteriors, we would expect these numbers to be around 95 per cent. Again, accounting for Poisson counting errors, this result is thus consistent with what we expect.

## 4.5 Composition and evolution

Having measured radii and masses for both planets allows to place them in a mass-radius (M-R) diagram, as shown in Fig. 7. Also shown on that figure are other planets with well-determined masses and radii and theoretical mass-radius relations for terrestrial and ocean worlds (Zeng et al. 2016), and for *evolved* planets with H/He envelopes (Zeng et al. 2019) accounting for 1–5 per cent of the total mass, at a temperature of 500 K, which is close to the equilibrium temperatures of both AU Mic b and c (593 ± 21 K and 454 ± 16 K, respectively, see Martioli et al. 2020 and Table 1). We note that although it is not ideal to directly compare these young planets with their older counterparts, this is a standard practice in the field.

Comparison to these theoretical composition curves for evolved, ~Gyr old planets indicates that planet c is compatible with a pure H$_2$O world, and planet b requires an atmosphere containing some H/He. However, it is difficult to imagine how a pure H$_2$O planet could form. Coupled with the fact planet b hosts an atmosphere containing H/He, we consider it more likely that both planets

possess voluminous H/He atmospheres. However, comparison to these evolved mass-radius relationships is indicative only. Young planets still contain considerable thermal energy left over from their assembly (e.g. Ginzburg, Schlichting & Sari 2016), resulting in significantly larger planetary radii compared to evolved planets with the same composition (e.g. Owen 2020).

Thus, we compare the measured masses and radii to evolutionary models for close-in exoplanets that include the impact of cooling and contraction of any H/He atmosphere and photoevaporative loss using MESA models (Owen & Wu 2013). The evolutionary calculations are identical to those used in Owen (2020) and Mann et al. (2021). In this manner, the measured masses and radii of planets b and c can be used to constrain their possible evolutionary pathways, and hence their composition and thermodynamic state both today and at their formation.

We find that both planets formed with voluminous, but low-mass H/He dominated atmospheres (~8 per cent and ~3 per cent by mass for planets b and c, respectively, as shown in the left-hand panel of Fig. 8), and as such are likely progenitors of the ubiquitous super-Earth/sub-Neptune population. This is further evidenced when we evolve these planets forward to an age of 5 Gyr finding that in the vast majority of cases these planets retain a fraction of their H/He atmosphere, reaching a radius of ~2.9 and ~2.7 $R_⊕$ for b and c, respectively (see Fig. 8, right-hand panel). However, there are two standout results. The first is that the amount of H/He these planets accreted from their parent discs is significantly less than predicted from standard models of core-accretion (e.g. Lee & Chiang 2015). Based on those models, we would have expected planet b, for example, to have accreted between 25 and 75 per cent of its mass in H/He. This fraction is high enough that it could plausibly have undergone runaway accretion. The envelope mass fraction we infer from the measured mass and radius of planet b
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Mass versus radius diagram for small exoplanets ($1 < R_p < 5R_⊕$ and $1 < M_p < 32M_⊕$). The grey points with error bars show planets with mass and radius measurements better than 50 per cent (As in the TEPCAT catalogue, https://www.astro.keele.ac.uk/jkt/tepcat/, Southworth, Wheatley & Sams 2007). AU Mic b and c are shown with a green and orange circle, respectively. The green and orange arrows indicate the expected evolutionary pathway for the planets and their final state after 5 Gyr. The green and orange squares show the expected position of AU Mic b and c after 5 Gyr, respectively. The solid lines represent two-layer models as given by Zeng, Sasselov & Jacobsen (2016) with a different colour corresponding to a different mixture of elements. The purple non-solid lines correspond to Earth-like rocky cores surrounded by a Hydrogen envelope of 1 per cent (dash-dotted line), 2 per cent (dotted line), and 5 per cent (dashed) Hydrogen mass for evolved exoplanets assuming an equilibrium temperature of 500 K (Zeng et al. 2019). **Note:** One should compare the expected evolved position of the planets with the composition curves, rather than their current position.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Initial envelope mass fractions and core masses (left – with contours shown at $0.5\sigma$ intervals) and radii after 5 Gyr (right) for AU Mic b and c (in green and orange, respectively) compatible with the measured masses and radii, based on evolutionary calculations including cooling and contraction and photoevaporation (see Section 4.5 for details).</figcaption></figure>

is 3–10 times smaller than this theoretical prediction. The second standout result is the trend between the two planets, where planet c has a higher inferred core-mass, but accreted less H/He than planet b. This is exactly the opposite of what is expected from core-accretion,

where more massive cores accrete more H/He (e.g. Pollack et al. 1996).

There have already been indications that close-in planets may accrete less H/He than predicted by core-accretion, based on a sample
of old individual planets (Jankovic, Owen & Mohanty 2019), and from inferences about the population of planets discovered by *Kepler* (Rogers & Owen 2021). However, this is the first time it has been directly identified in young planets that have not been significantly affected by atmospheric mass-loss. Possible mechanisms to resolve this tension with core-accretion include forming the planets late in the disc's lifetime, during the dispersal phase (e.g. Ikoma & Hori 2012; Lee & Chiang 2016), or the conclusion from 3D simulations that high-entropy gas is continually recycled into the forming atmosphere, preventing it from cooling and accreting more gas (e.g. Ormel, Shi & Kuiper 2015; Ali-Dib, Cumming & Lin 2020; Chen et al. 2020). Finally, additional mass-loss during disc dispersal can dramatically remove accreted material, as the disc's pressure confinement of the proto-atmosphere is rapidly removed (e.g. Ginzburg et al. 2016; Owen & Wu 2016).

However, none of these proposed mechanisms can reconcile core-accretion with the observation that planet c, with its more massive core, started off with a smaller H/He atmosphere. This is even more curious as planet c appears massive enough to have undergone runaway accretion, even at its current short-period orbit, and formed a giant planet (e.g. Rafikov 2006; Lee, Chiang & Ormel 2014). However, it only had an initial H/He mass fraction of a few per cent.

One possible explanation is that the two planets formed further from the star than their present-day locations, then migrated inwards. Standard core accretion models predict that the mass of the accreted atmosphere grows with core mass as ~ $M_{\text{core}}^{1.7}$ (as evidenced by Lee & Chiang 2015), and is approximately logarithmically sensitive to the disc pressure (e.g. Piso & Youdin 2014). Applying these scaling laws directly to the core and envelope masses we infer for AU Mic b and c, it would imply that the two planets formed in regions of the disc with extremely different densities (by a factor $\gtrsim \exp(10)$). Such a scenario would also require that planet c did not start migrating inwards until the disc dispersed (otherwise it would have accreted more gas). It is then difficult to see how the planet would have reached its current orbit by the present age of the system, in the absence of a gas disc. On the other hand, the critical core mass for runaway gas accretion depends on the rate at which the core is accreting solid material, which itself depends on location within the disc (Terquem 2014). This opens up the possibility of a scenario where the planets form at different locations in the disc, but less extremely so, and migrate inwards during the late stages of the disc's dispersal. Exploring this scenario further using detailed calculations is, however, beyond the scope of this work.

Another possible explanation would be that planet c started forming later, and thus acquired its envelope over a shorter time, than planet b. However, the accreted mass has a very weak dependence on accretion time ($M_{\text{env}} \propto t^{0.4}$, Lee & Chiang 2015). Following the same logic as the gas surface density arguments given above, this would require planet c to form over a time period $10^{10.4} \sim 300$ times shorter than planet b. Again, this would require a great deal of fine tuning.

A second, speculative scenario that would reconcile both c's large core-mass and lower H/He atmosphere mass compared to b is that c is the product of a post-disc-dispersal giant impact. Giant impacts are not an unlikely outcome for planets that find themselves on short-period orbits after disc-dispersal (e.g. Bonomo et al. 2019; Poon et al. 2020). Giant impacts can remove large amounts of H/He in any collision (e.g. Liu et al. 2015; Inamdar & Schlichting 2016) in addition to the fact that the collision inflates any residual H/He atmosphere resulting in enhanced photoevaporation (Biersteker & Schlichting 2019). Planets in multiple systems found by *Kepler* tend to have similar radii (e.g. Weiss et al. 2018) and masses (e.g.

Millholland, Wang & Laughlin 2017) to each other. This has been interpreted as evidence for uniformity in the underlying core-mass within each system (Millholland & Winn 2021). Thus, if AU Mic formed three planets with properties similar to b, and two underwent a giant impact, a planet with properties similar to c would be the natural outcome. One way to test this speculative scenario is to measure the obliquity of planet c. Any indication of difference between the spin-orbit angles of the two planets (the orbital plane of planet b having been found aligned with the stellar rotation axis; Hirano et al. 2020; Martioli et al. 2020; Addison et al. 2021; Szabó et al. 2021) would indicate that significant dynamical interaction between planets in the system took place after disc-dispersal. The amplitude of the Rossiter-McLaughlin (RM) effect induced by planet c is expected to be of order ~$10 \text{ m s}^{-1}$, which is feasible with current high-precision spectrographs.

## 5 CONCLUSIONS

We have analysed data from an intensive monitoring campaign using HARPS with multidimensional GP framework (B22), and measured the masses of AU Mic b and c to be $11.7 \pm 5.0 \ M_{\oplus}$ and $22.2 \pm 6.7 \ M_{\oplus}$. We detect planet c at $3.4\sigma$ confidence and planet b at $2.3\sigma$ confidence. Our results are insensitive to the choice of RV time-series (DRS, SERVAL, HARPS-TERRA) and to the specific activity indicators used, though this choice affects the degree of confidence in the detections slightly. Our mass measurements, combined with literature radius estimates, indicate that AU Mic b has a significant H/He envelope and may have an internal heat source. On the other hand, we find that AU Mic c is denser and compatible with either a rocky core surrounded by a ~2 per cent H/He envelope or a pure $H_2O$ composition. These results are in tension with current core-accretion models, which predict that planet c, with a more massive core, should have the larger H/He envelope of the two. We provide a possible, speculative explanation for this discrepancy involving a giant impact, but stress that additional observations are needed to further refine the mass measurements and constrain the obliquity of planet c.

The RV signals of the two planets around AU Mic are around 50 times smaller than the activity signals. Disentangling the former from the latter is very challenging, as highlighted by the discrepancies between our results and those of other teams working on the same system using different instruments and methods, which need to be better understood. Further observations of this system and independent analyses of the available data sets are needed to resolve these discrepancies. However, the fact that these detections were possible at all is very encouraging for the continued characterization of the increasing number of young transiting planets being discovered by surveys such as *K2* and *TESS*. The characterization of planets around young stars also provides an extreme test case for activity mitigation techniques in RVs. The methods being developed in this context will prove useful for future searches for smaller and cooler planets around less active stars, where the amplitude ratios between planetary and activity signals are similar.

## ACKNOWLEDGEMENTS

We gratefully acknowledge X. Dumusque and F. Bouchy for their coordination of the HARPS time-share, which allowed us to spread the observations over the two semesters with the appropriate time-sampling, and all the observers involved in the time-share (F. Bouchy, B. Canto, I. de Castro, G. Dransfield, M. Esposito, V. Van Eylen, F. Hawthorn, M. Hobson, V. Hodzic, D. Martin, J. McCormac,
H. Osborne, J. Otegi, A. Suarez, P. Torres) for carrying out the observations on our behalf. We also thank J. Patterson for managing the Oxford Astrophysics compute cluster, *glamdring*, which was used to carry out the data analysis and simulations. In addition, we thank Caroline Terquem for her insights and suggestions regarding the formation of the planets. This study is based on observations collected at the European Southern Observatory under ESO programmes 0105.C-02887 & 0106.C-0852 (PIs Aigrain and Zicher). This research has used NASA's Astrophysics Data System. NZ acknowledges support from the UK Science and Technology Facilities Council (STFC) under Grant Code ST/N504233/1, studentship no. 1947725. This publication is part of a project that has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 865624). LDN thanks the Swiss National Science Foundation for support under Early Postdoc.Mobility grant P2GEP2_200044. JEO is supported by a Royal Society University Research Fellowship. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 853022, PEVAP). Some of this work was performed using the DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations grant ST/R001014/1. DiRAC is part of the National e-Infrastructure. NZ, OB, and BK would like to thank the pizzeria and café in Oxford that gave us the energy to detect the two planetary signals during COVID-19 times.

## DATA AVAILABILITY

The HARPS spectra used in this study are available on the ESO archive.<fnref n="4" /> The radial velocities and activity indicators extracted with the DRS, HARPS-TERRA, and SERVAL are listed in Table 2, a machine-readable version of which is provided in the Supplementary Material. The code used to perform the evolutionary analysis is available at MESAPlanet https://github.com/jo276/MESAlanet.<fnref n="5" /> Pyaneti is available at pyaneti https://github.com/oscaribv/pyaneti.<fnref n="6" />

## REFERENCES

* Addison B. C. et al., 2021, AJ, 162, 137

* Afram N., Berdyugina S. V., 2019, A&A, 629, A83

* Aigrain S., Pont F., Zucker S., 2012, MNRAS, 419, 3147

* Ali-Dib M., Cumming A., Lin D. N. C., 2020, MNRAS, 494, 2440

* Anglada-Escudé G., Butler R. P., 2012, ApJS, 200, 15

* Astropy Collaboration, 2013, A&A, 558, A33

* Astropy Collaboration, 2018, AJ, 156, 123

* Barragán O. et al. 2022a, MNRAS, preprint (arXiv:2110.13069)

* Barragán O. et al., 2019b, MNRAS, 490, 698

* Barragán O., Aigrain S., Rajpaul V. M., Zicher N., 2022b, MNRAS, 509, 866

* Barragán O., Gandolfi D., Antoniciello G., 2019a, MNRAS, 482, 1017

* Batalha N. E., Lewis T., Fortney J. J., Batalha N. M., Kempton E., Lewis N. K., Line M. R., 2019, ApJ, 885, L25

* Bean J. L., Raymond S. N., Owen J. E., 2021, J. Geophys. Res., 126, e06639

* Beuzit J. L. et al., 2019, A&A, 631, A155

* Biersteker J. B., Schlichting H. E., 2019, MNRAS, 485, 4454

* Boccaletti A. et al., 2015, Nature, 526, 230

* Boccaletti A. et al., 2018, A&A, 614, A52

* Bonomo A. S. et al., 2019, Nat. Astron., 3, 416

* Cale B. L. et al., 2021, AJ, 162, 295

* Chambers J. E., 1999, MNRAS, 304, 793

* Chen Y.-X., Li Y.-P., Li H., Lin D. N. C., 2020, ApJ, 896, 135

* Claret A., 2018, A&A, 618, A20

* David T. J. et al., 2016a, AJ, 151, 112

* David T. J. et al., 2016b, Nature, 534, 658

* David T. J. et al., 2019a, AJ, 158, 79

* David T. J., Petigura E. A., Luger R., Foreman-Mackey D., Livingston J. H., Mamajek E. E., Hillenbrand L. A., 2019b, ApJ, 885, L12

* Deck K. M., Agol E., Holman M. J., Nesvorný D., 2014, ApJ, 787, 132

* Donati J. F. et al., 2016, Nature, 534, 662

* Donati J. F. et al., 2020, MNRAS, 498, 5684

* Gaia Collaboration, 2021, A&A, 649, A1

* Gaidos E. et al., 2017, MNRAS, 464, 850

* Gilbert E. A. et al., 2021, preprint (arXiv:2109.03924)

* Ginzburg S., Schlichting H. E., Sari R., 2016, ApJ, 825, 29

* Gomes da Silva J., Santos N. C., Bonfils X., Delfosse X., Forveille T., Udry S., 2011, A&A, 534, A30

* Grandjean A. et al., 2020, A&A, 633, A44

* Gregory P. C., 2005, ApJ, 631, 1198

* Hebb L., Petro L., Ford H. C., Ardila D. R., Toledo I., Minniti D., Golimowski D. A., Clampin M., 2007, MNRAS, 379, 63

* Hirano T. et al., 2020, ApJ, 899, L13

* Howell S. B. et al., 2014, PASP, 126, 398

* Ikoma M., Hori Y., 2012, ApJ, 753, 66

* Inamdar N. K., Schlichting H. E., 2016, ApJ, 817, L13

* Jankovic M. R., Owen J. E., Mohanty S., 2019, MNRAS, 484, 2296

* Kalas P., Liu M. C., Matthews B. C., 2004, Science, 303, 1990

* Klein B. et al., 2021, MNRAS, 502, 188

* Klein B. et al., 2022, MNRAS, Available at: https://academic.oup.com/mnras/advance-article/doi/10.1093/mnras/stac761/6552137

* Kossakowski D. et al., 2021, A&A, 656, A124

* Lannier J. et al., 2017, A&A, 603, A54

* Lee E. J., Chiang E., 2015, ApJ, 811, 41

* Lee E. J., Chiang E., 2016, ApJ, 817, 90

* Lee E. J., Chiang E., Ormel C. W., 2014, ApJ, 797, 95

* Liu S.-F., Hori Y., Lin D. N. C., Asphaug E., 2015, ApJ, 812, 164

* Lo Curto G. et al., 2015, The Messenger, 162, 9

* Lopez E. D., Fortney J. J., 2013, ApJ, 776, 2

* Mamajek E. E., Bell C. P. M., 2014, MNRAS, 445, 2169

* Mann A. W. et al., 2016a, AJ, 152, 61

* Mann A. W. et al., 2016b, ApJ, 818, 46

* Mann A. W. et al., 2017, AJ, 153, 64

* Mann A. W. et al., 2018, AJ, 155, 4

* Mann A. W. et al., 2020, AJ, 160, 179

* Mann A. W. et al., 2021, preprint (arXiv:2110.09531)

* Martioli E. et al., 2020, A&A, 641, L1

* Martioli E., Hébrard G., Correia A. C. M., Laskar J., Lecavelier des Etangs A., 2021, A&A, 649, A177

* Mayor M. et al., 2003, The Messenger, 114, 20

* Millholland S. C., Winn J. N., 2021, ApJ, 920, L34

* Millholland S., Wang S., Laughlin G., 2017, ApJ, 849, L33

* Newton E. R. et al., 2019, ApJ, 880, L17

* Newton E. R. et al., 2021, AJ, 161, 65

* Ormel C. W., Shi J.-M., Kuiper R., 2015, MNRAS, 447, 3512

* Owen J. E., 2019, Annu. Rev. Earth Planet. Sci., 47, 67

* Owen J. E., 2020, MNRAS, 498, 5030

* Owen J. E., Wu Y., 2013, ApJ, 775, 105

<footnote n="7">⁷The P105 observations were delayed because of the COVID-19 pandemic and taken in P107.</footnote>

<footnote n="4">⁴http://archive.eso.org/eso/eso_archive_main.html</footnote>

<footnote n="5">⁵https://github.com/jo276/MESAlanet</footnote>

<footnote n="6">⁶https://github.com/oscaribv/pyaneti</footnote>
Owen J. E., Wu Y., 2016, ApJ, 817, 107

Pepe F. A. et al., 2010, in Ground-based and Airborne Instrumentation for Astronomy III. Proceedings of SPIE - The International Society for Optical Engineering, p. 77350F

Pepper J. et al., 2017, AJ, 153, 177

Piso A.-M. A., Youdin A. N., 2014, ApJ, 786, 21

Plavchan P. et al., 2020, Nature, 582, 497

Plavchan P., Werner M. W., Chen C. H., Stapelfeldt K. R., Su K. Y. L., Stauffer J. R., Song I., 2009, ApJ, 698, 1068

Pollack J. B., Hubickyj O., Bodenheimer P., Lissauer J. J., Podolak M., Greenzweig Y., 1996, Icarus, 124, 62

Poon S. T. S., Nelson R. P., Jacobson S. A., Morbidelli A., 2020, MNRAS, 491, 5595

Quinn S. N. et al., 2012, ApJ, 756, L33

Quirrenbach A. et al., 2014, in Ramsay S. K., McLean I. S., Takami H., eds, Proc.SPIE Conf. Ser. Vol. 9147, Ground-based and Airborne Instrumentation for Astronomy V. SPIE, Bellingham, p. 91471F

Rafikov R. R., 2006, ApJ, 648, 666

Rajpaul V., Aigrain S., Osborne M. A., Reece S., Roberts S., 2015, MNRAS, 452, 2269

Rajpaul V., Aigrain S., Roberts S., 2016, MNRAS, 456, L6

Ricker G. R. et al., 2015, J. Astron. Telesc. Instrum. Syst., 1, 014003

Rizzuto A. C. et al., 2020, AJ, 160, 33

Rizzuto A. C., Vanderbilt A., Mann A. W., Kraus A. L., Dressing C. D., Agüeros M. A., Douglas S. T., Krolikowski D. M., 2018, AJ, 156, 195

Rogers J. G., Owen J. E., 2021, MNRAS, 503, 1526

Southworth J., Wheatley P. J., Sams G., 2007, MNRAS, 379, L11

Suárez Mascareño A. et al., 2021, Nat. Astron., 6, 232

Szabó G. M. et al., 2021, A&A, 654, A159

Terquem C., 2014, MNRAS, 444, 1738

Torres C. A. O., Quast G. R., da Silva L., de La Reza R., Melo C. H. F., Sterzik M., 2006, A&A, 460, 695

Van Eylen V. et al., 2019, AJ, 157, 61

Vanderburg A. et al., 2018, AJ, 156, 46

Weiss L. M. et al., 2018, AJ, 155, 48

White R., Schaefer G., Boyajian T., von Braun K., ten Brummelaar T., Bieryla A., Dupuy T., Latham D. W., 2019, in American Astronomical Society Meeting Abstracts #233. p. 259.41

Yu L. et al., 2017, MNRAS, 467, 1342

Zechmeister M. et al., 2018, A&A, 609, A12

Zeng L. et al., 2019, Proc. Natl. Acad. Sci., 116, 9723

Zeng L., Sasselov D. D., Jacobsen S. B., 2016, ApJ, 819, 127

## SUPPORTING INFORMATION

Supplementary data are available at MNRAS online.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

## APPENDIX A: CORRELATIONS PLOT

Fig. A1 shows the posterior and correlations plots for the sampled parameters of the described in Section 3. The parameters correspond to the SERVAL RVs and DLW model (Model 3 in Section 3.3).
<figure><img src="image_12.png" /><figcaption>**Figure A1.** Posterior and correlation plots for the sampled parameters (MCMC corner plot). The inferred parameters from these posteriors are given in Table 3.</figcaption></figure>

## APPENDIX B: ACTIVITY MODELS

In Table B1, we show the priors and inferred parameters of the three models described in Section 3.3.
<table><caption>**Table B1.** Priors and posterior median and confidence intervals for the different models.</caption><thead><tr><th>Parameter</th><th>Prior<sup>a</sup></th><th>Model 1<br>Posterior value<sup>b</sup></th><th>Model 2<br>Posterior value<sup>b</sup></th><th>Model 3<br>Posterior value<sup>b</sup></th></tr></thead><tbody><tr><td colspan="5"><strong>AU Mic b's parameters</strong></td></tr><tr><td>Orbital period P<sub>orb</sub> (d)</td><td>N[8.463000, 0.000002]</td><td>8.4629998<sup>+2.1e-06</sup><sub>-2.06</sub></td><td>8.463<sup>+2e-06</sup><sub>-2.1e-06</sub></td><td>8.463000 ± 0.000002</td></tr><tr><td>Transit epoch T<sub>0</sub> (BJD<sub>TDB</sub> −2450000)</td><td>N[8330.39051, 0.00015]</td><td>8330.39052<sup>+0.00015</sup><sub>-0.00016</sub></td><td>8330.39051<sup>+0.00016</sup><sub>-0.00015</sub></td><td>8330.39051 ± 0.00015</td></tr><tr><td>Eccentricity e</td><td>B[1.52, 29]<sup>c</sup></td><td>0.042<sup>+0.05</sup><sub>-0.027</sub></td><td>0.042<sup>+0.044</sup><sub>-0.026</sub></td><td>0.04<sup>+0.045</sup><sub>-0.025</sub></td></tr><tr><td>Angle of periastron ω (deg)</td><td>U[0, 360]</td><td>182.0<sup>+126.0</sup><sub>-126.0</sub></td><td>182.0<sup>+127.0</sup><sub>-129.0</sub></td><td>179<sup>+128</sup><sub>-125</sub></td></tr><tr><td>Doppler semi-amplitude variation K (m s<sup>-1</sup>)</td><td>U[0, 50]</td><td>5.89<sup>+2.98</sup><sub>-2.66</sub></td><td>5.1<sup>+2.32</sup><sub>-2.15</sub></td><td>5.8 ± 2.5</td></tr><tr><td colspan="5"><strong>AU Mic c's parameters</strong></td></tr><tr><td>Orbital period P<sub>orb</sub> (d)</td><td>N[18.859019, 0.000016]</td><td>18.859019<sup>+1.5e-05</sup><sub>-1.6e-05</sub></td><td>18.85902<sup>+1.6e-05</sup><sub>-1.5e-05</sub></td><td>18.859019 ± 0.0000016</td></tr><tr><td>Transit epoch T<sub>0</sub> (BJD<sub>TDB</sub> −2450000)</td><td>N[8342.2223, 0.0005]</td><td>8342.22227<sup>+0.00048</sup><sub>-0.00046</sub></td><td>8342.2223<sup>+0.00051</sup><sub>-0.00051</sub></td><td>8342.22231 ± 0.00050</td></tr><tr><td>Eccentricity e</td><td>B[1.52, 29]<sup>c</sup></td><td>0.039<sup>+0.045</sup><sub>-0.025</sub></td><td>0.04<sup>+0.047</sup><sub>-0.026</sub></td><td>0.041<sup>+0.047</sup><sub>-0.026</sub></td></tr><tr><td>Angle of periastron ω (deg)</td><td>U[0, 360]</td><td>182.0<sup>+126.0</sup><sub>-126.0</sub></td><td>143.3<sup>+141.4</sup><sub>-95.8</sub></td><td>153<sup>+124</sup><sub>-94</sub></td></tr><tr><td>Doppler semi-amplitude variation K (m s<sup>-1</sup>)</td><td>U[0, 50]</td><td>6.7<sup>+3.09</sup><sub>-3.06</sub></td><td>8.86<sup>+2.42</sup><sub>-2.51</sub></td><td>8.5 ± 2.5</td></tr><tr><td colspan="5"><strong>Multidimensional GP parameters</strong></td></tr><tr><td>GP Period P<sub>GP</sub> (d)</td><td>U[4.8, 4.9]</td><td>4.8604<sup>+0.0032</sup><sub>-0.0031</sub></td><td>4.8572<sup>+0.0037</sup><sub>-0.0025</sub></td><td>4.8571<sup>+0.0037</sup><sub>-0.0027</sub></td></tr><tr><td>λ<sub>p</sub></td><td>U[0.1, 2]</td><td>0.455<sup>+0.06</sup><sub>-0.048</sub></td><td>0.46<sup>+0.054</sup><sub>-0.045</sub></td><td>0.449<sup>+0.049</sup><sub>-0.043</sub></td></tr><tr><td>λ<sub>e</sub> (d)</td><td>U[10, 500]</td><td>142.7<sup>+35.1</sup><sub>-24.7</sub></td><td>111.7<sup>+17.3</sup><sub>-16.6</sub></td><td>108 ± 15</td></tr><tr><td>A<sub>V</sub> (m s<sup>-1</sup>)</td><td>U[0, 100]</td><td>21<sup>+24</sup><sub>-14</sub></td><td>10<sup>+15</sup><sub>-07</sub></td><td>7.5<sup>+10.9</sup><sub>-05.5</sub></td></tr><tr><td>B<sub>V</sub> (m s<sup>-1</sup> d)</td><td>U[-1000, 1000]</td><td>102<sup>+53</sup><sub>-26</sub></td><td>102<sup>+39</sup><sub>-23</sub></td><td>91<sup>+27</sup><sub>-18</sub></td></tr><tr><td>A<sub>DLW</sub> (1000 m<sup>2</sup> s<sup>-2</sup>)</td><td>U[-1000, 1000]</td><td>—</td><td>45<sup>+19</sup><sub>-10</sub></td><td>38.7<sup>+12.3</sup><sub>-08.5</sub></td></tr><tr><td>A<sub>FWHM</sub> (km s<sup>-1</sup>)</td><td>U[-1, 1]</td><td>0.149<sup>+0.076</sup><sub>-0.039</sub></td><td>—</td><td>—</td></tr><tr><td>A<sub>BIS</sub> (km s<sup>-1</sup>)</td><td>U[-1, 1]</td><td>-0.0048<sup>+0.0097</sup><sub>-0.0119</sub></td><td>-0.003<sup>+0.008</sup><sub>-0.0093</sub></td><td>—</td></tr><tr><td>B<sub>BIS</sub> (km s<sup>-1</sup> d)</td><td>U[-1, 1]</td><td>-0.06<sup>+0.016</sup><sub>-0.032</sub></td><td>-0.05<sup>+0.012</sup><sub>-0.02</sub></td><td>—</td></tr><tr><td colspan="5"><strong>Other parameters</strong></td></tr><tr><td>Offset RV (SERVAL, km s<sup>-1</sup>)</td><td>U[-1, 1]</td><td>—</td><td>-0.0972<sup>+0.0077</sup><sub>-0.005</sub></td><td>-0.0982<sup>+0.0047</sup><sub>-0.0043</sub></td></tr><tr><td>Offset RV (DRS, km s<sup>-1</sup>)</td><td>U[-5, -3]</td><td>-4.775<sup>+0.01</sup><sub>-0.013</sub></td><td>—</td><td>—</td></tr><tr><td>Offset FWHM (km s<sup>-1</sup>)</td><td>U[10, 12]</td><td>10.7<sup>+0.074</sup><sub>-0.085</sub></td><td>—</td><td>—</td></tr><tr><td>Offset DLW (m<sup>2</sup> s<sup>-2</sup>)</td><td>U[-1, 1]</td><td>—</td><td>-0.001<sup>+0.024</sup><sub>-0.021</sub></td><td>-0.003<sup>+0.018</sup><sub>-0.019</sub></td></tr><tr><td>Offset BIS (km s<sup>-1</sup>)</td><td>U[-1, 1]</td><td>0.0068<sup>+0.0076</sup><sub>-0.0073</sub></td><td>0.0072<sup>+0.0063</sup><sub>-0.007</sub></td><td>—</td></tr><tr><td>Jitter term σ<sub>j, RV</sub> (m s<sup>-1</sup>)</td><td>J[0, 1000]</td><td>13.14<sup>+2.26</sup><sub>-1.82</sub></td><td>10.32<sup>+1.68</sup><sub>-1.42</sub></td><td>10.3<sup>+1.6</sup><sub>-1.4</sub></td></tr><tr><td>Jitter term σ<sub>j, FWHM</sub> (m s<sup>-1</sup>)</td><td>J[0, 1000]</td><td>89.91<sup>+8.05</sup><sub>-7.0</sub></td><td>—</td><td>—</td></tr><tr><td>Jitter term σ<sub>j, DLW</sub> (1000 m<sup>2</sup> s<sup>-2</sup>)</td><td>J[0, 1000]</td><td>—</td><td>30.29<sup>+2.78</sup><sub>-2.39</sub></td><td>30.7<sup>+2.7</sup><sub>-2.4</sub></td></tr><tr><td>Jitter term σ<sub>j, BIS</sub> (m s<sup>-1</sup>)</td><td>J[0, 1000]</td><td>42.89<sup>+3.95</sup><sub>-3.52</sub></td><td>47.91<sup>+4.18</sup><sub>-3.73</sub></td><td>—</td></tr></tbody></table>

Note.ª, ᵇ, and ᶜ as defined in Table 3.

## APPENDIX C: MODEL COMPARISON FOR DIFFERENT PLANET CONFIGURATIONS

Table C1 shows a model comparison for different runs including different combinations of the AU Mic b and c planetary signals (we note that the activity model used corresponds to Model 3 in Appendix B).
<table><caption>Table C1. Model comparison.</caption><thead><tr><th>Model</th><th>N<sub>pars</sub></th><th>K<sub>b</sub> [m s<sup>-1</sup>]</th><th>K<sub>c</sub> [m s<sup>-1</sup>]</th><th>ln L</th><th>AIC<sup>a</sup></th><th>BIC<sup>b</sup></th></tr></thead><tbody><tr><td>Planet b signal only</td><td>15</td><td>6.1<sup>+2.9</sup><sub>-2.8</sub></td><td>–</td><td>336</td><td>-641</td><td>-595</td></tr><tr><td>Planet c signal only</td><td>15</td><td>–</td><td>8.4 ± 2.6</td><td>339</td><td>-648</td><td>-601</td></tr><tr><td>Both planetary signals</td><td>20</td><td>5.8 ± 2.5</td><td>8.5 ± 2.5</td><td>342</td><td>-643</td><td>-581</td></tr></tbody></table>

Note. $^{a}$Akaike Information Criterion ($AIC = 2N_{pars} - 2 \ln L$). $^{b}$Bayesian Information Criterion ($BIC = N_{pars} \ln N_{data} - 2 \ln L$).

## APPENDIX D: TTV PREDICTIONS

Fig. D1 shows how the amplitude of the expected TTVs changes when we vary the planet masses or eccentricities, keeping the

other parameters at the values reported in this work. Also shown are the TTVs expected for the nominal system parameters.

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

**Figure D1.** Expected TTV amplitudes for the AU Mic system. In the left-hand panel, we vary the eccentricities, while keeping the masses fixed. In the right-hand panel, we keep the eccentricities fixed but vary the masses. In both panels, the vertical and horizontal solid lines indicate the maximum a posteriori values from our modelling of the RVs. Note the difference in colour scale between the left- and right-hand panels. The explored range of eccentricities and masses extend slightly beyond the 1σ confidence limit in each direction. In the middle panel, the blue dots show the expected TTVs (time of transit centre relative to a linear ephemeris) when the planet parameters are set to their maximum a posteriori values. The x-axis corresponds to transit number. The green lines show a sinusoidal fit to those TTVs.

This paper has been typeset from a TEX/LATEX file prepared by the author.