Light scattering computation model for nonspherical aerosol particles based on multi-resolution time-domain scheme: model development and validation

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

doi:10.1364/OE.25.001463

1. Introduction

Light scattering process of aerosol particles and ice crystals is an important factor influencing radiative transfer since it can affect the reflection and absorption capabilities of atmosphere [1–3]. So far, the scattering properties of those particles are still with substantial uncertainties because of their complicated geometries and inhomogeneous compositions, thus further affecting the simulation accuracy of radiative transfer [4–6]. According to the IPCC reports [7], though the estimation accuracy of direct radiative force of aerosol is improved, an uncertainty about −0.9~-0.1W/m² still exists. The inadequate understanding of aerosol’s scattering properties also affects atmospheric optical remote sensing since they act as the input parameters of Radiative Transfer Models (RTM) [8–10]. In most of the RTM, the scattering characteristics are mainly calculated by Lorenz-Mie theory (taking aerosol particles as spherical ones), while particle shapes are not adequately taken into account, which will lead to the decreasing of the retrieval accuracy [9]. Therefore owning to its importance, obtaining the scattering properties of nonspherical aerosol has become a hotspot in atmospheric science [11,12].

In order to simulate the scattering properties of nonspherical aerosols, many scattering computational models are established. These models can be classified into two categories: scattering approximation models and exact computation models [11]. Scattering approximation models mainly apply for particles with extreme size parameters, namely, with their radii much larger or smaller than light wavelength. Typical models of this kind include Rayleigh approximation, Rayleigh–Gans–Stevenson approximation [11], bridging method [13,14], anomalous diffraction theory and GOA (Geometric Optics Approximation) [5], where GOA and its improved versions have been widely employed into the scattering simulation of ice crystals [5, 15]. For exact computation models, particle’s scattering properties are obtained by solving the Maxwell’s curl equations and Helmholtz equations either analytically or numerically with given boundary and initial conditions [11]. According to their calculation principles, exact computation models can be roughly divided into three types: field expansion models, models based on volume integral methods, and micro-element models. TMM (T-Matrix Method) [16, 17], SVM (Separation of Variables Method) [18] and PMM (Point Matching Method) [19] can be classified into the first type, where TMM has become a widely tested and explored model for the improvements by Mishchenko and other researchers [11,16,17,20,21]. The characteristics of these models is that the incident field and scattering field are firstly expanded by suitable vector spherical wave functions, and the transformation relationship of their expansion coefficients are established by solving differential equations or integral equations. These models have the advantages of high computational speed and low memory consumption, besides, high calculation precision can be easily achieved as well [11]. However due to the limitation of the boundary conditions for numerical calculation, most of these models can only apply for regular particles, such as ellipsoid, cylinder, etc. besides, the numerical simulation process becomes unstable for particles with large size parameters, refractive index or extreme shapes as well. Models based on volume integral method include Method of Moments (MoM) [22], Discrete Dipole Approximation (DDA) [23,24] and FIEM (Fredholm integral equation Method). In these models, scattering process can be obtained by numerically solving the Green-function solutions of Helmholtz equations. Though DDA and MoM can be applied for particles with any irregular shapes and inhomogeneous compositions, their calculation processes are usually unstable and time-consuming since large matrix equations are needed to be solved; besides, because their simulation precision closely depends on the number of the dipoles, the contradiction between calculation accuracy and memory consumption is evident as well. Owning to these reasons, DDA and MoM are mainly suitable for small particles, and their precision is not as high as that of field expansion models. With the development of computational electromagnetics, such micro-element models as Finite Difference Time Domain (FDTD), Finite Element Method (FEM) and their improved versions are gradually introduced into the atmospheric optics field, where FDTD has been widely used in the simulation of nonspherical aerosol’s scattering properties [25]. The main advantage of these methods is that they can overcome the singularity of volume integral methods, thus they are more suitable for particles with complex shapes and large sizes. However, similar to DDA and MoM, the simulation accuracy is contradicted to the computational efficiency and memory consumption since the higher precision depends on the smaller spatial resolution for calculation. To overcome this problem, Liu and Yang have introduced the Pseudo Spectral Time Domain (PSTD) into the numerical simulation of light scattering process [26,27]. Because this method can achieve high simulation accuracy with less girds per wavelength, it can decrease the computational time and memory dramatically. Now Liu has successfully calculated particle’s single-scattering properties for size parameters up to 200 [26]. The main disadvantage of PSTD is that the incident field can only be introduced by pure scattering field technique, which is time consuming and difficult for particles with large size and complex shapes [28].

Similar to PSTD, MRTD (Multi-Resolution Time-Domain) is a technique widely used in electromagnetic computational field, which is not only a higher order differential scheme of Maxwell’s equations, but also a combination of the superiorities of MoM and FDTD [29–31]. Owning to its excellent numerical dispersion performance and high precision, MRTD is widely applied to the scattering simulation of large targets, and now it is introduced into the light scattering simulation of nonspherical aerosol particles in atmospheric science field. In our MRTD scattering model, ADE-PML (Perfectly Matched Layer with Auxiliary Differential Equation) absorption boundary is improved by exponential difference method and the incident wave is introduced by Total Field-scattering Field (TF-SF) technique rather than pure scattering method. Since most aerosol particles are dielectric and the scattering parameters needed for atmospheric radiative transfer is different from that for electromagnetic computational field, the scattering parameters calculation models are derived from Maxwell’s curl equations directly in a similar way by Yang [32].

This paper is divided into three parts: in the section 2, the basic principles of MRTD scattering model and its implementation are described. In section 3 and 4, the accuracy of our model is validated against the well-tested models, and its parallel computational efficiency is analyzed as well. In the last part, a summary of our study is given in section 5.

2. Implementation of MRTD scattering model

As shown in Fig. 1, the MRTD scattering model is comprised of three modules: near field calculation module, near-to-far field transformation module, and scattering parameters computation module. The near field calculation module is designed to calculate the electromagnetic field in the near space containing the scatterers. There are four key techniques included in this module, i.e., the iterative technique for solving Maxwell’s equations, design of the absorption boundary conditions, and the introduction of incident wave by TF-SF technique. Near-to-far field transformation module consists of the processes as converting the electromagnetic field from time domain to frequency domain, and calculating the far electric field based on volume integral method, which is the basis for the calculation of scattering phase matrix. Aerosol scattering parameters computation module is to calculate the scattering properties such as absorption, scattering extinction cross section, and phase matrix based on the near or far electric field.

Fig. 1 The basic frame of MRTD scattering model for nonspherical aerosol.

Download Full Size | PDF

2.1 Basic principle of MRTD and its Courant stability condition

The propagation process of electromagnetic wave is fully described by Maxwell’s Equations. For space without any electric charge or current, the Maxwell’s curl Equations can be written as

\nabla \times H = ε \frac{\partial E}{\partial t} + σ E, \nabla \times E = - μ \frac{\partial H}{\partial t} - σ_{m} H;

where

E

and

H

refer to the intensity vector of electric and magnetic field; ε and μ are the permittivity and permeability of the medium;

σ

and

σ_{m}

represent the electric conductivity and magnetic conductivity, respectively.

In the field of atmospheric radiative transfer, the optical property of aerosol is described by refractive index m, where its real part and imaginary part denote the ability of scattering and absorption, respectively. Since the parameters needed for Maxwell’s equations are permittivity and electric conductivity, it is necessary to transform m into these electromagnetic parameters. The relationship of m and complex permittivity ε_c can be expressed as

ε_{c} = ε_{0} m^{2} = ε - \frac{j σ}{ω},

where,

ω

is the angular frequency,

ε_{0}

is the permittivity of vacuum, Assuming refractive index has the form as

m = m_{r} - j m_{i}

, then from Eq. (2),

ε

and

