# A Nonlocal Extension of Dispersion
Analysis for Closures in Reactive Flows

Omkar B. Shendea,∗, Ali Mania
aDepartment of Mechanical Engineering, Stanford University

###### Abstract

Nonlocal algebraic closure models associated with both unresolved advective
transport and nonlinear reaction terms in a Reynolds-averaged Navier-Stokes
context are presented in this work. In particular, a system of species subject
to binary reactions and transport by advection and diffusion are examined by
expanding upon analysis originally developed for binary reactions in the
context of Taylor dispersion of scalars. This work extends model forms from
weakly-nonlinear extensions of that dispersion theory and the role of
nonlocality in the presence of reactions is studied and captured by analytic
algebraic expressions. These expressions can be incorporated into an eddy
diffusivity matrix that explicitly capture the influence of chemical kinetics
on the closure operators. Furthermore, we demonstrate that the model form
derived in a laminar context directly translates to an analogous setup in
homogeneous isotropic turbulence. We show that this framework improves
prediction of mean quantities compared to previous results.

Keywords: Nonlocal Models; RANS Closures; Scalar Transport

##  1 Introduction

In considering the temporal evolution of a binary reaction system with two
independent scalars with corresponding concentrations labeled as $C_{1}$ and
$C_{2}$ in the presence of an imposed divergence-free velocity field, $U$, the
governing equations for the scalars can be formulated as

$\frac{\partial C_{i}}{\partial
t}+\nabla\cdot\left(UC_{i}\right)=D_{m}\nabla^{2}C_{i}-AC_{i}C_{j\neq i},$ (1)

where $A$ is a reaction coefficient and $D_{m}$ is a scalar diffusivity that
is identical for both indexed scalars. In examining the evolution of the
concentrations of passive scalars in such an incompressible flow, a one-way
coupling between the fluid momentum and the scalars can be assumed in writing
the transport equations given by Equation (1). This assumption is physically
realizable in the dilute species limit, where the effects of heat release from
the reaction are so small that neither the fluid density nor the reaction
kinetics are significantly affected.

While this form of the equation incorporates all effects present in the
system, in practice, the Reynolds-averaged analogue of Equation (1), written
as

$\frac{\partial\overline{C_{i}}}{\partial
t}+\nabla\cdot\left(\overline{U}_{\,}\overline{C}_{i}\right)+\nabla\cdot\left(\overline{U^{\prime}C_{i}^{\prime}}\right)=D_{m}\nabla^{2}\overline{C_{i}}\\\
-A\overline{C}_{i}\overline{C}_{j\neq
i}-A\overline{C^{\prime}_{i}C^{\prime}_{j\neq i}},$ (2)

is far less expensive to solve. This equation exists in the Reynolds-averaged
Navier-Stokes (RANS) space, and we will therefore refer to it as a RANS scalar
equation, even though we do not solve the RANS momentum equation herein.
Primed terms represent fluctuations about the averaged quantities, which are
denoted by over-bars. Finding models for the unclosed terms in Equation (2),
given by the scalar flux term, $\overline{U^{\prime}C_{i}^{\prime}}$, and a
reaction source term, $\overline{C^{\prime}_{i}C^{\prime}_{j\neq i}}$, is
essential to solving the full RANS problems and enabling faster predictions
for computational problem in combustion and other fields where similar
transport equations are solved, like atmospheric pollution.

In this work, a model form first proposed by [1] derived using a weakly-
nonlinear extension of dispersion analysis will be briefly introduced. In what
follows, we will first develop a model problem based on laminar flows that can
capture dynamics neglected in that work for a linear reaction setup that will
capture the effects of nonlocality. The resulting nonlocal model form will be
generalized using the reduced-order model (ROM) proposed in [1] such that it
can be applied to a binary reaction problem in a general turbulent flow. The
outcome of such analysis is an algebraic model in the RANS context for
nonlocal closures that clearly illustrates the role of chemical kinetics and
nonreactive flow parameters on the magnitude of transport and reaction
closures.

##  2 Local Formulation

In this study, we build off of the work of [1] to develop a model in the RANS
context that can provide algebraic closures to the scalar evolution equations
for a binary reactant setup. This is done in the spirit of dispersion
analysis, as first demonstrated by Taylor in [2, 3] and since extended by many
others, cf. [4, 5, 6, 7]. By considering the case of a parallel flow as an
analytical prototype, [1] extends Taylor’s analysis to write local model forms
for the case of dispersion of scalars undergoing binary reactions, which, in
particular, lead to scalar flux closures that can be written as

$\overline{u_{i}^{\prime}C_{j}^{\prime}}=-[D_{eff}]\left[\frac{\partial\overline{C}}{\partial
x_{i}}\right],$ (3)

where $D_{eff}$ is some generalized effective eddy diffusivity matrix and
gradients of both reactants play a role. If considering a single scalar and
using a standard gradient diffusion hypothesis, one would simply conclude that
$D_{eff}=D^{0}$, where $D^{0}$ is a turbulent eddy diffusivity that is set
purely by the flow. In the presence of reacting scalars, however, values of
$D_{eff}$ differ from those observed in non-reactive cases and these
differences can be explained though the use of a weakly-nonlinear extension of
dispersion analysis. In particular, the effective eddy diffusivity involves
the the superposition of two scalar gradients multiplied by a prefactor
dependent on both scalar fields. The equations for the scalar flux are written
as

$\begin{bmatrix}\overline{u^{\prime}C_{1}^{\prime}}\\\
\overline{u^{\prime}C_{2}^{\prime}}\end{bmatrix}=-\begin{bmatrix}D_{11}&D_{12}\\\
D_{21}&D_{22}\end{bmatrix}\frac{\partial}{\partial
x}\begin{bmatrix}\overline{C_{1}}\\\ \overline{C_{2}}\end{bmatrix},$ (4)

where $D_{kl}$ represents the diffusivity coefficient associated with the flux
of the $k$-th species due to a gradient in the species. These coefficients are
written explicitly as

$\displaystyle\begin{split}&D_{11}=\frac{D^{0}(1+A\overline{C_{1}}\tau_{mix})}{1+A\tau_{mix}\overline{C_{1}}+A\tau_{mix}\overline{C_{2}}}\\\
&D_{12}=-\frac{D^{0}A\tau_{mix}\overline{C_{1}}}{1+A\tau_{mix}\overline{C_{1}}+A\tau_{mix}\overline{C_{2}}}\\\
&D_{21}=-\frac{D^{0}A\tau_{mix}\overline{C_{2}}}{1+A\tau_{mix}\overline{C_{1}}+A\tau_{mix}\overline{C_{2}}}\\\
&D_{22}=\frac{D^{0}(1+A\overline{C_{2}}\tau_{mix})}{1+A\tau_{mix}\overline{C_{1}}+A\tau_{mix}\overline{C_{2}}}.\end{split}$
(5)

Here, $\tau_{mix}$ denotes the characteristic mixing time associated with a
particular turbulent flow. These two quantities, the mixing time and eddy
diffusivity, are measured from a non-reactive flow using methods such as the
macroscopic forcing method (MFM) of [6, 8].

The most novel implication of the derivation approach used to obtain this
model, however, is that it invokes an implied model for the fluctuating fields
themselves. As a result, it is possible to determine the actual local form of
the fluctuations and to close not only the scalar transport term, but also the
unclosed reaction terms in the standard RANS equations. This closure
expression can be written as

$\overline{AC_{1}^{\prime}C_{2}^{\prime}}=\frac{A}{u_{rms}^{2}}\\\
\left(D_{11}\frac{\partial\overline{C_{1}}}{\partial
x}+D_{12}\frac{\partial\overline{C_{2}}}{\partial
x}\right)\left(D_{21}\frac{\partial\overline{C_{1}}}{\partial
x}+D_{22}\frac{\partial\overline{C_{2}}}{\partial x}\right),$ (6)

where $u_{rms}$ is the measured root-mean-squared velocity field of the
underlying flow. As such, we can fully close the scalar transport equation for
both scalars using a single framework.

The basic weakly nonlinear model offers an alternative approach to
perturbation expansions based on directly linearizing the equations performed
by works like [9]. It recovers the scaling relationship between diffusivity
and Damköhler number derived in other numerical works and experiments, as in
[10, 11, 12, 7], and the model based on [13] explicated in [14].

This provides a description of the local model that solves the reduced-degree-
of-freedom system resulting from the RANS equations. The three algebraic
components describe an entire model form, and require only information about
the underlying flow, specifically $D^{0}$ and $u_{rms}$.

##  3 The model problem

In order to develop insights into a model form that captures nonlocal closure
terms in a binary mixture, we propose the following illustrative setup that
can be treated semi-analytically. We consider a flow field that is parallel,
steady, and two-dimensional and is contained in a domain that is an elongated
box. We place periodic boundary conditions on the two opposing elongated
sides, as pictured in Figure (1). By elongated, we mean that $L_{1}>>L_{2}$.

In such a two-dimensional domain, we define an averaging operator for the RANS
equations as

$\overline{f(x)}=\frac{1}{L_{2}}\int_{0}^{L_{2}}f(x,y)dy,$ (7)

which differs from the ensemble-averaged version implied earlier. In such a
domain, we assume a periodic structure for the flow-field with no imposed mean
velocity. Therefore, $u(x,y)=\overline{u}+u^{\prime}=U_{0}sin(ky)$, where
$k=2\pi/L_{2}$. At $x=-L1/2$ and $x=L_{1}/2$, the two domain boundaries, we
prescribe Dirichlet conditions for the scalar concentrations, as pictured in
Figure (1).

In addition, we consider a reaction of the form $C_{1}+C_{2}\rightarrow
C_{3}$, which is considered irreversible and fully activated, but is driven by
a finite rate law of mass action. While simpler than a three-dimensional fully
turbulent flow with reacting scalars, this model problem captures the
essential competing physics of mixing, transport, and reactions that govern
mean scalar concentrations in the more realistic case. However, we can now
formulate nonlocal closures with interpretable forms in this sandbox context.

Fig. 1: A schematic cartoon of the model problem, illustrating dimensions and
boundary conditions. The flow is depicted in black as $u^{\prime}(y)$, with a
representation of a reaction front given at the mid-plane. Axial will refer to
the $x$-direction, while spanwise will refer to the $y$-direction

In particular, past work in [1] and work on a similar problem in [6] has shown
that this setup allows for the formulation of solutions that translate
directly from a laminar flow topology to a turbulent context, so we can be
reasonably confident in the utility of these simplifications.

Fig. 2: Plot of macroscopic operator as a function of the wavenumber for
various magnitudes of the multi-physics coefficient with a fit of the proposed
model of effective eddy diffusivity as given by Equation (15). In this case,
$Z=|Z|e^{i\pi/4}$ and unsteady effects are included. Fig. 3: Plot of
macroscopic operator as a function of the wavenumber for various magnitudes of
the multi-physics coefficient with a fit of the proposed model of effective
eddy diffusivity as given by Equation (15). In this case, $Z$ is purely real
and we are in the steady limit. The imaginary portion in both model and MFM
results is zero.

###  3.1 Nonlocal formulation

To derive nonlocal insight, let us examine this problem in the context of a
linear reaction, as studied in works like [13]. Here, we are following the
basic steps of the MFM procedure that is fully explained in [6] and complete
details can be found in that reference.

For this particular setup, the concentration of one of the scalars, $C_{2}$,
is held constant in the entire domain at unity and does not evolve, while the
maximum value of $C_{1}$, at the left side of the domain, is held to be at
least an order of magnitude less than that of $C_{2}$. As $C_{2}$ is a
constant, we will denote the reaction rate as $A_{L}=AC_{2}$, denoting the
product of a standard binary species reaction coefficient times the
concentration of the constant species. The fluctuations for $C_{1}$ are now
governed by a transport equation as

$\frac{\partial C_{1}^{\prime}}{\partial
t}+\frac{\partial(u^{\prime}C_{1}^{\prime})^{\prime}}{\partial
x_{1}}+u^{\prime}(y)\frac{\partial\overline{C_{1}}}{\partial
x_{1}}=D_{m}\frac{\partial^{2}C_{1}^{\prime}}{\partial x_{i}\partial x_{i}}\\\
-A_{L}C_{1}^{\prime},$ (8)

