An In-depth Look at TOI-3884b: A Super-Neptune Transiting an M4Dwarf with Persistent Starspot Crossings

TOI-3884 (TIC 86263325; T_mag = 12.91; J_mag = 11.13) ²⁹ was flagged as an object-of-interest host in the TESS sector 22 (2020 February 19–2020 March 17) long-cadence (30 minutes) data by the TESS Quick Look Pipeline (QLP; Huang et al. 2020) during the Faint-Star Search (Kunimoto et al. 2022). The transit shape was noted to show an unusual shape by the TESS Follow-up Observing Program (TFOP). ³⁰ TOI-3884 was again observed by TESS in sector 46 (2021 December 4–2021 December 30) and sector 49 (2022 March 1–2022 March 25) with 2-minute exposures.

We used the lightkurve package (Lightkurve Collaboration et al. 2018) to download all three sectors assuming a "harder" quality flag, removing all NaNs and initial outliers from the Pre-search Data Conditioning Simple Aperture Photometry (PDCSAP; Figure 1; Jenkins et al. 2016; Caldwell et al. 2020). Folding the 2-minute short-cadence light curves in both sectors 46 and 49 on the expected 4.56-day period for TOI-3884b clearly shows an unusual transit shape (Figure 2)—an ingress, a bump that spans half the transit, and then the continuation of a normal transit shape through egress. This bump is also present in sector 22, though the long 30-minute cadence is too sparse to resolve any structure in transit.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We observed seven photometric transits/partial transits of TOI-3884b using three separate ground-based facilities using Bessell I, Sloan Digital Sky Survey (SDSS) ${i}^{{\prime} }$, and SDSS ${r}^{{\prime} }$ filters. We highlight each set of observations below and plot each individual transit, along with the folded TESS transits for sectors 22, 46, and 49 in Figure 2.

2.2.1. 0.3 m TMMT

2.2.2. 3.5 m ARC Telescope

2.2.3. 0.6 m RBO

We exclude potential background sources that may impact the overall transit signal (depth or shape) using the NN-EXPLORE Exoplanet Stellar Speckle Imager (NESSI; Scott et al. 2018) on the WIYN 3.5 m telescope at Kitt Peak National Observatory (KPNO) in Arizona. We took a 9-minute sequence of 40 ms exposures using NESSI's ${z}^{{\prime} }$ filter on 2022 April 18. These images were then combined and processed following the methods highlighted in Howell et al. (2011).

We plot the final contrast curve and speckle image in Figure 3. We detect no nearby sources with a Δz ${}^{{\prime} }$ magnitude brighter than 3.8'' from 0.2'' out to 0.8'' and magnitudes brighter than 5 from 0.8 out to 1.2''. We compliment this with archival Gaia DR3 data (Gaia Collaboration et al. 2021), which find no nearby sources within 20''. Gaia also assigns TOI-3884 a renormalized unit weight error (RUWE) equal to 1.25, which is consistent with a single star (Belokurov et al. 2020; Ziegler et al. 2020). TOI-3884 is a single star in a fairly sparse region of the night sky.

Zoom In Zoom Out Reset image size

We performed an intensive RV follow-up campaign of TOI-3884 using the Habitable-zone Planet Finder (HPF; Mahadevan et al. 2012, 2014) starting on 2021 December 1. HPF is a high-resolution (R ∼ 55,000), near-infrared (810–1280 nm), fiber-fed (Kanodia et al. 2018), stabilized (Stefansson et al. 2016) precision RV spectrograph on the 10 m Hobby-Eberly Telescope in Texas (Ramsey et al. 1998). Over the next 5 months, we observed TOI-3884 on 27 nights, with each night obtaining two 945 s exposure measurements. Each spectrum was analyzed using the HxRGproc package, which corrects for bias, nonlinearities, and cosmic rays and then calculates the flux and variance of the individual spectra as described in Ninan et al. (2018). We use barycorrpy (Kanodia & Wright 2018) to perform the barycentric correction on the individual spectra, which is the Python implementation of the algorithms from Wright & Eastman (2014). A wavelength solution was created by interpolating the wavelength over all other exposures in the same night of each observation, which was then applied to the respective TOI-3884 spectra.

We removed all nights (eight total) that possessed unbinned S/Ns less than 50% of the expected S/N of 74 at 1.04 μm calculated from the HPF Exposure Time Calculator. ³¹ These S/Ns ranged between 21 and 31. An inspection of these low-S/N observations determined that they were all obtained during less-than-optimal sky conditions (variable seeing >2'', background ${i}^{{\prime} }$-band magnitude was brighter than 16.5, transparency was <75%, and/or bad weather or clouds were noted in the night logs). Every other observation possessed an S/N > 43 and met our required observing conditions for transparency, seeing, and good weather conditions. We also removed the spectra from 2022 April 5, as these were observed during the transit spanning the large bump. As the planet is crossing an active region of the star, this may introduce potential contamination in the RV signal. This left 36 unbinned spectra taken over the course of 18 nights (Figures 4 and 5).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We applied a template-matching method (Anglada-Escudé & Butler 2012) using the SERVAL pipeline (Zechmeister et al. 2019) modified for HPF (Stefánsson et al. 2020). A master template was created by combining all spectra and masking tellurics and sky emission lines. This template was then shifted to match each individual spectrum by minimizing the χ² statistic before converting this shift into velocity space. We binned the two nightly individual RVs reported from SERVAL using a weighted average based on their respective S/Ns. The final binned RVs used for our analysis are listed in Table 1 and are plotted in Figure 5.

