An extreme test case for planet formation: a close-in Neptune orbiting an ultracool star

Funding: This work was partially supported by funding from the Center for Exoplanets and Habitable Worlds. The Center for Exoplanets and Habitable Worlds is supported by the Pennsylvania State University, the Eberly College of Science, and the Pennsylvania Space Grant Consortium. This work was supported by NASA Headquarters under the NASA Earth and Space Science Fellowship Program through grants 80NSSC18K1114. We acknowledge support from NSF grants AST-1006676, AST-1126413, AST-1310885, AST-1517592, AST-1310875, the NASA Astrobiology Institute (NAI; NNA09DA76A), and PSARC. We acknowledge support from the Heising-Simons Foundation via grant 2017-0494 and 2019-1177. G.S. acknowledges support provided by NASA through the NASA Hubble Fellowship grant HST-HF2-51519.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. Computations for this research were performed on the Pennsylvania State University’s Institute for Computational & Data Sciences (ICDS). The research was carried out, in part, at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). The Hobby-Eberly Telescope (HET) is a joint project of the University of Texas at Austin, the Pennsylvania State University, Ludwig-Maximilians-Universität München, and Georg-August-Universität Göttingen. The HET is named in honor of its principal benefactors, William P. Hobby and Robert E. Eberly. The Low Resolution Spectrograph 2 (LRS2) was developed and funded by the University of Texas at Austin McDonald Observatory and Department of Astronomy and by Pennsylvania State University. We thank the Leibniz-Institut für Astrophysik Potsdam (AIP) and the Institut für Astrophysik Göttingen (IAG) for their contributions to the construction of the integral field units.

Author Contributions: G.S. led the writing and interpretation of the result, including analyzing the RVs and host star properties. S.M. is PI of the HPF instrument, provided oversight, and assisted with analysis and interpretation. Y.M., in discussion with G.S. led the planet formation simulations. P.R. performed starspot simulations and aided in the detailed interpretation of the result. M.D. and G.S. performed the HPF RV sensitivity and occurrence analysis. C.C. performed the Gaia upper mass limit constraints. J. Winn provided expertise into the interpretation of the result and editing of the manuscript. G.Z., S.K., G.H. and G.S. reduced and interpreted the LRS2 spectra. B.Z. and E.F. contributed to the interpretation of planet formation. R.H. aided with the TESS photometric analysis. J.N., C.B., A.R., and R.T. extracted the HPF spectra and provided accurate wavelength calibration. S.D., A.Metcalf, and C.F. led the development of the HPF Laser Frequency Comb Calibrator. B.B. provided expertise into brown dwarf occurrence rates. W.C. and M.E. led observing programs contributing observing time. S.H., F.H., A.L., A.Monson, L.R., C.S., and J.Wright are part of the HPF team. All authors contributed to the interpretation of the result.

Competing interests: The authors declare no competing interests.

Data and materials availability: All time series data used (RVs, and activity indices) are available as machine-readable tables. The juliet code is available at https://github.com/nespinoza/juliet. The HPF-SpecMatch code is available at https://gummiks.github.io/hpfspecmatch/. The Panacea code is available at https://github.com/grzeimann/Panacea/.

Materials and Methods
Figs. S1 to S6
Tables S1 to S3
References (40-71)

HPF has a NIR Laser Frequency Comb (LFC) calibrator to provide a precise wavelength solution and track instrumental drifts, which has been shown to enable $\sim$20 $\mathrm{cm/s}$ RV calibration precision in 10 minute bins  (?). Following  (?), we elected not to use the simultaneous LFC calibration during the observations to minimize the risk of contaminating the science spectrum from scattered light from the LFC, and the drift correction was performed as described in  (?), which has been shown to enable precise wavelength calibration at the $\sim$$30~\mathrm{cm/s}$ level.