where $C_{1}=\overline{C_{1}}+C_{1}^{\prime}$ and so on. This equation is
derived by finding the full transport equation for $C_{1}$ and subtracting
from it the evolution equation for $\overline{C_{1}}$. Note that this equation
is linear in the scalar concentration, and so therefore can be described as
evolving via a linear operator acting on the scalar field. If we further
regard the term on the left-hand side that deals with fluctuations of
fluctuations as negligible and add an arbitrary forcing in time and the axial
$x-$ dimension to the equation, we can write the transport equation as

$\frac{\partial C^{\prime}_{1}}{\partial
t}+cos(x_{2})\frac{\partial\overline{C_{1}}}{\partial
x_{1}}=\frac{\partial^{2}C^{\prime}_{1}}{\partial x_{2}^{2}}\\\
+\frac{1}{Pe}\frac{\partial^{2}C^{\prime}_{1}}{\partial
x_{1}^{2}}-A_{L}C^{\prime}_{1}+s(x_{1},t),\\\ $ (9)

where $s=exp(i\omega t+ikx_{1})$ is a unit forcing term added to the right-
hand side and $Pe=uL/D_{m}$ is a Péclet number. We can then also expand the
concentration field as $C^{\prime}_{1}=\hat{c}(\omega,k;x_{2})exp(i\omega
t+ikx_{1})$ as suggested by the imposed velocity field to get a final
governing equation of

$\left[i\omega+\frac{k^{2}}{Pe}+A_{L}-\frac{\partial^{2}}{\partial
x_{2}^{2}}+ikcos(x_{2})\right]\hat{c}=1,$ (10)

where $\omega$ is a temporal frequency and the bracketed terms represent a
linear operator that acts on the scalar field. Note that we can now write this
equation simply as $L\hat{c}(\omega,k;x_{2})=1$, where $L$ denotes the linear
operator.

Following the process of MFM as described in [6], gives us the exact
macroscopic operator. We note that we can also rewrite this equation more
generally as

$\left[Z-\frac{\partial^{2}}{\partial
x_{2}^{2}}+ikcos(x_{2})\right]\hat{c}=1,$ (11)

where $Z$, which we call a “multi-physics coefficient,” is a complex number
that incorporates the effects of axial diffusion, chemical reaction, and
unsteadiness. If we were to apply the over-barred averaging operator to this
equation, only $Z$ and the advective term will remain, allowing us to find the
operator that acts on $\overline{\hat{c}}$.

We can solve this ordinary differential equation directly for each $\omega$
and $k$ and get $\hat{c}$, the scalar magnitude, using appropriate boundary
conditions. Now, after applying the averaging operator to the measured
solution field, we have found the operator that acts on the full Fourier
field. Thus, the macroscopic closure operator that captures purely unresolved
effects in Fourier space can be computed by subtracting the mean field
effects, so that we can write

$\hat{\overline{L^{\prime}}}=1/\overline{\hat{c}}-Z$ (12)

where the overbar represents averaging. This operator is now a function of the
real wavenumber, $k$, and $Z$, which can take complex values.

###  3.2 The closure model form

To showcase the real and imaginary parts of the operator, we will consider $Z$
to take the form of a scalar multiplied by $1+i$, and we can plot the real and
imaginary components of the macroscopic operator as in Figure (2). If we
consider purely the steady limit ($\omega=0$), we can get Figure (3), which
expresses the true steady macroscopic closure operator over a range of $Z$
values. In this limit, we see that the effect of increasing the reaction
coefficient in Equation (9) is to suppress the projection of perturbations to
the macroscopic space.

While we have calculated the values for the macroscopic operator directly, the
practical utility of this process is that one can now fit analytical curves to
the actual calculated MFM operator. We do this by examining the asymptotic
limits of said operator at zero and high $k$ and matching those limits for
intermediate wavenumbers. One analytic expression that matches the limit
values of $\hat{\overline{L^{\prime}}}$ is

$\hat{\overline{L^{\prime}}}=\left((1+Z)^{2}+k^{2}\right)^{1/2}-1-Z,$ (13)

which is a reasonable approximation to the true operator across decades of
wavenumbers. We want to consider the steady limit, and a fit that approximates
the MFM data can be expressed as

$\hat{\overline{L^{\prime}}}=\frac{k^{2}}{\sqrt{4(1+Z)^{2}+k^{2}}},$ (14)

which allows us to consider the steady RANS problem. However, we want to
modify this expression slightly in order to match previously explored model
forms in the local linear limit as well as in the nonlocal but nonreacting
limits.

For a linear reaction problem, [13] explicated a model form that we would like
to match with our new generalized nonlocal operator. In the limit of non-
reactive flow, we would also like to recover the non-reactive scalar transport
model derived in [6]. To do so, we can rewrite the operator derived here for
the linear reaction setup as

$\hat{\overline{L^{\prime}}}=\frac{\hat{\overline{L^{\prime}}}_{nr}}{1+Z\tau_{mix}}=\frac{0.5k^{2}/\sqrt{1+0.25k^{2}}}{1+Z/\sqrt{1+0.25k^{2}}},$
(15)

where the subscript nr denotes the non-reactive macroscopic operator. In the
absence of reaction, the only term that survives is
$\hat{\overline{L^{\prime}}}_{nr}$, which matches the nonreactive but nonlocal
scalar transport model derived in [6].

This model for the eddy diffusivty is plotted in Figures (2-3) against the
true, measured MFM operator. We see that the role of the multi-physics
coefficient is to suppress the magnitude of the macroscopic operator, which is
captured in Equation (15). In the specific steady case where $Z=A_{L}$ and
there is no axial diffusion, we recover the model form suggested by the local
work of [13] in the case of $k=0$.

Furthermore, Equation (15) isolates a nonlocal form for the mixing time scale,

$\tau=\left(1+0.25k^{2}\right)^{-1/2},$ (16)

and this assumed model form is plotted in Figure (4). The colored lines
represent the real measured mixing time as suggested by the middle expression
in Equation (15) from a range of $Z=A_{L}$ values. We see the mixing time is
highly scale-dependent, and is largely agnostic to the presence of a reaction.
While not covered here, the mixing time is sensitive to the steadiness of the
underlying flow.

Fig. 4: Plot of mixing time as a function of the wavenumber for various
magnitudes of the steady-state and real multi-physics coefficient with a fit
of the proposed mixing time.

From Figures (2-3), we can see that as the magnitude of reaction coefficient
or the axial diffusivity increases, the magnitude of the nonlocal behavior of
the closures rolls off.

##  4 Turbulent Connection

We have thus far written a Fourier space operator, but we must translate it
into physical space to draw conclusions about its utility. We start by noting
that each $ik$ term represents a spatial derivative, while each $i\omega$
denotes a time derivative. Now we can use the notation of $D^{0}$ and $l$, an
eddy mixing length, from [8] to transform our operator derived from a laminar
context into one suited for a fully turbulent flow. Such an operator for the
linear reaction problem is given by

$D_{eff}=\frac{D}{1+A_{L}\tau}\\\
=\frac{D^{0}/\sqrt{\mathcal{I}-l^{2}\frac{\partial^{2}}{\partial
x_{1}^{2}}}}{\mathcal{I}+A_{L}\tau_{mix}/\sqrt{\mathcal{I}-l^{2}\frac{\partial^{2}}{\partial
x_{1}^{2}}}},$ (17)

where $\mathcal{I}$ represents an identity operator and which assumes that
$Z=A_{L}$ only, as described in the previous section. We are concerned with
the steady-state limit only here, and at the high $Pe$ limit, where we can
neglect all terms in $Z$ except this scalar reaction coefficient.

It is important to note that while this model appears similar to a fractional
step operator, it is derived analytically by considering the asymptotic limits
of a governing equation solution. In addition, while this operator is
structurally similar to the work of [13], it is far more general and allows
for definitions for mixing time and eddy diffusivity that incorporate
nonlocality.

Fig. 5: Predictions of closures for the binary reaction problem subject to a
turbulent flow. $X=x/L_{1}$ and DNS results are presented with 95% confidence
intervals accounting for statistical fluctuations. Fig. 6: Results of _a
posteriori_ analysis, with DNS quantities indicated with overbars, the
nonlocal formalism by NL ROM, and the standard local model with ROM, with
$X=x/L_{1}$.

Based on the translation of the laminar model to a turbulent context, we can
write that

$D=\frac{D^{0}}{\sqrt{\mathcal{I}-l^{2}\frac{\partial^{2}}{\partial
x_{1}^{2}}}},$ (18)

and

$\tau=\frac{\tau_{mix}}{\sqrt{\mathcal{I}-l^{2}\frac{\partial^{2}}{\partial
x_{1}^{2}}}},$ (19)

using values for the eddy diffusivity, mixing time, and mixing length as
defined and tabulated in [8].

In previous work, [1] derived model forms for closures of the binary reaction
problem involving $\tau_{mix}$ and $D^{0}$. With the observation that we can
write independent expressions for eddy diffusivity and mixing time, an
extension of those parameters to the binary reaction case seems
straightforward by adapting the model forms of Equation (5) for the binary
reaction problem. The newly stated nonlocal eddy diffusivity, $D$, and
nonlocal mixing time, $\tau$, can be used in place of the previously used
local eddy diffusivity, $D^{0}$ and the local mixing time, $\tau_{mix}$. This
means, for example, we can write the first element of the diffusivity matrix
as

$D_{11}=\frac{D(1+A\overline{C_{1}}\tau)}{1+A\tau\overline{C_{1}}+A\tau\overline{C_{2}}}.$
(20)

Having conceived this model form, we seek to test its applicability for the
binary reactant problem. In particular, we replace the steady, parallel flow
examined in the previous section with three dimensional (3D) homogeneous
isotropic turbulence (HIT) in an elongated domain. To generalize the model
problem, the code of [15] in an incompressible mode was adapted for this work
to simulate HIT in a 3D domain of size $(2\pi)^{3}$. The resulting flow field
is periodically extended in the axial $x$-direction to generate a $20\pi\times
2\pi\times 2\pi$ computational domain for scalar transport. The scalar
transport equations are solved with $C_{1}=C_{ref}$ at $x=-L_{1}/2$ and
$C_{2}=C_{ref}$ at $x=L_{1}/2$. This provides a realistic reaction zone in the
middle of the domain, far from the boundaries.

For our flow simulations, we adopted parameters from [8] for a case of
$Re_{\lambda}=26$, which was obtained by setting the HIT box size of $L=2\pi$,
the turbulent forcing parameter following the prescription of [16] equal to
$0.2792$, and the kinematic viscosity to $\nu=0.0263$.

This setup allows adoption of the values for non-reactive eddy diffusivity.
For this specific case, we use values matching those in [1], which fall within
error bounds reported by [8]. These values, chosen to match constant scalar
flux values outside the reaction zone, are
$D^{0}=0.86u_{rms}^{4}/\epsilon=0.96$, $l=1.23u_{rms}^{3}/\epsilon=1.41$, and
$\tau_{mix}=0.86u_{rms}^{2}/\epsilon=1.02$, where $\epsilon=0.79$ is turbulent
kinetic energy dissipation rate and $u_{rms}=0.97$ is the single-component
root-mean-squared velocity.

For scalar transport, we consider molecular diffusivity $D_{m}=0.0263$,
matching $\nu$, $C_{ref}=1$ and a reaction coefficient of $A=100$. This leads
to $Pe\equiv D^{0}/D_{m}=37$ and $Da\equiv A\tau_{mix}C_{ref}=102$.

###  4.1 A priori analysis

First, let us use a priori analysis to examine the ensemble-averaged binary
reaction problem in the defined turbulent context. In Figure (5), we can see
the plotted closures for the case, as well as the mean scalar concentrations.
The nonlocal ROM recovers from some of the errors in capturing the true
transport closures that are incurred by the local ROM model of [1]. In
particular, we see that the nonlocal model performs far better than the local
model for the reaction closure term.

For these specific closure terms, the standard Gradient Diffusion model
appears to match the DNS data more closely than the local ROM and some regimes
of the nonlocal ROM. However, an overall assessment that considers both
reaction and transport closures reveals the advantage of the nonlocal ROM, as
the Gradient Diffusion transport closure model offers no answers to the
equally vital reaction closure question, and the local model does not
accurately match the peak of the DNS reaction closure term. Quantifying this
advantage can be accomplished by performing a posteriori analysis and
comparing one-dimensional model predictions against DNS data.

