\section{Introduction}
\label{sec:intro}
In image processing, spatial differentiation accomplishes image sharpening and edge-based segmentation, with broad applications ranging from microscopy and medical imaging to industrial inspection and object detection.\cite{Gonzales2008, Markham1963, Abramoff2004, Brosnan2004,Dalal2005} In these applications, of particular importance are isotropic derivative operators, whose response is rotationally invariant. The simplest and most widely used isotropic derivative operator is the two-dimensional Laplacian.\cite{Rosenfeld1982}

Spatial differentiation can certainly be carried out with conventional digital electronic computation. However, there are many big-data applications that require real-time and high-throughput image processing, for which digital computations become challenging.\cite{Pham2000,Holyer1989} Optical analog computing may overcome this challenge by offering high-throughput derivative operation with almost no energy consumption other than the inherent optical loss associated with the differential operators.\cite{Solli2015,Gorlitz1975,Eu1973,Reinhardt1978,Sirohi1977,Kasprzak1982,Yao1971,Nesteruk2001,Warde1983,Davis1991,Lancis1997} Moreover, the recent developments of nanophotonic structures such as meta-surfaces and plasmonic structures have offered the possibility of doing optical differentiation using compact devices.
However, most of the existing works on optical spatial differentiation using nanophotonic structures are restricted to one dimension, whereas  most images are two-dimensional objects.\cite{Zhu2017, Golovastikov2015, Golovastikov2014, Silva2014, Youssefi2016, AbdollahRamezani2015, Pors2015} To date, there has been no optical realization of the two-dimensional Laplacian at transmission mode with a compact nanophotonic device.
\begin{figure}
\centering
\includegraphics[width=0.9\columnwidth]{Fig1.png}
\caption{\label{fig1}(a) Geometry of the photonic crystal slab differentiator, which consists of a photonic crystal slab separated from a uniform dielectric slab by an air gap. For $\epsilon=12$, the geometry parameters are: $d=0.55a$, $r=0.111a$, $d_s=0.07a$, $d_g=0.21a$. The plane above and below shows the input and output image respectively. The transmitted image (e.g. ring) through the device is the Laplacian of the incident image (e.g. disk) illuminated by a normally incident light with frequency $\omega_0 = 0.47656 \times 2\pi c/a$. (b) Coordinate system. (c) The Brillioun zone of the system.}
\end{figure}

In this article, we show that a two-dimensional Lapacian can be implemented with the use of a photonic crystal slab, one of the most widely studied nanophotonic structures.\cite{Chow2000,Joannopoulos2011,Zhou2014} As shown in Figure~\ref{fig1}, the structure consists of a photonic crystal slab separated from a uniform dielectric slab by an air gap. The key in this design is to achieve an unusual band structure that is isotropic for the two polarizations.
\section{Theoretical Analysis}
\label{sec:analysis}
Our objective is to realize a Lapacian for a two-dimensional input field. Specifically, for a normally incident light beam along the $z$-axis with a transverse field profile $S_{in}(x,y)$, we would like to design an optical device for which the transmitted beam has a profile $S_{out}(x,y) \propto \nabla^2 S_{in}(x,y)$, whre $\nabla^2 = \partial_x^2 + \partial_y^2$ is the Laplacian. The task of realizing the Laplacian $\nabla^2$ in the real space is equivalent to designing an optical system with response function
\begin{equation}
t(\bm{k}) \propto (k_x^2+k_y^2) \label{eq:t-propto-k}
\end{equation} in the wavevector space ($\bm{k}$-space).\cite{Bracewell1986}
Equation~(\ref{eq:t-propto-k}) requires that $t(\bm{k})=0$ at $|\bm{k}|=0$. In a lossless photonic crystal slab, the transmission for normally incident light can vanish near a guided resonance.\cite{Fan2002} Thus we consider in more details the transmission coefficient of a photonic crystal slab near normal incidence. Consider a single photonic band of guided resonances, as characterized by $\bm{k}$-dependent resonant frequencies $\omega(\bm{k})$ and radiative linewidths $\gamma(\bm{k})$. (Here $\bm{k} = (k_x, k_y)$ refers to the in-plane wavevector.) Near the resonant frequencies,
 the transmitted amplitude $t$ is expressed as \cite{Fan2002}
\begin{equation}
t(\omega,\bm{k}) = t_d + f \frac{\gamma(\bm{k})}{i[\omega-\omega(\bm{k})]+\gamma(\bm{k})}~, \label{eq:t_origin}
\end{equation}
where $\omega$ is the incident light frequency, $t_d$ is the direct transmission coefficient, and $f$ is related to the complex decaying amplitude of the resonance to the transmission side of the slab.

In general, $f$ is constrained by the direct process due to energy conservation and time-reversal symmetry.\cite{Fan2002,Wang2013} In particular, if $t_d =1$, that is, the direct
pathway has a $100\%$ transmission coefficient, then
\begin{equation}
f = -t_d = -1 \label{eq:f_td}
\end{equation}
even for structure without z-mirror symmetry such as shown in Figure \ref{fig1}, as has been derived in Ref. \cite{Wang2013}. In this special case,
\begin{equation}
t(\omega,\bm{k}) = 1 - \frac{\gamma(\bm{k})}{i[\omega-\omega(\bm{k})]+\gamma(\bm{k})}~. \label{eq:t_special}
\end{equation}
We denote $\omega_0 = \omega(\bm{k}=\bm{0})$, $\gamma_0 = \gamma(\bm{k}=\bm{0})$. Using Equation~(\ref{eq:t_special}), zero transmission occurs at the $\Gamma$ point when the incident wave has the frequency:
\begin{equation}
	\omega = \omega_0 ~.
\end{equation}
At $\omega = \omega_0$, we perform an expansion of the transmission coefficient $t$ near $\bm{k}=\bm{0}$ :
\begin{equation}
	t(\omega_0, \bm{k}) = 0 + \frac{\partial t}{\partial \omega(\bm{k})}\Bigr|_{\Gamma} \delta \omega(\bm{k}) + \frac{\partial t}{\partial \gamma(\bm{k})}\Bigr|_{\Gamma} \delta \gamma(\bm{k}) ~, \label{eq:t}
\end{equation}
where
\begin{equation}
\delta \omega(\bm{k}) = \omega(\bm{k}) - \omega_0,\qquad \delta \gamma(\bm{k}) = \gamma(\bm{k})-\gamma_0 ~,
\end{equation}
\begin{equation}
\frac{\partial t}{\partial \omega(\bm{k})}\Bigr|_{\Gamma} = -\frac{i}{\gamma_0},\qquad \frac{\partial t}{\partial \gamma(\bm{k})}\Bigr|_{\Gamma} = 0~,
\end{equation}
therefore,
\begin{equation}
t(\omega_0, \bm{k}) = - \frac{i}{\gamma_0}\delta \omega(\bm{k})~.
\end{equation}
In this special case, $t(\bm{k})$ is simply proportional to the band dispersion $\delta\omega(\bm{k})$ near $\bm{k}=\bm{0}$. If $\delta\omega(\bm{k}) \propto |\bm{k}|^2$, then $t(\bm{k})\propto |\bm{k}|^2$ as well. Therefore the transmission through a photonic crystal slab can be used to achieve the Laplacian.
\begin{figure*}
\centering
\includegraphics[width=1.5\columnwidth]{Fig2.png}
\caption{\label{fig2} Band structure of the photonic crystal slab with the dielectric constant  $\epsilon=12$, the thickness $d=0.55a$, and the radius $r=0.111a$ of the holes. (a-c) Band structure near $\omega_0= 0.38749 \times 2\pi c/a$ as a typical example of conventional anisotropic bands that are doubly degenerate at $\Gamma$. (a) Band dispersion  diagram along $\Gamma\textrm{-}X$ and $\Gamma\textrm{-}M$. (b) Constant frequency contours of the lower band with respect to $(k_x,k_y)$. (c) Constant frequency contours of the upper band. (d-f) Band structure at $\omega_0= 0.47656 \times 2\pi c/a$, which is nearly isotropic. (d) Band dispersion along $\Gamma\textrm{-}X$ and $\Gamma\textrm{-}M$. (e) Constant frequency contours of the lower band. (f) Constant frequency contours of the upper band. The frequency $(\omega-\omega_0)$ is in the units of $10^{-4}\times2\pi c/a$, while $|\bm{k}|$, $k_x$ and $k_y$ are in the units of $10^{-3}\times2\pi/a$.}
\end{figure*}

The analysis above indicates that, to use a photonic crystal slab as a two-dimensional spatial differentiator, it is sufficient that the slab satisfies the following three conditions:
\begin{enumerate}
\item $t_d = 1$.
\item Only one guided resonance band is coupled.
\item The band satisfies the dispersion $(\omega(\bm{k})-\omega_0) \propto |\bm{k}|^2$.
\end{enumerate}

To satisfy the first condition above, we note that the direct transmission coefficient $t_d$ is related to the non-resonant transmission pathway \cite{Fan2002}. Hence it is possible to realize $t_d = 1$ by, e.g. changing the thickness of the slab. In the structure as shown in Figure~\ref{fig1}, we achieve $t_d = 1$ by placing a uniform dielectric slab in the vicinity of the photon crystal slab, and by tuning the distance between the slabs. This has the advantage that we can tune $t_d$ without significantly affecting the band structure of the photonic crystal slab.

