Grand Challenges in Protoplanetary Disc Modelling

Abstract

The Protoplanetary Discussions conference—held in Edinburgh, UK, from 2016 March 7th–11th—included several open sessions led by participants. This paper reports on the discussions collectively concerned with the multi-physics modelling of protoplanetary discs, including the self-consistent calculation of gas and dust dynamics, radiative transfer, and chemistry. After a short introduction to each of these disciplines in isolation, we identify a series of burning questions and grand challenges associated with their continuing development and integration. We then discuss potential pathways towards solving these challenges, grouped by strategical, technical, and collaborative developments. This paper is not intended to be a review, but rather to motivate and direct future research and collaboration across typically distinct fields based on community-driven input, to encourage further progress in our understanding of circumstellar and protoplanetary discs.

1 INTRODUCTION

For the first time in history, spatially resolved observations of the structures within protoplanetary discs are being obtained (see review by Casassus 2016). This has revealed a wealth of sub-structure, including rings and gaps (ALMA Partnership et al. 2015; Andrews et al. 2016; Canovas et al. 2016), spirals (e.g. Garufi et al. 2013; Benisty et al. 2015; Wagner et al. 2015), warps (e.g. Casassus et al. 2015), shadows (e.g. Stolker et al. 2016), cavities (e.g. Andrews et al. 2011), and dust traps (e.g. van der Marel et al. 2013, 2016). These recent observations, combined with the huge diversity of exoplanetary systems discovered over recent years (Winn & Fabrycky 2015), has stimulated a new wave of rapid development in the modelling of protoplanetary discs, to better understand their evolution, along with their connection to the planet formation process (e.g. Papaloizou & Terquem 2006).

Understanding the evolution of discs, the structures that we are observing within them and the planet formation process presents a substantial challenge to modellers. Discs are composed of non-primordial material spanning conditions ranging from cold, extremely dense, and molecular, through to diffuse, hot, and ionised. Densities and temperatures vary by ~ 10 and 3 orders of magnitude, respectively. The basic chemical composition of discs alone is the subject of at least four complex research fields distinguished by the local matter conditions and radiation field: dust–grains, gas–grain chemistry, photon-dominated chemistry, and photoionisation (e.g. Gorti & Hollenbach 2009; Thiabaud et al. 2015; Walsh, Nomura, & van Dishoeck 2015; Gorti, Hollenbach, & Dullemond 2015). The situation is even more challenging since the observational determination of a disc’s composition is often degenerate, making direct comparison between observations and theory (and thus validation of our models) difficult (e.g. Meijer et al. 2008; Woitke et al. 2016; Miotello et al. 2016; Boneberg et al. 2016; Kama et al. 2016).

The dynamics of protoplanetary discs are also extremely challenging. The gravitational potential from the parent star, self-gravity of the disc, hydrodynamic torques in the disc, radiation from the parent star or other nearby stars, dust, and (non-ideal) magnetohydrodynamics (MHDs) all play important roles (Bodenheimer 1995; Dullemond et al. 2007; Lodato 2008; Armitage 2011, 2015). Furthermore, the dynamical evolution of dust grains with moderate Stokes numbers St ≳ 0.01 must be solved in addition to the gas dynamics (for a recent review, see Testi et al. 2014). Discs are also not necessarily in a steady state, and can be subject to a range of instabilities, such as gravitational fragmentation (Durisen et al. 2007; Young & Clarke 2015; Forgan, Parker, & Rice 2015; Meru 2015; Takahashi, Tsukamoto, & Inutsuka 2016), the streaming instability (Youdin & Goodman 2005), Rossby wave instability (e.g. Lovelace et al. 1999; Tagger 2001; Lyra et al. 2008b, 2009), baroclinic and vertical shear instabilities, which can form and grow vortex structures (Lyra & Klahr 2011; Lesur & Papaloizou 2010; Nelson, Gressel, & Umurhan 2013; Richard, Nelson, & Umurhan 2016), the magneto-rotational instability (e.g. Balbus & Hawley 1991; Reyes-Ruiz et al. 2003), and dust-settling induced vortices (Lorén-Aguilar & Bate 2015, 2016). The local environment can also significantly modify disc evolution via mass transfer from the ambient medium onto the disc (Vorobyov, Lin, & Guedel 2015; Lomax, Whitworth, & Hubber 2015), nearby radiation sources (e.g. Bally, O’Dell, & McCaughrean 2000; Henney et al. 2002; Smith, Bally, & Morse 2003; Adams et al. 2004; Wright et al. 2012; Facchini, Clarke, & Bisbas 2016) and tidal encounters (e.g. Clarke & Pringle 1993; de Juan Ovelar et al. 2012; Rosotti et al. 2014; Vincke, Breslau, & Pfalzner 2015; Dai et al. 2015; Vincke & Pfalzner 2016). A summary of some of the key processes (local, not environmental) that modellers attempt to capture in discs is given in Figure 1.

Figure 1. A protoplanetary disc schematic highlighting some of the key disc mechanisms and physics we are required to model to capture them (in parentheses). These physical ingredients are hydrodynamics (HD), magnetohydrodynamics (MHD), radiation hydrodynamics (RHD), radiative transfer (RT), chemistry (CHEM), and dust dynamics (DD). The background image is a subset of a Hubble observation of R136, credit: NASA, ESA, and F. Paresce (INAF-IASF, Bologna, Italy), R. O’Connell (University of Virginia, Charlottesville), and the Wide Field Camera 3 Science Oversight Committee.

This physically rich environment is made even more complex given that most of these dynamic, magnetic, radiative, and chemical processes are interlinked. For example, the effect of magnetic fields depend upon the ion density, which in turn is determined by the composition, which in turn depends upon the radiation field (e.g. due to photoionisation of atoms, photodissociation of molecules and determination of the thermal properties through processes such as line and continuum cooling). Another distinct coupling is the interaction between the gravitational instability and the magnetorotational instability, which has been well-studied in the disc community using semi-analytic models as the cause of an accretion limit cycle causing protostellar outburst phenomena (Armitage, Livio, & Pringle 2001), but is only now being investigated with self-consistent hydrodynamic simulations (e.g. Bae et al. 2014). Another example is that the radiation field in a disc is sensitive to the distribution of small dust grains (the motions of which may also be influenced by the radiation field, e.g. Hutchison et al. 2016) which in turn is sensitive to dynamical effects such as shadowing caused by warping of the inner disc (Marino, Perez, & Casassus 2015; Stolker et al. 2016). Furthermore, radiative heating increases the gas sound speed, and hence the amount of turbulent motion transferred to dust grains via gas-dust coupling, which influences grain–grain collisions and therefore the growth and fragmentation of dust (e.g. Testi et al. 2014). As a final example, gravitational instability and fragmentation in discs is sensitive to radiation (e.g. Meru & Bate 2010; Forgan & Rice 2013) and magnetic fields (Price & Bate 2007; Wurster, Price, & Bate 2016), and can induce dramatic effects in the chemical composition of discs (see Section 3, Ilee et al. 2011; Evans et al. 2015).

Given the importance of these links, ultimately one wishes to identify which physical processes affect each other in a non-negligible fashion, and to model all of them simultaneously. The modelling of protoplanetary discs is therefore a daunting task—what might be termed a grand challenge. Each physical mechanism requires sufficient rigour and detail that modelling them constitutes an active field of protoplanetary disc research in their own right (for reviews of physical processes in protoplanetary discs, see e.g. Hartmann 1998; Armitage 2011; Williams & Cieza 2011; Armitage 2015). In practice, we have neither the numerical tools nor computational resources to achieve such multi-physics modelling of protoplanetary discs at present (nor in the immediate future). However, we can set out a roadmap towards this goal whilst outlining the more achievable milestones along the way.

In this paper, motivated by group discussion sessions at the ‘Protoplanetary Discussions’ conference in Edinburgh ¹ , we ultimately aim to stimulate progress in the multi-physics modelling of protoplanetary discs in order to deepen our understanding of them. This paper is presented in parallel with a second paper which focusses on the observations required to advance our understanding of discs (Sicilia-Aguilar et al., in preparation). Although our focus here is new numerical methods and the questions they might answer, it is important to remember that there are still many unsolved problems that can be tackled with existing techniques. Additionally, new numerical methods are likely to be computationally expensive so there will be many problems that are better tackled using existing techniques (e.g. parametric models used to interpret observations, Williams & Best 2014). Furthermore, this paper is not exhaustive, there will certainly be fruitful avenues of theoretical research into protoplanetary discs that are not discussed here (in particular regarding magnetic fields and the details of planet formation itself).

The structure of this paper is as follows—in Section 2, we provide an overview of some core ingredients of disc modelling. In Section 3, we then present a series of mid- and long-term challenges to motivate future development. Finally, in Sections 4–6, we discuss pathways towards meeting the challenges in terms of strategical, technical, and collaborative developments.

2 AN OVERVIEW OF CURRENT TECHNIQUES

