{'fraction_non_alphanumeric': 0.03408677800825439, 'fraction_numerical': 0.013878493680901568, 'mean_word_length': 4.366536607757632, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Ultra long spin decoherence times in graphene quantum dots with a small number of nuclear spins

We study the dynamics of an electron spin in a graphene quantum dot, which is interacting with a bath of less than ten nuclear spins via the anisotropic hyperfine interaction. Due to substantial progress in the fabrication of graphene quantum dots, the consideration of such a small number of nuclear spins is experimentally relevant. This choice allows us to use exact diagonalization to calculate the longtime average of the electron spin as well as its decoherence time. We investigate the dependence of spin observables on the initial states of nuclear spins and on the position of nuclear spins in the quantum dot. Moreover, we analyze the effects of the anisotropy of the hyperfine interaction for different orientations of the spin quantization axis with respect to the graphene plane. Interestingly, we then predict remarkable long decoherence times of more than 10ms in the limit of few nuclear spins.


I. INTRODUCTION
In recent years, spin qubits hosted in solid state nanostructures have been under extensive research due to their possible applications in quantum information processing and computation. Among the host materials of spin qubits, quite different approaches can be found, for instance, III-V-semiconductor and carbon nanotube quantum dots 1 (QD) as well as nitrogen vacancies in diamond 2 . These host materials show promising prospects but, unfortunately, also come with certain drawbacks.
A precise control of the qubit state is the major advantage of III-V-semiconductor QDs based on Al(Ga)As heterostructures. Preparation and readout of the qubit with high fidelity via electrostatical gates has been demonstrated in many ground-breaking experiments [3][4][5][6][7][8][9][10][11][12][13] . However, the disadvantage of this material system is the presence of many nuclear spins inherent to the atoms of group III and group V elements of the periodic table. These nuclear spins give rise to a fast decoherence of the electron spin 14,15 . Elaborate design of experiments including pulse sequences and methods to polarize the nuclear spins [16][17][18][19][20][21] such as dynamic nuclear polarization may help to compensate the effect of the spin bath. Nevertheless, it seems desirable to reduce the number of nuclear spins which can be achieved on the basis of other host materials.
Obvious candidates are carbon and silicon, since their spin carrying isotopes have only very low natural abundances of about 1% and 5%, respectively. In silicon, the qubits can be fabricated 22,23 either using donor impurities or by confining a single electron via electrostatical gates. However, a controlled localization of the donor impurities is still a challenging task and electrostatically confined Si QDs often involve nanostructures with other materials like Ge, which potentially introduce addi-tional nuclear spins 23 . Carbon based QDs can be realized by confining an electron spin in carbon nanotubes [24][25][26][27][28][29][30][31] (CNT) via electrical gates allowing for a good control in the few electron regime. However, the curvature of the CNTs gives rise to a sizable spin-orbit coupling, yet another intrinsic source of decoherence to the electron spin. A different approach to a carbon based QD is the use of nitrogen vacancies in diamonds 2,32-35 , which show tremendously long coherence times. Unfortunately, control and readout of the qubit have to be done optically, which is disadvantageous for the realization of future onchip electric circuits.
These examples illustrate a more general issue of designing qubits, where an easy (electric) control and scalability seem to compete with noiseless environments and, hence, long decoherence times. A system potentially providing the best of both worlds is a graphene QD 36,37 , which offers very interesting electronic properties 38 and a small spin-orbit coupling [39][40][41][42] , as well as the possibility to control the number of nuclear spins by isotopic purification 25,43,44 . Moreover, the hyperfine interaction between the remaining nuclear spins and the electron spin is much smaller than in GaAs or Si. Additionally, the hyperfine interaction in graphene is anisotropic which could provide interesting applications as we discuss at the end of this article.
Experimentally, graphene QDs are, for instance, realized by confining electrons with gates in bilayer graphene 45,46 and graphene nanoribbons 47,48 , respectively, or by etching the QD structure out of graphene flakes [49][50][51][52][53][54][55][56][57][58][59] . Typical diameters are of the order of tens to hundreds of nanometers resulting in K = 15 to K = 1500 nuclear spins assuming a natural abundance of spin carrying 13 C of 1%. Thus reducing the abundance of 13 C by only two orders of magnitude leads to very small spin baths even in the case of rather large QDs. Recently, ultra small graphene QDs with diameters in the 1 nm range were made by electroburning 60 . Altogether, these considerations show that the study of few nuclear spin models with K < 10 as considered in this work is highly relevant for future research in the field.
In this paper, we aim to set the basis for forthcoming investigations of the spin dynamics in graphene nanostructures. Besides quantum information theory, especially ongoing research on magnetism on edges [61][62][63][64][65][66] and vacancies 67,68 in graphene can benefit from a detailed knowledge of the properties of the anisotropic hyperfine interaction (AHI). Moreover, we intend to complement our previous analytic study 69 of the electron spin dynamics. Considering a large nuclear spin bath, we investigated the coherence of the electron spin in a non-Markovian approach using a generalized master equation. In this work, however, we were limited to large external magnetic fields in order to justify the perturbative treatment of the hyperfine interaction.
Since we restrict ourselves to less than ten nuclear spins in the present work, we can apply exact diagonalization to the hyperfine Hamiltonian, which offers a powerful tool to investigate the dynamics of the electron spin for a wide parameter regime 14,[70][71][72][73][74][75][76][77] . In particular, we analyze the role of the number of nuclear spins K, their position within the QD, as well as their initial spin state. Thereby, we use the long-time average S z T of the longitudinal electron spin and its decoherence time T D to quantify the influence of these different aspects. Moreover, we investigate the dependence of S z T and T D on the orientation of the spin quantization axis with respect to the graphene plane. For the long-time average, we find a continuous crossover from a initial state dominated regime for K < 5 to a regime more affected by the configuration of the nuclear spins for K > 6 where the relative positions of the nuclear spins with respect to each other matter. As we will show below, this behavior can be understood by an analysis of the Hilbert space dimensions as well as of the spatial distribution of the nuclear spins in the QD. Besides this regime change, a growing nuclear spin bath suppresses fluctuations around the long-time average more and more effectively, while the average itself is almost constant for all K < 9 with S z T ≈ /4 in the out-ofplane orientation and S z T ≈ 0 in the in-plane case. By resolving the orientation dependence in more detail, we find good agreement with S z T (β) = S z T (0) · cos(β) 2 , where β = 0 and β = π/2 correspond to the out-of-plane and in-plane orientation, respectively, cf. Fig. 1 (a).
Evidently, the decoherence time T D strongly depends on the number of nuclear spins K. We observe that the configuration of the nuclear spins is decisive even for very small numbers of nuclear spins, whereas the initial states play only a minor role. Depending on the relative positions of the nuclei, the decoherence times may deviate over several orders of magnitude. This behavior can be traced back to changes of the spectrum of eigenvalues of the full Hamiltonian. Moreover, the decoherence times significantly differ between out-of-plane and in-plane orientation. For K = 3 and β = 0, the majority of inves-tigated configurations show very long decoherence times above 10 ms where, in many cases, even no decoherence at all was found. For β = π/2, in contrast, we always find decoherence, which predominantly occurs within 500 µs. With increasing bath size, the decoherence times decrease for both orientations of the quantization axis. Then, decoherence times below 500 µs are most common.
The article is organized as follows. In Sec. II, we explain our model of the QD and discuss all relevant interactions of the spins with each other. Subsequently, in Sec. III, we present the method used to obtain both the long-time average of the electron spin and its decoherence times. All results are shown and analyzed in Sec. IV. Based on a summary in Sec. V, we give an outlook on possible applications of few nuclear spin graphene QDs and on interesting future projects in this field.

II. MODEL
We study the spin dynamics in a graphene quantum dot, where one electron spin is in contact with a bath of nuclear spins hosted by the 13 C atoms. Due to the confinement, the electron can occupy a discrete spectrum of bound states, with an energy splitting between different states 36,[78][79][80][81][82] . If the temperature is small compared to the level spacing ∆E of these bound-state energies, the electron resides in the ground state, which we describe by an envelope function φ( r). Hence, the probability to find the electron in a certain region of the QD can be described by its absolute square |φ( r)| 2 . In this paper we define the "center" of the dot as the region around r max , where the envelope function is maximal |φ( r max )| 2 . Far away from this center, the envelope function has to vanish: In this work, we model a graphene QD by the set of where C = 10 −6 is a constant cut-off. A plot of a QD realized in this way is shown in Fig. 1. The choice of this finite system of discrete sites r i imposes a normalization condition since we want to find the electron with probability 1 somewhere within the dot. Effectively, we ignore everything outside the barrier defined by the cut-off, which is justified by the vanishing probability to find the electron there. This choice will become more clear when we discuss the hyperfine interaction between the electron and a single nuclear spin below.
A possible choice for the envelope function is a Gaussian where r = | r| is the absolute value of the electron position and the norm φ 0 is chosen to satisfy the normalization condition in Eq. (3). This assumption is also in agreement with a recent experiment investigating the wave function of a graphene QD with soft confinement 83 . Note that the envelope function is not the exact electron wave function, but should give a good approximation to the precise solution. This can be seen, for instance, in graphene QDs based on semi-conducting armchair nanoribbons 36 . The most important aspects, which are captured by this specific choice, are the absence of nodes in the ground state, a peak of the wave function in the center as well as a strong decay at the edges of the dot, which is illustrated by Fig. 1 (c). Having defined the shape of our dot, we are now able to introduce the nuclear spins present in the system. Since we do not have further knowledge about the distribution of 13 C within the dot, we randomly place the nuclear spins on the sites defined by Eq. (2), where each site is chosen with equal probability. An example of a configuration of ten nuclear spins is shown in Fig. 1.
We now proceed to the interactions between the nuclear spins and the electron spin and between themselves. The most relevant spin-spin interaction in our system is the hyperfine interaction between the electron spin S and K nuclear spins I k located at sites r k , where the indices µ and ν run over spatial coordinates x, y, z. The energy scale of this interaction is given by A HI = 0.6 µeV and ← → A µν is a spherical tensor 84 of rank 2, which takes into account the anisotropy of the hyperfine interaction in graphene 85 . Remarkably, this interaction is strongly modulated by the envelope function, a fact which arises from the on-site nature of the hyperfine interaction. Thus, making the boundary of the dot smooth by taking a very small cutoff C → 0 in Eq. (2), would only add vanishingly small contributions to the interaction in Eq. (5). Therefore we choose a small, but finite cutoff for simplicity.
Besides the anisotropic hyperfine interaction (AHI), there is also a dipole-dipole interaction between pairs of nuclear spins. In the parameter regime considered in this work, this interaction, however, is about five orders of magnitude smaller than the AHI and, thus, neglected. We also proved its irrelevance by a numerical study, which we will not present here.
Since external magnetic fields allow to manipulate spins experimentally, we include a Zeeman-Hamiltonian to account for this: where we used the fact that the electron gyromagnetic ratio γ S = 1.76 × 10 11 s −1 T −1 is much larger than the gyromagnetic ratio of the nuclear spins γ13 C = 6.73 × 10 7 s −1 T −1 to justify the right-hand side of Eq. (6).
In the presence of an external magnetic field, an interplay of the spin orbit coupling with acoustic phonons can lead to spin relaxation times T 1 ranging from milliseconds to seconds [86][87][88][89] for small external magnetic fields. The exact value of T 1 significantly varies with the spectrum of the phonons which depends on the details of the dot nanostructure. Providing that the graphene flake is flat throughout the spatial extent of the QD, however, the spin orbit interaction should be small. Thus, this assumption justifies to neglect the influence of the spin orbit interaction on our problem.
In the following, we aim to simulate a model experiment consisting of a preparation of the spins and the actual measurement of the spin dynamics. For the preparation, one can think of two different scenarios. First, the states of both the electron spin and the nuclear spins can be prepared in the presence of a strong external magnetic field B 0 ≫ K · A HI |φ(r max )| 2 / γ S = K · |φ(r max )| 2 · 5.7 mT, which imprints a well defined quantization axis. At the beginning of the actual measurement, this field is turned off or reduced to a finite value and the time evolution of the quantity of interest is recorded. However, since the tuning of magnetic fields is typically slow, this preparation scheme might not be adequate for a real experiment and we have to look for other solutions. In this case, we can think of injecting a spin polarized electron via a spin dependent tunneling process from a normal lead or via a spin conserving tunneling process from a spin polarized lead.
Anyhow, in both considered scenarios our system features two natural reference frames, the first defined by the graphene plane and the second one by the quantization axis of the electron as depicted in Fig. 1 (a). In order to clarify the notation, the graphene coordinate system (GCS) is written as ṽ = (ṽ x1 ,ṽ x2 ,ṽ x3 ) with all objects being marked with a tilde, whereas a vector in the quantization axis coordinate system (QCS) is labeled by Fig. 1. If one chooses, without loss of generality, the x 2 -and the y-axis to coincide, both coordinate systems are connected via a rotationD(β) by an angle β around this common axis.
Due to the symmetries of the carbon orbitals 85 , the spherical tensor ← → A µν of the AHI in Eq. (5) takes its simplest form in the GCS, namely while the spin operators,Ŝ µ andÎ k,ν , and the spin states are most conveniently defined in the QCS.
In the main part of this work, we are interested in the time-dependent expectation value of an arbitrary opera-torÔ where |ψ 0 describes the initial state of the total spin system andÛ (t) = exp[−i −1 tĤ] is the time evolution operator determined by the total HamiltonianĤ = H HI +Ĥ ZE . Note, that this is the total Hamiltonian with respect to the QCS, where the Zeeman Hamiltonian is always diagonal and the AHI Hamiltonian is obtained from its simple GCS formĤ HI bŷ Likewise, we could keep the Hamiltonian fixed in its GCS form and instead transform the operatorsÔ andĤ ZE as well as the initial state |ψ 0 from the QCS to the GCS. For technical reasons, however, we choose to transform the Hamiltonian, while keeping the operator and the initial state fixed for arbitrary β.

III. METHOD
In order to numerically compute the time evolution in Eq. (8), we need a basis to represent the state of our system and the operators acting on it. A natural choice for N = 1 + K spins is given by the tensor product states of the electron spin and the nuclear spin eigenstates where the electron spin is represented by |m n S , m n S =⇓ , ⇑ and the nuclear spin states by |m n k , m n k =↓, ↑. For convenience, we have ordered the nuclear spins m n S m n K m n K−1 . . . according to the value of the envelope function at the corresponding site: Within this basis, an arbitrary state is given by a linear superposition of these 2 N states with complex coefficients α n , while all operators are represented by 2 N × 2 N matrices. By diagonalizing the total HamiltonianMĤM † = diag(λ 0 , λ 1 , . . . , λ 2 N −1 ) we are able to re-express the time evolution operator whereM is an unitary operator formed by the eigenvectors ofĤ and the λ n are the corresponding eigenvalues. Finally, we re-write Eq. (8) and find which we will evaluate for different parameter regimes in the following section. The numerical diagonalization is performed using the EIGEN 90 package for C++. As one can notice from the explanations above, we deal with a quite big parameter space in which we can analyze the outcome of Eq. (14). First, we control the shape of the dot by means of the envelope function |φ(r)| 2 , secondly the number of nuclear spins K is variable and finally these spins can have different positions or configurations C within the dot. All of these parameters change the AHI Hamiltonian in Eq. (5). Moreover, we will investigate different initial states |ψ 0 of the electron and the nuclear spins affecting Eq. (14). Additionally, the eigenvector matrixM , appearing in this equation is a function of the twisting angle β between the GCS and the QCS. Note that the spectrum of eigenvalues λ n is unaffected by a change of β. Finally, we can also modify the absolute value of the external magnetic field, which we will parametrize by the resulting Zeeman energy of the electron A ZE .

IV. RESULTS
In this section, we present our findings for the model system defined above. All calculations were carried out using an envelope function of the Gaussian type in Eq. (4) with R = 7 a N N and a cut-off C = 10 −6 . This corresponds to a dot with diameter D ≈ 7.2 nm containing N sites ≈ 10 3 carbon atoms, such that K = 9 atoms correspond to the natural abundance n I = 0.01 of 13 C. In order to investigate the impact of different initial states, we choose random complex (RC) initial states 14,70 . These states were created by drawing complex coefficients α n from Re[α n ], Im[α n ] ∈ [−1, 1] with equal probability and normalizing them according to Eq. (12). Moreover, we choose the electron spin always to point down resulting in initial states consisting of |⇓ . . . states only, which means that α n = 0 for n ≥ 2 K .
In order to determine qualitatively and quantitatively the impacts of the parameters, we investigate the time dependent expectation value S z (t) of the longitudinal electron spin component, which is calculated using Eq. (14). A typical time evolution of S z (t) is plotted in Fig.  2. Within the decoherence time T D , the initial amplitude of the electron spin of /2 decays to its long-time average value, where still finite oscillations and beatings occur. This can be traced back to the finite size of the spin bath considered here.
Its long-time average value is calculated by where we average over N T = (T max − T min )/∆T time steps of width ∆T . In order to investigate the oscillations of S z (t) quantitatively, we consider the standard deviation as well as the sample range .
(17) The latter quantity is a measure for the occurrence of oscillations with a big amplitude which originate from either recurrences of the signal, beatings or an entire lack of decoherence. While for beatings one expects rather small sample ranges ∆S z < S z (0), the former two cases should give values on the order of the initial amplitude, ∆S z ∼ O(S z (0)) in the out-of-plane case β = 0 and ∆S z ∼ 2 · O(S z (0)) in the in-plane case β = π/2.
Besides these quantities characterizing the long time average of the electron spin, we are also interested in the amount of time it takes to decohere the system. In order to be independent from specific models of the decay, such as exponential or power-law decoherence, and to account for the characteristics of the numerics, we find this decoherence time T D by the first minimum exceeding a certain threshold C S . For clarity, C S is also illustrated in Fig. 2. This approach is similar to the one used in Ref. 77 to find the decoherence times. Of course, the choice of this constant C S changes the value of T D . However, its order of magnitude and its dependence on the different parameters is rather independent from a specific choice as long as C S is not too close to S z T , which we confirmed for different values of C S .
In the following, we analyze both the decoherence time and the long-time average of the longitudinal electron spin component for different parameter sets. For each number K of nuclear spins many initial states and configurations are created and labeled by numbers 0, 1, 2 . . . for later comparison of the results. Note, that for different nuclear spin numbers K these labels describe different initial states and configurations. Moreover, we concentrate on two orientations of the quantization axis, namely out-of-plane orientation for β = 0 and in-plane orientation with β = π 2 . We investigated the effect of finite magnetic fields for exemplary initial states, configurations and K = 2, 4, 6 nuclear spins, where we varied the resulting Zeeman constant from A ZE /A HI ≪ 1 to A ZE /A HI ≫ 1. For increasing A ZE , we find a continuous crossover to a perfect alignment of the electron spin in the case of a very strong magnetic field. In the following, we put A ZE = 0 because we would like to better understand the low-magnetic field behavior of the spin dynamics in the presence of the AHI.
A. Dependence of the long-time average on different initial states, configurations, and the number of nuclear spins First, we investigate the consequences of both different RC initial states and different configurations of the nuclei within the dot. We calculated S z T , σ Sz and, ∆S z for different parameter sets and found stable results for T min = 0.5 × 10 9 τ HI , T max = 1.5 × 10 9 τ HI , and  Fig. 3 but for inplane orientation β = π/2. For all parameters the long-time average of the longitudinal electron spin component sharply saturates at Sz T ≈ 0. In contrast to the out-of-plane case, the results seem independent of the initial state even in the few spin regime. For certain configurations, however, we find a non negligible dependence on the initial state. ∆T = 10 4 τ HI with τ HI = /A HI ≈ 1 ns. In Fig. 3, we plot the long-time average S z T as a function of different RC states and configurations for K = 3 and K = 6, respectively, in out-of-plane orientation. The color map in Fig. 4 was created for the same parameters with in-plane orientation.
For a small number of nuclear spins K = 3 and β = 0, we observe strong fluctuations for both different RC states and different configurations around an average value of S z T ≈ −0.22 as depicted in the color map of Fig. 3. The horizontal stripes dominate over the vertical structures indicating, that the choice of the RC initial states has a greater influence on the results than the spatial configuration of the nuclear spins within the dot.
Moreover, we find large oscillations around this longtime average value for many configurations and initial states. This results in both sizable sample ranges ∆S z and standard deviations σ Sz . By averaging over all 51 × 51 results, we find S z T = (−0.22 ± 0.06) , σ Sz = (0.13 ± 0.04) , and ∆S z = (0.52 ± 0.11) , which is also shown in Fig. 5. The large average value of the sample range ∆S z shows that for most cases analyzed, there was at least one big change in amplitude. However, no total spin flip to + /2 is achieved. The occurrence of sizable standard deviations indicates that there are on average many of these events. Thus in the few nuclear spin regime, coherent oscillations of the electron spin are the dominant dynamics, where recurrences of the initial value take place with a period T P = / max i (|λ i |) ∼ /(A HI · |φ(r K )| 2 ) ≥ 100 ns.
If we consider a larger environment of nuclear spins as presented in Fig. 3 with K = 6, the behavior of the long-time average changes. First of all, the result is much more uniform with respect to both the RC initial states and the configurations. In addition, the remaining differences in S z T depend on the configurations rather than on the initial states, which is obvious from the vertical lines present in this color map. Averaging over all 51×51 results gives S z T = (−0.22 ± 0.02) , which is essentially the same as for K = 3. However, the standard deviation σ Sz = (0.06 ± 0.03) and the sample range ∆S z = (0.37±0.06) clearly decrease. We confirmed this trend of decreasing fluctuations by repeating the above averaging procedure for other numbers of nuclear spins. These results are presented as a function of K in Fig. 5. While the long-time average value is constant, both the standard deviation and the sample range become smaller. Especially, the pronounced decay of the sample range clearly indicates that recurrences occur much less and, hence, that the corresponding recurrence times are increasing with more nuclear spins.
Thus, the major effect of an increased number of nuclear spins is to suppress the oscillations around the longtime average and changing the system from an initial state dominated regime to a regime where the configuration of the nuclear spins is important. This behavior can be understood by analyzing the impact of the nuclear spin number on the dimension of the Hilbert space and on the strength of the hyperfine interaction.
For a small number of nuclear spins, the dimension of the corresponding Hilbert space D = 2 K+1 is small and, hence, we draw our RC initial states from a rather limited set, where individual single product states |n lead to very different dynamics of the electron spin. Due to the combination of only 2 K states |n to a RC initial state, it is not unlikely that one of these states dominates over the rest leading to rather diverse results.
By increasing K and, thus, the Hilbert space dimension this situation is changed. Since the individual state |⇓ . . . ↑↓ of nuclear spins at the border of the dot is al-most irrelevant due to a small |φ( r k )| 2 , groups of effectively equivalent states are superposed. Thus, a more effective averaging is achieved suppressing the dependence on a specific initial state. As a consequence, it is very unlikely for a single state to dominate over the rest.
The coupling strength of nuclei is the key in understanding the dependence of the results on the configuration. Its energy scale is given by the product A HI |φ(r K )| 2 of the hyperfine coupling constant and the maximal value of the envelope function at the sites of the nuclear spins. For a small number K, the probability to find two or more nuclear spins, which couple almost equally with the electron spin, is low due to the large gradient of the envelope function. Hence, effectively only one nuclear spin strongly interacts with the electron leading to simple oscillations. This behavior can be also easily derived by diagonalizing the resulting, effective 4 × 4 matrix of the AHI Hamiltonian given in Eq. (5). In doing so, one finds a discrete spectrum of frequencies given by the degenerate eigenvalues This fact is responsible for the rather uniform dynamics with respect to different configurations in a small K regime.
This situation can of course also occur for larger nuclear spin environments, as shown in Fig. 6 (a) for K = 6. It is, however, rather the exception from the more probable case of several nuclei coupling comparably to the electron, where a almost continuous spectrum is found as depicted in Fig. 6 (b). If we characterize these spectra quantitatively by counting the number of distinct eigenvalues, i.e., eigenvalues which differ significantly, we can map the configuration of the nuclei to the spectra as depicted in Fig. 6 (c).
For the in-plane case, our findings are quite different from the former ones. The electron spin saturates around S z T = 0 for both K = 3 and K = 6 as shown in Fig.  4. Interestingly, we find already for K = 3, that this average is reached very precisely with smaller fluctuations than in the out-of-plane case. This fact becomes also clear from averaging the longitudinal electron spin S z T = (0.000 ± 0.004) over all results. Moreover, the results are independent from the choice of the RC initial state. Some single configurations, however, give rise to deviations from this, where also a dependence on the initial state is restored. It seems, that this is the case, where several nuclear spins couple comparably to the electron spin explaining the sensitivity on initial states. The size of the fluctuations is on average given by σ Sz = (0.15 ± 0.02) . The mean value of the sample range of ∆S z = (0.92 ± 0.07) close to 1 indicates, that in most cases, the electron spin is at least once almost completely flipped to + total spin flips is also strongly decreased for K = 6 as is clear from the sample range ∆S z = (0.51 ± 0.09) . Analyzing this observable as a function of the number of nuclear spins, we observe again a prominent suppression of the fluctuations for growing K as is apparent in Fig.  5.
In order to understand the differences between the inplane and out-of-plane dynamics of the electron spin in more detail, an analytic analysis of the dynamics in the case of only one nuclear spin is very useful. Calculating the long-time average analytically for K = 1 yields Sz T (β) on the orientation of the quantization axis with respect to the graphene plane for K = 6 nuclear spins. The example shown here was calculated for the configuration C = 3, whose continuous spectrum is presented in Fig. 6 (b). The only parameter used to fit the numerical values to the analytic curve Sz T (β) = Sz T (0) · cos 2 (β) is the out-of-plane value Sz T (0).
where the initial density matrix is only non-zero for the electron spin pointing down as for the RC initial states. For more nuclear spins involved, the resulting equations become much more complicated. However, for the special case of only one strongly coupling nuclear spin, the structure of the AHI Hamiltonian remains the same and Eq. (18) still holds.
We investigated the β dependence of the long-time average numerically for some configurations and initial states and K = 2, 4, 6 and 9 nuclear spins, where we find good agreement of our results with S z T (β) = S z T (0) cos 2 (β) with increasing K. Particularly, we observed this behavior also for configurations with several nuclear spins coupling almost equally to the electron spin. As an example, we plot in Fig. 7 the β dependence of S z T for K = 6 and configuration C = 3. Its spectrum is shown in Fig. 6 (b). This fact is also supported by our results presented in Figs. 3 and 4, where we find on average S z T ≈ −0.22 for β = 0 and S z T ≈ 0 for β = π/2. The deviation of S z T (β = 0) from − /4 originates from the finite time window [T min , T max ] used in the numerical calculations, which misses recurrences of the full initial value of S z (t = 0) = − /2.
From this numerical findings and Eqs. (18) and (19), we suppose that contributions from the off diagonal parts cancel each other almost completely and that the elements of the diagonal parts of the density matrix ρ ↓↓ , ρ ↑↑ have approximately equal weight of 1/2 K , which seems reasonable for random complex initial states.

B. Decoherence times
In this section, we want to investigate the decoherence times of the longitudinal electron spin S z for different initial states and different configurations. We chose the threshold to be always about 0.1 below the obtained long-time average, which gives C S = −0.325 for the out-of-plane case β = 0 and C S = −0.1 for the inplane case β = π/2. Moreover, we used exactly the same initial states and configurations for all K as for the calculation of the long-time average. The decoherence times were estimated for times up to 10 7 τ HI ≈ 10 ms with a time resolution ∆T = 10 2 τ HI , which yields at least P = 2π/(∆T ·λ max ) ≈ 20 points per period of the highest absolute frequency max i (|λ i |). For K = 6 we extended the investigated time regime to 10 8 τ HI ≈ 100 ms using the same time resolution ∆T .
As it turns out, the decoherence times obtained by this method are rather independent from the initial states. Several factors are important for this fact. First of all, for larger numbers of nuclei of course the same arguments concerning the Hilbert space dimensions as for the longtime average hold. However, we also find for small K only little dependence on the initial states. One reason for this is probably, that our method is robust against small changes of the longitudinal electron spin caused by different initial states, since we measure when the signal is above a certain threshold, but not how much. Finally, as we show below, the decoherence seems strongly related to the presence of many incommensurate frequencies. These frequencies are proportional to the eigenvalues of the hyperfine Hamiltonian and, hence, independent from the initial state.
Therefore, we focus in the following on the consequences of different configurations on the decoherence times for different numbers of nuclear spins. In principle, there are two relevant aspects concerning the positions of the nuclei, the absolute value of the envelope function |φ( r K )| 2 at the site of the strongest coupling nuclear spin and the relative position of the nuclei with respect to each other. The importance of the former is obvious, since the envelope function sets the maximal energy scale of the AHI in Eq. (9) to A HI ·|φ(r K )| 2 and, consequently, rescales all times by a factor |φ(r K )| −2 . Therefore, if we want to analyze the influence of the relative positions, we have normalize the decoherence times according to We begin our discussion with investigating these normalized decoherence times for K = 6 nuclei in more detail and then turn to absolute decoherence times as a function of K afterwards.
A color map of the normalized decoherence times for 51 initial states and 51 configurations is shown in Fig. 8. For the out-of-plane case, we find that the decoherence times are almost independent of the initial state, but vary over several orders of magnitude for different configurations. If we plot the normalized times as a function of the number of distinct eigenvalues, cf. Fig. 6, we find a direct connection between these times and the configuration of the nuclei in the dot. As is clear from Fig.  9, long decoherence times can be only found for the discrete spectra, which are realized if only one nuclear spin strongly interacts with the electron. The configurations without any decoherence, which are indicated by white spaces in Fig. 8, exhibit discrete spectra with the minimal number of distinct eigenvalues of 3. An example of such a spectrum is shown in Fig. 6 (a). In these cases the dynamics of the longitudinal electron spin are coherent oscillations, where recurrences appear with a period of T P = |φ(r K )| −2 · τ HI ≥ 100 ns. In contrast to this, short normalized decoherence times are a consequence of continuous spectra as presented in Fig. 6 (b). Thus, by the configurations studied, we can proof a direct relation between the relative positions of the nuclear spins and their relative coupling strengths, respectively, and the order of magnitude of the decoherence times.
For the in-plane case, the qualitative picture is similar, however, with shorter normalized decoherence times over all, such that we find decoherence within the investigated times for all configurations. In contrast to the out-ofplane case, also discrete spectra can show rather short decoherence times for specific configurations. Altogether, this demonstrates a much faster decoherence due to the broken symmetry in the in-plane orientation.
Turning from normalized times to absolute decoherence times, the value of the envelope function |φ(r K )| 2 T D > 0.1s In the out-of-plane case, long decoherence times are clearly dominating for few nuclear spins. Increasing K leads to a quick decay of the decoherence times such that short times TD are common. Very short decoherence times start to become relevant for K > 6. (b): In the in-plane case, even for few nuclear spins short relaxation times are the rule. For increasing K, the percentage of short decoherence times is growing further.
at the site of the strongest coupling nuclear spins additionally becomes relevant, since it sets the order of magnitude of all times. Putting a larger and larger number of nuclear spins on a QD of constant area increases the average value of |φ(r K )| 2 , since it is more likely to find a spin very close to the center. Moreover, as we discussed above, an increased K makes it much more probable to have several nuclear spins coupling almost equally to the electron spin. Altogether, this lets us expect a prominent decay of long decoherence times as a function of growing K, which is confirmed by Fig. 10. For β = 0 and very few nuclear spins K = 3, we find that the majority of decoherence times is longer than 10 ms, whereas very short T D are almost completely irrelevant. For K = 8 the percentages of short and long times are inverse. Now, only less than 7% of the decoherence times are longer than 10 ms, while most of the decay of the electron spins takes place within 500 µs. However, for K = 6, surprisingly, still about one-fifth of the cases shows ultra long decoherence times.
In the in-plane orientation, long decoherence times make up only a small fraction even for few nuclear spins. Short decoherence times in the range of 5 µs to 500 µs are significantly increasing for more spins. Notably, ultra short times below 5 µs do not become much more important.
In summary, typical decoherence times are on the order of ms under ideal conditions of small nuclear spin numbers and out-of-plane orientation. In the case of such long decoherence times, of course, other effects like spin orbit coupling could become relevant. In the presence of acoustic phonons and small external magnetic fields, this spin orbit coupling 86,89 can lead to spin relaxation times of T 1 ∼ 1 ms below the decoherence times found here.
For larger numbers of nuclear spins and, generally, for in-plane orientation, decoherence times are smaller, but still above 5 µs. Typical decoherence times of GaAs QDs under spin echo 4 lie in the T 2,echo ∼ 1µs regime, whereas the current record of T 2,CP MG ≈ 200 µs was measured using the Carr-Purcell-Meiboom-Gill (CPMG) pulse sequence 10 . Pure dephasing times T * 2 are below 50 ns for GaAs. Although all our estimates for the decoherence times are done for a model without any effort to improve the coherence of the electron spin like pulse sequences or strong magnetic fields, in almost all considered cases, we are above the GaAs spin echo time T 2,echo . For smaller nuclear spin numbers, graphene even outperforms the CPMG time, which lets us expect very long decoherence times in graphene QDs when using pulse sequences.

V. SUMMARY AND CONCLUSION
Starting from a generic model of a graphene QD, we studied the dynamics of the electron spin caused by the hyperfine interaction with the nuclear spins present in the dot. The number of nuclei was varied from K = 2 to K = 9, where the upper limit corresponds to the natural abundance of spin carrying 13 C for the dot size considered in this article. Besides the role of the number of nuclei, we also investigated the influence of the initial conditions as well as the impact of different configurations of the nuclei in the dot. Moreover, we explored the consequences of the orientation of the spin quantization axis with respect to the graphene plane. In order to characterize and quantify these effects, we analyzed both the long-time average S z T of the longitudinal electron spin component and its decoherence time T D .
Since nuclear spins are usually very hard to control in the envisioned experiments, we chose the initial states to be random complex (RC) superpositions of single product states. For this class of initial states, we found an appreciable effect on the long-time average only in the case of very few nuclear spins with K < 5. Upon increasing the number of nuclear spins the effects of quantum parallelism and amplitude averaging 14,70 reduce the differences between individual RC initial states more and more effectively. In this parameter regime, the results are dominated by the configuration of the nuclear spins within the dot, i.e., by their relative positions with respect to each other. For different configurations, the spectrum of eigenvalues of the hyperfine interaction varies from a highly discrete one with many degenerate eigenvalues to a continuous spectrum with many incommensurate frequencies.
For all K, a pronounced dependence of the long-time average on the orientation angle β between the spin quantization axis and the normal vector of the graphene plane was found. It saturates at approximately one-half of its initial value of S z T ≈ − /4 for β = 0 and at S z T ≈ 0 in the in-plane case with β = π/2. While the long-time average of the electron spin is surprisingly almost constant with respect to K, we observed a strong reduction of fluctuations around it for larger nuclear spin baths.
In contrast to the long-time average, the decoherence times T D never showed a recognizable dependence on the initial states. Instead, the decoherence times depended decisively on the configuration of the nuclear spins in the dot. Long decoherence times were observed for only one nuclear spin strongly interacting with the electron spin, while several almost equally coupled nuclei lead to a very fast decoherence. Moreover, the decoherence times showed a strong dependence on the number of nuclear spins as well as on the orientation of the quantization axis. In the out-of-plane case, about 75% of our results experienced decoherence times longer than T max = 10ms for K = 3. For K = 8, instead, less than 10% showed no decoherence within this time frame while, in most cases, the electron spin decayed in less than 500 µs. Considering the in-plane orientation, already for K = 3 the majority of investigated initial state / configuration sets decohere within 500 µs.
Although our results were obtained for a specific model of the graphene quantum dot using a Gaussian envelope function, they could be generalized quite naturally. In our model, the QD was comparably small with a sharp boundary. This choice resulted in a steep envelope function. Physically, this situation corresponds approximately to an etched QD. Thinking of larger QDs with smoother boundaries, we expect a flatter envelope function which gives rise to more nuclear spins interacting comparably with the electron spin. Consequently, it becomes more likely to end up with rather low fluctuations around the long-time average and to find quite short decoherence times. In contrast, the realization of even smaller dots 60 with diameters of about 1 nm causes a very steep envelope function. This case should result in, at most, one nuclear spin interacting with the electron spin.
Both scenarios seem experimentally interesting in order to engineer QDs for different applications. A 13 C enriched QD could potentially be used to prepare the electron spin very precisely in a certain superposition of spin up and down for subsequent experiments. A very small QD, in contrast, could serve as a storage for the electron spin where very long decoherence times are to be expected.
Besides these technical points, graphene QDs could also serve as a rich playground to test fundamental aspects of quantum mechanics and quantum information theory in an interesting system-bath setup. If we consider the hosted electron spin as the system, we are able to control both its spatial size and its state electrostatically. Moreover, the electron spin can be straightforwardly addressed via external magnetic fields or in an optical way. In contrast, a direct preparation of the state of the nuclear spin bath seems challenging. However, the size of the nuclear spin bath can be modified systematically by isotopic purification of either 12 C or 13 C. Finally, as we argued above, the design of the QD enables the experimentalist to manipulate the strength as well as the nature of the system-bath interaction. Thus, these considerations render graphene QDs to be a flexible system-bath realization offering a controllable, fermionic bath of spin 1/2 nuclei.
Given these opportunities, our setup not only seems very promising for studying a quantum to classical crossover as a function of the bath size in a fermionic environment, but also for more advanced concepts of quantum information theory such as quantum Darwinism 91 . With the notion of "quantum Darwinism," W. H. Zurek summarizes his ideas of emergent classicality in a pure quantum universe, where the formation of classical and quantum system bath correlations, as well as the accompanied information exchange matter. Since measurements are most often indirect, relying on the environment as a mediator, it seems to be evident that the environment plays the crucial role in the creation of objective properties. However, as far as we know, this theory is up to now tested only in few experiments 92,93 in a rather indirect way. In our opinion, taking advantage of the controllable spin bath in graphene QDs, this system offers a unique opportunity to reveal new insights in the role of the environment in the classical world we experience.