Fourier modal method for the description of nanoparticle lattices in the dipole approximation
Ilia M. Fradkin, Sergey A. Dyakov, Nikolay A. Gippius
TL;DR
This paper introduces a fast and accurate method combining discrete dipole approximation with RCWA to analyze optical properties of nanoparticle lattices, enabling detailed study of plasmonic resonances in layered structures.
Contribution
It develops a novel approach integrating DDA with RCWA for efficient analysis of nanoparticle arrays, improving speed and accuracy over existing methods.
Findings
Successfully models localized surface plasmon resonances.
Accurately predicts lattice plasmon resonances and their hybridization.
Demonstrates high accuracy and fast convergence compared to other computational methods.
Abstract
Rigorous coupled-wave analysis (RCWA) is a very effective tool for the studying optical properties of multilayered vertically invariant periodic structures. However, it fails to deal with arrays of small particles because of high gradients in a local field. In this thesis, we implement discrete dipole approximation (DDA) for the construction of scattering matrices of arrays of resonant nanoparticles. This strongly speeds up the calculations and therefore provides an opportunity for thorough consideration of various layered structures with small periodic inclusions in terms of the RCWA. We study in detail three main stages of the method: calculation of polarizability tensor of a single nanoparticle, effective polarizability of this particle in a lattice and corresponding scattering matrix of the layer for further integration in the conventional RCWA approach. We demonstrate the…
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Fourier modal method for the description of nanoparticle lattices in the dipole approximation
Ilia M. Fradkin [email protected] Skolkovo Institute of Science and Technology, Nobel Street 3, Moscow 143025, Russia
Moscow Institute of Physics and Technology, Institutskiy pereulok 9, Moscow Region 141701, Russia
Sergey A. Dyakov
Skolkovo Institute of Science and Technology, Nobel Street 3, Moscow 143025, Russia
Nikolay A. Gippius
Skolkovo Institute of Science and Technology, Nobel Street 3, Moscow 143025, Russia
(March 15, 2024)
Abstract
Rigorous coupled-wave analysis (RCWA) is a very effective tool for the studying optical properties of multilayered vertically invariant periodic structures. However, it fails to deal with arrays of small particles because of high gradients in a local field. In this thesis, we implement discrete dipole approximation (DDA) for the construction of scattering matrices of arrays of resonant nanoparticles. This strongly speeds up the calculations and therefore provides an opportunity for thorough consideration of various layered structures with small periodic inclusions in terms of the RCWA. We study in detail three main stages of the method: calculation of polarizability tensor of a single nanoparticle, effective polarizability of this particle in a lattice and corresponding scattering matrix of the layer for further integration in the conventional RCWA approach. We demonstrate the performance of the proposed method by considering plasmonic lattices embedded in a homogeneous ambiance and placed inside and onto optical waveguides and compare our results with experimental papers. Such phenomena as localized surface plasmon resonances (LSPRs) and lattice plasmon resonances (LPRs) are observed as well as their hybridization with photonic guided modes. High accuracy and fast convergence of our approach are shown by a comparison with other computational approaches. Typical limits of applicability of our approximate method are determined by an exploration of the dependence of its error on the parameters of the structure.
This paper is an extended version of our article [1].
Contents
-
A.3 Resonant approximation for polarizability of indivisual particle
-
B.4 Analytical calculation of convolution for a homogeneous environment
Chapter 1 Introduction
Periodic photonic structures are one of the most important low-level components in modern photonics since they form the basis for plenty of optical elements and devices. They include bandgap materials [2, 3, 4], diffraction gratings [5, 6, 7, 8], frequency selective surfaces [9], antennas [10, 11], waveguides [12, 13, 14], metasurfaces [15, 16, 17, 18, 19], biosensors[20, 21] etc. Spatial periodicity naturally suggests that the electromagnetic field in periodic structures can be expanded into the Fourier series, that appears to be an effective tool for calculating of optical properties of such structures. Indeed, this is a basic idea of the rigorous coupled-wave analysis (RCWA) [22, 14] proved itself to be an extremely fast and efficient computational method. However, inclusions much smaller than a wavelength give rise to high gradients in the near field, which forces to take lots of harmonics into account and can significantly slow down this approach. This fundamental drawback of all Fourier-modal methods, which originates from the Gibbs phenomenon [23, 24], is most pronounced for high-contrast inclusions, such as plasmonic nanoparticles. Moreover, occurring of localized surface plasmon resonances makes the contribution of these inclusions determinative for the optical properties of the whole structure.
To overcome this problem, several approaches have been developed. One of them is the Li factorization rules [25], which solve the problem of poor convergence at concurrent jump discontinuities. Yet another approach is the use of an adaptive spatial resolution by choosing specially designed coordinate transformation which increases the resolution around the material boundaries [26, 27]. These methods significantly improve the convergence of the RCWA numerical scheme. However, if the size of metal inclusions is much smaller than the structure spatial period, practically unrealizable number of Fourier harmonics is required for the solution of Maxwell’s equations to be converged. As a result, the RCWA fails to describe optical properties of a periodic array of metallic particles in a dielectric.
At the same time, the light scattering by relatively small particles in most cases might be approximately described by substituting them with ideal electric dipoles, which brings us to discrete dipole approximation [28] (DDA).
In this thesis, we report the approach for calculation of a scattering matrix of a plasmonic lattice, which is based on a determination of an effective polarizability tensor of nanoparticles arranged in periodic arrays. We combine three widespread numerical methods: finite element method (FEM), DDA, and RCWA in order to implement each of them on a specific stage of computation and obtain precise results in a fast way. It helps to study any layered structures with inclusions of plasmonic lattices for any angles of incident light. Such an approach makes it possible to observe dispersion of lattice waves and phenomena of out-of-plane polarization of plasmonic particles [29, 30]. To illustrate the feasibility of the proposed method, we consider the same plasmonic lattice in a homogeneous ambience, on a photonic waveguide and inside it, the formation of photonic band structure, strong coupling of photonic guided modes with both LSPRs and LPRs. Also we reproduce several experiments and compare them with our results. High accuracy and fast convergence of our approach are shown by a comparison with conventional RCWA, RCWA enhanced by adaptive spatial resolution and FEM. Typical limits of applicability of our approximate method are determined by an exploration of the dependence of its error on the parameters of the structure.
Chapter 2 Effective polarizability
Dipole approximation makes it possible to split the problem of determination of particles’ polarization into two parts, which deal with a problem in different scales. Indeed, when a particle is incorporated in a certain structure of permittivity so that this permittivity is constant over the volume of a particle then it can be conveniently substituted with a solitary dipole (see Appendix A). The dipole moment of the -th particle, , is determined as a polarizability tensor, , acting on a background low-gradient electric field, , which excites electron oscillations in a plasmonic nanoparticle:
[TABLE]
In this paper, we define dipole moment as a dipole moment of free charges, which generates the same far-field as a real particle. To give an example, such a definition makes tensor times larger in comparison with a tensor conventionally defined for a particle in a homogeneous medium [31]. Computation of this polarizability tensor is the first subproblem, which can be easily solved numerically by any near-field computational method for a particle of any shape (see Appendix A) or even analytically for particles of trivial shape [32, 33]. The second subproblem is a determination of the background field itself, , which is a sum of an incident waves’ field, , which would have been in the absence of the lattice and electric field scattered by all the neighboring particles of a lattice at the coordinate of the considered particle, :
[TABLE]
where is the dyadic Green’s function showing electric field induced at the point by a dipole at the coordinate in a considered, not necessarily homogeneous environment. It should be emphasized that in this expression Green’s function acts on a dipole moment of free charges and is defined accordingly.
Background electric field can be found immediately by solving this linear algebraic system consisting of equations ( is a number of particles) via any specialized method [34], which is a general approach for DDA method. However, occurrence of periodicity makes it possible to express the solution of algebraic system in a simple form.
Indeed, let us consider an infinite lattice consisting of the same particles. In this case according to Bloch theorem the Floquet periodicity occur, which means that . Substituting this expression into Eqn. 2.2 gives us:
[TABLE]
Application of Eqn. 2.1, finally results in a following expression, which variations has been observed in literature [34, 35, 36, 37, 38, 39, 40, 41]:
[TABLE]
where, , is the identity matrix, is the in-plain component of photon wavevector and tensor is a so-called dynamic interaction constant [40], which is a sum of dyadic Green’s function over the lattice (see Appendix B for details of practical calculation):
[TABLE]
Relation 2.4 allows to introduce effective polarizability tensor, , which connects dipole moment with an incident electric field :
[TABLE]
As it can be seen from this expression, effective polarizability, , has resonances of two types. The first one occurs, when the ordinary individual polarizability experiences wide LSPR. Another resonance with non-trivial dispersion (LPR) is associated with collective oscillations of the lattice and occurs when the condition is fulfilled.
This approach allows to describe the structure in a very simple way but naturally has inherent limitations. Indeed, dipole approximation works until background field does not change on the dimensions of a particle, which results in a requirement for a particle to be much less than a wavelength (far-field limit) and a period of a structure (near-field limit). Also, an electric field of an ideal dipole and a real particle match each other for distances larger than several particle sizes. This means that all the in the vicinity of particles should be accounted in tensor by an appropriate choice of Green’s function and polarizability of a single particle .
Dipole approximation allows considering a wide range of structures. However, the proposed approach can be naturally extended by taking into account higher order multipole moments. Even accounting for quadrupole moment makes it possible to consider larger particles, place them closer to each other and observe quadrupole resonances. For instance, this can be potentially profitable for the description of dense plasmonic metasurfaces or periodic structures of large dielectric particles.
Chapter 3 Scattering matrix calculation
Calculation of the effective polarizability, , of a particle in a lattice is an important stage. However, our final goal is to obtain the total scattering matrix of the entire structure containing resonant particles. For this hereafter, we will use the formalism of RCWA which was conveniently developed for calculating of optical properties of vertically invariant periodic structures (e.g., multilayered gratings). The total scattering matrix of a multilayered structure is calculated iteratively [42]; at each iteration step the scattering matrix of the -th layer is connected with the total scattering matrix of layers 1 to found at previous iteration step. Usually, the scattering matrix of each layer is found by decomposing the electromagnetic field into spatial Fourier harmonics.
As mentioned above, such structures as small metallic particles cause high-gradient fields, which requires taking into account a very large number of spatial harmonics to resolve them in RCWA. However, strongly evanescent behavior of high- harmonics leads to their confinement inside a layer of several dozens nanometers thickness, which contains plasmonic nanoparticles (blue region in Fig. 3.1). This means, that construction of the scattering matrix of this layer as a whole, , via DDA gives us a possibility to work only with low- harmonics (see Fig. 3.1 (a)). In such an approach, all high- effects are described by the effective , while low- effects are treated conventionally by RCWA. Once the local scattering matrix is calculated, the initial RCWA problem of calculation of the total scattering matrix takes on the task of coordination of adjacent layers, which is consistent with the original spirit of RCWA.
By definition, the scattering matrix of the considered layer, , connects amplitudes of the incoming and outgoing waves on the boundaries of the considered layer (see Fig. 3.1 (a)) [14]:
[TABLE]
Hereinafter, we use the symbols and for the amplitudes of positively and negatively propagating waves taken at positions specified by the subscripts as shown in Fig. 3.1. Please note that axis is directed from top to bottom.
Elaborating the idea discussed in the section 2 we represent this matrix as a sum of two terms . The first term corresponds to a matrix calculated in the assumption of the absence of a lattice inclusion (see Fig. 3.1 (b)), whereas the second one (see Fig. 3.1 (c)) accounts for the radiation of the currents induced in the particles. Finally, we substitute the complex layer, which contains nanoparticles lattice by a black box, which is fully described by the matrix, , which has a small number of non-zero elements, describing low- harmonics.
We start with the consideration of an empty layer without plasmonic inclusions. Scattering matrix connects amplitudes on the boundaries:
[TABLE]
However, for our purposes, it is very important to know the vector of amplitudes at the plane, which will further contain the dipole lattice. Therefore, we introduce scattering matrices and of upper and lower parts of the layer respectively. These matrices act as follows:
[TABLE]
Moreover, they are obviously connected as:
[TABLE]
where the operand denotes the combination of two adjacent scattering matrices [42] (see Appendix C for details of its calculations).
It is convenient to introduce an auxiliary matrix , which allows to determine vectors directly from the incoming waves amplitudes. Equations 3.3 allows us to represent this matrix via and components:
[TABLE]
where
[TABLE]
Vector help us to determine subsequently the field, , induced at the position of the particle and corresponding current . According to Maxwell’s equations, the presence of currents leads to an appearance of discontinuities of the horizontal components of electric and magnetic fields [43, 44]. As a result, the vectors of amplitudes taken at the coordinates infinitesimally above and below the dipole’s plane (see Fig. 3.1 (c)) are connected as follows:
[TABLE]
where is the amplitude discontinuity vector determined by the induced current , which results in a connection with incoming waves’ amplitudes via a special tensor, , (see Appendix C).
The last thing, which is left to do is the calculation of outgoing from the layer waves’ amplitudes, determined by this emission. Similarly to , it is very convenient to introduce another matrix connecting these discontinuities with outgoing waves’ amplitudes associated with them. This matrix can be expressed through scattering matrices and components from the condition of the absence of incoming waves (see Fig. 3.1 (c)):
[TABLE]
It should be especially emphasized, that although there are no incoming waves in the subproblem of dipoles’ radiation, their amplitudes depend on the incoming waves, which allows finding the correction to the scattering matrix, defined as
[TABLE]
Summing up, in order to do this, we just have to (I) calculate the incoming waves’ amplitudes at a plane of the particles, (II) apply operator in order to find the amplitude discontinuity vector and finally (III) determine the correction to the outgoing from the layer waves.
Combining equations (3.5 – 3.10), we obtain the following expression for the correction to the scattering matrix of the layer:
[TABLE]
When the local scattering matrix of a plasmonic layer is known, it can be easily inserted in any structure. Its total scattering matrix is then found as
[TABLE]
These rather general expressions can be potentially used for a description of rather complex structures, for example, plasmonic lattice inside a photonic crystal or a lattice on a metamaterial. However, in this paper, we implement our approach only for two most practical cases and it makes sense to specify matrices for them. When the lattice is embedded in a homogeneous medium, matrices take a very simple form:
[TABLE]
where and are the thicknesses of higher and lower layers included into the local environment (see Fig. 3.2) and satisfies the relation
[TABLE]
In Eqns. (3.13) and (3.14) , and are the dimensionless diagonal operators [14] and
[TABLE]
where and are hypervectors of - and - projections of reciprocal lattice vectors. Matrix also takes a simple form discussed in details in Appendix C.
The second most practical case is a lattice placed in close proximity of an interface between two homogeneous media. It should be noted, that even if the real particles lay right on the boundary, effective dipole lattice is placed at a level of the centers of these particles. Thus, for a dipole lattice situated at a distance of above an interface matrices have the following form:
[TABLE]
where is the scattering matrix of the interface and operators and are calculated above and below the interface correspondingly, and are defined according to the Fig. 3.2 b.
The matrix has the same form as for homogeneous environment (whereas used for its calculation of course differs) since the plane of particles’ centers is fully inside a homogeneous medium.
It should be emphasized once again that the dyadic Green’s function , which is used for a calculation of effective polarizability tensor , accounts for the structure of the layer (for instance, presence of an interface) in the near field of metallic particles (see Appendix B). Indeed, this fact allows us not to account for the self-influence of dipoles by means of RCWA and escape from operating with all the high- harmonics. Also, it should be noted, that, although the matrices of the same form are implemented for both considered cases they differ in tensors applied for their calculation, which, in turn, strongly depend on the presence of an interface and a distance to it.
There is an interesting fact that such splitting parameters as and can be chosen arbitrarily and even set to zero. However, the method performance is determined by the distance to the first boundary not accounted by the local scattering matrix. Indeed, the closer the boundary, the larger is the number of required harmonics and the less is our advantage over the conventional RCWA. Moreover, for distances of the order of particle’s size our approach becomes inapplicable and an additional boundary should be included in a local environment to be accounted by the polarizability tensor.
Chapter 4 Example: hybrid resonances
Plasmonic metals play an important role in modern photonics, because of unique ability to enhance light-matter interaction via strong confinement of light in subwavelength dimensions. Very well-known localized surface plasmon resonances (LSPR) and surface plasmon resonances (SPR) have low Q-factors but can compete with conventional resonances because of extremely small mode volume. In some cases, an intermediate regime of middle Q-factor and mode volume is needed, which suggests the use of hybrid plasmon-photonic resonators [45, 46]. An elegant way of implementation of this idea is a construction a regular lattice of resonant plasmonic nanoparticles. Such a structure allows employing both advantages and features of plasmonic and periodic structures. Lattice plasmon resonances that occur in them are actively implemented in biosensors [47, 48, 20], sources of light [49, 39, 50], stretchable devices [51]. They might be used for lasing [52, 53], strong coupling with emitting systems [54, 55] and other purposes [56, 57, 58].
Plasmonic lattices have already been thoroughly studied theoretically. Such universal near-field methods as FEM or finite-difference time-domain method [59] (FDTD) can be implemented for a description of any periodic structure. However, they are too computationally expensive for observation of angle-dependent spectra with a good resolution, which requires at least thousands of computations. RCWA [14, 22] specialized for periodic structures as it was mentioned before also fails in this case. Plasmonic lattices not only consist from small particles but generate high contrast in permittivities of adjacent materials as well. These factors together with exciting physical phenomena make them the most promising candidate for the first application of our approach.
For this reason, in order to show the performance of our method and illustrate the physical phenomena, which can be investigated, we consider plasmonic lattice in different environments (see Fig. 4.1 (a–c)). We have chosen a square lattice of silver nanodisks with a period of nm. Disks have the radius of nm, the height of nm and are described by Johnson-Christy optical constants [60]. In the first case, this lattice is embedded in an infinite surrounding of silica (). In the second variant of the structure, the lattice is incorporated directly in the middle of this waveguide. And in the latter one, we deposit nanoparticles on an nm membrane silica waveguide.
4.1 Homogeneous environment
To start with we have calculated extinction spectra (, where is a transmission in the main diffraction order) of the plasmonic lattice in a bulk silica for both polarizations of light incident on a lattice along its crystallographic axis (see. Fig 4.1 d,g). Several specific phenomena are observed in this structure.
So-called Rayleigh anomalies [61] (magenta dashed lines in Fig. 4.1 (d,g)) correspond to openings of different diffraction channels. In other words, for this relation of and all the particles radiate light coherently along the plane of a lattice, which results in in-phase contribution into the tensor, its divergence and subsequent vanishing of the effective polarizability . Hence, in dipole approximation lattice becomes effectively transparent under Rayleigh condition, which is clearly seen in Fig. 4.1 (d,g,j). At the same time, when we go to the energies lower than Rayleigh anomalies, interaction constant slightly decreases, which leads to a fulfillment of the resonant condition for LPRs, observed in Fig. 4.1 (d,g,j) as well. These resonances are so-called Fano-Wood anomalies [61], which occur when one of the diffraction orders matches both frequency and an in-plane component of wavevector of the guided mode. LPRs in a homogeneous or almost homogeneous ambience have been thoroughly studied both theoretically and experimentally [62, 56, 38, 63, 64, 48]. Their dispersion spectra have been observed experimentally and explained in [62].
Extinction spectrum of the plasmonic lattice in a homogeneous environment under the normal incidence is shown in Fig. 4.1(j), whereas the reflection, transmission, absorption and diffraction spectra are depicted in Fig. 4.2.
4.2 Lattice in a waveguide
When the lattice is embedded in the middle of the waveguide of finite thickness (see Fig. 4.1 (b)) photonic guided modes come into play, while the positions of LSPR and LPR almost do not change. An interaction between plasmonic resonances with guided modes leads to an appearance of hybridized resonances (see Fig. 4.1 (e,h,k)), which represents the main difference with the lattice in the bulk space. Since the lattice is located strictly in the middle of the symmetric membrane waveguide the hybridization is determined by the parity of the guided modes. Therefore, half of the modes, which have zero electric fields in the center of the waveguide are optically active.
The charachteristic example is observed for -polarized wave at and (see Fig. 4.1 (h)). In this point almost horizontal dispersion curve of lattice plasmon resonance is crossed by a second TM guided mode, which results in Rabi splitting.
Reflection, transmission, absorption and diffraction spectra for the normal incidence of light are depicted in Fig. 4.3.
4.3 Lattice on a waveguide
Despite the apparent similarity of the structures, the appearance of an interface in proximity of a lattice (see Fig. 4.1 (c)) strongly changes optical properties. The difference comes from the fact that dipole located near an interface between two media almost does not radiate light along the boundary. It can be explained by the destructive interference of waves with their counterparts totally reflected in an opposite phase. Therefore, in this case, LSPRs are able to couple to each other only via near-field, which strongly suppresses any collective phenomena. However, the feedback required for observation of Fano-Wood anomalies can emerge again from the coupling of LSPRs to an external photonic guided mode.
In this way, we observe coupling between strongly blue-shifted LSPRs (see appendix A) and guided photonic modes of a membrane silica waveguide (see Fig. 4.1 (c)). Since the lattice is located on the surface of the waveguide, guided modes of both parities are excited in the structure (see Fig. 4.1 (f,i,l)), unlike the previous case.
Reflection, transmission, absorption and diffraction spectra for the normal incidence of light are depicted in Fig. 4.4.
4.4 Optical modes in waveguides of various thickness
Although reflection, transmission and extinction spectra are convenient for the comparison with experimental measurements, absorption spectra are the most convenient for the theoretical analysis of the band structure. First of all, absorption directly indicates excitation of the plasmonic lattice, since it is the only source of dissipation. Secondly, guided modes of the waveguide lead to a resonant suppression of the effective polarizability, , and subsequent disappearance of absorption. Thus, absorption spectra allow us to observe resonant modes and see their plasmonic nature.
In this section, we consider the same lattice of silver nanodisks inside a silica membrane waveguide but vary its thickness. First of all, we compare spectra of a lattice in a bulk silica (infinite thickness) with a case of very wide waveguides ( and ). As can be seen from Fig. 4.5 spectra in both and polarizations are very similar and does not have significant discrepancies. In a waveguide finite thickness results in very dense lines of the waveguide and Fabry-Perot modes. Complicated patterns observed in the spectra, including characteristic concentric ellipses, are a result of Moire fringes originating from the crossing of very high- resonances and computational grid. After the blurring, which is peculiar to any experiment, all these fringes will disappear and we will see the spectrum very similar to the case of bulk silica.
When the waveguide has the thickness of Fabry-Perot and waveguide resonances become much more sparse and pronounced. They ’cut’ both localized and lattice plasmon resonances, but the general picture does not change.
Next, we consider and transition between ultra-wide waveguides and conventional ones. In Fig. 4.6 extinctions spectra of and waveguides are presented. Now, there is a countable, but still a large number of resonances and their hybridizations. Lattice plasmon resonance is still distinguishable, whereas localized resonance is strongly coupled to plenty of guided modes.
For relatively thin waveguides ( and , see Fig. 4.7) several other features emerge. First of all, it is seen that Rayleigh anomalies and resonances associated with them slightly move upwards because of a slight increase of an average permittivity. Another feature, which becomes observable is that hybridized plasmon-photonic resonances are located below the eigenmodes of the waveguide, when (red zone), whereas in the blue zone of the LSPR () they swap places. It is another demonstration of well known resonances’ splitting.
Finally, we observe typical waveguides of and thickness (see Fig. 4.8). The main feature, which is observed immediately is that lattice resonance became almost dispersionless. It is explained by the fact, that the action of Rayleigh anomaly is strongly suppressed by the existence of two interfaces in close proximity of a lattice. That is why plasmonic particles are coupled to each other only via near field and waveguide modes and as a result, there is no conventional lattice plasmon resonance and only localized resonance, as well as resonances hybridized with photonic guided modes, are observed. An interesting deflection of LSPR on an air Rayleigh anomaly is observed in polarization. This resonance can be described according to the paper [65].
Chapter 5 Comparison with experimental papers
Theoretical calculations are very interesting by themselves. However, the final aim of any calculation is to explain experimental results or predict them. For this reason, here we compare our calculations with two experimental papers.
5.1 Homogeneous environment
The first article, which we observe [62], considers lattices with different unit cells in an optically homogeneous medium. Experimentally such lattices are fabricated on glass substrates () covered by 2 nm titanium as an adhesive layer. Lattices themselves are constituted by silver nanodisks of (nm) diameter and (nm) height and are embedded into the index-matching oil (see Fig. 5.1,5.2).
Angle-resolved extinction spectra of both square and hexagonal lattices (see Fig. 5.1,5.2) show Rayleigh anomalies, which correspond to dips in the extinction spectra. Slightly below them, lattice plasmon resonances are observed. Generally speaking, both the structure and its spectra are very similar to the case considered in the previous chapter and does not require any additional analysis.
In our calculations, we did not take into account the adhesive layer. Although, titanium is very dissipative material and even such a thin layer can strongly modify LSPR it can not qualitatively change spectra and their inherent features. In this way, experimental and theoretical spectra turn out to be very similar (see Fig. 5.1). All the resonances and peculiarities are observed at the same frequencies. The discrepancy in magnitude of extinction is probably associated with non-ideal periodicity, the difference between the particles or diffraction in an experiment.
5.2 Lattice on a waveguide
As it was discussed previously, although close proximity of an interface destroys conventional lattice plasmon resonances, a similar coupling can be provided not by the photons propagating along the lattice in a homogeneous medium, but by a photonic guided mode. Such an effect for golden nanoparticles on an ITO waveguide was experimentally studied in a paper [66] (see Fig. 5.3 top-left panel). Two branches of Fano-Wood resonances are observed in an extinction spectrum of such a structure for different angles of incidence (see Fig. 5.3 top-right panel). Experimental and theoretical peaks and dips match but also differ in amplitude, which most likely is a result of non-ideality of an experiment. However, the dispersion of dips extracted from the experiment also coincides with theoretical estimations (see Fig. 5.3 bottom panel). It is also worth to mention, that energy splitting and so-called bound states in a continuum are observed in this case.
Chapter 6 Convergence and accuracy
6.1 Convergence
In order to compare convergence rate with conventional RCWA method and RCWA enhanced with adaptive spatial resolution (ASR) [27], we have considered extinction emerging in the lattice in a homogeneous silica under the normal-incident light of eV energy. It is seen in Fig. 6.1 that DDA enhanced method converges almost immediately, which is actually determined by the convergence rate of the dynamical interaction constant (see appendix B for details). Concurrently, original RCWA method converges very slowly and does not provide reliable results even for harmonics, which is the maximum available value in our calculations. RCWA+ASR calculations conducted by prof. Thomas Weiss converges much faster, however does not reach the performance of dipole approximation.
6.2 Accuracy
In order to verify the considered method, we compare our results with FEM calculations conducted in COMSOL Multiphysics. Since the calculation of extinction coefficients for each and every frequency and in-plain wavevector component takes of the order of minute it is possible to conduct computations varying only one parameter in a reasonable time. We have considered spectra corresponding to a normal incidence of light on the same structures. As it can be seen from Figs. 4.1 (j–l), our results almost perfectly match with FEM calculations.
However, the accuracy of results might strongly depend on the structure. We know, that our main approximation is dipole one. The smaller is the characteristic size of the particle the more valid is the method. In this way, the ratio of the particle size and other spatial parameters, such as period, wavelength and distance to the closest interface, act as small parameters.
Here, we conduct a brief study of the accuracy dependence on these parameters in order to understand roughly limits of its applicability. We leave out of account study of the dependence on wavelength for several reasons. Plasmonic structures, which we consider, have strong dispersion and support bright resonance, which makes it hard to distinguish the impact of real physical effects and the drawbacks of approximation. For example, low precision can remain unnoticed, when the particle itself has an extremely small cross section. Also, the issue of the wavelength dependence refers not to lattice effects, but to a single particle and its own peculiarities. Therefore, it is better to leave this question for the study of the properties of specific particles, but not the method.
6.2.1 Dense lattice
If we consider near field of a particle then we can qualitatively assume that the wavelength is very large and we deal with the electrostatic problem. For such type of problems, it is well-known that quadrupolar and other high multipole fields decay in space much faster then dipolar one. Therefore, we can utilize dipole approximation for relatively large periods, which are smaller than a wavelength. However, when we approach the close proximity of the particle high multipole fields become comparable to a dipole contribution.
In order to estimate the minimum distance on which dipolar approximation is still valid, we have considered lattice of silver nanodisks in bulk silica and on the air/silica interface for normally incident light of (eV) energy. Their dependences of reflection, transmission, absorption, and diffraction on the period are depicted in Fig. 6.2. In both cases, a very good match is observed for almost the whole range of periods. Significant, but still not crucial deviations arise only for a period less than (nm). Considering the fact that the diameter of silver nanodisks is (nm) this case corresponds to a very dense lattice.
6.2.2 Lattice near an interface
When we consider plasmonic lattice near an interface there are two possible approaches to conduct calculations. The first one is to include an interface in the local layer. In this approach, we account for reflection from an interface in the calculation of effective polarizability and local scattering matrix. Polarizability of a single particle is calculated near the interface as well. The main advantage of this method is the accuracy since we account for the shape of the particle and its influence on the near-field in close proximity of an interface. The main drawback is the necessity to recalculate polarizability of a particle for each distance from an interface.
Another approach is to consider plasmonic lattice in a homogeneous medium and couple it to an external interface via a combination of scattering matrices in the spirit of RCWA. This method is more convenient since we can easily change the distance between the lattice and the interface. It is accurate for large distances, since in this case interface indeed <<sees>> particles as dipoles. However, when we bring them closer, the interface should know something about the shape of the particle even though we work in the dipole approximation in the end. Local field of the particle is complicated and can not be described in dipole approximation. Another problem is the rapid increase in the required number of harmonics with approaching an interface.
From this perspective, we understand that for large distances both methods are acceptable, but the second one is much more convenient. At the same time for small distances, the first approach is not only faster than the second one, but the latter is simply incorrect. However, these speculations do not give us any estimations on the transitional distance at which we have to switch from one approach to another.
In order to reveal this issue, we compare the results of two approaches for different distances. We have chosen dense square lattice of (nm) period in order to avoid lattice resonances and effects and put them into two environments depicted in Fig. 6.3. In Fig. 6.4 we compare absorption spectra under the normal incidence for these cases. Left column corresponds to the lattice embedded into the silica and right column to the lattice hanging in the air. In both cases, structures are illuminated from the same side of the interface from which they are located.
As it is seen from panels (a) and (b) of Fig. 6.4 both methods give the same results for relatively large . Frankly speaking, in this case, the polarizability of a particle in a homogeneous space and near the interface are almost the same. Therefore, results coincide not because of the match of the contribution of the interface calculated via different approach, but because the fact that in both cases this contribution is negligible. At the same time, FEM calculations generally match our DDA+RCWA approaches. The most significant deviation is observed for the side resonances. Most probably there are two main reasons for that. Firstly, this resonance is associated with highly-confined edge plasmons and the required dense mesh was not provided in 3D FEM calculations. Secondly, the electric field of this resonance is very heterogeneous, which makes dipole approximation less valid. Indeed, hot spots on the edges make a huge contribution in polarizability, but when we account for their interaction we should probably take the distance between the edges ( nm), but not between the centers ( nm). In other words, for such dense lattice and complex resonances accounting for higher multipole harmonics is required for better accuracy.
For slightly smaller distance (panels (c,d)) results of two approaches no longer match. Moreover, as we predicted the 1st method, which includes an interface into the local layer provide more precise data. When the lattice is further closer to the boundary 2nd method totally breaks down. However, the first one works much worse as well - amplitude of the main resonance no longer matches FEM results. I do not have a clear explanation for this phenomenon, but the trends observed in Fig. 6.4 allow us to make several assumptions. We see that that the bias of a lattice even on several nanometers leads to significant deviations in the spectrum. However, when we bring a particle closer to the interface fields inside it becomes less homogeneous. From this perspective, it might be more correct to put the ideal dipole not in the geometric center of the particle, but in a slightly deviated point, which will change the spectra. In other words, we return to the discussion of the validity of the dipole approximation. Most probably accuracy is lowering near an interface because of the increase of the gradients in the field and corresponding higher order multipole moments.
At the same time, we should remember, that the period of the lattice is extremely small and is comparable with the size of the particles. Therefore, we expect that more sparse structures should provide much better results. Summing up, the obtained results confirm that we can calmly exclude interface from the local layer if it is removed from the lattice by at least several dozens of nanometers.
Chapter 7 Conclusions
We have proposed the method for the implementation of the DDA method in RCWA. Although our approach is very general and allows considering lattices of nanoparticles in a complex environment, in this thesis, we have implemented it for lattices in a homogeneous ambiance and on a boundary between two homogeneous media, which are the most practical cases. We have demonstrated its operational feasibility by calculation of spectra for solitary plasmonic lattice and a lattice coupled with an optical waveguide. An occurrence of Fano-Wood anomalies, lattice plasmon resonances and strong coupling between them is observed via the proposed approach. The accuracy of our approach was verified by comparison with FEM calculation, whereas the convergence rate was shown to be much faster than in conventional and adaptive-spatial-resolution enhanced RCWA calculations. Also, limits of applicability of the method, determined by the dipole approximation, were estimated. Comparison of calculations with experimental results shows a good match, which can be considered as an additional confirmation.
Fast speed of calculation (typically several dozens of milliseconds for one computation at a fixed and ) and high accuracy of results makes this approach convenient and perspective for both interpretation of experimental results and fundamental analysis of phenomena occurring in plasmonic lattices. Although we formulate our method for a single particle in the unit cell, it can be easily generalized for the case of multiple particles, which is to be done in further research.
Appendix A Polarizability tensor calculation
A.1 Total/scattering field formulation
A.1.1 General formulation
In order to find polarizability, , of an isolated nanoparticle, it is very convenient to switch to total/scattered field formulation. Let us consider a certain background structure, which does not include the particle itself and is defined by a background permittivity . We suppose, that in the absence of a particle structure is rather simple and solution of the Maxwell equations can be obtained analytically. To give an example can represent a plane wave in a homogeneous medium or Fresnel reflection from an interface between two homogeneous media. Since is a solution of the Maxwell equations, it fulfills one of them in the following way:
[TABLE]
[TABLE]
Now, we will consider the total structure, including the particle itself, defined by permittivity. In a similar way, we define total field as a solution of the following equation:
[TABLE]
In contradistinction to the background field, total field typically can be determined only numerically. For a practical calculations it is very convenient to express total field as a sum of background and scattered field . In this way, we can reformulate the problem to a calculation of scattered field. Subtraction of Eqn. A.2 from Eqn. A.3 gives us a rule for its calculation:
[TABLE]
[TABLE]
where is a perturbation in permittivities introduced by a particle. As we can see from Eqns. (A.1, A.5), scattered field can be found as a field radiated by a current induced by a background electric field [67, 68]. It should be especially noted, that this current is distributed over the volume of a particle and that radiation is considered in the environment defined by . Therefore, if a particle supports a certain eigenmode (e.g. LSPR) then a Purcell effect will be observed, which leads to a resonance in scattered field. In this way, a connection between emission in resonant mode and scattering of resonator becomes obvious. Thus, in general case the problem of scattering can be considered as an emission problem, which is easier for an implementation by any near-field calculation method.
We have found out a convenient way to determine electromagnetic fields in the process of scattering, but according to the DDA, our final aim is the determination of a dipole moment, which can be inserted in a background environment in order to substitute a real particle. Eqn. A.5 can be easily reorganized for these purposes:
[TABLE]
which means that a current is excited in a particle [69].
All the calculations above are very general and can be applied in a very wide range of problems. However, small particles are considered in this paper and therefore, we are only interested in a dipole moment of the induced polarization density. It can be easily found by integrating over a volume of a particle:
[TABLE]
Since, in this paper we consider small particles in media with constant permittivity, , over their volume (boundaries between materials does not intersect particles), we assume (not ) to be constant in space on the dimensions of a particle . This means that the primary induced current is just a vector with components, not dependent on coordinates and makes it possible to introduce as a simple tensor, connecting additional dipole moment of a particle, , with a background field :
[TABLE]
where is the volume of a particle and is a scattered field (determined as a field generated by the current density ) averaged over the volume of a particle. Such an approach clearly shows that dipole approximation can be applied even for a particle in a very complicated environment and provides an easy-to-apply procedure of tensor calculation for the cases considered in this paper.
Thus, the only thing that should be done to calculate polarizability tensor, , is the calculation of scattered field, which is a field radiated by an effective current density . Since, the connection between dipole moment and electric field is linear a unitary background field (e.g. ) is taken in numerical simulations. In general case, independent calculations are required to determine the response of a particle on each polarization of a background field. However, in the presence of additional symmetry, there can be only 2 or even 1 independent components, which simplifies calculations additionally. We conduct the described calculations in COMSOL Multiphysics, whereas they might be potentially realized via any near-field calculation methods.
A.1.2 Axially symmetrical case
In case when the structure has additional axial symmetry, the calculation can be significantly simplified. Indeed, for such structures electric field can be expanded as a sum over angular harmonics:
[TABLE]
Each harmonic is absolutely independent on the others and can be determined in a 2D calculation, which is extremely faster than common 3D ones. Moreover, for the calculation of polarizability tensor we need only harmonics, which strongly simplifies the general problem.
Indeed, the presence of axial symmetry means that the axis is one of the main axes and the two latter are degenerate and lay in the plane. For the sake of simplicity, we can choose and axes as the main ones.
As it is seen from the Eqn. A.8 component of the tensor corresponds to the averaged value of field over the volume of a particle. Obviously, only harmonic provides a non-zero value of the integral. In such a way, the computation of component of polarizability tensor almost does not differ from the general 3D case. We just excite -directed current in a simulation of harmonic and then compute the corresponding averaged field.
However, the computation of in-plain components of the tensor is slightly more difficult. Only harmonics make non-zero contribution in and components of averaged fields. However, each harmonic corresponds to a certain circular polarization of electromagnetic fields and therefore we can not just excite pure -polarized current. Without loss of generality we can consider harmonic and excite corresponding circularly polarized current:
[TABLE]
which radiates corresponding electric field . Since the polarizability tensor is diagonal then we can consider and components independently. In other words, we can say that component of current radiates electric field, which has only component of averaged field and the same with components. -polarized component of the current, which can be easily found as:
[TABLE]
Similarly averaged component of corresponding field can be determined as follows:
[TABLE]
where averaging means the following:
[TABLE]
A.2 Example of calculations
In Fig. A.1 we show energy dependence of in-plain component of a polarizability tensor of the silver nanodisk in bulk silica, air and in the air on the air/silica interface. It is seen, that wide resonances of LSPRs correlate with the obtained spectra (see Fig. 4.1 (d–i)). Since silica has permittivity higher than air, corresponding plasmonic resonance in a disc is strongly red-shifted. Resonance of a particle laying on a boundary is obviously located between them. Almost perfect matching of extinction spectra with FEM calculations (see Fig. 4.1 (j–l)) indirectly proves that our approach provides precise results and substantiates our choice.
In this way, there are several other important advantages of the method applied for the calculation of tensor. First of all, it is very universal so that allows considering particles of any, even very complex shape, in any environment. Secondly, it provides high accuracy, which as a result leads to high-precision spectra.
A.3 Resonant approximation for polarizability of indivisual particle
A.3.1 General formulation
The proposed method of calculation of polarizability tensor is accurate and fast enough. However, the most interesting objects for consideration are resonant nanoparticles including plasmonic ones and they allow to simplify the computational problem additionally.
Indeed, theory of radiation in low- resonators [70, 68, 45] states that in our case, scattering field can be expressed as a sum of the fields of quasinormal modes with definite frequency-dependent amplitudes . In most cases, when resonances’ spectra do not intersect or they have orthogonal polarizations, then the amplitude of each mode can be approximately determined independently on the others:
[TABLE]
where the field of the mode is normalized as follows:
[TABLE]
There are several methods to normalize quasinormal modes, which exponentially diverge at infinity [45]. However, in this work, we apply a convenient method described in the paper [68]. The main idea is that if we excite a certain mode with a current on a complex frequency very close to the eigenfrequency , then the total field is very close to the field of this mode:
[TABLE]
Therefore
[TABLE]
and finally we obtain:
[TABLE]
The pole itself is found by an iterative search via the Pade approximation according to the paper. Another important issue, which raises during the search of the eigenmodes, is the extension of permittivity to the plane of complex frequencies. One way to tackle this problem is to fit the data for real frequencies by some analytical model (for instance, Drude-Lorentz one) and apply it for the calculation on a complex frequency. Although this approach is very powerful it has several drawbacks as well. Sometimes it is hard to fit the data onto a wide spectral range, in some cases the choice of the fitting model is not evident and not unique. That is why in our calculations we apply more straightforward approach - we just employ Taylor decomposition. Since the typical -factor of plasmonic resonance is approximately , we estimate first and second derivatives via the divided differences on a spectral range of approximately width. Finally, Taylor series together with Cauchy–Riemann equations give us an estimation on permittivity at the complex frequency.
A.3.2 Polarizability tensor calculation
Accounting for the fact that background current, as well as background field, are constant over the particle in dipole approximation, we simplify the expression for the amplitude of the mode:
[TABLE]
Therefore, dipole moment of a particle can be expressed as follows:
[TABLE]
Finally, polarizability tensor can be found as follows:
[TABLE]
This is a general expression, which can potentially deal with an arbitrary number of modes. However, in most cases, small particles have only 3 pronounced dipole resonances. Moreover, in the presence of symmetry, they are degenerate. For instance, if all the 3 modes are symmetric, then:
[TABLE]
where is an -polarized mode. If modes are degenerate only in plain, then:
[TABLE]
where and are and -polarized modes correspondingly.
A.3.3 Axially symmetrical case
As it was discussed for the straightforward calculations, axially symmetrical structures can be treated by much simpler simulations. Calculation of modes in symmetrical resonators is easier as well. Resonance along axis corresponds to axial harmonic whereas in-plain resonances along and axis are superpositions of circularly-polarized modes with .
Circularly polarized modes have the following fields:
[TABLE]
This means that if the current
[TABLE]
excites fields of corresponding angular dependences, then current
[TABLE]
will excite field
[TABLE]
which corresponds to the -polarized mode.
Thus, the integral required for the calculation of the normalized field of the mode can be found as follows:
[TABLE]
The field averaged over the volume is found correspondingly:
[TABLE]
A.4 Analytical approach
In some simple cases, polarizability of a particle can be easily found analytically. One of the most well-known examples is the calculation polarizability of an ellipsoid in a homogeneous medium, which is described in many books and papers (e.g. [32, 71, 39]). When the silver disk is embedded in homogeneous ambiance it can be approximated by an oblate spheroid with semiaxes of and in order to provide the same volume and axes ratio. In-plain polarizability of such ellipsoid can be found analytically in quasistatic approximation [32] as follows:
[TABLE]
where is the depolarization coefficient, which is determined by an asymmetry of an ellipsoid.
However, this approximation is valid only for very small particles and we apply so-called modified long-wavelength approximations (MLWA) [39, 71], which accounts for the first dynamical corrections to obtain more precise results:
[TABLE]
where is the wavevector in a medium.
A.5 Rotating tensor
In the case, when a particle does not have rotational symmetry (for instance, nanorod) an issue of particle rotation arise. Conventionally, polarizability, , is calculated in a basis of main axes and therefore should be recalculated for a fixed basis.
First, we consider only in-plain rotation. A particle is rotated on an angle clockwise when we are looking in the direction of propagation of -axis. In this way:
[TABLE]
Therefore:
[TABLE]
A.6 Comparison of different approaches
Optical properties of plasmonic nanodisks are mainly determined by its in-plain component and particularly localized plasmon resonance, which shapes the spectrum. Here, we compare calculated via different approaches and estimate their accuracy in such a way.
We start with the consideration of the same silver disks in homogeneous silica and air, which is the most simple case and might be tackled with any approach. As can be seen from Fig. A.2, MLWA model generally describes the localized plasmon resonance. Although it is possible to slightly adjust and semiaxes and make the match almost perfect, it does not make sense. Indeed, if we need to conduct FEM calculations in order to find parameters for an analytical model, then this model does not provide any benefits not in time consumption nor in accuracy.
As can be seen from Fig. A.2 (c,d) there are several side resonances, which are too complex to be described by any analytical model. At the same time, the theory of quasinormal modes is developed specifically to work with them and therefore provides us such an opportunity. Each peculiarity in the polarizability is associated with corresponding plasmonic resonance. We have found several resonances and calculated their contribution to polarizability (see Fig. A.3). It is seen that the curve calculated in resonant approximation matches the precise data in a good way and follows all the peculiarities. The only significant deviation is observed at (eV). The reason for this phenomenon remains unveiled, but probably it is related to a significant change of permittivity on a width of the resonance. The main resonance ((eV), ) is a classical LSPR, whose field is almost homogeneous inside the metal (see Fig. A.4 (a)). Second and third resonances are slightly red-shifted and are much narrower ((eV), and (eV), correspondingly). Distribution of their electric fields (see Fig. A.4 (b,c)) shows that most probably they are formed as a hybridization of LSPR with so-called edge plasmons. Although this result is clearly observed theoretically and helps us to show the power of resonant approximation, in practice it will be almost impossible to observe this phenomenon. The edge rounding by a few nanometers fillet (which will obviously occur in the experiment) almost destroys these modes. The low and wide peak at (eV) constitutes from many closely-pitched resonances. It does not make sense to resolve all the modes and consider the peak in resonant approximation, therefore we just found several modes to confirm their existence and nature of the peak.
Silver disk in the air has very similar optical properties. The wonderful fact is that although resonances are strongly red-shifted and have completely different quality factors, distributions of electric fields are almost unchanged. They are so similar that it is impossible to notice any difference by eye.
A.7 Crossing the interface
As it was discussed above, two most commonly occurring cases are a particle in a homogeneous ambiance or placed onto an interface. The first one is the most simple for theoretical consideration and provides several unique phenomena, whereas the second one is the most natural for experimental realization. However, from the theoretical point of view, it is very interesting to observe the transition between these states and to study the dependence of the polarizability on the distance from the interface, at least on a certain example.
We have already considered silver nanodisks in silica and air environments, so here, we place the particle on different distances from their interface (see Fig. A.8). Figure A.9 shows how the polarizability tensor transforms while crossing the interface. The distance specified in the legend shows the spacing between the center of the disk and the interface. In this way lines related to (nm) distance correspond to the particles touching an interface from different sides.
The most important fact is that the particle moved away from the border by only 5 nanometers behaves almost as in a bulk medium. This means that polarizability is mainly determined by the local environment and therefore it is not possible to predict the polarizability of a particle on an interface, knowing only bulk polarizabilities, since it strongly depends on the shape of a certain particle.
The behavior of the side resonances depicted in Fig. A.9 (c-d) is too complex to be analyzed from this graph.
Appendix B Sum calculation
B.1 Green’s function filtering
First of all, since in this paper we consider only examples of environments, which have translational symmetry in the plane, and all the particles are located in the same plane, for a fixed coordinate, Green’s function depends only on the difference of coordinates. Therefore, sum 2.5 transforms to a slightly simpler form:
[TABLE]
where is the -th translational vector of a lattice in real space.
Typically this sum converges very slowly. For instance, in a homogeneous medium dyadic Green’s function (where is a wavevector in a medium) decays as , which makes the sum calculation not a trivial problem. Various methods have been developed for an efficient calculation of a lattice sum for homogeneous ambience [72, 73, 74, 75, 40]. Many of them, starting from the classical Ewald’s method are based on the following idea: if a lattice sum of a function converges slowly in real space, then this function can be represented as a sum of two auxiliary functions. The first one accounts for high gradients and decays very fast at infinity, whereas the second one, on the contrary, should be very smooth. In this way, the sum of the first function converges fast in real space and the sum of the second one can be calculated efficiently in reciprocal space (via the Poisson formula [76]).
Following this idea, we implement a very similar method, which helps us to obtain fast convergence for lattices both in homogeneous ambiance and on an interface between two media and has very clear physical sense. Green’s function has high gradients only in the proximity of the zero point, but in our case, this point is not included in the sum, since self-action of a particle is already accounted for by the polarizability tensor, . This fact gives us an opportunity to split the Green’s function in such a way that the contribution of the first, high-gradient function is negligible. In order to do that we multiply the original Green’s function by a special filter, which is equal to zero at the zero point and tends to unity for distances larger or equal to the period of a structure (see Fig. B.1). In this way the summation can be conducted over all the nodes of a lattice, having exactly the same result:
[TABLE]
where , is the filtered dyadic Green’s function in reciprocal space, is an area of a unit cell in real space and is a -th vector of reciprocal lattice. If the filtered Green’s function, , is a smooth enough that it has derivatives and all of them are absolutely integrable, then its Fourier image for , which provides fast asymptotic convergence.
This simple approach not only provides fast convergence but also has a clear physical meaning. Indeed, if we apply a filter, which is smooth, then each summand in reciprocal space has a sense of a contribution of corresponding diffraction order in the interaction between particles. Therefore, in most cases, it is possible to account for the same number of diffraction orders in sum calculation and in the scattering matrix , which is very convenient.
However, calculation of Fourier image of a filtered Green’s function is not always an easy task. When the homogeneous space is considered, then it is possible to compose such a filtering function, that Fourier transform can be calculated analytically. But, when a particle near an interface is considered, then expressions become too complex. Since in real space filtering is just a multiplication then in reciprocal space we deal with a convolution , where and . However, since filtering function tends to unity at infinity, it is convenient to introduce complementary filtering function, , for further derivations. In this way we introduce , which helps to find the originally filtered Green’s function as .
B.2 Green’s function near an interface
Here, we consider Green’s function for an emitter located in the proximity of an interface between two media since it might be easily reduced to a case of a homogeneous medium. It is known, that this function can be expressed as a sum of two parts, [77]. The first term is just Green’s function of the homogeneous medium and the second one accounts for the field reflected from the boundary, which is naturally calculated in reciprocal space because it requires just multiplication of amplitudes of plane waves by Fresnel coefficients.
Green’s function for a homogeneous ambience in reciprocal space can be expressed as a sum of - and - polarized contributions:
[TABLE]
where is the component of the wavevector (positive imaginary part should be chosen for ). Sign corresponds to upper/lower half-spaces, however, when we consider plane , the choice of the sign does not matter and all the alternating-sign components finally do not make any contribution.
Splitting of the tensor into two terms corresponding to different polarizations is very convenient for an accounting of reflection from an interface. We just have to multiply each term by corresponding reflection coefficient, choose appropriate signs keeping in mind the fact that the reflected wave propagates in the direction opposite to an incident one and take into account an additional phase, which is gained during the propagation to the interface and back.
Without loss of generality, we consider a lattice of particles placed above an interface, whose centers are located at a distance of from a boundary (see Fig. 3.2 (b)). The reflected part of Green’s function has the following form:
[TABLE]
where and are Fresnel reflection coefficients for corresponding polarizations, and is a phase delay.
B.3 Convenient representation of matrices
Although, formal expressions and representations for Green’s functions such as B.3,B.4 are very well known and can be used in any purposes, it is very useful to study them in a rigorous way. Understanding the structure of these tensors as well as their possible representations helps to deal with them efficiently both analytically and numerically.
We start from the definition of Green’s function for a homogeneous ambiance :
[TABLE]
The latter integral can be easily calculated:
[TABLE]
where for and for . It is important to understand, that in Eqn. B.5 emerges from the differentiation and therefore has different signs depending on the half-space, which is considered. It is the only place in the whole text in which we do not state the sign of .
We know, that electromagnetic wave, which propagates in a homogeneous medium are perpendicular, which means that , where is a vector of polarization (not necessarily unit). Therefore, for a fixed direction of wave propagation subspace of its polarization is two-dimensional, which means that it makes sense to switch -basis to a basis of some accessible polarizations. In this way, the following expression acts as the unity operator on an electric field (and as a consequence):
[TABLE]
where denotes the polarization. Application of this operator to matrix gives its elegant representation [78]:
[TABLE]
where is a polarization vector defined by polarization , in-plain component of a wavevector and sign defining the direction of propagation, is a reciprocal vector, which is orthogonal to vectors of another polarization and gives unity under scalar product with vector of the same polarization.
It should be noted, that any basis can be implemented in this formula. Most commonly used is the basis of - and -polarized waves.
[TABLE]
Explicit representation of B.8 sum in this basis results in well-known expression B.3. However, in some cases another basis might be more convenient. For instance, one is commonly used in RCWA:
[TABLE]
[TABLE]
Basis of left- and right- hand polarized waves is also convenient for some applications:
[TABLE]
[TABLE]
Now, we are ready to consider Green’s function of a dipole located near an interface. Its reflected part, , can be easily constructed in an analogous way. According to the previous section we consider a particle, which is located in half-space. According to the choice of the -axis direction depicted in several figures, this corresponds to a lattice beneath an interface, although it does not matter. In order to find reflected part, we should take, , tensor, multiply it by a phase factor, reflection coefficient and substitute polarization vector by a polarization vector of wave propagating in an opposite direction. This results in the following formula:
[TABLE]
where is a reflection coefficient from polarization to one. It should be noted, that each multiplier in Eqn. B.14 has very clear meaning and therefore it can be easily reconstructed "mnemonically" in similar situations. The rule is that reciprocal polarization vector accounts for emission and therefore corresponds to the direction of initially radiated wave ( in our case), whereas common polarization vector corresponds to the wave in its eventual state ( in our case). Since Eqn. B.14 is derived in rather general form it can be potentially implemented for any basis and reflection from a photonic crystal, which can even mix polarizations.
B.4 Analytical calculation of convolution for a homogeneous environment
We have already described the process of the filtering of Green’s function. However, when we deal with a case of homogeneous space, this general idea might be simplified. Indeed, operator acts locally in real space. Therefore, if some function has a plateau on which it is equal to [math], then obviously, this plateau will remain after an application of the operator. This gives us an opportunity to filter not a dyadic Green’s function itself, which has lots of components, but so-called scalar Green’s function of the vector potential, , which is an argument of the discussed operator.
This approach is much easier since the simple form of the internal part of Green’s function allows to find such a filtering function that has an analytical Fourier image. Together with the technique presented in the previous section, we understand that this leads to the fact that calculation of filtered dyadic Green’s function in reciprocal space is reduced to a substitution of the value of the integral B.6 with a specialized , which is found below.
Let be a filter as it was defined above. However, we should still remember, that now we filter another function.
[TABLE]
We have found a class of functions of the following shape that allows analytical calculation of the latter integral:
[TABLE]
where is a famous gamma function. Width of this function can be roughly estimated as unity:
[TABLE]
Therefore, filtering function of width can be found as . Finally, we return to :
[TABLE]
The wonderful fact is that the expression can be found as a solution of recursive equation (i.e. iteratively):
[TABLE]
which brings us to the final expression:
[TABLE]
As it was announced above, filtered Green’s function might be easily found as:
[TABLE]
The interesting fact is that even though the filter is a smooth function and provides a narrow space spectrum of the filtered vector potential, its space spectrum of an electric field is significantly wider because of specific peculiarities of the operator connecting vector potential with an electric field. This phenomenon results in necessity to take into account much more diffraction orders than it makes sense physically (typically several hundred instead of a dozen), however it is still much less than is needed to describe local field of plasmonic resonance (if we apply RCWA directly) and require just summation over them, which has a complexity of instead of standard matrix inversion in RCWA. At the same time, the opportunity to filter only one function instead of at least 6 and feasibility to do it analytically makes such an approach extremely efficient. The main advantage is that we need just to calculate analytically several terms of the sum, which is done extremely fast, much faster than any other required operations.
B.5 Fast calculation of convolution
2D convolution, which is required for filtering of is a rather expensive operation. In order to tackle this problem, we boost it in several ways.
It can be easily noticed, that components of of all the tensor summands have some trivial angular dependence. Here, we consider axial symmetric filtering functions , which means that also depends only on the absolute value of the wavevector in-plain component. This makes it possible to reduce the convolution to a much simpler form.
Indeed, let us consider a convolution for a certain component of the tensor . From the explicit expressions for tensors’ components it can be easily noticed, that any of them can be represented in the following way:
[TABLE]
where angular part takes one of the following forms:
[TABLE]
In this way convolution can be calculated as follows:
[TABLE]
For each and every element of a set a decomposition of can be conducted. To give an example, a case of is considered below. Since we obtain:
[TABLE]
where
[TABLE]
[TABLE]
[TABLE]
In this way, for a definite frequency, -dependence might be calculated on a grid and then interpolated. The fact that angular dependence is determined analytically and all pre-calculations are conducted for a 1D grid as well as the possibility not to calculate functions many times in closely-spaced points increase the speed of computations drastically. The strongest speedup is observed when the angle dependence of a spectrum is considered. Moreover, for some filters, angular part of integrals (B.29-B.31) and similar ones can be calculated analytically, which additionally speeds the calculations up.
Hereinafter, we observe a specific filter of the following shape (see Fig. B.3), which we use in our calculations. This filter is convenient for a practical utilization because of a high rate of decay and the possibility to operate with it analytically, which is demonstrated below. Its width can be roughly estimated as . Therefore, in order to obtain the filter of the provided width, (see Fig. B.1), we should use . Since original Green’s function diverges as at the zero point and filter decays as , then the filtered Green’s function has to have at least 2 derivatives and therefore for .
Fourier image of the filter can be easily calculated analytically:
[TABLE]
At the same time, Fourier image of the original filter is subsequently derived as
The following expressions are valid for the particular filter utilized in this paper:
[TABLE]
[TABLE]
[TABLE]
[TABLE]
[TABLE]
where , and are the modified Bessel functions of the first kind of corresponding orders, and .
To conclude, in order to find a filtered Green’s function in reciprocal space, , at any point, we should calculate just several functions on a 1D grid and continue them analytically to the whole plane. Pre-calculation of this auxiliary functions on a grid and subsequent interpolation makes these computations much faster. Moreover, for specific filtering functions calculation of the angular part of an integral can be conducted analytically, which boosts calculations additionally.
Finally, this approach, allows us to calculate angle-dependent spectra in a few minutes on a regular laptop, which is more than enough for practical utilization. Moreover, calculation of convolution is even not a bottleneck in our calculations and does not limit the performance of the whole program.
Appendix C Details on RCWA matrix calculations
C.1 Combination of scattering matrices
The elements of the scattering matrix which is a combination of two scattering matrices and denoted as is given by the following formula
[TABLE]
C.2 Method of oscillating currents
The RCWA formalism allows to calculate the emission of oscillating dipoles embedded in an arbitrary layer. The basic principle to do that is to construct the amplitude discontinuity vector which connects the vectors of amplitudes at coordinates infinitesimally above and below the dipole’s plane [43]:
[TABLE]
In this paper, we consider only examples, when the currents are located in a section of a homogeneous layer (which does not exclude the existence of other layers in near-field slightly higher and lower), therefore all the derivations below are conducted under this assumption. Application of the material matrix, , [14] to the equation C.5 gives us jumps of tangential components of electric and magnetic fields:
[TABLE]
where is the vector of discontinuities of Fourier components of electric and magnetic fields. The elements of this vector are found from the Fourier components of the surface current [43]:
[TABLE]
where and are defined by formulas (3.15).
In this way, in order to find the matrix , we should just express harmonics of this current through the vector . The first step is the application of the material matrix, which gives us fields:
[TABLE]
This vector of electric and magnetic fields harmonics is consequently transformed to the vector of Fourier harmonics of all the components of electric field:
[TABLE]
In turn, this vector is used for calculation of electric field in real space at the position of the -th particle:
[TABLE]
where hypervector is set by the following expression:
[TABLE]
where is the coordinate of a nanoparticle in the unit cell in a real space.
Thus, the current of the -th particle is determined as:
[TABLE]
Finally, according to the Poisson formula a grid of point currents in real space corresponds to the following harmonics in Fourier space:
[TABLE]
where is a surface of a unit cell in a real space and dagger denotes Hermitian conjugate.
Combining equations (C.5-C.13) we obtain the following expression for tensor:
[TABLE]
where
[TABLE]
Acknowledgements
This work was supported by the Russian Foundation for Basic Research (Grant No. 18-29-20032).
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] Ilia M Fradkin, Sergey A Dyakov, and Nikolay A Gippius. Fourier modal method for the description of nanoparticle lattices in the dipole approximation. Physical Review B , 99(7):075310, 2019.
- 2[2] Vladimir P Bykov. Spontaneous emission in a periodic structure. Soviet Journal of Experimental and Theoretical Physics , 35:269, 1972.
- 3[3] Eli Yablonovitch. Photonic band-gap structures. JOSA B , 10(2):283–295, 1993.
- 4[4] Kazuaki Sakoda. Optical properties of photonic crystals , volume 80. Springer Science & Business Media, 2004.
- 5[5] Robert Williams Wood. Xlii. on a remarkable case of uneven distribution of light in a diffraction grating spectrum. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science , 4(21):396–402, 1902.
- 6[6] John William Strutt. On the dynamical theory of gratings. Proc. R. Soc. Lond. A , 79(532):399–416, 1907.
- 7[7] Ugo Fano. The theory of anomalous diffraction gratings and of quasi-stationary waves on metallic surfaces (sommerfeld’s waves). JOSA , 31(3):213–222, 1941.
- 8[8] Ugo Fano. Effects of configuration interaction on intensities and phase shifts. Physical Review , 124(6):1866, 1961.