Evidence that the Occurrence Rate of Hot Jupiters around Sun-like Stars Decreases with Stellar Age

The CLS catalog presented by R21 contains 179 planets detected around 719 nearby FGKM dwarfs based on the radial velocity (RV) data from the California Planet Search team (Howard et al., 2010). R21 made various data related to the CLS survey available, including spectroscopic parameters of the stars, properties of detected planets derived from RVs, and planet detection efficiencies (completeness),¹¹1https://github.com/leerosenthalj/CLSI which makes the catalog an ideal resource for occurrence rate studies. In the following subsections, we will define a subset of Sun-like stars from the CLS catalog for our occurrence analysis. The stars in our sample are shown in Figure 1 with filled orange circles. We will use the spectroscopic parameters provided in the catalog for stellar isochrone fitting in Section 2.2, and utilize the planetary parameters and detection efficiencies for inferring the planet occurrence rate in Section 4.

We fit the observed effective temperature $T_{\rm eff}$, iron metallicity [Fe/H], $K_{\rm s}$-band magnitude $K_{\rm s}$, and parallax $\varpi$ of the CLS stars with the MIST models (Paxton et al., 2011, 2013, 2015; Choi et al., 2016), and derive physical parameters of the stars including mass and age. The outcome is the samples drawn from the joint posterior distribution for the stellar parameters of each star, obtained using the jaxstar software²²2https://github.com/kemasuda/jaxstar (Masuda, 2022, 2023). The samples will be used for inferring the occurrence rate of planets in Section 3.

We use the measurements of $T_{\rm eff}$ and [Fe/H] from the CLS catalog, which were obtained in R21 by applying the Specmatch (Petigura, 2015) algorithm to the Keck-HIRES high-resolution spectra and using the Isoclassify (Morton, 2015; Huber et al., 2017; Berger et al., 2020) software package.³³3Strictly speaking, it is ideal here to use the direct outputs of Specmatch that do not rely on stellar models. We implicitly assume that this minor inconsistency does not affect the inferred age dependence given the errors we assigned for these parameters. We assume Gaussian errors of 100 K and 0.1 dex for $T_{\rm eff}$ and [Fe/H], respectively. We exclude five stars without $T_{\rm eff}$ measurements and eight stars with [Fe/H]$<-1$ in the CLS catalog. We collect the parallax measurements ($\varpi$) from Gaia DR3 (Gaia Collaboration et al., 2016, 2022) by cross-referencing using their coordinates and magnitudes. We correct the zero-points of the parallax values following the recipe of Lindegren et al. (2021) and their error bars following El-Badry et al. (2021). We obtain $K_{\rm s}$-band magnitudes from the Two Micron All Sky Survey (2MASS; Skrutskie et al., 2006) using the Gaia DR3 source IDs and correct the measurements for extinction using the dust map of Bayestar17 (Green et al., 2019), although the corrections turned out to be negligible in most cases due to their proximity and near-infrared magnitudes. We removed stars without the 2MASS IDs. Six stars in the CLS catalog do not have $K_{\rm s}$-band measurements, so we instead use $H$-band 2MASS magnitudes for these stars. In the end, all the necessary information was found for 697 stars.

We then perform isochrone fitting for the 697 stars. By fitting, we mean drawing samples from the joint posterior probability density function (PDF) of stellar parameters $\bm{\theta}$, $p_{0}(\bm{\theta}|D)$, conditioned on the data $D$ and a certain prior PDF $p_{0}(\bm{\theta})$ for “bookkeeping.”⁴⁴4The inference of the occurrence rate will not depend explicitly on the prior adopted here; see Section 3. Given the measured values of $\bm{y}^{\rm obs}=\{T_{\rm eff},\mathrm{[Fe/H]},\varpi,K_{\rm s}\}$ and their errors $\bm{\sigma}^{\rm obs}$, the likelihood function $\mathcal{L}(D|\bm{\theta})$ based on the assumption that errors follow Gaussian distributions and are independent is:

where $y_{i}(\bm{\theta})$ is computed (in a deterministic manner) by linearly interpolating the MIST grids of $(t_{\star},\mathrm{[Fe/H]},e)$ and by $\varpi=1/d$, where $t_{\star}$, $e$, and $d$ represent the stellar age, equivalent evolutionary phase (EEP; Dotter, 2016), and distance to the star, respectively. We then sample from:

