ket vectors represent the Fock states of the KPO and the output mode, respectively, and $N_l$ ($l=0, 1, \ldots , 6$) is the number at which the photon number in the KPO is truncated when the photon number in the output mode is $l$. In the present simulations, we set $N_0=N_1=N_2=6$, $N_3=5$, $N_4=4$, $N_5=3$, and $N_6=2$. The second ket vector is defined with the creation operators, e.g., as follows: $$\begin{aligned}
|j_1,j_2,j_3\rangle
=\mathcal{N}(j_1,j_2,j_3) \tilde{b}^{\dagger}_{j_1} \tilde{b}^{\dagger}_{j_2} \tilde{b}^{\dagger}_{j_3} |0\rangle,\end{aligned}$$ where the normalization factor $\mathcal{N}$ is defined as $$\begin{aligned}
\mathcal{N}(j_1,j_2,j_3)
=
\left\{
\begin{matrix}
1 & \cdots & j_1>j_2>j_3, \\
1/\sqrt{2!} & \cdots & j_1=j_2>j_3, \\
1/\sqrt{2!} & \cdots & j_1>j_2=j_3, \\
1/\sqrt{3!} & \cdots & j_1=j_2=j_3. \\
\end{matrix}
\right.\end{aligned}$$ Using this representation, we numerically solved the Schrödinger equation with the Hamiltonian in Eq. (8) in the main text. Since the Hamiltonian includes only one of $\{ \tilde{b}_j \}$, the corresponding Schrödinger equation is simple. Moreover, at time $t \in [z_{j-1}, z_j]$, $|\psi \rangle$ includes only the output-mode photons satisfying $j_1 \le j$ in Eq. (\[eq-psi\]). This enables fast implementation of the simulation.

Since the photon number in the KPO is large (small) in the first (second) half of the whole process, we set correspondingly 