Sparse‑Lagrangian PDF Modelling of Silica Synthesis from Silane Jets in Vitiated Co‑flows with Varying Inflow Conditions

This paper presents a comparison of experimental and numerical results for a series of turbulent reacting jets where silica nanoparticles are formed and grow due to surface growth and agglomeration. We use large-eddy simulation coupled with a multiple mapping conditioning approach for the solution of the transport equation for the joint probability density function of scalar composition and particulate size distribution. The model considers inception based on finite-rate chemistry, volumetric surface growth and agglomeration. The sub-models adopted for these particulate processes are the standard ones used by the com-munity. Validation follows the “paradigm shift” approach where elastic light scattering signals (that depend on particulate number and size), OH- and SiO-LIF signals are computed from the simulation results and compared with “raw signals” from laser diagnostics. The sensitivity towards variable boundary conditions such as co-flow temperature, Reynolds number and precursor doping of the jet is investigated. Agreement between simulation and experiments is very good for a reference case which is used to calibrate the signals. While keeping the model parameters constant, the sensitivity of the particulate size distribution on co-flow temperature is predicted satisfactorily upstream although quantitative differences with the data exist downstream for the lowest coflow temperature case that is considered. When the precursor concentration is varied, the model predicts the correct direction of the change in signal but notable qualitative and quantitative differences with the data are observed. In particular, the measured signals show a highly non-linear variation while the predictions exhibit a square dependence on precursor doping at best. So, while the results for the reference case appear to be very good, shortcomings in the standard submodels are revealed through variation of the boundary conditions. This demonstrates the importance of testing complex nanoparticle synthesis models on a flame series to ensure that the physical trends are correctly accounted for.


