Investigating the Origins of Hot Neptunes from Radial Velocity Data

Observations of extrasolar planets have revealed a variety of populations and enabled statistical studies to learn about their formation and evolution (e.g., Zhu & Dong 2021). Hot Neptunes are extrasolar planets that are similar in size (about 2–6 Earth radii) to Neptune in our solar system but are much closer to their stars, completing an orbit in 10 days or less (Dong et al. 2018), in comparison to 165 yr for our Neptune. The origin of hot Neptunes is not fully understood. In this work, we aim to investigate the origins of hot Neptunes based on radial velocity (RV) data.

Hot Neptunes share remarkable similarities with hot Jupiters (Dong et al. 2018). Both populations more likely reside in metal-rich stars and in systems with single-transiting planets. Furthermore, the dependence of their frequencies on host metallicity is similar, increasing by a factor of approximately 10 from subsolar to supersolar metallicity regimes. For each metallicity subsample, the frequencies of hot Neptunes and hot Jupiters are similar within a factor of 2. On the origin of hot Jupiters (see a review by Dawson & Johnson 2018), one theory suggests that the gravity of a massive third outer body brings a long-period Jupiter onto an orbit of high eccentricity (e.g., through the Lidov-Kozai mechanism; Lidov 1962; Kozai 1962 ) and the subsequent tidal dissipation leads to an orbit much closer to the star.

The first-order effect caused by the acceleration induced by the gravity of the putative third body would show up as a linear trend in the long-term RV variation of the planet system. Studying the linear trend in RV data of hot-Jupiter systems (e.g., Knutson et al. 2014) and searching for long-period companions (e.g., Bryan et al. 2016) lend support to the above dynamical origin mechanism for hot Jupiters.

Given the similarities between hot Jupiters and hot Neptunes mentioned above, it is natural to ask whether the same dynamical process can also explain the origin of hot Neptunes. The accumulated RV observations over the past decades (e.g., Trifonov et al. 2020) for hot-Neptune systems have made such an investigation possible.

In this work, we analyze the RV data for a sample of hot Neptunes, aiming to shed light on the origins of such a population of exoplanets. We present the sample of hot-Neptune systems, the RV data, and the method in Section 2. In Section 3, based on modeling the RV data, we compare the linear trend distributions of hot-Neptune and hot-Jupiter systems and constrain the distribution of the potential third bodies. We summarize and discuss our results in Section 4.

We construct a sample of hot-Neptune systems based on the NASA Exoplanet Archive. ¹

We apply various filters by defining hot Neptunes to have orbital periods of 10 days or less and sizes of 2–6 Earth radii, and by requiring that they were discovered through the transit method. In total, 580 exoplanets are identified. The period and radius distribution of these selected extrasolar planets are shown as blue points in Figure 1.

Zoom In Zoom Out Reset image size

Next, to constrain the linear trends, the selected systems need to have RV observations. The HARPS RV Bank (Trifonov et al. 2020) provides convenient RV data, with precise measurements from a uniform and improved spectroscopic pipeline. We use the HARPS RV Bank as our initial source of RV data. We search the HARPS RV Bank (Trifonov et al. 2020) and keep only those systems that have RV data in the Bank. We end up with 34 hot-Neptune systems for our investigation, which are marked as red points in Figure 1. Among them, 12 (22) appear to be single-planet (multiplanet) systems. For each system, in order to better constrain the RV trend, we further supplement the RV data with those found in literature (as listed in Table 1).

For the sample of the 34 hot-Neptune systems, their properties are similar to those of the parent sample. For example, the mean and standard deviation of the periods of the parent sample and selected sample are 6.0 ± 2.3 days and 5.9 ± 2.4 days, respectively. For radius, they are 2.9 ± 0.8 and 3.3 ± 1.0 earth radii; and for host star metallicity, 0.02 ± 0.16 and 0.07 ± 0.17.

