Diving Beneath the Sea of Stellar Activity: Chromatic Radial Velocities of the Young AU Mic Planetary System

Our analyses make use of new and archival high-resolution echelle spectra from a variety of facilities, which are summarized in Table 1. We briefly detail new spectroscopic observations and the corresponding RVs from observing programs primarily intended to characterize the AU Mic planetary system.

2.1.1. CARMENES

2.1.2. CHIRON

2.1.3. HIRES

2.1.4. Minerva-Australis

2.1.5. TRES

2.1.6. IRD

2.1.7. iSHELL

The NASA TESS mission (Ricker et al. 2015) observed AU Mic in Sectors 1 (2018 July 25–August 22) and 27 (2020 July 4–30). We download the light curves from the Mikulski Archive for Space Telescopes (MAST; Swade et al. 2018). We use the Science Processing Operations Center (SPOC; Jenkins et al. 2016) "Presearch Data Conditioning" light curves utilizing "Simple Aperture Photometry" (PDCSAP; Smith et al. 2012; Stumpe et al. 2012, 2014) to inform our model in Section 3.3.1.

We primarily seek to utilize a global (joint) Gaussian process model with multiple realizations that give rise to the data we observe with all of the above instruments simultaneously. To implement our desired framework, we have developed two Python packages. We leave the description of optimize—a high-level Bayesian inference framework to Appendix B.

To provide RV-specific routines, we extend the optimize package within the orbits submodule of the pychell (Cale et al. 2019) package. ⁴⁹ We define classes specific for RV data, models, and likelihoods, with much of the "boilerplate" code handled through optimize. A top-level "RVProblem" further defines a pool of RV-specific methods for pre- and post-optimization routines, such as plotting phased RVs, periodogram tools, model comparison tests, and propagation of Markov-Chain Monte Carlo (MCMC) chains for deterministic Keplerian parameters (e.g., planet masses, semimajor axes, and densities).

A Gaussian process kernel is defined through a square matrix, K (also called the covariance matrix), where each entry describes the covariance between two measurements. ⁵⁰ We introduce two GP kernels as extensions of the quasi-periodic (QP) kernel, which has been demonstrated in numerous cases to model rotationally modulated stellar activity in both photometric and RV observations (see Section 1): ⁵¹

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{K}}}_{\mathrm{QP}}({t}_{i},{t}_{j}) & = & {\eta }_{\sigma }^{2}\exp \left[-\displaystyle \frac{{\rm{\Delta }}{t}^{2}}{2{\eta }_{\tau }^{2}}-\displaystyle \frac{1}{2{\eta }_{{\ell }}^{2}}{\sin }^{2}\left(\pi \displaystyle \frac{{\rm{\Delta }}t}{{\eta }_{p}}\right)\right]\\ \mathrm{where}\ \ {\rm{\Delta }}t & = & | {t}_{i}-{t}_{j}| .\end{array}\end{eqnarray} \tag{ 1 }$$

Here, η_P typically represents the stellar-rotation period; η_τ the mean spot lifetime; η_ℓ is the relative contribution of the periodic term, which may be interpreted as a smoothing parameter (larger is smoother); and η_σ is the amplitude of the autocorrelation of the activity signal.

We seek to use a fully inclusive QP-like kernel that accounts for the wavelength-dependence of the stellar activity present in our multiwavelength data set. In this work, we only modify the amplitude parameter, η_σ ; we leave further chromatic modifications (namely, convective blueshift and limb-darkening; see Section 1) as subjects for future work. To first order, we expect the amplitude from activity to be linearly proportional to frequency (or inversely proportional to wavelength). This approximation is a direct result of the spot-contrast scaling with the photon frequency (or inversely with wavelength) from the ratio of two blackbody functions with different effective temperatures (Reiners et al. 2010).

We first reparameterize the amplitude through a linear kernel as follows:

$$\begin{eqnarray}&&{{\boldsymbol{K}}}_{{\boldsymbol{J}}1}({t}_{i},{t}_{j})={\eta }_{\sigma ,{\rm{s}}(i)}{\eta }_{\sigma ,{\rm{s}}(j)}\times \exp [...].\end{eqnarray} \tag{ 2 }$$

Here, η_σ,s(i) η_σ,s(j) are the effective amplitudes for the spectrographs at times t_i and t_j , respectively, where s(i) represents an indexing set between the observations at time t_i and spectrograph s. ⁵² Each amplitude is a free parameter.

We also consider a variation of this kernel that further enforces the expected inverse relationship between the amplitude with wavelength. We rewrite the kernel to become:

$$\begin{eqnarray}&&{{\boldsymbol{K}}}_{{\boldsymbol{J}}2}({t}_{i},{t}_{j},{\lambda }_{i},{\lambda }_{j})={\eta }_{\sigma ,0}^{2}{\left(\displaystyle \frac{{\lambda }_{0}}{\sqrt{{\lambda }_{i}{\lambda }_{j}}}\right)}^{2{\eta }_{\lambda }}\times \exp [...].\end{eqnarray} \tag{ 3 }$$