###  4.2 A posteriori analysis

In this section, having developed confidence in our nonlocal ROM, we can now
solve Equation (2) directly by invoking the closure models heretofore
presented. The problem setup and parameters used in this section are identical
to those for the _a priori_ analysis.

The only exception is that for the RANS equation, instead of Dirichlet
boundary conditions for the scalars, Neumann boundary conditions with slopes
matching the DNS profiles are used. This is done because the DNS develops
axial boundary layers near the Dirichlet conditions due to local outflow
advection. For the DNS, we have ensured these artificial effects do not
pollute the reaction zone results by ensuring the boundaries are far from the
reaction zone. These boundary layers are absent in the RANS case as the mean
velocity is zero. Appropriate matching of RANS solutions to the DNS ones
should consider concentration profiles outside of these artificial boundary
layers. We have done so by matching the slopes of the RANS concentration
profiles to those of the DNS outside of the boundary layers, but far from the
reaction zones.

In Figure (6), we can see the a posteriori results. In particular, we compare
the predicted mean profile of $C_{1}$ to the results derived from the DNS
described in the previous section. As the problem exhibits mirror symmetry in
the ensemble-averaged sense about the midplane, results for $C_{2}$ are
omitted as they provide the same insight.

It is clear that when all closures are considered, the local ROM outperforms
the standard Gradient Diffusion model, while the nonlocal model outperforms
the local model, especially for $X=x/L_{1}\approx 0$, where the mean reaction
zone, the “flame,” exists.

##  5 Conclusions

In this work, we have introduced a framework to extend analysis created to
study the dispersion of a single passive tracer to analysis of system of
species undergoing binary reactions with nonlocal closures. The nonlocality
for the binary reaction problem is derived from solving the corresponding
linear reaction problem exactly to extract scale-dependent measures of eddy
diffusivity and mixing time. The resulting model parameters and their
incorporation into the model form of [1] introduces closures to both advective
flux and reaction terms with nonlinear interactions between the mean state of
all species, even though the primary set of equations being solved undergo a
linearization procedure.

A particularly salient highlight from the a posteriori analysis is that there
is reduced scalar leakage of $C_{1}$ to the right of the flame-zone, where we
expect it to exponentially decay to zero concentration. As this leakage often
misrepresents the mean reaction zone width, it is particularly important to
capture.

It is also salient to reiterate that while the capturing of nonlocality
produces model parameters that appear notionally similar to a fraction power
operator, they are analytically explainable as fitting the high and low
wavenumber limits of diffusion.

A future avenue of further inquiry is applying the proposed ROM to a Large-
Eddy Simulation context as a subgrid-scale model, where mean space variables
are replaced with filtered variables and quantities like $u_{rms}$ are
evaluated at the grid scale. In additon, while this work examines the steady-
state limit of the nonlocal operator, the use of a generalized multi-physics
coefficient in the model derivation allows the definition of an operator
capable of modeling transient behaviour.

Acknowledgements

Support for this work was provided by the National Science Foundation Graduate
Research Fellowship Program under Grant No. 1656518, the Stanford Graduate
Fellowships in Science and Engineering, and National Science Foundation
Extreme Science and Engineering Discovery Environment resources under Grant
No. CTS190057.

## References

  * [1] O. B. Shende, A. Mani, Closures for multi-component reacting flows based on dispersion analysis (2021).  arXiv:2111.10551.
  * [2] G. I. Taylor, Dispersion of soluble matter in solvent flowing slowly through a tube, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 219 (1137) (1953) 186–203.
  * [3] G. I. Taylor, The dispersion of matter in turbulent flow through a pipe, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 223 (1155) (1954) 446–468.  doi:10.1098/rspa.1954.0130.
  * [4] R. Aris, On the dispersion of a solute in a fluid flowing through a tube, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 235 (1200) (1956) 67 – 77.
  * [5] M. J. Lighthill, Initial Development of Diffusion in Poiseuille Flow, IMA Journal of Applied Mathematics 2 (1) (1966) 97–108.  doi:10.1093/imamat/2.1.97.
  * [6] A. Mani, D. Park, Macroscopic forcing method: A tool for turbulence modeling and analysis of closures, Phys. Rev. Fluids 6 (2021) 054607.  doi:10.1103/PhysRevFluids.6.054607.
  * [7] C. J. Prend, G. R. Flierl, K. M. Smith, A. K. Kaminski, Parameterizing eddy transport of biogeochemical tracers, Geophysical Research Letters 48 (21) (2021) e2021GL094405.  doi:https://doi.org/10.1029/2021GL094405.
  * [8] Y. Shirian, A. Mani, Eddy diffusivity in homogeneous isotropic turbulence (2019).  arXiv:1905.08379.
  * [9] I. Battiato, D. Tartakovsky, Applicability regimes for macroscopic models of reactive transport in porous media, Journal of Contaminant Hydrology 120-121 (2011) 18–26.  doi:https://doi.org/10.1016/j.jconhyd.2010.05.005.
  * [10] T. Elperin, N. Kleeorin, M. Liberman, I. Rogachevskii, Turbulent diffusion of chemically reacting gaseous admixtures, Phys. Rev. E 90 (2014) 053001.  doi:10.1103/PhysRevE.90.053001.
  * [11] T. Elperin, N. Kleeorin, M. Liberman, A. N. Lipatnikov, I. Rogachevskii, R. Yu, Turbulent diffusion of chemically reacting flows: Theory and numerical simulations, Phys. Rev. E 96 (2017) 053111.  doi:10.1103/PhysRevE.96.053111.
  * [12] T. Watanabe, Y. Sakai, K. Nagata, O. Terashima, Turbulent schmidt number and eddy diffusivity change with a chemical reaction, Journal of Fluid Mechanics 754 (2014) 98–121.  doi:10.1017/jfm.2014.387.
  * [13] S. Corrsin, The reactant concentration spectrum in turbulent mixing with a first-order reaction, Journal of Fluid Mechanics 11 (3) (1961) 407–416.  doi:10.1017/S0022112061000615.
  * [14] N. Peters, Turbulent Combustion, Cambridge Monographs on Mechanics, Cambridge University Press, 2000.  doi:10.1017/CBO9780511612701.
  * [15] H. Pouransari, M. Mortazavi, A. Mani, Parallel variable-density particle-laden turbulence simulation (2016).  arXiv:1601.05448.
URL<EMAIL_ADDRESS>
  * [16] C. Rosales, C. Meneveau, Linear forcing in numerical simulations of isotropic turbulence: Physical space implementations and convergence properties, Physics of Fluids 17 (9) (2005) 095106.  doi:10.1063/1.2047568.


# Two-Dimensional Melting of two- and three-component mixtures

Yan-Wei Li<EMAIL_ADDRESS>Key Laboratory of Advanced Optoelectronic
Quantum Architecture and Measurement (MOE), School of Physics, Beijing
Institute of Technology, Beijing, 100081, China  Yugui Yao<EMAIL_ADDRESS>Key Laboratory of Advanced Optoelectronic Quantum Architecture and Measurement
(MOE), School of Physics, Beijing Institute of Technology, Beijing, 100081,
China  Massimo Pica Ciamarra<EMAIL_ADDRESS>Division of Physics and
Applied Physics, School of Physical and Mathematical Sciences, Nanyang
Technological University, Singapore 637371, Singapore  CNR–SPIN, Dipartimento
di Scienze Fisiche, Università di Napoli Federico II, I-80126, Napoli, Italy
CNRS@CREATE LTD, 1 Create Way, #08-01 CREATE Tower, Singapore 138602

###### Abstract

We elucidate the interplay between diverse two-dimensional melting pathways
and establish solid/hexatic and hexatic/liquid transition criteria via the
numerical simulations of the melting transition of two- and three-component
mixtures of hard polygons and disks. We show that a mixture’s melting pathway
may differ from its components and demonstrate eutectic mixtures that
crystallize at a higher density than their pure components. Comparing the
melting scenario of many two- and three-component mixtures, we establish
universal melting criteria: the solid and hexatic phases become unstable as
the density of topological defects respectively overcomes $\rho_{d,{\rm
s}}\simeq 0.046$ and $\rho_{d,{\rm h}}\simeq 0.123$.

The nature of the melting of two-dimensional (2d) solids is a long-standing
and fascinating problem in statistical physics Kosterlitz and Thouless (1973);
Halperin and Nelson (1978); Young (1979); Brinkman  _et al._ (1982); Heiney
_et al._ (1983); Strandburg (1988); Barker  _et al._ (1981); Bernard and
Krauth (2011); Thorneywork  _et al._ (2017); Li and Ciamarra (2020); Komatsu
and Tanaka (2015); Li and Pica Ciamarra (2018); Deutschländer  _et al._
(2013); Downs  _et al._ (2021); Deutschländer  _et al._ (2014); Abutbul and
Podolsky (2022). According to the Kosterlitz-Thouless-Halperin-Nelson-Young
(KTHNY) theory Kosterlitz and Thouless (1973); Halperin and Nelson (1978);
Young (1979), 2d solids melt via a continuous solid to hexatic transition
driven by the unbinding of dislocation pairs, followed by a continuous hexatic
to liquid transition driven by the further unbinding of isolated dislocation
into disclinations. Two alternative melting scenarios are possible. In the
mixed case, the solid-to-hexatic transition is continuous, while the hexatic-
to-liquid transition is discontinuous; Finally, the solid transforms into a
liquid via a first-order transition in the discontinuous case. Recent
computational and methodological advances have allowed establishing the
melting scenario of pure two-dimensional substances Bernard and Krauth (2011);
Kapfer and Krauth (2015); Anderson  _et al._ (2017); Deutschländer  _et al._
(2013); Prestipino  _et al._ (2011); Zu  _et al._ (2016); Komatsu and Tanaka
(2015); Li and Pica Ciamarra (2018); Deutschländer  _et al._ (2013); Downs
_et al._ (2021); Deutschländer  _et al._ (2014); Abutbul and Podolsky (2022);
Thorneywork  _et al._ (2017); Li and Ciamarra (2020). It has been clarified,
for instance, that shape controls the melting scenario in hard-particle
systems Bernard and Krauth (2011); Thorneywork  _et al._ (2017); Li and
Ciamarra (2020) with, e.g., hexagons, disks and pentagons following the KTHNY,
mixed and discontinuous melting Anderson  _et al._ (2017); Li and Ciamarra
(2020). In soft repulsive systems, the softness of the inter-particle
interaction influences the melting pathway Kapfer and Krauth (2015), while
attractive interactions promote discontinuous melting Li and Ciamarra (2020).

Phase diagrams are more complex in many-component systems. The relative
fractions become an essential control variable, novel phases may emerge, and
melting may interplay with phase separation. These aspects have been
thoroughly investigated in three-dimensional systems but are almost unexplored
in two dimensions Russo and Wilding (2017); Guo  _et al._ (2021); Sampedro
Ruiz  _et al._ (2019). Can a mixture melt via a scenario different from its
components? Is there a eutectic mixture whose melting temperature is lower
than its components or, for hard particles, a mixture that melts at a density
higher than its components? Do topological defects correlate with the phase
behaviour as in pure systems Digregorio  _et al._ (2022)? Investigating the
melting of mixtures further offers the possibility of validating proposed
melting criteria based on a generalization of Lindemann’s approach Zahn and
Maret (2000); Khrapak (2020); Dillmann  _et al._ (2012) or topological defects
Digregorio  _et al._ (2022); Guo  _et al._ (2021).

In this Letter, we investigate the interplay between diverse melting
mechanisms by tuning the relative composition of two- and three-component
mixtures of hard particles. In hexagons+disks and pentagons+hexagons mixtures,
the melting pathway smoothly crossovers from that of a pure substance to that
of the other. On the contrary, disks+pentagons mixtures may follow the KTHNY
melting scenario, not occurring in pure disks or pentagons. The study of a
large catalogue of two-dimensional systems further demonstrates that the
density of topological defects $\rho_{d}$ controls the stability of the solid
and hexatic phases. As the density decreases, the solid phase becomes unstable
at $\rho_{d,{\rm s}}\simeq 0.046$, and the hexatic one at $\rho_{d,{\rm
h}}\simeq 0.123$ Guo  _et al._ (2021). Our predictions open new routes to
self-assembly in two dimensions Zong and Zhao (2022) and can be experimentally
verified by changing a mixture’s composition.

