A Comprehensive Review on Study Methods of Aerosol Optical Properties in Different Dimensions

A comprehensive understanding of the optical properties of atmospheric aerosol is essential for a variety of applications, such as optical imaging, optical communication, and remote sensing. In recent years, many theories and numerical simulation methods have been developed to connect aerosol physical-chemical properties to their intrinsic and integral optical properties. Usually, simulations and measurements are intertwined to synergistically attain the retrieval of aerosol optical properties and mitigate or even eliminate the adverse impacts of aerosol during imaging, sensing, or communication. This review covers the fundamental theories of aerosol optical properties, the development of numerical simulations, the instrument-based sampling measurements, and the cutting-edge techniques of remote sensing. Numerical simulations have been progressing from symmetric particles to asymmetric particles over the past two decades, although any simulation method is limited by specific shape and a restricted size parameter range. Thus, this review also examines the most typical advances in aerosol instrumentation that are frequently used to measure the intrinsic optical properties of unknown aerosols. Such obtained properties validate simulations and constitute the basis of integral optical properties. In terms of practical applications, integral optical properties are the most critical knowledge about atmospheric aerosol. Remote sensing measurements, be it ground-based, airborne, or satellite-based, all retrieve integral optical properties of atmospheric aerosol from various perspectives, which are elaborated upon in this review. In conclusion, this review provides an all-encompassing comprehension of aerosol optical properties.


I. INTRODUCTION
Atmospheric aerosols are generally defined as solid or liquid airborne particles with size distributions ranging from 0.001 µm to 100 µm. These aerosols can be classified into various types, based on their wide size range spanning four orders of magnitude and complex formation processes. Consequently, they are divided into primary aerosols, secondary aerosols, and aerosol precursors according to their formation processes, and urban aerosols, desert aerosols, volcanic aerosols, and stratospheric aerosols according to the environments they inhabit [1]. Despite the diversity of atmospheric The associate editor coordinating the review of this manuscript and approving it for publication was Jon Atli Benediktsson . aerosols, their sizes typically follow a log-normal distribution, particularly in the size range spanning several orders of magnitude [2]. The aerosol optical properties (AOPs) are some of the most important physical properties, which are determined by macroscopic statistical parameters such as number concentration, mass concentration, and surface area concentration, as well as microscopic intrinsic parameters, such as spectral refractive index, dimension, shape, and composition. At present, aerosol optical thickness is the most widely used macroscopic parameter. Moreover, the AOPs, including absorption spectrum, optical scattering, and polarization degradation of aerosols, have a significant impact on optical remote sensing, optical imaging, optical measurement, and optical communication.
The weak absorption and strong scattering of solar radiation by aerosols can significantly influence remote sensing measurements by altering the planetary albedo of the Earth's atmospheric system [3]. Due to the similarity between the particle size of aerosols (fog, smoke, dust) and the optical wavelengths in free-space optical (FSO) communications, visible or infrared beams are strongly scattered, resulting in an attenuation of the optical signal [4]. In the field of imaging through aerosols, the absorption and scattering of light by aerosols can reduce the signal light from a target and increase the background light, leading to a blurred image [5]. If the AOPs have been obtained, the detrimental effects of aerosols on optical applications can be suppressed or completely avoided. For instance, selecting an optical channel with the lowest absorption according to the absorption spectrum of the aerosols can enhance the reliability of the communication link of the FSO system [4]. Accurately measuring the optical thickness of the atmospheric aerosols can enable the compensation of light intensity in optical remote sensing through digital operation [6]. In the field of optical imaging, blurred images can be reconstructed by numerically integrating the optical scattering along the optical link [7], [8]. Furthermore, measuring and analyzing the relationship between wavelength and aerosol light absorption can also help to evaluate the apportionment of aerosol sources, quantify the relative contribution of different types of light-absorbing aerosols, thus simulating aerosol radiative forcing (RF), and providing improved parameters for remote sensing [9]. This paper does not aim to provide an exhaustive overview of AOPs, but rather focuses on the characteristics that have a significant impact on optical remote sensing, imaging, and optical communication. Previous works have summarized the AOPs obtained from the development of experimental methods and instruments [10], [11], [12]. With the advancement of computer technology, numerical simulation has become a popular approach to studying the AOPs due to its adjustable parameters and universality. In this review, we first introduce a theoretical model for the analysis of AOPs and some physical parameters affecting the AOPs in detail. Second, we describe the principles and processes of numerical simulation for the AOPs developed so far, and compare the strengths and weaknesses of various simulation methods and working conditions. Third, we present the measurements of AOPs, including sampling measurements and remote sensing measurements. As illustrated in Fig. 1, this review covers three types of research methods for the AOPs. Numerical simulations consider both the microscopic properties of individual aerosol particles and the macroscopic statistical properties of atmospheric aerosols, thus providing a link between the aerosol intrinsic parameters and the AOPs. Sampling measurements focus on the intrinsic properties of aerosols that are independent of the environment in which they are located, while remote sensing measurements focus on the macroscopic AOPs in specific environments. Through the analysis of different dimensions, this review serves as a guide for future research on AOPs.

II. THEORETICAL OPTICAL CHARACTERISTICS OF AEROSOLS
A significant factor influencing the process of atmospheric radiative transfer is the light scattering of aerosols. Variations in the scattering properties of aerosols can alter the planetary albedo of the Earth's atmospheric system, thus impacting the energy balance of the Earth's atmospheric system [13]. Examining the Intergovernmental Panel on Climate Change (IPCC) reports over the years [14], it can be observed that, despite advances in the assessment accuracy of aerosol direct radiative forcing, there remain considerable uncertainties (−0.9 W·m3 -−0.1 W·m3). One of the key elements that restrict the accuracy of the assessment is the uncertainty of aerosol scattering properties, particularly the scattering properties of non-spherical particles such as ice crystals and cirrus clouds.
Despite the efforts of many researchers in constructing databases to describe the light scattering properties of aerosols [15], the increasing diversity of aerosol types in radiative transfer models has yet to be addressed. Owing to the complexity of microphysical parameters such as aerosol shapes, chemical components and scale spectra, these databases are unable to encompass the scattering properties of all aerosols, thus inevitably resulting in errors in the simulation of radiative transfer processes [16], [17]. Consequently, a comprehensive understanding of the light scattering properties of aerosol particles and their impacts on radiative transfer simulations is of paramount importance at present [18].
In 1908, Mie and Lorenz independently proposed a scattering model of homogeneous spherical particles, which enabled them to accurately calculate the light scattering properties of the particles by solving Maxwell's equations under spherical boundary conditions. As a result, the scattering from homogeneous spheres is also known as Mie scattering or Lorenz-Mie scattering [19], [20]. Nevertheless, experiments have revealed that the Mie scattering theory is not suitable for simulating the scattering properties of non-spherical particles, as it tends to underestimate the transverse scattering effects compared to other scattering theories. For instance, Hoyningen-Huene et al. found that the Mie scattering model can lead to an error of up to 60% in the scattering of large aerosols with size parameters exceeding 80 [21]. Similarly, Koepke et al. reported that the discrepancies between spherical and non-spherical particles can reach up to 60% in the transverse and backward scattering regions of the phase function in the solar spectral range [22]. Moreover, Sorribas highlighted that if spherical particles were used to assess the AOPs, the correction must be reconsidered to reduce the uncertainty of scattering in rough mode [23].
In this section, the optical properties of non-spherical scatterers will be discussed first, followed by a comparison of these definitions with those of the generally simpler spherical scatterers, thus providing a theoretical foundation for the subsequent introduction to simulations and measurements of light scattering.

A. THEORY OF ELECTROMAGNETIC SCATTERING
The electromagnetic radiation by any particle can be characterized by its electric vector ⃗ E and magnetic vector ⃗ H . It is generally assumed that plane waves can accurately represent radiation from very far-away point sources.
where k = 2π/λ is the wave number, λ is the wavelength, the vector ⃗ n in represents the direction of wave propagation, and ⃗ r is the radius vector. The right-angle coordinate system to be used subsequently is the same as the origin of the spherical coordinate system (r, θ, ϕ) adopted here, with ⃗ i θ ⃗ i ϕ as the unit vector. E in θ and E in ϕ denote the radiation of the incident field in the θ and ϕ directions, respectively.
As Fig. 2 illustrates, both reflective and refractive fields can be calculated using Fresnel's Formula. Furthermore, we only consider time-dependent harmonic fields, since any field can be decomposed into a Fourier series, and the components of the series can be treated independently.
l ) inside the particle and scattered field ( ⃗ E s , ⃗ H s ) in the medium around the particle.
For time-harmonic fields, the wave equation can be derived from Maxwell's electromagnetic equations and material equations. In the case of a spherical symmetric scatterer, the variables in the equation and boundary conditions can be separated to obtain the exact solution, i.e., Mie theory. However, for non-spherical scatterers, such separation is not possible, making the solution of the scattering problem highly challenging [24]. In this context, the appropriate selection of the axis orientation of the Cartesian coordinate system is of particular importance.

