2.1. Classical model

Laser amplification is often treated as a parametric process, where the signal and idler waves grow by consuming a pump wave. When the pump energy is being replenished, or when the pump energy dominates, one may approximate the pump amplitude $a_1$ as a constant, in which case the seed amplitude $a_2$ and the idler amplitude $a_3$ grow exponentially. However, when the pump amplitude is not held constant, the three-wave nature of the underlying interaction becomes apparent. The interaction is often described by the three-wave equations (Davidson Reference Davidson2012)

(2.1a)\begin{gather} d_t a_1 = g a_2 a_3, \end{gather}

(2.1b)\begin{gather}d_t a_2 ={-}g^* a_1 a_3^{\dagger}, \end{gather}

(2.1c)\begin{gather}d_t a_3 ={-}g^* a_1 a_2^{\dagger}, \end{gather}

where $d_t$ is the advective derivative, $g$ is the coupling coefficient, $g^*$ is its complex conjugate and $a^{\dagger}$ denotes the complex conjugate of $a$ in the classical model. The advective derivative is specific for each wave and is defined as $d_t=\partial _t + \boldsymbol {v}_g\boldsymbol {\cdot }\boldsymbol {\nabla } + \mu$, where $\boldsymbol {v}_g=\partial \omega /\partial \boldsymbol {k}$ is the group velocity and $\mu$ is the damping rate of the wave. The complex-valued amplitude $a$ is the slowly varying envelope of the classical wave. The amplitude is normalized such that $n=|a|^2$ is proportional to the wave action density, which is proportional to the number of photons in the wave.

Three-wave interactions satisfy a number of conservation laws. First, the equations describe resonant interactions where both energy and momentum are conserved. For waves with narrow bandwidth, their 4-momentum density is proportional to $k^\mu =({\omega }/{c},\boldsymbol {k})$. Energy and momentum are conserved because the interaction satisfies the resonance conditions $k_1^\mu =k_2^\mu +k_3^\mu$. Moreover, the interaction has two independently conserved actions $S_2=n_1+n_2$ and $S_3=n_1+n_3$. These actions are constants of motion because their advective derivatives are zero. The physical meaning of action conservation is that the three-wave interaction splits each pump photon into a seed and an idler photon. So, whenever $n_1$ decreases by one, both $n_2$ and $n_3$ increase by one, and vice versa.

The three-wave equations are partial differential equations that describe how the wave amplitudes evolve in both space and time. If all amplitudes are initially uniform in space, then they will remain uniform as they evolve in time. However, if the amplitudes are not uniform, two competing effects change the envelopes of the waves. First, wave advection transports the wave envelope in the direction of the group velocity. The advection is a linear process, and the envelope remains unchanged in the co-moving frame. Second, the three-wave interaction changes amplitudes locally. Since the interaction is nonlinear, the change is faster where the amplitudes of the other two waves are larger.

Laser pulse compression is a special scenario where the competition between the two effects leads to the amplification and shortening of the seed wave. At later stages of pulse compression, the intensity of the seed far exceeds that of the pump. The large $a_2$ induces an additional relativistic nonlinearity. The nonlinearity originates from the fact that, in plasmas, photons are massive particles due to their interactions with free charges. In unmagnetized plasmas, photons satisfy the dispersion relation $\omega ^2=\omega _p^2+c^2k^2$, where the photon mass can be identified with the plasma frequency $\omega _p=(e^2n_e/\epsilon _0 m_e)^{1/2}$. Here, $e$ is the electron charge and $n_e$ is the electron density. Because electrons oscillate in the laser's electric field, the effective electron mass $m_e$ is replaced by $\gamma m_e$ when the electron quiver speed $v_q$ becomes comparable to the speed of light $c$, where $\gamma =1/\sqrt {1-v_q^2/c^2}$. As the seed pulse propagates, at places where the pulse is more intense, the photon mass $\omega _p \propto \gamma ^{-1/2}$ becomes smaller. A smaller $\omega _p$ leads to a larger $k$ at a fixed $\omega$, which means a larger group velocity $\boldsymbol {v}_g=c^2\boldsymbol {k}/\omega$. Consequently, the more intense part of $a_2$ moves at a higher group velocity. If the envelope of $a_2$ has initial modulations, then they will pile up and grow. This process is known as relativistic modulational instability.

The modulational instability may be understood as a four-wave process, where the laser beats with its modulations to produce side bands. In the weakly relativistic limit, we can expand $1/\gamma \simeq 1-v_q^2/(2c^2)$. Since electrons respond primarily to the electric field of the laser pulse, Newton's equation $m_e \,{\rm d}\boldsymbol {v}_q/{\rm d} t\simeq e\boldsymbol {E}$ becomes $\tilde {\boldsymbol {v}}_q\simeq \textrm {i}e\tilde {\boldsymbol {E}}/(m_e\omega )$ in Fourier space. Denoting the normalized laser amplitude by $\boldsymbol {\alpha }=e\tilde {\boldsymbol {E}}/(m_e\omega c)$, then the average quiver speed is $\langle v_q^2/c^2\rangle =\frac {1}{2} e^2|\tilde {E}|^2/(m_e\omega c)^2= \frac {1}{2}|\alpha |^2$. Replacing $\omega _p^2\rightarrow \omega _p^2/\gamma \simeq \omega _p^2(1-\frac {1}{4}|\alpha |^2)$, the photon dispersion relation is approximated. The dispersion relation is derivable from the wave equation $[\partial _t^2-c^2\nabla ^2+\omega _p^2(1-\frac {1}{4}|\alpha |^2)]\boldsymbol {E}=0$. In the Wentzel–Kramers–Brillouin approximation, the complex wave is $\boldsymbol {E}(x,t)=\boldsymbol {A}(x,t)\,{\rm e}^{\textrm {i}\theta }$, where $\theta =\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {x}-\omega t$, and the envelope $\boldsymbol {A}$ varies slowly in the sense that $|\partial _t\boldsymbol {A}/\boldsymbol {A}|\ll \omega$ and $|\boldsymbol {\nabla } \boldsymbol {A}/\boldsymbol {A}|\ll k$. Then, to leading order, the wave equation is approximated by $[\partial _t + \boldsymbol {v}_g\boldsymbol {\cdot }\boldsymbol {\nabla } -\textrm {i}\omega _p^2|\alpha |^2/(8\omega )] \boldsymbol {A}=0$. When focusing on the scalar amplitude $A$, because $A\propto \alpha \propto a$, the equation that describes the modulational instability of the seed pulse is

(2.2)\begin{equation} d_t a_2 = \textrm{i} R a_2^{\dagger} a_2 a_2, \end{equation}

where $d_t$ is again the advective derivative and $R=\omega _p^2/(8\omega$) is the coupling coefficient. For relativistic modulational instability, $R>0$ is a real number, which means ${\rm i}R |a_2|^2$ is purely imaginary. The above equation thus modulates the phase of the complex $a_2$ in such a way that a larger $|a_2|$ leads to a faster phase evolution.

Laser pulse compression is described by the combined equations (2.1) and (2.2). The modulational instability only affects $a_2$ because during laser pulse compression, the peak values satisfy $|a_2|\gg |a_1|, |a_3|$. Since the three-wave interaction is a phase-sensitive process, the modulational instability spoils the amplification process by introducing a phase mismatch.

2.2. Quantum model

In the classical model, the amplitude $a$ is a complex-valued function, and $n=|a|^2$ is proportional to the number of photons. This set-up naturally admits canonical quantization for bosonic quantum fields $[a_i(\boldsymbol {x}), a_j^{\dagger} (\boldsymbol {y})] = \delta _{ij}(2{\rm \pi} )^3\delta ^{(3)}(\boldsymbol {x}-\boldsymbol {y})$, where the indices $i, j=1,2,3$, and the operators have spatial dependencies. Since we will later implement the model on quantum hardware, which has a limited number of qubits, in this paper we focus on the temporal problem with no spatial dependence. In this case, when damping is negligible, the advective derivative $d_t\rightarrow \partial _t$ is reduced to a partial derivative in time, and the operators satisfy the canonical quantization conditions

(2.3)\begin{equation} [a_i, a_j^{\dagger}] = \delta_{ij}. \end{equation}

The quantization promotes normalized amplitudes to creation and annihilation operators, and the Kronecker delta distinguishes the three types of waves in the system. For each wave type, the number operator is $n_i=a_i^{\dagger} a_i$. The eigenstates of $n_i$ are the Fock states $|m_i\rangle$, namely $n_i|m_i\rangle = m_i|m_i\rangle$, where $m_i=0,1,2\dots$ and $|m_i\rangle =(a_i^{\dagger} )^{m_i}|0_i\rangle /\sqrt {m_i!}$. Here, $|0_i\rangle$ is the ground state of wave $i$, which is annihilated by $a_i|0_i\rangle =0$. More generally, these quantum harmonic oscillators have the usual matrix elements $a_i|m_i\rangle =\sqrt {m_i}|m_i-1\rangle$ and $a_i^{\dagger} |m_i\rangle =\sqrt {m_i+1}|m_i+1\rangle$. Since we have three types of waves, it is convenient to abbreviate the tensor-product state $|m_1\rangle \otimes |m_2\rangle \otimes |m_3\rangle$ as $|m_1, m_2, m_3\rangle$. This number basis is natural for the quantized problem.

While the Schrödinger picture $\textrm {i}\partial _t|\psi \rangle =H |\psi \rangle$ is more convenient for quantum simulations, the connection between the quantum and classical models is more transparent in the Heisenberg picture $d_ta=\textrm {i}[H, a]$. For three-wave interactions, equations (2.1) are the Heisenberg equations from the cubic Hamiltonian

