Minimum Trotterization Formulas for a Time-Dependent Hamiltonian

When a time propagator eδtA for duration δt consists of two noncommuting parts A=X+Y, Trotterization approximately decomposes the propagator into a product of exponentials of X and Y. Various Trotterization formulas have been utilized in quantum and classical computers, but much less is known for the Trotterization with the time-dependent generator A(t). Here, for A(t) given by the sum of two operators X and Y with time-dependent coefficients A(t)=x(t)X+y(t)Y, we develop a systematic approach to derive high-order Trotterization formulas with minimum possible exponentials. In particular, we obtain fourth-order and sixth-order Trotterization formulas involving seven and fifteen exponentials, respectively, which are no more than those for time-independent generators. We also construct another fourth-order formula consisting of nine exponentials having a smaller error coefficient. Finally, we numerically benchmark the fourth-order formulas in a Hamiltonian simulation for a quantum Ising chain, showing that the 9-exponential formula accompanies smaller errors per local quantum gate than the well-known Suzuki formula.


Introduction
Let us consider an initial value problem where A(t), S(t, t ′ ) ∈ C d×d and I is the identity operator in C d×d .Equation (1) frequently appears in physics.For instance, the Schrödinger equation corresponds to the case where A(t) is anti-Hermitian (i.e., A(t) = −iH(t) with a Hermitian Hamiltonian H(t)).The Lindblad equation, which is the time evolution equation for Markovian open systems, corresponds to the case where A(t) is the Liouvillian, which is not Hermitian nor anti-Hermitian in general.Therefore, solving Eq. (1) is one of the fundamental problems in physics.
The formal solution of Eq. (1) is given by the following time-ordered exponential, S(t, t ′ ) = T exp t t ′ A(s)ds . ( When the operator A(t) is time-independent, i.e., A(t) = A ∀t, Eq. (3) reduces to S(t, t ′ ) = exp((t − t ′ )A).The time-ordered exponential is defined as an infinite sum of multiple integrals and is not easy to calculate since A(t) consists of non-commuting operators in general.Thus, a decomposition that expresses (3) as a product of exponentials of simpler operators is important in applications such as a digital quantum simulation.
For the time-independent case of A = X + Y , decompositions of S(δt, 0) = e (X+Y )δt into a product of e αX and e βY have been extensively studied.The lowest-order decomposition S(δt, 0) = e δtX e δtY + O(δt 2 ) was first obtained by Lie and its applicability was extended by Trotter [24] and Kato [12].Also, various higher-order decompositions have been proposed both in classical mechanics [6] and in quantum mechanics [7].Here, Trotterization means the form of with real numbers a j 's and b j 's, where 2q or 2q+1 denotes the number of exponentials in each formula.When a j 's and b j 's are chosen so that the error of the approximation is of the order of δt p+1 , i.e., we call T (⃗ a, ⃗ b) the p-th order Trotterization formula of S(µ + δt 2 , µ − δt 2 ).Unlike other methods [14], Trotterization allows us to implement Schrödinger equation in digital quantum computers without exploiting ancillary qubits and has been widely used in noisy intermediate-scale quantum computers.
However, much less is known about Trotterization for the time-dependent case [7,8,14,19], which includes more variety of nonequilibrium quantum phenomena.For simplicity, we focus, throughout this paper, on the following type of a time-dependent generator [1] A(t) = x(t)X + y(t)Y, (6) where x(t) and y(t) are smooth real functions and X, Y ∈ C d×d .Two different approaches are known for the Trotterization formula of a timedependent operator.The first is to take x(t) and y(t) at fine-tuned discrete points t k (n = 1, • • • , N ), and the algorithm is constructed using exp(a k Xδt) and exp(b k y(t k )Y ) [7,26].This approach systematically generates higherorder formulas based on the formulas for timeindependent cases, but its efficiency, especially the number of exponentials, is not necessarily optimal.The second is an approach based on a different approximation initiated by Huyghebaert and de Readt [8,19], where the time-ordered exponential (3) is approximated by a product of normal exponentials such as exp t t ′ x(s)Xds and exp t t ′ y(s)Y ds .However, its generalization to higher orders is not straightforward and has not been systematically developed.
Here we focus on the number of exponentials in the formula as its usefulness.The number is directly related to the calculation cost for, e.g., a digital quantum computer and a matrix-productstate-(MPS-)based classical simulation.Therefore, we look for the Trotter formula of the minimum number of exponentials.For k = 2, the well-known formula so-called midpoint rule (see Eq. (12) below) is a three-exponential formula and is the minimum formula.This is because it is readily shown that there does not exist a 2exponential second-order formula in general.For k = 4, by contrast, there are some known formulas.For example, Suzuki [7,23,26] gave a 15-exponential formula for the more general case X(t) + Y (t), and this reduces to an 11exponential one for the present class of problem given in Eq. (3).However, it has remained unsolved whether there exists another fourth-order formula having fewer exponentials, and the minimum number of exponentials has not been found yet, to the best of the authors' knowledge.
In this paper, we derive minimum Trotterization formulas for time-dependent A(t) in the form of Eq. (6).In particular, we obtain the fourth-order (sixth-order) explicit formulas consisting of 7 (15) exponentials and show that 7 (15) exponentials are the minimum among all the fourth-order (sixth-order) formulas in the form of Eq. (4).Thus we call our formulas the minimum fourth-order and sixth-order Trotterization (MFT and MST) formulas.We also find a fourth-order 9-exponential formula, whose error is less than that of the MFT.We numerically compare the error between the MFT, 9exponential, and Suzuki's fourth-order formulas in 1-qubit and many-spin models.We confirm that the MFT (9-exponential) formula has a slightly larger (smaller) error per exponential than Suzuki's does.In this way, the 9-exponential formula offers an efficient way for simulating dynamics on, e.g., a digital quantum computer and an MPS-based classical simulator.

The formulas and errors
We begin with the Magnus expansion [3] As is well known, Ω n is given by a multiple integral over n variables t 1 , . . ., t n on the domain , and the integrand consists of nested commutators of A(t 1 ), . . ., A(t n ) (see Appendix A for their concrete forms).Therefore, hold in general.Also, our parameterization ensures that Ω n does not involve even-order terms in δt when Taylor-expanded.The Magnus expansion, truncated at some order, has been utilized for classical and quantum simulations (see, e.g., Refs.[10,21]).Here we use it to derive Trotterization formulas, i.e., decompositions into X and Y as in Eq. (4).
For A(t) given in Eq. ( 6), the Magnus expansion reduces to the continuous BCH formula [3] giving (11) and so on, where we introduced useful notations For brevity, we shall omit (µ, δt) from β's in the following.
We remark that Eq. (7) becomes a trivial identity in the time-independent case, where x(t) = x and y(t) = y.In this case, the exact evolution does not need the time-ordered exponential, giving S(µ Interestingly, Ω n 's are actually smaller than the obvious estimate O(δt n ).
For instance, To prove this, one can, e.g., Taylor-expand x(t j ) and y(t j ) at t j = µ, finding that the coefficient of δt 2 vanishes.Furthermore, one can show Ω 3 = O(δt 5 ) and Ω 4 = O(δt 5 ) (see Appendix A for proof).This rapid increase of orders in Ω n enables us to construct efficient Trotterization formulas as shown below.
Using the lowest-order continuous BCH formula, we can reproduce the well-known secondorder Trotterization formula S(µ also known as the mid-point rule.To do this, we notice S(µ 3 ) and use the minimum secondorder formula for a time-independent problem 3 ) and that for β 2 .
To obtain a fourth-order formula, we invoke the following key result. 7).( 14) 7 ), where we used u = O(δt 2 ) ensured by our assumption.On the other hand, the continuous BCH formula gives Comparing these, we obtain Eq. (13).
Three remarks are in order regarding Theorem 1.The first remark is on the role of the assumption β 12 /β 2 = O(δt 2 ).Although β 2 = O(δt) always holds true, β 2 can happen to be as small as β 2 = Θ(δt 3 ).Here, β = Θ(δt n ) stands for 0 < lim δt→+0 (β/δt n ) < ∞, and we used the fact that β 2 consists of odd-order terms of δt (see Appendix B).In such a case, if β 12 = Θ(δt 3 ), u = Θ(1) and Eq.(13) does not hold.However, even when we can still use Theorem 1 by interchanging the roles of x(t)X and y(t)Y .Recall that our target S(µ + δt 2 , µ − δt 2 ) is invariant under the interchange of x(t)X and y(t)Y .If both β 12 /β i (i = 1, 2) are not O(δt 2 ), Theorem 1 can be used by setting u = 0 as follows.This case implies The second remark is that Theorem 1 can also be derived using (a non-Hermitian version of) the Schrieffer-Wolff Transformation (SWT) [4,20].When we regard Ω 1 and Ω 2 as the unperturbed and perturbation terms, respectively, the SWT aims to eliminate Ω 2 within O(δt 5 ) errors by an appropriate similarity transformation (S is not necessarily Hermitian below) whose exponential form reads Assuming that S = O(δt 2 ), we obtain the condition for S satisfying Eq. ( 16) as when O(δt 5 ) errors are neglected.Recalling that we obtain a solution S = (β 12 /β 2 )X = uX for Eq. ( 18), for which Eq. ( 17) coincides with Eq. (13).We remark that the transformation is unitary if X is Hermitian.
The final remark is that Theorem 1 becomes trivial for the time-independent case, where we have u = 0 and Υ 5 .This follows from Ω n = 0 for n > 1, as remarked above.In this sense, Υ 5 is the characteristic of time-dependent cases, which we will call the time-dependent component of the Trotter error.
Having given those three remarks, we now return to the Trotterization formula's derivation.Even though we are considering a time-dependent problem, Eq. (13) tells us that the Trotterization is achieved by Trotterizing e β 1 X+β 2 Y , which is done by invoking a conventional formula for time-independent problems.For example, we can apply the Forest-Ruth-Suzuki formula [6,22] where where s = (2 − 2 1/3 ) −1 and Γ 5 = O(δt 5 ) is the error [17,18], whose properties are generally investigated [5].Substituting Eq. ( 19) into Eq.( 13), we obtain the following corollary giving the MFT.

