Quantum simulation of antiferromagnetic Heisenberg chain with gate-defined quantum dots

Quantum-mechanical correlations of interacting fermions result in the emergence of exotic phases. Magnetic phases naturally arise in the Mott-insulator regime of the Fermi-Hubbard model, where charges are localized and the spin degree of freedom remains. In this regime, the occurrence of phenomena such as resonating valence bonds, frustrated magnetism, and spin liquids is predicted. Quantum systems with engineered Hamiltonians can be used as simulators of such spin physics to provide insights beyond the capabilities of analytical methods and classical computers. To be useful, methods for the preparation of intricate many-body spin states and access to relevant observables are required. Here, we show the quantum simulation of magnetism in the Mott-insulator regime with a linear quantum-dot array. We characterize the energy spectrum for a Heisenberg spin chain, from which we can identify when the conditions for homogeneous exchange couplings are met. Next, we study the multispin coherence with global exchange oscillations in both the singlet and triplet subspace of the Heisenberg Hamiltonian. Last, we adiabatically prepare the low-energy global singlet of the homogeneous spin chain and probe it with two-spin singlettriplet measurements on each nearest-neighbor pair and the correlations therein. The methods and control presented here open new opportunities for the simulation of quantum magnetism benefiting from the flexibility in tuning and layout of gate-defined quantum-dot arrays.


I. INTRODUCTION
Analog quantum simulations of magnetism [1] have been performed with a rich variety of experimental platforms, ranging from ultracold atoms in optical lattices [2][3][4][5][6][7] to trapped ions [8,9], scanning tunneling microscopy of atoms on metallic surfaces [10], and superconducting circuits [11]. A recent addition is the use of gate-defined quantum dots as a platform for quantum simulation of the Fermi-Hubbard model. The abilities to independently control the filling of the array, the local electrochemical potentials, and the hopping energy between sites are complemented with methods to probe the charge configuration across an array, spin states, and the electrical susceptibility as well as transport through the system [12]. These already enabled the observation of the transition from Coulomb blockade to collective Coulomb blockade [13] and the observation of Nagaoka ferromagnetism [14], a form of purely itinerant ferromagnetism which occurs at doping with a single hole.
In the Mott-insulator regime, where all sites are occupied by one electron, magnetism is governed by the Heisenberg exchange interaction [15], which favors antiferromagnetic spin alignment. Earlier studies on quantum-dot arrays in this regime demonstrate the sequential control of exchange couplings enabling coherent state transfer [16] and a method to handle cross talk in simultaneous control of exchange couplings [17]. In order to study the many-body properties of this system, novel methods are needed for state preparation in the presence of disorder and temperature as well as for probing spin correlations. Preparation of the Heisenberg ground state is both a useful and exciting goal because of its potential applications such as quantum information transfer [18][19][20][21][22] and quantum simulation of magnetic phases [23][24][25].
In this work, we simulate the antiferromagnetic Heisenberg spin chain in a gate-defined quadruple quantum dot. For this purpose, we develop experimental techniques based on energy spectroscopy and coherent oscillations of the global spin state. These include methods for many-body spin-state preparation and singlet-triplet correlation measurements, which form a powerful probe for the characterization of a many-body spin state [26]. We use these TABLE I. Four-spin shared eigenstates ofŜ 2 andŜ z expressed in a basis of two-spin singlets and triplets on either the left and right or middle and outer pair. States in the rightmost column are the same as states on the same row in the column to the left of it. The four-spin states shown here are, in general, not eigenstates of the Heisenberg Hamiltonian, but the Hamiltonian does operate within a specific ðS; m S Þ subspace.  the total spin operatorŜ 2 , with eigenvalues SðS þ 1Þ, and the spin operator in the z direction,Ŝ z , with eigenvalues m S ∈ ½−S; …; S. The dimensions for these subspaces can be obtained from the Clebsch-Gordan decomposition. For four spins, this results in two global singlet states for which S ¼ 0, nine global triplet states for which S ¼ 1, and five states with S ¼ 2, which form a quintuplet. The triplets are separated into three three-dimensional subspaces and are denoted by T α k , with the subspace magnetization α ∈ ½−; 0; þ and k labeling the energy level, where k ¼ 0 for the state with lowest energy in the respective subspace. Similarly, the singlets are denoted by S k , and the quintuplets by Q β with β ∈ ½−−; −; 0; þ; þþ.
The global spin states can be characterized in terms of the probabilities to measure either two-spin singlets and j indicate the site. The simultaneous eigenstates ofŜ 2 andŜ z for four-spin states can be expressed in this pairwise singlet-triplet basis as shown in Table I. Appendix A discusses the limitations of singlet-triplet measurements to distinguish spin states.
Alternatively, we can characterize the Heisenberg spin chain via its energy spectrum. Based on the symmetries of the Hamiltonian, for four spins the global singlet states form a two-dimensional subspace. For this subspace, the Heisenberg Hamiltonian is with basis states j0 S i and j1 S i from Table I. This subspace has been proposed as a singlet-only exchange-only qubit implementation, which offers increased coherence due to the reduced influence of nuclear spins [29].
For the global triplet states, the three three-dimensional subspaces are identical in terms of energy splittings. The Heisenberg Hamiltonian for each of these triplet subspaces is with basis j0 T þ i, j1 T þ i, and j2 T þ i from Table I. The quintuplet states have zero energy contribution from the Heisenberg Hamiltonian but can be energetically split with an external magnetic field. The energy differences in the subspaces reveal information about the exchange coupling strengths, and characteristic features can be identified. For the singlet subspace, the energy splitting is 1 . It follows that, given J 12 ¼ J 34 , the energy splitting is minimized when J 23 ¼ J 12 ¼ J 34 and, thus, for homogeneous exchange couplings. For the triplet subspace, the energy difference between the two lowest-energy states is minimized if J 12 ¼ J 34 . If J 12 ¼ J 23 ¼ J 34 , the triplet states are equally spaced in energy (see Appendix B for simulated energy diagrams). These characteristic features for the energy spectrum of the Heisenberg Hamiltonian will be experimentally identified, but first we introduce the quantum-dot device and the experimental operation.

