# Conditional moment methods for polydisperse cavitating flows

Spencer H. Bryngelson<EMAIL_ADDRESS>Rodney O. Fox  Tim Colonius  School of
Computational Science & Engineering, Georgia Institute of Technology, Atlanta,
GA 30332, USA  Department of Chemical and Biological Engineering, Iowa State
University, Ames, IA 50011, USA  Center for Multiphase Flow Research and
Education, Iowa State University, Ames, IA 50011, USA  Department of
Mechanical and Civil Engineering, California Institute of Technology,
Pasadena, CA 91125, USA

###### Abstract

The dynamics of cavitation bubbles are important in many flows, but their
small sizes and high number densities often preclude direct numerical
simulation. We present a computational model that averages their effect on the
flow over larger spatiotemporal scales. The model is based on solving a
generalized population balance equation (PBE) for nonlinear bubble dynamics
and explicitly represents the evolving probability density of bubble radii and
radial velocities. Conditional quadrature-based moment methods (QBMMs) are
adapted to solve this PBE. A one-way-coupled bubble dynamics problem
demonstrates the efficacy of different QBMMs for the evolving bubble
statistics. Results show that enforcing hyperbolicity during moment inversion
(CHyQMOM) provides comparable model-form accuracy to the traditional
conditional method of moments and decreases computational costs by about ten
times for a broad range of test cases. The CHyQMOM-based computational model
is implemented in MFC, an open-source multi-phase and high-order-accurate flow
solver. We assess the effect of the model and its parameters on a two-way
coupled bubble screen flow problem.

###### keywords:

Bubbly flow, cavitation, population balance modeling, quadrature-based moment
methods

††journal: J. Comp. Phys.

##  1 Introduction

Cavitating flows arise near ship propellers [1], in hydraulic machinery [2],
over spillways [3], and during medical therapies like lithotripsy and
histotripsy [4, 5, 6]. In these flows, polydisperse bubble clouds are
generated directly through nucleation or break-up and shedding larger vapor
regions. Directly simulating the flow of these dispersions is challenging
since the associated bubble dynamics frequently occur at smaller length and
time scales than the background flow [7].

In the dilute limit, a strategy for overcoming this scale separation
recognizes that the dynamics of each bubble are unimportant compared to
statistics of the bubble population. Ensemble phase-averaging [8, 9], for
example, results in equations for the continuous liquid face forced, in a two-
way-coupled manner, by statistics of (typically) the bubble radius radial
velocity. The associated bubble variables become stochastic, Eulerian fields.
This contrasts against Lagrangian models that track, model, and average over
samples of individual particles [10, 11]. Unless the bubble populations are
such that any small volume (compared to the mixture-averaged flow field
variations) contains a sufficient number of bubbles to describe the statistics
faithfully, these models simulate at best a single realization of a random
process.

The population balance equation (PBE) represents the evolution of a
dispersion’s probability (or number density function) according to a set of
internal variables [12, 13]. The internal variables must provide sufficient
information to close each particle’s dynamic and thermodynamic evolution. The
PBE is commonly used to model the particle size distribution. It depends
solely on particle motion, coagulation, and breakup [14, 15, 16], wherein the
internal variables are the spatial coordinates and velocities. PBE-based
models have successfully modeled flowing soot during combustion processes
[17], aerosol sprays [18], and more. For oscillating bubbles, the internal
variables must include the bubble radius, bubble radial velocity, and
sufficient information regarding the bubble contents, which can often
characterize with a single equilibrium radius [7]. When appropriate, these
three internal variables can be appended to those associated with relative
motion, coalescence, and break up [19, 20].

The current work formulates a PBE-based model for the distribution of bubble
radius, radial velocity, and equilibrium radius. Once formulated, the PBE is a
PDE in independent coordinates associated with each internal variable (as well
as spatiotemporal coordinates associated with the continuous flow in which the
bubbles reside). Solving this six-dimensional PDE via the method of lines is
intractable. Instead, the method of moments represents and evolves the number
density function via its statistical moments [21]. Other approaches to solving
PBEs exist, like class [22, 23], particle [24], and Monte Carlo [25, 26]
methods, though these are comparatively inefficient for multivariate PBEs or
large simulations [27].

The ensemble phase-averaging methods determine the governing flow equations
for oscillating bubble populations. As described in section 2, we develop an
Euler–Euler approach to represent the flow and the averaged bubble dynamics at
the sub-grid level [28, 10]. This method is two-way-coupled between the
dispersed bubbles and suspending liquid. In particular, moments of the
oscillating bubble dynamics alter the effective pressure and void fraction
evolution equation.

Next, in section 3, we develop and verify the conditional moment methods
discussed above to close these equations. This requires a
dynamical/thermodynamical model for each bubble. Even under the assumption of
spherical bubbles, the most general form of such a model consists of a set of
PDEs for the balance of mass, momentum, energy. Under further assumptions,
however, these can be reduced to a set of ODEs for the bubble radius and
radial velocity, with equilibrium radius as a parameter [29]. The set of ODEs
for each particle determines the derivatives of the internal variables that
close the moment transport equations [30]. Quadrature-based moment methods
then invert the moment set for a quadrature rule approximating the quantities
required to close the moment transport equations [31, 32]. With multiple
independent variables, conditional quadrature moment methods are
computationally preferable [33]. This work implements the conditional
quadrature method of moments (CQMOM [33]) and its hyperbolically constrained
version, CHyQMOM [34, 35, 36].

Section 4 describes our interface-capturing numerical algorithm for solving
the resulting system of equations. We demonstrate the methodology in section 5
by considering an acoustically excited bubble screen. Lastly, section 6
summarizes our results and potential next steps for extending and analyzing
the PBE-based method.

##  2 Model formulation

###  2.1 Compressible flow equations

A dilute suspension of dynamically evolving bubbles flow in a compressible
liquid is considered. For simplicity, we assume no slip between the bubbles
and the surrounding liquid and that the gas density is much smaller than the
liquid density. Under these assumptions, the mixture-averaged form of the
compressible flow equations are

$\displaystyle\frac{\partial\rho}{\partial t}+\nabla\cdot({\rho\bm{u}})=$
$\displaystyle\,0,$ $\displaystyle\frac{\partial{\rho\bm{u}}}{\partial
t}+\nabla\cdot(\rho\bm{u}\bm{u}+p\bm{I})=$ $\displaystyle\,0,$ (1)
$\displaystyle\frac{\partial{E}}{\partial t}+\nabla\cdot({E}+p)\bm{u}=$
$\displaystyle\,0,$

where $\rho$, $\bm{u}$, $p$, and $E$ are the density, velocity vector,
pressure, and total energy. Terms associated with the bubbles modify these
quantities and transport them in space according to ensemble phase-averaging,
discussed next.

###  2.2 Ensemble phase-averaging

The ensemble-averaged equations follow from Zhang and Prosperetti [8] and
Bryngelson et al. [10]. The disperse phase has a void fraction $\alpha$ and is
assumed to be a dilute ($\alpha\ll 1$) population of spherical bubbles. The
bubbles are represented statistically via random variables $R$, $\dot{R}$, and
$R_{o}$ corresponding to the instantaneous bubble radius, time derivative, and
equilibrium bubble radius. Section 2.3 presents this bubble dynamics model in
detail. The mixture-averaged pressure field is computed as

$\displaystyle p(\bm{x},t)=(1-\alpha)p_{\ell}+\alpha\left(\frac{\mkern
1.5mu\overline{\mkern-1.5muR^{3}p_{bw}\mkern-1.5mu}\mkern 1.5mu}{\mkern
1.5mu\overline{\mkern-1.5muR^{3}\mkern-1.5mu}\mkern 1.5mu}-\rho\frac{\mkern
1.5mu\overline{\mkern-1.5muR^{3}\dot{R}^{2}\mkern-1.5mu}\mkern 1.5mu}{\mkern
1.5mu\overline{\mkern-1.5muR^{3}\mkern-1.5mu}\mkern 1.5mu}\right),$ (2)

where $p_{bw}$ are the associated bubble wall pressure and
$p_{\ell}(\bm{x},t)$ is the liquid pressure according to the stiffened-gas
equation of state [37]:

$\displaystyle\Gamma_{\ell}p_{\ell}+\Pi_{\infty,\ell}=\frac{1}{1-\alpha}\left(E-\frac{1}{2}\rho\bm{u}^{2}\right).$
(3)

The coefficients of (3) represent water, with specific heat ratio
$\gamma_{\ell}=1+1/\Gamma_{\ell}=7.15$ and stiffness
$\Pi_{\infty,\ell}=$356\text{\,}\mathrm{MPa}$$ [38].

The bubble number density per unit volume $n(\bm{x},t)$ is conserved in the
absence of coalescence or breakup:

$\displaystyle\frac{\partial n}{\partial t}+\nabla\cdot(n\bm{u})=0.$ (4)

For the spherical bubbles considered here, $n$ is related to the void fraction
$\alpha$ via

$\displaystyle\alpha(\bm{x},t)=\frac{4}{3}\pi{\mkern
1.5mu\overline{\mkern-1.5muR^{3}\mkern-1.5mu}\mkern 1.5mu}n(\bm{x},t),$ (5)

and thus the void fraction $\alpha(\bm{x},t)$ transports as

$\displaystyle\frac{\partial\alpha}{\partial
t}+\bm{u}\cdot\nabla\alpha=3\alpha\frac{\mkern
1.5mu\overline{\mkern-1.5muR^{2}\dot{R}\mkern-1.5mu}\mkern 1.5mu}{\mkern
1.5mu\overline{\mkern-1.5muR^{3}\mkern-1.5mu}\mkern 1.5mu},$ (6)

where the right-hand-side represents the change in void fraction due to bubble
growth and collapse.

The over-barred terms appearing in (2), (5), and (6),

$\displaystyle\mkern
1.5mu\overline{\mkern-1.5muR^{3}\dot{R}^{2}\mkern-1.5mu}\mkern 1.5mu,\;\mkern
1.5mu\overline{\mkern-1.5muR^{3}\mkern-1.5mu}\mkern 1.5mu,\;\mkern
1.5mu\overline{\mkern-1.5muR^{2}\dot{R}\mkern-1.5mu}\mkern 1.5mu,\;\text{and
}\mkern 1.5mu\overline{\mkern-1.5muR^{3}p_{bw}\mkern-1.5mu}\mkern 1.5mu.$ (7)

denote average quantities of the bubble dispersion. In particular, they are
raw moments $\mu_{lmn}$ with respect to a bubble number density function
$f(R,\dot{R},R_{o})$,

$\displaystyle\mu_{lmn}=\mkern
1.5mu\overline{\mkern-1.5muR^{l}\dot{R}^{m}R_{o}^{n}\mkern-1.5mu}\mkern
1.5mu=\int_{\Omega}R^{l}\dot{R}^{m}R_{o}^{n}f(R,\dot{R},R_{o})\,\text{d}R\text{d}\dot{R}\text{d}R_{o},$
(8)

which are computed via the methods of section 3.

###  2.3 Bubble dynamics model

Even for spherical models, dozens of available models employ differing
assumptions about the bubble contents and simplifications of the physics.
These models range from a system of two or more ODEs up to a set of PDEs (in a
radial coordinate and time) describing the balance of mass, momentum, and
energy of each bubble. PDE-based models are deemed intractable (at present),
and we focus on ODE-based modeling. We choose one of the simplest possible
models to demonstrate our methodology but note that our framework can be
extended to arbitrarily complex ODE-based models.

We assume the spherical bubbles are filled with noncondensible gas and that,
insofar as their dynamics are concerned, the liquid may be assumed
incompressible. We assume that the gas undergoes a polytropic process during
compression and expansion. Under these assumptions, the bubble radius is
governed by a Rayleigh–Plesslet-like equation

$\displaystyle
R\ddot{R}+\frac{3}{2}\dot{R}^{2}+\frac{4}{\mbox{Re}}\frac{\dot{R}}{R}=\left(\frac{R_{o}}{R}\right)^{3\gamma}-\frac{1}{C_{p}}-\frac{2}{\mbox{We}\,R_{o}}\left[\frac{R_{o}}{R}-\left(\frac{R_{o}}{R}\right)^{3\gamma}\right],$
(9)

which is dimensionless via the reference bubble size $R_{o}^{\ast}$ and
ambient liquid pressure $p_{0}$, and density $\rho_{0}$. The polytropic index
is $\gamma$; in examples below we use $\gamma=1.4$. In (9), $C_{p}\equiv
p_{0}/p_{\ell}$ is the forcing pressure ratio and Reynolds and Weber numbers
correspond to viscous and surface tension effects as