To satisfy the second and third conditions above, one will need to design the band structure for the photonic crystal slab. The design here is in fact quite non-trivial due to the vectorial nature of  electromagnetic waves. Since the Laplacian is isotropic in $\bm{k}$-space, it is natural to consider a photonic crystal slab structure that has rotational symmetry. As an illustration, here we consider a slab structure with a square lattice of air holes that has $C_{4v}$ symmetry. For such a slab, it is known that  at the $\Gamma$ point, which corresponds to $|\bm{k}|=0$, the only modes that can couple to external plane wave must be two-fold degenerate, belonging to a two-dimensional irreducible representation of the $C_{4v}$ group.\cite{Ochiai2001,Fan2002} Near such  modes, in the vicinity of $\Gamma$ point, the band structure in general can be described by the following $2 \times 2$ effective Hamiltonian: (See Supplementary Materials)
\begin{equation} \label{eq:H-k}
\hat{\mathcal{H}}(\bm{k}) = (\omega_0 - i \gamma_0 + a |\bm{k}|^2)\hat{\mathbf{I}} + b (k_x^2-k_y^2) \hat{\bm{\sigma}}_z+ c k_x k_y \hat{\bm{\sigma}}_x~,
\end{equation}
where $a,b,c$ are three complex coefficients, the $\bm{\sigma}$'s are the Pauli matrices. This Hamiltonian has two eigenvalues of
\begin{equation}
\omega_{\pm}(\bm{k}) - i \gamma_{\pm}(\bm{k}) = \omega_0 - i \gamma_0 + a|\bm{k}|^2 \pm \sqrt{b^2(k_x^2-k_y^2)^2+ c^2 k_x^2k_y^2} ~.\label{eq:eigenvalue}
\end{equation}

The band structure as described by Equation~(\ref{eq:eigenvalue}) in general does not satisfy the conditions for ideal spatial differentiation as outlined above. In this paper for concreteness we will fix the dielectric constant of the material for the slabs to be  $\epsilon = 12$, which approximates that of Si or GaAs in the infrared wavelength range. In Figure~\ref{fig2}(a-c), we plot the band structure for a  slab with a thickness $d=0.55a$, and a radius  $r=0.111a$  of the holes, in the vicinity of a guided resonance at the $\Gamma$ point with a frequency $\omega_0= 0.38749 \times 2\pi c/a$. The band structure is numerically determined using guided-mode expansion method \cite{Andreani2006,Minkov2014}, and it agrees excellently with the analytic expression of Equation~(\ref{eq:eigenvalue}) (See Supplementary Materials).  At the $\Gamma$ point, the bands are two-fold degenerate. Away from the $\Gamma$ point, the degeneracy is lifted, and the two bands are strongly anisotropic, as can be seen in Figure~\ref{fig2}(a), where the effective masses along the $\Gamma \textrm{-}X$ and the $\Gamma \textrm{-}M$ direction are drastically different. The band anisotropy can also be visualized in the constant frequency contour plotted in Figure~\ref{fig2}(b) and (c), for the two bands. In general, it can be shown that, for a general choice of parameters $a$, $b$ and $c$, the constant frequency contour as described by Equation (\ref{eq:eigenvalue}) is not a circle even at the $|\bm{k}|\to0$ limit.

On the other hand, Equation~(\ref{eq:eigenvalue}) also indicates that when $c=\pm2b$, both bands will be isotropic:
\begin{equation}
\omega_{\pm}(\bm{k}) - i \gamma_{\pm}(\bm{k}) = \omega_0 - i \gamma_0 + (a\pm b)|\bm{k}|^2 ~.\label{eq:eigenvalue_iso}
\end{equation}
To achieve such an isotropic band structure requires detailed tuning of the parameters. Numerically, we found that such an isotropic band structure can be approximately achieved with the same slab as indicated above, but near the guided resonance at $\Gamma$ with a different frequency $\omega_0 = 0.47656 \times 2\pi c/a$, where the complex coefficients $a=0.68-0.14i,b=1.11-0.12i,c=2.24+0.15i$,  thus $c/2b=0.99 + 0.17i\approx1$. (See Supplementary Materials for details.) The resulting isotropic band structures are shown in Figure~\ref{fig2}(d-f). In Figure~\ref{fig2}(d), we see that the bands along the $\Gamma \textrm{-}M$ and $\Gamma \textrm{-}X$ direction have almost identical effective masses. And in Figure~\ref{fig2}(e) and (f), we see that the constant frequency contours are almost completely circular.

As we mentioned above, for a structure with $C_{4v}$ symmetry, the modes that a normally-incident plane wave can couple to at the $\Gamma$ point are always two-fold degenerate. Consequently, near the $\Gamma$ point, there are always two bands of guided resonance present. Moreover, off the normal direction, if the direction of incident waves is away from the high symmetry planes, both $S$ and $P$ polarized lights may couple to both bands, leading to complex polarization conversion effects. Thus,  in general, in addition to having an anisotropic band structure near $\Gamma$, a photonic crystal slab structure also does not satisfy the condition 2 above regarding the excitation of a single guided resonance band. Remarkably, however, below we show that once the condition for an isotropic band is satisfied, each of the two bands in fact only  couple to one single polarization, \textit{for every direction of incidence}. Below, we refer to this effect, where each polarization excites only a single band away from normal incidence, as the effect of \textit{single-band excitation}.

To illustrate the effect of single-band excitation for every direction when the bands are isotropic, we first show that single-band excitation always occurs when the incident direction is in a mirror symmetry plane of the structure. Using the group theory notation in Ref.~\citep{Sakoda2005}, as our system possesses $C_{4v}$ symmetry, the doubly degenerate states at $\Gamma$ are $E$ modes. Along the $\Gamma\textrm{-}X$ direction, this pair of $E$ modes splits into two singly degenerate states of $A$ and $B$ modes which are even and odd, respectively, with respect to the reflection operator associated with the mirror plane $y=0$. On the other hand, since the $S$ and $P$ polarized lights are also odd and even with respect to the same mirror plane, respectively, along the $\Gamma\textrm{-}X$ direction, the $S$($P$) polarized light can only couple to the $B$($A$) modes. Thus, in general, we have the effect of single-band excitation when the direction of incidence is in a high symmetry plane such as the $y=0$ plane.\cite{Tikhodeev2002}

Next, we prove the following statement: if the two-band Hamiltonian is isotropic, that is,
\begin{equation}\label{eq:hamilton}
 \quad \hat{\mathcal{R}}(\varphi) \hat{\mathcal{H}}(\bm{k}) \hat{\mathcal{R}}^{-1}(\varphi) = \hat{\mathcal{H}}(\bm{R}(\varphi)\bm{k})~,
\end{equation}
for every $\varphi \in (0, 2\pi)$, then we have the effect of single-band excitation along all directions. Here, $\bm{R}(\varphi)$ and $\hat{\mathcal{R}}(\varphi)$ are the rotation operators that describe the rotation around the $z$-axis by an angle $\varphi$, in the $\bm{k}$ space and the two-dimensional Hilbert space, respectively. To prove this, we denote the eigenstates of the two bands as $\ket{\bm{k}, A}$ and $\ket{\bm{k}, B}$, which connect to the $A$ and $B$ modes along the $\Gamma \textrm{-}X$ direction respectively. Equation~(\ref{eq:hamilton}) implies that:
\begin{equation}\label{eq:eigenmodes} \begin{split}
	\ket{\bm{R}(\varphi)\bm{k}, A} &= \hat{\mathcal{R}}(\varphi) \ket{\bm{k}, A}~, \\
    \ket{\bm{R}(\varphi)\bm{k}, B} &= \hat{\mathcal{R}}(\varphi) \ket{\bm{k}, B}~.
\end{split}
\end{equation}
We denote the $S$ and $P$ polarized modes as $\ket{\bm{k}, S}$ and $\ket{\bm{k}, P}$ respectively. By definition,
\begin{equation}\label{eq:SP} \begin{split}
	\ket{\bm{R}(\varphi)\bm{k}, S} &= \hat{\mathcal{R}}(\varphi) \ket{\bm{k}, S}~, \\
    \ket{\bm{R}(\varphi)\bm{k}, P} &= \hat{\mathcal{R}}(\varphi) \ket{\bm{k}, P} ~.
\end{split}
\end{equation}
For $\bm{k}=\bm{k_x}$ along the $\Gamma \textrm{-}X$ direction, $A$($B$) mode does not couple to the $S$($P$) polarizations, i.e.
\begin{equation} \begin{split}
\braket{\bm{k_x}, S|\bm{k_x}, A}&=0~, \\
\braket{\bm{k_x}, P|\bm{k_x}, B}&=0~.
\end{split}
\end{equation}
Then for $\bm{k} = \bm{R}(\varphi)\bm{k_x}$ along any direction, using Equations (\ref{eq:eigenmodes}) and (\ref{eq:SP}), we also have
\begin{equation} \begin{split}
\braket{\bm{k}, S|\bm{k},A}= \bra{\bm{k_x}, S}\inv{\hat{\mathcal{R}}}(\varphi)\hat{\mathcal{R}}(\varphi)\ket{\bm{k_x}, A} = 0 ~, \\
\braket{\bm{k}, P|\bm{k},B}= \bra{\bm{k_x}, P}\inv{\hat{\mathcal{R}}}(\varphi)\hat{\mathcal{R}}(\varphi)\ket{\bm{k_x}, B} = 0 ~.
\end{split}
\end{equation}

