Multipole Excitations and Nonlocality in 1d Plasmonic Nanostructures

Efficient simulation methods for taking nonlocal effects in nanostructures into account have been developed, but they are usually computationally expensive or provide little insight into underlying physics. A multipolar expansion approach, among others, holds promise to properly describe electromagnetic interactions in complex nanosystems. Conventionally, the electric dipole dominates in plasmonic nanostructures, while higher order multipoles, especially the magnetic dipole, electric quadrupole, magnetic quadrupole, and electric octopole, can be responsible for many optical phenomena. The higher order multipoles not only result in specific optical resonances, but they are also involved in the cross-multipole coupling, thus giving rise to new effects. In this work, we introduce a simple yet accurate simulation modeling technique, based on the transfer-matrix method, to compute higher-order nonlocal corrections to the effective permittivity of 1d plasmonic periodic nanostructures. In particular, we show how to specify the material parameters and the arrangement of the nanolayers in order to maximize or minimize various nonlocal corrections. The obtained results provide a framework for guiding and interpreting experiments, as well as for designing metamaterials with desired dielectric and optical properties.


Introduction
The electric dipole approximation is commonly applied to study optical properties of various materials and structures. Generally, any nonmagnetic structure can be considered as an array of coupled discrete electric dipoles, while its electromagnetic properties can be explained in terms of the interference of the related electric dipole (ED) fields. This idea is at the heart of the well-known simulation technique, known as the coupled dipole (or discrete-dipole) approximation [1]. However, being computationally efficient, the coupled dipole approximation may provide no insight into the underlying physics. At the same time, electric and magnetic multipoles in nanomaterials are known to be responsible for a variety of optical phenomena, including artificial magnetism, additional extraordinary waves, negative refraction, giant electromagnetic chirality, enhanced spontaneous emission, superscattering, cloaking, etc.
Involving higher multipoles in materials and nanostructures for the description of their optical properties often becomes crucial to achieve desired functionalities [2]. As examples of more recent contributions to the field, the works on the generalized Kerker effect [3], topological Weyl metamaterials (MMs) [4], and Pancharatnam-Berry metasurfaces [5] can be mentioned. In fact, various nanostructures with a complex arrangement could give rise to enhanced or suppressed particular multipole excitations, offering new techniques to control their optical properties. Recently, a good review on multipole concepts with emphasis on achieving high directivities from radiating and scattering systems was given by Ziolkowski [6].
In homogeneous media, Raab and de Lange considered symmetry constraints for the constitutive relations, taking multipolar effects into account [7], and recapitulated and further developed the multipole theory [8,9]. They classified long-wavelength radiation effects that involve different multipole orders in the expansions of the displacement D and of the magnetic field H. Their hierarchy of the magnitudes of multipole contributions to the radiation-matter interaction is as follows: electric dipole electric quadrupole + magnetic dipole electric octopole + magnetic quadrupole and so on. The literature on the role of multipole excitations in plasmonic nanosystems is abundant. So, higher-order multipolar modes are known to interfere with the broad dipole mode giving rise to high-order Fano resonances with asymmetric profile [10]. Cho et al. [11] considered the nonlocal effective response of MM consisting of a pair of coupled metal bars using the formalism of the effective permittivity and effective permeability. As they found, the contribution of the electric quadrupole (EQ) to the effective permeability can be comparable to that of the magnetic dipole (MD). Alu and Engheta [12] dealt with infinite linear chains of plasmonic nanoparticles, taking their quadrupolar resonance into account. They showed that under specific conditions, the quadrupolar response of a plasmonic nanoparticle can be dominant. Besides, they concluded that despite the narrow bandwidth of the individual quadrupolar resonance, the overall bandwidth of chain guidance can be large due to strong coupling among the nanoparticles. Grahn et al. [13], dealing with arrays of nanoscatterers, introduced a multipole decomposition of the electric currents in the scatterers and connected it to the classical multipole expansion of the scattered field. They found, in particular, that completely different multipoles can produce indistinguishable electromagnetic fields. Silveirinha [14] derived the boundary conditions for the fields at interfaces of media using the so-called quadrupolar media approximation, which takes into consideration only electric dipolar, magnetic dipolar, and electric quadrupolar excitations. Wells et al. [15] analyzed the electric field enhancement in a plasmonic conical nanorod array with an account for second-order nonlocal corrections to its effective permittivity tensor.
Lim et al. [16] applied a multipolar decomposition to the analysis of natural modes in complimentary MM structures and showed that the magnetic quadrupole (MQ) can give rise to the linewidth narrowing of the Fano resonance. Dirdal et al. [17] showed that the electric octopole (EO) -MQ term can be of the same order of magnitude as the MD and EQ terms in MMs. Babicheva and Evlyukhin [18] considered optical properties of subwavelength spherical nanoparticles and their periodic arrays with an emphasis on the MQ response. In particular, they concluded that under certain conditions the effective polarizabilities of the ED and MQ, taking multipole coupling into account, can be significantly reduced. Within the multipolar framework, Mun et al. [19] introduced the higher-order dynamic polarizability tensors, which can accurately represent anisotropic meta-atoms with higher-order multipole transitions. Liu et al. [20] showed that in double-Dirac-cone PCs with the electric field along the optical axis, quadrupole resonances contribute to the effective permittivity. Rahimzadegan et al. [21] dealt with periodic metasurfaces (2d arrays) and derived expressions for their optical response using a multipole expansion for fields including dipolar, quadrupolar, and octopolar terms. As an illustration, they applied the developed approach to a fully diffracting metagrating, a polarization filter, and Huygens' metasurface.
One-dimensional plasmonic gratings attract the attention of physicists, as they can exhibit a number of interesting phenomena, such as, e.g., multiband absorption [22] or plasmonically induced transparency [23,24]. At the same time, despite a few efforts, the effect of multipole excitations on the nonlocal permittivity tensor remains poorly studied even in the simplest case of 1d MMs. So, the issue of the relationship between the unit cell geometry and the nonlocal response has not been addressed yet. In the present study, we take a step forward to assess the relative contributions of various nonlocal corrections to the nonlocal effective permittivity. As the corrections cannot be obtained in analytical form, our approach relies primarily on numerical simulations. In particular, it involves calculating the nonlocal permittivity with the use of the transfer matrix method (TMM), which is especially efficient for 1d nanostructures.

