On the parametrization of optical particle counter response including instrument-induced broadening of size spectra and a self-consistent evaluation of calibration measurements

Optical particle counters (OPCs) are common tools for the in situ measurement of aerosol particle number size distributions. As the actual quantity measured by OPCs is the intensity of light scattered by individual particles, it is necessary to translate the collected distribution of scattering signals into the desired distribution of particle sizes. A crucial part in this challenge is the modeling of OPC response and the calibration of the instrument, i.e. to establish the relation between particle scattering cross-section and measured signal amplitude. To date, existing methods lack a comprehensive parametrization of 5 this instrument response, particularly regarding the instrument-induced broadening of size distribution widths. We introduce an approach overcoming the present shortcomings by implementing a simple parametrization of the broadening effect and a self-consistent way to evaluate calibration measurements using a Markov chain Monte Carlo (MCMC) method. We further outline how to obtain realistic uncertainty estimates for OPC size distributions within this new framework. Measurements of particle standards for two OPCs, the Grimm model 1.129 (SkyOPC) and the DMT Passive Cavity Aerosol Spectrometer Probe 10 (PCASP), show substantial reduction in residuals between measured and modeled response compared to existing methods. For the presented set of measurements only the new approach yields results that are conform with the true size distributions within the range of model uncertainty. The offered approach will help to improve the accuracy of derived size distributions and the assessment of their precision for OPC measurements in general.


Introduction
The size distribution of aerosol particles is a key property to understand the impact of aerosols on human health and Earth's climate.To measure aerosol size distributions, optical particle counters (OPCs) are widely used in air quality programs and atmospheric studies.However, several studies directly comparing size distributions from different OPC instruments (e.g., Belosi et al., 2013;Renard et al., 2016) and OPCs with other sizing methods (e.g., Reid et al., 2003;Müller et al., 2012) find significant disagreements and in some cases OPCs show systematic mis-sizing and artificial broadening of size spectra.This highlights that, although OPCs allow for a fast assessment of qualitative size information, the task to gain proper particle number size distributions can be challenging.One reason for this is the measurement principle itself, as particle size is only indirectly inferred from scattered light intensity.This intensity, in general, is a non-monotonic function of particle size and depends also on particle intrinsic properties, such as complex refractive index and shape (Szymanski and Liu, 1986;Szymanski et al., 2009).Especially for particle sizes that are comparable or larger than the wavelength of the incident light, the size dependency of scattered intensity tends to be flat and occasionally ambiguous, so that uncertainties in the particle's intrinsic properties can introduce large sizing uncertainties (Reid et al., 2003;Formenti et al., 2011).Another reason lies in the existing methods for OPC calibration and response parametrization.The available approaches (e.g., Cerni, 1983;Bemer et al., 1990;Rosenberg et al., 2012;Cai et al., 2013) are not consistent with each other.Further, they do not allow for a comprehensive description of instrument response and a satisfactory quantification of corresponding uncertainties.Figure 1 summarizes the major sources of uncertainty adjunct to OPC measurements.They can be divided into the last-mentioned uncertainty in the instrument response and calibration, the uncertainty in the particle intrinsic properties and the uncertainty in the measured concentrations themselves, e.g., arising from counting statistics, eventual particle losses.In order to allow for inter-comparability between different OPC instruments and the comparison with other measurement techniques, it is necessary to correct for systematic errors and to quantify all uncertainties as good as possible, i.e., to improve OPC data accuracy and assess its precision (Formenti et al., 2011;Mahowald et al., 2014).
In the following paper we focus on the central aspect of OPC response modeling and calibration and present a new approach that allows for a more accurate description of OPC instrument response and yields realistic associated uncertainty estimates.
We discuss the advantages of the new approach against the background of the prevailing concepts and present its superiority by means of measurement results for two optical particle counters that were used during the Saharan Aerosol Long-range TRansport and Aerosol-Cloud-Interaction Experiment (SALTRACE) (Weinzierl et al., 2017).Moreover, we outline a possible way to obtain adequate uncertainties for OPC size distributions within the new framework.

OPC measurement principle
The basic principle behind OPC measurements is that particles passing through a sampling volume illuminated by a light source -usually a monochromatic laser -scatter light into a photosensitive detector.The amplitudes of the detected scattering signal pulses are a function of particle size.An OPC counts these pulses and typically sorts them into different bins according to their amplitudes.Therefore the resulting measurement data are histograms of scattering signal amplitudes.The mathematical problem of retrieving number size distributions from recorded scattering signal amplitude histograms is of inverse nature and is described by a set of so-called Fredholm integral equations of the first kind: with the number of particles N i counted in bin i, a term N i accounting for potential counting errors, the corresponding kernel function κ i (D) giving the probability for each particle diameter D to be sorted into bin i and the number size distribution F (D) (Kandlikar and Ramachandran, 1999;Fiebig et al., 2005).Connecting the OPC output, i.e., the particle count histograms, and the desired information, i.e., particle number size distribution, the kernel functions are the key aspect of every OPC measurement.Deriving the kernel functions requires knowledge of the scattering signal amplitude threshold values defining the bin limits, the instrument-specific relationship between scattering signal amplitude and particle scattering cross section and the theoretical relationship between scattering cross section and particle size.The latter is subject to intrinsic particle properties such as complex refractive index and shape.For given intrinsic properties the sizedependent particle scattering cross section C scat (D) with respect to the incident light and OPC scattering geometry, i.e., the solid-angle range covered by the detector, can be calculated.In case of an homogeneous sphere, Mie-Lorenz theory (Mie, 1908) provides an analytical solution.For more complex particle shapes, complementary frameworks like the Tmatrix method (Waterman, 1965) or the discrete dipole approximation (Purcell and Pennypacker, 1973) can be applied.
Bridging the gap between theoretical calculations and the instrument output, i.e., finding the instrument-specific parameters linking C scat (D) with the measured scattering signal amplitude, is the purpose of an OPC calibration.The set of instrument-specific parameters resulting from the calibration in combination with scattering theory allows us to predict the OPC output, i.e., to determine the kernel functions, for any other material with the given optical properties.

