Crossing of the branch cut: the topological origin of a universal 2{\pi}-phase retardation in non-Hermitian metasurfaces

Full wavefront control by photonic components requires that the spatial phase modulation on an incoming optical beam ranges from 0 to 2{\pi}. Because of their radiative coupling to the environment, all optical components are intrinsically non-Hermitian systems, often described by reflection and transmission matrices with complex eigenfrequencies. Here, we show that Parity-Time symmetry breaking -- either explicit or spontaneous -- moves the position of Zero singularities of the reflection or transmission matrices from the real axis to the upper part of the complex frequency plane. A universal 0 to 2{\pi}-phase gradient of an output channel as a function of the real frequency excitation is thus realized whenever the discontinuity branch bridging a Zero and a Pole, i.e a pair of singularities, is crossing the real axis. This basic understanding is applied to engineer electromagnetic fields at interfaces, including, but not limited to, metasurfaces. Non-Hermitian topological features associated with exceptional degeneracies or branch cut crossing are shown to play a surprisingly pivotal role in the design of resonant photonic systems.


Abstract
Full wavefront control by photonic components requires that the spatial phase modulation on an incoming optical beam ranges from 0 to 2π. Because of their radiative coupling to the environment, all optical components are intrinsically non-Hermitian systems, often described by reflection and transmission matrices with complex eigenfrequencies. Here, we show that Parity-Time symmetry breaking -either explicit or spontaneous-moves the position of Zero singularities of the reflection or transmission matrices from the real axis to the upper part of the complex frequency plane. A universal 0 to 2π-phase gradient of an output channel as a function of the real frequency excitation is thus realized whenever the discontinuity branch bridging a Zero and a Pole, i.e. a pair of singularities, is crossing the real axis. This basic understanding is applied to engineer electromagnetic fields at interfaces, including, but not limited to, metasurfaces. Non-Hermitian topological features associated with exceptional degeneracies or branch cut crossing are shown to play a surprisingly pivotal role in the design of resonant photonic systems.

I. INTRODUCTION
Over the last decades, the development of artificial materials and artificial interfaces to address arbitrarily output-from-inputs signals has considerably modernized the fields of optics and optical design. A particularly impressive amount of research in the field of metamaterials has been devoted to the design and realization of ultra-thin artificial optical surfaces for wavefront engineering and control. These optical surfaces, also dubbed metasurfaces, rely on the coherent scattering of light by a sizable distribution of nanoscatterers of various shapes and material compositions.
The list of optical effects achieved using metasurfaces is extensive, ranging from anomalous reflection and refraction [1][2][3][4][5][6][7][8], all the way to the design of utterly complex vectorial holographic surfaces [9][10][11][12][13][14][15]. The interested reader could refer to reviews and references therein for detailed descriptions and applications [16][17][18][19]. In practice, light control is achieved by varying the geometries of adjacent elements so as to provide spatially-varying phase retardation on the incoming wavefront. Surprisingly, among all physical mechanisms of interest for the design of these phase building blocks, including the Pancharatnam-Berry phase in anisotropic nanoparticles [20,21] and effective index waveguide modes in nanopillars [2], the fundamental origin of the 2π-phase modulation associated with resonant scattering of Mie nanoparticles has not been deeply investigated.
Ultimately, and despite all the efforts in understanding which physical mechanisms are leading to optimal designs, the most advanced patterns often require witless numerical parameter searches.
Here, we discovered that full 2π-phase retardation commonly used for designing metasurfaces exclusively relies on the relative spectral positions of topological singularities of the metasurface response functions. These manifest as Zero-Pole pairs in the complex frequency plane. Resonantphase metasurfaces are designed to connect N input E in -modes with N output E out -modes. They behave as regular linear open wave systems mathematically described by N × N response matrices F (ω) such that, E out = F (ω)·E in [22,23], where F(ω) represents either the S-matrix, the reflection matrix R or the transmission matrix T [24] . Zeros and Poles of response matrices or functions behave as phase singularities around which the phase is spiraling in a vortex-like way [25,26].
Metasurfaces, like most optical devices, are intrinsically non-Hermitian devices, suffering from scattering losses and potentially from intrinsic loss or gain. Therefore, their phase singularities generally occur in the complex frequency plane. Singularities have long been known to play an important role in optics [27][28][29][30]. Here, we show how phase singularities occurring in the complex plane greatly influence the optical response of resonant metasurfaces. More precisely, we show that 2π phase accumulation needed for wavefront engineering of reflected or transmitted output channel has a deep topological origin. It requires two complementary singularities, known as Poles and Zeros, to be disposed in the complex frequency plane on either side of the real axis. This way, the branch-cut discontinuity bridging these two singularities crosses the real axis, providing 2π phase accumulation as a function of the real frequency excitation. We demonstrate that one can rely on symmetry-breaking in order to have two singularities of the same pair located on each side of the real axis which could induce 2π circulation for any path in the complex plane encircling a singularity. We propose two approaches inducing Parity-time symmetry breaking to move the Zeros from the real axis to the upper part of the complex plane.
The first solution is relatively straightforward and relies on explicit symmetry-breaking arguments. The second solution, which is linked to the regime known as Huygens metasurfaces [31], is more subtle as it occurs in symmetric systems featuring spontaneously symmetry-broken states.
We thus show that the physical origin of Huygens metasurfaces is deeply rooted in topological concepts. Finally, looking at the analytical expressions of metasurface boundary conditions, namely the "generalized sheet transition conditions" (GSTC), we identified the symmetry conditions of the electromagnetic modes to be considered for spontaneous symmetry-breaking.

