Simultaneously simulating the scattering properties of nonspherical aerosol particles with different sizes by the MRTD scattering model

Shuai Hu; Taichang Gao; Hao Li; Ming Chen; Feng Zhang; Bo Yang

doi:10.1364/OE.25.017872

1. Introduction

No matter in climate modeling or remote sensing applications, scattering properties of aerosol particles are extremely important since they act as the basic parameters of Radiative Transfer numerical Models (RTM) [1–4]. However, because of the complicated shapes and inhomogeneous compositions of aerosols, there are still substantial uncertainties existed in their scattering properties [5–7], thus further limiting the accuracy of radiative transfer calculations and remote sensing [8–10]. For polarized remote sensing, the influence of aerosol’s non-spherical shape is much more remarkable since the polarized components of electromagnetic wave are much more sensitive to the particles’ geometries [11]. Owning to these reasons, the scattering characteristics of the aerosol particles with irregular shapes have been being widely investigated at present.

To calculate light scattering by particles with arbitrary shapes, many numerical approaches have been developed, such as the T-Matrix Method [12,13], Invariant Imbedding T-matrix method [14,15], SVM (Separation of Variables Method) [16], PMM (Point Matching Method) [17], Method of Moments [18] and Discrete Dipole Approximation (DDA) [19]. With the development of computational electromagnetics, such methods as Finite Difference Time Domain (FDTD) [20] and Finite Element Method (FEM) [21] are also introduced into the light scattering simulation of non-spherical particles. Compared with other models, their numerical calculation are not limited by the boundary and initial conditions, and the singularity of volume integral methods can be avoided, thus they are more suitable for particles with complex shapes and large sizes. However, owning to the constraint of numerical dispersion, the spatial resolution for these models can’t be too large [22], so their computational amount would be huge for large particles. To alleviate this difficulty, Pseudo Spectral Time Domain (PSTD) is introduced into the numerical simulation of light scattering process [23,24]. Since PSTD can achieve high simulation accuracy with less girds per wavelength, thus less computational time and memory are needed [23–25]. Similar to PSTD, MRTD is also an excellent computational technique widely used in electromagnetic computational field, and now it has been applied to the light scattering simulation of non-spherical aerosol particles by our team [26]; to improve the computational efficiency, the MRTD model is further parallelized by MPI repeated non-blocking communication technique [27].

Compared with FDTD, though MRTD can save computational time to some extent, it can only calculate the scattering properties of one particle size (with specific orientation) at the given wavelength in one wave-particle interaction simulation [26], therefore, if a large number of particle sizes are needed to simulate, the computational process will still take a long time. This problem is indeed encountered in the field of atmospheric radiative transfer. Since the scattering parameters in RTMs are volume-averaged values (i.e., the scattering properties of aerosol should be integrated and normalized over their size distribution and orientations [1,28]), the computational process should be repeated for many times in order to obtain the scattering parameters of the randomly-orientated particles with given size distributions, which is a very time-consuming process. So it will be meaningful if the MRTD model can be improved so that the scattering parameters of particles with different sizes can be simultaneously simulated in one wave-particle interaction simulation. According to the scattering theories, particle radius r in the theoretical formulation of the light scattering problem is only encountered as a ratio to the wavelength (i.e., size parameter x, x = 2πr/λ, λ is the wavelength) and similar results can be achieved by increasing particle radius or reducing the light wavelength by the same ratio [29,30], therefore the light scattering simulation for particles with different sizes can be transformed into the scattering calculation of a size-fixed particle at different wavelengths once the shape, orientation and refractive index of the particle are given. Based on the consideration above, the concept of scattering transfer function and the method to transform transient field to time-harmonic field in computational electromagnetics are introduced into MRTD scattering model, and the Multi-size Synchronous-Computational Scheme (MSCS) is proposed, by using this method, the scattering process for incident light with a wide spectrum can be simulated by the IMRTD model, and the scattering parameters of particles with different sizes can be obtained by normalizing the scattering properties at different wavelengths.

Our work in this paper is divided into 4 parts, in section 2, the basic principles and frame of MSCS and IMRTD scattering model are introduced, and the method to design the model’s input parameters, such as the spatial resolution and discrete time interval, is proposed as well. In section 3, the precision of IMRTD model is validated against well-test scattering models, and the influence of spatial resolution on the simulation accuracy is discussed, besides, model’s computational efficiency is also discussed at the end of this part. In section 4, a brief summary of this study is given.

2. Principles and methodology

2.1 Basic principle of the MSCS

According to the theory of computational electromagnetics, aerosol particles can be regarded as linear scattering medium, and the light scattering process by aerosols can be taken as a linear transfer of the incident light [22]. Once the refractive index, size, shape and orientation of a particle are given, then the signal transfer function G(λ) of the scatterer in frequency domain is also fixed, and the transformation relationship of incident field E_inc and scattering field E_sc can be written as

E_{s c} (λ) = E_{i n c} (λ) \cdot G (λ) .

From the equation above, it can be found that if the incident light with a wide spectrum is scattered by aerosol particles, the scattered light should be the linear synthesis of the time-harmonic scattering wave at the frequencies in this spectrum correspondingly; therefore, if the time-harmonic components of the scattering and incident field at different frequencies are extracted, the scattering properties at different wavelengths can be obtained as well. Based on the idea above, the Multi-size Synchronous-Computational Scheme (MSCS) is proposed by introducing the transformation technique of transient to time-harmonic field into MRTD model, and its basic principle is shown in Fig. 1. As presented in the figure, the pulse plane wave with a wide spectrum is injected into the computational domain, and the transient near electromagnetic field is calculated by MRTD scheme; then, the time-harmonic field components of different frequencies are extracted by Fourier Transform (FT), and the scattering properties at different wavelengths are computed; at last, the scattering parameters of particles with different sizes are calculated by normalizing the obtained results.

Fig. 1 The principle of the Multi-size Synchronous Computational Scheme for the IMRTD scattering model.

Download Full Size | PDF

The detailed computational flow of the MSCS is further given in the Fig. 2. As described in the figure, the calculation process of the MSCS is divided into the following four steps:

Fig. 2 The flowchart of the Multi-size Synchronous Computational Scheme for IMRTD model.

Download Full Size | PDF

