Variational quantum algorithm for ergotropy estimation in quantum many-body batteries

Quantum batteries are predicted to have the potential to outperform their classical counterparts and are therefore an important element in the development of quantum technologies. Of particular interest is the role of correlations in many-body quantum batteries and how these can affect the maximal work extraction, quantified by the ergotropy. In this work we simulate the charging process and work extraction of many-body quantum batteries on noisy-intermediate scale quantum (NISQ) devices, and devise the Variational Quantum Ergotropy (VQErgo) algorithm which finds the optimal unitary operation that maximises work extraction from the battery. We test VQErgo by calculating the ergotropy of a many-body quantum battery undergoing transverse field Ising dynamics following a sudden quench. We investigate the battery for different system sizes and charging times, and analyze the minimum number of ansatz circuit repetitions needed for the variational optimization using both ideal and noisy simulators. We also discuss how the growth of long-range correlations can hamper the accuracy of VQErgo in larger systems, requiring increased repetitions of the ansatz circuit to reduce error. Finally, we optimize part of the VQErgo algorithm and calculate the ergotropy on one of IBM's quantum devices.


I. INTRODUCTION
The allure of modern quantum technologies relies on leveraging quantum effects such as coherence and entanglement to out-perform their classical counterparts.In recent years this has been motivated by rapid experimental advances which has increased control over quantum states and has allowed to explore fundamental concepts in these devices.In particular, quantum thermal machines allow to explore the foundations of quantum thermodynamics, with devices such as quantum heat engines and refrigerators designed to control work output and heat flow with quantum media [1][2][3].Energy can also be stored in quantum batteries to be extracted at a later time [4][5][6][7][8][9][10][11][12], which have the potential to outperform their classical counterparts in terms of total stored energy [13,14], charging speed [8,[15][16][17][18][19][20] and energy extraction [21,22].The maximum amount of energy that can be extracted from quantum systems through unitary processes is given by the ergotropy [23] which relies on finding the optimal unitary operation which transforms the system to its lowest energy state, known as its passive state.This can be a difficult task as the ergotropy can be sensitive to correlations which can also affect device performance, notably improving efficiency in quantum heat engines coupled to squeezed baths [24][25][26][27], while impairing energy extraction from many-body batteries [6,7,[28][29][30][31][32][33][34][35][36].Simulation of the latter problem will be the focus of our work.
Simulating the dynamics of many-body quantum systems in itself can be a complex problem due to the nonnegligible role of quantum correlations which arise from finite couplings between particles.Furthermore, by today numerical calculations carried out on classical hardware are limited to small numbers of particles.This is in contrast to algorithms based on quantum hardware, which promise to alleviate some of this complexity by simulating quantum wave-functions in the Hilbert space of quantum bits, rather than numerically in classical registers.In addition to quantum physics and other fundamental sciences [37][38][39][40][41][42][43], quantum computers promise applications in various technological sectors including chemistry [44][45][46][47][48][49] and materials design and research [50][51][52][53][54].It is for this reason that they have seen unprecedented growth in recent years: reported milestones include simulation of dynamics and calculations of accurate expectation values on a 127 qubit device [55], demonstration of fast converging quantum-enhanced Markov chain Monte-Carlo simulations [56] and generation of large-scale cluster states on superconducting qubit devices [57].While fault tolerant quantum computation (FTQC) based on error corrected qubits is still not technically possible, currently noisy intermediate-scale quantum (NISQ) processors are available.However, they only have short-lived qubits which are not protected from decoherence [58][59][60].In this NISQ era, quantum algorithms rely on shallow circuits where qubits are measured quickly [61,62].There is therefore a need for quantum algorithms which can simulate quantum systems within a limited time-span while still solving problems which exceed the capabilities of their classical counterparts.
Variational quantum algorithms [63][64][65] (VQA) have been deemed particularly promising for NISQ devices.These are a class of algorithms which can be used to find variational approximate solutions to problems of interest.A famous example is the variational quantum eigensolver (VQE) [61] which is used to determine the ground state of a Hamiltonian through repeated sampling of an ansatz wave-function in the eigenbases of a set of observables.Amongst others, VQAs have been developed to solve the max-cut problem via the quantum approximate optimization algorithm [66], to find numerous chemical properties of molecules [67][68][69], or to perform machine learning tasks like the classification of symmetry protected topological phases [70].While the performance of these algorithms is limited by barren plateaus [71], i.e., the problem of exponentially vanishing gradients with the system size, in recent years tools have been developed to study and mitigate this phenomenon [72,73].In addition, error mitigation has been shown to offer significant improvements when noise is an issue [74][75][76], and approaches inspired by FTQC have resulted in partial error correction schemes developed for NISQ devices [77,78].
In this work we propose a VQA called variational quantum ergotropy (VQErgo) to calculate the ergotropy of a quantum battery on NISQ computers.We use the transverse field Ising spin-chain model to benchmark our algorithm whereby the battery is charged by a sudden quench of the interaction among nearest-neighbor spins.The interactions will create correlations between the spins which results in non-trivial dynamics of the ergotropy.In order to simulate the dynamics we use projected -Variational Quantum Dynamics (p-VQD) [79] to find the timeevolved quantum state and then a variational optimization is carried out to obtain the optimal unitary which prepares the passive state.Although VQErgo has some elements in common with other variational quantum algorithms such as VQE, VQErgo does not find a variational ground state and its energy or other static quantities relating to a Hamiltonian.VQErgo time-evolves the input wavefunction to then evaluate the ergotropy by considering a reduced state of the full wave-function.In contrast to VQE and its variants which variationally prepares the pure eigenstates of a system, VQErgo performs its state optimization on a mixed state which can also be strongly correlated with the rest of the system.
The performance of the algorithm is analyzed for different system sizes and number of ansatz circuit repetitions, and its accuracy is assessed by comparison with exact results.We also analyze how the creation of correlations between spins can negatively affect the ergotropy estimation, requiring an increase in the number of repetitions of the variational ansatz circuit.Finally, we evaluate the effectiveness of our scheme in the presence of noise using a noisy simulator as well as real hardware.Our work represents one of the first NISQ algorithms designed specifically to calculate the ergotropy of quantum systems, expanding the tools already available to describe quantum thermodynamics and associated devices on quantum computers [80][81][82][83].
The manuscript is organized as follows.In Sec.II, we briefly review the operation of quantum batteries and the concept of ergotropy.We also present the VQErgo algorithm which is used to calculate the ergotropy on quantum hardware, separated into four steps: initialization, time evolution/charging, mean energy calculation and passive energy optimization.Then, the transverse field Ising spin-chain Hamiltonian and the charging protocol for our quantum battery are introduced in Sec.III.We describe our main results in Sec.IV, including the dynamics of the system, the measurement of the total energy and the ergotropy from noise-free (state-vector simulations), noisy simulations, and from calculations run on IBM quantum devices [84].Finally in Sec.V, we draw our conclusions and discuss future prospects for this algorithm.

A. The maximal extractable work -ergotropy
We describe a quantum battery made from N identical quantum cells which are charged through unitary dynamics by suddenly switching on an external field V .Initially the battery is prepared in the ground state |Ψ (t = 0)⟩ of a local Hamiltonian H 0 , and during charging it evolves according to The state of the charged battery is therefore timedependent, |Ψ (t)⟩, and energy is discharged by removing the external field V with the total work stored in the battery at time t then given by While this is the total energy that is stored in the entire battery after the time t, it is not necessarily all extractable, especially when only considering subsystems of the device.This would correspond to extracting energy from M ≤ N cells of the battery, which could be required due to a restriction on accessing the full state of the system or in order to only partially discharge the battery.In this scenario energy can be locked in correlations between the M and N − M cells thereby reducing the amount of energy that can be extracted [4,6,7].The maximum amount of work that can be extracted from the M -cell state ρ M = tr N −M {|Ψ (t)⟩ ⟨Ψ (t)|} = j=1 λ j |φ j ⟩ ⟨φ j | (with λ j ≥ λ j+1 ) through unitary transformations is given by the ergotropy [23], which is found by optimizing over all possible unitaries such that the resulting state has the minimum energy with respect to the M -cell Hamiltonian This state is known as the passive state P ρ = i λ i |ψ i ⟩ ⟨ψ i | and no further work can be extracted from it by unitary transformations.The ergotropy can then be expressed in the well-known form [23] where is the projection of ρ M on the eigenstates of H M 0 .To extract energy from the battery we therefore require that p i ̸ = λ i .If the reduced state ρ M is mixed there can be a difference between the work and the ergotropy, W (t) ≥ E(t), which becomes an equality if ρ M is pure.
In classical simulations of quantum batteries, the ergotropy is conventionally calculated by solving Eq. ( 4), i.e., by diagonalizing the sub-system Hamiltonian and the reduced density matrix of the battery state and by computing the relevant overlaps.The analogous way to run this sequence of calculations using NISQ hardware would be to first obtain the full spectrum of the Hamiltonian and its eigenstates [85][86][87].Then estimates of the overlaps with the time-evolved wavefunction can be obtained by measurement in the Hamiltonian eigenbasis.However, we can instead consider the optimization problem in Eq. ( 3) which is naturally expressed in terms of expectation values that can be efficiently computed on a quantum computer.Furthermore, the optimization over unitary operators for the passive state can be naturally phrased in the language of variational quantum algorithms [63][64][65] as we detail in the following section.Hence, current state of the art quantum devices allow us to readily simulate and investigate the ergotropy and other properties of many-body quantum batteries.

B. The Variational Quantum Ergotropy (VQErgo) algorithm
In the following we describe our framework for simulating quantum batteries on quantum hardware and how to extract interesting properties, in particular the ergotropy.We therefore refer to the overall algorithm as the Variational Quantum Ergotropy (VQErgo) algorithm.VQErgo can be divided into 4 subroutines which are (i) battery initialization, (ii) battery charging, (iii) mean energy calculation and (iv) passive energy optimization as shown in Fig. 1.
Battery initialization.The battery starts off in the uncharged state, corresponding to the ground state of the local Hamiltonian H 0 .Any ground state preparation routine (e.g.VQE) can be employed for this task.Note that in the rest of this work we choose the local Hamiltonian to be of the form H 0 = − N i σ z i .Thus, the ground state |0⟩ ⊗N naturally coincides with the initial computational basis state of digital-based quantum computers allowing us to omit a state preparation circuit.
Battery charging.The battery is charged by time evolving the initial state with the Hamiltonian H 1 for a total time t.On digital quantum computers time evolution can be achieved by a Trotter-Suzuki decomposition of the global time evolution unitary into local gates.However, the number of Trotter steps (the number of gates) grows with time t and hence, generally, only short evolution times can be simulated on noisy hardware.To overcome this limitation, there have been several proposals for performing the time evolution variationally using parameterized circuits with a fixed, small number of gates [88][89][90][91].Here, we employ the projected-variational quantum dynamics (p-VQD) algorithm due to its efficiency [79].p-VQD iteratively evolves the parameters w(t + δt) = w(t) + dw of a state ansatz |ψ w(t) ⟩ = U (w(t))|0⟩ in short time increments δt by minimizing the infidelity between the ansatz state |ψ w(t)+dw ⟩ = U (w(t) + dw)|0⟩ and the true time evolved state |ϕ(t The unitary e −iH1δt is typically approximated using the Trotter-Suzuki decomposition with a single Trotter step given that the time step size δt is chosen sufficiently small.Importantly, the repetitions of the state ansatz circuit and of the circuits used for evaluating Eq. ( 6) do not grow with time t.For further details regarding the p-VQD optimization, we refer to Appendix A 1.
Let us note here that our quantum circuit framework for quantum battery simulation is highly modular and, in principle, any of the subroutine algorithms can be exchanged with other viable quantum algorithms for the respective tasks.Specifically, for time evolution, we mention the time-dependent variational algorithm (TDVA) [89,90], and subspace variational quantum simulation (SVQS) [92], which may be used in cases where the number of excited states populated during time-evolution does not exceed the number of qubits.Additionally, one could restrict to quenches with Hamiltonians composed of only commuting terms, e.g., . In this case, the time evolution operator is trivially decomposed into a single layer of two-qubit gates and eliminates the need for involved, approximate time evolution algorithms (see Appendix B).Finally, the recent advances in analog quantum computing provide yet another promising architecture for quantum battery simulations since time evolution is naturally implemented via global Hamiltonians [93,94].
Mean energy calculation.The ergotropy of the quantum battery is calculated as the difference of the mean and passive energies of the charged state |ψ(t)⟩ after tracing over N − M sites (c.f.Eq. 3).Specifically, the mean energy can be expressed as an expectation value We approximate the time evolution unitary UTE = exp(−iH1t) via a variational circuit UTE(ω) optimized using the projected-variational quantum dynamics (p-VQD) algorithm.We also consider the case in which the magnetic field is switched off (h = 0) and the time evolution unitary can be exactly decomposed into a single layer of twoqubit gates.Mean energy measurement: The ergotropy is the difference between the mean and passive energy (c.f.Eq. ( 3)).The former can be measured as an expectation value of H M 0 on M < N qubits of the time evolved state |Ψ(t)⟩.(b) Passive energy optimization: The passive energy is defined as the minimum attainable expectation value of H M 0 over all unitary transformations UE acting on the time evolved state.We express UE (θ) in terms of a variational circuit with parameters θ which are optimized using a typical classical-quantum feedback loop.
of the local Hamiltonian acting only on the subsystem of M < N qubits Passive energy optimization.The computation of the passive energy requires us to find the optimal unitary transformation U E acting on M qubits of the charged state that minimizes the expectation value of the local Hamiltonian within the subsystem We can efficiently perform the optimization over unitaries on current quantum hardware using the tools of variational quantum algorithms.In particular, we define a circuit ansatz U E (θ) composed of two-qubit and singlequbit gates with a set of parameters θ.The optimization of the passive state then amounts to finding the optimal parameters that minimize the expectation value in Eq. ( 8) which can be iteratively achieved by using a classical optimizer like gradient descent and the parametershift rule for evaluating gradients [95,96].
In order to limit the amount of noise during the quantum simulation to a minimum, we employ a hardwareefficient ansatz for the p-VQD and the passive state variational circuits.This means that the ansatz circuits are composed of several layers each containing arbitrary, parameterized single-qubit rotations (decomposed into R Y R Z R Y gates) followed by a series of CNOT gates applied only to neighboring qubits (see Fig. 1(b)).Note that the passive state circuit U E is only defined on the M < N subsystem qubits.Appendix A 2 contains additional details about the circuit optimization.

III. MODEL
To model the quantum battery we consider the paradigmatic transverse field Ising spin-chain [6,22,97,98].The competition between the nearest-neighbour interactions and an external field can result in strongly correlated reduced states ρ M which give rise to a non-trivial dependence of the ergotropy on the charging time t and subsystem size M as we show further below.Moreover, the system is amenable to simulations on currently available NISQ devices for a couple of qubits.The discharged battery is described by the non-interacting Hamiltonian where h > 0 is the external magnetic field, N is the number of spins (cells) in the battery, and σ k i with k = x, y, z denotes the spin-1/2 Pauli matrices.At t = 0 the battery is initialized in the spin polarized ground state ⊗N and thus, naturally coincides with the initial state of the quantum computer.
In order to charge the battery we implement a sudden quench H 0 → H 1 for t > 0 which switches on the nearestneighbor interaction where J is the coupling strength and we consider open boundary conditions.We simulate the time evolved state |Ψ(t)⟩ = exp(−iH 1 t) |Ψ(0)⟩ on a quantum computer using the aforementioned p-VQD algorithm.In Appendix B we provide results obtained with an alternative charging protocol in which the external field is turned off during the charging time (h = 0).In this case, the time evolution operator can be trivially decomposed into a single layer of local two-qubit gates and thus, be implemented without the need for a variational optimization.Throughout the remainder of this work we set ℏ = 1 and consider the transverse field Ising model with fixed parameters h = 0.6 and J = 2.The quench dynamics of the state can be computed directly through where |Φ F j ⟩ and E F j are the eigenstates and eigenvalues of H 1 .The stored work and ergotropy are then calculated using Eqs.( 5) and (4), respectively.Fig. 2 shows the work W and ergotropy E stored in M cells as a function of charging time for a total system size of N = 8 spins.We also plot their ratio E/W describing how efficiently the battery can be discharged which saturates for M = N as expected.One can see that immediately following the quench the work and the ergotropy rapidly increase as the quench drives the system far outof-equilibrium, while subsequent oscillations are the result of finite size effects.We note that the maximum ergotropy and injected work can be achieved at short charging times t ∼ 0.4 for M > 2 and the charging process subsequently stabilizes in the region 2 ≲ t ≲ 6 after which revivals induce further oscillations.
In general, the larger the sub-system M , the more work, ergotropy and efficiency can be achieved.However, for smaller cell size M the efficiency necessarily suffers as correlated cells in the rest of the battery (N − M ) are discarded.This is apparent for times t > 2 when the reduced state ρ M is sufficiently mixed.We also note that, for long intervals for the M = 1 system the ergotropy is exactly zero although the total injected work is non-zero (see panel (b) and (c)).This is related to the equivalence of the reduced and passive states when M = 1, and will be discussed in detail in the next section.

A. VQErgo state-vector simulation
First we simulate the quantum battery using our proposed VQErgo algorithm and analyze the variational optimization in an ideal, noise-free setting via statevector simulation.We restrict ourselves to charging times 0 < t < 1.4 which include the first two maxima of the work and ergotropy curves (see Fig. 2).We simulate the time evolution starting from the polarized product state using p-VQD which optimizes the variational circuit parameters iteratively in small time increments and hence we automatically obtain the evolved states at all intermediate times as well.The ansatz used for the time-evolution and the passive state optimization was the hardware efficient ansatz.All the details regarding the optimization and choice of ansatz are reported in Appendix A 1. In particular, for a given number of spins N we repeat the optimization with different ansatz circuit repetitions, resulting in a different number of variational parameters and compare the p-VQD states to the exact time evolved states in Fig. 9 of Appendix A 1. For each simulated system size we choose a final p-VQD ansatz circuit repetition number that gives rise to small errors with respect to the exact state.
Once the optimized time evolution circuit is obtained we can measure the mean energy on the circuit (see Eq. ( 7)) which allows us to calculate the stored work W . Next, we perform the VQErgo optimization to prepare the passive state on a subsystem of M < N qubits and measure the passive energy from which we can extract the ergotropy E. In Fig. 3 we show the results obtained for a total system size of N = 8 qubits and subsystem sizes M = 1, . . ., 7. We compare the exact ergotropy (orange line) to the ergotropies evaluated on optimized circuits U E of different repetition numbers.Each point represents an average over 100 different runs of the algorithm (i.e., using different random seeds for the circuit initialization).Overall, we find a good agreement of the variationally obtained ergotropies with their exact values given that the number of repetitions is chosen large enough.Note that some of the observed discrepancies have to be attributed to the preceding p-VQD optimization, which also introduces an error in the state.
To further understand the dynamics of the ergotropy in Fig. 3 we examine its constituent time dependent parts, namely the exact probabilities p i of the reduced state and the λ i of its passive state.These are shown in Fig. 4 for both M = 1 and M = 6.The simplest case is M = 1 as its dynamics is that of a two level system (there are only two accessible states), with λ i = p i for times t < 0.4 and therefore the ergotropy is zero.At t = 0.4 there is a crossing in the probability distribution with p 2 > p 1 and finite energy can now be extracted from the battery through reordering of these occupancies (see Fig. 3(a)).For t ≥ 1.2 a subsequent crossing restores the original ordering of the probabilties and thus the ergotropy again vanishes.
This behaviour is echoed in the dynamics of larger subsystems, albeit with more complexity, as the number of occupied states p i and λ i is increased.For instance, the dynamics of λ i for M = 6 possess a similar structure (see Fig. 4(c)) with contributions mainly from the two lowest energy eigenstates.On the contrary, the dynamics of p i is more complex and includes contributions from higher energy eigenstates of H M 0 .This results in a non-zero ergotropy at all times t > 0. In Figs.4(e) and (g) we show λ i and p i at t = 0.4 which corresponds to the maximum ergotropy for M = 6.The large difference between the reduced state and its passive state is readily apparent, as (a), (c) The passive state probabilities λi as a function of t for M = 1 and M = 6, respectively.In (b), (d) we show the corresponding reduced state probabilities pi.For M = 6 we plot λi at (e) t = 0.4 and (f) t = 0.8, and pi at (g) t = 0.4 and (h) t = 0.8.Data in all figures are obtained numerically via exact diagonalization, with the specific times t = {0, 0.2, 0.4, 0.6, 0.8, 1, 1.2, 1.4} denoted by square markers.p i is distributed over all possible states, while λ i is again concentrated around the two lowest eigenstates.However, at t = 0.8 the ergotropy has a local minimum as the p i 's occupy lower energy states owing to less work stored in the battery (see Fig. 4(h)).
Similarly to other variational quantum algorithms, the number of repetitions of the ansatz circuit is an important hyper-parameter of the optimization.Fig. 3 suggests that for VQErgo the required repetitions depends both on the subsystem size M and charging time t.For better visualization, we plot the error in the measured ergotropies as a function of the subsystem size in Fig. 5 at t = 0.4 and t = 0.8.In the case of a single quantum cell M = 1, we always only require one general single-qubit rotation to prepare the passive state.However, with increasing subsystem size M > 1 more layers of singleand two-qubit gates are needed to reduce errors.This is due to correlations that are spread over larger distances within the system which can be quantified through the C XX and C ZZ correlations between qubit i and a second qubit at i + l where ⟨•⟩ = ⟨Ψ(t)| • |Ψ(t)⟩ denotes an expectation value calculated with the time evolved state.We take qubit i = 4 at the center of the N = 8 spin chain as an example and plot its correlations with the other qubits as a function of time in Figs. 6 (a) and (b).The maximum ergotropy coincides with the maximum correlations in the x-directions, while correlations in the z-direction vanish.Furthermore, up to times t ≲ 0.6  5.The absolute error between the ergotropy calculated variationally via a statevector simulation and its exact value versus the battery subsystem size M .We show the error for two distinct charging times t = 0.4 (a) and t = 0.8 (b) for different numbers of repetitions of the passive state ansatz.Note the different order of magnitudes in the error for the two considered times which indicates that the required number of repetitions depends not only on the cell size M , but also on the charging time t.
the qubit is correlated only with its nearest neighbors at ℓ = ±1.We observe a lack of long-range correlations also for the other spins in the chain (not shown here) and can thus infer that a single layer of two-qubit gates (paired with parameterized single-qubit rotations) is sufficient to disentangle all qubits of the subsystem and prepare the exact passive state.However, this is not the case for times t > 0.6 as long-range correlations and entanglement are built up.We therefore require multiple layers of two-qubit gates to rotate the reduced state into the uncorrelated basis set of the passive state and hence, to increase the accuracy of the ergotropy estimation.This is apparent in Fig. 5(b) which shows a significant increase in error at t = 0.8 (note the difference in order of magnitudes between Figs. 5 (a) and (b)).However, we have found that for the particular Ising system considered here, the errors quickly decrease with the number of repetitions and that two repetitions are already often sufficient.Any extra repetitions provide only a small additional advantage which suggests that the number of repetitions of the circuit scales sub-linear with the battery cell size M , making the optimization less prone to barren plateaus [71,99].For charging times t ≲ 0.6 nearest neighbor correlations dominate while at later times t > 0.6 also long-range correlations appear.All data is from exact diagonalization calculations.

B. VQErgo quantum device experiments
Following the analysis of VQErgo under ideal, noisefree conditions, we now evaluate its performance on a real quantum device.To that end, we perform VQErgo on one of the freely accessible quantum computers provided through the IBM Quantum cloud.While the most recent state-of-the-art quantum computers operate on more than 100 qubits and feature small error rates [55], the freely available quantum devices are still small in size and very noisy.Hence, we restrict ourselves to quantum batteries comprised of only a handful of spins that can be simulated with shallow-repetition circuits and can be mapped to the device qubit layout without the need for long-range gates (or SWAP gates).We also substitute our real hardware experiments with noisy, classical simulations that mimic the device noise model.All our experiments are performed on the 7 qubit ibm perth device and its classical simulator analog FakePerth.
Full VQErgo results We start by running the full VQErgo pipeline (including the p-VQD and passive state optimization) on the noisy classical simulator using the SPSA optimizer with 250 optimization steps and 2048 shots per measurement.The training curves and any further technical details can be found in Appendix A 2. We report the final measured work and ergotropy as a function of the charging time for a system with N = 2 and M = 1 battery cells in Fig. 7.Each point is again an average over 100 independent runs of the algorithm.For this small battery system, the ergotropies are in good agreement with their exact values.Any discrepancies and the increased standard deviation compared to the statevector simulation can be attributed to various error sources, such as shot noise, state preparation and measurement (SPAM) errors, coherent and incoherent noise.Note that the observed error in the work and ergotropy increases slightly with time which is to be expected since p-VQD iteratively evolves the ansatz state in time and hence, errors naturally build up in the charged state.
Noise free p-VQD optimization Using the noisy simulator has not allowed us to obtain converged results for p-VQD with N = 4.This can be understood from the fact that p-VQD carries out a variational optimization for each time-step that is simulated.Therefore any error resulting from the noisy hardware or a simulation of it compounds with every iteration.As explained in Appendix A 1, when N = 4, the repetition numbers of the p-VQD ansatz circuit must be increased to twice what it is when N = 2 while still needing 14 iterations.Although circuits of this number of repetitions can be run without the resulting quantum state decohering totally, the accumulated error due to noise is too great to result in accurate time-evolution.However, with the optimized p-VQD parameters from state-vector simulations we are still able to show convergence of the passive state optimization using the noisy simulator and the actual device.In the remainder of this work we perform the p-VQD optimiza- The exact stored work W (dashed line) and ergotropy E (solid line) are plotted against the charging time t for a system with N = 4, M = 2.The results were obtained using state-vector optimized p-VQD parameters together with a noisy simulation of the ibm perth backend as well as actual device results.For the real-device experiments only a single run of the passive state optimization is carried out, while we display the average and standard deviation of 100 runs in the case of the noisy classical simulator.Due to the limited available quantum computing time, we performed the VQErgo optimization on ibm perth only on a subset of the shown times.
tion using the classical statevector simulation and only run the optimized time evolution circuit on the quantum device followed by the variational passive state optimization.Note that in Appendix B we consider a simplified charging protocol that does not require variational time evolution and thus, can be executed end-to-end on current hardware.
In Fig. 8 we display the real-device and noisy simulator results obtained for N = 4 and subsystem size M = 2.The measured injected work is in nearly perfect agreement with their theoretically computed values with the exception of the point at time t = 0.6 which is slightly lower.We believe that this outlier can be attributed to naturally occurring fluctuations in the experimental hardware over time (our simulations have been performed over several weeks).On the other hand, the evaluated ergotropies are consistently lower than their exact values for the noisy simulator and the real-device experiments.We expect this to be the result of decoherence as the passive state circuit contains two additional layers of CNOT gates compared to the circuit the work was measured on.It would be interesting to investigate whether quantum error mitigation such as zero-noise extrapolation can improve these results [100][101][102].However, despite the small discrepancies, the qualitative dependence of the ergotropy with time can be successfully inferred.Importantly, VQErgo allows us to determine the time at which the ergotropy becomes maximal which is crucial for designing many-body quantum batteries that perform optimally.

V. CONCLUSIONS
In this work, we have studied the ergotropy -the maximal extractable work -of quantum batteries.We have shown that the calculation of the ergotropy can be naturally phrased in terms of a variational quantum algorithm and thus, the ergotropy is readily amenable to current NISQ device computations.We have embedded the ergotropy calculation in an end-to-end variational simulation routine for quantum batteries called VQErgo that includes battery initialization, charging, and the ergotropy estimation.Note that due to the modularity of the presented framework, different algorithms for any of its subroutines like initial state preparation or time evolution can be chosen and adapted to the system at hand.
We tested VQErgo on a battery undergoing transverse field Ising dynamics.To that end, we investigated the passive state optimization and the required number of repetitions of the variational circuit with the battery cell size and charging time.In particular, we showed that more than one repetition is necessary beyond a critical charging time after which long-range correlations set in.Subsequently, we demonstrated VQErgo using a noisy classical simulator with noise characteristics from a current IBM quantum device, and demonstrated that the passive state optimization can be carried out on the actual physical device.In both cases, we were able to successfully measure the injected work and for different charging times.While the estimated ergotropies were slightly below their true exact values, the qualitative dependence of the ergotropy with time still matched the theoretical predictions.In particular, the results allow us to infer the optimal charging time of the quantum battery that leads to the maximal ergotropy.
While we only consider quench dynamics in this work, our algorithm could also be expanded to also optimize the charging protocol for a specific charging time, allowing to time-dependently tune coupling terms in the Hamiltonian to maximize the ergotropy on short times.Our algorithm is also not restricted to the simulation of quantum batteries, but is also amenable to other thermodynamic devices where correlations can affect work statistics and therefore depend on accurate estimations of the ergotropy to account for all energy contributions.For instance, in quantum heat engines coupled to nonclassical environments, such as squeezed baths [24][25][26][27], both passive thermal energy and ergotropy can be transferred to the system.Work can also be extracted from quantum heat engines by coupling to quantum batteries or so called quantum flywheels [103,104], where again it is important to distinguish between the amount of useful work quantified by the ergotropy and thermal fluctuations that can degrade efficiency.Finally, due to its sensitivity to correlations the ergotropy can also be used to measure genuine multipartite entanglement [105], as the mixed state ρ M is entangled with the remaining N − M cells leading a non-trivial gap between the ergotropy and the average work.As the passive state is diagonal in H M 0 the entanglement spectrum λ i can be recovered, provided that the eigenbasis can be resolved on quantum hardware [106].
In this work we have shown a viable path towards studying many-body quantum batteries using quantum hardware.It is remarkable that the non-trivial dynamics of the transverse field Ising model can be probed even on the relatively noisy 7 qubit ibm perth device as this is a device with a quantum volume of only 32 [84,107].In contrast, the state of the art Falcon r10 device is reported to have a quantum volume of 256 [108].The complexity of the types of systems that can currently be interrogated with VQErgo is therefore expected to exceed the capabilities we have shown here.As a concrete example, it is feasible that the p-VQD optimization, which we had to carry out using a state-vector simulator for N = 4, could be carried out on a quantum volume 256 device.
We optimize p-VQD using the BFGS optimizer for the state-vector simulations and SPSA [110] for the noisy simulations.While SPSA only performs approximate gradient descent, it also only requires two circuit evaluations per optimization step independent of the number of parameters and can thus be efficiently executed on quantum devices.Moreover, the stochasticity in the perturbation directions make it robust to noise.The fidelity in Eq. ( 6) is evaluated via sampling and can be replaced by a local cost function to make the optimization less prone to barren plateaus [79,99].As a termination condition for the BFGS optimizer we set a precision goal of 10 −6 in the cost function, i.e., the infidelity in Eq. ( 6).When using SPSA we set the number of optimization steps per time step to 1000 instead.We run several state-vector simulations of the evolution up to a total time t = 1.4 with different time increments δt, ansatz repetitions, and 3 random seeds to determine the optimal hyper-parameters that minimize the infidelity with respect to the exact state computed using the QuSpin package [111].In Fig. 9 we plot the infidelity as a function of time for 4 different system sizes showing only the best out of each run.As expected, we find that the infidelity on average increases with time as errors build up in the time evolved state.Furthermore, the infidelity also grows with system size and as such we require more ansatz repetitions to faithfully represent the increasingly correlated quantum states.For the N = 8 state-vector simulations of Section IV A we use 5 repetitions of the hardware efficient gate-fabric, for the N = 2 noisy simulations of Fig. 7 we set the number of repetitions to 1 and for the N = 4 quantum device experiments of Fig. 8 the number of repetitions is equal to 2. Note that in the latter case, we used the pre-optimized p-VQD parameters from the state-vector simulation to time-evolve the state on the noisy hardware.Additionally, we experimented with different numbers of time-steps and found that the optimal number of time-steps needed is 7 for N = 6 and N = 8 while it is 14 for N = 2 and N = 4.

Passive state optimization
Analogous to p-VQD, we use the BFGS optimizer for the statevector simulations and SPSA for the noisy and real-device simulations.With the exception of the hardware experiments for which only a single data point per time is collected, we repeat each classical simulation with 100 random seeds and take the average.Fig. 10(a) shows the average number of required BFGS optimization steps to reach a precision of 10 −6 in the cost as a function of the subsystem size M for an example of a charged state Qubit T1 (µs) T2 (µs) f (GHz) δ (GHz) ϵ readout p(0|1) p(1|0) Single-qubit error CNOT error tgate (ns) 0 188.6   at t = 0.4.The number of repetitions was fixed to 2. Note that the number of parameters in the circuit ansatz grow linearly in M .We also display the observed standard deviation of the ergotropy versus the subsystem size (see Fig. 10(b)).
In Fig. 11 we show two training curves of the passive state optimization using SPSA that were collected for the two exemplary charging times t = 0.5 and t = 0.9 of Fig. 7 in the main text.The dark (bright) color corresponds to the mean (std) over 100 independent noisy simulations using the FakePerth backend while the dashed line indicates the theoretically exact ergotropy.The optimization usually converged within the first 50 ∼ 100 steps.However, the final value can deviate from its ex-act prediction due to noise coming from the mean and passive energy measurements as well as errors arising in the p-VQD time evolution.
We also show two typical training curves for the realdevice experiment results (see Fig. 8 in the main text) performed on ibm perth in Fig. 12. Rather than measuring the passive energy after each optimization step which would require additional circuit evaluations, we instead estimate its value by averaging the two expectation values used by SPSA at each iteration.

Appendix B: Alternative charging protocol
In this section we provide results obtained with a simplified battery charging scheme that does not require trotterization or variational optimization and thus can be easily implemented on NISQ hardware.Instead of evolving the system with the transverse field Ising Hamiltonian of Eq. ( 10) we turn off the magenetic field during the quench and only evolve with the term containing the nearest-neighbor coupling Note that H 1 is composed of only commuting terms.Therefore, the time evolution operator e −iH1t can be exactly decomposed into a single layer of two-qubit gates acting only on neighboring spins leading to where x and θ = −2Jt.
Note that the charging time t is encoded as the angle of the R XX gates.

Statevector simulation
We run VQErgo on charged states of an N = 10 spin system and show the achieved final ergotropies for different subsystem sizes M in Fig. 13.Interestingly, we find that a single repetition is sufficient to achieve a high accuracy with the exactly computed values (orange line) for all considered cell sizes and charging times.This suggests that the evolved state only contains nearest neighbor correlations irrespective of the charging time t.Moreover, we observe that the stored work and ergotropy coincide at times t = (k + 1/2)π/J, with k = 0, 1, 2, . . .reaching an energy W = E = 2h.At these times the charged state is in a fully disentangled product state.The ergotropy in this case depends solely on the number of qubits M < N from which we want to extract the energy instead of the full system size N .This is shown in Fig. 14 for N = 6 whereby the ergotropy and work are the same as in the N = 10 case in Fig. 13.Note that for the case of M = 1, it is also possible to obtain an analytical expression for the ergotropy E(t) = 0, if tan 2 (Jt) ≤ 1 2h sin 2 (Jt) − cos 2 (Jt) , if tan 2 (Jt) > 1 .
(B3) Overall, and unsurprisingly, the dynamics in this case is more trivial than the one generated by the full transverse field Ising Hamiltonian.The time evolved state is entangled only over short distances and thus, the passive state can be prepared with at most a single layer of nearestneighbor two-qubit gates.It would be interesting to study the quantum battery with charging protocols inter-polating the simplified case discussed here and the Ising dynamics from the main text by applying a small number of quenches with alternating non-commuting generators.

Noisy simulations
We have also tested VQErgo with the simplified charging protocol on noisy simulators and hardware.We choose 6 qubits of the 7 qubit ibm perth device and plot an average of the measured ergotropies for different subsystem sizes M obtained on the noisy classical simulator in Fig. 14.The number of repetitions of the passive state ansatz circuit is 1.However, extra SWAP gates are required to map the full circuit to the underlying topology of the real device which ultimately introduces more noise.The error between the exact and variationally obtained ergotropies grows with the battery cell size.On the other hand, the error is independent of the charging time which is in contrast to the p-VQD based simulation where errors naturally built up over time.Despite the noisy values, we can successfully infer the qualitative dependence of the ergotropy with the charging time.
Finally, we also report two results obtained on the ibm perth quantum device for a system with N = 2, M = 1 in Fig. 15.Again, we find an overall good agreement between the measured values and their exact prediction.

FIG. 1 .
FIG. 1.(a) Schematic depiction of the Variational Quantum Ergotropy (VQErgo) algorithm.Initialization: The uncharged battery state |ψ(0)⟩ given by the ground state of the local Hamiltonian H0 is prepared on the quantum device e.g. using the variational quantum eigensolver (VQE).Time evolution/charging: The battery is charged by time evolving the state with another Hamiltonian H1.As an example, in this work we consider quenches with the transverse field Ising Hamiltonian H1 = −J N −1

FIG. 2 .
FIG. 2. (a)The total injected work W , (b) the ergotropy E and (c) the efficiency of the battery E/W as a function of the charging time t.The total system is comprised of N = 8 spins while each line corresponds to a different sub-system size M from which energy is extracted.

4 FIG. 3 .
FIG. 3.The ergotropy E as a function of the charging time t for different subsystem sizes M and a total system size of N = 8.The grey dashed line denotes the work W stored in the battery cells.The orange line corresponds to the value of the ergotropy computed numerically from exact diagonalization (ED) calculations while the markers show the variationally obtained ergotropies for different repetitions of the passive state ansatz.Each point is an average over 100 optimization runs using a statevector simulation (the standard deviation is indicated by the error bars).
FIG.5.The absolute error between the ergotropy calculated variationally via a statevector simulation and its exact value versus the battery subsystem size M .We show the error for two distinct charging times t = 0.4 (a) and t = 0.8 (b) for different numbers of repetitions of the passive state ansatz.Note the different order of magnitudes in the error for the two considered times which indicates that the required number of repetitions depends not only on the cell size M , but also on the charging time t.

FIG. 6 .
FIG.6.The Pauli-XX (a) and ZZ (b) correlations in the charged state between the 4th qubit at the center of the chain and a qubit ℓ sites apart as a function of the charging time t.The insets show the respective correlations as a function of the distance ℓ for two specific times t = 0.4 and t = 0.8.For charging times t ≲ 0.6 nearest neighbor correlations dominate while at later times t > 0.6 also long-range correlations appear.All data is from exact diagonalization calculations.

FIG. 7 .
FIG. 7. The total amount of stored work W (dashed line) and the ergotropy E (solid line) as a function of the battery charging time t computed via exact techniques.We also show their values extracted via VQErgo using an ideal statevector simulation and a noisy classical simulation with FakePerth which mimics the ibm perth quantum device.The results are an average over 100 independently run optimizations.The error bars indicate the corresponding standard deviation.

8 FIG. 9 .
FIG. 9. Infidelity between the p-VQD optimized state (via ideal, noise-free statevector simulations) and the exact time-evolved state for system sizes N = 2 (a), N = 4 (b), N = 6 (c), N = 8 (d), and different number of repetitions of the variational circuit.The optimizations have been performed using BFGS for the Ising chain dynamics defined in the main text.

FIG. 10 . 9 FIG. 11 .
FIG. 10.(a) The average number of BFGS iterations required to achieve a final precision of 10 −6 in the cost of the passive state optimization as a function of the subsystem size M .(b) The corresponding standard deviation in the ergotropy over 100 runs.The optimizations were performed using the noise-free statevector simulator, a total system size N = 8, number of repetitions of 2, and a charging time t = 0.4.

4 FIG. 12 .
FIG.12.The estimated ergotropy E at each SPSA iteration for charging times t = 0.2 (a) and t = 0.4 (b) measured on the ibm perth quantum hardware (see also Fig.8in the main text).The black dashed line corresponds to the exact result calculated using ED methods.N = 4, M = 2.

9 FIG. 13 . 5 FIG. 14 .FIG. 15 .
FIG.13.The ergotropy E as a function of the charging time t for different subsystem sizes M (total system size N = 10).The transverse magnetic field is turned off during the charging quench.The orange line indicates the exact values while the markers denote average values of 100 optimization runs using a statevector simulation with different repetition numbers.The total injected work is shown as a gray dashed line.

TABLE I .
Further information about the ibm perth device used in this work.The table shows from left to right column: relaxation time T1, dephasing time T2, frequency f , anharmonicity δ, readout assignment error ϵ readout , the probability p(0|1) of measuring 0 if the qubit was prepared in the |1⟩ state, the probability p(1|0) of measuring 1 if the qubit was prepared in the |0⟩ state, single-qubit gate error, average CNOT gate error, average gate times tgate.The final row shows the median computed over all qubits.Note that our experiments were performed over the course of several weeks, while the calibration data is only a snapshot of the device properties at a single point in time (data from Nov 13, 2023).