Coagulation, non-associative algebras and binary trees

We consider the classical Smoluchowski coagulation equation with a general frequency kernel. We show that there exists a natural deterministic solution expansion in the non-associative algebra generated by the convolution product of the coalescence term. The non-associative solution expansion is equivalently represented by binary trees. We demonstrate that the existence of such solutions corresponds to establishing the compatibility of two binary-tree generating procedures, by: (i) grafting together the roots of all pairs of order-compatibile trees at preceding orders, or (ii) attaching binary branches to all free branches of trees at the previous order. We then show that the solution represents a linearised flow, and also establish a new numerical simulation method based on truncation of the solution tree expansion and approximating the integral terms at each order by fast Fourier transform. In particular, for general separable frequency kernels, the complexity of the method is linear-loglinear in the number of spatial modes/nodes.


Introduction
Herein we consider the classical Smoluchowski coagulation equation with a general frequency kernel.We show there exists a natural deterministic solution expansion in a non-associative algebra equivalent to an expansion on binary trees.The algebra product is generated by the coagulation convolution term and is non-associative for general frequency kernels.We demonstrate how the existence of solutions to the coagulation equation is equivalent to the compatibility of generating binary trees at each order, by the two following natural procedures, by: (i) grafting together the roots of all pairs of order-compatible trees at preceding orders, and (ii) attaching binary branches to all the possible free branches of trees at the previous order.We examine the classical constant, additive and multiplicative frequency kernels as special cases within the binary tree context.Lastly, we also establish that the flow has a linearised flow formulation, and outline how the flow can be approximated to develop new numerical simulation methods.
The Smoluchowski coagulation equation has the form, ∂ t g(x; t) = 1 2 x 0 K(x − y, y)g(x − y; t)g(y; t) dy − g(x; t) ∞ 0 K(x, y)g(y; t) dy, where g = g(x; t) denotes the density of molecular clusters of mass x, and K = K(x, y) is the given symmet-Email address: S.J.A.Malham@hw.ac.uk (Simon J.A. Malham) 1 SJAM was supported by an EPSRC Mathematical Sciences Small Grant EP/X018784/1 ric frequency kernel.We assume the initial data is g 0 so g(x; 0) = g 0 (x).All the variables concerned here are positive.The exact form of the frequency kernel depends on the application considered.Three special cases are often the focus of investigation, namely the constant, K = 1, additive, K = x + y, and multiplicative, K = xy, kernel cases.This is because the Smoluchowski equation can be solved explicitly for these cases.See Aldous [2] and Menon and Pego [77] for more details.In the multiplicative kernel case a phase transition occurs at a gelation time, beyond which the solution can be extended, see Leyvraz and Tschudi [65].Herein, we focus on the general K = K(x, y) kernel case.It is usual to consider the Laplace or Bernstein transform of the Smoluchowski's equation; see Menon and Pego [77].For the moment let us consider a general linear transformation H analogous to such transforms.We define the general linear transformation of a function g = g(x) with support on [0, ∞) by, where 0 h 1.If h := e −sx then H is the Laplace transform.If h(s, x) := 1 − e −sx , then H is the Bernstein transform.It is straightforward to show that if we consider the general transform of the Smoluchowski equation (1) and g := H g, then g = g(s; t) satisfies, ∂ t g(s; t) = g(y; t) H(s, y, z)K(y, z) g(z; t) dy dz, (3) where the double integral shown is over [0, ∞) 2 , and we set H(s, y, z) := 1 2 (h(s, y + z) − h(s, y) − h(s, z)).If we set g = H −1 g in both places where g appears on the right in (3), the resulting vector field involves a quadruple integral of a quadratic form involving g with a weight H(s, y, z)K(y, z)v(y, •)v(z, •).Here v is the kernel of the operator H −1 .So, for example, in the instance of the Laplace transform v(s, x) = e xs , and the quadruple integral additionally involves two Bromwich contour integrals.The resulting evolution equation for g has the form, where the product '⋆' is based on the quadratic form of the vector field in the evolution equation (3) with the quadruple integral.There are many ways to represent Smoluchowski's coagulation equation as an abstract evolution equation of the form (4) with the quadratic vector field shown.We consider another analogous representation for Smoluchowski's coagulation equation in Section 3. In all of these representations which have the abstract form (4), the product '⋆' is in general non-associative.That this is the case is straightforwardly checked.The constant frequency kernel case is the exception, the product '⋆' is associative in this singular instance.Naturally the form of f ⋆ g for general transforms of any pair of functions f and g with support on [0, ∞), for any of the representations mentioned, is straightforwardly implied.Hence, with '⋆' a non-associative product, our goal is to solve the abstract equation ( 4), or equivalently, for g = g(s; t), where ξ = ξ(s) is the general transform of the data g 0 = g 0 (x).The solution can be formally derived by iteration.This generates the solution expansion, In the context of Smoluchowski's equation, the product '⋆', though non-associative, is commutative.This means that we can simplify the expansion by combining some like terms.For example the terms of order t 2 combine, as do all the terms of order t 3 except the symmetric term with the real factor '2'.We discuss such symmetries in detail in Section 3.However, in general we assume the product '⋆' is non-commutative as we can complete all our analysis in this more general context, and specialise to the commutative context as and when we require.Keeping track of the bracketed terms due to the non-associativity at higher order is most easily accomplished using the representation of the terms in the expansion by rooted planar binary trees.
Besides the empty tree ∅, the first set of rooted planar binary trees, up to and including grade 3, are, , Respectively, in order, these trees encode the successive terms in our formal solution expansion (6) for g.Each vertex indicates a '⋆' product while each free branch indicates the instance of a ξ factor.Hence the first term 'ξ' in (6) is represented by the ' ', which, say, we express as ξ( ).The second term 'ξ ⋆ ξ' is represented by ' ' and we express this using the shorthand notation ξ( ), and so forth.Indeed a representation for g in ( 6) is thus, We observe, if we drop the ξ's, we can represent this solution expansion for g as an expansion in the algebra of planar binary trees over the field R.
The purpose of this paper is to: (i) Prove the full tree expansion (6) is the solution to the abstract non-associative and non-commutative quadratic differential equation ( 4) with data ξ.The proof is equivalent to establishing the compatability, with mutliplicity, of the grafting and branching tree generating procedures; (ii) Give closed form expressions for the solution; (iii) Show that the solution can be expressed as a linearised flow; (iv) Establish a practical numerical simulation method for the Smoluchowksi coagulation equation, with a general frequency kernel, to evaluate the solution at any time t, up to gelation.The method is particularly efficient for separable frequency kernels.
There are many comprehensive surveys on Smoluchowksi coagulation, more recent examples include Aldous [2], giving a general overview, da Costa [25], focusing on deterministic aspects and Hammond [50], focusing on stochastic aspects.Conditions for well-posedness of the Smoluchowski equation, the gelation time, as well as the phenomenon of instantaneous gelation can be found therein, as well as in Dubovskii [31].The work of Menon and Pego [77] and Iyer et al. [57] in particular, establishes a general rigorous foundation for the Smoluchowski equation and its extension to multiple mergers.On the theory of differential equations in the context of non-associative algebras, early work includes that by Markus [72], as well as Röhrl [84].
The algebra of rooted planar trees is exceedingly rich and has an illustrious history dating back to Cayley.More recently, in the work of Butcher [14], they arose in the context of Runge-Kutta methods in numerical analysis in the form of the Butcher group and in the work of Connes and Kreimer [23] on perturbative quantum field theory in the form of a Hopf algebra of trees.These two objects are equivalent, as demonsrtated by Brouder [12].A Hopf algebra is an algebra endowed with a compatabile co-product, and an antipode, see Foissy [42] for an introduction.In the algebra context, the Hopf algebra of planar binary trees was established by Loday and Ronco [66].Also see Hairer et al. [49] for a detailed account of the use of trees in numerical analysis.The connection between Smoluchowski coagulation and binary trees is not new, and indeed binary trees appear, for example, in Aldous [2] in the context of Galton-Watson processes, in the combinatorial approach in Spouge [89], and in Sheth and Pitman [88] as merger history trees.Also see Marckert and Wang [71].
Our paper is structured as follows.In Section 2 we introduce the algebra of planar binary trees.We solve the abstract equation ( 4) in this algebra and present three different equivalent solution forms.We apply our abstract results from Section 2 to coagulation examples in Section 3. In Section 4 we demonstrate that the solution flow can be represented as a linearised flow, and we outline how the binary tree expansion solution can be approximated and used for numerical simulation.Finally, in Section 5 we discuss the many future directions to which our results herein extend.