(1) The design of pulse wave and discrete grid in time and space domain: to obtain a wide frequency spectrum, the incident light is usually taken as Gauss pulse plane wave, modulated Gauss pulse and raised cosine pulse. The key step to design the pulse incident wave $E_{G P}^{i} (t)$ is to determine its upper cut-off frequency f_max since it decides the particle size range that the scattering model can effectively simulate. Assuming that we want to obtain particles’ scattering properties at the wavelength λ₀ (called “simulation wavelength” in this paper), and the radius of the particles needed to simulate is ranging from r_min to r_max, then f_max can be calculated by

f_{\max} = \frac{r_{\max}}{r_{0}} \cdot \frac{c}{λ_{0}},

where, r₀ is the radius of the size-fixed particle set in the computational space (called “original particle radius”), and c is the velocity of light. Once f_max is given, then other parameters of the pulse wave (like the pulse width, upper-limit frequency, etc.) can be determined as well (as shown in section 2.3).

Discrete grid size in MSCS is not only a parameter determining the accuracy of particle construction and iterative process, but also a factor influencing the upper cut-off frequency f_max of the effective signal, therefore its design is also an important step in the MSCS. Under the constraint of f_max determined by Eq. (2) and the requirement of MRTD scheme, the discrete grid size in time and space domain can be derived as well, which will also be introduced in section 2.3 in detail.

(2) Iterative solution of the near electromagnetic field: having determined the basic parameters of IMRTD model, the pulse plane wave $E_{G P}^{i} (t)$ is introduced into computational domain and the transient near electromagnetic field, $E (t, r)$ and $H (t, r)$ is calculated by MRTD technique. In this process, the incident wave is introduced by 1D-MRTD scheme and space-projection method [31]; in the iterative process, the incident electric field $E_{G P, r_{0}}^{i} (t)$ at the original point r₀ should be recorded at each iterative step.

(3) The extraction and normalization of the time-harmonic components of near field: based on the near field obtained above, then $E (t, r)$ and $E_{G P, r_{0}}^{i} (t)$ are transformed into the frequency domain by FT, and their time-harmonic electric components at the required frequencies ( $E (f, r)$ and $E_{G P, r_{0}}^{i} (f)$ , λ = c/f, the frequencies should be set according to particle sizes needed to be simulated) can be easily extracted from the obtained results; then, we further normalize the electric field $E (f, r)$ by $E_{G P, r_{0}}^{i} (f)$ , and the scattering transfer equation of the scatterer $G (f, r)$ can be calculated, i.e., $G (f, r) = E (f, r) / E_{G P, r_{0}}^{i} (f)$ . It can be found that $G (f, r)$ is independent of the incident wave, and it is only decided by particle’s size, shape and refractive index at λ₀. Since the near electric field $E (t, r)$ is discrete in space and time domain, its Fourier transform should be carried out by the following way:

E (f, r) = \int_{- \infty}^{+ \infty} E (t, r) \exp (- j 2 π c t / λ) d t = Δ t \sum_{n = 0}^{N} E (n Δ t, r) \exp (- j 2 π c n Δ t / λ),

where, t denotes the time,

Δ t

is time interval,

r

is space position vector, N is total iterative step. To save computer memory, the FT of

E (t, r)

can be implemented with iteration process.

(4) The calculation of the scattering parameters and results correction: based on the scattering transfer equation $G (f, r)$ , the far electric field at different frequencies can be calculated by volume integral principle, and the scattering properties, such as absorption cross section C_ab(r₀, λ), extinction cross section C_ext(r₀, λ) and scattering phase matrix F(r₀, λ) at different wavelengths can be computed by the models proposed in literature [26] (where r₀ is the original radius, and λ is the wavelength).

To obtain the scattering parameters of particles with different sizes, the results at different wavelengths should be further normalized to λ₀ (the simulation wavelength). Though the scattering properties of particles with the same size parameters x and refractive indices are similar, another scalar correction is required when the scattering parameters at a specific wavelength are transformed into those of another, and the scalar ratios for the correction are closely related to particle size and light wavelength. For particle with a radius of r (r satisfies r₀/λ₀ = r/λ = r/(c∙f), where r₀ is the original particle radius), the corrected scattering phase matrix F(r), absorption cross section C_ab(r) and extinction cross section C_ext(r) at λ₀ can be calculated by

F (r) = \frac{λ_{0}^{2}}{λ^{2}} F (r_{0}, λ);

C_{a b} (r) = \frac{r^{2} C_{a b} (r_{0}, λ)}{r_{0}^{2}} = \frac{r^{2} Z_{0} ∭_{V} σ (r) G (f, r) \cdot G^{*} (f, r) d^{3} r}{r_{0}^{2} {| E_{i n c} (f, r) |}^{2}};

C_{e x t} (r) = \frac{r^{2} C_{e x t} (r_{0}, λ)}{r_{0}^{2}} = - \frac{Z_{0} r^{2}}{r_{0}^{2} {| E_{i n c} (f, r) |}^{2}} {∭_{V} Im [ω (ε - ε_{0} - j \frac{σ}{ω}) E_{i n c} (f, r) \cdot G_{}^{*} (f, r)] d^{3} r};

where, Z₀ is the wave impedance,

σ

and

ε

are the permittivity and electric conductivity of the particle, ε₀ is the permittivity of vacuum, and

ω = 2 π f

, is the angular frequency.

Compared with traditional MRTD scattering model, IMRTD improved by the MSCS can reduce the computational time dramatically. In traditional MRTD model, the time-harmonic wave is generally taken as the incident wave, so in order to guarantee the simulation accuracy, the iteration process can’t be stopped until the electromagnetic field becomes stable. In this case, the time needed should be at least longer than the time for 3~5 round-trips along the diagonals of the computational space with the velocity of light. But for IMRTD, since the interaction between the particle and pulse plane wave is transient, the field stability is not required any more, therefore the time needed can be reduced by nearly 50% for the same space resolution.

2.2 Implementation of MRTD scattering model