We begin by providing an overview of some of the core ingredients of protoplanetary disc modelling, to introduce concepts and provide context for the rest of the paper. This is by no means intended to be a comprehensive review, rather it should provide some basic platform from which a reader unfamiliar with certain concepts can proceed through the rest of the paper. Figure 2 illustrates the four core disciplines that comprise the majority of protoplanetary disc modelling: gas and dust dynamics, magnetic fields, radiative transfer, and chemistry. As shown, these topics are all fundamentally linked. It is this interdependence that raises the possibility that multi-physics modelling will be important and is hence a key focus of this paper.

Figure 2. An illustration of the core disciplines in protoplanetary disc modelling: gas and dust dynamics, magnetic fields, radiative transfer, and chemistry. Each discipline is a field in its own right, subject to intensive study. However, they are all closely interlinked, affecting each other in a number of ways, of which we illustrate a few representative examples. It is this interdependence between fields that necessitates the drive towards multi-physics modelling of protoplanetary discs.

2.1. (Magneto-) hydrodynamics

Solving for the motion of fluids as a function of time is a key ingredient for understanding the evolution of protoplanetary discs. Numerical hydrodynamics is a relatively mature field. Numerical solvers are either Eulerian or Lagrangian in character. Eulerian solvers trace flows across fixed discrete spatial elements, whilst Lagrangian solvers follow the motion of the flow. In protostellar disc simulations, the majority of hydro solvers are either Eulerian/Lagrangian grid-based simulators, or the fully Lagrangian Smoothed Particle Hydrodynamics.

Depending on the resolution requirements, solvers are either global, in that the entire disc extent is simulated together, or local, where a region in the disc is simulated at high resolution, with appropriate boundary conditions to reflect the surrounding disc environment. Which construction is best used is dependent upon the problem being studied, as we discuss below.

2.1.1. Global disc simulations

Historically, the primary challenge for global simulations of protoplanetary discs with Eulerian codes was the Keplerian flow—advection of material at supersonic speeds across a stationary mesh is a recipe for high numerical diffusion. This has now been overcome with, for example, the fargo algorithm (Masset 2000), implemented in both the fargo (Masset 2000; Baruteau & Masset 2008; Benítez-Llambay & Masset 2016) and pluto (Mignone et al. 2007) codes. Eulerian codes perform best when the flow is aligned with the grid. This means that cylindrical or spherical grids are preferable which, when applicable, offer the best accuracy currently possible of any technique for a given level of computational expense or resolution. However, this means that adaptive mesh refinement (Berger & Colella 1989), which is mainly (but not exclusively) developed for Cartesian meshes, is not typically used (an example of an exception is Paardekooper & Mellema 2004). Furthermore, simulating warped, twisted, or broken discs remains difficult (e.g. Fragner & Nelson 2010).

Lagrangian schemes such as smoothed particle hydrodynamics (SPH, for reviews, see e.g. Monaghan 1992; Price 2012) are well suited to more geometrically complex global disc simulations because advection is computed exactly, angular momentum can be exactly conserved (e.g. an orbit can be correctly simulated with one particle) and there is no preferred geometry. Numerical propagation of warps using SPH has been shown to closely match the predictions by Ogilvie (1999) of α-disc theory (Lodato & Price 2010). In particular, a generic outcome of discs that are misaligned with respect to the orbits of central binaries or companions is that the disc ‘tears’ (Nixon et al. 2012; Nixon, King, & Price 2013; Nealon, Price, & Nixon 2015) or breaks (Nixon & King 2012; Facchini, Lodato, & Price 2013; Doğan et al. 2015). Such behaviour is well modelled by SPH codes, and appears to be relevant to observed protoplanetary discs, including HK Tau (Stapelfeldt et al. 1998), KH15D (Lodato & Facchini 2013), and HD142527 (Casassus et al. 2015). A limitation of the SPH approach is that the particles adaptively trace the densest regions, low density components of the disc, e.g. gaps and the disc upper layers, can therefore be under-resolved (e.g. de Val-Borro et al. 2006).

2.1.2. Local simulations

The most common technique utilised for local simulations of discs is the Cartesian shearing box (Hawley, Gammie, & Balbus 1995; Guan & Gammie 2008). This imposes the shear flow in a subset of a disc and allows for high-resolution simulations of disc microphysics in a Cartesian geometry, well suited to most Eulerian codes. This means that all the sophistication of modern Godunov-based hydrodynamics can be applied (there are many textbooks covering grid-based hydrodynamics, e.g. Toro 2013). This approach has been used almost exclusively for simulating the magnetorotational (see Balbus 2003, and references within) and other instabilities—in particular the streaming instability (e.g. Youdin & Goodman 2005; Youdin & Johansen 2007; Johansen et al. 2007; Bai & Stone 2010b)—in discs. Though other applications include the study of magnetically driven disc winds (e.g. Suzuki & Inutsuka 2009; Suzuki, Muto, & Inutsuka 2010).

By contrast, at present, there is no particular advantage to use Lagrangian schemes for local disc simulations. The cost for comparable results in cartesian boxes is up to an order of magnitude higher in SPH compared to Eulerian codes (e.g. Tasker et al. 2008; Price & Federrath 2010), mainly due to the additional costs associated with finding neighbouring particles, and the algorithms tend to be more dissipative than their grid-based counterparts, particularly when the flow is well matched to the grid geometry. However, Lagrangian techniques can accommodate open boundary conditions more naturally, so may offer advantages for certain problems in the future.

2.1.3. Other codes

In recent years, several new hydrodynamic solver methods have appeared. This broad class of Arbitrary Lagrangian Eulerian methods (ALE) offer the user the ability to switch between Lagrangian and Eulerian formalisms smoothly, in some cases during simulation runtime. Such ALE solvers include moving mesh codes (Springel 2010, 2011; Duffell & MacFadyen 2011) and meshless codes (Maron, McNally, & Mac Low 2012; McNally, Maron, & Mac Low 2012; Hopkins 2015). This extreme flexibility in approach appears to offer highly conservative schemes and adaptive resolution whilst capturing mixing and shear instabilities with high fidelity. The relative youth of these techniques (at least, in their application to computational astrophysics) means the full extent of weaknesses and strengths in these approaches remains to be seen (e.g. the ‘grid noise’ encountered during mesh regularisation; Mocz et al. 2015) although early applications to protostellar discs appear to be promising (see e.g. Muñoz et al. 2014).

Another recent development in numerical astrophysical fluid dynamics is the use of discontinuous Galerkin methods (which have a long history of application in the mathematical community). These grid-based techniques offer accurate, high-order solutions in a manner that is readily applied to adaptive meshes, and that scale efficiently on modern high performance computing facilities. In the astrophysical community, discontinuous Galerkin algorithms have now been implemented in both Cartesian (e.g. the tenet code; Schaal et al. 2015) and moving Voronoi mesh (e.g. the arepo code; Mocz et al. 2014) frameworks.

2.1.4. Magnetic fields

The above hydrodynamic solvers are able to include the evolution of the magnetic field in their fundamental equations. This has been most easily incorporated in Eulerian solvers, with mature MHD implementations in, for example, the athena (Stone et al. 2008), enzo (Bryan et al. 2014a), fargo (Benítez-Llambay & Masset 2016), pluto (Mignone et al. 2007), and pencil (Brandenburg & Dobler 2002) codes. SPH and other meshless codes can now also incorporate MHD (see review by Price 2012), provided that the ∇ · B = 0 condition can be sustained, for example, using divergence cleaning techniques (Tricco & Price 2012). Note however, that MHD with SPH is not a mature approach and is therefore somewhat less robust than Eulerian MHD at present (e.g. Lewis, Bate, & Tricco 2016).

Whilst ideal MHD disc simulations have been conducted for some time (see Balbus 2003, and references within), particularly important for protostellar discs is the role of non-ideal MHD, ever since the idea of a ‘dead zone’ was proposed by Gammie (1996). More recently, the interplay between the Hall effect, ambipolar diffusion, and Ohmic diffusion is yielding new turbulent behaviour (Sano & Stone 2002; Simon et al. 2015b), new forms of instability, and zonal flows in both MRI-active and ‘dead zone’ regions (e.g. Kunz & Lesur 2013; Bai & Stone 2014), not to mention addressing the so-called magnetic braking catastrophe that suppresses disc formation in ideal MHD (Tsukamoto et al. 2015; Wurster et al. 2016) (see recent review by Tsukamoto 2016, this volume).

For more general modelling of young stellar systems, global simulations are particularly important for modelling the launching of magnetised winds from the star and/or disc, and jets from the central star (e.g. Casse, Meliani, & Sauty 2007; Bai 2014; Lovelace et al. 2014; Suzuki & Inutsuka 2014; Staff et al. 2016).

2.1.5. Remarks on hydrodynamics

In summary, there are a number of options available as how to model the (magneto-)hydrodynamical evolution of a disc—the problem one is addressing determines which method is most appropriate. This ‘horses for courses’ approach is important, and is likely to extend to efforts which hope to further include elements from the other disciplines of disc modelling such as chemistry and radiation transport.

2.2. Dust–gas dynamics

The dynamics of small dust grains (Stokes number ≪ 1) is typically well coupled to that of the gas. For larger grains, however, the dust and gas dynamics can be decoupled. Properly modelling these decoupled motions is important both for disc dynamics, but also for interpreting observations. This latter point is particularly prudent given that some of the most important disc observations in recent years are millimetre continuum observations (i.e. of dust). For example, decoupled dust and gas dynamics is apparently important for understanding the symmetric gaps observed in discs (e.g. Dipierro et al. 2015b; Jin et al. 2016; Rosotti et al. 2016).

