Non-resonant exceptional points as enablers of noise-resilient sensors

Exceptional point degeneracies (EPDs) in the resonant spectrum of non-Hermitian systems have been recently employed for sensing due to the sublinear response of the resonance splitting when a perturbant interacts with the sensor. The sublinear response provides high sensitivity to small perturbations and a large dynamic range. However, the resonant-based EPD sensing abides to the resolution limit imposed by the resonant quality factors and by the signal-to-noise ratio reduction due to gain-elements. Moreover, it is susceptible to local mechanical disturbances and imperfections. Here, we propose a passive non-resonant (NR) EPD-sensor that is resilient to losses, local cavity variations, and noise. The NR-EPD describes the coalescence of Bloch eigenmodes associated with the spectrum of transfer matrices of periodic structures. This coalescence enables scattering cross-section cusps with a sublinear response to small detunings away from an NR-EPD. We show that these cusps can be utilized for enhanced noise-resilient sensing. Sensing schemes based on exceptional point degeneracies in the resonant spectrum of non-Hermitian systems have been recently criticized for being susceptible to various sources of noise and structural imperfections. Here, the authors study non-resonant exceptional point degeneracies in the spectrum of transfer matrices of periodic structures, based on which passive sensors with enhanced noise resilience are then proposed.

E xceptional point degeneracies (EPDs) are spectral singularities corresponding to points in the parameter space of a non-Hermitian operator at which its eigenvalues and the associated eigenvectors coalesce [1][2][3] . A prominent example includes EPDs in the resonant spectrum of non-Hermitian systems [3][4][5][6] . In their proximity, a small perturbation ε ( 1 leads to a sublinear response (SLR) in the resonant splitting Δω / ffiffi ε m p ) ε due to a fractional Puiseux expansion of the perturbed frequencies around an m th order EPD 1,2 . Such SLR provides an enhanced sensitivity to small perturbations ε [7][8][9] ; while also offering an additional advantage over other sensing schemes relying on high-Q resonances 10 : an enhanced dynamic range, which is the ability to measure both small and large perturbationrelated detunings. This observation has recently generated a substantial research effort in developing appropriate platforms where resonant EPDs are realized and their SLR is harvested for enhanced sensing application [11][12][13][14][15][16][17][18][19] . That said, the implementation of resonant-based EPD sensing has triggered an ongoing debate regarding the resolution limit and the signal-to-noise ratio (SNR) efficiency of such schemes [20][21][22] . Specifically, the EPD sensing schemes based on purely lossy systems have been hampered by the broadening of the resonance linewidths. The addition of gainelements can offset the losses, thereby, improving the resolution limit of the EPD sensing, however, they also introduce additional noise which becomes enhanced in the vicinity of the EPD and leads to a degradation of the SNR performance of the sensor. The noise can be intrinsic (e.g., due to amplification) or fundamental (due to the eigenbasis collapse at the EPD) and in some EPD-platforms, might offset the enhanced signal sensitivity, thus leading to an SNR which is not exceptional, but rather conventional [19][20][21] . Importantly, resonant EPD sensors, alike all resonant-based schemes, are susceptible to local mechanical disturbances (e.g., temperature variations, vibrations, etc.) and cavity imperfections that hamper their sensitivity.
Here, we propose a non-resonant sensing protocol based on EPDs occurring in the spectrum of operators other than the effective Hamiltonian of a resonant system. In our paradigm, the formation of an m th order non-resonant EPD (NR-EPD) occurs in the spectrum of transfer matrices of Hermitian periodic structures. Their existence enforces a stationary point in the Bloch dispersion relation ω k ð Þ $ ω SP þ ðk À k SP Þ m ; where m ≥ 2 Bloch modes coalesce and the group velocity v g ∂ω ∂k $ ðω À ω SP Þ mÀ1 m of a propagating wave inside the structure vanishes [23][24][25][26][27][28][29][30][31][32][33][34][35][36] . As a result, the differential scattering cross-section jσ T j 2 of a suitably designed incident wavefront shows a SLR with respect to small global parameter variations X SP ! X SP ± ν occurring in the proximity of the stationary point, i.e., where W T $ ν Àα is the energy density of the excited (slow) propagating mode. The exponent α dictates the formation of the cusp in Eq. (1) and its value is controlled by the incident wave via wavefront shaping techniques [37][38][39][40] . Here, we propose to utilize the sublinear response Eq. (1) of the differential cross-section, near stationary points, as a protocol for hypersensitive sublinear sensing. The special case of m ¼ 3 NR-EPDs, known as stationary inflection points (SIPs) 23,24,[27][28][29][30][31][32][33] , monopolizes our attention because of its resilience to common mechanical disturbances, structural imperfections 33,41 and losses 33,42 . Importantly, we show that the proposed SIP-sensing protocol exhibits an enhanced noise-resilient performance, as opposed to existing resonant based-schemes.
Results and discussion Implementation of SIP platforms. Unlike resonant EPDs, which rely on local perturbations done to the EPD-based platform, the proposed sensing platform (see "Methods") consists of two distinct units-the probing element and the SIP sensing element (see Fig. 1a). The probing element can be any type of optical or microwave high-Q frequency selective filter (see Fig. 1b, c) whose resonant frequency depends on the applied local perturbation (i.e., acceleration, rotation, particle, etc.). Such arrangement allows to make the probing element extremely compact-thus granting access to measurements on the microscopic scale in a potentially cramped environment, without affecting the measured physics. After being probed, the perturbation is then transduced to the sensing SIP element, which may be placed away from the probe. As a result, the SIP structure can be made sensitive and robust enough without major concern of its size, making such arrangement useful for microscale sensing applications. A proposed experimental implementation of stationary-point-based, on-chip sensors (e.g., accelerometers, gyroscopes, inclinometers), using an optics/microwave framework, is shown in Fig. 1 (see "Methods" for elaboration). Importantly, the proposed platform is scalable and, hence, applicable for various wavelengths ranging from optical to millimeter-wave and radio frequency. The SIP protocol can be utilized in a variety of applications ranging from avionics and temperature variation sensing to bio-and chemical sensing.
Transfer matrices and Bloch dispersion relation. The description of a wave propagating in an M-channel periodic structure with periodicity L 0 is typically done using the 2M-dimensional unit cell transfer matrix M ω ð ÞMðz ¼ z 0 þ L 0 ; z 0 ¼ 0; ωÞ. This non-normal operator connects the wave amplitudes Φ of a monochromatic wave at two different spatial positions of the structure at z ¼ z 0 þ L 0 and z 0 , through the relation where the real part of the Floquet-Bloch wavenumber ReðkÞ 2 ½Àπ=L 0 ; Àπ=L 0 is defined up to a multiple of 2π=L 0 (first Brillouin zone). At the right-hand-side of Eq. (2) we have imposed the Bloch theorem that requires Φ k z ¼ z 0 þ L 0 À Á ¼ e ikL 0 Φ k ðz 0 Þ, reflecting the periodicity of the underlying structure. Finally, the Bloch dispersion relation ω ¼ ωðkÞ is evaluated by calculating the 2M eigenvalues λ ω ð Þ of the transfer matrix by solving the following secular equation:  23,43 . Below, we will refer to the propagating Bloch modes with positive (negative) group velocity v g > 0 (v g < 0) and to the evanescent Bloch modes with Im k ð Þ > 0 (I m k ð Þ < 0) as forward (backward) waves.
EPDs in the spectrum of transfer matrices and their relation to stationary points in the Bloch dispersion relation. The nonnormal nature of the transfer matrix M allows for the formation of NR-EPDs in systems with non-trivial topology. Such degenerate points occur in the spectrum of MðXÞ by appropriately tuning one or a few of its parameters X (e.g., frequency ω; magnetic field, etc.). In this case, a complete basis is formed by augmenting the Bloch eigenvectors with generalized eigenvectors Φ r q . The latter are found by implementing the standard Jordan chain procedure 1,2,23,33 , defined by the set of vectors that satisfy the recursive equations: where m (2 ≤ m < 2M) is the order of the NR-EPD (NR-EPDÀm), eigenvector while the remaining m À 1 eigenvectors Φ r q > 1 are the generalized eigenvectors. One can show that the generalized eigenvectors diverge along the propagation direction as Φ q z ð Þ / z qÀ1 e ikz Φ q ð0Þ: At an NR-EPDÀm, the transfer matrix M X NR-EPD À Á is similar to a Jordan form or a matrix containing Jordan blocks, i.e., M X NR-EPD À Á ¼ VΛV À1 , and therefore, it is not diagonalizable 1,2,23,33 . The similarity transformation matrix V has columns consisting of the regular and the generalized eigenvectors. In the case when M has only one occurrence of an EPD in its spectrum, the matrix Λ consists of a Jordan block of size m and a diagonal matrix of size 2M À m with diagonal elements being the eigenvalues λ n with geometric multiplicity 1. Generally, each Jordan block of size m will include one propagating Bloch mode with zero group velocity and m À 1 generalized eigenvectors with algebraically diverging amplitude with respect to the propagation distance z.
An important feature of an NR-EPD is that any perturbation X NR-EPD ! X NR-EPD þ ν that detunes the system away from the degenerate point results in a fractional power series (Puiseux series) of the eigenvalues with respect to the perturbation parameter ν. In other words, when an NR-EPDÀm transfer matrix M X NR-EPD À Á ; which is similar to a matrix containing at least a Jordan block of size m m, is perturbed as where Δ is a constant perturbation matrix, then the degenerate eigenvalues λ, will obey the fractional expansion The perturbed eigenvalue λ q in Eq. (5) can be also written as c n δk n which allows us to deduce that the perturbed Bloch wavevector (to the first order correction in detuning ν) is δk $ ffiffiffi ν m p . When the associated detuning ν from X NR-EPD is identified with the frequency variation from the NR-EPD frequency ω NR-EPD ¼ ω SP , we get: Equation (6) signifies the formation of a stationary point in the Bloch dispersion relation and connects an NR-EPDÀm (and the associated size of the Jordan block) with the order of the stationary point. From Eq. (6) we deduce the presence of a slowlight mode with a vanishing group velocity v g The most familiar example of a stationary point is the regular band-edge (RBE) corresponding to m ¼ 2. Higher order Fig. 1 Implementation of an SIP-based sensing protocol. a The schematic illustration of the sensing platform. The protocol involves a source which emits a broadband signal to the sensing platform which acts as a high-Q frequency selective filter, with frequency that depends on the applied perturbation. b An optomechanical Fabry-Perot cavity whose one wall is acting as a test-mass. c A whispering-gallery-mode resonator whose resonant frequency is changed due to the Sagnac effect. Different platforms supporting a stationary inflection point (SIP) singularity: d a multilayer structure 25,26,32 , e a coupled microstrip-waveguide [27][28][29]35 , where H indicates an external magnetic field strength, or f a serpentine waveguide 36 . stationary points with m > 2 need to be specifically engineered and are divided into two categories: The stationary points corresponding to even order NR-EPDÀm 0 s that occur at the band-edges, and the odd order NR-EPDÀm 0 s that occur inside the band and form an inflection point in the Bloch dispersion relation. The former are known as degenerate band-edges (DBEs) while the latter as stationary inflection points 23,24 . As opposed to the RBE, the higher-order stationary points include the presence of degenerate evanescent modes in addition to the propagating modes 23,24,31 . These evanescent modes contribute significantly in the formation of the cusp anomaly in the differential scattering cross-sections. Their engineered excitation allows us to control the type of divergence that the energy density of the excited (slow) propagating mode W T demonstrates, see Eq. (1).
Although these cusp anomalies can occur for both odd and even order-m NR-EPDs, our focus will be on SIPs. Among all odd m NR-EPDs, the case of m ¼ 3 will be monopolizing our attention. Various reasons led us to this decision: First, an SIP can be engineered in a way that shows a symmetric cusp anomaly with respect to a left/right detuning from the NR-EPD conditions when it is designed to occur in the middle of the band. Second, RBEs and DBEs are often overwhelmed with much more powerful, giant slow wave resonances when the system is turned to a scattering set-up 24 . Such high-Q resonances destroy the formation of cusps and result in sensing protocols with a small dynamic range 10 . Third, SIPs are not particularly sensitive to the size and shape of the underlying photonic structure and, when compared to Fabry-Perot or transmission band-edge resonances (where the whole photonic structure acts as a resonator), they are more resilient to absorption and structural imperfections 33,41,42 . Finally, the SIP of order m ¼ 3 has been already implemented experimentally using various photonic platforms (see Fig. 1d-f), for efficient slow-light conversion [27][28][29] . While in these latter studies, the SIP of order m ¼ 3 was promoted due to its high conversion prospects 23-36 , our sensing scheme takes advantage of the opposite scenario i.e., the possibility of total decoupling between the incident light and the slow mode.
Cusp anomalies in the differential reflectance near a stationary point. Next, we analyze the differential reflectance ΔRðνÞ when a control parameter X (e.g., the frequency ω of the incident wavefront) is detuned away from the NR-EPD conditions by ν ¼ X À X NR-EPD . For the sake of the argument, we assume only one Jordan block of size m in the similarity matrix Λ. Furthermore, we assume the simplest possible scattering scenario of a lossless semi-infinite SIP structure of order m ¼ 3 occupying the positive semi-infinite space z ≥ 0. The validity of our conclusions in the case of finite structures (in the presence of losses), has been tested in the next section.
Our presentation below follows closely the study of Figotin and Vitebskiy 24 . At the interface, the incident Ψ I ω; z ð Þ; reflected Ψ R ω; z ð Þ and transmitted Ψ T ω; z ð Þ waves, must satisfy the boundary condition where the reflection coefficient has been absorbed in Ψ R . When the system is detuned away from the NR-EPD (i.e., ν ≠ 0), the transmitted wave Ψ T ω; z ð Þ inside the slow-light structure can be decomposed into a superposition of the M forward Bloch modes Φ þ n (see "Methods"). Specifically, Furthermore, the energy conservation condition requires that the transmitted and reflected waves satisfy the relation where S I , S T X ð Þ, and S R X ð Þ are the energy fluxes of the incident, transmitted and reflected waves. We further assume that the incident wave is normalized as S I ¼ 1. Since evanescent modes do not contribute to the energy flux, S T is where v ðnÞ g indicates the group velocities of the forward propagating Bloch modes inside the structure. In the case that the incident monochromatic wavefront that does not excite any fast-propagating modes inside the medium, Eq. (10) simplifies to is the energy density carried by the slow mode and v slow g is the corresponding group velocity given by Eq. (7). Therefore, the analysis of the exact form of the S T singularity collapses to the investigation of the divergence of the energy density W T $ ν Àα with respect to the detuning from the NR-EPD.
We first recall that for ν ≠ 0; in the specific case of an SIP of order m ¼ 3, the excitation Ψ T z ð Þ can be written as a superposition of a forward slow-propagating Bloch mode Φ þ slow ω; z ð Þ and of a forward evanescent mode Of course, the decaying evanescent contribution to the transmitted wave becomes negligible at a certain distance z c / 1=I m k ð Þ from the interface and this explains the approximation i.e., Nevertheless, the excitation of an evanescent mode is detrimental in the formation of a cusp as it can lead to a finite energy . Such a field intensity divergence is compatible with the boundary condition Eq. (8) (which dictates that Ψ T z ð Þ at the boundary z ¼ 0 is finite), provided that the incident wave excites both the slow forward propagating and the forward evanescent modes in a way that they are interfering destructively at the interface i.e., The effect of the dramatic amplitude growth of the transmitted wave in the vicinity of an SIP is a feature of the frozen mode regime. Past efforts [23][24][25][26][27][28][29][30][31] utilized this generic property of the frozen mode regime for enhancing light-matter interactions. Here, however, our interest lies in the opposite scenario where S T / ν j j 2=3 and, therefore, the transmittance (or the differential reflectance) forms a cusp. From the above discussion, it becomes clear that the formation of the cusp Eq. (1) requires the design of an incident wavefront which does not excite a Bloch evanescent mode. In this case, the differential reflectance in the proximity of an SIP of order m ¼ 3 behaves as where we have assumed that in the lossless scenario at the SIP the , and T ω ð Þ S T ðωÞ=S I and R ω ð Þ ¼ S R ðωÞ=S I are the transmittance and reflectance of the incident wavefront. Equation (11) signifies a sublinear response of the differential reflectance to frequency detuning and can be used as a measurand for enhanced SLR sensing protocols that also demonstrate a large dynamic range. In fact, the above SLR is still applicable for any global parameter variation ν ¼ X À X NRÀEPD which preserves the form Eq. (7) of the slow-light group velocity. Below, we will demonstrate that one such case is associated with variations of the external magnetic field applied to a slow-light structure. Let us finally remark, that a similar argument that led to Eq. (11) for the case of a stationary point of order m ¼ 3, can also apply for the stationary point of order m ¼ 2 (RBE) leading to a square root cusp i.e., ΔR ν ð Þ $ ffiffiffi ν p : Wave simulations and Coupled Mode Theory (CMT) modeling. We have confirmed our proposed sensing protocol by performing full-wave simulations with a realistic photonic platform that has been experimentally shown to demonstrate an SIP of order m ¼ 3 [27][28][29] . We consider a semi-infinite periodic system whose unit cell consists of a pair of coupled perfectly conductive microstrip waveguides on a ferrite magnetic substrate (see Fig. 1e and "Methods" for the design details). The unit cell contains one straight waveguide and one meandered waveguide. The waveguides are coupled together in the spatial domain, where they are closely separated from one another. We assumed that there is an applied out-of-plane constant magnetic field H % 86 mT, resulting in a violation of time-reversal symmetry, necessary for achieving an NR-EPD of order m ¼ 3 in the Bloch modes of the 4 4 transfer matrix M which describes the unit cell of the system. Using COMSOL's finite element method (FEM) solver, we have calculated the S-parameters of the unit cell by probing the system with 50Ω impedance ports and retrieved the associated transfer matrix M. This allowed us to calculate the dispersion relation and the Bloch modes, which are required for the analysis of the semi-infinite structure (see "Methods" for details).
In Fig. 2a, we show the calculated dispersion relation which displays an SIP of order m ¼ 3 (SIP-3) at the frequency f SP ¼ 2:00645 GHz (see the dashed horizontal line). In the inset of Fig. 2b, we also show the reflectance R as a function of frequency f . The SIP-3 frequency f SP , where the reflectance develops a cusp, is shown in the inset by the vertical dotted line. In the main panel of the same figure, we report the differential reflectance ΔR as a function of frequency detuning ν from the SIP frequency f SP . The incident wavefront is designed following the criteria specified in the previous section to guarantee the SLR of Eq. (11). The best fit (see black dotted line in Fig. 2b) indicates that the differential reflectance varies with the frequency detuning ν as ΔR $ ν j j 0:7 , which is consistent with the theoretical expectations of Eq. (11). We have also checked the validity of the SLR of ΔR with respect to other (global) parameter detunings. In Fig. 2c, we report the response of the differential reflectance ΔR with respect to magnetic field variations ν ΔH from the value of the applied magnetic field strength H SP at the SIP. The FEM simulations reveal a scaling ΔR $ ΔH j j 0:66 in perfect agreement with Eq. (11). The scattering properties of the simulated photonic circuit of Fig. 1e in the proximity of the SIP can be analyzed using CMT modeling (see "Methods"). Due to its generality, CMT modeling allows us to extend our conclusions beyond the specific platform of Fig. 1e, to any system that demonstrates stationary points in its Bloch dispersion relation. The geometry of the model is shown in Fig. 3b, while its dispersion relation is reported in the inset of Fig. 3a.
In the main panel of Fig. 3a, we report the differential reflectance ΔR ν ð Þ versus the detuning ν from the SIP frequency. The analysis required the evaluation of the eigenmodes of the transfer matrix of the slow-light structure and the decomposition of the incident wavefront Ψ I in this basis, see Eq. (9). Furthermore, a small imaginary part ðI m ε 0;1 $ 10 À5 Þ has been introduced (see "Methods") to simulate natural losses occurring in the structure and to allow us to lift the NR-EPD degeneracy in a controllable manner, in order to implement the decomposition process Eq. (9) numerically. We have made sure that such a prepared wavefront satisfies the boundary condition Eq. (8) together with the requirements for the appearance of the cusp anomaly in the reflectance, as outlined in the previous section. The data confirms nicely the validity of Eq. (11) for at least three orders of magnitude. The same analysis has been performed using, as a detuning parameter, the variations of the Peirels' phase from its NR-EPD value ϕ SP and for fixed incident frequency ω ¼ ω SP : These results are also presented in the main panel of Fig. 3b. Furthermore, we have tested that these results are robust for large, but finite, slow-light samples when a small amount of losses are present in the system. These losses are required to suppress Fabry-Perot resonances in the reflectance spectrum which can mask the otherwise robust SLR scaling. Achieving successful scaling in finite samples further indicates the viability of the proposed NR-EPD sensing platform. However, the optimal local loss strength γ with respect to the size of the system L is expected to be model dependent and determined by the condition that the absorption length ξ γ in the vicinity of the stationary point is less than the system size. If the loss strength is too high, the impedance mismatch will inevitably deteriorate the desired SLR response. In Fig. 3c, we also show how the finite size of the structure affects the scaling of the differential reflectance ΔR. From the figure, it is seen that the fractional response is already evident when the number of unit cells is twenty. Larger system sizes result in an increased dynamic range due to reduction of the lower bound of the sublinear scaling.
Noise analysis. While enhanced sensitivity to small perturbations is an important metric, the precision of the measurement is another important characteristic of the efficiency of a sensor. This is defined as the smallest measurable change of the input quantity given by the noise of the sensor output. Noise can stem from a variety of sources, including mechanical vibrations and mesoscopic fluctuations due to environmental thermal fluctuations, signal noise generated by the coupling to the interrogating channels, quantum uncertainty or fundamental detector resolution limits, and it can never be fully eliminated. It is therefore vital to study both-sensitivity and noise-in tandem. Such analysis allows one to estimate how noise in the observable (e.g., ΔR) is translated to uncertainty in the measurand (e.g., ν), which defines the actual precision of a sensor. The precision of our stationary point sensing protocol due to noise is scrutinized by the computational simplicity that the CMT modeling offers and is reconfirmed by COMSOL simulations for the platform shown in Fig. 1e. Below, we distinguish between two types of noise that might affect the performance of a sensor 9 : (a) multiplicative noise associated with classical noise sources describing a noisy NR-EPD, and (b) additive noise associated with noisy input channels. For a better assessment of the NR-EPD sensing efficiency we also compare the precision provided by SIP-3 and RBE sensing protocols. A quantitative description for the precision of the stationarypoint-based sensors is provided by analyzing the detuning error σðνÞ defined as where is the standard deviation of ΔR ν ð Þ due to noise fluctuations and χ is the sensitivity of the sensor, where Á h i indicates a noise averaging. The detuning error Eq. (12) provides an estimation of the uncertainty in the evaluation of the detuning ν (which is the signature that the perturbing agent leaves when it interacts with the sensing platform) via the output measurement associated with the differential reflectance.
Multiplicative noise. We first analyze the influence of practically inevitable fabrication imperfections or slow fluctuations of the environment associated with, temperature or pressure variations that affect the constituent parameters of the materials, and other types of mesoscopic fluctuations in the system parameters. Such fluctuations are characterized by a very large correlation time implying (quasi-)static disorder and lead to the existence of additional parasitic degrees of freedom of the system 44 . Therefore, it is instructive to study how these parametric fluctuations are translated to the error in the measured detuning. This situation is modeled by weakly fluctuating frequencies ε 0 þ δε 0 and ε 1 þ δε 1 of each resonant mode of the CMT model (see Eq. (22) in "Methods"). For demonstration purposes, we have assumed that both δε 0 and δε 1 are independent random variables taken from a uniform distribution ÀW; W ½ . We first consider an SIP-3 scenario in a finite-size, disordered system. To mimic the behavior of a semi-infinite structure, we have introduced losses, in such a way that the associated absorption length is smaller than the size of the system (see Methods for details). In Fig. 4a, we report the CMT results for the differential reflectance ΔR ν ð Þ (gray circles) versus the frequency detuning ν from ω SP . The ensembled average differential reflectance ΔR ν ð Þ is also shown with a blue solid line. We find that the power law scaling ΔR ν ð Þ $ ν α with the best-fitting value of α % 2=3 persists for three orders of magnitude in detuning, despite the presence of the disorder. Nevertheless, a smearing of the SIP cusp is unavoidable leading to the formation of a plateau ΔR ν ! 0 ð Þ $ W β for very small detunings ν c $ W 3β=2 (see inset of Fig. 4a). These detunings define the sensitivity bound of our SIP-sensors as far as the mesoscopic fluctuations are concerned. The numerical analysis gives the best-fit of the exponent, which is β % 2=3 (see "Methods" for further analysis).
Further statistical processing of ΔR ν ð Þ allows us to evaluate σ ΔR ν ð Þ and χ and from there, via Eq. (12), the detuning error σ ν . To get a better appreciation of the robustness of the SIP-3 based sensing with respect to disorder, we have also calculated σ ν associated with an RBE. The comparison is shown in Fig. 4b for three values of the disorder strength W. In this figure, we have highlighted (red region) the domain where the error in the measured detuning is larger than the detuning itself and, therefore, the precision of the sensor has completely deteriorated. In all cases, σ ν is decreasing as we are approaching the corresponding resolution limit (black dashed line). However, for the same disorder strength W, the detuning error of the SIP-3 based sensing is smaller by an order of magnitude than the detuning error of RBE-based sensing protocol indicating its superiority as far as (long-correlation time) multiplicative noise is concerned.
Additive noise. Next, we analyze the effect of noise due to lowfrequency thermal fluctuations in the input channels. For this Fig. 3 Differential reflectance near SIP frequency obtained from two-chain CMT model. a Logarithmic plot of the differential reflectance ΔR (blue line) as a function of the frequency detuning ν from the stationary inflection point (SIP) frequency ω SP . The black dashed line indicates a variation of the differential reflectance which is ΔR / ν j j 0:66 . The inset shows the dispersion relation of the Coupled Mode Theory (CMT) model. The SIP frequency ω SP is indicated by a horizontal dashed line. b Logarithmic plot of the differential reflectance ΔR (blue line) as a function of the Peirels' phase detuning ν ¼ Δϕ from the critical phase ϕ SP . The black dashed line indicates a variation of the differential reflectance which is ΔR / ν j j 0:66 . In the background we show a schematic of the CMT model Eq. (12). c Logarithmic plot of the differential reflectance ΔR (color lines) as a function of the frequency detuning ν from the SIP frequency ω SP in case of a finitesize system for various number of unit-cells L. The black dashed line indicates a variation of the differential reflectance which is ΔR / ν j j 0:66 . In order to optimize the scaling performance, we have used a ramped loss γ n ð Þ ¼ γ max purpose, we have performed Monte-Carlo simulations in a semiinfinite structure associated with the CMT Hamiltonian (see Eq. (22) in "Methods"). From these simulations we have extracted the standard deviation σ ΔR of ΔR ν ð Þ due to the presence of the additive noise and the sensitivity χ of the sensor, which allowed us to evaluate the detuning error σ ν via Eq. (12).
In Fig. 5a, we report some typical results of the SIP-based sensing scheme, accounting for the influence of input signal noise on the measured differential reflectance ΔR ν ð Þ. The blue solid line indicates the mean value of ΔR ν ð Þ for each specific value of frequency detuning ν. The height of the gray domain surrounding the blue line represents the standard deviation σ ΔR of ΔR ν ð Þ at each value of ν. When comparing the evaluated σ ΔR for two distant values of ν (see black vertical double-sided arrows) it is found that σ ΔR does not experience noticeable variations as a function of ν and remains approximately constant. At the same time, Eq. (12) implies that for σ ΔR ν ð Þ % const:, the detuning error will scale inversely proportional to the sensitivity i.e., σ ν $ χ À1 ðνÞ. Therefore, in the limit of small detuning values where the sensitivity is higher, we expect a decreasing detuning error σ ν . Figure 5b shows the probability density PðνÞ of the simulated detuning measurements performed for two distant values of ν (blue lines). The corresponding values of standard deviation σ ν are indicated by double-sided blue arrows. As expected, when the simulated measurements are performed close to the stationary point  singularity, the standard deviation σ ν decreases, indicating an enhanced precision of the sensor. A panorama of σ ν versus the detuning ν, for an SIP-based (blue line) and an RBE-based (red line) sensing protocol is shown in Fig. 5c. In both cases, we have found that σ ν $ 1=χ (see the blue and red dashed lines corresponding to σ ν $ ν 1=3 and σ ν $ ν 1=2 , respectively) which is a consequence of the fact that the standard deviation σ ΔR remains approximately constant with respect to the detuning from the stationary point frequency ω SP : The smaller detuning error of the RBE sensor compared to the SIP sensor in case of small detuning, does not imply that the RBE sensor is superior to the SIP sensor, since the former is extremely fragile to disorder and losses. In the same subfigure, we report COMSOL results for the case of the coupled microstrip waveguide of Fig. 1e (light blue line). It is instructive to point out that the results of the noise analysis of the Monte-Carlo simulations for the RBE sensing protocol agree with the recent experimental findings of enhanced precision accelerometers operating in the vicinity of a Wigner Cusp Anomaly (WCA) 45 . The WCAs are square root cusps in the differential scattering cross-section of processes around the frequency threshold of a newly open channel, and they can be seen as a generalization of RBEs.