We use exostriker (Trifonov 2019) to model the available RV data for each system in our sample and constrain the orbital parameters (e.g., period, semimajor axis, and eccentricity) of its planets and the linear trend. The initial inputs about planets' orbits are based on values found in literature. ² In particular, the number of planets and the orbit period of each planet are our adopted prior in performing the modeling. The jitter term is included in the error budget to ensure that the best fit has a reduced χ² value of unity. We find that the jitter values range from 0 to 50 m s⁻¹, in broad agreement with those found in stellar jitter studies (e.g., Bastien et al. 2014) and in hot-Jupiter studies (e.g., Knutson et al. 2014).

As an example, in Figure 2, we show the RV observation and the model fit for HD285968. The band in the left panel is in fact the sine-like curve with a period of 8.77 days from the model. When folded with the hot Neptune's orbital period of 8.77 days, the RV variation of the host star as a function of the phase can be seen in the right panel. In the left panel, the RV fit shows a clear slope. With RV data covering a baseline of ∼10 yr, the linear trend for this system is detected with high significance, (−2.57 ± 0.80) × 10⁻³ ms⁻¹d⁻¹.

Zoom In Zoom Out Reset image size

In the example shown in Figure 2, there are two sets of data (represented by blue and red points) with a clear RV offset. Between the ∼2552 day gap, an upgrade of the optical fibers for HARPS occurred in 2015 May (Trifonov et al. 2020). It may raise the question of whether the linear trend is an artifact caused by the instrument systematics in the RV measurements. However, this does not appear to be the case. First, the relative offset of the two sets of data has been accounted for in the model fitting, since it is included as a free parameter. Additionally, the model parameters (including the linear trend) remain consistent when fitting the first set of data (blue) alone. The second data set (red) has a shorter baseline, leading to weaker model constraints, but they are again consistent with those from fitting both data sets.

Figure 3 shows the RV fitting results for GJ1265, with observational data also covering a baseline of about 10 yr. This system is posed as an example of no significant detection of the linear trend, (−3.4 ± 3.9) × 10⁻⁴ ms⁻¹d⁻¹.

Zoom In Zoom Out Reset image size

HD1461 is a system with two known planets, and its RV fitting results are plotted in Figure 4. From the RV data, the effect of the two planets, with one being the hot Neptune, clearly shows up. With observation covering a baseline of ∼5000 days, the quadratic trend in the RV data is detected.

Zoom In Zoom Out Reset image size

In Table 1, the linear trend constraints for the 34 hot-Neptune systems are listed. To derive the constraints, the model only includes RV trend to the linear order. We then extend the model to have both linear and quadratic terms in the RV trend. Table 2 lists the RV trend constraints for systems with significant quadratic term.

We first analyze the derived linear trend constraints for our sample of hot-Neptune systems and compare them with those for hot-Jupiter systems. Then, we perform further analysis of the RV data to study the implications for the distributions of masses and distances of the perturbing objects and make comparisons to those for hot-Jupiter systems.

With our hot-Neptune systems, we test for correlation between the absolute value of the linear trend and planet parameters, including planet radius, orbit period, and orbit eccentricity. We find no clear dependence of the linear trend constraints on these planet parameters.

Among the 34 hot-Neptune systems, six systems have linear trend detection. In addition, four systems have quadratic trend detection. These numbers are 12 and 3 for the hot-Jupiter systems. The overall fraction of hot-Neptune systems with acceleration detection (10 out of 34) is similar to that of hot-Jupiter systems (15 out of 51).

In the top panel of Figure 5, we show the linear trend constraints as a function of planet mass for hot-Neptune systems from this work and for hot-Jupiter systems in Knutson et al. (2014). The sign of the linear trend is of no importance for our discussion and is neglected. For consistency, four hot-Neptune systems with significant quadratic trends and three hot-Jupiter systems with clear curvature in the fits are not shown. Filled circles with error bars represent systems with significant (>3σ) linear trend detection, while open circles show those without 3σ detection.

Zoom In Zoom Out Reset image size

The distributions of linear trends of hot-Neptune and hot-Jupiter systems appear to be different, if we treat the best-fit values plus and minus error bars as the whole constraints. For hot-Jupiter systems, the linear trend constraints are clustered around a few times 10⁻²ms⁻¹d⁻¹. While some hot-Neptune systems have linear trend constraints at such a level, a large faction of them have linear trends about more than one order of magnitude lower, down to a few times 10⁻⁴ms⁻¹d⁻¹.

