$S(\tilde{z},\tau)$ and the field $\mathcal{E}(\tilde{z},\tau)$ are unchanged, the definitions of the storage, retrieval, and total efficiencies remain the same. Furthermore, the boundary conditions for storage are the same as before except $P_2(\tilde{z},0)=0$ is added.

We can also find the four-level version of $J$ $$\begin{aligned}
J &=\int_0^1 d\tilde{z}\,S(\tilde{z},T)S^*(\tilde{z},T)+\int_0^T d\tau \int_0^1 d\tilde{z}\,\Big\{ \bar{\mathcal{E}}^* \Big[ -\partial_{\tilde{z}} \mathcal{E} \nonumber\\
&+ i \sqrt{d \gamma}\left(\tilde{\mu}_{1g} P_1 +\tilde{\mu}_{2g} P_2\right)\Big] +c.c.\Big\} \nonumber \\
&+\int_0^T d\tau \int_0^1 d\tilde{z}\, \Big\{ \bar{P}_1^*\Big[(-\partial_{\tau}+i \delta_g-\gamma) P_1+i \tilde{\mu}_{1s} \Omega(\tau) S  \nonumber \\
&+ i \tilde{\mu}_{1g} \sqrt{d \gamma} \mathcal{E} \Big] +c.c. \Big\} \nonumber \\
&+\int_0^T d\tau \int_0^1 d\tilde{z}\, \Big\{ \bar{P}_2^*\Big[ \left(-\partial_{\tau}+i\delta_g-i\Delta_e-\gamma  \right) P_2  \nonumber \\
&+i \tilde{\mu}_{2s} \Omega(\tau) S+ i \tilde{\mu}_{2g} \sqrt{d \gamma} \mathcal{E} \Big] +c.c. \Big\} \nonumber \\
&+\int_0^T d\tau \int_0^1 d\tilde{z}\,\Big\{ \bar{S}^*\Big[-\partial_{\tau}S+ i \left(\delta_g - \delta_s \right) S \nonumber \\
&+ i\Omega^*(\tau)\left( \tilde{\mu}_{1s}P_1+\tilde{\mu}_{2s}P_2  \right) \Big] +c.c. \Big\}, \label{eqn:4LJ}\end{aligned}$$ in order to obtain the optimal control and efficiencies. By requiring that $J$ is stationary with respect to variations in $\mathcal{E}$, $P_1$, $P_2$, and $S$, we find the equations of motion $$\begin{aligned}
\partial_{\tilde{z}} \bar{\mathcal{E}} &= i\sqrt{d \gamma}\left( \tilde{\mu}_{1g}\bar{P}_1+\tilde{\mu}_{2g}\bar{P}_2 \right)\label{eqn:4LEoMbE}\\
\partial_{\tau}\bar{P}_1 &= \left(i \delta_g+\gamma \right) \bar{P}_1+i \tilde{\mu}_{1s} \Omega(\tau) \bar{S} +i \tilde{\mu}_{1g} \sqrt{d \gamma} \bar{\mathcal{E}} \label{eqn:4LEoMbP1}\\
\partial_{\tau}\bar{P}_2 &= \left(i \delta_g-i\Delta_e+\gamma \right) \bar{P}_2+i \tilde{\mu}_{2s} \Omega(\tau) \bar{S} +i \tilde{\mu}_{2g} 