\section{Introduction}

Earthquake catalogs are a key seismological tool for problems involving location, precursors, and earthquake source parameters \citep{Gutenberg56,Scholz68,AkiRichards,Meier2017}. Catalogs are fundamental to determining where, when, and how a fault has ruptured \citep{AkiRichards,Scholz2002} and they are central to studies of earthquake precursors and seismic hazard based on both field observations \citep{Marsan2014,Bouchon2008,Bouchon2013,Huang2017,Meng2018} and theory \citep{AkiRichards,Ferdowsi2013,Kazemian}. Laboratory studies of earthquake precursors also rely on event catalogs \citep{Lei2014,JohnsonPrecursor,riviere2018}. Such studies include the temporal evolution of the Gutenberg-Richter b-value, which is a fundamental parameter relating laboratory and field studies of earthquake physics \citep{Smith1980,Wang2016,Scholz2015,Gulia,riviere2018,Kwiatek2014}. The seismic b-value has emerged as a key parameter for estimating seismic hazard and for connecting laboratory and field studies of seismicity \citep{Scholz2002}.

Laboratory-based studies of seismic hazard and earthquake forecasting have traditionally relied on earthquake catalogs. However, recent work shows that lab earthquakes can be predicted based on continuous acoustic data \citep{rouet17ML,BRLFriction,Hulbert2017SlowSlip}. These works show that statistical characteristics of the continuous seismic signal emanating from lab fault zones can predict the timing of future failure events as well as the frictional state of the fault zone. Moreover, recent work has extended this approach to field observation by showing that statistical characteristics of continuous seismicity recorded in Cascadia contain a fingerprint of the megathrust displacement rate that can be used to predict the timing of episodic tremor and slip events \citep{RouetCascadia18}. The methods are based on machine learning (ML) and rely exclusively on continuous measurements of acoustic emission (AE). Thus, a key question involves whether this approach can be extended to catalog based measurements of seismic activity.

Applications of ML to seismology and geosciences are becoming increasingly common \citep{lary2016,khawaja18,kenta18,zefeng18,Holtzman18}. Here, we apply the ML Random Forest method \citep{Breiman2001} to study laboratory earthquake catalogs.  We find that as the catalog degrades, so does our ability to infer fault physics, as quantified the ability to predict shear stress and the time to and since failure.  We begin by describing the biaxial shearing experiment, our methods to construct a catalog from the continuous waveform data, the statistical descriptors we extract from the catalog, and the random forest method used to make inferences from these descriptors. This results in successful models for estimating the fault physics variables using local-in-time catalog data.  Then, from the original catalog we construct a series of progressively more degraded catalogs, each of which contains only the events that exceed a variable magnitude of completeness.  We examine the performance of this method as a function of the magnitude of completeness, and observe that the learned signatures of fault physics degrade dramatically if small enough events are not detected. However, the performance plateaus for low magnitudes of completeness; the smallest 80\% events in the catalog can be dropped without sacrificing performance. Our results show that ML methods based on event catalogs can successfully model laboratory faults, with the caveat that predictive power requires a sufficiently complete catalog.

\section{Data Collection}

\subsection{Biaxial Shear Experiment}

We used data from laboratory friction experiments conducted with a biaxial shear apparatus \citep{marone98,JohnsonPrecursor} pictured in Fig.~\ref{figstress}a. Experiments were conducted in the double direct shear configuration in which two fault zones are sheared between three rigid forcing blocks. Our samples consisted of two 5 mm-thick layers of simulated fault gouge with a nominal contact area of 10 x 10 cm$^2$. Gouge material consisted of soda-lime glass beads with initial particle size between 105-149 {\textmu}m. Additional details about the apparatus and sample construction can be found in \citet{Anthony05}, \citet{riviere2018}, and Text S1.  Prior to shearing, we impose a constant fault normal stress of 2 MPa using a servo-controlled load-feedback mechanism and allow the sample to compact. Once the sample has reached a constant layer thickness, the central block is driven down at constant rate of 10~\textmu m/s. In tandem, we collect an AE signal continuously at 4~MHz (Fig.~\ref{figstress}b, yellow curve) from a piezoceramic sensor embedded in a steel forcing block $\approx 22$ mm from the gouge layer.

Figure~\ref{figstress}c shows the shear stress measured on the fault interface for a full experiment.  Experiments begin with a run-in stage where the shear stress increases and macroscopic shearing of the fault zone transitions from stable sliding to stick-slip failure. The repetitive cycles of loading and failure represent laboratory seismic cycles and transition from periodic to aperiodic as a function of load-point displacement in our experiments \citet{Anthony05,JohnsonPrecursor}. Here, we focus on aperiodic slip cycles and measure the slip history for a set of 23 failure events that occur over 268~seconds~\ref{figstress}c. We define large failure events as times for which stress drop exceeds 0.05 MPa within 1 ms.

\begin{figure}[h]
 \centering
 \includegraphics[width=32pc]{fig1.pdf}

 \caption{{\bf a)} Diagram of biaxial shear apparatus. The blue arrow indicates the direction of shear, and the red arrow locates the piezoceramic sensor that records AE. {\bf b)} Example cataloged events (grey circles) derived from continuously recorded AE (yellow) using the statistics of the smoothed signal envelope (black). {\bf c)} Observed shear stress for an experiment at fixed strain rate. Sharp drops in stress correspond to failure events. Inset: The data segment analyzed in this paper. An event catalog is constructed from AE data sampled at 4 MHz. Our models are constructed using event data from the training segment (blue) and evaluated on data from the test segment (green).
 }
 \label{figstress}
 \end{figure}

\subsection{Event Catalog}
\label{sec:events}

\begin{figure}[h]
 \centering
 \includegraphics[width=32pc]{fig2.pdf}
 \caption{ {\bf a)} Normalized count of training events with magnitude exceeding $M$. The dashed line represents ideal Gutenberg-Richter scaling with a $b$-value of 1.66. {\bf b)} Event rates vs. TTF stacked over all 23 failure cycles. {\bf c)} Stacked event statistics for median event magnitude vs. time to failure.}
 \label{figstack}
 \end{figure}

Our event catalog is built from the acoustic signal (example in Fig.~\ref{figstress}b, yellow) during the time interval noted in Fig.~\ref{figstress}c, inset. The catalog is composed of event amplitudes  $A(t_k)$ and times $t_k$, where $k$ is an event index, built following methods described in \citet{riviere2018} and Text S2. We define laboratory magnitudes $m_k = \log_{10}{A(t_k)}$. The catalog contains  $N_\textrm{tot} \approx 3.8 \times 10^5$ slip events with laboratory magnitudes that range from about 1.0 to 4.2. Our lab magnitudes are based only on acoustic signal amplitude, and although they do not connect directly to earthquake magnitude, previous works have demonstrated the connection between seismic b-value and lab frequency magnitude statistics \citep{Scholz2002}.

To visualize the distribution of magnitudes, we construct the cumulative event distribution
\begin{equation}
N(M) = \sum_{m_k \geq M} 1. \label{eq:n_M}
\end{equation}That is, $N(M)$ is the count of all events with magnitude at least $M$. We plot $\log N(M) / N_\textrm{tot}$ in Fig.~\ref{figstack}a and observe reasonable Gutenberg-Richter scaling with a $b$-value of 1.66.

Figure~\ref{figstack}b shows the rate of events as a function of time to failure (TTF), averaged over all slip cycles.  The rate of events increases significantly as failure is approached, similar to what has been observed in prior works \citep{JohnsonPrecursor,weeks1978,Ponomarev1997,riviere2018,Goebel2013}.  Figure~\ref{figstack}c shows that the median event magnitude changes steadily but slowly from approximately $M_{\textrm{median}}=1.4$ to $M_{\textrm{median}}=1.5$  as a function of time during the lab seismic cycle. The trends of the catalog event count and size reflected in Fig.~\ref{figstack}b and Fig.~\ref{figstack}c motivate the creation of machine learning features that describe the evolution event count and size in more detail.

\section{Machine Learning Methods}

\subsection{Regression with Random Forests}

Abstractly, a regression model maps input points $\mathbf{X}_i$ to continuous output labels  $\hat y_i$. Input points are represented by a feature vectors $\mathbf X_i = (X_{i,1}, X_{i,2},\dots X_{i, N_\textrm{features}})$ that are designed to capture the relevant information about $\hat y_i$. The training dataset $\{(\mathbf X_i, y_i)\}$ consists of many points $\mathbf X_i$ associated with true labels $y_i$. From this training set, one aims to learn a model that is successful at predicting labels $\hat y_i$ for {\em new} data points that were unseen within the training process. In this work, the features will be statistics extracted from the catalog, the labels will describe the macroscopic state of the fault (such as shear stress); both are indexed by time window $i$.

We  use the random forest (RF) algorithm \citep{Breiman2001} as it is implemented in the scikit-learn package \citep{pedregosa2011scikit}. A single RF model contains many decision trees, each of which is a simple regression model constructed stochastically from the training data. To make a prediction, the decision tree begins at its root node and asks a series of yes/no questions. Each question has the form: {\em Is some input feature above or below some threshold value?} Eventually a unique leaf node is reached, from which the decision tree makes its prediction. The RF makes a prediction by averaging over many decision trees, and tends to be more robust than any individual tree. Additional details about our training procedure are available in Text S3. We investigated other regression methods \citep{Friedman2001,hui2005,tibshirani96,hoerl1970} besides the RF, and although all were successful to some extent, none offered definitive improvement over the RF (see Table~S1).

A regression model is successful if its predicted labels $\hat y_i$ agree with the true labels $y_i$ for data points in the {\em testing dataset}, which is separate from the training dataset. We use the dimensionless $R^2$ value to quantify the performance of the model,

\begin{equation}
R^2 = 1 - \frac{\sum_i (\hat{y_i}-y_i)^2}{\sum_i (\bar{y}-y_i)^2},
\end{equation}
where $\bar{y}$ is the mean output label for the testing dataset. Let us review the significance of the $R^2$ metric, which is maximized at $R^2=1$, and can take on arbitrarily low values. Using the mean label as a constant predictor for all data points, $\hat y_i = \bar y$, yields $R^2=0$. Thus, $R^2<0$ indicates that a model performs worse than a constant learned only from the labels, without utilizing the features; if $R^2 \leq  0$, a model has failed to discover useful correlations between the features and the labels. The condition $R^2=1$ would indicate that a model makes perfect predictions for all labels in the test set.

\subsection{Feature and Label Creation}

Our task is to build a regression model that uses the event distribution within a {\em local} time window to understand and make predictions about the stress state and times to and from large failure events. Each time window is labeled, and predictions of the labels will depend only on features from within the window. Thus, each test prediction is independent of prior predictions, and the model cannot exploit the quasi-periodicity of the fault state.

There is significant interest in predicting the TTF for an upcoming earthquake event. \citet{rouet17ML} and \citet{BRLFriction} showed that the continuous AE from a fault is remarkably predictive of the TTF and instantaneous friction. In both previous works, $\mathbf X_i$ is a collection of statistical features extracted from the {\em continuous} AE signal for the time window. Here, we build features for each time window using only the event catalog, rather than the continuous acoustic signal.

We divide our full event catalog (cf. Sec.~\ref{sec:events}) into sub-catalogs that contain only the events occurring in non-overlapping time windows of length $\Delta T = 1$ second. The feature vector contains information about events in the sub-catalog for this time window. Our analysis region consists of 270 seconds, which we divide into 1 second windows. We train on the first $60\%$ (159 windows) and test on the remaining $40\%$ (111 windows).

We construct regression models for three labels: (a) the shear stress, (b) the time to the next failure event (TTF), and (c) the time since the last failure event (TSF). The shear stress label is assigned using the mean shear stress over that window. TTF is measured following the end of time window, and TSF is measured preceding the beginning of the time window. For windows including failures, both TTF and TSF are zero.

