Complex refractive index of crystalline quartz particles from UV to thermal infrared

Abstract

Crystalline quartz is an important element of many materials and one of the major compounds of mineral dust. Thus, knowledge of its optical properties is essential for many applications, in particular for remote sensing techniques. However, despite the fact that crystalline quartz has been an object of study for many years using various experimental approaches, its optical properties and in particular its complex refractive indices (CRI) remain uncertain. The focus of this investigation is the determination of a new set of CRI of crystalline quartz particles retrieved from a methodology which has already demonstrated its efficiency on pure amorphous materials and volcanic ashes. The experimental set-up allows recording simultaneously size distribution (from 15 nm to 20 µm) of airborne particles and extinction spectra over a wide continuous spectral range from thermal infrared (650 cm⁻¹) to UV (40,000 cm⁻¹) with high spectral resolution (up to 0.5 cm⁻¹). By associating these measurements and a numerical procedure coupling scattering theories, the single subtractive Kramers-Kronig relation, and an optimal estimation method in an iterative process, we precisely determine the real and imaginary parts of the CRI of crystalline SiO₂ with mean values of the uncertainties of 2.5% and 1.7%, respectively. Moreover, the use of these CRIs allows to simulate an extinction spectrum of another quartz sample very well, which shows their suitability as well as the good reproducibility of the whole process.

EDITOR:

• Hans Moosmüller

1. Introduction

Silicon dioxide (SiO₂) exists in different forms: crystalline, amorphous and glassy. Due to its optical and physico-chemical characteristics, it is an important component in many technical and scientific applications (thermal, electronic or optical). It is naturally present in many different astrophysical environments (Petrov, Zhuzhulina, and Savushkin 2021 and references therein), and it is one of the main oxides in the Earth’s crust (Heaney, Prewitt, and Gibbs 1994). It is therefore not surprising that silicon dioxide is one of the major components of atmospheric mineral dust emitted during desert storms and volcanic eruptions. Given the large quantities annually injected into the atmosphere, mineral dust aerosols play an important role in the radiative budget equilibrium (Granados-Muñoz et al. 2019; Saidou Chaibou, Ma, and Sha 2020) and atmospheric chemistry (Usher, Michel, and Grassian 2003). Accurate modeling of the effect of mineral dust aerosol on the different atmospheric processes requires information about its chemical/mineralogical composition, concentration, size distribution (SD), or height and thickness of the aerosol plume. These information can be obtained using different active (Sugimoto and Huang 2014; Winker et al. 2013) or passive (Alalam et al. 2022; Clarisse et al. 2010; DeSouza-Machado et al. 2006; Ginoux et al. 2007; Wang et al. 2019) remote sensing techniques. In particular, due to very strong resonance absorption peaks, mainly from the Si-O vibrations which leads to a typical spectral signature between 8 and 12 µm, mineral dusts (e.g., Ackerman 1997; DeSouza-Machado et al. 2010; Pierangelo et al. 2004; Klüser et al. 2012) as well as volcanic ash (Prata 1989; Grainger et al. 2013) can be detected in the thermal infrared (TIR) region.

However, regardless of the measurement method or spectral range used, it is essential to have a thorough knowledge of their optical properties. Indeed, the extinction of radiation produced by aerosol particles depends mainly on the particles size distribution, and complex refractive indices m($\overline{\nu })$ = n($\overline{\nu })$ + i $\kappa (\overline{\nu }),$ where $\overline{\nu }$ denotes spectral wavenumber. The predominant part of the uncertainty in modeling the effect of aerosols on atmospheric radiative transfer arises from the lack of knowledge of these parameters, which requires assumptions to be made about them leading to uncertainties in remote sensing and radiative forcing estimates (Boucher et al. 2013).

