Generalized Elastodynamic Model for Nanophotonics

A self-consistent theory for the classical description of the interaction of light and matter at the nano-scale is presented, which takes into account spatial dispersion. Up to now, the Maxwell equations in nanostructured materials with spatial dispersion have been solved by the introduction of the so-called Additional Boundary Conditions. In this paper, we derive an approach where non-local effects are studied in a precise and uniquely defined way, thus allowing the treatment of all solid-solid interfaces (metals, semiconductors or insulators), as well as solid-vacuum interfaces in the same framework. The theory is based on the derivation of a potential energy for an ensemble of electrons in a given potential, where the deformation of the ensemble is treated as in a solid, including both shear and compressional deformations, instead of a fluid described only by a bulk compressibility like in the hydrodynamical approach. The derived classical equation of motion for the ensemble describes the deformation vector and the corresponding polarization vector as an elastodynamic field, including viscous forces, from which a generalized non-local constitutive equation for the dielectric constant is derived. Boundary conditions are identical to that of elastodynamics and they emerge in a natural way, without any physical hypothesis outside the current description, as it is commonly required in other non-local approaches. This description does not require the discontinuity of any component of the electric, magnetic or polarization fields and, consequently, no bounded currents or charges are present at the interface, which is a more suitable description from the microscopic point of view. It is shown that the method converges to the local boundary conditions in the low spatial dispersion limit for insulators and conductors, quantified by means of a parameter defined as the"characteristic length".


I. INTRODUCTION
The interaction of light and matter at the nanoscale has been a topic of intense research in recent years [1,2], due to the extraordinary advances in the manipulation capabilities of materials at the nanoscale [3][4][5]. The range of applications of this science is extremely broad and continues growing up [6], therefore theoretical and numerical tools for their accurate description are of primary interest [7][8][9].
The nanoscale is a complex size limit for the study of the interaction of light and matter, since typical structures are big enough to consider the problem from the classical point of view although some quantum effects can be observable. However, it is still possible to use a full classical description, as long as we find a proper constitutive equation relating the electric and magnetic fields with the induced currents and charges, which can take into account quantum corrections [10].
The simpler form of the constitutive equations used in electrodynamics are linear and local in both space and * dtorrent@uji.es time, and they define the dielectric permittivity, magnetic permeability and electrical conductivity [11]. Although nonlocality in time is commonly assumed, resulting in frequency-dependent constitutive parameters, spatial non-locality is in general left behind to more refined models of matter, since they become important only at the nanoscale [12][13][14][15]. It has to be pointed out that in the domain of metamaterials spatial dispersion has been a topic of intense research as well, since the distance between the "meta-atoms" is only one order of magnitude smaller than the operating wavelength of the field [16][17][18][19].
Despite the great success of the spatially local description at the macro or even micro scales, the theory fails in the accurate description of nanomaterials, since spatial dispersion becomes more relevant and it has to be included in the constitutive parameters [20][21][22][23][24][25].
When the constitutive equations become spatially nonlocal additional modes emerge in the solution of the wave equation, and the boundary conditions derived within the framework of Maxwell's equations are insufficient to match all the excited fields at an interface between two materials. This problem has derived in a countless number works proposing the so-called "additional boundary conditions", which are different depending on the dif-ferent response models [25][26][27][28], what in turn means that spatial dispersion is not generally treated in the same way in insulators, semiconductors or metals. Despite the fact that these additional boundary conditions are deduced with more or less reasonable physical assumptions, they are not derived within the framework of Maxwell's equations complemented with the constitutive equations, what means that the description is not "closed" and it can result in an inaccurate description of the electrodynamic problem, as can be seen for the fact that the problem still remains unsolved [29][30][31].
In this work we will show that, by means of an elastodynamic model of the induced current, we can find a self-consistent description of the light-matter interaction which accounts for spatial dispersion. We will show that when an ensemble of electrons is deformed its potential energy is identical to that of an elastic body, so that dynamically it moves as a continuous elastodynamic field [32]. In this context, we can develop a non-local theory for electrodynamics, where the required boundary conditions arise in a natural way. Moreover, we will define the local limit by means of a "characteristic length" parameter in which we recover the local description, reinforcing therefore the generality of this model.
The paper is organized as follows. After this introduction, in section II a discussion about the problem of non-locality and its possible solutions is presented, then section III presents the notions of quantum pressure and stress, and the elastic energy of the ensemble of electrons is derived. Section IV presents the elastodynamic formulation for nanophotonics, with a discussion about boundary conditions. Section V analyzes the solution of the wave equations in isotropic materials and section VI shows several numerical examples concerning the solidsolid and vacuum-solid interfaces. Finally, section VII summarizes the work.

II. LOCAL AND NON-LOCAL MATERIALS
The evolution of the electromagnetic field inside a material where no sources are present is described by means of Maxwell's equations [11], where J i and ρ i are the induced current and charge densities in the material, respectively. The continuity equation for the current density is implicit in Maxwell's equations, and it is derived taking the divergence of equation (2) and using equation (4), Maxwell's equations with the induced currents are not enough to solve the full electrodynamic problem, and we need a relationship between the induced current J i and the electromagnetic field. This relationship is complex to obtain, since it implies a many-body problem solved in the framework of classical, semiclassical or quantum theories, however a phenomenological approach is commonly considered in which the classical permittivity, conductivity and permeability are defined. Then, in classical electrodynamics, for non-magnetic and non-chiral materials, the induced current is expressed as with σ E and χ E being the electric conductivity and susceptibility, respectively. It is common to define the polarization vector P as so that, assuming a harmonic time dependence of the fields of the form exp(−iωt), equation (6) is which defines the well known complex dielectric constant ε by means of the dielectric displacement D, defined as so that If either σ E or χ E are discontinuous at an interface, so it is the current J i and, by means of equation (5), a surface charge appears. In classical electrodynamics these discontinuities and surface fields are usual, and they do not imply any non-physical description of the fields.
However, there are two points which make equation (6) unpleasant from a deeper physical insight. First, this local relationship in both time and space implies an instantaneous response of matter (locality in time) and a point-to-point response (locality in space), that is, the induced current in the material is instantaneous and it depends on the field at a given point only. Clearly this is not a true physical situation, since the response of the charges will have some inertial response and will be influenced by the surrounding material. Another point, less discussed in the literature, is that current discontinuities and surface charges are not really possible at the microscopic level, since charges and currents are quantum entities described by wave functions which are in general continuous across the interfaces. Then, while it is true that the classical description of electrodynamics assumes that these discontinuities are just idealizations, we would like to find a description in which these do not occur, so that the theory would allow us to reduce more and more the scale of validity.
In this work we will show how, surprisingly, a nonlocal version of equation (6) is possible which additionally implies the continuity of all the fields involved in the interaction.

A. Non-local constitutive equation
When elementary models of the light-matter interaction are considered, the derived constitutive parameters χ E and σ E are found to be frequency-dependent, but this dependence has no consequences in the nature of the fields, since equations (1) to (4) are usually solved assuming time-harmonic dependence.
However, the constitutive parameters can depend on the wavenumber as well, and this dependence is not equivalent to the frequency-dependence since, usually, we have to work on bounded materials. Then, let us assume for instance that the dielectric constant is both frequency and wavenumber dependent. We will have (for the purposes of this section we will consider only a scalar-like relation between fields, the full vectorial theory is developed later on) while we can continue working on the frequency-domain, if our material is bounded we need to work on the space domain, and Fourier transform the above equation, thus we have where we have used the convolution theorem. The above equation clarifies the term non-local for wavenumberdependent constitutive paramaters: the dielectric displacement at point x depends linearly on the electric field at point x , distant from x a length d = x − x . The dielectric function ε(d, ω), as a non-local response function, weights the contribution of electric fields at different x , and will decay as d increases. In the local limit, this function is proportional to the Dirac delta function δ(x − x ). It is obvious that, if we introduce this relationship of the dielectric displacement with the electric field in equations (1) to (4), even in the time-harmonic regime we obtain a more complicated equation. As the convolution theorem applies only to a spatially invariant medium, for a bounded material the dielectric constant will not be in general a function of d = x − x and we get the most general linear relation between two scalar fields.
which complicates even more the problem. In next subsection we will discuss how this problem could be solved in Fourier space by means of the so called "additional boundary conditions".

B. The need of additional boundary conditions
The usual method to solve the electrodynamic equations in bounded materials is by means of the application of boundary conditions, roughly speaking: knowing the solution of the fields in two materials, we apply boundary conditions at the interface between them and we find the solution of the problem. Maxwell's equations (1) to (4) provide us of the required number of boundary conditions, once the constitutive equation for the induced current has been found, but only if this constitutive equation is local in space, otherwise we require of additional boundary conditions, and these additional boundary conditions has been (and continues being) a topic of a great discussion in the literature. The origin of this need is found in the number of solutions that a wave-number dependent dielectric constant provides, as will be explained below.
Let us assume the problem of reflection and transmission at a flat interface. If the dielectric constant is frequency-dependent only, the dispersion relation in the material is typically of the form which means that, for a given frequency, we have two waves propagating through oposite directions (at normal incidence and in this scalar version, the vectorial case is more complex, as will be explained later). For a classical reflection and transmission problem we will have therefore the incident wave, the reflected wave propagating backwards and the transmitted field at the other material. We will need to determine the amplitude of the reflected and transmitted waves and Maxwell equations provide us indeed two boundary conditions: the problem is perfectly defined. However, if the dielectric constant is also wavenumberdependent, for a given frequency our dispersion relation is which, due to the dependence on k of ε, can give more than two modes propagating in opposite directions. Then, a given incident field will excite, in the simpler of the situations, two reflected and two transmitted modes (assuming both materials non-local).
Maxwell equations provides only two boundary conditions, but we have four modes to determine, so that the problem is not well defined. We need the so-called "additional boundary conditions" to completely solve the problem.
The amount of works about these additional boundary conditions is huge and actually experts have developed the acronym ABC to refer to specific examples or generally, to the very issue presented above. There is no consensus about which ones are the correct ones. It has been assumed then that the right boundary conditions depend on the micro-structure of the materials and the interfaces, so that there is not a unique solution for the problem of additional boundary conditions at the macroscopic level. However, we believe there is a well-defined macroscopic solution to this problem based on a right definition of the constitutive parameters, whose specific values of course depend on the micro-structure of the material but whose macroscopic behaviour can be universally defined, as we do in the electrodynamics of local materials, where obviously we don't apply the same boundary conditions for metals and for lossless dielectrics, but it is due to the fact that their constitutive parameters have different values. Next subsection will present the approach we will follow in this work to overcome the problem of the additional boundary conditions.

C. Solution from additional field equations
The additional boundary conditions have been discussed in uncountable ways, and these are suggested by more or less acceptable physical arguments or microscopical models. Our objective is to find a set of macroscopic equations from which obtain the required boundary conditions, with these equations being functions of a set of (local) constitutive parameters whose numerical values define the different type of materials. It is obvious that the additional boundary conditions cannot be obtained from Maxwell equations, and physical arguments can be enough for simple situations (vacuum-solid interface at normal incidence), but more general situations require of a more rigorous solution.
The motivation to search for the field equation of the induced current comes from equation (13), which is equivalent to a relationship between the polarization P and the electric field E of the form P (r, ω) = ε 0 χ E (r, r , ω)E(r , ω)dV.
The above equation suggests indeed that the function χ E (r, r , ω) is the Green's function of differential equation of the form and with L being a partial differential operator. If we are able to find the operator L we will be able to find the right boundary conditions in a similar way as we find them from Maxwell equations. It has to be pointed out that the objective of this work is not to model nor discuss the consequences of the frequency-dependence of the complex dielectric constant, our objective is to discuss the consequences of the wavenumber dependence, which will enter into the model through spatial derivatives in L, and how this affects boundary conditions. Then, it has to be assumed that all the constants appearing in the model equations could present a more or less complex frequency dependence, even more if our material is artificially nano-structured, but we will not take care of this dependence.

III. QUANTUM STRESS-STRAIN RELATION
The most elementary microscopic model of the dielectric constant assumes that the electron is bounded to the atom by means of a spring-like restoring force [11], so that when the electric field interacts with it the equation of motion is from which we can solve for the contribution to polarizability of one single electron as p = −er. If the material is a conductor the restoring force −κr is set to zero and we recover Drude's model for the free electron. While these models are quite elementary, they allow to explain some aspects of light-matter interaction, and in both models we recover a local and complex dielectric function.

A. Quantum pressure and quantum stress
Sommerfeld's model applies quantum statistics to the physics of the free electron in the solid, which is described as an ensemble of non-interacting spin 1/2 particles. It is found that an ensemble of N particles in a box of volume Ω has a total energy E given by where E F is Fermi's energy and is the maximum energy level occupied by the electrons. This energy is proportional to n 2/3 , with n = N/Ω, so that we can obtain the so-called quantum pressure P of the gas from the thermodynamical definition of pressure, and the compressibility or bulk modulus The electron gas is then described as a fluid material and the linearized equation of motion is that of the acoustic field subject to a body pressure due to the external electric field. This is the so-called hydrodynamic model for plasmonics [10] and, although widely used, we will show here that a more accurate description is needed, since the electron gas is found to have not only a bulk but also a shear modulus, so that it is better described as a solid material.
Other models also include in this equation Ohmnic losses or diffusion [33], what essentially changes the frequency dependence of the constants involved in the model. Including other sources of dissipation at the quantum level, inter or intra band transitions for example [34], might also change this frequency response, but for the purpose of this work it will be enough to assume that all the parameters appearing in the models are complex and frequency-dependent, since our objective is not to model the interaction, but to understand the role and nature of the additional spatial derivatives appearing in the equation of motion for the induced current.
The hydrodynamic model assumes that the gas of electrons is a fluid, and that only volumetric changes are possible. Moreover, it assumes as well that the internal restoring force is due to the gradient of a scalar pressure field. We will assume here a more general deformation of the ensemble of electrons to demonstrate that they behave actually as a solid.
Let us assume we have an ensemble of electrons in equilibrium in a periodic potential V (r). Let us assume now that some external perturbation (like an electric field) modifies this wave function so that a deformation u(r) is applied, in such a way that the coordinates are transformed as If a strain is defined in the usual way the stress of the deformed system is found as with E being the elastic energy density. It is important to remark that, if rotations are excluded from the dynamics, the strain tensor is symmetric and the following property holds Let us assume now that the Hamiltonian of the system is given by where p α = −i ∇ α , r α and m α are the momentum, position and mass of the α particle, respectively. The term V (r α ) stands for the periodic potential of the lattice. The term U (r αβ ), whith r αβ = |r α − r β |, is a two-body interaction between electrons. Once the transformation is applied, each of these terms will contribute in a different way to the energy of the system. If a small deformation is assumed, the energy of the deformed system will be a function of the deformation u and the strain ∇u. The difference between the energy of the deformed and undeformed systems is called the elastic energy, and it will allow us to define the elastic constants of the ensamble and to deduce the equation of motion.

B. Kinetic Term
The kinetic energy of the deformed system, using the transformation ∂ i = ∂ i − ∂ j u i ∂ j , is given by The first term of the right hand side part of the above equation is the kinetic energy of the undeformed system, and the second term will be where we have used the equipartition energy theorem, Finally, the last term is so that the contribution to the elastic energy of the kinetic part is given by If rotations are ignored, the energy density can be expressed in terms of the strain as The volume Ω is also a function of the strain, and is with Ω 0 being the volume of the undeformed system. The energy density is therefore given by from which we can obtain the stress of the system using equation (25) where and If we consider only the kinetic part of the energy we obtain the free electron gas, and it is shown that the unperturbed system is at a pressure P defined as which, according to equation (38), is identical to that of the free electron gas derived from Fermi-Dirac statistics. However, equation (39) shows that the electron gas is not a fluid, but a solid with Lamé coefficients λ S = µ S = P, which actually gives the same bulk modulus of the free electron gas although the nature of the gas is not that of a fluid but of a solid, due to the presence of the shear modulus µ S . This result, already obtained in references [35,36] using kinetic arguments, suggest that the hydrodynamic model, which ignores this shear modulus, should be revisited. The electron gas moves actually as an isotropic solid, and its equation of motion should be that of elastodynamics. In the following two subsections the interaction terms will be added which will consist in a body force due to the single particle potential and a generalized stiffness tensor due to the two-body interaction.

C. One-Particle Term
Let us consider the contribution of the periodic potential to the elastic energy. Expanding the potential energy up to second order in the deformation we obtain The first term of the above expression corresponds to the unperturbed system, the second and the third ones contribute therefore to the elastic energy. If we take the expected value of the energy we obtain, where we have integrated by parts for the first term but not in the second one. The energy density is therefore where we have considered as well the variation of the volume Ω with the strain. We see then that the effect of the single particle potential is to add a quantity V /Ω 0 to the equilibrium pressure P and to the Lamé parameter λ S , while the shear modulus is not affected. Additionally, a term proportional to the square of the deformation appears, which will be the responsible of a body force density f = −∇ u E, being This term is clearly the responsible of the local dielectric function, and we would recover the classical oscillator model if we neglect the contribution of the non-local strain contribution.

D. Two-Body Interaction
The two-body interaction is traditionally introduced in the theory of phonons to derive the acoustic equation of motion, however in our case the two-body interaction takes place between electrons, since we are interested in the optical regime where the nucleus will remain at rest. If we assume that the deformation is small we have that Then, using r αβ = |r α − r β |, the two body potential is and The two-body interaction contributes therefore with a fourth-rank tensor to the stiffness of the system. Since this contribution has to be averaged through the unit cell and electrons are assumed to be Bloch wave functions, we expect this tensor to have the symmetries of the lattice.

IV. ELASTODYNAMIC FORMULATION
In the previous section we have shown that the elastic energy density of a deformed ensemble of electrons is given by where the intrinsic stress of the system is and the stiffness and body force tensors are given by and respectively.
We have assumed however a timeindependent deformation, for which the above energy is uniquely the potential energy of the system. To derive the classical equation of motion fo the deformation we have to assume that a kinetic energy term due to the time-variation of u will appear, so that the Lagrangian density of the ensemble can be expressed as where ρ M is the mass density of the electron solid. Assuming a constant intrinsic stress, results in the following equation of motion [37] Since the induced polarization is −ρ e u and the external force is due to the electric field, then F e i = ρ e E i , where ρ e is the charge density due to electrons, the above equation can be expressed as the system of equations where we have used the definition of the plasma frequency ε 0 ω 2 P = ρ 2 e /ρ M . Finally, dissipation can be added phenomenologically to this model. The problem of dissipation is complex and an accurate description is beyond the objective of the present work. It will be however considered in a phenomenological way. The local mechanism for dissipation is usually introduced by means of Ohms law, in which a finite conductivity appears so that an induced current proportional to the electric field is excited.
However, dissipation in elastodynamics occurs in a different way, since for both fluids and solids it is due to forces proportional to the time derivatives of the strain. In our model, due to the presence of energy terms proportional to both the square of the strain and the deformation, we will assume that the two types of dissipation could appear, i.e., due to time derivatives of the strain and the deformation, what means that a local dissipative force F γ = −γ∂ t P will appear in equation (56) and a viscous stress σ η = η : ∂ t ∇P will be added in equation (57).
In summary, the equation of motion of the electromagnetic field inside matter when no sources are present can now be described by means of the following set of equations, where we have used in equation (61) the continuity equation We can combine equations (58) and (59) in the usual way to have a second order partial differential equation relating P and E. Similarly, we can combine equation (63) and (62) and we will have another second order partial differential equation relating P and E. Combining these two we will have a fourth order partial differential equation, which in turn means that the eigenvalue problem for the fields will be a 3 × 3 matrix equation with a fourth power in the wavenumber, so that in principle we will have three polarizations times the four solutions for the wavenumbers, i.e., a total of twelve modes, but the restrictions due to the equations (60) and (61) reduces by two the number of modes, then the most general solution of equations (58) to (63) in a homogeneous material will consist in ten propagating modes, but with double degeneracy due to reciprocity. This means that, at an interface between two different materials, we will need to match five modes at each side, so that we will need ten boundary conditions. Equations (58) to (63) will allow us to derive these ten boundary conditions in a natural way, as will be discussed in the following section.
However, in some practical situations one of the materials can be a local material, as is the case of vacuum, local dielectrics or metals. In these situations, we will have less modes than boundary conditions, so that we will need to find the boundary conditions that are no longer satisfied. This is the opposite to the traditional approach in non-local electrodynamics, in which a great effort has been done trying to find additional boundary conditions. Clearly this approach is more efficient, since going from the most general situation to particular cases is always easier than beginning with a particular case and trying to find the most general one, as will be seen in next section.

A. Boundary Conditions
The set of equations (58) to (63) define the evolution of the electromagnetic and polarization fields. Boundary conditions arise in a natural way in this description, since at an interface (flat or infinitesimally flat) the parallel wavenumber is a conserved quantity, and we can make in these equations the substitution ∇ → n∂ n +ik t . Since no singularities are allowed in any of the fields (no surface currents or fields), we impose the continuity of all the fields in front of the normal derivative ∂ n . It is easy to deduce then that the following conditions have to be satisfied at the boundary, The same conclusion could be reached with the traditional pill-box and circulation arguments. The first two equations are the well-known continuity equations of electrodynamics, the only difference is that in magnetic materials the microscopic field B t has to be replaced by H t . The last two equations are identical to elastodynamics, where the continuity of the normal components of the stress tensor and the displacement vector are required. The electrodynamic boundary conditions provide four equations (there are two components of each transverse field), while the elastodynamic ones provide us of six (the normal component of a second rank tensor is a threecomponent vector, and the polarization vector has three components), therefore the above equations perfectly defines the general boundary value problem at the interface between two materials, where ten equations were needed to match the five modes excited at each material, as explained before.
It is interesting to point out that within the elastodynamic description the required boundary conditions are the continuity of the transverse components of the electric and magnetic fields, however, from equation (3) the continuity of the normal component of the magnetic field also holds, as usual, but also the continuity of the normal component of the electric field, since in equation (4) the continuity of the normal component of the polarization field also implies the continuity of the normal component of the electric field. The absence of discontinuities in the fields suggests that the elastodynamic description is more suitable for the study of nanostructured materials, where the continuous nature of the wavefunctions of the different polarization carriers has also to be imposed.

B. Vacuum and Local Materials
The boundary conditions derived in the previous subsection emerge in a natural way at the interface between two nonlocal materials. However, in many practical situations, one of the materials might be vacuum or a local material, in the sense that the effects of the stiffness tensor might be neglected. It will be also interesting to determine the conditions for the consideration of non-local effects, that is, the conditions under which the material cannot be considered a local material and this more complex theory has to be applied. We expect this limit happens for low frequencies, but these conditions will be derived later on. In this subsection we will just consider what happens at an interface between a non-local material (solid) and a local one, which can be either vacuum or a local dielectric with finite conductivity, since, in terms of the number of solutions, all these materials are identical.
When one of the interfaces is vacuum or a local material the number of modes to match reduces from ten to seven, since now we have the five modes of the solid material plus the two polarizations allowed in local electrodynamics. Therefore, we will need as well seven boundary conditions for this interface.
Let us consider the local material first. We can assume that this situation will happen when no free or nearly free electrons are allowed (insulator), so that the stiffness component due to the kinetic term cancels, and also when the two-body interaction between electrons is negligible, due to the fact that the restoring force constant κ is so high that electrons remain bounded around the nucleus. The material is then polarizable only locally. At the interface, the electrons traveling through the solid find a very tinny potential barrier, so that the electronic wavefunction will be continuous across the interface and, consequently, the polarization vector. However, it is clear that the free electrons from the solid will be able to apply a force at the boundary, but not the electrons at the local side, since the stress tensor is zero in the local material. In the case of local metals, these can be seen as those materials where the local dissipation term γ is high, so that the mean free path of the electrons is small. The effect is that the behaviour of the material is local, since electrons are similarly "trapped" in a finite region, so that we will have a similar situation as described before.
Thus, the stress tensor can be discontinuous at this interface, and it is this boundary condition which is released here. We have then the following seven boundary conditions at the local-solid interface: It is tempting now to define vacuum as a "material" in which κ is zero and also the stiffness tensor C. However, this would be a wrong picture of the vacuum-solid interface. Vacuum is indeed a region non-accesible for the electrons of the solid, in the sense that these are bounded on its volume, we can consider therefore that it is the same situation as before but with the potential barrier found by the electrons in the solid being infinite, then it is the case when κ → ∞ so that we have to impose the cancelation of the wavefunction at its boundary and, equivalently, the cancelation of the polarization vector. This is indeed the same boundary condition as (66c), since in vacuum the polarization vector is zero and its continuity implies the cancelation at the boundary. This is equivalent to the mechanically rigid body, in which we impose that the displacement of the surface be zero, so that P = 0, while we are able to apply a force on its surface, meaning that σ will be different than zero there. Then, boundary conditions will be The above boundary conditions are the ones used by Pekar [26], using the argument that the excitonic wavefunctions should be zero at a vacuum interface. In section VI we will numerically show the above boundary conditions as a limiting case of κ → ∞, and we will show numerically the transition from a non-local material to a local one. The properties of plasmonic materials, when described by means of the hydrodynamical model [27], implies only one additional mode, the longitudinal one, thus it is not possible to cancel the three polarization components, and only the normal one is imposed . However, the justification of the third "additional boundary condition" is typically done by suggesting that since no charge is leaving the surface we need that the normal component of the current and, therefore, of the polarization vector, has to be zero there. However this is a tricky argument, since in the local description of metals with finite conductivity the normal component of the electric current is different than zero, and it does not mean that charges are leaving the surface, it means that we have charges on the surface. This is what happens indeed during the non-stationary regime of the charge of a capacitor.
In this unified picture of spatial dispersion, these boundary conditions arise in a more natural and general way, and we can assume that the normal component of the polarization is continuous in general and it cancells at the solid-bacuum interface, because vacuum is a rigid body and then its surface cannot be displaced.

C. Final remarks about generalized boundary conditions
It is the current view within the community working on non-local effects that boundary conditions cannot be uniquely established, mostly because they depend on the microscopic properties of each material and interface. Although this is certainly true for local materials, we propose to deal with the "microscopic" dependence of boundary conditions by including a few macroscopic quantities which help to establish the type of boundary in each interface. This approach essentially includes the strain tensor σ, and deduces a local relationship with the induced strain ∇P . With this simple hypothesis several response models are found. For instance, the different boundary conditions discussed in [28] or [31] can be derived by just changing the values of the stiffness tensor C, local resonance ω R , plasma resonance ω P and inertia ρ M on one side and the other of the interface, however here boundary conditions are derived within the framework of the field equations, which is a more rigorous procedure. The statement "boundary conditions depend on the microstructure of the material" acquires a precise meaning in this context, since, obviously, the different values and symmetries of these parameters depend on the microstructure of the material. Limiting situations, like perfect conductors, are perfectly defined as particular situations where the boundary conditions are different.
The present formulation renounces to the design of the constitutive parameters, it can be left behind to other domains of physics, and within the framework of this theory we can just study the physical properties of these materials assuming at least the order of magnitude of these parameters, but also imagine new phenomena and devices and guess which values of these constants would be required. Obviously this is not the most general theory possible but, as we have found in the literature, it is a big step in the domain of nanophotonics since it unifies in just one model several materials widely studied.
The elastodynamic model is "closed", in the sense that we do not need additional physical considerations to properly define boundary conditions, while it is true that boundary conditions when spatial dispersion is present requires of the knowledge of the microstructure of matter, it is also true that phenomenological theories have been derived in all domains of physics, where the different limiting values of the constitutive parameters can define different materials and responses, and this description seems a good step forward towards the general understanding of spatial dispersion at the nanoscale.

V. SOLUTION IN ISOTROPIC MEDIA
Although our approach is general and can be applied to any anisotropic material, we will illustrate the application of this method with the isotropic case.
Let us assume that the fields have time-harmonic dependence of the form exp(−iωt), in this case the constitutive equation (63) defines a complex stiffness tensor which in isotropic media has only two independent components, then we have where α and β are complex parameters of the form α = α R − iωα I . The stress tensor is in this case given by and the equation of motion for P becomes (α + β)∇∇ · P + β∇ 2 P − ΓP = −ρ M ε 0 ω 2 P E, (70) where we have defined the frequency-dependent Γ constant as Γ = ρ M (ω 2 R − ω 2 ) − iγω, with κ = ρ M ω 2 R . The above equation is equivalent to which in Fourier space gives Similarly, Fourier transforming equations (59) and (58) we get from which we can solve for the longitudinal wavenumber and the transverse one whose solutions are which shows that for real β and Γ, i.e., with no dissipation, the squared wavenumber k 2 T is always real, but can be negative or positive, therefore we will have either a propagating wave or a non-dissipative evanescent one.
In summary, we have a longitudinal wave with wave number k L and two transverse modes with wavenumbers k T given by the two solutions of (79). Taking into account that each transverse mode is decomposed in two polarizations, we have a total of five modes propagating through the bulk material in each direction, as discussed before.

A. The quasi-local limit
It is interesting to see that, in the absence of dissipation, if ω 2 R + ω 2 P − ω 2 < 0 the longitudinal mode is propagative, and also the two transverse modes k T , since the product of the two solutions of equation (79) is However, if ω 2 R + ω 2 P − ω 2 > 0, only one of the transverse modes is propagative, and the longitudinal mode is also evanescent. We have therefore only two propagative modes and three evanescent ones. This is the situation that will be analyzed in this work, since in terms of propagative modes the number of solutions is the same as in local materials, however the role of the evanescent modes is to keep the continuity of the fields at the interface. Also, this situation will allow us to recover the local limit, as will be seen later, therefore we call this limit the "quasi-local" limit.
Then, in the limit k 2 0 β → 0, that is, the low frequency and low spatial dispersion limit, we can expand the solution for k 2 T as which gives two physically different solutions, the propagating one (taking the plus sign in the ± term), and hereafter labeled a1, and the evanescent one, labeled as a2, where we have defined the static electric susceptibility χ E as and the characteristic length 0 as (85) Therefore, in this limit, the propagative solution is identical to that in a dielectric material with dielectric constant 1 + χ E , and quickly evanescent modes will be excited up to a distance 0 from the surface. The case of conductors is slightly different. A Drude metal is found in the limit of ω R → 0, with γ = 0, in the quasi-local limit we would have, assuming that ω 2 P > ω 2 , while in this case so that even below the plasma frequency this model predicts the existence of a propagating mode at the speed of the electronic shear wave c 2 S = β/ρ M . This velocity can be estimated in the following way: since β if of the order of the elastic constants of a solid, while ρ M is about two thousand times smaller than its mass density (the ratio between the proton and electron mass) we can estimate the speed of this wave of being about forty times faster than the speed of sound in solids. This additional mode will have a ver short wavelength, so that it will be also very quickly attenuated.
Indeed, if we consider the dissipation factor γ, we recover from χ E the conductivity predicted by Drude model with σ 0 = ρ 2 e τ /ρ M and τ = ρ M /γ, as expected, while for the second mode we have which gives so that the length 0 is now approximately given by the distance covered by the electrons traveling at the speed c S , i.e., the shear speed of the ensemble. It will be clearly a short distance, and the local limit will be recovered in this case taking the limit of this distance to zero, with the boundary conditions derived before, as will be demonstrated later on.
B. The characteristic length 0 The term "non-local" has been used through this work paying attention specially to the consequences in terms of wave propagation and the number of modes to match. Although it is this analysis the relevant one in terms of the solution of boundary value problems, it is also interesting the analysis of the constitutive equations in real space, since it will help us to better understand the physical interpretation of the stiffness tensor C. We will limit our discussion to the isotropic case and neglecting the longitudinal response, but a deeper analysis can be found in [38].
For simplicity we will assume that α = 0, so that equations (73) and (74) defines a frequency and a wavenumber dependent electrical susceptibility which in real space relates the polarization and the electric field as where the non-local susceptibility is or, using the radial symmetry of the Fourier transform, Therefore, the polarization at a given point r is a linear combination of the electric field through all the space weighted by the non-local susceptibility χ E (r, ω). For both Γ and β reals and positive (the quasi-dielectric limit), we have (see [39], section 273, equation 3) where we have defined the characteristic length 0 as The meaning of the parameter 0 is therefore clear in this limit: a region of radius 0 has to be taken into account to compute the polarization P , since outside this region the susceptibility χ N L is negligible.
Once defined the characteristic length 0 , the effects of spatial dispersion will be more or less important depending on the wavelength of the field. As will be seen later, when λ >> 0 the influence region is averaged and the behaviour of the fields is similar to that of a local material, while when λ ≈ 0 the effects of spatial dispersion are more relevant and cannot be ignored.
The above discussion is valid only in the quasi-local limit, that is, when all the additional modes are evanescent and we still have only one propagative solution. When either Γ and β are complex, the residue theorem has to be applied to evaluate equation (94), and we find that the exponential decay rate depends on the absolute value of the real part of Γ/β. Then, let us assume for instance that β is real, and Γ = ρ M (ω 2 R − ω 2 ) − iωγ, we will have .
The above result shows that, as long as ω 2 R −ω 2 > 0 the characteristic length is of the order of 0 , as before, however when ω 2 R − ω 2 < 0 the imaginary and real parts are exchanged, and the characteristic length becomes proportional to 1/δ, which can be very high, depending on the value of this parameter. The same effect applies as well when β has both real and imaginary parts, but not when it is purely imaginary. A deeper analysis of this complex situation is beyond the objective of the present work, which is only to understand how the effects of spatial dispersion affects boundary conditions, therefore in the numerical calculations we will consider only the quasidielectric limit, in which the characteristic length is small and we can compare results with the local limit.

VI. NUMERICAL EXAMPLES: PLANAR INTERFACE
In this section we will consider the planar interface between two solids, labeled a and b, in order to analyze the consistency of the boundary conditions derived previously, as well as the conditions in which we recover the local limit. We will assume that material a is at the left hand side of the interface and material b is at the right hand side, as shown in figure 1. Let us assume that we are in the quasi-dielectric limit, in which one transverse wavenumber is propagative and the other transverse wavenumber is evanescent, as well as the longitudinal one. It is convenient now to derive the relationship between the different fields involved in boundary conditions and the electric field. Thus, after performing the substitution ∇ → ik t + ∂ n , from equation (58) we have that the transverse magnetic field is given by while the polarization vector is given by The ij element of the stress tensor σ in an isotropic material is given in by equation (69), so that the normal and transverse components, which are to be continuous, are σ nn = (α + 2β)∂ n P n + iαk t · P (100) σ nt = β∂ n P t + iβk t P n .
In the next subsection we will analyze the behaviour of these fields at normal incidence, and later on we will derive the general reflection coefficients at the interface between vacuum and a solid material considering oblique incidence.

A. Normal incidence
Let us asume that a propagating transverse wave is excited in material a and propagates along the x axis, so that n =x. Without loss of generality, we will assume that the interface is placed at x = 0 and that the field is polarized along the z direction. If the excited mode has wavevector k = k a1x , after reflection two additional modes are excited in material a, corresponding to the two solutions of the dispersion equation (79). Since there are no normal components of the fields there is no excitation of the longitudinal mode, whose discussion is left for next subsection. Then we have, for x < 0, The transverse component of the magnetic field is and the polarization vector is parallel to the electric field, therefore The normal component of the stress tensor σ is zero, and the transverse one is For x > 0 we have two modes excited as well, thus In summary, we will need to determine the value of four coefficients, and we have indeed four equations for this geometry, the continuity of E z , B y , P z and σ zx , then the solution for the coefficients can be found after inversion of the system The above equations show that the interface problem is well defined once we have defined the two transverse wavenumbers (at normal incidence there is no excitation of the longitudinal mode) and the β coefficient. We will see now how the above equations, derived using equations (65), can degenerate in equations (66) as the limit of low characteristic length 0 . Figure 2 shows the polarization P (x) (upper panel) and the stress tensor σ xz as a function of x along an interface between two solid materials. Material a is chosen so that β a = 3, k 2 a1 = (1 + χ Ea )k 2 0 and k 2 a2 = −1/ 2 a , with χ Ea = 2 and a = λ/10. Similarly, material b is selected as β b = 1, k 2 b1 = (1 + χ Eb )k 2 0 and k 2 b2 = −1/ 2 b , with χ Eb = 3 and b = λ/10, λ/50 and λ/100. We see that both solids are non-local but we have analyzed the evolution of material b towards a local material, showing therefore that in this transition the polarization vector remains continuous while the stress tensor presents a discontinuity, as required by boundary conditions (66). . We see how in the transition b → 0 material b becomes a local material and the stress tensor presents a discontinuity, as predicted by the theory (see text for further details). Figure 3 shows the same situation as before but now we set χ Eb equal to 0, so that in the limit of b → 0 material a converges towards vacuum. We see that, as expected for vacuum, both the polarization and the stress fields cancels, however the polarization remains continuous (zero) at the interface while the stress presents a discontinuity.
Once more this agrees with boundary conditions (67).  Figure 4 shows the same situation but now we set χ Eb = −100, which corresponds to a metalic material. The meaning of the characteristic length b is similar, as shown in equation (90), but now we add the contribution of the short wave to k b2 , so that we have k b2 = k S + i/ b , with k S = 10k 0 and b having the same range as before.
Finally, figure 5 shows the behaviour of the polarization vector at a vacuum-solid interface, when the solid has a dielectric (upper panel) or a conductor (lower panel) character, with parameters β b = 2, χ Eb = ±3, with the +(-) sign corresponding to the dielectric (conductor). For the conductor case we also set k S = 10k 0 and results are shown for b = λ/10, λ/25, λ/50 and λ/100. Boundary conditions imply the cancelation of the polarization at the interface, but we see how, as the characteristic length becomes smaller, we recover the local behaviour for both types of materials, and the polarization presents a step discontinuity for the dielectric while it is concentrated at the surface for the conductor. It has to be pointed out that, within the hydrodynamic description, the behaviour of the conductor could not be described in this way, since at normal incidence it is identical to a local metal, so that this description is clearly more correct if microscopic arguments are to be included.
The presented approach is clearly consistent in the limiting situations, in which we recover the "traditional" boundary conditions as a progressive limiting situation, something that is not given in other approaches, which  ). We see how, although the polarization cancels at the surface, as we reduce the characteristic length the polarization tends towards a step discontinuity in the case of the non-local dielectric and to concentrate all the current at the surface for the conductor, as predicted by local eletrdodynamics but also in agreement with the present theory.
just define a set of boundary conditions but do not specify when these are valid or not.

B. Oblique incidence
In the previous subsection we considered the different interfaces at normal incidence, and this analysis allowed us to study the behaviour of the fields when material's discontinuities appear. In most of the real problems the fields do not incide normally at a surface, as for example in cylindrical or spherical objects, therefore the analysis of the behaviour of the fields when a given angle θ 0 is formed with the normal to the surface will help us to understand the influence of spatial dispersion in more general problems. As will be seen below, the main difference is that the longitudinal mode is excited when a normal component of the electric field appears at the interface. This longitudinal mode will be required to satisfy the continuity of the fields, and it will play a similar role as the evanescent mode studied before.
For an isotropic material the possible solutions of the electromagnetic field can be decomposed in one longitudinal (L) and two transverse (T ) modes, and each of the transverse modes has two components, which we define as the S 1 , S 2 and P 1 , P 2 polarizations.
Since the component of the wavevector parallel to the interface is a conserved quantity and, consequently, identical to all the polarizations, we define the wavevector as where the index i indicates in which material the wave propagates (i = a, b in our case) and σ = L, S 1 , S 2 , P 1 , P 2 . The vector k t is the component of the wavevector parallel to the surface, therefore the component of the wavevector normal to the surface is given by with k iσ = |k ± iσ | being k L for σ = L and k T for the different transverse modes, i.e., the solutions of equations (77) and (79).
The unit vectors parallel to these polarizations are given by and the following relationships are found It has to be pointed out that the u S and u P are orthogonal for the same root of equation (79), but not for different roots. Also, u L is not in general orthogonal to these vectors.
We assume that material a is vacuum, and that the incident electric field is a plane wave with a general polarization state, thus E 0 = σ=S,P A σ e iqaσn e ikt·r u + σ (119) which excites a reflected field given by Material b is a non-local solid, therefore we will have two solutions for each of the S and P polarizations plus the longitudinal mode L, therefore the transmitted electric field will be The inputs of the system are the coefficients of the incident electric field A S and A P , and we have to obtain the two reflected amplitudes B S and B P and the five transmitted amplitudes C S1 , C S2 , C P 1 , C P 2 and C L . We need therefore, as discussed in section IV A, seven equations, which correspond to the seven boundary conditions (67a)-(67b).
It is easy to see that the S polarization uncouples from the P and L, which in turn means that the incident A S field excites only the B S reflected field and the C S1 and C S2 transmitted fields. The cancelation of the polarization vector parallel to u S implies while the continuity of the transverse components of the electric and magnetic fields gives A S − B S = q b1 q a0 C S1 + q b2 q a0 C S2 from which we can obtain the reflection coefficient of the S mode R S = B S /A S as We can proceed similarly to obtain the reflection coefficient of the P polarization R P = B P /A P , using then the continuity of the transverse electric and magnetic field is while the cancelation of the polarization vector gives From the above two equations we can obtain the relationship of the C P1 with the C P2 and C L coefficients, giving q bL q b1 + |k t | 2 q bL q b2 + |k t | 2 so that finally we obtain the R P coefficient as 1 − R P 1 + R P = k 2 0 q a0 q b1 /k 2 b1 + q b2 /k 2 b2 χ P + 1/k L χ L 1 + χ P .
At the interface between vacuum and a local dielectric material the normal component of the electric field is not continuous, and a bounded surface charge exists. This is due to the discontinuity of the polarizability of the material which induces a discontinuity in the polarization vector and, from the continuity equation, ρ = −∇ · P = −ik · P − ∂ n P n ∼ δ(n) this discontinuity implies that it is actually the normal component of dielectric displacement D = (1 + χ)E the quantity continuous through the interface. In the elastodynamic description however, due to the continuity of the polarization vector P , there are no bounded surface charges, consequently the normal component of the electric field is also continuous at an interface.

VII. SUMMARY
In summary, a self-consistent theory for the interaction of the classical electromagnetic field and matter has been presented. The theory, developed within the framework of elastodynamics, is based on a polarization vector whose equation of motion is identical to that of the elastodynamic field, and a set of boundary conditions arise in a natural way whose number is consistent with the number of electrodynamic modes. Elementary considerations about the limiting value of the parameters of the model allow us to define vacuum, local dielectrics, real conductors and hydrodynamic plasmas, recovering in each situation the boundary conditions employed in the literature. This description however includes the possibility of more advanced interfaces, composites and symmetries, since a full anisotropic description has been considered.
Therefore, the theory presented here can be a starting point for any model of matter at the nanoscale, where the effects of spatial dispersion and continuity of the mi-croscopic fields can be more relevant, and more refined models of matter are required.
This unified theory for the description of nanomaterials will doubtless provide a new and generalized description of nanophotonic structures which will be fundamental for either the characterization of nanomaterials and the accurate design of new nanophotonic devices.