B. AMPLITUDE SCATTERING MATRIX AND MUELLER MATRIX
If only far-field scattering (kr ≫ 1) is considered, the scattered field can be approximately treated as a transverse field, which can be depicted in the form of spherical waves: where E sca θ,0 and E sca ϕ,0 denote the propagation of scatter radiation along the direction of the z-axis and the direction of the x-axis, respectively. For any single scatterer, the linearity of Maxwell's equations guarantees that the linear invariance of incident, internal and scattered fields. Assuming that the scattered field is represented by a spherical wave, the amplitude scattering matrixŜ can be used to express the linear relationship between the incident field and scattered field.
In the particle-connected coordinate system, selection of TE or TM mode (utilizing S 2 , S 4 or S 3 , S 1 , respectively) yields two solutions, ⃗ E sca TE and ⃗ E sca TM , for the scattered radiation.
The degree of polarization and amplitude scattering matrixŜ are both determined by the direction of scattering, which can be described by the Stokes vector in terms of the relationship between the scattered and incident fields.
where (I sca , Q sca , U sca , V sca ) T and I in , Q in , U in , V in T are the Stokes vectors of scattered and incident radiation, respectively. I is equal to the total intensity, Q and U fully describe the linear polarization of light, and V denotes the left-or right-circularly polarized light component. The 4 × 4 matrix is a scattering matrixẐ , which is also known as the Mueller matrix of individual particles. The Mueller matrix elements can be accurately represented by the amplitude scattering matrix elements. For spherical scatterers, the Mueller matrix has 10 elements equal to zero; however, the Mueller matrix is always complex for non-spherical scatters. VOLUME 11, 2023 36765 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.

C. EXTINCTION, SCATTERING, AND ABSORPTION PROPERTIES
The electromagnetic radiation theory under far-field conditions facilitates the calculation of certain optical properties of scattering objects. For instance, the extinction cross section C ext indicates the proportion of incident light flux that is absorbed and scattered by particles. When C ext is multiplied by the irradiance of the electromagnetic wave incident on an object, the total amount of radiation scattered and absorbed by the object is obtained. The C ext can be expressed as: where θ is the scattering angle, i.e., the angle between the incident and scattered directions, and ⃗ i is the unit vector denotes the imaginary part of the projection of the inner product of ⃗ E sca and ⃗ i in the positive z-axis direction.
The scattering cross section C sca is determined by the direction and polarization state of the incident light,which can be defined as follows: The absorption cross section C abs can be defined by the law of conservation of energy as follows: The extinction cross section for a polarized incident light can be expressed as: (8) where φ denotes the angle of polarization. I in , Q in , U in are the elements constituting the Stokes vector of incident radiation. The single scattering albedo (SSA) ϖ , widely utilized in radiative transfer theory, is defined as the ratio of the scattering cross section to the extinction cross section ratio, i.e., ϖ = C sca C ext ≤ 1. The use of the single-scattering albedo (SSA) can be utilized to explain the likelihood that photons will be scattered rather than absorbed by particles. The calculation of the SSA for a polarized incident light can be accomplished as follows: The asymmetry parameter g (or ⟨cos θ⟩) is defined as the average cosine of the scattering angle. The phase function of a spherical particle exhibits rotational symmetry with respect to the propagation direction of the incident radiation; however, for a non-spherical particle, such symmetry is usually absent. Thus, the asymmetry parameter is a scalar for spherical particles, but a vector for non-spherical particles. In addition, the asymmetry factor is positive if the scattering is forward (θ = 0), and negative for backward scattering (θ = π). Furthermore, the value of the asymmetry factor is zero when the scattering is symmetric with respect to a plane perpendicular to the incident direction. The asymmetry parameter can be calculated as follows: where cos θ is scattering angle's cosine and p is the symbol of phase function. A comprehensive explanation of the procedure for computing the asymmetric parameter can be found in the reference [25].

D. OPTICAL THICKNESS AND ANGSTROM EXPONENT
The extinction, scattering, and absorption efficiencies (or efficiency factors) of the particles are quantified as: where G is the cross-sectional area projected onto a particle by a plane perpendicular to the incident beam (e.g., the crosssectional area G = πa 2 for a sphere with radius a). Utilizing the intuitive concept of geometric optics, it can be determined that the extinction efficiency of ordinary particles cannot exceed 1. However, there are numerous particles that are able to absorb and scatter more light than the incident light; thus, the efficiency can be viewed as a dimensionless cross section [26]. The optical properties of aerosol particle populations are closely related to their particle size distributions n(r), thus necessitating the integration of the scattering and absorption cross sections (in m −1 ) over all particles in order to calculate the scattering and absorption coefficients.
The extinction, scattering, and absorption efficiencies of aerosols are atmospheric local properties that vary with position and altitude. The integral of the extinction coefficient along the vertical direction is referred to as the aerosol optical thickness (AOT) or aerosol optical depth (AOD) and can be calculated as follows: The alteration of wavelength may have a significant effect on the AOPs; thus, it is essential to link the wavelength to a particular AOD. The Angstrom exponent (also referred to as the Angstrom parameter or Angstrom index) is utilized to depict the spectral dependence (or extinction coefficient) of the aerosol associated with the AOD. It is calculated using the following equation: 36766 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
where τ 1 and τ 2 denote the AODs at wavelength λ 1 and λ 2 , respectively, the Angstrom exponent varies with wavelength, but in the visible range, the variation is usually smaller than in infrared range. Practically, the AOD at wavelength λ can be estimated with the reference value (the reference AOD τ r at the reference wavelength λ r and the Angstrom exponent α) using the following equation: There is an unclear relationship between the Angstrom Exponent (AE) and aerosol size due to an array of factors, such as the aerosol refractive index, influencing the spectral dependence of light extinction. However, light absorption is independent with aerosol size and refractive index, but depends on aerosol types. Nevertheless, a more precise estimation of the fine mode aerosol proportion can be attained by taking into account the spectral curvature [27]. Futhermore, the absorption Angstrom exponent β is specified as follows: where τ abs and τ abs r mean the aerosol absorption optical depths (AAOD) at wavelengths λ andλ r , respectively. The Angstrom exponent for absorption is dependent on the aerosol type and can be leveraged for the classification of aerosol types [28].

III. NUMERICAL SIMULATION OF AEROSOL LIGHT SCATTERING
Electromagnetic scattering simulation is a confluence of applied physics and computational mathematics. The scattering of light from uniform spherical particles is one of the most general problems in this area. Several light scattering theories have been developed, each with their respective advantages and disadvantages. These theories have been described in detail by Wriedt et al., Kahnert et al.,and others [29], [30], [31]. For example, Wriedt et al. explored the practical applications of different algorithms in aerosol science and examined the advances in light scattering theory and computational procedures up to 2009 [31]. Similarly, Kahnert et al. explored the theoretical basis, pros and cons of various numerical calculations based on a rigorous electromagnetic theory, and the fundamentals of using particle models of various complexities in light scattering calculations [28], [32]. The theories (and applications) of light scattering have been continually developed in two directions: (i) from plane waves to beams of arbitrary shapes, and (ii) from uniform spherical particles to particles of arbitrary shapes. Gouesbet et al. offered a review to provide an integrated and unified understanding of these two trends [31].
Numerical methods for calculating light scattering of nonspherical particles can be divided into two categories: exact solution models and approximate calculation models. Vector wave function expansion, volume integral equations, and differential equations are the exact solution models, while the approximate computational models are not further subdivided. Table 1 outlines the applicability, particle types, algorithm complexity, and the advantages and disadvantages of the various models.

A. METHODS BASED ON VECTOR WAVE FUNCTION EXPANSION
Vector-based wave function expansion involves the expansion of the scattered, internal, and incident fields by vector wave functions such as spherical wave, cylindrical wave, and ellipsoidal wave. Subsequently, a differential or integral equation is solved, taking into account specific boundary conditions, in order to determine the relationship between the incident and scattered field, thereby enabling numerical calculation of the scattering. Commonly used models in this field include the T-matrix method (TMM), the point matching method (PMM), the Nystrom Method (NM), and the separation of variables method (SVM) [32].

