Nanoparticle size distribution from inversion of wide angle X-ray total scattering data

An increasingly important issue in nanoscience and nanotechnology is the accurate determination of nanoparticle sizing. Wide angle X-ray total scattering (WAXTS) data are frequently used to retrieve the Particle Size Distributions (PSDs) of nanocrystals of highly technological relevance; however, the PSD shape typically relies on an a-priori assumption. Here, we propose a modified version of the classical iterative Lucy-Richardson (LR) algorithm, which is simple, fast and highly reliable against noise, and demonstrate that the inversion of WAXTS data can be profitably used for recovering accurate PSD regardless of its shape. Computer simulations based on the use of the Debye Scattering Equation (DSE) modelling WAXTS data show that the algorithm is capable of recovering accurate PSDs even when the sample is made of a mixture of different polymorphs and/or exhibits microstrain effects. When applied to the inversion of WAXTS data taken on real samples, the method requires accurate modelling of the nanoparticle crystal structure, which includes structural defects, microstrain and surface induced distortions. Provided that this information is correctly fed to the program, the inversion algorithm reconstructs the WAXTS data with high accuracy and recovers highly robust (against noise) PSDs. Two examples reporting the characterization of Magnetite-Maghemite and commercial P25-Titania nanopowders, are discussed. We demonstrate that pre-assumption of wrong PSD shape leads to inaccurate number-based average sizes in highly polydisperse samples.


1) Inversion Algorithm: optimization of the smoothing parameter
As explained in the main text, the novelty of our algorithm, with respect to the original Lucy-Richardson algorithm, is the introduction of a smoothing procedure of the recovered distributions between successive iterations. This procedure consists in convolving the recovered distribution (for each phase p) with a 3points symmetric triangular operator Λ , in which the amplitudes of the lateral points are set by the parameter . Such a convolution acts as a regularization scheme that produces smooth distributions, but the price to pay for such regularization is a dampening of the ultimate resolution of the method. Indeed, the higher , the broader the distribution that can be reliably recovered. In other words, narrow distributions are artificially broadened if a too large value for is used, dampening in this way the resolving power of the method.
An example of these effects is shown in Figure 1, where the noisy data (~300) deriving from two Log-Normal distributions of anatase TiO2 nanocrystals with the same (weight) average diameter < 1 > = < 2 > = 15.6 , but quite different widths 1 = 6.2 ( Fig.S1(a)) and 2 = 0.50 [ Fig.S1(a)], were inverted with different values of (including = 0) varying in the range 10 −6 − 10 −1 . As one can immediately notice, the recovered distributions of Fig.S1(a) are fairly dependent on , passing from highly spiked curves at small to nicely smoothed curves for large 's where, however, the reconstruction is not so accurate. In between, there is an optimal value of ~10 −3 for which the matching between the recovered and the input distribution is reasonably good, as also witnessed by the residual plots reported in Fig.S1(b). When the input distribution is quite narrow as in Fig.S1(c), the optimal value tend to be ~0 and as larger and larger 's are used, the recovered distribution becomes increasingly over-smoothed, spoiling the resolution of the inversion algorithm. Thus, it is quite evident that, for this simulation, the optimal value of is highly dependent on distribution width and must be optimized. Another factor which is expected to influence the optimal value of to be used in the inversion procedure is the actual noise present in the data. To quantitatively investigate this effect, we repeated the same simulations of Fig.S1(a) and Fig.S1(b) by generating independent noisy WAXTS-DSE signals, all of them characterized by the same level of statistical Poisson noise (~300), which is equal to the typical noise encountered in total scattering experiments of solid (i.e. dry) nanoparticles performed at dedicated synchrotron beamlines. Each input signal was inverted with different values of and the inversion procedure was stopped as described below in the second paragraph of this SI. The agreement between the distributions recovered at the different and the input one, was ascertained by using the Relative Restoration Error (RRE) parameter defined as (S1) where the mass distributions are computed as ( ) = ( ) (note that the phase index p has been omitted for simplicity). Eq.(S1) corresponds to the relative average mean square deviations between the retrieved and input mass distribution. Thus, the optimal value of is, in principle, the one that minimizes RRE. Relative Reconstruction Error (RRE) between the input and recovered distributions; each curve exhibits a minimum at * whose average value is 〈 * 〉 = 1.04 × 10 −3 ± 5.4 × 10 −4 (see horizontal bar). (b) Goodness of Fit (GOF) parameter (see Eq.5 main text); (c) Relative error [(rec-inp)/inp] between the recovered and input average diameters. (d) Relative error between the recovered and input standard deviations. The horizontal bar indicates the range 〈 * 〉 ± . Vertical bar in (b) indicates the value ~10 〈 * 〉 at which the GOF starts to deviate from 1. Vertical bars in (c) and (d) show that both errors start to deviate from 0 at values well beyond ~10 〈 * 〉 .
In Fig.S2(a) we report the behavior of RRE as a function of for the same distribution of Fig.S1(a) corresponding to 10 statistically independent configurations of Poisson noise. As one can notice, the curves change remarkably for different configurations and exhibit minima, * , that are very broad and spread over a pretty large range of 's. Their average value is 〈 * 〉 = 1.04 × 10 −3 ± 5.4 × 10 −4 . In Fig.S2(b) we report the behavior of GOF as a function of corresponding to the same curves shown in Fig.S2(a). This figure suggests that there is an extremely large range of 's (from ≪ 〈 * 〉 to ~10 〈 * 〉, indicated by the vertical bar) where, regardless of the fact that RRE might be quite higher than its minimum value, the signal reconstruction is always excellent with values of GOF~1. Similarly, within this range, the distribution recovery is always accurate, as evidenced by Fig.S2(c) and S2(d), where the relative errors [(rec.-inp)/inp] between the recovered and input average diameters and standard deviations are reported as a function of . Both figures exhibit the same behaviors, showing that, as long as ≤ 10 〈 * 〉, both parameters are recovered quite accurately, with relative errors always smaller than ~0.1% (average diameter) and ~1% (standard deviation). Figure S3 reports the same analysis of Fig.S2 for the distribution of Fig.S1(b), which is much narrower. In this case, the minima * are even more scattered and, on average, are much smaller, i.e. 〈 * 〉 = 3.2 × 10 −7 ± 3.2 × 10 −7 Similarly to Fig.S2, as long as ≤ 10 〈 * 〉 (vertical bar), GOF~1 and the relative errors on the average diameter remain always smaller than ~0.1%. As to the standard deviation, Fig.S3(d) suggests that the optimal value of → 0, but also in this limit the errors on remain always Figure S3 -Behaviors, as a function of , of various parameters characterizing the distribution of Fig.S1(b) for 10 different configurations of statistical Poisson noise of the same level ( = 300). (a): Relative Reconstruction Error (RRE) between the input and recovered distributions; each curve exhibits a minimum at * whose average value is 〈 * 〉 = 3.2 × 10 −7 ± 3.2 × 10 −7 (see horizontal bar). (b) Goodness of Fit (GOF) parameter (see Eq.5 main text); (c) Relative error [(rec-inp)/inp] between the recovered and input average diameters. (d) Relative error between the recovered and input standard deviations. The horizontal bar indicates the range 〈 * 〉 ± . Vertical bar in (b) indicates the value ~10 〈 * 〉 at which the GOF starts to deviate from 1. Vertical bars in (c) and (d) show that both errors start to deviate from 0 at values well beyond ~10 〈 * 〉 .
~10%, a limitation due to the intrinsic finite resolution of the inversion procedure. On the other side, for values of > 〈 * 〉, the figure shows that, up to ~5 〈 * 〉, the errors on are always ≤ 20%.
Summarizing, Figs.2 and 3 show that, although the data associated to each single noise configuration should be inverted with their own optimal * , if we use a unique single optimal = 〈 * 〉, the errors introduced in the recovery of the distribution parameters are quite negligible (except for the parameter of the narrow distribution, which remain always ~10%, due to the intrinsic finite resolution of the procedure).
Clearly, we expect that is also dependent on many other parameters such as distribution shapes, crystal phases and the absolute level on noise present on the data. For this reason, we repeated the same test of Figs.2 and 3 with different levels of Poisson noise, by using crystals of anatase and rutile, characterized by Log-Normal and Gaussians distributions with various average diameters and standard deviations. Our findings show that, regardless of all the other parameters, is mainly dependent on only two parameters: the width of the (mass) recovered distribution and the noise level present in the data. Figure 4 reports the behavior of (symbols) as a function of for a variety of different distributions and four noise levels ( = 30, 100, 300, 1000). As one can notice, for each , the values of span over about two decade, passing from ~10 −7 − 10 −5 for very narrow distributions (~0.5 − 1 ) up to values which tend to saturate around ~10 −3 − 10 −1 for broad distributions with ≥ 10 . This range of variability of about two decade is mostly due to different noise levels, with noisiest data ( = 30) requiring values of ~ 100 times higher than those required for least noisy data ( = 1000). Figure S4 -Behaviors, as a function of , of the optimal smoothing parameter for a variety of different distributions and four noise levels. The lines are guide to the eye.
In conclusion, Fig.S4 represents a clear indication for the choice of when and SNR are known. However, when these two parameters are not (or only roughly) known prior the inversion, a rule of thumb for choosing is the following: (a) carry out the inversion with = 0; [which is not affected by → 0 as shown in Figs.S2(d) and S3(d)]; (c) by using Fig.S4, choose a mid-range value for corresponding to ; (d) repeat the inversion with the new value of and compute the new ′ ; (e) if ′ > decrease (for example by a factor 2) and go back to point (d); (f) if ′ ~ and the recovered distribution is still too spiky (see below), increase (for example by a factor 2) and go back to point (d); (g) if ′ ~ and the recovered distribution is sufficiently smooth (see below), accept the result and the procedure is over.
As a final comment about the spikiness (or smoothness) of the recovered distribution, we found that, although possible, a quantification of this feature was not necessary. Indeed, as long as < (see panels (c) and (d) of Figs.S2 and S3), both < > and are unaffected by .Thus, distributions with different levels of spikiness provide the same recovered parameters and selecting the one which is sufficiently smooth is a simple criterion of "good sense" based on a visual inspection of the curve.

2) Inversion Algorithm: stopping criterion
The iterative procedure was stopped according to the following criteria: (a) First of all, we impose a minimum number of iterations, , which is necessary for the algorithm to work properly, i.e. to reconstruct accurately the expected distribution under ideal conditions (noiseless data). This is necessary because the starting uniform distribution is (obviously) very different from the expected one and the LR algorithm attains convergence quite slowly. The parameter was estimated by finding the number of iterations necessary for retrieving the expected (mass) distributions with high accuracy (RRE ~10 −3 ). We tested many distributions of a single phase (anatase) TiO2 nanocrystals with different shapes (Log-Normals and Gaussians) and different average diameters and standard deviations. All the inversion were carried out by imposing = 0, which is the optimal value for noiseless data. The results indicate that, regardless of average diameters and distribution shapes, is mainly dependent on the (mass) standard deviations the input distributions as described in Table S1. Notice that scales approximately as the inverse of the distribution width ( ∝ 1/m) and for the narrowest distributions ( m = 1 nm) tends to be rather high (> 10 4 ).
(b) We also impose a maximum number of iterations = 10 6 , which ensures that the inversion stops even when the criteria below reported are not met.
(c) For any < < , the procedure is stopped as soon as the parameter GOF(r), computed for each iteration, attains a minimum and continues to increase for at least 10 consecutive iterations.
(d) Additionally, whenever, for < < , condition (c) is not met but the decrease of GOF(r) becomes increasingly slow, the procedure is stopped when the variation of GOF(r) is below a given threshold. Numerically, we monitored the parameter ( ) = [ GOF(r)/dr]/GOF(r) (equal to the (relative) first derivative of GOF with respect to r) and stopped the procedure when ( ) ≤ 10 −9 .
Finally, we checked that, when the procedure would prefer to stop at a number of iterations < either because conditions (c) or (d) are met, forcing it to continue up to does not jeopardize significantly the quality of the recovered distribution. Table S2 reports the comparison between the number input and recovered distribution parameters relative to the TiO2 simulation described in Fig.2 of the main text.  Table S3 reports the comparison between the number input and recovered distribution parameters relative to the Fe5Te4 simulation described in Fig.3 of the main text.

4) Inversion Algorithm: stability against noise
One of the main advantages of having introduced our smoothing procedure in the original LR algorithm is the remarkable stability of our algorithm against noise. Indeed, whereas for the original LR procedure the recovered distributions tend to become spiky when too many iterations are processed, in our case, they always remain nicely smooth without renouncing to accuracy (provided that is properly chosen). We have ascertained the stability of our algorithm against noise by repeating the same simulation test of Fig.2 of the main text (a TiO2 sample made of a mixture of nanocrystals of the three common polymorphs rutile, anatase and brookite) with three different noise levels corresponding to ~ 300, 100, 30.
For each of them, we accumulated statistics by generating (and inverting) 100 independent noisy ( ) scattering profiles obtained by adding 100 statistical independent configurations of Poisson noise to the ideal (noiseless) profile. Then we evaluated the algorithm performances in terms of accuracy on the recovered distribution parameter. The results of this test are shown in Figs.5 and 6.  Fig.S5(b), where the (systematic) errors are always smaller than 1%. On the right side of the figure, we report for all the three noise levels, as a function of the noise configuration number, the GOF (c), the stopping iteration number (d), and the relative error on the recovered background amplitude (e). As one can notice, the GOF is always around unity (±0.02), the number of iterations varies between ~1 − 5 × 10 4 , and the relative error on the background amplitude is always lower than ~1%.   Figure 6 reports for all the three noise levels, as a function of the noise configuration number and for all the three phases, the relative errors between the recovered and expected (mass) average diameter < > (first row), standard deviation (second row) and concentration (third row). As expected, the errors increase with the noise level (lower SNRs), but for common experimental conditions where typically ≥ 100, they are always much less than 1% for < > and , and less than a few percents for .
Therefore, we can conclude that the stability of our algorithm against noise is remarkably high.

5) Inversion algorithm: artefacts in the recovered distributions deriving from imperfect modeling
In this section we report some examples of the artefacts that might arise in the recovered distribution when there is some discrepancy between the modelled and actual structure of the nanocrystals or there are some errors in the (shape of the) background signal. We anticipate that these inaccuracies introduce systematic errors in the kernel functions that show up as spurious peaks in the recovered distributions. Figure 7 shows an example of what happens when the background signal is not properly taken into account. The input data have been generated by summing the WAXTS-DSE data ( Fig.S7(b), black curve) corresponding to a distribution of anatase TiO2 nanocrystals (Log Normal, < > = 24 , = 6.0 ) to a linear background (blue curve). The latter one is quite small with respect to the scattering peaks, but not so small with respect to the diffuse scattering in between the peaks. The inversion has been carried out by using either no background or a constant background. The corresponding recovered mass distributions (colored symbols) are shown in Fig.S7(a), together with the input distribution (black solid curve). As evident, the recovered distributions match the input one rather accurately over almost the entire diameter range, except for the very narrow sizes where two spurious peaks are present. These peaks are due to the fact that, since the background is not available or not properly shaped for the signal reconstruction, this missing contribution is attributed to the smaller sizes that are the ones whose scattering profiles resemble more closely the background signal (see magenta curve in Fig.S6(b), which corresponds to the second lowest size used in the nanocrystal distribution). It is worth noticing that the two spurious peaks are very well separated from the main bell shaped distribution; thus, they can be easily trimmed out when computing the distribution average parameters. For completeness, we report in Fig.S7(c) also the input and recovered number distributions, in which the spurious peaks are enormously amplified. Nevertheless, as for the mass distributions, they are still well isolated and can be easily removed from the analysis.
The second example shows what happens when the Debye−Waller thermal parameters are not correct. Figure 8a compares the input distribution of anatase TiO2 nanocristyals (Log Normal, < > = 12 , = 4.4 ) with the distributions recovered when the input data have been generated with = 0.7 Å, but inverted with = 0.5 Å (red circles) and vice versa (green squares). The effects of varying between 0.5 Å and 0.7 Å is shown in the figure inset [ Fig.S8(b)], where is quite evident that the main difference between the two profiles shows up in the diffuse scattering between the peaks. In these regions, the sample with higher scatters slightly more (~10%) than the sample with smaller . Thus, the situation is similar to that discussed in Fig.S7, where an incorrect background signal was used. It is indeed well known that thermal factors and the background level, even in conventional powder diffractometry, highly correlate, and manifest themselves in additional scattering, attributable to diffuse or extra-sample contributions, respectively. Correspondingly, the distribution of Fig.S8(a), which was recovered by using kernel functions with not enough scattering in between the peaks (very much as in the case of an underestimated background level) shows a spurious peak at small sizes (red circles). On the contrary, had the input data generated with = 0.5 Å, and inverted with = 0.7 Å, an excess diffuse scattering appears, which can be dampened only by reducing the population of the smallest classes belonging to the input distribution (green squares). In the latter case the occurrence of this artefact is less evident, and care must be taken into interpreting the inversion results. For completeness, we report in Fig.S8(c) also the input and recovered number distributions, in which the two artefact are remarkably amplified with respect to the ones shown in the mass distributions. Thus, the recovery of number distributions in the presence of this kind of artefacts becomes highly critical.

6) Inversion algorithm: comparison with DEBUSSY analysis.
In this section we report two examples concerning the comparison between our algorithm and the DEBUSSY analysis. As known, the DEBUSSY method pivots on the strong assumption that the PSD of the sample is supposed to be known (typically a LogNormal distribution) and the various distribution parameters (average size, standard deviation, Debye-Waller thermal parameters, background amplitude, etc.) are retrieved by standard X 2 minimization. Thus, in those cases where the PSD shapes are fairly different from a LogNormal, this wrong assumption could affect significantly the results, and the fitted distribution parameters might be highly inaccurate.
In Fig.9 we report this comparison for the case of a Weibull (number) input distribution of anatase TiO2 nanocrystals with 〈 〉 = 10 and = 5 . Figure 9(a) and (b) compare the input number and mass distributions (black solid curves) with the corresponding averaged distributions recovered with the DEBUSSY analysis (blue squares) and with our algorithm (red circles). Statistics was accumulated by averaging 100 distributions retrieved by processing 100 noisy data generated with = 300. As one can notice, the error bars are always remarkably small (not visible on the scale of the figure) except for the small sizes of the number distributions recovered with our algorithm (where, anyway, there is statistical consistency between input and recovered distributions). Thus, Figures 9(a) and 9(b) show that, whereas our algorithm is capable of retrieving both number and mass distributions with high accuracy, the DEBUSSY method recovers with a somewhat accuracy only the mass distribution, but wildly fails in the reconstruction of for the number distribution. The different performances of the two methods are also witnessed by the different qualities of signal reconstruction: with our algorithm we obtain a GOF ~ 1.00 ± 0.01 whereas with the Debussy analysis we get GOF ~ 3.16 ± 0.02. Although such differences are not appreciable in Fig.9(c) where both reconstructed signals are indistinguishable from the input signal (black circles), the residuals plot clearly shows that the DEBUSSY reconstruction (blue curve) exhibits systematic deviations that, in correspondence of the peaks, are much higher than the (non systematic) deviations associated to our method (red curve). A summary of the results of this test together with more tests carried out with = 100 and = 30 are reported in Table S4. Regardless of the , our method recovers both number (〈 〉 , ) and mass (〈 〉 , ) parameters that are always consistent (within the statistical accuracy) with the input ones. Conversely, for the DEBUSSY analysis, 〈 〉 and are both wrong by ~20% whereas the mass parameters are somewhat more accurate with only off by ~8%. It is also worth noticing that the GOF recovered with our method is always around unity, whereas for DEBUSSY depends sensitively on the levels. As expected at low = 30, the GOF is around unity because the systematic deviations between input and reconstructed data are smaller than noise (data not shown), but at high = 300, the opposite takes place (see Fig.9d) and we get GOF ~ 3.16 ± 0.02. and corresponding (averaged) recovered distributions obtained with the DEBUSSY analysis (blue squares) and with our inversion algorithm (red circles). Statistics was accumulated by processing 100 noisy data with ~ 300.; (b) corresponding input and recovered mass distributions; (c) Simulated input WAXTS (black circles) and reconstructed data obtained with the DEBUSSY analysis (blue line, not visible)and with our inversion algorithm (red line); (d) absolute residuals (recovered-input) for the data of panel c. DEBUSSY residuals are systematic and much higher than inversion residuals. The second test refers to an exponential decay (number) input distribution of anatase TiO2 nanocrystals with nominal 〈 〉 = = 5 . As above, we processed 100 noisy input data signals with different levels and, for = 300, we obtained the results reported in Fig.10. Once again the distributions obtained with the DEBUSSY analysis are partially (mass) and highly (number) inaccurate, whereas the ones obtained with our algorithm match quite nicely the expected ones (both number and mass). The differences in signal reconstruction are very similar to the ones obtained in the previous case, with similar systematic residuals [see Fig.10(d)] and a higher GOF. A summary of the results of this test carried out also with = 100 and = 30 are are reported in Table S5. All the comments done for Table S4  apply also to Table S5.
In conclusion, we have shown two examples (a Weibull and an Exponential distribution) of the different performances between our inversion algorithm and the DEBUSSY analysis. In both cases the number distributions to be recovered were quite broad ( / 〈 〉 ≥ 0.5) and their shapes quite different from that of a LogNormal distribution, which is characterized by long decaying tails toward large sizes. Conversely, the Weibull distribution exhibits long tails toward small sizes and the Exponential distribution exhibits no tails at all at small sizes. For these two examples, the DEBUSSY method wildly fails in recovering the number distributions with errors on 〈 〉 and of ~20%, but recovers the mass distribution with a somewhat satisfactory accuracy so that, in spite of the fact the PDS shape does not match accurately the expected one, all the mass parameters 〈 〉 and are retrieved with rather high accuracy (a few percents) except for of the Weibull distribution (~8%). Therefore, the DEBUSSY analysis appears to be reliable when dealing with mass distributions, but any comparison with other techniques based on the analysis of number PSDs (such as TEM or other optical microscopy methods) must be taken with care.  and corresponding (averaged) recovered distributions obtained with the DEBUSSY analysis (blue squares) and with our inversion algorithm (red circles). Statistics was accumulated by processing 100 noisy data signals with ~ 300.; (b) corresponding input and recovered mass distributions; (c) Simulated input WAXTS (black circles) and reconstructed data obtained with the DEBUSSY analysis (blue line, not visible)and with our inversion algorithm (red line); (d) absolute residuals (recovered-input) for the data of panel c. DEBUSSY residuals are systematic and much higher than inversion residuals.
Clearly, had we used in the DEBUSSY analysis the correct shape for the PSD, we would had found much more accurate results, consistent with the ones obtained with the inversion algorithm. It should be also pointed out that, regardless of the shape, when the PSDs are rather narrow ( / 〈 〉 ≤ 0.1) the differences between our inversion method and the DEBUSSY analysis become more and more negligible. Under these circumstances, the DESUSSY analysis in more convenient because much faster than the inversion algorithm.
High-resolution Wide Angle X-ray Total Scattering (WAXTS) measurements were performed at the MS-X04SA Powder Diffraction Beamline of the Swiss Light Source (Paul Scherrer Institute, Villigen, CH). 1 Two different beam energies of 15 KeV (Fe2O3) and 17 KeV (TiO2) were set and the operational wavelengths (λ15KeV = 0.82712 Å, λ17KeV = 0.70880 Å) accurately determined using a silicon powder standard (NIST 640d, a0 = 0.543123(8) nm at 22.5°C). Data were collected in the 0.5°-130° 2θ range using a single-photon counting silicon microstrip detector (MYTHEN II). 2 The spatial coherence length of the X-ray beam of the MS-X04SA beamline is claimed to be, in the longitudinal direction, of the order of 10 5 's, i.e., a few microns, and, in the transversal plane, up to 0.1 mm. Such coherence is much larger than the sizes of nanoparticles (< 100 nm) treated with DSE equation and the inversion algorithm (Eq.s 1 and 2 of the main text). Thus, the impinging field does not suffer of significant spatial variations and does not affect the analysis.
He/air background and empty glass capillaries were independently collected under the same experimental conditions. Additionally, an amorphous ferrihydrite sample was measured in the same experimental conditions used for the MM nanocrystals (NCs), to be added as background curve during the modelling.
Angle-dependent intensity corrections 3 were applied to the raw data to account for signal attenuations due to absorption effects; sample absorption curves were determined by measuring the transmitted beam from the filled capillaries, while for the empty capillaries the X-ray attenuation coefficient was computed using their nominal composition. Angular calibrations were applied to the zero angle and to x, y capillary offsets, derived from the certified silicon powder standard (NIST 640d) using locally developed procedures. Air and (absorption-corrected) capillary scattering contributions were subtracted from the signals of the samples.

