Modelling of artefacts in estimations of particle size of needle-like particles from laser diffraction measurements

Manufacturing of particulate products across many industries relies on accurate measurements of particle size distributions in dispersions or powders. Laser diffraction (or small angle light scattering) is commonly used, usually off-line, for particle size measurements. The estimation of particle sizes by this method requires the solution of an inverse problem using a suitable scattering model that takes into account size, shape and optical properties of the particles. However, laser diffraction instruments are usually accompanied by software that employs a default scattering model for spherical particles, which is then used to solve the inverse problem even though a significant number of particulate products occur in strongly non-spherical shapes such as needles. In this work, we demonstrate that using the spherical model for the estimation of sizes of needle-like particles can lead to the appearance of artefacts in the form of multimodal populations of particles with size modes much smaller than those actually present in the sample. This effect can result in a significant under-estimation of the mean particle size and in false modes in estimated particles size distributions.


Introduction
Particle size measurements are crucial across many industries in the manufacturing of particulate products, such as pharmaceuticals, agrochemicals, detergents, pigments and food. Particle size and shape have a profound influence on downstream processing as well as on final product properties through a variety of product attributes, such as solubility, dissolution kinetics, flowability, etc. There are various particles sizing methods commonly employed in practice, some of them online and others offline [1][2][3]. One of the widely used techniques for measuring of particle size distribution (PSD) in dispersions is laser diffraction [4] which is typically arranged in a flow-through setting but is most often used off line as there are limits on dispersion densities due to multiple scattering. Laser diffraction measurement involves the collection of scattered light from a dilute dispersion of particles by an array of detectors placed at different spatial locations so that they cover a certain span of scattering angles θ. Since the angular dependence of the scattered light intensity originating from a particle is a function of the size and shape of the particle, as well as the orientation of the particle with respect to the incident laser beam, the particle size and shape can be inferred from the corresponding scattering intensity pattern. However, as there is typically a distribution of particle sizes across a population, the intensity pattern measured by the detectors will be a convolution of the intensity patterns from all the particles of different sizes in the dispersion.
The estimation of the PSD from the measured scattering intensity pattern (scattering intensity as a function of scattering angle) involves solving an inverse problem using a suitable scattering model which describes the scattering intensity for particles of a given shape, size and optical properties. The inversion is implemented in the software accompanying laser diffraction instruments, typically using the Mie scattering model [5] for spherical particles as a default, regardless of the shape of particles in the measured sample. This can lead to various artefacts, such as apparently multimodal distributions in PSD estimates (e.g., [6,7]) when the shape of the particles in the sample deviates significantly from spherical. This could result in misleading estimates of mean particles sizes with severe consequences for applications where the process is very sensitive to the particle sizes. This is particularly important in the pharmaceutical industry where many of the active pharma-ceutical ingredients are crystallised in needle-like habits.
In this paper, we demonstrate that multimodal PSDs can be obtained from the inversion process even when the true particle size is monodisperse, i.e., all particles are of the same size. We simulate the scattering intensity patterns for monodisperse population of needle-like particles to explain how multimodal PSD artefacts arise due to solving the inverse problem using a scattering model for spherical particles. We compute the angularly dependent scattering intensity for needle-like particles of specified optical properties using a model for infinitely long cylinders with diameters ranging from 1 to 100 micrometres. Then we solve the inverse problem of estimating PSD from the angularly dependent scattering intensity pattern using the Mie theory for spherical particles, mimicking the analysis performed when using commercial laser diffraction instruments. Since the PSD of needle-like particles is exactly specified here, the estimated PSD from the inversion can be directly compared with the actual PSD. We consider two limiting cases, when needle-like particles are assumed to perfectly aligned with flow which is perpendicular to the incident laser beam and when they adopt random orientations. We show that in either case the inversion results in estimated PSD that are multimodal where smaller modes are mathematical artefacts due to inversion and we explain how these are related to different shapes of intensity scattering patterns for needle-like particles compared to those of spheres.