We derived the spectroscopic stellar parameters, effective temperature (T_eff), metallicity ([Fe/H]), and log g, of TOI-3884 by applying the template-matching methodology on the HPF spectra as outlined in Stefánsson et al. (2020). Using the HPF-SpecMatch package (Stefánsson et al. 2020), we apply the spectral-matching technique to the HPF Order 5 spectra (853–864 nm), which has little to no telluric contamination. We list the spectroscopically derived stellar parameters for TOI-3884 from this analysis in Table 2.

Table 2. Summary of Stellar Parameters for TOI-3884

Parameter	Description	Value	Reference
TOI	TESS object of interest	3884	TESS mission
TIC	TESS Input Catalogue	86263325	Stassun
2MASS	⋯	J12061746+1230249	2MASS
Gaia DR3	⋯	3919169687804622336	Gaia DR3
Equatorial coordinates, proper motion, and visual extinction:
αJ2016	Right Ascension (R.A.)	181.571808	Gaia DR3
δJ2016	Declination (decl.)	+12.507030	Gaia DR3
d	Distance in pc	43.1	Gaia DR3
${A}_{V,\max }$	Maximum visual extinction	0.04	Green
Optical and near-infrared magnitudes:
B	Johnson B mag	17.46 ± 0.23	APASS
V	Johnson V mag	15.74 ± 0.01	APASS
${g}^{{\prime} }$	Sloan ${g}^{{\prime} }$ mag	16.62 ± 0.09	Pan-STARRS1
${r}^{{\prime} }$	Sloan ${r}^{{\prime} }$ mag	15.17 ± 0.06	Pan-STARRS1
${i}^{{\prime} }$	Sloan ${i}^{{\prime} }$ mag	13.58 ± 0.06	Pan-STARRS1
J	J mag	11.13 ± 0.02	2MASS
H	H mag	10.55 ± 0.02	2MASS
Ks	Ks mag	10.24 ± 0.02	2MASS
W1	WISE1 mag	10.16 ± 0.02	WISE
W2	WISE2 mag	9.99 ± 0.02	WISE
W3	WISE3 mag	9.76 ± 0.05	WISE
SpecMatch spectroscopic parameters:
Teff	Effective temperature in K	3180 ± 88	This work
[Fe/H]	Metallicity in dex	0.04 ± 0.12	This work
$\mathrm{log}(g)$	Surface gravity in cgs units	4.97 ± 0.05	This work
Model-dependent stellar SED and isochrone fit parameters:
Ms	Mass in M⊙	0.298 ± 0.018	This work
Rs	Radius in R⊙	0.302 ± 0.012	This work
ρs	Density in g cm−3	15.26 ± 2.04	This work
Other stellar parameters:
$v\sin {i}_{s}$	Rotational velocity in km s−1	3.59 ± 0.92	This work
Prot(is = 90°)	Nontilted maximum rotational period in days	4.22 ± 1.09	This work
ΔRV	"Absolute" radial velocity in km s−1	3.16 ± 2.89	Gaia DR3

Note.

References are as follows: Stassun (Stassun et al. 2018), 2MASS (Cutri et al. 2003), Gaia DR3 (Gaia Collaboration et al. 2021), Green (Green et al. 2019), APASS (Henden et al. 2018), WISE (Wright et al. 2010), Pan-STARRS1 (Chambers et al. 2016; Magnier et al. 2020).

Download table as:  ASCIITypeset image

With the isochrones package (Morton 2015), we create a spectral energy distribution (SED) fit using the combination of the derived stellar spectroscopic values; the g, r, i, z, and y magnitudes reported from Pan-STARRS1 (PS1; Chambers et al. 2016; Magnier et al. 2020); the W1, W2, and W3 Wide-field Infrared Survey Explorer (WISE) band magnitudes (Wright et al. 2010); the J, H, and K magnitudes reported by the Two Micron All Sky Survey (2MASS; Cutri et al. 2003); and the parallax from Gaia DR3 (Gaia Collaboration et al. 2021). We utilize Gaussian priors for all parameters except for a flat prior on the A_V extinction and flat-log age prior up to 2 Gyr (see Section 3.1). We utilize the relations in Green et al. (2019), calculated for the Gaia-reported distance of 43 pc, to place an upper A_V extinction limit of 0.1. We determine a stellar mass and radius of 0.298 ± 0.018 M_⊙ and 0.302 ± 0.012 R_⊕, respectively (stellar density: 15.26 ± 2.04 g cm⁻³), for TOI-3884. We verify these values by repeating the same fits using the ExoFASTv2 package (Eastman et al. 2019), deriving masses and radii within 1σ. We verify this stellar density using the high-precision 2022 April 5 APO transit, where we obtain a best-fit density of 15.43 ± 0.39 g cm⁻³ (assuming a circular orbit).