Geometry and Assumptions
The nanostructures under consideration are assumed to be infinite in the x-y plane and periodic in the z-direction (see Figure 1). This allows us to avoid dealing with the violation of Maxwell boundary conditions resulted from nonlocality [25]. Furthermore, the following assumptions were made. (1) We deal with nonmagnetic, homogeneous, isotropic dielectric, and plasmonic layers.
Although multipolar excitations can occur in all-dielectric nanostructures due to retardation effects, this usually takes place at sizes where the long-wavelength multipole expansion is no longer accurate. In contrast, in plasmonic multilayered nanostructures strong optical nonlocality can occur for moderate size unit cells as a result from excitation of surface plasmon polaritons (SPPs). Currently, plasmonic MMs offer rich opportunities to control the light intensity and polarization, as well as the near-field heat transfer and local density of states on the subwavelength scale [26,27]. (2) We neglect intrinsic nonlocality of constituent layers. It means that the vectors of the electric displacement and magnetic fields D and H at any point of space can be written in terms of the spatially averaged electric and induction fields, E and B, at the same point. It is necessary to keep in mind that in some cases, especially for thin metal layers and gaps, the local response model for their permittivity may become inaccurate [28][29][30].
To reduce the number of material parameters, we deal only with two-phase MMs. It implies that if the filling factors of the constituents 1 and 2 are f 1 and f 2 , respectively, it is always f 1 + f 2 = 1. At the same time, multiphase MMs with arbitrary number of different phases within unit cell can be studied using our technique. (4) For simplicity, only microgeometries with two-layer and four-layer unit cells are considered as shown in Figure 1. We note that any three-layer microgeometry cannot be realized when dealing with 1d two-phase MMs. (5) Herein, we limit ourselves to considering only the fundamental mode, which propagates along the interfaces (k z = 0) with the electric field, oriented along the optical axis (TM polarization). This means that we consider only one diagonal component ( ) of the nonlocal effective permittivity tensor˜ ij = diag( , ⊥ ). Thus, in what follows we implicitly assume that˜ ≡ . This choice is not accidental. It is motivated by the fact that this component can vary in wide limits, that gives rise to numerous potential applications [31][32][33][34].
Due to the simplifying assumptions made, the determination of the nonlocal effective permittivity becomes straightforward. In principle, this allows one to evaluate the nonlocal corrections of any order. In this study, we restrict ourselves by magnetic dipole and electric quadrupole, as well as magnetic quadrupole, electric octopole, and higher-order terms, up to the eighth order term in the wavenumber expansion.

