\section{Anharmonic energy level structure of artificial atoms}

The cavity, which is designed for phase-matching condition that enables
non-degenerate $\chi^{(2)}$ interaction $\omega_{a}+\omega_{b}\rightarrow\omega_{c}$,
has a Hamiltonian of
\begin{eqnarray}
H & = & \omega_{\mathrm{a}}a^{\dagger}a+\omega_{\mathrm{b}}b^{\dagger}b+\omega_{\mathrm{c}}c^{\dagger}c+g_{\mathrm{nd}}\left(a^{\dagger}b^{\dagger}c+abc^{\dagger}\right),
\end{eqnarray}
where $a,b,c$ represents the bosonic operator for modes, $g_{\mathrm{nd}}$
is the nonlinear coupling strength. In the strong coupling regime
$g_{\mathrm{nd}}\gtrsim\kappa_{\mathrm{a}},\kappa_{\mathrm{b}},\kappa_{\mathrm{c}}$,
the system shows strong anhamonicity. For example, the state $|1_{\mathrm{a}}1_{\mathrm{b}}0_{\mathrm{c}}\rangle$
will couple to $|0_{\mathrm{a}}0_{\mathrm{b}}1_{\mathrm{c}}\rangle$
with a strength of $g_{\mathrm{nd}}$, and these two states will be
hybridized with each other, leading to a splitting of $2g_{\mathrm{nd}}$.
In this nonlinear cavity system, the eigenstates and eigenvalues in
each $N$-photon subspace could be solved independently. The energy
conversation requires
\begin{eqnarray}
N & = & N_{a}+N_{b}+2N_{c},
\end{eqnarray}
where $N_{i}$ is the number of excitation in mode $i$ ($i\in\left\{ \mathrm{a},\mathrm{b},\mathrm{c}\right\} $).
Therefore, the basis $\left\{ |\Psi_{N,m,n}\rangle=|m-1,n-1,N-m-n+2\rangle\right\} _{N}$
spans the $N$-photon subspace $H^{\left(N\right)}$, where $m,n\leq N$
are integers with $m+n-2<N$. The Hamiltonian is written in the matrix
form as
\begin{eqnarray}
H^{\left(N\right)} & = & \sum H_{m,n;m',n'}^{\left(N\right)}|\Psi_{N,m,n}\rangle\langle\Psi_{N,m',n'}|,
\end{eqnarray}
with the matrix element
\begin{eqnarray}
H_{m,n;m',n'}^{\left(N\right)} & = & \langle\Psi_{N,m',n'}|H|\Psi_{N,m,n}\rangle.
\end{eqnarray}
The eigenstates and eigenenergies are solved by diagonalizing the
matrix. Specially, we only consider the interaction between two photons,
thus truncate the dimension of the Fock state space to $2$ is enough.
The Hamiltonian in $2$-photon excitation is spanned by the eigenstate
$|\tilde{j}\rangle\in\left\{ \begin{array}{cc}
|1,1,0\rangle, & |0,0,1\rangle\end{array}\right\} $. For perfect phase matching, an appropriate rotating frame can be
chosen to transform $H$ to $H=g\left(a^{\dagger}b^{\dagger}c+abc^{\dagger}\right),$
\begin{eqnarray}
H_{ij}^{2} & = & \langle\tilde{j}|H|\tilde{i}\rangle_{2}\\
 & = & g_{\mathrm{nd}}\left(\delta(i-1,j)+\delta(i,j-1)\right),
\end{eqnarray}
where $\delta(x,y)$ is the Kronecker delta function. The energy levels
are obtained by diagolizing the matrix of the Hamiltonian. Therefore,
the degeneracy of the two-excitation states is removed by the nonlinear
interaction with energy levels split by $2g_{nd}$. If the cavity
is initially prepared in the state $|1_{a}0_{b}0_{c}\rangle$, the
photon resonant with mode $b$ can no longer efficiently induce the
transition from $|1_{a}0_{b}0_{c}\rangle$ to $|0_{a}0_{b}1_{c}\rangle$,
which is known as the photon-blockade effect. The energy levels of
second-harmonic generation $g_{\mathrm{d}}\left(a^{\dagger2}c+a^{2}c^{\dagger}\right)$
can be obtained following a similar calculation.

\section{Two-photon scattering by the artificial atom}

In this section, we investigate the two-photon transport through a
waveguide side-coupled to a cavity supporting $\chi^{(2)}$ interaction.
Figure 2(a) in the main text schematically illustrates a doubly-resonant
microcavity coupled with a bus waveguide. The two resonant modes in
the cavity are the fundamental (FM) mode $a$ and the second-harmonic
(SH) mode $c$ with frequency $\omega_{i}$ ($i\in\{a,c\}$). The
waveguide supports continuum states both in the fundamental and SH
frequency bands. Inside the cavity, modes $a$ and $c$ couple with
each other via the process of SH generation. The Hamiltonian of whole
system includes, the continuum states in the waveguide
\begin{eqnarray}
H_{W} & = & \int dxf^{\dagger}(x)\frac{1}{i}\frac{d}{dx}f(x)+\int dxs^{\dagger}(x)\frac{1}{i}\frac{d}{dx}s(x),
\end{eqnarray}
the localized states in the cavity
\begin{eqnarray}
H_{C} & = & \left(\omega_{a}-i\kappa_{a,0}\right)a^{\dagger}a+\left(\omega_{c}-i\kappa_{c,0}\right)c^{\dagger}c,
\end{eqnarray}
the nonlinear coupling inside the cavity
\begin{eqnarray}
H_{nl} & = & g_{\mathrm{d}}\left(a^{\dagger2}c+a^{2}c^{\dagger}\right),
\end{eqnarray}
the linear coupling between the waveguide and cavity
\begin{eqnarray}
H_{l} & = & V_{a}\int dx\delta(x)[f^{\dagger}(x)a+f(x)a^{\dagger}]+V_{c}\int dx\delta(x)[s^{\dagger}(x)c+s(x)c^{\dagger}],
\end{eqnarray}
where $x$ is coordinate along the waveguide, $f(x)$ ($s(x)$) and
$f^{\dagger}(x)$ ($s^{\dagger}(x)$) are the annihilation and creation
operators at position $x$ in the waveguide, $g_{\mathrm{d}}$ is
the nonlinear coupling strength of SH generation, $V_{a}$ ($V_{c}$)
is the coupling strength between waveguide and cavity in the fundamental
(SH) frequency band, $\kappa_{a,0}$ ($\kappa_{c,0}$) is the intrinsic
decay rates of the cavity mode. Since the coupling near the narrow
resonant frequency window is considered, the dispersions of the coupling
strengths $g_{\mathrm{d}}$, $V_{a}$ and $V_{c}$ are neglected.
The input state is initially prepared in a two-photon state at the
fundamental frequencies. A general form of two-excitation state can
be written as
\begin{eqnarray}
|E_{2}\rangle & = & \int dx_{1}dx_{2}g(x_{1},x_{2})f^{\dagger}(x_{1})f^{\dagger}(x_{2})|0\rangle+\int dxm(x)f^{\dagger}(x)|0\rangle+\nonumber \\
 &  & \int dxh(x)s^{\dagger}(x)|0\rangle+f_{a}\frac{1}{\sqrt{2}}a^{\dagger}a^{\dagger}|0\rangle+f_{c}c^{\dagger}|0\rangle.