II. COMPLEX FREQUENCY ANALYSIS AND PHASE INTEGRATION
We begin our analysis by studying the analytical properties of a non-Hermitian metasurface represented by a response function F (ω), e.g. its reflection R or its transmission T coefficients.
Expanding F (ω) on the basis of its Poles ω p,n and Zeros ω z,n using Weierstrass product expansion, we have [25,26]: where A and B are constants depending on the properties of the metasurface and the considered response function. By definition, the phase is given by Arg(F ) = Im (log(F )) and, as a consequence, Zeros and Poles of F correspond to logarithmic branch points at which the phase becomes singular [32]. The phase is thus spiraling as vortices around Zeros and Poles of F in the same way as the phase of a complex number in C rotates around the origin of the complex plane [25]. Zeros and Poles are thus topological defects [33]. While the Poles of S, R, and T are identical and are associated with resonances or quasi-normal modes solutions of the wave equation, their complex Zeros are different and provide insights into the physical effects governing the metasurface response.
In particular, they greatly influence the real-frequency response of metasurfaces. Zeros of the Smatrix are complex conjugate of Poles for passive lossless metasurfaces and can reach the real axis in lossy structures where they are linked to perfect coherent absorption [34] . Zeros of R and T coefficients are linked to Reflectionless and Transmissionless states, respectively [35] . Let us consider a toy-model F (ω) = (ω − ω z )/(ω − ω p ), derived from Eq. (1) for A = 1 and B = 0 and considering a single Zero-Pole pair, featuring only one Zero ω z and one Pole ω p located on each side of the real axis, as indicated by the blue and red regions respectively in Fig. 1(a.). We define as topological charge the quantity q calculated by evaluating the winding number along a counterclockwise contour encircling the singularities as [33,36]: equalling to q = +1 for a Zero and q = −1 for a Pole (see SI). The phase varies by 2π around the Zeros and the Poles with an opposite sign. As a result, the phase variation around both is null, see Fig. 1(b.). Moreover, Zeros and Poles occur in pairs and are connected by a phase discontinuity, also called branch cut. This can be seen in the panel in Fig. 1(b.) for a single Zero-Pole pair. In our example, the complex Zero and Pole values are purposely chosen such that the branch cut is crossing the real axis. The sufficient condition for a metasurface to support a 2π resonant phase accumulation as a function of the real frequency is thus to possess a Zero-Pole pair separated by the real axis, as illustrated in Fig. 1(c.). After unwrapping the phase discontinuity resulting from crossing the branch cut at ω = 1, we obtain a 2π phase difference ∆Arg (F ) = Arg (F (ω 2 )) − Arg (F (ω 1 )) as ω 1 → 0 and ω 2 → 2 (see Fig. 1(c)). Cauchy's residue theorem is used to link the 2π phase accumulation on the real axis with the existence of a phase singularity in the upper part of the complex plane. To integrate ∆Arg (F ) = ω 2 ω 1 dArg(F (ω)) dω dω in the complex plane, we consider the integral along a contour C l containing one singularity of charge q such as the one displayed in Fig. 1(b.). It consists of the line segment [ω 1 , ω 2 ] closed by a semi-circle C sc of radius R in the upper part of the complex plane. As a consequence, dω. However, it can be shown that Csc dArg(F (ω)) dω dω → 0 when the radius R increases, see SI for the exact case of the transmission coefficient. Consequently, ∆Arg (F ) → 2πq directly linking the phase accumulation on the real axis to the existence of a phase singularity above the real axis.