Formalism and Techniques Used
According to the classical electrodynamics theory [35], for a nonlocal medium, the nonlocal permittivity (ω, k) can be introduced to relate the averaged electric field E av and the generalized electric displacement vector D g = E av + 4πP g via D g = (ω, k)E av , where P g is the generalized polarization density vector. In what follows, we assume a harmonic time dependence exp(iωt) for all the fields. The generalized polarization density can be written as where M is the magnetization vector. The polarization density P g at any point of space depends on the distribution of the macroscopic electric field in a neighborhood of this considered point. As a result, the vector of P g and the spatially averaged induced current contain the effect of the dipolar moment, as well as the effect of all higher-order multipolar moments. The starting point of our analysis is the power expansion of the nonlocal effective permittivity˜ (ω, k) in terms of the small parameter η (the period-to-wavelength ratio) as [36,37] (ω, k) = 0 (ω) + ηβ 1 (2) where 0 is the quasistatic (local) permittivity in the limit of k → 0, Equation (2) can be derived within the framework of the Landau-Lifshitz approach [38]. It assumes that the induced current density is described in terms of the electric polarization density, while the magnetization is zero. Furthermore, because of the symmetry,˜ (ω, k) =˜ (ω, −k), that results in the vanishing of the odd-order terms in the expansion (2) [39].
We deal with the fundamental mode, which propagates along the x-direction with the propagation constant k x = k 0ñ = k 0 √˜ , whereñ is the nonlocal effective index of refraction. This mode is in essence a coupled SPP excited at individual interfaces of plasmonic and dielectric layers [40]. The determination of˜ is one of the key points of our work.
Generally, the nonlocal effective permittivity can be found using purely numerical techniques, like rigorous coupled-wave analysis, modal method, finite element method, and finite-difference time domain method. Here, we invoke the analytical transfer matrix method [41]. The TMM is a remarkably powerful and versatile tool that can be adapted for studying various microstructures, but it is especially convenient, fruitful, and easy to operate when applied to layered systems, such as superlattices and multilayered waveguides. It can be also applied to finite and infinite periodic multilayers, as well as to 1d graded-index waveguides [42,43].
According to the TMM, if the z-component of the wave vector is zero, the dispersion relation for the the fundamental TM 0 mode reads where is the transfer-matrix of the ith layer, i and d i are its permittivity and thickness, respectively, x are the z-components of the wave vector in the corresponding layers, and n is the number of layers in the unit cell. From Equation (3), the nonlocal effective permittivity can be found as˜ Equation (3) is, generally speaking, the complex transcendental equation, which can be numerically solved using standard mathematical techniques for both the real (˜ ) and imaginary (˜ ) parts of the nonlocal effective permittivity.
In the particular case of n = 2, Equation (3) for the TM case reduces to the well-known Rytov equation [44] cos In its turn, Equation (5) can be split into two equations, and Equations (6) and (7) describe the so-called even and odd plasmonic modes, respectively, which correspond to the symmetric and antisymmetric distributions of the local electric field with respect to the middle of layers. The quasistatic term 0 for this case can be derived in the limit of d 1 , d 2 → 0. This term can be associated with the ED, and it is well known to be 0 = 1 2 where f 1 = d 1 /d and f 2 = d 2 /d. As η n → 0 in the limit of n → ∞, the expansion (2) may be truncated. Having determined˜ at several values of the lattice period d, this opens a possibility to determine unknown functions β i and B i (see Appendix A).
It is commonly accepted that the nonlocal corrections B i are associated with multipole excitations. Their relationship, however, is not straightforward.

