C Yuan

F Laurent

R O Fox

An extended quadrature method of moments for population balance equations

Keywords: Population balance equation, extended quadrature method of moments (EQMOM), aggregation, breakage, condensation, evaporation

Population balance equations (PBE) for a number density function (NDF) arise in many applications of aerosol technology. Thus, there has been considerable interest in the development of numerical methods to find solutions to PBE, especially in the context of spatially inhomogeneous systems where moment realizability becomes a significant issue. Quadrature-based moment methods (QBMM) are an important class of methods for which the accuracy of the solution can be improved in a controlled manner by increasing the number of quadrature nodes. However, when a large number of nodes is required to achieve the desired accuracy, the moment-inversion problem can become ill-conditioned. Moreover, oftentimes pointwise values of the NDF are required, but are unavailable with existing QBMM. In this work, a new generation of QBMM is introduced that provides an explicit form for the NDF. This extended quadrature method of moments (EQMOM) approximates the NDF by a sum of non-negative weight functions, which allows unclosed source terms to be computed with great accuracy by increasing the number of quadrature nodes independent of the number of transported moments. Here, we use EQMOM to solve a spatially homogeneous PBE with aggregation, breakage, condensation, and evaporation terms, and compare the results with analytical solutions whenever possible. However, by employing realizable finite-volume methods, the extension of EQMOM to spatially inhomogeneous systems is straightforward.

Introduction

Populations of discrete particles in a carrier fluid can be described by a population balance equation (PBE) [START_REF] Ramkrishna | Population Balances: Theory and Applications to Particulate Systems in Engineering[END_REF]. Technologically important examples include the general dynamic equation [START_REF] Friedlander | Smoke, Dust, and Haze: Fundamentals of Aerosol Dynamics, 2nd Edition[END_REF] used in aerosol reactor [START_REF] Pratsinis | Simultaneous nucleation, condensation, and coagulation in aerosol reactors[END_REF] and atmospheric models (Seinfeld and Pandis, ⇤ Corresponding author Email addresses: canshengyuan@hotmail.com (C. Yuan), frederique.laurent@ecp.fr (F. Laurent), rofox@iastate.edu (and R. O. Fox) 2006), and the spray equation [START_REF] Sirignano | Fluid Dynamics and Transport of Droplets and Sprays, 2nd Edition[END_REF] used to describe the evolution of fuel droplets in spray combustion devices. In the simplest case where the mean particle velocity is the same as the fluid, the PBE describes the evolution of a number density function (NDF) for the number of particles with given set of internal coordinates (e.g., volume, chemical composition, etc.) whose values lie in a high-dimensional phase space. The NDF depends on time t, spatial location x, and the values of the internal coordinates ⇠. A typical PBE contains terms for spatial transport (e.g., advection and di↵usion), source terms for the formation of new particles from the surrounding fluid, loss terms due to evaporation, growth terms on individual particles, and aggregation and breakage [START_REF] Vigil | Destructive aggregation: Aggregation with collision-induced breakage[END_REF] terms involving multiple particles. The mathematical form of the PBE is therefore quite complex with both hyperbolic (e.g., growth) and integral (e.g., aggregation) terms in phase space, in addition to the di↵erential terms for advection and di↵usion in real space. The complexity of the PBE, and the high-dimensional phase space, make the direct solution of the PBE intractable for use in time-dependent transport codes seeking to model the spatial dependence of the NDF, such as those used in climate [START_REF] Seinfeld | Atmospheric Chemistry and Physics: From Air Pollution to Climate Change, 2nd Edition[END_REF] and aerosol reactor [START_REF] Mehta | Multiscale modeling of TiO 2 nanoparticle production in flame reactors: E↵ect of chemical mechanism[END_REF][START_REF] Pratsinis | Simultaneous nucleation, condensation, and coagulation in aerosol reactors[END_REF][START_REF] Sung | Large-eddy simulation based multiscale modeling of TiO 2 nanoparticle synthesis in turbulent flame reactors using detailed nucleation chemistry[END_REF] modeling.

In order to reduce the dimensionality of the PBE, a number of di↵erent strategies have been proposed.

