Assessment of lidar depolarization uncertainty by means of a polarimetric lidar simulator

. Lidar depolarization measurements distinguish between spherical and non-spherical aerosol particles based on the change of the polarization state between the emitted and received signal. The particle shape information in combination with other aerosol optical properties allows the characterization of different aerosol types and the retrieval of aerosol particle microphysical properties. Regarding the microphysical inversions, the lidar depolarization technique is becoming a key method since particle shape information can be used by algorithms based on spheres and spheroids, optimizing the retrieval procedure. Thus, the identiﬁcation of the depolarization error sources and the quantiﬁcation of their effects are crucial. This work presents a new tool to as-sess the systematic error of the volume linear depolarization ratio ( δ) , combining the Stokes–Müller formalism and the complete sampling of the error space using the lidar model presented in Freudenthaler (2016a). This tool is applied to a synthetic lidar system and to several EARLINET lidars with depolarization capabilities at 355 or 532 nm. The lidar systems show relative errors of δ larger than 100 % for δ values around molecular linear depolarization ratios ( ∼ 0 . 004 and up to ∼ 10 % for δ = 0 . 45 ) . However, one system shows only relative errors of 25 and 0.22 % for δ = 0 . 004 and δ = 0 . 45, respectively, and gives an example of how a proper identiﬁcation and reduction of the main error sources can drastically reduce the systematic errors of δ . In this regard, we provide some indications of how to reduce the systematic errors.


Introduction
The lidar depolarization technique is a useful tool for different applications in atmospheric science, such as the identification of the thermodynamic phase of clouds (e.g. Ansmann et al., 2005;Reichardt et al., 2003;Schotland et al., 1971) and the aerosol typing (e.g. Bravo-Aranda et al., 2013;Groß et al., 2011Groß et al., , 2012Groß et al., , 2013Groß et al., , 2015Navas-Guzmána, b et al., 2013). Additionally, the lidar depolarization technique is very important for improving the retrieval of microphysical aerosol properties (e.g. Ansmann et al., 2011;Chaikovsky et al., 2002;Granados-Muñoz et al., 2014;Wagner et al., 2013;Samaras et al., 2015), becoming crucial for inversion algorithms based on modelling aerosol particles as spheres and spheroids. Unfortunately, the reliability of the lidar depolarization technique is limited due to the complexity of the depolarization calibration.
A relative polarization calibration using molecular depolarization as a reference introduces a high uncertainty due to possibly residual depolarizing aerosol in the assumed clean air range and due to the low signal-to-noise ratio. Absolute calibration methods are more robust, but often do not take into account all polarization effects of the lidar optics, as e.g. the polarization-dependent receiver transmission (Mattis et al., 2009). Many authors have focussed their effort on the improvement of the lidar polarization calibration and on the determination of the depolarization uncertainties (e.g. Álvarez et al., 2006;David et al., 2012;Bravo-Aranda et al., 2013;Hayman and Thayer, 2012;Freudenthaler et al., 2009;Freudenthaler, 2016a). The identification of the lidar optics' influence on the depolarization measurements is relevant (i) for an appropriate assessment of depolarization ratio errors, (ii) for the prioritization of the technical improvements of the lidars according to the optical parameters with the larger impact on depolarization ratio errors, and (iii) for the development of future lidar generations. To this end, the Stokes-Müller formalism is used to provide the theoretical framework and the formulae for the polarization calibration factor covering the different calibration techniques and lidar set-ups used in this work (Freudenthaler, 2016a).
This work quantifies the volume linear depolarization ratio uncertainty ( δ) due to the unknown systematic errors caused by the uncertainties of the polarization-relevant parameters of the lidar optics and assesses the contribution of each lidar functional block to the total uncertainty. A software tool called Polarimetric Lidar Simulator (PLS) has been developed, based on the theoretical framework presented by Freudenthaler (2016a), which performs a complete search of the error space. The PLS is applied to several lidar systems in order to show the dependence of the systematic error on their design features. Random errors due to signal noise are neglected in this work. Their contribution to the uncertainty can be derived, in a similar way to Pappalardo et al. (2004) and Guerrero-Rascado et al. (2008), by simulating the com- plete error space of δ due to the uncertainties of the lidar parameters.
In Sect. 2 we describe the functional blocks. In Sect. 3, the PLS performance is explained in detail. In Sect. 4 a synthetic lidar set-up is used to identify the most important error sources. In Sect. 5, the total systematic error of δ is determined for several EARLINET lidar systems, the error sources are analysed, and possible ways to reduce the uncertainties are pointed out. Finally, conclusions are reported in Sect. 6.

