Simple and accurate method for central spin problems
SSimple and accurate method for central spin problems
Lachlan P. Lindoy and David E. Manolopoulos
Department of Chemistry, University of Oxford, Physical and TheoreticalChemistry Laboratory, South Parks Road, Oxford, OX1 3QZ, UK
We describe a simple quantum mechanical method that can be used to obtain accurate numericalresults over long time scales for the spin correlation tensor of an electron spin that is hyperfinecoupled to a large number of nuclear spins. This method does not suffer from the statistical errorsthat accompany a Monte Carlo sampling of the exact eigenstates of the central spin Hamiltonianobtained from the algebraic Bethe ansatz, or from the growth of the truncation error with time inthe time-dependent density matrix renormalization group (t-DMRG) approach. As a result, it canbe applied to larger central spin problems than the algebraic Bethe ansatz, and for longer timesthan the t-DMRG algorithm. It is therefore an ideal method to use to solve central spin problems,and we expect that it will also prove useful for a variety of related problems that arise in a numberof different research fields.
The central spin Hamiltonianˆ H = B ˆ S z + N X j =1 a j ˆ S · ˆ I j (1)has been widely studied in the condensed matter physicsliterature because of its relevance to the hyperfine-induced decoherence of electron spins on quantumdots. It also arises as an important ingredient inthe Hamiltonians that govern the dynamics of the elec-tron spins in radical pairs of chemical and biochem-ical interest, and the spin dynamics of polaronpairs of interest to the organic semiconducting devicecommunity.
There is therefore considerable currentinterest in being able to calculate the exact quantum me-chanical central (ˆ S ) spin dynamics of the Hamiltonian inEq. (1).In the context of an unpaired electron on a quantumdot, the first term in Eq. (1) is the Zeeman interactionbetween the electron spin and an external magnetic field,and the second term is a sum of hyperfine interactions be-tween the electron spin and the spins of the nuclei in thedot. The Zeeman interactions of the nuclear spins havebeen eliminated by defining B = ( g e − g n ) h , where g e and g n are the electron and nuclear spin g-factors and h is the magnetic field strength. Since ˆ J z = ˆ S z + P Nj =1 ˆ I jz is a constant of the motion, the difference g n h ˆ J z betweenEq. (1) and the Hamiltonian that includes nuclear Zee-man interactions is a trivial energy shift within each ˆ J z symmetry block. The dipolar interactions between thenuclear spins in the dot have been neglected, but sincethey are much weaker than the terms that have been re-tained in Eq. (1) they will only affect the electron spindynamics at very long times.The difficulty of the problem lies in the fact that thesize of the Hilbert space increases exponentially with N .Assuming for simplicity that all of the nuclear spins have I = 1 /
2, the dimension of the Hilbert space is 2 N +1 .This makes solving the problem by numerical diagonal-ization impractical for N &
20. However, quantum dotstypically have many more than 20 nuclear spins with non- negligible hyperfine interactions, and the same is true ofmany interesting radicals and polarons.In the condensed matter physics literature, severalhighly sophisticated techniques have been developed toovercome the exponential scaling. Perhaps the most el-egant of these exploits the fact that the Hamiltonian inEq. (1) has the form of a Gaudin magnet in an externalmagnetic field, which is well known to be an integrableproblem.
In other words, there are N constants of themotion in addition to ˆ H , of which ˆ J z is just one. Theeigenvalues and eigenstates of all N + 1 conserved opera-tors, and therefore of ˆ H itself, can be found by solving asystem of N + 1 coupled algebraic Bethe equations. The solution of these equations is obvious when 1 /B → /B . But once this has been done, one is still facedwith the problem of summing over all of the eigenstatesof ˆ H to calculate the dynamics of the central electronspin, with a numerical effort that again scales exponen-tially with N . Faribault and Schuricht have suggested aMonte Carlo sampling of the eigenstates to alleviate thisproblem, and demonstrated that this implementation ofthe algebraic Bethe ansatz is capable of solving centralspin problems with up to N = 48 nuclear spins. How-ever, the statistical errors from the Monte Carlo samplingare quite noticeable in their results, and one would expectthese errors to become even larger for larger N .A contrasting approach has been to adapt the time-dependent density matrix renormalization group (t-DMRG) method to the topology of the central spinproblem. Since the t-DMRG algorithm is bestsuited to a one-dimensional chain of spins with nearest-neighbour interactions, this requires some ingenuity. Itcan however be done, as Uhrig and co-workers haveshown in a pair of papers that include well convergedshort-time results for a number of central spin problemswith up to 999 hyperfine-coupled nuclear spins.
Thetrouble with this method is that it is limited to shorttime scales. From Eq. (1), the coupling between any twonuclear spins via the electron spin involves only two suc-cessive Hamiltonian interactions. This implies that theentanglement of the exact time-evolved wavefunction will a r X i v : . [ qu a n t - ph ] A ug grow rapidly with time, and that it will become increas-ingly difficult to represent it with the tensor product formthat is assumed in t-DMRG. Hence the truncation errorin this algorithm is bound to grow with time. Uhrig andco-workers have shown that this error becomes unaccept-able beyond 40 τ for a problem with N = 99 nuclear spins,where τ = N X j =1 a j − (2)is the characteristic time scale of the electron spin pre-cession in the nuclear hyperfine field. The t-DMRG trun-cation error is also expected to grow with N .Motivated by these issues with existing algorithms, weshall now present an entirely different approach to theproblem. Rather than attempting to calculate the centralspin dynamics of the Hamiltonian in Eq. (1), we shallinstead construct a sequence of simpler Hamiltonians ˆ H M for M = 1 , , . . . in such a way that their central spindynamics converges to that of the original Hamiltonianwith increasing M .The Hamiltonians ˆ H M we shall consider areˆ H M = B ˆ S z + M X j =1 A j N j X i =1 ˆ S · ˆ I ij , (3)in which the modified hyperfine coupling constants A j and the numbers N j of equivalent nuclear spins in eachof the M blocks are chosen to ensure that the first M + 1moments of the modified hyperfine distribution coincidewith those of the original distribution: µ k = N X j =1 a kj = M X j =1 N j A kj for k = 0 , , . . . , M. (4)In order to obtain these Hamiltonians, we first use adiscrete procedure of Stieltjes to construct a Gaussianquadrature rule with non-integer weights W j and nodes ¯ A j such that µ k = N X j =1 a kj = M X j =1 W j ¯ A kj for k = 0 , , . . . , M − . (5)We then search through all floor and ceiling possibilitiesfor the integers N j = b W j c or d W j e in Eq. (4) that areconsistent with the k = 0 moment constraint M X j =1 N j = N, (6)and use Newton’s method to solve the remaining momentequations in Eq. (4) for the M unknowns { A j } , startingfrom the initial guess { A j } = (cid:8) ¯ A j (cid:9) . Finally, we selectfrom among the < M possible solutions for the set ofintegers { N j } the one that minimises the error in µ M +1 . The reason for choosing simplified Hamiltonians ˆ H M ofthe form in Eq. (3) is that the presence of equivalent nu-clei dramatically simplifies the spin dynamics calculation.For example, a calculation with M = 1 and N = 100equivalent spin-1/2 nuclei can be reduced by symmetryto 51 separate calculations, each of which involves a sin-gle resultant spin ˆ I = P N i =1 ˆ I i with I between 0 and50. The Hilbert space of the largest of these calculationshas dimension 2 ×
101 = 202, whereas the dimension ofthe full Hilbert space is 2 . A calculation with M = 2and N = N = 50 can be reduced in the same way to26 ×
26 = 676 separate calculations with I and I be-tween 0 and 25. The Hilbert space of the largest of thesecalculations has dimension 2 × ×
51 = 5 , N , the effort increases with M , but itremains many orders of magnitude smaller than that ofthe full calculation for all M (cid:28) N .The choice of model Hamiltonians ˆ H M that have hy-perfine distributions with the same first moments asthose of the exact Hamiltonian ˆ H is physically motivated.The early semiclassical theory of Schulten and Wolynes only involves the second moment µ of the hyperfine dis-tribution, and yet it becomes exact for the central spindynamics in the limit as N → ∞ . By constructing ourmodel Hamiltonians on the basis of Eq. (4), we are there-fore guaranteed to obtain the correct result in the large N limit for any M ≥
2. And since Eq. (3) reduces to Eq. (1)when M = N (with A j = a j and N j = 1 for all j ), we arealso guaranteed to obtain the correct result in the small N limit in which a calculation based on Eq. (1) is feasi-ble. The only remaining question, therefore, is how wellthis method actually works for the intermediate values of N that are of interest in physical applications.In order to answer this question, we have used themethod to calculate the infinite temperature central spincorrelation tensors R αβ ( t ) = 1 Z tr h ˆ S α (0) ˆ S β ( t ) i , (7)with Z = 2 N +1 , for two model problems that havebeen studied previously using the t-DMRG and algebraicBethe ansatz methods. The first of these has a uniformdistribution of hyperfine coupling constants a j τ = r N N + 3 N + 1 N + 1 − jN , (8)and the second has an exponential distribution that arisesfrom the Fermi contact interactions in a two-dimensionalquantum dot with a Gaussian electronic wavefunction a j τ = s − e − / ( N − − e − N/ ( N − e − ( j − / ( N − , (9)where N is the number of nuclear spins within one Bohrradius of the Gaussian wavefunction and τ ≡ p /µ isthe timescale in Eq. (2). M = 20.030.060.090.12 M = 30.030.060.090.12 M = 40 20 40 60 80 1000.000.030.060.090.12 M = 5 t / R zz ( t ) FIG. 1. Convergence of R zz ( t ) for Model I with N = 49, as afunction of M . The solid black curve in each panel shows thefully-converged result obtained with M = 7, and the dashedred curve shows the result obtained with the specified valueof M . In what follows, we shall refer to these two modelsas Model I and Model II, respectively. Uhrig and co-workers used the t-DMRG method to study the low-field( B = 0) limit of Model I for N = 49, 99, and 999, andFaribault and Schuricht used the algebraic Bethe ansatzto study the low field limit of Model II with N = 48 and N = 24 and 36. These are the model problems we shallconsider. Since both have B = 0, we have that R xx ( t ) = R yy ( t ) = R zz ( t ) and R xy ( t ) = R yz ( t ) = R zx ( t ) = 0.It therefore suffices to consider just R zz ( t ), from whichall other properties of the central spin dynamics can berecovered. We have calculated R zz ( t ) for both of these models us-ing the method outlined above. Rather than numericallydiagonalising the matrix representations of the Hamil-tonians ˆ H M in Eq. (3), we found it more efficient touse a symplectic time-dependent wavepacket propagationalgorithm. The traces in Eq. (7) were evaluated deter-ministically for each symmetry block of ˆ H M containingless than 1000 states, and stochastically using a recently-developed spin coherent state algorithm for the largersymmetry blocks. Note that this is quite different from N = 490.050.100.150.200.25 R zz ( t ) N = 990 20 40 60 80 100 t /0.000.050.100.150.200.25 N = 999
10 20 30 400.0650.0750.085 10 20 30 40 500.0650.0750.085
FIG. 2. Comparison of the present results for Model I(dashed red curves) with the t-DMRG results of Uhrig andco-workers (solid black curves), for N = 49, 99, and 999. stochastically sampling the eigenstates of the full Hamil-tonian ˆ H as was done by Faribault and Schuricht. Ourstatistical errors were found to be negligible with just1000 Monte Carlo samples of the initial nuclear spin co-herent states in each of the larger symmetry blocks.Fig. 1 shows the convergence of the correlation function R zz ( t ) obtained from the present method with increasing M , for Model I with N = 49 nuclear spins. One seesthat the method remains accurate for longer times as M increases, and that it is fully converged over the timeinterval considered (up to t = 100 τ ) by the time M = 5.Similar convergence plots for N = 99 and 999 are givenin the supplementary material. Full convergence out to100 τ is obtained with M = 4 for N = 99, and withjust M = 3 for N = 999. Thus the method becomesincreasingly efficient with increasing N . Presumably thisis because N = 999 is approaching the large- N limit inwhich the central spin dynamics is determined entirelyby the second moment of the hyperfine distribution.Fig. 2 compares our converged results for Model Iwith the t-DMRG results of Uhrig and co-workers. The good agreement over the timescales for which the t-DMRG results are available confirms that the presentmethod converges on the correct quantum mechanical N = 240 20 40 60 80 100 t /0.000.050.100.150.200.25 N = 36
60 70 80 900.0580.0610.064 60 70 80 900.0680.0710.074 R zz ( t ) FIG. 3. Comparison of the present results for Model II(dashed red curves) with the algebraic Bethe ansatz resultsof Faribault and Schuricht (solid black curves), for N = 48. central spin dynamics of this model. It also avoids thegrowth of the truncation error with time in the t-DMRGalgorithm, which is especially apparent in the inset ofthe N = 99 panel in the figure. The t-DMRG results for N = 999 are only available out to t = 10 τ because thetruncation error increases with increasing N .The convergence tests we have performed for Model II(with N = 48 nuclear spins) are provided in the supple-mentary material. M = 5 was again found to be suffi-cient to converge R zz ( t ) out to t = 100 τ , in much thesame way as for Model I with N = 49 nuclear spins.Fig. 3 compares the resulting correlation functions withthose obtained by Faribault and Schuricht using the al-gebraic Bethe ansatz. The agreement between the twosets of calculations is again excellent, and confirms thatboth methods are giving the correct correlation functionsfor this model. The only difference between the two setsof results is that our correlation functions are smootherthan those of Faribault and Schuricht, which have no-ticeable stochastic errors associated with the incompleteMonte Carlo sampling of the eigenstates of ˆ H in theirmethod. As we have already mentioned, these stochasticerrors in the algebraic Bethe ansatz method are expectedto become even more pronounced for larger N .Finally, let us return to Model I and the question ofhow close N = 999 is to the large N limit in which theearly semiclassical theory of Schulten and Wolynes isexpected to become exact. Fig. 4 compares the presentquantum mechanical results for this model with thosegiven by their theory, and also with those of an im-proved semiclassical theory suggested by Manolopoulosand Hore. N = 490.050.100.150.200.25 R zz ( t ) N = 990 20 40 60 80 100 t /0.000.050.100.150.200.25 N = 999 SWQMSC10 20 30 400.0650.0750.085 10 20 30 40 500.0650.0750.085
FIG. 4. Comparison of the present quantum mechanical (QM)results for Model I with the results of the early semiclassicaltheory of Schulten and Wolynes (SW) and the more recentsemiclassical theory of Manolopoulos and Hore (SC). Thelong-time behaviours of these correlation functions are dis-cussed in the supplementary material. The quantum mechanical results in Fig. 4 are suffi-ciently well converged, and for sufficiently long times, toprovide a stringent test of these semiclassical approxi-mations. One sees from the figure that, although bothsemiclassical theories are qualitatively reasonable, nei-ther is quantitatively accurate, even for N = 999. TheSchulten-Wolynes theory misses the long-time decay ofthe central spin correlation function, and the improvedsemiclassical theory predicts too rapid a long-time decay.This is even more apparent in Fig. 5, which shows R xx ( t )over a longer time scale for Model I with N = 49 in amagnetic field of strength B = 1 / τ .In view of this, it would be interesting to apply thepresent method to a variety of problems that have previ-ously only been studied semiclassically. One such prob-lem is the spin dynamics of a photoexcited carotenoid-porphyrin-fullerene radical pair that has been shown tobe sensitive to an Earth strength magnetic field. Sincethis field is so weak ( ∼ µ T), and the experiment thatwas used to detect it measured a tiny field-on minus field-off difference signal, it is conceivable that the semi-classical calculations that have been used to study the t /0.050.000.050.100.150.200.25 R xx ( t ) SWQMSC
FIG. 5. Comparison of SW, QM, and SC results for R xx ( t )for Model I with N = 49 in a magnetic field of strength B =1 / τ . Note that the time scale available from a t-DMRGcalculation ( t < τ ) would not be enough to predict the long-time behaviour of this correlation function, which is discussedfurther in the supplementary material. problem might not have been sufficiently accurate. The carotenoid radical in the pair contains 45 protons withsignificant hyperfine interactions, which have so far pre-cluded an exact quantum mechanical calculation. Theresults we have presented here show that the presentmethod could easily handle this problem, and many otherinteresting spin dynamics problems as well. ACKNOWLEDGMENTS
We are grateful to Pranav Singh for a helpful commenton the first draft of this manuscript. Lachlan Lindoy issupported by a Perkin Research Studentship from Mag-dalen College, Oxford, an Eleanor Sophia Wood Post-graduate Research Travelling Scholarship from the Uni-versity of Sydney, and by a James Fairfax Oxford Aus-tralia Scholarship. The authors would like to acknowl-edge the use of the University of Oxford Advanced Re-search Computing (ARC) facility in carrying out thiswork. See http://dx.doi.org/10.5281/zenodo.22558. I. A. Merkulov, A. L. Efros and M. Rosen, Phys. Rev. B , 205309 (2002). A. Khaetskii, D. Loss and L. Glazman, Phys. Rev. B ,195329 (2003). S. I. Erlingsson and Y. V. Nazarov, Phys. Rev. B ,205327 (2004). P. F. Braun, X. Marie, L. Lombez, B. Urbaszek, T. Amand,P. Renucci, V. K. Kalevich, K. V. Kavokin, O. Krebs,P. Voisin and Y. Masumoto, Phys. Rev. Lett. , 116601(2005). K. A. Al-Hassanieh, V. V. Dobrovitski, E. Dagotto and B.N. Harmon, Phys. Rev. Lett. , 037204 (2006). G. Chen, D. L. Bergman and L. Balents, Phys. Rev. B ,045312 (2007). E. Barnes, L. Cywinski and S. Das Sarma, Phys. Rev. Lett. , 140403 (2012). U. E. Steiner and T. Ulrich, Chem. Rev. , 51 (1989). C. T. Rodgers, Pure Appl. Chem. , 19 (2009). P. J. Hore, Proc. Natl. Acad. Sci. USA , 1357 (2012). B. Koopmans, W. Wagemans, F. Bloom, P. A. Bobbert,M. Kemerink and M. Wohlgenannt, Philos. Trans. R. Soc.London Ser. A , 3602 (2011). M. Wohlgenannt, Phys. Status Solidi (RRL) , 229 (2012). E. Ehrenfreund and Z. V. Vardeny, Isr. J. Chem. , 552(2012). A. Faribault and D. Schuricht, Phys. Rev. B , 085323(2013). M. Gaudin, J. Phys. (Paris) , 1087 (1976). M. Gaudin,
La fonction d’onde de Bethe (Masson, Paris,1983);
The Bethe wavefunction (Cambridge UniversityPress, 2014). A. Faribault and D. Schuricht, Phys. Rev. Lett. , 040405 (2013). D. Stanek, C. Raas and G. S. Uhrig, Phys. Rev. B ,155305 (2013). D. Stanek, C. Raas and G. S. Uhrig, Phys. Rev. B ,064301 (2014). W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P.Flannery,
Numerical recipes in fortran (2nd Edition, Cam-bridge University Press, 1992), pp. 151-153. K. Schulten and P. G. Wolynes, J. Chem. Phys. , 3292(1978). For example, the observable considered by Faribault andSchuricht in Ref. 14 is Re [ h S + ( t ) i ] / h S + (0) i with ρ (0) =(2 S x + 1) /Z . It is easy to show that this is the same as4 R xx ( t ), which coincides with 4 R zz ( t ) when B = 0. S. K. Gray and D. E. Manolopoulos, J. Chem. Phys. ,7099 (1996). A. M. Lewis, T. P. Fay and D. E. Manolopoulos, J. Chem.Phys. , 244101 (2016). D. E. Manolopoulos and P. J. Hore, J. Chem. Phys. ,124106 (2013). The Schulten-Wolynes theory considers the electron spinprecession in a static nuclear hyperfine field, as discussedin more detail in the context of the central spin problem byMerkulov et al. (Ref. 1). The improved semiclassical the-ory combines the electron and nuclear spin precessions ina classical vector model with p S ( S + 1) and p I j ( I j + 1)vector lengths, in a way that is more consistent with New-ton’s third law. K. Maeda, K. B. Henbest, F. Cintolesi, I. Kuprov, C. T.Rodgers, P. A. Liddell, D. Gust, C. R. Timmel and P. J.Hore, Nature , 387 (2008). A. M. Lewis, D. E. Manolopoulos and P. J. Hore, J. Chem.Phys. , 044111 (2014). simple and accurate method for central spin problems:Supplementary material
Lachlan P. Lindoy and David E. Manolopoulos
Department of Chemistry, University of Oxford, Physical and TheoreticalChemistry Laboratory, South Parks Road, Oxford, OX1 3QZ, UK
This supplementary material contains four sections. The first provides the computer programthat we used to generate the simplified Hamiltonians ˆ H M in Eq. (3) of the text, along with exampleinput and output files. This information is provided so as to make our results reproducible, andto enable others to apply our method to other spin dynamics problems. The second contains moredetails on how to exploit the presence of equivalent spins in spin dynamics calculations, and thethird contains plots showing the convergence of our method with increasing M for Model I with N = 99 and 999, and for Model II with N = 24 and 36. The final section discusses semi-log plots( R αα ( t ) vs log t ) of the data in Figs. 4 and 5 of the manuscript, which provide a clearer picture ofthe long-time behaviour of the various (Schulten-Wolynes, semiclassical, and quantum mechanical)central spin correlation functions than the linear plots given in the manuscript. a r X i v : . [ qu a n t - ph ] A ug I. HOW TO CONSTRUCT THE HAMILTONIANS ˆ H M The text of the paper has simply outlined how we have constructed the simplified Hamiltonians ˆ H M in Eq. (3).Here we provide the computer code that we have actually used to do this, so that others can both reproduce ourresults and apply the same method to other spin dynamics problems.The following fortran program reads in 3 input parameters, N , N , and M . It then constructs ˆ H M , and outputsthe optimised hyperfine coupling constants { A j } M and numbers { N j } M of equivalent nuclear spins in each of itssymmetry blocks. It also outputs the percentage errors in the first M + 1 moments of the hyperfine distribution thatresult from replacing the original Hamiltonian with the simplified Hamiltonian.The results for the first central spin problem (Model I) considered in the text were obtained by inputting N = 0, N = 49, 99, and 999, and M = 5, 4, and 3, respectively. The results for the second problem (Model II) were obtainedby inputting N = 24 and 36, N = 48, and M = 5. The convergence tests we shall present below were obtained usingother values of M between 2 and 7. program momentsimplicit double precision (a-h,o-z)cc ------------------------------------------------------------------c This program constructs a simplified central spin Hamiltonainc with M sets of equivalent nuclei from a given Hamiltonian withc N inequivalent nuclei.c ------------------------------------------------------------------c allocatable :: a(:),abar(:),nbar(:)cc Setup:c read (5,*) n0,n,mallocate (a(n),abar(m),nbar(m))if (n0 .eq. 0) then! uniform hyperfine distributionfac = sqrt((6.0d0*n)/(2.0d0*n*n+3.0d0*n+1.0d0))do j = 1,na(j) = fac*(n-(j-1.0d0))/nenddoelse if (n0 .gt. 1) then! exponential hyperfine distributiontop = 1.d0-exp(-2.d0/(n0-1.d0))btm = 1.d0-exp(-(2.d0*n)/(n0-1.d0))fac = sqrt(top/btm)do j = 1,na(j) = fac*exp(-(j-1.d0)/(n0-1.d0))enddoelsestop ’moments 1’endifcc Calculation:c call shrink (a,n,abar,nbar,m)cc Output:c write (6,600) n0,n,m600 format(1x,’ N0 = ’,i5,’ N = ’,i5,’ M = ’,i5/1x,/1x,+ ’ j N_j A_j’/1x)ntot = 0do j = 1,mwrite (6,601) j,nbar(j),abar(j)
601 format(1x,i6,i9,f20.12)ntot = ntot+nbar(j)enddoif (ntot .ne. n) stop ’moments 2’write (6,602)602 format(/1x,’ k % error in mu_k’/1x)do k = 1,m+1exact = 0.d0do i = 1,nexact = exact+a(i)**kenddoapprox = 0.d0do j = 1,mapprox = approx+nbar(j)*abar(j)**kenddoerror = 100.d0*abs(exact-approx)/abs(exact)write (6,603) k,error603 format(1x,i6,f29.6)if (k .eq. m) print*enddodeallocate (a,abar,nbar)stopendsubroutine shrink (a,n,abar,nbar,m)implicit double precision (a-h,o-z)cc ------------------------------------------------------------------c Optimally approximates N inequivalent nuclei with hyperfinec coupling constants a(j) by M < N sets of equivalent nucleic with hyperfine constants abar(j), with nbar(j) nuclei in set j.c ------------------------------------------------------------------cc We use integer*8 so that we can deal with all M < 64:c integer*8 ib,nsetdimension a(n),as(n),awbar(m),abar(m),nbar(m),nfloor(m)dimension atemp(m),ntemp(m)dimension w(n),wbar(m)cc It only makes any sense to call this subroutine with M < N:c if (m .gt. n) stop ’M > N in shrink?’if (m .eq. n) stop ’M = N in shrink?’cc Optimum solution with non-integer weights:c do i = 1,nw(i) = 1.d0enddocall qrule (n,w,a,m,wbar,awbar)nftot = 0do j = 1,mnfloor(j) = int(wbar(j))nftot = nftot+nfloor(j)enddondiff = n-nftotc c Optimum solution with integer weights:c dmom = 1.d10do ib=0,LSHIFT(1,m)-1call count_set_bits(ib,nset)if(nset .eq. ndiff) thendo i=1,Mnind = RSHIFT(IAND(ib,LSHIFT(1,i-1)),i-1)ntemp(i) = nfloor(i)+nindatemp(i) = awbar(i)enddocall newton (awbar,wbar,atemp,ntemp,m,icode)if (icode.eq.0) thenemom = 0.0d0do k = 1,m+1exact = 0.0d0do i = 1,nexact = exact+a(i)**kenddoapprox = 0.0d0do j = 1,mapprox = approx+ntemp(j)*atemp(j)**kenddoemom = emom+abs(approx-exact)enddoif (emom.lt.dmom) thendmom = emomdo i = 1,mnbar(i) = ntemp(i)abar(i) = atemp(i)enddoendifendifendifenddoreturnendsubroutine count_set_bits (i,nset)implicit nonecc ------------------------------------------------------------------c Determines the total number of non-zero bits in a 64 bit integerc using a 64 bit implementation of the Hamming Weight algorithm.c ------------------------------------------------------------------c integer*8 i,nset,vinteger*8 m1,m2,m3,m4DATA m1 /Z’5555555555555555’/, m2/Z’3333333333333333’/DATA m3 /Z’f0f0f0f0f0f0f0f’/, m4/Z’101010101010101’/c v = iv = v-IAND(RSHIFT(v,1),m1)v = IAND(V,m2)+IAND(RSHIFT(V,2),m2)v = IAND(v+RSHIFT(v,4),m3)nset = RSHIFT(v*m4,(SIZEOF(nset)-1)*8)returnend subroutine newton (awbar,wbar,abar,nbar,m,icode)implicit double precision (a-h,o-z)cc ------------------------------------------------------------------c Uses Newton’s method to solve M non-linear moment equations in Mc unknowns (the hyperfine constants of M sets of equivalent nuclei),c starting from an appropriate initial guess.c ------------------------------------------------------------------c parameter (maxit = 30)dimension awbar(m),wbar(m),abar(m),nbar(m)dimension f(m),g(m,m),indx(m)c do iter = 1,maxitdo k = 1,mf(k) = 0.d0do j = 1,mf(k) = f(k)+nbar(j)*abar(j)**k-wbar(j)*awbar(j)**kg(k,j) = k*nbar(j)*abar(j)**(k-1)enddoenddocall rgfac (g,m,m,indx,ierr)if (ierr .ne. 0) thenicode = -1returnendifcall rgsol (g,m,m,indx,f)absa = 0.d0absd = 0.d0do j = 1,mabsa = absa+abs(abar(j))absd = absd+abs(f(j))g(j,1) = abar(j)abar(j) = abar(j)-f(j)enddoda = absd/absaif (da .gt. 1.d0) thenicode = 1returnendifif (1.d0+da .eq. 1.d0) thenicode = 0returnendifif (iter.gt.1 .and. da.ge.dap) thendo j = 1,mabar(j) = g(j,1)enddoicode = 0returnendifdap = daenddoicode = 2returnend subroutine qrule (np,wp,xp,n,w,x)implicit double precision (a-h,o-z)cc -----------------------------------------------------------------c Uses a discrete Stieltjes procedure to construct an n-pointc contracted quadrature rule from an np-point primitive quadraturec rule with the same (non-negative) weight function.c -----------------------------------------------------------------c dimension wp(np),xp(np),w(n),x(n)dimension p(np,n),q(np),v(n)c do i = 1,npq(i) = sqrt(wp(i))enddodo k = 1,nqq = 0.d0do i = 1,npqq = qq+q(i)**2enddoqq = sqrt(qq)if (qq .eq. 0.d0) stop ’qrule 1’do i = 1,npp(i,k) = q(i)/qqq(i) = xp(i)*p(i,k)enddodo j = k,1,-1pq = 0.d0do i = 1,nppq = pq+p(i,j)*q(i)enddoif (j .eq. k) thenv(k) = 0.d0w(k) = qqx(k) = pqendifdo i = 1,npq(i) = q(i)-p(i,j)*pqenddoenddoenddoweight = w(1)v(1) = 1.d0ldv = 1call rstqlv (x,w,n,v,ldv,ierr)if (ierr .ne. 0) stop ’qrule 2’do j = 1,nw(j) = (weight*v(j))**2enddoreturnendsubroutine rstqlv (d,e,n,v,ldv,ierr)implicit double precision (a-h,o-z)cc -----------------------------------------------------------------c Eigenvalues and eigenvectors of a real symmetric tridiagonalc matrix. Based on the Numerical Recipes routine tqli but c modified so as to calculate just the first ldv rows of thec eigenvector matrix. (Only the first row is required in thec discrete Stieltjes procedure.)cc Note that v must be initialized to (the first ldv rows of)c a unit matrix before entry.c -----------------------------------------------------------------c dimension d(n),e(n),v(ldv,n)c nv = min(ldv,n)do i = 2,ne(i-1) = e(i)enddoe(n) = 0.d0do l = 1,niter = 01 do m = l,n-1dd = abs(d(m))+abs(d(m+1))ee = abs(e(m))if (ee+dd .le. dd) goto 2enddom = n2 if (m .ne. l) thenif (iter .eq. 30) thenierr = 1returnendifiter = iter+1g = (d(l+1)-d(l))/(2.d0*e(l))r = sqrt(1.d0+g**2)if (abs(g-r) .gt. abs(g+r)) r = -rg = d(m)-d(l)+e(l)/(g+r)s = 1.d0c = 1.d0p = 0.d0do i = m-1,l,-1f = s*e(i)b = c*e(i)r = sqrt(f**2+g**2)e(i+1) = rif (abs(r) .eq. 0.d0) thend(i+1) = d(i+1)-pe(m) = 0.d0goto 1endifs = f/rc = g/rg = d(i+1)-pr = (d(i)-g)*s+2.d0*c*bp = s*rd(i+1) = g+pg = c*r-bdo k = 1,nvf = v(k,i+1)v(k,i+1) = s*v(k,i)+c*fv(k,i) = c*v(k,i)-s*fenddo enddod(l) = d(l)-pe(l) = ge(m) = 0.d0goto 1endifenddodo j = 1,n-1k = jdo i = j+1,nif (d(i) .lt. d(k)) k = ienddoif (k .ne. j) thenswap = d(k)d(k) = d(j)d(j) = swapdo i = 1,nvswap = v(i,k)v(i,k) = v(i,j)v(i,j) = swapenddoendifenddoierr = 0returnendsubroutine rgfac (a,lda,n,indx,ierr)implicit double precision (a-h,o-z)cc -----------------------------------------------------------------c LU decomposition routine for real matrices.cc Uses the same pivoting strategy as Numerical Recipes ludcmp,c but with a re-organization of the inner loops to exploitc sparsity and reduce the loop overhead.c -----------------------------------------------------------------c dimension a(lda,n),indx(n)dimension vv(n)c do i = 1,naamax = 0.d0do j = 1,nif (abs(a(i,j)) .gt. aamax) aamax = abs(a(i,j))enddoif (aamax .eq. 0.d0) thenierr = 1returnendifvv(i) = 1.d0/aamaxenddodo j = 1,ndo k = 1,j-1if (a(k,j) .ne. 0.d0) thendo i = k+1,na(i,j) = a(i,j)-a(i,k)*a(k,j)enddo endifenddoaamax = 0.d0do i = j,natest = vv(i)*abs(a(i,j))if (atest .ge. aamax) thenimax = iaamax = atestendifenddoif (j .ne. imax) thendo k = 1,nswap = a(imax,k)a(imax,k) = a(j,k)a(j,k) = swapenddovv(imax) = vv(j)endifindx(j) = imaxif (a(j,j) .eq. 0.d0) thenierr = 2returnendifif (j .ne. n) thenpivot = 1.d0/a(j,j)do i = j+1,na(i,j) = a(i,j)*pivotenddoendifenddoierr = 0returnendsubroutine rgsol (a,lda,n,indx,b)implicit double precision (a-h,o-z)cc -----------------------------------------------------------------c Uses the factorized A from rgfac to solve the linearc equations A*X = B, overwriting the solution X on B.c Like the Numerical Recipes routine lubksb, but withc a slightly different calling sequence.c -----------------------------------------------------------------c dimension a(lda,n),indx(n),b(n)c ii = 0do i = 1,nll = indx(i)sum = b(ll)b(ll) = b(i)if (ii .ne. 0) thendo j = ii,i-1sum = sum-a(i,j)*b(j)enddoelse if (sum .ne. 0.d0) thenii = iendif b(i) = sumenddodo i = n,1,-1do j = i+1,nb(i) = b(i)-a(i,j)*b(j)enddob(i) = b(i)/a(i,i)enddoreturnend This program can be compiled with the gfortran compiler, using the -Os compiler option. It generates the followingoutput file for Model I with N = 999 and M = 3: N0 = 0 N = 999 M = 3j N_j A_j1 278 0.0062115350152 444 0.0274385274273 277 0.048627307400k % error in mu_k1 0.0000002 0.0000003 0.0000004 0.000005
And the following output file for Model II with N = 24 and M = 5: N0 = 24 N = 48 M = 5j N_j A_j1 12 0.0458916723302 16 0.0888052602063 11 0.1618163898044 6 0.2318338801575 3 0.281681731186k % error in mu_k1 0.0000002 0.0000003 0.0000004 0.0000005 0.0000006 0.000567 II. EXPLOITING EQUIVALENT SPINS
The number of ways W ( N, I ) in which N equivalent spin-1/2 nuclei can be combined to give a resultant spin withangular momentum quantum number I is summarised in the following table: N I = 0
41 12 1 13 2 14 2 3 15 5 4 16 5 9 5 17 14 14 6 18 14 28 20 7 1(etc.)
The entries in this table are straightforward to generate on a computer using the recurrence relation W ( N, I ) = W ( N − , I + 1 / , I = 0 W ( N − , I − /
2) + W ( N − , I + 1 / , < I < N/ W ( N − , I − / , I = N/ , they are given explicitly by W ( N, I ) =
NN/ I ! (2 I + 1)( N/ I + 1) , and one can show that they satisfy X I W ( N, I )(2 I + 1) = 2 N . I.e., the 2 N states in the uncoupled representation | σ , . . . , σ N i (where σ I = ± / i -th nuclearspin on the z axis) can be combined to give the same number of states | I, M I i in the coupled representation, where I ranges from mod( N, / N/ M I ranges from − I to I in integer steps.It follows from the above table that a central spin problem with N = 4 equivalent spin-1/2 nuclei can be reduced to3 separate calculations, each of which involves a single nuclear spin with I = 0, 1, or 2 coupled to the central electronspin. The results of these calculations are simply multiplied by the weight factors W ( N, I ) = 2, 3, 1, and then addedtogether and divided by Z = 2 N +1 to obtain the central spin correlation tensor R αβ ( t ) = 1 Z tr h ˆ S α (0) ˆ S β ( t ) i . A central spin problem with one set of N = 4 equivalent spin-1/2 nuclei and another set of N = 3 equivalent spin-1/2nuclei can be reduced in the same way to 6 separate calculations, each of which involves pair of nuclear spins with( I , I ) = (0 , / W ( N , I ) W ( N , I ) = 4, 2, 6, 3, 2, 1, respectively. And so on. One can automate this procedure foran arbitrary number of sets of equivalent nuclei, and also generalise it to the case where the equivalent nuclei havespins other than 1/2.2 III. ADDITIONAL CONVERGENCE TESTS
Figure 1 of the paper shows the convergence of R zz ( t ) with increasing M for Model I with N = 49. The followingfigures show similar convergence tests for Model I with N = 99 and 999, and for Model II with N = 24 and 36. M = 20.030.060.090.12 M = 30 20 40 60 80 1000.000.030.060.090.12 M = 4 t / R zz ( t ) FIG. 1. Convergence of R zz ( t ) for Model I with N = 99, as a function of M . The solid black curve in each panel is the fullyconverged result obtained with M = 5, and the dashed red curve is the result obtained with the specified value of M . t / R zz ( t ) FIG. 2. Comparison of the results obtained with M = 2 (dashed red line) and M = 3 (solid black curve) for R zz ( t ) of Model Iwith N = 999. In this case, because of the large value of N , we did not have the computational resources to go up to M = 4.However, it is already clear from this comparison and the other convergence tests we have presented that the M = 3 resultsfor N = 999 are likely to be converged to graphical accuracy all the way out to t = 100 τ . M = 20.030.060.090.12 M = 30.030.060.090.12 M = 40 20 40 60 80 1000.000.030.060.090.12 M = 5 t / R zz ( t ) FIG. 3. Convergence of R zz ( t ) for Model II with N = 24, as a function of M . The solid black curve in each panel is the fullyconverged result obtained with M = 7, and the dashed red curve is the result obtained with the specified value of M . IV. LONG TIME CORRELATIONS
Figure 4 in the manuscript compares the Schulten-Wolynes (SW), improved semiclassical (SC), and quantummechanical (QM) correlation functions R zz ( t ) for Model I with N = 49, 99, and 999 nuclear spins. That figure wasplotted with linear axes ( R zz ( t ) versus t ) to emphasise that the SW theory misses the long time decay of the centralspin correlation function and the SC theory predicts too rapid a long time decay (see especially the insets in the N = 49 and N = 99 panels of the figure).The present Fig. 5 plots the same data on a semi-log plot ( R zz ( t ) versus log t ). This makes it clearer that, at leastfor N = 49 and 99, the SC and QM R zz ( t )’s have reached a plateau value by the time t = 100 τ . If the (cheaper) SCcalculation is extended to t = 200 τ , the computed R zz ( t ) remains at this plateau value, and since the SC and QMcalculations agree at t = 100 τ we expect that this would also be the case in the QM calculation. This suggests that,even with a finite number of nuclear spins in the central spin problem (here with a uniform distribution of hyperfinecoupling constants, and in the absence of an applied magnetic field), the central spin retains some information about itsinitial state in the long time limit. This is especially relevant to the quantum dot problem because it is a prerequisitefor being able to use a quantum dot as a qubit in a quantum computer (although of course in a real quantum dot thedipolar coupling between the nuclear spins – which we have ignored in the present calculations – will eventually playa role on a sufficiently long time scale).The SC and QM results for N = 999 in Fig. 5 have not yet reached their long-time limit at t = 100 τ , but we suspecton the basis of the N = 49 and 99 results that these R zz ( t )’s are also tending to a non-zero plateau value. We are notin a position to predict this value because the plateau values for N = 49 and 99 are both the same (0 . ± . M = 20.030.060.090.12 M = 30.030.060.090.12 M = 40 20 40 60 80 1000.000.030.060.090.12 M = 5 t / R zz ( t ) FIG. 4. Convergence of R zz ( t ) for Model II with N = 36, as a function of M . The solid black curve in each panel is the fullyconverged result obtained with M = 7, and the dashed red curve is the result obtained with the specified value of M . the nuclear spin precession by taking the limit as N → ∞ before the limit as t → ∞ , is R zz ( t → ∞ ) = 1 / ∼ . R xx ( t ) versus t ) to bring outthe long-time behaviour of the various (SW, SC, and QM) correlation functions. This plot is shown in the presentFig. 6. The SW theory is qualitatively wrong in this case – R xx ( t ) with a finite magnetic field in the z direction –in predicting a finite plateau value in the long time limit. The SC and QM curves have not converged to their longtime limits by the end of the plot ( t = 400 τ ), but are both seen to be oscillating around zero. When we extend the(cheaper) SC calculation to longer times, we find that the amplitude of the oscillation decays to zero, and we wouldexpect the same to be the case in the QM calculation. Combining this with the results for the other components ofthe spin correlation tensor (not shown in the figure), we find that a magnetic field of 1 / τ in the z direction causescomplete decoherence of the central spin in the xy plane ( R xx ( t → ∞ ) = R xy ( t → ∞ ) = R yy ( t → ∞ ) = 0), but notin the z direction ( R zz ( t → ∞ ) ∼ . ± . N = 49 SWQMSC R zz ( t ) N = 991 10 100 t /0.000.050.100.150.200.25 N = 999 FIG. 5. As in Fig. 4 of the paper, but with the data plotted on a semi-log plot ( R zz ( t ) versus log t ) to emphasise the long-timebehaviour of the Schulten-Wolynes (SW), semiclassical (SC) and quantum mechanical (QM) correlation functions. t /0.050.000.050.100.150.200.25 R xx ( t ) SWQMSC
FIG. 6. As in Fig. 5 of the paper, but with the data plotted on a semi-log plot ( R xx ( t ) versus log tt