While the optical properties of amorphous silica seem to be well understood (Hubert et al. 2017), this is not the case for crystalline silica (Quartz) for which the latter, in the infrared, remain largely uncertain. The optical constants of crystalline SiO₂ (Longtin et al. 1988 and references therein) were obtained from two different methods: (1)(1) $$D_m=D_a\sqrt{\chi \frac{C\left(D_a\right){\rho }_0}{C\left(D_m\right){\rho }_m}}$$(1) using KBr pellet sample (Steyer, Day, and Huffman 1974), (2)(2) $$\frac{dN}{d{\mathit{LogD}}_{}}=\frac{N}{\sqrt{2\pi }\mathrm{Log}\sigma }\exp \hspace{0.17em}\left[-\frac{1}{2}{\left(\frac{\mathrm{Log}\left(\raisebox{1ex}{$D$}\!\left/ \!\raisebox{-1ex}{$\overline{D}$}\right.\right)}{\mathrm{Log}\sigma }\right)}^2\right]$$(2) exploiting reflectance spectra and dispersion theory of solid crystal (Spitzer and Kleinman 1961) coupled with emissivity spectra and Fresnel’s law (Wenrich and Christensen 1996). Furthermore, it is now accepted that the optical constants obtained from these techniques are not suitable for aerosols (McPeath et al. 2002) on the one hand because the use of bulk material underestimates the scattering effect and on the other hand because the use of pellets modifies the size, the shape and the vibrational modes. From an overview of the existing literature, it is impossible to obtain consistent values for the optical constants on broad spectral range with a spectral resolution suitable for the latest remote sensing instruments. Previous work has shown that extinction spectra of quartz cannot be simulated from the literature indices, either using the spectral averaging method or the optical constant averaging method (Mogili et al. 2007). In this way, some studies have attempted to investigate the effects of particle shape for quartz such as: disk, needle, or ellipsoids associated with different scattering theories such as Continuous Distribution of Ellipsoids (CDE) (Hudson et al. 2008), or T-matrix (Kleiber et al. 2009; Reed et al. 2017). However, even if the comparisons between experimental and calculated spectra are better for the CDE and T-matrix methods than Mie theory, they remain very unsatisfactory in particular in the thermal infrared resonance peak region.

The present work aims to provide a set of consistent refractive indices of crystalline quartz over a wide spectral range (from IR to UV) with a spectral resolution adapted to the capacities of the most recent remote sensing techniques. The next section details the experimental setup and protocol. The size distribution determination and the extinction spectrum are also exposed. Section 3 describes the new aspects of the numerical methodology compared to that used previously (Herbin et al. 2017; Hubert et al. 2017, Deguine et al. 2020) and the newly retrieved refractive indices and their uncertainties are discussed. Finally, Section 4 summarizes our results and presents perspectives for future applications.

2. Experimental setup and results

2.1. Aerosol generation

The aerosol suspension process has already been presented in previous articles (Hubert et al. 2017, Deguine et al. 2020) and is only briefly recalled here. The quartz sample (S1) comes from Sigma-Aldrich (CAS number 14808-60-7, particle size < 63 μm) and has a density provided by the manufacturer of 2.65 g.cm⁻³. The sample in a petri dish is placed in an oven at 110 °C for 12 h to remove any water residues. Then, the particles in powder form are introduced into a glass reactor where the fine particles are resuspended and carried by a flow of nitrogen gas (N₂, purity 99.996%), under mechanical agitation from magnetic stirring. The concentration can be adjusted by varying the speed of the magnetic stirrer between 0 and 1500 rpm and the inlet flow rate from 4.3 to 6.5 L.min⁻¹ managed by a mass flow controller. A buffer volume is placed at the outlet of the reactor to decouple the resuspension and sampling parts. This volume also ensures a homogeneous mixture of aerosols and avoids overpressure in the device thanks to an exhaust. At fixed speed rotation and N₂ flow rate, the concentration remains stable over several minutes, allowing good quality recordings during the acquisition time of the spectra. Finally, tests carried out by placing particle counters upstream and then downstream of the spectrometers show that the loss in number of particles in the cells is negligible (<1%).

2.2. Size distribution determination

In order to accurately determine the entire size distribution, we use two instruments simultaneously. The first is an Aerodynamic Particle Sizer (APS, TSI 3321) measuring the time of flight of the particles, allowing the determination of their aerodynamic diameter (D_a) in the range 0.523 − 20 μm. A measured size distribution by the APS (32 channels per decade, 52 channels in total) of a quartz sample is represented by black bins on Figure 1.

Figure 1. Quartz particles size distribution measured by APS (black bins), SMPS (grey bins), and the corresponding lognormal function with $\overline{\mathit{D}}$ = 0.24 µm and σ = 2.24 in red curve.

The second device is a Scanning Mobility Particle Sizer (SMPS TSI) used to measure finer particles. This instrument couples a Differential Mobility Analyzer (DMA TSI 3080) to select particles according to their electrical mobility and a downstream Condensation Particle Counter (CPC TSI 3775) that counts the previously selected particles. The SMPS records the particle size distribution with an electrical mobility diameter D_m between 10 and 800 nm. The SMPS size distribution measurement is given by grey bins on Figure 1.

The two size distributions from SMPS and APS overlap between the 500 nm and 800 nm geometric diameters where both instruments have good counting efficiencies. As detailed in Khlystov 2004, by measuring the same size distribution based on different physical properties it is possible to estimate the dynamic shape factor χ, according to the equation (Hinds 1999): (1) $$D_m=D_a\sqrt{\chi \frac{C\left(D_a\right){\rho }_0}{C\left(D_m\right){\rho }_m}}$$(1) where ρ₀ and ρ_m are densities of pure water (1 g.cm⁻³) and particle under study, respectively, C(D_a) and C(D_m) are the Cunningham slip factors for the aerodynamic and mobility diameters, respectively, according to Kim et al. (2005).