III. DEVICE AND EXPERIMENTAL OPERATION
The prototype for the simulation of an antiferromagnetic Heisenberg spin chain consists of a quadruple dot and a sensing dot, which are formed in a device nominally identical to that shown in Fig. 1(a) (see Appendix C 1 for details). The device is based on a GaAs=AlGaAs heterostructure, since this is the technology in which we are able to fabricate high-quality and well-controlled quantum-dot arrays. The exchange couplings are induced by electron wave function overlap, which we here control by detuning the potentials of neighboring dots, such that one electron shifts toward the other [30] (we note that independent control of the exchange couplings can also be achieved by adjusting the tunnel couplings [13,17,31]). In order to control the detuning between one pair of dots without affecting the detuning between other pairs, we define the detunings ε ij as 0 B @ ε 12 where ε i is the negative local energy offset for site i and ε ij ¼ 0 at the interdot transition between charge occupations (1111) and (0211), (1201), and (1102) for ε 12 , ε 23 , and ε 34 , respectively. The ε i are independently controlled using virtual plunger gates, which are linear combinations of the voltages applied to the gates P i [13] (in the figures below, we express ε ij in units of mV; in Appendix C 2, we specify the conversion factor between energy and the applied voltage). The dependence of exchange couplings on detunings can be modeled as [32] with t ij the tunnel coupling between dots i and j. In this way, the exchange couplings can be set independently with the detunings, and increasing detuning results in increasing exchange strength. This method can be extended to larger chains based on alternating the detuning directions such as defined in Eq. (5), which prevents unwanted charge transitions. In addition, the amount of change in chemical potentials could, in the future, be further reduced, by leveraging the fact that the effect of detuning on exchange is weaker closer to the charge symmetry point. The experimental sequence used to operate the quantumdot spin chain is schematically depicted in Fig. 1(b) and is described step by step here (charge-stability diagrams and sequence details are provided in Appendix C 1). Initially, the quadruple dot is tuned in either the (0202) or (1102) charge occupation. For the (0202) case, we load local singlets in the second and fourth dots, by allowing tunneling or cotunneling between dots and the reservoirs. For the (1102) configuration, we load a thermal mixture of two-spin states on the left pair, postselect for triplet loading (see Appendix C 4), and load a singlet on dot four. Next, the electrons are separated to obtain (1111) charge occupation. The global spin state remains a product of local spin pairs, because the left-and right-pair exchange remain large compared to the hyperfine field from the nuclear spins, and the middle-pair exchange coupling is kept small. Then, during the manipulation stage, the exchange couplings are diabatically or adiabatically changed by applying gate voltage pulses with variable rise time. In this work, the pulses are always diabatic with respect to anticrossings between states with different magnetization. Subsequently, the four-spin state evolves under the newly set exchange couplings. Finally, the spin pair(s) to be measured is(are) diabatically isolated from the other spins and measured with single-shot singlet-triplet readout based on Pauli spin blockade [33]. Here, either the left and right pair are sequentially read out, while parking the other pair to avoid capacitive cross talk [34], or the middle pair is read out.
In the remainder of this work, we focus on realizing homogeneous exchange couplings throughout the spin chain. In principle, this could be achieved by calibrating the exchange couplings one at a time [17] and extrapolating to the required tuning while compensating for cross talk. Instead, we develop a two-step spectroscopy method from which we identify directly when the condition of homogeneous exchange couplings is met.