Calculating scattering intensity
The scattering intensity patterns for needle-like particles will be simulated using the scattering theory for infinitely long cylinders [5]. Even though the theory was developed for infinitely long cylinders, it is applicable for needlelike particles with lengths significantly larger than their diameters [8]. Such long thin particles of approximately cylindrical shape are similar to needlelike particles often encountered in pharmaceutical and chemical manufacturing. For example, in Fig. 1 we show typical particles of cellobiose octaacetate, benzoic acid and metformin hydrochloride. As needle-like particles are modelled as infinitely long cylinders in this work, the circular cross-sectional diameters of these cylinders will be used to represent the size of needle-like particles as their length cannot be specified. The scattering intensity pattern for spherical particles will be simulated using the Mie theory (see [5] and section 1 of the supplementary information for details). The procedure for performing the calculations is described below.  Figure 1: Images of typical needle-like crystals of (a) cellobiose octaacetate, (b) benzoic acid and (c) metformin hydrochloride. Consider a detector system (sketched in Fig. 2(a)) in which a monochromatic light with wave vector k i is incident on a particle of arbitrary size and shape. The scattered light with wave vector k s is then collected at different angles θ to the direction of propagation of the incident light by an array of detectors as depicted in Fig. 2(a). Both the incident and scattered light have components parallel and perpendicular to the scattering plane (the plane containing the incident and scattered light) [5]. The scattering wave vector q is the difference between the incident and scattered wave vectors as sketched in Fig. 2(b). The magnitude of the scattering wave vector is a function of the scattering angle θ, and it is given by [9] where λ is the wavelength of the incident light. The wavelength of λ = 0.633µm was used in all calculations in this work consistent with the wavelength of red light in typical commercial laser diffraction instruments. The intensity of the scattered light I (which is a function of the scattering angle θ) is related to the intensity of the incident light I 0 by means of the phase or scattering matrix Z as [5,10] (see also section 1 of the supplementary information for details) where the magnitude of the wave vector q has been discretised 1 into j = 1, 2 . . . , M angular positions. The particle size (represented by diameter of 1 The values of scattering angles θ = 0.0015 • to θ = 180 • were used in all the calculations in this work. This is to ensure that the entire forward direction for scattering is covered. spherical or circular cross section of needle-like particles modelled as infinitely long cylinders in this work) D i , i = 1, 2, . . . , N has been discretised into N particle size classes 2 such that the characteristic diameter D i of particles whose diameters lie between D i and D i+1 is given by The number of particles whose diameters lie between D i and D i+1 is given by X(D i ). Finally, the quantity Z 11 (q j , D i ) is the first component (see section 1 of the supplementary information) of the phase matrix averaged over particles whose sizes lie between D i and D i+1 .
Writing I(q j ) as I j , Z 11 (q j , D i ) asZ j,i and X(D i ) as X i , then Eq. (2) can be written as a matrix equation as The scaling factor in Eq. (3) can be removed if the scattering intensity is rescaled by the scattering intensity I 1 measured at an angular position close to zero. That is, the zero q limit of the scattering intensity. The choice of I 1 for rescaling the scattering intensity is reasonable since the quantity Z 11 assumes a flat profile at the zero q limit for particles of different shapes and sizes and refractive indices [5]. If the first detector in the array of detectors (sketched in Fig. 1) is placed at an angular position sufficiently close to zero, then the scattering intensity at the zero q limit for the particles in the N size The quantity 1/q has the dimension of length, and it sets the length scale accessible in a light scattering experiment [9]. The smallest and largest values of θ correspond to length dimensions 1/q of order 3800µm and 0.05µm respectively. These are the same order of magnitude as the largest and smallest particle diameters of D N +1 = 3500µm and D 1 = 0.01µm (respectively) in the particle size grid used in the calculations. The angular positions θ were discretised on a geometrically spaced grid with 100 bins. This then gives a size difference of order L 2 ∆q µm (where L is an arbitrary particle size and ∆q is the q spacing). As low q (∆q ≈ 10 −4 µm −1 ) values correspond to large particle sizes, then using L = 3800µm gives a size difference of 144µm, which is about 4% of D N +1 . This is the resolution of the largest particle size. In the large q (∆q ≈ 1µm −1 ) region, which corresponds to small particles, the size difference is of order 10 −4 µm (using L = 0.01µm) which is about 1% of D 1 , and this is the resolution of the smallest particle size.
2 A particle size grid with N = 200 size classes was used in the calculations here. The particle diameters on the grid run from D 1 = 0.01µm to D N +1 = 3500µm. This covers the entire range of particle diameters in typical commercial laser diffraction instruments [11].
classes can be constructed from the first row of matrixZ as Then the rescaled scattering intensityĨ j can be constructed as The rescaled scattering intensity in Eq. (5) is the form in which the scattering intensity data is reported in typical commercial laser diffraction instruments. Hence this form of the scattering intensity will be used in subsequent analysis in this work. The components of the phase matrix take different forms depending on the shape of the particles (see section 1 of the supplementary information).

