Simulating Static and Dynamic Properties of Magnetic Molecules with Prototype Quantum Computers

Abstract

Magnetic molecules are prototypical systems to investigate peculiar quantum mechanical phenomena. As such, simulating their static and dynamical behavior is intrinsically difficult for a classical computer, due to the exponential increase of required resources with the system size. Quantum computers solve this issue by providing an inherently quantum platform, suited to describe these magnetic systems. Here, we show that both the ground state properties and the spin dynamics of magnetic molecules can be simulated on prototype quantum computers, based on superconducting qubits. In particular, we study small-size anti-ferromagnetic spin chains and rings, which are ideal test-beds for these pioneering devices. We use the variational quantum eigensolver algorithm to determine the ground state wave-function with targeted ansatzes fulfilling the spin symmetries of the investigated models. The coherent spin dynamics are simulated by computing dynamical correlation functions, an essential ingredient to extract many experimentally accessible properties, such as the inelastic neutron cross-section.

1. Introduction

Molecular nanomagnets (MNMs) are finite-size spin systems that have attracted great interest, both for their potential technological applications in quantum computing [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27] and high-density data storage [28,29,30,31,32,33] and as model systems to investigate fundamental quantum mechanical phenomena, such as the tunneling of the magnetization [34,35,36,37], or of the Neel vector [38], frustration [39,40,41,42] or decoherence [43]. In particular, one-dimensional arrays of spins, coupled by anti-ferromagnetic exchange interactions, have been widely studied to evidence peculiar behaviors [44] and to shed light on finite-size [45] and parity [39,46] effects. The former distinguishes between chains (characterized by open boundary conditions) [45,47,48], as opposed to ring-shaped molecular structures, in which periodic boundary conditions yield new symmetries and a peculiar structures of the eigenstates [39,40,41,42,44,49,50,51].

The theoretical modeling of MNMs is often challenging, due to the large size of the Hilbert space, which is difficult to handle by classical computers. Indeed, the dimension of the Hilbert space (roughly scaling as the number of classical bits of information needed to store the system wave-function) increases exponentially with the number of spins in the examined system. For instance, in the case of Mn${}_{12}$, which can be considered the forefather of MNMs, such a dimension is ${10}^8$ [29], even by including only spin degrees of freedom. This makes the exact determination of the ground state wave-function and of its dynamical properties very hard for classical computers. To solve this issue, physicists have so far resorted to different kinds of approximations, which however often require a detailed knowledge of the investigated system, or at least of the hierarchy between the interactions at play.

The advent of quantum computers (QCs) [52,53,54] offers a completely new perspective on this hurdle. Indeed, an intrinsically quantum device can simulate the static and dynamic properties of another quantum system, object of our investigation (called target in the following), using a number of quantum-bits scaling only linearly with the number of units in the target system [55,56,57]. This provides an impressive gain compared to the exponential scaling of resources of a classical device and would allow a brute-force solution of many quantum problems, currently intractable even by using the best available super-computers, without needing an in-depth knowledge of the target system. Hence, although current QCs are still noisy [53], the incessant improvement of their capabilities could make them, in the near future, the ideal architecture to theoretically describe MNMs. In fact, a steady growth in the size of prototypical programmable digital quantum devices was achieved over the last few years with the currently most advanced technological platforms, namely superconducting circuits [58,59] and trapped ions [52]. In parallel, their performances as measured, e.g., by the quantum volume metric [60,61], have also been consistently improving, thus pointing towards promising new avenues of research. Furthermore, several error mitigation strategies have been recently put forward [62,63,64,65,66,67,68,69,70,71], showing great potential in significantly enhancing the quality of the results for quantum computing applications in chemistry and physics.

Here, we focus, in particular, on anti-ferromagnetically coupled chains and rings of spins 1/2, and we show that QCs can be used to efficiently determine the ground state of these systems, together with related observable quantities, and to investigate the dynamics. Besides demonstrating interesting physics, these finite-size spin systems are classically hard to simulate with increasing the system size [72]. Hence, they are also ideal test-beds of current available quantum chips, such as the IBM Quantum superconducting processors [59].

In order to find the ground state of the target system, we use the variational quantum eigensolver (VQE) method [73,74,75]. This belongs to the family of variational quantum algorithms [76,77,78], which collectively represent promising near-term quantum solutions for a wide range of problems, from optimization [79] to machine learning [80,81]. Although the VQE was originally introduced [73] and is nowadays widely applied [75,82,83,84] to treat quantum chemistry problems, it can easily be extended to deal with many-body [85] and spin models [75]. The latter, in particular, can be naturally mapped on a qubit-based hardware.

The specific problem under investigation is made non-trivial by the high degree of entanglement which is typically associated to the ground states of MNMs. In the following, we describe the use of tailored quantum circuits leveraging the intrinsic symmetries of the system, such as the rotational symmetry of the Heisenberg model Hamiltonian [86,87], to realize effective and scalable computations.

