Terahertz and infrared nonlocality and field saturation in extreme-scale nanoslits

: With advances in nanofabrication techniques, extreme-scale nanophotonic devices with critical gap dimensions of just 1-2 nm have been realized. The plasmonic response in these extreme-scale gaps is signiﬁcantly aﬀected by nonlocal electrodynamics, quenching ﬁeld enhancement and blue-shifting the resonance with respect to a purely local behavior. The extreme mismatch in lengthscales, ranging from millimeter-long wavelengths to atomic-scale charge distributions, poses a daunting computational challenge. In this paper, we perform computations of a single nanoslit using the hybridizable discontinuous Galerkin method to solve Maxwell’s equations augmented with the hydrodynamic model for the conduction-band electrons in noble metals. This method enables the eﬃcient simulation of the slit while accounting for the nonlocal interactions between electrons and the incident light. We study the impact of gap width, ﬁlm thickness and electron motion model on the plasmon resonances of the slit for two diﬀerent frequency regimes: (1) terahertz frequencies, which lead to 1000-fold ﬁeld amplitude enhancements that saturate as the gap shrinks; and (2) the near- and mid-infrared regime, where we show that narrow gaps and thick ﬁlms cluster Fabry-Pérot (FP) resonances towards lower frequencies, derive a dispersion relation for the ﬁrst FP resonance, in addition to observing that nonlocality boosts transmittance and reduces enhancement.

