# Casimir-Lifshitz force for graphene-covered gratings

Youssef Jeyar, Minggang Luo, Brahim Guizal, Ho Bun Chan, Mauro Antezza

► To cite this version:

Youssef Jeyar, Minggang Luo, Brahim Guizal, Ho Bun Chan, Mauro Antezza. Casimir-Lifshitz force for graphene-covered gratings. 2024. hal-04603644

HAL Id: hal-04603644

https://hal.science/hal-04603644

Preprint submitted on 19 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Casimir-Lifshitz force for graphene-covered gratings

Youssef Jeyar, Minggang Luo, Brahim Guizal, Ho Bun Chan, Mauro Antezza

► To cite this version:

Youssef Jeyar, Minggang Luo, Brahim Guizal, Ho Bun Chan, Mauro Antezza. Casimir-Lifshitz force for graphene-covered gratings. Physical Review B, 2024, 110 (4), pp.045446. 10.1103/Phys-RevB.110.045446. hal-04603644

HAL Id: hal-04603644

https://hal.science/hal-04603644v1

Submitted on 19 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Casimir-Lifshitz force for graphene-covered gratings

Youssef Jeyar,¹, * Minggang Luo,¹ Brahim Guizal,¹ H. B. Chan,², ³, ⁴ and Mauro Antezza¹, ⁵, †

¹Laboratoire Charles Coulomb (L2C), UMR 5221 CNRS-Université de Montpellier, F-34095 Montpellier, France

²Department of Physics, The Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong, China

³William Mong Institute of Nano Science and Technology,
The Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong, China

⁴Center for Metamaterial Research, The Hong Kong University of
Science and Technology, Clear Water Bay, Kowloon, Hong Kong, China

⁵Institut Universitaire de France, 1 rue Descartes, Paris Cedex 05, F-75231, France

(Dated: May 24, 2024)

We study the Casimir-Lifshitz force (CLF) between a gold plate and a graphene-covered dielectric grating. Using a scattering matrix (S-matrix) approach derived from the Fourier Modal Method (FMM), we find a significant enhancement in the CLF as compared to a mere dielectric slab coated with graphene, over a wide range of temperatures. Additionally, we demonstrate that the CLF depends strongly on the chemical potential of graphene, with maximal effects observed at lower filling fractions. Finally, we analyse the Casimir force gradient between a gold sphere and a graphene-coated dielectric grating, highlighting potential avenues for experimental measurements.

# I. INTRODUCTION

The Casimir-Lifshitz force (CLF), resulting from quantum vacuum fluctuations, manifests in the interactions between polarizable bodies. In 1948, Casimir [1] first developed a theoretical framework for this force between two perfectly conducting plates at zero temperature. Later, Lifshitz expanded this framework to encompass bodies with diverse optical properties and at finite temperatures [2]. Since its experimental verification [3], the CLF has attracted significant research interest, with investigations focused on exploring different geometries and materials [4–11].

Among diverse geometrical systems, the grating configuration is particularly important in the study of the CLF and has been explored both theoretically [12, 13] and experimentally [14, 15], in terms of the non-trivial interplay between geometry and dielectric properties of the interacting bodies. Analyzing such structures necessitates precise and efficient methods such as the Fourier Modal Method (FMM) or its improved version, the FMM with adaptive spatial resolution (FMM-ASR), which significantly enhances efficiency for metallic gratings [16].

The dielectric properties of the interacting objects can affect the scattering details, which are relevant to both CLF in and out of thermal equilibrium and radiative heat transfer. In particular, there has been much progress in graphene-based structures [17–28], that utilizes the unique optical properties of this material [29–33]. A recent study [34] delved into the CLF between graphene gratings, utilizing another enhanced version of the FMM approach known as FMM with local basis functions [35]. This study shows the potential for modulating the force by manipulating geometry parameters and/or graphene

properties. The investigation is further extended to analyzing heat transfer within these structures for both aligned [36] and misaligned [37] cases. Patterning the graphene sheet into grating will result in a topology transition of the modes for energy transfer from the circular one to the hyperbolic one, which allows a significant enhancement in the heat transfer.

In another example [38], graphene is integrated with one-dimensional plasmonic nanogratings with high aspect ratios and narrow grooves, the good agreement between theory and experiment suggests promising prospects for high-performance electrically tunable graphene-based infrared photodetectors and biological sensors.

In this work, we explore the CLF between graphene coated dielectric grating and gold semi-finite plane, using the scattering matrix approach obtained by an adapted FMM to include the graphene sheet. This configuration combines the features of both the grating structures and the special dielectric properties of graphene, which is expected to provide more freedom in controlling the CLF.

In section II, we describe the physical system, in section III we introduce the reflection operators obtained by the FMM. Finally, in section IV, we present and discuss the numerical results. The calculation details of the scattering matrix are provided in the appendix.

# II. PHYSICAL SYSTEM

The physical system depicted in Fig. 1 consists of a dielectric grating coated with a graphene sheet (body 1) facing a gold half-space (body 2). These two bodies are separated by a distance *d* and maintained at a temperature *T*. The dielectric grating, composed of fused silica (SiO₂), exhibits periodicity in the *x* direction with a period *D* (we use a fixed value of *D* = 1µm in this study), each period contains an SiO₂ region of width *a*, giving a filling fraction *f* = *a*/*D*. The grating has a thickness *h*₁

<footnote>* youssef.jeyar@umontpellier.fr</footnote>

<footnote>† mauro.antezza@umontpellier.fr</footnote>
Casimir-Lifshitz force for graphene-covered gratings

Youssef Jeyar,¹, * Minggang Luo,¹ Brahim Guizal,¹ H. B. Chan,², ³, ⁴ and Mauro Antezza¹, ⁵, †

¹Laboratoire Charles Coulomb (L2C), UMR 5221 CNRS-Université de Montpellier, F-34095 Montpellier, France

²Department of Physics, The Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong, China

³William Mong Institute of Nano Science and Technology,
The Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong, China

⁴Center for Metamaterial Research, The Hong Kong University of
Science and Technology, Clear Water Bay, Kowloon, Hong Kong, China

⁵Institut Universitaire de France, 1 rue Descartes, Paris Cedex 05, F-75231, France

(Dated: May 24, 2024)

We study the Casimir-Lifshitz force (CLF) between a gold plate and a graphene-covered dielectric grating. Using a scattering matrix (S-matrix) approach derived from the Fourier Modal Method (FMM), we find a significant enhancement in the CLF as compared to a mere dielectric slab coated with graphene, over a wide range of temperatures. Additionally, we demonstrate that the CLF depends strongly on the chemical potential of graphene, with maximal effects observed at lower filling fractions. Finally, we analyse the Casimir force gradient between a gold sphere and a graphene-coated dielectric grating, highlighting potential avenues for experimental measurements.

# I. INTRODUCTION

The Casimir-Lifshitz force (CLF), resulting from quantum vacuum fluctuations, manifests in the interactions between polarizable bodies. In 1948, Casimir [1] first developed a theoretical framework for this force between two perfectly conducting plates at zero temperature. Later, Lifshitz expanded this framework to encompass bodies with diverse optical properties and at finite temperatures [2]. Since its experimental verification [3], the CLF has attracted significant research interest, with investigations focused on exploring different geometries and materials [4–11].

Among diverse geometrical systems, the grating configuration is particularly important in the study of the CLF and has been explored both theoretically [12, 13] and experimentally [14, 15], in terms of the non-trivial interplay between geometry and dielectric properties of the interacting bodies. Analyzing such structures necessitates precise and efficient methods such as the Fourier Modal Method (FMM) or its improved version, the FMM with adaptive spatial resolution (FMM-ASR), which significantly enhances efficiency for metallic gratings [16].

The dielectric properties of the interacting objects can affect the scattering details, which are relevant to both CLF in and out of thermal equilibrium and radiative heat transfer. In particular, there has been much progress in graphene-based structures [17–28], that utilizes the unique optical properties of this material [29–33]. A recent study [34] delved into the CLF between graphene gratings, utilizing another enhanced version of the FMM approach known as FMM with local basis functions [35]. This study shows the potential for modulating the force by manipulating geometry parameters and/or graphene

properties. The investigation is further extended to analyzing heat transfer within these structures for both aligned [36] and misaligned [37] cases. Patterning the graphene sheet into grating will result in a topology transition of the modes for energy transfer from the circular one to the hyperbolic one, which allows a significant enhancement in the heat transfer.

In another example [38], graphene is integrated with one-dimensional plasmonic nanogratings with high aspect ratios and narrow grooves, the good agreement between theory and experiment suggests promising prospects for high-performance electrically tunable graphene-based infrared photodetectors and biological sensors.

In this work, we explore the CLF between graphene coated dielectric grating and gold semi-finite plane, using the scattering matrix approach obtained by an adapted FMM to include the graphene sheet. This configuration combines the features of both the grating structures and the special dielectric properties of graphene, which is expected to provide more freedom in controlling the CLF.

In section II, we describe the physical system, in section III we introduce the reflection operators obtained by the FMM. Finally, in section IV, we present and discuss the numerical results. The calculation details of the scattering matrix are provided in the appendix.

# II. PHYSICAL SYSTEM

The physical system depicted in Fig. 1 consists of a dielectric grating coated with a graphene sheet (body 1) facing a gold half-space (body 2). These two bodies are separated by a distance *d* and maintained at a temperature *T*. The dielectric grating, composed of fused silica (SiO₂), exhibits periodicity in the *x* direction with a period *D* (we use a fixed value of *D* = 1µm in this study), each period contains an SiO₂ region of width *a*, giving a filling fraction *f* = *a*/*D*. The grating has a thickness *h*₁

<footnote>* youssef.jeyar@umontpellier.fr</footnote>

<footnote>† mauro.antezza@umontpellier.fr</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. Schematic of the system: a dielectric grating covered with graphene sheet and is placed at a distance d from a gold half-space. The grating is periodic in the x direction with period D, characterized by width a, a depth h₁ and filling fraction f = a/D.</figcaption></figure>

of 300 nm and is positioned on top of an SiO₂ slab with a thickness h₂ of 100 nm, which in turn sits on an silicon (Si) half-space.

For this configuration, the Casimir-Lifshitz pressure (CLP) acting on body 1 along the z positive direction can be expressed as [13, 34, 39]

$$P(d, T, \mu) = \frac{k_B T}{4\pi^2} \sum_{m=0}^{+\infty}{}' \int_{-\frac{\pi}{D}}^{\frac{\pi}{D}} dk_x \int_{-\infty}^{+\infty} dk_y \mathrm{Tr}(\gamma' \mathcal{M}), \quad (1)$$

where $\gamma' = \mathrm{diag}(k'_{zn}, k'_{zn})_n$, $k'_{zn} = \sqrt{\xi_m^2/c^2 + \mathbf{k}_n^2}$, $\mathbf{k}_n = (k_{xn}, k_y)$, $k_{xn} = k_x + n\frac{2\pi}{D}$, $k_x$ is in the first Brillouin zone $[-\frac{\pi}{D}, \frac{\pi}{D}]$ and $k_y$ is in $\mathbb{R}$. The notation $\mathrm{diag}(u_n, u_n)_n$ represents a diagonal matrix of size $2(2N+1) \times 2(2N+1)$, where the diagonal elements are $u_{-N}, u_{-N+1}, \dots, u_N, u_{-N}, u_{-N+1}, \dots, u_N$. The expression for $\mathcal{M}$ is given by

$$\mathcal{M} = (U^{(12)}\mathcal{R}^{(1)+}\mathcal{R}^{(2)-} + U^{(21)}\mathcal{R}^{(2)-}\mathcal{R}^{(1)+}). \quad (2)$$

where $\mathcal{R}^{(1)+}$ and $\mathcal{R}^{(2)-}$ represent the reflection operators corresponding to the two bodies in the (TE, TM) basis obtained using the scattering matrix approach, as will be discussed hereafter. As for the inter-cavity operators, they are defined as follows ($\mathcal{I}$ denoting the operator identity):

$$\begin{aligned} U^{(12)} &= (\mathcal{I} - \mathcal{R}^{(1)+}\mathcal{R}^{(2)-})^{-1}, \\ U^{(21)} &= (\mathcal{I} - \mathcal{R}^{(2)-}\mathcal{R}^{(1)+})^{-1}. \end{aligned} \quad (3)$$

In Eq. (1), the sum is over the Matsubara frequencies ($\xi_m = 2\pi mk_B T/\hbar$) and the prime on this sum means that

the first term ($m=0$) is divided by 2. $k_B$ is the Boltzmann constant, and $\hbar$ is the reduced Planck's constant.

For fused silica, the dielectric response at imaginary frequencies, is obtained through the Kramers-Kronig relation: $\varepsilon(i\xi_m) = 1 + 2\pi^{-1} \int_0^\infty \omega\varepsilon''(\omega)/(\omega^2 + \xi_m^2)d\omega$ and relies on the data extracted from the dielectric function of SiO₂, as sourced from [40], extrapolated along the real frequency axis.

We adopt the following model for the dielectric function of silicon [15]:

$$\varepsilon_{\mathrm{Si}}(i\xi_m) = 1.035 + \frac{11.87 - 1.035}{1 + \xi_m^2/\omega_0^2} + \frac{\omega_p^2}{\xi_m(\xi_m + \Gamma)} \quad (4)$$

where $\omega_0 = 6.6 \times 10^{15}$ rad s⁻¹, $\omega_p = 2.37 \times 10^{14}$ rad s⁻¹, and $\Gamma = 6.45 \times 10^{13}$ rad s⁻¹.

We employ the Drude permittivity model for the dielectric function of gold (Au)

$$\varepsilon_{\mathrm{Au}}(i\xi_m) = 1 + \frac{\omega_p'^2}{\xi_m(\xi_m + \gamma)} \quad (5)$$

with $\hbar\omega_p' = 9$ eV and $\hbar\gamma = 35$ meV.

The contribution of graphene is introduced through its conductivity, which varies with temperature T and chemical potential $\mu$. We model the conductivity as the sum of intraband and interband components: $\sigma_g = \sigma_{\mathrm{intra}} + \sigma_{\mathrm{inter}}$. More precisely, their expressions on the imaginary frequency axis are given by [26,29,30,31]

$$\begin{aligned} \sigma_{\mathrm{intra}}(i\xi_m) &= \frac{8\sigma_0 k_B T}{\pi(\hbar\xi_m + \hbar/\tau)} \ln \left[ 2 \cosh \left( \frac{\mu}{2k_B T} \right) \right], \\ \sigma_{\mathrm{inter}}(i\xi_m) &= \frac{\sigma_0 4\hbar\xi_m}{\pi} \int_0^{+\infty} \frac{G(x)}{(\hbar\xi_m)^2 + 4x^2} dx, \end{aligned} \quad (6)$$

where, $\sigma_0 = e^2/(4\hbar)$, e is the electron charge, $G(x) = \sinh(x/k_B T)/[\cosh(\mu/k_B T) + \cosh(x/k_B T)]$, $\tau$ the relaxation time (we use $\tau = 10^{-13}$s).

### III. REFLECTION OPERATORS

In this section, we present the derivation of the reflection operators for the interacting bodies. We begin with body 1, which comprises a dielectric grating covered with graphene, as depicted in Fig. 1. The reflection operator can be directly obtained, as sub-blocs, from the scattering matrix, which has dimensions of $4(2N+1) \times 4(2N+1)$ and can represented as follows

$$S = \begin{pmatrix} \mathcal{R}^{-} & \mathcal{T}^{-} \\ \mathcal{T}^{+} & \mathcal{R}^{+} \end{pmatrix}. \quad (7)$$

This matrix is obtained through the use of the FMM (see detailed calculations in the appendix) and takes the form:

$$S = \begin{pmatrix} 1 & 0 \\ 0 & \Phi \end{pmatrix} \tilde{S} \begin{pmatrix} 1 & 0 \\ 0 & \Phi \end{pmatrix}, \quad (8)$$
<figure><img src="image_2.png" /><figcaption>FIG. 1. Schematic of the system: a dielectric grating covered with graphene sheet and is placed at a distance d from a gold half-space. The grating is periodic in the x direction with period D, characterized by width a, a depth h₁ and filling fraction f = a/D.</figcaption></figure>

of 300 nm and is positioned on top of an SiO₂ slab with a thickness h₂ of 100 nm, which in turn sits on an silicon (Si) half-space.

For this configuration, the Casimir-Lifshitz pressure (CLP) acting on body 1 along the z positive direction can be expressed as [13, 34, 39]

$$P(d, T, \mu) = \frac{k_B T}{4\pi^2} \sum_{m=0}^{+\infty}{}' \int_{-\frac{\pi}{D}}^{\frac{\pi}{D}} dk_x \int_{-\infty}^{+\infty} dk_y \mathrm{Tr}(\gamma' \mathcal{M}), \quad (1)$$