\end{eqnarray}
Following the Schrodinger equation $H|E_{2}\rangle=E_{2}|E_{2}\rangle$
($E_{2}=k_{1}+k_{2}$, $k_{i}$ is the frequency of the $i$-th photon),
one gets
\begin{eqnarray}
\left(-i\frac{\partial}{\partial x_{1}}-i\frac{\partial}{\partial x_{2}}-E_{2}\right)g(x_{1},x_{2})+\frac{V_{a}}{2}[\delta(x_{1})m(x_{2})+\delta(x_{2})m(x_{1})] & = & 0,\label{eq:two-ea-1}\\
\left(-i\frac{d}{dx}+\omega_{a}-i\kappa_{a,0}-E_{2}\right)m(x)+V_{a}[g(x,0)+g(0,x)]+\sqrt{2}V_{a}f_{a}\delta(x) & = & 0,\\
\left(-i\frac{d}{dx}-E_{2}\right)h(x)+V_{c}f_{c}\delta(x) & = & 0,\label{eq:two-ea}\\
\left(2\omega_{a}-2i\kappa_{a,0}-E_{2}\right)f_{a}+\sqrt{2}g_{\mathrm{d}}f_{c}+\sqrt{2}V_{a}m(0) & = & 0,\\
\left(\omega_{c}-i\kappa_{c,0}-E_{2}\right)f_{c}+\sqrt{2}g_{\mathrm{d}}f_{a}+V_{c}h(0) & = & 0.\label{eq:two-ea-5}
\end{eqnarray}
The discontinuous functions fulfills
\begin{eqnarray}
g(x,0) & = & g(0,x),\\
g(x,0) & = & [g(x,0^{+})+g(x,0^{-})]/2,\\
m(0) & = & [m(0^{+})+m(0^{-})]/2,\\
h(0) & = & [h(0^{+})+h(0^{-})]/2,\\
\frac{d}{dx}m(x)|_{x=0} & = & [m(0^{+})-m(0^{-})]\delta(x).
\end{eqnarray}
Then Eqs.(\ref{eq:two-ea-1})-(\ref{eq:two-ea-5}) are equivalent
to
\begin{eqnarray}
\left(-i\frac{\partial}{\partial x_{1}}-i\frac{\partial}{\partial x_{2}}-E_{2}\right)g(x_{1},x_{2}) & = & 0,\\
\frac{2i}{V_{a}}[g(0^{+},x)-g(0^{-},x)] & = & m(x),\label{eq:2}\\
\left(-i\frac{d}{dx}+\omega_{a}-i\kappa_{a,0}-E_{2}\right)m(x)+V_{a}[g(x,0^{+})+g(0^{-},x)] & = & 0,\label{eq:3}\\
\frac{i}{\sqrt{2}V_{a}}[m(0^{+})-m(0^{-})] & = & f_{a},\\
\left(-i\frac{d}{dx}-E_{2}\right)h(x) & = & 0,\\
\frac{i}{V_{c}}[h(0^{+})-h(0^{-})] & = & f_{c},\\
\left(2\omega_{a}-2i\kappa_{a,0}-E_{2}\right)f_{a}+\sqrt{2}g_{\mathrm{d}}f_{c}+\sqrt{2}V_{a}m(0) & = & 0,\\
\left(\omega_{c}-i\kappa_{c,0}-E_{2}\right)f_{c}+\sqrt{2}g_{\mathrm{d}}f_{a}+V_{c}h(0) & = & 0.
\end{eqnarray}
Combine Eqs.$\,$(\ref{eq:2}) and (\ref{eq:3}), we get
\begin{eqnarray}
\left(-i\frac{d}{dx}+\omega_{a}-i\kappa_{a,0}-i\kappa_{a,1}-E_{2}\right)g(0^{+},x) & = & \left(-i\frac{d}{dx}+\omega_{a}-i\kappa_{a,0}+i\kappa_{a,1}-E_{2}\right)g(0^{-},x).\label{eq:region}
\end{eqnarray}
Adopting the method in Ref.$\,$\cite{Zheng2013}, $g(x_{1},x_{2})$
is divided to three regions: $g_{1}(x_{1},x_{2})$ for $x_{1}\leq x_{2}<0$,
$g_{2}(x_{1},x_{2})$ for $x_{1}<0<x_{2}$ and $g_{3}(x_{1},x_{2})$
for $0<x_{1}\leq x_{2}$. For two-photon monochromatic wave input
state
\begin{eqnarray}
g_{1}(x_{1},x_{2}) & = & \frac{1}{2\text{\ensuremath{\sqrt{2}\pi}}}\left(e^{ik_{1}x_{1}}e^{ik_{2}x_{2}}+e^{ik_{2}x_{1}}e^{ik_{1}x_{2}}\right),
\end{eqnarray}
$g(x_{1},x_{2})$ in the other two regions are derived from Eq.$\,$(\ref{eq:region})
as
\begin{eqnarray}
g_{2}(x_{1},x_{2}) & = & \frac{1}{2\text{\ensuremath{\sqrt{2}\pi}}}\left(t_{k_{2}}e^{i(k_{1}x_{1}+k_{2}x_{2})}+t_{k_{1}}e^{i(k_{2}x_{1}+k_{1}x_{2})}\right),\\
g_{3}(x_{1},x_{2}) & = & \frac{1}{2\text{\ensuremath{\sqrt{2}\pi}}}t_{k_{1}}t_{k_{2}}\left(e^{i(k_{1}x_{1}+k_{2}x_{2})}+e^{i(k_{2}x_{1}+k_{1}x_{2})}\right)+Be^{-i(\omega_{a}-i\kappa_{a,0}-i\kappa_{a,1})(x_{2}-x_{1})}e^{i(k_{1}+k_{2})x_{2}},
\end{eqnarray}
with
\begin{eqnarray}
t_{k} & = & \frac{k-\omega_{a}+i\kappa_{a,0}-i\kappa_{a,1}}{k-\omega_{a}+i\kappa_{a,0}+i\kappa_{a,1}}.
\end{eqnarray}
Eliminating $f_{a}$ and $f_{c}$ ,
\begin{eqnarray}
\frac{2i}{V_{a}}[\alpha_{a}^{-}g(0^{+},0^{+})-(\alpha_{a}^{-}+\alpha_{a}^{+})g(0^{-},0^{+})+\alpha_{a}^{+}g(0^{-},0^{-})]+2g_{\mathrm{d}}\frac{V_{a}}{V_{c}}[h(0^{+})-h(0^{-})] & = & 0,\\
\alpha_{b}^{-}h(0^{+})-\alpha_{b}^{+}h(0^{-})+\frac{2ig_{\mathrm{d}}V_{c}}{V_{a}^{2}}[g(0^{+},0^{+})-2g(0^{-},0^{+})+g(0^{-},0^{-})] & = & 0,
\end{eqnarray}
with $\alpha_{a}^{\pm}=2\left(\omega_{a}-i\kappa_{a,0}\pm i\kappa_{a,1}\right)-k_{1}-k_{2}$,
$\alpha_{c}^{\pm}=\left(\omega_{c}-i\kappa_{c,0}\pm i\kappa_{c,1}\right)-k_{1}-k_{2}$,
$\kappa_{i,1}=V_{i}^{2}/2$. Using the values at the discontinuous
points, 
\begin{eqnarray}
g(0^{+},0^{+})=g_{3}(0,0) & = & \frac{1}{\sqrt{2}\pi}t_{k_{1}}t_{k_{2}}+B,\\
g(0^{-},0^{+})=g_{2}(0,0) & = & \frac{1}{2\sqrt{2}\pi}(t_{k_{1}}+t_{k_{2}}),\\
g(0^{-},0^{-})=g_{1}(0,0) & = & \frac{1}{\sqrt{2}\pi},\\
h(0^{+}) & = & \frac{1}{\sqrt{2\pi}}C,\\
h(0^{-}) & = & 0,
\end{eqnarray}
one finally obtains
\begin{eqnarray}
C & = & \frac{ig_{\mathrm{d}}\sqrt{2\kappa_{c,1}}\left(\alpha_{a}^{+}-\alpha_{a}^{-}\right)\left(2-t_{k_{1}}-t_{k_{2}}\right)}{2\sqrt{\pi}\kappa_{a,1}\left(\alpha_{a}^{-}\alpha_{c}^{-}-2g_{d}^{2}\right)}\nonumber \\
 & = & -\frac{2g_{\mathrm{d}}\sqrt{2\kappa_{c,1}}\left(2-t_{k_{1}}-t_{k_{2}}\right)}{\sqrt{\pi}\left(\alpha_{a}^{-}\alpha_{c}^{-}-2g_{\mathrm{d}}^{2}\right)},\\
B & = & \frac{4g_{\mathrm{d}}^{2}\left(t_{k_{1}}-1\right)\left(t_{k_{2}}-1\right)+\alpha_{a}^{+}\alpha_{c}^{-}\left(t_{k_{1}}+t_{k_{2}}-2\right)+\alpha_{a}^{-}\alpha_{c}^{-}\left(t_{k_{1}}+t_{k_{2}}-2t_{k_{1}}t_{k_{2}}\right)}{2\sqrt{2}\pi\left(\alpha_{a}^{-}\alpha_{c}^{-}-2g_{\mathrm{d}}^{2}\right)}\nonumber \\
 & = & \frac{\sqrt{2}g_{\mathrm{d}}^{2}\left(1-t_{k_{1}}\right)\left(1-t_{k_{2}}\right)}{\pi\left(\alpha_{a}^{-}\alpha_{c}^{-}-2g_{\mathrm{d}}^{2}\right)}.
\end{eqnarray}
The expressions can be simplified as
\begin{eqnarray}
C & = & -\frac{2g_{\mathrm{d}}\sqrt{2\kappa_{c,1}}}{\sqrt{\pi}\left[2\omega_{a}-k_{1}-k_{2}-2i\kappa_{a,tot}\right]\left[\omega_{c}-k_{1}-k_{2}-i\kappa_{c,tot}\right]-2g_{d}^{2}}\left(\frac{2i\kappa_{a,1}}{k_{1}-\omega_{a}+i\kappa_{a,tot}}+\frac{2i\kappa_{a,1}}{k_{2}-\omega_{a}+i\kappa_{a,tot}}\right),\\
B & = & \frac{\sqrt{2}g_{d}^{2}}{\pi\left[2\omega_{a}-k_{1}-k_{2}-2i\kappa_{a,tot}\right]\left[\omega_{c}-k_{1}-k_{2}-i\kappa_{c,tot}\right]-2g_{\mathrm{d}}^{2}}\frac{2i\kappa_{a,1}}{k_{1}-\omega_{a}+i\kappa_{a,tot}}\frac{2i\kappa_{a,1}}{k_{2}-\omega_{a}+i\kappa_{a,tot}}.
\end{eqnarray}
where $\kappa_{i,tot}=\kappa_{i,0}+\kappa_{i,1}$. The wave-functions
of two-photon input and output states are computed by the Lippmann-Schwinger
equations$\:$\cite{lipp}, 
\begin{eqnarray}
\phi_{in}(x_{1},x_{2}) & = & \langle x_{1},x_{2}|E_{2}\rangle-\langle x_{1},x_{2}|\frac{1}{k_{1}+k_{2}-H_{0}+i0^{+}}H_{int}|E_{2}\rangle\nonumber \\
 & = & g_{1}(x_{1},x_{2})\theta(x_{2}-x_{1})+g_{1}(x_{2},x_{1})\theta(x_{1}-x_{2}),\\
\phi_{out}(x_{1},x_{2}) & = & \langle x_{1},x_{2}|E_{2}\rangle-\langle x_{1},x_{2}|\frac{1}{k_{1}+k_{2}-H_{0}-i0^{+}}H_{int}|E_{2}\rangle\nonumber \\
 & = & g_{3}(x_{1},x_{2})\theta(x_{2}-x_{1})+g_{3}(x_{2},x_{1})\theta(x_{1}-x_{2}).
\end{eqnarray}
Then the two-photon scattering matrix is obtained by calculating the
wave-function overlap between the output state and the two photon
plain wave state,
\begin{eqnarray}
S(p_{1},p_{2};k_{1},k_{2}) & = & \int dx_{1}dx_{2}\frac{1}{2\sqrt{2}\pi}[e^{-i(p_{1}x_{1}+p_{2}x_{2})}+e^{-i(p_{2}x_{1}+p_{1}x_{2})}]\phi_{out}(x_{1},x_{2}).
\end{eqnarray}
The transmission plane wave term is
\begin{eqnarray}
S_{plane} & = & \frac{1}{8\pi^{2}}\int dx_{1}dx_{2}[e^{-i(p_{1}-k_{1})x_{1}}e^{-i(p_{2}-k_{2})x_{2}}+e^{-i(p_{1}-k_{2})x_{1}}e^{-i(p_{2}-k_{1})x_{2}}\nonumber \\
 &  & +e^{-i(p_{1}-k_{2})x_{2}}e^{-i(p_{2}-k_{1})x_{1}}+e^{-i(p_{2}-k_{2})x_{1}}e^{-i(p_{1}-k_{1})x_{2}}]\nonumber \\
 & = & t_{k_{1}}t_{k_{2}}[\delta(p_{1}-k_{1})\delta(p_{2}-k_{2})+\delta(p_{1}-k_{2})\delta(p_{2}-k_{1})]\nonumber \\
 & = & S_{p_{1}k_{1}}S_{p_{2}k_{2}}+S_{p_{1}k_{2}}S_{p_{2}k_{1}}.
\end{eqnarray}
The bound state term is
\begin{eqnarray}
S_{bound} & = & \frac{B}{2\sqrt{2}\pi}\int dx_{1}dx_{2}[e^{-i(p_{1}x_{1}+p_{2}x_{2})}+e^{-i(p_{2}x_{1}+p_{1}x_{2})}]e^{-i(\omega_{a}-i\kappa_{a,0}-i\kappa_{a,1})(x_{2}-x_{1})}e^{i(k_{1}+k_{2})x_{2}}\theta(x_{2}-x_{1})+\nonumber \\
 &  & \frac{B}{2\sqrt{2}\pi}\int dx_{1}dx_{2}[e^{-i(p_{1}x_{1}+p_{2}x_{2})}+e^{-i(p_{2}x_{1}+p_{1}x_{2})}]e^{-i(\omega_{a}-i\kappa_{a,0}-i\kappa_{a,1})(x_{1}-x_{2})}e^{i(k_{1}+k_{2})x_{1}}\theta(x_{1}-x_{2}).
\end{eqnarray}
The integral of the first term is simplified by replacing $x_{2}$
by $x_{3}+x_{1}$ and similar for the second term. We get
\begin{eqnarray}
S_{bound,1} & = & \frac{B}{2\sqrt{2}\pi}\int dx_{1}dx_{3}e^{-(\kappa_{a,0}+\kappa_{a,1})x_{3}}[e^{-i(p_{1}+p_{2}-k_{1}-k_{2})x_{1}}e^{-i(p_{2}+\omega_{a}-k_{1}-k_{2})x_{3}}\nonumber \\
 &  & +e^{-i(p_{1}+p_{2}-k_{1}-k_{2})x_{1}}e^{-i(p_{1}+\omega_{a}-k_{1}-k_{2})x_{3}}]\theta(x_{3})\nonumber \\
 & = & \frac{B}{2\sqrt{2}\pi}\times2\pi[\frac{1}{(\kappa_{a,0}+\kappa_{a,1})+i(p_{2}+\omega_{a}-k_{1}-k_{2})}\delta(p_{1}+p_{2}-k_{1}-k_{2})\nonumber \\
 &  & +\frac{1}{(\kappa_{a,0}+\kappa_{a,1})+i(p_{1}+\omega_{a}-k_{1}-k_{2})}\delta(p_{1}+p_{2}-k_{1}-k_{2})]\nonumber \\
 & = & \frac{B}{\sqrt{2}}[\frac{i}{p_{1}-\omega_{a}+i(\kappa_{a,0}+\kappa_{a,1})}+\frac{i}{p_{2}-\omega_{a}+i(\kappa_{a,0}+\kappa_{a,1})}]\delta(p_{1}+p_{2}-k_{1}-k_{2}).
\end{eqnarray}
Due to the symmetry of the wavefunction
\begin{eqnarray}
S_{bound} & = & 2S_{bound,1}\nonumber \\
 & = & \sqrt{2}B[\frac{i}{p_{1}-\omega_{a}+i(\kappa_{a,0}+\kappa_{a,1})}+\frac{i}{p_{2}-\omega_{a}+i(\kappa_{a,0}+\kappa_{a,1})}]\delta(p_{1}+p_{2}-k_{1}-k_{2})\nonumber \\
 & = & \sqrt{2}B[\frac{i}{p_{1}-\omega_{a}+i(\kappa_{a,0}+\kappa_{a,1})}+\frac{i}{p_{2}-\omega_{a}+i(\kappa_{a,0}+\kappa_{a,1})}]\delta(p_{1}+p_{2}-k_{1}-k_{2})\nonumber \\
 & = & \frac{B}{\sqrt{2}\kappa_{a,1}}(2-t_{p_{1}}-t_{p_{2}})\delta(p_{1}+p_{2}-k_{1}-k_{2}).
\end{eqnarray}
The total scattering matrix is 
\begin{eqnarray}
S(p_{1},p_{2};k_{1},k_{2}) & = & t_{k_{1}}t_{k_{2}}[\delta(p_{1}-k_{1})\delta(p_{2}-k_{2})+\delta(p_{1}-k_{2})\delta(p_{2}-k_{1})]\\
 &  & +\frac{B}{\sqrt{2}\kappa_{a,1}}(2-t_{p_{1}}-t_{p_{2}})\delta(p_{1}+p_{2}-k_{1}-k_{2}).
\end{eqnarray}
According to the scattering matrix, one can compute the output spectrum
of the two-photon state as
\begin{eqnarray}
\alpha_{out}(p_{1},p_{2}) & = & \int\int\alpha_{in}(k_{1},k_{2})S(p_{1},p_{2};k_{1},k_{2})dk_{1}dk_{2},
\end{eqnarray}
where $\alpha_{in}(k_{1},k_{2})$ is the input spectrum.

