# Transit Timing Variations for AU Microscopii (catalog ) b & c

Justin M. Wittrock Department of Physics & Astronomy, George Mason University, 4400 University Drive MS 3F3, Fairfax, VA 22030, USA<EMAIL_ADDRESS>Stefan Dreizler Institut für Astrophysik, Georg-August-Universität, Friedrich-Hund-Platz 1, 37077 Göttingen, Germany  Michael A. Reefe Department of Physics & Astronomy, George Mason University, 4400 University Drive MS 3F3, Fairfax, VA 22030, USA  Brett M. Morris Center for Space and Habitability, University of Bern, Gesellschaftsstrasse 6, 3012 Bern, Switzerland  Peter P. Plavchan Department of Physics & Astronomy, George Mason University, 4400 University Drive MS 3F3, Fairfax, VA 22030, USA  Patrick J. Lowrance IPAC-Spitzer, California Institute of Technology, MC 314-6, 1200 E. California Blvd., Pasadena, California 91125, USA  Brice-Olivier Demory Center for Space and Habitability, University of Bern, Gesellschaftsstrasse 6, 3012 Bern, Switzerland  James G. Ingalls IPAC-Spitzer, California Institute of Technology, MC 314-6, 1200 E. California Blvd., Pasadena, California 91125, USA  Emily A. Gilbert Department of Astronomy & Astrophysics, University of Chicago, 5640 S. Ellis Ave, Chicago, IL 60637, USA  University of Maryland, Baltimore County, 1000 Hilltop Circle, Baltimore, MD 21250, USA  The Adler Planetarium, 1300 South Lakeshore Drive, Chicago, IL 60605, USA  NASA Goddard Space Flight Center, 8800 Greenbelt Road, Greenbelt, MD 20771, USA  GSFC Sellers Exoplanet Environments Collaboration  Thomas Barclay University of Maryland, Baltimore County, 1000 Hilltop Circle, Baltimore, MD 21250, USA  NASA Goddard Space Flight Center, 8800 Greenbelt Road, Greenbelt, MD 20771, USA  Bryson L. Cale Department of Physics & Astronomy, George Mason University, 4400 University Drive MS 3F3, Fairfax, VA 22030, USA  Karen A. Collins Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Kevin I. Collins Department of Physics & Astronomy, George Mason University, 4400 University Drive MS 3F3, Fairfax, VA 22030, USA  Ian J. M. Crossfield Department of Physics & Astronomy, University of Kansas, 1251 Wescoe Hall Dr, Lawrence, KS 66045, USA  Diana Dragomir Department of Physics & Astronomy, University of New Mexico, Albuquerque, NM, USA  Jason D. Eastman Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Mohammed El Mufti Department of Physics & Astronomy, George Mason University, 4400 University Drive MS 3F3, Fairfax, VA 22030, USA  Dax Feliz Department of Physics & Astronomy, Vanderbilt University, 6301 Stevenson Center Ln., Nashville, TN 37235, USA  Jonathan Gagné Planétarium Rio Tinto Alcan, Espace pour la vie, 4801 av. Pierre-De Coubertin, Montréal, QC H1V 3V4, Canada  Institute for Research on Exoplanets, Université de Montréal, Département de Physique, C.P. 6128 Succ. Centre-ville, Montréal, QC H3C 3J7, Canada  Eric Gaidos University of Hawai‘i at M$\bar{a}$noa, 1680 East-West Road, Honolulu, HI 96822, USA  Peter Gao Earth and Planets Laboratory, Carnegie Institution for Science, 5241 Broad Branch Rd NW, Washington, DC 20015, USA  Claire S. Geneser Mississippi State University, 75 B. S. Hood Road, Mississippi State, MS 39762, USA  Leslie Hebb Physics Department, Hobart and William Smith Colleges, Geneva, NY 14456, USA  Christopher E. Henze  NASA Ames Research Center, MS 244-30, Moffett Field, CA 94035, USA  Keith D. Horne SUPA Physics and Astronomy, University of St. Andrews, Fife, KY16 9SS Scotland, UK  Jon M. Jenkins NASA Ames Research Center, MS 244-30, Moffett Field, CA 94035, USA  Eric L. N. Jensen Department of Physics & Astronomy, Swarthmore College, Swarthmore PA 19081, USA  Stephen R. Kane University of California, Riverside, 900 University Ave., Riverside, CA 92521, USA  Laurel Kaye Sub-department of Astrophysics, Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford, OX1 3RH, UK  Eder Martioli Sorbonne Université, CNRS, UMR 7095, Institut d’Astrophysique de Paris, 98 bis bd Arago, 75014 Paris, France  Laboratório Nacional de Astrofísica, Rua Estados Unidos 154, Itajubá, MG 37504-364, Brazil  Teresa A. Monsue NASA Goddard Space Flight Center, 8800 Greenbelt Road, Greenbelt, MD 20771, USA  Enric Pallé Instituto de Astrofísica de Canarias, Vía Láctea s/n, E-38205 La Laguna, Tenerife, Spain  Departamento de Astrofísica, Universidad de La Laguna, E-38206 La Laguna, Tenerife, Spain  Elisa V. Quintana NASA Goddard Space Flight Center, 8800 Greenbelt Road, Greenbelt, MD 20771, USA  Don J. Radford Brierfield Observatory, Bowral, New South Wales, Australia  Veronica Roccatagliata Dipartimento di Fisica “Enrico Fermi”, Universita’ di Pisa, Largo Pontecorvo 3, 56127 Pisa, Italy  INFN, Sezione di Pisa, Largo Bruno Pontecorvo 3, 56127 Pisa, Italy  INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy  Joshua E. Schlieder NASA Goddard Space Flight Center, 8800 Greenbelt Road, Greenbelt, MD 20771, USA  Richard P. Schwarz Patashnick Voorheesville Observatory, Voorheesville, NY 12186, USA  Avi Shporer Department of Physics and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA  Keivan G. Stassun Department of Physics & Astronomy, Vanderbilt University, 6301 Stevenson Center Ln., Nashville, TN 37235, USA  Christopher Stockdale Hazelwood Observatory, Hazelwood South, Victoria, Australia  Thiam-Guan Tan Perth Exoplanet Survey Telescope, Perth, Western Australia, Australia  Curtin Institute of Radio Astronomy, Curtin University, Bentley, Western Australia 6102  Angelle M. Tanner Mississippi State University, 75 B. S. Hood Road, Mississippi State, MS 39762, USA  Andrew Vanderburg Department of Physics and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA  Laura D. Vega Department of Astronomy, University of Maryland, College Park, MD 20742, USA  NASA Goddard Space Flight Center, 8800 Greenbelt Road, Greenbelt, MD 20771, USA  Center for Research and Exploration in Space Science & Technology, NASA/GSFC, Greenbelt, MD 20771, USA  Songhu Wang Department of Astronomy, Indiana University, Bloomington, IN 47405, USA

###### Abstract

We explore the transit timing variations (TTVs) of the young (22 Myr) nearby
AU Mic (catalog ) planetary system. For AU Mic (catalog ) b, we introduce
three Spitzer (4.5 $\mu$m) transits, five TESS transits, 11 LCO transits, one
PEST transit, one Brierfield transit, and two transit timing measurements from
Rossiter-McLaughlin observations; for AU Mic (catalog ) c, we introduce three
TESS transits. We present two independent TTV analyses. First, we use
EXOFASTv2 to jointly model the Spitzer and ground-based transits and to obtain
the midpoint transit times. We then construct an O–C diagram and model the
TTVs with Exo-Striker. Second, we reproduce our results with an independent
photodynamical analysis. We recover a TTV mass for AU Mic (catalog ) c of
10.8${}^{+2.3}_{-2.2}$ M⊕. We compare the TTV-derived constraints to a recent
radial-velocity (RV) mass determination. We also observe excess TTVs that do
not appear to be consistent with the dynamical interactions of b and c alone,
and do not appear to be due to spots or flares. Thus, we present a
hypothetical non-transiting “middle-d” candidate exoplanet that is consistent
with the observed TTVs, the candidate RV signal, and would establish the AU
Mic system as a compact resonant multi-planet chain in a 4:6:9 period
commensurability. These results demonstrate that the AU Mic (catalog )
planetary system is dynamically interacting producing detectable TTVs, and the
implied orbital dynamics may inform the formation mechanisms for this young
system. We recommend future RV and TTV observations of AU Mic (catalog ) b and
c to further constrain the masses and to confirm the existence of possible
additional planet(s).

††facilities: Brierfield:0.36 m (Moravian G4-16000 KAF-16803), CFHT (SPIRou),
ExoFOP, Exoplanet Archive, IRSA, LCOGT (SAAO:1 m & SSO:1 m; Sinistro), MAST,
PEST:0.30 m (SBIG ST-8XME), Spitzer (IRAC), TESS, VLT:Antu
(ESPRESSO)††software: AstroImageJ (Collins et al., 2017), astropy (Astropy
Collaboration et al., 2013, 2018), batman (Kreidberg, 2015), bayesflare
(Pitkin et al., 2014), celerite (Foreman-Mackey et al., 2017), celerite2
(Foreman-Mackey et al., 2017; Foreman-Mackey, 2018), emcee (Foreman-Mackey et
al., 2013), EXOFASTv2 (Eastman et al., 2019), exoplanet (Foreman-Mackey et
al., 2021), Exo-Striker (Trifonov, 2019), fleck (Morris, 2020), ipython (Pérez
& Granger, 2007), lightkurve (Lightkurve Collaboration et al., 2018),
matplotlib (Hunter, 2007), mercury6 (Chambers, 1999), numpy (Harris et al.,
2020), rebound (Rein & Liu, 2012; Rein & Spiegel, 2015), scipy (Virtanen et
al., 2020), TAPIR (Jensen, 2013), xoflares (Gilbert et al., 2021)

##  1 Introduction

Exoplanetary sciences have been expanding over the past few decades, with its
fields increasingly diversifying thanks in large part to several successful
and diligent missions, including Kepler (Borucki et al., 2010), K2 (Howell et
al., 2014), and the Transiting Exoplanet Survey Satellite (TESS, Ricker et
al., 2015). TESS has detected 4 548 transiting candidates (TOIs) as of 2021
October 7 and 159 confirmed planets as of 2021 October
4111https://exoplanetarchive.ipac.caltech.edu. Many discoveries have
challenged our theories of planet formation, such as hot Jupiters – e.g., 51
Pegasi b (Mayor & Queloz, 1995), HD 209458 b (Henry et al., 2000), TOI-628 b
(Rodriguez et al., 2021) – planets in highly eccentric orbit – e.g., 16 Cygni
B b (Cochran et al., 1997), BD+63 1405 b (Dalal et al., 2021), HD 26161 b
(Rosenthal et al., 2021) – and compact systems – e.g., HD 108236 (Daylan et
al., 2021; Bonfanti et al., 2021), TOI-178 (Leleu et al., 2021), TRAPPIST-1
(Gillon et al., 2016, 2017a). One way to investigate how these systems form
and evolve can be done by probing young stellar systems, when their
characteristics and orbital dynamics are still undergoing progression. Several
young exoplanet systems have recently been discovered by the TESS and K2
missions – e.g., DS Tucanae A (Newton et al., 2019), K2-25 (Mann et al.,
2016), K2-33 (David et al., 2016), V1298 Tauri (David et al., 2019) – and
other exoplanet detection methods including direct imaging and radial
velocities (RVs) – e.g., HD 47366 (Sato et al., 2016), HR 8799 (Marois et al.,
2008, 2010), PDS 70 (Keppler et al., 2018; Haffert et al., 2019), $\beta$
Pictoris (Lagrange et al., 2009, 2019). Further probing of certain systems
such as PDS 70 revealed a potentially moon-forming circumplanetary disk around
PDS 70 c (Benisty et al., 2021). This new and growing population of transiting
young exoplanets has recently enabled a new frontier in the study of planet
formation and evolution. Among the nearby young exoplanet systems, the nearest
one is AU Microscopii (catalog ) (Tables 1 & 1).

Table 1: Stellar properties for host star AU Mic (catalog ).

Property | Unit | Quantity | Ref
---|---|---|---
Spectral Type | … | M1Ve | …
mV | … | 8.81 $\pm$ 0.10 | …
mTESS | … | 6.755 $\pm$ 0.032 | …
$\alpha_{\rm J2000}$ | h:m:s | 20:45:09.53 | 1
$\delta_{\rm J2000}$ | deg:am:as | -31:20:27.24 | 1
$\mu_{\rm\alpha}$ | mas/yr | 281.424 $\pm$ 0.075 | 1
$\mu_{\rm\delta}$ | mas/yr | -359.895 $\pm$ 0.054 | 1
Distance | pc | 9.7221 $\pm$ 0.0046 | 2
Parallax | mas | 102.8295 $\pm$ 0.0486 | 1
M⋆ | M☉ | 0.50 $\pm$ 0.03 | 3
R⋆ | R☉ | 0.75 $\pm$ 0.03 | 4
Teff | K | 3 700 $\pm$ 100 | 5
L⋆ | L☉ | 0.09 | 5
Age | Myr | 22 $\pm$ 3 | 6
Prot | days | 4.863 $\pm$ 0.010 | 3
v $\sin{i}$ | km/s | 8.7 $\pm$ 0.2 | 7

References. — (1) Gaia Collaboration (2018); (2) Bailer-Jones et al. (2018);
(3) Plavchan et al. (2020); (4) White et al. (2019); (5) Plavchan et al.
(2009); (6) Mamajek & Bell (2014); (7) Lannier et al. (2017)

Table 2: Planetary properties for AU Mic (catalog ) system.

Property | Description | Unit |  AU Mic (catalog ) b |  AU Mic (catalog ) c | Ref
---|---|---|---|---|---
Porb | Orbital Period | days | 8.4630004${}^{+0.0000058}_{-0.0000060}$ | 18.858982${}^{+0.000053}_{-0.000050}$ | 1
a | Semi-Major Axis | au | 0.0645 $\pm$ 0.0013 | 0.1101 $\pm$ 0.0022 | 2
e | Eccentricity | … | 0.12${}^{+0.16}_{-0.08}$ | 0.13${}^{+0.16}_{-0.09}$ | 1
i | Inclination | deg | 89.5 $\pm$ 0.3 | 89.0${}^{+0.5}_{-0.4}$ | 2
$\omega$ | Argument of Periastron | deg | -0.3${}^{+2.4}_{-2.3}$ | -0.3${}^{+2.5}_{-2.2}$ | 1
Mp | Planetary Mass | MJ | 0.054 $\pm$ 0.015 | 0.007 $<$ Mc $<$ 0.079 | 2
M⊕ | 17 $\pm$ 5 | 2 $<$ Mc $<$ 25
Rp | Planetary Radius | RJ | 0.374${}^{+0.021}_{-0.020}$ | 0.249${}^{+0.028}_{-0.027}$ | 1
R⊕ | 4.19${}^{+0.24}_{-0.22}$ | 2.79${}^{+0.31}_{-0.30}$
$\rho_{\rm p}$ | Planetary Density | g/cm3 | 1.4 $\pm$ 0.4 | 0.4 $<$ $\rho_{c}$ $<$ 4.1 | 2
K | RV Semi-Amplitude | m/s | 8.5${}^{+2.3}_{-2.2}$ | 0.8 $<$ Kc $<$ 9.5 | 2
TC $-$ 2 458 000 | Time of Conjunction | BJD | 330.39051 $\pm$ 0.00015 | 342.2223 $\pm$ 0.0005 | 2
tduration | Transit Duration | hours | 3.50 $\pm$ 0.08 | 4.5 $\pm$ 0.8 | 2
Rp/R⋆ | … | … | 0.0512 $\pm$ 0.0020 | 0.0340${}^{+0.0034}_{-0.0033}$ | 1
a/R⋆ | … | … | 19.1 $\pm$ 0.3 | 29 $\pm$ 3 | 2
b | Impact Parameter | … | 0.16${}^{+0.13}_{-0.11}$ | 0.30${}^{+0.21}_{-0.20}$ | 1

References. — (1) Gilbert et al. (2021); (2) Martioli et al. (2021)

AU Microscopii (catalog ) (TOI-2221, TIC 441420236, HD 197481, GJ 803) is a
young ($22\pm 3$ Myr, Mamajek & Bell, 2014), nearby (9.7 pc, Bailer-Jones et
al., 2018) BY Draconis variable star with spectral type M1Ve and relative
brightness mV=8.81 mag. It is known to be fairly active, with numerous flares
having been observed and studied at various wavelengths (Butler et al., 1981;
Kundu et al., 1987; Cully et al., 1993; Tsikoudi & Kellett, 2000; Gilbert et
al., 2021). Kalas et al. (2004) observed the presence of a large dust disk
having a radius between 50 and 210 au from the young star, having first been
detected as a mid-infrared flux excess with IRAS (Fajardo-Acosta et al., 2000;
Zuckerman, 2001; Song et al., 2002; Liu et al., 2004; Plavchan et al., 2005).
Later, Plavchan et al. (2020) discovered a Neptune-sized transiting planet AU
Mic (catalog ) b interior to a spatially-resolved debris disc and with orbital
period of 8.46 days. Recently, Gilbert et al. (2021) and Martioli et al.
(2021) confirmed the existence of another planet AU Mic (catalog ) c with
orbital period of 18.86 days, which put the planets near a 4:9 orbital
commensurability. The aforementioned traits of AU Mic (catalog ) and its
planets make this system a unique, viable laboratory for studying the stellar
activity of a young M dwarf, the planetary formation, the evolution of
exoplanet radii as a function of age, orbital architectures of young giant
planet systems, atmospheric characteristics of young exoplanets, and the
interplay between planets and disks.