Corollary 1.1 (Minimum 7-exponential formula). It follows
Remarkably, the 7-exponential formula (20) is minimum among fourth-order formulas.To prove this, we recall that our time-dependent problem involves the time-independent problems as a special case of x(t) = x and y(t) = y.For the special case, Eq. ( 20) reduces to the fourth-order Forest-Ruth-Suzuki formula [6,22] as one can check easily.Recall that the 7-exponential fourth-order Forest-Ruth-Suzuki formula is minimum among time-independent problems.Thus every fourthorder formula for general time-dependent problems requires at least 7 exponentials.Theorem 1 shows that such a formula actually exists despite a reasonable inference that more general problems demand more exponentials.
The error parts Γ 5 and Υ 5 of the total error Λ 5 are time-independent and -dependent components, respectively.As we discussed, Υ 5 = 0 holds for the special case that x(t) and y(t) are time-independent.Since Γ 5 is the well-known Trotter error for time-independent case [5], Υ 5 is the additional error source particular to the time dependence.
We remark on the time-reversal character of the MFT.The exact evolution satisfies then we obtain The keys to confirming Eq. ( 24) are the relations ) and β ′ 12 = −β 12 and their consequence u ′ = u.The time-reversal character (24) implies that the error is of odd order [17].
Although we have substituted the minimum Forest-Ruth-Suzuki formula (19) into Eq.( 13), we can use any fourth-order Trotterization formulas instead (see, e.g., Ref. [18] for a list of them).One can use another to reduce the error Γ 5 by using more exponentials than seven.For instance, using Omelyan's Forest-Ruth formula (see Refs. [17,18] for detail and real numbers a j 's and b j 's) in Eq. (13), we obtain the following corollary giving the 9-exponential fourth-order formula.

