Exploring the formation by core accretion and the luminosity evolution of directly imaged planets: The case of HIP 65426 b

TL;DR

This paper investigates the formation and luminosity evolution of the directly imaged planet HIP 65426 b, using population synthesis, formation modeling, and scattering simulations to constrain its properties and possible formation pathways.

Contribution

It combines observational constraints with population synthesis and dynamical models to explore formation scenarios of wide-orbit giant planets like HIP 65426 b.

Findings

HIP 65426 b's mass is estimated around 10 Jupiter masses.

A formation pathway involving core accretion followed by outward scattering is plausible.

The planet is predicted to have a high eccentricity and possibly companions at smaller orbits.

Abstract

A low-mass companion to the two-solar mass star HIP65426 has recently been detected by SPHERE at around 100 au from its host. Explaining the presence of super-Jovian planets at large separations, as revealed by direct imaging, is currently an open question. We want to derive statistical constraints on the mass and initial entropy of HIP65426b and to explore possible formation pathways of directly imaged objects within the core-accretion paradigm, focusing on HIP65426b. Constraints on the planet's mass and post-formation entropy are derived from its age and luminosity combined with cooling models. For the first time, the results of population synthesis are also used to inform the results. Then, a formation model that includes N-body dynamics with several embryos per disc is used to study possible formation histories and the properties of possible additional companions. Finally, the outcomes of two- and three-planet scattering in the post-disc phase are analysed, taking tides into account. The mass of HIP65426b is found to be Mp = 9.9 +1.1 -1.8 MJ using the hot population and Mp = 10.9 +1.4 -2.0 MJ with the cold-nominal population. Core formation at small separations from the star followed by outward scattering and runaway accretion at a few hundred AU succeeds in reproducing the mass and separation of HIP65426b. Alternatively, systems having two or more giant planets close enough to be on an unstable orbit at disc dispersal are likely to end up with one planet on a wide HIP65426b-like orbit with a relatively high eccentricity (>~ 0.5). If this scattering scenario explains its formation, HIP65426b is predicted to have a high eccentricity and to be accompanied by one or several roughly Jovian-mass planets at smaller semi-major axes, which also could have a high eccentricity. This could be tested by further direct-imaging as well as radial-velocity observations.

Figures

Tables

Table 1. Table 1: Stellar and disc parameters used for the N 𝑁 N -body simulations.

Parameter	Value
Stellar mass	2 $M_{\odot }$
Stellar radius	2 $R_{\odot }$
Stellar temperature	10,000 K
Disc inner boundary	0.1 au
Disc outer boundary	200 au
Initial surface density exponent	$-1.5$
Initial surface density ${\Sigma }_0=\Sigma \left(1\mathrm{au}\right)$	8655 g cm-2
Disc metallicity	1 $\times$ solar
Photoevaporation factor ${\Phi }_{41}$	1000
Background viscous $\alpha$	$5\times {10}^{-3}$

Equations

$$L_{\mathrm{pf}}^{\text{hottest}}$$

$$L_{\mathrm{pf}}^{\text{cold-nom.}}$$

$$L_{\mathrm{pf}}^{\text{cold-class.}}$$

$$L_{\mathrm{pf}}^{\text{coldest}}$$

$$\frac{{\mathrm{d}}^{2}N}{\mathrm{d}{m}_{\mathrm{p}}\mathrm{d}{s}_{\text{pf}}}=p(m_{\mathrm{p}},s_{\text{pf}})=p_{s_{\text{pf}}}(m_{\mathrm{p}},s_{\text{pf}})\times p_{m_{\mathrm{p}}}(m_{\mathrm{p}}).$$