Approaches for modelling the dynamics of dust grains that are decoupled from the motions of the gas are often distinguished by whether they use a single or two-fluid approach, both of which we discuss below.

2.2.1. Two-fluid or ‘hybrid’ schemes

In an SPH framework, the two-fluid approach sees the dust and gas as separate particle populations, the dynamics for which are solved separately (Monaghan & Kocharyan 1995; Barrière-Fouchet et al. 2005; Laibe & Price 2012a, 2012b; Lorén-Aguilar & Bate 2014; Booth, Sijacki, & Clarke 2015). In grid-based methods, the dust is typically simulated as a particle population, with the hydrodynamics computed on the grid (e.g. Paardekooper 2007; Lyra et al. 2008a; Miniati 2010; Bai & Stone 2010a; Flock et al. 2015; Baruteau & Zhu 2016; Yang & Johansen 2016)—hence usually referred to as a ‘hybrid’ approach. The ‘hybrid’ or ‘two-fluid’ approaches are best suited to decoupled grains with Stokes number ≳ 1, where the interaction can be computed explicitly.

The traditional difficulty when dust is modelled by a separate set of particles is that short timesteps are required for small grains (Stokes numbers ≪ 1), requiring implicit timestepping schemes (Monaghan 1997; Bai & Stone 2010a; Miniati 2010; Laibe & Price 2012b). However, Laibe & Price (2012a) showed that simulating tightly coupled grains this way leads to ‘overdamping’ of the mixture, becoming increasingly inaccurate for small Stokes numbers, caused by the need to spatially resolve the ‘stopping length’ l ~ c _s t _s (where c _s is the sound speed and t _s is the stopping time). A similar issue was noted by Miniati (2010) in the context of grid based codes, finding only first-order convergence in the ‘stiff’ regime when the stopping time is shorter than the Courant timestep. However, by making use of the analytical solutions for the motion under drag forces that respect the underlying problem, this dissipation can be substantially reduced (or entirely avoided in the limit of negligible dust mass, Lorén-Aguilar & Bate 2014).

2.2.2. Single-fluid schemes

In the single-fluid approach, the dust parameters (dust to gas ratio, relative velocity) are properties of the ‘mixture’. In SPH, this means that a single population of SPH particles is used, representing the total fluid mass, with dust properties updated on each ‘mixture’ particle (Laibe & Price 2014a, 2014b, 2014c; Price & Laibe 2015; Hutchison et al. 2016). The same approach on a grid means evolving the dust density on the grid (called a ‘two-fluid’ approach by Miniati 2010—though not to be confused with the two-fluid approach mentioned above—to distinguish it from the ‘hybrid’ grid-plus-particles method). This is sometimes achieved using the approach suggested by Johansen & Klahr (2005) based on the ‘short friction time’ or ‘terminal velocity approximation’ for small grains. Here, the dust continuity equation is solved and the dust velocity is set equal to the gas velocity plus the stopping time times the differential forces between the gas and dust mixture. This is similar to the ‘diffusion approximation for dust’ derived by Laibe & Price (2014a) and implemented in SPH by Price & Laibe (2015) with an important caveat—that this formulation is only valid when the dust fraction is small (since it assumes that the gas velocity equals the barycentric velocity of the mixture). This assumption can easily be relaxed, at no additional computational expense, as shown by Laibe & Price (2014a).

An attraction of fluid-based dust models is that within their domain of validity, they provide a high degree of accuracy for their computational cost, whilst particle approaches typically suffer from sampling noise. However, the fluid approach is equivalent to use a moment-based method for solving the radiative transfer equations (see Section 2.3) where all moments of order greater than unity (or even zero in the short-friction time approach) are discarded. This means that in cases where the dust velocity becomes multi-valued the result may converge to the wrong answer. Possible examples of when this can occur include settling (for S_t > 1), turbulent motion (S_t ≳ R ^{− 1/2} _e ~ 10⁻⁴ in astrophysical flows, although Reynolds nymbers, R_e ≳ 10³, are rarely achieved numerically, Falkovich, Fouxon, & Stepanov 2002; Ormel & Cuzzi 2007), strong gravitational scattering, and in convergent flows at curved shocks. By including higher order moments, the fluid approximation could be extended to support multi-valued flows and thus support both large and small grains (Chalons, Kah, & Massot 2012; Chalons, Fox, & Massot 2010; Yuan & Fox 2011; Yuan, Laurent, & Fox 2012).

2.2.3. Dust post-processing approaches

Whilst the dynamical evolution of discs is clearly of importance to many problems, there are many cases in which the dynamic timescales are very different to other processes (see also Section 5.3 of this paper). For example, the short radiative timescale in discs has led to the standard approach of treating them as isothermal. Similarly, since dust growth often occurs on much longer timescales (> 10⁴ yr), the approach of post-processing the dust evolution according to some average over the short-term dynamics can be viable. For example, Brauer, Dullemond, & Henning (2008) and Birnstiel, Dullemond, & Brauer (2010) evolve the gas disc until a steady state is reached and then evolve the dust against this steady gas background.

This approach has also been applied to transition discs and discs with massive planets embedded, in particular following the growth of large particles trapped inside pressure maxima (Pinilla et al. 2015, 2016). Similarly, Dipierro et al. (2015a) applied this approach to self-gravitating discs in order to predict scattered light images. Miyake, Suzuki, & Inutsuka (2016) have also studied the motions of dust grains against a fixed gas background for the scenario of magneto-rotationally driven winds. However, we note that this approach can be fraught with difficulty, since it is difficult to know a priori what the representative average of the disc should be within which to evolve the dust. For example, particles with St ~ 1 can become trapped in the spiral arms of self-gravitating discs (or other pressure maxima), making azimuthal averaging unreliable. Similarly, although Rosotti et al. (2016) showed that azimuthal averaging works well for transition discs formed by planets of order, a Jupiter mass or less, ignoring the gas dynamics completely would predict an incorrect surface density profile and thus also incorrect growth rates. However, when the effects of combined dust–gas dynamics are taken properly into account (e.g. the short-friction time approximation can be used with hydrodynamic models to predict the evolution of dust grains 1 mm or smaller in transition discs), the post-processing approach will undoubtedly continue to provide important insights.

Conversely, coupling to live simulations of the dust/gas dynamics may prove to be essential for understanding some phenomena. For example, Gonzalez et al. (2015b) showed that by incorporating grain growth, radial drift, and feedback that self-induced dust traps may arise (to be explored in more detail in Gonzalez et al. in preparation). There will be many other important cases that likely require live simulations, for example, understanding whether planet formation can occur via the streaming instability in dust traps will require models that can show whether grains can grow to the required sizes without destabilizing the trap (e.g. Kato, Fujimoto, & Ida 2012; Taki, Fujimoto, & Ida 2016).

2.2.4. Remarks on dust dynamics

To date, there are virtually no simulations where both small and large grains are directly simultaneously evolved alongside the gas, in 3D, including the backreaction on the gas [though considerable progress towards this has been made by Paardekooper (2007), Lyra et al. (2008a), Gonzalez et al. (2015a, 2015b)]. Such a combination is important, because the grains, particularly when the dust-to-gas ratio becomes high, exert a backreaction on the gas, which in turn modifies the dynamics of the other grain species. For example, Laibe & Price (2014c) showed that under certain conditions, effects from the dynamics of multiple grain species could lead to the outward rather than inward migration of pebble-sized grains in discs. Whilst the large grain populations with S _t ≳ 1 are more interesting dynamically because they are more decoupled from the gas, modelling the small grains is necessary for coupling with radiative transfer and thus for comparison with observations. Paardekooper (2007) and Lyra et al. (2008a) do model a distribution of grain sizes using a particle appraoch, but not in regimes where the backreaction on to the gas is accounted for. Another often used approach is to perform a series of single grain-size simulations, and merge the results (e.g. Gonzalez et al. 2012; Dipierro et al. 2015b). Whilst these approaches neglect any feedback that the grain species have on the gas dynamics, they have proved a useful tool for direct comparison with observations.

From the perspective of dust dynamics, a long-term goal would be to model the dynamics of the whole grain population in discs simultaneously, in 3D, including the effects of the dust on the gas dynamics. Some progress towards this was made by Bai & Stone (2010a), Laibe & Price (2014c), showing how multiple grain species can be treated simultaneously within a one-fluid approach, but this is not yet implemented in any numerical code. Modelling the entire grain population would open the possibility of coupling the dust dynamics directly to the radiative transfer. In turn, the radiative transfer could then be used to set the gas temperature profile in the disc, allowing for thermodynamic feedback between the grain dynamics and the gas and the coupling to chemistry.

2.3. Radiative transfer

