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

### 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 2*M*-dimensional unit cell transfer matrix \({{{{{\mathcal{M}}}}}}\left(\omega \right)\,{{\equiv }}\,{{{{{\mathcal{M}}}}}}(z={z}_{0}+{L}_{0},{z}_{0}=0;\omega )\). This non-normal operator connects the wave amplitudes \({{{{{\boldsymbol{\Phi }}}}}}\) of a monochromatic wave at two different spatial positions of the structure at \(z={z}_{0}+{L}_{0}\) and \({z}_{0}\), through the relation

$${{{{{\mathcal{M}}}}}}\left(z,{z}_{0};\omega \right){{{{{{\boldsymbol{\Phi }}}}}}}_{k}\left({z}_{0}\right)={{{{{{\boldsymbol{\Phi }}}}}}}_{k}\left(z={z}_{0}+{L}_{0}\right)=\lambda \left(\omega \right){{{{{{\boldsymbol{\Phi }}}}}}}_{k}\left({z}_{0}\right);\:{{{{{\rm{where}}}}}}\;\lambda \left(\omega \right)\equiv {{{{{\rm{e}}}}}}^{{{{{{{\rm{i}}}}}}}k{L}_{0}},$$

(2)

where the real part of the Floquet–Bloch wavenumber \({{{{{\rm{Re}}}}}}(k)\in [-{{\pi }}/{L}_{0},-{{\pi }}/{L}_{0}]\) is defined up to a multiple of \(2{{\pi }}/{L}_{0}\) (first Brillouin zone). At the right-hand-side of Eq. (2) we have imposed the Bloch theorem that requires \({{{{{{\boldsymbol{\Phi }}}}}}}_{k}\left(z={z}_{0}+{L}_{0}\right)={{{{{\rm{e}}}}}}^{{{{{{{{\rm{i}}}}}}}}k{L}_{0}}{{{{{{\boldsymbol{\Phi }}}}}}}_{k}({z}_{0})\), reflecting the periodicity of the underlying structure. Finally, the Bloch dispersion relation \(\omega =\omega (k)\) is evaluated by calculating the \(2M\) eigenvalues \(\lambda \left(\omega \right)\) of the transfer matrix by solving the following secular equation:

$${{\det }}\left({{{{{\mathcal{M}}}}}}\left(z,{z}_{0};\omega \right)-{\lambda }_{n}\left(\omega \right)\cdot {\hat{1}}_{2M}\right)=0,\:n=1,\cdots ,2M,$$

(3)

where \({\hat{1}}_{2M}\) is a 2*M*-dimensional identity matrix. When the eigensystem of Eq. (2) has 2\(M\) independent eigensolutions (and under the assumption that the eigensolutions of Eq. (2) can be written in the Bloch form) the transfer matrix \({{{{{\mathcal{M}}}}}}\left(\omega \right)\) is diagonalizable and can be written in the form \({{{{{\mathcal{M}}}}}}\left(\omega \right)={{{{{\mathcal{V}}}}}}\Lambda {{{{{{\mathcal{V}}}}}}}^{-1}\). The nonsingular similarity transformation matrix \({{{{{\mathcal{V}}}}}}\) has columns consisting of the eigenvectors of \({{{{{\mathcal{M}}}}}}\left(z,{z}_{0};\omega \right)\) while \({\Lambda }_{{nm}}={\lambda }_{n}{\delta }_{{nm}}\).

By imposing current conservation, we can show that the transfer matrix \({{{{{\mathcal{M}}}}}}\left(\omega \right)\) satisfies the relation \({{{{{\mathcal{M}}}}}}{\left(\omega \right)}^{{{\dagger}} }\Sigma {{{{{\mathcal{M}}}}}}\left(\omega \right)=\Sigma\) where \(\Sigma ={\Sigma }^{{{\dagger}} }\) is an invertible matrix with \({\Sigma }^{2}=1\), i.e., it belongs to the pseudo-unitary group \(U\left(M,M\right)\). Consequently, \(\left|{{\det }}{{{{{\mathcal{M}}}}}}\left(\omega \right)\right|=1\) while for any frequency \(\omega\) the eigenvalues \({\lambda }_{n}\) satisfy the relations \(\left\{{\lambda }_{n}^{-1}\right\}=\left\{{\lambda }_{n}^{* }\right\}\) corresponding to Bloch wavevectors which are either real (corresponding to propagating waves), or occur in complex conjugate pairs, (corresponding to evanescent modes)^{23,43}. Below, we will refer to the propagating Bloch modes with positive (negative) group velocity \({v}_{{{{{{\rm{g}}}}}}} \, > \, 0\) (\({v}_{{{{{{\rm{g}}}}}}} \, < \,0\)) and to the evanescent Bloch modes with \({{{{{\mathcal{I}}}}}}{{{{{\rm{m}}}}}}\left(k\right) \, > \, 0\) (\({{{{{\mathcal{I}}}}}}{{{{{\rm{m}}}}}}\left(k\right) \, < \,0\)) as forward (backward) waves.