We assume that the prior PDF $p_{0}(\bm{\theta})$ is separable as $p_{0}(\bm{\theta})=p_{0}(t_{\star},\mathrm{[Fe/H]},e)\,p_{0}(d)$ and set:

in order that the prior PDF is constant in the parameter space of $(t_{\star},\mathrm{[Fe/H]},M_{\star})$ where valid stellar isochrone models exist (Morton, 2015; Masuda, 2022). The priors for age, [Fe/H], and EEP are bounded within $(0.1,13.8)$ Gyr, $(-0.5,0.5)$, and $(0,600)$, respectively. For the distance prior $p_{0}(d)$, we adopt an exponentially decreasing volume density prior with a length scale of 1.35 kpc (Bailer-Jones, 2015; Astraatmadja & Bailer-Jones, 2016). We obtain 10,000 posterior samples for each of the 697 stars in the CLS catalog (black crosses in Figure 3). For all the parameters, the Gelman–Rubin convergence diagnostics $\hat{R}$ (Gelman & Rubin, 1992) computed by splitting a single chain were $\hat{R}<1.01$ in most stars and $\hat{R}<1.2$ in all the stars.

In this work, we focus on planets around Sun-like stars, for which models are more accurately calibrated than for lower- and higher-mass stars (Soderblom, 2010b). We first removed 25 poorly-fitted stars whose marginalized posterior PDFs do not contain any of the measured values within 90% credible intervals, resulting in 672 stars. We construct the subsample of Sun-like stars by choosing the stars whose mass and [Fe/H], as evaluated using the median of the posterior samples, satisfy $0.8<M^{\rm med}_{\star}/M_{\odot}<1.2$ and $-0.4<\mathrm{[Fe/H]}^{\rm med}<0.4$, respectively, and then we are left with 382 stars. Typical $68\%$ credible intervals for each posterior are approximately $0.05\,M_{\odot}$ in mass, $0.09\,{\rm dex}$ in [Fe/H], and $2.8\,{\rm Gyr}$ in age. A previous comparison with asteroseismic stars demonstrated no significant bias larger than these uncertainties for the stars in the above mass-metallicity range (Masuda, 2022).

The Sun-like sample contains 124 planets shown in Figure 2. We focus on the following two classes of planets:

• •

  Hot Jupiter (HJ): Planets with orbital period $P$ ranging from 1–10$\;{\rm days}$ and minimum mass $M\sin{i}$ between 0.3–10$\;M_{\rm Jup}$.

• •

  Cold Jupiter (CJ): Planets with orbital period $P$ between 1–10$\;{\rm years}$ (corresponding to approximately 1–5 au in semi-major axis) and minimum mass $M\sin{i}$ in the range of 0.3–10$\;M_{\rm Jup}$.

The corresponding regions of the parameter space are shown by red and cyan rectangles in Figure 2. Also shown by the background color map and white contours are the mean detection efficiencies as computed using the results of injection-recovery simulation reported in R21. The minimum masses of the planets $M\sin{i}$ and the detection efficiencies are recalculated using the host masses $M_{\star}$ derived from isochrone fitting. For most of the HJs and CJs as defined here, the uncertainties in the minimum mass and orbital period are smaller than the bin size adopted for displaying the detection efficiency.

In Figure 2, each planet is color-coded based on the age of the host star, here evaluated as the median of the marginal posterior. Albeit with large uncertainties (see also Figures 4), this figure already hints that stars hosting HJs (all having reddish colors) tend to be younger than those with CJs (mixture of blue and red).⁵⁵5One may think that the lower limit of the HJ’s mass, $M\sin{i}=0.3M_{\rm Jup}$, can impact our result because a relatively older planet-host exists at $(P,M\sin{i})=(6\;{\rm day},70M_{\oplus})$ in Figure 2. We confirmed our conclusions remained unchanged when the planets down to $M\sin{i}=0.1\,M_{\rm Jup}$ were considered as HJs. The same trend is also seen in Figure 3, which displays the distribution of the posterior medians of $(\mathrm{[Fe/H]},M_{\star},t_{\star})$ in the Sun-like sample: the stars hosting HJs (red stars) are in the lower parts of the figures at a given [Fe/H] (left panel) or mass (right panel), compared to stars with CJs (cyan circles) or stars without known planets (black crosses).

