The paper presents derivation of analytical components of S-matrices for an arbitrary planar diffraction grating in the Fourier domain. General formulas for S-matrices attained can be applied within both formulations in the Cartesian and curvilinear metric. A numerical method based on these results can benefit from all previous improvements of the Fourier method, possess the same convergence rates and provides the same result as the FMM. As the new method does not require one to solve the eigenvalue problem it appears to be slightly more efficient than the FMM, though both methods have the same asymptotic numerical complexity proportional to the cubed number of Fourier harmonics.
Analytic S-matrix components for planar grating structures \authorAlexey A. Shcherbakov, Yury V. Stebunov, Denis F. Baidin\\ Thomas Kämpfe, and Yves Jourlin \\\\ Moscow Institute of Physics and Technology, Dolgoprudniy, Russia \\ University of Lyon, Laboratory Hubert Curien, Saint-Etienne, France \\firstname.lastname@example.org
Periodic optical structures ranging from conventional one-dimensional diffraction gratings to metasurfaces attract a lot of attention due to wide optimization possibilities in design of optical response functions. The wavelength scale nature of patterns of the most complex diffractive structures and metasurfaces requires an effort to rigorously solve Maxwell’s equations. Among various methods capable to suit this task the Fourier space methods  are among most popular and important ones. They are widely used owing to their relatively simple formulation, versatility in structures that can be modelled, and an output in form of S-matrices. The latter property features these methods with a yield that can be directly used in simulation of measurable quantitites. Not only are S-matrices physically important entities, but also they bring stability in numerical calculations [2, 3, Li1996-1, 4].
Fourier methods, however, conventionally operate with T-matrices, and to ensure stability of calculations an additional effort is needed, e.g., proposed in [5, 3]. In this paper we propose a way to overcome these complications and demonstrate how S-matrix components of and arbitrary diffractive structure can be derived analytically in the Fourier space. Our derivation is based on the integral equation solution of the Maxwell’s equations, albeit the distribution formalism given in  can be used to get the same result. Current work can be considered as generalization of results obtained in the latter reference.
In addition to derivation of S-matrices of bulk material gratings it is shown that similar expressions can be attained for corrugated layers of 2D materials, which provide possibility to efficiently simulate complex metasurfaces with 2D. This widens applicability of the previous models of [7, 8] where an electrodynamic response of a graphen sheet placed on top of a corrugated substrate was simulated by means of the Rayleigh-type method. Our method does not rely on the Rayleigh hypothesis and hence is free of convergence issues which the Rayleigh-type methods face [9, 10].
2 Volume Integral Equation
To start consider a planar structure, which can be either a metasurface of diffractive optical element, in Cartesian coordinates , with unit vectors , such that axis is orthogonal to structure plane. The structure is supposed to be doubly periodic along two non-collinear directions in plane . Denote unit vectors in the directions of the periods as , and let periods be . Reciprocal lattice vectors then read
The paper refers to linear phenomena only. Thus, we consider Maxwell’s equations for time-harmonic fields and sources with implicit time dependence exponential factor , which will be omitted:
Magnetic source is essential for a curvilinear coordinate formulation as will be clear from further derivations. Eqs. (2) yield Helmholtz equations providing that the dielectric permittivity and the magnetic permeability are constants. We will refer to these quantities as basis ones and denote them as , . Helmholtz equations then read
where is the wavenumber of the homogeneous basis space. A well-known solution of Eqs. (3) in form of the volume integral equation relies on the free-space Green’s functions. As we treat planar structures the Green’s function should be decomposed in the plane wave basis. Given a plane wave wavevector with sign distinguishing waves propagating upwards and downwards relative to axis , unit vectors of the TE and the TM polarized waves can be taken as
When dealing with problems in periodic domain the sources can be decomposed into the Fourier series in the plane of periodicity, so that source Fourier component is
where the integration is performed over one structure period, and plane harmonic wavevector depends on incident plane wave vector and two-dimensional index . Corresponding solution in the periodic domain is
A similar expression holds for the magnetic field. Eq. (9) shows that the electric field is a superposition of plane harmonics and a source ‘’delta‘’-term. Thus, introduction of a modified field as , makes plane wave decomposition of this modified field
These equations can be used to get a formulation of the Generalized Source Method either in Cartesian  or in curvilinear  coordinates by introducing generalized sources related to the fields and performing implicit numerical integration. Instead, in the next section Eqs. (11) are used to obtain analytical S-matrix components of an infinitely thin grating layer. The approach can also be seen as generalization of a theory of .
3 S-matrix of a thin grating slice
Consider sources inside a layer of thickness . If integration in Eqs. (11) reduces to multiplication of integrands by . Denote coordinate of layer center as . Ones incoming wave amplitudes at layer boundaries are known, Eqs. (11) yield diffracted amplitudes in form
These equations are accurate up to terms proportional to .
To get a closed form equations sources and are to be related with the fields. In the simplest case , and , as is done within the Generalized Source Method for index gratings  or in other Volume Integral Equation methods (e.g., ). When a formulation includes a correct treatment of the boundary conditions to a corrugation interface in the Fourier domain  or the generalized metric sources, which appear in a curvilinear formulation of the problem , the dependence of the sources from the field can be more complicated. Generally, such dependence can be taken in form
with matrices . In order to operate with plane harmonic amplitudes only introduce matrices composed of column vectors given by Eq. (4):
Additionally, introduce amplitude vectors of upward and downward propagating harmonics . Then, Eqs. (12) transform to
where are components of the Fourier block-Toeplitz matrices obtained by the Fourier transform of corresponding matrices evaluated at coordinate . Self-consistent field amplitudes in the right-hand part of the equation can be substituted by the incoming filed amplitudes according to the chosen accuracy of the second term, which is proportional to .
Derived Eq. (15) directly provide relations between incoming and outgoing wave amplitudes for diffraction on an arbitrary thin grating layer. Associating incident amplitudes with local fields at boundaries of the layer one can recompile these equations to compose a corresponding -matrix. Given a planar structure bounded by planes and , will refer to the -matrix of this structure as a block matrix that relates plane wave amplitudes at the boundaries in the following way
To apply Eqs. (15), (16) a grating layer should be divided into a number of slices analogous to the FMM and other Fourier methods. However, calculation of eigen modes is no more needed since -matrix components for each slice are readily available. Once -matrix of each slice is attained a corresponding matrix for the whole grating layer can be calculated by the -matrix propagation algorithm known to be numerically stable . The algorithm is based on the following multiplication rule. Given two -matrices of planar structures occupying plane layers and respectively, components of -matrix that relates amplitudes at interfaces and are
Due to presence of matrix inversions in Eq. (17) complexity of the algorithm can be estimated as , where is a number of Fourier harmonics, and is a number of slices.
Matrices generally depend on a particular implementation of the Fourier approach and follow from numerous results attained by different authors on the Fourier Modal Method, Differential Method, C-method, and the Generalized Source Method together with other volume integral implementations (e.g., [16, 18, 19, 20, 21]). For reader to get acquainted with possible output of the previous section we provide here two illustrative examples of widely used gratings followed by a derivation of an S-matrix for corrugated layer of 2D material.
4.1 1D lamellar grating
A 1D lamellar grating is one of the simplest but hereby among the most practically important examples of periodic corrugations. Consider a formulation in the Cartesian coordinates. Then, all components would be zero, and would be diagonal. Correct factorization of the products of discontinuous functions [16, 22] yields
Here the periodicity along is assumed, and denotes inverted truncated Fourier-matrix of periodic permittivity function. This permittivity function explicitly writes , , and , with , , and . This gives elements of the Fourier matrix
In the collinear diffraction case when -matrix splits into two independent parts for the TE and TM polarization separately. We get for the TE polarization
and for the TM polarization
Independence of function of coordinate allows calculating grating -matrix with complexity by sequentially multiplicating thin layer -matrix by itself.
4.2 1D sinusoidal grating in curvilinear coordinates
Necessity in the magnetic sources in the above derivations is due to the previously developed curvilinear coordinate Generalized Source Method with effective generalized metric sources [17, 14]. The core idea is to fit a grating corrugation profile with a suitable cuvilinear coordinate transformation similar to what is done in the C-method , while the chosen curvilinear metric should continuously become Cartesian in a region near the grating. This curvilinear coordinate approach exploits the similarity between Maxwell’s operators in the Cartesian and the curvilinear metric, and a possibility to treat metric contributions as source terms. Under this rationale the generalized sources derived from local metric variations read
where are local metric tensor components, , and summation over the repeated indices is implied. The approach yields matrices that depend on local metric and medium properties :
with , . -matrices for the two polarization then explicitly read:
These equations have a symmetric form relative TE/TM polarizations contrary to the previous example, which is a consequence of the symmetry of the generalized metric sources in Eq. (21). In case of sinusoidal corrugation coordinate transformation is defined as , and , where is the corrugation depth, and is a region with curvilinear metric (see  for details). Components of Fourier matrices for sinusoidally corrugated gratings are found analytically:
with and .
The two considered examples demonstrate the possibility to derive -matrix components explicitly, though resulting expressions can be rather bulky. Therefore here we restrict ourselves by only these two examples for gratings of bulk materials as -matrix components in other various cases can be attained with aids of the vast literature on the Fourier methods and general equations of the previous section. The method has the second order polynomial convergence relative to the slice number, and typical convergence plots are quite similar to those presented for the GSM in [13, 17, 14]. The convergence rate relative the number of Fourier orders for profiled gratings and S-matrices written in Cartesian coordinates is polynomial, and for continuously differentiable profiles with curvilinear coordinate S-matrices is exponential, again, similarly to the GSM , and the GSMCC  respectively.
4.3 Corrugated 2D material
We proceed with modification of the results of the previous subsection and consider a layer of 2D material, e.g., graphene, on top of a corrugated substrate. Here focus on 1D holographic gratings, which profiles can be well approximated by sinusoidal functions introduced above.
Electric current in such materials depends only on tangent electric field components. Therefore, within the rationale of curvilinear coordinate transfomation relation between an effective current and the field via surface conductivity should be written as follows
where upper and lower indices distinguish contravariant and covariant vector components, is surface conductivity, and . It is seen that normalization factor makes effective conductivity in the curvilinear metric be periodic even when the effect of corrugation on conductivity is neglected. Possible impact of the periodicity on conductivity  can be included in our method straightforwardly, though we do not account for it in the following examples. Analogously to the previous section substitution of the sources into Eq. (11) yields
Note that after we take limit the modified field no more differs from the real one due to the absence of the third field component in Eq. (25). Electric current at layer location depends on continuous tangential electric field. In case of 1D corrugation along direction and collinear diffraction with Eq. (26) reduces to two systems for the TE and TM polarization, which read in the S-matrix form as follows
where , . Inversion implies that matrices should be composed of corresponding truncated Fourier matrices and then numerically inverted. In case of 1D grating and Eq. (27) simplifies to known reflection and transmission coefficients for the TE polarization . Similar does Eq. (28) for the TM polarization in the absence of the corrugation. Also it is important to note that Eqs. (27,28) should be applied together with the results of the previous subsection.
Summarizing the results we derived a general form of S-matrix in the Fourier basis Eq. (15) from volume integral solutions of Maxwell’s equations. Provided examples include derivation of explicit analytical S-matrix components for 1D lamellar and sinusoidally corrugated gratings of bulk materials. S-matrices for other types of gratings including complex 2D periodic diffractive optical elements and metasurfaces can be attained in a similar manner on the basis of all previous results for the Fourier Modal Method (e.g., [18, 19, 20, 21]), and there is no more need in solution of the eigenvalue problem in each slice . Moreover the derived analytical S-matrix components can be directly used for resonant analysis of grating structures proposed in . The obtained S-matrix components are accurate up to terms in spatial discretization of the grating layer, which is similar to the other Fourier methods . In addition we expressed the S-matrix of a corrugated layer of 2D material with a given conductivity, which can be useful for further research of optical response functions of metasurfaces covered with such materials.
The work was supported in by the Russian Scientific Foundation, Grant No. 17-79-20345.
-  E. Popov, ed., Gratings: Theory and Numeric Applications. Institut Fresnel, AMU, 2012.
-  D. Y. K. Ko and J. C. Inkson, “Matrix method for tunneling in heterostructures: Resonant tunneling in multilayer systems,” Phys. Rev. B, vol. 38, pp. 9945–9951, 1988.
-  N. P. K. Cotter, T. W. Preist, and J. R. Sambles, “Scattering-matrix approach to multilayer diffraction,” J. Opt. Soc. Am. A, vol. 12, pp. 1097–1103, 1995.
-  N. A. Gippius and S. G. Tikhodeev, “The scattering matrix and optical properties of metamaterials,” Phys. Usp., vol. 52, p. 967â971, 2009.
-  M. G. Moharam, D. A. Pommet, and E. B. Grann, “Stable implementation of the rigorous coupled-wave analysis for surface-relief gratings: enhanced transmittance matrix approach,” J. Opt. Soc. Am. A, vol. 12, pp. 1077–1086, 1995.
-  J. E. Sipe, “New Green-function formalism for surface optics,” J. Opt. Soc. Am. B, vol. 4, pp. 481–489, 1987.
-  Y. V. Bludov, A. Ferreira, N. M. R. Peres, and M. I. Vasilevskiy, “A primer of surface plasmon-polaritons in graphene,” Int. J. Mod. Phys. B, vol. 27, pp. 1341001–74, 2013.
-  T. M. Slipchenko, M. L. Nesterov, L. Martin-Moreno, and A. Y. Nikitin, “Analytical solution for the diffraction of an electromagnetic wave by a graphene grating,” J. Opt., vol. 15, pp. 114008–11, 2013.
-  J. Wauer and T. Rother, “Considerations to Rayleigh’s hypothesis,” Opt. Comm., vol. 282, pp. 339–350, 2009.
-  A. V. Tishchenko, “Rayleigh was right: Electromagnetic fields and corrugated interfaces,” Opt. Express, vol. 21, pp. 50–54, 2010.
-  L. B. Popov and N. Marcuvitz, Radiation and Scattering of Waves. Prentice-Hall, 1972.
-  B. Amorim, P. A. D. Goncalves, M. I. Vasilevskiy, and N. M. R. Peres, “Impact of graphene on the polarizability of a neighbour nanoparticle: A dyadic Greenâs function study,” Appl. Sci., vol. 7, p. 1158, 2017.
-  A. A. Shcherbakov and A. V. Tishchenko, “New fast and memory-sparing method for rigorous electromagnetic analysis of 2d periodic dielectric structures,” J. Quantitative Spectrosc. Radiat. Transfer, vol. 113, pp. 158–171, 2012.
-  A. A. Shcherbakov and A. V. Tishchenko, “Generalized source method in curvilinear coordinates for 2d grating diffraction simulation,” J. Quantitative Spectrosc. Radiat. Transfer, vol. 187, pp. 76–96, 2017.
-  O. J. F. Martin and N. B. Piller, “Electromagnetic scattering in polarizable backgrounds,” Phys. Rev. E, vol. 58, pp. 3909–3915, 1998.
-  L. Li, “Use of Fourier series in the analysis of discontinuous periodic structures,” J. Opt. Soc. Am. A, vol. 13, pp. 1870–1876, 1996.
-  A. A. Shcherbakov and A. V. Tishchenko, “Efficient curvilinear coordinate method for grating diffraction simulation,” Opt. Express, vol. 21, pp. 25236–25247, 2013.
-  E. Popov and M. Nevi’ere, “Maxwell equations in Fourier space: fast-converging formulation for diffraction by arbitrary shaped, periodic, anisotropic media,” J. Opt. Soc. Am. A, vol. 18, pp. 2886–2894, 2001.
-  G. Granet and J.-P. Plumey, “Parametric formulation of the fourier modal method for crossed surface-relief gratings,” J. Opt. A: Pure Appl. Opt., vol. 4, pp. S145–S149, 2002.
-  T. Schuster, J. Ruoff, N. Kerwien, S. Rafler, and W. Osten, “Normal vector method for convergence improvement using the rcwa for crossed gratings,” J. Opt. Soc. Am. A, vol. 24, pp. 2880–2890, 2007.
-  M. C. van Beurden and I. D. Setija, “Local normal vector field formulation for periodic scattering problems formulated in the spectral domain,” J. Opt. Soc. Am. A, vol. 34, pp. 224–233, 2017.
-  P. Lalanne and G. M. Morris, “Highly improved convergence of the coupled-wave method for TM polarization,” J. Opt. Soc. Am. A, vol. 13, pp. 779–784, 1996.
-  J. Chandezon, D. Maystre, and G. Raoult, “A new theoretical method for diffraction gratings and its numerical application,” J. Opt. Paris, vol. 11, p. 235, 1980.
-  A. Isacsson, L. M. Jonsson, J. M. Kinaret, and M. Jonson, “Electronic superlattices in corrugated graphene,” Phys. Rev. B, vol. 77, pp. 035423–6, 2008.
-  T. Stauber, N. M. R. Peres, and A. K. Geim, “Optical conductivity of graphene in the visible region of the spectrum,” Phys. Rev. B, vol. 78, pp. 085432–8, 2008.
-  M. Liscidini, D. Gerace, L. C. Andreani, and J. E. Sipe, “Scattering-matrix analysis of periodically patterned multilayers with asymmetric unit cells and birefringent media,” Phys. Rev. B, vol. 77, pp. 035324–11, 2008.
-  D. A. Bykov and L. L. Doskolovich, “Numerical methods for calculating poles of the scattering matrix with applications in grating theory,” J. Lightwave Technol., vol. 31, pp. 793–801, 2013.
-  E. Popov, M. Nevi‘ere, B. Gralak, and G. Tayeb, “Staircase approximation validity for arbitrary-shaped gratings,” J. Opt. Soc. Am. A, vol. 19, p. 33â42, 2002.