1) T MATRIX METHOD (TMM)
Since Waterman first proposed the T-Matrix Method (TMM) in 1965 [33], this method has been continuously improved, developed and widely used. Currently, TMM is considered one of the most effective and extensively used approaches for numerical simulation of light scattering of resonant nonspherical particles. The elements of the T-matrix are independent of the incident and scattered fields and only depend on the shape, size parameter and refractive index of the scattered particles, as well as the direction relative to the reference frame. As a consequence, the T-matrix only has to be computed once and can be used to calculate incident and scattered light in any direction [34]. The method of calculating the Tmatrix under a single uniform scattering condition is referred to as the extended boundary condition method (EBCM). Mishchenko et al. provided a comprehensive overview of the most extensive peer-reviewed T-matrix publishing theme database from 2004 to 2019 [35].
Numerical stability (or convergence) has long been a key concern in TMM calculations. Ding et al. investigated the numerical iterative convergence of TMM and found that Mishchenko's ''mathematical convergence method'' sometimes yielded unreasonable results [36], [37], with the extinction efficiency factor Q ext being smaller than the scattering efficiency factor Q sca . Dallas proposed a condition operator for examining the convergence of TMM based on numerical analysis [38]. Kyurkchan et al. demonstrated that TMM could provide correct and stable results if the surfaces covering the scatterer's interior satisfied the null field criterion [39]. Zhai et al., building on the Invariant Imbedding T-matrix Method (II-TM), utilized Gauss Legendre (GL) quadrature to generate nodes and weights, which significantly enhanced convergence and computational efficiency [40]. The computational region of the II-TM-based non-spherical scatterer is illustrated in Fig. 3. The spherical core region is a homogeneous sphere whose T matrix can be calculated by VOLUME 11, 2023  Lorentz-Mie's theory. In the spherical mantle region, the Invariant Embedding technique is applied to solve the T matrix by setting its initial iteration value of the T matrix of a spherical core [41]. Yang et al. discussed the convergence of various methods (including II-TM) for solving electromagnetic scattering of non-spherical particles [42]. Kahnert et al. derived sufficient conditions for the convergence of the iterative T-matrix algorithm by utilizing the Banach Fixed Point Theorem [43].
TMM is used for calculating light scattering from a variety of uniform ellipsoids, including smooth ellipsoids, Chebyshev particles, generalized Chebyshev particles, and hyper ellipsoids, as well as finite cylindrical particles with sharp edges, polyhedral, eccentric inclusions, lamellar particles, and collections of multiple spheres. Furthermore, it can also account for light scattering from nonspherical particles such as aerosol particles and ice crystal particles.

2) POINT MATCH METHOD (PMM)
Unlike SVM and TMM, the point matching method expands the electromagnetic fields using truncated vector spherical wave functions. Through the selection of a limited number of matching points on the particle's surface, it is possible to determine the expansion coefficients of the electromagnetic fields, with the constraint of boundary conditions applied. Subsequently, the corresponding algebraic equations can be solved, with the number of matching points equating to the number of unknowns in the equation.
The accuracy of the results obtained from the PMM is largely contingent upon the selection of matching points. Consequently, inaccurate results or unstable numerical computations may occur when using the PMM [44]. To address this issue, Oguchi et al. proposed a generalized point match method (GPMM), which doubles the number of matching points as compared to the number of unknowns, thus causing the linear equation solving process to be overdetermined. The expansion coefficients can then be solved by least squares [45]. In comparison to PMM, GPMM is numerically more stable as the matching points of GPMM fulfill the boundary conditions [46]. Furthermore, Kahnert provided an in-depth description of the particular differences between the PPM and GPPM. [47]. Ohnuki et al. applied a region decomposition technique, which used PMM to solve the electromagnetic scattering from various polygonal conductor columns, and investigated the influence of scatterer shape variation in the shadow region on electromagnetic scattering [48]. Additionally, Nieminen, et al. applied the PMM to calculate the T-matrix for asymmetric particles, thereby avoiding the large time consumption in computing the surface integral when using EBCM [49]. Farafonov et al. replaced the relevant integrals in the Particle Mesh Method (PMM) with a summation of points on the particle surface, thereby increasing the accuracy of the calculation by several orders of magnitude in comparison to the commonly used Extinction-Based Cluster Method (EBCM) for particles with a ratio of maximum size to minimum size of less than 1.5 [50]. Subsequently, Loke et al. proposed a technique combining Diffraction-Discrete-Approximation (DDA) and Point-Particle-Method (PPM) to develop a T-matrix for simulating the light scattering of micro-sized particles of arbitrary shapes, utilizing the discrete rotation symmetry and mirror symmetry of the micro-sized particles in order to reduce memory usage and computation time by several orders of magnitude [51]. Yannopapas employed the Point Matching Module to amalgamate the DDA and Layer-Multiple-Scattering (LMS) methods in order to acquire the T-matrices of arbitrary scatterers for investigating the optical properties of 2D and 3D photonic lattices [52]. Ohnuki et al. furthermore extended the PMM to solve 3D electromagnetic scattering problems [53].

3) THE NYSTROM METHOD (NM)
The NM and PMM differ in terms of their selection of matching points: the former uses quadrature rules to select points, which involves a summation of discrete area elements rather than a continuous integration over small domains, while the latter selects points arbitrarily, but the matrix must be evaluated by numerical integration. Compared to the TMM and PMM, the NM has certain operational advantages: discretization of integral equations and removal of primary and test functions are straightforward, the requirements for grid quality are low, and higher-order integration can be used to improve numerical accuracy. However, it necessitates the accurate and effective treatment of hyper-singularity [54]. In recent years, various improvements to the NM have been put forward by researchers. For example, Tong et al. applied the NM to the solution of electromagnetic scattering from three-dimensional penetrable objects [55]. Yang et al. developed an efficient NM for solving electromagnetic scattering problems with non-uniformly anisotropic media [56]. Chen et al. developed a new high-precision NM with few matching points, which was verified through numerical calculations on multiple scatterers [57], and Zhang et al. derived effective analytic formulas to facilitate the practical application of the NM [58]. Furthermore, taking into account the similarities and differences between the NM and PMM, numerous hybrid strategies have been proposed, such as the one by He et al., which combines the NM and PMM to solve the electromagnetic scattering problems of open structures with no thickness, and which was verified by simulation calculations [59].

4) SEPARATION OF VARIABLES METHOD (SVM)
Oguchi and Asano et al. initially proposed the SVM to calculate the light scattering from uniform prolate spheroidal particles or oblate ellipsoid particles of arbitrary size using the separation of vector wave equations in the ellipsoidal coordinate system [60], [61]. Subsequently, Voshchinnikov et al. and Farafonov et al. further improved the SVM in order to attain a higher computational efficiency, especially for extremely prolate spheroidal particles or oblate ellipsoid particles [62]. Vinokurov, A. subsequently applied the SVM to calculate the scattering and absorption cross sections of two-layer confocal ellipsoidal particles and compared the results to those obtained from other algorithms, such as EBCM, finding that SVM was more suitable for confocal multilayer particles [63]. Additionally, Farafonov et al. used SVM to calculate non-spherical scatterers of up to 100 layers and compared it to the EBM (or EBCM), finding that SVM exhibited better performance [64].
The truncation error of the Support Vector Machine (SVM) is derived from the limitation of the vector ellipsoidal wave function on the expansion order of the electromagnetic field. When the particle radius or refractive index is large, the computation procedure of the associated linear algebraic equation system may also be subject to errors [47]. With regard to the SVM and the Extinction Boundary Condition Method (EBCM) algorithms, disparate boundary conditions result in disparate linear equations. The EBCM contains two systems of equations, whereas the SVM produces only one system of equations with a number greater than the square of its dimension, and thus, the EBCM is generally favored for practical applications.
Farafonov et al. conducted an analytical and numerical analysis of the applicability of TMM, PMM, and SVM, three vector wave function expansion-based scattering simulation models, in both near-field and far-field regions, while taking into consideration the Rayleigh assumption [65]. The results of this investigation suggested that these methods are interrelated and equivalent.

B. METHODS BASED ON VOLUME INTEGRAL EQUATIONS
The proposed method based on the volume integral equation is employed to realize the discretization of the integral equation. This method utilizes the Green's function of the surrounding medium, thereby automatically satisfying the radiation condition of the scattered field, and resulting in a full matrix for solution. Popular models for this purpose include the Method of Moments (MoM), the Discrete Dipole Approximation (DDA), and the Fredholm Integral Equation Method (FIEM).