The transport of radiation through matter is important for three primary reasons. First, radiation can modify the composition and thermal properties of matter. For example, changing the composition and heating through mechanisms such as photoionisation and photodissociation and cooling it through the escape of line emission. Radiation can also set the dust temperature, which is determined by radiative equilibrium between thermal emission from the grains and the local radiation field [there are a number of textbooks with extensive discussion of these topics, such as Spitzer (1978), Rybicki & Lightman (1979), Osterbrock & Ferland (2006)]. This impact on the composition and thermal structure drives many macroscopic processes in discs (see e.g. Section 1, Figure 1). Second, radiation pressure can directly impart a force upon matter, altering the dynamics. Finally, radiation is what is actually observed. Radiative transfer is therefore required to make the most meaningful and robust comparisons between theoretical models and observations.

Since radiative transfer is fundamentally coupled to matter (influencing the composition and temperature, which in turn modifies opacities and emissivities), the coupling of radiation transport and chemistry is already an established field, which will be discussed further in Section 2.4.

For purely dynamical applications, the only quantities of interest from radiative transfer are a temperature/pressure estimate and/or a radiation pressure estimate. To this end, popular techniques are flux limited diffusion (FLD) and similar moment methods, owing to their relatively minimal computational expense compared with more detailed radiative transfer methods (e.g. Levermore & Pomraning 1981; Whitehouse & Bate 2004; Whitehouse, Bate, & Monaghan 2005). In FLD schemes, the directional properties of the radiation field are replaced by angle averaged ones and the radiative transfer problem is solved using a diffusion equation. FLD has long been applied in optically thick regimes without sharp density contrasts, but can generate spurious results where this is not the case (Owen, Ercolano, & Clarke 2014; Kuiper & Klessen 2013). Most modern applications of FLD account for this failure at low optical depth by using boundary conditions (e.g. Mayer et al. 2007), or using hybrid methods to allow the system to radiate energy away from optically thin regions (e.g. Boley et al. 2007; Forgan et al. 2009). Other approximate temperature prescriptions have also been developed that are tailored to model the effect of higher energy extreme ultraviolet (EUV) and X-ray photons from the parent star on the disc evolution (e.g. Alexander, Clarke, & Pringle 2006a, 2006b; Owen et al. 2010; Owen, Ercolano, & Clarke 2011; Owen, Clarke, & Ercolano 2012; Haworth, Clarke, & Owen 2016b).

More rigorous radiation transport methods have historically typically been confined to compute synthetic observables, where the density structure is based on snapshots from dynamical models, hydrostatic equilibrium in a simple disc, or a parametric model. Perhaps, the most popular method in this context is Monte Carlo radiative transfer (Lucy 1999), which is used by the well-known codes radmc-3d (Dullemond 2012), mcmax (Min et al. 2009), hyperion (Robitaille 2011), mcfost (Pinte et al. 2006), and torus (Harries 2015, also discussed below). Monte Carlo radiation transport typically involves breaking the energy from radiative sources into discrete packets, which are propagated through space in a random walk akin to the propagation of real photons through matter (e.g. including scattering and absorption/re-emission events). This provides an estimate of the mean intensity everywhere which can be used, for example, to solve for the ionisation state of a gas, the dust radiative equilibrium temperature, or to generate synthetic observations. The Monte Carlo approach naturally accounts for the processed radiation field (scatterings, recombination photons), works in arbitrarily geometrically complex media and also treats multi-frequency radiation transport (conversely FLD approaches typically assume that the opacity is frequency independent).

In addition to the Monte Carlo approach, other well-known methods are also the pure (e.g. Abel & Wandelt 2002) and short characteristic (e.g. Davis, Stone, & Jiang 2012) ray tracing schemes. Recently, intermediate expense hybrid-methods have been developed which combine FLD and other (e.g. ray-tracing) methods to offer a better balance between the accuracy of a more sophisticated scheme and the speed of FLD for dynamical applications (Kuiper & Klessen 2013; Owen et al. 2014; Ramsey & Dullemond 2015).

2.4. Chemistry

Molecular line observations play a central role in determining both the conditions within, and kinematics of, protoplanetary discs. In particular, CO and its isotopologues are popular tracers which are relatively abundant, have a permanent dipole moment and estimates of canonical abundances in the interstellar medium (ISM). CO synthetic observations can therefore be generated relatively easily in discs by assuming the canonical abundance and that local thermodynamic equilibrium (LTE) applies, in which case the level populations are set analytically by the Boltzmann distribution (e.g. Williams & Best 2014). However, such a simple approach is not always valid. For example, in discs there is evidence of departure from the canonical CO abundance (e.g. Favre et al. 2013) and the relative abundance of isotoplogues does not necessarily scale as in the ISM (Miotello, Bruderer, & van Dishoeck 2014b). Furthermore, dust grain evolution and dynamical processes such as instabilities and planet–disc interactions can also affect the chemistry (e.g. Boley et al. 2007; Ilee et al. 2011; Evans et al. 2015; Öberg et al. 2015a, 2015b; Cleeves, Bergin, & Harries 2015; Huang, Öberg, & Andrews 2016). Although simple CO parameterisations yield useful insights into the global properties of discs (such as the disc mass, e.g. Miotello et al. 2016; Williams & McPartland 2016), they are substantially more limited when it comes to probing the local properties. Given the above, more substantial chemical models will play an important role in the interpretation of modern protoplanetary disc observations. Furthermore, such models would support observations using molecules other than CO that are less easily parameterised, but could be better suited for probing certain components of a disc. In addition to interpret observations, understanding the chemical evolution of discs will also have astrobiological implications in the connection to the chemical composition of planets themselves.

To date, almost 200 molecules have been detected in interstellar or circumstellar environments ² . The abundances of these molecules can be subject to change via a large number of chemical reactions (see Caselli 2005; Henning & Semenov 2013, for reviews). In order to accurately model the evolution of even a small number of these molecules, complex computational networks of chemical reactions are needed in the form of coupled ordinary differential equations (ODEs). Several research groups have compiled publicly available databases of both these chemical reaction networks, and data on the rates of individual chemical reactions themselves—including the UMIST Database for Astrochemistry ³ (UDfA; Millar, Farquhar, & Willacy 1997; Woodall et al. 2007; McElroy et al. 2013), the Ohio State University networks ⁴ , and the Kinetic Database for Astrochemistry ⁵ (KIDA; Wakelam et al. 2012). Databases either contain these rates explicitly, or include how such a rate depends on local properties in the form of a parameterised expression (often via the Arrhenius–Kooij equation, Arrhenius 1889; Kooij 1893).

Chemical reactions fall into several categories and can involve a variety of reactants. Table 1 lists the common types of astrophysical reactions. Whilst the majority of reactions are concerned with gas-phase species or their interaction with photons, dust grain surfaces provide a location for further chemistry to occur. Gas-phase molecules attach themselves to the surfaces of dust grains (a process known as adsorption) via two mechanisms: physisorption (involving weak van der Waals forces) or chemisorption (due to chemical valence bonds). Once species are adsorbed, they produce layers of ices on the surface of dust grains, which allows more complex surface chemistry to occur (Herbst & van Dishoeck 2009). An example of this is the process of hydrogenation, by which hydrogen reacts quickly with other surface species (including itself) to produce saturated molecules such as methane. Of particular interest for this paper is that the composition of ices on dust grains (e.g. CO-coated versus H₂O-coated) can also affect the subsequent evolution of the dust by affecting the sticking efficiency and coagulation and fragmentation efficiencies [not discussed in detail here, but see e.g. Kouchi et al. (2002), Blum & Wurm (2008), Johansen et al. (2014); Musiolik et al. (2016), for further information]. Regions that are well shielded from incident stellar radiation (such as the disc midplane) might be thought to be chemically inert, as there is not sufficient energy to overcome reaction activation barriers. However, in such regions, ionisations caused by cosmic rays can induce ion–molecule reaction sequences that dominate much of the gas-phase chemistry, including the production of secondary cosmic-ray-induced photons. Increased densities in the disc midplane also mean that three-body reactions in the gas phase will begin to have an important effect on the chemistry. In these cases, a third body (M, the most abundant species, often molecular hydrogen) acts as a non-reacting catalyst.

Table 1. Common gas–grain reactions in astrophysical environments. Species are all considered to be in the gas phase, unless shown as X_gr, which are considered to be located on the ice mantles of dust grains. Photons are shown as γ and cosmic rays are shown as γ_cr. Adapted from Caselli (2005).

In addition to (closely coupled to) the computation of abundances is the computation of the temperature. This is determined by the heating and cooling rates, which are themselves set by, to name just a few: radiative processes (e.g. photoionisation heating and line cooling), dust/PAH’s (e.g. PAH heating and grain radiative cooling), chemical processes (e.g. exothermic reactions), hydrodynamic work/viscous heating, and cosmic rays (a review is given by Woitke 2015). Many of these heating/cooling terms are linked to the composition of the gas, requiring chemical and thermal calculations to be solved iteratively. In principle, since the heating and cooling is also set by the dust and radiation field, it might also be necessary to iterate over the (decoupled dust–gas) dynamics and radiative transfer.