Conclusions
We have identified a paradigm of enhanced sensing protocols that utilize non-resonant exceptional point degeneracies occurring in the spectrum of transfer matrices of periodic structures. Their presence enforces a stationary point in the Bloch dispersion relation, which leads to a cusp in the reflectance from these structures. We have shown that the reflectance variation with respect to the frequency detuning (or any other global parameter detuning) from the NR-EPD configuration is sublinear and can be used as a measurand for enhanced sensing with a large dynamic range. In this respect, our proposed SLR sensing protocol is complimentary to the resonant-based EPD schemes that detect local perturbations imposed on the system from a perturbing agent. Furthermore, the enhanced sublinear sensitivity is resilient to noise generated by the input channels as opposed to resonant EPDs 5,7-9,11-22 where (fundamental) noise in the proximity of an EPD, due to the resonant eigenbasis collapse, is enhanced and, in some cases (depending on the underlying platform), might offset the enhanced sensitivity [19][20][21] . Sensors that utilize stationary inflection point NR-EPDs show additional robustness to mesoscopic fluctuations due to environmental temperature variations and fabrication imperfection. Our sensing proposal does not involve active elements, and therefore, it does not suffer from excess noise effects. Although the analysis performed here has been confined to the first two stationary points of order m ¼ 2; 3, it can be extended to higher-order NR-EPDs where one can achieve an even higher sensitivity depending on the number of evanescent modes excited by the incident wavefront. It will be interesting to extend these studies in this direction and analyze the robustness of these schemes to noise.

