Quantum computing Floquet energy spectra

Quantum systems can be dynamically controlled using time-periodic external fields, leading to the concept of Floquet engineering, with promising technological applications. Computing Floquet energy spectra is harder than only computing ground state properties or single time-dependent trajectories, and scales exponentially with the Hilbert space dimension. Especially for strongly correlated systems in the low frequency limit, classical approaches based on truncation break down. Here, we present two quantum algorithms to determine effective Floquet modes and energy spectra. We combine the defining properties of Floquet modes in time and frequency domains with the expressiveness of parametrized quantum circuits to overcome the limitations of classical approaches. We benchmark our algorithms and provide an analysis of the key properties relevant for near-term quantum hardware.


Introduction
The interaction between the electromagnetic field and matter is one of the basic principles, which is used to probe and manipulate electronic and magnetic degrees of freedom. Quantum systems that are subjected to time-periodic irradiation show intriguing phenomena such as light-induced surface states [1], topological phases of matter [2][3][4][5][6], high harmonic generation [7] and light-induced superconductivity [8]. With advances in highpower ultra-fast spectroscopy, non-linear phenomena, including multi-photon processes, have opened up new prospects of dynamically controlling quantum dynamics, engineering atomic [9], molecular [10][11][12] and solid-state [13,14] properties on demand. On the theoretical side the study of non-equilibrium dynamics in highly entangled many-body systems is at the forefront of current research [15][16][17][18][19][20][21][22][23]. Light-driven quantum systems are described by time-dependent Hamiltonians H(t), which include the interaction between the external electromagnetic field and the electronic or magnetic degrees of freedom. If the field is periodic in time, i.e., H(t) = H(t + T ) for some period T , one can use the concept of Floquet engineering to control the microscopic degrees of freedom [24]. Floquet engineering is based on the Floquet theorem [25,26]. It states that there exists a complete and orthogonal set of solutions |Ψ(t)⟩ to the time-dependent Schrödinger equation taking the form Here we have chosen the reduced Planck constanth = 1, |Φ α (t)⟩ is the Bloch amplitude periodic in time, and ϵ α denotes the Floquet quasi-energy. The quasi-energies depend on the eigenstates of the non-driven system as well as the specific form of the driving field. They are uniquely defined up to the fundamental frequency Ω = 2π/T . Changing the modulation of the driving field, such as shape, intensity and frequency allows to modify the ϵ α 's and thereby the quasi-energy or Floquet energy spectrum, inducing novel topological phases not present in equilibrium [27]. While the Floquet theorem is a strong statement, computing the quasi-energy spectrum explicitly, depending on the microscopic degrees of freedom and the time-dependence of the external field, remains a difficult problem due to the exponentially scaling Hilbert space. Various classical techniques, such as time-dependent dynamical mean-field theory (t-DMFT) [28][29][30][31][32], time-dependent density matrix renormalization group (t-DMRG) [33], kinetic equations [34,35], perturbative high-frequency expansions [36][37][38][39][40][41][42] and exact diagonalization [43][44][45] have been employed, but unfortunately most of them either are not universally applicable or scale exponentially in system size. Especially the computation of the whole quasi-energy spectrum requires more than a simple time-evolution. Another problem in Floquet-engineering comes from the limited theoretical modeling, in which only a few band model is considered and the perturbative infinite frequency expansion can be applied. In real materials however the situation is more difficult, as higher energy bands become important once the frequency is sufficiently large to induce optical transitions.
In this work, we propose to use quantum computers to overcome the limitations of previous classical approaches. We present two algorithms that in principle can calculate all Floquet eigenstates and quasi-energies. Our quantum-classical hybrid approach takes the limitations of modern Noisy Intermediate-Scale Quantum (NISQ) hardware [46] into account and allows to either improve accuracy by increasing the number of auxiliary qubits or by deepening the quantum circuits. Both algorithms make use of the enhanced expressiveness [47] of parametrized quantum circuits to find approximate Floquet eigenstates and compute the corresponding quasi-energies. While the first algorithm works in the original Hilbert space of the time-dependent problem, the second algorithm makes use of the extended Floquet Hilbert space in frequency domain. After a detailed explanation of the basic principles, we evaluate the performance of both algorithms by investigating the simplest system that is not analytically solvable: a spin-1 2 in a linearly polarized periodic magnetic field.