We perform a simple statistical test for the age distribution as an attempt to quantify this difference further and to highlight its issues as well. Figure 5 shows the cumulative distribution functions of ages (median of marginal posterior) for all stars in the sample (gray), stars hosting HJs (red), and stars with CJs (blue). We tested the null hypothesis that two of the subsamples are drawn from the same population using the Anderson–Darling test for $k$ samples and found the $p$-values of 0.022, 0.18, and $1.6\times 10^{-3}$ when comparing all-stars vs. HJ hosts, all-stars vs. CJ hosts, and HJ hosts vs. CJ hosts, respectively. These results confirm the above visual impression that HJ hosts are younger than CJ hosts as well as all surveyed stars.

We note, however, that interpreting this kind of statistical test is by no means straightforward. First, this does not take into account heterogeneous uncertainties in the age estimates. In general, the uncertainty varies from star to star for a variety of reasons; this may be due to the difference in the precision of the spectroscopic parameters that depends on the signal-to-noise of the observed spectrum, or due to the difference in other stellar parameters (e.g., mass) that affect the age precision. When this is the case, the above test may give low $p$-values even for two populations of stars with the same true age distributions; the distribution of the medians of the posteriors conditioned on a certain prior (i.e., the estimated ages) could still be different.⁶⁶6Consider an extreme case where there exist two populations of stars with the same age distributions and the ages can be inferred extremely precisely for one population, and only very poorly for the other. Then the distribution of the estimated ages will be very similar to the true distribution in the former but will be totally different for the latter, depending on how the estimate is given, and so the distributions of the estimated ages will look very different. Second, the apparent difference in the age distributions could be caused by the difference in other stellar parameters, such as the mass, that are correlated with age. If the occurrence rate of HJs increases with the stellar mass more rapidly than that of CJs, for example, then the HJ hosts will be younger than CJs even if their occurrences have the same age dependence at a given mass, because more massive stars are on average younger. One conceptually straightforward solution is to control the other stellar parameters in the sample by, for example, focusing on subsets of stars with similar masses and metallicities and computing the occurrence rates separately. When this is not practical — as in our case — one needs to consider the distribution of all the relevant stellar parameters in the sample and infer the dependence of the occurrence rate on these parameters simultaneously. Such an analysis must take into account the fact that the constraints on stellar parameters are often given in a degenerate manner; for example, the isochronal mass and age are strongly degenerate (Figure 4). The degeneracy in the inferred stellar parameters should not be confused with the correlation that actually exists in the stellar sample.

The rest of this paper describes our attempts to analyze the age difference taking into account the above issues associated with large and heterogeneous uncertainties in the degenerate stellar parameters. The following Section 3 presents a framework for doing so.

Two types of “occurrence rates” have been mainly discussed in the literature: (1) the number of planets per star (NPPS; e.g. Youdin, 2011; Fulton et al., 2021), and (2) the fraction of stars with planets (FSWP; e.g. Fischer & Valenti, 2005; Johnson et al., 2010). In this work, we discuss the NPPS as a function of both planet and stellar properties.

Let $\bm{z}$ and $\bm{x}$ represent the physical parameters of stars (mass, metallicity, age, etc.) and planets (mass, orbital period, radius, etc.), respectively. We define the NPPS function $f(\bm{x}|\bm{z})$ so that, around a star with given $\bm{z}$, the expected number of planets that fall within a small volume of the parameter space $d\bm{x}$ around $\bm{x}$ is (e.g., Tabachnik & Tremaine, 2002; Youdin, 2011):

Thus $f(\bm{x}|\bm{z})$ denotes the expected number of planets per star and per $\mathrm{d}\bm{x}$, but not per $d\bm{z}$.

By integrating over a certain region $\mathcal{P}$ in the $\bm{x}$-space, which defines the “planet” of consideration, we obtain the expected number of planets with $\bm{x}\in\mathcal{P}$ per star, in a manner that depends on the stellar properties:

Further integrating over $\bm{z}$, the expected NPPS around $N_{\star}$ stars whose properties $\bm{z}$ follow the probability distribution $p_{\star}(\bm{z})$ is:

where $\mathcal{S}$ denotes the domain of $\bm{z}$.

We emphasize again that $f(\bm{x}|\bm{z})$ does not have the unit of $1/\bm{z}$, and neither does $\overline{n}_{\bm{x}\in\mathcal{P}}(\bm{z})$. The PDF $p_{\star}(\bm{z})$ does, and satisfies $\int_{\mathcal{S}}p_{\star}(\bm{z})\,d\bm{z}=1$. One cannot integrate out $\bm{z}$ in $f(\bm{x}|\bm{z})$ without specifying $p_{\star}(\bm{z})$. This means the properties of all the surveyed stars need to be specified to relate $f$ with the observed number of planets. To put it the other way, the parameter distribution of the surveyed stars necessarily needs to be specified/inferred for inferring the NPPS as a function of stellar parameters.

Consider a planet survey in which the total number of stars is $N_{\star}$. We assume that the survey provides the following data sets:

• •

  $\bm{D}\equiv\{D_{j}\}^{N_{\star}}_{j=1}$: This data has information about the stellar properties $(\bm{z}_{j})$ for each star. In our case, this is the stellar data used for isochrone fitting.

• •

  $\bm{H}=\{H_{j}\}^{N_{\star}}_{j=1}$: This data set includes the parameters of planets around each star. For the $j$-th star with detected planets, the set $H_{j}$ contains the parameters for all the detected planets, represented as $H_{j}=\{\bm{x}^{(1)}_{j},\bm{x}^{(2)}_{j},...,\bm{x}^{(n_{j})}_{j}\}$, where $n_{j}$ is the number of detected planets orbiting the $j$-th star. If no planets are detected, the set will be empty. In our application, $\bm{x}$ consists of the minimum mass and orbital period of the planet.

We parameterize the stellar distribution $p_{\star}(\bm{z})$ and the NPPS function $f(\bm{x}|\bm{z})$ by the sets of parameters $\bm{\alpha}$ and $\bm{\gamma}$, respectively. We then aim to determine the joint posterior PDF of these parameters, given the data $\bm{D}$ and $\bm{H}$:

where $\mathcal{L}(\bm{\alpha},\bm{\gamma})$ is the likelihood function and $p(\bm{\alpha},\bm{\gamma})$ is the prior PDF for $\bm{\alpha}$ and $\bm{\gamma}$. We assume that the data for each star are independent when conditioned on $\bm{\alpha}$ and $\bm{\gamma}$. Then the likelihood is represented as the product of the probabilities of the data for each star:

For the $j$-th star, the likelihood may be evaluated as:

In the second line, we assume that $D_{j}$ and $H_{j}$ are independent when conditioned on $\bm{z}_{j}$, and that $D_{j}$ does not depend on $\bm{\alpha}$ nor $\bm{\gamma}$ when conditioned on $\bm{z}_{j}$. We compute $p(H_{j}|\bm{\gamma},\bm{z}_{j})$ assuming that the number of detected planets in a given small partition $\Delta_{l}$ of the parameter space follows the Poisson distribution, and that the detections in different partitions are made independently. From Equation (4), the expected number of planets in a given partition is $\eta_{j}(\bm{x}_{l})f(\bm{x}_{l}|\bm{z}_{j})\Delta_{l}$, where $\eta_{j}(\bm{x})$ is the detection efficiency of planets as a function of $\bm{x}$ for the $j$-th star. Note that $\eta_{j}(\bm{x})$ is not a density function per $\bm{x}$. Then $p(H_{j}|\bm{\gamma},\bm{z}_{j})$, the probability to find one planet in partitions around $\bm{x}_{j}^{(1)},\dots,\bm{x}_{j}^{(n_{j})}$, but none elsewhere, is (Tabachnik & Tremaine, 2002; Youdin, 2011):

Substituting this expression into Equation (3.2), we find

where we defined $\widetilde{f}_{j}(\bm{x}|\bm{z})=\eta_{j}(\bm{x})f(\bm{x}|\bm{z})$.

In this work, we model the distribution of stellar parameters in the sample, $p_{\star}(\bm{z})$, as follows:

Here, $k$ represents the index for bins that divide the domain of $\bm{z}$, and $\bm{1}_{k}(\bm{z})$ is a step function that returns 1 if the $k$-th bin contains $\bm{z}$, and 0 otherwise. The term $\exp{(\alpha_{k})}$ represents the height of the $k$-th bin, and must satisfy $\sum_{k}\exp{(\alpha_{k})}\Delta_{k}=1$ with $\Delta_{k}$ being the volume of the $k$-th bin, so that $p_{\star}$ is a normalized PDF. In this form, we assume that $p_{\star}(\bm{z}|\bm{\alpha})$ is constant within a given bin, and $\alpha_{k}$ represents the log probability density in the $k$-th bin. Using this expression for $p_{\star}(\bm{z}|\bm{\alpha})$, Equation (11) can be rewritten as:

where

Following Hogg et al. (2010), we evaluate the integral of $L_{j,k}$ via an importance sampling using the samples from the posterior PDF $p_{0}(\bm{z}|D_{j})$ conditioned on a certain bookkeeping prior $p_{0}(\bm{z})$, as we obtained in Section 2.2. Namely, considering Bayes’ theorem as

we approximate $L_{j,k}$ as:

Here, $M$ represents the number of samples drawn from the posterior PDF, and we can drop the constant factor $p_{0}(D)$ that is irrelevant to the inference. We note that the inference does not depend explicitly on the choice of the bookkeeping prior $p_{0}(\bm{z})$. This is because the evaluation of $L$ only involves $p_{0}(\bm{z}|D)/p_{0}(\bm{z})\propto p(D|\bm{z})$, which is the likelihood function that depends on the data alone, but not on $p_{0}(\bm{z})$ used to draw posterior samples in the isochrone fitting of individual stars.

Finally, the joint posterior PDF in Equation (7) is evaluated as

In this work, we perform posterior sampling using Hamiltonian Monte Carlo and No-U-Turn Sampler (Duane et al., 1987; Betancourt, 2017) as implemented in NumPyro (Phan et al., 2019; Bingham et al., 2019). The code is implemented using JAX (Bradbury et al., 2021).

Throughout this paper, the stellar distribution $p_{\star}(\bm{z})$ is modeled as a histogram, as shown in Equation (12). In this subsection, we model $\mathcal{G}(\bm{z},\bm{\gamma}_{s})$ as a histogram too, using the common bins in the $\bm{z}$ space as adopted for $p_{\star}(\bm{z})$. This is meant to serve as an analysis with minimal assumptions on the $\bm{z}$ dependence of the NPPS, which guides the interpretation of the subsequent results from other more constrained parametric models. We set up $4\times 4\times 5$ bins for [Fe/H], mass, and age spanning $[-0.4,0.4]$, $[0.8,1.2]$, and $[0,14]$, respectively, resulting in the total bin number of $K=80$. The bin widths are determined to be comparable to the typical uncertainty of these parameters (see also Section 4.3); we do not expect that the data provide useful information on finer parameter dependencies. Then we represent the expected NPPS values in each bin (assumed to be constant within the bin) by $\bm{\gamma}_{s}=\{f_{1},f_{2},...,f_{K}\}$. Consequently, we have

where $\Delta_{m}\Delta_{p}$ denotes the integral of the power law part of Equation (19) in the HJ/CJ domain, and depends both on $m$ and $p$. This factor makes $f_{k}$ to be the number of planets in the domain per star with parameters $\bm{z}$; see Equation (5).

The prior $p(\bm{\alpha},\bm{\gamma})$ is assumed to be separable as $p(\bm{\alpha},\bm{\gamma})=p(\bm{\alpha})\,p(\bm{\gamma})$. For $\bm{\alpha}$ that describes the stellar distribution $p_{\star}(\bm{z}|\bm{\alpha})$, we assume

where $\delta(x)$ denotes the Dirac delta function, and $\mathcal{U}(x;a,b)$ denotes the uniform distribution for $x$ between $a$ and $b$. The first term on the right side comes from the normalization requirement of the PDF. In the second line, we assign uniform priors on most $\alpha_{k}$, but assign essentially zero probability density $\alpha_{k}=-10$ at bins in the set $\mathcal{K}$ that satisfies