By using a similar method as described by Khlystov 2004, which consists in determining χ by fitting the Equation (1)(1) $$D_m=D_a\sqrt{\chi \frac{C\left(D_a\right){\rho }_0}{C\left(D_m\right){\rho }_m}}$$(1) from the experimental data D_m from the SMPS and D_a from the APS in the area of diameters overlap, we estimated the shape factor χ = 1.33. This value is typical for irregular particles and is consistent with these obtained by Alexander et al. 2016 for similar particle size ranges. By taking χ into account, the aerodynamic diameters D_a measured by the APS are converted into electric mobility diameters D_m.

Both distributions can then be fitted with one lognormal function, which is symmetrical with respect to the mean diameter described as: (2) $$\frac{dN}{d{\mathit{LogD}}_{}}=\frac{N}{\sqrt{2\pi }\mathrm{Log}\sigma }\exp \hspace{0.17em}\left[-\frac{1}{2}{\left(\frac{\mathrm{Log}\left(\raisebox{1ex}{$D$}\!\left/ \!\raisebox{-1ex}{$\overline{D}$}\right.\right)}{\mathrm{Log}\sigma }\right)}^2\right]$$(2) where N is the concentration in particles per cm⁻³, D is the diameter of the size distribution in μm, $\overline{D}$ is the mean diameter, and σ is the geometric standard deviation. Thus, this function generates 182 (107 identical to SMPS and 75 to APS) size classes starting from a diameter of 14.6 nm up to 10.2 μm and where the base-10 logarithms of the diameters are regularly spaced (see red curve on Figure 1).

2.3. Extinction measurements

The infrared part of the spectrum is recorded from a Fourier Transform InfraRed (FTIR) Spectrometer (iS50 from Thermo Scientific) equipped with a multi-pass cell with ZnSe (TIR) and BaF₂ (NIR) windows, and an optical path of 10 meters which allows to increase the sensitivity, i.e., to reduce the particle concentration in order to verify the single scattering criterion. Moreover, the cell is placed horizontally in order to minimize the deposition of particles on the mirrors. The spectrometer is also equipped with a set of beam splitters (KBr, CaF2) and detectors (MCT, InGaAs) to cover the 600–9000 cm⁻¹ spectral range with a spectral resolution of 0.5 cm⁻¹.

The UV-Visible part of the spectrum is obtained using a spectrometer composed of a source containing two specific lamps: a deuterium which emits in the UV (210–400 nm) and a halogen one covering the visible and near infrared parts (360 and 1500 nm). The light beam is injected into an optical fiber with a core diameter of 200 μm. At the output of the fiber a collimator delivers a parallel beam of 8 mm diameter through the cell. At the other end, the beam is collimated onto the input of an optical fiber which is connected to the Maya2000-Pro detector device (Ocean Optics). The latter consists of a grating (HC-1) and of a 2064 × 70 pixel CCD (Hamamatsu S10420) which collects the diffracted light in a spectral band between 200 and 1100 nm (i.e., 9000 and 40,000 cm⁻¹). A slit placed at the entrance of the device, with an opening of 25 μm, allows spectra to be obtained with a spectral resolution of 1 nm. The spectrometer has a single pass cell of 1 meter optical path and 48 mm inner diameter. Aerosols are injected at 45° into the cell to limit the accumulation of particles on the inner surfaces, and the whole is inside an opaque black box to limit parasite light.

The measured extinction spectrum of quartz corresponding to the size distribution discussed in Section 2.2 is shown in Figure 2. The latter is a combination of the three spectral ranges covered by the MCT (620–4000 cm⁻¹), the InGaAs (3500 − 9000 cm⁻¹) and CCD (9000–40,000 cm⁻¹) detectors. The spectrum obtained is very close to those already obtained previously (Mogili et al. 2007; Reed et al. 2017). However, we can notice the good signal-to-noise ratio, with mean values for each spectral range covered being 200 between 620–4000 cm⁻¹, 400 in the 3500–9000 cm⁻¹ spectral range and 100 in the UV-visible region. Despite the purge of several hours, two rovibrational bands of residual water vapors are present around 1600 and 3600 cm⁻¹. In the infrared region between 650 and 1300 cm⁻¹, three absorption bands caused by vibrations of the silicon-oxygen bonds are observed. The main absorption band is located around 1100 cm⁻¹, and is attributed to the asymmetric elongation vibrations of the Si-O-Si bridges (Kitamura, Pilon, and Jonasz 2007). The absorption band between 782 and 802 cm⁻¹ is due to Si-O-Si symmetric elongation vibrations, involving the displacement of the oxygen atom perpendicular to the plane in which the silicon atoms lie. Finally, the band at 696 cm⁻¹ corresponds to the angular deformation of the Si-O-Si bridge. The positions of each peak determined in our extinction spectra are in good agreement with literature values (Hlavay et al. 1978; He et al. 2019).