$\displaystyle\mbox{Re}\equiv\sqrt{\frac{p_{0}}{\rho_{0}}}\frac{R_{o}^{\ast}}{\nu_{0}}\quad\text{and}\quad\mbox{We}\equiv\frac{p_{0}R_{o}^{\ast}}{S},$
(10)

where $S$ is the air–water surface tension coefficient and $\nu_{0}$ is the
liquid kinematic viscosity. Under these conditions, the bubble wall pressure
of (2) simplifies to $p_{bw}=\left(R_{o}/R\right)^{3\gamma}$ and the last
moment of (7) reduces to $\mkern
1.5mu\overline{\mkern-1.5muR^{3}(R_{o}/R)^{3\gamma}\mkern-1.5mu}\mkern 1.5mu$.

##  3 Population balance formulation

Figure 1: Illustration of the quadrature-based moment method for a fully
coupled bubbly flow.

The population balance equation (PBE)

$\displaystyle\frac{\partial f}{\partial t}+\frac{\partial}{\partial
R}(f\dot{R})+\frac{\partial}{\partial\dot{R}}(f\ddot{R})=0$ (11)

governs the conditional PDF $f(R,\dot{R}|R_{o})$ in the absence of bubble
coalescence or breakup, though this approach can naturally accommodate these
effects if required. Figure 1 summarizes the PBE–quadrature approach.

###  3.1 Method of moments

Following the usual procedure, a finite set of raw moments
$\vec{\boldsymbol{\mu}}$ represent $f$ per (8) [14]. The specific moments that
make up $\vec{\boldsymbol{\mu}}$ depend on the inversion algorithm and are
defined in the appendices A and B. These moments transport on the grid and
evolve according to the bubble dynamics of section 2.3 as

$\displaystyle\frac{\partial n\vec{\boldsymbol{\mu}}}{\partial
t}+\nabla\cdot(n\vec{\boldsymbol{\mu}}\bm{u})=n{\dot{\vec{\boldsymbol{\mu}}}}=n\vec{\bm{g}}$
(12)

where

$\displaystyle
g_{lmn}=l\mu_{l-1,m+1,n}+m\iiint_{\Omega}\ddot{R}R^{l}\dot{R}^{m-1}R_{o}^{n}f(\vec{\boldsymbol{\mu}})\,\text{d}R\text{d}\dot{R}\text{d}R_{o}$
(13)

and
$\Omega=\Omega_{R}\times\Omega_{\dot{R}}\times\Omega_{R_{o}}=(0,\infty)\times(-\infty,\infty)\times(0,\infty)$.
The integrand of (13) is closed via the bubble dynamics model (9) for
$\ddot{R}$, and the integral is approximated via quadrature, which is
discussed next.

###  3.2 Conditional quadrature moment inversion

Since $R_{o}$ is not a dynamic variable, the number density function is split
as

$\displaystyle f(R,\dot{R},R_{o})=f(R,\dot{R}|R_{o})f(R_{o}).$ (14)

The raw moments (8) are then

$\displaystyle\mu_{lmn}$
$\displaystyle\equiv\int_{\Omega_{R_{o}}}f(R_{o})R_{o}^{m}\mu_{lm}(R_{o})\,\text{d}R_{o}$
(15)
$\displaystyle\approx\sum_{i=1}^{N_{R_{o}}}w_{i}\widehat{R}_{o,i}^{n}\,\mu_{lm}(\widehat{R}_{o,i}),$
(16)

with $N_{R_{o}}$ time-independent weights $w_{i}$ and nodes
$\widehat{R}_{o,i}$. Simpson’s rule computes these nodes and weights and the
accuracy in the approximation of (16) depends on $N_{R_{o}}$. Other numerical
integration methods are suitable and will be discussed in section 5.3. The
$R_{o}$-conditioned moments are

$\displaystyle\mu_{lm}(\widehat{R}_{o,i})$
$\displaystyle\equiv\iint_{\Omega_{R,\dot{R}}}f(R,\dot{R}|\widehat{R}_{o,i})R^{l}\dot{R}^{m}\,\text{d}R\text{d}\dot{R}$
(17)
$\displaystyle\approx\sum_{j=1}^{N_{R}}\sum_{k=1}^{N_{\dot{R}}}\left[\widehat{w}_{j,k}\widehat{R}_{j}^{l}\widehat{\dot{R}}_{k}^{m}\right]_{\widehat{R}_{o,i}}.$
(18)

The moment indices comprising the moment set of (18) are determined by the
conditional quadrature moment method used to invert those moments for a set of
quadrature points and weights (and the number of points desired) [33, 35]. In
particular, $\mu_{lm}(\widehat{R}_{o,i})$ is traded for quadrature points
$\\{\widehat{R}_{j},\widehat{\dot{R}}_{k}\\}(\widehat{R}_{o,i})$ and weights
$\widehat{w}_{j,k}(\widehat{R}_{o,i})$ for each $i=1,\dots,N_{R_{o}}$ (with
$j=1,\dots,N_{R}$; $k=1,\dots,N_{\dot{R}}$). The total moments (16) are then
approximated by substituting (18) into (16) as

$\displaystyle\mu_{lmn}=\sum_{i=1}^{N_{R_{o}}}w_{i}\widehat{R}_{o,i}^{n}\sum_{j=1}^{N_{R}}\sum_{k=1}^{N_{\dot{R}}}\left[\widehat{w}_{j,k}\widehat{R}_{j}^{l}\widehat{\dot{R}}_{k}^{m}\right]_{\widehat{R}_{o,i}}.$
(19)

The CHyQMOM and CQMOM algorithms are described in appendices A and B. Their
implementation is verified by ensuring that the error in closing a linear set
of moment transport equations is comparable to a finite-precision round-off
[39].

###  3.3 Model performance

We consider the ability of the QBMMs to represent the statistics of the bubble
dynamics. For this example, we take a monodisperse population with the same
equilibrium radius $R_{o}=1$ (and thus $N_{R_{o}}=1$). $R$ and $\dot{R}$ are
initialized with independent log-normal and normal distributions with
variances $\sigma_{R}^{2}$ and $\sigma_{\dot{R}}^{2}$, respectively. The
bubbles are forced by a step-change in pressure at $t=0^{+}$ from equilibrium
to a value $C_{p}$ (which is varied). This represents a stringent test of
model fidelity as high values of $C_{p}$ will induce strongly nonlinear bubble
dynamics.

Figure 2:  Example first- and second-order moments $\mu$ as labeled for
bubbles with Reynolds and Weber numbers corresponding to (a)
$R_{o}^{\ast}=$1\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ and (b)
$R_{o}^{\ast}=$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$$. The dimensionless
initial shape parameters are $\sigma_{R}=\sigma_{\dot{R}}=0.2$. The CHyQMOM
implementation uses two quadrature nodes in each internal coordinate direction
($2\times 2$), and “Gaussian” corresponds to Gaussian closure [10]. Monte
Carlo serves as an exact solution. The final time is $T=13.9$ and represents
10 natural periods of the largest bubble.

Figure 2 shows the evolution of the carried moment set for $C_{p}=0.3$. To
assess the accuracy, the CHyQMOM results are compared against Monte-Carlo
solutions that used $10^{4}$ samples to ensure that the sampling error was at
least 10-times smaller than the QBMM model-form error. Thus the Monte-Carlo
solutions can be regarded as exact solutions for these comparisons. The radial
moments $\mu_{10}$ and $\mu_{20}$ show how the bubbles oscillate in response
to the change in external pressure, eventually reaching a statistical
equilibrium [40]. These oscillations are damped rapidly for smaller bubbles
due to stronger viscous effects (smaller Re). Thus, larger bubbles entail
larger model-form errors as they build over time. This is most prominent for
the $\dot{R}$ moments, such as $\mu_{02}$, and smaller bubbles ((b)
$R_{o}^{\ast}=$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$$) as the velocity
statistics change more quickly than the radial ones.

Figure 3: Relative model-form errors for varying pressure ratio $C_{p}$ and
closure methods as labeled. All cases have $\mbox{Re}=100$ and
$\mbox{We}=13.9$.

A QBMM model-form relative $L_{2}$ error is

$\displaystyle\varepsilon\equiv\frac{1}{N_{t}}\sqrt{\sum_{i=1}^{N_{t}}\left[\frac{\mu^{\mathrm{(QBMM)}}(t_{i})-\mu^{\mathrm{(MC)}}(t_{i})}{\mu^{\mathrm{(MC)}}(t_{i})}\right]^{2}},$
(20)

where Monte-Carlo (MC) simulations serve as surrogate truth data and $t_{i}$
are $N_{t}=10^{3}$ uniformly spaced times in the time interval $t\in[0,T]$.
Figure 3 shows $\varepsilon$ for select first- and second-order moments $\mu$.
We see that $\varepsilon$ is smaller for CHyQMOM and CQMOM than Gaussian
closure, though the difference is modest and more strongly associated with
$C_{p}$. The larger errors for smaller $C_{p}$ are associated with stronger
bubble dynamics and the formation of non-Gaussian statistics, like skewness
and kurtosis, that these closures do not represent [41]. One can represent
higher-order statistics by carrying higher-order moments and thus inverting
for a higher-order quadrature rule. However, this is computationally
cumbersome and numerically unstable in the cases tested here.

Figure 4: Total simulation times for the cases of figure 3. A single-core of a
$2.8\text{\,}\mathrm{GHz}$ Intel i7-7700HQ performed the computations. Only
relative times are relevant, so the cheapest simulation normalizes the scale
(CHyQMOM, $C_{p}=0.8$).

Model accuracy was similar for all of the closure methods, but the
computational time to solution is not. Figure 4 shows the total relative
simulation time for the cases of figure 3 under the same adaptive time-step
tolerance. CHyQMOM simulations are about 10-times cheaper than those using
CQMOM or Gaussian closures, even though the accuracy is the same. This is
attributed to larger time step sizes (given the same error constraint) and a
smaller carried moment set than CQMOM. Due to its low cost for the same
accuracy, the fully-coupled simulation algorithm employs CHyQMOM and is
discussed next.

##  4 Interface-capturing flow solution algorithm

The QBMM approach of the previous section is solved using interface-capturing
numerics. MFC, an open-source flow solver, is the basis for the implementation
[42]. A brief description of the algorithm follows.

The governing equations (1), (6), and (12) combine as

$\displaystyle\frac{\partial\vec{\bm{q}}_{c}}{\partial
t}+\nabla\cdot\vec{\bm{F}}=\vec{\bm{r}}$ (21)

where
$\vec{\bm{q}}_{c}\equiv\\{\rho,\rho\bm{u},E,\alpha,n\vec{\boldsymbol{\mu}}\\}$
are the conservative variables, $\vec{\bm{F}}$ are the advective fluxes, and
$\vec{\bm{r}}$ are diffusive source terms. Finite volumes with uniform size
$\Delta$ discretize the domain.

###  4.1 Flow state initialization

We initialize the flow as follows. The number density function
$f(R,\dot{R},R_{o})$ is independently distributed with log-normal, normal, and
log-normal shapes in the $R$, $\dot{R}$, and $R_{o}$ directions with expected
values $\mathbb{E}[R]=\mathbb{E}[R_{o}]=1$ and $\mathbb{E}[\dot{R}]=0$ and
shape parameters $\sigma_{\cdot}$.

The NDF $f$ initializes the moment set $\vec{\boldsymbol{\mu}}$ via
integration. The primitive variables
$\vec{\bm{q}}_{p}\equiv\\{\rho,\bm{u},p,\alpha,\vec{\boldsymbol{\mu}}\\}$ are
initialized as appropriate. The bubble number density $n$ follows from (5),
$\vec{\bm{q}}_{p}$, and $\mu_{300}$. Lastly, the known primitive variables are
used to compute the conservative ones as
$\vec{\bm{q}}_{c}\equiv\\{\rho,\rho\bm{u},E,\alpha,n\vec{\boldsymbol{\mu}}\\}$.

1:$f(R,\dot{R},R_{o})\leftarrow$ Presumed form, uncorrelated,
$\\{\mu_{R},\mu_{\dot{R}},\mu_{R_{o}}\\}$ and
$\\{\sigma_{R},\sigma_{\dot{R}},\sigma_{R_{o}}\\}$

2:$\vec{\boldsymbol{\mu}}\leftarrow f$ and (8)

3:$\vec{\bm{q}}_{p}\equiv\\{\rho,\bm{u},p,\alpha,\vec{\boldsymbol{\mu}}\\}\leftarrow$
Patches

