Riemannian quantum circuit optimization for Hamiltonian simulation

Hamiltonian simulation, i.e. simulating the real time evolution of a target quantum system, is a natural application of quantum computing. Trotter-Suzuki splitting methods can generate corresponding quantum circuits; however, a faithful approximation can lead to relatively deep circuits. Here we start from the insight that for translation invariant systems, the gates in such circuit topologies can be further optimized on classical computers to decrease the circuit depth and/or increase the accuracy. We employ tensor network techniques and devise a method based on the Riemannian trust-region algorithm on the unitary matrix manifold for this purpose. For the Ising and Heisenberg models on a one-dimensional lattice, we achieve orders of magnitude accuracy improvements compared to fourth-order splitting methods. The optimized circuits could also be of practical use for the time-evolving block decimation algorithm.


I. INTRODUCTION
Hamiltonian simulation is a natural and promising application of quantum computing [1,2].For example, quantum time evolution gives access to the dynamical behavior of strongly correlated quantum systems, can quickly generate entanglement, and is a central ingredient of several quantum algorithms, like the HHL algorithm and quantum phase estimation.
A detailed numerical analysis of Trotter-Suzuki splitting methods [3] shows that they can approximate the time evolution with circuit depth scaling essentially linearly in simulated time.Recently, the authors of [4][5][6] have proposed and implemented the idea of optimizing variational circuit Ansätze inspired by Trotterized time evolution for the purpose of Hamiltonian simulation (using parametrized circuit gates).Here, we build upon the same idea and adapt a tensor network perspective as in [6], but take a step further by regarding the circuit gates as general unitary matrices, analogous to [7,8].Our main technical innovation is a derivation of how to employ the Riemannian trust-region algorithm [9, chapter 7] for the purpose of optimizing the quantum circuit.
Methods for approximating quantum time evolution on a quantum computer have also been explored in [10][11][12][13][14].The authors of [13,14] focus on the best-approximation of the time-dependent quantum state, instead of the overall time evolution operator considered here.Approaches based on quantum signal processing [10,11] might provide a complexity theoretic advantage in certain situations, but have the disadvantage of requiring additional auxiliary qubits and a block encoding of the Hamiltonian, which can incur a large overhead in practice.