\begin{figure}
\centering
\includegraphics[width=1.0\columnwidth]{Fig3.png}
\caption{\label{fig3}Transmittances of the two-slab structure as illustrated in Figure~\ref{fig1}. The photonic crystal slab has the same structure parameters as those of Figure~\ref{fig2}. The thickness of the uniform dielectric slab $d_s=0.07a$. The air gap width $d_g=2.1a$. (a,b) Transmittance near $\omega_0= 0.38749 \times 2\pi c/a$ as a function of $|\bm{k}|$ and frequency $(\omega-\omega_0)$ along a general direction $\varphi=14^\circ$ for (a) $S$ light and (b) $P$ light. Both bands couple to $S$ and $P$ light. (c,d) Transmittance near $\omega_0= 0.47656 \times 2\pi c/a$ along $\varphi=14^\circ$ for (c) $S$ light and (d) $P$ light. $S$ light only couples with the upper $B$ band while $P$ with the lower $A$ band, showing the single-band excitation effect. Frequency $(\omega-\omega_0)$ is in the units of $10^{-4}\times2\pi c/a$, while $|\bm{k}|$ is in the units of $10^{-3}\times2\pi/a$.}
\end{figure}
As a numerical demonstration of the connection of the single-band excitation effect with the isotropic band structure, we consider the two-slab structure as illustrated in Figure~\ref{fig1}. The photonic crystal slab structure has the same parameters as those of Figure~\ref{fig2}. The uniform dielectric slab has a thickness of $d_s=0.07a$. The air gap between the two slabs has a width of $d_g=2.1a$. Figure~\ref{fig3} shows the transmission coefficients through the two slab system, as we vary the in-plane wavevector, while maintaining $\varphi = 14^\circ$, and hence the direction of the incident light is away from any symmetry plane. Near  the frequency of $\omega = 0.38749 \times 2\pi c/a$, which is near the guided resonances as shown in Figure~\ref{fig2}(a-c), the band structure is strongly anisotropic, and both the $S$ and $P$ polarizations excite both bands as shown in Figure~\ref{fig3}(a-b). On the other hand, near the frequency of $\omega = 0.47656 \times 2\pi c/a$, which is near the guided resonances as shown in Figure~\ref{fig2}(d-f), each polarization excites only one band as shown in Figure~\ref{fig3}(c-d).

In this structure, the thickness of the uniform dielectric slab and the gap size are chosen such that $t_d = 1$ near $\omega = 0.47656 \times 2\pi c/a$. Therefore, we have designed a structure that satisfies all the sufficient conditions for achieving an ideal Laplacian spatial differentiator.

In Figure~\ref{fig4}(a) and (b), we plot the transmission coefficients, as a function of in-plane wavevector $\bm{k}$, for the $S$ and $P$ polarized incident light, at the frequency $\omega = 0.47656 \times 2\pi c/a$, for the same structure used in Figure~\ref{fig3}. Both polarizations exhibit a transmission coefficient that is proportional to $|\bm{k}|^2$. In Figure~\ref{fig4}(c), we plot the transmission coefficient for unpolarized light at the same frequency. The transmission coefficient for the unpolarized light can be derived as (See Supplementary Materials for a derivation)
\begin{equation}
|t_u| = \sqrt{\frac{|t_s|^2+|t_p|^2}{2}}~.
\end{equation}
Here $t_s$ and $t_p$ are the transmission coefficients for $S$ and $P$ light respectively. The transmission coefficient for the unpolarized light also shows an isotropic response. In Figure~\ref{fig4}(d), we plot $|t_u|$ as a function of $|\bm{k}|$,  which shows a quadratic dependency. To show this, we fit the curve with quadratic function $|t_u|=\alpha|\bm{k}|^2$, where $\alpha$ is the only fitting parameter. The fitting is almost perfect for $|\bm{k}|$ up to $6\times10^{-3}\times2\pi/a$, which as we will show provides a sufficient range of wavevector for image differentiation. Thus, the transmission coefficients have all the required $\bm{k}$-dependency in order to demonstrate a Laplacian.
\begin{figure}
\centering
\includegraphics[width=1.0\columnwidth]{Fig4.png}
\caption{\label{fig4}(a,b,c) Contour plot of transmittance $|t|$ as a function of $k_x$ and $k_y$ at frequency $\omega_0 = 0.47656 \times 2\pi c/a$ for (a) $S$, (b) $P$ and (c) unpolarized light, which are all isotropic near $\Gamma$. (d) $|t_u|$ as a function of $|\bm{k}|$ along $\varphi=14^\circ$ direction, and the quadratic fitting $t_u = \alpha |\bm{k}|^2$ where $\alpha$ is the only fitting parameter. $k_x$ and $k_y$ are in the units of $10^{-3}\times2\pi/a$.}
\end{figure}
\section{Numerical Demonstration}
\label{sec:demo}

\begin{figure}
\includegraphics[width=1.0\columnwidth]{Fig5.png}
\caption{\label{fig5}(a) Incident Stanford emblem with a size of $2610 a\times 1729 a$. (b) Calculated transmitted image, which clearly shows the edges with different orientations. (c) Incident slot patterns with length $500 a$ and width $100, 50, 30$ and $20 a$ respectively. (d) Calculated transmitted images, which show that  the spatial resolution of our design is around $30 a$.}
\end{figure}

We now show that the structure used in Figure~\ref{fig3} indeed operates as a two-dimensional Laplacian differentiator. Figure~\ref{fig5}(a) is the Stanford emblem as the incident image. Figure~\ref{fig5}(b) is the calculated transmitted image, which clearly shows all the edges with the same intensity despite different edge orientations. The ability to detect edges with different orientations is a significant advantage of an isotropic differential operator such as the two-dimensional Lapacian, as compared with the one-dimensional differential operators. The Laplacian also highlights the sharp corners, exhibiting its higher sensibility to fine details as a second-order derivative operator.\cite{Gonzales2008}

We also test the spatial resolution of such spatial differentiator. Figure~\ref{fig5}(c) is a series of slot patterns and Figure~\ref{fig5}(d) is the calculated differentiated image. The performance of edge detection degrades as the slot width decreases. The spatial resolution of our differentiator, which is the minimum separation between the two edges that can be resolved, is around $30 a$. Considering $a=0.67\mu m$ corresponding to the resonant wavelength $\lambda_0=1.4 \mu m$, the spatial resolution is $20 \mu m$, which is sufficient for most image processing application.
\section{Discussion and Conclusion}
In previous works, the two-dimensional Laplace operator has been implemented using holograms \cite{Guo2008,Ryle2010} and phase-shifted Bragg grating \cite{Bykov2014a}. These implementations either rely upon bulky optical components, or operate in reflection modes. In contrast, here we show this operation can be achieved with a compact optical structure in transmission mode, which is more compatible for image processing applications.

In conclusion, we have shown that a Laplace differentiator can be implemented using  a photonic crystal slab. Such a simple photonic device
may have various applications involving image processing. Future research will be worthy to extend the Laplace differentiator to multi-frequency with multi-layer structures, enabling differentiation of color images.

\section*{Funding Information}
This work is supported in part by Samsung Electronics, and by the U. S. Air Force Grant No.  FA9550-17-1-0002.

\section*{Acknowledgments}

The authors thank Yu Guo and Dr. Alexander Cerjan for helpful discussions.

\section*{Supplementary Materialsal Documents}

\bigskip \noindent See \href{link}{Supplementary Materials} for supporting content.

\section{Derivation of the effective Hamiltonian}\label{sec:Hamiltonian}
In this section, we derive the $2\times2$ effective Hamiltonian (Equation~(10) in the primary text) near the $\Gamma$ point. We assume that the system has $C_{4v}$ symmetry. And at the $\Gamma$ point the system supports a pair of doubly degenerate states, denoted as $\ket{x}$ and $\ket{y}$, respectively. With these states as bases, the $2\times2$ Hamiltonian in the vicinity of $\Gamma$ has the following general form:
\begin{equation}
\hat{\mathcal{H}}(\bm{k}) = \hat{\mathcal{A}}(\bm{k}) - i \hat{\mathcal{B}}(\bm{k})
\end{equation}
where $\hat{\mathcal{A}}(\bm{k})$ and $\hat{\mathcal{B}}(\bm{k})$ are both Hermitian and
\begin{align}
\hat{\mathcal{A}}(\bm{k}) &= f(\bm{k})\hat{\bm{\sigma}}_{+}+f^*(\bm{k})\hat{\bm{\sigma}}_{-}+g(\bm{k})\hat{\bm{\sigma}}_{z}+h(\bm{k})\hat{I}+\omega_0\hat{I} \notag \\
\hat{\mathcal{B}}(\bm{k}) &= r(\bm{k})\hat{\bm{\sigma}}_{+}+r^*(\bm{k})\hat{\bm{\sigma}}_{-}+s(\bm{k})\hat{\bm{\sigma}}_{z}+t(\bm{k})\hat{I}+\gamma_0\hat{I}
\end{align}
We are interested in the lowest-order non-vanishing terms in the functions $f, g, h, r, s$ and $t$. By virtual of the two-fold degeneracy, we have
\begin{equation}
f(\bm{0})=g(\bm{0})=h(\bm{0})=r(\bm{0})=s(\bm{0})=t(\bm{0})=0
\end{equation}
Also, $g,h,s$ and $t$ are real since $\hat{\mathcal{A}}(\bm{k})$ and $\hat{\mathcal{B}}(\bm{k})$ are Hermitian.