(2.4)\begin{equation} H_T = \textrm{i}g a_1^{\dagger} a_2 a_3 - \textrm{i}g^* a_1 a_2^{\dagger} a_3^{\dagger}. \end{equation}

A degenerate form of $H_T$, where $a_1=a_2$, commonly arises in optomechanical systems (Aspelmeyer & Schwab Reference Aspelmeyer and Schwab2008; Kong et al. Reference Kong, Li, You, Xiong and Wu2018; Lake et al. Reference Lake, Mitchell, Sanders and Barclay2020; Shang Reference Shang2023). In our case, $a_1\ne a_2$. The first term of $H_T$ annihilates a seed and an idler photon to create a pump photon, while the second term of $H_T$ is the reverse process where a pump photon decays into a seed and an idler photon. Although the Heisenberg equations for $a_j$ are formally identical to the classical three-wave equations, the differences between the quantum and classical systems become apparent when one calculates higher-order cumulants. For example, because $a_j$ and $a^{\dagger} _j$ do not commute, the Heisenberg equation for $n_i$ is different from its classical counterpart (Shi et al. Reference Shi, Castelli, Wu, Joseph, Geyko, Graziani, Libby, Parker, Rosen and Martinez2021a). Similarly, for the four-wave interaction, (2.2) is the Heisenberg equation from the quartic Hamiltonian

(2.5)\begin{equation} H_F ={-}\frac{R}{2} a_2^{\dagger} a_2^{\dagger} a_2 a_2, \end{equation}

which is also known as the self-Kerr nonlinearity in the quantum literature (Kerr Reference Kerr1875). Since $R>0$ for the modulational instability, the negative sign in $H_F$ means that photons tend to condense together, which leads to a lower energy of the system. In (2.5), the operators are normal ordered $a_2^{\dagger} a_2^{\dagger} a_2 a_2=n_2^2-n_2$, which is different from other orderings such as $a_2^{\dagger} a_2 a_2^{\dagger} a_2=n_2^2$. Different orderings differ by factors of the number operator $n_2$. For our purposes, the ordering does not matter, because the total Hamiltonian also includes the energy of the non-interacting harmonic oscillators $H_0=\sum _j\omega _ja_j^{\dagger} a_j=\sum _j\omega _jn_j$. In the interaction picture of $H_0$, these quadratic terms are trivially removed, and (2.4) and (2.5) should be understood as three- and four-wave interactions in the interaction picture of $H_0$. Then, different orderings in (2.5) are equivalent up to a re-definition of $\omega _2$. Hence, it is sufficient to consider normal-ordered operators.

To use quantum Hamiltonian simulations to solve the quantized wave–wave interaction problems, we focus on the Schrödinger picture and use a basis that respects action conservation. In classical wave–wave interactions, $S_2=n_1+n_2$ and $S_3=n_1+n_3$ are known as the conserved wave actions. In the quantized model, $[H_T, S_2]=[H_T, S_3]=0$, and $H_F$ also commutes with $S_2$ and $S_3$. Therefore, it is convenient to use eigenstates of $S_2$ and $S_3$ as the computational basis, which we call the action basis. For the laser pulse compression problem, since we are primarily interested in the seed wave $a_2$, we label the action basis by

(2.6)\begin{equation} |\phi_j^{s_2, s_3}\rangle = |s_2 - j, j, s_3-s_2+j\rangle, \end{equation}

where $j$ is the number of photons in the seed wave and is bounded within the range $j_{\min } \le j\le s_2$ with $j_{\min }=\max (0, s_2-s_3)$. In other words, the $|\phi _j^{s_2, s_3}\rangle$ basis spans a $D=\min (s_2, s_3)+1$ dimensional subspace of $|m_1, m_2, m_3\rangle$, where $m_1=s_2-j$, $m_2=j$ and $m_3=s_3-s_2+j$. The bounds for $j$ come from the fact that $m_i\ge 0$. The non-negative integers $s_2$ and $s_3$ are eigenvalues of $S_2$ and $S_3$, namely $S_2|\phi _j^{s_2, s_3}\rangle =s_2 |\phi _j^{s_2, s_3}\rangle$ and $S_3|\phi _j^{s_2, s_3}\rangle =s_3 |\phi _j^{s_2, s_3}\rangle$. In other words, the infinite-dimensional Hilbert space can be decomposed as a direct sum of finite-dimensional subspaces, where each subspace is labelled by a pair of quantum numbers $(s_2, s_3)$.

In the action basis, the Hamiltonian becomes block diagonal. The total Hamiltonian that governs the mixed three- and four-wave interaction problem is

(2.7)\begin{equation} H = H_T +H_F, \end{equation}

where $[H_T, H_F]\ne 0$. The matrix elements of $H_T$ are $H_T |\phi _j^{s_2, s_3}\rangle = \textrm {i}g \eta _{j-{1}/{2}}^{s_2, s_3} |\phi _{j-1}^{s_2, s_3}\rangle - \textrm {i}g^* \eta _{j+{1}/{2}}^{s_2, s_3} |\phi _{j+1}^{s_2, s_3}\rangle$, where the reduced matrix element $\eta _{j-{1}/{2}}^{s_2, s_3}=\sqrt {(s_2+1-j)j(s_3-s_2+j)}$ is inherited from the creation and annihilation operators. Notice that $H_T$ only couples nearest neighbours in the action basis. Moreover, since the annihilation operator terminates at $m=0$, the reduced matrix element vanishes for both bottom and top values of $j$. The bottom values of $j$ are either $j=0$ or $j=s_2-s_3$, and in both cases $\eta _{-{1}/{2}}^{s_2, s_3}=\eta _{s_2-s_3-{1}/{2}}^{s_2, s_3}=0$. The top value is $j=s_2$, for which $\eta _{s_2+{1}/{2}}^{s_2, s_3}=0$. In other words, in the $|\phi _j^{s_2, s_3}\rangle$ basis, the matrix of $H_T$ is block tridiagonal, where each block is finite-dimensional, with no need for artificial truncation. Within each block, $H_F |\phi _j^{s_2, s_3}\rangle = -R\zeta _j |\phi _j^{s_2, s_3}\rangle$. The reduced matrix element is $\zeta _j=\frac {1}{2}j(j-1)$, which is non-zero only when $j\ge 2$. This is intuitive because the self-Kerr effect is a photon–photon nonlinearity, so at least two photons are needed to see the effect. Because the matrix for $H_F$ is diagonal, one can easily exponentiate the matrix analytically. Efficient quantum simulation algorithms also exist for diagonal Hamiltonian matrices with a smooth structure (Welch et al. Reference Welch, Greenbaum, Mostame and Aspuru-Guzik2014).

Because $H$ is block diagonal in the action basis, we can perform Hamiltonian simulations separately in each $(s_2, s_3)$ subspace. We span the wavefunction by

(2.8)\begin{equation} |\phi(t)\rangle = \sum_{s_2, s_3=0}^{\infty} \sqrt{p^{s_2, s_3}} \sum_{j=j_{\min}}^{s_2} c^{s_2, s_3}_j(t) |\phi^{s_2, s_3}_j\rangle. \end{equation}

The time-independent $p^{s_2, s_3}\ge 0$ is the probability that the state is within the $(s_2, s_3)$ subspace, and the total probability $\sum _{s_2, s_3}p^{s_2, s_3} = 1$. Because different subspaces decouple, the wavefunction is normalized separately within each subspace $\sum _{j} |c^{s_2, s_3}_j|^2=1$. After performing Hamiltonian simulations for the subspaces, which can be done in parallel, one can sum their contributions using weights that are predetermined by initial conditions.

For plasma problems, the initial states are often tensor-product states $|\phi (t=0)\rangle =|\psi _1\rangle \otimes |\psi _2\rangle \otimes |\psi _3\rangle$, where the three waves are initially unentangled. For each wave, $|\psi _i\rangle =\sum _m \beta ^i_m|m\rangle$ can be expanded in its Fock basis, where $p^i_m=|\beta ^i_m|^2$ is the probability that the $i$th wave occupies its $m$th state. For example, for classical laser pulse compression, the pump and seed lasers are initially in a coherent state $|\xi \rangle _c$, where $\xi =r \,{\rm e}^{\mathrm {i}\theta }$ is a complex number. For the coherent state, the probability amplitudes are $\beta _m=\exp (-\frac {1}{2}r^2) \xi ^m/\sqrt {m!}$. In this case, $p_m$ is a Poisson distribution, with $\langle n \rangle = r^2$ and $\langle n^2\rangle = \langle n \rangle ^2 + \langle n \rangle$, and $p_m$ peaks near $m\sim r^2$. When $m\gg 1$, $p_m\simeq \exp (-r^2) (er^2/m)^m/\sqrt {2{\rm \pi} m}$ decays faster than exponential. As another example, for interactions between quantum light, the initial state may be a squeezed vacuum state $|\xi \rangle _q$, for which $\beta _{2m+1}=0$ for odd states and $\beta _{2m}=(\textrm {sech}\, r)^{1/2} (-\frac {1}{2} \,{\rm e}^{\mathrm {i}\theta } \textrm {tanh}\, r)^m \sqrt {(2m)!}/m!$ for even states. In this case, $\langle n \rangle = \textrm {sinh}^2 r$ and $\langle n^2\rangle = 3\langle n \rangle ^2 + 2\langle n \rangle$. The probability $p_{2m}$ monotonically decreases, and $p_{2m}\simeq \textrm {sech}\, r (\textrm {tanh}\, r)^{2m}/\sqrt {{\rm \pi} m}$ decays much slower than a coherent state but faster than an exponential. In both examples, due to the superexponential decays of $p_m$, it is sufficient to keep track of a finite number of states in the Fock basis. Expanding the initial tensor-product state in the action basis gives