In classical moment methods [START_REF] Barrett | A comparison of some approximate methods for solving the aerosol general dynamic equation[END_REF], the PBE is multiplied by test functions (e.g., integers powers of the internal coordinates) and integrated over phase space. Because the resulting moment transport equations are not closed in terms of a finite set of moments, this approach leads to a moment closure problem. Broadly speaking, there are two ways to achieve closure: (i) provide a functional dependence of the unknown moments using the transported moment set (e.g., an interpolative closure such as MOMIC [START_REF] Frenklach | Aerosol dynamics modeling using the method of moments[END_REF]), or (ii) reconstruct the NDF from the transported moments (i.e., the truncated Hausdor↵/Stieljes/Hamburger moment problem [START_REF] Hausdor↵ | Momentenprobleme für ein endliches intervall[END_REF][START_REF] Shohat | The Problem of Moments[END_REF] from which the unclosed terms can be evaluated. In either case, a very important consideration for problems involving spatial transport is the realizability of the moment set [START_REF] Wright | Numerical advection of moments of the particle size distribution in Eulerian models[END_REF]. In other words, even if a suitable closure can be found for the moment transport equation, numerical advection and di↵usion schemes can lead to moment sets that do not correspond to a realizable NDF (i.e., the NDF must be non-negative on the support of the internal coordinates). An ad hoc approach to ensuring realizability can be formulated in terms of a moment correction algorithm using, for example, the Hankel-Hadamard determinants [START_REF] Shohat | The Problem of Moments[END_REF]. However, such approaches are di cult to apply to multivariate moment sets, and are generally unsatisfactory because they attempt to correct errors introduced during the numerical approximation of the moment transport equations that do not exist at the level of the PBE. In recent works [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF][START_REF] Vikas | Realizable high-order finite-volume schemes for quadrature-based moment methods[END_REF], promising alternative approaches that guarantee realizable moments have been introduced and make use of the reconstructed NDF to ensure realizability for high-order finite-volume reconstruction schemes on unstructured grids.

In light of realizability, in this work we restrict our attention to moment closures that reconstruct the NDF from the moments. The simplest moment-inversion algorithms use a parameterized NDF where the unknown parameters are found from a few lower-order moments (e.g., a log-normal distribution) [START_REF] Hutton | Particle size distribution reconstruction: The moment surface method[END_REF]. However, such methods do not provide enough flexibility to describe a complex multi-modal NDF.

A more flexible formulation, introduce by Grad [START_REF] Grad | On the kinetic theory of rarefied gases[END_REF] in kinetic theory, consists of expanding the NDF in a family of known orthogonal polynomials in terms of the internal coordinates. For a univariate PBE with size or volume as the internal coordinate, the Legendre polynomials can be employed (e.g., FCMOM [START_REF] Strumendo | Solution of PBE by MOM in finite size domains[END_REF]), or any other family defined on a subset of the real line (e.g., DuQMoGeM [START_REF] Lage | On the representation of QMOM as a weighted-residual method -The dual-quadrature method of generalized moments[END_REF]). The principal shortcoming of the Grad approach is that because the number of polynomials that can be determined from a finite set of moments is finite, the reconstructed NDF will almost always be negative for some values of the internal coordinate. Thus, while small negative values can be tolerated for approximating the integral terms in the PBE, they can lead to instabilities in the spatial transport terms and, almost inevitably [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF][START_REF] Vikas | Realizable high-order finite-volume schemes for quadrature-based moment methods[END_REF][START_REF] Wright | Numerical advection of moments of the particle size distribution in Eulerian models[END_REF], will yield unrealizable moments. For this reason, only moment-inversion methods that guarantee a nonnegative NDF are acceptable for approximating solutions to the PBE in the context of spatially inhomogeneous systems.

In aerosol science, a now widely used closure for the moment transport equation is the quadrature method of moments (QMOM) [START_REF] Mcgraw | Description of aerosol dynamics by the quadrature method of moments[END_REF]. Strictly speaking, QMOM does not reconstruct the NDF. However, given the one-to-one relationship between the QMOM moment closure and point distribution function [START_REF] Dette | The Theory of Canonical Moments with Applications in Statistics, Probability, and Analysis[END_REF], it is natural to interpret QMOM with n nodes to be an n-point distribution function (i.e., the sum of n weighted Dirac delta functions). Indeed, from the theory of moments [START_REF] Dette | The Theory of Canonical Moments with Applications in Statistics, Probability, and Analysis[END_REF][START_REF] Shohat | The Problem of Moments[END_REF], the QMOM closure for the moment m 2n is the smallest value possible given the moment set (m 0 , . . . , m 2n 1 ). In other words, the QMOM reconstructed NDF has moments of order higher than 2n that lie on the boundary of moment space, and all other possible NDF with the same moments (m 0 , . . . , m 2n 1 ) lie in the interior of moment space. Thus, QMOM uniquely chooses one of an infinite number of NDF from the transported moment set.

The representation of the NDF as an n-point distribution is exploited in the direct quadrature method of moments (DQMOM) [START_REF] Marchisio | Solution of population balance equations using the direct quadrature method of moments[END_REF] to formulate transport equations for the weights and abscissas. Theoretically, an exact and regular solution to the DQMOM equations should be identical to the solution of the moment transport equations. However, since the moments are the natural "conserved" quantities needed to represent the PBE, and the relationship between the moments and the abscissas is nonlinear, conservative finite-volume schemes applied to the weights and abscissas cannot guarantee conservation of the moments. A even more serious shortcoming associated with the n-point NDF arises in the treatment of evaporation terms [START_REF] Fox | Numerical simulation of spray coalescence in an Eulerian framework: direct quadrature method of moments and multi-fluid method[END_REF]. The key problem here is that evaporation leads to a zero-order moment equation containing a term corresponding to the loss of particles of zero size. In order to evaluate this term, the value of the NDF at ⇠ = 0 is required, but not available in QMOM. A partial (but impractical) solution would be to increase n to a large enough value (e.g., n > 100) so that phase space is adequately "discretized". However, the moment-inversion algorithm used in QMOM is not accurate for n larger than about ten [START_REF] Gautschi | Orthogonal Polynomials: Computation and Approximation[END_REF][START_REF] Mcgraw | Description of aerosol dynamics by the quadrature method of moments[END_REF][START_REF] Wheeler | Modified moments and Gaussian quadratures[END_REF]. Thus, the only feasible alternative to capture accurately the e↵ect of evaporation is to reconstruct a continuous NDF that can be evaluated at ⇠ = 0 [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF].

As mentioned above, there exists an infinite number of continuous NDF with moments in the interior of moment space. Thus, an additional criteria is needed to choose one of them to represent the NDF. For example, the entropy maximization (EM) method [START_REF] Mead | Maximum entropy in the problem of moments[END_REF][START_REF] Tagliani | Hausdor↵ moment problem and maximum entropy: A unified approach[END_REF] chooses the NDF that minimizes a functional subject to moment constraints. In practice, the numerical implementation of the EM method requires the solution of a constrained multi-variate minimization problem whose dimension depends on the number of transported moments. An obvious advantage of EM over the Grad method is that the NDF is guaranteed to be non-negative. However, the extension of EM to the boundaries of moment space is ill-conditioned [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF], and the treatment of multiple internal coordinates increases the numerical di culties considerably. Another promising method for reconstructing a realizable NDF is the kernel density element method (KDEM) [START_REF] Athanassoulis | The truncated Hausdor↵ moment problem solved by using kernel density functions[END_REF] and related methods [START_REF] Diemer | A moment methodology for coagulation and breakage problems: Part 2. Moment models and distribution reconstruction[END_REF], which uses a weighted sum of known non-negative kernel density functions to represent the NDF. As in the EM method, KDEM fixes the unknown parameters by solving a constrained minimization problem wherein only a few of the lowest-order moments are exactly reproduced by the reconstructed NDF.

Thus, in comparison to QMOM and EM where all 2n moments are exactly reproduced, KDEM introduces a quadrature error into the moment closures.

In this work we develop an alternative moment-inversion algorithm that combines the most desirable properties of QMOM and KDEM, while eliminating their weaknesses. The basic idea [START_REF] Chalons | A multi-Gaussian quadrature method of moments for gas-particle flows in a LES framework[END_REF] behind the extended quadrature method of moments (EQMOM) is to choose a kernel density function, depending on a single parameter , for which the QMOM moment-inversion algorithm can be applied directly to find the weights and abscissas. The parameter is then determined by forcing one additional transported moment (i.e., m 2n ) to agree with the reconstructed NDF. Thus, in place of a multi-variate minimization algorithm, is determined using a one-dimensional root-finding method with relatively modest computational cost. Moreover, in the limit ! 0, EQMOM reduces to QMOM in a well-conditioned manner, making it possible to easily reconstruct the NDF all the way to the boundary of moment space.

(See Appendix A for a comparison of EQMOM with EM methods.) As we shall show in Sec. 4, another important advantage of EQMOM is the ability to construct a second Gaussian quadrature [START_REF] Gautschi | Orthogonal Polynomials: Computation and Approximation[END_REF] with respect to the kernel density function that is independent of the number of transported moments.

Thus, as with DuQMoGeM [START_REF] Lage | On the representation of QMOM as a weighted-residual method -The dual-quadrature method of generalized moments[END_REF], the accuracy of the moment closures used for integral terms can be greatly increased as very low computational cost.

The focus of this paper is on the formulation and validation of EQMOM in the context of a univariate NDF for an internal coordinate (such as particle volume or size) defined on either a semi-finite or finite subset of the real line. In Sec. 2 we briefly describe the PBE and the moment transport equations for a spatially homogeneous system. In Sec. 3 we introduce the mathematical formulation of EQMOM and describe the numerically robust moment-inversion algorithm (referred to hereinafter as the first quadrature) used in later sections. Then, in Sec. 4, we describe how moment closures are constructed using the second Gaussian quadrature. Next, in Sec. 5, EQMOM is applied to 13 test cases, many of which have analytical solutions, in order to investigate the accuracy of the moment closures and to illustrate how the EQMOM predictions depend on the parameters used in the first and second quadratures. Finally, in Sec. 6 we summarize the principal conclusions and briefly mention how EQMOM can be generalized to multiple internal coordinates and spatially inhomogeneous systems.

Moment Methods for Population Balance Equations

In this work we consider a PBE with growth (i.e., evaporation and condensation), aggregation, breakage and source terms. For clarity, we consider only spatially homogeneous systems. However, the methods developed here can be applied to a spatially inhomogeneous PBE using realizable finite-volume schemes [START_REF] Vikas | Realizable high-order finite-volume schemes for quadrature-based moment methods[END_REF]. In order to guarantee realizability of the moments sets, such schemes make use of the underlying connection between the PBE and the moment equations. The same is true for the moment method developed in this work.

Population balance equation

Consider the following PBE for the spatially homogeneous, univariate NDF f (t, ⇠):

@f (t, ⇠) @t + @ @⇠ [g(t, ⇠)f (t, ⇠)] = B agg (t, ⇠) D agg (t, ⇠) + B br (t, ⇠) D br (t, ⇠) + Q(t, ⇠)f (t, ⇠) (1)
where ⇠ 2 [0, ⇠ max ] (with ⇠ max possibly infinite) depends on what physical processes are investigated (e.g, size, surface area, or volume of particles) and Q(t, ⇠) is a source term used in this work to obtain a known form for f (t, ⇠). The aggregation terms are given by

B agg (t, ⇠) = 1 2 Z ⇠ 0 a(t, ⇠ ⇠ 0 , ⇠ 0 )f (t, ⇠ ⇠ 0 )f (t, ⇠ 0 ) d⇠ 0 (2)
and

D agg (t, ⇠) = Z ⇠max 0 a(t, ⇠, ⇠ 0 )f (t, ⇠)f (t, ⇠ 0 ) d⇠ 0 (3)
where a(t, ⇠, ⇠ 0 ) is the aggregation kernel. The breakage terms can be expressed as

B br (t, ⇠) = Z ⇠max ⇠ #(t, ⇠ 0 )b(t, ⇠ 0 )P (t, ⇠|⇠ 0 )f (t, ⇠ 0 ) d⇠ 0 (4)
and

D br (t, ⇠) = b(t, ⇠)f (t, ⇠) (5)
where #(t, ⇠ 0 ) is the average number of particles formed by breakage of a particle of type ⇠ 0 at time t, P (t, ⇠|⇠ 0 ) is probability density function for particle type ⇠ 0 at time t to produce a daughter particle with type ⇠, and b(t, ⇠) is the breakage kernel. In the examples in Sec. 5, the source term is used to force the exact solution of the PBE to have a chosen form for the NDF (i.e., a so-called manufactured solution). The accuracy of the proposed moment method can then be evaluated against the moments of the known NDF.

Moment equations

The integer moments of f (t, ⇠) are defined by

m k (t) ⌘ Z ⇠max 0 ⇠ k f (t, ⇠) d⇠ (6)
with ⇠ max finite or infinite. The integro-di↵erential equation for the moments can be found starting from Eq. ( 1):

dm k dt = g(t, ⇠)f (t, ⇠)⇠ k ⇠max 0 + Z ⇠max 0 k⇠ k 1 g(t, ⇠)f (t, ⇠) d⇠ + Z ⇠max 0 ⇠ k [B agg (t, ⇠) D agg (t, ⇠)] d⇠ + Z ⇠max 0 ⇠ k [B br (t, ⇠) D br (t, ⇠)] d⇠ + Z ⇠max 0 ⇠ k Q(t, ⇠)f (t, ⇠) d⇠. (7)
In general, this system of moment equations is not closed in terms of a finite set of moments due to the nonlinear dependencies of the growth, aggregation, and breakage terms on ⇠, which introduce higher-order (and non-integer) moments. In order to attain closure, a novel quadrature-based moment method (QBMM) is introduced in Sec. 3. However, we should note that the mathematical forms of the unclosed terms in Eq. ( 7) can make closure very challenging.

For example, the growth terms in Eq. ( 1) are hyperbolic in nature (i.e., solutions evolve along characteristics in ⇠-phase space), while the aggregation and breakage terms lend themselves to QBMM. In particular, when the growth term g(t, ⇠) is negative (i.e., evaporation [START_REF] Fox | Numerical simulation of spray coalescence in an Eulerian framework: direct quadrature method of moments and multi-fluid method[END_REF][START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF] and fines dissolution [START_REF] Grosch | Generalization and numerical investigation of QMOM[END_REF]), the zero-order moment m 0 requires a closure for the boundary flux term g(t, 0)f (t, 0), i.e., a pointwise value of the NDF. In standard QMOM [START_REF] Mcgraw | Description of aerosol dynamics by the quadrature method of moments[END_REF], pointwise values of the NDF are not available, and thus an important property of our proposed extension of QMOM is the ability to accurately evaluate the functional form of f (t, ⇠) at a discrete point ⇠. A related class of problems for which QMOM performs poorly is a PBE with selective removal of particles in a small interval of size space [START_REF] Grosch | Generalization and numerical investigation of QMOM[END_REF]) (e.g., fines removal or separation of large particles). For such problems, one needs a good estimate of the integral of f (t, ⇠) over a finite interval ⇠ 1  ⇠  ⇠ 2 , which -for reasonable accuracy -requires knowledge of f (t, ⇠) at several points in the interval.

Extended Quadrature Method of Moments

