No-slip boundary conditions for electron hydrodynamics and the thermal Casimir pressure

We derive modified reflection coefficients for electromagnetic waves in the THz and far infrared range. The idea is based on hydrodynamic boundary conditions for metallic conduction electrons. The temperature-dependent part of the Casimir pressure between metal plates is evaluated. The results should shed light on the"thermal anomaly"where measurements deviate from the standard fluctuation electrodynamics for conducting metals.


Introduction
The Universe is mainly filled with matter and radiation. While the former makes up a mass of roughly one Hydrogen atom per cubic meter, averaged over cosmological scales, quantum theory predicts since about one century ago that the zero-point energy of radiation, one half photon energy per electromagnetic mode, sums up to an energy per unit volume of (Adler & al., 1995) where the factor 2 accounts for two transverse polarizations, and Λ is a cutoff of the k-space available to the modes. Choosing this at the Planck length 1/Λ = ( G/c 3 ) 1/2 , one obtains an energy density a factor ∼ 10 123 above the energy equivalent of the matter content: the "wrongest formula in physics". It is intriguing to handle this discrepancy with an argument familiar from renormalization theory (and also used by Casimir (1948)): the vacuum energy scales in leading order with the volume of the system and can be subtracted by comparing two situations with the same volume -but differing in boundary conditions. However, what should set boundary conditions for the Universe as a whole and where? -perhaps the cosmological horizon, effectively considering the Universe as a bubble with Hubble radius c/H. Subtracting the energies, a term remains that scales with the surface of the bubble. Taking the short-wavelength cutoff Λ for k-vectors parallel to the horizon, the vacuum energy density in our Universe bubble becomes u → c Λ 2 (c/H) 2 ∼ infrared frequencies, the relevant length scales are below 100 nm, being determined by d. This calls for a reappraisal of methods that have been developed over the last century. Those coming from the context of infrared spectroscopy do not address the range of parameters k ω/c, since one is dealing with the response to long wavelengths λ = 2πc/ω. More relevant is work on electron energy loss spectroscopy (Verbeeck & al., 2005) where the fields correspond to the Coulomb potential of a moving charge. In that context, however, the focus has been on rather high frequencies (energies), even the surface plasmon resonance (in the visible or UV) being considered a low-energy feature (Zangwill, 1988). The situation is exacerbated by experiments addressing, in the distance range d ∼ 10 . . . 300 nm, the Casimir force (Chen & al., 2016;Bimonte & al., 2016) and non-contact heat transfer (Kloppstech & al., 2017;Cui & al., 2017): they give results that disagree with the standard theory of the fluctuating electromagnetic field (Lifshitz & Pitaevskii, 1980;Rytov & al., 1989).
We outline in this paper an improved, hydrodynamic approximation for the electromagnetic response of conduction electrons at a metallic surface. It is shown in particular that the classical Fresnel formulas based on a local dielectric function apply only in a specific range in the kωplane, as it happens also with other models including spatial dispersion García-Moliner & Flores (1979); Dressel & Grüner (2002). In our context, a kind of conspiration of scales has to be addressed. To fix a relevant set of parameters, consider the commonly used local Drude dielectric function and conductivity where ε b (possibly frequency-dependent, too) describes the response of bound electrons, the plasma frequency Ω p scales with the root of the conduction electron density, and τ is the scattering time. The latter can be determined from the DC conductivity σ 0 = ε 0 Ω 2 p τ . Typical parameters for gold at room temperature are Ω p = 9.1 eV and /τ = 27 meV (wavelength 46 µm, in the far infrared). The scattering rate separates the classical Hagen-Rubens regime (low frequencies) from the so-called relaxation regime 1/τ ω Ω p (see Sievers (1978) and Dressel & Grüner (2002, Appendix E) for more details). The first coincidence is that the typical thermal frequency is quite close to the Drude scattering rate k B T / = 0.94/τ (gold at room temperature). (For a detailed study of the behaviour of Casimir pressure and entropy at low temperatures, see Intravaia & al. (2010); Reiche & al. (2020).) The second coincidence is one of length scales. Recall that an electromagnetic field in the thermal frequency band penetrates into a metal in a diffusive way, leading to the characteristic length l m = ( D m /k B T ) 1/2 where D m = 1/(µ 0 σ 0 ) (µ 0 being the permeability) is the diffusion constant for magnetic fields (Jackson, 2014). The scale l m ∼ 20 nm explains the "thermal anomaly" of the Casimir pressure, namely that temperature-dependent corrections appear already at distances d much shorter than λ T (Boström & Sernelius, 2000;Intravaia & Henkel, 2009). The hydrodynamic model introduces another length in the same range, namely the electronic mean free path . From the Fermi velocity v F (gold: 1.4 × 10 6 m/s), we have = v F τ = 34 nm. This scale appeared already in the anomalous skin effect (Dressel & Grüner, 2002) that occurs when is larger than the classical penetration depth (c/ω)/ Im √ ε m . Its impact on the Casimir pressure has been studied using generalisations of the Fresnel reflection amplitudes (Esquivel & al., 2003;Svetovoy & Esquivel, 2006), although the modifications were found to occur only for p-polarized field modes. This polarization contains an electric field component perpendicular to the surface and probes the charge density profile at a metallic surface, whose characteristic scale is the Fermi wavelength 2π/k F = 5.2Å or the Thomas-Fermi screening length v F /Ω p = 1.0Å (Feibelman, 1982;Wegner & Henkel, 2020). We focus our hydrodynamic approach on the response to s-polarized fields whose electric field is parallel to the surface. To understand the basic idea, recall the Maxwell matching conditions for the tangential electric field in vacuum and metal, E (vac) = E (m). Using Ohm's law in local form j (m) = σ E (m), a nonzero current density right at the metal surface appears. This is not consistent with the no-slip boundary condition typical for the hydrodynamics of viscous fluids (Landau & Lifshitz, 1987). For electrons moving parallel to a surface, the no-slip condition takes into account, on length scales much larger than the Fermi wavelength, the scattering by surface roughness and by amorphous reconstructions of the sub-surface material. Note that this would not apply for atomically clean samples, but rather for metals kept in ambient conditions, as also suggested by experiments in the THz range (Laman & Grischkowsky, 2007). Our hydrodynamic calculations indeed predict that the electronic current density varies significantly in the sub-surface region on the scale of the mean free path .
A key parameter is the shear viscosity of the electron fluid. We use the observation of Conti & Vignale (1999) that the hydrodynamic (Navier-Stokes) equations can be phrased in the language of a visco-elastic medium. An elastic modulus and the viscosity are related, at finite frequencies, to the real and imaginary parts of the same mechanical response. It turns out that this response is encoded in the longitudinal and transverse dielectric functions of the charged Fermi gas. These are well-known within Lindhard theory (Lindhard, 1954) in the self-consistent field (or randomphase) approximation. By matching the long-wavelength expansion of these functions, extended to take into account collisions (Kliewer & Fuchs, 1969;Mermin, 1970;Conti & Vignale, 1999), we find a kinematic shear viscosity that scales, below the collision rate 1/τ , with v F = 2 /τ . This is formally a diffusion coefficient comparable in magnitude to the magnetic one D m = λ 2 p /τ because of the coincidence between the (reduced) plasma wavelengthλ p = c/Ω p = 22 nm and the mean free path .
An overview of our results is shown in Fig. 1 where the Casimir pressure (left) and the heat transfer (right) due to s-polarized modes is plotted. These modes give a sizeable thermal correction to the Casimir force (Boström & Sernelius, 2000;Torgerson & Lamoreaux, 2004;Intravaia & Henkel, 2009) and reduce the attraction prevailing between ideal reflectors. The p-polarisation does not contribute because of the efficient charge build-up at the metal surfaces (in the infrared, the dielectric screening is 1/|ε m | 1). Note in particular that the s-polarized modes alone account for nearly the entire difference between the measured Casimir force and the theory based on the local Drude approximation (open dots with error bars (Decca & al., 2005)). In the no-slip model, their repulsive contribution drops significantly compared to the local Drude model, so that the difference to observations gets smaller (Klimchitskaya & Mostepanenko, 2020a). The radiative heat current in Fig. 1 (right) is also reduced, but the data are well below the levels observed by Kloppstech & al. (2017); Cui & al. (2017). The outline of the paper is as follows. In Sec. 2, we motivate the formulas for the Casimir pressure between metallic plates at nonzero temperature and for the radiative heat transfer. After an introduction to the hydrodynamic approximation (Navier-Stokes equation) in Sec. 3.1, we solve the reflection/transmission problem at a conducting surface using the no-slip boundary condition (Sec. 3.2) and a modified surface current density (Sec. 3.3). A discussion of the reflection coefficients and the impact on the Casimir pressure is given in Sec. 4. Appendix A presents the derivation of the viscosity for conduction electrons based on the wave vector-and frequencydependent longitudinal and transverse dielectric functions.