Material Parameters
In what follows, we take the phase 1 as a transparent isotropic dielectric with the dispersionless permittivity 1 = 2.5 and the phase 2 as a plasmonic constituent material with the frequency-dependent permittivity For definiteness, we have set here ∞ = 4, the plasma frequency of ω p = 1.75 eV, and the collision rate (damping factor) of γ = 0.045 eV, that approximately corresponds to zinc oxide heavily doped with aluminum (2 wt%) [45]. On the one hand, such a choice of the Drude model parameters allows one to describe the plasmonic response of ZnO:Al in the telecommunication range of 1.5-2 µm, and on the other hand, if we are dealing with the normalized frequency ω/ω p and damping factor γ/ω p , it can be considered as typical for doped semiconductors. As to the permittivity of the dielectric phase, its specific choice is not a critical issue because the results of simulations depend on the ratio 2 / 1 , but not on the particular values of the permittivities 1 and 2 . Our focus on heavily doped semiconductors is motivated by the fact that, due to lower carrier density, they usually demonstrate a plasmonic resonance peak at lower frequencies (in the IR) as compared to typical metals. This opens up new possibilities for various photonic applications in the IR.

Two-Layer Microgeometry
Let us first consider the simplest case, when the unit cell consists of two layers. An analytical expression for the second-order term B 2 can be derived using an expansion tan(z) ≈ z + z 3 /3 for tan( k 1 d 1 2 ) and tan( k 2 d 2 2 ) in Equation (6). It is easy to make sure that it is of the form In the particular case of equal layers ( f 1 = f 2 = 1/2), Equation (10) reduces to The determination of higher-order terms involves cumbersome calculations. Nevertheless, in the case of f 1 = f 2 = 1/2, we have managed to derive an analytical expression for B 4 . It is of the form For further analysis, it is useful to study the behavior of the function B 2 . Qualitatively, its dependence on ( 2 ) is not critically important, and we take here for simplicity ( 2 ) = 0.3. The isolines of the absolute value of (B 2 ), (B 2 ) and B 2 (logarithmic scale) in coordinates of ( f 2 , ( 2 / 1 )), computed in accordance with Equation (10) are shown in Figure 2. The logarithmic scale allow us to clearly see all zeros and extrema of (B 2 ) and (B 2 ). As Figure 2 shows, the behavior of the frequency dependence of the second-order correction B 2 to the nonlocal effective permittivity, despite its simple analytical form, can be nontrivial. The behavior of the higher-order corrections can be even more complicated. Besides, it is necessary to keep in mind that as the plasmonic phase relative content f 2 becomes large, nonlocality can become strong. If so, the expansion (2) can become inaccurate, or at least more terms in Equation (2) should be retained.
In the following, we fix the period d at d = 150 nm, that yields η = d/λ approximately within the range of 0.05-0.1, and k 0 d approximately within the range of 0.3-0.6. In Figures 3 and 4, we display the nonlocal corrections η 2 B 2 , η 4 B 4 , η 6 B 6 , and η 8 B 8 for f 2 = 0.3 and f 2 = 0.5, respectively. The corresponding geometrical parameters are shown in Table 1.   It should be emphasized that the results, shown in Figures 3 and 4, do not depend significantly on the choice of the extraction methods, outlined in Appendix A. In addition, the results for the second-order correction B 2 and fourth-order correction B 4 , obtained with the use of the above numerical methods, are in good agreement with those obtained with the use of analytical calculations, Equations (11) and (12), respectively. This verifies both the consistency and high accuracy of the numerical procedures used.

