<figure><img src="image_1.png" /></figure>

HAL
open science

# Two-fluid dusty shocks: simple benchmarking problems and applications to protoplanetary discs

Andrew Lehmann, Mark Wardle

► To cite this version:

Andrew Lehmann, Mark Wardle. Two-fluid dusty shocks: simple benchmarking problems and applications to protoplanetary discs. Monthly Notices of the Royal Astronomical Society, 2018, 476 (3), pp.3185-3194. 10.1093/mnras/sty450. hal-02345624

HAL Id: hal-02345624

https://hal.science/hal-02345624

Submitted on 15 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Two-fluid dusty shocks: simple benchmarking problems and applications to protoplanetary discs

Andrew Lehmann<fnref n="1" />,<fnref n="2" />* and Mark Wardle<fnref n="1" />

Accepted 2018 February 15. Received 2018 February 15; in original form 2016 December 29

## ABSTRACT

The key role that dust plays in the interstellar medium has motivated the development of numerical codes designed to study the coupled evolution of dust and gas in systems such as turbulent molecular clouds and protoplanetary discs. Drift between dust and gas has proven to be important as well as numerically challenging. We provide simple benchmarking problems for dusty gas codes by numerically solving the two-fluid dust–gas equations for steady, plane-parallel shock waves. The two distinct shock solutions to these equations allow a numerical code to test different forms of drag between the two fluids, the strength of that drag and the dust to gas ratio. We also provide an astrophysical application of J-type dust–gas shocks to studying the structure of accretion shocks on to protoplanetary discs. We find that two-fluid effects are most important for grains larger than 1 µm, and that the peak dust temperature within an accretion shock provides a signature of the dust-to-gas ratio of the infalling material.

**Key words:** shock waves – protoplanetary disc.

## 1 INTRODUCTION

Dust plays a key role in numerous astrophysical environments. In the interstellar medium (ISM), dust is heavily involved in controlling the thermodynamics by being a major coolant, collisional partner, and source of opacity. It allows us to probe the magnetic field by measuring the polarization in thermal dust emission (Planck Collaboration XXXIII 2016), and is crucial to the formation of H₂ in molecular clouds by providing a catalytic surface and by attenuating the dissociating ultraviolet radiation field (Glover & Clark 2012). Thermal dust emission is observed with telescopes such as *Spitzer* (e.g. Stephens et al. 2014), *Herschel* (e.g. Launhardt et al. 2013), and ALMA (e.g. ALMA Partnership 2015), and these observations are used to obtain properties of the gas. Thus it is crucially important to understand not only the properties of dust grains, but also their coupled evolution with the gas phase in order to rigorously relate dust emission to properties of the ISM. In protoplanetary discs, the radial pressure gradient causes gas to flow at below the Keplerian velocity. The pressureless dust component is closer to Keplerian, so the two fluids drift with respect to each other and accurate modelling requires a careful treatment of the drag.

The importance of coupled gas–dust modelling has motivated the development of numerical codes designed to simulate gas and dust in various astrophysical systems. For example, radial gradients of

gas pressure in protoplanetary discs induce dust clumping leading to planetesimal formation (Bai & Stone 2010b), large dust-to-gas variations occur in turbulent molecular clouds and dust filaments do not necessarily correlate with gas filaments (Hopkins & Lee 2016), and dust-gaps are cleared more easily than gas-gaps in protoplanetary discs (Paardekooper & Mellema 2006). Both grid-based (Johansen, Henning & Klahr 2006; Bai & Stone 2010a) and particle-based smoothed particle hydrodynamics codes (Laibe & Price 2012a; Lorén-Aguilar & Bate 2014) have been used to study dusty gas flows in the ISM. However, Laibe & Price (2011) highlight a lack of simple analytic solutions to benchmark dusty gas codes in astrophysical conditions. The main goal of this work is to provide such a simple solution by computing the structure of steady-state, planar two-fluid dusty gas shock waves. Unlike the standard shock-tube tests, steady shocks comprise only one hydrodynamic component with a structure that can be computed by simply integrating the governing ordinary differential equations, as we do in Section 3. The numerical simulation is also simple: drive a piston represented by reflective boundary conditions into a uniform medium. This simple test can be used to benchmark how numerical codes behave with different dust-to-gas ratios, or e.g. linear, quadratic, or Epstein forms of the drag coefficients.

Two-fluid dusty shocks are not just ideal benchmarks for numerical codes. Supersonic flows occur ubiquitously in astrophysical systems. For example, in the inside-out collapse model of protostellar cores, material becomes thermally unsupported and free-falls on to the protoplanetary disc at a few km s⁻<fnref n="1" />. The sound speed in the

<footnote n="1">¹Department of Physics and Astronomy, and Research Centre for Astronomy, Astrophysics & Astrophotonics, Macquarie University, Sydney, NSW 2109, Australia</footnote>

<footnote n="2">²Sorbonne Université, Observatoire de Paris, Université PSL, École normale supérieure, CNRS, LERMA, F-75005 Paris, France</footnote>

<footnote>* E-mail: andrew.lehmann@ens.fr</footnote>
gas is only ~0.2 km s⁻¹, and so a shock wave forms as the material decelerates to settle on to the disc. In Section 4 we provide an application of our two-fluid shock solutions to study this type of accretion shock.

## 2 THEORY

In this section we outline the set of equations that describe the two-fluid dust-gas system. We use these equations to derive the dispersion relation for linear waves in the combined fluid, and discuss how this affects the possible dust-gas shock structures. We characterize the initial stationary states to outline the criteria for J- and C-type shocks to occur and discuss their structure.

### 2.1 Fluid equations

For a fluid with gas density $ρ$, velocity $v$, and pressure $P$ the equations of continuity and conservation of momentum can be written as

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0, \qquad (1) $$

$$ \rho \frac{\partial \mathbf{v}}{\partial t} + \rho (\mathbf{v} \cdot \nabla) \mathbf{v} = -\nabla P - \mathbf{F}_{\text{drag}}, \qquad (2) $$