Figure 2. Extinction spectrum of quartz between 650 and 40,000 cm⁻¹. The spectral range 8500–9500 cm⁻¹ is very noisy, the plot in this gap corresponds then to linear interpolated data. The two strong rovibrational absorption bands around 1600 and 3600 cm⁻¹ are due to residual water vapor and they are therefore removed in the CRI calculation. The inset shows a zoomed portion of the extinction spectrum in the TIR (between 650 and 1350 cm⁻¹) where the o- and e-rays are pointed out by vertical lines.

Quartz is a birefringent material whose optic axis correspond to the crystallographic axis (a, b and c). Electromagnetic wave propagating in any direction through a uniaxial mineral such as quartz is divided into two polarized rays, the ordinary (o) ray and the extraordinary (e) ray, traveling at different velocities. So, the o-ray energy can propagate along both the c and the a axes, however, the e-ray energy may not propagate along the c axis because the electric vector of the extraordinary ray is parallel to the c axis. Thus, the energy traveling parallel to the c axis is strictly the o-ray and the energy traveling parallel to the a axis is composed of equivalent amounts of o- and e-rays. Then, the accepted method for calculating the optical properties for a distribution of randomly oriented quartz particles is to calculate individual spectra for both o-ray and e-ray indices, and then to take a weighted average of the spectra in the form [(2/3) o-ray + (1/3) e-ray] (Bohren and Huffman, 1983). Due to the experimental method used for particle suspension and transportation, the spectrum in Figure 2 corresponds to an average over different particles orientations, which we checked using a polarizer. Indeed, no changes of the intensity in the extinction spectrum were detected when using an IR polarizer at different orientations. It is also interesting to notice that often, the tabulated optical constants are at fairly low resolution and the fine details of the resonance line shape may not be well modeled, but here thanks to a spectral resolution of 0.5 cm⁻¹ the absorption bands of o- and e-rays are quite well resolved in particular around 800 cm⁻¹.

3. Complex refractive indices

3.1. Methodology

The formulation of the extinction spectrum (cf. Figure 2) can be written under the single scattering criterion: (3) $$E\left(\overline{\nu }\right)=-\log \hspace{0.17em}\frac{I\left(\overline{\nu }\right)}{I_0\left(\overline{\nu }\right)}={\int }_{r_{\min }}^{r_{\max }}l\mathrm{ }N\left(r\right){\mathrm{ }\sigma }_{\mathit{ext}}(r,\overline{\nu })$$(3) where E($\overline{\nu }$) is the measured extinction spectrum, l is the optical path length (cm), N(r) is the number concentration of particles with radius r (part.cm⁻³), and $\sigma$_ext is the extinction cross-section (cm²).

The formulation of $\sigma$_ext depends on the scattering regime, but whatever its form it always depends on the size distribution and the real and imaginary parts of the refractive index. The simple idea is to use experimental measurements of extinction and size distribution to determine the CRIs.

For the extinction spectrum, the only manipulation of the experimental data is to smooth the UV-vis. spectral range to remove high frequency noise and, in the IR region the residual water vapor bands to avoid including them in the calculation of the indices. Then, the whole spectrum is interpolated onto a grid of wavelengths regularly spaced by 1 cm⁻¹.

From an optical point of view, we have to consider the volume equivalent diameter D_ve of the particles, we must therefore convert the size distribution in Figure 1 from the following relationship: (4) $$D_{ve}=D_m\frac{C(D_{ve)}}{\chi C\left(D_m\right)}$$(4) with C(D_ve) the Cunningham slip factor of volume equivalent diameter.

From Equation (4)(4) $$D_{ve}=D_m\frac{C(D_{ve)}}{\chi C\left(D_m\right)}$$(4) it can be seen that for irregular particles (χ > 1), the mobility diameter (D_m) is always larger than the volume equivalent diameter (D_ve). As we do not know the value of C(D_ve) and the value of the shape factor χ is only an average estimation for the whole SD with a large uncertainty, the lognormal size distribution parameters are adjusted using a least squares method close to that described in Mogili et al. 2007, but only in the spectral range 4000–40,000 cm⁻¹, where the extinction spectrum is dominated by the real part of the CRI, n, whose values are correctly retrieved regardless of the experimental method used (pellets, bulk material, suspended particles). The new SD parameters values are D_ve = 0.18 µm and σ = 2.22 and the simulated extinction spectrum using these fitted parameters is compared with the experimental one in Figure 3. As expected, the Mie theory using the indices from the literature reproduces quite well the experimental extinction spectrum on the spectral region: 4000–40,000 cm⁻¹.

Figure 3. Comparison between experimental extinction spectra of quartz (grey line) in the NIR and UV-visible regions, calculated from Mie theory using the CRIs from the literature (black line) and calculated from geometric optic (grey dashed line).

