© E. Martioli et al. 2021

1 Introduction

Young planetary systems represent an opportunity to observe planets in the early stages of planetary formation when gravitational interactions have not significantly changed the initial configuration of the system. AU Microscopii (AU Mic) is a particularly interesting young system, with an estimated age of 22 ± 3 Myr (Mamajek & Bell 2014), which is located at a distance of only 9.7248 ± 0.0046 pc (Gaia Collaboration 2018). Table 1 summarizes the star parameters of AU Mic that are relevant for this work. The AU Mic host is an M1 star with a spatially resolved edge-on debris disk (Kalas et al. 2004) and at least one transiting planet, AU Mic b (Plavchan et al. 2020). This Neptune-size planet is in a 8.5-d prograde orbit aligned with the stellar rotation axis (Martioli et al. 2020; Hirano et al. 2020; Palle et al. 2020). While Plavchan et al. (2020) only reported an upper limit on the mass of AU Mic b, Klein et al. (2021) measured ${17.1}_{-4.5}^{+4.7}$ M_⊕ thanks to infrared observations secured with the SPIRou spectropolarimeter (Donati et al. 2020). Plavchan et al. (2020) has also reported the detection of an isolated transit event of a second possible candidate planet, hereafter referred to as AU Mic c.

The notion that AU Mic could be a host for several planets is not a surprising one. Indeed, the occurrence rate of planets with a radius between 0.5 R_⊕ and 4.0 R_⊕ and a period between 0.5 and 256 d is estimated as ${8.4}_{-1.1}^{+1.2}$ planets per M dwarf (Hsu et al. 2020), although the occurrence rate of more massive planets in M dwarfs decreases significantly with mass (e.g., Bonfils et al. 2013). Moreover, provided the fact that AU Mic has an edge-on debris disk and that AU Mic b is a transiting planet with an aligned orbit, the chances that other planets also reside in co-planar orbits increase and, therefore, possible close-in additional planets are also likely to transit the star.

AU Mic is a magnetically active star with strong flaring activity (Robinson et al. 2001). Its surface is largely filled by starspots, producing a BY Draconis-type light curve with a quasi-periodic rotational modulation and a period of 4.863± 0.010 d (Plavchan et al. 2020). As in other active stars, the strong flaring and magnetic activity of AU Mic makes it more difficult to detect planetary transits in a photometric time series, especially for small planets where the occurrence rate is higher. In this work, we present an analysis of the Transiting Exoplanet Survey Satellite (TESS) data of AU Mic, including the new observations obtained in July 2020, where we implement a multi-flare model combined with the starspots model, improving the constraints on the planetary parameters from transit modeling and allowing for the detection of two additional transits of the second candidate planet AU Mic c.

2 TESS light curve

AU Mic was observed by TESS (Ricker et al. 2015) in Sector 1 from 2018-Jul.-25 to 2018-Aug.-22, in cycle 1, camera 1. Then it was observed again in Sector 27 from 2020-Jul.-04 to 2020-Jul.-30, in cycle 3, camera 1. The second visit was observed as part of the TESS Guest Investigator Programs G03263 (PI: P. Plavchan), G03141 (PI: E. Newton), and G03273 (PI: L. Vega) in fast mode with a time sampling of 20 s compared to the 2 min time sampling of the first visit. We obtained the de-trended photometrictime series from the Mikulski Archive for Space Telescopes (MAST) using the MAST astroquery tool¹. Figure 1 shows the AU Mic PDCSAP² flux in units of electrons-per-second (e⁻ s⁻¹) as a function of time given in TESS Barycentric Julian Date (TBJD = BJD – 2 457 000.0) for the two visits. Plavchan et al. (2020) have analyzed the same TESS data from the first visit only, and in the present paper we report for the first time an analysis including the more recent 2020 TESS data. The predicted times of the transits of AU Mic b and c as calculated from our ephemerides (see Sects. 5 and 6) are marked with vertical lines.

3 Starspots

As illustrated in Fig. 1, the TESS light curve of AU Mic shows a typical BY Draconis type of quasi-periodic variation due to starspots modulated by stellar rotation. The starspots evolve and, therefore, this modulation cannot be accurately modeled by a single periodic function. Thus, we treat this quasi-periodic modulation as a baseline “continuum” signal, where we model it by fitting a piece-wise fourth-order polynomial with iterative sigma-clipping. This approach works well when the ranges are carefully selected and inspected as follows. First, we can be certain that the polynomial function is a sufficiently accurate approximation for the local baseline variation within each range. Then we avoid the edges of the ranges to include either a transit or a flare event. Finally, we avoid the inclusion of large gaps (when TESS did not deliver data) within the same range. The edges for all selected ranges are represented by the blue vertical dashed lines in Fig. 1.

This starspot model is removed from the data to model flares and planetary transits as explained in Sects. 4–6, which are further removed from the original data to fit the starspot modulation again. We repeat this procedure iteratively until the standard deviation of residuals is not improved by more than 1%. We typically meet this criterion upon three to five iterations. The final starspot model is also shown in Fig. 1. In Appendix A, we present an independent measurement of the rotational period of AU Mic from the starspot modulation in the 2018 TESS data.

4 Flares