This represents our prior knowledge that stars do not exist in these regions. The value of $\alpha_{\rm max}$ is chosen so that the normalized $p_{\star}(\bm{z}|\bm{\alpha})$ is positive (i.e., $\exp(\alpha_{k})\Delta_{k}$ never exceeds unity for any $k$). For the parameters $\bm{\gamma}$ in the NPPS function $f(\bm{x}|\bm{\gamma},\bm{z})$, we assume independent and uniform priors:

We obtain 10,000 samples from the posterior distribution in Equation (16) and confirm $\hat{R}<1.01$ for all the parameters. Figure 6, Figure 7, Figure 8, and Table 1 present a summary of the result for the joint inference of the stellar parameter distribution $p_{\star}(\bm{z})$ and the NPPS function $f(\bm{x}|\bm{\gamma},\bm{z})$. Figure 6 shows the stellar parameter distribution $p_{\star}(\bm{z})$ inferred from the whole data. The left panel displays the joint distribution:

The right panels show the marginalized distributions of this joint PDF $p(\bm{z}|\bm{D},\bm{H})$ for [Fe/H], mass, and age. Figure 7 shows the inferred NPPS functions for HJs and CJs. What is shown in the right panels is essentially Equation (5) conditioned on the data, but for 1D visualization as a function of each of the mass, metallicity, and age (which we denote by $\hat{z}$), it has been marginalized over the other two stellar parameters ($\bm{z}_{\backslash\hat{z}}$). Namely, shown here by the violin plots are the distributions of

(see also Equation (5)) computed for the samples of $(\bm{\alpha},\bm{\gamma})$ drawn from $p(\bm{\alpha},\bm{\gamma}|\bm{D},\bm{H})$. From these plots, we see that the NPPS of both HJs and CJs is likely higher for more metal-rich stars. This is a known trend, but we recover this even considering the correlations with other stellar parameters. Furthermore, the NPPS of HJs shows a decreasing trend with increasing age, to a similar degree to the metallicity dependence. We further investigate these trends with the parametric models below. The left panel of Figure 7 shows the posterior PDFs for the NPPS of HJs and CJs in the Sun-like sample, averaged over the stellar parameters: the histograms show the distributions of Equation (6), $\overline{n}_{\bm{x}\in\mathcal{P},\bm{z}\sim p_{\star}(\bm{z})}(\bm{\alpha},\bm{\gamma})$, computed for the samples of $(\bm{\alpha},\bm{\gamma})$ drawn from $p(\bm{\alpha},\bm{\gamma}|D,H)$. The posterior median and equal-tail 68% intervals are $4.0^{+1.9}_{-1.2}\,\%$ for HJs (orange) and $18.0^{+5.1}_{-3.3}\,\%$ for CJs (blue), respectively. They are both consistent with the values reported by Zhu (2022), who analyzed a similar subset of Sun-like stars from the CLS catalog. In Figure 8, we show the posterior PDFs for the NPPS as a function of both mass and age, now marginalizing only over the metallicity. Although the results for separate stellar mass bins are more uncertain, the NPPS of HJs still shows a decreasing trend with increasing age at each stellar mass. This confirms that the age dependence seen in Figure 7 is not solely due to the mass dependence (cf. Section 2.4); even if we control the stellar mass, the NPPS of HJs tends to be higher at younger ages.

In this section, we model $\mathcal{G}(\bm{z},\bm{\gamma}_{s})$ as a parametric function of the form:

where $f_{0}$ is a normalization factor for the NPPS, and the function $\mathcal{G}^{\prime}(t_{\star},\bm{\gamma}^{\prime})$ that accounts for the age dependence will be specified below. The metallicity and mass dependencies are modeled as power-law functions with indices $\beta$ and $\kappa$, respectively, following the prior studies (e.g., Johnson et al., 2010). We set uniform priors for $\beta$ and $\kappa$ within $[-20,20]$ and $[-30,30]$, respectively, and a log-uniform prior for $f_{0}$ within $[10^{-5},1]$; see Table 2.