MRTD scattering model is a rigorous method for simulating light scattering by arbitrarily shaped aerosol particles. As demonstrated in Fig. 2, there are three core modules in the MRTD model, i.e., the near field computational module, near-to-far transformation module and scattering parameters calculation module. The near field computational module is designed to calculate the near electromagnetic field, in this module, the pulse plane wave is introduced into the computational space by TF/SF technique, and the Convolution Perfectly Matched Layer (CPML) is applied to truncate the computational domain. In the near-to-far transformation module, time-harmonic field components at different frequencies are extracted by FT, and the volume integral method is employed to transform the near electric field to the far field. After the far electric is obtained, the scattering parameters calculation module is applied to calculate the scattering properties of aerosol particles. In order to obtain the phase matrix, the simulation process should be implemented for the horizontally and vertically polarized incident wave independently.

In MRTD scheme, the Maxwell’s equations are discretized by Galerkin principle rather than the central differential method. Taking E_x for example, the field component is firstly expanded by the Haar scaling function $h_{n} (t)$ for time domain and Daubechies scaling function $ϕ_{i} (x)$ for space, written as

E_{x} (r, t) = \sum_{i, j, k, n = - \infty}^{+ \infty} E_{i + 1 / 2, j, k}^{ϕ x, n} ϕ_{i + 1 /2} (x) ϕ_{j} (y) ϕ_{k} (z) h_{n} (t) .

Substituting Eq. (7) into the Maxwell’s equations, and using Galerkin scheme for simplification, the iterative equation of E_x can be derived as [26]

E_{i + 1 / 2, j, k}^{ϕ x, n + 1} = \frac{2 ε - σ Δ t}{2 ε + σ Δ t} E_{i + 1 / 2, j, k}^{ϕ x, n} + \frac{2 Δ t}{2 ε + σ Δ t} \sum_{l = - L s}^{L s - 1} a (l) (H_{i + 0.5, j + l + 0.5, k}^{ϕ z, n + 1 / 2} / Δ y - H_{i + 0.5, j, k + l + 0.5}^{ϕ y, n + 1 / 2} / Δ z),

where, m = (i + 1/2, j, k), denotes the discrete coordinate of E_x, ∆x, ∆y and ∆z are the discrete interval;

a (l)

is called connection coefficient, where the variation range of l is -Ls~Ls-1, Ls = 2N-1 is the effective support size of scaling functions (N is vanishing wavelet moment).

a (l)

can be calculated by the inner products of scalar functions numerically in the Fourier domain, expressed as

a (l) \equiv 〈 \frac{d ϕ_{j + 1 / 2} (x)}{d x}, ϕ_{j - l} (x) 〉 = {\int_{0}^{\infty} ω | \hat{ϕ} (ω) |}^{2} \sin (ω (l + 1 / 2)) d ω,

where,

| \hat{ϕ} (ω) |

is the Fourier transform of

ϕ

.

2.3 Design of discrete grid and incident wave

The grid mesh of IMRTD scattering model is similar to the lowpass filter in signal processing field, and once the grid size is given, the upper cut-off frequency f_max of “the lowpass filter” is determined as well. If the aerosol particle is illuminated by the light with a wide spectrum, then the time-harmonic field components with their frequencies higher than f_max will be filtered in the iteration process and the simulation results at these frequencies will deviate from the real scattering properties of aerosol particles. Owning to this reason, the design of the spatial grid mesh and incident wave is extremely important in the implementation of the MSCS.

The method to determine the computational grid and incident pulse wave is presented in Fig. 3. This process can be roughly divided into three steps: at first, the pulse width and upper-limit frequency of incident light are determined according to the type of the pulse wave and the upper cut-off frequency f_max; then, under the constraint of the requirements of IMRTD scheme and discrete accuracy of the pulse wave, the discrete time interval is designed; at last, the method for the determination of spatial grid size can be derived based on the Courant stability condition and discrete time interval. To illustrate this process in detail, here we take the Gauss pulse wave as the example to introduce its implementation.

Fig. 3 The steps to design the spatial grid size and discrete time interval.

Download Full Size | PDF

(1) The determination of pulse width and upper-limit frequency of incident light. In the time domain, the Gauss pulse can be expressed as

E (t) = \exp [- 4 π {(t - t_{0})}^{2} / τ^{2}] .

Its expression in frequency domain can be derived by using Fourier transform, written as

E (f) = τ / 2 \cdot \exp [- j 2 π f t_{0} - π f^{2} τ^{2} / 4],

in which, t denotes the propagation time, t₀ is the time of the pulse peak (called ‘central time’ in this paper),

τ

is the pulse width. From Eq. (11), it can be found that the spectral amplitude of Gauss pulse decreases with the increasing of the frequency, and its maximum amplitude is located at f = 0Hz, therefore, when designing the pulse wave, what we should pay attention to is its relative amplitude at high frequencies. Assuming that the upper-limit frequency of the incident wave is f_c (f_c should be equal or larger than the upper cut-off frequency f_max (see Eq. (2)), so that all the frequencies we require can be included in [0, f_c]), and the relative spectral amplitude corresponding to f_c is k_c, then the following equation can be established easily, written as

| E (f = f_{c}) / E (f = 0) | = k_{c},

in which, k_c is called the cut-off relative amplitude, and in order to assure that the intensity of the scattering signal is large enough at high frequencies, k_c is generally set as a value larger than 10%. Once k_c is determined, the pulse width

τ_{c}

can be obtained by solving Eq. (11) and Eq. (12), presented as

τ_{c} = \frac{2}{\sqrt{π} f_{c}} \sqrt{\ln (1 / k_{c})}, f_{c} = \frac{2}{\sqrt{π} τ_{c}} \sqrt{\ln (1 / k_{c})} .

From Eq. (13), it can be found that, with the decreasing of the pulse width, the upper-limit frequency is becoming higher, thus the light scattering process can be simulated at a much shorter wavelength, and the scattering parameters of particles with larger sizes can be obtained.

(2) The design of discrete time interval. The design of the discrete time interval ∆t should be implemented under the constraint of the requirements of MRTD scheme, such as the Courant stability condition, numerical dispersion performance and the requirement of the anisotropy errors of phase velocity caused by discretization.

Supposing that the minimum wavelength at which light scattering process should be effectively simulated by IMRTD model is $λ_{\min} = c / f_{c}$ , then in order to guarantee the performance of IMRTD model, the spatial grid size δ should satisfy the following inequality according to computational electromagnetics theory, expressed as

