Collective modes and quantum effects in two-dimensional nanofluidic channels

Nanoscale fluid transport is typically pictured in terms of atomic-scale dynamics, as is natural in the real-space framework of molecular simulations. An alternative Fourier-space picture, that involves the collective charge fluctuation modes of both the liquid and the confining wall, has recently been successful at predicting new nanofluidic phenomena such as quantum friction and near-field heat transfer, that rely on the coupling of those fluctuations. Here, we study the charge fluctuation modes of a two-dimensional (planar) nanofluidic channel. Introducing confined response functions that generalize the notion of surface response function, we show that the channel walls exhibit coupled plasmon modes as soon as the confinement is comparable to the plasmon wavelength. Conversely, the water fluctuations remain remarkably bulk-like, with significant confinement effects arising only when the wall spacing is reduced to 7 Å. We apply the confined response formalism to predict the dependence of the solid–water quantum friction and thermal boundary conductance on channel width for model channel wall materials. Our results provide a general framework for Coulomb interactions of fluctuating matter under nanoscale confinement.


I. INTRODUCTION
Fluids confined at the nanometer scale underly many technologically important processes [1,2], including filtration, seawater desalination [3,4], blue energy harvesting [5,6] and electrochemical energy storage [7].Yet, they started to be fundamentally investigated not more than 20 years ago, and their initial theoretical description was largely inherited from macroscopic hydrodynamics, with generic walls imposing the same boundary conditions regardless of their material composition [8].The first nanofluidic effects emerged from the realization that, at the nanoscale, one may not neglect the wall's surface charge [9], which results in coupled ion-fluid transport phenomena such as electro-osmosis and streaming currents [10].There has been, however, accumulating evidence in recent years that surface charge is not a sufficient descriptor for the nanofluidic solid-liquid interface.From fluids near conducting surfaces [11,12] to strongly interacting ions due to dielectric contrast [13][14][15], several studies pointed to the need of describing the solid walls at the level of their electronic properties.
It may indeed be expected that, close enough to a solid wall, the Coulomb potentials produced by charged particles in a liquid are screened by the dielectric response of the wall material: this effect has been termed interaction confinement [15].Charged particles in a liquid are, first and foremost, ions: interaction confinement produces effective Coulomb interactions between ions in nanochannels that are modified compared to bulk Coulomb inter- * nikita.kavokine@mpip-mainz.mpg.deactions, leading to a wealth of correlation effects [13,14].But a polar liquid such as water, even though electrically neutral after time-avergaging, has a molecular-level charge structure: water thus exhibits thermal charge fluctuations at terahertz frequencies and on a wide range of length scales [16] (termed hydrons [17]).The corresponding Coulomb fields are also subject to interaction confinement: they are dynamically screened by the thermal and quantum fluctuations of the electrons in the solid wall [18,19].This solid-liquid coupling has been shown to result in a "quantum" contribution to hydrodynamic friction, and in direct near-field energy transfer between the liquid and the solid's electrons [19,20].These effects bridge the gap between fluid dynamics and condensed matter physics, opening the way to engineering nanoscale flows with the confining walls' electronic properties [17,21].
Fluctuation-induced effects in nanofluidics have so far been studied at the level of a single planar interface.The relevant many-body electrostatics were conveniently described in terms of surface response functions: surface analogues of the dielectric function that had been widely used, for instance, in the field of plasmonics [22,23].Here, we introduce confined response functions, that generalize surface response functions to a 2D nanochannel geometry (Fig. 1a), providing a general tool for the treatment of Coulomb interactions in 2D confinement.As an illustration, we study the confined response of water, and of solid walls described as either graphene sheets or jellium slabs [24,25]: this allows us to predict the confinement dependence of solid-water quantum friction and thermal boundary conductance.
This paper is organised as follows.In Sec.II, we introduce confined response functions and link them to eigen-modes of the Coulomb potential in the 2D nanochannel geometry.In Sec.III, we compute the confined responses of specific media.We take the examples of a graphene sheet and a semi-infinite jellium for the solid; for the liquid, we study water in the framework of both force field and ab initio molecular dynamics (MD) simulations.In Sec.IV, we describe the influence of confinement on fluctuation-induced effects, specifically quantum friction and near-field heat transfer.To do so, we carry out the field-theory derivation of quantum friction directly in the confined geometry, which leads to natural emergence of the confined response functions.Finally, Sec.V establishes our conclusions.
Units and conventions.We set the Boltzmann constant k B = 1 (that is, we express the temperature in energy units), but otherwise use SI units throughout the text.In real space, we use the cylindrical coordinates r = (ρ, z).The interfaces are at z = 0 for a singleinterface and at z = ±h/2 for a confined channel.We use Fourier transforms for both ρ and the time but never for the z-direction.We use the following convention for the d-dimensional Fourier transform: (2π) d+1 F (q, ω)e iq•r−iωt .The charge densities are expressed in units of e and the electrical potentials include an additional factor e. We denote V (r) = e 2 /(4πϵ 0 r) the Coulomb potential which becomes V (q, z) = e 2 /(2ϵ 0 q)e −q|z| in Fourier space.