AU Mic is an active young star with intense flaring activity. We carry out an analysis of the flares in the TESS data to improve the constraints on the transit events. We consider the starspot-subtracted residuals as shown in Fig. 2. The flares are typically clusters of points lying above the noise level. First we detect peaks as possible candidate flares using the scipy.signal.find_peaks³ routine, where it finds local maxima via a simple comparison of neighbouring values. We set a minimum peak amplitude of 2.5σ and a minimum horizontal distance between peaks of about ~100 min. Some flares in AU Mic are quite complex, while other flares occur before a previously initiated flare has ended. Therefore, we visually inspected all detected peaks to identify any possible additional peaks that have been missed by the find-peaks algorithm. After a few iterations, we identified a total of 324 flares in the whole TESS time series.

The times and amplitudes of detected peaks were adopted as initial values for the basis of performing a least-squares fit using the multi-flare model of Davenport et al. (2014), where each flare is represented by a two-phase model with a polynomial rise and an exponential decay. Figure 3 shows, as an example, the residual light curve and fit flares model for three days of flaring activity in AU Mic in both TESS visits. For each flare in the model, we fit the flare amplitude (f_max), the full width at half maximum (FWHM), and the time of maximum t_peak. The fit parameters for all flares are presented in Appendix B.

The total continuous monitoring time of AU Mic performed by TESS was estimated as 49.75 d, where we have discounted the large gaps in the data. We detected a total of 324 flares, from which 316 have fit amplitudes above F_max > 1σ, for σ = 161.5 e⁻s⁻¹. Thus, the AU Mic global occurrence rate of flares with amplitude above f_max ≳ 0.06% of the median stellar flux (279544 e⁻s⁻¹) is estimatedas 6.35 flares per day or about 1 flare every 3.8 h. Figure 4 shows the distribution of flares as function of flare amplitudes. We fit an exponential to the measured distribution and obtained an empirical model for the occurrence rate as function of flare amplitude of f_r = 0.17e^−4.0f_max d⁻¹, with f_max in units of 10² e⁻s⁻¹. We notice that the strong flares, which are more important energetically, have a much lower occurrence rate and therefore a low statistical significance in our analysis. Thus, this empirical relationship is only expected to hold for flares above the detection limit imposed by the noise and for flares with a statistically significant number of events observed in the TESS data, that is, those with amplitudes in the range of 0.06% < f_max < 1.5% of the star flux.

The transit duration of AU Mic b calculated by Plavchan et al. (2020) is about 3.5 h with transit depth of 0.26%. Therefore, considering that a transit observation of AU Mic b requires an out-of-transit baseline of about 2 h, it is most likely that any transit observation in AU Mic should be affected in some way by flare events. This shows the importance of modeling flares adequately to improve the constraints on the planetary parameters obtained by the transits.

Fig. 1 TESS light curve of AU Mic. Top panel: data for the first visit (from 2018-Jul.-25 to 2018-Aug.-22). Bottom panel: data for the second visit (from 2020-Jul.-04 to 2020-Jul.-30). Black points are the PDCSAP flux and vertical solid lines show the times of planetary transits, as obtained in our analysis, for AU Mic b (in orange) and AU Mic c (in green). Blue vertical dashed lines show the knots of selected ranges for the piece-wise polynomial fit and the red line shows the starspot model.

5 Transits of AU Mic b

Three transits of AU Mic b occurred during the first TESS visit, but only two have been observed due to a communication problem with the spacecraft during the second transit (see Plavchan et al. 2020). For the second visit, three other transits were expected, and indeed we were able to identify all of them. To analyze these five transits together, westart with the parameters of AU Mic b from Plavchan et al. (2020), referred to as “prior values” in Table 2, to remove the transit signal from the TESS light curve before fitting the starspot and flare models. Our transit model is calculated using the BATMAN toolkit by Kreidberg (2015), where we assume a circular orbit (e = 0). This assumption is in agreement with Plavchan et al. (2020), who did not detect any significant eccentricity from their analysis of the transits only, where the eccentricity would be slightly constrained by the duration of the transit. We adopt as priors the quadratic limb-darkening coefficients (LDC) from Claret (2018) with an arbitrary error of 0.1, where we obtained the values calculated for the photometric system of TESS and those matching the closest stellar parameters to those of AU Mic (see Table 1).

The transits of AU Mic b, flares, and starspots are fit simultaneously using the non-linear least squares optimization (OLS) fit tool scipy.optimize.leastsq. As explained in Sect. 3 this procedure is run iteratively, where we first obtain an independent fit for each component of the model and then we set these values as initial guess to run an OLS analysis with all free parameters in the three components of the model. Then we consider the data in the ranges around the transits of AU Mic b to sample the posterior distributions of the transit parameters using the emcee Markov chain Monte Carlo (MCMC) package (Foreman-Mackey et al. 2013). The chain is set with 100 walkers and 30 000 MCMC steps of which we discard the first 5000. The MCMC samples are presented in Fig. C.1, which show the chains reaching stability before the first 5000 discarded steps. The posterior distributions are illustrated in Fig. C.3. The best-fit values of the transit parameters are calculated as the medians of the posterior distribution with error bars defined by the 34% on each side of the median, all presented in Table 2.

