Magnetic‐Electric Metamirror and Polarizing Beam Splitter Composed of Anisotropic Nanoparticles

The emergence of new materials and fabrication techniques provides progress in the development of advanced photonic and communication devices. Transition metal dichalcogenides (e.g., molybdenum disulfide, MoS2) are novel materials possessing unique physical and chemical properties promising for optical applications. In this paper, a metasurface composed of particles made of bulk MoS2 is proposed and numerically studied considering its operation in the near‐infrared range. In the bulk configuration, MoS2 has a layered structure being a uniaxial anisotropic crystal demonstrating an optical birefringence property. It is supposed that the large‐scale and uniform MoS2 layers are synthesized in a vertical‐standing morphology, and then they are patterned into a regular 2D array of disks to form a metasurface. The natural anisotropy of MoS2 is utilized to realize the splitting of electric and magnetic dipole modes of the disks while optimizing their geometric parameters to bring the desired modes into overlap. At the corresponding resonant frequencies, the metasurface behaves as either an electric or a magnetic mirror, depending on the polarization of incident light. Based on the extraordinary reflection characteristics of the proposed metasurface, it can be considered an alternative to traditional mirrors and optical splitters when designing compact and highly efficient metadevices, which provide polarization and phase manipulation of electromagnetic waves on a subwavelength scale.


Introduction
Polarization-selective mirrors, optical switches, and beam splitters are the key components in photonic integrated circuits used for optical signal processing, communication, and sensing. Conventional polarization/selective components are usually made DOI: 10.1002/andp.202300111 of prisms or birefringent substrates and suffer from inherently bulky and voluminous configurations, which do not satisfy the current trend of integration and miniaturization in photonics. A modern approach is to construct these components on a metasurfacebased platform. [1][2][3] They are low-loss planar nanostructures composed of arrays of dielectric resonators with a subwavelength spatial resolution whose production is compatible with available on-chip fabrication and integration technologies.
The optical features of dielectric metasurfaces depend on the excitation conditions of corresponding Mie-type modes [4][5][6] arising in resonators. In particular, in the spectra of dielectric metasurfaces, there are two basic resonant states that correspond to the excitation of the magnetic dipole (MD) and electric dipole (ED) modes, whose position on the frequency scale depends on the aspect ratio of the particles. Tuning mode conditions by changing the aspect ratio of dielectric particles [7][8][9] or their symmetry properties [10,11] makes it possible to flexibly manipulate the resonant states of nanostructures.
Exploiting Mie-type modes in dielectric nanoparticles also allows one to create reflective metasurfaces (metamirrors). It is reported [12][13][14] that the dielectric metamirrors can reach near-unity reflection performance surpassing the reflectance of conventional metallic mirrors. [15] Remarkably, since nanoparticles in dielectric metasurfaces simultaneously support modes of the electric and magnetic types, both electric and magnetic mirrors can be implemented. [16][17][18] The difference between these two mirrors is related to their reflection phase: [19] the tangential component of the electric field of the wave reflected from the array of particles bearing the ED mode experiences a phase shift at the resonance, whereas no phase shift occurs after reflection from the array of particles bearing the MD mode. Moreover, by providing fine-tuning of the mode interaction in dielectric nanoparticles, full control over the reflected wavefront can be achieved. [20,21] This ability to flexibly tune both reflection and transmission characteristics in dielectric metasurfaces opens the prospect of realizing polarization beam splitters and polarization converters with improved characteristics. [22][23][24][25][26] Optical anisotropy plays a crucial role in light manipulation in many optical systems owing to the birefringence phenomena. [27,28] It was previously reported [29,30] that in diskshaped resonators, the birefringence property removes the degeneracy for modes, resulting in their splitting. Moreover, the anisotropy breaks the symmetry of the resonators, and additional leakage may appear for the axially symmetric modes, whereas the quality factor of the other modes may increase. Utilizing the birefringence property can introduce an additional degree of freedom and pave the way for a new generation of optical compact planar devices and systems based on anisotropic particles. In particular, such novel materials as transition metal dichalcogenides (TMDs) [31] have the required anisotropy. Moreover, with the use of TMDs, one can expect to make the metasurface-based platforms manageable.
Atomically thin TMDs are formed by a layer of transitionmetal atoms sandwiched between two layers of chalcogenide atoms. In the bulk form, TMDs appear as a layered structure composed of monolayer nanosheets. Therefore, bulk TMD is an anisotropic crystal that can be macroscopically characterized with a tensor-valued permittivity possessing in-plane ∥ and out-ofplane ⟂ terms. [32,33] With the use of the chemical vapor deposition method, the synthesis of large-scale and uniform monolayers with either horizontally [34] or vertically [35,36] standing morphology is possible. These extraordinary properties of TMDs have already been utilized in the development of anisotropic resonators for lasing systems [37] and dielectric metasurfaces. [38][39][40] The anisotropy property expands the number of possible effects that can be observed in metasurfaces composed of TMD anisotropic particles and options for their practical use. In particular, we have already shown earlier that anisotropy changes the properties of modes of TMD (MoS 2 ) resonators and can turn their operation from dark to bright states and vice versa. [41,42] To develop this concept further, in this paper, we combine the anisotropic properties of bulk MoS 2 that is considered for particle fabrication, and geometric parameters of these particles to construct a reflective metasurface that simultaneously exhibits the properties of either an electric or a magnetic mirror and beam splitter, depending on the polarization of incident light. In particular, in our metasurface, we utilize disk-shaped resonators due to their amenability to modern lithography techniques [43] used for producing metasurfaces operating in the near-infrared range (within the third telecommunication band). Our approach is based on the analysis of the eigenwaves of such anisotropic res-onators and the search for optimal parameters for the period of the metasurface and aspect ratio for its constitutive nanoparticles to obtain the specific crossing conditions for the corresponding MD and ED modes. To describe the effect of birefringence and parameters of the polarization transformations in the metasurfaces, we apply the multipole decomposition method developed for the transmission and reflection coefficients. [44] 2. Results and Discussion