II. ELECTRIC RESPONSE OF INTERFACIAL SYSTEMS
A. Single-interface: surface response function We first briefly recall the widely-used concept of surface response function [22,23].Consider a semi-infinite medium occupying the half-space (z < 0).Given the electrostatic potential ϕ ext applied by an external source (an appropriate charge distribution) inside the medium, we wish to determine the potential ϕ ind induced by the medium in the half-space z > 0. The potential ϕ ext must solve the Laplace equation for z < 0. The physicallymeaningful (non-diverging) solutions are given by the evanescent plane waves where we have introduced the surface weight function F 0 : this seemingly cumbersome notation will be useful upon generalization to a confined geometry.Assuming that the medium has a linear charge density response function χ, the induced potential is given by (see Fig. 1b) where we have introduced the surface response function g 0 (q, ω) = − e 2 2ϵ 0 q dzdz ′ F 0 (q, z)χ(q, z, z ′ , ω)F 0 (q, z ′ ).
(3) It is worth noting that g 0 (q, ω) is a scalar.Since we are considering a linearly responding medium, the response to an evanescent plane wave at (q, ω) is an evanescent plane wave at (q, ω), so that the induced potential has the same z-dependence as the external potential, given by the weight function F 0 .In other words, the evanescent plane waves form an eigenbasis for the surface response.
Let us illustrate the role of the surface response function with a simple example.We consider a static point charge e at at distance z 0 from the interface.This charge produces an "external" potential where we have introduced a Fourier decomposition into evanescent plane waves.Accounting for the response of the medium, the total potential is The surface response function gives the contribution of the medium's polarization to the total potential.

B. Double interface: confined response functions
We now generalize this approach to a two-dimensional nanochannel geometry.We consider two interfaces at z = ±h/2 that define the channel, with the outside medium (|z| > h/2) being distinct from the inside medium (|z| < h/2).Later, we will specify that we consider water inside the channel, but we remain general at this point.
Contrary to the previous case, there is no longer a symmetry between the two media, and we need therefore to distinguish the responses of the inner and outer medium.All the applied and induced potentials still satisfy the Laplace equation, but, in both media, the subspace of harmonic functions with wavevector q and frequency ω is now of dimension 2: the response function at a given (q, ω) is then in principle given by a 2 × 2 matrix.We shall, however, express the potentials in a basis of even (symmetric) and odd (antisymmetric) harmonic functions, where the matrix turns out to be diagonal.We call these basis functions confined weight functions and denote them F s/a i/o , where s/a stands for symmetric/antisymmetric function and i/o for inner/outer medium (see Table 1).The amplitude of the basis functions is in principle arbitrary: it is chosen so that the confined response functions defined in the following reduce to the conventional surface response functions in the non-confined case.The responding medium can be either in the outer space (c&d) or in the inner space (e&f).The external potential can be either symmetric (c&e) or antisymmetric (d&f).We observe that the inter-solid interactions bring corrections by lowering the induced potential in the symmetric case (c) and increasing the induced potential in the antisymmetric case (d).

Single interface
Half-space TABLE I. Weight functions used in the definitions of the surface and confined response functions (Eq.( 3) and ( 8)).
Let us start with the response of the inner medium.We consider a generic external potential ϕ ext (q, z, ω) = ϕ s (q, ω)F s i (q, z) + ϕ a (q, ω)F a i (q, z) applied on the inner medium of charge susceptibility χ i .The induced potential is then In the outer space |z| > h/2, taking advantage of the definition of the confined weight functions (see Table 1), this reduces to where we have defined a generic confined response function: where the interaction takes place over the domain of definition of the weight function F c m ; m = i, o and c = s, a.To summarise, the inner medium responds to a perturbation of the form F s/a i by a potential of the form F s/a o in the outer medium with an amplitude −g s/a i (see Fig. 1c-d).
We proceed similarly for the response of the outer medium (with charge susceptibility χ o ).It responds to an external potential ϕ ext (q, z, ω) = ϕ s (q, ω)F s o (q, z) + ϕ a (q, ω)F a o (q, z) with an induced potential in the inner space |z| < h/2: The two components of the outer medium response are shown in Fig. 1e-f.To summarize, we have generalized surface response functions to a 2D nanochannel geometry by identifying the new eigenmodes of the Coulomb potential, given by the confined weight functions.The response to a symmetric (antisymmetric) weight function is a symmetric (antisymmetric) weight function with amplitude given by the confined symmetric (antisymmetric) response function, in the same way that the response to an evanescent wave is an evanescent wave in the single interface case.

C. Fluctuation-dissipation theorem and physical interpretation
To obtain a physical interpretation of the confined response functions, it is useful to relate them to equilibrium charge density fluctuations using the fluctuationdissipation theorem (FDT).For the charge susceptibility χ, the FDT reads [19] where f (ω) = 2T /ω for classical dynamics and f (ω) = h cotanh(hω/(2T )) for quantum dynamics.The structure factor S is defined as where n m is the charge density of the medium.Adapted to the confined response function defined in Eq. ( 8), the FDT becomes where S c is a confined structure factor defined as where A is the area of the interface.In Eq. ( 13), if the weight function F c m is symmetric, the charge density n m (z) can be replaced with (n m (z) + n m (−z))/2: the structure factor only counts the symmetric charge fluctuations.Similarly, if the weight function is antisymmetric, the charge density n m (z) can be replaced with (n m (z) − n m (−z))/2: the structure factor only counts the antisymmetric charge fluctuations.Thus, for the channel walls, the symmetric (antisymmetric) response function accounts for in-phase (out-of-phase) coupled modes.For the inner medium, the symmetric (antisymmetric) response function accounts for monopolar (dipolar) charge fluctuations.