Figure 5 shows the results of our analysis for each transit range separately, and the bottom right panel shows the flux normalized by the starspot and flare models for all transits together, and our best fit model along with the previous model of Plavchan et al. (2020) for comparison. The dispersion of residuals are 358, 433, 633, 566, and 588 ppm, for each respective epoch. The global dispersion is 573 ppm. Notice that our measured parameters of AU Mic b agree within 3σ with the previous measurements by Plavchan et al. (2020) but with improved accuracy. Our measured planet-to-star radius ratio is slightly larger than that measured by Plavchan et al. (2020; see lower, right panel of Fig. 5). However, as described in Sect. 7, our derived effective planetary radius of 1.05 ± 0.04 R_Nep is slightly smaller than the value 1.08 ± 0.05 R_Nep reported by Plavchan et al. (2020).

Fig. 2 TESS residual light curve of AU Mic. Top panel: residual light curve for the first TESS visit (2018-Jul.-25 to 2018-Aug.-22). Bottom panel: residual curve for the second visit (2020-Jul.-04 to 2020-Jul.-30). Blue circles show AU Mic TESS flux data after subtracting the starspot model and black dotted lines show ± 2.5σ range around the starspot fit model, for σ = 162.5 e−s−1. The flaring activity can be seen as positive groups of points above + 2.5σ and the transits as negative groups of points below − 2.5σ. Our predicted times of transits for AU Mic b (in orange) and AU Mic c (in green) are marked with arrows.

Fig. 3 TESS residual light curve of AU Mic between TBJD = 1344.3 and TBJD = 1347.3 (top panel) and between TBJD = 2050.8 and TBJD = 2053.8 (bottom panel). Blue circles show the flux after subtracting the starspot model, brown line shows the multi-flare best fit model. Green points show the residuals with an arbitrary offset for better visualization. The dispersion of residuals are 107 e− s−1 and 178 e−s−1 for data presented in top panel and bottom panel, respectively.

Fig. 4 Number of flares normalized by the total observed time of 49.75 d as function of flare amplitude. Dashed line shows the fit to an exponential decay as shown in the legend, where variable x represents the flare amplitude.

6 Confirmation and characterization of planet AU Mic c

An isolated candidate transit event at T_c = 1342.22 ± 0.03 TBJD was identified in the 2018 TESS data by Plavchan et al. (2020), which was interpreted as a possible second planet in a 30 ± 6-d orbit. We have searched for other transit-like events in the TESS data by looking at regions with systematic flux dips below the noise level. We identified three regions as potential transits including the one previously identified by Plavchan et al. (2020) (see Fig. 2). We notice that there is an additional region presenting a significant flux dip around TBJD = 1347.8, which did not include as a possible transit. The 2018 TESS data after TBJD ~ 1347.7 contain several gaps, which have an impact on the modeling of the baseline modulation by starspots as well as on the detection and modeling of flares. Therefore, we presume that the apparent flux dip in this region is more likely due to a misfit other than a transit. In addition the shape of that dip is different from that of the three features above.

We hypothesize that the three selected events are transits of the same candidate planet AU Mic c. In order to test this hypothesis, we first obtain an orbital period that is consistent with the three transits, that is, a period given by the slope of a linear fit to the ephemeris equation, T_c = T₀ + P × E, with T_c as the times of conjunction measured independently for each event and E as the corresponding epochs to each event, assumed to be 0, 37, and 38. We find an orbital period of P = 18.85895 ± 0.00003 d. However, we note that with this particular periodicity, only three events would occur during the TESS observations (see Fig. 1). The possibility of shorter alias periods is ruled out by the fact that TESS would have observed more transits of this object and we did not detect them. To test further our hypothesis, we performed an independent analysis of each event assuming an orbital period, first with an uniform prior of $P=\mathcal{U}(1,500)$ d to explore a broad range of periods, and then with a normal prior of $P=\mathcal{N}(18.86,0.01)$ d, which explores a solution constrained by the previous knowledge imposed by our hypothesis. Analogously, the normalized semi-major axes are also first set with an uniform prior of $a/R_{\star }=\mathcal{U}(1,500)$ and then with a normal prior estimated using the Kepler’s law, that is, $a/R_{\star }=\mathcal{N}(31.5,1.4)$. The priors for the times of conjunction are measured locally on each event, for example, for the first event, we obtained $T_0=\mathcal{N}(1342.23,0.01)$. The prior for the planet-to-star radius ratio is estimated from the square root of the average flux depth of the three events, assuming a conservative error, that is, $R_{\text{p}}/R_{\star }=\mathcal{N}(0.04,0.01)$. The orbital inclination is considered with an uniform prior distribution of $i=\mathcal{U}({85}^{°},{90}^{°})$ and the limb darkening coefficients are fixed to the literature values presented in Table 1. In Appendix D, we present the results for this independent analysis of each transit. The posterior of all the transit parameters agree within 3σ (see Table D.1), and the dispersion of residuals improves when we adopted the fit period of $P=\mathcal{N}(18.86,0.01)$ d as prior, which supports our hypothesis that these three events were caused by the transits of the same planet AU Mic c. Finally, we perform a joint analysis of the three transits simultaneously, where we adopt the same priors as in the independent analysis above, except for the limb darkening coefficients, where we adopted a normal prior with the literature values and an arbitrary error of 1.0. We call attention to the fact that the limb darkening coefficients obtained in Sect. 5 could have been used as priors since both planets are transiting the same star. However, the two planets may transit the stellar disk in different regions, and since AU Mic is largely filled by starspots, the limb darkening may also be affected by the different temperatures of the transited regions in the photosphere. Table 3 presents the priors and fit parameters, while Fig. 6 shows the TESS data for the three transits and the best-fit model obtained from our analysis. As in Sect. 5, the MCMC samples and posterior distributions are illustrated in Figs. C.2 and C.4.