IV. ENERGY SPECTROSCOPY
For gate-defined quantum dots, information about the energy-level spectrum can be obtained from the degeneracies between spin states with different magnetization. This so-called spin funnel method has been used extensively in quantum-dot arrays of various lengths [30,35,36]. Here, we use the same underlying principles in a novel method for simultaneously characterizing multiple exchange coupling strengths in the spin chain. In addition, since the system size is small enough to allow classical numerical computation of its energy-level spectrum, we can validate the quantum simulator by comparing the measured energy spectrum to the numerically computed spectrum.
For the energy spectroscopy measurements, we prepare spin singlets on the left and right dot pairs [see Fig. 1(b)] (i.e., we prepare in the low-energy global singlet j0 S i), diabatically pulse the exchange couplings, allow the system to evolve for 100 ns, and read out the left and right pair. The duration of 100 ns is chosen to allow the coherent time evolution kick started by the pulse to largely damp out. This measurement gives access to correlations in the singlettriplet occupations, P ST , P TT , P TS , and P SS , where the left (right) subindex corresponds to the left-(right-) pair outcome. Decreased P SS indicates mixing of the low-energy global singlet with one of the triplet or quintuplet states. Such mixing occurs most manifestly at anticrossings between the low-energy global singlet state and the polarized states with m S ¼ 1, 2, induced by the gradients of the hyperfine field and the spin-orbit interaction. Depending on which of the other probabilities increases, we can infer information on the nature of the polarized state involved in that specific anticrossing.
We now examine and interpret the spectra in detail. For the measurement shown in Fig. 2(a), the left-and right-pair detunings during the manipulation stage are varied in the presence of a 40 mT magnetic field. The middle-pair detuning is kept fixed and such that the middle exchange coupling is small compared to the outer exchange couplings; thus, the low-energy global singlet state remains almost fully Figure 2(b) shows the result of a corresponding numerical simulation, which helps to interpret the data. The detunings for the anticrossing between the low-energy singlet and the T þ , over a large range of detunings. This demonstrates the independent control of J 12 and J 34 with the detunings ε 12 and ε 34 , respectively. For higher ε 12 (ε 34 ), the horizontal (vertical) T þ 0 line bends away toward lower ε 34 (ε 12 ), which is a manifestation of the capacitive coupling between the left and right pair of dots: The singlet energy on one pair is lowered when the charge occupation for a singlet on the other pair becomes more (02)-like [34]. The capacitive coupling is modeled by adding −DJ 12 J 34 to the diagonal matrix element for jS 12 S 34 i ¼ j0 S i in Eq. (3), with D ¼ 0.015 μeV −1 a prefactor for the interaction strength of the singlet dipoles [34,37]. At the left and right detunings for which the anticrossings with T þ 0 and T þ 1 are closest together, the condition with E Z the Zeeman splitting set by the magnetic field. The coupling between the singlet state and the quintuplet states is of second order in the hyperfine gradients; hence, the mixing between them is less efficient. The corresponding lines, such as the blue line in Fig. 2(b), are most visible when the quintuplet state energy is closest to a triplet state energy (see the white arrow), since the triplet states mediate the second-order coupling (see Appendix D).  Figure 2(c) shows the corresponding numerical simulation. As ε 23 is increased, the energy splitting between T þ 0 and T þ 1 increases due to increased middle-pair exchange. We experimentally reach the condition J 12 ¼ J 23 ¼ J 34 ¼ E Z at the middle-pair detuning for which the energy-level spacing between T þ 2 and T þ 1 is equal to that for T þ 1 and T þ 0 , which is indicated by the two equal-length double arrows in the bottom-left panel.