In general, the  Hamiltonian $\hat{\mathcal{H}}(\bm{k})$ is constrained by the symmetric group $G$ of the system:
\begin{equation} \label{eq:constraint}
\forall g\in G,\qquad \hat{\mathcal{D}}(g) \hat{\mathcal{H}}(\bm{k}) \hat{\mathcal{D}}^{-1}(g) = \hat{\mathcal{H}}(\bm{D}(g)\bm{k})
\end{equation}
where $\hat{\mathcal{D}}(g)$ and $\bm{D}(g)$ are the representations of $g$ in the Hilbert space and the $\bm{k}$ space, respectively.

For our system, $G=C_{4v}=\{E, 2C_4, C_2, 2\sigma_v, 2\sigma_d\}.$ As we choose $\ket{x}$ and $\ket{y}$ as bases for the Hilbert space and $\bm{k_x}$ and $\bm{k_y}$ for the $\bm{k}$ space, the representations $\hat{\mathcal{D}}(g)$ and $\bm{D}(g)$ have the same matrix forms. Now, we consider the constraint on $\hat{\mathcal{H}}(\bm{k})$ from each symmetry element in $G$.
\begin{enumerate}
\item {\bf Inversion symmetry $C_2$.}
\begin{equation}
\hat{\mathcal{D}}(C_2) = \bm{D}(C_2) \equiv \bm{\Pi}= \begin{pmatrix}
-1 & 0 \\
0 & -1
\end{pmatrix}
\end{equation}
From Equation~(\ref{eq:constraint}) we have:
\begin{align}
&f(\bm{k})=f(-\bm{k}),\quad g(\bm{k})=g(-\bm{k}),\quad h(\bm{k})=h(-\bm{k}) \notag \\
&r(\bm{k})=r(-\bm{k}),\quad s(\bm{k})=s(-\bm{k}),\quad t(\bm{k})=t(-\bm{k})
\end{align}
Thus, the lowest order terms in $f, g, h, r, s$ and $t$ must all be quadratic.
\item {\bf Four-fold rotational symmetry $C_4$.}
\begin{equation}
\hat{\mathcal{D}}(C_4) = \bm{D}(C_4) \equiv \bm{R} = \begin{pmatrix}
0 & 1 \\
-1 & 0
\end{pmatrix}
\end{equation}
From Equation~(\ref{eq:constraint}) we have:
\begin{align}
&f^*(\bm{k})=-f(\bm{R}\bm{k}),\quad g(\bm{k})=-g(\bm{R}\bm{k}),\quad h(\bm{k})=h(\bm{R}\bm{k}) \notag \\
&r^*(\bm{k})=-r(\bm{R}\bm{k}),\quad s(\bm{k})=-s(\bm{R}\bm{k}),\quad t(\bm{k})=t(\bm{R}\bm{k})
\end{align}
\item {\bf Mirror symmetry $\sigma_v$.}
\begin{equation}
\hat{\mathcal{D}}(\sigma_v) = \bm{D}(\sigma_v) \equiv \bm{M_y} = \begin{pmatrix}
1 & 0 \\
0 & -1
\end{pmatrix}
\end{equation}
Here we consider the mirror plane perpendicular to the y-axis. From Equation~(\ref{eq:constraint}) we have:
\begin{align}
&f(\bm{k})=-f(\bm{M_y}\bm{k}),~ g(\bm{k})=g(\bm{M_y}\bm{k}),~ h(\bm{k})=h(\bm{M_y}\bm{k}) \notag \\
&r(\bm{k})=-r(\bm{M_y}\bm{k}),~ s(\bm{k})=s(\bm{M_y}\bm{k}),~ t(\bm{k})=t(\bm{M_y}\bm{k})~
\end{align}
\end{enumerate}

The other mirror symmetry $\sigma_d$ doesn't provide new constraints as it can be realized with a combination of $\sigma_v$ and  $C_{4v}$.

Combining all the symmetry requirements, we determine the forms of the lowest order terms in all the functions. As an example, we consider the function $f(\bm{k})$, which can be expanded as
\begin{equation}
f(\bm{k}) = C_x k_x^2 + C_y k_y^2 + C_{xy} k_x k_y
\end{equation}
where all the $C$ coefficients are in general complex. Since $f(\bm{k})=-f^*(\bm{R}\bm{k})=-f(\bm{M_y}\bm{k})$, we have $C_x = C_y = 0$, $C_{xy}=C^*_{xy}\equiv C$. Thus $f(\bm{k})=C k_x k_y$, where $C$ is real. Repeating the procedures for all the functions, we summarize the results below:
\begin{align}
&f(\bm{k}) = C k_x k_y,\qquad &r(\bm{k})= C' k_x k_y \notag\\
&g(\bm{k}) = B (k_x^2-k_y^2),\qquad &s(\bm{k}) = B' (k_x^2-k_y^2) \notag\\
&h(\bm{k}) = A (k_x^2+k_y^2),\qquad &t(\bm{k}) = A' (k_x^2+k_y^2)
\end{align}
where all the coefficients are real.

Therefore, we obtain the Hamiltonian in Equation~(10) of the main text.:
\begin{equation} \label{eq:H-k}
\hat{\mathcal{H}}(\bm{k}) = (\omega_0 - i \gamma_0 + a |\bm{k}|^2)\hat{\mathbf{I}} + b (k_x^2-k_y^2) \hat{\bm{\sigma}}_z+ c k_x k_y \hat{\bm{\sigma}}_x
\end{equation}
where $a=A-iA',~b=B-iB',~c=C-iC'$ are the complex coefficients.

\section{Validation of the effective Hamiltonian}\label{sec:Validation_Hamiltonian}
In this section we provide quantitative validation of the effective Hamiltonian. Figure (\ref{fig:SI_model_anisotropy}) plots the band structure of the photonic crystal slab near frequency $\omega_0=0.38749 \times 2\pi c/a$. Figure (\ref{fig:SI_model_anisotropy}) (a-d) exhibit contour plots of $(\omega-\omega_0)$ and $(\gamma-\gamma_0)$ for the two bands numerically calculated using guided-mode expansion method, while (e-h) show the corresponding analytical results from the effective Hamiltonian (Equation (\ref{eq:H-k})). The complex coefficients $a=0.33-0.16i,b=0.24-0.16i,c=8.78-0.08i$ are fitted from the band dispersions along $\Gamma\textrm{-}X$ and $\Gamma\textrm{-}M$ direction. Figure (\ref{fig:SI_model}) plots the numerical and analytical results of the band structure of the same photonic crystal slab near frequency $\omega_0=0.47656 \times 2\pi c/a$.  The fitted complex coefficients are $a=0.68-0.14i,b=1.11-0.12i,c=2.24+0.15i$. In both cases, the numerical results agree with the analytical expressions excellently, which provide the validations of the effective Hamiltonian.

\begin{figure}
\includegraphics[width=0.9\columnwidth]{FIG_S2_anisotropy.png}
\centering
\caption{\label{fig:SI_model_anisotropy}Band structure of the photonic crystal slab near frequency $\omega_0=0.38749 \times 2\pi c/a$. The structure has the dielectric constant  $\epsilon=12$, the thickness $d=0.55a$, and the radius $r=0.111a$ of the holes. (a-d) Contour plots of the band structure numerically calculated using guided-mode expansion method. (a) $(\omega-\omega_0)\textrm{-}(k_x,k_y)$ for the lower band. (b) $(\omega-\omega_0)\textrm{-}(k_x,k_y)$ for the upper band. (c) $(\gamma-\gamma_0)\textrm{-}(k_x,k_y)$ for the lower band. (d) $(\gamma-\gamma_0)\textrm{-}(k_x,k_y)$ for the upper band. $\gamma_0= 2.0\times10^{-4}\times 2\pi c/a$. The $(\omega-\omega_0)$ and $(\gamma-\gamma_0)$ contours for both bands are anisotropic. (e-h) Corresponding contour plots of the band structure analytically calculated from the effective Hamiltonian (Equation (\ref{eq:H-k})). The complex coefficients $a=0.33
-0.16i,b=0.24-0.16i,c=8.78
-0.08i$ are fitted from the band dispersions along $\Gamma\textrm{-}X$ and $\Gamma\textrm{-}M$ direction. The numerical results agree with the analytical expressions excellently. In all the plots $k_x$ and $k_y$ are in the units of $10^{-3}\times2\pi/a$, $(\omega-\omega_0)$ and $(\gamma-\gamma_0)$ are in the units of $10^{-4}\times2\pi c/a$.
}
\end{figure}
\begin{figure}
\includegraphics[width=0.9\columnwidth]{FIG_S2.png}
\centering
\caption{\label{fig:SI_model}Band structure of the photonic crystal slab near frequency $\omega_0=0.47656 \times 2\pi c/a$. The structure has the dielectric constant  $\epsilon=12$, the thickness $d=0.55a$, and the radius $r=0.111a$ of the holes. (a-d) Contour plots of the band structure numerically calculated using guided-mode expansion method. (a) $(\omega-\omega_0)\textrm{-}(k_x,k_y)$ for the lower band. (b) $(\omega-\omega_0)\textrm{-}(k_x,k_y)$ for the upper band. (c) $(\gamma-\gamma_0)\textrm{-}(k_x,k_y)$ for the lower band. (d) $(\gamma-\gamma_0)\textrm{-}(k_x,k_y)$ for the upper band. $\gamma_0= 1.3\times10^{-4}\times 2\pi c/a$. The $(\omega-\omega_0)$ contours for both bands are almost completely circular, while the $(\gamma-\gamma_0)$ contours are anisotropic. Nonetheless, $(\gamma-\gamma_0)$ are much smaller than $(\omega-\omega_0)$ and thus don't affect the isotropy of the band structure much. (e-h) Corresponding contour plots of the band structure analytically calculated from the effective Hamiltonian (Equation (\ref{eq:H-k})). The complex coefficients $a=0.68-0.14i,b=1.11-0.12i,c=2.24+0.15i$ are fitted from the band dispersions along $\Gamma\textrm{-}X$ and $\Gamma\textrm{-}M$ direction. The numerical results agree with the analytical expressions excellently. In all the plots $k_x$ and $k_y$ are in the units of $10^{-3}\times2\pi/a$, $(\omega-\omega_0)$ and $(\gamma-\gamma_0)$ are in the units of $10^{-4}\times2\pi c/a$.
}
\end{figure}
\section{Transmittance of arbitrary polarized states}\label{sec:DensityMatrix}
In this section we compute the transmission of an arbitrarily polarized state,  given the transmission response of the system for $S$ and $P$ polarized light.

