Discrete dipole approximation for black carbon-containing aerosols in arbitrary mixing state: A hybrid discretization scheme

Discrete dipole approximation for black carbon-containing aerosols in arbitrary mixing state: A hybrid discretization scheme

Author’s Accepted Manuscript Discrete dipole approximation for black carboncontaining aerosols in arbitrary mixing state: A hybrid discretization sche...

935KB Sizes 0 Downloads 28 Views

Author’s Accepted Manuscript Discrete dipole approximation for black carboncontaining aerosols in arbitrary mixing state: A hybrid discretization scheme Nobuhiro Moteki www.elsevier.com/locate/jqsrt

PII: DOI: Reference:

S0022-4073(15)30098-4 http://dx.doi.org/10.1016/j.jqsrt.2016.01.025 JQSRT5200

To appear in: Journal of Quantitative Spectroscopy and Radiative Transfer Received date: 4 September 2015 Revised date: 20 January 2016 Accepted date: 20 January 2016 Cite this article as: Nobuhiro Moteki, Discrete dipole approximation for black carbon-containing aerosols in arbitrary mixing state: A hybrid discretization scheme, Journal of Quantitative Spectroscopy and Radiative Transfer, http://dx.doi.org/10.1016/j.jqsrt.2016.01.025 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

JQSRT_2015_103

Discrete dipole approximation for black carbon-containing aerosols in arbitrary mixing state: A hybrid discretization scheme Nobuhiro Moteki*

Graduate School of Science, Department of Earth and Planetary Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, 113-0033 Tokyo, Japan

* Corresponding author Email address: [email protected]

Last modified: January 19, 2016

1

Abstract An accurate and efficient simulation of light scattering by an atmospheric black carbon (BC)-containing aerosol—a fractal-like cluster of hundreds of carbon monomers that is internally mixed with other aerosol compounds such as sulfates, organics, and water—remains challenging owing to the enormous diversities of such aerosols’ size, shape, and mixing state. Although the discrete dipole approximation (DDA) is theoretically an exact numerical method that is applicable to arbitrary non-spherical inhomogeneous targets, in practice, it suffers from severe granularity-induced error and degradation of computational efficiency for such extremely complex targets. To solve this drawback, we propose herein a hybrid DDA method designed for arbitrary BC-containing aerosols: the monomer-dipole assumption is applied to a cluster of carbon monomers, whereas the efficient cubic-lattice discretization is applied to the remaining particle volume consisting of other materials. The hybrid DDA is free from the error induced by the surface granularity of carbon monomers that occurs in conventional cubic-lattice DDA. In the hybrid DDA, we successfully mitigate the artifact of neglecting the higher-order multipoles in the monomer-dipole assumption by incorporating the magnetic dipole in addition to the electric dipole into our DDA formulations. Our numerical experiments show that the hybrid DDA method is an efficient light-scattering solver for BC-containing aerosols in arbitrary mixing states. The hybrid DDA could be also useful for a cluster of metallic nanospheres associated with other dielectric materials.

Keywords: light scattering, discrete dipole approximation, aerosol, black carbon, mixing state

2

1. Introduction Black carbon (BC) is a light-absorbing aerosol produced by the combustion of fossil fuels and biomass. It may be the second-most important positive climate-forcing agent after carbon dioxide [1, 2, 3]. The typical morphology of pure BC aerosols is a fractal-like cluster of tens to hundreds of spherical monomers, where the monomer diameter ranges from 20 to 100 nm [4]. BC aerosols can be internally mixed with non-absorbing, or weakly absorbing, compounds such as sulfates, organics, or water through aging processes, forming BC-containing aerosols with various morphologies. The absorption cross section of a BC-containing aerosol per unit BC mass can change by several-folds, depending on its morphology [5, 6, 7]. Therefore, an efficient and reliable light scattering solver applicable to arbitrary BC-containing aerosols is necessary for modeling the optical properties of aerosols as a function of aging time [8] and characterizing the morphology of BC-containing aerosols using light-scattering techniques [9, 10, 11, 12]. The discrete dipole approximation (DDA), wherein a target is represented by a collection of interacting point electric dipoles, is a numerically exact solver of light scattering by arbitrary non-spherical inhomogeneous particles [13, 14, 15]. In research on atmospheric BC-containing aerosols, publicly available DDA codes have been used to investigate the effects of aerosols’ morphologies on their optical properties [16, 17]. Almost every DDA code, by default, assumes that the position of each point dipole is restricted to lie on the grid points of a cubic lattice (cubic-lattice DDA) as a necessary condition for applying an efficient algorithm using fast-Fourier transform (FFT) [18, 19, 20]. Without this efficient algorithm, the DDA would be prohibitively slow for particles with sizes comparable to or larger than the wavelength.

3