where $\gamma' = \mathrm{diag}(k'_{zn}, k'_{zn})_n$, $k'_{zn} = \sqrt{\xi_m^2/c^2 + \mathbf{k}_n^2}$, $\mathbf{k}_n = (k_{xn}, k_y)$, $k_{xn} = k_x + n\frac{2\pi}{D}$, $k_x$ is in the first Brillouin zone $[-\frac{\pi}{D}, \frac{\pi}{D}]$ and $k_y$ is in $\mathbb{R}$. The notation $\mathrm{diag}(u_n, u_n)_n$ represents a diagonal matrix of size $2(2N+1) \times 2(2N+1)$, where the diagonal elements are $u_{-N}, u_{-N+1}, \dots, u_N, u_{-N}, u_{-N+1}, \dots, u_N$. The expression for $\mathcal{M}$ is given by

$$\mathcal{M} = (U^{(12)}\mathcal{R}^{(1)+}\mathcal{R}^{(2)-} + U^{(21)}\mathcal{R}^{(2)-}\mathcal{R}^{(1)+}). \quad (2)$$

where $\mathcal{R}^{(1)+}$ and $\mathcal{R}^{(2)-}$ represent the reflection operators corresponding to the two bodies in the (TE, TM) basis obtained using the scattering matrix approach, as will be discussed hereafter. As for the inter-cavity operators, they are defined as follows ($\mathcal{I}$ denoting the operator identity):

$$\begin{aligned} U^{(12)} &= (\mathcal{I} - \mathcal{R}^{(1)+}\mathcal{R}^{(2)-})^{-1}, \\ U^{(21)} &= (\mathcal{I} - \mathcal{R}^{(2)-}\mathcal{R}^{(1)+})^{-1}. \end{aligned} \quad (3)$$

In Eq. (1), the sum is over the Matsubara frequencies ($\xi_m = 2\pi mk_B T/\hbar$) and the prime on this sum means that

the first term ($m=0$) is divided by 2. $k_B$ is the Boltzmann constant, and $\hbar$ is the reduced Planck's constant.

For fused silica, the dielectric response at imaginary frequencies, is obtained through the Kramers-Kronig relation: $\varepsilon(i\xi_m) = 1 + 2\pi^{-1} \int_0^\infty \omega\varepsilon''(\omega)/(\omega^2 + \xi_m^2)d\omega$ and relies on the data extracted from the dielectric function of SiO₂, as sourced from [40], extrapolated along the real frequency axis.

We adopt the following model for the dielectric function of silicon [15]:

$$\varepsilon_{\mathrm{Si}}(i\xi_m) = 1.035 + \frac{11.87 - 1.035}{1 + \xi_m^2/\omega_0^2} + \frac{\omega_p^2}{\xi_m(\xi_m + \Gamma)} \quad (4)$$

where $\omega_0 = 6.6 \times 10^{15}$ rad s⁻¹, $\omega_p = 2.37 \times 10^{14}$ rad s⁻¹, and $\Gamma = 6.45 \times 10^{13}$ rad s⁻¹.

We employ the Drude permittivity model for the dielectric function of gold (Au)