σ

can be expressed as

ε = ε_{0} (m_{r}^{2} - m_{i}^{2})

,

σ = 2 ω m_{r} m_{i} ε_{0}

. Base on the two equations above, the electromagnetic parameters of aerosol particles can be easily calculated from the given refractive index.

The basic principle of MRTD is first put forward by Krumpholz and Katehi [33]. Since MRTD is a high order differential form of Maxwell equations, its implementation is more complicated than FDTD, its principle is introduced as follows. First, the electromagnetic field components are expanded by the Haar scaling function for time domain and orthonormal scaling function for space, respectively, shown as follows [34]

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);

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

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

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

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

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

where,

r = (x, y, z)

is position vector of the field components.

E_{i + 1 / 2, j, k}^{ϕ x, n}

,

E_{i, j + 1 / 2, k}^{ϕ y, n}

and

E_{i, j, k + 1 / 2}^{ϕ z, n}

are the coefficients for the electric field expansions in terms of the scaling functions along x, y and z axis.

H_{i, j + 1 / 2, k + 1 / 2}^{ϕ x, n + 1 /2}

,

H_{i, j + 1 / 2, k + 1 / 2}^{ϕ y, n + 1 /2}

and

H_{i + 1 / 2, j + 1 / 2, k}^{ϕ z, n + 1 /2}

denote the expansion coefficients of magnetic field. i, j and k refer to the discrete spatial coordinate of the field components, and n is the iterative step.

h_{n} (t)

is Haar scaling function, defined as

h_{n} (t) = h (t / Δ t - n) .

ϕ_{k} (x)

is the Daubechies orthonormal scaling function, which has the form as

ϕ_{k} (x) = ϕ (x / Δ x - k + M_{1}),

where, M₁ is the first-order moment of the scaling function. In Daubechies’ scaling functions, the shifted interpolation property is almost satisfied, namely,

ϕ_{i} (x) \approx δ_{i, 0}

, where

δ_{i, i^{'}}

is the Dirac delta function.

For the sake of simplicity and generality, we only take E_x for example to introduce the implementation of MRTD scheme. Substituting Eq. (3) into Eq. (1) and using the Galerkin scheme for simplification, the iterative equation for E_x can be expressed as

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

where,

a (l)

is called connection coefficient. Ls = 2N-1, denotes the effective support size of the scaling functions (N is vanishing wavelet moment, which is set as 2 in our model).

a (l)

can be calculated by the inner products of scalar functions numerically in the Fourier domain, written 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 ω,

in which

| \hat{ϕ} (ω) |

is the Fourier transform of

ϕ (x)

. Then simplifying Eq. (11), the iterative equation for electric field component E_x can be obtained, given as

E_{i + 1 / 2, j, k}^{ϕ x, n + 1} = C A (m) E_{i + 1 / 2, j, k}^{ϕ x, n} + C B (m) \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),

in which, m = (i + 1/2, j, k), denotes the discrete coordinate of E_x component, and CA and CB are the iterative coefficients, which can be calculated by

C A (m) = \frac{2 ε - σ Δ t}{2 ε + σ Δ t}, C B (m) = \frac{2 Δ t}{2 ε + σ Δ t} .

The iterative equations for other electromagnetic field components can be obtained by the similar way straightforwardly. Because of the shifted interpolation property of $ϕ_{i} (x)$ , the actual value of the field component at $r_{0}$ equals to its expansion coefficient approximately, i.e., $E_{x} (r_{0}, t_{0}) = E_{i + 1 / 2, j, k}^{ϕ x, n}$ . Therefore the field construction process by integrating over space and time domain can be avoided, which is a great advantage of Daubechies scalar functions over others.

Because MRTD scheme is an explicit difference scheme of Maxwell’s equations, the time discrete interval $Δ t$ and spatial discrete intervals should satisfy the Courant condition to guarantee the stability of the iterative process. For MRTD scheme with Daubechies scalar functions, time discrete interval $Δ t$ should satisfy the following equation [34]

Δ 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}}}},

where, c is the velocity of light in vacuum,

Δ x

,

Δ y

and

Δ z

are the discrete intervals along x, y and z axis, respectively.

2.2 Design of 3D ADE-PML absorption boundary condition

The ADE-PML is employed as the boundary condition to truncate the spatial domain for the near-field simulation. The ADE-PML for MRTD scheme is first proposed by Liu et al for 2D space [35], now in our model, it has been generalized for 3D case and improved by exponential differential scheme since light attenuates rapidly in ADE-PML. Similar to section 2.1, for simplicity, we just taken E_x for example to introduce its calculation principle.

In frequency domain, the Maxwell’s equations in the stretched coordinate space can be expressed as

j ω ε E + σ E = \nabla_{s} \times H, - j ω μ H - σ_{m} H = \nabla_{s} \times E;

where

\nabla_{s} = \frac{1}{s_{x}} \frac{\partial}{\partial x} \hat{x} + \frac{1}{s_{y}} \frac{\partial}{\partial y} \hat{y} + \frac{1}{s_{z}} \frac{\partial}{\partial z} \hat{z}

, s_i (i = x, y or z) is the stretched-coordinate metric, which is proposed to be

s_{i} = κ_{i} + \frac{σ_{i}}{α_{i} + j ω ε_{0}} .

According to the equation above, the reciprocal of

s_{i}

can be written as