V. GLOBAL COHERENT OSCILLATIONS
Before we further characterize the homogeneously coupled spin chain, we first demonstrate the coherent nature of the coupled four-spin system. Figure 3 shows global coherent oscillations, during which the full four-spin system evolves, along with Fourier transforms of those oscillations. Because of the symmetries of the Heisenberg Hamiltonian, time evolution occurs within the subspaces of fixed total spin and magnetization. Since we initialize in either only local singlets or at most one local triplet, the subspaces here consist of global singlet states or triplet states, respectively. The insets show numerical simulations based on time evolution under a single-band Fermi-Hubbard model [13] without decoherence effects (see Appendix E). We note that the condition of homogeneous exchange couplings can be extracted from the coherent oscillations as well.
To observe global coherent oscillations, the spin chain is again operated as depicted in Fig. 1(b). A magnetic field of 200 mT is applied here and in the subsequent measurements, to suppress leakage to states with different magnetization during the manipulation stage. For the data shown in Figs. 3(a)-3(d), a triplet state is initialized on the left pair and a singlet on the right pair. The detunings are rapidly pulsed such that the exchange couplings diabatically change and the system evolves coherently under the Heisenberg Hamiltonian during the manipulation stage. This evolution results in oscillations in the singlet-triplet probabilities of the left-and right-pair readout. Figure 3(a) shows a chevron pattern from coherent oscillations between the global triplet states for varying differences between the left-and right-pair detuning and fixed middle-pair detuning. The Fourier transform of these oscillations is shown in Fig. 3(b). In line with the discussion in Sec. II, for fixed J 23 the energy difference between the two lowest-energy triplet states, and thus the oscillation frequency, is minimized if J 12 ¼ J 34 . We point out that all three exchanges are activated in this measurement, so all four spins coherently evolve together. Figure 3(c) shows global coherent oscillations in the triplet subspace for varying middle-pair detuning and with the left-and right-pair detuning fixed such that J 12 ¼ J 34 , as obtained from Figs. 3(a) and 3(b). In Fig. 3(d), the Fourier transform of these oscillations is shown. The middle-pair detuning for which J 12 ¼ J 23 ¼ J 34 , indicated by the white dotted line, can be identified from this Fourier transform as the point where the faint vertical line meets the other more visible line. Here, the triplets with identical magnetization are equidistant in energy (J hom = ffiffi ffi 2 p ), thus reaching the condition of a spin chain with homogeneous exchange couplings, as described in Sec. II. In both the experimental and numerical Fourier transform data shown in Figs. 3(b) and 3(d), one frequency component is typically much more visible than the others. This is caused by the fact that the initial state overlaps mostly with just two of the three eigenstates; hence, the energy difference between these two eigenstates dominates the time evolution.
For the measurements in Fig. 3(e), a product of singlets is initialized, and the middle pair is read out. As described in Sec. II, the energy splitting in the singlet subspace, given Figure 3(f) shows the Fourier transform of the oscillations in the singlet subspace. The ε 23 value for the frequency minimum cannot be precisely identified due to the limited frequency resolution, but an approximate identification is in agreement with the value of ε 23 corresponding to homogeneous exchange couplings from Fig. 3(d). Also, the ratio of the observed oscillation frequencies in the triplet and the singlet subspace at this value of ε 23 is consistent with theory.
The observation of global coherent oscillations demonstrates the coherent nature of the four-spin system. The coherence is limited by hyperfine and charge noise, of which the first can be strongly reduced by working with (isotopically purified) silicon or germanium as host materials [38]. The latter can be largely mitigated when simulating spin models at half filling (one electron per site) by operating at a so-called sweet spot [39,40]. Magnetic field gradients, such as due to hyperfine fields, can also induce leakage out of the fixed total spin and magnetization subspace. For the evolution, this can result in damping of the oscillations toward an offset that corresponds to the leakage state(s) [41], but this effect is strongly suppressed when exchange couplings dominate the hyperfine fields. The optimal visibility of the oscillations is, in general, lower than one, because the eigenstates partially overlap with the readout basis. The measured visibility is further lowered by relaxation during readout (which can be accounted for; see Appendix C 6), leakage, and partial adiabaticity of state preparation and transition to the readout configuration, which are further discussed in the next section.