4:$n\leftarrow\alpha,\mu_{300}$ and (5)

5:$\vec{\bm{q}}_{c}\leftarrow n,\vec{\bm{q}}_{p}$

Algorithm 1 Flow initialization procedure

###  4.2 Flux divergence computation

A fifth-order-accurate WENO [43] scheme reconstructs the primitive variables
$\vec{\bm{q}}_{p}$ and the HLLC approximate Riemann solver [44] computes the
fluxes. Algorithm 2 describes this process in detail. High-order WENO
reconstructions do not guarantee that the reconstructed moments are
realizable, though the moment sets remained invertible in the subsequent
simulations.

1:$\vec{\bm{q}}_{p}\leftarrow\vec{\bm{q}}_{c}$ and (3)

2:$\widehat{\bm{q}}_{p}\leftarrow\texttt{WENO}(\vec{\bm{q}}_{p})$

3:$\\{\widehat{R},\widehat{\dot{R}},\widehat{w}\\}_{jk}\leftarrow\texttt{QBMM}(\widehat{\boldsymbol{\mu}})$

4:$\widehat{\bm{g}},\\{\mkern
1.5mu\overline{\mkern-1.5muR^{3}\dot{R}^{2}\mkern-1.5mu}\mkern 1.5mu,\mkern
1.5mu\overline{\mkern-1.5muR^{3}\mkern-1.5mu}\mkern 1.5mu,\mkern
1.5mu\overline{\mkern-1.5muR^{2}\dot{R}\mkern-1.5mu}\mkern 1.5mu,\mkern
1.5mu\overline{\mkern-1.5muR^{3}p_{bw}\mkern-1.5mu}\mkern
1.5mu\\}\leftarrow\\{\widehat{R}_{o},w\\}_{i},\\{\widehat{R},\widehat{\dot{R}},\widehat{w}\\}_{jk,i}$

5:$\vec{\bm{F}}\leftarrow\texttt{HLLC}(\widehat{\bm{q}}_{p},\widehat{\bm{g}},\mkern
1.5mu\overline{\mkern-1.5mu\\{\cdot\\}\mkern-1.5mu}\mkern 1.5mu)$

6:$\nabla\cdot\vec{\bm{F}}\leftarrow\vec{\bm{F}}$

Algorithm 2 Algorithm for computing the flux divergence.

###  4.3 Time stepping

The conservative variables are integrated in time using third-order-accurate
SSP–RK3 time integration [45]. Once the spatial derivatives have been
approximated, (21) becomes a semi-discrete system of ordinary differential
equations in time. We treat the temporal derivative using a Runge–Kutta time-
marching scheme for the state variables. To achieve high-order accuracy and
avoid spurious oscillations, we use the third-order-accurate total variation
diminishing scheme of Gottlieb and Shu [46]:

$\displaystyle\vec{\bm{q}}_{c}^{\,(1)}$
$\displaystyle=\vec{\bm{q}}_{c}^{\,\tilde{n}}+\Delta
t\bm{L}(\vec{\bm{q}}_{c}^{\,\tilde{n}}),$
$\displaystyle\vec{\bm{q}}_{c}^{\,(2)}$
$\displaystyle=\frac{3}{4}\vec{\bm{q}}_{c}^{\,\tilde{n}}+\frac{1}{4}\vec{\bm{q}}_{c}^{(1)}+\frac{1}{4}\Delta
t\bm{L}(\vec{\bm{q}}_{c}^{\,(1)}),$ (22)
$\displaystyle\vec{\bm{q}}_{c}^{\,\tilde{n}+1}$
$\displaystyle=\frac{1}{3}\vec{\bm{q}}_{c}^{\,\tilde{n}}+\frac{2}{3}\vec{\bm{q}}_{c}^{\,(2)}+\frac{2}{3}\Delta
t\bm{L}(\vec{\bm{q}}_{c}^{\,(2)}),$

where superscripts $(1)$ and $(2)$ indicate intermediate time-step stages,
$\bm{L}$ represents the portion of (21) that does not include the time
derivative, and $\tilde{n}$ is the time-step index.

##  5 Application to polydisperse bubble screens

###  5.1 Problem setup

Figure 5:  Schematic of the acoustically excited bubble screen test problem.
The screen is centered at $x=0$.

We assess the statistics of bubble dynamics in an idealized model problem
consisting of an acoustically excited dilute bubble screen. The bubble screen
parameterization matches that of Bryngelson et al. [10], with a length of
$L_{s}=$5\text{\,}\mathrm{mm}$$, initial void fraction $\alpha_{o}=10^{-4}$,
median bubble equilibrium size
$R_{o}^{\ast}=$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ and log-normal
variance $\sigma_{R_{o}}$. The domain is $L_{d}=5L_{s}$ long, and its
boundaries are non-reflective via characteristic-based boundary conditions
[47]. The one-way (positive $\hat{x}$-direction) sound wave $p_{\infty}(t)$ is
generated via source terms in the governing equations (21) according to
Bryngelson et al. [10]. Its form is a single period of a sinusoid with peak
amplitude $0.3p_{0}$ and frequency $300\text{\,}\mathrm{kHz}$.

###  5.2 Bubble screen behavior

Figure 6: Bubble-screen-centered pressure for varying $R_{o}$ log-normal
distributions with shape parameter $\sigma_{R_{o}}$ and fixed
$\sigma_{R}=\sigma_{\dot{R}}=0.2$.

We start by considering a screen with fixed dynamic coordinate distributions
$\sigma_{R}=\sigma_{\dot{R}}$, but varying distributions of equilibrium sizes
$\sigma_{R_{o}}$. Polydispersity in $R_{o}$ is integrated via Simpson’s rule
$61$ quadrature points for all cases. Figure 6 shows the bubble screen
pressure for these cases as they evolve in time. For larger $\sigma_{R_{o}}$
(or broader distributions or bubble equilibrium sizes), the pressure is less
oscillatory in time. A similar observation was made by Bryngelson et al. [10]
for cases with no $R$ or $\dot{R}$ distributions. In the case of figure 6, we
instead observe high-frequency oscillations in addition to the long-wavelength
behaviors associated with the impinging pressure wave $p_{\infty}$. The origin
of these oscillations is discussed next.

Figure 7: Bubble-screen-centered pressure before, during, and after excitement
due to an acoustic wave. The bubbles are polydisperse with log-normal $R_{o}$
distribution ($\sigma_{R_{o}}=0.2$) and $\mbox{Re}=10^{3}$. Variations in (a)
$\sigma_{R}$ and (b) $\sigma_{\dot{R}}$ are shown about a
$\sigma_{R}=\sigma_{\dot{R}}=0.2$ representative state.

Figure 7 shows the dynamics associated with a bubble screen in varying degrees
of statistical disequilibrium, represented via different $\sigma_{R}$ and
$\sigma_{\dot{R}}$. Figure 7 (a) fixes $\sigma_{\dot{R}}$ and varies
$\sigma_{R}$. We observe the shorter-wavelength oscillatory behavior, observed
in figure 6, becoming more prominent for larger $\sigma_{R}$. These
wavelengths are commensurate with the mean bubble natural frequencies, which
superimpose the longer wavelength acoustics associated with the impinging
$p_{\infty}$ wave. Figure 7 (b) shows a smoother pressure profile for larger
$\sigma_{\dot{R}}$. Phase-cancellation between the larger waves associated
with broader $\sigma_{\dot{R}}$ distributions and those of the $\sigma_{R}$
distributions may account for this behavior. Notably, these behaviors are
qualitatively similar to those associated with varying $R_{o}$ distribution
widths. Thus, parameterizing an $R_{o}$ distribution based on single-probe
pressure measurements is insufficient.

###  5.3 Closure errors

We quantify the moment closure error, $\varepsilon_{c}$, via the mismatch in
bubble screen pressure $p(t,x=0)$ due to truncated $R_{o}$ integration.
Otherwise, the definition of $\varepsilon_{c}$ has the same form as
$\varepsilon$ of (20), though the truth values are changed from Monte-Carlo
simulations (which are unavailable) to high-resolution $N_{R_{o}}=401$
simulations as

$\displaystyle\varepsilon_{c}\equiv\frac{1}{N_{t}}\sqrt{\sum_{i=1}^{N_{t}}\left[\frac{p^{\mathrm{(QBMM)}}(t_{i},x=0)-p^{(N_{R_{o}}=401)}(t_{i},x=0)}{p^{(N_{R_{o}}=401)}(t_{i},x=0)}\right]^{2}}.$
(23)

This choice is made for two reasons. First, extending the QBMM method to
additional $R$ and $\dot{R}$ quadrature points (or moments, equivalently) was
found numerically unstable for most bubble cavitation problems. One approach
to addressing this specific problem is introducing a recurrent neural network
to correct the quadrature points and weights [48], though we do not discuss it
further here. Second, we will see that the closure errors are most strongly
associated with $N_{R_{o}}$, and thus focus on the influence of this
parameter.

Figure 8:  Relative closure error $\varepsilon_{c}$ (defined in (23)) for
increasing number of $R_{o}$-direction quadrature points $N_{R_{o}}$.
Variations in (a) $\sigma_{R}$, (b) $\sigma_{\dot{R}}$, and (c)
$\sigma_{R_{o}}$ are shown. Unless labeled otherwise, cases have the baseline
$\sigma_{R}=\sigma_{\dot{R}}=\sigma_{R_{o}}=0.2$.

Figure 8 shows the $N_{R_{o}}$ closure errors associated with variations in
all three PDF directions (panels a–c). For varying $\sigma_{R}$ (a) and
$\sigma_{R_{o}}$ (c) we see that increasing variance $\sigma_{\cdot}$ results
in larger closure errors $\varepsilon_{c}$. This appears to be associated with
the larger pressure oscillations $p(x=0,t)$ observed for such cases. For
figure 8 (b), we see the reverse trend, with larger $\sigma_{\dot{R}}$
corresponding to smaller closure errors, though this effect is small. Indeed,
this effect matches that of the $R$\- and $R_{o}$-direction effects, where
larger $\sigma_{\dot{R}}$ results in smoother pressure histories. We also
investigated the effect of Re, which results in smoother pressure profiles for
smaller Re, and found the same trend.

Figure 9: Closure error associated with
$\sigma_{R}=\sigma_{\dot{R}}=\sigma_{R_{o}}=0.2$ for different
$R_{o}$-direction quadrature rules as labeled.

Compelled by the closure error associated with $N_{R_{o}}$ dominating the
total simulation error, we also implemented Gaussian–Hermite and
Gauss–Legendre quadrature methods in the $R_{o}$-direction in an attempt to
reduce the error. Figure 9 shows these closures errors $\varepsilon_{c}$ for
increasing $N_{R_{o}}$. We see that these alternative quadrature rules do not
significantly change the closure error, despite being optimal for a given
$N_{R_{o}}$. This is because the integrand becomes highly oscillatory without
viscous damping, as bubbles of different equilibrium radii $R_{o}$ become out
of phase [40]. Long-time simulations without significant viscous effects,
broadly polydisperse bubble populations (large $\sigma_{R_{o}}$) will require
treatment of these developing oscillations. For example, a time-stationary
analysis could be used to approximate the $R_{o}$-moments, so long as the
bubble dynamics are fast compared to the fluid flow [40]. We will investigate
this in future work.

##  6 Conclusion

A fully coupled numerical method for simulating sub-grid cavitating bubble
dispersions was presented. The approach represents and evolves the statistics
of the bubble dynamic variables. This work built upon the well-established
framework of quadrature-based moment methods: a governing set of moment
transport equations were derived, effectively representing the underlying
statistics. Those moments were inverted for a quadrature rule to close the
fully coupled disperse flow equations.

Our results showed that only four quadrature points, two in each of the bubble
dynamics coordinates, can be sufficient to represent the statistics of the
monodisperse bubbles. This model was particularly accurate for weak pressure
forcing, and thus nearly linear bubble dynamics. It was also accurate for
highly damped dynamics, characteristic of low Reynolds numbers, as non-
Gaussian statistics cannot present themselves. This result contrasts against a
previous approach that assumed Gaussian statistics, which had demonstrably
worse performance and higher computational costs [41]. The method would
require more quadrature points for cases with more complicated dynamics, such
as larger forcing.