Here, η_σ,0 is the effective amplitude at λ = λ₀, and η_λ is an additional power-law scaling parameter with wavelength to allow for a more flexible nonlinear (with frequency) relation. λ_i and λ_j are the "effective" wavelengths for observations at times t_i and t_j , respectively. For both Equations (2) and (3), the expression within square brackets is identical to that in Equation (1).

To make predictions from K _{J 2} (Equation (3)), we follow Rasmussen & Williams (2006; their Equations (2.23) and (2.24)). We construct the matrix K _{J 2}(t_i,*, t_j , λ_*, λ_j ), which denotes the n_* × n matrix of the covariances evaluated at all pairs of test points and training points (the data). Wavelengths in the * dimension are identical, and therefore each realization corresponds to a unique wavelength. This formulation allows us to realize the GP with high accuracy for all wavelengths so long as at least one wavelength is sampled near t_i,*. Predictions with kernel K _{J 1} (Equation (2)) are found in a similar fashion, where each realization corresponds to a particular spectrograph.

We first bin out-of-transit RV observations from each night (per spectrograph). While not negligible, we expect changes from rotationally modulated activity to be small within a night, so we choose to mitigate activity on shorter timescales our model is not intended to capture (e.g., p-mode oscillations, granulation). The median RV for each spectrograph is also subtracted. We choose to ignore poorly sampled regions with respect to our adopted mean spot lifetime η_τ (100 days; see Section 3.3.1). Each instance of a covariance matrix represents a family of functions, and therefore the GP regression may be too flexible (and thus poorly constrained) in regions of low-cadence observations. We also ignore regions with only low-precision measurements (median errors ≳10 m s⁻¹). This limits our analyses to all observations between 2019 September and 2020 December, and the spectrographs HIRES, TRES, CARMENES-VIS and NIR, SPIRou, and iSHELL. We do not include six binned IRD or 13 binned Minerva-Australis observations in our primary analyses, as we expect the offsets to be poorly constrained in the presence of stellar activity. Finally, we discard three CARMENES-VIS and 13 CARMENES-NIR measurements from our analyses, primarily near the beginning of each season, due to residuals >100 m s⁻¹ that are inconsistent with our other data sets. We suspect that telluric contamination, which is further exacerbated by the high airmass of the observations, may have degraded the CARMENES observations. For completeness, we present fit results including all spectrographs in Appendix D. A summary of measurements is provided in Table 1.

Our RV model first consists of two Keplerian components for the known transiting planets, a GP model for stellar activity, and per-instrument zero points. The zero points are each assigned to 1 m s⁻¹ with a uniform prior of ±300 m s⁻¹. We further adopt a normal prior of ${ \mathcal N }(0,100)$ to make each offset well-behaved. When using multiple priors, the composite prior probability for such a parameter will not integrate to unity. For the combination of a uniform + normal prior, this is not a concern; the normal prior is properly normalized and takes on a continuous range of values, whereas the uniform prior will either result in a constant term added to the likelihood function if the parameter is in bounds, or − ∞ if out of bounds.

Analyses of the TESS transits in M21 found P_b = 8.4629991 ± 0.0000024 days, TC_b = 2458330.39046 ±0.00016, P_c = 18.858991 ± 0.00001 days, and TC_c =2458342.2243 ± 0.0003. For all of our analyses, we fix P and TC for planets b and c; the uncertainties in these measurements are insignificant even for our full baseline of ≈17 yr. The semi-amplitudes of each planet start at K_b = 8.5 m s⁻¹ and K_c = 5 m s⁻¹, and are only enforced to be positive. Preliminary analyses of a secondary eclipse observed in Spitzer observations support a moderately eccentric orbit for AU Mic b, with e_b = 0.189 ± 0.04 (K. I. Collins et al. 2021, in preparation), which is somewhat larger than the eccentricity determined from the duration of the primary transits observed with TESS (e_b = 0.12 ± 0.04; Gilbert et al. 2021). We assume a circular orbit for AU Mic c, and further examine eccentric cases in Section 5.1. The Keplerian component of our RV model in pychell is nearly identical to that used in RadVel (Fulton et al. 2018). Kepler's equation is written in Python and makes use of the numba.@njit decorator (Lam et al. 2015) for optimal performance. We exclusively use the orbit basis {P, TC, e, ω, K}.

Our optimizer seeks to maximize the natural logarithm of the a posteriori probability (MAP) under the assumption of normally distributed errors:

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }=-\displaystyle \frac{1}{2}\left[{{\boldsymbol{r}}}^{{\mathsf{T}}}{{\boldsymbol{K}}}_{{\boldsymbol{o}}}^{-1}{\boldsymbol{r}}+\mathrm{ln}| {{\boldsymbol{K}}}_{{\boldsymbol{o}}}| +N\mathrm{ln}(2\pi )\right]+\sum _{i}\mathrm{ln}{{ \mathcal P }}_{i}.\end{eqnarray} \tag{ 4 }$$