First Algorithm
In the first algorithm, Fauseweh-Zhu-1, we use a parametrized quantum circuit U Θ , with variational parameters Θ, in combination with the defining properties of the time evolution operator to determine a good approximation for Floquet eigenstates. Without loss of generality, we set t 0 = 0 and denote U t = U (t, 0). The Floquet theorem implies that for a full time period t = T the time evolution operator has complex eigenvalues of the form We represent the Floquet eigenstate |Ψ α (0)⟩ using a parametrized quantum circuit where |0⟩ is the initial state of the quantum computer. If U Θ |0⟩ is a Floquet eigenstate, then holds for the overlap between the time-evolved state and the initial state. Note that overlaps of the form |⟨0|U † V |0⟩| 2 can be efficiently computed on a quantum computer by interpreting them as the measurement of the computational initial state projector Thus we devise a hybrid algorithm that maximizes the overlap to determine a Floquet eigenstate. An optimal solution Θ α has a maximal overlap of 1. To obtain the complete set of eigenstates, we modify the target function using a Lagrange multiplier λ > 0 to project out solutions Θ β that have previously been computed: Here λ has to be choosen sufficiently large for the algorithm to find a new solution, see also [48]. Once a Floquet eigenstate has been determined, we use iterative quantum phase estimation (IQPE) [49][50][51] to compute the complex phase e −iϵαT of U T upon application to the eigenstate. The algorithm is sketched in Fig. 1 a). We turn to a high-level analysis of the algorithm and its requirements. The basis for the algorithm is the optimization of a parametrized quantum circuit. Thus one of the fundamental requirements is that the manifold spanned by all the possible quantum circuits within the ansatz [52] contains the Floquet eigenstates or is at least close to it. Thus carefully choosing a circuit ansatz is important for the applicability of our algorithm, and is currently a subject of intensive research [53][54][55][56][57][58]. Note that the target function L(Θ) allows for a gradient-based optimization using parameter shift rules [59]. Choose initial parameters Θ 3: Update parameters Θ to increase target L(Θ)

7:
Iterative quantum phase estimation on optimized state U Θα |0⟩ for ϵ α 8: The algorithm optimizes a global observable and is therefore, in principle subject to the problem of barren plateaus [60]. This can be avoided by replacing the global observable with a local observable with identical extremum, see also [61].
We do not specify how the time evolution in the first part of the algorithm is performed. Trotterization is in principle applicable to k-local Hamiltonians [62] and is the most accurate. However, it is disadvantageous with respect to circuit depth [63][64][65]. Hence it is best applied if T is the smallest time scale, i.e., in the large frequency limit. Various other methods have been developed [66][67][68][69][70] that make use of shallower circuits, including approaches that work directly within the variational manifold [71,61,72]. All of these methods can be combined with our algorithm, as long as they keep track of the global phase alignment.
The final part of the algorithm uses a single ancillary qubit to perform IQPE. It applies controlled (U T ) 2 n gates to compute the n-th bit of the phase 2πϕ = ϵ α T . The iteration depths n max of IQPE determines the phase resolution and thereby desired precision of the quasi-energies. Unless powers of U 2 T can be obtained within the same circuit depth, e.g., with variational time evolution, IQPE leads to an increasing circuit depth requirement for increasing precision.