Casimir pressure and boundary conditions
The famous formula by Casimir for the force per unit area between two ideally reflecting plates ( Fig. 2) reads where d is the distance and the negative sign denotes an attractive force. This formula ignores the physical properties of the plates, although its derivation requires that they become transparent in the far UV -to regularise the UV divergent vacuum energy. A more complete description is provided by the theory of dispersion forces. It describes electromagnetic fluctuations, both in the quantum and thermal regime, in and between macroscopic bodies and the ensuing interactions. If we restrict the discussion to distances d 1Å, it seems appropriate to use a continuum description and to describe the objects with the help of material equations, using the framework of macroscopic electrodynamics. This has been developed over the last twenty years into macroscopic quantum electrodynamics (Rytov & al., 1989;Buhmann, 2012a,b;Volokitin & Persson, 2017;Sernelius, 2018). Its basic idea is that the macroscopic response functions also determine the strength of the fluctuations produced by the bodies. By the very construction of this approach, there is no distinction to be made between virtual or real fields: the body's material is responding to an external field.
We focus on the simple geometry of two thick parallel plates a distance d apart, with the z-axis normal to the surfaces. They are kept at temperature T so that for a given (angular) frequency ω, the mean energy per photon mode is given by 1 2 ω coth( ω/2k B T ) = ω[ 1 2 +n(ω, T )] defining the Bose-Einstein distributionn. Due to translational and rotational symmetry parallel to the plates, the electromagnetic modes may be labelled by a two-dimensional k-vector k = (k x , k y ). Between the plates (vacuum), these modes vary with the wave vector In the second case, the modes are called evanescent, and they are localised to the vicinity of their sources. (Klimchitskaya & Mostepanenko (2020a) use the words "on-shell" ("off-shell") for propagating (evanescent) modes, respectively.) There are two transverse polarisations, usually called p (or TM) and s (TE). When a wave with polarisation µ is incident on the metal plate, it is reflected with amplitude r µ = r µ (k, ω). Multiple reflections between plate 1 and 2 can be represented by a geometric series At each reflection, a propagating photon imparts a recoil momentum of order k z onto the plate. Taking into account the reflection amplitudes and summing over all photon modes and their thermal occupation numbers, the electromagnetic stress normal to the surfaces yields the Lifshitz formula for the Casimir force per unit area (Lifshitz, 1956) This contains also the contribution of evanescent waves via the L-shaped path of the k z integral: it runs along the imaginary axis from i∞ to 0 and then to ω/c (recall the convention: negative P gives an attractive force).
Formula (5) is not suitable for calculating the pressure because its zero-temperature limit is plagued by the rapid oscillating factor e 2ikzd at high frequencies (along the real "leg" of the k zintegral). (The integral is physically cut off around the the plasma frequency where tabulated optical data rather than the Drude permittivity must be used.) Lifshitz shifted the ω-integration in the complex frequency plane to the imaginary axis ω = iξ which is possible because the integrand is built from response functions that are analytic in the upper half plane. The k z -integral is then taken from i∞ to iξ/c so that all exponentials become real and decay at large ξ. For finite temperatures, the integration is replaced by a summation over the Matsubara frequencies ω = iξ n = 2πink B T / , the poles of coth( ω/2k B T ), with the term n = 0 counting only one half: We have used here that the integrand f (ω) becomes a real function along the imaginary axis. Similar considerations have led Polder & Van Hove (1971) and Loomis & Maris (1994) to a formula for the heat current between two planar bodies of temperatures T 1 > T 2 separated by a vacuum gap of width d (Volokitin & Persson, 2017). For small gaps, the important contribution comes from evanescent waves, i.e. imaginary k z = iκ: where n(ω, T ) is the difference of Bose distributions. The contribution from real k z looks a bit different, Physical properties of the reflection coefficients, namely energy conservation for real k z (|r iµ | 2 ≤ 1) and passivity for imaginary k z (Im r iµ ≥ 0), ensure that the heat current is always oriented from hot to cold, consistent with the Second Law of thermodynamics. Note that in this approach, the concept of temperature has shifted from the field itself to its sources, namely currents and charges in the two bodies. A detailed discussion when the simple picture of two uniform temperatures T 1 = T 2 is applicable, has been given by Eckhardt (1982). The basic idea is that the material's heat capacity and thermal conductivity are sufficiently large so that the absorption of electromagnetic energy does not change its temperature. Similar arguments have been used to model transport in semiconductors at high fields (Ancona, 1995). Additional bodies in thermal contact are obviously also instrumental in maintaining the non-equilibrium setting.
In the following, our focus will be on the temperature-dependent part of the Casimir pressure and the radiative heat transfer. This is why we do not use the Matsubara sum: the thermal correction would be hidden in the difference between sum and integral [see Eq. (6)]. Using the argument principle rather than the Euler-MacLaurin formula to evaluate that difference, brings us back to the real-frequency integral (5). This makes one essential difference with respect to Klimchitskaya & Mostepanenko (2020a) where a modified surface response was also proposed, but the focus was on the behaviour of the zero'th Matsubara frequency ξ 0 = 0. Further comparison to that paper will be drawn in the Conclusion. Since we also focus on distances d much smaller than the thermal wavelengthλ T , the k z -integral is dominated by its imaginary leg (evanescent modes), thus providing numerically tractable expressions. A typical scale for the imaginary wavenumber is set by the inverse distance κ ∼ 1/d.
It remains in the following to analyse the reflection coefficients. If the plates are characterized by a local dielectric function ε m (ω) (a conductivity σ m (ω)), the Fresnel formulas can be used where is the decay constant inside the metal. (We take Re κ m ≥ 0.) These expressions arise from the dispersion relation k 2 −κ 2 m = (ω/c) 2 ε m in the metal and by matching the tangential components of the electric and magnetic fields at the vacuum-metal interface. In the following section, we derive a generalisation of these expressions using a hydrodynamic picture where the conduction electrons are modelled as a charged, visco-elastic medium.