Here, r is the vector of residuals between the observations and model, K _o is the covariance matrix sampled at the same observations, N is the number of data points, and $\{{{ \mathcal P }}_{i}\}$ is the set of prior knowledge. We maximize Equation (4) using the iterative Nelder–Mead algorithm described in Cale et al. (2019), which is included as part of the optimize package. We also sample the posterior distributions using the emcee package (Foreman-Mackey et al. 2013) for a subset of models to determine parameter uncertainties, always starting from the MAP-derived parameters. In all cases, we use twice the number of chains as varied parameters. We perform a burn-in phase of 1000 steps followed by a full MCMC analysis for ≈50× the median autocorrelation time (steps) of all chains.

3.3.1. Estimation of Kernel Parameters

We briefly analyze both sectors of TESS photometry in order to estimate the GP kernel parameters η_τ , η_ℓ , and η_P . We note that the rotationally modulated structure in both sectors is consistent (Figure 1). If we assume spots are spatially static in the rest frame of the stellar surface (i.e., spots do not migrate), this suggests a similar spot configuration and contrast for each sector. We first determine η_P by qualitatively analyzing both TESS sectors phased up to periods close to that used in M21 (4.862 ± 0.032 days) with a step size of 0.001 day (see Figure 1). We find η_P ≈ 4.836 or η_P ≈ 4.869 days from our range of periods tested; no periods between these two values are consistent with our assumption of an identical spot configuration. The difference in these two periods further corresponds to one additional period between the two sectors (i.e., ∣1/η_P,1 −1/η_P,2∣ ≈ 1/700 day⁻¹). The smaller of these two values implies AU Mic b is in a 7:4 resonance with the stellar rotation period (Szabó et al. 2021), potentially indicating tidal interactions between the planet and star. We adopt ${\eta }_{P}\sim { \mathcal N }(4.836,0.001)$ in all our analyses where the uncertainty is a conservative estimate determined by our step size.

Zoom In Zoom Out Reset image size

Although the TESS light curve itself can provide insight into η_τ and η_ℓ , we instead try to estimate these values directly from the predicted spot-induced RV variability via the ${FF}^{\prime}$ technique (Aigrain et al. 2012):

$$\begin{eqnarray}&&{\rm{\Delta }}{\mathrm{RV}}_{\mathrm{spots}}(t)=-F(t)F^{\prime} (t){R}_{\star }/f.\end{eqnarray} \tag{ 5 }$$

Here, F is the photometric flux and f represents the relative flux drop for a spot at the center of the stellar disk. To compute F and $F^{\prime}$ (the derivative of F with respect to time), we first fit the TESS light curve via cubic spline regression (scipy.interpolate.LSQUnivariateSpline; Virtanen et al. 2020) for each sector individually, with knots sampled in units of 0.5 day (≈10% of one rotation period), to average over transits and the majority of flare events (Figure 1). The nominal cubic splines are then used to directly compute both F and $F^{\prime}$ on a downsampled grid of 100 evenly spaced points for each sector. We then divide the resulting joint-sector curve by its standard deviation for normalization; we do not care to directly fit for the chromatic parameter f(TESS, ...). We further assume f to be constant in time (i.e., spots are well-dispersed on the stellar surface). We then perform both MAP and MCMC analyses for this curve using a standard QP kernel (Equation (1)) with loose uniform priors of ${\eta }_{\tau }\sim { \mathcal U }(10,2000)$ (days) and ${\eta }_{{\ell }}\sim { \mathcal U }(0.05,0.6)$. We set the intrinsic error bars of the curve to zero but include an additional "jitter" (white noise) term in the model with a Jeffreys prior (Jeffreys 1946) distribution with the knee at zero to help keep the jitter well-behaved by discouraging larger values unless it significantly improves the fit quality through an inversely proportional penalty term. The amplitude of the model is drawn from a wide uniform distribution of ${ \mathcal U }(0.3,3.0)$. The posterior distributions are provided in Figure 2.

Zoom In Zoom Out Reset image size

A fit to the ${FF}^{\prime}$ curve suggests a mean activity timescale ${\eta }_{\tau }\approx {92}_{-23}^{+29}$ days. Although our interpretation implies η_τ should be comparable to the gap between the two sectors, (∼700 days), we do not have photometric measurements between the two sectors, and therefore cannot speak to evolution that will be important for our 2019 observations. We further note that the TESS Sector 27 light curve exhibits moderate evolution whereas Sector 1 appears more stable (Figure 1). The posterior distributions are also consistent with a relatively smooth GP with the period length scale η_ℓ ≈ 0.45 ± 0.06.