Second algorithm
The second algorithm, Fauseweh-Zhu-2, uses Fourier analysis to map the problem onto an extended Hilbert space that is then solved using variational quantum eigensolver approaches. From the definition of the Floquet modes in (2) one can derive that the Hermitian operator Applying a Fourier expansion we can map the time-dependent evolution problem to a time-independent eigenvalue problem Note that the state Φ k α is part of an extended Hilbert space R⊗T containing the original Hilbert space R onto which H(t) acts and the Hilbert space T of square-integrable periodic functions of period T . This extended Hilbert space is not physical, but it has a clear connection to the original Hilbert space: the eigenvalues of the time-independent problem are identical to the original time-dependent problem, up to multiples of jΩ. A graphical interpretation of this procedure is shown in Fig. 1 b), as the original system is expanded into infinitely many copies. The Hamiltonian H 0 acts within the plane, while H ±γ , γ > 0 introduces hopping in the j direction. In the literature, this is understood as extending the system in a new dimension [24]. This is certainly true for non-interacting systems, but not in the strict many-body sense, as with each additional layer in the quantum number j the Hilbert space dimension increases by dim(R), but in a many-body system it would be multiplied by dim(R).
The eigenvalue problem in (11) is the starting point for our second algorithm. We parameterize the Floquet eigenstates Φ k α using a parametrized quantum circuit Notice that we used the notation |0⟩ T |0⟩ R to mark the extended Hilbert space computational basis. Here the first number marks the quantum number j of the T Hilbert space part, while the second number refers to the original Hilbert space R. In the following we neglect the indices T and R of the states. We define the effective Hamiltonian Now we want to compute the eigenstates of H eff to obtain the Floquet energy spectrum. We use the variational principle to optimize the parameters Θ. As the extended Hilbert space is infinite dimensional, we truncate H eff at ±j max . This truncation introduces finite size errors in our calculation. We expect that finite size errors are smallest in the center of the energy spectrum. We therefore minimize H 2 eff instead of H eff , giving us the squared eigenvalues ϵ 2 α . The sign of the eigenvalues can then be determined after optimization by measuring the expectation value of H eff . As in the first algorithm we define our target function using a Lagrange multiplier λ > 0 to project out solutions Θ β that have previously been computed which corresponds to an excited state variational quantum eigensolver [48] for the squared effective Hamiltonian. An overview of the Algorithm is sketched in Fig. 1 b). By analyzing the second algorithm in comparison to the first algorithm, we immediately identify the increased auxiliary qubits that are required due to the extension of the Hilbert space. While this increases the width of the quantum circuit, its depth now depends purely on the depth of the ansatz. This is in striking contrast to the first algorithm, where the maximum depth is determined by the IQPE part of the algorithm, and hence by the required quasi-energy resolution. The second algorithm also has the property of benefiting from a mixed qudit-qubit architecture, as the truncated T part of the Hilbert space naturally leads to states such as |±j⟩ |ϕ⟩, with |ϕ⟩ being a state in the R Hilbert space. This could be useful for quantum computer architectures that have access to more than two states per fundamental quantum building block. The accuracy of the approach also directly depends on the maximum truncation j max , that needs to be determined, depending on the driving and the original Hamiltonian. Choose initial parameters Θ 3: Update parameters Θ to increase target L(Θ)

Benchmark
In the following we evaluate the applicability, performance and scalability of both algorithms. We start with a simple system, that is not analytically solvable, the linearly driven spin- 1 2 , Here σ i are the Pauli matrices. We fix the energy spacing ∆ = 1 and the driving frequency Ω = 2.5, and investigate the Floquet energy spectrum as a function of the amplitude A of the external field. To investigate how the approach scales for larger systems with strong correlations we then perform simulations of the Floquet states for the circularly driven spin-1 2 Heisenberg chain with periodic boundary conditions In equilibrium the spin-1 2 chain is gapless and exhibits fractionalized spinon excitations [73] in the thermodynamic limit. In a finite size systems elementary excitations are best described by gapped S = 1 triplet excitations.

Linearly driven spin-
We use the U 3 gate as a generic single qubit rotation for the parametrized quantum circuit: The time evolution U T = T exp(−i dtH(t)) is performed using Trotterization with 100 time steps. Optimization is performed using conjugate gradient descent on a quantum computer simulator using the IBM qiskit package [74]. The eigenstates are optimized with

Algorithm Fauseweh-Zhu-2
As the second algorithm uses the extended Hilbert basis we transfer the problem in (15) to frequency space with a truncation value of j max = ±1 where the matrix S x is defined in the appendix A. To account for this, we must also enlarge the ansatz in the parametrized quantum circuit and cannot use a single qubit gate. We use a variational Hamiltonian approach [54] to effectively reach the target states, We specify the hermitian generators K l in appendix A. The parameters are optimized using the conjugate gradient descent method. We used 10 4 samples for optimization and evaluation of observables on a quantum computing simulator. For comparison we also computed the exact eigenvalues of the Floquet matrix. The results are shown in Fig. 3. We see a very good agreement between the simulated results and the exact quasi-energy structure of the truncated Hamiltonian. Naturally the truncation leads to an increasing error with amplitude A and only qualitatively captures the exact quasi-energy spectrum to a certain point. Increasing the truncation from j max = ±1 to j max = ±2 significantly