Existing concepts for size assignment and calibration evaluation
Though C scat (D) and, hence, the scattering signal amplitude generally are non-monotonic functions of particle size (see Fig. 2), the most popular approach of OPC bin size assignment is to assume or establish monotonicity in order to simplify Eq. ( 1) by allowing for a one-to-one mapping between particle diameter and bin threshold values.One way to achieve monotonicity is to replace the correct C scat (D) by a smoothed monotonic approximation (Cerni, 1983;Osborne et al., 2008).Another option is to simply merge bins in size regions affected by ambiguities in C scat (D), accepting a reduction in resolution (Pinnick et al., 1981).Following these concepts OPC manufacturers usually provide their instruments with a table of predefined (polystyrene latex (PSL)  equivalent) diameter bin threshold values.Mathematically, this means expressing the kernel functions as sharp, adjacent step functions in diameter space: with delta functions at D i and D i+1 , i.e., the lower and upper diameter threshold values of bin i.In doing so, Eq. ( 1) simplifies to and the size distribution can be directly represented by the measured counts in a discrete way as This simplification, however, has fundamental shortcomings: -Even if quasi-monotonicity between particle size and (discretized) scattering signal amplitude can be established for particles of certain intrinsic properties (e.g., polystyrene latex spheres) by a smart choice of OPC collecting optics (Barnard and Harrison, 1988) and/or bin threshold values, this does not automatically hold for particles of different intrinsic properties (e.g., different refractive index or shape) (Szymanski et al., 2009).
-Due to the involved approximations (e.g., a smoothing of C scat (D)) nominal manufacturer values can significantly deviate from reality for certain parts of the instrument size range.Such deviations are regularly reported (Szymanski and Liu, 1986;Rosenberg et al., 2012;Ryder et al., 2013).
-The instrument response can change over time, e.g., due to degradation of OPC light source intensity, pollution or misalignment of optical elements.Such changes usually do not induce a uniform shift in the apparent size distributions but rather cause a complicated deformation.
-No uncertainty estimates are provided for the nominal diameter threshold values.This lack entails an underestimation of size distribution uncertainties.
Some studies stick to the simplified concept of adjacent bins in diameter space but try to reduce possible sizing deviations.Lance et al. (2010) and Cai et al. (2013) use an empirical diameter offset to uniformly shift the manufacturer threshold values in order to yield best agreement between measured histogram modes and nominal diameter values of reference particles.A more universal calibration approach commonly used is to find the parameters for the linear relationship between measured mean scattering signal amplitudes and theoretical mean scattering cross sections for reference particles (Cerni, 1983;Bemer et al., 1990).Still assuming a monotonic C scat (D) they use the resulting linear fit parameters, i.e., slope m and intercept c, to derive a size-dependent scattering signal amplitude U (D) = m • C scat (D) + c and calculate the bin diameter threshold values D i from their predefined scattering signal amplitude counterparts U i .Rosenberg et al. (2012) presented another way of size assignment that avoids workarounds for the non-monotonic A. Walser et al.: On optical particle counter response parametrization and calibration behavior of C scat (D).Their main new concept is to use the same calibration parameters m and c and the unmodified/unsmoothed C scat (D) to define the kernel functions as diameter projections of the scattering signal amplitude bins (see Fig. 2a): C scat,i and C scat,i+1 denote the lower and upper scattering cross-section threshold values of bin i that are a linear function of the actual thresholds given by the scattering signal amplitude values U i and U i+1 .This means that particle diameters will be sorted into bin i if C scat (D) falls within the limits defined by U i and U i+1 scaled with the linear coefficients m and c.This can be further expressed as To simplify the inverse problem of Eq. ( 1) and, again, directly gain size distribution information from OPC histogram data, they use the kernel functions to calculate so-called "perfect" (mean) diameters D p,i and widths W i to characterize all bins: With these values a discrete representation for the size distribution is given by The uncertainties in the calibration parameters m and c are used to derive instrument-related uncertainties for D p,i , W i and, therewith, the resulting size distribution values F D p,i .
Although this approach supersedes workarounds for the ambiguities in C scat (D), it still has shortcomings.One conceptual inadequacy is that representing the bins by their perfect diameter D p,i is ultimately not appropriate, as it will only match with the real mean diameter of particles sorted into bin i in the unrealistic case of a flat size distribution.If, for instance, the size distribution is (strongly) dropping towards larger particles, the occurrence of smaller particle diameters is more likely, meaning that the real mean diameter of particles falling into bin i would be (much) smaller than D p,i .As a result, this causes a sizing bias between the real and calculated size distribution.Spiegel et al. (2012) offer yet another approach to directly estimate size distributions from measured histograms.They translate the range of possible scattering geometries seen by individual particles1 into a range of possible U (D) ∝ C scat (D), a so-called "Mie band".From this Mie band they calculate the a priori probabilities for discrete (equidistant) particle diameter intervals to contribute to the count rate of each bin.This approach potentially allows particle diameters to be sorted into more than one bin, i.e., overlapping OPC bins in diameter space.According to the derived contribution probabilities they then distribute the measured bin counts to the discrete diameter intervals.The major shortcoming of this method is similar to the one discussed above.Even if the a priori probabilities for two diameter intervals (of equal width) to contribute to a certain bin's count rate are the same, the true particle abundance in these intervals will not be the same for the realistic case of a non-flat size distribution.Therefore, the inverse direction, i.e., to equally distribute the counted particles to the two intervals, is generally incorrect.
In summary, none of the existing concepts for OPC calibration and bin size assignment prove completely satisfactory.The simplifications to the inverse problem of Eq. ( 1) and the attempts to directly gain size information from OPC histogram data are always accompanied by systematic errors.Further, some approaches, especially the use of the manufacturer-provided set of bin diameter threshold values, do not offer instrument-related (sizing) uncertainty estimates.