EQMOM is conceptually equivalent to the bi-Gaussian quadrature introduced elsewhere [START_REF] Chalons | A multi-Gaussian quadrature method of moments for gas-particle flows in a LES framework[END_REF], and shares many similarities (but important di↵erences) with the KDEM (Athanassoulis and Gavril-iadis, 2002).

A non-negative approximation for the NDF

The starting point of EQMOM is to represent the distribution function by a weighted sum of non-negative functions:

p n (⇠) = n X ↵=1 w ↵ (⇠; ⇠ ↵ ) (8)
where w ↵ are non-negative weights, ⇠ ↵ the corresponding abscissas, and (⇠; ⇠ ↵ ) is a kernel density function with a finite (or infinite) support determined by the parameter . Recall that in QMOM, the NDF is represented as

p n (⇠) = n X ↵=1 w ↵ (⇠ ⇠ ↵ ). (9)
Thus, for consistency, we require that (⇠; ⇠ ↵ ) be chosen to satisfy the following condition:

lim !0 (⇠; ⇠ ↵ ) = (⇠ ⇠ ↵ ), (10)
which is not di cult to obtain for positive kernel density functions [START_REF] Athanassoulis | The truncated Hausdor↵ moment problem solved by using kernel density functions[END_REF].

Notice that the same is shared for all ↵ 2 1, . . . , n in Eq. ( 8).

The reader interested in the mathematical properties of Eq. ( 8) as an approximation of f (t, ⇠) can consult the literature [START_REF] Athanassoulis | The truncated Hausdor↵ moment problem solved by using kernel density functions[END_REF]. In particular, it can be shown that Eq. ( 8) can approximate any continuous density function as closely as desired with some finite n.

The moment-inversion problem

For the EQMOM approximation, 2n + 1 moments are employed to compute the 2n + 1 unknowns: w ↵ , ⇠ ↵ for ↵ = 1, . . . , n, and . In fact, the key improvement of EQMOM as compared to KDEM is the ability to compute these unknowns exactly from the moment set

M 2n = {m 0 , m 1 , . . . , m 2n }
using a robust moment-inversion algorithm based on the adaptive Wheeler algorithm [START_REF] Yuan | Conditional quadrature method of moments for kinetic equations[END_REF].

The reader should note that in the limit = 0, the moment set

M 2n 1 = {m 0 , m 1 , . . . , m 2n 1 }
can be used to construct the QMOM representation of the NDF in Eq. ( 9). Thus, in essence, the moment m 2n is used to determine the value of in Eq. ( 8).

In the literature [START_REF] Athanassoulis | The truncated Hausdor↵ moment problem solved by using kernel density functions[END_REF], a detailed discussion is available of the numerical issues faced when trying to invert Eq. ( 8) directly to yield a prescribed set of moments (i.e., the Haus-dor↵ truncated moment problem [START_REF] Hausdor↵ | Momentenprobleme für ein endliches intervall[END_REF]). In the present work, the adaptive Wheeler algorithm [START_REF] Wheeler | Modified moments and Gaussian quadratures[END_REF][START_REF] Yuan | Conditional quadrature method of moments for kinetic equations[END_REF] (see also [START_REF] Dette | The Theory of Canonical Moments with Applications in Statistics, Probability, and Analysis[END_REF][START_REF] Mcgraw | Description of aerosol dynamics by the quadrature method of moments[END_REF]) is used for moment inversion in order to minimize the numerical di culties as much as possible (i.e., the recursion coe cients are computed from the moments directly instead of solving a minimization problem as done elsewhere [START_REF] Athanassoulis | The truncated Hausdor↵ moment problem solved by using kernel density functions[END_REF]). Even so, the maximum value of n that can be used with reasonable numerical accuracy is on the order of five to ten, and hence our main focus is on accurately predicting the moments up to order 2n + 1, as opposed to predicting the (unknown) exact NDF.

For a Gaussian kernel density function defined for ⇠ 2 ( 1, +1), the moment-inversion method used in EQMOM is described elsewhere [START_REF] Chalons | A multi-Gaussian quadrature method of moments for gas-particle flows in a LES framework[END_REF] for n = 2. For this kernel density function, the abscissas ⇠ ↵ correspond to the mean and to the standard deviation of the Gaussian. In this work, we consider kernel density functions that are appropriate for cases where ⇠ 2 [0, 1) and ⇠ 2 [0, 1]. In these cases, the abscissas and no longer correspond to the mean and standard deviation of the kernel density function and, in fact, we shall see that their definitions are critical for defining a well-posed moment-inversion algorithm. Note that by a linear change of variables, any semi-infinite interval [a, 1) or finite interval [a, b] can be transformed to [0, 1) and [0, 1], respectively. Thus, the EQMOM algorithms developed in this work can be applied to any semi-infinite or finite interval used to define the range of ⇠.

Gamma EQMOM

For ⇠ 2 [0, 1), a gamma distribution is a good choice for (⇠; ⇠ ↵ ):

(⇠; ⇠ ↵ ) = ⇠ ↵ 1 e ⇠/ ( ↵ ) ↵ , (11)
where ↵ = ⇠ ↵ / . Then f (t, ⇠) can be approximated by

p n (t, ⇠) = n X ↵=1 w ↵ ⇠ ↵ 1 e ⇠/ ( ↵ ) ↵ . ( 12
)
It remains to show how the moments are related to the parameters ↵ and .

The first step is to compute the integer moments of (⇠; ⇠ ↵ ) with respect to ⇠ analytically. For the gamma EQMOM, these integer moments are

m (↵) k = ( ↵ + k) ( ↵ ) ✓ ⇠ ↵ ↵ ◆ k , (13)
and thus

m k = n X ↵=1 w ↵ ( ↵ + k) ( ↵ ) k = n X ↵=1 w ↵ G k (⇠ ↵ , ) (14)
where

G k (⇠ ↵ , ) = 8 > < > : 1 i f k = 0, Q k 1 i=0 (⇠ ↵ + i ) if k 1. ( 15
)
Noting that G k can be written as

G k (⇠ ↵ , ) = ⇠ k ↵ + P k 1 (⇠ ↵ , ) (16)
where P k 1 (x, y) is a homogeneous polynomial of order k 1 in x and y, we can rewrite the integer moments as

m k = m ⇤ k + n X ↵=1 w ↵ P k 1 (⇠ ↵ , ) (17)
where

m ⇤ k = n X ↵=1 w ↵ ⇠ k ↵ . ( 18
)
The right-hand side of Eq. ( 17) can therefore be expressed in terms of and the moments m ⇤ 0 , . . . , m ⇤ k . In other words, Eq. ( 17) forms a lower-triangular linear system that can be inverted using forward substitutions.

As an example, up to k = 4 (which corresponds to n = 2), we have

m 0 = m ⇤ 0 , m 1 = m ⇤ 1 , m 2 = m ⇤ 2 + m ⇤ 1 , m 3 = m ⇤ 3 + 3 m ⇤ 2 + 2 2 m ⇤ 1 , m 4 = m ⇤ 4 + 6 m ⇤ 3 + 11 2 m ⇤ 2 + 6 3 m ⇤ 1 . (19)
These equations can be recast in matrix form:

m = A( )m ⇤ . ( 20
)
This system of equations can be solved iteratively to find using the methods described in Sec. 3.5 below.

It is worth mentioning that only by choosing ↵ as in Eq. (11) will A( ) be a closed lower-triangular matrix. This is a key technical point because, for a given , it allows us to find the moment set (m ⇤ 0 , . . . , m ⇤ 2n 1 ) from (m 0 , . . . , m 2n 1 ), and thus to find the weights w ↵ and abscissas ⇠ ↵ using the adaptive Wheeler algorithm [START_REF] Yuan | Conditional quadrature method of moments for kinetic equations[END_REF]. In this manner, the moments set (m 0 , . . . , m 2n 1 ) is exactly recovered for any value of , and we are free to fix the value of to agree with moment m 2n . This is accomplished using the scalar function

J n ( ) = m 2n m ⇤ 2n n X ↵=1 w ↵ P 2n 1 (⇠ ↵ , ) (21)
by finding the smallest for which J n ( ) = 0.

Beta EQMOM

If ⇠ is in the bounded interval [0, 1], (⇠; ⇠ ↵ ) can be set to a beta distribution. The NDF for ⇠ 2 [0, 1] is then approximated by

p n (⇠) = n X ↵=1 w ↵ ⇠ ↵ 1 (1 ⇠) µ↵ 1 B( ↵ , µ ↵ ) (22)
where the two parameters are defined by ↵ = ⇠ ↵ / and µ ↵ = (1 ⇠ ↵ )/ . For the beta distribution, the integer moments of (⇠; ⇠ ↵ ) can be found from a recursion formula:

m (↵) k = ⇠ ↵ + (k 1) 1 + (k 1) m (↵) k 1 for k > 0, (23)
and m

(↵) 0

= 1. We can thus express the integer moments of the distribution function in Eq (22) as

m k = n X ↵=1 w ↵ G k (⇠ ↵ , ) (24)
where

G k (⇠ ↵ , ) = 8 > < > : 1 i f k = 0, Q k 1 i=0 ( ⇠↵+i 1+i ) if k 1, (25)
which is very similar to Eq. ( 15).

Because ⇠ ↵ appears only in the numerator in Eq. ( 25), the product will yield polynomials in ⇠ ↵ of order k. Thus, we can express the integer moments of the distribution function in Eq. ( 22) as

m k = k m ⇤ k + k 1 m ⇤ k 1 + • • • + 1 m ⇤ 1 ( 26
)
where the non-negative coe cients k depend only on .

For example, up to k = 4, we have

