Figure 1 shows the configuration of the studied heterostructure, where the hBN layer is sandwiched between two graphene sheets and the composite structure is located on SiO2-Si layers. The whole structure is illuminated by a TM-polarized beam with an incident angle of θ. There is no layer below the Si layer. The conductivity of each graphene sheet can be modeled by the following relation [93]:

$${\sigma _{1,2}}\left( {\omega ,{\mu _c},\Gamma ,T} \right)=\frac{{ - j{e^2}}}{{4\pi \hbar }}Ln\left[ {\frac{{2\left| {{\mu _{c,1,2}}} \right| - (\omega - j2\Gamma )\hbar }}{{2\left| {{\mu _{c,1,2}}} \right|+(\omega - j2\Gamma )\hbar }}} \right]+\frac{{ - j{e^2}{K_B}T}}{{\pi {\hbar ^2}(\omega - j2\Gamma )}}\left[ {\frac{{{\mu _{c,1,2}}}}{{{K_B}T}}+2Ln\left( {1+{e^{ - {{{\mu _{c,1,2}}} \mathord{\left/ {\vphantom {{{\mu _{c,1,2}}} {{K_B}T}}} \right. \kern-0pt} {{K_B}T}}}}} \right)} \right]$$

1

In (1), \(\varGamma\) is the scattering rate, \(T\) is the temperature, and \({\mu }_{c,\text{1,2}}\) is the chemical potential of each graphene. Furthermore, \(h\) is the reduced Planck’s constant, \({K}_{B}\)is Boltzmann’s constant, ω is radian frequency, and \(e\)is the electron charge in this relation.

hBN is a polar dielectric, supporting two phonon modes related to hyperbolicity, with the following permittivity tensor [74]:

$${\varepsilon _m}\left( \omega \right)={\varepsilon _{\infty ,m}}+{\varepsilon _{\infty ,m}}.\frac{{{{\left( {{\omega _{LO,m}}} \right)}^2} - {{\left( {{\omega _{TO,m}}} \right)}^2}}}{{{{\left( {{\omega _{TO,m}}} \right)}^2} - {\omega ^2} - j\omega \,{\Gamma _m}}}$$

2

In (2), \(m=\Vert or\perp\) is related to the transverse and z-axis, respectively. Moreover, \({\omega }_{LO}, {\omega }_{TO }\)show the LO and TO phonon frequencies, respectively, in which each frequency has two values in the upper and lower Reststrahlen bands:\({\omega }_{LO,\perp }=24.9 THz,{\omega }_{TO,\perp }=23.4 THz, {\omega }_{LO,\Vert }=48.3 THz, {\omega }_{TO,\Vert }=41.1 THz\). In (2), \({\varGamma }_{m}\)is a damping factor (\({\varGamma }_{\perp }=0.15 THz,{\varGamma }_{\Vert }=0.12 THz\)) and \({\varepsilon }_{m}\)is related to the high-frequency permittivity (\({\varepsilon }_{\infty ,\perp }=4.87,{\varepsilon }_{\infty ,\Vert }=2.95\)) [74]. In Fig. 2, the dielectric function of hBN permittivity is depicted, which shows the lower and upper Reststrahlen bands.

To calculate the reflection coefficient and the reflected group delay, the transfer matrix method can be utilized. For various regions, TM-polarized waves (p-polarized waves) can be written as follows:

$$\begin{gathered} {H_{y,1}}=\left( {{a_1}{e^{i{k_{z,1}}z}}+{b_1}{e^{ - i{k_{z,1}}z}}} \right){e^{i{k_x}x}}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,z<0 \hfill \\ {H_{y,2}}=\left( {{a_2}{e^{i{k_{z,2}}z}}+{b_2}{e^{ - i{k_{z,2}}z}}} \right){e^{i{k_x}x}}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,0<z<t \hfill \\ {H_{y,3}}=\left( {{a_3}{e^{i{k_{z,3}}z}}+{b_3}{e^{ - i{k_{z,3}}z}}} \right){e^{i{k_x}x}}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,t<z<t+d \hfill \\ {H_{y,4}}=\left( {{a_4}{e^{i{k_{z,4}}z}}+{b_4}{e^{ - i{k_{z,4}}z}}} \right){e^{i{k_x}x}}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,t+d<z \hfill \\ \end{gathered}$$

3

Thus, the transfer matrix of the whole structure is obtained:

$$\left[ {\begin{array}{*{20}{c}} {{a_1}} \\ {{b_1}} \end{array}} \right]=M.\left[ {\begin{array}{*{20}{c}} {{a_4}} \\ {{b_4}} \end{array}} \right]$$

4

Where