Bulk
The key concepts for a hydrodynamic description are the density n and the velocity field v for carriers with mass m and charge e. The dynamics of the latter is given by nm(∂ t + v · ∇)v = f with the Navier-Stokes force density (Landau & Lifshitz, 1987) Here, 1/τ is the Drude scattering rate (see discussion in Appendix A), the compressibility is expressed via the velocity β, and the kinematic shear and bulk viscosities are η, ζ with ζ = ζ + 1 3 η. (This force density can be interpreted as a gradient expansion, assuming that n and v vary on large scales only. The hydrodynamic description does not resolve a microscopic scale like the Fermi wavelength.) We are interested in the linear response of conduction electrons to the electric field that splits naturally into a longitudinal (L) and transverse (T ) part E = −∇φ − ∂ t A (Coulomb gauge: ∇ · A = 0). Using the equation of continuity ∂ t n + ∇ · (nv) = 0, assuming all fields to evolve at a given frequency with exp(−iωt), and dropping second-order terms from the Navier-Stokes equation (11), we obtain The complex combination β 2 −iωζ can be interpreted as a dynamic modulus, using the language of visco-elastic media (Conti & Vignale, 1999). In a homogeneous medium where the fields vary with the wave vector q, these equations produce the longitudinal and transverse conductivities in the hydrodynamic approximation according to j L,T = n e e v L,T = σ L,T (q, ω)E L,T : where σ 0 = n e e 2 τ /m is the DC conductivity and n e e the equilibrium charge density. Note from the poles of these expressions for τ → ∞ how β determines the speed of longitudinal sound waves, while the transverse current behaves in a diffusive way with diffusion constant η.
The expressions (14, 15) provide a framework to actually find the hydrodynamic parameters for the electron fluid. In this paper, we focus on the seminal results of Lindhard (1954) for the conductivities in the self-consistent field (or random-phase) approximation, leaving a detailed study of exchange-correlation effects for later work (see Conti & Vignale (1999) for this more general approach). We incorporate collisions into the Lindhard functions in such a way that charge excitations relax to local equilibrium set by the electrochemical potential E F +eφ and their static limit is correctly reproduced (Kliewer & Fuchs, 1969;Mermin, 1970;Conti & Vignale, 1999). The Lindhard conductivities are not restricted by the hydrodynamic approximation and can be expanded for small q-vectors. As outlined in Appendix A, this procedure gives The expression (16) recovers the longitudinal speed of sound derived by Halevi (1995). Its real part crosses over from β = v F / √ 3 at low frequencies to 3 5 v F at high frequencies (isentropic limit). Its imaginary part is attributed here to the viscosities ζ, η of the electron fluid. The shear viscosity (17) turns out to be larger than the quantum scale /m and is plotted in Fig. 3 (limit q 1, black dashed lines). Its low-frequency limit is set by the diffusion constant . It scales β 2 T ∼ ω 2 at low frequencies, consistent with the picture that a liquid does not support low-frequency shear waves. Around the collision frequency, the response to shear changes from viscous to elastic and β T ≈ v F / √ 5 for ω 1/τ . The gray dashed lines correspond approximately to the longitudinal elastic parameters (see discussion of de Andrés & al. (1986) in Appendix A.2.3). The other curves in Fig. 3 illustrate corrections beyond hydrodynamics that appear when the dimensionless parameter qv F / (ω + i/τ ) ∼ 1. They become relevant on scales shorter than the mean free path, i.e., q > 1.
We make the following curious observation from Eqs. (16, 17): when the longitudinal speed of sound β and the bulk viscosity ζ are computed by subtracting the part involving the shear viscosity η, one obtains β 2 − iωζ = 1 3 v 2 F , a real constant. The bulk viscosity ζ thus vanishes (as also mentioned by Conti & Vignale (1999)). The dispersion of longitudinal sound waves found by Halevi (1995) originates entirely from the nonzero shear modulus 4 3 Re[−iωη(ω)]. We finally note a close coincidence of parameters. The penetration of s-polarised evanescent fields into a metal follows the decay constant κ m ≈ [k 2 − iωµ 0 σ 0 /(1 − iωτ )] 1/2 where the DC conductivity σ 0 can be expressed as a diffusion constant 1/µ 0 σ 0 =λ 2 p /τ that governs the spatio-temporal behaviour of low-frequency magnetic fields (Jackson, 2014). The kinematic shear viscosity η has also the dimension of a diffusion coefficient [see Eq.(15)] η ∼ v 2 F τ = 2 /τ . For the noble metal parameters, we have ∼λ p , and the diffusive behaviour of both types overlaps in space. The only way to formally isolate the local Drude model is the "dirty limit" where τ → 0 at fixedλ p , v F . A discussion of the opposite case, that is typical for low temperatures, is provided by Intravaia & al. (2010); Reiche & al. (2020).