where $F_{\text{drag}}$ is the rate at which momentum is added to the gas via drag from the dust fluid. The energy equation can be written as

$$ \frac{\partial}{\partial t} \left( \frac{1}{2} \rho v^2 + u + P \right) + \nabla \cdot \left( \left( \frac{1}{2} \rho v^2 + u + P \right) \mathbf{v} \right) = \Gamma - \Lambda, \quad (3) $$

where $u$ is the internal energy per unit volume, $\Gamma$ is the heating rate per volume, and $\Lambda$ is the cooling rate per volume. The analogous equations for the coextensive dust fluid with density $\rho_d$ and velocity $v_d$ are

$$ \frac{\partial \rho_d}{\partial t} + \nabla \cdot (\rho_d \mathbf{v}_d) = 0, \qquad (4) $$

$$ \rho_d \frac{\partial v_d}{\partial t} + \rho_d (v_d \cdot \nabla) v_d = F_{\text{drag}}, \qquad (5) $$

where we have assumed the dust to be pressureless. Finally, we use the ideal equation of state

$$ P = \frac{\rho k_B T}{\mu}, \qquad (6) $$

where $\mu$ is the mean mass per gas particle, $T$ is the gas temperature, and $k_B$ is the Boltzmann constant.

In Section 4, we introduce heating and cooling functions appropriate to astrophysical applications so that the gas temperature varies. Here, for the purposes of computing simple benchmarking problems, we assume the gas is isothermal so that equation (3) is ignored, and the ideal gas law becomes

$$ P = \rho c_s^2, \qquad (7) $$

where $c_s$ is the isothermal sound speed.

The form of the drag term has been thoroughly discussed by Laibe & Price (2012b) for various regimes of astrophysical interest. It is generally proportional to a power law of the drift velocity between the two fluids. If we consider the linear drag regime, the drag term on the total fluid can be written as

$$ F_{\text{drag}} = K (\mathbf{v} - \mathbf{v}_d) \qquad (8) $$

for drag coefficient $K$. A linear analysis gives the dispersion relation for waves with angular frequency $\omega$ and wavenumber $k$

$$ (\omega^2 - k^2 c_s^2) + \frac{i}{\omega \tau_s} (\omega^2 - k^2 \tilde{c}_s^2) = 0, \qquad (9) $$

where

$$ \tau_s = \frac{\rho_0 \rho_{d0}}{K (\rho_0 + \rho_{d0})} $$

is the drag stopping time – the characteristic time to damp the differential velocity between the dust and gas fluids – for unperturbed gas and dust density $\rho_0$ and $\rho_{d0}$, respectively. The combined fluid sound speed

$$ \tilde{c}_s = c_s (1 + D)^{-1/2}, \qquad (10) $$

with the dust-to-gas ratio

$$ D = \frac{\rho_{d0}}{\rho_0}. \qquad (11) $$

In the limit of weak coupling between the dust and gas ($K \to 0$, $\tau_s \to \infty$) we recover the dispersion relation for ordinary sound waves in a gas with phase velocity $\omega/k = c_s$. In the strong coupling limit ($K \to \infty$, $\tau_s \to 0$) the second term of equation (9) dominates, and so waves travel at the combined sound speed $\tilde{c}_s$.

The two signal speeds in the system, $c_s$ and $\tilde{c}_s$, determine the possible structures of dust-gas shocks. The dust-gas mixture will behave as a single fluid far ahead and behind a shock and so the combined sound speed $\tilde{c}_s$ is the relevant signal speed that, in the frame of reference comoving with the shock, the fluid velocity must transition across. As $\tilde{c}_s$ is necessarily less than $c_s$, we will see that two distinct classes of shocks arise depending on whether the shock speed is greater or less than the gas signal speed $c_s$.

For a supersonic shock (shock velocity $v_s > c_s$) the pre-shock fluid is overrun by high density gas in a thin shock front a few mean free paths wide that resembles an ordinary gas dynamic shock. The dust particles cannot respond quickly, and so there is a relaxation zone wherein the dust particles are accelerated until the two fluids flow at the same velocity. This structure is qualitatively sketched in Fig. 1.

When the shock speed is between the two signal speeds, sound waves in the gas fluid can travel ahead of the shock front and compress the gas and dust in such a way that all the fluid variables remain continuous through the shock. We will call these two classes J-type and C-type shocks in analogy to the kinds of magnetized two-fluid shocks outlined by Draine (1986), where ion-magnetosonic waves can travel ahead of a jump front to form a 'magnetic precursor'. In the next section we will characterize these classes in further detail.

### 2.2 Shock classification

Assuming a steady state, one-dimensional structure varying in the z-direction and a power law drag term with index $\beta$, the gas fluid equations (1) and (2) in the frame of reference comoving with the shock reduce to

$$ \frac{d}{dz} (\rho v) = 0, \qquad (12) $$

$$ \frac{d}{dz} (\rho v^2 + P) = K |v_d - v|^\beta, \qquad (13) $$

and the dust fluid equations (4) and (5) reduce to

$$ \frac{d}{dz} (\rho_d v_d) = 0, \qquad (14) $$
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Sketch of a J-type dust–gas two-fluid shock in the frame of reference where the pre-shock fluid is stationary. The red, smaller circles refer to gas particles whereas the larger black circles refer to dust particles. The black dashed line marks the jump transition in the gas fluid that takes place of the order of a few mean free paths.</figcaption></figure>

$$ \frac{d}{dz} (\rho_d v_d^2) = -K |v_d - v|^\beta. \qquad (15) $$

The pre-shock state, defined by the gas density $\rho_0$, dust density $\rho_{d0}$, shock velocity $v_s$, and initial pressure $P_0$ is a stationary point. In order to classify this state, we perturb the initial state and let the perturbation grow (or decay) exponentially as follows:

$$ v = v_s + \delta v e^{\lambda z}, \qquad (16) $$

$$ v_d = v_s + \delta v_d e^{\lambda z}, \qquad (17) $$

$$ \rho = \rho_0 + \delta \rho e^{\lambda z}, \qquad (18) $$

$$ \rho_d = \rho_{d0} + \delta \rho_d e^{\lambda z}, \qquad (19) $$

$$ P = P_0 + \delta P e^{\lambda z} = \rho_0 c_s^2 + c_s^2 \delta \rho e^{\lambda z}. \qquad (20) $$

Substituting these perturbed variables back into equations (12)–(15), assuming linear drag ($\beta = 1$), ignoring terms of second order in perturbations and solving for $\lambda$ gives the eigenvalue

$$ \lambda = - \frac{K}{\rho_{d0} v_s} \frac{1 + D}{D} \frac{v_s^2 - \tilde{c}_s^2}{v_s^2 - c_s^2}. \qquad (21) $$

This expression implies that for supersonic shocks ($v_s > c_s > \tilde{c}_s$), the eigenvalue $\lambda < 0$ and hence the perturbation decays. That is, the initial state is a stable stationary point. Hence it requires an initial discontinuity (jump) to get across the sound speed. In this kind of shock, the gas fluid is highly compressed over a few mean free paths. The dust cannot respond quickly, and so this jump is determined by the hydrodynamic jump conditions.

By replacing $v_s$ with the post-shock solution $v_{post}$ we can explore the final state. A shock in the total gas–dust fluid is a transition across the total speed $\tilde{c}_s$, and so $v_{post} < \tilde{c}_s$. Thus the eigenvalue near the post-shock state $\lambda < 0$, which defines a stable point, so jumping near this state will settle on to it. Two-fluid dust–gas J-type shocks

were first discussed by Carrier (1958), and have been thoroughly studied for various non-astrophysical applications (see review by Igra & Ben-Dor 1988, and references therein).

For C-type shock solutions to exist we require a positive eigenvalue in equation (21), so that perturbations smoothly grow away from the pre-shock state. The shock speed must then be in the range

$$ \tilde{c}_s < v_s < c_s \qquad (22) $$

and hence the Mach number is necessarily below unity for this type of shock. Without cooling, this shock will smoothly settle on to the post-shock state. However, with cooling there is the potential for the gas sound speed to drop (as the gas compresses in the shock) faster than the velocity. If $c_s = v$ in the shock somewhere, a jump will be required. Otherwise, the fluid variables will remain continuous throughout the shock. This kind of shock was first investigated by Kriebel (1964) and further developed by Miura (1972). In the next section, we discuss how these two kinds of shocks could be used to benchmark dusty gas numerical codes.

## 3 BENCHMARK PROBLEMS

In this section we describe the first order differential equation that we solve to investigate the structure of two-fluid dust–gas shocks. These isothermal shock solutions are ideal tests for benchmarking numerical codes that wish to simulate dusty gas. Our `PYTHON` code that returns the shock solutions described in this section is publicly available on the Python Package Index<fnref n="1" /> and `BitBucket`<fnref n="2" />.

In the standard shock-tube problem (Sod 1978) the simple setup breaks up into a shock wave, rarefaction wave and a contact discontinuity. In the two-fluid dust–gas version of this problem there is no known analytic solution, but Saito, Marumoto & Takayama (2003) find that a steady-state shock solution fits one of the components well. Unlike the Sod shock-tube, the steady shocks we compute here are very simple, comprising just one hydrodynamic structure.

<footnote n="1">¹ https://pypi.python.org/pypi/DustyShock</footnote>

<footnote n="2">² https://bitbucket.org/AndrewLehmann/dustyshock</footnote>
A numerical code can then be tested against the steady solution by setting up a reflective boundary representing a piston as the driver of the shock, as in Toth (1994).

## 3.1 Numerical integration

We get the dimensionless derivative of the dust velocity from combining equations (14) and (15):

$$ \frac{dw_d}{d\xi} = -|w_d - w|^\beta \quad (23) $$

for the normalized velocities and position defined by

$$ w_d = \frac{v_d}{v_s}, \quad (24) $$

$$ w = \frac{v}{v_s}, \quad (25) $$

$$ \xi = \frac{K}{\rho_{d0} v_s^{2-\beta}} z, \quad (26) $$

where $v_s$ is the shock velocity and $\rho_{d0}$ is the pre-shock dust mass density. Combining equations (13) and (15), and making use of isothermality, gives

$$ \frac{d}{dz} (\rho v^2 + \rho c_s^2 + \rho_d v_d^2) = 0, \quad (27) $$

which states that the sum of the gas and dust ram pressures and the thermal pressure is conserved through the shock. The sum of pressures at any point in the shock retains the pre-shock value, so that

$$ \rho v^2 + \rho c_s^2 + \rho_d v_d^2 = \rho_0 v_s^2 + \rho_0 c_s^2 + \rho_{d0} v_s^2. \quad (28) $$

From this equation, we derive a quadratic equation in the gas velocity

$$ w^2 + [D(w_d - 1) - 1 - M^{-2}] w + M^{-2} = 0, \quad (29) $$

where we have used the sonic Mach number $M \equiv v_s/c_s$, assuming that initially the gas and dust flow together at the shock velocity $v_s$. The two roots of the quadratic represent supersonic and subsonic (with regards to $c_s$) gas velocities, and closes equation (23).

The fluid variables defining the shock structure are obtained by integrating the first-order ordinary differential equation (ODE) defined by equation (23). In this paper we use the open source PYTHON module scikits.odes.<fnref n="3" /> This module solves initial value problems for ODEs using variable-order, variable-step, multistep methods.

### 3.1.1 J-type shock solutions

For J-type shocks, the pre-shock fluid is in a stable stationary state with gas fluid velocity equal to the supersonic root of equation (29). As discussed in Section 2.2, a hydrodynamic jump in the gas fluid is required to cross the sound speed $c_s$. For the isothermal shocks computed here, the density immediately jumps from $\rho_0$ to

