# Hydrodynamical simulations of proton ingestion flashes in Type I X-ray
Bursts

Simon Guichandut  Department of Physics and Trottier Space Institute, McGill
University, 3600 rue University, Montreal, QC, H3A 2T8, Canada  Michael
Zingale  Department of Physics and Astronomy, Stony Brook University, NY
11794-3800, USA  Andrew Cumming  Department of Physics and Trottier Space
Institute, McGill University, 3600 rue University, Montreal, QC, H3A 2T8,
Canada

###### Abstract

We perform the first multidimensional fluid simulations of thermonuclear
helium ignition underneath a hydrogen-rich shell. This situation is relevant
to Type I X-ray bursts on neutron stars which accrete from a hydrogen-rich
companion. Using the low Mach number fluid code MAESTROeX, we investigate the
growth of the convection zone due to nuclear burning, and the evolution of the
chemical abundances in the atmosphere of the star. We also examine the
convective boundary mixing processes which cause the evolution to differ
significantly from previous one-dimensional simulations that rely on mixing-
length theory. We find that the convection zone grows outwards as penetrating
fluid elements cool the overlying radiative layer, rather than directly from
the increasing entropy of the convection zone itself. Simultaneously, these
flows efficiently mix composition, carrying carbon out of, and protons into
the convection zone even before contact with the hydrogen shell. We discuss
the implications of these effects for future modeling of these events and
observations.

X-ray bursts — Neutron stars — Convection — Hydrodynamics

††software: MAESTROeX (Nonaka et al., 2010), AMReX (Zhang et al., 2019),
pynucastro (Smith et al., 2023), NumPy (Harris et al., 2020), SciPy (Virtanen
et al., 2020), Matplotlib (Hunter, 2007), and yt (Turk et al.,
2011).††facilities: NERSC

##  1 Introduction

Type I X-ray bursts are the result of thermonuclear runaways on the surface of
accreting neutron stars in low-mass X-ray binaries (Lewin et al., 1993;
Galloway & Keek, 2021; Strohmayer & Bildsten, 2006). As the most frequent
transient in the high-energy sky with over 7000 bursts catalogued (Galloway et
al., 2020), these “bursts” can help us understand the properties of the
neutron star surface. The challenge is to properly model the physics of
nuclear burning, fluid motions, radiation, and their interaction.

A critical mechanism which controls the evolution of these bursts is
convection. Even though ignition happens only a few tens of meters below the
star’s surface, the density of the envelope is so large that the fluid cannot
cool by radiation only. The burning fuel will therefore generate a convection
zone which will grow towards the surface in a time shorter than the thermal
time. Many studies have described this convection zone and its evolution. Joss
(1977) derived the basic timescales for radiation and convection, and showed
that convection could not reach all the way to the photosphere. Initial time-
dependent calculations made the simplifying assumption that regions unstable
to convection are perfectly adiabatic and instantaneously mix entropy (Joss,
1978; Hanawa & Sugimoto, 1982). Many simulations of bursts were performed
using stellar evolution codes such as KEPLER (Wallace et al., 1982; Woosley et
al., 2004), SHIVA (José et al., 2010), and MESA (Paxton et al., 2011; Yu &
Weinberg, 2018). These codes have in common that they treat convection with
some variation of mixing-length theory (MLT, e.g. Cox & Giuli, 1968; Henyey et
al., 1965), in which it is assumed that convective parcels travel a fixed
fraction of the local pressure scale height, at the velocity required to
transport the heat flux.

Meanwhile, few studies have been dedicated to modeling the multidimensional
hydrodynamic nature of convection within bursts. Numerically, this is a
challenge of timescales. Burning and convection take place over many scale
heights, and one therefore needs to account for the compressibility of the
fluid. This introduces sound waves into the system, which travel on a sound
crossing time of ${\sim}10\,\mu$s, much shorter than the tens of seconds burst
duration, and the computational cost of the simulations easily becomes
prohibitive. But since bursts proceed as subsonic deflagrations (Wallace et
al., 1982), an alternative approach is to model the fluid in the low Mach
number approximation (see Almgren et al., 2006, and references therein for the
origins of this method). This was first applied to X-ray bursts by Lin et al.
(2006) who performed 2D simulations of a pure helium burning layer and growing
convection zone, and confirmed that the Mach number remained small ($\lesssim
15$%) throughout, justifying the use of the low Mach method. Malone et al.
(2011) used the dedicated low Mach number hydrodynamics code MAESTRO (Nonaka
et al., 2010) to model a similar type of burst. Malone et al. (2014) then
modeled a burst from a uniform mixture of hydrogen and helium, a more common
scenario for accreting neutron stars in bursting sources. Both works found
that a very high spatial resolution (few centimeters per zone or even less,
depending on the temperature sensitivity of the nuclear reactions) was needed
to obtain numerical convergence.

Another possibility for the composition of the envelope prior to the outburst
is a layered structure, where a hydrogen-rich shell sits on top of a pure
helium layer. In sources where the neutron star is accreting solar-like
material (primarily H with mass fraction $X_{0}\sim 0.7$, then He, and some
CNO elements with $Z_{\rm CNO}$ of a few percent), at high enough accretion
rates $\dot{M}_{\rm acc}$ such that the temperature $T$ exceeds $8\times
10^{7}$ K, hydrogen will contiunously burn via the hot CNO cycle, and run out
at the depletion depth

$y_{\rm d}=2.7\times 10^{7}\,\rm{g}\;\rm{cm}^{-2}\left(\frac{\dot{M}_{\rm
acc}}{0.01\dot{M}_{\rm Edd}}\right)\left(\frac{0.02}{Z_{\rm
CNO}}\right)\left(\frac{X_{0}}{0.7}\right)\,,$ (1)

as shown by Cumming & Bildsten (2000). Here $\dot{M}_{\rm Edd}$ is the
Eddington accretion rate for a neutron star with a 12 km radius accreting gas
with $X_{0}=0.7$. The column depth
$y(r)\equiv\int_{r}^{\infty}\rho(r^{\prime})dr^{\prime}$ measures the amount
of mass above a radial coordinate $r$.

Recently, in Guichandut & Cumming (2023), we suggested that the lightcurves of
bursts which ignite at $y>y_{d}$ could be used to test the radial extent of
the convection zone. In particular, we showed that the peculiar observation of
a ${\approx}\,0.7\,$s pause in the rise of a burst from SAX J1808.4–3658 (Bult
et al., 2019) was consistent with the ejection of a hydrogen shell. For the
pause to be so short, this shell would have had to be eroded by convection. We
ran one-dimensional calculations using the stellar evolution code MESA (Paxton
et al., 2011) and verified that convection could indeed extend to the required
column depths of $y_{\rm conv}{\sim}10^{4}{-}10^{5}$ g cm${}^{-2}\ll y_{\rm
d}$. However, our results were sensitive to assumptions made about convection
in the code. One reason for this is that MLT is not appropriate for situations
where the thermodynamics of the gas are changing on timescales similar to the
convective turnover time. In these bursts, a convection zone in the helium
layer will grow and eventually “collide” with the hydrogen-rich shell above,
which will introduce new fuel in the form of free protons which quickly burn
and rapidly change the nature of the convection. This is an example of a
“proton-ingestion flash”, which also occur for example in evolved metal poor
stars (see e.g. Herwig et al., 2011, and references therein). The limitations
of MLT have been noted for these events as well, and recent works have instead
focused on fluid simulations (Herwig et al., 2014).

To understand how convection proceeds in X-ray bursts with proton-ingestion
events, we must model the multidimensional fluid. In this paper, we use the
low Mach number hydrodynamics code MAESTROeX (the successor to MAESTRO; Fan et
al., 2019) to evolve the same burst which we previously studied in Guichandut
& Cumming (2023). As in Malone et al. (2011, 2014), we consider a two-
dimensional box in the plane-parallel neutron star atmosphere on the verge of
thermonuclear ignition. We briefly describe the numerical method and how our
initial model is created in Section 2, and present the results of the
simulations in Section 3.

##  2 Methodology

###  2.1 Numerical method

The basic idea of the low Mach number method is to decompose the pressure into
base-state plus perturbative components,

$p=p(\bm{x})=p_{0}(r)+\pi(\bm{x})\,,$ (2)

where $\bm{x}$ is the cartesian position vector and $r$ is the radial
coordinate in the star. The base-state is assumed to be in hydrostatic
equilbrium, i.e.

$\nabla p_{0}=-\rho_{0}g\bm{\hat{r}}\,,$ (3)

where $\rho_{0}(r)$ is the base state density. It can be shown (Almgren et
al., 2006) that if the Mach number $M=|\bm{U}|/c_{s}$ ($c_{s}$ being the local
speed of sound) is small, then the pressure perturbations $\pi/p_{0}$ are of
order $M^{2}$. Thus, we can ignore $\pi$ in linking the pressure to the
density via the equation of state, which effectively decouples the two and
filters out sound waves. For full details on the equations being solved and
the algorithmic implementation, see Fan et al. (2019).

As in previous papers using MAESTRO to model bursts, we use the Helmholtz
equation of state (Timmes & Swesty, 2000), which has pressure contributions
from ideal gas, radiation, and degenerate electrons. Our reaction network is
equivalent to MESA’s “CNO_extras” network (Paxton et al., 2011), with inert
${}^{56}{\rm Fe}$ added. The complete set of reactions linking 20 isotopes
from ${}^{1}{\rm H}$ to ${}^{24}{\rm Mg}$ is assembled using pynucastro (Smith
et al., 2023).

###  2.2 Initial model

As in Malone et al. (2014), we create an initial model based on a 1D stellar
evolution calculation, in this case our previous MESA study of this type of
burst in Guichandut & Cumming (2023). There, we considered a neutron star with
a mass of $1.4M_{\odot}$, radius $12$ km, and surface gravity of $1.29\times
10^{14}$ cm s-2 (ignoring relativistic corrections). We first accreted a large
column of inert iron to form the bottom of the atmosphere, then let it accrete
a solar mixture of fuel (${}^{1}{\rm H}$, ${}^{4}{\rm He}$ and ${}^{12}{\rm
C}$ with mass fractions $X=0.7$, $Y=0.28$ and $Z=0.02$) until ignition. This
MESA model is plotted with solid lines in Figure 1, in terms of its
temperature and composition.

Figure 1: Initial model from MESA (solid lines), and the initial model fed
into the fluid simulations (dot-dashed lines). The top panel shows the
temperature profile as a function of height. The bottom panel shows the
composition profiles of the main chemical species, where the green color
combines all CNO elements. The offset in the position of the solar layer is
due to the temperature kick given to the initial fluid model.

In practice, it is difficult to directly port a MESA model to a fluid
simulation code, because the former’s grid is by mass coordinate while the
latter is spatial, and due to small differences in the EOS of both codes. This
is why we instead opt to construct a toy model which is defined by a small set
of physical parameters.

