Exploring quantum chaos with a single nuclear spin

TL;DR

This paper proposes using a single nuclear spin in silicon to experimentally observe quantum chaos, enabling detailed study of the quantum-classical boundary and dynamical tunneling with long coherence times.

Contribution

It introduces a realistic method to realize a quantum chaotic system using a nuclear spin in silicon, facilitating direct observation of quantum chaos phenomena.

Findings

Signatures of chaos are expected with feasible parameters.

The system allows exploration of quantum decoherence and classical chaos.

Potential to observe dynamical tunneling in a single-particle quantum system.

Abstract

Most classical dynamical systems are chaotic. The trajectories of two identical systems prepared in infinitesimally different initial conditions diverge exponentially with time. Quantum systems, instead, exhibit quasi-periodicity due to their discrete spectrum. Nonetheless, the dynamics of quantum systems whose classical counterparts are chaotic are expected to show some features that resemble chaotic motion. Among the many controversial aspects of the quantum-classical boundary, the emergence of chaos remains among the least experimentally verified. Time-resolved observations of quantum chaotic dynamics are particularly rare, and as yet unachieved in a single particle, where the subtle interplay between chaos and quantum measurement could be explored at its deepest levels. We present here a realistic proposal to construct a chaotic driven top from the nuclear spin of a single donor atom in silicon, in the presence of a nuclear quadrupole interaction. This system is exquisitely measurable and controllable, and possesses extremely long intrinsic quantum coherence times, allowing for the observation of subtle dynamical behavior over extended periods. We show that signatures of chaos are expected to arise for experimentally realizable parameters of the system, allowing the study of the relation between quantum decoherence and classical chaos, and the observation of dynamical tunneling.

Tables

Table 1. Table 1: Parameters of group-V donors in silicon. Ground 1 ​ s 1 𝑠 1s state binding energies taken from Ref.  Ramdas and Rodriguez, 1981 . Hyperfine interaction A 𝐴 A taken from Ref.  Feher, 1959 . Nuclear gyromagnetic ratio γ n = μ ​ μ n / I subscript 𝛾 n 𝜇 subscript 𝜇 n 𝐼 \gamma_{\mathrm{n}}=\mu\mu_{\mathrm{n}}/I is calculated from the nuclear magnetic moment μ 𝜇 \mu given in units of nuclear magneton μ n = subscript 𝜇 n absent \mu_{\mathrm{n}}= 7.62   MHz / T times 7.62 divide megahertz tesla 7.62\text{\,}\mathrm{MHz}\text{/}\mathrm{T} in Ref.  Stone, 2014 . Minimum and maximum values of nuclear quadrupole moment Q n subscript 𝑄 n Q_{\mathrm{n}} are given, based on the range of values reported in Ref.  Stone, 2014 .

donor	$I$	1s binding energy	$A$	${\gamma }_n$	$Q_n$
		($\mathrm{meV}$)	($\mathrm{MHz}$)	($\mathrm{MHz}\text{/}T$)	$\left({10}^{-28}{m}^2\right)$
${}_{}{}^{31}P$	$1/2$	45.59	117.53	17.26	-
${}_{}{}^{75}\mathrm{As}$	$3/2$	53.76	198.35	7.31	0.314
${}_{}{}^{121}\mathrm{Sb}$	$5/2$	42.74	186.80	10.26	[-0.36 , -0.54]
${}_{}{}^{123}\mathrm{Sb}$	$7/2$	42.74	101.52	5.55	[-0.49 , -0.69]
${}_{}{}^{209}\mathrm{Bi}$	$9/2$	70.98	1475.4	6.96	[-0.37 , -0.77]

Table 2. Table 2: Comparison of equivalent classical and quantum Hamiltonian parameters.

classical		quantum
original	dimensionless	original	dimensionless
$𝐋$	${𝐋}^{\prime }=\frac{𝐋}{\left|𝐋\right|}$	$𝐈$	${𝐈}^{\prime }=\frac{𝐈}{I}$
$\frac{\alpha }{2\pi }$	1	${\gamma }_n{B}_0$	1
$\frac{\beta }{2\pi }$	${\beta }^{\prime }=\frac{\beta \left|𝐋\right|}{\alpha }$	$Q$	$Q^{\prime }=\frac{QI}{{\gamma }_n{B}_0}$
$\frac{\gamma }{2\pi }$	${\gamma }^{\prime }=\frac{\gamma }{\alpha }$	${\gamma }_n{B}_1$	$B_1^{\prime }=\frac{B_1}{B_0}$
$f$	$f^{\prime }=\frac{2\pi f}{\alpha }$	$f$	$f^{\prime }=\frac{f}{{\gamma }_n{B}_0}$
$t$	$t^{\prime }=\frac{\alpha }{2\pi }t$	$t$	$t^{\prime }={\gamma }_n{B}_{0}t$

Equations

$${\mathcal{H}}_{\mathrm{classical}}=\alpha L_{\mathrm{z}}+\beta L_{\mathrm{x}}^2+\gamma \cos \left(2\pi ft\right)L_{\mathrm{y}},$$

$$\mathcal{H}=$$

$$+\left({\gamma }_{\mathrm{e}}{S}_{\mathrm{y}}-{\gamma }_{\mathrm{n}}{I}_{\mathrm{y}}\right)B_1\cos \left(2\pi ft\right),$$

$${\mathcal{H}}_{\mathrm{Q}}=Q\left(I_{{\mathrm{z}}^{\prime }}^2-\frac{I^2}{3}+\frac{\eta }{3}\left(I_{{\mathrm{x}}^{\prime }}^2-I_{{\mathrm{y}}^{\prime }}^2\right)\right),$$

$${\mathcal{H}}_{\mathrm{quantum}}=\left({\gamma }_{\mathrm{n}}{B}_0\pm \frac{1}{2}A\right)I_{\mathrm{z}}+Q{I}_{\mathrm{x}}^2+{\gamma }_{\mathrm{n}}{B}_0\cos \left(2\pi ft\right)I_{\mathrm{y}}.$$

$${\mathcal{H}}_{\mathrm{quantum},\mathrm{RF}}={\gamma }_{\mathrm{n}}{B}_0{I}_{\mathrm{z}}+Q{I}_{\mathrm{x}}^2+$$

$$\left[{\gamma }_{\mathrm{n}}{B}_{1,\mathrm{I}}\cos \left(2\pi f_{\mathrm{RF}}t\right)+{\gamma }_{\mathrm{n}}{B}_{1,\mathrm{Q}}\cos \left(2\pi ft\right)\sin \left(2\pi f_{\mathrm{RF}}t\right)\right]I_{\mathrm{y}}.$$

$${\mathcal{H}}_{\mathrm{quantum},\mathrm{RWA}}=$$

$$-\frac{1}{2}Q{I}_{\mathrm{z}}^2+\frac{1}{2}{\gamma }_{\mathrm{n}}{B}_{1,\mathrm{Q}}\cos \left(2\pi ft\right)I_{\mathrm{x}},$$

$$\frac{{\mathcal{H}}_{\mathrm{c}}}{2\pi |\mathbf{L}|}$$

$$=\frac{\alpha }{2\pi }{L}_{\mathrm{z}}^{\prime }+\frac{\beta }{2\pi }|\mathbf{L}|{L_{\mathrm{x}}^{\prime }}^2+\frac{\gamma }{2\pi }\cos \left(2\pi ft\right)L_{\mathrm{y}}^{\prime },$$

$${\mathcal{H}}_{\mathrm{c}}^{\prime }=\frac{{\mathcal{H}}_{\mathrm{c}}}{\alpha |\mathbf{L}|}$$

$$=L_z^{\prime }+{\beta }^{\prime }{L_{\mathrm{x}}^{\prime }}^2+{\gamma }^{\prime }\cos \left(2\pi f^{\prime }{t}^{\prime }\right)L_{\mathrm{y}}^{\prime }$$

$$\frac{{\mathcal{H}}_{\mathrm{q}}}{hI}={\gamma }_{\mathrm{n}}{B}_0{I}_{\mathrm{z}}^{\prime }+QI{I_{\mathrm{x}}^{\prime }}^2+{\gamma }_{\mathrm{n}}{B}_1\cos \left(2\pi ft\right)I_{\mathrm{y}}^{\prime },$$

$${\mathcal{H}}_{\mathrm{q}}^{\prime }=\frac{{\mathcal{H}}_{\mathrm{q}}}{hI{\gamma }_{\mathrm{n}}{B}_0}$$

$$=I_{\mathrm{z}}^{\prime }+Q^{\prime }{I_{\mathrm{x}}^{\prime }}^2+B_1^{\prime }\cos \left(2\pi f^{\prime }{t}^{\prime }\right)I_{\mathrm{y}}^{\prime }$$

$$\frac{\mathrm{d}\mathbf{L}}{\mathrm{d}t}=\{\mathbf{L},\mathcal{H}\}+\frac{\partial \mathbf{L}}{\partial t}$$

$${\dot{L}}_{\mathrm{x}}=$$

$${\dot{L}}_{\mathrm{y}}=$$

$${\dot{L}}_{\mathrm{z}}=$$

$${\dot{L}}_{\mathrm{y}}=$$

$$=$$

$${\dot{L}}_{\mathrm{x}}=$$

$${\dot{L}}_{\mathrm{y}}=$$

$${\dot{L}}_{\mathrm{z}}=$$

$${\mathcal{H}}_{\mathrm{quantum},\mathrm{RF}}={\gamma }_{\mathrm{n}}{B}_0{I}_{\mathrm{z}}+Q{I}_{\mathrm{x}}^2+$$

$$\left({\gamma }_{\mathrm{n}}{B}_{1,\mathrm{I}}\cos \left(2\pi f_{\mathrm{RF}}t\right)+{\gamma }_{\mathrm{n}}{B}_{1,\mathrm{Q}}\cos \left(2\pi ft\right)\sin \left(2\pi f_{\mathrm{RF}}t\right)\right)I_{\mathrm{y}}$$

$${\mathcal{H}}_{\mathrm{quantum},\mathrm{RF}}^{\prime }=R\left(-{\omega }_{\mathrm{RF}}t\right){\mathcal{H}}_{\mathrm{quantum},\mathrm{RF}}R\left({\omega }_{\mathrm{RF}}t\right)-f_{\mathrm{RF}}{I}_{\mathrm{z}}$$

$$R\left(\varphi \right)=e^{-i\varphi I_{\mathrm{z}}}=e^{-i{\omega }_{\mathrm{RF}}t{I}_{\mathrm{z}}}$$

$$R\left(-\varphi \right)I_{\mathrm{z}}R\left(\varphi \right)$$