Basis of the Polarimetric Lidar Simulator (PLS):
Stokes-Müller formalism applied to lidar functional blocks As described by Freudenthaler (2016a), we subdivide our lidar system model into functional blocks with the corresponding Stokes vector and Müller matrices: the laser I L , the emitting optics M E (beam expander, steering mirrors), the receiving optics M O (telescope, collimator, dichroic mirrors . . . ), and the polarizing beam splitter M S including the detectors. M S is further split into the transmitted path M T and the reflected path M R with the corresponding Stokes vectors I T and I R . Furthermore, the polarization calibrator C can be considered as an additional functional block. Figure 1 shows a scheme of the lidar based on the functional blocks. We also include a rotator R to enable different lidar set-ups as described below. The next sections explain each functional block, describe the assumptions, and show their Stokes vectors or Müller matrices.

Laser, I L
The specified "polarization purity" of lasers commonly used in lidars is of the order of 100 : 1, if it is specified at all. Already the terminology indicates that such specifications are rather vague. Actually, laser manufacturers do not measure the state of polarization of the laser beams, and seem to give values which are on the safe side under all circumstances. Theoretically, non-linear crystals as second and third harmonic generators should provide very clean linear polariza-tion, just depending on the quality and accuracy of alignment of the crystals. Due to a lack of detailed information we neglect this source of errors in this work. However, in order to remove this uncertainty, in some lidar systems high-quality polarizing beam splitters are used to improve the degree of linear polarization of the emitted laser beams. In both cases, the plane of polarization of the laser beam can be rotated by angle α with respect to the incident plane of the polarizing beam splitter in the receiver optics, which results in the Stokes vector I L of the emitted laser beam: where I L is the laser energy.

Laser emitting optics, M E
The emitting functional block can consist of a set of dichroic beam splitters and steering mirrors, which align the laser beam with the telescope axis, and optionally includes beam expanders to decrease the emitted beam divergence and for eye safety reasons. The depolarizing effect of beam expanders is neglected, since we cannot estimate the uncertainties introduced by possible birefringence as in the case of CaF2 lenses of apochromatic beam expanders. In addition, stress birefringence in windows in the transmitting part (e.g. roof window) is neglected due to its complex analysis and lack of information. The effect of these optical devices has to be investigated in the future. The general Müller matrix of dichroic beam splitters and steering mirrors is (Lu and Chipman, 1996;Freudenthaler, 2016a) This functional block, formed by the telescope and dichroic beam splitters, leads the received signal to the photomultipliers and, in case of multiwavelength lidar, separates the received signal by wavelength. In the same way as the emitting optics, M O can be described by a unique effective Müller matrix as follows: where T O , D O , and O are the effective transmittance, diattenuation, and retardance of M O , respectively, and γ describes the rotational misalignment angle of M O with respect to the incident plane of the PBS. The polarization effects of telescopes with small incidence angles of the light beam are neglected in this work. This approximation is valid for Cassegrain telescopes, but possibly not for Newtonian telescopes with 90 • fold mirrors (Clark and Breckinridge, 2011;Di et al., 2015;Lam and Chipman, 2015) as in the case of the PollyXT lidars . This source of error has to be investigated further. Furthermore, windows commonly used to protect the telescope can also affect the polarization due to intrinsic and stress birefringence, but this is also not considered in the simulator since the effect is very difficult to evaluate due to the lack of information about the properties and the time-dependent behaviour.