Using ESPRESSO, Almenara et al. (2022) suggest that TOI-3884 is a slowly rotating star with a slow $v\sin i$ of 1.1 km s⁻¹. They note that this slow rotation suggests an inactive star—in contrast with the large spot-crossing event. We use our HPF spectra in an attempt to verify the slow rotator scenario by constraining the rotational broadening of TOI-3884 using two separate methods. First, during the spectral-matching process, HPF-SpecMatch performs an optimization for the optimal rotational broadening (see Stefansson et al. 2020, for discussion). This results in a $v\sin i=3.6\pm 0.9$ km s⁻¹. Second, we compare the widths of cross-correlation functions (CCFs) of TOI-3884 to the CCF widths of the artificially broadened slowly rotating reference star of a similar spectral type. The HPF-SpecMatch analysis highlights Ross 128 as an excellent spectral match to TOI-3884b with a T_eff = 3192 ± 60 K (Mann et al. 2015), which matches well with the effective temperature of TOI-3884 of T_eff = 3180 ± 80 K (Table 2). Further, Bonfils et al. (2018) demonstrate that Ross 128 is an inactive slowly rotating M dwarf with a long rotation period of >100 days, suggesting minimal rotational broadening.

Figure 6 compares the CCFs of TOI-3884 to the CCFs of Ross 128 from six HPF orders clean of tellurics, suggesting that a $v\sin i\gt 3$ km s⁻¹ is warranted, and we derive a $v\sin i=3.2\pm 0.9$ km s⁻¹ estimate from the average and the standard deviation values from the six HPF orders, respectively. We perform the same check by using other slowly rotating standard stars for comparison and retrieve comprable  $v\sin i$ values for TOI-3884. We elect to formally adapt the $v\sin i$ value derived from HPF-SpecMatch, as through its χ² minimization process of the full spectra it can better account for differences in normalization offsets that could lead to differences in the CCFs. We were unsuccessful in resolving the slower 1.1 km s⁻¹ $v\sin i$ as originally published for this star.

Zoom In Zoom Out Reset image size

The relatively rapid $v\sin i$ from this work suggests that TOI-3884 should be active and relatively young (<1 Gyr; Newton et al. 2016). We support this conclusion with the LAMOST spectra, which cover Hα. LAMOST reports an Hα equivalent width (EW) of −3.86 ± 0.02 Å in emission. From Equation (1) in Newton et al. (2017), an inactive star with TOI-3884's properties should have an Hα EW of 0.18 Å in absorption. TESS observes three large flaring events in the two short-cadence sectors, and the HPF spectra also show clear Ca IR triplet (Ca IRT) excess in emission.

We apply pyHammer (Roulston et al. 2020) to the archival Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST; Yuan et al. 2015; Xiang et al. 2017) spectra assuming the metallicity derived from the HPF spectra. Using a template-matching routine of empirical M-dwarf spectra, we determine the best-fit spectral type to be either an M4 or M5 dwarf. We adopt M4 as the spectral type for this work.

Significant spot coverage can affect the measured photospheric temperature of the star and influence the SED-derived stellar mass/radius. However, we detect no significant deviation from a single-star SED fit to the observed magnitudes using ExoFASTv2. Moreover, we calculate that a spot covering fraction of 15% with a spot temperature of 2900 K will impact the actual stellar temperature by 50–100 K for TOI-3884. This is within the T_eff uncertainty reported by HPF-SpecMatch. Therefore, the derived stellar parameters in Table 2 are minimally affected by the large spot (see Section 4).

We perform a joint analysis of the transit and RV observations to measure TOI-3884b's mass and radius. However, the lack of a pristine non-spot-crossing transit light curve of TOI-3884b creates a challenge for transit analysis. We create an individualized starspot mask based on visual inspection of each ground-based transit and of the three folded TESS sectors. We then fit a transit model to the unmasked points. We calculate a ${\chi }_{r}^{2}$ of the residuals along with a by-eye examination. Points that still demonstrate bump structure are masked, and we repeat the procedure until we have minimized the ${\chi }_{r}^{2}$. Masked duration and location vary slightly between data sets, though all fell between 39 minutes prior to mid-transit (T₀) and 25 minutes after T₀, i.e., ∼67% of the transit duration. Figure 2 plots the ground-based and folded TESS light curves, with the masked points denoted in blue.

We perform a joint fit with exoplanet (Foreman-Mackey et al. 2021) using both the masked TESS and ground-based transits and the HPF RVs. We fit for a single transit ephemeris, period, impact parameter, a/R_s , and transit depth using the combination of the three masked TESS sectors, three TMMT transits, and three APO observations. We include a dilution term to the transit depth for each of the three TESS sectors fixed on the ARCTIC SDSS ${i}^{{\prime} }$ transit. exoplanet uses the built-in starry (Luger et al. 2019) package to model the quadratic limb-darkening parameters. As each instrument uses a different broadband filter, we fit for quadratic limb-darkening terms specific to the various wavelength coverage (TESS bandpass, Bessell I, SDSS ${i}^{{\prime} }$, and SDSS ${r}^{{\prime} }$). Lastly, we included a jitter term added in quadrature to the flux errors and a flux offset to each transit observation. Neither the TESS out-of-transit baseline (PDCSAP) nor the ground-based observations required additional detrending. We plot our best-fit transit model in red for each transit in Figure 2.

We assume a Keplerian model for the RVs allowing eccentricity and the argument of periastron (ω) to float, as well as the RV semiamplitude. Similar to the photometric transits, we include a jitter and RV offset terms, as well as a general trend line. Including a Gaussian Process (GP) had no effect on the RV results; thus, we do not include an activity-dependent GP for the RV orbit. We determine that TOI-3884b is a super-Neptune with a mass of ${32.59}_{-7.38}^{+7.31}\,{M}_{\oplus }$ and radius of 6.43 ± 0.20 R_⊕. Figure 5 plots the best-fit RV model, along with the 1σ contours. We report the best-fit properties from this joint analysis and the final planetary properties for TOI-3884b in Table 3. We note that the derived mass of TOI-3884b is based on the model's assumption that the planet is the main source of the RV variation. Periodograms of the Ca IRT, differential line widths, and chromatic index show no peaks with false-alarm probabilities <10% at the planet's period (nor any other period), indicating that the RV signal is not dominated by stellar activity.