II. NOTATION AND SETUP
Consider the unitary time evolution operator (in units of ℏ = 1) of a quantum system governed by a (time-independent) quantum Hamiltonian H defined on a lattice.We denote the local dimension of each lattice site by d, i.e., the local Hilbert space is C d .In our numerical simulations we will set d = 2, but the method works for general d.Translation invariance of the Hamiltonian is assumed throughout.Our goal is to approximate U by a quantum circuit.We designate the overall unitary transformation effected by the circuit as W (G 1 , . . ., G n ), with G ℓ ∈ U(d 2 ), ℓ = 1, . . ., n, to-be optimized quantum gates forming the circuit, as illustrated in Fig. 1 for a one-dimensional lattice and n = 3.Here and in the following, U(m) denotes the set of unitary m × m matrices.The circuit has a brick wall layout, and due to translation invariance, the gates within a layer are by construction all the same.We follow the mathematical convention of matrix chain ordering from right to left, i.e., the gates which are applied first are in the rightmost layer.
The Ansatz is motivated by the well-studied Trotterized time evolution approximation [3], which assumes that the Hamiltonian is a sum of "simpler" terms, H = Γ γ=1 H γ , such that each e −iHγ t can be exactly realized as quantum circuit; a basic example is the even-odd splitting of a Hamiltonian with nearest-neighbor interactions on a one-dimensional lattice, H = H even + H odd , and the Strang splitting approximation e −iHt = e −iHevent/2 e −iH odd t e −iHevent/2 +O(t 3 ). ( A benchmark comparison of our optimized circuits with such splitting methods is presented in Sect.VI.We will devise a numerical method for optimizing the gates G ℓ in sections IV and V.Each gate G ℓ is regarded as general unitary matrix (instead of being parametrized), which opens up the broader range of Riemannian optimization techniques.
The time parameter t is set to a fixed (model dependent) numerical value of order 1.The optimized gates can then be used on a quantum computer to reach times which are integer multiples of t by concatenating copies of the circuit.One expects that the approximation error increases only linearly with the final time [3].

III. LIGHT-CONE CONSIDERATIONS AND GENERALIZATION TO LARGER SYSTEMS
For practical reasons, we will perform the circuit optimization for rather small system sizes; nevertheless, it turns out that the circuit is a faithful representation of the time evolution operator for larger systems as well (assuming translation invariance), as already noted in [4,15], cf. the detailed mathematical analysis in [12].To provide an intuitive argument why this works, we consider the light-cone picture shown in Fig. 2. The causal correlations spread with a finite velocity, and cannot exceed the Lieb-Robinson bounds [16,17].
One observes that the following conditions have to be satisfied to arrive at a faithful approximation of the exact time evolution operator: (i) The evolution time t has to be small enough such that the extent of the light cone at t is smaller or equal to the system size.This property ensures that the light cone does not interfere with itself given the periodic boundary conditions.
(ii) The causal range of influence of the circuit gates (red thick lines in Fig. 2) has to enclose the physi- cal light cone, to be able to represent the physical information spreading.
Assuming these prerequisites hold, it becomes possible to use the circuit also for larger systems, simply by extending it with copies of the same gate G ℓ in layer ℓ.We will test this idea in our numerical experiments in Sect.VI.A practical use-case scenario consists of performing the optimization on a classical computer and small system size, and then using the optimized gates for larger systems on a quantum computer.

IV. MATHEMATICAL FORMALISM OF OPTIMIZATION ON THE MANIFOLD OF UNITARY MATRICES
In this section we briefly review the mathematical formalism for optimization under unitary constraints [7,9,18], following the notation in [9].This formalism forms the foundation for the numerical method proposed in the next section.
For fixed integer m, the set of unitary m × m matrices, with I m the m×m identity matrix, forms a mathematical manifold.U(m) can be interpreted as Riemannian submanifold of C m×m (with metric described below), and is a special case of the complex Stiefel manifold consisting of isometries.For notational simplicity, we omit the parameter m from U(m) in the following.A central concept is the tangent space of a manifold.The tangent space at a given V ∈ U is parametrized by the set of complex anti-Hermitian matrices [9]: where A † denotes the adjoint (conjugate transpose) of A. By construction, V † X is anti-Hermitian for any X ∈ T V U.
We introduce the following Riemannian metric on T V U, as in [7]: and thus the trace is real-valued since V † X and V † Y are anti-Hermitian.
Viewing U as embedded into C m×m , the corresponding projection onto the tangent space at V ∈ U reads [7,9]: with skew(A) = 1 2 (A − A † ) the anti-Hermitian part of a matrix.
We define the gradient of a smooth function f : C → R at point z = x + iy with x, y ∈ R as composed of the derivatives with respect to the real and imaginary components of z: This definition is straightforwardly generalized to functions depending on several complex numbers, e.g., f : , by applying the definition entrywise.
We will encounter the situation that the to-be optimized target function f : U → R is the restriction of a function f defined on C m×m .In this case, the gradient vector of f results from projecting the gradient vector of f onto the tangent space: For the purpose of computing second derivatives and Hessian matrices, we need some additional concepts.Let X(U) denote the set of smooth vector fields on U, following the notation of [9].We will use the unique Riemannian (Levi-Civita) connection ∇, which is formally defined as a map which is symmetric and compatible with the Riemannian metric.Intuitively, ∇ η is the derivative of a vector field in direction η.
As before, we interpret U as Riemannian submanifold of C m×m .Let ξ ∈ X(U) be a vector field.Then the derivative of ξ in gradient direction X ∈ T V U (V ∈ U) is given by [9, Eq. (5.15)] where Dξ(V )[X] is the gradient of ξ in direction X at point V .
A retraction on U [9, chapter 4] is a mapping from the tangent bundle of the unitary matrix manifold into the manifold.We have found it convenient to use the polar decomposition (V ∈ U) as retraction: where q polar (A) denotes the unitary matrix Q ∈ U from the polar decomposition of A ∈ C m×m as A = QP , with P a Hermitian positive semi-definite matrix of the same size as A. As a remark, alternative retraction methods have also been studied in the literature, based on QRdecompositions, projections, the Cayley transform and generally geodesic-like schemes [8,[19][20][21].We have found the polar decomposition to work well in practice, and leave a thorough exploration of alternative methods for future work.
In our numerical calculations, we will optimize the quantum gates G 1 , . . ., G n simultaneously (instead of one after another).Matching this procedure with the mathematical formalism requires a generalization to target functions depending on several unitary matrices: Formally, f is a function from the product manifold U ×n to the real numbers.The corresponding tangent space is the direct sum of the individual tangent spaces, cf.[7].In practice, the overall gradient vector is thus a concatenation of the individual gradient vectors in Eq. ( 8), and the overall retraction results from applying the retraction in Eq. ( 11) to the individual isometries and tangent vectors.
For minimizing f in ( 12), we will use the Riemannian trust-region algorithm [9, chapter 7].The central idea consists of a quadratic approximation of the target function in the neighborhood of a point G ∈ U ×n : for X ∈ T G U ×n , with the Riemannian Hessian The specific details for computing the gradient and Hessian will be discussed in the next section; with that, we have collected all ingredients for implementing the Riemannian trust-region algorithm [9, Algorithm 10], using the truncated conjugate-gradient method for the trustregion subproblem, see [9,Algorithm 11] and [22].

V. NUMERICAL METHOD FOR BRICK WALL CIRCUIT OPTIMIZATION
As just mentioned, we will use the Riemannian trustregion algorithm [9, chapter 7] for the optimization.Here we describe the details for the specialization to the brick wall circuit Ansatz.
We quantify the approximation error by the Frobenius norm distance ∥W − U ∥ 2 F between U = e −iHt and the brick wall circuit W . (To shorten notation, we omit the explicit t dependence of U .)We minimize this distance with respect to G = (G 1 , . . ., G n ), where G ℓ ∈ U(m) (m = d 2 ) for all ℓ = 1, . . ., n: Since both U and W (G) are unitary by construction, an expansion of the distance leads to where I denotes the identity matrix.Thus we may equivalently minimize the following target function: A graphical tensor diagram representation of Tr[U † W (G)] is shown in Fig. 3.
The gradient of f can be obtained as in Eq. ( 8).The straightforward extension of f reads f : (18) i.e., inserting general matrices into the brick wall diagram.In the present setting, since f depends on several unitary matrices G 1 , . . ., G n , the corresponding individual projections P G ℓ (ℓ = 1, . . ., n) have to be applied.
We make use of the Wirtinger formalism, summarized in appendix A, to obtain the gradient of f .The Wirtinger derivative of −Re Tr[U † W ] with respect to an entry of W is Next, applying the chain rule (A2), with W regarded as function of G ℓ , leads to Here we have used that ∂ G * ℓ W = 0, such that the second term of the chain rule vanishes.The derivative of W with respect to G ℓ can be expressed as graphical diagram (shown for three layers) as The uncontracted legs of the dashed "holes" in the network form the gradient tensor.The summation is due to the product rule.The tensor network diagram of ] then results from combining Fig. 3 with Eq. ( 21).Note that Tr[U † ∂ G ℓ W (G)] has the same dimensions as G ℓ (instead of a scalar quantity, as the trace might suggest).Finally, we can use relation (A4) to obtain the gradient of f with respect to the unitary matrices (G 1 , . . ., G n ): The gradient of f is then, according to Eq. ( 8), In practice, we have found it convenient to work with a real-valued representation of the gradient.For that purpose, we first parametrize the tangent spaces T G ℓ U, ℓ = 1, . . ., n: let denote the set of anti-Hermitian m×m matrices (which is a vector space over the real numbers).Then the following map defines an isometry between this space and the realvalued m × m matrices:  4) and ( 6), we may thus isometrically map the gradient of f to a list of real-valued matrices: For convenience, we finally reshape this gradient into a real vector of length nm 2 in our calculations.
For calculating the Hessian appearing in Eq. ( 14), note that X ∈ T G U ×n consists of a list of tangent vectors: Combined with the formula (23), we have to evaluate for all ℓ, ℓ ′ = 1, . . ., n.We use Eq.(10) for that purpose.In case ℓ ̸ = ℓ ′ , the derivative in direction X ℓ takes a similar form as in Eq. ( 21), but with the "holes" in layer ℓ filled by X ℓ .In case ℓ = ℓ ′ , there is a contribution from the trace in Eq. (27), formed by replacing one of the remaining G ℓ matrices by X ℓ and summing over all occurrences.Another contribution stems from the projector P G ℓ : for this, we first rewrite the definition in Eq. ( 6) as Thus the gradient of P G ℓ Z in direction X ℓ (Z regarded constant) equals Analogous to the gradient in Eq. ( 26), we have found it convenient to use the map s to parametrize the tangent vectors X ℓ in terms of real m × m matrices, and accordingly represent the Hessian as real symmetric nm 2 ×nm 2 matrix.Note that the Hessian matrix is (in general) not positive semidefinite.
As mentioned above, based on the gradient in Eq. ( 26) and the real-valued Hessian matrix, we now employ the Riemannian trust-region algorithm [9, Algorithm 10] combined with the truncated conjugate-gradient method for the trust-region subproblem [9, Algorithm 11], [22] to minimize the target function (17) with respect to the unitary matrices G = (G 1 , . . ., G n ) ∈ U(m) ×n .The hyperparameters of the algorithm are chosen as: initial radius ∆ 0 = 0.01, maximum radius ∆ = 0.1, and ρ ′ = 1  8 .The optimization is sensitive to the initial brick wall unitaries used as starting point, and is not always converging to the global optimum according to our numerical experiments.We have found the following two strategies useful for obtaining expedient starting values: (i) A "bootstrapping" approach, using the optimized gates from a circuit with fewer layers (typically two less) and padding identity layers (i.e., containing identity matrices) on the left and right.All the gates are then optimized simultaneously.(ii) Employing a splitting method from the literature as starting point.Given a Hamiltonian with two-body interactions, a splitting method provides twoqubit gates which form the same brick wall topology as our optimization Ansatz.By using these gates as starting point, the optimized circuit performs at least as good as the splitting method.For comparison, we demonstrate the effect of a more simplistic starting point, namely all identity matrices, for the Heisenberg model (see below).
The Python/NumPy source code of our implementation, including the numerical experiments of the following section, is available at [23].We have tested the gradient and Hessian computation by comparison with finite difference approximations of derivatives.