$$ \rho_j = M^2 \rho_0. \quad (30) $$

The velocity switches from the supersonic to the subsonic root of equation (29). Integration of equation (23) then gives the post-jump structure.

An example of the velocity structure of an isothermal J-type shock in the frame of reference comoving with the shock is shown in Fig. 2. The upper panel shows shock structures computed with Mach

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** J-type gas–dust shocks with Mach number $M = 2$ and initial dust-to-gas ratios $D = 1$ (solid lines), $D = 0.1$ (dashed lines) and $D = 0.01$ (dotted lines). The red lines give the gas velocity and the black lines give the dust velocity, both normalized to the shock velocity. The upper panel are solutions computed with a linear drag term while the lower panel uses a quadratic drag term. Note that the z-scale normalization differs according to equation (26).</figcaption></figure>

number $M = 2.0$ using a linear drag term ($\beta = 1$) and initial dust-to-gas ratios varying from 0.01 to 1. After the initial hydrodynamic jump in the gas, the dust lags behind the gas, but both eventually settle to the same velocity (below $\tilde{c}_s$). The combined fluid velocities far ahead ($v_s$) and behind ($v_2$) the shock front are related by

$$ \frac{v_2}{v_s} = \left( \frac{\tilde{c}_s}{v_s} \right)^2 = (1 + D)^{-1} M^{-2}. \quad (31) $$

Hence reducing $D$ changes the shock structure by increasing the post-shock velocity that the solution settles to. This effect allows this steady state solution to test how a numerical code behaves with different dust-to-gas ratios.

The lower panel of Fig. 2 shows isothermal shock solutions with the same conditions as the upper panel except that they are computed using a quadratic drag term

$$ F_{\text{drag}} = K |v - v_d|^2. \quad (32) $$

The structure is qualitatively similar to the linear drag case, however the shock thickness is an order of magnitude larger (in the dimensionless position variable $\xi$). From equation (23) the shock thickness

$$ \Delta\xi \sim (\Delta w)^{1-\beta}. \quad (33) $$

<footnote n="3">³ https://github.com/bmcage/odes</footnote>
$\Delta w$ is necessarily between 0 and 1, and so the shock thickness increases with the index of the power-law drag. For this reason, the solution with quadratic drag has an extended tail of very small but finite $\Delta w$. In real units, however, the shock thickness

$$ \Delta z = \frac{\rho_{d0} v_s^{2-\beta}}{K} |\Delta w|^{1-\beta}. \quad (34) $$

If the same value of K and $\rho_{d0}$ is used in the different drag cases, the size of the shock in the linear regime ($\beta = 1$) is larger than the shock in the quadratic regime ($\beta = 2$) by the factor $v_s|\Delta w|$. For astrophysical shocks of the order of a few km s$^{-1}$, this means shocks with linear drag will be ~10$^4$ times larger than shocks with quadratic drag. Note that the drag coefficient K has different dimensions in the different drag regimes.

These stark differences allow a numerical code's implementation of different drag coefficients to be tested by the shock problem. Note that regardless of the form of the drag term, the shock solution settles on to the same post-shock velocity. This is because the jump conditions relate the velocities of the combined fluid on either side of the shock.

The different shock structures resulting from different drag terms could be used to test a numerical code's implementation of the drag. However, in many codes artificial viscosity is used to smear a discontinuity over a finite distance. The scaling factor in equation (34) can be adjusted until the relaxation tail of the J-type analytic shock solution spreads over many computational cells, so that the behaviour of a numerical code in certain regimes of the drag coefficient (small K) and/or dust density (large $\rho_{d0}$) can be tested. The C-type shock solution in the following section is intrinsically smooth and so does not require any treatment of a discontinuity, isolating the dust-gas drag as the agent responsible for mediating the shock transition. A numerical code's implementation of the dust-gas interaction is therefore isolated and tested by comparing to these C-type solutions.

### 3.1.2 C-type shock solutions

C-type shocks are necessarily subsonic, and so the shock structure can be obtained by integrating equation (23) using the subsonic root of equation (29).

An example of the velocity structure of isothermal C-type shocks is shown in Fig. 3 using a linear drag term ($\beta = 1$). The pre-shock variables are the same as for the J-type shock discussed above, except that the Mach number $M = 0.95$. Sound waves in the gas fluid stream ahead of the shock and gradually compress the pre-shock medium. As the two fluids smoothly settle on to the post-shock state the drift velocity remains small compared to the drift velocities reached in the J-type shocks. When computed with a quadratic drag term ($\beta = 2$, lower panel), the structure is an order of magnitude wider (in the dimensionless variable $\xi$) than in the linear case, explained by equation (33). Similarly to the J-type shock with quadratic drag, there is an extended tail of very small but finite $\Delta w$.

Recall that the shock velocity $v_s$ for C-type dusty gas shocks is restricted to the range

$$ \tilde{c}_s = \frac{c_s}{\sqrt{1+D}} < v_s < c_s. \quad (35) $$

This means that when the initial dust-to-gas ratio becomes small – such as the typical interstellar value of 0.01 – $v_s$ cannot be much larger than $\tilde{c}_s$, resulting in a very weak shock. For this reason, C-type dusty shocks may not be relevant to the general ISM. However, strong variations of D have been found in simulations of turbulent molecular clouds when dust-gas decoupling has been modelled

<figure><img src="image_6.png" /><figcaption>Figure 3. Same as Fig. 2 but for a C-type shock with Mach number $M = 0.95$ and $D = 1$.</figcaption></figure>

(Hopkins & Lee 2016), and values as high as unity have been used in protostellar discs (e.g. Dipierro et al. 2015) to account for dust migration to the inner parts of the disc.

We have presented two potential benchmarking problems for numerical codes seeking to simulate dusty gas systems. These problems allow a code to test the implementation of different forms of drag, the strength of the drag (K) and the dust-to-gas mass density ratio (D). In the following section we provide an astrophysical application of two-fluid dust-gas shocks.

