© G.-D. Marleau et al. 2019

Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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_p = 8–12 M_J dusty L6 ± 1 companion to the m_⋆ = 1.96 ± 0.04 M_⊙ fast rotator HIP 65426, which has an equatorial velocity v_⋆ sini = 299 ± 9 km s⁻¹. Its projected separation is 92.0 ± 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 has a 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 Sect. 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 Sect. 3 we use detailed planet formation models following the disc evolution and N-body interactions. Then in Sect. 4 we use N-body integrations to look in detail at interactions between several companions once the disc has cleared. Finally, in Sect. 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 κ And b and β 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 Sect. 2.4.

2.1 Luminosity and age

Firstly, we discuss the input quantities for the MCMC. The adopted bolometric luminosity is log L∕L_⊙ = −4.06 ± 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 ± 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 ± 4 Myr. We note that at 2 M_⊙, 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 Δt_form (Fortney et al. 2005; Bonnefoy et al. 2014b), which we now briefly discuss.

While the dependence of the formation time Δt_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 , where Σ_m and ρ_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 Σ_m and ρ_gas are expected to increase with stellar mass, the formation time should decrease with planet mass. Also, observationally, the formation time Δt_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 Δt_form ~ 3–5 Myr given the lifetimes of protoplanetary discs (Haisch et al. 2001). Finally, population synthesis calculations for a 2 M_⊙ central star (Mordasini et al., in prep.) indicate that most ~10 M_J planets (approximately the mass of HIP 65426 b, as we show later) have reached their final mass after roughly Δt_form ~ 2 Myr, and the simulations presented in Sect. 3 using the Coleman & Nelson (2016a) models for a 2 M_⊙ star yield Δt_form ≈ 2.5–4 Myr. Therefore, we adopt Δt_form = 2 Myr, and thus t_cool = 12 ± 4 Myr as the fiducial age. We note that this Δt_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_cool = 14 ± 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. 2019)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 Δτ_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 afew percentage points. This coupling approach was applied recently to low-mass planets in Linder et al. (2019).

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

• (i)

  the Eddington assumption;

• (ii)

  AMES-COND (Baraffe et al. 2003);

• (iii)

  Burrows et al. (1997);

• (iv)

  petitCODE (Mollière et al. 2015);

• (v)

  HELIOS (Malik et al. 2017, 2019).

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) grids¹. 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_⊙ and 2 M_⊙ population syntheses, objects need a mass of 16 M_J (20 M_J) to have consumed even only ≈ 30% (≈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 mM_⊙ 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_p ≈ 8–11 M_J, which is notrare 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 basedon 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_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_p ≈ 0.3 to ≈ 12 M_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

respectively, for the hottest, cold-nominal, cold-classical, and coldest planets. Briefly, traces the brightest planet at every mass; corresponds to the cold-nominal population, in which gas is assumed to accrete cold; 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, 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 ), 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_KH ~ 100 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_KH ~ 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:

• (i)

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

• (ii)

  a hot population, which differs from the first only by the assumption that the entire accretion luminosity is brought into theplanet;