Forward and inverse problem
In a dispersion of particles of different sizes, the scattering intensityĨ * that is measured will be a convolution of the scattering intensities of the individual particles in the dispersion and the PSD X given by Eq. (2). If the PSD X of the particles in the dispersion is given, then the scattering intensity of the population can be calculated by solving the forward problem in Eq. (5). The calculated scattering intensityĨ can then be compared with the measured scattering intensity.
In reality the PSD X of the particles in a dispersion will not be known. Instead the challenge will be to estimate the PSD corresponding to a measured scattering intensity. This will involve solving an inverse problem. As mentioned earlier, the purpose of this work is to examine the effect of the model employed in solving the inverse problem on the solutions obtained. To achieve this objective, the experimentally measured scattering intensitỹ I * will be simulated with the model for infinitely long cylinders. The scattering intensityĨ calculated using the model for spherical particles (see section In this work, the inversion required to obtain the PSD X will be carried out by minimising the objective function f given as This is a weighted least square problem which is an unconstrained optimisation problem [12]. The weight function w j is given as where the quantities C 1 and C 2 are optimisation parameters with initial values C 1 = C 2 = 0. The weighting function in the objective function in Eq. (6) is necessary as the values of the scattering intensity cover several orders of magnitude over the entire q range of interest. A similar weighting function was employed previously [13] for the calculation of intensity for anti-Stokes Raman scattering but with fixed values of C 1 and C 2 .

Number and volume based PSD
The PSD defined in Eq. (5) which is calculated by solving the least square problem in Eq. (6) is number based. The number based PSD X i is defined as an exponential function of the parameter γ i as 3 Then the optimisation problem given in Eq. (6) is solved by searching for γ i (using the Levenberg-Marquardt algorithm as implemented in Matlab) which minimises the objective function f given in Eq. (6). As the Levenberg-Marquardt method is gradient based and the objective function in Eq. (6) contains local minima, then a multi-start strategy [14] is used to search for a solution which is close to the global minimum. Once a solution close to the global minimum is found, then it is used to estimate the volume based PSD X v i as described in section 2 of the supplementary information. Finally, the volume based PSD is normalised aŝ whereX v i is the volume fraction of spherical particles of size D i .

Mean particle size
The mean size of the particles in a population can be represented by various metrics depending on the application [15]. The volume weighted mean diameter D 43 is commonly reported by commercial laser diffraction instruments. The D 43 value is defined as [15] which upon using the substitution The D 43 value will coincide with the diameter of spherical particles in a monodispersed population.

Results and Discussion
Suspended needle-like particles experience various kinds of flow conditions depending on the vessel shape, agitation arrangement, solid loading and so on. Depending on these conditions, the particles may be able to perform completely random or partially restricted rotations or they could become fully aligned with the flow field. Flow-through cells typically used in laser diffraction instruments are a few millimetres wide where the incident laser beam is perpendicular to the flat transparent windows enclosing the cell and to the direction of flow through the cell. Depending on the shape, size and solid loading of particles, their rotations in such flow through cells may be restricted and in the limit of thin needles or thin platelets at high The needle-like particles were modelled as infinitely long cylinders solid loadings they can become fully aligned with the flow with their main axis being perpendicular to the incident laser beam. Motivated by this, we consider two limiting cases: case I, where the incident beam is perpendicular to the main axis of the particle, and case II, where the incident beam is at a random angle to the main axis of the particle. These two cases will be examined in detail below. When the particles are free to perform random rotations, the incident angle of the incoming monochromatic light will take all possible values from grazing to normal incidence with equal probability. However, when the particles are aligned with the flow field, the incident laser beam will be perpendicular to the particle main axis and hence the incident angle is fixed at 90 • (see Fig. 1 of the supplementary information for a schematic for light scattering by an infinitely long cylinder.).
The refractive indices of materials typically encountered in organic crystals (e.g., in pharmaceutical manufacturing) are around N r = 1.50 with zero absorption. For example, the cellobiose octaacetate, benzoic acid and metformin hydrochloride crystals shown in Fig. 1 have refractive indices of N r = 1.51, 1.50 and 1.58 respectively. The three materials have poor solubility in methanol which has a refractive index N r = 1.33. Hence the refractive index of N r = 1.50 was used in the simulation of the scattering intensities of needle-like particles (modelled as infinitely long cylinders), and the particles were assumed to be suspended in a non-solvent medium with refractive index N r = 1.33. The wavelength of the incident light was fixed at λ = 0.633µm in all the calculations which is consistent with the wavelength of red light in typical commercial laser diffraction instruments.

Case I: restricted rotations
The calculated scattering intensities for monodisperse needle-like particles with the diameter of circular cross-section of D = 1µm, 10µm, and 100µm are shown in Fig. 3 (black diamonds). The blue circles in Fig. 3 show the calculated scattering intensity for monodisperse spherical particles with diameter equal to the diameter of circular cross-section of needle-like particles. It can be seen in Fig. 3(a) that for D = 1µm both scattering patterns are quite similar for lower q values, while in the high q (q 5µm −1 ) region the scattering intensity for spheres is much lower than that for the needle-like particles. However, we can see that the scattering intensity pattern for monodisperse needle-like particles can be reasonably well fitted using scattering from a multimodal distribution of spherical particles. The fitted curve shown by the green crosses in Fig. 3(a) is obtained by solving the optimisation problem in Eq. (6).
In order to quantify the degree of fit, we compute the L 2 norm (L f it 2 ) of the difference between the scattering intensities of the needle-like particles (modelled as infinitely long cylinders) and that of the fitted curve from the spherical model. This gives L f it 2 = 0.142. This value is nearly one third of the L 2 norm (L sphere 2 ) of the difference between the scattering intensities of the needle-like particles and that of the spherical particles whose diameters are equal to those of the circular cross-sectional diameter of the infinitely long cylinder. The value of L sphere 2 is 0.415 for the case of D = 1µm. The multimodal PSD (in terms of volume fractionsX v ) of spherical particles is shown in Fig. 4(a). In order to achieve a reasonable fit (quantified by the L f it 2 value), the estimated PSD needs to contain at least two modes, with In the low q region, the normalised scattering intensityĨ is independent of q regardless of particle shape. In the Guinier region, the scattering intensity starts decreasing with increasing q and can be approximated asĨ(q)/Ĩ(0) = 1 − q 2 R 2 g /3 [16], where R g is the radius of gyration of the particles andĨ(0) is the zero q limit ofĨ which is indicated as I 1 in Eq. (4). The Guinier region extends towards qR g ≈ 1 where a noticeable downturn of the scattering  Fig. 4(c)). The black diamonds are the scattering intensity for needle-like particles of circular cross-sectional diameter D = 100µm. The needle-like particles were modelled as infinitely long cylinders intensity can be observed. For spherical particles R g = D (3/20) so that the downturn in the scattering intensity can be seen when q ≈ (2/D) (5/3) (i.e., at q ≈ 3µm −1 for D = 1µm). At large q values the scattering intensity transits into the Porod region where it decays with an exponent depending on the nature of the particle surface (for smooth surfacesĨ ∼ q −4 [9]).
It can be seen in Fig. 3(a) that the scattering intensity pattern for the monodisperse needle-like particles (modelled as infinitely long cylinders) of circular cross-sectional diameter D = 1µm (black diamonds) turns down at slightly lower q values than for spheres of the same diameter. Therefore a slightly larger sphere diameter is needed to reproduce the Guinier region of a thin cylinder and the estimated PSD (green crosses in Fig. 4(a)) contains a peak at a particle diameter slightly larger than 1µm. The scattering intensity pattern for thin cylinders of circular cross-sectional diameter D = 1µm initially decays faster than that of spheres of the same diameter at higher q values ( Fig. 3(a)). However, the decrease of the scattering intensity pattern for the thin cylinders slows down considerably as q increases, whereas the corresponding scattering intensity for spheres continue to decrease with increasing q. Hence, the scattering pattern of a thin cylinder cannot be fitted with just that of a single sphere. In order to get a reasonable fit (judging by the L f it 2 value), it is necessary to include additional scattering at larger q values, which can only come from much smaller spherical particles and the resulting PSD will become multimodal. This can be seen as a second peak at D ≈ 0.2µm in the estimated PSD in Fig. 4(a). This peak is a pure artefact in the estimated PSD as it does not correspond to any physical length scale related to the particle size and arises solely from the difference of scattering patterns between cylinders and spheres. The same phenomenon can be seen for other cylinder sizes as shown in Figs  To further understand that using the spherical model requires a multimodal PSD to fit the scattering intensity of the needle-like particles, consider  Fig. 4(c). A scattering intensity corresponding to the fitted curve shown by the green crosses in Fig. 3(c) can be constructed by performing a superposition of the scattering intensities for monodisperse spherical particles of sizes D = 0.8µm, D = 1.2µm, D = 1.9µm, D = 21µm and D = 140µm shown in Fig. 5. The superposition is performed by weighting the scattering intensities by their respective number fractions estimated asX s =X v s /D 3 s where subscript s represents particle size andX v s is the volume fraction corresponding to each size as shown in Fig. 4(c). The normalised scattering intensity obtained by this superposition procedure is shown by the green crosses in Fig. 5, and it has a good fit (L f it 2 = 0.182) to the original scattering intensity (black open diamonds in Fig. 5) for the needle-like particles (modelled as infinitely long cylinders) of circular cross-sectional diameter D = 100µm. This fit cannot be obtained if the weighted scattering intensities from the smaller particle sizes of D = 0.8µm, D = 1.2µm, D = 1.9µm and D = 21µm are not included in the superposition.

