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

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.


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. E-mail address: gabriel.marleau{@space.unibe.ch, @uni-tuebingen.de}. Current affiliation: Institut für Astronomie und Astrophysik, Eberhard Karls Universität Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany.
CHEOPS Fellow The first discovery of the SPHERE instrument at the VLT (Beuzit et al. 2008(Beuzit et al. , 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 sin i = 299 ± 9 km s −1 . 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 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 Article number, page 1 of 17 arXiv:1902.01869v2 [astro-ph.EP] 20 Feb 2019 A&A proofs: manuscript no. ms_neu evolution and N-body interactions. Then in Section 4 we use Nbody 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.

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 Section 2.4.

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 phasespace 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Ṁ ∝ m 1/6 Σ m ρ 2/5 gas , 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 Section 3 us-ing 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.

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 ∆τ 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).
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 1 . 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 Time (yr) 0.5, 1, 3, 5, 10, 16, 20 mM ⊙ Baraffe et al. (2003) Burrows et al. (1997 completo 21 + Eddington completo 21 + AMES-COND completo 21 + Burrows et al.  Baraffe et al. 2003, Burrows et al. 1997; see legend). The Bern evolution (cooling) code completo 21 is used and compares very well to the original models. Units of milli-solar masses (1 mM = 1.05 M J ) are used to reproduce as closely as possible the tracks of Burrows et al. (1997) and Baraffe et al. (2003). 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). The faint grey cross shows β Pic b (Bonnefoy et al. 2014b) as an example error bar. 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).  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 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.