Under the push of an external electric field, electrons in a metal system tend to accumulate at its surface. In the classical local-response description, electrons do not interact with each other so that they can occupy identical states, generating an infinitesimally thin surface distribution layer, which in turn leads to a discontinuity in the macroscopic normal field component. However, Coulomb interactions and, more importantly, the Pauli exclusion principle prevent electrons from occupying the same state, leading to a strong electron-electron repulsion. Thus an electron gas resists compression by the applied electromagnetic field, giving rise to a smeared charge distribution. An external field applied at a point r 0 then causes a response not only at that point, but also in the neighboring volume, producing microscopically a nonlocal response effect. Because abrupt spatial variations in the materials properties at the interface do not lead to equally abrupt variations in charge/field distributions, the nonlocal effect poses one of fundamental barriers to the realization of ultimate field confinement and enhancement in the ultra-narrow gaps [27].
Previous studies on nonlocal effects mainly focused on visible and near-infrared regimes. Extending such investigations to longer wavelength regimes, in particular, millimeter-wave or terahertz (THz) and infrared (MIR) regime (wavelengths 0.75-20 µm), is important for two reasons. First, researchers have demonstrated substantial field enhancements in metal-insulator-metal (MIM) structures that are resonant in the MIR [20,33] and THz regimes [34][35][36][37][38][39]. Naturally, one could ask how the nonlocal effect limits field confinement and enhancement in such cases. Second, besides its fundamental importance, nonlocality is also critical for the accurate design and modeling of antennas, metamaterials, and sensors for practical applications. Despite the rapid growth of MIR and THz nanophotonics [40], the enormous mismatch between the long free-space wavelengths and the atomic-scale nature of nonlocality poses severe challenges to both experimental and numerical studies. As the wavelength increases, computational modeling of nonlocality in the MIR or THz regimes, where the mismatch between the wavelength (several microns to millimeters) and the Thomas-Fermi length (∼0.1 nm) can be larger than eight orders of magnitude, exceeding the capability of many commercial finite-difference time-domain (FDTD) or finite-element method (FEM) solvers.
The hybridizable discontinuous Galerkin (HDG) method is a high-order accurate finite-element numerical scheme that can address aforementioned challenges. Our previous work showed the capability of the HDG method for efficiently modeling long-wavelength MIR and THz radiation in free space and plasmonic resonators [15,[41][42][43]. We have also equipped the HDG method with the hydrodynamic model for metals [44,45]. In this paper, we leverage this novel framework to address the daunting challenge of modeling THz nonlocality, which requires resolving the nonlocal charge distribution in the metal and the interaction of a millimeter-long wave within single-digit nanometer gaps. Here we apply the HDG method to study light transmission through a single slit in the MIR and THz regimes. We quantify the impact of nonlocality in field enhancement, localization, as well as field saturation phenomena [46,47].

Results and discussion
The structure analyzed in this article consists of an infinitely extended gold film of thickness T suspended in free space, with an infinitely long aperture or gap of width G, illuminated from below by a Gaussian beam orthogonally polarized to the slit, see Fig. 7 (in the Appendix). For the simulations we consider film thicknesses between 150 and 500 nm and perfectly smooth gaps ranging from 1 to 50 nm, which will allow us to identify the influence of nonlocal electrodynamics without accounting for purely quantum phenomena, such as electron tunneling, that gain importance below ∼1 nm [23,31]. In addition, two different frequency regimes are considered: the low terahertz regime (0.35-1.5 mm wavelengths) and the near to mid infrared regime (0.75-20 micron wavelengths), where Fabry-Pérot (FP) resonances are excited. Simulation parameters and results exhibit differences between these two regimes, hence they are treated separately in the discussion.

Simulation methods
The results presented in this work couple two sets of equations: (1) the time-harmonic Maxwell's equations are solved in free space; and (2) the time-harmonic Maxwell's equations augmented with the hydrodynamic model for metals [27,48] are solved in the gold film. For illustration purposes, the influence of the nonlocal parameter β is accounted through an auxiliary b factor described as β 2 = bβ 2 0 , where β 2 0 = 3/5v 2 F is the widely used value of the nonlocal parameter in terms of the Fermi velocity v F [49]. Note that b = 0 corresponds to a purely local description of electron motion in metals based on the Drude model.
The numerical experiments in this article have been carried out with the hybridizable discontinuous Galerkin (HDG) method, a fast and high-order accurate numerical method that has been developed to solve the 3-D time-harmonic Maxwell's equations for both the local [41] and the nonlocal [45] model. The HDG method is specifically tailored to resolve tight field localizations among multiple length scales, and has therefore been used to successfully simulate ultra-narrow coaxial apertures periodically patterned on optically thin gold films [15,42,43,45]. The HDG method belongs to the class of discontinuous Galerkin (DG) methods, which are unstructured, high accurate, locally conservative, exhibit low dissipation and dispersion and are high-order, meaning the numerical error in the approximation can be made insensitive to the mesh discretization, as opposed to other finite-difference time-domain (FDTD) or finite element method (FEM) commercial solvers. For the aforementioned reasons, DG methods have been extensively used to simulate plasmonics phenomena [50][51][52][53][54]. The HDG method differs from other DG methods in the introduction of auxiliary variables defined on the faces of the discretization -also known as hybrid variables, hence the name hybridizable-, which enable the solution of smaller linear systems of equations and attain optimal convergence rates. Thus, the HDG method is a strong candidate for solving extreme nanophotonics phenomena [44,45,52]. For a more thorough review on the HDG method and its implementation to the infinite slit problem, we refer the reader to Appendix A.
Leveraging the symmetry of the structure, we simulate only the 2-D cross-section of half of the slit, thus prescribing perfect electric conductor conditions on the boundary along the symmetry plane. Perfectly matched layers (PMLs) [55,56] are used to truncate the computational domain on all non-symmetry boundaries, precluding the outgoing waves from re-entering the domain. The structure is illuminated at normal incidence from below with a monochromatic Gaussian beam of waist w 0 = 3λ. Since we are interested in different frequency regimes, meshes of distinct resolution and degree of anisotropy have been used. In all cases, the solutions are approximated with cubic polynomials on quadrilateral elements. The computational domain is chosen such that the paraxial approximation of Gaussian beam remains valid and that the outgoing waves are totally attenuated in the PMLs region. The quantities of interest are the amplitude enhancement of the electric field along the upper tip of the aperture, hereinafter referred to as Q, and the transmittance. The expressions used to compute these outputs, along with the discussion on the meshes, governing equations, illumination settings, boundary conditions and parameter values are presented in the Appendix A.

Nonlocal electrodynamics in the terahertz regime
The single slit exhibits extraordinary optical transmission and huge field enhancements as the frequency approaches zero [11,46,57], and the enhancement has been shown to saturate for shrinking gap size [47,58], with a limiting saturation enhancement Q lim that is inversely proportional to the film thickness and the wave frequency. These studies have been carried out for a slit on a micrometer-thick film described with both the Drude model and as a perfect conductor, and here we extend the study to metal films thinner than a micron and for gaps ranging between 1 and 50 nm and modeled with both the Drude and the hydrodynamic model. The simulation of electromagnetic wave propagation for low terahertz frequencies is computationally challenging due to the huge mismatch in length scales -mm long wavelengths are squeezed into nanometer-wide apertures. In addition, simulations with the nonlocal model require a very fine resolution near the metal-vacuum interfaces so as to capture the extreme localization of the electron density and decay of the electric field. For each wavelength, we simulate a computational domain of 12λ × 6λ surrounded by a 2λ-thick PMLs in each direction, ensuring the relevant plasmonic phenomena are properly represented. The longest wavelength considered for this regime is 1.5 mm (0.2 THz).
The simulated field enhancement |E x |/|E inc | results for the 2 nm aperture, 300 nm film at 0.2 THz are shown in Fig. 1(a). The enhancement is almost constant along the slit with an extreme localization on the upper corners, thus we focus there to better appreciate the differences between the local and nonlocal model, namely a lower enhancement, a more diffused plasmonic hotspot and a smoother metal-dielectric transition. The induced electron density is depicted in Fig. 1(b) in logarithmic scale, where we observe the decay of several orders of magnitude within few Ångstroms, forming a boundary-layer structure at the metal-dielectric interface. A detail of the mesh used to compute these solutions is also shown in Fig. 1(b), and we refer the reader to Fig. 8 (in the Appendix) for additional images on the mesh dimensions, characteristics and resolution. To further investigate the nonlocal effects at the metal-dielectric interface, we show the field enhancement profile along the middle of the slit in Fig. 1(c), for 300 nm film at 0.2 THz as before, according to the axes at the inset. The b = 0 case corresponds to describing the metal with the Drude model, where the metal acts like a mirror and thus allows minimal penetration of the electric field, resulting in an abrupt decay of the enhancement at the interface. For positive b the hydrodynamic model allows the propagation of a new longitudinal mode that relaxes the field enhancement decay inside the metal, forming a boundary-layer structure characterized by [59,60]. Note that this is much smaller than the usual (transverse mode) skin depth. For γ = 0, ε ∞ = 1 and ω ω p , the depth at which the field has decayed by a factor of 1/e equals δ L β/ω p [61]. Our simulation results coincide with the analytic decay, for instance with b = 1 and G = 1 nm, we have δ L,exact = 0.786 Å and δ L,simulation = 0.783 Å. The depth where the field enhancement |E x |/|E inc | plateaus is roughly equivalent to 9, 12.5 and 15.5 Å for b = 0.5, 1, 1.5 respectively. Consequently, the effective aperture width is enlarged [62] compared to the classical local approximation, where the electric field decay occurs sharply at the interface, and larger values of β translate into increased skin depth inside the metal.
To study the electric field structure we monitor the field enhancement on the upper tip of the slit Q -details available in Appendix section A.2. The enhancement values for all gaps are depicted in Fig. 2(a) only for the nonlocal model (b = 1) and the 150 and 500 nm film -results are qualitatively identical for the local model-where we see that narrower gaps increase field enhancement. Nonetheless, we also observe that for a fixed thickness and frequency, shrinking the gap does not lead to arbitrarily large enhancements, but rather approaches a limiting value. This phenomenon is known as field enhancement saturation [47], and the limiting enhancement has been described as Q lim = λ∆H 2πT , where T is the film thickness and ∆H is the difference in magnetic enhancement between the upper and lower boundaries of the film under the same illumination conditions when there is no gap, that is the structure is an infinite gold film without an aperture [47]. In our simulations, we have found the limiting enhancement to be almost constant across models Q b=0 lim /Q b=1 lim ≈ 1, mainly because the nonlocal effects are not noticeable in the absence of a gap where light is channeled.
For each film thickness and frequency, we evaluate its limiting enhancement and compare it with the actual enhancement Q observed when simulating the slit for various gap sizes. The results for f = 0.2 THz are shown in Fig. 2(b) for both the local and the nonlocal model with b = 1. Note that as gap decreases and film thickness increases, the ratio Q/Q lim approaches 1, but at a different rate amongst electron models. For the local model, Q/Q lim is between 95-97% for the 1 nm gap for all thicknesses, whereas for the 10 nm gap saturation ranges from 66% (T = 150 nm) to 86% (T = 500 nm), thus differences among thicknesses vanish as the gap reduces. This effect is not as pronounced for the nonlocal model, see for instance the 1 nm with saturation values between 83-94%, whereas the 10 nm gap exhibits saturation values in the range of 64-85%, very similar to those of the local model. Due to the electron density profile that smears inside the metal, the gap is effectively enlarged in the nonlocal calculations, therefore resulting in lower enhancements and lower saturation values.
The numerical results are therefore in good agreement with the theoretical behavior of the perfectly conducting slit, and show that for this quasi-static regime, extraordinary optical transmission and field enhancement saturation may also be observed if a more realistic physical model for electron motion is used, albeit achieving lower saturation values than for perfect electric conductor and Drude model representations of the metal. Note the rate at which enhancement approaches Q lim for shrinking gap width depends on film thickness. (c) Saturation field enhancement spectra for different film thicknesses, where we observe that the nonlocal model predicts lower saturation enhancements than the local model since the gap is effectively at enlarged. Differences between models are more acute for thinner films.

Nonlocal electrodynamics in the infrared regime
For a given slit gap and film thickness, increasing the frequencies towards the infrared regime excites several resonant Fabry-Pérot (FP) modes along the vertical direction of the slit, that lead to extraordinary optical transmission and large field enhancements [63]. In this section, we investigate the impact of nonlocality, gap width and film thickness on the FP resonances, devoting special attention to the first Fabry-Pérot resonance (FP1). Since the simulations span frequencies between 30 and 600 THz, for each wavelength we prescribe a computational domain of 12λ × 6λ, that together with a PMLs of width 2λ in each direction ensures a faithful representation of the underlying physics.
Firstly, we identify the separate influence of the aperture width and the film thickness on the location and distribution of the subsequent FP resonances in the mid IR, using the nonlocal b = 1 model. Setting the thickness to 500 nm, we compute the field enhancement profile for the 1 and 10 nm gap, see Fig. 3 (solid and dash). Narrower gaps squeeze FP resonances to lower frequencies and exhibit larger enhancements, see that below 300 THz the first seven FP resonances of the 1 nm gap are excited compared to only the first two FP resonances of the 10 nm gap. Then, we set the gap to 1 nm and evaluate the field enhancement profile for a 150 and 500 nm gold film, see Fig. 3 (solid and dotted). Thicker films squeeze FP resonances to lower frequencies and exhibit lower enhancements, see that below 300 THz the first seven FP resonances of the 500 nm film are excited compared to only the first two FP resonances of the 150 nm film. Hence, we can conclude that narrower slits and thicker films cluster the FP resonances towards the lower end of the spectrum.
We now focus solely on FP1, whose peak field enhancement Q and transmittance as a function of the film thickness are depicted in Fig. 4 for several gaps using only the nonlocal model (b = 1), since the results for the Drude model are qualitatively similar. For a given thickness, smaller gaps exhibit more enhancement (since there is more confinement) and less transmittance (since there is less open area). If we account for thickness variation, the simulation results show that thinner films increase the field enhancement for gaps below 3 nm, and that the effect is reversed for gaps larger than 3 nm. On the other hand, transmittance is determined by both absorption -that augments as we consider thicker films-and field enhancement -whose behavior has just been described-and interact as follows: for gaps below 6 nm, film thinning from 500 to 150 nm leads to an increase of field enhancement (1.6-fold for the 1 nm gap) that is positively combined with absorption reduction, resulting in a transmittance boost (8-fold for the 1 nm gap); however, for gaps larger than 6 nm we observe that as the film thickens from 150 to 500 nm the increased absorption is outbalanced by the boost in field enhancement (3.3-fold for the 50 nm gap), which leads to more transmittance (2.5-fold for the 50 nm gap) for thicker films. To conclude, the combination of field enhancement and absorption, which interact constructively for narrow gaps and destructively for large gaps, drive transmittance. Incorporating nonlocal electrodynamics in the simulation of the slit induces a blue-shift in the FP resonances with respect to local simulations, and the shift does not depend on the thickness of the gold film, see Figs. 5(a) and 5(b) for FP1. The resonance blue-shift in noble metals has already been identified in previously published work [27,31,64] as a relevant consequence of nonlocal effects. Furthermore, nonlocality is responsible of reducing the field enhancement and boosting the transmission, see Figs. 5(c) and 5(d). This behavior is directly linked to the penetration depth δ shown in Fig. 1(b). Indeed, the sharp decay of the electric field inside the metal is smoothed by nonlocality (b > 0), effectively enlarging the gap seen by the incident wave. This enlargement leads to lower enhancements (less confinement of the wave) and increased transmissions (larger open-area to propagate through), and since δ does not depend on the original aperture width, its effects become less noticeable as the gap widens -5 nm seems to be the cutoff for pronounced nonlocal effects. If instead of representing the resonant wavelength against the gap size as in Fig. 5(a), we use the film thickness T, a linear dependence λ = α 0 + α 1 T arises. To that end, we use the resonant wavelength values from the simulations to regress the coefficients (α 0 , α 1 ), and show both the data (cross) and the linear fit (solid line) in Fig. 6(a), with coefficient of determination R 2 > 0.9999 for all gap widths and b values. The coefficient of determination is the proportion of variance in the resonant wavelength that may be predicted from the thickness value. Since only a tiny percentage of variance is not explained by the linear models, we may safely conclude the relationship between resonant wavelength and thickness is linear. The values of the linear fit parameters for several gaps and b values are collected in Table 1 (in the Appendix), along with information on how R 2 is computed. Moreover, the dispersion relation for FP1 can be derived following the results of [4] for a metal-insulator-metal waveguide. Namely, for each unique gap the wavelength of the propagating plasmon λ p = 2π/k along the slit at FP1 is twice the thickness, hence the wave vector reads k = π/T. The dispersion relation is therefore expressed as λ = α 0 + α 1 π/k, and is shown in Fig. 6(b) for several gaps and electron models.  Table 1. Linear fit coefficients for λ = α 0 + α 1 T , where T is film thickness in µm and λ is the FP1 resonant wavelength in µm.

Conclusions
This computational study shows that nonlocal description of electron motion in metals has a significant effect on the optical properties of an ultra-narrow slit in the THz and MIR regimes. First, for frequencies in the low-THz regime we extend the field saturation studies in [47] to both the local and nonlocal model, as well as to thinner films and narrower gaps. We have numerically verified that, for a given frequency and film thickness, field enhancement saturates as the gap is reduced and approaches the limiting enhancement value proposed in [47]. We have demonstrated that nonlocal saturation is less pronounced, since the theoretical limiting enhancement values are similar in both cases whereas nonlocal enhancement is significantly weaker due to the effective gap enlargement. In addition, we have numerically shown the impact of the parameter β on the field penetration at the metal-dielectric interface that characterizes nonlocality, and verified that the numerical values for the skin depth coincide with the analytical ones.
Furthermore, for MIR frequencies we observe a blue-shift of the gap plasmon for nonlocal calculations with respect to Drude model predictions. In addition, the nonlocal model predicts a decrease in enhancement and a boost in transmittance, due to the effectively increased gap width that stems from the nonzero penetration depth of the electric field into the metal. These results are consistent with previously published work [27,31], whereby a nonlocal representation of electron motion is indispensable to achieve accurate predictions for gaps below 5 nm.
On the numerical perspective, we have shown that the HDG method equipped with the hydrodynamic model for methods is a powerful scheme to simulate plasmonic structures. Indeed, it enables us to resolve the squeezing of millimeter-long wavelengths into nanometerwide cavities, along with the atomic-scale charge distribution layers that characterize nonlocal electrodynamics. This is even more dramatic for the THz regime, where both high order solvers and adaptive anisotropic meshes are required to accurately represent the solution in the presence of such lengthscale mismatch. The results in this paper showcase the strength of the HDG methods towards more efficient simulations of extreme-scale nanophotonic phenomena -terahertz nonlocality. Future applications of the HDG method include investigations of THz optoelectronics [39], nonlocality in the nonlinear regime [65] and extreme field confinement via polaritons in two-dimensional materials [66][67][68].

A.1. Governing equations
The hybridizable discontinuous Galerkin (HDG) method [41,43,45] is particularly suited to simulate wave propagation problems for a number of reasons: (1) it can be used on unstructured and anisotropic meshes, thus enabling the simulation of complex geometries; (2) it has low dissipation and dispersion, and the numerical error may be reduced to the level where it is insensitive to the geometry discretization; (3) it gives rise to a linear system that contains a reduced number of degrees of freedom defined on the faces of the discretization; and (4) enables the natural treatment of both boundary conditions and material contrasts at interfaces.
In this section, we describe the equations used in the simulation of the infinite 2-D slit shown in Fig. 7(a). The HDG method was used to simulate the full 2-D electromagnetic response on a computational domain Ω ⊂ R 2 that contains the slit. The domain Ω can be split according to whether the material is free space Ω f or a metal Ω m , thus Ω = Ω f ∪ Ω m , with boundary ∂Ω = Γ.
In Ω f the scattering is governed by the time-harmonic Maxwell's equations for a single frequency ω, that is In the metal Ω m , the time-harmonic Maxwell's equations are coupled with the hydrodynamic model for metals, namely ∇ × E − iωH = 0, where J is the internal current and ρ the internal charge density. Parameters (ε ∞ , ω p , γ, β) represent the metal permittivity, plasma frequency, the electron collision rate and the nonlocal parameter respectively. We use ε ∞ = 1, ω p = 9.02 eV, γ = 0.02678 eV from Ordal [69], and for the nonlocal parameter we have β 2 = b 3/5v 2 F , with v F = 1.39 · 10 6 m/s being the Fermi velocity of gold and b an auxiliary parameter to control the "amount" of nonlocality. The complete implementation details of the HDG method to solve (2) and (3) for metallic nanostructures may be found in [41,43,45].

A.2. Simulation details
The slit is illuminated from below with a Gaussian beam of wavelength λ that propagates upwards. The beam waist radius w 0 , which is the radius in the x-direction at which the beam intensity has dropped to 1/e 2 of its peak, is taken as 3λ, ensuring the paraxial approximation is valid. The length L is then taken as 4w 0 to guarantee the Gaussian beam intensity is negligible on the rightmost boundary. All in all, we prescribe h = 3λ, L = 12λ.
For computational purposes, we capitalize the symmetry of the structure and restrict the computational domain Ω to be half the slit, see Fig. 7(b). For a horizontally-polarized incident Gaussian beam (E inc , H inc ), mirror symmetry is imposed by prescribing a PEC condition (E × n = 0) at the symmetry axis Γ x , and perfectly matched layers (PMLs) [55,56] are appended to attenuate the outwards-propagating waves and simulate unboundedness. The PML dimensions are h 0 , h 1 = 2λ. Due to the attenuation effect of PMLs, we can safely employ PEC conditions on the remaining boundaries for simplicity. For the metal-dielectric interfaces, we need an additional condition to preclude the electrons from escaping the metal (no spill-out condition), given by J · n = 0.
The mesh used for the simulations is a structured anisotropic quadrilateral mesh consisting of 31K (resp. 42K) cubic elements for the low-THz (resp. mid-IR) regime. The meshes are highly anisotropic, and exhibit a significant concentration of elements near the tips of the slit and the metal-dielectric interface in order to capture the tightly localized features of the phenomenon, namely the atomically-thin charge distribution layer. The histogram of the mesh elements' areas, as well as several mesh details for 0.2 THz (the one with largest discrepancy between lengthscales) are depicted in Fig. 8. It can be observed that the mesh elements' areas span twelve orders of magnitude, thus highlighting the difficulties of the simulations and the necessity of high-order methods and anisotropic meshes to properly capture the solution. Transmittance is computed as the ratio between the transmitted power across a y-constant half-plane A located 1 micron above the gold film with x ∈ [0, 500] nm and the incident power across a y-constant plane A 0 located 1 micron below the gold film with x ∈ [0, 500] nm, namely The amplitude enhancement Q for on the upper tip of the gap S t similarly computed as the following ratio

Appendix B. Coefficients for linear fit of the dispersion relation
In this section, we report the coefficients α 0 , α 1 that govern the dispersion relation (λ, k) for the slit FP1 resonance, arising from the linear fit between resonant wavelength and film thickness, see Fig. 6 in the manuscript. The coefficients are collected in Table 1 for several gaps and electron models, b = 0, 0.5, 1, 1.5. As mentioned in the manuscript, the R 2 coefficient of the fit is above 0.9999 for all cases. To calculate this coefficient, for each gap and electron model b, we numerically obtain the resonant wavelength for each of the five thickness values considered (150, 200, 300, 400 and 500 nm), which can be expressed as five tuples (T i , λ i ) 5 i=1 . The coefficient of determination is therefore computed as where λ is the mean value of resonant wavelength for the particular gap and b, that is λ = 5 i=1 λ i /5, and L(T i ) is the model that we use, in this case a linear fit L(T i ) = α 0 + α 1 T i . Note the R 2 value ranges between −∞ (the fit can be arbitrarily bad) to 1 (perfect fit), and that R 2 = 0 corresponds to using a model that always predicts the mean value λ. We can thus conclude that our linear fit it is almost perfect.