Abstractly, our the features for each time window describe the cumulative statistics of event counts and amplitudes from the sub-catalog evaluated at magnitude thresholds extracted from the full catalog. More precisely, the feature vector $\mathbf X_i$ for time window $i$ will have components $X_{i,1}, X_{i,2}, \dots X_{i,N_\textrm{features}}$. Each feature component $X_{i,j}$  is associated with a characteristic magnitude $M_j$. For each $j$, we define $M_j$ to be the largest observed magnitude that satisfies
\begin{equation}
\frac{N(M_j)}{N_\textrm{tot}} \geq \alpha^j,
\end{equation}
where the $j$ in right hand side is applied as an exponent; $\alpha$ is a parameter of range $ 0 < \alpha <1$ that controls the fineness of the magnitude bins and the total number of bins, $N_\textrm{features}$. Recall from Eq.~\eqref{eq:n_M} that $N(M)$ is the cumulative event count of our {\em entire} training catalog, and that $N_\textrm{tot} \approx 3.8\times 10^5$. In this work, we select $\alpha=0.7$, and this leads to $N_\textrm{feature} = 35$ non-empty bins. We emphasize that the characteristic magnitudes $M_j$ are derived from the entire training catalog, independent of window index $i$.

We now define the feature vector $\mathbf{X}_i$ for each window $i$ from the sub-catalog of event magnitudes $m^{(i)}_k$, where $k$ indexes the events within the window. We test two simple schemes for $\mathbf X_i$, sensitive to the counts and amplitudes in the catalog. Mathematically, the $j$th components of these two feature vectors are
\begin{align}
X^{\textrm{count}}_{i,j} &=\sum_{m^{(i)}_k \geq M_j} 1 \\
X^{\textrm{ampl}}_{i,j} &= \sum_{m^{(i)}_k \geq M_j} 10^{m^{(i)}_k}  \label{eq:ampl_features}
\end{align}
where the sums run over all sub-catalog events $k$ with magnitude $m^{(i)}_k$ at least $M_j$. In words,  $X_{i,j}^\textrm{count}$ is the count of events in the $i$th window with magnitude exceeding $M_j$ (analogous to $N(M_j)$, but restricted to sub-catalog $i$). The features $X_{i,j}^\textrm{ampl}$ measure the total amplitude of all events in window $i$ above the magnitude $M_j$.

\subsection{Catalog Ablation Test}

To understand how our catalog-RF model performs with decreasing information, we vary the magnitude of completeness of our catalog by discarding all events with magnitude below a cutoff,  $M_\textrm{cut}$. In practice, we implement this ablation by removing features  $X_{i,j}$ for indices  $j$ that correspond to magnitudes $M_j$ below the cutoff $M_{\textrm{cut}}$. We can then quantify the performance as a function of the ablation cutoff $M_\textrm{cut}$.

\section{Results}
\subsection{Performance on Full Catalogs}

Table \ref{tabregression} shows the $R^2$ performance of our catalog-RF models for shear stress, TSF, and TTF, trained using the full catalogs.  The TTF models ($R^2 = 0.551, R^2=0.612$) do not capture fine details. The TTF prediction is worse than the TSF and shear stress models ($R^2 > 0.8$)---in other words, it is harder to predict the future than to estimate aspects of the past or present. Table \ref{tabregression} also indicates that the cumulative counts  $\mathbf{X}^\textrm{count}$ are better for predicting the shear stress, and the cumulative amplitudes  $\mathbf{X}^\textrm{ampl}$ are slightly better for TTF and TSF. The count-based catalog features predict shear stress with an accuracy of $R^2 = 0.898$, nearly equal to the continuous acoustic approach of \citet{BRLFriction}, who report $R^2 = 0.922$ using the same data (that is, same train and test segments on the same experiment). We note that although the same data are analyzed, the methods cannot be compared exactly; for example, the time window in \citet{BRLFriction} is $1.33$ seconds, and each time window has $90\%$ overlap with the previous window. (Recall that in our work the time window length is $1$ second, and windows do not overlap.)  Figure~\ref{figregression} shows the predictions compared to the true labels over time for the amplitude-based features $\mathbf{X}^\textrm{ampl}$, as well as scatter-plots showing a comparison of the true and predicted labels. Figure~S1 shows that models perform nearly as well for shear stress and TSF when analyzing a larger region of 573 s of data (44 failure cycles), but with reduced performance for TTF. This fact can be accounted for by drift in cycle length over time; TTF performance is restored by randomly assigning cycles to testing data from throughout the entire data analysis region (Fig.~S2).

\begin{table}[h]
 \centering
 \begin{tabular}{l  r  r  r }
 \toprule
   & Shear Stress & TSF  & TTF \\
 \midrule
  $\mathbf{X}^\textrm{count}$ & 0.898 & 0.842 & 0.551 \\
  $\mathbf{X}^\textrm{ampl}$ & 0.848 & 0.882 & 0.612 \\
  \bottomrule
 \end{tabular}
 \caption{$R^2$  performance for regression of shear stress, TSF, and TTF for catalog-RF models on complete catalogs. The choice of count or amplitude based features ( $\mathbf{X}^\textrm{count}$ or  $\mathbf{X}^\textrm{ampl}$) affects performance slightly.}
 \label{tabregression}
\end{table}

\begin{figure}[h]
 \centering
 \includegraphics[width=32pc]{fig3.pdf}

 \caption{Regression performance for the catalog-ML model using amplitude-based features. The rows show performance for shear stress, TSF, and TTF, respectively. The left panel of each row shows predictions on training and testing data over time. True values are shown as black squares, and predicted values as colored circles. Time windows for features are contiguous and non-overlapping. Each point is plotted at the end of the time window represented.  The right panel of each row plots the prediction vs. true values for each label. Testing data is shown with dark markers, and training data is shown with faint markers.}
 \label{figregression}
 \end{figure}

\subsection{Performance on Ablated Catalogs}

 \begin{figure}[h]
 \centering
 \includegraphics[width=32pc]{fig4.pdf}

 \caption{The accuracy of the catalog-RF model decreases as training sub-catalogs are artificially limited to varying magnitudes of completeness, $M_\textrm{cut}$. Left: $R^2$ performance vs. $M_{\textrm{cut}}$. The drop of $R^2$ for $M_\textrm{cut}$ between 2.0 and 3.0 indicates that this range of event magnitudes plays a key role in our RF models. The dashed vertical line at $M_{\textrm{cut}} = 3.67$ denotes the smallest event magnitude associated with a large stress drop. Right: $R^2$ is plotted against $N(M_\textrm{cut})/N_\textrm{tot}$, i.e. the fraction of events remaining in the full catalog after ablation. One observes that most events contribute very little to the catalog-RF performance.}
 \label{figablation}
 \end{figure}

Next we build our catalog-RF models on ablated catalogs that include events only above an imposed magnitude of completeness, $M_\textrm{cut}$. For this study, we use the amplitude features $\mathbf{X}^\textrm{ampl}$ of Eq.~\eqref{eq:ampl_features}.  Figure~\ref{figablation}, left side, shows the $R^2$ performance as a function of the imposed magnitude cutoff.  As one might expect, catalog-RF models generally perform worse when provided less information. However, it is interesting to note that $R^2$ decays relatively slowly until $M_\textrm{cut}$ reaches about 2.0, after which the performance drops precipitously, plateauing near magnitude 3.0. In this regard we say that $M_{\textrm{cut}} = 2.0$ provides a {\em sufficient magnitude of completeness} for the catalog-RF model with the given window size.  In Fig.~\ref{figablation}, right side, we plot $R^2$ as a function of the remaining event fraction $N(M_\textrm{cut})/N_\textrm{tot}$. This highlights that the performance drop at $M_{\textrm{cut}} \approx 2.0$ corresponds to ablating $\approx 80\%$ of the events from the catalog. Likewise, the sharp decay in performance indicates that events $M \lesssim 2.5$ are {\em required} to make an accurate determination of the fault state. The performance of the method for $M_{\textrm{cut}} \gtrsim 3 $ plateaus at a non-zero value of $R^2$ because the algorithm is still able to differentiate between windows that contain a large failure and those that do not, but is otherwise unable to discern useful differences between sub-catalogs.

A study of the ablation curve as a function of window size (Fig.~S3) shows that the sharp decay in performance shifts to smaller cutoff magnitudes as the window size is decreased; while regression can be performed with small window sizes, the sufficient magnitude of completeness is lower. This leads us to hypothesize that the sharp decay in performance may be related to the number of events within a window; to keep the number of events in a window fixed, one must lower the magnitude of completeness along with the window size.

\section{Conclusions}

We analyze a dense acoustic event catalog obtained from biaxial shearing friction experiment.  The device produces dozens of aperiodic stick-slip events during an experiment. Stacking catalog statistics over many cycles shows that as failure approaches, events are more common and typically larger (Fig.~\ref{figstack}).  We then use a machine learning workflow to show that physical characteristics of the fault (shear stress, TSF, and TTF) can be obtained from catalog statistics over a short time window using a Random Forest regression model. Previous work employed the continuous AE, and showed remarkable ability to forecast failure \citep{rouet17ML}, as well as to infer the instantaneous shear stress and friction \citep{BRLFriction}.  Here, we demonstrate that the continuous acoustic waveform is not needed if a catalog {\em with sufficient completeness} is available. Our models achieve similar accuracy to the continuous approach.

We note, however, that our catalogs are extraordinarily well resolved. To create these catalogs, we recorded AE at very high frequency (4 MHz) using a sensor very near the fault (distance $\approx 22$ mm). One would not expect such fidelity in real earthquake catalogs, except perhaps if the sensors happened to be located very close to the slip patch.  We find that the smallest events in our catalog contribute little to the prediction accuracy; increasing $M_\textrm{cut}$ from 1.0 to 2.0 only slightly diminishes the catalog-RF model performance (visible in Fig.~\ref{figablation}).

The fact that the smallest events are not {\em required} offers hope that a catalog-based ML approach could eventually contribute to the understanding of natural seismicity. We note a few potential differences between laboratory conditions and natural seismicity that may pose challenges. Laboratory data corresponds to a single isolated fault. Tectonic fault zones are typically comprised of many sub-faults, and predictions may be vastly more difficult. While precursors are frequently observed in laboratory studies, they are not as reliably observed in natural faults. Moreover, large seismic cycles in the earth are far slower \citep{nelson2006,Ellsworth2013}, and thus fewer cycles are available to train ML models.

We note that \citet{RouetCascadia18} have applied the continuous approach developed in the laboratory to slow slip in Cascadia with surprisingly good results.  Thus we are hopeful that the catalog-based approach described here may prove fruitful to the study of real earthquakes, given sufficient catalog quality. It may have particular importance in regions where continuous data are not yet available, where the continuous data is quite noisy, or for examining archival catalogs.

\acknowledgments
Work supported by institutional support (LDRD) at the Los Alamos National Laboratory.  We gratefully acknowledge the support of the Center for Nonlinear Studies. CM was supported by NSF-EAR1520760 and DE-EE0006762. We thank Andrew Delorey, Robert A. Guyer, Bertrand Rouet-Leduc, Claudia Hulbert, and James Theiler for productive discussions. The data used here are freely available from the Penn State Rock Mechanics lab at \url{http://www3.geosc.psu.edu/~cjm38/}.

\section*{Contents}

\begin{enumerate}
\item Text S1: Experimental details
\item Text S2: Catalog generation from acoustic emissions
\item Text S3: Random forest details
\item Table S1: Performance of various regression algorithms
\item Figure S1: Regression for long time segment with contiguous train/test split
\item Figure S2: Regression for long time segment with random cycle train/test split
\item Figure S3: Ablation analysis for varying time-window sizes
\end{enumerate}

\section*{Text S1: Experimental details}

We constructed samples for the double direct shear (DDS) assembly following the procedures detailed in previous works \citep{Anthony05,riviere2018}. Steel guide plates are mounted on the side blocks of the DDS to prevent gouge loss on the front and back edges. The faces of the forcing blocks in contact with the gouge layers are grooved to ensure that shear occurs within the layer; grooves are 0.8 mm deep and spaced every 1mm perpendicular to the direction of shear. In addition, a small rubber jacket covers the bottom half of the sample to minimize material loss throughout the experiment. Forces on faults are measured with strain-gauge load cells placed in series with the horizontal and vertical pistons. Displacements on the fault are measured using direct current displacement transducers that are coupled to the horizontal and vertical pistons. Stress and displacement data are measured continuously at 1kHz throughout the experiment from a 24-bit recording system. All measurements of shear stress and displacement for the DDS configuration are relative to one fault. In addition, we have instrumented 36 broad-band piezoceramic sensors (~0.02-2 MHz) into steel 10 x 10 cm\textsuperscript{2} blocks. The sensors are in blind holes 2-mm from the block face and the block is placed adjacent to the fault zone to monitor acoustic emission activity \citep{riviere2018}. The acoustic emission data analyzed here are recorded continuously at 4 MHz from a single piezoceramic sensor using a 14-bit Verasonics data acquisition system.

\section*{Text S2: Catalog generation from acoustic emissions}

We use a simple event detection algorithm to construct a catalog from the acoustic emission time series signal.  The algorithm has been designed to account for both the response of the sensor and the background noise level. The first step is to compute the envelope of the signal and use a moving average to smooth it over the time scale $\tau_\textrm{smooth}$. From this smoothed envelope we examine peaks and label them as events if they satisfy the following criteria: 1) The peak is above a minimum amplitude threshold, $A_\textrm{min}$, related to the average background noise. 2) The peak is separated from previous peaks by a minimum inter-event time threshold, $\tau_\textrm{min}$. 3) The peak does not occur during in the coda of an event -- this is defined using an exponential ring-down function computed over the ten previous events. Mathematically, event $k$ is not part of a coda if the peak amplitude $A(t_k)$ satisfies
 \begin{equation}