In addition, a spectrum calculated from geometric optics has been performed because at high wavenumber the extinction mainly depends on the SD. Indeed, the asymptotic representation of the Mie series gives the following approximate expression for the extinction efficiency (Nussenzveig and Wiscombe 1980): (5) $$Q_{e\mathrm{x}t}=2(1+x^{-\raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$3$}\right.})$$(5) where x = $\raisebox{1ex}{$2\pi r$}\!\left/ \!\raisebox{-1ex}{$\lambda $}\right.$ is the size parameter, with r the radius of the particles and λ the wavelength, both in µm.

The Equation (5)(5) $$Q_{e\mathrm{x}t}=2(1+x^{-\raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$3$}\right.})$$(5) shows that the extinction of micrometric particles does not depend on the refractive index values in the UV-Visible part. As shown in Figure 3, the coherence between the three curves for $\overline{\nu }$ > 30,000 cm⁻¹ validates the SD parameters retrieved

In practice, even knowing extinction and size distribution, obtaining CRIs is not trivial, because at each wavenumber corresponds two unknowns (n and κ). The problem is therefore mathematically ill-conditioned and under-constrained. We therefore use the Optimal Estimation Method (OEM) (Rodgers, 2000) which is perfectly adapted to weakly nonlinear problems and widely used in the remote sensing community to retrieve physical parameters from experimental measured spectra. Moreover, the use of Kramers-Kronig relations which link the real and imaginary parts of the CRIs adds a correlation between n and κ which allows on the one hand to increase the constraints on the restitution but also to ensure the respect of the principle of causality for the couple n and κ at the end of the iterative process on the whole covered spectral domain. In short, the goal of the OEM is to obtain a meaningful solution by regularizing the retrieval with a priori information about the variables. This information is composed of an a priori vector ${\mathbf{x}}_a$ formed by n_a $({\overline{\nu }}_a)$ and κ_a$({\overline{\nu }}_a)$ values, and an a priori variance-covariance matrix ${\mathbf{S}}_a$ (Rodgers, 2000). In this case, the goal of the retrievalis to find the real and imaginary parts of the refractive index $\stackrel{\mathbf{⁁}}{\mathbf{x}},$ which are the approximation of the true state x, which best agrees with both the measurement and the a priori information. The solution $\stackrel{\mathbf{⁁}}{\mathbf{x}}$ is found by an iterative process (Rodgers, 2000): (6) $${\stackrel{\mathbf{⁁}}{\mathbf{x}}}_{i+1}={\mathbf{x}}_a+{({\mathbf{K}}_i^T{\mathbf{S}}_{\epsilon }^{-1}{\mathbf{K}}_i+{\mathbf{S}}_a^{-1})}^{-1}{\mathbf{K}}_i^T{\mathbf{S}}_{\epsilon }^{-1}\left[\mathbf{y}-\mathbf{F}({\stackrel{\mathbf{⁁}}{\mathbf{x}}}_i)+{\mathbf{K}}_i({\stackrel{\mathbf{⁁}}{\mathbf{x}}}_i-{\mathbf{x}}_a)\right],$$(6)

where y is the measured extinction spectrum, F is the forward model used (here it corresponds to the scattering theory), $K=\frac{\partial F}{\partial x}$ is the Jacobian matrix, and K^T is its transposed. Previously ${\mathbf{S}}_{\epsilon }{=\mathrm{\sigma }}_{\epsilon }^2\mathbf{I}$ (Hubert et al. 2017; Deguine et al. 2020), was estimated only from the signal-to-noise ratio. We have here another procedure which determines an average noise in the two spectral regions including both the instrumental noise and the error made by the smoothing of the experimental spectrum. Sε is thus a diagonal matrix so each element is calculated, such that: (7) $${\mathrm{\sigma }}_{\epsilon }^{}\left(\overline{\nu }\right)=\left|\frac{1}{2}\left(Natural\mathrm{ }exp.\mathrm{ }spec.\left(\overline{\nu }\right)-Smoothed\mathrm{ }exp.\mathrm{ }spec.\left(\overline{\nu }\right)\right)\right|$$(7)