$$R\left(-\varphi \right)I_{\mathrm{y}}R\left(\varphi \right)$$

$$R\left(-\varphi \right)I_{\mathrm{x}}R\left(\varphi \right)$$

$$R\left(-\varphi \right)I_{\mathrm{z}}^{2}R\left(\varphi \right)$$

$$R\left(-\varphi \right)I_{\mathrm{y}}^{2}R\left(\varphi \right)$$

$$-\frac{1}{2}\cos 2\varphi \left(I_{\mathrm{x}}^2-I_{\mathrm{y}}^2\right)+\frac{1}{2}I\left(I+1\right)$$

$$R\left(-\varphi \right)I_{\mathrm{x}}^{2}R\left(\varphi \right)$$

$$+\frac{1}{2}\cos 2\varphi \left(I_{\mathrm{x}}^2-I_{\mathrm{y}}^2\right)+\frac{1}{2}I\left(I+1\right)$$

$${\mathcal{H}}_{\mathrm{quantum},\mathrm{RF}}^{\prime }=-\frac{1}{2}Q{I}_{\mathrm{z}}^2-\frac{1}{2}Q\sin 2{\omega }_{\mathrm{RF}}t\{I_{\mathrm{x}},I_{\mathrm{y}}\}$$

Full text

I Introduction

I.1 Quantum chaos

The correspondence principle, as formulated by the Copenhagen school of quantum mechanics, states that the dynamics of quantum systems should converge towards classical dynamics, in the limit where the system becomes large. Appealing (and, for simple cases, often correct) as it may sound, this point of view is afflicted by a plethora of complications and controversies around the precise nature of the quantum-classical transition Leggett (2002), such as decoherence Zurek (2003) and the quantum measurement problem Schlosshauer (2005). Another key aspect of the quantum-classical transition concerns reconciling the chaotic dynamics of certain classical systems with the unitary evolution of their quantum-mechanical counterparts.

Classical chaos is ubiquitous and well understood. It arises from nonlinear terms in the equations of motion, and from the lack of a sufficient number of constants of motion compared to the number of degrees of freedom of the system Gutzwiller (2013). The hallmark of chaotic dynamics is the extreme sensitivity to initial conditions, whereby the trajectory of a system prepared in two infinitesimally different states evolves along two trajectories that diverge exponentially. Chaos plays a fundamental role in, for example, establishing the validity of classical statistical mechanics and thermodynamics, and a practical role in a wide range of applications, from weather forecasting to the design of tokamaks for nuclear fusion.

The usual description of quantum systems, in terms of states vectors that evolve according the Schrödinger equation, can appear puzzling when examined in the context of the chaotic behavior of the equivalent classical Hamiltonian. Consider for example two slightly different quantum states at time $t=0$, $|\psi_{1}(0)\rangle$ and $|\psi_{2}(0)\rangle$ having an initial overlap $\Omega(0)=|\langle\psi_{1}(0)|\psi_{2}(0)\rangle|^{2}=1-\delta^{2}$, with $\delta\ll 1$. As time progresses, these states evolve according to the time evolution operation $U(t)$. The overlap at later times is thus $\Omega(t)=|\langle\psi_{1}(0)|U^{\dagger}(t)U(t)|\psi_{2}(0)\rangle|^{2}$. Since the time evolution is unitary, $U^{\dagger}(t)=U^{-1}(t)$, we find that $\Omega(t)=\Omega(0)=1-\delta^{2}$; i.e., the overlap remains constant at all times. Since the exponential divergence of trajectories typical of classical systems appears ruled out, does this mean that there cannot be chaos in quantum dynamics?

A more appropriate and illuminating comparison between classical and quantum dynamics is obtained by describing the classical system in terms of a density $f$ in phase space, and calculating its time evolution using the Liouville equation $i\frac{\partial f(t)}{\partial t}=\mathcal{L}f(t)$, where $\mathcal{L}$ is the Liouville operator Peres and Terno (1996). One then finds that, given two initially overlapping densities $f_{1}(0)$ and $f_{2}(0)$, the Liouville equation for a conservative Hamiltonian system ensures that their overlap remains constant at all times Koopman (1931). This property mirrors the quantum behavior described earlier, so now the question may be reversed: in what way, if at all, does classical chaos differ from the dynamics of quantum systems?

The answer to this question can be rather subtle. At its heart, quantum mechanics requires that the classical phase space is coarse-grained into volumes of size $\hbar^{N}$ (with $N$ the number of degrees of freedom), and forbids specifying the state of the system to a precision finer than that. An illuminating example of how this affects the dynamics of chaotic systems was provided by Korsch and Berry Korsch and Berry (1981), who analyzed a classically chaotic iterative map while varying the value of $\hbar$. In the classical limit ($\hbar\rightarrow 0$), the map ‘shreds’ the initially smooth distribution into thin chaotic-looking ‘tendrils’. Conversely, when $\hbar$ becomes sizable on the scale of the map’s effect, one finds that the distribution remains smooth and seems to lose its chaotic features, displaying instead a slow spread from the initial shape. More generally, since bounded quantum systems have a discrete spectrum, their dynamics exhibits quasiperiodic features that are at odds with the ‘true chaos’ seen in classical systems Casati and Chirikov (1995). Therefore, it is often stated that ‘quantum chaos’ constitutes a new and unique type of dynamics Berry (1987). Finally, quantum systems allow for interference effects that result in peculiar dynamical features, such as the development of structures in phase space at a scale smaller than the Planck constant, which are reached most quickly when the classical system is chaotic Zurek (2001).

The issue of how to observe and interpret signatures of chaos in quantum mechanics has profound repercussions on many important topics in physics. For example, classical chaos underpins the ergodic hypothesis in statistical mechanics, and it is expected that its quantum equivalent plays a fundamental role in the thermalization of isolated quantum systems Deutsch (1991); Srednicki (1994); Rigol et al. (2008); Eisert et al. (2015). Chaos is also thought to be related to the issue of decoherence Zurek and Paz (1994), which is crucial in the modern topic of quantum information science. There, one must answer the delicate question of whether an onset of chaos may harm the operation of a large-scale quantum computer Georgeot and Shepelyansky (2000); Silvestrov et al. (2001); Braun (2002). On the other hand, it has been suggested that the inherent ability of chaotic systems to quickly explore a vast configuration space can be used for the purpose of demonstrating ‘quantum supremacy’ in multiqubit devices without error correction Boixo et al. (2018).

I.2 Experimental tests of quantum chaos

Despite its broad and deep importance, experimental progress in ‘quantum chaos’ is rare. Early work focused on the study of static and statistical properties of chaotic systems Haq et al. (1982); Leviandier et al. (1986); Held et al. (1998), such as the energy spectra of chaotic billiards implemented in semiconductor quantum dots Marcus et al. (1992); Wilkinson et al. (1996); Folk et al. (1996). Even more rare is the ability to experimentally observe dynamical chaos, i.e. signatures of chaos in the time evolution of quantum systems. Crudely speaking, this is because most quantum systems decohere and randomize for trivial reasons (noise, uncontrolled environments, etc.) over time scales that are too short for signatures of chaotic behavior to reveal themselves. Conversely, systems with long coherence times, such as ensembles of nuclear spins in liquids, can show signatures of chaos in the dynamics of their macroscopic magnetization Lin et al. (2000). The experimental state of the art for truly quantum chaotic dynamics is found in ensembles of cold gases Hensinger et al. (2001); Steck et al. (2001); Chaudhury et al. (2009); Manai et al. (2015), whereas only very recently an experiment on three superconducting qubits has provided experimental insight into the link between chaos and thermalization in a small-scale quantum system Neill et al. (2016).

What is still missing is an experimental study of the quantum signatures of chaos in an individual quantum system. Such a study will be an important complement and extension to experiments conducted on ensembles of particles, since an individual quantum system allows a much broader choice of measurement strategies. Although the chaotic dynamics we will describe in this paper are the result of the Hamiltonian evolution alone, the use of an individual quantum system will allow us in the future to explore the interplay between the emergence of chaos and the measurements performed on the system. Theoretical studies Sanders and Milburn (1989); Eastman et al. (2017) predict that the measurement strength can be used as an additional experimental knob to tune the chaoticity of the system’s dynamics. The measurement strength on a single object can be tuned continuously Katz et al. (2006); Hatridge et al. (2013) from projective single-shot readout Blatt and Zoller (1988); Elzerman et al. (2004); Vijay et al. (2011) to arbitrarily weak measurements, partial wave function collapse Guerlin et al. (2007) and even measurement reversal Katz et al. (2008). Strings of individual measurement outcomes could be analyzed with sophisticated statistical techniques to extract the most accurate information on the trajectory of the quantum object Gammelmark et al. (2013), providing unprecedented insights into the chaotic dynamics of a monitored quantum system. Variable-strength measurements have already been experimentally demonstrated in the 31P donor system Muhonen et al. (2017).

Moreover, in the realm of classical computation, it has been recently shown that a network of individually chaotic electronic components can solve computationally hard problems faster than one using nonchaotic elements Kumar et al. (2017). The construction of a single quantum chaotic system amenable to networking and controlled interactions could provide insights into how the equivalent quantum circuit would perform in complex computational problems.

I.3 Quantum chaotic driven top

Here we present a detailed and quantitative proposal to experimentally realize a single-atom version of one of the best studied quantum chaotic systems, the ‘kicked top’ Haake (2000). For experimental convenience, we will focus on the case where the top is periodically driven, instead of kicked with $\delta$ functions. This system becomes chaotic in the presence of a term in the Hamiltonian that is quadratic in the angular momentum, and has a classical Hamiltonian of the form

[TABLE]

with angular momentum $\mathbf{L}=\left(L_{\mathrm{x}}\hskip 5.0ptL_{\mathrm{y}}\hskip 5.0ptL_{\mathrm{z}}\right)^{\mathrm{T}}$ ($|\mathbf{L}|=$ constant), $\alpha$, $\beta$, and $\gamma$ proportionality constants, and $f$ the frequency of the drive. The size and shape of the regular and chaotic regions in classical phase space are determined by the Hamiltonian parameters $\alpha$, $\beta$, $\gamma$, and $f$ (Fig. 1). A sizable region of chaos is found with linear and quadratic interactions of similar strength ($\beta\approx\alpha$), a sufficiently strong periodic drive ($\gamma\gtrsim 0.02\alpha$), and a drive frequency close to the range of resonance frequencies of the nonlinear system ($f\approx 1.4\alpha$). These are the conditions we seek to reproduce in the quantum driven top.