Fig. 5 Each panel shows the range around one of the three transits of AU Mic b observed by TESS. The top part of each panel shows the TESS data (blue circles) and the best-fit model for each component, as indicated in the legend. The bottom partin each panel shows the TESS data normalized by the starspot and flare model (blue points) and the best-fit transit model (red line). Vertical dashed lines show the central time of each transit. Bottom right panel: phase-folded light curve for all the data from the five transits (light blue points) and the binned data (dark blue points) with a bin size of 0.005 d, our best fit model (red line), and the previous fit model of Plavchan et al. (2020) for comparison (orange dashed line). The residuals (green points) are displayed with an arbitrary offset.

7 Discussion

We combinethe transit parameters of AU Mic b and c presented in Tables 2 and 3 with the star parameters from Table 1 to calculate a number of derived parameters, which are presented in Table 4.

The planet-to-star radius ratio obtained in Sects. 5 and 6 are likely overestimated due to the presence of spots covering a significant fraction of the stellar photosphere. This is because the observed transit depth is likely smaller than the nominal depth that would be obtained for a completely unspotted photosphere. Therefore, an effective planet-to-star radius ratio can be obtained by applying a correction factor to the observed transit depth as in Eq. (1) of Rackham et al. (2018), that is, $${\left(\frac{R_{\text{p}}}{R_{\star }}\right)}_{\text{eff}}^2=\left[1-f_{\text{spot}}\left(1-\frac{F_{\text{spot}}}{F_{\text{phot}}}\right)\right]{\left(\frac{R_{\text{p}}}{R_{\star }}\right)}_{\text{obs}}^2$$(1)

for f_spot being the spot filling factor estimated at 20% based on magnetic activity (Klein et al. 2021) and F_spot∕F_phot being the flux fraction between the spots and the photosphere. The latter can be estimated assuming a blackbody flux for each of the two components, where the photospheric temperature is T_phot = 3700 K and the spots’ temperature is estimated at 86% of the photospheric value in M1 stars (Rackham et al. 2018). The latter is in agreement with spectropolarimetric measurements of AU Mic by Berdyugina (2011). Our estimation for the spot-to-photosphere flux fraction integrated over the TESS band pass (between 60 and 1000 nm) is F_spot∕F_phot ~ 46%. Therefore, the effective radius for both transiting planets in the AU Mic system should be given by $R_{\text{p}}^{\text{eff}}=0.94R_{\text{p}}^{\text{obs}}$, that is, about 6% smaller than the measured radius.

In addition to the correction above, there is an uncertainty in the planet radius due to the rotational modulation by starspots and other photospheric heterogeneities, which implies a variable baseline flux. To account for this additional uncertainty, we considered the amplitude of the photospheric modulation in the AU Mic TESS light curve, which is about 5%. This variability propagates an uncertainty on the order of $~\sqrt{0.05}$ to the planet radius measurements, which is included in the uncertainty values reported in Table 4.

Klein et al. (2021) recently measured a semi-amplitude of the radial velocity of the reflex motion caused by AU Mic b of $K_{\text{b}}={8.5}_{-2.2}^{+2.3}$ m s⁻¹, providing an important constraint on the mass of this planet. We recalculated the planet density based on our radius measurement and obtained ρ_b = 1.4 ± 0.4 g cm⁻³.

To estimate a plausible mass range for AU Mic c, we considered a population of transiting exoplanets⁴ with radii within 1σ of the measured radius of AU Mic c, as presented in Fig. 7. The median and median absolute deviation of the mass of exoplanets in this population is 0.043 ± 0.036 M_Jup. Therefore, the mass of AU Mic c is likely to lie in the range of 0.007 M_Jup < M_c <  0.079 M_Jup, implying a planet density in the range of 0.4 g cm⁻³ < ρ_c <  4.1 g cm⁻³. Several different scenarios are possible for the internal structure of this planet depending on its mass. We estimate the semi-amplitude of the induced radial velocity caused by AU Mic c in the star motion to be in the range of 0.8 m s⁻¹ < K_c <  9.5 m s⁻¹. A number of world-class spectrometers can currently achieve the necessary precision to detect such reflex motion. However, the intense stellar activity in AU Mic may represent a challenge for such measurement.

In order to analyze the stability of the orbital solution (Table 4), we performed a global frequency analysis (Laskar 1990, 1993) in the vicinity of the best fit, in the same way as achieved for other planetary systems (e.g., Correia et al. 2005, 2010). The system is integrated on a regular 2D mesh of initial conditions, with varying semi-major axis and eccentricity of planet c, while the other parameters are retained at their nominal values. We used the symplectic integrator SABA1064 of Farrés et al. (2013), with a step size of 5 × 10⁻³ yr and general relativity corrections. Each initial condition is integrated over 5 kyr, and a stability indicator is computed to be the variation in the measured mean motion over the two consecutive 2.5 kyr intervals of time (for more details see Couetdic et al. 2010). For regular motion, there is no significant variation in the mean motion along the trajectory, while it can vary significantly for chaotic trajectories.