1) MOM AND DDA
The Method of Moments (MoM) transforms volume integrals into finite sums of elements through discrete methods, resulting in the simplification of the volume integral equation into a linear system of equations, which can be solved using the conjugate gradient method or Gaussian elimination. Gibson provides an in-depth description of how to obtain numerical solutions to electromagnetic integral equations using MoM [66]. Furthermore, certain literature reviews have evaluated the relationship between MoM and the Discrete Dipole Approximation (DDA) [67], as both methods derive a linear system of equations that can be solved using standard numerical algorithms. DDA has been further improved on the basis of MoM, and the primary distinction lies in the calculation of self-induced electric fields.
The Discrete Dipole Approximation (DDA) method, first proposed by Purcell and Pennypacker [68] -also known as the Coupled Dipole Method (CDM)-was improved by Draine to reduce its computational complexity by utilizing the Conjugate Gradient Method together with the Fast Fourier Transform (CGM-FFT) method to solve the linear matrix equations [69]. A comprehensive review of the theoretical and numerical development of DDA up until 2007 can be found in literature [70]. Yurkin and Hoekstra investigated the open source code ADDA for spheres with a large size parameter (χ = 320) [71]. Huang et al. combined the electronic structure of materials to transform aerosol condensed particles into a series of discrete dipoles [72], and then used DDA to numerically calculate the optical characteristic values of scattering, absorption, and extinction efficiency factors at different wavelengths. Vartia et al. quantitatively studied the simulation accuracy of DDA for particles of varying shapes [73]. Shabaninezhad et al. employed a graphics processing unit (GPU) in MATLAB to accelerate the calculation process, as illustrated in Fig. 4, and realize the calculation of optical properties of arbitrary shaped plasmonic nanoparticles (NPs). According to the results, a GUP can perform calculations approximately 10 times faster than a CPU for 1 million dipoles [74]. Due to the efforts of Draine, the Discrete Dipole Approximation (DDA) has now become a widely used scattering algorithm. However, the algorithm is characterized by low numerical accuracy and a computationally complex process. What's more, properly simulated media are limited to particles with a large size parameter and a small complex refractive index. With the advancement of computer technology, the optical properties of particles with large size parameters can be obtained through parallel computing. Yurkin et al. conducted a thorough convergence analysis of the DDA algorithm and implemented general simulations of particles with extremely high refractive indices using the ADDA code [75], [76]. Furthermore, Flatau and Draine proposed a fast DDA algorithm for gridded data, which requires less computation time for near-field light scattering calculations of large particles [77]. Groth et al. accelerated the computation of DDA by using a preprocessing method that made use of a circulant matrix to approximate the system matrix, allowing for an effective inversion of the FFT. This accelerated DDA was verified by simulation to reduce computation time by several orders of magnitude in most cases [78]. DDA has been widely applied to scatterers of different shapes, such as triaxial ellipsoids [79], hexagonal prisms [80] and decahedra [81].

2) FREDHOLM INTEGRAL EQUATION METHOD (FIEM)
The Finite Integration Element Method (FIEM) was first proposed by Holt et al. [82]. Through the Fourier transform of the internal field, this method transforms a volume integral into an integral over the wave number coordinate system, thereby circumventing the singular kernel problem encountered in the integral equation method. However, the FIEM algorithm necessitates that the elements of the solved matrix be resolvable, making it only applicable to particles with shapes of ellipsoids and finite cylinders. Moreover, the corresponding algorithms vary greatly according to the shape of the particle. Subsequent to its initial proposal, Papadakis et al. extended the applicability of the FIEM to anisotropic dielectric ellipsoidal particles and three-dimensional ellipsoid particle chains [83], while Ngobigha et al. improved the method for scattering calculations of irregular inhomogeneous dielectric particles [84].
One of the great advantages of FIEM over other models is its integration process, which is only related to the particle size parameter and shape and thus does not require repeated calculation when the incident light state changes, surpassing the abilities of the DDA and other models based on differential ideas (discussed further in the next section). However, the disadvantage of FIEM is that it is only applicable to smallscale particles with a small complex refractive index, and the model is computationally intensive [85]. To improve its computational efficiency, many researchers have proposed and tested methods. Alipanah and Dehghan suggested a method based on radial basis functions (RBF) interpolation to approximate the solution of Fredholm nonlinear integral equation, and verified its effectiveness and adaptability through numerical arithmetic examples [86]. Hosseinzadeh et al. studied the stability of the RBF-based method, applied it to the interpolation of scattered data, and mathematically proved its positive definiteness [87]. Zhang and Deng further improved the FIEM by employing the Schaubert-Wilton-Glisson and edge (SWG-Edge) hybrid basis to solve multi-boundary inhomogeneous dielectric objects; the unknowns were only 71% of those used in the general SWG basis and the memory was halved [88].

C. METHODS BASED ON DIFFERENCE IDEA
The difference-based methods discretize a given space using a series of small grids and simulate the electromagnetic field information of each grid by solving the electromagnetic propagation equation. This, in turn, enables the calculation of the external scattered field of a particle. To approximate electromagnetic waves propagating to infinity, a strong absorption boundary is typically used to limit the computational region. These methods mainly include the Finite Difference Time Domain (FDTD), Pseudo Spectral Time Domain (PSTD) and Multi-resolution Time-Domain (MRTD) approaches.

1) FDTD METHOD
The Finite-Difference Time-Domain (FDTD) method, initially proposed and developed by Yee and Chen and summarised by Yang and Liou [89], [90], is the most straightforward approach for solving Maxwell's curl equations. Solving Maxwell's equations in the frequency domain is equivalent to solving a boundary-value problem, which is generally more difficult than solving an initial-value problem [91]. The FDTD method is suitable for calculating the electromagnetic scattering of small, complex-shaped and inhomogeneous particles, directly from solving Maxwell's curl equations via the time domain method. However, the limitation of the FDTD method is that a small grid size is needed to increase the accuracy of the calculation, but this significantly increases the runtime and computational requirements [92].
Zheng et al. proposed an unconditionally stable threedimensional (3-D) FDTD method, which was analytically proven to be unconditionally stable and its effectiveness was verified by numerical experiments. The number of iterations was at least four times fewer than a common FDTD method at the same accuracy [93]. Cole reported a high-precision Yee algorithm based on second-order nonstandard finite differences (NSFDs) and numerically demonstrated its accuracy [94]. For a fixed frequency and grid space, NSFDs increased the maximal time step by 20% in two-dimensional space and 3% in three-dimensional space compared to the Yee FDTD method that employed standard finite differences (SFDs). Ahmed et al. proposed a threedimensional unconditionally-stable locally-one-dimensional finite-difference time-domain (LOD-FDTD) method, which, when compared to the conventional three-dimensional alternating direction implicit FDTD (ADI-FDTD) method, resulted in a runtime reduction of approximately 20%, owing to the decreased arithmetic operations [95]. Eng reviewed several one-dimensional FDTD methods and demonstrated that the simulation could be performed quickly and more efficiently using the fundamental alternating-direction-implicit finite-difference time-domain (FADI-FDTD) method [96]. Although the FDTD method is generally robust, it is susceptible to dispersion errors. In addition to the traditional methods of reducing dispersion, i.e., increasing the sampling rate or using higher order of accuracy, several new methods have been proposed. For example, Shlager and Schneider compared the accuracy of several 2-D low-dispersion FDTD methods, all of which had reduced dispersion errors when compared to the classical Yee FDTD [97]. Finkelstein and Kastner proposed a simplified scheme to weaken the FDTD dispersion errors based on a modification of the characteristic equation [98]. Kong et al. proposed a three-dimensional one-step leapfrog hybrid implicit-explicit finite-FDTD (HIE-FDTD) method, which was demonstrated to reduce the numerical dispersion error without augmenting the computational expense [99].

2) PSTD
Post-Staggered Time Domain (PSTD) is a variation of the Finite Difference Time Domain (FDTD) method. FDTD utilizes the finite difference technique to calculate the spatial derivatives of electric and magnetic fields, whereas PSTD employs the spectral method. More specifically, PSTD utilizes the Fourier transformation and its inverse transformation to represent the spatial derivatives as the integral transformation of the spectral domain [100]. There is no differential approximation for the PSTD method, thus making it more accurate in principle than the FDTD method. During the algorithm design of a PSTD method, two important problems need to be resolved: first, the expansion of the trigonometric basis function of the point circle during the fast Fourier transform is inaccurate (known as the Gibbs phenomenon); and second, the discontinuity of the spatial medium will induce some arithmetic errors of the FFT. Liu et al. and Li et al. were able to successfully address these two problems by combining the PSTD with a random wave superposition model of roughness and applying the model to the simulation of scattering properties of ice crystal particles, allowing a maximum particle scale parameter of 200 [101], [102]. Based on this model, they established a database of scattering properties of ice crystal particles from Rayleigh scattering to geometric optics, which is widely used in various radiative transfer modes. The traditional TF/SF technique is not applicable to PSTD, as it is too time consuming to use the scattering field technique to introduce the incident field. To this end, Hu et al. [103] proposed a weighted TF/SF technique that could accurately introduce incident light by modifying the electromagnetic elements in the connection region between the total field and the scattered field regions. The scattering phase matrix and integral scattering parameters measured by the improved PSTD method were found to be in excellent agreement with the well-tested particle model.
The two main methods used to calculate the light scattering of dielectric particles, PSTD and DDA, can be challenging to differentiate in terms of which is more suitable for particular applications. Liu et al. conducted a numerical simulation of both PSTD and DDA, and determined that the PSTD method was more efficient when the refractive index was greater than 1.4. This finding applied to a wide range of size parameters [104]. Podowitz et al. further extended the comparison between the two methods to the aspect of absorption [105]. Their results indicated that the PSTD method was more applicable when the imaginary part of the refractive index was less than 10 −3 and the size parameter was greater than 40. However, when the particles exhibited strong absorption, the calculation speed of the DDA method was approximately twice as fast as the PSTD method for particles with size parameters up to 100.