The cubic-lattice DDA suffers from an error caused by the surface granularity (granularity-induced error) [14], the magnitude of which is larger for a smaller number of dipoles N. Numerical convergence of this error is sluggish under a high refractive index |m| >> 1: the number of dipoles N scales to approximately |m|3 to attain a desired fractional error [14]. The granularity-induced error persists to the electrostatic limit because the internal field induced by an applied static field is sensitive to surface granularity [14, 21]. Thus, the granularity-induced error is especially problematic for a large cluster of small spheres with |m| >> 1 because a very large N is required to suppress the granularity-induced error in each spherical monomer. Moreover, the cubic-lattice DDA would no longer be efficient for a target of highly non-spherical shapes with a low packing density because the cubic-lattice DDA fully uses a three-dimensional domain of a rectangular cuboid that encloses the entire target volume [19, 20]. Therefore, the cubic-lattice DDA has difficulties in computing light scattering by a large fractal-like cluster of small spheres with |m| >> 1, such as a BC aerosol or an aggregate of metallic nanospheres. A different DDA approach that can be more effective for a large fractal-like cluster of small spheres is to approximate each monomer by a point electric dipole (monomer-dipole DDA) [22, 23, 24, 25, 26]. The monomer-dipole DDA is free from the granularity-induced error, and its computational costs do not vary according to the aggregate’s shape. Despite these advantages, monomer-dipole DDA is subject to another error as a result of its assumption that each monomer is an electric dipole, thereby ignoring contributions from the magnetic dipole and higher-order multipoles (multipole-truncation error) [27, 28]. This multipole-truncation error persists to the electrostatic limit, as illustrated by Mackowski [28], who demonstrated that

4

the electrostatic field is highly inhomogeneous around the points of contact between neighboring monomers, especially when |m| >> 1. The multipole-truncation error in the monomer-dipole DDA can be mitigated to some extent by extending the DDA formulation from the coupled electric dipoles (CEDs) to the coupled electric and magnetic dipoles (CEMDs), as demonstrated by Mulholland et al. [29]. On the basis of the abovementioned review, the monomer-dipole DDA or the more exact superposition T-matrix method [27] is recommended for pure BC aerosols, whereas the cubic-lattice DDA must be used for non-BC materials, e.g., sulfates, organics, and water, of arbitrary shape. An open issue is that neither of these DDA schemes is suitable for modeling general atmospheric BC-containing aerosols with enormous diversities in their size, shape, and mixing state. To fill the gap between the monomer-dipole and cubic-lattice DDAs, we introduce herein a hybrid DDA that is suitable for general atmospheric BC-containing aerosols. In Section 2, we explain the theory of DDA that forms the background of the subsequent sections. In Section 3, we construct the numerical algorithm of hybrid DDA as a combination of the algorithms of monomer-dipole and cubic-lattice DDAs. In Section 4, first, we discuss the selection of polarizability prescriptions for hybrid DDA. Second, we clarify the sources of DDA error by separately evaluating the granularity-induced error and multipole-truncation error for a cluster of BC monomers. Third, we compare the hybrid DDA with conventional cubic-lattice DDA in terms of accuracy and efficiency for a cluster of BC monomers associated with non-BC materials. Finally, conclusive remarks are given in Section 5.

5

2. Theory In our DDA, we adopt either of the two theoretical methods: the CEDs or the CEMDs. Because the former is a special case of the latter, we explain here only the CEMDs. In CEMDs, the scattering target is assumed to be a collection of N subvolumes with center position ri (i = 1,…,N). Under a given incident field (Einc,i, Hinc,i) (i = 1,…,N), unknown electric and magnetic dipole moments (pi, mi) (i = 1,…,N) satisfy a system of linear equations

(

where

)

∑[

(

is Dirac’s delta and

)

)(

(

)] (

is a 3 × 3 unit tensor. Terms

)

and

are electric and

magnetic scalar polarizabilities of ith subvolume, respectively. We assume here that the materials composing the target are isotropic. The polarizability formulas depend on the physical properties of the subvolume as well as on our choice of method for approximating the singular self-integral around the subvolume [15]. The 3 × 3 dyadic

relates an electric (resp. magnetic) dipole moment to a

radiating electric (resp. magnetic) field. The 3 × 3 dyadic

relates a magnetic (resp. electric)

dipole moment to a radiating electric (resp. magnetic) field. The formulations of

and

are derived from the multipole expansion of a vector potential produced by a localized electric current [31]. After the system of equations (1) is solved with respect to (pi, mi) (i = 1,…, N), the principal quantities of interests—the differential scattering cross section dCsca/dΩ, extinction

6

cross section Cext, and the absorption cross section Cabs—are computed by the formulas reported in the literature [32]. For example, the Cabs is given by

|

|

∑ {[

(

)

](

)

[

(

)

](

)}

where k denotes the wavenumber in the surrounding medium and the superscript * denotes the complex conjugate. The corresponding CED formulations can be obtained by omitting the magnetic dipole moment m and related terms from the CEMD formulations.

3. Numerical Methods The numerical algorithms for the CEMDs are explained here. In the case of the CEDs, the three symbols 6N, 6Nm, and 6Nc should be read as 3N, 3Nm, and 3Nc, respectively.

3.1. Conventional DDAs To obtain a numerical solution of the complex linear system (1) of size 6N × 6N, iterative methods are generally preferred over direct methods for large values of N. An iterative method for a general complex matrix typically involves two matrix–vector products in each iterative cycle. In monomer-dipole DDA, each monomer-dipole can be placed at an arbitrary location in 3D space

and the resulting matrix does not possess any

attractive structure in terms of numerical algorithm. In this case, an operation of matrix–vector product requires

computation time and

storage.

7