Sub-surface region
We now apply the Navier-Stokes equations to the response of a metallic half-space to an spolarized field and compute the reflection amplitude r s . We assume that all fields vary ∼  Conti & Vignale (1999). We plot the kinematic viscosity Re η(ω) and the square of the shear wave sound velocity β 2 T = ω Im η(ω). The thick black dashed line corresponds to the hydrodynamic limit (small q). The gray dashed line is the result of de Andrés & al. (1986) (see Appendix A.2.3). The solid colored lines are obtained from the inverse conductivity σ 0 /σ T (q, ω) computed according to Kliewer & Fuchs (1969) and Conti & Vignale (1999), by subtracting the local limit 1 − iωτ and dividing by q 2 τ [see Eq. (15) and Appendix A.2.2]. The kinks appearing at ω ∼ v F q signal the onset of Landau damping (creation of electron-hole pairs). exp i(kx − ωt) with the x-axis parallel to the surface and the metal occupying the region z ≥ 0. It is easy to see that the s-polarisation gives transverse fields senkrecht (orthogonal) to the xzplane, we denote by v = v(z) and A = A(z) the corresponding components of v T and A. The no-slip boundary condition v(0) = 0 thus pertains to the tangential velocity of the electron fluid, while A, proportional to the tangential electric field, is actually continuous across the surface. This rule has the advantage of not needing a dimensional parameter (apart from the bulk viscosity fixed from the bulk behaviour). It nevertheless provides an "additional boundary condition" in the language of optics in spatially dispersive media (Dressel & Grüner, 2002). The approach of Klimchitskaya & Mostepanenko (2020a) is very different since they modify the q-dependence of the bulk dielectric function. Compared to Eqs. (14, 15), their corrections are linear in k rather than quadratic, the anisotropy being justified by the presence of the surface. Our approach is quite the opposite, since an explicit boundary condition enters into the description of the surface, while (deep) inside the metal, the bulk dielectric functions are applied, keeping their spatial dispersion within the scope of hydrodynamics.
Eq. (13) and the Maxwell equations yield the following equations of motion where we used the link between electron density and plasma wavelengthλ p = c/Ω p to re-write the current density µ 0 n e ev [first term of Eq. (19)]. In κ 2 b = k 2 − ε b (ω/c) 2 , the displacement current will give a negligibly small contribution. Since the fields must decay deep into the bulk metal, this system can be solved with the Ansatz and similar for A. The decay constants are given by The ratio between the eigenmode amplitudes v l , A l (l = 1, 2) is mv l /eA l = (κ 2 b − κ 2 l )λ 2 p [Eq. (19)]. It is essential to have two decay modes here, otherwise the boundary conditions v(0) = 0 = v(∞) would make the velocity vanish everywhere. For typical good conductors and ωτ ∼ 1, the two terms under the root are comparable. If the mean free path is not resolved, κ 1 ≈ (1−iωτ ) 1/2 / diverges (thin boundary layer) and κ 2 ≈ κ m = (κ 2 b −iωµ 0 σ 0 /(1−iωτ )) 1/2 is the decay constant in the local Drude model and the Fresnel equation (9).
The no-slip boundary condition fixes from Eq. (20) the ratio v 1 = −v 2 , so that only one free parameter remains. It is fixed by the amplitude of the field incident from the vacuum side. A convenient quantity is the ratio Z = −A/(dA/dz) (a length) evaluated at the surface. Since E y = iωA is a tangential electric field and B x = −dA/dz a tangential magnetic field, the ratio Z is actually proportional to the surface impedance of the metallic half-space. A quick calculation gives The "impedance" Z matches with the incident and reflected fields on the vacuum side, A(z) = A 0 e ikzz + r s e −ikzz . This yields the beautiful formula for the s-polarised reflection amplitude the main result of this section. (The second form becomes exact for ε b = 1.) In the local limit (vanishing viscosity), κ 1 diverges [Eq. (21)], and r s goes into the Fresnel formula (9). Corrections to this thus depend on the ratio κ 1 /k z . The hydrodynamic description is illustrated by the results in Fig. 4 where the current profile v(z) is shown for different choices of parameters. The values of the reflection coefficients are also given and compared to the local (Drude-Fresnel) result.