Instrumental broadening of size spectra
A shortcoming common to all available methods is that they do not consider the artificial broadening of size distributions in the basic parametrization of OPC response.One primary cause for the increase in apparent size distribution width is the nonuniformity of light intensity inside the OPC sampling volume (e.g., Wendisch et al., 1996).An inhomogeneity of incident light intensity leads to differences in scattering signal amplitudes for particles passing the sampling volume at slightly varying locations.This means that due to the intrinsic nature of real OPCs, even spherical and homogeneous identical particles appear to vary in size.So far this signal broadening has, if at all, been treated separately from the basic instrument calibration, although it is an instrument-specific property, meaning it is always present in OPC measurements.Depending on the degree of instrument-induced signal broadening in relation to the actual size distribution width, this ef-Figure 2.An example subset of kernel functions for the SkyOPC describing the probabilities for particle diameters and corresponding scattering cross sections to be sorted into the predefined OPC scattering signal amplitude histogram bins, visualized by the different colors.The theoretical relationship between particle diameter and scattering cross section for non-absorbing PSL -with a refractive index of n r = 1.585 at the SkyOPC wavelength of 655 nm -is represented by the black curve.The upper graph (a) shows an ideal case without instrumental broadening of size spectra, whereas the lower graph (b) shows a more realistic case where the effect of signal broadening is considered.The broadening is parametrized by a constant relative Gaussian uncertainty on the scattering cross-section bin threshold values.
fect may lead to significant measurement biases being most pronounced for narrow size distributions (or size distribution modes).In atmospheric research, such narrow size spectra are, for example, met during ice residual measurements in (contrail) cirrus (Voigt et al., 2017) or aerosol chamber experiments (e.g., Schnaiter et al., 2012).
Signal broadening can be further enhanced by other effects such as varying orientation of aspherical particles with respect to the direction of the incident light (Reid et al., 2003) and coincident count events (Baumgardner et al., 1985;Cooper, 1988).The latter becomes relevant for very high particle concentrations when average inter-particle distances are not larger than the size of the sampling volume anymore.In such a case, the probability of erroneously interpreting the sum of several scattering signals from multiple particles as a single particle's signal increases.In addition to an artificial deformation of the size distribution towards larger sizes this entails an underestimation of total particle number concentration.
To correct for artificial broadening of size spectra the common procedure is to define a matrix that contains the probabilities (associated with the broadening effect) to find a particle of a certain size class in adjacent size classes in its elements (Cooper, 1988;Baumgardner and Spowart, 1990;Wendisch et al., 1996;Brenguier et al., 1998).The resulting inverse matrix equation is then solved for the true size distribution.One disadvantage of such methods based on empirical matrices is that their elements might not be universally valid, as for instance the magnitude of broadening that is related to varying particle orientations depends on the degree of particle asphericity.Moreover, the number of uncertaintyafflicted parameters becomes quite large.Assuming an OPC with K bins, the number of parameters required to describe signal broadening is K 2 .
Therefore, for the inversion of OPC histogram data it is advantageous to treat signal broadening in a more universal way.In Sect.2.5.1 we present a new approach that includes the instrument-specific part of signal broadening within the basic parametrization of OPC response.In addition, signal broadening resulting from different orientations of aspherical particles can be included in the inversion process via a set of possible size-to-scattering cross-section relations as outlined in Sect.2.5.3.

Uncertainty in particle properties
So far, we discussed the inverse nature of the OPC measurement principle, challenges and shortcomings of the parametrization of basic OPC response and the artificial broadening of size spectra.An aspect that further complicates OPC measurements is that, in most situations, the (sizedependent) optical properties of the aerosol particles are a priori unknown or at least subject to a considerable degree of uncertainty.Externally or internally mixed individual particles can be combinations of different non-homogeneously distributed materials (e.g., Kandler et al., 2011), making it difficult to find representative complex refractive indices for the bulk aerosol.In any case, the quality of OPC-derived size distributions depends on the quality of information on the optical particle properties.
In order to derive a size distribution uncertainty estimate from uncertainties in the particle properties, however, most studies follow the pragmatic approach and report the maximum impact on the size distribution as a conservative esti-mate (e.g., Osborne et al., 2008;Weinzierl et al., 2009;Brock et al., 2011;Ryder et al., 2013;Hermann et al., 2016).However, the size distribution uncertainty induced by the uncertainties in the particle properties can be substantially size dependent.To yield improved size distribution uncertainty estimates one needs realistic estimates for the set of possible size-to-scattering cross-section relations and a proper way to propagate these estimates (as, for example, outlined in Sect.2.5.3).

New approach
In this section, we introduce an approach to the parametrization of OPC response that involves instrument-specific signal broadening and overcomes the shortcomings of existing methods.We further propose a way to evaluate calibration measurements and to obtain aerosol particle number size distributions with realistic uncertainty estimates from OPC data.

Parametrization of the instrument response
Let a particle of intrinsic properties ϑ (e.g., complex refractive index) and diameter D have the scattering cross section C scat,ϑ with respect to the incident light and OPC scattering geometry.Usually, the scattering signal amplitude in the detector U ϑ is assumed to be a linear function of C scat,ϑ .Hence, for an ideal instrument U ϑ is completely defined by C scat,ϑ and the linear coefficients m and c.However, as explained in Sect.2.3, instrument-induced signal broadening causes the single signal amplitude U ϑ to be replaced by a probability density function (PDF) for a set of possible {U ϑ }.Under the assumption that the nonuniformity of light intensity is a primary reason for the broadening, the relative variance of this PDF is independent of the absolute value of U ϑ .In a simplified approach, the U ϑ PDF is, thus, approximated by a Gaussian distribution with a constant relative standard deviation b.This is equivalent to a "blurring" of the initially sharp OPC scattering cross-section bin threshold values (resulting from the predefined scattering signal amplitude bin thresholds) with the same relative standard deviation b.Replacing the delta functions in Eq. (3) (i.e., Gaussian functions of vanishing standard deviation) by Gaussian functions of constant relative standard deviation yields the new kernel functions: with the new instrument-specific parameter triplet (b, m, c) and the scattering signal amplitude threshold values U i and U i+1 defining bin i. Figure 2 illustrates the difference in OPC kernel functions between an ideal instrument that follows Eq. ( 3) and an instrument with a finite relative Gaussian signal broadening.While the assumption of a Gaussian broadening is an adequate approximation for the OPCs used in this study, it is possible to customize the shape of the broadening effect for other OPC types (e.g., open-path geometries), if necessary.

Calibration evaluation
Taking account of signal broadening, the new parametrization allows for an extension of the classical OPC calibration evaluation approach that is restricted to the determination of the linear coefficients m and c.Given a set of particle standards with known intrinsic properties and size distributions, the forward solution of Eq. (1) using Eq. ( 4) for the kernel functions yields the model count histograms, i.e., the parametrized theoretical instrument response: with the model counts M ij for OPC bin i and particle standard j and the corresponding number size distribution F j (D).With the real measured particle counts N ij the task of a calibration within the new framework is now to inversely find the values for the parameters that bring M ij and N ij into best agreement.For stable measurement conditions, i.e., constant OPC volumetric sample flow, the uncertainties of the measured particle counts follow the Poisson counting statistics.With increasing number of counts, the relative uncertainty hence decreases with N −1/2 ij .Naturally, the simplified model will not be able to reproduce the calibration measurements perfectly because there will be additional deviations that are not parametrized.Provided sufficiently high numbers of counts in the course of the sampling, the relative bin count uncertainties due to Poisson counting statistics will become negligible compared to these additional deviations.As a consequence, bringing model and measurement into agreement corresponds to maximizing the probability of the model counts M ij afflicted with a priori unknown uncertainties σ ij , which cover the additional model deviations, to occur given the measured N ij .To ensure that the modeled instrument response later agrees with reality within its margin of uncertainty, it is necessary to find a good representation of the unknown uncertainties σ ij and quantify them in the course of the calibration, too.
Comparing measured particle standard histograms with the corresponding model results for a suitable instrument parameter tuple (b, m, c) reveals that remaining deviations be-tween the two mainly appear as (nonuniform) small shifts, meaning that compared to the model histograms some measured histograms are shifted to smaller while others are shifted to larger scattering signal amplitudes.It thus seems natural to treat these deviations as a remaining uncertainty of the modeled scattering signal amplitudes.Apart from the experimental finding, there are also theoretical explanations for the observed shifts.OPC light source intensity fluctuations around a temporal average induce time-dependent scattering signal amplitude variations, causing measured histograms to move up and down slightly with time.In contrast to the instantaneous signal broadening discussed in Sect.2.3, these fluctuations act on greater timescales, leading to histograms shifts that are differently pronounced for samples recorded with certain time lags.Other possible sources for such timedependent shifts are changes in detector sensitivity or signal background noise.Due to the linear relationship, timedependent signal amplitude fluctuations can equivalently be thought of as relative fluctuations in the scattering cross sections, assuming a fixed light source intensity.Therefore, we express the model count uncertainties σ ij in terms of a relative uncertainty of the theoretical particle scattering cross sections C scat,ϑ (D).
For a given instrument parameter tuple (b, m, c) the set of model bin counts M ij results from Eqs. ( 5) and ( 4) with the (best estimate) theoretical function C scat,ϑ (D).A relative shift in the theoretical scattering cross sections corresponding to a multiplication of C scat,ϑ (D) by a factor ε = 1 leads to a different set of model bin counts M ij,ε .Assuming the PDF of the possible relative shifts ε to be a Gaussian function centered at 1 and having a standard deviation of σ ε one can derive the respective PDFs for the model bin counts M ij,ε .For the sake of convenience and simplicity the resulting model bin count PDFs can themselves be approximated by Gaussian PDFs, which is usually an adequate approximation.This leads to the following expression for the unknown model bin count uncertainties: with M ij,ε defined by Eqs. ( 5) and ( 4), replacing D).In summary, the new calibration evaluation should yield the set of model parameters (b, m, c, σ ε ) composed of the instrument-specific parameter tuple (b, m, c) and, according to the above considerations, the remaining relative uncertainty of the theoretical scattering cross sections σ ε .
A way to meet the challenge of model parameter probability maximization under initially unknown model uncertainties is to make use of Bayesian statistics and Markov chain Monte Carlo (MCMC) methods (e.g., Goodman and Weare, 2010).Following Bayes' theorem (Bayes and Price, 1763) the (posterior) probability P for a set of model bin counts M ij to occur under a set of measured bin counts N ij can A. Walser et al.: On optical particle counter response parametrization and calibration be expressed as i.e., the product of the likelihood function determining the probability of the N ij to occur given the M ij and the so-called prior probability P b, m, c, σ ij , including all prior knowledge on the model parameters for instance from physical constraints or invariance considerations (e.g., Jaynes, 1968).The proportionality factor equating both sides of Eq. ( 7) can be thought of as a normalization constant.Upon the assumption of Gaussian model bin count PDFs the likelihood function can be expressed as with M ij and σ ij defined by Eqs. ( 5) and ( 6) respectively.MCMC methods allow us to efficiently sample the model parameter space utilizing the forward solution to the problem to find the region of maximum probability according to Eq. ( 7).This way, the PDFs for the instrument parameters (b, m, c) and the relative uncertainty of the theoretical particle scattering cross sections σ ε are obtained together with all correlations between the individual parameters.In this study we utilize the Python-based sampler tool emcee (Foreman-Mackey et al., 2013).

