Hot Jupiter formation in dense stellar clusters: A Monte Carlo model applied to 47 Tucanae

Abstract

We study the efficiency of high-e migration as a pathway for Hot Jupiter formation in the dense globular cluster 47 Tuc. Gravitational N-body simulations are performed to investigate the orbital evolution of star-planet systems due to dynamical stellar perturbations. Planetary systems that have been scattered into orbits of sufficiently high eccentricity can undergo tidal circularisation, with Hot Jupiter formation being one possible stopping condition. We also account for the possibility of (i) ionisation due to high-energy encounters, (ii) tidal disruption of the planet by tidal forces inside the Roche limit and (iii) Warm Jupiter formation. The orbital evolution of a population of cold Jupiter progenitors, with initial semi-major axes between 1-30 au, is simulated over 12 Gyr using a simplified dynamical model of 47 Tuc. Our computational treatment of dynamical encounters yields an overall HJ occurrence rate of F_HJ = 5.9 x 10^-4 per cluster star (a 51 per cent enhancement relative to the analytic baseline). The probability of Hot Jupiter formation is highest in the core and falls off steeply beyond a few parsecs from the centre of the cluster, where the stellar density is too low to drive efficient eccentricity diffusion. The code can be found here: https://github.com/James-Wirth/HotJupiter.

Figures (12)

• Figure 1: The mean error $\Delta(\xi_{\mathrm{crit}})$ due to our truncation of the perturbing orbit at $\xi = \xi_{\mathrm{crit}}$ (Equation \ref{['equation:xi']}) is shown alongside the mean improvement $t_{\mathrm{int}}/t_{\mathrm{bm}}$ in the integration time. The mean is taken over a randomised sample of $500$ encounters with $\Pi \sim$47Tuc.
• Figure 2: The eccentricity excitation $\epsilon$ and semi-major axis excitation $\alpha$ are plotted against time for $v_{\infty} = 6 \ \mathrm{km} \ \mathrm{s}^{-1}$ and $v_{\infty} = 24 \ \mathrm{km} \ \mathrm{s}^{-1}$, with the remaining encounter parameters set to $b = 15 \ \mathrm{au}$, $\Omega = i = \omega = 1 \ \mathrm{rad}$, $e_0 = 0.3$, $a_0 = 1 \ \mathrm{au}$, $m_{\star} = m_{\mathrm{pert}} = M_{\odot}$ and $m_{\mathrm{p}} = M_{\mathrm{J}}$. The simulated results are averaged over the initial planetary phase. The experiments show good convergence to the benchmark value ($\xi_{\mathrm{bm}} = 10^{-10}$, blue) well within the restricted integration time. The encounter with $v = 24 \ \mathrm{km} \ \mathrm{s}^{-1}$ is secular, causing a permanent change in eccentricity but no permanent change in semi-major axis.
• Figure 3: The mean binned relative error $\Delta := |(\epsilon_{\mathrm{sim}} - \epsilon_{\mathrm{hr}})/\epsilon_{\mathrm{hr}}|$ of the analytic approximation compared to the simulated values, obtained from a sparse sample of $10^5$ strictly non-ionising encounters with $\Pi \sim$47Tuc, as a function of $(\mathcal{T}, \mathcal{S})$. We display the results for velocity dispersion $\sigma = 6 \ \mathrm{km} \ \mathrm{s}^{-1}$. The mean relative error is of order $\Delta \lesssim 10\%$ in the region $\mathcal{D}$ where $\mathcal{T} \gtrapprox \mathcal{T}_{\mathrm{min}} = 15$ and $\mathcal{S} \gtrapprox \mathcal{S}_{\mathrm{min}} = 300$.
• Figure 4: The simulated eccentricity excitations ($\epsilon_{\mathrm{sim}}$) are compared against the analytic results ($\epsilon_{\mathrm{sim}}$) over a range of encounter speeds $v_{\infty}/ \mathrm{au} \ \mathrm{yr}^{-1} \in [10^{-1}, 10^{3}]$ and an arbitrary selection of orientations. The remaining free parameters were set to $b = 50 \ \mathrm{au}$, $e_0 = 0.3$, $a_0 = 1 \ \mathrm{au}$, $m_{\star} = m_{\mathrm{pert}} = M_{\odot}$ and $m_{\mathrm{p}} = M_{\mathrm{J}}$. The relative error (lower panel) is small in the region $\mathcal{D}$ of parameter space where $T/T_{\mathrm{min}}>1$ and $S/S_{\mathrm{min}}>1$
• Figure 5: The eccentricity distribution of a statistical ensemble of planetary systems after $3$ Gyr of dynamical evolution due to stellar encounters. The initial eccentricity distribution is a sharply peaked Gaussian with $e_0 \sim \mathrm{Normal}(\mu=0.3, \sigma=10^{-3})$ and the remaining free parameters are set to $a_0 = 1 \ \mathrm{au}$, $m_{\star} = m_{\mathrm{pert}} = M_{\odot}$, $m_{\mathrm{p}} = M_{\mathrm{J}}$. The analytic MC curves (dashed histograms) rely solely on the Heggie1996 approximation, whereas the hybrid MC curves (solid histograms) employ direct N-body simulations in the domain $\mathcal{D}^{\mathrm{C}}$ of parameter space where the approximation breaks down (Section \ref{['section:MonteCarlo_Simulation']}). The solid black curve is the solution of the PDE in Equation \ref{['equation:eccentricity_diffusion']}. The MC simulations are computed for $b_{\mathrm{max}}/\mathrm{au} \in \{50,75,100\}$ to demonstrate the good convergence for $b_{\mathrm{max}} \gtrapprox 50 \ \mathrm{au}$. The hybrid MC results exhibit an enhanced rate of eccentricity diffusion toward high eccentricities compared to the analytic prediction.