Correlated quantum dynamics of graphene clusters

.


I. INTRODUCTION
Since its discovery [1], graphene has become an extremely rich arena for research, both from fundamental science as well as for practical applications [2,3].The electronic structure of graphene (and graphite) predates the experimental observation of graphene with a broad margin [4,5], where the initial efforts were based on tight-binding electronic structure theory.Since these early efforts many studies have been published (see e.g.Ref. [6]), with results of the electronic structure that do not deviate significantly from the early results [4,5].One of the most interesting aspects of the electronic structure of graphene is the linear dispersion relation around the so called Dirac point, K, at the Brillouin zone boundary, with its peculiar consequences for Klein tunneling [7,8] and half-integer quantum Hall effect [1].The unique electronic properties of graphene around the Fermi level has opened for possible applications in electronics and spintronics (see e.g.Refs.[9,10]).
The valence band states that have attracted most attention are the so called π and π * states, that represent occupied and unoccupied electron states of undoped graphene.These states are composed of p z orbitals centered at each C atom.These states are weakly bonding, and the strong chemical bonds of graphene come instead from sp 2 hybrids (composed of s, p x and P y orbitals of each C atom) that build up a strong network of σ bonds (see e.g.Ref. [11]).The energy bands corresponding to these σ bonds are however far below the Fermi level, and are from a transport point of view rather uninteresting.This produces a rather interesting scenario, where the basic electronic structure close to the Fermi level (a few electron volts on either side of the Fermi level) of the bipartite graphene system can be described by tightbinding theory with one orbital (p z ) per atomic site.We will utilize the simplicity of the orbital structure of the π and π * states in this paper, by investigating the electronic structure using tight-binding theory, including on-site correlations, as described by the Hubbard model (see Section II).This allows to study a typical many-body, model Hamiltonian, with the advantage of its ability to describe a real, physical system.
The dynamical properties of the electronic structure of graphene is the main focus of this investigation, and we have employed several approximations to do this, they are: the Hartree-Fock (HF) method, Exact Diagonalisation (ED) as well as the fermionic Truncated Wigner Approximation (fTWA) [12].Of these the ED method is exact but only tractable for small systems because of its exponential scaling, while the other two represent approximations described and analyzed below.It is noteworthy that the fTWA method is the least frequently investigated approximation, when it comes to electronic structure theory, and we will for this reason put special emphasis on this method.The basic equation describing the quantum dynamics here is the Liouville-von Neumann equation of the density matrix, and we detail below the different technical aspects of its solution, comparing in particular the time evolution of the occupation numbers and second-order correlation functions.

II. THE MICROSCOPIC HAMILTONIAN
In this work, we consider the Hubbard-Fermi Hamiltonian, written in the second quantization formalism: where j ij is the hopping interaction of p z orbitals between sites i and j, σ is the particle's spin, while u ij is the interaction between two particles of opposite spins situated on sites i and j.Furthermore, ĉ is an annihilation operator and ĉ † a creation operator.In this article we only numerically consider hopping between neighbour sites and with equal potentials: j ij = J if i and j are neighbours, and in addition we limit the work to a system with on-site interaction: In the present work we chose equal strengths of J and U .The Hamiltonian ( 1) is a good representation of the electronic structure of graphene around the Fermi level, according to the discussion of the introduction.Hence this Hamiltonian has been used in several instances to simulate the energy dispersion of a graphene layer (see e.g.Refs.[13] and [14]).