δ \leq \min {λ_{\min} / N_{D}, λ_{\min} / N_{A}},

where, N_D is the minimum grids per wavelength required to guarantee the numerical dispersion performance of IMRTD; N_A is the minimum grids per wavelength to assure that the simulation errors caused by the anisotropy of the phase velocity are small enough for scattering simulation.

To guarantee the stability of the iterative process, the discrete time interval ∆t should satisfies the Courant stability condition, given as

Δ t \leq \frac{1}{c \sum_{l = 0}^{L s - 1} | a (l) | \sqrt{\frac{1}{{(Δ x)}^{2}} + \frac{1}{{(Δ y)}^{2}} + \frac{1}{{(Δ z)}^{2}}}},

If the computational domain is discretized by cubic cells, i.e.,

Δ x = Δ y = Δ z = δ

, then substitute Eq. (14) into Eq. (15), and the requirement for the discrete time interval

Δ t

in the MSCS can be derived as

Δ t \leq \frac{λ_{\min}}{\sqrt{3} c N \sum_{l = 0}^{L s - 1} | a (l) |} \leq \frac{\sqrt{π} τ_{c}}{2 \sqrt{3} N \sum_{l = 0}^{L s - 1} | a (l) | \sqrt{\ln (1 / k_{c})}} = \frac{τ_{c}}{N_{τ}},

where, N is the minimum grids number in the wavelength λ_min, and it should satisfy the inequality,

N \geq \max {N_{D}, N_{A}}

;

N_{τ}

denotes the number of the sampling points in the pulse width

τ_{c}

. By using Eq. (16), the discrete time interval ∆t can be designed reasonably.

(3) The determination of the spatial grid size. Based on the discrete time interval we have obtained, the method to calculate the spatial grid size δ can be derived under the constraint of Courant stability condition, written as

δ \leq \sqrt{3} c \sum_{l = 0}^{L s - 1} | a (l) | Δ t \leq \frac{\sqrt{3} c \sum_{l = 0}^{L s - 1} | a (l) | τ_{c}}{N_{τ}} .

If the grid size δ is smaller enough to construct the scatter’s shape accurately, then δ will be set as the spatial resolution for IMRTD calculation.

3. Model simulation, validation and discussion

The validation of the IMRTD scattering model is carried out in three aspects: (1) first, the precision of our model is tested by comparing its simulation results with those of Mie and T-Matrix theory for spherical and cylindrical particles. (2) Considering the importance of discrete spatial resolution for IMRTD, the influence of spatial grid size on the simulation accuracy is discussed. (3) At last, the computational efficiency of the IMRTD model is analyzed and compared with traditional MRTD model.

3.1 Validation of the precision of IMRTD Model

3.1.1 Spherical particle case

In this test, the simulation wavelength λ₀ is 0.5μm, the original radius of the particle is set as r₀ = 0.5μm, and its refractive index is taken to be m = 1.532-0.008i; according to the method proposed in section 2.3, the spatial resolution for space discretization is set as δ = λ₀/60, and the discrete time interval is taken to be ∆t = δ/3c; the Gauss pulse light is introduced as the incident wave, where its pulse width is set as τ = 60∆t (i.e., τ = 5.56 × 10⁻⁸ns), and in order to assure that the intensity of the incident wave at t = 0 equals to 0 approximately, the central time t₀ is set as 0.8τ. To validate the simulation accuracy of our model, the 4 × 4 scattering phase matrix F and integral scattering properties (i.e., extinction, scattering and absorption cross sections) of five particles with different sizes are calculated by IMRTD in one computation, and the results are compared with those of Mie scattering theory and traditional MRTD model (for the integral scattering parameters), where the radii of these particles are 0.25μm, 0.3571μm, 0.5μm, 0.833μm and 1.0μm, respectively.

The validation of scattering phase matrix. The scattering phase matrices simulated by different models are presented in Fig. 4 and Fig. 5. In each figure, the variation of the phase matrix elements with scattering angle is shown in the left panel, and in the right panel, their computational errors are demonstrated, where for element F₁₁ (i.e., the scattering phase function), its simulation errors are evaluated by relative errors; for F₁₂, it is normalized by F₁₁, and its deviations are described by absolute errors.

Fig. 4 The comparison between the scattering phase functions simulated by the IMRTD model and Mie theory.

Download Full Size | PDF

Fig. 5 The comparison between the scattering phase matrix element F₁₂ simulated by the IMRTD scattering model and Mie theory.

Download Full Size | PDF

As demonstrated in the figures, an excellent agreement is achieved between the results of different scattering models, where for the particles with r = 0.25μm and 0.5μm, the maximum errors of F₁₁ are 6.4% and 11.3%, respectively, indicating that our model can effectively simulate the scattering properties of particles with different sizes in one computation. It can also be found that the simulation errors of F₁₁ and F₁₂/F₁₁ increase with the increasing of particle radius, and the simulation precision for small particles is notably higher than that for large ones. This phenomenon can be explained by the principle of MSCS. In our computational model, the scattering properties of large particles are normalized from those of a size-fixed particle at short wavelengths, namely, the electromagnetic field for their scattering calculation is the high-frequency component. According to the spectral distribution of Gauss pulse wave, the relative spectral amplitude of the signal is exponentially attenuated with the increasing of frequency, so the intensity of the signal is relatively weak at short wavelengths; therefore the Signal-Noise Ratio (SNR) of the electromagnetic field extracted at these wavelengths should be lower than that at longer ones correspondingly, which further results in the decreasing of their scattering simulation accuracy. From the angular distribution of the simulation errors, it can be found that the simulation errors in forward directions are almost as large as those at large scattering angles, where for particle with r = 0.8333μm, the maximum errors of F₁₁ in forward directions and backward scattering angles are 18.31% and 18.53%, respectively, which is different from the angular distribution of traditional MRTD model (for traditional MRTD model, the simulation errors in forward directions are much smaller than those at large scattering angles, readers can refer to ref [26] for detail).