Case II: random rotations
The Case II where the needle-like particles are able to make random rotations is analysed in this section. As in Case I, the scattering intensity pattern of the needle-like particles will be simulated with the model for infinitely long cylinders, while the inverse problem will be solved with the model for spherical particles. In this case, the scattering intensities are computed for various values of incident angle from 1 • to 90 • . Then the overall scattering intensity is obtained by averaging these scattering intensities over all angles.
Compared to Case I, the Guinier region of the scattering intensity (blue circles in Fig. 6) for spherical particles of diameter D = 1µm matches that of needle-like particles (black diamonds in Fig. 6) of the same circular crosssectional diameter much more closely. This leads the spherical model to estimate the largest particle size with the diameter which is very close to the circular cross-sectional diameter of the needle-like particles (modelled as infinitely long cylinders) as shown by the green crosses in Fig. 7. However, similar to Case I, the decay of the scattering intensity pattern for needle-like particles slows down at higher q values when compared with those of spheres of the same cross-sectional diameter as seen in Fig. 6. Hence as in Case I, in order to get a reasonable fit (based on the L f it 2 value), it is necessary to include additional scattering at larger q values, which must come from particles with much smaller diameters and the resulting PSD thus become multimodal, although the smaller size peaks are artefacts in the estimated PSD. As before, the consequence of these additional peaks introduced into the estimated PSD by the spherical model is that the estimated D 43 values can be significantly lower that the actual cross-sectional diameter of needlelike particles as seen in Fig. 7(d).

