Interference corrections to light scattering and absorption by metal nanoparticles

Interference corrections to light scattering and absorption by metal nanoparticles V. V. Maksimenko1,2, V. A. Zagaynov1,2, S. Yu. Krylov3, and I. E. Agranovski1,4,∗ 1National Research Nuclear University MEPhI (Moscow Engineering Physics Institute), 31, Kashirskoe shosse, 115409, Moscow, Russia 2Karpov Institute of Physical Chemistry, Vorontsovo Pole 10, 105064, Moscow, Russia 3Institute of Physical Chemistry and Electrochemistry, RAS, Leninsky prospect 31, 119071, Moscow, Russia 4School of Engineering, Griffith University, Brisbane, 4111 QLD, Australia ∗E-mail: I.Agranovski@griffith.edu.au


Introduction
New interest in the light localization problem was ignited by recent investigations of so-called random lasers where scattering powders with stimulated emissions were used [1]. A standard theoretical approach to attack the problem of localization is presented in [2][3][4][5][6]. Within their frameworks, the Bethe-Salpeter (BS) equation for an averaged two-photon propagator was reduced to a kinetic equation. The diffusion coefficient for electromagnetic energy, the transport velocity, and other parameters are indispensable attributes of this description.
In this paper we suggest some alternative approaches to attacking the problem. For Rayleigh's particles, within the same approximations used in the standard approach (summation of the most-crossed diagrams), the BS equation can be resolved directly in coordinate representation. Summation of these diagrams reduces to some corrections of classical results related to electromagnetic properties of dispersed media. In this case, there is no necessity to introduce the EM-energy diffusion coefficient, and the main localization criterion related to this vanishing coefficient becomes unnecessary.
The BS equation kernel,which is composed of the most-crossed diagrams, describes the interference of the probability amplitudes corresponding to two alternative ways of closed loop passing on trajectories of virtual photons. A main goal of this paper is to investigate the influence of this interference on elastic scattering and effective absorption of the light by a random dense system of non-absorbing small metallic particles. We propose the name "localization" for this interference PTEP 2015, 013I01 V. V. Maksimenko et al. phenomenon, intentionally avoiding the common term "strong localization," as the latter is related to the vanishing of the diffusion coefficient. Our "localized" photon is described by a pole of the BS equation, and consequently could be considered as a bound state of two virtual photons passing the closed loop on two alternative trajectories (clockwise and counterclockwise) [7,8].