$$\varepsilon_{\mathrm{Au}}(i\xi_m) = 1 + \frac{\omega_p'^2}{\xi_m(\xi_m + \gamma)} \quad (5)$$

with $\hbar\omega_p' = 9$ eV and $\hbar\gamma = 35$ meV.

The contribution of graphene is introduced through its conductivity, which varies with temperature T and chemical potential $\mu$. We model the conductivity as the sum of intraband and interband components: $\sigma_g = \sigma_{\mathrm{intra}} + \sigma_{\mathrm{inter}}$. More precisely, their expressions on the imaginary frequency axis are given by [26,29,30,31]

$$\begin{aligned} \sigma_{\mathrm{intra}}(i\xi_m) &= \frac{8\sigma_0 k_B T}{\pi(\hbar\xi_m + \hbar/\tau)} \ln \left[ 2 \cosh \left( \frac{\mu}{2k_B T} \right) \right], \\ \sigma_{\mathrm{inter}}(i\xi_m) &= \frac{\sigma_0 4\hbar\xi_m}{\pi} \int_0^{+\infty} \frac{G(x)}{(\hbar\xi_m)^2 + 4x^2} dx, \end{aligned} \quad (6)$$

where, $\sigma_0 = e^2/(4\hbar)$, e is the electron charge, $G(x) = \sinh(x/k_B T)/[\cosh(\mu/k_B T) + \cosh(x/k_B T)]$, $\tau$ the relaxation time (we use $\tau = 10^{-13}$s).

### III. REFLECTION OPERATORS

In this section, we present the derivation of the reflection operators for the interacting bodies. We begin with body 1, which comprises a dielectric grating covered with graphene, as depicted in Fig. 1. The reflection operator can be directly obtained, as sub-blocs, from the scattering matrix, which has dimensions of $4(2N+1) \times 4(2N+1)$ and can represented as follows

$$S = \begin{pmatrix} \mathcal{R}^{-} & \mathcal{T}^{-} \\ \mathcal{T}^{+} & \mathcal{R}^{+} \end{pmatrix}. \quad (7)$$

This matrix is obtained through the use of the FMM (see detailed calculations in the appendix) and takes the form:

$$S = \begin{pmatrix} 1 & 0 \\ 0 & \Phi \end{pmatrix} \tilde{S} \begin{pmatrix} 1 & 0 \\ 0 & \Phi \end{pmatrix}, \quad (8)$$
where

$$ \Phi = \left( \begin{array}{cc} \mathrm{diag} \left( e^{ik_{zn}^{(4)}(h_1+h_2)} \right)_n & \mathbb{0} \\ \mathbb{0} & \mathrm{diag} \left( e^{ik_{zn}^{(4)}(h_1+h_2)} \right)_n \end{array} \right), \quad (9) $$

$k_{zn}^{(i)} = \sqrt{\varepsilon_i(\omega)\frac{\omega^2}{c^2} - k_n^2}$ for medium i (see details in the appendix), and

$$ \tilde{S} = \tilde{S}_1 \circledast \tilde{S}_2 \circledast \tilde{S}_3. \quad (10) $$

Here, the star product operation $\mathcal{A} = \mathbb{B} \circledast \mathbb{C}$ is defined as [16]

$$ \begin{aligned} \mathcal{A}_{11} &= \mathbb{B}_{11} + \mathbb{B}_{12}(1 - \mathbb{C}_{11}\mathbb{B}_{22})^{-1}\mathbb{C}_{11}\mathbb{B}_{21}, \\ \mathcal{A}_{12} &= \mathbb{B}_{12}(1 - \mathbb{C}_{11}\mathbb{B}_{22})^{-1}\mathbb{C}_{12}, \\ \mathcal{A}_{21} &= \mathbb{C}_{21}(1 - \mathbb{B}_{22}\mathbb{C}_{11})^{-1}\mathbb{B}_{21}, \\ \mathcal{A}_{22} &= \mathbb{C}_{22} + \mathbb{C}_{21}(1 - \mathbb{B}_{22}\mathbb{C}_{11})^{-1}\mathbb{B}_{22}\mathbb{C}_{12}, \end{aligned} \quad (11) $$

and

$$ \begin{aligned} \tilde{S}_1 &= \left( \begin{array}{cc} \mathbb{K}'_1 & -\mathbb{P} \\ \mathbb{L}'_1 + \mathbb{M}'_1 & -\mathbb{P}' \end{array} \right)^{-1} \left( \begin{array}{cc} \mathbb{K}_1 & \mathbb{P} \\ \mathbb{L}_1 + \mathbb{M}_1 & -\mathbb{P}' \end{array} \right), \\ \tilde{S}_2 &= \left( \begin{array}{cc} \Phi_1 & \mathbb{0} \\ \mathbb{0} & 1 \end{array} \right) \left( \begin{array}{cc} -\mathbb{P} & -\mathbb{K}_3 \\ \mathbb{P}' & -\mathbb{L}_3 \end{array} \right)^{-1} \left( \begin{array}{cc} \mathbb{P} & -\mathbb{K}'_3 \\ \mathbb{P}' & -\mathbb{L}_3 \end{array} \right) \left( \begin{array}{cc} \Phi_1 & \mathbb{0} \\ \mathbb{0} & 1 \end{array} \right), \\ \tilde{S}_3 &= \left( \begin{array}{cc} \Phi_2 & \mathbb{0} \\ \mathbb{0} & 1 \end{array} \right) \left( \begin{array}{cc} \mathbb{K}'_3 & \mathbb{K}_4 \\ \mathbb{L}'_3 & \mathbb{L}_4 \end{array} \right)^{-1} \left( \begin{array}{cc} \mathbb{K}_3 & \mathbb{K}'_4 \\ \mathbb{L}_3 & \mathbb{L}'_4 \end{array} \right) \left( \begin{array}{cc} \Phi_2 & \mathbb{0} \\ \mathbb{0} & 1 \end{array} \right). \end{aligned} \quad (12) $$

In these matrices, we have introduced

$$ \begin{aligned} \mathbb{L}_i &= \sqrt{\varepsilon_i} \begin{pmatrix} \mathbb{B}_{xi} & \mathbb{A}_y \\ \mathbb{B}_{yi} & -\mathbb{A}_x \end{pmatrix}, & \mathbb{L}'_i &= \sqrt{\varepsilon_i} \begin{pmatrix} \mathbb{B}_{xi} & -\mathbb{A}_y \\ \mathbb{B}_{yi} & \mathbb{A}_x \end{pmatrix}, \\ \mathbb{K}_i &= \begin{pmatrix} \mathbb{A}_y & -\mathbb{B}_{xi} \\ -\mathbb{A}_x & -\mathbb{B}_{yi} \end{pmatrix}, & \mathbb{K}'_i &= \begin{pmatrix} -\mathbb{A}_y & -\mathbb{B}_{xi} \\ \mathbb{A}_x & -\mathbb{B}_{yi} \end{pmatrix}, \\ \mathbb{M}_i &= \sigma_g Z_0 \begin{pmatrix} -\mathbb{A}_x & -\mathbb{B}_{yi} \\ -\mathbb{A}_y & \mathbb{B}_{xi} \end{pmatrix}, & \mathbb{M}'_i &= \sigma_g Z_0 \begin{pmatrix} \mathbb{A}_x & -\mathbb{B}_{yi} \\ \mathbb{A}_y & \mathbb{B}_{xi} \end{pmatrix}, \end{aligned} \quad (13) $$

with

$$ \begin{aligned} \mathbb{A}_x &= \mathrm{diag} \left( \frac{k_{xn}}{k_n} \right)_n, & \mathbb{A}_y &= \mathrm{diag} \left( \frac{k_y}{k_n} \right)_n, \\ \mathbb{B}_{xi} &= \frac{c}{\sqrt{\varepsilon_i \omega}} \mathrm{diag} \left( \frac{k_{xn}}{k_n} k_{zn}^{(i)} \right)_n, & \\ \mathbb{B}_{yi} &= \frac{c}{\sqrt{\varepsilon_i \omega}} \mathrm{diag} \left( \frac{k_y}{k_n} k_{zn}^{(i)} \right)_n, & \end{aligned} \quad (14) $$

and

$$ \begin{aligned} \Phi_1 &= e^{\mathbb{D}h_1} = \left( \begin{array}{cc} e^{\mathbb{D}^{(11)}h_1} & \mathbb{0} \\ \mathbb{0} & e^{\mathbb{D}^{(22)}h_1} \end{array} \right), \\ \Phi_2 &= \left( \begin{array}{cc} \mathrm{diag} \left( e^{ik_{zn}^{(3)}h_2} \right)_n & \mathbb{0} \\ \mathbb{0} & \mathrm{diag} \left( e^{ik_{zn}^{(3)}h_2} \right)_n \end{array} \right). \end{aligned} \quad (15) $$

$Z_0$ denotes the impedance of vacuum.

The reflection matrix of body 1, denoted as $\mathcal{R}^{(1)+}$, can be derived directly from the matrix $\mathcal{R}^{-}$. Specifically, due to the distinct orientations along the z-axis, $\mathcal{R}^{(1)+}$ shares identical characteristics with $\mathcal{R}^{-}$ within the two diagonal blocks. However, a sign difference is observed within the two off-diagonal blocks, detailed in [16], as follows

$$ \mathcal{R}^{(1)+} = \begin{cases} \mathcal{R}_{p,p'}^{-} & p = p' \\ -\mathcal{R}_{p,p'}^{-} & p \neq p' \end{cases}. \quad (16) $$

For body 2, positioned at a distance $d$ from the origin, a phase shift is induced in the reflection matrix, as outlined in [39]. This reflection matrix can be formulated as

$$ \begin{aligned} & \langle p, \mathbf{k}, n | \mathcal{R}^{(2)-}(\omega) | p', \mathbf{k}', n' \rangle \\ &= e^{i(k_{zn} + k'_{zn'})d} \langle p, \mathbf{k}, n | \tilde{\mathcal{R}}^{(2)-}(\omega) | p', \mathbf{k}', n' \rangle. \end{aligned} \quad (17) $$

Here, $\tilde{\mathcal{R}}^{(2)-}$ is diagonal with dimensions $2(2N + 1) \times 2(2N + 1)$ and can be easily obtained using the Fresnel coefficients:

$$ \tilde{\mathcal{R}}^{(2)-} = \mathrm{diag}(\rho_n^{\mathrm{TE}}, \rho_n^{\mathrm{TM}}) \quad (18) $$

where

$$ \rho_n^{\mathrm{TE}} = \frac{k_{zn} - k_{zn}^{\mathrm{Au}}}{k_{zn} + k_{zn}^{\mathrm{Au}}} \quad (19) $$

$$ \rho_n^{\mathrm{TM}} = \frac{\varepsilon_{\mathrm{Au}}(\omega)k_{zn} - k_{zn}^{\mathrm{Au}}}{\varepsilon_{\mathrm{Au}}(\omega)k_z + k_z^{\mathrm{Au}}} \quad (20) $$

In these expressions, $k_{zn} = \sqrt{\frac{\omega^2}{c^2} - k_n^2}$.

Finally, it is worth noting that this calculation is also valid for imaginary Matsubara frequencies by replacing $\omega$ with $i\xi_m$.

# IV. RESULTS AND DISCUSSION

## A. The effect of graphene coating

Now that we have all the necessary ingredients to calculate the Casimir-Lifshitz pressure for the configuration described previously (see section II), let us discuss some numerical results. To begin with, we explore the impact on the CLP of covering the dielectric grating with graphene. The ratio of the pressure with graphene to that without graphene across various values of $\mu$ and $f$, for separation distances ranging from 60 nm to 8 $\mu$m, at temperatures $T = 300$ K and $T = 10$ K is depicted in Fig. 2.

At $T = 300$ K and for the slab configuration ($f = 1$), we observe in Fig. 2(a), that at short distances, below 60
where

$$ \Phi = \left( \begin{array}{cc} \mathrm{diag} \left( e^{ik_{zn}^{(4)}(h_1+h_2)} \right)_n & \mathbb{0} \\ \mathbb{0} & \mathrm{diag} \left( e^{ik_{zn}^{(4)}(h_1+h_2)} \right)_n \end{array} \right), \quad (9) $$

$k_{zn}^{(i)} = \sqrt{\varepsilon_i(\omega)\frac{\omega^2}{c^2} - \mathbf{k}_n^2}$ for medium i (see details in the appendix), and

$$ \tilde{\mathcal{S}} = \tilde{\mathcal{S}}_1 \circledast \tilde{\mathcal{S}}_2 \circledast \tilde{\mathcal{S}}_3. \quad (10) $$

Here, the star product operation $\mathcal{A} = \mathbb{B} \circledast \mathbb{C}$ is defined as [16]

$$ \begin{aligned} \mathcal{A}_{11} &= \mathbb{B}_{11} + \mathbb{B}_{12}(1 - \mathbb{C}_{11}\mathbb{B}_{22})^{-1}\mathbb{C}_{11}\mathbb{B}_{21}, \\ \mathcal{A}_{12} &= \mathbb{B}_{12}(1 - \mathbb{C}_{11}\mathbb{B}_{22})^{-1}\mathbb{C}_{12}, \\ \mathcal{A}_{21} &= \mathbb{C}_{21}(1 - \mathbb{B}_{22}\mathbb{C}_{11})^{-1}\mathbb{B}_{21}, \\ \mathcal{A}_{22} &= \mathbb{C}_{22} + \mathbb{C}_{21}(1 - \mathbb{B}_{22}\mathbb{C}_{11})^{-1}\mathbb{B}_{22}\mathbb{C}_{12}, \end{aligned} \quad (11) $$

and

$$ \begin{aligned} \tilde{\mathcal{S}}_1 &= \left( \begin{array}{cc} \mathbb{K}'_1 & -\mathbb{P} \\ \mathbb{L}'_1 + \mathbb{M}'_1 & -\mathbb{P}' \end{array} \right)^{-1} \left( \begin{array}{cc} \mathbb{K}_1 & \mathbb{P} \\ \mathbb{L}_1 + \mathbb{M}_1 & -\mathbb{P}' \end{array} \right), \\ \tilde{\mathcal{S}}_2 &= \left( \begin{array}{cc} \Phi_1 & \mathbb{0} \\ \mathbb{0} & 1 \end{array} \right) \left( \begin{array}{cc} -\mathbb{P} & -\mathbb{K}_3 \\ \mathbb{P}' & -\mathbb{L}_3 \end{array} \right)^{-1} \left( \begin{array}{cc} \mathbb{P} & -\mathbb{K}'_3 \\ \mathbb{P}' & -\mathbb{L}_3 \end{array} \right) \left( \begin{array}{cc} \Phi_1 & \mathbb{0} \\ \mathbb{0} & 1 \end{array} \right), \\ \tilde{\mathcal{S}}_3 &= \left( \begin{array}{cc} \Phi_2 & \mathbb{0} \\ \mathbb{0} & 1 \end{array} \right) \left( \begin{array}{cc} \mathbb{K}'_3 & \mathbb{K}_4 \\ \mathbb{L}'_3 & \mathbb{L}_4 \end{array} \right)^{-1} \left( \begin{array}{cc} \mathbb{K}_3 & \mathbb{K}'_4 \\ \mathbb{L}_3 & \mathbb{L}'_4 \end{array} \right) \left( \begin{array}{cc} \Phi_2 & \mathbb{0} \\ \mathbb{0} & 1 \end{array} \right). \end{aligned} \quad (12) $$

In these matrices, we have introduced

$$ \begin{aligned} \mathbb{L}_i &= \sqrt{\varepsilon_i} \begin{pmatrix} \mathbb{B}_{xi} & \mathbb{A}_y \\ \mathbb{B}_{yi} & -\mathbb{A}_x \end{pmatrix}, & \mathbb{L}'_i &= \sqrt{\varepsilon_i} \begin{pmatrix} \mathbb{B}_{xi} & -\mathbb{A}_y \\ \mathbb{B}_{yi} & \mathbb{A}_x \end{pmatrix}, \\ \mathbb{K}_i &= \begin{pmatrix} \mathbb{A}_y & -\mathbb{B}_{xi} \\ -\mathbb{A}_x & -\mathbb{B}_{yi} \end{pmatrix}, & \mathbb{K}'_i &= \begin{pmatrix} -\mathbb{A}_y & -\mathbb{B}_{xi} \\ \mathbb{A}_x & -\mathbb{B}_{yi} \end{pmatrix}, \\ \mathbb{M}_i &= \sigma_g Z_0 \begin{pmatrix} -\mathbb{A}_x & -\mathbb{B}_{yi} \\ -\mathbb{A}_y & \mathbb{B}_{xi} \end{pmatrix}, & \mathbb{M}'_i &= \sigma_g Z_0 \begin{pmatrix} \mathbb{A}_x & -\mathbb{B}_{yi} \\ \mathbb{A}_y & \mathbb{B}_{xi} \end{pmatrix}, \end{aligned} \quad (13) $$

with

$$ \begin{aligned} \mathbb{A}_x &= \mathrm{diag} \left( \frac{k_{xn}}{k_n} \right)_n, & \mathbb{A}_y &= \mathrm{diag} \left( \frac{k_y}{k_n} \right)_n, \\ \mathbb{B}_{xi} &= \frac{c}{\sqrt{\varepsilon_i \omega}} \mathrm{diag} \left( \frac{k_{xn}}{k_n} k_{zn}^{(i)} \right)_n, & \\ \mathbb{B}_{yi} &= \frac{c}{\sqrt{\varepsilon_i \omega}} \mathrm{diag} \left( \frac{k_y}{k_n} k_{zn}^{(i)} \right)_n, & \end{aligned} \quad (14) $$

and

$$ \begin{aligned} \Phi_1 &\equiv e^{\mathbb{D}h_1} = \begin{pmatrix} e^{\mathbb{D}^{(11)}h_1} & \mathbb{0} \\ \mathbb{0} & e^{\mathbb{D}^{(22)}h_1} \end{pmatrix}, \\ \Phi_2 &\equiv \begin{pmatrix} \mathrm{diag} \left( e^{ik_{zn}^{(3)}h_2} \right)_n & \mathbb{0} \\ \mathbb{0} & \mathrm{diag} \left( e^{ik_{zn}^{(3)}h_2} \right)_n \end{pmatrix}. \end{aligned} \quad (15) $$

$Z_0$ denotes the impedance of vacuum.

The reflection matrix of body 1, denoted as $\mathcal{R}^{(1)+}$, can be derived directly from the matrix $\mathcal{R}^{-}$. Specifically, due to the distinct orientations along the z-axis, $\mathcal{R}^{(1)+}$ shares identical characteristics with $\mathcal{R}^{-}$ within the two diagonal blocks. However, a sign difference is observed within the two off-diagonal blocks, detailed in [16], as follows

$$ \mathcal{R}^{(1)+} = \begin{cases} \mathcal{R}_{p,p'}^{-} & p = p' \\ -\mathcal{R}_{p,p'}^{-} & p \neq p' \end{cases}. \quad (16) $$

For body 2, positioned at a distance $d$ from the origin, a phase shift is induced in the reflection matrix, as outlined in [39]. This reflection matrix can be formulated as

$$ \langle p, \mathbf{k}, n | \tilde{\mathcal{R}}^{(2)-} (\omega) | p', \mathbf{k}', n' \rangle = e^{i(k_{zn} + k'_{zn'})d} \langle p, \mathbf{k}, n | \tilde{\mathcal{R}}^{(2)-} (\omega) | p', \mathbf{k}', n' \rangle. \quad (17) $$

Here, $\tilde{\mathcal{R}}^{(2)-}$ is diagonal with dimensions $2(2N + 1) \times 2(2N + 1)$ and can be easily obtained using the Fresnel coefficients:

$$ \tilde{\mathcal{R}}^{(2)-} = \mathrm{diag}(\rho_n^{\mathrm{TE}}, \rho_n^{\mathrm{TM}}) \quad (18) $$

where

$$ \rho_n^{\mathrm{TE}} = \frac{k_{zn} - k_{zn}^{\mathrm{Au}}}{k_{zn} + k_{zn}^{\mathrm{Au}}} \quad (19) $$

$$ \rho_n^{\mathrm{TM}} = \frac{\varepsilon_{\mathrm{Au}}(\omega)k_{zn} - k_{zn}^{\mathrm{Au}}}{\varepsilon_{\mathrm{Au}}(\omega)k_z + k_z^{\mathrm{Au}}} \quad (20) $$

In these expressions, $k_{zn} = \sqrt{\frac{\omega^2}{c^2} - \mathbf{k}_n^2}$.

Finally, it is worth noting that this calculation is also valid for imaginary Matsubara frequencies by replacing $\omega$ with $i\xi_m$.

# IV. RESULTS AND DISCUSSION

## A. The effect of graphene coating

Now that we have all the necessary ingredients to calculate the Casimir-Lifshitz pressure for the configuration described previously (see section II), let us discuss some numerical results. To begin with, we explore the impact on the CLP of covering the dielectric grating with graphene. The ratio of the pressure with graphene to that without graphene across various values of $\mu$ and $f$, for separation distances ranging from 60 nm to 8 $\mu$m, at temperatures $T = 300$ K and $T = 10$ K is depicted in Fig. 2.

At $T = 300$ K and for the slab configuration ($f = 1$), we observe in Fig. 2(a), that at short distances, below 60
nm, the CLP with graphene is approximately 10% higher than without graphene for $\mu = 0$ and $\mu = 0.2$ eV (solid lines with * and + markers), and 17% higher for $\mu = 0.5$ eV (solid line with square marker). This enhancement peaks at $d = 500$ nm before declining as $d$ increases. Beyond this separation, the CLP enhancement is almost identical for the three values of $\mu$. Figure 2(b) shows the ratio when the temperature is reduced to 10 K for the same range of $d$. The enhancement remains < 20%. For $\mu = 0.2$ eV and 0 eV, no peaks are observed. The ratio decreases to 5% at $d \approx 1\mu$m.

As the system cools down to $T = 10$ K, we observe a behavior similar to that at $T = 300$ K at short separation distances with minor variations. However, as the separation distance increases, the influence of graphene diminishes and becomes negligible, with this effect being less than 5% for $\mu = 0$ eV and $\mu = 0.2$ eV.

<figure><img src="image_3.png" /><figcaption>FIG. 2. Ratio of the CLP with graphene to CLP without graphene for different values of $\mu$ and $f$ at $T = 300$ K (a) and $T = 10$ K (b). Here, $D = 1 \mu$m, $h_1 = 300$ nm, and $h_2 = 100$ nm.</figcaption></figure>

On the other hand, the effect of graphene becomes more pronounced for the grating configuration. Considering $T = 300$ K and $f=0.2$, as shown Fig. 2(a), we clearly see significant changes on the CLP, with approximately a 100% increase for $\mu = 0$ eV (dotted line with *), 110% for $\mu = 0.2$ eV (dotted line with +), and 130% for $\mu = 0.5$ eV (dotted line with square) at short separation distances (i.e. below $d = 60$ nm). Similar trends persist at $T = 10$ K (Fig. 2(b)), where the relative change remains considerable, around 85% for $\mu = 0$ eV, 105% for $\mu = 0.2$ eV, and 130% for $\mu = 0.5$ eV. This enhancement decreases with increasing separation distance, reaching approximately 40% at $1 \mu$m for $T = 300$ K, irrespective of the chemical potential. However, at $T = 10$ K, the impact of graphene at this separation distance varies with the chemical potential, with enhancements of 13%, 20%, and 32% observed for $\mu = 0$ eV, $\mu = 0.2$ eV, and $\mu = 0.5$ eV, respectively. Additionally, the enhancement of the CLP by graphene continues to diminish with increasing separation distance at both temperatures.

<figure><img src="image_4.png" /><figcaption>FIG. 3. Casimir-Lifshitz pressure at $T = 300$ K (a) and $T = 10$ K (b), for $\mu = 0$ eV and different filling fractions $f = 0, 0.2, 0.5, 0.8$, and 1. The Lifshitz limit for metals, $P_{Lif}$, is shown as a black dashed line.</figcaption></figure>

To gain more insight into the order of magnitude of the CLP and how it varies with the geometric parameters and temperature, we plot the pressure $P$ as a function of $d$ on a log-log scale for $\mu = 0$ eV and at $T = 300$ K [Fig. 3(a)] and at $T = 10$ K [Fig. 3(b)]. These plots cover distances from 60 nm to $8 \mu$m for different filling fractions $f = 0, 0.2, 0.5, 0.8$, and 1. Notably, when $f = 0$, the configuration corresponds to the graphene sheet located on top of a vacuum gap 300 nm from the top surface of the $SiO_2$ layer with thickness $h_2 = 100$ nm, supported by a Silicon half-space. We observe that the CLP increases with the filling fraction $f$ regardless of temperature. Moreover, at $T = 300$ K, the CLP reaches the Lifshitz limit for metals ($P_{Lif} = -\frac{k_B T \zeta(3)}{8\pi d^3}$) at nearly 6 $\mu$m (for all filling fractions), as illustrated clearly in Fig. 3(a). However, at $T = 10$ K, depicted in Fig. 3(b), this limit remains unattained within the explored separation distances due to the larger thermal wavelength, $\lambda_T = \frac{hc}{k_B T}$, which is approximately 229 $\mu$m.

## B. Chemical potential and geometry effects

Next, we investigate the impact of the chemical potential combined with the geometry of the grating. For that, we analyze the ratio of the CLP at $\mu = 0.2$ eV and 0.5 eV to that at $\mu = 0$ eV, for various filling fractions $f = 0, 0.2, 0.5, 0.8$, and 1, and spanning distances from 60 nm to $8 \mu$m.

At $T = 300$ K and $\mu = 0.2$ eV, Fig. 4(a) shows that having a non-zero chemical potential always increases the pressure for any value of $f$. This augmentation is larger for smaller filling fraction. The ratio peaks at around 80 nm and decreases with distance. This effect is most pronounced when $f = 0$, reaching a maximum increase of ≈9% and for $f = 1$ (representing a slab configuration) the maximum increase is reduced to only 1%.

As the chemical potential increases to $\mu = 0.5$ eV, as shown in Fig. 4(b), the effect becomes more significant,
nm, the CLP with graphene is approximately 10% higher than without graphene for $\mu = 0$ and $\mu = 0.2$ eV (solid lines with * and + markers), and 17% higher for $\mu = 0.5$ eV (solid line with square marker). This enhancement peaks at $d = 500$ nm before declining as $d$ increases. Beyond this separation, the CLP enhancement is almost identical for the three values of $\mu$. Figure 2(b) shows the ratio when the temperature is reduced to 10 K for the same range of $d$. The enhancement remains < 20%. For $\mu = 0.2$ eV and 0 eV, no peaks are observed. The ratio decreases to 5% at $d \approx 1\mu$m.

As the system cools down to $T = 10$ K, we observe a behavior similar to that at $T = 300$ K at short separation distances with minor variations. However, as the separation distance increases, the influence of graphene diminishes and becomes negligible, with this effect being less than 5% for $\mu = 0$ eV and $\mu = 0.2$ eV.

<figure><img src="image_3.png" /><figcaption>FIG. 2. Ratio of the CLP with graphene to CLP without graphene for different values of $\mu$ and $f$ at $T = 300$ K (a) and $T = 10$ K (b). Here, $D = 1 \mu$m, $h_1 = 300$ nm, and $h_2 = 100$ nm.</figcaption></figure>

On the other hand, the effect of graphene becomes more pronounced for the grating configuration. Considering $T = 300$ K and $f=0.2$, as shown Fig. 2(a), we clearly see significant changes on the CLP, with approximately a 100% increase for $\mu = 0$ eV (dotted line with *), 110% for $\mu = 0.2$ eV (dotted line with +), and 130% for $\mu = 0.5$ eV (dotted line with square) at short separation distances (i.e. below $d = 60$ nm). Similar trends persist at $T = 10$ K (Fig. 2(b)), where the relative change remains considerable, around 85% for $\mu = 0$ eV, 105% for $\mu = 0.2$ eV, and 130% for $\mu = 0.5$ eV. This enhancement decreases with increasing separation distance, reaching approximately 40% at $1 \mu$m for $T = 300$ K, irrespective of the chemical potential. However, at $T = 10$ K, the impact of graphene at this separation distance varies with the chemical potential, with enhancements of 13%, 20%, and 32% observed for $\mu = 0$ eV, $\mu = 0.2$ eV, and $\mu = 0.5$ eV, respectively. Additionally, the enhancement of the CLP by graphene continues to diminish with increasing separation distance at both temperatures.

<figure><img src="image_4.png" /><figcaption>FIG. 3. Casimir-Lifshitz pressure at $T = 300$ K (a) and $T = 10$ K (b), for $\mu = 0$ eV and different filling fractions $f = 0, 0.2, 0.5, 0.8$, and 1. The Lifshitz limit for metals, $P_{Lif}$, is shown as a black dashed line.</figcaption></figure>

To gain more insight into the order of magnitude of the CLP and how it varies with the geometric parameters and temperature, we plot the pressure $P$ as a function of $d$ on a log-log scale for $\mu = 0$ eV and at $T = 300$ K [Fig. 3(a)] and at $T = 10$ K [Fig. 3(b)]. These plots cover distances from 60 nm to $8 \mu$m for different filling fractions $f = 0, 0.2, 0.5, 0.8$, and 1. Notably, when $f = 0$, the configuration corresponds to the graphene sheet located on top of a vacuum gap 300 nm from the top surface of the $SiO_2$ layer with thickness $h_2 = 100$ nm, supported by a Silicon half-space. We observe that the CLP increases with the filling fraction $f$ regardless of temperature. Moreover, at $T = 300$ K, the CLP reaches the Lifshitz limit for metals ($P_{Lif} = -\frac{k_B T \zeta(3)}{8\pi d^3}$) at nearly 6 $\mu$m (for all filling fractions), as illustrated clearly in Fig. 3(a). However, at $T = 10$ K, depicted in Fig. 3(b), this limit remains unattained within the explored separation distances due to the larger thermal wavelength, $\lambda_T = \frac{hc}{k_B T}$, which is approximately 229 $\mu$m.

## B. Chemical potential and geometry effects

Next, we investigate the impact of the chemical potential combined with the geometry of the grating. For that, we analyze the ratio of the CLP at $\mu = 0.2$ eV and 0.5 eV to that at $\mu = 0$ eV, for various filling fractions $f = 0, 0.2, 0.5, 0.8$, and 1, and spanning distances from 60 nm to $8 \mu$m.

At $T = 300$ K and $\mu = 0.2$ eV, Fig. 4(a) shows that having a non-zero chemical potential always increases the pressure for any value of $f$. This augmentation is larger for smaller filling fraction. The ratio peaks at around 80 nm and decreases with distance. This effect is most pronounced when $f = 0$, reaching a maximum increase of ≈9% and for $f = 1$ (representing a slab configuration) the maximum increase is reduced to only 1%.

As the chemical potential increases to $\mu = 0.5$ eV, as shown in Fig. 4(b), the effect becomes more significant,
with an increase of 32% for $f = 0$, 16% for $f = 0.2$ at $d \approx 100$ nm. The maximum ratio also becomes larger for other values of $f$.

<figure><img src="image_5.png" /><figcaption>FIG. 4. The ratio of CLP (calculated using Eq. (1)) at $\mu = 0.2$ eV (a) and $\mu = 0.5$ eV (b) to that at $\mu = 0$ eV, at $T = 300$ K, for various filling fractions $f = 0, 0.2, 0.5, 0.8$, and 1.</figcaption></figure>

Figure 5 shows the results when the temperature is reduced to 10 K. The most notable change is that the effect of $\mu$ becomes stronger. The ratio $P(\mu)/P(\mu=0)$ becomes larger for the chosen parameters. In particular, the maximum ratio increases for all values of $f$ compared to 300 K. For example, the peak ratio for $f = 0$ is 1.42 and 1.85 for $\mu = 0.2$ eV [Fig. 5(a)] and 0.5 eV [Fig. 5(b)] respectively at 10 K. Furthermore, it is apparent that the distance at which the ratio attains peak value is shifted towards larger values compared to 300K.

<figure><img src="image_6.png" /><figcaption>FIG. 5. The ratio of CLP (calculated using Eq. (1)) at $\mu = 0.2$ eV (a) and $\mu = 0.5$ eV (b) to that at $\mu = 0$ eV, at $T = 10$ K, for various filling fractions $f = 0, 0.2, 0.5, 0.8$, and 1.</figcaption></figure>

To further study the impact of geometry, we analyze the ratio of the CLP for a graphene-covered dielectric grating to that of both a graphene-covered slab and a bare slab (without graphene). In Fig. 6(a), where $T = 300$ K and $\mu = 0$ eV, the CLP between the graphene-covered grating and the gold half-space remains lower than that between a graphene-covered slab and a gold

<figure><img src="image_7.png" /><figcaption>FIG. 6. Normalized CLP for different filling fractions $f = 0, 0.2, 0.5, 0.8$, and 1 (if $\sigma_g \neq 0$) relative to the CLP for $f = 1$ with graphene (a) and without graphene (b) at $T = 300$ K and for $\mu = 0$ eV.</figcaption></figure>

half-space when the separation distance is less than 4 $\mu$m. Beyond 4 $\mu$m, the effect of $f$ diminishes, with the ratio approaching 1, since the two structures have the same asympthotique Lifshitz limit. Similarly, at $T = 10$ K, as depicted in Fig. 7(a), we observe a comparable trend with a greater saturation distance.

On the other hand, the CLP for the graphene-coated dielectric grating ($f \neq 1$) is initially lower than that for the bare slab at small separation distances, as illustrated in Fig. 6(b). However, as the separation distance $d$ increases, this trend reverses. For $f = 0.8$, the CLP exceeds that of the bare slab beyond 300 nm, while for $f = 0.5$, this occurs beyond 1 $\mu$m. As $f$ increases, the cross-over distance shifts to larger values. The crossover no longer occurs for $T = 10$ K, as it is clear from Fig. 7(b) where the CLP for the graphene-coated grating ($f \neq 1$) is lower than that of the bare slab. It is noteworthy that when the slab is covered with graphene ($f = 1$), the pressure is always higher compared to that without graphene whether at $T = 300$ K or $T = 10$ K.

### C. Thermal effect and Casimir force gradient

Now we explore the thermal impact of graphene by examining the ratio of the CLP at $T = 300$ K to that at $T = 10$ K. Initially, for $\mu = 0$ eV and $f = 0$, a notable thermal effect emerges around 300 nm, where the ratio of approximately 1.7 (Fig. 8(a)). As we increase the filling fraction, this peak gradually diminishes in height and shifts towards larger separation distances. For instance, the local peak decreases to approximately 1.4 for $f = 0.2$ (at $d=500$ nm) and 1.15 at $d = 1\mu$m for $f = 1$. After this peak (for every $f$), the ratio decreases, before eventually experiencing a sharp increase for distances $\approx 6\mu$m. The sharp rise is expected because the pressure at $T = 300$ K reaches the Lifshitz limit (see Fig. 3(a)) while remaining
with an increase of 32% for $f = 0$, 16% for $f = 0.2$ at $d \approx 100$ nm. The maximum ratio also becomes larger for other values of $f$.

<figure><img src="image_5.png" /><figcaption>FIG. 4. The ratio of CLP (calculated using Eq. (1)) at $\mu = 0.2$ eV (a) and $\mu = 0.5$ eV (b) to that at $\mu = 0$ eV, at $T = 300$ K, for various filling fractions $f = 0, 0.2, 0.5, 0.8$, and 1.</figcaption></figure>

Figure 5 shows the results when the temperature is reduced to 10 K. The most notable change is that the effect of $\mu$ becomes stronger. The ratio $P(\mu)/P(\mu=0)$ becomes larger for the chosen parameters. In particular, the maximum ratio increases for all values of $f$ compared to 300 K. For example, the peak ratio for $f = 0$ is 1.42 and 1.85 for $\mu = 0.2$ eV [Fig. 5(a)] and 0.5 eV [Fig. 5(b)] respectively at 10 K. Furthermore, it is apparent that the distance at which the ratio attains peak value is shifted towards larger values compared to 300K.

<figure><img src="image_6.png" /><figcaption>FIG. 5. The ratio of CLP (calculated using Eq. (1)) at $\mu = 0.2$ eV (a) and $\mu = 0.5$ eV (b) to that at $\mu = 0$ eV, at $T = 10$ K, for various filling fractions $f = 0, 0.2, 0.5, 0.8$, and 1.</figcaption></figure>

To further study the impact of geometry, we analyze the ratio of the CLP for a graphene-covered dielectric grating to that of both a graphene-covered slab and a bare slab (without graphene). In Fig. 6(a), where $T = 300$ K and $\mu = 0$ eV, the CLP between the graphene-covered grating and the gold half-space remains lower than that between a graphene-covered slab and a gold

<figure><img src="image_7.png" /><figcaption>FIG. 6. Normalized CLP for different filling fractions $f = 0, 0.2, 0.5, 0.8$, and 1 (if $\sigma_g \neq 0$) relative to the CLP for $f = 1$ with graphene (a) and without graphene (b) at $T = 300$ K and for $\mu = 0$ eV.</figcaption></figure>

half-space when the separation distance is less than 4 $\mu$m. Beyond 4 $\mu$m, the effect of $f$ diminishes, with the ratio approaching 1, since the two structures have the same asympthotique Lifshitz limit. Similarly, at $T = 10$ K, as depicted in Fig. 7(a), we observe a comparable trend with a greater saturation distance.

On the other hand, the CLP for the graphene-coated dielectric grating ($f \neq 1$) is initially lower than that for the bare slab at small separation distances, as illustrated in Fig. 6(b). However, as the separation distance $d$ increases, this trend reverses. For $f = 0.8$, the CLP exceeds that of the bare slab beyond 300 nm, while for $f = 0.5$, this occurs beyond 1 $\mu$m. As $f$ increases, the cross-over distance shifts to larger values. The crossover no longer occurs for $T = 10$ K, as it is clear from Fig. 7(b) where the CLP for the graphene-coated grating ($f \neq 1$) is lower than that of the bare slab. It is noteworthy that when the slab is covered with graphene ($f = 1$), the pressure is always higher compared to that without graphene whether at $T = 300$ K or $T = 10$ K.

### C. Thermal effect and Casimir force gradient

Now we explore the thermal impact of graphene by examining the ratio of the CLP at $T = 300$ K to that at $T = 10$ K. Initially, for $\mu = 0$ eV and $f = 0$, a notable thermal effect emerges around 300 nm, where the ratio of approximately 1.7 (Fig. 8(a)). As we increase the filling fraction, this peak gradually diminishes in height and shifts towards larger separation distances. For instance, the local peak decreases to approximately 1.4 for $f = 0.2$ (at $d=500$ nm) and 1.15 at $d = 1\mu$m for $f = 1$. After this peak (for every $f$), the ratio decreases, before eventually experiencing a sharp increase for distances $\approx 6\mu$m. The sharp rise is expected because the pressure at $T = 300$ K reaches the Lifshitz limit (see Fig. 3(a)) while remaining
<figure><img src="image_8.png" /><figcaption>FIG. 7. Normalized CLP for different filling fractions $f = 0, 0.2, 0.5, 0.8$, and $1$ (if $\sigma_g \neq 0$) relative to the CLP for $f = 1$ with graphene (a) and without graphene (b) at $T = 10$ K and for $\mu = 0$ eV.</figcaption></figure>

below it at $T = 10$ K (see Fig. 3(b)). In this limiting case, using the expression $P_{\text{Lif}} = k_B T \zeta(3)/(8\pi d^3)$ yields a ratio of 30. Consequently, it is expected that the ratio continues to increase until reaching the thermal wavelength, approximately 229 µm, beyond which the value levels off.

Moreover, this thermal effect decreases with the chemical potential, as demonstrated in Fig. 8(b) and 8(c) for $\mu = 0.2$ and $0.5$, respectively. This is due to the fact that the graphene conductivity (6) only depends weakly on the temperature $T \approx 10^{-3} \text{eV}$ when $\mu$ is larger than $T$.

<figure><img src="image_9.png" /><figcaption>FIG. 8. Ratio of the CLP at $T = 300$ K to that at $T = 10$ K for $\mu = 0$ eV (a), $\mu = 0.2$ eV (b), and $\mu = 0.5$ eV (c), for different filling fractions $f = 0, 0.2, 0.5, 0.8$, and 1.</figcaption></figure>

Finally, we consider the gradient of the Casimir force between a gold sphere and a dielectric grating covered with graphene. According to the Proximity Force Approximation (PFA) [11, 14], the Casimir pressure $P$ between two parallel plates is proportional to the gradient of the Casimir force in the sphere-plate geometry, and it can be expressed as $F' = 2\pi RP$, where $R$ is the sphere

radius. The latter configuration is used in experiments to avoid the difficulty in keeping two planar bodies parallel at short separations. We considered a gold sphere with a radius of 40 µm at $T = 300$ K and $T = 10$ K, and calculated the gradient over distances ranging from 60 nm to 1 µm for two configurations: a slab ($f = 1$) and a grating with $f = 0.2$ and $\mu = 0$ eV.

Our analysis revealed that the effect of the graphene on the force gradient is significantly more pronounced for a grating compared to solid SiO$_2$, particularly at short separation distances for both $T = 10$ K and 300 K. the force gradient at $d= 200$ nm increases by 94% and 55% respectively for $T = 300$ K and 10K when graphene is added to an SiO$_2$ grating with $f = 0.2$.

<figure><img src="image_10.png" /><figcaption>FIG. 9. Gradient of the Casimir force between a gold sphere with radius $R = 40$ µm and a graphene-coated fused silica grating with the following parameters: period $D = 1$ µm, depth $h_1 = 300$ nm, and slab thickness $h_2 = 100$ nm, at $T = 300$ K (a) and $T = 10$ K (b). for different values of $\mu$ and $f$.</figcaption></figure>

V. CONCLUSION

In conclusion, we have investigated the Casimir-Lifshitz pressure between a graphene-coated dielectric grating and a gold half-space using an adapted FMM. Our results reveal that covering the dielectric grating with graphene leads to a significant increase in pressure, with enhancements up to 130% for $\mu = 0.5$ eV and $f = 0.2$, compared to only 17% when covering a slab with graphene, over a wide range of temperatures from 10K to 300K. Additionally, we have shown that the CLP depends strongly on the chemical potential of graphene, with maximal effects observed at lower filling fractions ($f = 0.2$ in our case) for both $T = 300$ K and 10 K. Notably, our findings indicate that the pressure for a graphene-coated dielectric grating can surpass that of a bare slab at $T = 300$ K at certain separation distances even when $f$ is small. Furthermore, we have identified a thermal effect
<figure><img src="image_8.png" /><figcaption>FIG. 7. Normalized CLP for different filling fractions $f = 0, 0.2, 0.5, 0.8$, and $1$ (if $\sigma_g \neq 0$) relative to the CLP for $f = 1$ with graphene (a) and without graphene (b) at $T = 10$ K and for $\mu = 0$ eV.</figcaption></figure>

below it at $T = 10$ K (see Fig. 3(b)). In this limiting case, using the expression $P_{\text{Lif}} = k_B T \zeta(3)/(8\pi d^3)$ yields a ratio of 30. Consequently, it is expected that the ratio continues to increase until reaching the thermal wavelength, approximately 229 µm, beyond which the value levels off.

Moreover, this thermal effect decreases with the chemical potential, as demonstrated in Fig. 8(b) and 8(c) for $\mu = 0.2$ and $0.5$, respectively. This is due to the fact that the graphene conductivity (6) only depends weakly on the temperature $T \approx 10^{-3} \text{eV}$ when $\mu$ is larger than $T$.

<figure><img src="image_9.png" /><figcaption>FIG. 8. Ratio of the CLP at $T = 300$ K to that at $T = 10$ K for $\mu = 0$ eV (a), $\mu = 0.2$ eV (b), and $\mu = 0.5$ eV (c), for different filling fractions $f = 0, 0.2, 0.5, 0.8$, and $1$.</figcaption></figure>

Finally, we consider the gradient of the Casimir force between a gold sphere and a dielectric grating covered with graphene. According to the Proximity Force Approximation (PFA) [11, 14], the Casimir pressure $P$ between two parallel plates is proportional to the gradient of the Casimir force in the sphere-plate geometry, and it can be expressed as $F' = 2\pi RP$, where $R$ is the sphere

radius. The latter configuration is used in experiments to avoid the difficulty in keeping two planar bodies parallel at short separations. We considered a gold sphere with a radius of 40 µm at $T = 300$ K and $T = 10$ K, and calculated the gradient over distances ranging from 60 nm to 1 µm for two configurations: a slab ($f = 1$) and a grating with $f = 0.2$ and $\mu = 0$ eV.

Our analysis revealed that the effect of the graphene on the force gradient is significantly more pronounced for a grating compared to solid SiO₂, particularly at short separation distances for both $T = 10$ K and 300 K. the force gradient at $d= 200$ nm increases by 94% and 55% respectively for $T = 300$ K and 10K when graphene is added to an SiO₂ grating with $f = 0.2$.

<figure><img src="image_10.png" /><figcaption>FIG. 9. Gradient of the Casimir force between a gold sphere with radius $R = 40$ µm and a graphene-coated fused silica grating with the following parameters: period $D = 1$ µm, depth $h_1 = 300$ nm, and slab thickness $h_2 = 100$ nm, at $T = 300$ K (a) and $T = 10$ K (b). for different values of $\mu$ and $f$.</figcaption></figure>

V. CONCLUSION

In conclusion, we have investigated the Casimir-Lifshitz pressure between a graphene-coated dielectric grating and a gold half-space using an adapted FMM. Our results reveal that covering the dielectric grating with graphene leads to a significant increase in pressure, with enhancements up to 130% for $\mu = 0.5$ eV and $f = 0.2$, compared to only 17% when covering a slab with graphene, over a wide range of temperatures from 10K to 300K. Additionally, we have shown that the CLP depends strongly on the chemical potential of graphene, with maximal effects observed at lower filling fractions ($f = 0.2$ in our case) for both $T = 300$ K and 10 K. Notably, our findings indicate that the pressure for a graphene-coated dielectric grating can surpass that of a bare slab at $T = 300$ K at certain separation distances even when $f$ is small. Furthermore, we have identified a thermal effect
between 300 nm and 400 nm, which diminishes with increasing filling factor and chemical potential. Lastly, we have presented the Casimir force gradient between a gold sphere and a graphene-coated dielectric grating, a configuration that is used for experimental measurements.

## ACKNOWLEDGMENTS

The work described in this paper was supported by a grant "CAT" from the ANR/RGC Joint Research Scheme sponsored by the French National Research Agency (ANR) and the Research Grants Council (RGC) of the Hong Kong Special Administrative Region, China (Project No. A-HKUST604/20).

## APPENDIX: SCATTERING MATRIX CALCULATION FOR GRAPHENE-COVERED GRATING USING THE FOURIER MODAL METHOD

In this appendix, we will describe the calculation of the scattering matrix using the Fourier Modal Method for dielectric gratings coated with graphene, employing a zero-thickness model. This approach simplifies the representation of graphene and facilitates its direct inclusion into the boundary conditions of the FMM. The structure under study is depicted in Fig. 10, where a graphene sheet is positioned at $z = 0$. Four distinct zones can be distinguished: (i) zone 1 ($z < 0$) is characterized by $\varepsilon_1$, (ii) zone 2 ($0 < z < h_1$) contains a lamellar grating with a period $D$ along the x-axis, a relative dielectric permittivity defined denoted $\varepsilon_2(x)$ and a thickness $h_1$, (iii) zone 3 ($h_1 < z < h_1+h_2$), is a homogeneous slab with permittivity $\varepsilon_3$ and (iv) zone 4 ($z > h_1+h_2$) is a semi-infinite medium characterized by $\varepsilon_4$.

### A. Electromagnetic Fields

Due to periodicity along the x-axis, the electromagnetic fields in each homogeneous zone ($i = 1, 3, 4$) can be expanded into the so-called Rayleigh expansion as follows:

<figure><img src="image_11.png" /><figcaption>FIG. 10. Geometry of the FMM calculation: A dielectric grating characterized by dielectric permittivity $\varepsilon_2(x)$, with period $D$, width $a$, and depth $h_1$, covered with a graphene sheet and placed on a top of a slab with dielectric permittivity $\varepsilon_3$ and thickness $h_2$, which in turn rests on a dielectric half-space with dielectric permittivity $\varepsilon_4$.</figcaption></figure>

$$ egin{cases} 
\mathbf{E}^{(i)}(\mathbf{R}, \omega) = \sum_{p, \phi} \int_{-\frac{\pi}{D}}^{\frac{\pi}{D}} \frac{dk_x}{2\pi} \sum_{n \in \mathbb{Z}} \int_{-\infty}^{+\infty} \frac{dk_y}{2\pi} \\ 
\times e^{i\mathbf{K}_n^{(i)}\phi \cdot \mathbf{R}} \hat{\mathbf{e}}_p^{(i)}\phi(\mathbf{k}_n, \omega) E_p^{(i)}\phi(\mathbf{k}_n, \omega), \\[2ex] 
\mathbf{B}^{(i)}(\mathbf{R}, \omega) = \frac{\sqrt{\varepsilon_i(\omega)}}{c} \sum_{p, \phi} \int_{-\frac{\pi}{D}}^{\frac{\pi}{D}} \frac{dk_x}{2\pi} \sum_{n \in \mathbb{Z}} \int_{-\infty}^{+\infty} \frac{dk_y}{2\pi} \\ 
\times e^{i\mathbf{K}_n^{(i)}\phi \cdot \mathbf{R}} (-1)^p \hat{\mathbf{e}}_{S(p)}^{(i)}\phi(\mathbf{k}_n, \omega) E_p^{(i)}\phi(\mathbf{k}_n, \omega), 
\end{cases} 
\quad (21) $$

where $p$ denotes the index of polarization ($p=1,2$ for TE and TM respectively), $\mathbf{R} = (\mathbf{r}, z)$, $S(p)$ is a function ($S(1) = 2$ and $S(2) = 1$), and

$$ egin{cases} 
\hat{\mathbf{e}}_{\text{TE}}^{(i)\phi}(\mathbf{k}_n, \omega) = \frac{1}{k_n}(-k_y \hat{\mathbf{e}}_x + k_{xn} \hat{\mathbf{e}}_y), \\[2ex] 
\hat{\mathbf{e}}_{\text{TM}}^{(i)\phi}(\mathbf{k}_n, \omega) = \frac{c}{\omega\sqrt{\varepsilon_i(\omega)}}(-k_n \hat{\mathbf{e}}_z + \phi k_{zn}^i \hat{\mathbf{k}}_n). 
\end{cases} 
\quad (22) $$

Here, $\hat{\mathbf{e}}_x$, $\hat{\mathbf{e}}_y$ and $\hat{\mathbf{e}}_z$ are unit vectors in the $(x, y, z)$ Cartesian basis, $\hat{\mathbf{k}}_n = \mathbf{k}_n/k_n$, $\phi$ is the direction of propagation of the waves $(+, -)$ along the z-axis for the incident and the reflected fields, respectively, and $\mathbf{K}_n^{(i)\phi} = (\mathbf{k}_n, \phi k_{zn}^i)$.

For the different amplitudes in these zones, we use the notations:

$$ egin{cases} 
E^{(1)+} = I, & E^{(1)-} = R, \\ 
E^{(3)+} = C, & E^{(3)-} = C', \\ 
E^{(4)+} = T, & E^{(4)-} = J, 
\end{cases} 
\quad (23) $$

In the periodic region (zone 2), the electric field assumes the following form:

$$ \mathbf{E}^{(2)}(\mathbf{R}, \omega) = \int_{-\frac{\pi}{D}}^{\frac{\pi}{D}} \frac{dk_x}{2\pi} \sum_{n \in \mathbb{Z}} \int_{-\infty}^{+\infty} \frac{dk_y}{2\pi} e^{i\mathbf{k}_n \cdot \mathbf{r}} \mathbf{E}^{(2)}(z, \mathbf{k}_n, \omega). \quad (24) $$
between 300 nm and 400 nm, which diminishes with increasing filling factor and chemical potential. Lastly, we have presented the Casimir force gradient between a gold sphere and a graphene-coated dielectric grating, a configuration that is used for experimental measurements.

## ACKNOWLEDGMENTS

The work described in this paper was supported by a grant "CAT" from the ANR/RGC Joint Research Scheme sponsored by the French National Research Agency (ANR) and the Research Grants Council (RGC) of the Hong Kong Special Administrative Region, China (Project No. A-HKUST604/20).

## APPENDIX: SCATTERING MATRIX CALCULATION FOR GRAPHENE-COVERED GRATING USING THE FOURIER MODAL METHOD

In this appendix, we will describe the calculation of the scattering matrix using the Fourier Modal Method for dielectric gratings coated with graphene, employing a zero-thickness model. This approach simplifies the representation of graphene and facilitates its direct inclusion into the boundary conditions of the FMM. The structure under study is depicted in Fig. 10, where a graphene sheet is positioned at $z = 0$. Four distinct zones can be distinguished: (i) zone 1 ($z < 0$) is characterized by $\varepsilon_1$, (ii) zone 2 ($0 < z < h_1$) contains a lamellar grating with a period $D$ along the x-axis, a relative dielectric permittivity defined denoted $\varepsilon_2(x)$ and a thickness $h_1$, (iii) zone 3 ($h_1 < z < h_1+h_2$), is a homogeneous slab with permittivity $\varepsilon_3$ and (iv) zone 4 ($z > h_1+h_2$) is a semi-infinite medium characterized by $\varepsilon_4$.

### A. Electromagnetic Fields

Due to periodicity along the x-axis, the electromagnetic fields in each homogeneous zone ($i = 1, 3, 4$) can be expanded into the so-called Rayleigh expansion as follows:

<figure><img src="image_11.png" /><figcaption>FIG. 10. Geometry of the FMM calculation: A dielectric grating characterized by dielectric permittivity $\varepsilon_2(x)$, with period $D$, width $a$, and depth $h_1$, covered with a graphene sheet and placed on a top of a slab with dielectric permittivity $\varepsilon_3$ and thickness $h_2$, which in turn rests on a dielectric half-space with dielectric permittivity $\varepsilon_4$.</figcaption></figure>

$$ egin{cases} 
\mathbf{E}^{(i)}(\mathbf{R}, \omega) = \sum_{p, \phi} \int_{-\frac{\pi}{D}}^{\frac{\pi}{D}} \frac{dk_x}{2\pi} \sum_{n \in \mathbb{Z}} \int_{-\infty}^{+\infty} \frac{dk_y}{2\pi} \\ 
\times e^{i\mathbf{K}_n^{(i)}\phi \cdot \mathbf{R}} \hat{\mathbf{e}}_p^{(i)}\phi(\mathbf{k}_n, \omega) E_p^{(i)}\phi(\mathbf{k}_n, \omega), \\[2em] 
\mathbf{B}^{(i)}(\mathbf{R}, \omega) = \frac{\sqrt{\varepsilon_i(\omega)}}{c} \sum_{p, \phi} \int_{-\frac{\pi}{D}}^{\frac{\pi}{D}} \frac{dk_x}{2\pi} \sum_{n \in \mathbb{Z}} \int_{-\infty}^{+\infty} \frac{dk_y}{2\pi} \\ 
\times e^{i\mathbf{K}_n^{(i)}\phi \cdot \mathbf{R}} (-1)^p \hat{\mathbf{e}}_{S(p)}^{(i)}\phi(\mathbf{k}_n, \omega) E_p^{(i)}\phi(\mathbf{k}_n, \omega), 
\end{cases} 
\quad (21) $$

where $p$ denotes the index of polarization ($p=1,2$ for TE and TM respectively), $\mathbf{R} = (\mathbf{r}, z)$, $S(p)$ is a function ($S(1) = 2$ and $S(2) = 1$), and

$$ egin{cases} 
\hat{\mathbf{e}}_{\text{TE}}^{(i)\phi}(\mathbf{k}_n, \omega) = \frac{1}{k_n}(-k_y \hat{\mathbf{e}}_x + k_{xn} \hat{\mathbf{e}}_y), \\[2em] 
\hat{\mathbf{e}}_{\text{TM}}^{(i)\phi}(\mathbf{k}_n, \omega) = \frac{c}{\omega\sqrt{\varepsilon_i(\omega)}}(-k_n \hat{\mathbf{e}}_z + \phi k_{zn}^i \hat{\mathbf{k}}_n). 
\end{cases} 
\quad (22) $$

Here, $\hat{\mathbf{e}}_x$, $\hat{\mathbf{e}}_y$ and $\hat{\mathbf{e}}_z$ are unit vectors in the $(x, y, z)$ Cartesian basis, $\hat{\mathbf{k}}_n = \mathbf{k}_n/k_n$, $\phi$ is the direction of propagation of the waves $(+, -)$ along the z-axis for the incident and the reflected fields, respectively, and $\mathbf{K}_n^{(i)\phi} = (\mathbf{k}_n, \phi k_{zn}^i)$.

For the different amplitudes in these zones, we use the notations:

$$ egin{cases} 
E^{(1)+} = I, & E^{(1)-} = R, \\ 
E^{(3)+} = C, & E^{(3)-} = C', \\ 
E^{(4)+} = T, & E^{(4)-} = J, 
\end{cases} 
\quad (23) $$

In the periodic region (zone 2), the electric field assumes the following form:

$$ \mathbf{E}^{(2)}(\mathbf{R}, \omega) = \int_{-\frac{\pi}{D}}^{\frac{\pi}{D}} \frac{dk_x}{2\pi} \sum_{n \in \mathbb{Z}} \int_{-\infty}^{+\infty} \frac{dk_y}{2\pi} e^{i\mathbf{k}_n \cdot \mathbf{r}} \mathbf{E}^{(2)}(z, \mathbf{k}_n, \omega). \quad (24) $$
In the actual computation, we will truncate the sum ranging from $-\infty$ to $\infty$, retaining only $2N + 1$ Fourier harmonics, $N$ is termed the truncation order. In order to obtain the amplitudes $\mathbf{E}^{(2)}(z, \mathbf{k}_n, \omega)$, we have to solve Maxwell's equations within this zone:

$$ \left\{ \begin{array}{l} \partial_y E_z - \partial_z E_y = ik_0 \tilde{H}_x \\ \partial_z E_x - \partial_x E_z = ik_0 \tilde{H}_y \\ \partial_x E_y - \partial_y E_x = ik_0 \tilde{H}_z \end{array} \right. \left\{ \begin{array}{l} \partial_y \tilde{H}_z - \partial_z \tilde{H}_y = -i\varepsilon k_0 E_x \\ \partial_z \tilde{H}_x - \partial_x \tilde{H}_z = -i\varepsilon k_0 E_y \\ \partial_x \tilde{H}_y - \partial_y \tilde{H}_x = -i\varepsilon k_0 E_z, \end{array} \right. \quad (25) $$

with $k_0 = \omega/c$, and $\tilde{H} = Z_0 H$.

Eliminating the $z$ components from Eq. (25), we obtain two systems of equations involving only the parallel components over which the boundary conditions hold:

$$ \partial_z \begin{pmatrix} E_x \\ E_y \end{pmatrix} = \begin{pmatrix} -\frac{i}{k_0} \partial_x \frac{1}{\varepsilon(x)} \partial_y & ik_0 + \frac{i}{k_0} \partial_x \frac{1}{\varepsilon(x)} \partial_x \\ -ik_0 - \frac{i}{k_0} \partial_y \frac{1}{\varepsilon(x)} \partial_y & \frac{i}{k_0} \partial_y \frac{1}{\varepsilon(x)} \partial_x \end{pmatrix} \times \begin{pmatrix} \tilde{H}_x \\ \tilde{H}_y \end{pmatrix} \quad (26) $$

$$ \partial_z \begin{pmatrix} \tilde{H}_x \\ \tilde{H}_y \end{pmatrix} = \begin{pmatrix} \frac{i}{k_0} \partial_x \partial_y & -ik_0 \varepsilon(x) - \frac{i}{k_0} \partial_x \partial_x \\ ik_0 \varepsilon(x) + \frac{i}{k_0} \partial_y \partial_y & -\frac{i}{k_0} \partial_x \partial_y \end{pmatrix} \times \begin{pmatrix} E_x \\ E_y \end{pmatrix} \quad (27) $$

Subsequently, following the method elaborated in more detail, for instance, in [41], we transform these equations into Fourier space and get:

$$ \left\{ \begin{array}{l} \partial_z \mathcal{E} = \mathbb{F} \tilde{\mathcal{H}} \\ \partial_z \tilde{\mathcal{H}} = \mathcal{G} \mathcal{E} \end{array} \right. \quad (28) $$

where $\mathcal{E} = [\mathcal{E}_x, \mathcal{E}_y]^t$ and $\tilde{\mathcal{H}} = [\tilde{\mathcal{H}}_x, \tilde{\mathcal{H}}_y]^t$ are vectors containing the Fourier components of the electric and magnetic parallel components of the field in the grating zone. Moreover,

$$ \mathbb{F} = \begin{pmatrix} \frac{i\beta}{k_0}\alpha[\varepsilon]^{-1} & ik_0\mathbb{1} - \frac{i\alpha}{k_0}[\varepsilon]^{-1}\alpha \\ -ik_0\mathbb{1} + \frac{i\beta^2}{k_0}[\varepsilon]^{-1} & -\frac{i\beta}{k_0}[\varepsilon]^{-1}\alpha \end{pmatrix}, \quad (29) $$

and

$$ \mathcal{G} = \begin{pmatrix} -\frac{i\beta}{k_0}\alpha & -ik_0[\varepsilon] + \frac{i\alpha^2}{k_0} \\ ik_0[\varepsilon^{-1}]^{-1} & -\frac{i\beta^2}{k_0} + \frac{i\beta}{k_0}\alpha \end{pmatrix}. \quad (30) $$

Here, $[\varepsilon]$ is the Toeplitz matrix such that $[\varepsilon]_{mn} = \varepsilon_{m-n}$, $\beta$ is a scalar, and $\alpha = \operatorname{diag}(k_{xn})_n$.

Using the coupled equations (28), we deduce the equation for $\mathcal{E}$ alone:

$$ \partial_z^2 \mathcal{E} = \mathbb{F}\mathcal{G}\mathcal{E} = \mathbb{P}\mathbb{D}^2\mathbb{P}^{-1}\mathcal{E}, \quad (31) $$

where $\mathbb{P}$ and $\mathbb{D}^2$ contain, respectively, the eigenvectors and eigenvalues of the matrix $\mathbb{F}\mathcal{G}$, each having dimensions of $2(2N + 1) \times 2(2N + 1)$ and take the following form:

$$ \mathbb{P} = \begin{pmatrix} \mathbb{P}^{(11)} & \mathbb{P}^{(12)} \\ \mathbb{P}^{(21)} & \mathbb{P}^{(22)} \end{pmatrix}, \quad \mathbb{D} = \begin{pmatrix} \mathbb{D}^{(11)} & \mathbb{0} \\ \mathbb{0} & \mathbb{D}^{(22)} \end{pmatrix}. \quad (32) $$

The solution of the set of equations (28) can then be written:

$$ \left\{ \begin{array}{l} \mathcal{E} = \mathbb{P}(e^{\mathbb{D}z}\mathcal{A} + e^{-\mathbb{D}z}\mathcal{B}) \\ \tilde{\mathcal{H}} = \mathbb{P}'(e^{\mathbb{D}z}\mathcal{A} - e^{-\mathbb{D}z}\mathcal{B}) \end{array} \right. \quad (33) $$

where $\mathcal{A} = [A_{xn}, A_{yn}]^t$ and $\mathcal{B} = [B_{xn}, B_{yn}]^t$ ($n \in [-N, N]$) are arbitrary constant vectors, and $\mathbb{P}' = \mathbb{F}^{-1}\mathbb{P}\mathbb{D}$.

## B. Boundary conditions

After expressing the fields in all zones, we are going to enforce the boundary conditions at their interfaces. For both analytical and numerical convenience, we introduce an additional phase factor in the expressions of the fields within zones 3 and 4. Specifically in zone 3, we take the $z$ dependence of the fields under the form $e^{ik^{(i)}\phi(z-h_1)}$, while in zone 4, we take it under the form $e^{ik^{(i)}\phi(z-h_1-h_2)}$. These adjustments facilitate the calculations and can be easily reverted at the conclusion of the analysis. For this reason, we will first obtain the scattering matrix with this simplified phase $\tilde{\mathcal{S}}$, and then multiply it by the accumulated phase to obtain the actual scattering matrix $\mathcal{S}$ (see Eq. (8)).

Using the zero thickness model for the graphene sheet and applying the boundary conditions for the electric field components at the interface ($z=0$) between zones 1 and 2, in the cartesian coordinates system, we obtain (after projection on the Fourier basis):
In the actual computation, we will truncate the sum ranging from $-\infty$ to $\infty$, retaining only $2N + 1$ Fourier harmonics, $N$ is termed the truncation order. In order to obtain the amplitudes $\mathbf{E}^{(2)}(z, \mathbf{k}_n, \omega)$, we have to solve Maxwell's equations within this zone:

$$ \left\{ \begin{array}{l} \partial_y E_z - \partial_z E_y = ik_0 \tilde{H}_x \\ \partial_z E_x - \partial_x E_z = ik_0 \tilde{H}_y \\ \partial_x E_y - \partial_y E_x = ik_0 \tilde{H}_z \end{array} \right. \left\{ \begin{array}{l} \partial_y \tilde{H}_z - \partial_z \tilde{H}_y = -i\varepsilon k_0 E_x \\ \partial_z \tilde{H}_x - \partial_x \tilde{H}_z = -i\varepsilon k_0 E_y \\ \partial_x \tilde{H}_y - \partial_y \tilde{H}_x = -i\varepsilon k_0 E_z, \end{array} \right. \quad (25) $$

with $k_0 = \omega/c$, and $\tilde{H} = Z_0 H$.

Eliminating the $z$ components from Eq. (25), we obtain two systems of equations involving only the parallel components over which the boundary conditions hold:

$$ \partial_z \begin{pmatrix} E_x \\ E_y \end{pmatrix} = \begin{pmatrix} -\frac{i}{k_0} \partial_x \frac{1}{\varepsilon(x)} \partial_y & ik_0 + \frac{i}{k_0} \partial_x \frac{1}{\varepsilon(x)} \partial_x \\ -ik_0 - \frac{i}{k_0} \partial_y \frac{1}{\varepsilon(x)} \partial_y & \frac{i}{k_0} \partial_y \frac{1}{\varepsilon(x)} \partial_x \end{pmatrix} \times \begin{pmatrix} \tilde{H}_x \\ \tilde{H}_y \end{pmatrix} \quad (26) $$

$$ \partial_z \begin{pmatrix} \tilde{H}_x \\ \tilde{H}_y \end{pmatrix} = \begin{pmatrix} \frac{i}{k_0} \partial_x \partial_y & -ik_0 \varepsilon(x) - \frac{i}{k_0} \partial_x \partial_x \\ ik_0 \varepsilon(x) + \frac{i}{k_0} \partial_y \partial_y & -\frac{i}{k_0} \partial_x \partial_y \end{pmatrix} \times \begin{pmatrix} E_x \\ E_y \end{pmatrix} \quad (27) $$

Subsequently, following the method elaborated in more detail, for instance, in [41], we transform these equations into Fourier space and get:

$$ \left\{ \begin{array}{l} \partial_z \mathcal{E} = \mathbb{F} \tilde{\mathcal{H}} \\ \partial_z \tilde{\mathcal{H}} = \mathcal{G} \mathcal{E} \end{array} \right. \quad (28) $$

where $\mathcal{E} = [\mathcal{E}_x, \mathcal{E}_y]^t$ and $\tilde{\mathcal{H}} = [\tilde{\mathcal{H}}_x, \tilde{\mathcal{H}}_y]^t$ are vectors containing the Fourier components of the electric and magnetic parallel components of the field in the grating zone. Moreover,

$$ \mathbb{F} = \begin{pmatrix} \frac{i\beta}{k_0} \alpha \, [\varepsilon]^{-1} & ik_0 \mathbb{1} - \frac{i\alpha}{k_0} \, [\varepsilon]^{-1} \alpha \\ -ik_0 \mathbb{1} + \frac{i\beta^2}{k_0} \, [\varepsilon]^{-1} & -\frac{i\beta}{k_0} \, [\varepsilon]^{-1} \alpha \end{pmatrix}, \quad (29) $$

and

$$ \mathcal{G} = \begin{pmatrix} -\frac{i\beta}{k_0}\alpha & -ik_0[\varepsilon] + \frac{i\alpha^2}{k_0} \\ ik_0[\varepsilon^{-1}]^{-1} & -\frac{i\beta^2}{k_0} + \frac{i\beta}{k_0}\alpha \end{pmatrix}. \quad (30) $$

Here, $[\varepsilon]$ is the Toeplitz matrix such that $[\varepsilon]_{mn} = \varepsilon_{m-n}$, $\beta$ is a scalar, and $\alpha = \operatorname{diag}(k_{xn})_n$.

Using the coupled equations (28), we deduce the equation for $\mathcal{E}$ alone:

$$ \partial_z^2 \mathcal{E} = \mathbb{F} \mathcal{G} \mathcal{E} = \mathbb{P} \mathbb{D}^2 \mathbb{P}^{-1} \mathcal{E}, \quad (31) $$

where $\mathbb{P}$ and $\mathbb{D}^2$ contain, respectively, the eigenvectors and eigenvalues of the matrix $\mathbb{F}\mathcal{G}$, each having dimensions of $2(2N + 1) \times 2(2N + 1)$ and take the following form:

$$ \mathbb{P} = \begin{pmatrix} \mathbb{P}^{(11)} & \mathbb{P}^{(12)} \\ \mathbb{P}^{(21)} & \mathbb{P}^{(22)} \end{pmatrix}, \quad \mathbb{D} = \begin{pmatrix} \mathbb{D}^{(11)} & \mathbb{D} \\ \mathbb{D} & \mathbb{D}^{(22)} \end{pmatrix}. \quad (32) $$

The solution of the set of equations (28) can then be written:

$$ \left\{ \begin{array}{l} \mathcal{E} = \mathbb{P}(e^{\mathbb{D}z} \mathcal{A} + e^{-\mathbb{D}z} \mathcal{B}) \\ \tilde{\mathcal{H}} = \mathbb{P}'(e^{\mathbb{D}z} \mathcal{A} - e^{-\mathbb{D}z} \mathcal{B}) \end{array} \right. \quad (33) $$

where $\mathcal{A} = [A_{xn}, A_{yn}]^t$ and $\mathcal{B} = [B_{xn}, B_{yn}]^t$ ($n \in [-N, N]$) are arbitrary constant vectors, and $\mathbb{P}' = \mathbb{F}^{-1}\mathbb{P}\mathbb{D}$.

## B. Boundary conditions

After expressing the fields in all zones, we are going to enforce the boundary conditions at their interfaces. For both analytical and numerical convenience, we introduce an additional phase factor in the expressions of the fields within zones 3 and 4. Specifically in zone 3, we take the $z$ dependence of the fields under the form $e^{ik^{(i)}\phi(z-h_1)}$, while in zone 4, we take it under the form $e^{ik^{(i)}\phi(z-h_1-h_2)}$. These adjustments facilitate the calculations and can be easily reverted at the conclusion of the analysis. For this reason, we will first obtain the scattering matrix with this simplified phase $\tilde{\mathcal{S}}$, and then multiply it by the accumulated phase to obtain the actual scattering matrix $\mathcal{S}$ (see Eq. (8)).

Using the zero thickness model for the graphene sheet and applying the boundary conditions for the electric field components at the interface ($z=0$) between zones 1 and 2, in the cartesian coordinates system, we obtain (after projection on the Fourier basis):
$$ \left( \begin{matrix} -\frac{k_y}{k_n}(I_{1n} + R_{1n}) + \frac{c}{\sqrt{\varepsilon_1 \omega}} k_{zn}^{(1)} \frac{k_{xn}}{k_n} (I_{2n} - R_{2n}) \\ \frac{k_{xn}}{k_n}(I_{1n} + R_{1n}) + \frac{c}{\sqrt{\varepsilon_1 \omega}} k_{zn}^{(1)} \frac{k_y}{k_n} (I_{2n} - R_{2n}) \end{matrix} \right) = \left( \begin{matrix} \mathbb{P}_{nm}^{(11)}(A_{xm} + B_{xm}) + \mathbb{P}_{nm}^{(12)}(A_{ym} + B_{ym}) \\ \mathbb{P}_{nm}^{(21)}(A_{xm} + B_{xm}) + \mathbb{P}_{nm}^{(22)}(A_{ym} + B_{ym}) \end{matrix} \right). \quad (34) $$

Doing the same for the magnetic field components, we have:

$$ \left( \begin{matrix} \mathbb{P}_{nm}^{(11)}(A_{xm} - B_{xm}) + \mathbb{P}_{nm}^{(12)}(A_{ym} - B_{ym}) + \frac{c}{\omega} k_{zn}^{(1)} \frac{k_{xn}}{k_n} (I_{1n} - R_{1n}) + \sqrt{\varepsilon_1} \frac{k_y}{k_n} (I_{2n} + R_{2n}) \\ \mathbb{P}_{nm}^{(21)}(A_{xm} - B_{xm}) + \mathbb{P}_{nm}^{(22)}(A_{ym} - B_{ym}) + \frac{c}{\omega} k_{zn}^{(1)} \frac{k_y}{k_n} (I_{1n} - R_{1n}) - \sqrt{\varepsilon_1} \frac{k_{xn}}{k_n} (I_{2n} + R_{2n}) \end{matrix} \right) = \sigma_g Z_0 \left( \begin{matrix} \frac{k_{xn}}{k_n} (I_{1n} + R_{1n}) + \frac{c}{\sqrt{\varepsilon_1 \omega}} k_{zn}^{(1)} \frac{k_y}{k_n} (I_{2n} - R_{2n}) \\ \frac{k_y}{k_n} (I_{1n} + R_{1n}) - \frac{c}{\sqrt{\varepsilon_1 \omega}} k_{zn}^{(1)} \frac{k_{xn}}{k_n} (I_{2n} - R_{2n}) \end{matrix} \right). \quad (35) $$

At the interface $z = h_1$ we have the following equations for the electric field

$$ \left( \begin{matrix} -\frac{k_y}{k_n}(C_{1n} + C'_{1n}) + \frac{c}{\sqrt{\varepsilon_3 \omega}} k_{zn}^{(3)} \frac{k_{xn}}{k_n} (C_{2n} - C'_{2n}) \\ \frac{k_{xn}}{k_n}(C_{1n} + C'_{1n}) + \frac{c}{\sqrt{\varepsilon_3 \omega}} k_{zn}^{(3)} \frac{k_y}{k_n} (C_{2n} - C'_{2n}) \end{matrix} \right) = \left( \begin{matrix} \mathbb{P}_{nm}^{(11)}(e^{\mathbb{D}_{mm}^{(11)}h_1} A_{xm} + e^{-\mathbb{D}_{mm}^{(11)}h_1} B_{xm}) + \mathbb{P}_{nm}^{(12)}(e^{\mathbb{D}_{mm}^{(22)}h_1} A_{ym} + e^{-\mathbb{D}_{mm}^{(22)}h_1} B_{ym}) \\ \mathbb{P}_{nm}^{(21)}(e^{\mathbb{D}_{mm}^{(11)}h_1} A_{x,m} + e^{-\mathbb{D}_{mm}^{(11)}h_1} B_{xm}) + \mathbb{P}_{nm}^{(22)}(e^{\mathbb{D}_{mm}^{(22)}h_1} A_{ym} + e^{-\mathbb{D}_{mm}^{(22)}h_1} B_{ym}) \end{matrix} \right), \quad (36) $$

and the boundary conditions for the magnetic field are as follows:

$$ \left( \begin{matrix} -\frac{c}{\omega} k_{zn}^{(3)} \frac{k_{xn}}{k_n} (C_{1n} - C'_{1n}) - \sqrt{\varepsilon_3} \frac{k_y}{k_n} (C_{2n} + C'_{2n}) \\ -\frac{c}{\omega} k_{zn}^{(3)} \frac{k_y}{k_n} (C_{1n} - C'_{1n}) + \sqrt{\varepsilon_3} \frac{k_{xn}}{k_n} (C_{2n} + C'_{2n}) \end{matrix} \right) = \left( \begin{matrix} \mathbb{P}_{nm}^{(11)}(e^{\mathbb{D}_{mm}^{(11)}h_1} A_{xm} - e^{-\mathbb{D}_{mm}^{(11)}h_1} B_{xm}) + \mathbb{P}_{nm}^{(12)}(e^{\mathbb{D}_{mm}^{(22)}h_1} A_{ym} - e^{-\mathbb{D}_{mm}^{(22)}h_1} B_{ym}) \\ \mathbb{P}_{nm}^{(21)}(e^{\mathbb{D}_{mm}^{(11)}h_1} A_{xm} - e^{-\mathbb{D}_{mm}^{(11)}h_1} B_{xm}) + \mathbb{P}_{nm}^{(22)}(e^{\mathbb{D}_{mm}^{(22)}h_1} A_{ym} - e^{-\mathbb{D}_{mm}^{(22)}h_1} B_{ym}) \end{matrix} \right). \quad (37) $$

Likewise, at the final interface, $z = h_1 + h_2$, we obtain

$$ \left( \begin{matrix} -\frac{k_y}{k_n} \left( e^{ik_{zn}^{(3)}h_2} C_{1n} + e^{-ik_{zn}^{(3)}h_2} C'_{1n} \right) + \frac{c}{\sqrt{\varepsilon_3 \omega}} k_{zn}^{(3)} \frac{k_{xn}}{k_n} \left( e^{ik_{zn}^{(3)}h_2} C_{2n} - e^{-ik_{zn}^{(3)}h_2} C'_{2n} \right) \\ \frac{k_{xn}}{k_n} \left( e^{ik_{zn}^{(3)}h_2} C_{1n} + e^{-ik_{zn}^{(3)}h_2} C'_{1n} \right) + \frac{c}{\sqrt{\varepsilon_3 \omega}} k_{zn}^{(3)} \frac{k_y}{k_n} \left( e^{ik_{zn}^{(3)}h_2} C_{2n} - e^{-ik_{zn}^{(3)}h_2} C'_{2n} \right) \end{matrix} \right) = \left( \begin{matrix} -\frac{k_y}{k_n} (T_{1n} + J_{1n}) + \frac{c}{\sqrt{\varepsilon_4 \omega}} k_{zn}^{(4)} \frac{k_{xn}}{k_n} (T_{2n} - J_{2n}) \\ \frac{k_{xn}}{k_n} (T_{1n} + J_{1n}) + \frac{c}{\sqrt{\varepsilon_4 \omega}} k_{zn}^{(4)} \frac{k_y}{k_n} (T_{2n} - J_{2n}) \end{matrix} \right) \quad (38) $$

and

$$ \left( \begin{matrix} -\frac{c}{\omega} k_{zn}^{(3)} \frac{k_{xn}}{k_n} \left( e^{ik_{zn}^{(3)}h_2} C_{1n} - e^{-ik_{zn}^{(3)}h_2} C'_{1n} \right) - \sqrt{\varepsilon_3} \frac{k_y}{k_n} \left( e^{ik_{zn}^{(3)}h_2} C_{2n} + e^{-ik_{zn}^{(3)}h_2} C'_{2n} \right) \\ -\frac{c}{\omega} k_{zn}^{(3)} \frac{k_y}{k_n} \left( e^{ik_{zn}^{(3)}h_2} C_{1n} - e^{-ik_{zn}^{(3)}h_2} C'_{1n} \right) + \sqrt{\varepsilon_3} \frac{k_{xn}}{k_n} \left( e^{ik_{zn}^{(3)}h_2} C_{2n} + e^{-ik_{zn}^{(3)}h_2} C'_{2n} \right) \end{matrix} \right) = \left( \begin{matrix} -\frac{c}{\omega} k_{zn}^{(4)} \frac{k_{xn}}{k_n} (T_{1n} - J_{1n}) - \sqrt{\varepsilon_4} \frac{k_y}{k_n} (T_{2n} + J_{2n}) \\ -\frac{c}{\omega} k_{zn}^{(4)} \frac{k_y}{k_n} (T_{1n} - J_{1n}) + \sqrt{\varepsilon_4} \frac{k_{xn}}{k_n} (T_{2n} + J_{2n}) \end{matrix} \right). \quad (39) $$

### C. Scattering matrix

$$ \left( \begin{matrix} \mathcal{R} \\ \mathcal{A} \end{matrix} \right) = \tilde{\mathcal{S}}_1 \left( \begin{matrix} \mathcal{I} \\ \mathcal{B} \end{matrix} \right), \quad \left( \begin{matrix} \mathcal{B} \\ \mathcal{C} \end{matrix} \right) = \tilde{\mathcal{S}}_2 \left( \begin{matrix} \mathcal{A} \\ \mathcal{C}' \end{matrix} \right), \\ \text{and} \quad \left( \begin{matrix} \mathcal{C}' \\ \mathcal{T} \end{matrix} \right) = \tilde{\mathcal{S}}_3 \left( \begin{matrix} \mathcal{C} \\ \mathcal{J} \end{matrix} \right), \quad (40) $$

Finally, we can now write the previous equations in matrix compact form as

where the explicit expressions of $\tilde{\mathcal{S}}_1$, $\tilde{\mathcal{S}}_2$, and $\tilde{\mathcal{S}}_3$ are already given in section III.
$$ \left( \begin{matrix} -\frac{k_y}{k_n} (I_{1n} + R_{1n}) + \frac{c}{\sqrt{\epsilon_1 \omega}} k_{zn}^{(1)} \frac{k_{xn}}{k_n} (I_{2n} - R_{2n}) \\ \frac{k_{xn}}{k_n} (I_{1n} + R_{1n}) + \frac{c}{\sqrt{\epsilon_1 \omega}} k_{zn}^{(1)} \frac{k_y}{k_n} (I_{2n} - R_{2n}) \end{matrix} \right) = \left( \begin{matrix} \mathbb{P}_{nm}^{(11)} (A_{xm} + B_{xm}) + \mathbb{P}_{nm}^{(12)} (A_{ym} + B_{ym}) \\ \mathbb{P}_{nm}^{(21)} (A_{xm} + B_{xm}) + \mathbb{P}_{nm}^{(22)} (A_{ym} + B_{ym}) \end{matrix} \right). \quad (34) $$

Doing the same for the magnetic field components, we have:

$$ \left( \begin{matrix} \mathbb{P}_{nm}^{(11)} (A_{xm} - B_{xm}) + \mathbb{P}_{nm}^{(12)} (A_{ym} - B_{ym}) + \frac{c}{\omega} k_{zn}^{(1)} \frac{k_{xn}}{k_n} (I_{1n} - R_{1n}) + \sqrt{\epsilon_1} \frac{k_y}{k_n} (I_{2n} + R_{2n}) \\ \mathbb{P}_{nm}^{(21)} (A_{xm} - B_{xm}) + \mathbb{P}_{nm}^{(22)} (A_{ym} - B_{ym}) + \frac{c}{\omega} k_{zn}^{(1)} \frac{k_y}{k_n} (I_{1n} - R_{1n}) - \sqrt{\epsilon_1} \frac{k_{xn}}{k_n} (I_{2n} + R_{2n}) \end{matrix} \right) = \sigma_g Z_0 \left( \begin{matrix} \frac{k_{xn}}{k_n} (I_{1n} + R_{1n}) + \frac{c}{\sqrt{\epsilon_1 \omega}} k_{zn}^{(1)} \frac{k_y}{k_n} (I_{2n} - R_{2n}) \\ \frac{k_y}{k_n} (I_{1n} + R_{1n}) - \frac{c}{\sqrt{\epsilon_1 \omega}} k_{zn}^{(1)} \frac{k_{xn}}{k_n} (I_{2n} - R_{2n}) \end{matrix} \right). \quad (35) $$

At the interface $z = h_1$ we have the following equations for the electric field

$$ \left( \begin{matrix} -\frac{k_y}{k_n} (C_{1n} + C'_{1n}) + \frac{c}{\sqrt{\epsilon_3 \omega}} k_{zn}^{(3)} \frac{k_{xn}}{k_n} (C_{2n} - C'_{2n}) \\ \frac{k_{xn}}{k_n} (C_{1n} + C'_{1n}) + \frac{c}{\sqrt{\epsilon_3 \omega}} k_{zn}^{(3)} \frac{k_y}{k_n} (C_{2n} - C'_{2n}) \end{matrix} \right) = \left( \begin{matrix} \mathbb{P}_{nm}^{(11)} (e^{\mathbb{D}_{mm}^{(11)}h_1} A_{xm} + e^{-\mathbb{D}_{mm}^{(11)}h_1} B_{xm}) + \mathbb{P}_{nm}^{(12)} (e^{\mathbb{D}_{mm}^{(22)}h_1} A_{ym} + e^{-\mathbb{D}_{mm}^{(22)}h_1} B_{ym}) \\ \mathbb{P}_{nm}^{(21)} (e^{\mathbb{D}_{mm}^{(11)}h_1} A_{x,m} + e^{-\mathbb{D}_{mm}^{(11)}h_1} B_{xm}) + \mathbb{P}_{nm}^{(22)} (e^{\mathbb{D}_{mm}^{(22)}h_1} A_{ym} + e^{-\mathbb{D}_{mm}^{(22)}h_1} B_{ym}) \end{matrix} \right), \quad (36) $$

and the boundary conditions for the magnetic field are as follows:

$$ \left( \begin{matrix} -\frac{c}{\omega} k_{zn}^{(3)} \frac{k_{xn}}{k_n} (C_{1n} - C'_{1n}) - \sqrt{\epsilon_3} \frac{k_y}{k_n} (C_{2n} + C'_{2n}) \\ -\frac{c}{\omega} k_{zn}^{(3)} \frac{k_y}{k_n} (C_{1n} - C'_{1n}) + \sqrt{\epsilon_3} \frac{k_{xn}}{k_n} (C_{2n} + C'_{2n}) \end{matrix} \right) = \left( \begin{matrix} \mathbb{P}_{nm}^{(11)} (e^{\mathbb{D}_{mm}^{(11)}h_1} A_{xm} - e^{-\mathbb{D}_{mm}^{(11)}h_1} B_{xm}) + \mathbb{P}_{nm}^{(12)} (e^{\mathbb{D}_{mm}^{(22)}h_1} A_{ym} - e^{-\mathbb{D}_{mm}^{(22)}h_1} B_{ym}) \\ \mathbb{P}_{nm}^{(21)} (e^{\mathbb{D}_{mm}^{(11)}h_1} A_{xm} - e^{-\mathbb{D}_{mm}^{(11)}h_1} B_{xm}) + \mathbb{P}_{nm}^{(22)} (e^{\mathbb{D}_{mm}^{(22)}h_1} A_{ym} - e^{-\mathbb{D}_{mm}^{(22)}h_1} B_{ym}) \end{matrix} \right). \quad (37) $$

Likewise, at the final interface, $z = h_1 + h_2$, we obtain

$$ \left( \begin{matrix} -\frac{k_y}{k_n} \left( e^{ik_{zn}^{(3)}h_2} C_{1n} + e^{-ik_{zn}^{(3)}h_2} C'_{1n} \right) + \frac{c}{\sqrt{\epsilon_3 \omega}} k_{zn}^{(3)} \frac{k_{xn}}{k_n} \left( e^{ik_{zn}^{(3)}h_2} C_{2n} - e^{-ik_{zn}^{(3)}h_2} C'_{2n} \right) \\ \frac{k_{xn}}{k_n} \left( e^{ik_{zn}^{(3)}h_2} C_{1n} + e^{-ik_{zn}^{(3)}h_2} C'_{1n} \right) + \frac{c}{\sqrt{\epsilon_3 \omega}} k_{zn}^{(3)} \frac{k_y}{k_n} \left( e^{ik_{zn}^{(3)}h_2} C_{2n} - e^{-ik_{zn}^{(3)}h_2} C'_{2n} \right) \end{matrix} \right) = \left( \begin{matrix} -\frac{k_y}{k_n} (T_{1n} + J_{1n}) + \frac{c}{\sqrt{\epsilon_4 \omega}} k_{zn}^{(4)} \frac{k_{xn}}{k_n} (T_{2n} - J_{2n}) \\ \frac{k_{xn}}{k_n} (T_{1n} + J_{1n}) + \frac{c}{\sqrt{\epsilon_4 \omega}} k_{zn}^{(4)} \frac{k_y}{k_n} (T_{2n} - J_{2n}) \end{matrix} \right) \quad (38) $$

and

$$ \left( \begin{matrix} -\frac{c}{\omega} k_{zn}^{(3)} \frac{k_{xn}}{k_n} \left( e^{ik_{zn}^{(3)}h_2} C_{1n} - e^{-ik_{zn}^{(3)}h_2} C'_{1n} \right) - \sqrt{\epsilon_3} \frac{k_y}{k_n} \left( e^{ik_{zn}^{(3)}h_2} C_{2n} + e^{-ik_{zn}^{(3)}h_2} C'_{2n} \right) \\ -\frac{c}{\omega} k_{zn}^{(3)} \frac{k_y}{k_n} \left( e^{ik_{zn}^{(3)}h_2} C_{1n} - e^{-ik_{zn}^{(3)}h_2} C'_{1n} \right) + \sqrt{\epsilon_3} \frac{k_{xn}}{k_n} \left( e^{ik_{zn}^{(3)}h_2} C_{2n} + e^{-ik_{zn}^{(3)}h_2} C'_{2n} \right) \end{matrix} \right) = \left( \begin{matrix} -\frac{c}{\omega} k_{zn}^{(4)} \frac{k_{xn}}{k_n} (T_{1n} - J_{1n}) - \sqrt{\epsilon_4} \frac{k_y}{k_n} (T_{2n} + J_{2n}) \\ -\frac{c}{\omega} k_{zn}^{(4)} \frac{k_y}{k_n} (T_{1n} - J_{1n}) + \sqrt{\epsilon_4} \frac{k_{xn}}{k_n} (T_{2n} + J_{2n}) \end{matrix} \right). \quad (39) $$

### C. Scattering matrix

$$ \left( \begin{matrix} \mathcal{R} \\ \mathcal{A} \end{matrix} \right) = \tilde{\mathcal{S}}_1 \left( \begin{matrix} \mathcal{I} \\ \mathcal{B} \end{matrix} \right), \quad \left( \begin{matrix} \mathcal{B} \\ \mathcal{C} \end{matrix} \right) = \tilde{\mathcal{S}}_2 \left( \begin{matrix} \mathcal{A} \\ \mathcal{C}' \end{matrix} \right), \\ \text{and} \quad \left( \begin{matrix} \mathcal{C}' \\ \mathcal{T} \end{matrix} \right) = \tilde{\mathcal{S}}_3 \left( \begin{matrix} \mathcal{C} \\ \mathcal{J} \end{matrix} \right), \quad (40) $$

Finally, we can now write the previous equations in matrix compact form as

where the explicit expressions of $\tilde{\mathcal{S}}_1$, $\tilde{\mathcal{S}}_2$, and $\tilde{\mathcal{S}}_3$ are already given in section III.
Furthermore, to derive the actual scattering matrix, we

incorporate a phase correction, as previously mentioned.
This adjustment is implemented in Eq. (8).

[1] H. B. G. Casimir, On the attraction between two perfectly conducting plates, Indag. Math. **10**, 261 (1948).

[2] E. L. I.E. Dzyaloshinskii and L. Pitaevskii, The general theory of van der waals forces, *Advances in Physics* **10**, 165 (1961).

[3] D. Dalvit, P. Milonni, D. Roberts, and F. da Rosa, *Casimir Physics*, Lecture Notes in Physics (Springer Berlin Heidelberg, 2011).

[4] S. K. Lamoreaux, Demonstration of the casimir force in the 0.6 to 6µm range, *Phys. Rev. Lett.* **78**, 5 (1997).

[5] G. Bressi, G. Carugno, R. Onofrio, and G. Ruoso, Measurement of the casimir force between parallel metallic surfaces, *Phys. Rev. Lett.* **88**, 041804 (2002).

[6] M. Antezza, H. B. Chan, B. Guizal, V. N. Marachevsky, R. Messina, and M. Wang, Giant casimir torque between rotated gratings and the $\theta = 0$ anomaly, *Phys. Rev. Lett.* **124**, 013903 (2020).

[7] R. S. Decca, D. López, E. Fischbach, and D. E. Krause, Measurement of the casimir force between dissimilar metals, *Phys. Rev. Lett.* **91**, 050402 (2003).

[8] A. Canaguier-Durand, P. A. Maia Neto, I. Cavero-Pelaez, A. Lambrecht, and S. Reynaud, Casimir interaction between plane and spherical metallic surfaces, *Phys. Rev. Lett.* **102**, 230404 (2009).

[9] R. O. Nunes, B. Spreng, R. de Melo e Souza, G.-L. Ingold, P. A. Maia Neto, and F. S. S. Rosa, The casimir interaction between spheres immersed in electrolytes, *Universe* **7**, 156 (2021).

[10] T. Schoger, B. Spreng, G.-L. Ingold, P. A. Maia Neto, and S. Reynaud, Universal casimir interaction between two dielectric spheres in salted water, *Phys. Rev. Lett.* **128**, 230602 (2022).

[11] G. Bimonte, Beyond-proximity-force-approximation casimir force between two spheres at finite temperature, *Phys. Rev. D* **97**, 085011 (2018).

[12] A. Lambrecht and V. N. Marachevsky, Casimir interaction of dielectric gratings, *Phys. Rev. Lett.* **101**, 160403 (2008).

[13] A. Noto, R. Messina, B. Guizal, and M. Antezza, Casimir-lifshitz force out of thermal equilibrium between dielectric gratings, *Phys. Rev. A* **90**, 022120 (2014).

[14] H. B. Chan, Y. Bao, J. Zou, R. A. Cirelli, F. Klemens, W. M. Mansfield, and C. S. Pai, Measurement of the casimir force between a gold sphere and a silicon surface with nanoscale trench arrays, *Phys. Rev. Lett.* **101**, 030401 (2008).

[15] M. Wang, L. Tang, C. Y. Ng, R. Messina, B. Guizal, J. A. Crosse, M. Antezza, C. T. Chan, and H. B. Chan, Strong geometry dependence of the casimir force between interpenetrated rectangular gratings, *Nature Communications* **12**, 600 (2021).

[16] R. Messina, A. Noto, B. Guizal, and M. Antezza, Radiative heat transfer between metallic gratings using fourier modal method with adaptive spatial resolution, *Phys. Rev. B* **95**, 125404 (2017).

[17] R. Messina, P. A. Maia Neto, B. Guizal, and M. Antezza, Casimir interaction between a sphere and a grating, *Phys. Rev. A* **92**, 062504 (2015).

[18] Y. Jeyar, K. Austry, M. Luo, B. Guizal, H. B. Chan, and M. Antezza, Casimir-lifshitz force between graphene-based structures out of thermal equilibrium, *Phys. Rev. B* **108**, 115412 (2023).

[19] M. Antezza, N. Emelianova, and N. Khusnutdinov, The normal Casimir-Lifshitz force for laterally moving graphene, *Nanotechnol.* **35**, 235001 (2024), arXiv:2303.03115.

[20] J.-S. Wang and M. Antezza, Photon mediated transport of energy, linear momentum, and angular momentum in fullerene and graphene systems beyond local equilibrium, *Phys. Rev. B* **109**, 125105 (2024).

[21] V. B. Svetovoy, P. J. van Zwol, and J. Chevrier, Plasmon enhanced near-field radiative heat transfer for graphene covered dielectrics, *Phys. Rev. B* **85**, 155418 (2012).

[22] O. Ilic, M. Jablan, J. D. Joannopoulos, I. Celanovic, H. Buljan, and M. Soljačić, Near-field thermal radiation transfer controlled by plasmons in graphene, *Phys. Rev. B* **85**, 155422 (2012).

[23] Z. H. Zheng, X. L. Liu, A. Wang, and Y. M. Xuan, Graphene-assisted near-field radiative thermal rectifier based on phase transition of vanadium dioxide (VO₂), *Int. J. Heat Mass Transf.* **109**, 63 (2017).

[24] A. Volokitin, Casimir friction and near-field radiative heat transfer in graphene structures, *Z. Naturforsch. A* **72**, 171 (2017).

[25] M. Liu, Y. Zhang, G. L. Klimchitskaya, V. M. Mostepanenko, and U. Mohideen, Demonstration of an unusual thermal effect in the casimir force from graphene, *Phys. Rev. Lett.* **126**, 206802 (2021).

[26] C. Abbas, B. Guizal, and M. Antezza, Strong thermal and electrostatic manipulation of the casimir force in graphene multilayers, *Phys. Rev. Lett.* **118**, 126101 (2017).

[27] M. Liu, Y. Zhang, G. L. Klimchitskaya, V. M. Mostepanenko, and U. Mohideen, Experimental and theoretical investigation of the thermal effect in the casimir interaction from graphene, *Phys. Rev. B* **104**, 085436 (2021).

[28] P. Rodriguez-Lopez, D.-N. Le, M. J. Calderón, E. Bascones, and L. M. Woods, Twisted bilayered graphenes at magic angles and casimir interactions: correlation-driven effects, *2D Mater.* **10**, 014006 (2022).

[29] L. A. Falkovsky, Optical properties of graphene, *Journal of Physics: Conference Series* **129**, 012004 (2008).

[30] L. A. Falkovsky and A. A. Varlamov, Space-time dispersion of graphene conductivity, *The European Physical Journal B* **56**, 281 (2007).

[31] S. A. Awan, A. Lombardo, A. Colli, G. Privitera, T. S. Kulmala, J. M. Kivioja, M. Koshino, and A. C. Ferrari, Transport conductivity of graphene at rf and microwave frequencies, *2D Materials* **3**, 015010 (2016).

[32] T. Zhu, M. Antezza, and J.-S. Wang, Dynamical polarizability of graphene with spatial dispersion, *Phys. Rev. B* **103**, 125421 (2021).

[33] P. Rodriguez-Lopez and M. Antezza, Graphene conductivity: Kubo model versus QFT-based model, arXiv
Furthermore, to derive the actual scattering matrix, we

incorporate a phase correction, as previously mentioned.
This adjustment is implemented in Eq. (8).

[1] H. B. G. Casimir, On the attraction between two perfectly conducting plates, Indag. Math. **10**, 261 (1948).

[2] E. L. I.E. Dzyaloshinskii and L. Pitaevskii, The general theory of van der waals forces, *Advances in Physics* **10**, 165 (1961).

[3] D. Dalvit, P. Milonni, D. Roberts, and F. da Rosa, *Casimir Physics*, Lecture Notes in Physics (Springer Berlin Heidelberg, 2011).

[4] S. K. Lamoreaux, Demonstration of the casimir force in the 0.6 to 6µm range, *Phys. Rev. Lett.* **78**, 5 (1997).

[5] G. Bressi, G. Carugno, R. Onofrio, and G. Ruoso, Measurement of the casimir force between parallel metallic surfaces, *Phys. Rev. Lett.* **88**, 041804 (2002).

[6] M. Antezza, H. B. Chan, B. Guizal, V. N. Marachevsky, R. Messina, and M. Wang, Giant casimir torque between rotated gratings and the $\theta = 0$ anomaly, *Phys. Rev. Lett.* **124**, 013903 (2020).

[7] R. S. Decca, D. López, E. Fischbach, and D. E. Krause, Measurement of the casimir force between dissimilar metals, *Phys. Rev. Lett.* **91**, 050402 (2003).

[8] A. Canaguier-Durand, P. A. Maia Neto, I. Cavero-Pelaez, A. Lambrecht, and S. Reynaud, Casimir interaction between plane and spherical metallic surfaces, *Phys. Rev. Lett.* **102**, 230404 (2009).

[9] R. O. Nunes, B. Spreng, R. de Melo e Souza, G.-L. Ingold, P. A. Maia Neto, and F. S. S. Rosa, The casimir interaction between spheres immersed in electrolytes, *Universe* **7**, 156 (2021).

[10] T. Schoger, B. Spreng, G.-L. Ingold, P. A. Maia Neto, and S. Reynaud, Universal casimir interaction between two dielectric spheres in salted water, *Phys. Rev. Lett.* **128**, 230602 (2022).

[11] G. Bimonte, Beyond-proximity-force-approximation casimir force between two spheres at finite temperature, *Phys. Rev. D* **97**, 085011 (2018).

[12] A. Lambrecht and V. N. Marachevsky, Casimir interaction of dielectric gratings, *Phys. Rev. Lett.* **101**, 160403 (2008).

[13] A. Noto, R. Messina, B. Guizal, and M. Antezza, Casimir-lifshitz force out of thermal equilibrium between dielectric gratings, *Phys. Rev. A* **90**, 022120 (2014).

[14] H. B. Chan, Y. Bao, J. Zou, R. A. Cirelli, F. Klemens, W. M. Mansfield, and C. S. Pai, Measurement of the casimir force between a gold sphere and a silicon surface with nanoscale trench arrays, *Phys. Rev. Lett.* **101**, 030401 (2008).

[15] M. Wang, L. Tang, C. Y. Ng, R. Messina, B. Guizal, J. A. Crosse, M. Antezza, C. T. Chan, and H. B. Chan, Strong geometry dependence of the casimir force between interpenetrated rectangular gratings, *Nature Communications* **12**, 600 (2021).

[16] R. Messina, A. Noto, B. Guizal, and M. Antezza, Radiative heat transfer between metallic gratings using fourier modal method with adaptive spatial resolution, *Phys. Rev. B* **95**, 125404 (2017).

[17] R. Messina, P. A. Maia Neto, B. Guizal, and M. Antezza, Casimir interaction between a sphere and a grating, *Phys. Rev. A* **92**, 062504 (2015).

[18] Y. Jeyar, K. Austry, M. Luo, B. Guizal, H. B. Chan, and M. Antezza, Casimir-lifshitz force between graphene-based structures out of thermal equilibrium, *Phys. Rev. B* **108**, 115412 (2023).

[19] M. Antezza, N. Emelianova, and N. Khusnutdinov, The normal Casimir-Lifshitz force for laterally moving graphene, *Nanotechnol.* **35**, 235001 (2024), arXiv:2303.03115.

[20] J.-S. Wang and M. Antezza, Photon mediated transport of energy, linear momentum, and angular momentum in fullerene and graphene systems beyond local equilibrium, *Phys. Rev. B* **109**, 125105 (2024).

[21] V. B. Svetovoy, P. J. van Zwol, and J. Chevrier, Plasmon enhanced near-field radiative heat transfer for graphene covered dielectrics, *Phys. Rev. B* **85**, 155418 (2012).

[22] O. Ilic, M. Jablan, J. D. Joannopoulos, I. Celanovic, H. Buljan, and M. Soljačić, Near-field thermal radiation transfer controlled by plasmons in graphene, *Phys. Rev. B* **85**, 155422 (2012).

[23] Z. H. Zheng, X. L. Liu, A. Wang, and Y. M. Xuan, Graphene-assisted near-field radiative thermal rectifier based on phase transition of vanadium dioxide (VO₂), *Int. J. Heat Mass Transf.* **109**, 63 (2017).

[24] A. Volokitin, Casimir friction and near-field radiative heat transfer in graphene structures, *Z. Naturforsch. A* **72**, 171 (2017).

[25] M. Liu, Y. Zhang, G. L. Klimchitskaya, V. M. Mostepanenko, and U. Mohideen, Demonstration of an unusual thermal effect in the casimir force from graphene, *Phys. Rev. Lett.* **126**, 206802 (2021).

[26] C. Abbas, B. Guizal, and M. Antezza, Strong thermal and electrostatic manipulation of the casimir force in graphene multilayers, *Phys. Rev. Lett.* **118**, 126101 (2017).

[27] M. Liu, Y. Zhang, G. L. Klimchitskaya, V. M. Mostepanenko, and U. Mohideen, Experimental and theoretical investigation of the thermal effect in the casimir interaction from graphene, *Phys. Rev. B* **104**, 085436 (2021).

[28] P. Rodriguez-Lopez, D.-N. Le, M. J. Calderón, E. Bascones, and L. M. Woods, Twisted bilayered graphenes at magic angles and casimir interactions: correlation-driven effects, *2D Mater.* **10**, 014006 (2022).

[29] L. A. Falkovsky, Optical properties of graphene, *Journal of Physics: Conference Series* **129**, 012004 (2008).

[30] L. A. Falkovsky and A. A. Varlamov, Space-time dispersion of graphene conductivity, *The European Physical Journal B* **56**, 281 (2007).

[31] S. A. Awan, A. Lombardo, A. Colli, G. Privitera, T. S. Kulmala, J. M. Kivioja, M. Koshino, and A. C. Ferrari, Transport conductivity of graphene at rf and microwave frequencies, *2D Materials* **3**, 015010 (2016).

[32] T. Zhu, M. Antezza, and J.-S. Wang, Dynamical polarizability of graphene with spatial dispersion, *Phys. Rev. B* **103**, 125421 (2021).

[33] P. Rodriguez-Lopez and M. Antezza, Graphene conductivity: Kubo model versus QFT-based model, *arXiv*
[34] Y. Jeyar, M. Luo, K. Austry, B. Guizal, Y. Zheng, H. B. Chan, and M. Antezza, Tunable nonadditivity in the casimir-lifshitz force between graphene gratings, Phys. Rev. A 108, 062811 (2023).

[35] Y. Jeyar, M. Antezza, and B. Guizal, Electromagnetic scattering by a partially graphene-coated dielectric cylinder: Efficient computation and multiple plasmonic resonances, Phys. Rev. E 107, 025306 (2023).

[36] M. Luo, Y. Jeyar, B. Guizal, J. Zhao, and M. Antezza, Effect of graphene grating coating on near-field radiative heat transfer, Applied Physics Letters 123, 253902 (2023).

[37] M. Luo, Y. Jeyar, B. Guizal, and M. Antezza, Near-field radiative heat transfer between shifted graphene gratings (2024), arXiv:2401.14357.

[38] S. Ogawa, M. Iwakawa, M. Shimatani, and S. Fukushima, Hybrid graphene-high-aspect ratio plasmonic nanograting systems, Journal of Applied Physics 135, 153107 (2024).

[39] R. Messina and M. Antezza, Scattering-matrix approach to casimir-lifshitz force and heat transfer out of thermal equilibrium between arbitrary bodies, Phys. Rev. A 84, 042102 (2011).

[40] W. J. Tropf and M. E. Thomas, Handbook of Optical Constants of Solids (edited by E. Palik New York, 1998).

[41] G. Granet and B. Guizal, Efficient implementation of the coupled-wave method for metallic lamellar gratings in tm polarization, J. Opt. Soc. Am. A 13, 1019 (1996).
preprint arXiv:2403.02279 (2024).

[34] Y. Jeyar, M. Luo, K. Austry, B. Guizal, Y. Zheng, H. B. Chan, and M. Antezza, Tunable nonadditivity in the casimir-lifshitz force between graphene gratings, Phys. Rev. A 108, 062811 (2023).

[35] Y. Jeyar, M. Antezza, and B. Guizal, Electromagnetic scattering by a partially graphene-coated dielectric cylinder: Efficient computation and multiple plasmonic resonances, Phys. Rev. E 107, 025306 (2023).

[36] M. Luo, Y. Jeyar, B. Guizal, J. Zhao, and M. Antezza, Effect of graphene grating coating on near-field radiative heat transfer, Applied Physics Letters 123, 253902 (2023).

[37] M. Luo, Y. Jeyar, B. Guizal, and M. Antezza, Near-field radiative heat transfer between shifted graphene gratings

(2024), arXiv:2401.14357.

[38] S. Ogawa, M. Iwakawa, M. Shimatani, and S. Fukushima, Hybrid graphene-high-aspect ratio plasmonic nanograting systems, Journal of Applied Physics 135, 153107 (2024).

[39] R. Messina and M. Antezza, Scattering-matrix approach to casimir-lifshitz force and heat transfer out of thermal equilibrium between arbitrary bodies, Phys. Rev. A 84, 042102 (2011).

[40] W. J. Tropf and M. E. Thomas, Handbook of Optical Constants of Solids (edited by E. Palik New York, 1998).

[41] G. Granet and B. Guizal, Efficient implementation of the coupled-wave method for metallic lamellar gratings in tm polarization, J. Opt. Soc. Am. A 13, 1019 (1996).