An acoustically excited dilute bubble screen problem demonstrated the model’s
behavior in a coupled flow. Modeling the statistics in the bubble dynamics
variables resulted in qualitatively different behavior in the screen region.
For example, increasing the distribution breadth in the bubble radius
coordinate resulted in short-wavelength pressure oscillations of increasing
magnitude superimposing the background response to the impinging acoustic
wave. Thus, modeling the $R$–$\dot{R}$ distributions is potentially critical
to modeling cavitating bubble clouds and certainly important if bubbles are
ever in such statistical disequilibrium.

Finally, we found that for broadly polydisperse populations (large
$\sigma_{R_{o}}$), the closure error in the equilibrium radius coordinate
$R_{o}$ was most important. Our results also show that phase-cancellation can
modestly reduce some computational costs associated with resolving the $R_{o}$
coordinate. Still, $R_{o}$-direction quadrature dominates the solution cost of
these polydisperse bubble populations, even in the face of more sophisticated
quadrature rules like Gauss–Hermite (corresponding to the underlying log-
normal distribution). However, one can apply existing approaches if time-
separation is achieved, or a novel method based on, e.g., neural networks, can
also treat this. We will investigate these in-depth in future work.

## Acknowledgements

We thank Alexis Charalampopoulos and Themis Sapsis for invaluable discussions
regarding the presented method. The US Office of Naval Research supported this
work under grant numbers N0014-17-1-2676 and N0014-18-1-2625. Computations
were performed via the Extreme Science and Engineering Discovery Environment
(XSEDE) under allocations TG-CTS120005 (PI Colonius) and TG-PHY210084 (PI
Bryngelson), supported by National Science Foundation grant number
ACI-1548562.

## References

  * Cook [1928] S. S. Cook,  Erosion by water-hammer,  Proc. Royal Soc. London A 119 (1928) 481–488.
  * Arndt [1981] R. Arndt,  Cavitation in fluid machinery and hydraulic structures,  Annu. Rev. Fluid Mech. 13 (1981) 273–326.
  * Falvey [1990] H. T. Falvey, Cavitation in chutes and spillways, US Department of the Interior, Bureau of Reclamation Denver, 1990.
  * Cleveland et al. [2000] R. O. Cleveland, O. A. Sapozhnikov, M. R. Bailey, L. A. Crum,  A dual passive cavitation detector for localized detection of lithotripsy-induced cavitation in vitro,  J. Acoust. Soc. Am. 107 (2000) 1745–1758.
  * Pishchalnikov et al. [2003] Y. A. Pishchalnikov, O. A. Sapozhnikov, M. R. Bailey, J. C. Williams, R. O. Cleveland, T. Colonius, L. A. Crum, A. P. Evan, J. A. McAteer,  Cavitation bubble cluster activity in the breakage of kidney stones by lithotripter shockwaves,  J. Endourol. 17 (2003) 435–446.
  * Maxwell et al. [2011] A. D. Maxwell, T.-Y. Wang, C. A. Cain, J. B. Fowlkes, O. A. Sapozhnikov, M. R. Bailey, Z. Xu,  Cavitation clouds created by shock scattering from bubbles during histotripsy,  J. Acoust. Soc. Am. 130 (2011) 1888–1898.
  * Brennen [1995] C. E. Brennen, Cavitation and bubble dynamics, Oxford University Press, 1995\.
  * Zhang and Prosperetti [1994a] D. Z. Zhang, A. Prosperetti,  Ensemble phase-averaged equations for bubbly flows,  Phys. Fluids 6 (1994a).
  * Zhang and Prosperetti [1994b] D. Z. Zhang, A. Prosperetti,  Averaged equations for inviscid disperse two-phase flow,  J. Fluid Mech. 267 (1994b) 185–219.
  * Bryngelson et al. [2019] S. H. Bryngelson, K. Schmidmayer, T. Colonius,  A quantitative comparison of phase-averaged models for bubbly, cavitating flows,  Int. J. Mult. Flow 115 (2019) 137–143.
  * Capecelatro and Desjardins [2013] J. Capecelatro, O. Desjardins,  An Euler–Lagrange strategy for simulating particle-laden flows,  J. Comp. Phys. 238 (2013) 1–31.
  * Randolph and Larson [1987] A. D. Randolph, M. A. Larson, Theory of particulate processes, San Diego, Academic Press, 1987.
  * Ramkrishna [2000] D. Ramkrishna, Population Balances, Academic Press, New York, USA, 2000.
  * Fox [2003] R. O. Fox, Computational models for turbulent reacting flows, Cambridge University Press, 2003\.
  * Rigopoulos [2010] S. Rigopoulos,  Population balance modelling of polydispersed particles in reactive flows,  Prog. Energy Combust. Sci. 36 (2010) 412–443.
  * Lee et al. [2015] K. F. Lee, R. I. Patterson, W. Wagner, M. Kraft,  Stochastic weighted particle methods for population balance equations with coagulation, fragmentation and spatial inhomogeneity,  J. Comp. Phys. 303 (2015) 1–18.
  * Mueller et al. [2009] M. E. Mueller, G. Blanquart, H. Pitsch,  A joint volume-surface model of soot aggregation with the method of moments,  Proc. Combust. Inst. 32 I (2009) 785–792.
  * Sibra et al. [2017] A. Sibra, J. Dupays, A. Murrone, F. Laurent, M. Massot,  Simulation of reactive polydisperse sprays strongly coupled to unsteady flows in solid rocket motors: Efficient strategy using Eulerian multi-fluid methods,  J. Comp. Phys. 339 (2017) 210–246.
  * Carrica et al. [1999] P. Carrica, D. Drew, F. Bonetto, R. Lahey,  A polydisperse model for bubbly two-phase flow around a surface ship,  Int. J. Mult. Flow 25 (1999) 257–305.
  * Heylmun et al. [2019] J. C. Heylmun, B. Kong, A. Passalacqua, R. Fox,  A quadrature-based moment method for polydisperse bubbly flows,  Comp. Phys. Comm. (2019).
  * Hulburt and Katz [1964] H. M. Hulburt, S. Katz,  Some problems in particle technology. A statistical mechanical formulation,  Chem. Eng. Sci. 19 (1964) 555–574.
  * Hounslow et al. [1988] M. J. Hounslow, R. L. Ryall, V. R. Marshall,  A discretized population balance for nucleation, growth, and aggregation,  AIChE J. 34 (1988) 1821–1832.
  * Vanni [2000] M. Vanni,  Approximate population balance equations for aggregation breakage processes,  J. Colloid Interface Sci. 221 (2000) 143–160.
  * Patterson et al. [2011] R. I. Patterson, W. Wagner, M. Kraft,  Stochastic weighted particle methods for population balance equations,  J. Comp. Phys. 230 (2011) 7456–7472.
  * Zhao et al. [2007] H. Zhao, A. Maisels, T. Matsoukas, C. Zheng,  Analysis of four monte-carlo methods for the solution of population balances in dispersed systems,  Powder Technol. 173 (2007) 38–50.
  * Rosner et al. [2003] D. E. Rosner, R. McGraw, P. Tandon,  Multivariate population balances via moment and Monte Carlo simulation methods: an important sol reaction engineering bivariate example and “mixed” moments for the estimation of deposition, scavenging, and optical properties for populations of nonspherical suspended particles,  Indust. Eng. Chem. Res. 42 (2003) 2699–2711.
  * Zucca et al. [2007] A. Zucca, D. L. Marchisio, M. Vanni, A. A. Barresi,  Validation of bivariate DQMOM for nanoparticle processes simulation,  AiChE J. 53 (2007) 918–931.
  * Ando et al. [2011] K. Ando, T. Colonius, C. E. Brennen,  Numerical simulation of shock propagation in a polydisperse bubbly liquid,  Int. J. Mult. Flow 37 (2011) 596–608.
  * Plesset and Prosperetti [1977] M. S. Plesset, A. Prosperetti,  Bubble dynamics and cavitation,  Annu. Rev. Fluid Mech. 9 (1977) 145–185.
  * Frenklach and Harris [1987] M. Frenklach, S. Harris,  Aerosol dynamics modeling using the method of moments,  J. Colloid Interface Sci. 118 (1987) 252–261.
  * McGraw [1997] R. McGraw,  Description of aerosol dynamics by the quadrature method of moments,  Aerosol Sci. Technol. 27 (1997) 255–265.
  * Marchisio and Fox [2013] D. L. Marchisio, R. O. Fox, Computational models for polydisperse particulate and multiphase systems, Cambridge University Press, 2013\.
  * Yuan and Fox [2011] C. Yuan, R. O. Fox,  Conditional quadrature method of moments for kinetic equations,  J. Comp. Phys. 230 (2011) 8216–8246.
  * Fox et al. [2018] R. O. Fox, F. Laurent, A. Vié,  Conditional hyperbolic quadrature method of moments for kinetic equations,  J. Comp. Phys. 365 (2018) 269–293.
  * Patel et al. [2019] R. G. Patel, O. Desjardins, R. O. Fox,  Three-dimensional conditional hyperbolic quadrature method of moments,  J. Comp. Phys. X 1 (2019) 100006.
  * Fox and Laurent [2021] R. Fox, F. Laurent,  Hyperbolic quadrature method of moments for the one-dimensional kinetic equation,  arXiv:2103.10138 (2021).
  * Menikoff and Plohr [1989] R. Menikoff, B. J. Plohr,  The Riemann problem for fluid-flow of real materials,  Rev. Mod. Phys. 61 (1989) 75–130.
  * Maeda and Colonius [2018] K. Maeda, T. Colonius,  Eulerian–Lagrangian method for simulation of cloud cavitation,  J. Comp. Phys. 371 (2018) 994–1017.
  * Bryngelson et al. [2020] S. H. Bryngelson, T. Colonius, R. O. Fox,  QBMMlib: A library of quadrature-based moment methods,  SoftwareX 12 (2020) 100615.
  * Colonius et al. [2008] T. Colonius, R. Hagmeijer, K. Ando, C. E. Brennen,  Statistical equilibrium of bubble oscillations in dilute bubbly flows,  Phys. Fluids 20 (2008).
  * Bryngelson et al. [2020a] S. H. Bryngelson, A. Charalampopoulos, T. P. Sapsis, T. Colonius,  A Gaussian moment method and its augmentation via LSTM recurrent neural networks for the statistics of cavitating bubble populations,  Int. J. Mult. Flow 127 (2020a) 103262.
  * Bryngelson et al. [2020b] S. H. Bryngelson, K. Schmidmayer, V. Coralic, J. C. Meng, K. Maeda, T. Colonius,  MFC: An open-source high-order multi-component, multi-phase, and multi-scale compressible flow solver,  Comp. Phys. Comm. (2020b) 107396.
  * Jiang and Shu [1996] G.-S. Jiang, C.-W. Shu,  Efficient implementation of weighted eno schemes,  J. Comp. Phys. 126 (1996) 202–228.
  * Toro et al. [1994] E. Toro, M. Spruce, W. Speares,  Restoration of the contact surface in the HLL-Riemann solver,  Shock waves 4 (1994) 25–34.
  * Gottlieb et al. [2001] S. Gottlieb, C. W. Shu, E. Tadmor,  Strong stability-preserving high-order time discretization methods,  SIAM Rev. (2001).
  * Gottlieb and Shu [1998] S. Gottlieb, C.-W. Shu,  Total variation diminishing Runge–Kutta schemes,  Math. Comput. 67 (1998) 73–85.
  * Thompson et al. [1986] P. O. Thompson, W. C. Cummings, S. J. Ha,  Sounds, source levels, and associated behavior of humpback whales, Southeast Alaska,  J. Acoustic. Soc. Am. 80 (1986) 735–740.
  * Charalampopoulos et al. [2021] A.-T. Charalampopoulos, S. H. Bryngelson, T. Colonius, T. P. Sapsis,  Hybrid quadrature moment method for accurate and stable representation of non-Gaussian processes and their dynamics,  arXiv:2110.01374 (2021).
  * Wheeler [1974] J. C. Wheeler,  Modified Moments and Gaussian Quadratures,  Rocky Mt. J. Math. 4 (1974) 287–296.

##  Appendix A CHyQMOM moment inversion

This appendix describes the 2D, $2\times 2$ (4) node CHyQMOM algorithm of
[34]. We note that Patel et al. [35] provides the 3D version of this
algorithm, though it is not necessary for our bubble dynamics problem.
Algorithm 3 is the full 2D CHyQMOM algorithm, which references the 1D 2-node
HyQMOM algorithm 4. The optimal moment set for this case is

$\displaystyle\vec{\boldsymbol{\mu}}=\left\\{\mu_{00},\mu_{10},\mu_{01},\mu_{20},\mu_{02},\mu_{11}\right\\},$
(24)