VI. PROBING THE LOW-ENERGY SINGLET
We finally turn to the preparation and characterization of the Heisenberg spin chain with homogeneous exchange couplings. The ground state of the Heisenberg spin chain, in the absence of an external magnetic field, is the low-energy singlet eigenstate. For homogeneous exchange couplings, J ij ¼ J hom , the low-energy singlet eigenstate S 0 , written in the singlet-triplet basis for the left and right pair, is q . Upon measurement in the two-spin singlet-triplet basis for the left and right pair, we thus have a 1 4 ð2 þ ffiffi ffi 3 p Þ ≈ 0.93 singlet-singlet probability and approximately 0.07 triplet-triplet probability. Alternatively, the same global singlet state written in a basis given by the middle and outer pair is which indicates a 50∶50 probability to measure a singlet or a triplet on the middle pair. When quasistatic hyperfine and charge noise is included (see Appendix E), then numerical simulations result in probabilities of P M S ¼0.50, P M T ¼ 0.50, P SS ¼ 0.91, P ST ¼ 0.01, P TS ¼ 0.01, and P TT ¼ 0.07, which indicates that the noise in the device does not form a direct bottleneck for the quantitative characterization of the spin-chain ground state. Quantitative two-spin singlet-triplet characterization of the low-energy singlet eigenstate S 0 is facilitated by state preparation that is adiabatic with respect to the exchange couplings. Starting from singlets on the left and right pair of dots, the detunings are slowly varied using voltage ramps to increase the middle-pair exchange while reducing the leftand right-pair exchange. Ideally, the singlet product state evolves to the instantaneous low-energy singlet eigenstate of the Hamiltonian at the manipulation stage. Conversely, projection in the singlet-triplet basis for readout requires a diabatic transition from the manipulation to the readout stage. Figure 4(a) shows the singlet probability for the middlepair readout, and Fig. 4(b) shows the triplet-triplet probability for the left-and right-pair readout as a function of ramp-in time and evolution time. As expected, we observe oscillations with a decreasing visibility as the ramp-in time increases, as a result of the more adiabatic character of the state preparation. We note that, for the measurements shown in Figs. 4(a)-4(c) during the manipulation stage, the middle exchange J 23 ≈ 175 MHz is made larger than the outer exchanges, J 12 ¼ J 34 ≈ 85 MHz, in order to increase the visibility of the oscillations, so we can best evaluate the adiabaticity of the state preparation. For shorter ramp-in time, the oscillations bend toward longer evolution time, which is due to the evolution during the ramp-in; thus, at the start of the evolution stage, the state has already evolved further for longer ramp-in time. Figure 4(c) shows the numerically simulated overlap with the low-energy singlet eigenstate at the manipulation point as a function of ramp-in time, indicating >92% overlap for 12-36 ns rampin time. For longer ramp-in time, leakage to quintuplet and triplet states starts to dominate.
We next study how to maximize the degree of diabaticity for the readout pulses with respect to the exchange couplings, in order to acquire the singlet-triplet probabilities for the state at the end of the manipulation stage. To increase the diabaticity given the finite rise time of the arbitrary waveform generator, an isolation stage is added between the manipulation stage and the readout stage for the case where we aim to read the middle pair. In the isolation stage, the voltages are pulsed deep into the (1201) charge region, such that the voltage step is steeper, which makes the pulse more diabatic. Figure 4(d) shows the singlet and triplet probability for the middle-pair readout as a function of ε 23;iso , the middlepair detuning (relative to the readout position) for the isolation stage. The state is prepared with t ramp ¼ 25 ns and with homogeneous exchange couplings at the manipulation stage. As discussed previously, for this condition the low-energy singlet eigenstate ideally has equal two-spin singlet and triplet probability on the middle pair. We see in Fig. 4(d) that the measured two-spin singlet and triplet probabilities gradually approach each other as ε 23;iso is made more negative, pulsing deeper in the (1201) charge region, indicating that the isolation becomes more diabatic in exchange couplings. Pushing ε 23;iso even further, the singlet-triplet probabilities for the middle pair pass slightly beyond 50=50%, which we can trace back to an artifact from the digital filter in the arbitrary waveform generator. Specifically, we measure the detailed rising flank of the arbitrary waveform generator, which shows an undershoot just before the rising flank and ringing after the rising flank, and numerically simulate its effect on the measured singlettriplet probabilities (see Appendix E). Based on these findings, we set in what follows the ramp-in time to 26 ns when aiming to adiabatically prepare the low-energy singlet eigenstate. Numerical simulations similar to Fig. 4(c) show that an overlap of 95.5% with the low-energy singlet state is expected for the homogeneous spin chain. For readout of the middle pair, the middle-pair isolation detuning is set to ε 23;iso ¼ 10 mV. Figure 5 shows the singlet-triplet probabilities for the spin chain as a function of middle-pair detuning (with fixed J 12 ¼ J 34 ≈ 85 MHz) during the manipulation stage. The exchange couplings are calibrated to be homogeneous at ε 23 ¼ −1.8 mV based on measurements of coherent oscillations similar to Fig. 3. The exchange coupling favors twospin singlet formation on every pair, but, due to the monogamy of entanglement, there cannot simultaneously exist such singlets on overlapping pairs. Away from homogeneous exchange couplings, we see in Fig. 5(a) that, for increasing ε 23 and, thus, increasing J 23 , the singlet probability on the middle pair increases, as expected. Conversely, in Fig. 5(b), we observe that as ε 23 decreases (J 23 decreases), P SS (P TT ) increases (decreases), because the singlets on the left and right pair become energetically increasingly favorable compared to the middle-pair singlet. For higher ε 23 , also P TS and P ST increase, which is caused by leakage out of the singlet subspace, due to decreasing energy splitting between the global singlet states and other m S ¼ 0 states.
The experimentally measured numbers are in good agreement with the predicted probabilities based on Eqs. (7) and (8) and numerical simulations including noise. Both the qualitative trend and quantitative comparison between experiment and simulation indicate high-fidelity preparation of the low-energy singlet eigenstate for the homogeneously coupled spin chain and a high-fidelity characterization method based on singlet-triplet readout and correlations therein.