Figure 1:  Melting of disks+pentagons mixtures. (a) Isothermal equation of
state. Data are vertically shifted for clarity, and the fraction of disks
increases from top to bottom. (b) probability distribution of the local
density, and (c) the translational and (d) the bond-orientational correlation
functions for selected state points marked in (a). Symbols and colours
identify different phases: liquid ( ); hexatic ( ); solid ( ); liquid/solid
coexistence (✩); liquid/hexatic coexistence ( ).

We construct pentagons, hexagons and disks with the same circumscribing circle
of diameter $d$, by lumping together $N_{d}=40$ or $42$ beads equally spaced
along the perimeter, as detailed in the Supplemental Material (SM) SM  and in
Li and Ciamarra (2020). Beads of different particles with separation distance
$r_{\rm b}$ interact via the Weeks, Chandler, and Andersen (WCA) potential
Weeks  _et al._ (1971): $V_{\rm WCA}(r_{\rm
b})=4\epsilon\left[\left(\frac{\sigma}{r_{\rm
b}}\right)^{12}-\left(\frac{\sigma}{r_{\rm b}}\right)^{6}+c\right]$ for
$r_{\rm b}\leq r_{\rm cut}=2^{1/6}\sigma$. The constant $c$ enforces $V_{\rm
WCA}(r_{\rm cut})=0$. We set $\sigma\simeq 0.14d$. We carry out molecular
dynamics simulations in the canonical ensemble, under periodic boundary
conditions within a rectangular box with a side length ratio of $2:\sqrt{3}$
to accommodate the triangular lattice, fixing the temperature via the
Nos$\rm{\acute{e}}$-Hoover thermostat. We fix $T=20\epsilon/k_{\rm b}$ if not
otherwise specified and demonstrate in Fig. S4 SM  that the choice of the
temperature values has no effect on the melting scenario. All simulations are
performed using the graphics processing unit-accelerated GALAMOST package Zhu
_et al._ (2013). We ensure thermal equilibration as detailed in Fig. S2 SM ,
and report results for systems of $N=20521$ particles, where finite size
effects are negligible as detailed in Fig. S3 SM .

Most previous investigations of the equilibrium phase diagram of two-
dimensional mixtures of particles interacting via purely repulsive forces
considered the specific case of hard disks. These studies focused on
identifying the many possible crystalline phases determined by the size ratio
and the relative fraction, e.g. Likos and Henley (1993); Uche  _et al._
(2004); Perera-Burgos  _et al._ (2016); Fayen  _et al._ (2020), and on the
stability of the hard-disk melting scenario Russo and Wilding (2017); Guo  _et
al._ (2021); Sampedro Ruiz  _et al._ (2019). In these systems, geometric
frustration stemming from the size disparity prevents the triangular lattice
from being the ground state. We depart from these previous works as we
consider mixtures of particles with similar sizes but different shapes. By
focusing on mixtures of particles with similar sizes, all crystallizing in the
triangular lattice with the same lattice constant, we avoid phase separation
Prajwal and Escobedo (2021), glass formation Speedy (1999) and the emergence
of complex crystalline phases on increasing the density. Yet, the shapes we
consider melt via different pathways. Hence, while previous works considered
the competition between different possible ground states, we explore the
competition between alternative crystallization pathways to the same ground
state.

To determine the melting scenario, we first consider if the isothermal
equation of state (EOS) possesses a Mayer-Wood loop indicating the presence of
a first-order transition Mayer and Wood (1965). If so, we fit the EOS to a
fifth-order polynomial to identify the coexistence boundaries and further
investigate the local density distribution in the coexisting region 111We
associate to each particle a local density defined as $\phi_{\rm
loc}(\vec{r}_{i})=\frac{\sum_{j=1}^{N}A_{j}H(r_{c}-|\vec{r}_{i}-\vec{r}_{j}|)}{\pi
r_{c}^{2}}$, where $A_{j}$ is the area of particle $j$, $H$ is the Heaviside
step function, and $r_{c}=21d$.. Outside the coexistence region, or in its
absence, we determine the pure phases by investigating the correlation
functions $c(r)$ and $g_{6}(r)$ of the translational and bond-orientational
order Kosterlitz and Thouless (1973); Bernard and Krauth (2011); Li and
Ciamarra (2019). In the liquid phase, both functions decay exponentially. In
the hexatic phase Kosterlitz and Thouless (1973), the system possesses quasi-
long-range bond-orientational order and short-range translational order, so
that $g_{6}(r)\propto r^{-\eta_{6}}$ with $\eta_{6}<1/4$, while $c(r)$ decays
exponentially. In the solid phase, $g_{6}(r)$ does not decay as the bond-
orientational order is extended, while $c(r)\propto r^{-\eta}$ with $\eta\leq
1/3$ decay algebraically as Mermin-Wagner fluctuations makes the translational
order quasi-long-ranged Mermin and Wagner (1966).

Figure 2:  (a): Phase diagram of hexagons+disks, disks+pentagons, and
pentagons+hexagons mixtures, in the relative fraction $f$ and area fraction
$\phi$ plane. (b): the same diagrams are illustrated in the $f$-$\rho_{d}$
plane, with $\rho_{d}$ the density of defective particles that do not have six
immediate neighbours. We evaluate $\rho_{d}$ in $8$ sub-blocks of the system
to estimate its error. The liquid/hexatic, solid/hexatic, and
solid/coexistence boundaries occur at approximately constant defect densities.

Figure 1 illustrates the dependence of the melting scenario of disks+pentagons
mixtures on the fraction of disks, $f_{d}$. The equation of state (a) reveals
a first-order transition at both small and large $f_{d}$ that corresponds to
the discontinuous solid-liquid transition of pure pentagons and the
discontinuous liquid-hexatic transitions of pure disks. When the transition is
discontinuous, the local density distribution for state points inside the
coexistence region (panel b) shows clear bimodality. Consistently with the
absence of segregation phenomena, the coexisting phases have the same disk
fraction $f_{d}$.

The translational and bond-orientational correlation functions illustrated in
Figs. 1(c) and 1(d) clarify that on increasing the disk fraction, the
discontinuous hard-pentagons melting scenario first transforms into the
continuous KTHNY scenario and then into that of hard-disk. Indeed, we find
melting is continuous for disk fractions around $f_{d}\simeq 0.82$. This
result demonstrates that a two-component mixture may melt via a scenario
different from its components. We notice that the crossover from a mixed to
discontinuous melting may not involve an intermediate KTHNY melting scenario,
as observed in attractive hard disks Li and Ciamarra (2020) as the temperature
decreases.

We have also investigated the melting scenario of hexagons+disks and
pentagons+hexagons mixtures and provide details in the SM SM . We summarize
our investigations in Fig. 2 by illustrating how the melting scenario of
hexagons+disks, disks+pentagons and pentagons+hexagons varies with the
fraction of disks, $f_{d}$, and that of hexagons, $f_{h}$. The three diagrams
match at their boundaries corresponding to the melting behaviour of pure
disks, pure pentagons and pure hexagons, which we find to follow the mixed,
discontinuous and KTHNY scenario, consistently with previous results Bernard
and Krauth (2011); Anderson  _et al._ (2017); Li and Ciamarra (2020).

Hexagons+disks mixtures (panel a, left) crossover from the KTHNY to the mixed
scenario as the fraction of disks overcomes $f_{d}\simeq 0.5$, an intermediate
value suggesting that these two shapes do not frustrate each other
considerably. Conversely, pentagons+hexagons mixtures (panel a, right)
crossover from the discontinuous to the KTHNY scenario as the fraction of
hexagons overcomes $f_{h}\simeq 0.9$: adding a small fraction of pentagons
disrupts hexagons’ KTHNY melting scenario. Similarly, a small amount of
pentagons disrupts the melting scenario of hard disks (panel a, middle).
However, in this case, on increasing the pentagons fraction, the melting
scenario first becomes of KTHNY type and then becomes discontinuous, as we
previously noticed.

In hexagons+disks and pentagons+hexagons mixtures, the stability limit
$\phi_{L}(f)$ of the liquid phase weakly depends on $f$. Conversely, in
disks+pentagons mixtures, for $f\simeq 0.5$, the liquid phase is stable up to
$\phi_{L}\simeq 0.715$, a value sensibly higher than that characterizing the
pure phases, $\phi_{L}\simeq 0.70$. This mixture is thus one of the few known
examples of eutectic two-dimensional colloidal systems Toyotama  _et al._
(2016).

Melting in two dimensions has historically been related to the evolution of
the density of topological defects, clusters of particles not having $6$
neighbours. In the KTHNY scenario Kosterlitz and Thouless (1973); Halperin and
Nelson (1978); Young (1979), melting occurs through a continuous solid-hexatic
transition driven by the unbinding of dislocation pairs (clusters of 4
particles with 5+7+5+7 neighbours), followed by a continuous hexatic-liquid
one driven by the further unbinding of isolated dislocation (5+7) into
disclinations (5 or 7). While the observation of extended clusters of
defective particles Digregorio  _et al._ (2022) suggests a rather complex
relation between phase behaviour and defects, Guo et al. Guo  _et al._ (2021)
have found that in systems following the KTHNY scenario, at the hexatic/liquid
boundary, the density of defective particles acquires a universal value,
$\rho_{d,{\rm h}}=0.12$. This number then acts as an upper bound for the
defects’ density in the hexatic phase.

We explore the relationship between defects’ density and phase behaviour by
investigating the phase diagrams in the $\rho_{d}-f$ plane in Fig. 2(b). We
find that the defects’ density at the hexatic/liquid boundary (red circles)
attains the universal value of $\rho_{d,{\rm h}}=0.123\pm 0.006$, in agreement
with the earlier speculation Guo  _et al._ (2021). In addition, we also find
the defects’ density equals $\rho_{d,{\rm s}}=0.046\pm 0.005$ at the
solid/hexatic and solid/coexistence boundaries.

Figure 3:  (a) Ternary phase diagram for mixtures of pentagons, hexagons and
disks. Each point within the triangle corresponds to a composition determined
by the intersections of lines passing through that point and the sides of the
triangle. For example, in point $A$, the fraction of disks, pentagons, and
hexagons are $0.50$, $0.25$ and $0.25$. (b) Phase diagram at fixed hexagon
fraction $f_{h}=0.25$, in the disk fraction, defect density plane. $\rho_{d}$
acquires constant values at the solid/hexatic and hexatic/liquid boundaries.

To test the universality of these findings, we examine ternary mixtures of
pentagons, hexagons and disks. For each composition, we investigate the phase
behaviour as a function of the volume fraction to determine the melting
scenario SM . The ternary phase diagram of Fig. 3(a) summarizes the result of
this investigation. Each point in the diagram represents the melting scenario
observed at a given composition $f_{d},f_{p},f_{h}$. A point’s composition is
read by connecting it to the triangle’s sides, as illustrated with point ‘A’
for which $f_{d}=0.50$, $f_{p}=0.25$ and $f_{h}=0.25$. Hence, the three
corners correspond to the three pure systems and the three sides to the three
binary diagrams, e.g., the bottom left corner corresponds to pure pentagons
and the ‘pentagon’ side to the disks+pentagons mixture.

Pentagon’s discontinuous melting scenario extends over a large part of the
diagram: melting becomes discontinuous as soon as the pentagon’s fraction
overcomes $f_{h}\simeq 0.25$. The KTHNY scenario always separates the mixed
and discontinuous transitions regardless of the hexagon fraction.

In bidisperse systems, the defects’ density controls the stability of the
solid and hexatic phases. We have found the same occurs in the ternary
mixtures. To demonstrate this, we focus here on ternary mixtures with a fixed
fraction of hexagons, $f_{h}=0.25$, and report analogues results for other
$f_{h}$ values in the SM SM . For the studied fraction of disks, $f_{d}\leq
1-f_{h}$, we plot in Fig. 3(b) the value of $\rho_{d}$ at the boundaries
between the different phases the system traverses as the area fraction
increases. Regardless of the $f_{d}$ values, the solid and hexatic phases
become unstable at $\rho_{d,{\rm s}}\simeq 0.046$ and $\rho_{d,{\rm h}}\simeq
0.123$, respectively, as observed in the binary mixtures and the pure phases.