Before making use of our joint kernels, we first assess the performance of the standard QP kernel (Equation (1)) for each instrument individually. Here, each spectrograph makes use of a unique QP kernel and amplitude term, but the remaining three GP parameters are shared across all kernels. Each amplitude is drawn from a normal prior with mean equal to the standard deviation of the data set, and a conservative width of 30 m s⁻¹. The expected semi-amplitudes for AU Mic b and c (≲10 m s⁻¹) will negligibly affect this estimation. We also apply a Jeffreys prior with the knee at zero to help keep the amplitude well-behaved. ⁵³ For η_τ and η_ℓ , we first make use of the same priors used to model the ${FF}^{\prime}$ curve. We further include a fixed jitter term at 3 m s⁻¹ added in quadrature along the diagonal of the covariance matrix K _o for the HIRES observations only; HIRES observations provide the smallest intrinsic uncertainties, but are most impacted by activity (largest in amplitude), so we choose to moderately downweight the HIRES observations. Given the flexibility of GP regression with a nightly cadence, we choose not to fit for (nor include) jitter terms for other spectrographs, and further discuss this decision in Section 5.5. This is the most flexible model we employ to the RVs, and we therefore use these results to flag the aforementioned CARMENES-VIS and CARMENES-NIR measurements.

We find normally distributed posteriors for η_τ and η_ℓ (Figure 14), but the reduced χ² statistic of 0.32 indicates the model over fits the data. The per-spectrograph amplitudes are reasonably consistent with their respective priors, so we assert this is a result of η_τ (≈43 days) and/or η_ℓ (≈0.23) taking on too small of values, indicating our RV model is insufficient to constrain these values from the RV observations, due to insufficient cadence and/or an inadequate model. We therefore again fix η_τ = 100 days to let each season have mostly distinct activity models, while minimizing the flexibility within each season, which is consistent with what the ${FF}^{\prime}$ curve suggests. As a compromise between the ${FF}^{\prime}$ and RV analyses, we also fix η_ℓ = 0.28. Our adopted value of η_τ is larger than that used in K21 (≈70 days ⁵⁴ ), while η_ℓ is nearly identical. We further explore these decisions and their impact on our derived semi-amplitudes in Section 5.2. With fixed value for η_τ , we rerun MAP and MCMC fits with disjoint kernels, yielding a reduced χ² of 0.86, indicating the model is now only slightly overfit.

3.3.2. Joint-kernel RV Fitting

We use results from the disjoint case to inform our primary joint-kernel models. Although the different GPs appear similar (Figure 3), each still exhibits unique features, suggesting a simple scaling is not valid and/or insufficient sampling for each kernel individually. Regardless, our two joint kernels will enforce a perfect scaling between any two spectrographs. The RVs are shown in Figures 4 and 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We run MAP and MCMC fits using the joint kernel K _{J 1} (Equation (2)) again making use of the same normal and Jeffreys priors for each amplitude. We then fit the resulting set of best-fit amplitudes using our proposed power-law relation (see Equation (3)): ${\eta }_{\sigma }{(\lambda )={\eta }_{\sigma ,0}({\lambda }_{0}/\lambda )}^{{\eta }_{\lambda }}$ with scipy.optimize.curve_fit (Jones et al. 2001; Figure 6). We arbitrarily anchor λ₀ at λ = 565 nm. The effective mean wavelength of each spectrograph should consider the RV information content (stellar and calibration), and ignore regions with dense telluric features. For gas-cell-calibrated spectrographs (HIRES, CHIRON, and iSHELL), we do limit the the range to regions with gas-cell features. For all other spectrographs, we take the effective RV information content to be uniform over the full spectral range as a "zeroth-order" approximation (Reiners & Zechmeister 2020). We further do not consider regions of tellurics that may have been masked (e.g., CARMENES RVs generated with serval). Although these estimations are imperfect, they are only relevant to kernel K _{J 2} (Equation (3)). The adopted wavelengths for each spectrograph are listed in Table 1.

Zoom In Zoom Out Reset image size

We find η_σ,0 ≈ 221 m s⁻¹, and η_λ ≈ 1.17. This amplitude is significantly larger than the intrinsic scatter of our observations (namely HIRES) suggests, so we adopt a tight normal prior of ${ \mathcal N }(221,10)$ to restrict it from getting any larger. We only apply a loose uniform prior for ${\eta }_{\lambda }\sim { \mathcal U }(0.2,2)$. We then run corresponding MAP and MCMC fits with kernel K _{J 2} (Equation (3)). A summary of all parameters is provided in Table 2. We present and discuss fit results from both joint kernels in Section 4.

Table 2. The Model Parameters and Prior Distributions Used in Our Primary Fitting Routines