The obvious quantum equivalent of a classical spinning top is a spin. The challenge here is to find a spin system whose Hamiltonian maps onto that of the chaotic driven top. This requires in particular a quadratic term in the Hamiltonian, which is only possible for a spin quantum number $I>1/2$. Moreover, a larger spin is crucial in comparing its dynamics to the structure of the classical phase space, as its smaller relative uncertainty spread allows for better localization of the quantum state in a certain area of interest in classical phase space (see Appendix E for further details). Using a single spin with access to high-fidelity single-shot state readout and/or variable strength weak measurements opens up the largely unexplored area of the interplay between chaos and quantum measurements. Lastly, it is of paramount importance that this system does not lose coherence for trivial reasons, unrelated to chaos, on time scales short compared to the chaotic dynamics. This requires a long intrinsic quantum coherence time of the system.

I.4 Experimental platform

Our proposed system meeting these requirements is the nuclear spin of a heavy group-V substitutional donor in isotopically enriched 28Si Itoh and Watanabe (2014). The lightest group-V donor in silicon, 31P, has been extensively studied in the context of quantum information processing Kane (1998), since it naturally contains two quantum bits, the electron (with spin $S=1/2$) and the 31P nucleus (with spin $I=1/2$). High-fidelity single-shot readout Morello et al. (2010); Pla et al. (2013), coherent operation Pla et al. (2012, 2013), mutual entanglement Dehollain et al. (2016a) and variable-strength measurements Muhonen et al. (2017) have been experimentally demonstrated. When implanted Van Donkelaar et al. (2015) in isotopically enriched 28Si, these single-atom spins exhibit outstanding coherence times Muhonen et al. (2014) (up to 35 seconds for the nuclear spin) and control fidelities Muhonen et al. (2015); Dehollain et al. (2016b). This suggests that donor spin systems would be ideal platforms to study the subtle effects of dynamical chaos and its interplay with quantum measurement, if it were possible to engineer a suitable spin Hamiltonian. This is not the case with 31P, since its spin value of 1/2 forbids the presence of quadratic terms in the Hamiltonian. Heavier donors, such as 75As, 121Sb, 123Sb and 209Bi, all have nuclear spins $I>1/2$, which allows the existence of a nuclear quadrupole interaction, scaling quadratically in the spin operators. Below we show that, under realistic conditions of quadrupole interaction and periodic drive, a heavy group-V donor can become a single-atom solid-state implementation of a chaotic driven top.

II Results

II.1 A large nuclear spin donor as quantum driven top

The spin Hamiltonian of group-V donors in silicon, in the presence of a static magnetic field $B_{0}$ in the $z$ direction and an oscillating magnetic field $B_{1}$ at frequency $f$ in the $y$ direction, reads

[TABLE]

where $\gamma_{\mathrm{e}}$ and $\gamma_{\mathrm{n}}$ are the electron and nuclear gyromagnetic ratios (their magnetic moments have opposite sign), the electron spin $S=1/2$ is described by the vector of operators $\mathbf{S}=\left(S_{\mathrm{x}}\hskip 5.0ptS_{\mathrm{y}}\hskip 5.0ptS_{\mathrm{z}}\right)^{\mathrm{T}}$, the nuclear spin $I$ is described by $\mathbf{I}=\left(I_{\mathrm{x}}\hskip 5.0ptI_{\mathrm{y}}\hskip 5.0ptI_{\mathrm{z}}\right)^{\mathrm{T}}$, $A$ is the hyperfine interaction between electron and nuclear spin (assumed to be isotropic), and $\mathcal{H}_{\mathrm{Q}}$ accounts for the nuclear quadrupole interaction (discussed below).

Experimentally implemented via a broadband on-chip antenna, the oscillating $B_{1}$ field allows for coherent control of the spins through electron spin resonance (ESR) and nuclear magnetic resonance (NMR). This enables the application of numerous techniques for tomography and characterization of the spin system Wolfowicz and Morton (2007).

The hyperfine interaction $A$ couples the electron and the nuclear spins, and can be approximated by an effective interaction $AS_{\mathrm{z}}I_{\mathrm{z}}$ in the device operating regime where $\gamma_{\mathrm{e}}B_{0}\gg A$. This introduces a dependence of the ESR frequency on the state of the nucleus. In turn, this allows the measurement (and consequent initialization by measurement) of the nuclear spin state by observing at what frequency the electron spin responds to a resonant microwave excitation using electron spin read-out via a standard spin-to-charge conversion technique Pla et al. (2013) (Fig. 2).

The nuclear quadrupole moment is caused by the non-spherical charge distribution of the nucleus. This quadrupole interacts with electric field gradients to introduce a new term in the spin Hamiltonian Kaufmann and Vianden (1979). In general, the tensor describing the electric field gradient can be diagonalized to $\operatorname{diag}\left(V_{\mathrm{x^{\prime}x^{\prime}}},V_{\mathrm{y^{\prime}y^{\prime}}},V_{\mathrm{z^{\prime}z^{\prime}}}\right)$ by an appropriate choice of coordinate frame $(x^{\prime},y^{\prime},z^{\prime})$, where $V_{\mathrm{ij}}=\tfrac{\partial^{2}V}{\partial{i}\partial{j}}$, $\left(i,j\in x^{\prime},y^{\prime},z^{\prime}\right)$ are the partial second derivatives of the electrostatic potential $V$, and $|V_{\mathrm{x^{\prime}x^{\prime}}}|\leq|V_{\mathrm{y^{\prime}y^{\prime}}}|\leq|V_{\mathrm{z^{\prime}z^{\prime}}}|$. The expression for the nuclear quadrupole interaction is then simplified to Slichter (1990)

[TABLE]

where $Q$ is the effective quadrupole interaction strength, which scales linearly with both $V_{\mathrm{z^{\prime}z^{\prime}}}$ and the nuclear quadrupole moment $Q_{\mathrm{n}}$, and the asymmetry parameter $\eta=\frac{V_{\mathrm{x^{\prime}x^{\prime}}}-V_{\mathrm{y^{\prime}y^{\prime}}}}{V_{\mathrm{z^{\prime}z^{\prime}}}}$ quantifies the deviation from axial symmetry of the electric field gradient ($0\leq\eta\leq 1$, $\eta=0$ corresponds to axial symmetry; see Appendix H for further details on quadrupole interaction).

The important features of the quadrupole interaction are that it is quadratic in the spin operators and has a preferred quantization axis. It is often the case that the electric field gradient tensor has approximately axial symmetry ($\eta=0$) Slichter (1990). If, in addition, the static magnetic field $B_{0}$ can be oriented in an arbitrary direction (for example using a 3-axis vector magnet), the linear and quadratic terms in the spin Hamiltonian can be made orthogonal. Ignoring the static energy offset $QI^{2}/3$, and assuming for simplicity that the symmetry axis of the electric field gradient is orthogonal to the periodic driving field $B_{1}$, the Hamiltonian of a donor with nuclear spin $I>1/2$ takes the form

[TABLE]

Therefore, $\mathcal{H}_{\mathrm{quantum}}$ represents the quantum equivalent of the Hamiltonian of a classical periodically driven top [Eq. (1)].

II.2 Donor parameters and preferred operation regimes

II.2.1 Realising a quantum driven top in the laboratory frame

We now estimate the parameters of the spin Hamiltonian of the quantum driven top, in order to compare them to the parameters that are known to lead to chaotic dynamics in the equivalent classical case.

The key parameter values of the group-V donors are summarized in Table 1. The value of $I$ increases with atomic mass, while the hyperfine interaction $A$ has a non-monotonic behavior, with a significant jump for the heavy 209Bi donor. Large $I$ values are desirable to reduce the relative quantum uncertainty of the spin state (see Appendix E for details), whereas from the analysis of the classical driven top we know that interesting chaotic dynamics arises when linear and quadratic terms in the spin Hamiltonian are of comparable strength.

For operation in the neutral charge state, the linear term has strength $\gamma_{n}B_{0}\pm A/2$, which becomes very large in the case of ${}^{209}\mathrm{Bi}$ where $A=$1475.4\text{,}\mathrm{MHz}$$. The $\pm A/2$ contribution can be removed by operating in the ionized state Pla et al. (2013). The nuclear Zeeman term is minimized by operating at low $B_{0}$, with the caveat that reducing $B_{0}$ affects the electron readout and initialization fidelity. Using a minimum attainable value of $B_{0}=0.5$ T and ionized donors, the linear term in the spin Hamiltonian thus takes values of order 3 MHz.

Next, we wish to obtain a comparable value for the quadratic term, which is achieved by maximizing the strength of the quadrupole interaction, and hence the electric field gradient [Eq. (3)]. Recent experiments Thorbeck and Zimmerman (2015); Park et al. (2016); Pla et al. (2018) found that aluminum gates on top of a Si/SiO2 stack cause strain in the silicon; finite-element models indicate that strain values approaching 0.1% can be expected [Fig. 2 (c)]. This strain-induced disruption of the cubic lattice symmetry under the metal results in an electric field gradient at the donor site. Such an electric field gradient is further enhanced at the nucleus due to rearrangement of the closed electronic shells of the donor (the Sternheimer anti-shielding effect Kaufmann and Vianden (1979); Feiock and Johnson (1969)). These effects are present in both the neutral and the ionized charge state, albeit with different strengths.

Measurements of the quadrupole interaction strength of group-V donors in strained silicon, especially near a nanostructure interface, are limited and have only been conducted very recently Franke et al. (2015, 2016); Mortemousque et al. (2016). For ionized group-V donors, the available data only allow an order-of-magnitude estimate of $Q$ in the hundreds of kHz range (see Appendix H for an extensive review).

Upon comparing the classical and quantum systems, the coefficient $\beta$ in the classical Hamiltonian must be compared to $QI$ in the quantum case (Appendix A). This makes the hundreds of kHz range for $Q$ rather promising, since for large $I$ the corresponding $QI\sim$1\text{,}\mathrm{MHz} is within an order of magnitude of the linear interaction strength ($\gamma_{\mathrm{n}}B_{0}\sim$3\text{\,}\mathrm{MHz}). If this proves to be insufficient, strain engineering could be deployed to further increase the electric field gradient to reach the target value of $QI\sim$3\text{,}\mathrm{MHz}$$. This can be achieved either in MOSFET structures Thompson et al. (2006) or in Si/SiGe devices, where the ability to electrically detect a single dopant atom coupled to a quantum dot has also been recently demonstrated Foote et al. (2015). All the above options will deliver a fixed value of strain, set by the thermal expansion of the metallic electrodes and/or the built-in strain in the substrate. As a next step in experimental sophistication and control, one could consider fabricating on-chip piezoactuators to dynamically control the strain Dreher et al. (2011), allowing the study of the chaotic dynamics of a nuclear spin as a function of its Hamiltonian parameters. Finally, another option to tune in situ the quadrupole splitting could be to distort the electron wave function using strong voltages on gates placed above the donor, to the point where the electron wave function is significantly displaced from the nuclear site, potentially generating a substantial electric field gradient. This type of electron wave function distortion has been discussed in numerous papers Calderón et al. (2006); Lansbergen et al. (2008); Tosi et al. (2017), but no calculation of the resulting nuclear quadrupole splitting in the case of a $I>1/2$ nucleus has been performed to date.