3) MRTD
The Modified Resolution Time Domain (MRTD) approach, originally suggested by Krumpholz and Katehi [106], is a variant of the Finite-Difference Time-Domain (FDTD) method. This method entails expanding the electromagnetic field components in the spatial domain utilizing Battle-Lemarie scale functions and wavelet functions, followed by expanding in the time domain utilizing rectangular pulse functions. The Galerkin method from the Method of Moments (MoM) is then utilized to discretize and solve the Maxwell curl equation, allowing for the derivation of the iterative formula of the MRTD method. Although the MRTD method is capable of accelerating the computation speed in comparison to the FDTD method, it is limited in that it can only calculate the scattering properties of particles for a specific size (or orientation) at a given wavelength at one time [107]. This results in an immense amount of calculations if the simulated particles present a wide size distribution. Hu et al. proposed an approach for multi-size simultaneous computation (MSCS) in multi-resolution time domain, whose fundamental principle is depicted in Fig. 5. This approach calculates the transient near electromagnetic field by means of the MRTD method and simultaneously computes the scattering properties of particles of different sizes in a single wave-particle interaction simulation [108]. The uncertainty of a non-spherical aerosol Mueller scattering matrix is an important factor that affects the accuracy of polarization sensing. Therefore, determining a suitable nearfar field transformation approach is a critical step for simulating the MRTD scattering model with greater accuracy. Shuai et al. conducted a systematic comparison of the simulation accuracy of the Mueller scattering matrixes obtained from two near-far field transformation approaches: surface integration based on Huygens' principle and medium volume integration based on the Helmholtz equation, respectively [109]. The results revealed that the near-far field transformation approach based on the volume integration principle had better performance in terms of the error distribution of the Mueller matrix calculation. Additionally, to address the problems of long running time and memory consumption of serial MRTD scattering models, Shuai H. proposed a parallel computational model for non-spherical aerosols based on the message passing interface (MPI) technique. The results indicated that the MRTD model was able to simulate the scattering properties of a non-spherical particle more accurately and the parallel computing technique significantly improved the computation efficiency [107].

4) FEM
First proposed by Morgan and Mei the Finite Element Method (FEM) has become applicable to the simulation of the scattering properties of non-homogeneous particles with arbitrary shape [110]. This method spatially discretizes the Helmholtz equations in order to solve the electromagnetic scattering problem numerically as a side-value problem. Additionally, in order to reduce the number of variables and memory consumption, the FEM requires being performed in a limited area. Similar to the FDTD method, strong absorption boundary conditions need to be applied in order to suppress the electromagnetic wave reflections. To overcome the limitations of absorption boundaries and obtain an approximate solution of electromagnetic wave propagation to infinity, Volakis et al. added a surface integral equation to the FEM method, although this destroyed the sparsity of the coefficient matrix of the linear equations [111]. Ru-Shan et al. applied the symmetric successive overrelaxation (SSOR) preprocessing scheme to the conjugate-gradient (CG) method to solve a large system of linear equations generated when performing the edge-based FEM [112]. This method, when employed in solving large-scale time-harmonic electromagnetic field problems using edge finite elements, requires only one-fifth of the CPU time needed by the CG method to achieve convergence. Furthermore, many researchers have extended the application of FEM to the simulation of electromagnetic scattering from non-spherical objects. Xu et al. used a two-dimensional scattering model of the FEM to simulate microwave scattering from a layered medium [113]. Khodapanah modified the basis vectors of the FEM for solving the electromagnetic scattering problem of continuously varying nonuniform spherical particles [114]. As the spherical harmonics waves display an orthogonality in the angular direction, the FEM can be applied to generate extremely sparse (multi-diagonal) matrices, which effectively improves the computational efficiency. Additionally, the FEM avoids the problem of singular kernel in the integral equation method. However, it generates a large computation amount as its simulated domain is not limited to the scatter itself, and its numerical simulation accuracy is uncontrolled given the influence of the absorption boundary conditions [115].

D. APPROXIMATE METHODS
The approximate calculation model is primarily suited to scenarios in which the incident light wavelength significantly differs from the particle size. When the incident light wavelength is substantially greater than the particle size, the classical models such as the Rayleigh approximation and the Rayleigh-Gans-Stevenson approximation can be employed. Conversely, when the incident wavelength is significantly smaller than the particle size, the anomalous diffraction approximation, the geometric optical approximation, and the physical geometric optical approximation are applicable.

1) RAYLEIGH APPROXIMATION
The Rayleigh scattering approximation (RA) model, first proposed by Rayleigh in 1897, concentrates on particles with size much smaller than the wavelength of incident light. To solve this model, it is necessary to assume that the internal and incident fields of the particles can be approximately treated as electrostatic fields and the internal fields of the particles are uniform. Although complete analytical solutions exist for the RA model for some particles with simple shapes (e.g., ellipsoidal), it is still necessary to consider the integral equations of polarization vectors for particles with complex shapes. Il'in et al. extended the Rayleigh approximation to non-spherical scatterers and validated it using Chebyshev particles [116]. Athanasiadis et al. utilized the low-frequency formula of the Rayleigh approximation to solve the inverse electromagnetic scattering problem using the near field data for dielectric ellipsoids [117]. Recently, in 2020, the method was expanded to laminar ellipsoids [118] and elastic ellipsoids [118] and elastic ellipsoids [119]. The Rayleigh approximation is applicable to the calculation of electromagnetic scattering from uncharged particles whose radii are much smaller than the incident wavelength; however, it is not as efficient as Mie theory in calculating electromagnetic scattering from small charged particles [120].

2) RAYLEIGH-GANS-STEVENSON APPROXIMATION
The Rayleigh-Gans-Stevenson (RGA) approximation proposed by Stevenson is an improvement on the Rayleigh approximation [121], which increases the expansion order of the basis function for the incident and scattered field to the fourth power of the size parameter (in contrast to the quadratic expansion of the Rayleigh approximation). To obtain an analytical solution, it is necessary to assume that the scattering from elements in the scatterer is only excited by the incident field and independent of other factors in the model. Currently, the model is mainly applicable to optical ''soft'' particles with small particle size (where the refractive index m satisfies |m − 1| = 1, i.e., the refraction and reflection of light at the particle boundary is almost negligible), and the range of its applicable size parameter is wider than that of the Rayleigh approximation. Lu et al. improved the algorithm of the GGA, which was tens of times faster in calculating the scattering properties for different incidence angles than the A-DDA [122]. Hogan and Westbrook proposed the Self-Similar Rayleigh-Gans Approximation (SSRGA) to describe the ''self-similar'' structure of snowflakes by means of a power law to calculate the backscattering cross-sectional equations of microwaves for snowflakes [123]. In 2017, the SSRGA method was developed to calculate the scattering phase functions for immature aggregates with a size of 1 cm [124]. Leinonen et al. evaluated the applicability of RGA and SSRGA for edge snow microwave scattering and demonstrated that the SSRGA was able to eliminate the small constant bias in RGA with good accuracy [125]. It is acknowledged that the use of RGA always leads to a systematic underestimation of the scattering and the absorption, as well as an inability to predict polarization properties. McCusker et al. proposed a new method that saved much time compared to the DDA method by involving intra-monomer coupling but ignoring the inter-monomer coupling, and had a much better accuracy than the conventional RGA [126]. Following continuous developments by numerous researchers, the RGA model can now be applied to calculate the scattering of circularly symmetric particles, anisotropic particles, and particles with Gaussian random surfaces [127], [128].

3) ANOMALOUS DIFFRACTION APPROXIMATION
The Anomalous Diffraction Approximation (ADA), which takes into account refraction in addition to diffraction, was initially theorized by Hulst and van de Hulst and was mainly utilized to solve light scattering problems for large, optically soft particles [129]. This method can provide fast and fairly accurate approximations for such particles. In subsequent years, the ADA has been extended to encompass more complex particle shapes and to the semi-empirical calculation of light scattering from particles with a smaller size parameter and larger refractive index. The applicability of ADA was examined by Kokhanovsky [130]. Jacquier and Gruy then presented a chord length distribution to calculate the light scattering cross section of aggregates. This approach was applicable to a broad range of ordered and disordered aggregates and was more than 100 times faster than the conventional ADA [131]. Subsequently, Mitchell et al. proposed the Modified Anomalous Diffraction Approximation (MADA), which offers a straightforward ADA expression for an arbitrary particle shape, allowing for simpler and more accurate calculations in comparison with ADA [132]. In geometric optics, MADA is utilized for the simplified calculation of extinction and absorption properties of particles in the atmosphere. Brendel et al. discussed the applicability of a ADA method to calculate scattering from infinite cylinders [133]. Thomas then derived an improved ADA formula by relaxing some approximations and considering the polarizability response of particles to an applied magnetic field [134]. This enhanced method enables more precise prediction of the extinction and absorption efficiencies for particles covering the entire range of size parameters.