Parameter [units]	Initial Value (P0)	Priors	Citation
Pb [days]	8.4629991		Primary transit; M21
TCb [days]	2458330.39046		Primary transit;M21
eb	0.189	${ \mathcal N }({P}_{0},0.04)$	Secondary eclipse; K. L. Collins et al. (2021, in preparation)
ωb [rad]	1.5449655	${ \mathcal N }({P}_{0},0.004)$	Secondary eclipse; K. L. Collins et al. (2021, in preparation)
Kb [m s−1]	8.5	Positive	K21
Pc [days]	18.858991		Primary transit; M21
TCc [days]	2458342.2243		Primary transit; M21
ec	0		⋯
ωc [rad]	π		⋯
Kc [m s−1]	5	Positive	M21
ησ,0 [m s−1]	216	${ \mathcal J }(1,600)$, ${ \mathcal N }({P}_{0},10)$	RVs; this work
ηλ	1.18	${ \mathcal U }(0.3,2)$	RVs; this work
${\eta }_{\sigma ,\mathrm{HIRES}}$ [m s−1]	130	${ \mathcal J }(1,600)$, ${ \mathcal N }({P}_{0},30)$	RVs; this work
${\eta }_{\sigma ,\mathrm{TRES}}$ [m s−1]	103	${ \mathcal J }(1,600)$, ${ \mathcal N }({P}_{0},30)$	RVs; this work
ησ,CARM‐VIS [m s−1]	98	${ \mathcal J }(1,600)$, ${ \mathcal N }({P}_{0},30)$	RVs; this work
ησ,CARM‐NIR [m s−1]	80	${ \mathcal J }(1,600)$, ${ \mathcal N }({P}_{0},30)$	RVs; this work
ησ,SPIRou [m s−1]	42	${ \mathcal J }(1,600)$, ${ \mathcal N }({P}_{0},30)$	RVs; this work
ησ,iSHELL [m s−1]	40	${ \mathcal J }(1,600)$, ${ \mathcal N }({P}_{0},30)$	RVs; this work
ητ [days]	100		TESS light curve and RVs; this work
ηℓ	0.28	⋯	TESS light curve and RVs; this work
ηp [days]	4.836	${ \mathcal N }({P}_{0},0.001)$	TESS light curve; this work
γ (per spectrograph) [m s−1]	1	${ \mathcal U }(-300,300)$, ${ \mathcal N }(0,100)$	RVs; this work
${\sigma }_{\mathrm{HIRES}}$ [m s−1]	3		⋯
${\sigma }_{\mathrm{TRES}}$ [m s−1]	0		⋯
σCARM‐VIS [m s−1]	0		⋯
σCARM‐NIR [m s−1]	0		⋯
σSPIRou [m s−1]	0		⋯
σiSHELL [m s−1]	0		⋯
M⋆ [M⊙]	${0.5}_{-0.03}^{+0.03}$	⋯	P20
Rb [R⊕]	${4.38}_{-0.18}^{+0.18}$	⋯	P20
Rc [R⊕]	${3.51}_{-0.16}^{+0.16}$	⋯	M21

Note. Here, indicates the parameter is fixed. We run models utilizing K _{J 1} and K _{J 2}. We list the radii of AU Mic b and c measured in M21 that we use to compute the corresponding densities of each planet.

Download table as:  ASCIITypeset image

We compute periodograms to further assess the relative statistical confidence of the two transiting planets and to search for other planets in the system. We first compute a series of generalized Lomb–Scargle (GLS; Zechmeister & Kürster 2018; Czesla et al. 2019) periodograms out to 500 days after removing the nominal zero points, appropriate GPs, and the two planets, all generated using parameters from our nominal two-planet model (Table 3) with kernel K _{J 1} (Equation (2)) to model the stellar activity. We also compute an activity-filtered periodogram from a planet-free model to assess how much the GP model will absorb planetary signals, and inform our interpretation of other peaks present in the periodogram. We further plot the normalized power levels for false-alarm probabilities (FAPs) of 10%, 1%, and 0.1%.

We also compute "brute-force" periodograms by performing MAP fits for a wide range of fixed orbital periods for a user-defined "test planet" with various assumptions for other model parameters (see Addison et al. 2021). Given the time complexity of GP regression, we only consider periods out to 100 days. We first run two searches with no other planets in the model, first allowing for the test planet's TC to float, and second fixing TC to the nominal value for AU Mic b (Table 3). We then run searches for a second planet, this time including a planetary model to account for the orbit of AU Mic b, with ${K}_{b}\sim { \mathcal N }(8.5,2.5)$, consistent with the semi-amplitude found in K21. We again consider the case of letting the test planet's TC float, then run three cases with fixing the test planet's TC to each time of transit for AU Mic c from the TESS Sector 1 and 27 light curves (M21). Last, we perform a search for a third planet, letting its TC float, and including models for AU Mic b and c (${K}_{b}\sim { \mathcal N }(8.5,2.5)$, K_c > 0).