Reduced boundary layer conductivity
It is remarkable in Fig. 4 (bottom left) how the current distribution is "missing" a sub-surface sheet a few thick, when comparing to the local model that does not apply the no-slip boundary condition (light dashed lines). We outline in this section how this can be included into a modified boundary condition for the electromagnetic fields, using the excess field technique developed by Bedeaux & Vlieger (2002). This is actually a paradigmatic example of boundary layer approximations or multiple-scale expansions (Nayfeh, 1981;Bender & Orszag, 1978). The response of the charge density at a conducting surface has been recently analyzed in the same spirit by Mortensen & al. (2021). The excess field approach lumps the details about the behaviour of fields and currents in the surface region into a small number of response functions. The idea is based on a separation of scales where deviations from a homogeneous bulk material only occur in a thin region near the surface (sometimes called the selvedge (Sipe, 1980)). For simplicity, we focus on non-magnetic materials and neglect spatial dispersion in the bulk (far away from the interface). In the following, we provide a closer look at p-polarized waves because the calculations are more involved.
The central concept of an excess field is based on taking the difference between a smooth, microscopic field F (z), say, and its approximation F loc (z) that extrapolates the local-medium values down to the surface, Bedeaux & Vlieger (2002) define the "total excess" as the integral of this difference, F s = dz [F (z) − F loc (z)]. The integral typically converges even before the variation with depth of F loc (z) sets in, a manifestation of a separation of length scales. The excess may still depend on the coordinates x, y in the surface. In Fig. 4 (bottom left), the excess current j s = n e ev s would be the shaded area between the no-slip hydrodynamic and the local current profiles.
It now remains to connect the excesses to the fields outside the surface layer. Excess quantities play similar roles as surface charges and currents in macroscopic electrodynamics and determine jumps F | := F loc (z → 0 + ) − F loc (z → 0 − ) of the coarse-grained electromagnetic fields. For a non-magnetic system and a fixed frequency, an integration of the macroscopic Maxwell equations across the interface yields the boundary conditions (Bedeaux & Vlieger, 2002): where D is the displacement field, tangential components carry the index , andn is the unit normal pointing into the metal. For simplicity, the surface coordinates x, y have been suppressed everywhere. Note certain jumps that are absent from the ordinary Maxwell boundary conditions. Finally, material relations specific to the surface are needed to express the surface excesses by the bulk fields (Bedeaux & Vlieger, 2002). We focus here on the surface conductivity σ s and the surface resistivity R s , in order to capture the conductive properties of the selvedge. Introducing F = 1 2 (F loc (z → 0 + ) + F loc (z → 0 − )) as the average on both sides, one obtains the relations The time derivative of D s gives, of course, the excess current tangential to the surface, while E s z expresses a potential drop due to the normal displacement current. We proceed to solving the reflection and transmission problem for a p-polarized wave incident from the vacuum side. A plane-wave Ansatz ∼ exp i(kx − ωt) as in Sec. 3.2 leads to complex amplitudes E x (z) and E z (z) that away from the surface vary according to where E 0 is the amplitude of the incident field, n m = √ ε m the complex refractive index inside the metal, and κ m given in Eq. (10). The coefficients r p and t p are the reflection and transmission amplitudes. The corresponding magnetic fields can be obtained by using the relation ωB = q×E with the Snell-Descartes law giving the wave vectors q on both sides of the surface. The excess boundary conditions (29), (30) yield the set of equations which is solved for the reflection amplitude Compared with the Fresnel formula (9), this expression features three additional terms. The familiar Fresnel terms are paired with a mixed term containing surface conductivity and resistivity. Both contribute additional corrections as well. Interestingly, the terms including the resistivity are paired with the parallel component k of the wave vector, suggesting that this may be a minor correction in the long-wavelength limit (1/k much longer than the selvedge thickness). The same calculations can be done for the s-polarisation and lead to that depends only on the surface conductivity, since there is no normal electric field component in this case. For the surface quantities σ s and R s , we propose a Drude-like model with a relaxation time τ s and a length 0 that captures the thickness of the selvedge region: By comparing to the no-slip viscous model, we estimate 0 ∼ − , the negative sign translating the "missing" current due to the boundary condition (Fig. 4). The behaviour of the reflection coefficients in the far infrared is illustrated in Fig. 5: we note a reduction compared to the local approximation for r s . This may be attributed to a better impedance matching when the jump of the current density at the surface is reduced. A good agreement with the hydrodynamic model is found in the long-wavelength limit for the parameter combination 0 ≈ −0.36 and τ s ≈ 2τ . At larger k-vectors, some discrepancies occur. The p-polarization does not contribute significantly in the evanescent wave sector (away from the light line ω = ck).