The linear trend constraints depend on the baseline of observation. Since constraints on hot-Jupiter systems used in the comparison typically have shorter baselines than hot-Neptune systems in this work, for fair comparison, we perform a test by using shorter baselines for hot-Neptune systems. We take the average baseline (∼1900 days) for hot-Jupiter systems. For hot-Neptune systems with observation baseline shorter than 1900 days, we keep all the RV data, while for those with baseline longer than 1900 days, we only keep the observed data up to the first 1900 days. With the shorter-baseline data, we redo the model fits to obtain constraints. The number of systems with significant acceleration detection becomes seven (five with linear trend detection and two with quadratic trend detection). The linear trend constraints from the shorter baseline data are shown in the bottom panel of Figure 5. While the uncertainties in the constraints become larger, they are on average still smaller than those for hot-Jupiter systems, indicating the contribution of factors like signal-to-noise of the data to the constraints. While the overall trend in the bottom panel is similar to that in the top panel, it becomes less different from the hot-Jupiter constraints (especially when those with significant detections are compared).

To put the linear trend values into context, we can express the expected linear trend caused by a third body of mass m at distance a from the system as

$$\begin{eqnarray}&&\dot{\gamma }=5.2\times {10}^{-3}\left(\displaystyle \frac{m}{{M}_{J}}\right){\left(\displaystyle \frac{a}{10\mathrm{au}}\right)}^{-2}\cos \theta \,{\rm{m}}{{\rm{s}}}^{-1}{{\rm{d}}}^{-1},\end{eqnarray} \tag{ 1 }$$

or

$$\begin{eqnarray}&&\dot{\gamma }=5.2\times {10}^{-2}\left(\displaystyle \frac{m}{{M}_{\odot }}\right){\left(\displaystyle \frac{a}{100\mathrm{au}}\right)}^{-2}\cos \theta \,{\rm{m}}{{\rm{s}}}^{-1}{{\rm{d}}}^{-1},\end{eqnarray} \tag{ 2 }$$

where θ is the angle between the line of sight and the line connecting the star and the third body.

The linear trend constraints, if taken at face values, seem to hint that the third bodies for hot-Neptune systems on average have lower mass and/or are farther from the host stars. The six hot-Jupiter systems with their hot-Jupiter mass starting to overlap with hot-Neptune masses (≲0.3M_J ) seem to follow the linear trend constraints seen in the hot-Neptune systems. However, if we focus on those with significant detections, linear trends for hot-Neptune systems and hot-Jupiter systems do not differ substantially. In what follows, we perform more quantitative comparisons based on Kolmogorov–Smirnov tests and maximum likelihood analysis.

3.1. Linear Trend Comparisons from Kolmogorov–Smirnov Tests

3.2. Maximum Likelihood Analysis of the Third-body Distribution

To fully account for the cases of both detections and nondetections of linear RV trend, we implement a maximum likelihood method to constrain the mass and distance distributions of the potential third body and compare the results for hot-Neptune and hot-Jupiter systems.

As the first step, for each hot-Neptune system, we add a putative third body to fit the RV data (Wright et al. 2007), together with the existing planet(s). For the i-th hot-Neptune system, we find the best-fit χ² value, ${\chi }_{i}^{2}(m,a)$, at each location on a grid of semimajor axis a and third-body mass $m\sin i$ (assuming zero eccentricity; for simplicity, hereafter we use m for $m\sin i$). Similar to Knutson et al. (2014), we set the range of m to be 0.2–500M_J and that of a to be 1–100 au. The likelihood for the third body to have the corresponding m and a is computed as

$$\begin{eqnarray}&&{{ \mathcal L }}_{i}(m,a)\propto \exp \left\{-[{\chi }_{i}^{2}(m,a)-{\chi }_{i}^{2}]/2\right\},\end{eqnarray} \tag{ 3 }$$

