\section{Introduction}
\IEEEPARstart{I}{n} their seminal work, Vetterli et al. \cite{vetterli} developed a sampling framework for a certain class of nonbandlimited signals that have a finite rate of innovation (FRI) or finite degrees of freedom per unit time/space. The sampling process consists of passing the signal through a suitable kernel, followed by sampling the resulting signal at specific instants. Under certain conditions, the samples thus obtained are sufficient to completely characterize the signal. The degrees of freedom or parameters of the signal are estimated using a suitable reconstruction technique that is coupled to the sampling process. One of the major aspects of the FRI sampling framework is the design of suitable sampling kernels that are realizable and applicable to a larger class of FRI signals.

Consider the 2-D signal
\begin{equation}
\label{eq:2-Dfri_time}
    f(x,y) = \sum_{\ell = 1}^{L}{\gamma_\ell}\,h(x - x_\ell, y - y_\ell),
\end{equation}
which is a sum of scaled and shifted versions of a known function $h(x,y)$, and $\{\gamma_\ell, x_\ell, y_\ell\}_{\ell=1}^{L}$ are the parameters that specify the signal $f(x,y)$. The signal in \eqref{eq:2-Dfri_time} is an FRI signal with $3L$ degrees of freedom. This signal model is frequently encountered in several { \it imaging applications} such as localization microscopy \cite{betzig2006imaging,folling2008fluorescence}, astronomical imaging \cite{ImageAstronomy_Molina, PantinEtAl07}, deflectometry \cite{sudhakar2013compressive}, etc. In this paper, we address the problem of designing suitable sampling kernels for signals of the form \eqref{eq:2-Dfri_time}.  Specifically, we develop a generalized kernel design methodology.  Before proceeding with the design framework, we briefly recall the various 1-D and 2-D sampling kernels proposed in the literature so far.

\subsection{Sampling Kernels}
In the 1-D case, Vetterli et al. \cite{vetterli} proposed infinitely supported sinc and Gaussian kernels, whereas Dragotti et al. \cite{fri_strang} designed a class of compactly supported kernels that reproduce polynomials or exponentials. Tur et al. \cite{eldar_sos} developed alias cancellation conditions and proposed the sum-of-sincs kernel in the frequency domain. Seelamantula and Unser \cite{css}, and Olkkonen and Olkkonen \cite{olkkonen} employed practiaclly realizable kernels derived from resistor-capacitor circuits to sample and reconstruct a stream of Diracs. Recently, Mulleti and Seelamantula \cite{satish_paley_wiener} developed a generalized method for designing 1-D sampling kernels based on the Paley-Wiener theorem, and as a specific construction, they focussed on the class of sum-of-modulated-spline (SMS) kernels.

Following the 1-D sampling framework \cite{vetterli}, Maravi\'c and Vetterli \cite{maravic_2D} developed a sampling and reconstruction framework for 2-D FRI signals with 2-D sinc and Gaussian functions as sampling kernels. Shukla and Dragotti \cite{shukla2007sampling} considered multidimensional FRI signals and proposed polynomial reproducing kernels for sampling convex shapes and polygons. As an application to step-edge detection in images, Baboulaz et al. \cite{baboulaz} developed a local reconstruction scheme with the B-spline kernel. Improving upon \cite{baboulaz}, Hirabayashi et al. \cite{hirabayashi2010spline} designed exponential reproducing kernels and showed that they perform better than polynomial reproducing kernels. Chen et al. \cite{chen20122D} generalized the B-spline sampling kernels for step-edge detection to polygon signal reconstruction using practical sampling kernels. Pan et al. \cite{FRI_Curves_Pan} employed the 2-D sinc sampling kernel and developed a reconstruction scheme for a certain class of parameterizable 2-D curves. Depending upon the number of parameters, the signal model efficiently represents a wider class of curves that are more complex than polygons. Recently, De and Seelamantula \cite{anindita_2D_fri} developed the separable extension of the 1-D non-repeating sum-of-sincs (NR-SoS) kernel \cite{mulleti2014ultrasound} to arrive at 2-D NR-SoS and 2-D SMS kernels.

\subsection{This Paper}
The main contribution of this paper is a generalized framework for designing compactly supported 2-D sampling kernels. The framework allows for the design of separable or nonseparable kernels. To the best of our knowledge, this is the first methodology for the generic design of nonseparable kernels. Starting from a set of alias cancellation conditions that have to be satisfied by a sampling kernel in the frequency domain, admissible sampling kernels with compact support are developed. The characterization of a kernel starting from alias cancellation conditions to enforcing compact support is based on the Paley-Wiener theorem for functions in $\mathbb{R}^d$, $d \geq 1$. We leverage the 1-D kernel design framework proposed in \cite{satish_paley_wiener} and extend it to 2-D. To begin with, we show that a certain separable extension of the 1-D case results in 2-D SMS kernels \cite{anindita_2D_fri}, which have the attractive feature of reproducing separable 2-D polynomial-modulated exponentials. Further, to demonstrate the generalizability of the proposed design framework, we develop a nonseparable sampling kernel. We present simulation results demonstrating exact recovery of a 2-D Dirac stream using one such nonseparable sampling kernel.

\section{The 2-D FRI Sampling and Reconstruction Problems}
A schematic of the kernel-based sampling framework is shown in Fig.~\ref{fig:block}, where the input signal $f$ is passed through a suitable kernel $g$. The resulting signal $\psi(x,y)=(f*g)(x,y)$ is sampled with the sampling intervals $T_{sx}$ and $T_{sy}$ along $x$ and $y$ axes, respectively to get the measurements $\{\psi(n_1T_{sx},n_2T_{sy})\}$, $\{(n_1, n_2)\}\in \mathbb{Z}^2$.

Let $F(s_x,s_y)=\int_{\mathbb{R}^2}f(x,y) e^{-(xs_x+ys_y)} \mathrm{d}x \mathrm{d}y$ be the Laplace transform of $f$, where $s_x=\sigma_x+\mathrm{j}\Omega_x$ and $s_y=\sigma_y+\mathrm{j}\Omega_y$. The 2-D continuous-time Fourier transform (CTFT) of the signal $f$  in~\eqref{eq:2-Dfri_time} is
\begin{equation}
    \label{eq:2-Dfri_freq}
    F(\mathrm{j}\Omega_x,\mathrm{j}\Omega_y) = H(\mathrm{j}\Omega_x,\mathrm{j}\Omega_y)\,\sum_{\ell = 1}^{L}\gamma_\ell e^{-\mathrm{j}(\Omega_x x_\ell + \Omega_y y_\ell)},
\end{equation}
where $H(\mathrm{j}\Omega_x,\mathrm{j}\Omega_y)$ is the 2-D CTFT of $h(x,y)$. Let  $\mathcal{S}$ denote the set of frequency locations in the $\Omega_x-\Omega_y$ plane defined as $\mathcal{S}=\{(k_1\Omega_{0x}, k_2\Omega_{0y})\}_{k_1\in\mathcal{K}_1, k_2\in\mathcal{K}_2}$ for some non-zero $\Omega_{0x}$ and $\Omega_{0y}$, where $\mathcal{K}_1$ and $\mathcal{K}_2$ are sets of contiguous integers chosen suitably based on the model order $L$ and the noise statistics. Throughout the paper, unless specified otherwise, $k_1 \in \mathcal{K}_1$ and $k_2 \in \mathcal{K}_2$. Now, the measurements of $P(\mathrm{j}\Omega_x,\mathrm{j}\Omega_y) = \frac{F(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y)}{H(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y)}$ evaluated on $\mathcal{S}$ are given by
\begin{align}
 P(\mathrm{j}k_1\Omega_{0x},\mathrm{j}k_2\Omega_{0y}) =  \sum_{\ell = 1}^{L}\gamma_\ell e^{-\mathrm{j}(k_1\Omega_{0x} x_\ell + k_2\Omega_{0y} y_\ell)}
 \label{eq:FRI_kwo}.
\end{align}
To avoid singularities, the set $\mathcal{S}$ is chosen such that $H(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y) \neq 0$ on $\mathcal{S}$. The right-hand side of~\eqref{eq:FRI_kwo} is in the form of sum-of-weighted-complex exponentials (SWCEs). The estimation of the parameters $\{x_\ell, y_\ell\}_{\ell = 1}^{L}$ from the measurements of the SWCE form in~\eqref{eq:FRI_kwo} is performed by employing high-resolution spectral estimation (HRSE) techniques~\cite{moses} such as annihilating filter~\cite{prony} applied suitably for 2-D case or 2-D subspace methods~\cite{2DESPRIT_rouquette2001,2DHarmonicRetrieval_vanpoucke1994,2DparameterPaulraj1998,2DESPRIT_haardt1995}. For more details about the application of these techniques, conditions on the minimum number of measurements required for exact recovery in the absence of noise, etc., the reader is referred to~\cite{maravic_2D, maravic2005sampling}.
\begin{figure}[t]
\centering
{\includegraphics[width=4.0in]{block.pdf}}
\caption{Schematic of a 2-D kernel-based sampling scheme.}
\label{fig:block}
\end{figure}

As $h$ is assumed to be known a priori, given the non-aliased samples of $F$ on $\mathcal{S}$, one could exactly recover the parameters $\{ x_\ell, y_\ell\}_{\ell = 1}^{L}$. The sampling kernel $g$ has to be designed in such a way that the non-aliased samples of $F$ on $\mathcal{S}$ can be obtained by the spatial-domain measurements $\{\psi(n_1T_{sx},n_2T_{sy})\}$.  Recently, we \cite{anindita_2D_fri} derived the conditions on the frequency response of the sampling kernel $g$, which are necessary to obtain non-alised samples of $F$ on $\mathcal{S}$:
\begin{align}
\centering
    &G(\mathrm{j}k_1\Omega_{0x},\mathrm{j}k_2\Omega_{0y}) \neq 0, \, \forall\, k_1\in\mathcal{K}_1, k_2\in\mathcal{K}_2,\,\text{and} \label{eq:alias_cancle_1} \\
    &G(\mathrm{j}k_1\Omega_{0x} + \mathrm{j}m_1\Omega_{sx}\, ,\, \mathrm{j}k_2\Omega_{0y} + \mathrm{j}m_2\Omega_{sy}) = 0,\,\quad  \forall\,k_1\in\mathcal{K}_1, k_2\in\mathcal{K}_2, \, \text{and}\, \forall\, m_1, m_2 \in\mathbb{Z}\backslash\{0\}, \label{eq:alias_cancle_2} \\
    & \text{with} \,\, \Omega_{sx}\geq |\mathcal{K}_1|\Omega_{0x}\,\, \text{and}\,\, \Omega_{sy}\geq |\mathcal{K}_2|\Omega_{0y}.
    \label{eq:alias_cancle_3}
\end{align}
The spatial-domain sampling intervals are given by $T_{sx}=\frac{2\pi}{\Omega_{sx}}$ and $T_{sy}=\frac{2\pi}{\Omega_{sy}}$. Also, it was shown that if the kernel $g$ satisfies  (\ref{eq:alias_cancle_1})-(\ref{eq:alias_cancle_3}), then the 2-D discrete-time Fourier transform of $\{\psi(n_1T_{sx},n_2T_{sy})\}$ evaluated on $\mathcal{S}$ gives the non-aliased samples of $F$ on $\mathcal{S}$.

\section{Generalized Design of 2-D Sampling Kernels}
Recently, a generalized method for designing compactly supported sampling kernels for 1-D FRI signals was developed by Mulleti and Seelamantula \cite{satish_paley_wiener}. We extend that framework to the design of 2-D sampling kernels (both separable and nonseparable varieties). The framework consists of two steps: first, to design kernels that satisfy the alias-cancellation conditions (\ref{eq:alias_cancle_1})-(\ref{eq:alias_cancle_2}), and second, enforcing the kernels to be compactly supported.

Consider a sampling kernel $g(x,y)$, that has a rational 2-D Laplace transform of the form
\begin{equation}
    \displaystyle G(s_x,s_y)= \frac{V(s_x,s_y)\,Z(s_x,s_y)}{U(s_x,s_y)},
    \label{eq:gen_samp_kernel}
\end{equation}
where $Z(s_x,s_y)$ is a function that has zeros at $\{n_1\Omega_{0x},n_2\Omega_{0y}\}$, $\forall\,\{(n_1,n_2)\}\, \in\, \mathbb{Z}^2$, $U(s_x,s_y)$ is a function that introduces poles in $G(s_x, s_y)$ such that $G(s_x,s_y)$ satisfies the alias-cancellation conditions (\ref{eq:alias_cancle_1})-(\ref{eq:alias_cancle_2}), and $V(s_x,s_y)$ is a polynomial function that does not have zeros on $\mathcal{S}$. Next, we choose appropriate $Z(s_x,s_y)$ such that the designed kernels have a compact support. To this end, we invoke the 2-D Paley-Wiener theorem \cite{paley_wiener}, which gives the relation between the growth of entire functions of the exponential type (EFET) in the $s$-domain and the support of their time-domain counterparts. Although in higher dimensions, there are several versions of the theorem, we specifically recall the one given by Gel'fand and Shilov \cite[Chapter 4]{gelfand_generalized_functions}, that fits perfectly in the 2-D design framework.

\begin {theorem}
\label{theo:paley_multi}
A function $g(x,y) \in L^2(\mathbb{R}^2)$ is compactly supported over the domain $\mathcal{T}=\{|x| \leq \tau_x, |y| \leq \tau_y\}$ if and only if its Laplace transform $G(s_x, s_y)$ is an EFET, that is, there exist a constant $C$ such that $|G(s_x,s_y)| < Ce^{\tau_x |s_x|+\tau_y |s_y|}$, and $G(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y) \in L^2(\mathbb{R}^2)$.
\end{theorem}
In the case of 1-D, the design of a specific class of 1-D SMS kernels for a particular choice of $Z(s_x)$, $V(s_x)$, and $U(s_x)$ was demonstrated in \cite{satish_paley_wiener}. In the 2-D case, we consider two particular choices of $Z(s_x,s_y)$, one each for separable and nonseparable cases, that result in compactly supported kernels.

\subsection{Separable Kernels}
The design of separable 2-D sampling kernels is a straightforward extension of the 1-D result and is summarized in the following proposition.
\begin{prop}
\label{prop:separable}
Let $Z(s_x,s_y) = \sinh^{r_1}\left(\frac{\pi s_x}{\Omega_{0x}} \right) \sinh^{r_2}\left(\frac{\pi s_y}{\Omega_{0y}}\right)$ and $U(s_x,s_y) = \hspace{-1.0mm} \prod \limits_{\substack{k_1\in\mathcal{K}_1 \\ k_2\in\mathcal{K}_2}} \hspace{-1.0mm} (s_x-\mathrm{j}k_1\Omega_{0x})^{r_1}(s_y-\mathrm{j}k_2\Omega_{0y})^{r_2}$. Then $\displaystyle G(s_x,s_y)$ satisfies the alias-cancellation conditions \eqref{eq:alias_cancle_1} and \eqref{eq:alias_cancle_2}, and $g(x,y)$ is compactly supported to $\Big[-\frac{r_1\,T_{0x}}{2}, \frac{r_1\,T_{0x}}{2}\Big] \times \left[-\frac{r_2\,T_{0y}}{2}, \frac{r_2\,T_{0y}}{2}\right]$, where $T_{0x}=\frac{2\pi}{\Omega_{0x}}$ and $T_{0y}=\frac{2\pi}{\Omega_{0y}}$.
\end{prop}
Using tools such as the partial fraction decomposition and the binomial theorem, it can be shown that the impulse response of the kernel proposed in Proposition \ref{prop:separable} is a sum of modulated $r^{\text{th}}$-order separable polynomial B-splines (denoted by $\beta^r(\cdot)$). The proof of the proposition is given in Appendix \ref{appendix:prop_separable}.

For the choice of
\begin{align}
    V(\mathrm{j}\Omega_{0x},\mathrm{j}\Omega_{0y})= \sum \limits_{\substack{p_1\in\mathcal{K}_1 \\ p_2\in\mathcal{K}_2}}\hspace{-1.0mm}d_{p_1,p_2}\tilde{V}(\mathrm{j}\Omega_{0x},\mathrm{j}\Omega_{0y}),
    \label{eq:separable_V}
\end{align}
where $\tilde{V}(\mathrm{j}\Omega_{0x},\mathrm{j}\Omega_{0y})\hspace{-0.25mm}=\hspace{-1.5mm}\prod \limits_{\substack{k_1\in\mathcal{K}_1\backslash\{p_1\} \\ k_2\in\mathcal{K}_2\backslash\{p_2\}}} \hspace{-2.0mm} (\mathrm{j}\Omega_x-\mathrm{j}k_1\Omega_{0x})^{r_1} (\mathrm{j}\Omega_y-\mathrm{j}k_2\Omega_{0y})^{r_2}$ and with appropriate constants $d_{p_1,p_2}$, we arrive at the special class of compactly supported 2-D SMS kernels proposed in  \cite{anindita_2D_fri}, whose frequency and impulse responses are
\begin{align}
G_{S}(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y) =\hspace{-0.9mm}\sum \limits_{\substack{k_1\in\mathcal{K}_1 \\ k_2\in\mathcal{K}_2}} \hspace{-1.0mm}\text{sinc}^{r_1}\hspace{-1.0mm}\left(\hspace{-0.7mm}\frac{\Omega_x \hspace{-0.5mm}- \hspace{-0.5mm}k_1\,\Omega_{0x}}{\Omega_{0x}}\hspace{-0.7mm}\right)\hspace{-0.7mm}\text{sinc}^{r_2}\hspace{-1.0mm}\left(\hspace{-0.7mm}\frac{\Omega_y \hspace{-0.5mm}-\hspace{-0.5mm} k_2\,\Omega_{0y}}{\Omega_{0y}}\hspace{-0.7mm}\right), \quad \text{and}
 \label{eq:2-D_SMS_freq}
 \end{align}
\begin{align}
g_S(x,y)= \beta^{(r_1-1)}\hspace{-0.7mm}\left(\hspace{-0.7mm}\frac{x}{T_{0x}}\hspace{-0.7mm}\right) \hspace{-0.3mm} \beta^{(r_2-1)}\hspace{-0.7mm}\left(\hspace{-0.7mm}\frac{y}{T_{0y}}\hspace{-0.7mm}\right)
 \hspace{-1.5mm} \sum \limits_{\substack{k_1\in\mathcal{K}_1 \\ k_2\in\mathcal{K}_2}}\hspace{-1.5mm}e^{\mathrm{j}(k_1\Omega_{0x} x+ k_2\Omega_{0y} y)},
 \label{eq:2-D_SMS}
 \end{align}
respectively (Appendix \ref{appendix:kernel_separable_expression}). One could chose different values for $r_1$ and $r_2$ in $g_S(x,y)$ that would result in kernels with different shapes and spatial-domain supports. Figure~\ref{fig:time_dom_separable} shows impulse responses of two such 2-D SMS kernels.\\
\subsubsection{Polynomial and Exponential Reproducing Kernels}
Since the kernel $g_S(x,y)$ in (\ref{eq:2-D_SMS}) is separable, and it was shown in \cite{satish_paley_wiener} that the 1-D SMS kernels satisfy the generalized Strang-Fix conditions \cite{uriguen2013fri}, it is readily seen that the kernel $g_S(x,y)$ could be used to generate a particular class of polynomials/exponentials. Specifically, for separable 2-D SMS kernels, there exist constants $\{c_{k_1,k_2,n_1,n_2}^{i,j}\}$ such that
$\displaystyle\sum_{n_1,n_2\in \mathbb{Z}}c_{k_1,k_2,n_1,n_2}^{i,j}\,g_S \left( \frac{x-n_1T_{sx}}{T_{sx}}, \frac{y-n_2T_{sy}}{T_{sy}} \right) =$
$~\left(\displaystyle\frac{x^i y^j}{T_{sx}^i T_{sy}^j}\right)\, e^{\mathrm{j}k_1\Omega_{0x} x+\mathrm{j}k_2\Omega_{0y} y}$,
for $i \in \llbracket 0, r_1 \rrbracket$ and $j \in \llbracket 0, r_2 \rrbracket$, where $\llbracket a, b \rrbracket$ denotes the set of contiguous integers from $a$ to $b$, both included.

The support of the kernel $g_S(x,y)$ depends on $r_1$ and $r_2$, and is independent of $|\mathcal{K}_1|$ and $|\mathcal{K}_2|$. This is an attractive feature of the 2-D SMS kernels in the sense that they can reproduce exponentials $\{e^{\mathrm{j}k_1\Omega_{0x} x+\mathrm{j}k_2\Omega_{0y} y}\}_{k_1\in \mathcal{K}_1, k_2 \in \mathcal{K}_2}$, wherein the support of the kernel $g_S(x,y)$ is independent of the order.
\begin{figure}[t]
\begin{minipage}[b]{.48\linewidth}
  \centering
  \centerline{\includegraphics[width=0.65\linewidth,trim={2.4cm 0 1cm 0},clip]{time_domain_kernel_L_3_order_3_0.pdf}}
  \centerline{(a)}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.48\linewidth}
  \centering
  \centerline{\includegraphics[width=0.65\linewidth,trim={2.4cm 0 1cm 0},clip]{time_domain_kernel_L_3_order_2_4.pdf}}
  \centerline{(b)}\medskip
\end{minipage}
\hfill
\caption{[Color online] Impulse response of separable 2-D SMS kernel:  $\mathcal{K}_1=\mathcal{K}_2 = \llbracket -L,L\rrbracket$ for $L = 2$, $\Omega_{0x}=\Omega_{0y} = 2\pi$, and (a) $r_1 = 4, r_2 = 1$ and (b) $r_1 = 3, r_2 = 5$.}
\label{fig:time_dom_separable}
\end{figure}
\subsection{Nonseparable Kernels}
We consider a nonseparable $Z(s_x,s_y)$ and a suitable $U(s_x,s_y)$  that results in a nonseparable $G(s_x,s_y)$ and consequently, a compactly supported $g(x,y)$. The following proposition summarizes the result.
\begin{prop}
\label{prop:non_separable}
Let $\displaystyle Z(s_x,s_y)=\sinh\left(\frac{\pi s_x}{\Omega_{0x}} + \frac{\pi s_y}{\Omega_{0y}}\right)\sinh\left(\frac{\pi s_y}{\Omega_{0y}} - \frac{\pi s_x}{\Omega_{0x}}\right)\, \text{and}\\ \,$ $ U(s_x,s_y) \hspace{-0.25mm}= \hspace{-1.5mm}\prod \limits_{\substack{k_1\in\mathcal{K}_1 \\ k_2\in\mathcal{K}_2}}\hspace{-1.5mm}\left(\hspace{-0.7mm}\frac{s_x}{\Omega_{0x}} \hspace{-0.3mm}+\hspace{-0.3mm} \frac{s_y}{\Omega_{0y}} \hspace{-0.3mm}-\hspace{-0.3mm} \mathrm{j}k_1 \hspace{-0.3mm}-\hspace{-0.3mm} \mathrm{j}k_2\hspace{-0.7mm} \right)  \left(\hspace{-0.7mm}\frac{s_y}{\Omega_{0y}} \hspace{-0.3mm}-\hspace{-0.3mm} \frac{s_x}{\Omega_{0x}} \hspace{-0.3mm}-\hspace{-0.3mm} \mathrm{j}k_2 \hspace{-0.3mm}+\hspace{-0.3mm} \mathrm{j}k_1 \hspace{-0.7mm}\right)$. Then $\displaystyle G(s_x,s_y)$ satisfies the alias-cancellation conditions \eqref{eq:alias_cancle_1} and \eqref{eq:alias_cancle_2}, and $g(x,y)$ is supported over $\left[-{T_{0x}}, {T_{0x}}\right] \times \left[-{T_{0y}}, {T_{0y}}\right]$, where $T_{0x}=\frac{2\pi}{\Omega_{0x}}$ and $T_{0y}=\frac{2\pi}{\Omega_{0y}}$.
\end{prop}
The proof involves two steps: showing that (i) $G(s_x,s_y)$ is an EFET; and (ii) $G(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y)\in L^2(\mathbb{R}^2)$, and is provided in Appendix \ref{appendix:prop_nonseparable}. Substituting the functions $Z(s_x,s_y)$ and $U(s_x,s_y)$ of Proposition \ref{prop:non_separable} in (\ref{eq:gen_samp_kernel}), and using the partial fraction decomposition of $\frac{V(s_x,s_y)}{U(s_x,s_y)}$ and rotation property of the 2-D CTFT, it can be shown that (see Appendix \ref{appendix:kernel_nonseparable_expression}) the frequency response of the sampling kernel in Proposition \ref{prop:non_separable} is given by
\begin{align}
G_{NS}(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y)= \,\, \pi^2 \sum_{k_1\in\mathcal{K}_1}\sum_{k_2\in\mathcal{K}_2}q_{k_1,k_2} \text{sinc}\left(\frac{\Omega_x - k_1\,\Omega_{0x}}{\Omega_{0x}}+\frac{\Omega_y - k_2\,\Omega_{0y}}{\Omega_{0y}}\right)
\text{sinc}\left(\frac{\Omega_y - k_2\,\Omega_{0y}}{\Omega_{0y}}\,-\, \frac{\Omega_x - k_1\,\Omega_{0x}}{\Omega_{0x}}\right),
\label{eq:G_omega_non_sepa}
\end{align}
where $ q_{k_1,k_2}$ are the partial fraction decomposition coefficients. The corresponding spatial-domain kernel is
\begin{align}
g_{NS}(x,y)=\frac{\Omega_{0x}\Omega_{0y}}{8}\,\,\text{rect}\left( \frac{\Omega_{0x}x + \Omega_{0y}y}{4\pi}\right)\text{rect}\left( \frac{\Omega_{0y}y - \Omega_{0x}x}{4\pi}\right)  \sum_{k_1\in\mathcal{K}_1}\sum_{k_2\in\mathcal{K}_2}q_{k_1,k_2}\, e^{\mathrm{j}(k_1\Omega_{0x}x + k_2\Omega_{0y}y)},
\label{eq:s_non_sepa}
\end{align}
where $\text{rect}(\cdot)$ is defined as  $\text{rect}(\frac{x}{T_x}) = 1$, if $|x| \leq T_x/2$, and $0$ otherwise. The sampling kernel derived in (\ref{eq:s_non_sepa}) is for a generic polynomial function $V(s_x,s_y)$, except that it does not have zeros on $\mathcal{S}$. For the particular choice \\
$V(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y) =\hspace{-1.5mm}\sum \limits_{\substack{p_1\in\mathcal{K}_1 \\ p_2\in\mathcal{K}_2}}\hspace{-1mm} d_{p_1,p_2}
\hspace{-1mm} \prod \limits_{\substack{k_1\in\mathcal{K}_1\backslash\{p_1\} \\ k_2\in\mathcal{K}_2\backslash\{p_2\}}}
\hspace{-1mm}\left(\frac{\mathrm{j}\,\Omega_x}{\Omega_{0x}} + \frac{\mathrm{j}\,\Omega_y}{\Omega_{0y}} - \mathrm{j}k_1 - \mathrm{j}k_2 \right) ~ \left(\frac{\mathrm{j}\,\Omega_y}{\Omega_{0y}} - \frac{\mathrm{j}\,\Omega_x}{\Omega_{0x}} - \mathrm{j}k_2 + \mathrm{j}k_1 \right)$,
we get $q_{k_1,k_2} = d_{k_1,k_2}$. \\
\subsubsection{Discussion}
Even though the nonseparable kernel $g_{NS}(x,y)$ in (\ref{eq:s_non_sepa}) appears to be a rotated version of the separable kernel $g_S(x,y)$ in (\ref{eq:2-D_SMS}) with $r_1=r_2=1$, a closer observation reveals something more. If we rotate the kernel $g_S(x,y)$ or equivalently $G_S(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y)$, then the zeros of the kernel will also shift in the 2-D plane. Hence, the alias cancellation conditions specified in (\ref{eq:alias_cancle_1}) and (\ref{eq:alias_cancle_2}) are no more satisfied on the 2-D rectangular grid, but are valid on the rotated 2-D grid. This means, to counter the effect of rotation, the sampling mechanism and the reconstruction techniques have to be suitably modified. On the other hand, in the case of the proposed nonseparable kernel $G_{NS}(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y)$ in (\ref{eq:G_omega_non_sepa}), the alias cancellation conditions are met on the 2-D grid, and the usual sampling and reconstruction techniques that are applicable in the case of separable kernels could be deployed.

Unlike the separable kernels, in the case of a nonseparable kernel $g_{NS}(x,y)$, we have considered only polynomial B-splines of zeroth order. Developing nonseparable kernels that are a sum of modulated-splines of higher orders needs more investigation and makes an interesting case for future work. Design of such kernels might result in their impulse responses being non-isotropic, and could be employed to approximate a more wider class of point-spread functions in the imaging modalities such as localization microscopy, radio astronomy, etc. Another aspect that is worthy for future investigation is the analysis and design of nonseparable  kernels that reproduce exponentials of the form $e^{\mathrm{j}\alpha xy}$ for some $\alpha$.

\section{Simulation Results}
\begin{figure}[t]
\begin{minipage}[b]{.48\linewidth}
  \centering
  \centerline{\includegraphics[width=0.65\linewidth,trim={2.4cm 0 1cm 0},clip]{time_domain_kernel_L_4.pdf}}

  \centerline{(a)}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.48\linewidth}
  \centering
  \centerline{\includegraphics[width=0.65\linewidth,trim={2.4cm 0 1cm 0},clip]{freq_domain_kernel_L_4.pdf}}

  \centerline{(b)}\medskip
\end{minipage}
\begin{minipage}[b]{0.48\linewidth}
  \centering
  \centerline{\includegraphics[width=0.65\linewidth,trim={2.4cm 0 1cm 0},clip]{convolution_op_L_4.pdf}}

  \centerline{(c)}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.48\linewidth}
  \centering
  \centerline{\includegraphics[width=0.65\linewidth,trim={2.4cm 0 1cm 0},clip]{Estimated_ground_truth_diracs.pdf}}

  \centerline{(d)}\medskip
\end{minipage}
\hfill
\caption{[Color online] (a) Impulse response; (b) Frequency response of the nonseparable sampling kernel $g_{NS}(x,y)$ in (\ref{eq:s_non_sepa}) with $\mathcal{K}_1 =\mathcal{K}_2 = \llbracket -L,L\rrbracket$ for $L=4$, and $\Omega_{0x} =\Omega_{0y} = \frac{\pi}{0.99}$; (c) Spatial-domain signal obtained by convolving the Diracs with the sampling kernel; and (d) The ground truth and reconstructed Dirac locations.}
\label{fig:dirac_sim}
\end{figure}

\begin{figure}[t]
\begin{minipage}[b]{.48\linewidth}
  \centering
  \centerline{\includegraphics[width=0.65\linewidth,trim={2.4cm 0 1cm 0},clip]{convolution_op_spline_L_3_Q_5_15_dB_SNR.pdf}}
 \vspace{-2mm}
  \centerline{(a)}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.48\linewidth}
  \centering
  \centerline{\includegraphics[width=0.65\linewidth,trim={1.4cm 0 1cm 0},clip]{Estimated_ground_truth_splines.pdf}}
  \vspace{-2mm}
  \centerline{(b)}\medskip