The full description of our model is presented in Appendix A. Briefly, our
atmosphere is composed of solar fuel at the top, then pure helium and CNO past
the hydrogen depletion depth, then pure iron. At the bottom of the fuel layer
(just above the iron) is an isentropic zone which becomes convectively
unstable once the simulation begins. The temperature profile above is in
radiative equilibrium, and depends on the flux from nuclear burning. We set
the temperature of the iron layer to an artificially small value so as to
create a “buoyancy wall” that suppresses downward penetrating flows and
inhibits mixing of the iron into the fuel above. Lastly, we “kick” the model
by increasing the temperature of the isentropic zone in order to accelerate
the evolution towards the thermonuclear runaway. This model is shown in Figure
1.

Our model is simpler than the MESA model in two main aspects. First, we assume
that the initial CNO species are ${}^{14}{\rm O}$ and ${}^{15}{\rm O}$ only.
This is valid for $y<y_{\rm d}$ where hot CNO burning is taking place and the
rate-limiting step is the oxygen $\beta$-decays. However for $y<y_{\rm d}$,
there are no more protons to sustain the reactions and the oxygens should
decay to ${}^{15}{\rm N}$, ${}^{14}{\rm N}$, then ${}^{13}{\rm C}$ products
(see Figure 2 in Guichandut & Cumming, 2023). We do not expect this
simplification to affect the results because these species only represent a
few percent of the mass of the fuel, and because burning in the helium layer
is dominated by triple-$\alpha$. Second, we start at a model where the helium
is fully un-burned, whereas in reality triple-$\alpha$ burning would have
already converted some of the helium into carbon at the moment of ignition, as
can be seen with the MESA model in Figure 1. This effectively means that more
fuel is available at ignition, leading to a more explosive burst.

As in Malone et al. (2014), we configure a “sponge” region at the top of the
model. Below a density $\rho=\rho_{\rm cutoff}=10^{3}$ g cm-3, we hold the
density constant. Within this region, MAESTROeX solves the anelastic velocity
constraint, which effectively supresses large spurious velocities which would
otherwise make the calculation diverge. The sponge initially starts at
$z\approx 22$ m, but moves upwards, following the density, as the atmosphere
expands during the burst.

##  3 Results

We ran simulations, starting from the initial model presented in the last
section, at resolutions of 12, 6, 3 and 1.5 cm per zone. The grid is square
with a width of 3072 cm. This being much less than the radius of the star, the
atmosphere is effectively plane-parallel, with a constant gravity
$g=1.29\times 10^{14}$ cm2 s-1. The boundary conditions are periodic on the
sides, outflow at the top, and slip wall (vanishing vertical velocity and
tangential velocity gradient) at the bottom of the box. The density and
enthalpy are held constant at the bottom and at the top. We do not include
thermal diffusion as it is much slower than heat transport via convection at
the depths we are considering111Malone et al. (2011) also showed that thermal
diffusion only had minor effects on the evolution of convection, which is what
we are focusing on..

In Section 3.1, we compare the general evolution of the different simulations
and discuss convergence. In Section 3.2 we consider a single simulation to
analyze the evolution of convection and burning throughout the burst. In
Section 3.3, we delve into the details of the mixing near the convective
boundary.

###  3.1 Resolution dependence and convergence

In order to assess convergence, we plot the instantaneous maximum values in
the grid of the temperature and Mach number in Figure 2. Both peaks are
tracers of the evolution of the burst in time. The temperature curve clearly
shows a transition from a slow increase to a full runaway. As we will see
later, this transition is related to the collision of the convection zone with
the overlying hydrogen-rich shell.

The evolution of the burst in time clearly depends on the resolution. At
greater resolutions (smaller grid spacings), the burst evolves faster. A few
factors contribute to this, which were previously observed and discussed in
Malone et al. (2011) and Malone et al. (2014). First, since the nuclear
burning is very sensitive to temperature, under-resolving the peak of the
temperature profile will underestimate the energy-generation. Second, lower
resolution results in spurious large velocities which overshoot the convection
zone, taking heat away from it and inhibiting its growth. In our case, this
second factor is particularily important as convection is needed to bring
additional fuel into the mixture.

Figure 2: Peak values of the temperature (left axis, solid lines) and mach
number (right axis, dashed lines) as a function of time for the same burst at
different resolutions.

The simulations appear to converge at a grid spacing of 3 cm, as decreasing it
to 1.5 cm does not significantly change the evolution. For comparison, Malone
et al. (2011) obtained convergence at a resolution of 0.5 cm, however the
X-ray burst under study had different initial conditions (notably a much
larger initial temperature). Despite convergence demonstrating the robustness
of the MAESTROeX calculations, our simulations do reach appreciable fractions
of the mach number ($\gtrsim$15%). Shortly thereafter, the calculation becomes
unstable222When Mach numbers are large, finding solutions to the velocity
divergence constraint becomes more difficult., and the simulations stop. At
this point, the thermonuclear runaway is still underway and the flux has not
had time to escape the atmosphere. Therefore, the full evolution of Type I
X-ray bursts in the low Mach approximation remains intractable. We discuss
this problem further in Section 4.

###  3.2 Evolution of convection

We now look in depth at the evolution of the 3 cm model. We identify three
“special” times for the three stages which we explore in this section. In
Figure 3, we plot snapshots of the fluid at these times in terms of its
vorticity, nuclear burning, hydrogen and carbon mass fractions.

Growth in the He layer. We first start the simulation by applying a gaussian
random noise perturbation of 1 part in 1000 to every temperature in the grid.
This triggers initial flows, which rapidly merge into a convective region in
the initial isentropic zone. In the top-left panel of Figure 3, convection can
be seen by the presence of large vortices, a distinct feature of 2D
convection. These vortices travel up and down the convective region, entrained
by the greater overturning convective motions. With convective velocities of
$\sim$km s-1, the convective turnover times at this point are on the order of
milliseconds. Burning is initially confined to the bottom of the convection
zone. Triple-$\alpha$ reactions convert helium into carbon and release heat
into the convective region, causing it to expand. After ${\sim}10$ ms, carbon
becomes fully mixed in the convective region, uniformily increasing in mass-
fraction as helium burns.

Figure 3: Flow, burning and composition at different times in the 3 cm
MAESTROeX simulation. From top to bottom: vorticity, burning rate, mass
fraction of hydrogen, mass fraction of carbon. Note the different scales in
the bottom two rows. The left column is during the initial growth of the
convection zone in the helium layer, the middle column is roughly at the
collision with the hydrogen layer, and the right column is later during the
expansion into the hydrogen. An animated version of this figure is available
at this link. Figure 4: Evolution of the temperature as a function of column
depth (note that the x-axis is reversed), at intervals of 10 ms. The colors
get darker as time increases. The labeled and thicker lines are the three
times from Figure 3. The circle and square markers show the convective
boundaries according to the Schwarzschild and Ledoux criteria respectively.
The vertical dotted line is the depletion depth, i.e. the depth of the
hydrogen layer. As convection approaches this line, the Schwarzschild and
Ledoux boundaries begin to separate. The inset zooms into the evolution near
the hydrogen boundary.

In Figure 4, we show the evolution of the temperature profile333Note that here
and in upcoming figures, the profiles of a quantity with depth are obtained
via horizontal averaging at each $z$ level in the grid.. Throughout this
initial evolution, the temperature gradient of the convection zone remains
adiabatic ($\nabla=d\ln T/d\ln p$ is close to $\nabla_{\rm ad}$, as indicated
in the figure). As a result, the top of the convection zone is roughly located
at the intersection between this gradient and the initial radiative profile
(but not exactly, as we will discuss in 3.3). Because the convection zone and
overlying stable He zone have similar mean molecular weights, the
Schwarzschild and Ledoux444We evaluate the Ledoux gradient $\nabla_{\rm L}$ as
in MESA, see Section 3.3 of Paxton et al. (2013). criteria agree on the
location of the convective boundary, which can be seen in Figure 4 with the
square and circle markers lying on top of each other.

Collision with the H shell. After about 50 ms, overshooting motions above the
convective boundary reach $y_{\rm d}$, the hydrogen depletion depth (Equation
1). This is the beginning of the event which we have dubbed the “collision”.
The second column of Figure 3 shows the state of the fluid 30 ms later, once
the convection zone has fully crossed $y=y_{\rm d}$. In the collision, fresh
fuel in the form of free protons become available and re-ignite the CNO cycle.
A thin burning layer appears at $y_{\rm d}$ where protons are capturing onto
${}^{12}{\rm C}$ previously made from helium burning, and then onto the fresh
${}^{13}{\rm N}$. The protons that do not burn at this boundary are able to
travel deeper into the convection zone, but not very far because the
temperature is so high. Even though carbon is burning away at the boundary,
more is being supplied by helium burning below, and mixing remains efficient
enough such that the carbon fraction is uniform throughout the convective
region.

Figure 5: Top: Evolution of the nuclear luminosity (Equation 4) with time, for
helium burning and metal burning reactions, and the total. The production of
the most important CNO rates is plotted in dashed lines. Bottom: Size of the
convective and semiconvective regions as a function of time. The dotted line
shows the depletion depth.

With the addition of protons into the mix, CNO burning overtakes helium-
burning as the main energy source. One way to track the total energy
production in the grid, using horizontally averaged quantities, is with the
“nuclear luminosity”,

$L_{\rm nuc}(t)=4\pi R^{2}\int_{z}\epsilon_{\rm
nuc}(z^{\prime},t)\rho(z^{\prime},t)dz^{\prime}\,,$ (4)

where the integral is performed over the vertical direction in the grid. This
quantity gives the total energy produced by the burst, at all depths, per unit
time555The actual radiative luminosity coming out of the atmosphere during the
burst is much smaller, both because of the thermal time to the surface, and
because much of the energy will get used to drive mass-loss.. We plot $L_{\rm
nuc}$ in the top panel of Figure 5. We see that after 50 ms, protons initially
capture onto ${}^{12}{\rm C}$ to produce ${}^{13}{\rm N}$, which subsequently
burns either via proton or $\alpha$-capture. Burning to heavier species such
as ${}^{17}{\rm F}$ begins ${\sim}50$ ms later. This sequence of reactions is
very similar to what we have previously obtained in 1D simulations.

The evolution of the convection following the collision in these simulation is
signicantly different than what we previously found in 1D. Instead of having
burning energy being deposited locally, resulting in a drastic expansion of
the convection in less than a millisecond (see Figure 3 of Guichandut &
Cumming, 2023), we instead see that protons are able travel into the
convection zone and burn anywhere within it. As seen in the bottom panel of
Figure 5, while the rate of expansion of the convection does increase
following the collision, it does so in a smooth fashion. This confirms our
speculations in Guichandut & Cumming (2023) that this violent expansion was
likely due to an improper treatment convection with MLT.