The transmission response of the system is described by the $S$ matrix:
\begin{equation}
S = \begin{pmatrix}
t_{ss} & t_{sp} \\
t_{ps} & t_{pp}
\end{pmatrix}
\end{equation}
where $t_{ss}, t_{sp}, t_{ps}$ and $t_{pp}$ are complex. Here we adopt $\ket{s}$ and $\ket{p}$ as the basis states. For a general input state $\ket{e}$, the output state is $S\ket{e}$.

Consider an input beam with its polarization described by a  polarization density matrix $\rho$, normalized such that $I = \operatorname{tr} \rho$ is the incident light intensity. The transmitted state is then described by $\rho' = S \rho S^\dagger$,
with the transmission defined as:
\begin{equation} \label{eq:T}
T = \operatorname{tr} \rho' / \operatorname{tr} \rho
\end{equation}

Now let's consider some specific  examples of incident states:
\begin{enumerate}
\item {\bf S polarized light}
\begin{align}
\rho &= I \begin{pmatrix}
1 & 0 \\
0 & 0
\end{pmatrix}\notag \\
\rho' &= I \begin{pmatrix}
|t_{ss}|^2 & t_{ss} t_{ps}^* \\
t_{ss}^* t_{ps} & |t_{ps}|^2
\end{pmatrix}\notag\\
T &\equiv |t_s|^2 = |t_{ss}|^2+|t_{ps}|^2 \label{eq:ts}
\end{align}
Similarly, for P polarized light,
\begin{equation}
|t_p|^2 = |t_{sp}|^2+|t_{pp}|^2 \label{eq:tp}
\end{equation}

\item {\bf Unpolarized light}
\begin{align}
\rho &= \frac{I}{2} \begin{pmatrix}
1 & 0 \\
0 & 1
\end{pmatrix}\notag\\
\rho' &= \frac{I}{2}  \begin{pmatrix}
|t_{ss}|^2 + |t_{sp}|^2 & t_{ss} t_{ps}^* + t_{sp}t_{pp}^* \\
t_{ss}^* t_{ps}+ t_{sp}^*t_{pp} & |t_{ps}|^2 + |t_{pp}|^2
\end{pmatrix}\notag\\
T &\equiv |t_u|^2 = \frac{1}{2}(|t_{ss}|^2+|t_{ps}|^2+|t_ps|^2+|t_{pp}|^2) \label{eq:tu}
\end{align}

From Equation~(\ref{eq:ts}-\ref{eq:tu}) we have:
\begin{equation}
|t_u| = \sqrt{\frac{|t_s|^2+|t_p|^2}{2}}
\end{equation}
which is Equation~(18) in the main text.

\item {\bf Left Circularly polarized light}
\begin{align}
\rho &= \frac{I}{2} \begin{pmatrix}
1 & -i \\
i & 1
\end{pmatrix}\notag\\
\rho' &= \frac{I}{2}  \begin{pmatrix}
|t_{ss}|^2 + |t_{sp}|^2 + i(t_{sp} t_{ss}^* - t_{sp}^* t_{ss}) & t_{ss} t_{ps}^* + t_{sp}t_{pp}^* + i(t_{sp} t_{ps}^*-t_{ss}t_{pp}^*) \\
t_{ss}^* t_{ps}+ t_{sp}^*t_{pp} - i(t_{sp}^* t_{ps}-t_{ss}^*t_{pp}) & |t_{ps}|^2 + |t_{pp}|^2 +i(t_{ps}^* t_{pp} - t_{ps} t_{pp}^*)
\end{pmatrix}\notag\\
T & \equiv |t_{l}|^2 = \frac{1}{2}[|t_{ss}|^2+|t_{ps}|^2+|t_ps|^2+|t_{pp}|^2 + i (t_{sp}t_{ss}^*-t_{sp}^*t_{ss}+t_{ps}^*t_{pp}-t_{ps}t_{pp}^*)]\label{eq:tl}
\end{align}
which differs from Equation~(\ref{eq:tu}) only in the interference term.
\end{enumerate}

Note that for our system with an isotropic band structure, $S$ is diagonal:
\begin{equation}
S = \begin{pmatrix}
t_{ss} & 0 \\
0 & t_{pp}
\end{pmatrix}
\end{equation}
Then the interference term in Equation~(\ref{eq:tl}) disappears and we have $|t_l|=|t_u|$. So the differentiation performance will be the same under illumination of unpolarized or circularly polarized light.

\section{Transmittance at larger wavevectors}\label{sec:t_k}
\begin{figure}
\includegraphics[width=0.9\columnwidth]{SI-tmap.png}
\centering
\caption{\label{fig:t_map}Transmittance $|t|$ as a function of $k_x$ and $k_y$ at frequency $\omega_0 = 0.47656 \times 2\pi c/a$ for (a) $S$, (b) $P$ and (c) unpolarized light, which are all isotropic near $\Gamma$, which is the major relevant wavevector region for image processing, but anisotropic at larger $\bm{k}$. $k_x$ and $k_y$ are in the units of $10^{-2}\times2\pi/a$.}
\end{figure}
Figure~(\ref{fig:t_map}) plots the transmittance $|t|$ as a function of $k_x$ and $k_y$ for $S$, $P$ and unpolarized light in a larger wavevector region as compared to Figure 4 of the primary text. All $|t|$ are isotropic near the $\Gamma$ point. They become anisotropic at larger wavevector $\bm{k}$. The background transmittances are all $1$. The occurrence of transmittance dips at larger $\bm{k}$ are due to other guided resonance modes.

We calculate the performance of our differentiator based on the transmittance shown here, which covers all the relevant wavevector regions for the incident images that we assumed.

\end{document}


\section{Introduction}

Coherent optical frequency combs have
revolutionized
precision  measurement
with light
\cite{Udem1999_1,Diddams2000,Jones2000}.
However,
the complex setups needed to generate
and stabilize
these systems
limits their
relevance
in field-testing scenarios.
Combs generated in a
micro-cavity
\cite{DelHaye2007, Kippenberg2011}
(``microcombs”)
have been evolving rapidly,
and
provide a pathway to miniaturize frequency comb
systems. While early microcombs were subject to
instabilities and lacked the critical ability to form
femtosecond pulses, an important advancement has been
the realization of temporal soliton mode locking in optical
microcavities
\cite{Herr2014}.
First observed in optical fibers
\cite{Leo2010},
this form of mode locking has been demonstrated
across many microcavity platforms using
materials such as
silica (SiO$_2$) \cite{Yi2015}
and
silicon nitride (Si$_3$N$_4$)
\cite{Brasch2016_1, Huang2016, Wang2016, Joshi2016}.
The potential to fully integrate frequency comb systems is
driving interest in semiconductor platforms that leverage
complementary metal-oxide-semiconductor (CMOS)
fabrication infrastructures
\cite{Heck2013, Roelkens2015}.
These are also
of interest
as soliton stabilization
\cite{Yi2016, Brasch2016_2}
necessitates
certain electronic components.

Temporal
optical
solitons
\cite{Hasegawa1973_1}
were
first observed in
optical fibers
\cite{Mollenauer1980}.
These
nonlinear
waves
balance
second-order
dispersion
using the
optical Kerr effect.
For
the new
temporal solitons in
microcavities, an
second balance
also occurs
\cite{Herr2014}:
optical loss
is
compensated by
third-order parametric gain
\cite{Kippenberg2004}
provided by
an external pump.
The resulting dissipative Kerr soliton system is
therefore a regenerative mode-locked oscillator.
Soliton
mode locked microcavities have now been applied to
demonstrate dual-comb
spectroscopy \cite{Suh2016,YuM2017,Pavlov2017},
optical
frequency synthesis
\cite{Spencer2017},
distance ranging
\cite{Suh2017,Trocha2017}
and optical communications
with tremendous bandwidth
\cite{Marin-Palomo2017}.