In Fig. 8, we show the wide vicinity of the nominal solution for two different configurations, one with M_c = 25.0 M_⊕ (top) and another with M_c = 2.2 M_⊕ (bottom), corresponding to the maximum and minimum masses of the outer planet, respectively. The stability indicator is reported using a color index, where “red” represents the strongly chaotic trajectories and “dark-blue” shows the extremely stable ones. We observe that there are several islands of mean motion resonances nearby, however, the main difference is that when the mass of the outer planet is large, these resonances become unstable. The best fit nominal solution (vertical dotted line) is close to the 9:4 resonance, but outside, and so the nominal solution is in the stable zone for both mass values, as long as e_c < 0.2. We hence conclude that the AU Mic planetary system is stable and nearly circular.

The proximity to the 9:4 resonance could cause significant transit-timing variation (TTV). To measure the TTV for both planets, we run an MCMC analysis for each individual transit observed by TESS, where we obtained an independentmeasurement of the times of transits. The transit parameters are fixed to their best fit values, except T_c, which is set with an uniform prior distribution of $\mathcal{U}(T_{\text{c}}-0.005,T_{\text{c}}+0.005)$, for T_c being the time of conjunction calculated from the fit ephemeris. We subtract T_c from the measured times of conjunction to obtain the TTVs as presented in Table 5. We detect no significant TTV in the TESS data. TTVs with amplitudes of one minute or more would have been likely detected in this dataset.

Using the nominal solution from Table 4, we generated the synthetic TTVs as in Hébrard et al. (2020), for the minimum and the maximum masses of the outer planet. In Fig. 9, we show the variations corresponding to each solution superimposed with the observational data (Table 5). For the outer planet, the amplitude of the TTV is around one minute for both mass choices. However, for the inner planet, the minimum mass produces TTV with an amplitude of only a few seconds, while for the maximum mass the amplitude can reach up to three minutes. Thus, the TTV of the inner planet can place constraints on the mass of the outer planet.

The precision and the number of photometric measurements currently available for the AU Mic system (Table 5) do not allowus to run an exhaustive search for a best fit solution, but they allow us to reduce the uncertainty in the mass of the outer planet. The non-detection of TTVs that are larger than one minute place the mass of the outer planet on the low side as the likelier assumption. In Fig. 9, we show the TTV corresponding to an outer companion with M_c = 7.0 M_⊕. This mass generates TTVs with an amplitude around one minute corresponding to our upper limit. The observation of additional transits should help to resolve the present ambiguity in the outer planet mass.

Furthermore, we calculated the habitable zone (HZ) for AU Mic using the equations and data from Kopparapu et al. (2014), which gives an optimistic lower limit (recent Venus) at 0.25 au, and an upper limit (early Mars) at 0.64 au, with the runaway greenhouse limits (M_p = 1 M_⊕) ranging between 0.32 and 0.61 au. AU Mic b and c being at an orbital distance of 0.0645 ± 0.0013 au and 0.1101 ± 0.0022 au are not in the HZ. We estimate the equilibrium temperatures as in Heng & Demory (2013) assuming an uniform heat redistribution and an arbitrary geometric albedo of 0.1, which gives T_{eq, b} = 593 ± 21 K and T_{eq, c} = 454 ± 16 K, for AU Mic b and c, respectively, showing that these warm Neptunes are indeed too hot to sustain water in liquid state. Moreover, the actual temperature for these planets could be much higher due to a possible greenhouse effect depending on their atmospheric composition and also due to an elevated internal temperature provided the young age of the system.

Fig. 6 Each panel shows the range around one of the three transits of AU Mic b observed by TESS. Details are the same as in Fig. 5, but here we present the results for the three TESS transits of AU Mic c. The dispersion of residuals are 370, 646, 609 ppm, for each respective epoch. The global dispersion is 609 ppm.

Fig. 7 Mass-radius diagram of planets. Blue circles show the known transiting exoplanets compiled from exoplanet.eu. The orange point represents AU Mic b and the light green horizontal dashed and dotted lines represent the AU Mic c radius measurement and ± 1σ uncertainties. Dark green point represents the median of planets with radius within 1σ of the measured radius of AU Mic c. Dashed orange line shows an iso-density line for the radius and mass of AU Mic b and the olive lines correspond to the iso-densities obtained for the radius of AU Mic c and for the lower limit (dotted), median value (dashed), and upper limit (dot-dashed) of the mass of AU Mic c, based on the exoplanet population with the same radii.

8 Conclusions

We presented an analysis of the photometric TESS observations of the active young M1 star AU Mic, where we model the rotational modulation by starspots, the flaring activity, and the planetary transits simultaneously. Our analysis delivered anestimation of the flare occurrence rate of 6.35 flares per day for flares with amplitudes in the range of 0.06% < f_max < 1.5% of the star flux. With such a flare rate it is important to model flaring activity to improve the constraints on the planetary parameters from transits. In our simple multi-flare model, we did not explore a detailed physical modeling of flares (e.g., as in Tilley et al. 2019) that could, in principle, be studied more extensively given the high quality and broad time coverage of these TESS observations. Our aim was to find an empirical description of flares that improves the constraints on the planetary parameters of AU Mic b and that increases the sensitivity for detecting and characterizing the second planet AU Mic c.