The first step of the iterative process consists to determine an initial guess for the imaginary component, κ_a$({\overline{\nu }}_a),$ from the Rayleigh theory. The Single Subtractive Kramers-Kronig (SSKK) integration is then performed to calculate the real component, n_a $({\overline{\nu }}_a)$. The latter requires an anchor point, n_a (${\overline{\nu }}_a$), which is the real component of the refractive index at a wavenumber $({\overline{\nu }}_a)$ chosen in the UV-visible spectral region of the spectrum, where the sensitivity of the measurement is maximum for n and minimum for κ (see Herbin et al. 2017). S_a is estimated from the difference between the experimental spectrum and the calculated one using the n_a $({\overline{\nu }}_a)$ and κ_a$({\overline{\nu }}_a)$ estimations, typically set to 2%. These a priori x_a and S_a are used by the forward model F to simulate an extinction spectrum and modify at each iteration according to Equation (6)(6) $${\stackrel{\mathbf{⁁}}{\mathbf{x}}}_{i+1}={\mathbf{x}}_a+{({\mathbf{K}}_i^T{\mathbf{S}}_{\epsilon }^{-1}{\mathbf{K}}_i+{\mathbf{S}}_a^{-1})}^{-1}{\mathbf{K}}_i^T{\mathbf{S}}_{\epsilon }^{-1}\left[\mathbf{y}-\mathbf{F}({\stackrel{\mathbf{⁁}}{\mathbf{x}}}_i)+{\mathbf{K}}_i({\stackrel{\mathbf{⁁}}{\mathbf{x}}}_i-{\mathbf{x}}_a)\right],$$(6) to improve the comparison between measured and calculated spectra F(x).

The retrieved state is obtained after convergence, that is to say when there is no further improvement in the difference between the intensities of the calculated and measured spectra. This is generally the case when this difference is of the order of magnitude of the instrumental noise.

3.2. CRIs retrieval

In previous works we have used Mie theory to calculate the extinction. On the one hand, because Mie theory is well understood and easy to apply, given the size distribution and complex refractive indices of the particles (Bohren and Huffman, 1983), and on the other hand because it is classically used for modeling radiative transfer through the scattering atmosphere in both climate modeling calculations and satellite data retrieval (Bond et al. 2013; Nishizawa et al. 2008; Clarisse et al. 2010). The quartz retrieval using the method described in Section 3.1 give acceptable results in the NIR and UV-visible parts of the spectrum as suggested in the Section 3.1 and the Figure 3. However, the use of Mie theory for irregular particle shapes failed in the TIR spectral range. In practical terms, this means that despite the good quality of the measurements and the a priori (values and variabilities), the algorithm failed to converge satisfactory to reproduce the measured extinction spectrum in the thermal infrared region. This can be explained by the fact that here the shape factor χ is too significant to be neglected. Mie theory is derived only for homogeneous spherical particles and it is well known that mineral aerosols can be inhomogeneous and irregular in shape (Dick et al. 1998). This can have a significant impact on the optical properties (Kalashnikova and Sokolik 2004) and then introduces inaccuracies and errors in the aerosol characterization (Pierangelo et al. 2004).

To overcome this difficulty, some methods for calculating the scattering by nonspherical particles, such as T-matrix or discrete dipole approximation (DDA) (Drain and Flatau 1994; Mischenko et al. 2002) have been tested (Dubovik et al. 2006; Kalashnikova et al. 2005). These methods offer more flexibility in constructing model particles that are inhomogeneous and irregularly shaped, but they are more computationally demanding and they still involve significant approximations and assumptions. Thus, it is clear that further experimental studies are needed to assess the relevance of these different approaches, but at present it is still very difficult to accurately characterize the shape of the particles. For example, SEM appears to be an interesting tool, but it suffers from a lack of statistical representation of the sample and only gives information about the shape in two dimensions.

Furthermore, in the particular case of quartz, several works have tested the use of different methods to calculate the extinction (Mogili et al. 2007, Hudson et al. 2008, Reed et al. 2017) from literature CRI. In particular, Reed et al. (2017) have recently attempted to reproduce their quartz infrared spectrum from a minimization method using T-matrix by varying the form factor $\xi ,$ which offers more degrees of freedom. However, even using a wide form factor distribution, the fit remains unsatisfactory. This suggests that although the consideration of particle shape is important, it is not sufficient to explain the remaining discrepancies which are attributable to incorrect refractive index values.

Previous studies (Hudson et al. 2008; Reed et al. 2017) have shown that the continuous distribution of ellipsoids (CDE) model, corresponding to the case where all forms are equiprobable, better reproduce the resonance peak positions and band shapes than Mie theory and provides remarkably good fits to the infrared extinction data for some mineral components (especially calcite and quartz), it therefore seems appropriate to replace in our algorithm the Mie extinction cross section calculation by the CDE expression (Bohren and Huffman, 1983): (8) $${\mathit{\sigma }}_{\mathit{CDE}}^{\mathit{ext}}=\mathit{kV}\text{Im}\left[\frac{2{\mathit{m}}^2}{{\mathit{m}}^2-1}\ln \hspace{0.17em}\left({\mathit{m}}^2\right)\right]+\frac{{\mathit{k}}^4{\mathit{V}}^2|{\mathit{m}}^2-1{|}^2}{3\mathit{\pi }\text{Im}\left({\mathit{m}}^2\right)}\text{Im}\left[\frac{{\mathit{m}}^2}{{\mathit{m}}^2-1}\ln \hspace{0.17em}\left({\mathit{m}}^2\right)\right]$$(8) where k = 2 π/λ, and V = 4/3 π$r_{ve}$³ is the equivalent volume of particles.