• (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, and 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 nowturn to the total distribution function, which we write as (2)

Mordasini et al. (2017) showed that there is spread of post-formation entropies of approximately Δs_pf ≈ 1 k_B baryon⁻¹ 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_pf: (3)

The following functions, dropping the usual entropy units k_B baryon⁻¹, closely fit the envelope of points s_pf(m_p) in Mordasini et al. (2017). For the cold-nominal population, the lower and upper edges are given respectively by

where masses m_p are implicitly in Jupiter masses in these equations, while for the hot population, the bounds are

This holds down to 2 M_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_J, approximately given by (6)

i.e. the distribution is nearly flat in logm_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 fit of Marcy et al. (2005). We note, however, that Cumming et al. (2008) found but for periods < 2000 days, while Brandt et al. (2014) obtained from direct imaging at distances ≳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_pf ≥ 9.2 k_B baryon⁻¹ but it is nototherwise constrained. This lower limit holds independently of the formation pathway and for masses up to m_p ≈ 15 M_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_p = 9.6 ± 1.7 M_J. For high values of s_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_pf ≳ 14 to approximate what is usually thought of as hot starts, we find for a cooling age t_cool = 12 ± 4 Myr. This agrees well with the m_p = 10 ± 2 M_J reported by Chauvin et al. (2017) for the DUSTY models². As expected from Marleau & Cumming (2014), the relative uncertainty on the hot-start mass is .

Next we fold in the outcome of the population syntheses into the analysis. If we take only the mass prior (Eq. (6)) into account, we obtain , which is lower by Δm_p ≈ 0.2 M_J than the case without priors. Applying the s_pf priors (Eq. (3)) as well, we obtain for the cold population (Eq. (4)) and for the hot population (Eq. (5)).

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 Δm_p ≲ 0.2 M_J. These differences represent only a modest fraction of the error bars. However, the s_pf priors are mildly important, leading to a difference Δm_p ≈ 1 M_J between the hot and the cold populations and even Δm_p = 1.5 M_J between theflat-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_pf prior despite the scaling sincethe 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_pf ≳ 9.2 obtained with the uniform prior does not change, but the population-synthesis priors lead to the determination of an upper bound, yielding in the case of the hot population and 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 aconsequence of the core-mass effect (CME) as discussed by Mordasini et al. (2017, Sect. 5.2.1). These values s_pf ≈ 10.3 are clearly lower than classical (arbitrarily) hot starts (s_pf ≈ 13), with an initial Kelvin–Helmholtz time t_KH ~ 10 Myr as opposed to t_KH ≲ 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_pf posteriors.

2.6 Discussion

For comparison, with a shorter cooling age t_cool = 10 ± 4 Myr (i.e. coming from a longer formation period), we obtain with uniform priors and with only the mass prior , whereas using the mass and s_pf priors from hot-start (cold-start) populations yields (). Instead, taking a cooling age t_cool = 14 ± 4 Myr, i.e. the age of HIP 65426 as might correspond to formation by gravitational instability, with only mass priors we obtain ; instead, using the s_pf and mass priors from the hot-start (cold-start) populations yields (). This is somewhat higher than, but still consistent with, the mass derived by Cheetham et al. (2019). Using an age of 14 ± 4 Myr and the AMES-COND models, they found m_p = 7.5 ± 0.9 M_J based on magnitudes in individual bands and m_p = 8.3 ± 0.9 M_J based on their bolometric luminosity, which had an uncertainty σ_{log L} = 0.03 dex half as largeas the value used here. 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_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 d logL∕ds 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_pf, the posterior was also relatively flat, and a small distortion will not change the nature of the weak constraints on s_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 willnot 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 Sect. 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_p ≈ 13 ± 2 M_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 with priors from the cold population and 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_p ≲ 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 α and the photoevaporation factor Φ₄₁ to give the disc an appropriate lifetime. Using the values presented in Table 1, the initial disc had a total mass equivalent to ~ 8% of the mass of HIP 65426 (i.e. around 150 M_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 Sect. 2.3.3). This radial structure increases the local α 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 M_⊕≤ m_p ≤ 20 M_⊕ at a_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 2016a,b). 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 coreswhich 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 ~ 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_p ≈ 220 au with an eccentricity e_p ≈ 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_p ≈ 2 M_J, it has a close encounter (label D in Fig. 4) with the most massive giant in the system (with a mass m_p ≈ 4 M_J, shown by the blue line). This lowers the semi-major axis and eccentricity of the former to a_p ≈ 125 au and e_p ≈ 0.7 respectively (see the drops in semi-major axis and eccentricity at ~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_p ≈ 9.8 M_J and has migrated in to having a semi-major axis of a_p ≈ 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_p ≈ 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_p = 9.4, 15.6, and 25.8 au, and masses m_p = 0.45, 2.8, and 7.8 M_J, respectively, on nearly circular orbits (e_p ≲ 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_p > 50 au and mass m_p > 1 M_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_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_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 (~ 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_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 (Sect. 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 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σ the presence of further companions more massive than 16 M_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_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_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 ~ 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) (7)

their orbits will be unstable on short timescales, typically of order (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_⋆ 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 Eq. (7) after disc dispersal requires one of two scenarios:

• (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;

• (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 (Sect. 3).

On the other hand, a system with three or more planets does not have such a sharp stability boundary as in Eq. (7). 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 M_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 planets³.

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 × 10⁶ yr, which is comparable to the age of the system since disc dispersal (see Sect. 2.1). Alternatively, integrations were stopped when only one planet remained in the system. The mass of the star was set to m_⋆ = 2 M_⊙, while the mass of each planet was randomly picked in the interval m_k = 5–15 M_J. The radius R_k of each planet k was set to , which roughly fits the non-accreting hot population of Mordasini et al. (2012a) at 3–5 Myr⁴. Planets that entered the Roche limit of the star were removed from the simulation, with the Roche limit given by R_⋆,Roche ≈ 2.2 R_⊙ (for comparison, the stellar radius is R_⋆ = 1.77 R_⊙; 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_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 ~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°]. 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_p = 10 au, while the outer one was positioned slightly inside the instability domain (Eq. (7)), at . We integrated 500 systems with this set of initial conditions; the final outcomes are plotted in Fig. 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⁵ yr, 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_r lies between the semi-major axes of the initial ones (Ford & Rasio 2008): (8)

This is the most common outcome (≈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

where the subscripts “r” and “ej” refer to the remaining and ejected planet, respectively, and where q_ej is the ejected planet’s minimal distance to the star on its parabolic orbit. This scenario represents almost all other cases (≈ 36% of the systems). The range of possible orbits for the typical values m₁ + m₂ = 20 M_J and q_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 (≈ 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 mostof 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 ≈ 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) ~ R_⋆ ~ 0.01 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⁵ yr (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 (ΔL_tide) and energy (ΔE_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

where and are the angular momentum and energy scales, , ω₀ is the dimensionless frequency of the fundamental mode (normalised by the internal dynamical frequency ), 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:

We translate either of these two expressions for ΔL_tide and ΔE_tide into migration and eccentricity damping timescales using respectively (Nagasawa et al. 2008):

with 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 Eq. (11), and for strong tides, using the set of Eq. (10). The results are presented in the top and bottom panel of Fig. 8, respectively. We assumed that ω₀ = 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 M_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₂ at 10 au, and m₁ and m₃ at four mutual Hill radii inside and outside the orbit of m₂, respectively (this corresponds to spacings ≈5–7 au). This initial spacing does not necessarily ensure instability in the system within the 5 × 10⁶ yr 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). However, 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 × 10⁶ yr. 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.

Systemsthat 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. (7)) 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 (13)

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) (14)

where (15)

For our considered range of masses, this criterion (Eq. (14)) is valid when both eccentricities e_k ≳ 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 (Eq. (10)). The result is displayed in Fig. 10. The effect on the final e_p –a_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 ≈ 18% of the systems ended up with a planet whose apocentre is above 90 au after 5 × 10⁶ yr in the conservative case (7% of the single planets and 21% of the two-planet systems), against ≈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_p ≳ 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 ~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_p and e_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 M_J outside of 15 au (Chauvin et al. 2017), or a planet more massive than 16 M_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_⊙, 14 ± 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_p = 6–12 M_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_⋆ ≈ 300 km s⁻¹), 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 (Sect. 2). Assuming it formed by core accretion (CA), we argued that HIP 65426 b should be roughly Δt_form ≈ 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 d ² N∕(dm_p ds_pf) in Eqs. (2)–(6) may be useful in other work. Flat priors yielded , whereas thepriors from the hot and cold population from the population synthesis lead to and , respectively. Independently of the priors, the minimal post-formation entropy could be constrained to s_pf ≳ 9.2 k_B baryon⁻¹. Using the population synthesis priors made a large difference, providing an upper bound and yielding in the hot-population case and 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 (Sect. 4). Systems with two planets usually (~ 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_p ~ 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_p –e_p plane. For our choice of initial conditions, ~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:

• 1.

  We estimate a mass of using the hot population and with the cold-nominal population for HIP 65426 b.

• 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.

  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 Δt_form is subdominant to this.

• 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.

  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 ≤ e_p ≤ 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_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_p ≲ 0.2 for β Pic b (Wang et al. 2016; Lagrange et al. 2019), which might suggest it did not form by the scenario shown here. While Dupuy et al. (2019) recently excluded e_p < 0.1 at > 2σ, their derived eccentricity was only e_p = 0.24 ± 0.06. Independently of the (modest) eccentricity of β Pic b, however, its low semi-major axis a_p ≈ 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.

  If directlyimaged 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, and in prep.).

• 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 ≈ 400–2000 yr (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 dedicate this work to the memory of Bruno Ramsauer. 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.

References

---

All Tables

All Figures