which is also the input to algorithm 3.

1:procedure CHyQMOM4($\boldsymbol{\mu}$)

2: $D_{ij}=\mu_{ij}/\mu_{00}$

3: $C_{20}\leftarrow D_{20}-D_{10}^{2}$

4: $C_{02}\leftarrow D_{02}-D_{01}^{2}$

5: $C_{11}\leftarrow D_{11}-D_{10}D_{01}$

6: $\boldsymbol{\rho},\bm{x}^{\prime}=\textsc{HyQMOM2}(\\{1,0,C_{20}\\})$

7: $\bm{y}^{\prime}=\bm{x}^{\prime}C_{11}/C_{20}$

8: $\mu_{\mathrm{avg}}^{2}=C_{02}-\sum_{j}\rho_{j}{y_{j}^{\prime}}^{2}$

9:
$\widetilde{\boldsymbol{\rho}},\widetilde{\bm{x}}^{\prime}=\textsc{HyQMOM2}(\\{1,0,\mu_{\mathrm{avg}}^{2}\\})$

10: for $i,j\in[1,2]$ do

11: $w_{ij}\leftarrow\mu_{00}\rho_{i}\widetilde{\rho}_{j}$

12: $x_{ij}\leftarrow D_{10}+x_{j}^{\prime}$

13: $y_{ij}\leftarrow D_{01}+y_{j}^{\prime}+\widetilde{x}_{i}^{\prime}$

14: end for

15:return $\bm{w},\bm{x},\bm{y}$

16:end procedure

Algorithm 3 CHyQMOM $2\times 2$. $\boldsymbol{\mu}$ are the input moments, $D$
are normalized moments, $C$ are central moments, and $\bm{w}$, $\bm{x}$ and
$\bm{y}$ are the weights and node locations in the first and second coordinate
directions (corresponding to $R$ and $\dot{R}$ in the main text).

1:procedure HyQMOM2($\boldsymbol{\mu}$)

2: $C_{2}\leftarrow(\mu_{0}\mu_{2}-\mu_{1}^{2})/\mu_{0}^{2}$

3: $w_{1}=w_{2}\leftarrow\mu_{0}/2$

4: $x_{1}\leftarrow\mu_{1}/\mu_{0}+\sqrt{C_{2}}$

5: $x_{2}\leftarrow\mu_{1}/\mu_{0}-\sqrt{C_{2}}$

6:return $\bm{w},\bm{x}$

7:end procedure

Algorithm 4 The two-node 1D HyQMOM algorithm. Nomenclature follows algorithm
3, with $\bm{w}$ and $\bm{x}$ serving as dummy weights and node locations.

##  Appendix B CQMOM moment inversion

This appendix has the same form as appendix A, though it instead describes the
2D CQMOM algorithm of [33]. Algorithm B is the full 2D CQMOM algorithm,
referencing Wheeler’s method for computing optimal 1D quadrature nodes and
weights 6. The optimal moment set for this case is

$\displaystyle\vec{\boldsymbol{\mu}}=\left\\{\mu_{00},\mu_{10},\mu_{01},\mu_{20},\mu_{02},\mu_{11},\mu_{30},\mu_{03},\mu_{12},\mu_{13}\right\\},$
(25)

which is also the input to algorithm 5.

1:procedure CQMOM($\boldsymbol{\mu}$)

2:
$\bm{w}^{(x)},\bm{x}\leftarrow\textsc{Wheeler}(\\{\mu_{00},\mu_{10},\dots,\mu_{2N_{x}-1,0}\\})$

3: for $i,j\in[1,\dots,N_{x}]$ do

4: $V_{ij}=x_{j}^{i-1}$

5: $W_{ii}=w_{i}^{(x)}$

6: end for

7: for $i\in[0,\dots,2N_{y}-1]$ do

8:
$\bm{x}_{i}^{\prime}\leftarrow\bm{V}\bm{W}\,\bm{x}_{i}^{\prime}=\mu_{0:N_{x}-1,i}$

9: $w_{ij}^{(y)},y_{ij}\leftarrow\textsc{Wheeler}(\bm{x}_{i}^{\prime})$

10: end for

11: $w_{ij}=w_{i}^{(x)}w_{ij}^{(y)}$

12:return $\bm{w},\bm{x},\bm{y}$

13:end procedure

Algorithm 5 2D CQMOM algorithm with dummy coordinate direction $y$ conditioned
on $x$. $\bm{V}$ is a Vandermonde matrix and $\bm{W}$ is a diagonal weight
matrix.

1:procedure Wheeler($\boldsymbol{\mu}$)

2: for $i\in[1,\dots,2N]$ do

3: $\sigma_{1,i}\leftarrow\mu_{i-1}$

4: end for

5: $a_{0}\leftarrow\mu_{1}/\mu_{0}$

6: $b_{0}\leftarrow 0$

7: for $i\in[2,\dots,N]$ do

8: for $j\in[i,\dots,2N-i+1]$ do

9: $\sigma_{i,j}\leftarrow\sigma_{i-1,j+1}-a_{i-2}\sigma_{i-1,i-1}$

10: end for

11:
$a_{i-1}\leftarrow\sigma_{i,i+1}/\sigma_{i,i}-\sigma_{i-1,i}/\sigma_{i-1,i-1}$

12: $b_{i-1}\leftarrow\sigma_{i,i}/\sigma_{i-1,i-1}$

13: end for

14: for $i\in[2,\dots,N-1]$ do

15: $J_{i,i}\leftarrow a_{i-1}$

16: $J_{i+1,i}=J_{i,i+1}\leftarrow\sqrt{b_{i}}$

17: end for

18: $\bm{Q}\boldsymbol{\Lambda}\bm{Q}^{-1}\leftarrow\textsc{EVD}(\bm{J})$

19: for $i\in[1,\dots,N]$ do

20: $x_{i}\leftarrow\Lambda_{ii}$

21: $w_{i}\leftarrow\mu_{0}Q_{i,1}^{2}$

22: end for

23:return $\bm{w},\bm{x}$

24:end procedure

Algorithm 6 Wheeler’s algorithm for optimal quadrature points $\bm{x}$ and
weights $\bm{w}$ given a moment set $\boldsymbol{\mu}$ [49, 32]. $\bm{J}$ is a
Jacobi matrix, and EVD refers to the eigenvalue decomposition.


# How Multilingual Are Large Language Models
Fine-Tuned for Translation?

Aquia Richburg
AMSC
University of Maryland
College Park, MD 20742, USA
<EMAIL_ADDRESS>
&Marine Carpuat
Computer Science
University of Maryland
College Park, MD 20742, USA
<EMAIL_ADDRESS>
###### Abstract

A new paradigm for machine translation has recently emerged: fine-tuning large
language models (Llm) on parallel text has been shown to outperform dedicated
translation systems trained in a supervised fashion on much larger amounts of
parallel data (Xu et al., 2024a; Alves et al., 2024). However, it remains
unclear whether this paradigm can enable massively multilingual machine
translation or whether it requires fine-tuning dedicated models for a small
number of language pairs. How does translation fine-tuning impact the MT
capabilities of Llms for zero-shot languages, zero-shot language pairs, and
translation tasks that do not involve English? To address these questions, we
conduct an extensive empirical evaluation of the translation quality of the
Tower family of language models (Alves et al., 2024) on 132 translation tasks
from the multi-parallel Flores-200 data. We find that translation fine-tuning
improves translation quality even for zero-shot languages on average, but that
the impact is uneven depending on the language pairs involved. These results
call for further research to effectively enable massively multilingual
translation with Llms.

##  1 Introduction & Background

Machine translation (MT) systems have long been built for a specific pair of
input and output languages, and trained on a parallel corpus representing this
language pair and translation direction (Brown et al., 1990). End-to-end
neural models made it easier to train a single model on multiple language
pairs, enabling multilingual machine translation models (Johnson et al., 2017;
Arivazhagan et al., 2019b; Dabre et al., 2020; Team et al., 2022), where low-
resource languages benefit from transfer learning from high-resource
languages, including in zero-shot translation directions. While
transformative, this paradigm still assumes access to large amounts of
parallel data which is unrealistic for many of the world’s languages – less
than 5% of languages have labeled data useful enough for applications (Joshi
et al., 2020).

Over the past year, large language models (Llms) have started to exhibit
competitive translation quality compared to dedicated supervised MT systems
(Kocmi et al., 2023). While Llms have been used for translation in zero-shot
and few-shot settings (Vilar et al., 2023; Hendy et al., 2023; Zhang et al.,
2023a; Bawden & Yvon, 2023; Lin et al., 2022; Zhu et al., 2023, among others),
the performance of Llms often lagged behind dedicated MT systems trained in a
supervised fashion on parallel text, particularly for low-resource languages
Instruction fine-tuning of a multilingual Llm for MT has emerged as a strong
contender to build state-of-the-art systems by providing a straightforward
mechanism to combine the strengths of multilingual pre-training on primarily
monolingual data with the controlled injection of parallel data – rather than
only relying on incidental bilingualism Briakou et al. (2023). Xu et al.
(2024a) showed that fine-tuning LLaMA2 models in two stages – first on
monolingual data followed by small amounts of high-quality parallel data –
could outperform dedicated supervised MT systems (Team et al., 2022) and large
closed models, like GPT3.5, with 7B or 13B parameter open LLaMA2 models.
Translation quality has been further improved with contrastive preference
optimization (Xu et al., 2024b). Meanwhile, Alves et al. (2024) recently
introduced a family of Llm for translation-related tasks which is also based
on LLaMA2. It involves continued pre-training on monolingual and parallel
text, as well as fine-tuning on parallel text (and related tasks such as post-
editing). The resulting models exhibit impressive machine translation
performance, outperforming dedicated multilingual models (Team et al., 2022)
and Xu et al. (2024b)’s best models, when translating between English and the
non-English languages covered in the training set. In the same vein, recent
models such as BigTranslate (Yang et al., 2023) or BayLing (Zhang et al.,
2023b) employ continued pre-training on parallel data or translation
instruction fine-tuning to improve the translation capabilities of Llms.

However, it remains unclear how well Llms fine-tuned for MT translate
languages that they have not been trained on. Multilingual Llms are imbalanced
in their language coverage and are affected by the curse of multilinguality
during pre-training (Chang et al., 2023; Liang et al., 2023); this is
compounded by the fact that exposure to fine-tuning parallel data is
necessarily limited to a small number of language pairs raising the same
questions about maximizing the benefits of transfer learning for low-resource
languages without degrading performance for high-resource languages that have
been extensively studied in multilingual MT (Arivazhagan et al., 2019a; Liu &
Niehues, 2022; Chang et al., 2023; Yuan et al., 2023, among many others).

To fill this gap, this work seeks to characterize how Llms fine-tuned for
translation perform across a diverse set of translation tasks beyond the ones
they were fine-tuned on. Addressing this question is crucial to determine
whether using Llms for MT will require fine-tuning models for individual (or
small sets of) language pairs, or whether they can truly fulfill the promise
of multilingual Llm and enable massively multilingual MT. We present an
extensive empirical evaluation of the Tower models on 132 machine translation
tasks with varying degrees of supervision. Inspired by Choudhury & Deshpande
(2021), who critique the use of average performance over a set of languages as
the sole model selection rationale for massively multilingual Llm, we
contribute an experimental design and detailed analysis of results which aims
to balance the need to cover an extensive set of languages to understand the
generalization ability of the models, with the desire to understand
performance beyond a single score which averages potentially disparate
behaviors across languages. Our selection of languages is driven by both
linguistic typology and data size criteria. Furthermore, we include
translation between all language pairs, and do not limit ourselves to
translation into and out of English. While MT research naturally does not
focus on English as much as other Natural Language Processing research
(Bender, 2019; Joshi et al., 2020), the vast majority of MT evaluations
involve English as either the source or target language. Here we also seek to
assess to what degree systems trained on English-centric parallel data
generalize to translation tasks that do not involve English.111This goal is
shared with Gao et al. (2024), however, their evaluation is limited to
language pairs that involve languages seen at fine-tuning time, while we also
consider generalization to zero-shot languages in addition to zero-shot
language pairs.

To be explicit, this work does not introduce any new modeling techniques.
Instead, we contribute:

  * •

an extensive empirical evaluation of the generalization ability of Llms fine-
tuned for translation, using Tower models on MT tasks involving 132
translation directions between 12 languages, including German (De), English
(En), Korean (Ko), Dutch (Nl), Russian (Ru) and Chinese (Zh) – the supervised
languages seen during fine-tuning – and Czech (Cs), Icelandic (Is), Japanese
(Ja), Polish (Pl), Swedish (Sv) and Ukrainian (Uk) – the zero-shot languages
not seen during fine-tuning (Section 2).

  * •