One method that serves as a useful tool for probing the exoplanetary systems
is Transit Timing Variations (TTVs). Compared to other detection methods, TTVs
can detect terrestrial-mass planets with greater ease (Holman & Murray, 2005).
The planets that are in orbital resonance with each other can amplify the TTV
signals (Agol et al., 2005), so TTVs can be used to search for and measure the
masses of other planets within a given stellar system (e.g., including
noteworthy systems presented in Mazeh et al., 2013; Becker et al., 2015;
Gillon et al., 2017a; Grimm et al., 2018). Many systems have been
characterized with TTVs, such as HIP 41378 (Bryant et al., 2021), K2-146 (Lam
et al., 2020), TOI-216 (Dawson et al., 2021), TOI-1266 (Demory et al., 2020),
TrES-3 (Mannaday et al., 2020), and many Kepler systems (Lithwick et al.,
2012; Mazeh et al., 2013; Hadden & Lithwick, 2014). Martioli et al. (2021)
searched for the TTVs of AU Mic (catalog ) transits from TESS light curves but
did not identify any significant TTVs. Szabó et al. (2021) did a TTV joint
model with TESS and CHEOPS data; they found AU Mic (catalog ) b’s
$\sim$3.9-minute variation across 33 days and attributed AU Mic (catalog ) c
as the potential source of this perturbation. Gilbert et al. (2021) performed
an independent analysis of AU Mic (catalog ) transits from TESS light curves
and were able to detect the TTVs on the order of $\sim$80 seconds.

For this paper, we examine the TTVs of AU Mic (catalog ) planets by
incorporating additional ground and space observations to our analysis. We
present the TTVs of AU Mic (catalog ) b and c to recover constraints on the
mass for AU Mic (catalog ) c and which indicates the presence of TTV excess
that cannot be accounted by both planets b and c alone. In $\S$2, we list the
light curve data we include for TTV analysis and elaborate on some of the
processes that were involved in data reduction. $\S$3 covers the two critical
developments: joint-modeling both the ground-based photometric and Spitzer
light curves and extracting the midpoint transit times from these sets using
the EXOFASTv2 package (Eastman et al., 2019), and constructing the O–C diagram
using the extracted midpoint times from the observations. Then, as explained
in $\S$4, we model the extracted TTVs using the Exo-Striker package (Trifonov,
2019). Next, we attempt to reproduce our results with an independent and
direct photodynamical analysis as described in $\S$5. Lastly, we discuss the
results in $\S$6 and close this paper in $\S$7.

##  2 Data from Observations

We obtained 23 AU Mic (catalog ) b transits and 3 AU Mic (catalog ) c transits
from three years worth of observations with multiple telescopes and have
included them in the analysis (Table 2). In addition to space-based
observations, with original transit observations from TESS and follow-ups from
Spitzer, we have utilized several ground-based facilities in conducting
follow-ups of AU Mic (catalog ), including Brierfield, LCO SAAO, LCO SSO, &
PEST for photometric observations and CFHT equipped with SPIRou, IRTF equipped
with iSHELL, & VLT equipped with ESPRESSO for Rossiter-McLaughlin (R-M)
observations (Tables 2 & 2). The TESS transits and one of the Spitzer transits
have been previously presented in Plavchan et al. (2020); Gilbert et al.
(2021); Martioli et al. (2021), and the R-M observation in Martioli et al.
(2020); Palle et al. (2020). The following subsections detail each telescope
and the methodology employed upon its respective data sets.

Table 3: List of AU Mic (catalog ) observation data incorporated for TTV
analysis. All ground-based photometric observations listed here were organized
via TESS Follow-up Observing Program (TFOP) Working Group
(WG)aahttps://tess.mit.edu/followup.

Planet | Telescope | Date (UT) | Filter | Exposure | No. of | Obs. Dur. | Transit | Ref
---|---|---|---|---|---|---|---|---
Time (sec) | Images | (min) | Coverage |
b |  Brierfield 0.36 m | 2020-08-13 | I | 16 | 398 | 379 | full | …
b |  CFHT (SPIRou) | 2019-06-17 | 955-2 515 nm | 122.6 | 116 | 302.8 | egress | 1
b |  IRTF (iSHELL) | 2019-06-17 | 2.18-2.47 nm | 120 | 47 | 105.2 | egress | 1
b | LCO SAAO 1.0 m | 2020-05-20 | Pan-STARRS Y | 35 | 99 | 262 | egress | …
2020-05-20 | Pan-STARRS zs | 15 | 333 | 266 | egress
2020-06-06 | Pan-STARRS zs | 15 | 266 | 218 | egress
2020-06-23 | Pan-STARRS zs | 15 | 223 | 183 | egress
2020-09-07 | Pan-STARRS zs | 15 | 211 | 172 | ingress
2020-10-11 | Pan-STARRS zs | 15 | 311 | 266 | ingress
b | LCO SSO 1.0 m | 2020-04-25 | Pan-STARRS Y | 35 | 40 | 104 | egress | …
2020-04-25 | Pan-STARRS zs | 15 | 212 | 172 | egress
2020-08-13 | Pan-STARRS zs | 15 | 379 | 312 | full
2020-09-16 | Pan-STARRS zs | 15 | 408 | 340 | full
2020-10-03 | Pan-STARRS zs | 15 | 248 | 219 | egress
b |  PEST 0.30 m | 2020-07-10 | V | 15 | 1 143 | 556 | full | …
b | Spitzer (IRAC) | 2019-02-10 | 4.5 $\mu$m | 0.08 | 3 020 | 475.7 | full | …
2019-02-27 | 4.5 $\mu$m | 0.08 | 3 377 | 475.7 | egress
2019-09-09 | 4.5 $\mu$m | 0.08 | 6 002 | 990.9 | full
b | TESS bb$\sim$12-hour snippets of the $\sim$27-day duration TESS Cycle 1 and 3 light curves were extracted for our analysis, centered approximately on each transit. | 2018-07-26 | TESS | 120 | 329 | 718.0 | full | 2
2018-08-12 | TESS | 120 | 296 | 708.0 | full
2020-07-10 | TESS | 20 | 2 132 | 719.7 | full
2020-07-19 | TESS | 20 | 2 137 | 719.7 | full
2020-07-27 | TESS | 20 | 2 120 | 719.7 | full
c | TESS bb$\sim$12-hour snippets of the $\sim$27-day duration TESS Cycle 1 and 3 light curves were extracted for our analysis, centered approximately on each transit. | 2018-08-11 | TESS | 120 | 342 | 718.0 | full | 2
2020-07-09 | TESS | 20 | 2 138 | 719.7 | full
2020-07-28 | TESS | 20 | 2 133 | 719.7 | full
b |  VLT (ESPRESSO) | 2019-08-07 | 378.2-788.7 nm | 200 | 88 | 359 | full | 3

References. — (1) Martioli et al. (2020); (2) Gilbert et al. (2021); (3) Palle
et al. (2020)

Table 4: List of facilities utilized for photometric and Rossiter-McLaughlin
follow-up observations of AU Mic (catalog ).

Telescope | Instrument | Location | Aperture | Pixel Scale | Resolution | FOV | Ref
---|---|---|---|---|---|---|---
(m) | (arcsec) | (pixels) | (arcmin)
Brierfield | Moravian 16803 | Bowral, New South Wales | 0.36 | 0.732 | 4 096$\times$4 096 | 50$\times$50 | 1
CFHT | SPIRou | Maunakea, Hawai‘i | 3.58 | … | … | … | 2
IRTF | iSHELL | Maunakea, Hawai‘i | 3.2 | … | … | … | 3
LCO SAAO | Sinistro | Sutherland, South Africa | 1.0 | 0.389 | 4 096$\times$4 096 | 26.5$\times$26.5 | 4
LCO SSO | Sinistro | Mount Woorut, New South Wales | 1.0 | 0.389 | 4 096$\times$4 096 | 26.5$\times$26.5 | 4
PEST | SBIG ST-8XME | Perth, Western Australia | 0.3048 | 1.23 | 1 530$\times$1 020 | 31$\times$21 | 5
Spitzer | IRAC | … | 0.85 | 1.22 | 256$\times$256 | 5.2$\times$5.2 | 6
VLT | ESPRESSO | Cerro Paranal, Chile | 8.2 | … | … | … | 7

References. — (1) https://www.brierfieldobservatory.com; (2)
https://www.cfht.hawaii.edu; (3) http://irtfweb.ifa.hawaii.edu; (4)
https://lco.global/observatory; (5) http://pestobservatory.com; (6)
https://www.spitzer.caltech.edu; (7) https://www.eso.org/public/teles-
instr/paranal-observatory/vlt

Table 5: Specifications of instruments used for ground-based Rossiter-
McLaughlin follow-up observations of AU Mic (catalog ).

Instrument | Telescope | Observing Mode |  $\lambda$ Range | Resolving | Aperture | Average | Ref
---|---|---|---|---|---|---|---
(nm) | Power | (arcsec) | SNR
ESPRESSO | VLT | HR (1-UT) | 378.2-788.7 | 140 000 | 1.0 | 93.9 | 1
iSHELL | IRTF | Kgas | 2.18-2.47 | 75 000 | 0.125 | 65 | 2, 4
SPIRou | CFHT | Stokes $V$ Spectropolarimetric | 955-2 515 | 70 000 | 1.29 | 242 | 3, 4

References. — (1) Donati et al. (2020), Palle et al. (2020); (2) Rayner et al.
(2016); (3) Pepe et al. (2021); (4) Martioli et al. (2020)

###  2.1 TESS Photometry

TESS 222https://tess.mit.edu
https://heasarc.gsfc.nasa.gov/docs/tess is a space-based telescope designed to
scan nearby bright F5-M5 stars for transiting exoplanets (Ricker et al.,
2015). Since its launch on 2018 April 18 and the start of its primary mission
on 2018 July 25, TESS has been probing the sky for $\sim$3 years as of this
writing and has made numerous groundbreaking contributions to planetary
detection – e.g., DS Tuc A (Newton et al., 2019), TOI-700 (Gilbert et al.,
2020), TOI-1338 (Kostov et al., 2020). Its two-year primary mission divided
the sky into the Southern and the Northern Ecliptic Hemispheres, with each
being divided further into 13 sectors. TESS began its search in the Southern
Ecliptic Hemisphere and probed each sector for 28 days. Within each 28-day
span, a subset of primary exoplanet transit search target stars in a given
sector were monitored at 2-minute cadence, and Full Frame Images (FFIs) were
collected at 30-minute cadence. The data collected by TESS are then processed
by the Science Processing Operations Center (SPOC), which functions to
generate the calibrated images, perform aperture photometry, remove systematic
artifacts, and searches the light curves for transiting planet signatures
(Jenkins et al., 2016). TESS successfully completed its two-year primary
mission and is now in its extended mission by repeating its observation in
each of 26 sectors, with some notable differences: TESS is probing or will
probe new targets along with the old targets, the 2-minute cadence for 20 000
targets per sector is boosted with 20-second cadence for 1 000 targets per
sector, and FFIs are retrieved at a shorter 10-minute cadence.

Figure 1: TESS photometry of AU Mic (catalog ). The blue points are TESS
measurements. The top plot is from Cycle 1 (Sector 1, 2018 July 25 19:00:27 UT
$-$ 2018 August 22 16:14:51 UT), 2-minute cadence. The bottom plot is from
Cycle 3 (Sector 27, 2020 July 05 18:31:16 UT $-$ 2020 July 30 03:21:15 UT),
20-second cadence.

TESS observed AU Mic (catalog ) (Figure 1) at 2-minute cadence during Cycle 1
(Sector 1, 2018 July 25 19:00:27 UT to 2018 August 22 16:14:51 UT)333The
following Guest Investigator (GI) proposals were awarded for AU Mic (catalog
)’s Cycle 1 observations: G011176/PI Czekala, G011185/PI Davenport, G011264/PI
Davenport, G011180/PI Dressing, G011239/PI Kowalski, G011175/PI Mann, &
G011266/PI Schlieder. using Camera 1 CCD 4. This set has a 1.13-day gap due to
data downlink. During its observation between 2018 August 16 16:00 UT and 2018
August 18 16:00 UT, the Fine Pointing mode calibration was not configured
optimally, culminating in poorer data quality due to excessive spacecraft
pointing jitter. Gilbert et al. (2021) employed the data quality flags to
filter out the problematic part of the data set, resulting in some additional
gaps in data. TESS observed AU Mic (catalog ) again during Cycle 3 (Sector 27,
2020 July 05 18:31:16 UT to 2020 July 30 03:21:15 UT)444The following GI
proposals were awarded for AU Mic (catalog )’s Cycle 3 observations: G03272/PI
Burt, G03227/PI Davenport, G03063/PI Llama, G03228/PI Million, G03205/PI
Monsue, G03141/PI Newton, G03202/PI Paudel, G03263/PI Plavchan, G03226/PI
Silverstein, & G03273/PI Vega., this time at both 20-second and 2-minute
cadences, with the latter constructed by co-added six 20-second exposures. The
data downlink during this period led to a 1.02-day gap in the data.

We use the AU Mic (catalog ) TESS Cycle 1 and 3 transit light curves from
Gilbert et al. (2021) for our primary TTV analysis ($\S$4), and herein we
summarize their analysis. In our independent photodynamical analysis ($\S$5),
we reanalyze the TESS light curves directly. The AU Mic (catalog ) TESS Cycle
1 and 3 light curves were retrieved from the Mikulski Archive for Space
Telescopes (MAST)555https://archive.stsci.edu archive using the lightkurve
package (Lightkurve Collaboration et al., 2018) while setting its bitmask
filter to “default”. The Presearch Data Conditioning (PDC SAP) light curves
were chosen since they addressed crowding and instrumental systematics (Smith
et al., 2012; Stumpe et al., 2012, 2014). After filtering the NaNs out of the
data sets, 25.07 days of Cycle 1, 23.29 days of 2-minute Cycle 3, and 22.57
days of 20-second Cycle 3 are left. Next, the Savitzky-Golay filter is applied
to the light curves to smooth out AU Mic (catalog )’s spot modulation. Since
the flares are abundant in TESS light curves, especially during most of AU Mic
(catalog ) b’s and c’s transits, the bayesflare (Pitkin et al., 2014) and
xoflares packages were used to extract and model the flares instead of
trimming the flares out as in Plavchan et al. (2020). Then, the celerite2
package (Foreman-Mackey et al., 2017; Foreman-Mackey, 2018) was applied to
model the stellar variability of AU Mic (catalog ), and the exoplanet package
(Foreman-Mackey et al., 2021) was utilized to model the transits of AU Mic
(catalog ) b and c.

###  2.2 Spitzer (IRAC) Photometry

The Spitzer Space Telescope666https://www.spitzer.caltech.edu was constructed
as part of NASA’s Great Observatories Program’s final mission to probe various
astrophysical objects at infrared wavelengths (Werner et al., 2004). Launched
in 2003 August 25, Spitzer carried out its primary mission along with NASA’s
Astronomical Search for Origins Program for 5.75 years until the liquid helium
coolant was depleted on 2009 May 15; afterward, it continued under several
extended missions – starting with the Spitzer Warm mission, then the Spitzer
Beyond mission, and finally the Spitzer Final Voyage mission – for the next
10.5 years starting from 2009 July 27 until its decommission on 2020 January
30. Spitzer fulfilled an indispensable role in characterizing exoplanets
(Deming & Knutson, 2020) – e.g. including some benchmark systems observed, HD
189733 b (Grillmair et al., 2007; Todorov et al., 2014), HD 209458 b (Zellem
et al., 2014), HD 219134 b & c (Gillon et al., 2017b), TRAPPIST-1 b through h
(Morris et al., 2018; Zhang et al., 2018; Ducrot et al., 2018, 2020), WASP-26
b (Mahtani et al., 2013).

To collect more data on the planetary object detected by TESS, Spitzer
Director’s Discretionary Time (DDT) observations were proposed and awarded
(PID 14214 & 14241) in the final year of operations for observations of AU Mic
(catalog ) with the Infrared Array Camera (IRAC, Fazio et al., 2004) due to
possible calculated transits on 2019 February 10 & 27 and September 10, which
are presented in this paper. Spitzer observed AU Mic (catalog ) with IRAC on
five occasions during the Spitzer Beyond and Final Voyage missions (2019
February 10 10:58:58 to 18:54:38 UT, 2019 February 27 09:46:23 to 17:42:02 UT,
2019 September 9 13:50:57 to 21:46:45 UT, 2019 September 9 23:26:02 UT to 2019
September 10 06:21:54 UT, & 2019 September 14 03:40:07 to 12:33:36 UT). The
first two observations were originally considered to be eclipses from the
initially assumed orbital period for AU Mic (catalog ) b from the TESS mission
Cycle 1 observations. However, these observations detected additional transits
of AU Mic (catalog ) b, establishing a true period to be half as long as
originally thought (Plavchan et al., 2020). The third observation is of a
transit search for an originally estimated but incorrect period for AU Mic
(catalog ) c, and the fourth observation is of a third transit of AU Mic
(catalog ) b; these two observations have been combined into one light curve
set for this analysis. The final observation is a secondary eclipse search of
AU Mic (catalog ) b which will be described and analyzed in a separate paper
and is not included in this work.

All of the observations were taken using the 32$\times$32 pixel sub-array mode
with an exposure time of 0.08 seconds to avoid saturation on the star
(measurement cadence is 0.1 seconds). After placing the star on the “sweet-
spot” pixel, using the pointing calibration and reference sensor (PCRS) peak-
up mode (Ingalls et al., 2012), we exposed with continuous staring (no
dithers). The observations were all taken at 4.5 $\mu$m, as this channel has a
lower systematics due to the intra-pixel sensitivity. The coordinates were
adjusted for the high parallax and proper motion of AU Mic (catalog ) for the
proposed observation dates. Each observation set consisted of a 30-minute pre-
stare dither pattern, an 8-hour stare, and a 10-minute post-staring dither
pattern. All data were calibrated by the Spitzer pipeline S19.2 and can be
accessed using the Spitzer Heritage Archive
(SHA)777http://sha.ipac.caltech.edu.

For each of the three transit observations discussed here, the following data
reduction steps were performed on each AU Mic (catalog ) stellar image
measured on 0.1-second intervals. We used the IDL routine
box_centroider888https://irsa.ipac.caltech.edu/data/SPITZER/docs/irac/calibrationfiles/pixelphase/box_centroider.pro
supplied by the Spitzer Science Center to measure the location of AU Mic
(catalog ) on the pixel. We then performed aperture photometry on each image
using the IDL Astronomy Users Library routine
aper999https://idlastro.gsfc.nasa.gov/ftp/pro/idlphot/aper.pro with a fixed
aperture of 2.25 pixels and subtracting a sky annulus of 3-7 pixels about the
centroid.