The HPF 1D spectra were reduced using the HPF pipeline, following the procedures in  (?),  (?), and  (?). Following the 1D spectral extraction, we reduced the HPF radial velocities using an adapted version of the SERVAL (SpEctrum Radial Velocity Analyzer) pipeline  (?), which is described in  (?). In short, SERVAL uses the template matching algorithm to derive RVs, which has been shown to be particularly effective at producing precise radial velocities for M-dwarfs. SERVAL uses the barycorrpy package  (?), which uses the methodology of  (?) to calculate accurate barycentric velocities. Following  (?), we only use the 8 HPF orders that are cleanest of tellurics, covering the wavelength regions from 8540-8890Å, and 9940-10760Å. We subtracted the estimated sky background from the stellar spectrum using the dedicated HPF sky fiber. We explicitly masked out telluric lines and sky-emission lines to minimize their impact on the RV determination.

We modeled the RVs using juliet  (?), which uses the radvel code  (?) to calculate the Keplerian model. We used the dynesty nested sampler  (?) available in juliet to perform a dynamic nested sampling of the posteriors. We modeled the Keplerian using the usual prescription of orbital period ($P$), time of inferior conjunction ($T_{\mathrm{conj}}$), eccentricity ($e$), argument of peristron ($\omega$), and RV semi-amplitude ($K$). We account for any additional noise not accounted for the errors with a white noise jitter term. To account for any potential long-period outer companions, we additionally experimented with adding a slope to the RV model. The resulting posterior on the slope parameter was fully consistent with no slope over the observing baseline, and we therefore adopt a model without a slope. Table S3 summarizes the priors and the posteriors from the RV fit.

Additionally, we obtained an $R=2,500$ optical spectrum obtained with LRS2  (?) on the Hobby-Eberly Telescope to measure the pseudo-equivalent width of the H$\alpha$ line, pEW(H$\alpha$), which we then use to estimate the fractional H$\alpha$ luminosity to the overall bolometric luminosity $\log(L_{H\alpha}/L_{bol})$. We obtained a spectrum with the LRS2-R unit in the ‘red’ channel from 6450-8400Å. Prior to measuring pEW(H$\alpha$), we shift the spectra to zero radial velocity, accounting for the barycentric velocity, and absolute velocity of the star. The LRS2 spectrum was reduced using the Panacea pipeline  (?), which performs basic CCD reduction tasks, wavelength calibration, fiber extraction, sky subtraction and flux calibration. We measure the pEW(H$\alpha$) using the following equation:

where we integrate over the limit from $\lambda_{1}=6560$Å and $\lambda_{2}=6566$Å. $F_{pc}$ is the average of the median flux in the pseudo-continuum in the ranges from $6545-6559$Å, and $6567-6580$Å after removing a linear slope fit to that range seen in the pseudo-continuum surrounding the H$\alpha$ line for late M dwarfs.

Figure S2 shows the spectrum of LHS 3154 surrounding the H$\alpha$ line compared to an LRS2 spectrum of the active M8 star VB 10  (?). VB 10 shows clear H$\alpha$ emission, while we see no H$\alpha$ emission in the LHS 3154 spectrum. For LHS 3154, we measure $\mathrm{pEW(H\alpha)}=-0.31$Å. Accounting for the S/N of the observation, we obtain a statistical uncertainty of $0.1$Å.

To estimate $\log(L_{H\alpha}/L_{bol})$, we use the following equation,

where $\log\chi$ is the ratio of the flux in the continuum near H$\alpha$ to the bolometric flux. We use the $\chi(T_{\mathrm{eff}})$ function from  (?), who estimated the functional form of $\chi$ as a function of $T_{\mathrm{eff}}$ using PHOENIX spectra. We obtain a $\log(\chi)=-4.8$, suggesting a $\log(L_{H\alpha}/L_{\mathrm{bol}})=-5.3$. This agrees with the value reported by  (?) of $\log(L_{\mathrm{H\alpha}}/L_{\mathrm{bol}})\leq-5.21$ for LHS 3154. Judging from the sample presented in  (?), these values of $\log(L_{\mathrm{H\alpha}}/L_{\mathrm{bol}})$ are suggestive of an exceptionally quiet ultracool star whose rotation period should be in excess of 100 days.