a fine-grained analysis of results that goes beyond averaging scores across
languages and considers best and worst behavior per model and per source and
target language to provide a more complete picture of model performance across
languages(Choudhury & Deshpande, 2021) (Section 3).

  * •

findings that indicate that fine-tuning on translation-related tasks improve
Llm’s ability to handle the task of translation itself beyond the specific
language pairs seen during fine-tuning (Section 3).

  * •

however, the fine-grained analysis reveals an uneven picture with higher
variance in translation quality for language pairs that are not fully
supervised, and low quality in outlier languages whether they are seen during
fine-tuning (Korean) or not (Icelandic) (Section 3).

##  2 Evaluation Design

### Evaluation Languages

Since we are interested in measuring the supervised and zero-shot translation
capabilities of Llms before and after fine-tuning, we use a subset of the
fine-tuning languages from Tower as our supervised set and a collection of
related languages seen as some proportion of LLaMA2’s pre-training data as our
zero-shot set, as summarized in Table 1. The supervised language set consists
of German (De), English (En), Korean (Ko), Dutch (Nl), Russian (Ru) and
Chinese (Zh). The zero-shot language set consists of Czech (Cs), Icelandic
(Is), Japanese (Ja), Polish (Pl), Swedish (Sv) and Ukrainian (Uk). The
languages in the zero-shot set are chosen to represent a range in resource
support and to be related to languages in the supervised set by language
family, typological properties or orthography. This yields $\binom{12}{2}$ =
66 pairs from the set of languages times, and thus 132 translation tasks when
considering both translation directions. We divide these 132 tasks into 3
subcategories based on the amount of supervision available in Tower:

  * •

Fully supervised translation directions: These consist of De, Ko, Nl, Ru and
Zh to and from En. This totals 10 pairs.

  * •

Partially supervised translation directions: These are the pairs where at
least one language is from the supervised set, but the exact language pair was
not seen during fine-tuning. For example, De-Zh belongs to this set because no
direct De-Zh parallel data was used during fine-tuning. This totals 92 pairs.

  * •

Unsupervised translation directions: These are the pairs where both languages
are in the zero-shot set. This totals 30 pairs.

Language | Script | LLaMA-2 support | Similarity Groups
---|---|---|---
Czech | Latin | 0.03% | West Slavic
Polish | Latin | 0.09% | West Slavic
Russian | Cyrillic | 0.13% | East Slavic
Ukrainian | Cyrillic | 0.07% | East Slavic
German | Latin | 0.17% | West Germanic
English | Latin | 89.70% | West Germanic
Icelandic | Latin | possibly unseen | North Germanic
Dutch | Latin | 0.12% | West Germanic
Swedish | Latin | 0.15% | North Germanic
Japanese | Kana | 0.10% | Kanji from Hanzi, SOV order
Korean | Hangul | 0.06% | SOV order
Chinese | Hanzi | 0.13% | Hanzi to Kanji, loanwords to Ja and Ko
Table 1: Evaluated languages with rationales for similarity grouping. The
languages marked in bold belong to the supervised set for both Tower models

### Test Data & Metrics

For all experiments, we evaluate on the multilingual, multi-parallel
Flores-200 (Team et al., 2022) devtest set, which consists of English language
articles manually translated into multiple languages. We report reference-
based Comet-22 (Rei et al., 2022a) as our main metric for translation quality.
We also computed the commonly used Bleu (Papineni et al., 2002) metric,
however in addition to its well-documented weaknesses, it exhibits a higher
variance over diverse languages due to the differences in tokenization.

###  Llms Translation Models

Our zero-shot baseline model is MetaAI’s LLaMA2 (Touvron et al., 2023), an LLM
which has been shown to be a strong baseline on a number of English language
tasks. For fine-tuning, we use the Tower (Alves et al., 2024) family of
models, which fine-tunes LLaMA2 for a range of translation-related tasks.
TowerBase uses a two-third/one-third mixture of monolingual and parallel data
to continue pre-training LLaMA2. The monolingual data is sourced from mC4 (Xue
et al., 2021), a multilingual webcrawl corpus, and the parallel data comes
from several public sources ranging in domain from news to medical to
Wikimedia. TowerInstruct further optimizes TowerBase for translation through
instruction fine-tuning, formatting the translation task as natural language
instructions. A diverse collection of data resources are formulated into
question-answer templates including some zero-shot directions. TowerInstruct
expands with task diversity to include paraphrasing, dialog data and coding
instructions. We consider both the 7B and 13B parameter model variants for
each of the Llms. At inference time, for the Tower models, we use the same
prompt format as presented during fine-tuning. Since LLaMA2 was not trained
for translation, we default to a basic natural language prompt to translate
the source sent from Language A to Language B. We use beam search with a beam-
width of 5 and a maximum of 128 new tokens throughout.

### Dedicated MT system

We also include the state-of-the-art No Language Left Behind (Nllb) model Team
et al. (2022).222We use the distilled 600M to maximize inference speed. On a
more limited set of language pairs, Alves et al. (2024) report consistent
results with the larger 54B variant. Nllb is a supervised multilingual NMT
system trained to translate between a diverse variety of 200 languages,
including using parallel data that does not include English. Among the
languages we evaluated, all language pairs that involve English and Korean
$\leftrightarrow$ Japanese are seen at training time. At inference time, we
use language specific tokens marking the source and target consistent with
Nllb training and use the same beam search setting as for the Llms.

### Limitations

While our empirical study is extensive, it comes with inherent limitations.
Our use of the Flores-200 test data means that we use text that was originally
written in English and translated into other languages, no matter what MT
direction is evaluated, which might introduce translationese effects (Team et
al., 2022). We accept this limitation as the multi-way parallel nature of
Flores-200 has the advantage of letting us make more controlled comparisons
across languages. Our evaluation is based on automatic metrics that are known
to correlate highly with human judgments Rei et al. (2022b), it remains to be
seen how the translation quality differences we observe are perceived by
people.

##  3 Results

Figure 1: Distribution of Comet scores for each translation approach over
language pair supervision type: TowerInstruct-13B models are competitive with
Nllb on average across settings, suggesting that transfer learning benefits
many unsupervised languages, but not all based on the increased variance in
the least supervised conditions.

### Impact of Supervision Type

In Figure 1, we plot the distribution of Comet scores for each model according
to translation supervision type. The dedicated Nllb MT model represent a
strong baseline across supervision types, with mean Comet scores of 0.8 or
higher. As can be expected, the Llm systems reach higher mean Comet scores for
the fully supervised language pairs than others. When comparing Llm with the
dedicated Nllb system, we find that the zero-shot LLaMA2 models lag behind by
a wide margin. Interestingly, the smaller LLaMA2-7B model yields higher mean
Comet scores than the larger LLaMA2-13B model. For supervised language pairs,
all Tower models are on par or better than Nllb. However, their performance is
more uneven in the other supervision conditions: continued pre-training in the
TowerBase models improves Comet scores compared to the zero-shot LLaMA2
results but the mean Comet scores lag behind those of the dedicated Nllb
systems. Translation-related fine-tuning in TowerInstruct models improves
translation quality further, but only the larger TowerInstruct-13B yields
higher mean Comet scores than Nllb. These results suggest the translation
supervision at pre-training and fine-tuning time does improve Llm’s
translation abilities beyond the languages seen in parallel text, and thus
that transfer learning benefits other languages.

However, the mean scores only provide a partial picture of the translation
quality across the 132 language pairs considered. Figure 1 clearly shows, for
all models, that the variance of Comet scores across languages also increases
in the partial supervision and unsupervised settings compared to the fully
supervised settings, which motivates us to look at results in a finer-grained
fashion. The variance increases more for Llm than Nllb models indicating that
the dedicated Nllb models have a more reliable behavior across language pairs.
TowerInstruct has a large tail of low quality outliers, some of which are
lower than the mean zero-shot LLaMA2 scores. This shows that translation-
related fine-tuning does not benefit all language pairs uniformly.

Figure 2: Spread of translation quality as measured by Comet score for each
target language for the dedicated Nllb model (top) and the overall top-
performing TowerInstruct-13B (bottom). For each data point, the color
identifies the source language, and we mark the best/worst dissimilar paths
(D+/D-), the similar path (S), and translation from English (E).

### Best/Worst Translation Paths per Model

To dig deeper, we examine differences in translation quality across languages
for fixed models. Figure 2 represents the spread of COMET scores for each
target language for the dedicated Nllb model (top) and the overall top-
performing TowerInstruct-13B (bottom).333The complete set of graphs for all
models and source/target language groupings can be found in Appendix C.
Results show that while on average TowerInstruct-13B is the top performing
model, its worst performance is consistently worse than that of Nllb, for both
supervised and unsupervised target languages. Translation from Icelandic
almost always yields the worst COMET score no matter the target language (red
points). The best performing dissimilar translation paths (D+ points)
interestingly tend to overlap with other top performing languages, indicating
that dissimilarity between source and target language is not necessarily a
challenge for TowerInstruct-13B.

Figure 3: Models that achieve the absolute best and worst scores per language
pair. The group outlined in dotted-blue corresponds to the fully supervised
set, the group outlined in dot-dash-red corresponds to the unsupervised set
and the remaining pairs are partially supervised. The circle and cross symbols
denote whether the language pair consists of similar or dissimilar language,
respectively (as defined in Table 1). Symbols marked with asterisks denote
statistical significance according to a paired $t$-test at a significance
level of 0.05.

### Best/Worst Models per Language Pair

Next, we show the models with the absolute best and worst translation quality
per language pair in Figure 3. TowerInstruct-13B emerges as the top performing
model across language pairs (left graph), while the zero-shot LLaMA2 models
often yield the worst performance (right graph); however there are interesting
exceptions in both cases.

Among the best performing systems, we see that for translation into English
and for supervised languages, the largest TowerInstruct-13B is not necessarily
best: the smaller TowerInstruct-7B yields higher COMET scores for translation
from 5 language pairs into English (top graph, top row), and surprisingly the
TowerBase-13B model outperforms the fine-tuned TowerInstruct models for
translation from German into English. So bigger models with more supervision
are not uniformly better. Furthermore, Llm all struggle compared to Nllb with
translation involving Icelandic and Korean. When Korean is the source
TowerInstruct wins 9 out of 11 times but loses completely when Korean is
switched to the target. In the case of Icelandic, TowerInstruct loses in both
directions. Icelandic is not seen at fine-tuning time, accounts for 1% at most
of LLaMA2 pre-training data, and is a fusional language, so it is perhaps
unsurprising that is challenging to translate for the Llms. The Korean results
are more surprising since it is included in fine-tuning. However, its
typography (unique Hangul) and grammar (Subject-Object-Verb) differ vastly
from the other included supervised languages (mainly Subject-Verb-Object)
which might make it harder to model on the target side and explain why Nllb
remains stronger. This suggests that when selecting a set of fine-tuning
languages, more typological diversity might support better generalization to
new languages.

Interestingly, the worst performing systems are often the larger LLaMA2-13B
zero-shot models, except for translation into Korean and Icelandic where the
smaller LLaMA2-7B zero-shot models are worse. This suggests that the increased
capacity in the underlying model benefits the representation of these outlier
languages, but not for the more similar higher-resource languages.
Additionally, for a small number of unsupervised language pairs (English-
Swedish, German-Ukrainian, Dutch-Ukrainian and Dutch-Swedish), continued pre-
training hurts translation quality and TowerBase-7B models underperform the
zero-shot LLaMA2 models. Similar versus dissimilar translation directions do
not show obviously different patterns.

To understand these results better, we turn to examining the properties of the
inputs and outputs of Nllb versus Llm translations.

### Off-Target Outputs

A known failure mode of Llm is to generate translations in the wrong target
language: for instance, Xu et al. (2024a) note that this behavior for LLaMA2
when translating out of English. One possible explanation is that when faced
with translation on languages with lower support the model defaults to a
higher-resource language, usually English. We report the average percentage of
translations that are found not to match the specified target language for our
Llm’s using the lingua-py444https://github.com/pemistahl/lingua-py language
identification tool in Figure 4. The analysis of off-target translation across
the multiple language paths considered here are consistent with prior work,
zero-shot LLaMA2 models are often off target. Interestingly the larger
LLaMA2-13B is more often off-target than the smaller LLaMA2-7B model, which
contributes to its worse overall performance (Figure 1). Continued pre-
training and fine-tuning substantially improve the on-target behavior of
supervised languages (Figure 4 left), however off-target outputs remain a
problem even after fine-tuning when translating into unsupervised languages
(Figure 4 right). These results suggest that Llm translations into Korean and
Icelandic are bad for different reasons: translation into Icelandic is only on
target about 40% of the time, while translation into Korean is generally on
target (more than 90% of the time).