Non-associative algebras and binary trees
We have seen that the solution (6) to the abstract evolutionary quadratic equation ( 4), representing Smoluchowski's coagulation equation, is formally given as a series in non-associative products '⋆' of the data ξ.The coefficients of the series have a simple time-dependence.We are thus lead to consider the non-commutative, nonassociative algebra with a single generator, and product '⋆'.Such algebras have a representation in the real algebra of planar binary trees.Hence our primary focus herein is on expansions in such algebras.
Let H = H(ξ, ⋆) denote the real non-commutative, nonassociative algebra generated by the function ξ, with the product '⋆'.We do not require H to be unital.The elements of H are polynomials or series of the form, where all the coefficients ĝ1 , ĝ11 , ĝ1 (11) and so forth, are real.Our focus in this section is on the combinatorial structure of such series.In Section 3 we outline in what sense such series converge in H, corresponding to a restriction on the coefficients ĝ1 , ĝ11 , etc.An immediate insight from the form of ( 7) is that the elements of H are characterised by the binary parenthesisation of strings, and linear combinations thereof.The binary parenthesisation of strings is isomorphically mapped to planar binary trees; see Stanley [90,Ch. 6].See Tables 1 and 2 for examples of such trees.As a consequence H is isomorphically mapped to the real algebra of binary trees R T , which we outline in detail presently.The non-commutative, non-associative product '⋆' on H is replaced by the grafting of pairs of trees at their roots-see Procedure 1 just below.Indeed, as already highlighted in the introduction, we can represent series such as (7) in the form, We can thus regard ξ as a homomorphic map from R T to H such that for any pair of trees τ 1 and τ 2 , we have, In the abstract context here, we are seeking solutions of the evolutionary quadratic equation ( 4) in H. However with this homomorphic property of ξ in mind, if we pullback equation ( 4) from H to R T , the result is the following Smoluchowski tree equation on R T .
Definition 1 (Smoluchowski tree equation).We define the Smoluchowski tree equation as the evolutionary quadratic equation for g = g(t) in R T given by, We augment (8) with the data g(0) = .
Remark 1.With a slight abuse of notation we use the same label 'g' to represent the solution to the Smoluchowksi equation in its different guises.It is the solution to (4) when g = g(s; t) is a function, which we consider in more detail in Section 3, or the solution g ∈ R T to the Smoluchowski tree equation (8).It will always be clear from the context, which solution form we refer to.
Our goal herein is to find a solution to initial value problem for the Smoluchowski tree equation ( 8) in Definition 1.For now, we focus on the structure of the real algebra of planar binary trees R T .Here T denotes the set/forest of all rooted planar binary trees.The root is the bottom vertex and the trees grow upwards from there.
Definition 2 (Grade).The grade |τ | of any binary tree τ ∈ T, is the number of vertices it possesses.
Remark 2 (Vertices).We nominate the number of vertices of any tree as the number of internal branchings or 's the tree possesses.Other authors include the number of free branches in the vertex count as well, depending on the context.See Tables 1 and 2 for further clarification.
Two natural procedures for generating rooted planar binary trees are as follows.
Procedure 1 (Grafting).We can generate all the trees of grade n from the trees of grade 0 through n − 1 by rootgrafting as follows.We graft at the root all trees τ 1 and τ 2 such that |τ 1 | + |τ 2 | = n − 1 with |τ 1 | cycling through 0, 1, . .., n − 1.This procedure exhaustively generates all the possible trees of grade n.For example, we can generate all the binary trees of grade 4 in Table 2, using all the trees in Table 1, by: grafting the grade 0 tree ' ' on the left with all the grade 3 trees on the right; grafting the grade 1 tree ' ' on the left with both the grade 2 trees on the right; grafting both grade 2 trees on the left with the grade 1 tree ' ' on the right; and then finally grafting all the grade 3 trees on the left with the grade 0 tree ' ' on the right.
Remark 3 (Grafting and concatentation).As indicated above, one natural product on the real algebra of rooted planar trees R T is grafting in the manner described in Procedure 1; see for example Lundervold and Munthe-Kaas [67].We endow R T with this product.It is akin to a concatenation product of parenthesised strings.
Remark 4 (General grafting).In principle we can graft a planar binary tree τ ′ to any free branch or even vertex of another planar binary tree.However, we do not require this procedure herein.

Remark 5 (Catalan numbers).
There is an intimate link between Catalan numbers and planar binary trees.Indeed, Stanley [90,Ex. 6.19] gives 66 characterisations of Catalan numbers.For n ∈ N 0 , the nth Catalan number is, The Catalan number C n counts the number of planar binary trees with n vertices, in other words, the number of binary trees generated at each grade by the grafting procedure above.Planar binary trees, or variants thereof, are often nominated as Catalan trees.
Procedure 2 (Branching).We can generate all the trees of grade n from the trees of grade n − 1 as follows.
To each tree τ at grade n − 1 we attach a single branch ' ', successively, to each free end of τ .This procedure also exhaustively generates all the possible trees of grade n.However some trees at grade n are multiply generated by those at grade n − 1 by this procedure.The multiplicity of any tree generated by this procedure corresponds to the weight character of that tree; see Definition 3 and Lemmas 1 and 2 below.For example, in Table 1, if we successively add a branch ' ' to each free end of the trees of grade 2, we generate the five binary trees of grade 3 shown, however the middle symmetric tree is generated with multiplicity 2.
Remark 6 (Composition and budding).The procedure for branching outlined is a simple version of the composition of trees; see eg.Lundervold and Munthe-Kaas [67,68].More general trees and linear combinations of trees than just ' ' can be attached to free branches in the manner outlined in Procedure 2. Branching is also related to "budding".See Definition 7 in Section 4 and items (iii) and (v) in Section 5.
Remark 7 (Empty tree).Note, the empty tree ∅ doesn't play a role in either of the tree generating procedures above.Indeed, the empty tree ∅ plays the role of a unit in R T , though we do not explicitly use this here.
We define the weight character associated with any given tree in R T as follows.
Definition 3 (Weight character).We define the weight character ω : R T → Q of any tree τ ∈ R T recursively as follows (by convention ω( ) := 1): and then Thus in essence, to compute the weight character of a given binary tree, we split the tree at its root and compute the product of the weight characters of the two subtrees τ 1 and τ 2 generated and include the Leibniz factor The weight characters of τ 1 and τ 2 are computed by propagating this process.Example weight characters for trees up to grade 4 are given in Tables 1 and 2.
Therein we abbreviate products of Leibniz coefficients by, The weight character naturally arises in the process of tree generation by branching, as described above, and outlined below.See Lemmas 1 and 2.
Definition 4 (Exponential series).We define an exponential tree series in R T as an expansion of the form, where the fτ are the real coefficients associated with each tree τ ∈ T\∅, and ω(τ ) and |τ | are respectively the weight character and grade of τ .
Two further operators are useful to our development.
Definition 5 (Branching operator).We define the linear branching operator B ∨ : R T → R T as the operator that acts on any tree τ ∈ T by successively, additively attaching a branch ' ' to each free end of the tree τ , thus generating a sum of the corresponding trees at the next grade.
Table 1: We list all the rooted planar binary trees up to grade 3.For the grades in the left column, we list each tree τ ∈ T in the second column, its corresponding weight ω(τ ) in the third column, the number of symmetries σ(τ ) it possesses (see Remark 13) in the fourth column and its word-coding in terms of levels in the final column.
grade tree weight symm.levels Remark 8. We can extend the branching operator to B τ ′ , for any τ ′ ∈ T, so that it becomes the operator that successively, additively attaches the tree τ ′ to each free end of the tree τ .We can also linearly extend B τ ′ so that Example 1.Some examples of the action of the branching operator are as follows, we observe, Definition 6 (Grade operator).The grade operator G on R T is given for any τ ∈ T by, Naturally, G is an endomorphism on R T .
Example 2. The grading operator normalises the given tree by its grade, so that for example, Example 3. It is also natural to combine the grading and branching operators, and so we set, B∨ := G • B ∨ .Thus, adapting Example 1, we observe, Naturally, since the non-associative product '⋆' is bilinear, so is the grafting operator.In other words, for example, for any three trees τ 1 , τ 2 , τ 3 ∈ T, we have, The branching operator also acts like a derivation in the following sense.For any two trees τ 1 , τ 2 ∈ T, we observe, A crucial component of our main result herein is the following identity.For convenience, for all n ∈ N 0 , let g n denote the weighted sum of all trees of grade n, so, Lemma 1 (Grafting identity).The following grafting identity holds for all n ∈ N 0 : Proof.Using that the grafting operator is bilinear and the properties of the weight character given in Definition 3, we observe the right-hand side equals, giving the result.
As a consequence, we have the following expansion for any power of the branching operator B ∨ acting on the tree ' '.
Proof.We prove the result by induction.Assume that the result holds for all powers B ℓ ∨ with ℓ = 0, 1, . . ., n.Note that for any of the values ℓ = 0, 1, . . ., n − 1, the statement above is equivalent to the statement, or in other words, B ∨ (g ℓ ) = g ℓ+1 .First, we observe that, using our induction assumption.Second, by direct computation, using Lemma 1, we observe that, Third, using the grafting bilinear property (10) and the derivation property of B ∨ from (11), we have, .
Fourth, substituting for g n−1−k and g k and using the induction assumption (13), we observe , where we shifted the k-summation label in the second term by one.Fifth, carefully considering of the k = 0 case in the first term and the k = n case in the second term in the last line above, and using the identity, we observe that, using Lemma 1.This gives the result for ℓ = n + 1.
We now solve the initial value problem (8) for g = g(t) in R T , using Lemmas 1 and 2. Indeed, the solution to the Smoluchowski tree equation has three useful formulations.
Theorem 3 (Main result: solution).The solution g to the initial value problem (8) is given by, Further, with data g(0) = , g = g(t) solves the equation, Proof.Using (12), the first form of the solution can be re-written as, Now if we use Lemma 1, we observe, On the other hand using ( 18) again, we observe, which matches our expression for ∂ t g just above.The initial condition is naturally attained as t → 0. The second solution form shown follows from Lemma 2. We observe, The third solution form shown follows from the identification B∨ := G • B ∨ .For any n ∈ N 0 we have, giving the result.Finally, using ( 18) and ( 12) we observe, which also follows from the second solution form.