Both the GLS (Figure 8) and brute-force (Figure 9) periodograms exhibit clear aliasing with a frequency of ≈0.00281 days⁻¹ (or 356 days), which we attribute to having two seasons of observations separated by ≈200 days. Given the respective power of AU Mic b in both the GLS and brute-force planet-free periodograms, we briefly explore other peaks with similar power, even though all other peaks are below all three FAPs after removing the nominal two-planet model (Figure 8, row 3; Figure 9, row 7). Both two- and zero-planet periodograms (as well as GLS and brute-force) show power between the orbits of AU Mic b and c near 12.72 and 13.19 days, as well as power near 66.7 days for the residual RVs. Although these peaks are comparable in power to AU Mic b in both planet-free periodograms, they may be spurious. We further discuss the confirmation of AU Mic b and c as well as the validation of such additional potential candidates in Section 5.3. A mass–radius diagram is shown in Figure 10 to place the mass and radius of all AU Mic b and c in context with other known exoplanets, including a subset of young sample of exoplanets shown in P20. The plotted masses for AU Mic b and c are from our nominal two-planet model using kernel K _{J 1} (Equation (2)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Here, we briefly explore eccentric orbits for the two-transiting planets b and c. For each planet, we take $e\sim { \mathcal U }(0,0.7)$ and $\omega \sim { \mathcal U }(0,2\pi )$. We only use kernel K _{J 1} (Equation (2)) to model the stellar activity. Posterior distributions are presented in Figure 18. We find e_b = 0.30 ± 0.04, which is ≈50% larger than our prior informed by a secondary eclipse event indicates. The corresponding finding of ω_b = 3.01 ± 0.27 is also inconsistent with our adopted prior for ω_b . The posterior distribution for e_c is concentrated at the upper bound (0.7), implying an overlapping orbit with AU Mic b. Orbital stability calculations presented in M21 indicate e_c < 0.2, so we assert our model is unable to accurately constrain its eccentricity. The behavior of e_c further indicates our detection of K_c may not be significant.

Our analyses in Section 3.3.1 make use of a fixed mean spot lifetime η_τ = 100 days and smoothing parameter η_ℓ = 0.28. Here, we determine how sensitive the recovered semi-amplitudes of AU Mic b and c are to these two parameters. We consider η_τ ∈ {40, 70, 100, 200, 300} (days), and η_ℓ ∈ {0.15, 0.2, 0.25, 0.3, 0.35}. We perform MAP and MCMC fits for all pairs of these two fixed parameters using K _{J 1} (Equation (2)) for a two-planet model. All other parameters adopt initial values and priors from Table 2. Results are summarized in Table 6.

Table 6. MCMC Results with Different Assumptions for the Mean Spot Lifetime η_τ and η_ℓ Using Kernel K _{J 1} (Equation (2))

ητ (days)	ηℓ	Kb (m s−1)	Kc (m s−1)	${\chi }_{\mathrm{red}}^{2}$
40	0.15	8.79 ± 1.47	7.38 ± 1.65	1.58
40	0.2	8.84 ± 1.30	8.51 ± 1.33	2.10
40	0.25	8.23 ± 1.17	9.05 ± 1.17	2.45
40	0.3	7.41 ± 1.08	9.13 ± 1.13	2.68
40	0.35	6.95 ± 0.98	9.23 ± 1.05	2.85
70	0.15	8.74 ± 1.24	6.88 ± 1.30	1.99
70	0.2	10.32 ± 1.13	5.90 ± 1.07	2.69
70	0.25	10.45 ± 1.04	4.76 ± 0.95	3.25
70	0.3	9.61 ± 0.91	4.16 ± 0.89	3.65
70	0.35	9.18 ± 0.82	3.94 ± 0.83	3.90
100	0.15	9.28 ± 1.17	5.88 ± 1.08	2.46
100	0.2	10.85 ± 1.00	4.73 ± 0.98	3.20
100	0.25	10.78 ± 0.95	3.78 ± 0.87	3.76
100	0.3	9.81 ± 0.85	3.63 ± 0.80	4.12
100	0.35	9.22 ± 0.80	3.60 ± 0.77	4.39
200	0.15	9.38 ± 0.98	4.01 ± 0.96	3.44
200	0.2	11.04 ± 0.89	3.35 ± 0.84	4.16
200	0.25	11.09 ± 0.87	3.38 ± 0.77	4.73
200	0.3	10.14 ± 0.84	4.51 ± 0.76	5.28
200	0.35	9.06 ± 0.73	4.77 ± 0.66	5.59
300	0.15	9.32 ± 0.92	3.84 ± 0.87	3.82
300	0.2	10.60 ± 0.88	3.78 ± 0.81	4.68
300	0.25	10.42 ± 0.80	4.99 ± 0.74	5.59
300	0.3	10.51 ± 0.76	4.52 ± 0.68	5.96
300	0.35	9.89 ± 0.75	4.41 ± 0.68	6.27

Note. For each row, we fix the values of η_τ and η_ℓ . All other model parameters take on the initial values and priors from Table 2 for a two-planet model. We perform a MAP fit followed by MCMC sampling for each case. We report the nominal values and uncertainties for the semi-amplitudes of AU Mic b and c from the MCMC fitting, as well as the reduced chi-squared statistic, ${\chi }_{\mathrm{red}}^{2}$, using the MAP-derived parameters. Uncertainties reported here for K_b and K_c are the average of the upper and lower uncertainties.

Download table as:  ASCIITypeset image

We find K_b is only moderately sensitive to the values of each hyperparameter, ranging from ∼7 to 11 m s⁻¹. With a larger spot lifetime, K_b tends toward larger values, indicating the GP is likely absorbing power from planet b with a more flexible model (smaller η_τ ). However, K_b is relatively insensitive to the value of η_ℓ . The range of values for K_c is larger, changing by nearly a factor of three. Unlike K_b , K_c is more unstable and tends toward larger values when using a more flexible (smaller) spot lifetime. The reduced chi-squared statistic indicates the model is not overfit in any of the cases performed, but is also larger than unity by a several factors in most cases, indicating our modeling is inadequate.

Here, we assess the fidelity of our RV model applied to the AU Mic system through planetary injection and recovery tests. We first inject planetary signals into the RV data with well-defined semi-amplitudes, periods, and ephemerides (TC). We arbitrarily choose TC = 2457147.36589 for all injected cases. We consider 40 unique periods between 5.12345 and 100.12345 days, uniformly distributed in log space. For the semi-amplitude K, we consider values from 1 to 10 m s⁻¹ with a step size of 1 m s⁻¹, as well as values between 10 and 100 m s⁻¹ that are uniformly distributed in log space (20 total values). In all cases, we include a model for AU Mic b with fixed P and TC such that ${K}_{b}\sim { \mathcal N }(8.5,2.5)$. We first assess our recovery capabilities using a Gaussian prior for P such that $P\sim { \mathcal N }({P}_{\mathrm{inj}},{P}_{\mathrm{inj}}/50)$ and a uniform prior for TC such that $\mathrm{TC}\sim { \mathcal U }({\mathrm{TC}}_{\mathrm{inj}}\pm {P}_{\mathrm{inj}}/2)$. For each injected planet (one at a time), we run our MAP and MCMC analyses to determine the recovered K and corresponding uncertainty. The starting values for K and TC are always the injected values. We also consider the same injection and recovery test but with P and TC fixed to the injected value. We finally determine how susceptible our RV model is to pick out "fake" planets by running these same two trials with no injected planets. Although there are no injected planets, we still run the same trials as for the injected case with different initial values for K. A two-dimensional histogram of the recovered K as a fraction of the injected K, as well as the associated uncertainty (also as a fraction of the injected K) for each case are shown in Figures 11 and 12 for the injected and non-injected cases, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the case of injected planets, we find our RV model is able to confidently recover semi-amplitudes down to a few m s⁻¹ in this data set with a relative precision of ≳4σ. However, a closer inspection reveals the recovered semi-amplitudes are typically larger than the injected K, in particular for smaller injected values (1–5 m s⁻¹) that include our measured semi-amplitude AU Mic c. When the ephemeris is known, we tend to poorly measure the smallest values of K, indicating the recovered TC in the non-fixed case is unlikely what we have injected. In the 6–13 m s⁻¹ range, which covers the recovered semi-amplitude of AU Mic b, we find that the accuracy of the recovered semi-amplitudes are ∼50%. So, while we quote a formal precision on the mass of AU Mic b to be ${M}_{b}={20.12}_{-1.57}^{+1.72}\ {M}_{\oplus }$ (∼9% precision), our injection and recovery tests indicate that the accuracy on the mass of AU Mic b is only known to a factor of two.

Unfortunately, attempts to recover non-injected planets are "unsuccessful," in that our modeling finds strong evidence for planets we did not inject (Figure 12) in the case of allowing P and TC to float. A deeper investigation into the posteriors of such fits indicates certain parameters (primarily P and TC) are typically not well-behaved and yield non-Gaussian distributions. When fixing P and TC to "nominal" values, our modeling does not tend to find such nonexistent planets (Figure 12).

The confident recoveries of "fake" planets in our tests indicate our GP model is flexible enough to find relatively (quasi)-stable islands in probability space with high confidence for K specifically. Although several peaks stand out in our periodogram analyses (Figures 8 and 9), more observations and/or more sophisticated modeling are needed to robustly claim these periods as statistically validated planets. We further note that the recovered values of K for the smallest injected values are inaccurate, indicating our measurement of K_c = 3.68 m s⁻¹ is also moderately unconvincing, and is likely an overestimate given the behavior of all recoveries at this level of K. We finally note this analysis is limited by planets we do not account for in the model, which may impact our ability to recover certain combinations of P and TC. Further tests using several values for the injected TC may also yield different results. With these limitations in mind, we also provide an estimation of the upper limit to the mass of AU Mic c. We find a 5σ upper limit to the semi-amplitude of AU Mic c of ≤7.68 m s⁻¹, corresponding to a mass of ≤20.13 M_⊕.

Our chromatic kernel used in this work is an initial step to exploit the expected correlation of stellar activity versus wavelength by introducing a scaling relation between wavelengths (Equations (2) and (3)). Here, we examine the "RV color" for our multiwavelength data set in order to further assess the correlation between our RVs with expected activity:

$$\begin{eqnarray}&&{\mathrm{RV}}_{\mathrm{color}}(t,\lambda ,\lambda ^{\prime} )=\mathrm{RV}(t,\lambda )-\mathrm{RV}(t,\lambda ^{\prime} ).\end{eqnarray} \tag{ 6 }$$

We first determine which nights contain nearly simultaneous measurements at unique wavelengths. We require observations to be within 0.3 day (≈6% of one rotation period) of each other, in order to minimize differences from rotationally modulated activity but increase the number of pairs for our brief use. For each nearly simultaneous chromatic pair, we compute the "data color" directly from the measured RVs as well as the "GP color" by computing the differences between the two measurements and two GPs sampled at the identical times, respectively (such that $\lambda ^{\prime} \gt \lambda$). This calculation requires knowledge of the parameters in order to remove the per-instrument zero points and realize each appropriate GP, so we make use of the MAP-derived parameters in Table 3 with kernel K _{J 2} (Equation (3)). The correlation between the data and GP color is shown in Figure 13. The agreement between the data and the model (weighted R² ≈ 0.71) indicates that our chromatic GP technique is doing a good job of reproducing the RV-color phenomenon for multiple wavelength pairs.

Zoom In Zoom Out Reset image size

With a sufficient model for stellar activity, we expect the data and GP RV color to match (up to white noise). Therefore, the "RV color" between the data and GP may be used to further constrain (in future analyses) the model (and therefore prevent overfitting) by including an effective L2 regularization penalty as follows:

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }\ {\rm{+=}}\ -{\rm{\Lambda }}\sum _{t}{r}_{\mathrm{col}}{(t)}^{2}.\end{eqnarray} \tag{ 7 }$$