### EPDs in the spectrum of transfer matrices and their relation to stationary points in the Bloch dispersion relation

The non-normal nature of the transfer matrix \({{{{{\mathcal{M}}}}}}\) allows for the formation of NR-EPDs in systems with non-trivial topology. Such degenerate points occur in the spectrum of \({{{{{\mathcal{M}}}}}}(X)\) by appropriately tuning one or a few of its parameters \(X\) (e.g., frequency \(\omega ,\) magnetic field, etc.). In this case, a complete basis is formed by augmenting the Bloch eigenvectors with generalized eigenvectors \({{{{{{\boldsymbol{\Phi }}}}}}}_{q}^{{{{{{\rm{r}}}}}}}\). 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:

$$\left({{{{{\mathcal{M}}}}}}\left({X}_{{{{{{\rm{NR}}}}}}{\mbox{-}}{{{{{\rm{EPD}}}}}}}\right)-{\lambda }_{{{{{{\rm{NR}}}}}}{\mbox{-}}{{{{{\rm{EPD}}}}}}}\cdot {\hat{1}}_{2N}\right){{{{{{\boldsymbol{\Phi }}}}}}}_{q}^{{{{{{\rm{r}}}}}}}\left(z\right)={{{{{{\boldsymbol{\Phi }}}}}}}_{q-1}^{{{{{{\rm{r}}}}}}},\:q=1,\cdots ,m,$$

(4)

where \(m\) (\(2\le m \, < \,2M\)) is the order of the NR-EPD (NR-EPD\(-m\)), \({{{{{{\boldsymbol{\Phi}}}}}}}_{0}^{{{{{{\rm{r}}}}}}}=0\), \({{{{{{\boldsymbol{\Phi }}}}}}}_{1}^{{{{{{\rm{r}}}}}}}={{{{{{\boldsymbol{\Phi }}}}}}}_{{{{{{\rm{NR}}}}}}-{{{{{\rm{EPD}}}}}}}\) is the regular (Bloch) eigenvector while the remaining \(m-1\) eigenvectors \({{{{{{\boldsymbol{\Phi }}}}}}}_{q \, > \, 1}^{r}\) are the generalized eigenvectors. One can show that the generalized eigenvectors diverge along the propagation direction as \({{{{{{\boldsymbol{\Phi }}}}}}}_{q}\left(z\right)\propto {z}^{q-1}{{{{{\rm{e}}}}}}^{{{{{{{{\rm{i}}}}}}}}{kz}}{{{{{{\boldsymbol{\Phi }}}}}}}_{q}(0).\) At an NR-EPD\(-m\), the transfer matrix \({{{{{\mathcal{M}}}}}}\left({X}_{{{{{{\rm{NR}}}}}}{\mbox{-}}{{{{{\rm{EPD}}}}}}}\right)\) is similar to a Jordan form or a matrix containing Jordan blocks, i.e., \({{{{{\mathcal{M}}}}}}\left({X}_{{{{{{\rm{NR}}}}}}{\mbox{-}}{{{{{\rm{EPD}}}}}}}\right)={{{{{\mathcal{V}}}}}}\Lambda {{{{{{\mathcal{V}}}}}}}^{-1}\), and therefore, it is not diagonalizable^{1,2,23,33}. The similarity transformation matrix \({{{{{\mathcal{V}}}}}}\) has columns consisting of the regular and the generalized eigenvectors. In the case when \({{{{{\mathcal{M}}}}}}\) has only one occurrence of an EPD in its spectrum, the matrix \(\Lambda\) consists of a Jordan block of size \(m\) and a diagonal matrix of size \(2M-m\) with diagonal elements being the eigenvalues \({\lambda }_{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}_{{{{{{\rm{NR}}}}}}{\mbox{-}}{{{{{\rm{EPD}}}}}}}\to {X}_{{{{{{\rm{NR}}}}}}{\mbox{-}}{{{{{\rm{EPD}}}}}}}+\nu\) 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 \(\nu\). In other words, when an NR-EPD\(-m\) transfer matrix \({{{{{\mathcal{M}}}}}}\left({X}_{{{{{{\rm{NR}}}}}}{\mbox{-}}{{{{{\rm{EPD}}}}}}}\right),\) which is similar to a matrix containing at least a Jordan block of size \(m\, \times m\), is perturbed as \({{{{{\mathcal{M}}}}}}\left({X}_{{{{{{\rm{NR}}}}}}{\mbox{-}}{{{{{\rm{EPD}}}}}}}+\nu \right){{{{{\mathcal{\approx }}}}}}{{{{{\mathcal{M}}}}}}\left({X}_{{{{{{\rm{NR}}}}}}{\mbox{-}}{{{{{\rm{EPD}}}}}}}\right)+\nu \cdot \bar{\Delta },\) where \(\bar{\Delta }\) is a constant perturbation matrix, then the degenerate eigenvalues \(\lambda\), will obey the fractional expansion