Basic equation
The propagation of a photon in a system of small spherical metallic particles randomly distributed in a space is considered. The particle radius is R λ, where λ is the wavelength of incident radiation. We assume that the average distance between the particles is of the order of R. We further assume that electrodynamic properties of a separated particle are determined by an electron gas interacting via the Coulomb potential, and u(r) is the uniform ion potential. The electron density is defined as n(r) = n 0 η(r), where n 0 is the bulk electron density of the particle, η(r) = 1 inside and η(r) = 0 on the outside of the particle. The system of units where = c = 1 is used.
The general Hamiltonian of the system "metal particle + electromagnetic field" has the form where the Hamiltonian describes the noninteracting conduction electrons and free transversal electromagnetic field, is the electron density operator, ψ + (r) and ψ(r) are the electron creation and annihilation operators, and a + k,λ , a k,λ are the creation and annihilation operators for a photon with wave vector k, polarization λ (λ = 1, 2), and frequency ω = ck. The interaction Hamiltonian has the form where Q (r) = |r| −1 , δρ (r) = ρ (r) − n (r), is the electron-momentum density operator, A is the vector potential operator V is the normalization volume, and e k,λ is the unit polarization vector obeying the relations ke k,λ = 0, e k,λ e k,λ = δ λ,λ .
Localization manifests itself both in the elastic channel and in absorption. In this section, we propose quantum mechanical calculation of these corrections. Let us start from the absorption. The adsorption probability of a quant with frequency ω is determined by the Fermi "golden rule": where the amplitude corresponds with the S-matrix by the relation the summation is carried out over all excited states |s of electron gas, and ω s correspond to excitation energies. The S-matrix is determined in the standard way: where T is the chronological-ordering operator, and where superscript I indicates that the corresponding value is written in interaction representation.
In order to calculate W 1 the following identity for commutators is used: On the other hand, Using the operator e k,α ∇ k,α on the right-hand sides of (2) and (3), the following relation can be obtained: − ω s + k 2 2m e k,α ∇ k,α 0 ρ(r)e ikr dr s The last term on the right-hand side equals zero. Using equation (4), it is possible to relate the Fourier component of the irreducible polarization operator momentum-momentum P αβ (k, k |ω) to the Fourier component of the irreducible polarization operator density-density P 00 (k, k |ω), and P 00 (r, r |ω) is related to P 00 (k, k |ω) by Fourier transformation. The second term in brackets is r 0 ω/c times smaller than the first one (r 0 is the Compton electron wavelength) and can be omitted in calculations. Using the following Lehmann-type spectral representation for P 00 (k, k |ω) [9], and comparing (1) and (5) yields: The summation of the PT series is reduced to a change of the irreducible polarization operator P 00 by the corresponding reducible operator 00 : The density operators here are taken in the Heisenberg representation: ρ(r, τ ) = e i Hτ ρ(r)e −i Hτ .
The reducible operator 00 is calculated by the following procedure [10]: 00 r, r |ω = P 00 r, r |ω + e m 2 P 0α (r, r 1 |ω ) D αβ (r 1 , r 2 |ω ) P β0 r 2 , r |ω dr 1 dr 2 , where D αβ is the photon propagator obeying the Dyson equation D a αβ r, r |ω = D 0 αβ r, r |ω + D 0 αγ (r, r 1 |ω ) π a γ ν (r 1 , r 2 |ω ) D a νβ r 2 , r |ω dr 1 dr 2 , (6) D 0 is the free photon propagator in gauge with zero scalar potential, and π a αβ r, r |ω = − is the interaction "potential" of a photon with a particle characterized by radius vector a, ε(ω) = 1 − ω 2 0 /ω 2 is the dielectric permeability of the electron gas, θ(x) is the unit step function, and ω 0 is the classical plasma frequency of unbounded electron gas [11]. The irreducible mixed polarization operators momentum-density and density-momentum are determined by the following relations [10]: In the high frequency approximation (ω ε F /N 1/3 , where ε F is the Fermi energy and N is the total number of conducting electrons in a particle) [10,12]: As a result, In order to calculate the scattering cross section, the two-photon propagator is averaged over the positions of all particles. The latter are determined by the Bethe-Salpeter equation. In graphical representations of this equation, the unconnected diagrams describe coherent scattering in the forward direction, and are beyond our consideration. The incoherent scattering is described by the so-called irreducible four-point vertex which is a sum of all connected diagrams. The differential scattering cross section is obtained by joining wave functions of incident and scattered photons with the block , as shown in Fig. 1. Here, the vertical and inclined dotted lines represent interaction potentials and the expressions correspond to the wavy lines. Normalization with one photon in the volume is used: f k (r) = e k exp(ikr). The factor g = f /(1 − f ) is associated with any circle (the appearance of the packing factor f in the denominator is related to the particle volume).
The total scattering cross section in a system of particles is obtained by averaging the following expression over the particle coordinates [13][14][15]: where the summation is carried out over all particles of the system, and n f is the unit polarization vector along the direction of the scattered quantum. Indexes a and b mean that the initial and final particles on the photon trajectory could be different. If a = b, we deal with the classical cross section which is the square of the modulus of the scattering amplitude aa on this particle. If a = b, some spacing between the entrance and exit points appears where localization becomes possible. As seen from the structure of (11), this cross section could be a complex value.  The scattering amplitude ab on an arbitrary pair of particles a and b is determined by the following perturbation series [13][14][15]: The origin of the "off-diagonal" elements of σ ab (a = b) in the cross section is related to the contribution of the so-called most-crossed diagrams contained in block in Fig. 1. These two-storey diagrams resemble those describing multiple interactions of two real photons. In reality, we deal with the effective interaction of the pair of virtual photons passing the closed loop clockwise and counterclockwise. This interaction is related to their scattering on the same particles. From this point of view, the localization may be considered as the generation of a bound state of this pair.
We choose the block composed of the most-crossed diagrams presented by the second line in Fig. 1 as the irreducible vertex. As is known, the Ward identity determines some relation between the irreducible four-point vertex and the mass operator (third, fourth, and fifth lines) in the Dyson Maksimenko et al. equation for the one-photon propagator D in the system (the sixth line). In order to satisfy this identity, the photon line of an arbitrary diagram for is cut in the middle and its right and left parts are separated on the lower and upper banks of the corresponding diagram for . The point is that the amplitudes of the processes corresponding to the banks with different numbers of particles do not interfere. Regarding the t-matrix of a photon scattering on a separate particle, the latter is calculated within the coherent potential approximation (the last line) [16]. The sequence order of the tensor indexes and arguments is explained by the second line of Fig. 1.
The first diagram of the series describes typical incoherent elastic scattering. Subsequent diagrams of the magnetic circular dichroism (MCD) series describe the interference of a pair of virtual photons passing the closed loop by two alternative ways, clockwise and counterclockwise.
We suppose that the one-photon propagator in a heterogeneous medium looks like the propagator in some homogeneous medium with an unknown effective dielectric permeability. Therefore, the first of our problems is reduced to determination of this permeability.
The interaction amplitude is the solution of the following equation: The photon propagator in the system of particles D satisfies the Dyson equation where is the photon mass operator determined by the expression The equation for the t-matrix of the photon scattering on a separate particle has the form t a αβ r, r = π a αβ r, r + π a αγ (r, We assume that the photon propagator D in the system is similar to the propagator in a uniform effective medium characterized by the longitudinal ε l and transverse ε t dielectric permeabilities: Substituting expression (15) for D into Eqs. (12)- (14), one can obtain equations for ε l and ε t based on a self-consistency condition. Notice that 00 (r, r |ω) could be presented as follows: where N is determined by the expression where brackets designate integration over the internal coordinate and summation over the tensor indexes. The series for the mass operator could be created inside the brackets. Finally,

Effective dielectric permeability
Firstly, Eq. (14) can be rewritten in the form It is assumed that a particle is surrounded by some effective medium with unknown dielectric permeabilityε αβ (r, r ). Then the interaction potential has the following form (compare with expression (6)): where ε αβ (r, r ) is the effective dielectric tensor, related to the longitudinal ε l and transverse ε t dielectric permeabilities of the medium by the expression [17]ε Using the identity the explicit form of ε αβ (r, r ) can be obtained: Now, the t-matrix inside a particle could be calculated. For R λ the term containing the small parameter ω R/c characterizing retardation of electromagnetic interaction in the propagator D in Eq. (17) could be neglected. Thus D is presented in the following form: Now the solution of Eq. (17) is presented in the form Substituting expression (22) into Eq. (17) yields the following equation for the function M: where a = (ε − ε t )/4π , b = (ε t − ε l )/4π , and the operator L is defined as The solution of Eq. (23) is obtained in the form of an expansion in scalar spherical harmonics. In the dipole approximation: where Using expression (23), Eq. (12) is rewritten as follows: αβσ ν r, r ,r ,r = δ αβ δ σ ν δ(r − r )δ(r −r )ϕ(r −r ) where and t = R − |r −r| /2. The substitution αβσ ν r, r ,r ,r = δ αβ δ σ ν δ(r − r )δ(r −r )ϕ(r −r ) + ∂ 4 ∂r α ∂r β ∂r σ ∂r ν r, r ,r ,r gives the following equation for the function : ∂ 4 ∂r α ∂r β ∂r σ ∂r ν r, r ,r ,r = ϕ(r −r )ϕ(r −r)D αβ (r, r )D * σ ν (r ,r) Block describes the photon motion along the closed loop. Therefore, its upper and lower coordinates are separated by a distance of the order of R. In other words, this block should be local with 9/19 PTEP 2015, 013I01 V. V. Maksimenko et al. respect to the upper and lower coordinates, which essentially simplifies the subsequent calculations. Substituting (22) into Eq. (29), using integration by parts using (20), and taking into account only singular terms related to double gradients of the functions Q, the following expression for can be obtained: αβσ ν r, r ,r ,r ≡ αβσ ν r −r , r −r, r − r = αβσ ν (r 1 , r 2 , r 3 ) where γ = 4π c 2 /(|ε l | ω 2 ). Now, using (13), let us calculate the mass operator αβ (r, r ) and its Fourier transform. Then, using the Dyson equation in impulse representation, one could calculate the Fourier transform of the propagator D(k, ω). If the form of Eq. (15) of the propagator was correctly chosen, this procedure could provide explicit expressions for ε l and ε t . As a result,ε where x is the solution of the following equation: At f → 0,ε =ε l =ε t satisfies the classical equation of effective medium approximation [18]:

Elastic scattering and absorption
It is useful to present in the following form: αβσ ν (r, r ,r ,r) ≡ αβσ ν (r 1 , r 2 , r 3 ) = δ αβ δ σ ν δ(r 3 )δ(r − r 2 )ϕ(r 1 ) where r 1 = r −r , r 2 = r −r, and r 3 = r − r . In accordance with Fig. 1, the averaged differential cross section related to MCD is expressed via in the form where V 0 is the system volume and k i and k f are the wave vectors of incident and scattered photons. As shown below, the integral in (34) is a nonanalytic function of the scattering direction. Carrying out routine integration over r 3 yields (31) Obvious symmetrization transforms the integral in the expression (35) to the following form: where q = k + m and m = 1 2 k i + k f . It is obvious that Hence, the expression (35) takes the form The last expression is not defined for k i = −k f . This corresponds to strictly backward scattering. For k i = −k f we ought to attack the following integral: where ψ(r) = ϕ(r)/(1 − γ 2 ϕ(r)). If ω R/c 1, the last expression can be transformed to the following form: δ αβ δ σ ν + δ ασ δ βν + δ αν δ βσ dx.
The total cross section integrated over the angles can be conveniently represented as a sum of three components, The first of these corresponds to the first diagram of the MCD series and describes the usual incoherent scattering by the system of particles. The terms in the brackets are related to interference corrections of the cross section. As expected, these components are complex values. What does it mean?
In order to answer this question let us use the continuity equation where ρ is the density of photons inside an arbitrary domain containing the particles and j is the flux density of radiation. Integrating over the volume of this domain, we have where J = jdS is the full radiation flux through the surface and N is the number of photons inside the domain. Taking into account that the flux scattered along the unit vector n f is related to the differential cross section dσ/dn f by the relation where I 0 is the flux density of incident radiation, the relation (42) can be rewritten in the form where σ = dσ/dn f dn f is the total scattering cross section. If the cross section was purely imaginary (σ = iσ 2 ), the only possibility to satisfy equation (43) would be the introduction of the two following assumptions: 1) introduction of a certain number N 0 of photons oscillating inside this volume, 2) the flux density of incident radiation I 0 stipulated namely by these photons also oscillates, The mean number N 0 of these photons is determined by separation of the imaginary part of Eq. (43): Now we can consider the real part σ 1 of the cross section. In the case of the appearance of σ 1 , the number of virtual photons in the volume still oscillates, but their number decreases with time.  The same decrease is also related to I 0 : Substituting N 0 and I 0 in Eq. (43) and separating the real and imaginary parts of this equation, we can write another equation determining an exponent δ, Thus, the "mean" number of oscillating photons could be presented as: The characteristic lifetime of such photons is τ ∼ (σ 2 /σ 1 )T , where T is the period of the incident electromagnetic wave. Note, the photons under consideration are virtual and both N 0 and I 0 describe just those which correspond to neither the source nor the detector of radiation and are not observable. These photons differ from real ones. The latter are either escaping from the source or could be detected by an observer. At a given time, the virtual photon could be located at any position, with no information about coordinates available. Being spread across the space, its wave function is rapidly transformed to a point solely when the photon is registered (so-called reduction of the wave function). The time t and coordinate r for the virtual photon are two independent parameters, which is also valid for their Fourier transforms ω and k. Therefore, the virtual photon is rather described by the propagator D 0 (k, ω), and not by the plane electromagnetic wave.
The angular dependence of the imaginary part of dσ/dn f for both possible polarizations is shown in Fig. 2. Then one could neglect with s-scattering because of its negligibly small magnitude in comparison with p-scattering. Figure 3 presents the ratio σ 2 /σ 1 determining the characteristic lifetime of the localized photon.   It is crucial to calculate the group light velocity v within the frequency range where σ 2 appears. Using the dispersion low ω = kc/ ε t , the following expression for v could be written: The frequency dependence of Imε t , whose derivative determines the effective light velocity, is shown in Fig. 4. In the region where |d(Imε t )/dω| 1, the group velocity of electromagnetic waves tends to zero (here, as it can be shown that Reε t Imε t ). Although the values σ 1 and σ 2 describe properties of the localized (unobservable) photons, their behavior requires some attention. The angular dependence of dσ 1 /dn f shown in Fig. 5 is of special 15  interest. On one hand, this value is an addition to the "real" scattering. On the other hand, it determines the behavior of already localized oscillating photons. As is known, the amplitudes of both alternative ways of loop passing are always subjected to constructive interference. The constructive interference ought to increase the probability of looping, and scattering to the back hemisphere should increase respectively. This is seen in Fig. 5, where the scattering indicatrix for both polarizations per one particle under the condition of dense packing is presented. The Rayleigh indicatrix for a separate particle is also presented for comparison. The angular width of the slit and cog in the indicatrices of p-and s-polarized light, respectively, are equal to zero (in the first case the polarization vector belongs to the scattering plane, in the second case this vector is orthogonal to it).
Some calculation in accordance with formula (16) give the following expression for absorption scattering per one particle: Terms of the order of (ω R/c) 3 are omitted here. The frequency dependence of the absorption cross section related to the localization is presented in Fig. 6. According to traditional views, this "absorption" is very large. It takes place in the same σ 2 frequency interval. The absorption spectrum has the characteristic two-headed structure. The reason for that is as follows. At the frequency corresponding to the condition dω/dk = 0 (see Fig. 4), the light group velocity becomes zero, and the idea of motion of the virtual photon along closed loops loses any sense. The two-peak structure of the absorption spectrum of silver hydrosol was observed in [18]. The localization depends critically on the packing factor f . As is seen from Fig. 7, there are two localization thresholds-upper and lower exist. Below and above these thresholds, the scattering becomes insufficient in order to produce closed loops on the photon trajectories.
The localization phenomenon sheds a new light on the reason for the well-known problems of numerous effective medium approximations [19]. Within these theories, the effective dielectric permeability of a collective of small non-absorptive metallic particles, in a certain interval of frequencies and packing factors, becomes a complex value. This means that some mysterious absorption appears there. If the difference between the longitudinal and transverse effective dielectric permeability is not taken into account, the system is characterized byε as determined by Eq. (33), which can be considered as one of the variants of effective medium approximations. Here, a region of ω and f where 16  ε is a complex value also exists. However, the suggested approach clarifies the issue-the effective absorption is related to the localization.