However, expression 7 is only valid in the Rayleigh approximation, i.e., for size parameters x≪1, which is incompatible with our method of retrieving the indices, which requires the use of the whole extinction spectrum. In order to overcome this difficulty we added a step to the procedure which consists in determining the wave number limit of the validity of the Rayleigh approximation. In the method described in Herbin et al. (2017) the first step is to obtain an a priori set of indices using the Rayleigh approximation and a fit from the Mie theory. Here we jointly perform a fit from Mie and CDE calculations. The limit wavenumber corresponds to the one for which we obtain an identical value of the imaginary part of the index from both calculations. For this study this value is 1333 cm⁻¹. In addition, compared to previous work Deguine et al. 2020; Herbin et al. 2017; Hubert et al. 2017), we have included in the optimal estimation routine the fact that the fit of the CRIs by the CDE between 650 cm⁻¹ and 1333 cm⁻¹ is constrained by the values obtained from Mie between 1333 cm⁻¹ and 40,000 cm⁻¹. Note that the use of the whole spectrum and the KK relation ensures the homogeneity and continuity of the retrieved CRIs.

The latter are shown in Figure 4 and given as supplementary material.

Figure 4. Real (left panel) and imaginary (right panel) parts of the retrieved complex refractive index of quartz (in black lines). The uncertainties are illustrated by a grey band delimited by black dash.

The characterization of the retrieved quantities in terms of error sources is essential to estimate the quality of the results. A major advantage of the OEM is to provide an efficient way for characterizing the retrieved state, which is given by (Rodgers, 2000): The total error variance-covariance matrix can then be written as the sum of three individual contributions: (9) $${\mathbf{S}}_{\text{total}}={\mathbf{S}}_{\text{smoothing}}+{\mathbf{S}}_{\text{meas}}+{\mathbf{S}}_{\text{modparam}}$$(9)

In Equation (9)(9) $${\mathbf{S}}_{\text{total}}={\mathbf{S}}_{\text{smoothing}}+{\mathbf{S}}_{\text{meas}}+{\mathbf{S}}_{\text{modparam}}$$(9) , the first term ${\mathbf{S}}_{\text{smoothing}}=(\mathbf{A}-\mathbf{I}){\mathbf{S}}_a{(\mathbf{A}-\mathbf{I})}^T$ is the smoothing error, which accounts for the sensitivity of the measurements to the retrieved parameters; the second term ${\mathbf{S}}_{\text{meas}.}=\mathbf{G}\mathrm{ }{\mathbf{S}}_{\epsilon }\mathrm{ }{\mathbf{G}}^T$ is the measurement error (see Section 3.1); and the third one ${\mathbf{S}}_{\text{mod}.\text{param}.}=\mathbf{G}{\mathbf{K}}_b\mathrm{ }{\mathbf{S}}_b\mathrm{ }{(\mathbf{G}{\mathbf{K}}_b)}^T$ is the model parameter error, which represents the imperfect knowledge of the size distribution parameters. These covariance matrices are calculated from the identity matrix I; the gain matrix G, whose rows are the derivatives of the retrieved state with respect to the spectral points and it is defined by: $\mathbf{G}=\partial \stackrel{\mathbf{⁁}}{\mathbf{x}}/\partial \mathbf{y};$ and the averaging kernel matrix A, gives a measure of the sensitivity of the retrieved state to the true state, which is defined by: $\mathbf{A}=\partial \stackrel{\mathbf{⁁}}{\mathbf{x}}/\partial \mathbf{x}=\mathbf{GK}.$ S_b is the variability on the non-retrieved parameters of the model, corresponding here to the aerosol concentration and size distribution. S_b values are given by the uncertainties estimated by the SD retrieval procedure, i.e., 1% for the concentration, 0.5% for D_ve and 1% for σ. The trace of S_total matrix give the total uncertainty at each wavenumber, represented by dash lines in Figure 4, and they are on average for the whole spectral range of 2.5% for the real part n, and 1.7% for the imaginary part κ of the refractive indices.

The calculated extinction spectrum using the retrieved CRI is compared to the experimental one in Figure 5. We can notice the very good agreement quantified by an average root mean square (RMS) of 1.47 10⁻³ calculated for the whole spectral range and 3.68 10⁻⁵ in the infrared between 650 and 1333 cm⁻¹. These RMS are much lower than the measurement noise (see Section 2.3).

Figure 5. Comparison between the experimental quartz spectrum (grey line) and the calculated one at the end of the retrieval process (black line). The two strong rovibrational absorption bands visible around 1600 and 3600 cm⁻¹ on the grey line are due to residual water vapor and they are therefore not considered in the calculated extinction spectrum (black line). The inset shows a zoomed portion of the extinction spectra in the TIR (between 650 and 1350 cm⁻¹).