A(t_k) > A(t_{k-i})\exp{\frac{-{(t_k-t_{k-i})}}{\tau_\textrm{ring-down}}}
\label{eq:A_t0}
\end{equation}
for $i = 1...10$.  The parameters of this algorithm are selected by evaluating the catalog quality on random 5~$\mu s$ snapshots of the acoustic signal using trial parameters and visual evaluation of catalog quality. A high-quality catalog maximizes the number of distinct visual phenomena selected as events, while at the same time minimizes the assignment of multiple events to a single phenomenon. In the end, we selected $A_\textrm{min} = 15$, $\tau_\textrm{smooth} = 7.75$ {\textmu}s,  $\tau_\textrm{min} = 150$ {\textmu}s, and $\tau_\textrm{ring-down} = 156.7$ {\textmu}s.

\section*{Text S3: Random forest details}
We employed 100 decision trees per random forest, each of which sees a bootstrap resampling of the training set, with a mean-squared-error cost function. Before training a random forest model (i.e., designing each tree's yes/no questions from the data), one may specify the maximum depth of the decision trees, which controls the allowable complexity of the random forest model. (If the maximum depth is larger, the random forest can fit more complicated datasets, but also becomes more prone to overfitting irrelevant details of the data.) For each regression task, we used 10-fold cross-validation of the training set to select the maximum depth; the possible depths are 1, 2, 4, 6, 8, and 12. The folds for cross-validation are selected by randomly partitioning the training set. Cross-validation uses the $R^2$ score as the criterion for the best model.  After identifying the tree depth which cross-validates most successfully, we retrain a forest with that maximal depth to the entire training set. Other regularization hyperparameters such as the minimal leaf size and random feature sub-selection were not utilized.

\begin{table}[h]
 \centering
\begin{tabular}{llrrrrr}
\toprule
           &       &    RFR &    GBR &  ElasticNet &  Lasso &  Ridge \\
Regression task & Feature vector &        &        &               &          &          \\
\midrule
Shear stress & $\mathbf X^\textrm{count}$ &  {\bf 0.901} &  0.893 &         0.715 &    0.715 &    0.800 \\
           & $\mathbf X^\textrm{ampl}$ &  {\bf 0.860} &  0.827 &         0.817 &    0.817 &    0.796 \\
TSF & $\mathbf X^\textrm{count}$ &  0.843 &  \bf{0.886} &         0.630 &    0.634 &    0.743 \\
           & $\mathbf X^\textrm{ampl}$ &  0.888 &  \bf{0.891} &         0.779 &    0.779 &    0.771 \\
TTF & $\mathbf X^\textrm{count}$ &  0.547 &  \bf{0.579} &         0.533 &    0.533 &    0.571 \\
           & $\mathbf X^\textrm{ampl}$ &  \bf{0.618} &  0.578 &         0.588 &    0.588 &    0.471 \\
\bottomrule
\end{tabular}
\caption{$R^2$ performance for various regression models. We compare random forest regression (RFR) \citep{Breiman2001}, gradient boosting regression (GBR) \citep{Friedman2001}, and linear regression with ElasticNet \citep{hui2005}, lasso \citep{tibshirani96}, and ridge \citep{hoerl1970} regularization penalties as implemented in scikit-learn \citep{pedregosa2011scikit}. The best performing model for each regression task and feature vector type is shown in bold. RFR and GBR produce models of similar quality and both outperform linear regression schemes.}
 \label{methodchoices}
\end{table}

\begin{figure}[h]
 \centering
 \includegraphics[width=32pc]{figS1.pdf}

 \caption{Regression performance for the catalog-ML model using amplitude-based features using 573 seconds of data. The rows show performance for shear stress ($R^2 =0.869$) , TSF ($R^2 =0.813$), and TTF ($R^2 =0.284$), respectively. The left panel of each row shows predictions on training (region shaded blue) and testing (region shaded green) data over time. True values are shown as black squares, and predicted values as colored circles. The right panel of each row plots the prediction vs. true values for each label. Testing data is shown with bold markers, and training data is shown with faint markers.}
 \label{figregression}
 \end{figure}

\begin{figure}[h]
 \centering
 \includegraphics[width=32pc]{figS2.pdf}

 \caption{Regression performance for the catalog-ML model using amplitude-based features using 573 seconds of data, where each failure cycle is individually assigned to either training or testing data. The rows show performance for shear stress ($R^2 =0.881$), TSF ($R^2 =0.869$), and TTF ($R^2 =0.511$), respectively. The left panel of each row shows predictions on training (region shaded blue) and testing (region shaded green) data over time. True values are shown as black squares, and predicted values as colored circles. The right panel of each row plots the prediction vs. true values for each label. Testing data is shown with bold markers, and training data is shown with faint markers.}
 \label{figregression}
 \end{figure}
\begin{figure}[h]
 \centering
 \includegraphics[width=20pc]{figS3.pdf}

 \caption{Ablation analysis for time since failure (TSF) showing $R^2$ performance vs. artificial magnitude of completeness for different window sizes. As the window size decreases, performance decreases. Furthermore, the curve shifts left; the magnitude of completeness necessary to achieve a given $R^2$ performance lowers along with window size. This figure demonstrates that while regression can be performed with small window sizes, the process relies on having information about progressively smaller events.}
 \label{fig:windowsizeabalation}
 \end{figure}


\section{\label{sec:level1}Introduction}

In a previous paper \cite{pendry2017compacted} we commented on the curious mathematical structure of the spectra of singular surfaces which are characterized by three quantum numbers despite the two-dimensional nature of the surface. Transformation optics \cite{Ward1996,Pendry2006} shows how the third dimension is hidden within the singularity \cite{pendry2017compacted}. This paper will explore in detail the nature of the plasmonic excitations, their ability to generate large concentrations of optical energy, and the transition from the discrete excitation spectrum of a non-singular surface to the continuous spectrum of a singular metasurface. Sub-wavelength metal gratings couple external light into surface plasmons, efficiently localizing the electromagnetic energy and finding applications in optical bio-sensing and photovoltaics \cite{Wen2014,Munday2011,Dhawan2011,Bog2012,Kuo2012}. These can be used to control external radiation using metasurfaces \cite{holloway2012overview, lin2014dielectric}. Here we consider a singular metasurface. Singularities in plasmonic systems, such as sharp edges or touching points, concentrate the electromagnetic fields even to sub-nanometric volumes yielding huge energy densities \cite{kelly2003optical, prodan2003hybridization, hao2004electromagnetic, nordlander2004plasmon, lu2005nanophotonic, romero2006plasmons, bukasov2007highly, romero2008plasmon, wu2009optical, aubry2010plasmonic, luo2010surface, benz2016single, galiffi2018broadband}.

Transformation optics takes advantage of the coordinate invariance of Maxwell's equations to give a prescription of how the electromagnetic parameters $\epsilon$ and $\mu$ change under geometrical transformations. For the case of two-dimensional conformal transformations, $\epsilon$ and $\mu$ are left unchanged in the plane, which can be exploited for solving complex plasmonics problems by transforming them to a frame where geometry is simpler \cite{Pendry2012}. This is particularly useful when considering systems with singularities \cite{luo2010surface} as they give rise to divergences that cannot be treated exactly with numerical methods.

In this paper we use transformation optics to derive an analytical theory of the optical response of singular plasmonic metasurfaces. In previous works we have applied this framework to the study of metasurfaces with a smooth shape \cite{Kraft2015, Huidobro2017a}. Here we use the conformal transformation introduced in Ref. \cite{pendry2017compacted} for the design of singular metasurfaces, and present an analytical derivation to completely characterize its optical properties.

The text is structured as follows. First, in Section \ref{sec:singularmetasurface} we present the conformal transformation that generates a metasurface with grooves or wedges forming a sharp angle. The external fields incident on the metasurface are considered in Sec. \ref{sec:source}. Next, in Sec. \ref{sec:Fields} we derive the analytical expressions for the electromagnetic fields and absorption cross section in the metasurface. We then introduce in Sec. \ref{sec:conductivitymodel} a flat surface model from which we obtain an effective surface conductivity that allows us to unambiguously determine the optical response of the metasurface through its reflectivity. Finally, the results for the singular metasurface at normal incidence are presented in Sec. \ref{sec:reflection}, and for two cases which break the symmetry in Sec. \ref{sec:breakingsymmetry}: a symmetric metasurface under oblique incidence and an asymmetric metasurface.

\section{\label{sec:singularmetasurface}Creating a singular metasurface}

We start by describing the transformations that result in the singular metasurfaces shown in Fig. \ref{Transformation}. Let us first consider an array of metal slabs with periodicity along the vertical direction and translational invariance along the horizontal direction, placed in the slab frame ($z_1=x_1+iy_1$). The period of the array is $d$, the thickness of the slabs is $d_3$ and we take $d_1+d_2$ to be the thickness of the dielectric region. An exponential transformation maps the slab array into either a wedge when $d_1+d_2 > d_3$ [panel (a)], or a groove when $d_1+d_2 < d_3$ [see panel (b)], in frame $z_2$. Then an inverse transformation is carried out to get the two-touching-circular segments shown in the $z_3$ frame. As a last step, a logarithmic transformation is used to generate a surface with a periodic set of sharp wedges/grooves (see $z_4$ frame): these are the singular metasurfaces under consideration here.

\begin{figure}
\includegraphics[width=1\columnwidth]{Transformation.png}
\centering
\caption{A series of transformations to generate a singular metasurface with (a) concave shape ($d_1+d_2>d_3$) and (b) convex shape ($d_1+d_2<d_3$). The different coordinate frames are labelled as $z_i$, with $z_i=x_i+i y_i$. In the slab frame, $z_1$, the period of the slab array is $d$, $d_1 + d_2$ is the thickness of the air region and $d_3$ is the thickness of each metal slab. This is successively transformed to a single wedge/groove (frame $z_2$), two-touching-circular segments (frame $z_3$) and finally a singular metasurface (frame $z_4$). Note that x-axis and y-axis do not have the same scale. }
\label{Transformation}
\end{figure}

On the basis of this series of transformations, a one-step transformation from the slab frame to the metasurface frame can be written as
\begin{equation}
\begin{split}
z_4=\frac{T}{2\pi} \mathrm{ln} \bigg( \frac{1}{a(e^{2\pi z_1 /d}-1)}+\frac{1}{2a} \bigg)
\end{split}
\label{transformation_eq}
\end{equation}
Here, $T$ defines the size of the metasurface by fixing its period, and $a$ is chosen as 0.5 so that the singular point in the metasurface frame is located on the y-axis. As mentioned previously, $d=d_1+d_2+d_3$ is the periodicity in the slab frame, and the choice of $\{d_1, d_2, d_3\}$ determines the shape of the metasurface. For $d_1+d_2=d_3$ the transformation generates a flat surface, while $d_1+d_2 > (<) d_3$ results in concave (convex) singular metasurfaces formed by periodic sharp wedges (grooves). In addition, setting $d_1=d_2$ generates a metasurface that is symmetric with respect to the horizontal axis.

In the following we detail our analytical derivations to calculate the optical response of these singular metasurfaces. Throughout the paper we take the metal to be gold with permittivity approximated by the Drude model, $\varepsilon = 1 - \frac{\omega_p^2}{\omega(\omega+i\Gamma)}$, with plasma frequency $\omega_p = 8.95 $ eV/$\hbar$ and damping $\Gamma = 65.8 $ meV/$\hbar$ \cite{novotny2012principles}.

\section{\label{sec:source}Transforming the source}

The problem we set out to solve is that of a p-polarized plane wave (magnetic field out of the plane, $H_z$) incident on the singular metasurface, as any more complex wave front can be expressed as a superposition of plane waves. Since the transformation not only transforms the geometry but also the form of the source, we have to derive the representation of the source in the slab frame. For this purpose, we generate an incident wave using a periodic array of magnetic current line sources (``monopoles'') in the right hand-side of the metasurface frame, as depicted in Fig. 2(b). Taking the sources to be at infinity such that their near fields can be safely neglected \cite{pendry1975chain, tretyakov2003analytical}, their radiated field is a  plane wave incident on the surface. We assume that the period of the metasurface is much less than the free space wavelength. When the incident wave impinges on the singular surface, there will be reflected and transmitted waves. These three sorts of waves all participate in the excitation of surface plasmons at the singular surface. In addition to the source currents, there needs to be a sink at infinity to receive the reflected waves, and another sink at minus infinity to receive the transmitted waves. Then, the source in the slab frame can be obtained by recognising that: (i) a magnetic current line is conserved under the transformation and, (ii) sources at $+\infty$ in the metasurface frame are mapped to the point $z_1 = i n d$ while sources at $-\infty$ in the metasurface frame are mapped to $z_1 = i (n+\frac{1}{2})d$ (here $n$ is an integer). Hence, the monopole sources generating the incident wave and the monopoles receiving the reflected waves are located in the air region in the slab frame, while the monopoles receiving the transmitted wave are placed in the metal region [see Fig. 2(a)].

We start by writing the magnetic field of the incident, reflected and transmitted waves in the metasurface frame as,
\begin{equation}
\left\{
{\begin{array}{l}
H_z^{inc} = H_0 e^{-i k_{0x} x_4 + i k_{0y} y_4}\\
H_z^{ref} = r H_0 e^{i k_{0x} x_4 + i k_{0y} y_4}\\
H_z^{tra} = t H_0 e^{-i k_{0x}^{'}  x_4 + i k_{0y} y_4}
\end{array}}
\right. \label{eq:source}
\end{equation}
where $H_0$ is the wave amplitude, $r$ and $t$ are reflection and transmission coefficients, $k_{0x}= \sqrt{k_0^2-k_{0y}^2}$ and $k_{0x}^{'}= \sqrt{\varepsilon k_0^2-k_{0y}^2}$.

\begin{figure}
\includegraphics[width=0.6\columnwidth]{Source_transformation.png}
\centering
\caption{Sketch showing the source field in both frames. In the metasurface frame (b) the source is a plane wave incident on the surface, which we take to be generated by an array of magnetic line currents located at infinity. The source is mapped into an array of magnetic line currents in the slab frame (a). }
\label{Source transformation}
\end{figure}

In order to write the source field in the slab frame, we transform Eq. \ref{eq:source} using the mapping Eq. (\ref{transformation_eq}) to obtain,
\begin{equation}
\left\{
{\begin{array}{l}
H_z^{inc} = H_0 \bigg(1+i \frac{k_{0x} T}{2\pi} \mathrm{ln} \big( \frac{\pi}{d} \big) \bigg) +  \int\limits_{-\infty}^{\infty} a_a \frac{e^{-|k_x||y_1|}}{|k_x|} e^{i k_x x_1} d k_x +  \int\limits_{-\infty}^{\infty} a_s \frac{e^{-|k_x||y_1|}}{sgn(y_1)k_x} e^{i k_x x_1} d k_x\\
H_z^{ref} = r H_0 \bigg(1-i \frac{k_{0x} T}{2\pi} \mathrm{ln} \big( \frac{\pi}{d} \big) \bigg) -  \int\limits_{-\infty}^{\infty} r a_a \frac{e^{-|k_x||y_1|}}{|k_x|} e^{i k_x x_1} d k_x +  \int\limits_{-\infty}^{\infty} r a_s \frac{e^{-|k_x||y_1|}}{sgn(y_1)k_x} e^{i k_x x_1} d k_x\\
H_z^{tra} = t H_0 \bigg(1 - i \frac{k_{0x}^{'} T}{2\pi} \mathrm{ln} \big( \frac{\pi}{d} \big) \bigg) - \int\limits_{-\infty}^{\infty} t \frac{k_{0x}^{'}}{k_{0x}}  a_a \frac{e^{-|k_x||y_1|}}{|k_x|} e^{i k_x x_1} d k_x  -  \int\limits_{-\infty}^{\infty} t a_s  \frac{e^{-|k_x||y_1+\frac{d}{2}|}}{sgn(y_1+\frac{d}{2})k_x} e^{i k_x x_1} d k_x
\end{array}}
\right.
\end{equation}
The detailed derivation for these source representations is included in Appendix A. In writing the above equations we have assumed that the period of the metasurface is subwavelength ($T\ll |x_4| \ll \lambda$). Also, we have written the fields as a Fourier series and we have identified a symmetric and an antisymmetric component to the source with amplitudes $a_a = -i \frac{k_{0x}T}{4\pi}H_0$ (anti-symmetric source) and $a_s = \frac{k_{0y}T}{4\pi}H_0$ (symmetric source). Note that we define the symmetry of the modes by considering the $E_x$ component in the slab frame: the anti-symmetric (symmetric) mode has odd (even) symmetry of $E_x(y)$. From the above, we can write the source field in k-space for the antisymmetric mode as,
\begin{equation}
    H_z^a(k_x)=
        \left\{
        {\begin{array}{lr}
         (1-r)a_a \frac{e^{-|k_x||y|}}{|k_x|},&{ - {d_2} < y < {d_1}} \\
         -t \frac{k_{0x}^{'}}{k_{0x}}  a_a \frac{e^{-|k_x||y+\frac{d}{2}|}}{|k_x|},&{ - {d_2+d_3} < y < {-d_2}}
        \end{array}}
        \right.
\end{equation}
and for the symmetric mode as,
\begin{equation}
       H_z^s(k_x)= \left\{
        {\begin{array}{lr}
         (1+r)a_s \text{sgn}(k_x) \frac{e^{-|k_x||y|}}{\text{sgn}(y)|k_x|},&{ - {d_2} < y < {d_1}} \\
         -t a_s \text{sgn}(k_x) \frac{e^{-|k_x||y+\frac{d}{2}|}}{\text{sgn}(y+\frac{d}{2})|k_x|},&{ - {d_2+d_3} < y < {-d_2}}
        \end{array}}
        \right.
\end{equation}
in which the constant field components are ignored since they do not contribute to the excitation of surface plasmon polaritons. Finally, note that the source representation includes terms $\sim\frac{e^{-|k_x||y|}}{|k_x|}$, which is just the Fourier transformation of a Hankel function in the quasi-static limit \cite{chew1995waves}. Indeed, a Hankel function is the source representation of a line current, which further confirms our source representation.

\section{\label{sec:Fields}Electromagnetic fields and absorption cross section}

\subsection{\label{sec:dispersion}Surface plasmons dispersion relation}
Once we have the source fields, we can calculate the excited field components. When a SPP mode ($H_z \sim e^{i k_x x}$) on the boundary between metal and air is excited, the total field distribution in the slab frame can be written as
\begin{equation}
H_z(k_x)=\left\{{\begin{array}{lr}
(1-r)a_a \frac{e^{-|k_x||y|}}{|k_x|} + (1+r)a_s \frac{e^{-|k_x||y|}}{\text{sgn}(y)k_x}+ {b_ +}e^{-|k_x|y} + {b_ -}e^{|k_x|y},&{ -d_2 < y < d_1}\\
-t \frac{k_{0x}^{'}}{k_{0x}}  a_a \frac{e^{-|k_x||y+\frac{d}{2}|}}{|k_x|} -t a_s \frac{e^{-|k_x||y+\frac{d}{2}|}}{\text{sgn}(y+\frac{d}{2})k_x} + {c_ +}e^{-|k_x|y} + {c_ -}e^{|k_x|y},&{ - ({d_2} + {d_3}) < y <  - {d_2}}
\end{array}} \right.
\end{equation}
Here, $b_+$, $b_-$, $c_+$ and $c_-$ are the excited mode amplitudes, which have to be determined from the boundary conditions. This requires matching the tangent components of the fields, $H_z$ and $E_x$, at $y=d_1$, $y=-d_2$ and $y=-(d_2+d_3)$ \cite{luo2010surface} and details are given in Appendix B.

\begin{figure}[h]
\includegraphics[width=0.8\columnwidth]{Dispersion_relation.png}
\centering
\caption{Dispersion relations of surface plasmons in singular metasurfaces. (a), (c) SPP dispersion relation in a wedge singular metasurface at oblique (a) and normal incidence (c). The parameters are $d_3=0.1d$, $d_1=d_2=(d-d_3)/2$.(b), (d) SPP dispersion relation in a groove singular metasurface at oblique (b) and normal incidence (d). The parameters are $d_3=0.9d$, $d_1=d_2=(d-d_3)/2$.}
\label{Dispersion relation}
\end{figure}

Once the excited field mode amplitudes have been determined (see Eqs. B2-B5), the dispersion relation of surface plasmons can be obtained by looking at the poles of these coefficients. Neglecting the pole $k_x = 0$ as it is a branch point which corresponds to a localized virtual excitation rather than SPPs \cite{aubry2010broadband, aubry2010conformal}, we concentrate on the plasmon pole. From the amplitude of the anti-symmetric mode (Eqs. B2 and B3), we obtain
\begin{equation}
\begin{split}
(\varepsilon-1)(e^{|k_{x}|(d_1+d_2)} - e^{|k_{x}|d_3}) + (\varepsilon+1)(e^{|k_{x}|d}-1)=0
\end{split}
\end{equation}
While for the symmetric mode (Eqs. B4 and B5) we have,
\begin{equation}
\begin{split}
(\varepsilon-1)(e^{|k_{x}|(d_1+d_2)} - e^{|k_{x}|d_3}) - (\varepsilon+1)(e^{|k_{x}|d}-1)=0
\end{split}
\end{equation}

Figure \ref{Dispersion relation} shows the calculated dispersion relations for the singular wedge (a, c) and groove (b, d) metasurfaces. For the general case, at an oblique incidence, both the symmetric and antisymmetric bands can be excited, as shown in panels (a) and (b). On the other hand, at normal incidence only the antisymmetric band is excited [panels (c) and (d)]. It should be noted that the bands have lower and upper cut-off frequencies different from 0 and $\omega_p$.
For the wedge (groove) metasurface, the symmetric (antisymmetric) band spans a frequency range $\omega_{c1}\le\omega<\omega_{sp}$ while the antisymmetric (symmetric band) spans $\omega_{sp}<\omega\le \omega_{c2}$. Here, $\omega_{sp}$ is the surface plasmon frequency, and $\omega_{c(1,2)}$ are the lower and upper cut-off frequencies, which were shown to be $\omega_{c1}= \omega_{p} \sqrt{\frac{\theta}{2\pi}}$ and $\omega_{c2} = \omega_{p} \sqrt{\frac{2\pi-\theta}{2\pi}}$, where $\theta = 2\pi \frac{d_3}{d}$ for the wedge case ($d_1+d_2>d_3$), $\theta = 2\pi \frac{d_1+d_2}{d}$ for the groove case ($d_1+d_2<d_3$) \cite{luo2010surface}. The appearance of these cut-offs is a result of plasmon hybridization in the infinite periodic array of the slab frame \cite{luo2010surface}. Note that when analyzing the dispersion relation, we have assumed a lossless metal, while we use a finite damping for the rest of the paper.

In the following we focus on the groove singular metasurface, for which the antisymmetric band excited at normal incidence exists below $\omega_{sp}$, where the Drude model gives a more accurate description of the metal.

\subsection{\label{sec:Fieldsrealspace}Electric and magnetic fields in real space}

Next, we calculate the electromagnetic fields in real space. We write the SPP mode field distribution by taking the Fourier transform of the excited terms in the magnetic field given by Eq. (6),
\begin{equation}
    H_z(x,y) = \int\limits_{-\infty}^\infty H_z(k_x,y)e^{i{k_x}x}d{k_x} = 2\pi i \text{Res}\left[H_z(k_x,y)e^{i{k_x}x} \right]|_{k_x=k_{px}}
\end{equation}
where the residue theorem is applied at the plasmon pole, $k_x = k_{px}$. Also, from Eq. (6) we have $H_z (k_x,y) = {b_ +}e^{-|k_x|y} + {b_ -}e^{|k_x|y}$ in the dielectric region and  $H_z (k_x,y) = {c_ +}e^{-|k_x|y} + {c_ -}e^{|k_x|y}$ in the metal region. This yields the following field distribution in the slab frame:
\begin{equation}
    H_z(x,y)= \left\{
        {\begin{array}{lr}
        i 2\pi a ({\Gamma _ + }{e^{ - {\sqrt{k_{px}^2}}y}} + {\Gamma _ - }{e^{{\sqrt{k_{px}^2}}y}})e^{i{k_{px}}|x|},&{ - {d_2} < y < {d_1}}\\
        i 2\pi a ({\Lambda _ + }{e^{ - {\sqrt{k_{px}^2}}y}} + {\Lambda _ - }{e^{{\sqrt{k_{px}^2}}y}})e^{i{k_{px}}|x|},&{ - {d_2+d_3} < y < -{d_2}}
        \end{array}}
        \right.
\end{equation}
where $a$ stands for $a_{(a,s)}$ for the antisymmetric/symmetric mode and $|k_x|$ has been written as $\sqrt{k_x^2}$ in the complex integration. All the field coefficients ($\Gamma_+$, $\Gamma_-$, $\Lambda_+$ and $\Lambda_-$) are given in Appendix B (Eqs. B6-B9). The electric field can be derived from Maxwell's equations using $E_x= \frac{i}{\omega \varepsilon} \frac{\partial{H_z}}{\partial{y}}$,
\begin{equation}
    E_x(x,y) = \left\{
        {\begin{array}{lr}
        \frac{2\pi a \sqrt{k_{px}^2}}{\omega \varepsilon_0}({\Gamma _ + }{e^{ - {\sqrt{k_{px}^2}}y}} - {\Gamma _ - }{e^{{\sqrt{k_{px}^2}}y}})e^{i{k_{px}}|x|},&{ - {d_2} < y < {d_1}}\\
        \frac{2\pi a \sqrt{k_{px}^2}}{\omega \varepsilon_0 \varepsilon}({\Lambda _ + }{e^{ - {\sqrt{k_{px}^2}}y}} - {\Lambda _ - }{e^{{\sqrt{k_{px}^2}}y}})e^{i{k_{px}}|x|},&{ - {d_2+d_3} < y < -{d_2}}
        \end{array}}
        \right.
\end{equation}
and $E_y = -\frac{i}{\omega \varepsilon} \frac{\partial{H_z}}{\partial{x}}$,
\begin{equation}
    E_y(x,y) = \left\{
        {\begin{array}{lr}
        i \text{sgn}(x) \frac{2\pi a k_{px}}{\omega \varepsilon_0}({\Gamma _ + }{e^{ - {\sqrt{k_{px}^2}}y}} + {\Gamma _ - }{e^{{\sqrt{k_{px}^2}}y}})e^{i{k_{px}}|x|},&{ - {d_2} < y < {d_1}}\\
        i \text{sgn}(x) \frac{2\pi a k_{px}}{\omega \varepsilon_0 \varepsilon}({\Lambda _ + }{e^{ - {\sqrt{k_{px}^2}}y}} + {\Lambda _ - }{e^{{\sqrt{k_{px}^2}}y}})e^{i{k_{px}}|x|},&{ - {d_2+d_3} < y < -{d_2}}
        \end{array}}
        \right.
\end{equation}

\begin{figure}
\includegraphics[width=0.6\columnwidth]{modeplot.png}
\centering
\caption{Mode plot of $E_x$ in the metasurface frame for (a) anti-symmetric mode at $\omega=0.6\omega_p$ and (b) symmetric mode at $\omega=0.8\omega_p$. From left column to right column are: charge distribution (sketch), phase and amplitude along the singular surface, and field distribution in the unit cell. The parameters are: $d_3=0.9d$, $d_1=d_2=(d-d_3)/2$. In the color map, red stands for positive values of the field magnitude and blue for negative ones.}
\label{modeplot}
\end{figure}

Once we have the fields in the slab frame, the fields in the metasurface frame can be calculated by mapping them following the rules of transformation optics\cite{Ward1996,Pendry2006}. The obtained mode profile is plotted in Fig. \ref{modeplot} for the antisymmetric (a) and symmetric (b) modes at two frequencies of choice below and above the surface plasmon frequency. We show a sketch of the charge distribution for each mode (antibonding and bonding, respectively), together with the calculated phase and amplitude of the $E_x$ component along the singular surface, and a field plot for $E_x(x,y)$ in the metasurface frame. It is clear the phase oscillates very rapidly, and the in-plane electric field diverges at the singularity. This can be understood from the compression of an infinite dimension hidden at the singular point \cite{pendry2017compacted}, making the derivative of $H_z$ with respect to $x$ and $y$ to be infinite at the singularity. On the other hand, due to its invariance, the magnetic field keeps a finite value in all the frames shown in Fig. \ref{Transformation}, and in particular in the singular metasurface frame.

The presence of loss will attenuate the field in any real system. In this case, the amplitude of the excited SPP wave in the slab frame will be attenuated as it travels along the slab and away from the sources. As was shown previously, there is a critical angle of the groove/wedge where the enhancement by compression and attenuation by loss are balanced \cite{luo2010surface}. The electric field converges for $\theta \leq \theta_c$, while it diverges otherwise, where $\theta_c = \mathrm{Im}[ \mathrm{ln} \left( \frac{\varepsilon-1}{\varepsilon+1} \right) ]$ for the lower band, and $\theta_c = -\mathrm{Im}[ \mathrm{ln} \left( \frac{1-\varepsilon}{\varepsilon+1} \right) ]$ for the upper band. The geometrical parameters chosen in Fig. \ref{modeplot} ($d_3=0.9d$, $d_1=d_2=(d-d_3)/2$) yield an angle at the singularity of $\theta = 2\pi \frac{Min(d_1+d_2,d_3)}{d} = 0.2 \pi$. At the frequencies of choice, $\omega=0.6\omega_p$ for the antisymmetric mode and $\omega=0.8\omega_p$ for the symmetric one, the critical angles are $\theta_c = 0.03$ and $\theta_c = 0.04$, respectively. Therefore, the electric field at singular point of both anti-symmetric and symmetric modes diverges despite the presence material losses, as shown in the right column of Fig. \ref{modeplot}.

\subsection{\label{sec:dissipation}Energy dissipation by SPPs}
Once the field distribution is obtained, we calculate the energy dissipated by the excited SPP mode in the slab frame. Dissipation is due to loss in the metal, so we calculate the absorbed power as the following integral on the slab volume,
\begin{equation}
\begin{split}
    P^{(a,s)}_{abs} &= \int\limits_{slab} \frac{1}{2} \omega \varepsilon_0 \mathrm{Im}[\varepsilon] |E|^2 dxdy\\
    &= \frac{4\pi^2|a|^2|k_{px}|^2 \mathrm{Im}[\varepsilon]}{\omega \varepsilon_0 |\varepsilon|^2} (\frac{|\Lambda_{(a,s)+}|^2}{-2\mathrm{Re}[\sqrt{k_{px}^2}]} (e^{2 \mathrm{Re}[\sqrt{k_{px}^2}]d_2}-e^{2 \mathrm{Re}[\sqrt{k_{px}^2}](d_2+d_3}) \\
    &+ \frac{|\Lambda_{(a,s)-}|^2}{2\mathrm{Re}[\sqrt{k_{px}^2}]}(e^{-2 \mathrm{Re}[\sqrt{k_{px}^2}]d_2}-e^{-2 \mathrm{Re}[\sqrt{k_{px}^2}](d_2+d_3)}) \big) \frac{1}{\mathrm{Im}[k_{px}]}
\end{split}
\end{equation}
where the coefficients $\Lambda_{(a,s)\pm}$ stand for the antisymmetric and symmetric mode coefficients given in Eqs. (B7) and (B9). Using $\Lambda_{a\pm}$ ($\Lambda_{s\pm}$)  yields the power absorbed by the antisymmetric (symmetric) mode. Finally, the absorption cross section of the structure can be obtained by normalizing to the input energy on the system in one period,
\begin{equation}
    \sigma^{(a,s)}_{abs} = \frac{P^{(a,b)}_{abs}}{\frac{1}{2} \sqrt{\frac{\mu_0}{\varepsilon_0}} H_0^2 T \cos{\theta_{in}}},
\end{equation}
where $\theta_{in}$ is incident angle of the plane wave and $H_0$ is the wave amplitude. Thus, we have for the anti-symmetric mode,
\begin{equation}
\begin{split}
    \sigma_{abs}^a
    &= \frac{k_0 T \cos{\theta_{in}}}{2}\frac{|k_{px}|^2 \mathrm{Im}[\varepsilon]}{|\varepsilon|^2} (\frac{|\Lambda_{a+}|^2}{-2\mathrm{Re}[\sqrt{k_{px}^2}]} (e^{2 \mathrm{Re}[\sqrt{k_{px}^2}]d_2}-e^{2 \mathrm{Re}[\sqrt{k_{px}^2}](d_2+d_3}) \\
    &+ \frac{|\Lambda_{a+}|^2}{2\mathrm{Re}[\sqrt{k_{px}^2}]} ( e^{-2 \mathrm{Re}[\sqrt{k_{px}^2}]d_2}-e^{-2 \mathrm{Re}[\sqrt{k_{px}^2}](d_2+d_3)}) \big) \frac{1}{\mathrm{Im}[k_{px}]}
\end{split}
\end{equation}
and for the symmetric one,
\begin{equation}
\begin{split}
    \sigma_{abs}^s
    &= \frac{k_0 T \sin^2{\theta_{in}}}{2\cos{\theta_{in}}}\frac{|k_{px}|^2 \mathrm{Im}[\varepsilon]}{|\varepsilon|^2} (\frac{|\Lambda_{s+}|^2}{-2\mathrm{Re}[\sqrt{k_{px}^2}]} (e^{2 \mathrm{Re}[\sqrt{k_{px}^2}]d_2}-e^{2 \mathrm{Re}[\sqrt{k_{px}^2}](d_2+d_3}) \\
    &+ \frac{|\Lambda_{s+}|^2}{2\mathrm{Re}[\sqrt{k_{px}^2}]} ( e^{-2 \mathrm{Re}[\sqrt{k_{px}^2}]d_2}-e^{-2 \mathrm{Re}[\sqrt{k_{px}^2}](d_2+d_3)}) \big) \frac{1}{\mathrm{Im}[k_{px}]}
\end{split}
\end{equation}

\subsection{\label{sec:blunt}Metasurface with blunt singularities: From a continuous to a discrete spectrum}

Any fabricated metasurface will not show a perfect singularity but a blunt one, so we now move on to treat blunt singularities \cite{luo2012broadband}. A metasurface with rounded singularities maps into a truncated slab array (or truncated cavity array for the groove geometry), as shown in Fig. \ref{Transformation_discretespectrum}. This has an important consequence: when the slabs/cavities are not infinite, the excited SPP modes travelling along the slab will be reflected at its terminals, as depicted in Fig. \ref{Transformation_discretespectrum}(a), resulting in a quantization of the SPP modes supported by the slab structure. Hence, the continuous spectrum of a singular metasurface turns into a discrete spectrum when the singularities are blunt.

\begin{figure}
\includegraphics[width=0.6\columnwidth]{Transformation_discretespectrum.png}
\centering
\caption{Metasurface with blunt singularities. A periodic array of truncated slabs in the slab frame (a) maps to a metasurface with blunt singularities in the metasurface frame (b). }
\label{Transformation_discretespectrum}
\end{figure}

In order to calculate the field distribution and energy dissipation in the truncated slab/cavity arrays we first need to calculate the reflection coefficient of SPPs at the slab/cavity terminal. For conciseness, in the following we will refer only to the cavity array shown in Fig. \ref{Transformation_discretespectrum} (a) and we point out here that the same derivation can be applied to the truncated slab array. The radius of the rounded singularity in the metasurface frame determines the length of the cavities in the slab frame, $L$, as well as the exact shape of the terminal. In order to calculate the reflection coefficient of the SPP modes at the end of the cavity we assume that the terminal is a flat vertical air/metal interface. Then we consider the field of the SPP modes ($H_z^{sp}$) in the cavity, $-L/2<x_1<L/2$, and the field ($H_z^{out}$) in the region outside the cavity, $x_1<-L/2$ and $x_1>L/2$, and impose the continuity condition of the tangential fields and the power flow at $x_1=L/2$ \cite{chandran2012metal, gordon2006light}. This yields the equations,
\begin{eqnarray}
(1+r_{sp})H_z^{sp} &=& H_z^{out} \label{eq:contH}\\
(1-r_{sp})E_y^{sp} &=& E_y^{out} \label{eq:contE}\\
\int\limits_{-(d_2+d_3)}^{d_1}(1-r_{sp}^*)(1+r_{sp})E_y^{sp*} H_z^{sp} dy &=& \int\limits_{-(d_2+d_3)}^{d_1} E_y^{out*} H_z^{out} dy \label{eq:contS}
\end{eqnarray}
where, $H_z^{sp}$ and $E_y^{sp}$ are given in Sec. \ref{sec:Fieldsrealspace}.

In the region outside the periodic cavity array we expand the magnetic field as a series of Bloch waves,
\begin{equation}
\begin{split}
H_z^{out} = \sum_{g} h(g) e^{i g y}
\end{split}
\end{equation}
where $g=n \frac{2\pi}{d}$ and $n$ is an integer. The coefficients $h(g)$ can be expressed as
\begin{equation}
\begin{split}
h(g) &= \frac{1}{d} \int\limits_{-(d_2+d_3)}^{d_1} H_z^{out} e^{-i g y} dy
= \frac{1+r_{sp}}{d} \int\limits_{-(d_2+d_3)}^{d_1} H_z^{sp} e^{-i g y} dy
= \frac{1+r_{sp}}{d} I(g) \label{eq:hg}
\end{split}
\end{equation}
where we have made use of Eq. \ref{eq:contH} and we have named the integral in the last equality as $I(g)$. On the other hand, we can derive the electric field in this region as $E_y^{out} = - \frac{i}{\omega \varepsilon \varepsilon_0} \frac{\partial H_z^{out}}{\partial x}$. This yields
\begin{equation}
\begin{split}
E_y^{out} =  \frac{1}{\omega \varepsilon \varepsilon_0} \sum_{g} \sqrt{\varepsilon k_0^2 - g^2} h(g) e^{i g y}
\end{split}
\end{equation}
where we have used $h(g) \propto e^{i \sqrt{\varepsilon k_0^2 - g^2}x}$ and $\varepsilon$ should be replaced by $1$ for the wedge system. Replacing the expression for $h(g)$, Eq. \ref{eq:hg}, in the above equation we obtain
\begin{equation}
\begin{split}
E_y^{out} = \frac{1+r_{sp}}{\omega \varepsilon \varepsilon_0 d} \sum_{g} \sqrt{\varepsilon k_0^2 - g^2} I(g) e^{i g y}
\end{split}
\end{equation}
Substituting the expressions for $H_z^{out}$ and $E_y^{out}$ into Eq. \ref{eq:contS} and after some algebra we arrive to
\begin{equation}
\begin{split}
\frac{1-r_{sp}}{1+r_{sp}} = \frac{1}{ \omega \varepsilon \varepsilon_0 d} \frac{\sum_{g}  \left( \sqrt{\varepsilon k_0^2 - g^2} \right) |I(g)|^2}{\int\limits_{-(d_2+d_3)}^{d_1} E_y^{sp} H_z^{sp^*} dy} \equiv G
\end{split} \label{eq:G}
\end{equation}
where we take $k_0 = 0$ in the quasi-static limit. From Eq. \ref{eq:G}, we obtain the reflection coefficient as
\begin{equation}
    r_{sp} = \frac{1-G}{1+G} \label{eq:rsp}
\end{equation}

Since in the quasistatic limit radiative loss is small, the amplitude of the reflection coefficient is $|r_{sp}|\approx 1$ (which is indeed confirmed from our calculations using Eq. \ref{eq:rsp}). Hence, we consider only the reflection phase, $r_{sp} = e^{i \phi}$, and write the field of the SPP mode in the cavity of length $L$,
\begin{equation}
    H_z(x,y)= \left\{
        {\begin{array}{lr}
        \begin{split}
            & i 2\pi a ({\Gamma_+}e^{-|k_{px}|y}+ {\Gamma_-}e^{|k_{px}|y}) \\
            & \times (e^{i k_{px} |x|} + e^{-i k_{px} |x|+ i k_{px} L + i\phi}) \frac{1}{1 \mp e^{i k_{px} L + i \phi}}
        \end{split}
        ,&{ - {d_2} < y < {d_1}}\\
        \begin{split}
            & i 2\pi a ({\Lambda_+}e^{-|k_{px}|y}+ {\Lambda_-}e^{|k_{px}|y}) \\
            & \times (e^{i k_{px} |x|} + e^{-i k_{px} |x|+ i k_{px} L + i\phi}) \frac{1}{1 \mp e^{i k_{px} L + i \phi}}
        \end{split}
        ,&{ - {d_2+d_3} < y < -{d_2}}
        \end{array}}
        \right.
\end{equation}
in which $\mp$ stands for $-$ for the anti-symmetric source, and $+$ for the symmetric source. From $H_z$ the electric field is obtained using $E_x(x,y) = \frac{i}{\omega \varepsilon} \frac{\partial{H_z}}{\partial{y}} $ and $E_y(x,y) = -\frac{i}{\omega \varepsilon} \frac{\partial{H_z}}{\partial{x}} $.

Finally, the absorption cross section of the blunt singular metasurface is derived from the electric field following the same procedure detailed in Sec. \ref{sec:dissipation}. We have for the antisymmetric mode,
\begin{equation}
\begin{split}
    \sigma_{abs}^a
    &= \frac{k_0 T \cos{\theta_{in}} }{2} \frac{|k_{px}|^2 \mathrm{Im}[\varepsilon]}{|\varepsilon|^2} \frac{1}{|1-e^{i k_{px} L + i \phi}|^2}
    \bigg[ (\frac{|\Lambda_{a+}|^2}{-2 \mathrm{Re}[\sqrt{k_{px}^2}]} (e^{2 \mathrm{Re}[\sqrt{k_{px}^2}]d_2}-e^{2 \mathrm{Re}[\sqrt{k_{px}^2}](d_2+d_3)}) \\
    & + \frac{|\Lambda_{a-}|^2}{2 \mathrm{Re}[\sqrt{k_{px}^2}]} (e^{-2 \mathrm{Re}[\sqrt{k_{px}^2}]d_2}-e^{-2 \mathrm{Re}[\sqrt{k_{px}^2}](d_2+d_3)})) \frac{1-e^{-2\mathrm{Im}[k_{px}]L}}{\mathrm{Im}[k_{px}]} \\
    &- (\frac{\Lambda_{a+}^* \Lambda_{a-}}{i 2 \mathrm{Im}[\sqrt{k_{px}^2}]} (e^{-i 2 \mathrm{Im}[\sqrt{k_{px}^2}]d_2} - e^{-i 2 \mathrm{Im}[\sqrt{k_{px}^2}](d_2+d_3)}) \\
    &+ \frac{\Lambda_{a+} \Lambda_{a-}^*}{-i 2 \mathrm{Im}[\sqrt{k_{px}^2}]} (e^{i 2 \mathrm{Im}[\sqrt{k_{px}^2}]d_2} - e^{i 2 \mathrm{Im}[\sqrt{k_{px}^2}](d_2+d_3)})) \frac{2 e^{-\mathrm{Im}[k_{px}]L}}{\mathrm{Re}[k_{px}]}(\sin(\mathrm{Re}[k_{px}]L+\phi)-\sin(\phi)) \bigg]
\end{split}
\end{equation}
and for the symmetric one,
\begin{equation}
\begin{split}
    \sigma_{abs}^s
    &= \frac{k_0 T \sin^2{\theta_{in}} }{2 \cos{\theta_{in}}} \frac{|k_{px}|^2 \mathrm{Im}[\varepsilon]}{|\varepsilon|^2} \frac{1}{|1+e^{i k_{px} L + i \phi}|^2}
    \bigg[ (\frac{|\Lambda_{s+}|^2}{-2 \mathrm{Re}[\sqrt{k_{px}^2}]} (e^{2 \mathrm{Re}[\sqrt{k_{px}^2}]d_2}-e^{2 \mathrm{Re}[\sqrt{k_{px}^2}](d_2+d_3)}) \\
    & + \frac{|\Lambda_{s-}|^2}{2 \mathrm{Re}[\sqrt{k_{px}^2}]} (e^{-2 \mathrm{Re}[\sqrt{k_{px}^2}]d_2}-e^{-2 \mathrm{Re}[\sqrt{k_{px}^2}](d_2+d_3)})) \frac{1-e^{-2\mathrm{Im}[k_{px}]L}}{\mathrm{Im}[k_{px}]} \\
    &- (\frac{\Lambda_{s+}^* \Lambda_{s-}}{i 2 \mathrm{Im}[\sqrt{k_{px}^2}]} (e^{-i 2 \mathrm{Im}[\sqrt{k_{px}^2}]d_2} - e^{-i 2 \mathrm{Im}[\sqrt{k_{px}^2}](d_2+d_3)}) \\
    &+ \frac{\Lambda_{s+} \Lambda_{s-}^*}{-i 2 \mathrm{Im}[\sqrt{k_{px}^2}]} (e^{i 2 \mathrm{Im}[\sqrt{k_{px}^2}]d_2} - e^{i 2 \mathrm{Im}[\sqrt{k_{px}^2}](d_2+d_3)})) \frac{2 e^{-\mathrm{Im}[k_{px}]L}}{\mathrm{Re}[k_{px}]}(\sin(\mathrm{Re}[k_{px}]L+\phi)-\sin(\phi)) \bigg]
\end{split}
\end{equation}

\section{\label{sec:conductivitymodel}Modeling a singular metasurface with an effective surface conductivity}

\begin{figure}
\includegraphics[width=0.6\columnwidth]{General_Flat_surface_model.png}
\centering
\caption{Flat surface model: (a) when SPP is not excited, the singular metasurface behaves as a flat surface; (b) when the anti-symmetric SPP mode is excited, the singular metasurface is modeled as a flat surface with an electric surface conductivity $\sigma_e$; (c) when the symmetric SPP mode is excited, the singular metasurface is modeled as a flat surface with a magnetic surface conductivity $\sigma_m$.}
\label{Flat surface model}
\end{figure}

The next step in our analytical treatment is to calculate reflection off the singular metasurface. For this purpose, we have developed an effective surface conductivity model. The model is based on the assumption of subwavelength periodicity and is illustrated in Fig. \ref{Flat surface model}. If there is no SPP mode excited on the metasurface, an incident wave will not see it and as a consequence the singular surface behaves effectively as a flat surface (a). If on the other hand SPPs are excited, energy will be dissipated at the metasurface, which we model through an effective surface conductivity. Due to their different charge distributions at the metasurface, the antisymmetric SPP mode will induce an electric surface current (b), while the symmetric SPP mode will induce a magnetic surface current (c). Hence, we model the antisymmetric mode with an effective electric conductivity, $\sigma_e=\sigma_{er} + i \sigma_{ei}$, and the symmetric mode with an effective magnetic conductivity, $\sigma_m=\sigma_{mr} + i \sigma_{mi}$

Let us first consider the antisymmetric mode. The effective electric surface current generated by the excited SPP modes yields the discontinuity of the tangential magnetic field, while the tangential electric field is continuous,
\begin{eqnarray}
    E_y^{inc}+E_y^{ref}-E_y^{t} &=& 0     \\
    H_z^{inc}+H_z^{ref}-H_z^{t} &=& -\sigma_e E_y^{loc}
\end{eqnarray}
Here, $E_y^{loc}$ is the local tangent electric field in the homogenized system, $E_y^{loc} = \frac{1}{2}(E_y^{inc}+E_y^{ref}+E_y^{tra})$. From the first of the above equations we have $t = \frac{k_{0x} \varepsilon}{k_{0x}^{'}}(1-r)$, so we write
\begin{equation}
    \begin{split}
      E_y^{loc} &= -(1-r)\frac{k_{0x}}{\omega \varepsilon_0} H_0  \label{eq:Eyloc}
    \end{split}
\end{equation}
where we see that the amplitude of the anti-symmetric mode is $\propto 1-r$. In fact, using the proportionality between $t$ and $1-r$ for this mode, we find that the mode coefficients for this mode can be re-written as $\Lambda_{a\pm}=(1-r)\Lambda'_{a\pm}$ and $\Gamma_{a\pm}=(1-r)\Gamma'_{a\pm}$, where the primed coefficients are independent of $r$ and $t$. Detailed expressions for these normalized coefficients are given in Appendix C.

Similarly, for the symmetric mode we have an effective magnetic surface current which yields,
\begin{eqnarray}
    E_y^{inc}+E_y^{ref}-E_y^{t} &=& -\sigma_m H_z^{loc}     \\
    H_z^{inc}+H_z^{ref}-H_z^{t} &=& 0
\end{eqnarray}
where the local tangent magnetic field in the homogenized system is $H_z^{loc} = \frac{1}{2}(H_z^{inc}+H_z^{ref}+H_z^{tra})$. In this case we have $t=1+r$, so
\begin{equation}
    \begin{split}
      H_z^{loc} &= (1+r) H_0  \label{eq:Hzloc}
    \end{split}
\end{equation}
and the amplitude of the symmetric mode is $\propto 1+r$. Normalized mode coefficients can be defined in this case by writing $\Lambda_{s\pm}=(1+r)\Lambda'_{s\pm}$ and $\Gamma_{s\pm}=(1+r)\Gamma'_{s\pm}$. We give detailed expressions in Appendix C.

From the above derivation it is clear that the anti-symmetric mode amplitude is proportional to the local electric field $E_y^{loc}$, while the symmetric mode amplitude is proportional to the local magnetic field $H_z^{loc}$. This justifies the introduction of two kinds of surface conductivities, electric and magnetic ones, in our model in order to mimic energy dissipation by the excited SPP wave. With this, the complex singular metasurface has been simplified as an easy boundary value problem. The reflection coefficient of the singular metasurface can be written straightforwardly as
{\scriptsize
\begin{equation}
\begin{split}
        r=\frac{-4 \varepsilon  \sigma_m + 4 \sigma_e Z_0^2 \cos \theta_{in} \sqrt{\varepsilon -\sin ^2\theta_{in}}-\sigma_e \sigma_m Z_0 \sqrt{\varepsilon -\sin ^2\theta_{in}}+\varepsilon  \sigma_e \sigma_m Z_0 \cos \theta_{in}-4 Z_0 \sqrt{\varepsilon -\sin ^2\theta_{in}}+4 \varepsilon  Z_0 \cos \theta_{in}}{4 \varepsilon  \sigma_m + 4 \sigma_e Z_0^2 \cos \theta_{in} \sqrt{\varepsilon -\sin ^2\theta_{in}}+\sigma_e \sigma_m Z_0 \sqrt{\varepsilon -\sin ^2\theta_{in}}+\varepsilon  \sigma_e \sigma_m Z_0 \cos\theta_{in}+4 Z_0 \sqrt{\varepsilon -\sin ^2\theta_{in}}+4 \varepsilon  Z_0 \cos\theta_{in}}
\end{split}
\label{reflection}
\end{equation}
}
where $\theta_{in}$ is the angle of incidence and $Z_0 = \sqrt{\frac{\mu_0}{\varepsilon_0}}$ is the impedance of free space. Note that for $\sigma_{(e,m)}=0$,  Eq. \ref{reflection} reduces to the reflection coefficient for a flat surface.

The problem then reduces to finding the effective conductivities. We first derive their real parts by using energy conservation, as this is the term that takes away energy. The energy absorbed in the metasurface, $\sigma_{abs} P_{inc}$, must equal energy dissipated by the excited SPP, so we write for the symmetric and antisymmetric modes,
\begin{equation}
\left\{
\begin{split}
\frac{1}{2} \sigma_{er} |E_y^{loc}|^2 T &= \sigma_{abs}^{a} P_{inc}\\
\frac{1}{2} \sigma_{mr} |H_z^{loc}|^2 T &= \sigma_{abs}^{s} P_{inc}
\end{split}
\right.
\label{energy conservation}
\end{equation}
where $P_{inc}=\frac{1}{2} Z_0 H_0^2 T \cos{\theta_{in}}$. By substituting the expression of the local tangent fields, Eqs. \ref{eq:Eyloc} and \ref{eq:Hzloc}, we arrive to
\begin{equation}
\left\{
\begin{split}
    \sigma_{er}  &= \frac{\omega^2 \varepsilon_0^2}{k_{0x}^2} \frac{\sigma^a_{abs}}{|1-r|^2} Z_0 \cos{\theta_{in}} \\
    \sigma_{mr}  &= \frac{\sigma_{abs}^{s}}{|1+r|^2} Z_0 \cos{\theta_{in}}
\end{split}
\right.
\end{equation}
Now, noting that for the antisymmetric mode $t = \frac{k_{0x} \varepsilon}{k_{0x}^{'}}(1-r)$, we have that $\sigma^a_{abs}\propto|1-r|^2$ (see Appendix C for detailed derivations), and we can eliminate the reflection coefficient in the equation for $\sigma_{er}$. Similarly, for the symmetric mode we have $t=1+r$, so $\sigma^s_{abs}\propto|1+r|^2$ and we can also eliminate the reflection coefficient in the equation for $\sigma_{mr}$. We can then write,
\begin{equation}
\left\{
\begin{split}
    \sigma_{er}  &= \frac{\sigma_{abs}^{a'}}{Z_0} = \sigma_{abs}^{a'} \sigma_{e0} \\
    \sigma_{mr}  &= \sigma_{abs}^{s'}Z_0 \sin^2{\theta_{in}} = \sigma_{abs}^{s'} \sigma_{m0} \sin^2{\theta_{in}} \label{eq:sigmar}
    \end{split}
\right.
\end{equation}
where $\sigma_{e0} = Z_0^{-1}$ and $\sigma_{m0} = Z_0$ are the free space electric and magnetic conductivity and we have introduced the intrinsic absorption cross sections (denoted with primes), which do not depend on the incident angle. Hence, the real parts of the effective electric and magnetic conductivities are given by the intrinsic absorption cross sections of the anti-symmetric and symmetric modes, respectively. These read as,
\begin{equation}
\begin{split}
\sigma_{abs}^{(a,s)'} &= \frac{k_0 T}{2} \frac{|k_{px}|^2 \mathrm{Im}[\varepsilon]}{|\varepsilon|^2}
    (\frac{|\Lambda_{(a,s)+}^{'}|^2}{-2\mathrm{Re}[\sqrt{k_{px}^2}]} (e^{2 \mathrm{Re}[\sqrt{k_{px}^2}]d_2}-e^{2 \mathrm{Re}[\sqrt{k_{px}^2}](d_2+d_3}) \\
    &+ \frac{|\Lambda_{(a,s)-}^{'}|^2}{2\mathrm{Re}[\sqrt{k_{px}^2}]} ( e^{-2 \mathrm{Re}[\sqrt{k_{px}^2}]d_2}-e^{-2 \mathrm{Re}[\sqrt{k_{px}^2}](d_2+d_3)}) \big) \frac{1}{\mathrm{Im}[k_{px}]}
\end{split} \label{eq:sigap}
\end{equation}
where all the normalized excited field coefficients are given in Appendix C.

Eqs. \ref{eq:sigmar}-\ref{eq:sigap} determine the real parts of the electric and magnetic conductivities unambiguously. That is, $\sigma_{er}$ and $\sigma_{mr}$ are just functions of the frequency and independent of $r$ and $t$. Furthermore, it should be noted that the electric and magnetic conductivities are defined in the frequency ranges where the antisymmetric and symmetric modes are supported, respectively. For a symmetric metasurface, the symmetric and antisymmetric modes exist over different frequency ranges as we discussed in Sec. \ref{Dispersion relation}. Specifically, we may write for the electric conductivity,
\begin{equation}
    \sigma_{er} = \left\{
        {\begin{array}{lr}
         \sigma_{abs}^{a'} \sigma_{e0}, &{\omega_{c1} < \omega < \omega_{sp}}\\
        0, &{other}
        \end{array}}
        \right.
\label{electric surface conductivity}
\end{equation}
where $\omega_{c1} < \omega < \omega_{sp}$ corresponds to the frequency range when the anti-symmetric mode is excited. Likewise, for the magnetic conductivity we have
\begin{equation}
    \sigma_{mr} = \left\{
        {\begin{array}{lr}
         \sigma_{abs}^{s'} \sigma_{m0} \sin^2{\theta_{in}}, &{\omega_{sp} < \omega < \omega_{c2}}\\
        0, &{other}
        \end{array}}
        \right.
\label{magnetic surface conductivity}
\end{equation}
where $\omega_{sp} < \omega < \omega_{c2}$ corresponds to the frequency range when the symmetric mode is excited.

Finally, in order to fully determine the conductivities, we also need their imaginary parts. These can be obtained through Kramers-Kronig relations as the conductivity must satisfy causality \cite{jackson1999classical, dressel2005electrodynamics}. The imaginary part is thus calculated using
\begin{equation}
\begin{split}
\sigma_{(e,m)i} = -\frac{1}{\pi}P\int\limits_{-\infty}^{\infty}\frac{\sigma_{(e,m)r}(s)}{s-\omega} ds = \frac{1}{\pi}  P \int\limits_{-\infty}^{\infty}\mathrm{ln}\big|s-\omega\big|\frac{d \sigma_{(e,m)r}(s)}{ds} ds
\end{split}
\label{Kramers-Kronig}
\end{equation}
Using the above equation the complex surface conductivities are fully determined and the reflection coefficient is finally obtained by substituting $\sigma_{e}$ and $\sigma_{m}$ into Eq. \ref{reflection}. Furthermore, we note that Eqs. \ref{eq:sigmar} and \ref{Kramers-Kronig} also hold for the metasurface with blunt singularities with the appropriate expressions for $\sigma_{abs}^{(a,b)'}$. These are given in Appendix \ref{sec:appAbsblunt}.

\section{\label{sec:reflection} Reflection spectrum of the singular metasurface}

Using the analytical framework presented in the previous sections, we now discuss the spectrum of a singular groove metasurface of period $T=10$ nm. The metasurface is defined with the parameters $d_3=0.9d$, $d_1=d_2=(d-d_3)/2$, such that it is symmetric with respect to $y=0$. We first consider a normally incident plane wave as source, such that only the anti-symmetric band is excited. Correspondingly, in this scenario the metasurface is modelled by only an effective electric conductivity. The calculated conductivity is presented in Fig. \ref{continuous spectrum} (a). Its real part, shown as a solid blue line, is non-zero only within the anti-symmetric band (between $\omega_{c1}$ and $\omega_{sp}$), as given by Eq. \ref{electric surface conductivity}. Outside the band, no SPPs are excited so the real part of electric surface conductivity is zero, and the metasurface acts effectively as a flat surface without a surface current. The imaginary part of the conductivity is plotted as a dashed red line. In this case, $\sigma_{ei}$ is non-zero also outside the band, where $\sigma_{er}=0$. This is necessary to satisfy Kramers-Kronings relations and represents a phase shift of the reflected wave at the metasurface.

\begin{figure}
\includegraphics[width=0.5\columnwidth]{continuous_spectrum_groove.png}
\centering
\caption{Continuous spectrum of a singular metasurface at normal incidence. (a) Effective electric surface conductivity: real (solid blue line) and imaginary (dashed red) parts. (b) Reflectivity of the metasurface (blue solid line) and of a flat surface with the same permittivity (red dashed). The parameters taken are $T=10$ nm, $d_3=0.9d$, $d_1=d_2=(d-d_3)/2$.}
\label{continuous spectrum}
\end{figure}

From the effective conductivity we obtain the reflection coefficient and we plot the reflectivity ($|r|^2$) in panel (b). The solid blue line corresponds to the calculated reflectivity of the singular metasurface, and the red dashed line to the reflectivity of a flat metal surface calculated from Fresnel coefficients. It is clear that outside of the SPP band the optical response of the singular metasurface is the same as that of a flat surface. On the other hand, between $\omega_{c1}$ and $\omega_{sp}$ the reflectivity presents a continuous spectrum where the reflectivity is smaller than 1. This corresponds to the excitation the anti-symmetric SPP band. Note that we only present analytical results here as full wave simulations using commercial software cannot calculate the spectrum of an exactly singular metasurface because the electric field diverges at the singularity.

\begin{figure}
\includegraphics[width=0.5\columnwidth]{reflection_parameterplot_thetaVSomega_T10nm.jpg}
\centering
\caption{Reflectivity as a function of frequency, $\omega$, and angle of the singular grooves, $\theta$, for a singular metasurface of period $T=10$ nm. }
\label{reflection_parameterplot_thetaVSomega}
\end{figure}

\begin{figure}
\includegraphics[width=0.5\columnwidth]{discrete_spectrum_groove.png}
\centering
\caption{Spectrum of a blunt singular metasurface at normal incidence, showing a discrete rather than a continuum. (a) Phase picked-up by SPPs at the terminal of the truncated cavity. (b) Effective electric surface conductivity, real (solid blue line) and imaginary (red dashed) parts. (c) Reflectivity: analytical (solid blue line) and numerical (red dashed). The metasurface parameters are the same as in Fig. \ref{continuous spectrum}, and the truncation length is $L=d$.}
\label{discrete spectrum}
\end{figure}

Next we explore the dependence of the normal incidence reflectivity spectrum on the sharpness of the singularity. Figure \ref{reflection_parameterplot_thetaVSomega} shows the reflectivity as a function of frequency and of the singularity angle, $\theta = 2\pi\frac{d_1+d_2}{d}$. The blue region corresponds to lower values of reflectivity, where the antisymmetric SPP mode is excited, below $\omega_{sp}$. The low reflectivity band is broad for small angles, and it becomes narrower as $\theta$ gets closer to $180^{\circ}$. The reason for this is that the cut-off frequency for this band, $\omega_{c1} = \omega_{p} \sqrt{\frac{\theta}{2\pi}}$, approaches $\omega_{sp} = \frac{\omega_{p}}{\sqrt{2}}$. On the other hand, we can see that the maximum reflection is reached at a finite angle, $\sim30^{\circ}$.

As discussed previously, real singular metasurfaces will present blunt singularities, which results in a discrete rather than a continuum spectrum. Figure \ref{discrete spectrum} presents results for a metasurface with the same parameters as in Fig. \ref{continuous spectrum} but with blunt singularities. Panel (a) shows the phase, $\phi$, acquired by SPPs reflected off the blunt edge, which we calculate in the slab frame for SPPs reflecting at the truncated end of the cavity. From the calculated reflection phase we obtain the electric surface conductivity through Eqs. \ref{electric surface conductivity}, which is plotted in panel (b). For frequencies within the anti-symmetric band, the conductivity develops resonances which result from the quantization of SPP modes in the periodic array of finite cavities. Then, using Eq. \ref{reflection}, we calculate the corresponding reflectivity [see panel (c), solid blue line]. The reflectivity presents a discrete set of peaks, in contrast to the continuous spectrum of the singular metasurface. We also present the reflectivity obtained from full wave simulations (using the commercial finite element method solver Comsol Multiphysics) as a red dashed line, which shows an excellent agreement with the analytical one. This confirms the validity of our analytical modelling.

\section{\label{sec:breakingsymmetry}Breaking the symmetry}

In Sec. \ref{sec:reflection} we have considered a symmetric metasurface illuminated with a normally incident plane wave, such that the anti-symmetric band is excited below the surface plasmon frequency while the symmetric band (above $\omega_{sp}$) is dark. In this Section we discuss how both bands can be excited by breaking the symmetry. Two ways of breaking the symmetry will be considered. First, we study a symmetric metasurface under oblique incidence, such that the source breaks the symmetry and both bands can be excited. Second, we consider an asymmetric metasurface. Different from the first case, an asymmetry in the geometry mixes the anti-symmetric and symmetric modes so that they become coupled with each other.

\subsection{\label{sec:obliqueinc}Symmetric metasurface under oblique incidence}
When the incident field is coming at an oblique angle of incidence, the source is no longer symmetric with respect to the $y=0$ plane. As a consequence, the dark mode becomes bright and both the anti-symmetric and the symmetric bands are excited. The results obtained with our analytical model for the singular metasurface studied above but under oblique incidence (incident angle $\theta_{in} = 0.4\pi$) are shown in Fig. \ref{continuous_spectrum_obliqueincidence}. Since both anti-symmetric and symmetric modes are excited, there will be both an electric surface conductivity $\sigma_e$ and a magnetic surface conductivity $\sigma_m$. As discussed above, for the groove metasurface the lower band ($\omega_{c1}<\omega<\omega_{sp}$) is the anti-symmetric mode, and its energy dissipation is modelled as a complex electric surface conductivity $\sigma_e$, which is shown in panel (a). Since the conductivity only depends on the intrinsic absorption cross section of SPPs, and not on the incidence angle, it is the same as for the normal incidence case considered above. On the other hand, the upper band is the symmetric mode ($\omega_{sp}<\omega<\omega_{c2}$), which is equivalent to a magnetic surface conductivity $\sigma_m$, shown in panel (b). Finally, panel (c) shows the metasurface reflectivity at oblique incidence. Different from Fig. \ref{continuous spectrum}, where only the lower band is excited, we see how under oblique incidence there are two continuous bands corresponding to the excitation of the anti-symmetric and symmetric modes.

\begin{figure}
\includegraphics[width=0.5\columnwidth]{continuous_spectrum_obliqueincidence_groove.png}
\centering
\caption{Continuous spectrum of the singular metasurface at oblique incidence ($\theta_{in} = 0.4\pi$). (a,b) Effective electric (magnetic) surface conductivity: real (solid blue line) and imaginary (dashed red) parts. (c) Reflectivity of the metasurface (blue solid line) and of a flat surface with the same permittivity (red dashed). The parameters for the metasurface are the same as in Fig. \ref{continuous spectrum}. }
\label{continuous_spectrum_obliqueincidence}
\end{figure}

We next consider the metasurface with blunt singularities under oblique incidence. In this case we first need to calculate the phase picked-up by the SPP modes on reflection at the truncated end in the slab frame. Since both the lower and upper bands are now excited, we calculate $\phi$ not only for the anti-symmetric mode (the same as for normal incidence) but also for the symmetric mode. The SPP phase is shown in Fig. \ref{discrete_spectrum_obliqueincidence}(a), which ranges from the lower cut-off frequency $\omega_{c1}$ to higher cut-off frequency $\omega_{c2}$, and $\omega_{sp}$ clearly separates the two bands. Panels (b) and (c) present the calculated effective electric and magnetic conductivities, respectively, where the quantized modes in the truncated slab yield a discrete spectrum. Finally, we obtain the reflectivity of the metasurface which is shown in panel (d) as a blue solid line, together with results from full wave simulations (red dashed line).

\begin{figure}
\includegraphics[width=1\columnwidth]{discrete_spectrum_obliqueincidence_groove.png}
\centering
\caption{Spectrum of a blunt singular metasurface at oblique incidence, showing a discrete rather than a continuum. (a) Phase picked-up by SPPs at the terminal of the truncated cavity, of length $L=d$. (b,c) Effective electric (magnetic) surface conductivity, real (solid blue line) and imaginary (red dashed) parts. (c) Reflectivity: analytical (solid blue line) and numerical (red dashed). The metasurface parameters are the same as in Fig. \ref{continuous_spectrum_obliqueincidence}.}
\label{discrete_spectrum_obliqueincidence}
\end{figure}

\subsection{\label{sec:level2}Asymmetric metasurface}

We now turn to discussing metasurfaces which are not symmetric with respect to $y=0$, as the one sketched in Fig. \ref{continuous_spectrum_AsymmetricMetasurface}. This asymmetric metasurface can be generated from the slab array by choosing $d_1 \neq d_2$ and, in contrast to the symmetric metasurface, it supports modes that cannot be classified as symmetric and anti-symmetric. In order to treat this case, we proceed in exactly the same way as for the symmetric metasurface and derive the mode coefficients given in Appendix E. The calculated conductivity and reflectivity spectra for a metasurface of the same period as considered previously are shown in Fig. \ref{continuous_spectrum_AsymmetricMetasurface}. Under normal incidence, the asymmetric metasurface can be modelled with an electric conductivity [see panel (a)], while $\sigma_{mr}=0$. However, as a difference with the symmetric metasurface, in this case $\sigma_{er}$ is non-zero both for the lower and upper bands. As a consequence, the calculated reflectivity spectrum (b) shows two continuous bands below and above the surface plasmon frequency, corresponding to the excitation of both the lower and upper bands. It is interesting to note that the reflectivity of the antisymmetric metasurface below $\omega_{sp}$ is very similar to that of the symmetric metasurface, specially for low frequencies. The reason for this is that the reflectivity is mainly determined by the singularity, and the groove angle is the same in both cases.

\begin{figure}
\includegraphics[width=0.5\columnwidth]{continuous_spectrum_AsymmetricGrating_groove.png}
\centering
\caption{Asymmetric singular metasurface at normal incidence. (a) Effective electric surface conductivity: real (solid blue line) and imaginary (dashed red) parts. (b) Reflectivity of the metasurface (blue solid line) and of a flat surface with the same permittivity (red dashed). The parameters used are $T=10nm$, $d_3=0.9d$, $d_1=\alpha (d-d_3)$ and $d_2=(1-\alpha)(d-d_3)$, where the asymmetry factor $\alpha$ is chosen as $0.7$.}
\label{continuous_spectrum_AsymmetricMetasurface}
\end{figure}

Finally, we also present results for an asymmetric metasurface with blunt singularities in Fig. \ref{discrete_spectrum_AsymmetricMetasurface}. Again, when the singularities are not perfect the continuous spectrum of the singular metasurface turns into a discrete spectrum. As before, we first calculate the phase acquired by SPPs in the truncated array in the slab frame [given in panel (a)], and from $\phi$ we calculate the effective electrical surface conductivity (b), which has a non-zero real part both for the lower and upper bands in this case. With this we obtain the reflectivity, which we plot in panel (c) as a solid blue line. The good agreement with results from full wave simulations, shown as a red dashed line, confirms our analytical modelling.

\begin{figure}
\includegraphics[width=0.5\columnwidth]{discrete_spectrum_AsymmetricGrating_groove.png}
\centering
\caption{Spectrum of a blunt asymmetric singular metasurface at normal incidence. (a) Acquired phase by the SPP mode at the terminal of the truncated cavity. (b) Effective surface conductivity: real (solid blue line) and imaginary (dashed red) parts. (c) Reflectivity: analytical (solid blue line) vs numerical (dashed red). The metasurface parameters are the same as in Fig. \ref{continuous_spectrum_AsymmetricMetasurface} and the truncation length is $L=d$. }
\label{discrete_spectrum_AsymmetricMetasurface}
\end{figure}

\section{Conclusions}
In this paper we have presented an analytical theory to study the optical response of singular plasmonic metasurfaces. By means of transformation optics, we have shown that subwavelength metasurfaces with singular grooves (or wedges) are spectrally equivalent to a more symmetric system, a simple periodic array of metal slabs. This has allowed us to obtain analytical expressions for the dispersion relation of SPP modes in the metasurface, as well as for the absorption cross section. Then, by introducing a flat surface model we have derived effective surface conductivities to model energy dissipation by the SPPs, which has enabled us to derive analytical expressions for the metasurface reflectivity. We have shown how singular plasmonic metasurfaces have continuous spectra, with a reduced reflectivity over a broad range of frequencies. On the other hand, realistic metasurfaces will not have perfect singularities, first due to the difficulties in nanofabricating sharp angles, and ultimately to non-locality. We have shown how blunt singularities have a striking effect in the spectrum of metasurfaces. Since blunt singular metasurfaces map to truncated slab arrays, the SPP modes are quantized resulting in a discrete set of peaks in the spectrum. We have also discussed how under normal incidence only the anti-symmetric band is excited, which lies below the surface plasmon frequency for the case of grooves (above for wedges). Breaking the symmetry with an incident wave at an oblique angle allows for the excitation of the symmetric band above the surface plasmon frequency (below for wedges). Finally, we have also discussed how both bands are excited at normal incidence in asymmetric metasurfaces.

For metasurfaces with blunt singularities we have been able to compare our analytical results with full wave electrodynamic simulations, showing that they are in excellent agreement. At the same time, our analytical treatment gives a great physical insight in terms of singularities and allows for an easy optimization of the metasurface shape for a given purpose.

\section{Acknowledgements}
We are grateful for fruitful discussions with Yu Luo. This work was supported by the Gordon and Betty Moore Foundation. F.Y. acknowledges the Lee Family Scholarship for financial support. P.A.H. acknowledges funding from a Marie Sklodowska-Curie Fellowship.