VI. NUMERICAL SIMULATIONS
As demonstration, we apply the numerical method in section V to several quantum lattice models.To describe the physical setup, first consider a Hamiltonian on the one-dimensional lattice Z /(L) with L sites (enumerated as 0, 1, . . ., L − 1) and periodic boundary conditions: Here the subscripts indicate that the operator ĥ acts locally on sites j, j + 1.
A natural benchmark is a comparison with Trotterized even-odd splitting schemes.On a one-dimensional lattice, the Hamiltonian is partitioned into an even and odd part, namely H = H even + H odd with By construction, the summands commute pairwise since the operators act on disjoint sites.A time step ∆t of H even/odd is exactly realized by a quantum circuit layer consisting of copies of the two-qubit quantum gate e −i ĥ∆t .A step of the Strang splitting approximation, Eq. ( 2), can then be implemented using three circuit layers with the same layout as in Fig. 1.The construction works analogously for other splitting methods.We fix a time t and quantify the approximation error by the spectral norm distance between the numerically exact time evolution operator e −iHt and the unitary matrix resulting from the splitting method or optimized circuit W (G 1 , . . ., G n ), respectively.Note that the target function (17) actually minimizes the Frobenius norm distance, since this has a canonical and simpler tensor network representation, as shown in Fig. 3.
For existing splitting schemes we consider the approximation error for an increasing number of steps r ∈ N ≥1 .Thus a single time step has size ∆t = t r and consists, e.g., of three substeps for the Strang method.We choose the convention that the "cost" of an integration method is the number n of required circuit layers.Since we can always merge the last substep with the first substep of the next time step, the number of layers for the Strang method is n = 2r + 1.In general, for a method with s substeps, we have n = (s − 1)r + 1.
Besides the Strang method, we also apply the Suzuki formula of order 4 by M. Suzuki [24], the method of order 4 by H. Yoshida [25], the Runge-Kutta-Nyström method of order 4 by R. I. McLachlan [26], and the partitioned Runge-Kutta method S 6 of order 4 by S. Blanes and P. C. Moan [27].