D. Confined response function of the outer medium
We now provide an expression of the outer medium confined response function in terms of the usual surface response functions.We distinguish the contributions of the two solid walls to the potential induced in the inside medium: we denote them ϕ T ind and ϕ B ind for the "top" and "bottom" solids, respectively.We need to account for the solid-solid interactions: each solid responds to the external potential and to the potential induced by the other solid.Let us consider an external potential ϕ ext (q, z, ω) = ϕ o (q, ω)F c o (q, z) where the weight function The top solid induces a potential Combining these two equations, we obtain the total induced potential as in the inner space |z| < h/2.Comparing to Eq. ( 9), we deduce We find that the confined response functions reduce to surface response functions at wavelengths 1/q ≪ h: the deviation from the surface response function originates from the Coulomb interaction between the two walls.As expected on physical grounds, the inter-wall interaction reduces the in-phase (symmetric) response and enhances the out-of-phase (antisymmetric) response.However, the induced potential does not exceed the applied potential in confinement if it does not for a single interface: there is no confinement-induced overscreening (see ESI.1).

E. Interaction confinement
A first consequence of the wall electric response in a 2D channel is a modification of the effective Coulomb interactions between charged particles inside the channel: this effect has been termed interaction confinement [15].Two confined charges (for example, ions in water) interact not only directly, but also indirectly, through the polarization charges induced in the channel walls.Coulomb interactions near polarizable walls have been computed in various geometries and within different models for the wall polarizability.
In ref. [15], interaction confinement in a 2D channel geometry was addressed in the framework of surface response functions, which allowed for evaluation of effective Coulomb interactions between two ions in the channel mid-plane, for various models of the channel wall material.For example, with walls described by a Thomas-Fermi model, the ion-ion interaction is reinforced with respect to bulk water at small distances, and exponentially screened at large distances.
Our confined response formalism generalizes the results of ref. [15], allowing for the evaluation of the effective Coulomb interactions between arbitrary charge distributions.Any charge distribution can indeed be decomposed into an even and an odd part.The wall response to the even (odd) part is then given by the symmetric (antisymmetric) response function.In ESI.2, we provide expressions in terms of the confined response functions for the effective Coulomb potential produced by prototypical even (point charge) and odd (dipole) charge distributions.

III. CONFINED RESPONSE OF MODEL MEDIA
We now specialise to a 2D nanofluidic channel of height h filled with water, and we investigate the effect of confinement on the electric response of water and of model channel wall materials.

A. Outer medium: coupled plasmon modes
We consider two models for the channel wall material: a graphene monolayer and a semi-infinite jellium.The graphene surface response function is computed as detailed in ref. [19], starting from the charge susceptibility in the Dirac cone and zero-temperature approximations, at a chemical potential µ = 100 meV.The semi-infinite jellium is treated in the specular reflection approximation, as detailed also in ref. [19].In the jellium model, electrons are free in a uniform positive background, and are completely characterised by their chemical potential µ and effective mass m.We use µ = 180 meV and m = 0.1 m e as a model for a doped semi-conductor: this corresponds to an electron density parameter r s = 5.
For both systems, the surface response functions feature a sharply-defined surface plasmon mode and a broad particle-hole continuum (Fig. 2).The effect of confinement is most clearly visible on the plasmon mode: its energy is increases in the confined symmetric response and decreased in the confined asymmetric response.This is consistent with the physical interpretation outlined above: as the two solid walls face each other, in-phase (out-of-phase) charge density oscillations have an increased (decreased) energy cost due to the Coulomb interactions between the two solids.These interactions are significant only for charge fluctuations whose wavelength is longer than the confinement width h, so that the confined response functions differ from the surface response function only at small enough momenta q (see Eq. ( 17)).
We anticipate that the formation of coupled plasmon modes between the walls of a 2D nanofluidic channel will affect transport inside the channel, and particularly fluctuation-induced effects (see Sec. IV).

B. Inner medium: confined water spectra from simulations
We now turn to the confined response functions of water.Confinement may impact water charge fluctuations in two ways.First, the interaction between the two interfaces of the water slab is expected to result in a difference between the symmetric and antisymmetric responses.Second, the confinement-induced modifications of the water structure may intrinsically affect its charge fluctuations.
We determine the water response functions in the framework of molecular dynamics (MD) simulations.We carry out both force-field (FF) and ab intio density functional theory-based (DFT) simulations of water confined between two frozen graphene sheets, for various separations h between the graphene sheets.From the simulation trajectories, we compute the charge structure factor of water, integrated along z after multiplication by the weight functions summarized in Table 1, and determine the corresponding response functions through the fluctuation-dissipation theorem (Sec.II C).DFT simula- tions are required to capture intra-molecular modes: we use them to obtain the spectra at frequencies above 150 meV.At lower frequencies, we use the FF simulations to capture the contribution of inter-molecular modes, inaccessible with the short simulation times of DFT.Details of the numerical parameters and procedures are given in ESI.3.The results are presented in Fig. 3. Fig. 3b shows the water surface response function (obtained from the simulation at weakest confinement) at a fixed wavevector q 0 = 0.67 Å −1 and in the frequency range 0 − 200 meV.In this range, the water charge response can be decomposed into four modes: in order of increasing frequency, the Debye mode, the hydrogen-bond stretching mode, the libration mode and the OH-bond bending mode [16,26].We note that the OH-stretch mode (at around 450 meV) falls outside the studied frequency range.The confinement dependence of these modes could in principle be analysed in terms of their molecular origin; this is, however, beyond the scope of this article, and we restrict ourselves to a phenomenological description.
Overall, the water response functions are remarkably robust to confinement.Down to h = 1.4 nm, the lowerfrequency intermolecular modes remain unaffected.We observe, however, a slight red shift, and an increased oscillator strength for the bending mode.A significant effect on the intermolecular modes is visible only at 7 Å confinement.In the symmetric response, the Debye and hydrogen-bond-stretch modes are amplified, while the libration mode is suppressed; in the antisymmetric response, the libration peak is strongly amplified while the other modes are unaffected.As illustrated in Figs.3e-f, in the symmetric case, the perpendicular component of the applied electric field changes direction across the channel, while it maintains a constant sign in the antisymmetric case.This likely indicates that the libration mode is mostly excited by the electric field perpendicular to the interface, while the lower frequency modes are excited by the parallel component.