III. FERMIONIC TRUNCATED WIGNER APPROXIMATION
Phase-space representations are a family of methods that have already demonstrated their ability to model the dynamic of many-body bosonic systems.They work by mapping the system's density matrix to a quasi-probability density and the Liouville-von Neumann equation of the Hamiltonian to a corresponding density differential equation for the probability.More recently, phase-space methods have been adapted to model fermionic dynamics.They are especially useful for 2D and 3D systems for which DMRG (density matrix renormalisation group), and alike methods, are less successful.We focus here on the computational efficiency of one approximate phase-space representation, called the fermionic Truncated Wigner Approximation (fTWA), and apply it to the Fermi-Hubbard model Eq.(1).
We start by providing a brief introduction to fTWA.In phase-space representations, we formulate the problem via an expansion of the density operator ρ over an over-complete operator basis Λ(λ) [15][16][17]: where the expansion 'coefficients' W (λ) constitute a quasi-distribution over generalised complex phasespace variables λ.The Liouville-von Neumann equation that describes the density operator dynamic is then mapped into a partial differential equation (PDE) of the Wigner function W (for more details see Appendix A 3).In the Wigner-Weyl representation, phase-space variables are mapped to symmetrized operators, which leads to a PDE containing only oddorder derivatives.In particular, the PDE for the Wigner function has no diffusion term (second-order derivative).The time dependent quantum operators Ô(t) are mapped to their Weyl symbols O W and evaluated in the Heisenberg representation [18]: The Truncated Wigner Approximation (TWA) is the practical implementation of the Wigner method.It here relies on two approximations: The high-order derivative terms (third-order and above) of the PDE are truncated, and the initial density is chosen to represent the two first moments (average and covariance).We use a gaussian distribution, similar to what was done in Ref. [18].Hence the Wigner function dynamics is found by computing trajectories of deterministic differential equations whose initial conditions are drawn from a gaussian probability density.
Because of the anti-commuting property of fermionic operators, they cannot directly be represented by complex numbers, so to adapt the Wigner representation to fermions, we choose to map phasespace variables to bilinear operators [12,18]: Here, letters i and j label site indices and σ i labels the particle spin.In a system with constant particle number, we only use Êiσi jσj , and we call ρ iσi,jσj its corresponding complex phase-space variable.We also consider in this article that electrons will not flip their spin, which is a reasonable assumption since the spinorbit coupling in graphene is very weak [8].For this reason we can limit the representation to only samespin phase-space variables ρ ijσ .The observable values are recovered using their Weyl symbol and Eq.(3), for example, the occupation operator and the doublon operators are linked to statistical averages of phase-space variables: see Appendix A 3 for details.The line above the symbols in Eq. ( 5) denotes stochastic averages.The time evolution of the quantum system is then given by a first order PDE for W , which provides differential equations for ρ ijσ [18]: where the hopping, j jk , and Coulomb repulsion, u ki , are defined in Eq. ( 1).From a practical point of view, W (ρ, t) is represented by a set of independent realizations of Eq. ( 6), called trajectories, whose initial condition are distributed respecting the moments in (5).
If the initial condition is a thermal state in a diagonal basis with site occupations n iiσ , we can compute the first and second moments (mean and covariance) between phase-space variables to generate an initial gaussian density with the same moments, which gives: and if where ξ ijσ is from a complex normal distribution, with ξ jiσ = ξ * ijσ .Again, we outline details of the calculations further in Appendix A 3.
In this article the Hamiltonian (1) is timeindependent.However, time-dependence in e.g. an external potential is straightforward to implement in fTWA, and appears as explicit time-dependent terms in Eq. ( 6), just as they appear in the corresponding mean-field method (below).We have numerically evaluated fTWA against the other computational methods used in this article also for time-dependent Hamiltonians (not presented here).

A. Exact Diagonalisation and Hartree-Fock
For the evaluation of the ED dynamics, we compute step by step the general solution of the timedependent Schrödinger equation, using the Krylov subspace projection technique implemented in the expokit package [19].For details of these calculations we refer to Appendix B 1.
To obtain the mean-field dynamics, we use the Heisenberg equation of motions to compute the dynamic of the number operators nijσ .Then we map their factorised averages to the variables of the meanfield method nijσ → n ijσ , see Appendix A for details.We recover the same differential equations that are found for the phase-space variables ρ ijσ , as shown in [20].Without the initial noise (ξ ijσ = 0), Eqs. ( 6) and ( 7) provide the same results as the Hartree-Fock mean-field method.

IV. RESULTS
Below we present the results for two examples, we start with a small systems where a comparison between different theoretical methods can be made (Exact Diagonalization, fTWA and Hartree-Fock), then we study a larger system on which the Exact Diagonalisation method is unable to give a result.

A. Small graphene systems
We first study the accuracy of the fTWA method on a small graphene-like few-body system.The system is composed of ten sites organised in two hexagonal cells, see Fig. 1, which makes it large enough to be interesting and small enough to compute numerical solutions with the Exact Diagonalisation method.The system is assumed to be electronically half-filled, which means that there are as many particles as there are sites.We also consider a system with equal amount of spin-up as spin-down electrons.For the initial condition, we choose the p-particles Fock-space vector with the largest overlap to the ground state, where each site is filled with either a spin-up or spin-down particle, see Fig. 1.We motivate this choice further in Appendix B 2. The time evolution with the fTWA method has been computed with 10 5 trajectories, and until t = 5, a choice that was made so that one can see when all methods considered here start to deviate from each other.

B. Evolution of occupations, first-order moments
The first measure used to compare the methods is the time evolution of site occupations, n iiσ (t).We show in Fig. 2 the occupation of spin-up particles on sites 1 and 5, respectively n 11↑ and n 55↑ .We see that the occupations computed with the Hartree-Fock method deviate from the exact solution (the ED method) at t ≃ 1 and that the occupations computed with fTWA starts to deviate later, at t ≃ 3. We stress that the only computational cost to achieve this improvement is that the calculation is repeated, in parallel, for the different trajectories.Note that the computation times of both HF and fTWA methods scale quadratically with the system size, and the fTWA method is here correct for approximately three times longer.Here, using 10 5 trajectories in the fTWA method, the statistical error is negligible, hidden behind the width of the line.Hence, the deviation comes entirely from the truncation possible from the approximation of the initial density [18] in the formalism of the method, and is not due to statistical uncertainty.We also note an improvement on the long term dynamics.When t → ∞ the occupation stabilizes at n iiσ ≃ 0.5 for all sites and spins (data not shown).This result is recovered by the fTWA method, but not by the HF one that oscillates uncontrollably.

C. Evolution of correlations, second-order moments
One advantage of phase-space methods, like fTWA, is their ability to give information on the dynamics of higher-order moments, like second-order correlation functions, even for large system.The correlation function between two sites i, j and spins σ 1 , σ 2 , here denoted g iσ1,jσ2 , can be seen as the effect the presence of a spin-σ 1 particle on site i has on the probability to have a spin-σ 2 particle on site j.The explicit formulas are, in the Schrödinger picture for ED: and in the Heisenberg picture with the fTWA phasespace variables: iσ1,jσ2 (t) = For the mean-field Hartree-Fock approximation we use Wick's theorem, see Eq. (A4) in Appendix A 2. In Fig. 3, we plot the correlation between the spinup particle on site 1 and the spin-up particles on the two neighboring sites, 2 and 4, see Fig. 1 for the geometry.The correlation computed with fTWA is essentially the same as that of ED, until approximately t ≃ 2 where it starts to deviate visibly.This happens earlier than for the occupation, shown in Fig. 2. At t ≃ 0, we observe that the correlation functions computed with fTWA does not tend precisely to its theoretical value, Eq. (C7), of g 1↑,2↑ = 0.5.The quantitative difficulties fTWA has to compute correlations involving initial empty sites are known for bosons [21], the calculations made in Appendix C expose the same problem for the present Hamiltonian.In principle, 1↑,2↑ and g 1↑,4↑ .The different curves represent the same methods and parameter values as described in Fig. 2.
these results can be improved by adding more trajectories or using projection methods.
In Fig. 4, we plot the correlations between the spinup particle on site 1 and the spin-up particles on distant sites, 6 and 10, see again Fig. 1 for the geometry.We observe that the fTWA results starts to deviate from data obtained by ED, at a time between t ≃ 2 and t ≃ 3. We also observe in Fig. 4 that the time for a correlation to appear between two sites, i.e. when g (2)  deviate from unity, depends on the increasing distance between those two sites.
1↑,2↑ , for systems of different size, n = 4, 6, 10, 198.The horizontal (black) lines represent the limit (t → ∞) values of correlation function gathered in table I, we recover the values of formula (C11) for all cases.
In the case of different spins correlations, like e.g.g i↑,i↓ , the second term in Eq. (A4) is zero.As a consequence HF gives a constant, g (2) ≡ 1, and cannot be used for comparisons.However, fTWA gives accurate results for short times, as we have explored numerically in comparisons with ED for small systems.