The validation of integral scattering parameters. The integral scattering properties obtained by IMRTD, traditional MRTD and Mie scattering theory are also compared with each other, where, the spatial resolution for the traditional MRTD is set as δ = λ₀/40, and the results are shown in Table 1. As can be seen, the integral scattering parameters calculated by IMRTD shown a great consistency with those simulated by Mie theory, where for the particles with r≤0.8333μm, the relative simulation errors of the extinction cross section C_ext , absorption cross section C_ab and scattering cross section C_sc are less than 4% for most cases, validating the accuracy of our model. Compared with the results of the traditional MRTD, the precision of the integral scattering properties simulated by IMRTD is relatively lower; the phenomenon can be explained by the following reasons. In IMRTD model, both the incident and the scattering field are transient pulse signals, and the number of the sampling points in such a short time is limit (usually 40~80 points per pulse width τ), therefore, it is unavoidable that computational errors will be produced in the extraction of the harmonic components from the discretized pulse wave. Besides, the incident light in traditional MRTD model is monochromatic, and the energy are concentrated in single frequency, however for the IMRTD, signal’s energy is distributed in a wide frequency spectrum, and the energy at specific frequencies is very small, thus the SNR of the computational results of IMRTD model is relatively lower than that of the monochromatic light case, which further results in the decreasing of simulation accuracy.

Table 1. Comparison of the integral scattering properties simulated by the IMRTD, traditional MRTD and Mie scattering theory (Unit: μm²).

View Table | View all tables in this article

3.1.2 Cylindrical particle case

The precision of the IMRTD model is further validated against T-Matrix method for cylindrical particles. In this simulation, the diameter D and length L of the original cylindrical particle are all set as 1.0μm, the refractive index m, spatial resolution δ and simulation wavelength λ₀ are the same as those for spherical particle in section 3.1.1. The incident light is taken as the Gauss pulse plane wave, where its pulse width τ is taken as 60∆t and the central time t₀ is set as 0.8τ, its propagation direction is perpendicular to the bottom surface of the cylinder. The scattering properties of five particles with different sizes are simulated by IMRTD and T-Matrix method, and their differences are quantitatively given, where the equivalent radii of the particles are 0.2862μm, 0.4088μm, 0.5724μm, 0.9539μm and 1.1447μm, respectively (their equivalent radii are obtained by equal-volume principle).

The validation of scattering phase matrix. The scattering phase matrix elements calculated by two models are presented in Fig. 6 and Fig. 7, the patterns of the figures are the same as Fig. 4 and Fig. 5. From the figures, it can be found that the curves of the matrix elements, F₁₁ and F₁₂/F₁₁, obtained by different models are approximately coincided with each other, where for the particle with r_eq = 0.4088μm, the maximum simulation error of F₁₁ and F₁₂/F₁₁ are 6.84% and 0.0587, indicating that IMRTD can effectively simulate the scattering properties of non-spherical particles with different sizes in one calculation. Similar to the results of spherical particle, the simulation accuracy of the scattering phase matrix is reduced with the increasing of particle size, where for the particle with r_eq = 1.1447μm, the maximum simulation errors of F₁₁ and F₁₂/F₁₁ reach 91.25% and 0.29, respectively, which is much larger than those of the particle with r_eq = 0.2862μm, the phenomenon can be explained by the low signal intensity at high frequencies. On the whole, the angular distributions of the simulation errors of phase matrix elements are in accordance to their variation characteristics with scattering angle, i.e., if F₁₁ and F₁₂/F₁₁ fluctuate with the scattering angle dramatically, then their deviations will be notably larger than that at other directions.

Fig. 6 The comparison of the scattering phase functions simulated by the IMRTD model and T-matrix method.

Download Full Size | PDF

Fig. 7 The comparison of the scattering phase matrix element F₁₂ simulated by the IMRTD model and T-Matrix method.

Download Full Size | PDF

The validation of integral scattering parameters. The extinction, absorption and scattering cross sections of these particles obtained by IMRTD, traditional MRTD model and T-Matrix method are listed in Table 2. As can be seen, the relative simulation errors of the integral scattering parameters calculated by IMRTD are generally smaller than 4%, validating the effectiveness and precision of our model. it can also be found that the IMRTD achieves the best simulation accuracy for the particles with their equivalent radius comparable to the simulation wavelength λ₀, and the simulation errors increases with particle size becoming larger or smaller, where for particle with r_eq = 0.5724μm, the relative errors of C_ext, C_ab and C_sc are only 3.2112%, 3.5619% and 3.1552%, respectively, while for particle with r_eq = 1.1447μm, their simulation errors increase to 4.1063%, 5.7728% and 6.6897%. Similar to the phenomenon observed in section 3.1.1, for non-spherical particle case, the simulation accuracy of IMRTD model is also relatively lower than that of traditional MRTD scattering model, especially for particles with large particle size.

Table 2. Comparison of the integral scattering parameters simulated by the IMRTD, traditional MRTD and T-Matrix method (Unit: μm²).

View Table | View all tables in this article

3.2 Influence of grid size on simulation accuracy

Spatial grid size is an important factor influencing the performance of the IMRTD model since it can affect both the upper cut-off frequency of the model and the construction accuracy of particle shape. To this end, the influence of grid size on calculation accuracy is investigated in this section. To quantitatively evaluate the simulation results, the precision of the integral scattering parameters are evaluated by relative errors, and the simulation accuracy of the angular distribution of scattering light is estimated by mean square root error of the scattering phase function F₁₁, defined as

R S E F = {[\frac{1}{N} \sum_{i = 1}^{N} {(\frac{F_{11} {(θ_{i})}_{I M R T D} - F_{11} {(θ_{i})}_{M i e}}{F_{11} {(θ_{i})}_{M i e}})}^{2}]}^{1 / 2},

Light scattering processes are performed by IMRTD model for different spatial resolutions. In this test, the simulation is implemented for spherical particle, and its original radius, refractive index are the same as that in section 3.1.1, the spatial grid sizes δ are set as λ₀/40, λ₀/50, λ₀/60 and λ₀/80, respectively. Similar to the simulations above, the Gauss pulse is taken as the incident wave, where, its pulse width is τ = 20δ/c, and the central time is set as t₀ = 0.8τ. The simulation errors of particles with different radius are obtained by comparing the results of IMRTD with those of Mie theory, and the results are shown in Table 3.

Table 3. The influence of space resolution on the simulation accuracy of the scattering properties of aerosol particles with different sizes.

View Table | View all tables in this article