m 0 = m ⇤ 0 , m 1 = m ⇤ 1 , m 2 = 1 1 + (m ⇤ 2 + m ⇤ 1 ), m 3 = 1 (1 + 2 )(1 + ) (m ⇤ 3 + 3 m ⇤ 2 + 2 2 m ⇤ 1 ), m 4 = 1 (1 + 3 )(1 + 2 )(1 + ) (m ⇤ 4 + 6 m ⇤ 3 + 11 2 m ⇤ 2 + 6 3 m ⇤ 1 ), (27)
where m ⇤ k has the same definition as in gamma EQMOM. Again, the system of equations in Eq. ( 26) can be written as m = A( )m ⇤ where A( ) is a lower triangular matrix. For beta EQMOM, the scalar function is

J n ( ) = m 2n 2n m ⇤ 2n 2n 1 m ⇤ 2n 1 • • • 1 m ⇤ 1 , (28)
and we again find the smallest for which J n ( ) = 0. The system of moment equations is solved iteratively to find using the algorithm described next.

The EQMOM moment-inversion algorithm

To solve a system of equations such as Eq. ( 19) or Eq. ( 27), we propose the following algorithm. Given the 2n + 1 realizable moments m = M 2n , 1. guess , and compute the 2n moments m ⇤ k for k = 0, . . . , 2n 1 from m ⇤ = A( ) 1 m;

2. use the adaptive Wheeler algorithm with m ⇤ k for k = 0, . . . , 2n 1 to find n weights w ↵ and n abscissas

⇠ ↵ ;
3. compute m ⇤ 2n using w ↵ and ⇠ ↵ ; 4. compute J n ( ) from m ⇤ and ; 5. if J n ( ) 6 = 0, compute a new guess for and iterate until convergence.

In this algorithm, the adaptive Wheeler algorithm introduced in Yuan and Fox ( 2011) is applied so that n is chosen adaptively based on the moment error. Because of the highly nonlinear dependence of J n on , it is not convenient to find an analytical expression for the derivative J 0 n . Thus, we have successfully employed both a bounded secant method and the Ridder's method [START_REF] Press | Numerical Recipes in Fortran 77: The Art of Scientific Computing[END_REF] to update . For completeness, we should note that for a Gaussian kernel with n > 2, the same algorithm is applied with J n ( ) defined in terms of the Gaussian moments [START_REF] Chalons | A multi-Gaussian quadrature method of moments for gas-particle flows in a LES framework[END_REF].

From the theory of moments [START_REF] Dette | The Theory of Canonical Moments with Applications in Statistics, Probability, and Analysis[END_REF], we know that m 2n m ⇤ 2n , with equality implying that the NDF is exactly represented by the sum of at most n weighted delta functions. The scalar function has the properties J n (0) 0 and J 0 n (0) < 0. Thus, when J n (0) > 0, we wish to find the smallest value of for which (i) J n ( ) = 0 and (ii) the moment set (m ⇤ 0 , . . . , m ⇤ 2n 1 ) is realizable (i.e., all w ↵ > 0 and ⇠ ↵ 2 [0, 1) for gamma EQMOM or ⇠ ↵ 2 [0, 1] for beta EQMOM). In practice, we find that the smallest value of for which J( ) = 0 can sometimes yield an abscissa out of range. In such cases, we choose to be the largest value for which all abscissas are realizable, and redefine J n accordingly in the iteration algorithm.

(See Appendix B for details.) Another possible case is that J n ( ) = 0 but one (or more) w ↵ = 0. When this occurs, it implies that less than n kernel density functions are needed to exactly capture the first 2n moments (e.g., f (t, ⇠) is exactly a gamma/beta distribution so that n = 1 su ces to reproduce it). A final possible case is that the moment set (m ⇤ 0 , . . . , m ⇤ 2n 1 ) is unrealizable but J n ( ) > 0, which implies that we cannot reproduce all 2n + 1 moments for a given n. In such cases, we reduce n by one and repeat the search process using less moments.

In the literature [START_REF] Athanassoulis | The truncated Hausdor↵ moment problem solved by using kernel density functions[END_REF][START_REF] Diemer | A moment methodology for coagulation and breakage problems: Part 2. Moment models and distribution reconstruction[END_REF], the NDF has been approximated by a kernel density function whose coe cients are computed by solving a constrained, nonnegative, least-squares problem. In this manner, usually only one or two constraints are applied, which means that only one or two moments are exact and the rest are approximate. In comparison, with the EQMOM moment-inversion algorithm the moment set M 2n 1 is almost always exactly reproduced and moment m 2n is as accurate as possible. Furthermore, the iterative scheme described above is one-dimensional since only the value of is unknown. We are thus guaranteed to quickly converge to an optimal solution with relatively little computational cost, which is an important requirement for solving spatially inhomogeneous problems [START_REF] Vikas | Realizable high-order finite-volume schemes for quadrature-based moment methods[END_REF][START_REF] Yuan | Conditional quadrature method of moments for kinetic equations[END_REF].

Moment Closure with EQMOM

In this section, we assume that the parameters defining p n (⇠) have been successfully computed using the algorithm in Sec. 3.5, and address the question of how to compute the moment closures in Eq. ( 7). For this purpose, we follow the idea introduced elsewhere (Lage, 2011) of using a second Gaussian quadrature with respect to the kernel density function.

Choice of kernel density function

As done above, a convenient choice for univariate EQMOM is to define the kernel density function in terms of the weight function w(t) for a known family of orthogonal polynomials. For example, on the interval [ 1, 1] the associated Jacobi polynomials have the weight function

w(t) = (1 t) a (1 + t) b , (29)
where the parameters 1 < a, 1 < b defines a particular family of Jacobi polynomials [START_REF] Gautschi | Orthogonal Polynomials: Computation and Approximation[END_REF]).

As we can see, the beta distribution can be easily expressed in terms of this weight function by the change of variable ⇠ = (t + 1)/2. Likewise, for gamma EQMOM the corresponding weight function is

w(t) = t a e t , (30)
which is used with the associated Laguerre polynomials [START_REF] Gautschi | Orthogonal Polynomials: Computation and Approximation[END_REF], while for the Gaussian EQMOM the weight function corresponds to the Hermite polynomials. In principle, an EQMOM could be defined for any weight function with a corresponding set of orthogonal polynomials. As discussed in Appendix B, for fixed n the classical weight functions do not cover all of moment space (i.e., always up to m 2n 1 , but not always m 2n ). Nevertheless, as shown next, they allow us to accurately approximated integral terms at very low computational expense.

EQMOM approximation of integrals

The principal advantage of using the weight function for a known family of orthogonal polynomials is that the recursion coe cients are known in advance (i.e., we do not need to compute them from the moments

m (↵) k
of the kernel density function). Thus, the Jacobi matrix, which is used to find the N weights and N abscissas of the second quadrature [START_REF] Gautschi | Orthogonal Polynomials: Computation and Approximation[END_REF], can be computed with good accuracy for large values of N . In contrast, if the moments m (↵)

k were used with the Wheeler algorithm to find the second quadrature, the maximum value of N that could be employed with good accuracy would be on the order of ten.

As an example of the second quadrature, consider the beta kernel density function, which leads to the following quadrature formula for an arbitrary function g(⇠):

Z 1 0 g(⇠) (⇠; ⇠ ↵ ) d⇠ = 1 B( ↵ , µ ↵ ) Z 1 0 g(⇠)⇠ ↵ 1 (1 ⇠) µ↵ 1 d⇠ = ✓ 1 2 ◆ ↵ +µ↵ 1 1 B( ↵ , µ ↵ ) Z 1 1 g ✓ t + 1 2 ◆ (1 t) µ↵ 1 (1 + t) ↵ 1 dt ⇡ N↵ X =1 w ↵ g ✓ t ↵ + 1 2 ◆ . (31)
In the second line, we see the weight function in Eq. ( 29) with a ↵ = µ ↵ 1 and b ↵ = ↵ 1. In the last line, {w ↵ , t ↵ } are the N ↵ weights1 and N ↵ abscissas computed from the Jacobi polynomial recursion coe cients for a given set of parameters {a ↵ , b ↵ } using an algorithm described elsewhere [START_REF] Gautschi | Orthogonal Polynomials: Computation and Approximation[END_REF].

Using Eq. ( 31), we can see that the EQMOM closure defined by Eq. ( 8) approximates integrals with respect to the f (t, ⇠) as

Z ⇠1 ⇠0 g(⇠)p n (⇠) d⇠ = Z ⇠1 ⇠0 g(⇠) n X ↵=1 w ↵ (⇠; ⇠ ↵ ) d⇠ = n X ↵=1 N↵ X =1 ⇢ ↵ g(⇠ ↵ )I [⇠0,⇠1] (⇠ ↵ ) (32)
where Here, it is important to recognize that N ↵ in Eq. ( 32) can be chosen independently from the value of n. Therefore, we can choose N ↵ n to improve the accuracy of Eq. ( 32). As for any Gaussian quadrature [START_REF] Gautschi | Orthogonal Polynomials: Computation and Approximation[END_REF], the formula in Eq. ( 31) is exact if g(⇠) is a polynomial. For a 2n order polynomial g(⇠), since N ↵ quadrature nodes can yield an exact integration for polynomials of order 2N ↵ 1, the value of N ↵ should be chosen as 2N ↵ 1 2n so that Eq. ( 31) is exact. This is one important fact for choosing N ↵ in the test cases in Sec. 5 for cases where g(⇠) is a polynomial. Otherwise, N ↵ is chosen large enough to reduce the error in the second quadrature to a level that is smaller than the error in the first quadrature, which depends on n.

I [a,b] (x)

EQMOM approximation of hyperbolic terms

The approximation in Eq. ( 32) can be used to close the integral terms in Eq. ( 7). However, in order to close the hyperbolic terms (i.e., growth), we introduce an equivalent point representation of p n (⇠):