D. Large-time correlation functions
To investigate if the fTWA is also able to model the long term values of correlation functions, we computed the fTWA results of g 1↑,2↑ for systems of different size, n = 4, 6, 10, 198 sites, until t = 50 (note that we did not plot other correlations for clarity, but they all tend to the same limits).The results are shown in Fig. 5, and we see that the correlation tends to the limits that correspond to equal probability for particles, see the explanation and examples in Appendix C 3. Correlations between all the other pairs of different sites and same-spin particles have the same limit.For n = 4, 6, 10, the fTWA results were compared with Exact Diagonalisation.

V. FTWA ON A LARGE SYSTEM
Now that we have shown the ability of fTWA to efficiently model small Fermi-Hubbard systems, by a direct comparison to data from Exact Diagonalization, we study a substantially larger system, that is far out (Color online) Four different frames from the dynamics of the magnetic moment mi = nii ↑ − nii ↓ for times t = 0, 1, 2, 3. Sites in red have mi = 1, one spin-up particle, sites in blue have mi = −1, one spin-down particle.Sites in purple have mi = 0.For larger times we see that the magnetic moments all stabilize to mi ≃ 0, when the system tends to n ii↑ = n ii↓ ≃ 0.5.
of reach for Exact Diagonalisation.We compute the dynamics of fermions described by the Hamiltonian of Eq. ( 1) on a graphene-like system with a honeycomb structure involving 198 sites, see the overall geometry in Fig. 6.In this figure the value of the magnetisation, m i = n ii↑ − n ii↓ , is shown for four time frames of the calculation.At t = 0, the sites are filled with either a spin-up electron or a spin-down electron, represented by red and blue circles in Fig. 6, respectively.This represents a starting state similar to that considered in Fig. 1.For the subsequent times, the dynamics is such that the magnetisation at each site approaches zero (n ii↑ ≃ n ii↓ ≃ 0.5), representing an equal occupation of spin-up and spin-down electrons.This is seen most clearly in Fig. 6 by the purple color of all sites at t = 3.This result is consistent with experimental data of graphene at equilibrium conditions, that are known to reflect an equal occupation of spin-up and spin-down electrons [2].
For a more accurate view of the dynamics, we have plotted in Fig. 7 the occupation of spin-up electrons on sites 36 and 48, n 36↑ and n 48↑ , i.e. from sites in the middle of the system.Within the short time of the simulation, effects of the boundary of the 198 atom cluster play little role, and the system appears infinite, as reflected in the visible symmetry between n 36↑ and n 48↑ in Fig. 7.This symmetry breaks at later times, around t ≃ 4. Note that the data in Fig. 7 contains results from HF and fTWA calculations and that for short simulation times (t ≤ 3) we find similar occupations for the two approaches, in contrast to the example in Fig. 2 that considered a smaller system.The size of the system gives the site interactions a more predominant role, which we expect makes HF and fTWA better approximations of the real particle dynamics, but clearly ED is out of reach for comparisons.Also, for longer simulation times, the HF method gives highly oscillatory results, as shown in Fig. 7, which is not the case for the fTWA method.
Similar to the smaller system, we also study for the 198 atom system the correlation functions between different sites.We picked a central site, numbered 36, to be a representative one (see the arrow in the t = 0 subplot in Fig. 6), and followed the correlation functions with one of its nearest neighbours as well as with further distant sites, see the geometry in Fig. 8.In Fig. 9 we show the correlation functions of spin-up particles with a near neighbour site, g 36↑,45↑ (t).The correlation for the three neighbours are initially very similar because of local symmetries in the large system, where edge effects of the cluster play a lesser role.We recognise the short-time limit (t → 0), where g (2) (0) = 2/3, because each site has three neighbours, see the derivation in Appendix C, and the deviation from the initial value that we saw in the smaller system, in Fig. 3.For longer time scales the data in Fig. 9 approach a value close to one, that only depends on the total number of sites, see Appendix C 3. Again one may note large oscillations with the HF method for longer times.In Fig. 10 have plotted the correlation functions between the site numbered 36 and its neighbours at longer distance, see Fig. 8 for the geometry.Note that our choice of sites is such that they are initially filled with electrons of the same spin orientation, which means that at t = 0 the correlation function is one.We can observe from Fig. 10 that the further away two sites are, the longer it takes before g (2) starts to deviate from unity.