In addition, we show that QCs can be used to compute the dynamical properties of MNMs. We focus, in particular, on dynamical spin–spin correlation functions, which constitute the essential (and classically hard to compute) ingredient of many physically crucial observables. This is the case, for instance, of the inelastic neutron magnetic cross-section, which is probably the best experimental technique to probe spin excitations in MNMs [29,88]. As shown in Ref. [89], dynamical spin correlation functions can be obtained by a quantum simulation of the target system time evolution, which can already be reliably computed on small finite-size spin chains. We also discuss here the application of this approach to spin rings.

In the following, we present simulations using realistic noise models of the VQE algorithm to determine the ground state of spin chains and rings with a variable number of sites (from 4 to 6) and experiments on IBM Quantum hardware to compute dynamical correlation functions.

2. Determining the Ground State of Heisenberg Chains by VQE

In this work, we focus on MNMs which can be modeled as chains of spin 1/2 systems linked via pairwise Heisenberg interactions with a uniform strength J and subject to a uniform external magnetic field B. The target Hamiltonian is therefore expressed as

$$\mathcal{H}=2J\sum _i{\mathbf{s}}_i·{\mathbf{s}}_{i+1}+B\sum _{i=1}^N{s}_i^z$$

(1)

with $s_i^{\alpha }$ spin 1/2 operators and open or periodic boundary conditions. As mentioned in the Introduction, the first problem that we tackle is the reconstruction of the ground state wavefunction, and observable properties of these systems by implementing variational methods (VQE) on quantum computing platforms. These techniques are particularly well suited for the study of spin $1/2$ models, which can be directly mapped on quantum hardware by associating each spin to a single qubit, i.e., the elementary digital unit of quantum information.

2.1. Variational Quantum Eigensolver

In a VQE calculation, the exact ground state $|{\mathsf{\Psi }}_0\rangle$ of a target Hamiltonian $\mathcal{H}$ is approximated, in agreement with the variational theorem of quantum mechanics, by constructing a parametrized quantum state $\left|\psi \right(\theta )\rangle$ with the aim of minimizing the energy expectation value ${\langle \mathcal{H}\rangle }_{\theta }=\langle \psi \left(\theta \right)\left|\mathcal{H}\right|\psi \left(\theta \right)\rangle$. A quantum register is used to prepare the trial wavefunction $\left|\psi \right(\theta )\rangle$ through a sequence of unitary operations $U\left(\theta \right)$ (i.e., gates) and to directly measure ${\langle \mathcal{H}\rangle }_{\theta }$. An updated set of parameters ${\theta }^{\prime }$ is then obtained with the help of a classical minimization routine, and is fed back to the quantum circuit, leading to a new energy estimate. Iterating such quantum-classical procedure until convergence, an optimum point can eventually be reached, thus achieving $|\psi \left({\theta }_{opt}\right)\rangle \simeq |{\mathsf{\Psi }}_0\rangle$.

The performances of a VQE implementation, both in terms of accuracy and computational cost, crucially depend on the choice of the parametrized trial wavefunction to be employed. In practice, this is usually reflected in the structure and properties of the so-called variational ansatz, namely the unitary transformation $U\left(\theta \right)$ that prepares $\left|\psi \right(\theta )\rangle$ from some fixed reference state of the quantum register, i.e., it implements the transformation $\left|\psi \right(\theta )\rangle =U(\theta \left)\right|0\rangle$. This operation must ideally be designed in such a way that, for a suitable choice of the parameters $\theta$, the exact ground state $|{\mathsf{\Psi }}_0\rangle$ can either be represented exactly or, more often, be approximated with sufficiently good precision. It is worth noticing at this point that, although in principle, any unitary operation, and therefore any quantum state, can be realized with a sufficiently complex quantum circuit, the number of parameters and of elementary quantum operations is, in general, exponential with the number of qubits, i.e., with the size of the problem under study. A well designed variational ansatz is therefore one which is able to explore only the relevant portions of the Hilbert space, namely the ones where $|{\mathsf{\Psi }}_0\rangle$ or its approximations are located, by requiring only (at most) polynomial resources.

There are essentially two approaches that one can follow to construct a suitable $U\left(\theta \right)$. On the one hand, when no specific information or insight about the structure of the problem is available, an heuristic approach can be adopted. This approach is also particularly well suited for near-term implementations, where the required unitary gates must be adapted to the hardware constraints of the processors in terms of, e.g., native elementary operations or connectivity between qubits. On the other hand, as we will show in the following, one may adopt specialized solutions tailored for the specific task: in quantum chemistry applications, this is exemplified by the unitary coupled cluster (UCC) ansatz [82], and in more general cases, may be represented by symmetry adapted [86,87], optimal control [90] or Hamiltonian-inspired [91,92,93] ansatzes.