$$M={D_{1 \to 2}}.{P_2}.\,{D_{2 \to 3}}.\,{P_3}.\,{D_{3 \to 4}}$$

5

From Air to hBN, the transmission matrix is written as follows:

$${D_{1 \to 2}}=\frac{1}{2}\left( {\begin{array}{*{20}{c}} {1+{\eta _{TM}}+{\xi _{TM}}}&{1 - {\eta _{TM}} - {\xi _{TM}}} \\ {1 - {\eta _{TM}}+{\xi _{TM}}}&{1+{\eta _{TM}} - {\xi _{TM}}} \end{array}} \right)$$

6

In (6), the following parameters have been defined:

$${\eta _{TM}}=\frac{{{\varepsilon _0}\,{k_{2z}}}}{{{\varepsilon _ \bot }{k_{1z}}}}$$

7

$${\xi _{TM}}=\frac{{{\sigma _1}{k_{2z}}}}{{{\varepsilon _0}{\varepsilon _ \bot }\omega }}$$

8

Moreover, the wave number component of the incident beam in the z-direction can be obtained by (it is supposed that the direction of the incident angle is θ):

$${k_{1z}}={k_0}\cos \theta$$

9

$${k_{2z}}={k_0}\sqrt {\varepsilon _{ \bot }^{2} - \frac{{{\varepsilon _ \bot }{{\left( {\sin \theta } \right)}^2}}}{{{\varepsilon _\parallel }}}}$$

10

Now, the propagation matrix of the plasmonic wave inside the hBN layer is obtained by the following matrix (the thickness of the hBN medium is assumed to be *t*):

$${P_2}=\left( {\begin{array}{*{20}{c}} {{e^{ - j{k_{2z}}t}}}&0 \\ 0&{{e^{j{k_{2z}}t}}} \end{array}} \right)$$

11

When the beam inside the hBN medium reaches the surface of the second graphene sheet, the transmission matrix from the hBN to the SiO2 layer is expressed as:

$${D_{2 \to 3}}=\frac{1}{2}\left( {\begin{array}{*{20}{c}} {1+{{\eta ^{\prime}}_{TM}}+{{\xi ^{\prime}}_{TM}}}&{1 - {{\eta ^{\prime}}_{TM}} - {{\xi ^{\prime}}_{TM}}} \\ {1 - {{\eta ^{\prime}}_{TM}}+{{\xi ^{\prime}}_{TM}}}&{1+{{\eta ^{\prime}}_{TM}} - {{\xi ^{\prime}}_{TM}}} \end{array}} \right)$$

12

Where the following parameters have been utilized in (12):

$${\eta ^{\prime}_{TM}}=\frac{{{\varepsilon _ \bot }\,{k_{3z}}}}{{{\varepsilon _{Si{O_2}}}{k_{2z}}}}$$

13

$${\xi ^{\prime}_{TM}}=\frac{{{\sigma _2}{k_{3z}}}}{{{\varepsilon _0}{\varepsilon _{Si{O_2}}}\omega }}$$

14

Similar to the propagation matrix of the beam inside the hBN medium, the propagation matrix inside the SiO2 layer can be written as (the thickness of the SiO2 layer is assumed to be *d*):

$${P_3}=\left( {\begin{array}{*{20}{c}} {{e^{ - j{k_{3z}}d}}}&0 \\ 0&{{e^{j{k_{3z}}d}}} \end{array}} \right)$$

15

Finally, as the beam reaches the border of SiO2-Si layers, the transmission matrix can be expressed as:

$${D_{3 \to 4}}=\frac{1}{2}\left( {\begin{array}{*{20}{c}} {1+{{\eta ^{\prime\prime}}_{TM}}}&{1 - {{\eta ^{\prime\prime}}_{TM}}} \\ {1 - {{\eta ^{\prime\prime}}_{TM}}}&{1+{{\eta ^{\prime\prime}}_{TM}}} \end{array}} \right)$$

16

In (16), the following parameters have been defined:

$${\eta ^{\prime\prime}_{TM}}=\frac{{{\varepsilon _{Si{O_2}}}{k_{4z}}}}{{{\varepsilon _{Si}}{k_{3z}}}}$$

17

By calculating the elements of the transfer matrix, the reflection coefficient and the reflectance are derived by:

$$r=\frac{{{M_{21}}}}{{{M_{11}}}}=\left| {r\left( \omega \right)} \right|\exp \left( {j{\varphi _r}\left( \omega \right)} \right)$$

18

$$R={\left| {r\left( \omega \right)} \right|^2}$$

19

For the proposed structure, \({M}_{21}\)and \({M}_{11}\) are calculated:

$$\begin{gathered} {M_{11}}=\left( {1+{\eta _{TM}}+{\xi _{TM}}} \right).\left( {1+{{\eta ^{\prime}}_{TM}}+{{\xi ^{\prime}}_{TM}}} \right).\left( {1+{{\eta ^{\prime\prime}}_{TM}}} \right).{e^{ - j{k_{2z}}t}}.{e^{ - j{k_{3z}}d}}+ \hfill \\ \left( {1 - {\eta _{TM}} - {\xi _{TM}}} \right).\left( {1 - {{\eta ^{\prime}}_{TM}}+{{\xi ^{\prime}}_{TM}}} \right).\left( {1+{{\eta ^{\prime\prime}}_{TM}}} \right).{e^{+j{k_{2z}}t}}.{e^{ - j{k_{3z}}d}}+ \hfill \\ \left( {1+{\eta _{TM}}+{\xi _{TM}}} \right).\left( {1 - {{\eta ^{\prime}}_{TM}} - {{\xi ^{\prime}}_{TM}}} \right).\left( {1 - {{\eta ^{\prime\prime}}_{TM}}} \right).{e^{ - j{k_{2z}}t}}.{e^{+j{k_{3z}}d}}+ \hfill \\ \left( {1 - {\eta _{TM}} - {\xi _{TM}}} \right).\left( {1+{{\eta ^{\prime}}_{TM}} - {{\xi ^{\prime}}_{TM}}} \right).\left( {1 - {{\eta ^{\prime\prime}}_{TM}}} \right).{e^{+j{k_{2z}}t}}.{e^{+j{k_{3z}}d}} \hfill \\ \end{gathered}$$

20

$$\begin{gathered} {M_{21}}=\left( {1 - {\eta _{TM}}+{\xi _{TM}}} \right).\left( {1+{{\eta ^{\prime}}_{TM}}+{{\xi ^{\prime}}_{TM}}} \right).\left( {1+{{\eta ^{\prime\prime}}_{TM}}} \right).{e^{ - j{k_{2z}}t}}.{e^{ - j{k_{3z}}d}}+ \hfill \\ \left( {1+{\eta _{TM}} - {\xi _{TM}}} \right).\left( {1 - {{\eta ^{\prime}}_{TM}}+{{\xi ^{\prime}}_{TM}}} \right).\left( {1+{{\eta ^{\prime\prime}}_{TM}}} \right).{e^{+j{k_{2z}}t}}.{e^{ - j{k_{3z}}d}}+ \hfill \\ \left( {1 - {\eta _{TM}}+{\xi _{TM}}} \right).\left( {1 - {{\eta ^{\prime}}_{TM}} - {{\xi ^{\prime}}_{TM}}} \right).\left( {1 - {{\eta ^{\prime\prime}}_{TM}}} \right).{e^{ - j{k_{2z}}t}}.{e^{+j{k_{3z}}d}}+ \hfill \\ \left( {1+{\eta _{TM}} - {\xi _{TM}}} \right).\left( {1+{{\eta ^{\prime}}_{TM}} - {{\xi ^{\prime}}_{TM}}} \right).\left( {1 - {{\eta ^{\prime\prime}}_{TM}}} \right).{e^{+j{k_{2z}}t}}.{e^{+j{k_{3z}}d}} \hfill \\ \end{gathered}$$

21

Here, if we suppose that the incident pulse is a Gaussian beam with the central and half-width of\({ \omega }_{0},{\tau }_{0}\), respectively:

$${E_i}\left( {0,t} \right)={A_0}\exp \left( { - {{{t^2}} \mathord{\left/ {\vphantom {{{t^2}} {2\tau _{0}^{2}}}} \right. \kern-0pt} {2\tau _{0}^{2}}}} \right)\,\exp \left( { - i{\omega _0}t} \right)$$

22

It should be noted that the corresponding Fourier spectrum of (22) is:

$${E_i}\left( {0,\omega } \right)=\frac{{{A_0}\tau _{0}^{{}}}}{{2\sqrt \pi }}\exp \left( { - {{{{\left( {\omega - \omega _{0}^{2}} \right)}^2}\tau _{0}^{2}} \mathord{\left/ {\vphantom {{{{\left( {\omega - \omega _{0}^{2}} \right)}^2}\tau _{0}^{2}} 2}} \right. \kern-0pt} 2}} \right)$$

23

and thus the group delay is not a function of time (t) because all relations are written in the Fourier space (ω). Then, the reflected group delay is obtained:

$${\tau _r}={\left[ {\frac{{\partial {\varphi _r}\left( \omega \right)}}{{\partial \omega }}} \right]_{\omega ={\omega _c}}}$$

24

In (24), \({\omega }_{c}\) is the carrier frequency. Now, our model is completed for the proposed heterostructure. In what follows, we will investigate the analytical results of the above mathematical relations.