Scattering matrix for chiral harmonic generation and frequency mixing in nonlinear metasurfaces

We generalize the concept of optical scattering matrix (S-matrix) to characterize harmonic generation and frequency mixing in planar metasurfaces in the limit of undepleted pump approximation. We show that the symmetry properties of such nonlinear S-matrix are determined by the metasurface symmetries at the macroscopic and microscopic scale. We demonstrate that for description of degenerate frequency mixing processes such as optical harmonic generation, the multidimensional S-matrix can be replaced with a reduced two-dimensional S-matrix. We show that for metasurfaces possessing specific point group symmetries, the selection rules determining the transformation of the reduced nonlinear S-matrix are simplified substantially and can be expressed in a compact form. We apply the developed approach to analyze chiral harmonic generation in nonlinear metasurfaces with various symmetries including rotational, inversion, in-plane mirror, and out-of-plane mirror symmetries. For each of those symmetries, we confirm the results of the developed analysis by full-wave numerical calculations. We believe our results provide a new paradigm for engineering nonlinear optical properties of metasurfaces which may find applications in active and nonlinear optics, biosensing, and quantum information processing.


Introduction
Optical properties of metasurfaces composed of subwavelength planar metallic or dielectric elements have been attracting increased attention in the recent years [1,2,3,4,5].
Scattering properties of metasurfaces are substantially determined by the shape and arrangement of their meta-atoms.The description of the scattering processes can be conducted with the scattering matrix (S-matrix) formalism, which provides a relation between waves incoming at and outgoing from the metasurface.Breaking the symmetries of meta-atoms can be used for tailoring sharp resonant features in the scattering spectrum [6], polarization control [7], and other modifications of the transmitted and reflected wave amplitude and phase.The magnitude of these effects can be predicted by analysing the relative values of the elements of Smatrix.Numerically, S-matrix can be computed via different methods, including Fourier modal method [8], decomposition into quasi-normal modes [9,10,11], and temporal coupled-mode theory [12,13].
Enhancement of optical chirality and chiroptical effects in metasurfaces has been actively discussed in the past decade [14].In contrast to natural chirality associated with a microscopic material structure of chiral ions and molecules, induced chiroptical response in metasurfaces arises because of broken geometrical symmetries of meta-atoms and their relative position.Chiroptical phenomena, indicative of optical activity, can be analysed through two distinct optical properties, optical rotation and circular dichroism.
The latter describes of the unequal transmission or absorption of circularly polarized light within a chiral medium.Recently, the concept of optical S-matrix was utilized to derive the conditions of strong optical chirality and high circular dichroism depending on the symmetries of meta-atoms and their resonant properties [15,16].Specific properties of transmitted and reflected chiral light were linked to the metasurface symmetry, such as absence of polarization conversion in lattices with rotational symmetry [17].Very recently, S-matrix was used to design single handedness cavities [18], and to achieve maximal optical chirality in metasurfaces via the engineering of coupling to sharp optical resonances [17,19].
Lately, the primary emphasis in optical metasurface research has transitioned towards nonlinear metasurfaces [20].In the last decade, it was shown that the nonlinear signal generated due to microscopic nonlinearities of metals or dielectrics can be enhanced strongly by engineering optical resonances [21,22,23,24,25].In particular, special interest was attracted to sum-frequency and optical harmonic generation in the undepleted pump regime [26,27], which goes beyond the phase-matching approach typical for nonlinear optics of macroscopic structures [28].Very recently, it was shown that generation of higher-order harmonics can be achieved for low input power by virtue of sharp optical modes [29,30,31].
Chiroptical activity in nonlinear metastructures and metasurfaces has been a topic of very intense studies [32,33,34,35,36].Analogous to linear phenomena, the enhancement of nonlinear chiroptical effects beyond the constraints of inherent material responses can be realized through breaking macroscopic symmetries of metasurface lattice and meta-atoms.It was shown that the efficiency of generation of chiral harmonics can be increased by employing resonances in metasurfaces [37,38,39], and maximal nonlinear chirality can be achieved [40,41].
Pioneering studies of nonlinear chirality in optically active natural crystals and bulk media revealed that the nonlinear optical activity depends on the microscopic symmetry of susceptibility and permittivity tensors, and obeys a set of selection rules [42,43].Several later experimental studies showed that similar selection rules are also valid for nonlinear metasurfaces, meaning that only specific chiral harmonic generation processes are allowed with the restrictions imposed by the structure lattice symmetry [44,45,39].The question arises if a universal and simple approach to analyse nonlinear chiroptical response of metasurfaces can be constructed similar to the S-matrix calculations for linear metasurfaces.
In this paper, we generalize the concept of Smatrix to the description of nonlinear processes, such as frequency mixing, in the undepleted pump approximation [Eq.( 6)], schematically shown in Fig. 1.We derive the expression for the nonlinear S-matrix from Maxwell's equations [Eq.(A.14)] that shows how the S-matrix is connected to the incoming and outgoing fields and material susceptibilities.We demonstrate that S-matrix symmetry properties are connected to the microscopic and macroscopic symmetries of the metasurface [Eq.(7)].We further show that for degenerate sum-frequency generation processes such as optical harmonic generation, the exact multidimensional S-matrix can be replaced with a reduced two-dimensional matrix [Eq.(18)].We show that for metasurfaces possessing specific point symmetries [Eq.(20)], the selection rules determining the transformation of the reduced nonlinear S-matrix are simplified substantially and can be expressed via multiplication of a small number of two-dimensional matrices [Eq.(19)].
We apply the developed approach to analysis of selection rules for chiral harmonic generation in nonlinear metasurfaces excited in the normal incidence geometry.We analyse the fundamental and harmonic signals below the diffraction limit.We show that for metasurfaces with a rotational symmetry the harmonic signal generated by circularly polarized pump is allowed only in case of a specific relation between the harmonic order and symmetry order [Table 2].We demonstrate that this relation originates from the conservation of the component of the total angular momentum of light projected on the out-of-plane axis of the metasurface.We show that depending on the harmonic order the allowed harmonics in reflection and transmission can be either only co-or only cross-polarized relative to the incident beam.We further study how in-plane mirror symmetries of the metasurface affect the chiral harmonic generation.We show that presence of an in-plane mirror symmetry leads to zero nonlinear circular dichroism for all harmonic orders [Table 3].
We next study the nonlinear chiroptical effects in metasurfaces with and without out-of-plane mirror symmetry.We show that for suspended metasurfaces with up-down reflection symmetry, the nonlinear circular dichroism is nonzero despite the structure becomes geometrically achiral [Table 3].Finally, we derive the selection rules for the case of three-dimensional inversion symmetry [Table 3].For rotational and mirror refelction symmetries, we confirm the developed theory with full-wave numerical calculations.We believe the developed formalism paves a way towards a new generation of nonlinear metasurfaces which may find applications in active and nonlinear optics, biosensing and quantum information processing.