Conclusions
We have demonstrated the origin of artefacts which arise due to applying a model for scattering by spherical particles to solve the inverse problem for laser diffraction when the scattering intensity pattern comes from a population of needle-like particles. The scattering intensity patterns for the needlelike particles were simulated with the model for infinitely long cylinders using a monodisperse distribution (in terms of the circular cross-sectional diameters) of these cylinders. Our results show that using the scattering model for spherical particles it is possible to find a good fit for the scattering intensity patterns of needle-like particles, but the resulting estimated PSD is not necessarily representative of actual cross-sectional diameters of the needlelike particles. The estimated PSDs are typically multimodal with the largest size mode close to the actual cross-sectional diameter but additional smaller size modes are not physical. These modes are mathematical artefacts arising from different shapes of scattering patterns of spheres and thin cylinders. Needle-like particles of various lengths would be expected to give rise to similar effects which would depend on the aspect ratio of the particles. It can be expected that the same issue applies for particles of other shapes, such as platelets.
This situation is unavoidable as long as the scattering model for spherical particles is used to fit the data from particles of strongly non-spherical shapes as it is possible to fit essentially any scattering pattern with the scattering models for a polydiperse population of spheres. This approach is used in commercial laser diffraction instruments, and it can lead to misleading conclusions about the PSD of the particles under analysis, in terms of unrealistic multimodal distributions and underestimating the mean particle size in terms of volume weighted mean diameter values. Therefore, the way forward is to obtain information about particle shape and apply appropriate scattering models for particles to be analysed.