8) Modeling and scattering profiles of Magnetite-maghemite (Fe3O4 -γ-Fe2O3) nanocrystals
The kernels profiles ( , , ) used for the inversion algorithm were computed using the DEBUSSY Suite, a suite of programs implementing the Debye Scattering Equation (DSE) to model total scattering data of nanosized and disordered materials. 4 The Suite relies on a bottom-up approach that consists of two main steps. In the first one, a monovariate population of atomistic models of nanocrystals (NCs) with increasing size and desired shape are generated, with the set of multiplicities of sampled interatomic distances stored in suitable databases. In the second step, this stored information is used to compute, via the DSE equation, the kernels ( , , ).
Magnetite and maghemite are characterized by the same spinel-like crystal structure (space group Fd-3m), shown in Fig.S11, with oxygens forming a face centered cubic lattice and two crystallographic independent Fe atoms occupying octahedral and tetrahedral interstitial sites, respectively. In magnetite (Fe3O4, containing Fe 2+ and Fe 3+ in the 1:2 ratio) half Fe 3+ ions show a tetrahedral coordination, whereas the other half, together with Fe 2+ cations, occupy the octahedral sites. When magnetite is partially oxidized to maghemite (-Fe2O3), some additional iron vacancies are created. The crystal structure remains the same with a slight shrinking of the lattice ( = 8.397 Å for magnetite to = 8.346 Å for maghemite). 5 The DSE modelling strategy behind the calculation of the Magnetite-Maghemite (MM) kernel functions was similar to the one used in Ref. [5], but based on the use of constant (not size dependent) crystallographic parameters. In particular: (i) the cubic crystal structure was built with an average lattice parameter = 8.36052 Å, which is in between the ones characterizing the magnetite and maghemite structures. This figure (slightly expanded with respect to that of magnetite 5 ) accounts for surface expansion effects that are common in many oxides and are mainly due to repulsion between unsaturated ions in the NCs shell. 6 Based on this structure, a monovariate population of NCs clusters of spherical shape and increasing diameter ( = 0.65 ) were generated up to = 50 . The corresponding multiplicities of the sampled interatomic distances of each cluster were stored a suitable database.