IV. CONFINED QUANTUM FRICTION AND HEAT TRANSFER
In this section, we discuss the effect of 2D confinement on fluctuation-induced interfacial effects -solid-liquid quantum friction and near-field radiative heat transfermaking use of the confined response function formalism developed above.

A. Single-interface quantum friction
We start by briefly summarising the physics of quantum hydrodynamic friction.The classical friction between a liquid and a solid is usually determined by the solid's surface roughness [27,28].However, it was recently shown that the classical contribution is supplemented by a fluctuation-induced or "quantum" contribution, due to the coupling of water charged fluctuations (termed hydrons) to electronic excitations within the solid [19].When undergoing quantum friction, a liq- uid transfers momentum directly to the solid's electrons.Similarly, a liquid may transfer energy directly to the solid's electrons: this is near-field radiative heat transfer [20,29].
In the case of a single solid-liquid interface, if the liquid is flowing at velocity v, it is subject to a quantum friction force F = −λAv with the quantum friction coefficient λ given by [19] Anticipating the generalization to the confined case, we have introduced the notation where the g's (e corresponds to the solid, w to the liquid) are generalized response functions computed with the weight function F (see Table 1).Similarly, if there is a temperature difference ∆T between the solid and the liquid, the solid-liquid heat flux is given by Q = κA∆T , where the thermal boundary conductance κ reads [20]

B. Confined quantum friction
We now generalize the above results to the 2D confined geometry presented in Fig. 4b.We wish to compute the total quantum friction force applied by the solid walls on the flowing liquid, and the total heat transfer rate between the liquid and the two solid walls.The liquid and the solid's electrons are described by their fluctuating charge densities n w and n e , which have a Coulomb interaction of the form ⟨Q⟩ = dr w dr e V (r w − r e )∂ t ⟨n w (r w , t)n e (r e , t)⟩, (23) where the integration over r w (resp.r e ) runs over the space occupied by the liquid (resp.solid).The correlation functions appearing in Eqs. ( 22) and ( 23) may be computed in perturbation theory with respect to H int , as has been detailed in ref. [19].Since the system is subject either to liquid flow or to a temperature gradient, the perturbative expansion needs to be carried out in the non-equilibrium Keldysh formalism [30].Ultimately, both the friction force and the heat transfer rate can be obtained in terms of the solid-liquid charge density correlation function χ ew , which, upon resummation of the perturbation series, is found to satisfy the following Dyson equation: Here * stands for convolution in space and time, multiplication by the Coulomb potential, and contraction of the Keldysh indices.A Feynman diagram representation of Eq. ( 24) is given in Fig. 4d.The Keldysh indices carried by the χ's are not important for the geometrical discussion that follows.We will therefore not write them out explicitly and we refer the reader to ref. [19] for further details: once the space-time convolutions have been dealt with, the computation is completely analogous to ref. [19].
In the single interface case, upon Fourier transformation in time and in space parallel to the interface, Eq. ( 24) immediately becomes a scalar equation for surface response functions.In the channel geometry, we need to introduce confined response functions, first as 2 × 2 matrices in the indices ν, ξ = T,B: where a,b=e,w, ϵ w = +, ϵ e = −, σ T = + and σ B = −.The space has been divided into three regions: the central region Z w , and the bottom (B) and top (T) wall regions: e ∪ Z T e .In the convolution over z in Eq. ( 24), summing over the top and bottom solids corresponds to summing over the indices B, T. Thus, in terms of the confined response functions, Eq. ( 24) becomes where the dot represents the matrix product, and g a ≡ g aa .Neglecting the off-diagonal terms in Eq. ( 26), we would recover two copies of the single interface Dyson equation, one for the top and one for the bottom interface.The off-diagonal terms represent cross-talk between the walls (for example, the top solid wall responding to a water fluctuation near the bottom wall), which is expected to vanish for weak confinement.We assume in the following that the top and bottom wall materials are the same, so that the system is symmetric under the mirror transformation z → −z.As a consequence, g TT = g BB and g TB = g BT .Therefore, all the confined response matrices can be diagonalised in the form where we have used which satisfies P −1 = P.We thus recover the definition of the confined response function in its eigenbasis, as introduced in Sec. 2. Multiplying Eq. ( 26) by P on the left and on the right, we obtain two scalar equations for the symmetric and antisymmetric response functions: g s ew (q, ω) = −g s e g s w + g s e g s w g s ew g a ew (q, ω) = −g s a g a w + g a e g a w g a ew (29) Once the Dyson equation is reduced to a scalar equation, we may follow the steps of refs.[19,20] to obtain the friction coefficient and thermal boundary conductance in the confined geometry: The confined response functions thus emerge naturally in the theory of fluctuation-induced effects in a 2D channel.

