AB Aurigae: Possible evidence of planet formation through the gravitational instability

Recent observations of the protoplanetary disc surrounding AB Aurigae have revealed the possible presence of two giant planets in the process of forming. The young measured age of $1-4$Myr for this system allows us to place strict time constraints on the formation histories of the observed planets. Hence we may be able to make a crucial distinction between formation through core accretion (CA) or the gravitational instability (GI), as CA formation timescales are typically Myrs whilst formation through GI will occur within the first $\approx10^4-10^5$yrs of disc evolution. We focus our analysis on the $4-13$M$_{\rm Jup}$ planet observed at $R\approx30$AU. We find CA formation timescales for such a massive planet typically exceed the system's age. The planet's high mass and wide orbit may instead be indicative of formation through GI. We use smoothed particle hydrodynamic simulations to determine the system's critical disc mass for fragmentation, finding $M_{\rm d,crit}=0.3$M$_{\odot}$. Viscous evolution models of the disc's mass history indicate that it was likely massive enough to exceed $M_{\rm d,crit}$ in the recent past, thus it is possible that a young AB Aurigae disc may have fragmented to form multiple giant gaseous protoplanets. Calculations of the Jeans mass in an AB Aurigae-like disc find that fragments may initially form with masses $1.6-13.3$M$_{\rm Jup}$, consistent with the planets which have been observed. We therefore propose that the inferred planets in the disc surrounding AB Aurigae may be evidence of planet formation through GI.


