\section{The QCQP is equivalent to the original control problem}
In the main text, we derive the conservation laws of Eq.~(3) as necessary conditions that must be satisfied by any solution of the control problem. In this section, we show that the conservation laws are also \emph{sufficient} conditions: any solution of the ultimate quadratically constrained quadratic program (QCQP) must also be a solution of the original control problem. Hence the optimal solution of the QCQP must also be the optimal solution of the original problem. The two problems are equivalent in this sense. (A related observation was made in \citeasnoun{Angeris2021}.)

The conservation-law constraints of Eq.~(3) in the main text are:
\begin{align}
                    \int_{t_0}^T \int_{t_0}^T \Phi^{\dagger}(t) D_i(t) \left[ \frac{H_c^{-1}(t)}{\varepsilon_{\rm max}} \delta(t-t') + \frac{i}{\hbar} G_0^{+}(t,t') \right] \Phi(t')\D t \,\D t' = \int_{t_0}^T \Phi^{\dagger}(t) D_i(t) U_0(t,t_0) \d{t},
                    \label{eq:Dconstr}
\end{align}
The QCQP arises from optimizing a linear or quadratic objective $f(\Phi)$ subject to \eqref{Dconstr} being satisfied for all possible $D_i$. Let us consider a solution $\Phi(t)$ of the QCQP that satisfies \eqref{Dconstr} for all $D_i(t)$. In particular, we can assume that it satisfied \eqref{Dconstr} for all $D_i(t)$ of the form
\begin{align}
    D^{i}_{jk}(t) = \delta_{j\ell}\delta_{km} \delta(t-t'),
\end{align}
where we have momentarily moved the $i$ to a superscript to explicitly consider the $jk$ element of $D^{(i)}$ in the Hilbert space, and $t'$ can be any time in the interval of interest. Given this choice of $D^{(i)}$, we can rewrite \eqref{Dconstr} as
\begin{align}
    \Phi_{\ell i}^*(t') \left[\frac{H_c^{-1}(t')}{\varepsilon_{\rm max}}\Phi(t') + \frac{i}{\hbar}\int_{t_0}^T G_0^{+}(t',t'') \Phi(t'')\,{\rm d}t'' - U_0(t',t_0)\right]_{jm} = 0,
    \label{eq:prodform}
\end{align}
which must hold for all $i,j,\ell,m$ and $t'$. Since \eqref{prodform} must hold at all times, we can see that at any time $t'$ there are two possibilities for $\Phi(t')$: either $\Phi(t') = 0$, or $\Phi(t')$ satisfies the expression in square brackets. This dichotomy will dictate how to find a corresponding $U(t,t_0)$ that satisfied the original (Schrodinger) constraint. At times where $\Phi(t')=0$, we can take the control field $\varepsilon(t')=0$. We can also note that the integrand will be zero at such times, and that the domain of the integral can then be restricted to all times for which $\Phi$ is nonzero. \Eqref{prodform} can equivalently be written:
\begin{align}
    \Phi_{\ell i}^*(t') \left[\frac{H_c^{-1}(t')}{\varepsilon_{\rm max}}\Phi(t') + \frac{i}{\hbar}\int_{\left\{t''\in[t_0,T]|\Phi(t'')\neq0\right\}} G_0^{+}(t',t'') \Phi(t'')\,{\rm d}t'' - U_0(t',t_0)\right]_{jm} = 0,
    \label{eq:prodform2}
\end{align}
Now we consider all $t'$ for which $\Phi(t') \neq 0$. The term in square brackets must be zero at all such times. For such times, we can set $\varepsilon(t') = \varepsilon_{\rm max}$, and $U(t',t_0) = (H_c^{-1}(t')/\varepsilon_{\rm max}) \Phi(t')$. These $U(t',t_0)$ must satisfy the term in square brackets, which can now be written:
\begin{align}
    U(t',t_0) + \frac{i}{\hbar}\int_{\left\{t''\in[t_0,T]|\varepsilon(t'')\neq0\right\}} G_0^{+}(t',t'') \varepsilon_{\rm max} H_c(t'') U(t'',t_0)\,{\rm d}t'' - U_0(t',t_0) = 0.
    \label{eq:prodform3}
\end{align}
Finally, since $\varepsilon(t)$ is either $\varepsilon_{\rm max}$ or 0 at all times, we can rewrite the integral domain to take place over all possible times, with $\varepsilon_{\rm max}$ replaced by $\varepsilon(t)$, implying that $U(t,t_0)$ must satisfy
\begin{align}
    U(t',t_0) = U_0(t',t_0) - \frac{i}{\hbar}\int_{t_0}^T G_0^{+}(t',t'') \varepsilon(t'') H_c(t'') U(t'',t_0)\,{\rm d}t''.
    \label{eq:inteq}
\end{align}
\Eqref{inteq} is exactly the integral form of Eq.~(1) in the main text, which is equivalent to the Dyson equation and the differential Schrodinger equation! Hence we have shown that any solution $\Phi(t)$ that satisfies \eqref{Dconstr} for all $D_i$ implies a solution $U(t,t_0)$ of the original dynamical constraints. Therefore the optimal value of the QCQP will also be the optimal value of the original control problem.

\section{Three-level Hamiltonians}
The asymmetric double-well case is taken from Sec.~2.8 of \cite{Werschnik2007},
\begin{eqnarray}
	H &=& \bb{\mat{\omega_0 & 0&0\\ 0 & \omega_1 & 0\\ 0& 0& \omega_2 }} -\epsilon(t) \bb{\mat{\mu_{00} & \mu_{01}&\mu_{02}\\ \mu_{10} & \mu_{11} & \mu_{12}\\ \mu_{20}& \mu_{21}&\mu_{22}}}\nonumber\\
	 &=& \bb{\mat{0 & 0&0\\ 0 & 0.1568 & 0\\ 0& 0& 0.7022 }} -\epsilon(t) \bb{\mat{-2.5676 &0.3921& 0.6382\\ 0.3921 &2.3242& -0.7037\\  0.6382 & -0.7037 &-0.5988}}\,.
\end{eqnarray}
We consider a maximum control amplitude of $\magn{\epsilon} \le 0.15$ in this example.

For the transmon example, the exact Hamiltonian is
\begin{eqnarray}
	H &=& \bb{\mat{\omega_0 & 0&0\\ 0 & \omega_1 & 0\\ 0& 0& \omega_2 }} - \epsilon(t)  \bb{\mat{0 & \mu_{01}& 0 \\ \mu_{10} & 0 & \mu_{12}\\ 0& \mu_{21}& 0}}\nonumber\\
	&=&\bb{\mat{0 & 0&0\\ 0 & 1.9 & 0\\ 0& 0& 3.7 }} + \epsilon(t)  \bb{\mat{0 & 1 & 0 \\ 1 & 0 & \sqrt{2}\\ 0& \sqrt{2} & 0}}\,,
\end{eqnarray}
with a maximum control amplitude $\magn{\epsilon} \le 0.3$ in this example.

\section{Finding a binary pulse equivalent to a continuous one by local averaging}
Given a Hamiltonian $H(t) = H_0 (t) + \varepsilon(t)H_c(t)$ with some control $\varepsilon(t)$, designed controls might be of ``bang--bang'' type, taking only two discrete values, or they might be smooth, continuous controls, possibly bounded in magnitude. From a theoretical bound perspective, the latter case is subsumed by the former: any smooth, continuous control can be approximated to arbitrarily high accuracy with a particular bang--bang control. Intuitively, one can see that this might be true: if you oscillate a bang--bang control at high enough frequency, oscillating over smaller time scales than any transitions (real or virtual) induced by the Hamiltonian, then the wave function will not respond to the particular high-frequency details of the control but rather to a homogenized average (which will lie between the two extremes). Mathematically, there is a rich literature on the field of ``homogenization'' theory~\cite{Milton2002}, for applications from material science to optics. Here, we prove that any smooth, continuous control can be approximated by a bang--bang control, which then allows us to use bang--bang controls in the formulation of our bounds. (Note the key point that such bang--bang controls would not ever need to be implemented; they simply inform us of the generality of the bounds.)

Suppose that we are interested in the time evolution of the system from $t_0$ to $t_f$, with continuous control $\varepsilon_c(t)$ and propagator $U_c(t_f,t_0)$. We can divide $[t_0,t_f]$ into $N$ intervals of equal length $\tau = (t_f - t_0)/N$: $[t_0,t_1], [t_1,t_2], \ldots, [t_{N-1},t_N]$ ($t_f = t_N$), such that each interval is shorter than any timescale of the Hamiltonian. We can then write $U_c (t_f,t_0) = U_c (t_N,t_{N-1})\cdots U_c (t_2,t_1)U_c (t_1,t_0)$. Since each interval is smaller than any transition time scale, we can apply first-order perturbation theory to each $U(t_{i+1},t_i)$ to analyze the accuracy of a bang--bang-control approximation. First, we analyze the accuracy of using a first-order approximation to the time-evolution operator over each interval:
\begin{eqnarray}
	U_c(t_{i+1},t_i)& =& U_0 (t_{i+1},t_i) - \frac{\I}{\hbar}\int_{t_i}^{t_{i+1}}\D t\, U_0 (t_{i+1},t)H_I U_0 (t,t_i) \epsilon_c (t) + e'_i\\
	&&( \mbox{Taylor expand} \quad U_0 (t_{i+1},t)H_I U_0 (t,t_i) \quad \mathrm{around} \quad t = t_i )\nonumber \\
	&=& \equals{U_0 (t_{i+1},t_i) - \frac{\I}{\hbar}U_0 (t_{i+1},t_i)H_I \int_{t_i}^{t_{i+1}}\D t\, \epsilon_c (t)}{\tilde{U}_c(t_{i+1},t_i)} + e_{i}\label{eq:perturb}
\end{eqnarray}
where $U_0$ is the propagator of $H_0$ and $e'_{i}$ is the error term of the perturbation series of $\BigO{\tau^2}$. Together with the error due to Taylor expansion of $O(\tau^2)$, the overall local error $e_i$ is of $O(\tau^2)$. If we stitch these local approximations back to the global propagator, i.e. $U_c (t_f,t_0) = U_c (t_N,t_{N-1})\cdots U_c (t_2,t_1)U_c (t_1,t_0)$, we see that the global error is $\sum_{i=0}^{N-1} U_0 (t_N,t_{i+1}) e_i U_0 (t_{i},t_0)$, which is at most $N$ ($\sim 1/\tau$) copies of $O(\tau^2)$, which is of $O(\tau)$. Hence the global error goes to zero as $\tau \rightarrow 0$, and we can approximate $U_c$ with the local apprxomations $\tilde{U}_c$:
\begin{equation}
	U_c (t_f,t_0) \sim \tilde{U}_c (t_N,t_{N-1})\cdots \tilde{U}_c (t_2,t_1) \tilde{U}_c (t_1,t_0)\quad \mathrm{as} \quad \tau \rightarrow 0.
\end{equation}

The above perturbation analysis shows that if we can find a bang--bang control $\epsilon_b(t)$ with propagator $U_b(t_f,t_0)$ whose local approximation $\tilde{U}_b(t_{i+1},t_i)$ agrees with  $\tilde{U}_c(t_{i+1},t_i)$ up to the first order in $\tau$, then we will have
\begin{equation}
	U_c (t_f,t_0) \sim \tilde{U}_c (t_N,t_{N-1})\cdots \tilde{U}_c (t_1,t_0) \sim \tilde{U}_b (t_N,t_{N-1})\cdots \tilde{U}_b (t_1,t_0) \sim U_b (t_f,t_0)
\end{equation}
as $\tau \rightarrow 0$. This will show the global equivalence between $U_c (t_f,t_0)$ and $U_b (t_f,t_0)$. Such an $\epsilon_b (t)$ is easy to find: by the form of the perturbation in \eqref{perturb}, any choice of $\epsilon_b$ is valid as long as $\int_{t_i}^{t_{i+1}}\D t \, \epsilon_b (t) = \int_{t_i}^{t_{i+1}}\D t \, \epsilon_c (t)$ for all $i$, which can be achieve simply by requiring the average of the bang--bang control equaling the average of the continuous control. (This choice of $\epsilon_b(t) $ over $[t_i,t_{i+1}]$ has $\epsilon_b = \varepsilon_{\rm min}$ for a duration $t'$ and then switch to $\epsilon_b = \varepsilon_{\rm max}$ for the rest of the time to $t_{i+1}$, for $t' = \frac{\varepsilon_{\rm max} \tau - M_i}{\varepsilon_{\rm max} - \varepsilon_{\rm min}}$, where $M_i =  \int_{t_i}^{t_{i+1}}\D t \, \epsilon_c (t)$. Certainly, there are other choices as well.)

We illustrate the above proof numerically through the transmon example in the main text. We apply a continuous control $\varepsilon_c (t) = E \cos(\omega t)$ with $\omega = \omega_1$ and $E = 0.5$ from $t_0=0$ to $T = 12$. We want to find a binary pulse with $\varepsilon_{\rm max} = - \varepsilon_{\rm min} = E$,  such that the dynamics $\psi_{b}(t)$ approaches $\psi_{c}(t)$ under the continuous pulse when $\tau$ goes to zeros. We construct such $\varepsilon_b (t)$ according to the method described in the previous paragraph. In \figref{averaging}, we compare $\psi_b$ with $\psi_c$ for $\tau = 1.2$ for (a), $0.3$ for (c) and $0.15$ for (e) and show their respective pulses in (b), (d) and (f).  We see from In \figref{averaging} the the evolution under each $\varepsilon_b (t)$ (red) gets closer to the evolution under $\varepsilon_c (t)$ (black dash) as $\tau$ decreases. More explicitly, we plot in \figref{converge} the relative difference between evolution of $\varepsilon_c$ and $\varepsilon_b$, measured by $\frac{(\int_{t_0}^T |\psi_{b}(t) - \psi_{c}(t)|^2)^{\frac{1}{2}}}{(\int_{t_0}^T |\psi_{c}(t)|^2)^{\frac{1}{2}}}$, against $\frac{\tau}{T}$. We can see that the convergence demonstrates the $O(\tau)$ behavior proved above.
We emphasize again that the purpose of numerical results here is to show the equivalence between continuous and binary pulses when $\tau$ goes to zero. In the bound computation, we do not need to find such $\varepsilon_b(t)$ and the bound converges automatically and faster than $O(\tau)$ in the proof\,.

To conclude, we have proven that for any continuously valued control $\varepsilon_c(t)$ there is a corresponding binary control $\varepsilon_b(t)$ producing the same time evolution as $\epsilon_c(t)$, when the local averaging time $\tau$ approaches zero. Hence, a bound derived for bang--bang controls will include all possible bounded continuous controls as well.

\begin{figure}[thh]
    \centering
    \includegraphics[width=0.9\textwidth]{averaging}
    \caption{(a) Comparisons between probability time evolution under $\varepsilon_c (t)$(red) and $\varepsilon_b (t)$(black dash) for all three levels, from $t_0 = 0$ to $T = 12$. (b) $\varepsilon_c (t) = E\cos(\omega t)$ is plotted in red and the $\varepsilon_b (t)$ in black, obtained through averaging for $\tau = 1.2$. (c) and (d), (e) and (f) are similar to (a) and (b), but the former pair is for $\tau = 0.3$ and the latter is for $\tau = 0.15$\,. One can see that as $\tau$ gets smaller, the $\varepsilon_b(t)$ oscillates more rapidly and the evolution gets closer to the one produce by $\varepsilon_c (t)$. For $\tau = 0.15$, the effect of $\varepsilon_b(t)$ and $\varepsilon_c(t)$ are almost identical.}
    \label{fig:averaging}
\end{figure}

\begin{figure}[thb]
    \centering
    \includegraphics[width=0.5\textwidth]{conv_plot}
   \caption{The relative difference between evolution of $\varepsilon_c$ and $\varepsilon_b$, measured by $\frac{(\int_{t_0}^T |\psi_{b}(t) - \psi_{c}(t)|^2)^{\frac{1}{2}}}{(\int_{t_0}^T |\psi_{c}(t)|^2)^{\frac{1}{2}}}$ is plotted against $\frac{\tau}{T}$. It shows the $O(\tau)$ convergence for small $\tau$.}
    \label{fig:converge}
\end{figure}

\clearpage
\section{Local probability conservation laws as a subset of D-matrix constraints }
In this section, we show explicitly that Eq. 4 in the main text, which represents local probability conservation, is a subset of constraints  of Eq. 3  in the main text with $D_i (t) = \eta(t_1 - t) \mathcal{I} $ with $t\ge t_0$, i.e. identity operator $\mathcal{I}$ supported from $t_0$ to $t_1$.

Now consider Eq. 2 in the main text\,,
\begin{eqnarray}
                    \int_{t=t_0}^T \int_{t'=t_0}^T \Phi^\dagger (t) D_i(t) \sb{\frac{H^{-1}_c }{\varepsilon_{\rm max}} \delta(t-t') + \frac{i}{\hbar}G_0^+ (t,t')}\Phi (t') \D t \, \D t' = \int_{t_0}^T \Phi^\dagger (t)D_i(t) U_0 (t,t_0) \D t\,,
                    \label{eq:mt3}
\end{eqnarray}
where $G_0^+ (t,t') = U_0 (t,t') \eta(t-t')$.
Now with $D_i (t) = \eta(t_1 - t) \mathcal{I}$, it becomes
\begin{eqnarray}
               \int_{t_0}^{t_1} \Phi^\dagger (t) \frac{H^{-1}_c }{\varepsilon_{\rm max}}   \Phi (t) \D t +  \frac{i}{\hbar} \int_{t=t_0}^{t_1} \int_{t'=t_0}^t \Phi^\dagger (t) U_0(t,t')\Phi (t') \D t \, \D t' = \int_{t_0}^{t_1} \Phi^\dagger (t)U_0 (t,t_0) \D t\,.
                    \label{eq:mt3t1}
\end{eqnarray}
Taking Hermitian conjugation both side of \eqref{mt3t1}, and using $U^\dagger_0 (t,t') = U_0 (t',t)$ and $(H^{-1}_c)^\dagger = H^{-1}_c$\,, produces
\begin{eqnarray}
               \int_{t_0}^{t_1} \Phi^\dagger (t) \frac{H^{-1}_c }{\varepsilon_{\rm max}}   \Phi (t) \D t -  \frac{i}{\hbar} \int_{t=t_0}^{t_1} \int_{t'=t_0}^t \Phi^\dagger (t') U_0(t',t)\Phi (t) \D t \, \D t' = \int_{t_0}^{t_1} U^\dagger_0 (t,t_0) \Phi(t)\D t\,
                    \label{eq:mt3t1cg}
\end{eqnarray}
Note we can exchange the $t,t'$ integration in the the following integral
\begin{eqnarray*}
	\int_{t = t_0}^{t_1} 	\int_{t' = t_0}^{t} \ldots \D t' \D t &=&  \int_{t' = t_0}^{t_1} 	\int_{t = t'}^{t_1} \ldots \D t \D t' \,.
\end{eqnarray*}
As a result, the double integral in \eqref{mt3t1cg} becomes
\begin{eqnarray*}
	 && \int_{t=t_0}^{t_1} \int_{t'=t_0}^t \Phi^\dagger (t') U_0(t',t)\Phi (t) \D t \, \D t'  \quad (\mbox{exchange $t, t'$ integrations})\\
 	 &=& \int_{t' =t_0}^{t_1} \int_{t = t'}^{t_1} \Phi^\dagger (t') U_0(t',t)\Phi (t) \D t \, \D t' \quad (\mbox{relabeling by} \,t\rightarrow t', t' \rightarrow t)\\
 	 &=& \int_{t =t_0}^{t_1} \int_{t' = t}^{t_1} \Phi^\dagger (t) U_0(t,t')\Phi (t') \D t \, \D t'\,.
\end{eqnarray*}
Now we subtract \eqref{mt3t1} by \eqref{mt3t1cg} to obtain the imaginary part of \eqref{mt3t1}, and we get
\begin{equation}
	\frac{1}{\hbar} \int_{t_0}^{t_1} \int_{t_0}^{t_1}\Phi^\dagger (t) U_0(t,t')\Phi (t') \D t \, \D t' = 2\Im \int_{t_0}^{t_1} \Phi^\dagger (t)U_0 (t,t_0) \D t\,,
\end{equation}
which is exactly  Eq. 4 in the main text, since $U^\dagger U = \mathcal{I}$ by probability conservation. $ \int_{t_0}^{t_1} \Phi^\dagger (t) \frac{H^{-1}_c }{\varepsilon_{\rm max}}   \Phi (t) \D t$ disappears since its imaginary part is zero. This makes sense since $H_c$ is a Hermitian matrix and no dissipation is present.

Often $D_i (t) = \delta(t-t')D$ is localized at time $t'$ ($D$ is some constant matrix). To connect such $D_i(t)$ with the above result, we notice that $\int_{t_0} ^{t_1} \delta (t-t') \D t' = \eta(t_1-t)$ for $t \ge t_0$. As a result, we can interpret the above conservation law as the imaginary part of a linear combination  of constraints with $D_i (t) =  \delta(t-t') \mathcal{I}$ for $t_0 \le t' \le t_1$, a subset of all D-matrix constraints. More generally, $D$ is not necessarily $\mathcal{I}$ and the real part is also a valid D-matrix constraint. We thus can interpret the entire set of D-matrix constraints as generalized local probability conservation laws.

\newpage
\section{Bound formulation for problems with multiple controllable terms}
In the main text, we only consider problems with a single control. We show here that our framework extends naturally to problems with multiple controls. Consider a Hamiltonian $H(t) = H_0(t) + \sum_{i=1}^{N}  \varepsilon_i(t) H^i_{c}(t)$, where $H_0(t)$ is the non-controllable part and  $\sum_{i=1}^{N} \varepsilon_i(t) H^i_{c}$ is the controllable part. To avoid overloading the notation, we here formulate the bound problem for a two-term problem with $H(t) = H_0(t) + \varepsilon_1(t) H^1_{c}  +  \varepsilon_2(t) H^2_{c}$, which can be extended to more terms in a straightforward way.

The first step is to write the integral-equation form (analog of Eq.~(1) of the main text) for multiple controls. To do this, we will assume the two controls cannot be turned on at the same times. (Note that two controls that can be turned on at the same time is equivalent to three controls that cannot be turned on at the same time, so there is no loss of generality here.) With two independent controls, the integral equation for $U(t,t_0)$ is:
\begin{equation}
    U(t,t_0) = U_0(t,t_0) -\frac{i}{\hbar}\int_{t_0}^{T} G_0^+(t,t') H^1_c(t) \varepsilon_1(t')U(t',t_0) - \frac{i}{\hbar}\int_{t_0}^{T} G_0^+(t,t') H^2_c(t) \varepsilon_2(t')U(t',t_0)
\end{equation}
For times where $\varepsilon_2(t)$ is zero, we have
\begin{equation}
	U(t,t_0) = U_0(t,t_0) -\frac{i}{\hbar}\int_{\cb{\varepsilon_1(t) \ne  0}} G_0^+(t,t') H^1_c \varepsilon_1(t')U(t',t_0) \D t'\,,
	\label{eq:precon1}
\end{equation}
 and for times where $\varepsilon_1(t)$ is zero,
\begin{equation}
	U(t,t_0) = U_0(t,t_0) -\frac{i}{\hbar}\int_{\cb{\varepsilon_2(t) \ne  0}} G_0^+(t,t') H^2_c \varepsilon_2(t')U(t',t_0) \D t'\,.
	\label{eq:precon2}
\end{equation}
We then left multiply $U^\dagger(t,t_0)\varepsilon_1(t)H^1_c D_i(t)$ to \eqref{precon1} and $U^\dagger(t,t_0)\varepsilon_2(t)H^2_cD_j(t) $ to \eqref{precon2}, and integrate over $t$ from $t_0$ to $T$. Then \eqref{precon1}  and \eqref{precon2} become
\begin{eqnarray}
	&&\int_{t_0}^{T}\D t \, U^\dagger(t,t_0)\varepsilon_1(t)H^1_c D_i(t) U(t,t_0)
	= \int_{t_0}^{T}\D t \, U^\dagger(t,t_0)\varepsilon_1(t)H^1_c D_i(t)U_0(t,t_0) \nonumber\\
	&& -\frac{i}{\hbar}\int_{t_0}^{T}\D t \int_{\cb{\varepsilon_1(t) \ne 0}} \D t' \, U^\dagger(t,t_0)\varepsilon_1(t)H^1_cD_i(t)G_0^+(t,t') H^1_c \varepsilon_1(t')U(t',t_0)\,,
\end{eqnarray}
and
\begin{eqnarray}
	&&\int_{t_0}^{T}\D t \, U^\dagger(t,t_0)\varepsilon_2(t)H^2_c D_j(t)U(t,t_0)
	= \int_{t_0}^{T}\D t \, U^\dagger(t,t_0)\varepsilon_2(t)H^2_c D_j(t) U_0(t,t_0) \nonumber\\
	&& -\frac{i}{\hbar}\int_{t_0}^{T}\D t \int_{\cb{\varepsilon_2(t) \ne 0}} \D t' \, U^\dagger(t,t_0)\varepsilon_2(t)H^2_c D_j(t)G_0^+(t,t') H^2_c \varepsilon_2(t')U(t',t_0)\,,
\end{eqnarray}
We define $\Phi_1(t) = \varepsilon_1(t)H^1_c U(t,t_0)$ and $\Phi_2(t) = \varepsilon_2(t)H^2_c U(t,t_0)$.  The above equations are quadratic in $\Phi_1$ and $\Phi_2$, thus they hold for $t'$ from $t_0$ and $T$. Consequently, we can replace the integral range of $t'$ by $t_0$ and $T$ and obtain constraints
\begin{equation}
	\int_{t_0}^{T}\D t \int_{t_0}^{T}\D t' \,  \Phi^\dagger_1(t)D_i(t)\bb{ \delta(t-t')\frac{(H^1_c)^{-1}}{\varepsilon_{1\mathrm{max}}} + \frac{i}{\hbar}G_0^+(t,t')}  \Phi_1(t') = \int_{t_0}^{T}\D t \,  \Phi^\dagger_1(t) D_i(t)U_0(t,t_0) \,,
	\label{eq:con1}
\end{equation}
and
\begin{equation}
	\int_{t_0}^{T}\D t \int_{t_0}^{T}\D t' \,  \Phi^\dagger_2(t)D_j(t)\bb{ \delta(t-t')\frac{(H^2_c)^{-1}}{\varepsilon_{2\mathrm{max}}} + \frac{i}{\hbar}G_0^+(t,t')}  \Phi_2(t') = \int_{t_0}^{T}\D t \,  \Phi^\dagger_2(t)D_j(t) U_0(t,t_0) \,,
	\label{eq:con2}
\end{equation}
similar to Eq. (3) in the main text, but now with two different components. In addition to these two sets of constraints, we also need to impose that only one of $\varepsilon_1(t)$ and $\varepsilon_2(t)$ can be nonzero at any given time, which is represented in the constraint
\begin{equation}
    \Phi_1(t)\Phi_2^\dagger(t) = 0\, \quad \mbox{for all components},
	\label{eq:exclude}
\end{equation}
for all times $t$. Hence bounds for two controls can be formulated as
                \begin{equation}
                    \begin{aligned}
                        & \underset{\Phi_1,\Phi_2}{\text{max.}} & & f(\Phi_1,\Phi_2) \\
                        & \text{s.t.} & & \textrm{\Eqref{con1} satisfied for all } D_i(t)\,,\nonumber\\
					  &  & & \textrm{\Eqref{con2} satisfied for all } D_j(t)\,,\\
					  &&& \Phi_1(t)\Phi_2^\dagger(t)= 0 \quad \mbox{for all } t\,, \mbox{for all components}\,.
                    \end{aligned}
                    \label{eq:opt_prob}
                \end{equation}
                One can see immediately how to generalize to multiterm pulse designs. For $H(t) = H_0(t) + \sum_{i=1}^{N}  \varepsilon_i(t) H^i_{c}$, we will have $N$ different $\Phi_i$ with its corresponding sets of quadratic constraints. We also need to impose that $\Phi_i(t) \Phi_j^\dagger (t) = 0$ for all $i\ne j$.

Finally, note that if both $\varepsilon_1(t)$ and $\varepsilon_2(t)$ can be nonzero at the same time, it is equivalent to considering a problem with three controllable terms: $\varepsilon_1(t)H^1_c$, $\varepsilon_2(t)H^2_c$, and $\varepsilon_1(t)H^1_c + \varepsilon_2(t)H^2_c$, only one of which can be nonzero at any given time.

\section{Bound formulation for open systems governed by master equations}
In the main text, we only consider problems whose dynamics is governed by the Schr\"odinger equation. Now we show how to apply our bound framework to open systems described by differential equations for the density matrix $\rho$. In principle, this framework applies to any equation that is linear in $\rho$. Here we specifically consider the Lindblad master equation \cite{breuer2002theory}:

\begin{equation}
	\frac{\partial }{\partial t}\rho = -\frac{i}{\hbar}[H,\rho] + \sum_{i}\gamma_i \bb{A_i\rho A_i^\dagger - \frac{1}{2}\cb{A_i^\dagger A_i, \rho }}\,,
\end{equation}
where we also assume the Hamiltonian has the form $H = H_0 + \varepsilon(t)H_c$, and the second term is the dissipative part. A common trick is to vectorize the equation to form ``superoperators''~\cite{breuer2002theory},  through flattening the $n^2$ by $n^2$ density matrix $\rho$ into a column vector $\vec{r}$ by ordering the columns of $\rho$ in a single vector~\cite{horn2012matrix}. Then all operators are transformed according to this new ordering (basis):
\begin{eqnarray}
	&&[H,\rho] \longrightarrow \bb{I\otimes H - H^T\otimes I}\vec{r}\nonumber\\
	&&\quad\quad\quad\quad =  \bb{I\otimes H_0 - H_0^T\otimes I}\vec{r} +  \bb{I\otimes H_c - H_c^T\otimes I}\vec{r}\\
	&& A_i \rho A_i^\dagger \longrightarrow  \bb{(A_i^\dagger)^T \otimes A_i}\vec{r} \\
     && A_i A_i^\dagger \rho  \longrightarrow  \bb{I \otimes A_i^\dagger A_i}\vec{r} \\
     && \rho A_i^\dagger  A_i \longrightarrow  \bb{(A_i^\dagger A_i)^T \otimes I}\vec{r}\,.
\end{eqnarray}
If we define the vectorized  $-\frac{i}{\hbar}[H_0,\rho]$ by $L_0$, the vectorized dissipative term by $L_D$, and vectorized  $-\frac{i}{\hbar}[H_c,\rho]$ by $L_c$, we turn the original equation into
\begin{equation}
	\frac{\partial }{\partial t}\vec{r} = (L_0+ L_D)\vec{r} + L_c \vec{r}\,.
\end{equation}
This equation again can be turned into the integral form
\begin{equation}
	\vec{r}(t) = \vec{r}_0(t) + \int_{t_0}^T \D t' G^+(t,t')L_c(t’) \vec{r}(t’),
	\label{eq:int_r}
\end{equation}
where $\vec{r_0}(t)$ is the evolution of $\vec{r}$ without the presence of $L_c$, and $G^+(t,t')$ is the retarded Green's function of $\frac{\partial }{\partial t}\vec{r} =(L_0+ L_D)\vec{r}$. The dissipative nature of the problem is encoded in $G^+(t,t')$.

\eqref{int_r} is the starting point of the bound formulation, just as Eq. (1) in the main text for the Schr\"odinger equation. Carrying out the same manipulations as in the main text gives the constraints
\begin{eqnarray}
	&&\int_{t_0}^{T} \D t \,\vec{r}^\dagger(t) L_c (t) D_i(t) \vec{r}(t) + \int_{t_0}^{T} \D t \int_{t_0}^{T} \D t' \, \vec{r}^\dagger(t) L_c (t) D_i(t) G^+(t,t') L_c (t') \vec{r}(t') \nonumber \\
	&=&   \int_{t_0}^{T} \D t \vec{r}^\dagger(t) L_c (t) D_i(t) \vec{r}(t)\,,
\end{eqnarray}
which are analogous to Eq. (2) in the main text. The quadratic constraints can then be formed analogously.

\section{Numerical details}
In this section, we provide details for the numerical scheme used in the main text. Unlike conventional time-stepping methods to solve initial-value differential equations, we directly use the integral equation (Eq. (1) of the main text) to solve for the dynamics. We do this for three reasons: (1) In our bound method, constraints come from the integral equation, so the discretized form of the integral equation is directly available after solving the bound problem. (2) Pulse design is straightforward to perform with methods based on the integral equation. (3) Methods based on the integral equation are equally viable as those based on the differential equation.

\subsection{Nystr\"om discretization}
We consider the integral equation in Eq. (1). We consider systems with $L$ levels, modeling either discrete levels or a subset of continuous states, over times $[t_0,T]$. Here we isolate a single column of the time-evolution operator, essentially selecting an initial state $\ket{\psi_0}$ and the time-dependent state $\ket{\psi(t)} = U(t,t_0) \ket{\psi_0}$. Then the integral equation for $\ket{\psi(t)}$ is
\begin{equation}
	\ket{\psi(t)}= \ket{\psi_0(t)} -\frac{\I}{\hbar} \int_{t_0}^{t} U_0(t,t')H_c \varepsilon(t')\ket{\psi(t')}\D t'\,,
	\label{eq:VIE}
\end{equation}
where we have replaced the retarded Green's function with the time-evolution operator multiplied by a step function (changing the upper limit of the integral to $t$). This is in the form of a Volterra integral equation of the second kind.

The Nystr\"om method \cite{atkinson1967numerical} for discretizing this equation consistent of defining quadrature nodes $\{t_i\}_{i=0}^N\in[t_0,T]$ and corresponding weights\,, and enforcing
\eqref{VIE} to hold at $\{t_i\}_{i=0}^N$
\begin{equation}
	\ket{\psi(t_i)} = \ket{\psi_0(t_i)} -\frac{\I}{\hbar} \int_{t_0}^{t_i} U_0(t_i,t')H_c \varepsilon(t') \ket{\psi(t')} \D t'\,, \quad i = 0,1,2,\ldots, N\,,
\end{equation}
and then approximating the integral operator by
\begin{equation}
	\frac{\I}{\hbar} \int_{t_0}^{t_i} U_0(t_i,t')H_c \varepsilon(t')\ket{ \psi(t')} \D t' \approx \sum_{j = 1}^{t_i} w_{ij} \equals{\frac{\I}{\hbar}U_0(t_i,t_j)H_c \varepsilon(t_j)}{K_{ij}} \ket{\psi(t_j)}\,.
	 \label{eq:quad}
\end{equation}
for some weights $w_{ij}$. Then \eqref{VIE} becomes a linear system
\begin{equation}
	 \ket{\psi(t_i)} = \psi_0(t_i) -\sum_{j = 0}^{i} w_{ij} K_{ij}\ket{\psi(t_j)}\,, \quad i = 0,1,2,\ldots, N\,,
	 \label{eq:pre_linsys}
\end{equation}
or  more compactly as
\begin{equation}
	\vec{\psi}  = \vec{\psi}_0 - \vec{A}\vec{\psi}\,,
	\label{eq:linsys}
\end{equation}
where $\vec{\psi}$ and $\vec{\psi}_0$ are column vector of length $L(N+1)$ with the $i$th column block $\vec{\psi}_i = \ket{\psi(t_i)}$ and $(\vec{\psi}_0)_i = \ket{\psi_0(t_i)}$ respectively. And $\vec{A}$ is a matrix of size $L(N+1)\times L(N+1)$ whose explicit form we will show later. Then $\vec{\psi}$ can be solved as $(\vec{I} + \vec{A})^{-1}\vec{\psi}_0$\,.

We note that the convergence rate of $\vec{\psi}$ to the exact value $\psi$ is the same as the order of quadrature used in \eqref{quad}. The matrix norm of $\vec{A}$ is provably small, so $\vec{I} + \vec{A}$ is always well-conditioned. To conclude, the scheme we described above is stable and convergences quickly for good choice of $\{t_i\}_{i}^N$ and weights\,.

In this work, we use the trapezoidal rule
\begin{equation}
	\int_{a}^{b} f(x) \D x \approx \frac{h}{2} \bb{f(a) +2f(a+h) + \ldots + 2f(b-h)+f(b)} + O(h^2)
	\label{eq:trap}
\end{equation}
where $h$ is the spacing between the equally spaced nodes from $a$ to $b$. As a result, the method we are using has convergence rate equivalent to Runge Kutta of order 2. Without special techniques for the matrix inversion process, this method is not as efficient as Runge Kutta of order 2. Despite of this drawback and the fact that higher order methods are available, the method here is sufficient for the examples we consider. If we apply the trapezoidal rule to \eqref{pre_linsys}, then the nodes $\{t_i\}_{i=0}^N\in[t_0,T]$ are $t_0, t_1,\ldots,t_N = T$ with $t_i = t_0 + ih$ and $h = \frac{T-t_0}{N}$\,. \eqref{pre_linsys} then becomes
\begin{eqnarray}
	\ket{\psi(t_0)} &=& \ket{\psi_0(t_0)} \nonumber\\
	\ket{\psi(t_1)} &=& \ket{\psi_0(t_1)} - \frac{h}{2}(K_{10} \ket{\psi(t_0)}+ K_{11}\ket{\psi(t_1)} ) \nonumber\\
	\ket{\psi(t_2)} &=& \ket{\psi_0(t_2)} - \frac{h}{2}(K_{20}\ket{ \psi(t_0)} + 2K_{21}\ket{\psi(t_1)} + K_{22} \ket{\psi(t_2)}  ) \nonumber\\
	&\vdots& \nonumber\\
	\ket{\psi(t_N)} &=& \ket{\psi_0(t_N)} - \frac{h}{2}(K_{N0} \ket{\psi(t_0)} + 2\sum_{j =1}^{N-1}K_{Nj}\ket{\psi(t_j)} + K_{NN}\ket{\psi(t_N)}  )\nonumber\,,\\
	\label{eq:linsys_exp}
\end{eqnarray}
which is \eqref{linsys} in component form with the trapezoidal rule.

To solve for $U$\,, we simply replace  the vector $\vec{\psi}$ by the matrix $\vec{U} = \bb{\mat{\vec{\psi}^1 & \vec{\psi}^2 & \ldots} \vec{\psi}^L }$ and $\vec{\psi}_0$ by the matrix $\vec{U}_0 = \bb{\mat{\vec{\psi}_0^1 & \vec{\psi}_0^2 & \ldots \vec{\psi}_0^L}}$\,, both of size $L(N+1)\times L$\,, where $\vec{\psi}^i $ represents the discretized version of  different states wave function $\ket{\psi^i (t)} = U(t,t_0)\ket{\psi_0^i (t_0)}$.
By comparing \eqref{linsys_exp} and \eqref{pre_linsys}, we can see that the weight $w_{ij} = h$ for $i \ne j$ and $w_{ij} = h/2$ for $i=j$ and $w_{ij} = 0$ for $i<j$.
Then the discretized version of Eq. (1) in the main text is given by
\begin{equation}
	\vec{U}  = \vec{U}_0 - \vec{A}\vec{U}\,,
	\label{eq:linsys_U}
\end{equation}
and $\vec{U} = (\vec{I} + \vec{A})^{-1}\vec{U}_0$\,.
\subsection{Discretization for quadratic constraints}
With the understanding of how the integral equation is discretized, we are ready to give the discretization form of Eq. (3), which constitutes the constraints for computing the bound.
Constraints in Eq. (3) are given by
\begin{eqnarray}
	&& \int_{t_0}^{T}\D t \int_{t_0}^{T}\D t' \,  \Phi^\dagger (t)D_i(t)\bb{ \delta(t-t')\frac{H^{-1}_c}{\varepsilon_{\mathrm{max}}} + \frac{i}{\hbar}G_0^+(t,t')}  \Phi(t')
	\nonumber\\
	&& = \int_{t_0}^{T}\D t \,  \Phi^\dagger(t)D_i(t) U_0(t,t_0) \,,
\end{eqnarray}
or
\begin{eqnarray}
		&& \int_{t_0}^{T}\D t \int_{t_0}^{t}\D t' \,  \Phi^\dagger (t)D_i(t) \frac{i}{\hbar} U_0(t,t') \Phi(t') +  \int_{t_0}^{T}\D t \, \Phi^\dagger (t)D_i(t)\frac{H^{-1}_c}{\varepsilon_{\mathrm{max}}}   \Phi(t)
	\nonumber\\
	&& = \int_{t_0}^{T}\D t \,  \Phi^\dagger(t)D_i(t) U_0(t,t_0) \,,
	\label{eq:pre_dis_con}
\end{eqnarray}
We have seen in the previous section on how to discretize the $t'$ integral. The $t$ integral is in fact much easier with fixed integration limit, and we can apply the trapezoidal rule directly to the integral $\int_{t_0}^T \D t f(t) \approx \sum_{i=0}^{N} w_i f(t_i)$ with the same equally spaced nodes $\{t_i \}_{i=0}^N$ and weights $w_i$ in \eqref{trap}\,.
Then \eqref{pre_dis_con} becomes
\begin{equation}
	\vec{\Phi}^\dagger \vec{D}_i \vec{W}\vec{G}\vec{\Phi} + \vec{\Phi}^\dagger \vec{D}_i \vec{W}\frac{\vec{H}^{-1}_c}{\varepsilon_{\rm max}}\vec{\Phi} = \vec{\Phi}^\dagger \vec{D}_i \vec{W}\vec{U}_0\,,
	\label{eq:QC_con}
\end{equation}
where $\vec{\Phi}$ and $\vec{U}_0$ are matrix of size $L(N+1)\times L$.
$\vec{G}$ is the discretized version of the integral operator $\int_{t_0}^t \D t' U(t,t')$. Similar to the way we discretized \eqref{quad}\,, we have the $ij$ block of $\vec{G}$
\begin{equation}
	\vec{G}_{ij} = w_{ij} \frac{\I}{\hbar} U_0(t_i,t_j)
	\label{eq:G}
\end{equation}
with $w_{ij}$ the weight in \eqref{linsys_exp}.
$\vec{W}$ is the matrix contains the effect of weights $w_i$ from the $t$ integral, and it is of the from
\begin{equation}
	\vec{W} = \bb{\mat{w_0 &&&\\&w_2&&\\ &&\ddots&\\&&&w_N}}\otimes \vec{I}_{L\times L}\,,
\end{equation}
where $\otimes$ is the kroncker product for matrices and $ \vec{I}_{L\times L}$ represents an identity matrix of dimension $L$.
Here $\vec{H}^{-1}_c$ is
\begin{equation}
	\vec{H}^{-1}_c = \vec{I}_{(N+1)\times(N+1)}\otimes H^{-1}_c
\end{equation}
Also, since $D_i = \delta(t-t_i)D$ and $D$ is some constant matrix of dimension $L$ by $L$, we have
\begin{equation}
	\vec{D}_i =  \bb{\mat{0 &&&&\\&\ddots&&&\\ &&1 & &\\&&&\ddots&\\&&&&0}}\otimes D\,,
\end{equation}
where the first diagonal matrix is the discretized $\delta(t-t_i)$ and is nonzero only in the $i$th diagonal entry.
$D$ of dimension $L$ by $L$ consists of $L^2$ unit basis matrices. For a two level system, all possible independent $D$ consists of $\bb{\mat{1 &0\\ 0& 0}}$, $\bb{\mat{0 &1\\ 0& 0}}$, $\bb{\mat{0 &0\\1& 0}}$ and $\bb{\mat{0 &0\\0& 1}}$\,. Higher dimension systems can be generalized accordingly.

By varying $D_i$ over all possible $t_i$ and $D$, we obtain all possible constraints in our problem in the form of a $L$ by $L$ matrix equality, which must hold componentwise.
$\vec{\Phi}$ can be flattened column by column into a vector, and the equality constraints for each components can then be transformed  into a quadratic constraint of the flattened $\vec{\Phi}$.
Then semidefinite relaxation can be carried out to form the semidefinite program, which can be solved by standard convex optimization package.

\subsection{Pulse designs via local optimization}
In this section, we describe how to efficiently compute the gradient  $\frac{\partial f(U)}{\partial \varepsilon(t)}$ of the objective $f(U)$ for all $t$ from $[t_0, T]$ by solving for the dynamics twice only. Now, we return to \eqref{linsys_U}, the discretized version of Eq. (1) in the main text. To make the variables  $\varepsilon(t_i)$ here more prominent, we separate  $\varepsilon(t_i)$ from $\vec{A}$ in \eqref{linsys_U} and write it as
\begin{equation}
	\vec{U} = \vec{U}_0 - \equals{\vec{G}\vec{H}_c \vec{\varepsilon}}{\vec{A}}\vec{U}\,,
	\label{eq:VIE_design}
\end{equation}
where
\begin{equation}
	\vec{H}_c = \vec{I}_{(N+1)\times(N+1)}\otimes H_c\,,
\end{equation}
\begin{equation}
	\vec{\varepsilon} =  \bb{\mat{\varepsilon(t_0) &&&\\&\varepsilon(t_1) &&\\ &&\ddots&\\&&&\varepsilon(t_N) }}\otimes \vec{I}_{L\times L}\,,
\end{equation}
and $\vec{G}$ is given in \eqref{G}. We also denote the column vector $\bb{\mat{\varepsilon(t_0) & \varepsilon(t_1)&\ldots&\varepsilon(t_N)}}^T$ by $\vec{p}$\,.

The goal here is to compute the discretized gradient $\frac{\partial f(\vec{U})}{\partial \vec{p}}$ in a way that the work required is independent of the number of variable $\varepsilon(t_i)$\,. To do so, we first use the CR calculus to formally treat $\vec{U}$ and $\vec{U}^*$ as independent variables \cite{kreutz2009complex}; the chain rule yields
\begin{equation}
		\frac{\partial f(\vec{U})}{\partial \vec{p}} = \frac{\partial f(\vec{U})}{\partial \vec{U}} \frac{\partial \vec{U}}{\partial \vec{p}} + \frac{\partial f(\vec{U})}{\partial \vec{U}^*} \frac{\partial \vec{U}^*}{\partial \vec{p}}  = 2\Re{\frac{\partial f(\vec{U})}{\partial \vec{U}} \frac{\partial \vec{U}}{\partial \vec{p}} }\,.
		\label{eq:grad}
\end{equation}
Then we differentiate \eqref{VIE_design} with respect to $\vec{p}$ and obtain
\begin{equation}
	\vec{G}\vec{H}_c \frac{\partial \vec{\varepsilon}}{\partial \vec{p}}\vec{U} = -(\vec{I} + \vec{A})\frac{\partial \vec{U}}{\partial \vec{p}}\,,
\end{equation}
and we can solve for $\frac{\partial \vec{U}}{\partial \vec{p}} = -(\vec{I} + \vec{A})^{-1}\vec{G}\vec{H}_c \frac{\partial \vec{\varepsilon}}{\partial \vec{p}}\vec{U} $\,. Substituting this into \eqref{grad} gives
\begin{equation}
		\frac{\partial f(\vec{U})}{\partial \vec{p}}  = -2\Re{\frac{\partial f(\vec{U})}{\partial \vec{U}}(\vec{I} + \vec{A})^{-1}\vec{G}\vec{H}_c \frac{\partial \vec{\varepsilon}}{\partial \vec{p}}\vec{U} }.
\end{equation}
Since $\vec{G}$, $\vec{H}_c$ are known $\frac{\partial \vec{\varepsilon}}{\partial \vec{p}}$ can be calculated analytically, then if we define the so-called ``adjoint'' solution \cite{strang2007computational,miller2013photonic}
\begin{equation}
    (\vec{I} + \vec{A})^T\vec{U}_{\rm adj} =\bb{ \frac{\partial f(\vec{U})}{\partial \vec{U}}}^T\,,\
	\label{eq:adjoint}
\end{equation}
we can obtain the gradient as
\begin{equation}
    \frac{\partial f(\vec{U})}{\partial \vec{p}} = -2\Re{\vec{U}_{\rm adj}^T \vec{G}\vec{H}_c \frac{\partial \vec{\varepsilon}}{\partial \vec{p}}\vec{U} }\,,
\end{equation}
by solving for $\vec{U}_{\rm adj} = (\vec{I} + \vec{A}^T)^{-1}\bb{ \frac{\partial f(\vec{U})}{\partial \vec{U}}}^T$  from \eqref{adjoint} and $\vec{U} = (\vec{I}+\vec{A})^{-1} \vec{U}_0$ from \eqref{VIE_design} only.
This allows rapid computation of $\frac{\partial f(\vec{U})}{\partial \vec{p}}$. Then standard optimization methods, such as gradient ascent, can be applied to maximize or minimize $f(\vec{U})$\,. All designed pulses in the main text are obtained by this approach described here.

\section{Guide to computing prior-literature bounds}
In this section, we describe how we computed the bounds from prior literature in Fig.~1 of the main text: Mandelstam--Tamm (MT), Margolus--Levitin (ML), and the bounds from Refs.~\cite{Arenz2017,Lee2018}. We start with the MT bound: the minimum time, which we denote $\tau_{\rm MT}$, for a given evolution is given by
\begin{align}
    \tau_{\rm MT} = \frac{\pi}{2\Delta H},
\end{align}
where the standard deviation $\Delta H$ is given by $\Delta H = \sqrt{\langle H^2 \rangle - \langle H \rangle^2}$. We do not know the mean of the Hamiltonian over the optimal trajectory, but can bound its square below by zero, which gives $\Delta H \leq \sqrt{\langle H^2 \rangle}$. We don't know the value of $\langle H^2 \rangle$ over the optimal trajectory either, but it can be bounded above by the square of the largest possible eigenvalue of the Hamiltonian, which we denote $\lambda_{\rm max}$. This eigenvalue can be found by a grid-search-based optimization over all possible instantaneous Hamiltonians (assuming bounded controls), and we are left with the bound:
\begin{align}
    \tau_{\rm MT} \geq \frac{\pi}{2\lambda_{\rm max}}.
\end{align}
The ML bound, which we denote $\tau_{\rm ML}$, follows similar reasoning. Now given by $\tau_{\rm ML} = \pi / 2\langle H \rangle$, we again do not know the mean energy, but can now bound it above by $\lambda_{\rm max}$, which implies that
\begin{align}
    \tau_{\rm ML} \geq \frac{\pi}{2\lambda_{\rm max}}.
\end{align}
By this reasoning, the MT and ML bounds coincides when presented with uncertainty about the optimal trajectory and bounded in this way.

Next, we consider the bound of \citeasnoun{Arenz2017}. In this case, one defines two functions incorporating information about the controllable and non-controllable Hamiltonians, $H_c$ and $H_0$, respectively:
\begin{align}
    C(U_g,H_c) = \frac{\sqrt{2\left(d - \sum_{j=1}^d \left| \langle \phi_j^{(c)} | U_g | \phi_j^{(c)} \rangle \right| \right)}}{2}, \\
    C(U_g,H_0) = \frac{\sqrt{2\left(d - \sum_{j=1}^d \left| \langle \phi_j^{(0)} | U_g | \phi_j^{(0)} \rangle \right| \right)}}{2},
\end{align}
where $U_g$ is a target evolution operator, $d$ is the dimensionality, and $\phi_j^{(c)}$ and $\phi_j^{(0)}$ are the eigenfunctions of the controllable and non-controllable Hamiltonians, respectively. From these functions, the time bounds of \citeasnoun{Arenz2017} can be written:
\begin{align}
    T \geq {\rm max} \left\{ \frac{2C(U_g,H_c)}{\|H_0\|}, \frac{2C(U_g,H_0)}{|f_{\rm max}|\|H_c\|} \right\},
    \label{eq:Tmax}
\end{align}
where $|f_{\rm max}|$ is the maximum amplitude of the bounded control parameter. \Eqref{Tmax} is straightforward to compute, though hard to apply to the maximum transition probability in a multilevel system, as there is no specific target unitary evolution. In this case, to compute a meaningful bound at all (for the sake of comparison), we instead consider the best numerically-optimized control sequence, and use its evolution operator evaluated at the final time as the target unitary. In this way, we take a known achievable propagator and can generate a meaningful bound to compare it to. It should be noted that if one were to want to derive a strict bound in this way, it would be required to find the smallest bound over all possible unitaries that correspond to unity probability in the desired state; hence, our computation of a bound is actually an over-estimate of the ``true'' bound using this approach, which is necessarily smaller (and thereby looser).

The bound of \citeasnoun{Lee2018} is similarly for a target unitary, but now uses a different expression. Now, the minimum time is given by:
\begin{align}
    T \geq \max_{V\in {\rm Stab}(iH_c)} \frac{\|[U_g,V]\|}{\|[H_0,V]\|}.
    \label{eq:Tmax2}
\end{align}
This bound arises from maximizing the ratio of norms on the right-hand side of \eqref{Tmax2} over all matrices that commute with $H_c$. Generically, there is no simple way to optimize over all such matrices. To compute a reasonable bound, then, we choose $H_c$ for $V$, as it is guaranteed to commute with itself. Then we evaluate the ratio in \eqref{Tmax2}, again using a single unitary as described above (instead of minimizing over all possible unitaries). In both \eqref{Tmax} and \eqref{Tmax2}, the Frobenius norm is the norm that is used.