Our analysis of the five transits of AU Mic b provided measurements of the orbital period and time of conjunction, giving an improved ephemeris ofT_c = 2 458 330.39051 ± 0.00015 + E × 8.463000 ± 0.000002 JD for an accurate prediction of the times of future transits. Our detection of two new transits in the recent TESS observations have provided substantial evidence to confirm the second planet, that is, AU Mic c. We analyzed the three transits and obtained consistent results for the planetary parameters, which supports our hypothesis that these transits come from the same planet. Our MCMC analysis provided a determination of the planetary parameters for AU Mic c, delivering an ephemeris of $T_{\text{c}}=2458342.2223\pm 0.0005+E\times {18.859019}_{-0.000015}^{+0.000016}$ JD for the times of transits with a transit duration of 4.5 ± 0.8 hours. The derived radius of 3.24 ± 0.16 R_⊕ indicates that AU Mic c is a Neptune-size planet with an expected mass in the range of 2.2 M_⊕ < M_c <  25.0 M_⊕, estimated from the population of exoplanets of a similar size.

AU Mic b and c with an orbital period of 8.46 d and 18.86 d are near the 9:4 mean-motion resonance. Such a configuration is stable and could cause significant TTV, which we do not detect with the current dataset. Still, that non-detection provides an upper limit to the mass of AU Mic c of M_c < 7 M_⊕, implying an upper limit to the planet density of ρ_c <1.3 g cm⁻³. This suggests that both planets are likely inflated, as it should be expected for young planets (e.g., Helled et al. 2020). Therefore, these planets are interesting targets for atmospheric characterization by transmission spectroscopy.

The interactions between the planets coupled in a near 9:4 mean-motion resonance and the debris disk surrounding the star are to be the subject of subsequent studies. High-resolution images have revealed that the disk has a complex structure with small-sized substructures, such as “feature A” at a dozen astronomical units showing a “loop-like” morphology (Wisniewski et al. 2019). As for the case of β Pictoris, it has been shown that the planets can sculpt the disk morphology through the gravitational interaction between planets and the debris disk parents bodies at far distances (Lecavelier Des Etangs et al. 1996; Augereau et al. 2001), particularly through resonance mechanisms in young systems that are still evolving (Lecavelier Des Etangs 1998). The impact of the presence of AU Mic b and AU Mic c on the dust disk deserves further analysis.

Fig. 8 Stability analysis of the AU Mic planetary system, around the best fit solution (Table 4), with Mc = 25.0 M⊕ (a), and Mc = 2.2 M⊕ (b). The phase space of the system is explored by varying the semi-major axis ac (and thus the period Pc) and eccentricity ec of the outer planet; ec is plotted as a function of Pc and not ac, as the uncertainty on ac is larger. For the initial conditions, the step size is 0.01 in eccentricity, and 2 × 10−5 au in semi-major axis. For each initial condition, the system is integrated over 5 kyr and a stability criterion is derived with the frequency analysis of the mean longitude. The chaotic diffusion is measured by the variation in the mean motion. The color scale corresponds to the decimal logarithm of the variation of the mean motion (Correia et al. 2010). The “red” zone corresponds to highly unstable orbits, while the “dark-blue” region canbe assumed to be stable on a billion-years timescale. The main resonances in the vicinity of the solution (13:6, 11:5, 9:4, 7:3) are labeled. It should be noted that for the lower value of the mass Mc (b), the libration island is mostly stable while it is not the case for the larger value of Mc (a).

Fig. 9 Transit timing variations for AU Mic b (top) and AU Mic c (bottom). The colored lines correspond to synthetic data obtained with the best fit solution (Table 4) using different values for the mass of the outer planet: Mc = 2.2 M⊕ (blue), Mc = 7.0 M⊕ (green), and Mc = 25.0 M⊕ (red). The dots correspond to the observed TTVs. (Table 5).

Acknowledgements

We acknowledge funding from the French National Research Agency (ANR) under contract number ANR-18-CE31-0019 (SPlaSH). This paper includes data collected with the TESS mission, obtained from the MAST data archive at the Space Telescope Science Institute (STScI). In particular, the TESS data secured in 2020 are part of the TESS Guest Investigator Programs G03263 (PI: P. Plavchan), G03141 (PI: E. Newton), and G03273 (PI: L. Vega). Together with other data and analyses of that particularly interesting system, that will be the subjects of forthcoming papers by Cale et al., Collins et al., El Mufti et al., Gilbert et al., Wittrock et al. Funding for the TESS mission is provided by the NASA Explorer Program. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. We thank the TESS Team members for making available the extremely accurate photometric data they obtained. We acknowledge support by CFisUC projects (UIDB/04564/2020 and UIDP/04564/2020), GRAVITY (PTDC/FIS-AST/7002/2020), ENGAGE SKA (POCI-01-0145-FEDER-022217), and PHOBOS (POCI-01-0145-FEDER-029932), funded by COMPETE 2020 and FCT, Portugal. This work benefited from HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir.

Appendix A Star rotation

In this appendix, we present an independent measurement of the star rotation period of AU Mic.