Discussion of results
In Fig. 6, we plot in the kω-plane a "spectral representation" of the Casimir pressure (left) and the radiative heat transfer (right). Only the thermal contribution of the s-polarization is shown. The data are normalized to the maximum value (in the chosen domain) of the local (Drude-Fresnel) approximation. One notes for both quantities an upper limit ω < k B T / ∼ 1/τ , as expected from the Bose-Einstein distribution. The heat transfer data are shifted upwards in frequency due to the additional factor ω under the integrals (7, 8). The maximum in the kωplane is set by the magnetic diffusion constantλ 2 p /τ (dashed gray lines). A reduction of the pressure appears notably in the range set by the kinematic viscosity η(0) = 1 To summarize, in this paper we have extended the classic Fresnel formulas for the reflection of electromagnetic waves by a metal surface. Two methods have been used: a hydrodynamic description that captures the spatial dispersion of the metal's dielectric function, and a boundary layer technique introducing surface layers of charges and currents. Both methods build on the assumption that the electric current density right at the surface vanishes. This boundary condition corresponds to the behaviour of a viscous fluid, and mirrors the impact of surface roughness on the few-nm scale. The viscosity of conduction electrons was derived from a modification of the well-known Lindhard dielectric functions, taking into account collisions with impurities, but neglecting exchange-correlation effects (Conti & Vignale, 1999). Note that the no-slip condition is needed to solve the hydrodynamic Navier-Stokes equation that involves higher derivatives of the electronic velocity field. It corrects the spatial profile of the current density right below the surface, and the "missing current" is mapped in the boundary layer technique onto a tangential surface current sheet. An interesting consequence is a modification of the s-polarized reflection coefficient such that repulsive contributions to the Casimir pressure between metallic plates are reduced in both models. This brings theoretical predictions closer to the observed values, possibly pointing towards a physically motivated solution of the so-called "plasma vs. Drude" controversy.
Among similar attempts to modify the reflection amplitudes by taking spatial dispersion into account, we mention Reiche & al. (2020) and Klimchitskaya & Mostepanenko (2020a). For both, the starting point are the surface impedances of a metallic half-space based on a specular reflection boundary condition (García-Moliner & Flores, 1979;Ford & Weber, 1984) where the longitudinal and transverse dielectric functions appear. Reiche & al. (2020) use the nonlocal Lindhard theory corrected for collisions as in Appendix A, but also focus on the impact of Landau damping and the low-temperature behaviour of Casimir interactions. These authors have stressed as well that surface roughness on the scale of the mean free path may conflict with the specular reflection assumption. Klimchitskaya & Mostepanenko (2020a) invoke the breaking of translational symmetry due to the surface to introduce an anisotropic q-dependence into the dielectric function. This correction uses the same small parameter as our hydrodynamic approximation, but is otherwise quite different in form. The analysis presented here complements both approaches. The Navier-Stokes model allows to resolve spatial nonlocality on scales larger than the mean free path and predicts nontrivial variations in the sub-surface current when the no-slip boundary condition is applied. The excess field (boundary layer) technique collects some of the nonlocality into a modified surface response, while allowing for a simpler, local description of the bulk. This illustrates that the term 'surface' depends on the choice of length scales implicit in the formulation of fields and boundary conditions.
We conclude with a few remarks. -The hydrodynamic model has been used in metals long before, but the focus was almost exclusively on the longitudinal response (charge density waves). The corresponding speed of sound β(ω) was derived by Halevi (1995). Our analysis links it to complex visco-elastic moduli (Conti & Vignale, 1999) and provides an additional interpretation. The dispersion of β(ω) is actually due to the complex shear modulus of the collisional electron gas in the Navier-Stokes equation, while the bulk viscosity vanishes completely in the hydrodynamic approximation (Stokes hypothesis).
Finally, we expect that the boundary conditions considered here will also modify the surface plasmon dispersion relation (that appears as a peak in Im r p in Fig. 5). This is probably irrelevant to radiative heat transfer because it appears in the frequency range ω ∼ Ω p where thermal occupation is negligible. The plasmon dispersion has been studied since a long time (García-Moliner & Flores, 1979;Halevi, 1995) and depends on the spatial profile of the surface charge on the Thomas-Fermi scale v F /Ω p . We thus do not expect large modifications since the no-slip condition changes the current density on the much longer scale of the mean free path v F τ , but a quantitative analysis is beyond the scope of this paper. tics down to Atomic Scales (NOA)".