Somewhat distinct from gas–grain chemistry are the photoionisation and photon-dominated region (PDR) regimes, where the radiation field plays a significant role in setting the composition and temperature of a medium. Photoionised gases are composed exclusively of atoms and ions and are typically modelled more in a radiative transfer context than a chemical one. Photoionisation models are usually concerned with the transfer of EUV photons and X-rays to solve for the ionisation balance and thermal structure of a gas of assumed gas and dust abundances. Despite not requiring chemical networks, this can include a variety of processes that are not trivially captured such as resonant line transfer and inner shell ionisations of atoms by X-rays (the liberation of multiple electrons by a single photon). Some examples of famous photoionisation codes are cloudy (Ferland et al. 2013) and mocassin (Ercolano et al. 2003). The photoionised regime only applies to disc winds, the very surface layers/inner edge of discs and, if the disc is externally irradiated by high energy photons (e.g. from a nearby O star), components of the flow from the disc outer edge.

The PDR regime applies at the transition between photoionisation and gas–grain dominated regimes; between predominantly ionised and molecular gasses. For example, in surface layers of the disc, but generally wherever matter is not optically thick to far ultraviolet (FUV) radiation. PDR modelling, like the gas–grain regime, requires a chemical network to be solved. It is also further complicated because cooling by line photons can be very important. This means that although the local radiation energy density (exciting the gas) is a single parameter, the escape probability of the line photons depends upon the extinction in all directions, i.e. it depends on the 3D structure of the surrounding space. Many PDR codes therefore compute this escape probability in one direction only, either working in 1D (e.g. models such as those in Röllig et al. 2007) or making some assumption about the dominant direction (e.g. vertically in the disc). Of the latter type, so called 1+1D models are particularly popular, which assume that at any given radial distance from the sta, the disc is in hydrostatic balance and escaping photons only consider the vertical distribution of gas at that radius (e.g. Gorti, Dullemond, & Hollenbach 2009; Woitke et al., 2016). Recently, multi-dimensional numerical approaches to solving PDR chemistry have appeared that do compute the 3D escape probabilities (Bisbas et al. 2012, 2015b) which they do efficiently using healpix (Górski et al. 2005).

2.4.1. Remarks on chemistry and radiative transfer

Chemical networks are used in conjunction with radiative transfer models to compute chemical abundances in various astrophysical environments. In general, the abundances are functions of temperature, density, and local radiation field, though many other parameters can play a role (in particular in the regime where line cooling is important, a measure of the extinction in all directions is ideally required). Often, the chemical networks are integrated to equilibrium in regions where the physical conditions are not thought to change significantly with time. However, in many cases, the microphysical conditions are functions of both space and time and are therefore not independent of dynamical processes (an example of this is given in Figure 3, see also Boley et al. 2007; Ilee et al. 2011; Evans et al. 2015; Drozdovskaya et al. 2016).

Figure 3. Left: The three-dimensional evolution of a tracer particle in a self-gravitating disc, colour coded with temperature changes, overlaid on the final column density snapshot of the disc. Right: The corresponding chemical evolution of particle, showing gas-phase CO and H₂CO, and CO ice (gCO). The shocks induced by the self-gravity of the disc have a significant impact on the chemical composition of the disc material (see Boley et al. 2007; Ilee et al. 2011; Evans et al. 2015).

Recent work has seen an increase in performing chemical evolution calculations alongside the radiative transfer calculations (e.g. Bruderer, Doty, & Benz 2009; Woitke, Kamp, & Thi 2009). Furthermore, chemistry is now being coupled directly with hydrodynamic calculations: In the context of star-forming regions, there are the full hydrodynamic models of Glover et al. (2010) and in a 1+1D disc framework, there are models such as those by Gorti et al. (2009) which also include radiative transfer. Such coupling is particularly important in the regions of the discs where the gas is not thermally coupled to the dust (i.e. in the upper layers of the disc, or within the dust sublimation radius), since the gas temperature, gas abundances, and level populations are strongly correlated. Unfortunately, it is in these regions of importance that 1+1D models become less applicable due to deviations from hydrostatic equilibrium (for example, thermally driven winds are not hydrostatic, e.g. Clarke & Alexander 2016). Dynamically, some chemical regimes (in particular, the PDR regime) are definitely important for understanding certain processes. For example, PDR physics is required to model FUV driven photoevaporative flows from the outer edge of discs (Adams et al. 2004; Facchini et al. 2016; Haworth et al. 2016a). The dynamical importance of gas–grain chemistry in cooler regions of the disc is currently yet to be determined, for example, presently unidentified chemically induced dynamical instabilities could potentially arise (see the burning questions, Section 3.1).

Aside from the coupling of chemistry with new physics such as dynamics, it is very important to stress that our base understanding of astrochemistry is constantly and rapidly evolving, with new species, reactions, and regimes being identified that can only be studied in a dedicated manner [for example, Penteado et al. (in preparation) use 10 000 models to study the sensitivity of single point chemical models to bind energies]. It is important that such focussed study continues.

Considering again the dust, there is no obvious consensus at present as to the best way to perform self-consistent dusty radiation hydrodynamics calculations of protoplanetary disc evolution. Schemes such as the short characteristics Variable Eddington Tensor (VET) method implemented in the Athena code by Davis et al. (2012), or the hybrid approach by Kuiper & Klessen (2013) show promise for bridging the gap between FLD and ray-tracing, but still require accurate modelling of the small grain dust population to determine the opacities before they can be applied in the context of protoplanetary discs (see Section 2.2).

With respect to magnetic fields, there are now also some approaches capable of modelling both radiation and MHDs (e.g. Flock et al. 2013; Tomida, Okuzumi, & Machida 2015).

Based on the above, we are already making excellent progress in cross-disciplinary modelling of discs, but most of this progress is very recent. There are still a number of highly coupled processes that cannot yet be modelled. As we will now discuss, there is a long, but fruitful journey ahead of multi-physics disc modellers.

3 BRIDGING THE GAPS—CHALLENGES

The interconnectedness of different processes in discs means that to be able to answer many of the outstanding theoretical and observational questions regarding protoplanetary discs, we will require a combination of 3D, global, multi-phase simulations with radiation hydrodynamics, dust dynamics, and size evolution, and chemistry computed self-consistently (see Figure 2).

3.1. Burning questions

Some examples of ‘burning’ science questions raised either during our discussion sessions, or by members of the community commenting on this paper, which might motivate improved multi-physics modelling of discs, included:

• • What are the main drivers of global disc evolution? In particular, what is the main driver of the mass accretion rate in protoplanetary discs?

• • Alongside magnetic fields, what other processes govern or control the launching of jets and outflows?

• • What is the effect of environment on protoplanetary disc evolution? For example, discs close to O stars are clearly heavily disrupted by high energy photons (we observe such systems as proplyds), but what is the role of comparatively modest radiation fields?

• • Do chemical–dynamical instabilities exist, i.e. is there a chemical reaction that feeds back into the dynamics (e.g. thermally) but responds to the dynamical change with a faster reaction rate?

• • What happens to small grains at the surface of the disc or in outflows/winds?

• • What happens at high dust to gas ratios? How important are streaming instabilities, or other instabilities? How important are self-induced dust traps? What happens to dust in shocks?

• • How do magnetic fields in the disc affect the behaviour of charged dust grains, and how do the dynamics and ionisation chemistry of the grain population in turn affect the magnetic field evolution?

• • What are the conditions under which pebble accretion (e.g. Ormel & Klahr 2010; Lambrechts & Johansen 2012; Morbidelli & Nesvorny 2012) might operate, and how will this impact the diversity of planetary systems formed in protoplanetary discs (e.g. Bitsch, Lambrechts, & Johansen 2015; Chambers 2016; Ida, Guillot, & Morbidelli 2016)?

• • What is the nature of fragmentation in self-gravitating discs? Is there a well-defined parameter space where fragmentation occurs (c.f. Meru & Bate 2011; Michael et al. 2012; Rice, Forgan & Armitage 2012; Rice et al. 2014), or can it occur stochastically through rare high-amplitude density perturbations over long enough timescales (Paardekooper 2012; Young & Clarke 2016)?

• • What is the origin of rings, gaps, horseshoes, and cavities observed in mm-continuum emission? How common are these features?

• • How can the masses and properties of embedded protoplanets be constrained from observations?

• • How do planets affect observations of chemical tracers?

• • How do planets and circumplanetary discs affect the evolution of the protoplanetary disc (e.g. through thermal feedback or increased radiative heating in gaps). Conversely, how does the disc affect an embedded planet (e.g. the planetary atmosphere).

• • Will dust discs fragment?

• • What determines the scale height of the dust layer? How is this set by different processes, for example, coagulation (e.g. Krijt & Ciesla 2016)?

• • Under which conditions do warps develop in discs? Can radiation pressure drive warping?

• • What are the possible initial conditions of class I/II/III discs and how do they influence the subsequent evolution? In particular, how does the early evolution of discs affect the chemistry and grain distribution (e.g. Miotello et al. 2014a)? What is inherited from the star-formation process?

• • The vertical component of the magnetic field controls the mass flux of winds and the saturation level of MRI-driven turbulence. How does the competition between accretion (drawing the vertical field in towards smaller radii) and diffusivity (pushing it outwards towards larger radii) cause this component of the field to vary with time? In particular, what is the magnitude of the diffusivity term, which is set by microphysics (e.g. Lubow, Papaloizou, & Pringle 1994; Rothstein & Lovelace 2008; Takeuchi & Okuzumi 2014)?