In the heuristic case, the variational ansatz is usually constructed in a modular way, by decomposing $U\left(\theta \right)$ in a series of layers. Each layer introduces additional parameters and, in principle, extends the range of states which can be reached. A very convenient choice is obtained by combining single qubit rotations, whose angles represent the variational parameters, and two-qubit operations which introduce entanglement. The most general form of a single-qubit operation acting on the q-th qubit is represented by the matrix

$$U_3^{\left(q\right)}(\theta ,\varphi ,\lambda )=\left(\begin{array}{cc}cos(\theta /2)& -e^{i\lambda }sin(\theta /2)\\ e^{i\varphi }sin(\theta /2)& e^{i(\lambda +\varphi )}cos(\theta /2)\end{array}\right).$$

(2)

The entangling blocks are instead constructed as a product of controlled-NOT operations on different pairs of qubits. Notice that the use of $\mathrm{CNOT}$ operations is particularly appropriate for IBM Quantum processors, where these represent the native entangling gate that is realized at the hardware level. This is in fact one example of the potential practical advantages brought by heuristic ansatzes, which can typically reduce the overheads associated with the implementation of VQE calculations on real devices. These overheads may be represented by, e.g., the additional circuit depth required to decompose a given unitary transformation into the native universal set of hardware gates. However, these convenient features must be balanced with the disadvantages associated with the blindness of the ansatz to the properties of the target Hamiltonian. This typically results in an increased number of variational parameters, which can quickly become hard to handle in the classical optimization step [94,95].

2.2. Heisenberg Spin Chains and Adapted Ansatz

We begin our investigation by considering the case of an MNM taking the shape of a closed Heisenberg square. We thus specialize the general target Hamitonian in Equation (1) by assuming periodic boundary conditions.

In this regime, it is easy to show that the following commutation relations hold, independently from the value of the external field B:

$$[\mathcal{H},S^2]=0\mathrm{and}[\mathcal{H},S_z]=0.$$

(3)

Thus, the usual quantum numbers S and $m=-S,-S+1,\dots ,S$, associated to the total spin ($\overrightarrow{S}={\sum }_i\overrightarrow{s_i}$) and to its z projection respectively, can be assigned to the Hamiltonian eigenstates. As shown, for example, in Figure 1a,c, for increasing values of $B/J$, several crossings between energy levels are present in the spectrum of $\mathcal{H}$. In particular, the ground state changes around $B/J=2$ and after $B/J=4$, thus identifying three different regions of interest, corresponding to $S=0,1,2$, respectively (and $m=-S$). For VQE applications, we will focus only on low- and intermediate-field conditions, since at high values of B, all the spins align in a trivially factorized ground state of the form $|\downarrow \downarrow \dots \downarrow \rangle$.

For the heuristic ansatz (HA) approach, we make use of parametrized quantum circuits $U_{HA}\left(\overrightarrow{\theta }\right)={\prod }_{j=1}^L{V}_j\left({\overrightarrow{\theta }}_j\right)$ whose building blocks are of the so called $R_y$-CNOT form

(4)

Here, only nearest neighbours entangling CNOT operations are allowed and $R_y\left(\theta \right)=e^{-i{s}_y\theta }=U_3^{\left(q\right)}(\theta ,0,0)$ is a Pauli rotation around the y axis. Note that, since the target Hamiltonian is real symmetric, the ground state wave-function must also be real. Hence, compared to the general form of Equation (2), we can remove the parameters $\varphi$ and $\lambda$ (which control relative phases) from each $U_3$ and only leave a single variational parameter per rotation. Moreover, all qubits are intialized in state $|0\rangle$ at the beginning of the computation, and a final set of parametrized $R_y\left(\theta \right)$ is also appended at the end of the circuit.

In Figure 1, we report the results achieved by optimizing the set of free parameters with a numerical simulation performed with the Qiskit statevector_simulator, i.e., by evaluating the results of quantum circuits directly, via algebraic matrix multiplication and by using the COBYLA [96] classical minimization routine. An ansatz depth $L=2$, corresponding to 12 independent classical parameters and 6 CNOT operations, allows us to obtain good results in the intermediate $B/J$ region. At low field, a more complex heuristic circuit ($L=3$, with 16 free parameters and 9 CNOTs) is required to exactly reconstruct the true ground state. This is a direct consequence of the higher degree of entanglement present in the ground state when the Heisenberg interaction is the dominant contribution. We see here that, when working with heuristic ansatzes, the only natural way to enrich the space of possible trial wavefunctions (e.g., when complex ground states must be approximated) is to increase the number of layers L. However, this strategy may not be optimal in view of scaling the problem size, particularly when the classical cost of the optimization problem associated with VQE implementations is also taken into account. On the one hand, only relatively shallow circuits may be run in practice in a near term scenario, where all quantum operations are subject to environment noise and limited fidelity (see also a discussion on noisy simulations below). On the other hand, the number of classical optimization steps, and hence of circuit executions required to estimate the variational energy ${\langle \mathcal{H}\rangle }_{\theta }$, may become quickly impractical when many classical parameters, intially set at random values, are present in the ansatz.