First analysis of the mass
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, L hottest pf traces the brightest planet at every mass; L cold-nom. pf corresponds to the cold-nominal population, in which gas is assumed to accrete cold; L cold-class. pf 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 coldest pf 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 pf = 1.2×10 −5 L (m p /M 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 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.

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 the planet; (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., 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 Mordasini et al. (2017) showed that there is spread of postformation entropies of approximately ∆s pf ≈ 1 k B 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 pf : The following functions, dropping the usual entropy units where masses m p are implicitly in Jupiter masses in these equations, while for the hot population, the bounds are s hot pf, min = 10.00 + 0.12 (m p − 10) (5a) s hot pf, max = 11.300 + 0.116 (m p − 10).
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(Paxton et al. , 2013(Paxton et al. , 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 i.e. the distribution is nearly flat in log m 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 dN/dm p ∝ m p −1.05 fit of Marcy et al. (2005). We note, however, that Cumming et al. (2008) found dN/dm p ∝ m p −1.3±0.2 but for periods < 2000 days, while Brandt et al. (2014) obtained from direct imaging dN/dm p ∝ m p −0.7±0.6 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. Figure 3 shows the joint constraints on the mass and postformation (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 postformation entropy s pf 9.2 k B baryon −1 but it is not otherwise 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 m p = 9.0 +1.3 −1.5 M J 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 2 . As expected from Marleau & Cumming (2014), the relative uncertainty on the hot-start mass σ m p /m p ≈ 0.2 is ≈ 1 2 σ t cool /t cool ≈ 0.3. 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 m p = 9.4 +1.5 −2.0 M J , 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 m p = 10.9 +1.4 −2.0 M J for the cold population (Eq. (4)) and m p = 9.9 +1.1 −1.8 M J 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 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 pf prior despite the p m p ∝ m p −1 scaling since the mass interval is small.

Results: Mass-entropy constraints
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 s pf = 10.4 +0.7 −0.2 in the case of the hot population and s 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 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 The value m p = 7 +2 −1 M 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).

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 m p = 9.0 +1.9 −1.7 M J and with only the mass prior m p = 8.7 +1.9 −2.0 M J , whereas using the mass and s pf priors from hot-start (cold-start) populations yields m p = 9.3 +1.3 −1.8 M J (m p = 10.3 +1.6 −1.8 M J ). 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 m p = 9.9 +1.4 −1.7 M J ; instead, using the s pf and mass priors from the hot-start (cold-start) populations yields m p = 10.4 +1.0 −1.6 M J (m p = 11.3 +1.0 −1.7 M J ). 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 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 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 log L/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 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 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 m p = 10.9 +1.4 −2.0 M J with priors from the cold population and m p = 9.9 +1.1 −1.8 M J using the hot population.

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 protoplane-Article number, page 5 of 17  tary 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 semimajor axis could be similar to those of directly imaged plan-ets 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).

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 Φ 41 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 Notes. For the meaning of Φ 41 see Dullemond et al. (2007).
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)  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 α 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. (2010Paardekooper et al. ( , 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 selfconsistently 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 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. 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 eccentric- ity 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).

Example HIP 65426-like simulated system
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' semimajor 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).

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 distanceeccentricity correlation.
Also seen in Fig. 5 are numerous giant planets with semimajor 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 (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. 2016Collaboration et al. , 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.

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 planetplanet 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) their orbits will be unstable on short timescales, typically of order τ ∼ ((m 1 + m 2 )/m ) −1/3 P 1 , 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 . Being in the unstable area given by Equation (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 ).
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 (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 3 .
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 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 = 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 R k = (1.1 + 0.06 × m k /M J ) R J , which roughly fits the non-accreting hot population of Mordasini et al. (2012a) at 3-5 Myr 4 . Planets that entered the Roche limit of the star No tides Fig. 7. Outcome of the scattering of two giant planets initially on circular orbits with semi-major axes near 10 au in the conservative (i.e. non-dissipative) case. The few systems that kept their two planets over 5 × 10 6 years are displayed in grey ( 1 % of the systems). Other systems evolved into one-planet systems typically after 10 5 years, either through planet-planet collision (≈ 63 %; red), ejection of the other planet (≈ 36 %; blue), or collision of the other planet with the star (≈ 1 %; orange). The yellow line is the predicted orbit of the remaining planet after an ejection for typical values (see Eq. 9). To the right of the black line are planets whose apocentre is above 90 au (projected distance of HIP 65426 b) and the grey line shows the orbits whose pericentres are at 15 au, which is the detection limit of an eventual companion of mass m p 5 M J to HIP 65426 b (Chauvin et al. 2017).
were removed from the simulation, with the Roche limit given by R ,Roche ≈ 2.2 R (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 semimajor 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.

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 a 2 = a 1 (1 + 1.42((m 1 + m 2 )/m ) 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 r lies between the semi-major axes of the initial ones : This is the most common outcome (≈ 63 % of all systems), leading to the red clump between 10 and 15 au in Fig. 7. These semimajor 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 1 + m 2 = 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 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.

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 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 cir-cularisation, 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 L k = m k (Gm k R k ) 1/2 and E k = Gm 2 k /R k are the angular momentum and energy scales, ζ k = (m k q 3 k ) 1/2 (m R 3 k ) −1/2 , ω 0 is the dimensionless frequency of the fundamental mode (normalised by the internal dynamical frequency [Gm k p/R 3 k ] 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: 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 γ k = (1 − e 2 k )/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 (11), and for strong tides, using the set of equations (10). The results are presented in the top and bottom panel of Fig. 8, respectively. We assumed that ω 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 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 semimajor 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.

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 ≈ 5-7 au). This initial spacing does not necessarily ensure instability in the system within the 5 × 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×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. 7) is not valid. Instead, we check if the system are angular momentum deficit (AMD) stable Petit et al. 2017). For coplanar orbits, the AMD of a two-planet  Fig. 9. Outcome of the scattering of three giant planets initially on circular orbits with semi-major axes near 10 au in the conservative (tidefree) case. Only the systems that lost at least one planet (≈ 46 % of the initial 1000) are represented. The black and grey lines are as in Figs. 7 and 8. Top panel: Systems that ended up with a single planet at the end of the run (≈ 23% of the systems that underwent a strong instability). Blue dots represent planets whose two companions were ejected, purple for the systems that underwent both planet-planet collision and ejection, orange for those that underwent both ejection and collision with the star, and red when two planet-planet collisions occurred. Bottom panel: Systems that ended up with two planets at the end of the run (≈ 77% of the systems that underwent a strong instability). The colourcoding is the same as in Figs. 7 and 8, with open circles for the inner planet and filled circles for the outer one. system is given by A given system is AMD stable if the orbits of the two planets cannot intersect through free exchange of AMD between the two planets . 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. where 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 (Equations 10). The result is displayed in Figure 10. The effect on the final e p -a p distribution is clearly less important Article number, page 13 of 17  Fig. 11. Comparison of the orbits of the inner and outer planets of each system represented in Figure 9. In the top panel, the colour indicates which planet is the more massive: the inner one (red dots) or the outer one (blue dots). Axes show the ratio of the semi-major axes and of the eccentricity, respectively. The bottom panel displays the distribution of e out /e int . We note the logarithmic horizontal axis.
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 6 years in the conservative case (7 % of the single planets and 21 % of the two-planet systems), against ≈ 16 % when strong tides are modelled.

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 semimajor 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).

Summary
The planet imager SPHERE (Beuzit et al. 2008(Beuzit et al. , 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 −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 ∆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 2 N/(dm p ds pf ) in Eqs.
(2)-(6) may be useful in other work. Flat priors yielded m p = 9.8 +1.5 −2.0 M J , whereas the priors from the hot and cold population from the population synthesis lead to m p = 9.9 +1.1 −1.8 M J and m p = 10.9 +1.4 −2.0 M J , respectively. Independently of the priors, the minimal post-formation entropy could be constrained to s pf 9.2 k B baryon −1 . Using the population synthesis priors made a large difference, providing an upper bound and yielding s pf = 10.4 +0.7 −0.2 in the hot-population case and s 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 blike 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 (∼ 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).

Discussion
The main implications of our study are the following: 1. We estimate a mass of m p = 9.9 +1.1 −1.8 M J using the hot population and m p = 10.9 +1.4 −2.0 M J 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. 2018), 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 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 semimajor 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 planetarymass 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: