system of differential equations. In this work, we instead discretize the position, which results in simpler equations as shown below.

We first introduce the annihilation operator with respect to the position $z$: $$\begin{aligned}
\tilde{b}(z)
&=
\frac{1}{\sqrt{2\pi}}
\int_{-\infty}^{\infty} \!
b(\omega) e^{i \omega z} \, d\omega.\end{aligned}$$ Then, we move to the interaction picture with the unitary operator $U(t)=e^{-iH_{\mathrm{out}}t}$ as follows: $$\begin{aligned}
\tilde{b}_I(z,t)
&=
U^{\dagger}(t) \tilde{b}(z) U(t)
=
\tilde{b}(z-t),
\nonumber
\\
H_I(t)
&= 
U^{\dagger}(t) \left[ H_{\mathrm{KPO}}(t)
+ H_c \right] U(t)
\nonumber
\\
&=
H_{\mathrm{KPO}}(t)
+
i \sqrt{\kappa_{\mathrm{ex}}}
\left[
\tilde{b}^{\dagger}(-t) a
-
a^{\dagger} \tilde{b}(-t)
\right].
\nonumber\end{aligned}$$ A pulse-mode operator for the interval $[0,T]$ is defined as $$\begin{aligned}
b_p = \int_0^{T} \! f_p(z) \tilde{b}_I(z,T) \, dz
=
\int_0^{T} \! f_p(T - z) \tilde{b}(-z) \, dz,
\nonumber\end{aligned}$$ where $f_p(z)$ is the normalized envelope function for the output pulse satisfying $\displaystyle \int_{0}^T \! |f_p(z)|^2 dz =1$. In the present work, we define $f_p(z)$ as ${f_p(z) \propto \sqrt{\langle \tilde{b}^{\dagger}_I(z,T) \tilde{b}_I(z,T) \rangle}}$, where $\langle \tilde{b}^{\dagger}_I(z,T) \tilde{b}_I(z,T) \rangle$ is the spatial distribution of photons in the output mode at the final time $T$. Note that in experiments, we can find such $f_p(z)$ from the measurement of the output power from the KPO.

Next, we divide the interval $[0,T]$ into $J$ 