A more viable solution arises when the ansatz design is adapted to the properties of the underlying physical problem. In our case, we recall that the eigenstates of the target Heisenberg Hamiltonian have well defined values of total spin quantum numbers (S and m). This allows us to restrict the variational search to the appropriate sector of the Hilbert space by starting from a state with well defined S and m and using parametrized unitary operators that preserve these symmetries. A suitable choice is represented by pairwise Heisenberg interaction terms of the form

$$W_{ij}\left(\theta \right)=e^{-i\theta {\mathbf{s}}_i·{\mathbf{s}}_j}.$$

(5)

These are inspired by the real-time quantum evolution induced by the problem Hamiltonian and essentially correspond to the eSWAP operations proposed in Ref. [87], that preserve the values of S and m. Each 2-qubit operation $W_{ij}\left(\theta \right)$ can be decomposed into a sequence of elementary single-qubit gates combined with 3 CNOTs by using the identity [57]

(6)

where $R_{\alpha }\left(\varphi \right)=e^{-i{s}_{\alpha }\varphi }$ denote the usual single-qubit Pauli rotations around the coordinate axes. A physically motivated ansatz (PMA), based on the Heisenberg Hamiltonian and preserving spin symmetries $U_{PMA}\left(\overrightarrow{\theta }\right)={\prod }_{k=1}^K{V}_k\left({\overrightarrow{\theta }}_k\right)$ can be built by alternating $W_{ij}\left(\theta \right)$ operations:

(7)

Notice that, as shown in the circuit above, even and odd bonds, which do not share common qubits, can be operated in parallel. Moreover, we find that nearest neighbours connectivity is already sufficient to guarantee the reconstruction of the ground state, while limiting the amount of free parameters.

Since $U_{PMA}\left(\overrightarrow{\theta }\right)$ preserves S and m, the choice of the symmetry sector can be done at the level of the initial state of the qubits. This should, in fact, respect the same symmetries of the ground state with respect to total spin, while being relatively straightforward to prepare. In particular, in the presence of anti-ferromagnetic nearest-neighbors interactions, the ground state at low field minimizes S [44]. Hence, for even-numbered spin 1/2 chains or rings, this yields an $S=0$ ground state. By increasing the magnetic field, the Zeeman contribution lowers the energy of the $S=1,m=-1$ state, which becomes the ground state after the first crossing. At high field, where the Zeeman contribution to the total energy dominates, the ground state becomes ferromagnetic $|\downarrow \downarrow \downarrow \downarrow \rangle$. As an example, since the ground state for the Heisenberg square at low $B/J$ has $S=0$, $m=0$, a possible initial state is represented by local singlet pairs [87]   

(8)

The parametrized trial wavefunction then becomes $|\psi \left(\overrightarrow{\theta }\right)\rangle =U_{PMA}\left(\overrightarrow{\theta }\right)|{\psi }_0\rangle$. Similarly, a $S=1$, $m=-1$ initial state of the form

(9)

can be used at intermediate $B/J$ values. Notice that both proposed initial states break the spatial symmetry associated with the target systems, due to the uniform strength of the Heisenberg and Zeeman interaction (with periodic boundary conditions). Although such symmetry is recovered sufficiently well through $U_{PMA}\left(\overrightarrow{\theta }\right)$, which is also not translationally invariant due to the independence of the parameters associated to each 2-qubit bond, the symmetry breaking effect induced by the initial state may, in general, worsen the performances of the ansatz, or require a larger number of repetitions K to (approximately) restore the correct structure. In our case, we observe that the most effective strategy is to apply the parametrized $W_{ij}\left(\theta \right)$, first on the bonds which connect the local singlet or triplet (at higher $B/J$) states, as this balances the structure of the initial state and enables a faster recovery of the correct symmetry properties. As an example, the complete trial wavefunction for $K=1$ and $S=m=0$ designed according to this expedient can be expressed as (notice the reversed order of the $W_{ij}$ links compared to Equation (7))

(10)

As shown in Figure 1, the PMA ansatz achieves very good approximations of the exact ground state. For both the low and medium $B/J$ region, we use $K=1$, resulting in only 4 variational parameters and a total of 12 CNOT gates (which can always be operated in parallel on pairs of non-overlapping bonds). Notice that, although the circuit complexity is comparable or slightly higher in terms of the number of quantum gates with respect to $U_{HA}\left(\overrightarrow{\theta }\right)$, the classical optimization converges much faster compared to the heuristic case, as shown in Figure 1b.

In Figure 1c,d, we extend both the HA and PMA strategies to a larger Heisenberg ring with $N=6$ spins. Here, due to the increased complexity of the system, $L=5\left(4\right)$ layers are required for the Heuristic Ansatz, resulting in 24(20) classical parameters in the low and intermediate field regions, respectively. The circuit depth is also significantly increased, including 15(12) CNOT gates not parallelizable. The Heisenberg-inspired ansatz can also be generalized in a straightforward way by adding two additional nearest neighbours $W_{ij}\left(\theta \right)$ operations per layer compared to the $N=4$ case. Similarly, the initial state $|{\psi }_0\rangle$ can be constructed with a combination of three local singlet pairs (low $B/J$) or of two singlets and one $S=1,m=-1$ state (after the first level crossing). Ideal performances for the PMA are obtained with $K=2$, namely 8 free parameters and 36 (parallelizable) CNOT operations. Notice that, consistently with the $N=4$ ring, the minimal PMA required to achieve an exact reconstruction of the ground state has a slightly more complex circuit compared to the HA providing comparable accuracy, but can be optimized in a significantly lower number of iterations with respect to the latter. In summary, as we will also see below in the presence of simulated hardware noise, there exists a trade-off between shallower heuristic circuits with more parameters and more involved but physically motivated ones, which must be carefully balanced when assessing the overall computational efficiency and performance of the VQE applications under study.

Effect of Noise

After assessing the performances of the VQE algorithm in the ideal scenario where no errors affect the execution of the quantum circuits, a crucial question concerns the robustness of these solutions to hardware noise. This is particularly relevant in the near term perspective, where small to medium-sized noisy quantum processors will become available. The most important errors expected on present day real devices and included in our simulations are the following:

• finite relaxation ($T_1$) and coherence ($T_2$) times of the physical qubits, which typically lead to amplitude and phase damping effects;
• single- and 2-qubit gate errors (the latter being usually much higher), acting during the implementation of each quantum operation, due to both imperfections of the coherent qubit manipulation and additional incoherent effects (e.g., depolarizing Pauli noise);
• readout errors, associated with imperfect measurements and erroneous assignment of the outcome, which can be modeled, for example, as bit flip channels.

It is worth pointing out explicitly that all the above noise processes lead to an overall non-unitary dynamics of the quantum states under study, thus requiring a density matrix formalism for an adequate description (see Materials and Methods, Section 5). In practice, noise influences the VQE execution by affecting the precise estimation of the required observables, most importantly, the variational energy ${\langle \mathcal{H}\rangle }_{\theta }$ which serves as cost function for the classical optimization. As a result, on top of lowering the overall accuracy for the reconstruction of the ground state properties (even in the case of exact preparation of the ground state wavefuction), noise can in principle affect the solution of the minimization problem itself, e.g., by significantly increasing the number of iterations required for convergence. In fact, although variational and adaptive approaches, such as the VQE algorithm, may show some resilience to moderate levels of hardware errors [97], it has also been theoretically demonstrated that noise-induced barren plateaus [98] arise in general for deep variational ansatzes.

In Figure 2, we report the solution of the 4-spin Heisenberg square in the presence of typical levels of hardware noise. These are simulated via the Qiskit qasm_simulator [99], making use of the NoiseModel.from_backend function (see Materials and Methods, Section 5). The noise parameters are derived from the calibration data of the ibmq_kolkata 27-qubit quantum processor, for which a quantum volume of 128 is reported. For simplicity, we construct a uniform noise model by using average parameters across the qubit register ($T_1\sim 135$$\mathsf{\mu }$s, $T_2\sim 125$$\mathsf{\mu }$s, single-qubit gate error $\sim 2.5\times {10}^{-4}$, 2-qubit gate error $\sim 8\times {10}^{-3}$, readout error $\sim {10}^{-2}$). In all noisy simulations, we use the SPSA classical optimizer [100], which, due to its stochastic nature, is particularly effective in presence of fluctuations and errors affecting the evaluation of the variational energy. Moreover, we increase the stability and accuracy of all energy evaluations and observations made on the noisy quantum circuits by applying the measurement error mitigation protocol provided in the Qiskit Ignis framework [99]. This is based on a calibration matrix, which is first constructed via preparation and measurement experiments on computational basis states and later used as a filter on noisy measurement outcomes [66].

As can be seen, although some quantitative inaccuracies arise (together with an increase in the number of optimization iterations for convergence, Figure 2c), the energy profiles, the fidelity and ground state observables are well recovered by both ansatzes. We consider, in particular, the expectation values of the magnetization $M={\sum }_i{s}_i^z$, of the total spin $S^2$ and of local one- and two-body spin operators $s_i^z$ and $s_i^z{s}_{i+1}^z$. Results are shown in panels (b) and (d). The average fidelity to the exact solution, computed from the simulated density matrices before measurement, is around 0.9. Both HA and PMA reconstruct the correct expectation value of different global (b) and local (d) observables, thus demonstrating that the correct symmetries of the target model (i.e., rotational and translational invariance) are recovered. It is worth noting that PMA, thanks to the much smaller number of optimization parameters, converges much faster than HA, taking one order of magnitude less iterations (see panel c).