GULARITIES IN THE COMPLEX PLANE
As a consequence of causality, Poles are always restricted to the lower part of the complex plane in passive systems [37] (the convention used for the time dependence throughout this article is e −iωt ). It is known that in active systems, adding gain can in principle move Poles upwards, fulfilling thereby the lasing condition when they reach the real axis, but this is not a simple and technologically relevant solution for metasurfaces. This is, therefore, outside of the scope of the present study. A question then arises: how to design a passive structure that can support a phase singularity in the upper part of the complex plane thus displaying a 2π resonant phase accumulation over the real axis? Discarding active systems with complex Pole manipulation, we propose to move Zeros to the upper part of the complex plane for the reflection or the transmission coefficients. Let us point out that for lossless structures or at least for structures featuring a small amount of loss, the Zero-Pole pairs of the S-matrix coefficients are usually separated by the real axis. Therefore, the S-matrix coefficients display a 2π phase accumulation for each of its Zero-Pole pairs. This can be seen in [26] for a plasmonic dolmen metasurface. This is however of little practical interest for designing metasurfaces since the S-matrix links purely incoming fields to purely outgoing fields that are difficult to produce experimentally. In the remaining of this article, we will consequently focus on ways to move transmission and reflection Zeros to the upper half of the complex plane.
To this end, one first needs to better understand the constraints imposed on the locations of the Zeros in the complex plane by the symmetries of the metasurface.
Transmissionless states of a system supporting Time-reversal symmetry. Any dielectric c.  [23]. And the same is true for Reflectionless states.
Reflectionless states of a system supporting Parity-Time symmetry.
The study of the underlying symmetry properties of Reflectionless states necessitates considering the reflection symmetry with respect to the median plane (z=0) of the metasurface in addition to time reversal [22]. Let us consider a metasurface possessing a reflection zero for light impinging from one given side at a possibly complex-valued frequency ω RZ . If the metasurface is invariant under time-reversal, then it also possesses a reflection zero at ω * RZ for light impinging on the other side as illustrated in Fig. 1(d.). In addition, the metasurface may remain invariant under the symmetry transformation which consists of a reflection with respect to its median plane. For the latter condition to be fulfilled, not only should the permittivity of the substrate sub be equal to the permittivity of the top medium sup but the scatterers composing the metasurface should all remain invariant under this reflection transformation as well (which is the case for the nanocubes considered in the following of the article). A metasurface invariant under time reversal and under the latter reflection symmetry (that will be called parity symmetry in the following) possesses a reflection zero at both frequencies ω RZ and ω * RZ for light impinging from the same side as shown Fig. 1

(d.).
Real versus Complex-conjugate Zeros of symmetry-preserving systems. These considerations indicate that the time-reversed counterpart of a Transmissionless state is also a Transmissionless state, and similarly, the parity-time symmetric of a Reflectionless state is also a Reflectionless state, as illustrated in Fig. 1(d.). Upon conservation of symmetries, transmission or respectively reflection Zeros should thus occur either at real frequencies (ω TZ is real) or as conjugated pairs in the complex plane (i.e. ω TZ and its complex-conjugated counterpart are located symmetrically with respect to the real axis). The former case is generic and could happen even for a single Zero-Pole pair, but to understand the latter case it is necessary to remind that Eq. (1) associates one Zero to only one Pole. Therefore, symmetry-preserving systems having conjugated Zeros positioned symmetrically with respect to the real axis imply necessarily the interaction of at least two Zero-Pole pairs.