Figure 4: Ratio of on-target translation averaged over the source sentences
from the supervised set. The TowerInstruct models are generally on target for
supervised languages, but off-target translations are frequent for
unsupervised languages. By contrast, Nllb translations are on target for
97-100% of outputs.

### Input Tokenization

Based on the typological property of each language and its pre-training
representation, parallel texts can be encoded into subword sequences of vastly
different lengths across languages, leading to discrepancies in the quality
and cost of inference across languages Limisiewicz et al. (2023); Ahia et al.
(2023); Petrov et al. (2023). The tokenizers for Nllb and our Llm models are
all built using SentencePiece (Kudo & Richardson, 2018). Nllb trains a 256k
vocabulary model by temperature sampling 100M bitext pairs to downsample
higher resource languages and upsample lower resource languages (Team et al.,
2022). The LLaMA2 tokenizer (also used by the Tower models) has a vocabulary
size of 32k trained on the mostly English data mixture (Touvron et al., 2023).
Figure 5 plots the average number of source subwords per segment (top) in
parallel with the average COMET scores for all models for translation out of
each source language (bottom).555We provide the symmetric figure for output
length in Appendix Figure 21. The Nllb and LLaMA2 tokenizers both favor
English, with inputs in all languages being longer than English indicating
over-segmentation. The Llm tokenizer over-tokenizes at a higher rate than the
Nllb tokenizer, and has wider discrepancy across languages: Korean text is
more than 3 times longer than when translating into English, which might be a
factor in the lower translation quality of TowerInstruct-13B into Korean
compared to other language pairs, despite its supervised status. Icelandic
inputs are more than twice longer than English, but are somewhat shorter than
Japanese inputs which are better translated by the Llm models. Input length
thus does not entirely explain translation quality.

Figure 5: Average input length per language after subword segmentation for
each model (top) and average translation quality out of each source language
per model (bottom).

##  4 Conclusion

This work sought to characterize how Llm fine-tuned for translation perform
across a diverse set of translation tasks beyond the ones they were fine-tuned
on. By conducting an extensive empirical evaluation of the Tower family of
Llms on 132 translation tasks representing diverse degrees of supervision, we
find that fine-tuning on translation-related tasks improves Llm’s ability to
handle the task of translation itself beyond the specific language pairs seen
during fine-tuning. This encouragingly suggests that the fine-tuning paradigm
has the potential to enable massively multilingual MT.

However, the fine-grained results show that transfer learning has an uneven
impact: outlier languages remain hard to handle whether they are seen during
fine-tuning (Korean) or not (Icelandic). Analysis suggest that Korean might be
harmed by oversegmentation compared to other fine-tuning languages, while
translation involving Icelandic often results in off-target outputs.
Furthermore, the worst-case behavior of the best Llm (TowerInstruct-13B) is
consistently worse than that of the dedicated Nllb model. This highlights the
need for future work on improving instruction tuning techniques to benefit the
harder language pairs. Initial efforts in that direction are underway. For
instance, Gao et al. (2024) have ported a cross-lingual consistency
regularization method developed for dedicated multilingual MT (Gao et al.,
2023) to instruction fine-tuning – while promising when tested on languages
seen during fine-tuning, its impact on translation involving languages unseen
during fine-tuning remains to be seen. Our results further suggest that
ensembling methods that back-off to dedicated MT models might be useful to
mitigate off-target translations, and also call for addressing tokenization
fairness as a potential cause for discrepancies in translation quality across
languages.

While we have studied overall translation quality as measured by Comet, other
aspects of the outputs would be worth studying in future work including how
translationese effects (Dutta Chowdhury et al., 2020; Vanmassenhove et al.,
2021) and hallucination patterns (Guerreiro et al., 2023) are impacted by
different degrees of instruction-tuning supervision.