• • How turbulent are protoplanetary discs (e.g. Flaherty et al. 2015; Simon et al. 2015a; Teague et al. 2016)?

• • What is the process by which a protoplanetary disc becomes a debris disc? Transition discs; those with inner holes, are typically attributed to the action of photoevaporation by the host star (see e.g. Owen 2016), or planets (e.g. Zhu et al. 2011). But which, if either, of these is the dominant process [examples of models including both are Alexander & Armitage (2009), Rosotti, Ercolano, & Owen (2015)]? Are there other processes that contribute significantly to disc dispersal, such as magneto-thermal winds (Bai et al. 2016)? What are the initial conditions of debris disc models (e.g. Takeuchi, Clarke, & Lin 2005; Thilliez & Maddison 2015)?

Some of these questions might only be addressed by combining all of the physical ingredients of protoplanetary disc modelling. However, several will only require consideration of a smaller fraction. These smaller steps will be extremely valuable in bridging the gaps between fields, and will undoubtedly inform the production of a fully comprehensive modelling approach. We manifest these steps as a series of challenges, outlined below.

3.2. Grand challenges for gas modelling

C1: Model the pressure and temperature effects of photochemistry in multi-dimensional, fully hydrodynamic models

This challenges us to account for the (non-hydrostatic) dynamical impact of gas whose composition and temperature is set by photodissociation region processes. Specifically, the temperature should be accurately computed to within ~ 15% of a standard PDR network (which is the level of accuracy typically attained by reduced networks, see Section 5.1.1).

C2: Model the pressure and temperature effects of gas–grain chemistry in multi-dimensional, fully hydrodynamic models

Similar to challenge C1, this challenges us to account for the (non-hydrostatic) dynamical impact of chemical processes in optically thick regions of discs. There is a caveat to this challenge in that the dynamical importance of gas–grain chemistry is currently unknown. This therefore also (first) challenges us to determine what features of gas–grain chemistry might actually be dynamically important—such as chemically induced dynamical instabilities (see also the burning questions; Section 3.1).

C3: Incorporate the radiation field self-consistently whilst computing a multi-dimensional hydrodynamic model which satisfies challenges C1/C2

Challenges C1 and C2 are likely to be met by making simplifying assumptions about the incident radiation and cosmic ray background. The next step is then to properly account for the radiation field: Set by the central protostar, the disc material and any surrounding environment (e.g. the envelope or neighbouring stars/clouds/associations). This challenge will play a crucial role in understanding environmental influences on disc lifetimes.

C4: Model magnetic fields that can couple self-consistently to a realistic population of participating species

Models constructed to meet challenges C1–C3 that directly compute the composition of matter will deliver self-consistent populations of electrons, ions, and neutral species. The formation and evolution of magnetically active and dead zones, and the activation of MRI, is fundamental to the disc’s ability to accrete onto the star, as well as the launching of jets and outflows. We must therefore be able to couple the magnetic field evolution to the above gas–grain chemistry (see also challenge C9). Typically, MHD simulations that model the principal non-ideal processes (the Hall effect, Ohmic dissipation, and ambipolar diffusion) use simplified models for ion/grain mass and charge, often assuming single values for these properties. In practice, ion masses and charges will vary tremendously depending on the gas composition and the ambient radiation field.

In this challenge, non-ideal MHD models must be made flexible enough to accept arbitrary populations of a wide variety of ions (and grains, see C9) as input for computing subsequent magnetic field structure (c.f. the recent use of a reduced network by Tomida et al. 2015).

C5: Assemblage of gas modelling challenges

This essentially challenges us to model all components of the gas phase, i.e. to couple both C1 and C2, whilst incorporating C3 and C4. This challenge has two tiers. The lower tier involves accounting for all of the dynamical effects, without necessarily directly modelling the composition. Conversely, the higher tier does involve direct computation of the dynamically (and observationally) relevant chemical species.

3.3. Grand challenges for dust–gas modelling

Simultaneously compute the dynamics and size evolution of the entire grain population, coupled to self-consistent modelling of the gas and radiation field in the disc in global, 3D simulations. This can be broken into a series of smaller challenges, as follows:

C6: Model the dynamics of the entire grain population in a global disc simulation

Develop the means to accurately and efficiently model the dynamics of solids spanning an entire grain size distribution in global, 3D, disc simulations, including the effect of embedded companions and with feedback from the dust grains to the gas.

C7: Model the growth and fragmentation of solids

Develop an accurate prescription for growth and fragmentation of grains and incorporate it into 3D dynamical models of dust and gas evolution in global disc, with feedback from the dust grains to the gas.

C8: Radiative equilibrium and radiation pressure

Compute the radiative equilibrium temperature, as well as the radiation pressure force, in global 3D dynamical protoplanetary disc simulations, using multi-frequency radiative transfer.

C9: Coupling to MHD

Allow the dust–grain population, along with the radiation field, to determine the ionisation chemistry in the disc and use this to self-consistently model the development of jets, outflows, and MRI turbulence in both local and global disc models

C10: Assemblage of dust modelling challenges

Similar to C5, this challenges us to combine C6–C9. That is, to have a method of computing the motions of a whole grain distribution, including the evolution of grain sizes and the effects of radiation and magnetic fields.

C11: The grandest challenge (in this paper)

Develop a single model capable of reproducing multi-tracer, resolved, observations of a given protoplanetary disc. That is, perform a global disc simulation that solves for the gas and dust dynamics, as well as the dust and chemical evolution of the disc, that then predicts (to within a reasonable degree of accuracy) all observed properties of a given disc at a resolution comparable to that of current observational instrumentation. The model should retrieve the continuum morphology and intensity for wavelengths probing a range of grain sizes, whilst also reproducing molecular line observations of different tracers (for example, C¹⁸O, HCO⁺, ¹²CO, which probe different components of the disc and can be sensitive to different chemical effects).

Doing so will require simultaneous completion of many of the above challenges. It is therefore a long-term goal, but one which should be achievable given progress made on the other challenges stated above.

4 DISCUSSION—STRATEGIC STEPS TOWARDS THE FUTURE

The grand challenges discussed in the previous section are in a sense a strategic pathway towards long-term future development. In practice, models of discs are currently much more focussed, but could still be improved by the integration of previously uncoupled physics. In this section, we discuss broad strategy for the immediate future of more general disc modelling. More specific technical developments are discussed in the next section.

4.1. Which problems are the most pressing to solve and what physics is required to solve them?

It is inefficient to develop new software, or exhaust substantial CPU hours on an intensive state of the art multi-physics calculation, if the results have no value. A key strategic step, therefore, is to identify the combination of physics required to answer well motivated, well formulated, key problems.

Table 2 provides an example of a strategic overview. Such an overview can guide/motivate the development of numerical methods to include all of the physics essential to solve a given problem. It would also motivate us to understand whether the uncertain features really do play an important role.

Table 2. A qualitative summary of the effect of different components of disc modelling on the intrinsic physical properties of protoplanetary discs – ‘✓’ implies that an ingredient is identified as important, ‘?’ implies that the importance is uncertain, ‘✗’ implies that an ingredient is likely unimportant. It is our hope that such a summary would eventually become more quantitative, with the relative importance of different processes more formally assessed.

In addition to identifying the processes that might contribute to a problem (such as in Table 2), one could possibly then order the contributing physical processes in a hierarchy of importance to determine which are the most important features to include in a model (similar to the way that the dynamical importance of microphysics on H ii region expansion was categorised by Haworth et al. 2015). For example, consider the generation of synthetic molecular line observations. At the most basic level radiative transfer is required, as it is photons that are observed by astronomers, as well as some estimate of the density, temperature, molecular abundance, and molecular level populations. This can initially be done assuming some simple static disc structure, assuming an abundance of molecules and level populations determined analytically by the Boltzmann distribution. This could then be improved with proper non-local thermodynamic equilibrium (NLTE) statistical equilibrium calculations, which could be improved upon by using chemical networks/direct abundance calculations, which can be improved upon by further solving the dynamics/thermal balance, decoupled dust transport, and so on. In order to do this, one would first need to define some measure of importance. For example, if interested in accretion, a hierarchy of importance would place processes resulting in the largest contribution to the accretion rate at the top.

Deciding which problems are most pressing to address should also be informed by recent and upcoming observations. For example, which questions might be addressed by models in tandem with data from the Square Kilometer Array [SKA, which amongst other things will probe grain growth and disc chemistry (Testi et al. 2015)], James Webb Space Telescope (JWST) or the European-Extremely Large Telescope (E-ELT, e.g. Hippler et al. 2009)?

Another key strategic point is to stress that on the path towards multi-physics modelling of significantly interdependent physics, it is essential that all individual fields continue to develop as they are presently. Integration should be complementary to our current approaches. There are (at least) two key reasons for this. One reason is that an integrated approach is likely to be substantially more computationally expensive, which limits the parameter space of a given problem that can be studied. This also strongly limits the ability to quickly interpret observations (e.g. with parametric models). The other reason is that each field is continuing to evolve, with the development of new algorithms and the identification of new important mechanisms. This focussed field-by-field progression will likely answer a number of the burning questions and the techniques developed will ultimately feed back into multi-physics models of the future.