S-matrix for linear scattering
We briefly recall the basic properties and definitions used for the conventional S-matrix describing linear scattering in planar metasurfaces.The S-matrix provides a relation between incoming and outgoing channels of the metasurface, as shown in Fig. 1A.A channel is identified as a solution to Maxwell's equations within the medium external to the metasurface.We treat this medium as two half-spaces composed of uniform and Concept of nonlinear optical S-matrix and its reduction for degenerate nonlinear processes.A Schematic of the conventional S-matrix for linear optical scattering in metasurfaces.Columns ⃗ a (ω) and ⃗ b (ω) describe the field amplitudes in the input and output scattering channels, respectively.In case the metasurface has a specific point symmetry, the S-matrix satisfies Eq. ( 4).
B Top panel: Schematic of the nonlinear S-matrix for the specific case of sumfrequency generation with two input frequencies.We note that we assume the undepleted pump approximation.The nonlinear S-matrix is a matrix of the third rank in this case.We show that in case the metasurface has a specific point symmetry, the nonlinear S-matrix satisfies Eq. (7).Bottom panel: we introduce a concept of reduced nonlinear S-matrix described by a matrix of the second rank for degenerate sum-frequency generation, e.g. for (ω 1 = ω 2 ), independently on the order of harmonic.For specific symmetry transformations, its elements obey selection rules in Eq. (5).isotropic materials, separated by a flat interface.For such cases, Maxwell's equations provide complete sets of orthogonal plane waves functions that can be used to expand an arbitrary solution of Maxwell's equations outside of metasurface and on its surface.We provide explicit channel functions definition in Appendix A.
For the transversal waves propagating normal to the metasurface plane, the S-matrix Ŝ(ω) , schematically shown in Fig. 1A, can be defined as Here, i, j are the indices of the radiation channels, and a i and b i represent a set of incoming and outgoing amplitudes, respectively.Below the diffraction limit, the number of channels is limited to 4, including two orthogonal porarizations and two directions of incidence (from upper side and from bottom side).Thus, the S-matrix can be written as a (4 × 4) matrix of the second rank ( We assume that the metasurface possesses a specific point symmetry simultaneously at the macroscopic scale for the unit cell and meta-atom, and at the microscopic scale for the permittivity.We assume that the point symmetry is described by the operator T .
In the selected coordinate basis, the transformation is characterized by the matrix T ij , and The crystallographic symmetry principle (also called Neumann's principle) [46] imposes that the S-matrix must be invariant with respect to the symmetry operations of the point group of the metasurface unit cell.We can write that the corresponding operators commute Ŝ(ω) T = T Ŝ(ω) .Thus , Applying Eq. ( 4) for specific symmetry transformations T , one can find relations between different elements of Ŝ(ω) .This approach is widely used to analyse the conditions required for observation of strong chiroptical effects, e.g. in Refs.[15,17].

S-matrix for nonlinear sum-frequency generation
In this section, we aim to show that for specific symmetry transformations T and selected nonlinear processes of weak intensity in metasurfaces, we can introduce a reduced (4 × 4) nonlinear S-matrix, as schematically shown in Fig. 1B.The reduced matrix Ŝ(Nω) fully describes the nonlinear process and obeys simple selection rules resembling Eq. (4), We further derive Eq. ( 5) and show the limits of its applicability.We start with the general concept of S-matrix for nonlinear frequency mixing processes of weak intensity.We assume that the metasurface is pumped with N multiple frequencies ω 1 , . . ., ω N and generates a signal at the output frequency ω 1 + . . .+ ω N , as shown in the left panel of Fig. 2. We also assume the perturbative regime and undepleted pump approximation for the frequency conversion, which is widely used for dielectric and plasmonic metasurfaces due to low conversion efficiency.Then, the nonlinear S-matrix, schematically shown in the top panel of Fig. 1B, can be defined as b Here, the nonlinear S-matrix is a matrix of the nth rank, and i, j, . . ., k are the channel indices as in Eq. (1).
We can analyse the symmetry properties of Ŝ(ω1+...+ω N ) for metasurfaces characterized with a certain point group symmetry element T by combining Eqs. ( 3) and ( 6), resulting in Importantly, T is an element of the point group symmetry of the nonlinear S-matrix, only if it is included in the point symmetry group elements of the metasurface lattice, meta-atom, and material susceptibility tensors (Fig. 2).The proof is given in Appendix A. We also notice that the proposed nonlinear S-matrix is non-reciprocal, and does not obey energy conservation and time-reversal symmetry by definition, because we consider the undepleted pump approximation.One can see that the analysis of Eq. ( 7) is challenging due to a large dimensionality of the nonlinear S-matrix.Our goal is to show that for selected simple sum-frequency generation processes, such as generation of optical harmonics with a polarized pump beam, only specific elements of S (ω1+...+ω N ) ij...k are needed to describe the process, and Eq. ( 6) can be simplified to a form similar to Eq. ( 4).
Second-harmonic generation.We next consider the second-harmonic generation (SHG) process, for the sake of clarifity of the derivation.We note that the model developed below can be generalized for higher-order harmonic generation.For SHG with the pump and second-harmonic (SH) wavelength below the diffraction limit, the number of open radiative channels is 4, harmonic order N = 2, and Material susceptibility

Meta-atom
Figure 2. Nonlinear S-matrix for sum-frequency generation in metasurfaces.We consider that the metasurface is excited in the normal incidence geometry.For each input and output frequency, we consider four channels numbered by two orthogonal polarizations and two directions of incidence (only two channels per frequency are displayed), provided we analyse only the zero-order diffraction normal to the metasurface plane.The symmetry of nonlinear S-matrix is determined by the point symmetry achieved at the macroscopic (metasurface lattice point symmetry, meta-atom point symmetry) and microscopic (permittivity and nonlinear susceptibility tensor point symmetry) scales.
The nonlinear S-matrix Ŝ(2ω) for SHG is a matrix of the third rank with the shape of (4 × 4 × 4).Its definition can be written as [see Eq. ( 6)] We note that from here on we omit the frequency superscript for amplitudes a i and b i .In Eq. ( 9), the terms with S (2ω) ijk and S (2ω) ikj can be combined because of Eq. ( 8).Then, we can re-write Eq. ( 9) in a simpler form Here, we introduce an auxiliary S-matrix Ŝ(2ω) , which is a matrix of the second rank with the shape of (4×10) and elements defined as From here on, we use Greek symbols α, β, . . .for the indices ranging from 1 to 10, and Roman symbols i, j, k, . . .for the indices in the range from 1 to 4.
The correspondence between indexing with α and (j, k) is shown in Table 1.For example, S The S-matrix Ŝ(2ω) obeys the selection rules governed by [see Eq. ( 7)] We can re-write Eq. ( 12) in terms of transformation of Ŝ(2ω) by combining the corresponding elements of T , Here, T(2) is the auxiliary symmetry transformation operator with the matrix of (10 × 10) shape defined as The correspondence between β = 1 . . . 10 and p, s = 1 . . . 4 is the same as between α and (j, k) (see Table 1).For example, T (2) 78 = (T 12 T 43 + T 13 T 42 )/2.
We focus on incident beams with well-defined polarization, such as linearly-or circularly polarized plane waves.With a proper choice of coordinate basis for b and a matching the incident wave polarization, we can impose a i a j = 0 for (i ̸ = j) and further simplify Eq. (13).In this case, the input amplitude has a single nonzero component, so the output amplitudes b i are defined only by the elements S (2ω) iα with α = 1 . . .4, as can be seen from Eq. (10).We can then define this part of Ŝ(2ω) as a reduced nonlinear S-matrix Ŝ(2ω) The reduced S-matrix Ŝ(2ω) is a matrix of the second rank with the shape of (4 × 4) similar to the S-matrix Ŝ(ω) used for description of linear scattering.The selection rules for the reduced nonlinear S-matrix are defined via β,α(j,j) .( 16) The first term in Eq. ( 16) is defined via the reduced S-matrix only, resembling Eq. ( 4), while the second term cannot be expressed via Ŝ(2ω) .The transformation matrix in the second term can be written as T β(p,s),α(j,j) = T pj T sj with j, p, s = 1 . . . 4 and p ̸ = s due to β(p, s) = 5 . . .10.
We focus on specific symmetry transformations obeying T pj T sj = 0 for all p ̸ = s in the given coordinate basis.This means that the each column of the matrix T pj has no more than one nonzero element, as schematically shown in Fig. 3.Such transformations do not mix multiple channels, according to Eq. (3).As shown below in Sec. 3, this condition is valid for a large set of symmetries, especially, in the basis of circularly polarized waves with defined helicity.For this condition, the second term in Eq. ( 16) is zero, and, using T Optimal transformation More than one element � 0 Figure 3. Optimal and non-optimal transformation matrices.A General form of a non-optimal transformation matrix T .B General form of an optimal transformation matrix T .The optimal criteria are satisfied the transformation matrix does not mix multiple channels, i.e. if each column of T has only one non-zero element.The selection rules analysis is greatly simplified for optimal matrices, as Eq. ( 16) is transformed to Eq. ( 17) that contains two-dimensional matrices only.In the basis of circularly polarized waves with defined helicity, the optimal criteria are satisfied for all relevant point group symmetries, including the rotational, in-plane mirror, out-ofplane mirror and inversion symmetries.
Arbitrary order of harmonic.For harmonic generation of an arbitrary order N , we can introduce a reduced (4×4) nonlinear S-matrix in analogy to Eq. ( 15), Following the steps above (see Appendix B for details), one can show that Ŝ(Nω) obeys under specific symmetry transformations T that follow Equations ( 19) and ( 20) are the central result of our study.They represent the generalization of selection rules given by Eq. ( 4) for harmonic generation processes in the undepleted pump approximation.We note that for (N = 1) Eq. ( 19) transforms to Eq. ( 4), as expected.Importantly, instead of analysing the multidimensional matrix S (N ω) ij...k , the selection rules can be predicted from multiplication of three twodimensional (4 × 4) matrices, T −1 ij , (T ij ) N and S(Nω) ij .
3. Application to chiral harmonic generation and numerical validation

Model and definitions
We consider a metasurface shown schematically in the left panel of Fig. 2, located in the (x − y) plane.We assume exp(−iωt) time dependence for all fields.We consider that the nonlinear S-matrix is expressed in the basic of plane waves with defined helicity, i.e. polarized along the complex unit vectors For waves propagating along the positive z direction, e + and e − correspond to the right circular polarization (RCP) and left circular polarization (LCP), respectively.For waves propagating along the negative z direction, the sign of polarizations is flipped and e + and e − correspond to the LCP and RCP, respectively.The generation of the N -th order optical harmonic can be described by the reduced S-matrix equation relating the outgoing wave amplitudes b ±,u(d) at the harmonic frequency (N ω) with the incident wave amplitudes a ±,u(d) , Here, the parameters t ij,u(d) determine forward and backward scattered N -th harmonic signal, respectively, with i = L, R polarization of the outgoing waves for j = L, R polarized pump, and (u, d) define "up" (z < 0) and "down" (z > 0) direction of incoming wave.From here on, we omit the subscript ± for the reduced S-matrix for the sake of clarity.
We will focus on the key optical parameter describing chiroptical effects, circular dichroism (CD).We extend its definition to nonlinear harmonic generation.In general, due to nonlinear circular birefringence, an LCP wave at (ω) generates both LCP and RCP waves at (N ω).A The point group symmetry of the nonlinear S-matrix consists of identical elements of the point symmetry groups of rotational symmetry of the metasurface lattice, meta-atom, and material susceptibility tensors.B Schematic of the total angular momentum projection conservation.The total angular momentum σ inc z = N (in units of ℏ) of the incident beam composed of N right-polarized degenerate waves is summed up with the ∆σz arising from the interaction with the metasurface.∆σz is equal to an integer number of m because the metasurface has an angular period related to the rotational symmetry of m-th order Cm and conserves total angular quasi-momentum.The total angular momentum per one emitted photon for the output beam at harmonic frequency must be equal to ±1.If this condition could not be fulfilled, the chiral harmonic generation of specific order is forbidden.
We can define co-and cross-polarized nonlinear CD for p = u, d, CD (25) To calculate the matrices T± for specific symmetry transformations in the helical amplitude basis, we use a unitary transformation of corresponding matrices Txy defined in the Cartesian amplitude basis.The matrices Txy for various symmetry transformations are defined in Appendix C.

Rotational symmetry of m-th order
We assume that the structure has a rotational symmetry of m-th order (C m ) at both microscopic and macroscopic scale, as schematically shown for the case of C 3 in Fig. 4A.The matrix corresponding to C m transformation is where φ m = 2π/m.The matrix is unitary, thus T (φm) . We can write the selection rules Eq. ( 19) as Applying the selection rules to Eq. ( 27), we get four distinctive cases.First, if we can find an integer s 1 that for any integer s 2 we get the constraints for the coefficients, and the nonlinear S-matrix takes a simplified form, The CD calculation is simplified as Second, if we can find an integer s 1 that for any integer we get the constraints for the coefficients, RR,u(d) = 0, (33) and the nonlinear S-matrix takes a simplified form, The CD calculation is simplified as Third, if we can find integers s 1 and s 2 that we do not have any additional constraints.Finally, if for any integer s all S-matrix elements are zero Equations N + ms = ±1 represent the conservation law for the projection of total angular momentum of light on the z axis for chiral harmonic generation processes, shown schematically in Fig. 4B.We can interpret them as follows: N photons with the total momentum projection +1 excite the metasurface.Since the metasurface lattice, meta-atom and molecular lattice all possess rotational symmetry of m-th order, the incident photons can acquire a multiple integer number s of total momentum projection quanta m.The output photon at harmonic frequency can carry the total momentum projection ±1, thus the emitted photons should obey the conservation law.If these conditions cannot be fulfilled, chiral harmonic generation is prohibited due to inconsistency with the total angular momentum projection conservation.
We can apply this to analysis of specific harmonic processes.In the linear regime (N = 1), the crosspolarized conversion is prohibited for m > 2, thus incident LCP (RCP) light generates only LCP (RCP) harmonic signal in transmission.For SHG (N = 2), only m = 3 allows to fulfill the conservation law with s 1 = −1 and any s 2 .Only harmonic photons with the opposite polarization can be generated in transmission, as shown in Fig. 5A, e.g.incident LCP light generates only RCP transmitted harmonic signal and vise versa, in agreement with Eq. (33).For the  2. We confirm the developed theory with numerical calculations.More details about calculations can be found in Methods, see Appendix D. We first consider a dielectric metasurface made of the transition metal dichalcogenide MoS 2 in 3R phase, characterized with the C 3 symmetry of meta-atom, hexagonal metasurface lattice, and C 3 symmetry of the susceptibility tensor, shown in Fig. 5A.Therefore, for metasurface and susceptibility tensor lattices aligned, the symmetry of the nonlinear S-matrix is m = 3.We consider the structure pumped with a circularly polarized plane wave propagating along the out-ofplane axis.Figure 5B shows the linear transmission, SH signal, linear CD and nonlinear CD in the nearinfrared range of incident wavelengths.One can see that in the linear scattering, cross polarization conversion is prohibited, so RCP (LCP) light generates only RCP (LCP) transmitted signal.This agrees with the prediction made in the linear regime (N = 1) for m > 2 [Table 2].The linear CD is nonzero and is governed by resonances of the structure.In the forward generated SH signal, the co-polarized components are zero as predicted by the selection rules [Table 2].
The nonlinear SH CD is large and is governed by resonances.
Next, we consider a dielectric metasurface made of Si, characterized with the C 4 symmetry of meta-atom, square metasurface lattice, and cubic symmetry of the susceptibility tensor, shown in Fig. 5C.The symmetry of the nonlinear S-matrix is m = 4. Figure 5D shows the linear transmission, SH signal, linear CD and nonlinear CD in the near-infrared range of incident wavelengths.Similar to Fig. 5B, in linear transmission the conversion of polarization is prohibited, while for nonlinear signal it is opposite and co-polarized thirdharmonic (TH) signal is not generated in the forward direction.This confirms once again the prediction made by analysing the selection rules for the reduced nonlinear S-matrix.
These second-harmonic and third-harmonic examples give some insights that the reduced S-matrix with the simple selection rules can be helpful in designing a nonlinear chiral metasurface for experiments.Additionally, the knowledge of which circular polarization is expected to be transmitted or reflected at a given harmonic order will be significantly useful in constructing an optical measurement setup including different types of waveplates and polarizers.

In-plane mirror symmetry
We next analyse metasurfaces with an in-plane mirror symmetry, defined with respect to an inclined in-plane axis rotated by the angle θ from the x-axis, as shown in Fig. 6A.We note that θ has no direct physical meaning, because x, y axes can be redefined by rotating the coordinate system.The corresponding transformation matrix in the helical basis is The matrix is unitary, and T (θ) We can write the selection rules Eq. ( 19) as, For δ = 0 the meta-atom has L-shape with two equal arms of length a.The meta-atoms are arranged in a square lattice.The metasurface lattice has 4 distinct in-plane mirror symmetry planes.For δ ̸ = 0, 0.5, 1 the meta-atom does not have in-plane mirror symmetry planes, thus the metasurface as a whole does not have an in-plane mirror symmetry.C Calculation of linear and TH CD.The pump wave is incident from the upper side.The CD vanishes for linear and nonlinear TH signal for δ = 0, 0.5, 1 due to metasurface retaining one of in-plane mirror symmetries.Parameters of the simulation: meta-atom arm length a = 900 nm, height h = 400 nm, square lattice with period 1200 nm, substrate permittivity 2, material permittivity 16, nonlinear tensor χ (3) is defined by the space group 227 (e.g.Si).
Applying the selection rules to Eq. ( 40), we get Therefore, the nonlinear CD in all definitions reaches zero for an arbitrary N and θ, This is in agreement with the definition of geometrical chirality as absence of all mirror symmetries of the object.
The selection rules for different N are summarized in Table 3.We confirm the developed theory with numerical calculations, for details see Methods in Appendix D. We consider a dielectric Si metasurface characterized with the shift parameter δ for the meta-atom, square metasurface lattice, and cubic symmetry of the susceptibility tensor, shown in Fig. 6B. Figure 6C shows linear CD and nonlinear TH CD for pump in the near-infrared range.One can see that for parameter δ = 0, 0.5, 1, both linear and nonlinear CD are zero, as expected, because the meta-atom and metasurface lattice mirror symmetry planes coincide and the structure becomes geometrically achiral.

Out-of-plane mirror symmetry
We finally analyse metasurfaces with the out-of-plane (up-down) mirror symmetry (σ h ), defined with the matrix The matrix is unitary, and . We can write the selection rules Eq. 19 as We conclude from Eq. 44 that the selection rules are independent on the harmonic order N .Therefore, we can write that the transmission and reflection amplitudes are equal on both sides of the metasurface upon changing from the left-to right-polarization and vise versa Therefore, for all definitions of CD in Eq. (24), Linear scattering and reciprocity.For linear processes (N = 1), the fundamental principle of the Lorentz reciprocity adds additional constraints on the Smatrix.The selection rules imposed by the reciprocity in-plane out-of-plane 3D inversion CD (1)  co = CD (1)  cross = 0 CD (1)  co = 0 CD (1)  co = 0 CD (N )  co = 0 CD (N ) co ̸ = 0 CD (N ) co ̸ = 0 CD (N )  cross = 0 CD (N ) cross ̸ = 0 CD (N ) cross ̸ = 0 can be written in the helical basis as where T (R) ± is the corresponding transformation matrix (see Appendix C).The S-matrix coefficients are LL,u = t (48) Applying Eq. ( 48), we find that the co-polarized linear CD is zero, independently of the excitation side p = u, d, CD (1)  co = 0, CD (1) tot = CD (1)  cross .
(49) Thus, we can design structures with out-of-plane mirror symmetry, that demonstrate a zero linear CD reaching, but nonzero nonlinear CD, in case the inplane mirror symmetries are not present.The selection rules for different N are summarized in Table 3.
We confirm the developed theory with numerical calculations, for details see Methods in Appendix D. We consider the same dielectric Si metasurface as in Fig. 5C, with C 4 rotational symmetry.We compare the linear and nonlinear TH signal for the metasurface with and without the substrate, as shown in Fig. 7A.Figures 7B,C show that the linear CD in transmission vanishes in the absence of the substrate, as expected from the reciprocity and symmetry considerations, while the nonlinear CD is non-zero.

Other symmetries
Inversion symmetry.We assume that the structure has the three-dimensional inversion symmetry C i that represents a product of C 2 and σ h .
Then, the transformation matrix is The matrix is unitary, and T (I) (52) Thus, generation of both co-and cross-polarized chiral harmonics is allowed for circularly polarized incident beam, and Eq. ( 46) is fulfilled.For N = 1, Eqs. ( 48) and ( 49) are valid due to the reciprocity.The selection rules for different N are summarized in Table 3.

Conclusion and outlook
We have demonstrated that a large class of nonlinear frequency mixing processes in planar metasurfaces can be described within the framework of the scattering matrix approach similar to the case of the conventional linear scattering processes.We have introduced the definition and described properties of the nonlinear scattering matrix.We have shown that for optical harmonic generation the scattering matrix can be reduced to a simpler two-dimensional form.We have derived the simplified selection rules for the reduced scattering matrix in the case of metasurfaces obeying specific point group symmetry transformations, and have applied them to the analysis of chiral harmonic generation in nonlinear metasurfaces with rotational, mirror reflection and inversion symmetries, see Tables 2  and 3. We have verified our analytical results with numerical calculations.We note that the developed framework is frequency independent, as it relies on symmetries only, and can be applied to analysis of both resonant and non-resonant nonlinear response.Moreover, the method can be applied for nonlinear metasurfaces with misaligned microscopic and macroscopic lattices, as long as each of point symmetry groups (metaatom geometry, metasurface lattice arrangement and susceptibility tensor) includes the specific symmetry element.
We emphasize that the presented approach can be generalized to describe a variety of nonlinear frequency mixing processes, including but not limited to difference frequency generation, four-wave mixing, and self-induced nonlinear effects.The described formalism can be also generalized to the analysis of oblique incidence and analysis of diffraction effects with a larger number of scattering channels.
The presented nonlinear scattering matrix approach is not only theoretically insightful and robust but also has practical usage.Nonlinear chiral metasurfaces can be designed considering the selection rules at the desired harmonics prior to performing full-wave simulations and fabrication of samples.The material and meta-atom symmetries as well as the lattice structure and the substrate can be selected for a desired function of the metasurface.This provides a wellestablished engineering guideline strongly supported by solid theory.We believe our results can find many important applications for engineering nonlinear metasurfaces with the properties on demand for polarization control, filtering, and nonlinear wave-front shaping.
where k l = ω l /c, ε is the metasurface permittivity, ε bg is the background permittivity, E bg (ω l ) is the incident background field at the input frequency ω l that satisfies Here, ĜEE is the electric-electric Green function of the resonant metasurface at frequency ω l , defined as We are interested in the frequency of the input and nonlinear output beams below the diffraction order, thus we consider 4 input and 4 output channels.We define the S-matrix scattering channels as plane waves with specific polarization that carry power towards and away from the metasurface.As an example, we consider linearly polarized plane waves with polarizations along x and y axis, respectively.The input channels are at the input frequencies, and the output channels are at the output sum frequency.We define the electric fields in the channels as where Êx(y) are the unitary polarization vectors, the waves are propagating from the up (z < 0) and bottom (z > 0) direction, axis z is defined as in Fig. 4B, and H θ (±z) is the Heaviside theta function.
The magnetic fields can be reconstructed from the Maxwell equations.We assumed the fields in each channels are normalized such as they have a unity outgoing/incoming power, with their inner product defined via the Poynting flux.
We then define auxiliary functions E inc 1...4 (ω l ) as propagating waves composed of incoming and outgoing fields that are free of singularities at the input frequency ω l , [13] Any incident field can be expanded over input and output amplitudes as The auxiliary symmetry transformation operator T(N) matrix of (α max × α max ) shape can be defined as  20), thus the reduced nonlinear S-matrix approach cannot be applied for these transformations.
To calculate the corresponding matrices in the helical basis T± , we use matrices Txy and unitary transformation T± = Tc Txy For numerical simulations of the transmittance and harmonic spectra, we use the finite-element-method solver in COMSOL Multiphysics in the frequency domain.All calculations were realized for a metasurface placed on a semi-infinite substrate surrounded by a perfectly matched layer mimicking an infinite region.The simulation area is the unit cell extended to an infinite metasurface by using the Bloch boundary conditions.The material properties are listed in the corresponding figure caption.The incident field is a plane wave in the normal excitation geometry with left or right circular polarization incident from the upper side of the metasurface.The nonlinear simulations of chiral harmonic response are employed within the undepleted pump approximation, using two steps to calculate the intensity of the radiated nonlinear signal.First, we simulated the linear scattering at the pump wavelength, and then obtained the nonlinear polarization induced inside the meta-atom.Then, we employed the polarization as a source for the electromagnetic simulation at the harmonic wavelength to obtain the generated harmonic field.The nonlinear susceptibility χ(N) was considered as a constant tensor corresponding to the specific crystallographic point group listed in the corresponding figure caption.
Figure 1.Concept of nonlinear optical S-matrix and its reduction for degenerate nonlinear processes.A Schematic of the conventional S-matrix for linear optical scattering in metasurfaces.Columns ⃗ a(ω) and ⃗ b (ω) describe the field amplitudes in the input and output scattering channels, respectively.In case the metasurface has a specific point symmetry, the S-matrix satisfies Eq. (4).B Top panel: Schematic of the nonlinear S-matrix for the specific case of sumfrequency generation with two input frequencies.We note that we assume the undepleted pump approximation.The nonlinear S-matrix is a matrix of the third rank in this case.We show that in case the metasurface has a specific point symmetry, the nonlinear S-matrix satisfies Eq.(7).Bottom panel: we introduce a concept of reduced nonlinear S-matrix described by a matrix of the second rank for degenerate sum-frequency generation, e.g. for (ω 1 = ω 2 ), independently on the order of harmonic.For specific symmetry transformations, its elements obey selection rules in Eq. (5).

2 .
We note the number of columns of S (2ω) iα is defined as α max = d(d+1)/2, where d is the number of open channels (see more details in Appendix B).In our case, d = 4 and α max = 10.

Figure 4 .
Figure 4.General concept for metasurfaces with rotational symmetries.A The point group symmetry of the nonlinear S-matrix consists of identical elements of the point symmetry groups of rotational symmetry of the metasurface lattice, meta-atom, and material susceptibility tensors.B Schematic of the total angular momentum projection conservation.The total angular momentum σ inc

Figure 5 .
Figure 5. Numerical analysis for metasurfaces with rotational symmetries.A,C Schematic of linear and nonlinear transmission in metasurfaces with C 3 meta-atoms of C 3 material in a hexagonal lattice A and for C 4 meta-atoms of C 4 material in a square lattice C. B,D Calculated linear and nonlinear transmission signals and CD.The pump wave is incident from the upper side, index p = u is omitted.Parameters of the simulations for the metasurface in A, B: L = 530 nm, D = 180 nm, W = 430 nm, H = 640 nm, period 1440 nm, substrate permittivity 3.13, material permittivity 16.7 at pump and 22.5 at SH wavelength, respectively, χ(2) is defined by the point group 3m (e.g.MoS 2 ).Parameters of the simulations for the metasurface in C, D: L = 550 nm, D = 250 nm, W = 230 nm, H = 400 nm, period 1300 nm, substrate permittivity 2, material permittivity 16 at pump and 16 at TH wavelength, respectively, χ(3) is defined by the space group 227 (e.g.Si).

Figure 6 .
Figure 6.Metasurfaces with in-plane mirror symmetries.A Definition of in-plane mirror symmetry plane (purple) and relative angle θ.B Schematic of metasurface inverse T-shaped meta-atom with shift parameter δ controlling the horizontal position of the upper bar.For δ = 0 the meta-atom has L-shape with two equal arms of length a.The meta-atoms are arranged in a square lattice.The metasurface lattice has 4 distinct in-plane mirror symmetry planes.For δ ̸ = 0, 0.5, 1 the meta-atom does not have in-plane mirror symmetry planes, thus the metasurface as a whole does not have an in-plane mirror symmetry.C Calculation of linear and TH CD.The pump wave is incident from the upper side.The CD vanishes for linear and nonlinear TH signal for δ = 0, 0.5, 1 due to metasurface retaining one of in-plane mirror symmetries.Parameters of the simulation: meta-atom arm length a = 900 nm, height h = 400 nm, square lattice with period 1200 nm, substrate permittivity 2, material permittivity 16, nonlinear tensor χ(3) is defined by the space group 227 (e.g.Si).

Figure 7 .
Figure 7. Metasurfaces with out-of plane mirror symmetry.A Illustration of the vertical mirror symmetry σ h plane and symmetry breaking scenario.(B,C) Linear and TH nonlinear CD for a Si metasurface with C 4h meta-atoms in a square lattice in the near-infrared range.The pump wave is incident from the upper side.Parameters of the structure are the same as in the Fig. 5 C. B Linear CD is prohibited by the combination of up-down symmetry and reciprocity.C Both linear and nonlinear CD are allowed by symmetry.

. (A. 7 )
l )E in i + a out i (ω l )E out iAt the same time, from Eq. (A.6)E bg = 4 i=1 a inc i (ω l )E inc i .(A.8)Thus, a inc i = a in i , and we can define the incoming amplitude for nonlinear S-matrix asa i (ω l ) ≡ a in i (ω l ), i = 1 . . .4.(A.9)Nω) is a matrix of the second rank with the shape of (d × α max ), α(j, k, . . ., l) = 1, . . ., α max , and d is the number of radiation channels (we used d = 4 in the main text).The number of columns α max can be calculated as α max = d(d + 1) . . .(d + N − 1) N ! .(B.2)