In cubic-lattice DDA, the centroid position r of every subvolume is restricted to lie on a cubic lattice inside a rectangular cuboid sites

, a discrete space with

lattice

, where nx, ny, and nz denote the number of grid points along the x-,

y-, and z-axes of the Cartesian coordinate system, respectively (see Figure 1 for an example of rectangular cuboid

illustrated in a particle cross-section perpendicular to a Cartesian

coordinate axis). In cubic-lattice DDA, the matrix consisting of non-diagonal blocks, i.e., the second term inside the square brackets in equation (1), takes the form of a multilevel block-Toeplitz (MBT) matrix [18]. The matrix–vector product for an MBT matrix is equivalent to a discrete convolution that can be accelerated using the FFT [19, 20]. As this fact is exploited in cubic-lattice DDA, an operation of matrix–vector product requires only computation time and

storage. We use the novel method of Barrowes et al. [20]

for efficient implementation of the fast algorithm in cubic-lattice DDA.

3.2. Hybrid DDA In the hybrid DDA, the monomer-dipole DDA is applied to a cluster of BC monomers, whereas the cubic-lattice DDA is utilized for non-BC coating materials (Figure 1). To illustrate the numerical algorithm, we rewrite the matrix equation (1) in a compact form, ̃ ̃ ̃

, ̃

(

rewritten in block-partitioned form: ̃ ̃

, where

) , and so on. After an

appropriate reordering of the subvolume index, the matrix equation ̃ ̃

(

̃

̃ ̃ )( ) ̃ ̃

8

(

̃ ̃

)

̃

can be

where the vectors with superscripts m and c denote the properties of the monomers and coating materials, respectively. For instance, the ̃

vector has 6Nm complex elements, where Nm is

the number of BC monomers. The operators with superscripts mc and cm both describe the dipolar interactions between coating volumes and BC monomers. The operators with superscripts mm and cc describe the dipolar interactions within BC monomers and those within coating volumes, respectively. For example, ̃

is a rectangular matrix of size of 6Nm × 6Nc,

where Nc denotes the number of coating subvolumes. Other symbols in equation (3) should be interpreted in the same manner. In an iterative method for solving equation (3), we repeatedly update the vector ̃ to ̃ according to the matrix–vector product (

̃ ) ̃

(

̃ ̃

̃ ̃ ) ( ) ̃ ̃

In operation (4), the partial matrix–vector product ̃

̃ is computed using the FFT-based

algorithm utilizing the MBT structure of ̃ . The other three partial matrix–vector products in (4) are directly manipulated. Thus, operation (4) requires, in total, computation time and

storage.

By definition, the number of lattice sites NCL is almost always larger than the number of coating subvolumes Nc. Moreover, for the majority of atmospheric BC-containing aerosols, a large number of coating subvolumes Nc is required such that Nc >> Nm to mimic an actual shape of the coating materials. Under the situations of NCL > Nc >> Nm, we can estimate the computation time and storage requirements for the matrix–vector product (4) to be ~

and

, respectively.

9

Cubic-lattice DDA

BC-containing aerosol BC-monomer

Discritization of particle volume

3

Rectangular cuboid R CL

Hybrid DDA (= monomer-dipole DDA + cubic-lattice DDA) non-BC coating material

3

Rectangular cuboid R CL

Figure 1. A cross-sectional view of a BC-containing aerosol before (left) and after (right) a discretization in the conventional cubic-lattice DDA (top right) and in the new hybrid DDA (bottom right).

4. Results and Discussion In our numerical experiments, the complex refractive index of BC is always assumed to be 1.75 + 0.44i, a value given in the database of the Optical Properties of Aerosols and Clouds (OPAC) [30] for λ = 550 nm.

4.1. Polarizability of DDA subvolume We compare several polarizability formulas for the DDA subvolume that have been proposed on the basis of different ideas and discuss which is the most appropriate for application to our hybrid DDA. First, Draine [14] derived the

formula of Clausius–Mossotti with radiation

10

correction (CMR), where the electrostatic polarizability of a spherical subregion in a continuous material is corrected by incorporating the reaction of radiation from an oscillating electric dipole. In the CMR method, the

for each DDA subvolume is equivalent to the corrected

electrostatic polarizability of a spherical subregion of equal volume. Second, Peltoniemi [33] started from the volume integral form of Maxwell equations and then derived an

formula

for each DDA subvolume on the basis of an analytical approximation of singular self-integral over a spherical subregion in a continuous material; the resulting formula is hereafter abbreviated as PEL. In the PEL approach, the spherical subregion is assumed to be identical to the volume-equivalent sphere of the DDA subvolume. Third, we consider an the a1 term of Mie coefficients [29, 34, 35] (A1T). In this approach, the

formula using

is defined such that

the radiation fields from an electric dipole moment p are identical to the first transverse-magnetic (TM) mode radiated from the volume-equivalent sphere of a DDA subvolume excited by a plane electric field of incidence. In contrast to the CMR and PEL methods, the A1T method does not assume that the DDA subvolume is embedded in a continuous material. Fourth, Mulholland et al. [29] extended the idea of A1T to formulate both and

using the a1 and b1 terms of Mie coefficients, respectively (A1B1T). In this