External-cavity diode lasers (ECDLs)
are
critical
for many applications
and are frequently used to optically pump microcomb systems.
They can provide single mode and narrow linewidth,
which is crucial for telecommunications.
In addition, their emission frequency can be tuned,
which is needed
for spectroscopy \cite{Schilt2004}
and frequency synthesis \cite{Spencer2017}.
However, they
typically emit
modest output power.
Consequently,
all soliton
microcomb
systems to date
require
non-integrated components
to boost the pump laser power.
Additionally, other functions to enable
rapid control of laser power and frequency
(acoustic/electro-optic modulators)
\cite{Brasch2016_1,Yi2016,Stone2017arXiv}
have also been
necessary.
In this work,
we report an ECDL
with extremely narrow linewidth,
low relative intensity noise (RIN),
high output power, and
useful spectral tunability.
As a very first application for this prototype,
a temporal soliton
is directly generated and stabilized
in a
high-$Q$
silica micro-resonator,
without using an optical amplifier or an external modulator.
This drastically reduces the complexity, size,
and cost of the system,
and
is expected to precede
a fully integrated source of solitons.

\end{spacing}

\section{Methods}
\label{SEC_METHOD}

\subsection{Summary of temporal soliton theory}
\label{SEC_2A}

Optical frequency combs
can be generated by pumping
a micro-resonator with
a single-mode
continuous-wave
(CW) laser
\cite{DelHaye2007}.
A
non-linear Schr\"odinger-like equation,
in its extended form
\cite{Lugiato1987}
with driving, damping and detuning terms,
has proved remarkably successful
in
modeling
the dynamics of these
combs
\cite{Herr2014}.
Soliton solutions to
this differential equation
can be
characterized by
the following temporal width
\cite{Herr2014}:
\begin{align}
\label{Eq_tau}
\Delta\tau
= \mathtt{C}_1 / \sqrt{ f_0 - f_{\scriptsize \mbox{p}} } ,
\end{align}
where
$f_{\scriptsize \mbox{p}}$ is the laser frequency
and
the resonator mode is centered at $f_0$.
The parameter
$\mathtt{C}_1 \equiv
\sqrt{\frac{\vert \partial_{f} n_{\mathrm{g}} \vert}
{2 n_{\mathrm{g}}}} / (2 \pi)$
depends on the group index $n_{\mathrm{g}}$,
and its dispersion.
Note that
a necessary condition for
the generation of solitons
is for the pump to be slightly red-detuned
relative to the resonance
($f_{\scriptsize \mbox{p}} < f_0$)
\cite{Coen2013-2}.
However,
due to
the thermo-optic effect,
stable comb generation
can be
obtained only by
approaching the pump frequency
from the blue side of the resonance
\cite{Carmon2004}.
The upper envelope of the soliton
spectral power density
can be approximated by
\cite{Yi2016_2}:
\begin{align}
\label{F_1}
\mathfrak{I}\left( f \right)
\propto 1 /
\cosh^{2}  \left[ \pi^2
\left( f - f_{\scriptsize \mbox{p}} + \delta \right)
\Delta\tau
\right] ,
\end{align}
where
$\delta$ is the soliton self-frequency shift
\cite{Yi2015}.
The spectral
maximum occurs at $f = f_{\scriptsize \mbox{p}} - \delta$.
Indeed, the maximum of the soliton power
does not necessarily coincide with the pump frequency.
Integrating Eq.~(\ref{F_1})
and using Eq.~(\ref{Eq_tau}),
the time-averaged soliton power
follows
\cite{Yi2015}:
\begin{align}
\mathfrak{I}_{\mathrm{sol}}
\propto
\sqrt{f_0 - f_{\scriptsize \mbox{p}} } .
\label{I_sol}
\end{align}
This one-to-one relation between soliton power and the
pump-resonance detuning has been used to servo lock the
pump laser to the cavity by
fixing the soliton power
\cite{Yi2016}.
From Eq.~(\ref{Eq_tau}),
this
procedure
also
sets
the pulse temporal width.

\subsection{Ultra-low-noise, high-power ECDL}
\label{SEC_HL}

A schematic of the
prototype
ultra-low noise (ULN)
ECDL
\cite{Morton_website_short}
is shown in Fig.~\ref{FIG_M_02}(a).
It combines
a
high performance
semiconductor
gain chip (GC) and
a polarization-maintaining (PM)
optical fiber
with an integrated
custom designed
fiber Bragg grating
(FBG)
\cite{Bird1991},
which forms one end of
the ECDL.
The GC
consists of
a multi-quantum well (MQW) active region
grown on an InP substrate.
Dielectric layers
are deposited
on one
facet
to form a high-reflectivity (HR)
coating, defining
the other end
of
the ECDL
\cite{Morton1994}.
The opposite end of the GC has both an angled
waveguide and an anti-reflection (AR) coated facet
to provide extremely low optical reflectivity,
suppressing any parasitic Fabry-Perot reflections
and supporting ULN operation.
The light emitted from the GC
is coupled
to the PM fiber via
an AR coated
lensed fiber.

\renewcommand{50mm}{55mm}
\renewcommand{3cm}{45mm}
\begin{figure}[h!]
\begin{center}
\raisebox{17mm}{(a)}
\includegraphics[width=75mm]{M_01A}
\\
\vspace{2mm}
\raisebox{43mm}{(b)}
\includegraphics[height=46mm]{SPECTR_}
\\
\vspace{2mm}
\raisebox{47mm}{(c)}
\includegraphics[height=50mm]{FN_}
\caption{
(a) Top-view schematic of the ECDL used in this work.
Pads of the butterfly package drawn in red
are for current injection to the GC,
and those in green are for temperature control.
(b) Optical spectrum measured at 525 mA and
with $T_{\mathrm{GC}} = T_{\mathrm{FBG}} =$ 20\,$^{\circ}$C.
The inset
is close-up view of the
amplified spontaneous emission.
(c) Corresponding frequency noise.
Red lines are proportional to $1/f$ and $1/f^2$,
respectively,
and the green line
indicates the white noise floor.
}
\label{FIG_M_02}
\end{center}
\end{figure}

The
ECDL
design is optimized to provide
ULN performance together
with extremely stable single-mode operation
\cite{Morton2017}.
The temperature
of the GC and of the FBG
are independently controlled,
which allows for coarse frequency tuning.
They are
packaged
in an extended butterfly
mount.
The present ECDL
does not include any moving parts
and its footprint is greatly reduced
compared to
benchtop
ECDLs
with a piezo-controlled frequency tuning mechanism
\cite{Cook2012}.
These are inherently less robust,
as their
cavity alignment
is very sensitive to
temperature fluctuations and vibrations.
A PM
optical isolator
(OZ Optics) with $\sim$60-dB isolation
is spliced
to the output fiber
to preserve the stability of this device
against
optical reflections.
Note that it has recently been suggested
that
back reflections
from the micro-resonator
could be leveraged
\cite{Taheri2017},
so that
an isolator may not be necessary
for soliton generation.
This ECDL has a threshold near 40 mA,
and an output power of 70 mW
was measured
after the isolator
with a drive current of 525 mA.
At this current level, the
signal-to-noise ratio (SNR)
is $>$70 dB, as seen
in Fig.~\ref{FIG_M_02}(b).
The inset of Fig.~\ref{FIG_M_02}(b)
is a close-up view
of
the amplified spontaneous emission (ASE).
On the red side of the
lasing peak,
the ASE is $\sim$1.5 dB higher than on the blue side.
This asymmetry
can be
understood from
the non-linear interaction
between
the optical field
and the carrier density in the QWs
\cite{Bogatov1975,Kalagara2013}.
Specifically,
the beating
between the intense lasing mode
and the ASE
modulates the (complex) refractive index,
which induces extra gain at longer wavelengths.

The frequency noise was measured
at 525 mA
with an automated system
(OEwaves, model OE4000).
As seen in Fig.~\ref{FIG_M_02}(c),
a white noise floor of
$\sim$20 Hz$^2$/Hz
is obtained for frequencies
greater than $\sim$10$^6$ Hz.
This relates \cite{Kikuchi1985} to
a Lorentzian 3-dB linewidth of
$\sim$63 Hz.
The
low
linewidth
of
this laser is
due to its long external cavity length,
high storage of photons
within this cavity \cite{Henry1983},
and operation on the long-wavelength side
of the FBG reflection
\cite{Kazarinov1987, Vahala1984}.
In the frequency range 30 Hz-2 kHz,
where the $1/f^2$-noise
and the $1/f$-noise dominate,
several distinct peaks
can be seen in Fig.~\ref{FIG_M_02}(c).
They are inherent to the noise from
the electronics of the OEwaves measurement system.
On the other hand, contributions
in the range 1-10 Hz
come
from the temperature controllers.

RIN measurements were also
performed using the OEwaves OE4000 system.
Over the entire
offset frequency
range from 10 Hz to 100 MHz,
the measured RIN is equal
to the noise floor of the equipment,
thus only providing an upper limit
to the RIN of the ECDL.
The noise floor of the OE4000 unit varies
from $-130$ dBc/Hz at 10 Hz,
down to almost $-160$ dBc/Hz at 100 kHz,
and is then
flat out to 100 MHz.
This RIN performance is comparable to
state-of-the-art low-noise ECDLs
in the C band
\cite{Loh2011},
while the linewidth demonstrated
by the device in
the present
experiments
is
$\sim$15 times narrower.