(2.9)\begin{equation} \sqrt{p^{s_2, s_3}}c^{s_2, s_3}_j(t=0) = \beta^1_{s_2-j} \beta^2_j \beta^3_{s_3-s_2+j}. \end{equation}

The simplest case is when only the pump wave is initially excited, which means that the seed and idler waves are initially in the ground state $\beta ^2_m=\beta ^3_m=\delta _{m,0}$. In this case, $p^{s_2, s_3}=|\beta ^1_{s_2}|^2\delta _{s_2, s_3}$ and $c^{s_2, s_3}_j(t=0) = \delta _{j,0}$. At later time, the wavefunction is spanned by $|\phi (t)\rangle = \sum _{s=0}^\infty \beta ^1_s\sum _{j=0}^s c^{s,s}_j(t)|\phi ^{s,s}_j\rangle$. Another special case is when the idler is initially in the ground state, namely $\beta ^3_m=\delta _{m,0}$. In this case, $p^{s_2, s_3}=|\beta ^1_{s_3}|^2 |\beta ^2_{s_2-s_3}|^2\varTheta _{s_2,s_3}$, where the step function $\varTheta _{i,j}=0$ when $i< j$ and $\varTheta _{i,j}=1$ when $i\ge j$. The initial condition is $c^{s_2, s_3}_j(t=0)=\delta _{j, s_2-s_3}$, and at a later time the wavefunction is spanned by $|\phi (t)\rangle = \sum _{s_2=0}^\infty \sum _{s_3=0}^{s_2} \beta ^1_{s_3} \beta ^2_{s_2-s_3} \sum _{j=s_2-s_3}^{s_2} c^{s_2,s_3}_j(t)|\phi ^{s_2,s_3}_j\rangle$. In both examples, the wavefunction initially occupies the ground state, namely the $j=j_{\min }$ state, within each $(s_2, s_3)$ subspace.

For plasma simulations, the observables of interest are often the wave amplitudes, namely the expectation values $\langle n_i \rangle$ for the three waves. Due to action conservation, the three amplitudes are not independent. For given initial conditions, the conserved actions are $\langle S_2 \rangle = \sum _{s_2, s_3 = 0}^\infty p^{s_2, s_3} s_2$ and $\langle S_3 \rangle = \sum _{s_2, s_3 = 0}^\infty p^{s_2, s_3} s_3$. Because the probabilities $p^{s_2, s_3}$ are time-independent, the expectation values $\langle S_2 \rangle$ and $\langle S_3 \rangle$ are constants of motion. The amplitudes of the pump and idler waves are related to the amplitude of the seed wave by

(2.10a)\begin{gather} \langle n_1(t) \rangle = \langle S_2 \rangle - \langle n_2(t)\rangle , \end{gather}

(2.10b)\begin{gather}\langle n_3(t) \rangle = \langle S_3 \rangle - \langle S_2 \rangle + \langle n_2(t) \rangle. \end{gather}

To compute the seed amplitude, we use the orthonormality of the action basis, which gives $\langle n_2 \rangle = \langle \phi |n_2|\phi \rangle = \langle \phi | \sum _{s_2, s_3, j} \sqrt {p^{s_2, s_3}} c^{s_2, s_3}_j j |\phi ^{s_2, s_3}_j\rangle = \sum _{s_2, s_3, j} p^{s_2, s_3} |c^{s_2, s_3}_j|^2 j$. The amplitude can be computed in two steps:

(2.11a)\begin{gather} \langle n_2^{s_2, s_3} (t) \rangle = \sum_{j=j_{\min}}^{s_2} j |c^{s_2, s_3}_j(t)|^2, \end{gather}

(2.11b)\begin{gather}\langle n_2(t) \rangle = \sum_{s_2, s_3 = 0}^\infty p^{s_2, s_3} \langle n_2^{s_2, s_3} (t) \rangle, \end{gather}

where the first step computes the expectation value in each $(s_2, s_3)$ subspace and the second step performs weighted sums using predetermined probabilities. For interactions involving quantum light, one may also be interested in observing higher-order cumulants, which can be measured in similar ways.

2.3. Exact dynamics

In the remaining part of this paper, we focus on a single subspace and drop the $s_2$ and $s_3$ superscripts to keep notations compact. When $s_2=s_3=0$, the subspace is trivial because $D=1$ is one-dimensional. For non-trivial dynamics, either $s_2$ or $s_3$ must be positive. Notice that although the underlying three- and four-wave interactions couple three and four photons, the dimension of the subspace $D=\min (s_2, s_3)+1$ can take any integer value. The smallest non-trivial problem is $D=2$, which requires a single qubit.

The exact dynamics involves two fundamental frequency scales $g$ and $R$, from $H_T$ and $H_F$, respectively (equations (2.4) and (2.5)). When $g=0$, the dynamics is trivial because $H_F$ is diagonal: under the influence of $H_F$ alone, the occupation of $|\phi _j\rangle$ remains unchanged, and the dynamics is a pure phase precession. To change occupation numbers, a non-zero $g$ is needed. Hence, for non-trivial dynamics, we can always normalize time by $\tau =|g|t$ and normalize the four-wave coupling by $\rho =R/|g|$. The Schrödinger equation becomes $\textrm {i}\partial _\tau \boldsymbol {c} = \boldsymbol{\mathsf{H}}\boldsymbol {c}$, where the vector $\boldsymbol {c}$ is the expansion coefficients such that $|\phi \rangle =\sum _{l=0}^{D-1} c_l|\phi _{k}\rangle$. The index $k=j_{\min }+l$, so $c_0$ is the probability amplitude that the seed wave is in the lowest occupied state and $c_{D-1}$ is the probability amplitude that the seed wave is in the highest occupied state. The normalized matrix elements $H_{k',k}=({1}/{|g|})\langle \phi _{k'}|H|\phi _k\rangle$ are non-zero only along the diagonal and first off-diagonal bands:

(2.12a)\begin{gather} H_{k,k} ={-}\tfrac{1}{2}\rho k(k-1), \end{gather}

(2.12b)\begin{gather}H_{k-1,k} = {\rm e}^{\textrm{i}\theta}\sqrt{k(s_2+1-k)(s_3-s_2+k)}, \end{gather}

where ${\rm e}^{\textrm {i}\theta } = \textrm {i} g/|g|$ and $H_{k,k-1}=H^*_{k-1,k}$. Notice that for the amplitude $c_l$, the index $k=j_{\min }+l$ may be shifted. For example, when $s_2=3$ and $s_3=4$, we have $D=4$ and $j_{\min }=0$, so $k=l=0,1,2,3$. On the other hand, when $s_2=4$ and $s_3=3$, even though we still have $D=4$, now $j_{\min }=1$, so $k=l+1=1,2,3,4$. In terms of $l$, we can rewrite $H_{k-1,k} ={\rm e}^{\textrm {i}\theta }[l(D-l)(|s_3-s_2|+l)]^{1/2}$, which is symmetric under $s_2\leftrightarrow s_3$. The three-wave interaction is a phase-sensitive process because when $\varphi \ne \varphi '$, $|{\rm e}^{\textrm {i}\varphi } + {\rm e}^{\textrm {i}\varphi '}|<2$. Here, the phase $\varphi$ contains a contribution from $\theta$, as well as a dynamical phase, which accumulates in time due to the diagonal four-wave interaction. If the relative phase between two adjacent levels changes, then the population transfer between them is reduced.

Since the Hamiltonian $\boldsymbol{\mathsf{H}}$ is time-independent, the exact dynamics is described by the unitary evolution operator $\boldsymbol{\mathsf{U}}(\tau )=\exp (-\textrm {i}\boldsymbol{\mathsf{H}}\tau)$, which can be obtained by direct diagonalization and exponentiation, at least for small problem sizes. For size $D=4$, the exact behaviours of three examples are shown in figure 1. In all three examples, the coupling phase $\theta =0$ and the constants of motion are $s_2=4$ and $s_3=3$, which means $j_{\min }=1$ and $k=l+1$, so $|\phi _{k}\rangle =|3-l,1+l,l\rangle$. All examples start the time evolution from $c_0=1$ and $c_l=0$ for $l=1,2,3$ at $\tau =0$. This ground state of the computational basis is by no means the ground state of the dynamical system. In fact, at $\tau =0$, there are three photons in the pump wave $a_1$ and one photon in the seed wave $a_2$, while the plasma wave $a_3$ is initially unexcited. At small $\tau$, the three-wave interaction consumes $a_1$ to produce $a_2$ and $a_3$, so the probabilities cascade from the ground state to higher states of the computational basis. At larger $\tau$, the pump is depleted, so the inverse process dominates, where $a_2$ and $a_3$ merge into $a_1$, and the probabilities cascade back to lower states. As can be seen from figure 1, the upward and downward cascades of probabilities repeat, but the dynamics is not periodic.

Figure 1. Exact dynamics of mixed three- and four-wave interaction problems in a $D=4$ dimensional Hilbert space with constants of motion $s_2=4$ and $s_3=3$. Starting from the ground state, the probability amplitudes $\boldsymbol {c}$ are evolved in time, and the occupation probabilities $P_l=|c_l|^2$ (a–c) as well as the expected quanta in the three waves $\langle n_i\rangle$ (d–f) are computed on a classical computer. When $\rho =R/|g|=0.1$, three-wave interaction dominates; when $\rho =2$, three- and four-wave interactions compete; when $\rho =10$, four-wave interaction dominates.