Table 3. Derived Parameters for TOI-3884b

Parameter	Units	Value a
Orbital parameters:
Orbital period	P (days)	4.5445828 ± 0.0000098
Eccentricity	e	${0.06}_{-0.04}^{+0.06}$
Argument of periastron	ω (rad)	−1.96${}_{-0.04}^{+4.28}$
Semiamplitude velocity	K (m s−1)	28.03${}_{-6.23}^{+6.06}$
RV trend ...	dv/dt (m s−1 yr−1)	${0.58}_{-4.92}^{+4.78}$
RV jitter ...	σHPF (m s−1)	${7.86}_{-5.11}^{+5.68}$
Transit parameters:
Transit midpoint	TC (BJDTDB)	2459556.51669±0.00025
Scaled radius	Rp /Rs	0.197 ± 0.002
Scaled semimajor axis	a/Rs	${25.90}_{-0.71}^{+0.96}$
Orbital inclination	i (deg)	${89.81}_{-0.18}^{+0.13}$
Impact parameter	b	0.089${}_{-0.061}^{+0.082}$
Transit duration	T14 (days)	${0.0666}_{-0.0024}^{+0.0019}$
Dilution b	DTESS S22	0.98 ± 0.12
	DTESS S46	0.86 ± 0.03
	DTESS S49	0.84 ± 0.03
Planetary parameters:
Mass	Mp (M⊕)	${32.59}_{-7.38}^{+7.31}$
Radius	Rp (R⊕)	6.43 ± 0.20
Density	ρp (g cm−3)	${0.67}_{-0.16}^{+0.18}$
Semimajor axis	a (au)	0.0361 ± 0.0008
Planetary insolation	S (S⊕)	6.29 ± 0.84
Equilibrium temperature c	Teq (K)	441 ± 15

Notes.

^a The reported values refer to the 16%–50%–84% percentile of the posteriors. ^b Dilution due to presence of background stars in TESS aperture, not accounted for. ^c We assume the planet to be a blackbody with zero albedo and perfect energy redistribution to estimate the equilibrium temperature.

Download table as:  ASCIITypeset image

We apply the program STarSPot (STSP; Morris et al. 2017; Schutte et al. 2022) to the high-precision 2022 April 05 APO observation with the procedure outlined in Schutte et al. (2022). STSP is specifically designed to model the light curves of transiting systems in which starspots and/or faculae create localized surface brightness variations on the host star. Using an affine-invariant Markov Chain Monte Carlo (MCMC; Morris et al. 2017) optimizer, a single run of STSP samples different radii (R_spot/R_s ), latitudes, and longitudes (θ and ϕ, respectively) for every spot but applies a fixed spot contrast (defined by its temperature and the filter of the observed transit) in order to break the known degeneracies between these properties.

For TOI-3884, we assume a photospheric temperature of 3200 K derived from the HPF spectroscopic parameters and perform STSP runs with the following set of spot temperatures: 2600, 2700, 2800, 2900, 3000, and 3100 K. We calculate the contrasts defined by these temperatures by first interpolating PHOENIX synthetic (Husser et al. 2013) spectra at both the stellar and spot temperatures. We then integrate over the specific filter response curve of the observed light curve and sum the integrated flux in that filter. Finally, we calculate the contrast for each spot temperature by dividing the integrated spot flux by the integrated photospheric flux (M. Schutte et al. 2023, in preparation).

We first consider the star–planet orientation in which the spin axis of the star is in the plane of the sky and aligned with the orbital axis of the planet. In this case, the measured $v\sin i$ provides a rotation period of P_rot = 4.22 ± 1.09 days, which is consistent with the orbital period (P_orb = 4.54 days). Adopting a rotation period for the star of P_rot = 4.54 days that is synchronous with the orbital period and applying the STSP program to the APO ${i}^{{\prime} }$-band transit, we quickly find that a single circular starspot is insufficient to describe the structure of the feature, but a three-spot model, with one large spot surrounded by two smaller spots, produces the lowest χ², regardless of spot contrast. Even after forcing the model to have two medium-sized spots in the place of the large central spot, the optimization preferred one large pole spot combined with the two smaller nearby spots.

