Hostname: page-component-848d4c4894-cjp7w Total loading time: 0 Render date: 2024-06-19T21:09:24.601Z Has data issue: false hasContentIssue false

Grand Challenges in Protoplanetary Disc Modelling

Published online by Cambridge University Press:  21 October 2016

Thomas J. Haworth*
Institute of Astronomy, Madingley Road, Cambridge, CB3 0HA, UK Astrophysics Group, Imperial College London, Blackett Laboratory, Prince Consort Road, London SW7 2AZ, UK
John D. Ilee
Institute of Astronomy, Madingley Road, Cambridge, CB3 0HA, UK
Duncan H. Forgan
School of Physics and Astronomy, University of St Andrews, North Haugh, St Andrews KY16 9SS, UK
Stefano Facchini
Max-Planck-Institut für Extraterrestrische Physik, Giessenbachstrasse 1, 85748 Garching, Germany
Daniel J. Price
Monash Centre for Astrophysics and School of Physics and Astronomy, Monash University, Clayton, VIC 3800, Australia
Dominika M. Boneberg
Institute of Astronomy, Madingley Road, Cambridge, CB3 0HA, UK
Richard A. Booth
Institute of Astronomy, Madingley Road, Cambridge, CB3 0HA, UK
Cathie J. Clarke
Institute of Astronomy, Madingley Road, Cambridge, CB3 0HA, UK
Jean-François Gonzalez
Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230, Saint-Genis-Laval, France
Mark A. Hutchison
Centre for Astrophysics and Supercomputing, Swinburne University of Technology, Hawthorn, VIC 3122, Australia
Inga Kamp
Kapteyn Astronomical Institute, Postbus 800, 9700 AV Groningen, The Netherlands
Guillaume Laibe
School of Physics and Astronomy, University of St Andrews, North Haugh, St Andrews KY16 9SS, UK
Wladimir Lyra
Department of Physics and Astronomy, California State University Northridge, 18111 Nordhoff St, Northridge, CA 91330, USA Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA, 91109, USA
Farzana Meru
Institute of Astronomy, Madingley Road, Cambridge, CB3 0HA, UK
Subhanjoy Mohanty
Astrophysics Group, Imperial College London, Blackett Laboratory, Prince Consort Road, London SW7 2AZ, UK
Olja Panić
Institute of Astronomy, Madingley Road, Cambridge, CB3 0HA, UK School of Physics and Astronomy, University of Leeds, Leeds, LS2 9JT, UK
Ken Rice
Institute for Astronomy, University of Edinburgh, Blackford Hill, Edinburgh EH9 3HJ, UK
Takeru Suzuki
School of Arts & Sciences, University of Tokyo 3-8-1, Komaba, Meguro, Tokyo, 153-8902, Japan Department of Physics, Nagoya University, Furo-cho, Chikusa, Nagoya, Aichi, 464-8602, Japan
Richard Teague
Max-Planck-Institut fr Astronomie, Knigstuhl 17, 69117 Heidelberg, Germany
Catherine Walsh
Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA, Leiden, The Netherlands
Peter Woitke
Astrophysics Group, Imperial College London, Blackett Laboratory, Prince Consort Road, London SW7 2AZ, UK
Rights & Permissions [Opens in a new window]


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.

Research Article
Copyright © Astronomical Society of Australia 2016 


For the first time in history, spatially resolved observations of the structures within protoplanetary discs are being obtained (see review by Casassus Reference Casassus2016). This has revealed a wealth of sub-structure, including rings and gaps (ALMA Partnership et al. 2015; Andrews et al. Reference Andrews2016; Canovas et al. Reference Canovas, Caceres, Schreiber, Hardy, Cieza, Ménard and Hales2016), spirals (e.g. Garufi et al. Reference Garufi2013; Benisty et al. Reference Benisty2015; Wagner et al. Reference Wagner, Apai, Kasper and Robberto2015), warps (e.g. Casassus et al. Reference Casassus2015), shadows (e.g. Stolker et al. Reference Stolker2016), cavities (e.g. Andrews et al. Reference Andrews, Wilner, Espaillat, Hughes, Dullemond, McClure, Qi and Brown2011), and dust traps (e.g. van der Marel et al. Reference van der Marel2013, Reference van der Marel, van Dishoeck, Bruderer, Andrews, Pontoppidan, Herczeg, van Kempen and Miotello2016). These recent observations, combined with the huge diversity of exoplanetary systems discovered over recent years (Winn & Fabrycky Reference Winn and Fabrycky2015), 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 Reference Papaloizou and Terquem2006).

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 Reference Gorti and Hollenbach2009; Thiabaud et al. Reference Thiabaud, Marboeuf, Alibert, Leya and Mezger2015; Walsh, Nomura, & van Dishoeck Reference Walsh, Nomura and van Dishoeck2015; Gorti, Hollenbach, & Dullemond Reference Gorti, Hollenbach and Dullemond2015). 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. Reference Meijer, Dominik, de Koter, Dullemond and Waters2008; Woitke et al. Reference Woitke2016; Miotello et al. Reference Miotello, van Dishoeck, Kama and Bruderer2016; Boneberg et al. Reference Boneberg, Panić, Haworth, Clarke and Min2016; Kama et al. Reference Kama2016).

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 Reference Bodenheimer1995; Dullemond et al. Reference Dullemond, Hollenbach, Kamp and D’Alessio2007; Lodato Reference Lodato2008; Armitage Reference Armitage2011, Reference Armitage2015). 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. Reference Testi2014). 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. Reference Durisen, Boss, Mayer, Nelson, Quinn and Rice2007; Young & Clarke Reference Young and Clarke2015; Forgan, Parker, & Rice Reference Forgan, Parker and Rice2015; Meru Reference Meru2015; Takahashi, Tsukamoto, & Inutsuka Reference Takahashi, Tsukamoto and Inutsuka2016), the streaming instability (Youdin & Goodman Reference Youdin and Goodman2005), Rossby wave instability (e.g. Lovelace et al. Reference Lovelace, Li, Colgate and Nelson1999; Tagger Reference Tagger2001; Lyra et al. Reference Lyra, Johansen, Klahr and Piskunov2008b, Reference Lyra, Johansen, Klahr and Piskunov2009), baroclinic and vertical shear instabilities, which can form and grow vortex structures (Lyra & Klahr Reference Lyra and Klahr2011; Lesur & Papaloizou Reference Lesur and Papaloizou2010; Nelson, Gressel, & Umurhan Reference Nelson, Gressel and Umurhan2013; Richard, Nelson, & Umurhan Reference Richard, Nelson and Umurhan2016), the magneto-rotational instability (e.g. Balbus & Hawley Reference Balbus and Hawley1991; Reyes-Ruiz et al. Reference Reyes-Ruiz, Pérez-Tijerina, Sánchez-Salcedo, Reyes-Ruiz and Vázquez-Semadeni2003), and dust-settling induced vortices (Lorén-Aguilar & Bate Reference Lorén-Aguilar and Bate2015, Reference Lorén-Aguilar and Bate2016). The local environment can also significantly modify disc evolution via mass transfer from the ambient medium onto the disc (Vorobyov, Lin, & Guedel Reference Vorobyov, Lin and Guedel2015; Lomax, Whitworth, & Hubber Reference Lomax, Whitworth and Hubber2015), nearby radiation sources (e.g. Bally, O’Dell, & McCaughrean Reference Bally, O’Dell and McCaughrean2000; Henney et al. Reference Henney, O’Dell, Meaburn, Garrington and Lopez2002; Smith, Bally, & Morse Reference Smith, Bally and Morse2003; Adams et al. Reference Adams, Hollenbach, Laughlin and Gorti2004; Wright et al. Reference Wright, Drake, Drew, Guarcello, Gutermuth, Hora and Kraemer2012; Facchini, Clarke, & Bisbas Reference Facchini, Clarke and Bisbas2016) and tidal encounters (e.g. Clarke & Pringle Reference Clarke and Pringle1993; de Juan Ovelar et al. Reference de Juan Ovelar, Kruijssen, Bressert, Testi, Bastian and Cánovas2012; Rosotti et al. Reference Rosotti, Dale, Ovelar, Hubber, Kruijssen, Ercolano and Walch2014; Vincke, Breslau, & Pfalzner Reference Vincke, Breslau and Pfalzner2015; Dai et al. Reference Dai, Facchini, Clarke and Haworth2015; Vincke & Pfalzner Reference Vincke and Pfalzner2016). 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 Reference Armitage, Livio and Pringle2001), but is only now being investigated with self-consistent hydrodynamic simulations (e.g. Bae et al. Reference Bae, Hartmann, Zhu and Nelson2014). 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. Reference Hutchison, Price, Laibe and Maddison2016) which in turn is sensitive to dynamical effects such as shadowing caused by warping of the inner disc (Marino, Perez, & Casassus Reference Marino, Perez and Casassus2015; Stolker et al. Reference Stolker2016). 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. Reference Testi2014). As a final example, gravitational instability and fragmentation in discs is sensitive to radiation (e.g. Meru & Bate Reference Meru and Bate2010; Forgan & Rice Reference Forgan and Rice2013) and magnetic fields (Price & Bate Reference Price and Bate2007; Wurster, Price, & Bate Reference Wurster, Price and Bate2016), and can induce dramatic effects in the chemical composition of discs (see Section 3, Ilee et al. Reference Ilee, Boley, Caselli, Durisen, Hartquist and Rawlings2011; Evans et al. Reference Evans, Ilee, Boley, Caselli, Durisen, Hartquist and Rawlings2015).

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 Reference Hartmann1998; Armitage Reference Armitage2011; Williams & Cieza Reference Williams and Cieza2011; Armitage Reference Armitage2015). 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 EdinburghFootnote 1 , 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 Reference Williams and Best2014). 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 46, we discuss pathways towards meeting the challenges in terms of strategical, technical, and collaborative developments.


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 Reference Masset2000), implemented in both the fargo (Masset Reference Masset2000; Baruteau & Masset Reference Baruteau and Masset2008; Benítez-Llambay & Masset Reference Benítez-Llambay and Masset2016) and pluto (Mignone et al. Reference Mignone, Bodo, Massaglia, Matsakos, Tesileanu, Zanni and Ferrari2007) 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 Reference Berger and Colella1989), which is mainly (but not exclusively) developed for Cartesian meshes, is not typically used (an example of an exception is Paardekooper & Mellema Reference Paardekooper and Mellema2004). Furthermore, simulating warped, twisted, or broken discs remains difficult (e.g. Fragner & Nelson Reference Fragner and Nelson2010).