## References

  * Ahia et al. (2023) Orevaoghene Ahia, Sachin Kumar, Hila Gonen, Jungo Kasai, David R. Mortensen, Noah A. Smith, and Yulia Tsvetkov.  Do All Languages Cost the Same? Tokenization in the Era of Commercial Language Models, May 2023.
  * Alves et al. (2024) Duarte M. Alves, José Pombal, Nuno M. Guerreiro, Pedro H. Martins, João Alves, Amin Farajian, Ben Peters, Ricardo Rei, Patrick Fernandes, Sweta Agrawal, Pierre Colombo, José G. C. de Souza, and André F. T. Martins.  Tower: An open multilingual large language model for translation-related tasks, 2024.
  * Arivazhagan et al. (2019a) Naveen Arivazhagan, Ankur Bapna, Orhan Firat, Roee Aharoni, Melvin Johnson, and Wolfgang Macherey.  The Missing Ingredient in Zero-Shot Neural Machine Translation, March 2019a.
  * Arivazhagan et al. (2019b) Naveen Arivazhagan, Ankur Bapna, Orhan Firat, Dmitry Lepikhin, Melvin Johnson, Maxim Krikun, Mia Xu Chen, Yuan Cao, George Foster, Colin Cherry, Wolfgang Macherey, Zhifeng Chen, and Yonghui Wu.  Massively Multilingual Neural Machine Translation in the Wild: Findings and Challenges, July 2019b.
  * Bawden & Yvon (2023) Rachel Bawden and François Yvon.  Investigating the Translation Performance of a Large Multilingual Language Model: The Case of BLOOM.  In Mary Nurminen, Judith Brenner, Maarit Koponen, Sirkku Latomaa, Mikhail Mikhailov, Frederike Schierl, Tharindu Ranasinghe, Eva Vanmassenhove, Sergi Alvarez Vidal, Nora Aranberri, Mara Nunziatini, Carla Parra Escartín, Mikel Forcada, Maja Popovic, Carolina Scarton, and Helena Moniz (eds.), _Proceedings of the 24th Annual Conference of the European Association for Machine Translation_ , pp. 157–170, Tampere, Finland, June 2023. European Association for Machine Translation.
  * Bender (2019) Emily M. Bender.  The #BenderRule: On Naming the Languages We Study and Why It Matters.  https://thegradient.pub/the-benderrule-on-naming-the-languages-we-study-and-why-it-matters/, September 2019.
  * Briakou et al. (2023) Eleftheria Briakou, Colin Cherry, and George Foster.  Searching for Needles in a Haystack: On the Role of Incidental Bilingualism in PaLM’s Translation Capability, May 2023.
  * Brown et al. (1990) Peter Brown, John Cocke, Stephen Della Pietra, Vincent Della Pietra, Frederik Jelinek, John Lafferty, Robert Mercer, and Paul Rossin.  A statistical approach to machine translation.  _Computational Linguistics_ , 16:79–85, 1990.
  * Chang et al. (2023) Tyler A. Chang, Catherine Arnett, Zhuowen Tu, and Benjamin K. Bergen.  When is multilinguality a curse? language modeling for 250 high- and low-resource languages, 2023.
  * Choudhury & Deshpande (2021) Monojit Choudhury and Amit Deshpande.  How Linguistically Fair Are Multilingual Pre-Trained Language Models?  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 35(14):12710–12718, May 2021.  ISSN 2374-3468.  doi: 10.1609/aaai.v35i14.17505.
  * Dabre et al. (2020) Raj Dabre, Chenhui Chu, and Anoop Kunchukuttan.  A Survey of Multilingual Neural Machine Translation.  _ACM Computing Surveys_ , 53(5):99:1–99:38, September 2020.  ISSN 0360-0300.  doi: 10.1145/3406095.
  * Dutta Chowdhury et al. (2020) Koel Dutta Chowdhury, Cristina España-Bonet, and Josef van Genabith.  Understanding Translationese in Multi-view Embedding Spaces.  In _Proceedings of the 28th International Conference on Computational Linguistics_ , pp. 6056–6062, Barcelona, Spain (Online), December 2020. International Committee on Computational Linguistics.  doi: 10.18653/v1/2020.coling-main.532.
  * Gao et al. (2023) Pengzhi Gao, Liwen Zhang, Zhongjun He, Hua Wu, and Haifeng Wang.  Improving Zero-shot Multilingual Neural Machine Translation by Leveraging Cross-lingual Consistency Regularization.  In Anna Rogers, Jordan Boyd-Graber, and Naoaki Okazaki (eds.), _Findings of the Association for Computational Linguistics: ACL 2023_ , pp. 12103–12119, Toronto, Canada, July 2023. Association for Computational Linguistics.  doi: 10.18653/v1/2023.findings-acl.766.
  * Gao et al. (2024) Pengzhi Gao, Zhongjun He, Hua Wu, and Haifeng Wang.  Towards Boosting Many-to-Many Multilingual Machine Translation with Large Language Models.  https://arxiv.org/abs/2401.05861v1, January 2024.
  * Guerreiro et al. (2023) Nuno M. Guerreiro, Duarte Alves, Jonas Waldendorf, Barry Haddow, Alexandra Birch, Pierre Colombo, and André F. T. Martins.  Hallucinations in Large Multilingual Translation Models, March 2023.
  * Hendy et al. (2023) Amr Hendy, Mohamed Abdelrehim, Amr Sharaf, Vikas Raunak, Mohamed Gabr, Hitokazu Matsushita, Young Jin Kim, Mohamed Afify, and Hany Hassan Awadalla.  How Good Are GPT Models at Machine Translation? A Comprehensive Evaluation, February 2023.
  * Johnson et al. (2017) Melvin Johnson, Mike Schuster, Quoc V. Le, Maxim Krikun, Yonghui Wu, Zhifeng Chen, Nikhil Thorat, Fernanda Viégas, Martin Wattenberg, Greg Corrado, Macduff Hughes, and Jeffrey Dean.  Google’s Multilingual Neural Machine Translation System: Enabling Zero-Shot Translation.  _Transactions of the Association for Computational Linguistics_ , 5:339–351, October 2017.  ISSN 2307-387X.  doi: 10.1162/tacl˙a˙00065.
  * Joshi et al. (2020) Pratik Joshi, Sebastin Santy, Amar Budhiraja, Kalika Bali, and Monojit Choudhury.  The State and Fate of Linguistic Diversity and Inclusion in the NLP World.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_ , pp. 6282–6293, Online, July 2020. Association for Computational Linguistics.  doi: 10.18653/v1/2020.acl-main.560.
  * Kocmi et al. (2023) Tom Kocmi, Eleftherios Avramidis, Rachel Bawden, Ondřej Bojar, Anton Dvorkovich, Christian Federmann, Mark Fishel, Markus Freitag, Thamme Gowda, Roman Grundkiewicz, Barry Haddow, Philipp Koehn, Benjamin Marie, Christof Monz, Makoto Morishita, Kenton Murray, Makoto Nagata, Toshiaki Nakazawa, Martin Popel, Maja Popović, and Mariya Shmatova.  Findings of the 2023 Conference on Machine Translation (WMT23): LLMs Are Here but Not Quite There Yet.  In Philipp Koehn, Barry Haddow, Tom Kocmi, and Christof Monz (eds.), _Proceedings of the Eighth Conference on Machine Translation_ , pp. 1–42, Singapore, December 2023. Association for Computational Linguistics.  doi: 10.18653/v1/2023.wmt-1.1.
  * Kudo & Richardson (2018) Taku Kudo and John Richardson.  SentencePiece: A simple and language independent subword tokenizer and detokenizer for neural text processing.  In Eduardo Blanco and Wei Lu (eds.), _Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing: System Demonstrations_ , pp. 66–71, Brussels, Belgium, November 2018. Association for Computational Linguistics.  doi: 10.18653/v1/D18-2012.  URL https://aclanthology.org/D18-2012.
  * Liang et al. (2023) Davis Liang, Hila Gonen, Yuning Mao, Rui Hou, Naman Goyal, Marjan Ghazvininejad, Luke Zettlemoyer, and Madian Khabsa.  XLM-V: Overcoming the vocabulary bottleneck in multilingual masked language models.  In Houda Bouamor, Juan Pino, and Kalika Bali (eds.), _Proceedings of the 2023 Conference on Empirical Methods in Natural Language Processing_ , pp. 13142–13152, Singapore, December 2023. Association for Computational Linguistics.  doi: 10.18653/v1/2023.emnlp-main.813.  URL https://aclanthology.org/2023.emnlp-main.813.
  * Limisiewicz et al. (2023) Tomasz Limisiewicz, Jiří Balhar, and David Mareček.  Tokenization Impacts Multilingual Language Modeling: Assessing Vocabulary Allocation and Overlap Across Languages.  In Anna Rogers, Jordan Boyd-Graber, and Naoaki Okazaki (eds.), _Findings of the Association for Computational Linguistics: ACL 2023_ , pp. 5661–5681, Toronto, Canada, July 2023. Association for Computational Linguistics.  doi: 10.18653/v1/2023.findings-acl.350.
  * Lin et al. (2022) Xi Victoria Lin, Todor Mihaylov, Mikel Artetxe, Tianlu Wang, Shuohui Chen, Daniel Simig, Myle Ott, Naman Goyal, Shruti Bhosale, Jingfei Du, Ramakanth Pasunuru, Sam Shleifer, Punit Singh Koura, Vishrav Chaudhary, Brian O’Horo, Jeff Wang, Luke Zettlemoyer, Zornitsa Kozareva, Mona Diab, Veselin Stoyanov, and Xian Li.  Few-shot Learning with Multilingual Language Models, October 2022.
  * Liu & Niehues (2022) Danni Liu and Jan Niehues.  Learning an artificial language for knowledge-sharing in multilingual translation, 2022.
  * Papineni et al. (2002) Kishore Papineni, Salim Roukos, Todd Ward, and Wei-Jing Zhu.  Bleu: a method for automatic evaluation of machine translation.  In _ACL_ , 2002.
  * Petrov et al. (2023) Aleksandar Petrov, Emanuele La Malfa, Philip Torr, and Adel Bibi.  Language Model Tokenizers Introduce Unfairness Between Languages.  In _Thirty-Seventh Conference on Neural Information Processing Systems_ , November 2023.
  * Rei et al. (2022a) Ricardo Rei, José G. C. de Souza, Duarte Alves, Chrysoula Zerva, Ana C Farinha, Taisiya Glushkova, Alon Lavie, Luisa Coheur, and André F. T. Martins.  COMET-22: Unbabel-IST 2022 submission for the metrics shared task.  In Philipp Koehn, Loïc Barrault, Ondřej Bojar, Fethi Bougares, Rajen Chatterjee, Marta R. Costa-jussà, Christian Federmann, Mark Fishel, Alexander Fraser, Markus Freitag, Yvette Graham, Roman Grundkiewicz, Paco Guzman, Barry Haddow, Matthias Huck, Antonio Jimeno Yepes, Tom Kocmi, André Martins, Makoto Morishita, Christof Monz, Masaaki Nagata, Toshiaki Nakazawa, Matteo Negri, Aurélie Névéol, Mariana Neves, Martin Popel, Marco Turchi, and Marcos Zampieri (eds.), _Proceedings of the Seventh Conference on Machine Translation (WMT)_ , pp. 578–585, Abu Dhabi, United Arab Emirates (Hybrid), December 2022a. Association for Computational Linguistics.  URL https://aclanthology.org/2022.wmt-1.52.
  * Rei et al. (2022b) Ricardo Rei, José G. C. de Souza, Duarte Alves, Chrysoula Zerva, Ana C Farinha, Taisiya Glushkova, Alon Lavie, Luisa Coheur, and André F. T. Martins.  COMET-22: Unbabel-IST 2022 Submission for the Metrics Shared Task.  In Philipp Koehn, Loïc Barrault, Ondřej Bojar, Fethi Bougares, Rajen Chatterjee, Marta R. Costa-jussà, Christian Federmann, Mark Fishel, Alexander Fraser, Markus Freitag, Yvette Graham, Roman Grundkiewicz, Paco Guzman, Barry Haddow, Matthias Huck, Antonio Jimeno Yepes, Tom Kocmi, André Martins, Makoto Morishita, Christof Monz, Masaaki Nagata, Toshiaki Nakazawa, Matteo Negri, Aurélie Névéol, Mariana Neves, Martin Popel, Marco Turchi, and Marcos Zampieri (eds.), _Proceedings of the Seventh Conference on Machine Translation (WMT)_ , pp. 578–585, Abu Dhabi, United Arab Emirates (Hybrid), December 2022b. Association for Computational Linguistics.
  * Team et al. (2022) NLLB Team, Marta R. Costa-jussà, James Cross, Onur Çelebi, Maha Elbayad, Kenneth Heafield, Kevin Heffernan, Elahe Kalbassi, Janice Lam, Daniel Licht, Jean Maillard, Anna Sun, Skyler Wang, Guillaume Wenzek, Al Youngblood, Bapi Akula, Loic Barrault, Gabriel Mejia Gonzalez, Prangthip Hansanti, John Hoffman, Semarley Jarrett, Kaushik Ram Sadagopan, Dirk Rowe, Shannon Spruit, Chau Tran, Pierre Andrews, Necip Fazil Ayan, Shruti Bhosale, Sergey Edunov, Angela Fan, Cynthia Gao, Vedanuj Goswami, Francisco Guzmán, Philipp Koehn, Alexandre Mourachko, Christophe Ropers, Safiyyah Saleem, Holger Schwenk, and Jeff Wang.  No language left behind: Scaling human-centered machine translation, 2022.
  * Touvron et al. (2023) Hugo Touvron, Louis Martin, Kevin Stone, Peter Albert, Amjad Almahairi, Yasmine Babaei, Nikolay Bashlykov, Soumya Batra, Prajjwal Bhargava, Shruti Bhosale, Dan Bikel, Lukas Blecher, Cristian Canton Ferrer, Moya Chen, Guillem Cucurull, David Esiobu, Jude Fernandes, Jeremy Fu, Wenyin Fu, Brian Fuller, Cynthia Gao, Vedanuj Goswami, Naman Goyal, Anthony Hartshorn, Saghar Hosseini, Rui Hou, Hakan Inan, Marcin Kardas, Viktor Kerkez, Madian Khabsa, Isabel Kloumann, Artem Korenev, Punit Singh Koura, Marie-Anne Lachaux, Thibaut Lavril, Jenya Lee, Diana Liskovich, Yinghai Lu, Yuning Mao, Xavier Martinet, Todor Mihaylov, Pushkar Mishra, Igor Molybog, Yixin Nie, Andrew Poulton, Jeremy Reizenstein, Rashi Rungta, Kalyan Saladi, Alan Schelten, Ruan Silva, Eric Michael Smith, Ranjan Subramanian, Xiaoqing Ellen Tan, Binh Tang, Ross Taylor, Adina Williams, Jian Xiang Kuan, Puxin Xu, Zheng Yan, Iliyan Zarov, Yuchen Zhang, Angela Fan, Melanie Kambadur, Sharan Narang, Aurelien Rodriguez, Robert Stojnic, Sergey Edunov, and Thomas Scialom.  Llama 2: Open foundation and fine-tuned chat models, 2023.
  * Vanmassenhove et al. (2021) Eva Vanmassenhove, Dimitar Shterionov, and Matthew Gwilliam.  Machine Translationese: Effects of Algorithmic Bias on Linguistic Complexity in Machine Translation.  In Paola Merlo, Jorg Tiedemann, and Reut Tsarfaty (eds.), _Proceedings of the 16th Conference of the European Chapter of the Association for Computational Linguistics: Main Volume_ , pp. 2203–2213, Online, April 2021. Association for Computational Linguistics.  doi: 10.18653/v1/2021.eacl-main.188.
  * Vilar et al. (2023) David Vilar, Markus Freitag, Colin Cherry, Jiaming Luo, Viresh Ratnakar, and George Foster.  Prompting PaLM for Translation: Assessing Strategies and Performance, June 2023.
  * Xu et al. (2024a) Haoran Xu, Young Jin Kim, Amr Sharaf, and Hany Hassan Awadalla.  A paradigm shift in machine translation: Boosting translation performance of large language models, 2024a.
  * Xu et al. (2024b) Haoran Xu, Amr Sharaf, Yunmo Chen, Weiting Tan, Lingfeng Shen, Benjamin Van Durme, Kenton Murray, and Young Jin Kim.  Contrastive preference optimization: Pushing the boundaries of llm performance in machine translation, 2024b.
  * Xue et al. (2021) Linting Xue, Noah Constant, Adam Roberts, Mihir Kale, Rami Al-Rfou, Aditya Siddhant, Aditya Barua, and Colin Raffel.  mT5: A massively multilingual pre-trained text-to-text transformer.  In Kristina Toutanova, Anna Rumshisky, Luke Zettlemoyer, Dilek Hakkani-Tur, Iz Beltagy, Steven Bethard, Ryan Cotterell, Tanmoy Chakraborty, and Yichao Zhou (eds.), _Proceedings of the 2021 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , pp. 483–498, Online, June 2021. Association for Computational Linguistics.  doi: 10.18653/v1/2021.naacl-main.41.  URL https://aclanthology.org/2021.naacl-main.41.
  * Yang et al. (2023) Wen Yang, Chong Li, Jiajun Zhang, and Chengqing Zong.  BigTranslate: Augmenting Large Language Models with Multilingual Translation Capability over 100 Languages, November 2023.
  * Yuan et al. (2023) Fei Yuan, Yinquan Lu, WenHao Zhu, Lingpeng Kong, Lei Li, Yu Qiao, and Jingjing Xu.  Lego-mt: Learning detachable models for massively multilingual machine translation, 2023.
  * Zhang et al. (2023a) Biao Zhang, Barry Haddow, and Alexandra Birch.  Prompting Large Language Model for Machine Translation: A Case Study, January 2023a.
  * Zhang et al. (2023b) Shaolei Zhang, Qingkai Fang, Zhuocheng Zhang, Zhengrui Ma, Yan Zhou, Langlin Huang, Mengyu Bu, Shangtong Gui, Yunji Chen, Xilin Chen, and Yang Feng.  BayLing: Bridging Cross-lingual Alignment and Instruction Following through Interactive Translation for Large Language Models, June 2023b.
  * Zhu et al. (2023) Wenhao Zhu, Hongyi Liu, Qingxiu Dong, Jingjing Xu, Shujian Huang, Lingpeng Kong, Jiajun Chen, and Lei Li.  Multilingual Machine Translation with Large Language Models: Empirical Results and Analysis, May 2023.

##  Appendix A Appendix: Other Metrics for Translation Quality

Figure 6: Distribution of ChrF scores for each translation approach over
language pair supervision type.

We plot the distribution of ChrF scores in Figure 6. We generally see a
similar trend where fine-tuning shows improvement across all supervision
types. The range in scores is wider compared to COMET. This indicates that
translation outputs vary to a degree from the reference.

##  Appendix B Appendix: Textual diversity

Figure 7:

Figure 8:

##  Appendix C Appendix: Translation Quality per Model and Source/Target
Language

We illustrate the spread of translation quality scores for all models and
source or target languages in the figures 9 to 20 below, to complement Figures
2 in the main paper which show this spread per target languages, for the Nllb
and TowerInstruct-13B models only.

Figure 9: TO TARGET; S: similar, E: English, D+: best dissimilar, D-: worst
dissimilar  Figure 10: TO TARGET; S: similar, E: English, D+: best dissimilar,
D-: worst dissimilar  Figure 11: TO TARGET; S: similar, E: English, D+: best
dissimilar, D-: worst dissimilar  Figure 12: TO TARGET; S: similar, E:
English, D+: best dissimilar, D-: worst dissimilar  Figure 13: TO TARGET; S:
similar, E: English, D+: best dissimilar, D-: worst dissimilar  Figure 14:
FROM SOURCE; S: similar, E: English, D+: best dissimilar, D-: worst dissimilar
Figure 15: FROM SOURCE; S: similar, E: English, D+: best dissimilar, D-: worst
dissimilar  Figure 16: FROM SOURCE; S: similar, E: English, D+: best
dissimilar, D-: worst dissimilar  Figure 17: FROM SOURCE; S: similar, E:
English, D+: best dissimilar, D-: worst dissimilar  Figure 18: FROM SOURCE; S:
similar, E: English, D+: best dissimilar, D-: worst dissimilar  Figure 19:
FROM SOURCE; S: similar, E: English, D+: best dissimilar, D-: worst dissimilar
Figure 20: FROM SOURCE; S: similar, E: English, D+: best dissimilar, D-: worst
dissimilar

##  Appendix D Appendix: Impact of Target Tokenization

Figure 21: Average translation quality into target and average number of
generated tokens in target.