# Dynamics of cold circumstellar gas in debris discs

Can Cui, Sebastian Marino, Quentin Kral, Henrik Latter

► To cite this version:

Can Cui, Sebastian Marino, Quentin Kral, Henrik Latter. Dynamics of cold circumstellar gas in debris discs. Monthly Notices of the Royal Astronomical Society, 2024, 530, pp.1766-1780. 10.1093/mn-ras/stae924. insu-04853447

HAL Id: insu-04853447

https://insu.hal.science/insu-04853447v1

Submitted on 23 Dec 2024

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Dynamics of cold circumstellar gas in debris discs

Can Cui<sup>1,★</sup> Sebastian Marino<sup>1,2</sup> Quentin Kral<sup>1,3</sup> and Henrik Latter<sup>1</sup>

Accepted 2024 March 25. Received 2024 March 13; in original form 2023 October 23

## ABSTRACT

Mounting observational evidence indicates that cold circumstellar gas is present in debris disc systems. This work focuses on various dynamical processes that debris disc gas may undergo. We review five mechanisms that can transport angular momentum and their applications to debris discs. These include molecular viscosity, hydrodynamic turbulence, magnetohydrodynamic (MHD) turbulence, magnetized disc winds, and laminar magnetic stress. We find that molecular viscosity can result in $\alpha$ as high as $\lesssim 0.1$ for sufficiently low densities, while the Rossby wave instability is a possible source of hydrodynamic turbulence and structure formation. We argue that the vertical shear instability is unlikely due to the long cooling times. The onset of the magnetorotational instability (MRI) is dichotomous: for low-density discs the MRI can be excited at the mid-plane, while for high-mass discs it may only be operating at $z > 2 - 3H$, if at all. The MHD wind and laminar magnetic stress mechanisms rely on the configuration and strength of any background large-scale magnetic field, the existence of which is uncertain and possibly unlikely. We conclude that the dominant mechanism and its efficiency in transporting angular momentum vary from one system to the other, depending especially closely on the gas density. More detailed analyses shall be performed in the future focusing on representative, nearby debris discs.

**Key words:** hydrodynamics – MHD – turbulence – methods: analytical – planetary systems.

## 1 INTRODUCTION

Main-sequence stars are commonly orbited by solid circumstellar material, such as planets and debris (Wyatt 2008; Matthews et al. 2014; Hughes, Duchêne & Matthews 2018). The latter is made up of planetesimals as well as the dust and gas derived from them. Debris discs are tenuous, optically thin (with respect to dust), and may persist over gigayears. The nearest debris disc is the Solar system’s asteroid and Kuiper belts. Extrasolar systems also commonly harbour debris discs: to name a few, Fomalhaut, $\beta$ Pictoris, Vega, and $\epsilon$ Eridani (e.g. Holland et al. 1998).

Originally, the absence of gas was a distinguishing characteristic of debris discs. However, mounting observational evidence in recent years indicates that cold circumstellar gas is present around a large fraction of stars hosting debris discs, identified mainly via CO gas emission at millimetre wavelengths (see Strøm et al. 2020; Marino 2022, for recent reviews). Large quantities of CO gas were found mainly around young A-type stars initially (e.g. Zuckerman & Song 2012; Kóspál et al. 2013; Dent et al. 2014; Lieman-Sifry et al. 2016; Moór et al. 2017), while later and deeper observations detected CO gas at lower levels around later spectral types (Marino et al. 2016; Matrà et al. 2019a; Kral et al. 2020b). CO is known to be short-lived due to photodissociation driven by the ultraviolet (UV) interstellar radiation field, unless shielded, and its dissociation products have also been found in some of these systems: C I, C II, and O I (e.g.

Roberge et al. 2006; Cataldi et al. 2014, 2018; Kral et al. 2019). Atacama Large Millimeter Array - ALMA Radio Telescope (ALMA) observations of these species have shown that the gas and dust spatial distributions overlay, with large inner cavities of tens of au in radius (e.g. Kóspál et al. 2013; Dent et al. 2014; Marino et al. 2016; Moór et al. 2017). Some of the observations also suggest that the gas is radially spreading relative to the dust (e.g. Kóspál et al. 2013; Moór et al. 2013).

The origin of the circumstellar gas still remains unclear, particularly in very massive gas discs. Both primordial and secondary origin scenarios have been proposed to explain some of the observed features. The primordial origin idea proposes that the circumstellar gas is inherited from earlier Class II protoplanetary discs, and hence is expected to be $H_2$-dominated (e.g. Moór et al. 2017; Nakatani et al. 2020; Klusmeyer et al. 2021; Smirnov-Pinchukov et al. 2022). Meanwhile, the secondary origin idea suggests that the gas is released by solid bodies in the planetesimal belt, via the grinding down of larger volatile-rich solid bodies in a collisional cascade, sublimation by external heating, or radionuclides; the gas is hence $H_2$-poor (e.g. Czechowski & Mann 2007; Moór et al. 2011; Zuckerman & Song 2012; Dent et al. 2014; Matrà et al. 2017a).

Though tenuous, circumstellar gas has been invoked in a number of physical processes in debris disc systems. Gas can influence dust dynamics and shape its spatial distribution (e.g. Takeuchi & Artymowicz 2001; Krivov et al. 2009; Olofsson et al. 2022). The photoelectric instability, whose onset depends on gas radial pressure gradients, can concentrate dust into annular structures when the disc is optically thin to starlight (Lyra & Kuchner 2013). Circumstellar gas

<footnote n="1"><sup>1</sup>DAMTP, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK</footnote>

<footnote n="2"><sup>2</sup>Department of Physics and Astronomy, University of Exeter, Stocker Road, Exeter EX4 4QL, UK</footnote>

<footnote n="3"><sup>3</sup>LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université Paris Cité, 5 place Jules Janssen, F-92195 Meudon, France</footnote>

<footnote>* E-mail: can.cui@astro.utoronto.ca</footnote>
can also accrete on to neighbouring planets and hence contribute to the chemical composition of secondary planetary atmospheres (Kral, Davoult & Charnay 2020a). Finally, it has been proposed that gas dynamics might be the origin of the dust (and gas) clump detected in β Pic, via a gas vortex triggered by the Rossby wave instability (RWI; Skaf et al. 2023).

Despite its importance in many physical processes, the understanding of gas dynamics in debris discs is relatively underdeveloped (but see Kral & Latter 2016, for low gas density disc of β Pic). Perhaps, the most salient dynamical questions in this context concern the radial spreading and vertical mixing of gas. The former is closely related to angular momentum transport, and the latter strongly influences the degree of mid-plane shielding of CO by surface C1, which in turn helps determine the origin of the gas. To this aim, this paper seeks solutions to some fundamental questions regarding gas dynamics. These include the following: what are the nature and strength of the angular momentum transport? Are hydrodynamic mechanisms important, such as molecular viscosity and hydrodynamic turbulence? If the gas dynamics is controlled by magnetohydrodynamics (MHD), is the disc magnetorotational instability (MRI) turbulent or do magnetized disc winds dominate mass accretion? Do non-ideal MHD effects come into play? In short, what is the origin of any 'anomalous viscosity' in debris discs, and can a simple Shakura-Sunyaev α parameter describe the underlying physical processes adequately?

The paper is organized as follows. In Section 2, we start with a background introduction on gas observations, theoretical modelling, and a comparison of gas physical properties between protoplanetary and debris discs. In Section 3, our methods are outlined. In Section 4, we introduce hydrodynamic mechanisms that may shape gas structures. In Section 5, we investigate the MHD effects on the gas dynamics. Finally, we discuss the main results in Section 6 and conclude in Section 7.

# 2 BACKGROUND

## 2.1 Observations of cold gas

The main observed species at tens of au have been CO, C1, CII, and OI. CO has a longer photodissociation lifetime compared to other molecules, and its bright rotational lines at millimetre wavelengths are where ALMA is most sensitive. C1 most likely originates from the photodissociation of CO and is long-lived, allowing resolved observations with ALMA. Additionally, CII and OI have been detected by the Herschel far-infrared (far-IR) survey in a few systems. Of these observable species, the gas mass tends to be dominated by CO, except that β Pic and 49 Ceti may have slightly higher C1 masses (Cataldi et al. 2023).

There is a wide distribution of CO masses from ~0.1 M<sub>⊕</sub> down to marginally detected levels of 10<sup>-7</sup> M<sub>⊕</sub> (e.g. Matrà et al. 2017b). The high end of the distribution, M<sub>CO</sub> ~ 10<sup>-1</sup> to 10<sup>-4</sup> M<sub>⊕</sub>, is dominated by A-type stars with bright debris discs (Marino et al. 2020), while the lower end of the distribution contains systems with a wide range of spectral types, from B- to M-type stars (e.g. Marino et al. 2016; Matrà et al. 2019a; Kral et al. 2020a; Rebollido et al. 2022). Age and dust levels also play a role. Massive gas discs with M<sub>CO</sub> > 10<sup>-4</sup> M<sub>⊕</sub> are found mostly around A-type stars younger than 50 Myr, and with dust fractional luminosities above 10<sup>-4</sup> (Moór et al. 2017). However, this might be due to observational biases since no surveys have been dedicated to searching for CO around older and less dusty discs. Among the surveyed dust-rich debris discs around young A-type stars, CO is readily detected in 11/17, 10 of which have CO

gas masses above 10<sup>-4</sup> M<sub>⊕</sub> (Moór et al. 2017; Kennedy et al. 2018; Hales et al. 2022).

Despite the abundance of gas, its bulk composition is still largely unknown. In particular, there are no strong constraints on the presence of H or H<sub>2</sub> in the most massive systems, which would dominate the mass if the gas was of primordial origin. For lower mass systems, it is thought that the H<sub>2</sub> content is low because of the subthermal gas excitation (Matrà et al. 2017a). There have been multiple attempts to detect carbon- and oxygen-bearing molecules other than CO that are frequently detected in protoplanetary discs, such as CN, HCN, HCO<sup>+</sup>, CCH, CH<sub>3</sub>OH, CS, and SiO (Klusmeyer et al. 2021; Smirnov-Pinchukov et al. 2022). However, none have been found, ruling out abundances relative to CO similar to protoplanetary discs. These low abundances, however, can be simply explained by the optically thin nature of debris discs: the UV stellar and interstellar radiation can penetrate the disc and quickly dissociate most molecules, even if the gas is interstellar medium-like with a high abundance of H or H<sub>2</sub>.

## 2.2 1D modelling of secondary gas evolution

Previously, modelling of the gas evolution has focused on two issues: how gas viscously expands radially, forming an accretion disc, and how CO photodissociates, producing atomic carbon and oxygen.

To model the radial viscous expansion, two free parameters, the Shakura-Sunyaev α and the CO gas production rate, have been generically employed and varied, to reproduce CO gas surface densities and radial distributions. Kral et al. (2016, 2017) presented a one-dimensional (1D) radial viscous evolution model to predict the CO, C1, CII, and OI gas distribution and abundances. This model simulated how gas released at a fixed rate at the planetesimal belt location will viscously spread, using an α model, ionize, and heat or cool. Subsequently, it was realized that CO could be shielded by C1, and thus if enough C1 accumulates, the shielding could protect CO and explain the rich CO detected around many A-type stars (Kral et al. 2019). This realization meant that all gas found in debris discs could be of secondary origin. Marino et al. (2020) incorporated this into a 1D radial evolution model to produce synthetic populations of discs where the gas production decreases with age as discs lose mass. Comparing these against observations found that the distribution of CO masses is best explained by α = 0.1.

Besides the viscous expansion in radius, the gas can be blown out over short time-scales and persist beyond the planetesimal belt by stellar wind protons, creating belt winds of CO, CO<sup>+</sup>, C1, and OI (Kral et al. 2021, 2023). This occurs when the gas density is below a critical threshold, for example, in the least massive gas discs detected so far.

More recently, attention has been drawn to the vertical distribution of the gas (Olofsson et al. 2022). The gas' precise vertical structure significantly influences the lifetime of CO, and thus on how much CO gas a disc can accumulate. The vertical mixing is determined by the vertical turbulent diffusion. If the gas is well mixed, Cataldi et al. (2020) showed that C1 shielding is weak and the CO lifetime increases only linearly with the column density of C1. Instead, if the gas is not mixed, where C1 is distributed on a layer high above the mid-plane and above the CO gas, it can shield CO more effectively, thus increasing the lifetime of CO exponentially (Kral et al. 2019).

Marino et al. (2022) investigated the vertical distribution of gas in the secondary origin scenario, using a 1D model that resolves the gas disc vertically. They showed that both scenarios are possible depending on how strong vertical mixing is, parametrized by a
<table><caption>**Table 1.** Comparison of physical parameters between protoplanetary and debris discs at 100 au (in cgs units).</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Protoplanetary</th><th>Debris</th></tr></thead><tbody><tr><td>Gas number density</td><td>n</td><td>10<sup>5</sup>–10<sup>9</sup></td><td>10<sup>2</sup>–10<sup>6</sup></td></tr><tr><td>Mean free path</td><td>l/H</td><td>10<sup>-8</sup> to 10<sup>-4</sup></td><td>10<sup>-5</sup> to 1</td></tr><tr><td>Ionization</td><td>x<sub>e</sub></td><td>10<sup>-13</sup> to 10<sup>-10</sup></td><td>10<sup>-7</sup> to 1</td></tr><tr><td>Dust-to-gas ratio</td><td>ρ<sub>d</sub>/ρ<sub>g</sub></td><td>10<sup>-2</sup></td><td>1–10<sup>5</sup></td></tr><tr><td>Stokes number</td><td>St</td><td>&lt;1</td><td>&gt;1</td></tr></tbody></table>

vertical (dimensionless) diffusivity $\alpha_v$. If the vertical mixing is efficient, CO and C I would have similar vertical distributions and the shielding by C I would weaken. If the vertical diffusivity is much weaker than the radial, and gas is released at a high rate, then a layered structure can be formed. This is a consequence of the disc becoming optically thick to photodissociating UV photons and the vertical diffusion time-scale being shorter than the time-scale at which gas is removed or replenished. In this case, the CO dominates the mid-plane, the C I density peaks above the mid-plane, at a height that depends on the gas surface density, and a low-density layer of C II dominates the uppermost layers.

## 2.3 Protoplanetary versus debris discs

Most previous work, as discussed in the last section, explicitly includes the effects of dynamical processes via simple parametrizations (such as by dimensionless diffusivities, $\alpha$ and $\alpha_v$), but without specifying the nature of those mechanisms or interrogating the accuracy of the parametrizations. It is the goal of this paper to explore the poorly understood dynamics of debris disc gas in more detail. To make inroads here, it is almost unavoidable to first consider protoplanetary discs and their rich dynamical behaviour. To help facilitate this, it is beneficial to compare the key characteristic physical and dynamical parameters of the two systems. These parameters and respective values are listed in Table 1. To make easy comparison, the values are calculated at 100 au and vary from the mid-plane to the surface.

Protoplanetary discs are gas-rich and mainly composed of H$_2$. The gas surface density is about $10^2$ g cm$^{-2}$ and number density of $10^5$–$10^9$ cm$^{-3}$ at 100 au (Weidenschilling 1977; Henning & Semenov 2013). Debris discs are much more dilute. The surface density of CO typically peaks at 100 au and varies from $10^{-7}$ to $10^{-6}$ g cm$^{-2}$ for tenuous discs, such as $\beta$ Pic, to $10^{-6}$ to $10^{-4}$ g cm$^{-2}$ for dense discs, such as 49 Ceti (Kral & Latter 2016; Hughes et al. 2017; Higuchi et al. 2020; Cataldi et al. 2023). The gas number density is about $10^2$–$10^7$ cm$^{-3}$ (Section 5.1.1). Note that the high end of the range of CO surface densities could be much higher. This is because $^{12}$CO emission is optically thick and CO isotopologue-selective photodissociation is usually ignored when estimating CO gas masses. The low end, on the other hand, is determined by the current sensitivity limits.

The gas mean free path determines the collisionality and molecular viscosity. Protoplanetary discs are generically neutral and collisional. Taking the cross-section of H$_2$ to be approximately $10^{-15}$ cm$^2$, the mean free path is $l/H \sim 10^{-8}$ to $10^{-4}$, where $l$ is the gas mean free path and $H$ the pressure scale height. Debris discs possess a non-negligible fraction of ionized gas (C II). The mean free path for C II–C II collisions is short due to the Coulomb encounters, and $l/H \sim 10^{-8}$ to $10^{-4}$. On the other hand, neutral–neutral and neutral–ion collisions have longer mean free paths, resulting in $l/H \sim 10^{-5}$ to 1 (Kral et al. 2016; Section 4.1). The ionization fraction is a key parameter for MHD processes, defined as the ratio of electron number

density to neutral number density, $x_e = n_e/n$. The ionization fraction is extremely low for protoplanetary discs, $x_e \sim 10^{-13}$ to $10^{-10}$ (Lesur 2021), whereas it is substantially higher in debris discs, $x_e \sim 10^{-7}$ to 1 (Fig. 4; Kral et al. 2016; Marino 2022). Rough observational constraints have been placed on the ionization fraction in some debris disc systems and found to be $x_e \sim 0.1–1$ based on C I and C II (Cataldi et al. 2023).

Aerodynamic coupling between dust and gas is subject to a mutual drag force. The force is determined by the dust-to-gas density ratio and the dimensionless Stokes number. Protoplanetary discs are gas-dominated; the dust-to-gas ratio reaches only 1 per cent. On the other hand, debris discs are dust-dominated; the dust-to-gas (CO) ratio can vary from unity to $10^5$ (Rebollido et al. 2022), though this quantity is uncertain because it depends on modelling assumptions, such as the size of the largest solids. The Stokes number, $St$, is defined as the ratio of the stopping time to the dynamical time. The stopping time is the time-scale required for the dust to respond to differences in the gas’ and dust’s velocities; the dynamical time is the local orbital period of the gas. When $St$ is smaller (greater) than unity, the dust and gas are dynamically well (poorly) coupled. For micron- to millimetre-sized dust at the mid-plane and at 100 au in protoplanetary discs, $St < 1$. For the same sized dust in debris discs, generally $St > 1$, because of their considerably lower gas densities (Marino et al. 2020; Huang & Ormel 2023).

## 2.4 Angular momentum transport theory

To better understand angular momentum transport in debris discs, and thus how they spread and accrete, we briefly introduce its underlying mathematical theory. The momentum equation in conservative form is (Landau & Lifshitz 1959)

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

where the momentum flux tensor is defined as

$$ \mathbf{M} \equiv \rho \mathbf{v} \mathbf{v} - \frac{\mathbf{B} \mathbf{B}}{4\pi} + \left( \mathbf{P} + \frac{\mathbf{B}^2}{8\pi} \right) \mathbf{I} - \mathbf{\tau}, \qquad (2) $$

in which $\mathbf{v}$ and $\mathbf{B}$ are the gas velocity and magnetic field vector, $\mathbf{I}$ is the identity tensor, and $\mathbf{\tau}$ is the viscous stress tensor,

$$ \mathbf{\tau} = \eta[\nabla\mathbf{v} + (\nabla\mathbf{v})^T] - \frac{2}{3}\eta(\nabla \cdot \mathbf{v})\mathbf{I}. \qquad (3) $$

The dynamical viscosity is denoted by $\eta = \rho v$, and $v$ is the kinematic viscosity. In cylindrical coordinates, the viscous stress tensor is $\mathbf{\tau}_{R\phi} = \eta R \partial\Omega / \partial R$ for a Keplerian disc with angular speed $\Omega$. Note that we have omitted the dust drag on the gas in this model, and throughout the paper. This can be justified for a large class of debris discs (see Kral & Latter 2016; Kral et al. 2016), but when the dust-to-gas ratio is large the drag term must be reinstated.

Under axisymmetry, the $R\phi$ and $z\phi$ components of the total stress tensor $\mathbf{T}$ are

$$ T_{R\phi} = \rho \delta v_R \delta v_{\phi} - \frac{B_R B_{\phi}}{4\pi} + \frac{3}{2} \eta \Omega, \qquad (4) $$

$$ T_{z\phi} = \rho v_z v_{\phi} - \frac{B_z B_{\phi}}{4\pi}, \qquad (5) $$

where $\delta v_R$ and $\delta v_{\phi}$ are the radial and azimuthal velocity fluctuations, respectively.

From above, the angular momentum can be transported radially and vertically (see reviews by Armitage 2011; Weiss, Bai & Fu 2021; Lesur et al. 2023). The $T_{R\phi}$ term facilitates the radial transport of angular momentum, where part of the gas loses angular momentum
and accretes on to the central star, and the rest gains angular momentum and expands radially outwards. The α-disc model introduces a dimensionless parameter to quantify the strength of the Rφ stress,

$$ \alpha = \frac{T_R \phi}{P}, \qquad (6) $$

where P is the gas pressure and an overbar denotes a vertical average. The $T_{z\phi}$ term represents the vertical transport of angular momentum. Magnetized disc winds can drive such vertical transport. They extract angular momentum from the disc by exerting a torque at the disc surface, forcing the entire disc to accrete.

Five mechanisms can be identified that contribute to the angular momentum transport process: molecular viscosity, hydrodynamic turbulence, magnetic turbulence, magnetized disc winds, and laminar magnetic stress. The molecular viscosity can transfer angular momentum through the last term of the $T_{R\phi}$ stress. Turbulence contributes to transport via the $T_{R\phi}$ stress. Hydrodynamic turbulence only manifests in the first term of $T_{R\phi}$. One example might be that arises from the vertical shear instability (VSI; Section 4.2). Magnetic turbulence can produce both the first and second terms in $T_{R\phi}$. One well-known example is the MRI (Section 5.3). On the other hand, a large-scale laminar magnetic stress can contribute to the angular momentum transport through the second term in $T_{R\phi}$. The second term of $T_{z\phi}$ arises from large-scale poloidal magnetic fields, which may permit the launching of magnetized disc winds. In the following sections, we will assess the viability of each of these mechanisms in debris disc gas.

## 3 METHODS

To help us investigate the different hydrodynamic and MHD processes going on in debris disc gas, we require a disc model describing the spatial distribution of the gas species. Our reference discs in this paper have been calculated by Marino et al. (2022), the details of which we now describe.

The model of Marino et al. (2022) follows the vertical evolution of gas using a 1D grid that resolves the vertical axis. It considers a belt centred at 100 au orbiting a $1.5 M_{\odot}$ star. The model accounts for the release of CO gas from a planetesimal belt at different release rates, and it then follows how photodissociation and ionization by interstellar UV and radial viscous spreading and vertical mixing shape the vertical distribution of gas species. The viscous evolution and vertical turbulent mixing are parametrized by the Shakura–Sunyaev $\alpha$ parameter. The model adopted in our paper employs radial and vertical diffusion coefficients of $\alpha = 10^{-2}$ and $\alpha_v = 10^{-2}$, respectively. The disc aspect ratio is set to $H/R = 0.05$, equivalent to a temperature of ~50 K. The mean molecular weight $\mu$ in their calculation is assumed to be equal to 14, equivalent to a gas dominated by carbon and oxygen in equal proportions.

A key parameter in the model is the solids' gas release rate, which varies between $10^{-3}$ and $10^1 M_{\oplus} Myr^{-1}$. This is a wide range that covers the expected release rate in known systems with detected CO gas, from the bottom of the distribution (~$10^{-2} M_{\oplus} Myr^{-1}$, e.g. HD 181327; Marino et al. 2016) to the most massive gas-rich systems (~1 $M_{\oplus} Myr^{-1}$, e.g. HD 121617; Kral et al. 2019). Rates as high as $10 M_{\oplus} Myr^{-1}$ might be plausible for a few Myr depending on the mass of the belt and the processes that drive the release of CO from solids (Bonsor et al. 2023). These rates map straightforwardly on to the steady state total gas surface densities at 100 au, namely $\Sigma \sim 1.7 \times 10^{-6}$ to $1.7 \times 10^{-2}$ g cm$^{-2}$, once $\alpha$ is specified. Note that the high end of this range is a hundred times larger than the highest CO surface densities observed using ALMA (~2 $\times 10^{-4}$ g cm$^{-2}$;

Cataldi et al. 2023). However, those observational estimates only account for CO and ignore the possibility of CO isotopologue-selective photodissociation. This means that the gas surface densities could be orders of magnitude higher.

We note that Marino et al. (2022) treat the $\alpha$ and $\alpha_v$ as free parameters that are independent of mass production rate, density, temperature, radius, height, etc. However, as we will show in Sections 4.1, 4.2, and 5.3, molecular viscosity, the VSI, and the MRI all result in different values of $\alpha$ over density and space. This breaking of the connection between the dynamical parameters and the disc thermodynamic state may lead to inconsistencies and possibly obscure important trends in the modelling. Future global numerical simulations may be needed to self-consistently and thus reliably obtain the values of $\alpha$.

Finally, Marino et al. (2022) investigate various vertical diffusion coefficients $\alpha_v$, ranging from $10^{-5}$ to $10^{-2}$, to explore the vertical distributions of gas species as $\alpha_v/\alpha$ varies. While our chosen value for $\alpha$ and $\alpha_v$ of $10^{-2}$ may seem high and arbitrary, we note that our key results are not sensitive to this choice but rather to the surface density of gas (set by the gas release rate and $\alpha$). A weaker $\alpha_v$ relative to $\alpha$ does not quantitatively change our main conclusions.

## 4 HYDRODYNAMICS

In this section, we present hydrodynamic mechanisms of angular momentum transport that may be relevant in the context of debris discs. Molecular viscosity, which can be of consequence due to the low gas densities, is discussed first, in Section 4.1. Two hydrodynamic instabilities that have the potential to prevail in debris discs are also explored: the VSI (Section 4.2) and the RWI (Section 4.3).

### 4.1 Non-negligible molecular transport

As described in Section 2.4, angular momentum can be transported by molecular viscosity, manifesting in the last term of equation (4). The kinematic viscosity can be estimated as $\nu \sim lc_s$, where $l$ and $c_s$ are the mean free path and sound speed of the molecules, respectively. Hence, $\alpha \sim \Omega\rho v/P \sim l/H$. Transport by molecular viscosity requires some degree of collisionality, and the estimate for $\nu$ above holds only if particle trajectories between collisions are not overly influenced by body forces and thus follow straight lines. When $l < H$, these restrictions are satisfied. However, when $l > H$, gas particles will tend to undergo vertical and radial epicyclic oscillation with a characteristic length-scale $\sim H$, and a different estimate for $\nu$ is needed (see Goldreich & Tremaine 1978; Pringle 1981).

We first compute the mean free path for different gas species, including collisions among neutrals (CO, C I, and O I) and ions (C II). The collisions with electrons are ignored due to their negligible inertia compared to other gas species. The mean free path is $l = 1/(n\sigma)$, where $n$ is the number density of the impacted species, and $\sigma$ is the relevant cross-section of the interaction. A rigorous calculation of the mean free path requires the cross-section for each projectile–target pair. Here, we simplify the calculation by considering three cases: neutral–neutral collision, neutral–C II collision, and C II–C II collision. For C II–C II collisions, $\sigma = \pi b^2$ is taken as the Coulomb cross-section, where $b = e^2/(4\pi\epsilon_0 m\mu c_s^2)$, $m\mu$ the reduced mass of the two colliders, and $c_s$ the sound speed. For neutral–neutral and neutral–C II collisions, the cross-section is approximated as geometric, $\sigma \approx \pi (10^{-8} \text{ cm})^2$. As the Coulomb impact parameter $b$ is orders of magnitude larger than the atomic size of the gas species, the mean free path of C II–C II collisions is much shorter than the other two. Thus, neutral–neutral and neutral–C II collisions result in
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Neutral gas mean free path *l* divided by pressure scale height *H* as a function of height. Colours represent different gas surface densities. The horizontal dashed line denotes *l*/*H* = 1.</figcaption></figure>

more efficient angular momentum transport than C II–C II collisions, noting also that the neutral species never comprise less than half the gas mass.

In Fig. 1, we plot the neutral mean free path *l* divided by the gas pressure scale height *H* as a function of height. The pressure scale height is calculated by assuming a constant temperature *T* = 50 K with height throughout the paper. As is clear, in the bulk of the gas disc the mean free path is much less than the scale height and the assumptions made above are validated. Even the most dilute discs in our model sample are collisional towards the mid-plane, in agreement with previous work (Kral et al. 2016). Note that values of *l* ≳ *H* (above the dashed horizontal line) are inaccurate because of the dominance of body forces in this weakly collisional regime.

Next, by vertically averaging the dynamical viscosity for *l* < *H*, we obtain the molecular viscous *α* by equation (6). We find *α* ~ 3.2 × 10⁻⁵, 2 × 10⁻⁴, 1.6 × 10⁻³, 1.2 × 10⁻², and 9 × 10⁻² for the five surface densities selected in Fig. 1, ranging from high to low, respectively. Lower density results in longer mean free paths, and thus enhanced transport. However, note that, as mentioned above, in more dilute discs, which are weakly collisional at all heights, gas particles undergo confined epicyclic oscillation before colliding and thus molecular transport, and the viscosity, will be lower (section 4; Pringle 1981).

The molecular viscosity sets a non-negligible lower limit for angular momentum transport in debris disc gas. If the disc is completely laminar, exhibiting no other hydrodynamic or MHD process, then the disc will still spread and evolve under the action of its intrinsic molecular viscosity. Note that in such a situation the disc may be viscously unstable, due to the inverse proportionality between *ν* and surface density, and may possibly break up into rings and gaps (Schmidt et al. 2009). The time-scale for this process is of the order of (*α*Ω)⁻¹, and thus might be relevant in more dilute systems. Note also that this background low level of molecular viscosity could interfere with other mechanisms that can drive angular momentum transport, for example, the onset and development of the hydrodynamic instabilities described later. However, we exclude its impact in this work for simplicity.

Finally, molecular transport may also give rise to significant vertical mixing, and thus to a non-turbulent $\alpha_v$. Naive estimates from kinetic theory would suggest that $\alpha_v$ is of the same order of magnitude as the molecular $\alpha$. As a consequence, in more dilute system it could provide a non-trivial contribution to the vertical transport, in line with our adoption of the fiducial estimate $\alpha_v = 10^{-2}$.

## 4.2 Vertical shear instability

The VSI is a promising hydrodynamic instability that may be operating in debris disc gas (Nelson, Gressel & Umurhan 2013; Latter & Papaloizou 2018; Cui & Lin 2021; Cui & Latter 2022; Latter & Kunz 2022; Svanberg, Cui & Latter 2022; Dang, Cui & Barraza-Alfaro 2024). Two conditions must be met for the onset of the VSI: a vertical gradient in the angular velocity, i.e. $d\Omega/dz \neq 0$, and a sufficiently short thermal relaxation time-scale $t_c$. The former is readily satisfied in a thin locally isothermal disc if there is a radial temperature gradient (Nelson et al. 2013). Numerical simulations of the VSI conducted, in the protoplanetary disc context, yield a Shakura–Sunyaev $\alpha$ ranging from ~10⁻³ with locally isothermal equation of state (Nelson et al. 2013) down to ~10⁻⁵ to 10⁻⁴ when including realistic radiative transport (Stoll & Kley 2016; Flock et al. 2017). Non-ideal MHD simulations of VSI, applicable to the outer regions of protoplanetary discs, suggest $\alpha \sim 10^{-4}$ to 10⁻³ (Cui & Bai 2020, 2022).

The critical thermal relaxation time-scale for the onset of the VSI was derived by Lin & Youdin (2015) in a vertically global and radially local disc model. They show that the instability criterion can be framed as

$$t_c < \frac{|q|h}{(\gamma - 1)\Omega_K}, \qquad (7)$$

where *h* = *H*/*R* is the disc aspect ratio, *q* is the radial power-law index of the temperature, and *γ* is the adiabatic index. In the context of debris discs, using *h* ≈ 0.05 for *T* = 50 K at 100 au, *q* ≈ −0.5, and *γ* ≈ 1.4, we find $t_c \lesssim 0.01P_{\text{orb}}$, where $P_{\text{orb}} = 2\pi/\Omega_K$ is the local orbital time-scale.

To estimate $t_c$ in debris disc gas, we assume that the thermal relaxation rate in the disc is of the order of the heating rate. To estimate the latter, we assume that heating is dominated by the dust’s photoelectric effect (Zagorovsky, Brandeker & Wu 2010). This assumption will produce an upper bound on the thermal relaxation time-scale because, if enough carbon gas is in the system, photoionization heating may take over (Kral et al. 2016). We adopt the photoelectric heating rate of equation (23) in Zagorovsky et al. (2010), and divide it by the total internal energy of the gas $\int e dV \sim \int P_{\text{gas}} dV$, where $P_{\text{gas}} = \rho(z)\mathbb{R}T_{\text{gas}}/\mu$ is the gas pressure, $\mathbb{R} = 8.314 \times 10^7 \text{ erg K}^{-1} \text{ g}^{-1}$ is the gas constant, and $\mu = 14$ is the dimensionless mean molecular weight (equipartition between carbon and oxygen). The parameters taken are as follows: $T_{\text{gas}} = 50$ K, $R = 100$ au, $L_{\text{dust}}/L_* = 10^{-3}$, $e\phi = 1$ eV, and the number density of electrons are from Marino et al. (2022)’s model. We leave to future work a more accurate and detailed calculation of the disc thermal physics.

Fig. 2 shows the dimensionless thermal relaxation time-scale $\tau = t_c/P_{\text{orb}}$. In the bulk of the disc $|z| < H$, the time-scale is too long for the VSI to develop except for the most dilute disc in our sample. The densest discs possess mid-plane relaxation times many orders of magnitudes greater than the dynamical time-scale. In all cases, at sufficiently high altitudes the relaxation time reaches a minimum value of $\tau \sim 0.01$, where the VSI may be excited, though if it does appear it will be in a form different from that appearing in
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Thermal relaxation time-scale $\tau = t_c/P_{\text{orb}}$ as a function of height, assuming that cooling rate is balanced by photoelectric heating rate of dust grains. Dashed line is the VSI critical thermal relaxation time-scale for $h = 0.05$, $q = -0.5$, and $\gamma = 1.4$ at 100 au.</figcaption></figure>

protoplanetary discs, and unlikely to transport angular momentum effectively.

Finally, the above discussion is based on pure hydrodynamics, but magnetic fields will negatively impact the onset of the VSI directly, via magnetic tension, and indirectly, via competing MRI turbulence (Latter & Kunz 2022). The MRI will be discussed in Section 5. We also omit the effect of dust drag, which may be important in discs with large dust-to-gas ratios (Kral & Latter 2016) and could provide a strong stabilizing effect (cf. Lin & Youdin 2017; Lehmann & Lin 2022).

## 4.3 Rossby wave instability

The 1D modelling of radial gas spreading in debris discs typically exhibits a bump structure in the surface density. This provides a natural precondition for the RWI. The RWI (Lovelace et al. 1999) is considered one way to generate large-scale vortices, and associated spiral density waves, in the context of accretion discs (Lovelace & Romanova 2014). In barotropic discs, a necessary condition for instability is an extremum in the vortensity profile (Lovelace et al. 1999; Li et al. 2000):

$$ \zeta = \frac{\hat{z} \cdot (\nabla \times \mathbf{v})}{\Sigma_g} = \frac{1}{\Sigma_g R} \frac{\partial(R^2 \Omega)}{\partial R}, \quad (8) $$

where $\Sigma_g = \int \rho dz$ is the surface density.

One prominent feature of the RWI's non-linear evolution is the formation of anticyclonic vortices near the initial surface density extremum (Li et al. 2001). Depending on the most unstable azimuthal modes, multiple vortices will emerge but eventually merge into one single giant vortex (Ono et al. 2016, 2018). These vortices enclose a pressure maximum that can trap dust grains in azimuth. In addition, spiral density waves excited around the vortex can transport angular momentum, with $\alpha \sim 10^{-4}$ to $10^{-2}$ (Li et al. 2001).

Fig. 3 shows the vortensity and gas surface density at different simulation times, using data taken from 1D radial modelling of gas in Marino et al. (2022, their section 4.5), where they consider CO gas released at a rate of 0.1 M$_{\oplus}$ Myr$^{-1}$ in a belt centred at 100

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

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Top panel: vortensity $\zeta$ as a function of radius at different times as evolved forward in the 1D model of Marino et al. (2020). Bottom panel: gas surface density as a function of radius at different times.</figcaption></figure>

and 50 au wide (FWHM), with radial $\alpha = 10^{-2}$. The viscous time-scale is 1.5 Myr. It can be seen that all the $\zeta$ profiles possess a local minimum at $r = 100$ au, even for times much longer than the viscous time, and are hence potentially RWI-unstable. Note, however, that the RWI condition quoted above is often insufficiently stringent (Chang, Youdin & Krapp 2023), especially if there is a radial entropy gradient or the density bump is very wide and viscosity is present. It is likely that the very broad vortensity bumps at late times are stable. Direct stability calculations are needed to truly decide on the prevalence of the RWI in debris discs.

## 5 MAGNETOHYDRODYNAMICS

In this section, we take our fiducial disc models and calculate the ionization fractions and plasma $\beta$ in Section 5.1, as well as the strength of non-ideal MHD effects in Section 5.2. MHD mechanisms that can transport angular momentum include magnetic turbulence
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Top panel: neutral number density as a function of height. Bottom panel: ionization fraction as a function of height, calculated by the ratio of electron number density to neutral number density. As earlier, different colours represent different surface densities.</figcaption></figure>

(MRI) and laminar magnetic stresses, the latter usually associated with a magnetized disc wind. We discuss their application to debris disc gas in Section 5.3.

## 5.1 Ionization and plasma β

### 5.1.1 Ionization fraction

Utilizing the model described in Section 3, Fig. 4 (top) shows the gas density profile, and Fig. 4 (bottom) shows the ionization fraction $x_e = n_e/n$ as functions of vertical height. At low surface densities, the ionization is high even at the mid-plane, with $x_e \sim 10^{-2}$ to $10^{-1}$. For high surface densities, the ionization is low, with $x_e \sim 10^{-7}$ to $10^{-5}$ at the mid-plane; it is expected that non-ideal MHD effects will set in in this environment. The maximum ionization is achieved towards the disc surface, which is $x_e = 0.5$ because carbon atoms are entirely ionized but oxygen atoms remain neutral, due to the lower first ionization potential of carbon.

### 5.1.2 Plasma β

The disc magnetization is parametrized by the plasma β, which is defined as the ratio of the gas pressure to the magnetic pressure:

$$ \beta = \frac{P_{\text{gas}}}{P_{\text{B}}} = \frac{\rho \mathbb{R} T / \mu}{B_z^2 / 8\pi}, \qquad (9) $$

where $\mu = 14$ is the dimensionless mean molecular weight. Fig. 5 shows plasma β as a function of height for different gas surface densities. The left panel is with a magnetic field strength of $B_z = 10^{-6}$ G and right panel of $B_z = 10^{-5}$ G. As the magnetic field strength in debris discs is unconstrained (see discussion in Section 6.2), these values are not necessarily representative, though they do correspond to typical solar interplanetary values. Note that the excitation of linear MRI modes requires $\beta \gtrsim 1$ in ideal MHD (Latter & Kunz 2022).

## 5.2 Non-ideal MHD effects

The ideal MHD regime holds only when gas is fully ionized and the gas is frozen into the magnetic field. Non-ideal MHD physics arises when the gas is partially ionized, which is the case in debris discs (Fig. 4, bottom). The non-ideal MHD effects are manifested in the generalized Ohm's Law. The current density j associated with the electric field E writes (Parks 1991)

$$ j = \sigma \cdot E, \qquad (10) $$

where the conductivity tensor is

$$ \sigma = \begin{pmatrix} \sigma_{\text{P}} & \sigma_{\text{H}} & 0 \\ -\sigma_{\text{H}} & \sigma_{\text{P}} & 0 \\ 0 & 0 & \sigma_{\text{O}} \end{pmatrix}. \qquad (11) $$

The Ohmic, Hall, and Pederson conductivities are denoted by $\sigma_O$, $\sigma_H$, and $\sigma_P$, respectively. Then, the non-ideal diffusivities for Ohmic resistivity, Hall drift, and ambipolar diffusion are written as

$$ \eta_O = \frac{c^2}{4\pi} \left[ \frac{1}{\sigma_O} \right], \quad \eta_H = \frac{c^2}{4\pi} \left[ \frac{\sigma_H}{\sigma_H^2 + \sigma_P^2} \right], \qquad (12) $$

$$ \eta_A = \frac{c^2}{4\pi} \left[ \frac{\sigma_P}{\sigma_H^2 + \sigma_P^2} \right] - \eta_O. \qquad (13) $$

For a three-component fluid composed of electrons, ions, and neutrals, the diffusivities are expressed as (Maeda 1977; Bai 2017)

$$ \eta_O = \frac{m_e c^2}{4\pi e^2} \langle \sigma v \rangle_{\text{en}} x_e^{-1}, \qquad (14) $$

$$ \eta_H = \frac{cB}{4\pi en} x_e^{-1}, \qquad (15) $$

$$ \eta_A = \frac{B^2}{4\pi \langle \sigma v \rangle_{\text{in}} m_n n^2} x_e^{-1}, \qquad (16) $$

where $m_e$ and $m_n$ are the molecular masses of electrons and neutrals, $e$ is elementary charge, $x_e$ is the ionization fraction, and $\langle \sigma v \rangle_{\text{en}}$ and $\langle \sigma v \rangle_{\text{in}}$ are the momentum exchange rate coefficients in electron-neutral and ion-neutral collisions, respectively, where angle brackets represent an average over a Maxwellian distribution function for the relative velocity $v$.

### 5.2.1 Ohmic resistivity

Ohmic resistivity issues from electron-neutral collisions. It can erase magnetic fluctuations on length-scales shorter than $\sim \eta_O/v_A$, where $v_A = B/(4\pi \rho)^{1/2}$ is the Alfvén speed. To inspect the operation of
<figure><img src="image_8.png" /><figcaption>**Figure 5.** The plasma $\beta$ parameter as a function of height. Left: $B_z = 10^{-6}$ G. Right: $B_z = 10^{-5}$ G. Dashed horizontal lines denote the critical value for the onset of MRI turbulence.</figcaption></figure>

the MRI with resistivity, we define the Ohmic magnetic Reynolds number

$$Rm \equiv \frac{c_s H}{\eta_o}. \qquad (17)$$

This quantity is independent of the magnetic field strength and has a minimum $Rm \sim 10^8$ at the mid-plane for the largest surface density model. Local linear analysis finds the condition for the onset of MRI in resistive discs to be $1 \le \beta \le Rm^2$ (Latter & Kunz 2022). In Fig. 5, $B_z = 10^{-6}$ G gives rise to a maximum beta at the mid-plane of $\beta \sim 10^6 \ll Rm^2 \sim 10^{16}$. Thus, Ohmic resistivity is likely to play a negligible role in debris discs.

### 5.2.2 Hall drift

The ion–electron drift is known as the Hall effect. The Hall Elsässer number is defined as

$$Ha \equiv \frac{v_A^2}{\eta_H \Omega}. \qquad (18)$$

The Hall Lundqvist number, independent of the magnetic field strength, is defined as

$$L_H \equiv \frac{v_A H}{\eta_H}. \qquad (19)$$

The effect of Hall drift on disc dynamics depends on the magnetic polarity that co- or counter-aligns with the angular momentum of the disc. However, when $|Ha| > 1$, which is the case for $B_z = 10^{-6}$ G and above, the MRI instability criterion reduces to a slightly modified version of its standard form, $\beta \le 1 + Ha^{-1}$. This can be reframed as

$$\beta \lesssim \frac{1}{4} \left( \sqrt{2L_H^{-1}} + \sqrt{2L_H^{-2} + 4} \right)^2 \qquad (20)$$

(see equation 49 in Latter & Kunz 2022). Fig. 6 demonstrates that the Lundqvist number $L_H \gg 1$ for all our disc models. As a result, the criterion reduces to the ideal MHD limit, where the MRI linear modes are only suppressed when $\beta \le 1$. Observing $\beta$ in Fig. 5, the MRI will be stabilized towards the disc surface.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Hall Lundqvist number as a function of height for different surface densities (colours).</figcaption></figure>

### 5.2.3 Ambipolar diffusion

Charged particles are coupled to magnetic fields and move relative to the neutrals, but suffer from friction due to ion–neutral collisions. Ambipolar diffusion arises from the drift between ions and neutrals. The ambipolar Elsässer number is defined by

$$Am \equiv \frac{v_A^2}{\eta_A \Omega}. \qquad (21)$$

Assuming $m_i \approx m_n$, the Elsässer number can be written as

$$\begin{aligned} Am &\approx \frac{\gamma_i \rho_i}{\Omega} \\ &\approx 4.1 \times 10^2 x_e \left[ \frac{\langle \sigma v \rangle_{\text{in}}}{10^{-9} \text{ cm}^3 \text{ s}^{-1}} \right] \left[ \frac{n}{10^2 \text{ cm}^{-3}} \right], \end{aligned} \qquad (22)$$

where $\rho_i$ and $m_i$ are ion density and ion molecular mass, respectively, and $\gamma_i = \langle \sigma v \rangle_{\text{ni}}/(m_i + m_n)$ is the drag coefficient of neutral–ion collisions. The momentum exchange rate coefficients $\langle \sigma v \rangle_{\text{ni}}$ between
$C^+$ and CO, C I, and O I can be found in Appendix A. Equation (22) is interpreted as the number of collisions between neutrals and ions per dynamical time. Therefore, when Am is smaller (greater) than unity, the ambipolar diffusion is strong (weak). Note that Am is independent of the magnetic field strength. Fig. 7 plots Am versus z for different surface densities and for two different vertical diffusion coefficients. Observing that Am can take significantly low values, we conclude that ambipolar diffusion is the key non-ideal effect in debris disc gas, and one of the main impediments to onset of the MRI. We leave a detailed discussion of its effect to Section 5.3.

## 5.3 MHD mechanisms of angular momentum transport

MHD mechanisms that can transport angular momentum include MRI-induced turbulence, magnetized disc winds, and laminar magnetic stress. Here, we discuss their application to the gas disc based on the calculations carried out in the previous sections.

### 5.3.1 MRI turbulence

The MRI, if present, can drive vigorous turbulence and contribute significantly to the angular momentum transport. The MRI operates best in the ideal MHD limit, where the gas is fully ionized, but is weakened by non-ideal MHD effects. In the previous section, it was demonstrated that Ohmic resistivity and Hall effect are too weak to become dynamically important. Hence, whether the gas disc is MRI-active depends on the strength of ambipolar diffusion and the plasma $\beta$. While local linear theory indicates that the onset of the MRI is unimpeded when $Am > 1$ (e.g. Latter & Kunz 2022), non-linear simulations suggest that fully developed turbulence only appears when $Am \gtrsim 10$ and $\beta > 1$ (Bai & Stone 2011).

Fig. 7 shows ambipolar Elsässer number as a function of height. We find that the mid-plane gas in the various disc models falls into one or the other of two distinct scenarios, depending on the surface densities. For high surface densities, the small $Am \sim 0.1–10$ at the mid-plane can extinguish or at least weaken the MRI, limiting it to higher disc regions with $z > 2 – 3H$ where $Am \sim 100$ and fortunately the plasma $\beta$ is generically above unity. Note, however, that these regions are collisionless and MRI stability will be altered. Note also the influence of strong vertical diffusion in aiding the ionization fraction at the mid-plane (Marino et al. 2022). For low surface densities, the ambipolar Elsässer number is large at the mid-plane, $Am \sim 100$, and thus MRI turbulence is likely, though one still needs to be sure that the plasma $\beta$ is above unity.

It would be useful to give a critical $B_z$ strength that above which the MRI is quenched at all the heights and surface densities considered. As $\beta \propto B^{-2}$, it can be immediately seen from Fig. 5 that if the magnetic field strength is 1 magnitude higher, $B_z = 10^{-4}$ G, the plasma $\beta < 1$ for all z and for all gas densities $\le 1.7 \times 10^{-4}$ g cm$^{-2}$. Higher densities yield ambipolar Elsässer numbers below the critical value for the development of MRI turbulence, $Am < 10$ (Fig. 7; note that Am is independent of B-field strength). It follows that a magnetic field strength of $B_z = 10^{-4}$ G should be considered as the upper limit for the onset of MRI turbulence for the disc models employed.

Importantly, the existence or not of the MRI turbulence will feedback strongly on the input parameters of the disc model, such as vertical diffusion and the radial $\alpha$ parameter. This feedback can work to further enable the MRI if it is already present, or disable it if it is not. For example, if MRI turbulence has begun, vertical diffusion will be larger, thus helping to keep the ionization levels higher (and thus MRI-favourable) at the mid-plane; meanwhile,

radial diffusion of mass will be more efficient, leading to a more dilute disc, and generally higher ionization levels again. On the other hand, if the disc starts off ‘dead’ with low radial and vertical diffusion, ionization may remain below the MRI threshold and the disc stays ‘dead’. This dichotomous behaviour, for the same gas density, might be reflected in hysteresis-type behaviour, though we leave this possibility to future work.

### 5.3.2 Laminar stresses and magnetized disc winds

Numerical work in the protoplanetary disc context shows that under some circumstances a non-trivial laminar MHD state is possible, with strong $R\phi$ magnetic stresses (e.g. Bai 2014; Lesur, Kunz & Fromang 2014; Béthune, Lesur & Ferreira 2017). Interestingly, this accreting state can exist independently of a wind. However, it relies on a significant Hall effect, which we can exclude from debris disc gas (Section 5.2.2). We hence put this transport route to one side.

On the other hand, protoplanetary discs are known to transport angular momentum through magnetized disc winds (Bai & Stone 2013), mainly via their accompanying laminar $z\phi$ magnetic stresses. The operation of magnetized disc winds closely depends on the configuration of the background magnetic field, and can work in the presence of ambipolar diffusion (e.g. Gressel et al. 2015; Cui & Bai 2021). Large-scale, ordered poloidal magnetic fields are critical to launch magnetized disc winds. At the disc surface, a wind is launched and propagates along the magnetic field line, either by a magneto-centrifugal force or by a magnetic pressure gradient (Blandford & Payne 1982; Lynden-Bell 1996). Meanwhile, large-scale correlated azimuthal and poloidal fields exert a torque, by their $z\phi$ Maxwell stress, on the disc surface, which efficiently extracts angular momentum, bestowing it on the gas in the wind.

Determining the background magnetic configuration and strength in debris discs is challenging because these properties are closely linked to the largely uncertain origin of the magnetic field (see discussion in Section 6.2). As we argue later, it is unlikely that there is any large-scale poloidal field threading debris discs, and thus this transport route is probably absent.

# 6 DISCUSSION

## 6.1 Applications to real disc systems

### 6.1.1 Vertical mixing of gas

The theoretical understanding of gas dynamics can shed light on some key issues in debris disc systems. The most fundamental question is the gas origin. It can be either primordial, a leftover from early Class II protoplanetary discs, or secondary, released from solids in the planetesimal belt. The secondary origin is particularly promising for dilute gas discs, characterized by low CO abundance. To explain CO-rich systems by the secondary origin route, gas vertical mixing needs to be weak; that way, C I at the disc surface can effectively shield mid-plane CO. However, if significant vertical mixing is present, the C I shielding mechanism may be rendered ineffective, suggesting instead that CO-rich systems are of primordial origin.

Vertical mixing can be provided by some of the instabilities discussed in this paper, especially if they give rise to fully developed turbulence. The MRI and VSI are possible candidates. MRI turbulence leads to random motions of gas, and current numerical simulations of the VSI generically show large-scale coherent vertical oscillations that are effective at mixing gas in the vertical direction. Importantly, we find that it is the dense CO-rich systems that are most
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Ambipolar Elsässer number as a function of height at different gas surface densities (colours). Left: strong vertical diffusion $\alpha_v = 10^{-2}$. Right: weak vertical diffusion $\alpha_v = 10^{-5}$. Dashed horizontal lines denote the critical value for the onset of MRI turbulence.</figcaption></figure>

likely to be MRI-stable (Section 5.3.1) and thus exhibit the weakest vertical diffusion. In addition, such dense systems are unable to support significant molecular mixing (see Section 4.1). This situation is consistent with the well-shielded scenario and the idea of the secondary origin of CO (Kral et al. 2019). The disc, however, may still be VSI-unstable, but our preliminary conclusion is that this is unlikely for such dense system, though more detailed calculations of the heating and cooling rates need to be performed in order to better determine the onset of VSI.

### 6.1.2 Radial spreading of gas

If the gas evolves viscously, it will spread in radius and, ultimately, accrete on to the central star. Radial gas spreading could contaminate the atmospheres of embedded planets, enhancing their C/O ratios. It is predicted that an Earth's atmospheric mass of gas is readily accreted on terrestrial planets in tenuous gas discs (Kral et al. 2020a). The best constraints on the radial distribution of gas come from resolved ALMA observations of CO (and its isotopologues) and C I gas. The spatial distribution of these two species has only been well constrained for 49 Ceti and $\beta$ Pic.<fnref n="1" /> For 49 Ceti, it was found that CO is tentatively more compact than C I (Higuchi et al. 2019). Kral et al. (2019) used a model to fit the carbon gas distribution in HD 131835 and found tentative evidence that C I extends further inwards than CO.

There are mixed results when comparing the gas distribution (CO or C I) relative to the mm-sized dust, perhaps the best measure of potential spreading of gas from its place of production. In some systems, the gas distribution is more extended (either towards inner or towards outer radii, e.g. HD 21997 and HD 32297; Kóspál et al. 2013; Moór et al. 2013). For others, the gas distribution appears less extended [e.g. 49 Ceti and HD 138813; Hughes et al. 2018; Higuchi et al. 2019; Moór et al. (in preparation)]. In others, they are consistent with each other [e.g. $\beta$ Pic, HD 181327, HD 121617, HD 131488, and HD 131835; Marino et al. 2016; Matrà et al. 2017a, 2019b;

Moór et al. 2017; Kral et al. 2019; Pawellek et al. (in preparation)), suggesting that the gas is not subject to any angular momentum transport process and is quiescent.

These somewhat confounding observations perhaps reflect the complicated nature of angular momentum transport (or its lack) in debris disc systems. While observed gas spreading could issue from any of the transport mechanisms discussed in Section 2.4, we emphasize that the existence or not of these mechanisms is contingent on several, poorly constrained, properties of the system: transport could be magnetic or hydrodynamic, it could be efficient or completely absent. Moreover, there is the added complication that systems could support more than one state for the same set of parameters (see discussion of 'hysteresis' in Section 5.3.1). Though dissatisfying, it does highlight the need to self-consistently connect parameters such as $\alpha$ and vertical mixing to the ionization and other properties of the gas in future modelling.

### 6.1.3 Vortices and dust clumps

The RWI saturates via the formation of vortices (Section 4.3), which in the protoplanetary disc context is a well-studied route to the accumulation of dust grains. In this section, we explore whether a similar effect may occur in debris discs.

The critical issue to resolve here is the strength of the dynamical coupling between dust and gas, which is quantified by the dimensionless Stokes number $St$. When the dust size is less than the gas mean free path, the case for debris discs, the Epstein drag regime is applicable, and $St = \pi \rho_s a \Omega/(\rho c_s) \sim \pi \rho_s a/(2\Sigma_g)$, where $\rho_s$ and $a$ denote the internal density and radius of the dust grains (Weidenschilling 1977). Analytical and numerical studies indicate that dust grains of $St \lesssim 1$ are the most susceptible to trapping in a vortex (Johansen, Andersen & Brandenburg 2004; Surville, Mayer & Lin 2016), though effects peculiar to debris disc gas (e.g. radiation pressure forcing grains into eccentric orbits; see Skaf et al. 2023) may complicate the picture. While generally it is thought that $St > 1$ in debris discs, we now check how dense a disc must be for the dust to feel the effect of any vortex present. A rough criterion is

$$ \Sigma_g \gtrsim 2 \times 10^{-4} \left( \frac{\rho_s}{\text{g cm}^{-3}} \right) \left( \frac{a}{\mu\text{m}} \right) \text{g cm}^{-2}. \quad (23) $$

<footnote n="1">¹There are other systems with resolved observations of CO and C I (e.g. HD 32297), but unfortunately only CO or C I has been analysed or modelled to constrain its spatial distribution.</footnote>
Thus, only rather dense gas discs can couple effectively to the dust.

It so happens that $\beta$ Pic exhibits a prominent dust clump, in southwest (SW) at a stellocentric distance of ~55 au, observed in the mid-IR via continuum thermal emission (Telesco et al. 2005) and comprising µm (but not mm)-sized dust (Telesco et al. 2005; Birnstiel et al. 2018; Matrà et al. 2019b; Sierra & Lizano 2020). Molecular gas clumps have also been detected in CO and C I emission, likely co-located with the dust clump (Dent et al. 2014; Matrà et al. 2017a; Cataldi et al. 2018). Several physical mechanisms have been proposed in the literature to explain this feature, including an origin via a RWI-induced vortex (Skaf et al. 2023). However, given the low gas surface densities in $\beta$ Pic, which are estimated to be $\Sigma_g \sim 10^{-7}$ to $10^{-6}$ g cm$^{-2}$ (Kral et al. 2016), the Stokes number for µm grains is just too high. It should be emphasized that $\beta$ Pic is encircled by a rather tenuous gas disc; denser gas systems, subject to the RWI, may be more capable of interacting and imprinting clumpy structure on the dust component.

## 6.2 Origin of the background magnetic field

The presence, strength, and configuration of any background large-scale magnetic fields are critical for the dynamics of the gas, if it is sufficiently well ionized (see Sections 2.4 and 5). To begin, the MRI's onset, saturation, and angular momentum transport are linked to the existence or not of large-scale net magnetic flux threading the gaseous disc (e.g. Zhu & Stone 2018; Jacquemin-Ide, Lesur & Ferreira 2021). On the other hand, large-scale laminar fields can lead to quasi-steady states that drive significant angular momentum transport, usually alongside an MHD wind, which are favoured if non-ideal effects kill off the MRI (e.g. Bai & Stone 2013; Lesur et al. 2014; Gressel et al. 2015; Béthune et al. 2017; Cui & Bai 2021). Presently, we are only just beginning to constrain debris disc magnetism (Hull et al. 2022), and so there are no observational constraints on its magnetic field properties. In this section, we instead discuss plausible origin scenarios for a background field and what its properties might be.

First, we consider the interplanetary magnetic field permeating a debris disc system. In the Solar system, the Cassini and Voyager 1 and 2 reported an interplanetary magnetic field of ~1–10 nT, equivalent to 1–10 µG (e.g. Jackman et al. 2005; Burlaga et al. 2022; Tsurutani et al. 2023), which is frozen into the solar wind. While such interplanetary magnetic fields exhibit a large-scale structure (see Owens & Forsyth 2013), debris disc gas will not feel their effects unless they are well penetrated by the stellar wind plasma. This might be possible in dilute systems, but then the disc gas is likely to be simply blown away in a 'belt wind' arising from the interaction (Kral et al. 2023). In denser systems, it is unclear how well the stellar wind plasma can mix with the debris disc gas; while this may be facilitated by instabilities and turbulence at the disc's bow shock, it seems unlikely that the wind's large-scale magnetic structure can be imprinted on the disc in any coherent way, and certainly not in the form demanded by current protoplanetary disc wind modelling (see previous paragraph). That said, the stellar wind should be able to provide a fund of smaller scale magnetic fluctuations that might serve as the seeds for local MHD activity.

Another possibility is that a background magnetic field might be inherited from the preceding protoplanetary disc stage, with the field 'carried over' to the later debris stage by residual ionized gas, if it exists, or any charged grains that evade coagulation into pebbles and planetesimals (Umebayashi & Nakano 1980; Ciolek & Mouschovias 1993). The viability of this process, however, is yet to be demonstrated. In any case, the final stages of a protoplanetary disc might witness the effective evacuation of most or even all of its

large-scale magnetic flux via ambipolar-assisted disc winds (see e.g. Lesur 2021). One is hence tempted to conclude that debris discs will not possess a coherent large-scale magnetic field of any significance via this route.

In summary, it is likely that any magnetic field within debris disc gas is internally generated from a small seed, via a dynamo process of some kind, most likely via a form of the MRI. We assume that this dynamo is unable to construct a large-scale magnetic field capable of producing meaningful wind accretion, though this may need to be checked (Lesur, private communication). On the other hand, we expect an internally generated field to provide further protection from penetration/erosion by the stellar wind, an effect that may work against the launching of belt winds (Kral et al. 2023).

## 6.3 Ionization sources

In this section, we flag a caveat arising from the calculation of the ionization fractions in our debris disc gas model. Previous work on gas evolution includes only the photoionization of gas by an external stellar and interstellar UV radiation field (e.g. Kral et al. 2016; Marino et al. 2022). We note that other processes may be important, e.g. the photoelectric effect of dust grains, cosmic rays, and stellar winds, and give them a brief discussion below.

Dust grains are an essential player in the ionization of protoplanetary discs; free electrons are produced by photoionization and are absorbed on to grains' surfaces, leading to negatively charged dust (Umebayashi & Nakano 1980). Tiny submicron-sized grains are more efficient in capturing electrons because of their larger total surface density. However, in debris discs the smallest grains (~0.1–10 µm, depending on the stellar mass) are blown out due to radiation pressure (Krivov, Löhne & Sremčević 2006), thus significantly diminishing the capture of electrons by this route. Moreover, in debris discs, electrons are more prone to recombine on to ionized carbon atoms rather than to be absorbed by grains when the gas ionization fraction is relatively high, as tested by Cloudy simulations of $\beta$ Pic (Kral et al. 2016).

Photoelectric charging could also play a role as the dust-to-gas density ratio is greater than unity in debris discs. It is found that photoelectric charging, when balanced by the thermal electron collection current, results in positively charged grains (e.g. Fernández, Brandeker & Wu 2006; Besla & Wu 2007; Zagorovsky et al. 2010). For massive gas discs, we expect that both photoelectric and photoionization processes will weaken, because they both depend on the same UV penetrating photons. A detailed study of ionization, including dust, in debris discs shall be conducted in the future, on a case-by-case basis, as it strongly depends on the gas versus dust abundances and radiation sources.

Direct and indirect cosmic ray ionization of carbon can be important for denser discs, in which UV photons are less able to penetrate. We present a detailed discussion in Appendix B. Its main conclusion is that, unless there are significant levels of primordial $H_2$, cosmic rays increase the ambipolar Elsasser number by a factor of 5 at most, and then only in the densest discs and with efficient vertical mixing. Lastly, high-velocity stellar wind protons can turn a small fraction of CO to $CO^+$ (Kral et al. 2021, 2023).

## 6.4 Future MHD simulations

Global numerical simulations should be performed in the context of debris discs in the future, in order to understand the contribution of MRI to the angular momentum transport, turbulence, radial spreading, and vertical mixing of gas. The initial magnetic field
configuration is key to understanding the gas dynamics, especially when it comes to determining the launching of magnetized disc winds. Initial configuration of large-scale open magnetic fields ensures winds due to flux freezing at and above the disc surface, while poloidal magnetic loops would hinder the launching of winds. In the context of protoplanetary discs, global simulations incorporating non-ideal MHD effects and magnetized disc winds have been conducted, employing large-scale net poloidal magnetic fields (e.g. Bai 2017; Béthune et al. 2017; Wang, Bai & Goodman 2019; Gressel et al. 2020; Cui & Bai 2021), as inherited from the primordial molecular cloud (e.g. Galli & Shu 1993; Girart, Rao & Marrone 2006; Girart et al. 2009). Future simulations for debris discs could test distinct initial magnetic field set-ups, since it is demanding to constrain the magnetic field strength and configuration. Future observational constraints on the gas distribution and kinematics would be helpful to probe the magnetic field properties of debris discs.

## 7 CONCLUSION

We review five mechanisms that can transport angular momentum in the context of debris disc gas: molecular viscosity, hydrodynamic turbulence, MHD turbulence, magnetized disc winds, and laminar magnetic stress. Disc models from Marino (2022) are employed to describe the spatial distribution of the gas species, including CO, C I, C II, and O I (Section 3), and thus provide us with quantitative estimates. We summarize the main findings as follows:

(i) The molecular viscosity sets a non-negligible lower limit for angular momentum transport in debris disc gas. The molecular-driven, vertically averaged $\alpha$ ranges from ~$10^{-5}$ to ~0.1 from high to low gas surface densities (Section 4.1). In addition, molecular diffusion may contribute significantly to vertical mixing in more dilute systems, as we expect the kinetic $\alpha_v$ to be ~$\alpha$.

(ii) The VSI is unlikely to drive hydrodynamic turbulence. Only in especially dilute gas regions is the thermal relaxation time sufficiently rapid, but the gas is reasonably well ionized and the MRI and/or magnetic tension will intervene in the VSI's development (Section 4.2).

(iii) In the secondary origin scenario, the gas disc naturally exhibits a bump structure in the surface density, at least early in its evolution. This may generate gas vortices via the RWI, but a sufficiently large gas surface density is needed to ensure $St \lesssim 1$ and hence the formation of dust clumps (Sections 4.3 and 6.1.3).

(iv) The MRI is controlled by the ambipolar diffusion strength, while Ohmic resistivity and Hall effect are too weak to make a significant impact in our disc model. More dilute gas discs should be MRI turbulent at all heights, while in denser discs the MRI may be confined to the better ionized upper layers $z > 2 - 3H$, though these conclusions rely on a background magnetic field $B_z = 10^{-6}$ to $10^{-5}$ G, which is poorly constrained. However, we do argue that above $B_z = 10^{-4}$ G the MRI will be quenched in all disc models we consider (Section 5.3.1).

(v) Laminar $R\phi$ magnetic stresses can transport angular momentum but rely on a significant Hall effect, and are hence ruled out in the debris disc context. Laminar $z\phi$ magnetic stresses rely on a magnetized disc wind, which in turn requires a large-scale vertical magnetic field (Section 5.3.2) threading the disc. We are sceptical that such a B-field configuration is possible in debris discs, though future investigation is necessary to decide on this point (Section 6.2).

(vi) Of the angular momentum transport mechanisms reviewed, MRI turbulence appears the most robust, but it is only effective in more dilute, and hence better ionized, discs. In denser systems,

it is likely that there is little transport at all. This suggests that gaseous debris discs are dichotomous. If a disc begins dilute and MRI turbulent, then, via effective vertical mixing and radial spreading, it will remain dilute, well ionized, and thus MRI turbulent. If a disc begins dense, poorly ionized, and quiescent, it will remain so; in fact, mass may continue to build up, as it will be inefficiently spread radially (Section 5.3.1). This dependence on the disc's dynamical history (and on its poorly constrained physical parameters) perhaps explains the extant, somewhat confounding, results on gas spreading (Section 6.1.2).

## ACKNOWLEDGEMENTS

We thank the anonymous referee for the detailed comments that improved the presentation of this paper. We thank Daniele Galli for pointing out where we could obtain the ionization rate of C I by cosmic rays. We also thank Mark Wyatt and Ya-Ping Li for useful discussions at the earlier stage of this project. CC and HL acknowledge funding from STFC grant ST/T00049X/1. SM was supported by a Royal Society University Research Fellowship (URF-R1-221669).

## DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

## REFERENCES

Armitage P. J., 2011, ARA&A, 49, 195

Bai X.-N., 2014, ApJ, 791, 137

Bai X.-N., 2017, ApJ, 845, 75

Bai X.-N., Stone J. M., 2011, ApJ, 736, 144

Bai X.-N., Stone J. M., 2013, ApJ, 769, 76

Besla G., Wu Y., 2007, ApJ, 655, 528

Béthune W., Lesur G., Ferreira J., 2017, A&A, 600, A75

Birnstiel T. et al., 2018, ApJ, 869, L45

Blandford R. D., Payne D. G., 1982, MNRAS, 199, 883

Bonsor A., Wyatt M. C., Marino S., Davidsson B. J. R., Kral Q., Thebault P., 2023, MNRAS, 526, 3115

Burlaga L. F., Ness N. F., Berdichevsky D. B., Jian L. K., Kurth W., Park J., Rankin J., Szabo A., 2022, ApJ, 932, 59

Cataldi G. et al., 2014, A&A, 563, A66

Cataldi G. et al., 2018, ApJ, 861, 72

Cataldi G. et al., 2020, ApJ, 892, 99

Cataldi G. et al., 2023, ApJ, 951, 111

Chang E., Youdin A. N., Krapp L., 2023, ApJ, 946, L1

Ciolek G. E., Mouschovias T. C., 1993, ApJ, 418, 774

Cui C., Bai X.-N., 2020, ApJ, 891, 30

Cui C., Bai X.-N., 2021, MNRAS, 507, 1106

Cui C., Bai X.-N., 2022, MNRAS, 516, 4660

Cui C., Latter H. N., 2022, MNRAS, 512, 1639

Cui C., Lin M.-K., 2021, MNRAS, 505, 2983

Czechowski A., Mann I., 2007, ApJ, 660, 1541

Dang Y., Cui C., Barraza-Alfaro M., 2024, MNRAS, 529, 918

Dent W. R. F. et al., 2014, Science, 343, 1490

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

Fernández R., Brandeker A., Wu Y., 2006, ApJ, 643, 509

Flock M., Nelson R. P., Turner N. J., Bertrang G. H. M., Carrasco-González C., Henning T., Lyra W., Teague R., 2017, ApJ, 850, 131

Galli D., Shu F. H., 1993, ApJ, 417, 220

Girart J. M., Rao R., Marrone D. P., 2006, Science, 313, 812

Girart J. M., Beltrán M. T., Zhang Q., Rao R., Estalella R., 2009, Science, 324, 1408
Goldreich P., Tremaine S. D., 1978, Icarus, 34, 227

Gressel O., Turner N. J., Nelson R. P., McNally C. P., 2015, ApJ, 801, 84

Gressel O., Ramsey J. P., Brinch C., Nelson R. P., Turner N. J., Bruderer S., 2020, ApJ, 896, 126

Hales A. S. et al., 2022, ApJ, 940, 161

Heays A. N., Bosman A. D., van Dishoeck E. F., 2017, A&A, 602, A105

Henning T., Semenov D., 2013, Chem. Rev., 113, 9016

Henry R. J. W., Burke P. G., Sinfailam A.-L., 1969, Phys. Rev., 178, 218

Higuchi A. E. et al., 2019, ApJ, 883, 180

Higuchi A. E., Kóspál Á., Moór A., Nomura H., Yamamoto S., 2020, ApJ, 905, 122

Holland W. S. et al., 1998, Nature, 392, 788

Huang H., Ormel C. W., 2023, MNRAS, 522, 2241

Hughes A. M. et al., 2017, ApJ, 839, 86

Hughes A. M., Duchêne G., Matthews B. C., 2018, ARA&A, 56, 541

Hull C. L. H. et al., 2022, ApJ, 930, 49

Itikawa Y., 2015, J. Phys. Chem. Ref. Data, 44, 013105

Itikawa Y., Ichimura A., 1990, J. Phys. Chem. Ref. Data, 19, 637

Jackman C. M., Achilleos N., Bunce E. J., Cecconi B., Clarke J. T., Cowley S. W. H., Kurth W. S., Zarka P., 2005, J. Geophys. Res.: Space Phys., 110, A10212

Jacquemin-Ide J., Lesur G., Ferreira J., 2021, A&A, 647, A192

Johansen A., Andersen A. C., Brandenburg A., 2004, A&A, 417, 361

Kennedy G. M., Marino S., Matrà L., Panić O., Wilner D., Wyatt M. C., Yelverton B., 2018, MNRAS, 475, 4924

Klusmeyer J. et al., 2021, ApJ, 921, 56

Kóspál Á., et al., 2013, ApJ, 776, 77

Kral Q., Latter H., 2016, MNRAS, 461, 1614

Kral Q., Wyatt M., Carswell R. F., Pringle J. E., Matrà L., Juhász A., 2016, MNRAS, 461, 845

Kral Q., Matrà L., Wyatt M. C., Kennedy G. M., 2017, MNRAS, 469, 521

Kral Q., Marino S., Wyatt M. C., Kama M., Matrà L., 2019, MNRAS, 489, 3670

Kral Q., Davoult J., Charnay B., 2020a, Nat. Astron., 4, 769

Kral Q., Matrà L., Kennedy G. M., Marino S., Wyatt M. C., 2020b, MNRAS, 497, 2811

Kral Q. et al., 2021, A&A, 653, L11

Kral Q., Pringle J. E., Matrà L., Thébault P., 2023, A&A, 669, A116

Krivov A. V., Löhne T., Sremčević M., 2006, A&A, 455, 509

Krivov A. V., Herrmann F., Brandeker A., Thébault P., 2009, A&A, 507, 1503

Landau L. D., Lifshitz E. M., 1959, Fluid Mechanics. Pergamon Press, Oxford

Latter H. N., Kunz M. W., 2022, MNRAS, 511, 1182

Latter H. N., Papaloizou J., 2018, MNRAS, 474, 3110

Lehmann M., Lin M.-K., 2022, A&A, 658, A156

Lesur G., 2021, J. Plasma Phys., 87, 205870101

Lesur G., Kunz M. W., Fromang S., 2014, A&A, 566, A56

Lesur G. et al., 2023, in Inutsuka S., Aikawa Y., Muto T., Tomida K., Tamura M., eds, ASP Conf. Ser. Vol. 534, Protostars and Planets VII. Astron. Soc. Pac., San Francisco, p. 465

Li H., Finn J. M., Lovelace R. V. E., Colgate S. A., 2000, ApJ, 533, 1023

Li H., Colgate S. A., Wendroff B., Liska R., 2001, ApJ, 551, 874

Lieman-Sifry J., Hughes A. M., Carpenter J. M., Gorti U., Hales A., Flaherty K. M., 2016, ApJ, 828, 25

Lin M.-K., Youdin A. N., 2015, ApJ, 811, 17

Lin M.-K., Youdin A. N., 2017, ApJ, 849, 129

Lovelace R. V. E., Romanova M. M., 2014, Fluid Dyn. Res., 46, 041401

Lovelace R. V. E., Li H., Colgate S. A., Nelson A. F., 1999, ApJ, 513, 805

Lynden-Bell D., 1996, MNRAS, 279, 389

Lyra W., Kuchner M., 2013, Nature, 499, 184

Maeda K.-I., 1977, J. Atmos. Terr. Phys., 39, 1041

Marino S., 2022, preprint (arXiv:2202.03053)

Marino S. et al., 2016, MNRAS, 460, 2933

Marino S., Flock M., Henning T., Kral Q., Matrà L., Wyatt M. C., 2020, MNRAS, 492, 4409

Marino S., Cataldi G., Jankovic M. R., Matrà L., Wyatt M. C., 2022, MNRAS, 515, 507

Maroulis G., 1996, J. Phys. Chem., 100, 13466

Matrà L. et al., 2017a, MNRAS, 464, 1415

Matrà L. et al., 2017b, ApJ, 842, 9

Matrà L., Öberg K. I., Wilner D. J., Olofsson J., Bayo A., 2019a, AJ, 157, 117

Matrà L., Wyatt M. C., Wilner D. J., Dent W. R. F., Marino S., Kennedy G. M., Milli J., 2019b, AJ, 157, 135

Matthews B. C., Krivov A. V., Wyatt M. C., Bryden G., Eiroa C., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. University of Arizona Press, Tucson, p. 521

Moór A. et al., 2011, ApJ, 740, L7

Moór A. et al., 2013, ApJ, 777, L25

Moór A. et al., 2017, ApJ, 849, 123

Nakatani R., Liu H. B., Ohashi S., Zhang Y., Hanawa T., Chandler C., Oya Y., Sakai N., 2020, ApJ, 895, L2

Nelson R. P., Gressel O., Umurhan O. M., 2013, MNRAS, 435, 2610

Olofsson J. et al., 2022, MNRAS, 513, 713

Ono T., Muto T., Takeuchi T., Nomura H., 2016, ApJ, 823, 84

Ono T., Muto T., Tomida K., Zhu Z., 2018, ApJ, 864, 70

Owens M. J., Forsyth R. J., 2013, Living Rev. Sol. Phys., 10, 5

Parks G. K., 1991, Physics of Space Plasmas: An Introduction. Addison-Wesley Publishing Co., Redwood City, CA

Prasad S. S., Huntress W. T. J., 1980, ApJS, 43, 1

Pringle J. E., 1981, ARA&A, 19, 137

Rebollido I. et al., 2022, MNRAS, 509, 693

Roberge A., Feldman P. D., Weinberger A. J., Deleuil M., Bouret J.-C., 2006, Nature, 441, 724

Schmidt J., Ohtsuki K., Rappaport N., Salo H., Spahn F., 2009, in Dougherty M. K., Esposito L. W., Krimigis S. M., eds, Saturn from Cassini-Huygens. Springer, Dordrecht, p. 413

Schwerdtfeger P., Nagle J. K., 2019, Mol. Phys., 117, 1200

Sierra A., Lizano S., 2020, ApJ, 892, 136

Skaf N. et al., 2023, A&A, 675, A35

Smirnov-Pinchukov G. V., Moór A., Semenov D. A., Ábrahám P., Henning T., Kóspál Á., Hughes A. M., di Folco E., 2022, MNRAS, 510, 1148

Stoll M. H. R., Kley W., 2016, A&A, 594, A57

Størm P. A. et al., 2020, PASP, 132, 101001

Surville C., Mayer L., Lin D. N. C., 2016, ApJ, 831, 82

Svanberg E., Cui C., Latter H. N., 2022, MNRAS, 514, 4581

Takeuchi T., Artymowicz P., 2001, ApJ, 557, 990

Telesco C. M. et al., 2005, Nature, 433, 133

Thomas L. D., Nesbet R. K., 1975, Phys. Rev. A, 12, 1729

Tsurutani B. T. et al., 2023, IEEE Trans. Plasma Sci., 51, 1595

Umebayashi T., Nakano T., 1980, PASJ, 32, 405

Umebayashi T., Nakano T., 1981, PASJ, 33, 617

Wang L., Bai X.-N., Goodman J., 2019, ApJ, 874, 90

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

Weiss B. P., Bai X.-N., Fu R. R., 2021, Sci. Adv., 7, eaba5967

Wyatt M. C., 2008, ARA&A, 46, 339

Zagorovsky K., Brandeker A., Wu Y., 2010, ApJ, 720, 923

Zhu Z., Stone J. M., 2018, ApJ, 857, 34

Zuckerman B., Song I., 2012, ApJ, 758, 77

## APPENDIX A: MOMENTUM EXCHANGE RATE COEFFICIENTS

The momentum exchange rate coefficients of electron-neutral and ion-neutral collisions are calculated here in order to compute Ohmic and ambipolar Elsässer number.

Electron-neutral scattering limits the electrical conductivity in low-ionization regions. The two-body collision rate coefficient is (Draine 2011)

$$ \langle \sigma v \rangle = \left[ \frac{8k_B T}{\pi \mu} \right]^{1/2} \int_0^\infty \sigma(E) \frac{E}{k_B T} e^{-E/k_B T} \frac{dE}{k_B T} \quad (A1) $$

where $k_B$ is the Boltzmann constant, $\mu$ the reduced mass of the two collision species, $T$ the thermal equilibrium temperature, $E$ the centre-
<figure><img src="image_11.png" /><figcaption>Figure A1. The integration part of equation (A1) for C, O, and CO.</figcaption></figure>

of-mass energy, and $\sigma$ the reaction cross-section. The momentum exchange cross-sections $\sigma_{\text{mt}}(E)$ for e$^-$–O and e$^-$–CO collisions are tabulated in Thomas & Nesbet (1975), Itikawa & Ichimura (1990), and Itikawa (2015). The elastic scattering cross-section for e$^-$–C is read in Henry, Burke & Sinfailam (1969). Note that the cross-sections measured in laboratory experiments are usually given as a function of energy in the laboratory frame. For electron-neutral collisions, $E \approx E_{\text{lab}}$.

In Fig. A1, the integration part of equation (A1) $\tilde{\sigma} = \int_0^\infty \sigma(E) \frac{E}{k_B T} e^{-E/k_B T} \frac{dE}{k_B T}$ is shown for C, O, and CO. When 10 K $\le T \le 1000$ K, $\tilde{\sigma}_{e^-,C} = 2 \times 10^{-15} \text{ cm}^2$, $\tilde{\sigma}_{e^-,O} = 2 \times 10^{-16} \text{ cm}^2$, and $\tilde{\sigma}_{e^-,CO} = 6.5 \times 10^{-16} \text{ cm}^2$. The corresponding rate coefficients can then be computed by

$$ \langle \sigma v \rangle_{e^-,C} = 4.14 \times 10^5 T^{1/2} \tilde{\sigma} \text{ cm}^3 \text{ s}^{-1}, \quad (A2) $$

and at $T = 50$ K,

$$ \langle \sigma v \rangle_{e^-,C} = 5.85 \times 10^{-9} \text{ cm}^3 \text{ s}^{-1}, \quad (A3) $$

$$ \langle \sigma v \rangle_{e^-,O} = 5.85 \times 10^{-10} \text{ cm}^3 \text{ s}^{-1}, \quad (A4) $$

$$ \langle \sigma v \rangle_{e^-,CO} = 1.9 \times 10^{-9} \text{ cm}^3 \text{ s}^{-1}. \quad (A5) $$

The ion-neutral collision rates determine the ambipolar diffusion strength. The momentum transfer rate coefficients for ion-neutral collisions between C$^+$ and C, O, and CO are calculated by equations (2.34) and (2.39) of Draine (2011), with the tabulated dipole polarizabilities (Maroulis 1996; Schwerdtfeger & Nagle 2019). This yields

$$ \langle \sigma v \rangle_{C^+,C} = 9.57 \times 10^{-10} \text{ cm}^3 \text{ s}^{-1}, \quad (A6) $$

$$ \langle \sigma v \rangle_{C^+,O} = 1.49 \times 10^{-9} \text{ cm}^3 \text{ s}^{-1}, \quad (A7) $$

$$ \langle \sigma v \rangle_{C^+,CO} = 1.35 \times 10^{-9} \text{ cm}^3 \text{ s}^{-1}. \quad (A8) $$

## APPENDIX B: COSMIC RAY IONIZATION OF CARBON

For tenuous debris discs such as β Pic, stellar and interstellar UV should dominate the ionization processes. Meanwhile, denser discs may become optically thick to ionizing UV photons at some height. In these circumstances, cosmic ray ionization might be significant and may be the leading ionization source towards the mid-plane (Umebayashi & Nakano 1981). To test this, we focus on the direct and indirect cosmic ray ionization of carbon as it has a lower first ionization potential than oxygen.

We implemented this new source of ionization for carbon in exogas<fnref n="2" /> (Marino et al. 2022), and evaluate two cases. First, we treat the case where the gas is of secondary origin and thus there is no H₂. Here, the dominant effect of cosmic rays is via the direct ionization of carbon at a rate of 1.94 relative to the standard value for H₂ (i.e. a C I ionization rate of $1.94 \times 10^{-16} \text{ s}^{-1}$; Prasad & Huntress 1980; Heays, Bosman & van Dishoeck 2017). Secondly, we treat the case where the gas is of a primordial origin and there are orders of magnitude more H₂ than C I. Now, the dominant effect of cosmic rays on C I is the indirect ionization via the excitation of H₂ and H followed by the emission of UV photons capable of ionizing C I. In this case, we use the tabulated value in Heays et al. (2017) of $2.6 \times 10^{-14} \text{ s}^{-1}$. Note that this value depends on the abundance of dust and C I and it assumes that the disc is optically thick to ionizing UV photons (valid in the primordial scenario). While the abundance of dust and C I could be different than in the context studied by Heays et al. (2017, C I abundance of $10^{-5}$ to $10^{-4}$), in the absence of a better estimate we use their value.

Fig. B1 shows the ambipolar Elsässer number Am for direct and indirect cosmic ray ionization. For strong vertical mixing ($\alpha_v = 10^{-2}$), the impact of cosmic rays is noticeable (left panel). Direct cosmic ray ionization yields a maximum amplification of ionization fraction by a factor of ~5 for dense discs at the mid-plane. Indirect cosmic ray ionization (due to the excitation of H₂) has a stronger effect, a factor of ~30–40 increase of Am at the mid-plane. It indicates that the entire disc can be MRI turbulent. Of course, indirect ionization requires the presence of large amounts of H₂, which would only be the case if the gas were of primordial origin.

For weak vertical mixing ($\alpha_v = 10^{-5}$), the impact of cosmic ray is limited. As mixing is inefficient, the amount of C I in the mid-plane is low, which lowers the ionization fraction and Am in the mid-plane. For dilute discs, the ionization by UV photons dominates over the cosmic ray, and hence the three curves almost overlap for both weak and strong vertical diffusion.

<footnote n="2">²https://github.com/SebaMarino/exogas</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure B1.** Ambipolar Elsässer number *Am* for direct and indirect cosmic ray ionization. Solid line: no cosmic ray ionization. Dashed line: direct cosmic ray ionization. Dotted line: indirect cosmic ray ionization of carbon by the excitation of H₂ and H. Left: strong vertical diffusion $\alpha_v = 10^{-2}$. Right: weak vertical diffusion $\alpha_v = 10^{-5}$. Different colours indicate different surface densities.</figcaption></figure>

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