A similar set of results is shown in Figure 3 for the $N=6$ ring: here, the degradation of performances in terms of fidelity and energy values suggests that the complexity of the required circuits (for both the HA and PMA cases) is closer to the limits of current technology. Nevertheless, it is immediately clear that the PMA still requires a significantly smaller number of iterations to converge to a similar accuracy, given the much smaller number of optimization parameters. Of course, the reduction in the number of parameters of the ansatz is often associated to more complex quantum circuits required to build the (highly entangled) trial state. In a near term perspective, where only noisy quantum processors are available, this makes PMA advantageous if the noise level is not too high. In fact, similar limitations in terms of sensitivity to noise (affecting also circuit trainability) may apply to other well known application-motivated ansatzes which require polynomial-depth quantum circuits, such as some instances of the Quantum Alternating Operator Ansatz (QAOA, commonly employed for combinatorial optimization problems) or the unitary coupled cluster (UCC) ansatz for quantum chemistry (see e.g., [98]). In addition, while (consistently with the $N=4$ case) the HA provides slightly better estimates for the system energy, the total magnetization and the total spin of the ground state wavefunction are much better reproduced by the PMA. This highlights the benefit of imposing native symmetry constraints in the ansatz even when these are not preserved, in general, by the effect of noise. In fact, symmetry-based error mitigation techniques could be specifically designed to further improve the performances [63,65,70,101]. The symmetry breaking effect induced by hardware noise is evident in the intermediate field region by comparing the shape of the radar charts in Figure 3d obtained using HA and PMA ansatzes: only the latter preserves the star-shape shown by the exact results (black dots).

Overall, the noisy simulations indicate that, for the systems sizes of interest, the execution of VQE on quantum processors is still quite demanding in terms of the number of measurements and circuit depths. In addition, the simulations reported above were performed assuming at least circular connectivity between the qubits. When using real hardware, one has to face the limited (often linear, nearest-neighbor) connectivity of the chip, which requires one to introduce additional SWAP operations to implement two-qubit operations between distant qubits, thus significantly increasing the circuit depth. Nevertheless, state-of-the-art devices are quickly approaching the required levels of accuracy, making these calculations possible in the near future. In addition to the readout error mitigation strategy adopted in this study and the symmetry-based tools mentioned above, other techniques could also be applied to enhance the quality of the results, and may become a critical resource to enable real hardware experiments. Examples include zero-noise extrapolation methods [64] and systematic procedures designed to increase the effective overlap of (generally mixed) variational quantum states with the true Hamiltonian ground state [67]. All these protocols could ultimately give access to more faithful reconstructions of all ground state observable properties, like energy and magnetization.

2.3. Finite-Size and Parity Effects

Having demonstrated the improvement obtained by using PMA versus HA, we now use the PMA to study the effect on local observables of changing the topology of the spin structure. In particular, we compare expectation values of local spin operators, $\langle s_i^z\rangle$, on the ground state of open vs close and odd vs. even-numbered rings, computed by VQE. These quantities give access to finite-size and parity effects, already evidenced in experiments on molecular nanomagnets [45].

Results obtained by VQE simulations are reported in Figure 4, with symbols of different colors corresponding to a different number of sites and closed/open topology, compared with exact values (indicated by dashed lines). Closed rings display translational invariance symmetry, yielding uniform values of $\langle s_i^z\rangle \forall i$. At low field, the ground state is known to be $S=0$, leading to $\langle s_i^z\rangle =0$, while at intermediate field $M=-1$, and hence $\langle s_i^z\rangle =-1/N$. These values are perfectly reproduced at low field (panel a), while at intermediate field (b), the effect of noise and of the finite number of shots used to evaluate energy in VQE iterations slightly remove the symmetry of the $N=6$ case, which is nevertheless still substantially preserved.

Open chains are instead symmetrical with respect to the central site (if odd-numbered) or bond (if even). At low field [Figure 4a], even chains still display an $S=0$ ground state and uniform $\langle s_i^z\rangle =0$, while the odd $N=5$ chain is characterized by an alternating pattern of positive and negative $\langle s_i^z\rangle$, satisfying $\langle M\rangle =\langle {\sum }_i{s}_i^z\rangle =-1/2$. Both cases perfectly match our noisy VQE simulations. At intermediate field [Figure 4b], the mirror symmetry of $\langle s_i^z\rangle$ is apparent in our simulations, again in very good agreement with expected results, for both even- and odd-numbered open chains. In general, VQE simulations reproduce the expected behavior well, with results only slightly worsened by increasing N. This demonstrates that the PMA inspired by the Heisenberg Hamiltonian is able to correctly predict the structure and the symmetries of the system ground state.