Calculating scattering intensity
As discussed in section 2 of the main text, the intensity I of scattered light is related to that of the incident light I 0 by means of the phase matrix. The calculation is performed by considering the detector system (sketched in Fig. 2(a) of the main text) in which a monochromatic light with wave vector k i is incident on a particle of arbitrary size and shape. The scattered light with wave vector k s is then collected at different angles θ to the direction of propagation of the incident light by an array of detectors as depicted in Fig.  2(a) of the main text. Both the incident and scattered light have components parallel and perpendicular to the scattering plane (the plane containing the incident and scattered light) [5]. The scattering wave vector q is the difference between the incident and scattered wave vectors as sketched in Fig. 2(a) of the main text.
The general information concerning the intensity and state of polarisation of both the incident and scattered light is contained in the Stokes parameters I, Q, U, V [5,10], where I is the light flux or intensity, Q and U describe the state of linear polarisation and V describes the state of circular polarisation. The Stokes parameters of the incident light are related to those of the scattered light by means of the scattering matrix or phase matrix Z as [5,10]  where k = 2π/λ is the wave number and r is distance in the radial direction. The Stokes parameters of the scattered light are indicated by the subscript s, while those of the incident wave are indicated by the subscript i in Eq.
(1). The elements Z ij , ij = 1, 2, 3, 4 of the phase matrix Z are related to the elements of the amplitude matrix S. The amplitude matrix relates the incident electric field E i to the scattered electric field E s as [5,10] where the incident light is taken to propagate in the z direction. The elements S i , i = 1, 2, 3, 4 of the amplitude matrix are functions of the scattering angle θ, the size of the particle D (the size of a particle is taken to be the diameter of spherical particles or the diameter of circular cross section of needle-like particles which are modelled as infinitely long cylinders in this work) and the refractive index N R of the particle. Most commercial laser diffraction instruments only measure the intensity of scattered light from an incident unpolarised light. Hence it is assumed in this work that there is no polarisation of the scattered light so that the Stokes parameters Q = U = V = 0. Thus Eq. (1) reduces to where Z 11 is given by [5] In a population of particles of different sizes and shapes, the phase matrix for the population will be the sum of the phase matrices of the individual particles in the population. Then the phase matrix for the population is given as [5,10] where N is the number of particles in the population. The sum in Eq. (5) holds provided that the particles are randomly positioned and oriented, and that the particles are sufficiently spaced such that there is no multiple scattering. Hence the intensity of scattered light from the population of particles will be given as where I 0 and I are the intensities of the incident and scattered light respectively and Z 11,k is the contribution to the Z 11 component of the phase matrix from particle k. The conditions required for Eq. (5) to hold can be met in commercial laser diffraction instruments, hence the conditions will be assumed to apply in this work.
Here, the polydispersed population of particles shall be assumed to consist of particles of the same shape and refractive index but different sizes.
The distribution of particle sizes is characterised by a probability density function n(D) such that n(D)dD is the probability of finding particles with diameters between D and D + dD. If the diameters of the particles are discretised and grouped into N geometrically spaced size classes such that the characteristic diameter of particles whose diameters lie between D i and D i+1 (i = 1, 2, . . . , N ) is given by D i = √ D i D i+1 , then a discretised particle size distribution (PSD) X(D i ) can be defined such that X(D i ) is the number of particles whose diameters lie between D i and D i+1 .
Then the number of particles N i with characteristic diameter D i is given by N i = X(D i ). Then Eq. (6) can be rewritten as where the scattering angle θ has been discretised into j = 1, 2 . . . , M angular positions and Z 11 (θ j , D i ) is the averaged phase matrix component for particles whose sizes lie between D i and D i+1 . Since the magnitude of the scattering wave vector q is a function of angle (as defined in Eq. (1) of the main text), then the scattering intensity and component Z 11 of the phase matrix in Eq. (7) can be rewritten as functions of q (discretised) as Writing I(q j ) as I j , Z 11 (q j , D i ) asZ j,i and X(D i ) as X i , then Eq. (8) can be written as a matrix equation as The scaling factor in Eq. (9) is removed by rescaling the scattering intensity by the scattering intensity I 1 to give the rescaled scattering intensitỹ I j asĨ