VI. CONCLUSION
In this work we have studied the quantum dynamics of the electronic structure of graphene-like systems, using an electronic Hamiltonian that allows for hopping and one-site Coulomb repulsion.The analysis is focused on the electron states close to the Fermi level, and are hence limited to p z orbitals of spin-up or spin-down character centered on each site of a honeycomb lattice site.This allows to study the dynamics of an electronic Hamiltonian that includes the minimum interactions to represent a realistic system, i.e. the electron hopping and on-site Coulomb repulsion.
We have in this investigation compared three methods with which to solve the time evolution of the electronic system; the Exact Diagonalization technique, the Hartree-Fock (HF) approximaand the fermionic Truncated Wigner Approximation (fTWA).In comparing the three approaches for smaller graphene-like systems we conclude that fTWA reproduces the results of Exact Diagonaliza-  6, at t = 0.The sites in red starts with a spin-up particle, e.g.n 36↑ (0) = 1, and the blue ones starts with a spin-down particle, e.g.n 48↓ (0) = 1.We will follow the correlation function of spin-up particles between site 36 and one of the three closest neighbours (blue lines) in Fig. 9. Then between site 36 and three distant sites (red lines) in Fig. 10.tion, for significantly longer times compared to HF, and for this reason we have focused on fTWA for larger systems.Previous works of fTWA have focused on long-range interactions [12], but mean-field and phasespace representation methods have larger difficulties to model on-site interactions because of the more pre-dominant role of quantum effects.Under those con- ditions, fTWA demonstrates a net improvement over mean-field methods.As shown here, the evolution of site occupations agrees well with exact results, for a period that is three times longer for fTWA compared to HF (see Fig. 2), and its long-time behaviour is also quantitatively recovered.The second-order correlation functions are also found to be well approximated, both on short-time dynamics (Figs. 3, 4) and longtime dynamics (Fig. 5).When comparing large and small systems, the results here are consistent with previous results; that smaller systems exhibit larger fluctuations in e.g. the correlation function, compared to larger ones.The improvements of fTWA over mean-field methods come with an acceptable computational cost.A fTWA computation scales as O(n 2 ) with the number of sites n, similar to the HF computation.It needs however multiple repetitions of computations (trajectories) to average upon.This cost is manageable on a single computer for the systems studied here and embarrassingly parallelizable for larger systems.
Which gives the differential equation for the singleparticle density matrix: (A5) In the Hamiltonian (1), no term allows the particle's spin to flip, (j αβ = j ab δ σaσ b ) so variables representing spin-flip are always null and are neglected.Moreover, the electron-electron interaction is only between opposite-spin particles, u αβ = u ab δ σa σb .We denote n ijσ = n , and Eq.(A5) then becomes These are the Hartree-Fock differential equations, as presented in Ref. [20].

Fermionic Truncated Wigner approximation
The fTWA is one of the phase-space representation methods, where we use a distribution to describe the electrons density matrix and map operators averages to the distribution's moments.For the fTWA in systems with constant number of particles, we can choose the symmetrically ordered one-body operator Êα β we defined in Eq. ( 4) and map it to the phase-space variables ρ αβ .First-order and second-order moments of the phase-space variable distribution are linked respectively to one-body and two-body operators: From those equations we can find the links between first-and second-order stochastic averages of phasespace variables and quantum operators in e.g.Eq. (5).
To obtain the equations of motion of the phase-space distribution, we follow the work of Polkovnikov [18,23].Here a Jordan-Schwinger mapping from the fermionic operators Ê to pair-bosonic operators was introduced, then the bosonic Truncated Wigner Approximation formalism can be applied.It was shown that the equation of motion for the phasespace variables ρ αβ are determined by the Poisson brackets with f being the structure constants of the bilinear operators; (A9) In the expression above, H W is the Hamiltonian in the ρ αβ variables: (A10) which gives us the differential equation (6).
We may not be able to reproduce all exact initial Wigner function with all high-order moments, but within the accuracy of the truncation, the two first moments are sufficient.So we approximate W (ρ, 0) with a gaussian distribution having the same first-and second-order moments (mean and covariance).From the relations (A7) we compute the covariance between phase-space variables: (A12) Hence we arrive at Eq. ( 7) for the random starting point of trajectories.