We have studied the two-dimensional melting of mixtures of similar-sized hard
disks, hexagons and pentagons. The size similarity suppresses geometric
frustration and ensures crystallization in the triangular lattice for all
compositions. As such, changes in the relative composition allow us to
investigate the interplay between the discontinuous melting scenario of hard
pentagons, the mixed scenario of hard disks, and the KTHNY scenario of hard
hexagons. Hexagons+disks mixtures smoothly crossover from the KTHNY to the
mixed melting scenario; similarly, hexagons+pentagons mixtures crossover from
the KTHNY to the discontinuous melting scenario. On the contrary, in
disks+pentagons mixtures, the KTHNY melting may result from the competition
between hard disks’ mixed melting scenario and pentagons’ discontinuous one:
two dogs strive for a bone, and the third runs away with it. The exhaustive
investigation of the melting of two- and three-component mixtures further
corroborates speculated Guo  _et al._ (2021) melting criteria and establishes
new ones based on the density of defects. The solid phase becomes unstable for
$\rho_{d}>\rho_{d,{\rm s}}$, the hexatic phase for $\rho_{d}>\rho_{d,{\rm
h}}$. The connection of these criteria with previously speculated Lindemann-
like ones Zahn and Maret (2000); Khrapak (2020); Dillmann  _et al._ (2012) and
the possibility of relating the free energy of the phases to the defect’s
density are fascinating future research directions that could lead to a better
understanding of the mechanisms driving the different melting pathways.

Y.-W.L. acknowledges support of the start-up funding of Beijing Institute of
Technology and support from the National Natural Science Foundation (NSF) of
China (Grants No. 12105012). Y.Y. is supported by the NSF of China (Grants
Nos. 11734003, 12061131002). M.P.C. acknowledges support by the Ministry of
Education, Singapore, under its Academic Research Fund Tier 2
(MOE-T2EP50221-0016) and Tier 1 (MOE-T1-RG56/21).

## References

  * Kosterlitz and Thouless (1973) J. M. Kosterlitz and D. J. Thouless, J. Phys. C 6, 1181 (1973).
  * Halperin and Nelson (1978) B. I. Halperin and D. R. Nelson, Phys. Rev. Lett. 41, 121 (1978).
  * Young (1979) A. P. Young, Phys. Rev. B 19, 1855 (1979).
  * Brinkman  _et al._ (1982) W. F. Brinkman, D. S. Fisher, and D. E. Moncton, Science 217, 693 (1982).
  * Heiney  _et al._ (1983) P. A. Heiney, P. W. Stephens, R. J. Birgeneau, P. M. Horn, and D. E. Moncton, Phys. Rev. B 28, 6416 (1983).
  * Strandburg (1988) K. J. Strandburg, Rev. Mod. Phys. 60, 161 (1988).
  * Barker  _et al._ (1981) J. A. Barker, D. Henderson, and F. F. Abraham, Physica A 106, 226 (1981).
  * Bernard and Krauth (2011) E. P. Bernard and W. Krauth, Phys. Rev. Lett. 107, 155704 (2011).
  * Thorneywork  _et al._ (2017) A. L. Thorneywork, J. L. Abbott, D. G. A. L. Aarts, and R. P. A. Dullens, Phys. Rev. Lett. 118, 158001 (2017).
  * Li and Ciamarra (2020) Y.-W. Li and M. P. Ciamarra, Phys. Rev. Lett. 124, 218002 (2020).
  * Komatsu and Tanaka (2015) Y. Komatsu and H. Tanaka, Phys. Rev. X 5, 031025 (2015).
  * Li and Pica Ciamarra (2018) Y.-W. Li and M. Pica Ciamarra, Phys. Rev. Mater. 2, 045602 (2018).
  * Deutschländer  _et al._ (2013) S. Deutschländer, T. Horn, H. Löwen, G. Maret, and P. Keim, Phys. Rev. Lett. 111, 098301 (2013).
  * Downs  _et al._ (2021) J. Downs, N. Smith, K. Mandadapu, J. P. Garrahan, and M. I. Smith, Phys. Rev. Lett. 127, 268002 (2021).
  * Deutschländer  _et al._ (2014) S. Deutschländer, A. M. Puertas, G. Maret, and P. Keim, Phys. Rev. Lett. 113, 127801 (2014).
  * Abutbul and Podolsky (2022) D. Abutbul and D. Podolsky, Phys. Rev. Lett. 128, 255501 (2022).
  * Kapfer and Krauth (2015) S. C. Kapfer and W. Krauth, Phys. Rev. Lett. 114, 035702 (2015).
  * Anderson  _et al._ (2017) J. A. Anderson, J. Antonaglia, J. A. Millan, M. Engel, and S. C. Glotzer, Phys. Rev. X 7, 021001 (2017).
  * Prestipino  _et al._ (2011) S. Prestipino, F. Saija, and P. V. Giaquinta, Phys. Rev. Lett. 106, 235701 (2011).
  * Zu  _et al._ (2016) M. Zu, J. Liu, H. Tong, and N. Xu, Phys. Rev. Lett. 117, 85702 (2016).
  * Russo and Wilding (2017) J. Russo and N. B. Wilding, Phys. Rev. Lett. 119, 115702 (2017).
  * Guo  _et al._ (2021) J. Guo, Y. Nie, and N. Xu, Soft Matter 17, 3397 (2021).
  * Sampedro Ruiz  _et al._ (2019) P. Sampedro Ruiz, Q.-l. Lei, and R. Ni, Commun. Phys. 2, 70 (2019).
  * Digregorio  _et al._ (2022) P. Digregorio, D. Levis, L. F. Cugliandolo, G. Gonnella, and I. Pagonabarraga, Soft Matter 18, 566 (2022).
  * Zahn and Maret (2000) K. Zahn and G. Maret, Phys. Rev. Lett. 85, 3656 (2000).
  * Khrapak (2020) S. A. Khrapak, Phys. Rev. Res. 2, 012040 (2020).
  * Dillmann  _et al._ (2012) P. Dillmann, G. Maret, and P. Keim, J. Phys. Condens. Matter 24, 464118 (2012).
  * Zong and Zhao (2022) Y. Zong and K. Zhao, Curr. Opin. Solid State Mater. Sci. 26, 101022 (2022).
  * (29) See Supplemental Material at http://… for additional information about numerical details, thermal equilibration, finite size effects, temperature dependence of the melting, the melting of bidisperse and ternary systems, and defects.
  * Weeks  _et al._ (1971) J. D. Weeks, D. Chandler, and H. C. Andersen, J. Chem. Phys. 54, 5237 (1971).
  * Zhu  _et al._ (2013) Y. Zhu, H. Liu, Z. Li, H. Qian, G. Milano, and Z. Lu, J. Comput. Chem. 34, 2197 (2013).
  * Likos and Henley (1993) C. N. Likos and C. L. Henley, Philos. mag., B 68, 85 (1993).
  * Uche  _et al._ (2004) O. Uche, F. Stillinger, and S. Torquato, Phys. A: Stat. Mech. 342, 428 (2004).
  * Perera-Burgos  _et al._ (2016) J. A. Perera-Burgos, J. M. Méndez-Alcaraz, G. Pérez-Ángel, and R. Castañeda-Priego, J. Chem. Phys. 145, 104905 (2016).
  * Fayen  _et al._ (2020) E. Fayen, A. Jagannathan, G. Foffi, and F. Smallenburg, J. Chem. Phys. 152, 204901 (2020).
  * Prajwal and Escobedo (2021) B. P. Prajwal and F. A. Escobedo, Phys. Rev. Mater. 5, 024003 (2021).
  * Speedy (1999) R. J. Speedy, J. Chem. Phys. 110, 4559 (1999).
  * Mayer and Wood (1965) J. E. Mayer and W. W. Wood, J. Chem. Phys. 42, 4268 (1965).
  * Li and Ciamarra (2019) Y.-W. Li and M. P. Ciamarra, Phys. Rev. E 100, 062606 (2019).
  * Mermin and Wagner (1966) N. D. Mermin and H. Wagner, Phys. Rev. Lett. 17, 1133 (1966).
  * Toyotama  _et al._ (2016) A. Toyotama, T. Okuzono, and J. Yamanaka, Sci. Rep. 6, 23292 (2016).


# Improving the fusion of acoustic and text representations in RNN-T

###### Abstract

The recurrent neural network transducer (RNN-T) has recently become the
mainstream end-to-end approach for streaming automatic speech recognition
(ASR). To estimate the output distributions over subword units, RNN-T uses a
fully connected layer as the joint network to fuse the acoustic
representations extracted using the acoustic encoder with the text
representations obtained using the prediction network based on the previous
subword units. In this paper, we propose to use gating, bilinear pooling, and
a combination of them in the joint network to produce more expressive
representations to feed into the output layer. A regularisation method is also
proposed to enable better acoustic encoder training by reducing the gradients
back-propagated into the prediction network at the beginning of RNN-T
training. Experimental results on a multilingual ASR setting for voice search
over nine languages show that the joint use of the proposed methods can result
in 4%–5% relative word error rate reductions with only a few million extra
parameters.

Index Terms—  ASR, RNN-T, fusion, gating, bilinear pooling

##  1 Introduction

In contrast to the traditional modular-based ASR system that consists of an
acoustic model, a language model (LM) and a rule-based decoder, the recent
end-to-end (E2E) approach aims at implementing the ASR process using a single
neural network model. In particular, both attention-based encoder-decoder [1,
2, 3, 4] and recurrent neural network transducer (RNN-T) [5, 6, 7, 8, 9, 10]
methods achieve a coherently integrate acoustic and text information using a
recurrent structure between the previous and current subword units in the
output text sequence in contrast to modular-based ASR systems. Recently, RNN-T
has become more prevalent due to its streaming benefits [11, 12, 13, 14, 15,
16, 17].

RNN-T was first proposed to extend a connectionist temporal classification
(CTC) acoustic encoder [18] with a prediction network serving as an LM.
Acoustic and text representations over phonemes are derived separately from
the acoustic encoder and prediction network, and fused using an addition
followed by a softmax function to produce the final output distributions [5].
Shortly afterward, an improvement was introduced to fuse more compressed
hidden representations using concatenation and a fully connected (FC) layer
with a hyperbolic tangent (tanh) function [6], which is termed as the joint
network. The fused representations are further transformed by the output
layer, which is another FC layer with the softmax function, into the output
distributions. Since then, there has been a research focus to improve the
RNN-T encoder structure, from using the long short-term memory (LSTM) model
[19, 20] to Transformer and Conformer [14, 21], and to improve their streaming
performance and on-device efficiency [11, 22, 23]. More recently, many studies
focus on improving the recognition accuracy on long tail words/phrases and
long-form utterances, which results in the use of extra model components [15,
24], novel prediction network structures and decoding algorithms [7, 25, 26],
alternative subword units to output [20, 11], test-time external language
model integration [27, 28, 29], and synthetic data augmentation and knowledge
distillation methods [11, 30, 31, 32].

Though joint network is the component closest to the output layer and is
important to RNN-T performance, there are only a few studies related to it [6,
28, 33]. In this paper, by viewing the function of the joint network as to
fuse the representations of the acoustic and text modalities, we propose to
improve the joint network implementation using different structures for
information fusion, including gating and a low-rank approximation of bilinear
pooling with shortcut connections and a tanh transform. A better-performing
structure is further proposed by stacking bilinear pooling on top of gating.
Furthermore, since text priors are often easier to learn than complex acoustic
patterns, the prediction network often converges much faster than the acoustic
encoder that causes the joint network overly biased towards the prediction
network. To alleviate this issue, a novel regularisation method is proposed to
penalise the gradients back-propagated into the prediction network in early
training stages, which can improve RNN-T performance without bringing any cost
to both training and test. Experiments were conducted on a large-scale
multilingual ASR setup for voice search, similar to the one used in [34]. As a
result, by jointly using these proposed methods, more than 4% relative word
error rate (WER) reductions were achieved by only increasing a few million
extra parameters.

The rest of the paper is organised as follows: Sec. 2 reviews RNN-T and the
work related to joint network. Sec. 3 gives details of the proposed joint
network structures and the regularisation method. Sec. 4 describes our
experimental setup, followed by the discussions on the results in Sec. 5. We
conclude in Sec. 6.

##  2 Background

###  2.1 RNN Transducer