The broadband antenna near the donor can be used to apply a radio-frequency periodic drive. Previous work has been conducted with drive strengths up to $B_{1}\sim$2\text{,}\mathrm{mT}$$ Pla et al. (2013), which correspond to radio-frequency powers of order 0.5 mW (at the chip). Those values were sufficient to achieve high-fidelity coherent control of the ${}^{31}\mathrm{P}$ nuclear spin qubit. Here, we wish to compare $\gamma_{\mathrm{n}}B_{1}$ to the classical parameter $\gamma$. Chaos arises when $\gamma\approx 0.02\alpha$ in the classical model. For $B_{0}\sim 0.5$ T and thus $\alpha\sim 3$ MHz, this implies $B_{1}\sim 10$ mT. Assuming the same setup and antenna as Ref. Pla et al., 2013, this value would require $\sim 10$ mW radio-frequency power at the chip. This is a very high value for operation at millikelvin temperatures, but we note that the broadband microwave antenna is terminated by a short circuit, constituting (ideally) a fully reflective load. Therefore, only a small fraction of the incident power is actually dissipated on the chip, while the rest is reflected and dissipated at stages of the setup with large cooling powers. Alternatively, high $B_{1}$ values with low incident power could be obtained by using $LC$ resonators.

II.2.2 Realizing a quantum driven top in the rotating frame

The Hamiltonian of the periodically driven top described so far was defined in the laboratory frame. An alternative approach is to describe it in the rotating frame, defined by the frequency of an oscillating field. This results in a system ‘dressed’ by a continuous radio-frequency field, at a frequency that matches the nuclear Zeeman interaction strength ($f_{\mathrm{RF}}=\gamma_{\mathrm{n}}B_{\mathrm{0}}$). This is a well-established method that originates from quantum optics Mollow (1969) and has recently been extended to microwave frequencies Timoney et al. (2011), including with the electron spin of the ${}^{31}\mathrm{P}$ donor Laucht et al. (2017). Here we analyze its application to the higher-dimensional nuclear spin of the heavier group-V donors. We consider a Hamiltonian for the ionized nucleus $\left(A=0\right)$ of the form:

[TABLE]

Here $B_{\mathrm{1,I}}$ and $B_{\mathrm{1,Q}}$ are the in-phase and quadrature amplitudes of an $\mathrm{IQ}$-modulated radio-frequency drive at frequency $f_{\mathrm{RF}}=\gamma_{\mathrm{n}}B_{\mathrm{0}}$, with an additional amplitude modulation applied to the in-phase component of the drive, at a frequency $f$. Switching the system to the rotating frame and applying the rotating wave approximation (RWA) will effectively remove the static Zeeman energy, while reintroducing a linear term which scales with the strength of the continuous drive at frequency $f_{\mathrm{RF}}=\gamma_{\mathrm{n}}B_{0}$ (see Appendix C for a derivation). Now the Hamiltonian reads

[TABLE]

which is, up to a trivial rotation, equivalent to the quantum driven top [Eq. (4)].

Engineering a dressed system and considering this in the rotating frame has several important benefits for the experimental feasibility of our proposal. Firstly, the linear term in the new Hamiltonian, $\alpha=\gamma_{\mathrm{n}}B_{\mathrm{1,I}}/2$, is continuously (and rapidly, if desired) tunable all the way to zero, and up to a maximum set by the strongest attainable oscillating field strength ($B_{\mathrm{1}}\sim 10$ mT as assumed earlier), corresponding to $\alpha\sim$30\text{,}\mathrm{kHz}. This means that we can access a strongly chaotic regime, where the quadratic term $\beta=QI$ is comparable to the linear term, with quadrupole interaction strengths of only $Q\sim$10\text{\,}\mathrm{kHz}.

Second, $\mathrm{IQ}$ modulation combined with amplitude modulation is a standard microwave control technique, which allows full and independent control over the strength of the periodic drive at frequency $f$ between 0 and the maximum strength of the linear interaction $\tfrac{1}{2}\gamma_{\mathrm{n}}B_{\mathrm{1,I}}=0\sim$30\text{,}\mathrm{kHz}$$. This opens up a new parameter regime of very strong periodic drive, with increased size of classical chaotic regions (see Appendix C), which would be challenging to obtain in the laboratory frame.

Lastly, the constraint of orthogonality between axis of quadrupole interaction and direction of periodic drive, as imposed by the classical system [Eq. (1)], is relaxed in the rotating frame under the RWA (see Appendix C). This is important, since now the axis of the static Zeeman field $B_{0}$ only needs to be perpendicular to the plane defined by the directions of quadrupole interaction and periodic drive, regardless of the relative angle between the latter two. This is easily achievable using a 3D vector magnet.

Overall, moving to the rotating frame can allow exploring a wider parameter space, but the actual timescale of dynamical phenomena will be scaled down by a factor $\sim 100$ (for $\tfrac{1}{2}\gamma_{\mathrm{n}}B_{\mathrm{1,I}}\sim$30\text{,}\mathrm{kHz} vs $\gamma_{\mathrm{n}}B_{\mathrm{0}}=$2.8\text{\,}\mathrm{MHz} at $B_{\mathrm{0}}=$0.5\text{,}\mathrm{T}$$). For example, the period of dynamical tunneling [Fig. 4(c) for the case of the laboratory frame] will become $\sim 100\leavevmode\nobreak\ \mu$s for the same choice of relative parameter strengths. This remains several orders of magnitude faster than the expected intrinsic coherence time of the nuclear spin, noting also that the technique of dressing a spin with a driving field often yields an extra order of magnitude in coherence time Laucht et al. (2017).

II.2.3 Summary