Growth into the H shell. Later, as the convection zone, which is becoming
richer in heavy elements, keeps growing into the hydrogen-rich shell, it does
so fighting against a buoyancy barrier. At this stage, composition gradients
start playing an important role in setting the convective boundary. This can
be seen in Figure 4 with the Ledoux and Schwarzschild boundaries becoming
separate. In order to grow, the convection zone has to “wait” until it is hot
enough that it is buoyant in the light material above. Therefore, a strong
temperature gradient develops at the convective boundary. This was predicted
by Weinberg et al. (2006). The material in this gradient region is thermally
unstable, but overall convectively stable due to the composition gradient.
These are the right conditions for semiconvection to develop (Kippenhahn et
al., 2012), shown as the blue regions in the bottom panel of Figure 5. To
model this type of mixing accurately would require following thermal diffusion
and also diffusive mixing of species, which MAESTROeX does not include. Note
however that we do not expect either thermal diffusion nor semiconvection to
prevent the existence of this sharp temperature gradient. The diffusion
coefficient for radiation is $4acT^{3}/(3\kappa\rho^{2}c_{p})\sim(10^{3}$ cm2
s-1) $T_{8}^{3}/\rho_{5}^{2}$, and thus thermal time over the ${\sim}10$ cm
gradient is 0.1 s. The diffusion coefficient for semiconvection is similar or
even smaller by a factor of up to 1000 (Langer, 1991).

During this period, the flow is dominated by large vortices (see top right
panel of Figure 3), which is a consequence of the inverse energy cascade of 2D
turbulence. Such large vortices have been shown to disappear in 3D simulations
of X-ray bursts (Zingale et al., 2015). Interestingly, these vortices appear
to carry within them a larger amount of hydrogen than the surrounding mean
flow, and are therefore burning at a faster rate than the horizontal mean.
Whether the presence of large vortices meaningfully affects the transport of
hydrogen and thefore the depth at which burning is taking place is a question
which can only be answered by 3D simulations using the same setup. We leave
this to future work.

Towards the end of the simulation, carbon and nitrogen have been fully
converted to heavier elements, and are no longer available to capture protons,
which are still inflowing from the top. The main proton-capture reaction
happening at this point, ${}^{16}{\rm O}(p,\gamma)^{17}{\rm F}$, is competing
with its photodisintegration inverse. As a result, an equilibrium mass-
fraction of ${\sim}10^{-3}$ of hydrogen is left behind, which is seen in the
third row of Figure 3. We note however that this is likely just an effect of
the limited size of our network — these protons would in reality have other
avenues to burn.

Throughout the simulation, a secondary convection zone has been developing in
the iron substrate underneath the burning layer, as can clearly be seen in the
bottom panel of Figure 5. This is the convective undershoot problem,
previously noted by Malone et al. (2011), where down currents in the
convective region are able to penetrate the stable layer despite opposing
buoyancy (which was even artificially enhanced by setting the temperature of
the iron substrate to a small value, Section 2.2). Chemical mixing between the
layers occurs, in which helium and carbon are brought down and a large amount
of iron gets dredged up into the burning layer (up to a mass fraction of
${\sim}50$% towards the end of the simulation), which limits the burning. It
would be interesting to investigate in future fluid simulations how much of
this is physical versus numerical. If correct, this dredge up of heavy
elements (not created in this burst but presumably in previous ones) could
have implications for the detection of spectral features in some photospheric
radius expansion bursts (in ’t Zand & Weinberg, 2010; Kajava et al., 2017;
Strohmayer et al., 2019).

###  3.3 Convective boundary mixing

In Guichandut & Cumming (2023), the interaction between the convection zone
and stable layer at the moment of collision caused an explosive transient with
convective velocities at the boundary suddenly increasing, and the convection
zone itself extending by many scale heights in the space of a few ${\mu}$s. It
was clear that these effects were numerical artefacts, caused by the
approximate treatment of convection in the regime of rapid heating. This
convective-reactive interaction is an example of a convective boundary mixing
process, which are an active area of research (see Anders & Pedersen 2023 for
a recent review). Now, we have the necessary fluid simulations to describe
what really happens at the top of the convection zone.

Starting with a simple question: How does the convection zone (CZ) grow into
the radiative zone (RZ)? The simplest model goes as follows: at any point in
time, convection is restricted to the adiabatic region, whereas the RZ above
is unchanging. Therefore, the way to move the convective boundary upwards is
by raising the adiabat of the CZ, i.e. increase its entropy. This is of course
accomplished by nuclear reactions which release heat into the mixture. As
such, the convective boundary is always located at the intersection between
the initial entropy profile of the radiative zone, $s_{\rm RZ,0}$666The
entropy of the radiative zone does increase over time, but much slower than
the evolution of the burst, see Hanawa & Sugimoto (1982)., and the current
entropy of the convection zone, $s_{\rm CZ}=s_{\rm CZ,0}+\Delta s_{\rm nuc}$.
This is well illustrated in Figure 4 of Hanawa & Sugimoto (1982).

It turns out that this simple picture is not at all what is happening in our
simulations. Surprisingly, we find that the specific entropy of the CZ is in
fact slowly decreasing over time. On its own, this behavior is easily
explained by the dredge-up of iron from below, which increases the mean
molecular weight faster than the temperature increases from nuclear burning.
Nevertheless, it hints at another mechanism for the growth of the CZ. The
answer starts in the inset of Figure 4. We see that above the convective
boundary, the fluid is becoming cooler than the initial RZ. Somehow, heat is
being transported away from this region. This reduces the entropy of the RZ,
allowing it to match to the CZ at smaller and smaller column depths — this is
how the CZ grows.

In Figure 6, we look at different properties of the flow near the convective
boundary, at $t=45$ ms. In the top panel, we notice a temperature bump above
the convective boundary. This bump is also clearly seen in Figure 4. The
origin of this temperature excess is a spike in burning at the same location,
as seen by the red dashed line. This is the same spike which appears in the
second row, middle panel of Figure 3, and is a result of the proton captures
onto ${}^{12}{\rm C}$ and ${}^{13}{\rm N}$. Indeed, in the middle of Figure 6,
we see that the root mean square (rms) vertical velocity does not go to zero
at the convective boundary; despite strong deceleration, some amount of fluid
is crossing the boundary and reaching all the way to the hydrogen layer at
$y=y_{\rm d}$. This is how carbon, generated from helium burning at the bottom
of the CZ, is able to travel to the hydrogen layer. This is confirmed by the
third panel in which we show the horizontally and time averaged flux of carbon
$F_{C}\equiv\rho v_{z}X(^{12}{\rm C})$, which stays positive all the way to
$y_{\rm d}$. In the other direction, the returning downflows carry with them
hydrogen from the top, which results in a negative hydrogen flux $F_{H}$. As
mentioned in the previous section, hydrogen does not burn only at the
boundary; much of it is able to travel downwards many scale heights. We
emphasize that this is all before the CZ has actually reached the hydrogen
layer — due to these overshooting motions, the collision is initiated earlier
than prescribed by the Schwarzschild criterion.

Figure 6: Bulk properties of the flow around the interface between the
convective and stable zone prior to the collision, i.e. the moment where this
boundary is below the hydrogen layer (at a column depth $y>y_{\rm d}$, the
vertical blue line). In each panel, the solid black line and the color dashed
lines are scaled on the left-hand and right-hand side respectively. Every
quantity is horizontally averaged over the domain and each height. Overbarred
quantities are time-averaged over 2 ms. Top: Temperature profile and burning
rate. Middle: RMS velocity fluctuation and adiabatic excess. The crossing of
the dashed orange with 0 sets the location of the vertical orange line
(Schwarzschild boundary). Bottom: Convective heat flux and composition fluxes
of hydrogen and carbon-12.

Let us now return to the issue of cooling above the convective boundary. In
the middle panel of Figure 6, we see that $\nabla$ remains just under
$\nabla_{\rm ad}$ over a short distance (about $20\%$ of the local scale
height), before dropping, as the temperature gradient goes to the initial
radiative gradient $\nabla_{\rm rad}$. This structure is well understood in
theories of convective boundary mixing. The marginally sub-adiabatic region
above the Schwarzschild boundary is known as the penetrative zone (PZ), while
the radiative region above with non-zero vertical velocity is the overshoot
zone (OZ) (Anders & Pedersen, 2023). In both regions, upflows are expanding
adiabatically in a sub-adiabatic background, thereby cooling their
surroundings. As a result, the convective heat flux, $F_{\rm conv}\equiv\rho
c_{p}v_{z}\Delta T$, where $\Delta T$ is the temperature difference relative
to the horizontal mean, is negative, as seen in the third panel of Figure 6.
Negative heat fluxes in the PZ have been observed in various simulations of
convection (e.g. Hurlburt et al., 1986; Singh et al., 1994; Browning et al.,
2004), but this is to our knowledge the first time the effect is observed in
concert with a growing convection zone. This cooling reduces the entropy in
the PZ and allows the underlying CZ to expand outwards.

##  4 Summary and Discussion

We have run low Mach number simulations of X-ray bursts in the regime where
helium is igniting and the convection zone eventually reaches new fuel in the
form of protons, which accelerates the thermonuclear runaway. Our main
findings are as follows:

  1. 1.

The resolution required to reach convergence is 3 cm per zone or smaller
(Figure 2). At coarser resolution, the evolution of the burst is delayed due
to improper resolving of the temperature peak.

  2. 2.

Nuclear burning throughout the burst evolves in a similar fashion to 1D
simulations. The burst is initially driven by helium burning. Once protons
enter the mix, there is a build up of ${}^{13}{\rm N}$, which then quickly
burns to ${}^{14}{\rm O}$, sharply increasing the nuclear energy output
(Figure 5).

  3. 3.

However, the evolution of convection following the collision is very different
than that in 1D simulations using mixing-length theory. While the rate of
growth of the convection zone does increase at this moment, driven by the
accelerated burning, it is not as violent as in 1D.

  4. 4.

The collision is also not a single precise moment, but rather a gradual
encroaching of the convection zone into the hydrogen shell, which is mediated
by convective boundary mixing processes (Figure 6). Early on, convective
motions penetrate through the convective boundary, exchanging entropy and
composition with the stable layer above. A negative heat flux reduces the
entropy above the boundary, which allows the convection zone to grow.

Our simulations are not able to follow the full evolution of convection during
the burst. This is mainly a dynamic range issue; there is a $\sim$4 order of
magnitude difference between the density at the bottom of the fuel layer
($y\sim 10^{8}$ g cm-2, $\rho\sim 10^{6}$ g cm-3) and that of the theoretical
maximal extent of the convection zone ($y\sim 10^{4}$ g cm-2, $\rho\sim
10^{2}$ g cm-3, e.g. Weinberg et al. (2006)). It is difficult to track this
range in the low Mach approximation because of the need for a velocity damping
region (the sponge) at the top of the model. The smallest values of $\rho_{\rm
cutoff}$ which would allow us to evolve the burst for a long time was $10^{3}$
g cm-3. Following the growth of the convection zone into regions of lower
density therefore likely requires a fully compressible fluid simulation. These
will allow us to constrain the duration of the observed “pause” as in SAX
J1808.4–3658 (Bult et al., 2019; Guichandut & Cumming, 2023).

