1 Introduction

In macroscopic electrodynamics at near-IR or optical frequencies, linear material properties are commonly described by frequency-dependent permittivities \(\epsilon (\omega )\). Their experimental determination happens, e.g., by means of ellipsometry [1] and the so-obtained data can directly be fed into computational approaches that solve Maxwell’s equations in the frequency domain. Their incorporation into time-domain Maxwell solvers requires more effort because, in the time domain, the constitutive relation \(\vec {D}(\omega ) = \epsilon _0 \epsilon (\omega ) \vec {E}(\omega )\) translates to a convolution integral, thus turning Maxwell’s equations into a set of numerically demanding integro-differential equations. Here, \(\epsilon _0\) denotes the vacuum permittivity.

Within the framework of auxiliary differential equations (ADEs), this convolution integral is circumvented and replaced by a set of differential equations for auxiliary polarizations or auxiliary currents. This is accomplished by writing the permittivity as a sum of rational functions in \(\omega\) that respect appropriate Kramers–Kronig relations in order to ensure causality [2] (for concrete examples of ADEs, see Sect. 3.1 below). For instance, metals are often treated via the so-called Drude–Lorentz model

$$\begin{aligned} \epsilon (\omega ) = 1 - \frac{\omega _p^2}{\omega (\omega +i\gamma _{\textrm{D}})} + \sum _{m=1}^N \frac{f_m \omega _p^2}{\omega _m^2-\omega ^2 - i \omega \gamma _m}. \end{aligned}$$
(1)

Here, the 1st and 2nd term on the r.h.s. represent, respectively, the vacuum contribution and the (collective) response of the metal’s free electrons, the so-called Drude model with plasma frequency \(\omega _p\) and phenomenological damping constant \(\gamma _{\textrm{D}}\).

In dielectrics, the 3rd term on the r.h.s. of Eq. (1) models electrons that are elastically bound to their parent nuclei so that the \(f_m\) and \(\omega _m\) \((m=1,\dots , N)\) may be interpreted, respectively, as oscillator strengths and frequencies of corresponding chemical bonds or may be associated with certain transitions in the electronic band structure of semiconductors. This particular form of the linear response has been suggested on empirical grounds by Sellmeier in 1871 [3]. In addition, a set of phenomenological damping constants \(\gamma _m\) \((m=1,\dots , N)\) is introduced. As a result, each Lorentz pole features three independently adjustable parameters that are obtained from fits to experimental permittivity data, where the oscillator strengths obey a sum rule. Often, additional frequency dependencies are associated to the damping constants [4].

In metals, Lorentz poles are used to account for interband transitions (IBTs) so that \(f_m\), \(\omega _m\), and \(\gamma _m\) (\(m=1,\dots ,N\)) and the number N of Lorentz poles represent effective parameters without straightforward physical interpretation. Nonetheless, sophisticated and rather cumbersome methods [5] have been developed to fit the parameters of the Drude–Lorentz model, Eq. (1), to experimental data such as those provided by the classic measurements of Johnson & Christy for noble metals [6] or the data reported in the handbook of optical constants of solids edited by Palik [7].

In the present work, we adopt a more microscopic viewpoint on \(\epsilon (\omega )\) and demonstrate (i) how this leads to a considerable reduction in the number of fit parameters and (ii) how the remaining parameters may be interpreted physically. Further, our approach provides an alternative perspective on Eq. (1) and offers a systematic way of improving the accuracy of time-domain computations without introducing additional fit parameters. Our work is organized as follows: In Sect. 2, we review some of the available literature, develop our approach first in the frequency domain and apply it to the permittivity of gold, the perhaps most frequently used plasmonic material. We proceed in Sect. 3 to the adaptation of our approach for time-domain simulations and provide numerical assessments of its validity. In Sect. 4, we summarize our findings.

2 Theory of the frequency-dependent permittivity of metals