In the traditional statistical ASR framework, [35], speech is produced and
encoded via a noisy channel and the ASR system is to find the most probable
source text sequence $\mathbf{y}^{*}$ given the acoustic feature sequence
$\mathbf{x}_{1:T}$ of length $T$ observed as the output of the channel. Based
on Bayes' rule, decoding follows the maximum a posteriori rule to search over
each possible hypothesized text sequence $\mathbf{y}$ by

$P(\mathbf{y}|\mathbf{x}_{1:T})\propto
p(\mathbf{x}_{1:T}|\mathbf{y})P(\mathbf{y}),$ (1)

where $p(\mathbf{x}_{1:T}|\mathbf{y})$ is estimated by the acoustic model and
is the likelihood of generating $\mathbf{x}_{1:T}$ through the channel;
$P(\mathbf{y})$ is estimated by an LM, describing the underlying probabilistic
distribution of the source text.

Instead of modelling $p(\mathbf{x}_{1:T}|\mathbf{y})$ and $P(\mathbf{y})$ by
independent models in a modularised system, E2E methods, such as RNN-T,
directly models $P(\mathbf{y}|\mathbf{x}_{1:T})$ by a single model. Let
$\mathbf{y}={y}_{1:U}$ where $U$ is the number of subword units in
$\mathbf{y}$. For a streaming setting without any look ahead frame and time
reduction, $\mathbf{h}^{\text{enc}}_{t}$, the $D^{\text{enc}}$-dimensional
(-dim) acoustic representation extracted by the acoustic encoder at time $t$,
$\mathbf{h}^{\text{pred}}_{u}$, the $D^{\text{pred}}$-dim text representation
of the $u$-th subword unit by the prediction network, and
$\mathbf{h}^{\text{joint}}_{t,u}$, the $D^{\text{joint}}$-dim fused
representation generated by the joint network, are calculated as follows:

$\displaystyle\mathbf{h}^{\text{enc}}_{t}=\text{AcousticEncoder}(\mathbf{x}_{1:t}),$
(2)
$\displaystyle\mathbf{h}^{\text{pred}}_{u}=\text{PredictionNetwork}(y_{1:u-1}),$
(3)
$\displaystyle\mathbf{h}^{\text{joint}}_{t,u}=\text{JointNetwork}(\mathbf{h}^{\text{enc}}_{t},\mathbf{h}^{\text{pred}}_{u}),$
(4) $\displaystyle
P(\hat{y}_{i}=k|y_{0:u-1},\mathbf{x}_{1:t})=\text{Softmax}(\text{W}^{\text{out}}\mathbf{h}^{\text{joint}}_{t,u})|_{k},$
(5)

where $y_{0}$ refers to the special start of sentence symbol; $k$ and
$\mathbf{W}^{\text{out}}$ are the $k$-th node and weights of the output layer.
Regarding a set of subword units $\mathcal{V}$, the symbol that $k$ represents
belongs to $\mathcal{V}\cup\\{\varnothing\\}$, where $\varnothing$ is the
blank symbol indicating no subword is emitted. During training, let
$\hat{\mathbf{y}}=\\{\hat{y}_{1},\hat{y}_{2},\ldots,\hat{y}_{T+U}\\}$ be an
alignment sequence of $\mathbf{y}$ that can be converted into $\mathbf{y}$ by
removing all occurrences of $\varnothing$,
$\mathcal{A}(\mathbf{x}_{1:T},\mathbf{y})$ be the reference lattice including
all possible alignment sequences between $\mathbf{y}$ and $\mathbf{x}_{1:T}$,
$P(\mathbf{y}|\mathbf{x}_{1:T})$ can be computed efficiently with the forward-
backward procedure. There is

$\displaystyle
P(\mathbf{y}|\mathbf{x}_{1:T})=\sum\nolimits_{\hat{\mathbf{y}}\in\mathcal{A}(\mathbf{x}_{1:T},\mathbf{y})}\prod\nolimits^{T+U}_{i=1}P(\hat{y}_{i}|\mathbf{x}_{1:t_{i}},y_{1:u_{i}}),$
(6)

where $t_{i}$ and $u_{i}$ are the values of $t$ and $u$ corresponding to
$\hat{y}_{i}$ in $\hat{\mathbf{y}}$.

In practice, $\text{AcousticEncoder}(\cdot)$ can be Conformer with a fixed
number of look ahead frames and a fixed time reduction rate.
$\text{PredictionNetwork}(\cdot)$ is often a multi-layer LSTM. The joint
network is often defined as an FC layer [6] that

$\displaystyle\mathbf{h}^{\text{joint}}_{t,u}=\tanh(\mathbf{W}^{\text{joint}}_{1}\mathbf{h}^{\text{enc}}_{t}+\mathbf{W}^{\text{joint}}_{2}\mathbf{h}^{\text{pred}}_{u}),$
(7)

where $\mathbf{W}^{\text{joint}}_{1}$ and $\mathbf{W}^{\text{joint}}_{2}$ are
weight matrices. For simplicity, bias is ignored in Eqn. (7) and the rest of
the paper.

If $\mathbf{W}^{\text{joint}}_{2}$ is $\mathbf{0}_{D^{\text{joint}}\times
D^{\text{pred}}}$ and the joint network transforms only the acoustic
representation, apart from their difference in
$\mathcal{A}(\mathbf{x}_{1:T},\mathbf{y})$ [11, 28], RNN-T becomes CTC that
calculates $P(\hat{y}_{i}|\mathbf{x}_{1:t_{i}})$ by making an independence
assumption between any subword units in $\mathbf{y}$. This reveals the
importance of the joint network.

###  2.2 Related work

In Eqn. (7), by enforcing $\mathbf{h}^{\text{enc}}_{t}=\mathbf{0}$, the
prediction network, joint network, and output layer jointly form an LSTM LM
that is often referred to as the internal LM [28]. Studies showed that more
WER reductions can be found from shallow fusion with external LMs by
discounting the internal LM scores at test-time [27, 28, 29]. More recently,
stateless RNN-T has been proposed to truncate the LM history embedded in the
prediction network to $n$-gram [25].

The fusion of representations associated with different modalities plays a key
role in multimodal intelligence [36]: attention, gating, and bilinear pooling
are the most commonly used structures for the purpose. In RNN-T, the standard
joint network implemented based on Eqn. (7) can be viewed as to fuse the
acoustic representation $\mathbf{h}^{\text{enc}}_{t}$ and text representation
$\mathbf{h}^{\text{pred}}_{u}$ using an FC layer. Recently, an alternative
joint network structure is proposed to model the multiplicative interactions
between the two representations [33]:

$\displaystyle\mathbf{h}^{\text{joint}}_{t,u}=\tanh(\mathbf{W}^{\text{joint}}_{1}\mathbf{h}^{\text{enc}}_{t}\odot\mathbf{W}^{\text{joint}}_{2}\mathbf{h}^{\text{pred}}_{u}),$
(8)

where $\odot$ is the Hadamard product.

##  3 Fusing acoustic and text representations

Fusing acoustic and text representations is arguably a difficult task, and the
standard joint network simply uses one FC layer. Sec. 3.1 and 3.2 propose more
complex joint network structures, and Sec. 3.3 proposes to improve the balance
between the two modalities.

###  3.1 Gating mechanism

Gating has been widely used in recurrent and shortcut structures [19, 37],
whose most famous application is LSTM. It allows each element in each
representation vector to be scaled with a different dynamic weight, before
being integrated via vector addition. Specifically,

$\displaystyle\mathbf{g}_{t,u}$
$\displaystyle=\sigma(\mathbf{W}^{\text{gate}}_{1}\mathbf{h}^{\text{enc}}_{t}+\mathbf{W}^{\text{gate}}_{2}\mathbf{h}^{\text{pred}}_{u}),$
(9) $\displaystyle\mathbf{h}^{\text{joint}}_{t,u}$
$\displaystyle=\mathbf{g}_{t,u}\odot\tanh(\mathbf{W}^{\text{joint}}_{1}\mathbf{h}^{\text{enc}}_{t})+(1-\mathbf{g}_{t,u})\odot\tanh(\mathbf{W}^{\text{joint}}_{2}\mathbf{h}^{\text{text}}_{u}),$

where $\mathbf{g}_{t,u}$ is the gating vector, $\sigma(\cdot)$ is sigmoid
function, and $\mathbf{W}^{\text{gate}}_{1}$ and
$\mathbf{W}^{\text{gate}}_{2}$ are weight matrices of the gating layer.
Notably a different gating vector can be computed to replace
$1-\mathbf{g}_{t,u}$. However, we observed worse WERs when using two separate
gating vectors.

###  3.2 Bilinear pooling

Compared to gating, bilinear pooling is a more powerful and expensive method
for fusing multimodal representations [38], which combines two vectors using
the bilinear form (with bias ignored)

$\displaystyle{h}^{\text{joint}}_{t,u,d}=(\mathbf{h}^{\text{enc}}_{t})^{\text{T}}\mathbf{W}^{\text{bi}}_{d}\,\mathbf{h}^{\text{pred}}_{u},$
(10)

where $\mathbf{W}^{\text{bi}}_{d}$ is a $D^{\text{enc}}\times
D^{\text{pred}}$-dim matrix, and ${h}^{\text{joint}}_{t,u,d}$ is the $d$-th
element of $\mathbf{h}^{\text{joint}}_{t,u}$. Considering all elements in
$\mathbf{h}^{\text{joint}}_{t,u}$,
$[\mathbf{W}^{\text{bi}}_{1},\ldots,\mathbf{W}^{\text{bi}}_{D^{\text{joint}}}]$
is a $D^{\text{enc}}\times D^{\text{pred}}\times D^{\text{joint}}$-dim weight
tensor. Fusing $\mathbf{h}^{\text{enc}}_{t}$ and
$\mathbf{h}^{\text{pred}}_{u}$ into $\mathbf{h}^{\text{joint}}_{t,u}$ is
therefore equivalent to perform

$\displaystyle\mathbf{h}^{\text{joint}}_{t,u}=[\text{Vector}(\mathbf{W}^{\text{bi}}_{1}),\ldots,\text{Vector}(\mathbf{W}^{\text{bi}}_{D^{\text{joint}}})]^{\text{T}}\text{Vector}(\mathbf{h}^{\text{enc}}_{t}\otimes\mathbf{h}^{\text{pred}}_{u}),$

where $\text{Vector}(\cdot)$ and $\otimes$ are the vectorisation and outer
product operations. Compared to gating, bilinear pooling first computes the
outer product of the two vectors to capture the multiplicative interactions
between all possible element pairs in a more expressive $D^{\text{enc}}\times
D^{\text{pred}}$-dim space, and then projects it into a $D^{\text{joint}}$-dim
vector space.

To avoid the issues when estimating a high dimensional weight tensor, a low-
rank approximation
$\mathbf{W}^{\text{bi}}_{d}\approx\mathbf{W}^{\text{low}}_{1,d}(\mathbf{W}^{\text{low}}_{2,d})^{\text{T}}$
is suggested [39], where $\mathbf{W}^{\text{low}}_{1,d}$ and
$\mathbf{W}^{\text{low}}_{2,d}$ are $D^{\text{enc}}\times D^{\text{rank}}$-dim
and $D^{\text{pred}}\times D^{\text{rank}}$-dim matrices, and
$D^{\text{rank}}$ is the rank of $\mathbf{W}^{\text{bi}}_{d}$. Therefore Eqn.
(10) can be rewritten as
${h}^{\text{join}}_{t,u,d}\approx(\mathbf{h}^{\text{enc}}_{t})^{\text{T}}\mathbf{W}^{\text{low}}_{1,d}(\mathbf{W}^{\text{low}}_{2,d})^{\text{T}}\mathbf{h}^{\text{pred}}_{u}=\mathbf{1}^{\text{T}}((\mathbf{W}^{\text{low}}_{1,d})^{\text{T}}\mathbf{h}^{\text{enc}}_{t}\odot(\mathbf{W}^{\text{low}}_{2,d})^{\text{T}}\mathbf{h}^{\text{pred}}_{u})$.
It was proposed to tie all $\mathbf{W}^{\text{low}}_{1,d}$ as
$\mathbf{W}^{\text{low}}_{1}$ and all $\mathbf{W}^{\text{low}}_{2,d}$ as
$\mathbf{W}^{\text{low}}_{2}$, and to use a projection matrix
$\mathbf{W}^{\text{proj}}$ to distinguish the elements in
$\mathbf{h}^{\text{joint}}_{t,u}$ [40]. When a tanh function is used to
transform the vectors before the Hadamard product, there is

