## 2 Methods
### Distributed Dielectric Constant-Based RCWA  
For a volume grating whose refractive index distributed in the volume grating can be defined by an arbitrary periodic spatial function \(n(x)\) as shown Figure 1a, a RCWA model based on a distributed refractive index model can be discretized as Figure 1b. The geometric domain is divided into three regions, i.e. region I, II and the grating region, and the coordinate system is defined as shown in Figure 1. Angle of incident \(\theta\) and angle of azimuth \(\varphi\) are defined as shown in Figure 1a as well. If a specific function of refractive index distribution can be predetermined, such as the sinusoidal distribution shown in Figure 1a, the material optical properties can be quantified by some parameters, such as the pitch \(\Lambda\), the height of grating region \(H\), the initial refractive index \(n_{0}\), and the amplitude of the maximum refractive index variation \(n_{A}\). Regardless of magnetic material, the relationship between refractive index \(n\) and relative dielectric constant \(\varepsilon\) (i.e. permittivity) is that \(\varepsilon=n^{2}\). At first, the dielectric material in a period can be uniformly discretized into \(N\) units along \(x\) axis. Without losing generality, supposing the distribution can be described by a sinusoidal function, the refractive indices at an arbitrary position \(x\) can be achieved using interpolation. At the same time, if the optical property variates along \(z\) direction, a commonly used layer-by-layer discretization along \(z\) axis can be carried out as well. As shown in Figure 1b, the grating region is sliced into \(L\) layers in \(z\) direction and \(N\) units in \(x\) direction. Then, the \(x\)-dependent relative dielectric constant at \(l\)-th layer can be expanded as Fourier series  
\[\varepsilon_{l}(x)=\sum\limits_{g}\varepsilon_{l,g}\exp\left(j\frac{2\pi g}{ \Lambda}x\right), \tag{1}\]  
where \(\varepsilon_{l,g}\) is the \(g\)-th component of the Fourier expansion, which can be obtained as  
\[\varepsilon_{l,g}=\left\{ \begin{array}{c}\frac{1}{N}\sum\limits_{n=1}^{N} \varepsilon_{ln}\quad g=0\\ \frac{1}{-j\frac{2\pi g}{N}}\left[\sum\limits_{n=1}^{N-1}\left(\varepsilon_{ ln}-\varepsilon_{l(n+1)}\right)\exp\left(-j\frac{2\pi ng}{N}\right)+\varepsilon_{lN }-\varepsilon_{l1}\right]\quad g\
eq 0\end{array} \right., \tag{2}\]  
Figure 1: Scheme of the (**a**) spatial refractive index distribution in a volume grating; and (**b**) the corresponding distributed dielectric constant model used for rigorous coupled-wave analysis. where \(N\) is the number of units in one period. In the same way, the reciprocal of dielectric constant can be expanded as  
\[\frac{1}{\varepsilon_{l}(x)}=\sum\limits_{g}\left(\frac{1}{\varepsilon}\right)_{ l,g}\exp\left(j\frac{2\pi g}{\Lambda}x\right), \tag{3}\]  
where  
\[\left(\frac{1}{\varepsilon}\right)_{l,g}=\left\{ \begin{array}{c}\frac{1}{N} \sum\limits_{n=1}^{N}\left(\frac{1}{\varepsilon}\right)_{ln}\ \ \ g=0\\ \frac{1}{-j2\pi g}\left[\sum\limits_{n=1}^{N-1}\left(\left(\frac{1}{\varepsilon }\right)_{ln}-\left(\frac{1}{\varepsilon}\right)_{l(n+1)}\right)\exp\left(-j \frac{2\pi ng}{N}\right)+\left(\frac{1}{\varepsilon}\right)_{IN}-\left(\frac{1 }{\varepsilon}\right)_{l1}\right]\ \ \ g\
eq 0\end{array} \right. \tag{4}\]  
The electrical field expressions of region I (\(z<0\)) and II (\(z>H\)) can be achieved as Equations (5) and (6) by Rayleigh expansion  
\[\mathbf{E}_{I}=\mathbf{E}_{inc}+\sum\limits_{i}[R_{i}\exp(-jk_{R,i}\bullet r)], \tag{5}\]  
\[\mathbf{E}_{II}=\sum\limits_{i}\{T_{i}\exp[-jk_{T,i}\bullet(r-Hz)]\}, \tag{6}\]  
where \(\mathbf{E}_{inc}\) is normalized incident electrical field; \(H\) is total height of the volume grating; \(\mathbf{R}_{i}\) is the amplitude of \(i\)-th diffractive wave; \(\mathbf{T}_{i}\) is the amplitude of \(i\)-th transmit wave; \(\mathbf{r}\) is the position vector of an arbitrary point on the plane wave front; and \(\mathbf{k}_{R,i}\) is the wave vector of \(i\)-th transmit wave.