p n (⇠) = n X ↵=1 N↵ X =1 ⇢ ↵ (⇠ ⇠ ↵ ). ( 33
)
When this representation is substituted into the PBE, the hyperbolic terms on the left-hand side of Eq. ( 1) are described exactly by the ordinary di↵erential equations (ODE) for the characteristics:

d⇢ ↵ dt = 0 and d⇠ ↵ dt = g (⇠ ↵ ) . (34)
Thus, it is quite natural to employ a time-splitting method to solve the moment equations wherein the growth terms are updated using the method of characteristics (i.e., by solving Eq. ( 34) over each time step).

Because N ↵ can be chosen arbitrarily large, the number of characteristics can be made large enough to attain a specified numerical error in the evaluation of the growth terms. Moreover, because ⇢ ↵ is constant along a characteristic, the updated moments are guaranteed to be realizable. Finally, note that for evaporation with g(0) < 0, characteristics can leave the realizable domain (⇠ 2 [0, 1)). When this occurs during a time step, the corresponding weights are simply set to zero. In this manner, the flux of m 0 at ⇠ = 0 is captured in a numerically robust manner and can be computed with any desired accuracy by increasing N ↵ .

EQMOM approximation of the moment equations

In summary, the 2n + 1 moments M 2n are updated by solving Eq. ( 7) in the following manner:

1. Given M 2n at time t, use EQMOM to find the first quadrature: w ↵ , ⇠ ↵ and .

2. Given the first quadrature, find the second quadrature for each ↵: ⇢ ↵ and ⇠ ↵ .

3. For a given t, update M 2n due to growth by solving Eq. ( 34). For evaporation problems, the time step is fixed such that the two abscissas closest to the origin at time t (i.e., ⇠ 1 (t) and ⇠ 2 (t)) exactly straddle the origin at t + t (i.e., ⇠ 1 (t + t) = ⇠ 2 (t + t)).

4. Given the updated M 2n due to growth, recompute the first and second quadratures.

5. For the same t, update M 2n by solving

dm k dt = 1 2 n X ↵1=1 N↵ 1 X 1 =1 n X ↵2=1 N↵ 2 X 2 =1 ⇢ ↵1 1 ⇢ ↵2 2 h (⇠ ↵1 1 + ⇠ ↵2 2 ) k ⇠ k ↵1 1 ⇠ k ↵2 2 i a (t, ⇠ ↵1 1 , ⇠ ↵2 2 ) + n X ↵=1 N↵ X =1 ⇢ ↵ b(t, ⇠ ↵ ) ⇥ ⇠ k ↵ #(t, ⇠ ↵ ) ⇤ M k P (t, ⇠ ↵ ) + n X ↵=1 N↵ X =1 ⇢ ↵ Q(t, ⇠ ↵ ) (35)
where

M k P (t, ⇠) = R ⇠ 0 ⇠ 0k P (t, ⇠|⇠ 0
) d⇠ 0 can usually be found analytically. Note that the right-hand side of Eq. ( 35) was evaluated assuming that ⇠ is the particle volume.

For cases without growth, only Eq. ( 35) is solved using a realizable ODE solver [START_REF] Vikas | Realizable high-order finite-volume schemes for quadrature-based moment methods[END_REF].

In Sec. 5, we use a suite of test problems to assess the accuracy of the predicted moments as a function of the quadrature parameters. The reader should note that when functions such as a(t, ⇠, ⇠ 0 ) and b(t, ⇠) appearing in Eq. ( 35) are polynomials in ⇠, the value of N ↵ can be chosen large enough to make the summations over exact. We will return to this point in the specific examples given below.

Results and Discussion

In this section, we apply the EQMOM to a series of test cases taken from the literature [START_REF] Ernst | Coagulation processes with a phase transition[END_REF][START_REF] Lage | On the representation of QMOM as a weighted-residual method -The dual-quadrature method of generalized moments[END_REF][START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF][START_REF] Strumendo | Solution of PBE by MOM in finite size domains[END_REF]. These test cases were chosen either because they have analytical solutions or because they are known to be challenging to standard QMOM.

Evaporation

As mentioned earlier, problems with evaporation are particularly di cult to treat accurately with standard QMOM due to the lack of information concerning the pointwise flux [START_REF] Fox | Numerical simulation of spray coalescence in an Eulerian framework: direct quadrature method of moments and multi-fluid method[END_REF]. Thus, in our opinion, any successful extension of QMOM must be able to accurately reproduce the behavior of the zero-order moment (m 0 ) for arbitrary functional forms for g(t, ⇠). In pure evaporation problems, ⇠ is usually chosen to be the droplet surface area, and then, assuming that the evaporation rate is proportional to the droplet surface area (i.e., the d2 -evaporation law [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF]), g(t, ⇠) is a negative constant. More general evaporation expressions have a dependence on ⇠, but have g(t, 0) < 0, which implies that m 0 will decrease due to the loss of droplets at ⇠ = 0.

In this section, the following evaporation test cases, taken the literature [START_REF] Lage | On the representation of QMOM as a weighted-residual method -The dual-quadrature method of generalized moments[END_REF][START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF] 2 , are analyzed:

• Case 1: g(t, ⇠) = 1/2 with initial condition f (0, ⇠) = 60⇠ 2 (1 ⇠) 3 for ⇠ 2 [0, 1]. • Case 2: g(t, ⇠) = 1/2 with initial condition f (0, ⇠) = 1 for ⇠ 2 [0, 1]. • Case 3: g(t, ⇠) = (⇠ + 1/2) with initial condition f (0, ⇠) = (1 + 8⇠)(1 ⇠) 2 exp h 0.001 ⇣ 1 1 (1 ⇠) 2 ⌘i for ⇠ 2 [0, 1].
For the two cases with constant g, the analytical solution simply shifts the initial NDF to the left. For example, the analytical solution in Case 1 is

f (t, ⇠) = max[60(⇠ + t/2) 2 [1 (⇠ + t/2)] 3
, 0]. The evaporation test cases are simulated using beta EQMOM on the interval

m k,error (t) = |m k,EQMOM (t) m k,exact (t)| m k,exact (t)
for k = 0, . . . , 2n + 1 where m k,EQMOM is the kth-order moment predicted by EQMOM and m k,exact is the corresponding exact moment.

Case 1. For this case, the NDF predicted using beta EQMOM with n = 4 at di↵erent times is shown in Fig. 1. In Fig. 1(a), it can be seen that the EQMOM results found with n = 4 and N ↵ = 8 match quite well with the exact NDF. However, there are small deviations from the exact NDF close to the lower boundary at t = 0.4 and t = 0.8. The reason for these deviations is that N 1 and N 2 are too small. As discussed in Sec. 4.4, the time step is determined by the second quadrature nodes and N ↵ = 8 leads to a relatively large time step. Therefore, in Fig. 1(b), N 1 and N 2 are increased to 50, while leaving N 3 and N 4 set to 5. It can be observed that the resulting NDF falls almost exactly on the exact NDF. The reason why N 3 and N 4 can be set to 5 without incurring any errors is that as long as N ↵ n + 1, the second quadrature will not a↵ect the accuracy of the moment predictions.

[Figure 1 about here.]

The relative errors in the moments for Case 1 are shown in Fig. 2 for n = 4. From Fig. 2(a), the higherorder moment errors (k = 5, 6, 7, 8) start with smaller values as compared to m 0 , but gradually increase with time. However, at t = 1.2, they surpass m 0,error . This occurs mainly because for large times, the NDF tends to zero, leading to very small values of the higher-order moments m k,exact . Thus, the absolute moment errors remain very small as t ! 2 (i.e., complete evaporation). Similar to the NDF, the relative moment errors improve by increasing N 1 and N 2 . For this case, Lage (2011) uses an ODE solver with a fixed time step to solve the moment transport equations. Since the value of NDF at the origin (f (t, 0)), which is used to compute the mass flux out of domain, only depends on the first quadrature, increasing N ↵ cannot improve the accuracy of the flux. Therefore, if we compare our results for N 1,2 = 50 with Fig. 3(a) in Lage (2011), we see that EQMOM has similar accuracy for m 0 . For higher-order moments (k 1), EQMOM gives better predictions of the moments than DuQMoGeM by at least one-half order of magnitude.

[Figure 2 about here.]

Finally, the e↵ect of n on the relative moment errors is shown in Fig. 3. Figure 3(a) for n = 1 shows abnormal behavior compared with larger values of n. Indeed, the relative moment error reaches unity at about t = 1.3 for n = 1. The reason is that only one beta function does not approximate the NDF well enough, and leads to faster evaporation than expected. This example shows that using a simple NDF may not be a good choice to simulate evaporation. Comparing the other three figures, we observe that all relative moment errors decrease by about one order of magnitude by increasing the node number from 2 to 3.

Increasing n from 3 to 4, the relative moment errors show about one-half order of magnitude improvement.

As expected, these comparisons show that the error in the first quadrature dominates when n is small, but with three or four nodes the beta EQMOM can capture accurately the moments of the NDF for evaporation.

[Figure 3 about here.]

Case 2. Figure 4 compares the predicted to the exact NDF at t = 1 for Case 2 with n = 4 and di↵erent values of N ↵ . We should note that for Case 2, the exact NDF corresponds to a beta EQMOM with n = 1 ( = ⇠ 1 ), and thus our focus here is on the errors incurred by choosing n > 1. As can be seen from Fig. 4, a significant improvement in the predicted NDF is obtained by increasing N 1 and N 2 from 8 to 80. As in Case 1, N 3 and N 4 can be set as small as 5 without changing the accuracy of the predicted moments. If we compare Fig. 4(b) in [START_REF] Lage | On the representation of QMOM as a weighted-residual method -The dual-quadrature method of generalized moments[END_REF] to Fig. 4, we can clearly see that DuQMoGeM yields a negative NDF in the interval [0.5, 1] at t = 1, which is due to the polynomial approximation used in DuQMoGeM. In comparison, EQMOM always yields a realizable NDF. For this reason, the relative moment errors with DuQMoGeM are significantly larger than with EQMOM, as is obvious by comparing Fig. 5 with Fig. 4(a) in [START_REF] Lage | On the representation of QMOM as a weighted-residual method -The dual-quadrature method of generalized moments[END_REF]. For example, EQMOM with the same parameters (i.e. Fig. 5(a)) as DuQMoGeM gives one order of magnitude smaller error for m 0 , and more than two orders of magnitude smaller errors for the higher-order moments.