4) GEOMETRIC OPTICS APPROXIMATION
The Geometric Optics Approximation (GOA), with particular focus on the improved version based on ray-tracing techniques, is a common method for calculating light scattering from arbitrarily oriented particles with any shape at scales much larger than the incident wavelength. Cai et al. were the first to consider the polarization effects and phase interference in ray tracing techniques, and developed vectorbased geometric ray tracing algorithms for several different coordinate systems [135]. Although the conventional GOA disregards the particle absorption effect, this is questionable for strongly absorbing particles. To incorporate the absorption effect in simulations, Stratton and Born et al. expanded the traditional concept of refractive index to complex refractive index. As a result, the imaginary part of the refractive index can describe the absorption properties of the particles, and new Fresnel reflection and refraction coefficients were defined to construct a more general GOA method that takes the absorption effect into account [136], [137]. Takano et al. utilized the GOA method in combination with the RGA method to calculate the single scattering properties of eight black carbon fractal agglomerates consisting of 7-600 primary spherical particulates [138]. In comparison with the results obtained through the superposition TMM, the GOA method was found to be effective in calculating light scattering from agglomerates comprising many primary spherical particulates (>6000) with large size parameter (≫2). To address complex particle problems, the application of the Monte Carlo method (MC) in light scattering calculations can be combined with a geometric ray tracing algorithm to accelerate computation speed [139]. Although the GOA method is theoretically an approximate one, it is important to consider the applicability range of the minimum size parameter χ min when using it in practice.

5) PHYSICAL GEOMETRIC OPTICAL METHODS
A new algorithm combining the near-field analytical beamtracking technique and the far-field accurate mapping technique was named the Physical Geometric Optics Method (PGOM) to distinguish it from GOA [140]. The accuracy of PGOM, however, has not been well quantified because ray tracing techniques only provide approximate values. Bi, L. employed PGOM and II-TM+IGOM to calculate the optical properties of ice cloud bodies, with the calculated results of II-TM+IGOM serving as a benchmark to estimate the uncertainty of PGOM in remote sensing and radiative transfer simulations [141]. Subsequently, Sun et al. proposed a new PGOM that can be applied to particles with arbitrary surfaces, significantly reducing the computation time to 10-30 percent of the original PGOM, depending on the particle shape and refractive index. This was achieved by dividing the original beam into several sub-beams, so that each sub-beam would only be incident on a specific surface [142]. Zhu et al. further employed the beam-tracking technique from the GOA method to obtain the analytical expressions for the absorbance and scattering albedo of multilayer particles by quantifying the analytic radiation absorption and scattering of the particles [143]. Finally, Ding et al. determined the limit of geometric optical ray tracing technique by numerically solving the vector Kirchhoff integral equations [144].

IV. MEASUREMENTS OF AEROSOLS A. SAMPLING MEASUREMENTS OF AEROSOLS
The AOPs are characterized by the scale, shape, composition, number concentration, mass concentration and surface area concentration of aerosol particles. In the past two decades, the simulation algorithms of atmospheric aerosols have displayed rapid advances, though the simulation results are still typically treated as a reference for experimental observations and out-door measurements. The most crucial particle parameters for simulations are obtainable only through measurements; thus, only experimental observations and measurements can yield representative and reliable information about aerosols. For this reason, it is pertinent to gain a comprehensive understanding of the fundamental principles, measurement accuracy and applicable conditions of the measuring methods. However, a comprehensive review of all modern sampling measurement methods is far beyond the scope of this paper. This section provides an overview of direct measurement techniques of the scattering and absorption properties of atmospheric aerosols, emphasizing the in-situ measurements of AOPs. Additionally, some aerosol detection and characterization techniques are briefly discussed, because they are relative to the microscopic optical properties of aerosols.

1) MEASUREMENTS OF SCATTERING PROPERTIES
The single scattering albedo of aerosols can be determined experimentally by measuring the scattering coefficients of the aerosols, or it can be inverted through various remote sensing methods. However, direct measurement of the aerosols' scattering coefficients in a laboratory environment is more accurate than remote sensing methods, and is thus regularly employed to verify the results obtained from remote sensing inversion.
Integral turbidimeters are instruments that measure the forward or backward scattering of aerosols. However, the limitations of the instrument's geometry prevent the capture of the full range of scattered radiation, which is typically distributed between a forward scattering angle of 7 • -170 • and a backscattering angle of 90 • -170 • . Therefore, the missing scattering angles in the front and back directions need to be corrected, as well as other corrections to eliminate the effects of molecular scattering, instrument wall scattering, and parasitic light. When equipped with a photon counting detector, an integrated turbidimeter can measure an extremely small light scattering coefficient of aerosols, which is even less than 1% of that for clean air at normal atmospheric pressure. However, there are various sources of error associated with the use of an integral turbidimeter, such as: 1) an underestimation of the contribution of coarse particles (particle size > 5 µm) due to their tendency to be deposited at the inlet, with the inlet loss increasing with size; 2) the inability to measure nearly forward scattering light (from 0 • or 5 • to 10 • , depending on the instrument design); and 3) evaporation of droplets induced by aerosol heating, which may introduce a very large error, especially in high relative humidity (>90%) where water constitutes the majority of the particle volume.
The integral turbidimeter, initially proposed by Beuttell and Brewer [145], is a series of small-scale instruments designed to measure visibility by directly measuring the scattering coefficient of a relatively small volume of air. Heintzenberg and Charlson provided an in-depth review of the design philosophy, instrumentation principles, possible designs, calibration, systematic errors, application to scientific problems, and inherent limitations of integral turbidimeters [146]. The Umov effect was further studied quantitatively, revealing a negative correlation between the detected maximum of linear polarization and the reflectance near backscattering when light is scattered from the particle surface [147], [148]. Thus, the application of polarization properties for measuring the light scattering properties of aerosols has become a new technology, with polarization turbidimeters providing better portability and a wider measuring angle. Dolgos and Martins developed a portable polarization imaging turbidimeter (PI-Neph) with the optical system and sensor layout as shown in Fig. 6, which captures the P11 and P12 matrix elements by taking two consecutive pictures. Their measured results for artificial spherical aerosols were found to be consistent with the expectations from Mie theory [149]. In 2016, Espinosa et al. used the Polarized Imaging Nephelometer to achieve high-precision insite measurements of the phase functions and degree of polarization over a wide angle of 3 • -177 • at three visible wavelengths [150]. In addition, laser imaging turbidimeters also have great application prospects. In 2018, Manfred et al. used 375 nm and 405 nm diode lasers to measure the nonpolarized scattering phase function of aerosols, and their results exhibited a high angular resolution of 0.5 • over the scattering angle of 4 • -175 • [151]. Moreover, Qiu et al. proposed a machine learning method for correcting the scattering coefficient obtained from a three-wavelength turbidimeter in different relative humidity, with more than 85% of the corrected values displaying errors less than 2% in 2021 [152]. Teri et al. evaluated the effectiveness of angle correction for different turbidimeters through experiments and simulations, developing guidelines for the optimal angle correction with respect to aerosol types and their size ranges [153]. Cavity spectrometer: Cavity spectrometer is a type of photon counter and divided into active cavity and passive cavity according to the probe. By measuring the visible VOLUME 11, 2023 radiation caused by the scattering of aerosol particles, the equivalent size distribution of particles within a certain size range can be estimated after the refractive index of aerosol particles is obtained. Varma et al. compared the performance of three broadband cavity spectrometers used to measure the extinction coefficients of aerosols and found that the measured results were consistent with those calculated according to Mie theory, demonstrating the great potential of these instruments in investigating AOPs [154]. Bluvshtein et al. performed the first measurements of extinction coefficients in the UV range from 315 nm to 345 nm using cavity enhancement spectroscopy [155]. Min et al. applied a dualchannel broadband cavity enhanced absorption spectrometer (BBCEAS) mounted on an aircraft to measure the AOPs at a high altitude [156]. Wang et al. developed a portable cavityenhanced absorption spectrometer that enabled fast and stable operation for practical applications [157]. In contrast to earlier cavity spectrometers, which were limited to operating at fixed wavelengths or measuring a single optical parameter, Xu et al. designed a three-wavelength cavity-enhanced photometer, as illustrated in Fig. 7. The light exiting the cavity was focused into an optical fiber before being coupled into a CCD spectrometer with a 100 µm-wide slit, thus enabling multi-wavelength, multi-parameter measurements of AOPs [158]. Since most extinction measurements require a stable light source to attain high accuracy and precision, Chen et al. proposed a method to enhance the stability of the light source, thereby improving the accuracy and precision of the broadband cavity system [159].