3. Dynamical Correlation Functions

Having outlined a general strategy to find the system ground state using VQE, we now move to compute dynamical properties. A quantity of fundamental interest in this context is represented by dynamical spin-spin correlation functions, which are key ingredients for computing several observables, such as the magnetic inelastic neutron cross-section [89]. They are defined as follows:

$${\mathcal{C}}_{ij}^{\alpha \beta }\left(t\right)={\langle s_i^{\alpha }\left(t\right)s_j^{\beta }\rangle }_0=\sum _p\langle 0|s_i^{\alpha }|p\rangle \langle p|s_j^{\beta }|0\rangle e^{-i{E}_{p}t},$$

(11)

where $|p\rangle$ are eigenstates of the system with energy $E_p$. Their direct calculation constitutes a hard task for a classical computer. Here, we apply the scheme introduced in Refs. [57,89] for open spin chains to the case of a closed ring and we compute ${\mathcal{C}}_{ij}^{\alpha \beta }\left(t\right)$ using a quantum circuit of the form

(12)

where we consider an ancilla qubit a, in addition to the N qubits used to simulate the dynamics of the target system. To compute ${\mathcal{C}}_{ij}^{\alpha \beta }$ correlation functions, the ancilla is initially prepared in a $\left(\right|0\rangle +|1\rangle )/\sqrt{2}$ superposition via a Hadamard gate, then it is entangled with qubit j of the target system by a controlled-$\beta$ gate in which qubit j is the target. In the following step, we perform a quantum simulation on the target qubits, i.e., we implement a sequence of gates reproducing the time evolution operator $U\left(t\right)=e^{-i\mathcal{H}t}$. Finally, a controlled-$\alpha$ gate between a (control) and i (target), conditioned by the ancilla being in 0, precedes the measurement of $s_x$ or $s_y$, which are used to obtain the real and imaginary parts of ${\mathcal{C}}_{ij}^{\alpha \beta }$, respectively. Measurement in a basis different from ${\sigma }_z$ is achieved by implementing a proper rotation before the standard basis measurement, i.e., a H gate to measure $\langle s_x\rangle$ and a $R_x(\pi /2)$ rotation to measure $\langle s_y\rangle$.

The unitary operator $U\left(t\right)$ describing the system time evolution can in general be obtained by decomposing $U\left(t\right)=e^{-i\mathcal{H}t}$ into the product of elementary one- and two-qubit gates that we are able to implement on the hardware. Differently from the VQE calculations considered in the previous sections, this represents a direct implementation of a non-variational digital quantum simulation protocol [57]. We consider, for instance, a $N=4$ closed Heisenberg ring in the large-field region, where the Zeeman overcomes the exchange interaction. Here, the ground state is known to be $|\downarrow \downarrow \downarrow \downarrow \rangle$ and the time evolution operator can be decomposed as follows:

(13)

Notice the close analogy with the PMA, which features the same building-blocks. It is worth reminding that, although this is not required in the particular case presented here, computing the time evolution operator $U\left(t\right)$ would require, in general, the Suzuki–Trotter decomposition [57]. On a side note, even though we adopted in this work a digital approach for the simulation of quantum dynamics, we recall that variational methods have also been proposed for carrying out the same task. As an example, this can be done by leveraging the so called Dirac-Frenkel, McLachlan or time-dependent variational principles [102,103]. These methods could become a useful resource to scale up the size of the treatable systems in the near term.

As shown in Ref. [89], correlations functions computed on the real hardware are usually subject to phase and scale errors, compared to theoretical values. Independently from their origin, these errors can be mitigated by taking into account general properties of the examined observables. In particular, we exploit the fact that ${\mathcal{C}}_{ii}^{\alpha \alpha }\left(0\right)$ must be real to fix the phase between real and imaginary parts of ${\mathcal{C}}_{ii}^{\alpha \alpha }\left(t\right)$ and the sum rule ${\mathcal{C}}_{ii}^{xx}\left(0\right)+{\mathcal{C}}_{ii}^{yy}\left(0\right)+{\mathcal{C}}_{ii}^{zz}\left(0\right)=s_i(s_i+1)=3/4$ (or, equivalently for the isotropic Heisenberg model under study, ${\mathcal{C}}_{ii}^{\alpha \alpha }\left(0\right)=1/4$) in order to fix the scale factor. We then apply the same phase and scale (PaS) correction at all times. In the case of cross-correlations (${\mathcal{C}}_{12}$), we use average phase and scale values obtained from the analysis of ${\mathcal{C}}_{11}$ and ${\mathcal{C}}_{22}$.

Raw experimental results obtained from implementing the circuit (12) on the ibmq_bogota IBM Quantum processor are shown in Figure 5 (light colors), for both real and imaginary parts, compared with exact values (continuous lines). After applying the PaS correction (dark colors), the agreement with experiments on the real quantum hardware using 5 qubits is very good. Notice, in particular, how the error mitigation strategy helps to recover the correct amplitude and phase of the oscillations for ${\mathcal{C}}_{12}$ cross-correlations.