Our simulations did not include thermal diffusion. Although this is possible
with MAESTROeX, it is not expected to be a big effect in the early stages of
the burst, where convection is taking place at high densities and the
conductivity is small However, it is thermal diffusion which ultimately stops
the advancement of convection at lower depths, when it becomes more efficient
as a means to transport energy than convection. Future simulations will
therefore need to include it to describe the full evolution of convection.

In addition, although previous comparisons between 2D and 3D simulations of
bursts have demonstrated good overall agreement in the burst evolution
(Zingale et al., 2015), it is well known that the turbulent flow behaves very
differently. In this work, we noticed the presence of large vortices which
carry a different composition than the surroundings. These might completely
disappear in 3D, which might alter the mixing and burning at different
heights. We remark that a fully compressible 3D simulation of this burst is
within reach with current computational resources, and we have already started
some experiments with the CASTRO code (Almgren et al., 2010) which show a
similar early burst evolution as our MAESTROeX simulations. Furthermore,
rotation could be important in dictating the flow patterns, as the spin period
of accreting neutron stars can be on the order of milliseconds, comparable to
the convective turnover time. These effects have recently been explored in the
Boussinesq approximation by Garcia et al. (2018), who modeled burning driven
convection over the whole surface of the star.

Lastly, a complete treatment of Type I X-ray bursts requires not only the
vertical, but also the lateral heat transport around the surface of the
neutron star. This is the well studied flame-propagation problem, which has
been applied to pure He (e.g. Cavecchi et al., 2013; Eiden et al., 2020; Chen
et al., 2023) and most recently to mixed H/He atmospheres (Johnson & Zingale,
2024), but not yet to the H/He layered case of this paper. How does our
concept of “collision” fit into this alternate picture? Is the flame fast
enough to propagate before convection has had time to connect with the
hydrogen shell, or do these processes interact in some other way? We hope our
work motivates future explorations into this question.

MAESTROeX is open-source and under active development. The code as well as the
inputs required to run the simulations presented in this paper can be found at
https://github.com/amrex-Astro/Maestroex.

We thank J.R. Fuentes for useful discussions on convective boundaries, E.T.
Johnson and A. Smith Clark for technical assistance. This project was borne
out of a visit at Stony Brook University, thanks to funding by the Centre de
Recherche en Astrophysique du Québec (CRAQ). S.G. is also supported by an
NSERC scholarship. The work at Stony Brook was supported by DOE/Office of
Nuclear Physics grant DE-FG02-87ER40317. This work was supported by NSERC
Discovery Grant RGPIN-2023-03620. This research used resources of the National
Energy Research Scientific Computing Center (NERSC), a Department of Energy
Office of Science User Facility using NERSC award NP-ERCAP0027167.