2) MEASUREMENTS OF ABSORPTION PROPERTIES
The light absorption of aerosol particles and air molecules exert considerable influences on the radiative properties of aerosols and atmosphere, necessitating accurate measurement of the absorption coefficient. While standard methods enable the concentration of gas to be easily determined, there is no analogous access to determine the light absorption of aerosol particles. Furthermore, these properties vary depending on the particle type and size, making it difficult to measure the absorption coefficient of aerosol particles. Moreover, all aerosol particles are capable of scattering light, making it almost impossible to completely eliminate the influence of scattering when calculating the absorption coefficient of particles. Nonetheless, there are still several common methods for measuring aerosol absorption coefficients σ abs : Based on photoacoustic (PA) spectral effect: When light pulses exhibiting a repetition rate within the range of acoustic frequencies pass through aerosols, the absorbed light energy heats the air through photothermal conversion, producing acoustic waves. The light absorption coefficient of aerosol particles can be determined by measuring the amplitude of the corresponding acoustic waves. At present, the photoacoustic method is a commonly used, highly sensitive technique for real-time aerosol absorption coefficient measurements. Careful calibration of the chopping frequency is important, as the acoustic signal is affected by the thermal relaxation processes in the particles. Arnott et al. provided a detailed description of the principles of calibration instruments [160], while Li et al. reviewed the development of infrared photoacoustic spectroscopy (PAS) techniques in the field of gas-phase analysis and presented their applications [161]. Cremer et al. employed a single-photon photoacoustic spectrometer to make highly sensitive measurements of aerosol particle size [162]. One of the challenges when measuring aerosol photo absorption using the photoacoustic method is the light absorption by atmospheric gases, such as NO2 and O3. To accurately measure the light absorption of particulate matter from the UV to the IR bands, the influence from gas-phase absorption must be eliminated [163]. Additionally, the lack of aerosol light absorption standards necessitates the use of light-absorbing gases for the calibration of the response to photoacoustic signals [164]. Foster et al. calibrated a multichannel PAS using small particles with strong light absorption (SSA < 0.5), and the standard deviation of the calibrated slope was less than 2% at 660 nm and less than 5% at 405 nm [165]. Yu et al. proposed a three-wavelength differential photoacoustic spectrometer (RGB-DPAS), whose experimental setup is depicted in Fig. 8 [166]. This instrument measures the aerosol absorption coefficients directly from the normalized signal responses (NSR) without needing to calibrate the photoacoustic signal response. The three-wavelength photoacoustic signals of aerosols were recorded by two identical photoacoustic units, which reduced interference from atmospheric gas absorption and low frequency ambient acoustic noise. Cao et al. proposed a three-wavelength photoacoustic spectrometer (TW-PAS) to address the reduced sensitivity of multi-band PAS operating under non-resonant states. This instrument could achieve a high temporal resolution and a high sensitivity because it operated in the resonant mode of each wavelength [167]. Filter-based methods: Filter-based methods are the most widely used technique for the sampling of atmospheric aerosols due to their low cost and convenience. However, they tend to overestimate the absorption coefficient, which is caused by two factors. Firstly, the particulate light scattering reduces transmission through the loaded filter. Secondly, the filter itself causes the light scattering, allowing a single photon to traverse the absorbing particle layer multiple times. To improve the absorption coefficient, correction coefficients based on empirical measurements or radiative transfer considerations of the reference system have been applied [10], [168].
The integrating sphere technique is a method proposed by Fischer for measuring the mass absorption coefficient of atmospheric aerosols with diameters of 0.4-2.4 µm. The detector simultaneously collects the transmission light through the sample and the reflection light from the sample and the signal strength is proportional to the sum of these two types of light. Therefore, the absorption coefficient of a particle is derived by comparing the signal difference between a pure substrate and a substrate loaded with particles [169]. The integrating plate technique was a simplified method based on the integrating sphere technique, first developed by Lin et al. [170], which can be used to measure the absorption coefficients of granular deposit. Montilla et al. designed a variable-resolution integrating sphere spectroscopy system to measure the light absorption coefficients of atmospheric aerosol in the spectral range of 320-800 nm [171]. Wang et al. developed a two-sphere integration (TSI) technique, illustrated in Fig. 9, which combines two separate integrating spheres to reduce the scattering of light-absorbing particles and thus accurately determine the total light absorption of particles collected on a nuclear pore filter [172].
The extinction coefficient of aerosol is measured using an extinction meter or a far-light photometer, and the light scattering coefficient is measured using an integral turbidimeter. Subsequently, the difference between the two measured results is treated as the absorption coefficient of the aerosol [173], [174]. However, for aerosols with low absorption, the difference between the extinction and scattering coefficients is small, requiring both values to be measured with high precision. To achieve this, either improved truncation errors or special calibrations for the integral turbidimeter that have already been corrected for truncation errors is necessary. Dial et al. developed an enhanced albedometer by implementing a cavity ring-down (CRD) technique. This strategy reduced the sample volume to less than 1% of the original design and shortened the response time by more than 30 times [175]. Onasch et al. introduced a single scattering albedo monitor, which enabled the simultaneous measurement of the extinction and scattering coefficients of aerosol particles in air in the spectrum range from 450 nm to 780 nm [176]. As errors were easily amplified by a subtraction operation, every source of error must be carefully checked before using an EMS method. Modini et al. developed a new modularized framework for the single scattering albedo monitor (illustrated in Fig. 10), and used it to evaluate the calculations of truncation error correction [177]. Singh et al. studied the errors and uncertainties of AOPs from cavity oscillation spectroscopy, integral turbidimeter method and EMS method. The experimental results showed that the overall uncertainty in the extinction cross section was 10 -11%, and the error was primarily from the measuring error of the condensed particle counter (CPC) [178].

3) OPTICAL CHARACTERIZATION METHOD OF AEROSOLS
Many optical analysis methods can characterize the biochemical components of atmospheric aerosols according to the interaction between light and some typical chemical molecules. On the one hand, these optical methods can detect the components of aerosols, on the other hand, they provide some effective paths to obtain the microscopic optical properties of aerosols which construct the fundament of macroscopic AOPs (such as absorption spectrum, infrared radiation spectrum and Raman spectrum). Thus, this part describes some common optical characterization methods of aerosols.
Spectroscopy of inelastic scattering: Laser-Induced Fluorescence (LIF) and Raman Spectroscopy (RS), have been commonly used in the real-time detection of Primary Biological Aerosol Particles. LIF is mainly used to analyze the fluorescence characteristics of single particles, and has been extensively employed for the identification and characterization of biological and other organic carbon aerosol particles. RS, also referred to as spontaneous Raman scattering, is an ideal method for rapidly characterizing atmospheric aerosol particles, with advantages such as molecular specificity and high resolution of spectrum. RS has been applied to the detection of biological samples in aerosols, and used to analyze the characteristics of atmospheric aerosols at the single particle level. Numerous reviews and comprehensive studies have discussed the applications of these techniques in aerosol characterization. For instance, Pan et al. reviewed the use of LIF to detect and characterize the biological and other organic carbon aerosol particles in the atmosphere [179]. Estefany et al. gave a detailed introduction to the analysis methods and detecting principles of different aerosol components (inorganic salts, organic substances, smoke dust, etc.) using RS, and comprehensively reviewed the progress of Raman spectroscopy in the microscopic characterization of atmospheric aerosol [180].

B. REMOTE SENSING MEASUREMENTS OF AEROSOLS
Measuring the characteristics of aerosols over a wide space in situ is difficult, so synchronous development of measurements based on remote sensing techniques was needed to address these problems. The light scattering and absorption of aerosol particles can significantly modify the properties (e.g., intensity, polarization, and phase) of incident radiation. Therefore, it is mathematically possible to invert the aerosol optical properties (AOPs) by constructing a correlation model between the modified properties of incident radiation and the intrinsic properties of aerosols, which is the fundamental principle of remote sensing measurements of aerosols. Remote sensing methods used to measure AOPs can be divided into ground remote sensing and satellite remote sensing according to the working condition of equipment, as well as active remote sensing and passive remote sensing according to the detection mechanisms. Because active remote sensing techniques are mainly used in ground remote sensing, we further divided ground remote sensing into three parts: passive remote sensing, active remote sensing and joint measurements, while satellite remote sensing is focused on the passive remote sensing.