A. Ising model on a one-dimensional lattice
We first consider the transverse-field Ising model (TFIM) Hamiltonian on the one-dimensional lattice Z /(L) with L sites (enumerated as 0, 1, . . ., L − 1) and periodic boundary conditions: Here X j and Z j are the usual Pauli matrices acting on site j, and J, g, h ∈ R are parameters.Without loss of generality, we set J = 1.We will first consider the integrable case h = 0, and then investigate general parameter values for which the model becomes non-integrable.To demonstrate the behavior of the Riemannian trustregion algorithm, Fig. 4 visualizes the spectral norm distance and target function during the optimization iterations.For 7 circuit layers, the two curves almost overlap exactly.One also notices that a plateau is reached after around 160 iterations, such that we stop after 200 iterations.We have adopted the strategy of either using the circuit gates of an existing splitting method as starting point for the optimization, or start with the optimized gates from a circuit with less layers and pad identity gates (which are then subject to the optimization as well).
A benchmark evaluation of the numerical optimization for the integrable Ising model and t = 1 is shown in Fig. 5.For the even-odd splitting methods, we have used the two-qubit gate e −iJ(Z⊗Z+g 1 2 (X⊗I2+I2⊗X))∆t .The approximation error of the optimized circuit is represented by the thick blue curve.One observes a clear advantage: for example, the optimized circuit achieves the same or better accuracy using 9 layers as compared to the S 6 method by S. Blanes and P. C. Moan using 49 layers.The chosen parameter g = 0.75 corresponds to the ordered phase; as indication that the results are not parameter-specific, we have repeated the calculations for g = 1.5 corresponding to the disordered phase, with qualitatively same outcome (data not shown).
As discussed at the end of Sect.III, it is possible to use the optimized circuit gates also for larger system sizes L. The corresponding approximation error quantified by the spectral norm distance is shown in Fig. 5b.As expected from the light cone picture, the error hardly increases for larger L.
Naturally, one could attribute the large effectiveness of the circuit optimization to the integrable property of the TFIM.In this sense, our method could be regarded as numerical equivalent of circuit compression based on the integrable structure manifest in the Yang-Baxter equation, as studied in Refs.[28][29][30][31].Surprisingly, however, the optimization results remain qualitatively unaffected by turning on an integrability-breaking longitudinal field, i.e., setting h to a non-zero value, as shown in Fig. 5c.