$\displaystyle\hat{\mathbf{h}}^{\text{joint}}_{t,u}=\mathbf{W}^{\text{proj}}(\tanh((\mathbf{W}^{\text{low}}_{1})^{\text{T}}\mathbf{h}^{\text{enc}}_{t})\odot\tanh((\mathbf{W}^{\text{low}}_{2})^{\text{T}}\mathbf{h}^{\text{pred}}_{u})).$
(11)

We found using shortcut connections [41] and a final tanh transform are
important for bilinear pooling for RNN-T, and thus propose

$\displaystyle{\mathbf{h}}^{\text{joint}}_{t,u}=\tanh(\hat{\mathbf{h}}^{\text{joint}}_{t,u}+\mathbf{W}^{\text{joint}}_{1}\mathbf{h}^{\text{enc}}_{t}+\mathbf{W}^{\text{joint}}_{2}\mathbf{h}^{\text{pred}}_{u}),$
(12)

where $\mathbf{W}^{\text{joint}}_{1}\mathbf{h}^{\text{enc}}_{t}$ and
$\mathbf{W}^{\text{joint}}_{2}\mathbf{h}^{\text{pred}}_{u}$ are the shortcut
connections here.

At last, we propose a stack structure to combine gating and bilinear pooling
to leverage their complementarity. It is implemented by replacing Eqn. (11)
with

$\displaystyle\hat{\mathbf{h}}^{\text{joint}}_{t,u}=\mathbf{W}^{\text{proj}}(\tanh((\mathbf{W}^{\text{low}}_{1})^{\text{T}}\mathbf{h}^{\text{enc}}_{t})\odot\tanh((\mathbf{W}^{\text{low}}_{2})^{\text{T}}\mathbf{h}^{\text{gate}}_{t,u})),$
(13)

where $\mathbf{h}^{\text{gate}}_{t,u}$ refers to the joint representation
computed by Eqn. (9).

###  3.3 Prediction network regularisation

It has been observed that strong and prevalent text priors (e.g. ``bananas are
yellow'') often caused image-text multimodal systems to overfit to the
statistical biases and tendencies, and largely circumvents the need to
understand visual scenes [42]. Similary in RNN-T, since text priors are often
easier to learn than the acoustic patterns, it is possible that the faster
converging speed of the prediction network (than the acoustic encoder) makes
$\mathbf{h}^{\text{pred}}_{u}$ overly weighted in
$\mathbf{h}^{\text{joint}}_{u,t}$. In that situation, the acoustic encoder is
less well trained to handle the audio samples with high internal LM scores.

In order to resolve this issue, we propose a prediction network regularisation
method applied to the beginning of RNN-T training. It is implemented by
recomputing the text representation as

$\displaystyle\mathbf{h}^{\text{pred}}_{u}=\alpha_{m}\,\mathbf{h}^{\text{pred}}_{u}-\text{sg}((\alpha_{m}-1)\,\mathbf{h}^{\text{pred}}_{u}),$
(14)

where $m$ is the index of current training step, $\alpha_{m}$ is a scaling
factor, $\text{sg}(\cdot)$ is the stop gradient function whose input tensor
will have zero gradients. When $0\leqslant\alpha_{m}\leqslant 1$, the value of
$\mathbf{h}^{\text{pred}}_{u}$ will not be changed but its corresponding
gradients that back-propagate into the prediction network will be reduced by a
factor of $\alpha_{m}$. This slows down the convergence of the prediction
network and makes the joint network fuse $\mathbf{h}^{\text{enc}}_{t}$ and
$\mathbf{h}^{\text{pred}}_{u}$ in a more balanced way. In this paper, a piece-
wise linear scheduler is used to control the value of $\alpha_{m}$:

$\displaystyle\alpha_{m}=\left\\{\begin{array}[]{lll}0&\text{if }m<m_{1}\\\
1&\text{else if }m\geqslant m_{2}\\\
(m-m_{1})/(m_{2}-m_{1})&\text{otherwise}\\\ \end{array}\right.,$ (18)

where $m_{1}$ and $m_{2}$ are two pre-defined hyper-parameters. Notably, this
method is different from initialising RNN-T with a pre-trained CTC model even
when $\alpha_{m}=0$, since the prediction network serves as a random but fix-
valued projection, through which RNN-T is still able to obtain $y_{u-1}$. This
links to the stateless RNN-T [25].

##  4 Experimental Setup

###  4.1 Datasets

Experiments were conducted on a dataset with 9 language locales: US English
(en-US), UK English (en-GB), French (fr-FR), Italian (it-IT), Germany (de-DE),
US Spanish (es-US), ES Spanish (es-ES), Taiwan Chinese (zh-TW) and Japanese
(ja-JP). All data are anonymised and hand-transcribed. There are totally
214.2M utterances which correspond to 142.3K hours of speech data collected
from Google’s Voice Search traffic. en-US and en-GB take about 25% and 5% of
the training data, while each of the rest 7 languages takes about 10% of the
data. The SpecAugment method is used to improve ASR robustness against noisy
conditions [43]. The training data is mixed with all languages without using
any language id information. The test sets are kept distinct for each language
with each of them containing 3.3K$\sim$15.4K utterances. The testing
utterances are also sampled from Google’s Voice Search traffic with a maximum
duration constraint of 5.5 second long for each utterance. The test sets have
no overlapping from the training set for evaluation purpose.

###  4.2 Model setup

The 80-dim log-mel filter bank features are used, which are computed using a
32ms frame length and a 10ms shift. Acoustic features from 3 contiguous frames
are stacked and subsampled to form a 240-dim input representation with 30ms
frame rate, which are then transformed using a linear projection to 512-dim
and added with positional embeddings. Twelve Conformer [21] encoder blocks
with 8-head self-attention and a convolution kernel size of 15 are used to
further transform the stacked features. A concatenation operation is performed
after the 3rd block to achieve a time reduction rate of 2, and the resulted
1024-dim vectors are transformed by the 4th Conformer block and then projected
back to 512-dim using another linear transform. Afterwards comes with another
8 Conformer blocks followed by a final linear normalisation layer. These
layers combined together form the RNN-T acoustic encoder. The prediction
network consists of two layers of 2,048-dim LSTM with a 640-dim linear
projection to make $D^{\text{pred}}=$640\. The dimension of the fused
representation $D^{\text{joint}}$ is also set to 640. All models are trained
to predict 16,384 word-piece units [44]. As a result, the final RNN-T baseline
has 110M parameters in the acoustic encoder and 33M parameters in the rest of
the model. All models are trained in Tensorflow using the Lingvo toolkit [45]
on Google’s Tensor Processing Units V3 with a global batch size of 4,096
utterances. Models are optimized using synchronized stochastic gradient
descent based on the Adam optimiser with $\beta_{1}=0.9$ and
$\beta_{2}=0.999$. During test, the models are tested in a fully E2E fashion
without any external LMs.

##  5 Experimental Results

In this section, RNN-T systems with different joint network structures are
first compared at the 200K-th training step, whose details are listed in Table
1. Next, prediction network regularisation with different hyper-parameter
values are compared using S4 at the 500K-th step. Final results are presented
with 800K training steps.

Table 1: Details of systems with different joint network structures. ``#Params'' refers to the number of joint network parameters. ID | Structure | Equations | $D^{\text{joint}}$ | $D^{\text{rank}}$ | #Params
---|---|---|---|---|---
S0 | FC with add. |  (7) | 640 | – | 0.73M
S1 | FC with add. |  (7) | 790 | – | 3.36M
S2 | FC with mul. |  (8) | 640 | – | 0.73M
S3 | Gating |  (9) | 640 | – | 1.47M
S4 | Bilinear |  (11) & (12) | 640 | 640 | 1.88M
S5 | Bilinear |  (11) & (12) | 640 | 1280 | 3.03M
S6 | Combination |  (13) & (12) | 640 | 640 | 3.36M
Table 2: The 200K training step WERs with different joint network structures. ID | de-DE | en-GB | en-US | es-ES | es-US | fr-FR | it-IT | ja-JP | zh-TW | Avg.
---|---|---|---|---|---|---|---|---|---|---
S0 [6] | 15.4 | 7.1 | 7.9 | 8.1 | 8.1 | 13.4 | 9.7 | 13.9 | 6.1 | 9.97
S1 [6] | 15.2 | 7.4 | 7.8 | 8.1 | 7.8 | 13.3 | 9.7 | 14.1 | 6.1 | 9.94
S2 [33] | 15.3 | 7.4 | 7.9 | 8.1 | 7.9 | 13.1 | 9.7 | 14.0 | 6.2 | 9.96
S3 | 15.1 | 7.4 | 7.6 | 7.9 | 8.0 | 12.9 | 9.4 | 13.7 | 6.1 | 9.78
S4 | 15.1 | 7.2 | 7.5 | 7.8 | 7.8 | 13.1 | 9.2 | 13.7 | 5.9 | 9.71
S5 | 15.0 | 7.2 | 7.5 | 7.7 | 7.6 | 13.3 | 9.1 | 13.6 | 6.0 | 9.67
S6 | 14.6 | 7.1 | 7.4 | 7.8 | 7.7 | 13.1 | 9.4 | 13.6 | 6.0 | 9.63
Table 3: Final 800K training step WERs. Superscript ``reg'' means to use regularisation with $m_{1}=$25K and $m_{2}=$200K. ID | de-DE | en-GB | en-US | es-ES | es-US | fr-FR | it-IT | ja-JP | zh-TW | Avg.
---|---|---|---|---|---|---|---|---|---|---
S0 | 14.2 | 6.6 | 6.8 | 7.5 | 7.5 | 13.0 | 8.6 | 12.9 | 5.6 | 9.19
S0${}^{\text{reg}}$ | 13.9 | 6.5 | 6.4 | 7.3 | 7.4 | 12.7 | 8.1 | 12.5 | 5.6 | 8.93
S1 | 14.3 | 6.6 | 7.0 | 7.5 | 7.6 | 13.0 | 8.6 | 12.6 | 5.6 | 9.20
S1${}^{\text{reg}}$ | 14.0 | 6.5 | 6.8 | 7.2 | 7.5 | 12.7 | 8.3 | 12.4 | 5.4 | 8.98
S4 | 13.8 | 6.4 | 6.7 | 7.2 | 7.7 | 12.8 | 8.1 | 12.3 | 5.4 | 8.93
S4${}^{\text{reg}}$ | 13.7 | 6.4 | 6.6 | 7.0 | 7.3 | 12.7 | 7.9 | 12.2 | 5.4 | 8.80
S5 | 13.8 | 6.4 | 6.6 | 7.1 | 7.4 | 12.8 | 8.2 | 12.3 | 5.5 | 8.90
S6 | 13.8 | 6.2 | 6.6 | 7.1 | 7.4 | 12.7 | 8.2 | 12.3 | 5.4 | 8.86
S6${}^{\text{reg}}$ | 13.5 | 6.3 | 6.5 | 7.0 | 7.2 | 12.6 | 8.0 | 12.1 | 5.3 | 8.72

###  5.1 On joint network structures

The results of RNN-T models with different joint networks are presented in
Table 2. First, S0 and S2 result in similar WERs, which validates the finding
in [33]. Next, by increasing $D^{\text{joint}}$ from $640$ (S0) to $790$ (S1),
the averaged WER was only slightly reduced by 0.03%. We also tried to add more
FC layers to the joint network that resulted in worse training loss values and
higher WERs.

Both gating (S3) and bilinear pooling (S4) outperformed FC fusion systems S0
and S1, where S1 has the same amount of parameters as S3 and S4 indicating
that the lower WERs of S3 and S4 do not come from the extra model parameters.
Furthermore, S4 outperformed S3 meaning that bilinear pooling can produce more
expressive joint representations. Another advantage of bilinear pooling is
having the flexibility to control $D^{\text{joint}}$ and $D^{\text{rank}}$
separately. Increasing $D^{\text{rank}}$ from 640 to 1280 allows S5 to have
more parameters and lower WERs without increasing the input size of the big
output layer with 16,384 nodes. At last, by stacking gating and bilinear
pooling in the joint network, S6 leverages the complementarity of both
structures and outperforms any other systems. In particular, with the same
amount of parameters, S6 outperformed S1 by a 0.31% WER reduction. In
conclusion, all of our proposed systems (S3 – S6) improved WERs considerably
by increasing only a small percent of the model parameters (S0 has 144M
parameters as given in Sec. 4.2).

Fig. 1: S4 system WERs with different hyper-parameters for prediction network
regularsation.

###  5.2 On prediction network regularisation

The results of using the prediction network regularsation defined in Eqns.
(14) and (18) with different $m_{1}$ and $m_{2}$ values are shown in Fig 1.
S4, the bilinear pooling joint network with $D^{\text{rank}}=640$, is used in
this section. The best results were found with $m_{1}=$25K and $m_{2}=$200K,
which improved the averaged WER by 0.09%, 0.07%, and 0.13% absolute with 200K,
500K, and 800K steps separately.

###  5.3 Final results

The models are trained for 800K training steps and the final results are shown
in Table 3. Comparing S4, S5, and S6 to S1, better joint network structures
lead to lower WERs without requiring more parameters. The prediction network
regularisation improved averaged WERs by 0.26%, 0.13% and 0.14% for S0, S4,
and S6, although the improvements are not always consistent regarding each
individual language. Compared to baseline system S0, our best-performing
systems S4 and S6 with the regularisation method achieved 4.2% and 5.1%
relative reductions in averaged WER by increasing only 1.15M and 2.63M
parameters.

##  6 Conclusions

By viewing the function of the joint network of RNN-T as to fuse the acoustic
and text representations derived from the acoustic encoder and prediction
network, we propose in this paper to apply the structures widely used for
multimodal representation fusion, including gating and bilinear pooling, to
improve the joint network implementation. These structures are modified and
combined to fit into RNN-T. A novel prediction network regularisation method
is also proposed to make the fusion between acoustic and text representations
more balanced. When evaluated using a large-scale multilingual voice search
setup, 4.2% and 5.1% relative WER reductions were found by using the bilinear
pooling and the combination structure separately along with the
regularisation.

##  7 REFERENCES

## References

  * [1] J. Chorowski, D. Bahdanau, D. Serdyuk, K. Cho, and Y. Bengio,  ``Attention-based models for speech recognition,''  in Proc. NIPS, 2015.
  * [2] L. Lu, X. Zhang, K. Cho, and S. Renals,  ``A study of the recurrent neural network encoder-decoder for large vocabulary speech recognition,''  in Proc. Interspeech, Dresden, 2015.
  * [3] W. Chan, N. Jaitly, Q.V. Le, and O. Vinyals,  ``Listen, attend and spell: A neural network for large vocabulary conversational speech recognition,''  in Proc. ICASSP, Shanghai, 2016.
  * [4] S. Kim, T. Hori, and S. Watanabe,  ``Joint CTC-attention based end-to-end speech recognition using multi-task learning,''  in Proc. ICASSP, New Orleans, 2017.
  * [5] A. Graves,  ``Sequence transduction with recurrent neural networks,''  in Proc. ICML Representation Learning Workshop, Edinburgh, 2012\.
  * [6] A. Graves, A. Mohamed, and G. Hinton,  ``Speech recognition with deep recurrent neural networks,''  in Proc. ICASSP, Vancouver, 2013.
  * [7] R. Prabhavalkar, K. Rao, T. Sainath, B. Li, L. Johnson, and N. Jaitly,  ``A comparison of sequence-to-sequence models for speech recognition,''  in Proc. Interspeech, Stockholm, 2017.
  * [8] E. Battenberg, J. Chen, R. Child, A. Coates, Y. Gaur, Y. Li, H. Liu, S. Satheesh, D. Seetapun, A. Sriram, and Z. Zhu,  ``Exploring neural transducers for end-to-end speech recognition,''  in Proc. ASRU, 2017.
  * [9] J. Li, R. Zhao, H. Hu, and Y. Gong,  ``Improving RNN transducer modeling for end-to-end speech recognition,''  in Proc. ASRU, 2019.
  * [10] A. Zeyer, A. Merboldt, R. Schlüter, and H. Ney,  ``A new training pipeline for an improved neural transducer,''  in Proc. Interspeech, 2020.
  * [11] Y. He, T. Sainath, R. Prabhavalkar, I. McGraw, R. Alvarez, D. Zhao, D. Rybach, A. Kannan, Y. Wu, R. Pang, Q. Liang, D. Bhatia, Y. Shangguan, B. Li, G. Pundak, K.C. Sim, T. Bagby, S. Chang, R. Rao, and A. Gruenstein,  ``Streaming end-to-end speech recognition for mobile devices,''  in Proc. ICASSP, Brighton, 2019.
  * [12] S. Chang, R. Prabhavalkar, Y. He, T.N. Sainath, and G. Simko,  ``Joint endpointing and decoding with end-to-end models,''  in Proc. ICASSP, Brighton, 2019.
  * [13] B. Li, S. Chang, T.N. Sainath, R. Pang, Y. He, T. Strohman, and Y. Wu,  ``Towards fast and accurate streaming end-to-end ASR,''  in Proc. ICASSP, Barcelona, 2020.
  * [14] Q. Zhang, H. Lu, H. Sak, A. Tripathi, E. McDermott, S. Koo, and S. Kumar,  ``Transformer Transducer: A streamable speech recognition model with Transformer encoders and RNN-T loss,''  in Proc. ICASSP, Barcelona, 2020.
  * [15] D. Le, M. Jain, G. Keren, S. Kim, Y. Shi, J. Mahadeokar, J. Chan, Y. Shangguan, C. Fuegen, O. Kalinli, Y. Saraf, and M.L. Seltzer,  ``Contextualized streaming end-to-end speech recognition with Trie-based deep biasing and shallow fusion,''  in Proc. Interspeech, Brno, 2021.
  * [16] Y. Shi, Y. Wang, C. Wu, C.-F. Yeh, J. Chan, F. Zhang, D. Le, and M. Seltzer,  ``Emformer: Efficient memory transformer based acoustic model for low latency streaming speech recognition,''  in Proc. ICASSP, Toronto, 2021.
  * [17] X. Chen, Y. Wu, Z. Wang, S. Liu, and J. Li,  ``Developing real-time streaming Transformer transducer for speech recognition on large-scale dataset,''  in Proc. ICASSP, Toronto, 2012.
  * [18] A. Graves, S. Fernandez, F. Gomez, and J. Schmidhuber,  ``Connectionist temporal classification: Labeling unsegmented sequenece data with recurrent neural networks,''  in Proc. ICML, Pittsburgh, 2006.
  * [19] S. Hochreiter and J. Schmidhuber,  ``Long short-term memory,''  Neural Computation, vol. 9, pp. 1735–1780, 1997.
  * [20] K. Rao, H. Sak, and R. Prabhavalkar,  ``Exploring architectures, data and units for streaming end-to-end speech recognition with RNN-transducer,''  in Proc. ASRU, Okinawa, 2017.
  * [21] A. Gulati, J. Qin, C.-C. Chiu, N. Parmar, Y. Zhang, J. Yu, W. Han, S. Wang, Z. Zhang, Y. Wu, and R. Pang,  ``Conformer: Convolution-augmented transformer for speech recognition,''  in Proc. Interspeech, Shanghai, 2020.
  * [22] A. Narayanan, T.N. Sainath, R. Pang, J. Yu, C.-C. Chiu, R. Prabhavalkar, E. Variani, and T. Strohman,  ``Cascaded encoders for unifying streaming and non-streaming ASR,''  in Proc. ICASSP, Toronto, 2021.
  * [23] Y. Zhang, S. Sun, and L. Ma,  ``Tiny transducer: A highly-efficient speech recognition model on edge devices,''  in Proc. ICASSP, 2021.
  * [24] G. Sun, C. Zhang, and P.C. Woodland,  ``Tree-constrained pointer generator for end-to-end contextual speech recognition,''  in Proc. ASRU, Cartagena, 2021.
  * [25] M. Ghodsi, X. Liu, J. Apfel, R. Cabrera, and E. Weinstein,  ``RNN transducer with stateless prediction network,''  in Proc. ICASSP, 2020.
  * [26] G. Saon, Z. Tüske, and K. Audhkhasi,  ``Alignment-length synchronous decoding for RNN transducer,''  in Proc. ICASSP, Barcelona, 2020.
  * [27] E. McDermott, H. Sak, and E. Variani,  ``A density ratio approach to language model fusion in end-to-end automatic speech recognition,''  in Proc. ASRU, Signapore, 2019.
  * [28] E. Variani, D. Rybach, C. Allauzen, and M. Riley,  ``Hybrid autoregressive transducer (HAT),''  in Proc. ICASSP, Barcelona, 2020.
  * [29] Z. Meng, N. Kanda, Y. Gaur, S. Parthasarathy, E. Sun, L. Lu, X. Chen, J. Li, and Y. Gong,  ``Internal language model training for domain-adaptive end-to-end speech recognition,''  in Proc. ICASSP, 2021.
  * [30] J. Li, R. Zhao, Z. Meng, Y. Liu, W. Wei, S. Parthasarathy, V. Mazalov, Z. Wang, L. He, S. Zhao, and Y. Gong,  ``Developing RNN-T models surpassing high-performance hybrid models with customization capability,''  in Proc. Interspeech, Shanghai, 2020.
  * [31] X. Zheng, Y. Liu, D. Gunceler, and D. Willett,  ``Using synthetic audio to improve the recognition of out-of-vocabulary words in end-to-end ASR systems,''  in Proc. ICASSP, Toronto, 2021.
  * [32] T. Doutre, W. Han, M. Ma, Z. Lu, C.-C. Chiu, R. Pang, A. Narayanan, A. Misra, Y. Zhang, and L. Cao,  ``Improving streaming automatic speech recognition with non-streaming model distillation on unsupervised data,''  in Proc. ICASSP, Toronto, 2021.
  * [33] G. Saon, Z. Tüske, D. Bolanos, and B. Kingsbury,  ``Advancing RNN transducer technology for speech recognition,''  in Proc. ICASSP, 2021.
  * [34] B. Li, R. Pang, T.N. Sainath, A. Gulati, Y. Zhang, J. Qin, P. Haghani, W.R. Huang, and M. Ma,  ``Scaling end-to-end models for large-scale multilingual ASR,''  in Proc. ASRU, Cartagena, 2021.
  * [35] F. Jelinek,  Statistical Methods for Speech Recognition,  MIT Press, 1997.
  * [36] C. Zhang, Z. Yang, X. He, and L. Deng,  ``Multimodal intelligence: Representation learning, information fusion, and applications,''  IEEE Journal of Selected Topics in Signal Proc., vol. 14, pp. 478–493, 2020.
  * [37] Y. Zhang, G. Chen, D. Yu, K. Yao, S. Khudanpur, and J. Glass,  ``Highway long short-term memory RNNs for distant speech recognition,''  in Proc. ICASSP, Shanghai, 2016.
  * [38] J. Tenenbaum and W. Freeman,  ``Separating style and content with bilinear models,''  Neural computation, vol. 12, pp. 1247–1283, 2000.
  * [39] H. Pirsiavash, D. Ramanan, and C. Fowlkes,  ``Bilinear classifiers for visual recognition,''  in Proc. NIPS, Vancouver, 2009.
  * [40] J.H. Kim, K.W. On, W. Lim, J. Kim, J.W. Ha, and B.T. Zhang,  ``Hadamard product for low-rank bilinear pooling,''  in Proc. ICLR, Toulon, 2017.
  * [41] G. Sun, C. Zhang, and P.C. Woodland,  ``Combination of deep speaker embeddings for diarisation,''  Neural Networks, vol. 141, pp. 372–384, 2021.
  * [42] S. Ramakrishnan, A. Agrawal, , and S. Lee,  ``Overcoming language priors in visual question answering with adversarial regularization,''  in Proc. NeurIPS, Montreal, 2018.
  * [43] D.S. Park, Y. Zhang, C.-C. Chiu, Y. Chen, B. Li, W. Chan, Q.V. Le, and Y. Wu,  ``SpecAugment on large scale datasets,''  in Proc. ICASSP, Toronto, 2020.
  * [44] M. Schuster and K. Nakajima,  ``Japanese and Korean voice search,''  in Proc. ICASSP, Kyoto, 2012.
  * [45] J. Shen, P. Nguyen, et al.,  ``Lingvo: A modular and scalable framework for sequence-to-sequence modeling,''  arXiv:2005.08100, 2019.