Population of spherical particles
The components S 3 = S 4 = 0 in the amplitude matrix defined in Eq. (2) in the case of spherical particles [5]. The components S 1 and S 2 are defined as [5] S 1 = n 2n + 1 n(n + 1) (a n π n + b n τ n ) (11a) where the quantities π n and τ n are functions of the scattering angle θ. They can be obtained by the following recurrence formulas [5] where µ = cos(θ), π 0 = 0 and π 1 = 1. The series in Eqs. (11) and (12) are truncated after n c terms; the value of which is related to the size of the particle [5]. The scattering coefficients a n and b n are functions of the particle size and refractive index. They are obtained by solving the Maxwell equations for a spherical particle with appropriate boundary conditions 1 . Using Eq. (4), then the Z 11 component of the phase matrix can be calculated and subsequently the rescaled scattering intensity of the population can be calculated using Eq. (10).

Population of needle-like particles
The scattering intensity of needle-like particles is approximated with that for infinitely long cylinders in this work. The sketch of such an infinitely long cylinder is shown in Fig. 1. The axis of the cylinder lies along the z-axis. The incident light (with wave vector k i ) which is contained in the x − z plane makes an angle ζ with the cylinder axis, while the scattered light (with wave vector k s ) is contained in the x − y plane and makes an angle θ with the negative x-axis. The scattered light is measured in the x − y plane by an array of detectors as shown in Fig. 1. The scattering plane contains the cylinder axis and scattered light [5].
The amplitude matrix relates the incident electric field E i to the scattered electric field E s as [5] E s E ⊥s = e i3π/4 2 πkr sin(ζ) e ik(r sin(ζ)−z cos(ζ)) T 1 −T 3 The quantities T i , i = 1, 2, 3 in Eq. (13) correspond to the quantities S i , i = 1, 2, 3, 4 in Eq. (2) for the general case. The quantities T i , i = 1, 2, 3 are given as [5] T 1 =b 0I + 2 ∞ n=1b nI cos (nθ) (14a) nII cos (nθ) (14b) The scattering coefficientsb 0I ,b nI ,â 0II ,â nII ,â nI (similar to the case of a n and b n in Eq. (11)) are functions of particle size (the diameter of the infinitely long cylinder) and refractive index. They are obtained by solving the Maxwell equations for the cylindrical geometry with appropriate boundary conditions 2 .
Using the components T i , i = 1, 2, 3 of the amplitude matrix for the infinitely long cylinder, then the Z 11 component of the phase matrix can be calculated as Subsequently the rescaled scattering intensity in Eq. (10) can be calculated for the infinitely long cylinders that are used to represent needle-like particles in this work.