5 TECHNICAL STEPS TOWARDS THE FUTURE

We now discuss possible near-term developments of our numerical methods towards resolution of the grand challenges, focussing on the coupling of physical ingredients with a particular emphasis on chemistry.

5.1. Simplified chemistry for dynamics

We currently identify three possible approaches to include chemistry in dynamical simulations: direct calculation of a full network and heating/cooling rates, direct calculation of a reduced network, or implementation of pre-computed look-up tables. We discuss these further below.

5.1.1. Reduced chemical networks

Reduced chemical networks prioritise only the species and reactions of most importance to a given aim. For example, if prioritising dynamics, then an ideal reduced network would be one that yields a temperature/pressure to within an acceptable degree of accuracy (say 10–15%). The established method of generating a reduced network is to start with a comprehensive one and systematically remove components, checking that it does not have a substantial impact on the resulting quantity of interest. There are already codes available capable of computing chemistry based on very large networks, such as prodimo (Woitke et al. 2009), dali (Bruderer 2013), ucl-chem (Viti et al. 2004, 2011), ucl-pdr (Bell et al. 2005, 2006), and the models of Walsh et al. (2012). Any of these networks could be analysed to determine which processes are essential for dynamics, and then reduced accordingly. Additionally, it is also possible to optimise calculations of large networks (e.g. Grassi et al. 2013). It is likely that a combined approach of reduction and optimisation will yield the most efficient results.

PDR chemistry is important in surface layers and the disc outer edge if the disc is externally irradiated. Reduced PDR networks already exist (e.g. Röllig et al. 2007). Such a network is already used in dynamical models by torus-3dpdr (see Section 5.2). However, existing reduced PDR networks are predominantly motivated by studies of star-forming regions/the ISM. New reduced networks tailored for discs would be extremely valuable for future dynamical models including PDR chemistry.

In the same vein as reduced chemical networks, there are also some recent promising developments concerning the relatively computationally cheap determination of the ionisation state in dense, dusty, optically thick regions of discs (in particular where dust-phase recombination dominates over the gas phase) which is particularly important for MHD and evolution of the dust population (e.g. regarding coagulation). Ivlev, Akimkin, & Caselli (2016) present an analytic model that yields the ionisation state of such dusty media, which could be incorporated into non-ideal MHD codes—offering an imminently achievable significant advance.

5.1.2. Lookup tables and functional parameterisations

An alternative to direct computation of the chemistry/temperature using a ‘full’ or reduced network is to tabulate temperatures or heating/cooling rates as a function of local properties in a disc. For example, Owen et al. (2010) prescribe the temperature of gas optically thin to X-rays as a function of local ionisation parameter (i.e. the density, distance from the source, and stellar X-ray luminosity) where the function (itself only published in full in Haworth et al. 2016b) was computed by the dedicated photoionisation code mocassin (Ercolano et al. 2003, 2008). A similar approach to obtaining PDR or gas–grain chemistry temperatures, where lookup tables are computed prior to run-time, could vastly reduce the potential computational expense of dynamical models.

Unfortunately, chemistry (both gas–grain and PDR) is not generally so easily parameterised as a simple function of the local properties. In order to include all relevant effects of heating and cooling, such a look-up table could easily grow very large. Below we briefly list several example quantities that would need to be included, along with a typical dimensionality for each in parenthesis (I. Kamp, private communication):

• • The temperature of dust grains (1).

• • The dust grain size(s), including second moment of the size distribution for grain surface chemistry and collisional gas–grain coupling (2).

• • The dust grain density (1).

• • The gas density (1).

• • Column densities towards the central star of key species (H, C, CO) for evaluating the amount of shielding (3).

• • The cosmic ray ionisation rate (though this can perhaps be approximated as constant throughout the disc) (1).

• • The strength of the radiation field in several bands, including X-Ray, UV, and optical (10).

• • The optical depth of the dust in direction of closest escape (1).

• • The fractional abundance of polycyclic aromatic hydrocarbons (PAHs) and further details if also using PAHs as an opacity source (3). Again, these parameters may be constants throughout the disc.

• • Column densities of all species to be considered, both towards the radiation source, and the direction of closest escape (≳ 10).

Given that the above list is by no means exhaustive, it is easy to see that such a look-up table may reach a dimensionality of 30–40. One of the key factors accounting for this issue is that the local chemistry depends upon the 3D non-local density distribution, because this sets photon escape probabilities, i.e. the chemistry at some point in space cares about the gas distribution in all directions from that point. It is therefore not solely dependent upon local properties, even if the local radiation field from external sources has been computed.

However, many of these quantities are likely not entirely independent, and relations between them could be identified in a statistically robust manner using grids of simulations. This may allow a reduction in the number of dimensions required. Of further note is that a ‘simplified’ thermodynamic prescription based on chemical modelling was developed by Woitke, Krueger, & Sedlmayr (1996a), Woitke, Goeres, & Sedlmayr (1996b), Schirrmacher, Woitke, & Sedlmayr (2003) for application to pulsating stars, which might offer some guidance on how to streamline some of the aforementioned dependencies.

5.2. Direct hybridisation

Historically, the approach to include more physics in dynamical models is to use a hydrodynamical code as the foundation and incorporate simplified physics modules subsequently. For example, Glover et al. (2010) and Dzyurkevich et al. (2016) patch reduced chemical networks into zeus-mp and ramses, respectively. Flock et al. (2013) also present an extension of the pluto code that includes both magnetic fields and an FLD radiation transport scheme. There is another approach, which is to start with a state of the art chemistry/radiative transfer code and subsequently incorporate somewhat more simple hydrodynamics. An example of this latter approach is the torus radiation transport and hydrodynamics code. This code began its life solely as a Monte Carlo radiative transfer code (Harries 2000) but now includes hydrodynamics, so can perform radiation hydrodynamic simulations with all the features of a dedicated radiation transport code (e.g. detailed photoionisation, dust radiative equilibrium, and radiation pressure in arbitrarily complex system geometries, etc.; Haworth & Harries 2012; Harries 2015; Haworth et al. 2015). Furthermore, torus-3dpdr is an extension of torus that also includes PDR chemistry with 3D extinction and escape probabilities (Bisbas et al. 2015b). The UV radiation field everywhere is computed using the Monte Carlo radiation transport and the escape probabilities are estimated in 3D using an algorithm based on healpix (Górski et al. 2005). torus-3dpdr is capable of directly modelling the role of FUV photons dynamically in non-hydrostatic scenarios, such as the external irradiation of discs by FUV radiation that has only been possible semi-analytically in the past (Adams et al. 2004; Facchini et al. 2016; Haworth et al. 2016a). It could also be used to test the validitiy of escape probability methods that assume a single dominant trajectory (the 1+1D methods).

One argument in favour of adding hydrodynamics to a radiative transfer/chemistry code is development time, since a simple but effective hydrodynamics algorithm is usually much more straightforward to develop than a radiative transfer/chemistry algorithm (though of course care must be taken to ensure that the hydrodynamics algorithm is appropriate for any given application). The obvious argument against this coupling of state of the art physics models with hydrodynamics is that they are not necessarily well streamlined and can be very computationally expensive (though this is not necessarily a problem if the code is optimised and/or highly scalable, as is the case for torus, Harries 2015).

Constructing a dedicated photochemical-dynamical code from scratch is another possible option, but potentially requires a lot of development time (e.g. the recent PDR-dynamical code of Motoyama et al. 2015).

Another promising avenue is the development of diverse, flexible self-consistent physics libraries that can be ported into other numerical (and therefore potentially hydrodynamical) codes. The krome code is an excellent example of this approach, which quickly solves arbitrary chemical networks and can also calculate heating and cooling terms (Grassi et al. 2014). Spectral codes, which solve partial differential equations flexibly and efficiently, could also offer a powerful means of combining other physical ingredients in a relatively straightforward manner. Spectral codes appear not to have featured in multi-physics disc modelling to date, but options for doing so include the dedalus (Burns et al. 2016) and snoopy (Lesur 2015) codes.

5.3. Temporal and spatial resolution

A very specific problem is that (in particular for non-equilibrium chemical-dynamics) we have to determine what the spatial and temporal scales are that we have to resolve in a given scenario. As an example, chemical timescales in the disc upper layers (that is, in the PDR regions) are rather short, whereas timescales deeper in the disc are usually much longer (for example, the case of CO being converted into CH₄ on timescales even longer than protoplanetary discs lifetimes). The time steps required to model the upper layers may therefore eventually be limited by the chemical timescales (in non-equilibrium scenarios) rather than the dynamical timescales, which might drastically increase computational expense. In such a regime where the chemical timescale is very small (much smaller than the dynamical timescale), we may be able to alleviate the problem somewhat with chemical sub-stepping—running multiple chemical updates per hydrodynamic update. Conversely, if the chemical/thermal timescales (reaction/heating/cooling rates) are very long, many dynamical steps can be taken between the more expensive chemical updates, improving the computation time substantially.

Alternatively, if the system is expected to reach a steady state, and all that is desired is an accurate model of this steady state (rather than the pathway to reaching the steady state), it may be possible to run chemical calculations very infrequently even if the chemical timescale is very short.