VII. CONCLUSION AND OUTLOOK
In summary, we have implemented a quantum simulation of a quantum coherent antiferromagnetic Heisenberg spin chain. For this purpose, we have developed energy spectroscopy methods to identify the condition of homogeneous exchange couplings. Furthermore, we have devised methods to prepare the global ground state of the homogeneous Heisenberg spin-chain Hamiltonian and to probe this state via local measurements in the singlet-triplet basis and correlations of such measurements. We find both qualitative and quantitative agreement between experiment and numerical simulation. Finally, coherent oscillations of the full four-spin system indicate the coherent nature of the system, despite the presence of hyperfine noise in the GaAs host material.
Future quantum magnetism simulation experiments with quantum dots may leverage the recent developments of (isotopically purified) silicon and germanium as host materials, due to the lower concentration of nuclear spins, which further enhances coherence and facilitates highresolution spectroscopy. The demonstrated control of exchange couplings, as facilitated by the independent control with virtual gates, is a powerful technique for quantum simulations in larger systems. The techniques demonstrated here also pave the way for quantum magnetism simulations in other lattice configurations, such as square and triangular ladders for which simulations of, respectively, resonating valence bonds and frustrated magnetism are now within the capabilities of gate-defined quantum dots.
The data reported in this paper and scripts to generate the figures are uploaded to Zenodo [42]. The two-spin singlet and triplet projection operators are, respectively, The global spin-raising operatorŜ þ ¼ P

commutes with the singlet and triplet projection operators.
In addition, the spin-raising operator commutes with the Heisenberg Hamiltonian thus, the spin-raising and -lowering operators map between states in the same eigenspace of the Hamiltonian. From the commutativity, it follows that two-spin singlettriplet measurements cannot distinguish states in the same eigenspace for the Heisenberg Hamiltonian, and those states can be mapped onto one another with spin-raising or -lowering operators.

Device
The quadruple quantum dot and sensing dot are formed in a device designed for eight dots and two sensors. A scanning electron micrograph image of the active region of a device similar to the one used in the experiment is shown in Fig. 1(a). The device is mounted in a dilution refrigerator, which results in an electron reservoir temperature of about 100 mK (roughly 10 μeV). By applying voltages on the electrodes on the surface, we shape the potential landscape in a two-dimensional electron gas 90 nm below the surface, formed in a silicon-doped GaAs=AlGaAs heterostructure. The plunger gates, labeled P i [red in Fig. 1(a)] for the spin-chain dots and SDP for the sensing dot, control the electrochemical potentials, and the barrier gates [green in Fig. 1(a)] control the tunnel couplings between dots or between a dot and a reservoir. When an external magnetic field is applied, it is oriented in the plane of the 2D electron gas. Table II shows an overview of the pulse sequence, with the durations and descriptions of the pulse stages.