## 4 PROTOPLANETARY DISC ACCRETION SHOCK

Here we present a crude astrophysical application of two-fluid dust-gas shocks. At very early stages of star formation, the system is characterized by an embedded protostellar disc surrounded by an infalling envelope of dust and gas. The material falls, due to gravity, through an accretion shock and then eventually settles on to the disc. Here we model this accretion shock as a two-fluid dust-gas shock. We first derive the shock parameters by considering a non-tional disc model, the minimum-mass solar nebula (Weidenschilling 1977; Hayashi 1981), and assume vertical infall of material. We then detail the astrophysical drag, heating and cooling terms. Finally, we discuss the structure of the accretion shock. While this application is highly simplified, it serves as a starting point for more detailed analyses of protoplanetary disc accretion shocks. We finish this section with a discussion of these simplifications and indicate improvements to the model.
## 4.1 Shock parameters

In this section we consider typical physical properties of protoplanetary discs in order to constrain the appropriate pre-shock conditions to model the accretion shock. We then consider the shock geometry by equating the disc thermal pressure with the shock’s ram pressure.

The accretion rate on to a protoplanetary disc in the inside-out collapse model of a singular isothermal sphere can be approximated (Larson 2003) as

$$ \dot{M} \sim \frac{c_s^3}{G} = 2 \times 10^{-6} M_{\odot} \text{yr}^{-1}, \qquad (36) $$

where G is the gravitational constant, and we have assumed the typical interstellar medium sound speed $c_s = 0.2 \text{ km s}^{-1}$. This mass is spread over the area of the disc, so that the mass flux of the accretion shock is

$$ \rho_0 v_s = \frac{\dot{M}}{\pi R_d^2} \sim \frac{c_s^3}{G \pi R_d^2} \qquad (37) $$

$$ \sim 1.7 \times 10^{-11} \left( \frac{R_d}{100 \text{ au}} \right)^{-2} \text{g s}^{-1} \text{cm}^{-2}, \qquad (38) $$

where $R_d$ is the disc radius. The material accretes in free-fall and thus reaches the protoplanetary disc at the escape speed. This free-falling gas meets a disc in Keplerian motion and so the relative velocity, giving the shock velocity, is approximately

$$ v_s \sim \sqrt{\frac{GM_*}{r}}, \qquad (39) $$

where r is the distance from the central protostar with mass $M_*$. Using this velocity and equation (38) for a solar mass protostar we get a pre-shock density of

$$ \rho_0 \sim 4 \times 10^{-17} \left( \frac{R_d}{100 \text{ au}} \right)^{-2} \left( \frac{r}{50 \text{ au}} \right)^{1/2} \text{g cm}^{-3}. \qquad (40) $$

This corresponds to a total hydrogen number density, $n_H \approx n(H) + 2n(H_2)$, of

$$ n_0 = \frac{\rho_0}{1.4 m_H} \sim 2 \times 10^7 \left( \frac{R_d}{100 \text{ au}} \right)^{-2} \left( \frac{r}{50 \text{ au}} \right)^{1/2} \text{cm}^{-3}. \qquad (41) $$

The shock will be located where its ram pressure is balanced by the thermal pressure of the disc. That is, where

$$ n_0 v_s^2 = n_{\text{disc}}(r, z) c_{\text{disc}}(r)^2, \qquad (42) $$

where the disc density and sound speed are functions of radial distance from the star, r, and vertical distance from the disc, z. These functions can be approximated in the minimum-mass solar nebula model (Wardle 2007) as

$$ n_{\text{disc}}(r, z) \sim 5.8 \times 10^{14} \text{cm}^{-3} \left(\frac{r}{\text{au}}\right)^{-11/4} \exp\left(-\frac{z^2}{2h^2}\right), \qquad (43) $$

$$ c_{\text{disc}}(r) \sim 0.99 \text{ km s}^{-1} \left(\frac{r}{\text{au}}\right)^{-1/4}, \qquad (44) $$

with the scale height, h, is given by

$$ \frac{h}{r} \sim 0.03 \left( \frac{r}{\text{au}} \right)^{1/4}. \qquad (45) $$

Substituting these approximations into equation (42) and rearranging for the vertical height gives

$$ z_s^2 = 2h^2 \ln \left( 2.62 \times 10^5 \left( \frac{R_d}{100 \text{ au}} \right)^2 \left( \frac{r}{\text{au}} \right)^{-11/4} \right). \qquad (46) $$

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Location of shock (solid line), where the ram pressure of free-falling material balances the thermal pressure of the protoplanetary disc. The dotted line marks the scale height for comparison.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Sketch of model for accretion shock (dashed) on to protoplanetary disc.</figcaption></figure>

The location where the shock ram pressure balances the disc thermal pressure ($z_s$) is shown in Fig. 4 for a disc radius $R_d = 100 \text{ au}$. Beyond $r \sim 95 \text{ au}$, the disc density and sound speed has dropped so low that its thermal pressure never balances the shock ram pressure. In this region the shock will run up against the material freely falling from the other side of the disc. A sketch of the system is shown in Fig. 5. In the next section we model this accretion shock using representative values of the shock velocity and pre-shock density between 30 and 80 au, i.e. $v_s \sim 4 \text{ km s}^{-1}$ and $\rho_0 \sim 4 \times 10^{-17} \text{ g cm}^{-3}$, respectively. These parameters could be used to model an isothermal shock. However, heating and cooling impact the shock structure, and provide observational predictions through radiation. So in the next section we derive an ODE for the temperature profile of the shock, to be solved simultaneously with equation (23), and describe the appropriate heating and cooling terms.

## 4.2 Physical processes

At 4 km s$^{-1}$, this shock is highly supersonic so we use the J-type model described in Section 2.2. In this section we outline the modifications to the fluid equations to introduce astrophysically relevant drag, heating, and cooling.