## References

  * Almgren et al. (2006) Almgren, A. S., Bell, J. B., Rendleman, C. A., & Zingale, M. 2006, ApJ, 649, 927, doi: 10.1086/507089
  * Almgren et al. (2010) Almgren, A. S., Beckner, V. E., Bell, J. B., et al. 2010, ApJ, 715, 1221, doi: 10.1088/0004-637X/715/2/1221
  * Anders & Pedersen (2023) Anders, E. H., & Pedersen, M. G. 2023, Galaxies, 11, 1, doi: 10.3390/galaxies11020056
  * Browning et al. (2004) Browning, M. K., Brun, A. S., & Toomre, J. 2004, Eur. Sp. Agency, (Special Publ. ESA SP, 349
  * Bult et al. (2019) Bult, P., Jaisawal, G. K., Güver, T., et al. 2019, ApJ, 885, L1, doi: 10.3847/2041-8213/ab4ae1
  * Cavecchi et al. (2013) Cavecchi, Y., Watts, A. L., Braithwaite, J., & Levin, Y. 2013, MNRAS, 434, 3526, doi: 10.1093/mnras/stt1273
  * Chen et al. (2023) Chen, Z., Zingale, M., & Eiden, K. 2023, ApJ, 955, 128, doi: 10.3847/1538-4357/acec72
  * Cox & Giuli (1968) Cox, J. P., & Giuli, R. T. 1968, Principles of stellar structure (New York: Gordon and Breach)
  * Cumming & Bildsten (2000) Cumming, A., & Bildsten, L. 2000, ApJ, 544, 453
  * Eiden et al. (2020) Eiden, K., Zingale, M., Harpole, A., et al. 2020, ApJ, 894, 6, doi: 10.3847/1538-4357/ab80bc
  * Fan et al. (2019) Fan, D., Nonaka, A., Almgren, A. S., Harpole, A., & Zingale, M. 2019, ApJ, 887, 212, doi: 10.3847/1538-4357/ab4f75
  * Galloway & Keek (2021) Galloway, D. K., & Keek, L. 2021, in Timing Neutron Stars Pulsations, Oscil. Explos., ed. T. M. Belloni, M. Mendez, & C. Zhang (Springer, Berlin, Heidelberg), 209–262, doi: 10.1007/978-3-662-62110-3_5
  * Galloway et al. (2020) Galloway, D. K., in ’t Zand, J. J. M., Chenevez, J., et al. 2020, ApJS, 249, 32, doi: 10.26180/5e4a697d9b8b6
  * Garcia et al. (2018) Garcia, F., Chambers, F. R., & Watts, A. L. 2018, Phys. Rev. Fluids, 3, doi: 10.1103/PhysRevFluids.3.123501
  * Guichandut & Cumming (2023) Guichandut, S., & Cumming, A. 2023, ApJ, 954, 54, doi: 10.3847/1538-4357/ace43c
  * Hanawa & Sugimoto (1982) Hanawa, T., & Sugimoto, D. 1982, PASJ, 34, 1
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Henyey et al. (1965) Henyey, L., Vardya, M. S., & Bodenheimer, P. 1965, ApJ, 142, 841
  * Herwig et al. (2011) Herwig, F., Pignatari, M., Woodward, P. R., et al. 2011, ApJ, 727, 89, doi: 10.1088/0004-637X/727/2/89
  * Herwig et al. (2014) Herwig, F., Woodward, P. R., Lin, P. H., Knox, M., & Fryer, C. 2014, ApJL, 792, 1, doi: 10.1088/2041-8205/792/1/L3
  * Hunter (2007) Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90
  * Hurlburt et al. (1986) Hurlburt, N. E., Toomre, J., & Massaguer, J. M. 1986, ApJ, 311, 563, doi: 10.1086/164796
  * in ’t Zand & Weinberg (2010) in ’t Zand, J. J. M., & Weinberg, N. N. 2010, A&A, 520, A81, doi: 10.1051/0004-6361/200913952
  * Johnson & Zingale (2024) Johnson, E. T., & Zingale, M. 2024.  https://arxiv.org/abs/2402.02013
  * José et al. (2010) José, J., Moreno, F., Parikh, A., & Iliadis, C. 2010, ApJS, 189, 204, doi: 10.1088/0067-0049/189/1/204
  * Joss (1977) Joss, P. C. 1977, Nature, 270, 310, doi: 10.1038/270310a0
  * Joss (1978) —. 1978, ApJ, 225, L123, doi: 10.1086/182808
  * Kajava et al. (2017) Kajava, J. J. E., Nättilä, J., Poutanen, J., et al. 2017, MNRAS, 464, L6, doi: 10.1093/mnrasl/slw167
  * Kippenhahn et al. (2012) Kippenhahn, R., Weigert, A., & Weiss, A. 2012, Stellar Structure and Evolution, Astronomy and Astrophysics Library (Berlin, Heidelberg: Springer Berlin Heidelberg), doi: 10.1007/978-3-642-30304-3
  * Kondev et al. (2021) Kondev, F. G., Wang, M., Huang, W. J., Naimi, S., & Audi, G. 2021, Chinese Phys. C, 45, 15, doi: 10.1088/1674-1137/abddae
  * Langer (1991) Langer, N. 1991, A&A, 252, 669
  * Lewin et al. (1993) Lewin, W. H. G., van Paradijs, J., & Taam, R. E. 1993, SSRv, 62, 223, doi: 10.1007/BF00196124
  * Lin et al. (2006) Lin, D. J., Bayliss, A., & Taam, R. E. 2006, ApJ, 653, 545, doi: 10.1086/508863
  * Malone et al. (2011) Malone, C. M., Nonaka, A., Almgren, A. S., Bell, J. B., & Zingale, M. 2011, ApJ, 728, 118, doi: 10.1088/0004-637X/728/2/118
  * Malone et al. (2014) Malone, C. M., Zingale, M., Nonaka, A., Almgren, A. S., & Bell, J. B. 2014, ApJ, 788, 115, doi: 10.1088/0004-637X/788/2/115
  * Nonaka et al. (2010) Nonaka, A., Almgren, A. S., Bell, J. B., et al. 2010, ApJS, 188, 358, doi: 10.1088/0067-0049/188/2/358
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3, doi: 10.1088/0067-0049/192/1/3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4, doi: 10.1088/0067-0049/208/1/4
  * Singh et al. (1994) Singh, H. P., Roxburgh, I. W., & Chan, K. L. 1994, A&A, 281, L73
  * Smith et al. (2023) Smith, A. I., Johnson, E. T., Chen, Z., et al. 2023, ApJ, 947, 65, doi: 10.3847/1538-4357/acbaff
  * Strohmayer & Bildsten (2006) Strohmayer, T., & Bildsten, L. 2006, in Compact Stellar X-ray Sources, ed. W. Lewin & M. van der Klis (Cambridge University Press), 113–156, doi: 10.1017/CBO9780511536281.004
  * Strohmayer et al. (2019) Strohmayer, T. E., Altamirano, D., Arzoumanian, Z., et al. 2019, ApJL, 878, L27, doi: 10.3847/2041-8213/ab25eb
  * Timmes & Swesty (2000) Timmes, F. X., & Swesty, F. D. 2000, ApJS, 126, 501, doi: 10.1086/313304
  * Turk et al. (2011) Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, doi: 10.1088/0067-0049/192/1/9
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wallace et al. (1982) Wallace, R. K., Woosley, S. E., & Weaver, T. A. 1982, ApJ, 258, 696, doi: 10.1086/160119
  * Weinberg et al. (2006) Weinberg, N. N., Bildsten, L., & Schatz, H. 2006, ApJ, 639, 1018, doi: 10.1086/499426
  * Woosley et al. (2004) Woosley, S. E., Heger, A., Cumming, A., et al. 2004, ApJS, 151, 75, doi: 10.1086/381533
  * Yu & Weinberg (2018) Yu, H., & Weinberg, N. N. 2018, ApJ, 863, 53, doi: 10.3847/1538-4357/aad045
  * Zhang et al. (2019) Zhang, W., Almgren, A., Beckner, V., et al. 2019, J. Open Source Softw., 4, 1370, doi: 10.21105/joss.01370
  * Zingale et al. (2015) Zingale, M., Malone, C. M., Nonaka, A., Almgren, A. S., & Bell, J. B. 2015, ApJ, 807, 60, doi: 10.1088/0004-637X/807/1/60

##  Appendix A Building the initial model

Our model is an extension of the one used in Malone et al. (2014), in which
the main feature is an isentropic zone in the fuel layer which becomes
convectively unstable once the simulation begins. For this, we provide the
temperature of the “neutron star” $T_{\rm star}$, i.e. the temperature in the
iron substrate, the height $H_{\rm star}$ of this substrate, the temperature
and column depth of the base of the fuel layer, $T_{\rm base}$ and $y_{\rm
base}$. At $r=H_{\rm star}$, the composition transitions from pure iron to the
fuel composition with a hyperbolic tangent of width $\delta$ (see appendix of
Malone et al. (2014)). We set $\delta=5$ cm so that this transition is thin,
but well resolved at the largest grid resolutions.

In Malone et al. (2014), the fuel layer was a mixture of H and He at fixed
mass fractions. Here, the composition depends on the radial location. In
particular, the mass-fraction of hydrogen $X$ depends analytically on the
column depth $y$ (Cumming & Bildsten, 2000),

$X(y)=\rm{max}(0,\,X_{0}[1-y/y_{\rm d}])\,$ (A1)

where $X_{0}$, the initial hydrogen fraction, and $y_{\rm d}$, the hydrogen
depletion depth, are two more input parameters. At any point during the
building of the model, $y$ can be obtained from the pressure, since $p=gy$ in
HSE. For the CNO species, we set constant mass fractions of ${}^{14}{\rm O}$
and ${}^{15}{\rm O}$, with a sum of $Z_{\rm CNO}$, and a ratio of
70.621/122.266, which is the ratio of beta-decay half-lives of these isotopes
(Kondev et al., 2021). Finally, the He mass-fraction is simply $Y=1-X-Z_{\rm
CNO}$.

The initial temperature profile is also more complicated than in Malone et al.
(2014). Above the isentropic zone where the temperature gradient,$\nabla\equiv
d\ln T/d\ln y$, is the adiabatic gradient $\nabla_{\rm ad}$, the atmosphere is
in radiative equilibrium, $dT^{4}/dr\propto-F$, where $F$ is the radiative
flux, which depends on the depth and nuclear burning from steady-state energy
balance, $dF/dy=-\epsilon_{\rm nuc}$. In practice, instead of solving for an
exact temperature profile consistent with the burning, we divide the fuel
layer into three zones, and obtain the best-fitting value of $\nabla$ in each
zone from the MESA model. This is illustrated in Figure 7. On top of the three
$\nabla$’s, the column depth boundaries between the first and second zone,
$y_{2}$, and second and third zone, $y_{3}$, also must be specified as input
parameters.

Figure 7: Illustration of the temperature profile considered for the toy
model. In each zone demarcated by the vertical dotted lines (except for the
isothermal substrate at large depths $y>y_{1}$, the best-fitting temperature
gradient $\nabla=d\ln T/d\ln y$ is obtained from the MESA model. This
temperature profile is then fed into a hydrostatic solver which determines the
density at every step, while the composition is an analytical function of the
pressure with Equation A1. See text for further details.

In 1D calculations, ignition is manifested by a rapid reduction in the
timestep which is limited by the nuclear energy generation rate. However, in a
fluid simulation, the timestep is restricted from the start by the fluid
velocity (Courant condition). Therefore, a simulation starting at this point
will take a long time (many timesteps) to evolve and proceed into the
thermonuclear runaway. In order to accelerate it, we “kick” the initial model,
i.e. we artifically increase the temperature at the base of the fuel layer
($T_{\rm base}$ parameter), while keeping the temperature gradient in that
region $\nabla_{1}$ constant. As a result, the boundary of the adiabatic
region moves according to

$y_{2}^{\prime}=y_{2}\left(\frac{T_{\rm base}}{T_{\rm
base}^{\prime}}\right)^{(\nabla_{1}-\nabla_{2})^{-1}}\,,$ (A2)

where the prime (′) values are those of the kicked model.

The table below lists all of the parameters used to build the initial fluid
model.

Table 1: Input parameters for initial model

Parameter | Definiton | Value | Units
---|---|---|---
$T_{\rm star}$ | Temperature at the bottom of the atmosphere | $10^{8}$ | K
$H_{\rm star}$ | Height of the iron substrate | 1070 | cm
$\delta$ | Width of the hyperbolic tangent transition between iron and fuel | 5 | cm
$T_{\rm base}$ | Temperature at the bottom of the fuel layer | $3.5\times 10^{8}\,{}^{\rm(a)}$ | K
$\nabla_{1}$ | Temperature gradient in the He isentropic zone | $0.37\,^{\rm(b)}$ |
$\nabla_{2}$ | Temperature gradient in the inert He zone | 0.12 |
$\nabla_{3}$ | Temperature gradient in the H-rich zone | 0.22 |
$y_{\rm base}$ | Column depth at the base of the fuel layer | $5.6\times 10^{7}$ | g cm-2
$y_{2}$ | Column depth boundary between the $\nabla_{1}$ and $\nabla_{2}$ zones | $5.6\times 10^{7}\,{}^{\rm(a)}$ | g cm-2
$y_{3}$ | Column depth boundary between the $\nabla_{2}$ and $\nabla_{3}$ zones | $2\times 10^{7}$ | g cm-2
$X_{0}$ | Hydrogen fraction in the accreted material | 0.7 |
$Z_{\rm CNO}$ | Total CNO fraction in the accreted material | 0.02 |
$y_{\rm d}$ | Hydrogen depletion depth from stable hot CNO burning | $3.7\times 10^{7}$ | g cm-2
$T_{\rm cutoff}$ | Minimum temperature in the model | $10^{6}$ | K
$\rho_{\rm cutoff}$ | Minimum density in the model | $10^{3}$ | g cm-3

(a) These are the values used for the kicked model, i.e. the one shown in Fig.
1, not Fig. 7.
(b) This is approximately the correct value of $\nabla_{\rm ad}$ in that
region.


[affiliation=1]YuiSudo [affiliation=1]YosukeFukumoto
[affiliation=1]MuhammadShakeel [affiliation=2]YifanPeng
[affiliation=2]ShinjiWatanabe

# Contextualized Automatic Speech Recognition with Dynamic Vocabulary

###### Abstract

Deep biasing (DB) improves the performance of end-to-end automatic speech
recognition (E2E-ASR) for rare words or contextual phrases using a bias list.
However, most existing methods treat bias phrases as sequences of subwords in
a predefined static vocabulary, which can result in ineffective learning of
the dependencies between subwords. More advanced techniques address this
problem by incorporating additional text data, which increases the overall
workload. This paper proposes a dynamic vocabulary where phrase-level bias
tokens can be added during the inference phase. Each bias token represents an
entire bias phrase within a single token, thereby eliminating the need to
learn the dependencies between the subwords within the bias phrases. This
method can be applied to various architectures because it only extends the
embedding and output layers in common E2E-ASR architectures. Experimental
results demonstrate that the proposed method improves the performance of bias
phrases on English and Japanese datasets.

###### keywords:

speech recognition, contextualization, biasing, dynamic vocabulary

##  1 Introduction

End-to-end automatic speech recognition (E2E-ASR) [1]methods convert acoustic
feature sequences directly to token sequences without requiring the components
used in conventional ASR systems, e.g., acoustic models and language models
(LM). Various E2E-ASR methods have been proposed, including connectionist
temporal classification (CTC) [2, 3], recurrent neural network transducer
(RNN-T) [4, 5], attention mechanisms [6, 7], and their hybrids [8, 9, 10, 11].
However, the effectiveness of E2E-ASR models is strongly dependent on the
context of the training data, which can result in performance inconsistencies
for specific user contexts. However, retraining models for every possible
context is infeasible, which highlights the need for a method that allows
users and developers to contextualize the model efficiently without additional
training process.

To address the challenge of contextualization in ASR, a common strategy is
shallow fusion with an external LM [12, 13], which frequently involves using a
weighted finite state transducer (WFST) to create a specialized in-class LM
that improves the recognition of the target named entities. In addition,
neural LM fusion methods [14, 15] attempt to improve accuracy by integrating
an external neural LM with an E2E-ASR system. This integration involves
rescoring the output hypotheses of the E2E-ASR model; however, integrating an
external LM requires additional training, which increases the overall
complexity.

Deep biasing (DB) [16, 17, 18, 19, 20, 21, 22] provides an efficient strategy
for contextualizing E2E-ASR models without additional retraining processes.
This approach uses an editable list of phrases, referred to as a bias list, to
tailor the E2E-ASR model to specific contexts. A key feature of many DB
techniques is the use of a cross-attention layer to ensure accurate
recognition of bias phrases. To improve the effectiveness of this layer,
previous studies [18, 19, 20, 21, 22] have introduced an auxiliary loss
function specifically for bias phrase detection. This function facilitates the
optimization of the parameters of the cross-attention layer through multitask
learning. However, these modules are tightly integrated into individual
architectures, making their application to other architectures complex, and
multitask learning requires multiple experimental training phases to fine-tune
the learning weights, which is a time-consuming process.

In addition, most existing techniques process bias phrases as sequences of
subwords in a predefined static vocabulary, which can lead to inadequate
learning of the dependencies between the subwords, especially for rare
phrases. For example, the personal name ``Nelly'' might be segmented into a
subword sequence, e.g., ``N'', ``el'', ``ly''. However, this approach
frequently fails to capture the dependencies between these subwords. To
address this problem, prefix tree-based approaches [23, 24, 25] have been
proposed to efficiently organize bias phrases with common prefixes into a tree
structure to manage large bias lists effectively. This tree structure is
relatively lightweight to create; however, it requires additional processing.
Advanced strategies have incorporated additional text data during training to
improve the handling of large bias lists. For example, previous studies [25,
26] have demonstrated significant performance improvements in terms of bias
phrase recognition by training an LM or text encoder using additional text
data. However, this technique can increase the workload considerably, and in
some cases the acquisition of additional data may be infeasible. Other
attempts have been proposed to improve contextualization using extra
information, such as phonemes for the bias phrases [27, 28, 29], named entity
tags [30], and synthesized speech [31, 32]; however, these strategies also
increase the workload.

This paper proposes a simple but effective deep biasing method by introducing
a dynamic vocabulary where phrase-level bias tokens can be added during the
inference phase. Each bias token represents an entire bias phrase within a
single token, thereby bypassing the complex process of learning subword
dependencies within the bias phrases. This approach allows for effective
biasing without relying on additional resources such as text data and tree
structures, unlike the previous methods [25, 26]. In addition, the proposed
bias tokens are integrated seamlessly with the conventional tokens in the
embedding and output layers, which eliminates the need for auxiliary losses,
unlike the previous studies [18, 19, 20, 21, 22]. Furthermore, unlike the
previous methods [19, 20, 21, 22, 25, 26], the proposed method can be applied
to various architectures because it only extends the embedding and output
layers of the common E2E-ASR models.

##  2 End-to-end ASR

This section describes an attention-based E2E-ASR system, which is extended to
the proposed method.

###  2.1 Audio encoder

The audio encoder comprises two convolutional layers, a linear projection
layer, and $M_{\text{a}}$ Conformer blocks [5]. The Conformer encoder
transforms an audio feature sequence $\bm{X}$ to a $T$ length hidden state
vector sequence $\bm{H}=[\bm{h}_{1},...,\bm{h}_{T}]\in\mathbb{R}^{d\times T}$
where $d$ represents the dimension:

$\bm{H}=\mathrm{AudioEnc}(\bm{X}).$ (1)

###  2.2 Decoder

Given $\bm{H}$ generated by the audio encoder in Eq. (1) and the previously
estimated token sequence $y_{0:s-1}=[y_{0},...,y_{s-1}]$, the decoder
estimates the next token $y_{s}$ recursively as follows:

$P(y_{s}|y_{0:s-1},\bm{X})=\mathrm{Decoder}(y_{0:s-1},\bm{H}),$ (2)

where $y_{s}$ is the $s$-th subword-level token in the pre-defined static
vocabulary $\mathcal{V}^{\text{n}}$ of size $K$
($y_{s}\in\mathcal{V}^{\text{n}}$). Here, $y_{0}$ represents the $<$sos$>$
token. Specifically, the decoder comprises an embedding layer, $M_{d}$
transformer blocks, and an output layer. First, the embedding layer with
positional encoding converts the input token sequence $y_{0:s-1}$ to an
embedding vector sequence
$\bm{E}_{0:s-1}=[\bm{e}_{0},...,\bm{e}_{s-1}]\in\mathbb{R}^{d\times s}$ as
follows:

$\bm{E}_{0:s-1}=\mathrm{Embedding}(y_{0:s-1}).$ (3)

Then, $\bm{E}_{0:s-1}$ is input to the $M_{d}$ transformer blocks with the
hidden state vectors $\bm{H}$ in Eq. (1) to generate a hidden state vector
$\bm{u}_{s}$ as follows:

$\bm{u}_{s}=\mathrm{Transformer}(\bm{E}_{0:s-1},\bm{H}),$ (4)

Subsequently, a token-wise score
$\bm{\alpha}^{\text{n}}=[\alpha^{\text{n}}_{1},...,\alpha^{\text{n}}_{K}]$ and
corresponding probability are calculated as follows:

$\bm{\alpha}^{\text{n}}=\mathrm{Linear}(\bm{u}_{s}),$ (5)
$P\left(y_{s}\mid\bm{y}_{0:s-1},\bm{X}\right)=\mathrm{Softmax}(\bm{\alpha}^{\text{n}}),$
(6)

By repeating these processes recursively, the posterior probability is
formulated as follows:

$P(y_{0:S}\mid\bm{X})=\prod_{s=1}^{S}P\left(y_{s}\mid
y_{0:s-1},\bm{X}\right),$ (7)

where $S$ is the total number of tokens. Model parameters are optimized by
minimizing the negative log-likelihood described as follows:

$L=-\log P(y_{0:S}\mid\bm{X}).$ (8)

The embedding and output layers are extended to the proposed biasing method in
Section 3.2.

##  3 Proposed method

Figure 1 shows the overall architecture of the proposed method, which
comprises the audio encoder, bias encoder, and decoder with extended embedding
and output layers. These components are described in the following
subsections.

Figure 1: Overall architecture of the proposed method, including the audio
encoder, bias encoder, and bias decoder with extended embedding and output
layers.

###  3.1 Bias encoder

Similar to [22], the bias encoder comprises an embedding layer, $M_{\text{e}}$
transformer blocks, a mean pooling layer, and a bias list
$B=\\{b_{1},\cdots,b_{N}$}, where $b_{n}\in\mathcal{V}^{\text{n}}$ is the
$I$-lengnth subword token sequence of the $n$-th bias phrase (e.g., [$<$N$>$,
$<$el$>$, $<$ly$>$]). After converting the bias list $B$ into a matrix
$\bm{B}\in\mathbb{R}^{I_{\text{max}}\times N}$ through zero padding based on
the maximum token length $I_{\text{max}}$ in $B$, the embedding layer and the
$M_{\text{e}}$ transformer blocks extract the high level representation
$\bm{G}\in\mathbb{R}^{d\times I_{\text{max}}\times N}$ as follows:

$\bm{G}=\mathrm{TransformerEnc}(\mathrm{Embedding}(\bm{B})).$ (9)

Then, a mean pooling layer extracts phrase-level embedding vectors
$\bm{V}=[\bm{v}_{1},\cdots,\bm{v}_{N}]\in\mathbb{R}^{d\times N}$ as follows:

$\bm{V}=\mathrm{MeanPool}(\bm{G}).$ (10)

###  3.2 Bias decoder with dynamic vocabulary

(a) Embedding layer

(b) Output layer

Figure 2: Extended embedding and output layers. (a) If the input token is a
dynamic bias token, the corresponding phrase-level bias embedding $\bm{v}_{n}$
is extracted. (b) The bias phrase score $\bm{\alpha}^{\text{b}}$ is calculated
using the inner product.

To avoid the complexity of learning dependencies within the bias phrases, we
introduce a dynamic vocabulary
$\mathcal{V}^{\text{b}}=\\{<b_{1}>,\cdots,<b_{N}>\\}$ where the phrase-level
bias tokens represent the bias phrases with the single entities for the $N$
bias phrases in the bias list $\bm{B}$. Unlike Eq. (2), the bias decoder
estimates the next token
$y^{\prime}_{s}\in\\{\mathcal{V}^{\text{n}}\cup\mathcal{V}^{\text{b}}\\}$
given $\bm{H}$, $\bm{V}$ in Eqs. (1) and (10), and $y^{\prime}_{0:s-1}$ as
follows:

$P(y^{\prime}_{s}|y^{\prime}_{0:s-1},\bm{X},\bm{B})=\mathrm{BiasDecoder}(y^{\prime}_{0:s-1},\bm{H},\bm{V}).$
(11)

For example, if a bias phrase ``Nelly'' exists in the bias list, the bias
decoder outputs the corresponding bias token [$<$Nelly$>$] rather than the
normal token sequence [$<$N$>$, $<$el$>$, $<$ly$>$].

The bias decoder comprises an extended embedding layer, $M_{d}$ transformer
blocks, and an extended output layer. First, the input token sequence
$y^{\prime}_{0:s-1}$ is converted to the embedding vector sequence
$\bm{E}^{\prime}_{0:s-1}=[\bm{e}^{\prime}_{0},...,\bm{e}^{\prime}_{s-1}]\in\mathbb{R}^{d\times
s}$. Unlike Eq. (3), if the input token $y^{\prime}_{s-1}$ is a bias token,
the corresponding bias embedding $\bm{v}_{n}$ is extracted from $\bm{V}$
(Figure 2(a)), otherwise, the embedding layer is used with a linear layer as
follows:

$\bm{e}^{\prime}_{s-1}=\begin{cases}\mathrm{Linear}(\mathrm{Extract}(\bm{V},y^{\prime}_{s-1}))&(y^{\prime}_{s-1}\in\mathcal{V}^{\text{b}})\\\
\mathrm{Linear}(\mathrm{Embedding}(y^{\prime}_{s-1}))&(y^{\prime}_{s-1}\in\mathcal{V}^{\text{n}}).\end{cases}$
(12)

Then, the $M_{d}$ transformer blocks convert $\bm{E}^{\prime}_{0:s-1}$ to the
hidden state vector $\bm{u}^{\prime}_{s}$ as in Eq. (4). In addition to the
normal token score $\bm{\alpha}^{\text{n}}$ in Eq. (5), the bias token score
$\bm{\alpha}^{\text{b}}=[\alpha^{\text{b}}_{1},...,\alpha^{\text{b}}_{N}]$ is
calculated using an inner product with linear layers (Figure 2(b)):

$\bm{\alpha}^{\text{b}}=\frac{\mathrm{Linear}(\bm{u}^{\prime}_{s})\mathrm{Linear}(\bm{V}^{T})}{\sqrt{d}}.$
(13)

After $\bm{\alpha}^{\text{b}}$ is concatenated with $\bm{\alpha}^{\text{n}}$,
resulting in
$\bm{\alpha}=[\alpha^{\text{n}}_{1},...,\alpha^{\text{n}}_{K},\alpha^{\text{b}}_{1},...,\alpha^{\text{b}}_{N}]$,
the softmax function calculates the token-wise probability, including the bias
token, as follows:

$P\left(y^{\prime}_{s}\mid
y^{\prime}_{0:s-1},\bm{X},\bm{B}\right)=\mathrm{Softmax}(\mathrm{Concat}(\bm{\alpha}^{\text{n}},\bm{\alpha}^{\text{b}})).$
(14)

Similar to Eqs. (7) and (8), the posterior probability and the loss function
are formulated as follows:

$P(y^{\prime}_{0:S^{\prime}}\mid\bm{X},\bm{B})=\prod_{s=1}^{S^{\prime}}P\left(y^{\prime}_{s}\mid
y^{\prime}_{0:s-1},\bm{X},\bm{B}\right),$ (15) $L^{\prime}=-\log
P(y^{\prime}_{0:S^{\prime}}\mid\bm{X},\bm{B}).$ (16)

where $S^{\prime}$ represents the total number of tokens based on the proposed
bias tokens, respectively. Note that Eqs. (13) and (14) are flexible to the
bias list size $N$, and the proposed method is optimized using Eq. (16)
without auxiliary losses.

The proposed method can be applied to other E2E-ASR architectures, such as
CTC, RNN-T, and their hybrids, including streaming systems and multilingual
systems [8, 9, 10, 33, 34, 35, 36, 37], because the proposed method only
extends the embedding and output layers in addition to the bias encoder
(Figure 3).

(a) CTC-based biasing

(b) RNN-T-based biasing

Figure 3: Various architectures utilized in the proposed method.

###  3.3 Training

During the training phase, a bias list $\bm{B}$ is created randomly from the
reference transcriptions for each batch, where $N_{\text{utt}}$ bias phrases
are selected per utterance, each has $I$ tokens in length. This process yields
a total of $N$ bias phrases calculated as $N_{\text{utt}}\times$ batch size.
Once the bias list $\bm{B}$ is defined, the corresponding reference
transcriptions are replaced with the bias tokens. For example, if ``$<$N$>$,
$<$el$>$, $<$ly$>$'' is extracted from the complete utterance $y$ = [$<$Hi$>$,
$<$N$>$, $<$el$>$, $<$ly$>$], the reference transcription is modified to
$y^{\prime}$ = [$<$Hi$>$, $<$Nelly$>$].

###  3.4 Bias weight during inference

Considering practicality, we introduce a bias weight to Eq. (14) to avoid
over/under-biasing during inference as follows:

$\mathrm{WeightSoftmax}_{i}(\bm{\alpha},\bm{w})=\frac{w_{i}\mathrm{exp}(\alpha_{i})}{\sum_{j=1}^{(K+N)}w_{j}\mathrm{exp}(\alpha_{j})},$
(17)

where $\bm{w}=[w_{1}....,w_{(K+N)}]$ and $i$ represent a weight vector and its
index for
$\bm{\alpha}=[\alpha^{\text{n}}_{1},...,\alpha^{\text{n}}_{K},\alpha^{\text{b}}_{1},...,\alpha^{\text{b}}_{N}]$,
respectively. The same bias weight $\mu$ is applied to bias tokens as follows:

$w_{i}=\begin{cases}1.0&(i\leqq K)\\\ \mu&(i>K).\end{cases}$ (18)

##  4 Experiment

###  4.1 Experimental setup

To demonstrate the adaptability of the proposed method, we applied it to
offline CTC/attention [8] and streaming RNN-T models. The input features are
80-dimensional Mel filterbanks with a window size of 512 samples and a hop
length of 160 samples. Then, SpecAugment is applied. The audio encoder
comprises two convolutional layers with a stride of two and a 256-dimensional
linear projection layer followed by 12 conformer layers with 1024 linear units
and layer normalization. For the streaming RNN-T, the audio encoder is
processed blockwisely [38] with the block size and look ahead of 800 and 320
ms, respectively. The bias encoder has six transformer blocks with 1024 linear
units. Regarding the bias decoder, the offline CTC/attention model has six
transformer blocks with 2048 linear units, and the streaming RNN-T model has a
single long short-term memory layer with a hidden size of 256 and a linear
layer of 320 joint sizes for prediction and joint networks. The attention
layers in the audio encoder, the bias encoder, and the bias decoder are four-
multihead attentions with a dimension $d$ of 256. The offline CTC/attention
and streaming RNN-T models have 40.58 M and 31.38 M parameters, respectively,
including the bias encoders. The bias weight $\mu$ (Section 3.4) is set to 0.8
(this is discussed further in Section 4.4). During training, a bias list
$\bm{B}$ is created randomly for each batch with $N_{\text{utt}}$ = [2 - 10]
and $I$ = [2 - 10] (Section 3.3). The proposed models are trained 150 epochs
at a learning rate of 0.0025/0.002 for the CTC/attention-based and RNN-T-based
systems, respectively.

The Librispeech-960 corpus [39] is used to evaluate the proposed method using
ESPnet as the E2E-ASR toolkit [40]. The proposed method is evaluated in terms
of word error rate (WER), bias phrase WER (B-WER), and unbiased phrase WER
(U-WER) [25]. Note that insertion errors are counted toward the B-WER if the
inserted phrases are present in the bias list; otherwise, insertion errors are
counted toward the U-WER.

Table 1: WER results of offline CTC/attention-based systems obtained on
Librispeech-960 (U-WER/B-WER).

|  $N$ = 0 (no-bias) |  $N$ = 100 |  $N$ = 500 |  $N$ = 1000
---|---|---|---|---
Model | test-clean | test-other | test-clean | test-other | test-clean | test-other | test-clean | test-other
Baseline [8] | 2.57 | 5.98 | 2.57 | 5.98 | 2.57 | 5.98 | 2.57 | 5.98
(CTC/attention) | (1.5/10.9) | (4.0/23.1) | (1.5/10.9) | (4.0/23.1) | (1.5/10.9) | (4.0/23.1) | (1.5/10.9) | (4.0/23.1)
CPPNet [18] | 4.29 | 9.16 | 3.40 | 7.77 | 3.68 | 8.31 | 3.81 | 8.75
| (2.6/18.3) | (5.9/37.5) | (2.6/10.4) | (6.0/23.0) | (2.8/10.9) | (6.5/24.3) | (2.9/11.4) | (6.9/25.3)
Attention-based DB | 5.05 | 8.81 | 2.75 | 5.60 | 3.21 | 6.28 | 3.47 | 7.34
\+ BPB beam search [22] | (3.9/14.1) | (6.6/27.9) | (2.3/6.0) | (4.9/12.0) | (2.7/7.0) | (5.5/13.5) | (3.0/7.7) | (6.4/15.8)
Proposed | 3.16 | 6.95 | 1.80 | 4.63 | 1.92 | 4.81 | 2.01 | 4.97
| (1.9/13.8) | (4.6/27.5) | (1.7/2.9) | (4.3/7.1) | (1.8/3.1) | (4.5/7.9) | (1.9/3.3) | (4.6/8.5)

###  4.2 Results on offline CTC/attention-based system

Table 1 shows the results of the offline CTC/attention-based systems obtained
on the Librispeech-960 dataset for different bias list sizes $N$. With a bias
list size $N>0$, the proposed method improves the B-WER considerably despite a
small increase in the U-WER, which results in a substantial improvement in the
overall WER. While the B-WER and U-WER tend to deteriorate with larger $N$,
the proposed method remains superior to other DB techniques across all bias
list sizes. Although the proposed method does not perform as well as the
baseline when $N=0$, this is not considered a particularly limiting drawback
considering that users typically register important keywords in the bias list.

###  4.3 Analysis of proposed bias token

Figure 4 shows an example of the cumulative log probability described in Eq.
(15), where the blue and red lines indicate the results obtained with and
without the proposed bias tokens. Without the use of dynamic tokens, the model
struggles to capture the relationships between the subwords, which results in
significantly lower scores for each subword. Conversely, the proposed method
assigns high scores to the bias token ($<$Nelly$>$), resulting in a better
B-WER (Table 1). Interestingly, the log probabilities before and after the
bias token ($<$fresh$>$ and $<$is$>$) remain stable, which indicates that the
proposed method maintains the context successfully throughout the entire token
sequence.

Figure 4: Example of cumulative log probability

Figure 5: Effect of the bias weight $\mu$.

###  4.4 Effect of the bias weight during inference

Figure 5 shows the effect of the bias weights $\mu$ (Section 3.4) on the WER,
U-WER, and B-WER results when the bias list size $N$ = 2000. Increasing the
bias weights $\mu$ improves the B-WER but deteriorates U-WER due to its
tendency to overbias. Under this experimental condition, there is a slight
tendency to overbias when no bias weights are introduced; thus, when
$\mu=0.8$, over-bias can be suppressed. The degree to which the model should
be biased depends on the target user domain; thus, we believe that this
mechanism to adjust the bias weights easily during inference is effective.

###  4.5 Validation on Japanese dataset

We also validate the proposed method on a Japanese dataset containing the
Corpus of Spontaneous Japanese (581 h) [41], 181 hours of Japanese speech in
the database developed by the Advanced Telecommunications Research Institute
International [42], and 93 hours of Japanese speech data, including meeting
and morning assembly scenarios with the attention-based system described in
Section 4.1. Table 2 shows the results obtained when $N$ = 203 technical terms
are registered in the bias list $\bm{B}$. The proposed method improves the
B-CER significantly with a slight degradation in U-CER, thereby resulting in
the best overall CER.

Table 2: Experimental results obtained on Japanese dataset.

Model | CER | U-CER | B-CER
---|---|---|---
Baseline [8] | 9.85 | 8.17 | 21.76
BPB beam search [22] | 9.67 | 9.20 | 13.16
Proposed | 9.03 | 8.93 | 9.73

###  4.6 Validation on streaming RNN-T-based system

Table 3 shows the results of the streaming RNN-T-based system with a bias list
of size $N$ = 100. Consistent with the results from the offline CTC/attention-
based system (Tabel 1), the proposed method significantly improves the B-WER
compared to the conventional DB method without relying on additional
information such as phonemes, while maintaining a comparable overall WER to
conventional DB approaches (A1-2 vs. A3). Furthermore, the conventional DB
method [25, 26] significantly improves the B-WER using additional text data
(A1-2 vs. B1-2), whereas the proposed method achieves similar B-WER without
using such additional data (A3 vs. B1-2).

Table 3: WER results of streaming RNN-T-based systems on Librispeech-960 test-
clean.

ID | Model | Text data | WER/B-WER
---|---|---|---
A0 | Baseline (RNN-T) | None | 3.7 / 14.1
A1 | Trie-based DB + WFST [25] | None | 2.8 / 7.4
A2 | Phoneme-based DB [26] | None |  2.6 / 6.8
A3 | Proposed | None |  2.6 / 3.1
B1 | A1 + NNLM* [25] | *810M words |  2.0 / 5.7
B2 | A2 + USTR* + WFST [26] | *810M words | 2.1 / 2.0

##  5 Conclusion

In this paper, we present a simple but effective deep biasing method that
introduces a dynamic vocabulary where the phrase-level dynamic bias tokens
with the single entities. To enhance the usability of the proposed method, we
introduce a bias weight to adjust the bias intensity during inference. The
experimental results obtained by applying the proposed method to an offline
CTC/attention-based system and a streaming RNN-T-based system demonstrate that
it significantly improves bias phrase recognition on English and Japanese
datasets.

## References

  * [1] R. Prabhavalkar _et al._ , ``End-to-end speech recognition: A survey,'' _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 32, pp. 325–351, 2023.
  * [2] A. Graves, S. Fernández, F. Gomez, and J. Schmidhuber, ``Connectionist temporal classification: Labelling unsegmented sequence data with recurrent neural networks,'' in _Proc. ICML_ , 2006, pp. 369–376.
  * [3] A. Graves and N. Jaitly, ``Towards end-to-end speech recognition with recurrent neural networks,'' in _Proc. ICML_ , 2014, pp. 1764–1772.
  * [4] A. Graves, ``Sequence transduction with recurrent neural networks,'' in _Proc. ICML_ , 2012.
  * [5] A. Gulati _et al._ , ``Conformer: Convolution-augmented transformer for speech recognition,'' in _Proc. Interspeech_ , 2020, pp. 5036–5040.
  * [6] J. K. Chorowski, D. Bahdanau, D. Serdyuk, K. Cho, and Y. Bengio, ``Attention-based models for speech recognition,'' _Advances in neural information processing systems_ , vol. 28, 2015.
  * [7] W. Chan, N. Jaitly, Q. Le, and O. Vinyals, ``Listen, attend and spell: A neural network for large vocabulary conversational speech recognition,'' in _Proc. ICASSP_ , 2016, pp. 4960–4964.
  * [8] S. Watanabe, T. Hori, S. Kim, J. R. Hershey, and T. Hayashi, ``Hybrid ctc/attention architecture for end-to-end speech recognition,'' _IEEE Journal of Selected Topics in Signal Processing_ , vol. 11, no. 8, pp. 1240–1253, 2017.
  * [9] K. Hu, T. N. Sainath, R. Pang, and R. Prabhavalkar, ``Deliberation model based two-pass end-to-end speech recognition,'' in _Proc. ICASSP_ , 2020, pp. 7799–7803.
  * [10] Y. Sudo, M. Shakeel, B. Yan, J. Shi, and S. Watanabe, ``4D ASR: Joint modeling of CTC, attention, transducer, and mask-predict decoders,'' in _Proc. Interspeech_ , 2023, pp. 3312–3316.
  * [11] Y. Peng _et al._ , ``Reproducing whisper-style training using an open-source toolkit and publicly available data,'' in _Proc. ASRU_ , 2023, pp. 1–8.
  * [12] R. Huang, O. Abdel-Hamid, X. Li, and G. Evermann, ``Class lm and word mapping for contextual biasing in end-to-end asr,'' in _Proc. Interspeech_ , 2020, pp. 4348–4351.
  * [13] I. Williams, A. Kannan, P. Aleksic, D. Rybach, and T. Sainath, ``Contextual speech recognition in end-to-end neural network systems using beam search,'' in _Proc. Interspeech_ , 2018.
  * [14] A. Kannan _et al._ , ``An analysis of incorporating an external language model into a sequence-to-sequence model,'' in _Proc. ICASSP_ , 2018, pp. 5824–5828.
  * [15] A. Sriram, H. Jun, S. Satheesh, and A. Coates, ``Cold fusion: training seq2seq models together with language models,'' in _Proc. Interspeech_ , 2018, pp. 387–391.
  * [16] G. Pundak, T. N. Sainath, R. Prabhavalkar, A. Kannan, and D. Zhao, ``Deep context: End-to-end contextual speech recognition,'' in _Proc. SLT_ , 2018, pp. 418–425.
  * [17] M. Jain, G. Keren, J. Mahadeokar, and Y. Saraf, ``Contextual rnn-t for open domain asr,'' in _Proc. Interspeech_ , 2020, pp. 11–15.
  * [18] K. Huang _et al._ , ``Contextualized End-to-End Speech Recognition with Contextual Phrase Prediction Network,'' in _Proc. Interspeech_ , 2023, pp. 4933–4937.
  * [19] M. Han _et al._ , ``Improving end-to-end contextual speech recognition with fine-grained contextual knowledge selection,'' in _Proc. ICASSP_ , 2022, pp. 491–495.
  * [20] C. Huber, J. Hussain, S. Stüker, and A. Waibel, ``Instant one-shot word-learning for context-specific neural sequence-to-sequence speech recognition,'' in _Proc. ASRU_ , 2021, pp. 1–7.
  * [21] S. Zhou, Z. Li, Y. Hong, M. Zhang, Z. Wang, and B. Huai, ``Copyne: Better contextual asr by copying named entities,'' _arXiv preprint arXiv:2305.12839_ , 2023.
  * [22] Y. Sudo, M. Shakeel, Y. Fukumoto, Y. Peng, and S. Watanabe, ``Contextualized automatic speech recognition with attention-based bias phrase boosted beam search,'' in _Proc. ICASSP_ , 2024, pp. 10 896–10 900.
  * [23] S. Dingliwal, M. Sunkara, S. Ronanki, J. Farris, K. Kirchhoff, and S. Bodapati, ``Personalization of ctc speech recognition models,'' in _Proc. SLT_ , 2023, pp. 302–309.
  * [24] G. Sun, C. Zhang, and P. C. Woodland, ``Tree-constrained pointer generator for end-to-end contextual speech recognition,'' in _Proc. ASRU_ , 2021, pp. 780–787.
  * [25] D. Le, Jain _et al._ , ``Contextualized streaming end-to-end speech recognition with trie-based deep biasing and shallow fusion,'' in _Proc. Interspeech_ , 2021, pp. 1772–1776.
  * [26] J. Qiu, L. Huang, B. Li, J. Zhang, L. Lu, and Z. Ma, ``Improving large-scale deep biasing with phoneme features and text-only data in streaming transducer,'' in _Proc. ASRU_ , 2023, pp. 1–8.
  * [27] A. Bruguier, R. Prabhavalkar, G. Pundak, and T. N. Sainath, ``Phoebe: Pronunciation-aware contextualization for end-to-end speech recognition,'' in _Proc. ICASSP_ , 2019, pp. 6171–6175.
  * [28] Z. Chen, M. Jain, Y. Wang, M. L. Seltzer, and C. Fuegen, ``Joint grapheme and phoneme embeddings for contextual end-to-end asr.'' in _Proc. Interspeech_ , 2019, pp. 3490–3494.
  * [29] H. Futami, E. Tsunoo, Y. Kashiwagi, H. Ogawa, S. Arora, and S. Watanabe, ``Phoneme-aware encoding for prefix-tree-based contextual asr,'' in _Proc. ICASSP_ , 2024.
  * [30] Y. Sudo, K. Hata, and K. Nakadai, ``Retraining-free customized asr for enharmonic words based on a named-entity-aware model and phoneme similarity estimation,'' in _Proc. Interspeech_ , 2023, pp. 3312–3316.
  * [31] X. Wang _et al._ , ``Towards contextual spelling correction for customization of end-to-end speech recognition systems,'' _IEEE Trans. Audio, Speech, Lang. Process._ , vol. 30, pp. 3089–3097, 2022.
  * [32] X. Wang, Y. Liu, J. Li, and S. Zhao, ``Improving contextual spelling correction by external acoustics attention and semantic aware data augmentation,'' in _Proc. ICASSP_ , 2023, pp. 1–5.
  * [33] Y. Sudo, S. Muhammad, Y. Peng, and S. Watanabe, ``Time-synchronous one-pass Beam Search for Parallel Online and Offline Transducers with Dynamic Block Training,'' in _Proc. Interspeech_ , 2023, pp. 4479–4483.
  * [34] E. Tsunoo, H. Futami, Y. Kashiwagi, S. Arora, and S. Watanabe, ``Integration of Frame- and Label-synchronous Beam Search for Streaming Encoder-decoder Speech Recognition,'' in _Proc. Interspeech_ , 2023, pp. 1369–1373.
  * [35] A. Radford, J. W. Kim, T. Xu, G. Brockman, C. McLeavey, and I. Sutskever, ``Robust speech recognition via large-scale weak supervision,'' in _Proc. ICML_ , 2023, pp. 28 492–28 518.
  * [36] Y. Peng _et al._ , ``Owsm v3.1: Better and faster open whisper-style speech models based on e-branchformer,'' _arXiv preprint arXiv:2401.16658_ , 2024.
  * [37] Y. Peng, Y. Sudo, M. Shakeel, and S. Watanabe, ``Owsm-ctc: An open encoder-only speech foundation model for speech recognition, translation, and language identification,'' _arXiv preprint arXiv:2402.12654_ , 2024.
  * [38] E. Tsunoo, Y. Kashiwagi, T. Kumakura, and S. Watanabe, ``Transformer ASR with contextual block processing,'' in _Proc. ASRU_ , 2019, pp. 427–433.
  * [39] V. Panayotov _et al._ , ``Librispeech: an asr corpus based on public domain audio books,'' in _Proc. ICASSP_ , 2015, pp. 5206–5210.
  * [40] S. Watanabe _et al._ , ``ESPnet: End-to-end speech processing toolkit,'' in _Proc. Interspeech_ , 2018, pp. 2207–2211.
  * [41] K. Maekawa, ``Corpus of spontaneous Japanese: Its design and evaluation,'' in _ISCA & IEEE Workshop on Spontaneous Speech Processing and Recognition_, 2003.
  * [42] A. Kurematsu _et al._ , ``Atr japanese speech database as a tool of speech recognition and synthesis,'' _Speech Communication_ , vol. 9, no. 4, pp. 357–363, 1990.


####  D.3.2 The Practical Generalization Bound of
$\widetilde{\mathcal{F}}_{\beta,\nu}$

The following lemma bridges $(B,d)$-Lipschitz parameterized and covering
number, and Proposition 5 further provides the upper bound of the Gaussian
complexity in Eq.(17).

######  Lemma 5 (Upper bound of Covering Number [78]).

If the score function set $\mathcal{F}$ is $(B,d)$-Lipschitz parameterized,
then the metric entropy defined on the metric space
$(\mathcal{F},d_{\infty,\mathcal{S}})$ is bounded:

$\log\mathfrak{C}(\mathcal{F},\epsilon,d_{\infty,\mathcal{S}})\leq
d\log{\frac{3B}{\epsilon}}$ (35)

######  Proposition 5.

Given the score function $f\in\widetilde{\mathcal{F}}$, if
$\widetilde{\mathcal{F}}$ is $(B,d)$-Lipschitz parameterized, then we have:

$\mathfrak{G}_{\tilde{\mathcal{S}}}(\widetilde{\mathcal{F}})\precsim\mathcal{O}\left(\frac{d\log{BnC}}{\sqrt{nC}}\right),$
(36)

and

$\mathfrak{G}_{\tilde{\mathcal{S}}^{\prime}}(\widetilde{\mathcal{F}})\precsim\mathcal{O}\left(\frac{d\log{Bn}}{\sqrt{n}}\right),$
(37)

where $\mathcal{O}$ is the complexity.

###### Proof.

Applying Lem. 5 to Eq. (19), we have

$\displaystyle\mathfrak{G}_{\tilde{\mathcal{S}}}(\widetilde{\mathcal{F}})$
$\displaystyle\leq\inf_{\alpha\geq
0}\left(C_{1}\alpha+\frac{C_{2}}{\sqrt{nC}}\int_{\alpha}^{1}d\log{\frac{3B}{\epsilon}}d\epsilon\right)$
$\displaystyle\leq\inf_{\alpha\geq
0}\left(C_{1}\alpha+\frac{C_{2}}{\sqrt{nC}}(1-\alpha)d\log{\frac{3B}{\alpha}}\right)$

The last equation holds since $d\log{\frac{3B}{\epsilon}}$ is non-increasing
w.r.t. $\epsilon$. By setting $\alpha=\frac{1}{\sqrt{nC}}$, we have

$\mathfrak{G}_{\tilde{\mathcal{S}}}(\widetilde{\mathcal{F}})\leq\frac{1}{\sqrt{nC}}\left[C_{1}+C_{2}d(1-\frac{1}{\sqrt{nC}})\log{3B\sqrt{nC}}\right].$
(38)

Similarly, we have

$\mathfrak{G}_{\mathcal{S}}(\widetilde{\mathcal{F}})\leq\frac{1}{\sqrt{n}}\left[C_{1}+C_{2}d(1-\frac{1}{\sqrt{n}})\log{3B\sqrt{n}}\right].$
(39)

Then the proof ends. ∎

Finally, the practical generalization bound of the convolutional neural
networks we discuss is clear: See 7

###### Proof.

The proof ends by Proposition 2, Proposition 5 and Corollary 3. ∎