C. Diagrammatic approach to confined response function
Using the diagrammatic approach developed for the fluctuation-induced effects, we may interpret the confined response functions of the solid walls as surface response functions that have been renormalised by the inter-wall interactions within the random phase approximation (RPA).
We start from the intra-wall response function χ e (z, z ′ ), that has been renormalized by the intra-wall Coulomb interactions at the RPA level, according to the Dyson equation shown diagrammatically in Fig. 4c.These are identical in the top and bottom walls (χ BB e = χ TT e ), and, at this stage, there is no inter-wall response χ BT e , since we do not allow for electron tunneling between the walls.
In the presence of Coulomb interactions between the walls, we may introduce the exchange term Π = χ BB e * χ TT e (Fig. 4e).Still at the RPA level, the intra-wall response function is then renormalised according to (Fig. 4f) The inter-wall response is no longer vanishing, and satisfies (Fig. 4g) Fourier-transforming these Dyson equations as detailed above, we obtain relations between confined and surface response functions.Using in particular that (36) We thus recover Eq. ( 17), that we previously obtained from purely electrodynamic considerations.

D. Effect of confinement on the fluctuation-induced effects
We now investigate the effect of confinement on the quantum friction coefficient and thermal boundary conductance of water in a 2D channel, with walls made of either graphene or a semi-infinite jellium, with the parameters detailed in Sec.III.The results are presented in Fig. 5.
We first focus on the effect of inter-solid interactions and thus evaluate the fluctuation-induced effects using the single-interface surface response function for water (continuous lines in Fig. 5).Interestingly, we observe opposite trends for graphene and for jellium walls.In the case of graphene, for both friction and thermal conductance, the antisymmetric contribution is enhanced and the symmetric contribution is reduced with confinement.Indeed, both effects are governed by the coupled plasmon modes of the walls, and the out-of-phase mode has lower energy than the in-phase mode, thus making a larger contribution.For our jellium model, the plasmon is well above the thermal energy (around 300 meV), and the fluctuation-induced effects are governed by single-particle excitations: we find that, in this case, the confinement enhances the symmetric contribution and reduces the antisymmetric contribution.The antisymmetric contribution dominates the behaviour of the total friction and thermal conductance, but the overall confinement effect remains lower than 10%, except for the thermal conductance with graphene walls, where it reaches 50%.In general, the confinement effect is stronger for graphene walls than for jellium walls because the electronic fluctuations that mediate quantum friction and near-field heat transfer have a longer wavelength (smaller momentum q) in the case of graphene.
We now turn to the effect of the confinement-induced changes in the water fluctuations.Our simulations have shown that these changes become significant only at 7 Å confinement (Fig. 3), with an amplification of the Debye peak in the symmetric response and of the libration peak in the antisymmetric response.This translates into an enhancement of the friction coefficient and thermal conductance for both solid models, by up to a factor of 2 in the case of graphene.The study of fluctuationinduced effects specifically in 7 Å confinement is thus of particular interest.For instance, the thermal conductance of the interface between graphene and nanoconfined water may be probed with optical pump -terahertz probe spectroscopy: the optically-excited graphene electrons would be expected to cool faster than in the nonconfined case [20].

V. CONCLUSIONS
In this paper, we have developed a theoretical framework for studying confinement-and fluctuation-induced effects in two-dimensional nanofluidic channels: effects of (fluctuating) Coulomb interactions between the liquid and the solid.The key element of our framework is the description of the response of the solid walls and of the confined liquid to the Coulomb potentials that they apply to each other.In the case of a single solid-liquid interface, the surface response function -the reflection coefficient for evanescent plane waves -was found to be the most convenient descriptor.This convenience was due to the evanescent plane waves being eigenmodes of the Coulomb potential: the response to an evanescent wave is an evanescent wave.Generalizing this idea to the 2D channel geometry, we have introduced confined response functions, that play the role of reflection coefficients for the potential eigenmodes of the confined system.Our approach is systematic, and potentially extendable to more complex geometries.
The confined response functions reveal electrodynamic cross-talk between the walls of a 2D nanochannel.Investigating model materials that exhibit a surface plasmon mode, we found that the plasmons of the two walls couple as soon as the confinement is comparable to the plasmon wavelength.While the coupling of collective modes through Coulomb interactions is in principle a well-known phenomenon [31], our framework allows for the investigation of its effect on nanoscale fluid transport.From the fluid side, we have investigated confined water through molecular dynamics simulations.We found that the water confined response remains essentially bulk-like in the thermal frequency range down to 1.4 nm confinement, but undergoes significant changes when the confinement reaches 7 Å.
As an application of our framework, we have investigated quantum friction and near-field radiative heat transfer between water and the walls of a 2D nanochannel.We have generalized the derivation of refs.[19,20] to a the confined geometry and found that the confined response functions naturally emerge.In channels wider than 7 Å, the friction coefficient and thermal boundary conductance are modified compared to their bulk values when the confinement is comparable to the typical wavelength of the relevant charge fluctuations.At 7 Å confinement, a significant enhancement occurs for both effets, due to the drastic modification of the water response.Observing a confinement-induced modification of quantum friction or heat transfer thus appears most promising for systems where charge fluctuations are on longer wavelengths, yet these are also the systems where the effects are the weakest.Nevertheless, it has been shown that graphene-water heat transfer can be measured with ultrafast spectroscopy [20], and a small quantum friction coefficient can still result in a large quantum friction force (termed quantum feedback) if the electrons are driven by a phonon wind [17].Our results thus provide guidelines for engineering fluctuation-induced effects in nanoscale fluid transport.
Since the real part of the surface response functions consists of positive numbers smaller than one, Eq. ( 16) of the main text allows the antisymmetric response function to become larger than 1.This is in particular the case at zero frequency when the imaginary parts vanish.
Physically, it is easy to understand that the symmetric response function should be smaller than the surface response function while its antisymmetric counterpart should be larger.Indeed, let us focus on the top solid, which responds to the sum of the (positive) external potential and the potential induced by the bottom solid.In the symmetric case the bottom solid also responds to a positive external potential and then induces a negative potential which is damped exponentially.However, the surviving part reaching the top solid is still negative and then reduces the effective external potential to * Contact: nikita.kavokine@mpip-mainz.mpg.de.
which the top solid responds.The symmetric response is then smaller than the surface response because the external potential is screened by the bottom solid.On the contrary, in the antisymmetric case, the external potential is negative on the bottom solid which then responds by inducing a positive field.By the same process, the top solid now sees an effective external potential which is larger, generating a larger response.The antisymmetric response should then be larger than the surface response.
However, this cannot generate overscreening.Indeed, even if the antisymmetric response function is larger than 1, one should remember that the induced potential is not of the form of the external potential but of its "conjugate".The antisymmetric inner weight function is smaller than its outer counterpart: this solves the apparent problem.Let us formalise this intuition by computing the ratio of the induced field by the external field between the interfaces: where η = ±1 depending on whether we consider the symmetric or antisymmetric case.This ratio is maximal for z = ±h/2 where it reaches (2) Therefore, the induced potential is always smaller than the external potential, even when the confined response function is larger than 1.There is no overscreening.

