and Lagrange multipliers $\bar{\mathcal{E}}$, $\bar{P}_1$, and $\bar{S}$ have been introduced to include the equations of motion. The maximum storage efficiency is found when $J$ is stationary with respect to variations in all dynamical variables and $\Omega(\tau)$. Requiring stationarity with respect to variations in $\mathcal{E}$, $P_1$, and $S$ results in equations of motion and boundary conditions for the Lagrange multipliers $\bar{\mathcal{E}}$, $\bar{P}_1$, and $\bar{S}$. The equations of motion are $$\begin{aligned}
\partial_{\tilde{z}} \bar{\mathcal{E}} &= i \tilde{\mu}_{1g}\sqrt{d \gamma}\bar{P}_1 \label{eqn:3LEoMbE}\\
\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:3LEoMbP1}\\
\partial_{\tau}\bar{S} &= i \left(\delta_g - \delta_s \right) \bar{S} + i \tilde{\mu}_{1s}\Omega^*(\tau)\bar{P}_1 \label{eqn:3LEoMbS}\end{aligned}$$ with boundary conditions $\bar{S}(\tilde{z},T) = S(\tilde{z},T)$ and $\bar{\mathcal{E}}(1,\tau)=\bar{P}_1(\tilde{z},T)=0$. As pointed out in Ref. [@ref:Gorshkov_PRA4], these are exactly the equations of motion and boundary conditions for backwards retrieval. Thus, solving these equations will yield the retrieval efficiency, $\eta_r$, and the total efficiency $\eta_{tot} = \eta_s \eta_r$ for storage followed by backwards retrieval. Explicitly, these quantities can be determined using $$\begin{aligned}
\eta_{r} &= \frac{\int_0^T d\tau \, |\mathcal{E}_{\textrm{out}}(\tau)|^2}{\int_0^1 d\tilde{z}\left|S(\tilde{z},T)\right|^2} \\
\eta_{tot} &= \int_0^T d\tau \, |\mathcal{E}_{\textrm{out}}(\tau)|^2\end{aligned}$$ where $\mathcal{E}_{\textrm{out}}(\tau) = \bar{\mathcal{E}}(0,\tau)$ is the output quantum field.

The variation of $J$ with respect to variations in $\Omega(\tau)$ (the gradient) can be identified 