\frac{1}{s_{i}} = \frac{1}{κ_{i}} + \frac{σ_{i}^{'}}{α_{i}^{'} + j ω ε_{i}^{'}}

, where

σ_{i}^{'} = - σ_{i}^{}

,

ε_{i}^{'} = ε_{0} κ_{i}^{2}

,

α_{i}^{'} = κ_{i}^{2} α_{i} + κ_{i} σ_{i}

,

α_{i}

is a nonnegative number, and

κ_{i} \geq 1

. Substituting it into Eq. (16), then the Maxwell’s equations for component E_x can be written as

\frac{1}{κ_{y}} \frac{\partial H_{z}}{\partial y} - \frac{1}{κ_{z}} \frac{\partial H_{y}}{\partial z} + ψ_{e x y} - ψ_{e x z} = j ω ε E_{x} + σ E_{x};

ψ_{e x y} = \frac{σ_{y}^{'}}{α_{y}^{'} + j ω ε_{y}^{'}} \frac{\partial}{\partial y} H_{z}, ψ_{e x z} = \frac{σ_{z}^{'}}{α_{z}^{'} + j ω ε_{z}^{'}} \frac{\partial}{\partial z} H_{y};

where

ψ_{e x y}

and

ψ_{e x z}

are two auxiliary variables for simplification. Further, transforming Eq. (18) and Eq. (19) into the time domain and we can obtain the following equations, expressed as

\frac{1}{κ_{y}} \frac{\partial H_{z}}{\partial y} - \frac{1}{κ_{z}} \frac{\partial H_{y}}{\partial z} + ψ_{e x y} - ψ_{e x z} = ε \frac{\partial}{\partial t} E_{x} + σ E_{x};

ε_{y}^{'} \frac{\partial}{\partial t} ψ_{e x y} + α_{y}^{'} ψ_{e x y} = σ_{y}^{'} \frac{\partial}{\partial y} H_{z}, ε_{z}^{'} \frac{\partial}{\partial t} ψ_{e x z} + α_{z}^{'} ψ_{e x z} = σ_{z}^{'} \frac{\partial}{\partial z} H_{y};

Discretizing the differential equations in time domain by exponential differential method, then Eq. (20) and Eq. (21) can be modified as

E_{x}^{n + 1} = C A (m) \cdot E_{x}^{n} + C B (m) \cdot {(\frac{1}{κ_{y}} \frac{\partial H_{z}}{\partial y} - \frac{1}{κ_{z}} \frac{\partial H_{y}}{\partial z} + ψ_{e x y} - ψ_{e x z}) |}^{n + 1 / 2};

ψ_{e x y}^{n + 1} = ψ_{e x y}^{n} \exp (- \frac{α_{y}^{'}}{ε_{y}^{'}} Δ t) + \frac{σ_{y}^{'}}{α_{y}^{'}} [1 - \exp (- \frac{α_{y}^{'}}{ε_{y}^{'}} Δ t)] {\frac{\partial}{\partial y} H_{z} |}^{n + 1 / 2};

ψ_{e x z}^{n + 1} = ψ_{e x z}^{n - 1} \exp (- \frac{α_{z}^{'}}{ε_{z}^{'}} Δ t) + \frac{σ_{z}^{'}}{α_{z}^{'}} [1 - \exp (- \frac{α_{z}^{'}}{ε_{z}^{'}} Δ t)] {\frac{\partial}{\partial z} H_{y} |}^{n + 1 / 2} .

Similar to section 2.1, expanding E_x, $ψ_{e x y}$ and $ψ_{e x z}$ by Daubechies scalar functions and substituting them into Eq. (22) - Eq. (24), then these equations above can be further discretized in space domain by using Galerkin principle, written as

\begin{array}{l} E_{i + 1 / 2, j, k}^{ϕ x, n + 1} = C A_{i + 1 / 2, j, k} E_{i + 1 / 2, j, k}^{ϕ x, n} + C B_{i + 1 / 2, j, k} [\frac{1}{κ_{y} Δ y} \sum_{l = - L s}^{L s - 1} a (l) H_{i + 1 / 2, j + l + 1 / 2, k}^{ϕ z, n + 1 / 2} - \frac{1}{κ_{z} Δ z} H_{i + 1 / 2, j, k + l + 1 / 2}^{ϕ y, n + 1 / 2}] \\ + C B_{i + 1 / 2, j, k} (ψ_{e x y, i + 1 / 2, j, k}^{ϕ x, n + 1} - ψ_{e x z, i + 1 / 2, j, k}^{ϕ x, n + 1}); \end{array}

ψ_{e x y, i + 1 / 2, j, k}^{ϕ x, n + 1} = b_{y} ψ_{e x y, i + 1 / 2, j, k}^{ϕ x, n} + \frac{a_{y}}{Δ y} \sum_{l = - L s}^{L s - 1} a (l) H_{i + 1 / 2, j + l + 1 / 2, k}^{ϕ z, n + 1 / 2};

ψ_{e x z, i + 1 / 2, j, k}^{ϕ x, n + 1} = b_{z} ψ_{e x z, i + 1 / 2, j, k}^{ϕ x, n} + \frac{a_{z}}{Δ z} \sum_{l = - L s}^{L s - 1} a (l) H_{i + 1 / 2, j, k + l + 1 / 2}^{ϕ y, n + 1 / 2};

where

a_{i}

and

b_{i}

(

i = x, y, z

) have the forms as

b_{i} = \exp (- \frac{α_{i}^{'}}{ε_{i}^{'}} Δ t), a_{i} = \frac{σ_{i}^{'}}{α_{i}^{'}} [1 - \exp (- \frac{α_{i}^{'}}{ε_{i}^{'}} Δ t)] .

ADE-PML absorption boundary is truncated by PEC boundary in our model. Since the iteration of the field components near the PEC boundary needs the field components outside the computational domain, image principle is applied to obtain these electromagnetic components outside the simulation space.

2.3 Introducing incident wave by TF-SF Technique

Incident wave is introduced into the computational area by the TF-SF Technique proposed by Gao in our model [31]. In this scheme, computational domain is splited into the total field and scattering field domains by Connection Boundary (CB), where in the total field domain, the electromagnetic field is regarded as the superposition of scattering field and incident field [31], expressed as

E = E_{i n c} + E_{s}, H = H_{i n c} + H_{s};

in which,

E_{i n c}

and

H_{i n c}

denote the incident electromagnetic field,

E_{s}

and

H_{s}

are the scattering field. Similar to FDTD, TF-SF Technique for MRTD also introduces the incident field by modifying the field value in the regions around CB, however it is more complicated than FDTD scheme because of the expansion of the MRTD basis functions. According to the iterative equations of MRTD scheme, the regions needed to be corrected by the incident field are the space domains around the CB rather than the CB itself for FDTD. These “correction regions” are generally regularly distributed outside the scatterer, thus the contour of the scatterer needs not to be considered, which is a great advantage over PSTD.

Similar to section 2.1 and 2.2, we just take E_x as the example to introduce the implementation of TF-SF Technique. From the iterative equation of the field component E_x, it can be found that incident light should be introduced by adding or subtracting the incident field from their initial field value in 12 space regions around CB, as shown in Fig. 2, where regions 1 to 4 are the correction regions for scattering field domain, regions 5 to 12 are the correction regions for the total field. For simplicity, we only give the computational methods for region 1 and region 5 as representative.

Fig. 2 Correction regions of the field component E_x for the TF-SF Technique in MRTD scheme.

Download Full Size | PDF

The range of region 1 is $i \in [I_{\min}, I_{\max} - 1]$ , $j \in [J_{\min}, J_{\max}]$ , $k \in [K_{\max} + 1, K_{\max} + L s - 1]$ , since this region locates in scattering field domain, all the magnetic field components of Eq. (13) inside the CB should subtract the incident field value correspondingly. So the iterative equation for E_x in this region can be modified as

\begin{array}{l} _{s} E_{i + 1 / 2, j, k}^{ϕ x, n + 1} =_{s} E_{i + 1 / 2, j, k}^{ϕ x, n + 1} + \frac{Δ t}{ε Δ y} \sum_{l = - L s}^{L s - 1} a {(l)}_{s} H_{i + 1 / 2, j + l + 1 / 2, k}^{ϕ z, n + 1 / 2} \\ - \frac{Δ t}{ε Δ z} [\sum_{l = K_{\max} - k}^{L s - 1} a {(l)}_{s} H_{i + 1 / 2, j, k + l + 1 / 2}^{ϕ y, n + 1 / 2} + \sum_{l = - L s}^{K_{\max} - k - 1} a (l) (_{tol} H_{i + 1 / 2, j, k + l + 1 / 2}^{ϕ y, n + 1 / 2} -_{inc} H_{i + 1 / 2, j, k + l + 1 / 2}^{ϕ y, n + 1 / 2})], \end{array}

The range of region 5 is $i \in [I_{\min}, I_{\max} - 1]$ , $j \in [J_{\min} - L s + 1, J_{\min} - 1]$ , $k \in [K_{\max} - L s + 1, K_{\max}]$ _. Because the electric components of this region belong to the total field, the incident field should be added to the field components outside the CB correspondingly, then iterative equation of E_x in this correction region can be expressed as

\begin{array}{l} _{t o l} E_{i + 1 / 2, j, k}^{ϕ x, n + 1} =_{t o l} E_{i + 1 / 2, j, k}^{ϕ x, n + 1} + \frac{Δ t}{ε Δ y} [\sum_{l = - J_{\min} - j}^{L s - 1} a {(l)}_{tol} H_{i + 1 / 2, j + l + 1 / 2, k}^{ϕ z, n + 1 / 2} \\ + \sum_{l = - J_{\min} - j}^{L s - 1} a (l) (_{s} H_{i + 1 / 2, j + l + 1 / 2, k}^{ϕ z, n + 1 / 2} +_{inc} H_{i + 1 / 2, j + l + 1 / 2, k}^{ϕ z, n + 1 / 2})] \\ - \frac{Δ t}{ε Δ z} [\sum_{l = K_{\max} - k}^{L s - 1} a {(l)}_{tol} H_{i + 1 / 2, j, k + l + 1 / 2}^{ϕ y, n + 1 / 2} + \sum_{l = - L s}^{K_{\max} - k - 1} a (l) (_{s} H_{i + 1 / 2, j, k + l + 1 / 2}^{ϕ y, n + 1 / 2} +_{inc} H_{i + 1 / 2, j, k + l + 1 / 2}^{ϕ y, n + 1 / 2})] . \end{array}

In Eq. (30) and Eq. (31), the subscripts s, tol and inc represent the scattering field, total field and incident filed, respectively. Similar to FDTD, The incident electromagnetic field values in these correction regions are obtained by projecting the electromagnetic field calculated by 1D MRTD into 3D space.

2.4 Transformation of near field to far field

Since most aerosol particles are nonferromagnetic dielectric medium, so when they are illuminated by incident light, there is electric polarization vector $P (r)$ existing inside the particles, while the magnetic polarization vector $M (r) = 0$ [36]. Assuming that the incident light is a harmonic wave whose time dependence is given by $\exp (j ω t)$ , then the Helmholtz equations for aerosol particles can be derived as

(Δ + ω^{2} ε_{0} μ_{0}) E = - (ω^{2} ε_{0} μ_{0} I + \nabla \nabla) \cdot P (r),

where

Δ

is the Laplace operator,

P (r) = χ E (r) = (ε (r) / ε_{0} - 1) E (r)

,

χ

is the electric permittivity.

Solving Eq. (32) by Green-function method, then the far electric field can be derived as [11]

E (r) = E_{inc} (r) + \underset{V}{∭} \frac{\exp (- j k | r - r^{'} |)}{| r - r^{'} |} (k^{2} I + \nabla \nabla) \cdot P (r^{'}) d^{3} r^{'},

where,

k = ω \sqrt{ε_{0} μ_{0}}

, denotes the wave number;

r^{'}

is the position vector of the electric component in the scatterer. In the far field, taking the approximation of

| r - r^{'} |

=

r - r^{'} \cdot e_{r}

, and placing

\nabla

by

j k e_{r}

, then the scattered electric field can be modified as [11]

\begin{array}{l} E_{s c a} (r) = E (r) - E_{inc} (r) \\ = \frac{k^{2} \exp (- j k r)}{4 π r} \underset{V}{∭} [\frac{ε (r^{'})}{ε_{0}} - 1] {E (r^{'}) - e_{r} \cdot [e_{r} \cdot E (r^{'})]} \exp (j k e_{r} \cdot r^{'}) d^{3} r^{'}, \end{array}

where e_r is a unit vector in the scattering direction. Since the computational area is discrete, the continuous integral of Eq. (34) should be discretized as well. The discretization scheme applied in our model is shown as follows

\begin{array}{l} E_{s c a}^{x} (r) = \frac{k^{2} \exp (- j k r)}{4 π r} \sum_{i^{'}}^{} \sum_{j^{'}} \sum_{k^{'}} [\frac{\bar{ε} (i^{'}, j^{'}, k^{'})}{ε_{0}} - 1] {{\bar{E}}_{x}^{} (i^{'}, j^{'}, k^{'}) - \sin θ_{o b} \cos φ_{o b} \cdot T s} \\ \cdot \exp (j k_{x} \cdot i^{'} Δ x + j k_{y} \cdot j^{'} Δ y + j k_{z} \cdot k^{'} Δ z) Δ x Δ y Δ z; \end{array}

\begin{array}{l} E_{s c a}^{y} (r) = \frac{k^{2} \exp (- j k r)}{4 π r} \sum_{i^{'}}^{} \sum_{j^{'}} \sum_{k^{'}} [\frac{\bar{ε} (i^{'}, j^{'}, k^{'})}{ε_{0}} - 1] {{\bar{E}}_{y}^{} (i^{'}, j^{'}, k^{'}) - \sin θ_{o b} \sin φ_{o b} \cdot T s} \\ \cdot \exp (j k_{x} \cdot i^{'} Δ x + j k_{y} \cdot j^{'} Δ y + j k_{z} \cdot k^{'} Δ z) Δ x Δ y Δ z; \end{array}

\begin{array}{l} E_{s c a}^{z} (r) = \frac{k^{2} \exp (- j k r)}{4 π r} \sum_{i^{'}}^{} \sum_{j^{'}} \sum_{k^{'}} [\frac{\bar{ε} (i^{'}, j^{'}, k^{'})}{ε_{0}} - 1] {{\bar{E}}_{z}^{} (i^{'}, j^{'}, k^{'}) - \cos θ_{o b} \cdot T s} \\ \cdot \exp (j k_{x} \cdot i^{'} Δ x + j k_{y} \cdot j^{'} Δ y + j k_{z} \cdot k^{'} Δ z) Δ x Δ y Δ z; \end{array}

where,

(θ_{o b}, φ_{o b})

denotes the scattering direction,

\bar{ε}

is the averaged electric permittivity,

{\bar{E}}_{x}^{}

,

{\bar{E}}_{y}^{}

and

{\bar{E}}_{z}^{}

are the electric components averaged over adjacent Yee cells. k = ke_r = (k_x, k_y, k_z), is the wave vector, Ts is a temporary variable for simplification, which can be written as

T s = \sin θ_{o b} \cos φ_{o b} {\bar{E}}_{x}^{} (i^{'}, j^{'}, k^{'}) + \sin θ_{o b} \sin φ_{o b} {\bar{E}}_{y}^{} (i^{'}, j^{'}, k^{'}) + \cos φ_{o b} {\bar{E}}_{z}^{} (i^{'}, j^{'}, k^{'}) .

Having obtained the electric field in Cartesian coordinate, then we can transform them into the spherical coordinate system, expressed as

[\begin{matrix} E_{s c a}^{φ} (r) \\ E_{s c a}^{θ} (r) \\ E_{s c a}^{r} (r) \end{matrix}] = U^{T} [\begin{matrix} E_{s c a}^{x} (r) \\ E_{s c a}^{y} (r) \\ E_{s c a}^{z} (r) \end{matrix}] = [\begin{matrix} - \sin φ_{i} & \cos φ_{i} & 0 \\ \cos θ_{i} \cos φ_{i} & \cos θ_{i} \sin φ_{i} & - \sin θ_{i} \\ \sin θ_{i} \cos φ_{i} & \sin θ_{i} \sin φ_{i} & \cos θ_{i} \end{matrix}] [\begin{matrix} E_{s c a}^{x} (r) \\ E_{s c a}^{y} (r) \\ E_{s c a}^{z} (r) \end{matrix}] .

It should be noted that the result of $E_{s c a}^{r} (r)$ should approximately equal to 0, and it can be applied for the validation of the near-to-far transformation. Based on the far electric field we have obtained, the scattering amplitude matrix and phase matrix can be calculated subsequently, as illustrated in section 2.5.

2.5 Calculation models of particle’ s scattering properties

Derivation of absorption cross section. Absorption cross section describes the ability to absorb the electromagnetic wave and transform it into thermal energy. Similar to the derivation process proposed by Yang and Liou [32], the energy absorbed by the aerosol particles $P_{a b}$ can be calculated by integrating the Poynting vector $S (r)$ over a closed surface with scattering particle inside, expressed as

P_{a b} = - Re [\iint_{S} n_{s} \cdot S (r) d^{2} r] = - Re [∭_{V} \nabla \cdot S (r) d^{3} r] = - \frac{1}{2} Re [∭_{V} \nabla \cdot (E (r) \times H^{*} (r)) d^{3} r],

in which,

n_{s}

is the outward-pointing unit normal vector of the integral surface, Re(.) denotes the real part. Simplifying the equation above by Maxwell’s curl equations, then we can obtain

P_{a b} = \frac{1}{2} Re ∭_{V} j ω (ε (r) E (r) \cdot E^{*} (r) - μ (r) H (r) \cdot H^{*} (r)) - σ (r) E (r) \cdot E^{*} (r) d^{3} r .

From the equation above, it can be easily found that the energy absorbed by aerosol particles can be presented as

P_{a b} = \frac{1}{2} ∭_{V} σ (r) E (r) \cdot E^{*} (r) d^{3} r .

As can be seen, the energy absorbed by particles is transformed into Joule heat, which is in accordance with the definition of absorption cross section. After $P_{a b}$ is obtained, Then the absorption cross section $C_{a b}$ can be calculated by normalizing $P_{a b}$ by the incident flux $P_{i n c}$ , written as

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

where,

Z_{0}

is the wave impedance,

E_{i n c}

is the intensity vector of incident electric field.

Derivation of extinction cross section. In the frequency domain, the total Poynting vector can be decomposed into three parts: the Poynting vector of incident field $S_{i n c}$ , Poynting vector of scattering field $S_{s c a}$ and Poynting vector of extinction field $S_{e x t}$ , where the specific form of $S_{e x t}$ can be expressed as [37]

S_{e x t} (r) = \frac{1}{2} (E_{i n c} (r) \times H_{s}^{*} (r) + E_{s} (r) \times H_{i n c}^{*} (r)) .

Similar to the implementation introduced by Yang [32], integrating $S_{e x t}$ over a closed surface surrounding the scatterer and taking its real part, and then the extinction power $P_{e x t}$ can be presented as

P_{e x t} = - Re (\iint_{S} S_{e x t} (r) \cdot n_{s} d^{2} r) = - \frac{1}{2} Re {∭_{V} \nabla \cdot (E_{i n c} (r) \times H_{s}^{*} (r) + E_{s} (r) \times H_{i n c}^{*} (r)) d^{3} r} .

Using the vector equation $\nabla \cdot (A \times B) = B \cdot \nabla \times A - A \cdot \nabla \times B$ for simplification, then Eq. (45) can be rewritten as

P_{e x t} = - \frac{1}{2} Re {∭_{V} [\begin{array}{l} j ω μ_{0} (H_{i n c} (r) \cdot H_{s}^{*} (r) + H_{i n c}^{*} (r) \cdot H_{s} (r)) \\ + j ω (ε E_{i n c} (r) \cdot E_{s}^{*} (r) + ε_{0} E_{i n c}^{*} (r) \cdot E_{s} (r) - \frac{σ}{j ω} E_{i n c} (r) \cdot E_{s}^{*} (r)) \end{array}] d^{3} r} .

Because the items of magnetic field are in Eq. (46) conjugate, it can be further presented as

P_{e x t} = - \frac{1}{2} {∭_{V} Im [ω (ε E_{i n c} (r) \cdot E_{s}^{*} (r) - ε_{0} E_{i n c} (r) \cdot E_{s}^{*} (r) + \frac{σ}{j ω} E_{i n c} (r) \cdot E_{s}^{*} (r))] d^{3} r} .

Having obtained extinction power

P_{e x t}

, the extinction cross section

C_{e x t}

can be calculated by normalizing it by the incident flux, given as

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

Based on the extinction cross section $C_{e x t}$ and absorption cross section $C_{a b}$ we have obtained, the scattering cross section $C_{s c}$ can be calculated by

C_{s c} = C_{e x t} - C_{a b} .

The extinction, absorption and scattering efficiencies can also be calculated by normalizing the corresponding cross section by particle’s geometrical cross section [2].

Calculation of the scattering amplitude matrix and phase matrix. Scattering amplitude matrix $S (θ)$ is a parameter describing the relationship between the incident and scattering electric field, which is defined as [38]

[\begin{matrix} E_{s c a}^{⊥} \\ E_{s c a}^{/ /} \end{matrix}] = \frac{\exp j (- k r + z)}{j k r} S (θ) [\begin{matrix} E_{i n c}^{⊥} \\ E_{i n c}^{/ /} \end{matrix}] = \frac{\exp j (- k r + z)}{j k r} [\begin{matrix} S_{1} & S_{4} \\ S_{3} & S_{2} \end{matrix}] [\begin{matrix} E_{i n c}^{⊥} \\ E_{i n c}^{/ /} \end{matrix}],

in which, θ is the scattering angle,

E_{i n c}^{⊥}

and

E_{i n c}^{/ /}

denote the incident electric field perpendicular and parallel to the scattering plane.

E_{s c a}^{⊥}

and

E_{s c a}^{/ /}

are vertically-polarized and horizontally-polarized components of the scattering field, which satisfy the following relationship:

E_{s c a}^{⊥} = - E_{s c a}^{φ}

and

E_{s c a}^{/ /} = E_{s c a}^{θ}

[11]. According to Eq. (50),

S (θ)

can be obtained by the following way: First, setting the incident light perpendicular and parallel to the scattering plane (the intensity of the incident field satisfies

E_{i n c}^{⊥} = E_{i n c}^{/ /} = E_{0}

), respectively, then the scattering electric field

{(E_{s c a, ⊥}^{⊥}, E_{s c a, ⊥}^{/ /})}^{T}

and

{(E_{s c a, / /}^{⊥}, E_{s c a, / /}^{/ /})}^{T}

can be calculated by the methods described above. Based on Eq. (50), the relationship between the incident and scattering field can be further established by matrix operation, written as

[\begin{matrix} E_{s c a, ⊥}^{⊥} & E_{s c a, / /}^{⊥} \\ E_{s c a, ⊥}^{/ /} & E_{s c a, / /}^{/ /} \end{matrix}] = \frac{\exp j (- k r + z)}{j k r} [\begin{matrix} S_{1} & S_{4} \\ S_{3} & S_{2} \end{matrix}] [\begin{matrix} E_{i n c}^{⊥} & 0 \\ 0 & E_{i n c}^{/ /} \end{matrix}] .

Solving the equation above, then the scattering amplitude matrix can be obtained. According to the relationship between the amplitude matrix and phase matrix F(θ), the elements of F(θ) can easily calculated from the known amplitude matrix, readers can refer to literature [11] and [37] for detail.

3. Model simulation and validation

All the procedure of the MRTD model is implemented in Fortran, and MPI repeated non-blocking communication technique is applied to accelerate the calculation speed as well. The parallelization scheme is similar to that proposed by Liu et al [39], where, the computational domain is averagely divided into several small domains along z axis according to the number of the processors. The simulation processes of different processors are kept in step by data exchange, and the method for data exchange process is shown in Fig. 3. To validate the accuracy of our model, its simulation results are compared against Lorenz-Mie, Aden-Kerker scattering theory and T-matrix method for four cases: (1) sphere particles with its size approximately equal and much larger than light wavelength; (2) concentric sphere; (3) ellipsoid particle; (4) cylinder particle.

Fig. 3 The data exchange scheme for parallel MRTD scattering model.

Download Full Size | PDF

3.1 Model validation for spherical particles

Case 1: Model validation for small spherical particle. In this case, the phase matrix simulated by MRTD scattering model is compared against that obtained by Lorenz-Mie theory. In the simulation process, the wavelength of incident light is set as λ = 0.633μm, particle’s radius and refractive index are taken as 0.8μm and 1.33-10⁻⁹i, respectively. The computational domain is discretized by a grid mesh with a spatial resolution of λ/60.

The nonzero elements of the phase matrix obtained by MRTD and Lorenz-Mie scattering theory are presented in Fig. 4. As can be seen, an excellent agreement is achieved between the results of different scattering models, indicating our model is accurate and reliable for the small spherical particle. For element F₁₁ (i.e., scattering phase function), its relative simulation errors fluctuate with scattering angle notably, with their values ranging from near-zero in forward directions to 11.8% at the scattering angles near 160°. For F₁₂/ F₁₁, F₃₄/F₁₁ and F₄₄/F₁₁, the absolute simulation errors are less than 0.05 for scattering angles smaller than 130°, however for backward scattering directions, their simulation errors are much larger than that at small scattering angles though the curves obtained by different models are coincided with each other approximately, where the maximum absolute simulation errors of F₁₂/ F₁₁, F₃₄/F₁₁ and F₄₄/ F₁₁ approach 0.11, 0.13 and 0.15, respectively. The phenomenon can be explained by stair approximation error in particle construction process. Since the particle is approximated by cubic cells, its real shape is actually different from the actual sphere we have assumed, and because the scattering phase matrix is more sensitive to particle’s shape in backscattering directions than forward directions, larger simulation errors are caused at large scattering angles correspondingly.

Fig. 4 Comparison of the scattering phase matrices simulated by MRTD and Lorenz-Mie theory (r = 0.5μm). The left panel presents the variation of phase matrix elements with scattering angle, and the right panel shows the variation of the simulation errors, where the simulation errors of F₁₁ are evaluated by the relative error; for F₁₂/ F₁₁, F₃₄/F₁₁ and F₄₄/F₁₁, the errors are described by absolute errors.

Download Full Size | PDF

Case 2: Model validation for large spherical particle. The nonzero phase matrix elements of a sphere with a size parameter x = 100 is presented in Fig. 5. In this case, the incident light wavelength λ is set as 0.86μm, and the spatial resolution is set to be λ/16. The refractive index of the particle is taken as 1.52-0.8i.

Fig. 5 Comparison of the scattering phase matrices simulated by MRTD and Lorenz-Mie theory for the spherical particle with a size parameter of 100. The left panel presents the curves of phase matrix elements with scattering angle, and the right panel shows the simulation error.

Download Full Size | PDF

As shown in the figure, the simulation results of MRTD model agree with the exact solution of Lorenz-Mie theory well in terms of their variation patterns. In forward directions, the relative errors of F₁₁ are less 20% generally; however for scattering angles ranging from 90° to 180°, the simulation errors are much larger, and the maximum error reaches 25% approximately. Contrary to the results in case 1, F₁₂/ F₁₁ and F₃₄/ F₁₁ have larger simulation errors (the maximum absolute errors is up to 0.0822 and 0.0601, respectively) in forward scattering directions than backward directions. The simulation errors of F₄₄/ F₁₁ has a similar variation pattern to F₁₁, and its absolute simulation errors ranges from near zero in forward directions to near 0.1 at 78°and 176°.

Case 3: validation for integral scattering parameters. To further validate the accuracy of our model, the integral scattering properties as extinction, absorption and scattering efficiencies obtained by MRTD model are compared with Lorenz-Mie theory. The simulation is performed for spherical particles with a size parameter of 5, 10, 20, 40 and 80, where the light wavelength λ is set as 0.8μm, and the refractive index of the aerosol particle is taken as 1.33-0.005i. The relative simulation errors are calculated for different cases, and the results are represented in Table 1.

Table 1. Comparison of the integral scattering parameters simulated by Lorenz-Mie theory and parallelized MRTD scattering model.

View Table | View all tables in this article

As shown in the Table, the integral scattering properties simulated by MRTD and Lorenz-Mie theory show an excellent agreement, where the simulation errors of extinction and absorption efficiencies for the particle with a size parameter of 10 are −0.4891% and −1.6933% respectively, and the error for single scattering albedo is only 0.1168%, which validates the accuracy of our model. On the whole, the simulation errors of the integral scattering properties are smaller than those of phase matrix.

3.2 Model validation for particles with inhomogeneous compositions

Model validation for inhomogeneous particles is performed for a concentric sphere with a high absorptive spherical core and a pure-scattering spherical shell. In this test, the wavelength of incident light is set as λ = 0.55μm, the radius of the inner and coated sphere are taken to be 0.25μm and 0.5μm, and their refractive indices are set as 1.52-0.75i and 1.33-0.00i, respectively. The computational space containing the particle is discretized by cubic cells with an edge length of λ/68. Then the phase matrix is calculated by MRTD scattering model and is compared with that by Aden-Kerker theory, as shown in Fig. 6.

Fig. 6 Comparison of the scattering phase matrices simulated by MRTD and Aden-Kerker theory for the concentric sphere. The left panel presents the curves of phase matrix elements with scattering angle, and the right panel shows their simulation errors.

Download Full Size | PDF

As can be seen, the curves of the phase matrix elements simulated by different models are coincided with each other approximately, especially in the forward scattering directions. For most scattering angles, the relative simulation errors of F₁₁ are less than 5%, and the absolute errors of F₁₂/F₁₁, F₃₄/F₁₁ and F₄₄/F₁₁ are less than 0.05. Overall, simulation errors are gradually increases with the increasing of scattering angle in terms of their variation characteristics, which is similar to the results for the small spherical particle.

3.3 Model validation for spheroidal particle

For nonspherical particles case, the scattering phase matrix obtained by our model is firstly compared with T-matrix method for spheroidal particle, as shown in Fig. 7. In this test, the particle is set as an oblate spheroid with its horizontal semi-axis a = 1.0μm and its vertical semi-axis b = 0.5μm, the refractive index is taken as 1.414-0.00i. The direction of incident light is parallel to vertical axis, and the light wavelength is set as λ = 0.86μm. The computational space is discretized by a grid mesh with a spatial resolution of λ/68.

Fig. 7 Comparison of the phase matrix elements simulated by MRTD and T-matrix method for spheroidal particle. The left panel presents the variation of phase matrix elements with scattering angle, and the right panel shows their relative simulation errors.

Download Full Size | PDF

As can be seen, a high agreement is achieved between the results of MRTD model and those of T-matrix method, where the relative errors of F₁₁ are generally less than 5% except for the scattering angles larger than 160°; the averaged absolute simulation error of F₁₂/F₁₁, F₃₄/F₁₁ and F₄₄/F₁₁ over all scattering angles are less than −0.014, −0.010 and 0.008, respectively. Similar to the results of spherical particles, large simulation errors are mainly distributed at large scattering angles.

The results of the integral scattering parameters simulated for particles with different sizes and refractive indices by different scattering models are given in Table 2. In this computational process, the wavelength of incident light is taken as 0.5μm, and the aspect ratio of the spheroidal particle is set as a/b = 0.5, the equivalent radii of the particles are obtained by equal-volume principle. From the Table, it can be found that the results simulated by MRTD model show a high agreement with those of T-matrix theory, and for most cases, the relative errors of both the extinction and absorption efficiencies are smaller than 1%.

Table 2. Comparison of the integral scattering parameters of spheroidal particles simulated by T-matrix theory and parallelized MRTD scattering model.

View Table | View all tables in this article

3.4 Model validation for Model validation for cylindrical particles

The MRTD scattering model is further validated against T-Matrix method for circular cylindrical particle. In the simulation process, the spatial resolution is the same as that in section 3.3 for the spheroidal particle. The diameter and length of the cylinder are all set as 1.0μm, and the refractive index is taken as 1.33-0.000i; the light wavelength is taken as 0.532μm, and its direction is perpendicular to the undersurface of the cylinder. Then the particle’s phase matrix and its simulation errors are calculated and presented in Fig. 8.

Fig. 8 Comparison of scattering phase matrix elements simulated by MRTD and T-matrix method for circular cylindrical particle. The left panel presents the curves of phase matrix elements with scattering angle, and the right panel shows the variation of their simulation errors.

Download Full Size | PDF

As can be seen, the simulation results of MRTD model and T-matrix theory show a high consistency, where the relative simulation errors of F₁₁ are less than 5% for most scattering angles, and the maximum absolute errors of F₁₂/F₁₁, F₃₄/F₁₁ and F₄₄/F₁₁ in forward scattering direction (0°<θ<90°) are less than 0.04, 0.08 and 0.05, respectively. Similar to the results of spheroidal and spherical particles, as the scattering angle approaches 180°, the simulation errors increase notably due to stair approximation error in the particle construction process. On the whole, MRTD model can simulate the scattering process of cylinder with high precision.

The simulation accuracy of the integral scattering properties by MRTD is also investigated for cylindrical particles. In this test, light wavelength is set as 0.55μm, and ratio of diameter and length of the cylinders (D/L) are all taken to be 1.0, the simulation results are demonstrated in Table 3. It is seen that the simulation accuracy of extinction and absorption efficiencies is much higher than that of phase matrix, and the relative simulation errors are less than 1% for most cases. It can also be found that higher accuracy is achieved by MRTD model for particles with their size comparable to the light wavelength than that for other particles.

Table 3. Comparison of the integral scattering parameters of cylindrical particles simulated by T-matrix theory and parallelized MRTD scattering model.

View Table | View all tables in this article

4. Parallel computational efficiency analysis

The efficiency of parallel computation is usually evaluated by the parallel acceleration ratio S_p and computational efficiency for single process η, where the parallel acceleration ratio S_p is defined as

S_{p} = T_{1} / T_{n},

in which, T₁ and T_n are the computational time for one processor and n processors, respectively. η is a parameter describing the time ratio that one processor is effectively used for computation, which can be calculated by

η = S_{p} / n .

The evaluation of the parallel computational efficiency is performed for the spherical particle with a size parameter of 80. In simulation process, the light wavelength is taken as λ = 0.633μm, the computational area is discretized by a grid mesh of λ/20, and the refractive index is set to be 1.52-0.008i. Scattering process is simulated for different number of processors (32bit 3.1GHz), and their computational time, parallel acceleration ratio and parallel computational efficiency are presented in the Table 4.

Table 4. the parallel computational efficiencies for different number of processors.

View Table | View all tables in this article

As shown in the table, the computational time decreases notably with the increasing of the number of processors, where the simulation time for 8 processors decrease to 1/6 of that for 1 processor approximately. It also should be noted that η decreases as the number of processors increases; this can be explained by the increasing of the communication time between different processors.

5. Conclusion

Scattering process of nonspherical aerosols and ice crystal plays an important role in radiative transfer since it can alter the transmission, reflection and absorption ability of atmospheric system. Due to the inadequate understanding of their scattering properties, there still exist considerable uncertainties in the radiative transfer simulation. To this end, a new scattering computation model based on MRTD method is proposed, and to obtain the higher computational efficiency, it is further parallelized by MPI technique. The model can be roughly decomposed into three parts: Near field calculation module, Near-to-far transformation module and computational module for scattering properties. By using this model, such scattering parameters as phase matrix, extinction efficiency, absorption efficiency and single scattering albedo can be simulated. The precision of our model is validated against the well-tested scattering theories like Lorenz-Mie, Aden-Kerker and T-matrix scattering theories for three conditions: spherical particles, inhomogeneous particles and nonspherical particles. The parallel computational efficiency of our model is analyzed as well in the last section. Several conclusions are drawn as follows:

The phase matrix simulated by parallel MRTD scattering model shows an excellent agreement with those of Lorenz-Mie, Aden-Kerker and T-matrix scattering theories. For element F₁₁, the relative simulation errors are less than 5% for most scattering angles. For F₁₂/ F₁₁, F₃₄/ F₁₁ and F₄₄/ F₁₁, their absolute simulation errors are generally less than 0.05. At large scattering angles, the simulation errors are much larger than that in forward scattering directions, where the relative errors of F₁₁ can reach 10%~20% for several scattering angles. The reason is that these particles are approximated by cubic cells, and their actual shapes are different from the ideal shapes we have assumed.

The computational accuracy of the integral scattering parameters as extinction efficiency and absorption efficiency is higher than that of phase matrix overall, where the simulation errors of extinction and absorption efficiencies for the spherical particle with a size parameter of 10 are −0.4891% and −1.6933%, respectively.

With the increasing of processors, the parallel acceleration ratio is improved dramatically; on the contrary, the computational efficiency for single processor decreases to some extent, this phenomenon can be explained by the increasing of the communication time between the processors.

Compared with the conventional FDTD, MRTD can achieve high simulation accuracy with less girds per wavelength because of its excellent numerical dispersion property, thus less computational time and memory is needed, which is an advantage similar to PSTD. Since the incident field can only be introduced by TF-SF technique in MRTD scheme, the implementation of MRTD is much simpler for particles with large size and complex shapes than PSTD.

Funding

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

Acknowledgments

Thanks 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. The authors would like to acknowledge the reviewers for their helpful remarks as well.

References and links

1. K. N. Liou and Y. Takano, “Light scattering by nonspherical particles: Remote sensing and climatic implications,” Atmos. Res. 31(4), 271–298 (1994). [CrossRef]

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

3. H. Iwabuchi and P. Yang, “Temperature dependence of ice optical constants: Implications for simulating the single-scattering properties of cold ice clouds,” J. Quant. Spectrosc. Radiat. Transf. 112(15), 2520–2525 (2011). [CrossRef]

4. M. I. Mishchenko, “Electromagnetic scattering by nonspherical particles:A tutorial review,” J. Quant. Spectrosc. Radiat. Transf. 110(11), 808–832 (2009). [CrossRef]

5. 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]

6. H. Yi, X. Ben, and H. Tan, “Transient radiative transfer in a scattering slab considering polarization,” Opt. Express 21(22), 26693–26713 (2013). [CrossRef] [PubMed]

7. IPCC: Climate Change 2007, “Intergovernmental Panel of Global Climate Change,” (2007).

8. T. Cheng, X. Gu, Y. Wu, H. Chen, and T. Yu, “The optical properties of absorbing aerosols with fractal soot aggregates: Implications for aerosol remote sensing,” J. Quant. Spectrosc. Radiat. Transf. 125, 93–104 (2013). [CrossRef]

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

10. Y. Liu, W. P. Arnott, and J. Hallett, “Particle size distribution retrieval from multispectral optical depth:Influences of particle nonsphericity and refractive index,” J. Geophys. Res. 104(D24), 31753–31762 (1999). [CrossRef]

11. M. I. Mishchenko, J. W. Hovenier, and L. D. Travis, Light Scattering by Nonspherical Particles, Theory, Measurements, and Application (Academic Press, 2000).

12. H. Iwabuchi and T. Suzuki, “Fast and accurate radiance calculations using truncation approximation for anisotropic scattering phase functions,” J. Quant. Spectrosc. Radiat. Transf. 110(17), 1926–1939 (2009). [CrossRef]

13. 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]

14. J.-Q. Zhao, G. Shi, H. Che, and G. Cheng, “Approximations of the scattering phase functions of particles,” Adv. Atmos. Sci. 23(5), 802–808 (2006). [CrossRef]

15. P. Yang, H. Wei, H.-L. Huang, B. A. Baum, Y. X. Hu, G. W. Kattawar, M. I. Mishchenko, and Q. Fu, “Scattering and absorption property database for nonspherical ice particles in the near- through far-infrared spectral region,” Appl. Opt. 44(26), 5512–5523 (2005). [CrossRef] [PubMed]

16. 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]

17. 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]

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

19. 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]

20. L. Liu, M. I. Mishchenko, B. Cairns, B. E. Carlson, and L. D. Travis, “Modeling single-scattering properties of small cirrus particles by use of a size-shape distribution of ice spheroids and cylinders,” J. Quant. Spectrosc. Radiat. Transf. 101(3), 488–497 (2006). [CrossRef]

21. A. Quirantes, “A T-matrix method and computer code for randomly oriented, axially symmetric coated scatterers,” J. Quant. Spectrosc. Radiat. Transf. 92(3), 373–381 (2005). [CrossRef]

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

23. B. T. Draine, “Discrete-dipole approximation and its application to interstellar graphite grains,” Astrophys. J. 333, 848–872 (1988). [CrossRef]

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

25. P. Yang and K. N. Liou, “Light scattering by hexagonal ice crystals: comparison of finite-difference time domain and geometric optics models,” J. Opt. Soc. Am. A 12(1), 162 (1995). [CrossRef]

26. 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]

27. 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]

28. S. Y. Dai, “Application of Wavelet and Multiresolution Time domain (MRTD) on Electromagnetic scattering,” (Xidian University, 2007).

29. Y. W. Cheong, Y. M. Lee, K. H. Ra, J. G. Kang, and C. C. Shin, “Wavelet-Galerkin scheme of time dependent inhomogeneous electromagnetic problems,” IEEE Microw. Guided Wave Lett. 9(8), 297–299 (1999). [CrossRef]

30. L. Massy, N. Pena, and M. M. Ney, “Dispersion characteristics and stability for an arbitray MRTD scheme with variable mesh,” Int. J. Numer. Model.: Electron. Netw. Dev. Fields 23(6), 470–491 (2010). [CrossRef]

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.

32. 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), 2072–2085 (1996). [CrossRef]

33. M. Krumpholz and L. P. B. Katehi, “MRTD: New Time-Domain Schemes Based on Multiresolution Analysis,” IEEE Trans. Microw. Theory Tech. 44(4), 555–571 (1996). [CrossRef]

34. E. M. Tentzeris, A. Cangellaris, L. P. B. Katehi, and J. Harvey, “MultiresolutionTime-Domain(MRTD)Adaptive Schemes Using Arbitrary Resolutions of Wavelets,” IEEE Trans. Microw. Tech. 50(2), 501–516 (2002). [CrossRef]

35. Y. Liu, Y. Chen, X. Xu, and Z. Liu, “Implementation and analysis of the perfectly matched layer with auxiliary differential equation for the multiresolution time-domain method,” Wuli Xuebao 62, 034101 (2013).

36. D. J. Griffiths, Introduction to Electrodynamics, 3^rd ed. (Pearson Education, Inc, 2014).

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

38. H. C. van der Hulst, Light Scattering by Small Particles (Dover Publications, 1981).

39. Y. Liu, Y. Chen, and P. Zhang, “Parallel implementation and application of the MRTD with an efficient CFS-PML,” Prog. Electromagnetics Res. 143, 223–242 (2013). [CrossRef]

x	Grid size	Extinction efficiency Q_ex			Absorption efficiency Q_ab			Single scattering albedo
x	Grid size	Mie	MRTD	RE(/%)	Mie	MRTD	RE(/%)	Mie	MRTD	RE(/%)
5	λ/50	3.5362	3.4964	−1.1255	0.1018	0.0995	−2.2298	0.9712	0.9715	0.0350
10	λ/40	2.2284	2.2175	−0.4891	0.2067	0.2032	−1.6933	0.9073	0.9084	0.1168
20	λ/40	2.1296	2.0948	−1.6341	0.3701	0.3585	−3.1343	0.8262	0.8289	0.3219
40	λ/30	2.0478	1.9997	−2.3488	0.5372	0.5227	−2.6992	0.7377	0.7386	0.1234
80	λ/20	2.0815	2.0265	−2.6423	0.7526	0.7283	−3.2289	0.6384	0.6406	0.3462

m	Radius /(μm)	Grid size	Extinction efficiency Q_ex			Absorption efficiency Q_ab
m	Radius /(μm)	Grid size	T-Matrix	MRTD	RE(/%)	T-Matrix	MRTD	RE(/%)
1.33-0.08i	0.1	λ/50	0.39056	0.38853	−0.52	0.25476	0.25338	−0.542
	0.5	λ/40	3.12926	3.12911	−4.8E-3	1.01461	1.01418	−0.042
	1.0	λ/40	2.56958	2.55594	−0.531	1.22015	1.21131	−0. 725
	2.0	λ/30	2.40848	2.39054	−0.745	1.19413	1.17958	−1.218
1.53-0.0008i	0.1	λ/50	0.37475	0.37172	−0.809	0.00287	0.00279	−2.787
	0.5	λ/40	3.52982	3.52963	−5.4E-3	0.02545	0.02552	0.275
	1.0	λ/40	2.87662	2.86278	−0.481	0.05415	0.05385	−0.554
	2.0	λ/30	2.30260	2.28645	−0.701	0.09452	0.09294	−1.672

m	Radius /(μm)	Grid size	Extinction efficiency Q_ex			Absorption efficiency Q_ab
m	Radius /(μm)	Grid size	T-Matrix	MRTD	RE(/%)	T-Matrix	MRTD	RE(/%)
1.33-0.08i	0.1	λ/50	0.39672	0.39394	−0.701	0.25564	0.25398	−0.649
	0.5	λ/40	3.06473	3.06446	−8.8E-3	1.01168	1.01153	−0.015
	1.0	λ/40	2.43360	2.43289	−0.029	1.21785	1.21496	−0.237
	2.0	λ/30	2.42828	2.40743	−0. 859	1.19871	1.18943	−0.774
1.53-0.0008i	0.1	λ/50	0.38117	0.38004	−0. 296	0.00287	0.00276	−3.833
	0.5	λ/40	2.95258	2.95248	−3.4E-3	0.02601	0.02617	0.615
	1.0	λ/40	2.86020	2.85125	−0. 313	0.05364	0.05288	−1.417
	2.0	λ/30	2.61422	2.60143	−0. 489	0.09529	0.09409	−1.259

Processors	Computational time/s	S_p	η
1	1.0783 × 10⁶	1.000	100%
2	6.7142 × 10⁵	1.606	80.3%
3	4.7670 × 10⁵	2.262	75.4%
4	3.6380 × 10⁵	2.964	74.1%
8	1.8747 × 10⁵	5.752	71.9%

x	Grid size	Extinction efficiency Q_ex			Absorption efficiency Q_ab			Single scattering albedo
x	Grid size	Mie	MRTD	RE(/%)	Mie	MRTD	RE(/%)	Mie	MRTD	RE(/%)
5	λ/50	3.5362	3.4964	−1.1255	0.1018	0.0995	−2.2298	0.9712	0.9715	0.0350
10	λ/40	2.2284	2.2175	−0.4891	0.2067	0.2032	−1.6933	0.9073	0.9084	0.1168
20	λ/40	2.1296	2.0948	−1.6341	0.3701	0.3585	−3.1343	0.8262	0.8289	0.3219
40	λ/30	2.0478	1.9997	−2.3488	0.5372	0.5227	−2.6992	0.7377	0.7386	0.1234
80	λ/20	2.0815	2.0265	−2.6423	0.7526	0.7283	−3.2289	0.6384	0.6406	0.3462

Light scattering computation model for nonspherical aerosol particles based on multi-resolution time-domain scheme: model development and validation

Abstract

1. Introduction

2. Implementation of MRTD scattering model

2.1 Basic principle of MRTD and its Courant stability condition

2.2 Design of 3D ADE-PML absorption boundary condition

2.3 Introducing incident wave by TF-SF Technique

2.4 Transformation of near field to far field

2.5 Calculation models of particle’ s scattering properties

3. Model simulation and validation

3.1 Model validation for spherical particles

3.2 Model validation for particles with inhomogeneous compositions

3.3 Model validation for spheroidal particle

3.4 Model validation for Model validation for cylindrical particles

4. Parallel computational efficiency analysis

5. Conclusion

Funding

Acknowledgments

References and links

Cited By

Figures (8)

Tables (4)

Equations (53)

Optics Express