Circularly driven spin-1 2 Heisenberg chain
To evaluate the performance of our algorithms for the spin-1 2 chain we fix the Heisenberg interaction to J = 1 and the field amplitude A = 2. We perform simulation without shot or other noise on a quantum computer simulator. To evaluate the effect of hardware noise we perform simulations assuming a symmetric depolarizing noise channel with details given in Appendix B. As we are interested in Floquet states that still have a well defined energy spectrum, as opposed to essentially random states [75], we choose Ω = 5JN site , where N site is the number of sites in the chain. This scaling of the driving frequency with the number of sites avoids that U (T ) exhibits properties of random matrices belonging to circular ensembles. Of course physically interesting cases in the thermodynamic limit do not necessarily fulfill this criterion, but there are numerical indications that in many cases there exists a critical frequency Ω crit separating regimes of finite and infinite heating [76][77][78][79], so that Floquet spectra are well defined even in the thermodynamic limit.

Algorithm Fauseweh-Zhu-1
We choose an Ansatz |Ψ i (Θ)⟩ for the parametrized quantum circuit that has a layered structure, with each layer containing single qubit gates acting on all sites, which contain the variational parameters Θ, followed by entangling all qubits with CNOT gates. The details are shown in Appendix A. We use the first algorithm to optimize the variational parameters for all quasi-energies. To quantify the convergence of the approach we define the error as the following 2-norm where ⟨Ψ i (Θ)|U T Ψ i (Θ)⟩ measures the overlap between the initial state |Ψ i (Θ)⟩ of the approximate Floquet state with quasi-energie ϵ i and the time evolved state. The results are shown in Fig. 4 for up to 8 sites in the chain. As expected, increasing the number of layers improves on the error of the variational approximation. Interestingly we observe the buildup of a plateau in the reachable accuracy with increasing chain length. Before the breakdown of this plateau the error reduces only slowly before it converges exponentially to zero with the number of layers. This observation is similar to the recently observed scaling of variational quantum circuits representing ground states of systems at a quantum critical point [80]. Our numerical analysis suggests that approximating Floquet eigenstates with variational quantum circuits is a hard problem similar to variationally approximating gapless quantum phases. This indeed makes sense, as a quantum critical point is driven by high-energy fluctuations, while Floquet systems are coupled across the whole (quasi)-energy spectrum. For ground states this behaviour breaks down once the quantum circuit is deep enough to resolve the finite size of the system, in which a finite spectral gap prevails. For Floquet systems a similar mechanism, in which the circuit can resolve the size of Floquet quasi-gaps, could provide an explanation for the breakdown of the plateau scaling.

Algorithm Fauseweh-Zhu-2
The effective Hamiltonian of the driven spin- 1 2 Heisenberg chain in the extended Hilbert space reads, The Ansatz for the parameterized quantum circuit follows the definition in (19). The generators are given in Appendix A. We define the error as the difference between the exact quasi-energy and the variational quasi-energy averaged over all Floquet states, The results, shown in Fig. 5, reveal a similar behavior as the first algorithm but with a slightly increased prefactor. Thus the usage of the extended Hilbert space does not change the scaling of the representability of Floquet states with variational quantum circuits. Due to the increased simulation requirements coming with the qudit-qubit structure we only simulated spin chains with up to 5 sites in this case.