approach, the radiation fields from a magnetic dipole moment m are identical to the first transverse-magnetic (TE) mode radiated from the volume-equivalent sphere of a DDA subvolume excited by a plane magnetic field of incidence. The

defined by A1B1T does not

vanish even for non-magnetic materials because an eddy current excited by an incident magnetic field also contributes to the magnetic dipole moment of each subvolume. The A1B1T implicitly assumes the use of CEMDs.

11

The most fundamental assumption in DDA is that the size of every subvolume must be sufficiently small compared to the length scale of field variations such that any macroscopic responses of each subvolume can be approximated well by those of an electric dipole. Typically, the subvolume size must be less than 1/10th the wavelength in the surrounding medium. In the cubic-lattice DDA, the subvolume size can be controlled to be small according to the choice of lattice spacing. In the monomer-dipole DDA, by contrast, the fundamental assumption could be violated because monomer size is not always much smaller than the wavelength. Therefore, we recommend for the hybrid DDA choosing a polarizability formula that is less erroneous under the monomer-dipole assumption over a practical range of monomer sizes. Accuracy tests for the four polarizability formulas in monomer-dipole DDA were performed for a single BC monomer; the results are shown in Figure 2. As expected, the DDA solution always converges to an exact Mie solution as the monomer size parameter xpp approaches zero, irrespective of the polarizability formula. The results begin to diverge when xpp > 0.05 because of differences in the size-dependent artifact inherent to each polarizability formula. As shown in Figure 2, the CMR and PEL formulas always result in larger DDA errors compared to that obtained when using the A1T formula, possibly because the monomer-dipole DDA violates the basic assumption common to CME and PEL, i.e., the DDA subvolume is assumed to be embedded inside the continuous material. Furthermore, the A1B1T performs substantially better than the A1T when xpp > ~0.1 because the contributions of the magnetic dipole relative to the electric dipole are more noticeable for larger xpp.

12

On the basis of these results, for application of the hybrid DDA, we recommend the A1B1T formula with CEMDs, whereas the A1T formula might be chosen for use with CEDs under a strict limitation of computational costs.

(a)

0.20

CMR PEL A1T A1B1T

0.1

maximum of |S11/S11,exact -1|

Qabs/Qabs,exact -1

0.2

0.0

-0.1

-0.2 0.0

0.1

0.2

0.3

0.4

0.5

(b)

0.15

CMR PEL A1T A1B1T

0.10

0.05

0.00 0.0

Monomer-size parameter xpp

0.1 0.2 0.3 0.4 Monomer-size parameter xpp

0.5

Figure 2. Monomer-dipole DDA results for a single BC monomer: (a) Fractional error in the absorption efficiency factor Qab plotted versus monomer size parameter xpp. (b) Same as (a) but showing the maximum fractional error in the scattering matrix element S11 over the entire scattering range 0 ≤ θ ≤ 180°.

4.2. Experimental conditions Our Fortran 2003 codes implementing the cubic-lattice DDA and hybrid DDA were designed to be identical except for the algorithms for volume discretization (see Figure 1) and the matrix– vector product. We chose the biconjugate gradient stabilized (BICGSTAB) method with Jacobi preconditioning [39] as an iterative method for solving the matrix equation (1). The FFTW 3.3.4 package [40] was used for FFT and inverse FFT operations for the fast algorithm of the matrix– vector product in the cubic-lattice and hybrid DDAs. The convergence criterion in BICGSTAB

13

was set to ‖ ̃ ̃

̃

‖⁄‖ ̃

‖ < 10−6. All the numerical experiments were performed on an

Apple Mac Pro with a 3.5 GHz 6-core Intel Xeon E5 processor in single-thread mode, using the gcc-4.9.3 GFortran compiler. We performed numerical experiments of DDAs for the model BC-containing aerosol target shown in Figure 3. This target consists of a collection of spheres without surface cuts; its geometric parameters are detailed in Table 1. The radius and center position of every sphere were scaled by the monomer size parameter xpp specified in each experiment. Near-exact reference solutions were obtained using the multisphere T-matrix code (MSTM version 3.0) [36, 37], with the multipole expansion order L set to L = 12. In both the cubic-lattice DDA and hybrid DDA, two DDA runs assuming the CEDs (with A1T) and CEMDs (with A1B1T) were performed separately. The lattice spacing d for both the cubic-lattice DDA and the hybrid DDA was always fixed to d = 0.2Dpp, where Dpp is the monomer diameter. In repetitive DDA runs for atmospheric BC-containing aerosols on a desktop workstation, the lattice spacing of d ≈ 0.2Dpp likely approaches the lowest acceptable value in terms of computational costs.

6 4

z

2 0 -2 -4 -6 -6 -4 -2 0

x

2

4

6

Figure 3. A model BC-containing aerosol used for numerical experiments of our DDAs. In this figure, the actual 3D geometry of the particle has been projected onto the x–z plane of the

14

Cartesian coordinate system. The cluster of 21 BC monomers (black filled circles) is referred to in the text as pure BC aerosol. The cluster of 21 BC-monomers (black filled circles) wherein seven monomers are embedded in a coating sphere (gray open circle) is referred to as BC-containing aerosol in the text. The complex refractive index of the coating material is 1.5 + 0.0i. Geometric parameters for this model particle are shown in Table 1.

