Dynamic matrix theory for resonance response of plasmonic metamaterial lattice

In this paper, a lattice dynamics method, named M-K matrix method, is proposed to investigate the near-field resonance response of a plasmonic metamaterial lattice under an oblique incident field with an arbitrary incident angle. By considering the electric, magnetic and field-dipole interactions, we construct a dissipative many-body Lagrange model for a reference lattice. A collective forced vibration equation, with the degree of freedom equals to the number of nanoparticles in a cell, is introduced to describe the lattice resonance under a polarized field. The resonance frequencies can be conveniently obtained from the poles of transfer function matrix. Based on this elegant matrix differential equation, one can calculate the amplitude-frequency and phase-frequency responses of plasmonic lattice, and analysis the normal modes from dispersion relations. The analytical results, which are from three examples: simple square lattice, binary chain and chessboard lattice, are perfectly matched with numerical simulations in a large frequency band, proving it to be an effective tool to calculate the dynamic response of plasmonic lattice.


Introduction
A metallic nanoparticle can restrict the optical field into a small space surrounding it, known as localized surface plasmon (LSP) resonance. The resonance frequency (RF) of LSP is not only sensitive to the background refractive index, but also determined by the particle's geometric and physical properties such as size, shape and electromagnetic parameters [1][2][3][4][5][6]. When a number of particles are arranged into an array, the LSP modes interact with each other, presenting fascinating properties significantly different from the individual particles, which is so-called lattice resonance [7][8][9][10][11][12][13]. Compared with isolated units, lattice resonance in close proximity exhibits a high-Q-factor extinction peak and supports several propagational modes (lattice waves). The optical devices based on periodic plasmonic array are therefore widely applied in nanophotonics, such as plasmonic sensing [9,11,12], plasmonic waveguides [13][14][15] and nano-gratings [16][17][18].
The measurement and calculation to the RF of LSP is a fundamental issue for both individual and periodic structures. For a single metallic nanoparticle with radius R and permittivity e, the near-field response of LSP under an applied field can be described by its electric potential function, which is calculated under quasi-static approximation [4], q pe e F = -+ ⋅ E r r p r cos 4 , For metallic nanoparticles with small loss of external optical field, the imaginary part of e e + 2 m is ignored. We see that w is independent of the incident and polarization directions of external field, as it presents a symmetric shape, indicating that RF is an inherent property. Moreover, it tells the RF of a nanoparticle is also independent of its radius R (even though it may affect the resonance amplitude). Generally, this theory is applicable when a particle's size is much smaller than the wavelength of light in the surrounding medium, that is,  l R. This is always satisfied in the sub-wavelength structures.
A direct way to enhance resonance is to place many units together to construct an oligomer or lattice. The near field coupling among these units will cause mode hybridization, and the corresponding LSP resonance is much more complicated than the case of isolated ones [19][20][21][22][23][24]. Recent studies also show that lattice with asymmetric dimer cell or oligomer structure can exhibit stronger resonance than that of symmetric ones [25][26][27][28][29]. These structures, with at least two units in a cell, cannot be simply treated as Bravais lattice. As they are ubiquitous in plasmonic systems, the main problem is, how to calculate the resonance response of a plasmonic lattice (periodic array) once a single nanoparticle's RF is known. Despite there have been already many experimental studies about plasmonic lattice, the analytical model for the coupled resonance modes excited by an external optical field in a plasmonic many-body system, to the best of our knowledge, has not been paid much attention.
As known, it is complicated to deal with a general infinite-degree-of-freedom oscillation equation analytically. If using numerical method, the computational cost is high and the truncation error is inevitable, and it makes little sense in analyzing the intrinsic property of resonance. Due to the particularity of plasmonic lattice, if a metallic nanoparticle can be regarded as a 'capacitor', one can reckon its 'capacitance' by the geometric parameters, and thus infer its 'inductance' from the RF of LSP (as the circuit model gives w = LC 1 [23,24,30]). As mentioned above, computing the RF of LSP is more convenient than that of general artificial metamaterials. Therefore, the circuit model can be applied to conversely evaluating the storage coefficients of plasmonic metamaterials [9,[31][32][33] once the individual RF of LSP is known. In this paper, we will show that, by considering the electric, magnetic and field-dipole interactions of all nanoparticles, an N-degree-of-freedom forced vibration equation is introduced to describe the collective resonance of a plasmonic metamaterial array. From this elegant matrix differential equation, the dynamic behavior of each particle in a cell, including its resonance amplitude and phase, can be tracked accurately. One can also find the RF from the poles of the inputoutput transfer function matrix, which is physically significant. The remarkable advantage of this theory is, although the mode hybridization for an infinite array is complicated, the coupled-mode equation can be reduced into a cell equation by applying the translation-invariant properties. The degree of freedom is reduced from infinite to N, equals to the units in a cell. Instead of analyzing each particle separately, the contribution from all lattice points can now be composed into a dynamic matrix. This will greatly reduce the computational cost and be beneficial to observe the property of normal mode and dynamic response. This paper is organized as follows: In section 2, we establish a dissipative Lagrange model for a translationinvariant plasmonic metamaterial array, and deduce the forced vibration solution of this many-body system and the eigenvalue problem about RF. In section 3, we present the dispersion relations of simple square lattice, binary chain and chessboard lattice, and analyze the property of normal modes. In section 4, we calculate the dynamic responses from the solution of forced vibration equation of these three lattices, and test its accuracy by comparing with the numerical results. Section 5 summarizes this paper. Figure 1 shows the schematics of resonance mode excited by an applied optical field in a plasmonic metamaterial lattice. According to circuit model, the coupling energy of nanoparticles comes from three parts: magnetic (current), electric (charge) and field-dipole couplings [23,24]. If we choose the charge variables as generalized coordinates, the current should be regarded as generalized velocity. Therefore, it is clear that the magnetic and electric energy can be respectively regarded as generalized kinetic and elastic potential energy, while the dipole energy in the external field should be regarded as generalized gravitational potential energy. The Lagrange model about the resonance of two-body systems has been investigated in previous studies [21,22,30]. While for a many-body system discussed here, it is more suitable to begin with a matrix form. We construct the Lagrangian for the electromagnetic coupling of nanoparticle array with an N-particle cell as(

Theoretical model
where s is nanoparticle index in cell l, going from 1 to N, represents the charge variables of particles 1 to N.  ss s / respectively mean the 'selfinductance' and the inverse of 'self-capacitance' of particle s.
The storage coefficients can be understood as counterparts of atom's mass and elastic force coefficients in crystal, and the effective dipole length relies on the geometric and physical parameters of nanoparticle. Therefore they are inherent property of lattice and possible to be measured, fitted from theoretical and experimental results, or calculated from numerical integrals [24,30]. For spherically symmetric (or concentric) particles, they are isotropic and independent of polarization of applied field, as mentioned in section 1. This is in contrast to some non-symmetric particles with polarization-dependent RFs and storage coefficients [23,25,30,31]. Here we only discuss the particles with a symmetric shape. In the following discussion, they are always assumed known as characteristic constants of a given lattice.
To describe a non-conservative dynamic system, we should further define the Rayleigh dissipation function, which is related to the sum of Joule heat of each nanoparticle [23], expressed as Because the structure is periodic, all of the cell blocks are identical, One can derive a matrix differential equation by substituting Figure 1. Schematic diagram of a plasmonic metamaterial lattice. As an example of 4-unit cell, the nanoparticles is numbered with indices s=1, 2, 3 and 4. The external optical field is produced with a polarization beam with incident angle θ. The lattice resonance mode is excited when the nanoparticle array is properly illuminated by the incident field.
equations (1) and (3) into disspative Lagrangian equation The derivative formulas of function matrix are seen in appendix)˜̈˜˜( (4) is an infinite-degree-of-freedom coupled-mode forced vibration equation, which describes the collective resonance of plasmonic lattice. Formally, it is similar to the second-order LRC circuit equation. Since the mode hybridization between different cells is complicated, seeking for the solution in a direct way is laborious. However, the translation-invariant property of the lattice indicates that the degree of freedom can be reduced. According to Bloch theorem, we assume the vector block ( ) is the in-plane projection of k, that is, for an arbitrary in-plane vector r, ⋅ = ⋅ k r k r, // and the ( ) Now, we have decoupled the associated vibration equation into an independent vibration equation at site = l 0. The resonance of each lattice point is determined by the in-plane projection of the wave vector. Obviously, once Q 0 k , // is solved, the whole resonance response can be determined. Note that the lattice is translation-invariant, the summation goes from -¥ to +¥. Therefore, we can simply replace ¢ l l with l. Equation In equation (7), we define the N-order total magnetic storage coefficient matrix ( ) S M k , // and the N-order total electric storage coefficient matrix ( ) both of which are periodic functions of the reciprocal lattice vector. Therefore, we only need to consider k // in the first Brillouin zone (FBZ).
Through these procedures, we have greatly reduced the difficulty of solving this many-body problem. Actually, as seen in equation (7), the matrices ( ) S M k // and ( ) S K k // have taken all of the coupling among in-cell ( = l 0) and out-of-cell ( ¹ l 0) nanoparticles into account. Moreover, the degree of freedom N is clearly presented in equation (7), and it equals to the number of units in a cell.Equation (7) can be analogous to an Ndegree-of-freedom second-order circuit equation. The lattice-wave solution of equation (7) is Matrix plays a role of transfer function matrix from the input (external field) to output (charge response). Principally, the collective RFs are derived from its poles. But for weak damping, as the approximation introduced in section 1, we can obtain the resonance condition by discarding the imaginary part (10) shows that, for non-degenerate cases, a single normal mode with wave vector k // is corresponding to N resonance sub-modes with RFs given from the singular points of N-order matrix ( ) ( ) w - // and the number of resonance sub-modes also equals to the number of the units in a cell. Equation (10) is equivalent to the following generalized eigenvalue problem Generally the eigenvalue of equation (11) is not always real and positive, and if not, the corresponding resonance mode cannot be excited. Fortunately, the self-storage coefficients are usually far greater than the mutual ones, and in such a case ( ) S M k // and ( ) S K k // are positively definite. Thus equation (11) has real eigenfrequency. Theoretically, N normal sub-modes can be solved corresponding to a fixed k .

//
Equations (9) and (11) are the main results of this paper. We name it as M-K matrix method. From equation (9), one can directly solve the dynamic response under an applied optical field and find the resonance (extinction) peaks. From equation (11), one can obtain the dispersion relation and normal modes. In the following two sections, we will test the feasibility and accuracy of M-K matrix method in lattice resonance problems.

Free vibration-normal modes
Principally, the free vibration will soon get eliminated because of the damping. But it is also worthwhile to analysis the free vibration modes since it can give the dispersion relation, which includes the dependence of RF on the wave vector. From dispersion diagrams, one can analysis the characters of resonance modes excited by external optical fields with different wavelengths and incident directions.  (11) are reduced to scalars. As a sub-wavelength structure, one can estimate the wavelength is larger than the size of nanoparticle but smaller than the lattice constant. The first and second nearest nanoparticles are most likely located in the near-field zone, in which the electric and magnetic interaction will soon decay when the distance increase from zero [24,30]. On the other hand, the neighbors over the third nearest sites can be regarded as the far-field zone, and the couplings are much weaker than those in near-field zone. Therefore, we only consider the first and second nearest couplings. Equation (7)  x y and and a is lattice constant. The eigenfrequency can be obtained from equation ( x y x y x y x y  We depict the dispersion relation of a plasmonic square lattice in figure 3. The collective RF can be continuously modulated by the direction and wavelength of the applied field, and it fluctuates around the individual RF w . 0 If the electric couplings are stronger than magnetic ones, the minimum RF appears under normal incidence. If reverse, the maximum RF appears here. The rotational and inverse symmetry of square lattice is consequentially reflected in the FBZ of its reciprocal lattice. That is, the RF in each mode satisfies ( ) ). The corresponding 8 symmetric points in FBZ share a same RF under a certain normal mode. The symmetry property of vibration normal mode is similar to that of the photonic band structure [34], despite they are different issues.

Binary chain
As an example of non-Bravais lattice, we consider another plasmonic lattice with two nanoparticles (a dimer) in a cell, forming a 1D binary chain, as shown in figure 4. The binary chain is somewhat like the Su-Schrieffer-Heeger model [35] and Rice-Mele model [36] which supports topological states. However, because the lattice constant of optical superlattice is much larger than that of solid lattice, here the quantum effect and topological states are generally not observable. But this does not prevent us to analyze its dispersion properties under classical dynamics.
According to equation (8), 2 Taking the nearestneighbor approximation, Then   (18), the dispersion curve is depicted in figure 5. A binary chain contains two normal modes without any intersections, meaning the modes are non-degenerate. Analogously to the diatomic chain of crystal, the lower and upper branches can be referred to acoustic and optical branches [37]. However, this is quite different from crystal lattice wave because the generalized kinetic energy in equation (1) also contains the contribution of 'mutual inductance', which never happens in crystal model. As seen, a remarkable consequence cause by this fact is that ( ) w ¹ 0 0for acoustic mode. where ( ) M r is a 4×4 matrix. The approximation cuts off the coupling terms beyond the second-nearest sites, and therefore the matrix elements satisfy the following relationship (seen from figure 6)  The rest are zero. Equation (19)∼ (23) give the expression of ( ) The same goes for ( ) Then the eigenvalue problem equation (11) becomes  are respectively defined as nearest magnetic and electric coupling coefficients (subscript 1), second-nearest ones of particles 1 (or 4) (subscript 2), and second-nearest ones of particle 2 (or 3) (subscript 3).  The 4 dispersion surfaces have no intersection, which means these 4 eigenmodes are non-degenerate. Apart from some individual eigenmodes, generally , , , , , .
x y y x This is different from the simple square lattice because C C ,

Amplitude-frequency and phase-frequency responses
The near-field resonance response of a plasmonic array can be calculated from equation (9), once ( ) S M k // and ( ) S K k // are known. The dynamic response can be clearly seen in the resonance curves, that is, the amplitude- frequency response (AFR) and phase-frequency response (PFR) curves. To test the accuracy of M-K matrix method, we also use the full-wave simulation results computed by lumerical FDTD solutions as a benchmark for comparison. Let us first estimate the self-storage coefficients for a single nanoparticle theoretically. Provided a spherical particle whose capacitance can be expressed as ( ) and therefore ( ) =C K 1 0 fF. For a micro-nano structure with sub-wavelength characters, the magnitude of RF approximately ranges from THz to PHz (this interval covers the frequency range from near-infrared region to near-ultraviolent region). Therefore the 'selfinductance' ( ) =L M 0 fH.The nanoparticle tends to be like a 'capacitor' rather than an 'inductor'. For a metallic nanorod, the estimation value is roughly the same.
According to the above estimation, the storage coefficients of particle 1 are set as / during the following calculation. One can obtain the individual RF is The resistance of particles 1 and 2 are respectively set as 0.06Ω and 0.04Ω.
These parameters can ensure a weak damping in equation (9). To compare the properties of different lattices, the dynamic responses are all calculated under a unit external field, and adjacent nanoparticles of complex lattice (in-cell or out-of-cell) are assumed isometric.

Simple square lattice
For simple square lattice, the storage coefficient matrices and generalized force vector in equation (9)  In the simulation results, we consider the interaction from all neighbors included those over third-nearest ones. While in the analytical results, the interaction is only included up to the second-nearest neighbors, same as that in section 3. The analytical result is well coincident with numerical one, proving the interaction energy of far field can be neglected.  (9) and depicted in figure 9. As predicted, there are two resonance peaks appearing in AFR curves. The lower and higher ones are referred to the acoustic and optical modes. The two particles share same RFs, while the oscillation encounters different levels of enhancement near the RFs. The phase also has a step at these points. The analytical and numerical results are also well matched.  figure 10. As expected, when illuminated with an external optical field, the chessboard lattice produces 4 resonance peaks, corresponding to the above-mentioned 4 eigenmodes. The amplitudes of 4 particles present different levels of enhancement near their common RFs (some resonance modes are not obvious). It is notable that the dynamic behaviors of particles 1 and 4 (2 and 3) present a high proximity, only phase difference. This is reasonable because they are identical particles and occupy the equivalent sites. The numerical result also verifies the accuracy of M-K matrix method. By comparing the maximal resonance amplitudes with a same intensity of external field, we see that complex lattice (chessboard) can achieve stronger resonance than simple one (square). If both lattices are complex, 2D planar structure shows better performance than 1D chain. These clear results may provide some reference for the experimental design of plasmonic lattice.

Conclusions
In this paper, the M-K dynamic matrix method is proposed to deal with the collective near-field resonance problem of plasmonic metamaterial lattice. By constructing a dissipative Lagrange model for a nanoparticle lattice, the resonance response is entirely described by an N-degree-of-freedom forced vibration equation, in which the interactions among nanoparticles are totally composed into the dynamic matrix. Based on this matrix differential equation, we presented the detailed calculation of the normal modes and dynamic responses of square lattice, binary chain and chessboard lattice. The three examples indicate that complex lattice can achieve stronger resonance than Bravais lattice, and 2D array is superior to 1D chain. The results of this analytical method excellently agree with the numerical simulation, and it is feasible and accurate in a large range of frequency band. We believe M-K matrix is a powerful method in treating resonance problems of oligomer plasmonic lattice [27][28][29]38], and may also provide a reference theory for relative researches on harvesting light and photons by using nanoparticle arrays [39,40].

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors. is a block matrix, and x kl is element of matrix X.