A Derivation of hydrodynamic parameters
For the ease of comparison to other work, we first write down the dielectric functions that follow from the hydrodynamic expressions obtained with the Navier-Stokes model (11) of the electronic liquid: (For simplicity, we put the background dielectric constant ε b = 1 in this Appendix.) We expect the hydrodynamic description to be accurate in the semiclassical and long-wavelength limits. This requires at least the regime q k F , the Fermi momentum. It is also apparent that the inverses 1/(ε L,T (q, ω) − 1) are polynomials in q 2 . We shall fix their complex coefficients by a corresponding expansion of the dielectric functions of the electron gas. For simplicity, we focus on the degenerate case (temperature much smaller than the Fermi energy E F = 5.5 eV) and on the self-consistent field (or random phase) approximation where the dielectric functions are given by Lindhard theory (Dressel & Grüner, 2002;Lindhard, 1954).
A technical challenge is to take into account collisions, since we expect their rate ∼ 1/τ to be comparable to ω. We apply results from Kliewer & Fuchs (1969); Mermin (1970); Conti & Vignale (1999) who combined the relaxation-time approximation with the longitudinal and transverse Lindhard functions. This provides to make contact with the local Drude conductivity as well. The same dielectric functions have been used by Reiche & al. (2020). Let us recall that collisions may involve different scenarios. Our focus is on impurity scattering that does not conserve total momentum, while carrier-carrier collisions do. Electron-phonon scattering is probably some intermediate case due to Umklapp processes Ashcroft & Mermin (1976). These mechanisms lead to distinct temperature dependences of τ . For an overview of the implications for Casimir-Polder interactions, see Reiche & al. (2020); Bordag (2017).
We match in the following the hydrodynamic parameters to the power series for the inverse susceptibilities whereω = ω + i/τ . Multiplying these expressions with iτ /ω, we obtain the normalized inverse conductivities σ 0 /σ L,T (q, ω) given in Eqs. (14, 15).