Discussion
In this paper we have presented two quantum algorithms, Fauseweh-Zhu-1 and Fauseweh-Zhu-2, to compute Floquet eigenstates and quasi-energies. Both algorithms use parametrized quantum circuits in combination with a quantum-classical hybrid approach to find solutions to the defining Floquet properties in time-and frequency-domain, respectively. While the precision of the first algorithm depends on the depth of the quantum circuit due to the IQPE application, the precision of the second algorithm mainly depends on the frequency truncation j max and thereby on the width of the parametrized quantum circuit. There are no fundamental hurdles in applying both algorithms on NISQ devices, with the advantage of complementary requirements in circuit depth and width for increasing system size. We have tested both algorithms on a quantum computer simulator for the linearly driven spin-1 2 problem and showed the principal applicability of our approach. To investigate the performance of our algorithms for larger systems we simulated a circularly driven spin-1 2 chain with up to 8 sites. We have demonstrated the feasibility of our variational approach and observed a scaling behavior in the accuracy of both approaches which was previously seen for ground states of quantum critical systems. Our work has uncovered a connection between variationally approximating ground state properties of critical systems and Floquet modes. Exploring the universality of this observation is left for future work. We have also simulated the effect of noise in Appendix B, demonstrating the stability of our variational approach towards this perturbation. As with other variational hybrid algorithms, the ansatz for the parametrized quantum circuit is fundamental for the success of the approach. In this context, the qudit-qubit structure of the second algorithm calls for novel schemes that have so far not been explored. It will be interesting to explore more complicated driving schemes with our approach and to test the performance on real devices using advanced error mitigation methods [63] in the near future.

A Generators for variational Hamiltonian ansatz
A.1 Linearly driven spin- 1 2 To find an efficient variational ansatz for the Fauseweh-Zhu-2 algorithm we first separate the gates required for the T and the R part of the Hilbert space. For j max = ±1 the qudit space is 3 dimensional. A general unitary matrix in SU(3) can be parametrized using the Gell-Mann matrices, which span the corresponding Lie-Algebra. For our particular case we find, that a much more limited set is sufficient to find all Floquet eigenstates. For R a simple Ry gate is sufficient. Naturally for nonzero A the Floquet eigenstates do not separate into the two Hilbert spaces anymore and we need an entangling gate to account for this. We use the matrices to define our variational Hamiltonian ansatz where e iΘ 4 Sxσx is the aforementioned entangling gate. For j max = ±2 we use a similar Ansatz by replacing the matrices S 4 and S 5 with the following generators of the SU(5) In this case, the matrix S x is given as,  The part of the Ansatz acting only on the single qubit remains invariant. A single layer of this Ansatz is sufficient to produce the results shown in Fig. 3.

A.2 Circularly driven spin-1 2 Heisenberg chain
For the driven Heisenberg spin chain we use a layered variational circuit that consists of an initial layer of single qubit gates, followed by a varying number of layers containing a linear chain of CNOT entangling gates and general single qubit rotation gates. A simple example for a 4-site chain with a single layer is shown in Fig. 6. For the second algorithm we use a layered ansatz of the following form, whereΘ is an index function that assigns a unique parameter from the vector Θ to the unitary operator and Γ j are the Gell-Mann matrices [81]. The ansatz therefore consists of single qubit and qudit rotations followed by qubit-qubit entangling gates and finally qubit-qudit entangling gates.

B Noisy simulations
To evaluate the performance of our algorithms in the presence of noise we simulate the second algorithm in the presence of symmetric depolarising noise [82] in the driven spin-1 2 Heisenberg chain. We employ the same variational Ansatz used in Fig. 5 for a chain length of 2 and 3 sites and measure the error given in Eq. (22). The noise channel is implemented into the simulation using the operator Λ l , modifying the resulting unitary evolution in Eq.
The noise operators act only on the qubit part of the system, neglecting decoherence in the qudit subspace as well as correlated noise between the qubits. We choose a symmetric depolarising noise, so that the noise operators are randomly chosen σ x , σ y or σ z matrices acting on a random qubit i. We denote the probability of the noise operator deviating from the identity with p. We also take a finite shot noise into account by averaging over 10 5 simulations. The results for the error are shown in Fig. 7. For the low noise simulated here we see a linear connection between the probability of noise occurring with the error of the circuit independent from the system size. However the overall amplitude of the error indeed increases with the number of qubits. This behavior is similar to previous works investigating the effect of noise on variational quantum circuits [83], demonstrating the applicability of the approach for sufficiently low noise and the scaling behavior of the noise-induced error with the system size.