1) GROUND REMOTE SENSING
The most accurate aerosol characteristics can be acquired through ground remote sensing, as it does not require the dealing with reflectance issues of ground. Such techniques have been developed over the past few decades, providing valuable data to validate satellite-based observed results and models. According to different detection mechanisms, the techniques of ground remote sensing are divided into passive remote sensing and active remote sensing.
Passive remote sensing: Remote sensing techniques utilizing a photometer: Utilizing a solar photometer to detect AOD is currently the most precise technique in ground-based remote sensing, and is often utilized to corroborate the detection results acquired from satellites. This approach is based on inverting the AOD and the distribution of particle size according to the atmospheric extinction spectrum, which is acquired from a series of narrow-band filters (typically with a full width at half maximum of approximately 10 nm) operating in the visible to near-infrared wavelength range [183]. Fig. 11 illustrates a schematic diagram of various passive remote sensing techniques for observing atmospheric aerosols. In practical applications, errors may arise from multiple sources, such as the operating waveband of photometer, half-wave width, instrument calibration (Langley calibration method), atmospheric absorption, and observation conditions. To further explore the applications of photometer, several studies have been conducted. For example, in 2002, NASA Goddard Space Flight Center (GSFC) initiated research on multi-band solar photometer and calibrated  [195]; (b) A space-borne High-spectral-resolution Lidar [196].
the aerosol optical thickness from measured signals at four wavelengths (340, 440, 675, and 870 nm), yielding more reliable results than with Langley calibration [184]. Li et al. derived most key features of atmospheric aerosols, including single scattering albedo and scattering matrix, across a spectrum from ultraviolet to near infrared by means of polarized solar photometer [185]. Zhang et al. systematically evaluated the effect of photometer's filtering function on aerosol optical depth retrieval [186]. Moreover, since the singleangle detection technique is unable to differentiate between scattered and absorbed components, the multi-angle detection technique and multi-spectral polarization detection technique have become increasingly popular areas of research in remote sensing. Si et al. summarized the typical instruments used for atmospheric aerosol detection in different countries through aerosol characterization with remote sensing measurements. Taking the multi-angle imaging spectrometer (MISR) on Terra as an example, they reviewed the research progress of aerosol optical properties from different aspects, including data set accuracy, inverse algorithm improvement and product application [187]. To investigate techniques for characterizing the vertical distribution of aerosols in the planetary boundary layer (PBL) using passive remote sensing, Choi et al. analyzed the aerosol optical properties of the Los Angeles Basin by measuring radiation and polarization spectrum in the near-infrared band and applying the Observation System Simulation Experiment (OSSE) to estimate the reflected near-infrared solar radiation information from various sensors [188].
Several major solar photometer networks exist worldwide. For example, NASA's AERONET and the French PHOTONS networks can provide accurate original data and high-quality inverse results regarding Aerosol Optical Properties (AOPs) over a large area in cloud-free conditions. As an illustration, an inverse algorithm was developed for AERONET to infer the aerosol single scattering albedo and particle size distribution in combination with the radiation data measured by a solar photometer [189].
Another instrument widely used for the measurement of AOD in visible range is the Multi-Filter Rotating Shadowband Radiometer (MFRSR) [190]. MFRSR obtains the direct solar irradiance by subtracting the shadowed irradiance from the measured global irradiance. A new algorithm allowed MFRSR to partition the aerosol optical thickness into fine aerosol mode and coarse aerosol mode and it can also retrieve the effective radius and the Angstrom Exponent of the fine mode [191]. Alexandrov et al. provided a detailed review of the latest advances in the use of MFRSR for characterizing atmospheric aerosols and in the construction of MFRSR networks [192]. As a supplement of AERONET, MFRSR networks usually provide better spatial density of measurement sites. Chen et al. provided an exhaustive review of the calibration methods of MFRSR [193]. Di Sarra et al. combined MFRSR with Cimel sunphotometer to develop an MFRSR-Cimel dataset. After correction, the mean bias (MB) between simultaneous AOD determinations by MFRSR and Cimel was always below 0.004, with the root mean square difference ≤0.031 at all wavelengths [194].
Active remote sensing: All of above passive remote sensing techniques utilize solar radiation as the light source, while Lidar (Light Detection and Ranging) uses laser as a light probe to obtain the vertical distribution information of aerosols. Lidar recovers AOPs by measuring the intensity of the backscattered light following the collision between the laser pulses and aerosol particles. The recorded backscattering signal is expressible as a function of time, thus allowing for the precise estimation of the aerosol height from the flight time of light. Laser radar remote sensing can be divided into two categories: one is to place a laser radar on the ground for aerial detection, and the other is to load a laser radar on aircraft for ground detection. Most aerosol lidars operate in the visible and near-infrared wavebands, as the aerosol optical depths (AODs) in these wavebands are sufficiently large to generate intense backscattered light. Fig. 12 (a) illustrates a novel vibrational Raman lidar system that employs a neodymium-doped yttrium aluminum garnet (Nd:YAG) pulsed laser as the light source. The reflected lidar signals are collected by a telescope, after which they are coupled into a multimode fiber and redirected into a newly-developed spectrum system.
Commonly used lidars for aerosol detection include Mie scattering lidar, polarization lidar, Raman lidar, and high spectrum resolution lidar. Mie scattering lidar is the most widely-utilized lidar for aerosol detection. Polarization lidar adds a polarization prism to the receiving system, which is based on the Mie scattering lidar. Polarization remote sensing can acquire polarization data of a specific target and quantify the polarization characteristics of the target based on its intrinsic properties; thus, a polarization lidar can obtain more atmosphere information than a Mie scattering lidar. Yan et al. provided an extensive review of recent technological advances in optical polarization remote sensing and their new application areas [197]. Although the structure of the Mie scattering lidar is simple, its reversion error can be more than 30% due to inaccurate reference point and lidar ratio. To address this limitation, techniques for direct detection of aerosol extinction coefficients are modified by adding Raman channels or Rayleigh hyperspectral resolution channels. Ceolato and Berg gave a comprehensive review of these measurements' principles, modeling, and applications [198].
Joint measurement: It is worthwhile to pay attention to the joint measurement method which combines the active and passive remote sensing. It is well known difficult when lidar is used alone to retrieve the AOD of aerosols without any knowledge of lidar ratio (extinction to backscatter ratio). However, passive remote sensing exhibits some capacity to estimate the vertically average value of lidar ratio. Consequently, a joint measurement of active and passive remote sensing had become an important strategy to retrieve some passively measured optical parameters (such as AOD) and certain optical parameters (such as lidar ratio) [199]. In addition, the relationship between aerosol vertical profile information (acquired from lidar) and passive measured optical parameters is a critical issue for understanding the complex aerosol layering phenomena and the resulting atmospheric effects [200]. Hence, a joint measuring method combining passive and active remote sensing is the most useful approach to improve our knowledge about the distribution and variation of atmospheric aerosols [201].

2) SATELLITE REMOTE SENSING
Because of the short life, complex components and extensive distribution of aerosols, satellite remote sensing is the only way to observe aerosols on a global scale. Fig. 12(b) presents a new, space-borne High-Spectral-Resolution Lidar (ACHSRL) designed for aerosol and cloud optical property profiling, which can precisely retrieve particle optical properties. To effectively detect aerosol optical properties (AOPs), several advanced sensors were employed for global aerosol monitoring. For instance, Mishchenko et al. used the Advanced High-Resolution Radiometer (AVHRR) to retrieve the global monthly average optical thickness and Angstrom index of tropospheric aerosols, and dedicated years to monitor the aerosol optical thickness, thus rationalizing the overall decreasing trend in global tropospheric aerosol optical thickness from 1991-2005 [203]. Li et al. constructed a dataset consisting of long-time observations of AOD on land by intercalibrating AVHRR results with radiation data from the widefield-of-view SeaWiFS sensor [204]. Waquet et al. applied a scanning cyclotron for aerosol detection on land, achieving accurate polarization measurements over a wide spectrum (410 -2250 nm) and large angular range (±60 • ) [205]. Dubovik et al. provided a comprehensive summary of instruments, methods, results, and perspectives of polarization remote sensing for atmospheric aerosol [206]. NASA carried the Moderate Resolution Imaging Spectrometer (MODIS) on the Terra and Aqua satellites for daily global aerosol observations over a broad spectrum (0.41 -15 µm), leading to numerous achievements, including the AOT data at three visible wavelengths on land and seven wavelengths from 0.47 -2.13 µm above the ocean [202]. Kahn et al. statistically compared the aerosol data generated from the Multi-angle Imaging Spectroradiometer (MISR) with that from MODIS, successfully reverting AOPs at more than 75% of randomly selected regions. The comparison showed that the correlation coefficients between MISR and MODIS were about 0.9 above the ocean and 0.7 on land [207].

V. SUMMARY AND PROSPECTS
AOPs have been studied by numerical simulations, sampling measurements, ground-based observations, and satellite remote sensing, the effectiveness of which methods can be validated against each other. In the last few decades, more and more progresses have been made in the fields of numerical simulation to obtain the scattering properties of spherical particles and non-spherical particles. This paper reviews the exact and approximate simulation methods calculating the scattering properties of various particles. In terms of current computational power, strict numerical calculations are still practically confined to a specific range of dimensional parameters. Consequently, a more flexible approach is considered to combine accurate numerical simulation, such as FDTD method, with random algorithms or geometric-optical approximations. This expectation is prospective to realize in the near future. Besides that, this paper also describes the advances in instrument-based sampling measurements for aerosols, and the obtained accurate information from them will eventually enable theoretical models to better describe the intrinsic optical properties and predict the integral optical properties of atmospheric aerosols. Traditionally, sampling measurements use filter-based techniques to collect particles for subsequent analysis which leads to a striking disadvantage that the scale, amount, and morphology feature of aerosol particles are changed inevitably during sampling. The question, therefore, is how to use representative samples to characterize actual atmospheric aerosols.
While some random algorithms can give several limited integral optical properties, the most useful and accurate methods are remote sensing measurements. Most current remote sensing techniques hold the limitation of keeping the consistency of results collected by different inversion algorithms and sensors. Because amounts of sensors have come into service on the ground or in space, and more sensors are also under development, there are some urgent requirements for more accurate and robust inversion algorithms and more representative aerosol models to better infer the integral optical properties of aerosols. A feasible approach is the combination of multidimensional optical detection, such as broadband absorption and scattering spectral, polarization states and light intensity distribution, which is helpful to limit the freedom degree of atmospheric inversion model and expected to reduce the uncertainties of AOPs through data assimilation.