As can be seen, with the increasing of spatial resolution, the simulation accuracy of the scattering parameters is improved for all particles, however the improvement is gradually reduced as the spatial resolution becomes higher. From the table, it can also be found that the precision improvement for large particles is much more remarkable than that for small ones, where for particle with r = 1.5μm, RSEF decreases by a percentage of 13.375% with δ decreasing from λ₀/40 to λ₀/80, while RSEF is reduced only by 2.891% for particle with r = 0.3μm. The phenomenon above can be explained from two aspects: on one hand, as δ becomes smaller, the precision of particle’s shape becomes higher, thus the simulation errors caused by stair-effect are reduced. On the other hand, since the pulse width is determined by τ = 20δ/c, both the upper cut-off frequency of IMRTD model and the relative amplitude of transient field at high frequencies will be improved with the decreasing of δ (see Eq. (13)), thus the accuracy of the light scattering simulation for large particles can be improved notably. While for the small particles, its scattering properties are obtained from the scattering signal at low frequencies, and the relative spectral amplitude is large enough at these wavelengths, thus the improvement of its simulation precision is mainly resulted from the better-construction of particle shape, so the increment of its calculation accuracy is not as notable as that for large particles.

3.3 Evaluation of the model’s computational efficiency

The evaluation of model’s efficiency is performed by the following way: (1) 10 spherical particles with different sizes are simultaneously simulated by IMRTD model, and its computational time T₁ is recorded, where the radius of 10 particles are uniformly distributed in 0~r_max (r_max is maximum particle radius that can be effectively simulated by IMRTD, where the simulation accuracy for particles with a radius of r_max should satisfy RSEF<15%), the refractive index is set as 1.42-0.05i, the simulation wavelength λ₀ is taken as 0.5μm, and the original radius r₀ and spatial grid sizes are listed in Table 4. (2) Calculate the scattering properties of these particles by traditional MRTD model, and record the total time consumed, T₀, in this simulation, the spatial resolution for particles with r = 0~0.3μm, r = 0.3~0.8μm and r>0.8μm are set as λ₀/60, λ₀/48, and λ₀/36, respectively. (3) Then, the computational time of different models are compared with each other for different original radius r₀.

Table 4. The comparison between the computational efficiency of improved and traditional MRTD scattering model.

View Table | View all tables in this article

The computational processes are performed on the computer with 2.3GHz CPU and 8G memories, and the results are shown in Table 4. It can be found that the computational efficiency of IMRTD is much higher than that of traditional MRTD, where in the case that r₀ = 0.5μm, the computational time consumed by IMRTD is approximate 1/6 of the time needed by the traditional method. With the increasing of r₀, the maximum simulation radius r_max increase as well, but it should be noted that the ratio r_max/ r₀ is reduced to some extent as r₀ increases. This can be explained from two aspects. On one hand, similar to traditional MRTD model, the simulation precision will be reduced as the increasing of particles size; on the other hand, the SNR of the electromagnetic field simulated at the high frequencies is relatively low, which in turn aggravates the deterioration of the calculation accuracy of large particles.

The scattering properties needed for atmospheric radiative transfer should be integrated over a wide size range and different orientations, therefore two issues should be draw attention to: the first issue is the contradiction between the pulse width and maximum particle radius r_max, r_max will increase with the decreasing of τ, but smaller τ means high time and space resolution, so for large particles, its computational amount would be large. The second issue is that the number of orientations needed for orientation averaging process increases with respect to the size parameter, if the same number of orientations are performed for all sizes, it will not be efficient. To overcome these difficulties, the whole particle size range should be divided into several subsections, and the scattering simulation should be implemented in these subsections independently. In this scheme, the original radius r₀ is taken as the central value of the subsection, thus the pulse width τ needn’t to be too small, and the spatial resolution can be degraded as well (but the discrete time interval should be reduced to some extent to guarantee that there are enough sampling points in the pulse width). To assure the accuracy of orientation averaging, more orientations should be simulated for the subsections with large sizes; while for the subsections with small particles, the number of the orientations can be reduced. The determination of size subsections should be based on the size range and the performance of IMRTD. From the results above, 4~5 time-domain simulations are required to cover the size parameter range 0.1~50.

4. Conclusion

MRTD scattering model is a rigorous tool to simulate light scattering by nonspherical aerosol particles. In order to improve the efficiency of the MRTD model, MSCS is proposed and applied to light scattering simulation. By using this method, the scattering properties of particles with different sizes can be computed simultaneously in one wave-particle interaction simulation. From the analysis and discussion of the results above, several conclusions are drawn as follows:

The IMRTD model can simulate the light scattering process of particles with different sizes simultaneously and accurately, where for spherical particles with a radius of r = 0.5μm, the maximum relative errors of F₁₁ and F₁₂/F₁₁ are only 11.3% and 0.045; on the whole, the simulation precision of small particles are higher than that of large ones, and the simulation errors for the integral parameters, such as extinction, absorption cross sections, are smaller than that of phase matrix.

The grid size has a notable impact on the simulation accuracy. With the increasing of spatial resolution, the simulation precision of the scattering parameters is improved as well, and the improvement for large particles is much more remarkable than that for small ones. The efficiency of IMRTD is much higher than that of traditional one; for light scattering simulation of 10 particles (the original radius of particle is r₀ = 0.5μm), the computational time consumed by IMRTD is approximate 1/6 of the time needed by traditional one.

It should also be noted that in the MSCS, we only simultaneously simulate the scattering properties of particles with different sizes at one wavelength, thus the frequency dependence of particle’s refractive index is not needed to be considered in the computational process. “The scattering properties at different wavelengths” stated in the paper are just the transition results for the calculation (which isn’t the same as the actual scattering properties at different wavelengths), and they will be transformed into the scattering properties at the simulation wavelength finally. Since the near field components at different frequencies for near-to-far transformation should be stored in the memory, the memory consumption for IMRTD will much larger than that for traditional MRTD method.

Funding

Natural Science Foundation of China (NSFC) (41575025, 41575024, 41675003).

Acknowledgments

Thank for the help provided by Professor Chen Bin and Dr. Liu Yawen of National Key Laboratory on Electromagnetic Environment and Electro-Optical Engineering in the study of electromagnetic theory.

References and links