\renewcommand{50mm}{55mm}
\renewcommand{3cm}{45mm}
\begin{figure}[h!]
\begin{center}
\raisebox{47mm}{(a)}
\includegraphics[height=50mm]{fitting_}
\\
\vspace{2mm}
\raisebox{47mm}{(b)}
\includegraphics[height=50mm]{Efficiency_}
\caption{
(a)
Upper,
small-signal
current modulation
of the ECDL
for a certain value of
$I_{\mathrm{0p}}$,
and
lower,
intensity out of the MZI
fitted (in black) with Eq.~(\ref{EQ_4_all}).
(b)
Peak
frequency deviation
$\Delta f_{\max}$
obtained
for different current amplitudes
$I_{\mathrm{0p}}$.
The FM efficiency
is extracted from the slope of
the linear fit (in red).
}
\label{FIG_M_FM}
\end{center}
\end{figure}

To characterize the tuning range,
the ECDL output
was
sent
through
a fiber-based
unbalanced Mach-Zehnder interferometer
(MZI)
with a free spectral range (FSR) of $\sim$40 MHz.
By
slowly
ramping the current in the ECDL, a normalized
frequency change of $\sim$20 MHz/mA is measured,
and the widest
mode-hop-free tuning range is
2.28 GHz (or 18.3 pm).
This is much narrower than
specified for benchtop ECDLs,
but as will be shown below,
the tuning speed is at least 2 orders higher.
To characterize the possible
tuning speed,
a
time-harmonic
voltage
was
applied
to
the ECDL,
with a modulation frequency
$f_{\mathrm{m}} = \omega_{\mathrm{m}} / (2 \pi)$.
The current
through the ECDL
can be expressed as
$I(t) = I_0 + \Delta I(t)$,
with
a mean value $I_0$,
a modulation
$\Delta I(t)
= I_{\mathrm{0p}} \cos \left( \omega_{\mathrm{m}} t \right)$,
and a current peak amplitude
$I_{\mathrm{0p}}$.
This
causes
the laser frequency
to vary harmonically,
with a peak deviation $\Delta f_{\max}$.
For
$f_{\mathrm{m}} \ll \Delta f_{\max}$
and
small-signal modulation
($I_{\mathrm{0p}} \ll I_{0}$),
the
intensity out
of the MZI can be written as
\cite{Schilt2004}:
\begin{subequations}
\label{EQ_4_all}
\begin{align}
\mathfrak{I}_{\mathrm{MZI}}
=
\mathfrak{I}_{0}
+ \Delta \mathfrak{I}_{\mathrm{MZI}}
= \left( \mathfrak{I}_{0}
+ \Delta \mathfrak{I}_{\mathrm{IM}} \right)
\left(
1 + \Delta \mathfrak{I}_{\mathrm{FM}} / \mathfrak{I}_{0}
\right) ,
\end{align}
where
$\mathfrak{I}_{0}$ is the unmodulated intensity,
and:
\begin{align}
\Delta \mathfrak{I}_{\mathrm{IM}}
&= m \mathfrak{I}_{0}
\cos \left( \omega_{\mathrm{m}} t
+ \varphi_{\mathrm{IM}} \right)
\label{EQ_4b}
\\
\Delta \mathfrak{I}_{\mathrm{FM}}
&= R \mathfrak{I}_{0}
\cos \left[ 2 \pi \frac{\Delta f_{\max}}{\mathrm{FSR}}
\cos \left( \omega_{\mathrm{m}} t + \varphi_f \right)
+ \varphi_{\mathrm{FM}} \right] .
\label{EQ_4c}
\end{align}
\end{subequations}
Eq.~(\ref{EQ_4b}) represents
the spurious intensity modulation (IM)
as the laser frequency is chirped,
and $m$ is the IM index.
In contrast,
Eq.~(\ref{EQ_4c})
is used to extract the FM efficiency:
$\eta_{\mathrm{FM}}
\equiv \Delta f_{\max} / I_{\mathrm{0p}}$.
The parameter
$R$
is related to the MZI couplers,
and the $\varphi$'s are phase
constants.
In the measurement,
the amplitude $I_{\mathrm{0p}}$ is varied up to 61 mA
and the modulation frequency is set to
$f_{\mathrm{m}} = 1$ MHz.
Note that for $f_{\mathrm{m}}$ below 50 MHz,
temperature-induced FM effects
are expected to be significant
\cite{Kobayashi1982,Bowers1985}.
Specifically,
$\eta_{\mathrm{FM}}$
is expected to decrease as
$f_{\mathrm{m}}$ increases from DC to 50 MHz.
It is then usually relatively constant
up to the relaxation resonance.
An example is shown in
Fig.~\ref{FIG_M_FM}(a)
for $I_{\mathrm{0p}} = 42.7$ mA.
Fig.~\ref{FIG_M_FM}(b)
demonstrates frequency tuning speed
above
280 MHz/\textmu s.
Significantly,
a frequency deviation of 40 MHz
(required for soliton generation in a silica micro-disk),
can be generated with
a current amplitude of 40 mA.
Data
in
Fig.~\ref{FIG_M_FM}(b)
are well fitted with a straight line,
and the slope is
$\eta_{\mathrm{FM}} \approx 1.1$ MHz/mA.

\subsection{Silica micro-resonator}

The thermally-grown silica
micro-disk resonator
is fabricated on a Si substrate
with a technique reported
elsewhere \cite{Lee2012}.
Its diameter is 3 mm, which corresponds
to an FSR of
22 GHz.
The resonator
profile
is a wedge
with an angle near 30$^{\circ}$
and a thickness of 8 \textmu m.
This geometry can simultaneously
provide anomalous GVD
\cite{Li2012a,Yi2015},
high $Q$ and minimal mode crossing.
The latter effect
is known to
hinder soliton formation
\cite{Herr2014_2}.
Light is evanescently coupled
to
the micro-resonator via a tapered fiber
\cite{Cai2000,Spillane2003}.
From transmission measurements,
the resonance used
in this work
for soliton generation
has a
loaded,
full-width-at-half-depth
(FWHD)
linewidth
of
1.2 MHz,
and therefore
a total quality factor
$Q_{\scriptsize \mbox{tot}} \approx 160$ M.
The micro-resonator was under-coupled
in the measurement and had
an intrinsic value
$Q_{\scriptsize \mbox{i}} \approx 260$ M.

To overcome the thermo-optic effect so that the pump
laser frequency is red-detuned relative to the cavity
resonance for soliton generation requires rapid tuning over
frequency spans as large as 40 MHz at rates in the range of
0.1-2 MHz/\textmu s \cite{Yi2016}.
This range and rate are readily
achievable with the ECDL.
As an aside, rates in the range
of 0.3-3 MHz/\textmu s are required
in SiN combs \cite{Brasch2016_2} and
also of reach for the ECDL.

\subsection{Setup for soliton generation}

Fig.~\ref{FIG_setup}(a)-(b)
shows a schematic of the setup used
for soliton generation and
locking
with
the ECDL
described
in Section~\ref{SEC_HL}.
The ECDL is driven
by a low-noise current source
(Newport
LDX-3620B).
A signal generator
(SG, Keysight 33522B)
produces a
local oscillation (LO)
used to modulate
the voltage supplied to the ECDL,
and thereby its emission frequency.
The output of the ECDL is connected to
a polarization
controller (PC, Thorlabs FPC560),
and
the output of the tapered fiber
is connected to
a
FBG (AOS).
Its reflection port (R),
with a
0.1-nm
pass-band,
transmits the pump laser
light
which is sent to a photo-detector
(New Focus 1811, labeled PD$_1$)
and
monitored
on an oscilloscope.
In contrast,
the transmission port (T) of the FBG
acts as a notch filter and
suppresses
the pump
by 25 dB,
so that
only the light produced by
comb emission
and the ASE
are detected by
PD$_2$.
Part of this photo-current is
recorded on the oscilloscope,
and the
remainder
is sent to
a servo controller
(Vescent D2-125)
that
subtracts
a
constant
offset.
The resulting difference
serves
as the error signal
for a feedback-locking loop
\cite{Yi2015,Yi2016}
essentially based on Eq.~(\ref{I_sol}).

Once the pump frequency reaches
the red side of the resonance
(prerequisite for soliton generation,
see
Section~\ref{SEC_2A}),
the SG sends a digital signal
to
the servo
to engage
the locking loop.
The servo controller
includes
an op-amp
integrator.
To establish a stable loop,
it
forces a null
at its input by integrating the error signal
within a
10-kHz
bandwidth.
Its output is
a correction voltage
signal
that is fed back to
the ECDL current source,
combined with the LO.
By
adjusting
the offset value,
states with different soliton numbers
can be generated
\cite{Yi2016}.
As shown in the next Section,
the tuning and locking parameters
were
optimized to guarantee
the generation
and stabilization
of a \emph{single}-soliton state.

\begin{figure}[h!]
\begin{center}
\begin{minipage}[b]{15cm}
\centering
\includegraphics[width=80mm]
{Full_}
\caption{
(a)-(b)
Schematic of the setup
used in this work
for soliton generation and
locking
without an OA.
Black and green lines
represent optical and
electrical
paths,
respectively.
A cross-section schematic of the micro-disk
resonator and its wedge profile is shown in (b),
with intensity pattern
simulated for the fundamental mode.
Alternatively,
for soliton generation with
a piezo-controlled ECDL,
part (a)
is to be replaced
by part (c).
}
\label{FIG_setup}
\end{minipage}
\end{center}
\end{figure}