Remark 9 (Grafting and branching equilibrium).
An immediate consequence of Theorem 3 is that the solution form g given therein, satisfies, We can thus interpret the solution flow g in Theorem 3 as one for which, as time evolves, the processes of grafting and branching precisely match each other.Their actions are equivalent on such a flow, or exactly compatible.
Remark 10 (Exponential series solution).The solution form g in Theorem 3 is an exponential series of the form in Definition 4, with, for each τ ∈ T\∅: fτ = t |τ | .

Coagulation examples
Herein we explore some of the ramifications of the binary tree solution expansions in the last section, and in particular, make connections to cases with known explicit solutions.In this and the next section we are focused on the class of exponential tree expansions of the form (9) given in Definition 4. In our coagulation application we ally any such tree expansion with a positive function ξ representing the generalised transform of the initial data.Hence an exponential tree series of the form ( 9) in R T represents the following expansion in terms of the generalised transform function ξ, namely, Here, for any tree τ ∈ T\∅, the forms ξ(τ ) are to be interpreted as outlined in the introduction and at the beginning of Section 2, so that for example, and so forth.Suppose there exist positive constants c and Λ < 1 such that for all τ ∈ T\∅, we have, where here, • represents the absolute value-in this case of the real-valued coefficients fτ and functions ξ(τ ).
If condition (20) holds, the series ( 19) is convergent as, Here we used that for a fixed grade n, the sum, over all trees τ of grade n, of ω(τ ) is n!.We can prove this by induction.Assume the result to be true for k = 0, 1, . . ., n.
Recall the each of the trees at grade 'n+1' are constructed uniquely by the root grafting all the trees of grade 'k' to those of grade 'n − k' for k = 0, 1, . . ., n. Hence using Table 2: We list all the rooted planar binary trees of grade 4. For each tree τ ∈ T in the left column, we give its corresponding weight ω(τ ) in the second column, the total number of symmetries σ(τ ) it possesses (see Remark 13) in the third column and its word-coding in terms of levels in the final column.tree weight symm.levels the property of the weight character in Definition 3, we observe, which equals (n + 1)!, giving the result.Hence our algebra H = H(ξ, ⋆) in Section 2, denotes the class of convergent exponential tree series, convergent with respect to the data ξ, i.e. those exponential tree series for which the condition ( 20) is satisfied.
For the general frequency kernel case K = K(y, z), the solution to the coagulation equation ( 4) with the data ξ, represented by (14) in Theorem 3, is given by, We identify the exponential tree series coefficients in this case as ĝτ = t |τ | .This series is convergent provided t < 1 and ξ(τ ) c for some c > 0, for all τ ∈ T\∅.This brief, general analysis here, for suitable kernels for which ξ(τ ) c, thus only establishes existence of such a solution locally in time, indeed, for t < 1.
Presently we consider the three classical coagulation cases of the constant, additive and multiplicative frequency kernels.In preparation, as in Menon and Pego [77], we specialise the general transform H to the Bernstein transform.In this case the kernel of H is h(s, x) := 1−e −sx .Thus, the Bernstein transform represents a desingularised Laplace transform.In fact the Bernstein transform exists for any positive Radon measure ν t = ν t (dx) on (0, ∞) with scalar real parameter t, in the form, We assume the data ν 0 is a positive Radon measure on (0, ∞) and ξ is the Bernstein transform of ν 0 .Menon and Pego thus consider a weak formulation of (5) for which, where H(s, y, z) := 1 2 (1 − e −sy )(1 − e −sz ) and the double integral is over [0, ∞) 2 .For the precise details of the weak formulation setting, see Menon and Pego [77].For this special form of H, we can give an alternative characterisation for g ⋆ g as follows.Suppose we can expand K in the following separable form, for some constants c kℓ 0, where where for k ∈ N 0 , the moments M k := ∞ 0 x k g(x) dx.In this sum, when k or ℓ are zero, the added M 0 term should be taken to be zero.If the separable form expansion (24) is finite, then this definition for the non-associative product '⋆' could be computationally useful as we only have to compute a finite number of moments.Indeed, this is the case in the classical coagulation cases we now focus on.
Example 4 (Constant kernel).Consider the constant K = 1 frequency kernel case.The Bernstein transform g = g(s; t) of ν t satisfies, The product '⋆' is just the usual associative real product, with the − 1 2 factor.The solution form (16) in this case is given by g = (id − t B∨ ) −1 • ξ( ).The action of the operator B∨ in this associative case is simple.Indeed we observe that B∨ • ξ( ) = − 1 2 ξ 2 .Then the action of B∨ on this last form is to successively, additively replace each ξ factor on the right by − 1 2 ξ 2 , giving ( B∨ ) 2 • ξ( ) = 1 4 ξ 3 .Recall here that the grading operator is implicit in B∨ .We also deduce that ( B∨ ) 3 • ξ( ) = − 1 8 ξ 4 .and so forth.Hence we observe that, which is the solution in this constant kernel case.
Example 5 (Additive/multiplicative case).The additive K = y + z and multiplicative K = yz frequency kernel cases can be considered together for the following reason.The modified Bernstein transform g * of ν t is similar to the Bernstein transform except that the kernel h(s, x) := x(1 − e −sx ).Respectively in the additive and multiplicative cases, the Bernstein transform g of ν t and the modified Bernstein transform g * of ν t , satisfy, In the additive case we have normalised the constant first moment M 1 = 1.If we set h := e t g, then ∂h = e −t h∂h.If we set κ := 1 − e −t and h * (s; κ) := h(s, t), then h * satisfies, Hence it is sufficient to focus on the case, with g(s; 0) = ξ(s).It is well-known that for this case, corresponding to the multiplicative frequency kernel case for g * , that after a suitable normalisation, the gelation time occurs at t = 1.Herein we only consider the solution up to gelation, though a solution exists beyond that time, see item (ii) in Section 5. Note that for the rescaled additive case for h * , the gelation time corresponds to κ = 1 which translates back to t = ∞.See Menon and Pego [77] for more details.The solution to the inviscid Burgers equation ( 25) is uniquely obtained up to t = 1 via characteristics as, where the centre term is a compositional inverse.For the equation ( 25) with solution g and data ξ, we have, which is non-associative.We know the solution form in this context is any of the forms ( 14)-( 16) in Theorem 3, there is a close visual match to (16) in particular.By uniqueness the solution forms ( 26) and ( 14) are one and the same-one could simply iterate (25) as we did for (4).
Remark 11 (Explicit match).We can be more explicit about the solution match in Example 5 as follows.If we set h := (id − tξ) •(−1) , then the solution ( 26) can be expressed in the form , where the second factor is a reciprocal.Thus, since ∂ t g = (∂ t h)(∂ξ • h), we observe that, We can use this form to compute the Taylor series expansion for g in powers of t.Thus for example, we observe, and so forth.
In Doikou et al. [29], we show that the additive and multiplicative kernel cases can be described in terms of the Faà di Bruno composition algebra of analytic exponential power series.See Figueroa et al. [40] and Gessel [46] for more details.Indeed, in these cases, inspired by the result of Byrnes [15] and Byrnes and Jhemi [16] on Riccati partial differential equations, in Doikou et al. [29], we show that the solution g has a representation as a rank-one analytic Grassmannian flow.See Section 5 for more details on Grassmannian flows and the extension to the general frequency kernel case.
Remark 13 (Non-planar trees).With the exception of Examples 4 and 5, and Remark 11, for full generality, we assume throughout that the non-associative product '⋆' is non-commutative, and thus the rooted binary trees we consider are planar.In Smoluchowski's equation (1) the fields are scalar and thus the resulting product '⋆' is commutative.This means that many of the terms ξ(τ ), for trees τ of the same grade, are equivalent/collapse together.Indeed, we can represent any exponential tree series expansion on the set of non-planar binary trees as opposed to planar ones.Thus, for example, commutativity of '⋆' means that we do not distinguish between the trees, and , and the only two distinct rooted non-planar binary trees of grade 3 are, and .
We can think of the set of rooted non-planar binary trees as a set of equivalence classes of rooted planar binary trees, where we choose a single representative for the classes of planar binary trees that are equivalent in the non-planar setting.For example, in the case of the two grade 2 planar trees shown above, we might choose the left one to be the representative.We have given two possible representatives in the grade 3 case above.In general for example, we might choose the representative to be the planar binary tree with the lowest word-order code-see Section 4. We can determine which planar trees are equivalent in the non-planar setting by systematically, successively 'twisting' (or 'rotating') all the vertices present in a given planar tree, to see, if by doing so, they can be transformed into our given representative.For example, for the right-hand grade 2 planar above, twisting the top left vertex does not alter the planar tree-indeed this is true for any free vertex ends of any planar trees (i.e.vertices with no further higher attachments).However, if we twist/rotate the bottom vertex, we obtain the grade 2 tree on the left.In the case of the grade 3 planar tree on the right, we see that twisting any of the vertices it contains, does not alter it.However if we examine all the other trees of grade 3 in Table 1, we see that we can twist the vertices at the bottom as well as the first level up, to transform them into the tree of grade 3 shown on the left above.A careful examination of Table 2 reveals that there are three non-planar trees of grade 4, a further quick enumeration reveals there are six at grade 5, and so forth.This substantial reduction in the number of trees required to represent the exponential tree series solution (14) in Theorem 3, has important implications for numerical approximations, as we discuss in Section 4. In Tables 1 and 2 we give a column entitled 'symm.' which indicates the total number of symmetries σ(τ ) the planar trees τ ∈ T shown have.Thus for example, the trees of grade 2 in Table 1, can each be transformed into each other in the non-planar setting by vertex twisting.Thus for either of these trees τ , each generates 2 1 trees (including itself) in the non-planar setting and we set σ(τ ) = 1.The number of non-planar trees each generates is 2 σ(τ ) .Now consider the trees of grade 3 in Table 1.The middle tree in the list shown cannot be transformed into any other and thus σ(τ ) = 0 and the number of non-planar trees, including itself, it can generate is 2 σ(τ ) = 1.The four other trees τ of grade 3 shown can each be transformed into each other, and so for each of them, σ(τ ) = 2, and the number of non-planar trees each can generate by vertex twisting, including themselves, is 2 σ(τ ) = 4.The first tree τ in Table 2 has σ(τ ) = 3 as, including itself, it can generate 8 non-planar binary tree copies.And so forth.
Remark 14 (Positivity).An important property of solutions to Smoluchowski's equation ( 1) is that for positive initial data, the solutions should be positive thereafter while they exist.In the case of the constant, additive and multiplicative frequency kernel cases, Menon and Pego [77] (20) holds, there exists a convergent exponential tree series solution, we have not studied in detail when condition (20) does hold.For example, as mentioned in Example 5, it is well-known, see Menon and Pego [77], that in the multiplicative case a positive Radon measure-valued solution exists up to the gelation time which is governed by the initial second moment M 2 (0), which we can normalise to be one.This result can be found in McLeod [74,75,76].Comprehensive accounts of existing well-posedness results for the Smoluchowski equation can be found in da Costa [25] and Dubovskii [31].They include well-posedness for the case of a frequency kernel K(y, z) k 0 (1 + y + z) for some constant k 0 > 0, and a striking result by Carr and Da Costa [19] establishing instantaneous gelation if there exist constants α, β with β > α > 1 such that k 0 (y α + z β ) K(y, z) k 1 (yz) β for some positive constants k 0 and k 1 .Further general rigorous well-posedness results are given in Escobedo et al. [36] and Escobedo et al. [37] for coagulation-fragmentation models.Probabilistic approaches to solutions include that by Deaconu and Tanré [26] for the three classical kernels, and the characterisation of 'eternal' solutions, existing for t ∈ R, by Bertoin [10] for the additive kernel.Such well-posedness results are very kernel-specific.The frequency kernel K is an integral part of the product '⋆' and thus solution series expansion (21) through the terms ξ(τ ).A general investigation into well-posedness results that can be established through the analysis of the terms ξ(τ ) in ( 21) is a definitive worthwhile future endeavour.For example, we know from Example 5, that if g is the solution to the additive kernel case, then h * (s; t) = (1 − t) −1 g(s; − log(1 − t)) is the solution to the multiplicative kernel case.In this instance, in the multiplicative case, we observe there must be a natural renormalisation of the solution series expansion (21) generating a pre-factor (1 − t) −1 , thus explicitly exposing the gelation time.