Comparing Fig. 5(a) to Fig. 5(b), it is interesting to note that the errors in the predicted NDF in Fig. 4 do not cause large relative moment errors.

[Figure 4 about here.]

[Figure 5 about here.]

Case 2 was also considered in [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF] using the EM method to reconstruct f (t, ⇠) on the finite interval [0, 1]. Although the EM method always guarantees positivity of the NDF, from Fig. 6 (t = 0.5) in [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF], it can be seen that the predicted NDF does not agree very well with the exact NDF. For this reason, the relative moment errors found with EM methods are significantly larger than with EQMOM.

For this case, the EM method could be improved by reconstructing the NDF on the interval [0, b(t)].

Case 3. Figure 6 presents results for the NDF predicted using beta EQMOM with n = 4. In this case, the larger surface area droplets evaporate faster than smaller droplets, which yields the higher values of f (t, 0) at intermediate times as compared to the NDF found with constant g. From Fig. 6, it can be observed that the EQMOM solution matches very well with the exact NDF until very few droplets remain, such as at time t = 1 in Fig. 6. Nevertheless, the predicted NDF always remains realizable as t ! 1.

[Figure 6 about here.]

Figure 7 presents the relative moment errors for Case 3. However, in order to compare with Fig. 4(left)

in [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF], the relative moment errors in Fig. 7 are defined as

m k,error = |m k,EQMOM m k,exact | m k (0) for k = 0, . . . , 2n + 1; ( 36
)
where m k (0) are the initial moments. From Fig. 7(a) we can observe that, using the same number of nodes as in [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF], the relative moment errors found with EQMOM are one order of magnitude smaller than with the EM method before t = 0.6. For larger t, the higher-order moments become very small, and a beta kernel density function may not be the optimal choice for approximating the NDF (see Fig. 6). Nevertheless, beta EQMOM still achieves the same order of accuracy as the EM method. By using four nodes, as shown in Fig. 7(b), beta EQMOM achieves more than one and half orders of magnitude better accuracy as compared with the EM results in [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF] for all t.

[Figure 7 about here.]

In summary, the results for Cases 1-3 with pure evaporation clearly illustrate the favorable properties of EQMOM to handle purely hyperbolic problems with negative flux terms at the boundary of phase space.

Not only is the predicted NDF guaranteed to be realizable for all choices of the quadrature parameters, but

we have also seen that by suitably choosing n and N ↵ the errors due to the first and second quadratures can be reduced to levels that ensure accurate predictions for the moments of the NDF. In fact, in comparison to Lagrangian droplet tracking methods that treat evaporation by simulating an ensemble of droplets, EQMOM should be able to attain equivalent or better accuracy simply by increasing N 1 and N 2 (which represent "notional droplets" in the discretized NDF given in Eq. ( 33)).

Condensation problems

For condensation (positive g(t, ⇠)), the phase-space variable ⇠ is chosen to be particle volume. As in the literature [START_REF] Lage | On the representation of QMOM as a weighted-residual method -The dual-quadrature method of generalized moments[END_REF][START_REF] Strumendo | Solution of PBE by MOM in finite size domains[END_REF], we consider the following two cases with known analytical solutions.

• Case 4: g(t, ⇠) = ⇠/2 with initial condition f (0, ⇠) = 6⇠ 3 e ⇠ for ⇠ 2 [0, 1) [START_REF] Lage | On the representation of QMOM as a weighted-residual method -The dual-quadrature method of generalized moments[END_REF].

• Case 5: g(t, ⇠) = K/⇠ with initial condition f (0, ⇠) = (2 p /7 q ) (⇠ 1) q (15 ⇠) q for ⇠ 2 [1, 15] where p = 2, q = 8, and K = 0.78 [START_REF] Strumendo | Solution of PBE by MOM in finite size domains[END_REF].

For Case 4, gamma EQMOM is used to approximate the NDF, while beta EQMOM is used for the bounded interval in Case 5. As with pure evaporation, the boundaries of the interval [a(t), b(t)] are found using the ODEs a 0 = g(t, a) and b 0 = g(t, b). The reader can note that the initial NDF in Case 4 can be exactly represented by gamma EQMOM with n = 1 ( = 1 and ⇠ 1 = 4). Likewise, in Case 5, the initial NDF is a transformed beta distribution with parameters = µ = p 1.

Case 4. The exact NDF for Case 4 [START_REF] Lage | On the representation of QMOM as a weighted-residual method -The dual-quadrature method of generalized moments[END_REF] is

f (t, ⇠) = (⇠e t/2 ) 3 e ⇠e t/2 6e t/2
, which again can be represented exactly using gamma EQMOM with n = 1, = e t/2 , and ↵ = 4. We fix n = 4, which means N ↵ = 5, and set the quadrature parameters in the adaptive Wheeler algorithm [START_REF] Yuan | Conditional quadrature method of moments for kinetic equations[END_REF] to rmax(1) = 0, rmax(2) = 10 6 , rmax(3) = 10 5 , and rmax(4) = 10 4 . Thus, only one quadrature node will be used if the NDF approximate is exact. As shown in Fig. 8 for t = 10, the gamma EQMOM algorithm correctly reduces to n = 1 and sets w ↵ = 0 for ↵ = 2, 3, 4. In fact, for Case 4, gamma EQMOM yields the exact NDF for all times.

[Figure 8 about here.]

Case 5. For Case 5, the growth rate corresponds to di↵usion-controlled growth [START_REF] Strumendo | Solution of PBE by MOM in finite size domains[END_REF] and the exact NDF is

f (t, ⇠) = f ⇣ 0, p ⇠ 2 2Kt ⌘ ⇠ p ⇠ 2 2Kt
.

Note that the boundaries on the exact NDF are a(t) = p 1 + 2Kt and b(t) = p 15 2 + 2Kt. Figure 9 shows the predicted NDF using beta EQMOM with 7 moments (n = 3) as compared with the exact solution at t = 20. The agreement is excellent and the two solutions coincide within machine precision with each other. Next, we can compare the beta EQMOM result with Fig. 1 in [START_REF] Strumendo | Solution of PBE by MOM in finite size domains[END_REF], which uses the FCMOM method to reconstruct the NDF. With almost the same number of moments controlled (7 in EQMOM, 8 in FCMOM), beta EQMOM predicts the NDF much better. In fact, since FCMOM is just a special case of DuQMoGeM, unavoidably, it gives negative values in the tail of the NDF. In contrast, the NDF reconstructed with EQMOM is always nonnegative.

[Figure 9 about here.]

In summary, for pure growth problems EQMOM provides a very accurate closure for the NDF. In fact, since the growth is hyperbolic in phase space, the EQMOM solution method using characteristics is essentially exact. Thus, the error in the NDF is associated mainly with the first quadrature representation in terms of n kernel density functions. Obviously, for very complex shapes, larger values of n will be needed to accurately capture the NDF. Nevertheless, the EQMOM algorithm provides an accurate and robust method for advancing the moments in time for hyperbolic problems with complex growth rates.

Breakage problems

For pure breakage problems, by convention, the phase variable is chosen to be the particle volume. The moments evolve according to Eq. ( 35), which is solved numerically using the second-order Runge-Kutta method with fixed time step t = 0.001. For convenience, we assume that the NDF is nonzero for ⇠ 2 [0, 1] and use beta EQMOM to reconstruct the NDF. The following three cases with analytical solutions [START_REF] Lage | On the representation of QMOM as a weighted-residual method -The dual-quadrature method of generalized moments[END_REF][START_REF] Zi↵ | The kinetics of cluster fragmentation and depolymerization[END_REF] are analyzed:

• Case 6: #(t, ⇠) = 2 (binary breakage), P (t, ⇠|⇠ 0 ) = H(⇠ 0 ⇠)/⇠ 0 where H is the Heavyside function, b(t, ⇠) = ⇠ 2 and Q(t, ⇠) = [2⇠ 2 (2 e t ) 2(1 e t )]/f (t, ⇠) with initial condition f (0, ⇠) = 6⇠ 3 e ⇠ .
• Case 7:

#(t, ⇠) = 2, P (t, ⇠|⇠ 0 ) = H(⇠ 0 ⇠)/⇠ 0 , b(t, ⇠) = ⇠ 1/3 and Q(t, ⇠) = [7e t 12 + 7(2 e t )⇠ 1/3 ]/f (t, ⇠) with initial condition f (0, ⇠) = 6⇠ 3 e ⇠ .
• Case 8:

#(t, ⇠) = 2, P (t, ⇠|⇠ 0 ) = H(⇠ 0 ⇠)/⇠ 0 , b(t, ⇠) = ⇠ 2 and Q(t, ⇠) = 0 with initial condition f (0, ⇠) = (⇠ 1).
Note that the source term Q(t, ⇠) in Cases 6 and 7 is chosen to yield a known form for f (t, ⇠), and thus

Q(t, ⇠
) is a known function of t and ⇠.

Case 6. For this case, the analytical solution [START_REF] Lage | On the representation of QMOM as a weighted-residual method -The dual-quadrature method of generalized moments[END_REF]