Corollary 1.2 (9-exponential formula). It follows
Here, in exchange for adding two more exponentials, we have smaller coefficients in Γ ′ 5 than in Γ 5 .As we will see below in an example model in Sec. 5, the 9-exponential formula produces less error than the MFT and the Suzuki formula model when compared at the same number of elementary quantum gates.
Before closing this section, we remark on the generalization beyond the fourth-order formula.
Recall that the MFT (20) reduces to the Forest-Ruth-Suzuki for time-independent cases when we set u = 0 or eliminate the SWT.This implies that the MFT can be derived by applying an appropriate SWT to the time-independent formula.This approach is promising to derive higher-order minimum Trotterization formulas and, in fact, leads to the following minimum sixth-order Trotterization (MST) formula T 15,7 consisting of 15 exponentials (see Appendix C for the derivation).
Theorem 2 (Minimum sixth-order 15-exponential formula).Suppose that β 1 = Θ(δt) and β 2 = Θ(δt).Then it follows where a 1 , . . ., a We make some remarks on Eq. (27).First, as in the MFT, Eq. ( 27) reduces to Yoshida's minimum sixth-order formula [27] in the time-independent case, where Namely, these variables are introduced as a timedependent generalization.Second, the 15 exponentials in Eq. ( 27) are placed symmetrically around the central e b 4 β 2 Y , except for the asymmetric appearance of ±u j (j = 1, . . ., 4) due to the SWTs.Third, w and z appearing symmetrically are not SWTs but are newly introduced at the sixth order, enabling T 15,7 to reproduce Ω 3 in the Magnus expansion.These structures would be generalized to even higher orders: The SWTs and the symmetric decorations like w and z would produce the minimum Trotterization formulas for time-dependent cases out of those for time-independent ones.Finally, like the MFT, the MST is a minimum sixth-order formula because it has the same number of exponentials as Yoshida's formula does for time-independent generators.