Numerical simulation method
There are two parts to this section.In the first part we show how the solution flow ( 14)-( 16) can be formulated as a linearised flow.In the second part we establish a new practical numerical simulation method based on the binary tree expansion solution (14).
First, focusing on a linearised flow formulation, we derive two equivalent prescriptions.Though not strictly required here, it is useful to introduce the budding operator.We discuss it more generally in item (v) in Section 5.Here it is a re-interpretation of the branching operator.
Definition 7 (Budding operator).For any τ ′ ∈ T, we define the budding operator B τ ′ : R T → R T as the operator that acts on any tree τ ∈ T by successively, additively performing the following operation to each free end of τ .The operator B τ ′ extends the existing free end to a new branch on the right, and attaches, via a new branch, the tree τ ′ on the left.So a new vertex is created where the free end was, and the left branch of the new vertex has τ ′ attached, while the right branch is a new free end.Example 6.For example, for any τ ′ ∈ T we have, .
Importantly, we observe that B • ≡ B ∨ .Further, we can extend B τ ′ linearly so that Prescription 1 (Linearised flow).Suppose p ∈ R T , the linear operator Q : R T → R T , and the exponential tree series g ∈ R T satisfy, p(0) = , and Q(0) = id, and the linear system of equations, We can solve the first two equations giving, Naturally, we observe that the linearised flow Prescription 1, with (28) in hand, is just a re-writing of the solution flow ( 14)-( 16), however the interpretation of Prescription 1 is useful.Recall that ' ' in T represents the data function ξ = ξ(s) in the coagulation context.Note that Q = id − t B∨ is a linear endomorphism on R T , and thus the relation p = Qg in ( 27) is indeed a linear equation for g, whose solution is (16).To interpret this further we need to impose an ordering on the set of rooted binary trees T. Recall Tables 1 and 2 and focus on the final column labelled 'levels'.Roughly, we observe from the ordered listing of the trees shown in the tables, that we order them according to more foliage on the right-this being also interpreted at each vertex-through to more foliage on the left.The words in the final column 'levels' provide a coding for the trees represented, and, can either be read left to right, or from the root up.Each level refers to the vertical position of the vertices of the trees shown.
Procedure 3 (Tree order).An ordering for the set of rooted planar binary trees at each grade can be established as follows.Let us focus on the grade 4 trees in Table 2.
Interpreting left to right, consider the tree labelled 3212.Scanning from the left, the first vertex we encounter is at level 3, the next at level 2, then level 1 and then level 2 again.For the tree labelled 2341, the first vertex we encounter scanning from the left is at level 2. We then follow those connected vertices that head upwards, so we get 34, before, having exhausted that branch, we proceed to the next vertex at level 1.And so forth.The length of the word-code naturally corresponds to the grade of the tree.The interpretation of the 'levels' word-coding from the root up, proceeds as follows.Consider the tree labelled 3212.The 1 indicates the first bottom vertex.Each wordcode of any length corresponding to a tree with at least one vertex, has only a single 1, naturally.The two 2's straddling the '1' indicate that there are vertices attached to each of the ends of the bottom vertex.Each word-code of any length has either none, one or two 2's.The case of no '2' means that the tree in question is just the grade 1 tree.The case of only one '2' means that only one branch of the root vertex has a vertex attached.If the '2' occurs somewhere to the right of the '1' in the word-code, the vertex is attached to the right branch.If the '2' occurs somewhere to the left of the '1', the vertex is attached to the left branch.In the word-code 3212, the single 3-there can only be one further digit for a grade 4 tree-means that there is a single vertex at level 3.That the '3' is to the left of the left '2' means that the vertex is attached to the left branch of the left vertex at level 2. As another example consider the tree corresponding to the word-code 2431.The '1' denotes the root vertex.The single '2' to the left of the '1' indicates there is only one vertex attached to the left branch of the root vertex.The single '3' to the right of the '2' indicates there is a single vertex attached to the right branch of the single vertex at level 2. The single '4' to the left of the '3', indicates there is a single vertex at level 4 attached to the left branch of the vertex at level 3.And so forth.The set of rooted planar binary trees is then ordered by the numerical ordering, within the set of integers, of the word-codes.See the ordered lists in Tables 1 and 2, and Example 7 below.Remark 16.The word-code order of rooted planar binary trees above is closely related to their coding by permutations given in the original paper by Loday and Ronco [66].Also see Aguiar and Sottile [1], Arcis and Márquez [4] and Chatel and Pilaud [21].
Consider the class of exponential tree series f of the form ( 9), which we henceforth restrict ourselves to.We order the terms in any such series according to the tree ordering outlined in Procedure 3.With this in hand, we can represent any such series by the vector of the coefficients ĥτ := ω(τ ) fτ , associated with each tree τ ∈ T\∅, contained therein.In other words, we represent any such f by the vector, f = ( ĥ0 ; ĥ1 ; ĥ12 , ĥ21 ; ĥ123 , ĥ132 , ĥ212 , ĥ231 , ĥ321 ; ĥ1234 , where the sub-indices are the word-codings of the corresponding trees.We use semi-colons to separate coefficients associated with trees of the same grade.The basis elements are the corresponding trees τ with the reciprocal of the grade factorial factors '1/|τ |!'.Let us now consider the action of the branching operator B ∨ from Definition 5 in the context of such coefficient vectors f.Since B ∨ successively, additively attaches a branch ' ' to each free end of any tree τ ∈ T, its action on the word-coding of binary trees can be described as follows.For trees of lower grades, we observe: Here, in the second example, we separate the digits of the word-code by a '|', and the arrows '↑' indicate where a free branch of the corresponding tree is and thus where we can successively, additively add a next level integer.Consider, for example, the tree represented by the word-code 12.This contains a vertex at level 1 and attached to its right branch is a level 2 vertex.We can always attach a vertex ' ' to either end (left or right) of any trees as these are always free branches.For the word-code 12 let us scan free branches right-to-left.The form ↑ 2 ↑ indicates we can attach a vertex ' ' at level 3 to the right branch of the level 2 vertex, giving 123.And since the integer digit to the immediate left of the 2 is less in integer-order than the 2, there must be a free left-branch of the vertex at level 2 in 12, thus generating 132.Then the form ↑ 1 indicates, since the integer digit to the immediate right is greater than 1 in integer-order, that there is not a free branch to the right for this level 1 vertex.However, there is free branch on the left to which we can add a level 2 vertex, generating the tree 212.The interpretation of the form ↑ 2 ↑ | 1 ↑ is now straightforward.Note, starting with 12 + 21 in tree order, when we consider attaching vertices to free branches of 12 and 21, we can do so successively, additively from right-to-left.The order of the trees shown as the result of applying B ∨ to '12 + 21' above is then a sum of trees in tree-order.Though scanning right-to-left generates the trees at the next grade in tree order for a given word-code, this property does not generally propagate across different word-codes in tree order.Indeed, let us consider some more examples, as follows.We observe, For the example 32434512 we had above, we observe, Note that as we consider each digit in 32434512, free branches only occur when neighbouring digits are smaller in integer-order, as indicated by the ↑'s above.Now consider an exponential series represented by the vector (29).The action of the graded branching operator B∨ on f, i.e.B∨ ( f), is given by, B∨ ( ĥ0 ; ĥ1 ; ĥ12 , ĥ21 ; ĥ123 , ĥ132 , ĥ212 , ĥ231 , ĥ321 ; ĥ1234 , • • • ) T = (0; ĥ0 ; ĥ1 , ĥ1 ; ĥ12 , ĥ12 , ĥ12 + ĥ21 , ĥ21 , ĥ21 ; ĥ123 , • • • ) T .
Let V denote the vector space of exponential tree series (9) represented by vectors of the form (29).We can now give a second linearised flow prescription corresponding to that in Prescription 1, in the context of V, as follows.
Prescription 2 (Linearised flow: reprise).Suppose that the linear operators P, Q, G : V → V satisfy P(0) = Q(0) = id, and the linear system of equations, For this linearised flow, the first two equations imply, Then for sufficiently small t we have G = (id − t Bm ) −1 .
Relating this linearised flow to the linearised flow Prescription 1, we naturally observe that p = P( ) and the linearised flow solution g is given by g = G( ), where recall ' ' represents the data ξ( ) ∈ H.We can also see this from another perspective.The solution (14) for g is an exponential series of the form ( 9) with ĥτ = ω(τ ) t |τ | .Using the properties of the weight character ω from Definition 3, it is straightforward to show that the corresponding vector f of the form (29), satisfies ∂ t f = B m f.Note here, B m is the non-graded form of Bm .It is the same as Bm but a factor corresponding to the grade of the tree at that position in the vector is attached first, before applying Bm .
Second, we now explore new numerical simulation methods that can be constructed by approximating the general solution form (14) in Theorem 3, for general frequency kernels, as well as in particular, separable frequency kernels.Consider the exponential tree series solution (14), which is given in R T , transformed to the form (21) in H(ξ, ⋆), the algebra of convergent exponential tree series in the sense of condition (20).Suppose we truncate (21) to only include trees of grade, up to and including, N ∈ N. Denote this truncated series by, Let S N (ξ) := {ξ(τ ) : τ ∈ T and |τ | N } denote the set of tree-parametrised terms ξ(τ ) in the truncated expansion (30).We construct a new, simple numerical method, based on the truncated expansion (30), to integrate Smoluchowski's coagulation equation for a general frequency kernel K, over a global time interval [0, T ] for some T > 0. Herein, we suppose T precedes any gelation time.The basic aspects we need to consider are: 1. Time discretisation: If T is very small, we may not need to discretise at all.We might be able to simply evaluate the set S N (ξ) for the initial data ξ only, for sufficiently large N .Then the solution is given by g N , or in coagulation variables as H −1 g N .However, in general, we need to discretise the interval [0, T ] into, say, M equal subintervals (for simplicity for here) so that [0, T ] = ∪ M−1 m=0 [t m , t m+1 ], with t 0 = 0 and t M = T .Let dt := T /M denote the uniform time step.We need to compute the set S N anew at each time step t m , m ∈ {1, . . ., M − 1}-in addition to computing it at the initial time t 0 = 0. Since the computational effort for S N is large for large N , there is a natural trade-off between the sizes of M and N .Determining the optimal trade-off is very much of interest.

Evaluation of the expansion terms S N :
There are two components of this computational effort: (a) Non-planar trees: The '⋆' product associated with Smoluchowski's equation ( 1) for the scalar field g = g(x; t) is commutative and the exponential tree series solution ( 14) simplifies to an expansion over rooted nonplanar binary trees.See Remark 13 for the full details.This means that for example, the integrals represented by ξ ⋆ (ξ ⋆ ξ) and (ξ ⋆ ξ) ⋆ ξ are the same, and can be combined.At grade 3, there are only two distinct non-planar binary trees.Thus for the commutative case, which includes the original Smoluchowski equation for a general frequency kernel K, there is a very significant reduction in computation effort associated with computing the set S N at each time step t m , m ∈ {0, . . ., M − 1}.
(b)Generalised transform: We opted for the Bernstein transform in Section 3, but in particular circumstances, another linear transformation-such as the modified Bernstein cosine transform-might significantly improve the effort associated with computing S N .Indeed, in practice, the fast Fourier transform works extremely well.
With regards the time discretisation in item 1, the local error associated with the truncated approximation (30) which includes trees of grade N , across any time interval [t m , t m+1 ], is O (dt) N +1 .Since this local error accumulates linearly at leading order over the global time interval [0, T ], the global time discretisation error is O (dt) N .This of course assumes we have evaluated the S N sufficiently accurately at each time step.With regards item 2, let us consider part (a) first; part (b) is implicitly realised further below.Recall Remark 13 in Section 3. In Tables 1  and 2, we list all the planar binary trees τ up to, and including, those of grade 4, together with their associated weights ω(τ ) and symmetries σ(τ ).Recall that the symmetry σ(τ ) associated with a non-planar tree determines the number of non-planar copies, given by 2 σ(τ ) , of itself it generates.Let us consider the triples τ, ω(τ ), σ(τ ) of non-planar trees at each grade.Recall the word-coding representation we introduced in Procedure 3. At grade 1 there is only one non-planar tree ' ' equivalent to the triple (1, 1, 0); see Table 1.The triple of the non-planar tree of grade 2 is (12, 1, 1), while at grade 3 there are only Table 3: We list all the non-planar binary trees of grades 4, 5 and 6.For each word-coded tree in the left column, we give its weight ω(τ ) and symmetry σ(τ ) in the second and third columns, and how it is generated from lower grade non-planar trees in the final column.two non-planar trees with triples (123, 1, 2) and (212, 2, 0).Note, for non-planar trees, we choose the convention to use the planar tree with the lowest word-order coding as the non-planar representative tree.The triples for all the nonplanar trees of grades 4, 5 and 6 are given in Table 3.In the Table , we can compute the weight associated with each tree from weights of lower grade trees, using the formula in Definition 3. The symmetries σ(τ ) are easily computed in the manner outlined in Remark 13.Further, we note that ξ( ) = ξ(1) is given by ξ(1) = ξ(0) ⋆ ξ(0), where ξ(0) = ξ( ), or more abstractly 1 = 0 ⋆ 0. This is naturally just grafting.Thus we also have ξ(12) = ξ ⋆ ξ(1), or abstractly 12 = 0 ⋆ 1, and also we have, abstractly, 123 = 0 ⋆ 12 and 212 = 1 ⋆ 1.All the non-planar binary trees up to, and including, grade 3, are thus generated from lower grade non-planar trees in this way.In Table 3 in the final column, we indicate how all the non-planar binary trees of grades 4, 5 and 6 can be similarly generated.Let S np N denote the subset of S N consisting of non-planar binary trees, using the representation convention outlined above.The truncated binary tree expansion for non-planar binary trees has the form, Remark 17 (Algebra of word-codes).In general, for planar binary trees, we can define an algebra isomorphic to R T , based on the word-codes in Procedure 3. The product '⋆' is this algebra would be a ), with the first simple cases taking the form 0 ⋆ 0 = 1, 0 ⋆ 1 = 12, 1 ⋆ 0 = 21, 0 ⋆ 12 = 123, 0 ⋆ 21 = 132, 1 ⋆ 1 = 212, and so forth.
Let us now outline a practical approach to evaluating the set S np N at each time step t m , m ∈ {0, 1, . . ., M − 1}.Given the details just above, we essentially need to compute, H(s, y, z)K(y, z) g(y)f (z) dy dz, (32) for suitable pairs of functions ξ and η.In this product, ξ = Hg and η = Hf , where H is the generalised transform with kernel h, as outlined in the Introduction and discussed in Section 3, and H(s, y, z) := 1 2 (h(s, y +z)−h(s, y)−h(s, z)).Let us immediately restrict ourselves to the case when the frequency kernel is separable, in the the sense that, where k is a given function k : [0, ∞) → [0, ∞).One reason for this choice is to compare our results to those in Filbert and Laurençot [41].We outline at the end of this section our strategy for more general frequency kernels.
The Bernstein transform we considered in Section 3 is a natural choice for the separable frequency kernel case as , and in this case the double integral in (32) decouples to become the real product of two integrals.We wish to retain this property and take advantage of the computational efficiency of the fast Fourier transform.For convenience, we define the simple linear function ℓ : [0, ∞) → [0, ∞) to be, Given the domain of integration in (32) involves the positive quadrant, a natural choice is to suppose h generates the modified Bernstein cosine transform, i.e. we take h(x, s) = 2x cos(2πsx).With this choice, If we subsitute this form for H and the separable frequency kernel form (33) into the product (32), it decomposes into the sum of four terms, each one of which is the real product of two factors.The factors are either cosine, Bernstein cosine, or sine transforms of kg, kf , ℓkg or ℓkf .These arguments can be extended as even or odd functions on the whole real line and thus the transforms in the factors can be expressed in terms of Fourier transforms of the extended arguments.Hence we can in principle use the fast Fourier transform to compute all the factors described.Note, for example, for the first term in binary tree expansion (31), given initial data such as g 0 = exp(−x)/x, then in the first time step the expression ξ⋆ξ can be computed analytically.
In practice, using the fast Fourier transform to evaluate (32) for separable frequency kernels is even more straightforward.We use the following form for the Fourier transform f of the function f and its inverse, Let us also outline the discrete Fourier transform (DFT).Understanding how this works is crucial to the numerical simulation method we propose.See Press et al. [81] for more details.Suppose we are given a function f on a finite set of equispaced nodes x n := νh, ν = 0, 1, . . ., n−1.Here, h is the discrete 'spatial' scale.Set f ν := f (x ν ).Then the Fourier transform f is approximated by, where The fast Fourier transform (FFT) computes the sum on the right, dropping the h prefactor.The discrete inverse Fourier transform is given by, For the FFT, the set of n frequencies k are reordered.For convenience we denote the Fourier transform of f by F (f ).Now suppose h(x, s) := x exp(2πisx).Then in this case, If we substitute this form for H into (32), using the separable frequency kernel form (33), we find that, where F 0 (f ) denotes the Fourier transform evaluated at s = 0.This result is not correct unless we extend the arguments kg, ℓkg, kf and ℓkf as even or odd functions on the whole real line.Given one of its original applications, signal processing, the DFT/FFT is ideally suited to the situation we now find ourselves in.In practice we are given initial data g 0 = g 0 (x) for x ∈ [0, ∞), and it is natural to sample this at the nodes x ν := νh for ν = 0, 1, . . ., n − 1, where h := L/n.Here [0, L], for L > 0, is a sufficiently large truncation of the coagulation domain [0, ∞).Often we are given data that is singular at the origin, for example g 0 (x) = exp(−x)/x, and in this instance we take x ν = (ν + 1)h, for ν = 0, 1, . . ., n − 1. Indeed this is our modus operandi herein.However, note that the DFT in (34), naturally does not reference the actual nodal positions x ν , but utilises the n function nodal values f ν and that the nodes are equispaced.It generates a discrete transform at the frequencies k, indicated above.When computing the inverse DFT via (35), we only use the n transform values f k .We are thus free to chose the equispaced nodes we wish to evaluate our original function at, thus determining its locale on the real line.Naturally we choose the nodal points so the locale is in [0, ∞), in particular we choose x ν = (ν + 1)h, for ν = 0, 1, . . ., n − 1.Thus in practice, having sampled the functions g = g(x) and f = f (x) at the nodal points x ν , as well as k and ℓ, we can replace the Fourier transforms F of the arguments kg, ℓkg, kf and ℓkf in (37) by the DFT/FFT.It remains to describe the overall algorithm, implement it for some poignant examples, and demonstrate how our method can be improved and generalised to wider classes of frequency kernels K. To describe the overall algorithm, it suffices to show how the approximate solution is evaluated over the first time step, as subsequent time steps simply repeat the process.Let us consider the case when we use the truncated non-planar binary tree expansion (31) for N = 3, which generates a third order integrator in time.Given initial data g 0 , we thus need to compute, To compute ξ = ξ(0), we simply compute the FFT of g 0 , having sampled g 0 on the n nodes x ν = (ν + 1)h, for ν = 0, 1, . . ., n − 1.For the first term on the right above, we actually don't need to do this.At the end of each time step we evaluate the approxomate solution g = g(x, t m ) in the coagulation space, as explained presently.Then to compute ξ(1) = ξ ⋆ ξ we use (37) with the FFT.We simply need to compute h times the FFT of kg 0 and ℓkg 0 .The term F 0 (kg 0 ) is h times the first element in the FFT vector of kg 0 .We then compute the inverse fast Fourier transform (iFFT) of ξ(1), call this g(1), as we need this in the next step.To compute ξ(12) = ξ ⋆ ξ(1), we use (37) with g given by g 0 , and f given by g (1).We then compute the iFFT of ξ (12), nominate this as g (12).Finally, to compute ξ(123) = ξ ⋆ ξ(12), we use (37) with g given by g 0 , and f given by g (12).To compute ξ(212) = ξ(1)⋆ξ(1) we use (37) with g and f both given by g (1).We then compute the iFFT of 4 • ξ(123) + 2 • ξ(212).We can then evaluate g 3 , in coagulation space, by considering the linear combination of the corresponding terms in (38).The approximation g 3 is then the initial data corrsponding to g 0 for the next time step.Using Table 3 the procedure for computing higher time-order approximations is now straightforward.For the third order case we considered here, the total number of FFTs required, including the iFFTs, is ten.And of course the effort of each such FFT operation is proportional to n log n.For the separable frequency kernel case, this compares favourably with the numerical method of Filbet and Laurençot [41].We implement six numerical approximations based on the truncated non-planar binary tree expansion (31), as just described, for N ∈ {1, 2, 3, 4, 5, 6}, for different numbers of time steps M on a global time interval [0, T ].In all cases the truncated 'spatial' domain is [0, L] with L = 100.We used n = 2 20 modes to evaulate the required FFTs in all cases, except where stated otherwise.As in some of the examples outlined in Filbet and Laurençot [41], we set, We consider three separate cases, when the parameter λ is set to be λ = 2, λ = 3/2 and λ = 2/3.The first case, λ = 2, corresponds to the solvable multiplicative frequency kernel case.Indeed, for the initial data g 0 = exp(−x)/x there is a well-known closed form solution given by, for t 1, where I 1 is the modified Bessel function of the first kind.This solution extends beyond the gelation time t = 1; see (ii) in Section 5.In the top panel in Figure 1, we give a log-log plot of global error versus the number of steps M used to compute the approximation at time T = 0.5, for all six cases N ∈ {1, 2, 3, 4, 5, 6}.The L 2 error was computed by comparing the approximations to the exact solution (39).All the methods have the order anticipated.We used the same number of nodes n = 2 20  for all the methods.We observe that the higher order methods error curves flatten off around 10 −3 , this just an artifact of the 'spatial' discretisation error, determined by n, starts to exceed the time step error.
The second case, λ = 3/2, also corresponds to a gelation case.For the same initial data g 0 (x) = exp(−x)/x, we computed the solution up to time t = 0.9 for all six integrators.In the middle panel in Figure 1, we give a log-log plot of global error versus the number of steps M .The L 2 error was computed by comparing the approximations computed with n = 2 20 modes versus the sixth order approximation computed with the smallest step size and n = 2 21 'spatial' modes.This log-log error plot looks very much like that for the λ = 2 case, and again the flattenning off of the error curves for the higher order methods is just an artifact of the 'spatial' discretisation error exceeding the time step error.We give a log-log plot of the solution in the left panel in Figure 2.
The third case, λ = 2/3, is non-gelling.For the initial data g 0 (x) = exp(−x), in the bottom panel in Figure 1, we give a log-log plot of global error versus the number of steps M , computed to the time T = 1.5.The L 2 error was computed using the same procedure we outlined for the case λ = 3/2 just above.Again the log-log error plot has the same characteristics as the previous two cases above.
In the right panel in Figure 2, we give a log-log plot of the solution computed to the time T = 15 in this caseusing M = 2 7 time steps.This compares well with the corresponding plot in Filbet and Laurençot [41], namely in Figure 11, in the lower left panel therein.Overall, our numerical method based on the truncated non-planar binary tree expansion (31) and FFT computations appears robust and efficient.
The method we have outlined can be improved and generalised in several ways as follows.
1. Non-uniform fast Fourier transform: The FFT computations above relied on the nodal points x ν being equispaced across the domain [0, L], where L = 100.However, for all the cases we considered, for large x, the initial data decays like 'exp(−x)'.It would therefore seem sensible to distribute the nodal points to this 'distribution', so that more nodal points are concentrated around smaller values of x.The non-uniform fast Fourier transform (NUFFT) does indeed exist and the effort required is still proportional to n log n, though with a larger constant.In its implementation in our context here, for an exponential distributed set of nodes x ν , we would compute the NUFFT on the same set of frequencies as above.For a different initial data, we may wish to choose a differently distributed set of nodes x ν .The inverse NUFFT can also be computed, though this requires some additional steps.For more details, see Kircheis and Potts [62].In particular, once implemented, we could in principle, evaluate the solution for very small nodal points, eg. 10 −40 , as in Filbet and Laurençot [41].
2. General additive kernels K(x, y) = k(x) + k(y): The method we have outlined is straightforwardly adapted to this class of frequency kernels, which were also considered by Filbet and Laurençot [41].