The dynamics is controlled by the dimensionless parameter $\rho$. When $\rho \ll 1$, as shown in figures 1(a) and 1(d), three-wave interaction dominates, which causes population transfer between the three waves. In this case, the much weaker four-wave interaction slowly accumulates phase mismatches that reduces the efficiency of population transfer, which is manifested by the decreasing oscillation amplitudes of $\langle n\rangle$ in figure 1(d). Here, the expected number of quanta in the three waves are calculated from the probability amplitudes $\boldsymbol {c}$ by $\langle n_1 \rangle = \sum _l (s_2-j)|c_l|^2$, $\langle n_2 \rangle = \sum _l j|c_l|^2$ and $\langle n_3 \rangle = \sum _l (s_3-s_2+j)|c_l|^2$, where $j=j_{\text {min}}+l$, and the summation is over $l=0,\ldots,D-1$. From the above equations, it is clear that $\langle S_2\rangle =s_2$ and $\langle S_3\rangle =s_3$ are exact constants of motion, as marked by horizontal dashed lines in figure 1(d–f). In the opposite limit $\rho \gg 1$, as shown in figures 1(c) and 1(f), four-wave interaction dominates. In this case, the phases of different $|\phi _j\rangle$ precess at drastically different rates, which inhibits population transfer because three-wave interaction requires phase matching. In this example ($\rho =10$), the $0\leftrightarrow 1$ transfer is strongly suppressed, while the $1\leftrightarrow 2$ transfer, which accumulates phase mismatch at a greater rate, becomes nearly impossible. Finally, in the intermediate case $\rho \sim 1$, as shown in figures 1(b) and 1(e), three- and four-wave interactions compete, and the dynamics is more complicated. While four-wave interactions generate phase mismatches that suppress population transfer, three-wave interactions change the populations and affect how the phases are weighted. The intermediate cases are where simulations are most needed for predicting the behaviour of the system.

The expectation value $\langle n (\tau ) \rangle$ can be measured efficiently using $O(\log D)$ circuits. Due to action conservation, the three expectation values $\langle n_1 \rangle = s_2 - \langle n_2 \rangle$, $\langle n_2 \rangle = j_{\min } + \langle \boldsymbol{\mathsf{O}} \rangle$ and $\langle n_3 \rangle = s_3 - s_2 + \langle n_2 \rangle$ are derivable from

(2.13)\begin{equation} \langle \boldsymbol{\mathsf{O}} \rangle = \sum_{l=0}^{2^n-1} l |c_l|^2 = \frac{2^n-1}{2} - \sum_{j=1}^n 2^{n-1-j} \langle \phi|\boldsymbol{\mathsf{Z}}_j|\phi\rangle, \end{equation}

where we embed the $D$-level system into the Hilbert space of $n=\lceil \log _2(D) \rceil$ qubits. We identify $|l\rangle = |\phi _{j_{\min }+l}\rangle$ and map this state to $|q_1\cdots q_n\rangle$ when $l=(q_1\cdots q_n)_2$, where $q=0$ or $1$ and $(q_1\cdots q_n)_2 = 2^{n-1}q_1 + \cdots 2^0 q_n$ is the binary representation of $l$. When $N=2^n>D$, we pad the state vector with zeros, namely we set $c_l=0$ for $l\ge D$. On the right-hand side of (2.13), $|\phi \rangle = \sum _l c_l |l\rangle$ is the quantum state, and $\boldsymbol{\mathsf{Z}}_j$ is the Pauli $\sigma _z$ gate acting on the $j$th qubit. By measuring $n$ single-qubit Pauli strings $\langle \boldsymbol{\mathsf{Z}}_j\rangle$ for $j=1,\ldots, n$, (2.13) allows us to obtain $\langle \boldsymbol{\mathsf{O}} \rangle$ using $n=O(\log D)$ measurement circuits with $n$ steps of post-processing. To see why (2.13) holds, we start from $\boldsymbol{\mathsf{O}}=\sum _l l |l\rangle \langle l |$. In the computational basis, $\boldsymbol{\mathsf{O}}=\textrm {diag}(0,1,\ldots, N-1)$ is a diagonal matrix, which can be expanded in Pauli basis as $\boldsymbol{\mathsf{O}}=\sum _l r_l \sigma ^n_l$, where $\sigma ^n_l$ is a tensor product of $n$ Pauli matrices. Explicitly, $\sigma ^n_{(q_1\cdots q_n)_2}:=\bigotimes _{i=1}^n \sigma _{3q_i}$, where $\sigma _{3q}=\boldsymbol{\mathsf{I}}$ when $q=0$ and $\sigma _{3q}=\boldsymbol{\mathsf{Z}}$ when $q=1$. Using the fact that $\text {tr}(\sigma ^n_i\sigma ^n_j)=2^n\delta _{ij}$, where $\delta _{ij}$ is the Kronecker delta, the expansion coefficient equals $r_l=2^{-n}\text {tr}(\sigma ^n_l \boldsymbol{\mathsf{O}})$. To calculate the trace, we need diagonal elements of $\sigma ^n_l$, and we denote its $k$th diagonal element by $[\sigma ^n_l]_k$. By induction, one can show that $[\sigma ^n_{(q_1\cdots q_n)_2}]_{(b_1\cdots b_n)_2} = \prod _{i=1}^n (-1)^{q_ib_i}$. Then, $\text {tr}(\sigma ^n_{(q_1\cdots q_n)_2} \boldsymbol{\mathsf{O}}) = \sum _k [\sigma ^n_{(q_1\cdots q_n)_2}]_k [\boldsymbol{\mathsf{O}}]_k = \sum _{b_1,\ldots, b_n=0,1} [\sigma ^n_{(q_1\cdots q_n)_2}]_{(b_1\cdots b_n)_2} (b_1\cdots b_n)_2$ $= \sum _{b_1,\ldots, b_n=0,1} \prod _{i=1}^n (-1)^{q_ib_i} \sum _{j=1}^n 2^{n-j} b_j \,{=}\, \sum _{j=1}^n 2^{n-j} \xi _j$. To calculate the sum $\xi _j=\sum _{b_1,\ldots, b_n=0,1} \prod _{i=1}^n (-1)^{q_ib_i} b_j$, one can show by induction that if $q_i=0\ \forall i\ne j$, $\xi _j= 2^{n-1}(-1)^{q_j}$, otherwise $\xi _j=0$. Consequently, only $n+1$ traces are non-zero. The first non-zero trace is when $q_i=0$ for all $i$, in which case $\sigma ^n_{(0\cdots 0)_2} = \boldsymbol{\mathsf{I}}$ and $\text {tr}(\boldsymbol{\mathsf{O}})=2^{n-1}(2^n-1)$. Second, when $q_i=\delta _{ij}$ for a given $j$, $\sigma ^n_{(0\cdots 1_j \cdots 0)_2} = \boldsymbol{\mathsf{Z}}_j$ is a single-qubit gate and $\text {tr}(\boldsymbol{\mathsf{Z}}_j\boldsymbol{\mathsf{O}})=-2^{n-1} 2^{n-j}$. When more than one $q$ is non-zero, the trace is zero. Therefore, we obtain $\boldsymbol{\mathsf{O}}=2^{-1}(2^n-1)\boldsymbol{\mathsf{I}}-\sum _{j=1}^n 2^{n-1-j}\boldsymbol{\mathsf{Z}}_j$, from which (2.13) follows.

The exact quantum dynamics can be simulated efficiently using quantum Hamiltonian simulations (Berry et al. Reference Berry, Ahokas, Cleve and Sanders2007,  Reference Berry, Childs, Cleve, Kothari and Somma2014; Low & Chuang Reference Low and Chuang2017), exploiting the fact that our Hamiltonian matrix is 3-sparse. For example, the qubitization algorithm of Low & Chuang (Reference Low and Chuang2017) requires $O(1)$ ancilla qubits in addition to $O(\log D)$ qubits that encode the action basis. The query complexity of the qubitization algorithm, namely the number of terms one need to keep in the Jacobi–Anger expansion, is $O(\tau \|H\|_{\max } + \log (1/\epsilon )/\log \log (1/\epsilon ))$, where $\epsilon$ is the desired precision and $\|H\|_{\max }$ is the matrix element with the largest absolute value. For plasma pulse compression, one is typically interested in attaining the maximum seed wave intensity. As shown in figure 1, the first peak of $\langle n_2\rangle$ tends to be the highest. So, it is usually sufficient to simulate the dynamics for the first few oscillations. In other words, the maximum simulation time of interest is $\tau =O(1/|\lambda |_{\max })$, where $|\lambda |_{\max }$ is the eigenvalue of $H$ that has the largest absolute value. Therefore, the query complexity for the pulse compression problem is $O(\|H\|_{\max }/|\lambda |_{\max })$. In other problems, such as simulating the chaotic regime of wave–wave interactions, one might be interested in long-time dynamics, in which case the complexity may be higher. To estimate $\|H\|_{\max }$, using (2.12), the diagonal elements attain the maximum at the largest $k$ value, which is $k=s_2$. When $s_2\rightarrow \infty$, $\|H_{k,k}\|_{\max }\simeq \frac {1}{2}\rho s_2^2$. For off-diagonal elements, we express $H_{k,k-1}$ in terms of $l$. When $s_3\rightarrow \infty$ at fixed $s_2$, or $s_2\rightarrow \infty$ at fixed $s_3$, the maximum is attained at $l\simeq D/2$, where $\|H_{k,k-1}\|_{\max }\simeq (D/2)\sqrt {|s_2-s_3|+D/2}$. When $s_2=s_3\rightarrow \infty$, the maximum is attained at $l\simeq 2D/3$ where $\|H_{k,k-1}\|_{\max }\simeq 2(D/3)^{3/2}$. Although $\|H\|_{\max }$ grows with $D$, the complexity remains $O(1)$ because $|\lambda |_{\max }$ also grows with $D$ at a similar rate. To estimate the eigenvalues, consider two limits. (i) Consider the limit $\rho \rightarrow \infty$. When $s_2> s_3$, because $j_{\min }>0$, diagonal matrix elements dominate, and the eigenvalues are $\lambda _l\simeq H_{k,k}\rightarrow -\infty$, for $l=0,\ldots,D-1$. When $s_2\le s_3$, because $j_{\min }=0$, the first two diagonal elements are zero. The eigenvalues are $\lambda _0\simeq |H_{0,1}|, \lambda _1\simeq -|H_{0,1}|$ and $\lambda _l\simeq H_{k,k}$ for $l=2,\ldots,D-1$, where $|H_{0,1}|^2=(D-1)(|s_2-s_3|+1)\ll \rho$. In both cases, $|\lambda |_{\max } \simeq \|H\|_{\max }$. (ii) Consider the limit $\rho \rightarrow 0$. Due to the absence of diagonal elements, the eigenvalues are either $0$ or appear in $\pm \lambda _i$ pairs, where $i=1,\ldots, \lfloor D/2\rfloor$. From the characteristic polynomial, $\sum _{i=1}^{\lfloor D/2\rfloor } \lambda _i^2 = \sum _{l=1}^{D-1} |H_{k-1,k}|^2=\frac {1}{12}D(D^2-1)(D + 2|s_2-s_3|)$. Because $\sum _{i=1}^{\lfloor D/2\rfloor } \lambda _i^2\le \lfloor D/2\rfloor |\lambda |^2_{\max }$, we obtain a lower bound $|\lambda |_{\max }\gtrsim [\frac {1}{6}(D^2-1)(D+2|s_2-s_3|)]^{1/2}$. Comparing with expressions for $\|H_{k,k-1}\|_{\max }$, we see $\|H\|_{\max } = O(|\lambda |_{\max })$ when $D\rightarrow \infty$. For both limits, as well as for intermediate $\rho$, the complexity of quantum Hamiltonian simulation is $O(\|H\|_{\max }/|\lambda |_{\max })=O(1)$, so the pulse compression problem can be simulated efficiently.

3.1. Implementation on a superconducting device

As the first test of the quantum device, we use it to enact the exact unitary operator. In figure 2, the solid blue lines are the exact occupation probabilities of the four states in our computational basis, which are computed using exact exponentiation on classical computers. The test problem uses parameters $\rho =2$, $\theta =0$, $s_2=4$ and $s_3=3$, which are identical to those of figure 1(b,e). The exact solutions serve as references for results on the quantum device. This first test is the simplest task that a quantum hardware can perform: for each time $\tau =N\varDelta$, we compute the unitary exactly on a classical computer. The sequence of unitary operations for this experiment is thus constituted of just a single unitary, $\boldsymbol{\mathsf{U}}(N\varDelta )$, and the results are shown in figure 2 as the black dashed lines. As can be seen from the figure, even when enacting a single dense unitary on the device, the fidelity is far from perfect. The important point is that because each time step uses the same gate depth, the performance does not degrade with $\tau$, except perhaps at the very beginning where most population remains in the ground state and the unitary is near identity. In this test, because the gate sequence is so short, decoherence is not a leading cause of infidelity. Instead, most infidelity comes from coherent gate errors, in the sense that each gate realizes a slightly different unitary than what is intended.

Figure 2. Occupation probabilities in a test problem with parameters $\rho =2$, $\theta =0$, $s_2=4$ and $s_3=3$. The blue lines are exact solutions from a classical computer, and the coloured symbols with error bars are results from the quantum device. When asked to enact the final unitary (black), the device performance is acceptable but not ideal. However, when asked to perform time evolution (orange), results on the device degrade to noise level after a few oscillations. The results are significantly improved using error-mitigation techniques (red), but the error bars grow exponentially.

In this simplest test, another source of error is readout, for which we have already corrected using an iterative Bayesian unfolding technique (Nachman et al. Reference Nachman, Urbanek, de Jong and Bauer2020). On the device level, the dispersive readout process is a scattering experiment, where a microwave pulse, whose frequency is off-resonant from the qubit transition frequency, is injected to interact with the qubit. The amplitude and phase of the returned microwave pulse depend on the state of the qubit, and therefore enables a measurement of the qubit states. During the readout process, the state of the qubit may change, for example, due to intrinsic decoherence or induced interactions from the readout pulse. As a consequence, the true probability $p_i$ that qubits are in state $i$ may differ from the measured probability $m_i$. The vectors $\boldsymbol {m}$ and $\boldsymbol {p}$ are related by a response matrix $\boldsymbol{\mathsf{R}}$ by $\boldsymbol {m}=\boldsymbol{\mathsf{R}}\boldsymbol {p}$, where $\boldsymbol{\mathsf{R}}$ is also called the confusion matrix, which is ideally close to the identity matrix. In practice, the matrix $\boldsymbol{\mathsf{R}}$ is constructed during calibrations, which prepare known quantum states using single-qubit gates, assuming they are ideal, and then immediately measure the qubit states. Due to statistical noise, directly inverting the matrix to find $\boldsymbol {p}=\boldsymbol{\mathsf{R}}^{-1}\boldsymbol {m}$ often leads to artefacts, including $p\notin [0,1]$ and $\sum _i p_i\ne 1$. These artefacts can be removed using iterative Bayesian unfolding $p^{n+1}_i=\sum _j (m_jR_{ji}p_i^n/\sum _kR_{jk}p_k^n)$, where $n$ is the iteration number. When the iteration converges, the denominator $\sum _kR_{jk}p_k=m_j$ cancels with the numerator. Because the marginalized probability $\sum _j R_{ji}=1$ for all $i$, $p^{n+1}_i=p^n_i$ converges to a stationary true value. In practice, starting from the initial guess $\boldsymbol {p}^0=\boldsymbol {m}$, a few iterations are usually sufficient. Without readout error mitigation, results differ only slightly from the black dashed lines in figure 2, which suggests that the leading error is not readout error but rather coherent gate errors.

As the second test, we perform time evolution using the exact unitary $\boldsymbol{\mathsf{U}}(\varDelta )$, and the results are shown by the orange dashed lines in figure 2. In this set of experiments, $\boldsymbol{\mathsf{U}}(\varDelta )$ is compiled to native gates, and the gate sequence is repeated $N$ times to enact $\boldsymbol{\mathsf{U}}^N(\varDelta )$. Because of the repetition, as $\tau =N\varDelta$ increases, the gate depth increases linearly. The accumulation of errors leads to a degradation of fidelity, as can be seen from figure 2. The oscillation amplitudes decrease and $p(\tau )$ deviates further from the true solution as $\tau$ increases. At even larger $\tau$ values, the quantum states become fully scrambled, so $p\rightarrow 1/4$ approaches the fully mixed value for the four quantum states. Because $\boldsymbol{\mathsf{U}}^N(\varDelta )$ has a larger depth than $\boldsymbol{\mathsf{U}}(N\varDelta )$, the device performs worse in this test (orange lines) than in the previous test (black lines) as expected. The $\boldsymbol{\mathsf{U}}^N(\varDelta )$ results improve noticeably from Shi et al. (Reference Shi, Castelli, Wu, Joseph, Geyko, Graziani, Libby, Parker, Rosen and Martinez2021a) primarily because of the SQISW gate, which has significantly shorter duration and higher fidelity than the two-qubit gate used in our previous work.

3.2. Improving results with error suppression and mitigation

Because product formulas require even larger gate depth, we need to improve the results for the exact unitary before moving on to the next test. The dominant source of error on the quantum processor are two-qubit gate errors, which are typically an order of magnitude larger than single-qubit gate errors. On the Aspen architecture, this is not only due to a significantly longer two-qubit gate time, but also because activating the two-qubit gate requires tuning one of the qubits away from its optimal operating point so the qubit becomes more sensitive to flux noise (Didier Reference Didier2019), which leads to a higher dephasing rate for the qubit.

We take a multi-pronged approach to minimizing two-qubit gate errors. First, the Aspen-M-3 chip used in our experiment has $\sim$200 calibrated two-qubit gates available. We require only one of these for this circuit and are thus able to select high-performing candidates based on the reported fidelities. In our experiment, we perform this selection manually, but this optimization can be performed by compilers and is often referred to as addressing. Second, the Aspen chip offers both native CZ and XY($\theta$) gates, thus providing a choice of how to express our problem unitary. We observe that the XY family is particularly expressive (Peterson et al. Reference Peterson, Crooks and Smith2020), allowing the expression of our target unitary using two XY(${\rm \pi} /2)$ gates and single-qubit gates (Huang et al. Reference Huang, Wang, Wu, Ding, Ye, Kong, Zhang, Ni, Song and Shi2023). The parametric XY($\theta$) gate is composed of two pulses with a swap angle of ${\rm \pi} /2$ while the variable angle is achieved by a phase shift between the pulses (Abrams et al. Reference Abrams, Didier, Johnson, da Silva and Ryan2020). However, we find that by using SQISW as our native gate, we can discard the second pulse from the XY gate, cutting the duration in half and reducing the two-qubit gate error by around 40 %. Our native two-qubit gate is thus a 64 ns SQISW which is combined with single-qubit rotations to produce highly expressive native cycles. Our target unitary is expressible in either two or three of these clock cycles depending on parameters of the problem. Our combined approach to native gate selection, which takes into account both the fidelity and expressiveness of the native two-qubit gates, allows us to achieve sufficiently low error rates so that it is possible to reach relatively large depths.

