Open Access
Add to BibSonomyAbstract
A T-matrix method for scattering by particles with small-scale surface roughness is presented. The method combines group theory with a perturbation expansion approach. Group theory is found to reduce CPU-time by 4–6 orders of magnitude. The perturbation expansion extends the range of size parameters by a factor of 5 compared to non-perturbative methods. An application to optically hard particles shows that small-scale surface roughness changes scattering in side- and backscattering directions, and it impacts the single-scattering albedo. This can have important implications for interpreting remote sensing observations, and for the climate impact of mineral aerosols.
©2011 Optical Society of America
1. Introduction
Small-scale surface roughness is a morphological property that is encountered in many types of aerosols in planetary atmospheres, as well as in mineral particles in the interplanetary and interstellar medium. Modeling the optical properties of such particles is of high relevance in remote sensing, astrophysics, and in quantifying the radiative climate forcing effect of mineral dust. However, small-scale surface roughness is also among the most challenging morphological features in numerical electromagnetic scattering computations. For instance, mineral aerosols in planetary atmospheres typically have large size parameters in the visible part of the spectrum. (The size parameter is defined as x = 2πr/λ, where r is the particle radius, and λ is the wavelength of light.) However, geometric optics, which is an approximate method valid for large size parameters, cannot be applied to such particles owing to the small size-scale of the surface perturbations. On the other hand, computational methods based on rigorous electromagnetic theory are typically plagued by ill-conditioning problems and rapidly growing computation time for increasing size parameters.
Previous computational studies of particles with small-scale surface roughness have therefore been limited to rather moderate size parameters. For instance, a recent modeling study of hematite aerosols at a wavelength of λ = 633 nm considered Chebyshev particles up to r = 1.4μm [1], which corresponds to a size parameter of 14. In the terrestrial atmosphere, aerosols in the coarse mode (i.e. with radii of 1–5 μm) typically make the dominant contribution to the mass concentration of mineral dust aerosols under background conditions, while during dust storm events considerably larger particles can be suspended in air. Thus, current computational methods for particles with small-scale surface roughness are severely limited in the range of size parameters for which numerical computations are sufficiently stable and expedient. Also, to the best of our knowledge, all previous studies based on T-matrix methods have been limited to model particles with axisymmetric geometries (e.g. [1–3])
Our current understanding, although based on relatively few investigations, suggests a potentially high relevance of surface roughness for the optical properties of dielectric particles. For instance, it has been demonstrated that the phase function of high-order Chebyshev particles with a small perturbation amplitude can substantially differ from that of a size-equivalent homogeneous sphere [2]. At higher values of the imaginary part of the refractive index, differences between spheres and spheroids disappear, while differences between spheres and high-order Chebyshev particles become more pronounced [2]. A plausible explanation is that internal resonances inside the particle, which are mainly responsible for the differences between the phase functions of spheres and spheroids, are quenched in more strongly absorbing particles. On the other hand, the impact of small scale surface roughness in high-order Chebyshev particles is not reduced by high absorption inside the particle.
It has also been noted that the impact of small-scale surface roughness may become more pronounced for particles with larger real parts of the refractive index [1]. However, these results are still preliminary. The conditions under which small-scale surface roughness has a dominant impact on the optical properties may be a complex interplay of several physical properties, such as the real and imaginary parts of the refractive index, the roughness amplitude, and the particle size. An important prerequisite for further investigations is to overcome the present limitations of computational methods. This is the main purpose of the present study.
Before proceeding, we will attempt to give a formal definition for small-scale surface roughness. Consider a perturbation of the particle surface with a (mean) perturbation wavelength Λ and a (mean) perturbation amplitude A that are small compared to the characteristic size r 0 of the particle and small compared to the wavelength λ of the incident light. To take a specific example, consider an axisymmetric Chebyshev particle with the surface parameterization
where ℓ is the order of the Chebyshev polynomial Tℓ(θ) = cos(ℓθ), and ε is known as the deformation parameter. The perturbation wavelength and amplitude in this case are given by In a recent study on the optical properties of Chebyshev particles [2] it has been shown that for low polynomial orders ℓ the phase function changes with ℓ. However, for high-order Cheby-shev particles, the phase function becomes independent of ℓ, even though the phase function is distinctly different from that of the unperturbed sphere. So for a surface perturbation with a sufficiently small Λ, the optical properties become independent of Λ, even though they are sensitive to the perturbation amplitude. Based on this observation (and following [1]), we define small-scale surface roughness as follows.- The roughness wavelength Λ is sufficiently small so that any further decrease in Λ does not alter the optical properties.
- Λ ≪ 2πr 0, i.e., Λ is much smaller than the circumference of the particle.
- Λ ≪ λ, where λ is the wavelength of light.
- A ≪ r 0, where A is the roughness amplitude.
- A ≪ λ.
- A is sufficiently large so that the optical properties of a particle with a perturbed boundary surface differ from those of the corresponding unperturbed geometry.
In the following section we present our numerical approach for modeling electromagnetic scattering by particles with small-scale surface roughness. The method is based on the T-matrix formulation of the scattering problem. In Sect. 3 we present some illustrative results of our numerical implementation. Concluding remarks are given in Sect. 4.
2. Methods
Numerically exact electromagnetic scattering computations for irregular particles become prohibitively time consuming for size parameters much larger than unity. For this reason, expedient numerical methods need to rely on adequate simplifications. One approach is based on abandoning rigorous electromagnetic theory, and invoking approximations to the physics of the scattering problem, such as in ray tracing methods or Rayleigh-Debye-Gans theory. Possible disadvantages of such ad hoc approximations have been discussed in [4]. Another approach is to use numerically exact methods to solve the electromagnetic scattering problem in conjunction with certain symmetry assumptions about the geometry of the scatterer. The simplest example is Lorenz-Mie theory for scattering by homogeneous spheres. A much more flexible, but also computationally more costly approach is to consider particles with axial symmetry. For instance, Chebyshev particles as defined in Eq. (1) belong to this class of geometries.
Figure 1 (left) shows an example for a Chebyshev particle. The prize we pay for restricting the geometry to axisymmetric symmetry is that we can only account for surface roughness in the polar direction, while the particle surface is unperturbed in the azimuthal direction. A more general model is obtained by perturbing the sphere by Chebyshev polynomials in both the polar and the azimuthal direction, i.e.
Such a particle is shown in Fig. 1 (right). Since axisymmetric Chebyshev particles are described by the coordinates (θ,r(θ)), we will refer to them as “2D Chebyshev particles”. Non-axisymmetric Chebyshev particles are characterized by the coordinates (θ,ϕ,r(θ,ϕ)), so we will refer to them as “3D Chebyshev particles”. 3D Chebyshev particles have a lower symmetry than 2D Chebyshev particles. On the other hand, they appear to provide a more realistic model for particles with small-scale surface roughness. So we consider them here as a compromise between the computationally efficient 2D Chebyshev model, and a realistic model that would assume a completely irregular perturbation of the surface, such as a Gaussian random sphere [5].
Fig. 1 2D (left) and 3D (right) Chebyshev particle with polynomial order ℓ=45 and deformation parameter ε =0.05.
Although not axisymmetric, 3D Chebychev particles still have a high degree of symmetry, of which we will take advantage in the computations. The way to systematically exploit symmetries in electromagnetic scattering, and in fact in all disciplines of physics, is to use group theory.
2.1. Application of group theory
Systematic accounts of the use of group theory in electromagnetic scattering theory are given, e.g., in [6, 7]. Here we only summarize the most important points. Symmetry operations are coordinate transformations that bring a particle into a new orientation indistinguishable from the original one. Therefore the optical properties are invariant under such transformations. The set of all symmetry operations of a particle forms a so-called point group. All point groups are sub-groups of the orthogonal group 𝒪(3), so the elements of such groups consist of rotations and pseudo-rotations.
In electromagnetic scattering theory we represent the elements g of a point group 𝒢 by unitary matrices U(g) that operate on the vector space of the vectorial eigensolutions of the Helmholtz equation. The invariance of the optical properties is expressed by
or where the T-matrix T contains the complete information on the particle’s scattering and absorption properties for a given wavelength, and where [A, B]=A · B – B · A is know as the commutator of the two matrices A and B. The commutation relation of the T-matrix can, in fact, be derived from very general considerations about boundary symmetries in differential and integral equation problems [8].For each point group we can identify a minimum set of generators g 1,..., gr ∈ 𝒢, from which all other group elements can be obtained by combination of the generators. Only the generators provide us with independent commutation relations for the T-matrix [7].
Consider as an example a 3D Chebyshev particle of even order ℓ. The symmetries of such a particle are described by the point group 𝒟ℓh, which contains 4ℓ elements. The generators of that group are the elements Cℓ, , and σh, where Cℓ represents a rotation about the main symmetry axis by an angle 2π/ℓ, denotes a rotation by π about an axis perpendicular to the main symmetry axis, and σh stands for a reflection in a plane perpendicular to the main symmetry axis. For this group, the commutation relations (Eq. (6)) of the three generators become in explicit form (see [7] for details)
The indices n,n ′ = 1, 2,... are related to the degree, m = −n,...,n, m ′ = –n ′,...,n ′ to the order, and τ, τ ′ = 1, 2 to the mode of the vector spherical wave functions. The first commutation relation comes from the rotational symmetry operation Cℓ, the second comes from the dihedral symmetry , and the third relation originates from the reflection symmetry σh.In numerical calculations, symmetries can be exploited in three different ways.
- The commutation relations reduce the number of non-zero, independent T-matrix elements that need to be numerically evaluated by a factor of 1/M [7], where M denotes the order of the symmetry group (i.e. the number of the elements in the group). For instance, for 3D Chebyshev particles of even order ℓ this reduces the computation time by a factor of 1/(4ℓ).
- In the extended boundary condition method [9], the T-matrix elements are computed by numerically evaluating surface integrals over cross products of vector spherical wave functions, where the integration surface is the boundary surface of the particle. It can be shown by use of the commutation relations that the integration area can be reduced by a factor of 1/M [6]. In conjunction with the reduction in the number of T-matrix elements, this results in a total reduction in computation time by a factor of 1/M 2. For even-order 3D Chebyshev particles, this would reduce the computation time by a total of 1/(4ℓ)2.
- By use of the matrices U(g) and the so-called characters of the group one can construct a transformation matrix that brings all matrices that commute with U(g) into block diagonal form, where the number of block matrices is equal to the number of irreducible representations of the group. This method has been presented and tested in [7]. The merit of the method is to save additional computation time and, most importantly, to significantly reduce numerical ill-conditioning problems. The method is completely general and can be applied to any geometry with discrete symmetries. However, for particles with small-scale surface roughness there is an even more efficient method for circumventing the notorious ill-conditioning problems in numerical T-matrix computations, which will be discussed in the following subsection.
2.2. Perturbation expansion of the T-matrix
Most approaches for computing a T-matrix, such as the extended boundary condition method [9], the separation of variables method [10], or the generalized point-matching method [11] are based on computing two matrices Q and Rg Q, from which the T-matrix is obtained according to
The origin of the ill-conditioning problems is the numerical inversion of the Q-matrix. The ill-conditioning problems tend to drastically increase with increasing size parameters. For particles that deviate only mildly from a reference geometry, such as a sphere, one can avoid the ill-conditioning problem by performing a perturbation expansion. This has been proposed and tested for the scalar Helmholtz equation in [12]. Here we will use this approach for the vector Helmholtz problem.Suppose we have a reference geometry (such as a sphere) with Q-matrix Q 0, and suppose we perform a small perturbation of the reference geometry, resulting in a new particle (e.g. a Chebyshev particle with small-scale surface roughness) with matrices Q and Rg Q. We formally define ΔQ=Q – Q 0 and substitute this into Eq. (10), which yields after rearranging terms
We subtract T · ΔQ and multiply by , which gives While Eq. (10) involves the potentially ill-conditioned inversion of the matrix Q, Eq. (12) only requires us to invert the matrix Q 0. For instance, if the unperturbed geometry is a sphere, then the matrix Q 0 is diagonal, and computation of is trivial; the ill-conditioning problem has completely disappeared! The prize we have to pay for this is that Eq. (12) only provides us with an implicit equation for the T-matrix.Equation (12) is of the same form as the Lippmann-Schwinger equation for the Stückelberg-Feynman propagator in quantum electrodynamics (e.g. [13]). In practice, one solves this type of equation by performing a perturbation expansion. To this end, we obtain a zeroth-order approximation by setting T = 0 on the rhs of Eq. (12), i.e.
A first-order approximation T (1) is obtained by substituting T (0) in the rhs of Eq. (12). This can be continued iteratively. So, more generally, if we have an approximate solution T ( n −1) of order n − 1, then we obtain a solution of order n according to Equation (14) in conjunction with the initial value given in Eq. (13) is the sought perturbative T-matrix solution. As to the convergence of the method for n → ∞, we take a pragmatic point of view by numerically testing the accuracy of the results for increasing n. To this end, we exploit the reciprocity condition [12, 14]. In general, we expect the method to be most efficient and robust for geometries that deviate only mildly from the unperturbed geometry.3. Results
We test our approach by implementing the perturbation expansion method into the Tsym program, which is a T-matrix code for scattering by 3D targets that has been specifically made for accounting for point-group symmetries [6, 7]. This code has previously been applied to polyhedral prisms only (e.g. [15]). We now added 2D and 3D Chebyshev particles of arbitrary order to the code. We will here present comparisons with results computed with mieschka, which is a comprehensively tested T-matrix code for axisymmetric particles [14].
We thus follow the traditional way of testing the accuracy of newly developed numerical methods by comparing their performance with well-established existing codes. A main motivation of the work presented here is to develop a method for particles with small-scale surface roughness that goes significantly beyond the current state-of-the-art by extending the accessible range of size parameters. We therefore expect that direct comparisons with existing codes are only possible within a limited size range; the most interesting results are those obtained for larger size parameters, which are beyond the reach of existing codes. However, it is nevertheless possible to test the method at larger size parameters. There exists a highly sensitive method for testing the accuracy of electromagnetic scattering computations, namely, the reciprocity condition [12, 14]. We will apply the reciprocity condition to test our approach for particle size parameters that lie beyond the capabilities of existing methods.
3.1. Comparison to 2D Chebyshev computations with mieschka
We compute the polarized differential scattering cross sections for a 2D Chebyshev particle of order ℓ = 45, deformation parameter ε = 0.03, and refractive index m = 3 + 0.1i. The size of the unperturbed sphere is r 0=1.4 μm, and the wavelength is λ =0.6328 μm, so the size parameter is approximately x=14. Our choice of the refractive index is typical for hematite at visible wavelengths [16]. The incident field is taken to be in the positive z-direction, the z-axis is assumed to coincide with the particle’s main rotational symmetry axis, and the scattered field is computed in the xz-plane as a function of the scattering angle Θ. Figure 2 shows Sα,β (Θ)= k 2 (dσ/dΩ)α,β (Θ), where k = 2π/λ is the wavenumber, and where (dσ/dΩ)α,β is the polarized differential scattering cross section. α = h means that the incident field is polarized “horizontally” that is in the xz plane, while α = v means that the incident field is polarized “vertically”, i.e. perpendicularly to the scattering plane. Similarly, β = h and β = v refer to the polarization state of the scattered field. Figures 2a and 2b show (dσ/dΩ)h,h and (dσ/dΩ)v,v, respectively. Results obtained with Tsym in conjunction with the perturbation expansion approach are represented by a black line, while the results computed with mieschka are plotted in red. The perturbation expansion of the T-matrix has been carried out to third order. The Tsym and mieschka results are indistinguishable. The cross-polarization components (dσ/dΩ)h,v and (dσ/dΩ)v,h (not shown) are essentially zero in this case.
Fig. 2 Top row: hh (left) and vv (right) components of the polarized differential scattering cross sections of 2D Chebyshev particles in a fixed orientation, computed with the group theoretical/perturbative T-matrix approach (black) and with mieschka (red). Bottom row: hh (left) and vv (right) components computed with mieschka for two different orientations.
3.2. Reciprocity condition
In general, if the wavevector of the incident field points in the direction k̂ inc, and that of the scattered field in the direction k̂ sca, then the reciprocity condition states
For instance, if we choose k̂ inc=ẑ and k̂ sca=x̂ (i.e. Θ = 90°), then in the reciprocal case we need to take the incident field in the direction −x̂ and the scattered field in the direction −ẑ. Equivalently, we can keep the direction of the incident field fixed, rotate the particle by an angle θp = 90° around the y-axis, and take the scattered field in the direction −x̂, i.e. Θ = 270°. So the reciprocity condition in this example becomes We usually require that the reciprocity condition has to be satisfied with an error smaller than 2 % [14]. This has proven to be a highly sensitive necessary condition for the accuracy of T-matrix computations.Table 1 shows Sα,β (Θ,θp) computed with Tsym, where k denotes the wavenumber. Since the differential scattering cross section has units μm2 sr−1, Sα,β has units sr−1. We see that for both polarization components the error is less than 2 %, so the reciprocity condition is fulfilled with the required accuracy. In the following subsection we will present more reciprocity tests for an extended range of size parameters and for the case of 3D Chebyshev particles. But we first want to mention an interesting observation.
Table 1. Reciprocity Test for 2D Chebyshev Particles
Figure 2c shows a comparison of Sh,h(Θ;θp = 0°) (black) and Sh,h(Θ;θp = 90°) (red) computed with mieschka. Figure 2d shows a corresponding comparison for Sv,v. We see that the differential scattering cross sections for the two particle orientations are very similar, even though a 2D Chebyshev particle does not possess spherical symmetry. In fact, Sα,β averaged over particle orientations (not shown) is rather similar to that of particles in a fixed orientation. This fact may be exploited in simplifying orientational averaging in numerical computations.
3.3. Illustrative application to 3D Chebyshev particles
We performed computations for 3D Chebyshev particles at an optical wavelength of λ =0.6328 μm, assuming a refractive index of m = 3+0.1i, and considering particles sizes r 0 = 1, 2,...,7 μm. Thus the range of size parameters now extends up to x max=70. We choose a size-dependent Chebyshev order ℓ such that the perturbation wavelength Λ = 2πr 0 /ℓ is fixed at Λ = λ/4. This means that for r 0 = 1, 2,...,7 μm we use ℓ = 40, 80,...,280, respectively. For Λ ≲ λ/4, the optical properties do no longer depend on Λ, which was one of the essential characteristics in our definition of small-scale surface roughness. We experiment with two different cases for the deformation parameter. In the first case, we use a constant value of the relative amplitude ε = 0.01, in the second case we use a constant value of the absolute amplitude A = εr 0 = 0.11λ. For r = 7 μm, these two cases coincide. In all cases we carried out the perturbation expansion of the T-matrix to sixth order. Tables 2 and 3 show the reciprocity tests for the two choices of the perturbation amplitudes. In either case, the reciprocity condition is satisfied with high accuracy.
Table 2. Reciprocity Test for 3D Chebyshev Particles with ε = 0.01
Table 3. Reciprocity Test for 3D Chebyshev Particles with A = 0.11λ
Figure 3 shows the single scattering albedo ω (top left), the asymmetry parameter g (top right), and the backscattering cross section C bak (bottom left) as a function of particle size after averaging over particle orientations and polarization states. For comparison, corresponding results for size-equivalent unperturbed spheres are also shown (dashed line). For the case in which we keep ε fixed at 0.01 (blue line), ω, g, and C bak computed for spheres and 3D Chebyshev particles are similar for the smallest particles, but they already start diverging at about r 0 = 2 μm. For the case in which we keep the perturbation amplitude A fixed at 0.11λ (red line), ω computed for 3D Chebyshev particles is lower by about 0.1 than that computed for spheres over the entire size range, while g is higher for 3D Chebyshev particles by about 0.13 than the corresponding value for spheres for all sizes. These are surprisingly large differences that may even be important in radiative forcing computations. The results suggest that neglecting the effect of small scale surface roughness results in too high values of ω and too low values of g. In radiative transfer computations, this would result in too much total scattering in relation to absorption, and too much side- and backscattering, both resulting in too much aerosol cooling. Thus, these two sources of error would be additive; the homogeneous sphere model is expected to predict a larger radiative cooling effect than the 3D Chebyshev model.
Fig. 3 Optical properties of spheres (dashed line), and 3D Chebyshev particles with ε =0.01 (blue) and A=0.11λ (red): ω (top left), g (top right), and C bak (bottom left). Also shown is the CPU time (bottom right) for computations using group theory (circles) and not using group theory (pluses).
Perhaps the most remarkable result is the large difference in C bak computed for spheres and 3D Chebyshev particles, which increases with particle size. For the largest particles, C bak computed for spheres is almost 6 times larger than that computed for 3D Chebyshev particles. This can have important consequences for interpreting lidar observations of the backscattering coefficient of mineral dust particles. Our results suggest that model particles that neglect the effect of surface roughness may significantly overestimate the backscattering coefficient. If used in a retrieval method, the retrieval algorithm would interpret a lidar return signal backscattered on rough dust particles by underestimating the particle concentrations, so that the product of the low particle concentration and the high value of C bak would reproduce the observed backscattering coefficient.
Figure 4 compares the Mueller matrix elements F 11 (left column) and −F 12 /F 11 (right column) of spheres (blue) and randomly oriented 3D Chebyshev particles (red) for particle sizes of r 0 = 1 μm (top row) and 6 μm (bottom row). For both sizes, we see that the oscillation of both elements as a function of the scattering angles are qualitatively similar, but the amplitude is larger for 3D Chebyshev particles than for homogeneous spheres. Most importantly, we see that for scattering angles larger than about 30°, spheres predict considerably larger values of F 11 than 3D Chebyshev particles. In particular, this explains the differences in the backscattering cross section C bak.
Fig. 4 Mueller matrix elements F 11 (left column) and –F 12 /F 11 (right column) of 3D Chebyshev particles (red) and spheres (blue) with particle radius 1 μm (top row) and 6 μm (bottom row).
3.4. Reduction of CPU time requirements by the use of group theory
As mentioned earlier, the use of symmetries is expected to reduce CPU time requirements by a factor of about 1/(4ℓ)2. In the Tsym program, a more detailed estimate shows that the actual reduction is roughly on that order, but slightly smaller. This is due to the optimized T-matrix truncation scheme used in the Tsym program, which is identical with the truncation method used in mieschka — see [14] for details.
Figure 3 (lower right) shows the CPU time in seconds as a function of particle size. The curve marked with circles shows the actual CPU time used in the calculations that fully exploit the symmetries of the 3D Chebyshev particles. The curve marked with pluses shows the CPU time that would have been required without the use of symmetries. The values shown are theoretical values based on computing the number of extra numerical operations that would be required in Tsym without the use of group theory. To check the correctness of the theoretical predictions, we ran the T-matrix code for 3D Chebyshev particles of r 0 = 1 μm with all symmetries switched off. The CPU time of this calculation is indicated by the red square in the figure. It agrees well with the theoretical prediction. For larger particle sizes, performing the computations without symmetries would require large computational resources.
For r 0 = 1 μm the computation without the use of symmetries takes about 1.25 hours. The use of symmetries reduces the computations time to 0.5 seconds. So, group theory helps us to save about 4 orders of magnitude in computations time in this case. For r 0 = 7 μm using symmetries results in a CPU time of 7.25 minutes. Without symmetries, the same computation is estimated to take about 4.5 years! Thus the use of symmetries saves between 5–6 orders of magnitude of CPU time in this case.
The CPU times of many electromagnetic scattering methods scale with size parameter x according to a power law, i.e. CPU ∼ xL with some power L that depends on the method. A power-law fit of the curves in Fig. 3 (lower right) reveals that the CPU-time with symmetries scales like CPU ∼ x 3.5, while the calculations without the use of symmetries give CPU ∼ x 5.5. Thus exploitation of group theory does not only reduce the CPU-time for any given particle size; it actually reduced the size-scaling by 2 powers! This remarkable result is related to the way in which we treat small-scale surface roughness in this application. We keep the perturbation wavelength Λ constant at Λ = λ/4, which, as we discussed earlier, is a reasonable assumption for treating small-scale surface roughness. However, since for 3D Chebyshev particles Λ = 2πr 0 /ℓ, this means that ℓ has to increase linearly with r 0, so it increases linearly with size parameter, ℓ ∝ x. Further, since the order M of the symmetry group of Chebyshev particles scales like M ∼ ℓ, and since the CPU-time reduction is proportional to M −2, this means that the achieved reduction in CPU-time scales like x −2. Thus, if the computation time without symmetries scales like CPU ∼ xL, then with the use of symmetries we obtain CPU ∼ xL −2.
4. Summary and conclusions
We have presented an approach for modeling electromagnetic scattering by particles with small-scale surface roughness. The method is based on numerically exact electromagnetic scattering computations, which are often severely limited in the range of accessible size parameters. Usually, the problems are particularly severe for non-axisymmetric particles. The limitations are caused by (i) numerical ill-conditioning problems and (ii) CPU-time requirements that rapidly increase with size parameter. In our approach we combine two different ideas to address these problems.
The main approximation we make is to impose symmetry assumptions for the structure of the small-scale surface roughness. However, we do not limit our method to axisymmetric symmetry; so we are able to account for the effect of 3D surface roughness. The symmetry assumptions allow us to exploit group theory for making the computations sufficiently expedient. For Chebyshev particles with a fixed perturbation wavelength, computation times are reduced by 4–6 orders of magnitude, and the scaling of the CPU-time with size parameter x is reduced by two powers from CPU ∼ x 5.5 (without symmetries) to CPU ∼ x 3.5 (with symmetries).
To alleviate numerical ill-conditioning problems we use a perturbation expansion approach of the T-matrix. This approach is ideally suited for particles with small-scale surface roughness. For the geometries considered, it turned out that a perturbation expansion carried out to sixth order was sufficient to obtain numerically accurate results. The computation time required for performing the perturbation expansion of the T-matrix was only about 20–25 % of that needed for computing the matrices Q and Rg Q. Without the perturbation expansion, we obtained numerically stable results for size parameters up to x=14. With the perturbation approach, we performed computations up to x=70, which is an increase by a factor of 5!
The method was implemented into the Tsym code, and the accuracy of the results was tested by performing direct comparisons with mieschka within the range of size parameters accessible to a non-perturbative T-matrix code. In addition, the reciprocity condition for the polarized differential scattering cross section was used as a necessary condition for the accuracy of the results. These tests were performed for size parameters up to x=70, and the reciprocity condition was found to be satisfied with high accuracy. These results demonstrate that the use of the perturbation approach allows us to considerably extend the range of size parameters in our T-matrix calculations. We emphasize that this depends on the perturbation amplitude we chose in the test cases. For larger amplitudes, the range of accessible size parameters is likely to be smaller, while for smaller amplitudes the size range will be larger.
The computational results we showed were mainly meant to illustrate possible applications of the method. They underline the potentially high impact of small-scale surface roughness on the optical properties of dielectric particles. The single scattering albedo, the asymmetry parameter, and especially the backscattering cross section were strongly modulated by the presence of small-scale surface roughness. This can have important implications for the interpretation of lidar remote sensing measurements, and possibly even for the radiative forcing effect of mineral dust aerosols.
Our study was limited to Chebyshev particles with a spherical base geometry. In that case, the inversion of the unperturbed Q-matrix becomes a trivial task. However, one could equally well apply the method to other base geometries, such as Chebyshev spheroids, which have been considered in [2]. Since T-matrix computations tend to be much more well-conditioned for smooth spheroids than for Chebyshev spheroids, the perturbation method is expected to significantly improve the numerical stability of T-matrix computations for such particles.
We emphasize, once more, that the ideas on which this study was based were rather general. It would be possible to apply both the group theoretical method and the perturbation approach to model particles with surface perturbations other than Chebyshev polynomials. The group theoretical approach merely requires that the perturbations be symmetric, while the perturbation approach works best for small perturbation amplitudes and perturbation wavelengths. However, it will require careful examinations to study the size-parameter ranges and the CPU-time reductions that can be achieved when applying the method to other geometries. Further, we point out that the general ideas of this study may not be limited to traditional T-matrix computations with Waterman’s extended boundary condition method [9]. For instance, an extention of the T-matrix concept known as the shape matrix has recently been applied to Chebyshev spheres [17]. Group theoretical methods would lend itself easily for applications in shape matrix computations.
There are many open questions on the significance of small-scale surface roughness that need to be addressed in future studies. We will have to better understand the effect of small-scale surface roughness as a function of dielectric properties, particle size, and perturbation amplitude. Also, the differences between irregular surface roughness and regular 2D and 3D surface roughness needs to be studied comprehensively for particles of different sizes and refractive indices. Finally, the effect of surface roughness in relation to other morphological features, such as aggregation or perturbations with low Λ-values, need to be investigated. The main purpose of this work was to establish an accurate and expedient method that will be a useful tool in such future studies.
Acknowledgments
M. Kahnert acknowledges funding from the Swedish Research Council under contract 80438701.
References and links
1. M. Kahnert, T. Nousiainen, and P. Mauno, “On the impact of small-scale surface roughness and non-sphericity on the optical properties of hematite aerosols,” J. Quant. Spectrosc. Radiat. Transfer (to be published).
2. T. Rother, K. Schmidt, J. Wauer, V. Shcherbakov, and J.-F. Gaeyt, “Light scattering on Chebyshev particles of higher order,” Appl. Opt. 45, 6030–6037 (2006). [CrossRef] [PubMed]
3. J.-C. Auger, G. Fernandes, K. Aptowicz, Y.-L. Pan, and R. Chang, “Influence of surface roughness on the elastic-light scattering patterns of micron-sized aerosol particles,” Appl. Phys. B 99, 229–234 (2010). [CrossRef]
4. M. I. Mishchenko, V. P. Tishkovets, L. D. Travis, B. Cairns, J. M. Dlugach, L. Liu, V. K. Rosenbush, and N. N. Kiselev, “Electromagnetic scattering by a morphologically complex object: fundamental concepts and common misconceptions,” J. Quant. Spectrosc. Radiat. Transfer 112, 671–692 (2011). [CrossRef]
5. K. Muinonen, “Light scattering by stochastically shaped particles,” in Light Scattering by Nonspherical Particles, M. I. Mishchenko, J. W. Hovenier, and L. D. Travis, eds. (Academic Press, 2000), pp. 323–354. [CrossRef]
6. F. M. Kahnert, J. J. Stamnes, and K. Stamnes, “Application of the extended boundary condition method to homogeneous particles with point group symmetries,” Appl. Opt. 40, 3110–3123 (2001). [CrossRef]
7. M. Kahnert, “Irreducible representations of finite groups in the T matrix formulation of the electromagnetic scattering problem,” J. Opt. Soc. Am. A 22, 1187–1199 (2005). [CrossRef]
8. M. Kahnert, “Boundary symmetries in linear differential and integral equation problems applied to the self-consistent Green’s function formalism of acoustic and electromagnetic scattering,” Opt. Commun. 265, 383–393 (2006). [CrossRef]
9. P. C. Waterman, “Matrix formulation of electromagnetic scattering,” Proc. IEEE 53, 805–812 (1965). [CrossRef]
10. F. M. Schulz, K. Stamnes, and J. J. Stamnes, “Scattering of electromagnetic waves by spheroidal particles: a novel approach exploiting the T-matrix computed in spheroidal coordinates,” Appl. Opt. 37, 7875–7896 (1998). [CrossRef]
11. T. A. Nieminen, H. Rubinsztein-Dunlop, and N. R. Heckenberg, “Calculation of the T-matrix: general considerations and application of the point-matching method,” J. Quant. Spectrosc. Radiat. Transfer 79–80, 1019–1029 (2003). [CrossRef]
12. T. Rother and J. Wauer, “Case study about the accuracy behavior of three different T-matrix methods,” Appl. Opt. 49, 5746–5756 (2010). [CrossRef] [PubMed]
13. W. Greiner and J. Reinhardt, Quantum Electrodynamics (Springer, 2008).
14. T. Rother, Electromagnetic Wave Scattering on Nonspherical Particles (Springer, 2009). [CrossRef]
15. F. M. Kahnert, J. J. Stamnes, and K. Stamnes, “Application of the extended boundary condition method to particles with sharp edges: a comparison of two different surface integration approaches,” Appl. Opt. 40, 3101–3109 (2001). [CrossRef]
16. O. Muñoz, H. Volten, J. W. Hovenier, M. Min, Y. G. Shkuratov, J. P. Jalava, W. J. van der Zande, and L. B. F. M. Waters, “Experimental and computational study of light scattering by irregular particles with extreme refractive indices: hematite and rutile,” Astron. Astrophys. 446, 525–535 (2006). [CrossRef]
17. D. Petrov, Y. Shkuvatov, and G. Videen, “Analytical light-scattering solution for Chebyshev particles,” J. Opt. Soc. Am. A 24, 1103–1119 (2007). [CrossRef]