where ${\chi }_{i}^{2}$ is the minimum of all the ${\chi }_{i}^{2}(m,a)$ values.

The results of the constraints on the mass and distance of the putative third body (a.k.a. Wright plots) for all the systems are shown in the Appendix (Figure 17). Figure 6 shows two examples. In the first example (GJ876; left panel), the combination of mass and distance is constrained. The third body can be a low-mass one closer to the star or a high-mass one farther away from the star. This is largely driven by the fact that the linear trend is proportional to m/a², as shown in Equations (1) and (2), and the RV data does not show significant trend beyond the linear order. For comparison, the RV data for the system CoRoT-24 show significant quadratic trend (see Figure 14 in Appendix A), allowing good constraints on the mass and distance of the third body (∼2M_J at ∼4 au).

Zoom In Zoom Out Reset image size

With the third body's mass and distance constraints for each hot-Neptune system, we perform the maximum likelihood analysis to obtain the constraints on the mass and distance distribution of the third body. Following previous work (e.g., Tabachnik & Tremaine 2002; Knutson et al. 2014), we model both the mass (m) and the distance (a) distribution as power law,

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}P}{d\mathrm{ln}m\,d\mathrm{ln}a}={{Cm}}^{\alpha }{a}^{\beta },\end{eqnarray} \tag{ 4 }$$

where C is the normalization factor. Instead of leaving C as a free parameter (e.g., Knutson et al. 2014), we fix it for given values of α and β by requiring $\iint {{Cm}}^{\alpha }{a}^{\beta }d\mathrm{ln}{md}\mathrm{ln}a=1$. We then introduce a fraction parameter F to explicitly incorporate the possibility that only a fraction of the systems have a third body. The third-body distribution is then described as f(m, a) = F × Cm^α a^β . By definition, F is the fraction of systems with third bodies and has values between 0 and 1. We consider three sets of hot-Neptune systems: the whole sample, the sample of single-planet systems, and the sample of multiplanet systems. For each set of samples, we aim to constrain α, β, and F. Comparisons of the constraints for different samples allow us to tell how different their third body distributions are.

Given the observational data, for each distribution model represented by the model parameters (α, β, F), the probability to have a detection for system i is computed as

$$\begin{eqnarray}&&{{ \mathcal P }}_{i}(\alpha ,\beta ,F)=\iint f(m,a| \alpha ,\beta ,F){p}_{i}(m,a)d\mathrm{ln}m\,d\mathrm{ln}a,\end{eqnarray} \tag{ 5 }$$

where p_i (m, a) is the probability in the corresponding ($\mathrm{ln}m$, $\mathrm{ln}a$) cell in the normalized Wright plot, derived based on Equation (3). The probability of nondetection is then $1-{{ \mathcal P }}_{i}$. The total likelihood of the model for a given sample is (Tabachnik & Tremaine 2002; Knutson et al. 2014)

$$\begin{eqnarray}&&{ \mathcal L }(\alpha ,\beta ,F)=\displaystyle \prod _{i=1}^{{N}_{{\rm{d}}}}{{ \mathcal P }}_{i}(\alpha ,\beta ,F)\times \displaystyle \prod _{j=1}^{{N}_{\mathrm{nd}}}\left[1-{{ \mathcal P }}_{j}(\alpha ,\beta ,F)\right],\end{eqnarray} \tag{ 6 }$$

where N_d is the number of systems with a significant (>3σ) linear trend detection and N_nd is that for nondetection.

We apply the maximum likelihood analysis to all the hot-Neptune systems. Figure 7 shows the marginalized distribution in the α–F and β–F plane. The parameter F, representing the overall fraction of systems with third bodies, is consistent with unity. In fact, the maximum likelihood occurs at F = 1, with α = − 0.253 and β = 0.114. The marginalized distribution of F peaks at F = 1, and we find F > 0.88, >0.71, and >0.52 at the 68.3%, 95.4%, and 99.7% confidence level.

Zoom In Zoom Out Reset image size