Eigenwaves Analysis
In what follows, we consider that a metasurface is composed of MoS 2 anisotropic disk-shaped nanoparticles arranged in the array with a square unit cell of size P. The diameter and height of the disks are denoted as D and H, respectively. The disks are made from bulk MoS 2 whose nanosheets are deposited in a vertically-standing morphology. The use of multilayer MoS 2 for the fabrication of the nanoparticles has advantages over conventional semiconductors such as due to its high refractive index and low level of infrared losses. [33] Bulk MoS 2 is a nonmagnetic uniaxial anisotropic crystal whose permittivity is described by a tensor quantity . For the sake of definiteness, we assume that the anisotropic resonators are manufactured and arranged in the metasurface in such a way that their anisotropy axis coincides with the direction of the x-axis of the chosen Cartesian coordinate frame, as shown in Figure 1a (in this coordinate frame, the horizontal plane at z = 0 is related to the mid-height of the resonators). For such a geometry, the particles permittivity is defined in a diagonal form = [ ⟂ , ∥ , ∥ ], where ⟂ − ∥ < 0 (negative birefringence). Since in the chosen operating frequency range, material losses in bulk MoS 2 are negligibly small (see Appendix A.1), we exclude them from our consideration. Without loss of generality, in this study, we assume that the particles are located in a homogeneous surrounding space (air). In order not to overload our research, we ignore the presence of the substrate since it is insignificant to the effects discussed. Accounting for the substrate can be performed later in the engineering implementation of the metasurface by optimizing the geometric parameters of the resonators used.
At the initial stage, we perform an analysis of the eigenwaves of the given metasurface. To find complex eigenfrequencies of our metasurface, we use the commercial COMSOL Multiphysics finite-element electromagnetic solver. In this solver, we perform full-wave eigenfrequency analysis in RF Module, frequency domain. In our model, a single unit cell of the metasurface with the corresponding constitutive and geometrical parameters is implemented, whereas the 2D array of resonators is imitated by applying the periodic (Floquet) boundary conditions and repeating the unit cell infinitely in both transverse directions. The thick perfectly matched layers are applied in the vertical direction at some distance above and below the metasurface to circumvent the boundary reflections. In this study, we use the default physicscontrolled COMSOL settings for fine mesh generation.
From the solutions found, we select those that correspond to the excitation of the lowest-order (dipole) modes only. They are the hybrid EH 11 and HE 11 eigenwaves [45] of a cylindrical resonator (the scattering cross-section of a single anisotropic MoS 2 nanoparticle with the far-field radiation patterns of correspond- ing modes can be found in Appendix A.2). These eigenwaves manifest themselves like an electric dipole p and magnetic dipole m, respectively, lying in the metasurface plane. Since the constitutive nanoparticles are made of anisotropic material, the hybrid EH 11 and HE 11 eigenwaves lose their degeneracy and split apart. [42] They can be distinguished in terms of ordinary (o) and extraordinary (e) modes with the corresponding dipoles oriented perpendicular and parallel to the anisotropy axis, as accepted in the theory of anisotropic media [27] (i.e., the ordinary mode is related to the case when the electric (magnetic) dipole moment is oriented perpendicular (parallel) to the anisotropy axis and vice versa for the extraordinary one). In what follows, for brevity, we denote them as EH o(e) 11 ⇔ ED o(e) and HE o(e) 11 ⇔ MD o(e) modes. In particular, a pair of low-frequency ordinary modes (ED o and MD o ) and a pair of high-frequency extraordinary modes (ED e and MD e ) arise whose frequency separation depends on the difference in the values of ∥ and ⟂ of the MoS 2 crystal. The corresponding eigenfrequencies and eigenfields for these four modes of interest are presented in Figures 1b,c, respectively. One can conclude that for the given metasurface, the eigenfrequencies of ED o and MD e modes appear to be closely spaced, and further, we utilize this feature in our designs.
Thus, to construct the reflective metasurface, which simultaneously behaves as an electric and magnetic metamirror, these two closely spaced ED o and MD e modes must be tuned to exhibit the crossing conditions at the same frequency. [17] To find the corresponding parameters for crossing, we fix in our simulations the diameter D of the disks while varying their height H. The simulations are made for two particular values of the period P of the array providing the period is less than the wavelength of the incident radiation (subwavelength conditions). The results obtained indicate that changing the period affects all modes equally, and as the period decreases, the eigenfrequencies of all modes increase. Nevertheless, changing the height of the disks has a different effect on the ED and MD modes. [7,46,47] As the height of the disks increases, the eigenfrequencies of all four modes acquire downshifting, where the decrease in frequency is more significant for the MD modes. This can bring the ED o and MD e modes into overlap.
Therefore, the H variation at fixed D and P allows us to determine the specific height of the disks where the crossing of dispersion curves of the ED o and MD e modes occurs. With such a parameter tuning, an optimal geometrical configuration of the metasurface can be found, and both the crossing frequency (f cr ) www.advancedsciencenews.com www.ann-phys.org of the ED o and MD e modes and the corresponding disks' height (H cr ) can be determined.