Here, r _col is the vector of residuals between the GP and data RV color, Λ > 0 is a tunable hyperparameter whose value is directly correlated with the relative importance and confidence of the stellar activity model, and + = represents the standard "addition assignment" operator. The vector r _col may be computed for all pairs of wavelengths with (nearly) simultaneous measurements, and each pair can make use of identical or unique values of Λ. We finally note this regularization term is not limited to our assumption of a simple scaling relation, and could also be used in the case of disjoint kernels.

Kernels K _{J 1} (2) and K _{J 2} (Equation (3)) make use of a scaling relation for stellar activity models at different wavelengths (spectrographs) where each activity model is drawn from a Gaussian process characterized by a covariance matrix utilizing all observations. Using such joint kernels yields fits with larger scatter than cases using disjoint QP kernels (one per spectrograph; Equation (1)). In the latter case, we find that although each of the activity models appear to be "in-phase" with one another, each GP exhibits unique features that are inconsistent with a simple scaling relation (Figure 3). With nightly sampling, it is difficult to determine whether the observed differences between disjoint GPs are indicative of inadequate sampling or an inadequate RV model (activity + planets). Further, all activity models used in this work make use of identical kernel hyperparameters (excluding the amplitude), which may further be an inadequate assumption. We expect the stellar rotation period (η_P ) to be identical across wavelengths (or nearly so); however, it is not clear whether the mean activity timescale (η_τ ) or period length scale (η_ℓ ) in particular should be achromatic hyperparameters.