ACCUMULATION ON THE REAL AXIS
Given these symmetry constraints, two approaches can be considered to expel Zeros to the upper part of the complex plane: (i) Explicit breaking of a relevant symmetry, which relaxes the real-value constraint for an isolated Zero, thus allowing to move it into the complex plane.
This solution is quite straightforward and works both for reflection, and respectively transmission, depending if T -symmetry, or respectively P T -symmetry, breaking is engineered, for instance by introducing gain or losses or by putting the metasurface on a substrate. (ii) Avoided crossing between two Zero-Pole pairs: Considering two pairs, the system could remain T -or P T -symmetric as a whole, but decomposed into two states which are not symmetric anymore.
Spontaneous symmetry-breaking of the Zero states thus creates the condition for which two different Zeros bifurcate as a conjugated pair in the complex plane. The bifurcation of the once-coalesced Zeros forces them to move symmetrically into the upper and lower halves of the complex plane.
As the overall symmetries are preserved, the conjugated Zeros are positioned symmetrically with respect to the real axis. The trajectories of the Zero-Pole pairs as a whole, presented in the supplementary videos (without or with absorption losses), indicate that the trajectories of the two pairs are such that the branch cuts linking each pair avoid each other. To this end, one of the Zeros has to go down in the complex plane while the Pole from the other pair moves up so that the branch cuts from each pair do not cross. This will be further discussed in section V [38] .
As a result, both scenarios leverage on Zeros expelled to the upper part of the complex plane to achieve a 2π-phase accumulation caused by the branch cut of one zero-pole pair crossing the real axis. We illustrate both scenarios considering a metasurface composed of square nanocubes with a length L = 350 nm and various heights h, made of a dielectric medium with relative permittivity ε = 8.05 (which approximately corresponds to Sb 2 S 3 phase-change material in the amorphous state) arranged in a 2D square array with a period p = 500 nm. The substrate relative permittivity is equal to ε sub = 2.25 while the relative permittivity of the top medium will be varied in the following examples. The numerical calculations of the optical response of these metasurfaces for complexvalued frequencies are performed using JCMsuite, a software based on the Finite-element method [39]. For details on the numerical simulations, please see the Supporting Information section SIV.

A. Explicit Parity-time symmetry-breaking for phase engineering in reflection
To illustrate the first method, we start by considering a PT-symmetric version of the Sb 2 S 3 metasurface: it is lossless and the superstrate is a dielectric material with the same relative permittivity as that of the substrate: ε sub = ε sub = 2.25, making the overall system invariant under reflection with respect to its median plane at z = 0. For nanocubes of height 160 nm, it can be seen in Fig. 2 (a.) that a reflection Zero exists on the real axis for ω ≈ 2.15 × 10 15 rad · s −1 . The parity symmetry may be broken simply by changing the dielectric permittivity of the top material to ε sup = 2.12, keeping the substrate permittivity unchanged and preserving the time-reversal symmetry. The results for an incident field impinging from the substrate side are shown in Fig.   2 (b.). As a consequence of the parity symmetry breaking, the reflection Zero moves upwards. A drastic change of the behavior of the reflection coefficient phase in the ω-plane can thus be observed. The plot of the phase dependency for real frequencies in Fig. 2 (a.) shows that there is a phase discontinuity when the Zero is located on the real axis. On the other hand, the phase varies from 0 to beyond 2π when the Zero moves to the upper part of the complex plane as seen in Fig. 2 (b.). The reason why the phase accumulation exceeds 2π is the existence of another Zero-Pole pair in the lower part of the complex plane for larger frequencies (see SI). Transposing this to the transmission case, i.e. moving up the Transmissionless states, requires breaking the time-reversal symmetry by adding loss or gain in the metasurface [40].