Number and volume based PSD
The PSD defined in Eq. (7) which is calculated by solving the least square problem in Eq. (6) of the main text is number based. The number based PSD X i is defined as an exponential function of the parameter γ i as X i = e γ i , i = 1, 2, . . . , N.
Then the optimisation problem given in Eq. (6) of the main text is solved by searching for γ i (using the Levenberg-Marquardt algorithm as implemented in Matlab) which minimises the objective function f given in Eq. (6) of the main text. As the Levenberg-Marquardt method is gradient based and the objective function in Eq. (6) of the main text contains local minima, then a multi-start strategy [14] is used to search for a solution which is close to the global minimum. This involves using different random starting solutions for γ i , then the solution for which the L 2 norm given as is minimum is then chosen as the optimum (the scattering intensityĨ * is that of the needle-like particles which are modelled as infinitely long cylinders, while the scattering intensityĨ is that of spherical particles). Subsequently the number based PSD X i in Eq. (16) is calculated from this optimum solution for γ i . However, as commercial laser diffraction instruments typically report a volume based PSD, then it is necessary to calculate a corresponding volume based PSD. This can be achieved as follows. Consider the scattering intensity I j (obtained in a manner similar to the case of I j in Eq. (9)) given as where α = I 0 k 2 r 2 The scattering intensity I j is associated with the number based PSD X i by means of Eq. (18). However, the scattering intensity I j can also be associated with the volume based PSD X v i (the total volume of particles with diameters between D i and D i+1 ) by writing where Then the scattering intensity I can be normalised to remove the scaling factor α (as in the case of Eq. (10)) aŝ to the global minimum of the objective function f in Eq. (6) of the main text, then an initial estimate X v0 i of the volume based PSD which is close to the global minimum of the objective function f v in Eq. (27) can be constructed by the method of truncated singular value decomposition (TSVD) [14] using the scattering intensityÎ j in Eq. (25). This initially estimated volume based PSD X v0 can then be used to obtain an initial starting solution for γ v i and passed on to the Levenberg-Marquardt algorithm to solve the optimisation problem in Eq. (27). The TSVD method can only be used to generate an initial estimate for γ v i as it predicts PSDs with negative values because of the ill-conditioning of the inverse problem. However, the formulation of the volume based PSD X v i in Eq. (26) guarantees non-negative values. Finally, the volume based PSD is normalised aŝ