From the frequency conversion point of view, one might care about
the conversion efficiency of the two-photon second-harmonic generation.
The two-photon to one-photon conversion matrix is
\begin{eqnarray}
S(p;k_{1},k_{2}) & = & \int dxe^{-ipx}\frac{1}{\sqrt{2\pi}}C(k_{1}+k_{2})e^{i(k_{1}+k_{x})x}\nonumber \\
 & = & \frac{1}{\sqrt{2\pi}}C(k_{1}+k_{2})\int dxe^{-i(p-k_{1}-k_{2})x}\nonumber \\
 & = & C(k_{1}+k_{2})\delta(p-k_{1}-k_{2}).
\end{eqnarray}
The output spectrum of the second-harmonic photon is 
\begin{eqnarray}
\beta_{out}(p) & = & \int\int\alpha_{in}(k_{1},k_{2})S(p;k_{1},k_{2})dk_{1}dk_{2}\nonumber \\
 &  & \int\int\alpha_{in}(k_{1},k_{2})C(k_{1}+k_{2})\delta(p-k_{1}-k_{2})dk_{1}dk_{2}.
\end{eqnarray}
Then, we get the conversion efficiency
\begin{eqnarray}
\eta & = & \frac{\int|\beta_{out}(p)|^{2}dp}{2\int\int|\alpha_{in}(k_{1},k_{2})|^{2}dk_{1}dk_{2}}.
\end{eqnarray}


\section{Control phase gate via photonic molecule}

\subsection{Coupled microresonator}

In our scheme, an active antenna is used to couple the traveling photons
and the artificial atom, as shown by Fig.$\:$3 in the main text.
For artificial atom supporting mode $a,b$ (with resonant frequency
$\omega_{\mathrm{a}},\omega_{\mathrm{b}}$), drive field with frequency
$\omega_{\mathrm{EO}}=\omega_{A}-\omega_{0}$ can be used to induce
the coupling between photons of frequency $\omega_{0}$ and the mode
$A$ ($A\in\left\{ \mathrm{a},\mathrm{b}\right\} $). Such an antenna
can be physically realized with a microcavity and electro-optic (EO)
effect \cite{Zhang2019}. A microwave driving is applied on the two
cavities which shifts the two cavity resonant frequencies to opposite
directions, with the photonic molecular Hamiltonian as
\begin{eqnarray}
H & = & \omega_{m,1}m_{1}^{\dagger}m_{1}+\omega_{m,2}m_{2}^{\dagger}m_{2}+J\left(m_{1}^{\dagger}m_{2}+m_{1}m_{2}^{\dagger}\right)+G\cos(\omega_{\mathrm{EO}}t)\left(m_{1}^{\dagger}m_{1}-m_{2}^{\dagger}m_{2}\right),
\end{eqnarray}
where $\omega_{m,i}$ is the resonant frequency of the cavity mode
$m_{i}$, $J$ is the linear coupling strength between the two cavities,
$G$ is the coupling strength due to $\chi^{(2)}$ interaction between
photonic mode and electric field, i.e. the EO drive. Denoting $\omega_{0}=\frac{\omega_{m,1}+\omega_{m,2}}{2}$,
$\delta\omega=\frac{\omega_{m,1}-\omega_{m,2}}{2}$, $\tan\theta=\frac{J}{\delta\omega}$
and making the linear transformation 
\begin{eqnarray}
b & = & \cos\frac{\theta}{2}m_{1}+\sin\frac{\theta}{2}m_{2},\\
d & = & -\sin\frac{\theta}{2}m_{1}+\cos\frac{\theta}{2}m_{2},
\end{eqnarray}
the Hamiltonian transforms to
\begin{eqnarray}
H & = & \omega_{b}b^{\dagger}b+\omega_{d}d^{\dagger}d+G\cos(\omega_{\mathrm{EO}}t)\sin\theta\left(b^{\dagger}d+bd^{\dagger}\right)+G\cos(\omega_{\mathrm{EO}}t)\cos\theta\left(b^{\dagger}b-d^{\dagger}d\right),
\end{eqnarray}
where $\omega_{b}=\omega_{0}+\sqrt{\delta\omega^{2}+J^{2}}$, $\omega_{d}=\omega_{0}-\sqrt{\delta\omega^{2}+J^{2}}$.
For modulation frequency $\omega_{m}$ much larger than the cavity
decay rates $\kappa_{b},\kappa_{d}$ of the supermodes $b,d$, the
sidebands of mode $b$ and $d$ generated by the frequency modulated
term $\Omega\cos(\omega_{\mathrm{EO}}t)\cos\theta\left(b^{\dagger}b-dd^{\dagger}\right)$
are no longer on resonance, which can be neglected safely. Therefore,
the Hamiltonian reduces to
\begin{eqnarray}
H & = & \omega_{b}b^{\dagger}b+\omega_{d}d^{\dagger}d+G\cos(\omega_{\mathrm{EO}}t)\sin\theta\left(b^{\dagger}d+bd^{\dagger}\right).
\end{eqnarray}
 If only mode $m_{2}$ couples with the waveguide with rate $\sqrt{2\kappa_{m,2}}$,
mode $b$ and $d$ experiences effective coupling strength of $\sqrt{2\kappa_{m,2}}\cos\frac{\theta}{2}$
and $\sqrt{2\kappa_{m,2}}\sin\frac{\theta}{2}$, respectively. For
the case of large detuning $\delta\omega\gg\left\{ J,\:\kappa\right\} $,
$b\approx m_{1}$, $d\approx m_{2}$, mode $b$ is approximately decoupled
from the waveguide and its lifetime is only limited by the intrinsic
loss. In the rotating frame of $\omega_{b}b^{\dagger}b+\omega_{d}d^{\dagger}d$
and apply the rotating wave approximation, the Hamiltonian is
\begin{eqnarray}
H & = & \frac{1}{2}\Omega\sin\theta\left(b^{\dagger}de^{i\delta}+bd^{\dagger}e^{-i\delta}\right),
\end{eqnarray}
where $\delta=\omega_{d}-(\omega_{b}-\omega_{d})$. Transforming the
Hamiltonian further to the frame of $\frac{\delta}{2}c^{\dagger}c-\frac{\delta}{2}d^{\dagger}d$,
we obtain the time-independent form
\begin{eqnarray}
H & = & -\frac{\delta}{2}b^{\dagger}b+\frac{\delta}{2}d^{\dagger}d+\frac{1}{2}G\sin\theta\left(b^{\dagger}d+bd^{\dagger}\right).
\end{eqnarray}
In this way, we obtain the coupling between mode $b$ ($m_{1}$) in
the artificial atom cavity and the mode $d$ ($m_{2}$) in the antenna
cavity. The frequency difference between mode $d$ and $b$ is compensated
by the frequency of the EO drive. For input photon with pulse shape
$A_{in}\left(t\right)$, it is rapidly coupled into the antenna cavity
due to the large waveguide-cavity coupling rate and then coupled to
the artificial atom by a time-dependent $G\left(t\right)$.For two
incident photons with the same frequency, we use different $\omega_{\mathrm{EO}}$
to couple them with energy levels of different frequencies without
crosstalk.

\subsection{Quantum storage of control photon pulse}

We use the photonic molecule introduced in the last section for the
storage and readout of the control photon from the waveguide to the
artificial atom. The total Hamiltonian of the system is 
\begin{eqnarray}
H & = & \sum_{j}\omega_{j}j^{\dagger}j+g_{\mathrm{nd}}abc^{\dagger}+\Omega\left(t\right)bd^{\dagger}+h.c.
\end{eqnarray}
where $j\in\{a,b\}$ represents the mode in the atomic cavity in the
fundamental frequency band, $c$ is the mode in the atomic cavity
in the SH frequency band, $d$ is the mode in the antenna cavity.
The dynamics of the cavity field operator $A$ follow 
\begin{eqnarray}
\frac{d}{dt}A & = & -i\left[A,H\right]-\kappa_{A}A+\sqrt{2\kappa_{A,1}}A_{in}\left(t\right).
\end{eqnarray}
The information of the pulse shape is contained in the time dependent
term $A_{in}\left(t\right)$.

Here, we investigate the quantum storage of single photon with the
photonic molecule. Without the control photon, the strong $\chi^{(2)}$
interaction is forbidden, the Hamiltonian for photon storage is
\begin{eqnarray}
H & = & \delta_{b}b^{\dagger}b+\delta_{d}d^{\dagger}d+\Omega b^{\dagger}d+\Omega^{*}bd^{\dagger}.\label{eq:storage}
\end{eqnarray}
For single-photon excitation, the state can be described by a pure
state with a general form of the superposition of Fock states. The
single-photon state is expressed as
\begin{eqnarray}
|\psi\rangle_{db} & = & c_{10}|10\rangle+c_{01}|01\rangle,
\end{eqnarray}
where $|mn\rangle=|m\rangle_{d}\otimes|n\rangle_{b}$. The dynamics
follows
\begin{eqnarray}
\frac{d}{dt}c_{mn} & = & \frac{d}{dt}\langle mn|\psi\rangle\nonumber \\
 & = & \langle mn|\frac{d}{dt}\psi\rangle\nonumber \\
 & = & -i\langle mn|H_{\mathrm{eff}}|\psi\rangle\nonumber \\
 & = & \left(-i\delta_{b}-\kappa_{b}\right)\sum_{m'n'}c_{m'n'}\langle mn|b^{\dagger}b|m'n'\rangle+\left(-i\delta_{d}-\kappa_{d}\right)\sum_{m'n'}c_{m'n'}\langle mn|d^{\dagger}d|m'n'\rangle\nonumber \\
 &  & -i\left(\Omega\sum_{m'n'}c_{m'n'}\langle mn|b^{\dagger}d|m'n'\rangle+\Omega^{*}\sum_{m'n'}c_{m'n'}\langle mn|bd^{\dagger}|m'n'\rangle\right)\nonumber \\
 & = & \alpha_{b}nc_{mn}+\alpha_{d}mc_{mn}-i\Omega\sqrt{(m+1)n}c_{m+1,n-1}-i\Omega^{*}\sqrt{m(n+1)}c_{m-1,n+1},
\end{eqnarray}
where $\alpha_{i}=-i\delta_{i}-\kappa_{i}$, $H_{\mathrm{eff}}$ is
the effective Hamiltonian that includes the mode dissipation, which
is obtained by replacing $\delta_{i}$ by $\alpha_{i}$ in Eq.$\:$\ref{eq:storage}.
Under single-photon pulse excitation of a perfectly phase-matched
cavity $\alpha_{i}=\kappa_{i}$, the dynamics of the intracavity and
output fields follow 
\begin{eqnarray}
\frac{d}{dt}c_{10} & = & -\kappa_{d}c_{10}-i\Omega^{*}\left(t\right)c_{01}+\sqrt{2\kappa_{d,1}}\phi_{in}(t),\label{eq:cme1-1}\\
\frac{d}{dt}c_{01} & = & -\kappa_{b}c_{01}-i\Omega\left(t\right)c_{10},\label{eq:cme1-2}\\
\phi_{out} & = & \phi_{in}-\sqrt{2\kappa_{c,1}}c_{10},\label{eq:cme1-3}
\end{eqnarray}
where $\phi_{in}(t)$ ($\phi_{out}(t)$) represents the pulse shape
of the input (output) single photon.

\subsubsection{Figure of merit}

We assume that both cavities are in the ground state, so $c_{10}=c_{01}=0$
at $t=0$. We also assume that the input field excitation $\phi_{in}(t)$
only has nonzero values in time interval $[0,T]$. Our goal is to
store the input field to $|01\rangle$ and retrieve the photon out
of the cavity after. During the storage process, $c_{10}(0)=c_{01}(0)=0$
and the input field fulfills $\int_{0}^{T}dt|\phi_{in}(t)|^{2}=1$.
The storage efficiency is given by
\begin{eqnarray}
\eta_{s} & = & \frac{\mathrm{Population}\:\mathrm{of\:|01\rangle}}{\mathrm{Number\:\mathrm{of\:\mathrm{incoming}\:\mathrm{photons}}}}\nonumber \\
 & = & |c_{01}(T)|^{2}.
\end{eqnarray}
For the retrieval process, we set $c_{10}(t_{r})=0$, $c_{01}(t_{r})=1$
with $t_{r}$ being the start time of the retrieve. The retrieval
efficiency is given by
\begin{eqnarray}
\eta_{r} & = & \frac{\mathrm{Number}\:\mathrm{of\:\mathrm{retrieved}\:\mathrm{photons}}}{\mathrm{Population}\:\mathrm{of\:|01\rangle}}\nonumber \\
 & = & \int_{t_{r}}^{\infty}|\phi_{out}(t)|^{2}dt.\label{eq:retrive}