Birefringence Manifestation
Since the metasurface is composed of anisotropic nanoparticles, all peculiarities of the birefringence manifestation must be elucidated. For this study, we consider that the linear polarization of the incident field is oriented at a certain angle with respect to the x-axis. To reveal how the response of the metasurface composed of anisotropic nanoparticles is formed, we use a strategy based on the component multipole analysis in the coordinate system associated with the incident field polarization. Thus, we introduce an auxiliary coordinate system x ′ -y ′ rotated at an angle relative to the primary coordinate system x-y in such a way that the polarization of the incident field E 0 is oriented along the x ′axis (see Figure A3 in Appendix A.3). We analyze the polarization characteristics of the metasurface using a set of the perpendicular (⟂) and parallel (||) coefficients R ⟂(||) and T ⟂(||) for the reflected and transmitted waves, respectively, with the electric field perpendicular (parallel) to the x ′ -axis. When two coordinate sys- and T o(e) are the corresponding coefficients related to the reflected and transmitted ordinary (extraordinary) waves.
To explicitly show that the optical properties of the given metasurface are determined by the MD and ED moments of the disks, we apply the multipole decomposition of the metasurface's reflection and transmission coefficient. [44] We denote these analytical coefficients with the superscript 'an' (R an ⟂(∥) and T an ⟂(∥) , see Appendix A.3). Then we compare the results of our analytical calculations with full-wave numerical simulations made in the COM-SOL solver. To calculate these coefficients in the solver, we designate the radiating and receiving ports above and below the metasurface, respectively. The port boundary conditions are placed on the interior boundaries of the perfectly matched layers. This allows us to determine the reflection and transmission characteristics of the metasurface in terms of S-parameters related to the ports (R ⟂(∥) = |S11| 2 and T ⟂(∥) = |S21| 2 ). The results of our calculations are collected in Figure 2.
One can readily see that the curves for R an ⟂(∥) (T an ⟂(∥) ) approximate well R ⟂(∥) (T ⟂(∥) ), which indicates the negligible contribution of high-order multipole terms. At the same time, one can see that for the given metasurface, the resonant features of the reflection and transmission spectra around f ≈ 200 THz arise due to the MD and ED modes manifestation. The direct transmission coefficient T ∥ associated with the incident polarization depends on both R ⟂ and R ∥ reflection coefficients [see Equation ( A2) in Appendix A.3], as a result, its suppression shown in Figure 2a is a consequence of the resonant reflection R ⟂ with the orthogonal polarization (see Figure 2b). This resonant reflection, in turn, is a result of the "anti-Kerker" effect, [48] when both the MD ∥ and ED ⟂ moments, induced in each disk, emit electromagnetic waves in phase (out of phase) in the backward (forward) direction. The destructive interference between waves scattered in the forward direction (transmitted wave) by the MD ∥ and ED ⟂ moments also leads to the suppression of the transmission T ⟂ at f ≈ 200 THz (see Figure 2b). Note that T ⟂ is determined only by the radiation emitted by the disk's dipole (multipole) moments [see Equations ( A1) and ( A5) in Appendix A.3]. In contrast to R ⟂ , the reflection coefficient R ∥ shown in Figure 2a is significantly reduced due to the destructive interference between the contributions of MD ⟂ and ED ∥ in the reflection (the Kerker effect [48] ).
The most pronounced polarization transformation in the orthogonal component happens for H ≈ 340 nm and f ≈ 200 THz. At these parameters, the reflected field appears to be elliptically polarized (see Figure 2c), and the polarization ellipse is rotated anticlockwise at the angle = ∕4 = 45 • . Thus, it confirms that the maximal orthogonal polarization component is obtained for the incident polarization with = ± ∕4 (see Figure 2d).
The birefringence characteristics of the MoS 2 nanoparticles give us the ability to manipulate the polarization states of the reflected waves. Moreover, in the next section, we show how the resonant anisotropic properties of the metasurface's building blocks can be adopted for the realization of a magnetic-electric metamirror.