To check if TESS reveals periodicities close to the period of LHS 3154b, we used the Systematics-Insensitive Periodogram (SIP) package  (?). SIP is capable of creating periodograms while simultaneously accounting for TESS systematics within and across TESS sectors. To do so, SIP uses the TESS Target Pixel File data and apertures assigned by the TESS pipeline to build simple aperture photometry light curves of a given target while building an estimate of the background contribution. To build a periodogram, SIP simultaneously fits systematic regressors and sinusoidal components to the time-series data, including regularization terms to avoid overfitting. Figure S3 shows the resulting corrected TESS photometry along with corresponding SIP periodograms for all sectors together, along with periodograms for individual sectors. We see no periodicity at the period of the planet.

To create the SIP periodogram, we elected to use all data irrespective of quality flag, as during testing, we observed this resulted in smoother SIP-corrected light curves. No peaks were seen at the planet period when removing data that had non-zero TESS quality flags. Additionally, we also tried generating Generalized-Lomb Scargle periodograms of the PDCSAP light curves for each sector, which revealed no consistent peaks at the planet period.

Given the close-in orbit of the planet, the planet has a relatively high geometric transit probability of $R_{\star}/a=2.9\%$. Using the ephemerides from the HPF RVs, we looked for evidence of transits in the TESS data. Figure S3d shows the expected transit time using the RV ephemerides. The model in Figure S3d assumes a most likely radius of the planet of $R=3.65R_{\oplus}$ using the Forecaster package  (?). We see that the TESS data confidently rules out any transits.

Additionally, the Zwicky Transient Facility (ZTF) observed LHS 3154 in three different filters, $z_{g}$, $z_{r}$, and $z_{i}$, covering a baseline of 796, 788, and 814 days, respectively. All filters had a median exposure time of 30s. We removed 17, 6 and 4 points as high ($>1.9$) airmass points for the three filters, respectively. We additionally removed 12, 20, and 6 points as $>4\sigma$ outliers, leaving a total of 681, 695 and 152 points for periodogram analysis for the three filters, respectively. Figure S4 shows the ZTF photometry along with the associated generalized Lomb-Scargle periodograms for each individual filter. We see no peaks at the planet period in any of the three filters. In $z_{g}$, we see peaks at 28.4 days and half of that (14.3 days), which we attribute to the impact of the moon on the faint star in the bluest filter, where the impact of scattered light from the moon is expected to be the highest. In $z_{r}$ we see a peak at 38.7 days above the 0.1% false alarm level, which we attribute to be a third of the $P=116$-day peak seen in the $z_{i}$ filter. The 116-day peak seen in the $z_{i}$ filter, and the two broad peaks seen in the $z_{g}$ filter at $\sim$90-140 days, are fully consistent with the expected rotation period value of $P=114\pm 22$days as calculated by the relations from  (?) for inactive ultracool M dwarfs. Given the broad array of peaks seen in the $z_{g}$, and $z_{i}$ filters from $\sim$90-140 days, we adopt a rotation period of $114\pm 22$ days to encompass these possible peaks as being the true rotation period.

Our SOAP 2.0 models adopt the same stellar parameters listed in Table 1, except that the rotation period is assumed to be 3.7 days—matching the period of the RV signal—and the stellar inclination angle $i$ is set to $90^{\circ}$. We generated starspot signals for a single equatorial spot of varying size and spot-photosphere temperature difference ($\Delta T$).

Across a range of starspot size/temperature combinations, we found no model that reproduced the high-amplitude RV signal alongside a photometric modulation amplitude consistent with the TESS lightcurve. For a $\Delta T=600$ K—considered to be within the typical range for M dwarf spots  (?)—we find that a spot radius $R=0.15R_{*}$ is required to match the observed RV amplitude. However, assuming non-polar stellar inclinations such a spot should create photometric modulation with amplitude of $\sim$2%, which is easily ruled out by the TESS lightcurve. In combination with the collection of evidence suggesting LHS 3154 is especially old and slowly rotating, these results indicate that conventional starspot modulation cannot easily explain the combination of RV and photometric signals observed for this star.