One of the earliest discussions on the permittivity of metals dates back to two works of Darwin [8, 9]. There, on a classical level, it was explained that the bound and free electron response are additive and that there is no local field correction in a plasma. On the quantum mechanical level, Bohm and Pines [10,11,12,13] have elaborated the physical properties of metals. Regarding the metals’ response to electromagnetic fields [12], they have demonstrated that the complicated multi-particle interactions in a metal can be classified into two types which require rather distinct treatments. More precisely, there exist collective excitations of the electron system that are dominant at long wavelengths (low energies) and they can be treated via Fermi liquid theory (Landau–Silin theory [14]) where the periodic (i.e., atomic) nature of the ionic background is essentially irrelevant. Through a corresponding canonical transformation of the original Hamiltonian, the remaining short-range interactions are then mapped onto a collection of independent particles that move in an effective periodic potential. Consequently, this separation of spatial scales leads to a description of the electron system as a mixture of a Fermi liquid (FL) and a set of non-interacting effective particles (EPs).

It is to this mixture that the aforementioned arguments of Darwin may be applied so that the linear permittivity of metals may be decomposed into

$$\begin{aligned} \epsilon _{\textrm{Metal}}&= 1 + \chi _{\textrm{FL}} (\omega ) + \chi _{\textrm{EP}} (\omega ), \end{aligned}$$
(2)

where, \(\chi _{\textrm{FL}} (\omega )\) and \(\chi _{\textrm{EP}} (\omega )\) denote the linear susceptibilities of the Fermi liquid and the effective particles, respectively.

2.1 Drude model and extensions

The Fermi liquid contributions can be treated in numerous ways with varying degrees of sophistication. The aforementioned Drude model

$$\begin{aligned} 1 + \chi _{\textrm{FL}} (\omega ) \approx \epsilon _{\textrm{D}}(\omega ) = 1 + \chi _{\textrm{D}} (\omega ) = 1 - \frac{\omega _p^2}{\omega (\omega +i\gamma _{\textrm{D}})} \end{aligned}$$
(3)

represents the simplest approach that works rather well for simple metals and can even be derived from purely classical arguments [8] where the plasma frequency \(\omega _p^2 = n e^2/(\epsilon _0 m)\) is expressed in terms of the density n, the charge e, and mass m of conduction electrons. In an actual Fermi liquid, n and m need to be replaced by density and effective mass of the quasi-particles. Since, however, \(\omega _p\) is usually determined by fitting \(\epsilon _{\textrm{D}} (\omega )\) to experimental data (see [15] for a recent example of fitting the Drude permittivity of gold), in practice, this difference does not matter too much.

For strongly correlated systems, such as the ferromagnetic metals nickel, cobalt or iron, the Drude model must be extended to better account for the electronic correlations [16].

Further, if other relevant length scales such as the size of a metallic particle or the distance between two metallic particles become comparable to the electrons’ mean free path, nonlocal effects, such as Landau damping (which already occurs in classical theories) and the Fermi pressure (a purely quantum statistical effect), as well as the detailed behavior of the liquid near surfaces have to be taken into account [17]. Since the Fermi liquid contributions are not at the focus of our work, we would like to refer to the recent literature [17,18,19,20,21,22,23] and references therein for details.

2.2 Interband transitions

The effective particle susceptibility can be obtained from standard linear response theory. To do so, it is important to note that the effective particle moves in a periodic potential which leads to the formation of an effective energy band structure that is related to the actual electronic band structure in the following way. First, for a metal, the conduction band is partially filled and the electrons that occupy these states represent the free (conduction) electrons that are treated via Fermi liquid theory as described in the previous Sect. 2.1. Consequently, as the aforementioned canonical transformations facilitate a separate treatment of the Fermi liquid (i.e., the intraband contributions), the effective particle corresponds to transitions from bands that are completely filled to bands that are at most only partially filled. In other words, to the level of approximation discussed above, the effective particle contributions to the susceptibility are associated with interband transitions of an (effective) and potentially highly doped semiconductor. Wysin et al. [24] have adopted this viewpoint in order to determine the contribution of these interband transitions to the Faraday rotation in metallic nanostructures. The corresponding contribution of the hydrodynamic part has been developed in Ref. [16].