Retrieval of size distributions within the new framework
The new instrument parametrization, including instrumentspecific signal broadening and the parameter PDFs resulting from the MCMC-based calibration evaluation, now permits us to derive size distributions from OPC measurements in a self-consistent way.Propagating the parameter uncertainties yields improved estimates for the corresponding size distribution uncertainties.Figure 3 illustrates a possible workflow within the proposed framework to go from measured OPC count histogram data to PDFs in size distribution solutions.Similar to what has been proposed by Fiebig et al. (2005), the basic idea is to start with random Monte Carlo samples drawn from the model parameter PDFs and a set of possible theoretical particle diameter to scattering crosssection relationships C scat,ϑ (D) , e.g., given by the likely range of aerosol particle complex refractive indices, their shape and orientation.2In addition, a random relative shift from the chosen C scat,ϑ (D) is picked according to its potential (time-dependent) systematic deviation, e.g., induced by light source intensity fluctuations.This relative shift is drawn from the Gaussian PDF N 1, σ 2 ε parametrized by σ ε , which is derived as part of the calibration evaluation.With the resulting (shifted) C scat,ϑ,ε (D) and the instrument parameter tuple (b, m, c), the set of OPC kernel functions {κ i (D)} can be calculated following Eq.(4).Given this set of bin kernel functions the aerosol size distribution F (D) is adjusted such that the (uncertainty-weighted) deviations between modeled and measured bin counts are minimized.The result is then either one best solution for F (D) or an ensemble of possible solutions for each iteration, depending on the respective inversion algorithm (see, e.g., Kandlikar and Ramachandran, 1999;Fiebig et al., 2005, and the references herein).By repeating this procedure multiple times one finally acquires a collective solution ensemble, representing members of the size distribution solution PDF, considering all uncertainties in instrument-specific parameters and the theoretical diameter to scattering cross-section relationships.In this work (see Sect. 4) we use a parametrized size distribution and, again, a MCMC method for the inversion.The corresponding size distribution parameter solutions obtained for each Monte Carlo iteration are merged into a final size distribution parameter solution ensemble.
Apart from a thorough and transparent derivation of size distribution uncertainties, the proposed retrieval method has further advantages.For one, the Monte Carlo sampling enables a one-to-one mapping between each size distribution solution (ensemble member) and the corresponding initial parameter picks, thereby facilitating, for example, parameter sensitivity studies.This can help to identify dominating initial parameter uncertainties and even allow us to confine the initial parameter estimates by comparing the traceable solutions to size distribution results obtained by independent measurements.It should be clarified, however, that the proposed retrieval method alone simply propagates the initial parameter PDFs and cannot provide information on their adequacy; i.e., it cannot judge the value of individual parameter picks.A big advantage of the method is that the solution ensemble itself allows for a simple yet appropriate further propagation of size distribution uncertainties.This is explicitly useful when calculating quantities depending on the size distribution.For instance, PDFs for the effective particle diameter or the aerosol extinction coefficient can easily be (numerically) derived by collecting the individual results obtained for each size distribution solution (ensemble member).In light of all benefits, we recommend using the proposed retrieval method (or equivalent approaches) for every occasion, despite the additional effort involved.3 Experimental setup

Involved OPCs
The two central OPCs examined in this study are the Grimm model 1.129 (SkyOPC) and the DMT Airborne Passive Cavity Aerosol Spectrometer Probe (PCASP-100X with an upgraded signal processing package SPP-200, abbreviated PCASP hereafter).Both aerosol spectrometers were part of the airborne in situ instrumentation used in the SALTRACE campaign.The SkyOPCs were operated inside the cabin of the German Aerospace Center's Falcon research aircraft behind an isokinetic aerosol inlet, while the PCASP was mounted in one of the under-wing stations.Detailed descriptions of the instruments can be found in Bundke et al. (2015) for the SkyOPC and in Liu et al. (1992) and Strapp et al. (1992) for the PCASP.Both are closed-path spectrometers in which the aerosol particle beam is confined to the inner area of light source focus.The particles scatter light coming from a monochromatic laser of visible red wavelength (633 nm helium-neon laser for the PCASP and 655 nm diode laser for the SkyOPC), which is then detected in a sideways direction.The applied wide-angle collection of scattered light minimizes ambiguities in C scat (D) (see Heim et al., 2008, andRosenberg et al., 2012, for details on the scattering geometry).
During SALTRACE and the lab measurements presented here the SkyOPC was operated in the fast mode for smaller sizes, covering a nominal diameter range of 0.25 to about 3 µm.The corresponding scattering signal amplitude range is A. Walser et al.: On optical particle counter response parametrization and calibration separated into 16 preset bins defined by a set of digital threshold values.In standard configuration, the PCASP sorts scattering pulses into 30 bins over a nominal diameter range of 0.1 to 3 µm.It further allows for a custom selection of the digital bin threshold values.In this study, we only consider the PCASP low gain stage. 3In order to better study differences between the approaches discussed in Sect.2, we present results for a custom high-resolution binning of this gain stage in addition to the default binning.For the custom binning the gain stage's signal amplitude range is divided into bins with constant width.
The DMT Ultra-High Sensitivity Aerosol Spectrometer (UHSAS) (Cai et al., 2008) lab version covering a size range of about 0.06 to 1 µm in high resolution (99 bins) was utilized as a reference for total particle concentration during the SkyOPC calibration measurements and further served for a qualitative assessment of sizing to support the particle mobility filtering described in Sect.3.2.

Measurements
The calibration measurements were performed using monodisperse aerosols of PSL spheres.For the SkyOPC, the data set is complemented by di(2-ethylhexyl) sebacate (DEHS) aerosol samples.The complex refractive indices for PSL and DEHS are approximately 1.585 + i0 (Sultanova et al., 2009) and 1.45 + i0 (manufacturer data sheet) in the wavelength range of the SkyOPC and PCASP.PSL spheres dispersed in distilled water were mobilized via nebulization with the DMT portable aerosol generator running with aerosol-free carrier air.The resulting aerosol was subsequently dehumidified by an arrangement of silica gel dryers and diluted to avoid OPC measurement issues related to coincident count events.The DEHS aerosol was produced using a TSI model 3475 condensation aerosol generator based on the Sinclair-LaMer principle (Altmann and Peters, 1992): with nitrogen as the carrier gas an aqueous sodium chloride solution is nebulized and dried to yield a high-concentration condensation nuclei aerosol.Passing through a heated vessel filled with liquid DEHS, a reheater unit and a condensation chimney, the precursor particles allow for heterogeneous condensation of the supersaturated DEHS vapor.The mean particle size of the resulting (quasi-)monodisperse aerosol is a function of the ratio between vapor concentration and condensation nuclei number concentration.Again, the DEHS aerosol was diluted prior to the measurement in the OPC to circumvent counting coincidences.For all measurements we chose sampling interval times long enough to minimize relative uncertainties from counting statistics.
For mean particle diameters up to 800 nm the aerosol was additionally filtered with the aid of a differential mobility analyzer (Grimm Vienna-type L-DMA, abbreviated DMA hereafter; Reischl et al., 1997).DMA filtering substantially reduced the widths of the DEHS aerosol size distributions and allows us to obtain quantitative information on their mean diameters and widths.Using the UHSAS, the DMA transfer function was carefully centered to the middle of the initial DEHS generator aerosol size distribution to guarantee that the resulting size distribution is in good approximation, represented by the Gaussian DMA transfer function itself.The relative standard deviations of the DMA transfer functions are calculated by means of the formulae given in Reischl et al. (1997) and Stolzenburg (1988).The initial PSL particle size distributions, which are traceable via the United States National Institute of Standards and Technology (NIST), are of Gaussian shape with known mean and standard deviation.Although they are narrower than the width of the DMA transfer functions, additional DMA filtering helped to effectively remove the interfering background at smaller particle diameters than are caused by the nebulization (see Hermann et al., 2016).For mean particle diameters larger than 800 nm the presented counting histograms are empirically separated from this background.

