Inversion of electrical mobility measurements using bipolar or unipolar chargers for the arbitrary distribution of channels

Please cite this article in press as: Domat, M., et al. Inversion of electrical mobility measurements using bipolar or unipolar chargers for the arbitrary distribution of channels. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2014.08.007 ARTICLE IN PRESS G Model PARTIC 738 1–10 Particuology xxx (2014) xxx–xxx Contents lists available at ScienceDirect Particuology journal homepage: www.elsevier.com/locate/partic Inversion of electrical mobility measurements using bipolar or unipolar chargers for the arbitrary distribution of channels 1 2 M. Domata,∗,F.E. Kruisb,N.L. Azong-Warab,J.M. Fernandez-Diaza Q1 3 aDepartment of Physics, University of Oviedo, C/ Calvo Sotelo, s/n, E-33007 Oviedo, Spain4 bInstitute of Technology for Nanostructures (NST) and Center for Nanointegration Duisburg-Essen (CENIDE), University of Duisburg-Essen, Bismarckstr. 81, D-47057 Duisburg, Germany 5 6 7 article info 8 9 Article history: 10 Received 16 March 201411 Received in revised form 1 August 201412 Accepted 15 August 201413 14 Keywords:15 Tandem differential mobility analysis16 Inverse problem17 Regularization methods18 Size distribution19 Transfer function20 Unipolar charger21 abstract The inversion of the particle size distribution from electrical mobility measurements is analyzed. Three different methods are adapted for a dot-matrix approach to the problem, especially for non-square or singular matrices, and applied to electrical mobility measurements from fixed or scanning voltages. Multiply charged particles, diffusion losses, arbitrary voltage steps, and noise were considered, which results in non-adjoining and overlapping transfer functions. The individual contribution of the transfer functions in each size interval was geometrically estimated, which requires only its characteristic mobilities. The methodology is applied to mobility measurements from particles charged with unipolar and bipolar chargers. However, the method can be extrapolated to any charging method with a defined charge distribution, and retrieval of the singly charged particle distribution and mean charge from a tandem differential mobility analysis configuration was successfully demonstrated. © 2014 Published by Elsevier B.V. on behalf of Chinese Society of Particuology and Institute of Process Engineering, Chinese Academy of Sciences. 22 Introduction23 Q3 Knowledge of the particle size distribution (PSD) is a funda-24 mental requisite to characterize an aerosol population. Differential25 mobility analysis (DMA) is the most common principle to measure26 submicron PSDs, because it is able to characterize a wide range of27 sizes with high resolution in almost real time independently of the28 composition of the particles.29 However, the PSD cannot be directly measured and must be30 inferred from related mobility measurements using inversion tech31 niques. The first attempt at inversion was by Knutson and Whitby 32 (1975a, 1975b) based on integration of the particle trajectory equa-33 tion inside the DMA. The particle size was represented by the size34 Abbreviations: CN, condition number; CPC, condensation particle counter; DMA, differential mobility analyzer; DMPS, differential mobility particle sizer; FWHM, full width at half maximum; GSD, geometric standard deviation; GCV, generalized cross-validation; LC, L-Curve; LS, least squares; NNLS, non-negative least squares; NRMSE, normalized root mean square error; PSD, particle size distribution; SMPS, scanning mobility particle sizer; SVD, singular value decomposition; TDMA, tandem differential mobility analysis; TF, transfer function. ∗Corresponding author. Tel.: +34 659390180. Q2 E-mail address: [email protected] (M. Domat). of a sphere with the same mobility Zpas the particle, which is called 35 the mobility diameter (DZp): 36 Zp=qeCc(DZp) 3DZp ,(1) 37 where qis the number of elementary charges eon the particle, is 38 the dynamic viscosity of the gas flow, and Cc(DZp) is the Cunning39 ham slip correction factor. 40 The particle size distribution N(Dp) and particle mobility dis41 tribution M(DZp) are related through the non-negative Kernel 42 function Kq(Dp,D Zp) by a Fredholm-type integral equation that also 43 considers the uncertainty in the measurements ε(DZp): 44 M(DZp)=Dp max Dp min Kq(Dp,D Zp)N(Dp)dDp+ε(DZp).(2) 45 When the Kernel is linear with respect to the size distribution, 46 Eq. (2) can be expressed in matrix form as 47 Mm≈Kq mnNn+εm.(3) 48 The equation is approximate because of the error εm∈Rmnot 49 only in the data, but also in the approach of the integral (Eq. (2))toa 50 set of discrete linear equations. The indices mand nare the number 51 of discrete measurements and output channels, respectively. 52 http://dx.doi.org/10.1016/j.partic.2014.08.007 1674-2001/© 2014 Published by Elsevier B.V. on behalf of Chinese Society of Particuology and Institute of Process Engineering, Chinese Academy of Sciences. Please cite this article in press as: Domat, M., et al. Inversion of electrical mobility measurements using bipolar or unipolar chargers for the arbitrary distribution of channels. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2014.08.007 ARTICLE IN PRESS G Model PARTIC 738 1–10 2M. Domat et al. / Particuology xxx (2014) xxx–xxx Nomenclature Aarea of the transfer function CcCunningham slip correction factor dratio between diameters Dgap vertical distance from needle tip to output of ionizing region in corona Dpparticle diameter Dpg geometric mean diameter DZpmobility diameter eelementary unit of charge (1.602 ×10−19)C fcfraction of charged particles Hheight of the fitting TF triangle kBBoltzmann’s constant (1.380 ×10−23) J/K Kkernel matrix Llength of DMA mnumber of measurement intervals Mparticle mobility distribution nnumber of output intervals Nparticle size distribution Ppressure Pnpenetration efficiency Pe Peclet number qnumber of elementary charges on a particle Qiion dilution flow Qaaerosol flow Qex excess flow Qmaerosol monodisperse flow Qsh sheath flow Ri,R0inner and outer radius of the DMA electrodes tresidence time V0voltage applied to the DMA central rod Zccentroid mobility of transfer function Zl,Zulower and upper boundaries from the mobility TF intervals Z∗ l,Z∗ ulower and upper boundaries from the size intervals Zpelectrical particle mobility Greek letters ˇ0relative width of the TF εerror regularization parameter dynamic viscosity of gases ggeometric standard deviation scan time transfer function (math.) One of the main drawbacks of this technique is that large par-53 ticles with multiple charges can have the same mobility as small54 singly charged particles. It is an ill-posed problem in the sense that55 a stable solution may not exist or may not be unique.56 In general, the rank of the control matrix Kq mn determines the57 reliability of finding a solution. For determined systems (m=n), 58 there may be a unique solution. If the system is under-determined59 (m<n), it has many possible solutions, whereas if m>nit is an60 overdetermined system and there is no exact solution (Voutilainen,61 2001; Talukdar & Swihart, 2003).62 Data inversion technique63 The ill-posedness of the problem can be rectified by replac-64 ing it with an approximate well-posed problem whose solution is65 assumed to be close to the actual PSD. The main methodology to66 solve Eq. (3) is based on methods derived from the least squares 67 (LS): 68 NLS =arg min KNinv −M .(4) 69 In this work, only deterministic techniques were applied, such 70 as iterative (Alofs & Balakumar, 1982; Bazan & Francisco, 2009; 71 Collins, Flagan, & Seinfeld, 2002; Crump & Seinfeld, 1981; Fiebig, 72 Stein, Schröder, Feldpausch, & Petzold, 2005; Hagen & Alofs, 1983; 73 Pfeifer et al., 2013; Rojas & Steihaug, 2002; Talukdar & Swihart, 74 2003; Twomey, 1975), regularization (Hansen & O’Leary, 1993; 75 Talukdar & Swihart, 2003), or linear inversion methods (Hagen & 76 Alofs, 1983), although recently the statistical techniques have been 77 improved (Dubey & Dhaniyala, 2013; Voutilainen, 2001). 78 When the control matrix is strongly ill-conditioned, small errors 79 in the data can be greatly magnified in the solution, even allowing 80 negative values for the PSD. One way to avoid negative values is 81 to use the non-negative least squares (NNLS) algorithm (Lawson & 82 Hanson, 1974), which is a reformulation of the LS solution in which 83 a dual vector w=KT(M−KN) forces the solution to be positive: 84 NN=arg min KNinv −M ;NN≥0.(5) 85 The biggest drawback of this method is the convergence: there 86 is no single factorization and the results can converge to a different 87 local minimum. Moreover, it is only valid for systems where m≥n,88 excluding under-determined problems. 89 The Tikhonov regularization (Willoughby, 1979) is widely used 90 to solve problems of signal processing such as noise reduction, 91 image restoration, and data inversion. The method replaces Eq. (4) 92 by a minimization problem: 93 N␭=N∈Rnarg min KN−M  2+ N  2,(6) 94 where is the regularization parameter, which is a positive con95 stant that balances the fit and smoothness of the distribution and 96 is chosen such that Nbecomes as close as possible to the noise-free 97 solution. 98 The main advantage of the Tikhonov regularization is that diag99 nosis of the problem can be made through the singular value 100 decomposition (SVD). It achieves a pseudo-inverse of the kernel 101 matrix, K+, by disaggregating it into two orthogonal matrices, 102 U∈Rm×mand W∈Rn×n, and a diagonal matrix formed by non103 negative elements, the singular values,˙∈Rm×n:104 K+=W˙+UT= k wkuT k k .(7) 105 Expressing Eq. (6) in vectorized form, the inverted PSD is 106 N␭=(K+K+II)−1K+M= r  k=1 k uT kM k wk,(8) 107 where kare the filter factors or Wiener weights: 108 k=2 k 2 k+.(9) 109 The L-Curve (Hagen & Alofs, 1983; Hansen & O’Leary, 1993; 110 Lloyd, Taylor, Lawson, & Shields, 1997) is a powerful tool for esti111 mating the optimal regularization parameter. In brief, it consists of 112 minimizing Eq. (8) for each value of until the point of maximum 113 curvature is found (Johnston & Gulrajani, 2000). Another regular114 ization method is based on generalized cross-validation (GCV). It 115 Please cite this article in press as: Domat, M., et al. Inversion of electrical mobility measurements using bipolar or unipolar chargers for the arbitrary distribution of channels. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2014.08.007 ARTICLE IN PRESS G Model PARTIC 738 1–10 M. Domat et al. / Particuology xxx (2014) xxx–xxx 3 determines a reduced solution omitting one data point and recal-116 culating. The optimal value is the minimizer of the function:117 NG=||KN␭−M||2 [trace(II −K(K+K+II)−1K+)]2,(10)118 where trace()is the sum of the diagonal elements (Crump & 119 Seinfeld, 1981).120 The sensitivity of the solutions to inaccuracies in the data is121 identified by the condition number (CN), which is the ratio of the122 largest to the smallest singular value other than zero. The larger123 the condition number, the greater the ill-conditioning, which is124 infinite for singular matrices. Its logarithm is an estimate of how125 many digits are lost in the process.126 Geometrical approach to the DMA transfer function127 The probability of particles leaving the DMA with the classified128 flow is determined by the DMA transfer function (TF), ˝q mn. For an129 ideal case, the TF can be fitted by a triangular shape centered at130 the mean mobility Zc.Zcdepends on the voltage Vapplied to the131 DMA and on its geometry, which is assumed to be cylindrical with 132 Riand R0as the inner and outer radii of the DMA electrodes and L133 the length between the entrance slit of the aerosol flow Qaand the134 exit slit of the classified flow Qm:135 Zc=ln (R0Ri) 4L (Qsh +Qex) V.(11) 136 The relative width ˇ0of the TF only depends on the flow rates, 137 considering also the sheath Qsh and excess Qex flowrates, then138 ˇ0=Zc Zc=Qa+Qm Qsh +Qex .(12)139 For the particle detector (i.e., CPC Model 3775) used in the140 present experiments (TSI Inc., 2007), the aerosol-to-sheath air flow141 ratio was kept constant at ˇ0= 0.1, ensuring sufficient resolution142 ᏾=ˇ−1 0.143 The DMA voltage can be varied in stepping or scanning mode:144 DMPS or SMPS (differential or scanning mobility particle sizing).145 SMPS is exponentially ramped at fixed values with a constant ratio146 V(t)=V0e±t/, where is the scan time, tis the residence time, and147 V0the voltage applied to the DMA central rod. At slow scan rates148 (our case, with = 120 s and t= 0.1567 s), the resulting TF for the 149 SMPS system can be approximated to the fixed TF (He & Dhaniyala,150 2013), although it will be shown that the proposed procedure is151 equally valid for both methods.152 The diffusion of particles caused by Brownian motion enhances 153 particle losses, and thereby influences the shape of the TF, because154 it is broadened and resolution is reduced. An analytical expression155 for the relative full width at half maximum (FWHM) of a triangular 156 TF assuming diffusional broadening, , was obtained by Karlsson157 and Martinsson (2003) for any DMA type:158 =1+0.15 ˇ021 2 0−1−12 ,(13)159 where 0is the experimental broadening expression for a Vienna160 type DMA with ˇ0= 0.15 and Pe is the Peclet number:161 0=1.69 1−exp −Pe 105 0.33−0.69.(14)162 Stolzenburg (1988) proposed a realistic bell-shaped TF model, 163 while Stratmann, Kauffeldt, Hummes, and Fissan (1997) showed 164 that the triangular shape is a good approximation because the error165 introduced is below the experimental uncertainty. A geometrical TF 166 Fig. 1. Transfer function for singly charged particles calculated in this work. The mobility intervals mare represented by the central mobility of each triangular TF, and the size intervals nby the vertical lines. Both intervals can be independently chosen, giving a different coverage of the TF area and thus different precision in the inversion process. can be calculated simply knowing the FWHM (namely, Zc) and the 167 height of the fitting triangle H:168 Zc=ˇ0 Zc,(15) 169 H=Pn,(16) 170 where the penetration efficiency Pn(Dp) was described in detail by 171 Reist (1993). The maximum height of the transfer function is ide172 ally 1 when particle losses are ignored. However, as particle size is 173 reduced, diffusion enhances particle losses and the width of the TF 174 broadens as its height decreases, as shown in Fig. 1.175 Estimation of the TF area within every size region (numbered 176 from I to VII in Fig. 2) is calculated comparing the lower and upper 177 boundaries from the TF, Zl(q,m) and Zu(q,m), and from the size 178 intervals, Z∗ l(q, n) and Z∗ u(q, n). Table 1 shows all of the final values, 179 where calculation of section “Results and discussion” is shown for 180 q= 1 as an example. In this particular case, the mobilities follow 181 Zc(m)≥Z∗ u(n)>Z ∗ l(n)>Z l(m), and the corresponding TF area is cal182 culated from the polygonal shape among the total area: 183 A˝(m, n)=1 2a(b+c),(17) 184 Fig. 2. Example of the geometrical area estimation of the TF within the different size intervals. Please cite this article in press as: Domat, M., et al. Inversion of electrical mobility measurements using bipolar or unipolar chargers for the arbitrary distribution of channels. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2014.08.007 ARTICLE IN PRESS G Model PARTIC 738 1–10 4M. Domat et al. / Particuology xxx (2014) xxx–xxx Table 1 Geometrical calculation of the TF area within each size interval. The dependence (m) or (n) with the interval is omitted and the variable J=H/2Zcis defined to simplify the equations. Sections A–C originated when wider size intervals cover wider areas of the TF, described here as a sum of several sections of Fig.2. Section Bound TF area IZ∗ u≤Zl0 II Zc≥Z∗ u>Z l≥Z∗ lJ(Z∗ u−Zl)2 III Zc≥Z∗ u>Z ∗ l≥ZlZ∗ pJ(Z∗ u+Z∗ l−2Zl) IV Zu>Z ∗ u>Z c>Z ∗ l>Z lJ[(Zu−Z∗ u)(Z∗ u−Zc)+ (Z∗ l−Zl)(Zc−Z∗ l)+ Z∗ pZc] VZu>Z ∗ u>Z ∗ l≥ZcZ∗ pJ(2Zu−Z∗ u−Z∗ l) VI Z∗ u≥Zu>Z ∗ l≥ZcJ(Zu−Z∗ l)2 VII Z∗ l≥Zu0 A: (II + III + IV) Zu>Z ∗ u>Z c>Z l≥Z∗ lJ(2Z2 c−(Zu−Z∗ u)2) B:(IV+V+VI) Z∗ u≥Zu>Z c>Z ∗ l>Z lJ(2Z2 c−(Z∗ l−Zl)2) C:(II+III+IV+V+VI) Z∗ u≥Zu>Z c>Z l≥Z∗ lHZc where a=Z∗ p(n)=Z∗ u(n)−Z∗ l(n) and the unknown heights band185 care 186 b=H(m) Zc(m)[Z∗ u(n)−Zl(m)] and c=H(m) Zc(m)[Z∗ l(n)−Zl(m)].(18)187 Substituting Eq. (18) into Eq. (17), the TF area is obtained188 depending only on known mobilities:189 A˝(m, n)=H(m)Z∗ p(n) 2Zc(m)(Z∗ u(n)+Z∗ l(n)−2Zl).(19)190 In many instruments, such as in TSI SMPS-3080 series devices,191 the number of channels per decade for both input and output 192 data are fixed to a multiple of 4, indicating that m=n.Inthe193 present work, the number of mobility measurements and output 194 size channels are independent and can be chosen arbitrarily, giving 195 overlapping TF intervals that can cover more mobilities when the196 TF is broadened by diffusion.197 The number of intervals is estimated through the number of 198 channels per decade measured by the instrument; they are related199 by the expression of the ratio between diameters for the PSD log-200 distribution, d= ln(Dp(i)/Dp(i−1)):201 d=ln (DpUp/DpLow) num Decades =ln(DpMax/DpMin) num Intervals ,(20) 202 where DpUp and DpLow are usually in a ratio of DpUp/DpLow = 10 for203 the decade, and DpMax and DpMin are the boundary diameters of the204 PSD range. Modifying the number of intervals (mor n), the covered205 area of the TF will change accordingly. 206 The TF is crucial to control matrix estimation, which also 207 depends on the fraction of charge fc(q,Dp) and the input to output208 flow ratio. The variation of these values, along with N(Dp), is small209 across the width of the integration intervals. From the geometrical 210 TF estimation, ˝(q,Zp,Zc)=q˝(Zp,Zc), Eq. (2) becomes211 dM(DZp) dln DZp=Qa Qm q qfc(q, Dp)Pn(Dp)dN(Dp) dln Dp 212 ×Dp max Dp min ˝(Zp,Z c)dDp+ε(DZp).(21) 213 214 Changing the integration variable simplifies the calculation of215 the TF area: 216 Zp max Zp min ˝(Zp,Z c) dln Dp dZp dZp=A˝(m, n).(22)217 FromFig. 3, the proposed geometric method shows similar218 results to the bell-shaped TF for the total area, with a relative error219 Fig. 3. Comparison of the TF areas calculated with the bell-shaped model from Stolzenburg (1988) and the triangular shape for the ideal case (upper set of curves) and considering losses (lower set of curves), including the SMPS analysis estimated with the geometric TF method. Both methods cover practically the same area with a difference of about 1%. below 1% for particles >2 nm for both ideal (no losses) and real 220 approximations. 221 For SMPS measurements, the mobility shift depends on the ratio 222 between residence and scanning times, which was experimentally 223 determined by Dubey and Dhaniyala (2008) for a nano-DMA, where 224 the area of the triangular scanning TF is time-dependent: 225 A˝(t)=1 2tH =ˇ01+k1 t +k2t 2.(23) 226 The constants k1and k2depend on the flow ratio, and the height 227 H=2A˝(t)/tis the same for the temporal and mobility based TFs. 228 Combined with the expression for the voltage scan ratio, the tem229 poral and mobility-based areas of the TF are related by 230 A˝(Zp)=A˝(t)Zp ln(Zu/Zl).(24) 231 The results for the ideal and real cases are shown in Fig. 3. The 232 differences between the results from scanned and fixed voltage 233 areas, which are less than 2% for particle sizes above 15 nm when 234 losses are considered, are probably because of the dependence of 235 the temporal area A˝(t) on the scan time and the mobility shift in 236 the conversion to A˝(Zp). 237 Particle size and charge distribution retrieval from any 238 charging method 239 Calculation of the PSD requires accurate information about 240 the charge distribution of the incoming aerosol. Although bipo241 lar radioactive chargers have been traditionally used because of 242 their well-defined charge distributions, unipolar diffusion chargers 243 attain higher charging efficiencies while avoiding the severe safety 244 restrictions of bipolar radioactive chargers. Some recent unipo245 lar devices are based on X-rays (Kim, Kim, Kwon, & Park, 2011), 246 UV photoionization (Honta˜ non & Kruis, 2008), corona discharge 247 (Chien, Tsai, Chen, Lin, & Wu, 2011; Domat, Kruis, & Fernandez248 Diaz, 2014b; Intra, 2012; Li & Chen, 2011; Qi, Chen, & Pui, 2007; 249 Tsai, Lin, Chen, Huang, & Alonso, 2010; Vivas, Honta˜ non, & Schmidt250 Ott, 2008)and others (Kwon, Fujimoto, Kuga, Sakurai, & Seto, 2005; 251 Kwon, Sakurai, & Seto, 2007). 252 However, inversion methods are seldom applied to unipo253 lar chargers, possibly because these devices are often in the 254 experimental phase and there is no homogeneous expression for 255 their charge distribution, such as there is for bipolar radioactive 256 Please cite this article in press as: Domat, M., et al. Inversion of electrical mobility measurements using bipolar or unipolar chargers for the arbitrary distribution of channels. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2014.08.007 ARTICLE IN PRESS G Model PARTIC 738 1–10 M. Domat et al. / Particuology xxx (2014) xxx–xxx 5 Table 2 Coefficients of the log-normal fitting of the charge fractions from the corona charger reported by Domat et al. (2014b). Dgap (mm) qf 0(Dp)¯ Dpg (nm) gCR 2 3 1−0.02065 20.4725 1.6373 13.44060 0.98987 2 0.00192 33.6301 1.3967 12.86995 0.96662 3 0.00131 47.0735 1.2406 10.62972 0.99802 4 0.00125 69.9434 1.2823 16.75291 0.99868 5 6.69 ×10−5102.5040 1.3247 25.18727 0.99867 10 1−0.02375 66.9368 2.0618 48.98775 0.96978 2−2.44 ×10−484.9966 1.5394 26.87955 0.99990 3−7.78 ×10−5134.0929 1.5506 23.01109 0.99934 chargers (Wiedensohler et al., 2012). Furthermore, the ratio of mul-257 tiply charged particles in unipolar chargers is higher than in bipolar258 chargers, complicating the multiple charge correction in the decon259 volution.260 Once the charge distribution is known, it can be fitted by a261 reproducible expression (Kaminski et al., 2012). There are analytic262 models to estimate the charge fraction, such as the birth-and-263 death model (Adachi, Kousaka, & Okuyama, 1985; Biskos, Reavell, & 264 Collings, 2005; Boisdron & Brock, 1970; Fuchs, 1964). However, we 265 chose simpler empirical fitting (the log-normal distribution shown266 in Eq. (25)) because analytic models do not account for the loss of 267 particles and ions in a first approximation and these assumptions268 introduce more uncertainty in the inversion process when fitting269 to the experimental data. 270 fc(Dp)=f0(Dp)+C √2ln gDp exp −ln (Dp/¯ Dpg)2 2ln2 g,(25)271 where a possible offset is denoted by f0(Dp), Cis the probability of272 acquiring charges, ¯ Dpg is the geometric mean diameter, and gis273 the geometric standard deviation.274 Charge fractions from a previously developed unipolar corona275 charger (Domat et al., 2014b) were fitted. The charger had sep-276 arate ionization and particle charging zones, and its charging 277 performance was obtained experimentally with concentrations of278 monodisperse particles between 6 and 60 nm below 106cm−3. The279 discharge electrode was a needle mounted on a micrometer screw, 280 allowing the distance to the walls (namely, Dgap) to change, thereby281 varying the Nit-product in the charging zone between 2.4 ×1012 282 and 1.8 ×1013 s/m3(Domat, Kruis, & Fernandez-Diaz, 2014a). Thus,283 different charging levels can be attained, with the aim of improving 284 the concentration of singly charged particles. The fitting parame-285 ters are listed in Table 2 for 3 and 10 mm gap distances, along with286 the regression coefficient R2.287 Multiple charge correction288 One of the difficulties of accurate retrieval of the PSD is the289 multiple charging of particles, which causes a proportion of the290 particles measured under the same applied voltage to correspond291 to particles with different diameters. Hoppel (1978) developed a 292 technique in which Eq. (3) can be disaggregated depending on the293 charge contribution:294 Mq m≈ ntot  n=1(K1 mnN1 n+ε1 m)+ qtot  q>1 (Kq mnNq n+εq m)≈M1 m+Mq>1 m.295 (26)296 297 Therefore, retrieval of the singly charged PSD would consist of298 inverting M1 m=Mq m−Mq>1 mfollowing some initial steps. First, esti-299 mation of M1 massuming only singly charged particles, leading to a300 first approximation of N1 nthat will be used to calculate the contri301 bution of the penultimate channel to the mobility measurements, 302 Mq>1 m−1.. The singly charged distribution is re-calculated considering 303 N1 n≈[K1 mn]+(Mq m−Mq>1 m−1), and this process is repeated for every 304 mobility channel mto successively correct the presence of different 305 particle sizes. 306 The algorithm of Hoppel supposes that an impactor is placed to 307 remove the large particles that can contribute multiple charges and 308 distort the measurements. This algorithm was recently updated by 309 He and Dhaniyala (2013), who added an intermediate step after the 310 first approximation that increases the accuracy when no impactor 311 is applied. They fitted the large diameters of the zeroth-order singly 312 charged PSD by a Gumbel distribution function, and the resultant 313 fit was used to correct the multiply charged contribution. However, 314 this approach is only suitable for smooth, unimodal distributions, 315 otherwise the Hoppel algorithm is applied. 316 Retrieval of the charge distribution in the TDMA setup 317 Tandem differential mobility analysis (TDMA) (Liu et al., 1978; 318 Rader & McMurry, 1986) is a popular technique to characterize 319 changes in the shape, size, concentration, and charge of particles 320 or to calibrate devices. The arrangement consists of two identi321 cal DMAs in series with a particle conditioner in between. In our 322 method, we used the aforementioned corona charger, where a fixed 323 voltage applied on the first DMA classifies particles of a certain size, 324 and the resulting distribution is neutralized and then recharged 325 with the new charger to be measured with the second DMA. 326 The response after the second DMA is analogous to the response 327 after the first DMA (Eq. (21)). However, this has to be expressed 328 taking into account the TF and other characteristics of the previous 329 DMA system: 330 dM2(DZp) dln DZp=Qa1 Qm1 Qa2 Qm2  q qfc2(q, Dp)Pn2(Dp)dN2(Dp) dln Dp 331 ×Dp max Dp min ˝1(Zp,Z c1)˝2(Zp,Z c2)dDp+ε(DZp).(27) 332 333 Usually, the charging properties of the test charger are 334 determined through a TDMA arrangement by comparing the con335 centrations at different points of the setup (see Chien et al., 2011 336 or Domat et al., 2014b for details). However, these properties are 337 subjected to the uncertainties of the concentration measurements. 338 Considering that the mean charge can be estimated as 339 ¯ q2(Dp)= q qfc2(q, Dp).(28) 340 A new method to calculate the mean charge is proposed, that is, 341 through inversion of the response from the second DMA (Eq. (27))342 once the PSD input in the first DMA is known. 343 For proper interpretation of the inverted results, the known sys344 tematic error (up to 2%) between the mean mobility diameters 345 measured by DMA-1 and DMA-2 has to be considered (Rader & 346 McMurry, 1986). This error can be partially corrected by adjusting 347 the voltages of each DMA (Alonso & Kousaka, 1996), although about 348 1% of the shift is unavoidable (Domat et al., 2014b). In addition, vari349 ations of the concentration between DMAs because of diffusion and 350 electrostatic effects are produced by the charging-neutralization 351 process: 352 dN2(Dp) dln Dp=Pn0(Dp)dN1(Dp) dln Dp .(29) 353 Please cite this article in press as: Domat, M., et al. Inversion of electrical mobility measurements using bipolar or unipolar chargers for the arbitrary distribution of channels. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2014.08.007 ARTICLE IN PRESS G Model PARTIC 738 1–10 6M. Domat et al. / Particuology xxx (2014) xxx–xxx Fig. 4. Simulated and inverted PSD for ranges measurable with nano-DMA (a) and long-DMA (b), and fraction of singly charged particles among the total number of charged particles (c) charged for the bipolar radioactive distribution and recovered from the previous distributions with the multiple charge correction algorithm. The losses according to fittings of the experimental data for 3354 and 10 mm gap distances are represented by Pn0, respectively, as355 Pn0 =0.83 −0.89e−0.09Dp(R2=0.984),(30)356 Pn0=0.93 −0.85e−0.13Dp(R2=0.962).(31) 357 Because the tubing to the particle counter is the same for DMA-1358 and DMA-2, Pn1 =Pn2 (namely Pn), and Eq. (27) becomes 359 dM2(DZp) dln DZp=Qa1 Qm1 Qa2 Qm2 ¯ q2(Dp)Pn(Dp)Pn0(Dp)dN1(Dp) dln Dp 360 ×Dp max Dp min ˝1(Zp,Z c1)˝2(Zp,Z c2)dDp+ε(DZp).(32)361 362 Then, analogous to PSD inversion, the mean charge can be 363 deconvoluted:364 Mm2≈˜ Kmn ¯ qn+εm,(33) 365 where ˜ Kis the new definition of the kernel matrix according to Eq.366 (32).367 A similar procedure has been used by Li, Li, and Chen (2006) and368 Stolzenburg and McMurry (2008) to deconvolute the TF of particle369 sizing devices, and could also be applied here. However, because 370 no reports of charge distribution deconvolution were found, and 371 it is an important issue for PSD characterization, this method was372 chosen among the others.373 Results and discussion374 The proposed inversion theory was tested with two types of 375 data: simulated and experimental. Simulated data were generated376 to test the validity of the techniques by comparing the results with 377 reference distributions. Both the radioactive and corona charging 378 techniques were tested and inverted by the NNLS and SVD methods 379 comparing with the regularization parameters from Section “Data 380 inversion technique”. 381 The accuracy of the simulated inversion was estimated by the 382 normalized root-mean-squared error (NRMSE), which gives a single 383 value for the mean of the standard deviation of the residuals from 384 the inversion between the original (simulated, N) and calculated 385 data (N*). That is, 386 NRMSE =1 Nmax −Nmin     1 n n  i=l (Ni−N∗ i)2.(34) 387 On the other hand, in experimental measurements, both the 388 uncertainty and the original data are unknown. In these cases, the 389 precision of the inverted results was calculated by averaging the 390 results of several tests from the same set of data to give an estima391 tion of the deviation from the original data. 392 Results from the simulated data 393 Two representative sets of log-normal PSD (Eq. (25) without 394 offset) charged with a bipolar 85Kr charger were tested. One for a 395 geometric mean diameter of 50 nm with 1.3 standard deviation and 396 a total concentration of C= 3.5 ×1011 m−3, and the other for a larger 397 diameter range ( ¯ Dpg =300 nm, g=1.5,C=4.5×1015 m−3),398 each within the measurable ranges of the nanoand long-DMA, 399 respectively. Some uncertainty was added to imitate the mea400 surement process: around 5% for ambient aerosol and up to 10% 401 for generated aerosol (Gysel, McFiggans, & Coe, 2009). However, 402 Please cite this article in press as: Domat, M., et al. Inversion of electrical mobility measurements using bipolar or unipolar chargers for the arbitrary distribution of channels. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2014.08.007 ARTICLE IN PRESS G Model PARTIC 738 1–10 M. Domat et al. / Particuology xxx (2014) xxx–xxx 7 Fig. 5. (a) Inverted PSD for the unipolar corona charger at different distances between electrodes and (b) inverted mean charge for the unipolar and bipolar test chargers in a TDMA configuration. phenomena such as coagulation and agglomeration were not con-403 sidered.404 In Fig. 4(a) and (b), the inversion of the determined problem405 with 64 channels/decade and 5% noise is shown. In the smaller 406 diameter range (50 nm, Fig. 4(a)), both SVD and NNLS methods407 lead to the same indistinguishable result with identical deviation408 and an NRMSE below 0.5%. For the larger diameter range (300 nm,409 Fig. 4(b)), the inversion performed by SVD was slightly better than 410 that by NNLS (NRMSESVD = 0.8%, NRMSENNLS = 1.2%). The L-Curve411 (LC) and GCV regularization methods also gave very similar results,412 although the LC had a lower NRMSE than GCV even for the nois413 iest data. Because the performance with the different methods is 414 reasonable, only the SVD-LC will be shown in the following.415 The results of the multiple charge correction algorithm from416 Section “Multiple charge correction” are shown in Fig. 4(c). The 417 singly and multiply charged distributions are very similar in the418 small particle range, but become different as the particle size419 increases and acquires more charge. The theoretical fraction of 420 singly charged particles among the total number of charged parti421 cles (Wiedensohler, 1988; Wiedensohler et al., 2012) is compared422 with the results from previous plots in Fig. 4(c). There is good agree-423 ment between the predicted and obtained results, although for the 424 largest diameters our prediction is clearly underestimated, prob-425 ably because of the decrease in the concentration in the tails of426 the PSD and the possible contribution of large particles with more 427 charges. 428 The multimodal distributions are shown in Fig. 5(a) ( ¯ Dpg =429 6,10,and 23 nm; g= 1.18, 1.3, and 1.05; C=(3.0, 5.0, and430 1.2) ×1010 m−3, respectively). The data were simulated with431 unipolar charging at two different Nitvalues given by the432 Fig. 6. Precision of the results for (a) an increase in artificial uncertainty during the measurements and (b) a variable ratio of mobility measurement (input, m) intervals to size (output, n) intervals. The NRMSE and the logarithm of the condition number CN give an indication of the accuracy of the inversion process. gap distance. The inversion performed better with D=10mm 433 (NRMSE = 1.6%) than with D= 3 mm gap (NRMSE = 3%), because a 434 higher Nit-product, and thus a higher level of particle charging, 435 leads to more singularity in the estimation of the kernel matrix. 436 The mean charge of a test charger in a TDMA arrangement was 437 estimated from the recovered PSD charged in a first DMA sys438 tem with a known bipolar 85Kr charge distribution. In the second 439 DMA system, three types of charge distribution were tested: bipolar 440 85Kr and unipolar corona charging with gaps at 3 and 10 mm. The 441 recovered results are shown in Fig. 5(b). The bipolar mean charge 442 deconvoluted with a NRMSE of 2%, while the unipolar charging was 443 slightly underestimated with a NRMSE of 8% for Dgap =10mmand 444 up to 20% NRMSE for Dgap = 3 mm. It is important to notice from Eq. 445 (32) that there is a recovered function only for the inverted N1(Dp)446 values that are not equal to zero. Moreover, the uncertainties in the 447 inversion of the PSD will propagate to the inversion of the mean 448 charge. 449 In Fig. 6(a), the variation of the error in the inversion process 450 with the artificial noise is shown. It is expected that the greater 451 the added noise, the greater the error in the inversion, as occurs in 452 PSD inversion. However, although the inverted mean charge gener453 ally has higher error, it is almost insensitive to added noise, which 454 indicates that the main contributors to the residuals in the second 455 inversion are the errors from the first inversion. 456 While the CN seems to be generally insensitive to errors and 457 noise during measurement of the data, as shown in the lower 458 graph of Fig. 6(a), the number of intervals directly affects the size 459 Please cite this article in press as: Domat, M., et al. Inversion of electrical mobility measurements using bipolar or unipolar chargers for the arbitrary distribution of channels. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2014.08.007 ARTICLE IN PRESS G Model PARTIC 738 1–10 8M. Domat et al. / Particuology xxx (2014) xxx–xxx Fig. 7. Deconvolution in left axis of the mobility distributed data from right axis for the over-determined, under-determined, and determined problems. The data were acquired in a DMPS system for symmetric flows of QQ= 1.5 lpm and Qsh =15lpm with 64 channels per decade. of the matrix and thus the accuracy of the results, following the460 expected trend of being more precise for over-determined than 461 under-determined systems (Fig. 6(b)). The error in the inversion462 of the PSD decreases to <2% when m/n≥1. The same occurs for463 the inverted mean charge, but, as in previous cases, there is an464 accumulated uncertainty from the previous inversion that the sec-465 ond calculation cannot overcome, reaching a lower “limit” for the 466 NRMSE.467 It must be mentioned that, although extra equations in the468 over-determined problem can be used as additional information 469 to reduce the high levels of uncertainty in the data, they can over470 smooth the retrieved data, losing information about narrow peaks.471 Moreover, it is common in real measurements to have a limited472 number of channels to measure, and it is therefore important to be473 able to solve under-determined problems. However, in the lower474 plot of Fig. 6(b), the logarithm of the CN shows lower precision loss475 for the under-determined problems than for the determined prob-476 lems, and first increases and then decreases with increasing m/n.477 This can lead to confusion, but it is caused because the most under-478 determined case has such a lack of information that the difference479 between the smallest and largest singular values is minimal.480 Results from experimental measurements481 The results from simulated measurements are useful to deter-482 mine the applicability of the proposed method to experimental data 483 and to analyze the factors that affect the inversion process. The next 484 step is to test experimental data from bipolar or unipolar charging485 measured by stepping and slow-scanning voltages in the DMA.486 The inversion of a set of data measured in a DMPS system work487 ing with a bipolar radioactive 85Kr charger is shown in Fig. 7. Three488 ratios of intervals (m/n) were tested. The under-determined and489 determined problems gave results with similar shapes and areas to 490 the measured mobility distribution. However, the over-determined 491 problem showed a different shaped curve near the tails, with two492 small peaks probably because of multiply charged particles that493 entered with the selected mobility in the classifier. This indicates 494 the importance of being able to modify the ratio of intervals.495 Measurements in a TDMA configuration are always performed496 for monodisperse distributions, sometimes with very narrow peaks497 and low concentration of particles, which complicates the inver-498 sion. In Fig. 8, a set of monodisperse measurements used to499 calibrate the corona charger (with Dgap =10 mm) is shown in a500 Fig. 8. Measured response and inverted PSD of the set of monodisperse measurements used to calibrate the corona charger in the TDMA arrangement. TDMA configuration at slow-scan rates to avoid smearing effects. 501 The monodisperse distributions were treated together to form a 502 wider distribution and to use the information of each of the single 503 peaks in the inversion process of the response. The recovered result 504 from the PSD inversion is also shown in Fig. 8, and, although the 505 Fig. 9. (a) Fraction of singly charged particles among the total charged particles recovered with the multiple charge correction algorithm and (b) inverted mean charge for the unipolar corona charger with Dgap = 10 mm in a TDMA configuration. Please cite this article in press as: Domat, M., et al. Inversion of electrical mobility measurements using bipolar or unipolar chargers for the arbitrary distribution of channels. Particuology (2014), http://dx.doi.org/10.1016/j.partic.2014.08.007 ARTICLE IN PRESS G Model PARTIC 738 1–10 M. Domat et al. / Particuology xxx (2014) xxx–xxx 9 singly charged particle distribution was calculated, its representa-506 tion results are almost indistinguishable from the total.507 Separately representing the fraction of particles with charge508 unity (Fig. 9(a)) indicates that this fraction is overestimated509 compared with the expected value to a similar extent to the under-510 estimation of the doubly charged particles. A reason for this could511 be that it was only possible to measure up to three charges in the512 charge distribution characterization, and therefore higher particle513 charges are present in the measurements that are not being taken514 into account.515 The retrieval of the mean charge from the PSD recovered in516 Fig. 8 through the method presented in Section “Retrieval of the517 charge distribution in the TDMA setup” is shown in Fig. 9(b). Oscil-518 lations due to errors and uncertainties in the data appear below 519 20 nm, which overestimates the results. However, for large par-520 ticle sizes, the differences with the expected values are reduced,521 probably because of a higher concentration and level of charge of522 the particles. 523 The diameter range below 20 nm is always the most sensitive524 range because of the forces to which particles are subjected, such525 as diffusion and electrostatic effects. These phenomena are compa-526 rable to the fluctuations of the distribution, and the algorithm has527 to be able to distinguish between information from the real dis-528 tribution and noise and uncertainties. We have to also taken into529 account that these particle ranges usually have the lowest charging530 efficiencies, and thus have low concentrations. Therefore, inversion531 of size ranges where diffusion plays an important role always have 532 large uncertainty, which is even larger if the error is able to prop533 agate from the first to the second inversion procedures, and thus534 precise estimation requires more effort.535 Conclusions536 A comprehensive study of the requirements for obtaining accu-537 rate inversion of data from mobility measurements was conducted 538 based mainly on NNLS and SVD algorithms, with the SVD algorithm539 comparing two methods of regularization: L-Curve and GCV. The540 methodology was developed to be used with any type of charger,541 once its charge distribution is defined, and an approximation is 542 performed to use data from scanning or stepping voltage config-543 urations in the DMA. Great flexibility is allowed in choosing the544 number of intervals and the number of data channels, allowing545 arbitrary voltage steps and overlapping DMA transfer functions. 546 The transfer function is estimated by a simple geometrical cal-547 culation from the areas within each output interval, requiring only548 the mobilities that define the TF, which is an accurate approxima-549 tion compared with the bell-shaped approximation. The analysis550 of the kernel matrix is a key to understanding the origin of the551 instabilities, because its rank conditions the number of potential552 mathematical solutions to the problem. 553 From direct comparison of the SVD and NNLS methods, it was554 found that the errors were similar and usually depended on the555 characteristics of the distribution. However, the SVD method com-556 bined with the L-Curve regularization has the advantage of greater 557 flexibility in selecting the number of channels, resulting in better558 performance because it is capable of accurately solving even under-559 determined systems. Inversion performed in this way gave a good 560 fit with less than 1% error for noise levels of 5%, which could be561 even less when the intervals are properly chosen.562 A method to calculate the contribution of singly charged563 particles to the total distribution was successfully applied to dis-564 tributions from unipolar and bipolar charging. However, in the 565 tails of the distributions, where the concentration is low, the frac-566 tion of singly charged particles was misestimated. The technique567 was also applied to the TDMA arrangement with the possibility 568 of retrieving the mean charge per particle of the distribution. The 569 results gave a good approximation when the original PSD was accu570 rately estimated, although errors can propagate and distort the 571 inverted mean charge and further method development is required 572 to reduce the error propagation. 573 Acknowledgments 574 M. Domat and J.M. Fernandez-Diaz thank the Ministerio de EduQ4575 cacion y Ciencia (MEC) of Spain for support under the project 576 MEC05CGL2005-05244/CLI and grant BES-2006-12469. 577 F.E. Kruis is financially supported by the Deutsche Forschungs578 gemeinschaft (DFG) in the framework of the joint research 579 program “Multiparameter Characterization of Particle-based Func580 tional Materials by Innovative Online Measurement Technology” 581 (PAK688). 582 References 583 Adachi, M., Kousaka, Y., & Okuyama, K. (1985). Unipolar and bipolar diffusion charg584 ing of ultrafine aerosol particles. Journal of Aerosol Science,16(2), 109–123. 585 Alofs, D. J., & Balakumar, P. (1982). Inversion to obtain aerosol size distributions from 586 measurements with a differential mobility analyzer. Journal of Aerosol Science,587 13(6), 513–527. 588 Alonso, M., & Kousaka, Y. (1996). Mobility shift in the differential mobility analyzer 589 due to Brownian diffusion and space-charge effects. Journal of Aerosol Science,590 27(8), 1201–1225. 591 Bazan, F. S. V., & Francisco, J. B. (2009). An improved fixed-point algorithm for deter592 mining a Tikhonov regularization parameter. Inverse Problems,25(4), 045007. 593 Biskos, G., Reavell, K., & Collings, N. (2005). Unipolar diffusion charging of aerosol 594 particles in the transition regime. Journal of Aerosol Science,36(2), 247–265. 595 Boisdron, K., & Brock, J. R. (1970). On the stochastic nature of the acquisition of elec596 trical charge and radioactivity by aerosol particles. Atmospheric Environment,597 4(1), 35–50. 598 Chien, C. L., Tsai, C. J., Chen, H. L., Lin, G. Y., & Wu, J. S. (2011). Modeling and validation 599 of nanoparticle charging efficiency of a single-wire corona unipolar charger. 600 Aerosol Science and Technology,45, 1468–1479. 601 Collins, D. R., Flagan, R. C., & Seinfeld, J. H. (2002). Improved inversion of scanning 602 DMA data. Aerosol Science and Technology,36(1), 1–9. 603 Crump, J., & Seinfeld, J. (1981). A new algorithm for inversion of aerosol size distri604 bution data. Aerosol Science and Technology,1(1), 15–34. 605 Domat, M., Kruis, F., & Fernandez-Diaz, J. (2014a). Determination of the relevant 606 charging parameters for the modelling of unipolar chargers. Journal of Aerosol 607 Science,71, 16–28. 608 Domat, M., Kruis, F., & Fernandez-Diaz, J. (2014b). Investigations of the effect of 609 electrode gap on the performance of a corona charger having separated corona 610 and charging zones. Journal of Aerosol Science,68, 1–13. 611 Dubey, P., & Dhaniyala, S. (2008). Analysis of scanning DMA transfer functions. 612 Aerosol Science and Technology,42(7), 544–555. 613 Dubey, P., & Dhaniyala, S. (2013). Improved inversion of scanning electrical mobility 614 spectrometer data using a new multiscale expectation maximization algorithm. 615 Aerosol Science and Technology,47(1), 69–80. 616 Fiebig, M., Stein, C., Schröder, F., Feldpausch, P., & Petzold, A. (2005). Inversion of data 617 containing information on the aerosol particle size distribution using multiple 618 instruments. Journal of Aerosol Science,36(11), 1353–1372. 619 Fuchs, N. (1964). The mechanics of aerosols. New York: Pergamon Press. 620 Gysel, M., McFiggans, G. B., & Coe, H. (2009). Inversion of tandem differential mobility 621 analyser (TDMA) measurements. Journal of Aerosol Science,40(2), 134–151. 622 Hagen, D. E., & Alofs, D. J. (1983). Linear inversion method to obtain aerosol size 623 distributions from measurements with a differential mobility analyzer. Aerosol 624 Science and Technology,2(4), 465–475. 625 Hansen, P. C., & O’Leary, D. P. (1993). The use of the L-curve in the regularization of 626 discrete ill-posed problems. SIAM Journal on Scientific Computing,14, 1487–1503. 627 He, M., & Dhaniyala, S. (2013). A multiple charging correction algorithm for scanning 628 electrical mobility spectrometer data. Journal of Aerosol Science,61, 13–26. 629 Honta˜ non, E., & Kruis, F. E. (2008). Single charging of nanoparticles by UV photoion630 ization at high flow rates. Aerosol Science and Technology,42, 310–323. 631 Hoppel, W. A. (1978). Determination of the aerosol size distribution from the mobil632 ity distribution of the charged fraction of aerosols. Journal of Aerosol Science,9(1), 633 41–54. 634 Intra, P. (2012). Corona discharge in a cylindrical triode charger for unipolar diffusion 635 aerosol charging. Journal of Electrostatics,70(1), 136–143. 636 Johnston, P. R., & Gulrajani, R. M. (2000). Selecting the corner in the L-curve approach 637 to Tikhonov regularization. IEEE Transactions on Biomedical Engineering,47,638 1293–1296. 639 Kaminski, H., Kuhlbusch, T. A. J., Fissan, H., Ravi, L., Horn, H.-G., Han, H.-S., 640 et al. (2012). Mathematical description of experimentally determined charge 641