New constraints on the planetary system around the young active star AU Mic. Two transiting warm Neptunes near mean-motion resonance

AU Mic is a young, active star whose transiting planet was recently detected. We report our analysis of its TESS data, where we modeled the BY Draconis type quasi-periodic rotational modulation by starspots simultaneously to the flaring activity and planetary transits. We measured a flare occurrence rate of 6.35 flares per day for flares with amplitudes in the range of $0.06\%<f_{\rm max}<1.5\%$ of the star flux. We employed a Bayesian MCMC analysis to model the five transits of AU Mic b, improving the constraints on the planetary parameters. The planet radius of $4.07\pm0.17$~R$_{\oplus}$ and a mean density of $1.4\pm0.4$~g~cm$^{-3}$ confirms that it is a Neptune-size moderately inflated planet. While a single feature possibly due to a second planet was previously reported in the former TESS data, we report the detection of two additional transit-like events in the new TESS observations of July 2020. This represents substantial evidence for a second planet (AU Mic c) in the system. We analyzed its three transits and obtained an orbital period of $18.859019\pm0.000016$~d and a planetary radius of $3.24\pm0.16$~R$_{\oplus}$, which defines it as a warm Neptune-size planet with an expected mass in the range of 2.2~M$_{\oplus}$~$<M_{\rm c}<$25.0~M$_{\oplus}$. The two planets in the system are in near 9:4 mean-motion resonance. We show that this configuration is dynamically stable and should produce transit-timing variations (TTV). Our non-detection of significant TTV in AU Mic b suggests an upper limit for the mass of AU Mic c of $<7$~M$_{\oplus}$, indicating that this planet is also likely to be inflated. As a young multi-planet system with at least two transiting planets, AU Mic becomes a key system for the study of atmospheres of infant planets and of planet-planet and planet-disk dynamics at the early stages of planetary evolution.


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. (2020) measured 17.1`4 .7 4.5 M C thanks to infrared observations secured with the SPIRou spectropolarimeter . 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 C and 4.0 R C and a period between 0.5 and 256 d is estimated as 8.4`1 .2 1.1 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 edgeon debris disk and that AU Mic b is a transiting planet with an aligned orbit, the chances that other planets also reside in coplanar 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. 8630 .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 References.
(1) Plavchan et al. (2020); (2) White et al. (2015); (3) Mamajek & Bell (2014); (4) Gaia Collaboration (2018); (5) Claret (2018) 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 multiflare 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.

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 sec compared to the 2 min time sampling of the first visit. We obtained the de-trended photometric time series from the Mikulski Archive for Space Telescopes (MAST) using the MAST astroquery tool 1 . Figure 1 shows the AU Mic PDCSAP 2 flux in units of electrons-persecond (e´s´1) 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.

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, 5, and 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.

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 3 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 minutes. 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´1. Thus, the AU Mic global occurrence rate of flares with amplitude above f max Á 0.06% of the median stellar flux (279544 e´s´1) is estimated as 6.35 flares per day or about 1 flare every 3.8 hours. 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 .0 f max d´1, with f max in units of 10 2 e´s´1. 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.  Plavchan et al. (2020) is about 3.5 hours 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 hours, 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.

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 an-alyze these five transits together, we start 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 limbdarkening 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) 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. 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  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

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 0`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 " Up1, 500q d to explore a broad range of periods, and then with a normal prior of P " Np18.86, 0.01q 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 ‹ " Up1, 500q and then with a normal prior estimated using the Kepler's law, that is, a{R ‹ " Np31.5, 1.4q. The priors for the times of conjunction are measured locally on each event, for example, for the first event, we obtained T 0 " Np1342.23, 0.01q. 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 p {R ‹ " Np0.04, 0.01q. The orbital inclination is considered with an uniform prior distribution of i " Up85˝, 90˝q 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 " Np18.86, 0.01q 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 in Appendix C.

Discussion
We combine the transit parameters of AU Mic b and c presented in Tables 2 and 3 with the star parameters from Table 1 to cal- The bottom right panel shows the 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. culate 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, for f spot being the spot filling factor estimated at 20% based on magnetic activity (Klein et al. 2020) 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-tophotosphere flux fraction integrated over the TESS band pass (between 600 nm 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 eff p " 0.94R obs p , 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 " ? 0.05 to the planet radius measurements, which is included in the uncertainty values reported in Table 4. Klein et al. (2020) recently measured a semi-amplitude of the radial velocity of the reflex motion caused by AU Mic b of K b " 8.5`2 .3 2.2 m s´1, 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´3.
To estimate a plausible mass range for AU Mic c, we considered a population of transiting exoplanets 4 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 4 exoplanet parameters compiled from exoplanet.eu 0.4 g cm´3 ă ρ c ă4.1 g cm´3. 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´1 ă K c ă9.5 m s´1. 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(Laskar , 1993 in the vicinity of the best fit, in the same way as achieved for other planetary systems (e.g., Correia et al. 2005Correia et al. , 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´3 yr and general relativity corrections. References. : semi-major axis derived from the fit period and the Kepler's law; § radius values correspond to the effective radius, which is already corrected for spot coverage and contains an additional uncertainty due to rotational modulation, as explained in Sect. 7; § § Klein et al. (2020) 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 independent measurement 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 UpT c´0 .005, T c`0 .005q, 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 allow us 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 C ) ranging between 0.32 au 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.

Conclusions
We presented an analysis of the photometric TESS observations of the active young M1 star AU Mic, where we model the rota-  (Table 4), with M c " 25.0 M ' (a), and M c " 2.2 M ' (b). The phase space of the system is explored by varying the semimajor axis a c (and thus the period P c ) and eccentricity e c of the outer planet; e c is plotted as a function of P c and not a c , as the uncertainty on a c 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 . The "red" zone corresponds to highly unstable orbits, while the "dark-blue" region can be 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 M c (b), the libration island is mostly stable while it is not the case for the larger value of M c (a).
tional modulation by starspots, the flaring activity, and the planetary transits simultaneously. Our analysis delivered an estimation 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 of T c " 2 458 330.39051˘0.00015Èˆ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  (Table 5).
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 c " 2 458 342.2223˘0.0005`E1 8.859019`0 .000016 0.000015 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 C indicates that AU Mic c is a Neptune-size planet with an expected mass in the range of 2.2 M C ă M c ă25.0 M C , 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 C , implying an upper limit to the planet density of ρ c ă1.3 g cm´3. 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.

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 removed from the original data using the models computed in Sects. 4, 5, and 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 χ 2 " ř pF i´Fc q 2 {σ 2 i , 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 χ 2 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 by Plavchan 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´1, which is more than twice the dispersion of residuals from the fit to the time series, that is, 99 e´s´1 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. 2020), which could generate a dispersion in the phase diagram since starspots at different latitudes are modulated by different periods.
Assuming the radius of AU Mic measured by interferometry of R ‹ " 0.75˘0.03 R d (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´1. This gives a maximum value of the projected velocity of v eq sin pi ‹ q ă 8.1 km s´1, for i ‹ being the inclination angle of the rotation axis of the star with respect to the line of sight.