Polarizing beam splitters (M R and M T )
A polarizing splitter separates the received signal into reflected and transmitted signals. Consequently, two Müller matrices are required to describe the reflection (M R ) and transmission (M T ) processes. Due to the similarity in the shape of the matrices, both matrices can be described by one matrix with the subscript S = {R, T } with M S as follows: where S is the phase shift of M S , T S , D S are defined by with T p S and T s S as the parallel p and perpendicular s transmittance/reflectance, respectively, and Z S = √ 1 − D S . Parallel and perpendicular polarization is defined with respect to the incidence plane of the polarizing beam splitter. We attribute any polarization-independent absorption to the channel gain (see Sect. 2.6) and get therewith Ideal PBS split light in the two orthogonally polarized components without crosstalk, which means R p = T s = 0. But real PBS always transmit a fraction of the perpendicular polarization and reflect part of parallel polarization. This effect of the crosstalk on the linear depolarization ratio has been studied for example by Álvarez et al. (2006), and Snels et al. (2009. As shown in Freudenthaler (2016a), it is relatively easy to reduce the crosstalk by adding polarizer sheet filters after the PBS. The Müller matrix of the cleaned PBS, M # S , is where D # S = D # T = 1, D # R = −1 and the superscript "#" indicates that the PBS is "cleaned".
We assume that the extinction ratio of the cleaning polarizing sheet filters is sufficiently small and that they can be oriented with an accuracy much better than ±5 • with respect to the PBS, and that therefore the resulting error of the D # S can be neglected (See Appendix A).

Rotator, R
The parallel polarized component of the emitted laser beam can be detected either in the transmitted or in the reflected path behind the PBS. This depends on the orientation of the PBS with respect to the laser polarization. We consider this by means of a rotator R(ψ) (Eq. 14) before the PBS (see Fig. 1). For = 90 • , the reflected and transmitted signals correspond to the parallel and perpendicular polarized components, respectively, and vice versa for = 0 • .

Photomultipliers, η R , η T
The reflected and transmitted signals are detected by the photomultipliers which perform the light-to-electrical signal conversion. They affect the depolarization measurements as, in general, different photomultipliers have different gains. Regarding the Stokes-Müller formalism, we define the optoelectronic gains η R and η T for the photomultiplier gains of the transmitted and reflected signals including all optical attenuation of the lidar system in the transmitted and reflected path that is independent of polarization. We set them equal to 1 since we only investigate the polarization-dependent errors of the lidar optics.

Calibrator, C
The calibrator allows the determination of the polarization effect of the optical devices located behind the calibrator and of the gain ratio of the photomultipliers (PMTs). For the lidar set-up in Fig. 1, the calibration factor η includes the effects of the polarizing splitter (M R and M T ) and photomultiplier gains, η R and η T : Different calibration methods have been proposed either using the theoretical value of molecular depolarization (e.g. Cairo et al., 1999) or using additional optical devices as halfwave plates or polarization filters (e.g. Álvarez et al., 2006;Snels et al., 2009;Freudenthaler et al., 2009). Particularly, the 90 • calibration method has been implemented within EARLINET lidar systems (e.g. Freudenthaler et al. 2009;Nemuc et al., 2013;Mamouri and Ansmann, 2014;Bravo-Aranda et al., 2015). This method uses two measurements, rotating the polarizing plane of the received signal at angles φ 1 and φ 2 around the nominal axial rotation ( ) with the constraint |φ 2 − φ 1 | = 90 • (e.g. φ 1 = 45 • ,φ 2 = −45 • around the measurement position = 0 • ). These rotations allow the equalization of the reflected and transmitted signals and thus, any difference between the reflected and transmitted signals is only due to the polarizing effects of the optical devices between the calibrator and the photomultipliers. The equalization of the reflected and transmitted signals can be made by a physical rotation of the receiving optics including the photomultipliers, by rotating a half-wave plate placed before the PBS, or by using a linear polarizing filter rotated accordingly (see Freudenthaler, 2016a, for more details). The measured calibration factor, η * , is where the two rotation angles, φ 1 and φ 2 , are written as x45 • + ε, with x = ±1 indicating the rotational direction and ε taking the deviation from the nominal rotation into account. Then, their geometric mean, is calculated since it is less influenced by ε than η * ( , x45 • + ε), as indicated by Freudenthaler et al. (2009).
It is possible to show (Freudenthaler, 2016a) that the analytical expressions of η * √ ± ( , ε), corresponding to a different experimental set-up, are always in the form is the combined influence of all the parameters which are not corrected by η. Hereafter, the correction factor, will be noted by f (α, . . .).

Reflected and transmitted signals, I R and I T
According to the Stokes-Müller formalism, the reflected (I R ) and transmitted (I T ) signals can be obtained by multiplying the laser beam Stokes vector (I L ) by the subsequent Müller matrices, which represent the different functional blocks, and the atmosphere F, which describes the scattering matrix for randomly oriented particles (van de Hulst, 1957;Mishchenko and Hovenier, 1995), by where F 11 is the backscatter coefficient and a is the polarization parameter. Both parameters are range-dependent, but since this dependency has no influence on our error calculation, we omitted it in the following. Please note that, instead of the polarization parameter a, different but equivalent expressions are used in other publications as described in more detail in Freudenthaler (2016a). Probably most known is the depolarization parameter d = 1 − a used in Gimmestadt (2008).
With that, the Stokes vector of the lidar and calibration measurements is described by where the first element of the Stokes vectors is the energy detected by the photomultipliers. Based on the Stokes-Müller formalism presented, the detected energy value depends on the 18 lidar parameters summarized in Table 1. However, only some of them are considered by PLS for the assessment of the error sources (see "error source" column in Table 1).

Volume linear depolarization ratio, δ
The volume linear depolarization ratio δ which contains the contributions of particles and air molecules, is directly related to the polarization parameter a by (Mishchenko and  Hovenier, 1995) For further details regarding mathematical expressions of depolarization parameters, see Cairo et al. (1999). δ can be retrieved from lidar measurements by the following general equation given by Freudenthalter (2016a): where the parameters G T , G R , H T , and H R are determined by solving the matrix multiplication of Eq. (24) and separating the measured energy I S into terms with and without the polarization parameter, a, as follows: (Further details are given by Freudenthaler, 2016.) and δ * ( ) is the reflected-to-transmitted signal ratio divided by the calibration factor, η (Eq. 23): , where η has to be derived from the measured gain ratio, η * √ ± (Eq. 26), estimating f (α, . . .) either from information already available (e.g. from technical specifications) or from additional measurements performed for this purpose (see Freudenthaler, 2016a). Without any information, we have to assume f (α, . . .) ∼ = 1. Therefore, a better characterization of the lidar system through G S , H S , and f (α, . . .) decreases the systematic errors of the depolarization measurements.

Polarimetric Lidar Simulator (PLS) performance
To quantify the systematic error of the linear depolarization ratio, the Polarimetric Lidar Simulator (PLS) is developed based on the matrix equations resulting from the theoretical framework given by Freudenthaler (2016a) (see Sect. 2). Because the parameters involved in the simulation are not always independent, we use a complete grid sampling of the parameter error space to calculate all possibly measured signals by means of numerical matrix calculation (Eqs. 28 and 29) and correct each measurement for the polarization effects using the assumed real parameter values and the correction Eqs. (26) and (31). The steps of the PLS workflow are shown in Fig. 2 and explained in detail below.
1. Creation of a parametric model: the parameters are noted by x 1 , . . ., x n in Fig. 2 and listed in Table 1, and their uncertainties are either taken from the technical specifications of the optical devices, or assumed with reasonable values. Particularly, for the beam-splitter properties, the reflectance and transmittance coefficients require additional calculations according to the splitter 3. Generation of simulated values: Gaussian or uniform distributions with a large number of random values (∼ 100 or larger) are commonly used in order to obtain a reliable δ. In our case, the parameter uncertainties are not related to random variations but to a lack of knowledge of the real value, and therefore a uniform distribution with a complete grid sampling is used since all the combinations are equality probable. However, a large number of simulations for each parameter would lead to an unmanageable quantity of combinations, wherefore we adjust the number of iteration as follows.
a. For the parameters α, D E , β, γ , and ε, we use only three values: b. For the parameters E , D 0 , 0 , D T , and D R we use values between x i − x i and x i + x i which a fixed step calculated to provide around 10 6 combinations: Especially for E and 0 , the steps are chosen dense enough to avoid aliasing effects.
4. Evaluation of the model for each x i,j combination and the atmospheric parameters, F 11 and δ r : to determine the influence of the atmospheric δ on the systematic error, this procedure is performed at δ r values of 0.004 and 0.45 as representative values of the minimum and maximum atmospheric δ: a. calibration and measurement signals I i,j s ( ) and I i,j s ( , φ); b. calibration signals used to retrieve the calibration factor, η * ,j √ ± ; c. δ * j retrieved using η * ,j √ ± and f (α, . . .); d. δ s,j retrieved by means of Eq. (31).

Coloured squares highlight the workflow linked to
the calibration (red) and to the corrections performed (green) using the assumed real parameter values.
6. The analysis of the results is performed in three different ways.
a. The uncertainty propagation of each simulation parameter, x i , is analysed through the simulated-to-real δ difference, E δ (x i ) = δ s,j − δ r , varying x i within its uncertainty range [x i − x i , x i + x i ], while all the other parameters are kept. This method is used in Sect. 4.
b. The systematic errors of our lidar set-ups are not "randomly distributed", but are in a fixed state at the time of the measurements. We do not know this fixed state, but assume another fixed set of parameter values (assumed to be real) to correct the systematic errors. Therefore we cannot apply a standard error propagation as for randomly distributed error sources, but determine the depolarization uncertainty δ from the minimum and maximum of the simulated δ values: min δ i,j s , max δ i,j s . www.atmos-meas-tech.net/9/4935/2016/  , where a displacement towards being larger/smaller than δ r indicates an over-/underestimation of δ.

Depolarization uncertainties of the synthetic lidar
In this section, we define a synthetic lidar system in order to evaluate the error caused by each functional block. We choose the 90 • calibration method implemented by means of a mechanical rotator in front of the PBS. The synthetic lidar system is based on technical specifications of commercial optical devices (Table 2). It is worth noting the large uncertainty of the effective retardance of ±180 • , since this parameter is usually not specified at all by manufacturers and can take any value within this range.
The effect of each simulation parameter x i is analysed through the difference between the real and the simulated δ E δ (x i ) = δ s,j − δ r , varying x i within its uncertainty range [x i − x i , x i + x i ], while all the other parameters are kept. This is done for the two real linear depolarization ratios δ r = 0.004 and 0.45 (Figs. 4-7). Finally, the total systematic error δ is determined by max (E δ )−min (E δ ) from the complete set of δ simulations (Fig. 8).

Synthetic lidar: influence of the laser
The laser can introduce errors in the depolarization measurements due to a misalignment angle α of the plane of the laser polarization with respect to the incident plane of the PBS (see Eq. 1). Figure 3 shows E δ over α with different values of δ r . As can be seen in Fig. 3, E δ (α) increases with α in absolute terms. δ (α) ranges between [0, 0.031] and [0, 0.024] for δ values of 0.004 and 0.45, respectively, showing a low δ dependence. Since the smallest atmospheric δ m is of the order of 10 −3 , we assume that δ (α) can be neglected (< 1 × 10 −4 ) if α is in the range 0 • ± 0.6 • .

Synthetic lidar: influence of the emitter optics
M E is characterized by the effective diattenuation D E , the retardance E , and the rotational misalignment angle β of M E with respect to the incident plane of PBS. These parameters are interdependent, and thus Fig. 4 shows the E δ dependence on E and D E for different β. Additionally, the influence of the atmospheric depolarization is assessed for δ r values: 0.004 and 0.45 in Fig. 4 (top). As aforementioned, E varies in the range [−180 • , 180 • ] because the retardance of steering and dichroic mirrors is generally not provided in the majority of the technical specifications. As can be seen in Fig. 4 (top), D E introduces systematic errors, except for β = 0 • with δ (D E , β) = [0, 0.001]. Figure 4 (bottom) shows that the maximum of E δ ( E , β) is larger than 0.03 for β larger than 5 • in absolute terms, indicating that the lack of information of E can lead to huge uncertainties even larger than 0.13 for β larger than 10 • . Figure 4 shows also that E δ (D E , β) and E δ ( E , β) decrease with increasing δ values.
In summary, the total δ due to the systematic uncertainty of M E is [0, 0.13] for δ r = 0.004 and [0, 0.1] for δ r = 0.45. It is recommended that the laser beam is emitted directly to the atmosphere to avoid this error source. Otherwise, it is crucial to set β = 0 • to keep δ (M E ) as low as possible, independently of the effective diattenuation and retardance. For the synthetic lidar set-up, β = 0 • ± 2.5 • would lead to negligible E δ (D E , β) and E δ ( E , β) with values lower than 1 × 10 −4 .

Synthetic lidar: influence of the receiver optics
The parameters of the receiving optics (M O ) are the effective diattenuation D O , retardance O , and the misalignment angle γ between the receiving optics and the incident plane of the PBS. As in the case of M E (β, D E , E ), the influence of any of these parameters on E δ is not independent of the others. However, the variation of E δ with γ is very weak, and thus 3)| = 0.12 considering δ r = 0.25 in both cases) because the parallel signal is stronger than the perpendicular one. In order to make the systematic error caused by M O negligible, the uncertainty of D o should be lower than ±0.0010 for E δ (D O ) <10 −4 , which is a very small value compared to what is available on the market. Thus, we advise the use of accurate calibration methods which can correct for D O , or to measure this value accurately. Figure 5 (bottom) shows that the effect of the uncertainty in the retardance E δ ( O )is comparable to that of the emitter block E δ ( E ) with E δ ( O ) larger than 0.03 for γ = ±5 • . Therefore, it is recommended that γ is kept as small as possible with little uncertainty.
Summarising, the uncertainties of γ and D O have a huge impact on δ, and thus, it is very important to carefully determine these parameters of the receiving optics.