We continue to model the stellar distribution as a histogram as in Section 4.1, but given the fewer number of free parameters in $\mathcal{G}(\bm{z},\bm{\gamma}_{s})$, we increase the resolution. We set up 10 bins for [Fe/H] spanning $[-0.5,0.5]$ at 0.1 dex intervals, 5 bins for mass spanning $[0.75,1.25]\;M_{\odot}$ at $0.1\,M_{\odot}$ intervals, 7 bins for ages spanning $[0,14]$ Gyr at 2 Gyr intervals, which results in the total number of bins of $K=10\times 5\times 7=350$. The prior PDF for $\bm{\alpha}$ is again given by Equation (22).

Both the histogram analysis (Section 4.1) and the other parametric analysis adopting two different functional forms (Section 4.2) consistently suggest that the NPPS of HJs decreases with increasing age. The parametric models give tighter constraints on the NPPS than the histogram one, presumably because of the additional assumption that the NPPS changes smoothly as a function of age. The consistency of the results based on two different parametric models suggests that the inferred age dependence is not too sensitive to how this smoothness is implemented in the NPPS function.

The NPPS of HJs decreases in the latter half of the main-sequence lifetime, while the behavior in the first half is not well constrained partly due to the small number of stars searched (see Figure 9 left); the lack of young stars may be due to the selection of the Doppler survey against active stars that are not suited for precise RV measurements. Although the uncertainty in the inferred NPPS is still large due to the limited number of HJs, the result suggests that the change in the NPPS around older main-sequence stars could well be quite rapid: both of the parametric models suggest that the NPPS decreases by an order of magnitude over a few Gyrs.

In contrast, the NPPS of CJs does not show strong evidence for age dependency during the main sequence. Given that we do not know of physical processes that change their occurrence rate during the main sequence, this absence of age dependence supports that our inference framework is working properly. For evolved stars, on the other hand, the sigmoid model may hint at the decrease of NPPS; we are not fully sure if this is a real trend, given that this may partly be driven by the construction of the sigmoid model, and that the NPPS is not well constrained due to the lack of such old stars. The slight inconsistency in the inferred NPPS of CJs around the oldest stars based on the two parametric models implies that the NPPS of CJs may have a more complicated dependence on age than we assume here and that the available data may already provide such information. It is beyond the scope of the present work to perform a more detailed study of the NPPS evolution for CJs.

In this work, we adopted the MIST models for the stellar isochrone fitting, and have not examined the dependence of the results on the adopted stellar model. Tayar et al. (2022) demonstrated that the model-dependent offsets are typically $\sim 5\%$ in mass and $\sim 20\%$ in age for main-sequence and sub-giant stars. Considering the current precision of the inference, these model-dependent offsets do not significantly affect the above conclusions but may become more important in future analyses using larger samples. Because our analysis is based on the homogeneous analysis of stars using the same stellar model and the observing data with similar qualities, the conclusions based on relative comparisons, such as the NPPS difference between younger and older stars, and the difference between HJs and CJs, are more robust against the model-dependent systematics than those based on the absolute age (e.g., at what age the NPPS starts decreasing).

Tidal orbital decay appears to be the most plausible explanation for the decrease of HJ’s NPPS around older stars (Section 4). In principle, the NPPS as a function of age, as obtained in this work, enables the joint inference of the mass-period distribution of HJs before tidal evolution as well as the efficiency of tidal dissipation inside the star. Because the NPPS analysis uses the information for all the stars that have been surveyed, such an analysis can consider both HJs that currently survive and that have already been engulfed in a consistent manner, which is otherwise difficult. It is even possible, at least in principle, to include the dependence of the dissipation efficiency on the star’s mass and age, tidal forcing frequency, and planet’s mass. The information on the “initial” mass–period distribution may be useful to distinguish the origin scenarios of HJs. The empirical constraint on the tidal dissipation efficiency will be useful for improving the theory of tidal dissipation, and also for better understanding the orbital evolution of stellar binaries in general.