B. Heisenberg model on a one-dimensional lattice
As next example, we consider the Heisenberg-type Hamiltonian with (σ 1 , σ 2 , σ 3 ) = (X, Y, Z) the vector of Pauli matrices and ⃗ J ∈ R 3 , ⃗ h ∈ R 3 parameters.In our numerical simulations we set ⃗ J = (1, 1, − 1 2 ) and ⃗ h = ( 3 4 , 0, 0).As before, the model is defined on a lattice with L sites and periodic boundary conditions.The final time is set to t = 1 4 , which is smaller than for the Ising model to compensate for the faster spreading velocity of correlations (data not shown).
The approximation error is plotted in Fig. 6a, together with a benchmark comparison of splitting methods from the literature, as before.The advantage of the optimized circuit for the Heisenberg time evolution is less  pronounced than for the Ising model, but nevertheless, the error is sill more than an order of magnitude smaller as compared to the Suzuki method.The dashed curve in Fig. 6a shows the results for a more simplistic optimization protocol, namely choosing identity gates as starting points.The gap to the best results widens with increasing number of layers, which is likely due to the higher-dimensional optimization landscape.In general, one observes that the results can sensitively depend on the starting point.As described above, employing a splitting method as starting point guarantees that the optimized circuit performs at least as well as the splitting method.
As for the Ising model, we also probe the generalizability to larger systems, by showing the error depending on L in Fig. 6b.As expected, the error increases only slightly with L.