Four-Layer Microgeometry
If the unit cell contains 4 layers, there are four parameters (say, d 1 , d 2 , d 3 and d 4 ), which specify microgeometry. If we fix d = d 1 + d 2 + d 3 + d 4 , three of them can be considered as independent. It is obvious that dealing with the full parameter space involves difficulties. Instead of this, we introduce a new geometrical parameter, µ, as This parameter takes its minimal value, µ = 0, if d 1 /d = d 2 /d = d 3 /d = d 4 /d = 1/4, i.e., for fully symmetric two-layer microgeometry, and its maximum value, µ = 0.625, if the relative thickness of any of the layers, d i /d → 1. Thus, in a sense, the parameter µ sets a measure of asymmetry. We are interested in how it can affect the nonlocal effective permittivity in general and various nonlocal corrections B i in particular.
For definiteness, in the paper, we restrict ourselves to the imaginary part of the nonlocal corrections. Doing so, it is necessary to keep in mind that the peaks of the real part are usually correspond to zeros of the imaginary part, and vice versa. For completeness, the corresponding real parts are given in Supplementary Materials, Figures S1-S3.
By analogy to the two-layer case, we take d = 150 nm. As one can expect, three options can be of interest: (1) f 2 = 1/2, (2) f 2 < 1/2, and (3) f 2 > 1/2. First, in Figure 5    Next, in Figure 6, we show the imaginary parts of the nonlocal corrections at f 1 = 0.7, f 2 = 0.3, computed for three different arrangements. All the geometrical parameters are listed in Table 3.  Table 3. The geometrical parameters which correspond to Figure 6.

Discussion
First of all, from Equation (8), one can conclude that ( 0 ) exhibits a Lorentzian profile that is symmetric with respect to the ED resonance, taking place at f 1 2 + f 2 1 0, while ( 0 ) is antisymmetric about this point. At the same time, despite the multipolar origin of all nonlocal corrections B i , the peaks and zeros in the spectra of B i (ω) cannot be unambiguously associated with specific multipolar excitations.
From Equation (10), one can see that the second-order term B 2 contains two summands, which can be associated with EQ and MD contributions. The first summand, which originates from the first term in the square brackets of Equation (10), disappears when f 1 = f 2 = 1/2, i.e., for a completely symmetrical arrangement. In contrast, the second summand, which originates from the second term in the square brackets of Equation (10), is always nonzero. Moreover, at f 1 = f 2 = 1/2, it takes its maximum value, given by Equation (11). We identify the first summand as associated with the EQ, while the second summand as associated with the MD. While usually the MD is characterized by an electric displacement loop, it is always present when there is enough retardation (nonuniformity) of the local electric field within the unit cell. As to the EQ moment, it is expected to be zero as a result of high symmetry at f 1 = f 2 = 1/2.
It is obvious that the first term in the square brackets of Equation (10) can be neglected close to the point of f 1 = f 2 = 1/2. In contrast, the first term prevails, while the second term can be neglected well below or well above this point. In the general case, of course, both terms should be taken account. This gives rise to an asymmetry of the B 2 (ω) profile.
Looking at Figure 3 ( f 1 = 0.7, f 2 = 0.3), we can conclude that (B i ) are symmetric, while (B i ) are antisymmetric with respect to the point of ω/ω p 0.4439, which corresponds to ( f 1 2 + f 2 1 ) 0, i.e., to the ED resonance. Near this resonant frequency, the second term in the square brackets of Equation (10) can be neglected. Neglecting the frequency dispersion of 2 (ω) in this range, it becomes clear why the second-order term B 2 exhibits the above symmetry or antisymmetry. Furthermore, the terms B 2 and B 4 , as well as B 6 and B 8 have opposite signs and hence they compensate each other. Because of the compensation effect,˜ 0 + η 2 B 2 + η 4 B 4 might be a more accurate approximation for this case than˜ 0 + η 2 B 2 + η 4 B 4 + η 6 B 6 . Ultimately, although the particular terms B i can be rather large, resulted nonlocality is moderate.
According to Figure 4 ( f 1 = f 2 = 0.5), the zeros of (B 2 ), (B 6 ), (B 8 ), as well as the minimum of (B 4 ) take place at approximately the same point, ω/ω p 0.394, that is close to the ED resonance, which occurs at 2 − 1 . The frequency dependence of (B 2 ) near this point is negative and reaches its minimum, in accordance with Equation (11).
(B 2 ), (B 6 ), (B 8 ) and (B 4 ) exhibit antisymmetry, while (B 2 ), (B 6 ), (B 8 ) and (B 4 ) exhibit symmetry with respect to the above point. Overall, nonlocality is relatively weak for this case; that can be related to the disappearance of the EQ contribution.
When f 2 is well above 1/2 (this case is not shown here), (B 2 ) is antisymmetric with respect to the ED resonance frequency; the position of its zero corresponds to a zero of ( 0 ) or the maximum of ( 0 ). At the same time, B 2 then can take large values and nonlocality is strong, that restricts the use of the expansion (2).
Let us now consider the four-layer case. At f 1 = f 2 = 0.5 ( Figure 5), the behavior of all B i is qualitatively similar to that of the two-layer case (Figure 4). At the same time, nonlocality is weaker as a whole. All nonlocal corrections rise with µ; moreover, higherorder terms B 4 , B 6 and B 8 rise faster than B 2 .
By analogy with the previous case, at f 1 = 0.7, f 2 = 0.3 (Figure 6), the behavior of the nonlocal corrections is qualitatively similar to that in Figure 4. The fast rise of the higher-order corrections with the parameter µ should be noticed.
Finally, we invoke Figure 7 to analyze the behavior of the nonlocal corrections B i at f 1 = 0.45, f 2 = 0.55. As can be seen, zeros of (B i ) take place at approximately the same frequency, ω/ω p 0.38. It is easy to check that this is close to the ED resonance. At this point, the contribution of the nonlocal corrections to (˜ ) is about zero. In contrast, their contribution to (˜ ) (not shown here) is maximal. All B i (ω) spectra become strongly asymmetric with respect to the ED resonance frequency. As one can expect, this is because the taken value of f 2 is close to 1/2; hence, both terms should be accounted for in Equation (10). Again, all the nonlocal corrections tend to rise with the parameter µ.
In general, the higher the term B i is, the more oscillations appear. For example, when (B 2 ) has only one zero, (B 4 ) has two ones, and (B 6 ) has five zeros. This means that more multipoles are coming into the game for every next nonlocal correction B i .
Overall, if the unit cell period d is fixed, nonlocality is weaker for four-layer microgeometries as compared to two-layer ones. This is consistent with the fact that all layers within the unit cell become thinner; hence, MM becomes "more homogeneous".
All the nonlocal terms B i could be controlled by tuning both material and geometrical parameters. For example, according to Equation (12), at f 1 = f 2 = 1/2, the B 4 term becomes especially large near the ED resonance in the case of high dielectric contrast. Currently, various techniques are known to tune the permittivity of constituents, say, using applied electric field [46], magnetic field [47], or temperature [48].