is f (t, ⇠) = 2 e t , which is constant for ⇠ 2 [0, 1].
Here, the breakage kernel involves only integer powers of ⇠ and the highest-order is 2n + 2. Therefore, N ↵ must satisfy 2N ↵ 2n + 3, and we set n = 2 and N ↵ = 4. Sample results for Case 6 are shown in Fig. 10.

In Fig. 10(a), the reconstructed NDF is on top of the exact NDF. As shown in [START_REF] Lage | On the representation of QMOM as a weighted-residual method -The dual-quadrature method of generalized moments[END_REF], the abscissas must be constant for all time, but standard QMOM cannot preserve this property, even though the moments are predicted quite well. As can be seen from Fig. 10(b), having constant abscissas is not a problem for EQMOM.

[Figure 10 about here.]

Case 7. For this case, the analytical solution [START_REF] Lage | On the representation of QMOM as a weighted-residual method -The dual-quadrature method of generalized moments[END_REF]

is f (t, ⇠) = 2 e t , which is constant for ⇠ 2 [0, 1].
Here, the breakage kernel b(⇠) = ⇠ 1/3 cannot be represented by integer moments, so larger N ↵ has to be used to reduce the integration error. Figure 11 shows the e↵ect of N ↵ on the reconstructed NDF, and it can be observed that a much better prediction is obtained by increasing N ↵ from 25 to 100. The same dependence on N ↵ is found for the relative moment errors. From Fig. 12, we can see that the first-order moment error decreases as N ↵ becomes larger. Comparing this figure with Fig. 12(a) in [START_REF] Lage | On the representation of QMOM as a weighted-residual method -The dual-quadrature method of generalized moments[END_REF] for the same values of n and N ↵ , we observe that EQMOM gives more then one and a half order of magnitude smaller moment errors as compared to DuQMoGeM.

[Figure 11 about here.]

[Figure 12 about here.]

Case 8. For this case, the analytical solution [START_REF] Zi↵ | The kinetics of cluster fragmentation and depolymerization[END_REF] is

f (t, ⇠) = e t⇠ 2 [ (⇠ 1)+2tH(1 ⇠)]
for ⇠ 2 [0, 1]. The principal di↵erence with Case 6 is the initial condition, so N ↵ is again set to n + 2. (2011), it is clear that EQMOM gives much better predictions for the higher-order moments in Case 8.

[Figure 13 about here.]

[Figure 14 about here.]

In summary, the results for pure breakage confirm the excellent predictive capability of EQMOM for solving moment equations with linear integro-di↵erential terms. When the breakage kernel is a polynomial of finite order, N ↵ can be chosen large enough to exactly evaluate the second quadrature. Otherwise, N ↵ can be used as a parameter to reduce the error in the second quadrature to a level that is smaller than the error in the first quadrature (which is controlled by n).

Aggregation/coalescence problems

For aggregation and coalescence problems, the size variable is volume. The following three cases are analyzed [START_REF] Ernst | Coagulation processes with a phase transition[END_REF][START_REF] Gelbard | Numerical solution of the dynamic equation for particulate systems[END_REF][START_REF] Lage | On the representation of QMOM as a weighted-residual method -The dual-quadrature method of generalized moments[END_REF][START_REF] Vemury | Self-preserving size distributions of agglomerates[END_REF]:

• Case 9: a(t, ⇠, ⇠ 0 ) = ⇠ + ⇠ 0 (sum kernel) with initial condition f (0, ⇠) = e ⇠ .
• Case 10: a(t, ⇠, ⇠ 0 ) = ⇠⇠ 0 (product kernel) with initial condition f (0, ⇠) = 4e 2⇠ .

• Case 11: a(t, ⇠, ⇠ 0 ) = ⇠ 1/3 + ⇠ 01/3 1/⇠ 1/3 + 1/⇠ 01/3 (Brownian kernel) with initial condition f (0, ⇠) = e ⇠ .

Gamma EQMOM is employed to solve Eq. ( 35).

Case 9. The analytical solution [START_REF] Gelbard | Numerical solution of the dynamic equation for particulate systems[END_REF] for this case is

f (t, ⇠) = e t 2⇠+⇠e t ⇠ p 1 e t I 1 ⇣ 2⇠ p 1 e t ⌘
where I 1 is a modified Bessel function. For Case 9, the aggregation kernel involves only integer powers of ⇠ and the highest order is 2n + 1, so N ↵ is chosen such that N ↵ n + 1. The distribution profile at t = 3 is shown in Fig. 15 and EQMOM gives a reasonable approximation of this complex distribution function.

Table 1 shows the relative moment errors at t = 3. For aggregation, the first-order moment is constant, and m 1,error is zero to machine precision. The relative moment errors of the higher-order moments are slightly larger than for the lower-order ones, but all moment errors are quite small (⇡ 10 5 ).

[Figure 15 about here.]

[Table 1 about here.]

Case 10. For this case, to our knowledge, no analytical solution for the NDF exists, but an analytical solution can be found for the first three moments [START_REF] Ernst | Coagulation processes with a phase transition[END_REF]:

m 0,exact (t) = m 0 (0) t/2, m 1,exact (t) = 1, m 2,exact (t) = m 2 (0) 1 m 2 (0)t ,
where m k (0) is kth-order initial moment. Because the aggregation kernel is a polynomial, N ↵ is set to

N ↵ = n + 1.
For the initial condition f (0, ⇠) = 4e 2⇠ , m 2 (0) = 1, and the critical time is t = 1. Therefore, the simulation time is set as t = 0.99. Due to conservation of the first-order moment, the error for m 1 (t) is zero. For this reason, only two moment errors are plotted in Fig. 16. It can be seen that EQMOM gives very accurate predictions for the product kernel up to very close to the critical time, where m 2 reaches infinity.

[Figure 16 about here.]

Case 11. For the Brownian aggregation kernel, no analytical solution for the NDF exists, but a selfpreserving solution is found for large t. In order to get this solution, a dimensionless volume ⌘ and a dimensionless NDF are defined as

⌘ = m 0 m 1 ⇠, = m 1 m 2 0 f.
The aggregation kernel a(t, ⇠, ⇠ 0 ) cannot be represented by integer moments, so N ↵ = 100 is set to reduce the integration error. Figure 17 shows the EQMOM result compared with the self-preserving solution from [START_REF] Vemury | Self-preserving size distributions of agglomerates[END_REF]. It can be observed that the EQMOM result matches very well with the literature solution.

[Figure 17 about here.]

Coupled phenomena

The last cases involve either simultaneous breakage and aggregation [START_REF] Lage | On the representation of QMOM as a weighted-residual method -The dual-quadrature method of generalized moments[END_REF][START_REF] Mccoy | Analytical solution for a population balance equation with aggregation and fragmentation[END_REF], or simultaneous evaporation and coalescence:

• Case 12: Breakage and aggregation with (t, ⇠) = 2, P (t,

⇠|⇠ 0 ) = H(⇠ 0 ⇠)/⇠ 0 , b(t, ⇠) = 2 (1)⇠/2, a(t, ⇠, ⇠ 0 ) = 1, Q(t, ⇠) = 0 and (1) = 2 with initial condition f (0, ⇠) = e ⇠ .
• Case 13: Evaporation and coalescence with g(t,

⇠) = k e ⇠ 1/3 , a(t, ⇠, ⇠ 0 ) = k c ⇠ 1/3 + ⇠ 01/3 1/⇠ 1/3 + 1/⇠ 01/3 with initial condition f (0, ⇠) = ⇠ 2 e ⇠ .
Gamma EQMOM can be used to simulate these cases with ⇠ equal to particle volume.

Case 12. The analytical solution [START_REF] Mccoy | Analytical solution for a population balance equation with aggregation and fragmentation[END_REF] for this case is f (t, ⇠) = 2 (t)e (t)⇠ where (t) = (1) 1 + (1) tanh( (1)t/2) (1) + tanh( (1)t/2) .

This solution is just a gamma distribution with the parameter depending on t. Thus, it can be expected that one-node gamma EQMOM su ces to capture the exact NDF and, because the kernel is a polynomial, we set n = 1 and N 1 = 2. The EQMOM results for this case are shown in Fig. 18. As we can see, the first quadrature finds the correct parameters (i.e., ⇠ ↵ , w ↵ and ) so that the EQMOM solution is almost on top of the exact NDF for all t.

[Figure 18 about here.]

Case 13. In the final case, a simultaneous evaporation and coalescence case is used to illustrate that EQMOM also works well with strongly coupled phenomena. To our knowledge, no analytical solution exists for this case. The parameter for the second quadrature is set to N 1,2 = 200 in order to get a relatively small time step for evaporation using the methods introduced in Sec. 4.4. The same time step is adopted for the coalescence term. Figure 19 shows how evaporation and coalescence a↵ect the NDF at t = 10. It can be observed that the NDF for k c = 0.05 and k e = 1 lies mostly on the left-hand side of the initial NDF because evaporation is stronger than coalescence. Increasing k c to 0.25, the coalescence dominates and moves the NDF toward the pure coalescence NDF. Overall, the coupling between growth and coalescence poses no particular problems for EQMOM.

[Figure 19 about here.]

Conclusions

The extended quadrature method of moments (EQMOM) presented in this work represents a significant improvement over existing methods for solving population balances equations (PBE) in the context of moment methods. The proposed NDF reconstruction algorithm is robust, computationally e cient, and results in a continuous, non-negative NDF that always reproduces 2n moments, and in many cases 2n + 1. By choosing a kernel density function that is equivalent to the weight function of a known family of orthogonal polynomials, we have demonstrated that EQMOM leads to a second Gaussian quadrature that can be employed to greatly improve the accuracy with which the source terms in the PBE can be estimated. Furthermore, by using operator splitting, the second quadrature can be used as a point-particle representation to e ciently and accurately solve hyperbolic terms in the PBE (such as evaporation and condensation). Using thirteen test cases, we demonstrate unequivocally the applicability and accuracy of EQMOM for approximating the moments of the PBE for cases of interest in aerosol science, as well as its relative strengths and shortcomings vis-a-vis other methods.