$${\lambda }_{q}={\lambda }_{{{{{{\rm{NR}}}}}}-{{{{{\rm{EPD}}}}}}}+\mathop{\sum }\limits_{n=1}^{\infty }{a}_{n}^{(q)}{\nu }^{n/m}.$$

(5)

The perturbed eigenvalue \({\lambda }_{q}\) in Eq. (5) can be also written as \({\lambda }_{q}={{{{{\rm{exp }}}}}}\left[{{{{{\rm{i}}}}}}\left({k}_{{{{{{\rm{NR}}}}}}{\mbox{-}}{{{{{\rm{EPD}}}}}}}+\delta k\right){L}_{0}\right]\approx {\lambda }_{{{{{{\rm{NR}}}}}}{\mbox{-}}{{{{{\rm{EPD}}}}}}}+{\sum }_{n=1}{c}_{n}{\delta k}^{n}\) which allows us to deduce that the perturbed Bloch wavevector (to the first order correction in detuning \(\nu\)) is \(\delta k\sim \root {m} \of{\nu }\). When the associated detuning \(\nu\) from \({X}_{{{{{{\rm{NR}}}}}}{\mbox{-}}{{{{{\rm{EPD}}}}}}}\) is identified with the frequency variation from the NR-EPD frequency \({{\omega }_{{{{{{\rm{NR}}}}}}{\mbox{-}}{{{{{\rm{EPD}}}}}}}=\omega }_{{{{{{\rm{SP}}}}}}}\), we get:

$$\omega -{\omega }_{{{{{{\rm{SP}}}}}}}=\nu \sim {\left(k-{k}_{{{{{{\rm{SP}}}}}}}\right)}^{m}.$$

(6)

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 slow-light mode with a vanishing group velocity \({v}_{{{{{{\rm{g}}}}}}}\)

$${v}_{{{{{{\rm{g}}}}}}}=\frac{\partial \omega \left(k\right)}{\partial k}\sim {\left|\omega -{\omega }_{{{{{{\rm{SP}}}}}}}\right|}^{\frac{m-1}{m}}.$$

(7)