Our work further excluded per-spectrograph uncorrelated "jitter" terms. We suspect this may be the source of our model's ability to find planets we did not inject into the model (Section 5.3), which we defer to future work. The reduced χ² values in Tables 4 and 6 quantify the degree to which our models do not capture signals from possible additional planets, incorrect values for eccentricity and/or ω, per-spectrographsystematics not included in the formal measurement uncertainties, stellar activity such as p-mode oscillations, convection noise, or longer timescale variations. Therefore, although our specific likelihood function (Equation (4)) assumes normally distributed errors, we choose not to combine any remaining (i.e., unaccounted for by the provided error bars) potentially correlated noise into an additional uncorrelated jitter term to keep our model simple.

More accurately characterizing the masses and orbits of AU Mic b and c may require a more sophisticated stellar activity model and more intensive multiwavelength cadence. Our work further does not make use of activity indicators (e.g., Ca ii H and K, Hα) or asymmetries in the cross-correlation function (e.g., the bisector inverse slope (BIS) or differential line width dLW; Zechmeister et al. 2018) to help constrain the activity model (see Rajpaul et al. 2015). The serval pipeline in particular provides a measure of the chromaticity (CRX) for both the CARMENES-VIS and NIR data sets, which we do not use in our modeling. For AU Mic, we expect that each spectrograph is precise enough to resolve first-order chromatic effects within their respective spectral grasp's, which unfortunately will make the formal uncertainties of each spectrograph larger. Further, our QP-based kernels are primarily intended to capture rotationally modulated activity induced by temperature inhomogeneities on the stellar surface. Although the flexibility of disjoint GPs likely captures other rotationally modulated effects such as convective blueshift and limb-darkening, it will not capture short-term activity such as flares. We finally note that more seasons with high-cadence RVs will help mitigate the strong 1 yr alias present in our data set, and will help determine the correct periods for potential non-transiting planets.