All IRAC photometry at 4.5 $\mu$m contain instrumental systematics caused by
the coupling between spacecraft pointing fluctuations and drifts with intra-
pixel sensitivity variations. For the three transit observations analyzed
here, we take three approaches for detrending the instrument systematics and
compare the results. First, as in Plavchan et al. (2020), we detrended this
effect using an independent pixel mapping dataset measured for non-variable
star BD+67 1044 (Ingalls et al., 2018). Because this calibration star doesn’t
intrinsically vary, we take its photometric variations to reflect the pixel
sensitivity map. We estimated the relative pixel sensitivity at the $(x,y)$
centroid locations of each AU Mic (catalog ) observation using the K-Nearest
Neighbors with Kernel Regression technique described by Ingalls et al. (2018)
and divided all AU Mic (catalog ) measurements by the sensitivities. This
approach was published for the first Spitzer transit in Plavchan et al.
(2020). However, we noticed that additional high-frequency (shorter than the
transit duration) photometric variability remained in the light curve that
looked like astrophysical “hot spot” crossings. But we subsequently identified
a strong correlation of these light curve features with the Spitzer PSF FWHM.
Therefore, second, we detrend the light curve using both the trend time series
for the pixel centroid motion and the PSF FWHM; this is the systematic-
detrended time series that we adopt in this work for further analysis. Third,
we also tested the noise-pixel technique detailed by Lewis et al. (2013) and
achieve qualitatively similar results to our second approach.

We modeled the IRAC intra-pixel sensitivity (Ingalls et al., 2016) using a
modified implementation of the BLISS (BiLinearly-Interpolated Sub-pixel
Sensitivity) mapping algorithm (Stevenson et al., 2012). We used a modified
version of the BLISS mapping (BM) approach to mitigate the correlated noise
associated with intra-pixel sensitivity. In our photometric baseline model, we
complement the BM correction with a linear function of the Point Response
Function (PRF) Full Width at Half-Maximum (FWHM). In addition to the BM, our
baseline model includes the PRF’s FWHM along the $x$ and $y$ axes, which
significantly reduces the level of correlated noise as shown in previous
studies (e.g., Lanotte et al. 2014, Demory et al. 2016a, Demory et al. 2016b,
Gillon et al. 2017a, Mendonça et al. 2018). Our baseline model does not
include time-dependent parameters. Our implementation of this baseline model
is included in a Markov Chain Monte Carlo (MCMC) framework already presented
in the literature (Gillon et al., 2012). We run two chains of 200 000 steps
each and check for convergence and efficient mixing using the Gelman-Rubin
statistic (Gelman & Rubin, 1992); all of the chains have converged with their
GR statistic $<$ 1.01.

We next construct a 2nd-order polynomial model fit to account for the
rotational modulation of stellar activity present in the three Spitzer light
curves. AU Mic (catalog ) is active with its rotational modulation of stellar
activity, and on the timescale of a transit duration can be described by a
2nd-order polynomial (Addison et al., 2020); longer timescales would
necessitate a Gaussian Process or similar analysis as undertaken in Plavchan
et al. (2020) and Gilbert et al. (2021) for the TESS transits. These
polynomial coefficients are marginalized over in our TTV analysis to account
for the timing uncertainties introduced from the rotational modulation of
stellar activity. We cross-check our approach to that using a Gaussian Process
model in Plavchan et al. (2020) and derive consistent TTVs and corresponding
uncertainties.

The second Spitzer transit also has an unusual “jump” feature during the
middle of the transit that was thought to be caused by either a flare or a
transit egress of another planet; we do not identify any systematic indicator
that this jump coincides with. We explored the timing of AU Mic (catalog ) c’s
transits but found that none line up with the Spitzer’s transit of AU Mic
(catalog ) b. So instead we constructed and fit a custom flare model for this
feature, consisting of a linear rapid rise followed by an exponential decay.
The amplitude of the flare model is marginalized over in our analysis of the
TTVs to account for the impact it has on the transit times; however, the flare
rise, peak, and decay times are fixed in our analysis. The adopted flare rise
and decay times are informed by and consistent with the characterization of
the flares in the TESS light curve analyzed in Gilbert et al. (2021). Here,
the 2nd-order polynomial coefficients are degenerate with the flare times and
are marginalized over to account for the impact the flare has in our derived
transit time and uncertainty. Since the flare did not occur during ingress or
egress, it has minimal impact on our derived transit midpoint time and
corresponding uncertainty.

The third Spitzer light curve was additionally detrended with an ad hoc
Gaussian model given the presence of a low-level Gaussian-like trend
coincident with the transit midpoint time (note, not a Gaussian Process, but a
Gaussian change in brightness with time). Again, we do not identify a
systematic indicator correlated with this brightness variation in the light
curve and associate it with an astrophysical origin for AU Mic (catalog ). The
amplitude of the Gaussian model is marginalized over in our analysis to assess
its impact on the TTVs, but the width and peak time were fixed; in this case,
marginalizing over the 2nd-order polynomial coefficients in our model again
compensates for and is degenerate with any error in the fixed Gaussian
centroid time and width. The astrophysical origins of this brightness change
and its coincidence with the transit midpoint time, as well as other remaining
residuals present in the AU Mic (catalog ) Spitzer light curves as seen in
Figure 2, are beyond the scope of this work and are the subject of a future
publication.

For all but one detrending time series, the additive coefficients were set to
zero and the multiplicative coefficients were set to one. The exception is the
flare detrending from the second Spitzer set, with the additive coefficient
set to one and the multiplicative coefficient set to zero. Afterward, we did a
joint model of the Spitzer and ground-based photometric transits using the
EXOFASTv2 package (Eastman et al., 2019). This process to extract the midpoint
transit times from the Spitzer light curves is explained in more details in
$\S$3.

###  2.3 Rossiter-McLaughlin Spectroscopy

The Rossiter-McLaughlin (R-M) technique is an advantageous tool in detecting
and characterizing transiting exoplanets, including determining their spin-
orbit alignments (Ohta et al., 2005; Winn, 2007; Triaud, 2018). The R-M effect
is observed when the planet crosses the host star during the RV observation,
blocking a portion of the star’s rotational signal and generating a
characteristic feature on the time series RV profile (Holt, 1893; Rossiter,
1924; McLaughlin, 1924; Ohta et al., 2005; Winn, 2007). Many exoplanets have
been characterized with R-M’s – e.g. including some benchmark systems, CoRoT-3
b & HD 189733 b (Triaud et al., 2009), KELT-20 b (Rainer et al., 2021), K2-232
b (Wang et al., 2021), TOI-1431 b (Stangret et al., 2021), WASP-17 b (Anderson
et al., 2010).

R-M observation also offer an additional way in which to derive transit mid-
point times independent of photometric observations, a method that has not
previously been commonly used for TTV analysis because it is resource-
intensive with its use of high-resolution spectrometers on large aperture
telescopes. Today, however, TESS mission candidates are relatively brighter
and nearby compared to Kepler systems, and more amenable to R-M observations.
We include transit midpoint times derived from two R-M observations of AU Mic
(catalog ) b’s transits that were obtained at relatively important epochs
shortly after the Spitzer observations and between the two-year gap in the
TESS observations. The first was collected using the SPIRou and iSHELL
instruments, and the second the ESPRESSO instrument. We retrieved the two
transit midpoint times from Martioli (2020, private communication) and Pallé
(2020, private communication), respectively. The following sections summarize
the work done by Martioli et al. (2020) and Palle et al. (2020) on processing
the respective SPIRou \+ iSHELL and ESPRESSO data.

####  2.3.1 CFHT (SPIRou) & IRTF (iSHELL) Spectroscopy

The SpectroPolarimètre InfraRouge
(SPIRou)101010https://www.cfht.hawaii.edu/Instruments/SPIRou, mounted on the
3.6 m Canada-France-Hawai‘i Telescope (CFHT) located atop Maunakea, Hawai‘i,
is a high-resolution near-infrared (NIR) spectrometer that is capable of
imaging in YJHK bands (0.95-2.5 $\mu$m) with a resolving power of $\sim$70 000
and is equipped with a fiber-fed cryogenic high-resolution échelle
spectrograph that can perform high-precision velocimetry and
spectropolarimetry, which allows it to simultaneously observe magnetic
features and stellar activities of the host stars (Donati et al., 2020). The
iSHELL 111111http://irtfweb.ifa.hawaii.edu/~ishell, installed on the 3.2 m
NASA Infrared Telescope Facility (IRTF) also located atop Maunakea, Hawai‘i,
is a high-resolution 1.1-5.3 $\mu$m spectrometer with a resolving power of 75
000 and was designed to replace CSHELL as an instrument with enhanced
spectroscopic capabilities (Rayner et al., 2016). These qualities make SPIRou
and iSHELL useful tools to carry out follow-up observations on transiting
exoplanets of young, active M dwarfs (Morin et al., 2010; Afram & Berdyugina,
2019), such as AU Mic (catalog ).

As part of the SPIRou Legacy Survey’s Work Package 2 (WP2) (Donati et al.,
2020), Martioli et al. (2020) observed AU Mic (catalog ) on 2019 June 17
10:10:56 to 15:13:45 UT with SPIRou set in Stokes V spectropolarimetric mode
and captured an egress of AU Mic (catalog ) b. 116 spectra were collected that
night, each taken at 122.6 seconds exposure, with the average SNR of 242.
Martioli et al. (2020) also used iSHELL set at K${}_{\mbox{\scriptsize gas}}$
mode for simultaneous but shorter observation of AU Mic (catalog ) (2019 June
17 11:08:19 to 12:53:32 UT); 47 120-second spectra were collected with that
instrument, with SNR of $\sim$60-70. The typical seeing condition was 0.96
$\pm$ 0.13”, the initial and final airmass were 2.9 and 1.8, respectively,
with its minimum being 1.59, and the Moon was 99% illuminated and 40.3∘ from
the target.

Martioli et al. (2020) implemented the reduction pipeline APERO (A PipelinE to
Reduce Observations, Cook et al. in prep.) to reduce and process the SPIRou
data and to calculate the cross-correlation functions (CCFs). Next, the
“M2_weighted_RV_-5.mas” line mask was applied to the spectra, and the lines
were masked if their telluric absorption is deeper than 40%. The line mask
then underwent further refinement by removing lines not present in AU Mic
(catalog )’s Stokes-I spectrum using the technique from Moutou et al. (2020).
The CCFs were calculated from each spectral order and summed to achieve
greater precision, then the RVs were measured from the CCF by using the
velocity shift’s least-square fit. The iSHELL RVs were extracted using the
pychell pipeline, and yielded consistent RVs and precision with the SPIRou RVs
(Cale et al., 2019).

Next, Martioli et al. (2020) constructed the R-M model using the emcee Markov
chain Monte Carlo (MCMC) package (Foreman-Mackey et al., 2013) and the stellar
activity model using the approach from Donati et al. (1997). The stellar
activity model was then subtracted from the measured RVs, and the R-M model
was then applied as a correction to the subtracted RVs. Finally, Martioli
(2020, private communication) modeled AU Mic (catalog ) b’s midpoint time from
SPIRou \+ iSHELL’s best-fit subtracted and corrected RV model using TC =
$\mathcal{N}$(2458330.39153, 0.00070) as a prior.

####  2.3.2 VLT (ESPRESSO) Spectroscopy

Echelle Spectrograph for Rocky Exoplanet and Stable Spectroscopic Observations
(ESPRESSO)121212https://www.eso.org/sci/facilities/paranal/instruments/espresso.html
is a high-precision RV spectrometer situated in the Combined-Coudé Laboratory
(CCL) at the focus of the Very Large Telescope (VLT) atop Cerro Paranal, Chile
(Pepe et al., 2021). Its spectrograph probes the sky at 378.2-788.7 nm range,
and it can use either four 8.2 m telescopes (4-UT) with lower resolution of
$\sim$70 000 or only one of them (1-UT) with higher resolutions of $\sim$140
000 in the High-Resolution (HR) mode or $>$190 000 in the Ultra-High-
Resolution (UHR) mode.

Palle et al. (2020) observed AU Mic (catalog ) on 2019 August 7 3:24 to 9:23
UT with ESPRESSO set at the standard HR (1-UT) mode and captured a full
transit of AU Mic (catalog ) b. 88 spectra were collected that night, each
taken at 200 seconds exposure. The SNR averaged around 93.9, the initial and
final airmass were 1.03 and 2.37, respectively, with its minimum being 1.007,
and the sky was clear.

Palle et al. (2020) applied several separate approaches in reducing and
modeling the ESPRESSO data; however, we only highlight one of those approaches
that provided us the midpoint time for this paper. The SERVAL package
(Zechmeister et al., 2018) was implemented to extract and calibrate the
spectra and generate the RV profile of AU Mic (catalog ) b. Next, the R-M
effect was modeled using the combination of the celerite (Foreman-Mackey et
al., 2017) package’s Gaussian process (GP) and PyAstronomy package (Czesla et
al., 2019), and the stellar activity was modeled with GP described by a Matérn
3/2 kernel implemented by celerite. These models are then applied as
corrections to the ESPRESSO RV profile. Pallé (2020, private communication)
extracted the midpoint time from the SERVAL GP + PyAstronomy best-fit RV
profile using TC = $\mathcal{N}$(2458702.77277, 0.00189) as a prior.

###  2.4 Ground-Based Photometry

The TESS Follow-up Observing Program (TFOP) Working Group
(WG)131313https://tess.mit.edu/followup coordinated numerous ground-based
follow-ups for various TOIs, including AU Mic (catalog ). As a result, 13 AU
Mic (catalog ) follow-up photometric transit observations were made using
different observatories: one Brierfield 0.36 m, six LCO SAAO 1.0 m, five LCO
SSO 1.0 m, and one PEST 0.30 m. The light curves from these observations are
available online through ExoFOP-TESS141414https://exofop.ipac.caltech.edu/tess
(Akeson et al., 2013). The follow-up observation schedules were conducted with
the online version of the TAPIR package (Jensen, 2013). We utilized the
AstroImageJ (AIJ, Collins et al., 2017) to process the ground-based light
curves (except PEST, which was processed through their own pipeline) and then
create a subset table containing only BJD_TDB, normalized detrended flux, flux
uncertainty, and detrending columns from the ground-based light curves to
prepare them for EXOFASTv2 modeling and extraction of midpoint times ($\S$3).
The following sub-subsections describe the role each telescope played in
collecting and processing the light curves.

####  2.4.1 LCOGT (Sinistro) Photometry

We made use of two 1.0 m LCO Ritchey-Chretien Cassegrain telescopes, both
equipped with Sinistro, that are part of the Las Cumbres Observatory Global
Telescope network (LCOGT, Brown et al.,
2013)151515https://lco.global/observatory. Two of Sinistro’s filters used for
AU Mic (catalog ) observations were Pan-STARRS Y and Pan-STARRS zs, with their
central wavelength peaks at 1004.0 and 870.0 nm, respectively. The third
filter B was used for simultaneous observation with Y; however, the data
collected with B are omitted from this paper due to non-detection of AU Mic
(catalog ) b’s transits and more pronounced stellar activity in the bluer
band.

LCO SSO, located on Mount Woorut near Coonabarabran, New South Wales,
Australia, observed AU Mic (catalog ) on four separate nights (2020 April 25
16:03:03 to 17:49:38 UT at 35 seconds exposure with Y and 2020 April 25
16:00:47 to 18:52:43 UT, 2020 August 13 12:41:53 to 17:50:19 UT, 2020
September 16 09:12:31 to 14:41:47 UT, & 2020 October 3 09:17:36 to 12:36:33 UT
at 15 seconds exposure with zs). 40 images from the first night were collected
with Y, and 212, 379, 408, & 248 images from the respective first, second,
third, & fourth nights were collected with zs. An egress was captured on the
first and fourth nights while a full transit was captured on the second and
third nights.

LCO SAAO, located in Sutherland, South Africa, observed AU Mic (catalog ) on
five separate nights (2020 May 20 22:53:29 UT to 2020 May 21 03:15:29 UT at 35
seconds exposure with Y and 2020 May 20 22:51:11 UT to 2020 May 21 03:17:38
UT, 2020 June 6 21:44:20 UT to 2020 June 7 01:22:36 UT, 2020 June 23 20:37:29
to 23:36:35 UT, 2020 September 7 21:55:46 UT to 2020 September 8 00:43:57 UT,
& 2020 October 11 18:08:39 to 22:29:54 UT at 15 seconds exposure with zs). 99
images from the first night were collected with Y, and 333, 266, 223, 211, &
311 images from the respective first, second, third, fourth, & fifth nights
were collected with zs. The second night’s photometric quality was impacted by
a combination of clouds and a full Moon. An egress was captured on the first
through third nights while an ingress was captured on the fourth and fifth
nights.

All light curves from LCOGT were reduced and detrended with AIJ. For each
LCOGT night, the following detrending parameters were applied: AIRMASS for
UT2020-04-25 (Y), UT2020-05-20 (Y & zs), UT2020-06-06, & UT2020-06-23;
Width_T1 for UT2020-04-25 (zs), UT2020-09-16, UT2020-10-03, & UT2020-10-11;
AIRMASS + Width_T1 for UT2020-08-13; and Width_T1 + Sky/Pixel_T1 for
UT2020-09-07. We also used AIJ to generate a subset table for each light
curve.

####  2.4.2 PEST (SBIG ST-8XME) Photometry

Perth Exoplanet Survey Telescope (PEST)161616http://pestobservatory.com, based
in Perth, Western Australia, is a 12” (0.3048 m) Meade LX200 Schmidt-
Cassegrain Telescope that was equipped with SBIG ST-8XME camera at the time of
the AU Mic (catalog ) observation. PEST observed AU Mic (catalog ) on 2020
July 10 13:26:53 to 22:42:33 UT with V and captured a full transit. 1143
images were collected, each at 15 seconds exposure. The PEST light curve was
reduced and processed through the PEST
pipeline171717http://pestobservatory.com/the-pest-pipeline. We then use AIJ to
create a subset table that included PEST-generated detrending parameters
comp_flux + x_coord + y_coord + dist_center + fwhm + airmass + sky.