Synthetic lidar: influence of the polarization splitter
For the synthetic lidar, we consider a non-cleaned polarizing beam splitter which has T p T , T p R , T s T , and T s R values and uncertainties as shown in Table 2 The ratio E δ T s T / T s T is about 8 times the ratio E δ T p T / T p T (i) because the parallel intensity, detected in the reflected path, is much larger than the perpendicular intensity, and (ii) because T s T is small; a small change in T s T makes a big difference in the crosstalk from the parallel to the perpendicular signal. The influence of T p T on E δ T p T , T s T increases and that of T s T decreases with increasing δ. The systematic errors due to crosstalk can be avoided by means of a cleaned PBS (see Sect. 2.4).

Synthetic lidar: influence of the calibrator
Here the 90 • calibration with a rotator in front of the PBS is considered. This calibration method is described in Sect. 2.7 (see Freudenthaler, 2016a, for further details). It is worth noting that the uncertainties due to a rotation calibrator affect both calibration and normal measurements, whereas the use of a polarizing filter as calibrator only affects the calibration measurements because it has to be removed for the normal measurement, and hence there is no rotation error ε in that.
As can be seen in Fig. 7, E δ (ε) is [0, 0.0008] for δ r = 0.004 and [0, 0.0006] for δ r = 0.45. Considering the large uncertainty range (ε = ±5 • ), δ (ε) is really low. Since ε is usually smaller in real lidars than in the synthetic lidar, we can conclude that δ (ε) can be neglected. This happens because the 90 • calibration method is used, but for other calibration methods, the uncertainty of this parameter needs to be taken into account.
Despite the fact that ε does not introduce a considerable systematic error, the position of the calibrator within the lidar system has to be considered. For example, the 90 • calibration in front of the receiving optics corrects the calibration factor for diattenuation of the receiving optics D O but it does not if the calibrator is located behind the receiving optics, as is the case for the LB21, IPRAL, and MUSA lidars (see Table 5).