Lagrangian schemes such as smoothed particle hydrodynamics (SPH, for reviews, see e.g. Monaghan Reference Monaghan1992; Price Reference Price2012) 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 (Reference Ogilvie1999) of α-disc theory (Lodato & Price Reference Lodato and Price2010). 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. Reference Nixon, King, Price and Frank2012; Nixon, King, & Price Reference Nixon, King and Price2013; Nealon, Price, & Nixon Reference Nealon, Price and Nixon2015) or breaks (Nixon & King Reference Nixon and King2012; Facchini, Lodato, & Price Reference Facchini, Lodato and Price2013; Doğan et al. Reference Doğan, Nixon, King and Price2015). Such behaviour is well modelled by SPH codes, and appears to be relevant to observed protoplanetary discs, including HK Tau (Stapelfeldt et al. Reference Stapelfeldt, Krist, Ménard, Bouvier, Padgett and Burrows1998), KH15D (Lodato & Facchini Reference Lodato and Facchini2013), and HD142527 (Casassus et al. Reference Casassus2015). 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. Reference de Val-Borro2006).

2.1.2. Local simulations

The most common technique utilised for local simulations of discs is the Cartesian shearing box (Hawley, Gammie, & Balbus Reference Hawley, Gammie and Balbus1995; Guan & Gammie Reference Guan and Gammie2008). 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 Reference Toro2013). This approach has been used almost exclusively for simulating the magnetorotational (see Balbus Reference Balbus2003, and references within) and other instabilities—in particular the streaming instability (e.g. Youdin & Goodman Reference Youdin and Goodman2005; Youdin & Johansen Reference Youdin and Johansen2007; Johansen et al. Reference Johansen, Oishi, Mac Low, Klahr, Henning and Youdin2007; Bai & Stone Reference Bai and Stone2010b)—in discs. Though other applications include the study of magnetically driven disc winds (e.g. Suzuki & Inutsuka Reference Suzuki and Inutsuka2009; Suzuki, Muto, & Inutsuka Reference Suzuki, Muto and Inutsuka2010).

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. Reference Tasker, Brunino, Mitchell, Michielsen, Hopton, Pearce, Bryan and Theuns2008; Price & Federrath Reference Price and Federrath2010), 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 Reference Springel2010, Reference Springel2011; Duffell & MacFadyen Reference Duffell and MacFadyen2011) and meshless codes (Maron, McNally, & Mac Low Reference Maron, McNally and Mac Low2012; McNally, Maron, & Mac Low Reference McNally, Maron and Mac Low2012; Hopkins Reference Hopkins2015). 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. Reference Mocz, Vogelsberger, Pakmor, Genel, Springel and Hernquist2015) although early applications to protostellar discs appear to be promising (see e.g. Muñoz et al. Reference Muñoz, Kratter, Springel and Hernquist2014).

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. Reference Schaal, Bauer, Chandrashekar, Pakmor, Klingenberg and Springel2015) and moving Voronoi mesh (e.g. the arepo code; Mocz et al. Reference Mocz, Vogelsberger, Sijacki, Pakmor and Hernquist2014) 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. Reference Stone, Gardiner, Teuben, Hawley and Simon2008), enzo (Bryan et al. Reference Bryan2014a), fargo (Benítez-Llambay & Masset Reference Benítez-Llambay and Masset2016), pluto (Mignone et al. Reference Mignone, Bodo, Massaglia, Matsakos, Tesileanu, Zanni and Ferrari2007), and pencil (Brandenburg & Dobler Reference Brandenburg and Dobler2002) codes. SPH and other meshless codes can now also incorporate MHD (see review by Price Reference Price2012), provided that the ∇ · B = 0 condition can be sustained, for example, using divergence cleaning techniques (Tricco & Price Reference Tricco and Price2012). 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 Reference Lewis, Bate and Tricco2016).

Whilst ideal MHD disc simulations have been conducted for some time (see Balbus Reference Balbus2003, 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 (Reference Gammie1996). More recently, the interplay between the Hall effect, ambipolar diffusion, and Ohmic diffusion is yielding new turbulent behaviour (Sano & Stone Reference Sano and Stone2002; Simon et al. Reference Simon, Lesur, Kunz and Armitage2015b), new forms of instability, and zonal flows in both MRI-active and ‘dead zone’ regions (e.g. Kunz & Lesur Reference Kunz and Lesur2013; Bai & Stone Reference Bai and Stone2014), not to mention addressing the so-called magnetic braking catastrophe that suppresses disc formation in ideal MHD (Tsukamoto et al. Reference Tsukamoto, Iwasaki, Okuzumi, Machida and Inutsuka2015; Wurster et al. Reference Wurster, Price and Bate2016) (see recent review by Tsukamoto Reference Tsukamoto2016, 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 Reference Casse, Meliani and Sauty2007; Bai Reference Bai2014; Lovelace et al. Reference Lovelace, Romanova, Lii and Dyda2014; Suzuki & Inutsuka Reference Suzuki and Inutsuka2014; Staff et al. Reference Staff, Koning, Ouyed, Tanaka and Tan2016).

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. Reference Dipierro, Price, Laibe, Hirsh, Cerioli and Lodato2015b; Jin et al. Reference Jin, Li, Isella, Li and Ji2016; Rosotti et al. Reference Rosotti, Juhasz, Booth and Clarke2016).

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 Reference Monaghan and Kocharyan1995; Barrière-Fouchet et al. Reference Barrière-Fouchet, Murray, Humble and Maddison2005; Laibe & Price Reference Laibe and Price2012a, Reference Laibe and Price2012b; Lorén-Aguilar & Bate Reference Lorén-Aguilar and Bate2014; Booth, Sijacki, & Clarke Reference Booth, Sijacki and Clarke2015). In grid-based methods, the dust is typically simulated as a particle population, with the hydrodynamics computed on the grid (e.g. Paardekooper Reference Paardekooper2007; Lyra et al. Reference Lyra, Johansen, Klahr and Piskunov2008a; Miniati Reference Miniati2010; Bai & Stone Reference Bai and Stone2010a; Flock et al. Reference Flock, Ruge, Dzyurkevich, Henning, Klahr and Wolf2015; Baruteau & Zhu Reference Baruteau and Zhu2016; Yang & Johansen Reference Yang and Johansen2016)—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 Reference Monaghan1997; Bai & Stone Reference Bai and Stone2010a; Miniati Reference Miniati2010; Laibe & Price Reference Laibe and Price2012b). However, Laibe & Price (Reference Laibe and Price2012a) 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 (Reference Miniati2010) 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 Reference Lorén-Aguilar and Bate2014).

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 Reference Laibe and Price2014a, Reference Laibe and Price2014b, Reference Laibe and Price2014c; Price & Laibe Reference Price and Laibe2015; Hutchison et al. Reference Hutchison, Price, Laibe and Maddison2016). The same approach on a grid means evolving the dust density on the grid (called a ‘two-fluid’ approach by Miniati Reference Miniati2010—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 (Reference Johansen and Klahr2005) 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 (Reference Laibe and Price2014a) and implemented in SPH by Price & Laibe (Reference Price and Laibe2015) 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 (Reference Laibe and Price2014a).

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 St > 1), turbulent motion (St R − 1/2 e ~ 10−4 in astrophysical flows, although Reynolds nymbers, Re ≳ 103, are rarely achieved numerically, Falkovich, Fouxon, & Stepanov Reference Falkovich, Fouxon and Stepanov2002; Ormel & Cuzzi Reference Ormel and Cuzzi2007), 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 Reference Chalons, Kah and Massot2012; Chalons, Fox, & Massot Reference Chalons, Fox and Massot2010; Yuan & Fox Reference Yuan and Fox2011; Yuan, Laurent, & Fox Reference Yuan, Laurent and Fox2012).

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 (> 104 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 (Reference Brauer, Dullemond and Henning2008) and Birnstiel, Dullemond, & Brauer (Reference Birnstiel, Dullemond and Brauer2010) 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. Reference Pinilla, Ovelar, Ataiee, Benisty, van Dishoeck and Min2015, Reference Pinilla, Klarmann, Birnstiel, Benisty, Dominik and Dullemond2016). Similarly, Dipierro et al. (Reference Dipierro, Pinilla, Lodato and Testi2015a) applied this approach to self-gravitating discs in order to predict scattered light images. Miyake, Suzuki, & Inutsuka (Reference Miyake, Suzuki and Inutsuka2016) 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. (Reference Rosotti, Juhasz, Booth and Clarke2016) 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. (Reference Gonzalez, Laibe, Maddison, Pinte and Ménard2015b) 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 Reference Kato, Fujimoto and Ida2012; Taki, Fujimoto, & Ida Reference Taki, Fujimoto and Ida2016).

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 (Reference Paardekooper2007), Lyra et al. (Reference Lyra, Johansen, Klahr and Piskunov2008a), Gonzalez et al. (Reference Gonzalez, Laibe, Maddison, Pinte and Ménard2015a, Reference Gonzalez, Laibe, Maddison, Pinte and Ménard2015b)]. 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 (Reference Laibe and Price2014c) 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 (Reference Paardekooper2007) and Lyra et al. (Reference Lyra, Johansen, Klahr and Piskunov2008a) 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. Reference Gonzalez, Pinte, Maddison, Ménard and Fouchet2012; Dipierro et al. Reference Dipierro, Price, Laibe, Hirsh, Cerioli and Lodato2015b). 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 (Reference Bai and Stone2010a), Laibe & Price (Reference Laibe and Price2014c), 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 (Reference Spitzer1978), Rybicki & Lightman (Reference Rybicki and Lightman1979), Osterbrock & Ferland (Reference Osterbrock, Ferland, Osterbrock and Ferland2006)]. 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 Reference Levermore and Pomraning1981; Whitehouse & Bate Reference Whitehouse and Bate2004; Whitehouse, Bate, & Monaghan Reference Whitehouse, Bate and Monaghan2005). 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 Reference Owen, Ercolano and Clarke2014; Kuiper & Klessen Reference Kuiper and Klessen2013). Most modern applications of FLD account for this failure at low optical depth by using boundary conditions (e.g. Mayer et al. Reference Mayer, Lufkin, Quinn and Wadsley2007), or using hybrid methods to allow the system to radiate energy away from optically thin regions (e.g. Boley et al. Reference Boley, Durisen, Nordlund and Lord2007; Forgan et al. Reference Forgan, Rice, Stamatellos and Whitworth2009). 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 Reference Alexander, Clarke and Pringle2006a, Reference Alexander, Clarke and Pringle2006b; Owen et al. Reference Owen, Ercolano, Clarke and Alexander2010; Owen, Ercolano, & Clarke Reference Owen, Ercolano and Clarke2011; Owen, Clarke, & Ercolano Reference Owen, Clarke and Ercolano2012; Haworth, Clarke, & Owen Reference Haworth, Clarke and Owen2016b).

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 Reference Lucy1999), which is used by the well-known codes radmc-3d (Dullemond Reference Dullemond2012), mcmax (Min et al. Reference Min, Dullemond, Dominik, de Koter and Hovenier2009), hyperion (Robitaille Reference Robitaille2011), mcfost (Pinte et al. Reference Pinte, Ménard, Duchêne and Bastien2006), and torus (Harries Reference Harries2015, 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 Reference Abel and Wandelt2002) and short characteristic (e.g. Davis, Stone, & Jiang Reference Davis, Stone and Jiang2012) 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 Reference Kuiper and Klessen2013; Owen et al. Reference Owen, Ercolano and Clarke2014; Ramsey & Dullemond Reference Ramsey and Dullemond2015).

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 Reference Williams and Best2014). 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. Reference Favre, Cleeves, Bergin, Qi and Blake2013) and the relative abundance of isotoplogues does not necessarily scale as in the ISM (Miotello, Bruderer, & van Dishoeck Reference Miotello, Bruderer and van Dishoeck2014b). Furthermore, dust grain evolution and dynamical processes such as instabilities and planet–disc interactions can also affect the chemistry (e.g. Boley et al. Reference Boley, Durisen, Nordlund and Lord2007; Ilee et al. Reference Ilee, Boley, Caselli, Durisen, Hartquist and Rawlings2011; Evans et al. Reference Evans, Ilee, Boley, Caselli, Durisen, Hartquist and Rawlings2015; Öberg et al. Reference Öberg, Guzmán, Furuya, Qi, Aikawa, Andrews, Loomis and Wilner2015a, Reference Öberg, Furuya, Loomis, Aikawa, Andrews, Qi, van Dishoeck and Wilner2015b; Cleeves, Bergin, & Harries Reference Cleeves, Bergin and Harries2015; Huang, Öberg, & Andrews Reference Huang, Öberg and Andrews2016). Although simple CO parameterisations yield useful insights into the global properties of discs (such as the disc mass, e.g. Miotello et al. Reference Miotello, van Dishoeck, Kama and Bruderer2016; Williams & McPartland Reference Williams and McPartland2016), 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 environmentsFootnote 2 . The abundances of these molecules can be subject to change via a large number of chemical reactions (see Caselli Reference Caselli, Kumar, Tafalla and Caselli2005; Henning & Semenov Reference Henning and Semenov2013, 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 AstrochemistryFootnote 3 (UDfA; Millar, Farquhar, & Willacy Reference Millar, Farquhar and Willacy1997; Woodall et al. Reference Woodall, Agúndez, Markwick-Kemper and Millar2007; McElroy et al. Reference McElroy, Walsh, Markwick, Cordiner, Smith and Millar2013), the Ohio State University networksFootnote 4 , and the Kinetic Database for AstrochemistryFootnote 5 (KIDA; Wakelam et al. Reference Wakelam2012). 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 Reference Arrhenius1889; Kooij Reference Kooij1893).

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 Reference Herbst and van Dishoeck2009). 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 H2O-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. (Reference Kouchi, Kudo, Nakano, Arakawa, Watanabe, Sirono, Higa and Maeno2002), Blum & Wurm (Reference Blum and Wurm2008), Johansen et al. (Reference Johansen, Blum, Tanaka, Ormel, Bizzarro and Rickman2014); Musiolik et al. (Reference Musiolik, Teiser, Jankowski and Wurm2016), 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 Xgr, 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 (Reference Caselli, Kumar, Tafalla and Caselli2005).

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 Reference Woitke, Kamp, Woitke and Ilee2015). 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. Reference Ferland2013) and mocassin (Ercolano et al. Reference Ercolano, Barlow, Storey and Liu2003). 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. Reference Röllig2007) 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 Reference Gorti, Dullemond and Hollenbach2009; Woitke et al., Reference Woitke2016). Recently, multi-dimensional numerical approaches to solving PDR chemistry have appeared that do compute the 3D escape probabilities (Bisbas et al. Reference Bisbas, Bell, Viti, Yates and Barlow2012, Reference Bisbas, Haworth, Barlow, Viti, Harries, Bell and Yates2015b) which they do efficiently using healpix (Górski et al. Reference Górski, Hivon, Banday, Wandelt, Hansen, Reinecke and Bartelmann2005).

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. Reference Boley, Durisen, Nordlund and Lord2007; Ilee et al. Reference Ilee, Boley, Caselli, Durisen, Hartquist and Rawlings2011; Evans et al. Reference Evans, Ilee, Boley, Caselli, Durisen, Hartquist and Rawlings2015; Drozdovskaya et al. Reference Drozdovskaya, Walsh, van Dishoeck, Furuya, Marboeuf, Thiabaud, Harsono and Visser2016).

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 H2CO, 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. Reference Boley, Durisen, Nordlund and Lord2007; Ilee et al. Reference Ilee, Boley, Caselli, Durisen, Hartquist and Rawlings2011; Evans et al. Reference Evans, Ilee, Boley, Caselli, Durisen, Hartquist and Rawlings2015).