TT 1 )
jp T ks + . . .+ T js T kp N ! .(B.3)Appendix C. Transformation matrices in Cartesian basisIn this section we provide transformation matrices Txy in the in the Cartesian basis of amplitudes a x,y,u(d) and b x,y,u(d) for waves travelling in positive and negative direction along the z axis.The matrices can be constructed from analysis Eq. (3).The matrix corresponding to the rotational symmetry of m-th order (C m ) isT (φm)The matrix for the in-plane mirror symmetry, defined with respect to an inclined axis rotated by the angle θ from the x-axis, is T The matrix corresponding to out-of-plane up-down mirror symmetry (σ h ) isT (σ h )

5 ) 8 )
The selection rules for the linear S-matrix Ŝ(ω) xy imposed by the Lorentz reciprocity can be written in Cartezian basis as Ŝ(ω) xy = Ŝ(ω) xy T .(C.6)To transform them into the helical basis, we apply the unitary transformation to the S-matrix, Appendix D. Methods: numerical calculations.
S 11 S 12 S 13 S 14 S 12 S 22 S 23 S 24 S 13 S 23 S 33 S 34 S 14 S 24 S 34 S 44

Table 2 .
Selection rules for the reduced nonlinear Smatrix elements for chiral harmonic generation for structures with symmetry Cm.RR , t LL , r LR , r RL are allowed, only t RL , t LR , r LL , r RR are allowed, all elements are zero.

Table 3 .
Selection rules for CD for structures with in-plane mirror, out-of-plane mirror, and three-dimensional inversion symmetries.