INTRODUCTION
Most of the known exoplanets are believed to have formed in discs of gas and dust around young stars. Owing to recent advances in high resolution infrared (IR) imaging we are now capable of observing the planet formation process taking place. Observations of these discs have revealed substructures indicative of the presence of planetary companions, such as rings (ALMA Partnership et al. 2015;Andrews et al. 2016;Avenhaus et al. 2018;Bertrang et al. 2018;Dipierro et al. 2018;Huang et al. 2018a), gaps (Andrews et al. 2011;Perez et al. 2015;Ginski et al. 2016;van Boekel et al. 2017) and spirals (Garufi et al. 2013;Grady et al. 2013;Benisty et al. 2015;Pérez et al. 2016;Tang et al. 2017;Huang et al. 2018b;Dong et al. 2018), and recently it has even become possible to directly image giant protoplanets forming (Keppler et al. 2018;Müller et al. 2018;Hashimoto et al. 2011;Tang et al. 2012Tang et al. , 2017, bright inner spirals (Piétu et al. 2005), extended CO spirals (Tang et al. 2012), and the possible presence of multiple, planetary-mass companions (Piétu et al. 2005;Tang et al. 2012Tang et al. , 2017. Recent high resolution, scattered light observations of AB Aurigae performed by Boccaletti et al. (2020) using SPHERE provide some of the most spectacular images of a protoplanetary disc to date, revealing detailed spiral features, and placing new constraints on the properties of any potential companions. A kink in the inner spiral at ≈ 30 AU is found to be consistent with the presence of a protoplanet with mass of 4 − 13 M Jup (hereafter referred to as planet P1), which is also consistent with conclusions from previous authors (Piétu et al. 2005;Tang et al. 2012Tang et al. , 2017. The authors also report a point-source located at the outer edge of the inner disc, which is characterised by a gas and dust cavity at ≈ 140 AU, for which they tentatively derive a planetary mass of 3 M Jup (hereafter referred to as planet P2). Throughout this paper we aim to explore the likely formation history of planet P1.
In the core accretion (CA) model of giant planet formation (Mizuno 1980;Pollack et al. 1996), growth proceeds through the steady collisional accumulation of planetesimals onto a rocky core, which may eventually become massive enough for the onset of accretion of a gaseous envelope. Currently this model provides the most popular explanation for the formation of giant planets. However, it has been shown that formation timescales, which may be anywhere up to 10 Myr, may exceed typical disc lifetimes (Haisch et al. 2001), specifically in the case of giant planets on wide orbits where the planetesimal surface densities will be low. The discovery of systems such as HR 8799, where four ultra wide-orbit (15 AU < < 70 AU), super-Jupiter mass ( P > 5 M Jup ) planets have been directly imaged (Marois et al. 2008(Marois et al. , 2010, is an example of a particularly challenging system to explain through in-situ CA (Nero & Bjorkman 2009;Kratter et al. 2010).
CA also faces challenges when establishing how the first solids are able to grow up to and beyond metre sizes, as it is anticipated that grains will encounter growth barriers, such as the fragmentation (Birnstiel et al. 2012), bouncing (Zsom et al. 2010) and radial drift (Weidenschilling 1977) barriers. Mechanisms such as the streaming instability (Youdin & Goodman 2005; may be capable of generating local regions of extremely high particle densities, which may then undergo gravitational collapse to form the first 100 − 1000 km planetesimals . Dust trapping in the spiral regions of self-gravitating discs may also provide suitable conditions for accelerated growth (Rice et al. 2004), and possible fragmentation of the planetesimal disc (Rice et al. 2006). Multi wavelength infrared (IR) observations of discs may allow us to probe their grain size distributions (Draine 2006;Williams & Cieza 2011;Dutrey et al. 2014;Testi et al. 2014;Ilee et al. 2020), place constraints on the rate of grain growth, and investigate whether significant growth may occur very early in the disc's evolution when it is massive enough to be self-gravitating (Dipierro et al. 2015;Cadman et al. 2020a).
In the gravitational instability (GI) model of planet formation (Boss 1997), unstable regions of the disc may directly collapse to rapidly form giant gaseous protoplanets and brown dwarfs. In a differentially rotating disc, susceptibility to GI can be determined by considering the Toomre parameter (Toomre 1964), where s is the local sound speed, Ω is the orbital frequency in a rotationally supported disc, G is the gravitational constant and Σ is the local surface density. A disc may become susceptible to GI, the growth of spiral substructure, and potentially disc fragmentation, when 1.5 − 1.7 (Durisen et al. 2007).
Fragmentation of the disc will occur if unstable regions are able to cool at a faster rate than the thermal energy is generated during collapse (i.e. if the cooling rate is greater than the heating rate). If the disc is able to cool rapidly, the instability will continue to grow until the inevitable outcome of fragmentation ensues. This requirement for rapid disc cooling, which can be characterised by a critical cooling rate (Gammie 2001;Rice et al. 2005), demands that fragmentation will occur at large radii where the disc material is less optically thick, hence can cool more efficiently (Clarke 2009;Rice & Armitage 2009;Hall et al. 2017).
Calculation of the Jeans mass in a gravitationally unstable disc can be used to estimate the likely initial masses of fragments formed in this way. Using analytic approximations it has been shown that, with some dependence on the level of disc irradiation, fragmentation may initially form objects with masses between a few and a few 10s of Jupiter masses (Forgan & Rice 2011, 2013aCadman et al. 2020b). Dynamical evolution, migration, tidal stripping and growth will then follow, during which the fragment may contract to form a compact planetary/brown dwarf mass object, or be entirely torn apart and destroyed (Nayakshin 2010a(Nayakshin ,b, 2011Forgan & Rice 2013b;Nayakshin & Fletcher 2015;Forgan et al. 2018;Humphries et al. 2019).
It has also been shown that discs around higher mass stars ( * ≥ 2 M ) may be more susceptible to GI (Cadman et al. 2020b;Haworth et al. 2020), which is consistent with observations that show a higher occurrence rate of giant planets and brown dwarfs orbiting these systems (Johnson et al. 2007;Bowler et al. 2010;Nielsen et al. 2019). This suggested existence of two distinct populations of exoplanets is indicative of two modes of planet formation.
Although disc instability may not be a viable mechanism for directly forming many of the known exoplanets, it may play a role in the early growth of planet building material (Rice et al. 2004), as multi-fluid simulations of self-gravitating discs have shown significant enhancement of dust abundance present in spiral arms . It is possible, but still disputed, that it could also lead to the direct formation of the wide-orbit objects that are found via direct imaging (Vigan et al. 2017). We find ourselves in a unique position with AB Aurigae, as most of the exoplanets discovered to date have already undergone significant migration and dynamical evolution since their formation. The young age of AB Aurigae places strict time constraints on the possible formation histories of the observed planets, thus it is an ideal site for testing theories of planet formation.
In this paper we focus on the formation history of planet P1 (CA vs. GI), and whether it is possible that the AB Aurigae system could be evidence of planet formation through GI. This paper is organised as follows. In Section 2 we calculate the likely CA formation timescale of planet P1, and in Section 3 we evaluate the possibility that the planet may have formed directly through GI during AB Aurigae's early evolution. We determine the critical disc-to-star mass ratio for fragmentation in Section 3.1, and use viscous evolution models in Section 3.2 to predict whether the disc may have ever been massive enough to fragment at some point in the recent past. We place new constraints on the current mass of the disc in Section 3.2.2, and in Section 3.3 we calculate the Jeans mass in a gravitationally unstable, AB Aurigae-like disc. We discuss our results and draw conclusions in Sections 4 and 5 respectively.

Methods
To model the formation timescale of a gas giant planet through CA, we use a similar approach to that outlined in Ida & Lin (2004). We begin by assuming that either a core,init = 0.01 M ⊕ or a core,init = 0.1 M ⊕ core, with density core = 3.2 gcm −2 , has formed at a semi-major axis, , which we vary between 5 AU and 50 AU. For simplicity, we consider planet growth in-situ and neglect any migration through the disc, the effect of which is discussed in Section 4.
Core growth proceeds at a rate (Safronov 1969), where is the radius of the core, Σ p is the local planetesimal surface density, Ω is the angular frequency and is the gravitational enhancement factor, calculated using the equations from Greenzweig & Lissauer (1992). The local planetesimal surface density, Σ p , is defined as the surface density of dust within a radial annulus defined by the protoplanet's Hill radius, H , where, where * is the mass of the host star and p is the total planet mass, equal to the sum of the core and envelope masses. Whilst the core mass is still low, growth initially proceeds through planetesimal accretion, and we update p using Equation 2 at each timestep. A planet may begin to retain a gaseous envelope if the core exceeds the critical mass for the onset of gas accretion, crit , where (Ikoma et al. 2000), where is the planetesimal opacity, for which we use = 1 gcm −2 . We use a simple approach to calculate the accretion rate of a gaseous envelope onto the core, gas , based on the Kelvin-Helmholtz cooling timescale, KH , of the protoplanet, where, KH = 10 9 ( p /M ⊕ ) −3 years, and, This approximation is only valid provided that there is sufficient disc gas present for the planet to accrete, and envelope accretion will cease if the planet is able to deplete all the gas available within its feeding zone. This can be defined in terms of an upper mass limit for in-situ formation, known as the gas isolation mass, g,iso , where, g,iso = 50 where Σ g is the local gas density. We prevent further growth once p ≥ g,iso .
We set up the gas component of the disc with a total mass gas = 0.6 M , hence a disc-to-star mass ratio, = 0.25, and with Σ g ∝ −1 . The surface density profile of the gas disc is evolved using the one dimensional model outlined in Rice where we assume a radially constant, fixed value for the Shakura-Sunyaev viscous-of = 10 −3 (Shakura & Sunyaev 1973). The planetesimal component of the disc is set up as, where dust is a scale factor such that we set Σ p at 5 AU to be 2 gcm −2 , 3 gcm −2 , 5 gcm −2 and 10 gcm −2 . ice is a constant where, and ice is the ice line located at, In each case, we allow the planets to evolve in the disc for a maximum of 10 Myr. Figure 1 illustrates the resultant planet growth tracks using this formalism, considering the setups with core,init = 0.01 M ⊕ . Planet formation begins with a phase of core growth, which may either be slow or rapid depending on the local planetesimal surface density. This phase tends to plateau once the local planetesimal surface density is depleted, at which point the planet mass remains approximately constant. The critical core mass for the onset of gas accretion is proportional to the planetesimal accretion rate onto the core, and as the heating from accretion ceases the contraction of a gas envelope may ensue. Wide-orbit, giant planet formation is generally favoured near to, and just beyond the ice line due to the enhancement in the local planetesimal surface density. We calculate ice ≈ 15.6 AU for a star of mass 2.4 M . If the local planetesimal surface density is particularly high, for example near to the ice line in Figure 1c, the core mass may pass straight through the critical mass without plateauing. If the local planetesimal surface density is low, for example at a large semi-major axis in Figure 1d, the core may never experience significant growth.

Results
In Table 1 we show the results of 32 runs, where we measure the time for the core to have accreted a significant envelope, of mass equal to the core mass ( p > 2 core ), and to reach a total planetary mass of P1 = 4 M Jup , equal to the lower limit of the estimated mass for planet P1. We find it challenging to produce a planet of at least 4 M Jup in an AB Aurigae-like disc in 1 − 4 Myr. In the majority of setups considered here the planet will either reach its isolation mass before reaching the mass of planet P1, as seen in Figure 1a, or will not grow rapidly enough to reach P1 within the duration modelled here. To rapidly form a planet this massive generally requires a significant core has initially formed, in a disc with an extremely high planetesimal surface density, with a planet on a shorter orbit than where planet P1 is currently located.
The planetesimal surface densities considered here, with Σ P,5AU = 2 gcm −2 , 3 gcm −2 , 5 gcm −2 and 10 gcm −2 correspond to total planetesimal masses across the disc of 0.012 M , 0.024 M , 0.048 M and 0.072 M . These are equivalent to initial dust-to-gas ratios of 0.01, 0.02, 0.04 and 0.06 if we assume that the dust and gas in the disc decoupled when the gas mass was 1.2 M (therefore when = 0.5), and that all of this dust then went on to form planetesimals. However, if any of the dust was depleted by some other mechanism or did not go on to form planetesimals, which would likely be the case, then the planetesimal surface densities used here would demand much higher initial dust-to-gas ratio prior (1) Initial core mass.
(3) Planetesimal surface density at 5 AU. (4) Time before the planet reaches runaway growth, where the envelope mass exceeds the core mass. (5) Time before the planet mass reaches 4 M Jup to decoupling. Therefore given the generously high planetesimal surface densities used here, we would consider these core accretion timescales as optimistic lower limits to what we might expect in a realistic disc. Mechanisms which may be capable of speeding up these core accretion timescales and have not been included in the models here, such as pebble accretion, are discussed in Section 4.2.

THE GRAVITATIONAL INSTABILITY
Regions of a protoplanetary disc which become sufficiently gravitationally unstable may undergo a period of rapid collapse to directly form giant gaseous protoplanets and brown dwarfs. GI potentially offers an alternative formation mechanism for wide-orbit, giant planets whose formation timescales are difficult to explain in the CA paradigm. Currently it remains unclear whether GI may be a viable planet formation mechanism; it is uncertain whether discs may ever become sufficiently unstable to fragment, and if they are able to fragment it may be that only stellar and brown dwarf-mass companions are capable of forming in this way.

Methods
We use the P (Price et al. 2018) smoothed particle hydrodynamics (SPH) code to determine the critical mass limit for fragmentation in a disc around a 2.4 M star, analogous to AB Aurigae. SPH allows us to model the detailed hydrodynamics of a fluid, represented as pseudo-particles, each with an assigned mass, position, velocity and internal energy. A continuous fluid is approximated by calculating particle interactions through a Gaussian kernel function, with a characteristic smoothing length.
We represent the disc with = 1 × 10 6 SPH particles, distributed between in = 2.5 AU and out = 400 AU with a surface density profile Σ ∝ −1 and sound speed profile s ∝ −0.25 . We modify P such that we model radiative cooling using the hybrid radiative transfer method outlined in Forgan et al. (2009), which combines the polytropic cooling formalism from Stamatellos et al. (2007) and the flux-limited diffusion method (Bodenheimer et al. 1990;Cleary & Monaghan 1999;Mayer et al. 2007). In using a combination of these two cooling methods we model both the overall energy loss from the system and the detailed energy exchange between neighbouring particles at both high and low optical depths. We assume that disc irradiation leads to a constant background temperature, which we represent as, irr = 10 K. Artificial disc viscosity is modelled using the standard − viscosity prescription, where we use SPH = 0.1 and SPH = 0.2.
Each disc is allowed to evolve for a maximum of = 15, 550 yrs, equal to 3 orbital periods at out = 400 AU, or until fragments form and the computational timestep becomes prohibitively long for the simulations to continue. We calculate the thermalisation timescale, therm,i , from Forgan et al. (2009) for each of our disc final states, which represents the time for the disc material to reach thermal equilibrium. We find that in the discs that don't fragment within 15, 550 yrs, max( ttherm,i ) 1 kyr. It is therefore reasonable to assume that if these discs have not fragmented after 3 orbital periods, they will not do so in future.

Results
Final states of these SPH simulations are shown in Figure 2, where we vary the initial disc mass between 0.2 M − 0.35 M , which correspond to disc-to-star mass ratios = 0.08 − 0.15.
From the final states of these disc models, we expect a disc similar to AB Aurigae to fragment and form multiple clumps if disc ≥ d,crit = 0.3 M ( ≥ 0.125), and to display nonaxisymmetric substructure if disc ≥ 0.25 M ( ≥ 0.1). For disc ≤ 0.2 M ( ≤ 0.08), it is unlikely that the gravitational instability will lead to the growth of significant spirals and, in the absence of a perturber, it should be almost entirely axisymmetric. Therefore given the current low mass state of the AB Aurigae disc we predict that it should be gravitationally stable, as expected.
When also considering a set of discs with outer radii out = 300 AU and out = 500 AU we find that this critical disc-to-star mass ratio has some dependence on disc size, with more extended discs being more stable. When out = 500 AU we find the threshold for fragmentation at d,crit = 0.35 M ( crit = 0.15), and when out = 300 AU we find d,crit = 0.3 M ( crit = 0.125).  Figure 1. Evolution of planet core masses (solid line) and core + envelope masses (dashed lines) for in-situ CA planet formation at radii = 5 AU (1a), 10 AU (1b), 20 AU (1c) and 30 AU (1d) from the stellar host. In each case the models begin with an initial core mass core,init = 0.01 M ⊕ at = 0. We vary the planetesimal surface densities in the disc such that Σ P,5AU = 2 gcm −2 , 3 gcm −2 , 5 gcm −2 and 10 gcm −2 , which correspond to total planetesimal masses across the disc of 0.012 M , 0.024 M , 0.048 M and 0.072 M .

Subsequent migration of the clumps
Fragmentation will only occur if the disc is able to radiate energy away at a rate faster than the clump will collapse, hence primarily operates at large radii from the central star where the disc opacity is low thus it can cool efficiently. In the disc with = 0.125, the fragment forms at ≈ 200 AU, much further out than the current semi-major axis of planet P1. 2D hydrodynamical simulations indicate that once fragments form in a gravitationally unstable disc they will rapidly migrate to the inner regions within a few orbital periods (Baruteau et al. 2011). Computation times become prohibitively long for us to model the long-term migration of clumps in these simulations, as to resolve the high densities at the clump centres requires long integration times. Instead, typical migration of protoplanets can be approximated using the analytic calculations from Nayakshin (2010a). For type I migration, the time to move from radii out to in will be, where, and for type II migration, where, where is the disc scale height at = .
Whether a planet is in the type I or type II regime can be established in terms of a transition mass, , which roughly corresponds to the mass at which protoplanets become capable of gap-opening. For ≤ (lower-mass, faster migrating protoplanets) the planet will be in the type I regime, and for ≥ (higher-mass, slower migrating protoplanets) the planet will be in the type II regime, where, We can calculate the time for planet P1 to migrate from out = 200 AU to in = 30 AU, by substituting * = 2.4 M , p = 4 M Jup , = 0.06 for a saturated disc, and calculating the azimuthally averaged disc scale height, taken from the SPH disc where disc = 0.3 M . Integrating Equations 11 and 13 we calculate Δ mig,I = 6.9 kyr and Δ mig,II = 1.0 Myr. Note that the value of used here should be considered an upper limit as will decrease as the planet migrates. Thus the calculated mig,II would be a lower limit.
From Equation 15 we calculate the transition mass for gap opening to be = 2.4 M Jup , which would place planet P1 comfortably in the type II regime. Baruteau et al. (2011) however suggest that GI protoplanets will migrate inwards much faster than the gap opening timescale, and that their migration may be better explained in the type I regime. Δ mig,I and Δ mig,II are likely more representative of lower and upper limits on the migration timescale of planet P1, and the subsequent migration of a GI protoplanet will be best explained by a combination of both regimes. In either case, these simple calculations demonstrate that, to first approximation, it should be entirely possible for a fragment formed on a wide orbit to migrate inward to the current location of planet P1 within the lifetime of the AB Aurigae disc.

Viscous evolution models of AB Aurigae
Despite the system's disc mass being too low to be gravitationally unstable currently, it will likely have been much more massive in the past prior to depletion by stellar accretion and photoevaporative winds, as massive discs will rapidly evolve away from an initially high mass state (Hall et al. 2019). Viscous evolution models use analytic prescriptions to calculate the evolution history of a protoplanetary disc's surface density profile. Hence, we may use them to predict the mass evolution history of AB Aurigae.

Methods
Full details of the model used here to calculate the evolution of a disc whose primary source of viscosity is provided by self-gravity can be found in Rice & Armitage (2009). We also outline the basic equations here.
Viscous evolution of the surface density, Σ( , ), can be modelled using the one-dimensional prescription from Lynden- Bell & Pringle (1974); Pringle (1981), where Σ wind represents the photoevaporative mass loss due to radiation from the central star. Here we implement the x-ray photoionization model described in detail in Owen et al. (2011) and assume a moderate x-ray luminosity of 1 × 10 30 erg s −1 , noting that the rate of photoevaporative mass loss scales linearly with x-ray luminosity. Disc viscosity, ( , ), is modelled using the Shakura-Sunyaev viscous− prescription (Shakura & Sunyaev 1973), where the disc scale height, = s /Ω, and Ω = √︁ * / 3 in a rotationally supported disc. The sound speed, s , is calculated by solving Equation 1, where we force the disc to be in a marginally unstable state with = 1.5.
The volume density can then be calculated as = Σ/2 , and the temperature, , optical depth, , and ratio of specific heats, , can be determined by interpolation of the equation of state table from Stamatellos et al. (2007) using the Rosseland mean opacities from Bell & Lin (1994).
To calculate the viscous-term from Equation 17, we must first determine the disc cooling time, which requires that we calculate the radiative cooling term (Hubeny 1990), and determine the local cooling time as cool = /Λ, where the energy per unit surface area is, In a disc where the primary source of viscosity comes from self-gravity, the effective viscous-term can be calculated as, (20) We set a lower limit, min , below which we assume that GI is not the dominant source of viscosity but instead, in a sufficiently ionized disc, MRI may dominate, for example. If < min we set = min and recalculate the disc properties, now no longer requiring the disc to be gravitationally unstable with = 1.5. Equations 1 and 20 can then be solved to calculate s and for use in Equation 17, and Equation 16 can be integrated to determine the time evolution of the disc's surface density profile, hence its mass evolution.

A note on the current mass of the AB Aurigae disc
Protoplanetary disc masses are notoriously challenging to measure. They often rely on empirical conversions between a disc's flux density and its mass, which requires uncertain assumptions about the disc optical depth, metallicity, dust-to-gas ratio and grain size distribution. Combined with uncertainties in the flux measurement and distance toward the system, mass estimates may be uncertain by up to an order of magnitude, and are usually considered to represent lower bounds. Estimates of the disc mass surrounding AB Aurigae find a low mass disc, with d = 0.01 M and uncertainty up to a factor ≈ 10 (Andrews & Williams 2005;Corder et al. 2005;Piétu et al. 2005;Semenov et al. 2005).
The accretion rate onto the star may also provide us a with rough estimate of the disc mass, as it is indicative of the mass reservoir available to the star from the disc. A protoplanetary disc is expected to settle into a steady-state with a constant mass accretion rate (Pringle 1981 In a disc with sound speed profile, s = s,0 −0.25 , and surface density profile, Σ = Σ 0 −1 , the disc may have a radially constant viscous-given by, We can substitute in for s,0 by assuming a flattened disc with / = 0.1 at = 100 AU, and substituting = s /Ω, where is the local disc scale height. Similarly, we can substitute Σ 0 for the disc outer radius, out = 400 AU and disc mass to obtain an equation in terms of and d , We plot this equation in Figure 3   From Figure 3 we see that for a very-low mass disc ( d ≤ 0.1 M ) to have an accretion rate = 1.3 × 10 −7 M yr −1 would require a disc viscosity much higher than we would usually expect from a quasi-stable disc, with ≥ 0.1. If instead the disc is still massive, with d ≥ 0.1 M , the viscous-required to explain the high accretion rate drops significantly. In a quasi-stable disc we might typically expect ≈ 10 −2 − 10 −4 (Hartmann et al. 1998;Rafikov 2017).
We do not attempt to propose an exact disc mass for AB Aurigae here, but instead wish to highlight that in order to explain the system's high accretion rate may require that the disc is more massive than has previously been suggested, and that it is likely at least as massive as the upper bound of the current disc mass estimates. . Viscous− vs. disc mass for a steady-state disc with = 1.3 × 10 −7 M yr −1 , equal to the mass accretion rate measured in AB Aurigae. We calculate as a function of disc mass using Equation 23, which assumes that the disc has a radially constant viscous− .

Results
With this in mind we use these viscous evolution models to predict how long ago the AB Aurigae disc may have been massive enough to exceed the critical mass limit for fragmentation, where d = crit = 0.3 M , assuming the system to have a current disc mass approximately equal to the upper bound on the mass estimate, d = 0.1 M .
In order to do this, we set up discs with initial masses d = crit = 0.3 M and evolve them forward in time until their mass has been depleted to d = 0.1 M , assuming min values in the range 0.01 − 0.05. Discs are set up with initial parameters similar to what we might expect in a young AB Aurigae disc, with * = 2.4 M , out,init = 400 AU, surface density profile Σ ∝ −1 and temperature profile ∝ −0.75 . We assume again that irradiation leads to a constant background temperature irr = 10 K.
The results of these models are shown in Figure 4. To illustrate how long in the recent past the AB Aurigae disc may have been massive enough to exceed the fragmentation threshold, we have plotted the disc mass evolution in reverse order. Hence = 0 represents the disc in its current state, with d = 0.1 M , and the x-axis measures Myrs in the past. For example, in the case of the min = 0.05 model, we predict that the AB Aurigae disc may have been more massive than d,crit = 0.3 M approximately 1.3 Myrs ago. This approach is equivalent to if we had run the models in Section 3.2.1 backwards, beginning at d = 0.1 M .
Discs with higher viscous− values will evolve at a faster rate, hence the time between the disc mass being in its current state, with d = 0.1 M , and exceeding the critical mass limit, d,crit = 0.3 M , will be shorter.
These models again reiterate how it is challenging to reconcile the current low estimated disc mass with the high measured accretion rate, leading us to conclude that AB Aurigae is either currently more massive than observations suggest, or was almost certainly so in its recent past. In the highest accreting case, with min = 0.05, we find the accretion rate when d = 0.1 M to be = 4.80 × 10 −8 M yr −1 , and in the lowest accreting case with min = 0.01 we find = 9.0 × 10 −8 M yr −1 , both of which are significantly lower than the currently measured value of = 1.3 × 10 −7 M yr −1 (Salyk et al. 2013).
Crucially though, the plots in Figure 4 demonstrate how we can trace the AB Aurigae disc back to a previously higher mass state, and how the disc mass may have exceeded the fragmentation threshold in the recent past. When assuming a moderate background min we find that the disc mass may have exceeded d,crit within the past ≈ 1.25 − 4 Myr. Thus it is plausible that a young AB Aurigae disc may have fragmented to form one or multiple giant gaseous protoplanets during its early evolution.

Jeans mass in an AB Aurigae-like disc
The local Jeans mass in a self-gravitating disc can be used, in the case where a region of the disc fragments, to estimate the masses of the bound clumps that will form.
We can use the same approach as in Section 3.2.1 to calculate how the Jeans mass varies as a function of and out . We assume the disc to be marginally unstable, with = 1.5, and use equation 1 to obtain s , and solve equation 20 to obtain for use in equation 21, allowing us to calculate the Jeans mass for a range of disc outer radii and accretion rates.
In Figure 5 we plot equation 24, for a disc around a 2.4 M star, with between 1 × 10 −9 M yr −1 and 1 × 10 −4 M yr −1 , and out between 50 AU and 500 AU. We assume that disc irradiation Hence the x-axis measures Myrs in the past. We vary the value of min , which represents a background viscous-value generated by some process other than disc self-gravity.
leads to a constant background temperature, and consider two cases where irr = 10 K and irr = 50 K. Higher disc temperatures reduce the effective viscous− from Equation 20 for discs of the same mass, whilst also providing greater pressure support against direct collapse, thus stabilising the system against GI. Hence, for a given and out the Jeans mass increases as a function of irradiation. A gravitationally unstable disc may fragment if a collapsing clump is able to cool and radiate energy away at a rate faster than the local dynamical time. This condition can be expressed in terms of a critical value of the dimensionless cooling parameter, = cool Ω, which in turn can be expressed in terms of a critical viscous− (see equation 20). We typically expect this value to be somewhere between crit ≈ 0.06 − 0.1 (Gammie 2001;Rice et al. 2005;Baehr et al. 2017), thus we include contours of = 0.01 and = 0.1 in Figure 5 to indicate regions of the parameter space that may fragment.
These plots reiterate that at an earlier stage of AB Aurigae's evolution, when the mass accretion rate was likely higher than it currently is, it is entirely plausible that the disc may have been gravitationally unstable and may have fragmented, as these higher accretion rate states lie in an unstable region of parameter space.
For a given disc radius the minimum Jeans mass doesn't vary much whether we assume fragmentation can only occur for ≥ 0.01 or ≥ 0.1. Assuming that crit = 0.1 we find the minimum Jeans masses at = 200 AU, 300 AU and 400 AU to be 1.6 M Jup , 2.5 M Jup and 3.4 M Jup respectively when irr = 10 K, and to be 10.3 M Jup , 12.4 M Jup and 13.3 M Jup respectively when irr = 50 K, roughly coinciding with what we observe from the mass of planet P1.

Implications for formation through CA
Significant fine tuning of the model parameters is required in Section 2 to form planet P1 through CA within the strict time constraint of the system's measured age. To form a planet of 4 M Jup within Figure 5. The Jeans mass in a self-gravitating disc surrounding a 2.4 M star. We consider two cases of disc irradiation, one where it leads to constant background temperature of irr = 10 K (left) and one where it leads to constant background temperature of irr = 50 K (right). We expect a disc to be unstable against fragmentation for crit ≈ 0.06 − 0.1, thus we plot contours of = 0.01 and = 0.1 to indicate regions of parameter space which would likely be unstable against fragmentation. Higher temperatures act to stabilize the disc against the gravitational instability by reducing the effective-. Hence for a given and R out the Jeans masses will be higher when irr = 50 K compared to when irr = 10 K.
1 − 4 Myr generally requires a planetesimal surface density much higher than would usually be expected, with a total planetesimal mass across the disc ≥ 0.072 M when M core,init = 0.01 M , and ≥ 0.024 M when M core,init = 0.1 M . When Σ p,5AU = 2 gcm −2 , hence with a total planetesimal mass across the disc of 0.012 M , we generally see very slow planet growth. It is important to note however that we have only considered a simple formalism for our modelling of CA here, and that processes not included in our models, such as planet migration, pebble accretion and disc instabilities, may be capable of accelerating initial growth. We discuss the effect of these next.

Limitations of the CA models
Migration allows the planet to sample a wider region of the disc, therefore preventing the local planetesimal surface density becoming depleted as rapidly as when it grows in-situ. When we include core migration in Section 2 the planets generally grow at a faster rate. However we chose to only consider in-situ formation here, as including migration causes all the cores to migrate to the inner disc ( 3 AU) away from the location where we currently find planet P1, and toward the regions of higher planetesimal surface density where they accrete at a faster rate. Some other mechanism, such as planet-planet scattering, would then be required to explain planet P1's subsequent migration out to ≈ 30 AU. When modelling insitu formation at the current semi-major axis of planet P1, we see only slow growth when core,init = 0.1 M ⊕ , and almost no growth when core,init = 0.01 M ⊕ .
Instabilities in discs may be capable of generating large overdensities of solids, hence they have been suggested as possible mechanisms for accelerated planetesimal growth and, in extreme cases, fragmentation of the disc solids under their self-gravity. The spiral arms of young, GI discs have been shown to cause strong dust-trapping (Rice et al. 2004), whilst the gravitational collapse of filaments generated in the streaming instability (Youdin & Goodman 2005;) may form planetesimals of radii 100−1000 km (Johansen et al. , 2011(Johansen et al. , 2012, thus providing a possible mechanism for the initial formation of rocky cores.
Whilst refraining from modelling the detailed physics of dust trapping through disc instabilities, we can crudely represent local grain enhancements by simply increasing the total dust-to-gas ratio in the disc, which by default will increase the planetesimal surface density local to the accreting core. We account this by increasing the total planetesimal surface density by up to a factor of 6 in Section 2.
Mechanisms for accelerated growth and rapid core formation become necessary as CA faces challenges when establishing how the first planetesimals are able to grow beyond metre sizes. The initial stages of growth are believed to be slow, as dust grains may encounter growth barriers beyond metre sizes (Brauer et al. 2008;Mordasini et al. 2010). It has been shown, as consequence of intrinsic gas-dust drag in the disc, that grains of a critical size will radially migrate and be accreted onto the star within a fraction of the disc lifetime (the radial drift barrier, Weidenschilling 1977). Further, solids of millimetre to centimetre sizes, with Stokes number close to 1, are expected to have high relative azimuthal velocities, hence grain-grain collisions may become destructive, resulting in shattering (the fragmentation barrier, Birnstiel et al. 2012), or neutral and result in recoiling (the bouncing barrier, Zsom et al. 2010), both of which prevent a positive outcome of coagulation. In our model we assume that a core of mass 0.01 M ⊕ or 0.1 M ⊕ , with core,init = 1.6 × 10 3 km and core,init = 3.5 × 10 3 km respectively, has already formed at = 0, therefore avoiding the detailed physics of this initial phase of core growth. Note that these initial core sizes are consistent with, but slightly larger than, the planetesimals expected to form through direct collapse of the dust disc during the streaming instability (Johansen et al. , 2011(Johansen et al. , 2012. Possibly most importantly, we note that we do not include a prescription for pebble accretion in our model (for a review see Johansen & Lambrechts 2017). Accretion of millimetre to centimetre sized pebbles onto planetesimal cores may have the potential to generate significantly faster growth rates than the planetesimalplanetesimal accretion we consider here. Pebbles may be abundant in protoplanetary discs, since it is a natural outcome from the fragmentation and bouncing barriers. Pebbles of millimetre-centimetre sizes are coupled to the gas in the disc. The gas component orbits at sub-Keplerian velocities due to the outward gas pressure. The solids, which are orbiting at Keplerian velocities, will experience a drag force which, in a smooth, laminar disc, will cause them to radially drift inward. This migration of pebbles can lead to them being transported to within the path of the growing planetesimal core, constantly replenishing the pebbles within the planetesimal's feeding zone and preventing it from reaching its isolation mass as quickly as they do in Section 2. If the planetesimal is gravitationally massive and capable of perturbing the velocities of nearby solids, the pebbles may enter into complex trajectories, orbiting and eventually settling down into its gravitational potential well. If the planetesimal's gravitational cross section exceeds its geometric cross section, pebble accretion may become the dominant growth mechanism. In their review paper Johansen & Lambrechts (2017) show that pebble accretion may be capable of resolving many of the timescale problems associated with CA, whilst being able to explain the formation of all planet types.

Implications for formation through GI
In Section 3.1 we used SPH simulations to determine the critical mass limit for fragmentation in a disc surrounding a 2.4 M star, finding that for a out = 400 AU disc, d,crit = 0.3 M ( crit = 0.125). Whilst we have mostly focused our discussion on the case of single fragment formation from our SPH simulations, it is also likely that multiple clumps may form in a disc with a mass slightly higher than d,crit (see Figure 2). The initial formation of multiple protoplanets may then also provide an explanation for the widerorbit planet P2 which has also been inferred, located at a distance ≈ 140 AU from the parent star (Boccaletti et al. 2020). We have refrained from analysing the formation history of planet P2, due to its slightly more tentative detection, choosing instead to focus on planet P1. However it would seem that the formation of a 3 M Jup planet at ≈ 140 AU may be even more challenging to explain in the CA paradigm than is the case for planet P1, as the gas and dust surface densities in the disc will drop off as Σ ∝ −1 , hence will be exceedingly low at such a large radius. As we see only minimal core growth at = 30 AU in Figure 1d, it is likely that growth at = 140 AU would be near-negligible. It may then be the case that in fact planets P1 and P2 represent two survivors from several fragments which could have initially formed.
Despite the AB Aurigae disc being far too low mass to be gravitationally unstable currently, models of the system's viscous evolution in Section 3.2 suggest that it may have been much more massive when it was younger, potentially exceeding the critical mass limit for fragmentation. It seems reasonable to expect that the disc might have previously fragmented in an extended system such as AB Aurigae, as previous studies suggest that fragmentation is inevitable in GI discs at radii, 50 − 100 AU (Rafikov 2005;Whitworth & Stamatellos 2006;Clarke 2009;Forgan & Rice 2011). Further, Cadman et al. (2020b);Haworth et al. (2020) used hydrodynamic simulations to demonstrate that, whilst lower mass stars may support gravitationally stable massive discs, susceptibility to fragmentation increases as a function of stellar mass, and that discs around higher mass stars ( * ≥ 2 M ) may fragment for relatively low disc-tostar mass ratios. AB Aurigae being an extended disc around a higher mass star therefore seems to be an ideal candidate system to search for surviving products of GI.
If the disc had been able to fragment whilst it was young, it is not necessarily true that the clumps will have survived the 1 − 4 Myr lifetime of the AB Aurigae system. We find that fragments may initially form on wide-orbits with 200 AU, and use analytic calculations to predict initial clump masses 1.6−13.3 M Jup .
However subsequent evolution is inevitable, and the fragments will rapidly migrate through the disc (Baruteau et al. 2011).
In the tidal downsizing hypothesis of planet formation (Nayakshin 2010a,b, 2011) GI embryos will cool and contract as they migrate. Dust sedimentation may lead to the formation of a solid core, potentially of mass comparable to that of a terrestrial planet (Boss 1998). If the embryo's outer layers contract slowly whilst migration occurs rapidly then tidal stripping from the parent star may occur once the embryo reaches the inner disc, as its physical radius may exceed its Hill sphere (Nayakshin 2010a). It is possible that many of the initially formed fragments may be entirely destroyed during this tidal downsizing process (Nayakshin & Fletcher 2015;Humphries et al. 2019). Population synthesis calculations find this may be the true of ≈ 50% of GI protoplanets, with the remaining objects eventually residing at 20 AU (Forgan & Rice 2013b), although when including fragment-fragment scattering this survival fraction may be significantly less (Forgan et al. 2018). The initial formation of multiple clumps would then be necessary if any are to survive beyond this early phase of evolution. Accretion of material onto the protoplanets will also occur as they migrate through the disc. Kratter et al. (2010) showed that most GI fragments will grow well beyond the mass limit for Deuterium burning, and that any GI-born planets likely represent the low mass tail of the eventual GI fragment mass distribution. The Jeans mass estimates that we present in Section 3.3 therefore represent those shortly after collapse only, as dynamical evolution will significantly influence the embryo's eventual mass.
We also tentatively suggest that the previous disc mass estimates ( d ≈ 0.01 M ) (DeWarf et al. 2003;Andrews & Williams 2005;Corder et al. 2005;Semenov et al. 2005) appear too low to be consistent with the high stellar accretion rate (Salyk et al. 2013), which is indicative of the presence of a large mass reservoir. Assuming the disc to be in a quasi-steady state with a radially constant viscous− suggests a lower limit for the current disc mass as d 0.1 M (see Fig. 3). This rough lower limit is in fact consistent with the upper bound of the uncertainty on the current disc mass estimates. However even when assuming this slightly higher disc mass, we still find the calculated accretion rates from our viscous evolution models in Section 3.2.3 to be significantly lower than the accretion rate measured from the system. On the unusually high stellar accretion rate, Tang et al. (2012) suggest a possible explanation is the presence of an inner disc, characterised by a gas/dust cavity observed at ≈ 100 AU, which is being replenished through accretion from the remnant envelope above and below the disc midplane. This would suggest that the measured accretion rate does not represent that of a settled, out = 400 AU disc as we have assumed here, and would allow for the existence of a low mass disc whilst being consistent with a high accretion rate. We only attempt to further highlight this discrepancy between the measured disc mass and accretion rate, and note that the current mass of the disc does not significantly affect the overall conclusions from this paper in regards to the formation history of planet P1.

CONCLUSIONS
In this paper we have analysed the possible formation history of the 4 − 13 M Jup planet observed at ≈ 30 AU within the protoplanetary disc surrounding AB Aurigae (Piétu et al. 2005;Tang et al. 2012Tang et al. , 2017Boccaletti et al. 2020). The young age of the star-disc system places strict constraints on the CA formation timescale, which we find challenging to explain within 1 − 4 Myr. The planet's high mass and wide-orbit are indicative of a planet which may have instead formed through disc instability in the natal AB Aurigae disc.
The key results are as follows.
(i) Typical in-situ CA formation timescales for planet P1 exceed the system's measured age. Fine tuning of the model parameters is required in order to form a planet of 4 M Jup within 1 − 4 Myr, including significant enhancement of the planetesimal surface density in the disc, and, in most cases, that a large planetesimal core with core,init = 0.1 M ⊕ has already formed near to the snow line at = 0. At the current semi-major axis of planet P1 ( = 30 AU) we find extremely slow in-situ growth due to the low disc surface densities at wide orbits. We do not include a prescription for pebble accretion in our models here, but note that it may be capable of speeding up planet growth.
(ii) A disc surrounding a 2.4 M star, analogous to young AB Aurigae, would have fragmented if its initial mass exceeded d,crit = 0.3 M ( crit = 0.125). If the disc mass is slightly higher than d,crit several fragments may form. Formation of multiple fragments would allow margin for some fragment destruction, which is likely inevitable during their subsequent dynamical evolution of GI protoplanets.
(iii) Viscous evolution models of the AB Aurigae disc suggest that it may have been massive enough to exceed d,crit during its early evolution whilst the disc was still young and massive. We find that a 0.1 M disc may have exceeded d,crit = 0.3 M within the past ≈ 1.25 − 4 Myr when considering moderate background viscosity.
(iv) Fragments will initially form on wide orbits, where the disc material is cool, and then rapidly migrate inwards. Typical migration timescales of a GI protoplanet which formed at ≈ 200 AU within a young AB Aurigae disc are found to be shorter than the current age of the system. We use analytic calculations to determine type I and type II migration timescales, finding that for migration from out = 200 AU to in = 30 AU, Δ mig,I = 6.9 kyr and Δ mig,II = 1.0 Myr when considering disc conditions taken from our hydrodynamic simulations.
(v) Calculations of the Jeans mass in a moderately irradiated proto-AB Aurigae disc represent what the initial fragment masses might have been immediately after formation. We find that J = 1.6 − 13.3 M Jup , which is consistent with the masses of the planets P1 and P2 in the AB Aurigae disc.
(vi) Although we focus our discussion on the formation history of planet P1, we highlight that planet P2 found at ≈ 140 AU with an estimated mass P2 = 3 M Jup may be even more challenging to reconcile with formation through CA.
We therefore propose that planets P1 and P2 which have been inferred through scattered light observations of the AB Aurigae disc (Boccaletti et al. 2020) may stand as evidence of planet formation through GI.