Our current research is aimed at generalizing EQMOM in two directions: (i) a multi-variate version of EQMOM based on an extension of the CQMOM algorithm described in Yuan and Fox (2011), and (ii) the application of EQMOM to evaluate the spatial fluxes (i.e., advection and di↵usion) for cases where the advection velocity and di↵usion coe cient depend on the internal coordinates (⇠). For the latter, the realizable high-order finite-volume schemes described in [START_REF] Vikas | Realizable high-order finite-volume schemes for quadrature-based moment methods[END_REF] can be applied with the pointparticle representation from EQMOM to guarantee that the transported moment set is always realizable.

smallest for which J n ( ) = 0 can sometimes yield an abscissa out of range. In such cases, we choose to be the largest value for which the Hankel determinants (defined below) of the moment set (m ⇤ 0 , . . . , m ⇤ 2n 1 ) are nonnegative [START_REF] Dette | The Theory of Canonical Moments with Applications in Statistics, Probability, and Analysis[END_REF][START_REF] Wright | Numerical advection of moments of the particle size distribution in Eulerian models[END_REF]. In this manner, the moment set (m ⇤ 0 , . . . , m ⇤ 2n 1 ) is always realizable (i.e., all w ↵ > 0 and ⇠ ↵ 2 [0, 1] (beta) or ⇠ ↵ 2 [0, 1) (gamma)).

Hankel determinants for beta EQMOM. The Hankel determinants for beta EQMOM are defined by

H 2i = m ⇤ 0 . . . m ⇤ i . . . . . . m ⇤ i . . . m ⇤ 2i , H 2i+1 = m ⇤ 0 m ⇤ 1 . . . m ⇤ i m ⇤ i+1 . . . . . . m ⇤ i m ⇤ i+1 . . . m ⇤ 2i m ⇤ 2i+1 , H 2i+1 = m ⇤ 1 . . . m ⇤ i+1 . . . . . . m ⇤ i+1 . . . m ⇤ 2i+1 , H 2i = m ⇤ 1 i ⇤ 2 . . . m ⇤ i m ⇤ i+1 . . . . . . m ⇤ i m ⇤ i+1 . . . m ⇤ 2i 1 m ⇤ 2i . (B.3)
In practice, H and H are computed for i = 0, . . . , n 1; and J n ( ) is equal to the smallest Hankel determinant (such that all determinants are non-negative). We then choose to be the largest value for which J n ( ) = 0.

Hankel determinants for gamma EQMOM. The Hankel determinants for gamma EQMOM are defined [START_REF] Wright | Numerical advection of moments of the particle size distribution in Eulerian models[END_REF] by

i,j = m ⇤ i m ⇤ i+1 . . . m ⇤ i+j m ⇤ i+1 m ⇤ i+2 . . . m ⇤ i+j+1 . . . . . . . . . . . . m ⇤ i+j m ⇤ i+j+1 . . . m ⇤ i+2j . (B.4)
In practice, the i,j for i = 0, 1 and j = 0, . . . , n 1 have to be non-negative; and J n ( ) is equal to the smallest Hankel determinant.

  is the indicator function for the interval [a, b], and ⇢ ↵ = w ↵ w ↵ . The n weights w ↵ , n abscissas ⇠ ↵ and are found from the first 2n + 1 integer moments of f (t, ⇠) using EQMOM (which we refer to as the first quadrature).

  [0, b(t)], where b(t) is found by solving the ODE db dt = g(t, b) with b(0) = 1. The quadrature parameters n and N ↵ are varied to investigate the dependence of the moment errors on the quadrature reconstruction. Recall from Sec. 4.4 that the time step is determined by the value of N ↵ , as this parameter fixes the location of the two abscissas closest to the origin. Hence, by increasing N ↵ , the accuracy of the evaporative flux prediction should improve. Unless stated otherwise, the relative moment errors are defined by

Figure 13 (

 13 Figure 13(a) shows the reconstructed NDF at t = 10 for n = 3. The overall agreement between EQMOM and the exact NDF is very good for n = 3, except near the lower bound. By increasing n to 4 (Fig. 13(b)),

Figure 2 :Figure 3 :Figure 5 :Figure 6 :Figure 7 :Figure 8 :

 235678 Figure1: Beta EQMOM results for the time evolution of the NDF in Case 1 at selected times t with n = 4. Symbols: beta EQMOM. Lines: exact NDF.

Figure 9 :Figure 10 :Figure 11 :Figure 12 :Figure 14 :

 910111214 Figure 9: Beta EQMOM results for the NDF in Case 5 at t = 20 with n = 3 and N↵ = 5. Beta EQMOM (symbols). Exact NDF (line).

Figure 15 :

 15 Figure 15: Reconstructed NDF at t = 3 for Case 9. Symbols: gamma EQMOM with n = 4 and N↵ = 5. Line: exact NDF.

Figure 16 :Figure 17 :Figure 18 :

 161718 Figure16: Relative moment errors for Case 10 found using gamma EQMOM with n = 1 and N 1 = 2. m 0 : solid. m 2 : dashed.

  error 5 ⇥ 10 7 4 ⇥ 10 6 2 ⇥ 10 5 3 ⇥ 10 5 6 ⇥ 10 5 9 ⇥ 10 5 1 ⇥ 10 5 2 ⇥ 10 5

Table 1 :

 1 Relative moment errors at t = 3 for Case 9.

The weights w ↵ are normalized such that P w ↵ = 1.

In the literature[START_REF] Lage | On the representation of QMOM as a weighted-residual method -The dual-quadrature method of generalized moments[END_REF] a case with initial conditions given by the sum of two delta functions is also considered. With EQMOM, this case is trivial since it can be solved exactly with n = 2.

Acknowledgments

This work was partially supported by a grant from the US National Science Foundation (CCF-0830214).

One of the authors (R.O.F.) wishes to thank the Ecole Centrale Paris for hosting his stay in France during the 2010-11 academic year.

Appendix A. Comparison of EQMOM with EM method

Unlike EQMOM, the entropy maximization (EM) method [START_REF] Mead | Maximum entropy in the problem of moments[END_REF] can always reconstruct the NDF using the entire moment set (m 0 , m 1 , . . . , m 2n ) when it is realizable. In Fig. 20 we compare two example NDF reconstructed with beta EQMOM and EM, and it is evident that both methods generate similar shapes. In terms of computational e ciency, our reconstruction algorithm for EQMOM is approximately 100 times faster than EM for the same value of n. However, even if the computational times for both methods were similar, the evaluation of integrals using the second quadrature with EQMOM will be many times faster and more accurate than attempting to construct a Gaussian quadrature from the NDF found using EM.

[Figure 20 about here.] Appendix B. Realizability for beta and gamma EQMOM For fixed n, the moment set (m 0 , m 1 , . . . , m 2n ) is used by EQMOM to reconstruct the NDF. In this appendix, we briefly describe how the moment-inversion algorithm in EQMOM treats cases for which the chosen kernel density function can not recover the highest-order moment m 2n by varying with fixed n. By construction, all moments up to m 2n 1 are always reproduced by EQMOM and thus the ability to capture m 2n depends on the shape of the kernel density function. This is most easily understood by considering the canonical moments.

Appendix B.1. Canonical moments

Consider the set of moments (m 0 , m 1 , . . . , m 4 ) found from a realizable NDF on the interval [0, 1]. The 4-dimensional moment space associated with this moment set is convex, but its boundaries are not easily expressed in terms of the moments. However, a set of canonical moments (p 1 , . . . , p 4 ) can be defined [START_REF] Dette | The Theory of Canonical Moments with Applications in Statistics, Probability, and Analysis[END_REF] to map the moment space into the 4-dimensional hypercube [0, 1] 4 . Thus, every point in the hypercube corresponds to a realizable NDF.

With beta EQMOM, the moments (m ⇤ 0 , m ⇤ 1 , . . . , m ⇤ 3 ) depend on and can be expressed in terms of the canonical moments:

Using the bounds of the canonical moments, the expression for p ⇤ 2 places an upper bound on :  max = p 2 /(1 p 2 ). Likewise, the bounds on p ⇤ 3 introduce additional constraints:

Note that as the value of n in EQMOM is increased, the upper bound on becomes tighter.

By definition [START_REF] Dette | The Theory of Canonical Moments with Applications in Statistics, Probability, and Analysis[END_REF], the QMOM reconstruction has p ⇤ 2n = 0. Thus, the EQMOM reconstruction attempts to increase until p 2n ( ) equals the value of the 2n th canonical moment found with the true m 2n , which we denote by p 2n . Theoretically, any value 0  p 2n  1 corresponds to a realizable NDF. However, max 2[0, max ] p 2n ( ) is often less than unity, implying that there exists realizable moment sets for which beta EQMOM can not be made to match m 2n (even though all lower-order moments are exactly reproduced). One possible remedy to this problem would be to use a more general kernel density function for which 0  p 2n ( )  1. However, it would then be necessary to find explicit formulas for the family of orthogonal polynomials corresponding to this new weight function, and it is likely that the computational cost would dramatically increase. Alternatively, we can simply adapt the moment-inversion algorithm to minimize the di↵erence between p 2n and p 2n ( ), and increase n when more moments are needed to improve the accuracy of the reconstructed NDF.

Appendix B.2. Treatment of realizability with EQMOM

From the discussion above, it is inevitable that some values of will generate unrealizable moments.

Indeed, when we implement the EQMOM moment-inversion algorithm described in Sec. 3.5, we find that the