II. ELECTROSTATIC INTERACTION CONFINEMENT
A. Model and solid's static response function In this section, we focus on how confinement will modify the electrostatic interactions between charges inside a arXiv:2306.00837v1[cond-mat.mes-hall] 1 Jun 2023 .This leads to a modified ion-ion interactions.b) Interaction confinement of an electrostatic dipole in the center of the channel: the generated electric field is partly confined in the channel by the surrounding solid which develop counter-charges (pictured in dark brown).This leads to a modified dipole-dipole interactions.c) Ion-ion interaction energy shift (normalised by the bulk interaction energy) for two ions at a distance ρ in bulk (black dashed line) and in a channel of confinement h = 2 nm.Different models of solid are used.d) Dipole-dipole interaction energy shift (normalised by the bulk interaction energy) for two dipoles in the z-direction at a distance ρ in bulk (black dashed line) and in a channel of confinement h = 2 nm.Different models of solid are used.

Single interface
Half-space channel.While our formalism allows us to treat the electric interactions of any static or dynamic distribution of charge, we restrict here to basic cases.We consider either an ion or a dipole in the center of the channel, as represented in Fig. 1a-b, to deduce the ion-ion and the dipole-dipole electrostatic interactions in presence of confinement.
To apply our theory and to compare with the existing literature, we need a simple model of the solid' surface response function.For this, we use the work of Kavokine et al. [2].To start, let us duplicate their eye-opening derivation of the surface response function of an insulator of static dielectric constant ϵ m ≈ 2, taking into account the presence of water of static dielectric constant ϵ w ≈ 80 at the interface.
We consider an external electric potential ϕ ext (q, z) = ϕ 0 ext (q)e qz applied on the solid in z < 0. On the solid area the total electric potential is denoted ϕ m (q, z) = ϕ 0 m (q)e qz and the solid generates an induced electric potential ϕ ind (q, z) = ϕ 0 m (q)e −qz on the area z > 0. The boundary conditions at z = 0 are the continuity of both the total electric potential ϕ and the electric displacement field D = −ϵ∇ϕ, leading to Finally, we obtain the surface response function: This derivation is valid for an insulator, while the dielectric constant of the solid depends on the momentum q in general.A generalisation of this formula exists for a Thomas-Fermi model of solid, parametrised with a Thomas-Fermi wavelength q TF [2]: (6) The case q TF = 0 corresponds to an insulator.The case q TF → ∞, that is g 0 o (q) → 1, corresponds to a metal.Let us highlight that g 0 m is negative at large wavelength and positive at small wavelength.The wavelength of the transition depends on the Thomas-Fermi wavelength.When the surface response function is negative, which is expected to happen at small distances, the effect of the boundary solids is to confine the electric field inside the channel: this is interaction confinement.At larger distances, the surface response function becomes ultimately positive: the electric field is no longer confined.

B. Interaction confinement for an ion and ion-ion electrostatic interaction
Let us start with the case of a static ion of charge e (see Fig. 1a).Such a charge generates a potential where we have introduced a Fourier decomposition into independent modes.Taking into account the effect of the material, the total potential in the inner space now writes In particular, in the center of the channel, that is at z = 0, we obtain using Eq. ( 17) of the main text: consistently with the result of Kavokine et al. [2].An ion in the center of the channel is a symmetric distribution of charge and then its potential is corrected by the symmetric response function of the outer medium.The effect of the medium is to screen the charge and then to confine the electric potential, as pictured in Fig. 1a.Considering another ion of charge e in the center of the channel at a distance ρ, the ion-ion electrostatic interaction is ) where we have used isotropy and introduced the Bessel function of the first kind J 0 .The first term provides te bulk interaction energy while the second term provides a correction due to confinement.This correction can be computed numerically using Eq. ( 6).The resulting ion-ion interaction (normalised by E bulk ii (ρ)) is plotted in Fig. 1c.Consistently with [2], we find that the interaction is reduced at large distance and increased at small distance for a generic solid.