algorithm for unitary dynamics
The arguments in the previous sections have been so general that A(t) may or may not be anti-Hermitian.Nonetheless, the case where A(t) is anti-Hermitian is particularly important because it includes the Schrödinger equation, for which A(t) = −iH(t) with H(t) being a Hamiltonian.In this section, we focus on this case, summarizing the algorithm for the Hamiltonian simulation.Considering Eq. (6), we focus on the following Hamiltonian where f (t) and g(t) are real, and F and G are noncommuting Hermitian operators.To address these problems, one can utilize the above results with the replacements A(t) → −iH(t), x(t) → f (t), X → −iF , etc.As relevant applications, the MFT offers an efficient implementation in quantum dynamics simulation on a digital quantum computer and matrix-product-state(MPS-)based classical simulations.
For usage, we summarize the algorithm as a pseudocode in Fig. 1.Here we focus on the MFT since it works similarly with the other formulas.In the algorithm, we provide an ordered list of unitaries that brings an initial state |ψ(t i )⟩ at time t = t i to a final state |ψ(t f )⟩ according to the Hamiltonian (28) based on the MFT.We can also use the 9-exponential formula instead of the MFT in the algorithm by replacing the seven exponentials to append with 9 exponentials appearing in Eq. (26).
We note that each time step requires the oracle functions that evaluate integrals, which can be cumbersome.In such a case, we can do it approximately up to the fourth order of δt by using the Taylor or orthogonal-polynomial expansion for x(t) and y(t) (see also Appendix B).

Comparison to Other Time-Dependent Formulas
Here we compare the MFT and 9-exponential formulas with other Trotterization formulas for time-dependent Hamiltonians (28).The known second-order formulas are equivalent to the midpoint rule (12) within errors of O(δt 3 ).In fact, the second-order Suzuki [7,23] is the same as the mid-point rule, and another integrator T HdR =