1. K. N. Liou, An Introduction to Atmospheric Radiation (Academic Press, 2003).

2. F. Zhang, K. Wu, J. Li, Q. Yang, J.-Q. Zhao, and J. Li, “Analytical Infrared Delta-Four-Stream Adding Method from Invariance Principle,” J. Atmos. Sci. 73(10), 4171–4188 (2016). [CrossRef]

3. S. Hu, T. Gao, H. Li, L. Liu, X.-C. Liu, T. Zhang, T.-J. Cheng, W.-T. Li, Z.-H. Dai, and X.-J. Su, “Effect of atmospheric refraction on radiative transfer in visible and near-infrared band: Model development, validation, and applications,” J. Geophys. Res. Atmos. 121(5), 2349–2368 (2016). [CrossRef]

4. J.-Q. Zhao, G. Shi, H. Chen, and G. Cheng, “Approximation of scattering phase function of particles,” adv. Atmos. Sci. 23(5), 802–808 (2006). [CrossRef]

5. J.-Q. Zhao and Y.-Q. Hu, “Bridging technique for calculating the extinction efficiency of arbitrary shaped particles,” Appl. Opt. 42(24), 4937–4945 (2003). [CrossRef] [PubMed]

6. J. Li, C. Liu, Y. Yin, and K. R. Kumar, “Numerical investigation on the Ångström exponent of black carbon aerosol,” J. Geophys. Res. Atmos. 121(7), 3506–3518 (2016). [CrossRef]

7. C. Liu and Y. Yin, “Inherent optical properties of pollen particles: A case study for the morning glory pollen,” Opt. Express 24(2), A104–A113 (2016). [CrossRef] [PubMed]

8. O. Dubovik, A. Sinyuk, T. Lapyonok, B. N. Holben, M. Mishchenko, P. Yang, T. F. Eck, H. Volten, O. Munoz, B. Veihelmann, W. J. van der Zande, J. F. Leon, M. Sorokin, and I. Slutsker, “Application of spheroid models to account for aerosol particle nonsphericity in remote sensing of desert dust,” J. Geophys. Res. 111(D11), D11208 (2006). [CrossRef]

9. M. Herman, J. L. Deuzé, A. Marchand, B. Roger, and P. Lallart, “Aerosol remote sensing from POLDER//ADEOS over the ocean:Improved retrieval using a nonspherical particle model,” J. Geophys. Res. 110(D10), D10S02 (2005). [CrossRef]

10. P. Yang, K.-N. Liou, L. Bi, C. Liu, B. Yi, and B. A. Baum, “On the Radiative Properties of Ice Clouds: Light Scattering, Remote Sensing, and Radiation Parameterization,” Adv. Atmos. Sci. 32(1), 32–63 (2015). [CrossRef]

11. T. H. Cheng, X. F. Gu, T. Yu, and G. L. Tian, “The reflection and polarization properties of non-spherical aerosol particles,” J. Quant. Spectrosc. Radiat. Transf. 111(6), 895–906 (2010). [CrossRef]

12. M. I. Mishchenko and L. D. Travis, “Capabilities and Limitations of a Current Fortran Implementation of the T-Martrix Method for Randomly Oriented,Rotationally Symmetric Scatterers,” J. Quant. Spectrosc. Radiat. Transf. 60(3), 309–324 (1998). [CrossRef]

13. M. I. Mishchenko and L. D. Travis, “T-Martrix computations of light scattering by large spheriodal particles,” Opt. Commun. 109(1-2), 16–21 (1994). [CrossRef]

14. L. Bi and P. Yang, “Accurate simulation of the optical properties of atmospheric ice crystals with the invariant imbedding T-matrix method,” J. Quant. Spectrosc. Radiat. Transf. 138, 17–35 (2014). [CrossRef]

15. L. Bi and P. Yang, “Impact of calcification state on the inherent optical properties of Emiliania huxleyi coccoliths and coccolithophores,” J. Quant. Spectrosc. Radiat. Transf. 155, 10–21 (2015). [CrossRef]

16. N. V. Voshchinnikov and V. G. Farafonov, “Optical properties of spheroidal particles,” Astrophys. Space Sci. 204(1), 19–86 (1993). [CrossRef]

17. H. M. Al-Rizzo and J. M. Tranquilla, “Electromagnetic scattering from dielectrically coated axisymmetric objects using the generalized point-matching technique,” J. Comput. Phys. 119(2), 356–373 (1995). [CrossRef]

18. R. F. Harrington, Field Computation by Moment Methods (Macmillan, 1968).

19. B. T. Draine and P. J. Flatau, “Discrete dipole approximation for scattering calculations,” J. Opt. Soc. Am. A 11(4), 1491–1499 (1994). [CrossRef]

20. P. Yang and K. N. Liou, “Finite-difference time domain method for light scattering by small ice crystals in three-dimensional space,” J. Opt. Soc. Am. A 13(10), 2073–2085 (1996). [CrossRef]

21. M. A. Morgan and K. K. Mei, “Finite-element computation of scattering by inhomogeneous penetrable bodies of revolution,” IEEE Trans. Antenn. Propag. 27(2), 202–214 (1979). [CrossRef]

22. D.-B. Ge and Y.-B. Yan, Finite-Difference Time-Domain method for Electromagnetic Waves (Xidian University, 2011).

23. C. Liu, R. L. Panetta, and P. Yang, “Application of the pseudo-spectral time domain method to compute particle single-scattering properties for size parameters up to 200,” J. Quant. Spectrosc. Radiat. Transf. 113(13), 1728–1740 (2012). [CrossRef]

24. C. Liu, R. L. Panetta, and P. Yang, “The effects of surface roughness on the scattering properties of hexagonal columns with sizes from the Rayleigh to the geometric optics regimes,” J. Quant. Spectrosc. Radiat. Transf. 129, 169–185 (2013). [CrossRef]

25. G. Chen, P. Yang, and G. W. Kattawar, “Application of the pseudospectral time-domain method to the scattering of light by nonspherical particles,” J. Opt. Soc. Am. A 25(3), 785–790 (2008). [CrossRef] [PubMed]