Table 1. Geometric parameters for the model BC-containing aerosol shown in Fig. 3 Sphere

Radius

Center position (x, y, z)

BC-monomer 1 BC-monomer 2 BC-monomer 3 BC-monomer 4 BC-monomer 5 BC-monomer 6 BC-monomer 7 BC-monomer 8 BC-monomer 9

1 1 1 1 1 1 1 1 1

−1.3051 0.6949 −0.3327 −0.8503 0.57784 1.4625 −0.2472 2.3849 4.3849

0.2203 0.2203 −1.0497 0.6526 1.9444 −0.7680 −1.2198 3.9102 3.9102

−0.7628 −0.7628 0.4361 1.3529 0.2448 1.2761 −1.7844 2.9271 2.9271

BC-monomer 10 BC-monomer 11 BC-monomer 12 BC-monomer 13 BC-monomer 14 BC-monomer 15 BC-monomer 16 BC-monomer 17 BC-monomer 18 BC-monomer 19 BC-monomer 20 BC-monomer 21

1 1 1 1 1 1 1 1 1 1 1 1

3.3572 2.8396 4.2678 5.1525 3.4427 −1.3051 0.6949 −0.3327 −0.8503 0.5778 1.4625 −0.2472

2.6403 4.3425 5.6344 2.9219 2.4701 0.2203 0.22026 −1.0497 0.6526 1.9444 −0.7680 −1.2198

4.1261 5.0428 3.9348 4.9660 1.9056 −7.1548 −7.1548 −5.9559 −5.0391 −6.1472 −5.1159 −8.1764

Coating sphere

3.196

0

0

0

4.3. Numerical results for pure BC aerosol We first present the cubic-lattice DDA and monomer-dipole DDA results for the pure BC aerosol (i.e., a cluster of 21 BC-monomers, see Fig. 3). Figures 4a–c show fractional errors in the DDA results as a function of monomer size parameter xpp. The cubic-lattice DDA results show positive bias in both absorption Qabs and scattering Qsca, and this bias is persistent over the

15

entire investigated xpp range as a result of the granularity-induced error. In the case of the cubic-lattice DDA, the contributions of magnetic dipoles to observable quantities are negligible compared to those of electric dipoles because the size parameter of each subvolume (~0.2xpp) is much less than unity. Therefore, the CED and CEMD results are almost identical. By contrast, in monomer-dipole DDAs, the contributions of magnetic dipoles to observable quantities could be important because the size parameter of each subvolume—specifically, the monomer size parameter xpp—is not always substantially less than unity. Thus, the CEMD result remains reasonably stable, whereas the corresponding CED result diverges for larger xpp values. In Figures 4a–c, the magnitude of fractional error in monomer-dipole CEMDs is almost always less than that in cubic-lattice CEMDs, demonstrating that the use of CEMDs in monomer-dipole DDA successfully reduces its multipole-truncation error to almost always be less than the discretization-induced error in cubic-lattice DDA. In addition to the DDA results, Figures 4a–c also present the corresponding fractional error in MSTM with multipole expansion order L = 1 (MSTM (L = 1)). The good agreements between the results of monomer-dipole CEMDs and MSTM (L = 1) were expected on the basis of their similar theoretical assumptions. The slight difference between their results might stem from the assumption of plane-wave excitation used in formulating the A1B1T polarizability. To clarify the nature of the multipole-truncation error, we observed the convergence behavior of MSTM with respect to the multipole expansion order L at a particular value of monomer size parameter xpp (Figure 4d). In Figure 4d, we show the results at monomer size parameter xpp = 0.4. The corresponding plots of the results at other xpp values (0.01 ≤ xpp ≤ 0.4)

16

were similar to those presented in Figure 4d and were omitted. Figure 4d shows that the MSTM (L) results rapidly approach exact solutions as L increases: the magnitude of fractional error in MSTM (1) diminishes by approximately half in MSTM (2). On the basis of the results in Figures 4a–d, we attribute the multipole-truncation error in monomer-dipole CEMDs to an artifact of neglecting higher-order multipoles, mostly L = 2. Thus, we expect a further reduction of the multipole-truncation error by incorporating the effects of the electric quadrupole (second-order TM mode) into our DDA formulations. Lemaire [38] proposed a coupled multipole method that includes up to the electric quadrupole by extending the coupled dipole method in DDA. Unfortunately, the advantages of the coupled multipole method over its coupled dipole counterpart were not clearly illustrated herein [38] because his numerical experiments were performed for a sphere discretized using a cubic-lattice, wherein the granularity-induced error almost always dominates the total DDA error.

17

cubic-lattice CED cubic-lattice CEMD monomer-dipole CED monomer-dipole CEMD MSTM with L=1

(a)

0.1

0.0

0.2

Qsca/Q sca,exact -1

Qabs/Q abs,exact -1

0.2

-0.1

-0.2 0.1 0.2 0.3 0.4 Monomer-size parameter xpp

0.0

-0.1

0.5

0.20

0.1

0.2

0.3

0.4

0.5

Monomer-size parameter xpp

cubic-lattice CED cubic-lattice CEMD monomer-dipole CED monomer-dipole CEMD MSTM with L=1

(c)

0.0

0.15 0.10 0.05 0.00

0.06 Fractional error in MSTM result