\end{eqnarray}


\subsubsection{Photon Retrieval}

First, we investigate the retrieval process with $\phi_{in}(t)=0$
for $t>t_{r}$. In our case, we have assumed that the second cavity
has a quality factor much larger than the first one which couples
with a bus waveguide, thus it is reasonable to make a approximation
$\kappa_{b}\approx0$. From Eqs.~(\ref{eq:cme1-1})-(\ref{eq:cme1-2}),
\begin{eqnarray}
\frac{d}{dt}\left(|c_{10}|^{2}+|c_{01}|^{2}\right) & = & -2\kappa_{d}|c_{10}|^{2}.
\end{eqnarray}
Using Eq.~(\ref{eq:cme1-3}) and Eq.~(\ref{eq:retrive}), the retrieval
efficiency is derived as
\begin{eqnarray}
\eta_{s} & = & 2\kappa_{d,1}\int_{t_{r}}^{\infty}dt|c_{10}|^{2}\nonumber \\
 & = & -\frac{\kappa_{d,1}}{\kappa_{d}}\int_{t_{r}}^{\infty}dt\frac{d}{dt}\left(|c_{10}|^{2}+|c_{01}|^{2}\right)\nonumber \\
 & = & \frac{\kappa_{d,1}}{\kappa_{d,0}+\kappa_{d,1}}\left(|c_{10}(t_{r})|^{2}+|c_{01}(t_{r})|^{2}-|c_{10}(\infty)|^{2}-|c_{01}(\infty)|^{2}\right).
\end{eqnarray}
If there is no excitation in the cavity after the retrieval process,
$|c_{10}(\infty)|^{2}=|c_{01}(\infty)|^{2}=0$, the efficiency reduces
to
\begin{equation}
\eta_{s}=\frac{\kappa_{d,1}}{\kappa_{d,0}+\kappa_{d,1}}.
\end{equation}
Therefore, the retrieval efficiency is independent of the pulse shape
of the driving fields, provided that the driving field pumps all excitation
out of the cavity completely. To achieve high retrieval efficiency,
over-coupled condition of the waveguide-cavity system is required.

In the following, we derive the relationship between the pulse shapes
of the driving field and the output field. For adiabatic retrieval
process, the input and driving fields are smooth. It is reasonable
to adiabatically eliminate $c_{10}$, the population of state $|10\rangle$,
in the coupled equations, if $\kappa_{d}\gg g/T$. Then Eq.~(\ref{eq:cme1-1})
transforms to
\begin{eqnarray}
-\kappa_{d}c_{10}-i\Omega c_{01} & = & 0.
\end{eqnarray}
Combining it with Eq.~(\ref{eq:cme1-2}), we get 
\begin{eqnarray}
\frac{d}{dt}c_{01} & = & -\frac{|\Omega(t)|^{2}}{\kappa_{d}}c_{01}.
\end{eqnarray}
Solving this differential equation and using the boundary condition
$c_{01}(0)=1$ (the start time of the retrieval process is set to
be $0$ for simplicity),
\begin{eqnarray}
c_{01}(t) & = & e^{-\frac{1}{\kappa_{d}}\int_{0}^{t}|g(t')|^{2}dt'}.
\end{eqnarray}
Denote $h(t,t')=\int_{t}^{t'}|\Omega(t^{''})|^{2}dt^{''}$ and submit
is into Eq.~(\ref{eq:cme1-3}), the output field is derived as
\begin{eqnarray}
\phi_{out}(t) & = & i\frac{\sqrt{2\kappa_{d.1}}}{\kappa_{d}}\Omega(t)e^{-\frac{1}{\kappa_{d}}h(0,t)}.\label{eq:out}
\end{eqnarray}
Computing the retrieval efficiency using Eq.~(\ref{eq:out}),
\begin{eqnarray}
\eta_{r} & = & \int_{0}^{\infty}dt|\phi_{out}|^{2}\nonumber \\
 & = & \frac{2\kappa_{d,1}}{\kappa_{d}^{2}}\int_{0}^{\infty}dt|\Omega(t)|^{2}e^{-\frac{2}{\kappa_{d}}\int_{0}^{t}dt'|g(t')|^{2}}\nonumber \\
 & = & \frac{\kappa_{d,1}}{\kappa_{d}}\left(e^{-\frac{2}{\kappa_{d}}h(0,0)}-e^{-\frac{2}{\kappa_{d}}h(0,\infty)}\right)\nonumber \\
 & = & \frac{\kappa_{d,1}}{\kappa_{d,0}+\kappa_{d,1}}\left(1-e^{-\frac{2}{\kappa_{d}}h(0,\infty)}\right).\label{eq:out-prac}
\end{eqnarray}
For strong or long driving pulse, $h(0,\infty)\gg1,$\textbf{ }the
retrieval efficiency reduces to the previously derived form $\frac{\kappa_{d,1}}{\kappa_{d,0}+\kappa_{d,1}}$.

To be more practical, we consider the case that not all excitation
are driven out of the cavity, which means either $c_{01}(\infty)$
or $c_{10}(\infty)$ does not equal to zero. The efficiency can be
obtained by giving $h(0,\infty)$ an finite value. Here, we shape
the driving field $\Omega(t)$ to retrieve the photon to the desired
mode $A_{in}(t)$ (the input shpae). Equaling $\sqrt{\frac{\kappa_{d,1}}{\kappa_{d,0}+\kappa_{d,1}}}A_{in}(t)$
and eq.\ref{eq:out} and their norm square integrals,
\begin{eqnarray}
\sqrt{\frac{\kappa_{d,1}}{\kappa_{d,0}+\kappa_{d,1}}}A_{in}(t) & = & i\frac{\sqrt{2\kappa_{d.1}}}{\kappa_{d}}\Omega(t)e^{-\frac{1}{\kappa_{d}}h(0,t)},
\end{eqnarray}
so
\begin{eqnarray}
\int_{0}^{t}dt'|A_{in}(t)|^{2} & = & 1-e^{-\frac{2}{\kappa_{d}}h(0,t)}.
\end{eqnarray}
We derive the optimal retrieval pulse shape
\begin{eqnarray}
\Omega(t) & = & -i\sqrt{\frac{\kappa_{d}}{2}}e^{h(0,t)/\kappa_{c}}A_{in}(t)\nonumber \\
 & = & -i\sqrt{\frac{\kappa_{d}}{2}}\frac{A_{in}(t)}{\sqrt{\int_{t}^{\infty}dt'|A_{in}(t')|^{2}}}.