Introduction
A large variety of nano-sized particulates are produced by combustion processes, some as unwanted by-products in the energy sector such as soot, others as commodities such as catalysts or stabilising agents (Iler 1979). In all cases, the final product characteristics are significantly influenced by precursor reactions, interparticulate interactions and the gas-phase conditions (composition, temperature, fluid dynamics) under which they are formed. The numerical simulation of all these processes and the accurate prediction of the particulate size distribution (PSD) are vital if we are to understand and control flame synthesis. The evolution of the PSD is usually described by the population balance equation (PBE) that can be solved by Direct Monte Carlo methods (Kruis et al. 2000), moment methods (Pratsinis 1988) or discretised sectional methods (Friedlander 2000;Marchisio and Fox 2013). When the PBE is coupled with the Navier-Stokes equation, modelling is required for the multiscale, non-linear interactions of the particlulate dynamics and the turbulent flow. In the context of large eddy simulations (LES), Loeffler et al. (2011) and Wang and Garrick (2006) investigated the formation and growth of titanium dioxide nanoparticulates, but neglected particulate-turbulence interactions. Neuber et al. (2017) used the sectional method and showed that the error associated with omission of sub-grid turbulence interactions can be large for the highly non-linear nucleation and surface growth processes. Similar conclusions were drawn by Pesmazoglou et al. (2017) and Sewerin and Rigopoulos (2018) who modelled particulate aggregation and soot formation and growth, respectively.
In Neuber et al. (2017Neuber et al. ( , 2019a we presented the implementation of a stochastic Monte-Carlo method for the spatial and temporal evolution of the joint filtered density function (FDF) of the gaseous composition and number density of the particulates. The number density is discretised into sections representing particulates of different sizes. To avoid ambiguities, we strictly distinguish in this paper between "stochastic particles", which are computational elements for solving the FDF equation in the Lagrangian sense, and "particulates", which are the physical nanoparticulates produced during the flame synthesis process. Conventional stochastic FDF methods require stochastic particle numbers up to 20-50 per LES cell (You et al. 2017), and even though optimised methods-such as in-situ adaptive tabulation (Pope 1997)-exist for the calculation of reaction kinetics, the simulation of flame synthesis processes remains computationally expensive. When using a sectional method, the computing time for the particulate synthesis may even exceed the computing time for the chemical kinetics of the gas phase. For such problems it is therefore beneficial to use a so-called sparse-Lagrangian FDF approach, which requires far fewer stochastic particles than traditional intensive methods (Cleary and Klimenko 2011). The key enabler of the sparse FDF approach is the use of a novel mixing model called multiple mapping conditioning (MMC) (Klimenko and Pope 2003) which helps to correctly emulate molecular and turbulent diffusion despite the relatively large physical distances between the particles to be mixed. The MMC model has been validated widely for different turbulent flames, e.g. see (Cleary and Klimenko 2011;Vo et al. 2017a;Galindo-Lopez et al. 2018;Neuber et al. 2019b). Recent publications have shown that the method can also be applied to particulate nucleation and growth processes (Neuber et al. 2017;Vo et al. 2017a) where interactions with turbulence are important. We demonstrated the capability of the combined PBE-MMC-LES method to model particulate inception, surface growth and agglomeration of non-spherical agglomerates in a turbulent, reactive jet flow (Neuber et al. 2019a). That study included, however, only one specific setup. MMC could capture some features of the flame synthesis process but notable differences between computations and experiments were observed.
Complete model evaluation demands the computation of a series of cases with varying boundary conditions to establish the sensitivity of the model and its input parameters. The comparison with a series will also ensure the model's capacity to capture trends and thus to have incorporated all the important physics of the particulate synthesis process. Here, we apply the model to the flame synthesis of silica nanoparticulates from a silane doped nitrogen jet issuing into a vitiated hot co-flow and investigate the sensitivity towards variations in precursor loading and gas phase temperatures and assess the model's performance against experimental OH-LIF, SiO-LIF and elastic light scattering (ELS) signals.

Methodology
PBE-MMC-LES is a hybrid Euler/Lagrange approach where the mass and momentum transport equations are solved by an Eulerian LES solver, here employing standard closures for the sub-grid terms. The gas-phase composition and particulate size distribution are solved using a Lagrangian Monte Carlo approach whose details are provided now. For the solution of the particulate size distribution we use the nodal form of the sectional approach (Prakash et al. 2003) which discretises the particulate volume space into a finite number of sections. If the particulate matter is transported by a fluid of gas-phase composition Y , temperature T and velocity field u j (x i , t) , the nodal form of the population balance equation takes the form where n k is the number density of particulates in section k, D k is the diffusion coefficient of section k, and G(v, Y, T) is the volumetric growth term that can be modelled as a function of particulate volume, v, and Y . Changes of the PSD due to inception and agglomeration are included via the source term ̇s k . The rate at which primary particulates are formed is determined by the production rate of the precursor species, which in turn is given by the gas phase chemistry described by the underlying reaction mechanism. Incipient particulates that are newly formed by chemical reaction are added to the first section of the discretised particulate size distribution.
The volumetric surface growth is determined by two processes, each dominating in different size ranges. For small particulates the surface growth is determined by a collision process whereas for large agglomerates the surface growth is driven by diffusion processes. For the former, the growth rate is based on the free-molecular collision kernel. The collision diameter of the species which is depositing on the agglomerates surface is assumed to be the diameter of the molecule (Shekar et al. 2012). In our sectional approach we consider fractal-like structures with their morphology given by a power law (Friedlander 2000). The collision diameter of the agglomerate is given by d c,k = d p,0 N 1∕D f k , with N k being the number of primary particulates in the agglomerate, d p,0 being the diameter of the primary particulate and D f being the fractal dimension which is set to D f = 1.8 as suggested in Refs. Shekar et al. (2012); Schaefer and Hurd (1990). Hence, the volumetric growth rate is given by the collision rate of the depositing species with the agglomerate. For the larger particulates we apply a diffusion-limited growth method (Witten and Sander 1981) (1) because there the volumetric surface growth is mainly determined by the diffusion rate of the depositing species towards the agglomerates' surface. We use the harmonic mean of the two formulations which serves as a blending function between the two regimes. For agglomeration the Fuchs interpolation expression between the free molecular and continuum regimes was proposed by Seinfeld (1986). For arbitrarily shaped agglomerates (Kruis et al. 1993) proposed to replace the spherical particulate diameter in the Fuchs interpolation expression by the collision diameter of the agglomerate and this was also done here. The implemented discretised representation of the PSD requires a size-splitting operator for the agglomeration term which ensures particulate number and mass conservation by distributing the formed particulates into two adjacent sections (Loeffler et al. 2011). This procedure introduces a slight broadening of the PSD but this effect is small compared to the broadening due to growth and agglomeration. Due to the relatively low temperatures in regions where particulates are formed we assume that the agglomerates are not subject to any sintering processes. This assumption is justified as agglomerate samples in the exhaust do not show significant signs of neck formation between the primary particulates (see also Fig. 2 and discussion in Sect. 3).
We apply a Lagrangian scheme to approximate the evolution of Eq. (1) as part of the joint FDF, F , of the composition field = (Y 1 , … , Y s , h, Z, n 1 , … , n e )which includes the gas phase composition vector Y , the total enthalpy, h, the mixture fraction, Z, and the discretised particulate number density field n . The FDF method has the major advantage that the chemical source terms as well as the rates for nucleation, volumetric surface growth and agglomeration appear in closed form in its governing equations. Thus, no closures are required for these terms to incorporate unknown sub-grid effects due to turbulence. As we use a Monte-Carlo technique a fractional step approach is applied to describe the spatial dispersion, gas phase reaction, mixing and aerosol dynamics (Pope 1985). The time evolution is then given by the equivalent stochastic differential equations (Cleary and Klimenko 2011) for particle transport and the change of the composition field in time is governed by Cleary and Klimenko (2011) where is a Wiener process and standard notation is used for diffusion, density, velocity, space and time. The molecular diffusion coefficient, D, is set equal for all species and particulate sizes except in our discussion conducted in Appendix 3 and D t denotes the turbulent diffusivity. S accounts for source terms including chemical reactions, particulate inception, surface growth and agglomeration and M is a mixing operator which emulates the sub-filter scalar dissipation. A suitable reaction mechanism provides chemical reaction rates and thus determines the particulate inception (cf. Sect. 4). Models for volumetric surface growth and agglomeration have been discussed above.
The mixing operator appears in unclosed form and requires modelling. A mixing model like the interaction by exchange with the mean (IEM) (Villermaux and Devillion 1972), Curl's mixing model (Curl 1963), the modified Curl mixing model (Janicka et al. 1979) or the Euclidean minimum spanning tree (EMST) model (Subramaniam and Pope 1998) can be used to close the equation on mixing operator level. Subramaniam and Pope (1998) specified requirements for mixing models, such as linearity, independence and localness in composition space. A mixing model which provides such requirements is the multiple mapping conditioning (MMC) model (Klimenko and Pope 2003) as it enforces localisation in composition space indirectly by localisation in an independent reference space. MMC requires the additional solution of such a reference field. This reference field is the (Eulerian) LES solution of mixture fraction, f , and the mixing of the stochastic particles is then conditioned on f . The localisation in reference space allows to increase the spatial distance of the stochastic particles as long as localness in composition space is maintained. This has led to the development of the so-called sparse-Lagrangian MMC method, where there are fewer stochastic particles than LES grid cells. In contrast to conventional mixing models sparse particle methods use up to three orders of magnitude fewer particles resulting in significant computational savings. Here, we use Curl's mixing model where particles are pairwise mixed in combination with the MMC conditioning. The particle pairs are selected such that they are local in reference mixture fraction space and the mixing extent is controlled by a mixing time scale, L . For a detailed discussion on the modelling of the mixing time scale the reader is referred to Vo et al. (2017a, b) and Neuber et al. (2019b).

Experimental Configuration
The general experimental configuration is based on the Cabra burner (Cabra et al. (2002)) where a central turbulent jet (here nitrogen) issues into a hot vitiated co-flow of premixed hydrogen-air combustion products. The central jet has a diameter of D = 4.57 mm and the outer diameter of the vitiated co-flow is 210 mm . The cold jet bulk velocity is U j = 33.2 m/s resulting in a jet Reynolds number of Re D = 10,000 . The co-flow velocity is U c = 1 m/s based on cold conditions (i.e. no combustion of co-flow mixture). The nozzle exit plane extends 70 mm above the perforated plate of the co-flow. The experimental campaign consists of measured signals from parameter variations where the central jet is doped with different concentrations of silane (0, 300, 1000, 2500, 2700, 2900 and 3100 ppm). Also, the co-flow temperatures (1300, 1500 and 1800K ) and the jet Reynolds numbers (5000, 10,000 and 15,000) were varied. The case with no silane doping (0 ppm) is used as baseline for the OH-LIF signal and served for the validation of the flow and mixing field predictions (Neuber et al. 2019a). Preliminary measurements demonstrated that dopings below 2000 ppm do not lead to any detectable light scattering signal from the particulate matter while for a range between 2500 and 3100 ppm a large variation of the signal is observed (cf. Fig. 1). There is a clear temperature dependence of the signal over the temperature range between 1300 to 1800 K and the dependence is less pronounced for lower loadings such that here, the analysis with respect to silane loading considers cases with silane doping larger or equal 2500 ppm only. The temperature sensitivity studies focus on the case with 3100 ppm and the Reynolds number dependence is based on 2500 ppm doping due to restrictions of the mass flow rate controller in the high Reynolds number case.
Silica particulate synthesis is investigated experimentally by two laser diagnostic techniques, planar laser-induced fluorescence (PLIF) and elastic light scattering (ELS). The investigation is supplemented with extractive particulate sampling. OH-PLIF and SiO-LIF are used to validate the chemical kinetic model whereas ELS is used to validate both the particulate formation model and the temperature field. Model validation follows the "paradigm shift" approach described by Connelly et al. (2009), which involves computing 'synthetic signals' that are compared with experimentally-acquired signals as a strategy to avoid incurring assumptions needed by the experimentalist for conversion of the signal to physically meaningful quantities.
For ELS, a four-head Q-switched frequency-doubled Nd:YAG laser is used in conjuction with sheet-forming optics to illuminate a 23 mm × 400 μm region across the burner centerline. The ELS signal is collected using a CCD camera with a 105 mm f/16 lens perpendicular to the laser beam. Using three of the four available heads, the shot energy is approximately 3 × 450 mJ . Measurements of different regions of the jet are performed by vertical translation of the burner. Data post-processing includes background subtraction, shot energy correction, beam profile correction, spatial calibration and image de-warping. Besides corrections that stem from the particulars of the experimental setup such as background subtraction and laser profile and energy corrections, no physical interpretation of the signal is attempted, as the measurement process is already simulated within the numerical model.
The fourth head of the Nd:YAG laser is used to pump a frequency-doubled Rhodamine 6G dye laser, tuned to a wavelength of 283.6 nm to excite the Q 1 (8) transition ( = 1 ← 0) of the A 2 Σ ← X 2 Π electronic system of the OH radical. This transition has been chosen so that the LIF signal is only weakly sensitive to temperature. A CCD camera with intensifier (multi-alkali, P43) and f/2.8 UV lens is placed at 90 • to collect the fluorescence signals from the ( = 0 ← 0) and ( = 1 ← 1) branches of OH around 309-315nm . Laser energy and sheet thickness have been adjusted to ensure operation within the linear LIF regime (Seitzman and Hanson 1993). The area probed is 23 mm x 700 μm and the shot energy 1.15 mJ. Again, background subtraction, shot energy correction, beam profile correction and beam extinction along the direction of travel of the beam are applied to the experimental signal.
For the SiO-LIF experiments, the fourth head from the Nd:YAG laser was converted to third harmonic emission (355 nm) and the dye laser converted to run on a Coumarin dye mixture. 100 mJ of 355 nm radiation was used to pump the dye laser, generating around 15 mJ pulse energy, tuneable around 460-470 nm. The beam was then frequency doubled in a BBO crystal, producing up to 1.5 mJ of narrowband radiation which could be varied in the range 230-235 nm. The laser was tuned to the Q 11 (J = 32) A 1 Π ← X 1 Σ + (0, 0) absorption line in SiO (235.087 nm) and the fluorescence signal was detected and corrected in the same way as the OH experiments (extinction corrections were not made, since SiO was not found in the co-flow).
Additionally, particulate samples were obtained in two ways, using a thermophoretic sampling device (TPS) similar to that described in Dobbins and Megaridis (1987) and using a Dekati PM10 impactor. The TPS samples were collected at specific locations in the jet, whilst the impactor samples were collected far downstream ( ∼ 2 m), in the extraction hood positioned over the burner. The TPS device uses a double-action pneumatic cylinder to rapidly push a 3 mm perforated carbon TEM grid into the jet and remove it after a pre-set residence time of 50 ms. The grid is held vertically during sampling in order to minimize disruption to the flow. Both types of samples are analysed using a Jeol 2100+ transmission electron microscope (TEM) at an acceleration potential of 200 kV without additional preparation and images are exported to Image-J for analysis. Figure 2 shows typical TEM images of particulates extracted from within the flame and from the exhaust. TPS samples at z = 70 mm had a dwell time in the jet of 250 and 50 ms in the 1000 ppm and 3100 ppm cases, respectively. The downstream samples, collected at z = 2 m , resulted from drawing the sampled gas through the finest impactor filter for 5 min at 3100 ppm and 20 min at 1000 ppm. Only few and small particulates can be seen on the samples for a silane loading of 1000 ppm indicating the presence of nucleation and surface growth but only moderate rates of agglomeration at this silane doping concentration. The picture for the case with 3100 ppm is not that clear. Much larger silica structures can be observed for samples taken within the flame, but clear agglomerate structures are absent. This "coating" of the sample probe may be due to condensation of the silica on the surface of the sample probe forming a surface film. When sampling within the exhaust two meters above the burner, large particulates of fractal shape have deposited on the sample probe indicating that the agglomeration process is the predominant growth mechanism and that sintering effects are irrelevant at conditions investigated here, see Fig. 2d. These images do certainly not represent the agglomerates' PSD, however, they clearly demonstrate the moderate particulate nucleation in the case with 1000 ppm and significant agglomeration yielding very large cluster for loadings with 3100 ppm.

Numerical Configuration
The hybrid Euler/Lagrange approach has been implemented into a code package called mmcFoam (Galindo-Lopez et al. 2018) which is based on OpenFOAM-5.0. The computational domain extends 25D in the axial direction and 11D in the radial direction. A priori investigations of three different grids with 0.5, 1.5 and 4 million cells revealed that results for the latter two cases are similar. The 1.5 million cell grid is used for all results presented in this paper. The mesh is refined near the nozzle, which is resolved by 45 cells along the jet diameter giving smallest cell sizes of 50 μm in the center and the shear layer of the jet. The turbulent sub-grid viscosity is modelled by the -model (Nicoud et al. 2011) and a model constant of C = 1.5 is applied. Pipe flow simulations were conducted inside the nozzle to provide realistic turbulent inflow velocity boundary conditions and zero-gradient outflow conditions are used at all other boundaries. Second-order central difference and TVD schemes are used for discretization of momentum and species transport, respectively. For the Lagrangian scheme, 230 000 stochastic particles have been used to compute the subgrid distribution of composition and particulate number density, which corresponds to one particle for every 6.5 LES cells. All MMC modelling parameters are standard as defined in earlier publications (Neuber et al. 2017;Cleary and Klimenko 2011;Vo et al. 2017b;Neuber et al. 2019b).
A finite rate chemistry model for the precursor chemistry is applied to model silane combustion (Suh et al. 2001;Suh et al. 2002). We use the mechanism provided by Suh et al. (2001) including 63 species and 264 reactions at their reported rates for atmospheric pressure. The mechanism includes clustering of silicon oxides in the gas phase that leads to nucleation of the first incipient particulates once the molecules are large enough. Suh et al. (2001) selected ( SiO n ) 11 to represent the first solid particulate but acknowledged this choice to be arbitrary as it depended on a trade-off between increased detail for gas-phase precursors and increasing uncertainties in the reaction kinetics for larger gas-phase molecules. Assuming spherical symmetry of incipient particulates with a solid matter density of p = 2196 kg∕m 3 gives a diameter of d p,0 = 0.98 nm for all primary particulates. A constant primary particulate size appears restrictive and Fig. 2d indeed demonstrates that particulates of various sizes are present. However, a sectional method including a primary particulate size distribution (i.e. an effectively two-dimensional representation of the particulate characteristics parameterized by the radius of gyration, R g , and d p,0 ) is not yet feasible within a PBE-MMC-LES approach and is not attempted here. Therefore, the stochastic particles carry the information of 63 species, enthalpy, mixture fraction and a suitable number of sections of the PSD. This "suitable" number warrants some more discussion: The number should be as small as possible to avoid unnecessary computational overhead but Neuber et al. (2017) showed that an unsuitably coarse discretisation can lead to excessive numerical diffusion and to an overprediction of the number of larger particulates. Since only nucleation and condensation were considered therein, the sensitivity study should be repeated here to ensure independence of results in the presence of agglomeration. Figure 3 shows the time-averaged PSD on the centerline of the jet with a 3100 ppm silane loading and a co-flow temperature of 1500 K at two different downstream positions. The number of sections has been varied from 30 to 120. It can be seen that all PSDs look similar and have a bimodal character, which indicates that numerical diffusion is small and that agglomeration effects are already dominating the distribution's spread across particulate size space at z∕D = 10 . For a higher number of sections the solutions converge. We can clearly observe that for 30 and 60 sections numerical diffusion leads to deviations from the converged solution for the PSD. The results for 90 and 120 sections are very similar. Since calculations with 120 sections increase computational requirements by 25% but differences in predictions are small, we use 90 sections for the approximation of the PSD for all further results presented in this study.
In the next section, we compare measured with computed signals following the "paradigm shift" approach which was used for a variety of analyses in turbulent combustion (Torniainen et al. 1998;Floyd and Kempf 2011) including particulate nucleation and growth (Connelly et al. 2009). This strategy avoids incurring assumptions needed to convert experimental signals to physical meaningful quantities such as particulate number density. The procedures to calculate the ELS and LIF signals are described in Appendices 1 and 2 , respectively.

Results
The PBE-MMC-LES calculations are performed for different silane loadings of the jet stream, different jet Reynolds numbers and different co-flow temperatures. Statistics have been collected for a minimum of 15 flow through times after a statistically The general numerical setup of this configuration including the predicted flow and mixing fields were validated in Neuber et al. (2019a). There it was shown that good agreement of measured and computed mean and standard deviation of the elastic light scattering (ELS) signal can be achieved. For the case with no silane doping, this indicated good predictions of the temperature (and therefore the mixing) field and that the most important flow phenomena can accurately be reproduced. The agreement for the ELS signal also indicated good predictions for silica nucleation and growth in the case with 3100 ppm silane doping of the jet. However, as measured signals require calibration and calibration has been realized by matching the peak centerline values, a certain degree of agreement is expected to be observed. For a more critical analysis, a series of test cases needs to be assessed and the model's capability to predict trends and to capture the underlying physics needs to be investigated. We now show further validation for the case with the highest silane loading of 3100 ppm (Sect. 5.1) and take this case as reference. Sections 5.2 (variation of silane loading) and 5.3 (variation of coflow temperature and jet Reynolds number) then provide a critical assessment of the model's capability to predict changes in process and boundary conditions. The variation of model parameters allows the identification of reasons for remaining discrepancies between measurements and computations. We complement our analysis by sensitivity studies towards variation of model parameters. However, these variations are ad hoc, they shall serve as indication of sensitivities and identify needs for further model improvements, but do not necessarily provide new suitable modelling constants. Results are therefore deferred to Appendices 3 and 4 .

Evaluation of the Reference Case
In this section we present results of our reference case with a silane doping of 3100 ppm and a co-flow temperature of 1500 K at a Reynolds number of Re D = 10,000 . Figure 4 compares the ELS signals from the experiments with those from the simulations. The signals from the LES are calculated as described in Appendix 1. We have calibrated the experimental signal such that it is unity in the co-flow. The predicted signal has two contributions. The signal originating from the gas phase is normalised to match the mean co-flow values, and the signal associated with the particulate matter is normalised such that the peak values of experiment and simulation agree. Note that this procedure is followed for the reference case only. For all other cases reported below the ELS signals are normalised by the same constants and no further adjustments have been applied. Overall, the simulations agree well with the measurments when assessing the mean of the ELS signal. Most notably the positions of the signals' peaks are very similar.
There are, however, a few pronounced differences that warrant some discussion: 1. The spatial extent of high particulate matter concentrations is somewhat too wide. This might be caused by the unity Lewis number assumption used here for the particle composition space including the size sections. The particulates can have very large Schmidt-numbers that are up to several orders of magnitude larger than the gaseous species' Schmidt numbers, and their molecular diffusion will tend to zero. Larger particulates emitting a higher signal tend to follow the streamlines while equal diffusivity assumptions enhance lateral diffusion. Vo et al. (2017a) suggested a modification of the Lagrangian mixing time scale for the particulate matter and this effect is discussed in Appendix 3. In addition, the random walk model (represented by the last RHS term in Eq. 2) induces particle dispersion and thus diffusion independent of the species specific molecular diffusion coefficient. This "enhanced" diffusion of heavy particles is more pronounced in hot regions where molecular diffusion coefficients are large and of the same order of magnitude as the turbulent contribution. In the present setup, the co-flow is hot which contributes to the over-prediction of particulate dispersion. The mean drift model developed by McDermott and Pope (2007) may extenuate this artificial diffusion but its implementation in the context of a sparse particle method is unclear to date and beyond the focus of this paper. 2. A comparison of Fig. 4c, d reveals qualitative differences between experiments and computations. The experiments show a clear maximum along the centerline while simulations place the highest fluctuations at the edges of the jet right into the shear layers. Again, neglecting differential diffusion in the model may be the cause of these differences. The heavy particulates would cluster in the inner jet and the larger concentrations there would also lead to larger gradients and thus variances within the jet's center. Indeed, modification of the mixing time scale such that there is differential dissipation of variances of the number densities of small and large particulates enhances variances along the centreline as shown in Appendix 3. 3. The computations predict an early onset of the mean ELS signal and corresponding fluctuations in the shear layers between jet and co-flow (cf. Fig. 4b, d). This is not observed in the experiments. The appearance of the simulated ELS signal within the upstream shear layer can be explained by the (modelled) dynamics within this jet that dominate silica particulate formation and growth. Figure 5 shows the inception rate, volumetric surface growth, primary particulate number density and agglomerate number density on a plane at centerline position. The inception rate is very high in the shear layer close to the jet exit, where the hot surrounding co-flow mixes with the jet and mixing leads to oxidation of silane and particulate inception. Along the centerline the inception rate peaks at z∕D = 10 and decreases significantly further downstream, which is in line with the observation that the precursor species is completely consumed at z∕D = 20 (not shown). The volumetric surface growth differs qualitatively from the inception rate. It is very large in the shear layer close to the jet exit but rapidly decreases further downstream. This is consistent with the observation that strong surface growth occur in regions where many intermediates are present which then deposit on the particulates' surfaces. Consistent with the inception rate, the number of primary particulates increases first in the turbulent shear layer and reaches its highest value at about z∕D = 15 (cf. Fig. 5c). Further downstream, the primary particulate number density decreases again due to ceasing primary particulate inception and due to particulate dispersion. Comparison of agglomerate number density with the number density of primary particulates highlights the influence of agglomeration which is in line with the analysis of the TEM pictures in Sect. 3. In regions where the primary particulate number density is highest, the number density of the agglomerates decreases very quickly due to particulate collision events.
The early inception of particulates in the shear layer is likely to be related to the gas phase chemistry including silane oxidation and formation of gas phase precursors. The validation of the gas phase kinetics includes a comparison of OH and SiO signals. OH can be understood as an indication of the position of the shear layer and it is linked to the kinetics of the silane oxidation process that is typically relatively fast and intimately coupled with the underlying hydrogen-oxygen kinetics. Radial profiles of OH mean values and standard deviation are shown in Fig. 6a, where experimental and predicted signals are normalised with their respective co-flow values at z∕D = 3 . It is seen that the positions of silane decomposition and OH production are well detected but peak values of the mean OH-LIF signal is over-predicted by 27%. This is a good agreement given the uncertainties that can be associated with the kinetics of silane oxidation. In line with the trends observed for the mean values the standard deviations are also overestimated. Note that the experimental standard deviations in the co-flow do not tend to zero. This is mainly due to shot noise and to a lesser extent to readout noise and temperature fluctuations in the hot co-flow. Overall, we achieve a satisfactory agreement between the measured and predicted position of the reaction zone and between actual and predicted silane decomposition representing the fast reactions within this process. OH is not, however, a good indicator for the silica formation process. Instead, SiO may serve as a representative species as it is the key intermediate for the formation of gaseous SiO 2 and one key species for (i) surface growth and (ii) growth of the gas phase precursors leading to the first particulates. Figure 6b shows excellent agreement between measured and predicted SiO-LIF signals. The location of the peak values and the variation with downstream distance for mean and standard deviations agree very well. Deviations at larger radii can be associated with the experimental read noise and shot noise that have a more substantial impact for SiO than for OH due to the relatively small SiO signal intensity.
With respect to the chemical kinetics and possible turbulence-chemistry interactions we may conclude that MMC-LES captures silane oxidation and formation of SiO (and possibly of the first silicon dioxide molecules) in the gas phase accurately. High SiO concentrations are correctly predicted in the shear layer close to the jet exit. The apparent agreement of key species in the gas phase and our prediction of the early onset of particulate nucleation (with the latter not being observed in the experiments) suggests some inaccuracies related to the modelling of cluster growth (the formation of ( SiO 2 ) m and (SiO) m with m ∈ [2, 10] ) in the gas-phase leading to particulate inception. The relevant kinetics were developed for standard pressures at T= 773 K (Suh et al. 2002) and extrapolation of the respective rate constants to the range of temperatures relevant in the present configuration may introduce some unavoidable inaccuracy. Direct validation of the cluster growth process is difficult as direct measurements of SiO and SiO 2 clusters do not exist and are not easy to perform-in particular for SiO 2 due to its photophysical properties. We can, however, assess the sensitivity of predictions on precursor kinetics and this is discussed in more detail at the end of Sect. 5.2. Figure 7 depicts the means of the ELS signals along the centerline for cases with 2500, 2700, 2900 and 3100 ppm silane loading. For the reference case with 3100 ppm the predicted ELS signal increases with the same rate as the signal from the measurements, indicating that particulate number and size are well predicted. Also growth and agglomeration are likely to be suitably modelled. The measured signal reaches its maximum value at z∕D = 14 , while the predicted peak value is very close at z∕D = 15 . The simulations predict the correct trends, i.e. lower silane doping leads to lower agglomerate numbers and therefore lower ELS signals. However, the sensitivities towards changes in silane doping are moderate while experiments feature much larger differences in signal strength. Only a small signal is detected for 2500 ppm and the increase in signal strength is strongly nonlinear with increasing silane concentrations. The strong dependence is unexpected but consistent with the TEM images shown in Sect. 3: for the low silane loadings only few small agglomerates were captured while higher loadings led to significant increases in particulate production and agglomeration. In contrast, the relatively low sensitivity observed in the simulations is not unexpected and can be explained by the models used for nucleation, growth and agglomeration. Nucleation and growth are linearly dependent on silane loading, and this linear trend is observed in the predictions where the predicted peak ELS signal is (to a good approximation) proportional to the silane loading. Non-linearities that may explain the experimental trends appear in the current models in three expressions only:

Sensitivity Study for Varying Silane Loadings
(1) the agglomeration reveals a square dependence on particulate number, (2) the computed ELS signal increases non-linearly with primary particulate size, agglomerate size and fractal dimension, and (3) chemical kinetics are non-linear. Appendices 1 and 5 give Exp., 2500 ppm Exp., 2700 ppm Exp., 2900 ppm Exp., 3100 ppm Sim., 2500 ppm Sim., 2700 ppm Sim., 2900 ppm Sim., 3100 ppm details on the computation of the ELS signal, highlight the signals sensitity to parameters and demonstrate that neither agglomeration nor ELS signal computation can reasonably explain the large discrepancies. The analysis in App. 5 hints at possible changes in D f that might be caused by the different seedings and lead to the different signal strenghts. These changes are, however, not clearly identifiable nor quantifiable for seedings between 2500 and 3100 ppm as the TEM images are snapshots and do not provide adequate statistics of the agglomerates' morphologies. This leaves-next to differential diffusion effects (cf. Sect. 5.1)-chemical kinetics as the primary suspect. Precursor growth in the gas-phase follows the sequence for m ∈ [1, 9] and n ∈ [1, 2] . The cluster ( SiO n ) 11 represents the first incipient particulates (Suh et al. 2001). An adjustment of the rate constants of these cluster growth mechanisms will lead to non-linear changes in particulate properties. In addition, Fig. 4 indicates particulate nucleation to be predicted too early and a reduction in rate constants seems justifiable. Appendix 4 demonstrates that adjusted rate constants could yield a behaviour similar to the experimental trends observed in Fig. 7. Simple perfectly stirred reactor computations show that the ELS signal depends strongly non-linearly on silane doping at early times. A sufficient delay in particulate nucleation would then yield similar non-linearities at z∕D = 15 where the maxima of the measured ELS signals are observed. Appendix 4 also shows additional MMC-LES computations of the Cabra burner with reduced rate parameters. The results demonstrate that (i) nucleation and growth is shifted downstream leading to absence of a particulate signal in the shear layer and continuous growth along the centerline, but (ii) shifts the peak ELS signal significantly further downstream and overall agreement does not improve. It seems that particulate dispersion and lateral species diffusion is not strong enough to suppress a further increase of the ELS signal on the centerline beyond z∕D = 15 . We may conclude that this sensitivity study corroborates a strong influence of the kinetics on the predictions but simple scaling will not suffice and key to success will be detailed and validated chemical schemes that ensure the correct growth of precursor species in the gas-phase. Figure 8a shows the ELS signal's sensitivity towards co-flow temperature. The ELS signals for the upstream positions are in excellent agreement with experiments indicating an accurate numerical treatment of the varying boundary conditions. The ELS signal further downstream is correlated to nucleation and growth of the particulates indicating that the silane conversion is slower for the lowest co-flow temperature leading to lower particulate number densities. Measured and predicted ELS signals have a similar peak value for the two higher co-flow temperature cases and we may conclude that temperatures above 1500 K are high enough for a fast and complete conversion of the precursor. Again, the larger changes observed in the experiments cannot be fully reproduced by the simulations when decreasing the co-flow temperature from 1500 to 1300 K, but we note that the model correctly captures the highest ELS signal for the middle temperature of 1500 K.

Sensitivity Study for Varying Co-flow Temperature and Reynolds Number
In Fig. 8b the sensitivity of the ELS signal towards the jet Reynolds number is shown. Due to restrictions of the mass flow rate controller for the high Reynolds number case, the silane doping of the central jet is set to 2500 ppm for all three cases. Results for the high and medium Reynolds number cases are consistent with the observations made above: no significant signal could be measured, but MMC-LES predicts presence of some particulate matter. Notably different is the low Reynolds number case. The much slower convective velocity provides sufficient time for particulate nucleation to set-in. At the same time, (turbulent) dispersion is very low as the flow is quasi-laminar. The reduced dispersion leads to rather high particulate concentrations, fast agglomeration and a strongly non-linear dependence of the ELS signal on Reynolds number variations. MMC-LES cannot capture this re-laminarization of the flow and-similar to above-unity Lewis number assumptions lead to enhanced dispersion and an almost linear dependence of the peak ELS signal is predicted.

Conclusions
Sparse-Lagrangian PBE-MMC-LES calculations including detailed precursor chemistry, particulate inception, volumetric surface growth and agglomeration have been conducted for a series of silane doped nitrogen jets in a hot co-flow. The reference case with high silane loading is very well predicted and trends for variations in silane loading, jet Reynolds number and co-flow temperature can be captured. However, the measured sensitivities are rather strong and non-linear while the model (based on the assumptions inherent in the submodels) gives an almost linear or at most quadratic dependence on silane concentration. Further comparison with LIF measurements of key gas-phase species representing the fuel conversion process and intermediates for silica formation demonstrates the method's capability to represent gas-phase conversion and its interactions with turbulence. However, sensitivity studies indicate that the gas-phase precursor (cluster formation) chemistry can cause these strong non-linearities, but simple scaling of precursor (cluster) growth does not lead to success. Also, differential diffusion effects are likely to have significant effect on the specific locations in the flame where particulates nucleate and grow. The current work shall be understood as an attempt to identify the specific submodels' sensitivities, to provide guidelines for the need of future development in the area of modelling particulate flame synthesis and to highlight the need for proper model validation. The validation with one set of measurements only may not identify shortcomings of models or sub-models as the "right" choice of modelling constants may often allow for a good match between simulations and "the" one target experiment. As has been established within the combustion community, a series of measurements with parameter variation is needed to provide an unbiased assessment of modelling approaches. The-at certain conditions-modest agreement between measurements and simulations should not lead to dismissal of some of the sub-models used here, it rather highlights the complexity of the particulate formation process. Nanoparticulate flame synthesis in general and silica particulate formation in particular is similar to soot inception and growth. The quality of predictions presented here is comparable to the quality of predictions of soot in turbulent flames found in the literature. Also, the difficulties with respect to modelling soot formation are commonly associated with the complex kinetics leading to the incipient soot particulates and with differential diffusion. This is consistent with the current paper that identifies the precursor chemistry and differential mixing of heavy species and particulates as key modelling issues that can explain the strongly non-linear dependence of particulate properties on silane doping, and more work is required in these areas if simulations are to be truly predictive for particulate synthesis processes.
where T is the temperature, x i is the species mole fraction and i is the ELS (Rayleigh scattering) cross section of species i at 532 nm . Normalised ELS cross sections have been determined by Fuest et al. (2012). Since composition and temperature are fully defined in the simulation the molecular signal can be predicted within the model. The second contribution of the ELS signal originates from light scattered by agglomerates and is calculated according to the Rayleigh-Debye-Gans theory of light scattering by fractal agglomerates (RDG-FA) (Sorensen 2001). In the case that a sectional method with N s sections is used the contribution of each class of agglomerates is added to give the total signal where the C agg sca is the scattering cross section of the fractal agglomerate and n k is the number density within the section k. The former is a function of the agglomerate size (Sorensen 2001) and needs to be weighted with the computed PSD to be comparable to the measured signal.
For loose, spherical primary particulates several times smaller than the wavelength, the scattering cross section varies with the sixth power of the primary particulate diameter (Rayleigh's approximation). When large fractal aggregates form, the scattering cross section also depends on the radius of gyration. Following the conventions of Link et al. (2011) we can formulate a dependency of the scattering cross section on fractal dimension and size for three different size regimes: where N is the number of primary particulates within the aggregate, the scattering wave vector q = 4 sin( ) and k = 2 ∕ are constant throughout the experiment since = 532 nm and = ∕2 . The value of the constant C is discussed in Sorensen (2001) and is set to C = 0.77 . As discussed in Sect. 2 the fractal dimension and the primary particulate diameter are set to D f = 1.8 and d p 0 = 0.98 nm , respectively.
The normalisation factors in Eq. (4) are determined based on the results of our reference case with a silane doping of 3100 ppm and a co-flow temperature of T c = 1500 K at Re D = 10,000 . The normalisation factor for the molecular ELS signal is determined to a = 6.52374 ⋅ 10 −4 such that a value of unity is obtained for the co-flow. The ELS signal from agglomerates is normalised with b = 1.01 ⋅ 10 −4 to match the same maximum value as the experimental reference case.
Two implications of Eq. (7) are noted: 1. C agg sca scales with d 6 p,0 . We may be tempted to assume that the modelling uncertainty with respect to the initial primary particulate size and variations in growth rates may have a very significant non-linear effect on the ELS signal. Indeed, particulate diameters measured at 70 mm downstream can easily vary by a factor of 2. It is hypothesized, however, that particulate number density scales with d −3 p,0 and therefore S agg ELS computed by Eq. (6) is much more likely to show a square dependence on the size of the primary particulates. Also note that TEM pictures indicate a very similar mean for the primary particulate size for the entire seedings range. Mean diameters are around 11 nm for 300 ppm and around 14 nm for 3100 ppm, and it may be concluded that primary particulate size is not (primarily nor solely) responsible for the strong non-linear effect observed for the different silane loadings. These observations may not hold in case of significant sintering after collision. However, this is not observed in the current application where primary particulates can easily be identified in the samples collected in the exhaust (cf. Fig. 2d).

Appendix 3: The Influence of a Mixing Time Scale Modification for the Transported Number Densities of the Discretised Particulate Size Distribution
Here, a modification of the Lagrangian mixing time scale is discussed. If unity Lewis number assumptions are invoked, all transported scalars are mixed with a mixing extent calculated on the basis of the anisotropic mixing time scale model proposed by Vo et al. (2017b).
Here we show and discuss results obtained with a modification of the Lagrangian mixing time scale for the particulate number densities only, as proposed by Vo et al. (2017a). They showed that for transported particulate matter a modification of the Lagrangian mixing time scale can achieve a significant improvement of the conditional variance. The standard mixing time scale is scaled with the diffusivity of the individual sections in order to reduce the sub-grid mixing of the aerosol, L,k = D k ∕D L , where L,k is the mixing time scale of section k. This emulates reduced mixing rates for large agglomerates. The reference simulation uses the same Lagrangian mixing time scale for all transported scalars including the particulate number densities and has been presented in Sect. 5.1. There, the predicted mean of the ELS signal agrees well with measurements but qualitative differences are observed for the standard deviations (see Fig. 4). It is noted, however, that the location of the peak mean ELS signal moves 5 diameters downstream. Figure 9 now shows simulated mean and standard deviation of the ELS signal for adjusted mixing time scales mimicking differential diffusion of the particulate matter. All remaining (gas-phase) quantities are mixed with the standard model. Improvements can be observed as (i) particulate nucleation and growth are moved further downstream leading to a reduced ELS signal in the shear layers upstream and (ii) a much better qualitative agreement for the standard deviation is observed (compare Fig. 9b with Fig. 4c).
A more quantative comparison is provided by Fig. 10 where radial profiles of the ELS signal are shown for different downstream positions. Here, the presence of the ELS signal in the shear layer for the case with equal diffusivities is pronounced and predictions with a modified mixing time scale show improved predictions. Also, the peak values of the standard deviations are much better predicted with the reduced mixing frequencies for the aerosol matter, especially for centerline positions where the peak value are perfectly matched. However, overall agreement does not really improve: the predicted standard deviation profiles are much broader than the measured data and more importantly, the centerline dependence of the mean signal is not capture well. We should point out here that the modifications suggested by Vo et al. (2017a) do not provide a fully consistent differential diffusion model and more sophisticated models as derived e.g. by McDermott and Pope (2007) and Dialameh et al. (2014) should be extended to the present application for accurate predictions of differential diffusion effect. Here, however, we limit ourselves to identify the sensitivities of results on molecular diffusion processes and point out that differential diffusion may need to be properly modelled if agreement with measurements is to be achieved at all positions in the flow.

Appendix 4: Discussion of the Influence of Reaction Rate Parameters
In Sect. 5.1 we have shown results of the reference case with a silane doping of 3100 ppm and a co-flow temperature of T = 1500 K . It has been discussed that the computed ELS signal becomes large in the shear layer, whereas this characteristic cannot be observed for the experimental data. There, the ELS signal increases at centerline position in streamwise direction. One reason could be an incorrect precursor chemistry where chemical reaction rates are overpredicted and inception is too fast. To investigate the influence of the precursor chemistry we first investigate the time evolution of particulate inception and growth. Figure 11 (left) shows perfectly stirred reactor calculations with different initial silane concentrations. After a transient, the ELS signals are almost linearly dependent on the silane loading. However, during an initial stage (Fig. 11 (right)) the ELS signals are strongly non-linearly dependent and the sensitivity resembles the sensitivity observed in Fig. 7. The rate coefficient for the clustering reactions given in Suh et al. (2002) are estimates and some adjustments seem justified. A reduction in rates may delay particulate nucleation sufficiently to (i) reduce particulate nucleation in the shear layer and (ii) "freeze" particulate formation sufficiently early due to dispersion and mixing prior to completion of the entire conversion process from silane to silica particulates. Figure 12 shows the mean ELS signal of simulations where the pre-exponential factors of all clustering reactions are divided by 2 and 10, respectively. As expected, the particulate inception is shifted further downstream with negligible particulate formation in the upstream shear layers which is consistent with experiments. Thus, the precursor chemistry model-and uncertainties related with it-can have a strong impact on the particulate evolution processes. We note, however, that the peak ELS signals are shifted further downstream (with the maximum for the ̇i∕10 case likely to be located outside the computational domain). The overall agreement with measurements is therefore not improved, and we may postulate a need for improved precursor kinetics if simulations of silica production processes are to be predictive.

Appendix 5: Investigation of the Influence of the Particulates' Fractal Dimension on the ELS Signal
The frequency of agglomerate-agglomerate collision is given by the collision kernel. Kernels require the specification of collisional cross-sections that can be parameterized by the agglomerates' radii of gyration and fractal dimensions. A typical fractal dimension for silica flame synthesis reported in the literature is D f = 1.8 (Shekar et al. 2012;Schaefer and Hurd 1990) and is assumed to be constant in our sectional approach. However, the fractal dimension may change during the growth process due to collision of different agglomerates (Inci et al. 2017) or due to surface growth and neck formation. To asses the influence of different morphologies Fig. 13 shows the mean ELS signal along the centerline for four MMC-LES simulation of the reference case (as defined in Sect. 5.1). The only change relates to using different-but still constant-values of the fractal dimension. The signal originating from the gas molecules is not affected and signals are nearly equal down to z/ D=5 when the first particulates form on the centerline. Then, the ELS signal increases with different rates and the signal originating from the solid phase starts to dominate. It can be observed that the maximum signal strength shifts further downstream for larger fractal dimensions. This indicates-as expected-a decreased agglomeration rate for more compact agglomerates. It needs to be noted, however, that the respective maximum values do not reveal a clear dependence on D f as the maximum value is observed for D f = 2.0. This non-monotonic dependence on D f can be explained with the aid of Fig. 14. To generate Fig. 14a, the ELS signals were computed from agglomerates with a uniform particulate number, N, and fractal dimension, D f , keeping the total number of primary particulates, n 0 = 1 ⋅ 10 20 , fixed. The figure highlights changes in the ELS signal due to agglomeration and different morphologies. For small and medium sized agglomerates in the Rayleigh and Guinier regime, the signal increases almost linearly with N. This is consistent with Eqs. (6) and (7) as n k ∼ N −1 . For large agglomerates in the powerlaw regime, the scattering signal remains constant which is in line with our discussion in App. 1. As the scattered signal from large agglomerates is much higher, only the large agglomerates affect the ELS signal. We can also see that the fractal dimension has a strong and non-linear influence on the ELS signal. Here, the scattered light of compact agglomerates is much higher than that of agglomerates with low fractal dimensions. Figure 14b shows the temporal evolution of the ELS signal from a perfectly stirred reactor calculation. The initial composition is given by a burning solution for a mixture fraction of Z = 0.35 of the reactive jet simulations with 3100 ppm silane loading. The simulation time corresponds to the time a fluid element would need to travel along the centerline in our jet configuration described in the main body of the paper, and the PSR simulations thus approximate the time history of a burning notional particle in the absence of mixing. It can be seen that agglomerates with smaller fractal dimensions emit stronger signals at the beginning. After some time, however, this dependence is reversed. This behaviour can now easily be understood with data from Fig. 14a. Initially, agglomerates with small fractal dimension grow faster due to a larger collisional cross section giving faster increases in the ELS signal. After some time, however, scatter is dominated by agglomerates in the power law regime, further growth does not increase the ELS signal, and the increase in the ELS signal considerably slows  (as nucleation is ongoing and smaller agglomerates continue to exist). Since more compact agglomerates with larger fractal dimension tend to yield stronger ELS signals, simulations using larger D f will provide larger ELS signals and the initial dependence on D f is reversed.
The evolution of the ELS signal observed in Fig. 13 can now be explained: The ELS signals for larger D f rise more slowly as agglomeration is delayed. This trend is very quickly reversed for D f < 2.2 and the more compact agglomerates tend to give larger ELS signals. For the fractal dimension of D f = 2.2 the results do not quite follow the trend. Upstream, the computed ELS signal is in line with our expectations but further downstream, the signal does not surpass the values computed for D f = 2.0 . It can be hypothezised here that this is due to broadening effects further downstream. Radial (turbulent and molecular) diffusion and particulate and agglomerate dispersion decrease the local particulate number density and thus the collision frequency and counteract the expected increase in ELS signal.
The discussions with respect to the computed ELS signal show that an unambiguous comparison between measured and computed ELS signals may require a model for the solid phase that includes information not only on mean parameters for particulate number density, radius of gyration, primary particulate diameter and fractal dimension, but also on their distributions and correlations. Such a multi-variate characterisation of the nanoparticles using a sectional approach is currently beyond reach as (i) computational requirements will be huge and (ii) models describing the dynamics within this multi-dimensional space do not exist. Only two-dimensional models based on a multi-sectional moment method were proposed (e.g. Yang and Mueller 2019;Xiong et al. 1993) and may guide future research in the context of MMC-LES-PBE of nanoparticle flame synthesis for improved predictions of particulate characteristics.