C. Interaction confinement for a dipole and dipole-dipole electrostatic interaction
Let us now turn to an electrostatic dipole d oriented in the z direction (see Fig. 1b).Such a dipole generates a potential where we have introduced a Fourier decomposition into independent modes.Taking into account the effect of the material, the total potential in the inner space now writes 13) In particular, in the center of the channel, using Eq. ( 17) of the main text, the electric field writes An electrostatic dipole in the center of the channel is an antisymmetric distribution of charge and then its potential is corrected by the antisymmetric response function of the outer medium.Here again, the effect of the outer medium is to confine the electric potential, as pictured in Fig. 1b.Considering another dipole of same moment d in the center of the channel at a distance ρ, the dipole-dipole electrostatic interaction is The first term provides te bulk interaction energy while the second term provides a correction due to confinement.This correction can be computed numerically using Eq. ( 6).The resulting dipole-dipole interaction (normalised by E bulk dd (ρ)) is plotted in Fig. 1d.After a reduction of the interaction at small distances compared with the confinement h, the interaction energy is increased by a factor ≈ 2.6 independent on the model of solid.
To conclude this section, we have derived the correction of confinement to the electric potential.This approach generalises the results of interaction confinement [2] to the generic distribution of charge, or equivalently to a generic external potential, which we can decompose into a symmetric and an antisymmetric parts.

III. WATER RESPONSE FUNCTIONS AND SIMULATIONS
A. Numerical methods for water simulations We have carried out both force-field (FF) and ab initio density functional theory-based (DFT) molecular dynamics (MD) simulations of water confined between two frozen graphene sheets for various separations h between the graphene sheets.
FF-MD simulations are performed using GROMACS 2020 using the SPC/E water model and GROMOS 53A6 parameters for carbon atoms [3,4].Simulations are performed in the NVT ensemble using the CSVR thermostat at 300 K. We use a 2 fs time step and truncate Lennard-Jones interactions at 9 Å.Electrostatics employ a realspace cut-off at 0.9 nm while long-range interactions are handled with the particle mesh Ewald method.The considered systems have lateral dimensions of 42.6 Å × 44.3 Å and contain 210, 1872 and 3489 water molecules at h =7 Å, 34 Å and 60 Ågraphene sheet separations, respectively.After 1 ns of equilibration, we perform production runs of 20 ns.
DFT-MD simulations are carried out with the CP2K 6.1 software which performs Born-Oppenheimer MD [5].Electronic structures are calculated at every time step based on the BLYP exchange correlation functional with Grimme-D3 dispersion correction [6], while core electrons are represented by GTH pseudo potentials and valence electrons are expanded in the DZVP-SR-MOLOPT basis set using a plane wave cutoff of 400 Ry [7,8].DFT-MD simulations are performed in the NVT ensemble at 300K with a 0.5 fs time step using the CSVR thermostat.The simulated systems have lateral dimensions of 29.82 Å × 27.05 Å and contain 90 and 420 water molecules for graphene sheet separations of h =7 and 17.8 Å, respectively.After preequilibration in FF-MD, DFT-simulations are equilibrated for another 5 ps followed by production runs of 70 ps.

C. Evaluation of the time correlation function
For each time step, we first compute the Fouriertransformed charge density.From MD simulation, we have access to the center of each atom and use a partial charge Z = δ = 0.41e for the hydrogen atoms and Z = −2δ for the oxygen atoms.The density is then n(q) = j Z j F (q, z j )e iqxx+iqyy (21) where F is the weight function and q x , q y are multiples of 2π/L x (2π/L y respectively).We make this computation for N t = 10 5 time steps, and then obtain the fluctuating part of the density by subtracting the average δn(q, t) = n(q, t) − ⟨n(q, t)⟩ t (22) where ⟨A(t)⟩ t = 1 Nt t A(t).
We then carry out a fast Fourier transform (FFT): δn(q, ω) = FFT[n(q, t)] t (23) and only keep the positive frequencies.The response function is then computed as Im [g(q, ω)] = e 2 ω 4ϵ 0 qT A |δn(q, ω) where dt is the length of the time step.In practice, g only depends on the norm q.Finally, the spectra are smoothed by a Gaussian smoothing of width 0.2 meV for the FF simulations and 0.5 meV for DFT simulations.

D. Momentum dependence and position of the interface
The limited size of the simulation box does not allow us to access all momenta necessary for the computation of fluctuation-induced effects (typically below 1 Å −1 ).Therefore, we restrict ourselves to one momentum compatible with the size of the box and extrapolate the response functions to lower momenta.We use q 0 = 0.67 Å −1 for FF simulations and q 0 = 1 Å −1 for DFT simulations.
From the simulations carried out with a large box in [1], we know that the non-confined water surface response function has only a weak momentum dependence: we therefore approximate the surface response function at all momenta by the one computed at q 0 .For the confined response functions, however, the form of the weight functions imply that the static (ω = 0) response fulfills g x w ≤ 1 + ηe −qh , where η = ±1 for the symmetric (antisymmetric) component (see.ESI.1).In particular, the static antisymmetric response function goes to zero at q → 0. We thus infer the momentum dependence of the confined response function according to g x w (q, ω) = 1 + ηe −qh 1 + ηe −q0h g x w (q 0 , ω).
In practice, in the evaluation of friction coefficients, this procedure makes only a 1% difference compared to the assumption of no momentum dependence.Formally, however, such a procedure allows us to obtain wellbehaved response functions that produce no unphysical overscreening.
The evaluation of the response functions requires to precisely define the positions of the interfaces at z = ±h/2.In ref. [1], a procedure for positioning the interface based on imposing the compressibility sum rule for the surface response function was proposed.This procedure no longer formally holds in a confined geometry, but we expect the position of the interface to not be significantly affected by confinement.With our simulations at weakest confinement (FF at 60 Å and DFT at 34 Å) and following the procedure of ref. [1], we an effective distance between the plane of the carbon atoms and the water surface d ≈ 1.4 Å, which is compatible with d = 1.3 Åfound in [1].We then adopt the prescription d = 1.4 Å for the simulations under confinement.