Results and discussion
In this section, we present the results for the evaluation of the PSL calibration measurements following the new method proposed in Sect.2.5.We compare these results with the theoretical instrument response for nominal manufacturer diameter bin thresholds and results obtained for the approach of Rosenberg et al. (2012), representing a state-of-the-art conventional method.Hereafter, we abbreviate these approaches as MFR and R12, respectively.We further demonstrate the impact of method choice on size distribution inversion results for measurements of DEHS samples.
The measurements of PSL particle standards, carried out as described in Sect.3.2, are utilized to calibrate the OPCs following both the new and the R12 approach (introduced in Sect.2.2).Figures 4 and 5 contrast the resulting modeled relative bin count histograms and the measured relative histograms for the SkyOPC and the PCASP (low gain stage) respectively.The model histograms are calculated by means of Eq. ( 5) with the well-defined Gaussian PSL size distributions and the kernel functions given by Eq. ( 3) for the R12 (shown in red brown colors) and Eq. ( 4) for the new approach (shown in blue colors).The best estimate model histograms, i.e., the model histograms for the maximum probability model parameter tuple -(m, c) best for the R12 and (b, m, c, σ ε ) best for the new approach -are represented by the color-framed white histogram bars.For the SkyOPC, additionally the model histograms for the MFR approach following Eq.( 2) and using the manufacturer-supplied set of nominal values are displayed in golden colors.The underlying measured histograms are depicted by the gray bars.For the new and the R12 approach the parameter PDFs resulting from the evaluation of the calibration measurements (see Figs. A1 and A2) are sampled using a Monte Carlo method to yield the corresponding PDFs of the model histogram bin counts that are visualized by error bars spanning the range between the 16 and 84th percentiles.In each panel of Figs. 4 and 5 the mean diameter and standard deviation of the Gaussian PSL size distribution is displayed in the left upper corner.Figure 6 supplements the SkyOPC histogram comparisons with scatter plots showing all modeled and measured relative bin counts for the different approaches.Finally, Fig. 7 quantitatively compares the total sum of residuals, between the measured N and best estimate model relative bin counts M for the two instruments and the different approaches.The subscripts i and j represent the different OPC bins and used particle standards respectively.
The model histograms for the MFR approach (e.g., Fig. 4, golden colors) exhibit significant deviations from the underlying measured histograms.They offer much smaller widths than their measured counterparts.In addition, absolute offsets between the histogram modes are apparent for both Sky-OPC and PCASP (not shown).Deviations are largest for the SkyOPC because it was operating under dusty conditions during SALTRACE over a longer period previous to the presented measurements, presumably causing a pollution of optical elements.In consequence, the scatter plots for the MFR approach in the upper row of Fig. 6 show the largest discrepancy between model and measurements.This becomes also obvious for both instruments when looking at the total sums of residuals in Fig. 7.The residuals for the MFR approach are substantially enhanced compared to the others.
The R12 approach allows for the correction of the absolute shifts of the histogram modes.Nevertheless, instrumentspecific signal broadening is still ignored.The modeled histograms, thus, continue to underestimate the widths of the acwww.atmos-meas-tech.net/10/4341/2017/Atmos.Meas.Tech., 10, 4341-4361, 2017 Figure 5. Same as Fig. 4 but for the PCASP and a different set of PSL particle standards.Instead of the PCASP (low gain stage) default binning, a custom high-resolution linear partitioning is applied here to better highlight the differences between the approaches.
tually measured histograms, which is visible in the histogram plots in Figs. 4 and 5. Here, and especially in Fig. 6, it is also apparent that the R12 approach remains unable to reproduce the measurements within the margins of model uncertainty for most of the relative bin counts.Particularly for the smaller relative count values the absence of a parametrization of signal broadening leads to large model deviations.However, in comparison to the MFR approach total residuals for the model best estimates are reduced by 25 and 35 % for the SkyOPC and PCASP (low gain stage, default binning) respectively.Beyond that, an estimate for the model uncertainty is established.By introducing a simple parametrization of instrumentspecific signal broadening and a self-consistent way of evaluating OPC calibration measurements, the new method succeeds in modeling the measured histogram widths correctly (see Figs. 4 and 5 rightmost columns).As a result, the total residuals between measured and modeled relative bin counts for the model's best estimates decrease by 82 and 77 % compared to the MFR approach for the SkyOPC and PCASP (low gain stage, default binning) respectively.With respect to the R12 approach total residuals for the SkyOPC, the PCASP default binning and the finer PCASP custom binning are lowered by 77, 64 and 76 %.Further, Fig. 6 shows that the new approach proves capable to correctly reproduce the measured histograms within the margins of model uncertainty over the complete range of relative bin counts.
Figure 8 shows the SkyOPC counting efficiency curve obtained by parallel measurements with the UHSAS as a reference counter during the PSL calibration measurements.These measurements offer another perspective on the comparison between the two approaches.The mean total concen- are calculated from the respective total number of counts and the volumetric instrument flow rates ϕ for each particle standard j and are depicted by the red diamond markers.The associated 68 % confidence intervals (approximately corresponding to ± 1 standard deviation) that result from error propagation involving count rate scatter and instrument sample flow uncertainties are represented by the red error bars.
The modeled concentration fractions are derived from the bin kernel functions κ i as and are visualized by the solid lines for the model best estimates, again in red brown for the R12 and in blue for the new approach.The shaded areas show the range between the 16 and 84th percentiles derived from the model parameter PDFs.The R12 approach predicts a sharp drop-off to smaller particle diameters in contrast to the measurements.The new approach is able to correctly model both shape and absolute values of the observed sigmoidal behavior of the counting efficiency curve.Measurements of DEHS samples, as outlined in Sect.3.2, allow us to test the possible implication of the choice of method for size distribution inversion results using an independent material.As proposed in Sect.2.5.3 and illustrated in Fig. 3, the inversion of measured OPC histogram data is based on the parametrization of instrument response, the respective parameter PDFs derived from the calibration and C scat,ϑ (D) for the new material.The use of DEHS spherical droplets guarantees that this latter relationship is welldefined for the given scattering geometry as complex refractive index and shape of the aerosol particles are known, thus adding no further complexity to the retrieval.Moreover, the size distribution of the filtered DEHS samples approximately follows a Gaussian distribution, simplifying the inversion in this case to the determination of the size distribution parameters, mean diameter µ sd and standard deviation σ sd .The inversion algorithm used here to solve Eq. (1) for the parametrized size distribution is based on a MCMC method (Goodman and Weare, 2010;Foreman-Mackey et al., 2013).To obtain adequate size distribution parameter PDFs, 10 000 Monte Carlo samples are drawn from the corresponding instrument parameter PDFs. Figure 9 shows the inversion re-  sults for two DEHS samples and the two methods, i.e., the R12 approach (red brown) and the new one (blue).The theoretical (true) size distributions and the corresponding values for the means and standard deviations are depicted by the red lines and markers.For both the new and R12 approach the retrieved size distribution means agree with the theoretical values within their range of uncertainty, meaning that both methods allow for a correct (mean) sizing.This finding additionally proves the validity of the used C scat,ϑ (D) (for DEHS and the calibration material PSL).Upon closer inspection the retrieved means tend to slightly underestimate the true values, which could imply minor deviations between the true and the used OPC scattering geometry and/or refractive index values for PSL and DEHS.The parameter PDFs for the size distribution means µ sd are almost identical for the two methods concerning both PDF median values and widths, i.e., uncertainty ranges.This agreement disappears for the size distribution standard deviations σ sd .The new method again agrees with the theoretical values within the range of parameter uncertainty and, hence, successfully predicts the full shape of the size distribution.The R12 approach attributes the width of a measured histogram completely to the width of the size distribution, thus overestimating σ sd significantly.For the examples shown here, the R12 approach overestimates the true values for σ sd by 714 and 302 % with respect to the medians of the retrieved parameter PDFs.The widths of the σ sd parameter PDFs, i.e., the estimated range of uncertainty in this parameter, also differ for the two methods.With respect to the distance between the 16 and 84th percentiles the R12 approach yields 285 and 224 % higher PDF widths than the new method leading to greater overall uncertainties in the retrieved size distributions, which are nonetheless unable to encompass the true ones.
It should be noted, though, that the standard deviations of the DEHS size distributions used here are quite small.When size distributions become broader the impact of instrumentspecific signal broadening on the width of the recorded histograms decreases and, hence, differences between the methods will become less pronounced.Besides, uncertainties in aerosol properties like complex refractive index and shape might be the dominant source of size distribution uncertainty in many situations.However, this example demonstrates that the new method is able to retrieve even narrow size distributions correctly and, hence, to provide access to realistic uncertainty estimates for all situations.The results also imply that even for the same data and OPC instrument, calibrated with the same set of measurements, retrieved size distributions can be contradictory solely due to different instrument response parametrizations and calibration evaluation approaches.