The best-fit three-spot model is shown in Figure 7, consisting of one large central spot (R_spot/R_s = 0.44) surrounded by two smaller spots (R_spot/R_s = 0.10 and 0.07, respectively). This model has a reduced χ² of 2.14 and corresponds to a spot temperature of 2900 K (contrast = 0.5). In comparison, the best-fit one- and two-spot models have reduced χ² values of 2.25 and 2.16, respectively. If we compare the AIC_c (Sugiura 1978) values between the one-, two-, and three-spot models (−9.24, −6.12, and −2.86, respectively), the AIC_c favors the one-spot model, as it has the least amount of fitting parameters. It is important to note that both the one- and two-spot models do not fit the data points as well by eye (the one-spot model is shown in Figure 7). Therefore, even though the AIC_c favors the simplest one-spot model, the χ² statistic prefers the three-spot model. Additionally, the spot feature is asymmetric, which is hard to fit with fixed circles as is required by STSP, but if we instead have two spots close together, this can replicate an asymmetric feature, which further adds to the three-spot model being the best-fit model. Spot temperatures of 2700 and 2800 K produce equally good solutions (reduced χ² values of 2.19 and 2.17, respectively, which fall within ${\sigma }_{{\chi }^{2}}=34$) with similar spot configurations. While the reduced χ² values are larger than 1, the data are ground-based data with likely underestimated error bars. Hotter spot temperatures (3000 and 3100 K) cannot reproduce the amplitude of the feature and are therefore not possible solutions regardless of the spot configuration. The coolest spot temperature of 2600 K produces a similar spot configuration, but the reduced χ² falls just outside of the above variance. Given the discrete 100 K sampling of our models, we find the temperature of the large spot to be between 2700 and 2900 K, with a preference for the hotter 2900 K spot. From these fits we constrain the radius of the large spot to be R_spot/R_s = 0.44 ± 0.08.

Zoom In Zoom Out Reset image size