B. Spontaneous symmetry-breaking at EP for phase engineering in transmission
Let us now illustrate the second approach. The map of the phase of the transmission coefficient of the Sb 2 S 3 metasurface for nanocube height equal to 190 nm is displayed in Fig. 2 Fig. 3 (a.) and (b.). These maps reveal a behavior identical to the one reported for dielectric Huygens metasurfaces (HMS) [31] with the existence of a range of heights for which the amplitude of the transmission coefficient becomes large and the phase covers the full 2π interval. The real and imaginary parts of transmission Zeros as a function of the height of the metasurface nanocubes are shown in Fig. 3 (b.) and (c.). Outside the HMS operating region, the two Zeros occur for two different real parts but share the same imaginary part (which vanishes). However, the two Zero real EPs are a type of degeneracies specific to non-Hermitian systems [43,44]. At the degeneracy, not only do the eigenfrequencies become identical but their associated eigenstates coalesce. In the field of metasurfaces, the coalescence of eigenstates at an EP has been used for the design of polarization-dependent metasurfaces recently [45]. EP have also been proven useful in many other applications [46], in particular for sensing [47,48]  tuned to obtain these EPs, while these are usually found by tuning two parameters. This is a consequence of the underlying symmetries of the system as pointed out in [49]. An important consequence of the discovery of this mechanism is that it provides a fundamental explanation of the physics underlying the design of widely used Dielectric Huygens metasurfaces (HMS) [31,[50][51][52][53][54][55][56][57][58][59][60][61][62][63][64][65][66].
The original understanding of HMS properties relied exclusively on a spectral overlap of electric and magnetic dipolar resonances, which was derived following a theoretical modal analysis of arrays of spherical silicon particles [67]. By carefully studying the trajectories of Zeros, we prove that the boundaries of the 2π phase change regions correspond to EPs where two Transmissionless Zeros become degenerated. We also computed the electric and magnetic resonance frequencies of this Sb 2 S 3 metasurface, and the silicon nanodisk metasurface proposed in [31] (see SI), and we found that they never fully coincide. The 2π phase accumulation of the Huygens metasurfaces, originally attributed to the superposition of π phase of two coinciding resonances [31], is the consequence of the discovered 2π-phase mechanism reported herein, i.e. only one Zero-Pole branch-cut crossing the real axis.
Finally, it is worth mentioning that the symmetry considerations used throughout this article were determined for a field impinging on the metasurface at normal incidence. At non-normal incidence angles, additional spatial symmetries of the metasurface play an important role, in particular inplane symmetries, explaining the high sensitivity of Huygens metasurfaces response with respect to the incident angle [68].  expressed in terms of effective surface susceptibilities [69]. In this way, we will be able to link our previous results obtained using symmetry-based discussions to the electric and magnetic dipole responses of metasurfaces which are widely used as a design tool. The effective metasurface susceptibilities provide simple expressions for the reflection R or transmission T coefficients. For a normally propagating incident wave, they read as [69]: where χ e and χ m are the metasurface electric and magnetic isotropic susceptibilities and k = ω/c with c being the speed of light in the metasurface surrounding medium. It is important to notice that magnetic currents generate odd modes and electric currents generate even modes as seen in Fig.   4(a). Here, even and odd are related to the parity of the solutions with respect to the median plane of the structure z = 0. Because of modal symmetry, previous works on HMS rightfully explained that destructive interference between these modes on one or the other side of the metasurface may lead to the appearance of reflection and transmission Zeros, also related to the first and second Kerker conditions [31]. Here, we particularly focus on Transmissionless states with zero optical where A e/m is the amplitude, ω e/m is the resonance frequency and γ e/m is the damping factor.  Fig. 4 (e) shows the same configuration as previously but with additional damping in the GSTC formula (γ e = 0, γ m = 0). In this case, time-reversal symmetry is broken. It can then be noted that there is no EPs of transmission Zeros. In fact, it is well known that the tuning of two parameters is usually necessary for reaching an EP. The reason why only one was necessary in the case of the lossless metasurface considered previously was the underlying symmetry of the system [49], in our case time-reversal symmetry. As soon as this symmetry is broken, two parameters are again needed for reaching an EP. In Fig. 4 (e), the Zeros thus do not merge, but the avoided crossing behavior remains, leading again to the existence of one Zero in the upper half of the complex plane. As a result, the 2π-phase accumulation is still observed without Zeros degeneracy and EP, indicating that this approach is robust and would also apply to metasurfaces made of arbitrary materials.

VI. CONCLUSION
In summary, we have proposed a general method to address 2π phase accumulation of the reflected or transmitted light at interfaces. We discovered that a sufficient condition is to have  [2] P. Lalanne, S. Astilean, P. Chavel, E. Cambril, and H. Launois, Blazed binary subwavelength gratings with efficiencies larger than those of conventionaléchelette gratings, Optics Letters 23, 1081 (1998).

SII. TOPOLOGICAL CHARGE AND RESIDUE THEOREM:
Let F (ω) be a response function characterizing the optical response of a metasurface. In general, If ω z,i is a Zero of order q of F (ω), then F (ω) can be well approximated by the function (ω − ω z,i ) q g(ω), g(ω) being a regular function with no Zeros and Poles, for ω belonging to the vicinity of ω z,i . log (F (ω)) can then be approximated by the following expression for ω belonging to the vicinity of ω z,i : The residue of d log(F (ω)) dω associated with its Pole at ω z,i is by definition equal to lim ω→ω z,i (ω − Eq. (S2) allows us to show that this residue is equal to q. If now, a Pole ω p,i of F (ω) of order q, analogous derivations allow us to show that its residue for d log(F (ω)) dω is equal to −q.
Several results can be derived using Cauchy's residue theorem. Now if we consider a contour containing several Zeros ω z,i of order q z,i and Poles ω p,i of order q p,i , the residue theorem leads to the following result: As a result, the contour integral on C d log(F (ω)) dω dω vanishes when C contains both a simple Pole and a simple Zero for example. Eq. (S3) is more general.
Let us emphasize that the Pole expansion of the logarithmic derivative of a response function is the starting point for deriving the Weierstrass product expansion provided in Eq.

SIII. PHASE ACCUMULATION AND GAUSS THEOREM:
Here, we are providing some additional derivations allowing us to illustrate the link between positions of topological charges and the phase accumulation on the real axis. To this end, we consider the type of contours shown in Fig. 1 of the main text: a contour made of a segment on the real axis plus a semi-circle in the upper half of the complex plane.
We will use derivations for the transmission coefficient to illustrate these results. Since these derivations require determining the logarithmic derivative of the transmission coefficient, we will use here the analytical formulas provided by the GSTC method (see Eqs. (4a) and (5)  We want to compute the phase difference between two frequencies ω 2 and ω 1 : ∆Arg(T ) = Arg(T (ω 2 )) − Arg(T (ω 1 )). For the sake of simplicity, we choose ω 2 and ω 1 that are symmetric with respect to a center frequency ω 0 : ω 1 = ω 0 − R and ω 2 = ω 0 + R. Computing this phase difference can be done by performing the following integration: ∆Arg(T ) = dω dω . The evaluation of the latter integral will be done via a deformation of the contour in the complex plane.
Let us then consider the following integral: C d log(T (ω)) dω dω where C is the contour shown in Fig. S1(a.). Since there is only one Zero within this contour, this integral is equal to 2iπ (see results in the previous section). The contour is split up into a segment on the real axis [ω 0 − R; ω 0 + R], where ω 0 = 2.159 × 10 15 and R is varied, and the semi circle of radius R in the upper complex plane centered on ω 0 : Consequently, in order to have Arg(T ) → 2π, we need to show that Im The behavior of the semi-circle contour term on the right hand side of Eq. (S4) as R increases is shown in Fig. S1(b.). The discontinuity corresponds the radius at which the Zero enters the contour. Then, for larger radii, we see clearly that this contour term decreases towards 0 and eventually vanishes when R is of the order of 10 15 . As a consequence, the integral on the real axis in the right-hand side of Eq. (S4) tends toward 2π as seen in Fig. S1. This result allows to directly link the behavior of the phase on the real axis to the existence of a phase singularity in the upper half of the complex plane.

SIV. REAL AND COMPLEX FREQUENCY NUMERICAL CALCULATION
Results from Fig. 2(a.)-(d.) and Fig. 3(d.) were obtained using the Maxwell solver JCMsuite based on a higher-order finite element method [39]. It allows calculating a complex transmission coefficient for complex excitation frequencies from which we extract transmission amplitude and phase.
The trajectories of the transmission Zeros shown in Fig. 3 Here we show the trajectories of the electric and magnetic eigenfrequencies for the Sb 2 S 3 metasurface used in the main text and the silicon metasurface proposed in [31]. These results were derived by solving an eigenvalue problem where outgoing boundary conditions were imposed in the regions at the top and bottom of the metasurface. These calculations were performed using the eigenvalue solver implemented in JCMsuite [71].
We note that, in both cases, the two Poles never fully coincide. For the Sb 2 S 3 metasurface S2 (a.) and (b.), there is a crossing of the real parts of the electric and magnetic resonance frequencies in the Huygens metasurface operating region (between the dashed lines) but their imaginary parts never cross. Similarly, the resonances of the silicon nanodisk metasurface never fully coincide as can be seen in Fig. S2 (c.) and (d.).

SVI. ADDITIONAL RESULTS METASURFACE WORKING IN REFLECTION:
A. Second Zero-Pole pair contributing to the phase accumulation in Fig. 2 In Figs 2 (a.) and (b.), it was shown how the explicit parity symmetry breaking could be used to design a metasurface displaying a 2π resonant phase accumulation over the real axis. It was in fact noted that the phase accumulation is larger than 2π in Fig. 2 (b.). This is due to the existence of an additional Zero-Pole pair at larger frequencies as can be seen in Fig S3. Given that this second Zero remains under the real axis, it doesn't lead to an additional 2π accumulation on the real axis but it adds instead about π to the phase accumulation seen in Fig. 2 (b.). B. Time-reversed reflection zero in the structure considered in Fig. 2 (b.) Fig. 2 (b.) shows the existence of the reflection zero for a frequency ω RZ located in the upper half of the complex plane for an asymmetric structure with the permittivity of the substrate sub different from the permittivity of the top medium sup when light is impinging from the substrate. Time-reversal symmetry imposes the existence of a reflection zero at ω * RZ for light impinging from the top medium. In Fig. S4 (c.) and (d.), the complex frequency map of log (|R(ω)|) and Arg (R(ω)) are shown for light impinging from the top medium. A reflection zero is seen at about ω/ω 0 ≈ 2.16 − i0.005 which is indeed the conjugated value of the reflection zero observed in Fig. 2 (b.). Since the zero remains in the lower half of the complex plane, the phase accumulation on the real axis is smaller than 2π as seen in Fig. S4 (b.).

SVII. METASURFACE SUSCEPTIBILITY MODELING
For known scattering parameters, relations (4) may be transformed to compute the metasurface effective susceptibilities as [69] χ e = − 2i (T + R − 1) k (T + R + 1) , For the general case of a metasurface with damped Lorentzian responses, the complex frequencies corresponding to the Kerker condition (i.e., full transmission through the metasurface, R = 0) are obtained by using (5) The analytical expression of the complex frequencies corresponding to transmission Zeros of a damped metasurface being very lengthy is omitted here.
In this situation, the complex frequencies of the transmission Zeros read where We now focus on the trajectories of transmission zeros for the case of undamped metasurface discussed previously in Fig. 3. For this purpose, we use (S10) to predict the branching behavior already observed above where transmission Zeros bifurcate from real to complex values. This happens when the term under the square root in the expression of B is negative, satisfied at the condition: (ω e − ω m ) 2 < A e A m 4c 2 . (S12) The exceptional points then occur at the boundary of this region where B becomes null: (ω e − ω m ) 2 = A e A m 4c 2 . (S13) The exceptional point positions are then provided by the following expression: ω = ± √ ω e ω m .
Looking at the asymptotic behaviors of Poles and transmission Zeros, we can note in Fig. 4(c) that the real frequencies of the Poles and Zeros tend toward each other far from the region of interest, i.e., h > 250 nm and h < 100 nm or when condition (S12) is not satisfied. This is verified by neglecting the contributions from A e A m in (S9) and (S10), which yields