26. S. Hu, T. Gao, H. Li, B. Yang, F. Zhang, M. Chen, and L. Liu, “Light scattering computation model for nonspherical aerosol particles based on multi-resolution time-domain scheme: model development and validation,” Opt. Express 25(2), 1463–1486 (2017). [CrossRef] [PubMed]

27. S. Hu, T.-C. Gao, H. Li, B. Yang, M. Chen, L. Liu, and G. Li, “Design and validation of the Parallelized scattering model for nonspherical aerosol based on MRTD method,” Acta Opt. Sin. (accepted).

28. M. I. Mishchenko, J. W. Hovenier, and L. D. Travis, light scattering by nonspherical particles, Thoery, Measurements, and Application (Academic Press, 2000).

29. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (John Wiley & Sons, Inc, 1983).

30. H. C. van der Hulst, light scattering by small particles (Dover Publications, Inc, 1981).

31. Q. Gao, Q. Cao, and J. Zhou, “Application of Total-Field/Scattered-Field Technique to 3D-MRTD Scattering Scheme,” in International Forum on Information Technology and Applications (2009), pp. 359–362.

r /μm	C_ext			C_ab			C_sc
r /μm	Mie	MRTD	IMRTD	Mie	MRTD	IMRTD	Mie	MRTD	IMRTD
0.25	0.7156	0.7003	0.6925	0.0235	0.0229	0.0227	0.6921	0.6774	0.6698
0.3571	1.2313	1.2095	1.1746	0.0473	0.0482	0.0460	1.1838	1.1613	1.1286
0.5	1.6607	1.6803	1.7027	0.1941	0.1924	0.1934	1.4666	1.4879	1.5093
0.8333	6.4150	6.5021	6.5629	0.8619	0.8801	0.9030	5.5531	5.6220	5.6599
1.0	6.3776	6.4017	6.7685	1.1938	1.1738	1.2478	5.1838	5.2279	5.5207

r_eq /μm	C_ext			C_ab			C_sc
r_eq /μm	T-Matrix	MRTD	IMRTD	T-Matrix	MRTD	IMRTD	T-Matrix	MRTD	IMRTD
0.2862	1.0633	1.0387	1.0236	0.0365	0.0354	0.0382	1.0268	1.0033	0.9854
0.4088	1.8153	1.7904	1.7540	0.1093	0.1078	0.0972	1.7060	1.6826	1.6568
0.5724	2.1830	2.1552	2.1129	0.3004	0.2965	0.2897	1.8826	1.8587	1.8232
0.9539	7.4621	7.4223	7.2102	1.2606	1.2494	1.2046	6.2012	6.1729	6.0056
1.1447	9.6315	9.5761	9.2360	1.9973	2.0403	2.1126	7.6341	7.5358	7.1234

Grid	r	0.2μm	0.3μm	0.4μm	0.5μm	0.6μm	0.8μm	1.0μm	1.5μm
λ₀/40	C_ext (%)	3.8753	3.8405	3.6108	3.2795	4.4602	4.7589	8.9084	18.8753
	C_ab (%)	4.6794	4.5833	4.7324	3.6231	5.2003	5.9733	8.7723	14.7718
	RSEF (%)	5.9084	6.8993	5.7874	6.8875	10.0944	12.0931	19.4301	31.0094
λ₀/50	C_ext (%)	3.7364	3.4635	4.3145	2.7343	3.4957	3.6773	6.9245	13.4834
	C_ab (%)	4.7927	4.4176	4.0428	2.4322	4.0488	4.8931	10.6423	10.7312
	RSEF (%)	5.3273	5.4641	4.9445	6.4744	9.1334	11.1032	16.3335	27.2341
λ₀/60	C_ext (%)	3.4328	2.9235	3.2133	2.5200	3.0134	3.3432	6.1293	9.3289
	C_ab (%)	3.7215	3.8143	4.6214	0.3606	3.0452	3.6323	4.5323	7.7842
	RSEF (%)	4.4230	4.4624	4.8147	5.9639	8.3445	10.1529	12.3425	21.2332
λ₀/80	C_ext (%)	2.6673	2.5262	3.8721	3.8432	3.5487	2.6466	4.0734	7.7634
	C_ab (%)	3.0232	2.9775	4.4498	1.2364	3.6593	3.0873	4.8840	6.5329
	RSEF (%)	4.2635	4.0086	4.6705	5.2322	6.2552	8.6687	10.3425	17.6342

No.	r₀ /μm	Grid	r_max/μm	Computational time of the improved MRTD model /s	Computational time of the traditional MRTD model /s
1	0.25	λ₀/60	0.55	559	4.5303 × 10³
2	0.5	λ₀/60	1.00	3462	2.1542 × 10⁴
3	0.75	λ₀/60	1.40	1.3318 × 10⁴	4.1593 × 10⁴
4	1.0	λ₀/50	1.65	2.3726 × 10⁴	6.8595 × 10⁴

r /μm	C_ext			C_ab			C_sc
r /μm	Mie	MRTD	IMRTD	Mie	MRTD	IMRTD	Mie	MRTD	IMRTD
0.25	0.7156	0.7003	0.6925	0.0235	0.0229	0.0227	0.6921	0.6774	0.6698
0.3571	1.2313	1.2095	1.1746	0.0473	0.0482	0.0460	1.1838	1.1613	1.1286
0.5	1.6607	1.6803	1.7027	0.1941	0.1924	0.1934	1.4666	1.4879	1.5093
0.8333	6.4150	6.5021	6.5629	0.8619	0.8801	0.9030	5.5531	5.6220	5.6599
1.0	6.3776	6.4017	6.7685	1.1938	1.1738	1.2478	5.1838	5.2279	5.5207

Simultaneously simulating the scattering properties of nonspherical aerosol particles with different sizes by the MRTD scattering model

Abstract

1. Introduction

2. Principles and methodology

2.1 Basic principle of the MSCS

2.2 Implementation of MRTD scattering model

2.3 Design of discrete grid and incident wave

3. Model simulation, validation and discussion

3.1 Validation of the precision of IMRTD Model

3.1.1 Spherical particle case

3.1.2 Cylindrical particle case

3.2 Influence of grid size on simulation accuracy

3.3 Evaluation of the model’s computational efficiency

4. Conclusion

Funding

Acknowledgments

References and links

Cited By

Figures (7)

Tables (4)

Equations (18)

Optics Express