The marginalized distribution of model parameters α and β constrained from all the hot-Neptune systems is shown in the left panel of Figure 8 (blue shaded contours). We find that the contours do not vary much if F is fixed to unity. The 1σ, 2σ, and 3σ constraints for α are (−0.59, −0.17), (−0.84, 0.01), and (−1.14, 0.19), respectively. Those for β are (−0.20, 0.39), (−0.50, 0.72), and (−0.80, 1.10).

Zoom In Zoom Out Reset image size

If the hot-Neptune systems are divided into those of single- and multiplanet systems, the constraints on the model parameters appear to be consistent with each other (see the black and cyan dashed contours).

As a comparison, in the left panel of Figure 8 we also show the constraints on the third-body distribution for hot-Jupiter systems (red dashed contours; adopted from Knutson et al. 2014). The 1σ contours for hot-Neptune and hot-Jupiter systems are well separated. For parameter α, the 2σ ranges of the marginalized distributions for the hot-Neptune and hot-Jupiter systems are different, (−0.84, 0.01) versus (0.1, 2.3). The lower values of α for hot-Neptune systems suggest that the corresponding third bodies on average tend to have lower masses. For parameter β, the hot-Neptune systems have slightly lower values, indicating third bodies closer to the host stars. However, the constraints for β are mostly consistent with those from hot-Jupiter systems, e.g., with 1σ ranges of (−0.20, 0.39) versus (0.3, 1.0). As a whole, the third-body distributions for hot-Neptune and hot-Jupiter systems are consistent with each other (e.g., with 2σ contours overlapped), showing no significant differences.

In the middle panel of Figure 8, we show the comparison based on the analysis using the shorter baseline data for hot-Neptune systems. As expected, the constraints become loosened, and the difference from the hot-Jupiter results becomes even less significant.

Finally, as a test, in the right panel of Figure 8, we show the results based on an maximum likelihood analysis with weight assigned to each hot-Neptune system. The motivation is to approximately test how the constraints may change if the whole parent sample (blue points in Figure 1) is uniformly sampled. We divide the period-radius plane into 17 rectangular cells, with each cell including at least one hot-Neptune system in our sample (red points in Figure 1). We assign weight to each hot-Neptune system in our sample as the ratio of the number of blue points to that of red points in the corresponding cell and perform the maximum likelihood analysis. The resulting contours shift. For example, the contours from all hot-Neptune systems shift toward lower α and higher β. However, the constraints remain consistent with the original constraints in the left panel. They also remain consistent with the hot-Jupiter results.

Overall, while the results from the maximum likelihood analysis qualitatively agree with those indicated by the linear trend distributions and those from the KS tests, they clearly show that the difference in the third-body distributions for hot-Neptune and hot-Jupiter systems is not significant. Both the linear trend comparison and the KS tests have the degeneracy between mass and distance of the third body (e.g., Equation (1)). The maximum likelihood analysis helps break such a degeneracy, as the RV data can provide constraints on the orbital motion of the putative third body.

We analyze the acceleration trends in the RV data of hot-Neptune systems and compare their distribution to that of hot-Jupiter systems, in order to test the theory about the origin of the hot-Neptune systems under the gravitational influence of a third body.

For our hot-Neptune systems, approximately 29% (10/34) have a significant linear trend or quadratic trend, similar to the proportion in hot-Jupiter systems (15/51≃29%). For single-planet hot-Neptune systems, this number is 25% (3/12). The results strongly indicate the existence of external third bodies, supporting the dynamic origins of hot Neptunes.

In detail, the linear trend distribution simply based on the constraints for all the hot-Neptune systems and that for the hot-Jupiter systems appear to be different, with the former having more systems with lower linear trends. However, the significance of the difference cannot be readily told, as the constraints depends on observation baseline and many systems do not have linear trend detection. If we limit to systems with significant linear trend, based on KS tests, we find no evidence for difference in the linear trend distribution between hot-Neptune and hot-Jupiter systems. To fully account for acceleration trend constraints, we employ a maximum likelihood analysis to constrain the mass and distance distribution of the putative third body.