Conclusions
Retrieving aerosol particle number size distributions and associated uncertainties from OPC histogram data is a challenging task.Scattered light intensity (the measurand) generally is a non-monotonic function of particle size (the quantity of interest) and depends also on particle intrinsic properties such as complex refractive index.Besides, due to the non-ideal behavior of real OPCs, measured intensity distri- butions are artificially broadened.To realistically model OPC response, i.e., to find suitable OPC bin kernel functions defining the probabilities for particle diameters to be sorted into the instrument's discrete scattering signal amplitude bins, is thus a crucial requirement.
We have introduced a new approach to model OPC response and, within this framework, a self-consistent way for the evaluation of calibration measurements.Two OPCs involved in the SALTRACE campaign, the SkyOPC and the PCASP, and measurements of PSL particles have been utilized to compare the new approach with existing concepts.The results lead to the following conclusions.
The manufacturer-provided set of (PSL-equivalent) nominal diameter threshold values for the OPC bin borders should be treated with caution and the resultant size distributions should be considered as rather qualitative measures.Not only can the concept of adjacent continuous bins in diameter space be problematic given the non-monotonic relation between particle size and scattering signal amplitude, but the values are also material-dependent and drifts in size assignment, e.g., due to pollution of OPC optics or light source intensity drifts, can occur over time.We have shown that the corresponding size distributions can significantly deviate from reality, even for the reference material.Furthermore, no uncertainty estimates are provided for the nominal diameter values that could be used to infer instrument-related size distribution uncertainties.
Calibrating the instrument can remove absolute sizing offsets.The results for a state-of-the-art OPC calibration and response parametrization approach (Rosenberg et al., 2012) exhibit clear improvements in sizing and, therewith, a reduction in total residuals between modeled and measured bin histograms.The introduction of instrument parameter uncertainties that go along with the calibration evaluation alwww.atmos-meas-tech.net/10/4341/2017/Atmos.Meas.Tech., 10, 4341-4361, 2017 lows us to derive related size distribution uncertainty estimates.However, these estimates fail to explain the remaining differences between modeled and measured instrument response for the presented data.The main reason for this is the absence of a parametrization of instrument-induced signal broadening.This artificial increase in apparent size distribution width, which is stronger the narrower the actual size distribution is compared to the degree of broadening, may involve significant systematic OPC measurement biases (in atmospheric research) when disregarded.By introducing a simple (one parameter) approach to describe this ever-present broadening of size spectra, the new method leads to substantial improvements.Residuals between modeled and measured OPC response are considerably reduced compared to the other methods.The new method further correctly predicts the size dependency of OPC counting efficiency.Most importantly, the measurements are successfully reproduced within the range of model uncertainty.
In the context of the new method we have also outlined a self-consistent way to thoroughly propagate parameter uncertainties and gain realistic size distribution PDFs without avoiding to address the actual inverse problem underlying OPC measurements.Besides the advanced uncertainty assessment, a benefit of the proposed Monte Carlo retrieval procedure is the facilitation of subsequent uncertainty propagation for quantities calculated from the size distribution (e.g., the effective diameter).When this procedure is combined with the new OPC response model, exemplary results for measurements of DEHS samples demonstrate that even narrow size distributions are retrieved correctly.For the conventional method the same retrieval procedure, propagating the corresponding parameter uncertainties, yields larger size distribution uncertainties and significantly overestimated size distribution widths.
In summary, the new method has the following major advantages over existing concepts for OPC bin size assignment: -The inevitable instrument-specific broadening of measured size spectra is parametrized for the first time, leading to a more accurate modeling of OPC response.
-The model parameter PDFs resulting from the evaluation of calibration measurements allow for realistic uncertainty estimates for this response and, as a consequence, provide a basis for proper size distribution uncertainties.
Figure A1 and A2 show the calibration results for the PCASP low gain stage with the custom high-resolution binning following the new and the R12 approach, respectively.As explained in detail in Sect.2.5, the parameter solution ensemble for the new approach is obtained with the aid of the Pythonbased MCMC sampler tool emcee (Foreman-Mackey et al., 2013), using Eq. ( 8) for the likelihood function.
For the R12 approach the theoretical particle scattering cross-section means C scatj and standard deviations σ C scat j for each particle standard j are calculated from the known PSL size distributions and the theoretical C scat (D) following Eqs.( 4) and ( 5) in Rosenberg et al. (2012).Accordingly, for the scattering signal amplitudes the measured histogram mode bin centers U j and half widths (σ U ) j are calculated.Further adapting their procedure, for the linear fit between the two properties the scattering cross-section and signal amplitude PDFs for each particle standard are approximated by uncorrelated Gaussian distributions with the aforementioned values defining the respective means and standard deviations.In order to yield a parameter solution ensemble similar to the new approach, in this study the fitting is conducted by means of the same MCMC tool using the following logarithmic likelihood function for a linear relation between two properties with uncorrelated Gaussian uncertainties: where m and c represent the fit slope and intercept, K is a (neglectable) constant offset that has no influence on the maximization process and the resulting parameter solution PDFs.
The two approaches' results for m and c are consistent with respect to the corresponding uncertainty ranges (see Figs. A1 and A2) with the new method yielding smaller uncertainty ranges for both m and c.The median values for the new method's additional parameters b and σ ε , describing the relative standard deviation ("blurring") of the scattering cross-section bin threshold values and the remaining relative uncertainty of C scat (D), approximately amount to 22 and 10 %, respectively.The slight cutting of larger values in the PDF for c resulting for the R12 approach is caused by the physical constraint that the scattering cross-section values for the bin thresholds may not be negative.