Magnetic-Electric Metamirror
In this section, we consider the reflection properties of the metasurface at the ED o and MD e mode crossing conditions, when = 0 the metasurface is irradiated by either x-polarized or ypolarized waves. The result of our calculations of the R an o(e) coefficients based on the multipole decomposition method is shown in Figure 3a where the crossing point of the resonant reflection related to the excitation of the ED o and MD e modes in the disks can be revealed. Since these modes can be excited only by a wave with a proper orientation of its electric field vector, the resonances are polarization-dependent. For the chosen geometry of the problem, the excitation of the MD e and ED o modes is carried out by the E x and E y components of the incident wave, respectively.
In Figure 3b, the reflected spectra of the metasurface irradiated by the two orthogonal linearly polarized waves are presented considering three different values of the disks' height. The crossing point of the resonant reflection is highlighted in these figures with gray arrows. One can conclude that close to the eigenfrequencies of either MD or ED modes, the metasurface re-flects the incident wave with high efficiency. Since the obtained resonances have a Fano form rather than the Lorentzian, the frequencies of the resonant reflection are slightly different from the derived eigenfrequencies of modes (each eigenfrequency is shifted to a middle point of the slope of the corresponding Fano resonance). Nevertheless, as expected from the eigenwaves analysis, the corresponding MD and ED resonances appear in the higher frequency range for the reflected x-polarized (extraordinary) wave compared to those for the reflected y-polarized (ordinary) wave. For the particular height of disks (H cr = 340 nm), the frequencies of MD e and ED o resonances are coincident, and our metasurface resonantly reflects the waves of both orthogonal linear polarizations.
The appearance of the resonant reflected field is schematically shown in Figure 3c. It is plotted for both orthogonal linear polarizations of the incident wave for the parameters where the crossing of the MD e and ED o modes occurs. The arrangement of red arrows indicates the appearance of the electric field. The boundary of the mirror lies on the upper surface of the resonators, forming the metasurface (z = H∕2). One can conclude that at the resonant frequency, a standing wave appears from this boundary for waves of both orthogonal polarizations. For the incident y-polarized wave, the reflected ordinary wave formed by the ED o mode has a node at the mirror surface where there is a minimum of the electric field amplitude. This is similar to the wave reflection conditions for the perfect electric conductor, and, thus, our metasurface behaves as an electric mirror. Contrariwise, for the incident x-polarized wave, the reflected extraordinary wave produced by the MD e mode has an anti-node at the surface, a point where the amplitude of the standing wave is at maximum, and the behaviors of a magnetic mirror arise. The amplitude patterns show that the electric field is significantly perturbed near the boundary of the metasurface, while at some distance from this boundary, a plane front of the reflected wave is formed, where the minimum amplitude of the reflected ordinary wave corresponds to the maximum amplitude of the reflected extraordinary wave.
Since the metasurface is composed of anisotropic particles, its eigen-polarization state is elliptical. Therefore, we perform a study of the reflection feature of our metasurface with respect to all polarization states of the incident wave. For this study, the reflection coefficient magnitude of the metasurface for different polarization states is mapped onto the surface of the full Poincaré sphere presented in Figure 4a. In our calculations, the electric field vector E 0 of the incident wave is defined by components E x = [0.5(Ψ + 1)] 1∕2 and E y = [0.5(Ψ − 1)] 1∕2 exp(i ), where Ψ ∈ [−1, 1] and ∈ [− , ]. [49][50][51] The sphere is related to the reflection characteristics of the metasurface at the ED o and MD e mode crossing conditions. The axes of the Poincaré sphere are labeled in terms of the Stokes parameters. Linearly polarized states are arranged along the equator of the sphere, while right-handed (RCP) and left-handed (LCP) circularly polarized waves are located at their north and south poles, respectively. All other points of the sphere depict elliptically polarized states of the incident wave. The uniformly colored appearance of the sphere indicates that the metamirror exhibits total reflection for all possible polarization states of the incident wave.
In particular, the metasurface reflects an incident linearly polarized wave without changing its polarization, when the direction of the vector E 0 of the linearly polarized incident wave is orthogonal or collinear to the direction of the MoS 2 crystal axis; in other cases, the reflected wave is elliptically polarized. Nevertheless, there is a certain peculiarity in the reflection of waves with circular polarization. As Figures 4b,c suggest, a phase dif-ference of appears between the RCP and LCP waves at the crossing resonant frequency f cr . Therefore, the polarization helicity of the waves incident on the metasurface is preserved, and the metamirror reflects circularly polarized waves keeping their original handedness. This is a distinctive feature of the given metasurface, which cannot be achieved in either purely electric or purely magnetic mirrors (we should note here that previously this feature was obtained for a dielectric metamirror composed of silicon cuboid resonators [17] ).
One can also notice in Figure 3b that near the frequency of 240 THz, the interference of the ED e mode with either an electric quadrupole (EQ) mode or magnetic quadrupole (MQ) mode may occur. Figure 4b indicates that due to such an interference, the maximal polarization transformation between orthogonally polarized waves arises. This effect deserves separate consideration and may be the subject of future research.