C. Ising model on a ladder geometry
As last example, we again consider the Ising model, but now on a lattice with ladder geometry, i.e., of dimension L x × 2, with periodic boundary conditions along the ladder direction.The Ising Hamiltonian (without longitudinal field) on a general lattice reads where the first sum runs over all nearest neighbor lattice sites.For our simulation, we set J = 1, g = 3, and final time t = 1 4 .The ladder geometry requires a modification of the splitting scheme: besides the interactions in x-direction, which we can split using an even-odd scheme as before, there are now additional interactions along the steps of the ladder.Thus we split the Hamiltonian into three terms: H Ising ladder = H Ising x,even + H Ising x,odd + H Ising y .The Ansatz circuit layout for the optimization is modified accordingly, using three different gate topologies for the layers.Likewise, for the benchmark comparison we use splitting schemes supporting three terms.The Strang scheme and the methods by Suzuki and Yoshida can be adapted for this purpose.For example, the Strang splitting scheme for partitioning H = H a + H b + H c reads e −iHt = e −iHat/2 e −iH b t/2 e −iHct e −iH b t/2 e −iHat/2 + O(t 3 ).(35) In our setting, H a = H Ising x,even , H b = H Ising x,odd and H c = H Ising y .We additionally include the AY 15-6 method of order 6 by Auzinger et al. [32] in the comparison.A common feature of the splitting methods and our Ansatz is the sequential ordering of the matrix exponentials and circuit layouts, respectively, namely as "abcbabcba. . .".The approximation error after numerical optimization and a comparison with the existing splitting methods is shown in Fig. 7a.Similar to the Ising model on a onedimensional lattice, one observes a large advantage of the optimized circuits.
To probe the suitability of the optimized gates for larger systems, we use them for a 6 × 2 layout and record the approximation error in Fig. 7b.As expected, the error increases only slightly.We cannot reach even larger systems in this comparison due to the difficulty of computing and storing the exact matrix exponential.