In summary, these estimates suggest that the ${}^{123}\mathrm{Sb}$ and ${}^{209}\mathrm{Bi}$ donors in silicon are suitable candidates to implement the quantum driven top, due to their high spin quantum number ($I=7/2$ and $I=9/2$, respectively), low nuclear gyromagnetic ratio ($\gamma_{\mathrm{n}}=$5.55\text{,}\mathrm{MHz}\text{,}{\mathrm{T}}^{-1} and $\gamma_{\mathrm{n}}=$6.96\text{\,}\mathrm{MHz}\text{\,}{\mathrm{T}}^{-1}), and, in the case of ${}^{123}\mathrm{Sb}$ , low hyperfine coupling strength ($A=$101.52\text{,}\mathrm{MHz}). ${}^{123}\mathrm{Sb}$ has the additional advantage that its suitability for ion implantation is well documented Schenkel *et al.* ([2006](#bib.bib77)): after low-energy implantation and high-temperature rapid thermal anneal, the Sb atoms are fully activated, and the implantation damage to the silicon lattice is thoroughly repaired. Recent work Weis *et al.* ([2012](#bib.bib78)) suggests that the implantation damage can be efficiently repaired also in the case of ${}^{209}\mathrm{Bi}$ , although the electrical activation yield remains lower than that of Sb. The attainable quadrupole interaction is not well known, but recent work Franke *et al.* ([2015](#bib.bib62), [2016](#bib.bib63)) indicates that it can plausibly reach a comparable value to the nuclear Zeeman term in low ($B_{0}\sim$0.5\text{\,}\mathrm{T}) static magnetic field. Furthermore, dressing the system at the nuclear Zeeman frequency $\gamma_{\mathrm{n}}B_{\mathrm{0}}$ replaces the linear interaction strength by the drive strength. This lowers the linear interaction strength, and thereby the minimum quadrupole strength, by two orders of magnitude, which brings the parameter regime of similar strength linear and quadratic interactions within reach using current device technology only. When combined with a strong oscillating magnetic field ($B_{1}\sim$10\text{,}\mathrm{mT}$$), we conclude that the parameter range where the equivalent classical driven top behaves chaotic throughout sizable areas of its phase space is within reach. In what follows, we will concentrate our discussion on the use of ${}^{123}\mathrm{Sb}$ as the model system to study quantum chaos in a single spin.

II.3 Quantum versus classical dynamics: a comparison

To illustrate the applicability of the ${}^{123}\mathrm{Sb}$ system to the study of quantum chaos, we propose two types of experiments: one aimed at finding a correspondence between the classically chaotic driven top and its quantum counterpart, the other at demonstrating a violation of classical dynamics, exposing the true quantum nature of the system. In what follows, we focus on the system in the laboratory frame, as this puts the most stringent conditions on size and shape of the chaotic region in phase space of the classical equivalent system, however; the suggested experiments are equally well applicable to the system in the rotating frame.

II.3.1 Decoherence as a precursor of chaos

Classical chaos is characterized by an extreme sensitivity to perturbations, and as a consequence, neighboring trajectories with slightly different initial coordinates rapidly diverge. This is in stark contrast with quantum dynamics, where the discrete nature of the energy spectrum results in quasiperiodic behavior, leading to partial revivals of the initial quantum state instead (see Appendix D for details). However, quantum states evolving under slightly perturbed Hamiltonians do experience divergence. Since a quantum system is never truly isolated, interactions with its environment lead to unknown perturbations of the Hamiltonian, effectively entangling the system with its environment. The unknown nature of this process translates to the quantum state losing its purity, thus decohering into a mixed state.

Crucial to the quantum driven top, certain initial quantum states are more prone to decoherence, while others remain relatively unperturbed. This behavior appears related to the high sensitivity of certain classical states to perturbations, but is caused here by a varying sensitivity of the time-evolution operator’s eigenstates to perturbations (see Appendix D for details). States containing more of these high-sensitivity eigenstates are therefore more susceptible to phase errors. In line with this picture, Zurek and co-workers Zurek and Paz (1994) predict that the rate at which different initial quantum states decohere provides a mapping to the chaotic or nonchaotic nature of the corresponding classical system, with chaotic classical regions corresponding to more rapidly decohering initial quantum states.

The driven-top system can be used to verify this prediction, both through simulations and experiments. By evolving an initial state for a certain duration using the driven-top Hamiltonian [Eq. (4)], the resulting density matrix $\rho$ provides information about the degree of decoherence through its purity $\operatorname{Tr}(\rho^{2})$. Decoherence is simulated by randomly varying a Hamiltonian parameter during the state’s evolution, and calculating the state purity from the ensemble average of many final states, each obtained with a different randomized evolution (see App. K for details). By sampling over all spin coherent states, a ‘purity map’ of the quantum driven top is obtained, which we compare to its classical counterpart (Fig. 3). The simulations highlight a correspondence between the classically-chaotic regions and quantum regions of strong decoherence, and between classically-regular regions and quantum regions of weak decoherence. To experimentally verify these predictions, we aim to prepare spin coherent states (see App. E and App. J for details), evolve the system under the driven-top Hamiltonian, and finally reconstruct $\rho_{\mathrm{final}}$ using quantum state tomography. Repeating this for different initial spin coherent states allows for experimental reconstruction of the ‘purity map’, which can then be compared to the corresponding classical phase space.

II.3.2 Dynamical tunneling

In the absence of a periodic drive, trajectories of the classical top are closed orbits confined to distinct regions in a two-dimensional phase space (the surface of a sphere due to $\left|\mathbf{L}\right|$ being a constant of motion). Upon addition of a periodic drive, this behavior is largely kept intact, except near boundaries between the regular regions, where a chaotic behavior appears. The classically-separated regions of regular motion have an analogue in the corresponding quantum system, where they can be identified as regions of weak decoherence (Fig. 3). However, there is a fundamental difference between the quantum and the classical case. In the classical system, the Kolmogorov, Arnol’d and Moser (KAM) theorem ensures that a system initially prepared within one of the regular regions will remain on a periodic orbit within such a region. The quantum system, however, cannot be precisely localized within a certain region, due to the uncertainty principle. The ‘leakage’ of the quantum wave function into a different regular region results in the phenomenon of dynamical tunneling, i.e. the tunneling of the quantum state between separate regular regions, in violation of the KAM theorem.

Dynamical tunneling manifests itself in the quantum driven top as a periodic oscillation of a spin coherent state between the two classically-regular regions associated with the quadratic interaction [Fig. 4(a)]. In contrast, a spin coherent state prepared within a classically chaotic region rapidly spreads out and shows no apparent revival to a spin coherent state [Fig. 4(b), 4(c) and movies in Supplemental Material sup ].

Numerical simulations, conducted using Hamiltonian parameters appropriate for ${}^{123}\mathrm{Sb}$ , clearly show the appearance of dynamical tunneling for spin coherent states prepared initially within the regions of classically regular periodic orbits [Fig. 4(b), 4(c)]. The predicted period of dynamical tunneling is $\sim$3\text{,}\mathrm{\SIUnitSymbolMicro s} (see Appendix [F](#A6) for dependence of tunnel rate on system parameters); this will increase by a factor $\sim 100$ to $\sim$300\text{\,}\mathrm{\SIUnitSymbolMicro s} upon considering the system in the rotating frame (assuming $\gamma_{n}B_{1,\mathrm{I}}/2\leavevmode\nobreak\ \sim$30\text{,}\mathrm{kHz} and $Q$ is reduced by a factor 100 to $\sim$8\text{\,}\mathrm{kHz}). This is a crucial result, since this period is orders of magnitude shorter than the dephasing time of a nuclear spin in ${}^{28}\mathrm{Si}$ Muhonen et al. (2014), ensuring that the coherent dynamical tunneling oscillations can be observed in an experiment over unprecedented timescales.

III Conclusion and outlook

In this paper, we have quantitatively described a proposal for the realization of a single quantum chaotic system, based upon the nuclear spin of a substitutional group-V donor in silicon. In particular, we have shown that, with realistically achievable parameters, the $I=7/2$ nucleus of a 123Sb donor can exhibit the whole spectrum of features of interest in the study of experimental quantum chaos, from state-dependent decoherence to dynamical tunneling. The experimental verification of these predictions would constitute the first observation of quantum chaos in an individual physical system. Such an achievement would reinvigorate the fundamental study of quantum-classical boundaries by providing a well-defined and exquisitely controllable experimental test bed.

In terms of applications, one could envisage laying out and operating individually chaotic 123Sb nuclei in the same types of multispin architectures that are being extensively studied in the context of quantum information processing with 31P donors Hill et al. (2015); O’Gorman et al. (2016); Pica et al. (2016); Tosi et al. (2017). Substituting a simple $I=1/2$ spin with a multilevel system like 123Sb could allow the study of quantum information processing where the information is encoded in an intrinsically chaotic system. This would be a different and complementary approach to the one taken, e.g., in superconducting systems, where it is the nature of the interaction between multiple qubits that produces a chaotic dynamics Neill et al. (2016); Boixo et al. (2018). Rather, it could constitute the quantum version of a type of analog computation that has started to show promise in the context of classical neural networks, where having individually chaotic elements can speed up the solution of complex problems Kumar et al. (2017).

The supporting data for this work are available at Research Data Australia Mourik et al. (2018).

Acknowledgments

We thank A. Laucht, F.A. Mohiyaddin, V. Schmitt, and G. Tosi for suggestions and comments. This work was funded by the Australian Research Council Discovery Projects No. DP150101863 and No. DP180100969. V.M. acknowledges support from a Niels Stensen Fellowship. V.M. and S.A. performed the numerical simulations. H.F. wrote the C code used in the classical simulations. J.J.P. provided estimates of the quadrupole interaction strength. V.M., S.A., C.H., G.J.M., J.C.M., and A.M. conceived the experimental design and measurements proposed in this paper. V.M., S.A., and A.M. wrote the paper, with contributions from all authors. A.M. supervised the project.

Appendix A Classical and quantum Hamiltonian parameters

The classical simulation results shown in Fig. 1 use dimensionless parameters which eases a direct comparison to corresponding dimensionless parameters in the quantum case. Here, their relation to the parameters in the original Hamiltonian [Eq. (1)] and to the corresponding quantum parameters is given.

A.0.1 Rescaling of classical Hamiltonian parameters

Both the variables $\left(\mathbf{L},t\right)$ and the parameters $\left(\alpha,\beta,\gamma,f\right)$ need a dimensionless equivalent. Dividing the classical Hamiltonian $\mathcal{H}_{\mathrm{c}}$ by $|\mathbf{L}|$ (including a factor $2\pi$ to convert from $\mathrm{rad}\text{\,}{\mathrm{s}}^{-1}$ to $\mathrm{Hz}$) we obtain

[TABLE]

where the normalized angular momentum variable $\mathbf{L}^{\prime}=\mathbf{L}/|\mathbf{L}|$ is introduced. Next we divide by $\alpha/2\pi$ and time $t^{\prime}=(\alpha/2\pi)t$ is introduced:

[TABLE]

This makes the parameters $\beta^{\prime}$, $\gamma^{\prime}$ and $f^{\prime}$ dimensionless and relative to $\alpha$. $\alpha$ itself has units of $\mathrm{rad}\text{\,}{\mathrm{s}}^{-1}$ and time variable $t^{\prime}$ has units of $2\pi/\alpha$.

A.0.2 Correspondence between classical and quantum Hamiltonian parameters

The above implies that the same convention of a normalized angular momentum has to be followed quantum mechanically. Hence the quantum Hamiltonian $\mathcal{H}_{\mathrm{q}}$ is divided by $h$ and transformed with $\mathbf{I^{\prime}}=\mathbf{I}/I$, thus converting to frequency units and normalizing the spin operators:

[TABLE]

This assumes units of $\mathrm{Hz}\text{\,}{\mathrm{T}}^{-1}$ for $\gamma_{\mathrm{n}}$, units of $\mathrm{Hz}$ for $Q$, and dimensionless spin operators $\mathbf{I}$. Note that after introducing $\mathbf{I}^{\prime}$ an additional factor $I$ appears in the second term due to its quadratic nature.

Next, as in the classical case, dividing by $\gamma_{\mathrm{n}}B_{0}$ and introducing time variable $t^{\prime}=\gamma_{\mathrm{n}}B_{0}t$ results in the dimensionless Hamiltonian

[TABLE]

Table 2 gives an overview of the different parameters used throughout the classical and quantum simulations.

Appendix B Derivation of classical equations of motion

The general expressions for the equations of motion of the driven top are given here. The starting point is Hamilton’s equation of motion in Poisson bracket formulation:

[TABLE]

with $\mathbf{L}=\left(L_{\mathrm{x}}\quad L_{\mathrm{y}}\quad L_{\mathrm{z}}\right)^{\mathrm{T}}$ the angular momentum vector, $\mathcal{H}$ the Hamiltonian, and $\{\mathbf{L},\mathcal{H}\}$ the Poisson bracket relation between $\mathbf{L}$ and $\mathcal{H}$. Angular momentum conservation implies $\partial\mathbf{L}/\partial t=0$. Introducing $\dot{L}_{\mathrm{i}}\equiv\mathrm{d}L_{\mathrm{i}}/\mathrm{d}t$, the equations of motion are

[TABLE]

B.1 Equations of motion for the classical driven top

The equations of motion given above can be applied to the classical driven top Hamiltonian of the main text. As an example, we derive an expression for $\dot{L_{\mathrm{y}}}$:

[TABLE]

where the product rule for Poisson brackets is used in the third line and whenever Poisson brackets are computed, the relation $\{L_{\mathrm{i}},L_{\mathrm{j}}\}=\epsilon_{\mathrm{ijk}}L_{\mathrm{k}}$, with $\epsilon_{\mathrm{ijk}}$ the Levi-Civita symbol, is used. Similarly, equations for $L_{\mathrm{x}}$ and $L_{\mathrm{z}}$ can be derived, resulting in the system of equations

[TABLE]

Appendix C Quantum system in the rotating frame and rotating wave approximation

The technique of ‘dressing’ a quantum spin state relies on applying a microwave tone with a frequency matching the dominant linear Zeeman interaction term in the system. Upon transforming the system to the rotating frame, and applying the rotating wave approximation (RWA), the original linear Zeeman interaction term disappears, and an effective linear interaction with a strength set by the microwave amplitude appears. We derive the effective ionized nuclear spin Hamiltonian of a donor in silicon using this approach; the starting point is the Hamiltonian proposed in Eq. (5),

[TABLE]

C.1 Transforming spin operators to the rotating frame

The transformation of $\mathcal{H}_{\mathrm{quantum,RF}}$ to a frame rotating with angular velocity $\omega_{\mathrm{RF}}=2\pi f_{\mathrm{RF}}$ is given by

[TABLE]

where $R$ is the rotation operator corresponding to a basis rotation over an angle $\phi$ around the $z$ axis:

[TABLE]

The remaining task is to transform the (squared) spin operators of the original Hamiltonian to the rotating frame. Using the series expansion of the matrix exponent, commutator rules of spin operators and recognizing sine or cosine series in the expansion, one can derive the following identities for rotated spin operators:

[TABLE]

With the same strategy, albeit less trivially, one can derive for the squared spin operators the following identities:

[TABLE]

where $\{I_{\mathrm{x}},I_{\mathrm{y}}\}$ denotes the anticommutator of $I_{\mathrm{x}}$ and $I_{\mathrm{y}}$.

C.2 Hamiltonian under the rotating wave approximation

Using the results of the previous section, it is straightforward to arrive at the rotating frame version of Eq. 5, which is given by:

[TABLE]

with $\omega=2\pi f$ the drive frequency to create the quantum driven top, and $\omega_{\mathrm{RF}}=2\pi f_{\mathrm{RF}}$, with $f_{\mathrm{RF}}=\gamma_{\mathrm{n}}B_{0}$ the drive frequency of the rotating frame (exactly canceling the Zeeman interaction term $\gamma_{\mathrm{n}}B_{0}I_{\mathrm{z}}$). Applying the RWA now reduces to neglecting all terms involving oscillatory factors at frequency $2\omega_{\mathrm{RF}}$. Further ignoring the irrelevant static energy offset $\dfrac{1}{2}QI\left(I+1\right)$, the Hamiltonian reduces to Eq. (6):

[TABLE]

Upon applying two trivial rotations, first by an angle $\pi/2$ around the $y$ axis, followed by an angle $\pi/2$ around the $x$ axis, the original quantum driven top Hamiltonian of Eq. (4) is recovered, demonstrating equivalence between the laboratory and rotating frame approach to creating a quantum driven top.

C.3 Relative angle between quadrupole interaction and periodic drive

The Hamiltonian of Eq. (15) still contains strong constraints on the relative directions of the different terms, as the linear Zeeman interaction, quadratic quadrupole interaction, and periodic drive are all orthogonal to each other. In a realistic experiment, the angle $\theta$ between the principal axis of the quadrupole interaction and the direction of the periodic driving field will be an intrinsic and uncontrollable device property. Under the RWA, however, the constraint of orthogonality between quadrupole interaction and periodic drive may be relaxed. Defining the $x,y$ plane as the one containing the principal axis of the quadrupole coupling and the periodic driving field, the periodic drive term in Eq. (15) may be rewritten as

[TABLE]

where a phase shift by an angle $\varphi$ is included in the drive. Upon choosing $\varphi$ to be equal to the angle $\theta$, in the rotating frame, the only oscillatory factors containing the phase $\theta$ are all terms in $2\omega_{\mathrm{RF}}$, which are neglected under the RWA. One then recovers as the remaining terms the desired combination $-\tfrac{1}{2}\gamma_{\mathrm{n}}B_{\mathrm{1,I}}I_{\mathrm{y}}+\tfrac{1}{2}\gamma_{\mathrm{n}}B_{\mathrm{1,Q}}\cos\left(2\pi f\right)I_{\mathrm{x}}$. This implies that one may correct for the angle $\theta$ between quadrupole interaction and drive axis straightforwardly by including this angle as a phase shift in the periodic drive.

C.4 Classical chaos for rotating frame parameters

As mentioned in the main text, by working in the rotating frame new parameters regimes may be explored that are difficult to reach in the laboratory frame. In particular, the regime of quadrupolar interaction strength much larger than effective linear interaction strength may now be explored. Similarly, the regime of large periodic drive may be explored, as the technique introduced here allows for any periodic drive strength up to the effective linear interaction strength. We extended our analysis of the classical dynamics presented in the main text to explore these regimes as well. Not surprisingly, we find that the regime of large periodic drive strength allows for very large chaotic fractions of the classical phase space. Even in regimes where the linear interaction strength is an order of magnitude smaller compared to the quadratic interaction strength, we still find significant chaotic fractions. These findings are summarized in Fig. 5.

Appendix D Quantum dynamics and Floquet formalism

The evolution of time-dependent periodic Hamiltonians, such as the quantum driven top, can be transformed to time-independent evolutions using the Floquet formalism. The Floquet operator $\mathcal{F}$ is equal to the time evolution operator $\mathcal{U}$ over one full period $\tau$:

[TABLE]

As a result, the Floquet operator $\mathcal{F}$ has the property $\ket{\psi(\tau)}=\mathcal{F}\ket{\psi(0)}$, irrespective of the state $\ket{\psi(0)}$. A consequence is that $\ket{\psi(N\tau)}=\mathcal{F}^{N}\ket{\psi(0)}$, and so once $\mathcal{F}$ is known, any state can be straightforwardly evolved over a discrete number of periods by repeated application of $\mathcal{F}$.

The Floquet operator can be decomposed into eigenstates $\ket{\Phi_{i}}$ and corresponding eigenvalues $\lambda_{i}$, which all satisfy $\left|\lambda_{i}\right|=1$ ($\mathcal{F}$ is unitary). All eigenvalues are therefore of the form $\lambda_{i}=\exp{\left(-i\epsilon_{i}\tau/\hbar\right)}$, where the angular frequency $\epsilon_{i}$ is known as the quasienergy of the corresponding eigenstate. This enables decomposition of any initial state $\ket{\psi(0)}$ into the Floquet eigenstates, and straightforward calculation of its state after evolution of $N$ periods:

[TABLE]

Each Floquet eigenstate accumulates a phase determined by its respective quasienergy, and so superpositions of eigenstates will lead to interference effects. This offers an explanation for the evolution of spin coherent states (Fig. 4); whereas the spin coherent state with many Floquet components [orange (light gray)] has many interfering frequencies, the state that exhibits dynamical tunneling [blue (dark gray)] primarily consists of only two Floquet components. In the latter case, there are two main interfering frequencies, and their difference determines the dynamical tunneling frequency. The Floquet formalism also connects to the varying dephasing rates for different initial states (Fig. 3), as the individual Floquet quasienergies have a different degree of sensitivity to perturbations. More generally, the Floquet formalism clearly emphasizes the discrete nature of quantum mechanics, leading to quasiperiodicity that causes partial revivals of initially localized states, as opposed to the exponential divergence of trajectories in the case of classical chaos.

Appendix E Spin $I>1/2$: some properties and definitions

E.1 Spin coherent states

Spin coherent states, also known as Bloch states, form a subset of the possible states of a spin-$I$ system. The formulation and some relevant properties of the spin coherent states are reviewed here.

For a spin with spin quantum number $I$ and magnetic quantum number $m\in\left[-I,-I+1,...,I-1,I\right]$, states can be described by the basis states $\ket{I,m}$, which are eigenstates of the $I_{\mathrm{z}}$ spin operator with corresponding eigenvalue $m$.

We introduce spherical angles $\phi$ (azimuthal) and $\theta$ (polar), and an operator $R_{\mathrm{\theta,\phi}}$ corresponding to a rotation over an angle $\theta$ about an axis $\left(\sin\phi,-\cos\phi,0\right)$, given by

[TABLE]

The spin coherent states $\ket{\theta,\phi}$ are now defined as the state $\ket{I,I}$ rotated by $R_{\mathrm{\theta,\phi}}$ Arecchi et al. (1972):

[TABLE]

All spin coherent states share the property

[TABLE]

Furthermore, the spin coherent states form an overcomplete normalized basis, having nonzero overlap with each other. Only opposite spin coherent states with a $\pi$ difference in polar angle $\theta$ are orthogonal and have zero overlap, analogous to the orthogonality of $\ket{I,-I}$ and $\ket{I,I}$.

Rotated spin operators can be obtained by applying $R_{\mathrm{\theta,\phi}}$ to the original spin operators:

[TABLE]

in which case the spin coherent state $\ket{\theta,\phi}$ is an eigenstate of $I_{\mathrm{z^{\prime}}}$ with eigenvalue $I$. Spin coherent states are the only states for which the uncertainty relation ($\sigma_{I_{\mathrm{x^{\prime}}}}^{2}\sigma_{I_{\mathrm{y^{\prime}}}}^{2}\geq\tfrac{\hbar^{2}}{4}\braket{I_{\mathrm{z^{\prime}}}}^{2}$) becomes an equality, and are therefore also known as minimum-uncertainty states.

For $I=1/2$, all pure spin states are spin coherent states, leaving the concept unnecessary. However, this concept is very useful for $I>1/2$, since now the spin coherent states form a distinct subset of pure states for which the spin is maximally aligned in a certain direction $(\theta,\phi)$, with a minimum-uncertainty spread around it. As such, these states are the closest analog to classical angular momentum, and are therefore used in the quantum driven-top experiments as corresponding initial quantum states.

E.2 Husimi Q distribution

The Husimi Q distribution Husimi (1940) is a quasiprobability distribution that is used to represent quantum states. It is particularly useful here to provide a visualization of high-dimensional quantum states (Fig. 6). For a given density matrix $\rho$, the Husimi Q function is defined as

[TABLE]

where $\ket{\theta,\phi}$ is a spin coherent state. In the case of a pure state $\left(\rho=\ket{\psi}\bra{\psi}\right)$, the Husimi Q distribution simplifies to $Q(\theta,\phi)=\tfrac{1}{\pi}\left|\braket{\psi}{\theta,\phi}\right|^{2}$, the overlap-squared between $\ket{\psi}$ and $\ket{\theta,\phi}$.

The Husimi Q distribution satisfies certain properties required for a joint probability distribution, because the distribution is normalized and non-negative with values ranging between $0\leq Q(\theta,\phi)\leq 1/\pi$. However, since different spin coherent states are nonorthogonal, different coordinates $(\theta,\phi)$ do not represent distinct physical contingencies, and different values of $Q$ are not the probability of mutually exclusive states, a requirement for a joint probability distribution. This reflects the fact that quantum mechanics lacks a clear phase-space description, as opposed to classical mechanics. Although not providing a true phase-space description, quasiprobability distributions such as the Husimi Q distribution provide the closest quantum-mechanical proxy to it as an alternative and complete representation of a quantum state (invertible to the original density matrix representation).

The Husimi Q distribution is preferred here over the Wigner distribution as the phase space representation of quantum states because classical effects are emphasized. For instance, application of the Husimi Q function to a spin coherent state closely matches a classical point in phase space with the addition of an uncertainty spread (Fig. 7). A downside is that typical quantum phenomena such as interference are not clearly visible. To visualize these effects, other distributions such as the Wigner distribution are more attractive candidates.

E.3 Uncertainty and size of $I$

Spin coherent states are the set of minimum uncertainty states ($\sigma_{I_{\mathrm{x^{\prime}}}}^{2}\sigma_{I_{\mathrm{y^{\prime}}}}^{2}=\tfrac{\hbar^{2}}{4}\braket{I_{\mathrm{z^{\prime}}}}^{2}=\tfrac{\hbar^{2}}{4}I^{2}$). Its uncertainty defines a typical area over which a coherent state spreads in phase space, and should be compared to the total surface of the phase space ($4\pi I^{2}$). This leads to a concept of relative uncertainty $\sigma_{I_{\mathrm{x^{\prime}}}}\sigma_{I_{\mathrm{y^{\prime}}}}/4\pi I^{2}\propto 1/I$, which is a measure of how well a quantum state is localized in phase space. This underlines the importance of having large enough $I$, as illustrated by comparing the Husimi Q function of spin coherent states for different $I$ (Fig. 7).

Appendix F Dependence of dynamical tunneling rate on system parameters

We investigate how the dynamical tunneling rate is influenced by different system parameters. Dynamical tunneling arises naturally in the quantum system, since the uncertainty spread of a quantum state prevents it from being truly localized within a classical region of phase space. Therefore, a state prepared within one of the classically stable regions of the driven top will have a finite overlap of its wave function with the other classically stable region, and may tunnel back and forth between these two regions. Qualitatively, the dynamical tunneling rate can be understood as determined by the amount of such a wave function overlap.

To get a better understanding of the tunneling rate we have numerically studied its dependence on different system parameters. First, as the spin quantum number $I$ grows, the relative uncertainty spread of a spin coherent state on the sphere of radius $I$ shrinks, and the dynamical tunneling rate decreases accordingly [Fig. 8(a)]. Second, the parameters $B_{0}$ and $Q$ define the relative distance of the corresponding classically stable regions. Upon increasing $B_{0}$, the two stable regions come together, and the tunneling rate increases [Fig. 8(b)], whereas increasing $Q$ has the opposite effect.

Appendix G Quantum state purity vs fluctuations in $B_{0}$ and $B_{1}$

In addition to Fig. 3, where fluctuations in the parameter $Q$ are considered, fluctuations in the parameters $B_{0}$ and $B_{1}$ are studied here (Fig. 9). The approach used to obtain the results shown in Fig. 3 is repeated for this study and details of these simulations may be found in Appendix K. Similar behavior is found compared to considering fluctuations in $Q$. We note that in all cases the correspondence between classical phase space being regular or chaotic and quantum state purity decaying slow or fast is not exact. Rather the general structure of the underlying classical phase space is recovered; i.e., classically stable points correspond to quantum states being less sensitive to fluctuations.

Appendix H Nuclear quadrupole interaction

As pointed out in the main text, nuclei with a spin $I>1/2$ have a quadrupole moment $Q_{\mathrm{n}}$ due to a nonspherical charge distribution Kaufmann and Vianden (1979). This charge distribution has an axis of symmetry that aligns with the nuclear angular momentum and interacts with an electric field gradient (EFG). We estimate the strength of quadrupole interaction due to such an EFG for the different donors in silicon.

For group-V donors in silicon, the EFG is produced by external charges, such as the donor-bound electron or the crystal lattice. In a sample of unperturbed bulk silicon, the tetrahedral donor symmetry results in canceling EFG components and consequently a vanishing quadrupole interaction in the ground state. In a realistic micro- or nanoelectronic device, strain and/or electric fields can break the bulk symmetry of the donor wave functions Slichter (1990). Strain is typically present in devices as a result of the mismatch of the coefficient of thermal expansion of different materials, for example between the metal electrodes and semiconductor substrate Thorbeck and Zimmerman (2015); Pla et al. (2018). The subsequent rearrangement of the multivalley electron state generates electric field gradients which can couple to the electric quadrupole moment of the donor Pla et al. (2018); Franke et al. (2016). Strain also acts on the silicon crystal to produce an EFG through the gradient elastic tensor $\mathbf{S}$ Sundfors (1979), which results in a quadrupole interaction even in the absence of the donor-bound electron Franke et al. (2015).

H.1 Estimates of Nuclear Quadrupole Interaction

We define the quadrupole interaction strength as the factor that precedes the quadratic $I_{\mathrm{z^{\prime}}}$ term in Eq. (3) of the main text,

[TABLE]

Here $\gamma_{\mathrm{s}}$ is the Sternheimer antishielding factor, $e$ is the elementary charge, $Q_{\mathrm{n}}$ is the nuclear quadrupole moment, $V_{\mathrm{ij}}=\tfrac{\partial^{2}V}{\partial{i}\partial{j}}$, $\left(i,j\in x^{\prime},y^{\prime},z^{\prime}\right)$, are the partial second derivatives of electric potential $V$, and $h$ is Planck’s constant. This differs from the more conventional definition of the quadrupole interaction (by including multiplication factor $3/4I\left(2I-1\right)$), as this allows a direct comparison to the equivalent classical parameter $\beta$. Accurate calculation of $Q$ for donors is complicated by the multiplicative term $\gamma_{\mathrm{s}}$. This factor relates to the Sternheimer antishielding effect, a phenomenon that describes the rearrangement of the inner electron shells in response to an external EFG, effectively enhancing the EFG experienced by the nucleus Kaufmann and Vianden (1979). The Sternheimer antishielding factor $\gamma_{\mathrm{s}}$ can be considerable; theoretical calculations Feiock and Johnson (1969) for isolated As and Bi ions show an enhancement of about one order of magnitude for As and up to three orders of magnitude for Bi. To the best of our knowledge, no such calculations have been completed for Sb. Furthermore, it is unknown how the covalent bonding of the donor to the silicon lattice affects $\gamma_{\mathrm{s}}$. As a result of the uncertainty in $\gamma_{\mathrm{s}}$, it is difficult to make purely theoretical predictions of $Q$ for donors in silicon.

Recent experiments Franke et al. (2015); Mortemousque et al. (2016); Franke et al. (2016); Pla et al. (2018) on quadrupole effects in silicon devices have produced some quantitative results that can be used to estimate $Q$ for As, Sb and Bi donors in Si. We will present an analysis for each of the donors in sections below, predicting the quadrupole coupling in the ionized charge state D+, where the EFG is produced by the crystal lattice alone.

H.1.1 Arsenic

While arsenic has the lowest nuclear spin of the donors considered here – making it less suitable for comparison with classical dynamics – it is a relatively well-studied donor for its quadrupole properties. In Refs. 62 and 63, spectroscopy of As donors in a strained silicon sample has been performed (uni-axial strain $\epsilon_{\perp}\approx 3\times 10^{-4}$). For an ionized donor, the EFG is generated through the gradient elastic tensor $V=\mathbf{S}\cdot\mathbf{\epsilon}$ (where $\mathbf{\epsilon}$ is the strain tensor in Voigt notation), implying a linear relationship between the applied strain and quadrupole interaction. Measurement of the quadrupole shifts in two samples of different surface planes [(100) and (111)] enabled the extraction of the nontrivial gradient elastic tensor components, $S_{11}=$1.5\text{\times}{10}^{22}\text{,}\mathrm{V}\text{,}{\mathrm{m}}^{-2} and $S_{44}=$6.8\text{\times}{10}^{22}\text{\,}\mathrm{V}\text{\,}{\mathrm{m}}^{-2}. These components (which include the Sternheimer antishielding factor $\gamma_{\mathrm{s}}\approx-7$ for As Feiock and Johnson (1969)) can be used to provide a rough estimation of the D+ quadrupole coupling for Bi and Sb (see below). In a nanodevice, strains of order $10^{-3}$ are expected directly underneath the metallic surface electrodes Thorbeck and Zimmerman (2015) due to the mismatch in thermal expansion coefficients, similar in magnitude to those observed in the prestrained devices in Refs. 62 and 63. This allows us to estimate the quadrupole interaction strength in a device: $Q\approx$210\text{,}\mathrm{kHz} for a (111) surface and $Q=$60\text{\,}\mathrm{kHz} for a (100) surface.

H.1.2 Bismuth

In order to estimate $Q$ in the D+ state for Bi, we require the antishielding factor $\gamma_{\mathrm{s}}$. We take an order-of-magnitude estimate, only serving as a rough guide, of $\gamma_{s}\approx 100$ for the antishielding factor of Bi. This value is based on the simulations of measured data reported in Ref. 64.

Using the estimated magnitude of $\gamma_{s}$ for Bi, the measured gradient elastic tensor components $S_{11}$ and $S_{44}$ of As (converted to strain using the theoretical magnitude of $|\gamma_{s}|\approx 7$ for As), we estimate the quadrupole interaction strength achievable in a nanodevice for the D+ state to be $Q\approx$240\text{,}\mathrm{kHz} for a (111) silicon surface and $Q=$45\text{\,}\mathrm{kHz} for a (100) surface.

H.1.3 Antimony

Antimony is the least understood of the group-V donors for its quadrupole properties. There are no theoretical calculations for the Sternheimer antishielding factor, and no experimental data on the interplay between strain and quadrupole interaction. In the EDMR experiments of Ref. 64, $Q$ was found for the neutral 121Sb donor to be approximately half that of the measured value for 209B, with the caveat that the implantation conditions were different for the Sb and Bi samples, and no estimate was given for the likely separation between the donors and the readout centers (which strongly influences the EFG) in the Sb sample.

Appendix I System characterization via NMR spectroscopy

High-precision measurements of the nuclear-spin Hamiltonian parameters are needed to achieve an accurate comparison between the quantum system and its classical equivalent. Furthermore, this is a crucial requirement in enabling arbitrary state preparation as discussed later (Appendix J). This section describes how the different parameters can be extracted through NMR spectroscopy of the ionized nucleus.

I.1 NMR spectroscopy

Starting from a system where none of the NMR transition frequencies are known, the first step is to find the ESR frequency, which uniquely depends on the (unknown) nuclear spin state through the hyperfine coupling. To this end, the electron spin is initialized in the down state through spin-dependent tunneling from the SET island onto the ionized donor, after which a voltage is applied to ensure that the electron is unable to tunnel back. A pulse with linearly increasing frequency is applied that adiabatically inverts the electron state if its resonance frequency lies within the frequency range Laucht et al. (2014). The bias voltage is then modified to ensure that the electron can only tunnel back onto the SET island if its spin state has successfully been flipped, which is measured as a finite SET current.

Once the first ESR frequency is known, an NMR frequency of the nuclear spin state can be found. The measurement sequence is nearly identical, with the exception that before loading an electron, a voltage is applied to force the donor into the ionized state, followed by another pulse with linearly increasing frequency in the expected NMR frequency range. If this pulse sweeps over one of the active NMR frequencies (which have an unique value due to the quadrupole interaction), this will adiabatically invert the nuclear spin state. As a result, the ESR frequency will change, and so no subsequent adiabatic inversion of the electron spin state will occur. Repeating the measurement flips the nuclear spin state, resulting in an alternating electron spin flip probability, allowing determination of the NMR frequency. The previous two steps can then be repeated to iteratively determine all ESR and NMR frequencies.

The above scheme assumes that loading or unloading an electron leaves the nuclear-spin eigenstates intact, which may require an adiabatic electron-loading scheme (Appendix J).

I.2 Extraction of Hamiltonian parameters

Combining NMR spectroscopy with full control of both the direction and strength of $B_{0}$ allows extraction of all Hamiltonian parameters. This relies on the underlying assumption that each of the NMR transitions is individually addressable, a condition that can be satisfied by assuming a $Q$ larger than the NMR transition linewidth [Fig. 10(a)]. It eases analysis to operate in the large magnetic field limit ($Q\ll\gamma_{\mathrm{n}}B_{0}$), where the quadrupole interaction can be treated as a perturbation to Eq. (2) of the main text [Fig. 10(b)].

The first goal is to determine the orientation of the quadrupole coordinate system $(\hat{x}^{\prime},\hat{y}^{\prime},\hat{z}^{\prime})$ [Eq. (3) main text]. Once known, aligning the $B_{0}$-axis with $\hat{z}^{\prime}$ provides the strength of $Q$, since successive NMR transitions have a constant separation of $2Q$ [Fig. 10(c)]. This equidistant spacing is lifted by a nonzero asymmetry parameter ($\eta>0$), but is recovered in the high magnetic-field limit [Fig. 10(d)].

The quadrupole’s primary axis $\hat{z}^{\prime}$ can be found through successive spectroscopy measurements while rotating $B_{0}$. Without prior knowledge of the quadrupole’s coordinate system, $B_{0}$ is initially rotated around an arbitrary axis. It is guaranteed that such a spectroscopy will reveal two symmetry axes, one of which is perpendicular to $\hat{z}^{\prime}$. A second rotation of $B_{0}$ around this particular symmetry axis will align $B_{0}$ with $\hat{z}^{\prime}$ at some specific angle. This point has maximum and (nearly) equidistant separation of the spectral lines, thereby revealing the quadrupole’s primary axis $\hat{z}^{\prime}$ [Fig. 10(e)].

The two remaining unknown parameters of the quadrupole interaction, the asymmetry $\eta$ and its orientation, can be found through a final rotation of $B_{0}$ around the $\hat{z}^{\prime}$ axis. If $\eta=0$, spectral lines will be independent of this rotation, while for $\eta>0$ again two symmetry axes will be revealed. The symmetry axes with the largest separation of spectral lines corresponds to $B_{0}$ being parallel to the secondary $\hat{y}^{\prime}$ axis, and the strength of spectral line variation is determined by the size of $\eta$.

The sequence sketched here allows variation of a single experimental handle (NMR frequency, $B_{0}$ direction) to isolate the effect of each parameter. This allows an accurate determination of all the relevant Hamiltonian parameters and a detailed understanding of the system.

Appendix J State preparation and measurement

J.1 Arbitrary state preparation

The ability to create an arbitrary target state $\ket{\psi_{T}}$ requires addressability of individual transition frequencies and accurate knowledge of all Hamiltonian parameters. The procedure described in this section (see Fig. 11) can be used for arbitrary state preparation, provided that the above conditions are fulfilled.

Assuming the system to be initialized in the ground state $\ket{e_{1}}$, the necessary sequence of pulses that results in the target state $\ket{\psi_{T}}$ can be found by solving the problem in reverse: go backwards in time and find the pulse sequence to end up in $\ket{e_{1}}$ starting from $\ket{\psi_{T}}$. Starting at $t=$0\text{,}\mathrm{s}$$ with initial state $\ket{\psi_{T}}=\sum_{k}a_{k}\ket{e_{k}}$, where $a_{k}=\braket{e_{k}}{\psi_{T}}$, pulses are iteratively chosen that transfer population from the populated highest-energy eigenstate to a lower-energy eigenstate. The procedure is as follows:

Choose eigenstate $\ket{e_{k}}$ with highest eigenvalue $\lambda_{k}$ and nonzero population $\left|a_{k}\right|$. 2. 2.

Find eigenstate $\ket{e_{k^{\prime}}}$ with lower eigenvalue $\lambda_{k^{\prime}}$ that has the highest coupling $\left|\braket{e_{k^{\prime}}}{I_{y}}{e_{k}}\right|$. 3. 3.

Transfer the population from $\ket{e_{k}}$ to $\ket{e_{k^{\prime}}}$ using a pulse with the following properties:

• •

frequency $f_{k,k^{\prime}}$, which is the transition frequency between $\ket{e_{k}}$ and $\ket{e_{k^{\prime}}}$;

• •

duration $t_{p}=\frac{2}{\Omega_{k,k^{\prime}}}\arctan{\frac{|a_{k}|}{|a_{k}^{\prime}|}}$, where $\Omega_{k,k^{\prime}}$ is the Rabi frequency;

• •

phase is the relative phase difference between $a_{k}$ and $a_{k^{\prime}}$, minus the phase offset $2\pi f_{k,k^{\prime}}(t-t_{p})$. 4. 4.

Update $t\rightarrow t-t_{p}$ and repeat steps until all population is transferred to $\ket{e_{1}}$.

These steps can be followed to create an arbitrary pure state, and in particular the spin coherent states $\ket{\theta,\phi}$ used in the proposed quantum driven-top experiments. While we do not expect our state-preparation fidelity to be limited by effects such as decoherence (state preparation takes $<$1\text{,}\mathrm{ms}, $T_{2}^{*}\sim$1\text{\,}\mathrm{s}), the total pulse-sequence duration can be further reduced by using multi frequency pulses.

J.2 Measurement

The presence of a strong quadrupole interaction with axis perpendicular to that of the linear interaction inhibits defining a clear quantization axis. Upon adding an electron to the ionized donor, the accompanying hyperfine interaction can be approximated as an enhancement of the linear interaction by about an order of magnitude for ${}^{123}\mathrm{Sb}$ (two orders for ${}^{209}\mathrm{Bi}$ ), significantly altering the eigenbasis of the donor. The existing techniques for initialization and readout of the nucleus rely on electron tunneling events between donor and SET island. Here, these tunneling events are accompanied by a change of eigenbasis, resulting in a probabilistically modified state after every such event. Modified protocols may therefore be necessary for high-fidelity nuclear state initialization and readout.

The proposed solution is the adiabatic transferring of the electron from donor to the silicon/silicon-dioxide interface, which maps neutral-donor eigenstates to ionized-donor eigenstates and vice versa. This can be achieved by the addition of an electrostatic gate above the donor that can attract its outer electron, a technique that is currently being developed in the context of achieving electrically active transitions and long-range coupling of donor nuclear spin qubits Tosi et al. (2017). Once the electron is moved from the donor to the interface, spin-dependent tunneling of the electron from the interface to the SET island allows for readout without affecting the nucleus.

Appendix K Numerical methods and simulation details

Both the classical and quantum simulations were primarily performed using the commercial software package MATLAB Release 2016b, The MathWorks, Inc., Natick, Massachusetts, United States. In this section we will describe the simulation techniques for both the classical and the quantum simulations. The source code of our simulations is available online in the Supplemental Material sup .

K.1 Classical simulations

The most computationally expensive simulations were those to determine the percentage of phase space that is chaotic [Fig. 1(h), 1(i)). Both color maps consist of $25\times 25$ logarithmically spaced points, each of which corresponds to a particular parameter set. For each parameter set, a total of 2000 initial angular-momentum coordinates were chosen uniformly distributed over the phase space. To determine whether the dynamics of an initial coordinate is chaotic, a neighboring point with distance $10^{-8}$ was chosen, and both were evolved for a fixed duration of $100/\alpha$. Whether or not a trajectory is chaotic is determined by measuring the distance between the two points over time, and fitting this to an exponential curve. Chaos is characterized by an exponential sensitivity to perturbations, and so the trajectory is categorized as chaotic if its exponent is above a certain threshold. This procedure is repeated for each of the 2000 initial conditions, resulting in the percentage of phase space that is chaotic.

Some trajectories, especially near a chaotic-regular boundary, can display an initial exponential divergence but nevertheless behave regularly over sufficient evolution time. These cases, although uncommon, can result in the trajectory being wrongly categorized as chaotic, and we expect a small uncertainty in the percentages of Fig. 1(h), 1(i). It should furthermore be noted that the fitted exponent is not necessarily equal to the Lyapunov exponent, as the intertrajectory distance may have increased sufficiently to be limited by the finite size of the phase space. Although the exponent is then an underestimate of the Lyapunov exponent, it will certainly be above the chaotic threshold, thereby correctly characterizing the trajectory as chaotic.

The classical simulation results shown in Fig. 1 were obtained using the ordinary-differential-equation (ODE) solver SUNDIALS Hindmarsh et al. (2005). Communication between MATLAB and Sundials was through self-written C code that was optimized for the driven-top system. The combination of SUNDIALS and the intermediate C code resulted in a computational speedup of over an order of magnitude compared to the native MATLAB ODE solvers. As chaotic dynamics are highly sensitive to perturbations, stringent error tolerances were chosen to ensure a high degree of accuracy in the computation of the trajectories.

K.2 Quantum simulations

The quantum driven-top system is evolved using the Floquet operator $F$, which can be approximated through segmentation as

[TABLE]

which becomes an equality in the limit $N\rightarrow\infty$. In the simulations, a fixed value of $N=1000$ is used, as results showed that the Floquet operator did not significantly change upon further increasing $N$. Additionally, the SUNDIALS ODE solver was used to compute $F$, and was found to be nearly identical to $F$ computed using the above method.

Decoherence of spin coherent states under influence of the driven-top Hamiltonian [Eq. (4)] was simulated by randomly fluctuating a Hamiltonian parameter during its evolution and averaging over many such evolutions. To this end, Floquet operators were calculated for 30 values of $Q$ uniformly distributed within three standard deviations of the its mean value ($Q=800\pm 4$ $\mathrm{kHz}$). For the evolution, a sequence of Floquet operators was chosen through random sampling of this set using a Gaussian distribution. This sequence was then applied to all initial spin coherent states, and this process was repeated for 300 such sequences. For each spin coherent state, the final density matrices were averaged, resulting in a mixed state $\rho$, from which the purity $\operatorname{Tr}(\rho^{2})$ was determined.