Recent work has seen an increase in performing chemical evolution calculations alongside the radiative transfer calculations (e.g. Bruderer, Doty, & Benz Reference Bruderer, Doty and Benz2009; Woitke, Kamp, & Thi Reference Woitke, Kamp and Thi2009). 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. (Reference Glover, Federrath, Low and Klessen2010) and in a 1+1D disc framework, there are models such as those by Gorti et al. (Reference Gorti, Dullemond and Hollenbach2009) 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 Reference Clarke and Alexander2016). 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. Reference Adams, Hollenbach, Laughlin and Gorti2004; Facchini et al. Reference Facchini, Clarke and Bisbas2016; Haworth et al. Reference Haworth, Boubert, Facchini, Bisbas and Clarke2016a). 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. (Reference Davis, Stone and Jiang2012), or the hybrid approach by Kuiper & Klessen (Reference Kuiper and Klessen2013) 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. Reference Flock, Fromang, González and Commerçon2013; Tomida, Okuzumi, & Machida Reference Tomida, Okuzumi and Machida2015).

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.


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:

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. Reference Tomida, Okuzumi and Machida2015).

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 C6C9. 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, C18O, HCO+, 12CO, 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.


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. Reference Haworth, Harries, Acreman and Bisbas2015). 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. Reference Testi2015)], James Webb Space Telescope (JWST) or the European-Extremely Large Telescope (E-ELT, e.g. Hippler et al. Reference Hippler, Usuda, Tamura and Ishii2009)?

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.


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. Reference Woitke, Kamp and Thi2009), dali (Bruderer Reference Bruderer2013), ucl-chem (Viti et al. Reference Viti, Collings, Dever, McCoustra and Williams2004, Reference Viti, Jimenez-Serra, Yates, Codella, Vasta, Caselli, Lefloch and Ceccarelli2011), ucl-pdr (Bell et al. Reference Bell, Viti, Williams, Crawford and Price2005, Reference Bell, Roueff, Viti and Williams2006), and the models of Walsh et al. (Reference Walsh, Nomura, Millar and Aikawa2012). 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. Reference Grassi, Bovino, Schleicher and Gianturco2013). 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. Reference Röllig2007). 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 (Reference Ivlev, Akimkin and Caselli2016) 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. (Reference Owen, Ercolano, Clarke and Alexander2010) 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. Reference Haworth, Clarke and Owen2016b) was computed by the dedicated photoionisation code mocassin (Ercolano et al. Reference Ercolano, Barlow, Storey and Liu2003, Reference Ercolano, Young, Drake and Raymond2008). 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 (Reference Woitke, Krueger and Sedlmayr1996a), Woitke, Goeres, & Sedlmayr (Reference Woitke, Goeres and Sedlmayr1996b), Schirrmacher, Woitke, & Sedlmayr (Reference Schirrmacher, Woitke and Sedlmayr2003) 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. (Reference Glover, Federrath, Low and Klessen2010) and Dzyurkevich et al. (Reference Dzyurkevich, Commerçon, Lesaffre and Semenov2016) patch reduced chemical networks into zeus-mp and ramses, respectively. Flock et al. (Reference Flock, Fromang, González and Commerçon2013) 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 Reference Harries2000) 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 Reference Haworth and Harries2012; Harries Reference Harries2015; Haworth et al. Reference Haworth, Harries, Acreman and Bisbas2015). Furthermore, torus-3dpdr is an extension of torus that also includes PDR chemistry with 3D extinction and escape probabilities (Bisbas et al. Reference Bisbas, Haworth, Barlow, Viti, Harries, Bell and Yates2015b). 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. Reference Górski, Hivon, Banday, Wandelt, Hansen, Reinecke and Bartelmann2005). 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. Reference Adams, Hollenbach, Laughlin and Gorti2004; Facchini et al. Reference Facchini, Clarke and Bisbas2016; Haworth et al. Reference Haworth, Boubert, Facchini, Bisbas and Clarke2016a). 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 Reference Harries2015).

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. Reference Motoyama, Morata, Shang, Krasnopolsky and Hasegawa2015).

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. Reference Grassi, Bovino, Schleicher, Prieto, Seifried, Simoncini and Gianturco2014). 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. Reference Burns, Brown, Lecoanet, Oishi and Vasil2016) and snoopy (Lesur Reference Lesur2015) 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 CH4 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 Reference Clark, Glover and Klessen2012, 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 Reference Schive, Tsai and Chiueh2010; Bryan et al. Reference Bryan2014b). However, a GPU implementation of the fargo disc-modelling code was developed by Benítez-Llambay & Masset (Reference Benítez-Llambay and Masset2016), 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. Reference Nickolls, Buck, Garland and Skadron2008) or opencl (Stone, Gohara, & Shi Reference Stone, Gohara and Shi2010) 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 standardFootnote 6 .

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 Reference Jeffers and Reinders2013)—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.


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 dotAstronomyFootnote 7 or AstropyFootnote 8 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. Reference de Val-Borro2006; Röllig et al. Reference Röllig2007; Pinte et al. Reference Pinte, Harries, Min, Watson, Dullemond, Woitke, Ménard and Durán-Rojas2009; Iliev et al. Reference Iliev2009).

A good example of a successful comparison project is the recent StarBench code comparison workshopsFootnote 9 (Bisbas et al. Reference Bisbas2015a). 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. Reference Bisbas2015a). 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 Reference Bate and Burkert1997; Commerçon et al. Reference Commerçon, Hennebelle, Audit, Chabrier and Teyssier2008; Price & Federrath Reference Price and Federrath2010; Hubber, Falle, & Goodwin Reference Hubber, Falle and Goodwin2013). 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 GitHubFootnote 10 ). 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. Reference Grassi, Bovino, Schleicher, Prieto, Seifried, Simoncini and Gianturco2014) and lime (Brinch & Hogerheijde Reference Brinch and Hogerheijde2010) 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.


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.


We thank the anonymous referee for their comments, and positive review of the paper. We would like to acknowledge the Protoplanetary Discussions conference, including the members of the local and scientific organising committees. We also thank the attendees of the conference (Table 3) and those who took part in the discussion sessions chaired by John Ilee and Daniel Price, all of whom helped inform many of the statements made here. We also particularly thank Phil Armitage for comments on the paper.

Through most of this work TJH was funded by the STFC consolidated grant ST/K000985/1 and is now funded by an Imperial Junior research fellowship. JDI gratefully acknowledges support from the DISCSIM project, grant agreement 341137, funded by the European Research Council under ERC-2013-ADG. DHF acknowledges support from the ECOGAL project, grant agreement 291227, funded by the European Research Council under ERC-2011-ADG. DJP gratefully acknowledges funding via Future Fellowship FT130100034 from the Australian Research Council.OP is supported by the Royal Society Dorothy Hodgkin Fellowship.