We additionally investigated the possibility of planets in 2:1 (1.86 days) or 1:2 (7.43days) inner/our mean motion resonances (MMR) similar to the case seen in the M dwarf system GJ 876  (?). We see no significant evidence for RV signals at those periods, with joint fits of the known 3.7day signal and informative priors constrained at the 2:1 or 1:2 MMR periods returning RV semi-amplitudes that are consistent with a flat line at $2\sigma$. All signals were assumed to be circular. We place formal $3\sigma$ upper mass limits of $<2.5M_{\oplus}$ and $<5.0M_{\oplus}$ for the 1.86day and 7.43day MMR case, respectively. Given the difficulty of recognizing planetary systems in a 2:1 MMR from the present dataset, we recommend long-term RV monitoring of LHS 3154 to search for evidence of any possible orbital precession.

For each planetary system, we start the calculations assuming a disk of gas and solids and a number of initial embryos growing in the midplane. Each initial embryo has a radius of 500 km, corresponding to a mass of approximately $1.4\times 10^{-4}$ Earth masses and the total number of initial embryos in each system is a few hundred with the exact number depending on the mass of the host star and disk. For example, for the nominal case of M${}_{disk}=0.01$M_⋆, we have approximately 600 initial embryos in the disk, which adds 0.08 Earth masses (or $2.38\times 10^{-7}$ M_⊙) of solids to the system. We note that this is two orders of magnitude lower than the total mass of solids in the disk.

The embryos grow with a timescale that is determined by the velocity dispersion of the swarm of planetesimals in its surroundings. In order to model the formation of a Neptune-like planet, we included gas accretion in our calculations. Following  (?), we start the gas accretion once the planet has a mass larger than the critical mass,

where $\dot{M_{s}}$ is the accretion rate of solids. The gaseous envelope is accreted on a Kelvin–Helmholtz time-scale,

where M_p is the planet’s mass and $\tau_{KH}$ is the Kelvin–Helmholtz time-scale given by,

Planets stop the gas accretion process when gas dissipates from the disk or when the planet’s Hill radius becomes larger than two times the disk scale height, a result based on detailed numerical calculations (see  (?) and references therein). The planets do not remain in their initial orbits, and migrate with type I migration due to an imbalance in the corotation and Lindblad torques, and this switches to type II migration once the mass of the planet is big enough to allow the opening of a gap in the orbit. In our calculations, we slow down migration randomly by multiplying the rate of change of the semimajor axis by a constant value between 1 and 0.1, to take into account the uncertainties in these processes  (?, ?). Planets can additionally be trapped in resonances, and may undergo orbit crossings and possible mergers after the gas has dissipated from the disk. Each planet formation model is run for $10^{8}$ years.

The solid dust disks and gaseous disks evolve with time in our calculations. The solid disk evolves due to the accretion of the planetary embryos, which slowly decreases the amount of solids in the planets feeding zones. More importantly, the dust is also depleted globally due to the gas drag effect. Planetesimals orbit at a Keplerian speed and suffer from a drag force caused by the gas moving at sub-Keplerian speed. This causes the planetesimals to drift towards the star and deplete the disk of solids quite efficiently  (?, ?). The gas disk also evolves locally and globally: locally due to the gas accretion onto the embryos, and globally it decreases exponentially with a dissipation time-scale which we vary between $10^{6}$ and $10^{7}$ years.

The initial conditions of the protoplanetary disk strongly influence the outcomes of the planet formation simulations, in particular the final masses and architectures of the planetary systems  (?, ?). Observations of protoplanetary disks suggest a relation between the mass of the star and that of the protoplanetary disk of $M_{\mathrm{disk}}/M_{\star}=0.01$  (?). Subsequent studies showed that this relation also depends on the age of the disk  (?, ?), but for young disks of a few million years, the relation M_disk/M${}_{\star}=0.01$ is still valid considering the large dispersion in the observational data. Moreover, a recent study  (?) showed that for stars of 0.1 M_⊙, most of the observed disks have masses between 0.01 and 0.001 times that of the star. We use a relation of M_disk/M${}_{\star}=0.01$ as the nominal case, and explore extreme cases based on protoplanetary disk observations in order to explain the formation of LHS 3154b.