While our calibration and gate selection techniques suppress errors, setting the fundamental performance ceiling of the circuit, a family of error-mitigation techniques allows us to tailor the noise, and recover unbiased estimates of observables. As mentioned earlier, a major source of error is coherent gate error. Coherent errors may be the result of either control errors, which should have ideally been calibrated away, or crosstalk, which are coherent errors caused by the state of nearby idling spectator qubits. Such errors can be suppressed using dynamical decoupling (Tripathi et al. Reference Tripathi, Chen, Khezri, Yip, Levenson-Falk and Lidar2022; Zhou et al. Reference Zhou, Sitler, Oda, Schultz and Quiroz2023; Evert et al. Reference Evert, Izquierdo, Sud, Hu, Grabbe, Rieffel, Reagor and Wang2024). We do not attempt this here because our circuit does not have significant idle periods. Finally, we may suffer from unintended control errors from neighbouring control signals. All three classes of errors are major contributors to infidelity on current devices. Moreover, while long-term changes, such as temperature drifts that affect the control electronics, can be removed by routine calibrations, current devices also suffer from short-term and unpredictable changes. For example, chemical residues from fabrication process or material defects can interact with the qubits, which changes qubit frequencies and coherence times (Klimov et al. Reference Klimov, Kelly, Chen, Neeley, Megrant, Burkett, Barends, Arya, Chiaro and Chen2018; Müller, Cole & Lisenfeld Reference Müller, Cole and Lisenfeld2019; Cho et al. Reference Cho, Beck, Castelli, Wendt, Evert, Reagor and DuBois2023). As another example, when the device is struck by cosmic rays, whose energy is often higher than that of the superconducting gap, Cooper pairs are broken, creating quaisparticles, and thus changing the qubit frequency (Vepsäläinen et al. Reference Vepsäläinen, Karamlou, Orrell, Dogra, Loer, Vasconcelos, Kim, Melville, Niedzielski and Yoder2020). Due to these changes, a control pulse that is calibrated to enact a specific unitary for the original set of qubit parameters will now deviate from the intended unitary operation. The resultant coherent gate errors are particularly damaging because their behaviours are difficult to predict. Rather than adding up systematic errors in a simple way, quantum interference may cause the errors to transiently disappear, only to re-emerge at a later time. In comparison, errors due to decoherence, or more specifically due to depolarizing noise, behave in a much simpler and more predictable way: the depolarizing errors simply lead to an exponential decay of the Bloch vector, which can be easily modelled once the decay exponent is measured.