4. Discussion and Conclusions

In summary, we have shown how quantum computers can be used to efficiently simulate the static and dynamical properties of finite-size spin systems, such as molecular nanomagnets. We have demonstrated this by performing simulations, including a realistic model of noisy on typical devices, and experiments on real noisy quantum processors.

In particular, we have used the hybrid (quantum-classical) variational quantum eigensolver algorithm to find the ground state of Heisenberg spin chains and rings of variable length. For a successful determination of ground state observables on near-term noisy devices, we have highlighted the importance of restricting both the variational search and the circuit size by exploiting the spin symmetries of the investigated systems.

We have then computed dynamical spin-spin correlation functions on a closed Heisenberg square, extending the approach introduced in Ref. [89]. In this respect, we have demonstrated that dynamical properties can be reliably obtained by applying error-mitigation strategies (based on general properties of the computed observables) to the raw output data.

The scheme followed in this work can be extended to investigate anisotropic systems, characterized for instance by anisotropic or anti-symmetric exchange couplings, which could yield anti-crossings between $|S,m\rangle$ states in the energy-level diagram. Far from these crossings, S and m are still practically conserved, while in the vicinity of the anti-crossings (where anisotropic terms of the target Hamiltonian are expected to mix different S multiplets), one should generate a superposition of S states to create a variational ansatz or to investigate the dynamics.

Further improvements in the capabilities of quantum hardware will enable us, in the near future, to increase the complexity of simulations. For instance, although we have focused here on spin 1/2 systems, qubit-based architectures can also be used to simulate spin systems consisting of interacting $s_i>1/2$ ions [20,104,105], by using more than one qubit to encode each spin. In addition, more complex molecular spin structures, with increased connectivity between the spins, could be investigated, e.g., in a regime of competing spin-spin interactions [39]. Finally, the quantum simulation of dynamical properties could be appended to VQE in order to compute the system evolution starting from non-trivial ground states.

All these perspective applications will require a larger number of qubits and larger circuit depth. Hence, error mitigation strategies will still be required to improve the quality of the results on noisy devices. In the case of Heisenberg chains, one could, for instance, exploit other symmetries besides the total spin, such as those related to intermediate spin quantum numbers.

5. Materials and Methods

Simulations

All numerical simulations and circuit executions on IBM Quantum superconducting processors were performed using the Qiskit python library [99]. For noiseless simulations, the statevector_simulator backend was employed, thus solving for the exact output states and observables via linear algebraic manipulations. Noisy simulations were instead executed with the qasm_simulator backend, which includes the effect of statistical sampling on the observed results, mimicking the real measurement process. In all calculations, each circuit was repeated $n_{shot}=8192$ times to reconstruct output averages. The hardware noise was modelled with the NoiseModel.from_backend method [99], which can make use of calibration data from real IBM Quantum devices, implementing thermal relaxation and decoherence effects based on the real physical parameters ($T_1$, $T_2$, frequency and temperature) found on the processor. Gate errors are associated also with local depolarizing channels, whose strength is tuned in such a way that, when combined with thermal noise and taking into account the gate durations, the effective total noise matches the hardware gate infidelities reported in the calibration data. Moreover, readout errors are included with local bit flip channels. We also notice that fidelity tests were performed, in the case of noisy simulations, by  operating the qasm_simulator in density matrix mode, thus extracting the simulated noise-affected quantum state (before measurement) to be compared, e.g., with ideal ground states obtained from exact diagonalization.

Two different optimization algorithms have been used: COBYLA (Constrained Optimization By Linear Approximation) for noiseless simulations with statevector_simulator backend and SPSA (Simultaneous Perturbation Stochastic Approximation algorithm) for noisy simulations. The latter (which is a stochastic gradient-free optimizer) works better in the presence of noise [106]. Executions on real quantum hardware were performed on IBM Quantum devices, made available online for remote access via quantum-computing.ibm.com. These processors are based on superconducting technology, with fixed-frequency transmon qubits, connected to superconducting resonator waveguides for control and readout via microwave pulses [59,107]. The native 2-qubit entangling CNOT operation is obtained via microwave-activated cross resonance interactions [108], mediated by dedicated waveguides connecting neighbouring qubits in linear chains or (portions of) heavy hexagonal lattices.

To accelerate the simulation study and to perform all the calculations, we developed and relied on a general framework to exploit all the Qiskit functionalities. This framework can be used in many different formats and for different use cases. Starting from a general configuration file in which all the required variables and parameters are set, we were able to leverage on several HPC resources in parallel to perform calculations and post processing of data. The developed framework can be found at https://automatic-computation-framework.readthedocs.io/en/latest/.