T exp(
Huyghebaert-de Readt [8] satisfies T HdR = T 2 + O(δt 3 ) as one can check easily.Furthermore, those second-order 3-exponential formulas are minimum as they are known to be minimum for time-independent cases.Thus we take Suzuki's formula (i.e., mid-point rule) as the representative for the second-order formula.
The fourth-order formula can vary in terms of error and gate complexity depending on the choice of algorithm.The fourth-order Suzuki formula for time-dependent Hamiltonians [7,23,26] is made of the mid-point rule (12) (with x(t)X → −if (t)F and y(t)Y → −ig(t)G) as where 2 δt, and w = (4 − 4 1/3 ) −1 .Note that T 4 consists of 11 exponentials for Hamiltonians in Eq. (28) because four pairs of e −iαF (α ∈ R) can be combined into one between T 2 's.Note again that the MFT has only 7 exponentials.Now we implement each algorithm and numerically compare their accuracy.To quantify the accuracy, we define the error ϵ as where T is the time-evolution unitary achieved through an algorithm and || • || is the Frobenius norm of the operator.We have confirmed that the results do not change qualitatively when the spectral norm is used, as shown in Appendix D.
Figure 2 shows the simulation errors for the seminal Landau-Zener Hamiltonian confirming the expected error scaling (σ x and σ z are the Pauli matrices).Here we fix µ = 1, for instance, and plot the errors of each formula against the time step δt.Recall that we have two choices in assigning which of σ x and tσ z to f (t)F (and the other to g(t)G).In Fig. 2, the assignment σ x → f (t)F and tσ z → g(t)G is shown by filled symbols and the other assignment σ x → g(t)G and tσ z → f (t)F by open ones.For both assignments, the mid-point rule gives O(δt 3 ) errors, whereas the MFT, 9-exponential, and the fourthorder Suzuki formulas give O(δt 5 ) errors, as expected.
For µ = 1 in Fig. 2, we observe that the error coefficient of the MFT is larger than that of the fourth-order Suzuki formula.This is a penalty to pay with the MFT instead of its advantage of fewer exponentials.However, it is noteworthy that the 9-exponential formula has error coefficients as small as Suzuki's, even though it still has fewer exponentials.Here we recall that the errors of the MFT Λ 5 and the 9-exponential formula Λ ′ 5 consist of two parts as in Eqs.(20) and (26), including the common contribution Υ 5 inherent to the time dependence.Our observation about the difference in error coefficients derives rather from the time-independent formulas, Γ 5 and Γ ′ 5 .To investigate the µ-dependence systematically, we plot the errors in Fig. 3.Here we set δt = 0.1/∥H(µ)∥ = 0.1/ 1 + µ 2 for each µ since the magnitude of H(t) significantly changes in such a wide range.The figure shows the following two features.First, as µ → 0, the errors are divergent for the MFT and 9-exponential formulas in the assignment σ x → g(t)G and tσ z → f (t)F .This derives from the fact that, in that limit, β 1 = O(δt 2 ) and β 12 /β 1 .Hence, to obtain a better approximation, we must choose the appropriate assignment σ x → f (t)F and tσ z → g(t)G.
Second, the error of the 9-exponential formula, with the appropriate assignment, is as small as the fourth-order Suzuki's over the entire range of µ, whereas that of the MFT can be as good as them in two limiting cases µ → 0 and µ → ∞.
The worst case for the MFT is µ ∼ 1, for which the δt-scaling is shown in Fig. 2. We note again that its error comes dominantly from Γ 5 rather than the time-dependent part Υ 5 .In fact, the µ −1 -scaling and the difference between the two MFT results in the large µ regime can be understood by analyzing Γ 5 (see Appendix E for details).

Application to a quantum circuit
Having studied the basic properties of the MFT and the 9-exponential formula in a single qubit model, we here apply them to a more practical model.We consider a quantum Ising chain of length L under an oscillatory transverse field, whose Hamiltonian is in the form of Eq. (28) with  where σ i x and σ i z are the Pauli matrices acting on the site i.We impose the periodic boundary condition σ L+1 z = σ 1 z and set J = −1.0,h z = 0.2, and h x = −2.0.
Let us compare the cost and accuracy of each algorithm by considering a time evolution from t i = 0 to t f = π, for example.For a number N of steps and the step size δt = (t f − t i )/N , we introduce t k = t i + kδt (k = 0, 1, . . ., N ) to study the error where the Trotterization T (t k , t k−1 ) is taken for either the mid-point rule, MFT, 9-exponential, or Suzuki formulas.We focus on the number of 1-qubit and 2-qubit gates, N gates , in each Trotterization.For the mid-point rule, we have T (t, s) = e −i sin(µ)δtF/2 e −iδtG e −i sin(µ)δtF/2 (µ ≡ t+s 2 ), which is further decomposed into 1-and 2-qubit gates by e −i sin(µ)δtF/2 = L i=1 e −i sin(µ)δtσ i x /2 and e . Thus, we have N gates = 5LN for the mid-point rule.Likewise, we have N gates = 10LN for the MFT, 13LN for the 9-exponential formula, and 15LN for the Suzuki formula.
Figure 4 shows the error (34) of each Trotterization calculated for several N ranging from 5 to 400.Here, the horizontal axis is taken as N gates /L, and the exact solution S(t f , t i ) was obtained by integrating with sufficiently small step-size to guarantee asymptotic convergence below all other error scales.We observe expected error scalings when N gates ∝ δt −1 increases; ϵ ∝ δN −2  gates for the second-order mid-point rule and ϵ ∝ N −4  gates for the other fourth-order formulas.We note that these scaling exponents are larger by one than those for a single step because we successively operate N ∝ δt −1 ∝ N gates steps.
Importantly, the 9-exponential formula produces less error than the Suzuki formula at a given number of gates.When compared with the same N , the Suzuki formula has less error (data not shown), but it requires N gates = 15LN gates that are more than N gates = 13LN of the 9-exponential formula.This difference in the required number of gates results in the 9exponential formula obtaining the best performance, as shown in Fig. 4 In contrast, the MFT gives larger errors than the other fourth-order formula, even though it requires fewer exponentials due to the large error coefficient in Γ 5 .

Discussion and Conclusion
We studied the generalization of Trotterization formulas to time-dependent generators.While the second-order formula is well-known, we derived a 7-exponential fourth-order formula T 7,4 and a 15-exponential sixth-order formula T 15,7 as in Eqs.(20) and (27).Remarkably, we showed that the number seven (fifteen) of exponentials are the minimum to make a fourth-order (sixthorder) formula, and hence we named them the minimum fourth-order and sixth-order formulas (MFT and MST).Our approach can be generalized to higher-order formulas.An important application of Trotterization formulas is quantum dynamics simulation on digital quantum computers, where reducing the number of exponentials (i.e., quantum gates) is beneficial for suppressing errors accompanied by each gate operation.Although Trotterization-based algorithms do not have optimal error scalings of quantum-singularvalue-transformation-based ones [13,15,16,25], they have fewer overheads and are useful in nearterm quantum computers.As demonstrated in a single-qubit model, the MFT's error is as tiny as Suzuki's fourth-order formula, even though the MFT involves fewer exponentials.Also, the 9exponential formula has less error per gate than Suzuki's formula, as demonstrated in a many-spin model.These formulas would also be useful in MPS-based dynamics simulations, where fewer exponentials could speed up the simulation.In addition, the fourth-order formulas can be used for estimating the errors of lower-order Trotterizations [9].
Our result is limited in two ways.First, we assumed Eq. (6) and excluded more general time dependence A(t) = X(t) + Y (t).Equation (6) includes important applications like the adiabatic state preparation, but it would be desirable if one could generalize the assumption.Second, our result is currently limited to the decomposition into two parts X and Y .In some physically relevant models, such as the quantum Ising model with external fields and the Heisenberg model, the Trotterization for the two operators is sufficient since those Hamiltonians can be viewed as a sum of two operator sets, each of which consists of commuting operators.However, it would be useful if one could extend the MFT to more than two operators.
The error of the MFT was obtained as Λ 5 in Eq. (20) consisting of time-independent anddependent contributions.For a quantum manybody local Hamiltonian on a lattice, both contributions are proportional to the volume since they consist of (nested) commutators between local Hamiltonian parts.If we use some norm for the error operator as an error estimator, it will also be proportional to the volume.However, one is usually interested in the evolution of a given vector v, and the error is evaluated as ∥Λ 5 v∥.Such an idea was proposed to improve the error estimates by Jahnke et al. [11] for time-independent cases and generalized by An et al. [1] to the timedependent Suzuki formulas.A similar idea could be used to estimate the error of the MFT.Also, errors could be reduced by making the time step δt adaptive [28,29].We leave further investigations of the time-dependent error Υ 5 for future work.

A Magnus expansion and continuous BCH formula
The Magnus expansion [3] up to the fourth order in our notation (7) reads

C Derivation of the MST
Here we derive the minimum sixth-order Trotterization formula in Theorem 2. To begin with, we invoke the minimum sixth-order formula for time-independent operators [27], having where Note that this Yoshida formula is known to have a large error coefficient [18].To reduce the coefficient with adding more exponential, one can use the Blanes and Moan formula [2], like we derived the 9-exponential formula in the main text.Now we consider its time-dependent generalization, for which according to the Magnus expansion, where β 12 = O(δt 3 ), and β i12 and β ij12 are O(δt 5 ).We aim to decorate the right-hand side of Eq. (51) so that the left-hand side coincides with Eq. (53).
To reproduce the six commutator terms in the exponent of Eq. (53), we introduce six unknown variables to consider where we impose as working hypotheses that u 1 , . . ., u 4 are O(δt 2 ), and w and z are O(δt 3 ).Applying BCH-type formulas repeatedly and using Eq.(51), we obtain where c 12 , c 112 , . . ., c 1212 are second-order polynomials in terms of the unknowns u 1 , . . ., u 4 , w, and z.Using Eq. (52), these polynomials are given by  Thus our problem is reduced to solving the following set of six equations, for u 1 , . . ., u 4 , w, and z.Due to the linearity and independence, Eqs.(62) can be uniquely solved for u 1 , . . ., u 4 under our assumptions β 1 = Θ(δt) and β 2 = Θ(δt).Then, substituting these solutions to Eqs. (63), we have two linear and independent equations for w and z, which can also be uniquely solved.Thus we have proved that Φ with these solutions coincides with the exact propagator S(µ − δt/2, µ + δt/2) up to the sixth order of δt.

D Spectral norm
In the main text, the Frobenius norm in the error definition (30).To emphasize this norm choice, we rewrite Eq. (30) as where ∥A∥ = tr(A † A).Note that 2ϵ F represents the trace distance of S and T .We can also consider the spectral norm where ∥A∥ S coincides with the largest singular value of A. We confirm that both choices give qualitatively similar results in the benchmark calculation presented in Fig. 2. We illustrate, in Fig. 5, the error ratio obtained for the same simulation, where we observe that the different choice results only in a constant factor independent of the algorithm and time step.
E Details on Γ 5 and Γ ′

5
Here we supplement more details on the time-independent errors Γ 5 and Γ ′ 5 .According to Ref. [17], each symmetric Trotterization formula has the following relation, e (A+B)h+C 1 h+C 3 h 3 +C 5 h 5 +••• = e Aa 1 h e Bb 1 h • • • e Bbqh e Aa q+1 h , (66 where h is the time step and We have q = 3 for the Forest-Ruth-Suzuki formula and q = 4 for Omelyan's Forest-Ruth formula, for which we have different values of γ 1 , . . ., γ 6 .Omelyan's Forest-Ruth formula has more degrees of freedom that are fine-tuned to decrease an error quantifier 6 i=1 |γ i | 6 [17].The time-independent contributions Γ 5 and Γ ′ 5 are given, in the leading order of δt, by h 5 C 5 under the substitution of Ah → β 1 X and Bh → β 2 Y (or the other assignment Ah → β 2 Y and Bh → β 1 X).Now we explain why the filled and open symbols for the MFT in Fig. 3 deviate significantly in the large µ region.We recall that where we used where we used Eqs.(69) and (70) and the asymptotic behaviors of β 1 and β 2 in µ ≫ 1.Note that we obtain the observed µ −1 scaling for both cases.Interestingly, the difference between these cases comes from the imbalance of γ 4 and γ 1 , which are [17] γ 1 = −0.0004138,γ 4 = +0.0046844, and thus differ by a factor of 10.This is why we observed a magnitude difference between ■ and □ in Fig. 3.
For generic values of µ, µ-dependence of δt and more than one coefficients γ 1 , • • • , γ 6 come into play.This may cause a peak near µ ∼ 1 in Fig. 3, but we do not go into more detail since this should depend on the model and parameters.

Figure 1 :
Figure 1: Pseudocode for the Hamiltonian simulation based on the MFT.

Figure 2 :
Figure 2: Errors (30) of Trotterization formulas for the Hamiltonian (31) at µ = 1.Different symbols correspond to the fourth-order Suzuki (circle), MFT (square), 9-exponential (triangle), and mid-point rule (diamond) formulas.The filled symbols correspond to the assignment σ x → f (t)F and tσ z → g(t)G in Eqs.(28) and (31), whereas the open symbols to the other assignment σ x → g(t)G and tσ z → f (t)F in these equations.The solid line guides the eye for the power law ∝ δt 5 .

1 Figure 3 :
Figure 3: Errors (30) of Trotterization formulas for the time-dependent Hamiltonian (31) for various µ.The time step δt is taken, depending on µ, as δt = 0.1/∥H(µ)∥ = 0.1/ 1 + µ 2 .The symbols are the same as in Fig. 2: The filled symbols correspond to the assignment σ x → f (t)F and tσ z → g(t)G in Eqs.(28) and (31), whereas the open symbols to the other assignment σ x → g(t)G and tσ z → f (t)F .The solid line guides the eye for the power law ∝ µ −1 .

Figure 4 :
Figure 4: Errors (34) of Trotterization formulas for the time-dependent quantum Ising Hamiltonian [Eqs.(32) and (33)] at L = 6 plotted against N gates /L.The time evolution from t i = 0 to t f = π is calculated for steps N ranging from 5 to 400, which is converted to N gate for each Trotterization formula shown in the legend (see also text).

Figure 5 :
Figure 5: Error ratio between those evaluated by the spectral and Frobenius norm for Trotterization formulas.The parameters and legends are all shared by those in Fig. 2.