Conclusions
In this paper, we have presented a straightforward method to retrieve the nonlocal corrections to the effective permittivity of 1d plasmonic two-phase MMs with two-layer and four-layer unit cells under TM polarization. This opens up an avenue to analyze the spectrum of multipolar excitations in 1d MMs. Our study shows that MMs with such simple microgeometries are fraught with many surprises and offer many opportunities for tailoring their diverse functionalities. The key findings of our study can be formulated as follows: (1) When the dielectric phase content prevails ( f 1 > f 2 ), particular nonlocal terms B i can be large, but overall nonlocality can be moderate due to their mutual compensation. (2) When f 1 = f 2 = 1/2, symmetry of the B 4 term is opposite to that of other nonlocal contribution terms under consideration.
When the plasmonic phase content prevails ( f 1 < f 2 ), strong nonlocality can occur. (4) The spectra of multipolar excitations of the MMs with four-layer unit cell reproduce main features of MMs with two-layer unit cell, but nonlocality in the former case is weaker. (5) Higher asymmetry in four-layer arrangements enhances nonlocality.
This study provides information that could be useful for the design of nanomaterials and nanostructures with extraordinary optical properties. As we expect, it could have an impact on the quantitative modeling in typical applications and, in turn, could allow for developing nano-optical devices with unique functionalities.
Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/nano13081395/s1. Figure S1: The same as Figure 5, but for the corresponding real parts. Figure S2: The same as Figure 6, but for the corresponding real parts. Figure S3: The same as Figure 7, but for the corresponding real parts.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.