Virtual gates
For the independent control of site-specific offsets [43], the cross talk is compensated with the matrix transformation The lever arms for the single-particle energy offsets are ½76; 81; 87; 84 μeV mV −1 , which are obtained with photonassisted tunneling [31,44].
3. Charge-stability diagrams with sequence details Figure 7 shows charge-stability diagrams for the nearestneighbor pairs in the quadruple dot. Typical pulse voltage positions and the order of the pulse sequence are indicated. The compensation and the parking for the readout of the left and right pair are not displayed to preserve clarity. Figure 8 shows single-shot results of readout directly after initialization. The data for the readout after the evolution stage are postselected by thresholding the signal from the readout after initialization. The signal is also used for background subtraction to suppress the effect of lowfrequency charge noise on the signal for the final readout. Errors in initialization for both the singlet-singlet product state and the triplet-singlet product state can occur because . Such errors can be reduced by increasing the duration of the initialization stage or increasing the dot-reservoir tunnel rate. Additionally, errors in the triplet-singlet product initialization can be caused by thermal excitations, due to which a singlet-singlet product state in the (1102) charge occupation can be occupied, which results in the bottom-left peak in Fig. 8(b). This error can be reduced by increasing the magnetic field strength. Figure 9 shows separate exchange coupling measurements for each of the neighboring pairs. For the spin funnel measurements for the left and right pair, the middle exchange coupling is set to be small. The spin funnels are fitted with the exchange model in Eq. (6). For the middle pair, the Fourier transform of coherent oscillations in the triplet-subspace is used. The model for the Fourier transform follows from Eq. (4) and is 1 2

Separate exchange coupling measurements
the homogeneous exchange coupling as obtained from Fig. 3. From these three separate fits, the tunnel couplings are extracted as 8.5, 7.5, and 11.9 μeV for the left, middle, and right pair, respectively. These tunnel coupling values are used for the numerical simulations of the experiment.

Readout with relaxation and histogram models
Single-shot readout characteristics for each of the nearest-neighbor pairs are shown in Fig. 10. The data in Figs. 10(b) and 10(d) are from the data for Fig. 5(a), and the data from Fig. 10(e) is from part of the data for Fig. 2(a).
The effect on the probabilities of relaxation during the readout is taken into account by modeling the single-shot data with a histogram. For readout on a single pair, the model for the histogram is described in Refs. [33,45]. For with N tot the number of single-shot repetitions, P ij the average probability for outcome ij, n ij the probability density distribution, and w i the bin widths. The readout of the left and right pair of spins is performed with the same sensing dot; thus, noise on the sensor signal can induce correlations between the two readouts. Low-frequency noise, which causes signal differences between repetitions, is taken into account by the subtraction of the sensor signal from the readout directly after the initialization. Correlations in the sensor signals due to higher-frequency noise remain, which results in the diagonally elongated signal peaks in the two-dimensional readout histogram as shown in Fig. 10(e). The effect of the correlations is incorporated in the model for the twodimensional Gaussian histogram by modeling the Gaussian peaks with rotated ellipses.
The two-dimensional Gaussian is with μ the mean coordinates and where the parameters for the shapes of the Gaussian peaks are with θ the rotation angle of the ellipsoidal shape of the Gaussian peaks, and σ 1 and σ 2 describe the ellipse width and length. In the experiment, the first and second readout signals are integrated for equal durations, which sets θ ¼ π=4, which means that the one-dimensional histograms for the left-and right-pair readout have the same Gaussian widths. In the extreme when there are no noise correlations, then σ 1 ¼ σ 2 , and, in the other extreme where the noise would be fully correlated, the histogram effectively is a one-dimensional Gaussian. The probability distributions for the correlated singlettriplet outcomes are with α i ¼ t i =T 1;i , which is the ratio between the signal integration time t i and the relaxation time T 1;i for the left pair or the right pair. In n TT , the first term corresponds to states which do not decay during both readouts, the second term to states which decay during the left-pair readout but not the right-pair readout, the third term to states which do not during the left-pair readout but do decay during the right-pair readout, and the last term to states which decay during both readouts. For the fitting procedure, first the histogram of singleshot results for all pulse set points is fitted to obtain the Gaussian widths and peak positions, and then the histogram for each pulse set point is separately fitted to obtain the probability amplitudes.

Automatic pulse offset calibration
The spin-chain operation of the device is automatically recalibrated on a daily basis to correct for the effect of irregular charge jumps. The results of the calibration are incorporated as offsets for the pulse voltages. The static voltage on the sensing dot plunger is sometimes adjusted between measurements, but all other dc voltages are left untouched to avoid introducing instabilities of the device. Figure 11 shows the result of the daily automatic recalibration, which consists of five traces. These traces are a sweep of the sensing dot plunger gate voltage and four one-dimensional cuts through the charge-stability space of the quadruple quantum dot. The sensing dot scan is performed to calibrate the gate voltage for the position on the flank of the sensing dot Coulomb peak. The four one-dimensional cuts are performed to locate the gate voltages corresponding to specific charge transitions. The gate voltages for these transitions form a reference for the voltages used in the pulsed control of the spin chain, such that electrostatic shifts of the device are corrected for. The analysis procedures for the Coulomb peak and the dot-reservoir transition are based on Ref. [46], and the interdot transition fitting is from Ref. [44]. The automatic calibration routine takes in total approximately 30 s. Occasionally, large charge jumps require coarse manual tuning of the device, after which the automatic calibration routine is used for refinement.