* This paper was coordinated and written by the first five authors: Haworth, Ilee, Forgan, Facchini, and Price. The additional ‘community authors’, presented alphabetically, made valuable contributions that helped to inform the paper.
Royal Society Dorothy Hodgkin Fellow



Abel, T., & Wandelt, B. D. 2002, MNRAS, 330, L53 j.1365-8711.2002.05206.x 2002MNRAS.330L..53A CrossRefGoogle Scholar
Adams, F. C., Hollenbach, D., Laughlin, G., & Gorti, U. 2004, ApJ, 611, 360 10.1086/421989 2004ApJ. . .611..360A CrossRefGoogle Scholar
Alexander, R. D., & Armitage, P. J. 2009, ApJ, 704, 989 10.1088/0004-637X/704/2/989 2009ApJ. . .704..989A CrossRefGoogle Scholar
Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006a, MNRAS, 369, 216 10.1111/j.1365-2966.2006.10293.x 2006MNRAS.369..216A Google Scholar
Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006b, MNRAS, 369, 229 10.1111/j.1365-2966.2006.10294.x 2006MNRAS.369..229A Google Scholar
ALMA Partnership, et al. 2015, ApJ, 808, L3 10.1088/2041-8205/808/1/L3 2015ApJ. . .808L. . .3A Google Scholar
Andrews, S. M., Wilner, D. J., Espaillat, C., Hughes, A. M., Dullemond, C. P., McClure, M. K., Qi, C., & Brown, J. M. 2011, ApJ, 732, 42 10.1088/0004-637X/732/1/42 2011ApJ. . .732. . .42A Google Scholar
Andrews, S. M., et al. 2016, ApJ, 820, L40 10.3847/2041-8205/820/2/L40 2016ApJ. . .820L..40A Google Scholar
Armitage, P. J. 2011, ARA&A, 49, 195 10.1146/annurev-astro-081710-102521 2011ARA&A..49..195A Google Scholar
Armitage, P. J. 2015, preprint (arXiv:1509.06382)2015arXiv150906382A Google Scholar
Armitage, P. J., Livio, M., & Pringle, J. E. 2001, MNRAS, 324, 705 j.1365-8711.2001.04356.x 2001MNRAS.324..705A Google Scholar
Arrhenius, S. 1889, ZPC, 4, 96 Google Scholar
Bae, J., Hartmann, L., Zhu, Z., & Nelson, R. P. 2014, ApJ, 795, 61 10.1088/0004-637X/795/1/61 2014ApJ. . .795. . .61B CrossRefGoogle Scholar
Bai, X.-N. 2014, ApJ, 791, 137 10.1088/0004-637X/791/2/137 2014ApJ. . .791..137B CrossRefGoogle Scholar
Bai, X.-N., & Stone, J. M. 2010a, ApJS, 190, 297 10.1088/0067-0049/190/2/297 2010ApJS..190..297B Google Scholar
Bai, X.-N., & Stone, J. M. 2010b, ApJ, 722, 1437 10.1088/0004-637X/722/2/1437 2010ApJ. . .722.1437B CrossRefGoogle Scholar
Bai, X.-N., & Stone, J. M. 2014, ApJ, 796, 31 10.1088/0004-637X/796/1/31 2014ApJ. . .796. . .31B CrossRefGoogle Scholar
Bai, X.-N., Ye, J., Goodman, J., & Yuan, F. 2016, ApJ, 818, 152 10.3847/0004-637X/818/2/152 2016ApJ. . .818..152B CrossRefGoogle Scholar
Balbus, S. A. 2003, ARA&A, 41, 555 10.1146/annurev.astro.41.081401.155207 2003ARA&A..41..555B Google Scholar
Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 10.1086/170270 1991ApJ. . .376..214B CrossRefGoogle Scholar
Bally, J., O’Dell, C. R., & McCaughrean, M. J. 2000, AJ, 119, 2919 10.1086/301385 2000AJ. . ..119.2919B Google Scholar
Barrière-Fouchet, L., Gonzalez Murray, J. R., Humble, R. J., & Maddison, S. T. 2005, A&A, 443, 185 10.1051/0004-6361:20042249 2005A&A. . .443..185B Google Scholar
Baruteau, C., & Masset, F. 2008, ApJ, 678, 483 10.1086/529487 2008ApJ. . .678..483B Google Scholar
Baruteau, C., & Zhu, Z. 2016, MNRAS, 458, 3927 10.1093/mnras/stv2527 2016MNRAS.458.3927B Google Scholar
Bate, M. R., & Burkert, A. 1997, MNRAS, 288, 1060 10.1093/mnras/288.4.1060 1997MNRAS.288.1060B Google Scholar
Bell, T. A., Roueff, E., Viti, S., & Williams, D. A. 2006, MNRAS, 371, 1865 10.1111/j.1365-2966.2006.10817.x 2006MNRAS.371.1865B Google Scholar
Bell, T. A., Viti, S., Williams, D. A., Crawford, I. A., & Price, R. J. 2005, MNRAS, 357, 961 10.1111/j.1365-2966.2005.08693.x 2005MNRAS.357..961B Google Scholar
Benisty, M., et al. 2015, A&A, 578, L6 10.1051/0004-6361/201526011 2015A&A. . .578L. . .6B Google Scholar
Benítez-Llambay, P., & Masset, F. S. 2016, ApJS, 223, 11 10.3847/0067-0049/223/1/11 2016ApJS..223. . .11B Google Scholar
Berger, M. J., & Colella, P. 1989, JCoPh, 82, 64 10.1016/0021-9991(89)90035-1 1989JCoPh..82. . .64B Google Scholar
Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79 10.1051/0004-6361/200913731 2010A&A. . .513A..79B Google Scholar
Bisbas, T. G., Bell, T. A., Viti, S., Yates, J., & Barlow, M. J. 2012, MNRAS, 427, 2100 10.1111/j.1365-2966.2012.22077.x 2012MNRAS.427.2100B Google Scholar
Bisbas, T. G., Haworth, T. J., Barlow, M. J., Viti, S., Harries, T. J., Bell, T., & Yates, J. A. 2015b, MNRAS, 454, 2828 10.1093/mnras/stv2156 2015MNRAS.454.2828B Google Scholar
Bisbas, T. G., et al. 2015a, MNRAS, 453, 1324 10.1093/mnras/stv1659 2015MNRAS.453.1324B Google Scholar
Bitsch, B., Lambrechts, M., & Johansen, A. 2015, A&A, 582, A112 10.1051/0004-6361/201526463 2015A&A. . .582A.112B Google Scholar
Blum, J., & Wurm, G. 2008, ARA&A, 46, 21 10.1146/annurev.astro.46.060407.145152 2008ARA&A..46. . .21B Google Scholar
Bodenheimer, P. 1995, ARA&A, 33, 199 10.1146/annurev.aa.33.090195.001215 1995ARA&A..33..199B Google Scholar
Boley, A. C., Durisen, R. H., Nordlund, Å., & Lord, J. 2007, ApJ, 665, 1254 10.1086/519767 2007ApJ. . .665.1254B Google Scholar
Boneberg, D. M., Panić, O., Haworth, T. J., Clarke, C. J., & Min, M. 2016, MNRAS, 461, 385 10.1093/mnras/stw1325 2016MNRAS.461..385B Google Scholar
Booth, R. A., Sijacki, D., & Clarke, C. J. 2015, MNRAS, 452, 3932 10.1093/mnras/stv1486 2015MNRAS.452.3932B CrossRefGoogle Scholar
Brandenburg, A., & Dobler, W. 2002, CoPhC, 147, 471 10.1016/S0010-4655(02)00334-X 2002CoPhC.147..471B Google Scholar
Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859 10.1051/0004-6361:20077759 2008A&A. . .480..859B Google Scholar
Brinch, C., & Hogerheijde, M. R. 2010, A&A, 523, A25 10.1051/0004-6361/201015333 2010A&A. . .523A..25B Google Scholar
Bruderer, S. 2013, A&A, 559, A46 10.1051/0004-6361/201321171 2013A&A. . .559A..46B Google Scholar
Bruderer, S., Doty, S. D., & Benz, A. O. 2009, ApJS, 183, 179 10.1088/0067-0049/183/2/179 2009ApJS..183..179B Google Scholar
Bryan, G. L., et al. 2014a, ApJS, 211, 19 10.1088/0067-0049/211/2/19 2014ApJS..211. . .19B Google Scholar
Bryan, G. L., et al. 2014b, ApJS, 211, 19 10.1088/0067-0049/211/2/19 2014ApJS..211. . .19B Google Scholar
Burns, K., Brown, B., Lecoanet, D., Oishi, J., & Vasil, G. 2016, Dedalus: Flexible framework for spectrally solving differential equations, Astrophysics Source Code Library (ascl:1603.015)Google Scholar
Canovas, H., Caceres, C., Schreiber, M. R., Hardy, A., Cieza, L., Ménard, F., & Hales, A. 2016, MNRAS, 458, L29 10.1093/mnrasl/slw006 2016MNRAS.458L..29C CrossRefGoogle Scholar
Casassus, S. 2016, PASA, 33, e013 10.1017/pasa.2016.7 2016PASA. . .33. . .13C Google Scholar
Casassus, S., et al. 2015, ApJ, 811, 92 10.1088/0004-637X/811/2/92 2015ApJ. . .811. . .92C Google Scholar
Caselli, P. 2005, in Cores to Clusters: Star Formation with Next Generation Telescopes, eds. Kumar, M. S. N., Tafalla, M., & Caselli, P. (New York: Springer), 47 astro-ph/0504298 CrossRefGoogle Scholar
Casse, F., Meliani, Z., & Sauty, C. 2007, Ap&SS, 311, 57 10.1007/s10509-007-9561-1 2007Ap&SS.311. . .57C Google Scholar
Chalons, C., Fox, R., & Massot, M. 2010, in Studying Turbulence Using Numerical Simulation Databases (Stanford: Stanford University, Center for Turbulence Research), 347 Google Scholar
Chalons, C., Kah, D., & Massot, M. 2012, Comm. in Math. Sciences, 10, 1241 Google Scholar
Chambers, J. E. 2016, ApJ, 825, 63 10.3847/0004-637X/825/1/63 2016ApJ. . .825. . .63C Google Scholar
Clark, P. C., Glover, S. C. O., & Klessen, R. S. 2012, MNRAS, 420, 745 10.1111/j.1365-2966.2011.20087.x 2012MNRAS.420..745C Google Scholar
Clarke, C. J., & Alexander, R. D. 2016, MNRAS 10.1093/mnras/stw1178 2016MNRAS.tmp..880C Google Scholar
Clarke, C. J., & Pringle, J. E. 1993, MNRAS, 261, 190 10.1093/mnras/261.1.190 1993MNRAS.261..190C Google Scholar
Cleeves, L. I., Bergin, E. A., & Harries, T. J. 2015, ApJ, 807, 2 10.1088/0004-637X/807/1/2 2015ApJ. . .807. . ..2C Google Scholar
Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2008, A&A, 482, 371 10.1051/0004-6361:20078591 2008A&A. . .482..371C Google Scholar
Dai, F., Facchini, S., Clarke, C. J., & Haworth, T. J. 2015, MNRAS, 449, 1996 10.1093/mnras/stv403 2015MNRAS.449.1996D Google Scholar
Davis, S. W., Stone, J. M., & Jiang, Y.-F. 2012, ApJS, 199, 9 10.1088/0067-0049/199/1/9 2012ApJS..199. . ..9D Google Scholar
de Juan Ovelar, M., Kruijssen, J. M. D., Bressert, E., Testi, L., Bastian, N., & Cánovas, H. 2012, A&A, 546, L1 10.1051/0004-6361/201219627 2012A&A. . .546L. . .1D Google Scholar
de Val-Borro, M., et al. 2006, MNRAS, 370, 529 10.1111/j.1365-2966.2006.10488.x 2006MNRAS.370..529D Google Scholar
Dipierro, G., Pinilla, P., Lodato, G., & Testi, L. 2015a, MNRAS, 451, 974 10.1093/mnras/stv970 2015MNRAS.451..974D CrossRefGoogle Scholar
Dipierro, G., Price, D., Laibe, G., Hirsh, K., Cerioli, A., & Lodato, G. 2015b, MNRAS, 453, L73 10.1093/mnrasl/slv105 2015MNRAS.453L..73D Google Scholar
Doğan, S., Nixon, C., King, A., & Price, D. J. 2015, MNRAS, 449, 1251 10.1093/mnras/stv347 2015MNRAS.449.1251D Google Scholar
Drozdovskaya, M. N., Walsh, C., van Dishoeck, E. F., Furuya, K., Marboeuf, U., Thiabaud, A., Harsono, D., & Visser, R. 2016, MNRAS, 462, 977 2016MNRAS.462..977D Google Scholar
Duffell, P. C., & MacFadyen, A. I. 2011, ApJS, 197, 15 10.1088/0067-0049/197/2/15 2011ApJS..197. . .15D Google Scholar
Dullemond, C. P. 2012, RADMC-3D: A multi-purpose radiative transfer tool, Astrophysics Source Code Library (ascl:1202.015)Google Scholar
Dullemond, C. P., Hollenbach, D., Kamp, I., & D’Alessio, P. 2007, Protostars and Planets V, 5552007prpl.conf..555D Google Scholar
Durisen, R. H., Boss, A. P., Mayer, L., Nelson, A. F., Quinn, T., Rice, W. K. M. 2007, Protostars and Planets V, 6072007prpl.conf..607D Google Scholar
Dzyurkevich, N., Commerçon, B., Lesaffre, P., & Semenov, D. 2016, preprint (arXiv:1605.08032)2016arXiv160508032D Google Scholar
Ercolano, B., Barlow, M. J., Storey, P. J., & Liu, X.-W. 2003, MNRAS, 340, 1136 j.1365-8711.2003.06371.x 2003MNRAS.340.1136E Google Scholar
Ercolano, B., Young, P. R., Drake, J. J., & Raymond, J. C. 2008, ApJS, 175, 534 10.1086/524378 2008ApJS..175..534E Google Scholar
Evans, M. G., Ilee, J. D., Boley, A. C., Caselli, P., Durisen, R. H., Hartquist, T. W., & Rawlings, J. M. C. 2015, MNRAS, 453, 1147 10.1093/mnras/stv1698 2015MNRAS.453.1147E Google Scholar
Facchini, S., Clarke, C. J., & Bisbas, T. G. 2016, MNRAS, 457, 3593 10.1093/mnras/stw240 2016MNRAS.457.3593F Google Scholar
Facchini, S., Lodato, G., & Price, D. J. 2013, MNRAS, 433, 2142 10.1093/mnras/stt877 2013MNRAS.433.2142F Google Scholar
Falkovich, G., Fouxon, A., & Stepanov, M. G. 2002, Nature, 419, 151 10.1038/nature00983 2002Natur.419..151F CrossRefGoogle Scholar
Favre, C., Cleeves, L. I., Bergin, E. A., Qi, C., & Blake, G. A. 2013, ApJ, 776, L38 10.1088/2041-8205/776/2/L38 2013ApJ. . .776L..38F CrossRefGoogle Scholar
Ferland, G. J., et al. 2013, RMXAA, 49, 137 2013RMxAA..49..137F Google Scholar
Flaherty, K. M., Hughes, A. M., Rosenfeld, K. A., Andrews, S. M., Chiang, E., Simon, J. B., Kerzner, S., & Wilner, D. J. 2015, ApJ, 813, 99 10.1088/0004-637X/813/2/99 2015ApJ. . .813. . .99F Google Scholar
Flock, M., Fromang, S., González, M., & Commerçon, B. 2013, A&A, 560, A43 10.1051/0004-6361/201322451 2013A&A. . .560A..43F Google Scholar
Flock, M., Ruge, J. P., Dzyurkevich, N., Henning, T., Klahr, H., & Wolf, S. 2015, A&A, 574, A68 10.1051/0004-6361/201424693 2015A&A. . .574A..68F Google Scholar
Forgan, D., Parker, R. J., & Rice, K. 2015, MNRAS, 447, 836 10.1093/mnras/stu2504 2015MNRAS.447..836F Google Scholar
Forgan, D., & Rice, K. 2013, MNRAS, 430, 2082 10.1093/mnras/stt032 2013MNRAS.430.2082F Google Scholar
Forgan, D., Rice, K., Stamatellos, D., & Whitworth, A. 2009, MNRAS, 394, 882 10.1111/j.1365-2966.2008.14373.x 2009MNRAS.394..882F Google Scholar
Fragner, M. M., & Nelson, R. P. 2010, A&A, 511, A77 10.1051/0004-6361/200913088 2010A&A. . .511A..77F Google Scholar
Gammie, C. F. 1996, ApJ, 457, 355 10.1086/176735 1996ApJ. . .457..355G Google Scholar
Garufi, A., et al. 2013, A&A, 560, A105 10.1051/0004-6361/201322429 2013A&A. . .560A.105G Google Scholar
Glover, S. C. O., Federrath, C., Mac Low, M.-M., & Klessen, R. S. 2010, MNRAS, 404, 2 10.1111/j.1365-2966.2009.15718.x 2010MNRAS.404. . ..2G Google Scholar
Gonzalez, J.-F., Pinte, C., Maddison, S. T., Ménard, F., & Fouchet, L. 2012, A&A, 547, A58 10.1051/0004-6361/201218806 2012A&A. . .547A..58G Google Scholar
Gonzalez, J.-F., Laibe, G., Maddison, S. T., Pinte, C., & Ménard, F. 2015a, P&SS, 116, 48 10.1016/j.pss.2015.05.018 2015P&SS..116. . .48G Google Scholar
Gonzalez, J.-F., Laibe, G., Maddison, S. T., Pinte, C., & Ménard, F. 2015b, MNRAS, 454, L36 10.1093/mnrasl/slv120 2015MNRAS.454L..36G Google Scholar
Górski, K. M., Hivon, E., Banday, A. J., Wandelt, B. D., Hansen, F. K., Reinecke, M., & Bartelmann, M. 2005, ApJ, 622, 759 10.1086/427976 2005ApJ. . .622..759G Google Scholar
Gorti, U., Dullemond, C. P., & Hollenbach, D. 2009, ApJ, 705, 1237 10.1088/0004-637X/705/2/1237 2009ApJ. . .705.1237G Google Scholar
Gorti, U., & Hollenbach, D. 2009, ApJ, 690, 1539 10.1088/0004-637X/690/2/1539 2009ApJ. . .690.1539G CrossRefGoogle Scholar
Gorti, U., Hollenbach, D., & Dullemond, C. P. 2015, ApJ, 804, 29 10.1088/0004-637X/804/1/29 2015ApJ. . .804. . .29G Google Scholar
Grassi, T., Bovino, S., Schleicher, D., & Gianturco, F. A. 2013, MNRAS, 431, 1659 10.1093/mnras/stt284 2013MNRAS.431.1659G Google Scholar
Grassi, T., Bovino, S., Schleicher, D. R. G., Prieto, J., Seifried, D., Simoncini, E., & Gianturco, F. A. 2014, MNRAS, 439, 2386 10.1093/mnras/stu114 2014MNRAS.439.2386G CrossRefGoogle Scholar
Guan, X., & Gammie, C. F. 2008, ApJS, 174, 145 10.1086/521147 2008ApJS..174..145G Google Scholar
Harries, T. J. 2000, MNRAS, 315, 722 j.1365-8711.2000.03505.x 2000MNRAS.315..722H Google Scholar
Harries, T. J. 2015, MNRAS, 448, 3156 10.1093/mnras/stv158 2015MNRAS.448.3156H Google Scholar
Hartmann, L. 1998, Accretion Processes in Star Formation (New York: Cambridge University Press)Google Scholar
Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742 10.1086/175311 1995ApJ. . .440..742H CrossRefGoogle Scholar
Haworth, T. J., Boubert, D., Facchini, S., Bisbas, T. G., & Clarke, C. J. 2016a, MNRAS, 463, 3616 10.1093/mnras/stw2280 463/4/3616 Google Scholar
Haworth, T. J., Clarke, C. J., & Owen, J. E. 2016b, MNRAS, 457, 1905 10.1093/mnras/stv3016 2016MNRAS.457.1905H Google Scholar
Haworth, T. J., & Harries, T. J. 2012, MNRAS, 420, 562 10.1111/j.1365-2966.2011.20062.x 2012MNRAS.420..562H Google Scholar
Haworth, T. J., Harries, T. J., Acreman, D. M., & Bisbas, T. G. 2015, MNRAS, 453, 2277 10.1093/mnras/stv1814 2015MNRAS.453.2277H Google Scholar
Henney, W. J., O’Dell, C. R., Meaburn, J., Garrington, S. T., & Lopez, J. A. 2002, ApJ, 566, 315 10.1086/338055 2002ApJ. . .566..315H Google Scholar
Henning, T., & Semenov, D. 2013, ChRv, 113, 9016 10.1021/cr400128p 2013ChRv..113.9016H Google Scholar
Herbst, E., & van Dishoeck, E. F. 2009, ARA&A, 47, 427 10.1146/annurev-astro-082708-101654 2009ARA&A..47..427H Google Scholar
Hippler, S., et al. 2009, in AIP Conf. Ser., Vol. 1158, ed. Usuda, T., Tamura, M., & Ishii, M. (San Francisco: AIP), 333 2009AIPC.1158..333H Google Scholar
Hopkins, P. F. 2015, MNRAS, 450, 53 10.1093/mnras/stv195 2015MNRAS.450. . .53H Google Scholar
Huang, J., Öberg, K. I., & Andrews, S. M. 2016, ApJ, 823, L18 10.3847/2041-8205/823/1/L18 2016ApJ. . .823L..18H Google Scholar
Hubber, D. A., Falle, S. A. E. G., & Goodwin, S. P. 2013, MNRAS, 432, 711 10.1093/mnras/stt509 2013MNRAS.432..711H Google Scholar
Hutchison, M. A., Price, D. J., Laibe, G., & Maddison, S. T. 2016, MNRAS, 461, 742 10.1093/mnras/stw1126 2016MNRAS.461..742H Google Scholar
Ida, S., Guillot, T., & Morbidelli, A. 2016, A&A, 591, A72 10.1051/0004-6361/201628099 2016A&A. . .591A..72I Google Scholar
Ilee, J. D., Boley, A. C., Caselli, P., Durisen, R. H., Hartquist, T. W., & Rawlings, J. M. C. 2011, MNRAS, 417, 2950 10.1111/j.1365-2966.2011.19455.x 2011MNRAS.417.2950I Google Scholar
Iliev, I. T., et al. 2009, MNRAS, 400, 1283 10.1111/j.1365-2966.2009.15558.x 2009MNRAS.400.1283I Google Scholar
Ivlev, A. V., Akimkin, V. V., & Caselli, P. 2016, preprint (arXiv:1607.03701)2016arXiv160703701I Google Scholar
Jeffers, J., & Reinders, J. 2013, Intel Xeon Phi coprocessor high-performance programming. NewnesGoogle Scholar
Jin, S., Li, S., Isella, A., Li, H., & Ji, J. 2016, ApJ, 818, 76 10.3847/0004-637X/818/1/76 2016ApJ. . .818. . .76J Google Scholar
Johansen, A., & Klahr, H. 2005, ApJ, 634, 1353 10.1086/497118 2005ApJ. . .634.1353J Google Scholar
Johansen, A., Oishi, J. S., Mac Low, M.-M., Klahr, H., Henning, T., & Youdin, A. 2007, Nature, 448, 1022 10.1038/nature06086 2007Natur.448.1022J Google Scholar
Johansen, A., Blum, J., Tanaka, H., Ormel, C., Bizzarro, M., & Rickman, H. 2014, Protostars and Planets VI, 54710.2458/azu_uapress_9780816531240-ch024 2014prpl.conf..547J Google Scholar
Kama, M., et al. 2016, A&A, 592, A83 2016A%26A. . .592A..83K Google Scholar
Kato, M. T., Fujimoto, M., & Ida, S. 2012, ApJ, 747, 11 10.1088/0004-637X/747/1/11 2012ApJ. . .747. . .11K Google Scholar
Kooij, D. M. 1893, ZPC, 12, 155 Google Scholar
Kouchi, A., Kudo, T., Nakano, H., Arakawa, M., Watanabe, N., Sirono, S.-i., Higa, M., & Maeno, N. 2002, ApJ, 566, L121 10.1086/339618 2002ApJ. . .566L.121K Google Scholar
Krijt, S., & Ciesla, F. J. 2016, ApJ, 822, 111 10.3847/0004-637X/822/2/111 2016ApJ. . .822..111K Google Scholar
Kuiper, R., & Klessen, R. S. 2013, A&A, 555, A7 10.1051/0004-6361/201321404 2013A&A. . .555A. . .7K Google Scholar
Kunz, M. W., & Lesur, G. 2013, MNRAS, 434, 2295 10.1093/mnras/stt1171 2013MNRAS.434.2295K Google Scholar
Laibe, G., & Price, D. J. 2012a, MNRAS, 420, 2345 10.1111/j.1365-2966.2011.20202.x 2012MNRAS.420.2345L Google Scholar
Laibe, G., & Price, D. J. 2012b, MNRAS, 420, 2365 10.1111/j.1365-2966.2011.20201.x 2012MNRAS.420.2365L Google Scholar
Laibe, G., & Price, D. J. 2014a, MNRAS, 440, 2136 10.1093/mnras/stu355 2014MNRAS.440.2136L Google Scholar
Laibe, G., & Price, D. J. 2014b, MNRAS, 440, 2147 10.1093/mnras/stu359 2014MNRAS.440.2147L CrossRefGoogle Scholar
Laibe, G., & Price, D. J. 2014c, MNRAS, 444, 1940 10.1093/mnras/stu1367 2014MNRAS.444.1940L Google Scholar
Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32 10.1051/0004-6361/201219127 2012A&A. . .544A..32L Google Scholar
Lesur, G. 2015, Snoopy: General purpose spectral solver, Astrophysics Source Code Library (ascl:1505.022)Google Scholar
Lesur, G., & Papaloizou, J. C. B. 2010, A&A, 513, A60 10.1051/0004-6361/200913594 2010A&A. . .513A..60L Google Scholar
Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321 10.1086/159157 1981ApJ. . .248..321L Google Scholar
Lewis, B. T., Bate, M. R., & Tricco, T. S. 2016, Proc. 11th Int. SPHERIC Workshop, Munich, 13–16 June2016arXiv160606972L Google Scholar
Lodato, G. 2008, NewAR, 52, 21 10.1016/j.newar.2008.04.002 2008NewAR..52. . .21L Google Scholar
Lodato, G., & Facchini, S. 2013, MNRAS, 433, 2157 10.1093/mnras/stt878 2013MNRAS.433.2157L Google Scholar
Lodato, G., & Price, D. J. 2010, MNRAS, 405, 1212 10.1111/j.1365-2966.2010.16526.x 2010MNRAS.405.1212L Google Scholar
Lomax, O., Whitworth, A. P., & Hubber, D. A. 2015, MNRAS, 449, 662 10.1093/mnras/stv310 2015MNRAS.449..662L CrossRefGoogle Scholar
Lorén-Aguilar, P., & Bate, M. R. 2014, MNRAS, 443, 927 10.1093/mnras/stu1173 2014MNRAS.443..927L Google Scholar
Lorén-Aguilar, P., & Bate, M. R. 2015, MNRAS, 453, L78 10.1093/mnrasl/slv109 2015MNRAS.453L..78L Google Scholar
Lorén-Aguilar, P., & Bate, M. R. 2016, MNRAS, 457, L54 10.1093/mnrasl/slv206 2016MNRAS.457L..54L Google Scholar
Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805 10.1086/306900 1999ApJ. . .513..805L Google Scholar
Lovelace, R. V. E., Romanova, M. M., Lii, P., & Dyda, S. 2014, ComAC, 1, 3 10.1186/s40668-014-0003-5 2014ComAC. . .1. . ..3L Google Scholar
Lubow, S. H., Papaloizou, J. C. B., & Pringle, J. E. 1994, MNRAS, 267, 235 10.1093/mnras/267.2.235 1994MNRAS.267..235L Google Scholar
Lucy, L. B. 1999, A&A, 344, 282 1999A&A. . .344..282L Google Scholar
Lyra, W., Johansen, A., Klahr, H., & Piskunov, N. 2008a, A&A, 479, 883 10.1051/0004-6361:20077948 2008A&A. . .479..883L Google Scholar
Lyra, W., Johansen, A., Klahr, H., & Piskunov, N. 2008b, A&A, 491, L41 10.1051/0004-6361:200810626 2008A&A. . .491L..41L Google Scholar
Lyra, W., Johansen, A., Klahr, H., & Piskunov, N. 2009, A&A, 493, 1125 10.1051/0004-6361:200810797 2009A&A. . .493.1125L Google Scholar
Lyra, W., & Klahr, H. 2011, A&A, 527, A138 10.1051/0004-6361/201015568 2011A&A. . .527A.138L Google Scholar
Marino, S., Perez, S., & Casassus, S. 2015, ApJ, 798, L44 10.1088/2041-8205/798/2/L44 2015ApJ. . .798L..44M Google Scholar
Maron, J. L., McNally, C. P., & Mac Low, M.-M. 2012, ApJS, 200, 6 10.1088/0067-0049/200/1/6 2012ApJS..200. . ..6M Google Scholar
Masset, F. 2000, A&AS, 141, 165 10.1051/aas:2000116 2000A&AS..141..165M Google Scholar
Mayer, L., Lufkin, G., Quinn, T., & Wadsley, J. 2007, ApJ, 661, L77 10.1086/518433 2007ApJ. . .661L..77M Google Scholar
McElroy, D., Walsh, C., Markwick, A. J., Cordiner, M. A., Smith, K., & Millar, T. J. 2013, A&A, 550, A36 10.1051/0004-6361/201220465 2013A&A. . .550A..36M Google Scholar
McNally, C. P., Maron, J. L., & Mac Low, M.-M. 2012, ApJS, 200, 7 10.1088/0067-0049/200/1/7 2012ApJS..200. . ..7M Google Scholar
Meijer, J., Dominik, C., de Koter, A., Dullemond, C. P., van Boekel R., & Waters, L. B. F. M. 2008, A&A, 492, 451 10.1051/0004-6361:20077967 2008A&A. . .492..451M Google Scholar
Meru, F. 2015, MNRAS, 454, 2529 10.1093/mnras/stv2128 2015MNRAS.454.2529M Google Scholar
Meru, F., & Bate, M. R. 2010, MNRAS, 406, 2279 10.1111/j.1365-2966.2010.16867.x 2010MNRAS.406.2279M Google Scholar
Meru, F., & Bate, M. R. 2011, MNRAS, 411, L1 10.1111/j.1745-3933.2010.00978.x Google Scholar
Michael, S., Steiman-Cameron, T. Y., Durisen, R. H., & Boley, A. C. 2012, ApJ, 746, 98 10.1088/0004-637X/746/1/98 Google Scholar
Mignone, A., Bodo, G., Massaglia, S., Matsakos, T., Tesileanu, O., Zanni, C., & Ferrari, A. 2007, ApJS, 170, 228 10.1086/513316 2007ApJS..170..228M Google Scholar
Millar, T. J., Farquhar, P. R. A., & Willacy, K. 1997, A&AS, 12110.1051/aas:1997118 1997A&AS..121..139M Google Scholar
Min, M., Dullemond, C. P., Dominik, C., de Koter, A., & Hovenier, J. W. 2009, A&A, 497, 155 10.1051/0004-6361/200811470 2009A&A. . .497..155M Google Scholar
Miniati, F. 2010, JCoPh, 229, 3916 10.1016/ 2010JCoPh.229.3916M Google Scholar
Miotello, A., Bruderer, S., & van Dishoeck, E. F. 2014b, A&A, 572, A96 10.1051/0004-6361/201424712 2014A&A. . .572A..96M Google Scholar
Miotello, A., Testi, L., Lodato, G., Ricci, L., Rosotti, G., Brooks, K., Maury, A., & Natta, A. 2014a, A&A, 567, A32 10.1051/0004-6361/201322945 2014A&A. . .567A..32M Google Scholar
Miotello, A., van Dishoeck, E. F., Kama, M., & Bruderer, S. 2016, preprint (arXiv:1605.07780)2016arXiv160507780M Google Scholar
Miyake, T., Suzuki, T. K., & Inutsuka, S.-i. 2016, ApJ, 821, 3 10.3847/0004-637X/821/1/3 2016ApJ. . .821. . ..3M Google Scholar
Mocz, P., Vogelsberger, M., Pakmor, R., Genel, S., Springel, V., & Hernquist, L. 2015, MNRAS, 452, 3853 10.1093/mnras/stv1598 2015MNRAS.452.3853M Google Scholar
Mocz, P., Vogelsberger, M., Sijacki, D., Pakmor, R., & Hernquist, L. 2014, MNRAS, 437, 397 10.1093/mnras/stt1890 2014MNRAS.437..397M Google Scholar
Monaghan, J. J. 1992, ARA&A, 30, 543 10.1146/annurev.aa.30.090192.002551 1992ARA&A..30..543M Google Scholar
Monaghan, J. J. 1997, JCoPh, 138, 801 10.1006/jcph.1997.5846 1997JCoPh.138..801M Google Scholar
Monaghan, J. J., & Kocharyan, A. 1995, CoPhC, 87, 225 10.1016/0010-4655(94)00174-Z 1995CoPhC..87..225M Google Scholar
Morbidelli, A., & Nesvorny, D. 2012, A&A, 546, A18 10.1051/0004-6361/201219824 2012A&A. . .546A..18M Google Scholar
Motoyama, K., Morata, O., Shang, H., Krasnopolsky, R., & Hasegawa, T. 2015, ApJ, 808, 46 10.1088/0004-637X/808/1/46 2015ApJ. . .808. . .46M Google Scholar
Muñoz, D. J., Kratter, K., Springel, V., & Hernquist, L. 2014, MNRAS, 445, 3475 10.1093/mnras/stu1918 2014MNRAS.445.3475M Google Scholar
Musiolik, G., Teiser, J., Jankowski, T., & Wurm, G. 2016, ApJ, 818, 16 Google Scholar
Nealon, R., Price, D. J., & Nixon, C. J. 2015, MNRAS, 448, 1526 10.1093/mnras/stv014 2015MNRAS.448.1526N Google Scholar
Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610 10.1093/mnras/stt1475 2013MNRAS.435.2610N Google Scholar
Nickolls, J., Buck, I., Garland, M., & Skadron, K. 2008, Queue, 6, 40 10.1145/1365490.1365500 Google Scholar
Nixon, C., King, A., & Price, D. 2013, MNRAS, 434, 1946 10.1093/mnras/stt1136 2013MNRAS.434.1946N Google Scholar
Nixon, C., King, A., Price, D., & Frank, J. 2012, ApJ, 757, L24 10.1088/2041-8205/757/2/L24 2012ApJ. . .757L..24N Google Scholar
Nixon, C. J., & King, A. R. 2012, MNRAS, 421, 1201 10.1111/j.1365-2966.2011.20377.x 2012MNRAS.421.1201N Google Scholar
Öberg, K. I., Furuya, K., Loomis, R., Aikawa, Y., Andrews, S. M., Qi, C., van Dishoeck, E. F., & Wilner, D. J. 2015b, ApJ, 810, 112 10.1088/0004-637X/810/2/112 2015ApJ. . .810..112O Google Scholar
Öberg, K. I., Guzmán, V. V., Furuya, K., Qi, C., Aikawa, Y., Andrews, S. M., Loomis, R., & Wilner, D. J. 2015a, Nature, 520, 198 10.1038/nature14276 2015Natur.520..198O Google Scholar
Ogilvie, G. I. 1999, MNRAS, 304, 557 j.1365-8711.1999.02340.x 1999MNRAS.304..557O Google Scholar
Ormel, C. W., & Cuzzi, J. N. 2007, A&A, 466, 413 10.1051/0004-6361:20066899 2007A&A. . .466..413O Google Scholar
Ormel, C. W., & Klahr, H. H. 2010, A&A, 520, A43 10.1051/0004-6361/201014903 2010A&A. . .520A..43O Google Scholar
Osterbrock, D. E., & Ferland, G. J. 2006, Astrophysics of gaseous nebulae and active galactic nuclei, ed. Osterbrock, D. E. & Ferland, G. J. (Sausalito: University Science Books) . . . .O Google Scholar
Owen, J. E. 2016, PASA, 33, e005 10.1017/pasa.2016.2 2016PASA. . .33. . ..5O Google Scholar
Owen, J. E., Clarke, C. J., & Ercolano, B. 2012, MNRAS, 422, 1880 10.1111/j.1365-2966.2011.20337.x 2012MNRAS.422.1880O Google Scholar
Owen, J. E., Ercolano, B., & Clarke, C. J. 2011, MNRAS, 412, 13 10.1111/j.1365-2966.2010.17818.x 2011MNRAS.412. . .13O Google Scholar
Owen, J. E., Ercolano, B., & Clarke, C. J. 2014, Ap&SS, 36, 127 10.1007/978-3-319-03041-8_23 2014ASSP. . .36..127O Google Scholar
Owen, J. E., Ercolano, B., Clarke, C. J., & Alexander, R. D. 2010, MNRAS, 401, 1415 10.1111/j.1365-2966.2009.15771.x 2010MNRAS.401.1415O Google Scholar
Paardekooper, S.-J. 2007, A&A, 462, 355 10.1051/0004-6361:20066326 2007A&A. . .462..355P Google Scholar
Paardekooper, S.-J. 2012, MNRAS, 421, 3286 10.1111/j.1365-2966.2012.20553.x Google Scholar
Paardekooper, S.-J., & Mellema, G. 2004, A&A, 425, L9 10.1051/0004-6361:200400053 2004A&A. . .425L. . .9P Google Scholar
Papaloizou, J. C. B., & Terquem, C. 2006, RPPh, 69, 119 10.1088/0034-4885/69/1/R03 2006RPPh. . .69..119P Google Scholar
Pinilla, P., de Juan Ovelar, M., Ataiee, S., Benisty, M., Birnstiel T., van Dishoeck, E. F., & Min, M. 2015, A&A, 573, A9 10.1051/0004-6361/201424679 2015A&A. . .573A. . .9P Google Scholar
Pinilla, P., Klarmann, L., Birnstiel, T., Benisty, M., Dominik, C., & Dullemond, C. P. 2016, A&A, 585, A35 10.1051/0004-6361/201527131 2016A&A. . .585A..35P Google Scholar
Pinte, C., Harries, T. J., Min, M., Watson, A. M., Dullemond, C. P., Woitke, P., Ménard, F., & Durán-Rojas, M. C. 2009, A&A, 498, 967 10.1051/0004-6361/200811555 2009A&A. . .498..967P Google Scholar
Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797 10.1051/0004-6361:20053275 2006A&A. . .459..797P Google Scholar
Price, D. J. 2012, JCoPh, 231, 759 10.1016/ 2012JCoPh.231..759P Google Scholar
Price, D. J., & Bate, M. R. 2007, MNRAS, 377, 77 10.1111/j.1365-2966.2007.11621.x 2007MNRAS.377. . .77P Google Scholar
Price, D. J., & Federrath, C. 2010, MNRAS, 406, 1659 10.1111/j.1365-2966.2010.16810.x 2010MNRAS.406.1659P Google Scholar
Price, D. J., & Laibe, G. 2015, MNRAS, 451, 813 10.1093/mnras/stv996 2015MNRAS.451..813P Google Scholar
Ramsey, J. P., & Dullemond, C. P. 2015, A&A, 574, A81 10.1051/0004-6361/201424954 2015A&A. . .574A..81R Google Scholar
Reyes-Ruiz, M., Pérez-Tijerina, E., & Sánchez-Salcedo, F. J. 2003, in Revista Mexicana de Astronomia y Astrofisica Conf. Ser., Vol. 18, eds. Reyes-Ruiz, M. & Vázquez-Semadeni, E. (Mexico City: Instituto de Astronomía, Universidad Nacional Autónoma de México), 926 2003RMxAC..18. . .92R Google Scholar
Rice, W. K. M., Forgan, D. H., & Armitage, P. J. 2012, MNRAS, 420, 1640 10.1111/j.1365-2966.2011.20153.x Google Scholar
Rice, W. K. M., Paardekooper, S.-J., Forgan, D. H., & Armitage, P. J. 2014, MNRAS, 438, 1593 10.1093/mnras/stt2297 Google Scholar
Richard, S., Nelson, R. P., & Umurhan, O. M. 2016, MNRAS, 456, 3571 10.1093/mnras/stv2898 2016MNRAS.456.3571R Google Scholar
Robitaille, T. P. 2011, A&A, 536, A79 10.1051/0004-6361/201117150 2011A&A. . .536A..79R Google Scholar
Röllig, M., et al. 2007, A&A, 467, 187 10.1051/0004-6361:20065918 2007A&A. . .467..187R Google Scholar
Rosotti, G. P., Dale, J. E., de Juan Ovelar, M., Hubber, D. A., Kruijssen, J. M. D., Ercolano, B., & Walch, S. 2014, MNRAS, 441, 2094 10.1093/mnras/stu679 2014MNRAS.441.2094R Google Scholar
Rosotti, G. P., Ercolano, B., & Owen, J. E. 2015, MNRAS, 454, 2173 10.1093/mnras/stv2102 2015MNRAS.454.2173R Google Scholar
Rosotti, G. P., Juhasz, A., Booth, R. A., & Clarke, C. J. 2016, MNRAS, 459, 2790 10.1093/mnras/stw691 2016MNRAS.459.2790R Google Scholar
Rothstein, D. M., & Lovelace, R. V. E. 2008, ApJ, 677, 1221 10.1086/529128 2008ApJ. . .677.1221R Google Scholar
Rybicki, G. B., & Lightman, A. P. 1979, Radiative processes in astrophysics (New York: Wiley-Interscience) . . . .R Google Scholar
Sano, T., & Stone, J. M. 2002, ApJ, 570, 314 10.1086/339504 2002ApJ. . .570..314S Google Scholar
Schaal, K., Bauer, A., Chandrashekar, P., Pakmor, R., Klingenberg, C., & Springel, V. 2015, MNRAS, 453, 4278 10.1093/mnras/stv1859 2015MNRAS.453.4278S Google Scholar
Schirrmacher, V., Woitke, P., & Sedlmayr, E. 2003, A&A, 404, 267 10.1051/0004-6361:20030444 2003A&A. . .404..267S Google Scholar
Schive, H.-Y., Tsai, Y.-C., & Chiueh, T. 2010, ApJS, 186, 457 10.1088/0067-0049/186/2/457 2010ApJS..186..457S Google Scholar
Simon, J. B., Hughes, A. M., Flaherty, K. M., Bai, X.-N., & Armitage, P.J. 2015a, ApJ, 808, 180 2015ApJ. . .808..180S Google Scholar
Simon, J. B., Lesur, G., Kunz, M. W., & Armitage, P. J. 2015b, MNRAS, 454, 1117 10.1093/mnras/stv2070 2015MNRAS.454.1117S Google Scholar
Smith, N., Bally, J., & Morse, J. A. 2003, ApJ, 587, L105 10.1086/375312 2003ApJ. . .587L.105S Google Scholar
Spitzer, L. 1978, Physical processes in the interstellar medium, (New York: Wiley-Interscience) . . . .S Google Scholar
Springel, V. 2010, MNRAS, 401, 791 10.1111/j.1365-2966.2009.15715.x 2010MNRAS.401..791S Google Scholar
Springel, V. 2011, preprint (arXiv:1109.2218)2011arXiv1109.2218S Google Scholar
Staff, J. E., Koning, N., Ouyed, R., Tanaka, K., & Tan, J. C. 2016, in AAS Meeting Abstracts, 227th Meeting of the AAS, Kissimmee, 4–8 January, 236.072016AAS. . .22723607S Google Scholar
Stapelfeldt, K. R., Krist, J. E., Ménard, F., Bouvier, J., Padgett, D. L., & Burrows, C. J. 1998, ApJ, 502, L65 10.1086/311479 1998ApJ. . .502L..65S Google Scholar
Stolker, T., et al. 2016, preprint (arXiv:1603.00481)2016arXiv160300481S Google Scholar
Stone, J. E., Gohara, D., & Shi, G. 2010, CSE, 12, 66 Google Scholar
Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137 10.1086/588755 2008ApJS..178..137S Google Scholar
Suzuki, T. K., & Inutsuka, S.-i. 2009, ApJ, 691, L49 10.1088/0004-637X/691/1/L49 2009ApJ. . .691L..49S Google Scholar
Suzuki, T. K., & Inutsuka, S.-i. 2014, ApJ, 784, 121 10.1088/0004-637X/784/2/121 2014ApJ. . .784..121S Google Scholar
Suzuki, T. K., Muto, T., & Inutsuka, S.-i. 2010, ApJ, 718, 1289 10.1088/0004-637X/718/2/1289 2010ApJ. . .718.1289S Google Scholar
Tagger, M. 2001, A&A, 380, 750 10.1051/0004-6361:20011423 2001A&A. . .380..750T Google Scholar
Takahashi, S. Z., Tsukamoto, Y., & Inutsuka, S. 2016, MNRAS, 458, 3597 10.1093/mnras/stw557 2016MNRAS.458.3597T Google Scholar
Takeuchi, T., Clarke, C. J., & Lin, D. N. C. 2005, ApJ, 627, 286 10.1086/430393 2005ApJ. . .627..286T Google Scholar
Takeuchi, T., & Okuzumi, S. 2014, ApJ, 797, 132 10.1088/0004-637X/797/2/132 2014ApJ. . .797..132T Google Scholar
Taki, T., Fujimoto, M., & Ida, S. 2016, A&A, 591, A86 10.1051/0004-6361/201527732 2016A&A. . .591A..86T Google Scholar
Tasker, E. J., Brunino, R., Mitchell, N. L., Michielsen, D., Hopton, S., Pearce, F. R., Bryan, G. L., & Theuns, T. 2008, MNRAS, 390, 1267 10.1111/j.1365-2966.2008.13836.x 2008MNRAS.390.1267T Google Scholar
Teague, R., et al. 2016, A&A, 592, A49 10.1051/0004-6361/201628550 2016A&A. . .592A..49T Google Scholar
Testi, L., et al. 2014, Protostars and Planets VI, 33910.2458/azu_uapress_9780816531240-ch015 2014prpl.conf..339T Google Scholar
Testi, L., et al. 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14), 117 2015aska.confE.117T Google Scholar
Thiabaud, A., Marboeuf, U., Alibert, Y., Leya, I., & Mezger, K. 2015, A&A, 574, A138 10.1051/0004-6361/201424868 2015A&A. . .574A.138T Google Scholar
Thilliez, E., & Maddison, S. T. 2015, PASA, 32, e039 10.1017/pasa.2015.41 2015PASA. . .32. . .39T Google Scholar
Tomida, K., Okuzumi, S., & Machida, M. N. 2015, ApJ, 801, 117 10.1088/0004-637X/801/2/117 2015ApJ. . .801..117T Google Scholar
Toro, E. F. 2013, Riemann solvers and numerical methods for fluid dynamics: a practical introduction (Berlin: Springer Science & Business Media)Google Scholar
Tricco, T. S., & Price, D. J. 2012, JCoPh, 231, 7214 10.1016/ 2012JCoPh.231.7214T Google Scholar
Tsukamoto, Y. 2016, PASA, 33, e010 10.1017/pasa.2016.6 2016PASA. . .33. . .10T Google Scholar
Tsukamoto, Y., Iwasaki, K., Okuzumi, S., Machida, M. N., & Inutsuka, S. 2015, ApJ, 810, L26 10.1088/2041-8205/810/2/L26 2015ApJ. . .810L..26T Google Scholar
van der Marel, N., van Dishoeck, E. F., Bruderer, S., Andrews, S. M., Pontoppidan, K. M., Herczeg, G. J., van Kempen, T., & Miotello, A. 2016, A&A, 585, A58 10.1051/0004-6361/201526988 2016A&A. . .585A..58V Google Scholar
van der Marel, N., et al. 2013, Science, 340, 1199 10.1126/science.1236770 2013Sci. . .340.1199V Google Scholar
Vincke, K., Breslau, A., & Pfalzner, S. 2015, A&A, 577, A115 10.1051/0004-6361/201425552 2015A&A. . .577A.115V Google Scholar
Vincke, K., & Pfalzner, S. 2016, ApJ, 828, 48 2016ApJ. . .828. . .48V Google Scholar
Viti, S., Collings, M. P., Dever, J. W., McCoustra, M. R. S., & Williams, D. A. 2004, MNRAS, 354, 1141 10.1111/j.1365-2966.2004.08273.x 2004MNRAS.354.1141V Google Scholar
Viti, S., Jimenez-Serra, I., Yates, J. A., Codella, C., Vasta, M., Caselli, P., Lefloch, B., & Ceccarelli, C. 2011, ApJ, 740, L3 10.1088/2041-8205/740/1/L3 2011ApJ. . .740L. . .3V Google Scholar
Vorobyov, E. I., Lin, D. N. C., & Guedel, M. 2015, A&A, 573, A5 10.1051/0004-6361/201424583 2015A&A. . .573A. . .5V Google Scholar
Wagner, K., Apai, D., Kasper, M., & Robberto, M. 2015, ApJ, 813, L2 10.1088/2041-8205/813/1/L2 2015ApJ. . .813L. . .2W Google Scholar
Wakelam, V., et al. 2012, ApJS, 199, 21 10.1088/0067-0049/199/1/21 2012ApJS..199. . .21W Google Scholar
Walsh, C., Nomura, H., Millar, T. J., & Aikawa, Y. 2012, ApJ, 747, 114 10.1088/0004-637X/747/2/114 2012ApJ. . .747..114W Google Scholar
Walsh, C., Nomura, H., & van Dishoeck, E. 2015, A&A, 582, A88 10.1051/0004-6361/201526751 2015A&A. . .582A..88W Google Scholar
Whitehouse, S. C., & Bate, M. R. 2004, MNRAS, 353, 1078 10.1111/j.1365-2966.2004.08131.x 2004MNRAS.353.1078W Google Scholar
Whitehouse, S. C., Bate, M. R., & Monaghan, J. J. 2005, MNRAS, 364, 1367 10.1111/j.1365-2966.2005.09683.x 2005MNRAS.364.1367W Google Scholar
Williams, J. P., & Best, W. M. J. 2014, ApJ, 788, 59 10.1088/0004-637X/788/1/59 2014ApJ. . .788. . .59W Google Scholar
Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67 10.1146/annurev-astro-081710-102548 2011ARA&A..49. . .67W Google Scholar
Williams, J. P., & McPartland, C. 2016, preprint (arXiv:1606.05646)2016arXiv160605646W Google Scholar
Winn, J. N., & Fabrycky, D. C. 2015, ARA&A, 53, 409 10.1146/annurev-astro-082214-122246 2015ARA&A..53..409W Google Scholar
Woitke, P. 2015, in EPJ Web of Conferences, Summer School - Protoplanetary Disks: Theory and Modeling Meet Observations, ed. Kamp, I., Woitke, P., & Ilee, J. D. (Ameland, The Netherlands: EDP Sciences), 00011 10.1051/epjconf/201510200011 Google Scholar
Woitke, P., Goeres, A., & Sedlmayr, E. 1996b, A&A, 313, 217 1996A&A. . .313..217W Google Scholar
Woitke, P., Kamp, I., & Thi, W.-F. 2009, A&A, 501, 383 10.1051/0004-6361/200911821 2009A&A. . .501..383W Google Scholar
Woitke, P., Krueger, D., & Sedlmayr, E. 1996a, A&A, 311, 927 1996A&A. . .311..927W Google Scholar
Woitke, P., et al. 2016, A&A, 586, A103 10.1051/0004-6361/201526538 2016A&A. . .586A.103W Google Scholar
Woodall, J., Agúndez, M., Markwick-Kemper, A. J., & Millar, T. J. 2007, A&A, 466, 1197 10.1051/0004-6361:20064981 2007A&A. . .466.1197W Google Scholar
Wright, N. J., Drake, J. J., Drew, J. E., Guarcello, M. G., Gutermuth, R. A., Hora, J. L., & Kraemer, K. E. 2012, ApJ, 746, L21 10.1088/2041-8205/746/2/L21 2012ApJ. . .746L..21W Google Scholar
Wurster, J., Price, D. J., & Bate, M. R. 2016, MNRAS, 457, 1037 10.1093/mnras/stw013 2016MNRAS.457.1037W Google Scholar
Yang, C.-C., & Johansen, A. 2016, ApJS, 224, 39 10.3847/0067-0049/224/2/39 2016ApJS..224. . .39Y Google Scholar
Youdin, A., & Johansen, A. 2007, ApJ, 662, 613 10.1086/516729 2007ApJ. . .662..613Y Google Scholar
Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 10.1086/426895 2005ApJ. . .620..459Y Google Scholar
Young, M. D., & Clarke, C. J. 2015, MNRAS, 451, 3987 10.1093/mnras/stv1266 2015MNRAS.451.3987Y Google Scholar
Young, M. D., & Clarke, C. J. 2016, MNRAS, 455, 1438 10.1093/mnras/stv2378 Google Scholar
Yuan, C., & Fox, R. 2011, JCoPh, 230, 8216 10.1016/ Google Scholar
Yuan, C., Laurent, F., & Fox, R. 2012, JAerS, 51, 1 10.1016/j.jaerosci.2012.04.003 Google Scholar
Zhu, Z., Nelson, R. P., Hartmann, L., Espaillat, C., & Calvet, N. 2011, ApJ, 729, 47 10.1088/0004-637X/729/1/47 2011ApJ. . .729. . .47Z Google Scholar
Figure 0

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.

Figure 1

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.

Figure 2

Table 1. Common gas–grain reactions in astrophysical environments. Species are all considered to be in the gas phase, unless shown as Xgr, 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).

Figure 3

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 H2CO, 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).

Figure 4

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.

Figure 5

Table 3. Participants of protoplanetary discussions 2016.