maximum of |S11/S11,exact -1|

0.1

-0.2 0.0

0.25

cubic-lattice CED cubic-lattice CEMD monomer-dipole CED monomer-dipole CEMD MSTM with L=1

(b)

(d)

0.04

xpp = 0.4 Q abs Q sca S11 (absolute max.)

0.02 0.00 -0.02 -0.04 -0.06

0.0

0.1

0.2

0.3

0.4

0.5

0

Monomer-size parameter xpp

2 4 6 8 10 Multipole expansion order L

12

Figure 4. Results of the cubic-lattice DDA, monomer-dipole DDA, and MSTM for pure BC aerosol. (a) Fractional error in the absorption efficiency factor Qabs plotted versus the monomer size parameter xpp. (b) Same as (a) but showing the fractional error in the scattering efficiency factor Qsca. (c) Same as (a) and (b) but showing the maximum fractional error in the scattering matrix element S11 over the whole scattering angle 0 ≤ θ ≤ 180° at the azimuth angle of 45°. (d) Fractional error in the MSTM (L) results versus the multipole expansion order L for monomer size parameter xpp = 0.4.

4.4. Numerical results for BC-containing aerosol Finally, we present our cubic-lattice DDA and the hybrid DDA results for the BC-containing aerosol—a cluster of 21 BC monomers partially embedded inside a coating sphere (Figure 3,

18

Table 1). In the cubic-lattice DDA, the CED and CEMD results are almost identical for the same reason explained in Section 4.3. The cubic-lattice DDA suffers from positive bias in Qabs as a result of the granularity-induced error, mainly in the 14 BC monomers outside the coating sphere. In the hybrid DDA, the CEMD is considerably less erroneous than the CED for larger monomer sizes (xpp > 0.2) as a result of mitigating the multipole-truncation error in BC monomers by incorporating the magnetic dipole. The maximum fractional error in the S11 is sometimes greater in the hybrid CEMD compared to that in the cubic-lattice CEMD (Figure 5c). As shown in Figure 5d, the maximum fractional error in the scattering matrix element S11 occurs around the backward-scattering direction (170° < θ < 180°), where the magnitude of S11 approaches the minimum. Figure 5d shows that the S11 results for the hybrid CEMD and the cubic-lattice CEMD are comparably accurate over the whole scattering angle range, except for the range near the backward direction.

19

(a)

cubic-lattice CED cubic-lattice CEMD hybrid CED hybrid CEMD

0.05

0.00

-0.05

-0.10 0.0

0.1

0.2

0.3

0.4

0.10

(b)

0.00

-0.05

-0.10 0.0

0.5

(d)

cubic-lattice CED cubic-lattice CEMD hybrid CED hybrid CEMD

0.20

0.1

0.2

0.3

0.4

0.5

cubic-lattice CEMD hybrid CEMD exact (MSTM)

xpp = 0.49 0.01

S11

maximum of |S11/S11,exact -1|

(c)

0.25

0.1

Monomer-size parameter xpp

Monomer-size parameter xpp 0.30

cubic-lattice CED cubic-lattice CEMD hybrid CED hybrid CEMD

0.05

Qsca/Qsca,exact -1

Qabs/Qabs,exact -1

0.10

0.15 0.10

xpp = 0.25

0.001 0.05 0.00 0.0

0.0001 0.1

0.2

0.3

0.4

0.5

0

Monomer-size parameter xpp

40 80 120 160 Scattering angle q (degree)

Figure 5. Cubic-lattice DDA and hybrid DDA results for the model BC-containing aerosol. (a) Fractional error in the absorption efficiency factor Qabs plotted versus the monomer size parameter xpp. (b) Same as (a) but showing the fractional error in the scattering efficiency factor Qsca. (c) Same as (a) and (b) but showing the maximum fractional error in the scattering matrix element S11 over the whole scattering angle range 0 ≤ θ ≤ 180° at the azimuth angle of 45°. (d) The S11 versus scattering angle θ at the azimuth angle of 45°, evaluated at two different xpp values.

Table 2 shows a comparison of the computational costs of the four different DDAs for the BC-containing aerosol. The CPU time denotes the time required to solve matrix equation (1) once. The hybrid DDA is several times more efficient compared to the cubic-lattice DDA in

20

terms of both CPU time and memory. As in this example, the hybrid DDA is substantially more efficient than the conventional cubic-lattice DDA when the size of the rectangular cuboid is substantially smaller in the former DDA (see Figure 1).

Table 2. Computational costs in different DDAs for the BC-containing aerosol. DDA methods

Number of dipoles (N)

FFT length1 ~ f 2 2M NCL

CPU time (sec)

Memory (GB)

Cubic-lattice CED

2914

1,389,555

28

0.2

Cubic-lattice CEMD

2914

5,558,220

104

0.8

Hybrid CED

21+1720

268,119

3.8

0.07

Hybrid CEMD

21+1720

1,072,476

13

0.2

1

Memory requirement and CPU time are approximately proportional to the FFT length. Parameter f

denotes the dimension of smallest building block of the DDA matrix (f =3 for CED, f = 6 for CEMD), M is the dimension of space (= 3), NCL is the number of lattice sites in the

.