To mitigate coherent errors, we convert them into stochastic Pauli errors using a random compilation technique (Wallman & Emerson Reference Wallman and Emerson2016; Hashim et al. Reference Hashim, Naik, Morvan, Ville, Mitchell, Kreikebaum, Davis, Smith, Iancu and O'Brien2020). The technique exploits the fact that the decomposition of a target unitary into elementary gates is not unique. Suppose one has found a particular decomposition $\boldsymbol{\mathsf{U}}=(\prod _{k=1}^K \boldsymbol{\mathsf{C}}_{k}\boldsymbol{\mathsf{G}}_k)\boldsymbol{\mathsf{C}}_{0}$, where $\boldsymbol{\mathsf{C}}$ are single-qubit gates and $\boldsymbol{\mathsf{G}}$ are two-qubit gates. Then, a different but equivalent decomposition can be constructed in two steps. First, the single-qubit gates are replaced by $\tilde {\boldsymbol{\mathsf{C}}}_{k}=\boldsymbol{\mathsf{T}}_{k+1}\boldsymbol{\mathsf{C}}_{k} \boldsymbol{\mathsf{T}}^c_{k}$ for $k=0,\ldots,K$, where $\boldsymbol{\mathsf{T}}_{k+1}$ is chosen at random from a subgroup of $\boldsymbol{\mathsf{C}}$ called the twirling group. An important property of the twirling group is that $\boldsymbol{\mathsf{T}}\boldsymbol{\mathsf{G}}=\boldsymbol{\mathsf{G}}\boldsymbol{\mathsf{T}}'$. In other words, a single-qubit gate in $\boldsymbol{\mathsf{T}}$ can be commuted across two-qubit gates to become another single-qubit gate. Due to this property, once $\boldsymbol{\mathsf{T}}_k$ is chosen, if one takes $\boldsymbol{\mathsf{T}}^c_k=\boldsymbol{\mathsf{G}}_k\boldsymbol{\mathsf{T}}^{\dagger} _k\boldsymbol{\mathsf{G}}^{\dagger} _k$ for $k=1,\ldots,K$, then $\boldsymbol{\mathsf{T}}^c_k\boldsymbol{\mathsf{G}}_k\boldsymbol{\mathsf{T}}_k=\boldsymbol{\mathsf{G}}_k$. Taking $\boldsymbol{\mathsf{T}}^c_0=\boldsymbol{\mathsf{I}}$ to be the identity, one thus finds another decomposition $\boldsymbol{\mathsf{U}}=\boldsymbol{\mathsf{T}}^{\dagger} _{K+1}(\prod _{k=1}^K \tilde {\boldsymbol{\mathsf{C}}}_{k}\boldsymbol{\mathsf{G}}_k)\boldsymbol{\mathsf{C}}_{0}$. Second, to reduce unnecessary single-qubit gate depth, the new single-qubit operations $\tilde {\boldsymbol{\mathsf{C}}}_{k}$, for $k=1,\ldots,K-1$, and $\boldsymbol{\mathsf{T}}^{\dagger} _{K+1}\tilde {\boldsymbol{\mathsf{C}}}_{K}$ are compressed and simplified to elementary single-qubit gates. In this paper, the two-qubit gate we use is the SQISW gate $\boldsymbol{\mathsf{G}}=\exp [-{\rm i}({{\rm \pi} }/{8})(\boldsymbol{\mathsf{X}}\otimes \boldsymbol{\mathsf{X}} +\boldsymbol{\mathsf{Y}}\otimes \boldsymbol{\mathsf{Y}})]$, where $\boldsymbol{\mathsf{X}}$ and $\boldsymbol{\mathsf{Y}}$ are Pauli matrices. Notice that the twirling does not directly touch the two-qubit gates, which are hard to implement on hardware. The randomness of the compilation is introduced solely from the intermediate single-qubit gates, which are easy to adjust. Wallman & Emerson (Reference Wallman and Emerson2016) shows that if the Pauli twirling gates are chosen independently and if hardware Pauli errors are also independent, then random compilation transforms any gate errors into stochastic Pauli errors. In other words, as shown in the Appendix, suppose the errors of a quantum channel, when represented by the Pauli-transfer matrix, have off-diagonal components before twirling. Then, after twirling, the errors become purely diagonal, which means coherent interference of errors is removed. Because our native two-qubit gate is a non-Clifford gate, we cannot apply full Pauli twirling. Rather, we use a pseudo-twirling technique which tailors a smaller subset of coherent errors using the group of single-qubit rotations which can be successfully inverted by single-qubit gates. This twirling group is less powerful than Pauli twirling, but still tailors the noise effectively in most situations. The twirling is performed using the TrueQ software library (Beale et al. Reference Beale, Boone, Carignan-Dugas, Chytros, Dahlen, Dawkins, Emerson, Ferracin, Frey and Hincks2020). In our experiments, we construct the logical circuit and compute 50 random compilations. Each compilation has an identical pulse schedule, and thus an identical noise model. The randomization of single-qubit gates is performed by updating angles of our virtual $\boldsymbol{\mathsf{Z}}$ gates. Such updates can be made with high efficiency, allowing a large number of randomizations of the circuit to be executed in quick succession.

The final step of error mitigation is to compensate for the suppression of observables using a rescaling technique (Ville et al. Reference Ville, Morvan, Hashim, Naik, Lu, Mitchell, Kreikebaum, O'Brien, Wallman and Hincks2022). After twirling of a unitary operation $\boldsymbol{\mathsf{U}}$, the noise channel becomes approximately $\mathcal {E}(\rho )=\sum \lambda _P \boldsymbol{\mathsf{P}}\rho \boldsymbol{\mathsf{P}}^{\dagger}$, where the summation is over all tensor products of single-qubit Pauli operators $\boldsymbol{\mathsf{P}}$. The coefficient $\lambda _P$, called the Pauli decay constant, is specific to the Pauli operator $\boldsymbol{\mathsf{P}}$ but is independent of the unitary $\boldsymbol{\mathsf{U}}$ that is being performed. Because the $\lambda _P$ values are bounded by their mean $\bar {\lambda }$ as $2\bar {\lambda }-1\le \lambda _P\le 1$, Ville et al. (Reference Ville, Morvan, Hashim, Naik, Lu, Mitchell, Kreikebaum, O'Brien, Wallman and Hincks2022) propose to use a single $\bar {\lambda }$ value to correct for all Pauli errors. This approximation becomes exact when the errors are fully depolarizing, which means that all Pauli channels decay in the same way. In this case, the measured expectation value $\tilde {E}$ for any $E$ of interest is given by a simple rescaling $\tilde {E}=\bar {\lambda }E$, because the length of the Bloch vector, which measures the purity of the state, shrinks by $\bar {\lambda }$. For example, in our two-qubit problem, the expected occupation of a state beyond its fully mixed value is $(\tilde {p}-\frac {1}{4})=\bar {\lambda }(p-\frac {1}{4})$, where $p$ is the occupation probability for a pure state after one unitary operation. Then, after $N$ unitary operations, we can purify the probability by a rescaling $p=\frac {1}{4}+\frac{1}{\bar {\lambda }^{N}}(\tilde {p}-\frac {1}{4})$. In other words, after measuring the probability of a state, we subtract the noise $\frac {1}{4}$ and amplify the remaining signal exponentially by a rescaling factor $(1/\bar {\lambda })^{N}$. We provide justifications for using the simple rescaling technique in the Appendix. Notice that while amplifying the signal, this purification procedure also amplifies statistical error bars exponentially.

In practice, we estimate the mean value of Pauli decay constants $\bar {\lambda }$ using cycle error reconstruction (Erhard et al. Reference Erhard, Wallman, Postler, Meth, Stricker, Martinez, Schindler, Monz, Emerson and Blatt2019; Carignan-Dugas et al. Reference Carignan-Dugas, Dahlen, Hincks, Ospadov, Beale, Ferracin, Skanes-Norman, Emerson and Wallman2023). Cycle error reconstruction is a technique for measuring the Pauli infidelities of a dressed cycle with multiplicative precision. A Pauli-twirled cycle results in a noise channel which is diagonal in the Pauli basis. After the qubits are prepared in an initial state, the cycle of interest is repeated $M$ times and inverted to create an identity. By repeating the procedure at different values of $M$, we can we extract a state-preparation and measurement (SPAM) robust error rate for the initial state. The ensemble of error rates allow us to reconstruct the Pauli error rates of the cycle and their orbital averages, as shown in the Appendix. Together, this set of measurements provides not only an overall measurement of the fidelity, but structured information about the error channel which can be used in more advanced error-mitigation techniques. The technique we use is scalable with the number of qubits, provided that the noise is sparse, which is a reasonable assumption in superconducting qubit architectures. Empirically, the result converges with the square root of the sampling size. In our case, the gate of interest is the SQISW gate. For a system of two qubits, it is affordable to perform cycle error reconstruction over all possible Pauli channels. The SQISW is a non-Clifford gate, meaning that Pauli twirls cannot be inverted by a following pair of single-qubit rotations. However, for the purposes of characterization, this is not a problem. The inversions are propagated to the end of the benchmarking circuit, and then inverted with a final SQISW gate. This procedure of measuring the purification constant $\bar {\lambda }$ is performed each time we run a batch of experiments on the hardware.

After performing the above error-mitigation steps, the hardware results for our test problem are shown in figure 2 by the dotted red lines. The mitigated results of the second test now closely track the exact solutions, and the test performs even better than the first test (black lines), which does not use any mitigation. While the mitigation significantly improves the signals, without noticeably increasing the hardware overhead, the price we pay is exponentially growing error bars. At even larger simulation depth, the error bars will become comparable to the signals, beyond which the simulations need to stop.

4.1. Product formulas

For given problem parameters, we compile unitary matrices $\boldsymbol{\mathsf{U}}_T(\tau )=\exp (-\textrm {i}\boldsymbol{\mathsf{H}}_T\tau )$ and $\boldsymbol{\mathsf{U}}_F(\tau )=\exp (-\textrm {i}\boldsymbol{\mathsf{H}}_F\tau )$ to native gates using a Cartan decomposition (Smith et al. Reference Smith, Peterson, Skilbeck and Davis2020). Then, we use $\boldsymbol{\mathsf{U}}_T$ and $\boldsymbol{\mathsf{U}}_F$ to approximate the exact $\boldsymbol{\mathsf{U}}$. To first order, we approximate a single step with time step size $\tau =\varDelta$ using the formula

(4.1)\begin{equation} \boldsymbol{\mathsf{U}}_1(\varDelta) = \boldsymbol{\mathsf{U}}_T(\varDelta) \boldsymbol{\mathsf{U}}_F(\rho\varDelta), \end{equation}

where the normalized four-wave coupling $\rho$ appears as a scaling of time, which is convenient on quantum computing platforms that support parametric compiling. The error per step of the first order formula is $O(\varDelta ^2)$, whose prefactor is proportional to the norm of the commutator $[\boldsymbol{\mathsf{H}}_T, \boldsymbol{\mathsf{H}}_F]$. To second order, we use the symmetric formula

(4.2)\begin{equation} \boldsymbol{\mathsf{U}}_2(\varDelta) = \boldsymbol{\mathsf{U}}_T \left(\frac{\varDelta}{2}\right) \boldsymbol{\mathsf{U}}_F(\rho\varDelta) \boldsymbol{\mathsf{U}}_T\left(\frac{\varDelta}{2}\right), \end{equation}

whose error per step is $O(\varDelta ^3)$. One could instead place $\boldsymbol{\mathsf{U}}_F$ on the outside. With either placement, when repeating $\boldsymbol{\mathsf{U}}_2$, the adjacent $\boldsymbol{\mathsf{U}}$ of the same kind can be merged to reduce the required number of operations. For example, using (4.2), we can simplify $\boldsymbol{\mathsf{U}}^2_2(\varDelta )=\boldsymbol{\mathsf{U}}_T(\varDelta /2) \boldsymbol{\mathsf{U}}_F(\rho \varDelta )\boldsymbol{\mathsf{U}}_T(\varDelta ) \boldsymbol{\mathsf{U}}_F(\rho \varDelta )\boldsymbol{\mathsf{U}}_T(\varDelta /2)$. Implementing this unitary sequence thus requires compiling three different unitary matrices $\boldsymbol{\mathsf{U}}_T(\varDelta /2)$, $\boldsymbol{\mathsf{U}}_T(\varDelta )$ and $\boldsymbol{\mathsf{U}}_F(\rho \varDelta )$. In a similar spirit, higher-order product formulas require more types of unitary operations. For example, to third order, we use the product formula

(4.3)\begin{equation} \boldsymbol{\mathsf{U}}_3(\varDelta) = \boldsymbol{\mathsf{U}}_T\left(\frac{7\varDelta}{24}\right) \boldsymbol{\mathsf{U}}_F\left(\frac{2\rho\varDelta}{3}\right) \boldsymbol{\mathsf{U}}_T\left(\frac{3\varDelta}{4}\right) \boldsymbol{\mathsf{U}}_F\left(-\frac{2\rho\varDelta}{3}\right) \boldsymbol{\mathsf{U}}_T\left(-\frac{\varDelta}{24}\right) \boldsymbol{\mathsf{U}}_F(\rho\varDelta), \end{equation}

whose error per step is $O(\varDelta ^4)$. Notice that two of the six unitary operations above are evolving backward in time. This feature is common for high-order product formulas: in order to cancel higher-order errors, more operations are needed and the time step sizes become increasingly constrained such that negative values become necessary. To achieve even higher-order approximations, we use Suzuki's symmetric recurrence formula (Suzuki Reference Suzuki1990), which allows construction of $\boldsymbol{\mathsf{U}}_{k+2}$ from $\boldsymbol{\mathsf{U}}_k$ by

(4.4)\begin{equation} \boldsymbol{\mathsf{U}}_{k+2}(\varDelta) = \boldsymbol{\mathsf{U}}^2_k(p\varDelta) \boldsymbol{\mathsf{U}}_k[(1-4p)\varDelta] \boldsymbol{\mathsf{U}}^2_k(p\varDelta), \end{equation}

where $p=1/(4-4^{1/(k+1)})>1/4$, so the middle step is an evolution backward in time. Suzuki's formula has a self-similar structure once it is fully expanded in terms of elementary $\boldsymbol{\mathsf{U}}_T$ and $\boldsymbol{\mathsf{U}}_F$ operations. As the order $k$ increases, $\boldsymbol{\mathsf{U}}_k$ becomes increasingly accurate with $O(\varDelta ^{k+1})$ errors per step, but the required number of elementary operations increases exponentially as $O(5^{k/2-1})$, which incurs a significant computational cost. Thus, in practice, high-order Suzuki formulas are rarely used and low-order formulas offer the best compromise of accuracy versus speed.

With error-mitigation techniques, we are able to run experiments on the quantum device for up to about 200 two-qubit gates. The gate depth is deep enough that we can begin to compare results of different product-formula algorithms, which are shown in figure 3. In this set of experiments, different product formulas are employed to reach the same final simulation time $\tau _f=3$ at a fixed gate depth. We choose not to exhaust the maximum gate depth such that error bars at the final time remain small.

Figure 3. At a fixed gate depth, different product formulas are used for a test problem with parameters $\rho =4$, $\theta =0$, $s_2=3$ and $s_3=3$. For clarity, results from four cases are each split in two separate panels. The exact results are obtained by exponentiating the total Hamiltonian on a classical computer. Results of product formulas on a classical computer are shown by open symbols, and results obtained on quantum devices are shown by filled symbols with error bars. Because higher-order algorithms require more gates per step, for a fixed gate depth, they must use larger time step sizes to reach the same targeted final time $\tau _f=3$, leading to larger discretization errors.

With a fixed gate budget, because lower-order algorithms (equations (4.1) and (4.2)) require fewer unitary operations per step, they can afford to use smaller time step sizes. In figure 3, both first-order (blue) and second-order (orange) algorithms yield results that match closely with the exact solutions (dashed lines), which are obtained on a classical computer by exact diagonalization and exponentiation of the total Hamiltonian matrix. Moreover, results on the quantum device (filled symbols) closely match results when the same product formula is used on a classical computer (open symbols).

In contrast, higher-order algorithms require significantly more unitary operations per step, and thus can only afford to use a much larger $\varDelta$ for a fixed total gate depth. At third order (equation (4.3)), the time resolution is still sufficient to capture the oscillatory dynamics, and results on hardware are close to expected results (figure 3, green symbols). However, at fourth order (equation (4.4)), the coarse time step leads to a large discretization error. Although product-formula results (figure 3, red symbols) on a quantum device remain close to results on a classical computer, the product formula no longer provides a good approximation to the exact dynamics, as can be seen from the deviations of the open red diamonds from the black dashed lines. The approximation becomes worse at even higher orders (not shown), which require exponentially more gates per step.

It is worth emphasizing that product formulas indeed become more accurate at higher orders, provided that the time step size $\varDelta$ is fixed. In our tests, higher-order algorithms perform worse because $\varDelta$ is changed, such that the total gate depth does not exceed what is viable on the quantum device. An analogy here is the run time on classical computers. While higher order algorithms are more accurate at a fixed resolution, they require more operations and therefore longer run time. When given a fixed run time, one is forced to use a coarser resolution, in which case higher order algorithms may perform worse than lower order algorithms.

4.2. Optimal use of limited quantum resources

On current quantum devices, which do not yet have operational error correction, the maximum gate depth is limited. To make the best use of the limited quantum resources, we can adjust the choice of algorithms and resolutions for a given problem. In our case, the pulse compression problem seeks to determine the final seed laser intensity at the end of the interactions. The final time is set, for example, by the duration of laser pulses or the time to traverse the size of the mediating plasma. As a test problem, we fix $\tau _f=1$ with parameter values $\rho =4$, $\theta =0$, $s_2=3$ and $s_3=3$. We perform Hamiltonian simulation on the quantum device using product formulas to evolve quantum states, and the measured occupation probabilities are post-processed to compute the expectation values of the three waves using (2.13). We measure the error of the simulation by $\epsilon =\{({1}/{N}) \sum _{k=1}^N [n(k\varDelta )-\langle n(k\varDelta )\rangle ]^2\}^{1/2}$, where $n$ is the exact result on a classical computer and $\langle n\rangle$ is the expectation value obtained from the quantum hardware. In other words, we define the overall error of the simulation to be the 2-norm between the exact and measured time series, normalized by the number of time steps. The average error per step, $\epsilon$, receives higher contributions from later steps of the simulation. In the definition of $\epsilon$, it makes no difference whether we use $n_1$, $n_2$ or $n_3$, because $s_2=n_1+n_2$ and $s_3=n_1+n_3$ are exact constants of motion. Notice that the error of $\langle n\rangle$ is different from, albeit correlated with, the errors in the unitary $\boldsymbol{\mathsf{U}}$. The unitary error, which is also known as the process infidelity, gives a more complete characterization of the hardware performance. But the expectation-value error is of more interest to the pulse compression problem: it is the same type of error that would typically be determined using a classical algorithm, and is much easier to measure than full process tomography in experiments.

The overall error receives contributions from two fundamental sources. First, algorithmic errors often arise when simulating a dynamical system with knowledge of only the exact solutions of its non-commuting subsystems. Algorithmic errors are unavoidable even on classical computers. In our test problem, we assume that the separate three- and four-wave unitary can be implemented exactly, and then use product formulas to approximate the total unitary. In this case, any finite time step size introduces a discretization error $\epsilon _\varDelta$, which can be reduced either by using higher-order formulas at fixed $\varDelta$ or by using the formula at a fixed order but with decreasing $\varDelta$. When using the Suzuki formula, demanding errors to scale as $\varDelta ^{q+1}$ per step requires $M_q$ operations, which grows exponentially with $q$. On the other hand, to reach a target final time $\tau _f$, the number of steps $N=\tau _f/\varDelta$ increases only linearly when decreasing $\varDelta$. The total algorithmic error $\epsilon _1=O(\tau _f\varDelta ^q)$ can in principle be made arbitrarily small by increasing $q$ and decreasing $\varDelta$. However, in practice, given a limited run time, finite algorithm precision must be chosen. The trade-off between using a larger $q$ versus a smaller $\varDelta$ is strongly influenced by the second fundamental source of error: the hardware error. We measure hardware error by $\epsilon _Q$, the error per unitary operation, which is analogous to round-off errors on classical computers. On future error-corrected quantum computers, it will be possible to suppress $\epsilon _Q$ to arbitrarily small values. However, on current noisy devices, with only error mitigation rather than error correction, $\epsilon _Q$ is substantial. Using randomized compilation, we transform coherent errors into random Pauli errors, which contribute to depolarizing noise together with intrinsic quantum decoherence. After error mitigation, the error for one operation becomes independent from that for the previous operation, so the total hardware error $\epsilon _2=O(NM_q\epsilon _Q)$ accumulates linearly with the number of operations in the worst-case scenario. Because the algorithmic error $\epsilon _1$ is independent of the hardware error $\epsilon _2$, the overall error is $\epsilon =\epsilon _1+\epsilon _2$. Notice that $\epsilon _1$ decreases with $N$, whereas $\epsilon _2$ increases with $N$, so there is an optimal resolution $\varDelta$ at which $\epsilon$ is minimized.

The trade-off between hardware and algorithmic errors is demonstrated by a suite of experiments, whose results are shown in figure 4. We test the four product formulas using a common test problem, whose parameters are $\rho =4$, $\theta =0$, $s_2=3$, $s_3=3$ and $\tau _f=1$. For each order of the product formula, the overall error first decreases with $N$ due to the reduction of algorithmic errors at finer resolution. However, when $N$ exceeds an optimal value $N_*$, the error starts to increase due to the accumulation of hardware errors. At small $N$, the resolution is too coarse to resolve the dynamics, so the errors are $O(1)$, which are comparable to the signals. If $N_*$ had been larger, one would expect to see that $\epsilon$ for higher-order algorithms is smaller and decreases at a steeper slope. In our tests, because $N_*$ is not large enough, such a behaviour is not clearly observed. At large $N$, where the accumulation of hardware errors dominates, $\epsilon$ increases roughly linearly with $N$ for all orders. Because higher-order algorithms use more operations per step $M_q$, higher-order curves reside above lower-order ones in the log–log plots of $\epsilon$, except between orders $1$ and $2$. Notice that although $M_1=2$ and $M_2=3$, after merging adjacent unitary operations of the same type, as discussed after (4.2), the first-order sequence has $2N$ operations, while the second-order sequence has $2N+1$ operations, which is only slightly larger. The advantage of the second-order formula is not apparent because $\epsilon$ measures the error in $n$ rather than in $\boldsymbol{\mathsf{U}}$. On a classical computer, we observe a similar behaviour: at small $N$, while the second-order formula has a smaller error in $\boldsymbol{\mathsf{U}}$, it has a larger error in $n$; at large $N$, while the second-order formula has a larger error in $\boldsymbol{\mathsf{U}}$, it has a smaller error in $n$. These behaviours are perhaps peculiar to our mixed three- and four-wave interaction problem. At intermediate $N$, higher-order algorithms tend to have a smaller optimal $N_*$. This behaviour can be understood from $\epsilon _1\simeq C/N^q$ and $\epsilon _2\simeq AB^qN$, where $B^q$ comes from the exponential scaling of $M_q$. The minimum of $\epsilon (N)$ is reached at $N_*=({1}/{B})(R q)^{1/(q+1)}$, where $R=BC/A$. The function $N_*(q)$ is not monotonic: it increases to a maximum before decreasing as $q^{1/q}$. When the hardware error is small compared with the algorithmic error, $R$ is large, in which case the maximum of $N_*(q)$ is reached at a small $q$ value. This means that $N_*$ decreases with order, which is what we observe in our experiments.

Figure 4. For a given problem, minimum error is obtained at the trade-off between algorithmic errors and hardware errors. All test problems use common parameters $\rho =4$, $\theta =0$, $s_2=3$ and $s_3=3$. The targeted final time $\tau _f=N\varDelta =1$ is fixed, so a finer resolution $\varDelta$ requires more steps $N$, which means algorithmic errors decrease with $N$ at the expense of accumulating more hardware errors. An optimal resolution exists, where the overall error is minimized. At higher order, the optimal $N$ shifts towards lower resolution, and the minimum error does not improve with the algorithm order.

While hardware errors are negligible in classical computers, current quantum devices have significant hardware errors. These errors make the trade-off with algorithmic errors a relevant issue. It is worth pointing out that when considering the effect of round-off errors, such a trade-off also exists for classical computers, but, due to the high precision available on today's classical computers, this is far less consequential. In the small hardware error limit, $R=BC/A\to \infty$, so $N_*$ monotonically decreases for $q\ge 1$. At a fixed $q$, $N_*$ increases with $R$, which means a smaller hardware error can support the use of a finer resolution. On a typical classical computer, $A\sim 10^{-12}$ for our test problem, so that $N_*\sim 10^4$, at which the total error $\epsilon \sim 10^{-8}$ is negligible. In contrast, current quantum hardware has error $A\sim 10^{-2}$, which means in our test problem, the hardware can only support $N_*\sim 10^1$, for which the total error $\epsilon \sim 10^{-2}$ is noticeable. In such a scenario, there is no clear benefit of using higher-order product formulas because at a given hardware error, such formulas favour the use of lower resolution. But when the resolution becomes too low to resolve the dynamics, the overall error becomes worse with increasing algorithm order. In our test, higher-order algorithms do not perform better than a first-order algorithm at their respective optimal $N_*$. However, this may change on future quantum devices where $\epsilon _Q$ becomes even smaller. Similar conclusions have recently been reached for simulating the transverse-field Ising model and the XY model on noisy quantum computers (Avtandilyan & Pogosov Reference Avtandilyan and Pogosov2024).