In addition to the above timescale arguments, resolution also needs to be considered. For example, some chemical features may only arise if the spatial resolution (e.g. around shocks) is sufficiently high—capturing such processes will of course increase computational expense.

5.4. Scaling

A key technical consideration is the scaling of the various physical ingredients in terms of both elements (cells, rays, chemical species, reactions, etc.) and computational resources (number of cores), since a calculation is going to be limited by its least tractable component. Different numerical approaches to computing a given ingredient can scale very differently. For example, in the case of radiative transfer, Monte Carlo radiation transport and treecol (see Clark, Glover, & Klessen 2012, for details of the latter) scale much more efficiently than long characteristic ray tracing. There are therefore multiple scaling options per ingredient.

For applications comprising two or more ingredients that scale very differently, there will likely be idle cores/inefficient CPU usage in the components of the code that do not scale so well. Furthermore, some techniques have specific constraints on the number/configuration of cores which may vary for different calculation ingredients. For example, if the hydrodynamic component of a calculation were confined to i distributed memory MPI threads (plus an arbitrary number of shared memory openMP threads), but the radiative transfer to j > i MPI threads, there will be unused MPI threads during each hydrodynamics step. This is a situation where dynamically optimising between shared and distributed memory processes is worthwhile, setting the otherwise idle MPI threads to contribute to openMP or other useful tasks.

5.5. Hardware developments

It is also important to assess new and projected hardware developments. We are approaching a time in which access to large numbers of processors increasingly outweighs the developments in performance of the processors themselves. Efficiently scalable numerical methods, such as Monte Carlo radiation transport and discontinuous Galerkin hydrodynamics solvers, will therefore be extremely advantageous in the near future.

Another significant realisation (only recently for astronomers) is that graphics processing units (GPUs) can offer significant speedup per core. A relatively small (but growing) fraction of astrophysical codes have a GPU implementation, and those that do are often those used for cosmological applications (e.g. Schive, Tsai, & Chiueh 2010; Bryan et al. 2014b). However, a GPU implementation of the fargo disc-modelling code was developed by Benítez-Llambay & Masset (2016), where they quote a typical speedup per core of a factor 40. It is beyond the scope of this paper to discuss GPU programming in detail, but we note that GPUs are fundamentally different architectures to CPUs and are therefore programmed in a somewhat different manner (taking time to learn), typically using either the cuda (Nickolls et al. 2008) or opencl (Stone, Gohara, & Shi 2010) standards. The high speeds of GPUs make them a powerful tool for the future of astronomy, where applicable, and they are likely to feature much more frequently in astronomy in the coming years, especially with the advent of directive-based GPU acceleration using the OpenACC standard ⁶ .

A final example, mentioned here only in passing, is the introduction of new types of processor such as the Intel Xeon Phi (e.g. Jeffers & Reinders 2013)—which combines some of the performance advantages of GPUs with an easier programming framework.

In general, the writing of new codes, or adapting old ones, to take advantage of hardware developments will be important. Given that more specialised hardware might continue to appear over time, it would also be advantageous if codes could be developed in such a fashion that they are easily ported, but it is unclear (to us at least) exactly how this might work in practice. This is an area where increased collaboration between astrophysicists and computer scientists will be advantageous. Interaction with computer scientists could also lead to other benefits such as improved efficiency of our codes and the promotion of better coding practice.

6 COLLABORATIVE STEPS TOWARDS THE FUTURE

As already mentioned, the components that we want to couple in the future of disc modelling are themselves already established and complex fields. It is therefore clear that these challenges are a whole-community effort, and substantial progress will only be made via collaboration. To this end, we have identified several key collaborative steps that we discuss below.

6.1. Workshops

Workshops are likely to be essential for stimulating cross-disciplinary collaboration. Whilst a typical conference setting will be important for each sub-discipline to discuss their work generally, events with ample time for break-out sessions and collaborative spaces are likely to be very productive. Such events allow large-scale discussion, but also allow for specific problems to be tackled one-on-one or in small groups in an ‘unconference’ setting (for example, the dotAstronomy ⁷ or Astropy ⁸ conference series). The identification of key ingredients to be swapped between respective fields will be important to establish, e.g. heating and cooling rates are likely to be of interest to those running dynamic models, whilst detailed abundance results may not be required.

6.2. Benchmarking

In addition to workshops, it is important for each field to develop an agreed set of benchmark problems, with the aim of transparency and reproducibility. Code comparison projects are key, but can require a lot of work for a small number of publications (albeit high impact, e.g. de Val-Borro et al. 2006; Röllig et al. 2007; Pinte et al. 2009; Iliev et al. 2009).

A good example of a successful comparison project is the recent StarBench code comparison workshops ⁹ (Bisbas et al. 2015a). These workshops aimed at testing radiation hydrodynamics codes used to study problems in star formation, with an emphasis on doing so in a positive and friendly environment. The workshops involved attendees running tests before arrival, which spanned a range of complexity. In the first meeting at the University of Exeter in 2013, every code passed the purely hydrodynamic shock tests without issue. However, the instant that radiative transfer/photoionisation was introduced into the dynamical problem we generally had poor agreement, even for the simplest case of tracking the time evolution of the extent of an ionised region about a star in a uniform density medium composed solely of hydrogen. The origin of the inconsistency between codes was that they were all running slightly different models (e.g. inconsistent recombination rates) and, after extremely careful rewriting of the specifications of this simple test, were subsequently able to get truly excellent agreement between the codes. This process highlighted to the community all of the things that should be explicitly stated in a paper in order to make it truly reproducible. Last but not least, in the case of an expanding H ii region, we actually discovered that although the codes all agreed perfectly, they did not agree with the classic analytic solution that everyone would compare with in their numerical methods paper and suggest that they get ‘good enough’ agreement with—validating their approach. Following re-investigation, as a result of code comparison, a direct improvement in our understanding of this fundamental analytic problem has been established (Bisbas et al. 2015a). In summary, code comparison

• • verifies that codes are working as desired;

• • informs the community what needs to be specified in papers to make them reproducible—a key factor, especially since there are likely to be many more ingredients in disc models of the future than there were in the relatively straightforward StarBench tests;

• • improves our understanding of each other’s numerical methods, including relative strengths and weaknesses. This can be done in a friendly way;

• • highlights the importance of careful numerics (e.g. understanding resolution dependency and which techniques are appropriate for a given scenario);

• • results in high impact publications;

• • leads to an improvement in our understanding of the underlying more fundamental (even analytic) problems.

Key to a successful comparison is active feedback between participants and iteration towards understanding the origin of any differences encountered. This can often be achieved just as easily with a comparison involving just two or three codes performed by a relatively small team (e.g. Bate & Burkert 1997; Commerçon et al. 2008; Price & Federrath 2010; Hubber, Falle, & Goodwin 2013). Such an approach avoids much of the friction associated with large-scale comparison projects whilst achieving the same objectives.

6.3. Open source software

A more applied collaborative practice is to develop software in an open source format (e.g. using GitHub ¹⁰ ). This is potentially very useful for both transparency and distributed development (i.e. international contributors). Examples taking such an approach are the krome (Grassi et al. 2014) and lime (Brinch & Hogerheijde 2010) projects.

Although the open source mentality is desirable, it should not be imposed since there may be legitimate reasons to protect intellectual interests. For example, if an early-career researcher invests substantial time into code development, the current academic culture requires a period where they are able to see a return on their time investment, in terms of first author publications where they lead astrophysical research (in the current culture, this is more important than a number of co-authored publications). There is no reason that their code cannot be shared collaboratively during such a phase of research. More widespread access can subsequently be yielded once the developers have seen sufficient return.

7 SUMMARY

Protoplanetary discs are a key focus of modern astronomy, being subject to extensive modelling including the dynamics of gas and dust, magnetic fields, radiation transport, and chemistry. These facets of physics required to model discs are, however, not independent, so as we proceed into the future we must consider their coupling in multi-physics modelling of discs. In particular, we perceive that it will be important to self-consistently model decoupled gas and dust dynamics, with radiative transfer, dust growth/fragmentation, and different chemical regimes (gas–grain, PDR). This paper aims to stimulate this development and consisted of the following components.

First, to establish a platform from which to discuss the coupling of different disciplines, we provide an overview of each in isolation, as well as the progress made towards multi-physics modelling to date. Using this, we have identified a series of challenges for the future of protoplanetary disc modelling, which are supposed to act as milestones towards the ultimate goal of a self-consistent gas, dust, radiation transport, and chemistry model mentioned above. Our first category of challenges regards gas modelling, with a particular focus on composition (e.g. gas–grain and photochemistry) coupled with dynamics. Our second category of challenges regards dust, including modelling of an entire grain size distribution as well as growth and fragmentation of grains and any additional physics (such as radiation) that alters the dust dynamics. We also discuss pathways towards addressing these challenges, which are grouped by whether they are strategic (e.g. identifying what needs to be done), technical (e.g. working out how to do it) and collaborative (working together to do it).

Table 3. Participants of protoplanetary discussions 2016.

We finish by noting that, as a further motivational strategy, appropriate agents mights offer prize(s) for completing more rigorously defined versions of one or more of the challenges presented here.