9) Modeling and scattering profiles of commercial Titania (TiO2) nanocrystals.
Two spherical databases for the two polymorphs of TiO2 (anatase and rutile) were built, implementing the same strategy detailed for MM NCs and using the structural parameters available in literature, 7 further optimized by a Rietveld refinement using the Topas program: 8 a = 3.78596 Å and c = 9.50805 Å for anatase (space group I41/amd); a =4.59469 Å and c =2.95921 Å for rutile (space group P42/mnm).
The TiO2 kernel profiles ( , , ) for both phases were computed by using the DEBUSSY suite with the exact DSE (Eq.1) for spherical NCs up to diameters of ~80 . Above this size, the computation via the DEBUSSY Suite become rather impractical because of very long computational times.
While for the anatase phase sizes up to 80 are large enough to recover correctly the PSD, the recovered distribution of the rutile appears to be truncated and suggests the presence of larger sizes. Thus we resorted to an alternative approach based on (Rietveld-inspired) analytical pseudo-Voigt functions describing the shapes of the diffraction peaks and a polynomial description of the diffuse scattering hidden in the background baseline. All these parameters were derived upon calibration using the DSE signals for the smaller NCs (up to ca. = 80 ) as benchmarks (shown in Fig.S12).
Once derived through calibration, these parameters were used to calculate, using the TOPAS program, 8 all the kernels for rutile used in the inversion algorithm, up to ~ 200 .
The same Debye-Waller factor (0.6Å 2 ) and s.o.f.=1.0 were used for both atomic species (Ti and O) and both phases (anatase and rutile). For both phases, the inelastic Compton scattering contribution has been added as an additional model component at the final simulated patterns.