Synthetic lidar: total uncertainty analysis
The total systematic error δ, including all possible mutual dependencies, is determined with a complete grid search of the error space. In order to keep the number of samples around 10 6 , α, D E , β, γ , and ε are simulated with three values, as e.g. α = −10, 0, 10  Table 3. For both simulations, the values obtained for δ span over quite a large range reaching even impossible values below the molecular δ (0.004) or δ values larger than 1. Since the distribution of simulated δ is displaced to the right of δ r , an overestimation of δ is more probable than an underestimation for this synthetic lidar.

Systematic depolarization errors of seven EARLINET lidar systems
The PLS is applied to the seven EARLINET lidar systems listed in Table 4, of which POLIS measures the linear depolarization ratio at two wavelengths. Detailed information about the lidar systems is given by Wandinger et al. (2016), except for the upgraded POLIS (Freudenthaler et al., 2016b) and IPRAL (IPSL high-Performance multi-wavelength RAman Lidar for Cloud Aerosol Water Vapor Research) which has recently been deployed at SIRTA atmospheric research observatory (Haeffelin et al., 2005). IPRAL provides measurements at 355 (parallel and perpendicular polarized components), 532, and 1064 nm (elastic backscatter), and the Raman-shifted backscatter at 387 (from N 2 ), 408 (from H 2 O), and 607 nm (from N 2 ). Table 5 shows the values and uncertainties of the lidar parameters used for the simulation. The main differences among the lidars are (i) the use of steering optics in the emitter (e.g. MUSA and POLIS do not have this functional   ∼ 1 × 10 6 . Finally, for IPRAL, the uncertainties of O , D T , and D R are neglected, and thus α, D E , β, γ , and ε are simulated with three values, whereas D O and E are simulated with 65 values, resulting in 3 5 65 2 ∼ 1 × 10 6 combinations. Figures 9 and 10 show the histograms of the simulated δ at δ r = 0.004 and δ r = 0.45 for the EARLINET lidars in Table 4. The discontinued distributions for LB21, POLLY-XT SEA, and POLIS show that certain parameters with a big impact are under-sampled. Disregarding the discontinuities, the distributions in the histograms are very different between the lidar systems; some are almost centred Gaussian-like, some more top-hat, and others one-sided, skewed, and strongly peaked. Therefore we only list the minima and maxima of the simulated δ values for each lidar in Table 6.
With the plots in Fig. 11, it is possible to show the contribution of individual parameters to the histograms as coloured sub-histograms. For the case of LB21, the diattenuation parameter of the receiving optics D O is clearly the dominant error source, and the spread of each sub-histogram shows the combined contribution of the other parameters. In the plot for MULHACEN, with coloured sub-histograms for the three values of the rotation angle of the laser polarization α, we -0.500 0.015 0.225 -0.8 see that the shift between the sub-histograms due to the uncertainty of α accounts for about one-third of the total width of the histogram. The RALI plot, with coloured ε-histograms indistinguishably lying above each other, shows that the uncertainty of ε has a negligible impact on the total error. The small error of POLIS compared to the other lidars can be attributed to (i) the custom-made dichroic mirrors designed for negligible diattenuation, (ii) the absence of any emitter optics, (iii) the cleaned PBS, and (iv) the well-adjusted rotation of the laser polarization with small uncertainty (Table 5).
Despite δ still being large for some lidars, the comparison with the much larger error of the synthetic lidar shows the success of the effort carried out by the EARLINET community towards the identification of the main error sources and a better characterization of the individual optical elements. From the results of our study, some indications can be provided to reduce the polarization-dependent systematic errors.  The polarization purity of the laser beam can be improved by using a high-energy polarizer after the emitting optics. To reduce the uncertainty introduced by M E , it is highly recommended that emitter optics is avoided. The errors due to the emitter and receiver optics can be reduced by improving their rotational alignment (i.e. β and γ ) with respect to the polarizing beam splitter. Finally, the PBS cross-talk can be removed by using a cleaned PBS. Additionally, a good characterization of the parameters and consequential correction can drastically reduce the systematic error. For example, the α and D O values can be determined by experimental assessments, as indicated by Freudenthaler (2016a).

Conclusions
This work shows the numerical analysis of the polarizationrelated systematic errors of the linear depolarization ratio δ of one synthetic and seven EARLINET lidar systems, which all use one of the 90 calibration techniques described by Freudenthaler (2016a). It uses the lidar model described there to apply corrections to the signals and calibrations with assumed real parameter values. For the error estimation we use a complete grid scan over all relevant error sources.
From this simulation we conclude that the systematic error can be very large if the lidar system is not well characterized and aligned.
Within the ranges of the typical uncertainties, the diattenuation of the receiving optics and of the polarizing beam splitter have the biggest impact. The next important parameters are the retardance of the emitting and of the receiving optics, and the rotational misalignment between the plane of polarization of the laser and the incident plane of the PBS.
The simulations of the EARLINET lidars show possible measured values of δ in the range of [−0.012, 0.039] for a δ r of 0.004 and in the range of [0.399, 0.512] for a δ r of 0.45. Compared to the large error of the synthetic lidar that is not well-characterized, this is a big improvement.
The uncertainty of some parameters, as the retardance of dichroic mirrors, is still very large, often because the manufacturers of optics do not provide specific information. The numerical error analysis shown in this work is actually a sensitivity analysis, which can be used to identify the parameters that need more accurate characterization. A positive example for such an improvement is the receiving optics of IPRAL, which employs custom-made dichroic beam splitters with almost negligible diattenuation and 0 • retardance.
Finally, further investigations are still required for a better understanding of the polarization effects of large windows, special lenses, and Newtonian telescopes. Possible evidence of error sources not considered in this work is if δ values are measured outside of the simulated δ distribution. For example, elliptical polarization of the outgoing laser beam could strongly affect the δ determination.

Data availability
The polarimetric lidar simulator and simulations supporting this article are available upon request from the corresponding author.