This scenario produces a very large (radius of ∼44% the star's radius) starspot that will produce significant photometric out-of-transit variability of >1% if that is the only large feature on the star. Interestingly, we do not detect any clear photometric modulations in the two short-cadence TESS sectors or in the publicly available ground-based monitoring with the Zwicky Transient Facility (ZTF; Masci et al. 2019), the All-Sky Automated Survey for Supernovae (ASAS-SN; Kochanek et al. 2017), and the Asteroid Terrestrial-impact Last Alert System (ATLAS; Tonry et al. 2018). Figure 8 shows TESS sector 46 data as an example, with sector 49 showing the same pattern. To explore this scenario fully, we model the full light curve for both TESS sectors 46 and 49 with the in-transit starspots fixed but with an additional three spots allowed to vary. We found that in order to decrease the out-of-transit variability enough to be less than the noise level (<0.5%) of the TESS light curves, there must be additional spots such that as the star rotates there is always a near-equal fraction of spotted area (20%) rotating out of view as is rotating into view. Thus, it is possible that the photometric spot modulation of TOI-3884 is simply hidden within the noise assuming that the spots are configured such that they are uniformly spread across the surface of the star and cover a large fraction of the star.

Zoom In Zoom Out Reset image size

This leads to the concern that the RV-observed 4.56-day signal is partially due to the stellar rotation and not the planet. However, none of the HPF activity indicators show any periodic signal, which would be characteristic for large spots (Robertson et al. 2020).

While observations cannot formally exclude this scenario, the requirements are contrived: TOI-3884 must have a rotation exactly equal to its planet, possess a spotted surface such that the photometric variability is <0.5% over the two TESS sectors, and maintain the same starspot with very little evolution across the transit chord while keeping the rest of the transit chord nearly spot-free. We therefore disfavor this hypothesis.

Tilting the star such that the spot does not rotate in and out of view would minimize the out-of-transit variability (Jackson & Jeffries 2012). In this scenario, TOI-3884's large spot must be located on or near the pole of the star, with the star's spin axis inclined (i_s ) such that the pole of the star is pointed toward the observer. In order for the spot feature to occur at the same phase in the first half of the transit, the spin axis of the star and the planet's orbital axis must be misaligned (i.e., λ ≠ 0). Because STSP assumes that the planet's position and the star's tilt are well known, it is not designed to derive the optimal stellar inclination and λ. However, the fixed phase of the spot feature allows us to constrain the tilt of the star's spin axis and λ. For example, if the pole of the star is pointed exactly toward the observer (stellar inclination of 0°), the bump would occur in exactly the middle of the transit regardless of λ. Conversely, a tilt that is too close to the plane of the sky (i.e., inclination >60°) would produce spot crossings during ingress of the transit. Thus, we first performed a comprehensive search of every stellar inclination value from 60° down to exactly pole-on in increments of 5°. For our search, we fixed the rotation period to be exactly equal to the orbital period, and we assumed that one spot with a radius of 30% the radius of the star was directly on the southern pole of the star. From our search, we found the only possible i_s was 25° ± 5°. After determining the best stellar inclination for the star, we then performed a series of simulations that varied λ from 0° to 180° in increments of 10°. From our search, we determined that λ = 75° ± 10° provided the best fit to the APO SDSS ${i}^{{\prime} }$ light curve. It is important to note that the provided uncertainties were derived from an exploration of the possible stellar inclinations that fit the data, and then, the uncertainties for λ assume that the stellar inclination is constant. Since these parameters are actually entwined, a more formal determination of the error bars is left to future work.

Once we determined the stellar spin axis and λ for the misaligned scenario, we modeled the starspots in the same way as before with STSP where the radii and locations of the spots are allowed to vary. For this scenario, we assume the same number of spots and spot temperature as found for the best-fit aligned model (three spots with temperatures of 2900 K), as it is likely that the spot temperature and number of spots are the same no matter the tilt of the star. We found the best-fit stellar surface features for this scenario to be one large spot (R_spot/R_s = 0.29) that is slightly off-center to the pole, with two smaller spots on either side (R_spot/R_s = 0.16 and 0.09, respectively). This spot configuration is shown in Figure 9 and has a final reduced χ² of 2.6. We also fit this scenario with one large spot instead of three spots, though the reduced χ² of this model was 3.96. Similar to the results with the nontilted star, the AIC_c value for the one-spot model is lower than the three-spot model (−10.38 and −3.04, respectively) owing to the difference in parameters. The one-spot model does not fit the data well by eye, and again the χ² statistic favors the three-spot model. Therefore, we opt for the three-spot model.

Zoom In Zoom Out Reset image size

Finally, we calculate the out-of-transit variability for the TESS short-cadence data and find that this spot configuration produces variability well below the TESS noise level with no additional large spots needed (see Figure 8). Table 4 reports the best-fit values from this analysis.

We use SOAPv2 (Dumusque et al. 2014) to test the impact of a pole spot on the RVs. A pole spot under this scenario will inject a ∼10 m s⁻¹ signal, which is 30% of the RV semiamplitude. However, SOAPv2 assumes an optical bandpass, as it was designed specifically for the HARPS wavelength range (380–700 nm). Stellar activity decreases at longer wavelengths, where the contrast between the spot and photosphere temperatures decreases (Reiners et al. 2010). As HPF operates at near-infrared wavelengths, we expect the overall impact of the spot to be suppressed by ∼2 × (Reiners et al. 2010; Robertson et al. 2020). Thus, the pole spot's impact with this configuration is <5 m s⁻¹—within the 1σ semiamplitude uncertainty of 28.0 ± 6.3 m s⁻¹.

We extend our spot model derived from the SDSS ${i}^{{\prime} }$ transit to the APO SDSS ${r}^{{\prime} }$ and TESS short-cadence observations. We calculate new contrast values for the TESS and SDSS ${r}^{{\prime} }$-band filters for a spot temperature of 2900 K. Using the exact spot-complex configuration for the misaligned system, we model the ${r}^{{\prime} }$ band and TESS short-cadence transits using STSP. The results showed that the same spot configuration could not fit either the TESS or APO ${r}^{{\prime} }$-band transits.

A close inspection of individual TESS transits suggests slight changes in spot amplitude, duration, and location (though it always starts during ingress). However, the lack of precision within the individual TESS transits makes it nearly impossible to map spot evolution of subsequent transits.

For the APO SDSS ${r}^{{\prime} }$ transit, we chose to assume the same number of spots and spot temperature while allowing the location and spot radii to vary. We discover that the polar spot remains approximately the same radius but shifts slightly, while the other two spots slightly increase in area (blue spots in Figure 10). Thus, while the general locations of the features near the pole are consistent across 6 months, the individual starspots are most likely evolving from one transit to the next. This tentative evidence for small-scale spot changes suggests caution when directly comparing transits observed at different times. If we instead allow the spot temperature to change rather than the spot location and radii, we find that the required contrast to fit the APO SDSS ${r}^{{\prime} }$ transit is nearly perfectly dark (c = 0.90) if we assume the same spot configuration as found in Figure 9. Since the spot contrast required to fit the SDSS ${r}^{{\prime} }$ transit is unreasonably dark, it is more likely that there is small-scale spot evolution between transits.

Zoom In Zoom Out Reset image size

While our qualitative conclusions generally agree with those from Almenara et al. (2022), there are notable exceptions.

Stellar rotation period.—We determine a stellar $v\sin i$ of 3.59 ± 0.92 km s⁻¹ for TOI-3884, in contrast to Almenara et al. (2022), who find a $v\sin i$ of 1.1 km s⁻¹. We arrive at the more rapid rotational value via both the HPF-SpecMatch and CCF methods. We also attempted to duplicate their method by using a template star nearly identical to their comparison star, LHS 1140, ³² which is slightly cooler and has slightly lower mass than TOI-3884. However, we still obtain a $v\sin i$ of 3.5 km s⁻¹. Using the LAMOST-derived Hα EW, models from Newton et al. (2017) constrain TOI-3884's rotation period <10 days—faster than a $v\sin i$ of 1.1 km s⁻¹.

Almenara et al. (2022) do note significant variation in their CCF profiles between their two spectra, likely a result of the large spot. In turn, we do not observe significant variation in either the differential line width (dLW) or the chromatic index (CRX) across the 34 unbinned spectra. We attribute this to HPF's near-infrared wavelength range, as well as HPF's lower resolution of R ∼ 53,000 (compared to ESPRESSO's bluer wavelength coverage and higher resolution). Thus, it is also possible that the discrepancy between reported $v\sin i$ is due to stellar activity—which HPF is less affected by.

Planetary mass.—Almenara et al. (2022) derive a planetary mass of ${16.5}_{-1.8}^{+3.5}$ M_⊕ using two RV ESPRESSO points. This is a 2.2σ discrepancy from our higher mass measured with HPF. We checked for correlations between the HPF RVs and the dLW and CRX. We found no statistically significant correlation, with a ρ of −0.191 (p-value: 0.282) and 0.171 (p-value: 0.334) for the respective activity indicators. It is therefore unlikely that activity alone is amplifying the periodic planetary signal observed by HPF. We attempted to jointly fit both the HPF and ESPRESSO RV points; however, our model requires a 14 m s⁻¹ jitter term added to the ESPRESSO RVs. ESPRESSO's bandpass of 380–780 nm is more susceptible to stellar activity (Reiners et al. 2010). SOAPv2 approximates that our pole spot model should introduce a ∼10 m s⁻¹ signal into optical ESPRESSO RVs—similar to the required jitter of our model. With 17 near-IR HPF RVs, we robustly measure a 4σ planetary mass.

Planetary radius and transit depth chromaticity.—We measure a larger planet of 6.43 ± 0.20 R_⊕ compared to Almenara et al. (2022), who report two different radii: 6.31 ± 0.28 R_⊕ from GP fitting and 6.00 ± 0.18 R_⊕ from starry. We test the impact of our transit mask on the best-fit transit depth and planetary radius derived from the joint fit. To accomplish this, we fit a transit model to the 2022 April 5 APO observation without any mask deriving a planetary radius of 6.05 ± 0.14 R_⊕—comparable to Almenara et al. (2022) but a poor fit to the APO light curve. We then created 10 random spot masks of various positions and sizes, though we required the mask to include the most extreme spot-crossing event (between 30 minutes pre-transit and 7 minutes post-transit). A transit model was again fit to these light curves, where we derived planetary radii spanning 6.39–6.61 R_⊕, well within 1σ of our reported planetary radius, as well as in good agreement with the Almenara et al. (2022) GP-derived radius for TOI-3884b. Therefore, the discrepancy between the two reported radii is not dependent on our mask selection alone.

When further investigating this discrepancy, we discovered that their transits (see Figure 1 in Almenara et al. 2022) demonstrate significant chromatic variability even outside of the modeled pole spot. Chromatic transits can be explained via either a background eclipsing binary (Wang et al. 2021) or unocculted stellar activity (Rackham et al. 2018). As we rule out nearby companions, we explore the impact of stellar activity on our transit depth.

We fit a transit model to the two masked APO SDSS ${i}^{{\prime} }$ and SDSS ${r}^{{\prime} }$ transits holding a/R_s , impact parameter, transit ephemeris, and transit depth (R_p /R_s )² constant across both transits. We check for chromaticity in the masked APO transit depths between the SDSS ${i}^{{\prime} }$ and SDSS ${r}^{{\prime} }$ observations, which could indicate a contaminating background source. Using the two masked transits, we fit a transit using exoplanet (Foreman-Mackey et al. 2021) holding a/R_s , impact parameter, transit ephemeris, constant and including a dilution term multiplied to the SDSS ${r}^{{\prime} }$ transit depth. Assuming no chromaticity, the dilution term is equal to 1. We determine a dilution term of 1.01 ± 0.02. There is no significant transit depth difference between these two wavelength bandpasses. However, both these depths are slightly shallower than the bluer g' transit in Almenara et al. (2022) and deeper than the IR ExTRA observation. Differing spot contrasts compared to the hotter photosphere create deeper transits at bluer wavelengths (Rackham et al. 2018). We approximate by eye the depths of their individual transits. We fit both our and their wavelength-dependent depths using a simple unocculted starspot model from Rackham et al. (2018):

$$\begin{eqnarray}&&{D}_{\mathrm{obs}}=\displaystyle \frac{D}{1-{f}_{\mathrm{sp}}(1-\tfrac{{F}_{\lambda ,\mathrm{sp}}}{{F}_{\lambda },\mathrm{ph}})},\end{eqnarray} \tag{ 1 }$$

where D_obs is the wavelength-dependent observed transit depth, D is the true transit depth, f_sp is the spot coverage fraction, and F_λ is the wavelength-dependent flux of the spot (sp) and photosphere (ph), respectively. The unocculted spot model fits the four transit depths assuming a spot temperature of 2900 K, total unocculted spot coverage of 16%, and a true planet radius of ∼6.2 R_⊕. We do not include uncertainties in these numbers, as the Almenara et al. (2022) transit depths and uncertainties rely on by-eye approximations. However, this model demonstrates that the chromatic transit depth is explained by unocculted stellar activity that slightly impacts the measured radius of the planet (∼1σ discrepancy from our radius).

Stellar inclination and spot properties.—We fit a spot model based on the values reported in Almenara et al. (2022) to the SDSS ${i}^{{\prime} }$ transit, determining a best-fit ${\chi }_{r}^{2}$ of 6.09. Assuming that the spot is evolving over time, it is possible that the spot evolved between the two observations. However, their spot model generates significant (∼1%) out-of-transit variability that is not observed in the TESS or ground-based photometry. Accounting for the lack of baseline variability enabled us to constrain the stellar inclination to 25° ± 5°, which in turn impacted our overall best-fit spot model.

Super-Neptunes (4 R_⊕ < R_p < 8 R_⊕) represent a transitional population of planets between the rocky terrestrial planets and Jovian gas giants. TOI-3884b adds to the growing sample of well-characterized (with precise 3σ masses and radii) super-Neptunes orbiting M dwarfs. In particular, Figure 11(a) shows TOI-3884b's position in a planetary mass–radius plane with respect to other M-dwarf (T_eff < 4000 K) planets with known (>3σ) masses and radii. Figure 11(b) plots the same sample as a function of stellar effective temperature.

Zoom In Zoom Out Reset image size

The formation of Jovian planets around M dwarfs should be inhibited by the longer orbital timescales with respect to the disk lifetimes (Laughlin et al. 2004; Ida & Lin 2005). This is corroborated by empirical data from RV surveys (Endl et al. 2006; Johnson et al. 2010; Maldonado et al. 2019; Sabotta et al. 2021). However, Laughlin et al. (2004) predict that M dwarfs should host an abundance of Neptunes that fail to accrete a massive-enough core (∼10 M_⊕; Pollack et al. 1996) in a timely manner to initiate runaway gaseous accretion. Using models from Fortney et al. (2007) and propagating the uncertainties in planetary parameters using the Monte Carlo method, we predict TOI-3884b's core mass to be about 21 ± 4 M_⊕. It should therefore have experienced some runaway gaseous accretion. The fact that TOI-3884b did not accrete a Jovian-mass atmosphere suggests that its core was slow to form, or there was a lack of nearby gas/material for rapid accretion, or both.

While the spot portion of the transit may complicate the analysis, TOI-3884b has the highest transmission spectroscopic metric (TSM; Kempton et al. 2018) of any known planet with an equilibrium temperature <500 K (TSM: 230; Figure 12). Due to its bright host star and large transit depth, this planet also has one of the highest metrics of any known non–hot Jupiter planet. At 430 K, TOI-3884b's atmosphere likely contains methane as the carbon-dominant molecule, along with water and some ammonia—assuming equilibrium chemistry (Zahnle & Marley 2014; Kempton et al. 2017; Fortney et al. 2020; Hu 2021). Deriving the overall abundances of these molecules would provide an approximate C/N/O ratio, a useful measurement for constraining where this planet originally formed in its disk (Öberg & Wordsworth 2019; Turrini et al. 2021; Hobbs et al. 2022). Öberg et al. (2011) demonstrate the connection between C/O ratios and various molecular snowlines. Dash et al. (2022) expand on this study by noting that nitrogen provides information surrounding the disk's overall metallicity, as it is unaffected by the condensation of molecules such as water, carbon dioxide, and methane. Only the ammonia snowline at ∼100 K and N₂ snowline at ∼78 K significantly affect its overall ratio in the disk. TOI-3884b is therefore an extremely promising target to observationally test the link between nitrogen abundance and formation location.

Zoom In Zoom Out Reset image size

Due to the combination of its cool equilibrium temperature along with experiencing UV radiation from its active M-dwarf host, TOI-3884b's atmosphere is likely composed of photochemically created hazes such as tholins (e.g., Morley et al. 2013) or even soot (Gao et al. 2017). Photochemically created hazes (e.g., Tsai et al. 2022), and aerosols in general, are common in exoplanetary atmospheres, with several studies linking their presence to temperature (Crossfield & Kreidberg 2017; Dymont et al. 2021; Yu et al. 2021). Assuming that the trend highlighted in Yu et al. (2021) holds, we would expect TOI-3884b to possess a fairly hazy transmission spectrum in the near-infrared. However, Kawashima et al. (2019) show that hazes should become translucent at longer wavelengths assuming that the overall particle sizes are small. Thus, extending out to >3 μm should enable atmospheric characterization of TOI-3884b regardless of its expected hazy atmosphere.

We generate the expected transmission spectrum of TOI-3884b using ExoTransmit (Kempton et al. 2017) assuming a 100× solar metallicity atmosphere with no aerosols, aerosols at pressures of 100 μbar, and aerosols at pressures of 10 μbar (Figure 13). For these simulations we assume a gray-opacity aerosol layer that is wavelength independent. Using the cloud-free model, we used PandExo (Batalha et al. 2017) to simulate two transit observations with JWST/NIRSpec-Prism. We find that we should easily retrieve methane and water in both the cloud-free and 100 μbar cases. At 10 μbar we will still observe methane absorption features with tentative detection of water. It should be noted, however, that these simulations assume a typical transit shape that allows for easily derived uncontaminated transit depths. Assuming that the bump feature in TOI-3884b's transit remains long-lived, it is possible that the uncertainties presented in Figure 13 are underestimated, and stellar contamination will also need to be included in modeling the observed transmission spectrum.

Zoom In Zoom Out Reset image size

TOI-3884b presents a second unique opportunity: the impact of starspots on the transmission spectrum of a planet. Stellar activity due to an inhomogeneous photosphere may introduce spurious features into a planet's transmission spectrum (e.g., Rackham et al. 2018; Barclay et al. 2021). This is of particular concern with M dwarfs, whose spots are cool enough to host their own water absorption features (Jones et al. 1995). Rackham et al. (2023) emphasize the need for future in-depth studies into untangling star and planetary spectra, especially as we begin to probe the atmospheres of terrestrial worlds with JWST and future instruments. With half of its transit covered by a spot and the other over the photosphere, comparing the resulting transmission spectrum from the first half of the transit to that from the second may yield an unprecedented probe into the effects of cooler spots.

Assuming the pole spot hypothesis, TOI-3884b possesses a misaligned orbit with an obliquity of 75° ± 10°. TOI-3884b therefore joins the growing population of misaligned warm Neptunes (R > 4 R_⊕), which includes the two M-dwarf Neptunes: GJ 3470b (Stefànsson et al. 2022) and GJ 436b (Bourrier et al. 2018).

Neither Gaia nor the HPF RV residuals detect evidence of any outer massive companion in the TOI-3884, which could have been responsible for TOI-3884b's misaligned orbit (Petrovich et al. 2020). Of the four misaligned Neptunes orbiting K and M dwarfs, two (HAT-P-11b and WASP-107b; Yee et al. 2018; Piaulet et al. 2021) have a confirmed outer companion, while Stefànsson et al. (2022) do not exclude the existence of an outer planet in the GJ 3470 system. Giants around M dwarfs are uncommon; it is unlikely that TOI-3884 hosts an additional gas giant responsible for the misalignment of TOI-3884b. However, our RV observations are limited to <6 months. Continued RV monitoring is required to detect longer-period massive planets in this system.