When the grain size is small compared to the surrounding gas mean free path the Epstein drag regime applies (Epstein 1924) and the drag term, or momentum rate of change per volume due to
elastic scattering, can be expressed as (Draine 1986):

$$F = \frac{\sigma \rho \rho_d}{m_d} \sqrt{\frac{2k_B T}{\pi \mu}} (v_d - v) I(v, v_d, T), \quad (47)$$

where $\sigma$ is the dust cross-section, $m_d$ is the mass of a dust particle, $\mu$ is the mean mass per gas particle ((7/3)$m_H$ in molecular gas), $T$ is the gas temperature, and $k_B$ is the Boltzmann constant. The function $I(v, v_d, T)$ is well approximated by (Kwok 1975):

$$I(v, v_d, T) \approx \frac{8}{3} \left( 1 + \frac{9\pi}{64} \frac{\frac{1}{2}\mu |v - v_d|^2}{k_B T} \right)^{1/2}. \quad (48)$$

The derivative of the dust velocity through the shock can then be written as

$$\frac{dv_d}{dz} = \alpha v_s \sqrt{\frac{k_B T}{\mu}} I(v, v_d, T) \frac{v - v_d}{v_d v}, \quad (49)$$

where

$$\alpha = \frac{\rho_0 r_d^2 \sqrt{2\pi}}{m_d} \quad (50)$$

for pre-shock gas density $\rho_0$ and dust radius $r_d$. We have assumed that the dust cross-section is $\pi r_d^2$.

The energy equation (3) reduces, in the static plane-parallel case, to

$$v \frac{dP}{dz} + \gamma P \frac{dv}{dz} = (\gamma - 1)(\Gamma - \Lambda), \quad (51)$$

where we have assumed that the internal energy is proportional to the pressure

$$u = \frac{P}{\gamma - 1}. \quad (52)$$

For this work, we will assume the adiabatic index $\gamma = 7/5$, appropriate for molecular gas at the high temperatures (~1000 K) reached in post-shock regions. Combining equation (51) with equations (12) and (13), along with the ideal gas law $P = \rho k_B T / \mu$, we derive the gas temperature derivative

$$\frac{1}{T} \frac{dT}{dz} = \frac{1}{\rho_0 v_s^3} \frac{\gamma - 1}{w^2 - \gamma \tau} \left( \left( \frac{w^2}{\tau} - 1 \right) (\Gamma - \Lambda) + vF \right), \quad (53)$$

where $\tau = k_B T / \mu v_s^2$. We also modify equation (29) to account for the variable temperature

$$w^2 + [D(w_d - 1) - 1 - M^{-2}] w + \frac{k_B T}{\mu v_s^2} = 0. \quad (54)$$

Finally, the temperature jumps from the pre-shock value $T_{g0}$ to $T_j$ across the initial discontinuity following

$$\frac{T_j}{T_{g0}} = \left( 1 + \frac{2\gamma}{\gamma + 1} (M^2 - 1) \right) \frac{M^2 (\gamma - 1) + 2}{M^2 (\gamma + 1)}. \quad (55)$$

From Draine (1986), the rate of change of the thermal energy per unit volume of the gas due to elastic scattering by dust with a velocity-independent cross-section, $\sigma$, is

$$\Gamma_{\text{drag}} = \frac{\sigma \rho \rho_d}{m_d^2} \sqrt{\frac{8k_B T}{\pi \mu}} [k_B (T - T_d) I_2 + k_B T_d I_3], \quad (56)$$

where

$$I_2 \approx \left( 1 + \frac{9\pi}{64} \frac{\frac{1}{2}\mu |v - v_d|^2}{k_B T} \right)^{1/2} \left( 4 + \frac{8}{3} \frac{\frac{1}{2}\mu |v - v_d|^2}{k_B T} \right), \quad (57)$$

and

$$I_3 \approx \left( 1 + \frac{9\pi}{64} \frac{\frac{1}{2}\mu |v - v_d|^2}{k_B T} \right)^{1/2} \frac{8}{3} \frac{\frac{1}{2}\mu |v - v_d|^2}{k_B T}. \quad (58)$$

To calculate the dust temperature we assume that the frictional heating per grain, $\Gamma_{\text{drag}}/n_d$, is always balanced by the power radiated by a dust grain. Following Draine (2011), grains lose energy by infrared emission at a rate, per grain,

$$\Lambda_d = 4\pi r_d^2 \langle Q_{\text{abs}} \rangle \sigma_B T_d^4, \quad (59)$$

where $\sigma_B$ is the Stefan–Boltzmann constant and $\langle Q_{\text{abs}} \rangle$ is the Planck-averaged emission efficiency, which for carbon grains is

$$\langle Q_{\text{abs}} \rangle_C \sim 8 \times 10^{-7} \left( \frac{r_d}{0.1 \, \mu\text{m}} \right) \left( \frac{T_d}{K} \right)^2. \quad (60)$$

We include rotational line cooling from CO and H$_2$ using the cooling functions of Neufeld & Kaufman (1993) and Neufeld, Lepp & Melnick (1995). They give the cooling rate per volume $\Lambda(M) = n(M)n(H_2)L_M$ for a molecule M using a cooling rate coefficient $L_M$ obtained by fitting to four parameters of the form

$$\frac{1}{L_M} = \frac{1}{L_0} + \frac{n(H_2)}{L_{\text{LTE}}} + \frac{1}{L_0} \left( \frac{n(H_2)}{n_{1/2}} \right)^\alpha \left( 1 - \frac{n_{1/2}L_0}{L_{\text{LTE}}} \right). \quad (61)$$

The parameters $L_0$, $L_{\text{LTE}}$, $n_{1/2}$, and $\alpha$ are tabulated for temperatures up to a few thousand K, and depend on an optical depth parameter $\tilde{N}$. Modelling the shock as a plane-parallel slab of thickness $d$, this parameter is given as