Discussion
The common procedure for the description of the light localization is based on reduction of the Bethe-Salpeter equation in momentum representation to a transport equation. On this basis, an introduction of the diffusion constant of electromagnetic energy seems to be an unavoidable step. We proposed another approach to attack the problem. Within the framework of standard approximation (summation of maximally crossed diagrams under calculation of four-point irreducible vertex ) we resolved the latter equation directly in coordinate representation. This approach was made possible due to the locality of with respect to some of its arguments as describes the photon movement along the closed loop.
Calculation of integrals in the cross section and bypassing the pole of lying on the real axis of the complex plane causes the appearance of some imaginary parts of the cross section. One used to think that the scattering cross section is a purely positive value determined by the square of the 17/19 PTEP 2015, 013I01 V. V. Maksimenko et al. modulus of the corresponding amplitude. However, for a system of many particles this is not valid. In this case, the differential cross section is obtained from averaging over the positions of the particles of the value determined by expression (11), where ac is the amplitude of the photon transition from particle "a" to particle "c," and summation is carried out over all particles. The structure of this expression explains the complexity of the cross section.
The reason for the localization is related to specific interference effects existing even in an absolutely disordered system of particles. The photon is entangled in the system of loops on its trajectory. It is important to consider that the loop could be passed by two alternative ways, clockwise and counterclockwise. The phase shifts corresponding to the two possible ways are zero. Consequently the corresponding amplitudes interfere constructively. Therefore, the probability of the loop generation increases dramatically. Any loop means a return backward. This increased probability stimulates the generation of new loops, etc. This self-sustained process results in the photon localization.
As was shown, only p-polarized light (when the polarization vector belongs to the scattering plane) could be localized. The localized photon possesses the angular momentum related to its motion along the closed loop. If the photon scattering is accompanied by localization, the light changes the polarization vector direction. The p-polarized light always changes its polarization direction upon scattering to any non-zero angle. The s-polarized light never changes the direction of the polarization vector after scattering on spheres or a collective of spheres. Therefore, the p-polarized photon could be localized, whilst the s-polarized one could not. Consequently, the peculiarities of the σ 1 indicatrix for scattering to the back hemisphere (see Fig. 5) are related to the p-polarization only, while the s-indicatrix remains unchanged, with only one exception. The polarization plane becomes indefinite for light scattering in the forward and backward directions. In these cases there is no difference between s-and p-polarized light. Respectively, there is no difference between scattering in the forward and backward directions. This is illustrated in Fig. 5. The appearance of the slit and cog in the indicatrix of p-and s-polarized light, respectively, are the results of this effect.
As is known, one more type of interference correction to the elastic scattering in a system of many particles exists. The point is that the amplitudes of direct and opposite photon passage of an arbitrary chain of particles always interfere constructively, because the phase shift in both cases is the same (k → −k; r → −r). The ladder diagrams describe the interference of just these amplitudes. Here also the analogical pole arises. It can be shown that the contribution of ladder diagrams to the scattering would be considerably smaller than the MCD case.