We leave this “full” analysis for future work, given that the NPPS function based on the CLS sample still has a large uncertainty. Instead, here we give a simple order-of-magnitude estimate for the tidal dissipation efficiency based on our inference, assuming that the decrease of the HJ’s NPPS is real and solely due to the tidal orbital decay. We characterize the efficiency of tidal dissipation by the modified stellar tidal quality factor, $Q^{\prime}_{\star}=3Q_{\star}/2k_{2}$, where $Q_{\star}$ represents the ratio of the maximum energy stored in the tide to the energy dissipated in one cycle, and $k_{2}$ denotes the tidal Love number (equals 3/2 for a homogeneous sphere). Highly dissipative stars have small $Q^{\prime}_{\star}$ values, and the $Q^{\prime}_{\star}$ value, in general, depends on tidal forcing frequency and amplitude as well as on the internal structure of the star (e.g., Barker, 2020). If $Q^{\prime}_{\star}$ is assumed to be constant, the inspiral time for a planet with orbital period $P$ is analytically given by (e.g., Ogilvie, 2014):

where $q$ is the planet-to-star mass ratio, and $\bar{\rho}_{\star}$ and $\bar{\rho}_{\odot}$ are the mean densities of the host star and the Sun, respectively. The values of $q$ and $P$ in the above equation correspond to typical HJs in our sample (see Figure 2), which should be considered as planets that have not yet been engulfed by the stars and mostly preserve their initial orbits (because the change in the orbit is small except for the last moment of tidal engulfment). The results of the parametric inferences in Figure 9 imply that they will be lost during the first half of the main sequence, roughly within $\sim 6\,\mathrm{Gyr}$. Thus Equation (5.1) implies $Q^{\prime}_{\star}\sim 10^{6}$ for a typical HJ orbiting a typical Sun-like star in our sample. We note that the corresponding inspiral time in Equation (5.1) is also comparable to the $\sim 1\,\mathrm{Gyr}$ timescale inferred from our exponential model. This estimate of the tidal quality factor is essentially similar to $\log_{10}{Q^{\prime}_{\star}}<6.5^{+0.5}_{-0.6}$ given by Hamer & Schlaufman (2019) who required $t_{\rm in}$ to be shorter than the main-sequence lifetime of Sun-like stars. The estimate is also comparable to the theoretical estimate assuming the efficient non-linear dissipation of internal gravity waves for the tidal forcing frequency of $\approx 3\,\mathrm{days}/2=1.5\,\mathrm{days}$ (see Equation 54 of Barker, 2020).⁷⁷7Barker (2020) derived $Q^{\prime}_{\star}$ scales as $P_{\rm tide}^{8/3}$ with $P_{\rm tide}$ being tidal forcing frequency. If we modify Equation 5.1 considering this dependence and assuming $P_{\rm tide}=P/2$, and apply the same argument, we estimate $Q^{\prime}_{\star}$ at $P_{\rm tide}=0.5\,\mathrm{days}$ to be $\sim 10^{5}$ instead. This also agrees with his Equation 54.

In this work, we focused on $382$ Sun-like stars from the CLS survey, which includes nine HJs. Due to this small number, the timescale for the decline of the NPPS is still consistent with a wide range of values. Our framework can be extended to transit surveys, and will also be useful to investigate HJ samples that are more than an order of magnitude larger (e.g., Yee et al., 2023). Such applications would provide more precise age dependence of the NPPS, and might also reveal the age dependence of the mass–period distribution of HJs: in this work, we ignored possible stellar parameter dependence of the shape of $f(\bm{x}|\bm{\gamma},\bm{z})$ (i.e., $m$ and $p$ in Equation 19 do not depend on stellar parameters) but this assumption may be relaxed with a larger sample.

Our result suggests that the age dependence of the HJ occurrence could be as important as that on metallicity. Then it is crucial to consider the age distribution of the sample stars as well in interpreting the results of different surveys. A general tendency would be that surveys targeting a larger fraction of older, brighter, more evolved Sun-like stars tend to give lower HJ occurrence rates, all else being equal. The age dependency might turn out to be the key to explaining marginally different survey results as mentioned in Section 1 in a unified manner. Such consideration will be even more important as the Nancy Grace Roman Space Telescope will uncover a large number of HJs orbiting stars across the Milky Way, including those in the Galactic Center (Montet et al., 2017; Miyazaki et al., 2021; Wilson et al., 2023) where we expect a large difference in the age distribution from the solar neighborhood. It is also possible for our framework to take into account the NPPS dependence on stellar properties other than mass, age, and metallicity, such as the position within the Galaxy.