Polarizing Beam Splitter
The extraordinary features revealed by our anisotropic metasurface can be utilized to construct a polarizing beam splitter. For the given metasurface, all resonances of interest have an asymmetric (peak-and-trough) Fano profile with alternating total reflection and total transmission, as is typical for low-loss metasurfaces. [52] Due to the anisotropy of MoS 2 nanoparticles, the position of the Fano resonances associated with the excitation of the ED and MD modes is different on the frequency scale for waves of orthogonal polarizations. Therefore, by choosing an appropriate aspect ratio of the disks, particular frequencies can be defined where the coincidence of the peak and trough appears for the Fano resonances associated with the excitation of the closely spaced ED o and MD e modes (we should note here that previously this specific property of the total transmission and total reflection in metasurfaces bearing asymmetric Fano resonances has been utilized for performing both linear [53,54] and nonlinear bistable [55][56][57] optical switching).
To estimate the crossing parameters for the splitting, we again employ the multipole decomposition method. The corresponding parametric dependences of the reflection and transmission spectra on the disks' height are presented in Figures 5a,b in the form of surface plots. From these plots, one can conclude that polarization splitting is possible for the metasurface when the resonant frequency of the total reflection of either ED o mode (electric metamirror splitter) or MD e mode (magnetic metamirror splitter) coincides with the total transmission frequency of its mode counterpart. In this way, the electric metamirror splitter transmits and reflects the x-polarized and y-polarized waves, respectively, at the frequency f 1 sp that is below the crossing frequency f cr of the ED o and MD e modes. This is possible for the disk-shaped resonators with the higher height (H 1 sp > H cr ). Contrariwise, the magnetic metamirror splitter transmits and reflects the y-polarized and x-polarized waves, respectively, at the higher frequency f