####  2.4.3 Brierfield (Moravian 16803) Photometry

The Brierfield Observatory181818https://www.brierfieldobservatory.com, located
in Bowral, New South Wales, Australia, houses the 14” (0.36 m) Planewave
Corrected Dall-Kirkham Astrograph telescope mounted with the instrument
Moravian G4-16000 KAF-16803. Brierfield observed AU Mic (catalog ) on 2020
August 13 11:35:21 to 17:54:35 UT with I and captured a full transit. 398
images were collected, each at 16 seconds exposure. The Brierfield light curve
was reduced and detrended with AIJ; the detrending parameters were
Meridian_Flip + X(FITS)_T1 + Y(FITS)_T1 + tot_C_cnts. Afterward, a subset
table was generated from this light curve with AIJ.

##  3 TESS, Spitzer, & Ground-Based Photometric Joint Modeling

We use the EXOFASTv2 package (Eastman et al., 2019) to model the transits and
characterize our light curves. EXOFASTv2 estimates the posterior probabilities
through the Markov chain Monte Carlo (MCMC) to both determine the statistical
significance of our ground-based & Spitzer detections and the confidence in
the time of conjunction measurements to assess for the presence of detectable
TTVs. Five TESS transits, three Spitzer transits and 13 ground-based
photometric transits of AU Mic (catalog ) b and three TESS transits of AU Mic
(catalog ) c are included in the model. The following detrending parameters
are treated as additive: flare (Spitzer), sky (Spitzer & PEST), & Sky/Pixel_T1
(LCO SAAO); the remaining detrending parameters are treated as multiplicative.

Table 6: Stellar, planetary, and transit priors for EXOFASTv2 modeling with AU
Mic (catalog )’s eight TESS, three Spitzer, & 13 ground transits.

Prior | Unit | Input
---|---|---
AU Mic (catalog ) b |  AU Mic (catalog ) c
log${}_{10}\left(\frac{{\rm M}_{\rm\star}}{{\rm M}_{\rm\sun}}\right)$ | … |  $\mathcal{N}$(-0.301, 0.026)
R⋆ | R☉ |  $\mathcal{N}$(0.75, 0.03)
Teff | K |  $\mathcal{N}$(3700, 100)
Age | Gyr |  $\mathcal{N}$(0.022, 0.003)
Parallax | mas |  $\mathcal{N}$(102.8295, 0.0486)
TC | BJD_TDB |  $\mathcal{N}$(2458330.39051) |  $\mathcal{N}$(2458342.2223)
log${}_{10}\left(\frac{\rm Period}{\rm days}\right)$ | … |  $\mathcal{N}$(0.92752436) |  $\mathcal{N}$(1.2755182)
Rp/R⋆ | … |  $\mathcal{N}$(0.0512, 0.0020) |  $\mathcal{N}$(0.0340, 0.0034)
e | … |  $\mathcal{N}$(0.12, 0.16) |  $\mathcal{N}$(0.13, 0.16)
TTV Offset | days |  $\mathcal{U}$(-0.02, 0.02)
Depth Offset | … |  $\mathcal{U}$(-0.01, 0.01)

, while the logarithmic functions were calculated. TTV and depth offsets are
arbitrary and applied as constraints to all transits.

Note. — $\mathcal{N}$ denotes the Gaussian priors, and $\mathcal{U}$ denotes
the uniform priors. The Gaussian priors were taken from Tables 1 & 1

Table 7: Apparent magnitude priors for EXOFASTv2’s spectral energy
distribution (SED) fitting of AU Mic (catalog ). This is intended to place
constraints on AU Mic (catalog )’s MIST evolutionary models (Choi et al.,
2016; Dotter, 2016) and is applied only to the first of the two EXOFASTv2 runs

.  Band Apparent Magnitude Ref Gaia 7.84 $\pm$ 0.02 1 GaiaBP 8.94 $\pm$ 0.02 1
GaiaRP 6.81 $\pm$ 0.02 1 J2M 5.44 $\pm$ 0.02 2 H2M 4.83 $\pm$ 0.02 2 K2M 4.53
$\pm$ 0.02 3 B 10.06 $\pm$ 0.02 … V 8.89 $\pm$ 0.18 … gSDSS 9.58 $\pm$ 0.05 4
rSDSS 8.64 $\pm$ 0.09 4 iSDSS 7.36 $\pm$ 0.14 4

References. — (1) Gaia Collaboration (2018); (2) Cutri et al. (2003); (3)
Stauffer et al. (2010); (4) Zacharias et al. (2012)

The Gaussian priors from Table 6 were taken from Tables 1 and 1, while the
logarithmic functions were calculated; the logarithmic version of stellar mass
and orbital period were used because they are the fitted priors in EXOFASTv2.
The TTV and depth offset priors were implemented to place constraints on the
variation of transit timing and depth of all light curves; any transit depth
variability is not investigated further herein. Since both Pan-STARRS Y and
Pan-STARRS zs are not available among the filters in EXOFASTv2, y and z’
(Sloan z) were used as respective approximate substitutes.

Table 8: Detrending parameters incorporated into EXOFASTv2 modeling of AU Mic
(catalog ) b transits. The flare (Spitzer), sky (Spitzer & PEST), &
Sky/Pixel_T1 (LCO SAAO) were implemented as additives; the remaining
detrending parameters were implemented as multiplicative. See $\S$2 for
details on detrending parameters used for each observation. Since both Pan-
STARRS Y and Pan-STARRS zs are not available among the filters in EXOFASTv2, y
and z’ (Sloan z) were used as respective approximate substitutes.

Telescope | Date (UT) | Filter | Detrending Parameter(s) | Note
---|---|---|---|---
Spitzer | 2019-02-10 | 4.5 $\mu$m | x, y, noise/pixel, FWHM_x, FWHM_y, sky, linear, quadratic | 1
Spitzer | 2019-02-27 | 4.5 $\mu$m | x, y, noise/pixel, FWHM_x, FWHM_y, sky, linear, quadratic, flare | 1
Spitzer | 2019-09-09 | 4.5 $\mu$m | x, y, noise/pixel, FWHM_x, FWHM_y, sky, linear, quadratic, Gaussian | 1
LCO SSO | 2020-04-25 | z’ | Width_T1 | 2
LCO SSO | 2020-04-25 | y | AIRMASS | 2
LCO SAAO | 2020-05-20 | z’ | AIRMASS | 2
LCO SAAO | 2020-05-20 | y | AIRMASS | 2
LCO SAAO | 2020-06-06 | z’ | AIRMASS | 2
LCO SAAO | 2020-06-23 | z’ | AIRMASS | 2
PEST | 2020-07-10 | V | comp_flux, x_coord, y_coord, dist_center, fwhm, airmass, sky | 3
Brierfield | 2020-08-13 | I | Meridian_Flip, X(FITS)_T1, Y(FITS)_T1, tot_C_cnts | 2
LCO SSO | 2020-08-13 | z’ | AIRMASS, Width_T1 | 2
LCO SAAO | 2020-09-07 | z’ | Width_T1, Sky/Pixel_T1 | 2
LCO SSO | 2020-09-16 | z’ | Width_T1 | 2
LCO SSO | 2020-10-03 | z’ | Width_T1 | 2
LCO SAAO | 2020-10-11 | z’ | Width_T1 | 2

11footnotetext: See $\S$2.2 for details on the detrending parameters applied
to Spitzer data.

22footnotetext: Detrending parameters generated from AIJ (Collins et al.,
2017).