A.1 Expansion of Lindhard functions
We collect here, for the convenience of the reader, the Lindhard formulas with the dielectric functions of the homogeneous degenerate electron gas (Lindhard, 1954). Lindhard introduced the dimensionless variables The longitudinal dielectric function takes the form and the transverse is The (natural) logarithms are to be evaluated on their principal branchs, approaching the real frequency axis from above. This can also be denoted by u = (ω + i0)/qv F (hence the superscript 0 in Eqs. (43,45). The resulting imaginary parts appear in the domains u + z < 1 and |u − z| < 1 < u + z (and are positive there); they vanish for |u − z| > 1.
For the matching with the hydrodynamic expressions (40, 41), we perform a double expansion in the Lindhard variables: small z and large u. In this limit, the imaginary parts do not play a role, and we obtain a regular power series whose first few terms are Among the last two terms, the first one dominates if we restrict to frequencies ω v F k F = 2E F / . This is well justified for ω ∼ 1/τ and a collisional width /τ E F . This result is consistent with Lindhard's Eqs. (3.5, 3.10) apart from the order 1/u 4 which has been obtained, however, by Arista & Brandt (1984) in the same limit (see Table 1 there). Klimchitskaya & Mostepanenko (2020a) have also used the small parameter v F q/ω to add correction terms to the local dielectric function of the Drude model. Their correction is, however, of the first order and is anisotropic (the wavevector k parallel to the metal surface is used in place of q).
The corresponding expansion of the transverse dielectric function yields Lindhard's Eq. (3.19) contains a term z 2 /u 2 (which vanishes in our calculations) and stops before the order 1/u 4 . We conclude that the small parameter for corrections beyond the hydrodynamic approximation is (qv F /ω) 2 ∼ (q ) 2 if we focus on ω ∼ 1/τ . Fortunately enough, they appear with relatively small numerical coefficients.

A.2 Including collisions A.2.1 Longitudinal dielectric function
Kliewer and Fuchs (Kliewer & Fuchs, 1969), Mermin (Mermin, 1970), and Conti and Vignale (Conti & Vignale, 1999) have constructed the collisional form of ε L based on the requirement that the electron gas relaxes to a state defined by a shifted electrochemical potential E F + eφ where φ is computed self-consistently from the induced charge density. This argument can be carried out for a broad class of dielectric functions, and for the Lindhard function introduced above, it gives the formulaω We denote by the superscript τ the collisional form. In the first term on the rhs, the electric susceptibility is evaluated at the complex frequencyω, while the second one involves the static susceptibility that is responsible for the screening of a static charge density. The latter is evaluated from the Lindhard formula (43) by taking the limit u → 0 (approaching zero from the upper half of the complex plane), yielding Expanding for small q, one obtains Adding the two terms in Eq. (49), the hydrodynamic power series (40) becomes Note that this treatment of collisions is necessary to match the zero'th order term. The quadratic term yields the complex combination β 2 − iω(ζ + 4 3 η) spelled out in Eq. (16). Its low-frequency limit 1 3 v 2 F arises from the static susceptibility in Eq. (49). The second term of order q 4 is again negligible compared to the one before. Halevi (1995) found his formula for the speed of longitudinal sound waves by a similar hydrodynamic argument. The present formalism provides a visco-elastic view with a splitting into elastic moduli and viscosities (the latter appear as the imaginary part of Halevi's β 2 ). The connection between Navier-Stokes hydrodynamics and elasticity theory was also made by Conti & Vignale (1999).

A.2.2 Transverse dielectric function
To include collisions into the self-consistent field approximation, similar considerations are applied by Conti & Vignale (1999). If we assume that the total carrier momentum is not conserved (as it happens for impurity scattering), then the resulting susceptibility takes a form slightly simpler than (49) The power series for the comparison to the hydrodynamic form (41) is thus and its quadratic term yields the complex shear viscosity η(ω) given in Eq. (17). It is interesting to subtract the contribution of the shear viscosity η from β 2 − iω(ζ + 4 3 η). It turns out that only the static term survives The bulk modulus (as expressed by the real part β 2 ) is thus determined by the static density response, as expected for a compressibility. Its frequency dependence is negligible, consistent with the remark of Conti & Vignale (1999) that the relevant frequency scale is the much larger Fermi frequency E F / . The bulk viscosity ζ, however, vanishes: collisions do not contribute any losses when compressing the electron gas. This has been observed also for other models of the dielectric response by Conti & Vignale (1999). In particular, they went beyond the random phase approximation and included exchange correlation (XC) effects by dynamic local field factors. Apart from the self-consistent field, the Coulomb interaction between electrons is neglected when using Lindhard's transverse dielectric function. A residual footprint of XC effects may be encoded in the electronic lifetime τ , however, as soon as the Lindhard functions are generalized to a collisional model.

A.2.3 De Andrés & al
de Andrés & al. (1986) have argued that the construction of Kliewer & Fuchs (1969) for the collisional transverse dielectric function was in error because Eq. (53) could not reproduce, in the static limit, the weak diamagnetism of the electron gas. Their reasoning is based on the relation between the relative permeability µ and the dielectric functions. They apply Eq. (49) with the substitution ε → µ to derive a collisional permeability µ τ . The strict analogy between electric and magnetic response relies on the assumption that magnetic charges be conserved. It turns out that µ τ (q, ω) − 1 is for frequencies ω ∼ 1/τ close to its static value −Ω 2 p /(4k 2 F c 2 ) much smaller than unity. (The spin contribution to the magnetic response would only give small corrections (Lindhard, 1954;de Andrés & al., 1986).) Relation (56) then predicts that ε τ T,L are practically the same. We obtain a complex shear viscosity de Andrés & al.: η(ω) = iv 2 where the last term is negligible for ω, /τ E F . This is plotted in dashed gray in Fig. 3. Note that the 1/ω pole would yield a finite velocity for acoustic shear waves, quite unexpected for a liquid. Also the bulk viscosity ζ [see Eq. (55)] would be nonzero in disagreement with the general observations of Conti & Vignale (1999). We thus believe that the close analogy between electric and magnetic responses put forward by de Andrés & al. (1986) (magnetic charge conservation) is not warranted, at least in the low-frequency region.