Figure 1 .
Figure1.Sources of uncertainty for size distributions derived from OPC measurements.The main purpose of this study is to introduce an advanced description of OPC response and to offer improved estimates for the corresponding uncertainties (red-rimmed box).

Figure 3 .
Figure 3. Flow chart demonstrating a possible pathway for the retrieval of size distribution information from OPC histogram data within the new framework.

Figure 4 .
Figure 4. Comparison of modeled relative histograms (colored) and measured counterparts (gray, hatched) for the SkyOPC and different PSL particle standards (rows) for different approaches of OPC kernel function parametrization (columns).The colored histogram bars represent each model's best estimate and the error bars are the range between the 16 and 84th percentiles of the corresponding PDFs.Panel (a) shows the theoretical instrument response according to the manufacturer-provided set of nominal diameter threshold values in gold (MFR), panel (b) the results following the calibration and instrument parametrization approach by Rosenberg et al. (2012) in red brown (R12) and panel (c) the results of the new approach in blue.

Figure 6 .
Figure 6.Scatter plots of all modeled relative SkyOPC bin counts for the PSL standards versus their measured counterparts for the three different approaches (rows).The comparisons are shown on linear and logarithmic scales on the left-and right-hand side respectively.The markers represent the model best estimates and the error bars are the range between the 16 and 84th percentiles of the corresponding PDFs.The black lines follow the one-to-one relationship.Significant model underestimations, i.e., vanishingly small model values where nonvanishing bin counts are measured, occur in the two upper rows.The number fraction of significantly underestimated values is noted in the upper left corner of the logarithmic scale plots and the corresponding values are shown with triangular markers in the linear scale plots.

Figure 7 .
Figure 7.Total sum of residuals between measured relative bin counts and the corresponding model best estimates including all PSL calibration measurements for the SkyOPC and PCASP.In addition to the absolute residual values (solid bars), the arrows and percentage numbers demonstrate the relative reduction by changing the approach.

Figure 8 .
Figure 8.Comparison between modeled and measured SkyOPC (bins 1-15) counting efficiency.The measured mean counting efficiency values are plotted with red diamond markers and their associated 68 % confidence intervals with red error bars.The solid lines represent the model best estimates for the different approaches.The shaded areas correspond to the range between the 16 and 84th percentiles.

Figure 9 .
Figure 9. Parametrized size distribution retrieval results for two DEHS samples with mean diameters of 0.4 (upper row, graphs a1 to a3) and 0.5 µm (lower row, graphs b1 to b3).The (normalized) Gaussian size distributions are shown in graphs (a1) and (b1).For both the R12 and the new approach the size distribution retrieval PDFs are represented by their diameter-wise medians (solid lines) and 2nd, 16, 84 and 98th percentiles (shaded areas).The theoretical (true) size distributions are indicated by the red lines.The corner plots display the solution PDFs for the size distribution parameters, i.e., mean and standard deviation for the R12 approach in graphs (a2) and (b2) and the new approach in graphs (a3) and (b3).The dashed lines in the 1-D histograms represent the parameter PDF medians, 16 and 84th percentiles (in µm).The median values and their distances to the percentiles are noted on top of each histogram.The 2-D correlation plots show the solution scatter (black points) superposed with color-coded 2-D histograms and smoothed Gaussian contours at 0.5, 1.0, 1.5 and 2σ .The true parameter values are again indicated by the red lines and markers.
Figure A1.PCASP (low gain stage) calibration results following the new method.The model parameter solution ensemble is visualized in the same way as in Fig. 9.

Figure A2 .
Figure A2.PCASP (low gain stage) calibration results for the approach of Rosenberg et al. (2012).The black markers in panel (a) represent the centers of the scattering signal amplitude histogram modes measured for the PSL standards and the corresponding calculated mean scattering cross sections.The error bars represent the histogram modes' half widths and scattering cross-section standard deviations respectively.The red brown solid line shows the best fit and the shaded area the range between the 16 and 84th percentiles of the fit function PDFs.To allow for direct comparison with the results of the new method, panel (b) shows the parameter solution ensemble in the same way as in Fig. A1.