$$\tilde{N}(\text{CO}) = \frac{n(\text{CO})d}{9\Delta v}. \quad (62)$$

We have chosen a CO abundance $x(\text{CO}) = 1.24 \times 10^{-4}$ with respect to the total hydrogen density and molecular hydrogen abundance $x(H_2) = 0.5$, both constant throughout the shock. We thus use an optical depth parameter $\tilde{N}(\text{CO}) \sim 10^{15} \text{ cm}^{-2}/(\text{km s}^{-1})$, appropriate for a shock thickness of 1 au and $\Delta v = 4 \text{ km s}^{-1}$.

Finally, we choose a pre-shock gas and dust temperature $T_{g0} = T_{d0} = 10$ K, corresponding to an isothermal sound speed $c_s = 0.188 \text{ km s}^{-1}$, and we do not allow either temperature to fall below their initial value.

Here we estimate the thickness of the dust–gas drift region for the $v_s \sim 4 \text{ km s}^{-1}$ shock we are considering. The gas temperature immediately behind the shock will jump from 10 to $\sim 10^3$ K (equation 55). The velocity jumps as

$$\frac{v_j}{v_s} = \frac{\gamma - 1}{\gamma + 1} + \frac{2}{M^2 (\gamma + 1)}, \quad (63)$$

which gives $v_j \sim v_s/6$. If we estimate the drift velocity $\Delta v = v_d - v \sim v_s/2$, assume spherical carbon dust particles with average mass density $\sim 2.2 \text{ g cm}^{-3}$, then we the size over which the dust will drift with respect to the gas is approximately

$$\Delta z \sim \frac{v_d v}{\alpha v_s I} \left( \frac{k_B T}{\mu} \right)^{-1/2} \sim 0.05 \text{ au} \left( \frac{r_d}{\mu\text{m}} \right). \quad (64)$$

Hence the shock remains thin compared to the size of the system shown in Fig. 4.

## 4.3 Results and discussion

With the heating and cooling processes in place, we numerically integrate the coupled ODEs as discussed in Section 3.1 for J-type shocks. We first investigate the effect of different sizes of dust grains
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Accretion shock profiles for $v_s = 4$ km s$^{-1}$, $\rho_0 = 4 \times 10^{-17}$ g cm$^{-3}$, pre-shock dust-to-gas ratio $D = 0.01$ and three different dust sizes, $r_d = 0.1$ (dotted), 1 (dash-dotted), and 10 (dashed) µm. Note also that the z-axis increases from the shock front at $z = 0$ towards the mid-plane of the disc. Upper left panel shows the velocity profiles, lower left panel shows the density profiles, upper right shows the temperature profile and lower panel shows the cooling profiles. The gas velocity, density, and temperature profiles (solid red lines) and the H$_2$ and CO cooling rates cannot be distinguished in the three cases. Note that in the density profile, the gas and dust densities are normalized by their respective pre-shock densities, which differs by a factor of 100 in this case.</figcaption></figure>

by considering a constant initial dust-to-gas ratio $D = 0.01$. The resulting velocity, density, temperature, and cooling rate profiles computed for dust sizes $r_d = 0.1, 1,$ and $10$ µm are shown in Fig. 6.

When $r_d = 0.1$ µm, the region of the shock with any drift between the dust and gas velocities is negligible compared to the size of the shock, and so closely resembles a one-fluid shock. However, when $r_d = 10$ µm this region is about half the size of the shock and agrees with the estimate of equation (64). Hence we expect two-fluid effects to be more prominent in accretion shocks where the dust has coagulated into large grains ($r_d > 1$ µm). In this initial region of dust-gas drift, the dust to gas ratio decreases from its initial value to a minimum 6 times smaller than D immediately after the jump. The ratio returns to its pre-shock value after ~1 au in the $r_d = 10$ µm case. Note that the shock thickness ~6 au, which is approaching the limits of validity for this model.

We next consider the dependence on the initial dust-to-gas ratio. The canonical value in the interstellar medium is $D = 0.01$, however values as high as unity have been used in protostellar discs (e.g. Dipierro et al. 2015) to account for dust migration to the inner parts of the disc. The velocity, density, temperature and cooling profiles for accretion shocks with $D = 0.01$ and $D = 1$ are shown in Fig. 7 for the case where the dust size $r_d = 10$ µm. The profiles are similar, with the larger dust-to-gas ratio resulting in a more compressed structure. There is strong CO and H$_2$ emission in these shocks as well as those shown in Fig. 6, due to the high temperatures reached (~$10^3$ K). This suggests that observing rotational lines of CO or H$_2$ are good tracers of the presence of an accretion shock, but are not sensitive to the shock parameters. On the other hand, while CO is the dominant coolant when $D = 0.01$, dust cooling dominates after

~1 au when $D = 1$. This means dust emission probes the initial dust-to-gas ratio of the accretion shock. It must be noted, however, that we required a large dust-to-gas ratio and large dust radius in order to significantly deviate from an ordinary hydrodynamic shock. In addition, the infalling material will comprise dust grains with a distribution of sizes, rather than the single size assumed here. The relative drift between grains of different radii induces sticking and shattering in grain-grain collisions, and these effects must be carefully treated (e.g. Guillet, Jones & Pineau des Forêts 2009; Guillet, Pineau des Forêts & Jones 2011).