Analysis of the ground state
As initial condition for the dynamics, we are looking for a position eigenstate that has the largest overlap with the ground state of the Hamiltonian (1).In Fig. 11, we have plotted a representation of the ground state |Ψ gs = i b i |ψ i .On the x-axis are the indices of the N-particle Hilbert basis vectors and on the yaxis are the coordinates of the ground state in this basis, b i .We find two vectors that stand out with b i ≃ −0.05, the first one is the initial condition we have chosen, see Fig. 1, the other one is its spin symmetric state.
When we let the dynamics evolve, we have found that the density on all sites tends to 0.5 for large times.The corresponding t → ∞ wave-function is a broad distribution of all N-particle Hilbert space vectors, |Ψ (t) → i b i (∞) |ψ i with |b i (∞)| 2 ≃ 1/N .Before we converge the noises η to 0, the leading terms in the fraction (C1) are in O(ǫ).They determine the value of g (2) (0) when the number of trajectories is too small.Hence the difficulties to have an accurate value of g (2) (0), see Fig. 3, and the need for more trajectories.

Figure 1 .
Figure 1.(Color online) Illustration of one of the two dominating pure states of the ground state of the Hamiltonian (Eqn.1).Here we consider a 10 sites Fermi-Hubbard model with half filling.The geometry is composed of two perfect hexagons of equal side lengths.This state was used as initial condition for the quantum dynamics calculations.In blue, the spin-down sites, in red the spin-up sites.The second most dominating state is the spin symmetry of this one, i.e. a state where all electrons have flipped their spin (and blue and red colours have been interchanged in the figure).

Figure 2 .
Figure2.(Color online) Dynamics of site occupations for spin-up particles on sites 1 (n11) and 5 (n55).The solid (blue) curves show the fTWA dynamics, the dashed (black) ones show the exact diagonalisation dynamics, the dasheddotted (red) ones the Hartree-Fock dynamics.For the occupation, the HF dynamic starts to deviate from the ED at t ≃ 1, while the fTWA dynamic starts to deviate from the ED solution at around t ≃ 3.

Figure 3 .
Figure 3. (Color online) Dynamics of correlation functions between neighbour sites, here g

Figure 4 .
Figure 4. (Color online) Dynamics of correlation functions between distant sites, here g

Figure 5 .
Figure 5. (Color online) Long-time fTWA dynamics for correlations between spin-up particles on sites 1 and 2, g

Figure 6 .
Figure 6.(Color online) Four different frames from the dynamics of the magnetic moment mi = nii ↑ − nii ↓ for times t = 0, 1, 2, 3. Sites in red have mi = 1, one spin-up particle, sites in blue have mi = −1, one spin-down particle.Sites in purple have mi = 0.For larger times we see that the magnetic moments all stabilize to mi ≃ 0, when the system tends to n ii↑ = n ii↓ ≃ 0.5.

Figure 7 .
Figure 7. (Color online) Dynamics of site occupation for spin-up particles on sites 36 and 48 for a 198 atom cluster.Data obtained by HF given by dashed red curves and fTWA full blue curves.

Figure 8 .
Figure 8. (Color online) Zoom in on the 198-sites Fermi-Hubbard system, of Fig.6, at t = 0.The sites in red starts with a spin-up particle, e.g.n 36↑ (0) = 1, and the blue ones starts with a spin-down particle, e.g.n 48↓ (0) = 1.We will follow the correlation function of spin-up particles between site 36 and one of the three closest neighbours (blue lines) in Fig.9.Then between site 36 and three distant sites (red lines) in Fig.10.

Figure 9 .
Figure 9. (Color online) Dynamics of correlation functions between neighbour sites, here the site 36 and a neighbour, site 45.Data obtained by HF given by dashed red curves and fTWA full blue curves.

Figure 11 .
Figure 11.Ground state representation of the Fermi-Hubbard system with 10 sites, 5 spin-up and 5 spin-down particles.If |Ψgs = i bi |ψi , the x-axis represents the Nparticle Hilbert space vector indices i in a given order, and the y-axis represents the value of the coefficients of those basis vectors, bi.We have chosen the first basis vector with the larger absolute coordinate, bi ≃ −0.05, as initial condition for the dynamics reported in the main text.