$$\displaystyle p_{s_{\textrm{pf}}}({m_{\rm p}},s_{\textrm{pf}})=\left\{\begin{array}[]{l}1\;\;\;\;\mbox{if $s_{\textrm{pf, min}}({m_{\rm p}})<s_{\textrm{pf}}<s_{\textrm{pf, max}}({m_{\rm p}})$}\\ 0\;\;\;\;\mbox{otherwise}.\end{array}\right.$$

$$s_{\text{pf, min}}^{\text{cold}}$$

$$s_{\text{pf, max}}^{\text{cold}}$$

$$s_{\text{pf, min}}^{\text{hot}}$$

$$s_{\text{pf, max}}^{\text{hot}}$$

$$p_{m_{\mathrm{p}}}(m_{\mathrm{p}})=\frac{\mathrm{d}N}{\mathrm{d}{m}_{\mathrm{p}}}\propto {m_{\mathrm{p}}}^{-1},$$

$$\frac{a_1-a_2}{a_1}<1.46{\left(\frac{m_1+m_2}{m_{\star }}\right)}^{2/7},$$

$$a_{\mathrm{r}}=\frac{a_1{a}_2(m_1+m_2)}{m_1{a}_2+m_2{a}_1}.$$

$$a_{\mathrm{r}}$$

$$e_{\mathrm{r}}$$

$$\Delta L_{\mathrm{tide}}$$

$$\Delta E_{\mathrm{tide}}$$

$$\Delta L_{\mathrm{tide}}$$

$$\Delta E_{\mathrm{tide}}$$

$${\tau }_a$$

$${\tau }_e$$

$$\displaystyle\textrm{AMD}=\sqrt{Gm_{\star}}\Bigg{(}$$

$$\displaystyle+m_{2}\sqrt{a_{2}}\left(1-\sqrt{1-{e_{2}}^{2}}\right)\Bigg{)}.$$

$$\frac{a_{\mathrm{int}}}{a_{\mathrm{out}}}<1-1.36{\left(\frac{m_1+m_2}{m_{\star }}\right)}^{1/5}{c}_{\min }^{1/10},$$

$$c_{\min }<2\left(2-\sqrt{1-e_1^2}-\sqrt{1-e_2^2}\right).$$

Full text

*Aims. * We want to derive statistical constraints on the mass and initial entropy of HIP 65426 b and to explore possible formation pathways of directly imaged objects within the core-accretion paradigm, focusing on HIP 65426 b.

*Methods. * Constraints on the planet’s mass and post-formation entropy are derived from its age and luminosity combined with cooling models. For the first time, the results of population synthesis are also used to inform the results. Then a formation model that includes $N$-body dynamics with several embryos per disc is used to study possible formation histories and the properties of possible additional companions. Finally, the outcomes of two- and three-planet scattering in the post-disc phase are analysed, taking tides into account for small-pericentre orbits.

*Results. * The mass of HIP 65426 b is found to be ${m_{\rm p}}=9.9^{+1.1}_{-1.8}\leavevmode\nobreak\ {M_{\textnormal{J}}}$ using the hot population and ${m_{\rm p}}=10.9^{+1.4}_{-2.0}\leavevmode\nobreak\ {M_{\textnormal{J}}}$ with the cold-nominal population. We find that core formation at small separations from the star followed by outward scattering and runaway accretion at a few hundred astronomical units succeeds in reproducing the mass and separation of HIP 65426 b. Alternatively, systems having two or more giant planets close enough to be on an unstable orbit at disc dispersal are likely to end up with one planet on a wide HIP 65426 b-like orbit with a relatively high eccentricity ($\gtrsim 0.5$).

*Conclusions. * If this scattering scenario explains its formation, HIP 65426 b is predicted to have a high eccentricity and to be accompanied by one or several roughly Jovian-mass planets at smaller semi-major axes, which also could have a high eccentricity. This could be tested by further direct-imaging as well as radial-velocity observations.

1 Introduction

In the last decade, direct imaging efforts have revealed a population of super-Jovian planets at large separations from their host stars. It has been well established that these planets are rare; only a small percentage of stars possess such a companion (Bowler, 2016). What is not yet clear is whether the formation process is intrinsically inefficient there and how important post-formation architectural changes to the system (through migration or interactions between protoplanets) are. The formation mechanism that produces these planets has not yet been convincingly identified. The main contenders are the different flavours of core accretion (CA; with planetesimals or pebbles building up the core) and of gravitational instability (GI; with or without tidal stripping). Therefore, given the current low numbers of detections, every new data point can represent an important new challenge for planet formation.

The first discovery of the SPHERE instrument at the VLT (Beuzit et al., 2008, 2019), HIP 65426 b, is an ${m_{\rm p}}=8$–12 ${M_{\textnormal{J}}}$ dusty L$6\pm 1$ companion to the $m_{\star}=1.96\pm 0.04\leavevmode\nobreak\ M_{\odot}$ fast rotator HIP 65426, which has an equatorial velocity $v_{\star}\sin i=299\pm 9$ km s*-1*. Its projected separation is $92.0\pm 0.2$ au, and the star is seen close to pole-on (Chauvin et al., 2017). If the planet is not captured and its orbital plane is the same as the midplane of the star, the projected separation is very close to the true separation.

In this paper, we set out to explore how core accretion could lead to the objects observed in direct imaging. We take a closer look at HIP 65426 b because it is of low mass and is at a relatively large separation, while its host star is a fast rotator. Essentially, we are following up on the comment in Chauvin et al. (2017) that the ‘planet location would not favor a formation by core accretion unless HIP 65426 b formed significantly closer to the star followed by a planet–planet scattering event.’

This paper is structured as follows. In Section 2 we use planet evolution models and work backwards from the observations to derive joint constraints on the mass and initial (i.e. post-formation) entropy of HIP 65426 b, where ‘initial’ refers to the beginning of the cooling. We then switch to a forward approach and study the possible formation of HIP 65426 b. In Section 3 we use detailed planet formation models following the disc evolution and $N$-body interactions. Then in Section 4 we use $N$-body integrations to look in detail at interactions between several companions once the disc has cleared. Finally, in Section 5 we present our conclusions and a discussion.

2 Constraints on the mass and post-formation entropy

In this section we use the luminosity to derive, with planet evolution models, constraints on the mass and initial (post-formation) entropy of HIP 65426 b, following the approach of Marleau & Cumming (2014), as also applied to $\kappa$ And b and $\beta$ Pic b (Bonnefoy et al., 2014a, b). The idea is to explore the parameter space of mass and initial entropy through the Markov chain Monte Carlo (MCMC) method. We assume Gaussian error bars on the logarithm of the luminosity and on the linear age. Both flat and non-flat priors in linear mass and post-formation entropy are considered, as detailed in Section 2.4.

2.1 Luminosity and age

Firstly, we discuss the input quantities for the MCMC. The adopted bolometric luminosity is $\log L/L_{\odot}=-4.06\pm 0.10$ as derived by Chauvin et al. (2017). Contrary to estimates based on the photometry in individual bands, this quantity should be robust as it is based on the comparison to young L5–L7 dwarfs with a similar near-infrared spectrum.

The typical age of stars in the Lower Centaurus–Crux group around HIP 65426 is $14\pm 2$ Myr, but the placement of phase-space neighbours of HIP 65426 in a Hertzsprung–Russell diagram suggests an age of 9–10 Myr. This lead Chauvin et al. (2017) to adopt an age of $14\pm 4$ Myr. We note that at 2 $M_{\odot}$, HIP 65426 is predicted by stellar evolution models to have a pre-main sequence lifetime of approximately 15 Myr (see the overview as a function of stellar mass in fig. 1 of Dotter, 2016), so that it is approaching the main sequence or has only recently joined it. A point to consider is that if HIP 65426 b formed by core accretion (CA), its cooling age would be smaller by a not entirely negligible formation delay ${\Delta t_{\rm form}}$ (Fortney et al., 2005; Bonnefoy et al., 2014b), which we now briefly discuss.

While the dependence of the formation time ${\Delta t_{\rm form}}$ on stellar mass has not yet been studied in detail, it seems plausible that giants form more quickly around more massive stars. Since runaway accretion proceeds very quickly by construction, it is the oligarchic growth phase that dominates the total formation time. For instance, Thommes et al. (2003, their eq. (11)) found that in this regime, the planet growth rate scales as $\dot{M}\propto m_{\star}^{1/6}\Sigma_{m}\rho_{\mathrm{gas}}^{2/5}$, where $\Sigma_{m}$ and $\rho_{\mathrm{gas}}$ are respectively the surface density of planetesimals and the (midplane) gas density. This scaling reflects in part the fact that the core accretion rate is proportional to the Keplerian frequency, which at fixed orbital distance increases with stellar mass. Since both $\Sigma_{m}$ and $\rho_{\mathrm{gas}}$ are expected to increase with stellar mass, the formation time should decrease with planet mass. Also, observationally, the formation time ${\Delta t_{\rm form}}$ is unlikely much longer than 3 Myr since discs around more massive stars are shorter-lived (Kennedy & Kenyon, 2009; Ribas et al., 2015); already at solar masses, gas giants must form typically in at most ${\Delta t_{\rm form}}\sim 3$–5 Myr given the lifetimes of protoplanetary discs (Haisch et al., 2001). Finally, population synthesis calculations for a 2 $M_{\odot}$ central star (Mordasini et al., in prep.) indicate that most $\sim 10\leavevmode\nobreak\ {M_{\textnormal{J}}}$ planets (approximately the mass of  HIP 65426 b, as we show later) have reached their final mass after roughly ${\Delta t_{\rm form}}\sim 2$ Myr, and the simulations presented in Section 3 using the Coleman & Nelson (2016a) models for a 2 $M_{\odot}$ star yield ${\Delta t_{\rm form}}\approx 2.5$–4 Myr. Therefore, we adopt ${\Delta t_{\rm form}}=2$ Myr, and thus $t_{\rm cool}=12\pm 4$ Myr as the fiducial age. We note that this ${\Delta t_{\rm form}}$ is of the order of or smaller than the one-sigma error bar on the age. However, to address formation by gravitational instability, where we expect the planet to be approximately coeval with the star, we also study the case of $t_{\rm cool}=14\pm 4$ Myr.

2.2 BEX cooling curves

For the MCMC we use the Bern EXoplanet cooling curves (BEX) with the AMES-COND atmospheres. The BEX models use the Bern planet evolution (cooling) code completo 21, which includes the cooling and contraction of the core and envelope at constant mass (see sects. 3.2 and 3.8.3 of Mordasini et al. 2012b, sect. 2.3 of Mordasini et al. 2012a, and sect. 2 of Linder et al. 2018) as well as deuterium burning (Mollière & Mordasini, 2012). The boundary conditions are provided by atmospheric models. Previously, only the simple Eddington model had been implemented, but we can now use arbitrary atmospheric models, following the coupling approach of Chabrier & Baraffe (1997). This entails simply taking a pressure–temperature point in the adiabatic part of the deep atmosphere as the starting point of the interior structure calculation. Since the structure is adiabatic, the precise location (e.g. at a Rosseland optical depth $\Delta{\tau_{\mathrm{R}}}=100$, at a pressure $P=50$ bar, or at the top of deepest convection zone) will not matter, and it is easy to verify that in any case the error in the radius is at most of a few percentage points. This coupling approach was applied recently to low-mass planets in Linder et al. (2018).

Currently, the BEX models are available with boundary conditions provided by

1. (i)

the Eddington assumption; 2. (ii)

AMES-COND (Baraffe et al., 2003); 3. (iii)

Burrows et al. (1997); 4. (iv)

petitCODE (Mollière et al., 2015); 5. (v)

HELIOS (Malik et al., 2017, Malik et al., in review).

For flavours (ii) and (iii) we extracted the relevant information from the publicly available Burrows et al. (1997) tracks and the Baraffe et al. (2003) grids111See https://www.astro.princeton.edu/~burrows/dat-html/data/ and https://phoenix.ens-lyon.fr/Grids/AMES-Cond/STRUCTURES/, respectively.. The details will be described in a dedicated publication, but we note already that we can reproduce very well the Burrows et al. (1997) and the AMES-COND tracks (see Fig. 1).

By default, the BEX curves assume full ISM deuterium abundance at the beginning of cooling, while in fact in core accretion a mass-dependent fraction will be burnt during formation (Mollière & Mordasini, 2012; Mordasini et al., 2017). However, in both cold 1 $M_{\odot}$ and 2 $M_{\odot}$ population syntheses, objects need a mass of 16 ${M_{\textnormal{J}}}$ (20 ${M_{\textnormal{J}}}$) to have consumed even only $\approx 30$ % ($\approx 70$ %) of their initial D abundance by the end of formation. Given the masses we find later for HIP 65426 b, and since in GI it is likely that no deuterium is destroyed during formation, the use of full deuterium abundance at the beginning of cooling is inconsequential.

We display in Fig. 1 the different flavours of the BEX cooling curves compared to the classical models of Burrows et al. (1997) and Baraffe et al. (2003). The initial luminosities are set to the same values as in Burrows et al. (1997), except for the 20 m$M_{\odot}$ case for which we took a slightly lower initial luminosity to avoid non-monotonicities in the re-interpolation of the original Burrows et al. (1997) data. Otherwise, the BEX curves clearly follow the Burrows et al. (1997) models, including the ‘shoulder’ that occurs during deuterium burning. At very old ages (20 Gyr) the black lines diverge because the models are beyond the tabulated range of input atmospheric structures.

We also see that the choice of either of the three classic atmospheric models (Eddington, AMES-Cond, Burrows) as outer boundary conditions only has a small effect on the cooling, as expected (Baraffe et al., 2003; Chabrier et al., 2000). Furthermore, it should be noted that the starting luminosities of the original AMES-COND (Baraffe et al., 2003) tracks are apparently not quite the same as for Burrows et al. (1997).

2.3 First analysis of the mass

Figure 2a compares HIP 65426 b to direct detections and the ‘hottest start’ cooling tracks of Mordasini et al. (2017), which use the simple Eddington outer boundary condition. A direct comparison with these cooling curves suggests a mass ${m_{\rm p}}\approx 8$–11 ${M_{\textnormal{J}}}$, which is not rare for direct detections of young isolated brown dwarfs (in the sense of substellar-mass objects; see the mass histogram in fig. 18 of Gagné et al. 2015). As shown in Fig. 1, these simpler models quite closely match cooling tracks based on detailed atmospheric models such as Burrows et al. (1997) or Baraffe et al. (2003). However, the luminosity error bar (0.1 dex, also a typical size; see e.g. Bowler, 2016) is small enough for the derived mass to depend slightly on the choice of the cooling curves.

After this first estimate of the mass based on models with an arbitrarily high post-formation luminosity, we look at cooling curves whose post-formation (also termed initial) luminosity $L_{\rm pf}$ follows the four relations seen in the population syntheses of Mordasini et al. (2017, their sect. 5.2.2 and their fig. 13). For ${m_{\rm p}}\approx 0.3$ to $\approx 12$ ${M_{\textnormal{J}}}$ (i.e. for planets that are massive enough to undergo the detached phase during the presence of the nebula, but not massive enough for deuterium burning to occur), these relations are given by

[TABLE]

respectively, for the hottest, cold-nominal, cold-classical, and coldest planets. Briefly, $L_{\rm pf}^{\textrm{hottest}}$ traces the brightest planet at every mass; $L_{\rm pf}^{\textrm{cold-nom.}}$ corresponds to the cold-nominal population, in which gas is assumed to accrete cold; $L_{\rm pf}^{\textrm{cold-class.}}$ is the best fit to the cold-classical population (which however shows an appreciable spread in luminosity at a given mass), in which the core artificially stops growing in the runaway phase à la Marley et al. (2007); and finally, $L_{\rm pf}^{\textrm{coldest}}$ traces the coldest planets at a given mass, which come from the small-core (coldest-start) population. It should be noted that we defined here the cold-nominal relation (Eq. (1b)) not as the mean of the cold-nominal population (as in Mordasini et al. 2017, with $L_{\rm pf}=1.2\times 10^{-5}L_{\odot}\,({m_{\rm p}}/{M_{\textnormal{J}}})^{1.3}$), but as the approximate upper envelope of points of that population.

Cooling tracks from all four relations are shown in Fig. 2b. At this age and for this mass there is barely any difference in the cooling curves of the hottest and the cold-nominal starts. However, the luminosities in the cold-classical population are one order of magnitude lower, the initial cooling (Kelvin–Helmholtz) timescale being $t_{\textrm{KH}}\sim 100\leavevmode\nobreak\ \textrm{Myr}$, which is roughly ten times longer than the age of HIP 65426 b. The coldest starts, finally, are even several orders of magnitude fainter than the others, with at the lower masses an initial $t_{\textrm{KH}}\sim 500$ Myr. Since the initial luminosities are well below the observed luminosity, no coldest start can match the observed luminosity of HIP 65426 b. Deuterium burning in more massive objects would be required here to reproduce the observed luminosity. In any case, as argued from different points of view by Mordasini et al. (2017) and Berardo et al. (2017), the coldest starts are not expected to be realistic.

We conclude that in this simple analysis, only the hottest and cold-nominal populations can reproduce HIP 65426 b. In the next section, we revisit this analysis in a more systematic fashion and take the error bars on age and luminosity into account.

2.4 Inputs: priors on mass and luminosity

Recently, using the tool of population synthesis, Mordasini et al. (2017) presented the first discussion of the statistics of planetary luminosities as predicted by a planet formation model. They looked in particular at the core accretion paradigm (Pollack et al., 1996; Mordasini et al., 2012b) and considered three populations, differing in the assumed efficiencies of the accretional heating of gas and planetesimals during formation:

1. (i)

a cold-nominal population, in which the entire gas accretion luminosity is radiated away at the shock, as in Mordasini et al. (2012a); 2. (ii)

a hot population, which differs from the first only by the assumption that the entire accretion luminosity is brought into the planet; 3. (iii)

a cold-classical population, which assumes, as in the classical work by Marley et al. (2007), that planetesimal accretion stops artificially once a giant planet enters the disc-limited gas accretion (detached) phase, and also does not include planetary migration.

Since the cold-classical population serves rather for model comparisons, and given that first dedicated and systematic simulations of the accretion shock have been recently performed (Marleau et al., 2017; Marleau et al., in prep.) but not yet used to produce cooling curves, we consider in this work the cold-nominal and hot populations as more realistic extreme scenarios.

We now turn to the total distribution function, which we write as

[TABLE]

Mordasini et al. (2017) showed that there is spread of post-formation entropies of approximately $\Delta s_{\textrm{pf}}\approx 1\leavevmode\nobreak\ k_{\textrm{B}}\,\textrm{baryon}^{-1}$ at a given mass (see their fig. 12), coming mostly from the core-mass effect (Mordasini, 2013; Bodenheimer et al., 2013). Given that the distribution of entropies is rather uniform for a given mass, we fit simple mass-dependent top-hat functions to the probability distributions of $s_{\textrm{pf}}$:

[TABLE]

The following functions, dropping the usual entropy units $k_{\textrm{B}}\,\textrm{baryon}^{-1}$, closely fit the envelope of points $s_{\textrm{pf}}({m_{\rm p}})$ in Mordasini et al. (2017). For the cold-nominal population, the lower and upper edges are given respectively by

[TABLE]

where masses ${m_{\rm p}}$ are implicitly in Jupiter masses in these equations, while for the hot population, the bounds are

[TABLE]

This holds down to 2 ${M_{\textnormal{J}}}$. We point out that in the Bern planet formation code, as in most codes using the Saumon et al. (1995) equation of state, there is an entropy offset relative to the published Saumon et al. (1995) tables (see Appendix B in Marleau & Cumming, 2014). This difference has no physical meaning, but care must be taken when comparing to work using codes with other entropy reference points such as MESA (Paxton et al., 2011, 2013, 2015) as used by Berardo et al. (2017).

Marginalising over entropy, the mass function in both the cold- and hot-start populations is, for about 1 to 10 ${M_{\textnormal{J}}}$, approximately given by

[TABLE]

i.e. the distribution is nearly flat in $\log{m_{\rm p}}$. As mentioned by Mordasini et al. (2017), this is similar to the distribution found by Mordasini et al. (2009) for population synthesis planets detectable by radial velocity, which in turn agreed with the ${\rm d}N/{\rm d}{m_{\rm p}}\propto{m_{\rm p}}^{-1.05}$ fit of Marcy et al. (2005). We note, however, that Cumming et al. (2008) found ${\rm d}N/{\rm d}{m_{\rm p}}\propto{m_{\rm p}}^{-1.3\pm 0.2}$ but for periods $<2000$ days, while Brandt et al. (2014) obtained from direct imaging ${\rm d}N/{\rm d}{m_{\rm p}}\propto{m_{\rm p}}^{-0.7\pm 0.6}$ at distances $\gtrsim 10$ au. Larger numbers of log-period radial velocity and direct-imaging detections will be necessary to reduce the error bars on these exponents.

2.5 Results: Mass–entropy constraints

Figure 3 shows the joint constraints on the mass and post-formation (or initial) entropy using the different priors discussed above. Considering first the case of uniform priors (i.e. not using information from formation scenarios), we find that the post-formation entropy $s_{\textrm{pf}}\geqslant\leavevmode\nobreak\ 9.2\leavevmode\nobreak\ k_{\textrm{B}}\,\textrm{baryon}^{-1}$ but it is not otherwise constrained. This lower limit holds independently of the formation pathway and for masses up to ${m_{\rm p}}\approx 15\leavevmode\nobreak\ {M_{\textnormal{J}}}$ (a conservative assumption). Marginalising instead over entropy, the 68.3% confidence interval (which is used throughout this section despite the non-Gaussianity of the posteriors) on the mass is ${m_{\rm p}}=9.6\pm 1.7\leavevmode\nobreak\ {M_{\textnormal{J}}}$. For high values of $s_{\textrm{pf}}$, the BEX models using the AMES-COND boundary conditions closely match the Baraffe et al. (2003) hot-start cooling tracks for these masses. If we consider somewhat arbitrarily $s_{\textrm{pf}}\gtrsim 14$ to approximate what is usually thought of as hot starts, we find ${m_{\rm p}}=9.0^{+1.3}_{-1.5}\leavevmode\nobreak\ {M_{\textnormal{J}}}$ for a cooling age $t_{\rm cool}12\pm 4$ Myr. This agrees well with the ${m_{\rm p}}=10\pm 2\leavevmode\nobreak\ {M_{\textnormal{J}}}$ reported by Chauvin et al. (2017) for the DUSTY models222 The value ${m_{\rm p}}=7_{-1}^{+2}\leavevmode\nobreak\ {M_{\textnormal{J}}}$ quoted by Chauvin et al. (2017) for COND03 (Baraffe et al., 2003) does not come from a luminosity comparison and is therefore less robust. However, COND03 and DUSTY use by construction the same luminosity tracks (Baraffe et al., 2003).. As expected from Marleau & Cumming (2014), the relative uncertainty on the hot-start mass $\sigma_{m_{\rm p}}/{m_{\rm p}}\approx 0.2$ is $\approx\frac{1}{2}\sigma_{t_{\rm cool}}/t_{\rm cool}\approx 0.3$.

Next we fold in the outcome of the population syntheses into the analysis. If we take only the mass prior (Eq. (8)) into account, we obtain ${m_{\rm p}}=9.4^{+1.5}_{-2.0}\leavevmode\nobreak\ {M_{\textnormal{J}}}$, which is lower by $\Delta{m_{\rm p}}\approx 0.2\leavevmode\nobreak\ {M_{\textnormal{J}}}$ than the case without priors. Applying the $s_{\textrm{pf}}$ priors (Eq. (5)) as well, we obtain ${m_{\rm p}}=10.9^{+1.4}_{-2.0}\leavevmode\nobreak\ {M_{\textnormal{J}}}$ for the cold population (Eq. (6)) and ${m_{\rm p}}=9.9^{+1.1}_{-1.8}\leavevmode\nobreak\ {M_{\textnormal{J}}}$ for the hot population (Eq. (7)).

These masses are shown as points with horizontal bars in the main panel of Fig. 3. The difference between the mass inferred with and without the mass priors is small, with $\Delta{m_{\rm p}}\lesssim 0.2\leavevmode\nobreak\ {M_{\textnormal{J}}}$. These differences represent only a modest fraction of the error bars. However, the $s_{\textrm{pf}}$ priors are mildly important, leading to a difference $\Delta{m_{\rm p}}\approx 1\leavevmode\nobreak\ {M_{\textnormal{J}}}$ between the hot and the cold populations and even $\Delta{m_{\rm p}}=1.5\leavevmode\nobreak\ {M_{\textnormal{J}}}$ between the flat-prior and (with Eq. (2)) the cold-population cases. Finally, we note the distinctly asymmetrical shape of the confidence intervals when using the priors. This asymmetry comes mostly from the $s_{\textrm{pf}}$ prior despite the $p_{m_{\rm p}}\propto{m_{\rm p}}^{-1}$ scaling since the mass interval is small.

The posterior on the post-formation entropy changes dramatically when taking the population-synthesis priors into account, as visible in the right panel of Fig. 3. The lower bound $s_{\textrm{pf}}\gtrsim 9.2$ obtained with the uniform prior does not change, but the population-synthesis priors lead to the determination of an upper bound, yielding $s_{\textrm{pf}}=10.4^{+0.7}_{-0.2}$ in the case of the hot population and $s_{\textrm{pf}}=10.2^{+0.3}_{-0.7}$ for the cold. It should be noted that the probability maxima are rather flat. These values differ only marginally from each other, reflecting the large overlap between the post-formation entropies or luminosities of the cold- and hot-start populations, which is ultimately a consequence of the core-mass effect (CME) as discussed by Mordasini et al. (2017, sect. 5.2.1). These values $s_{\textrm{pf}}\approx 10.3$ are clearly lower than classical (arbitrarily) hot starts ($s_{\textrm{pf}}\approx 13$), with an initial Kelvin–Helmholtz time $t_{\textrm{KH}}\sim 10$ Myr as opposed to $t_{\textrm{KH}}\lesssim 1$ Myr for classical hot starts. Thus, HIP 65426 b would have just begun joining the hot-start cooling track (see Marleau & Cumming, 2014 for a general discussion of the shape of cooling tracks). We finally note that the mass prior barely changes the $s_{\textrm{pf}}$ posteriors.

2.6 Discussion

For comparison, with a shorter cooling age $t_{\rm cool}=10\pm 4$ Myr (i.e. coming from a longer formation period), we obtain with uniform priors ${m_{\rm p}}=9.0^{+1.9}_{-1.7}\leavevmode\nobreak\ {M_{\textnormal{J}}}$ and with only the mass prior ${m_{\rm p}}=8.7^{+1.9}_{-2.0}\leavevmode\nobreak\ {M_{\textnormal{J}}}$, whereas using the mass and $s_{\textrm{pf}}$ priors from hot-start (cold-start) populations yields ${m_{\rm p}}=9.3^{+1.3}_{-1.8}\leavevmode\nobreak\ {M_{\textnormal{J}}}$ (${m_{\rm p}}=10.3^{+1.6}_{-1.8}\leavevmode\nobreak\ {M_{\textnormal{J}}}$). Instead, taking a cooling age $t_{\rm cool}=14\pm 4$ Myr, i.e. the age of HIP 65426 as might correspond to formation by gravitational instability, with only mass priors we obtain ${m_{\rm p}}=9.9^{+1.4}_{-1.7}\leavevmode\nobreak\ {M_{\textnormal{J}}}$; instead, using the $s_{\textrm{pf}}$ and mass priors from the hot-start (cold-start) populations yields ${m_{\rm p}}=10.4^{+1.0}_{-1.6}\leavevmode\nobreak\ {M_{\textnormal{J}}}$ (${m_{\rm p}}=11.3^{+1.0}_{-1.7}\leavevmode\nobreak\ {M_{\textnormal{J}}}$). This is somewhat higher than, but still consistent with, the mass derived by Cheetham et al. (2019). Using an age of $14\pm 4$ Myr and the AMES-Cond models, they found ${m_{\rm p}}=7.5\pm 0.9\leavevmode\nobreak\ {M_{\textnormal{J}}}$ based on magnitudes in individual bands and ${m_{\rm p}}=8.3\pm 0.9\leavevmode\nobreak\ {M_{\textnormal{J}}}$ based on their bolometric luminosity, which had an uncertainty $\sigma_{\log L}=0.03$ dex half as large as the value used here. That the mass found by Cheetham et al. (2019) is lower than that derived here is not surprising since the AMES-Cond models they used correspond only to hotter starts, whereas here a range of $s_{\textrm{pf}}$ was considered.

In general, one could expect somewhat different results if using the logarithm of the post-formation luminosity instead of the post-formation entropy as an independent variable (along with the mass). Indeed, the luminosity $L$ and entropy $s$ are monotonic functions of each other at a given mass, but the slope ${\rm d}\log L/{\rm d}s$ depends on both mass and entropy (Marleau & Cumming, 2014). This means that a prior which is uniform in $s$ for all masses is not uniform in $\log L$ for all masses, and vice versa.

However, one can argue that this should be of negligible concern. In the case of a flat prior in $s_{\textrm{pf}}$, the posterior was also relatively flat, and a small distortion will not change the nature of the weak constraints on $s_{\textrm{pf}}$. The distortion should be small judging by the precise scalings identified in eq. (9) of Marleau & Cumming (2014), and while these hold specifically for their Eddington atmospheric models, the $L(s)$ relation will not be entirely different for AMES-COND. In the case of the hot or cold priors, the posteriors are non-zero over a relatively small region, so that in this case too there should not be any significant skew.

We finally note that, as mentioned in Section 2.1, Fig. 3 shows that HIP 65426 b is unlikely to have a mass for which a meaningful fraction of deuterium could be burnt (cf. Spiegel et al., 2011; Mollière & Mordasini, 2012). This justifies a posteriori the use of cooling curves that assume full deuterium abundance at the start. In the case of other detections close to the deuterium-burning limit (${m_{\rm p}}\approx 13\pm 2\leavevmode\nobreak\ {M_{\textnormal{J}}}$), however, this assumption would need to be revisited if they formed over a longer timescale than expected from gravitational instability.

To summarise, we find that the mass of HIP 65426 b is ${m_{\rm p}}=10.9^{+1.4}_{-2.0}\leavevmode\nobreak\ {M_{\textnormal{J}}}$ with priors from the cold population and ${m_{\rm p}}=9.9^{+1.1}_{-1.8}\leavevmode\nobreak\ {M_{\textnormal{J}}}$ using the hot population.

3 Forming HIP 65426 b in core-accretion models

We now switch from the study of HIP 65426 b’s post-formation thermodynamical evolution to numerical experiments concerning its formation. Core accretion models typically involve forming cores of giant planets and having them undergo runaway gas accretion at small orbital radii (${a_{\rm p}}\lesssim 20$ au) before they migrate in towards the central star, inconsistent with the location of HIP 65426 b. At larger orbital radii, the time taken to form a core through planetesimal accretion is longer than typical protoplanetary disc lifetimes, though forming a core through pebble accretion could be significantly faster (Lambrechts & Johansen, 2014; but see also Rosenthal & Murray-Clay, 2018). Even if a single planetary core is able to form at large orbital radii, either through pebble or planetesimal accretion, interactions with the local protoplanetary disc will force the planet to migrate through type I migration to small orbital separations on timescales shorter than that required for the core to accrete a significant gaseous envelope and undergo runaway gas accretion (Coleman & Nelson, 2016b). This fast migration poses the main problem for forming a planet that has properties consistent with that found for HIP 65426 b.

To overcome these problems for the core accretion model in forming planets such as HIP 65426 b, we ran numerous $N$-body simulations in which we placed a number of giant planet cores in a protoplanetary disc and allowed them to mutually interact, migrate throughout the disc, and accrete gaseous disc material. The idea is that as one giant planet core undergoes runaway gas accretion and rapidly increases its mass, the system of planets becomes dynamically unstable, leading to the scattering of one of the less massive cores. This core, once scattered out into the outer disc will then circularise its orbit and begin to migrate back in towards the central star. However the core will continue accreting gas from the surrounding disc and could then undergo runaway gas accretion in the outer disc, becoming a gas giant and transitioning to the slower type II migration regime. If the planet is scattered out far enough and has insufficient time to migrate back in towards the inner disc, its final mass and semi-major axis could be similar to those of directly imaged planets and HIP 65426 b. This process has been observed in population synthesis models when, in massive discs, multiple gas giants form as a first generation and subsequently destabilise the orbits of surrounding embryos, scattering them to larger orbits where they then grow into gas giants (Ida et al., 2013).

3.1 Simulation set-up

In order to run these simulations, we adapted the $N$-body and disc model of Coleman & Nelson (2016b) to be appropriate for a protoplanetary disc surrounding an A-type star such as HIP 65426. This model couples a 1D thermally evolving viscous disc model (Shakura & Sunyaev, 1973) to the Mercury-6 symplectic integrator (Chambers, 1999) and includes prescriptions for photoevaporation (Dullemond et al., 2007; Alexander & Armitage, 2009), both type I and II planet migration (Paardekooper et al., 2011; Lin & Papaloizou, 1986), and gas accretion from the surrounding disc (Coleman et al., 2017). Table 1 shows the disc parameters used for the simulations. We chose the values for the viscosity parameter $\alpha$ and the photoevaporation factor $\Phi_{41}$ to give the disc an appropriate lifetime. Using the values presented in Table 1, the initial disc had a total mass equivalent to $\sim 8\%$ of the mass of HIP 65426 (i.e. around $150\leavevmode\nobreak\ {M_{\textnormal{J}}}$) and a lifetime of 3.5 Myr. The lifetime of the disc in the simulations is always shorter since a significant fraction of the total gas mass is accreted onto the planets.

Since type I migration timescales for giant planet cores are shorter than the timescales for the cores to reach runaway gas accretion (Coleman & Nelson, 2016b), we require a mechanism to stall type I migration (see also Pudritz et al., 2018). To stall type I migration and counter the short timescales experienced by giant planet cores, we placed a radial structure in the disc that mimics the effects of a zonal flow. Zonal flows have been observed in both local (Johansen et al., 2009) and global (Steinacker & Papaloizou, 2002; Papaloizou & Nelson, 2003; Fromang & Nelson, 2006) numerical simulations of magnetised discs, including those incorporating non-ideal MHD effects (Bai & Stone, 2014; Béthune et al., 2016; Zhu et al., 2014). Radial structures which could also be reminiscent of zonal flows have also been seen in numerous observations of protoplanetary discs (ALMA Partnership et al., 2015; Andrews et al., 2016; van Boekel et al., 2017). The effect of zonal flows on a protoplanetary disc is to create a radial pressure bump in the disc, which results in a positive surface density gradient. This positive surface density gradient increases the strength of the vortensity component of a planet’s corotation torque, allowing it to balance the planet’s Lindblad torque, thus creating a planet trap that stalls type I migration (Masset et al., 2006; Hasegawa & Pudritz, 2011; Coleman & Nelson, 2016a). To account for a radial structure in the disc that mimics the effects of a zonal flow, we included a single radial structure in each simulation, following the approach used in Coleman & Nelson (2016a, see their section 2.3.3). This radial structure increases the local $\alpha$ parameter when calculating the viscosity, which results in a reduction in the local surface density, creating a positive surface density gradient that acts as a planet trap as described above. We assume that this structure remains at the same location in the disc, placed arbitrarily at either 15 or 20 au in our simulations, and has a lifetime equivalent to the disc lifetime. The lifetimes of zonal flows in MHD simulations are still unexplored due to long simulation run times, but since these structures are seen in both young and old protoplanetary disc observations (ALMA Partnership et al., 2015; Andrews et al., 2016), it seems reasonable to assume that the flows are long lived.

To account for planet migration we use the torque formulae of Paardekooper et al. (2010, 2011) whilst the planet is embedded in the disc, to simulate type I migration due to Lindblad and corotation torques. Our model accounts for the possible saturation of the corotation torque (Paardekooper et al., 2011), and also the influences of eccentricity and inclination on the disc forces (Cresswell & Nelson, 2008; Fendyke & Nelson, 2014). Once the planet has become massive enough to open a gap in the disc we use the impulse approximation to calculate the torques acting on the planet from the surrounding disc as it undergoes type II migration (Lin & Papaloizou, 1986). To calculate gas accretion on to the planet, we use the accretion routine presented in Coleman et al. (2017). In this model, whilst the planet is embedded in the disc we construct a 1D envelope structure model that self-consistently calculates the gas accretion rate taking into account local disc conditions. After the planet has opened a gap in the disc, we assume that the gas accretion is equal to the viscous supply rate. All gas that is accreted on to the planet is removed from the surrounding disc.

For each simulation we placed five planets of masses $15\,{\rm M}_{\oplus}\leq{m_{\rm p}}\leq 20\,{\rm M}_{\oplus}$ at ${a_{\rm p}}=15$–30 au, i.e. in the outer disc beyond the radial structure, in close proximity to each other (initial period ratios between neighbouring planets ranging from 1.08 to 1.7). We placed the planets in close proximity to each other to ensure that they were able to become trapped in resonant chains fairly quickly before a single core could undergo runaway gas accretion and thus destabilise the system. This configuration is frequently seen to arise in global planet formation simulations that include planet migration, planetesimal accretion, mutual interactions between planetary embryos, and evolution of the protoplanetary disc (Coleman & Nelson, 2016b, a). However, due to the chaotic nature of the formation processes (i.e. migration, planetesimal accretion rates, $N$-body interactions), we force this initial set-up onto the planets for these simulations so as to save on computational time. We also varied the location of the radial structure as described above and the formation time of the giant planet cores which ranged between 1.5 and 2.5 Myr. These different initial conditions led to the computation of 792 simulations.

3.2 Example HIP 65426-like simulated system

Figure 4 shows the mass versus orbital distance evolution (left panel) and the temporal evolution of planet masses, semi-major axes, and eccentricities (right panel) of a typical example of such a simulation. The mass versus orbital distance tracks of the planets are shown with solid lines indicating semi-major axes and dashed lines displaying the planets’ pericentres and apocentres. Black dots represent the final masses and semi-major axes of the planets with the red cross showing the mass and orbital distance of HIP 65426 b as discussed in Sect. 2. As the simulation starts, all of the planets begin to accrete gas and migrate inwards towards the radial structure. The planets’ migration stalls as they approach the radial structure (see label A in Fig. 4) due to the enhanced corotation torques arising from the radial structure’s effect on the local disc profile. The outer three cores (blue, green, and yellow lines) then undergo runaway gas accretion, opening a gap in the disc (see the sudden increase in mass for some of the planets at $\sim 2.8$ Myr in the top right panel of Fig. 4). The inner planets (purple and orange lines) are then starved of gas due to the opening of gaps in the disc, delaying their ability to transition to runaway gas accretion. After a further 50 kyr, the system becomes unstable with two of the giants impacting each other. Other cores are also scattered, with one core (purple line) being scattered to ${a_{\rm p}}\approx 220$ au with an eccentricity ${e_{\rm p}}\approx 0.85$ (label B in Fig. 4). This core then undergoes runaway gas accretion at this large orbital distance, accreting gas when it enters the disc on its eccentric orbit (label C in Fig. 4). When the planet’s mass has risen to ${m_{\rm p}}\approx 2\leavevmode\nobreak\ {M_{\textnormal{J}}}$, it has a close encounter (label D in Fig. 4) with the most massive giant in the system (with a mass ${m_{\rm p}}\approx 4\leavevmode\nobreak\ {M_{\textnormal{J}}}$, shown by the blue line). This lowers the semi-major axis and eccentricity of the former to ${a_{\rm p}}\approx 125$ au and ${e_{\rm p}}\approx 0.7$ respectively (see the drops in semi-major axis and eccentricity at $\sim 2.9$ Myr in the middle and bottom right panels of Fig. 4). This planet then continues to accrete gas and slowly migrate in towards the central star (label E in Fig. 4), with the disc gradually damping the planet’s eccentricity. We implement eccentricity damping for giant planets by setting the damping timescale to 100 local orbital periods. This timescale is consistent with eccentricity damping timescales found for eccentric planets in isothermal discs (Bitsch et al., 2013).

By the time the disc has fully dispersed, the planet has grown to ${m_{\rm p}}\approx 9.8\leavevmode\nobreak\ {M_{\textnormal{J}}}$ and has migrated in to having a semi-major axis of ${a_{\rm p}}\approx 77$ au. The temporal evolution of the planets’ semi-major axis, mass and eccentricity can be seen in the right panel of Fig. 4, with the shaded grey region indicating the times during which the disc is present. Due to the circularisation of its orbit, the planet’s eccentricity has dropped to ${e_{\rm p}}\approx 0.27$, resulting in the planet orbiting between 56 and 98 au, spending most of its orbit near apocentre with an orbital distance greater than 90 au. This range in orbital distance is shown by the horizontal black bar in Fig. 4 and is compatible with the observed position of HIP 65426 b. Also remaining in the system are three other giant planets with semi-major axes of ${a_{\rm p}}=9.4$, 15.6, and 25.8 au, and masses ${m_{\rm p}}=0.45$, 2.8, and 7.8 ${M_{\textnormal{J}}}$, respectively, on nearly circular orbits (${e_{\rm p}}\lesssim 0.05$ for all).

3.3 Overall results

The system described above, with a giant planet similar to HIP 65426 b and numerous planets with shorter periods, is a common outcome of these simulations. The left panel of Fig. 5 shows the final masses and semi-major axes from the simulations that had a similar outcome to that described above. We define these systems as containing at least one planet with semi-major axis ${a_{\rm p}}>50$ au and mass ${m_{\rm p}}>1\leavevmode\nobreak\ {M_{\textnormal{J}}}$. As can be seen, there are numerous giant planets that are similar to HIP 65426 b (shown by the red diamond) in terms of mass and semi-major axis (projected orbital distance for HIP 65426 b), but with a wide range of eccentricities, spanning essentially ${e_{\rm p}}=0$–1. All of these planets are accompanied by a number of interior giant companions that could be detectable in long-baseline radial-velocity surveys.

The right panel of Fig. 5 shows the final semi-major axes and masses of all surviving planets in the simulations. The colour of the marker indicates each planet’s eccentricity, and the red diamonds show the currently observed planets found in direct imaging surveys. There is good agreement between the observations and simulated giant planets in terms of semi-major axes and planet mass. The eccentricities cannot really be compared since for observed planets they are not well constrained due to insufficient time sampling of their long orbital periods (Bowler & Nielsen, 2018). The simulated giant planets have non-zero eccentricities which in many cases are significant (${e_{\rm p}}>0.5$), with most of these high-eccentricity giants having semi-major axes greater than 100 au. This is not surprising since for these planets to attain such large semi-major axes, they need to undergo significant scattering, which induces high eccentricities, and since the planets have had little time to migrate back in towards the central star, their orbits have also had insufficient time to circularise.

The plot also shows that for the distant giant planets, eccentricity and orbital distances are positively correlated, an imprint of that planet’s main scattering event. This is expected from the fact that the original formation region ($\sim 20$–30 au) remains, at least for a single scattering event, part of the orbit as its pericentre distance. In the simulations, however, eccentricity damping from the gas disc and minor interactions with other planets in the system can decrease the planet’s eccentricity over time, raising the pericentre away from the formation region. Since the distant planets that formed in the simulations had insufficient time to circularise fully, due to dispersal of the gas disc, this imprint of the main scattering event remains, explaining the distance–eccentricity correlation.

Also seen in Fig. 5 are numerous giant planets with semi-major axes ${a_{\rm p}}<10$ au. Quite often these giant planets were responsible for scattering a giant planet core into the outer system that could then undergo runaway gas accretion at hundreds of astronomical units, as was described in the example simulation above (Section 3.2). While these planets are typically too faint and too close to the star to be observed in direct imaging surveys, they could be observed in radial velocity surveys (Butler et al., 2017) or in astrometry surveys such as GAIA (Gaia Collaboration et al., 2016, 2018). Further observations of HIP 65426 using the radial velocity or astrometry technique could yield additional giant planets in the system closer to the star than HIP 65426 b. Very recently, Cheetham et al. (2019) ruled out to 5 $\sigma$ the presence of further companions more massive than 16 ${M_{\textnormal{J}}}$ down to orbital separations of 3 au. This is consistent with the simulated planets shown in Fig. 5, which all fall below the detection limits out to 30 au. Most of the inner companions to HIP 65426 b-like planets have masses well below 10 ${M_{\textnormal{J}}}$. Should these planets exist, and if HIP 65426 b were found to have an eccentric orbit, this could suggest the formation origin of HIP 65426 b as described here, and may indicate that other directly imaged giant planets should have giant planet companions closer to their host star than has been observed up to now.

For systems that contained two or more giant planets at the end of the disc lifetime, it is possible that dynamical instabilities between giant planets as the systems age lead to the planets having wider orbits, similar to HIP 65426 b. Figure 6 shows the planet mass versus semi-major axis evolution (left panel) and the temporal evolution of planet mass, semi-major axis, and eccentricity for such a scenario. Here the planets undergo a similar initial evolution to that described in Sect. 3.2, but as the disc fully disperses, all five giant planets have relatively stable orbits (given the eccentricity-damping effect of the gas) with ${a_{\rm p}}=8$–70 au. However these orbits are not stable on long timescales after disc dispersal, and within 0.2 Myr of the disc dispersing, the planets undergo significant dynamical instabilities, increasing eccentricities and scattering some of the giants to larger semi-major axes. Continued interactions resulted in three of the giant planets being ejected from the system, with the two most massive giant planets remaining. These surviving planets are shown by the black dots in the left panel of Fig. 6, where the solid black horizontal lines shows the extent of their orbit from pericentre to apocentre. A more detailed study of forming HIP 65426 b through giant–giant scattering is discussed in Sect. 4.

4 Post-formation scattering of giant planets

In this section we explore the formation of systems with giant planets on wide orbits, such as HIP 65426 b, through planet–planet scattering after disc dispersal. This mechanism is known to create highly eccentric planets (Rasio & Ford, 1996; Lin & Ida, 1997; Chatterjee et al., 2008; Marzari & Weidenschilling, 2002). We estimate here its efficiency in raising the apocentre of a giant planet above $\sim 100$ au without ejecting it. We examine two scenarios: two-planet scattering and three-planet scattering. A system with two planets behaves qualitatively differently than a system with three or more planets.

On initially coplanar circular orbits, if the initial semi-major axes of two planets are closer than (Wisdom, 1980; Deck et al., 2013)

[TABLE]

their orbits will be unstable on short timescales, typically of order $\tau\sim((m_{1}+m_{2})/m_{\star})^{-1/3}P_{1}$ (Petit et al., 2017), where $P_{k}$, $a_{k}$, and $m_{k}$ are the orbital periods, semi-major axes, and the masses of the two planets, and $m_{\star}$ is the mass of the central star. Once such a system undergoes an instability, the most probable outcome is a single-planet system (Ford & Rasio, 2008). Being in the unstable area given by Equation (9) after disc dispersal requires one of two scenarios:

1. (a)

the planets have migrated into a stable configuration such as a 1:2 mean-motion resonance (MMR) during the disc phase (Lee & Peale, 2002), and this stable configuration was disrupted after disc dispersal; 2. (b)

the planets were in the unstable area during their formation, but the significant disc mass postponed the instability, for instance through eccentricity damping (Ford & Rasio, 2008).

In the case (b), instability can ensue before the total dispersal of the disc, which might still affect the orbital evolution of the giant planets. The results of this section thus have to be compared to the occurrences of giant–giant scattering during the disc phase (Section 3).

On the other hand, a system with three or more planets does not have such a sharp stability boundary as in Equation (9). These systems can become unstable for much larger initial spacings. However, as the initial spacings increase, the timescale of the first close encounter increases as well (Chambers et al., 1996). Single giant planets or pairs are common outcomes of this instability (Chatterjee et al., 2008), as seen for example in Fig. 6.

The two- and the three-planet-scattering scenarios are both consistent with observational constraints. Out of the hundreds of giant planets of mass above $2\leavevmode\nobreak\ {M_{\textnormal{J}}}$ that have been observed with semi-major axes ranging from 1 to 20 au, tens are known to belong to multi-planetary systems containing at least two giant planets444See exoplanet.eu. However, these statistics being incomplete due to observational biases, the multiplicity of giant planet systems is probably underestimated..

To explore these two scenarios, we performed $N$-body simulations of HIP 65426-like systems after disc dispersal. We used the variable-step integrator DOPRI, whose behaviour for highly eccentric orbits was validated in a previous work (Leleu et al., 2018). We integrated the synthetic systems for $5\times 10^{6}$ years, which is comparable to the age of the system since disc dispersal (see Section 2.1). Alternatively, integrations were stopped when only one planet remained in the system. The mass of the star was set to $m_{\star}=2\leavevmode\nobreak\ M_{\odot}$, while the mass of each planet was randomly picked in the interval $m_{k}=5$–15 ${M_{\textnormal{J}}}$. The radius $R_{k}$ of each planet $k$ was set to $R_{k}=\left(1.1+0.06\times m_{k}/{M_{\textnormal{J}}}\right)\leavevmode\nobreak\ {R_{\textnormal{J}}}$, which roughly fits the non-accreting hot population of Mordasini et al. (2012a) at 3–5 Myr555This can be explored with data from the Data Analysis Centre for Exoplanets (DACE) platform at https://dace.unige.ch, in the ‘Evolution’ section.. Planets that entered the Roche limit of the star were removed from the simulation, with the Roche limit given by $R_{\star,\textrm{Roche}}\approx 2.2\leavevmode\nobreak\ R_{\odot}$ (the stellar radius is $R_{\star}=1.77\leavevmode\nobreak\ R_{\odot}$; Chauvin et al., 2017) for our considered range of planetary masses and radii. Collisions between planets were detected when the physical radii of the two objects intersect and were treated as completely inelastic, i.e. assuming perfect merging and conservation of total momentum and mass. Other collision models including possible hit-and-runs and energy dissipation might change the outcomes of the simulations slightly. However, they should not create a significant number of broader orbits as these collisions typically reduce the eccentricities of the bodies.

Initial eccentricities were set to ${e_{\rm p}}=0$ and will be excited though planet–planet interactions. As we are primarily interested in the feasibility of raising the apocentre of a planet above $\sim 100$ au, we restrict our study to the coplanar case and set all inclinations to $i=0$. All other angular orbital elements where chosen randomly within [0:$360^{\circ}$]. The initial distribution of semi-major axes depends on the considered scenario, but they are generally taken in the 10–15 au range as it is the upper limit for the typical formation of giants in the core accretion scenario as mentioned above. Taking large initial semi-major axis makes it possible, through angular momentum transfer with other planets, to raise the apocentre of a given planet to greater values without being too close to ejection.

4.1 Two-planet scattering

4.1.1 Conservative case

In the two-planet-scattering scenario, the inner planet was positioned at ${a_{\rm p}}=10$ au, while the outer one was positioned slightly inside the instability domain (Eq. 9), at $a_{2}=a_{1}\left(1+1.42((m_{1}+m_{2})/m_{\star}\right)^{2/7}$. We integrated 500 systems with this set of initial conditions; the final outcomes are plotted in Figure 7.

In the conservative case, two planets on intersecting orbits will continue to experience close encounters until one of the three following outcomes happens: planet–planet collision, planet–star collision, or planet ejection. These events occurred within a few $10^{5}$ years, which is significantly shorter than the estimated age of the system. We now discuss each in turn.

Planet–planet collisions tend to decrease the eccentricity that the planets acquired during their stay in the unstable domain, while energy conservation ensures that the semi-major axis of the resulting planet $a_{\rm r}$ lies between the semi-major axes of the initial ones (Ford & Rasio, 2008):

[TABLE]

This is the most common outcome ($\approx 63\,\%$ of all systems), leading to the red clump between 10 and 15 au in Fig. 7. These semi-major axes are too small to correspond to HIP 65426 b.

In the case of ejection, the escaping planet typically leaves the system with a very low (positive) energy (Moorhead & Adams, 2005). The orbit of the remaining planet is hence predictable using angular momentum and energy conservation, yielding

[TABLE]

where the subscripts ‘r’ and ‘ej’ refer to the remaining and ejected planet, respectively, and where $q_{\rm ej}$ is the ejected planet’s minimal distance to the star on its parabolic orbit. This scenario represents almost all other cases ($\approx 36\,\%$ of the systems). The range of possible orbits for the typical values $m_{1}+m_{2}=20\leavevmode\nobreak\ {M_{\textnormal{J}}}$ and $q_{\rm ej}=10$ au is shown by the yellow line in Figs. 7 and 8. It closely matches the distribution from the $N$-body integrations. This scenario again leads to planets with orbital distances too small in comparison to HIP 65426 b.

The last outcome, planet–star collisions, is less likely for our range of initial conditions ($\approx 1\,\%$), but yields a wider range of final configurations. To fall onto the star, a planet initially on a circular orbit at $10$ au needs to give most of its angular momentum to the other planet. Depending on the mass ratio of the planets, this might not be enough to eject the outer planet, which therefore would remain on a wider orbit. Figure 7 shows that most of these remaining planets are consistent with the observed projected distance of HIP 65426 b. In this set-up (two planets, no tides) HIP 65426 b would be on a highly eccentric orbit, and it would be the only body in the system. We show later, when tides are included in the two-body scenario, that the outcomes explaining HIP 65426 b can again only contain the scattered body alone (the other was sent into the star), but also configurations with two remaining bodies, with the second companion very close to the star, circularised by tides.

4.1.2 Effect of dynamical tides

In the conservative case, we saw that almost all systems with two planets initially at $\approx 10$ au underwent close encounters and ended up in the planets colliding or one of them being ejected. However, even after gas dispersal, a planet that gets close enough to the star on its orbit, typically with a pericentre $q=a(1-e)\sim R_{\star}\sim 0.01\leavevmode\nobreak\ \textrm{au}$, will undergo tidal circularisation (Ivanov & Papaloizou, 2004). This process has the main effect of lowering the apocentre of an eccentric planet, while keeping the pericentre roughly constant, which can stop the orbits from intersecting before the ejection of the outer planet.

Since the planets that experience tidal effects will be on wide eccentric orbits, we consider dynamical tides, which are a succession of tidal excitation (when the planet is close to its pericentre) and relaxation (during the rest of the orbit) (Ivanov & Papaloizou, 2004). For giant planets, the migration and eccentricity timescales of these tides can be below $10^{5}$ years (Nagasawa et al., 2008) and hence can be comparable to the lifetime of the two-planet systems integrated in the conservative case. To take these tides into account, we adopt the formula of Ivanov & Papaloizou (2004) who calculated analytically the strongest normal modes, the $l=2$ fundamental modes, of the tidal deformation. Depending on the rotation of the planet that undergoes tidal circularisation, they derived the tidally gained angular momentum ($\Delta L_{\rm tide}$) and energy ($\Delta E_{\rm tide}$) during a single pericentre passage in two extreme cases:

• •

an initially non-rotating planet, which tends to maximise the effect of the tides, with

[TABLE]

where $L_{k}=m_{k}(Gm_{k}R_{k})^{1/2}$ and $E_{k}=Gm_{k}^{2}/R_{k}$ are the angular momentum and energy scales, $\zeta_{k}=(m_{k}q_{k}^{3})^{1/2}(m_{\star}R_{k}^{3})^{-1/2}$, $\omega_{0}$ is the dimensionless frequency of the fundamental mode (normalised by the internal dynamical frequency $[Gm_{k}p/R_{k}^{3}]^{1/2}$), and $Q$ is a dimensionless overlap integral that depends on the planetary interior model;

• •

a planet that is spinning at the critical rotation rate, for which the passage at pericentre does not provide an increase in angular momentum. This minimises the effect of the tides:

[TABLE]

We translate either of these two expressions for $\Delta L_{\rm tide}$ and $\Delta E_{\rm tide}$ into migration and eccentricity damping timescales using respectively (Nagasawa et al., 2008)

[TABLE]

with $\gamma_{k}=(1-e_{k}^{2})/{e_{k}}^{2}$ and where $P$ is the orbital period of the planet.

We note that for these tidal models to be realistic, it is necessary that the normal modes arising near the pericentre passage be fully dissipated before the next pericentre passage, which is typically the case for a semi-major axis above a few astronomical units (Ivanov & Papaloizou, 2004). Moreover, the actual spin of the planet evolves over time, which causes the effectiveness of the tides to vary. As a result, we assume that this model allows us to correctly represent the evolution of the orbit of an inner planet during the early stages of its apocentre lowering, but does not represent correctly the final state of the inner planet. This is however not of concern as we are primarily interested in the final orbit of the outer planet.

To estimate the effect of the tides on the systems that were integrated in the conservative case, we re-ran the same set of initial conditions as in Fig. 7 in two cases: for weak tides, using the set of equations (13), and for strong tides, using the set of equations (12). The results are presented in the top and bottom panel of Fig. 8, respectively. We assumed that $\omega_{0}=1.2$ and $Q=0.56$ for all planets, as these dimensionless parameters tend to be independent of the radius of the planet for $m_{k}=5\leavevmode\nobreak\ {M_{\textnormal{J}}}$, and we assume that it remains the case for more massive planets (Ivanov & Papaloizou, 2004).

In both cases, the majority of the systems still evolve towards the two main outcomes of the conservative case: either ejection of one of the planets or planet–planet collision, leaving a single planet with a semi-major axis below 15 au. The relative occurrences are very similar to the conservative case. The similarity with the conservative case is easily understandable as the tides affect systems for which the pericentre of one of the planets goes below a few hundredths of an astronomical unit, which is relevant only for a few systems. Nevertheless, including tides increases the number of systems that exhibit a planet–star collision or that retain two planets until the end of the simulation. The ‘weak tides’ model produced more planets with large semi-major axes than did the ‘strong tides’ model. The reason is that the latter tends to lower the apocentre of the inner planet before it can exchange enough angular momentum with the outer planet. It is important to note that the number of planet–star collisions is considerably overestimated due to our continuous application of dynamical tides even when the apocentre of the inner planet goes below a few astronomical units. In that sense, most of these systems are more likely to retain a close-in circularised giant planet in addition to the outer ones displayed in orange in Fig. 8. Although the tides allow for a broader diversity of outcomes for the scattering of two giant planets, only a small fraction of the systems contain planets with apocentres above 90 au.

4.2 Three-planet scattering

As mentioned previously, in the three-planet case there is no sharp stability condition regarding the initial semi-major axes of the giant planets. Following Marzari & Weidenschilling (2002), we initially position $m_{2}$ at 10 au, and $m_{1}$ and $m_{3}$ at four mutual Hill radii inside and outside the orbit of $m_{2}$, respectively (this corresponds to spacings $\approx 5$–7 au). This initial spacing does not necessarily ensure instability in the system within the $5\times 10^{6}$ years of integration, and on the other hand these systems can become unstable even for much wider initial spacings, but the timescale of first encounter will increase as well (Chambers et al., 1996). We chose this spacing in order to have a good probability of close encounters within the age of the system (see next paragraph). We point out that our results are probably more general than our restricted set of initial conditions might suggest, since the time until instability for a particular set of initial conditions does not affect the statistical properties of final outcomes in this kind of study (Chatterjee et al., 2008).

Out of 1000 initial conditions, the planets strongly interacted in $46\,\%$ of the systems, which resulted in the loss of at least one planet within $5\times 10^{6}$ years. These systems are shown in Fig. 9. In the remaining systems, the planets oscillated around their initial semi-major axis without significant increase of eccentricity and will not be discussed further. Out of the systems that interacted, we separated those resulting in single-planet systems (top panel) and two-planet systems (bottom panel).

Single-planet systems generally underwent two ejections, a planet–planet collision and an ejection, or a planet–star collision and an ejection. Although the outcomes are less predictable than in the two-planet case, it is still planet–star collisions that tend to allow a single planet to remain on a wide orbit after the removal of its companions.

Systems that lost only one planet (bottom panel of Fig. 9) end up with two planets on well separated orbits, generally after an ejection or a planet–planet collision. As the eccentricity of these orbits is significant, the stability criterion used previously (Eq. 9) is not valid. Instead, we check if the system are angular momentum deficit (AMD) stable (Laskar & Petit, 2017; Petit et al., 2017). For coplanar orbits, the AMD of a two-planet system is given by

[TABLE]

A given system is AMD stable if the orbits of the two planets cannot intersect through free exchange of AMD between the two planets (Laskar & Petit, 2017). This criterion is valid as long as the two planets are not in mean-motion resonance. For completeness, we also check if the systems are in the chaotic area due to the overlap of first-order MMRs, which is given by (Petit et al., 2017)

[TABLE]

where

[TABLE]

For our considered range of masses, this criterion (Eq. 16) is valid when both eccentricities $e_{k}\gtrsim 0.2$. We find that more than 99 % of the resulting systems with two planets are AMD stable.

The two-planet systems represented in Fig. 9 are significantly more diverse than in the two-planet scattering case (cf. Fig. 7). They generally have an inner planet with a semi-major axis comparable to or lower than the initial innermost planet, while the outer planets (filled circles) have their pericentre distributed around 15 au (grey curve). This means that, roughly, their pericentre remains near their initial semi-major axis. However, the departure from this curve can be significant.

For comparison, we re-ran the same initial conditions with strong tides (Equations 12). The result is displayed in Figure 10. The effect on the final ${e_{\rm p}}$–${a_{\rm p}}$ distribution is clearly less important than in the two-planet scattering case (cf. Figs. 7 and 8), which implies that the ejections or planet–planet collisions tend to occur before the pericentre of the innermost planet reaches a few hundredths of an astronomical unit. In fact, only approximately 1%  of the inner planets see their pericentre drop below 0.1 au throughout their orbital evolution.

In total, a significant fraction of the systems ends up with a planet on a wide orbit, with an apocentre several times higher than the initial semi-major axis. If we compare this to the projected distance of HIP 65426 b (92 au), with our choice of initial conditions $\approx 18\,\%$ of the systems ended up with a planet whose apocentre is above 90 au after $5\times 10^{6}$ years in the conservative case ($7\,\%$ of the single planets and $21\,\%$ of the two-planet systems), against $\approx 16\,\%$ when strong tides are modelled.

4.3 Conclusion about giant planet scattering

Both two-planet and three-planet scattering scenarios are able to create systems with giant planets on wide orbits, with semi-major axes above 100 au, even starting with planets in the vicinity of 10 au. In both cases, these planets tend to be highly eccentric (${e_{\rm p}}\gtrsim 0.5$, and generally more) as they retain a pericentre close to their initial semi-major axis. However, the occurrence rate of these orbits, as well as the presence and properties of an eventual giant planet companion, greatly depend on the studied scenario:

• •

In two-planet scattering, we have found that at most a small fraction of the systems (depending on the chosen tidal model) end up with a planet on a stable orbit with an apocentre significantly raised with respect to their initial semi-major axis. However, the instability between two planets may occur while at least a partial disc is remaining, which may lead to a broader range of outcomes (see Sect. 3). Although most of the systems that ended up with a planet on a wide orbit (with an apocentre significantly larger than its initial one) were single-planet systems, a proper model of the tides can circularise an inner planet on a tight orbit instead of letting it migrate all the way into the star. This would cause HIP 65426 to have another giant planet on a much shorter orbital period, possibly observable using the radial-velocity method. However, it would not be observable with current direct-imaging techniques such as Sparse Aperture Masking, which push down to a few au for this system (Cheetham et al., 2019).

• •

In the three-planet case, the outcomes are much more diverse. In $\sim 3/4$ of the cases, two planets remain in the system on stable orbits. Most of these systems have a planet with a semi-major axis significantly higher than initially and an inner planet with a semi-major axis comparable to the initial one or lower (see Figs. 9 and 10). Figure 11 shows the ${a_{\rm p}}$ and ${e_{\rm p}}$ ratios between the inner and outer planets of Fig. 9, and which planet of the pair is more massive. There is no tendency for either the inner or the outer planet to be the more massive one, and both planets tend to have comparable eccentricities (histogram in Fig. 11, bottom panel). If directly imaged planets such as HIP 65426 b obtained their wide orbit through planet–planet scattering of three giant planets initially in the 3–20 au range, it is probable that these systems also contain an eccentric inner planet with a semi-major axis of or greater than a few astronomical units. Observations are not yet constraining enough to confirm or refute the existence of such a planet around HIP 65426  as the current limits only exclude a planet more massive than $5\leavevmode\nobreak\ {M_{\textnormal{J}}}$ outside of 15 au (Chauvin et al., 2017), or a planet more massive than $16\leavevmode\nobreak\ {M_{\textnormal{J}}}$ outside of 3 au (Cheetham et al., 2019), and we recall that in the three-planet scattering case the remaining inner planet would not necessarily be more massive than HIP 65426 b (see Fig. 11).

5 Conclusion

5.1 Summary

The planet imager SPHERE (Beuzit et al., 2008, 2019) recently revealed a companion to the 2 $M_{\odot}$, $14\pm 4$ Myr Lower Centaurus-Crux group member HIP 65426. The initial analysis by Chauvin et al. (2017) showed it to be of planetary mass, with ${m_{\rm p}}=6$–12 ${M_{\textnormal{J}}}$, while located at a separation from its host star (projected: 92 au) at which formation by core accretion is not expected to be efficient. Combined with the star’s unusually high rotation rate ($v_{\star}\approx 300$ km s*-1*), this motivated us to take a closer look at the system to (i) infer joint constraints on the mass and initial (post-formation) entropy, (ii) explore the formation of wide-orbit (directly imaged) planets by core accretion, and (iii) derive predictions about the presence of further companions in the system. While we focused on HIP 65426 b, it is an excellent representative for the relatively recent and modestly populated class of directly imaged exoplanets in terms of mass, age, and separation from its host star.

First, we derived constraints on the mass and initial entropy of HIP 65426 b from its age and luminosity (Section 2). Assuming it formed by core accretion (CA), we argued that HIP 65426 b should be roughly ${\Delta t_{\rm form}}\approx 2$ Myr younger than its host. We considered different priors on the mass and entropy, including for the first time the mass and post-formation entropy distribution of the Mordasini et al. (2017) population synthesis. The simple but robust 2D fits for ${\rm d}^{2}N/({\rm d}{m_{\rm p}}\,{\rm d}s_{\textrm{pf}})$ in Eqs. (2)–(8) may be useful in other work. Flat priors yielded ${m_{\rm p}}=9.8^{+1.5}_{-2.0}\leavevmode\nobreak\ {M_{\textnormal{J}}}$, whereas the priors from the hot and cold population from the population synthesis lead to ${m_{\rm p}}=9.9^{+1.1}_{-1.8}\leavevmode\nobreak\ {M_{\textnormal{J}}}$ and ${m_{\rm p}}=10.9^{+1.4}_{-2.0}\leavevmode\nobreak\ {M_{\textnormal{J}}}$, respectively. Independently of the priors, the minimal post-formation entropy could be constrained to $s_{\textrm{pf}}\gtrsim 9.2\leavevmode\nobreak\ k_{\textrm{B}}\,\textrm{baryon}^{-1}$. Using the population synthesis priors made a large difference, providing an upper bound and yielding $s_{\textrm{pf}}=10.4^{+0.7}_{-0.2}$ in the hot-population case and $s_{\textrm{pf}}=10.2^{+0.3}_{-0.7}$ for the cold one.

Next, we studied the formation of wide-orbit gas giants by core accretion (Sect. 3). The idea is to let a core that formed in the inner disc be scattered by a companion into the outer disc, where it can undergo runaway accretion. If this scattering happens late enough, the finite lifetime of the disc combined with the slower type II migration rate should allow the planet to stay at large semi-major axes. To counter the fast type I migration while the core forms, we included, as in Coleman & Nelson (2016a), a specific radial structure which acts as a planet trap and could be due to zonal flows.

This scenario was seen to work well, producing HIP 65426 b-like planets in a number of cases (Fig. 5). In almost all systems, they were accompanied by interior giant companions that could be detectable in long-baseline radial-velocity surveys. Another possibility is of instabilities after disc dispersal. This too was shown to be a possible origin for HIP 65426-like systems, again with the prediction of further interior companions

Finally, we focused on the post-disc phase with $N$-body integrations of two- or three-planet systems including tides (Section 4). Systems with two planets usually ($\sim 2/3$ of the time) featured a planet–planet collision, with almost all other cases ending up with a planet ejection. For both outcomes, the remaining planet still retained too small a semi-major axis (${a_{\rm p}}\sim 10$ au) to explain HIP 65426 b. In the case of three planets initially, roughly half of the systems did not interact significantly within 5 Myr. Of the others, about 1/4 lost two planets, with the remaining planet matching HIP 65426 b only a small fraction of the time. Systems with two remaining planets however had more diverse configurations in the ${a_{\rm p}}$–${e_{\rm p}}$ plane. For our choice of initial conditions, $\sim 1/5$ of the systems ended up with a planet with an apocentre above 90 au (HIP 65426 b’s projected separation). We also looked at the effect of tidal circularisation, which can affect the orbit of highly eccentric planets that pass close to the star. We showed, however, that in both the two- and three-planet scenarios the outcomes are changed only slightly (Figs. 8 and 10).

5.2 Discussion

The main implications of our study are the following:

We estimate a mass of ${m_{\rm p}}=9.9^{+1.1}_{-1.8}\leavevmode\nobreak\ {M_{\textnormal{J}}}$ using the hot population and ${m_{\rm p}}=10.9^{+1.4}_{-2.0}\leavevmode\nobreak\ {M_{\textnormal{J}}}$ with the cold-nominal population for HIP 65426 b. 2. 2.

As for almost all other directly imaged planets, we find that HIP 65426 b is not consistent with the extreme cold starts à la Marley et al. (2007). This is also in agreement with recent theoretical work (Berardo et al., 2017; Marleau et al., 2017). 3. 3.

A more precise mass determination is hindered here less by systematics between the different atmospheric models (see Fig. 1) than by the large relative uncertainty on the stellar age. The uncertain formation time ${\Delta t_{\rm form}}$ is subdominant to this. 4. 4.

Both runaway accretion at a large separation after outward scattering of the core as well as post-disc-phase scattering of inner gas giants were seen as viable scenarios to explain HIP 65426 b-like objects. 5. 5.

Our simulations show that if it formed through core accretion, HIP 65426 b likely has some eccentricity. This eccentricity arises from scattering with other planets in the system. If these scattering events occur before the end of the disc lifetime, damping with the gas disc can act to reduce the eccentricity. In this case the planet would have a modest eccentricity, $0\leq{e_{\rm p}}\leq 0.5$, where the time of scattering with respect to the end of the disc lifetime determines how much eccentricity can be damped. If the scattering event took place after the end of the disc lifetime, then the eccentricity can be higher, depending on the scattering conditions. Therefore, if future observations revealed the eccentricity to be ${e_{\rm p}}<0.5$, this would not rule out the scenario of scattering before disc dispersal. It would, however, make scattering at the end of or after the disc phase unlikely, unless we could invoke another kind of eccentricity-damping mechanism.

The high-eccentricity cases are in contrast with the very tentative result that directly imaged planets might tend to have low eccentricities. However, this is mostly based on a relatively small number of upper limits (Bowler & Nielsen, 2018), and the few cases with more robust determinations are not likely candidates for the scenario presented here.

For example, several authors have favoured ${e_{\rm p}}\lesssim 0.2$ for $\beta$ Pic b (Wang et al., 2016; Lagrange et al., 2018), which might suggest it did not form by the scenario shown here. While Dupuy et al. (2019) recently excluded ${e_{\rm p}}<0.1$ at $>2\,\sigma$, their derived eccentricity was only ${e_{\rm p}}=0.24\pm 0.06$. Independently of the (modest) eccentricity of $\beta$ Pic b, however, its low semi-major axis ${a_{\rm p}}\approx 12$ au makes it a somewhat unlikely candidate for formation by scattering. Also, the presence of a debris disc makes any speculation about its dynamical origin more challenging.

As for the specific case of the HR 8799 planets, it is unlikely that a scenario such as that studied here is responsible for their formation, irrespective of their exact eccentricities (which appear to be low to moderate; Wang et al., 2018). Indeed, this would require an unlikely series of interactions without ejections, for example four times in a row.

In any case, a longer coverage of the orbits will be necessary to improve the statistics of the eccentricity determinations. 6. 6.

If directly imaged planets such as HIP 65426 b obtained their wide orbit through planet–planet scattering of three giant planets initially in the 3–20 au range, it is probable that these systems also contain an eccentric inner planet with a semi-major axis equal to or greater than a few astronomical units.

Some previous studies have put upper limits to the existence of inner companions in systems with wide-orbit planetary-mass companions (see Bryan et al., 2016 and references therein). However, these limits typically reach down to several tens of astronomical units for several Jupiter masses, and thus leave open a parameter space consistent with our tentative prediction. This could be tested by future radial-velocity surveys or further direct-imaging observations.

We can wonder whether the inferred initial entropy will reveal clues to the location of the runaway gas accretion phase. In the absence of detailed studies of this question, the answer seems complex since several effects are relevant at the same time:

• •

For a given planet mass, radius (or entropy), and accretion rate, it is easy to show that the properties of the accretion shock onto the planet will not depend, at least not directly, on its location in the disc (Marleau et al., 2017, Marleau et al., subm.).

• •

The Core Mass Effect (Mordasini, 2013) predicts higher post-formation entropies for higher planetesimal surface density, as should be found closer in to the star.

• •

Berardo et al. (2017) showed that the important quantity determining the influence of the shock is the pre-runaway entropy of the protoplanet. This in turn might be different for planets formed at different locations, but the magnitude of the effect is challenging to assess without dedicated simulations.

Further factors might come into play, such as the metallicity of the gas. Depending on which way the different effects go, the location of runaway gas accretion may or may not be imprinted in the post-formation entropy. Clearly, a global dedicated study is warranted here.

Thanks to the large separation and super-Jupiter mass of the companion, HIP 65426 represents an important system to explore the dynamical interactions of (proto)planets and the limits of planet formation by core accretion. Recent studies have reached contrasting results about gravitational instability (GI), arguing that it must be an intrinsically rare process (Forgan & Rice, 2013; Vigan et al., 2017) or rather that it is common but associated with very fast migration of the clumps (see Nayakshin, 2017 and discussion therein as well as Vorobyov & Elbakyan, 2018). However, these studies looked mainly at FGK stars, whereas the planet occurrence rate seems to increase with stellar mass (Bowler, 2016). In any case, it would be interesting also to perform a study similar to the present one in the context of GI, following self-consistently the formation of the central star (see e.g. Nixon et al., 2018) and trying to explain its high spin frequency. Also, predictions of the post-formation entropies and luminosities in GI formation models (e.g. Forgan & Rice, 2013; Forgan et al., 2015) would be a welcome counterpart to those of core accretion (Mordasini et al., 2017). With an orbital period $P\approx 400$–2000 years (Cheetham et al., 2019), approximately five to ten years are required until the eccentricity of HIP 65426 b can be robustly determined if orbital curvature begins to be resolved (G. Chauvin 2019, priv. comm.). However, this – and more generally for other systems too – will be an important constraint on the formation model. Similarly, further radial-velocity monitoring of the host to reveal or rule out the presence of further companions would also help confirm or exclude some of the formation pathways discussed here.

Acknowledgements.

We thank G. Chauvin for helpful comments about the manuscript and answers to observational questions, A.-M. Lagrange for the interesting discussions, and the referee for the comments and questions that helped clarify the manuscript. G.-D.M. and C.M. acknowledge support from the Swiss National Science Foundation under grant BSSGI0_155816 ‘PlanetsInTime’. G.-D.M. acknowledges the support of the DFG priority program SPP 1992 ‘Exploring the Diversity of Extrasolar Planets’ (KU 2849/7-1). Parts of this work have been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation.