# Simulations of Protoplanetary Disk Dispersal: Stellar Mass Dependence of the
Disk Lifetime

Ayano Komaki Department of Physics, The University of Tokyo, 7-3-1 Hongo,
Bunkyo, Tokyo 113-0033, Japan<EMAIL_ADDRESS>Shuhei
Fukuhara  Department of Multi-Disciplinary Sciences, Graduate School of Arts
and Sciences, The University of Tokyo, 3-8-1 Komaba, Meguro, Tokyo 153-8902,
Japan  Takeru K. Suzuki  Department of Multi-Disciplinary Sciences, Graduate
School of Arts and Sciences, The University of Tokyo, 3-8-1 Komaba, Meguro,
Tokyo 153-8902, Japan  Komaba Institute for Science, The University of Tokyo,
3-8-1 Komaba, Meguro, Tokyo 153-8902, Japan  Department of Astronomy, Graduate
School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo
113-0033, Japan  Naoki Yoshida Department of Physics, The University of Tokyo,
7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan  Kavli Institute for the Physics
and Mathematics of the Universe (WPI), UT Institute for Advanced Study, The
University of Tokyo, Kashiwa, Chiba 277-8583, Japan  Research Center for the
Early Universe (RESCEU), School of Science, The University of Tokyo, 7-3-1
Hongo, Bunkyo, Tokyo 113-0033, Japan

###### Abstract

Recent infrared and submillimeter observations suggest that the protoplanetary
disk lifetime depends on the central stellar mass. The disk dispersal is
thought to be driven by viscous accretion, magneto-hydrodynamics (MHD) winds,
and photoevaporation by the central star. We perform a set of one-dimensional
simulations of long-term disk evolution that include all the three processes.
We vary the stellar mass in the range of 0.5–7${\rm\,M_{\odot}}$, and study
the mass dependence of the disk evolution. We show that a significant fraction
of the disk gas is lost by MHD winds in the early stage, but the later disk
evolution is mainly governed by photoevaporation. The disk radius decreases as
photoevaporation clears out the gas in the outer disk efficiently. The
qualitative evolutionary trends of the disk mass are remarkably similar for
the wide range of the central stellar mass we consider, and the time evolution
of the disk mass can be well fitted by a simple function. The dispersal time
is approximately ten million years for low mass stars with weak mass
dependence, but gets as short as two million years around a
$7{\rm\,M_{\odot}}$ star. In the latter case, a prominent inner hole is formed
by the combined effect of accretion and MHD winds within about one million
years. The strength of the MHD wind and viscous accretion controls the overall
mass-loss rate, but does not alter the dependence of the dispersal timescale
on the central stellar mass.

††journal: ApJ

##  1 Introduction

Recent observations discovered diverse planetary architectures around various
types of stars (Fulton et al., 2017; Zhang et al., 2018). Planets are formed
in protoplanetary disks (PPDs) out of disk materials within some limited time,
and thus the dynamical evolution of a PPD affects significantly planet
formation. In particular, disk dispersal process and its characteristic
timescale are thought to be critically important for setting the scene for the
formation of diverse planets.

Observations of several star forming regions suggest that PPDs have a lifetime
of a few million years (e.g., Haisch et al., 2001; Meyer et al., 2007;
Hernández et al., 2007; Mamajek, 2009; Bayo et al., 2012; Ribas et al., 2014).
It is also found that the disk lifetime may depend on the central stellar mass
(e.g., Carpenter et al., 2006; Lada et al., 2006; Allers et al., 2007; Dahm &
Hillenbrand, 2007; Kennedy & Kenyon, 2009; Fang et al., 2012; Yasui et al.,
2014; Ribas et al., 2015). For example, Ribas et al. (2015) sorted disks by
the system age and the central stellar mass, and calculated the disk fraction
of each sample. Based on the finding that the disk fraction decreased as the
stellar mass increased, they concluded that disks around high-mass stars had
shorter lifetimes.

The three major disk dispersal mechanisms proposed so far are viscous
accretion onto the central star, magnetohydrodynamic (MHD) winds, and
photoevaporation. In a PPD, efficient angular momentum transfer causes the gas
of the inner disk to fall onto the star. The accretion rate is characterised
by $\alpha$ parameter (Shakura & Sunyaev, 1973; Lynden-Bell & Pringle, 1974),
which can be estimated in several ways based on observations (Calvet et al.,
2005; Fedele et al., 2010; Mathews et al., 2012). Interestingly, a positive
correlation is suggested between the stellar mass and the accretion rate
(Muzerolle et al., 2003; Hartmann et al., 2006; Herczeg & Hillenbrand, 2008;
Fairlamb et al., 2015; Hartmann et al., 2016).

The role of MHD disk winds in the evolution of PPDs has attracted considerable
attention in recent years (e.g., Pascucci et al., 2022). In addition to the
direct mass loss (Suzuki & Inutsuka, 2009), MHD disk winds remove angular
momentum from the disk, which promotes mass accretion (Blandford & Payne,
1982; Bai, 2016; Béthune et al., 2017; Gressel et al., 2020). As a result, the
radial profile of surface density is significantly altered from the one
without the effect of MHD disk winds (Suzuki et al., 2016; Hasegawa et al.,
2017; Tabone et al., 2022). MHD-driven winds also cause a great impact on the
evolution of solid particles in PPDs (Taki et al., 2021; Arakawa et al., 2021)
and the formation and migration of planets (Ogihara et al., 2018; Kimmig et
al., 2020).

Photoevaporation is driven by high-energy radiation such as far-ultraviolet
(FUV; $6{\rm\,eV}\lesssim h\nu<13.6{\rm\,eV}$), extreme-ultraviolet (EUV;
$13.6{\rm\,eV}<h\nu\lesssim 100{\rm\,eV}$), and X-ray ($100{\rm\,eV}\lesssim
h\nu\lesssim 10{\rm\,keV}$) photons, which heat the gas on the disk surface to
launch photoevaporative flows. The effect of photoevaporation on the disk
evolution has been studied by a number of authors (Hollenbach et al., 1994;
Richling & Yorke, 1997; Gorti & Hollenbach, 2009; Ercolano et al., 2009; Owen
et al., 2010; Tanaka et al., 2013; Wang & Goodman, 2017; Nakatani et al.,
2018a, b; Picogna et al., 2019; Komaki et al., 2021; Picogna et al., 2021).
Gorti & Hollenbach (2009) performed 1+1 dimensional simulations of disk
photoevaporation by varying the central stellar mass in the range of
0.3–7.0${\rm\,M_{\odot}}$. They incorporated both accretion and
photoevaporation as major disk dispersal processes. Komaki et al. (2021) used
radiation hydrodynamics simulations to show that the mass-loss rate by
photoevaporation increased as the central stellar mass.

It is important to notice that disk dispersal may likely be caused by a
combination of multiple physical mechanisms, between which there can be
complicated interplay, and that the whole dispersal process can last over a
long time of several to ten million years. Clearly, it is necessary to study
the long-term disk evolution with incorporating all the proposed dispersal
mechanisms in a consistent manner. This motivated several theoretical and
numerical studies (Clarke et al., 2001; Gorti et al., 2015; Kunitomo et al.,
2020, 2021), but none of them consider accurate radial profiles of
photoevaporation that are derived from detailed radiation hydrodynamics
calculation with non-equilibrium chemistry.

In the present paper, we perform a set of long-term one-dimensional
simulations of PPD evolution considering realistic photoevaporation profiles.
Our calculations incorporate, for the first time, physical models of viscous
accretion, MHD winds, and photoevaporation for PPDs around a wide range of the
central stellar mass. We follow the disk evolution until the disk is dispersed
nearly completely, and derive the disk lifetime accurately. We also study in
detail the mass-loss processes at a variety of evolutionary phases.

The rest of the paper is organized as follows. In Section 2, we explain the
methods we apply. In Section 3, we show the main results. In Section 4, we
discuss detailed properties of our model calculations. Finally in Section 5,
we summarise the paper.

##  2 Numerical methods

We perform long-term disk evolution simulations varying the central stellar
mass in the range of 0.5–7.0${\rm\,M_{\odot}}$. We adopt the cylindrical
coordinates $(r,\phi,z)$ assuming the disk is axisymmetric about the z-axis.
We follow the evolution of the disk surface density, $\Sigma=\int\rho\,{\rm
d}z$, which is the integrated gas density in $z$-direction. We calculate the
time evolution of the gas temperature consistently until the disk disperses.
We incorporate accretion, MHD winds and photoevaporation in the following
manner. The governing equations are

$\begin{split}&\frac{\partial\Sigma}{\partial
t}+\frac{1}{r}\frac{\partial}{\partial r}\left(r\Sigma
v_{r}\right)+\dot{\Sigma}_{\textrm{w}}+\dot{\Sigma}_{\textrm{pe}}=0,\\\
&r\Sigma v_{r}=-\frac{2}{r\Omega}\left[\frac{\partial}{\partial
r}\left(r^{2}\Sigma\overline{\alpha_{r\phi}}c_{\rm
s}^{2}\right)+r^{2}\overline{\alpha_{\phi z}}\left(\rho c_{\rm
s}^{2}\right)_{\textrm{mid}}\right],\end{split}$

where $\dot{\Sigma}_{\textrm{w}}$ and $\dot{\Sigma}_{\textrm{pe}}$ are the
surface mass-loss rates by MHD winds and photoevaporation, respectively.
Details of the mass-loss profiles by MHD winds and by photoevaporation are
described later in this section. We incorporate the two mass-loss mechanisms
as a simple sum of $\dot{\Sigma}_{\textrm{w}}$ and
$\dot{\Sigma}_{\textrm{pe}}$. Although this might overestimate the total mass-
loss rates, we have confirmed that the disk evolution does not vary
significantly if we use a conservative setup with taking the larger value of
$\dot{\Sigma}_{\textrm{w}}$ and $\dot{\Sigma}_{\textrm{pe}}$, in which the
only dominant mechanism operates at a given time. This is simply because only
one effect dominates for most of the time during the disk evolution. In the
above equations, $v_{r}$, $\Omega$ and $c_{\rm s}$ are the velocity in
$r$-direction, the angular velocity and the sound speed respectively. The
subscript ‘mid’ expresses the value at the mid-plane. The dimensionless
parameters $\alpha_{r\phi}$ and $\alpha_{\phi z}$ express the efficiency of
the viscous and wind-driven accretion, respectively. The mass-weighted
averages are denoted as $\overline{\alpha_{r\phi}}$ and
$\overline{\alpha_{\phi z}}$ (see Suzuki et al. (2016) for the definition.)

We calculate the disk gas temperature by considering irradiation from the
central star and viscous heating following Suzuki et al. (2016). The
$z$-averaged temperature is calculated from the balance equation

$T^{4}=T_{\textrm{irr}}^{4}+T_{\textrm{vis}}^{4},$

where $T_{\textrm{irr}}$ and $T_{\textrm{vis}}$ are the temperatures
corresponding to the energy by irradiation and by viscous heating. The
irradiation temperature is given by

$T_{\textrm{irr}}=280{\rm\,K}\left(\frac{L_{*}}{{\rm\,L_{\odot}}}\right)^{1/4}\left(\frac{r}{1{\rm\,au}}\right)^{-1/2},$

where $L_{*}$ is the stellar bolometric luminosity and ${\rm\,L_{\odot}}$ is
the solar bolometric luminosity. Here we assume that the temperature reaches
equilibrium quickly so that the gas temperature is equal to the dust
temperature. The above equation effectively approximates that the dust
temperature is determined by the incident radiation from the central star. The
bolometric luminosity of the central star at the chracteristic age of $1{\,\rm
Myr}$ is adopted from Gorti & Hollenbach (2009). There are several stellar
evolution models for pre-main sequence stars, but there remains a discrepancy
especially for stars younger than $1{\,\rm Myr}$ (Tout et al., 1999; Siess et
al., 2000). For completeness of our model, we use the stellar luminosity at
the age of $1{\,\rm Myr}$ (Table 1).

The viscous heating temperature is determined by the equation

$2\sigma_{\textrm{SB}}T_{\textrm{vis}}^{4}=\left(\frac{3}{8}\tau_{\textrm{R}}+\frac{1}{2\tau_{\textrm{P}}}\right)F_{\textrm{rad}},$

where $\sigma_{\textrm{SB}},\tau_{\textrm{R}},\tau_{\textrm{P}}$, and
$F_{\textrm{rad}}$ are Stefan-Boltzman constant, the Rosseland mean optical
depth, the Planck mean optical depth, and the radiation flux from the disk
surface (Nakamoto & Nakagawa, 1994). This equation expresses the energy
transfer in $z$-direction under the assumption of thermodynamical equilibrium
in a geometrically thin disk. The disk gas is heated by viscous dissipation,
and the deposited thermal energy diffuses out to the disk surface. A part of
the energy liberated by accretion is finally released from the disk surface as
radiation. The Rosseland mean optical depth is expressed as

