Vectorization of the density matrix and quantum simulation of the von Neumann equation of time-dependent Hamiltonians

Based oh the properties of Lie algebras, in this work we develop a general framework to linearize the von Neumann equation rendering it in a suitable form for quantum simulations. Departing from the conventional method of expanding the density matrix in the Liouville space formed by matrices unit we express the von Neumann equation in terms of Pauli strings. This provides several advantages related to the quantum tomography of the density matrix and the formulation of the unitary gates that generate the time evolution. The use of Pauli strings facilitates the quantum tomography of the density matrix whose elements are purely real. As for any other basis of Hermitian matrices, this eliminates the need to calculate the phase of the complex entries of the density matrix. This approach also enables to express the evolution operator as a sequence of commuting Hamiltonian gates of Pauli strings that can readily be synthetized using Clifford gates. Additionally, the fact that these gates commute with each other along with the unique properties of the algebra formed by Pauli strings allows to avoid the use of Trotterization hence considerably reducing the circuit depth. The algorithm is demonstrated for three Hamiltonians using the IBM noisy quantum circuit simulator.


I. INTRODUCTION
One of the most important applications of quantum computers is the simulation of large quantum systems that are intractable using conventional classical computers [1][2][3].A vast number of quantum techniques to simulate quantum dynamics has been developed with this purpose in mind [4,5].Specially many body systems give rise to extremely large Hamiltonians [6] that are classically impossible to tackle even for a small number of particles.
A very common quantum simulation problem is the digital quantum simulation (DQS) of a closed quantum system [2,4,[6][7][8][9][10][11][12][13][14] which consists in solving the timedependent Schrödinger equation.In general terms, a DQS algorithm consists of two steps.First the time evolution operator of a given Hamiltonian is expressed as a series of unitary quantum gates through the Jordan-Wigner isomorphism [11], for example.Second, the time evolution operator is used to propagate an initial state [2,5], typically by shifting locally the wave function forward in time over discrete and sufficiently small time slices [2,4,5].The Hamiltonian is written as the sum over many simpler interactions H = m i=1 H i and the evolution operator can be approximated by U(t) = e −iH1∆t/ℏ e −iH2∆t/ℏ . . .e −iHm∆t/ℏ as long as the interval ∆t is small enough.Through the Trotter-Suzuki formula [15][16][17], that takes into account the non-commutativity of the H i 's, the accuracy of U (t) can then be improved but at the cost of increasing the circuit depth.
Nearly all realistic quantum systems are open systems whose evolution is non-unitary due to the decoherence induced by the unavoidable interaction with the environment.This presents a major difficulty for digital quantum computers that can only operate with a set of unitary gates.This obstacle has spawned a wide range of solutions for the quantum simulation of open quantum systems [18][19][20][21][22][23][24].The following are some examples.The open-system nature of a two-spin NMR ensemble has been employed to simulate the decoherence in a quantum simulation [18].Also, a set of ancilla qubits that are designed to have the same effect as the simulated environment have been used to represent the decoherence of an open quantum system [20].Quantum algorithms have been developed specifically for duality quantum computers that can perform non-unitary operations [21].More recently, eigenstate and thermal state quantum simulations have been demonstrated employing the quantum analogues of imaginary time evolution and Lanczos algorithms [22].Of direct relevance to this work is a new kind of open-quantum system simulation algorithm that solves the Lindblad master equation and that has successfully been tested in real digital quantum hardware [24].The algorithm mainly relies on the Kraus matrix representations [25][26][27] and an adaptation of the imaginary time evolution algorithm [22].In general, the Lindblad equation is first vectorized in the form of a Scrhödinger equation through the Kraus operator sum [25][26][27][28] and then, the anti-Hermitian component of the Hamiltonian is treated through quantum imaginary evolution.The elements of the vectorized density matrix are finally obtained by quantum tomography.The unitary evolution of the Hermitian component of the Hamiltonian is calculated by the standard techniques of DQS listed above.
In this work we are primarily interested in the linearization of the density matrix that is governed by the von Neumann equation where H and ρ(t) are the N × N matrices for the Hamiltonian and the density matrix.The linearization is typically done by expressing the density matrix using the Kraus operator sum and column stacking the density matrix elements.As we show below, this procedure yields the following linearized von Neumann equation [24,26] where I is the N × N identity matrix and |ρ⟩ denotes the vector resulting from stacking the columns of ρ(t) from left to right.We develop a method that generalizes the vectorization process of the von Neumann equation expanding the density matrix in terms of the elements of a Lie algebra.Moreover, we show that the column-stacked density matrix in Eq. ( 2) is in fact one of the many possible ways of representing ρ(t) as a vector using the elements of a very specific Lie algebra.Instead of using this algebra, whose elements are non Hermitian, we propose expanding the von Neumann equation in terms of Hermitian matrices.The advantage is that the thus vectorized density matrix has purely real coefficients that greatly simplify the quantum tomography of |ρ⟩.We present the special case of the algebra formed by Pauli strings and with it, implement an algorithm that solves the von Neumann equation to illustrate the method.The method, together with the special properties of Pauli strings, allows to avoid the use of Trotterization therefore reducing the circuit depth.The algorithm was tested for the von Neumann equation of the magnetic resonance Hamiltonian using Qiskit in the noisy quantum circuit simulator qasm_simulator [29] .The code can be downloaded from [30].

II. VECTORIZATION OF THE DENSITY MATRIX.
Any N × N square matrix W can be expressed as the linear combination of the elements of the finite Lie algebra g n = {g 1 , g 2 , . . ., g n }, as where g ⊤ = (g 1 , g 2 , . . ., g n ), w ⊤ = (w 1 , w 2 , . . ., w n ), g i ∈ C N ×N and w i ∈ C. If W is known, the coefficients w i can be calculated as provided that the elements of the algebra are orthogonal under the Frobenius inner product, namely Not all algebras have orthogonal elements, however, one can always find a linear combination of them that does meet (5).Such combinations might be generated, for example, through the Gram-Schmidt method.One of the fundamental properties of Lie algebras is that their elements form a closed commutator algebra described by where the Lie bracket of g n has been chosen to be the commutator [g i , g j ] = g i g j − g j g i .From this definition it is straight forward to show that the structure constants c i,j,k have the property The explicit form of the structure constants is obtained by combining Eqs. ( 5) and ( 6).In the case of an algebra h n = {h 1 , h 2 , . . ., h n } formed only by Hermitian matrices the structure constants acquire a new and useful property, the invariance with respect to cyclic permutation of the indices as long as Note that the properties (8) and ( 9) imply that c i,j,k is fully antisymmetric.
With these definitions at hand we can proceed to vectorize the von Neumann equation.In terms of the most general algebra g n the density matrix can be expanded as where ρ(t) = (ρ 1 (t), ρ 2 (t), . . ., ρ n (t)) and from Eq. (4) the coefficients are Multiplying the left and right-hand side terms of (1) by Because of ( 3) and ( 4), the Hamiltonian can be expressed as the decomposition where are the Hamiltonian coefficients and a ⊤ (t) = (a 1 (t), a 2 (t), . . ., a n (t)).
Equation ( 12) can be put in the more succinct form where is the Hamiltonian superoperator, and C i is the matrix whose elements are the structure constants (C i ) j,k = c i,j,k .Equation ( 15) is a general form of the vectorized von Neumann equation.

At this point it is natural to ask what is the connection between this version of the von Neumann equation and
(2).To address this question we introduce two auxiliary algebras.First consider the algebra q n = {q 1 , q 2 , . . ., q n } formed by the N × N sparse matrices q i that only have a 1 that progressively shifts from top to bottom and from left to right as i goes from 1 to n = N 2 .More explicitly Therefore, the vector of dimension n = N 2 formed by the coefficients of a given N × N matrix in terms of q n corresponds to the stacked columns thereof.Indeed, the vector formed by the coefficients of W in terms of q n are the matrix elements of W where i 1 and i 2 are given by ( 18) and (19), respectively.The second auxiliary algebra is q ⊗2 n = q n ⊗ q n = {q 1 ⊗ q 1 , q 1 ⊗ q 2 , . . ., q n−1 ⊗ q n , q n ⊗ q n } = {Q 1,1 , Q 1,2 , . . ., Q n,n−1 , Q n,n } where Q i,j is an N 2 × N 2 sparse matrix that has a 1 in the i, j position.It is important to note that even though in general Q i,j ̸ = q i ⊗ q j , it is possible to establish a one to one relation for Following a similar procedure as in Eq. ( 12) the expansion of the von Neumann equation in the q n altakes the form In Appendix A it is shown that if A and B are two N ×N matrices then With this, Eq. ( 24) becomes Tr q † j ρ(t) , (26) that, in full consistency with Eq. ( 2), can also be written as the vectorial expression Hence, the von Neumann equation that stems from the Kraus sum is a particular case of (15) when it is expanded in terms of the algebra q n .In this case ρ(t) corresponds to the stacked columns of the density matrix.The entries of the vector ρ(t) are, however, complex numbers since q n is comprised of non-Hermitian matrices.This substantially complicates the quantum tomography of ρ(t) whose elements may have non-vanishing phases.
To circumvent this difficulty one can use any Hermitian algebra h n that yields real entries for ρ(t).In this base, the Hamiltonian and the density matrix are expressed as the linear combination of the elements of the finite Lie algebra are real and, in general, time-dependent.Similarly, the density matrix coefficients are As before, a(t) = (a 1 (t), a 2 (t), . . ., a n (t)) and ρ(t) = (ρ 1 (t), ρ 2 (t), . . ., ρ n (t)).From ( 16) the Hamiltonian superoperator is given by where Note that H is Hermitian because the coefficients a i (t) are real and the structure constants are fully antisymmetric, i.e., and therefore

III. TIME EVOLUTION
At first glance it would seem reasonable to find the evolution operator by directly applying the standard techniques of DQS to (15).However, to do so efficiently, it would be convenient to leverage the structure and the symmetries of the linear operator H in (16).These are not obvious from the expression of H.But the contraction C i a i (t) (i = 1, 2, . . ., n) signals the fact that only the n structure constants are needed to generate the time evolution [17] and not n 2 = n × n as the dimension of H would suggest.Instead of following this line of reasoning we resort to the Lie algebraic method to prove that indeed only n elements are needed to generate the time evolution.
The Lie algebraic method enables the exact determination of evolution operators for time-dependent Hamiltonians of finite dimension or having a dynamical algebra.Though the general method is thoroughly discussed in [31] a quick review is provided here for completeness.
One possible way of expressing the evolution operator of the Hamiltonian (28) is where is the unitary transformation generated by the k-th element of h n .Each U k transforms h according to where the explicit form of the matrices M k (α k (t)) is given by The matrices M k (α k (t)) are unitary since, as mentioned earlier, C k = −C ⊤ k is skew-symmetric in accordance with (34).Therefore, by performing a time evolution over h we obtain where The explicit time-dependence of the α k (t) parameters is yet unknown but can be derived from the Schrödinger equation where p t = iℏ∂/∂t is the energy operator.The Hamiltonian and the energy operator transform according to where It then follows from Eqs. ( 41) and (42) that the Schrödinger equation under U (t) becomes or, equivalently we arrive at the result that p t U (t) |ψ(t)⟩ = 0. Since p t = iℏ∂/∂t, the previous equation will only hold if ⟩ is the initial state.Equation (47) consists of a system of n ordinary differential equations for α 1 (t), α 2 (t), . . .α n (t).The  explicit time dependence of the α(t) parameters comes from the solution of (47) along with the initial condition necessary to ensure that U (0) = U(0) = 1.As we mentioned above, since both ρ(t) and h are Hermitian matrices, the elements of ρ are real.This crucial feature significantly simplify the quantum tomography of the density matrix coefficients allowing to easily determine its elements through phase kickback [32] with only one auxiliary qbit.Using Eqs.(38) and ( 29) the evolution of the density matrix operator in the Schrödinger picture is expressed according to From the comparison of this result with Eq. ( 29) it follows that the time-dependent density matrix may be cast in the shape of a vector as Except for a normalization constant, notice that M ⊤ (α(t)) is associated to the superoperator of U † (t) acting on the vectorized density matrix ρ(t) in the Fock-Liouville space generated by h n .In other words, the evolution of the density matrix thus maps onto the evolution of a quantum state vector where M (α(t)) is a super-evolution operator.The density matrix can be expanded as where ρ i (t) are the normalized density matrix coefficients ( i ρ 2 i (t) = 1) and there is a one-to-one relation between the Fock-Liouville states |h i ⟩ and the matrices h i .
Equation (50) proves that the time evolution of the vectorized density matrix is generated by only n elements that correspond to the structure constants C k as was hinted at the beginning of this section.

IV. ALGORITHM
We now provide a quantum simulation algorithm for the density matrix ρ(t) obeying the von-Neumann equation (1) expanded in a general algebra h n .It illustrates the time evolution and the quantum tomography of the density matrix.
Equation ( 50) is the underpinning element of the algorithm.It allows for the computation of ρ(t) through a series of time steps t m of step size ∆t = t m − t m−1 as where, from (47) Each time step in (54) can be regarded as an independent differential time evolution so that, from the initial condition (48), we can make α(t m−1 ) = 0 without any loss of generality.Because of this, V −1 (α(t m−1 )) = V −1 (0) = I and M ⊤ (α(t m−1 )) = M ⊤ (0) = I and consequently the first order term in Eq. (54) gives The second order term is Using ( 43), (44) and the identity the first element of (56) reduces to The second element of (56) vanishes, i.e., because c i,j,k a j (t)a k (t) is fully antisymmetric.Thus, the second order term is given by Gathering all the contributions up to sencond order in ∆t yields where dα m i is the i-th component of dα m with i = 1, 2, . . .n. Equation (47) and consequently (61) take into account the ordering of the Hamiltonian gates prescribed by Eq. ( 41) hence enabling us to avoid Trotterization provided that we aim to errors of the order of ∆t 3 .More accurate approximations of dα(t m ) could iteratively be built up to a certain tolerated error of ∆t reducing the number of time steps and quantum gates.However, this would imply a heavy load of operations in the classical stage of the computation that could undermine the effectiveness of the algorithm.Hence, it should be established at what order of ∆t the costs of calculating dα m outweigh the gain in precision.

V. QUANTUM CIRCUIT
The fundamental stages of the quantum circuit are presented in Fig. 1.The Hamiltonian and M (α(t)) are considered to have arbitrary dimension n = 2 nq−1 where n q is the number of required qbits to perform the quantum simulation considering that an extra control qbit is needed for the quantum tomography of |ρ⟩.The control qbit and the remaining ones will be termed first and second registers, respectively.
The first stage is the initial state preparation.In it, the goal is to set the initial state to the superposition The second part of the previous state initializes the density matrix coefficients to and the first is used as an ancillary uniform superposition later on used to carry out the phase kickback.Conveniently, the density matrix coefficients ρ i (0) are arranged so that the binary string of i − 1 matches the second register.
The first Hadamard gate (H-gate) puts the control qbit in the superposition From this point on the part of the state within square braces refers to the second register.The role of the controlled U H gate is to set the initial condition for the normalized components of the density matrix vector |ρ(0)⟩.This section of the circuit can be represented by the unitary transformation where P 0 = |0⟩ ⟨0| and P 1 = |1⟩ ⟨1| are the standard projectors.The purpose of this gate is to initialize the second register, whose control qbit is |1⟩, to the normalized coefficients of |ρ(0)⟩ living the other unmodified.Given that the second register is initially set to |0⟩ by default, the first column vector of the U H must be (ρ 1 (0), ρ 2 (0), . . ., ρ n (0)).There are many variants of the U H matrix that meet this requirement and that of being unitary.One option is the Householder matrix given by where u 1 = (1, 0, . . ., 0) ⊤ .The singularity produced when ρ 1 = 1 and ρ 2 = .. = ρ n = 0 can be avoided by choosing u 2 ,... or u n instead of u 1 .The following controlled H-gates set the part of the second register that is proportional to the control qbit |0⟩ to a uniform superposition.This segment of the circuit can be cast in the form of the transformation After collecting the gate transformations and from some elementary algebra it follows that the initial state is given by which is consistent with (62).
The next portion of the circuit is devoted to the actual time evolution of |ρ(t)⟩.The action of the controlled M (α(t)) gate on the initial state |ψ(0)⟩ can be expressed as As it is shown in Fig. 2, the controlled M (α) gate is broken down into small time steps according to Eq. (53).
In agreement with Eqs. ( 37) and (39) each differential time step M (dα) is decomposed into n − 1 Hamiltonian gates of the form where the Hamiltonian operator is connected to the structure constants through as can be seen in Fig. 3.These gates can be synthetized for SU (2 N ) operations with arbitrary number of qbits by means of Cartan decomposition [33][34][35][36], and other alternative methods [37,38].Yet, the complexity and depth of the quantum circuits increases exponentially with the number of qbits producing fast fidelity decays.Fortunately, the structure constants C k acquire properties that facilitate the implementation of the Hamiltonian gates if we restrict ourselves to the algebra h n whose elements are the Pauli strings The structure constants are n × n matrices that can be expanded in terms of the elements of the algebra 2 , . . ., h Thereby the Hamiltonian gate becomes where i spans only the non-vanishing coefficients Tr C k h (2) i . Furthermore, in the Appendix B we prove that in the previous expansion the elements h (2) i with non-vanishing coefficients Tr C k h (2) i commute with each other.This is quite advantageous because it allows to build the Hamiltonian gate (71) as the sequence of Pauli gates without regard to the ordering thereof, and consequently without the need of Trotterization.Additionally the Hamiltonian gates of Pauli strings can be very efficiently simulated using Clifford gates [3] .Up to this point (after the last H-gate), the wave function takes the form Finally, the density matrix coefficients ρ i (t) are computed as where p 0,i and p 1,i are the probabilities corresponding to |0⟩ ⊗ [|h i ⟩] and |1⟩ ⊗ [|h i ⟩], respectively.These are given by and are computed through the counts from the circuit execution.

VI. EXAMPLE 1: ONE SPIN 1/2 PARTICLE SUBJECT TO A TIME VARYING MAGNETIC FIELD
As a toy example, in this section we examine a spin 1/2 particle in a varying magnetic field.The time evolution of the density matrix elements ρ 1 (t), ρ 2 (t), ρ 3 (t) and ρ 4 (t) were calculated by means of the classical and quantum algorithms.The Hamiltonian for this system is In order to illustrate how to calculate expected values using the density matrix coefficients and to further compare the classical and quantum algorithms Fig. 5 exhibits the expected values of the populations of the lowest and highest energy levels and the z projection of the spin given by All the results in this section where also confirmed by directly solving the von Neumann equation ( 1) for the matrix elements of ρ(t) and then projecting on to the elements of h n using Eq.(31).

VII. EXAMPLE 2: TWO SPIN 1/2 PARTICLES COUPLED THROUGH EXCHANGE INTERACTION SUBJECT TO A OSCILLATING MAGNETIC FIELD.
In order to test a larger system, in this section we deal with the time evolution of two spin 1/2 particles coupled through exchange interaction subject to a time-varying magnetic field.The Hamiltonian for two electron spins coupled through the exchange interaction characterized by A and subject to a varying magnetic field is given by where S 1 and S 2 are the spin operators for the first and second particles, respectively.Projecting onto the base formed by the Pauli strings the non vanishing coefficients of the Hamiltonian are a 2 (t) = a 5 (t) = ω 1 cos(ωt), a 3 (t) = a 9 (t) = −ω 1 cos(ωt + ϕ), a 4 (t) = a 13 (t) = −ω 0 and a 6 (t) = a 11 (t) = a 16 (t) = A/2.The initial condition, that corresponds to both electrons being in the lowest energy level, is given by The non vanishing coefficients corresponding to this condition are ρ 1 (0) = ρ 4 (0) = −ρ 13 (0) = −ρ 16 (0) = 1/2.In this example the Hamiltonian parameters are ω = 0.9, ω 0 = 1, ω = 22.0, ϕ = π/2 and A = 3.0.
In Fig. 6 we observe the 16 density matrix coefficients as a function of time.Both classical and quantum computations are consistent.To further compare the classical and quantum algorithms Fig. 7 presents the expected values of the population corresponding to the lowest and highest energy levels of both atoms P 1,L (t), P 1,H (t), P 2,L (t), P 2,H (t) where This figure also shows the expected values of the spin where I is the 2 × 2 identity matrix.As in the previous example, the coefficients of the density matrix where also confirmed by solving Eq. ( 1) and then projecting onto the elements of h n through Eq. (??).

VIII. CONCLUSIONS
In the present paper, a general framework to linearize the von-Neumann equation was developed.It mainly relies on the projection of the density matrix and the Hamiltonian on an operator base formed by the elements of a Lie algebra.For the particular case of q n the von-Neumann equation is mapped to the conventional Shrödinger-like equation (2), where the state vector |ρ(t)⟩ corresponds to the column stacked matrix elements of the density matrix ρ(t) and the Hamiltonian to the superoperator H = I ⊗ H(t) − H ⊤ (t) ⊗ I.It is shown that this is but one of the multiple ways of linearizing the von-Neumann equation.Other versions can be obtained by projecting onto different algebras.In the case of the Hermitian algebras h n the linearization yields a density matrix vector |ρ(t)⟩ with purely real entries which highly simplifies the quantum tomography   All these notions were used to implement a quantum algorithm that solves the von-Neumann equation.The quantum algorithm was tested against the classical solution of the von Neumann equation for two toy Hamiltonians using the QASM simulator of IBMQ [29] giving identical results.Even though these algorithms were devised for the algebra whose elements are the Pauli strings, they can readily be adapted for any other Hermitian algebra.
We have seen in Eqs. ( 71) and (72) that the structure constants generate the time evolution of the density matrix through a series of unitary quantum gates.This framework thus gives us the flexibility to engineer the structure constants and consequently the required quan-tum gates by choosing the elements of the algebra.
This same scheme can also be applied to the analysis of the time evolution of open-quantum systems through the linearization of the Lindblad-Von Neumann master equation.This problem will be addressed elsewhere.Q i,j of q n ⊗ q n can be related to the expansion in terms of q k ⊗ q l through Expanding A and B in the right-hand side of (A1) we find that Tr q † i q † a q j q b Tr q b ⊗ q a B ⊤ ⊗ A .(A2) In the last step we have used Tr q b B ⊤ = Tr q † b B .Substituting the expansion of q † b ⊗ q † a in terms of the Q r,s as into Eq.(A2) we obtain Only if a and b are given by ( 22) and ( 23), as stated by (21), the term Tr Q † k,l q b ⊗ q a = 1, otherwise Tr Q † k,l q b ⊗ q a = 0. Hence, the only non vanishing term from the sum in (A4) is Plugging ( 18) and ( 19) in to the previous equation, one gets Substituting the explicit expressions for a and b from ( 22) and ( 23) in the equation above yields Introducing this result in (A4) we finally obtain (A1).

Appendix B: Useful properties of the Pauli strings algebra
In this appendix we derive useful properties for the Pauli strings algebra mainly to obtain recurrence relations that make the classical stage of the algorithm more efficient.Pauli strings have the form (1/2) n/2 σ k1 ⊗ σ k2 • • • ⊗ σ kn where σ 0 is the 2 × 2 identity matrix, σ i (i = 1, 2, 3) are proportional to the standard Pauli matrices and the factor 1/ √ 2 is used to normalize the base to unity.Here we have adopted a notation where the superindex n of the element h n i indicates the number of qbits spanned by the algebra.The elements of algebras of higher dimension can be computed by sequentially taking the Kronecker product of the elements of lower dimensional algebras.For example, the elements of the algebra that spans two qbits may be obtained from the Kronecker product of the elements of the algebra corresponding to one qbit as where i = 4(i 1 − 1) + i 2 , i 1 , i 2 = 1, 2, 3, 4 and i = 1, 2, . . ., 16.More generally, the algebra of n + m qbits can be generated from the Kronecker product of the algebras of n and m qbits, i.e., where In particular, a base corresponding to an arbitrary number of qbits can be constructed by recursively applying In this notation, the commutation of two elements of the algebra is given in terms of the structure constants as Using the orthonormailty of this base, the explicit form of the structure constant is given as The structure constants b n i,j,k ensued from taking the anticommutator as the Lie braket allow to express the anticommutator as Although these do not play any role in the construction of the Hamiltonian gates, they will be very useful as auxiliary parameters in determining recurrence relations for c n i,j,k .
From the general commutator relations if follows immediately that the commutator and anticommutator of an algebra of dimension 2 2(n+m) can be expressed in terms of the commutator and anticommutator of algebras of smaller dimensions 2 2n and 2 2m as Now we move on to how to workout recurrence relations for the structure constants.Multiplying both sides of (B15) and (B16) by h n+m k = h n k1 ⊗ h m k2 and using the definitions (B10) and (B12) we obtain the following recurrence relations or more succintly where (C n k ) i,j = c n i,j,k and (B n k ) i,j = b n i,j,k .The use of these recurrence relations to obtain structure constants is far more efficient than the direct application of (B10) and (B12) because it avoids the tracing and all the matrix multiplications required by (B10) and (B12).For instance, one could initially calculate the structure constants for n = 1 (one qbit) through (B10) and (B12), and then, by setting m = 1 in (B19) and (B20) it is possible to recursively compute the commutors and anticommutors of progressively larger algebras.
In a similar fashion we are able compute recurrence relations for the coefficients Tr C k h

Figure 1 .
Figure 1.Main quantum algorithm to determine the dynamics of the density matrix coefficients ρ1(t),. . ., ρn(t).(b) Controlled M (α) gate expressed as a sequence of differential time steps.(c) One time step controlled differential M (dα) gate expressed as the series of Hamiltonian gates generated by the structure constants of the Pauli strings.

Figure 2 .
Figure 2. Controlled M (α) gate expressed as a sequence of differential time steps.

Figure 3 .
Figure 3.One time step controlled differential M (dα) gate expressed as the series of Hamiltonian gates generated by the structure constants of the Pauli strings.

Figure 4 .
Figure 4. Density matrix coefficients as functions of time for the magnetic resonance Hamiltonian.The plots from the classical (continuous lines) and quantum (circles) computations are shown.

Figure 5 .
Figure 5. Expected values of the lowest energy level population, highest energy level population and spin projection along the z axis.The classical and quantum computations are shown as solid lines and circles, respectively.

Figure 7 .
Figure 7. Expected values of the lowest energy level population, highest energy level population and spin projection along the z axis for (a) the first and (b) second atoms.The classical and quantum computations are shown as solid lines and circles, respectively.

2
,m h 2n m .(B24) Then, multiplying both sides by h 2n l and tracing, we haveTr C n k h 2n i 4n = 2 14nterms are needed to prove that.Even for a small number of qbits this is a challenging task.However, (B25) can be recast in the form of a recurrence relation by substituting in it (B19), (B21) and (B22) for which we find Figure 6.Density matrix coefficients as functions of time for the Hamiltonian corresponding to two spin 1/2 particles coupled through exchange interaction subject to a oscillating magnetic field.The classical and quantum computations are shown as solid lines and circles, respectively.
It only remains to prove that in the expansion of the Hamiltonian gate (75) the elements h 2n This can be easily proven by showing that the commutator of the projections of C n k over two generic elements h 2n i and h 2n j vanishes, namely