Starspots in AU Mic have a life time sufficiently large that the overall spot pattern did not change significantly during the 2018 TESS observations in the first visit. As one can see in the bottom panel of Fig. 1, there has been a more significant change in the spot pattern during the second visit of TESS, where the amplitude of the cyclic flux variations decreases with time. For this reason, we did not include the second visit data to measure the rotation period. Flares and planetary transits are removedfrom the original data using the models computed in Sects. 4–6. We apply the phase dispersion minimization (PDM) method inspired by the work of Stellingwerf (1978), where we calculate the phase curve for several values of trial rotation periods ranging from 4.85 d to 4.88 d in steps of 0.0001 d. The rotation period of P_r = 4.865 d measured by Torres et al. (1972) was considered as initial guess to define our search range. For each trial period, we fit a cubic spline with 15 knots and calculate the chi-square by ${\chi }^2={\displaystyle \sum {(F_{\text{i}}-F_{\text{c}})}^2/{\sigma }_{\text{i}}^2}$, for F_i being the flux after removing flares and transits, F_c being the fit model, and σ_i being the flux uncertainty given by the TESS pipeline. The results are presented in Fig. A.1, where we find a minimum χ² at P_r = 4.862 ± 0.032 d, which is consistent but less accurate than the rotation period of 4.863 ± 0.010 d obtained byPlavchan et al. (2020) from the same data set. The resulting phase curve and best fit model are presented in Fig. A.2. The final dispersion of residuals is 265 e⁻s⁻¹, which is more than twice the dispersion of residuals from the fit to the time series, that is, 99 e⁻s⁻¹ for the data from the first visit of TESS only. This discrepancy is likely due to the evolution or drifts of starspots during the several rotation cycles covered in the first visit and perhaps also due to differential rotation (Klein et al. 2021), which could generate a dispersion in the phase diagram since starspots at different latitudes are modulated by different periods.

Fig. A.1 Rotation period versus chi-square for the reduced light curve of AU Mic after removing flares and planetary transits. Red dashed line shows the minimum χ2 at P = 4.862 ± 0.032 d obtained by a parabolic fit (green line).

Fig. A.2 Phased light curve of AU Mic using our measured rotation period of 4.862 d. Black points show the 2018 TESS (first visit only) fluxes after removing the flares and transits best fit models, and red line shows the best fit cubic spline model to the phased data.

Assuming the radius of AU Mic measured by interferometry of R_⋆ = 0.75 ± 0.03 R_⊙ (White et al. 2015) and the rotation period of P_r = 4.862 ± 0.032 d, we calculate the velocity at the equator as v_eq = 2πR_⋆∕P_r = 7.8 ± 0.3 km s⁻¹. This gives a maximum value of the projected velocity of v_eq sin(i_⋆) < 8.1 km s⁻¹, for i_⋆ being the inclination angle of the rotation axis of the star with respect to the line of sight.

Appendix B Flares fit parameters

Appendix C MCMC samples and posterior distributions

Fig. C.1 MCMC samples (black dots) used to calculate the posterior distributions of the transit parameters of AU Mic b. The red dots show the discarded samples corresponding to the first 5000 steps of the chain. The blue lines show the median calculated for bins of 1000 steps. Dashed lines show the 16- and 84-percentiles and the solid horizontal lines shows the mediansof the valid samples.)

Fig. C.2 Same Fig. C.1 for the MCMC samples of the transit parameters of AU Mic c.

Fig. C.3 Pairs plot showing the MCMC samples and posterior distributions for the transit parameters of AU Mic b as presented in Table 2. The contours mark the 1, 2, and 3σ regions of the distribution. The blue crosses indicate the best fit values for each parameter and the dashed vertical lines in the projected distributions show the median values and the 1σ uncertainty (34% on each side of the median).

Fig. C.4 Pairs plot showing the MCMC samples and posterior distributions for the transit parameters of AU Mic b as presented in Table 2. Details are the same as in Fig. C.3 for the analysis of the three transits of AU Mic c.

Appendix D Independent fit to the transits of AU Mic c

This appendix presents the results obtained from an MCMC analysis (as in Sect. 5) performed on each transit of AU Mic c independently. Table D.1 shows the posterior of the transit parameters, where we present first the results using broad priors and then using the more restricted priors, as described in Sect. 6. We notice that the period and semi-major axis are not well constrained by the shape of the transits. As a consequence, when using broad priors, their posteriors turn out to have relatively large uncertainties. However, the posterior of all parameters agree within 3σ among the three events and they also agree with our derived period of P ~ 18.86-d. The fit parameters for the more constrained priors also agree within 3σ, and the dispersion of residuals are improved with respect to the unconstrained solution, which supports our hypothesis that these events are caused by the transits of the same planet, namely, AU Mic c.

References

All Tables

All Figures

Fig. 1 TESS light curve of AU Mic. Top panel: data for the first visit (from 2018-Jul.-25 to 2018-Aug.-22). Bottom panel: data for the second visit (from 2020-Jul.-04 to 2020-Jul.-30). Black points are the PDCSAP flux and vertical solid lines show the times of planetary transits, as obtained in our analysis, for AU Mic b (in orange) and AU Mic c (in green). Blue vertical dashed lines show the knots of selected ranges for the piece-wise polynomial fit and the red line shows the starspot model.

In the text

Fig. 2 TESS residual light curve of AU Mic. Top panel: residual light curve for the first TESS visit (2018-Jul.-25 to 2018-Aug.-22). Bottom panel: residual curve for the second visit (2020-Jul.-04 to 2020-Jul.-30). Blue circles show AU Mic TESS flux data after subtracting the starspot model and black dotted lines show ± 2.5σ range around the starspot fit model, for σ = 162.5 e−s−1. The flaring activity can be seen as positive groups of points above + 2.5σ and the transits as negative groups of points below − 2.5σ. Our predicted times of transits for AU Mic b (in orange) and AU Mic c (in green) are marked with arrows.

In the text