VII. CONCLUSIONS AND OUTLOOK
In this work we have constructed and explored a numerical scheme for optimizing the two-qubit gates (G 1 , . . ., G n ) in a quantum circuit Ansatz as general unitary matrices.To use these gates in physical quantum computers still requires their further decomposition into hardware-native gates; for this purpose, one can rely on existing algorithms from the literature [33,34].For example, if the hardware supports CNOT and arbitrary single-qubit gates, each G ℓ layer can be represented by (at most) seven layers of such gates [34].In general, a good decomposition strategy and the "cost" of applying the gates will depend on the specific hardware, cf.[35].We leave a detailed use-case study and comparison with parametrized gates as in [4,6] for future work.
In the numerical experiments, we have seen that the improvement due to the optimization compared to existing splitting methods strongly depends on the model.It would be interesting to gain a deeper mathematical understanding of these differences, in particular regarding the (ir-)relevance of integrability.A related question is how to find suitable starting gates for the optimization, to arrive at a global optimum in the best case.Conversely, further improvements of the optimized circuits shown in this work might be possible.
The periodic boundary conditions are a tool to avoid finite size effects in the numerical simulations.Using the optimized gates on actual quantum computers requires additional considerations regarding the boundary conditions.We note that a one-dimensional system with periodic boundary conditions (as studied in the present work) could be realized on a physical quantum computer having a two-dimensional topology and nearest-neighbor connectivity, by mapping the logical qubits to a loop formed by physical qubits.Nevertheless, endowing the logical system with open boundary conditions would require a dedicated optimization of the brick wall circuit gates without translation invariance.We leave this task for future work.
We remark that our method can also be applied to quantum models with longer-range interactions, for which an even-odd splitting of the Hamiltonian is not feasible, assuming that the conditions based on the lightcone picture in Sect.III are still satisfied.Related to that, our optimized gates could also be of practical use for the time-evolving block decimation (TEBD) algorithm, when simulating quantum systems on a classical computer.
Another natural application is the approximation of the time evolution on a two-dimensional lattice.The reasonably smallest dimension on a square lattice is a window of size 4 × 4, to avoid interference due to the periodic boundary conditions.Performing the numerical optimization for 16 sites is technically challenging, but we plan to tackle this scenario via a tailored implementation, possibly running on a compute cluster.One could avoid very large matrices via a matrix-free application of circuit gates and matrix exponentials to statevectors, or via the "exponential-free" approach investigated in [36].
The numerical method developed in our work can approximate a general translation invariant unitary target matrix.Thus another use-case could be the decomposition of multiple-qubit gates.
In principle, the Riemannian gate optimization framework is applicable for non-translation invariant systems as well when admitting independent gates within each layer.Since the number of to-be optimized gates then increases with the system size L (or number of "orbitals" in a chemistry setting), the Riemannian trust-region part of the algorithm becomes computationally more expensive in practice.More specifically, in this scenario the overall number of gates is then n = D • L/2 for D layers.

FIG. 1 .
FIG.1.Example of a quantum circuit with brick wall layout and periodic boundary conditions, for approximating the exact time evolution operator.

FIG. 4 .
FIG. 4. Progress of the Riemannian trust-region algorithm iteration, for the Ising model Hamiltonian (32) on a onedimensional lattice with L = 6 sites, h = 0 and time t = 1.The blue curves shows the approximation error quantified by the spectral norm distance, and the orange curves the shifted and rescaled target function f in Eq. (17).

FIG. 5 .
FIG. 5. Approximation error (quantified by the spectral norm distance) of the quantum time evolution operator for t = 1 governed by the Ising model Hamiltonian (32) on a onedimensional lattice.The circuit gates have been optimized for L = 6 lattice sites.(a) and (c) show a comparison with existing splitting methods from the literature, with the thick blue curves corresponding to the optimized circuits of the present work.Subplot (b) displays the error for larger systems, always using the same optimized brickwall circuit gates (G1, . . ., Gn) found in (a).

FIG. 7 .
FIG. 7. Approximation error of the quantum time evolution operator for t =1  4 governed by the Ising model Hamiltonian (34) with J = 1 and g = 3 on a ladder geometry.The circuit gates have been optimized on a 4 × 2 lattice.(a) shows a comparison with existing splitting methods from the literature, with the blue curves corresponding to the optimized circuit of the present work.Subplot (b) displays the error for system size 6 × 2 as well, using the same optimized brickwall circuit gates (G1, . . ., Gn).