General kernels:
The method we have outlined is easily adapted to any frequency kernels K = K(y, z) consisting of a finite linear combination of the forms k i (y)k j (z), for i, j = 1, . . ., l, and such that k i (z)k j (y) = k i (y)k j (z).However, when l is large, or if K cannot be expressed as such a finite linear combination, then we might be able to proceed as follows.We set h(x, s) = x exp(2πisx) as above.Then given g and f , computing the product (32), with H given by (36), amounts to computing the two-dimensional Fourier transform, or two-dimensional FFT, along the diagonal set of wavenumbers.In principle, we require the one-dimensional iFFT of the resulting expression for ξ ⋆ η, in order to implement our method, as outlined after (38).

Discussion
There are many connections and extensions to the material herein.A brief summary of examples is as follows.
(i) Grassmannian flow: In Doikou et al. [29] we outlined in detail how the constant, additive and multiplicative frequency kernels can be considered as Grassmannian  flows.In principle, we can extend these cases to the general frequency kernel case, i.e. the solution flow ( 14)-( 16) can formally be formulated as a Grassmannian flow as follows.Recall the vector space V from Section 4. We now construct the Grassmannian Gr(V ⊕ V ⊥ ; V), the set of all subspaces W of V ⊕ V ⊥ such that: (1) the orthogonal projection pr : W → V is a Fredholm operator, indeed a Hilbert-Schmidt perturbation of the identity; (2) the orthogonal projection pr : W → V ⊥ is a Hilbert-Schmidt operator.See Pressley and Segal [82].Herein we assume V ⊥ is isomorphic and isometric to V. Any such subspace W has a representation of the form, W = (Q; P), where Q is a Fredholm operator on V which is a Hilbert-Schmidt perturbation of the identity, and P is a Hilbert-Schmidt operator on V. Let V 0 denote the canonical subspace with the representation, V 0 = (id; O), where O is the infinite matrix of zeros.The projections pr : W → V 0 and pr : W → V ⊥ 0 respectively give, W := (Q; O) and W ⊥ := (O; P).These projections are possible provided det 2 Q = 0.The subspace spanned by the columns of W coincides with V 0 , and indeed, Q −1 transforms span{W } to V 0 .Under this transformation, the representation W for W becomes, (id; G), where G = P Q −1 .Any subspace that can be projected onto V 0 can be represented in this way and vice-versa.The Hilbert-Schmidt operators G parametrise all subspaces W that can be projected in this way.If det 2 Q = 0, a different coordinate patch must be chosen.For more details, see Pressley and Segal [82] or Doikou et al. [30], as well as Beck et al. [5,6].A possible Grassmannian flow prescription for our flow herein, corresponding to the linearised flow Prescription 1, is the linearised flow Prescription (2) for the linear operators P, Q and G.This approach needs further investigation.In particular, an important aspect of Grassmannian flows in a given coordinate patch, is that flow singularities correspond to a poor representative patch, where det 2 Q = 0, but the solution flow can still be represented and continued in a different coordinate patch.This is relevant to the following item on gelation.
(ii) Computing beyond gelation: Quoting from van Roessel and Shirvani [93]: "The phenomenon whereby conservation of mass breaks down in finite time is known as gelation and is physically interpreted as being caused by the appearance of an infinite 'gel' or 'superparticle'...".In turn, also see Ernst et al. [35].Herein, for frequency kernels for which gelation occurs, we have focused on computing the solution up to the time of gelation.During this interval of time, the quantity M t := ∞ 0 x g(x; t) dx is conserved by the flow.However, post-gelation, quoting from Normand and Zambotti [79], the superparticles "do not count in the computation of the mass so from the gelation time on, M t starts to decrease".Indeed the solution can be computed, see for example, Leyvraz and Tschudi [65], van Roessel and Shirvani [93] and Normand and Zambotti [79].Using the binary tree expansion to compute the solution g, analytically and/or numerically, post-gelation, is very much of interest.Indeed a Grassmannian flow context would seem to be a natural one to encapsulate the bahaviour before, during and after such a phase transition, via different coordinate patch representations.
(iii) Kingman coalescent and Galton-Watson processes: Smoluchowski coagulation models are naturally represented in terms of coalescent stochastic processes, in particular Kingman coalescent processes or Galton-Watson branching processes (reversed in time).See, for example, Aldous [2], Iyer et al. [57], Harris et al. [51] and Johnston et al. [58].Etheridge [38,Ch. 2] demonstrates the following insightful connection for the solution of a classical quadratically semilinear parabolic partial differential equation.The solution can be expanded deterministically by an iterative procedure akin to that we performed for the coagulation equation (5).The terms in the solution expansion are indexed by rooted planar binary trees.Etheridge then shows that, at a given time, the terms at each grade are given by the expectation across branches of a given branching process-by analogy with McKean's [73] solution of the Kolmogorov-Petrovski-Piskunov equation via branching Brownian motion.It is thus natural to try to establish such a connection between each term of the rooted planar binary tree expansion (14) and the branches of the underlying Galton-Watson process, reversed in time.
(iv) Multiple mergers: A natural extension of the Smoluchowski coagulation model is to coalescent phenomena involving multiple mergers.In principle the solution in such cases can be similarly expanded as an exponential trees series analogous to (14) indexed by planar trees, except now the set of indexing planar trees would be more general and could include the complete collection of n-ary planar trees reflecting the class of merger coalesence included.See, for example, Iyer et al. [57] for more details on Smoluchowski models with multiple coalescence.Also see Doikou et al. [29,Sec. 3].
(v) Decorated trees: We considered a non-commutative, non-associative algebra with one generator and the nonassociative product '⋆'.Naturally in general, we can construct such algebras with more than one generator.In this instance we can use the set of decorated binary trees to represent the monomials in such an algebra-or in the multiple merger case, just decorated trees.In the decorated binary tree context, we could extend the action of the budding operator B as follows.For example, suppose there are three generators ξ, η and ζ in H, respectively represented by '•', '•' and '⋄' in T. Then we could define, .
Also recall Remark 8, we could also consider extending this to include B τ , for any τ ∈ T, which might model multiple stage reactions depending on the tree attached.Note how B is left-budding-for each "free" branch we extend the existing bud to the right and attach the new one to the left.Naturally there is a right-budding version as well.
(vi) Species: The coagulation equation ( 5), whose solution ( 14) is expressed as an exponential series in rooted planar binary trees, can be expressed as an equation in the species of rooted planar binary trees.The extensions mentioned in (iii) above to more general coalescent processes are, in principle, examples of further species of structures.See Bergeron et al. [9] for more details on species.
(vii) Multi-component coagulation: There are multicomponent generalisations of Smoluchowski's coagulation model (1), in particular in the context of atmospheric sciences, where "clusters can be formed by different types of particles"; see Throm [92].Adapting our binary tree expansion approach to such generalisations appears to be straightforward, and again, very much of interest.
(viii) Hopf algebras of trees: We mentioned in the introduction the depth and wide ranging applications of algebras of planar trees.If the algebra of planar trees is endowed with the grafting product, one can define a compatabile co-product, and an antipode, and thus establish a Hopf algebra of such trees; see Loday and Ronco [66].We mention, due to their more direct relevance here, the use of rooted planar trees in Lie group methods and backward error analysis in Munthe-Kaas and Wright [78] and Lundervold and Munthe-Kaas [67,68].Indeed, there is more than one Hopf algebra of planar trees, see Calaque et al. [17].One co-product that would be useful in our analysis would be that which, for a given planar binary tree τ , it additively enumerates all the possible pairs of trees that when root grafted together, generate τ .Thus for example, for a given general tree series expansion with terms fτ •τ , such a coproduct helps keep track of the origin of coefficients that are generated by root grafting two such series together.This means that, in principle, there is an equivalent formulation of ( 14) as a co-product expansion.And, in principle, this could lead to a more abstract formulation of (14) in terms of Hopf algebra endomorphisms, by analogy with such expansions in Ebrahimi-Fard et al. [33] and Ebrahimi-Fard et al. [34].Lastly, we also mention here the work by Ishida [56] on the Lie algebra of rooted planar trees, Chapoton [20] on exponential-like series and Gerritzen [45] on non-associative exponential series.
(ix) Free pre-Lie algebra: Al Kaabi [3] considers the free pre-Lie algebra structure associated with rooted planar trees.The construction of numerical algorithms in the free pre-Lie algebra context is very much of interest.
(x) Branching Brownian motion and colloids: Smoluchowski diffusion models incorporate the spatial Brownian motion of clusters.They have applications in theory of colloids and sedimentation.In such models, the density g = g(x, ζ; t) is recorded at position ζ.A branching Brownian motion can be considered as a "Gaussian process indexed by the leaves of a Galton-Watson process"-Bovier [11].Viewed backwards in time we observe a diffusive coalescent, and Smoluchowski diffusion models can be interpreted in this light; see Harris et al. [51].Can we extend the connection in (iii) above, here between the deterministic interative solution expansion and the underlying diffusive coalescent process, in this case?See Etheridge [38] and Dynkin [32] for more details of the diffusive branching case, and Berestycki and Berestycki [7] and Berestycki et al. [8] for more details of the diffusive coalescent case.

Acknowledgement
The author would like to thank the referees for their positive reports and suggestion to implement the numerical scheme.This directly led to the material in the latter half of Section 4 and helped to significantly improve the original manuscript.The author is also very grateful to one of the referees for bringing [65] to his attention.

Funding and conflicts or competing interests
SJAM was supported by an EPSRC Mathematical Sciences Small Grant EP/X018784/1.There are no conflicts of interests or competing interests.

Data availability statement
No data was used in this work.All the Matlab codes are provided in the electronic supplementary material.
any a, b ∈ R. As in Example 3, we set B := G•B.Naturally we also have B• ≡ B∨ .We now present a linearised flow prescription in the context of R T .

Example 7 .
Two further illustrative examples of the word-coding are as follows,
[39,blish positive Radon measure-valued solutions provided the initial measure is positive.This achieved by establishing that the Bernstein transform of the measure solution, or its derivative, is completely monotone in the sense of Feller[39, XIII.4].Menon and Pego's results cover Examples 4 and 5 above.We do not pursue positivity in this sense further here, however a valuable future investigation would be to determine or characterise the conditions under which the general solution (14) for g possesses complete monotonicity or a similar property guaranteeing positive Radon measure-valued solutions.