The most familiar example of a stationary point is the regular band-edge (RBE) corresponding to \(m=2\). Higher order 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^{\prime} s\) that occur at the band-edges, and the odd order NR-EPD\(-m^{\prime} 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}_{{{{{{\rm{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,24,25,26,27,28,29,30,31,32,33,34,35,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 \(\Delta R(\nu )\) when a control parameter \(X\) (e.g., the frequency \(\omega\) of the incident wavefront) is detuned away from the NR-EPD conditions by \(\nu =X-{X}_{{{{{{\rm{NR}}}}}}{\mbox{-}}{{{{{\rm{EPD}}}}}}}\). For the sake of the argument, we assume only one Jordan block of size \(m\) in the similarity matrix \(\Lambda\). 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\ge 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 \({{{{{{\boldsymbol{\Psi }}}}}}}_{{{{{{\rm{I}}}}}}}\left(\omega ,z\right),\) reflected \({{{{{{\boldsymbol{\Psi }}}}}}}_{{{{{{\rm{R}}}}}}}\left(\omega ,z\right)\) and transmitted \({{{{{{\boldsymbol{\Psi }}}}}}}_{{{{{{\rm{T}}}}}}}\left(\omega ,z\right)\) waves, must satisfy the boundary condition

$${{{{{{\boldsymbol{\Psi }}}}}}}_{{{{{{\rm{T}}}}}}}\left(\omega ,z=0\right)={{{{{{\boldsymbol{\Psi }}}}}}}_{{{{{{\rm{I}}}}}}}\left(\omega ,z=0\right)+{{{{{{\boldsymbol{\Psi }}}}}}}_{{{{{{\rm{R}}}}}}}\left(\omega ,z=0\right),$$

(8)

where the reflection coefficient has been absorbed in \({{{{{{\boldsymbol{\Psi }}}}}}}_{{{{{{\rm{R}}}}}}}\). When the system is detuned away from the NR-EPD (i.e., \(\nu \; \ne \; 0\)), the transmitted wave \({{{{{{\boldsymbol{\Psi }}}}}}}_{{{{{{\rm{T}}}}}}}\left(\omega ,z\right)\) inside the slow-light structure can be decomposed into a superposition of the \(M\) forward Bloch modes \({{{{{{\boldsymbol{\Phi }}}}}}}_{n}^{+}\) (see “Methods”). Specifically,

$${{{{{{\boldsymbol{\Psi }}}}}}}_{{{{{{\rm{T}}}}}}}\left(\omega ,z\right)=\mathop{\sum }\limits_{n=1}^{M}{{{{{{\boldsymbol{\Phi }}}}}}}_{n}^{+}(\omega ,z)\,0\le z.$$

(9)

Furthermore, the energy conservation condition requires that the transmitted and reflected waves satisfy the relation \({S}_{{{{{{\rm{I}}}}}}}+{S}_{{{{{{\rm{R}}}}}}}\left(X\right)={S}_{{{{{{\rm{T}}}}}}}\left(X\right)\), where \({S}_{{{{{{\rm{I}}}}}}}\), \({S}_{{{{{{\rm{T}}}}}}}\left(X\right)\), and \({S}_{{{{{{\rm{R}}}}}}}\left(X\right)\) are the energy fluxes of the incident, transmitted and reflected waves. We further assume that the incident wave is normalized as \({S}_{{{{{{\rm{I}}}}}}}=1\). Since evanescent modes do not contribute to the energy flux, \({S}_{{{{{{\rm{T}}}}}}}\) is

$${S}_{{{{{{\rm{T}}}}}}}\propto \mathop{\sum}\limits_{n\in {{{{{\rm{prop}}}}}}}{\left|{{{{{{\boldsymbol{\Phi }}}}}}}_{n}^{+}\left(z\right)\right|}^{2}{v}_{{{{{{\rm{g}}}}}}}^{(n)},$$

(10)

where \({v}_{{{{{{\rm{g}}}}}}}^{(n)}\) 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 \({S}_{{{{{{\rm{T}}}}}}}\propto {W}_{{{{{{\rm{T}}}}}}}{v}_{{{{{{\rm{g}}}}}}}^{{{{{{\rm{slow}}}}}}}\), where \({{W}_{{{{{{\rm{T}}}}}}}\equiv {\left|{{{{{{\boldsymbol{\Phi }}}}}}}_{{{{{{\rm{slow}}}}}}}^{+}\left(z\right)\right|}^{2}\approx \left|{{{{{{\boldsymbol{\Psi }}}}}}}_{{{{{{\rm{T}}}}}}}\left(z\right)\right|}^{2}\) is the energy density carried by the slow mode and \({v}_{{{{{{\rm{g}}}}}}}^{{{{{{\rm{slow}}}}}}}\) is the corresponding group velocity given by Eq. (7). Therefore, the analysis of the exact form of the \({S}_{{{{{{\rm{T}}}}}}}\) singularity collapses to the investigation of the divergence of the energy density \({W}_{{{{{{\rm{T}}}}}}}\sim {\nu }^{-\alpha }\) with respect to the detuning from the NR-EPD.

We first recall that for \(\nu \,\ne\, 0,\) in the specific case of an SIP of order \(m=3\), the excitation \({{{{{{\boldsymbol{\Psi }}}}}}}_{{{{{{\rm{T}}}}}}}\left(z\right)\) can be written as a superposition of a forward slow-propagating Bloch mode \({{{{{{\boldsymbol{\Phi }}}}}}}_{{{{{{\rm{slow}}}}}}}^{+}\left(\omega ,z\right)\) and of a forward evanescent mode \({{{{{{\boldsymbol{\Phi }}}}}}}_{{{{{{\rm{ev}}}}}}}^{+}\left(\omega ,z\right)\) i.e., \({{{{{{\boldsymbol{\Psi }}}}}}}_{{{{{{\rm{T}}}}}}}\left(\omega ,z\right)={{{{{{\boldsymbol{\Phi }}}}}}}_{{{{{{\rm{slow}}}}}}}^{+}\left(\omega ,z\right)+{{{{{{\boldsymbol{\Phi }}}}}}}_{{{{{{\rm{ev}}}}}}}^{+}(\omega ,z)\). Of course, the decaying evanescent contribution to the transmitted wave becomes negligible at a certain distance \({z}_{{{{{{\rm{c}}}}}}}\propto 1/{{{{{\mathcal{I}}}}}}{{{{{\rm{m}}}}}}\left(k\right)\) from the interface and this explains the approximation i.e., \({\left|{{{{{{\boldsymbol{\Phi }}}}}}}_{{{{{{\rm{slow}}}}}}}^{+}\left(z\right)\right|}^{2}\approx {\left|{{{{{{\boldsymbol{\Psi }}}}}}}_{{{{{{\rm{T}}}}}}}\left(z\right)\right|}^{2}\). Nevertheless, the excitation of an evanescent mode is detrimental in the formation of a cusp as it can lead to a finite energy flux \({S}_{{{{{{\rm{T}}}}}}}\propto {\left|{{{{{{\boldsymbol{\Phi }}}}}}}_{{{{{{\rm{slow}}}}}}}^{+}\left(z\right)\right|}^{2}{v}_{{{{{{\rm{g}}}}}}}={{{{{\rm{const}}}}}}.\) under the condition that \({\left|{{{{{{\boldsymbol{\Phi }}}}}}}_{{{{{{\rm{slow}}}}}}}^{+}\left(z\right)\right|}^{2}\sim {\left|{{{{{\rm{\nu }}}}}}\right|}^{-2/3}\). At the same time, we recall that \({\left|{{{{{{\boldsymbol{\Psi }}}}}}}_{{{{{{\rm{T}}}}}}}\left(z\right)\right|}^{2}\approx \left|{{{{{{\boldsymbol{\Phi }}}}}}}_{{{{{{\rm{slow}}}}}}}^{+}\left(z\right)\right|^{2}\) leading to the conclusion that \(\left|{{{\boldsymbol{\Psi}}}}_{{{\mathrm{T}}}}(z)\right|^{2}\sim\left|{{{\mathrm{\nu}}}}\right|^{-2/3}\). Such a field intensity divergence is compatible with the boundary condition Eq. (8) (which dictates that \({{{{{{\boldsymbol{\Psi }}}}}}}_{{{{{{\rm{T}}}}}}}\left(z\right)\) 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., \({{{{{{\boldsymbol{\Phi }}}}}}}_{{{{{{\rm{slow}}}}}}}^{+}(z=0)\approx -{{{{{{\boldsymbol{\Phi }}}}}}}_{{{{{{\rm{ev}}}}}}}^{+}\left(z=0\right)\sim {\left|{{{{{\rm{\nu }}}}}}\right|}^{-1/3}.\) 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}_{{{{{{\rm{T}}}}}}}\propto {\left|{{{{{\rm{\nu }}}}}}\right|}^{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

$$\Delta R\left(\nu \right)\equiv \left|R\left({\omega }_{{{{{{\rm{SP}}}}}}}\right)-R\left({\omega }_{{{{{{\rm{SP}}}}}}}+\nu \right)\right|=T\left({\omega }_{{{{{{\rm{SP}}}}}}}+\nu \right)\sim {\left|\nu \right|}^{\frac{2}{3}},$$

(11)

where we have assumed that in the lossless scenario at the SIP the reflectance is \(R\left({\omega }_{{{{{{\rm{SP}}}}}}}\right)=1\) (since \({S}_{{{{{{\rm{T}}}}}}}\propto {\left|{{{{{\rm{\nu }}}}}}\right|}^{2/3}\mathop{\to }\limits^{\nu \to 0}0\)), and \(T\left(\omega \right)\equiv {S}_{{{{{{\rm{T}}}}}}}(\omega )/{S}_{{{{{{\rm{I}}}}}}}\) and \(R\left(\omega \right)={S}_{{{{{{\rm{R}}}}}}}(\omega )/{S}_{{{{{{\rm{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 \(\nu =X-{X}_{{{{{{\rm{NR}}}}}}-{{{{{\rm{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., \(\Delta R\left(\nu \right)\sim \sqrt{\nu }.\)

### 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\approx 86\) \({{{{{\rm{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 \,\times 4\) transfer matrix \({{{{{\mathcal{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\Omega\) impedance ports and retrieved the associated transfer matrix \({{{{{\mathcal{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}_{{{{{{\rm{SP}}}}}}}=2.00645\;{{{\mathrm{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}_{{{{{{\rm{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 \(\Delta R\) as a function of frequency detuning \(\nu\) from the SIP frequency \({f}_{{{{{{\rm{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 \(\nu\) as \(\Delta R\sim {\left|\nu \right|}^{0.7}\), which is consistent with the theoretical expectations of Eq. (11). We have also checked the validity of the SLR of \(\Delta R\) with respect to other (global) parameter detunings. In Fig. 2c, we report the response of the differential reflectance \(\Delta R\) with respect to magnetic field variations \(\nu \equiv \Delta H\) from the value of the applied magnetic field strength \({H}_{{{{{{\rm{SP}}}}}}}\) at the SIP. The FEM simulations reveal a scaling \(\Delta R\sim {\left|\Delta H\right|}^{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 \(\Delta R\left(\nu \right)\) versus the detuning \(\nu\) 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 \({{{{{{\boldsymbol{\Psi }}}}}}}_{{{{{{\rm{I}}}}}}}\) in this basis, see Eq. (9). Furthermore, a small imaginary part \(({{{{{\mathcal{I}}}}}}{{{{{\rm{m}}}}}}\left({\varepsilon }_{{{{{\mathrm{0,1}}}}}}\right)\sim {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 \({\phi }_{{{{{{\rm{SP}}}}}}}\) and for fixed incident frequency \(\omega ={\omega }_{{{{{{\rm{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 \(\gamma\) with respect to the size of the system \(L\) is expected to be model dependent and determined by the condition that the absorption length \({\xi }_{\gamma }\) 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 \(\Delta 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., \(\Delta R\)) is translated to uncertainty in the measurand (e.g., \(\nu\)), 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 stationary-point-based sensors is provided by analyzing the detuning error \(\sigma (\nu )\) defined as

$${\sigma }_{\nu }=\frac{{\sigma }_{\Delta R}(\nu )}{\chi (\nu )},$$

(12)

where \({\sigma }_{\Delta R}\left(\nu \right)\equiv \sqrt{\left\langle \Delta {R}^{2}\left(\nu \right)\right\rangle -{\left\langle \Delta R\left(\nu \right)\right\rangle }^{2}}\) is the standard deviation of \(\Delta R\left(\nu \right)\) due to noise fluctuations and \(\chi \equiv \frac{\partial \left\langle \Delta R\left(\nu \right)\right\rangle }{\partial \nu }\) is the sensitivity of the sensor, where \(\left\langle \cdot \right\rangle\) indicates a noise averaging. The detuning error Eq. (12) provides an estimation of the uncertainty in the evaluation of the detuning \(\nu\) (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 \(\varepsilon_0 +\delta\varepsilon_0\) and \(\varepsilon_1 +\delta\varepsilon_1\) of each resonant mode of the CMT model (see Eq. (22) in “Methods”). For demonstration purposes, we have assumed that both \(\delta\varepsilon_0\) and \(\delta\varepsilon_1\) are independent random variables taken from a uniform distribution \(\left[-W,W\right]\).

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 \(\Delta R\left(\nu \right)\) (gray circles) versus the frequency detuning \(\nu\) from \({\omega }_{{{{{{\rm{SP}}}}}}}\). The ensembled average differential reflectance \(\left\langle {\Delta} R\left(\nu \right)\right\rangle\) is also shown with a blue solid line. We find that the power law scaling \(\left\langle {\Delta} R\left(\nu \right)\right\rangle \sim {\nu }^{\alpha }\) with the best-fitting value of \(\alpha \approx 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 \(\left\langle \Delta R\left(\nu \to 0\right)\right\rangle \sim {W}^{\beta }\) for very small detunings \({\nu }_{{{{{{\rm{c}}}}}}}\sim {W}^{3\beta /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 \(\beta \approx 2/3\) (see “Methods” for further analysis).

Further statistical processing of \(\Delta R\left(\nu \right)\) allows us to evaluate \({\sigma }_{\Delta R}\left(\nu \right)\) and \(\chi\) and from there, via Eq. (12), the detuning error \({\sigma }_{\nu }\). To get a better appreciation of the robustness of the SIP-3 based sensing with respect to disorder, we have also calculated \({\sigma }_{\nu }\) 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, \({\sigma }_{\nu }\) 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 low-frequency thermal fluctuations in the input channels. For this purpose, we have performed Monte-Carlo simulations in a semi-infinite structure associated with the CMT Hamiltonian (see Eq. (22) in “Methods”). From these simulations we have extracted the standard deviation \({\sigma }_{\Delta R}\) of \(\Delta R\left(\nu \right)\) due to the presence of the additive noise and the sensitivity \(\chi\) of the sensor, which allowed us to evaluate the detuning error \({\sigma }_{\nu }\) 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 \(\Delta R\left(\nu \right)\). The blue solid line indicates the mean value of \(\left\langle \Delta R\left(\nu \right)\right\rangle\) for each specific value of frequency detuning \({{{{{\rm{\nu }}}}}}\). The height of the gray domain surrounding the blue line represents the standard deviation \({\sigma }_{\Delta R}\) of \(\Delta R\left(\nu \right)\) at each value of \({{{{{\rm{\nu }}}}}}\). When comparing the evaluated \({\sigma }_{\Delta R}\) for two distant values of \({{{{{\rm{\nu }}}}}}\) (see black vertical double-sided arrows) it is found that \({\sigma }_{\Delta R}\) does not experience noticeable variations as a function of \({{{{{\rm{\nu }}}}}}\) and remains approximately constant. At the same time, Eq. (12) implies that for \({\sigma }_{\Delta R}\left(\nu \right)\approx {{{{{\rm{const}}}}}}.\), the detuning error will scale inversely proportional to the sensitivity i.e., \({\sigma }_{\nu }\sim {\chi }^{-1}(\nu )\). Therefore, in the limit of small detuning values where the sensitivity is higher, we expect a decreasing detuning error \({\sigma }_{\nu }\). Figure 5b shows the probability density \({{{{{\mathcal{P}}}}}}(\nu )\) of the simulated detuning measurements performed for two distant values of \({{{{{\rm{\nu }}}}}}\) (blue lines). The corresponding values of standard deviation \({\sigma }_{\nu }\) are indicated by double-sided blue arrows. As expected, when the simulated measurements are performed close to the stationary point singularity, the standard deviation \({\sigma }_{\nu }\) decreases, indicating an enhanced precision of the sensor.

A panorama of \({\sigma }_{\nu }\) versus the detuning \(\nu\), 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 \({\sigma }_{\nu }\sim 1/\chi\) (see the blue and red dashed lines corresponding to \({\sigma }_{\nu }{\sim \nu }^{1/3}\) and \({\sigma }_{\nu }\sim {\nu }^{1/2}\), respectively) which is a consequence of the fact that the standard deviation \({\sigma }_{\Delta R}\) remains approximately constant with respect to the detuning from the stationary point frequency \({\omega }_{{{{{{\rm{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.