We have presented an astrophysical application of two-fluid dust-gas shocks by studying the accretion shock above a protoplanetary disc. We have simplified the system by not including any chemical reactions or evaporation of grain mantles. At the temperatures reached (~$10^3$ K) there is significant driving of neutral-neutral reactions, and coolants such as H$_2$O and OH could be produced. Cooling by these molecules could change the detailed structure of the shock and/or provide radiative signatures of the shock parameters. Our simple treatment has shown that a detailed analysis of dust-gas shocks is needed to investigate infalling material on to protoplanetary discs. In addition, we have only considered a shock velocity appropriate for a solar mass protostar. Lower mass protostars would lead to lower shock velocities (equation 39) and therefore lower post-shock temperatures (equation 55). Thus, detailed shock models at a range of velocities could provide radiative signatures of protostellar mass. In this work, we adopted a shock velocity and pre-shock density appropriate for the minimum-mass solar nebula. This could be improved by considering more appropriate physical conditions for earlier stages of star formation and accounting for material flowing in from all directions.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Accretion shock profiles for $v_s = 4$ km s$^{-1}$, $\rho_0 = 4 \times 10^{-17}$ g cm$^{-3}$, and dust size $r_d = 10$ µm. Note also that the z-axis increases from the shock front at $z=0$ towards the mid-plane of the disc. Upper left panel shows the velocity profiles, lower left panel shows the density profiles, upper right shows the temperature profile, and lower panel shows the cooling profiles. Red lines refer to shocks with initial dust-to-gas ratio $D = 0.01$ while black lines have $D = 1$. Dashed lines refer to dust variables and solid lines refer to gas variables in the velocity and density profiles. In the cooling profile, solid lines refer to CO, dotted lines refer to H$_2$ and dashed lines refer to dust. Note that in the density profile, the gas and dust densities are normalized by their respective pre-shock densities, which differs in the $D = 0.01$ case.</figcaption></figure>

## 5 CONCLUSION

We have numerically solved the two-fluid dust-gas equations assuming a steady-state, planar structure. Two distinct shock solutions exist where the gas fluid drags the dust fluid along through a discontinuity (J-type) or smoothly (C-type) until both fluids settle on to post-shock values. These shocks are ideal tests for benchmarking the behaviour of numerical codes seeking to simulate dusty gas with different expressions for the drag or dust-to-gas mass density ratios. Our **PYTHON** code that returns shock solutions for user defined parameters is publicly available on the Python Package Index<fnref n="4" /> and BitBucket.<fnref n="5" />

We used a J-type two-fluid dust-gas shock to study the accretion shock settling material on to a protoplanetary disc. We found that two-fluid effects are most likely to be important for larger grains ($r_d > 1$ µm). Dust emission from within the shock front was found to be a sensitive probe of the dust-to-gas ratio that eventually falls on to the protoplanetary disc. This work shows that a detailed analysis of two-fluid dust-gas shocks could be a fruitful avenue to investigating the composition of infalling material on to protoplanetary systems.

## ACKNOWLEDGEMENTS

The authors gratefully acknowledge discussions with James Tocknell, Shane Vickers, and Birendra Pandey, and especially thank Daniel Price for bringing our attention to this problem and for

comments on the manuscript. This research was supported by the Australian Research Council through Discovery Project grant DP130104873. AL was supported by an Australian Postgraduate Award. This research made use of Astropy, a community-developed core **PYTHON** package for Astronomy (Astropy Collaboration 2013).

## REFERENCES

* ALMA Partnership, 2015, ApJ, 808, L3

* Astropy Collaboration, 2013, A&A, 558, A33

* Bai X.-N., Stone J. M., 2010a, ApJS, 190, 297

* Bai X.-N., Stone J. M., 2010b, ApJ, 722, L220

* Carrier G. F., 1958, J. Fluid Mech., 4, 376

* Dipierro G., Price D., Laibe G., Hirsh K., Cerioli A., Lodato G., 2015, MNRAS, 453, L73

* Draine B. T., 1986, MNRAS, 220, 133

* Draine B. T., 2011, Physics of the Interstellar and Intergalactic Medium. Princeton Univ. Press, Princeton, NJ

* Epstein P. S., 1924, Phys. Rev., 23, 710

* Glover S. C. O., Clark P. C., 2012, MNRAS, 421, 9

* Guillet V., Jones A. P., Pineau des Forêts G., 2009, A&A, 497, 145

* Guillet V., Pineau des Forêts G., Jones A. P., 2011, A&A, 527, 123

* Hayashi C., 1981, Prog. Theor. Phys. Suppl., 70, 35

* Hopkins P. F., Lee H., 2016, MNRAS, 456, 4174

* Igra O., Ben-Dor G., 1988, Appl. Mech. Rev., 41, 379

* Johansen A., Henning T., Klahr H., 2006, ApJ, 643, 1219

* Kriebel A. R., 1964, J. Basic Eng., 86, 655

* Kwok S., 1975, ApJ, 198, 583

* Laibe G., Price D. J., 2011, MNRAS, 418, 1491

* Laibe G., Price D. J., 2012a, MNRAS, 420, 2345

* Laibe G., Price D. J., 2012b, MNRAS, 420, 2365

* Larson R. B., 2003, Rep. Prog. Phys., 66, 1651

<footnote n="4">⁴ https://pypi.python.org/pypi/DustyShock</footnote>

<footnote n="5">⁵ https://bitbucket.org/AndrewLehmann/dustyshock</footnote>
Launhardt R. et al., 2013, A&A, 551, A98

Lorén-Aguilar P., Bate M. R., 2014, MNRAS, 443, 927

Miura H., 1972, JPSJ, 33, 1688

Neufeld D. A., Kaufman M. J., 1993, ApJ, 418, 263

Neufeld D. A., Lepp S., Melnick G. J., 1995, ApJS, 100, 132

Paardekooper S.-J., Mellema G., 2006, A&A, 453, 1129

Planck Collaboration XXXIII, 2016, A&A, 586, A136

Saito T., Marumoto M., Takayama K., 2003, Shock Waves, 13, 299

Sod G. A., 1978, J. Comput. Phys., 27, 1

Stephens I. W., Evans J. M., Xue R., Chu Y.-H., Gruendl R. A., Segura-Cox
D. M., 2014, ApJ, 784, 147

Toth G., 1994, ApJ, 425, 171

Wardle M., 2007, Ap&SS, 311, 35

Weidenschilling S. J., 1977, Ap&SS, 51, 153

This paper has been typeset from a TEX/LATEX file prepared by the author.