\end{eqnarray}
For photon storage, it is the time reverse of the retrieve process
\cite{storage}. One can use driving pulse $\Omega^{*}(t_{r}-t)$
to store input pulse of shape $A_{in}\left(t\right)$ with maximal
efficiency.

\subsection{Transmission of the target photon}

The target photon ($a_{in}\left(t\right)$) can also be coupled to
the artificial atom mediated by the antenna under an EO drive. However,
when the two modes $a$ and $b$ in the artificial atom are designed
to fulfill phase-match condition with mode $c$, whether the target
photon can enter into mode $b$ of the artificial atom ($|1_{a}1_{b}\rangle$)
strongly depends on the quantum state of mode $a$. For sufficiently
large nonlinear coupling strength $g_{\mathrm{nd}}$, the degeneracy
of two-excitation states $|1_{a}1_{b}0_{c}\rangle$ and $|0_{a}0_{b}1_{c}\rangle$
will be removed, resulting in two new eigenstates with the form of
superposition of $|1_{a}1_{b}0_{c}\rangle$ and $|0_{a}0_{b}1_{c}\rangle$
and eigenenergiessplited by $2g_{\mathrm{nd}}$. The sum of the two
photon frequencies no longer equals to the resonant frequency of two-excitation
states. If the control photon has already been stored into mode $b$,
the target photon is no longer resonant with the transition from $|0_{a}1_{b}0_{c}\rangle$
to $|1_{a}1_{b}0_{c}\rangle$, thus will be blocked. Accordingly,
the phase of the target photon is controlled by the state of the stored
control photon in the artificial atom.

In this work, we take the artificial atom based on non-degenerate
three-wave mixing for constructing quantum phase gate, thus the total
system Hamiltonian is 
\begin{eqnarray}
H & = & \omega_{\mathrm{d}}d^{\dagger}d+\omega_{\mathrm{a}}a^{\dagger}a+\omega_{\mathrm{a}}b^{\dagger}b+\delta_{c}c^{\dagger}c+\Omega a^{\dagger}d+g_{nd}abc^{\dagger}+h.c.
\end{eqnarray}
 For EO driving strength $\kappa_{A}\ll\Omega\ll g_{nd}$, the waveguide
is effectively over-coupled with the atom. The target photon couples
into and out from the atom for the control photon in state $|0\rangle$,
whereas transmit directly without entering the atom for control photon
in state $|1\rangle$. In this case, it is no longer necessary to
design the shape of the EO drive, but requires carefully choosing
an appropriate strength.

In the two-excitation subspace, the state can be described by a pure
state with a general form of the superposition of Fock states. For
mode $a$ in state $|1\rangle$ and very long lifetime, the single-photon
state is expressed as
\begin{eqnarray}
|\psi\rangle_{dbac} & = & c_{1010}|1010\rangle+c_{0110}|0110\rangle+c_{0001}|0001\rangle,
\end{eqnarray}
where $|mnkl\rangle=|m\rangle_{d}\otimes|n\rangle_{a}\otimes|k\rangle_{b}\otimes|l\rangle_{c}$.
The dynamics follows
\begin{eqnarray}
\frac{d}{dt}c_{mnkl} & = & \frac{d}{dt}\langle mnkl|\psi\rangle\nonumber \\
 & = & \langle mnkl|\frac{d}{dt}\psi\rangle\nonumber \\
 & = & -i\langle mnkl|H_{eff}|\psi\rangle\nonumber \\
 & = & \alpha_{d}\sum_{m'n'k'l'}c_{m'n'k'l'}\langle mnkl|d^{\dagger}d|m'n'k'l'\rangle+\alpha_{b}\sum_{m'n'k'l'}c_{m'n'k'l'}\langle mnkl|b^{\dagger}b|m'n'k'l'\rangle\nonumber \\
 &  & +\alpha_{a}\sum_{m'n'k'l'}c_{m'n'k'l'}\langle mnkl|a^{\dagger}a|m'n'k'l'\rangle+\alpha_{c}\sum_{m'n'k'l'}c_{m'n'k'l'}\langle mnkl|c^{\dagger}c|m'n'k'l'\rangle\nonumber \\
 &  & -i\Omega\left(\sum_{m'n'k'l'}c_{m'n'k'l'}\langle mnkl|a^{\dagger}d|m'n'k'l'\rangle+\sum_{m'n'k'l'}c_{m'n'k'l'}\langle mnkl|ad^{\dagger}|m'n'k'l'\rangle\right)\nonumber \\
 &  & -ig_{nd}\left(\sum_{m'n'k'l'}c_{m'n'k'l'}\langle mnkl|a^{\dagger}b^{\dagger}c|m'n'k'l'\rangle+\sum_{m'n'k'l'}c_{m'n'k'l'}\langle mnkl|abc^{\dagger}|m'n'k'l'\rangle\right)\nonumber \\
 & = & \alpha_{d}mc_{mnkl}+\alpha_{a}nc_{mnkl}+\alpha_{b}kc_{mnkl}+\alpha_{c}lc_{mnkl}\nonumber \\
 &  & -i\Omega\sqrt{m(n+1)}c_{(m-1)(n+1)kl}-i\Omega\sqrt{(m+1)n}c_{(m+1)(n-1)kl}\nonumber \\
 &  & -ig_{nd}\sqrt{nk(l+1)}c_{m(n-1)(k-1)(l+1)}-ig_{nd}\sqrt{(n+1)(k+1)l}c_{m(n+1)(k+1)(l-1)}.
\end{eqnarray}
For our two-photon case
\begin{eqnarray}
\frac{d}{dt}c_{1010} & = & \alpha_{d}c_{1010}+\alpha_{b}c_{1010}-i\Omega c_{0110}-i\phi_{in}(t),\label{eq:cme2-1}\\
\frac{d}{dt}c_{0110} & = & \alpha_{b}c_{0110}+\alpha_{a}c_{0110}-i\Omega c_{1010}-ig_{nd}c_{0001},\label{eq:cme2-2}\\
\frac{d}{dt}c_{0001} & = & \alpha_{c}c_{0001}-ig_{nd}c_{0110}.\label{eq:cme2-3}
\end{eqnarray}
The dynamics of the system depends on the initial state of the system.
Initially, we set $c_{1010}=c_{0110}=c_{0001}=0$ to solve the dynamics
of the system.
\end{document}