5. Concluding Remarks We proposed a novel computational framework of DDA that largely resolves the drawbacks of conventional DDA schemes in applications involving atmospheric BC-containing aerosols in an arbitrary mixing state. In this framework, the monomer-dipole DDA is applied to a cluster of BC monomers, whereas the efficient cubic-lattice DDA is applied to non-BC coating materials. Thus we refer to this framework as hybrid DDA. The granularity-induced error inherent to the conventional cubic-lattice DDA is suppressed in the hybrid DDA by avoiding the discretization of BC monomers. The multipole-truncation error inherent to the monomer-dipole DDA is successfully mitigated through the use of the CEMDs instead of the CEDs. Owing to the flexibility in the hybrid DDA of choosing the size of the rectangular cuboid

21

, a severe

reduction of computational efficiencies in the cubic-lattice DDA, depending on the particle morphology, can be avoided in the hybrid DDA. In atmospheric modeling, the hybrid DDA is useful for quantitative investigations of the optical properties of atmospheric BC-containing aerosols with enormous variations in size, shape, and state of mixing with other compounds such as sulfates, organics, and water. With respect to aerosol measurements, the hybrid DDA is a new choice of light scattering simulator for designing or evaluating optical particle characterization techniques for morphologically complex BC-containing aerosols, as well as a cluster of metallic nanospheres associated with other dielectric compounds. Finally, we note that the multipole-truncation error in the hybrid DDA would be suppressed further by an extension of DDA to the coupled multipole method.

Acknowledgments I acknowledge a suggestion by M. Yurkin for analyzing the multipole-truncation error in monomer-dipole DDA using a superposition T-matrix method. I thank D. W. Mackowski for assisting my interpretations of the MSTM results. This work was supported by JSPS KAKENHI Grant No. 15H05465 and the global environment research fund of the Ministry of the Environment, Japan (2-1403). The DDA codes used herein are available upon request to the author.

References [1] M. Z. Jacobson, Strong radiative heating due to the mixing state of black carbon in atmospheric aerosols, Nature 409 (2001) 695-697.

22

[2] V. Ramanathan, G. Carmichael, Global and regional climate changes due to black carbon, Nature Geoscience 1 (2008) 221-227.

[3] T. C. Bond, S. J. Doherty, D. Fahey, P. Forster, T. Berntsen, B. DeAngelo, M. Flanner, S. Ghan, B. Kärcher, D. Koch, et al., Bounding the role of black carbon in the climate system: A scientific assessment, Journal of Geophysical Research: Atmospheres 118 (2013) 5380-5552.

[4] K. Adachi, P. R. Buseck, Internally mixed soot, sulfates, and organic matter in aerosol particles from Mexico City, Atmospheric Chemistry and Physics Discussions 8 (2008) 9179-9207.

[5] K. A. Fuller, W. C. Malm, S. M. Kreidenweis, Effects of mixing on extinction by carbonaceous particles, Journal of Geophysical Research 104 (1999) 15941-15954.

[6] M. Kahnert, T. Nousiainen, H. Lindqvist, M. Ebert, Optical properties of light absorbing carbon aggregates mixed with sulfate: assessment of different model geometries for climate forcing calculations, Optics Express 20 (2012) 10042-10058.

[7] M. Kahnert, T. Nousiainen, H. Lindqvist, Models for integrated and differential scattering optical properties of encapsulated light absorbing carbon aggregates, Optics Express 21 (2013) 7974-7993.

[8] C. He, K.-N. Liou, Y. Takano, R. Zhang, M. Levy Zamora, P. Yang, Q. Li, L. R. Leung, Variation of the radiative properties during black carbon aging: theoretical and experimental intercomparison, Atmospheric Chemistry and Physics 15 (2015) 11967-11980.

[9] R. Gao, J. Schwarz, K. Kelly, D. Fahey, L. Watts, T. Thompson, J. Spackman, J. Slowik, E. Cross, J.-H. Han, et al., A novel method for estimating light-scattering properties of soot aerosols using a modified single-particle soot photometer, Aerosol Science and Technology 41 (2007) 125-135.

23

[10] N. Moteki, Y. Kondo, K. Adachi, Identification by single-particle soot photometer of black carbon particles attached to other particles: Laboratory experiments and ground observations in Tokyo, Journal of Geophysical Research: Atmospheres 119 (2014) 1031-1043.

[11] J. Schwarz, A. Perring, M. Markovic, R. Gao, S. Ohata, J. Langridge, D. Law, R. McLaughlin, D. Fahey, Technique and theoretical approach for quantifying the hygroscopicity of black-carbon-containing aerosol using a single particle soot photometer, Journal of Aerosol Science 81 (2015) 110-126.

[12] T. Cheng, X. Gu, Y. Wu, H. Chen, T. Yu, The optical properties of absorbing aerosols with fractal soot aggregates: Implications for aerosol remote sensing, Journal of Quantitative Spectroscopy & Radiative Transfer 125 (2013) 93-104.

[13] E. M. Purcell, C. R. Pennypacker, Scattering and absorption of light by nonspherical dielectric grains, The Astrophysical Journal 186 (1973) 705-714.

[14] B. T. Draine, The discrete-dipole approximation and its application to interstellar graphite grains, The Astrophysical Journal 333 (1988) 848-872.