In total, we run four different simulation scenarios, all of which consist of 300 planetary systems that are considered stable against gravitational collapse by the Toomre Q stability criterion. The initial conditions of the four disks considered are summarized in Figure S6, and the outcomes of the planet formation simulations are summarized in Figure 3. For the first scenario, we adopt the nominal relation between the mass of the star and the mass of the disk, and we assume a $M_{\mathrm{dust}}/M_{\mathrm{gas}}=1\%$, typical of the value seen in the ISM  (?). This model is shown in Figure 3a. This model is capable of creating planets larger than Mars and up to a few Earth-mass planets, similar to the TRAPPIST-1 planets. However, it is not capable of creating planets as massive as LHS 3154b.

To create more massive planets, we explored increasing the mass of the protoplanetary disk. In Figure 3b, we increase the mass of the protoplanetary disk to the maximum dust-mass value consistent with observations of protoplanetary disks from  (?) for the lowest mass stars of $0.1M_{\odot}$ with $M_{\mathrm{dust}}=54M_{\oplus}$. Assuming a $M_{\mathrm{dust}}/M_{\mathrm{gas}}=1\%$, this corresponds to a disk mass of $M_{\mathrm{disk}}=0.14M_{\star}$. From Figure 3b, we see that this model is barely capable of creating planets that have similar masses and orbital distances as LHS 3154b.

In addition to exploring changing the protoplanetary disk mass, we explored changing the gas-to-dust ratio. A $M_{\mathrm{dust}}/M_{\mathrm{gas}}=1\%$ is typically assumed to derive total gas disk masses from measured dust disk masses. However, the mass of dust derived from protoplanetary disk observations is only sensitive to grain sizes up to $\sim$cm, and is based on the emission from regions with r$>$10 AU  (?, ?). Therefore, we take dust masses derived from observations as a minimum value in our calculations and explore different $M_{\mathrm{dust}}/M_{\mathrm{gas}}$ relations. In Figure 3c, we also show the limiting case at the upper limit of dust measurements in protoplanetary disks, but now assuming a $M_{\mathrm{dust}}/M_{\mathrm{gas}}=10\%$. This corresponds to a protoplanetary disk mass of $M_{\mathrm{disk}}=0.015M_{\star}$. From Figure 3b, we see that this model is also capable of creating planets that have similar masses and orbital distances as LHS 3154b.

We additionally note that for both the models in Figure 3b and 3c, when a massive planet is formed, it is the only massive planet formed in the system. During their formation, these planets accrete solids more efficiently and become an oligarch among the group of planet embryos, growing faster than the others and merge with the remaining embryos. This result obtained in the simulations also agrees with what we see in the LHS 3154 system: we see no evidence of other massive planets in the system.

Lastly, in Figure 3d, we keep the mass of the disk the same as in Figure 3b, while assuming the same $M_{\mathrm{dust}}/M_{\mathrm{gas}}=10\%$ as in Figure 3c. As expected, this model has the highest disk mass and the most amount of dust available ($M_{\mathrm{dust}}=500M_{\oplus}$), and is capable of creating the most massive planets. However, we note that this model would imply dust masses that are a factor of 10 larger than that those seen by protoplanetary disk observations  (?, ?). As such, we prefer the models in Figure 3b and 3c over the model in 3d, even if they marginally form planets such as LHS 3154b.

Disks are complex systems, and the disks around low mass stars are difficult to observe and are poorly understood. Along with detecting new planets around the lowest mass stars, further observations of protoplanetary disks at the lowest mass end will help improve our understanding of planet formation around the lowest mass stars.