Following the above prescription, the contributions of the effective particles to the susceptibility can be determined via linear response theory under dipole coupling. This gives

$$\begin{aligned} \chi _{\textrm{EP}} (\omega ) = \frac{2 n}{\hbar } \sum _{l \vec {k}} \sum _{l^\prime \vec {k}^\prime } \frac{p_{l \vec {k}} \, \omega _{l l^\prime \vec {k} \vec {k}^\prime } \vert \langle l \vec {k} \vert \hat{\vec {d}} \vert l^\prime \vec {k}^\prime \rangle \vert ^2}{ \omega _{l l^\prime \vec {k} \vec {k}^\prime }^2 - \gamma _{l l^\prime \vec {k} \vec {k}^\prime }^2 - 2i \gamma _{l l^\prime \vec {k} \vec {k}^\prime } \omega - \omega ^2 }. \end{aligned}$$
(4)

Here, \(p_{l \vec {k}}\) and \(\omega _{l l^\prime \vec {k} \vec {k}^\prime } = (\epsilon _{l \vec {k}}-\epsilon _{l^\prime \vec {k}^\prime })/\hbar\) denote, respectively, the occupation of a Bloch state (labeled by band index l and wave vector \(\vec {k}\)) of the effective particle with energy \(\epsilon _{l \vec {k}}\) and the transition frequency between two Bloch states. In addition, we have introduced phenomenological damping constants \(\gamma _{l l^\prime \vec {k} \vec {k}^\prime }\) for the respective transitions. Finally, \(\langle l \vec {k} \vert \hat{\vec {d}} \vert l^\prime \vec {k}^\prime \rangle\) is the matrix element of the electric dipole operator \(\hat{\vec {d}}\) between two Bloch states. For the details of the calculations leading to (4), we refer to the work of Wysin et al. [24].

2.3 Parabolic two-band model for gold

We now apply the framework outlined in Eqs. (2), (3), and Sect. 2.2 to an actual metal—we chose the arguably most important metal for applications in plasmonics: gold. To do so, it is important to recall the results of corresponding electronic band structure computations. Specifically, electronic band structure theory reveals that gold features a fully occupied 5d band that, energetically, lies below the partially filled 6sp band [25] where the energy difference for direct optical transitions is between 1.8 and 2.45 eV, thus covering the frequency range from the red to the green/cyan range of the optical spectrum. As a result, the electrons of the partially filled 6sp band correspond to the Fermi liquid discussed above. Symmetry considerations stipulate that the dipole matrix elements between the 5d- and the 6-sp band are strongly suppressed. However, the 5d band exhibits a strong van Hove singularity near the X-point and near the L point of the Brillouin zone which compensates the small dipole matrix elements. Consequently, if we limit our considerations of the permittivity to the optical frequency range, we can restrict the effective particle susceptibility in Eq. (4) to a two-band model that only includes direct optical transitions, \(\vec {k} = \vec {k}^\prime\). As the transitions are appreciable near the X- and L points, and given, that the 5d band is much flatter in \(\Gamma\)-X and \(\Gamma\)-L direction than in other directions, we may adopt a parabolic approximation to both bands and limit the remaining summation in Eq. (4) to integrations along the relevant parts of the \(\Gamma\)-X and the \(\Gamma\)-L lines. As a matter of fact, recent ultrafast pump-probe experiments on the dynamics of hot electrons [26] suggest that transitions near the X- and transitions near the L point contribute in a comparable manner. Within our effective particle framework, we are thus left with two options. On the simplest level, we can combine the contributions of the \(\Gamma\)-X and the \(\Gamma\)-L line into one average effective particle contribution and this is what we will pursue in the remainder of the manuscript. However, at this point, we would like to note that a more refined treatment could be employed, where the two aforementioned contributions would be treated separately, i.e., we would have two terms of the form described in Eq. (5) (see also the discussion in Sect. 4).

Irrespective of this choice, due to the treatment of the EPs as effective two-band semiconductor, we have that the effective Fermi energy represents a fitting parameter that lies between the top of the lower and the bottom of the upper band. Assuming a step function for the occupation of the effective two-band model and further assuming that the corresponding dipole matrix elements and damping constants do not depend on wave vector, Eq. (4) evaluates to (for details of the calculation using the same notation, see Ref. [24])

$$\begin{aligned} \chi _{\textrm{EP}} (\omega )&= Q \int _0^{s_{\textrm{U}}} \textrm{d}s \, \frac{2s^2}{\left( \omega _g + s^2 \right) \left( \left( \omega _g + s^2 \right) ^2 - \left( \omega + i \gamma \right) ^2 \right) } \end{aligned}$$
(5)
$$\begin{aligned} Q&= \frac{2ne^2 \vert M \vert ^2 \hbar }{\epsilon _0 m_e^2} \cdot \frac{2V s_{\textrm{U}}^2}{\pi ^3} \cdot \left( \frac{\hbar }{2m^*} \right) ^{1/2} \end{aligned}$$
(6)

Here, we have introduced the scaled wave vector \(s = (\hbar /2m^*)k\) (with dimension 1/time\(^{1/2}\)), the upper limit \(s_{\textrm{U}}\) of the scaled wave vector integration, and have denoted the effective dipole matrix elements with M. The energies are measured from the top of the effective valence band, so that the dispersion of the parabolic bands are \(\epsilon _{l \vec {k}} = - \hbar ^2 k^2/2m_v\) and \(\epsilon _{l^\prime k^\prime } = \hbar \omega _g + \hbar ^2 k^2/2 m_c\), where \(m_v\) and \(m_c\) represent the curvatures of the effective valence and the effective conduction band around the Fermi surface, respectively. The reduced mass is defined via \(1/m^* = 1/m_v + 1/m_c\). The above integral can be evaluated exactly in terms of a sum of arc tangents. However, as we shall discuss in the next section, this does neither provide further physical insight nor does it help with formulating ADEs for time-domain simulations.

Combining Eqs. (2), (3), and (5), we now have a model for the linear permittivity of gold with a total of 6 fitting parameters: \(\omega _p\), \(\gamma _{\textrm{D}}\), Q, \(\omega _g\), \(\gamma\), and the integration limit \(s_{\textrm{U}}\). We expect that this model covers the low frequencies all the way to the blue end of the visible frequency range. In addition, we expect that the gap frequency \(\omega _g\) of the effective two-band model is close to the aforementioned energetic separation between the 5d- and the 6-sp band. Similarly, we expect that the effective matrix element M in Eq. (6) is close to the dipole moment of the associated transition, thereby allowing to interpret the quantity Q in Eq. (6) as an effective response strength that consists of the strength of the dipole transition weighted by the associated density of states (see also Ref. [25]).

Upon carrying out the above-sketched fitting procedure to the experimental data for gold by Johnson & Christy [6] via a standard Monte Carlo approach, we obtain the values \(\omega _p = 1.32 \cdot 10^{16}\) Hz, \(\gamma _{\textrm{D}} = 1.23 \cdot 10^{14}\) Hz, \(Q = 2.72 \cdot 10^{24}\) 1/Hz\(^{3/2}\), \(\omega _g = 3.63 \cdot 10^{15}\) Hz, \(\gamma = 2.41 \cdot 10^{14}\) Hz, and \(s_{\textrm{U}} = 2.66 \cdot 10^{8}\) Hz\(^{1/2}\). In Fig. 1, we display a comparison between the experimental data and our fit using Eqs. (2), (3), and (5).

Fig. 1
figure 1

Johnson & Christy data (dashed lines) for the complex refractive index (real part: blue lines, imaginary part: red lines) of gold in the range between 1600nm and 200nm as well as our best fit (solid lines) of the parabolic two-band model via Eqs. (2), (3), and (5). The fit has been done for the range of 1937nm to 187nm. This is approximately the range considered in Ref. [5]

Full size image

The above gap frequency corresponds to an energy of 2.39 eV and, thus, is consistent with expectation. The fit shows very good agreement with the experimental data for wavelengths up to and including the violet part of the visible spectrum. For even shorter wavelengths, deviations between the experimental data and the fit become more pronounced. Obviously, the fit could be improved on and further extended into the high frequency range by considering (energetically higher-lying) transitions between, e.g., the 4d-band and the 6-sp band thereby extending the two-band model to a multi-band model. Nonetheless, we would like to emphasize that the above results have been obtained with only 6 fitting parameters while corresponding fits of somewhat comparable quality using Eq. (1) require at least 5 or more Lorentz terms, totalling at least 15 or more fitting parameters. [5].

3 Time-domain modeling of interband transitions

As described above, in time-domain numerical schemes of Maxwell’s equations frequency-dependent permittivities are usually treated via ADEs. Unfortunately, while Eq. (3) is amenable to an ADE approach, Eq. (5) is not. And neither is its closed-form solution in terms of a sum of arc tangents. Therefore, in this section, we will develop an appropriate scheme for handling Eq. (5) via ADEs to any desired accuracy without introducing additional fitting parameters.

3.1 Discretization of the parabolic two-band model

The specific form of Eq. (5) suggests that we discretize the integral via a Gauss quadrature rule according to

$$\begin{aligned} \int _a^b \textrm{d}x \, f(x)&\approx \frac{b-a}{2} \, \sum _{i=1}^{N_{\textrm{G}}} w_i \, f(u_i). \end{aligned}$$
(7)

Here, \(u_i = (a+b)/2 + x_i(b-a)/2\) with \(i=1,\dots ,N_{\textrm{G}}\) denote the Gauss quadrature nodes which are expressed in terms of the roots \(x_i\) of the Legendre polynomial \(P_{N_{\textrm{G}}}(x)\) of order \(N_\textrm{G}\) within the interval \([-1,1]\) [27]. In addition, the weights \(w_i\) are given by \(w_i = -2/((1-x_i^2) (P_{N_\textrm{G}}^\prime (x_i)^2)\) so that Eq. (7) integrates exactly polynomials of degree \(2N_{\textrm{G}}-1\) [27].

Upon applying the above Gauss quadrature to Eq. (5), we obtain

$$\begin{aligned} \chi _{\textrm{EP}} (\omega )&\approx \sum _{m=1}^{N_{\textrm{G}}} \frac{a_m^2}{c_m^2 - (\omega +i\gamma )^2}, \end{aligned}$$
(8)
$$\begin{aligned} a_m&= \sqrt{\frac{Q \, s_{\textrm{U}} \, s_m^2 w_m}{c_m}}, \quad c_m = \omega _g + s_m^2, \quad s_m = \frac{s_{\textrm{U}}}{2} \left( x_m +1 \right) . \end{aligned}$$
(9)

With the above discretization, Eq. (8), of the single particle contribution within the parabolic two-band model, Eq. (5), we end up with a set of effective Lorentz poles which can be included into time-domain simulations of the Maxwell equations via ADEs. The number \(N_{\textrm{G}}\) of these Lorentz poles controls the level of accuracy with which we are approximating Eq. (5). However, contrary to the standard Drude–Lorentz approximations based on Eq. (1) or Lorentz pole-only approximations [5], our approximation scheme, Eqs. (2), (3), and (8), does not introduce any additional fit parameters when increasing the number of Lorentz poles (or, equivalently, when increasing the accuracy of the approximation).

Fig. 2
figure 2

Complex refractive index \(\sqrt{\epsilon (\omega )} = n(\omega ) + i k (\omega )\) of gold for different Gauss quadratures of the parabolic two-band model in the range between 200nm and 900nm. Left panel: Real part of the refractive index. Right panel: Imaginary part of the refractive index. The free interband parameters Q, \(\omega _g\), \(\gamma\), and \(s_{\textrm{U}}\) have been determined by fitting Eqs. (2), (3), and (5) to the experimental data of Johnson & Christy (dashed lines). When the number \(N_{\textrm{G}}\) of Lorentz poles (i.e., Gauss quadrature points) increases, both, real and imaginary part, converge to the full model

Full size image

In Fig. 2, we depict details of the accuracy of our effective Lorentz pole discretization for gold. A single effective Lorentz pole works very well for wavelengths larger than about 600 nm, an approximation based on three effective Lorentz poles works well down to about 550 nm. For smaller wavelengths, a considerably larger number of Lorentz poles are required for obtaining accurate results. In particular, by its very nature, the Lorentz pole approximation leads to unphysical peaks both, in the real and imaginary part of the complex refractive index which become less pronounced when the number of effective Lorentz poles is increased. Therefore, considerable care must be exercised when designing gold-based plasmonic elements for operation below 600 nm. Nonetheless, we would like to point out that independent of how many effective Lorentz poles we use, we still have only the four interband fit parameters discussed above. Using only these four fit parameters (along with their physical meaning, see Sect. 2.3 above), we have in hand a systematic way of improving the accuracy of our modeling of the dielectric properties of gold via ’simply’ increasing the number of effective Lorentz poles, thus allowing us to avoid unsubstantiated or unphysical designs and results.

3.2 Alternative discretization of the parabolic two-band model

Based on the above considerations, we consider a variant of the above approach. For a given \(N_{\textrm{G}}\), we may regard Eqs. (2), (3), and (8) as a replacement of Eq. (1) with fixed Drude parameters \(\omega _p\) and \(\gamma _{\textrm{D}}\) and a set of four fit parameters Q, \(\omega _g\), \(\gamma\), and \(s_\textrm{U}\). We then proceed to fit the resulting approximation to the experimental data. Note, that in this approach, for different values of \(N_{\textrm{G}}\), we obtain (slightly) different values for these four fit parameters. In Fig. 3 we depict the results of such a strategy for different numbers \(N_{\textrm{G}}\) of Lorentz poles. The corresponding values of the fit parameters are summarized in Table 1.

Fig. 3
figure 3

Complex refractive index \(\sqrt{\epsilon (\omega )} = n(\omega ) + i k (\omega )\) of gold for the alternative approach to Gauss quadrature approximations of the parabolic two-band model in the range between 200 nm and 900 nm. Left panel: Real part of the refractive index. Right panel: Imaginary part of the refractive index. Contrary to the results depicted in Fig. 2, in this approach the IBT-integral Eq. (5) is first discretized with a Gauss quadrature scheme with \(N_{\textrm{G}}\) points (i.e. Lorentz poles) according to Eq. (8). Then, for this choice of \(N_{\textrm{G}}\), the free parameters Q, \(\omega _g\), \(\gamma\), and \(s_{\textrm{U}}\) are determined by fitting to the experimental data of Johnson & Christy (dashed lines). When the number \(N_{\textrm{G}}\) of Lorentz poles (i.e., Gauss quadrature points) increases, both, real and imaginary part, converge to the full model. An approximation with \(N_{\textrm{G}}=3\) Lorentz poles works well for wavelengths larger than 450 nm and a significantly larger number of Lorentz poles is required for obtaining accurate results below 450 nm. The values for the fit parameters are listed in Table 1

Full size image
Table 1 Values of the fit parameters for the alternative discretization of the parabolic two-band model for gold
Full size table

The results of Fig. 3 show that this approach leads to a further improvement of the approximation. The single effective Lorentz pole approximation now works well for wavelengths larger than 550 nm, while the approximation based on three effective Lorentz poles works well down to about 450 nm. For even smaller wavelengths, a significantly larger number of effective Lorentz poles is required for accurate results where the discretization converges to the results displayed in Fig. 1. However, compared with the original discretization discussed in Sect. 3.1, the unphysical peaks in the complex refractive index are less pronounced within this alternative discretization scheme of the IBT-integral, thus somewhat alleviating but certainly not eliminating the afore-discussed problem of unphysical peaks in the complex refractive index.

Finally, we would like to emphasize that, on very general grounds, our approach of using only 6 fitting parameters and systematically increasing the number of Lorentz poles for improving the accuracy is preferred over the traditional approach (see, for instance, Ref. [5]) that increases the number fitting parameters through using more and more Lorentz poles. While in terms of computational effort for comparable accuracy, our approach might be only slightly more efficient, the issue here is that models with ever larger numbers of fitting parameters have an inherent tendency to produce less reliable results [29].

4 Conclusions

In summary, we have developed an efficient, electronic band structure-informed model for the linear dielectric response of simple metals with an emphasis on treating interband transitions and the model’s incorporation into time-domain Maxwell solvers. We have illustrated this approach by way for arguably the most important plasmonic material, gold.

Our model is based on the notion of separating the metal electron’s hydrodynamic (intraband) characteristics from the interband transitions by way of an appropriate canonical transformation. The resulting effective modeling of interband transitions leads to an excellent agreement with experimental data by fitting only four free parameters.

In addition, we have demonstrated two closely connected approaches how this model can be efficiently discretized for time-domain Maxwell simulations via an arbitrary number of effective Lorentz poles for standard auxiliary differential equations. This provides a systematic way of improving the results without introducing additional fit parameters. With three effective Lorentz poles, these models work well for wavelengths larger than about 450 nm and significantly more effective Lorentz poles are required for accurate results below 450 nm. In this latter wavelength range, peaks associated with the effective Lorentz poles may lead to unphysical results. Consequently, the aforementioned systematic way of improving the discretization (which reduces the amplitudes of the Lorentz peaks) together with the fact that in all instances only six physically meaningful parameters are required (two for the Drude permittivity and four parameters for the interband transitions) represents a useful tool for numerical computations.

Furthermore, we should like to mention that our model naturally allows more advanced treatments of the hydrodynamic (intraband) contributions well beyond the Drude model. Specifically, the Drude permittivity may simply be replaced by any of its nonlocal extensions [17]. Due to the additivity of the different contributions to the linear dielectric constant, the extension of our approach to the treatment of frequency ranges and/or materials where the underlying electronic band structure suggests that, for a single band, different high-symmetry lines in the Brillouin zone have to be resolved and/or additional bands become relevant is straightforward. For instance and as alluded to Sect. 2.3, the accurate modeling of interband transitions in ultrafast pump-probe settings [26] suggests a separate treatment of the transitions along the \(\Gamma\)-L and the \(\Gamma\)-X line. Naturally, this also increases the number of fit parameters. Especially in the context of hot-electron dynamics in metals, a further extension of our model is conceivable. Specifically, considering a non-equilibrium occupation \(p_{l \vec {k}}\) in the general expression for the interband transitions, Eq. (4) would allow the treatment of hot-electron dynamics and thermalization effects. As a matter of fact, a similar approach has already been described by Saavedra et al. [30], albeit within a standard Drude–Lorentz fitting approach where the Lorentz-pole oscillator strengths \(f_m\) (see Eq. (1)) have been equipped with a temperature dependence.

As a result, our model will be very useful for accurate modeling of plasmonic nanostructures in regimes where both, nonlocal and interband transitions, are important.