Fig. 3 TESS residual light curve of AU Mic between TBJD = 1344.3 and TBJD = 1347.3 (top panel) and between TBJD = 2050.8 and TBJD = 2053.8 (bottom panel). Blue circles show the flux after subtracting the starspot model, brown line shows the multi-flare best fit model. Green points show the residuals with an arbitrary offset for better visualization. The dispersion of residuals are 107 e− s−1 and 178 e−s−1 for data presented in top panel and bottom panel, respectively.

In the text

	Fig. 4 Number of flares normalized by the total observed time of 49.75 d as function of flare amplitude. Dashed line shows the fit to an exponential decay as shown in the legend, where variable x represents the flare amplitude.
In the text

Fig. 5 Each panel shows the range around one of the three transits of AU Mic b observed by TESS. The top part of each panel shows the TESS data (blue circles) and the best-fit model for each component, as indicated in the legend. The bottom partin each panel shows the TESS data normalized by the starspot and flare model (blue points) and the best-fit transit model (red line). Vertical dashed lines show the central time of each transit. Bottom right panel: phase-folded light curve for all the data from the five transits (light blue points) and the binned data (dark blue points) with a bin size of 0.005 d, our best fit model (red line), and the previous fit model of Plavchan et al. (2020) for comparison (orange dashed line). The residuals (green points) are displayed with an arbitrary offset.

In the text

	Fig. 6 Each panel shows the range around one of the three transits of AU Mic b observed by TESS. Details are the same as in Fig. 5, but here we present the results for the three TESS transits of AU Mic c. The dispersion of residuals are 370, 646, 609 ppm, for each respective epoch. The global dispersion is 609 ppm.
In the text

Fig. 7 Mass-radius diagram of planets. Blue circles show the known transiting exoplanets compiled from exoplanet.eu. The orange point represents AU Mic b and the light green horizontal dashed and dotted lines represent the AU Mic c radius measurement and ± 1σ uncertainties. Dark green point represents the median of planets with radius within 1σ of the measured radius of AU Mic c. Dashed orange line shows an iso-density line for the radius and mass of AU Mic b and the olive lines correspond to the iso-densities obtained for the radius of AU Mic c and for the lower limit (dotted), median value (dashed), and upper limit (dot-dashed) of the mass of AU Mic c, based on the exoplanet population with the same radii.

In the text

Fig. 8 Stability analysis of the AU Mic planetary system, around the best fit solution (Table 4), with Mc = 25.0 M⊕ (a), and Mc = 2.2 M⊕ (b). The phase space of the system is explored by varying the semi-major axis ac (and thus the period Pc) and eccentricity ec of the outer planet; ec is plotted as a function of Pc and not ac, as the uncertainty on ac is larger. For the initial conditions, the step size is 0.01 in eccentricity, and 2 × 10−5 au in semi-major axis. For each initial condition, the system is integrated over 5 kyr and a stability criterion is derived with the frequency analysis of the mean longitude. The chaotic diffusion is measured by the variation in the mean motion. The color scale corresponds to the decimal logarithm of the variation of the mean motion (Correia et al. 2010). The “red” zone corresponds to highly unstable orbits, while the “dark-blue” region canbe assumed to be stable on a billion-years timescale. The main resonances in the vicinity of the solution (13:6, 11:5, 9:4, 7:3) are labeled. It should be noted that for the lower value of the mass Mc (b), the libration island is mostly stable while it is not the case for the larger value of Mc (a).

In the text

	Fig. 9 Transit timing variations for AU Mic b (top) and AU Mic c (bottom). The colored lines correspond to synthetic data obtained with the best fit solution (Table 4) using different values for the mass of the outer planet: Mc = 2.2 M⊕ (blue), Mc = 7.0 M⊕ (green), and Mc = 25.0 M⊕ (red). The dots correspond to the observed TTVs. (Table 5).
In the text

	Fig. A.1 Rotation period versus chi-square for the reduced light curve of AU Mic after removing flares and planetary transits. Red dashed line shows the minimum χ2 at P = 4.862 ± 0.032 d obtained by a parabolic fit (green line).
In the text

	Fig. A.2 Phased light curve of AU Mic using our measured rotation period of 4.862 d. Black points show the 2018 TESS (first visit only) fluxes after removing the flares and transits best fit models, and red line shows the best fit cubic spline model to the phased data.
In the text

	Fig. C.1 MCMC samples (black dots) used to calculate the posterior distributions of the transit parameters of AU Mic b. The red dots show the discarded samples corresponding to the first 5000 steps of the chain. The blue lines show the median calculated for bins of 1000 steps. Dashed lines show the 16- and 84-percentiles and the solid horizontal lines shows the mediansof the valid samples.)
In the text

	Fig. C.2 Same Fig. C.1 for the MCMC samples of the transit parameters of AU Mic c.
In the text

Fig. C.3 Pairs plot showing the MCMC samples and posterior distributions for the transit parameters of AU Mic b as presented in Table 2. The contours mark the 1, 2, and 3σ regions of the distribution. The blue crosses indicate the best fit values for each parameter and the dashed vertical lines in the projected distributions show the median values and the 1σ uncertainty (34% on each side of the median).

In the text

	Fig. C.4 Pairs plot showing the MCMC samples and posterior distributions for the transit parameters of AU Mic b as presented in Table 2. Details are the same as in Fig. C.3 for the analysis of the three transits of AU Mic c.
In the text