\end{minipage}
\hfill
\caption{[Color online] Localization of Gaussian blobs: (a) Result of convolution of three ($L=3$) truncated Gaussian functions with the nonseparable sampling kernel with $\mathcal{K}_1 = \mathcal{K}_2 = \llbracket -5L,5L\rrbracket$ and $\Omega_{0x} = \Omega_{0y} = \frac{\pi}{0.99}$; and  (b) Original 2-D FRI signal with the estimated locations (marked `+'). The SNR of the measurements used for reconstruction was $15$ dB and the MSE in the estimation of locations (over 50 realizations) was computed to be $-55$ dB.}
\label{fig:spline_sim}
\end{figure}

We conducted two experiments to demonstrate the applicability of the proposed non-separable sampling kernel. In the first experiment, we consider a 2-D signal having four Diracs $(L = 4)$ with parameters $\{\gamma_\ell,x_\ell,\, y_\ell\}_{\ell = 1}^{4}$ that are selected uniformly at random over $(0\,,\,1)$. A compactly supported nonseparable kernel as in (\ref{eq:s_non_sepa}) with the support restricted to $\left[-{T_{0x}}, {T_{0x}}\right] \times \left[-{T_{0y}}, {T_{0y}}\right]$ is simulated with the following parameters:  $\mathcal{K}_1 = \mathcal{K}_2 =\llbracket -4,4\rrbracket$,  $\Omega_{0x} = \Omega_{0y} =\frac{\pi}{0.99}$. The spatial domain samples are acquired at the critical sampling rate of  $\Omega_{sx} = (2L + 1)\Omega_{0x}$ and $\Omega_{sy} = (2L + 1)\Omega_{0y}$. The parameters were estimated using the algebraically coupled matrix pencil reconstruction method \cite{2DHarmonicRetrieval_vanpoucke1994}. The input signal along with the estimated Dirac locations, the impulse and frequency response of the sampling kernel, and convolution output of the signal and the sampling kernel are shown in Fig. \ref{fig:dirac_sim}. The mean-square error in the estimation of Dirac locations was computed to be $-280$ dB implying perfect reconstruction up to machine precision.

In the second experiment, we consider a 2-D FRI signal consisting of  three truncated Gaussian functions located at $\{x_\ell,\, y_\ell\}_{\ell = 1}^{3}$ that are selected uniformly at random over $(0\,,\,1)$. Zero-mean,  additive white Gaussian noise is added to the spatial-domain samples such that the signal-to-noise ratio (SNR) of the resulting samples is $15$ dB. We oversample the signal with the sampling rates $\Omega_{sx} = (10L + 1)\Omega_{0x}$ and $\Omega_{sy} = (10L + 1)\Omega_{0y}$. Figure \ref{fig:spline_sim} shows the ground truth signal with the accurately estimated locations of the Gaussian blobs. A comprehensive assessment of noise robustness of the kernels vis-\`a-vis the Cram\'er-Rao bounds will be addressed in a future work.

\section{Conclusions}
In this paper, we proposed a generalized framework for designing compactly supported sampling kernels for 2-D FRI signals. The first key idea in this generalization was to design the frequency response of the kernel, which satisfies a set of alias cancellation conditions; and second, to characterize admissible kernels with a compact spatial support by invoking the 2-D Paley-Wiener theorem.  The proposed framework allows for the design of both separable and nonseparable 2-D sampling kernels. As a particular case, we showed that a special case of the separable sampling kernel is the class of 2-D SMS kernels,  which has the attractive feature of reproducing a certain class of exponentials and  the support of the kernel is independent of the order. We also demonstrated the design of a nonseparable kernel
and validated it by performing simulations to extract the exact locations of Diracs in the 2-D plane. The design of higher-order nonseparable kernels and analysis of their exponential reproducing properties are some interesting aspects for further study.

\appendices
\section{Proof of Proposition \ref{prop:separable}}
\label{appendix:prop_separable}

The Proof of Proposition \ref{prop:separable} involves two steps: (i) showing $G_S(s_x,s_y)$ is an EFET; and (ii) showing $G_S(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y) \in L^2(\mathbb{R})$.
\subsection{Proof that $G_S(s_x,s_y)$ is an EFET}
\label{appendix:G_efet_separable}
\begin{proof}[Proof]

For the separable case, we have
\begin{equation}
    Z_{r_1,r_2}(s_x,s_y) = \sinh^{r_1}\left(\frac{\pi s_x}{\Omega_{0x}} \right) \sinh^{r_2}\left(\frac{\pi s_y}{\Omega_{0y}}\right). \nonumber
\end{equation}

We know that

\begin{align}
    \left|\sinh\left(\frac{\pi s_x}{\Omega_{0x}} \right) \sinh\left(\frac{\pi s_y}{\Omega_{0y}}\right)\right| \leq& \left|\sinh\left(\frac{\pi s_x}{\Omega_{0x}} \right)\right| \left|\sinh\left(\frac{\pi s_y}{\Omega_{0y}}\right)\right|, \nonumber \\
    =& \left| \frac{e^{\frac{\pi s_x}{\Omega_{0x}}} - e^{ - \frac{\pi s_x}{\Omega_{0x}}} }{2} \right| \left|\frac{e^{\frac{\pi s_y}{\Omega_{0y}}} - e^{ - \frac{\pi s_y}{\Omega_{0y}}} }{2} \right|, \nonumber \\
    \leq& \frac{1}{4} \left( \left| e^{\frac{\pi s_x}{\Omega_{0x}}}\right| + \left|e^{ - \frac{\pi s_x}{\Omega_{0x}}} \right| \right) \left(\left|e^{\frac{\pi s_y}{\Omega_{0y}}} \right| + \left| e^{-\frac{\pi s_y}{\Omega_{0y}}} \right|\right), \nonumber \\
    \leq& \frac{1}{4} \left(  e^{\left|\frac{\pi s_x}{\Omega_{0x}}\right|} + e^{\left|\frac{\pi s_x}{\Omega_{0x}}\right|}  \right) \left(e^{\left|\frac{\pi s_y}{\Omega_{0y}}\right|}  +  e^{\left|\frac{\pi s_y}{\Omega_{0y}}\right|} \right), \nonumber \\
    =& \frac{1}{4} \left( 2\, e^{\pi\left|\frac{s_x}{\Omega_{0x}}\right|}\right) \left( 2\, e^{\pi\left|\frac{ s_y}{\Omega_{0y}}\right|}\right), \nonumber \\
    =& \, e^{\left(\frac{\pi}{\Omega_{0x}}\left|s_x \right|  + \frac{\pi}{\Omega_{0y}}\left|s_y\right|\right)} . \nonumber

\end{align}

Consequently,
\begin{equation}
    \left|\sinh^{r_1}\left(\frac{\pi s_x}{\Omega_{0x}} \right) \sinh^{r_2}\left(\frac{\pi s_y}{\Omega_{0y}}\right) \right| \leq \,  e^{\left(\frac{r_1\pi}{\Omega_{0x}} \left|s_x \right|  + \frac{r_2\pi}{\Omega_{0y}}\left|s_y\right|\right)}. \nonumber
\end{equation}

Since $Z(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y)$  is a sum of entire functions, it is an EFET, and consequently $G_S(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y)$ is also an EFET (as shown in \cite{satish_paley_wiener}). The support of $g_S(x,y)$ is given by $\left[-\frac{r_1\pi}{\Omega_{0x}}, \frac{r_1\pi}{\Omega_{0x}}\right] \times \left[-\frac{r_2\pi}{\Omega_{0y}}, \frac{r_2\pi}{\Omega_{0y}}\right]$ i.e., $\Big[-\frac{r_1\,T_{0x}}{2}, \frac{r_1\,T_{0x}}{2}\Big] \times \left[-\frac{r_2\,T_{0y}}{2}, \frac{r_2\,T_{0y}}{2}\right]$.

\end{proof}

\subsection{Proof that $G_S(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y) \in L^2(\mathbb{R})$}
\label{appendix:G_energy_separable}
\begin{proof}[Proof]

Using the hyperbolic sine identity: $\sinh(\mathrm{j}\theta)=\mathrm{j}\sin(\theta)$, $\forall \,\theta\in \mathbb{R}$, $G_S(\mathrm{j}\Omega_{0x}, \mathrm{j}\Omega_{0y})$ can be expressed as
\begin{equation}
\centering
G_S(\mathrm{j}\Omega_{0x},\mathrm{j}\Omega_{0y}) = V(\mathrm{j}\Omega_{0x},\mathrm{j}\Omega_{0y})\,\mathrm{j}^{(r_1 (1 - |\mathcal{K}_1|) + r_2(1 - |\mathcal{K}_2|))}\frac{\sin^{r_1}\left(\frac{\pi \Omega_x}{\Omega_{0x}} \right) \sin^{r_2}\left(\frac{\pi \Omega_y}{\Omega_{0y}}\right)}{\prod \limits_{\substack{k_1\in\mathcal{K}_1 \\ k_2\in\mathcal{K}_2}} (\Omega_{x}-k_1\Omega_{0x})^{r_1}(\Omega_{y}-k_2\Omega_{0y})^{r_2}}. \nonumber
\end{equation}
Using partial fraction decomposition, $G_S(\mathrm{j}\Omega_{0x}, \mathrm{j}\Omega_{0y})$ takes the form
\begin{equation}
\centering
G_S(\mathrm{j}\Omega_{0x},\mathrm{j}\Omega_{0y})=\sum_{k_1\in\mathcal{K}_1}\sum_{k_2\in\mathcal{K}_2}\sum_{m_1=1}^{r_1}\sum_{m_2=1}^{r_2}\frac{q_{k_1,k_2,m_1,m_2}}{(\Omega_{x}-k_1\Omega_{0x})^{m_1}(\Omega_{y}-k_1\Omega_{0y})^{m_2}}\sin^{r_1}\left(\frac{\pi \Omega_x}{\Omega_{0x}} \right) \sin^{r_2}\left(\frac{\pi \Omega_y}{\Omega_{0y}}\right),
\label{gs_partial}
\end{equation}
where $q_{k_1,k_2,m_1,m_2}$ are the coefficients of the partial fraction expansion and are given by
\begin{equation}
\centering
q_{k_1,k_2,m_1,m_2}=\frac{\mathrm{j}^{(r_1 (1 - |\mathcal{K}_1|) + r_2(1 - |\mathcal{K}_2|))}}{(r_1-m_1)!(r_2-m_2)!}\frac{\text{d}^{(r_1-m_1)}}  {\text{d}\Omega_{x}^{(r_1-m_1)}} \frac{\text{d}^{(r_2-m_2)}}  {\text{d}\Omega_{y}^{(r_2-m_2)}} \left. \left( \frac{V(\mathrm{j}\Omega_{0x},\mathrm{j}\Omega_{0y})}{\prod \limits_{\substack{p_1\in\mathcal{K}_1\backslash\{k_1\} \\ p_2\in\mathcal{K}_2\backslash\{k_2\}}} (\Omega_{x}-p_1\Omega_{0x})^{r_1}(\Omega_{y}-p_1\Omega_{0y})^{r_2}} \right ) \right|_{\substack{\Omega_x=k_1\Omega_{0x} \\ \Omega_y=k_2\Omega_{0y}}},\nonumber
\end{equation}
for $m_1= 1, 2, \cdots, r_1$, and $m_2= 1, 2, \cdots, r_2$.

Using the Cauchy-Schwarts inequality, we have
\begin{align}
    \int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}|G_S(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y)|^2 \text{d}\Omega_x\text{d}\Omega_y = & \int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty} \left|\sum_{k_1\in\mathcal{K}_1}\sum_{k_2\in\mathcal{K}_2}\sum_{m_1=1}^{r_1}\sum_{m_2=1}^{r_2}q_{k_1,k_2,m_1,m_2}\frac{\sin^{r_1}\left(\frac{\pi \Omega_x}{\Omega_{0x}} \right) \sin^{r_2}\left(\frac{\pi \Omega_y}{\Omega_{0y}}\right)}{(\Omega_{x}-k_1\Omega_{0x})^{m_1}(\Omega_{y}-k_1\Omega_{0y})^{m_2}}\right|^2\text{d}\Omega_x\text{d}\Omega_y, \nonumber \\
    \leq& \,\,r_1 r_2 |\mathcal{K}_1||\mathcal{K}_2|\sum_{k_1\in\mathcal{K}_1}\sum_{k_2\in\mathcal{K}_2}\sum_{m_1=1}^{r_1}\sum_{m_2=1}^{r_2}\,|q_{k_1,k_2,m_1,m_2}|^2 \nonumber \\
    & \int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}\left|\frac{\sin^{r_1}\left(\frac{\pi \Omega_x}{\Omega_{0x}} \right) \sin^{r_2}\left(\frac{\pi \Omega_y}{\Omega_{0y}}\right)}{(\Omega_{x}-k_1\Omega_{0x})^{m_1}(\Omega_{y}-k_1\Omega_{0y})^{m_2}}\right|^2 \text{d}\Omega_x\text{d}\Omega_y, \nonumber \\
    = & \,\,r_1 r_2 |\mathcal{K}_1||\mathcal{K}_2|\sum_{k_1\in\mathcal{K}_1}\sum_{k_2\in\mathcal{K}_2}\sum_{m_1=1}^{r_1}\sum_{m_2=1}^{r_2}\,|q_{k_1,k_2,m_1,m_2}|^2  \int\limits_{-\infty}^{\infty}\left|\frac{\sin^{r_1}\left(\frac{\pi \Omega_x}{\Omega_{0x}} \right) }{(\Omega_{x}-k_1\Omega_{0x})^{m_1}}\right|^2 \text{d}\Omega_x \nonumber \\
    & \int\limits_{-\infty}^{\infty}\left|\frac{\sin^{r_2}\left(\frac{\pi \Omega_y}{\Omega_{0y}}\right)}{(\Omega_{y}-k_1\Omega_{0y})^{m_2}}\right|^2\text{d}\Omega_y. \nonumber
\end{align}

As shown in \cite{satish_paley_wiener}, $\int\limits_{-\infty}^{\infty}\left|\frac{\sin^{r_1}\left(\frac{\pi \Omega_x}{\Omega_{0x}} \right) }{(\Omega_{x}-k_1\Omega_{0x})^{m_1}}\right|^2 \text{d}\Omega_x < \infty$ and $\int\limits_{-\infty}^{\infty}\left|\frac{\sin^{r_2}\left(\frac{\pi \Omega_y}{\Omega_{0y}}\right)}{(\Omega_{y}-k_1\Omega_{0x})^{m_2}}\right|^2\text{d}\Omega_y < \infty$. Hence,  we have \\ $\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}|G_S(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y)|^2 \text{d}\Omega_x\text{d}\Omega_y < \infty$.

\end{proof}

\subsection{Expression for $g_S(x,y)$}
\label{appendix:kernel_separable_expression}

Consider
\begin{equation}
\centering
G_S(\mathrm{j}\Omega_{0x},\mathrm{j}\Omega_{0y})=\sum_{k_1\in\mathcal{K}_1}\sum_{k_2\in\mathcal{K}_2}\sum_{m_1=1}^{r_1}\sum_{m_2=1}^{r_2}\frac{q_{k_1,k_2,m_1,m_2}}{(\Omega_{x}-k_1\Omega_{0x})^{m_1}(\Omega_{y}-k_1\Omega_{0y})^{m_2}}\sin^{r_1}\left(\frac{\pi \Omega_x}{\Omega_{0x}} \right) \sin^{r_2}\left(\frac{\pi \Omega_y}{\Omega_{0y}}\right).
\end{equation}
Further simplifying the above equation, we get
\begin{align}
G_S(\mathrm{j}\Omega_{0x},\mathrm{j}\Omega_{0y})=&\sum_{k_1\in\mathcal{K}_1}\sum_{k_2\in\mathcal{K}_2}\sum_{m_1=1}^{r_1}\sum_{m_2=1}^{r_2} q_{k_1,k_2,m_1,m_2} (\Omega_x-k\Omega_{0x})^{r_1-m_1}(\Omega_y-k_2\Omega_{0y})^{r_2-m_2}\,\, \frac{\sin^{r_1}\left( \frac{\pi \Omega_x}{\Omega_{0x}} \right)}{(\Omega_x-k_1\Omega_{0x})^{r_1}}\frac{\sin^{r_2}\left( \frac{\pi \Omega_y}{\Omega_{0y}} \right)}{(\Omega_y-k\Omega_{0y})^{r_2}}, \nonumber \\
=&\sum_{k_1\in\mathcal{K}_1}\sum_{k_2\in\mathcal{K}_2}\sum_{m_1=1}^{r_1}\sum_{m_2=1}^{r_2} q_{k_1,k_2,m_1,m_2}(\Omega_x-k\Omega_{0x})^{r_1-m_1}(\Omega_y-k_2\Omega_{0y})^{r_2-m_2}(-1)^{r_1k_1 + r_2k_2} \nonumber \\
&\frac{\sin^{r_1}\left( \frac{ \Omega_x-k_1\Omega_{0x}}{\Omega_{0x}}\pi \right)}{(\Omega_x-k_1\Omega_{0x})^{r_1}}\frac{\sin^{r_2}\left( \frac{ \Omega_y-k_2\Omega_{0y}}{\Omega_{0y}}\pi \right)}{(\Omega_y-k_2\Omega_{0y})^{r_2}}, \nonumber\\
=& \left(\frac{\pi^{(r_1 + r2)}}{\Omega_{0x}^{r_1}\Omega_{0y}^{r_2}}\right)\sum_{k_1\in\mathcal{K}_1}\sum_{k_2\in\mathcal{K}_2}\sum_{m_1=1}^{r_1}\sum_{m_2=1}^{r_2} q_{k_1,k_2,m_1,m_2} (-1)^{r_1k_1 + r_2k_2} \sum_{c_1=0}^{r_1-m_1} {r_1-m_1 \choose c_1} (-k_1\Omega_{0x})^{r_1-m_1-c_1}\Omega_x^{c_1} \nonumber \\  & \text{sinc}^{r_1}\left( \frac{ \Omega_x-k_1\Omega_{0x}}{\Omega_{0x}} \right) \sum_{c_2=0}^{r_2-m_2} {r_2-m_2 \choose c_2} (-k_2\Omega_{0y})^{r_2-m_2-c_2}\Omega_y^{c_2} \text{sinc}^{r_2}\left( \frac{ \Omega_y-k_2\Omega_{0y}}{\Omega_{0y}} \right) .
\label{eq:G1}
\end{align}
Now, applying the inverse Fourier transform on (\ref{eq:G1}), $g_S(x,y)$ is obtained as
\begin{align}
g_S(x,y)=& \left(\frac{\pi^{(r_1 + r2)}}{\Omega_{0x}^{r_1}\Omega_{0y}^{r_2}}\right) \sum_{k_1\in\mathcal{K}_1}\sum_{k_2\in\mathcal{K}_2}\sum_{m_1=1}^{r_1}\sum_{m_2=1}^{r_2} q_{k_1,k_2,m_1,m_2} (-1)^{r_1k_1 + r_2k_2} \sum_{c_1=0}^{r_1-m_1} {r_1-m_1 \choose c_1} (-k_1\Omega_{0x})^{r_1-m_1-c_1} \nonumber\\
&{\frac{1}{\mathrm{j}^{c_1}} \frac{\text{d}^{c_1}}{\text{d}x^{c_1}}\left[e^{\mathrm{j}k_1\Omega_{0x}x} \frac{1}{T_{0x}}\beta^{r_1-1} \left(\frac{x}{T_{0x}}\right) \right]} \sum_{c_2=0}^{r_2-m_2} {r_2-m_2 \choose c_2} (-k_2\Omega_{0y})^{r_2-m_2-c_2} {\frac{1}{\mathrm{j}^{c_2}} \frac{\text{d}^{c_2}}{\text{d}y^{c_2}}\left[e^{\mathrm{j}k_2\Omega_{0y}y} \frac{1}{T_{0y}}\beta^{r_2-1} \left(\frac{y}{T_{0y}}\right) \right]},
\label{eq:g_x_y_sepa}
\end{align}
where $\beta^{r}(t)$ denotes the $r^{\text{th}}$-order polynomial B-spline.

Now, evaluating the above equation for $V(\mathrm{j}\Omega_{0x},\mathrm{j}\Omega_{0y})$ given in (\ref{eq:separable_V}), we arrive at the class of separable 2-D SMS kernels given in (\ref{eq:2-D_SMS}).

\section{Proof of Proposition \ref{prop:non_separable}}
\label{appendix:prop_nonseparable}

The Proof of Proposition \ref{prop:non_separable} involves two steps: (i) showing $G_{NS}(s_x,s_y)$ is an EFET; and (ii) showing $G_{NS}(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y) \in L^2(\mathbb{R})$.
\subsection{Proof that $G_{NS}(s_x,s_y)$ is an EFET }
\label{appendix:G_efet_nonseparable}
\begin{proof}[Proof]

For the nonseparable case, we have
\begin{equation}
    Z_{r_1,r_2}(s_x,s_y) = \sinh\left(\frac{\pi s_x}{\Omega_{0x}} + \frac{\pi s_y}{\Omega_{0y}}\right)\sinh\left(\frac{\pi s_y}{\Omega_{0y}} - \frac{\pi s_x}{\Omega_{0x}}\right). \nonumber
\end{equation}

We know that
\begin{align}
    \left|\sinh\left(\frac{\pi s_x}{\Omega_{0x}} + \frac{\pi s_y}{\Omega_{0y}}\right)\sinh\left(\frac{\pi s_y}{\Omega_{0y}} - \frac{\pi s_x}{\Omega_{0x}}\right)\right| \leq& \left|\sinh\left(\frac{\pi s_x}{\Omega_{0x}} + \frac{\pi s_y}{\Omega_{0y}}\right)\right| \left|\sinh\left(\frac{\pi s_y}{\Omega_{0y}} - \frac{\pi s_x}{\Omega_{0x}}\right)\right|, \nonumber \\
    =&\, \left| \frac{e^{\left(\frac{\pi s_x}{\Omega_{0x}} + \frac{\pi s_y}{\Omega_{0y}}\right)} - e^{ - \left(\frac{\pi s_x}{\Omega_{0x}} + \frac{\pi s_y}{\Omega_{0y}}\right)} }{2} \right| \left|\frac{e^{\left(\frac{\pi s_y}{\Omega_{0y}} - \frac{\pi s_x}{\Omega_{0x}}\right)} - e^{ - \left(\frac{\pi s_y}{\Omega_{0y}} - \frac{\pi s_x}{\Omega_{0x}}\right)} }{2} \right|, \nonumber \\
    \leq&\, \frac{1}{4} \left( \left| e^{\left(\frac{\pi s_x}{\Omega_{0x}} + \frac{\pi s_y}{\Omega_{0y}}\right)}\right| + \left|e^{ - \left(\frac{\pi s_x}{\Omega_{0x}} + \frac{\pi s_y}{\Omega_{0y}}\right)} \right| \right) \left(\left|e^{\left(\frac{\pi s_y}{\Omega_{0y}} - \frac{\pi s_x}{\Omega_{0x}}\right)} \right| + \left| e^{-\left(\frac{\pi s_y}{\Omega_{0y}} - \frac{\pi s_x}{\Omega_{0x}}\right)} \right|\right), \nonumber \\
    \leq& \,\frac{1}{4} \left(  e^{\left|\frac{\pi s_x}{\Omega_{0x}} + \frac{\pi s_y}{\Omega_{0y}}\right|} + e^{\left|\frac{\pi s_x}{\Omega_{0x}} + \frac{\pi s_y}{\Omega_{0y}}\right|}  \right) \left(e^{\left|\frac{\pi s_y}{\Omega_{0y}} - \frac{\pi s_x}{\Omega_{0x}}\right|}  +  e^{\left|\frac{\pi s_y}{\Omega_{0y}} - \frac{\pi s_x}{\Omega_{0x}}\right|} \right), \nonumber \\
    =& \,\frac{1}{4} \left( 2\, e^{\pi\left|\frac{ s_x}{\Omega_{0x}} + \frac{ s_y}{\Omega_{0y}}\right|}\right) \left( 2\, e^{\pi\left|\frac{ s_y}{\Omega_{0y}} - \frac{ s_x}{\Omega_{0x}}\right|}\right), \nonumber \\
    \le& \, \left(e^{\pi \left(\left|\frac{s_x}{\Omega_{0x}}\right| + \left|\frac{s_y}{\Omega_{0y}}\right|\right)} \right)\left( e^{\pi \left(\left|\frac{s_y}{\Omega_{0y}}\right| + \left|\frac{s_x}{\Omega_{0x}}\right|\right)} \right), \nonumber \\
    =& \, e^{\left(\frac{2\pi}{\Omega_{0x}}\left|s_x \right|  + \frac{2\pi}{\Omega_{0y}}\left|s_y\right|\right)} . \nonumber

\end{align}

Since $Z(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y)$  is a sum of entire functions, it is an EFET, and consequently $G_{NS}(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y)$ is an EFET as well. The support of $g_{NS}(x,y)$ is given by $\left[-\frac{2\pi}{\Omega_{0x}}, \frac{2\pi}{\Omega_{0x}} \right] \times \left[-\frac{2\pi}{\Omega_{0y}}, \frac{2\pi}{\Omega_{0y}} \right] $ i.e., $\left[-T_{0x}, T_{0x} \right] \times \left[-T_{0y}, T_{0y} \right] $.

\end{proof}

\subsection{Proof that $G_{NS}(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y) \in L^2(\mathbb{R})$}
\label{appendix:G_energy_nonseparable}
\begin{proof}[Proof]

Using the hyperbolic sine identity: $\sinh(\mathrm{j}\theta)=\mathrm{j}\sin(\theta)$, $\forall \,\theta\in \mathbb{R}$, $G_{NS}(\mathrm{j}\Omega_{0x}, \mathrm{j}\Omega_{0y})$ is expressed as
\begin{align}
G_{NS}(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y)
= \frac{\mathrm{j}^{(2 - |\mathcal{K}_1| - |\mathcal{K}_2|)}\,V(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y)\,\sin\left(\frac{\pi\,\Omega_x}{\Omega_{0x}}\,+\, \frac{\pi\,\Omega_y}{\Omega_{0y}}\right)\sin\left(\frac{\pi\,\Omega_y}{\Omega_{0y}}\,-\, \frac{\pi\,\Omega_x}{\Omega_{0x}}\right)}{\prod \limits_{\substack{k_1\in\mathcal{K}_1 \\ k_2\in\mathcal{K}_2}}\left(\frac{\Omega_x}{\Omega_{0x}} + \frac{\Omega_y}{\Omega_{0y}} - k_1 - k_2 \right)\left(\frac{\Omega_y}{\Omega_{0y}} - \frac{\Omega_x}{\Omega_{0x}} - k_2 + k_1 \right)}, \nonumber
\end{align}

which when expanded using partial fraction decomposition takes the form
\begin{equation}
\centering
G_{NS}(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y)=\sum_{k_1\in\mathcal{K}_1}\sum_{k_2\in\mathcal{K}_2}\frac{q_{k_1,k_2}}{\left(\frac{\Omega_x}{\Omega_{0x}} + \frac{\Omega_y}{\Omega_{0y}} - k_1 - k_2 \right)\left(\frac{\Omega_y}{\Omega_{0y}} - \frac{\Omega_x}{\Omega_{0x}} - k_2 + k_1 \right)}\sin\left(\frac{\pi\,\Omega_x}{\Omega_{0x}}\,+\, \frac{\pi\,\Omega_y}{\Omega_{0y}}\right)\sin\left(\frac{\pi\,\Omega_y}{\Omega_{0y}}\,-\, \frac{\pi\,\Omega_x}{\Omega_{0x}}\right).
\label{eq:non_sepa_partial}
\end{equation}
$q_{k_1,k_2}$ are the coefficients of the partial fraction expansion, and are given by
\begin{equation}
\centering
q_{k_1,k_2}= \left.\left( \frac{\mathrm{j}^{(2 - |\mathcal{K}_1| - |\mathcal{K}_2|)}\,V(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y)}{\prod \limits_{\substack{p_1\in\mathcal{K}_1\backslash\{k_1\} \\ p_2\in\mathcal{K}_2\backslash\{k_2\}}}\left(\frac{\Omega_x}{\Omega_{0x}} + \frac{\Omega_y}{\Omega_{0y}} - p_1 - p_2 \right)\left(\frac{\Omega_y}{\Omega_{0y}} - \frac{\Omega_x}{\Omega_{0x}} - p_2 + p_1 \right)} \right ) \right|_{\substack{\Omega_x=k_1\Omega_{0x} \\ \Omega_y=k_2\Omega_{0y}}}.\nonumber
\end{equation}

Using the Cauchy-Schwarts inequality, we have
\begin{align}
    \int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}|G_{NS}(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y)|^2 \text{d}\Omega_x\text{d}\Omega_y = & \int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty} \left| \sum_{k_1\in\mathcal{K}_1}\sum_{k_2\in\mathcal{K}_2}\,q_{k_1,k_2}\,\frac{\sin\left(\frac{\pi\,\Omega_x}{\Omega_{0x}}\,+\, \frac{\pi\,\Omega_y}{\Omega_{0y}}\right)\sin\left(\frac{\pi\,\Omega_y}{\Omega_{0y}}\,-\, \frac{\pi\,\Omega_x}{\Omega_{0x}}\right)}{\left(\frac{\Omega_x - k_1\,\Omega_{0x}}{\Omega_{0x}} + \frac{\Omega_y - k_2\,\Omega_{0y}}{\Omega_{0y}}\right)\left(\frac{\Omega_y - k_2\,\Omega_{0y}}{\Omega_{0y}} - \frac{\Omega_x - k_1\,\Omega_{0x}}{\Omega_{0x}}\right)}\right|^2\text{d}\Omega_x\text{d}\Omega_y, \nonumber \\
    \leq& |\mathcal{K}_1||\mathcal{K}_2|\sum_{k_1\in\mathcal{K}_1}\sum_{k_2\in\mathcal{K}_2}\,|q_{k_1,k_2}|^2 \nonumber \\ & \int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}\left|\frac{\sin\left(\frac{\pi\,\Omega_x}{\Omega_{0x}}\,+\, \frac{\pi\,\Omega_y}{\Omega_{0y}}\right)\sin\left(\frac{\pi\,\Omega_y}{\Omega_{0y}}\,-\, \frac{\pi\,\Omega_x}{\Omega_{0x}}\right)}{\left(\frac{\Omega_x - k_1\,\Omega_{0x}}{\Omega_{0x}} + \frac{\Omega_y - k_2\,\Omega_{0y}}{\Omega_{0y}}\right)\left(\frac{\Omega_y - k_2\,\Omega_{0y}}{\Omega_{0y}} - \frac{\Omega_x - k_1\,\Omega_{0x}}{\Omega_{0x}}\right)}\right|^2 \text{d}\Omega_x\text{d}\Omega_y . \nonumber
\end{align}

Let us consider
\begin{align}
    \zeta = \int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}\left|\frac{\sin\left(\frac{\pi\,\Omega_x}{\Omega_{0x}}\,+\, \frac{\pi\,\Omega_y}{\Omega_{0y}}\right)\sin\left(\frac{\pi\,\Omega_y}{\Omega_{0y}}\,-\, \frac{\pi\,\Omega_x}{\Omega_{0x}}\right)}{\left(\frac{\Omega_x - k_1\,\Omega_{0x}}{\Omega_{0x}} + \frac{\Omega_y - k_2\,\Omega_{0y}}{\Omega_{0y}}\right)\left(\frac{\Omega_y - k_2\,\Omega_{0y}}{\Omega_{0y}} - \frac{\Omega_x - k_1\,\Omega_{0x}}{\Omega_{0x}}\right)}\right|^2 \text{d}\Omega_x\text{d}\Omega_y .
    \label{eq:double_int}
\end{align}

Substituting $\widetilde{\Omega}_x = \frac{\Omega_x - k_1\,\Omega_{0x}}{\Omega_{0x}}$ and $ \widetilde{\Omega}_y = \frac{\Omega_y - k_2\,\Omega_{0y}}{\Omega_{0y}}  $ in (\ref{eq:double_int}) and simplifying it further, we get

\begin{align}
    \zeta =& \,\,\Omega_{0x}\Omega_{0y} \int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}\left|\frac{\sin\left(\pi\,\widetilde{\Omega}_x + k_1\pi\,+\, \pi\,\widetilde{\Omega}_y + k_2\pi \right)\sin\left(\pi\,\widetilde{\Omega}_y + k_2\pi\,-\, \pi\,\widetilde{\Omega}_x + k_1\pi\right)}{\left(\widetilde{\Omega}_x + \widetilde{\Omega}_y\right)\left(\widetilde{\Omega}_x + \widetilde{\Omega}_y\right)}\right|^2 \text{d}\widetilde{\Omega}_x\text{d}\widetilde{\Omega}_y, \nonumber \\
    =& \,\,\Omega_{0x}\Omega_{0y} \int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}\left|\frac{\sin\left(\pi\,(\widetilde{\Omega}_x + \widetilde{\Omega}_y) \right)\sin\left(\pi\,(\widetilde{\Omega}_y - \widetilde{\Omega}_x)\right)}{\left(\widetilde{\Omega}_x + \widetilde{\Omega}_y\right)\left(\widetilde{\Omega}_x + \widetilde{\Omega}_y\right)}\right|^2 \text{d}\widetilde{\Omega}_x\text{d}\widetilde{\Omega}_y, \nonumber \\
    =& \,\,\pi^4\Omega_{0x}\Omega_{0y} \int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}\left|\text{sinc}\left(\widetilde{\Omega}_x + \widetilde{\Omega}_y \right)\text{sinc}\left(\widetilde{\Omega}_y - \widetilde{\Omega}_x\right)\right|^2 \text{d}\widetilde{\Omega}_x\text{d}\widetilde{\Omega}_y. \label{eq:before_rot}
\end{align}

Consider rotation of a function $\vartheta(\widetilde{\Omega}_x,\widetilde{\Omega}_y) = \text{sinc}\left(\widetilde{\Omega}_x + \widetilde{\Omega}_y \right)\text{sinc}\left(\widetilde{\Omega}_y - \widetilde{\Omega}_x\right)$  by $\theta = 45^{\circ}$. Since energy of $\vartheta(\widetilde{\Omega}_x,\widetilde{\Omega}_y)$ is preserved by the rotation operator $R = \left[\begin{smallmatrix}\cos\theta&-\sin\theta\\\sin\theta&\cos\theta\end{smallmatrix}\right]$, we substitute the rotated $\vartheta(\widetilde{\Omega}_x,\widetilde{\Omega}_y)$ in \eqref{eq:before_rot}. Thus the right-hand side of (\ref{eq:before_rot}) is given by
\begin{align}
    \zeta =&\,\,\pi^4\Omega_{0x}\Omega_{0y} \int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}\left|\text{sinc}\left(\sqrt{2}\widetilde{\Omega}_x \right)\text{sinc}\left(\sqrt{2}\widetilde{\Omega}_y\right)\right|^2 \text{d}\widetilde{\Omega}_x\text{d}\widetilde{\Omega}_y, \nonumber \\
    =& \,\,\pi^4\Omega_{0x}\Omega_{0y} \int\limits_{-\infty}^{\infty}\left|\text{sinc}\left(\sqrt{2}\widetilde{\Omega}_x \right)\right|^2\text{d}\widetilde{\Omega}_x\int\limits_{-\infty}^{\infty}\left| \text{sinc}\left(\sqrt{2}\widetilde{\Omega}_y\right)\right|^2 \text{d}\widetilde{\Omega}_y,\nonumber\\
    =& \,\,\pi^4\Omega_{0x}\Omega_{0y}\left(\frac{1}{\sqrt{2}}\right)\left( \frac{1}{\sqrt{2}}\right),\nonumber\\
    =& \,\,\frac{\pi^4\Omega_{0x}\Omega_{0y}}{2}.\nonumber
\end{align}
Hence, we have $\zeta < \infty $, and consequently, $\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}|G_{NS}(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y)|^2 \text{d}\Omega_x\text{d}\Omega_y < \infty $.

\end{proof}

\subsection{Expression for $g_{NS}(x,y)$}
\label{appendix:kernel_nonseparable_expression}
Consider
\begin{equation}
\centering
G_{NS}(\mathrm{j}\Omega_x, \mathrm{j}\Omega_y)=\sum_{k_1\in\mathcal{K}_1}\sum_{k_2\in\mathcal{K}_2}\frac{q_{k_1,k_2}}{\left(\frac{\Omega_x}{\Omega_{0x}} + \frac{\Omega_y}{\Omega_{0y}} - k_1 - k_2 \right)\left(\frac{\Omega_y}{\Omega_{0y}} - \frac{\Omega_x}{\Omega_{0x}} - k_2 + k_1 \right)}\sin\left(\frac{\pi\,\Omega_x}{\Omega_{0x}}\,+\, \frac{\pi\,\Omega_y}{\Omega_{0y}}\right)\sin\left(\frac{\pi\,\Omega_y}{\Omega_{0y}}\,-\, \frac{\pi\,\Omega_x}{\Omega_{0x}}\right).
\end{equation}
Upon simplifying further, we get
\begin{align}
G_{NS}(\mathrm{j}\Omega_{0x},\mathrm{j}\Omega_{0y}) =& \sum_{k_1\in\mathcal{K}_1}\sum_{k_2\in\mathcal{K}_2}\,q_{k_1,k_2}\,\frac{\sin\left(\frac{\pi\,\Omega_x}{\Omega_{0x}}\,+\, \frac{\pi\,\Omega_y}{\Omega_{0y}}\right)\sin\left(\frac{\pi\,\Omega_y}{\Omega_{0y}}\,-\, \frac{\pi\,\Omega_x}{\Omega_{0x}}\right)}{\left(\frac{\Omega_x - k_1\,\Omega_{0x}}{\Omega_{0x}} + \frac{\Omega_y - k_2\,\Omega_{0y}}{\Omega_{0y}}\right)\left(\frac{\Omega_y - k_2\,\Omega_{0y}}{\Omega_{0y}} - \frac{\Omega_x - k_1\,\Omega_{0x}}{\Omega_{0x}}\right)}, \nonumber\\
=&  \sum_{k_1\in\mathcal{K}_1}\sum_{k_2\in\mathcal{K}_2}\,q_{k_1,k_2}\,(-1)^{(k_1 + k_2 + k_2 - k_1)}\frac{\sin\left(\frac{\pi\,\Omega_x}{\Omega_{0x}}\,+\, \frac{\pi\,\Omega_y}{\Omega_{0y}}  - \pi(k_1 + k_2)\right)\sin\left(\frac{\pi\,\Omega_y}{\Omega_{0y}}\,-\, \frac{\pi\,\Omega_x}{\Omega_{0x}}  - \pi(k_2 - k_1)\right)}{\left(\frac{\Omega_x - k_1\,\Omega_{0x}}{\Omega_{0x}} + \frac{\Omega_y - k_2\,\Omega_{0y}}{\Omega_{0y}}\right)\left(\frac{\Omega_y - k_2\,\Omega_{0y}}{\Omega_{0y}} - \frac{\Omega_x - k_1\,\Omega_{0x}}{\Omega_{0x}}\right)}, \nonumber\\
=&  \,\, \pi^2 \sum_{k_1\in\mathcal{K}_1}\sum_{k_2\in\mathcal{K}_2}\,q_{k_1,k_2}\,{\text{sinc}\left(\frac{\Omega_x - k_1\,\Omega_{0x}}{\Omega_{0x}}\,+\, \frac{\Omega_y - k_2\,\Omega_{0y}}{\Omega_{0y}}\right)\text{sinc}\left(\frac{\Omega_y - k_2\,\Omega_{0y}}{\Omega_{0y}}\,-\, \frac{\Omega_x - k_1\,\Omega_{0x}}{\Omega_{0x}}\right)}.
\label{eq:G_non_sepa}
\end{align}
Using the Lemma \ref{lemma:inverse_fourier_nonseparable_sinc} in Appendix \ref{appendix:add_proposition} and the frequency shift property, the inverse Fourier transform of (\ref{eq:G_non_sepa}) is given by
\begin{align}
g_{NS}(x,y)=&\,\,\,\,\frac{\Omega_{0x}\Omega_{0y}}{8}\,\,\text{rect}\left( \frac{\Omega_{0x}x + \Omega_{0y}y}{4\pi}\right)\text{rect}\left( \frac{\Omega_{0y}y - \Omega_{0x}x}{4\pi}\right)\, \sum_{k_1\in\mathcal{K}_1}\sum_{k_2\in\mathcal{K}_2}q_{k_1,k_2}\, e^{\mathrm{j}(k_1\Omega_{0x}x + k_2\Omega_{0y}y)}.
\label{eq:g_x_y_non_sepa}
\end{align}

\section{ }
\label{appendix:add_proposition}

\begin{lemma}
\label{lemma:inverse_fourier_nonseparable_sinc}
If $F(\Omega_x, \Omega_y) = \text{sinc}\left(\frac{\Omega_x}{\Omega_{0x}} + \frac{\Omega_y}{\Omega_{0y}}\right)\text{sinc}\left(\frac{\Omega_y}{\Omega_{0y}} - \frac{\Omega_x}{\Omega_{0x}}\right)$, then its inverse Fourier transform is given by $f(x,y) = \,\,\frac{\Omega_{0x}\Omega_{0y}}{8\pi^2}\,\,\text{rect}\left( \frac{\Omega_{0x}x + \Omega_{0y}y}{4\pi}\right)\text{rect}\left( \frac{\Omega_{0y}y - \Omega_{0x}x}{4\pi}\right)$.
\end{lemma}
\begin{proof}[Proof]

We know that the 1-D inverse Fourier transform of $F(\Omega) = \text{sinc}(\sqrt{2}\Omega)$  is given by

\begin{equation}
    f(t) = \frac{1}{2\sqrt{2}\pi}\text{rect}\left(\frac{t}{2\sqrt{2}\pi}\right).\nonumber
\end{equation}

For the separable two-dimensional counterpart, the inverse Fourier transform of the function
\begin{equation}
    \widetilde{F}(\Omega_x, \Omega_y) = \text{sinc}(\sqrt{2}\Omega_x)\text{sinc}(\sqrt{2}\Omega_y), \label{eq:foo}
\end{equation}
 is given by,
\begin{equation}
    \widetilde{f}(x,y) = \frac{1}{8\pi^2}\text{rect}\left(\frac{x}{2\sqrt{2}\pi}\right)\text{rect}\left(\frac{y}{2\sqrt{2}\pi}\right). \nonumber
\end{equation}

Applying rotation of $\theta = - 45^{\circ}$ on the function in (\ref{eq:foo}) and using the 2-D rotation property of the Fourier transform we have,
\begin{align}
F\left(\frac{\Omega_x}{\sqrt{2}} + \frac{\Omega_y}{\sqrt{2}}, - \frac{\Omega_x}{\sqrt{2}} + \frac{\Omega_y}{\sqrt{2}}\right)
=&\,\,\text{sinc}\left({\Omega_x} + {\Omega_y}\right)\text{sinc}\left({\Omega_y} - {\Omega_x}\right),\nonumber
\end{align}
and its corresponding inverse Fourier transform is given by
\begin{align}
\hspace{-.6in} f(x,y) =& \,\,\frac{1}{8\pi^2}\text{rect}\left(\frac{\frac{x}{\sqrt{2}} + \frac{y}{\sqrt{2}}}{2\sqrt{2}\pi}\right)\text{rect}\left(\frac{\frac{y}{\sqrt{2}} - \frac{x}{\sqrt{2}}}{2\sqrt{2}\pi}\right),\nonumber\\
\hspace{0.1in}=& \,\,\frac{1}{8\pi^2}\text{rect}\left(\frac{x + y}{4\pi}\right)\text{rect}\left(\frac{y - x}{4\pi}\right).
\label{eq:goo}
\end{align}

Using the two dimensional scaling property on (\ref{eq:goo}), the inverse Fourier transform of

\begin{equation}
    F(\Omega_x, \Omega_y) = \text{sinc}\left(\frac{\Omega_x}{\Omega_{0x}} + \frac{\Omega_y}{\Omega_{0y}}\right)\text{sinc}\left(\frac{\Omega_y}{\Omega_{0y}} - \frac{\Omega_x}{\Omega_{0x}}\right),\nonumber
\end{equation}
is given by,
\begin{equation}
    f(x,y) = \,\,\frac{\Omega_{0x}\Omega_{0y}}{8\pi^2}\,\,\text{rect}\left( \frac{\Omega_{0x}x + \Omega_{0y}y}{4\pi}\right)\text{rect}\left( \frac{\Omega_{0y}y - \Omega_{0x}x}{4\pi}\right).\nonumber
\end{equation}

\end{proof}

\bibliographystyle{IEEEbib}


\section{Introduction}

The auditory and vestibular systems display remarkable mechanical sensitivity.  The end organs are able to detect motion in the \AA\space regime, below the level of thermal fluctuations in the surrounding fluid. \cite{Hudspeth14}  Humans are able to resolve two clicks temporally separated by only 10 microseconds, \cite{Leshowitz71} with the stimulus waveforms presented into both ears.  The sensitivity and temporal resolution of the inner ear are not fully understood, and physics of its signal detection remains an active area of research. \cite{Reichenbach14}

Mechanical detection of sound and acceleration is performed by hair cells, sensory cell that were named after the organelle that protrudes from their apical surface.  This organelle consists of rod-like stereovilli that are organized in interconnecting rows and are collectively named the hair bundle. Transduction channels are embedded at the tips of these stereovilli and are coupled to the interconnecting tip links.  When an external signal deflects the hair bundle, it modulates the opening probability of its transduction channels. Thus, the hair bundle transduces the mechanical energy of the signal into an electrical potential change in the cell via the influx of ionic current through the channels. \cite{LeMasurier05}

Hair bundles of several species exhibit autonomous oscillations in the absence of an external stimulus. \cite{Martin03}  These spontaneous oscillations were shown to violate the fluctuation dissipation theorem, \cite{Martin01} indicating that they are a manifestation of an internal active mechanism. \cite{Hudspeth08}  We previously demonstrated in a series of experiments that these spontaneous oscillations exhibit chaotic dynamics. \cite{Faber18a}  By tuning the experimental parameters of the biological preparation, we were able to modulate the degree of chaos in the system and demonstrated that the weakly chaotic regime gives rise to enhanced sensitivity. \cite{Faber18b}.  Further, we showed that the temporal resolution of the detector increases with increasing levels of chaos.  These results were shown to be consistent with simulations of a numerical model for hair cell dynamics.

In the current work, we provide a theoretical treatment of the beneficial role of chaos in signal detection by the hair cells. We focus on noise-induced chaos, as biological sensory systems operate in thermal environments. The normal form equation for the Hopf bifurcation has  been used to model many dynamical systems. \cite{Marsden76} It was shown to capture a number of experimentally observed phenomena exhibited by the inner ear, including the sensitivity and frequency selectivity of hearing.  \cite{Eguiluz00, Kern03}  Furthermore, it provides us with the simplest model that reproduces the main characteristics of hair bundle dynamics. We consider here the nonisochronous Hopf oscillator, for which the angular frequency depends on the radius. \cite{Pikovsky01}  By varying the parameters of the system, we modulate the degree of nonisochronicity, and thus the level of chaos that would be observed in the presence of noise. \cite{Faber18b}  We calculate the sensitivity to a constant-force stimulus analytically and find the nonisochronous system to be both more sensitive and faster to respond than the isochronous one.

Experiments performed on live hair cells have shown that the bundles exhibit a non-monotonic response to large, step-function stimuli, which resembles the ringing of an underdamped oscillator. \cite{Benser96, Tinevez07}  This so-called ``twitch'' is observed only in living, active cells and is believed to be a manifestation of an internal active process, as a passive hair bundle exhibits an overdamped response. The phenomenon may play an important role in hearing, as it provides a possible mechanism by which the hair cell can amplify an incoming signal.  In the current work, we show that the Hopf oscillator can reproduce the twitch when the stimulus induces a saddle-node on invariant circle (SNIC) bifurcation resulting in a spiral-sink fixed point.  Further, we find the parameter conditions that allow such behavior, providing us with a very simple model for this biological phenomenon.

Lastly, we explore the routes to chaos in the system, by finding the exact analytic expressions for the stroboscopic maps in the presence of a periodic delta-function stimulus.  Using these maps, we find a period-doubling cascade to chaos upon increasing the forcing strength when the system is nonisochronous.  Further, we find the same route to chaos upon increasing the degree of nonisochronicity, while keeping the forcing strength constant.

\section{Theoretical Model}

We use the normal-form equation for the Hopf bifurcation, including terms up to cubic order, and introduce additive forcing:

\begin{ceqn}
\begin{align}
\frac{dz}{dt} = (\mu + i\omega_0)z - (\alpha + i\beta)|z|^2z + F_z(t),
\label{eq:Hopf}
\end{align}
\end{ceqn}

\noindent where $z(t) = x(t) + iy(t)$ and $F_z(t) = F_x(t) + iF_y(t)$.  Here, $x(t)$ represents the bundle position, while $y(t)$ reflects internal parameters of the bundle and is not assigned a specific measurable quantity.  Chaotic dynamics arise when this system is driven by common Gaussian white noise: \cite{Zhou02}  $F_z(t) = \eta_x(t) + i\eta_y(t)$, where $\langle \eta_x(t)\eta_x(t') \rangle  = \langle \eta_y(t)\eta_y(t') \rangle = 2D \delta (t-t')$, and $\langle \eta_x(t)\eta_y(t') \rangle  = 0$ (Fig. \ref{fig1}).  We previously demonstrated that, under these conditions, the Lyapunov exponent can be approximated as \cite{Faber18b}

\begin{ceqn}
\begin{align}
\lambda \approx \frac{D\beta}{\mu}.
\label{eq:Lyapunov}
\end{align}
\end{ceqn}

\noindent When $\beta \neq 0$, the oscillator is nonisochronous: the frequency is dependent on the amplitude of oscillations.  Prior work has demonstrated that nonisochronous systems exhibit noise-induced chaos \cite{Faber18b}.  For simplicity, we study this instability in a noiseless system, while noting that in the presence of noise, the degree of chaos would be proportional to $\beta$.

Auditory and vestibular stimuli induce lateral deflections on the hair bundle, so we consider forces in the $\hat{x}$ direction, which coincides with the direction of autonomous oscillation.  In polar coordinates, the model takes the form:

\begin{ceqn}
\begin{align}
\frac{dr}{dt} = \mu r - \alpha r^3 + F_x(t)\cos\theta
\label{eq:r_dot}
\end{align}
\end{ceqn}

\begin{ceqn}
\begin{align}
\frac{d\theta}{dt} = \omega_0 - \beta r^2 - F_x(t)\frac{\sin\theta}{r}.
\label{eq:theta_dot}
\end{align}
\end{ceqn}

\noindent We observe that, in the absence of forcing, there is a stable limit cycle of radius $r_0 = \sqrt{\frac{\mu}{\alpha}}$ for $\mu > 0$.  Further, there is a nullcline in the phase at radius $r_{nc} = \sqrt{\frac{\omega_0}{\beta}}$.  If the force causes the stable limit cycle to approach this radius, the angular frequency slows down, and the system approaches a SNIC bifurcation.  We show that the forcing can perturb the shape of the limit cycle, causing an intersection between the stable limit cycle and $r_{nc}$.  This SNIC bifurcation causes a pair of fixed points (one stable, one unstable) to appear on the limit cycle (see Appendix A).  When the system is poised in close proximity to this bifurcation, noise can most easily induce a chaotic response, and the system is very sensitive to external signals.

\section{Response to a Constant-Force Stimulus}

We consider forcing of the form $F_x(t) = f_0 \Theta(t-t_0)$, where $\Theta$ is the Heaviside step function, and look at the response of a system in the oscillatory regime ($\mu > 0$).  We assume that the forcing is weak compared to the stability of the limit cycle ($\frac{f_0}{\mu r_0} << 1$).  We can therefore assume that the forcing acts as a small perturbation on the shape of the limit cycle and ignore the transient in $r$ at the onset of the forcing.  Looking for solutions near the limit cycle, we let $r(t) = r_0 + \delta r(t)$, where $r_0 = \sqrt{\frac{\mu}{\alpha}}$ is the limit cycle radius in the absence of forcing.  Inserting this approximation into Eq. (\ref{eq:r_dot}) and keeping only terms that are linear in $\delta r$ yields

\begin{ceqn}
\begin{align}
\frac{d}{dt}\delta r = -2\mu\delta r + f_0 \Theta(t-t_0)\cos\theta,
\end{align}
\end{ceqn}

\noindent which has a steady-state solution $\delta r = \frac{f_0\cos\theta}{2\mu}$ for $t > t_0$.  Therefore, the first order perturbation to the shape of the limit cycle is given by:

\begin{ceqn}
\begin{align}
r(t) = r_0 + \frac{f_0\cos\theta(t)}{2\mu}.
\label{eq:r}
\end{align}
\end{ceqn}

\noindent Inserting this solution into Eq. (\ref{eq:theta_dot}) and keeping only the first-order forcing terms, we find that

\begin{ceqn}
\begin{align}
\frac{d\theta}{dt} = \Omega_0 - \frac{f_0\sin\theta}{r_0} - \frac{\beta r_0 f_0\cos\theta}{\mu},
\end{align}
\end{ceqn}

\noindent where $\Omega_0 = \omega_0 - \beta r_0^2$ is the natural frequency.  We integrate the equation to solve for $\theta (t)$:

\begin{ceqn}
\begin{align}
\int_{t_0}^{t} dt' = \int_{\theta_0}^{\theta} \frac{d\theta'}{\Omega_0 - \frac{f_0\sin\theta'}{r_0} - \frac{\beta r_0 f_0\cos\theta'}{\mu}},
\end{align}
\end{ceqn}

\noindent where $\theta_0$ is the phase of the oscillator at the onset of the step. Evaluating the integral yields

\begin{ceqn}
\begin{align}
-\frac{1}{2} \gamma (t - t_0)  =  \tanh^{-1}(Q(\theta)) - \tanh^{-1}(Q(\theta_0)),
\label{eq:tanh}
\end{align}
\end{ceqn}

\noindent where

\begin{ceqn}
\begin{align}
\gamma = \sqrt{ \Big(\frac{f_0}{r_0}\Big)^2 + \Big(\frac{\beta r_0 f_0}{\mu}\Big)^2 - \Omega_0^2}
\label{eq:gamma}
\end{align}
\end{ceqn}

\noindent and

\begin{ceqn}
\begin{align}
Q(\theta) = \frac{(\Omega_0 + \frac{\beta r_0 f_0}{\mu})\tan(\frac{\theta}{2}) - \frac{f_0}{r_0}}{\gamma}.
\end{align}
\end{ceqn}

\noindent Inverting Eq. (\ref{eq:tanh}) yields

\begin{ceqn}
\begin{align}
\theta(t) = 2 \tan^{-1} \Bigg( \frac{ \gamma \tanh \Big(-\frac{1}{2} \gamma (t - t_0) + \tanh^{-1}(Q(\theta_0))\Big) + \frac{f_0}{r_0} }{ \Omega_0 + \frac{\beta r_0 f_0}{\mu} } \Bigg).
\label{eq:theta}
\end{align}
\end{ceqn}

\noindent Using the identity $\tanh^{-1}(x) = \frac{1}{2}\log(1 + x) - \frac{1}{2}\log(1 - x)$, we can also write Eq. (\ref{eq:tanh}) in exponental form:

\begin{ceqn}
\begin{align}
e^{- \gamma (t - t_0)}  =  \frac{(1 + Q(\theta))(1 - Q(\theta_0))}{(1 - Q(\theta))(1 + Q(\theta_0))}.
\label{eq:exp}
\end{align}
\end{ceqn}

\begin{figure*}[t!]
\includegraphics[width=17cm]{Fig2}
\caption{\label{fig2} (a) Mean displacement of the Hopf oscillator in response to a constant-force stimulus found semi-analytically (black points) and fully analytically (solid curve).  (b) Response time of the oscillator characterized by the exponential decay time in response to a step force.  Black points represent the values found from numerical simulations, and the solid curve represents the analytic approximation.  For (a-b),  $\mu = \alpha = \Omega_0 = 1$ and $f_0 = 0.01$.  The dashed line indicates the point where the forcing is able to induce a SNIC bifurcation.  (c) An example of the twitch, a non-monotonic response (solid) to a step force (dashed), for $\mu = \alpha = \Omega_0 = 1$, $\beta = 50$, and $f_0 = 0.1$.  The response curve represents the average over 200 different initial conditions, evenly separated in phase.}
\end{figure*}

\subsection{Response Amplitude}

We use the mean displacement in $x$ to characterize the response amplitude:

\begin{ceqn}
\begin{align}
\langle \Delta x(t) \rangle = \langle r(t)\cos\theta(t) \rangle_{stim} - \langle r(t)\cos\theta(t)\rangle_{0},
\end{align}
\end{ceqn}

\noindent where the first and second terms on the right side are time averages in the presence and absence of the stimulus, respectively.  The second term yields zero.  Inserting Eq. (\ref{eq:r}), we obtain

\begin{ceqn}
\begin{align}
\langle \Delta x(t) \rangle = r_0 \langle \cos\theta(t) \rangle_{stim} + \frac{f_0}{2\mu} \langle \cos^{2}\theta(t) \rangle_{stim},
\label{eq:mean_disp}
\end{align}
\end{ceqn}

\noindent which can be solved semi-analytically by inserting Eq. (\ref{eq:theta}) and taking the time average over a long time series (Fig. \ref{fig2}).

Depending on the choice of parameters, $\gamma$ is either purely real or purely imaginary. Eq. (\ref{eq:exp}) indicates that the solutions are oscillatory for imaginary $\gamma$ and decay exponentially to a fixed point for real $\gamma$.  In the latter case, we can find the mean displacement analytically.  This fixed point corresponds to the intersection between the perturbed limit cycle, $r(t) = r_0 + \frac{f_0\cos\theta}{2\mu}$, and the angular nullcline $r_{nc} = \sqrt{\frac{\omega_0}{\beta}}$.  We set these two radii equal to each other and find that

\begin{ceqn}
\begin{align}
\cos\theta = \frac{2\mu}{f_0} \bigg(\sqrt{\frac{\omega_0}{\beta}} - r_0 \bigg).
\end{align}
\end{ceqn}

\noindent We insert this solution into Eq. (\ref{eq:mean_disp}) and substitute in the natural frequency, $\Omega_0 = \omega_0 - \beta r_0^2$, to obtain:

\begin{ceqn}
\begin{align}
\langle \Delta x(t) \rangle = \frac{2\mu r_0}{f_0} \bigg(\sqrt{\frac{\Omega_0}{\beta} + r_0^2} - r_0 \bigg) +  \frac{2\mu}{f_0} \bigg(\sqrt{\frac{\Omega_0}{\beta} + r_0^2} - r_0 \bigg)^2.
\label{mean_disp_analytic}
\end{align}
\end{ceqn}

\noindent This analytic curve is plotted in Fig. \ref{fig2}.  Note that at a fixed natural frequency $\Omega_0$ and in the limit of large $\beta$, the radius of the angular nullcline $r_{nc}$ approaches the limit cycle radius.  With these two circles close together, any small perturbation to the limit cycle can cause an intersection, resulting in a stable fixed point.

\subsection{Temporal Resolution}

We characterize the temporal resolution of the system by calculating the response time (decay time) to a step stimulus.  Our results are limited to real $\gamma$, as this regime shows an exponential decay to a stable fixed point upon the application of a stimulus.  For imaginary $\gamma$, the system would continue to oscillate in the presence of a step stimulus.  We perform a Taylor series expansion of the right side of Eq. (\ref{eq:exp}) in the vicinity of the stable fixed point (see Appendix B):

\begin{ceqn}
\begin{align}
\theta(t) \approx \frac{1-Q(\theta_c)^2}{2Q'(\theta_c)}\Big(e^{- \gamma (t - t_0)} - 1\Big) + \theta_0,
\end{align}
\end{ceqn}

\noindent where $\theta_c$ is the angle of the stable fixed point.  The decay time of this solution is

\begin{ceqn}
\begin{align}
\tau_{res} = \frac{1}{\gamma} = \frac{1}{ \sqrt{ (\frac{f_0}{r_0})^2 + (\frac{\beta r_0 f_0}{\mu})^2 - \Omega_0^2} }.
\end{align}
\end{ceqn}

\noindent We fit exponential functions to the numerical simulations of the response and compare the values to the analytic form of $\tau_{res}$ in Fig. \ref{fig2}.

\begin{figure*}[t!]
\includegraphics[width=17cm]{Fig3}
\caption{\label{fig3} Mean displacement of the Hopf oscillator in response to a brief, constant-force stimulus.  (a) and (b) represent the responses to a pulse of duration $\frac{1}{4}$ and $\frac{1}{16}$ of the natural period of oscillation, respectively.  $\mu = \alpha = \Omega_0 = 1$, $f_0 = 0.01$, and the vertical dashed line represents the value of $\beta$ at which the forcing is able to induce a SNIC bifurcation.}
\end{figure*}

\subsection{Non-monotonic response}

In the regime of large $\beta$ and/or strong forcing, the averaged response curve does not resemble an exponential decay.  Rather, the response is non-monotonic and resembles the observed phenomenon known in hair bundle mechanics as ``the twitch'' \cite{Benser96, Tinevez07} (Fig. \ref{fig2}). In our simple model, the non-monotonic behavior is a result of a spiral sink forming after the step force induces a SNIC bifurcation.  The trajectories spiraling to the fixed point give rise to the ringing behavior in $x(t)$.  In this section, we find the conditions necessary for the Hopf oscillator to produce this phenomenon.  We find the fixed point $(R_{c}, \theta_{c})$ using Eqs. (\ref{eq:r_dot}) and (\ref{eq:theta_dot}):

\begin{ceqn}
\begin{align}
0 = \mu R_{c} - \alpha R_{c}^3 + f_0\cos\theta_{c}
\label{eq:r_dot_fp}
\end{align}
\end{ceqn}

\begin{ceqn}
\begin{align}
0 = \omega_0 - \beta R_{c}^2 - f_0\frac{\sin\theta_{c}}{R_{c}}.
\label{eq:theta_dot_fp}
\end{align}
\end{ceqn}

\noindent The Jacobian at the fixed point takes the form:

\begin{ceqn}
\begin{align}
J =
\quad
\begin{pmatrix}
\mu - 3\alpha R_{c}^2 & -f_0\sin\theta_{c} \\
- 2\beta R_{c} + f_0\frac{\sin\theta_{c}}{R_{c}^2} & - f_0\frac{\cos\theta_{c}}{R_{c}}
\end{pmatrix}
\quad
\end{align}
\end{ceqn}

\noindent which can also be written as

\begin{ceqn}
\begin{align}
J =
\quad
\begin{pmatrix}
\mu - 3\alpha R_{c}^2 & \beta R_c^3 - \omega_0 R_c \\
 \frac{\omega_0}{R_c} - 3\beta R_c & \mu - \alpha R_c^2
\end{pmatrix}
\quad
\end{align}
\end{ceqn}

\noindent by using Eqs. (\ref{eq:r_dot_fp}) and (\ref{eq:theta_dot_fp}) to eliminate $\theta_c$.  A spiral sink is present if and only if the trace of J is less than zero, the determinant of J is greater than zero, and $4\det(J) - Tr(J)^2 > 0$. \cite{Strogatz94}  Since we are considering the limit cycle regime, we have $\mu > 0$.  We will also use the condition that $R_c \geq r_0 = \sqrt{\frac{\mu}{\alpha}}$.  Since this effect is most pronounced when $\beta$ is large, we will assume that the second term in Eq. (\ref{eq:theta_dot}) is larger in magnitude than the first ($\beta r_0^2 > \omega_0$).  With these assumptions, the first two conditions are always satisfied:

\begin{ceqn}
\begin{align}
Tr(J) = 2\mu - 4\alpha R_c^2 \leq 2\mu - 4\alpha r_0^2 = -2\mu < 0
\end{align}
\end{ceqn}

\noindent and

\begin{eqnarray}
\begin{split}
\det(J) = \mu^2 \bigg(3\Big(\frac{R_c}{r_0}\Big)^4 - 4\Big(\frac{R_c}{r_0}\Big)^2 + 1 \bigg)   \\
+ \omega_0^2 \bigg( 3\Big(\frac{\beta R_c^2}{\omega_0}\Big)^2 - 4\Big(\frac{\beta R_c^2}{\omega_0}\Big) + 1 \bigg)  >  0,
\end{split}
\end{eqnarray}

\noindent since $\frac{R_c}{r_0} \geq 1$ and $\frac{\beta R_c^2}{\omega_0} \geq \frac{\beta r_0^2}{\omega_0} > 1$.  The third and final condition is satisfied when

\begin{ceqn}
\begin{align}
\beta > \frac{2\omega_0 + \sqrt{3\mu^2 + \omega_0^2}}{3r_0^2}
\end{align}
\end{ceqn}

\noindent  and $\beta > \frac{\alpha}{\sqrt{3}}$ (see Appendix C).  Thus, we have found the conditions on the parameters of the Hopf oscillator such that it can reproduce the experimentally observed twitch.

\section{Response to a Short-Pulse Stimulus}

We now find the response amplitude to a brief, weak force of amplitude $f_0$ and duration $T$:

\begin{ceqn}
\begin{align}
F_x(t) = f_0 \Big(\Theta(t) - \Theta(t - T)\Big).
\end{align}
\end{ceqn}

\noindent The solution during the force step was already found in Eqs. (\ref{eq:r}) and (\ref{eq:theta}), provided that the initial conditions at the onset of the pulse are on the stable limit cycle.  We use these equations to find the displacement induced in $x(t)$.  The mean displacement is obtained by computing the average of $(x(T) - x(0))$ over many different initial conditions.  The mean displacement is used to quantify the responsiveness of the system.

In Fig. \ref{fig3}, we show that there is a local maximum in responsiveness at the SNIC bifurcation, as in the steady-state case.  However, we observe an additional local maximum as $\beta$ increases beyond the bifurcation point.  Increasing $\beta$ in this regime increases the speed at which the trajectories move across the circle to the fixed point.  This increase is beneficial for the responsiveness of the system.  However, the increase in $\beta$  also moves the fixed point in the $-\hat{x}$ direction, thereby reducing the mean displacement in $x$ (see Appendix A).  The competition between these two effects results in a local maximum.  The location of this maximum is dependent on the pulse duration, which determines the time needed for the trajectories to reach the fixed point, starting from the initial conditions.

\begin{figure*}[t!]
\includegraphics[width=17cm]{Fig4}
\caption{\label{fig4}  Bifurcation diagrams generated from the radial stroboscopic map at each impulse of the spike train.  The interval between spikes was set to be the natural period of oscillations ($\tau = \frac{2\pi}{\Omega_0}$, $\mu = \alpha = \Omega_0 = 1$).  (a) $\beta = 3$, while $f_0$ is varied.  (b) $f_0 = 0.1$, while $\beta$ is varied.}
\end{figure*}

\section{Effects of a Spike Train}

We now find the response to periodic, delta-function forcing.  We construct the stroboscopic Poincar\'e maps analytically \cite{Gonzalez83} and observe the transition to chaos.  First, we find the analytic solution in the absence of forcing by integrating Eqs. (\ref{eq:r_dot}) and (\ref{eq:theta_dot}) with $F_x(t) = 0$.  The first can be solved using partial fraction decomposition:

\begin{ceqn}
\begin{align}
\int_{t_i}^t dt' = \int_{r_i}^r \frac{dr'}{\mu r' - \alpha r'^3}  = \int_{r_i}^r  \Big( \frac{1}{\mu r'} - \frac{\alpha r'}{\mu (\alpha r'^2 - \mu)} \Big)dr',
\end{align}
\end{ceqn}

\noindent which yields

\begin{ceqn}
\begin{align}
r_{so}(t, r_i) = \sqrt{ \frac{\mu}{ \alpha - (\alpha - \frac{\mu}{r_i^2})e^{-2\mu t} } },
\end{align}
\end{ceqn}

\noindent where we have set the initial time $t_i$ to zero, and $r_i$ is the radius at this time.  We insert this solution into Eq. (\ref{eq:theta_dot}):

\begin{ceqn}
\begin{align}
\frac{d\theta}{dt} = \omega_0 - \frac{\beta\mu}{ \alpha - (\alpha - \frac{\mu}{r_i^2})e^{-2\mu t} }.
\end{align}
\end{ceqn}

\noindent Integrating this equation yields

\begin{ceqn}
\begin{align}
\theta_{so}(t, r_i, \theta_i) = \theta_i + \omega_0 t + \frac{\beta}{2\alpha} \log \Bigg( \frac{\mu}{\alpha r_i^2 (e^{2\mu t} - 1) + \mu}  \Bigg),
\end{align}
\end{ceqn}

\noindent where $\theta_i$ is the phase at $t=0$.  We now include the spike-train stimulus:

\begin{ceqn}
\begin{align}
F_x(t) = f_0 \sum_{n=1}^{n=\infty} \delta(t - n\tau).
\end{align}
\end{ceqn}

\noindent The first impulse occurs at time $t=\tau$.  Before this time ($t = \tau - \epsilon$), the solution is simply $r(\tau - \epsilon) = r_{so}(\tau - \epsilon, r_i)$ and $\theta(\tau - \epsilon) = \theta_{so}(\tau - \epsilon, r_i, \theta_i)$.  Projecting $F_x(\tau)$ onto polar coordinates, we find the solution at the time of the first impulse to be:

\begin{ceqn}
\begin{align}
r(t = \tau) = R_1 = r_{so}(\tau, r_i) + f_0 \cos\big( \theta_{so}(\tau, r_i, \theta_i)\big)
\end{align}
\end{ceqn}

\begin{ceqn}
\begin{align}
\theta(t=\tau) = \Theta_1 = \theta_{so}(\tau, r_i, \theta_i)  -  \frac{f_0 \sin\big( \theta_{so}(\tau, r_i, \theta_i)\big)}{ r_{so}(\tau, r_i) }.
\end{align}
\end{ceqn}

\noindent The solution between the first and second impulses ($t = \tau + \epsilon$) takes the form:

\begin{ceqn}
\begin{align}
r(\tau + \epsilon) = r_{so}(\epsilon, R_1)
\end{align}
\end{ceqn}

\begin{ceqn}
\begin{align}
\theta(\tau + \epsilon) = \theta_{so}(\epsilon, R_1, \Theta_1),
\end{align}
\end{ceqn}

\noindent where we treat the coordinates at the previous impulse as the initial conditions and reset the time.  We now construct the stroboscopic Poincar\'e maps using this recursive method.  We find the radius and phase at every spike, $R_n = r(t=n\tau)$ and $\Theta_n = \theta(t=n\tau)$, to be:

\begin{ceqn}
\begin{align}
R_{n+1} = r_{so}(\tau, R_n) + f_0 \cos\big( \theta_{so}(\tau, R_n, \Theta_n)\big)
\end{align}
\end{ceqn}

\begin{ceqn}
\begin{align}
\Theta_{n+1} = \theta_{so}(\tau, R_n, \Theta_n)  -  \frac{f_0 \sin\big( \theta_{so}(\tau, R_n, \Theta_n)\big)}{ r_{so}(\tau, R_n) }.
\end{align}
\end{ceqn}

\noindent These stroboscopic maps are valid for both the oscillatory and quiescent regimes.  We have made no assumptions on the choice of parameters or the strength of the forcing.  Upon increasing either $f_0$ or $\beta$, we observe a period-doubling cascade to chaos (Fig. \ref{fig4}).

\section{Discussion}

In a prior study, we demonstrated that nonisochronicity results in noise-induced chaos for even very weak levels of noise and any nonzero $\beta$.  In the present work, we demonstrated analytically that nonisochronicity is responsible for enhanced responsiveness to weak force steps, with the most responsive system being poised near the SNIC bifurcation. Further, we demonstrated that detection of a single constant-force pulse is highly rapid in a nonisochronous system. By varying $\beta$ in the Hopf oscillator equation, we control the degree of nonisochronicity and study its impact on response characteristics. We show that the oscillator can detect very short pulses for a wide range of $\beta$, with the ideal value depending on the duration of the pulse.  We also showed that the speed of the response, as measured by the exponential decay time to the constant force, increases with increasing $\beta$.

Although we considered the noiseless system in this work, upon including even a small amount of noise, $\beta$ can be directly mapped onto the Lyapunov exponent. We previously demonstrated enhancement of sensitivity and temporal resolution with numerical models, where we varied the Lyapunov exponent directly \cite{Faber18b}.  The results of this numerical study are consistent with that of the present work. Nonisochronicity in a noiseless system hence directly correlates with the degree of chaos in the presence of noise.

Numerical models of hair cell dynamics typically require many variables to reproduce the experimentally observed twitch.  Here, we demonstrated that this 2-dimensional model can capture the main features of this effect.  This has, thus far, been the simplest model of hair cell dynamics able to produce this behavior.  We have shown that the twitch can be explained through bifurcation theory: the stimulus yields a SNIC bifurcation with one of the fixed points being a spiral sink.  Further, we found the parameter conditions required to see the effect.

Finally, we found the exact analytic expressions for the stroboscopic Poincar\'e maps in $r$ and $\theta$, when the oscillator is subjected to periodic, impulsive forcing.  These maps are valid for any parameter regime and any degree of forcing.  Using these stroboscopic maps, we observe a period-doubling cascade to chaos upon increasing either the forcing strength or $\beta$.  The nature of this stimulus resembles that of efferent activity.  In hair cell dynamics, efferent activity is believed to be the biological feedback responsible for modulating the sensitivity of the detector and may serve as a protective mechanism in preventing damage caused by acoustic trauma. \cite{Patuzzi91}  In the present work, we demonstrated that an efferent-like stimulus can induce or modulate the degree of chaos, thus affecting the sensitivity of the system.  Future work entails measuring the change in chaoticity and sensitivity caused by efferent activity.

Chaos is often considered an unfavorable element in dynamical systems, as it limits their predictability and regularity.  However, chaos is likely present in a number of real-world systems, as they contain many degrees of freedom and exhibit nonlinearities. Here, we have shown analytically that the instabilities that give rise to chaotic dynamics are also responsible for enhanced sensitivity and temporal resolution in the Hopf oscillator, suggesting a beneficial role of chaos in the auditory and vestibular systems.  We propose that chaos may play a role in signal detection in other noisy biological systems where timing and sensitivity are essential.

\begin{acknowledgments}
The authors gratefully acknowledge support of NSF Physics of Living Systems, under grant 1705139.
\end{acknowledgments}