Optical properties are dependent of the chemical composition and of the sample fabrication process. However, to ensure that the complex refractive indices are those of the material and not of the medium, an additional work has been done. Using the same experimental conditions, the extinction spectrum of another quartz sample (S2) purchased from Sigma-Aldrich (Fluka, purity 99%, particle size < 10 μm, density of 2.6 g.cm⁻³) have been recorded as well as the corresponding size distribution (D_ve = 0.2 µm and σ = 2.4). Figure 6 shows the comparison between the experimental extinction spectrum obtained from S2 and the simulated one calculated by the CDE relation (6)(6) $${\stackrel{\mathbf{⁁}}{\mathbf{x}}}_{i+1}={\mathbf{x}}_a+{({\mathbf{K}}_i^T{\mathbf{S}}_{\epsilon }^{-1}{\mathbf{K}}_i+{\mathbf{S}}_a^{-1})}^{-1}{\mathbf{K}}_i^T{\mathbf{S}}_{\epsilon }^{-1}\left[\mathbf{y}-\mathbf{F}({\stackrel{\mathbf{⁁}}{\mathbf{x}}}_i)+{\mathbf{K}}_i({\stackrel{\mathbf{⁁}}{\mathbf{x}}}_i-{\mathbf{x}}_a)\right],$$(6) using the complex refractive index previously retrieved from S1 sample (Figure 4). Comparison has been done on the whole spectral range, but only an expanded view in the IR spectral region has been selected, where the CRIs values are the most discussed in the literature. This extinction spectrum covers the spectral region 650–1333 cm⁻¹ at a spectral resolution of 1 cm⁻¹. Simulated extinction spectrum from retrieved optical constants presented in Figure 4 shows an excellent agreement with the data in terms of band shapes, peak positions and intensities. The overlay is not perfect, but the comparison is very satisfactory considering the fact that the quartz sample is not exactly similar. This result illustrates not only the consistency of the simulated extinction from the retrieved CRIs with the measurement, but also the good reproducibility of the whole process, both experimental and numerical.

Figure 6. Comparison between the experimental quartz spectrum (grey line) corresponding to sample S2 and the calculated one (black line) using the CRIs retrieved for S1 (Figure 4) and the size distribution measurement of S2.

4. Conclusion

Previous studies demonstrated that Mie theory was suitable for reproducing the extinction of aerosols even if they were not spherical. However, this assumption failed to reproduce the infrared part of the spectrum in particular around the resonance peak around 10 µm for birefringent compounds such as crystalline quartz and calcite. Thus, whatever the scattering theory (Mie, CDE, T-Matrix, DDA) and the refraction indices of the literature used, it appears difficult to correctly simulate the spectrum of crystalline quartz in the thermal infrared. The latter being one of the major elements of mineral dust, this represents a significant source of uncertainty for their retrieval by remote sensing measurements and therefore for the quantification of their radiative impact. This highlights the necessity to improve the quartz optical properties. For this purpose, we used an experimental set-up allows recording simultaneously high spectral resolution extinction spectra and size distribution of homogeneous well-mixed aerosol flow. These data are then combined in a numerical procedure to obtain new quartz optical constants. As the retrieval process of the CRIs from the Mie theory does not lead to a satisfactory convergence in the TIR, we adapted our algorithm in order to replace it by the CDE formula between 650 and 1333 cm⁻¹. Since the latter is limited to small size parameters, the rest of the spectrum is adjusted simultaneously from Mie theory. We then obtain a continuous and coherent set of real n and imaginary $\kappa$ parts of the complex refractive index from 650 to 40,000 cm⁻¹. These new optical constants were then used to fit an extinction spectrum from another different sample of crystalline quartz. The general agreement shows that our approach is suitable to reproduce the experimental quartz extinction spectrum and validates the consistency of the new CRIs. The latter are available on request or in online supplementary materials and they can be introduced in the existing optical constants databases such as ARIA (Aerosol Refractive Index Archive). Finally, recent work exploiting these data has shown the interest of having homogeneous optical properties with high spectral resolution for a better physico-chemical characterization of aerosols from the IASI satellite instrument (Alalam et al. 2022). Current experimental efforts focus on extending the spectral range covered toward the far infrared (100–650 cm⁻¹), particularly in preparation for the FORUM space mission.

Acknowledgments

The authors thank Sylvie Gosselin and Valentine Bizet for the technical support.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work is a contribution to the LabEx CaPPA (Chemical and Physical Properties of the Atmosphere) project funded by the French National Research Agency (ANR) under contract “ANR-11-LABX-0005-01” and to the CPER research project CLIMIBIO funded by the French Ministère de l’Enseignement Supérieur et de la Recherche. The authors thank the Regional Council ‘Hauts-de-France’ and the European Regional Development Fund for their financial support for these projects.