FIG. 1 .
FIG. 1. Electric response of interfacial systems.a) Schematic of the confined geometry under consideration, emphasizing the role of the channel walls' electronic degrees of freedom.b-f): Electric response in different geometries.b) Single-interface case: the response function is the usual surface response function.c-f) In confined geometry there are four situations to distinguish.The responding medium can be either in the outer space (c&d) or in the inner space (e&f).The external potential can be either symmetric (c&e) or antisymmetric (d&f).We observe that the inter-solid interactions bring corrections by lowering the induced potential in the symmetric case (c) and increasing the induced potential in the antisymmetric case (d).

FIG. 2 .
FIG. 2. Surface and confined response functions for different solids.a-f) Interfacial response functions Im[g(q, ω)] of the solid as a function of momentum q and frequency ω.The first row (a-c) is for graphene with Fermi level µ = 100 meV.The second row (d-f) is for a jellium model with an effective mass m = 0.1 me and a Fermi level µ = 180 meV (electron density parameter rs = 5), corresponding to a doped semi-conductor.The first column (a&d) corresponds to the surface response functions.The second column (b&e) corresponds to the symmetric confined response functions with a confinement of 7 Å.The third column (c&f) corresponds to the anti-symmetric confined response functions with a confinement of 7 Å.

FIG. 3 .
FIG.3.Water response functions from simulation.a) Snapshot of a FF molecular dynamics simulation with 7 Å confinement b-d) Response functions Im[g(q, ω)] of water as a function of the frequency ω obtained from the simulations, at fixed wavevector q = 0.67 Å −1 for FF and q = 1 Å −1 for DFT.Different confinements are used: 7 Å, 14 Å (only FF), 18 Å (only DFT), 34 Å and 60 Å (only FF).b) Surface response function computed from simulations at weak confinement (FF at h =60 Å and DFT at h =34 Å). c) Symmetric response function for different confinements.The dashed black line is the non-confined surface response function.d) Antisymmetric response function for different confinements.The dashed black line is the non-confined surface response function.e) Schematic of a channel containing a single water monolayer on which a symmetric potential is applied.f) Same as (e) with application of an antisymmetric potential.

FIG. 4 .
FIG. 4. Fluctuation-induced effects.(a) Schematic of the interactions involved in single-interface quantum friction.Friction results from the Coulomb coupling of liquid fluctuations (hydrons) and electronic fluctuations.(b) Schematic of the interactions involved in confined quantum friction.Compared to the non-confined case, there is an additional Coulomb interaction between the charge fluctuations in the two solid walls.(c) Dyson equation for the renormalization of the bottom wall's charge density response function (polarization "bubble") by the intra-wall Coulomb interactions.The bare bubble is empty and the renormalized bubble is filled with gray.(d) Dyson equation for the electron-water susceptibility χew.(e) Definition of the exchange term Π (see text).(f ) Dyson equation for the renormalization of the bottom wall's charge density response function by the inter-wall Coulomb interactions.The stripes indicate a fully-renormalized bubble.(g) Dyson equation for the inter-wall charge density response function, which vanishes in the absence of inter-wall Coulomb interactions (Π = 0).

H
h−e (t) = dr w dr e n w (r w , t)V (r w − r e )n e (r e , t).(21)The dynamics of the systems are governed by the interaction Hamiltonian that also comprises the electronelectron Coulomb interactions: H int = H e−e + H h−e .The friction force and heat transfer rate are given by ⟨F⟩ = − dr w dr e ∇V (r w − r e )⟨n w (r w − vt, t)n e (r e , t)⟩,

FIG. 5 .
FIG.5.Effect of confinement on the fluctuation-induced phenomena.a) Quantum friction coefficient λ normalised by the single interface quantum friction coefficient λ 0 as a function of confinement between graphene walls.b) Thermal boundary conductance κ normalised by the single interface thermal boundary conductance κ 0 as a function of confinement between graphene walls.c) Quantum friction coefficient λ normalised by the single interface quantum friction coefficient λ 0 as a function of confinement between jellium walls with an effective mass m = 0.1 me and a Fermi level µ = 180 meV (electron density parameter rs = 5) d) Thermal boundary conductance κ normalised by the single interface thermal boundary conductance κ 0 as a function of confinement between jellium walls with an effective mass m = 0.1 me and a Fermi level µ = 180 meV (electron density parameter rs = 5).In all panels, the continuous lines are obtained using the single interface surface response function for water, while the crosses are obtained with the confined response function of water at the relevant confinement, as obtained from molecular dynamics simulations.

FIG. 1 .
FIG.1.a) Interaction confinement of an ion in the center of the channel: the generated electric field is partly confined in the channel by the surrounding solid which develop counter-charges (pictured in dark brown).This leads to a modified ion-ion interactions.b) Interaction confinement of an electrostatic dipole in the center of the channel: the generated electric field is partly confined in the channel by the surrounding solid which develop counter-charges (pictured in dark brown).This leads to a modified dipole-dipole interactions.c) Ion-ion interaction energy shift (normalised by the bulk interaction energy) for two ions at a distance ρ in bulk (black dashed line) and in a channel of confinement h = 2 nm.Different models of solid are used.d) Dipole-dipole interaction energy shift (normalised by the bulk interaction energy) for two dipoles in the z-direction at a distance ρ in bulk (black dashed line) and in a channel of confinement h = 2 nm.Different models of solid are used.

TABLE I .
Weight functions used in the definitions of the surface and confined response functions (Eq.(??) and (??)).