[15] M. A. Yurkin, A. G. Hoekstra, The discrete dipole approximation: an overview and recent developments, Journal of Quantitative Spectroscopy & Radiative Transfer 106 (2007) 558-589.

[16] K. Adachi, S. H. Chung, P. R. Buseck, Shapes of soot aerosol particles and implications for their effects on climate, Journal of Geophysical Research 115 (2010) D15206.

[17] B. Scarnato, S. Vahidinia, D. Richard, T. Kirchstetter, Effects of internal mixing and aggregate morphology on optical properties of black carbon using a discrete dipole approximation model, Atmospheric Chemistry and Physics 13 (2013) 5089-5101.

[18] P. J. Flatau, B. T. Draine, G. L. Stephens, Light scattering by rectangular solids in the discrete-dipole approximation: a new algorithm exploiting the block-Toeplitz structure, Journal of Optical Society of America A 7 (1990) 593-600.

24

[19] J. J. Goodman, B. T. Draine, P. J. Flatau, Application of fast-Fourier transform techniques to the discrete-dipole approximation, Optics Letters 16 (1991) 1198-1200.

[20] B. E. Barrowes, F. L. Teixeira, J. A. Kong, Fast algorithm for matrix-vector multiply of asymmetric multilevel block-toeplitz matrices in 3-D scattering, Microwave and Optical Technology Letters 31 (2001) 28-32.

[21] A. Rahmani, P. C. Chaumet, G. W. Bryant, On the importance of local-field corrections for polarizable particles on a finite lattice: application to the discrete dipole approximation, The Astrophysical Journal 607 (2004) 873-878.

[22] A. R. Jones, Electromagnetic wave scattering by assemblies of particles in the Rayleigh approximation, Proceedings of the Royal Society of London A 366 (1979) 111-127.

[23] A. R. Jones, Scattering efficiency factors for agglomerates for small spheres, Journal of Physics D: Applied Physics 12 (1979) 1661-1672.

[24] M. F. Iskander, H. Y. Chen, J. E. Penner, Optical scattering and absorption by branched chains of aerosols, Applied Optics 28 (1989) 3083-3091.

[25] T. Kozasa, J. Blum, T. Mukai, Optical properties of dust aggregates I. Wavelength dependence, Astronomy and Astrophysics 263 (1992) 423-432.

[26] T. Kozasa, J. Blum, H. Okamoto, T. Mukai, Optical properties of dust aggregates II. Angular dependence of scattered light, Astronomy and Astrophysics 276 (1993) 278-288.

[27] D. W. Mackowski, Calculation of total cross sections of multiple-sphere clusters, Journal of Optical Society of America A 11 (1994) 2851-2861.

[28] D. W. Mackowski, Electrostatics analysis of radiative absorption by sphere clusters in the Rayleigh limit: application to soot particles, Applied Optics 34 (1995) 3535-3545.

25

[29] G. W. Mullholland, C. F. Bohren, K. A. Fuller, Light scattering by agglomerates: coupled electric and magnetic dipole method, Langmuir 10 (1994) 2533-2546.

[30] M. Hess, P. Koepke, I. Schult, Optical properties of aerosols and clouds: The software package OPAC, Bulletin of the American Meteorological Society 79 (1998) 831-844.

[31] R. K. Wangsness, Electromagnetic Fields 2nd edition, John Wiley & Sons Hoboken, NJ (1986)

[32] P. C. Chaumet, A. Rahmani, Coupled-dipole method for magnetic and negative-refraction materials, Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 22-29.

[33] J. I. Peltoniemi, Variational volume integral equation method for electromagnetic scattering by irregular grains, Journal of Quantitative Spectroscopy & Radiative Transfer 55 (1996) 637-647.

[34] W. T. Doyle, Optical properties of a suspension of metal spheres, Physical Review B 39 (1989) 9852-9858.

[35] H. Okamoto, Light scattering by clusters: the a1-term method, Optical Review 2 (1995) 407-412.

[36] D. W. Mackowski, M. I. Mishchenko, A multiple sphere T-matrix Fortran code for use on parallel computer clusters, Journal of Quantitative Spectroscopy & Radiative Transfer 112 (2011) 2182–2192.

[37] D. W. Mackowski, A general superposition solution for electromagnetic scattering by multiple spherical domains of optically active media, Journal of Quantitative Spectroscopy & Radiative Transfer 133 (2014) 264-270.

[38] T. Lemaire, Coupled-multipole formulation for the treatment of electromagnetic scattering by a small dielectric particle of arbitrary shape, Journal of Optical Society of America A 14 (1997) 470-474.

26

[39] R. Barrett, M. W. Berry, T. F. Chan, J. Demmel, J. M. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine, H. Van der Vorst, Templates for the solution of linear systems: building blocks for iterative methods, SIAM, Philadelphia PA, (1994).

[40] M. Frigo, S. G. Johnson, The design and implementation of fftw3, Proceedings of the IEEE 93 (2005) 216-231.

Highlights 

Light scattering code is developed for atmospheric black carbon aerosols.



A hybrid discretization scheme in discrete dipole approximation (DDA) is proposed.



The addition of the magnetic dipole improves the accuracy of the hybrid DDA.



Improved computational efficiencies are achieved in the hybrid DDA.

27