APPENDIX D: HYPERFINE GRADIENT COUPLINGS
The effect of the nuclear spins in the material environment on the electron spins in the quantum dots can effectively be described by hyperfine fields [32] as with ⃗ h i the local hyperfine field for dot i. The hyperfine term breaks the conservation of total spin and spin z of the Heisenberg Hamiltonian and has been studied extensively experimentally and theoretically for two-spin systems [32,47]. For the energy spectroscopy of the fourspin system as shown in Fig. 2, the couplings between the singlet-subspace and the polarized states are given by the hyperfine matrix, which in the basis ðQ þþ ; Q þ ; j2 T þ i; j1 T þ i; j0 T þ i; j1 S i; j0 S iÞ is with dB þ;ij ¼dB x;ij þidB y;ij and dB −;ij ¼ dB x;ij − idB y;ij , where dB α;ij ¼ ðB α;i − B α;j Þ=2, and B α;sum ¼ P i B α;i . From the hyperfine matrix, it follows that, for any state in the singlet subspace jS k i, the first-order coupling to a quintuplet is hS k jH hf jQ þþ i ¼ hS k jH hf jQ þ i ¼ 0. The coupling between the singlet and the quintuplet(s), as visible in the energy spectroscopy in Fig. 2, can be explained by a second-order effect where the singlet state couples to a triplet state, which, in turn, couples to a quintuplet state. The second-order coupling becomes effective only when the singlet and quintuplet state energies are near the energy of a coupling-mediating triplet state. This is observed in the energy spectroscopy where the coupling of the singlet to the Q þþ state decreases as the singlet energy differs more from the T þ 0 and T þ 1 energies. Note that the spin-orbit interaction can induce couplings between states with different charge occupation [48,49], which can contribute to the signal for the anticrossings in the energy spectroscopy. where ϵ i is the negative single-particle energy offset, n i ¼ c † i c i is the dot occupation, c ð †Þ i is the annihilation (creation) operator, U i is the on-site Coulomb repulsion, V ij is the intersite Coulomb repulsion, and t ij is the interdot tunnel coupling. The parameter values for the simulations are U i ¼ 3 meV, V i;iþ1 ¼0.5meV, V i;iþ2 ¼ 0.1 meV, V i;iþ3 ¼ 20 μeV, and ½t 12 ; t 23 ; t 34 ¼ ½8.5; 7.5; 11.9 μeV. The tunnel coupling values are experimentally obtained from fits to spin funnels for the left and right pair and from a fit of the Fourier transform in Fig. 3(d) (see Appendix C 5). Differences between the experimental and simulation values for the interaction parameters are accounted for in the simulations with offsets in the single-particle energies. The Hamiltonian matrix for the simulation is generated with QuTip [50]. For the simulation, the same sequence is followed as the experimental sequence for the spin chain, which is shown in Fig. 1(b) and detailed in Table II. Time evolution of the states is computed using an in-house density matrix solver package [51]. The simulated probabilities shown in the figures in the main text are the probabilities at the isolation stage after the evolution. In order to simulate the adiabaticity of the experimental sequence, the voltage pulse shape for the simulation is based on the experimental voltage pulse shape from the arbitrary waveform generator (AWG). The simulated shape is obtained from the Fourier transform of the ideal pulse shapes and the subsequent inverse Fourier transform with the experimentally measured frequency response of the AWG. The Fourier component amplitudes are corrected for the finite sampling rate f s ¼ 1 GHz of the AWG by dividing with sincðf=f s Þ, where f is the frequency of the Fourier component.
The effects of charge noise and quasistatic hyperfine noise are included for the simulations in the main text, except for Fig. 3. The quasistatic hyperfine noise is considered by repeating the simulation and for each repetition taking a sample from a Gaussian distribution with a root mean square of 3.2 mT, which is experimentally obtained from the amplitude decay of the global exchange oscillations. The power spectral density of the charge noise is modeled with A=f α , with A ¼ 0.26 μeV 2 =Hz and α ¼ 0.79, which are obtained from a charge-noise measurement. The charge noise is measured from 1 Hz to 5 kHz. For the simulations, a charge-noise frequency range from 0.1 Hz to 100 GHz is used, where noise on a timescale longer than a single shot of the experimental sequence is integrated and added as quasistatic noise.