Conclusion
In conclusion, we have demonstrated that the metasurface made of an array of MoS 2 anisotropic nanoparticles behaves as either an electric or a magnetic mirror for the waves of orthogonal linear polarizations. These dual properties of the metasurface are associated with the characteristics of the electric and magnetic dipole modes, which are spectrally split in anisotropic nanoparticles. It is also revealed that such a magnetic-electric metamirror reflects the circularly polarized waves keeping their original handedness.
We have used the extraordinary reflection properties of the metasurface to implement a polarizing beam splitter. To do this, we match the corresponding frequencies of total reflection and total transmission of the Fano resonances for the orthogonal linearly polarized waves.
Although the theory presented here was obtained for resonators made from MoS 2 , it remains applicable for description metasurfaces composed of resonators made of any other uniaxial anisotropic material including monocrystalline ones. Never- Figure A1. Relative permittivity of MoS 2 extracted from experimental data of Ref. [33].
theless, due to such features of MoS 2 as the presence of direct and indirect bandgaps on the micrometer and nanometer scale, and its tunable dielectric constant, which can be controlled using the intrinsic carrier density and voltage bias, there are prospects to realize tunable highly efficient splitters based on our design.

A.1. MoS 2 Relative Permittivity
Dispersion characteristics of components of the tensor of the relative permittivity ( = ′ + i ′′ ) of bulk MoS 2 are shown in Figure A1. Since in our work, we focus on the spectral range near the telecom wavelength of 1500 nm, where negligibly weak dispersion exists (see the spectral region marked gray), we adopt ′ ∥ = 16.6, ′ ⟂ = 6.3, ′′ ∥ = ′′ ⟂ = 0.0 in all considered region.

A.2. Single Disk Response
We performed additional simulations on a single anisotropic MoS 2 nanoparticle scattering in order to prove that the discussed resonant states are influenced by the modes of an individual resonator rather than by the electromagnetic coupling (periodic effect) between them. It is supposed that the disk is illuminated by a linearly polarized plane wave propagating along the disk's axis (z-axis, frontal excitation). The normalized scattering cross-section (SCS) of a disk-shaped resonator is presented in Figure A2 for two orthogonal polarizations of the incident wave in the frequency band of interest. For this simulation, we implemented the multipole decomposition method [58] in the COMSOL Multiphysics electromagnetic solver. One can find a correspondence between the positions of the SCS resonances on the frequency scale, eigenfrequencies of the corresponding ED and MD modes of the single nanoparticle, and those existing in the unit cell of the metasurface under study (as noted above in the main body of the manuscript, all the discussed resonances are Fano-shaped and each eigenfrequency is shifted to the midpoint of the slope of the corresponding resonance on the SCS curve; see also Figures 1b and 3b).

A.3. Reflection and Transmission Coefficients
To analyze the polarization transformations for the incident wave and the wave transmitted through the metasurface, we proceed to the new coordinate system x ′ -y ′ rotated relative to the x-y coordinate system at the same  angle as the polarization of the incident wave (see Figure A3). In this case, all multipole components, as well as the polarization components of the coefficients R and T are denoted by subscript ∥ corresponding to the waves polarized along the x ′ -axis (parallel to the incident wave polarization E 0 ) and subscript ⟂ corresponding to the waves polarized along the y ′ -axis (perpendicular to E 0 ). Using the expansion up to the quadrupole order terms, the corresponding coefficients take the form: where the amplitudes of the corresponding coefficients can be calculated through the combination of the corresponding multipoles as follows: where k d = k 0 √ d is the wave number in a surrounding medium with permittivity d = 1, v d = 1∕ √ 0 0 d is the speed of light in the surrounding nonmagnetic medium, E 0 is the field calculated by direct numerical simulation in COMSOL Multiphysics for an empty box, S L = P 2 is the area of a lattice unit cell, P is the lattice constant, k 0 is the wave number in a vacuum, 0 and 0 are the vacuum permittivity and permeability, respectively. We take the corresponding formulas from Ref. [44] since their derivation is invariant to the rotation of the coordinate system. However, we calculate the multipole components p x , p y , m x , m y , and others through the components of the scattered field, [44] which we obtain from COMSOL Multiphysics in the usual coordinate system x-y. Therefore, it is necessary to recalculate these components. Taking into account that the z component does not change, the electric and magnetic dipoles d ′ = is the rotation matrix.