$\tau_{\textrm{R}}=\kappa_{\textrm{R}}\Sigma,$

where $\kappa_{\textrm{R}}$ is the opacity of dust grains and given by Hueso &
Guillot (2005) as

$\kappa_{\textrm{R}}=\left\\{\begin{array}[]{ll}4.5\left(\frac{T}{150{\rm\,K}}\right)^{2}{\rm\,cm}^{2}{\rm\,g}^{-1}&(T<150{\rm\,K})\\\
4.5{\rm\,cm}^{2}{\rm\,g}^{-1}&(150{\rm\,K}\leq T\leq 1500{\rm\,K})\\\
0{\rm\,cm}^{2}{\rm\,g}^{-1}&(T>1500{\rm\,K})\end{array},\right.$

assuming silicate and water ice as the dominant opacity agents (Pollack et
al., 1985). Note that the temperature dependence changes at around
$150{\rm\,K}$ because of water ice depletion. In our calculation, we adopt a
smooth function

$\begin{split}\kappa_{\textrm{R}}=2.25{\rm\,cm}^{2}{\rm\,g}^{-1}\,&\left[1-\tanh\left(\frac{T-1500{\rm\,K}}{150{\rm\,K}}\right)\right]\\\
&\times\textrm{min}\left[1,\left(\frac{T}{150{\rm\,K}}\right)^{2}\right]\end{split}$

which approximates the above values well. For low-temperature dust, the
Rosseland and Planck mean opacities can be expressed by a power-law, and the
relation $\tau_{\textrm{P}}=2.4\tau_{\textrm{R}}$ holds at low temperature
(Nakamoto & Nakagawa, 1994). We thus set the Planck mean optical depth as

$\tau_{\textrm{P}}=\textrm{max}(2.4\tau_{\textrm{R}},0.5),$

where a lower limit of 0.5 is adopted to reproduce the optically thin limit.

We consider a wide range of the effective viscosity $\alpha$ as a parameter.
Recent ALMA observations toward star-forming regions suggest that the $\alpha$
parameter varies over two orders of magnitude (Rafikov, 2017; Ansdell et al.,
2018). Ansdell et al. (2018) estimated the gas radius of each disk from
${}^{12}\ce{CO}$ line emission. They simulated the evolution of gas radius
with varying $\alpha$ parameter in the range of $10^{-4}$–$10^{-2}$ assuming
an initial disk model by Facchini et al. (2017) and compared the result with
the observations. They found a wide distribution of the value of $\alpha$ over
a few orders of magnitudes. Several observations have been conducted to
estimate the disk accretion rate using the H$\alpha$ equivalent width (Fedele
et al., 2010; Mathews et al., 2012; Hartmann et al., 2016). These studies also
have shown that there is a large variation in the accretion rate. The optical-
UV observations toward T Tauri stars have shown that there is a positive
relationship between the accretion rate and the stellar mass given as
$\dot{M}_{\textrm{acc}}\propto M_{*}^{2}$ (Muzerolle et al., 2003; Calvet et
al., 2004; Hartmann et al., 2006). We thus assume

$\overline{\alpha_{r\phi}}\propto M_{*}$

to match the observed trend. In general, viscous accretion is effective only
in the radial portions of the disk where magnetorotational instability (MRI
hereafter) (Velikhov, 1959; Chandrasekhar, 1961; Balbus & Hawley, 1991)
operates to generate strong turbulence. While MRI actively induces MHD
turbulence in the inner region $<1{\rm\,au}$ and the regions near the disk
surfaces owing to sufficient ionization, an MRI-inactive region called a dead
zone (e.g., Gammie, 1996) is supposed to occupy $1\lesssim r\lesssim$ a few
tens ${\rm\,au}$ of a PPD because of insufficient ionization near the midplane
(Mori & Okuzumi, 2016; Pinte et al., 2016; Flaherty et al., 2017). On the
other hand, a low but finite value of $\alpha_{r\phi}$ can still be sustained
by purely hydrodynamical processes such as vertical shear instability (Urpin &
Brandenburg, 1998; Nelson et al., 2013; Lin & Youdin, 2015; Flock et al.,
2020; Manger et al., 2020). Although the ionization degree is expected to be
varied over a large distance from the central star, the detailed radial extent
of the dead zone is not still well understood quantitatively. Considering
viscous accretion is dominant in the inner region, we assume a constant
$\overline{\alpha_{r\phi}}$ throughout the disk. We adopt
$\overline{\alpha_{r\phi}}=1.0\times 10^{-4}(M_{*}/1{\rm\,M_{\odot}})$ as our
fiducial value. In order to examine the impact of enhanced accretion in MRI-
active disks, we also perform simulations with
$\overline{\alpha_{r\phi}}=1.0\times
10^{-2}(M_{*}/1{\rm\,M_{\odot}}),1.0\times 10^{-3}(M_{*}/1{\rm\,M_{\odot}})$.

We parametrise the wind loss term $\dot{\Sigma}_{\textrm{w}}$, following
Suzuki et al. (2010):

$\dot{\Sigma}_{\textrm{w}}=(\rho c_{\rm
s})_{\textrm{mid}}C_{\textrm{w}}=\frac{\Sigma\Omega}{\sqrt{2\pi}}C_{\textrm{w}},$
(1)

where $C_{\textrm{w}}$ expresses a dimensionless mass flux given by

$C_{\textrm{w}}=\textrm{min}\left(C_{\textrm{w},0},C_{\textrm{w},e}\right).$

Here, the constant maximum value $C_{\textrm{w},0}$ is estimated from local
shearing box MHD simulations of Suzuki & Inutsuka (2009). We adopt
$C_{\textrm{w},0}=2.0\times 10^{-5}$ for the MRI-active case and
$C_{\textrm{w},0}=1.0\times 10^{-5}$ for the MRI-inactive case. We calculate
the energetics-constrained mass flux following Suzuki et al. (2016). We
consider two cases with strong and weak winds. In the strong wind case, all
the liberated gravitational energy is used to drive MHD winds. In this case,
the mass flux $C_{\textrm{w},e}$ and the energy flux $F_{\textrm{rad}}$ are
given by

$\begin{split}C_{\textrm{w},e}&=\textrm{max}\left[\frac{2}{r^{3}\Omega(\rho
c_{\rm s})_{\textrm{mid}}}\frac{\partial}{\partial
r}\left(r^{2}\Sigma\overline{\alpha_{r\phi}}c_{\rm s}^{2}\right)+\frac{2c_{\rm
s}}{r\Omega}\overline{\alpha_{\phi z}},0\right]\\\
F_{\textrm{rad}}&=\textrm{max}\left[-\frac{1}{r}\frac{\partial}{\partial
r}\left(r^{2}\Sigma\Omega\overline{\alpha_{r\phi}}c_{\rm
s}^{2}\right),0\right].\end{split}$

In the weak wind case, a relatively small fraction of the sum of the liberated
gravitational energy and the energy by viscous heating is spent to drive MHD
winds, and the rest is emitted as radiation. We define the fractional ratio of
energy used to launch winds as $\epsilon_{\textrm{rad}}$. We then set

$\begin{split}C_{\textrm{w}}&=(1-\epsilon_{\textrm{rad}})\left[\frac{3\sqrt{2\pi}c_{\rm
s}^{2}}{r^{2}\Omega^{2}}\overline{\alpha_{r\phi}}+\frac{2c_{\rm
s}}{r\Omega}\overline{\alpha_{\phi z}}\right]\\\
F_{\textrm{rad}}&=\epsilon_{\textrm{rad}}\left[\frac{3\sqrt{2\pi}(\rho c_{\rm
s}^{3})_{\textrm{mid}}}{2}\overline{\alpha_{r\phi}}+r\Omega\overline{\alpha_{\phi
z}}(\rho c_{\rm s}^{2})_{\textrm{mid}}\right].\end{split}$

For the main results presented in the following, we adopt the strong disk wind
case as our fiducial model. We discuss the effect of weak wind separately in
Section 4, where we set $\epsilon_{\textrm{rad}}=0.9$ for the weak wind case.
Note that this is a limiting case in which the radiation loss is maximally
evaluated.

Local shearing-box simulations show that the MHD wind torque satisfies
$\overline{\alpha_{\phi z}}\sim 10^{-5}$–$10^{-3}$ and has a positive
dependence on the strength of the net vertical magnetic field $B_{z}$ (Bai,
2013). The dependence is well approximated as

$\overline{\alpha_{\phi z}}\propto\left(\frac{B_{z}^{2}}{8\pi(\rho c_{\rm
s}^{2})_{\textrm{mid}}}\right)^{0.66}.$

The variation of $\overline{\alpha_{\phi z}}$ can be described by the ratio of
the surface density with respect to the initial value on the assumption that
the vertical magnetic flux stays constant during the disk evolution. Then
$\overline{\alpha_{\phi z}}$ is expressed as

$\overline{\alpha_{\phi
z}}=\textrm{min}\left[10^{-5}\left(\frac{\Sigma}{\Sigma_{\textrm{int}}}\right)^{-0.66},1\right],$
(2)

where $\Sigma_{\textrm{int}}$ is the initial surface density. Since the
evolution of the magnetic field is still poorly understood, we place a
conservative upper limit on the value of $\overline{\alpha_{\phi z}}$. We
should note that Armitage et al. (2013) reported that the vertical magnetic
field also diffuses away as the surface density decreases; in this case, the
dependence of $\alpha_{\phi,z}$ on $\Sigma$ would be weaker than in Eq. 2.

We adopt the profile of $\dot{\Sigma}_{\textrm{pe}}$ based on the calculations
by Komaki et al. (2021). In practice, we fit the results of the radiation
hydrodynamics simulations by a function that combines a quadratic function and
a function with negative power. The radiation hydrodynamics simulations
incorporate EUV, FUV and X-ray radiation from the central star. Komaki et al.
(2021) show that strong FUV radiation heats the disk gas effectively and
drives rapid mass loss around a high mass star. It is also found that the gas
in the outer disk is heated to generate flared structure.

Previous studies have shown that the mass-loss rate by photoevaporation,
$\dot{M}_{\textrm{pe}}$, does not strongly depend on the initial disk mass by
performing hydrodynamics simulations (Wölfer et al., 2019; Nakatani et al.,
2021). $\dot{M}_{\textrm{pe}}$ is primarily determined by the density at the
disk surface rather than the density at the mid-plane, which is directly
connected to the total mass. For the same reason, we also assume that
$\dot{\Sigma}_{\textrm{pe}}$ is constant throughout the calculation.

The initial surface density is configured following the minimum solar disk
model (Hayashi, 1981), which is defined as

$\Sigma_{\textrm{int}}=\Sigma_{1{\rm\,au}}\left(\frac{r}{1{\rm\,au}}\right)^{-3/2}\exp{\left(-\frac{r}{r_{\textrm{cut}}}\right)},$

where $r_{\textrm{cut}}$ is a cut-off radius. Motivated by the observation
that the disk radius increases proportionally to the stellar mass (Andrews et
al., 2018), we assume $r_{\textrm{cut}}\propto M_{*}$. We list the physical
parameters of our simulations in Table 1. Note that we normalize the disk gas
density (and hence mass) by $\Sigma_{\textrm{1{\rm\,au}}}$, the surface
density at $r=1{\rm\,au}$. We determine the value so that the total disk mass
satisfies $M_{\textrm{disk}}=0.117(M_{*}/1{\rm\,M_{\odot}}){\rm\,M_{\odot}}$
as Hayashi (1981) suggests.

Table 1: Fiducial stellar parameters in the model (adopted from Gorti & Hollenbach (2009)) $M_{*}$ (${\rm\,M_{\odot}}$) |  $L_{\textrm{bol}}$ (${\rm\,L_{\odot}}$) |  $M_{\textrm{disk}}$ (${\rm\,M_{\odot}}$) |  $r_{\textrm{cut}}$ (${\rm\,au}$) | log$L_{\textrm{FUV}}$ (${\rm\,erg}{\rm\,s}^{-1}$) | log$\phi_{\textrm{EUV}}$ (${\rm\,s}^{-1}$) | log$L_{\textrm{X-ray}}$ (${\rm\,erg}{\rm\,s}^{-1}$)
---|---|---|---|---|---|---
0.5 | 0.93 | 0.0585 | 15 | 30.9 | 40.1 | 29.8
0.7 | 1.72 | 0.0819 | 21 | 31.3 | 40.5 | 30.2
1.0 | 2.34 | 0.117 | 30 | 31.7 | 40.7 | 30.4
1.7 | 5.00 | 0.199 | 51 | 32.3 | 41.0 | 30.7
3.0 | 14.85 | 0.351 | 90 | 32.9 | 39.0 | 28.7
7.0 | 1687 | 0.819 | 210 | 36.5 | 44.1 | 30.8

We set the computational domain at
$r_{\textrm{in}}=10^{-2}{\rm\,au}<r<r_{\textrm{out}}=10^{4}{\rm\,au}$. We
impose the zero-torque condition on the inner and outer boundaries (Suzuki et
al., 2016; Kunitomo et al., 2020). We calculate the evolution until the disk
mass decreases to an extremely small value of
$M_{\textrm{disk}}<10^{-10}{\rm\,M_{\odot}}$. Since the disk mass decreases
rapidly in the last stage, the calculation time does not change even if a
different threshold mass is adopted, provided that it is
$M_{\textrm{disk}}<10^{-5}{\rm\,M_{\odot}}$.

##  3 Results

###  3.1 Around a solar-mass star

Figure 1:  Snapshots of surface density in simulations with
$M_{*}=1,3,7{\rm\,M_{\odot}}$. The initial surface density is shown in dark
blue and as time goes on, the surface density is shown in more yellowish line.

Figure 1 shows the snapshots from 0 to $8{\,\rm Myr}$ for the run with
$M_{*}=1,3,7{\rm\,M_{\odot}}$. In this section, we focus on the
$M_{*}=1{\rm\,M_{\odot}}$ case and we explain the other cases in Section 3.2.
The inner disk in the vicinity of the central star loses quickly a significant
amount of mass by MHD winds in a few million years. After $2{\,\rm Myr}$, the
mass-loss rate by MHD winds, $\dot{M}_{\textrm{w}}$, decreases, and a steep
density profile is found near the central star. This is partly a numerical
effect caused by the upper limit of $\overline{\alpha_{\phi z}}$ defined in
Eq. 2; the mass-loss rate by MHD winds is effectively limited. To examine the
numerical effect, we perform a set of simulations varying the upper limit of
$\overline{\alpha_{\phi z}}$ by a factor of 0.1 and 10. We have confirmed that
the mass-loss rate and the disk lifetime are not affected by this change.

Figure 2:  The surface mass-loss profiles by the two disk dispersal processes;
MHD winds(orange), photoevaporation(green) in the case of
$M_{*}=1{\rm\,M_{\odot}}$. The triangles express the disk radius at the age of
$0,1,3,5{\,\rm Myr}$.

Figure 2 shows $\dot{\Sigma}_{\textrm{w}}$ and $\dot{\Sigma}_{\textrm{pe}}$ at
the age of $t=1,3$ and $5{\,\rm Myr}$. We also indicate by triangles the disk
radii measured from our simulation outputs. Photoevaporation is the main
dispersal process in the outer disk with $r\gtrsim 20{\rm\,au}$ from the
initial stage, while the inner disk is mainly dispersed by MHD winds. As shown
in Figure 2, $\dot{\Sigma}_{\textrm{w}}$ decreases with time because it is
proportional to the surface density as given in Eq. 1. These features are
consistent with the results of Kunitomo et al. (2020).

An important result is that the disk radius decreases with time as indicated
by the triangles in Figure 2. Photoevaporation is the dominant dispersal
process in the outer region, which we expect to cause the disk radius to
decrease. We have run test simulations without photoevaporation to study its
effect clearly. We have found that viscous accretion re-distributes the
angular momentum and can indeed cause bloating of the disk even if
$\overline{\alpha_{r\phi}}$ is set to a very low value of $10^{-4}$. We thus
conclude that photoevaporation plays a crucial role in shaping the disk
morphology and its size through the effective mass-loss in the outer region.

At the age of $\sim 5{\,\rm Myr}$, the major mass-loss is caused by
photoevaporation in $r>1{\rm\,au}$. We explain the details of
$\dot{\Sigma}_{\textrm{pe}}$ and investigate the impact of distinct
photoevaporation profile on disk evolution later in Section 4.3.

Figure 3:  The mass-loss rates by the three dispersal processes; accretion
(blue), MHD winds (orange), and photoevaporation (green) in our fiducial run
with $M_{*}=1{\rm\,M_{\odot}}$. The instantaneous mass-loss rates are
calculated using Eq. LABEL:eq:accretionprocess and Eq. 4.

Figure 3 shows the mass-loss rate as a function of time. The mass-loss rate by
accretion is calculated as

$\begin{split}&\dot{M}_{\textrm{acc}}=-2\pi(rv_{r}\Sigma),\\\ \end{split}$ (3)

using the values at the inner boundary. The mass-loss rates by MHD winds and
photoevaporation are calculated as follows.

$\begin{split}&\dot{M}_{\textrm{w}}=\int_{r_{\textrm{in}}}^{r_{\textrm{out}}^{\prime}}\dot{\Sigma}_{\textrm{w}}\,{\rm
d}r\\\
&\dot{M}_{\textrm{pe}}=\int_{r_{\textrm{in}}}^{r_{\textrm{out}}^{\prime}}\dot{\Sigma}_{\textrm{pe}}\,{\rm
d}r,\\\ \end{split}$ (4)

where $r_{\textrm{out}}^{\prime}$ is the disk radius. Initially, strong MHD
winds disperse a large amount of mass, and photoevaporation becomes the
dominant dispersal process after a few Myrs, as has been discussed in the
above. In contrast, $\dot{M}_{\textrm{pe}}$ is almost constant with time. The
slight decrease in the photoevaporative mass-loss (Figure 3) is explained by
the decreasing outer radius of the disk.

We integrate the mass-loss rate over time to derive the total mass lost. We
calculate the fractional contributions from the three dispersal processes. The
fraction of accretion, MHD winds and photoevaporation are $24\%$, $62\%$ and
$13\%$, respectively. Clearly, more than a half of the disk mass is dispersed
by MHD winds. Since $\dot{\Sigma}_{\textrm{w}}$ scales proportional to the
surface density, the winds are blown with high densities from near the central
star in the initial phase. We note that previous studies suggest that
$\dot{M}_{\textrm{w}}$ may often be overestimated. The 3D MHD simulations of
an accretion disk show that a part of the gas launched as MHD winds falls back
onto the central star through funnel-walls (Takasao et al., 2018, 2020). In
practice, the partially re-accreted gas should be included in
$\dot{M}_{\textrm{acc}}$ rather than in $\dot{M}_{\textrm{w}}$. For a similar
reason, the upper limit of $C_{\textrm{w}}$ should be appropriately calibrated
in our future work.

Protostellar jets can be another mass dispersal process. The mass-loss rate
due to jets is estimated to be about $10\%$ of the accreted mass around low-
mass stars (Calvet et al., 2004; Natta et al., 2014). Since accretion itself
is a relatively minor mass-loss process according to our calculations, we
expect that the influence of jets on the disk dispersal would be limited.

###  3.2 Simulations with different stellar masses

We perform disk evolution simulations with varying the stellar mass in the
range of 0.5–7${\rm\,M_{\odot}}$. We compute the mass-loss rate by each
process in the same manner as in the 1${\rm\,M_{\odot}}$ case. In all cases,
the disk is dispersed dominantly by MHD winds in the early phase, and then
photoevaporation replaces to become the dominant disk dispersal process in the
last stage.

We find very similar disk evolution in the runs with lower mass stars
($<1{\rm\,M_{\odot}}$). Initially, MHD winds disperse the disk gas, and
photoevaporation becomes the main dispersal process at $t\sim 4$–$5{\,\rm
Myr}$. Similarly to our $1{\rm\,M_{\odot}}$ run, about $60\%$ of the total
mass-loss is due to MHD winds.

The case with $3{\rm\,M_{\odot}}$ shows a noticeable difference. The evolution
of the surface density is shown in the middle panel of Figure 1. The main disk
dispersal process is MHD winds and later switches to photoevaporation as in
our $1{\rm\,M_{\odot}}$ run, but the transition occurs later at the age of
$\sim 6{\,\rm Myr}$. Around a $3{\rm\,M_{\odot}}$ star,
$\dot{M}_{\textrm{pe}}$ is relatively low, because an intermediate-mass star
does not have a well developed convective layer and the magnetic activity on
the stellar surface is weak and hence generates lower X-ray emission. The low
X-ray luminosity results in the overall low mass-loss rate by photoevaporation
(Komaki et al., 2021).

The run with $M_{*}=7{\rm\,M_{\odot}}$ shows another difference; the disk is
dispersed efficiently by both MHD winds and photoevaporation. While MHD winds
dominate the mass-loss in the very initial phase, the primary mechanism shifts
to photoevaporation already at $t=1{\,\rm Myr}$ because the inner disk is
completely dispersed as shown in Figure 1. The contribution of accretion, MHD
winds and photoevaporation to the total mass loss is $23\%$, $45\%$ and
$32\%$, respectively. The relative contribution from the photoevaporation is
larger than the other lower mass cases.

The disk structure and its evolution can be seen more clearly by performing
the following visualization. We reconstruct two-dimensional density maps using
the surface density obtained in the simulation with $M_{*}=7{\rm\,M_{\odot}}$.
Figure 4 shows the snapshots of the reconstructed disk density maps. The
vertical scale height throughout the disk is calculated by assuming
hydrostatic equilibrium in the $z$-direction. We see clearly that the diffuse
gas in the outer disk disperses quickly in less than 1 Myrs. Photoevaporation
critically affects the overall disk shape even though its contribution to the
total mass loss is small in the early epoch. The inner disk has a small scale
height, and the gas density is high on and near the mid-plane. An inner ”hole”
opens at $\sim 10{\rm\,au}$ at $t\sim 1.4{\,\rm Myr}$ by the combined effect
of accretion and MHD winds, but the disk gas remains at $10{\rm\,au}\lesssim
r\lesssim$ 200 ${\rm\,au}$ where the effect is relatively weak.

Figure 4:  Snapshots of the reconstructed gas density in a simulation with
$M_{*}=7.0{\rm\,M_{\odot}}$ at $t=0,0.8,1.0,1.8{\,\rm Myr}$. The color map
shows the density and the black line expresses the surface density.

###  3.3 stellar mass dependence of the disk lifetime

Figure 5:  The disk mass evolution around a $1{\rm\,M_{\odot}}$ star. The blue
line shows the simulation data and the orange dotted line expresses the fit.

Based on our simulations, we propose an objective way of determining the disk
lifetime. Figure 5 shows the evolution of the disk mass in the case of
$M_{*}=1{\rm\,M_{\odot}}$. We find that the time evolution is accurately
fitted with a simple vertical sigmoid function given by

$\begin{split}\log_{10}M_{\textrm{disk}}&=a\log\left(\frac{1-x}{x}\right)+d\\\
x&=\frac{(t/1{\,\rm Myr})-b}{c}.\end{split}$

We treat $a,b,c,d$ in the equation as fitting parameters. With this functional
form, the characteristic dispersal time is given by $(b+c)$ Myr.

Figure 6:  The disk dispersal timescale of each stellar mass case. The blue
circles are the simulation data. The black points and green horizontal bars
represent inner-disk lifetimes estimated from observational disk fractions in
Bayo et al. (2012) and Ribas et al. (2015), respectively. See the main text
for more detail.

We do the fitting and estimate the disk dispersal times for all the cases we
simulated. In Figure 6 we compare the results of our simulations with
available observational estimates based on the disk fractions (Bayo et al.,
2012; Ribas et al., 2015). We assume that the disk fraction decreases
exponentially, following $\exp(-t/t_{\textrm{dis}})$, where $t_{\textrm{dis}}$
is a disk dispersal timescale. Overall our results show that the disk lifetime
is shorter for higher stellar mass, except in the case with
$3{\rm\,M_{\odot}}$ which we discuss separately in Section 4.2. The trend is
consistent with the observations.

It is worth noting here that observations toward protoplanetary disks are
conducted often in infrared wavelength, and thus thermal emission from dust
grains is observed primarily (Ribas et al., 2015; Andrews et al., 2018). Both
theoretical studies and observations suggest that the gas component and dust
are distributed differently (de Gregorio-Monsalvo et al., 2013; Ansdell et
al., 2018; Toci et al., 2021; Long et al., 2022). Recent high-resolution
observations have made it possible to observe the gas disk at the
$10{\rm\,au}$ scale using a variety of molecular lines (Öberg et al., 2021),
to show clear images of disk morphology around stars with
$1$–$2{\rm\,M_{\odot}}$ in the nearby star-forming regions. We expect that
future observations toward a number of disks will reveal evolution of disk
morphology around stars with different stellar masses. In our future work, we
compare our simulation results directly with the gas observations to study if
a disk loses its mass from outside to inside around low-mass stars and also
whether or a disk around a high-mass star has an inner hole at the later
stage.

##  4 Discussion

###  4.1 disk parameter

We have performed disk evolution simulations incorporating accretion, MHD
winds and photoevaporation with several assumptions and approximations. For
instance, we assume that all the energy liberated by accretion is transferred
to drive MHD winds in our fiducial run. Here we discuss possible variations of
the results due to our choice of various model parameters.

First we study cases when only a part of gravitational energy liberated by
accretion is used for winds. In order to examine the effects quantitatively,
we run a series of simulations assuming $10\%$ of liberated gravitational
energy is converted to launch MHD winds.

Figure 7:  Same as Figure 3 but for the weak wind case.

Figure 7 shows the mass-loss rate by each disk dispersal process. The
contribution by MHD winds is lower than in the fiducial case. At the early
phase, accretion and MHD winds are the main dispersal processes. At the later
stage, photoevaporation becomes the dominant process. In the case of
$M_{*}=7{\rm\,M_{\odot}}$, the disk is dispersed mainly by photoevaporation.
We estimate the disk lifetime in the same way as the fiducial case. In the
case of $1{\rm\,M_{\odot}}$, the disk lifetime is longer than that of the
fiducial case only by $\sim 1.4{\,\rm Myr}$, which corresponds to 1.2 times of
the fiducial case. Considering the weak wind case is an extreme limit of the
reduced mass-loss by MHD winds, we conclude that the mass dependence of the
disk lifetime is not affected by the choice of the strong wind or weak wind
setups.

Recent observations show that the viscous accretion efficiency, $\alpha$
parameter, varies by a few orders of magnitudes (Hartmann et al., 2016;
Ansdell et al., 2018). In order to investigate the evolution of disk with high
viscous accretion, we also run simulations with
$\overline{\alpha}_{r\phi}=10^{-2}(M_{*}/1{\rm\,M_{\odot}})$ and
$10^{-3}(M_{*}/1{\rm\,M_{\odot}})$. Both $\dot{M}_{\textrm{acc}}$ and
$\dot{M}_{\textrm{w}}$ increases accordingly. Figure 8 shows the resulting
disk lifetime with different values of $\overline{\alpha}_{r\phi}$. The
lifetime is reduced by 1/5 with an increase in $\overline{\alpha_{r\phi}}$ by
a factor of 100. Interestingly the stellar mass dependence remains essentially
the same. As $\overline{\alpha_{r\phi}}$ increases, both
$\dot{M}_{\textrm{acc}}$ and $\dot{M}_{\textrm{w}}$ become larger. At the same
time, $\dot{M}_{\textrm{pe}}$ also increases because the disk radius increases
by the strong angular momentum transport caused by the high viscosity. These
effects roughly cancell out, and the overall disk dispersal time is not
significantly affected by the choice of $\overline{\alpha_{r\phi}}$.

Figure 8:  The disk dispersal timescale in simulations varying
$\overline{\alpha_{r\phi}}$ parameter. The blue, brown and red circles are the
simulation data with
$\overline{\alpha_{r\phi}}=10^{-4}(M_{*}/1{\rm\,M_{\odot}}),10^{-3}(M_{*}/1{\rm\,M_{\odot}})$
and $10^{-2}(M_{*}/1{\rm\,M_{\odot}})$.

###  4.2 star and disk evolution

Kunitomo et al. (2021) performed long-term simulations of disk evolution
considering accretion, photoevaporation and stellar evolution. They adopted
the photoevaporation mass-loss rate depending on the time varying stellar
luminosity in order to study the effect of stellar evolution on disk
evolution. In particular, the X-ray luminosity of a $M_{*}=3{\rm\,M_{\odot}}$
star decreases by a factor of $\sim 3$ to $t\sim 1.5{\,\rm Myr}$. They
clarified that the reduction of X-ray luminosity resulted in the decrease of
the mass-loss rate. They also showed that a disk around a $\sim
3{\rm\,M_{\odot}}$ star had a long lifetime because of the low X-ray
luminosity as we have also found in our simulation.

In our calculations, we use $\dot{\Sigma}_{\textrm{pe}}$ derived from the
numerical simulations of Komaki et al. (2021). EUV, FUV and X-ray photons are
considered as heating sources, and it is shown that FUV radiation efficiently
heats the gas at the disk surface. According to Kunitomo et al. (2021), FUV
luminosity increases by at least a few orders of magnitudes for a star with
intermediate mass of 1.5–4${\rm\,M_{\odot}}$. The FUV luminosity increases
because the stellar surface becomes hot enough to radiate FUV photons
directly. In the lower-mass side of $M_{*}<2.5{\rm\,M_{\odot}}$, the
luminosity increases after the age of several ${\,\rm Myr}$, which is
comparable to the disk dispersal timescale. On the other hand, FUV luminosity
increases earlier at the age of $\sim 1{\,\rm Myr}$ around a
$3{\rm\,M_{\odot}}$ star. This is possibly early enough to affect the disk
evolution by photoevaporation. Komaki et al. (2021) also run a series of
photoevaporation simulations varying the stellar luminosity. In all the cases
of $M_{*}=0.5,1,3{\rm\,M_{\odot}}$, the mass-loss rates follow the
relationship, $\dot{M}_{\textrm{pe}}\propto L_{\textrm{FUV}}^{0.5}$. In the
light of this, we perform a disk evolution calculation around a
$3{\rm\,M_{\odot}}$ star by incorporating the time-dependent FUV luminosity.
We approximate the evolution of FUV luminosity given by Kunitomo et al. (2021)
with a sigmoid function. We obtain the surface mass-loss profile by
multiplying the same factor as $\dot{M}_{\textrm{pe}}$ from the fiducial
value. The disk around a $3{\rm\,M_{\odot}}$ star has a lifetime of $\sim
5.1{\,\rm Myr}$, which is half of the fiducial value. Photoevaporation plays a
crucial role in determining the disk lifetime. It is important to consider
evolution of the stellar luminosity around an intermediate-mass star.

Kunitomo et al. (2021) have shown that the FUV luminosity of a
$1{\rm\,M_{\odot}}$ star decreases gradually by a factor of 10 from $t=1{\,\rm
Myr}$ to $3{\,\rm Myr}$ when the main source of the FUV radiation changes from
accretion to radiation from the stellar chromosphere because of the rapid
decrease of the accretion rate. Since $\dot{M}_{\textrm{pe}}$ decreases as the
FUV luminosity decreases (Komaki et al., 2021), our calculation may
overestimate the effect of photoevaporation, and thus may underestimate the
disk lifetime. We expect that the stellar dependence of disk lifetime is
steeper if we consider evolution of stellar luminosity.

###  4.3 photoevaporation model

We constructed $\dot{\Sigma}_{\textrm{pe}}$ based on Komaki et al. (2021).
While previous theoretical studies incorporated EUV and X-ray radiation as a
heating source, Komaki et al. (2021) also considered FUV radiation. EUV
photons are absorbed by the gas near the central star and ionize hydrogen
atoms. As a result, Hii regions are formed in the vicinity of the central
star, and EUV photons do not reach the outer disk. FUV photons contribute to
the disk heating by photoelectric heating on dust grains. X-ray photons
penetrate into the deeper region of a disk and launch a dense gas flow.
Previous disk evolution simulations take the maximum value among the mass-loss
rates by each radiation and incorporate as $\dot{\Sigma}_{\textrm{pe}}$. The
photoevaporation simulations by Nakatani et al. (2018a, b) and Komaki et al.
(2021) showed that X-ray photons promoted photoelectric heating by ionization
of various species and the decrease of the positive charges of dust grains.
This suggests that it is necessary to perform disk photoevaporation
simulations considering FUV and X-ray radiation at the same time.

Previous calculations of disk evolution (e.g., Kunitomo et al., 2020, 2021,
and more if any) often use $\dot{\Sigma}_{\textrm{pe}}$ given by Owen et al.
(2010). Figure 9 shows the difference in photoevaporation mass-loss profile by
Owen et al. (2010) and Komaki et al. (2021) in the case of
$M_{*}=1{\rm\,M_{\odot}}$. We multiply by $3.1\times 10^{-10}$ the profile of
Owen et al. (2010) by $3.1\times 10^{-10}$ which is given in an arbitrary
unit. We determine the coefficient so that the mass-loss rate inside
$100{\rm\,au}$ becomes equal to that of our photoevaporation profile.

The profile of Komaki et al. (2021) has two peaks generated by EUV and FUV
heating, while the profile by Owen et al. (2010) has one peak. However, this
interesting difference at the inner region of several au does not affect
significantly the results of the simulations because MHD winds are the
dominant disk dispersal process in the vicinity of the central star.

The computational region of Komaki et al. (2021) is
$[0.89{\rm\,au},178{\rm\,au}]$, while that of Owen et al. (2010) is
$[0.82{\rm\,au},100{\rm\,au}]$. In the present paper, we extrapolate the mass-
loss profile for the outer disk. Since photoevaporation becomes the dominant
dispersal process at $r>20{\rm\,au}$, the high mass-loss profile results in
the fast dispersal in the outer region.

Figure 9:  The photoevaporation mass-loss profiles: fitted based on Komaki et
al. (2021)(blue), Owen et al. (2010)(orange) and a simple
approximation(green). Since the fitted data by Owen et al. (2010) is given in
an arbitrary unit, the mass-loss profile in the figure is multiplied by
$3.1\times 10^{-10}$.

In order to examine the dependence on the mass-loss profile of disk dispersal
process, we conduct a disk evolution simulation with a profile given by

$\dot{\Sigma}_{\textrm{pe,s}}=3.2\times
10^{-12}{\rm\,g}{\rm\,cm}^{-2}{\rm\,s}^{-1}\left(\frac{L_{\rm{X-ray}}}{10^{30}{\rm\,erg}{\rm\,s}^{-1}}\right)\left(\frac{r}{2.5{\rm\,au}}\right)^{-2},$

which is a simple fit for the profile by Owen et al. (2010). This simple
profile is shown by an orange line in Figure 9. Even though Owen et al. (2010)
only provides a profile within $100{\rm\,au}$, we extend the fit to
$r>100{\rm\,au}$ so that we incorporate photoevaporation in the outer disk. We
set the inner boundary to $0.14r_{\rm g}=1.2{\rm\,au}$ following Liffman
(2003), which calculate the gravitational radius theoretically. The dominant
disk dispersal process changes from MHD winds to photoevaporation at the age
of $\sim 2{\,\rm Myr}$.

Figure 10:  Snapshots of the reconstructed gas density in a simulation with
$M_{*}=1{\rm\,M_{\odot}}$ at $t=1,4,7{\,\rm Myr}$. The color map shows the
density and the black line expresses the surface density. The top row shows
the density distribution in the fiducial simulation and the bottom row shows
the density distribution in the simulation with the simple photoevaporation
profile .

Figure 10 shows the reconstructed 2D snapshots of the simulation. As shown in
Figure 9, $\dot{\Sigma}_{\textrm{pe,s}}$ decreases with increasing distance
from the central star following $r^{-2}$, while the profile by Komaki et al.
(2021) maintains $\sim 3.0\times 10^{-14}{\rm\,g}{\rm\,cm}^{-2}{\rm\,s}^{-1}$
in the outer region. The disk at a several tens of ${\rm\,au}$ remains at the
last stage of disk dissipation. Following these, we need to apply realistic
photoevaporation profile to obtain the disk density distribution.

###  4.4 Other effects

We assume $\dot{M}_{\textrm{pe}}$ is constant throughout our calculations. It
has been suggested that photoevaporation can be enhanced when an inner hole
opens (Alexander et al., 2014; Owen et al., 2010; Nakatani et al., 2021). In
our study, only the simulation with $M_{*}=7{\rm\,M_{\odot}}$ shows the
formation of a hole at $\sim 1{\,\rm Myr}$, before the disk disperses. The
time evolution of $\dot{M}_{\textrm{pe}}$, especially when a hole or a gap is
formed, may need to be taken into account to evaluate the disk lifetime more
accurately.

Disk observations have shown that disks in a strong radiation field have
specific shapes with long tails (O’dell et al., 1993; McCaughrean & O’dell,
1996; Winter et al., 2018). Theoretical studies suggested that the dense gas
was eroded by high radiation from the nearby high-mass stars (Richling &
Yorke, 1998; Haworth & Clarke, 2019). Haworth & Clarke (2019) performed 2D
hydrodynamics simulations. They assumed that the disk was exposed to the high
FUV radiation fields and solved radiation hydrodynamics to calculate the mass-
loss rate. As s result, the outer disk is effectively removed because the thin
disk is heated efficiently by the radiation from the nearby stars. The present
study does not consider radiation by nearby stars. We expect that the outer
disk would be dispersed more efficiently under an external radiation field,
and the disk radius would shrink rapidly.

##  5 Summary

Recent PPD observations in star forming regions suggested that the lifetime
decreases with increasing central stellar mass. Theoretical studies proposed
three disk dispersal mechanisms: accretion, MHD winds and photoevaporation. In
order to understand both the typical lifetime of a few million years and the
stellar mass dependence, it is necessary to conduct disk evolution simulations
throughout the disk lifetime considering all the disk dispersal mechanisms. We
performed 1D disk evolution simulations varying the central stellar mass in
the range of 0.5–7${\rm\,M_{\odot}}$. We showed the disk loses its mass mainly
by MHD winds in the early stage, and later by photoevaporation. Especially, in
the case with $M_{*}=7{\rm\,M_{\odot}}$, photoevaporation by the high
radiation contributes more than other cases, and a gap opens at $\sim 1{\,\rm
Myr}$. We find that photoevaporation shapes the disk morphology by clearing
out the gas in the outer disk efficiently. The time evolution of disk mass can
be described by a simple function that is given by transforming the sigmoid
function, which yields an accurate estimate of the disk dispersal time. The
dispersal time around a high-mass star is $\sim 2{\,\rm Myr}$, which is
shorter than the other cases by the factor of $\sim 5$. The trend is
consistent with recent observations. Finally, the mass dependence of the
dispersal timescale does not vary by the choice of strong wind or weak wind or
the choice of viscous parameter.

This research was supported by Forefront Physics and Mathematics Program to
Drive Transformation (FoPM), a World-leading Innovative Graduate Study (WINGS)
Program, the University of Tokyo. NY acknowledges financial support from
MEXT/JSPS KAKENHI 20H05844. T.K.S. is supported in part by Grants-in-Aid for
Scientific Research from the MEXT/JSPS of Japan, 17H01105, 21H00033 and
22H01263 and by Program for Promoting Research on the Supercomputer Fugaku by
the RIKEN Center for Computational Science (Toward a unified view of the
universe: from large-scale structures to planets; grant 20351188-PI J. Makino)
from the MEXT of Japan.

## References

  * Alexander et al. (2014) Alexander, R., Pascucci, I., Andrews, S., Armitage, P., & Cieza, L. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 475, doi: 10.2458/azu_uapress_9780816531240-ch021
  * Allers et al. (2007) Allers, K. N., Jaffe, D. T., Luhman, K. L., et al. 2007, ApJ, 657, 511, doi: 10.1086/510845
  * Andrews et al. (2018) Andrews, S. M., Terrell, M., Tripathi, A., et al. 2018, ApJ, 865, 157, doi: 10.3847/1538-4357/aadd9f
  * Ansdell et al. (2018) Ansdell, M., Williams, J. P., Trapman, L., et al. 2018, ApJ, 859, 21, doi: 10.3847/1538-4357/aab890
  * Arakawa et al. (2021) Arakawa, S., Matsumoto, Y., & Honda, M. 2021, ApJ, 920, 27, doi: 10.3847/1538-4357/ac157e
  * Armitage et al. (2013) Armitage, P. J., Simon, J. B., & Martin, R. G. 2013, ApJ, 778, L14, doi: 10.1088/2041-8205/778/1/L14
  * Bai (2013) Bai, X.-N. 2013, ApJ, 772, 96, doi: 10.1088/0004-637X/772/2/96
  * Bai (2016) —. 2016, ApJ, 821, 80, doi: 10.3847/0004-637X/821/2/80
  * Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214, doi: 10.1086/170270
  * Bayo et al. (2012) Bayo, A., Barrado, D., Huélamo, N., et al. 2012, A&A, 547, A80, doi: 10.1051/0004-6361/201219374
  * Béthune et al. (2017) Béthune, W., Lesur, G., & Ferreira, J. 2017, A&A, 600, A75, doi: 10.1051/0004-6361/201630056
  * Blandford & Payne (1982) Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883
  * Calvet et al. (2005) Calvet, N., Briceño, C., Hernández, J., et al. 2005, AJ, 129, 935, doi: 10.1086/426910
  * Calvet et al. (2004) Calvet, N., Muzerolle, J., Briceño, C., et al. 2004, AJ, 128, 1294, doi: 10.1086/422733
  * Carpenter et al. (2006) Carpenter, J. M., Mamajek, E. E., Hillenbrand, L. A., & Meyer, M. R. 2006, ApJ, 651, L49, doi: 10.1086/509121
  * Chandrasekhar (1961) Chandrasekhar, S. 1961, Hydrodynamic and hydromagnetic stability
  * Clarke et al. (2001) Clarke, C. J., Gendrin, A., & Sotomayor, M. 2001, MNRAS, 328, 485, doi: 10.1046/j.1365-8711.2001.04891.x
  * Dahm & Hillenbrand (2007) Dahm, S. E., & Hillenbrand, L. A. 2007, AJ, 133, 2072, doi: 10.1086/512156
  * de Gregorio-Monsalvo et al. (2013) de Gregorio-Monsalvo, I., Ménard, F., Dent, W., et al. 2013, A&A, 557, A133, doi: 10.1051/0004-6361/201321603
  * Ercolano et al. (2009) Ercolano, B., Clarke, C. J., & Drake, J. J. 2009, ApJ, 699, 1639, doi: 10.1088/0004-637X/699/2/1639
  * Facchini et al. (2017) Facchini, S., Birnstiel, T., Bruderer, S., & van Dishoeck, E. F. 2017, A&A, 605, A16, doi: 10.1051/0004-6361/201630329
  * Fairlamb et al. (2015) Fairlamb, J. R., Oudmaijer, R. D., Mendigutía, I., Ilee, J. D., & van den Ancker, M. E. 2015, MNRAS, 453, 976, doi: 10.1093/mnras/stv1576
  * Fang et al. (2012) Fang, M., van Boekel, R., King, R. R., et al. 2012, A&A, 539, A119, doi: 10.1051/0004-6361/201015914
  * Fedele et al. (2010) Fedele, D., van den Ancker, M. E., Henning, T., Jayawardhana, R., & Oliveira, J. M. 2010, A&A, 510, A72, doi: 10.1051/0004-6361/200912810
  * Flaherty et al. (2017) Flaherty, K. M., Hughes, A. M., Rose, S. C., et al. 2017, ApJ, 843, 150, doi: 10.3847/1538-4357/aa79f9
  * Flock et al. (2020) Flock, M., Turner, N. J., Nelson, R. P., et al. 2020, ApJ, 897, 155, doi: 10.3847/1538-4357/ab9641
  * Fulton et al. (2017) Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109, doi: 10.3847/1538-3881/aa80eb
  * Gammie (1996) Gammie, C. F. 1996, ApJ, 457, 355, doi: 10.1086/176735
  * Gorti & Hollenbach (2009) Gorti, U., & Hollenbach, D. 2009, ApJ, 690, 1539, doi: 10.1088/0004-637X/690/2/1539
  * Gorti et al. (2015) Gorti, U., Hollenbach, D., & Dullemond, C. P. 2015, ApJ, 804, 29, doi: 10.1088/0004-637X/804/1/29
  * Gressel et al. (2020) Gressel, O., Ramsey, J. P., Brinch, C., et al. 2020, ApJ, 896, 126, doi: 10.3847/1538-4357/ab91b7
  * Haisch et al. (2001) Haisch, Karl E., J., Lada, E. A., & Lada, C. J. 2001, ApJ, 553, L153, doi: 10.1086/320685
  * Hartmann et al. (2006) Hartmann, L., D’Alessio, P., Calvet, N., & Muzerolle, J. 2006, ApJ, 648, 484, doi: 10.1086/505788
  * Hartmann et al. (2016) Hartmann, L., Herczeg, G., & Calvet, N. 2016, ARA&A, 54, 135, doi: 10.1146/annurev-astro-081915-023347
  * Hasegawa et al. (2017) Hasegawa, Y., Okuzumi, S., Flock, M., & Turner, N. J. 2017, ApJ, 845, 31, doi: 10.3847/1538-4357/aa7d55
  * Haworth & Clarke (2019) Haworth, T. J., & Clarke, C. J. 2019, MNRAS, 485, 3895, doi: 10.1093/mnras/stz706
  * Hayashi (1981) Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35, doi: 10.1143/PTPS.70.35
  * Herczeg & Hillenbrand (2008) Herczeg, G. J., & Hillenbrand, L. A. 2008, ApJ, 681, 594, doi: 10.1086/586728
  * Hernández et al. (2007) Hernández, J., Hartmann, L., Megeath, T., et al. 2007, ApJ, 662, 1067, doi: 10.1086/513735
  * Hollenbach et al. (1994) Hollenbach, D., Johnstone, D., Lizano, S., & Shu, F. 1994, ApJ, 428, 654, doi: 10.1086/174276
  * Hueso & Guillot (2005) Hueso, R., & Guillot, T. 2005, A&A, 442, 703, doi: 10.1051/0004-6361:20041905
  * Kennedy & Kenyon (2009) Kennedy, G. M., & Kenyon, S. J. 2009, ApJ, 695, 1210, doi: 10.1088/0004-637X/695/2/1210
  * Kimmig et al. (2020) Kimmig, C. N., Dullemond, C. P., & Kley, W. 2020, A&A, 633, A4, doi: 10.1051/0004-6361/201936412
  * Komaki et al. (2021) Komaki, A., Nakatani, R., & Yoshida, N. 2021, ApJ, 910, 51, doi: 10.3847/1538-4357/abe2af
  * Kunitomo et al. (2021) Kunitomo, M., Ida, S., Takeuchi, T., et al. 2021, ApJ, 909, 109, doi: 10.3847/1538-4357/abdb2a
  * Kunitomo et al. (2020) Kunitomo, M., Suzuki, T. K., & Inutsuka, S.-i. 2020, MNRAS, 492, 3849, doi: 10.1093/mnras/staa087
  * Lada et al. (2006) Lada, C. J., Muench, A. A., Luhman, K. L., et al. 2006, AJ, 131, 1574, doi: 10.1086/499808
  * Liffman (2003) Liffman, K. 2003, PASA, 20, 337, doi: 10.1071/AS03019
  * Lin & Youdin (2015) Lin, M.-K., & Youdin, A. N. 2015, ApJ, 811, 17, doi: 10.1088/0004-637X/811/1/17
  * Long et al. (2022) Long, F., Andrews, S. M., Rosotti, G., et al. 2022, ApJ, 931, 6, doi: 10.3847/1538-4357/ac634e
  * Lynden-Bell & Pringle (1974) Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603, doi: 10.1093/mnras/168.3.603
  * Mamajek (2009) Mamajek, E. E. 2009, in American Institute of Physics Conference Series, Vol. 1158, American Institute of Physics Conference Series, ed. T. Usuda, M. Tamura, & M. Ishii, 3–10, doi: 10.1063/1.3215910
  * Manger et al. (2020) Manger, N., Klahr, H., Kley, W., & Flock, M. 2020, MNRAS, 499, 1841, doi: 10.1093/mnras/staa2943
  * Mathews et al. (2012) Mathews, G. S., Williams, J. P., Ménard, F., et al. 2012, ApJ, 745, 23, doi: 10.1088/0004-637X/745/1/23
  * McCaughrean & O’dell (1996) McCaughrean, M. J., & O’dell, C. R. 1996, AJ, 111, 1977, doi: 10.1086/117934
  * Meyer et al. (2007) Meyer, M. R., Backman, D. E., Weinberger, A. J., & Wyatt, M. C. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 573\.  https://arxiv.org/abs/astro-ph/0606399
  * Mori & Okuzumi (2016) Mori, S., & Okuzumi, S. 2016, ApJ, 817, 52, doi: 10.3847/0004-637X/817/1/52
  * Muzerolle et al. (2003) Muzerolle, J., Hillenbrand, L., Calvet, N., Briceño, C., & Hartmann, L. 2003, ApJ, 592, 266, doi: 10.1086/375704
  * Nakamoto & Nakagawa (1994) Nakamoto, T., & Nakagawa, Y. 1994, ApJ, 421, 640, doi: 10.1086/173678
  * Nakatani et al. (2018a) Nakatani, R., Hosokawa, T., Yoshida, N., Nomura, H., & Kuiper, R. 2018a, ApJ, 857, 57, doi: 10.3847/1538-4357/aab70b
  * Nakatani et al. (2018b) —. 2018b, ApJ, 865, 75, doi: 10.3847/1538-4357/aad9fd
  * Nakatani et al. (2021) Nakatani, R., Kobayashi, H., Kuiper, R., Nomura, H., & Aikawa, Y. 2021, ApJ, 915, 90, doi: 10.3847/1538-4357/ac0137
  * Natta et al. (2014) Natta, A., Testi, L., Alcalá, J. M., et al. 2014, A&A, 569, A5, doi: 10.1051/0004-6361/201424136
  * Nelson et al. (2013) Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610, doi: 10.1093/mnras/stt1475
  * Öberg et al. (2021) Öberg, K. I., Guzmán, V. V., Walsh, C., et al. 2021, ApJS, 257, 1, doi: 10.3847/1538-4365/ac1432
  * O’dell et al. (1993) O’dell, C. R., Wen, Z., & Hu, X. 1993, ApJ, 410, 696, doi: 10.1086/172786
  * Ogihara et al. (2018) Ogihara, M., Kokubo, E., Suzuki, T. K., & Morbidelli, A. 2018, A&A, 615, A63, doi: 10.1051/0004-6361/201832720
  * Owen et al. (2010) Owen, J. E., Ercolano, B., Clarke, C. J., & Alexand er, R. D. 2010, MNRAS, 401, 1415, doi: 10.1111/j.1365-2966.2009.15771.x
  * Pascucci et al. (2022) Pascucci, I., Cabrit, S., Edwards, S., et al. 2022, arXiv e-prints, arXiv:2203.10068, doi: 10.48550/arXiv.2203.10068
  * Picogna et al. (2021) Picogna, G., Ercolano, B., & Espaillat, C. C. 2021, MNRAS, 508, 3611, doi: 10.1093/mnras/stab2883
  * Picogna et al. (2019) Picogna, G., Ercolano, B., Owen, J. E., & Weber, M. L. 2019, MNRAS, 487, 691, doi: 10.1093/mnras/stz1166
  * Pinte et al. (2016) Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25, doi: 10.3847/0004-637X/816/1/25
  * Pollack et al. (1985) Pollack, J. B., McKay, C. P., & Christofferson, B. M. 1985, Icarus, 64, 471, doi: 10.1016/0019-1035(85)90069-7
  * Rafikov (2017) Rafikov, R. R. 2017, ApJ, 837, 163, doi: 10.3847/1538-4357/aa6249
  * Ribas et al. (2015) Ribas, Á., Bouy, H., & Merín, B. 2015, A&A, 576, A52, doi: 10.1051/0004-6361/201424846
  * Ribas et al. (2014) Ribas, Á., Merín, B., Bouy, H., & Maud, L. T. 2014, A&A, 561, A54, doi: 10.1051/0004-6361/201322597
  * Richling & Yorke (1997) Richling, S., & Yorke, H. W. 1997, A&A, 327, 317
  * Richling & Yorke (1998) —. 1998, A&A, 340, 508
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33
  * Siess et al. (2000) Siess, L., Dufour, E., & Forestini, M. 2000, A&A, 358, 593, doi: 10.48550/arXiv.astro-ph/0003477
  * Suzuki & Inutsuka (2009) Suzuki, T. K., & Inutsuka, S.-i. 2009, ApJ, 691, L49, doi: 10.1088/0004-637X/691/1/L49
  * Suzuki et al. (2010) Suzuki, T. K., Muto, T., & Inutsuka, S.-i. 2010, ApJ, 718, 1289, doi: 10.1088/0004-637X/718/2/1289
  * Suzuki et al. (2016) Suzuki, T. K., Ogihara, M., Morbidelli, A. r., Crida, A., & Guillot, T. 2016, A&A, 596, A74, doi: 10.1051/0004-6361/201628955
  * Tabone et al. (2022) Tabone, B., Rosotti, G. P., Cridland, A. J., Armitage, P. J., & Lodato, G. 2022, MNRAS, 512, 2290, doi: 10.1093/mnras/stab3442
  * Takasao et al. (2018) Takasao, S., Tomida, K., Iwasaki, K., & Suzuki, T. K. 2018, ApJ, 857, 4, doi: 10.3847/1538-4357/aab5b3
  * Takasao et al. (2020) —. 2020, IAU Focus Meeting, 30, 138, doi: 10.1017/S1743921319003892
  * Taki et al. (2021) Taki, T., Kuwabara, K., Kobayashi, H., & Suzuki, T. K. 2021, ApJ, 909, 75, doi: 10.3847/1538-4357/abd79f
  * Tanaka et al. (2013) Tanaka, K. E. I., Nakamoto, T., & Omukai, K. 2013, ApJ, 773, 155, doi: 10.1088/0004-637X/773/2/155
  * Toci et al. (2021) Toci, C., Rosotti, G., Lodato, G., Testi, L., & Trapman, L. 2021, MNRAS, 507, 818, doi: 10.1093/mnras/stab2112
  * Tout et al. (1999) Tout, C. A., Livio, M., & Bonnell, I. A. 1999, MNRAS, 310, 360, doi: 10.1046/j.1365-8711.1999.02987.x
  * Urpin & Brandenburg (1998) Urpin, V., & Brandenburg, A. 1998, MNRAS, 294, 399, doi: 10.1046/j.1365-8711.1998.01118.x
  * Velikhov (1959) Velikhov, E. P. 1959, Zhur. Eksptl’. i Teoret. Fiz., 36.  https://www.osti.gov/biblio/4232891
  * Wang & Goodman (2017) Wang, L., & Goodman, J. 2017, ApJ, 847, 11, doi: 10.3847/1538-4357/aa8726
  * Winter et al. (2018) Winter, A. J., Clarke, C. J., Rosotti, G., et al. 2018, MNRAS, 478, 2700, doi: 10.1093/mnras/sty984
  * Wölfer et al. (2019) Wölfer, L., Picogna, G., Ercolano, B., & van Dishoeck, E. F. 2019, MNRAS, 490, 5596, doi: 10.1093/mnras/stz2939
  * Yasui et al. (2014) Yasui, C., Kobayashi, N., Tokunaga, A. T., & Saito, M. 2014, MNRAS, 442, 2543, doi: 10.1093/mnras/stu1013
  * Zhang et al. (2018) Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJ, 869, L47, doi: 10.3847/2041-8213/aaf744


# An Experimental Study on Private Aggregation of Teacher Ensemble Learning
for End-to-End Speech Recognition

###### Abstract

Differential privacy (DP) is one data protection avenue to safeguard user
information used for training deep models by imposing noisy distortion on
privacy data. Such a noise perturbation often results in a severe performance
degradation in automatic speech recognition (ASR) in order to meet a privacy
budget $\varepsilon$. Private aggregation of teacher ensemble (PATE) utilizes
ensemble probabilities to improve ASR accuracy when dealing with the noise
effects controlled by small values of $\varepsilon$. We extend PATE learning
to work with dynamic patterns, namely speech utterances, and perform a first
experimental demonstration that it prevents acoustic data leakage in ASR
training. We evaluate three end-to-end deep models, including LAS, hybrid
CTC/attention, and RNN transducer, on the open-source LibriSpeech and TIMIT
corpora. PATE learning-enhanced ASR models outperform the benchmark DP-SGD
mechanisms, especially under strict DP budgets, giving relative word error
rate reductions between 26.2% and 27.5% for an RNN transducer model evaluated
with LibriSpeech. We also introduce a DP-preserving ASR solution for
pretraining on public speech corpora.

Index Terms—  privacy-preserving learning, automatic speech recognition,
teacher-student learning, ensemble training.

##  1 Introduction

Automatic speech recognition (ASR) [1] is widely used in spoken language
processing applications, such as smart device control [2], intelligent human-
machine dialogue [3], and spoken language understanding [4]. To build ASR
systems, a large collection of user speech [5] is often required for training
high-performance acoustic models. Protecting information privacy and measuring
numerical privacy loss, such as whether data from a specific user is used for
model training, are becoming critical and prominent research topics for on-
device speech applications [6, 7, 8, 9, 10, 11, 12].

Differential privacy (DP) introduces a noise addition scheme for information
protection characterized by measurable privacy budgets. The noise level is
defined by a privacy budget (e.g., controlled by a minimum $\varepsilon$
value) in $\varepsilon$-DP [13, 14] for a given data set. Machine learning
frameworks based on $\varepsilon$-DP have been shown effective against leakage
of training data (e.g., human faces), model inversion attacks [15] (MIA) and
membership inference, in which a query-free algorithm is used to generate
highly confident test data similar to its training set. Deploying
$\varepsilon$-DP in speech applications based on deep models has recently
attracted much interest [16]. However, directly applying $\varepsilon$-DP
perturbation on the training data could lead to severe performance (e.g.,
prediction accuracy) degradation [17]. Therefore, the noise-enabled protection
process needs to be designed carefully for incorporation into machine learning
for training speech models. Most published DP-based approaches are still
limited to recognition of isolated spoken commands [7]. Designing a continuous
speech recognition system with $\varepsilon$-DP protection needs further
investigation under a variety of privacy settings.

In this work we use $\varepsilon$-DP protection to show that acoustic features
used in ASR training data can been protected against model inversion attack
(MIA) [15]. We also demonstrate that $\varepsilon$-DP can prevent data leakage
from a pretrained ASR model.

Fig. 1: Proposed framework for utilizing private aggregation teacher ensemble
(PATE) to train end-to-end ASR with $\varepsilon$-differential privacy.
Teacher ASR models could also be combined with pretraining on public data (our
second case study).

Private aggregation of teacher ensembles (PATE) [18] learning is a recently
proposed framework that aims to avoid accuracy loss in large-scale visual
prediction models while guaranteeing $\varepsilon$-DP protection. PATE is
guided by a teacher-student learning process, where multiple teachers make up
an ensemble for knowledge transfer. The idea underlying PATE is to benefit
from aggregated noisy outputs of teacher models to compartmentalize
nonsensitive public data with $\varepsilon$-DP protections. PATE and related
approaches [19] were proven effective in mitigating model accuracy
deterioration by employing an ensemble of teacher models with independent
noise perturbations. Nonetheless, for ASR systems dealing with continuous
speech, label stream prediction (e.g., by sequence level modeling) is
required. In this study, we apply the PATE framework to continuous speech
recognition and design different ensemble strategies to ensure
$\varepsilon$-DP for end-to-end ASR models, including RNN transducer [20],
hybrid CTC/attention [21], and LAS (Listener, Attender, Speller) [22]
networks, as shown in the blue squares in Figure 1. Under a strict DP budget
($\varepsilon$=1), the PATE-trained end-to-end models maintain data privacy at
the expense of only a slight increase in word error rate (WER). It also shows
competitive advantages when compared to differentially private stochastic
gradient descent (DP-SGD) [16, 23] benchmarks.

##  2 Related Work

Recent research efforts to preserve data privacy in an ASR system can be
categorized into two groups: (i) systemic, such as federated learning [24],
data isolation [25], and data encryption [26], and (ii) algorithmic, mainly
differentially private machine learning [16]. In this section, we first
summarize some of the privacy-preserving solutions proposed for ASR. Next, we
describe the substratum of differential privacy devised for machine learning
applications and discuss their difference to our proposed approach, while
highlighting its key contributions.

###  2.1 Privacy-Preserving Automatic Speech Recognition

Federated machine learning algorithms [24, 27, 25] have been investigated in
the ASR community to improve information protection. For instance, the average
gradient method [27] was deployed to update the model in ASR training.
Vertical federated learning methods [25] show some other benefits from
isolated features extractors under a heterogeneous computation. However, these
system-level frameworks often rely on assumptions about the constrained
accessibility of the malicious threats and barely provide universal and
statistical justification for privacy guarantees. Cryptographic encryption
[26, 28] and computation protocols [29] are established techniques for
privacy-preserving speaker identification. Meanwhile, these encryption
algorithms and protocols do not consider privacy protection for training
samples, which is a crucial element in developing machine learning models on a
large scale. Lately, differentially-private stochastic gradient descent (DP-
SGD) [16, 23] was introduced to allow quantitative measurements and further
protect identity inference (e.g., of accent or speaking condition) by
introducing an additive distortion during model training. In the remainder of
the paper, a mathematical formulation of differential privacy is provided, as
well as a study of its effectiveness for ASR.

###  2.2 Differential Privacy Fundamentals

The differential privacy mechanism [13] is a standard method to deploy
algorithms with a target privacy guarantee.

Definition 1. A randomized algorithm $\mathcal{M}$ with domain $\mathcal{D}$
and range $\mathcal{R}$ is $(\varepsilon,\delta)$-differentially private if
for any two neighboring inputs (e.g., speech data)
$d,d^{\prime}\in\mathcal{D}$ and for any subset of output predictions (e.g.,
labels) $S\subseteq\mathcal{R}$, the following holds:

$\operatorname{Pr}[\mathcal{M}(d)\in S]\leq
e^{\varepsilon}\operatorname{Pr}\left[\mathcal{M}\left(d^{\prime}\right)\in
S\right]+\delta.$ (1)

The definition above produces a notion of privacy that can be expressed as a
measure of the probabilistic difference of a specific outcome by a
multiplicative factor, $e^{\varepsilon}$, and an additive amount, $\delta$.
The $DP$ mechanism with post-processing [16] (e.g., batch-wise training) is
under a general Renyi-divergence [14] measurement with order
$\alpha\in(1,\infty)$, called $\text{RDP}_{\alpha}$, for all neighboring data
sets $d,d^{\prime}$:

$\begin{array}[]{l}\text{RDP}_{\alpha}\left(\mathcal{M}(d)\|\mathcal{M}\left(d^{\prime}\right)\right)\\\
=\frac{1}{\alpha-1}\log
S_{\theta\sim\mathcal{M}\left(d^{\prime}\right)}\left[\left(\frac{p_{\mathcal{M}(d)}(\theta)}{p_{\mathcal{M}\left(d^{\prime}\right)}(\theta)}\right)^{\alpha}\right]\leq\varepsilon\end{array}$
(2)

As $\alpha\rightarrow\infty$, $\text{RDP}_{\alpha}$ converges to the standard
$(\varepsilon,0)$-DP. Both $\varepsilon$ and $\delta$ should be non-negative.
Considering $\delta\rightarrow$ 0 with only minor relaxation, a smaller value
of $\varepsilon$ indicates a stronger $(\varepsilon,0)$-differentially private
guarantee. In other words, nearly equal probabilities in Eq. (1) would be
given from the neighboring inputs $d$ and $d^{\prime}$, which makes data
identity much hard to infer. Moreover, learning from post-processing features
(e.g., Mel-spectrum [30, 31]) based on the speech data could also be
differentially private, which has been shown by the theorem given in [13, 16].

##  3 PATE Learning for ASR

###  3.1 Noisy Teacher Ensembles for Acoustic Modeling

We now describe PATE [18, 19] based acoustic modeling to enable
$\varepsilon$-DP. First, an ensemble of teacher models is built by
partitioning the training dataset into $I$ disjoint subsets:
$\mathcal{D}_{1},...,\mathcal{D}_{I}$. Next, these are used to train $I$
teacher models independently: $\mathcal{T}_{1},...,\mathcal{T}_{I}$. Each such
model is employed to generate frame-level acoustic model scores, which are
then aggregated by weighted average and used as a teacher for student model
training [32, 33, 34]. For each frame of audio $x$, the final aggregated
teacher model produces a vector of posteriors
$\mathcal{T}_{\mathrm{ens}}(s\mid x)$ over context-dependent states $s$
computed as follows:

$\mathcal{T}_{\mathrm{ens}}(s\mid x)=\sum_{i=1}^{I}w_{i}\mathcal{T}_{i}(s\mid
x),$ (3)

where $\mathcal{T}_{i}(s\mid x)$ is the posterior from the $i$th model, and
$w_{i}$ is its weighting coefficient. The corresponding states from individual
teacher models have the same dimension $J$ before the output alignment (e.g.,
considering a special silent character before alignment).

To ensure $\varepsilon$-DP under the PATE method, a random perturbation was
introduced into the individual teachers’ predictions ($\mathcal{T}_{i}$). We
revise Eq. (3) to obtain a final ensemble from noisy teachers:

$\mathcal{T}_{\operatorname{PATE}}(x,{\lambda})=\sum_{i=1}^{I}w_{i}\left(\mathcal{T}_{i}(x)+Y_{j}({\lambda})\right),$
(4)

where $Y_{1},...,Y_{m}$ are i.i.d. Laplacian or Gaussian random variables with
location $0$ and scale $\lambda^{-1}$. $\lambda$ refers to a privacy parameter
that influences $(\epsilon,\delta)$-DP guarantees and for which a bound has
been proven under composition theorems applicable to model aggregation [18,
19, 33].

As shown in Figure 1, the next PATE step is a process of knowledge transfer,
where the noisy ensemble output is used to relabel a nonsensitive data set,
with a total sample number $K$, which in turn is used to train a student
model, $\mathcal{S}$. Both the prediction outputs and the trained student
model’s internal parameters are free from querying requests, which allows the
only privacy cost to be associated with acquiring the training data for the
student model. We evaluate two noisy aggregation processes, Gaussian NoisyMax
(GNMax) and Laplacian NoisyMax (LNMax) presented in previous studies [18, 19].
Under this setup, the student model is $(\varepsilon,10^{-3})$-DP guaranteed
using $\lambda$ from an analysis in [18, 19]. According to Eq. (4), a large
$\lambda$ refers to a smaller $\varepsilon$, providing a stronger privacy
guarantee, but degrades the accuracy of the labels from the noisy maximum
prediction output of the PATE function.

###  3.2 Sequence-level Teacher Distillations

To transfer the label data for knowledge distillation (KD) on a privacy-
preserving student model, $\mathbf{h}_{\mathcal{T}}^{e}$ and
$\mathbf{h}_{\mathcal{S}}^{e}$ denote the hidden vector representations of the
teacher and student encoders respectively, and
$P_{\mathcal{T}}\left(v\mid\mathbf{h}_{\mathcal{T}}^{d}\right)$ and
$P_{\mathcal{S}}\left(v\mid\mathbf{h}_{\mathcal{S}}^{d}\right)$ the posterior
probabilities computed by teacher and student models, respectively, in regards
to the label $v$. $\mathcal{L}_{KD}$ is defined at the sequence-level as:

$\mathcal{L}_{KD}=-\sum_{S}\sum_{n=1}^{N}P_{\mathcal{T}}\left(\hat{\mathbf{y}}_{n}\mid\mathbf{h}_{\mathcal{T}}^{e}\right)\log
p_{\mathcal{S}}\left(\hat{\mathbf{y}}_{n}\mid\mathbf{h}_{\mathcal{S}}^{e}\right)$
(5)

with $\hat{\mathbf{y}}_{n}$ being the $n$th hypothesis from the set of
$N$-best hypothesis derived from beam-search (i.e., beam width $=N$ ) for the
teacher.

###  3.3 Three Evaluated Acoustic Models

Three deep ASR models are adopted to investigate the impact of PATE. We
consider homogeneously distributed settings, where the same network
architecture is used in both teacher ($\mathcal{T}$) and student
($\mathcal{S}$) models in teacher-student learning.

RNN transducer (RNN-T) [20] utilizes a joint network to combine current
acoustic observations from an encoder network and predictions based on
previously recognized tokens. The encoder network is an RNN that converts the
input acoustic features into a hidden representation for each frame of input.
The prediction network generates output from previous nonblank output labels.
The joint network computes output token logits. The blank symbol is treated as
a possible output to account for the length mismatch between input and output
token sequences, as in CTC [35] models. The loss function of RNN-T [20] is
computed as the negative log posterior of the output label sequence given the
input acoustic feature.

Hybrid CTC/attention (CTC/Att) [21] combines a sequence-to-sequence (seq2seq)
attention-based [36] model with a CTC loss [35]. The encoder processes an
input sequence and creates a hidden latent representation of the same length
as the target sequence. The training data set contains target utterances
composed of predefined vocabulary tokens (characters, tokens, or words). The
CTC loss is computed from the predictions obtained from the encoder with
sequence modeling.

LAS [22] consists of 3 components: (i) a listener (encoder), which is similar
to a standard acoustic model, takes a time-frequency representation of the
speech input and uses a group of neural network layers to map the input to
higher-level features; (ii) an attender, which takes encoder features as input
to learn an alignment between the input features and the predicted subword
units, where each subword is typically a grapheme or word piece; and finally,
(iii) speller (decoder), which takes the output of the attention module and
generates a probability for each hypothesized word.

##  4 Experiments

###  4.1 Experimental Setup

We consider two privacy-preserving conditions for ASR modeling with
$\varepsilon$-DP measurement and protection:

  1. 1.

Parts of the training data require $\varepsilon$-DP protection, and

  2. 2.

All training data requires $\varepsilon$-DP protection,

using two open-source speech corpora for training and evaluation: LibriSpeech
[5] and TIMIT [37].

In Condition 1 ($\mathcal{C}$1), we assume both public and nonpublic data are
from the same domain (e.g., same recording process and acoustic condition),
but part of the training speech is associated with information that is
potentially sensitive and needs protection. We select LibriSpeech and split
the “train-clean960” set into 800 hours of “sensitive” data for learning
teachers models, and 160 hours (as the nonsensitive data shown in Fig. 1) for
training a $\varepsilon$-DP preserving student models with open access to end-
users. We use the Librispeech “test-clean” portion to report word error rates
(WERs).

In Condition 2 ($\mathcal{C}$2), we assume that two data sets coming from
different domains can be accessed during acoustic modeling. One data set
(LibriSpeech in our case) is publicly available in its entirety without
privacy concerns, while the other data set (TIMIT in our case) is private and
needs to be protected. We select different ASR models with pretraining on 800
hours of LibriSpeech used in $\mathcal{C}$1, and adopt the models for fine-
tuning on the TIMIT training set (55 hours) as teacher models. The
$\varepsilon$-DP preserving student ASR model is trained with a knowledge
distillations process using the same unlabeled data as $\mathcal{C}$1\. We use
the 5-hour TIMIT validation data to report WERs.

###  4.2 Baseline and PATE Results

Ensemble baseline: Before incorporating the aggregated ASR system for noisy
teacher-student learning, we first evaluate how the size of the teacher
ensemble influences general ASR performance. In previous work on PATE,
researchers have demonstrated that having sufficient training data for each
teacher model is crucial to ensure a high-performance end-to-end model. We
start with the same number of teacher models reported in [18, 19] with $I$=20
in Eq. (3) and evaluate with the various ASR models discussed in Section 3.3.
The experimental results show that the data scale of the nonpublic training
set used for teacher models could affect an optimal ensemble number. As shown
in Fig. 2(a), the selected ASR models have lower WERs ($<10\%$) with $10$ to
$20$ teachers in $\mathcal{C}$1 (training with LibriSpeech 800 hours as
private). For $\mathcal{C}$2 (training with TIMIT 55 hours as private) shown
in Fig. 2(b), the desired number of teachers becomes smaller ($10$ to $12$) to
obtain competitive WER performance ($<14\%$). We fix the number of teachers at
$15$ in $\mathcal{C}$1 and $10$ in $\mathcal{C}$2, for all following PATE-
based experiments.

(a) $\mathcal{C}$1: LibriSpeech

(b) $\mathcal{C}$2: TIMIT

Fig. 2: Baseline performance under ensemble training for ASR.

Additive noise selection: to study how the $\varepsilon$-DP based noise
injection could impact the ASR performances, we compare two additive noise
mechanisms, GNMax (with Gaussian noise) and LNMax (with Laplacian noise), with
an empirical DP budget of $\varepsilon$=10, as has been reported for different
on-device applications [38]. We follow the PATE framework proposed in Section
3 to obtain the final student ASR models with $\varepsilon$-DP protection. As
shown in Fig. 3, both PATE-GNMax (blue) and PATE-LNMax (orange) show absolute
WER degradation of $2.97\%$ to $8.22\%$ ($\mathcal{C}$1) and $2.73\%$ to
$8.23\%$ ($\mathcal{C}$2) when compared to student model training with clean
aggregated teacher models (gray). RNN-T based ASR with GNMax achieves the best
performance under the $\varepsilon$=10 privacy budget, where PATE-GNMax shows
an average of $4.32\pm 0.12$% absolute WER reduction compared with PATE-LNMax
results.

(a) $\mathcal{C}$1: LibriSpeech

(b) $\mathcal{C}$2: TIMIT

Fig. 3: Student ASR model performance with two additive noisy mechanisms,
GNMax and LNMax, under PATE for teacher training.

PATE for ASR with different privacy budgets: We evaluate the proposed PATE-
based algorithms at different noise levels, which represent different target
privacy budgets ($\lambda=\frac{K}{2\varepsilon}$). When compared with the
benchmark $\varepsilon$-DP preserving DP-SGD [16] algorithm, the results in
Figures 4 (LAS) 5 (Hybrid CTC/Att), and 6 (RNN-T) demonstrate that PATE under
the GNMax mechanism can achieve good performance when combined with different
end-to-end architectures. They show $5.45\pm 0.23\%$ and $9.81\pm 0.43\%$
average WER reduction when comparing PATE-LNMax against DP-SGD [16] powered
$\varepsilon$-DP preserving ASR models. They also highlight that the deployed
models suffer major WER increases for privacy budgets $\varepsilon\leq 10$.
The RNN-T based PATE model outperforms both LAS and Hybrid CTC/attention based
ASR models when $\varepsilon\leq 100$. As an extreme case, the ASR model
starts to converge to its clean ensemble when $\varepsilon\rightarrow 1000$ in
both $\mathcal{C}$1 and $\mathcal{C}$2 conditions.

(a) $\mathcal{C}$1: LibriSpeech

(b) $\mathcal{C}$2: TIMIT

Fig. 4: LAS results with PATEs and DP-SGD.

(a) $\mathcal{C}$1: LibriSpeech

(b) $\mathcal{C}$2: TIMIT

Fig. 5: CTC/Attention results with PATEs and DP-SGD.

(a) $\mathcal{C}$1: LibriSpeech

(b) $\mathcal{C}$2: TIMIT

Fig. 6: RNN-T results with PATEs and DP-SGD.

###  4.3 A Case Study for $\varepsilon$-DP Protection Against MIA

As a first step toward incorporating $\varepsilon$-DP into ASR acoustic
modeling, we demonstrate how PATE can help protect user information in the
training set, as used for training the teacher models. The demonstration
considers the major privacy leakage method based on model inversion attack
(MIA) [15], which utilizes the maximum a posteriori principle and constructs
the input features that maximize the likelihood of observed model output
queries (e.g., using an on-device API). We select the RNN-T based ensemble
model trained under $\mathcal{C}$1 for a privacy-preserving case study against
MIA with 10,000 service queries. We use the MIA algorithm to maximize the
likelihood of a target output word of “stop” and observe its reconstructed
high-confidence input audio. As shown in Fig. 7, MIA can generate inverse Mel-
spectrum output as in Fig. 7(b), similar to an original utterance clip shown
in Fig. 7(a), presumably including speech characteristics unique to the
speaker. However, $\varepsilon$-DP shows its effectiveness in preventing this
information recovery attack with $\varepsilon\leq 10$, as shown Fig. 7 (c) and
Fig. 7 (d). This demonstrates how PATE-trained ASR models can successfully
protect user information.

(a) original (no MIA)

(b) MIA without DP ($\varepsilon$ = $\infty$)

(c) MIA with DP ($\varepsilon$ = 10)

(d) MIA with DP ($\varepsilon$ = 1)

Fig. 7: Model inversion attack (MIA) against ASR with a target word of “stop.”
(c) and (d) show effective $\varepsilon$-DP protection.

##  5 Conclusion

We have conducted a first study in privacy-preserving ASR applying the PATE
framework for $\varepsilon$-differential privacy to three popular ASR
architectures. For all the DP budgets tested, PATE greatly outperforms
previously proposed DP-SGD mechanisms as benchmarks, especially under strict
DP budgets ($\varepsilon$=1). Using the best RNN-T model, we also demonstrated
that the proposed $\varepsilon$-DP preserving PATE models can protect against
potential private data leakage from the training set using model inversion
attacks.

## References

  * [1] Lawrence R. Rabiner and Ronald W. Schafer,  Introduction to Digital Speech Processing, vol. 1,  Now Publishers Inc, 2007.
  * [2] Ian McGraw, Rohit Prabhavalkar, Raziel Alvarez, Montse Gonzalez Arenas, Kanishka Rao, David Rybach, Ouais Alsharif, Haşim Sak, Alexander Gruenstein, Françoise Beaufays, et al.,  “Personalized speech recognition on mobile devices,”  in Proc. IEEE ICASSP, 2016, pp. 5955–5959.
  * [3] Andreas Stolcke, Klaus Ries, Noah Coccaro, Elizabeth Shriberg, Rebecca Bates, Daniel Jurafsky, Paul Taylor, Rachel Martin, Carol Van Ess-Dykema, and Marie Meteer,  “Dialogue act modeling for automatic tagging and recognition of conversational speech,”  Computational Linguistics, vol. 26, no. 3, pp. 339–373, 2000.
  * [4] Janet M Baker, Li Deng, James Glass, Sanjeev Khudanpur, Chin-Hui Lee, Nelson Morgan, and Douglas O’Shaughnessy,  “Developments and directions in speech recognition and understanding, Part 1 [DSP Education],”  IEEE Signal Processing Magazine, vol. 26, no. 3, pp. 75–80, 2009\.
  * [5] Vassil Panayotov, Guoguo Chen, Daniel Povey, and Sanjeev Khudanpur,  “LibriSpeech: an ASR corpus based on public domain audio books,”  in Proc. IEEE ICASSP, 2015, pp. 5206–5210.
  * [6] Peng Cheng and Utz Roedig,  “Personal voice assistant security and privacy—a survey,”  Proceedings of the IEEE, vol. 10, no. 4, pp. 476–507, 2022.
  * [7] Chao-Han Huck Yang, Sabato Marco Siniscalchi, and Chin-Hui Lee,  “PATE-AAE: Incorporating adversarial autoencoder into private aggregation of teacher ensembles for spoken command classification,”  in Proc. Interspeech, 2021, pp. 881–885.
  * [8] Xiaodong Cui, Songtao Lu, and Brian Kingsbury,  “Federated acoustic modeling for automatic speech recognition,”  in Proc. IEEE ICASSP, 2021, pp. 6748–6752.
  * [9] Natalia Tomashenko, Xin Wang, Emmanuel Vincent, Jose Patino, Brij Mohan Lal Srivastava, Paul-Gauthier Noé, Andreas Nautsch, Nicholas Evans, Junichi Yamagishi, Benjamin O’Brien, et al.,  “The VoicePrivacy 2020 Challenge: Results and findings,”  Computer Speech & Language, vol. 74, pp. 101362, 2022.
  * [10] Chao-Han Huck Yang, Zeeshan Ahmed, Yile Gu, Joseph Szurley, Roger Ren, Linda Liu, Andreas Stolcke, and Ivan Bulyko,  “Mitigating closed-model adversarial examples with Bayesian neural modeling for enhanced end-to-end speech recognition,”  in Proc. IEEE ICASSP, 2022, pp. 6302–6306.
  * [11] Haibin Wu, Xu Li, Andy T Liu, Zhiyong Wu, Helen Meng, and Hung-Yi Lee,  “Improving the adversarial robustness for speaker verification by self-supervised learning,”  IEEE/ACM Transactions on Audio, Speech, and Language Processing, vol. 30, pp. 202–217, 2021.
  * [12] Yu-Chen Lin, Tsun-An Hsieh, Kuo-Hsuan Hung, Cheng Yu, Harinath Garudadri, Yu Tsao, and Tei-Wei Kuo,  “Speech recovery for real-world self-powered intermittent devices,”  in Proc. IEEE ICASSP, 2022, pp. 26–30.
  * [13] Cynthia Dwork,  “Differential privacy: A survey of results,”  in International Conference on Theory and Applications of Models of Computation. Springer, 2008, pp. 1–19.
  * [14] Ilya Mironov,  “Rényi differential privacy,”  in 2017 IEEE 30th Computer Security Foundations Symposium (CSF). IEEE, 2017, pp. 263–275.
  * [15] Matt Fredrikson, Somesh Jha, and Thomas Ristenpart,  “Model inversion attacks that exploit confidence information and basic countermeasures,”  in Proceedings of the 22nd ACM SIGSAC Conference on Computer and Communications Security, 2015, pp. 1322–1333.
  * [16] Martin Abadi, Andy Chu, Ian Goodfellow, H Brendan McMahan, Ilya Mironov, Kunal Talwar, and Li Zhang,  “Deep learning with differential privacy,”  in Proceedings of the 2016 ACM SIGSAC conference on computer and communications security, 2016, pp. 308–318.
  * [17] Arun Rajkumar and Shivani Agarwal,  “A differentially private stochastic gradient descent algorithm for multiparty classification,”  in Artificial Intelligence and Statistics. PMLR, 2012, pp. 933–941.
  * [18] Nicolas Papernot, Martín Abadi, Ulfar Erlingsson, Ian Goodfellow, and Kunal Talwar,  “Semi-supervised knowledge transfer for deep learning from private training data,”  in Proc. International Conference on Learning Representations, 2017\.
  * [19] Nicolas Papernot, Shuang Song, Ilya Mironov, Ananth Raghunathan, Kunal Talwar, and Ulfar Erlingsson,  “Scalable private learning with PATE,”  in Proc. International Conference on Learning Representations, 2018\.
  * [20] Alex Graves,  “Sequence transduction with recurrent neural networks,”  Representation Learning Workshop, ICML, 2012.
  * [21] Suyoun Kim, Takaaki Hori, and Shinji Watanabe,  “Joint CTC-attention based end-to-end speech recognition using multi-task learning,”  in Proc. IEEE ICASSP, 2017, pp. 4835–4839.
  * [22] William Chan, Navdeep Jaitly, Quoc Le, and Oriol Vinyals,  “Listen, Attend and Spell: A neural network for large vocabulary conversational speech recognition,”  in Proc. IEEE ICASSP, 2016, pp. 4960–4964.
  * [23] Cynthia Dwork, Guy N Rothblum, and Salil Vadhan,  “Boosting and differential privacy,”  in Proc. IEEE 51st Annual Symposium on Foundations of Computer Science, 2010, pp. 51–60.
  * [24] David Leroy, Alice Coucke, Thibaut Lavril, Thibault Gisselbrecht, and Joseph Dureau,  “Federated learning for keyword spotting,”  in Proc. IEEE ICASSP, 2019, pp. 6341–6345.
  * [25] Chao-Han Huck Yang, Jun Qi, Samuel Yen-Chi Chen, Pin-Yu Chen, Sabato Marco Siniscalchi, Xiaoli Ma, and Chin-Hui Lee,  “Decentralizing feature extraction with quantum convolutional neural network for automatic speech recognition,”  in Proc. IEEE ICASSP, 2021, pp. 6523–6527.
  * [26] Cornelius Glackin, Gerard Chollet, Nazim Dugan, Nigel Cannings, Julie Wall, Shahzaib Tahir, Indranil Ghosh Ray, and Muttukrishnan Rajarajan,  “Privacy preserving encrypted phonetic search of speech data,”  in Proc. IEEE ICASSP, 2017, pp. 6414–6418.
  * [27] Dimitrios Dimitriadis, Kenichi Kumatani, Robert Gmyr, Yashesh Gaur, and Sefik Emre Eskimez,  “A federated approach in training acoustic models,”  in Proc. Interspeech, 2020.
  * [28] Ferdinand Brasser, Tommaso Frassetto, Korbinian Riedhammer, Ahmad-Reza Sadeghi, Thomas Schneider, and Christian Weinert,  “VoiceGuard: Secure and private speech processing.,”  in Proc. Interspeech, 2018, pp. 1303–1307.
  * [29] Manas A Pathak and Bhiksha Raj,  “Privacy-preserving speaker verification and identification using Gaussian mixture models,”  IEEE Transactions on Audio, Speech, and Language Processing, vol. 21, no. 2, pp. 397–406, 2012.
  * [30] Steven Davis and Paul Mermelstein,  “Comparison of parametric representations for monosyllabic word recognition in continuously spoken sentences,”  IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 28, no. 4, pp. 357–366, 1980.
  * [31] John S. Bridle and Michael D. Brown,  “An experimental automatic word recognition system,”  JSRU Report, vol. 1003, no. 5, pp. 33, 1974,  Joint Speech Research Unit Ruislip, England.
  * [32] Yevgen Chebotar and Austin Waters,  “Distilling knowledge from ensembles of neural networks for speech recognition.,”  in Proc. Interspeech, 2016, pp. 3439–3443.
  * [33] Yan Gao, Titouan Parcollet, and Nicholas D Lane,  “Distilling knowledge from ensembles of acoustic models for joint CTC-attention end-to-end speech recognition,”  in Proc. IEEE Automatic Speech Recognition and Understanding Workshop, 2021, pp. 138–145.
  * [34] Chao-Han Huck Yang, Linda Liu, Ankur Gandhe, Yile Gu, Anirudh Raju, Denis Filimonov, and Ivan Bulyko,  “Multi-task language modeling for improving speech recognition of rare words,”  in Proc. IEEE Automatic Speech Recognition and Understanding Workshop, 2021, pp. 1087–1093.
  * [35] Alex Graves, Santiago Fernández, Faustino Gomez, and Jürgen Schmidhuber,  “Connectionist temporal classification: labelling unsegmented sequence data with recurrent neural networks,”  in Proc. 23rd International Conference on Machine Learning, 2006, pp. 369–376.
  * [36] Dzmitry Bahdanau, Jan Chorowski, Dmitriy Serdyuk, Philemon Brakel, and Yoshua Bengio,  “End-to-end attention-based large vocabulary speech recognition,”  in Proc. IEEE ICASSP, 2016, pp. 4945–4949.
  * [37] John S Garofolo, Lori F Lamel, William M Fisher, Jonathan G Fiscus, and David S Pallett,  “DARPA TIMIT acoustic-phonetic continuous speech corpus CD-ROM. NIST speech disc 1-1.1,”  vol. 93, pp. 27403, 1993.
  * [38] Wennan Zhu, Peter Kairouz, Brendan McMahan, Haicheng Sun, and Wei Li,  “Federated heavy hitters discovery with differential privacy,”  in International Conference on Artificial Intelligence and Statistics. PMLR, 2020, pp. 3837–3847.