10) Inversion from multimodal distributions.
In this section we show two examples of the capability of our method in recovering bi-modal distributions that were characterized by two peaks of different widths and relative heights. Figure S13 reports the case of a bi-Gaussian distribution of anatase TiO2 nanocrystals characterized by two broad peaks, namely < 1 > = 6.88 , 1 = 1.75 , mass fraction 1 = 0.15 and by < 2 > = 23.4 , 2 = 4.65 and mass fraction 2 = 0.85. The noisy data ( = 300) were inverted by setting = 2 × 10 −4 and the iterative procedure was stopped after ~7.2 × 10 4 iterations when the relative variation of GOF attained the threshold ( ) ≤ 10 −9 (see step (d) of section 4) and GOF~1.01. As one can notice, the matching between the input (black solid curves) and the recovered (symbols) distributions is excellent, for both number (Fig.S13a) and mass (Fig.S13b) distributions. Figure S14 reports the case of a bi-Gaussian distribution of anatase TiO2 nanocrystals characterized by narrow peaks, namely < 1 > = 20.0 , 1 = 1.0 , mass fraction 1 = 0.54 and by < 2 > = 30.0.6 , 2 = 1.0 and mass fraction 2 = 0.46. The noisy data ( = 300) were inverted by setting = 1 × 10 −6 and the iterative procedure was stopped at the maximum number of iterations 10 6 where GOF~0.99. As for the case of Fig.S13, the matching between the input (black solid curves) and the recovered (symbols) distributions is excellent, for both number (Fig.S14a) and mass (Fig.S14b) distributions.

11) Ill-posedness analysis of the WAXTS-DSE data inversion problem.
We show in this section that the inversion of WAXTS-DSE data is not a severely ill-posed problem. Indeed, as mentioned in the main text, the kernels ( , , ) associated to Eq.(1) are highly structured with the presence of a large number of relatively narrow and differently shaped peaks, whose amplitude, width, and positions depend sensitively on NC size and morphology. Thus, as long as the data ( ) are taken over a large Q-range with a high Q-resolution and high ratio, the difference of the Intensity profiles of two adjacent size classes is higher than noise and the inversion algorithm can recover the correct distribution without introducing artefacts. An example of the dependence of the intensity profiles for 7 anatase TiO2 NCs with sizes ranging between 0.5 and 83 nm, is reported in Fig.S15. As one can appreciate, the average intensity, peaks height and widths vary sensitively with the NCs sizes.
A quantitative analysis of this behaviors is reported in Fig.S16, where we show that the average intensity scales as 〈 ( )〉 ~ 3 [ Fig.S16(a)    As already mentioned above, the ill-posedness of an inversion problem is highly related to the noise level present on the data in comparison with the (r.m.s.) difference between the kernels of two adjacent sizes. Fig.S17 shows the relative difference between three couples of adjacent kernels spanning the entire range of the anatase TiO2 sizes used in the reconstruction. Before taking the difference, the intensity profile of each kernel has been normalized to ( j ) 3 , so that all the kernels are rescaled to the same average intensity. As one can notice, the relative difference is pretty large for the smaller sizes (kernels 4 and 5 with 4 = 2.02 and 5 = 2.53 ), whereas becomes increasingly smaller for large sizes (kernels 163 and 164 with 163 = 82.60 and 5 = 83.11 ). The equivalent associated to such a difference is = (∑ 〈 〉 2 =1 ∑ ∆ 2 =1 ⁄ ) 1/2 where 〈 〉 = [ 1 ( ) + 2 ( )]/2 and ∆ = 2 ( ) − 1 ( ) (the suffixes 1 and 2 indicating the members of the couple). Such figures, indicated in the legend of Fig.S17, are smaller than the of the data for each couple of adjacent kernels spanning the entire size range. Thus, the inversion problem is expected to be not severely ill-posed.
Finally, when the kernels of different polymorphs are compared, the peaks of the intensity profiles show up at rather different Q values (occasionally overlapping, particularly at high Q's), implying a very high (up to ~100%) relative differences. Correspondingly, the values of couple of kernels of distinct polymorphs are very low, and the problem is clearly not ill-posed.
Figure S17 -relative difference between three couples of adjacent kernels for the anatase TiO2 NCs of Fig.S15. Small sizes exhibit relative differences much higher than those associated to large sizes. Their corresponding (see text) are always smaller than typical present on the data.