Molecular spin qudits for quantum simulation of light-matter interactions

We show that molecular spin qudits provide an ideal platform to simulate the quantum dynamics of photon fields strongly interacting with matter. The basic unit of the proposed molecular quantum simulator can be realized by a simple dimer of a spin 1/2 and a spin $S$ transition metal ion, solely controlled by microwave pulses. The spin $S$ ion is exploited to encode the photon field in a flexible architecture, which enables the digital simulation of a wide range of spin-boson models much more efficiently than by using a multi-qubit register. The effectiveness of our proposal is demonstrated by numerical simulations using realistic molecular parameters, whose prerequisites delineating possible chemical approaches are also discussed.


I. INTRODUCTION
In the last few years, quantum computers have emerged as a disruptive technology that promises to solve a large class of problems much more efficiently than any classical machine. The first noisy quantum processors [1] are already available and enable the implementation of non-trivial algorithms targeted to specific tasks [2][3][4][5]. In particular, thanks to their intrinsically quantum logic [6], they could be used already in the short term to simulate the dynamics of classically intractable quantum systems. Understanding the behavior of matter at the nano-scale is a fundamental step to design new molecules, materials and devices. However, the "wonderful problem" of quantum simulation "doesn't look so easy" [7]. In fact, many examples of interest for Physics and Chemistry, such as atoms interacting with light or with thermal baths, are intrinsically difficult to be modeled on current qubit-based architectures [8].
In this respect, Chemistry offers the change of perspective which could overcome the aforementioned difficulties. Indeed, molecular spin systems characterized by a sizeable number of accessible levels can be used to encode multi-level logical units (qudits). Each molecular qudit could replace several distinct qubits in various algorithms [9], thus greatly simplifying manipulations of the register. Magnetic molecules are the ideal playground to implement this alternative architecture. Indeed, they are characterized by long coherence times [10][11][12][13][14][15][16], which can be even enhanced by chemically designing the molecular structure [17,18] or targeting protected transitions [19,20]. Moreover, the spin state of these systems can be easily manipulated by microwave or radio-frequency pulses [21], thus implementing single-and two-qubit gates in permanently coupled [22][23][24] or scalable architectures, for which different ways of switching the qubit-qubit interaction were proposed [25][26][27][28]. Recently, the idea of exploiting the additional levels typical of these systems for quantum error correction has been put forward [21,29,30].
In these works, a multi-level molecule is used to encode a protected qubit within a single object, in place of the many qubits required by standard block-codes [31].
Here we show how the qudit nature of magnetic molecules could simplify the practical implementation of important quantum simulation algorithms. We focus, in particular, on the simulation of light-matter interaction processes in the ultra-strong coupling regime, a problem that does not generally allow for a perturbative treatment and is therefore hard to be solved on a classical computer. This class of models, which are of crucial importance for many fundamental investigations ranging from cavity quantum electrodynamics to photochemistry [32,33], have mostly been tackled so far with analog [34][35][36] or digital-analog [37,38] simulation strategies. The necessity to describe radiation modes (characterized in principle by infinitely many degrees of freedom) represents a major challenge for digital approaches. Standard encodings, designed for multiqubit architectures [39], either employ an exponentially large Hilbert space (using a number of qubits equal to the number of simulated photons [40]) or reduce the number of qubits at the price of non-local qubit-qubit interactions and hence complex quantum circuits [39,41].
Conversely, here we reduce both the hardware overhead and the complexity of manipulations by mapping each photon mode to a single spin S qudit. Thanks to the power of coordination chemistry, different qudits can be linked together and, e.g., to spin 1/2 units encoding two-level atoms [42][43][44][45], in non-trivial molecular geometries. This, together with the capability of manipulating the state of the hardware by resonant and semi/resonant microwave pulses, would allow us to digitally simulate the atom-photon dynamics involving multi-mode fields and/or multiple atoms [40,46]. In particular, we show that very simple molecules consisting of dimers of transition metal ions (a spin 1/2 and a spin S ≥ 3/2) can be used to efficiently simulate atom-photon interactions in a non-trivial range of parameters up to strong and ultra-strong coupling [32,33]. The same approach can be extended to simulate, e.g., lattice gauge models involving many field excitations [41], by exploiting the remarkable capabilities of coordination chemistry in synthesizing multi-center molecules with very large total spin [44]. In the following, we design the sequence of pulses allowing us (i) to determine the ground state of the simulated system using the variational quantum eigensolver algorithm (VQE) [3,[47][48][49] and (ii) to follow the time evolution of the system prepared in an out-ofequilibrium initial state. The remarkable performance of the proposed hardware is demonstrated by numerical simulations with parameters corresponding to existing coordination compounds [12,20,[50][51][52][53][54], including the effect of decoherence and the full sequence of pulses needed to implement the algorithms. These results make the here-proposed molecular quantum simulator very promising and pave the way to forthcoming proof-ofprinciple experiments. We finally note that such qudit encoding can be easily extended to any other boson field, thus allowing one to simulate along the same lines many other important models, ranging from phonon vibrations [55][56][57], possibly interacting with spins [58], to lattice gauge theories [41] and complex quantum optical setups [59].

A. Molecular Quantum Simulator
The proposed molecular hardware for quantum simulation is sketched in Fig. 1 (left part). It is a dimer consisting of a spin S 1 ≥ 3/2 qudit that we exploit to encode the boson field, and an effective s 2 = 1/2, described by the following Hamiltonian: (1) Here, the first two terms represent the Zeeman interaction of the two spins with an external magnetic field B applied along z axis and µ B is the Bohr magneton. The third term is the zero-field splitting on the qudit (important to make all qudits transitions spectroscopically distinguishable) and the last one models an exchange or dipolar interaction between the two ions. To reduce our assumptions, we consider in the following axially anisotropic (J z = −2J x,y ) coupling, modeling a dipole-dipole interaction between the two centers. Different forms of the spin-spin interaction or of the singleion anisotropy do not hinder the implementation of our scheme. The only requirement concerns the hierarchy of interactions: the transverse component of the spin-spin coupling must be much smaller than the difference be-tween the excitation energies of the two spins S 1 and s 2 . This condition guarantees that the eigenstates of Hamiltonian 1 are practically factorized products of the states of the two spins, and can thus be labeled by S z1 and s z2 eigenvalues: |ψ m1m2 ≈ |m 1 |m 2 . These requirements are easily fulfilled in coordination compounds containing a spin 1/2 ion coupled to a spin S 1 transition metal ion. The latter provides the ideal qudit for the proposed architecture. As shown below, the relatively small number of levels of these qudits (d = 2S 1 + 1 ≤ 6) is already sufficient to simulate lightmatter interaction from strong to ultra-strong coupling regimes. In addition, transition metal ion complexes with quenched orbital angular momentum ensure significantly long coherence times [12,17,20,50,54], important to achieve a good simulation. We consider, in particular, two paradigmatic cases: Cr III and Fe III ions in distorted octahedral environment, yielding 3d 3 and 3d 5 electronic configurations with a single electron per orbital and thus S = 3/2 and 5/2, respectively [52,53]. Due to the practically complete quenching of the orbital angular momentum, the spectroscopic tensor g is isotropic and close to the free electron value, while single ion anisotropy is typically in the ∼ 0.2 − 0.3 cm −1 range [52,53]. As an illustrative example, in the simulations reported below we use D = 0.24 − 0.30 cm −1 and g = 1.98 for Cr III , as in [17,20] and D = −0.30 cm −1 and g = 2.00 for Fe III , as reported e.g. in Ref. [50,60]. [61] These single-ion qudits can be weakly coupled through bond or through space to a spin 1/2 ion, such as Cu II in distorted octahedral ligand cage [12,52,53], typically characterized by g ∼ 2.1 − 2.3 and in some cases also by remarkable coherence times [12]. In the following we assume g 2z = 2.3, significantly different from g 1z = 1.98−2 to ensure factorization of the system wave-function. For the dipolar interaction we assume J x,y = 0.008 cm −1 , which corresponds to a dipolar coupling (in the point dipole approximation) between ions at a distance ∼ 6Å. A more extensive discussion on possible physical implementations is provided in Sec. III. These parameters, combined with a static field of ∼ 0.3 − 0.5 T ensure that ∆m 1,2 = ±1 transitions needed to manipulate the state of the system fall within the 20 GHz range typically explored in coplanar microwave resonators [62,63].
Having described in detail the molecular hardware, we now switch our attention to the target model, object of our simulation, and on how to map it onto the hardware. The target Hamiltonian is the Rabi model [64][65][66]: Here σ z and σ x are the usual spin 1/2 operators, while a † (a) are bosonic creation (annihilation) operators, [a, a † ] = 1, G is the atom-photon coupling, Ω (ω a ) is the photon (atom) excitation energy, and we have assumed = 1. Hamiltonian (Eq. 2) describes the interaction between a radiation field and a two level system, such as an atom or a spin 1/2 particle. It has recently attracted a great interest in the context of quantum computing, with efforts devoted to achieve the strong coupling between superconducting qubits or spin systems and quantized photons within wave-guide resonators [67,68]. Behind its apparent simplicity, our target model can reveal interesting Physics and non-trivial behaviors associated to ultra or deep strong coupling regimes [33], in which light and matter strongly mix together and exchange excitations without conserving energy [32]. Such a regime can also give insights into fundamental principles of lattice gauge theories [69]. We fix in the following ω a = Ω/2 and study the model for increasing values of the G/Ω ratio, the threshold for the ultra-strong coupling regime usually being G/Ω 0.25. The molecular processor described by Hamiltonian 1 can be used to compute ground state properties and to mimic the dynamics of the target Hamiltonian 2. To achieve this, we first need to encode the boson field into the spin qudit. Notice that a very good approximation can be obtained by truncating the boson field to a relatively small number of levels. Hence, the d = 2S 1 + 1 levels of the qudit are sufficient to encode the radiation field with negligible error, by truncating the radiation field to a maximum number n M = 2S 1 of bosons. The mapping between S z1 eigenvalues and number of bosons (n = a † a) is shown in the bottom part of Fig. 1. In parallel, the two-state atom appearing in the Rabi Hamiltonian can be directly encoded on the hardware spin 1/2 degrees of freedom.
Complete control of the hardware is achieved via microwave pulses resonant (or semi-resonant) with specific excitations of the spin 1/2 or of the qudit. In particular, ∆m 2 = ±1 transitions allow us to rotate the state of the qubit, while ∆m 1 = ±1 pulses are used to excite the qudit. Moreover, the spin-qudit interaction enables conditioned (entangling) operations.

B. Variational quantum eigensolver
The starting point to derive many important properties of the examined systems is the determination of its ground state wave-function. This task can be achieved using the variational quantum eigensolver (VQE) approach [3,[47][48][49]. This is a hybrid quantum-classical algorithm, particularly resilient to noise and therefore well suited for near term quantum processors. It exploits the fact that the energy expectation value is minimum for the ground state of the system. The quantum hardware is used to generate an approximation of the ground state (also known as trial wavefunction or variational ansatz) for the target model, which depends on a set of free parameters θ i , and to evaluate the energy expectation value. Minimization of the evaluated energy by a classical subroutine allows us to explore the parameter space until convergence to the system ground state. It is worth noticing that this method is typically much less demanding, compared to the digital simulation of real time evolution, in terms of the complexity and length of the required sequences of quantum operations to be implemented.
Here we demonstrate an implementation of the VQE on the proposed qudit architecture applied to the target Rabi Hamiltonian, Eq. (2). We construct the trial wavefunction by designing some basic quantum operations achieved in practice via external control microwave pulses. In particular, as shown in Fig. 2a, we assume a (S 1 , s 2 ) = (3/2, 1/2) hardware platform and we combine pulses resonant with transitions of s 2 (green arrows), implementing rotations of the qubit, with ∆m 1 = ±1 pulses on the S 1 = 3/2 spin (black). To introduce entanglement in the approximate ground state, the operations on the S 1 = 3/2 spin are actually conditioned by the state of the spin 1/2, i.e. we rotate each pair of qudit levels by ±θ i depending on the sign of m 2 (see Fig. 2a). In total, the ansatz contains only 4 free parameters, and can be implemented with a sequence of microwave pulses that can be as fast as 100−200 ns. We also mention that such variational structure, which can be natively realized on our proposed qudit architecture, is closely related to the socalled polaron ansatz, which was recently implemented on superconducting quantum hardware [40] through nontrivial decompositions into elementary qubit operations.
In this demonstration, we combine a classical optimization routine, namely the Nelder-Mead simplex algorithm [70], with numerical simulations of the unitary transformations corresponding to every choice of the variational parameters. In fact, each sequence of microwave pulses can be seen as the series of quantum operations reported in the inset of Fig. 2b. Here the black thick (green narrow) line represents the qudit (qubit). Conditioned qubit-qudit operations are depicted with black boxes, while single qubit rotations are shown in green, in direct correspondence with Fig. 2a.
Simulations are performed according to a realistic hardware setup, including all the required external control pulses and molecular parameters discussed above. The effect of a finite spin coherence time T 2 is included by simulating the dynamics of the hardware density ma-trix ρ according to the Lindblad master equation [71] dρ dt where time dependent H 1 term in the Hamiltonian indicates the presence of external oscillating control fields.
For simplicity, we assume the same value for the T 2 for both spins in the hardware. In Fig. 2c, we report results of the VQE algorithm simulated by assuming a realistic value of the spin coherence time (10 µs, symbols), compared with exact values (lines) for both the ground state energy and some ground state properties of interest. Notice that, over a wide range of G/Ω values, the proposed ansatz achieves very good approximations of the exact ground state. The limiting factor is essentially the expressibility of the trial wavefunction, i.e. the fact that by using the set of operations reported in Fig. 2a we may not achieve the exact form of the true ground state. This limitation can in principle be overcome by repeating the same basic parametrized structure more than once. It is worth noting that a finite coherence time, similarly to small imperfections in the practical realization of quantum gates, only minimally affects the final results. In fact, consistently with the underlying variational principle, noisy ground state energy estimates sometimes converge to values slightly larger than the exact ones.

C. Digital quantum simulation of strong light-matter interaction
After investigating ground state properties, we now move on to show how the proposed molecular qudit-based processor can be used to simulate the dynamics of the Rabi model. The digital quantum simulation of the target Hamiltonian H S requires to implement the transformation: This can be approximated to the product of unitary terms e −iH S t ≈ (e −iωaσzt/N e −iΩaa † at/N e −i2Gσx(a+a † )t/N ) N by dividing the transformation into small time steps t/N , according to the Suzuki-Trotter decomposition. Each unitary step is then implemented by a sequence of micro-wave pulses. For instance, the effect of the diagonal operator a † a is obtained by pulses semi-resonant with ∆m 1 = ±1 transitions [72], while the term σ x (a + a † ) is simulated by resonant ∆m 1 = ±1 transitions conditioned by the state of the qubit and essentially correspond to the similar ones employed in the VQE above.
In Fig. 3 we show the digital quantum simulation of the Rabi model, Eq. (2), realized with the spin qudit encoding described above and for increasingly challenging choices of the G/Ω ratio. Large G/Ω values introduce peculiar features in the dynamics of the target system: the rotating wave approximation fails and the total number of excitations is not conserved. This non-trivial behavior emerges in our simulations below, where we report the time evolution of the average number of photons n photons in the radiation mode and of the atom population σ z , assuming an initial vacuum state with zero photons and the atom in its ground state. This vacuum state (with no excitations) would not be subject to any evolution for small G/Ω ratios. Hence, oscillations in n photons and σ z are a direct signature of the ultrastrong coupling regime. In all panels, we compare the reference curves, computed via exact matrix exponentiation, with numerical simulations of a realistic hardware obtained again by integrating Eq. (3). A quantitative assessment of the overall quality of the results can be obtained by computing the fidelity F = ψ id |ρ|ψ id between the hardware output ρ and the ideal result ψ id of a digital quantum simulation algorithm realized with the same number of Suzuki-Trotter time steps and the same size of the bosonic Hilbert space. The latter can be obtained by with standard matrix algebra.
In the first example, Fig. 3a-b, we show the results of the quantum simulation of the target Hamiltonian H S with G/Ω = 0.25 realized with N = 4 (for t ≤ 5) and N = 6 (for t > 5) Suzuki-Trotter steps. Here, the hardware setup is composed of a spin S 1 = 3/2, encoding a d = 4 photonic space, and a spin s 2 = 1/2 representing the atomic degrees of freedom. The longest pulse sequence requires 1.7 µs, resulting in large average fidelities: F 0.984 for T 2 = 50 µs and F 0.951 for T 2 = 10 µs.
Increasing values of the target G/Ω ratio, Fig 3c-f, yield larger oscillations in the average number of photons and atom populations. To capture these features we need, on the one hand to increase the number of digital steps (N ), on the other hand to enlarge the bosonic space (n M ). This last step is fundamental to correctly capture the system dynamics at significant G/Ω, as clearly shown in panels (g-h), where we compare the time evolution obtained by truncating the number of photons to 3 or 5, for G/Ω = 0.7. Indeed, by slightly increasing n M , we practically obtain the exact dynamics (continuous line). Given n M = 2S 1 , on the synthetic side, this simply translates in changing the qudit spin from 3/2 to 5/2. Conversely, increasing N (and hence the length of our manipulations) requires larger T 2 or faster pulses, e.g. by engineering the molecular spectrum to better resolve all transitions. In this respect, the large degree of chemical flexibility represents a valuable resource. In particular, it is helpful to replace the s 2 = 1/2 with a spin s 2 = 1 system. A promising candidate ion is for example Ni II , for which coherence times in the regime of microseconds were reported [54]. While only two consecutive levels, e.g. m 2 = 0, 1, are used for the actual encoding of the target model, the presence of an additional zero-field splitting term D s 2 z2 in the hardware Hamiltonian greatly improves the frequency resolution of the relevant transitions, thus allowing for larger operation fidelities with reasonably fast control pulses. For Ni II , D can be in the 0.1 − 1 cm −1 range (in octahedral ligand field) [52,53]. In Fig. 3c-d we report a digital simulation for G/Ω = 0.5, obtained with N = 7 on a (S 1 , s 2 ) = (3/2, 1) model hardware.
Here, the pulse sequences last approximately 0.9 µs  on average, resulting in average fidelities F ≥ 0.92 also for T 2 = 10 µs. Finally, we achieve in Fig. 3e-f a digital simulation well above the ultra-strong coupling threshold (G/Ω = 0.7, N = 8) with a model hardware (S 1 , s 2 ) = (5/2, 1) (i.e. with a bosonic space truncated at d = 6). More demanding pulse sequences are required in this case, with an average duration of ∼ 1.6 µs and average fidelity around F 0.84 for the shortest T 2 .

III. POSSIBLE PHYSICAL IMPLEMENTATIONS
Let us now explore potential realizations of molecular qudits displaying a set of properties consistent with the ones employed in our calculations. In many cases, we refer to chemical building blocks already discussed or characterized in the literature.
To identify a suitable molecular platform, we need to combine requirements on the different units discussed in the previous sections. As already illustrated, a prototypical hardware could consist of a dimer of transition metal ions, respectively with spins S 1 ≥ 3/2 and s 2 ≥ 1/2. In order to ensure factorization of the two-ion wave-function, the two ions should be weakly interacting through space or through bond and characterized by g factors significantly different along a given direction. Single-ion anisotropy on both S 1 and s 2 (if the latter is ≥ 1) could for example help to better resolve different transitions. Such single constraints do not appear so stringent. For instance, Cr III and Cu II have sufficiently different g values g Cr = 1.98, g Cu = 2.10 − 2.3 to allow factorization of the wavefunction. At the same time, the individual spin's resonance frequencies are both accessible in the same resonator. D values in the order of the tenth of cm −1 characterize ions that have half filled valence orbitals, like Mn II , Fe III , or Gd III , as well as half-filled t 2g orbitals in octahedral ligand field, such as Cr III . It must be said that in this case the rhombicity and principal directions of the magnetic anisotropy are difficult to predict and control synthetically, but they are not crucial for the feasibility of our scheme. More demanding is the control of the interaction between the spin qubit and the qudit. Dipolar interactions can be easily computed and controlled. The required strength, ca 0.01 cm −1 , is associated to a distance of about 6Å. Such a relatively short distance inside a molecular architecture is compatible with compact linkers like oxalate, cyanide, azide etc. These bridging ligands are very efficient in transmitting also exchange interactions, and thus unsuitable for single spin addressing. The optimal choice falls on very weak exchange interactions that are expected to be almost ubiquitous when the two spin centers are embedded in the same molecular scaffold. However, such weak interactions (of the same order of inter-molecular ones) have been poorly characterized through standard magnetometry techniques in concentrated solids or ab initio calculations. An elucidating example of the wide range of achievable interactions is the case of condensed Cu II porphyrin complexes. Complexes with conjugated macrocyclic ligands have been attracting increasing interest for the relatively long and robust coherence combined with semiconducting properties and convenient processability [13,[73][74][75]. Electron-electron double resonance has been recently used to investigate the spin-spin interactions in edge-fused coplanar Cu II dimers and in meso-meso singly linked dimers [76]. In the latter, the Cu-porphyrin rings are mutually orthogonal and exchange interaction fully suppressed, significantly smaller than the dipolar interaction, estimated to be 0.0028 cm −1 . On the contrary, the planarity imposed by the triple link between the two units boosts the antiferromagnetic exchange interaction to J = −2.64 cm −1 . The difficulty to predict the actual spin-spin interaction (based on ab-initio calculations or simple geometrical considerations) does not hinder the implementation of our scheme. Indeed, once the complex has been characterized, it is possible to tune the external field in order to ensure factorization of the wave-function. This could require to adapt the experimental setup to work at larger frequencies than commercial resonators, as demonstrated for instance in Ref. [63], where superconducting coplanar resonators operating up to 50 GHz were reported. These superconducting resonators could also employ high-T c superconductors to support large magnetic fields [77]. The choice of the linkers between the two magnetic ions should also fulfil other constraints. In particular, we need to control the decoherence of the system. A coherence time T 2 of 10 µs at low temperatures is often observed for S = 1/2 transition metal ions, especially if the first coordination sphere is nuclear spin free, e.g. oxygen, sulphur, or carbon donor atoms, and if total or partial deutera-tion of the ligand is affordable. This requires to eliminate nitrogen from the first coordination sphere and aliphatic CH n groups in the molecular scaffold, thus reducing the available library of molecular candidates. We should finally remind that an efficient operation of the simulator requires that the qudit-qubit pairs are well isolated, still retaining a control over the molecular orientation. An isostructural diamagnetic matrix is thus mandatory. While this is usually accessible for single qubits, in the case of a two-spins architecture the cocrystallization of the para-and dia-magnetic molecules must occur without metal scrambling. This can be easily avoided using inert d 3 /d 6 ions, as in the case of Cr III and low spin Co III . Metal scrambling is however much more common for labile d 1 /d 9 ions, such as Cu II , requiring the use of polydentate linkers in the design of the molecular architecture.

IV. DISCUSSION AND CONCLUSIONS
Summarizing, we have shown that magnetic molecules are very promising quantum simulators for complex physical systems, in particular for target Hamiltonians involving bosonic variables representing e.g. radiation fields. The many degrees of freedom present in this class of target Hamiltonians make their simulation with a multiqubit register very demanding, both in terms of number of qubits and sequence of operations. In contrast, the multi-level structure typical of magnetic molecules allow us to encode a boson into a single spin qudit, thus greatly simplifying the architecture of the register and its manipulation. The latter can be achieved solely by sequences of microwave pulses, resonant with specific transitions. As an example, we have reported both ground state calculations, performed with the VQE algorithm, and the digital quantum simulation of real time evolution for the Rabi model up to the ultra-strong coupling regime. In all cases, the outcomes obtained by considering realistic hardware parameters are in very good quantitative agreement with exact predictions. These results pave the way to proof-of-principle experiments demonstrating the effectiveness of our proposal. The scheme is flexible and allows one to simulate a wide range of interesting models, thanks to the chemical tunability of the proposed hardware. Indeed, although we have focused here on very simple single ions, much larger S 1 can be obtained by exploiting the total spin ground multiplet of multi-nuclear complexes with tailored interactions [44,78]. With larger S 1 , one could for example include more photons in the simulations, thus enabling the treatment of more exotic regimes such as the deep strong coupling for light-matter interactions [79] or fundamental models such as lattice gauge theories [41,80,81]. The latter require a large number of boson modes and excitations for a detailed description in arbitrary dimension, thus representing a challenging task for both classical devices and near time qubit-based architectures [41]. Ad-ditionally, models involving multiple two-level atoms or boson modes [46] can be simulated by chemically engineering the structures in order to link together several qudit and/or qubits [42,43].
In conclusion, it is worth stressing that a synthetic effort to achieve the conditions highlighted in this work would place molecular nanomagnets among the most promising platforms for the realization of effective quantum simulators.