33footnotetext: Detrending parameters generated from PEST pipeline
(http://pestobservatory.com/the-pest-pipeline).

Given that AU Mic (catalog ) is a low-mass red dwarf, we configured EXOFASTv2
to use MIST for evolutionary models (Choi et al., 2016; Dotter, 2016) and
ignore the Claret & Bloemen limb darkening tables (Claret & Bloemen, 2011).
Additionally, we incorporate the spectral energy distribution (SED) to place
constraint on MIST evolutionary models; the bands and their corresponding
magnitude priors are presented in Table 3. We also assume the orbit of both AU
Mic (catalog ) b & c to be non-circular. For the EXOFASTv2 modeling, each of
the 16 observations are detrended as indicated in Table 8; $\S$2 describes
additional details on detrending parameters used for each data set. We split
up the EXOFASTv2 modeling into two sequential MCMC runs. For the first run, we
integrate up to 15 000 steps while setting NTHIN = 12; the first run was also
configured to integrate the priors from Tables 6 & 3 and to invoke the
rejectflatmodel option for all light curves with NTEMPS = 8 to aid in faster
convergence. After the first run, EXOFASTv2 generates the new prior.2 file,
which we then implement while repeating the process to achieve better
convergence. For the second run, we integrate up to 20 000 steps while setting
NTHIN = 5, the rejectflatmodel option was turned off, and the MIST SED file
was omitted. After these runs were completed, EXOFASTv2 generated the
posteriors, including the transit models (Figure 2 and Tables 9, 10, 11, & 12)
and midpoint times (see Table 13). Of particular note are our eccentricity
posteriors of 0.079${}^{+0.160}_{-0.058}$ for AU Mic (catalog ) b and
0.114${}^{+0.120}_{-0.074}$ for AU Mic (catalog ) c, which exclude moderate to
high eccentricities. Additional analyses of the transits individually
indicates that this posterior is most constrained by the Spitzer transits
presented herein.

|  |
---|---|---
|  |
|  |
Figure 2: Two-panel plots of comparison between ground-based + Spitzer \+ TESS
transits (cyan) and EXOFASTv2’s best-fit model (black) for AU Mic b & c. The
modelings are in the upper panels, and the residuals are in the lower panels.
Herein epoch refers to the number of transits since the first transit of
either b or c. The Brierfield transit is obtained at epoch 88; LCO SAAO at 78,
80, 82, 91, & 95; LCO SSO at 75, 88, 92, & 94; PEST at 84; Spitzer at 23, 25,
& 48; and TESS b at 0, 2, 84, 85, & 86, relative to the first TESS b transit.
The TESS c transit is obtained at epoch 0, 37, & 38 relative to the first TESS
c transit.

Table 9: EXOFASTv2-generated median values and 68% confidence interval for AU
Mic system.

Posterior | Description | Unit | Quantity
---|---|---|---
AU Mic (catalog ) b |  AU Mic (catalog ) c
M⋆ | Stellar Mass | M☉ | 0.558${}^{+0.068}_{0.090}$
R⋆ | Stellar Radius | R☉ | 0.749${}^{+0.030}_{-0.037}$
L⋆ | Stellar Luminosity | L☉ | 0.094${}^{+0.014}_{-0.015}$
$\rho_{\star}$ | Stellar Density | g/cm3 | 1.89${}^{+0.24}_{-0.21}$
$\log{g}$ | Surface Gravity | … | 4.443${}^{+0.044}_{-0.048}$
Teff | Effective Temperature | K | 3701${}^{+100}_{-98}$
$[{\rm Fe/H}]$ | Metallicity | … | 0.17${}^{+0.27}_{-0.37}$
$[{\rm Fe/H}]_{0}$ | Initial MetallicityaaThe metallicity of the star at birth. | … | 0.12${}^{+0.24}_{-0.36}$
Age | … | Gyr | 0.0216${}^{+0.0034}_{-0.0031}$
EEP | Equal Evolutionary PhasebbCorresponds to static points in a star’s evolutionary history. See $\S$2 in Dotter (2016). | … | 164.9${}^{+5.4}_{-5.8}$
Porb | Orbital Period | days | 8.462993${}^{+0.000048}_{-0.000039}$ | 18.859005${}^{+0.000084}_{-0.000072}$
Mp | Planetary MassccUses measured radius and estimated mass from Chen & Kipping (2017). | MJ | 0.058${}^{+0.021}_{-0.016}$ | 0.0291${}^{+0.0110}_{0.0080}$
Rp | Planetary Radius | RJ | 0.373${}^{+0.017}_{-0.021}$ | 0.250${}^{+0.013}_{-0.015}$
TC | Time of ConjunctionddTime of Conjunction is commonly reported as the “transit time”. | BJD_TDB | 2458330.3916${}^{+0.0031}_{-0.0030}$ | 2458342.2238${}^{+0.0026}_{-0.0023}$
TT | Time of Minimum Projected SeparationeeTime of Minimum Projected Separation is a more correct “transit time”. | BJD_TDB | 2458330.3916${}^{+0.0031}_{-0.0030}$ | 2458342.2238${}^{+0.0026}_{-0.0023}$
T0 | Optimal Conjunction TimeffOptimal Time of Conjunction minimizes the covariance between TC and Period. | BJD_TDB | 2458457.3365${}^{+0.0031}_{-0.0029}$ | 2458455.3777${}^{+0.0027}_{-0.0022}$
a | Semi-Major Axis | au | 0.0669${}^{+0.0026}_{-0.0038}$ | 0.1141${}^{+0.0044}_{-0.0065}$
e | Eccentricity | … | 0.079${}^{+0.160}_{-0.058}$ | 0.114${}^{+0.120}_{-0.074}$
i | Inclination | deg | 89.72${}^{+0.22}_{-0.29}$ | 89.39${}^{+0.47}_{-0.28}$
$\omega$ | Argument of Periastron | deg | -90 $\pm$ 100 | -90 $\pm$ 100
Teq | Equilibrium TemperatureggAssumes no albedo and perfect redistribution. | K | 595${}^{+18}_{-21}$ | 456${}^{+14}_{-16}$
$\tau_{\rm circ}$ | Tidal Circularization Timescale | Gyr | 152${}^{+81}_{-55}$ | 16400${}^{+9000}_{-7500}$
K | RV Semi-AmplitudeccUses measured radius and estimated mass from Chen & Kipping (2017). | m/s | 8.6${}^{+3.2}_{-2.4}$ | 3.33${}^{+1.30}_{-0.92}$
Rp/R⋆ | … | … | 0.05137${}^{+0.00099}_{-0.00120}$ | 0.03429${}^{+0.00100}_{-0.00099}$
a/R⋆ | … | … | 19.28${}^{+0.79}_{-0.73}$ | 32.9 $\pm$ 1.3
$\delta$ | Transit Depth (Rp/R⋆)2 | … | 0.00264${}^{+0.00010}_{-0.00012}$ | 0.001176${}^{+0.000070}_{-0.000067}$
$\delta_{\rm I}$ | Transit Depth in I | … | 0.00366${}^{+0.00130}_{-0.00072}$ | 0.00152${}^{+0.00048}_{-0.00027}$
$\delta_{\rm z^{\prime}}$ | Transit Depth in z’ | … | 0.00409${}^{+0.00053}_{-0.00044}$ | 0.00169${}^{+0.00020}_{-0.00018}$
$\delta_{\rm 4.5\mu m}$ | Transit Depth in $4.5\mu m$ | … | 0.00265${}^{+0.00010}_{-0.00012}$ | 0.001182${}^{+0.000069}_{-0.000068}$
$\delta_{\rm TESS}$ | Transit Depth in TESS | … | 0.00339${}^{+0.00019}_{-0.00022}$ | 0.001462${}^{+0.000077}_{-0.000093}$
$\delta_{\rm V}$ | Transit Depth in V | … | 0.00341${}^{+0.0014}_{-0.00058}$ | 0.00146${}^{+0.00051}_{-0.00024}$
$\delta_{\rm y}$ | Transit Depth in y | … | 0.00357${}^{+0.00160}_{-0.00068}$ | 0.00149${}^{+0.00057}_{-0.00025}$
$\tau$ | Ingress/Egress Transit Duration | days | 0.00721${}^{+0.00022}_{-0.00025}$ | 0.00664${}^{+0.00180}_{-0.00084}$
T14 | Total Transit Duration | days | 0.14585${}^{+0.00031}_{-0.00036}$ | 0.1773${}^{+0.0017}_{-0.0015}$
TFWHM | FWHM Transit Duration | days | 0.13861${}^{+0.00020}_{-0.00023}$ | 0.1703${}^{+0.0011}_{-0.0013}$
b | Transit Impact Parameter | … | 0.098${}^{+0.092}_{-0.065}$ | 0.35${}^{+0.19}_{-0.26}$
bS | Eclipse Impact Parameter | … | 0.093${}^{+0.091}_{-0.063}$ | 0.346${}^{+0.095}_{-0.260}$
$\tau_{\rm S}$ | Ingress/Egress Eclipse Duration | days | 0.00715${}^{+0.00059}_{-0.00070}$ | 0.00647${}^{+0.00061}_{-0.00065}$
TS,14 | Total Eclipse Duration | days | 0.145${}^{+0.011}_{-0.014}$ | 0.173${}^{+0.020}_{-0.027}$
TS,FWHM | FWHM Eclipse Duration | days | 0.137${}^{+0.010}_{-0.014}$ | 0.166${}^{+0.020}_{-0.026}$
$\delta_{\rm S,2.5\mu m}$ | Blackbody Eclipse Depth at 2.5 $\mu$m | ppm | 0.63${}^{+0.18}_{-0.16}$ | 0.0147${}^{+0.0060}_{-0.0051}$
$\delta_{\rm S,5.0\mu m}$ | Blackbody Eclipse Depth at 5.0 $\mu$m | ppm | 24.9${}^{+3.3}_{-3.4}$ | 2.55${}^{+0.48}_{-0.49}$
$\delta_{\rm S,7.5\mu m}$ | Blackbody Eclipse Depth at 7.5 $\mu$m | ppm | 74.3${}^{+6.7}_{-7.4}$ | 12.2${}^{+1.5}_{-1.7}$
$\rho_{\rm p}$ | Planetary DensityccUses measured radius and estimated mass from Chen & Kipping (2017). | g/cm3 | 1.42${}^{+0.60}_{-0.35}$ | 2.28${}^{+0.94}_{-0.62}$
$\log{g_{\rm p}}$ | Surface GravityccUses measured radius and estimated mass from Chen & Kipping (2017). | … | 3.02${}^{+0.14}_{-0.12}$ | 3.07${}^{+0.14}_{-0.11}$
$\Theta$ | Safronov Number | … | 0.0373${}^{+0.0140}_{-0.0097}$ | 0.048${}^{+0.019}_{-0.013}$
$\langle F\rangle$ | Incident Flux | 109 $\frac{{\rm erg}}{{\rm s}~{}{\rm cm}^{2}}$ | 0.0282${}^{+0.0035}_{-0.0039}$ | 0.0096 $\pm$ 0.0013
TP | Time of Periastron | BJD_TDB | 2458323.3${}^{+2.7}_{-3.0}$ | 2458332.8${}^{+6.1}_{-6.2}$
TS | Time of Eclipse | BJD_TDB | 2458334.62${}^{+0.54}_{-0.50}$ | 2458332.8 $\pm$ 1.3
TA | Time of Ascending Node | BJD_TDB | 2458328.27${}^{+0.30}_{-0.48}$ | 2458337.41${}^{+0.72}_{-1.20}$
TD | Time of Descending Node | BJD_TDB | 2458332.51${}^{+0.50}_{-0.28}$ | 2458347.03${}^{+1.20}_{-0.72}$
Vc/Ve | h | … | 1.000${}^{+0.045}_{-0.038}$ | 1.009${}^{+0.140}_{-0.060}$
$\sqrt{\left(1-\frac{{\rm R}_{\rm p}}{{\rm R}_{\star}}\right)^{2}-\mbox{b}^{2}}$ | Transit Chord | … | 1.0466${}^{+0.0042}_{-0.0130}$ | 0.974${}^{+0.050}_{-0.092}$
sign | i | … | 1.40${}^{+0.40}_{-0.33}$ | 1.48${}^{+0.33}_{-0.34}$
e$\cos{\omega}$ | … | … | -0.001${}^{+0.100}_{-0.093}$ | $0.00^{+0.10}_{-0.11}$
e$\sin{\omega}$ | … | … | -0.004${}^{+0.037}_{-0.054}$ | -0.015${}^{+0.060}_{-0.14}$
Mp $\sin{i}$ | Minimum MassccUses measured radius and estimated mass from Chen & Kipping (2017). | MJ | 0.058${}^{+0.021}_{-0.016}$ | 0.0291${}^{+0.0110}_{-0.0080}$
Mp/M⋆ | Mass RatioccUses measured radius and estimated mass from Chen & Kipping (2017). | … | 0.000099${}^{+0.000039}_{-0.000028}$ | 0.000050${}^{+0.000021}_{-0.000014}$
d/R⋆ | Separation at Mid Transit | … | 19.2${}^{+1.6}_{-1.5}$ | 33.1${}^{+4.2}_{-2.7}$
PT | A Priori Non-grazing Transit Probability | … | 0.0495${}^{+0.0043}_{-0.0038}$ | 0.0291${}^{+0.0026}_{-0.0033}$
PT,G | A Priori Transit Probability | … | 0.0548${}^{+0.0047}_{-0.0041}$ | 0.0312${}^{+0.0028}_{-0.0035}$
PS | A Priori Non-grazing Eclipse Probability | … | 0.04958${}^{+0.00210}_{-0.00067}$ | 0.0299${}^{+0.0038}_{-0.0020}$
PS,G | A Priori Eclipse Probability | … | 0.05492${}^{+0.00240}_{-0.00073}$ | 0.0320${}^{+0.0042}_{-0.0021}$
u1,I | Linear Limb-darkening Coefficient in I | … | 0.57${}^{+0.39}_{-0.36}$
u${}_{\rm 1,z^{\prime}}$ | Linear Limb-darkening Coefficient in z’ | … | 0.72${}^{+0.13}_{-0.14}$
u1,4.5μm | Linear Limb-darkening Coefficient in 4.5 $\mu$m | … | 0.0086${}^{+0.0130}_{-0.0065}$
u1,TESS | Linear Limb-darkening Coefficient in TESS | … | 0.453${}^{+0.059}_{-0.064}$
u1,V | Linear Limb-darkening Coefficient in V | … | 0.47${}^{+0.45}_{-0.33}$
u1,y | Linear Limb-darkening Coefficient in y | … | 0.54${}^{+0.45}_{-0.36}$
u2,I | Quadratic Limb-darkening Coefficient in I | … | 0.01${}^{+0.42}_{-0.33}$
u${}_{\rm 2,z^{\prime}}$ | Quadratic Limb-darkening Coefficient in z’ | … | -0.28${}^{+0.17}_{-0.10}$
u2,4.5μm | Quadratic Limb-darkening Coefficient in 4.5 $\mu$m | … | 0.140${}^{+0.037}_{-0.038}$
u2,TESS | Quadratic Limb-darkening Coefficient in TESS | … | 0.207${}^{+0.100}_{-0.097}$
u2,V | Quadratic Limb-darkening Coefficient in V | … | 0.09${}^{+0.45}_{-0.37}$
u2,y | Quadratic Limb-darkening Coefficient in y | … | 0.02${}^{+0.40}_{-0.34}$

hhfootnotetext: The velocity at TC of an assumed circular orbit divided by the
velocity of the modeled eccentric orbit.

iifootnotetext: The sign of the solution to the quadratic mapping from Vc/Ve
to e.

Note. — See Table 3 in Eastman et al. (2019) for a detailed description of all
parameters. Since both Pan-STARRS Y and Pan-STARRS zs are not available among
the filters in EXOFASTv2, y and z’ (Sloan z) were used as respective
approximate substitutes. Additionally, the Claret & Bloemen limb darkening
tables (Claret & Bloemen, 2011) default option has been disabled since AU Mic
(catalog ) is a low-mass red dwarf.

=30mm

Table 10: EXOFASTv2-generated median values and 68% confidence interval for
follow-up observations of AU Mic transits (Part I).

Planet | Telescope | Date (UT) | Filter |  $\sigma^{2}$ (Added Variance) | TTVaaTransit Timing Variation (days) | T$\delta$ VbbTransit Depth Variation | F0ccBaseline Flux
---|---|---|---|---|---|---|---
b | TESS | 2018-07-26 | TESS | 0.0000000019${}^{+0.0000000100}_{-0.0000000088}$ | -0.0024${}^{+0.0028}_{-0.0031}$ | -0.0006${}^{+0.0013}_{-0.0012}$ | 1.000002${}^{+0.000024}_{-0.000023}$
c | TESS | 2018-08-11 | TESS | 0.0000000030 $\pm$ 0.0000000097 | -0.0005${}^{+0.0020}_{-0.0026}$ | -0.00006${}^{+0.00080}_{-0.00074}$ | 0.999994${}^{+0.000026}_{-0.000027}$
b | TESS | 2018-08-12 | TESS | 0.000000005${}^{+0.000000016}_{-0.000000013}$ | -0.0006${}^{+0.0029}_{-0.0030}$ | -0.0006${}^{+0.0014}_{-0.0012}$ | 0.999999 $\pm$ 0.000030
b | Spitzer | 2019-02-10 | 4.5 $\mu$m | 0.0000000600${}^{+0.0000000050}_{-0.0000000047}$ | 0.0039 $\pm$ 0.0030 | -0.00853${}^{+0.00130}_{-0.00100}$ | 1.000075 $\pm$ 0.000022
b | Spitzer | 2019-02-27 | 4.5 $\mu$m | 0.0000000374${}^{+0.0000000040}_{-0.0000000036}$ | 0.0049${}^{+0.0029}_{-0.0031}$ | -0.0055${}^{+0.0015}_{-0.0012}$ | 1.000220${}^{+0.000019}_{-0.000017}$
b | Spitzer | 2019-09-09 | 4.5 $\mu$m | 0.0000000933${}^{+0.0000000042}_{-0.0000000037}$ | 0.0020${}^{+0.0032}_{-0.0035}$ | -0.0039${}^{+0.0014}_{-0.0011}$ | 0.999893${}^{+0.000013}_{-0.000014}$
b | LCO SSO | 2020-04-25 | z’ | 0.00000234${}^{+0.00000029}_{-0.00000028}$ | -0.0028${}^{+0.0040}_{-0.0048}$ | -0.0054${}^{+0.0047}_{-0.0032}$ | 1.00000${}^{+0.00013}_{-0.00014}$
b | LCO SSO | 2020-04-25 | y | 0.00000104${}^{+0.00000100}_{-0.00000067}$ | -0.0037${}^{+0.0048}_{-0.0059}$ | 0.0037${}^{+0.0046}_{-0.0076}$ | 0.99964${}^{+0.00031}_{-0.00038}$
b | LCO SAAO | 2020-05-20 | z’ | 0.00000188 $\pm$ 0.00000018 | -0.0019${}^{+0.0044}_{-0.0049}$ | 0.0005${}^{+0.0025}_{-0.0027}$ | 1.00034 $\pm$ 0.00018
b | LCO SAAO | 2020-05-20 | y | 0.00000140${}^{+0.00000060}_{-0.00000047}$ | 0.0084${}^{+0.0057}_{-0.0066}$ | 0.0000${}^{+0.0047}_{-0.0052}$ | 0.99993${}^{+0.00044}_{-0.00039}$
b | LCO SAAO | 2020-06-06 | z’ | 0.00000677${}^{+0.00000081}_{-0.00000066}$ | -0.0053${}^{+0.0050}_{-0.0056}$ | -0.0018${}^{+0.0053}_{-0.0050}$ | 1.00005${}^{+0.00041}_{-0.00034}$
b | LCO SAAO | 2020-06-23 | z’ | 0.00000185${}^{+0.00000025}_{-0.00000021}$ | 0.0003${}^{+0.0040}_{-0.0049}$ | 0.0079${}^{+0.0016}_{-0.0031}$ | 0.99987${}^{+0.00018}_{-0.00021}$
c | TESS | 2020-07-09 | TESS | 0.000000001 $\pm$ 0.000000015 | -0.0024${}^{+0.0031}_{-0.0047}$ | 0.00023${}^{+0.00100}_{-0.00063}$ | 1.000001${}^{+0.000020}_{-0.000021}$
b | TESS | 2020-07-10 | TESS | 0.000000002${}^{+0.000000012}_{-0.000000011}$ | -0.0018${}^{+0.0038}_{-0.0047}$ | -0.0006${}^{+0.0013}_{-0.0011}$ | 0.999999${}^{+0.000015}_{-0.000018}$
b | PEST | 2020-07-10 | V | 0.00000600${}^{+0.00000074}_{-0.00000076}$ | -0.0038 $\pm$ 0.0061 | -0.0081${}^{+0.0029}_{-0.0014}$ | 0.99817${}^{+0.00018}_{-0.00015}$
b | TESS | 2020-07-19 | TESS | 0.0000000015${}^{+0.0000000100}_{-0.0000000099}$ | -0.0009${}^{+0.0036}_{-0.0050}$ | -0.0006${}^{+0.0013}_{-0.0011}$ | 1.000000 $\pm$ 0.000015
b | TESS | 2020-07-27 | TESS | 0.0000000011 $\pm$ 0.0000000097 | 0.0002${}^{+0.0037}_{-0.0050}$ | -0.0007${}^{+0.0013}_{-0.0011}$ | 0.999999 $\pm$ 0.000015
c | TESS | 2020-07-28 | TESS | 0.000000001${}^{+0.000000013}_{-0.000000012}$ | 0.0001${}^{+0.0032}_{-0.0039}$ | 0.00022${}^{+0.00097}_{-0.00058}$ | 1.000001${}^{+0.000017}_{-0.000019}$
b | Brierfield | 2020-08-13 | I | 0.0000117${}^{+0.0000011}_{-0.0000010}$ | 0.0085${}^{+0.0047}_{-0.0055}$ | 0.0070${}^{+0.0022}_{-0.0036}$ | 0.99935${}^{+0.00025}_{-0.00028}$
b | LCO SSO | 2020-08-13 | z’ | 0.00000406${}^{+0.00000039}_{-0.00000037}$ | -0.0018${}^{+0.0044}_{-0.0062}$ | -0.0013${}^{+0.0034}_{-0.0036}$ | 1.00009${}^{+0.00030}_{-0.00026}$
b | LCO SAAO | 2020-09-07 | z’ | 0.00000688${}^{+0.00000083}_{-0.00000070}$ | 0.0064${}^{+0.0092}_{-0.0068}$ | -0.0057${}^{+0.0042}_{-0.0030}$ | 0.99988${}^{+0.00045}_{-0.00035}$
b | LCO SSO | 2020-09-16 | z’ | 0.00000131${}^{+0.00000013}_{-0.00000011}$ | -0.0020${}^{+0.0041}_{-0.0054}$ | 0.0000${}^{+0.0021}_{-0.0023}$ | 0.99996 $\pm$ 0.00013
b | LCO SSO | 2020-10-03 | z’ | 0.00000125${}^{+0.00000017}_{-0.00000014}$ | 0.0011${}^{+0.0041}_{-0.0056}$ | 0.0020${}^{+0.0022}_{-0.0024}$ | 1.00216${}^{+0.00017}_{-0.00016}$
b | LCO SAAO | 2020-10-11 | z’ | 0.00000651${}^{+0.00000069}_{-0.00000058}$ | 0.0100${}^{+0.0070}_{-0.0074}$ | 0.0039${}^{+0.0030}_{-0.0033}$ | 1.00099${}^{+0.00048}_{-0.00049}$

Note. — See Table 3 in Eastman et al. (2019) for a detailed description of all
parameters. Since both Pan-STARRS Y and Pan-STARRS zs are not available among
the filters in EXOFASTv2, y and z’ (Sloan z) were used as respective
approximate substitutes.

=30mm

Table 11: EXOFASTv2-generated median values and 68% confidence interval for
follow-up observations of AU Mic b (Part II).

Planet | Telescope | Date (UT) | Filter | C0aaAdditive Detrending Coefficient | C1aaAdditive Detrending Coefficient | M0bbMultiplicative Detrending Coefficient | M1bbMultiplicative Detrending Coefficient | M2bbMultiplicative Detrending Coefficient
---|---|---|---|---|---|---|---|---
b | TESS | 2018-07-26 | TESS | … | … | … | … | …
c | TESS | 2018-08-11 | TESS | … | … | … | … | …
b | TESS | 2018-08-12 | TESS | … | … | … | … | …
b | Spitzer | 2019-02-10 | 4.5 $\mu$m | -0.000165${}^{+0.000043}_{-0.000041}$ | … | -0.00110${}^{+0.00034}_{-0.00033}$ | -0.0112 $\pm$ 0.0010 | 0.0108 $\pm$ 0.0010
b | Spitzer | 2019-02-27 | 4.5 $\mu$m | -0.000247${}^{+0.000044}_{-0.000042}$ | 0.000048${}^{+0.000041}_{-0.000032}$ | 0.00070${}^{+0.00033}_{-0.00037}$ | 0.00037${}^{+0.00037}_{-0.00041}$ | -0.00303${}^{+0.00081}_{-0.00076}$
b | Spitzer | 2019-09-09 | 4.5 $\mu$m | -0.000576${}^{+0.000093}_{-0.000110}$ | … | -0.00030 $\pm$ 0.00016 | 0.00194${}^{+0.00020}_{-0.00021}$ | -0.00022${}^{+0.00061}_{-0.00063}$
b | LCO SSO | 2020-04-25 | z’ | … | … | 0.00000${}^{+0.00024}_{-0.00026}$ | … | …
b | LCO SSO | 2020-04-25 | y | … | … | -0.00111${}^{+0.00074}_{-0.00079}$ | … | …
b | LCO SAAO | 2020-05-20 | z’ | … | … | 0.00102 $\pm$ 0.00035 | … | …
b | LCO SAAO | 2020-05-20 | y | … | … | -0.00007${}^{+0.00071}_{-0.00070}$ | … | …
b | LCO SAAO | 2020-06-06 | z’ | … | … | 0.00025${}^{+0.00081}_{-0.00074}$ | … | …
b | LCO SAAO | 2020-06-23 | z’ | … | … | -0.00021${}^{+0.00038}_{-0.00049}$ | … | …
c | TESS | 2020-07-09 | TESS | … | … | … | … | …
b | TESS | 2020-07-10 | TESS | … | … | … | … | …
b | PEST | 2020-07-10 | V | -0.0037${}^{+0.0050}_{-0.0049}$ | … | 0.1038 $\pm$ 0.0040 | -0.1074 $\pm$ 0.0043 | -0.0002${}^{+0.0017}_{-0.0016}$
b | TESS | 2020-07-19 | TESS | … | … | … | … | …
b | TESS | 2020-07-27 | TESS | … | … | … | … | …
c | TESS | 2020-07-28 | TESS | … | … | … | … | …
b | Brierfield | 2020-08-13 | I | … | … | -0.0038 $\pm$ 0.0052 | -0.0026${}^{+0.0045}_{-0.0047}$ | -0.0009${}^{+0.0022}_{-0.0023}$
b | LCO SSO | 2020-08-13 | z’ | … | … | -0.00002${}^{+0.00043}_{-0.00052}$ | 0.00028${}^{+0.00029}_{-0.00033}$ | …
b | LCO SAAO | 2020-09-07 | z’ | -$0.00020^{+0.00063}_{-0.00066}$ | … | -0.00035 $\pm$ 0.00096 | … | …
b | LCO SSO | 2020-09-16 | z’ | … | … | 0.00014${}^{+0.00033}_{-0.00029}$ | … | …
b | LCO SSO | 2020-10-03 | z’ | … | … | -0.00179 $\pm$ 0.00028 | … | …
b | LCO SAAO | 2020-10-11 | z’ | … | … | -0.0014 $\pm$ 0.0016 | … | …

Note. — See Table 3 in Eastman et al. (2019) for a detailed description of all
parameters. Since both Pan-STARRS Y and Pan-STARRS zs are not available among
the filters in EXOFASTv2, y and z’ (Sloan z) were used as respective
approximate substitutes. Also see $\S$2 of this paper for details on
detrending parameters used for each observation; the flare (Spitzer), sky
(Spitzer & PEST), & Sky/Pixel_T1 (LCO SAAO) were set as additive while the
remaining detrending parameters were set as multiplicative.

=30mm

Table 12: EXOFASTv2-generated median values and 68% confidence interval for
follow-up observations of AU Mic b (Part III).

Planet | Telescope | Date (UT) | Filter | M3aaMultiplicative Detrending Coefficient | M4aaMultiplicative Detrending Coefficient | M5aaMultiplicative Detrending Coefficient | M6aaMultiplicative Detrending Coefficient | M7aaMultiplicative Detrending Coefficient
---|---|---|---|---|---|---|---|---
b | TESS | 2018-07-26 | TESS | … | … | … | … | …
c | TESS | 2018-08-11 | TESS | … | … | … | … | …
b | TESS | 2018-08-12 | TESS | … | … | … | … | …
b | Spitzer | 2019-02-10 | 4.5 $\mu$m | -0.00623${}^{+0.00075}_{-0.00074}$ | 0.00496 $\pm$ 0.00048 | -0.000407${}^{+0.000041}_{-0.000042}$ | -0.000270${}^{+0.000057}_{-0.000053}$ | …
b | Spitzer | 2019-02-27 | 4.5 $\mu$m | 0.00326${}^{+0.00073}_{-0.00080}$ | 0.00118${}^{+0.00065}_{-0.00060}$ | 0.000873${}^{+0.000095}_{-0.000110}$ | -0.00088${}^{+0.00015}_{-0.00013}$ | …
b | Spitzer | 2019-09-09 | 4.5 $\mu$m | 0.00049 $\pm$ 0.00056 | -0.00174${}^{+0.00046}_{-0.00045}$ | -0.000003${}^{+0.000044}_{-0.000046}$ | -0.000403${}^{+0.000061}_{-0.000057}$ | -0.000519${}^{+0.000059}_{-0.000058}$
b | LCO SSO | 2020-04-25 | z’ | … | … | … | … | …
b | LCO SSO | 2020-04-25 | y | … | … | … | … | …
b | LCO SAAO | 2020-05-20 | z’ | … | … | … | … | …
b | LCO SAAO | 2020-05-20 | y | … | … | … | … | …
b | LCO SAAO | 2020-06-06 | z’ | … | … | … | … | …
b | LCO SAAO | 2020-06-23 | z’ | … | … | … | … | …
c | TESS | 2020-07-09 | TESS | … | … | … | … | …
b | TESS | 2020-07-10 | TESS | … | … | … | … | …
b | PEST | 2020-07-10 | V | 0.0004${}^{+0.0012}_{-0.0014}$ | -0.0023${}^{+0.0028}_{-0.0026}$ | 0.0013 $\pm$ 0.0011 | -0.00284${}^{+0.00088}_{-0.00087}$ | …
b | TESS | 2020-07-19 | TESS | … | … | … | … | …
b | TESS | 2020-07-27 | TESS | … | … | … | … | …
c | TESS | 2020-07-28 | TESS | … | … | … | … | …
b | Brierfield | 2020-08-13 | I | -0.0009${}^{+0.0015}_{-0.0014}$ | … | … | … | …
b | LCO SSO | 2020-08-13 | z’ | … | … | … | … | …
b | LCO SAAO | 2020-09-07 | z’ | … | … | … | … | …
b | LCO SSO | 2020-09-16 | z’ | … | … | … | … | …
b | LCO SSO | 2020-10-03 | z’ | … | … | … | … | …
b | LCO SAAO | 2020-10-11 | z’ | … | … | … | … | …

Note. — See Table 3 in Eastman et al. (2019) for a detailed description of all
parameters. Since both Pan-STARRS Y and Pan-STARRS zs are not available among
the filters in EXOFASTv2, y and z’ (Sloan z) were used as respective
approximate substitutes. Also see $\S$2 of this paper for details on
detrending parameters used for each observation; the flare (Spitzer), sky
(Spitzer & PEST), & Sky/Pixel_T1 (LCO SAAO) were set as additive while the
remaining detrending parameters were set as multiplicative.

##  4 TTV Analysis

In this section, we present our O–C diagram and TTV dynamical modeling with
Exo-Striker, considering both a two-planet model and an example three-planet
model to account for the observed TTVs.

###  4.1 O–C Diagram

We calculate the expected midpoint times using AU Mic (catalog ) b’s period
and TC from Gilbert et al. (2021); the period and TC from Martioli et al.
(2021) yield similar results. Using the measured and expected midpoint times,
we construct the O–C diagram in Figure 3. We make use of the combined transit
midpoint times (Table 13), now extracted from all light curves and R-M
observations through the processes described in the previous $\S$2 for data
from Gilbert et al. (2021) and the R-M observations, and from our own analysis
in $\S$3. With a $\chi^{2}\gg 1$, it is readily apparent by eye that the
precise Spitzer transit times are significantly deviant from those expected
from a linear ephemeris from the TESS transit times alone. Note, the Spitzer
transit times are in BJD and corrected for the relative light-time travel
delay between Spitzer and the Solar System barycenter. Additionally, the
ground-based R-M transit midpoints are similarly late and consistent with the
Spitzer transit times. The ground-based photometric transits show large
scatter and correspondingly larger timing uncertainties relative to the space-
based transit midpoint times.

Figure 3: O–C Diagram of AU Mic (catalog ) b, using the EXOFASTv2-generated
measured midpoint times and the calculated expected midpoint times for all 23
AU Mic (catalog ) b transit data sets from Table 2. AU Mic (catalog ) b’s
period and TC from Gilbert et al. (2021) were used for the calculation of
expected midpoint times. Table 13: AU Mic (catalog ) planets’ midpoint time
priors for Exo-Striker models.

Planet | Telescope | Transit N | Date (UT) | Band | T0 (BJD)
---|---|---|---|---|---
b | TESS | 1 | 2018-07-26 | TESS | 2458330.38920 $\pm$ 0.00041
TESS | 3 | 2018-08-12 | TESS | 2458347.31699 $\pm$ 0.00060
Spitzer | 24 | 2019-02-10 | 4.5 $\mu$m | 2458525.04439 $\pm$ 0.00017
Spitzer | 26 | 2019-02-27 | 4.5 $\mu$m | 2458541.97136 $\pm$ 0.00016
SPIRou \+ iSHELL | 39 | 2019-06-17 | (a) | 2458651.99020 $\pm$ 0.00180
ESPRESSO | 45 | 2019-08-07 | 378.2-788.7 nm | 2458702.77397 $\pm$ 0.00178
Spitzer | 49 | 2019-09-09 | 4.5 $\mu$m | 2458736.61730 $\pm$ 0.00015
LCO SSO | 76 | 2020-04-25 | Pan-STARRS Y | 2458965.11250 $\pm$ 0.00300
LCO SSO | 76 | 2020-04-25 | Pan-STARRS zs | 2458965.11300 $\pm$ 0.00140
LCO SAAO | 79 | 2020-05-20 | Pan-STARRS zs | 2458990.50270 $\pm$ 0.00240
LCO SAAO | 79 | 2020-05-20 | Pan-STARRS Y | 2458990.51350 $\pm$ 0.00430
LCO SAAO | 81 | 2020-06-06 | Pan-STARRS zs | 2459007.42580 $\pm$ 0.00340
LCO SAAO | 83 | 2020-06-23 | Pan-STARRS zs | 2459024.35740 $\pm$ 0.00140
PEST | 85 | 2020-07-10 | V | 2459041.27900 $\pm$ 0.00570
TESS | 85 | 2020-07-10 | TESS | 2459041.28120 $\pm$ 0.00030
TESS | 86 | 2020-07-19 | TESS | 2459049.74491 $\pm$ 0.00028
TESS | 87 | 2020-07-27 | TESS | 2459058.20901 $\pm$ 0.00026
LCO SSO | 89 | 2020-08-13 | Pan-STARRS zs | 2459075.13340 $\pm$ 0.00250
Brierfield | 89 | 2020-08-13 | I | 2459075.14340 $\pm$ 0.00240
LCO SAAO | 92 | 2020-09-07 | Pan-STARRS zs | 2459100.52910 $\pm$ 0.00980
LCO SSO | 93 | 2020-09-16 | Pan-STARRS zs | 2459108.98502 $\pm$ 0.00092
LCO SSO | 95 | 2020-10-03 | Pan-STARRS zs | 2459125.91400 $\pm$ 0.00110
LCO SAAO | 96 | 2020-10-11 | Pan-STARRS zs | 2459134.38820 $\pm$ 0.00990
c | TESS | 1 | 2018-08-11 | TESS | 2458342.22330 $\pm$ 0.00110
TESS | 38 | 2020-07-09 | TESS | 2459040.00432 $\pm$ 0.00082
TESS | 39 | 2020-07-28 | TESS | 2459058.86603 $\pm$ 0.00072

Note. — (a) 955-2 515 nm (SPIRou) and 2.18-2.47 nm (iSHELL)

###  4.2 TTV Exo-Striker Dynamical Modeling

Motivated by the apparent TTV variability deviating from a linear ephemeris in
$\S$4.1, we utilize the Exo-Striker package (Trifonov, 2019) to model the
variation in transit timings of AU Mic (catalog ) b and c. Like EXOFASTv2,
Exo-Striker utilizes a Markov chain Monte Carlo (MCMC) to assess the
statistical significance of the measured TTVs and the confidence in the
dynamical model posteriors that can be inferred from the TTVs. Additionally,
we incorporate the stellar priors from Table 14 and the planet priors from
Table 15.

Table 14: AU Mic (catalog )’s stellar priors for Exo-Striker best-fit and MCMC
modeling.

Prior | Unit | AU Mic
---|---|---
Mass | M☉ |  $\mathcal{N}$(0.50, 0.03)
Radius | R☉ |  $\mathcal{N}$(0.75, 0.03)
Luminosity | L☉ |  $\mathcal{N}$(0.0900, 0.0001)
Teff | K |  $\mathcal{N}$(3700, 100)
v sin i | km/s |  $\mathcal{N}$(8.7, 0.2)

Note. — These priors are taken from Plavchan et al. (2020).

Table 15: AU Mic (catalog )’s planetary priors for Exo-Striker best-fit and
MCMC modeling.

Prior | Unit |  AU Mic (catalog ) b |  AU Mic (catalog ) c |  AU Mic (catalog ) daafootnotemark:
---|---|---|---|---
K | m/s |  $\mathcal{N}$(8.5, 2.3) |  $\mathcal{U}$(0.8, 9.5) |  $\mathcal{U}$(0.0, 10000.0)
Porb | day |  $\mathcal{N}$(8.463, 0.001) |  $\mathcal{N}$(18.859, 0.001) |  $\mathcal{N}$(12.742, 0.020)
e | … |  $\mathcal{U}$(0.00000, 0.58038) |  $\mathcal{U}$(0.00000, 0.37308) |  $\mathcal{U}$(0.000, 0.999)
$\omega$ | deg |  $\mathcal{U}$(0.0, 360.0) |  $\mathcal{U}$(0.0, 360.0) |  $\mathcal{U}$(0.0, 360.0)
M0 | deg |  $\mathcal{U}$(0.0, 360.0) |  $\mathcal{U}$(0.0, 360.0) |  $\mathcal{U}$(0.0, 360.0)
i | deg |  $\mathcal{N}$(89.5, 0.3) |  $\mathcal{N}$(89.0, 0.5) |  $\mathcal{U}$(0.0, 180.0)
$\Omega$ | deg |  $\mathcal{U}$(0.0, 360.0) |  $\mathcal{U}$(0.0, 360.0) |  $\mathcal{U}$(0.0, 360.0)

Note. — M0 $\equiv$ mean anomaly, and $\Omega$ $\equiv$ longitude of ascending
node.

(a) For 3-planet model only.

We configured Exo-Striker to use the Simplex and Dynamical algorithms during
both the model fitting and MCMC runs; the justification for using Dynamical
instead of Keplerian is that the Keplerian algorithm does not work as well in
a compact system with orbital resonances (Fabrycky, 2010) like AU Mic (catalog
). The dynamical model time steps has been set to the lowest possible 0.01 day
given the short orbital periods of both planets. To find the best-fit TTV
model, we employ Exo-Striker’s built-in scipy minimizer algorithms; we use
truncated Newton algorithm
(TNC)191919https://docs.scipy.org/doc/scipy/reference/optimize.minimize-
tnc.html as a primary minimizer and Nelder-Mead
algorithm202020https://docs.scipy.org/doc/scipy/reference/optimize.minimize-
neldermead.html as a secondary minimizer, with both set at one consecutive
integration and 5 000 integration steps; the rest of the configurations for
those minimizers were left to default settings. After we find a best-fit
model, we compute an MCMC with 50 000 burn-ins, 200 000 integration steps, and
4 walkers; the other settings for MCMC were left to defaults, including
adopting 68.300% confidence intervals for estimating one-$\sigma$ posterior
uncertainties. These best-fit and MCMC configurations apply to both the
2-planet and 3-planet dynamical models presented herein.

We first explore the 2-planet scenario (AU Mic (catalog ) b and c). Then, we
explore a representative 3-planet scenario. The following sub-subsections
detail the process in exploring these cases.

####  4.2.1 2-Planet Dynamical Modeling

We explored a best-fit scenario for a 2-planet model with Exo-Striker. The
eccentricity posteriors from the EXOFASTv2 analysis (Table 9) provided us
$\sim$4.056$\sigma$ upper limits on both planets’ eccentricity. Analysis of
the transit light curves themselves exclude such high eccentricities as in
Gilbert et al. (2021) and Plavchan et al. (2020), but we are only modeling the
transit midpoint times herein. We also explored a best-fit scenario for a
“mass-less” no-TTV 2-planet model with Exo-Striker, as a control on testing
for the presence and statistical significance of TTVs.

####  4.2.2 3-Planet Dynamical Modeling

Cale et al. (2021) explores the analysis of RVs of AU Mic (catalog ), and
searches for additional candidate non-transiting planet signals. In
particular, Cale et al. (2021) identifies a candidate RV signal in-between the
orbits of AU Mic (catalog ) b and c with an orbital period of 12.742 days,
which to date is unconfirmed, which we call a hypothetical “d” planet. This
“middle-d” non-transiting planet scenario, if confirmed, would establish the
AU Mic (catalog ) planetary system in a 4:6:9 orbital resonant chain. This
would not be the first known non-transiting planet to exist between two
transiting planets; Christiansen et al. (2017), Buchhave et al. (2016), Sun et
al. (2019), and Osborn et al. (2021) identified similar exoplanet
configurations for the HD 3167, Kepler-20, Kepler-411, and TOI-431 planetary
systems respectively. Similarly to the 2-planet model, we imposed 4$\sigma$
upper limits on planets b’s and c’s eccentricity for this paper.

###  4.3 TTV Results

Here we present the results of both the no-TTV, 2-planet and 3-planet TTV
modeling using the Exo-Striker package. We also calculate the mass of AU Mic
(catalog ) c using the generated posteriors from Exo-Striker.

####  4.3.1 Results on Exo-Striker 2-Planet Modeling

|
---|---
Figure 4: 2-planet O–C diagram of AU Mic (catalog ) b (left) and AU Mic (catalog ) c (right), with comparison between TTVs (green) and Exo-Striker-generated best-fit models (black). Table 16: Exo-Striker-generated parameters for AU Mic (catalog ) 2-planet best-fit and MCMC modeling. Parameter | Unit | Best-fit | MCMC
---|---|---|---
AU Mic (catalog ) b |  AU Mic (catalog ) c |  AU Mic (catalog ) b |  AU Mic (catalog ) c
K | m/s | 4.03906 | 9.49903 | 4.10830 $\pm$ 1.47951 | 9.31429 $\pm$ 0.27661
Porb | day | 8.46255 | 18.86109 | 8.46257 $\pm$ 0.00003 | 18.86112 $\pm$ 0.00076
e | … | 0.00000 | 0.08329 | 0.00089 $\pm$ 0.00065 | 0.08220 $\pm$ 0.00473
$\omega$ | deg | 89.79621 | 216.71379 | 61.89388 $\pm$ 36.98216 | 215.56678 $\pm$ 1.90214
M0 | deg | 0.00000 | 0.00032 | 27.96061 $\pm$ 20.80132 | 1.12536 $\pm$ 0.80464
i | deg | 89.50262 | 88.99146 | 89.51412 $\pm$ 0.31257 | 88.95986 $\pm$ 0.51613
$\Omega$ | deg | 0.00000 | 0.00006 | 185.70005 $\pm$ 116.59810 | 184.90569 $\pm$ 116.58777
$-\ln\mathcal{L}$ | … | -61.93923 | -74.96317
$\chi^{2}$ | … | 428 | 455
$\chi^{2}_{\rm red}$ | … | 36 | 38
$\Delta AICc$ | … | -57.69664 | -83.74452

The best-fit O–C diagram (Figure 4), posteriors (Table 16), and MCMC corner
plot (Figure 12) showcase the best possible 2-planet model. The Exo-Striker’s
Angular Momentum Deficit (AMD, Laskar, 1997, 2000; Laskar & Petit, 2017)
criteria indicated that the 2-planet model is stable. However, it is clear
from Figure 4 that the model does not converge very well with the data points,
especially with TESS and Spitzer. This implies that we potentially need a
third planet to account for the observed TTV variability. This interpretation
is one of several hypotheses that are being explored and which are discussed
in greater details in $\S$6. For our control no-TTV scenario of “mass-less”
planets b and c, we present the best fit O-C diagram (Figure 5), posteriors
(Table 17), and corner plot (Figure 13).

|
---|---
Figure 5: Mass-less planets O–C diagram of AU Mic (catalog ) b (left) and AU Mic (catalog ) c (right), with comparison between TTVs (green) and Exo-Striker-generated best-fit models (black). Table 17: Exo-Striker-generated parameters for AU Mic (catalog ) mass-less planets best-fit and MCMC modeling. Parameter | Unit | Best-fit | MCMC
---|---|---|---
AU Mic (catalog ) b |  AU Mic (catalog ) c |  AU Mic (catalog ) b |  AU Mic (catalog ) c
K | m/s | 0.00000 | 0.00000 | 0.00000 $\pm$ 0.00000 | 0.00000 $\pm$ 0.00000
Porb | day | 8.46293 | 19.04753 | 8.46293 $\pm$ 0.00000 | 18.98613 $\pm$ 0.00002
e | … | 0.00004 | 0.00001 | 0.47442 $\pm$ 0.09381 | 0.37308 $\pm$ 0.00000
$\omega$ | deg | 223.78788 | 0.00007 | 230.38681 $\pm$ 18.24298 | 0.00000 $\pm$ 0.00000
M0 | deg | 225.89190 | 0.00010 | 193.13783 $\pm$ 7.76425 | 0.00000 $\pm$ 0.00000
i | deg | 90.73560 | 88.59738 | 89.47176 $\pm$ 0.29016 | 69.69803 $\pm$ 0.00001
$\Omega$ | deg | 0.00002 | 0.00000 | 189.38073 $\pm$ 126.18823 | 180.59530 $\pm$ 125.12807
$-\ln\mathcal{L}$ | … | -20683716.47978 | -6083619013.93697
$\chi^{2}$ | … | 41367738 | 12167264028
$\chi^{2}_{\rm red}$ | … | 2954838 | 869090288
$\Delta AICc$ | … | -41367384.95956 | -12167237979.87394

####  4.3.2 Results on Exo-Striker 3-Planet Modeling

Since the 2-planet circular orbit model does not adequately reproduce the
observed TTVs as evidenced in Figure 4, we explored a representative and non-
exhaustive hypothetical 3-planet scenario. We find a best-fit and MCMC model
for this representative 3-planet scenario that adequately accounts for the
TESS and Spitzer TTVs and is also consistent with the RV candidate signal in
Cale et al. (2021), presented in Table 18 and Figures 6 & 14. Moreover, the
3-planet case’s log likelihood, $\chi^{2}$, and reduced-$\chi^{2}$ are better
than those of the 2-planet case, the latter of which are better than the no-
TTV “mass-less” scenario. The delta log-likelihoods and corresponding $\Delta
AICc$ indicate that the 3-planet scenario is strongly favored and the two
planet scenario and the no-TTV scenario are statistically excluded. Exo-
Striker indicates that this three-planet solution fails the AMD criterion.
However, given the 4:6:9 orbital period commensurability for b, d, and c,
respectively, we investigate the dynamical stability of this system
configuration with two N-body codes, the latter of which is used for a
consistency check. We generate simulations with rebound (Rein & Liu, 2012;
Rein & Spiegel, 2015) and mercury6 (Chambers, 1999) to test the stability of
this representative 3-planet system; both indicate that this 3-planet
configuration is stable (see $\S$6 for more detailed discussions).

|
---|---
Figure 6: 3-planet O–C diagram of AU Mic (catalog ) b (left) and AU Mic (catalog ) c (right), with comparison between TTVs (green) and Exo-Striker-generated best-fit models (black). Table 18: Exo-Striker-generated parameters for AU Mic (catalog ) 3-planet best-fit and MCMC modeling. Parameter | Unit | Best-fit | MCMC
---|---|---|---
AU Mic (catalog ) b |  AU Mic (catalog ) c |  AU Mic (catalog ) d |  AU Mic (catalog ) b |  AU Mic (catalog ) c |  AU Mic (catalog ) d
K | m/s | 17.40198 | 7.65369 | 5.07363 | 1.91815 $\pm$ 1.27252 | 9.21910 $\pm$ 0.50712 | 1.06232 $\pm$ 0.30432
Porb | day | 8.46340 | 18.85872 | 13.48517 | 8.46329 $\pm$ 0.00027 | 18.86224 $\pm$ 0.00126 | 13.46591 $\pm$ 0.01022
e | … | 0.02348 | 0.00000 | 0.00000 | 0.07436 $\pm$ 0.00994 | 0.02508 $\pm$ 0.01531 | 0.00998 $\pm$ 0.00752
$\omega$ | deg | 89.96574 | 223.91754 | 70.01180 | 85.48567 $\pm$ 6.11803 | 210.12425 $\pm$ 16.99204 | 43.40723 $\pm$ 31.78809
M0 | deg | 0.00000 | 0.00000 | 0.00000 | 5.18482 $\pm$ 3.84571 | 11.58305 $\pm$ 8.87570 | 42.36104 $\pm$ 29.98678
i | deg | 89.47231 | 89.10159 | 115.68159 | 89.43632 $\pm$ 0.30223 | 88.99503 $\pm$ 0.54176 | 103.58128 $\pm$ 4.14661
$\Omega$ | deg | 0.00000 | 0.17063 | 0.00000 | 123.30784 $\pm$ 87.12008 | 162.09943 $\pm$ 91.36624 | 139.42421 $\pm$ 90.81975
$-\ln\mathcal{L}$ | … | 73.34798 | -848.73127
$\chi^{2}$ | … | 158 | 1991
$\chi^{2}_{\rm red}$ | … | 32 | 398
$\Delta AICc$ | … | 419.69596 | -1424.46254
Table 19: Comparisons of best-fit model fitting parameters among the mass-less planets, 2-planet, and 3-planet cases. Best-fit Model | $\chi^{2}$ | $\chi^{2}_{\rm red}$ | $-\ln\mathcal{L}$ | $\Delta AICc$
---|---|---|---|---
No TTVs | 41367738 | 2954838 | -20683716.480 | -41367384.960
2-planet | 428 | 36 | -61.939 | -57.697
3-planet | 158 | 32 | 73.348 | 419.696

####  4.3.3 The Mass of AU Mic c

We use equation (1) from Cumming et al. (1999) and the best-fit parameters &
MCMC uncertainties from Tables 16 and 18 to calculate the mass of AU Mic
(catalog ) c. To simplify the equation, we made an assumption that M⋆ $>>$ Mp.
For the best-fit 2-planet scenario, our calculation yields the mass of c to be
0.0781 $\pm$ 0.0039 MJ (or 24.8 $\pm$ 1.2 M⊕) at 20.7-$\sigma$ significance,
which makes it roughly a Neptune-sized planet. In the case of our
representative 3-planet scenario, our calculation yields the mass of c to be
0.0631 $\pm$ 0.0049 MJ (or 20.1 $\pm$ 1.6 M⊕) at 12.6-$\sigma$ significance,
implying that it again has roughly the mass of Neptune.

##  5 Photodynamical Analysis

Table 20: Parameters for the photodynamical models, with model a the
simultaneous fit of stellar activity variation & flares and model b the pre-
cleaned data.

Parameter | Unit | Model a | Model b
---|---|---|---
AU Mic (catalog ) b |  AU Mic (catalog ) c |  AU Mic (catalog ) b |  AU Mic (catalog ) c
Mp | M⊕ | 16${}^{+12}_{-9}$ | 10.8${}^{+2.3}_{-2.2}$ | 13${}^{+8}_{-2}$ | 13${}^{+5}_{-2}$
Porb | day | 8.4626${}^{+0.0001}_{-0.0001}$ | 18.8624${}^{+0.0015}_{-0.0011}$ | 8.4626${}^{+0.0001}_{-0.0002}$ | 18.860${}^{+0.002}_{-0.002}$
e | … | 0.022${}^{+0.005}_{-0.006}$ | 0.097${}^{+0.010}_{-0.013}$ | 0.012${}^{+0.003}_{-0.003}$ | 0.069${}^{+0.025}_{-0.014}$
$\omega$ | deg | 84${}^{+7}_{-9}$ | 99.1${}^{+3.7}_{-3.5}$ | 103${}^{+6}_{-7}$ | 112${}^{+5}_{-9}$
tperi | day | 4.87${}^{+0.16}_{-0.21}$ | 4.89${}^{+0.15}_{-0.16}$ | 5.3${}^{+0.2}_{-0.2}$ | 5.5${}^{+0.2}_{-0.4}$
Rp/R⋆ | … | 0.0529${}^{+0.0002}_{-0.0002}$ | 0.032${}^{+0.002}_{-0.001}$ | 0.0511${}^{+0.0002}_{-0.0001}$ | 0.034${}^{+0.006}_{-0.001}$
i | deg | 89.20${}^{+0.06}_{-0.06}$ | 90.8${}^{+0.1}_{-0.1}$ | 89.11${}^{+0.03}_{-0.03}$ | 89.4${}^{+0.08}_{-0.07}$

Note. — Orbital elements are given for BJD 2458300.

In a second independent analysis to validate our methods, we used a
photodynamical code to simultaneously fit the transit light curves and TTVs.
The core of this code is based on rebound (Rein & Liu, 2012; Rein & Spiegel,
2015) with which the N-body problem is integrated using the high accuracy non-
symplectic integrator with adaptive time-stepping IAS15. At the times of the
measured mid-transit times, the current orbital elements are used to calculate
the corresponding transit light curve using the python implementation of the
model from Mandel & Agol (2002) from Ian
Crossfield212121https://www.lpl.arizona.edu/~ianc/python/_modules/transit.html.
Using the Markov chain Monte Carlo sample emcee (Foreman-Mackey et al., 2013),
we sample the parameter space maximizing the $\log$ likelihood. In addition to
the transit light curves, the code also fits the radial velocity (RV) data of
AU Mic (catalog ) taken with SPIRou (Cale et al., 2021; Klein et al., 2021),
which at near-infrared wavelengths is least impacted by the stellar activity
of this young system. The RV model is calculated from the same rebound N-body
integration of the planetary system. As free parameters we have the planetary
masses, orbital periods, eccentricities, longitudes of periastron,
inclination, and the Rp/R⋆ ratio. The stellar mass and radius, hence the
stellar density, is fixed. The limb darkening for the various bands are taken
from (Claret et al., 2012; Claret, 2017) and are kept at the literature
values. Since the Spitzer transits show depth variations, we account for this
effect by a third light contribution as free parameter (to be explored further
in future work). All data sets have their individual offsets, while the SPIRou
data are modeled by adding a “jitter” white noise term.

The stellar activity is modeled using celerite (Foreman-Mackey et al., 2017),
which has the advantage to offer a fast and scaleable implementation of
Gaussian Process (GP) regression, especially important for large data sets.
The celerite package provides several built-in covariance kernels, one is
representing damped oscillations driven by white noise called SHO. It has
three parameters: the undamped oscillator frequency or period $\omega_{0}$ =
2$\pi$/P, the quality factor $Q$ of the oscillator (which is reversely
proportional to the damping time scale $\tau$), and the variance $S_{0}$. For
more details, see equations (19)-(24) in Foreman-Mackey et al. (2017). Star
spots typically manifest in variations at the rotation period as well as on
its first harmonic. In order to represent this in a rotational kernel, we
follow the idea presented in celerite2 (Foreman-Mackey et al., 2017; Foreman-
Mackey, 2018)222222https://celerite2.readthedocs.io/en/latest/api/python and
use the sum of two SHO kernels, where only one of the frequencies is a free
hyper-parameter and the other is fixed to its first harmonic. Both oscillators
are forced by a lower boundary of the quality factor to be in the underdamped
regime (i.e., $Q>1/2$). For the RV data and the TESS light curves, the same
rotation period and quality factors are enforced; however, the variance can
differ.

We performed the photodynamical modeling in two versions. In the first (model
a, using 45 000 steps after 10 000 burn-in steps and 400 walkers), we
simultaneously fit the flares in the TESS light curves using the flare model
aflare (Davenport et al., 2014). In the second version (model b, using 10 000
steps after 10 000 burn-in steps and 100 walkers, will be increased), we use
the cleaned TESS light curve from Gilbert et al. (2021), where also the
variation due to stellar activity has been corrected, i.e. no GP modeling is
required here for the transit light curves. The parameters of the
photodynamical fits are listed in Table 20 and are displayed in Figures 15 &
16.

###  5.1 Photodynamical Modeling Results

Figure 7: The first TESS transit of AU Mic (catalog ) b in Sector 1. From top
to bottom: TESS light curve with stellar activity variation and flares
overlayed with the best-fit model a from the photodynamical analysis, best-fit
model a with contribution of stellar activity variation and flares removed,
best-fit model b.

Figure 8: The first TESS transit of AU Mic (catalog ) c in Sector 27. From top
to bottom: TESS light curve with stellar activity variation and flares
overlayed with the best-fit model a from the photodynamical analysis, best-fit
model a with contribution of stellar activity variation and flares removed,
best-fit model b, which shows by eye a clear timing offset when the activity
is modeled separately from the photodynamical modeling.

Figure 9: The first Spitzer transit of AU Mic (catalog ) b, with the best-fit
model a (top) and best-fit model b (bottom) from the photodynamical analysis.

In this subsection, we present a summary of the key results of our
photodynamical modeling. Model a has the advantage to provide a simultaneous
fit of the transits, stellar activity modulation, and flares. Since the
timescales of ingress and egress are similar to those of the flares, a
simultaneous fit should provide a more robust result. The flexibility of the
GP kernels, however, may compensate for small deviations in time between
observed and modelled transits. Comparing the calculated transit midpoints of
a best fit to the a priori measured mid-transit times indeed show deviations
on the order of 10 minutes. Either the transit midpoint measurements are
affected by systematic errors due to the stellar activity and flaring, or the
GP model is too flexible and hides these discrepancies. In model b, this is
prevented due to a pre-cleaning of the light curve, however, with the caveat
that the measured mid-transit times may be affected by activity modulation and
flaring. A comparison between the two approaches shows that we can obtain a
2-planet model which fits all transits (TESS, Spitzer, and ground-based) if
the stellar variability and flares are fitted simultaneously, while all models
compared to the pre-cleaned data show noticeable discrepancies. As an example,
we show the first transit of AU Mic (catalog ) b in Sector 1 of TESS in Figure
7, on the top two panels from the simultaneous stellar variability and flare
fit and on the bottom panel from the pre-cleaned data. We also show the best-
fit for one of the Spitzer light curves in Figures 8 & 9.

While the best fit of the photodynamical model including stellar variability
and flares is promising, we note the above mentioned compensation of
deviations between models and observation due to the flexibility of the GP
model. Further transits will be required to check this model. The timing
discrepancies in the second version of the photodynamical fit corroborate the
conclusion from our TTV fitting (i.e., a 2-planet model is not able to
reproduce the small but detectable TTVs).

##  6 Discussion

In this section we present the key results of our analyses in $\S$6.1,
followed by the impact of the rotational modulation of stellar spots and
plages in $\S$6.2 and flares in $\S$6.3. In $\S$6.4, we assess the likelihood
of an additional middle non-transiting planet candidate to explain in the
observed TTVs, and in $\S$6.5 we compare our analyses’ constraints on the mass
of AU Mic (catalog ) c. Lastly, in $\S$6.6, we discuss the limitation of our
joint modeling.

###  6.1 Key Results of Analyses

In the preceding sections, we modeled the transit observations of AU Mic
(catalog ) b and c to ascertain whether or not transit timing variations were
present, and whether or not they can be accounted for by our existing
knowledge of this planetary system and be used to constrain the dynamical
orbits and masses of the two transiting planets. Given the inherent stellar
activity in this young system, we explored two different and independent
analyses: in $\S$4 a dynamical model of the transit midpoint times, and in
$\S$5 a full photodynamical model. We found in both scenarios indications of
transit timing variations deviating from a linear ephemeris on the order of
ten minutes, particularly when comparing the TESS and Spitzer transit times,
and generally consistent results, albeit with larger timing uncertainties,
from the ground-based transits and R-M observations (Figure 3).

We next explored whether the two transiting planets could account for the
observed TTVs. We first attempted a 2-planet dynamical model for the TTVs with
only the high precision TESS transits, which we were able to find a fit akin
to that derived in Martioli et al. (2021). These transit times were derived
with the joint stellar activity model in Gilbert et al. (2021) who find
evidence for TTVs on the order of a few minutes. However, when including the
Spitzer transits, which are less susceptible to stellar activity at 4.5
$\mu$m, we find that the Spitzer transits are incompatible with this model
within the measured timing uncertainties, showing a significant departure from
a circular two-planet model on the order of $\sim$10 minutes. We next
considered whether orbital eccentricity could account for the observed TTVs,
but derive high eccentricities that are incompatible with the light curves
themselves which exclude high eccentricity scenarios (Plavchan et al., 2020;
Gilbert et al., 2021).

With the full photodynamical analysis in $\S$5, we were able to jointly model
the activity and transits to reproduce the observed transit midpoint times.
However, this left open the possibility that the flexibility of the GP for the
rotational modulation of the starspots could also absorb any of the TTVs
present. In fact, with a second photodynamical model using Gilbert et al.
(2021) activity model detrended transit light curves (e.g. the stellar
activity in corrected serially rather than jointly with the photodynamical
analysis), we found excess TTVs not explained by the 2-planet model alone.

In order to explain the observed excess of TTVs than accounted for by our two-
planet modeling, we next investigate in turn the following possibilities.
First, that our modeling of the rotational modulation of star spots is
underestimating the impact on the TTVs. Second, we consider the impact the
flares have on our TTVs. Third, we consider the possibility of a third planet
in the AU Mic (catalog ) system not yet detected to transit nor through RVs.
Finally, we consider the limitations of our analysis and statistical methods.

###  6.2 Stellar Spots

Stellar spots can significantly impact the recovered stellar and planetary
model parameters (Pont et al., 2007; Czesla et al., 2009; Berta et al., 2011;
Désert et al., 2011; Ballerini et al., 2012) and can produce TTV-like signals
that potentially lead to false-positive detection of non-transiting planets
(Alonso et al., 2009; Sanchis-Ojeda et al., 2011; Sanchis-Ojeda & Winn, 2011;
Oshagh et al., 2012). For Rp/R⋆ $\approx$ 0.05 (as is the case for AU Mic
(catalog ) b), the expected maximum amplitude of the spot-induced TTV in
seconds at visible wavelengths is $AMP$ = 139 $\times$ $f$, where $f$ is the
stellar spot filling factor in percent (Oshagh et al., 2013). If $f=$ 0.25%,
1%, or 3%, the maximum amplitude would be 34.75, 139, or 417 seconds (0.6,
2.3, or 7.0 minutes), respectively. Given the long-lived spot lifetimes for AU
Mic (catalog ), it is not unreasonable to assume that the spot-filling
fraction is $>$3% and thus the impact of the rotational modulation of
starspots could be on the order of $>$7 minutes. Oshagh et al. (2013) added
that if the transit duration and depth priors were fixed (and not modeled),
the TTV amplitude would be smaller than if they were floating.

First, we consider the wavelength dependence of the TTVs. The Spitzer transit
times have greater photometric precision due its larger aperture, a much
higher cadence, and consequently more precise transit midpoint times (to
within 13 seconds). Further, the Spitzer observations are less impacted by the
rotational modulation of stellar activity at 4.5 $\mu$m due to the decreased
flux contrast of any spots; the rotational modulation of stellar activity
should be significantly decreased in amplitude. This decreased amplitude of
the impacts of stellar activity with increasing wavelength is also observed
for AU Mic (catalog ) in photometry (Hebb et al., 2007) and in radial
velocities (Cale et al., 2021). However, the Spitzer transits show the most
significant deviations from a linear ephemeris than those derived by the TESS
data alone, running counter-intuitive to the expected impact the rotational
modulation of starspots would have on the derived transit times.

Second, Szabó et al. (2021) measured AU Mic (catalog )’s rotation period with
TESS data and AU Mic (catalog ) b’s orbital period with both TESS and CHEOPS
light curves and found that the ratio between these two periods implies the
7:4 spin-orbit commensurability. Due to this resonance, for every fourth
transit of AU Mic (catalog ) b, the transit takes place over the same range of
stellar longitudes of AU Mic (catalog ). In other words, the same stellar
“side” is transited every fourth transit. Thus the impact of the stellar
rotational modulation of spots should produce the same systematic impact on
the TTVs every fourth transit. For example, if the 1st transit of AU Mic
(catalog ) b is 4 minutes “late” from the apparent effects of starspots, then
this will also be true for the 5th, 9th, 13th, 17th, …, Nth transits. Given AU
Mic (catalog )’s long spot lifetime, relatively unchanged between TESS sectors
1 and 27 two years later, this pattern will have persisted over the duration
of our observations.

|  |
---|---|---

Figure 10: O$-$C Diagram of AU Mic (catalog ) b from initial model (left),
Exo-Striker 2-planet model (center) and Exo-Striker 3-planet model (right)

, color-coded as a function of transit number modulo 4, using the
EXOFASTv2-generated measured midpoint times and the calculated midpoint times
for all 23 AU Mic (catalog ) b transit data sets from Table 2. Here, we
explore the impact the rotational modulation of starspots, given the 7:4 spin-
orbit commensurability, may have on the observed TTVs. The weighted averages
of the distribution of TTVs for each of the four stellar longitude crossing
phases are in Table 21; they show correlations between phases 1 & 3 and phases
2 & 4\. However, after subtracting the 2 and 3 planet models from the TTV
data, the weighted averages show no clear correlation of more than $\sim$2
minutes, indicating that stellar spots of AU Mic (catalog ) are not a
significant factor in accounting for the timing of transits for AU Mic
(catalog ) b. The epochs are relative to the first TESS transit of AU Mic
(catalog ) b.

Table 21: Weighted averages of the respective initial, 2-planet-model-subtracted, & 3-planet-model-subtracted TTVs as a function of transit number modulo 4. Phase | TTVi (min) | TTV2 (min) | TTV3 (min)
---|---|---|---
1 | 0.52 $\pm$ 1.17 | -0.90 $\pm$ 0.59 | -0.40 $\pm$ 0.59
2 | 4.61 $\pm$ 1.00 | 1.11 $\pm$ 0.52 | 0.06 $\pm$ 0.52
3 | 0.99 $\pm$ 1.28 | 0.61 $\pm$ 0.79 | 1.88 $\pm$ 0.79
4 | 4.89 $\pm$ 1.13 | -0.36 $\pm$ 0.54 | -0.40 $\pm$ 0.54

To explore the impact of spot modulations on TTVs, we calculated the weighted
average of the initial, 2-planet model, and 3-planet model TTVs for each of
the four phases (Table 21); from there we can see correlations between phases
1 & 3 and phases 2 & 4\. However, after subtracting the 2-planet and 3-planet
models from the respective TTVs, these correlations disappeared. Thus, no such
“modulo 4” pattern is observed in the AU Mic (catalog ) b transit times
(Figure 10), with transits being neither systematically late nor early in
phase with the stellar longitude crossings. The three Spitzer transits are
from three different sets of stellar longitude crossings, and one of the R-M
observation midpoint times is from the fourth set of stellar longitude
crossings; all four transit-midpoint times are systematically late w/r/t to
the TESS transit times by $\sim$5-10 minutes. The TESS transits themselves
also encompass three of the four stellar longitude crossing sets, and the
variation in transit times between them is less than a couple of minutes.

Third, Gilbert et al. (2021), in deriving the TESS transit midpoint times,
includes the joint modeling of the rotational modulation of the starspots in
deriving the transit midpoint times. This analysis approach thus mitigates any
impact the spot modulation has on the values of the derived transit midpoint
times as estimated in Oshagh et al. (2013). Additionally, because of the joint
modeling of the TESS transits, the impact of the rotational modulation is
reflected in the posterior uncertainties and thus precision in the derived
transit times. Our photodynamical analysis has shown that when modeling the
rotational modulation of starspots with a GP in parallel with deriving transit
midpoint times, as opposed to in series, the GP can absorb rather than
introduce additional TTVs on the order of a few minutes for AU Mic (catalog ).

Thus, in light of the modeling presented herein and Gilbert et al. (2021) the
rotational modulation of starspots can likely be ruled out as an explanation
for our observed TTVs between TESS and Spitzer. More detailed synthetic
simulations will be needed in the future to further quantify the impact of the
rotational modulation of AU Mic (catalog ) starspots on the derived transit
times, but that is beyond the scope of this work. We next turn to assess the
impact of flares on the derived transit times.

###  6.3 Stellar Flares

Like stellar spots, stellar flares can impact the midpoint timing of the
transits, particularly when they occur during ingress or egress, which
happened several times during the TESS transits. For the ground-based transit
data, it is difficult to resolve the flares due to the relatively lower
photometric precision, and additionally AIJ and EXOFASTv2 do not include joint
modeling of flares during transits. The ground-based transit times also
possess larger timing uncertainties, which may mask the impact of flares.

Several flares occurred during the TESS observations, but Gilbert et al.
(2021) jointly modeled them with the bayesflare (Pitkin et al., 2014) and
xoflares packages, some of which do occur during egress and thus could impact
the derived transit midpoint times were they not modeled jointly.
Additionally, the O–C diagram from our photodynamical analysis demonstrates
that the derived transit times are sensitive to the methods employed for
accounting for the stellar flares. However, the analysis done by Gilbert et
al. (2021) shows no dependence of transit timing on activity after
marginalizing over models for the flares and spot modulation, and the TESS
transit times of AU Mic (catalog ) b are fairly constant to within $\sim$4
minutes.

The second partial Spitzer transit of AU Mic (catalog ) b contains an obvious
flare during transit, which we marginalize over in our modeling, thus
impacting our derived timing uncertainties. Additionally, the first and and
third Spitzer transits may show some small flares pre-ingress which we do not
account for in our modeling. As with spots, however, flares are smaller in
amplitude at 4.5 $\mu$m than at visible wavelengths and thus should have a
relatively smaller impact on the derived transit times. However, we find that
the first two Spitzer transits times are consistent with one another, and the
third Spitzer transit is most contemporaneous with the R-M transit
observations.

For the SPIRou \+ iSHELL R-M transit observation, the analysis took into
account the impact of magnetic activity and flares by constructing a model
using a similar approach from Donati et al. (2020), which is then subtracted
from the RV observations. Similarly, for the ESPRESSO data, the stellar
activity was modeled using the celerite package’s GP that is described by a
Matérn 3/2 kernel and then subtracted from the RV data.

All five R-M and Spitzer transits show significantly deviant and late transit
times; it would be difficult and randomly unlucky to have randomly timed
flares during the different transits all impact the derived transit midpoint
times in the same way – late, as opposed to early. While these results are
based on a relatively small number of transits, it seems unlikely that flares
can account for the observed TTV excesses. More detailed simulations in the
future will be needed to assess the impact of flares on the transit times of
AU Mic (catalog ), but that is beyond the scope of this work. We next turn to
consider the possibility of additional planets in the AU Mic (catalog )
system.

###  6.4 Existence of a planet d candidate?

In the previous subsections, we explored the stellar activity of AU Mic
(catalog ) and the significance of its impact on the TTVs, both through the
rotational modulation of starspots and flares. It is possible that there is
some unaccounted-for effect in the derived TTV uncertainties, but we deem this
scenario unlikely given the above marginalization over our activity models in
the derived transit time posteriors. Thus, we turn to another possibility: the
presence of a third planet.

Cale et al. (2021) explored additional candidate RV planet signals when
modeling the RVs of the AU Mic (catalog ) system. One candidate period that
was explored in particular was the presence of a planet in-between b and c, a
“middle-d” non-transiting planet with a period of 12.742 days. Such a planet
would establish that the AU Mic (catalog ) system is in a 4:6:9 orbital
commensurability and result in significant TTVs. In $\S$4.2.2, we modeled this
3-planet configuration with Exo-Striker based on this candidate RV signal from
Cale et al. (2021). Since the AMD criterion indicated that this configuration
is unstable, we tested its stability with N-body packages, including rebound
(Rein & Liu, 2012; Rein & Spiegel, 2015) and mercury6 (Chambers, 1999), the
latter of which was used as a consistency check. The mercury6 simulation ran
for 10 Myrs while the rebound simulation ran for 2 Myrs; both indicated that
the 3-planet configuration from $\S$4.2.2 is stable (see Figure 11 for outputs
from rebound). Additional possible configurations, including having a planet
beyond c, is beyond the scope of this work. While we find that our TTVs are
consistent with a possible third “middle d” transiting planet, we neither
confirm nor rule out its existence. If this third planet does exist, our
dynamical modeling implies the impact on the observed TTVs will be more
readily apparent with additional TTV measurements in 2021 and the next several
years, as the TTVs for AU Mic (catalog ) b would deviate further from a linear
ephemerides, and given the “curvature” in the TTVs from the first three years
of TTVs presented herein.

Figure 11: rebound model of the stability of 3-planet system for AU Mic
(catalog ). The top three plots span 2 Myrs while the bottom three plots are
zoomed in between 0.95 and 1.05 Myrs.

###  6.5 Implication on the Mass of AU Mic c

With only three transits of AU Mic (catalog ) c, all with relatively larger
transit timing uncertainties, we do not derive a TTV mass constraint for AU
Mic (catalog ) b. However, from the perturbations of the transit times of AU
Mic (catalog ) b, we are able to place constraints on the mass of AU Mic
(catalog ) c. Table 22 summarizes our Exo-Striker and photodynamical TTV
masses with the RV masses from Cale et al. (2021). In both our 2-planet and
3-planet Exo-Striker models, the mass of AU Mic (catalog ) c is approximately
that of Neptune. Both photodynamical models determine the mass of AU Mic
(catalog ) c to be more comparable to that of Uranus instead of the more
massive Neptune. The densities of AU Mic (catalog ) c using the masses from
Exo-Striker analyses are quite high. We instead adopt as our final mass
estimate for AU Mic (catalog ) c of $10.8^{+2.3}_{-2.2}$ M⊕, the value from
our “a” photodynamical model which jointly models the stellar activity with
the transit dynamics with no third planet. This value is consistent with the
mass for AU Mic (catalog ) c as determined from RVs (as are all the masses we
derive). Note that we only list the masses of AU Mic (catalog ) c; we
currently do not have enough transits nor enough precision on the timing of
transits for AU Mic (catalog ) c to place a meaningful constraint on the mass
of AU Mic (catalog ) b.

Table 22: Characteristics of AU Mic (catalog ) c. Property | Unit | Quantity | Ref
---|---|---|---
Mass | M⊕ | 9.60${}^{+2.07}_{-2.31}$ | KJ1 (Cale et al., 2021)
14.12${}^{+2.48}_{-2.71}$ | KJ2 (Cale et al., 2021)
24.8 $\pm$ 1.2 |  Exo-Striker 2-planet Model
20.1 $\pm$ 1.6 |  Exo-Striker 3-planet Model
10.8${}^{+2.3}_{-2.2}$ | Photodynamical Model a
13${}^{+5}_{-2}$ | Photodynamical Model b
Radius | R⊕ | 2.79${}^{+0.31}_{-0.30}$ | Gilbert et al. (2021)
Density | g/cm3 | 2.44${}^{+0.97}_{-0.98}$ | KJ1 (Cale et al., 2021)
3.58 $\pm$ 1.35 | KJ2 (Cale et al., 2021)
6.28${}^{+2.12}_{-2.05}$ |  Exo-Striker 2-planet Model
5.08${}^{+1.74}_{-1.69}$ |  Exo-Striker 3-planet Model
2.75${}^{+1.09}_{-1.05}$ | Photodynamical Model a
3.3${}^{+1.7}_{-1.2}$ | Photodynamical Model b

###  6.6 Limitation of Our Joint Modeling & Future Work

In this subsection, we discuss some of the limitations of our analysis. First,
we do not perform joint RV + TTV modeling because the RVs require a custom
treatment of stellar activity that is an order of magnitude larger in
effective amplitude relative to the Keplerian orbital reflex motion of the
star (Cale et al., 2021). We also sought to conduct an independent TTV
modeling to compare with the RV analysis, given the complexity of the RV
model. Additionally, to our knowledge, no prior work exists in performing a
joint RV + TTV modeling in the presence of such stellar activity as exhibited
by AU Mic (catalog ).

Second, in RVs it is common to adopt a “jitter” term representing the white-
noise measurement error that is not captured by the formal measurement
uncertainties and added in quadrature, without which Keplerian orbital fits
can yield reduced-$\chi^{2}\gg 1$. A similar approach may be employed for AU
Mic (catalog ) in a future study for the TTV modeling – including a “jitter”
term that accounts for timing uncertainties not captured by our modeling of
spots, flares and orbital dynamics. We do not undertake such an analysis
herein for several reasons. First, Exo-Striker does not include a “jitter” TTV
error term, and thus we cannot marginalize over this parameter; e.g., we would
need to assume a “jitter” timing term and inflate our timing measurement
uncertainties. Second, our timing data is fairly in-homogeneous, spanning a
range of precision and wavelength from ground and space facilities; we would
need multiple independent jitter terms that are not well constrained; our
Spitzer data that most deviates from a two-planet TTV model is least impacted
by stellar activity. Third, as shown in Szabó et al. (2021), due to the 7:4
period commensurability of the stellar rotation to the orbital period of AU
Mic (catalog ) b, the TTV “jitter” from spots will not be a white noise term;
given the known stellar rotation period, the timing impact is not simply
random and fairly deterministic, as we have undertaken our “mod-4” transit
analysis to identify whether or not the transit times correlation with the
stellar longitude crossings. That being said, the flares are random in time
and will randomly impact the measured transit times. Therefore, a “jitter”
term would be appropriate for accounting for the impact of flares on measured
the TTVs that are not captured by our flare modeling of the light curves.

Third, the uncertainties in the midpoint times from the R-M observations are
likely to be underestimated. As allowing certain parameters, such as the
orbital period and eccentricity, to remain free will result in degeneracy in
the R-M model, setting such parameters fixed would avoid this issue but would
cause the uncertainties to be underestimated.

As a consequence of these limitations in our analysis, we do not perform any
statistically robust model comparisons between the 2- and 3-planet TTV models.
Future observations and analyses will be necessary to reach a definitive
conclusion for the hypothetical third planet candidate.

##  7 Conclusion

AU Mic (catalog ) hosts a young nearby exoplanet system that serves as a
useful laboratory for probing and characterizing young exoplanetary systems.
We have collected 23 transits from AU Mic (catalog ) b and 3 transits from AU
Mic (catalog ) c over the course of three years. We model the observed
transits and derive new transit model posteriors. We have run two independent
methods (Exo-Striker and photodynamical) in modeling the transits and the
timing of those transits. Our observations and analyses of the transits of AU
Mic (catalog ) b and c show that AU Mic (catalog ) b is exhibiting TTVs,
consistent with Szabó et al. (2021), Martioli et al. (2021), and Gilbert et
al. (2021). Our photodynamical model yields the mass for AU Mic (catalog ) c
of 10.8${}^{+2.3}_{-2.2}$ M⊕, consistent with the RV-determined mass in Cale
et al. (2021). However, going beyond the work of Szabó et al. (2021) and
Martioli et al. (2021), our TTVs show timing excess of $>$5 minutes that
appear to discrepant with a 2-planet model alone, particularly when comparing
the Spitzer and TESS derived transit times, the former of which is presented
for the first time herein. Further, we marginalize our TTV models over our
models for the rotational modulation of stellar activity and flares.
Consequently, stellar activity, while not excluded through statistically
robust model comparison, does not appear likely to be able to account for the
observed TTV excess. We mapped the O–C diagram taking into account AU Mic
(catalog )’s stellar spin 7:4 commensurability with the orbital period of AU
Mic (catalog ) b from Szabó et al. (2021), and we did not identify that the
spot modulation results in a significant effect on the observed TTVs as a
function of stellar longitude crossings. We explore a possible,
representative, and non-exhaustive 3-planet configuration scenario that is:
consistent with the identified non-transiting “middle-d” RV candidate signal
in Cale et al. (2021), is dynamically stable, would establish the AU Mic
(catalog ) system of planets in a 4:6:9 orbital period commensurability, and
can account for the observed TTV excess of our data.

Nonetheless, given the high level of stellar activity for AU Mic (catalog ),
we cannot ignore the possibility that our modeling and marginalization over
our stellar activity models does not fully account for some effects that
significantly impact the observed TTVs. Thus, additional ground- and space-
based TTVs in the next few years are needed to further vet the impact of
stellar activity on the observed TTVs, or to confirm the possibility that the
excess TTVs are due to the RV candidate highlighted in Cale et al. (2021), and
to enable a more thorough search of the possible orbital periods for possible