Methods
Proposed experimental implementation of stationary-point-based avionic sensors. In case of avionic sensing (e.g., acceleration), the probing element can be an optomechanical Fabry-Perot cavity whose one wall is acting as a test-mass (see Fig. 1b). The same concept can be used in microwaves where now the variable cavity consists of an LC circuit with a variable capacitor plate playing the role of the test-mass. When an acceleration is applied, the test-mass is displaced due to the inertia; thus, inducing a variation in the length of the cavity and, consequently, in its resonant frequency. Another possible implementation of the sensing protocol is as a hypersensitive gyroscope. In the optical framework, the probing platform will be a whispering-gallery-mode whose resonant frequency is changed due to the Sagnac effect, which is a shift of the resonant mode frequency by an amount proportional to the angular velocity of the rotating platform (see Fig. 1c). A microwave implementation of the whispering-gallery-mode sensing platform might involve a combination of two microwave accelerometers (see Fig. 1b) which are normal to one another. In both cases (see Fig. 1b, c), the variable resonators act as filtering devices that turn a broadband incident signal into a continuous-wavelike signal with a narrow spectral width around a detuned resonant frequency induced by the motion of the sensing platform. The continuous-wave-signal transmits through the splitter and is reflected by the sensing element which supports an NR-EPD of order m ¼ 3 corresponding to a stationary (inflection) point in its Bloch dispersion relation. The reflected signal is measured by a detector and shows a cusp anomaly in the proximity of the stationary point frequency f SP . The reflectance variations ΔR ¼ R À R SP with respect to the frequency detuning ν ¼ f À f SP demonstrate a sublinear response ΔR $ ν 2=3 which is used as a measurand for the proposed sensing, see Eq. (1). Our stationary point sensing scheme utilizes intensity variation measurements (i.e., transmittance/reflectance) near non-resonant EPDs as opposed to resonant shift measurements. The latter might be masked by a linewidth broadening of the resonances that appear in the transmission (or the reflection) spectrum or by the generation of additional noise due to the presence of gain elements. Furthermore, resonant shift sensing requires broadband frequency sweeps, which impose an upper bound on the dynamic range. Such measurements are done in two ways: (i) signal (intensity) measurements, by probing the system at individual frequencies within the scanned domain, (ii) Fourier transform of the temporal signal emitted by the system. In the former case one needs to perform a large number of intensity measurements for each respective frequency, which results either in longer sampling time compared to the direct intensity measurements or worse signal-tonoise ratio. The latter measuring scheme is applicable only to active systems operating at or above the lasing threshold, which, besides adding extra complexity to the platform, also results in the addition of the quantum noise generated by the system.
Decomposition into Bloch modes and wavefront shaping protocol. The analysis of the transport characteristics of a semi-infinite structure is typically based on a scattering matrix approach. The latter can be directly extracted from the analysis of the following three 2M 2M transfer matrices: (a) the matrix M S that determines the wave propagation from a unit cell to a unit cell inside the semi-infinite structure; (b) the transfer matrix M L that dictates the wave evolution in free space (or the leads), and; (c) the transfer matrix M I which describes the boundary conditions Eq. (8) that must be satisfied by the field at the interface z ¼ 0. The latter can be expressed in a compact form as: where Ψ j ðj ¼ 1; 2Þ denote the fields before and after the interface respectively. While the scattering matrix connects incoming to outgoing waves and, therefore, provides easy access to the transmission/reflection coefficients associated with the scattering process, the transfer matrix formalism allows us to identify the associated Bloch modes of the semi-infinite structure and, via this identification, establish an appropriate basis where the conditions which lead to Eq. (1), for the formation of a cusp in the differential reflectance, can be formulated.
We start with the description of the process that allows us to decompose the fields into the Bloch eigenmodes (see Eq. (9)). First, we diagonalize the transfer matrices M S , M L in order to extract the Bloch modes and classify them as forward/backward and propagating/evanescent according to their wavenumbers and group velocities as described in the main text. This mode-basis is then arranged as columns of 2M m x j matrices B x j , where m x j denotes the number of each type x of the Bloch eigenmodes at a respective frequency. For example, x ¼ f denotes forward propagating modes, x ¼ b backward propagating modes, and x ¼ e evanescent modes which decay away from the interface. Bloch modes that grow away from the interface are excluded from the decomposition based on the physical requirement that the fields must not diverge as jzj ! 1. The decomposition proceeds as, where α x j contain the expansion coefficients in the basis of the Bloch eigenmodes. Inserting Eq. (14) in Eq. (13) leads to the following relation which, after appropriate re-arrangement, can be written as Equation (16) relates the incoming propagating and the outgoing propagating and evanescent modes on the right and left sides, respectively.
Next, we want to rewrite the above relation in terms of the incoming ϕ þ (defined in an (m To this end, we first introduce the projection operators ε f 1;2 , ε e 1;2 ; and ε b 1;2 defined as ε Substituting in Eq. (16) the expressions for α f ;b;e 1;2 from Eq. (17), we find the relation that connects the incoming with the outgoing (including the evanescent) modes, Equation (18) defines the non-unitary scattering matrix e S in terms of the Knowledge of e S allows us to extract further information about the expansion coefficients α x j and, consequently, identify conditions under which an incident wavefront can (cannot) excite specific Bloch modes inside the semi-infinite stationary point structure.
Let us, for example, analyze the expansion coefficients of a generic incident wavefront that is injected into a semi-infinite structure that supports a stationary point from two single-mode transmission lines. We assume that the wavefront can be written as a linear combination of propagating waves ϕ ðnÞ þ ¼ ðδ 1;n ; δ 2;n Þ T (n ¼ 1; 2) in each of the transmission lines i.e., ϕ þ ¼ ϕ ð1Þ þ þ βϕ ð2Þ þ . The complex amplitude β contains information about the relative magnitudes and phases of the two propagating waves. From superposition principle, we know that the resulting Bloch modal excitation in the semi-infinite structure can be written as a linear combination of the modes that are excited by each individual incident propagating wave ϕ ðnÞ þ : Since there is only one forward evanescent mode in the SIP structure, the amplitude of the expansion coefficients associated with the evanescent mode α e 2 ðnÞ are evaluated using Eq. (18). We get α e 2 n ð Þ ε e 2 À Á y e Sϕ n ð Þ þ : The requirement of zero evanescent mode excitation now translates to the condition α e 2 ð1Þ þ βα e 2 ð2Þ ¼ 0 which allows us to extract the wavefront shaping parameter β: The exact value of β might depend on the frequency of the incident wavefront or other parameters X of the semi-infinite structure. However, our detailed numerical analysis indicated that the extracted value of β at X ¼ X SP maintains the sublinear scaling of the differential reflectance.
As a point of caution, Eq. (18) does not assume a flux normalization of the incident wavefront. The latter is crucial when one needs to evaluate transmission coefficients that describe scattering processes between channels that support different group velocities. In our case, we focus on reflectances/transmittances to channels with the same group velocity defined by the dispersion relation of the single mode transmission lines of equal impedance.
Geometry of the unit cell and COMSOL simulations. The full-wave simulations have been performed with a periodic structure which has a unit cell that is shown schematically in Fig. 6. The unit cell consists of a ferrite board (green) with a perfectly conductive grounded bottom surface. The board hosts a pair of perfectly conductive microstrip waveguides indicated with a brown color. The unit cell consists of two spatial domains of length l 1 ¼ 4:318 mm, and l 2 ¼ 5:588 mm, where the two waveguides are separated by distances s 1 ¼ 3:5052 mm, and s 2 ¼ 0:0889 mm respectively. Within the first spatial domain, the width of the bent waveguide is W 1 ¼ 1:524 mm, while the straight one has a constant width W 2 ¼ 0:762 mm within the whole unit cell. Within the second spatial domain, the width of the bent waveguide is reduced to W 3 ¼ 1:016 mm. Finally, the board is made from a magnetic garnet G-810 which has thickness T ¼ 1:524 mm. The dielectric permittivity of the board is ε R ¼ 14:6, while the magnetization is 4πM s ¼ 800 G.
Coupled Mode Theory modeling. The lower chain consists of coupled modes with resonant frequencies ε 0 while the n.n. coupling between them is V 0 2 R. The upper chain consists of coupled modes with resonant frequencies ε 1 ¼ 0 and n.n. coupling between them, which is V 1 2R. The same n.n. coupling constant describes the coupling between the modes of the first chain with the modes of the second chain. Finally, the diagonal coupling between modes of the first and second chain is V 2 e iϕ ðV 2 2 RÞ. The Peirels' phase ϕ models a magnetic flux which is responsible for violating the time-reversal symmetry of the system. The effective CMT Hamiltonian of the model has a block-tri-diagonal form, H nl ¼ δ n;l H n;n þ δ nÀ1;l H y n;nþ1 þ δ nþ1;l H n;nþ1 ; ð22Þ where the 2 2 block matrices have elements, H m;m 0 n;n ¼ δ m;2 δ m;m 0 ε 0 þ δ m;m 0 ± 1 V 1 and H m;m 0 The corresponding eigenvalue problem can be written in the following form, where ω is the eigenfrequency and hm ψ n represents the complex amplitude of the field occupying the site in the m th array of the n th unit cell. Substituting where ϵ 1 ðkÞ ¼ 2V 1 cosðkÞ; ϵ 2 ðkÞ ¼ ϵ 0 þ 2V 0 cosðKÞ and uðkÞ ¼ V 1 þ V 2 e iðkþϕÞ . The dispersion relation ωðkÞ is obtained by solving the associated secular equation det H k ð Þ À ω k ð Þ Á1 2 À Á ¼ 0. We get the following two-band dispersion curve which for an appropriate choice of the model parameters ε 0 ; V 0 ; V 1 ; V 2 ; ϕ À Á SP % ð4; 1; 2; 3; 3:32Þ demonstrates an SIP-3 at ω SP % 6:328 (k SP % À0:507), see the inset of Fig. 3a.
The system turns to a scattering setup by coupling each of the left-most resonant modes of the n ¼ 1 unit cell of the sample with transmission lines that consist of coupled resonant modes of resonant frequency ε L ¼ 0 and n.n. coupling V L ¼ À3:5. Each of the transmission lines supports propagating waves with dispersion relation ω k ð Þ ¼ ε L þ 2V L cosðkÞ.
Transfer matrix for the CMT model. The calculations reported in Fig. 3 were performed using a CMT formalism where we considered a semi-infinite SIP structure with two single-mode transmission lines attached to each one of the resonant modes of the first unit cell. Our analysis follows closely the steps presented in the second subsection of the "Methods". The Bloch modes of the leads have a trivial plane waveform with dispersion relation ω ¼ ε L þ 2V L cos k a ð Þ, where a ¼ 1 and ε L ; V L À Á ¼ 0; À3:5 ð Þ . The Bloch modes of the semi-infinite structure Eq. (23) are extracted from the transfer matrix of the unit cell, which takes the form, where1 M is the M M identity matrix,0 M is the M M matrix of zeroes, jΨ n i ¼ ðjψ nþ1 i; jψ n iÞ T , ξ 1;n ¼ H À1 n;nþ1 ðω1 M À H n;n Þ and ξ 2;n ¼ ÀH À1 n;nþ1 H y n;nþ1 : Finally, the transfer matrix of the interface is given by Eq. (26) by substituting the 2 2 submatrix H y n;nþ1 with the hopping matrix V y that connects the transmission lines with the first unit cell and has matrix elements V Coalescence parameter of the transfer matrix eigenvectors and EPD. A detailed analysis of the eigenmodes of the transfer matrices associated with the CMT model and of the microstrip waveguide of Figs. 1e and 6 allows us to establish convincing evidence that the SIP occurring in the Bloch dispersion relation is associated with the formation of the EPD in the spectrum of the corresponding transfer matrices. We start our analysis with the evaluation of the spectrum of the transfer matrix (see Eq. 26) associated with the CMT model. The four eigenvalues of the unit-cell transfer matrix are expressed in terms of the wavevectors k q as λ q ω ð Þ e ik q L 0 (see Eqs. 2, 3) and can be calculated via direct diagonalization of M n . Two of the four eigenvalues correspond to forward and backward propagating modes and have real k values (red and green solid lines in Fig. 7a, b respectively). When ω is plotted as a function of k, it provides the Bloch dispersion relation of the system ωðkÞ, which demonstrates an SIP (see Fig. 7a). The remaining two eigenvalues are associated with one forward and one backward evanescent mode, and they are characterized by a complex k vector (blue and orange solid lines in Fig. 7a, b respectively). From Fig. 7a, b, we see a third order EPD associated with the coalescence of the slow forward propagating mode (red line) and the forward/backward evanescent modes (blue/orange lines).
In Fig. 7c, we further report the eigenvector coalescence parameter defined as 46 , where Φ m refer to the coalescing right eigenvectors of the transfer matrix and j Á j and hÁ; Ái indicate the norm and inner product respectively. When D H ! 0 all vectors involved in the summation are parallel. From the plots, we can immediately recognize the formation of the SIP at frequency ω SP (see horizontal dashed line). In Fig. 7d-f, we again present the same set of plots, which now are associated with the microstrip system of Fig. 1e. The scattering parameters have been extracted from COMSOL and transformed to the transfer matrix of the unit-cell. A behavior similar to the one found in the CMT model is evident, confirming that the SIP (indicated by a horizontal black dashed line) corresponds to an EPD.
Reflectance in the presence of disorder. The computation of the reflectance in the presence of disorder (see Fig. 4) is performed by considering scattering from a large, but finite, sample. Such a system, in principle is modeled by the temporal CMT equations, where Ψ j i represents the field inside the system and jS ± i represents the timedependent incoming/outgoing signals to/from the system. H eff is the effective Hamiltonian of the system, while H 0 is the Hamiltonian of the isolated system when it is decoupled from the transmission lines. Finally, Δ ¼ Λ À i 2 DD T is a diagonal matrix describing the self-energy term associated with the presence of the transmission lines which support plane waves with dispersion characteristics ω ¼ ε L þ 2V L cos k. Here, Λ is the purely real normalization matrix accounting for resonant shifts, while the extra loss in the system, due to coupling to the transmission lines, is accounted for by the purely real coupling matrix D. The matrix elements of the self-energy term are Δ ½ jl ¼ V L e ik δ jl ðδ j1 þ δ j2 þ δ j2LÀ1 þ δ j2L Þ.
The stationary solutions of Eq. (28) admit the following form Ψ j i ¼ e Àiωt ψ where we have assumed monochromatic incident waves jS ± i ¼ e Àiωt s ± j i.
Substituting these expressions back in Eq. (28) yields the following equations: By solving the first equation for ψ and then applying the result to the second one, we have where S is the scattering matrix that describes the scattering process, G ¼ ω1 2N À H eff À Á À1 is the Green's function of the system and1 2M is the 2M 2M identity matrix. This formalism allows us to evaluate the reflectance for any incident wavefront. In the case of disordered systems, we have injected an appropriately chosen wavefront that produces the reflectance cusp of Eq. (11) in the corresponding perfect (semi-infinite) case. To mimic the response of the semi-infinite system, a controllable amount of loss is introduced to the model. The amount of loss was chosen in such a way that the corresponding absorption length is smaller than the length of the system. At the same time, we made sure that the loss was gradual along the sample in order to avoid smearing of the cusp due to losses. This scheme allows to prevent reflections from the opposite interface of the finite system and, therefore, leads to the suppression of Fabry-Perot resonances.
In the case of SIP-structures, we distributed the loss following a quadraticallyramped spatial profile, resulting in a loss-growth from γ ¼ 0 to γ ¼ 1:5 10 À2 over L ¼ 1000 unit cells. Our analysis indicated that the RBE is more susceptible to local losses than the SIP. For this reason, in the case of the RBE, we have distributed the quadratically-ramped loss in a way that it increases with a much slower rate over a larger total distance (L ¼ 5000 unit cells) from γ ¼ 0 to γ ¼ 1:5 10 À3 .
Phenomenological description of sensitivity bound due to multiplicative noise. The scaling of the sensitivity bound calls for a general argument for its explanation. The following heuristic argument provides some understanding of the power-law for ν c $ W. To this end, we consider a photon propagating inside the semi-infinite structure with velocity v g . During time τ W ; related to the Wigner-Smith delay time, the photon will be propagating a distance ξ 1 ¼ v g τ W inside the structure where ξ 1 is the localization length due to the presence of disorder. On the other hand, the presence of the SIP singularity is resolved in a scattering experiment, whenever the incident wave interacts with the sample for times that are, at least, inversely proportional to the detuning from ω SP i.e., τ W $ 1=ν. Substituting this estimation, together with the expression Eq. (7) for the group velocity (for m ¼ 3) we get that ν c $ ð 1 ξ 1 Þ 3 . At the same time, previous scaling analysis for the localization scaling in the proximity of an SIP-3 indicated that ξ 1 $ W À0:4 ± 0:02 41,42 which eventually gives ν c $ W 1:2 ± 0:06 in reasonable agreement with our numerical findings (see inset of Fig. 4a).
Additive noise analysis. The results presented in Fig. 5 require considering random fluctuations in the incident signal s þ , generated by thermal noise from the transmission lines. The noise at a channel n has been characterized by the autocorrelation function of its amplitude f n where k B is the Boltzmann constant, and T n ¼ T is the temperature characterizing the noise fluctuations which in our simulations has been taken to be the same for all channels. The angular frequency bandwidth B over which the noise is collected by the detector during the data acquisition, along with the continuous-wave signal associated with a carrier frequency, is inversely proportional to the sampling time.
In the Monte-Carlo simulations we have accounted for 10 7 noise realizations for a specific value of the correlation strength σ 2 þ . In our CMT (and COMSOL) modeling, this additional noise term f has been accounted for by assuming that the i th components f i have real and imaginary parts randomly selected from a gaussian distribution with zero mean and variance σ 2 þ ¼ hReðf i Þ 2 i ¼ hI mðf i Þ 2 i ¼ 10 À5 ð5 10 À4 Þ. The noisy outgoing field is evaluated using the scattering matrix formalism of the second subsection of the "Methods", where the S-matrix appearing in Eq. (32) is unitary after having projected out the evanescent modes and renormalizing the remaining elements of the non-unitary scattering matrix e S (see Eq. (18)) with the associated group velocities of the respective modes 43 . The reflectance is calculated as R ¼ s À j j 2 s þ j j 2 . For better statistical processing, we have simulated the reflected signal over an ensemble of 10 7 realizations of the noise. The probability distribution of the extracted detuning PðνÞ in Fig. 5b is calculated by mapping the ensemble of differential reflectances to the corresponding detuning values ν f g, according to the best power-law fit of the mean value hΔRðνÞi.
Finally, in the COMSOL simulations, we have acquired a smooth curve for the reflectance data by performing an additional moving averaging window over a length of 25 data points. At the same time, the resulting detuning error σ ν was averaged over a length of 10 data points.

Data availability
The data used to plot the figures are available at Zenodo (https://doi.org/10.5281/zenodo. 6626113). All other data are available from the corresponding author upon reasonable request.