The conventional methods for generation of solitons using
piezo-controlled ECDLs can be understood
by replacing
part (a)
with part (c) in Fig.~\ref{FIG_setup}.
Specifically,
two common methods
have been
developed to overcome the thermal
instability
and generate solitons.
They are based on
abrupt changes of either the pump power
\cite{Yi2015, Yi2016, Brasch2016_1, Brasch2016_2}
or the pump frequency \cite{Stone2017arXiv},
and referred to as
``power kicking" or ``frequency kicking",
respectively.
These changes occur
over
timescales that are
faster than the thermal
time constant of the micro-resonator.
Referring to
Fig.~\ref{FIG_setup}(c),
the pump is
typically
provided
by a benchtop ECDL
with a piezo-controlled frequency.
The servo feedback
is applied to the piezo-controller,
and not to the current source like
in
the
procedure described
above.
For power kicking,
the modulator
in Fig.~\ref{FIG_setup}(c)
consists of a
combination of
acousto-optic and electro-optic modulators
(with modulation frequencies up to $\sim$100 MHz).
In contrast, the
frequency kicking protocol
is achieved
by manipulating the single sideband
from a quadrature phase shift keying (QPSK) modulator driven
by a fast-tuning voltage-controlled oscillator (VCO).
In both protocols,
the modulators
are controlled by an extra SG
and
have optical insertion loss of a few dB.
An
optical amplifier (OA)
is thus inevitably required
to boost
the pump power.
A band-pass filter (BPF)
is also preferred
to suppress the undesired ASE
from the OA.
Note that AOMs and QPSKs
require a specific polarization
for optimal operation,
so that an extra PC is
also needed.
These instruments dramatically increase the footprint,
power consumption,
and cost
of the micro-comb system.
In this work, we are able to eliminate these bulky components
by
implementing a frequency-kicking protocol
directly with an ECDL.

\section{Results and discussion}

\newcommand{171mm}{170mm}

\renewcommand{3cm}{45mm}

\begin{figure*}[h!]
\begin{center}
\includegraphics[width=171mm]{Soliton_Locking_}
\caption{
Frequency
tuning leading to
the generation of a single soliton,
and its subsequent locking.
The green and the red curves
are the transmitted power
from the pump and from the comb, respectively.
The inset a close-up view of
the thermal drift and
the soliton locking.
The black line is a moving average,
and the horizontal scale is logarithmic.
}
\label{FIG_soliton_locking}
\end{center}
\end{figure*}

The ECDL is
driven
at 350 mA where its output power is 42 mW,
and the power
coupled into
the tapered fiber
is 34 mW.
Fig.~\ref{FIG_soliton_locking}
shows the evolution of
the pump power $P_1$
measured after the micro-resonator,
as the pump frequency $f_{\scriptsize \mbox{p}}$
is linearly tuned
from the blue to the red side of
the resonance.
Also shown is the
generated
comb power $P_2$.
Up to $t_1$,
$f_{\scriptsize \mbox{p}}$
is so far from resonance
that it does not couple significantly to the micro-disk
and
$P_1$
is relatively constant.
This situation changes
near $t_1$
where
a
reduction of
$P_1$
is observed
as more light couples to the micro-disk.
By
$t_{2.1}$,
the power coupled to the resonator
eventually reaches
threshold for parametric oscillation,
leading to
the onset of a primary comb
and a sharp increase of $P_2$.
From $t_{2.1}$ to
$t_{3}$,
$P_1$
further decreases
as more pump power is
coupled to the resonator,
and between $t_{2.2}$
and $t_{3}$,
PD$_2$ is saturated.
Heating of the
resonator induces the overall triangular-shaped profile of
the pump power transmission
\cite{Carmon2004}.
Upon arriving at the edge of the blue-detuned side,
the pump frequency
is kicked a few MHz to the red-detuned side for a few \textmu s,
which induces the soliton state
in
the resonator.
Finally, at $t_{4}$,
\textit{i.e.} a few hundred \textmu s after the scanning,
the servo loop is engaged to adjust $f_{\scriptsize \mbox{p}}$,
so as to
lock the
soliton state indefinitely against thermal drifting.
Variations in
$P_2$
are further suppressed,
suggesting from Eq.~(\ref{I_sol})
that the
pump-resonance
detuning is
indeed stabilized.

\renewcommand{3cm}{45mm}
\begin{figure}[h!]
\begin{center}
\begin{minipage}[b]{15cm}
\centering
\raisebox{49mm}{(a)}
\includegraphics[height=52mm]{soliton1_}
\\
\vspace{2mm}
\raisebox{37.5mm}{(b)}
\includegraphics[height=40.5mm]{RF_}
\raisebox{37.5mm}{(c)}
\includegraphics[height=40.5mm]{PhN_}
\caption{
(a) Optical spectrum measured
in the single-soliton regime.
The red line is a fit according to Eq.~(\ref{F_1}).
(b)
Corresponding
RF electrical spectrum,
with $f_{\mathrm{c}} \approx 22.02$ GHz,
fitted with a Lorentzian (in red).
The bottom green line indicates the noise level.
(c)
Phase noise
spectral density.
The red line is proportional to $1/f^3$
and black squares indicate the instrument sensitivity.
}
\label{FIG_soliton}
\end{minipage}
\end{center}
\end{figure}

To confirm soliton generation,
the
light output from the T-port of the FBG
in Fig.~\ref{FIG_setup}(b)
is sent to an optical spectrum analyser
(OSA, Yokogawa AQ6370).
It is seen in
Fig.~\ref{FIG_soliton}(a)
that the spectral envelope
is relatively smooth,
without significant
mode crossings.
However, an intensity spike
with an amplitude 17 dB above the ASE level
is centered
near 1529 nm
(corresponding to
a relative mode number $\mu = 120$).
This spectral feature
is a signature of a dispersive wave
\cite{Yi2017},
which can occur from the
interaction
of the soliton with
different transverse modes in the micro-resonator.
Eq.~(\ref{F_1})
provides an excellent fit for
the
remainder
of this spectrum,
and
leads
$\Delta\tau \approx 196$ fs,
and $\delta \approx 60$ GHz.
The obtained temporal width
is within
the range of values (125-215 fs) reported in
\cite{Yi2016_2}.
The
value obtained for $\delta$
indicates a red-shifting
of the soliton
from the pump
by nearly 3 FSRs.
This
self-frequency shift
is attributed to
the spectral recoil
caused by the dispersive wave,
and possible
stimulated Raman scattering
\cite{Yi2016_2, Karpov2016, Yi2017}.
Here,
a single-soliton state is
directly
generated,
in contrast to \cite{Guo2017}.
No step-like features are visible
in the oscilloscope traces
of Fig.~\ref{FIG_soliton_locking},
as these would occur
at the transition
between different soliton states.

The soliton is further confirmed
by
assessing
the
coherence
of the generated comb
\cite{Erkintalo2014}.
Indeed,
a low-noise and narrow RF signal is
known to be
a necessary signature of stable soliton formation
\cite{Herr2014}.
The OSA
is replaced by a PD
(50-GHz bandwidth)
connected to
an electrical spectrum analyser (ESA,
Rohde \& Schwarz, FSUP26)
that
also measures
the phase noise.
Fig.~\ref{FIG_soliton}(b)
shows the radio-frequency (RF) electrical spectrum
recorded
with a 100-Hz resolution bandwidth.
The
carrier frequency
$f_{\mathrm{c}} \approx 22$ GHz
corresponds to the
comb repetition frequency,
\textit{i.e.}
to the beating between neighboring comb lines.
These data
have an SNR
$>$80 dB.
They are fitted with a Lorentzian,
and the 3-dB linewidth is $\sim$25 Hz.
Note that
the central part is Gaussian
with a 3-dB linewidth of $\sim$1.1 kHz.

Finally,
Fig.~\ref{FIG_soliton}(c)
shows the
phase noise spectral
density
of the repetition beat note
plotted versus the offset frequency $f$.
Between 10 Hz and 100 kHz,
\textit{i.e.} over 4 decades,
the phase noise is
approximately
proportional to $1/f^3$.
This
is believed to be contributed by a
combination of pump laser RIN and frequency noise
\cite{Lee2013, Liang2015}.
The detected phase noise is smaller than $-100$ dBc/Hz
for offset frequencies higher than
28 kHz,
which
is comparable to the data
from a previous report
\cite{Yi2017},
where the silica micro-resonators
were pumped
with a benchtop ECDL.
Beyond
1 MHz,
the detected phase noise reaches
its
floor value
of $-129$ dBc/Hz,
attributed to the PD shot noise
\cite{Liang2015}.

\section{Conclusion}

Soliton generation was demonstrated in a micro-resonator
with an
ultra-low-noise diode laser.
A single-soliton state was successfully stabilized
by locking the frequency of the pump laser
to the power level of the
soliton comb.
To our knowledge, it is the first time
a temporal soliton
is realized
in a
chip-based
micro-resonator
without an optical amplifier.
This demonstration reduces
the complexity and cost of soliton experiments,
representing a significant
step towards
soliton comb systems fully
integrated on a chip.

\section*{Funding Information}
This research is supported by a DARPA MTO
DODOS contract (HR0011-15-C-055).
The views and conclusions contained in this document are those of the authors and should not be interpreted as representing official policies of the Defense Advanced Research Projects Agency or the U.S. Government.
N.V. acknowledges support from the Swiss National Science Foundation (SNSF).

\section*{Acknowledgments}

The authors thank
Rui-Lin Chao,
Sarat Chandra Gundavarapu,
and Xinbai Li
for experimental assistance.

\color{black}

\tolerance=1
\emergencystretch=\maxdimen
\hyphenpenalty=10000
\hbadness=10000

\bibliographystyle{lpr}