* Three calls of the unitary channels A and one call of unitary channels B, in the order AAAB.

In order to show these results, we numerically evaluate the maximum probability of success of simulating the quantum switch when the input channels are unitary. We do so in the case where the input channels are acting on qubit systems, i.e., when $d = 2$, and in the restricted simulation scenario where the input control and target systems are fixed, and the output target system is discarded [see Eq. (D1)]. We remark once again that the impossibility of a deterministic simulation in the restricted, fixed-dimension case implies the impossibility of a deterministic simulation in general.

In this case, we do not construct an explicit basis analytically for the subspace spanned by $k$ identical copies of a qubit unitary channel, but instead, we randomly sample a set of $d_I$ linearly independent qubit unitary channels $\{U_i^{\otimes k}\}_{i=1}^{d_I}$ to form our basis. We repeat here that for the case where $k = 1$, $d_I = 10$, for $k = 2$, $d_I = 35$, and for $k = 3$, and $d_I = 84$. This implies that in the $(k_A, k_B) = (2, 1)$ case, one needs $35 \cdot 10 = 350$ pairs of qubit unitary channels and for $(k_A, k_B) = (3, 1)$, one needs $84 \cdot 10 = 840$ pairs of qubit unitary channels as input for the primal and dual SDPs [Eqs. (27) and (28) of the main text]. We numerically obtain the values of

$$ \text{order} = \text{AAB (qubit, unitary), } p \approx 0.600 \qquad (\text{E1}) $$

$$ \text{order} = \text{BAA (qubit, unitary), } p \approx 0.851 \qquad (\text{E2}) $$

$$ \text{order} = \text{AAAB (qubit, unitary), } p \approx 0.708. \qquad (\text{E3}) $$

## Appendix F: Efficient certification that a matrix is positive semidefinite

A self-adjoint linear operator $A \in \mathcal{L}(\mathbb{C}^d)$ is positive semidefinite if and only if there exists an operator $L \in \mathcal{L}(\mathbb{C}^d)$ such that $A = LL^\dagger$. Moreover, the operator $L$ can be taken to be a lower triangular matrix, i.e., a matrix in which all entries above the main diagonal are zero. The decomposition of a positive semidefinite matrix as $A = LL^\dagger$ is referred to as the Cholesky decomposition, and finding such a decomposition can be done efficiently [66].

In the final step of the algorithm used for computer-assisted proofs presented in the Methods section of the main text, one is required to certify that an operator $A$ is positive semidefinite. When $A$ is stored as a symbolic matrix, due to the way computers manipulate symbolic variables, if the matrix $A$ is not sparse enough, ensuring that $A \ge 0$ may be a prohibitively time-consuming task even when using the Cholesky decomposition. This was the case, for instance, when trying to ensure that the matrix $\overline{\Gamma}(\Gamma^{\text{sym}}) + \eta \mathbb{1} - \sum_{i,j} R_{ij}^{\text{OK}} \otimes (J_i^{A^{\otimes k_A}} \otimes J_j^{B^{\otimes k_B}})^T$ in our computer-assisted proofs algorithm involving $k_A + k_B = 4$-slot quantum combs. Below, we describe an algorithm that can be used to ensure that a symbolic matrix $A$ is positive semidefinite, which is considerably faster than performing Cholesky decomposition on a symbolic matrix $A$.

The algorithm we present below is based on three key ideas.

1. It is possible to rigorously certify that a matrix is positive definite using floating-point arithmetic quickly. One way to attain this goal is to use the methods presented in Ref. [67]. It shows how to efficiently certify that a matrix $A$ is positive definite using a rigorous algorithm that accounts for all possible computational and rounding errors and remains valid in the presence of underflow.

2. If a matrix $A$ is “close” to another matrix $A'$, and $A'$ is “far” from the set of non-positive semidefinite matrices, then $A$ has to be positive semidefinite.

3. If $A$ is matrix with symbolic entries, we can obtain a floating-point variable matrix $A'$ that is guaranteed to be close to $A$. This can be done via arbitrary-precision arithmetic [68].

### Algorithm to prove that a symbolic self-adjoint matrix $A$ is positive semidefinite:

1. Construct a self-adjoint matrix $A'$ that is equal to $A$ up to $n$ decimal digits. This step can be accomplished using arbitrary-precision arithmetic [68] with a precision of $n$ decimal digits. In this way, all matrix elements of $A' - A$ are between $-10^{-n+1}$ and $10^{-n+1}$, ensuring that

$$ -J \cdot 10^{-n+1} \le A' - A \le J \cdot 10^{-n+1}, \qquad (\text{F1}) $$

where $J \ge 0$ is a matrix in which all entries are the number one.
2. Prove that $A' - J \geq 0$ using the algorithm presented in Ref. [67].
Since $A' - J \cdot 10^{-n+1} \leq A$ holds, if we ensure that $A' - J \geq 0$, then by transitivity it follows that $A \geq 0$.