The maximum likelihood analysis of the mass and distance distribution of the external third bodies shows that the overall fraction of hot-Neptune systems with third bodies is consistent with unity, higher than 71% at a 2σ confidence level. On average, the external bodies for hot-Neptune systems are lower in mass than those for hot-Jupiter systems, but the difference is only at the 2σ level. The third bodies for hot-Neptune systems also appear to be closer to host stars than those of hot-Jupiter systems, but again the difference is not significant. In the α–β plane, the 2σ contours for hot-Neptune and hot-Jupiter systems overlap, suggesting that no significant difference is found in the distribution of third bodies.

The constraints suggest that Neptune-like planets, being less massive than Jupiter-like planets, can be affected by external bodies similar to those for hot-Jupiter systems to migrate to become hot Neptunes.

The constraints on mass and distance of the third bodies of individual hot-Neptune systems in this work can guide the search for them, e.g., through imaging observations with adaptive optics (e.g., Bryan et al. 2016). For example, nine of the hot-Neptune systems (CoRoT-24, GJ480, GJ536, HD1461, HD39091, HD77338, HD181433, HD285968, and K2-27) are interesting targets for observation, as the distance and mass of the putative third body have been tightly constrained from the RV analysis (see Figure 17). The observational search for third bodies would provide strong tests on the dynamic origins of hot Neptunes.

With RV observations of more hot-Neptune/Jupiter systems and in a longer time span, the constraints on the third bodies are expected to become tighter, which will help further test the origins of these extrasolar planet populations.

S.Z. thanks Subo Dong for his helpful guidance and insightful discussions. S.Z. also thanks the referee for their comments that helped improve the paper.

This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.

Facility: Exoplanet Archive - .

Software: exostriker (Trifonov 2019).

Model fits to the RV data and the corresponding planet phase plots for two hot-Neptune systems with significant acceleration detection are shown in the main text (Figures 2 and 4). For completeness, we show in this appendix those for the other eight hot-Neptune systems with significant acceleration detection, five with significant linear trend detection (Figures 9–13), and three with significant quadratic trend detection (Figures 14–16).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A few systems are worth mentioning. HD47186 has a significant linear trend detection (∼19σ; see Table 1). However, the model fit to its RV data (Figure 12) shows clear residuals. If we add quadratic term when fitting the RV data, the detection on the quadratic trend is barely significant (∼2.96σ). Therefore we formally put this system into the category of having a significant linear trend detection. More accurately, such residuals are a manifestation of a third body at work that will be characterized in our maximum likelihood analysis (Figure 17).

Zoom In Zoom Out Reset image size

K2-138 has as many as six planets (Figure 13). Although the RV data only span about 70 days, given the period priors we adopt for the six known planets, the fit is able to show the effects of the six planets in the folded RV data.

Zoom In Zoom Out Reset image size

For CoRoT-24 (Figure 14), although we use the period priors for the two known planets, the phase-folded RV data do not show clear trends because of the large uncertainties in the data. The quadratic trend in the RV data, however, is strong enough to be detected, leading to good constraints on the mass and distance of the putative third body (right panel of Figure 6).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For each hot-Neptune system, we fit the RV data by adding a putative body and obtain the constraints on the mass m and distance a (semimajor axis) of the third body (Wright et al. 2007). The results are shown in Figure 17, which are used in our maximum likelihood analysis. The contours correspond to Δχ² = 2.30, 6.17, and 11.80, to approximate the 1σ, 2σ, and 3σ confidence levels for two parameters.

Zoom In Zoom Out Reset image size

For CoRoT-7 and GJ9827, the fits do not converge at a large fraction of grid points. These two systems are excluded in our maximum likelihood analysis, although they are put as place holders in Figure 17.

The constraints in most systems show the expected trend of a larger mass at a larger distance. The RV data of several systems allow tighter constraints on the distribution of the third body. For the five systems with a significant quadratic trend detection, three of them (CoRoT-24, HD1461, and HD181433) have good constraints on the mass and distance of the third body. For systems with good linear trend constraints (≳2σ), some have good third-body distribution constraints, such as GJ536, HD77338, and HD285968, while GJ3634 and HD47186 show tight constraints along the m/a² degeneracy direction.