Numerical Solver for the out-of-equilibrium time dependent Boltzmann Collision operator: Application to 2D materials
NNumerical Solver for the out-of-equilibrium Boltzmann Collision operator: 2Dsecond-order momentum discretisation
I. Wadgaonkar and M. Battiato ∗ Nanyang Technological University, 21 Nanyang Link, Singapore, Singapore
M. Wais
Technical University Wein, Vienna, Austria (Dated: February 17, 2021)Quantum Boltzmann equation (QBE) is a viable option to study strongly out-of-equilibrium ther-malization dynamics which are becoming increasingly critical for many novel physical applicationslike Ultrafast thermalization, Terahertz radiation etc. However its applicability is greatly limited bythe impractical scaling of the solution to the collision integral in QBE. In our previous work[1] wehad proposed a numerical solver for the solution of the collision integral in QBE and then improvedon it[2] to include second degree momentum discretisation and adaptive time stepping, thus makingit fully compatible to the standard numerical solvers for the transport part of QBE. The improvedsolver is numerically efficient and extremely robust against inherent numerical instabilities.Here weshowcase the applications of this improved solver to a simple 2D system,doped graphene, and analysethermalizations of the introduced out-of-equilibrium populations.
I. INTRODUCTIONII. THE QUANTUM BOLTZMANNEQUATION: COLLISION OPERATOR
The Collision operator in QBE is given by the Quan-tum Fokker-Planck equation. Using first order timedependent perturbation theory applied to a system ofweakly interacting quasiparticles, the Quantum Fokker-Planck equation calculates the time evolution of the mo-mentum resolved population, f n ( t, (cid:126) k ) (index n includes both, quasiparticle type and the quasiparticle band la-bel, while index (cid:126) k is for the crystal momentum vector)for each band for each quasiparticle involved in the scat-tering. Every combination of bands with an active inter-action (called as ’scattering channel’ from hence forth)contributes to the time evolution.For eg. the four leg fermionic scattering process and itstime reversed process depicted in fig.1 represent a singlescattering channel. The change in the population in eachleg involved in the scattering is obtained by the QuantumFokker-Planck equation for each leg separately. So thetime evolution of the population in the first leg or thefirst involved band, n , is, (cid:32) ∂f ∂t (cid:33) | n k > + | n k > ⇐⇒| n k > + | n k > = (cid:88) G (cid:90) (cid:90) (cid:90) V BZ d d (cid:126) k d d (cid:126) k d d (cid:126) k w e − e δ ( (cid:15) n ( (cid:126) k ) + (cid:15) n ( (cid:126) k ) − (cid:15) n ( (cid:126) k ) − (cid:15) n ( (cid:126) k )) δ ( (cid:126) k + (cid:126) k − (cid:126) k − (cid:126) k + G ) [(1 − f )(1 − f ) f f − f f (1 − f )(1 − f )] (cid:124) (cid:123)(cid:122) (cid:125) P (1)Here, f is short hand for f n ( t, (cid:126) k ) while (cid:15) n ( (cid:126) k ) isthe quasiparticle dispersion which is taken to be known. w e − e , which we call scattering amplitude, is, in principle,a function of the momenta (cid:126) k n of all the involved states inthe scattering process. Each scattering channel will havea corresponding value of w e − e . We remind that, in thiswork, (cid:15) n ( (cid:126) k ) and w e − e are inputs to eq.1 and their cal-culation is not addressed here. G is the reciprocal latticevector and the summation over G in eq.(1) is required ∗ [email protected] to account for Umklapp scatterings. The triple integralover the volume of brillouin zone, V BZ , in conjunctionwith the two Dirac deltas ensures that all combinations ofmomenta, which satisfy energy and momentum conserva-tion ( up to a reciprocal lattice vector G ), are accountedfor. d is the dimensionality of the system described andso, for our simulations of doped graphene we use d=2in eq.1 to get the appropriate expression of the Collisionoperator.Broadly speaking, eq.(1) represents scattering of elec-tron population into and out of a state at | n = n (cid:126) k = (cid:126) k > and so the phase factor, denoted as P , is sim-ply the total simultaneous probability for i) the states at a r X i v : . [ phy s i c s . c o m p - ph ] F e b FIG. 1. The 4-Leg scattering channel denoted as (cid:16) | n (cid:126) k > + | n (cid:126) k > ⇐⇒ | n (cid:126) k > + | n (cid:126) k > (cid:17) . An electronin state | n (cid:126) k > and another in state | n (cid:126) k > scatter witheach other and they end in state | n (cid:126) k > and state | n (cid:126) k > .b) the time reversed process | n (cid:126) k > and | n (cid:126) k > being unoccupied and the statesat | n (cid:126) k > and | n (cid:126) k > being occupied and ii) the timereversed process. Eq.(1) is completely general and canbe easily extended to any type of scatterings with differ-ent number of legs ( i.e. involved states) and betweendifferent quasi particles (eg. fermions and bosons), whilekeeping the mathematical structure unaltered. For thisreason we only show the numerical solution to eq.1.However, the solution of eq.(1) presents serious compu-tational challenges. i) As all the populations appearinginside the high dimensional integral in eq.1 are time de-pendent and known only at run time, the integral is aquartic operator. This results in an impractical scalingof the computational and storage cost. ii) The presenceof multiple Dirac deltas (depending on the dimensional-ity of the system considered), one of which has a highlynon-trivial form, make the integration domain highly dis- continuous. As such, the use of integration methods likeMonte Carlo Integration is unreliable. iii) Conservationof extensive quantities like particle number,momentumand energy is challenging and any non-conformity in thisregard leads to spurious results.These computational challenges have largely motivatedthe use of strong approximations and significantly limitedthe applicability of QBE to close-to-equilibrium scenar-ios. We have proposed a numerical algorithm[1] for thesolution of eq.1 and further extended it to include secondorder momentum discretisation and adaptive time step-ping [2] thereby making our numerical solver compatibleto the standard numerical solvers for the transport partof QBE which require a minimum of second order basisfunctions. While in our previous work we applied theimproved numerical solver to 1D materials, in this workwe showcase the applications of our solver to 2D ma-terials. In particular, we analyse thermalization of theintroduced excitations in doped graphene. A. Scattering Rates
Eq.(1) gives the change of the particular populationdistribution function f n ( t, (cid:126) k ) in a scattering process.However as outlined before, the numerical solution iscomputationally taxing and so it becomes important togenerate a visual aid,under suitable assumptions,into thebehaviour of eq.(1). It can be proved (refer [1]) thatfor a small and localized excitation in the state | n (cid:126) k > over the equilibrium distribution (appropriate Fermi-Dirac distribution in our case), the population decaysback exponentially to the equilibrium distribution. Theinverse of the time constant of this exponential decay isthe k-resolved scattering rate, λ n ( (cid:126) k ) which is given as( See [1] for details), λ n ( (cid:126) k ) = (cid:88) G (cid:90) (cid:90) (cid:90) V BZ d (cid:126) k d (cid:126) k d (cid:126) k w e − e δ ( (cid:15) n ( (cid:126) k ) + (cid:15) n ( (cid:126) k ) − (cid:15) n ( (cid:126) k ) − (cid:15) n ( (cid:126) k )) δ ( (cid:126) k + (cid:126) k − (cid:126) k − (cid:126) k + G )[(1 − f eq ) f eq f eq − f eq (1 − f eq )(1 − f eq )] (2)where f eqi ( shorthand for f eqn i ( (cid:126) k i ) ) is the relevantequilibrium distribution (Fermi-Dirac distribution in thiscase). Eq.2 appears very similar in structure to eq.1. Butthere is a very crucial difference. Notice that the phasefactor in eq.2 is composed of equilibrium population dis-tributions instead of the time dependent population dis-tributions. Therefore, the integral in eq.2 is not a quarticoperator and, as we will show in the next section, it canbe estimated at a fraction of the computational cost ofthe complete time propagation in eq.1. III. DISCRETIZATION
We highlight here only the changes made to our previ-ous solver (described in reference [1]) in order to increasethe order of convergence of momentum space. First, webriefly summarise the structure of the method which isindependent of the basis functions and then we describethe required modifications in the calculation of the scat-tering tensor elements owing to the second degree poly-nomial basis functions. Importantly, we focus on a sin-gle 2D band since we intend to apply ( for the sake ofthis study) the modified solver to excitations in dopedgraphene. However all the details are individually appli-cable to each band and to any number of bands. Also it is
FIG. 2. Sample chosen basis functions for a particular bandn. From left: Ψ A n ( (cid:126) k ) , Ψ A n ( (cid:126) k ) and Ψ B n ( (cid:126) k ) understood that the dimensionality of the system, d=2,implies that the crystal momenta, (cid:126) k n , are two dimen-sional vectors with components { k nx , k ny } in the carte-sian coordinate system.The band and momentum resolved population distri-bution, f n ( t, (cid:126) k ) , and dispersion relation, (cid:15) n ( (cid:126) k ) (which isassumed to be known before hand) are assumed to bedefined over a domain which can be a compact subset ofthe Brillouin zone. Note that each band can have its owndefined domain. This allows for a selective resolution ofthe Brillouin zone based on actual dynamics and can leadto large computational savings. A. General Structure of the method
First, we split the domain of the band in non-overlapping elements thereby forming a mesh. For ourchosen 2D system (doped graphene) the mesh consists ofequal rectangular elements. Now we project the solutionand the dispersion on a set of band specific momentumbasis functions Ψ Aan ( (cid:126) k ) . The basis functions are zero ev-erywhere in the domain except on the element identifiedby the index A and they are continuous everywhere ex-cept on the edges of the element A. Each element in the mesh can have several basis functions which are labelledby sub index a, such that the basis functions are linearlyindependent. In this study we assume that the basisfunctions are orthonormal: (cid:90) d (cid:126) k Ψ Bbn ( (cid:126) k ) Ψ Aan ( (cid:126) k ) = δ AB δ ab (3)where δ is the Kronecker delta. The method is equallyvalid for a non-orthonormal basis set with the RHS ineq.3 replaced by a mass matrix. A few sample basis func-tions from the basis set chosen for this study are shownin fig.2. The detailed derivation of the exact functionalform of the chosen basis set is given in Appendix A.Once the basis functions are defined, we can expressthe solution, f n ( t, (cid:126) k ) and the dispersion, (cid:15) n ( (cid:126) k ) as a linearcombination of these basis functions: f n ( t, (cid:126) k ) = (cid:88) Bb f Bbn ( t ) Ψ Bbn ( (cid:126) k ) (cid:15) n ( (cid:126) k ) = (cid:88) Bb (cid:15) Bbn Ψ Bbn ( (cid:126) k )1 = (cid:88) Bb Bbn Ψ Bbn ( (cid:126) k ) (4)Here we have, for later convenience, also written thediscretized representation of a constant function whichis equal to 1 on the whole k space. f Bbn ( t ) , (cid:15) Bbn and Bbn are the coefficients of the discretised representations ofthe population distribution function, dispersion and theconstant function respectively.Projecting eq.1 on the chosen orthonormal basis func-tions, using eq.4 and using the fact that the basis func-tions are non-zero only over a single element ( See ref-erence [1] for details) we get the final expression for thesemi-discretised form of the collision integral as, (cid:88) Aa (cid:48) df Aa (cid:48) n ( t ) dt = (cid:88) Aa (cid:88) Bb (cid:88) Cc (cid:88) Dd S a (cid:48) abcdAABCDn n n n (cid:32) (1 Aan − f Aan ( t ))(1 Bbn − f Bbn ( t )) f Ccn ( t ) f Ddn ( t ) − f Aan ( t ) f Bbn ( t )(1 Ccn − f Ccn ( t ))(1 Ddn − f Ddn ( t )) (cid:33) (5)with S a (cid:48) abcdAABCDn n n n = (cid:88) G (cid:90) An (cid:90) Bn (cid:90) Cn (cid:90) Dn w e − e d (cid:126) k d (cid:126) k d (cid:126) k d (cid:126) k Ψ Aa (cid:48) n ( (cid:126) k ) Ψ Aan ( (cid:126) k ) Ψ Bbn ( (cid:126) k ) Ψ Ccn ( (cid:126) k ) Ψ Ddn ( (cid:126) k ) δ ( (cid:126) k + (cid:126) k − (cid:126) k − (cid:126) k + G ) δ (cid:16) (cid:88) α (cid:15) Aαn Ψ Aαn ( (cid:126) k ) + (cid:88) β (cid:15) Bβn Ψ Bβn ( (cid:126) k ) − (cid:88) γ (cid:15) Cγn Ψ Cγn ( (cid:126) k ) − (cid:88) ξ (cid:15) Dξn Ψ Dξn ( (cid:126) k ) (cid:17) (6)Eq.5 is still semi-discrete as the time variable has not been discretised yet. We refer to S a (cid:48) abcdAABCDn n n n as the ’Scat-tering tensor’ as it practically contains all the informationabout the scattering. However, note that the Scatteringtensor calculation must be repeated for each scatteringchannel. In this study we consider only one band (whichmeans only one scattering channel) and hence we requireonly one calculation of the scattering tensor.Once the scattering tensor is calculated the time propa-gation of the population distribution function is straight-forward using eq.5. We highlight here that the integralsin eq.6 are no longer over the whole Brillouin zone butrather only on single elements owing to our choice of piecewise continuous polynomials as basis functions. This has far reaching consequences on the scaling of the overallcomputational cost. We study this aspect in further de-tail in sec. VI A.
1. Discretized form of Scattering rates
To obtain the discretised form of the scatteringrates, λ Aa (cid:48) n , we project the equation for the scatteringrates ,eq. 2, on the chosen basis functions and followa similar procedure as above (See ref.[1] for details) toobtain: λ Aa (cid:48) n = (cid:88) Aa (cid:88) Bb (cid:88) Cc (cid:88) Dd S a (cid:48) abcdAABCDn n n n (cid:32) Aan (cid:104) Bbn − f eqBbn (cid:105) f eqCcn f eqDdn − Aan f eqBbn (cid:104) Ccn − f eqCcn (cid:105)(cid:104) Ddn − f eqDdn (cid:105)(cid:33) (7) IV. CALCULATION OF THE SCATTERINGTENSOR ELEMENTS
In this section we detail the calculation of the scat-tering tensor elements using eq.6. Due to the presenceof multiple Dirac deltas the integration domain in eq.6is an extremely discontinuous hyper-surface and so astraightforward use of Monte Carlo Integration methodis impractical. However, our choice of polynomial ba-sis functions converts the expressions in both the Diracdeltas to simple polynomials as we have shown in eq.6,making it possible to analytically invert them. But theDirac delta for momentum conservation is unsymmet-rical as a result of which the final expression after an- alytical reduction of both the Dirac deltas is depen-dent on the order of variables we choose for the inver-sion. This disparity in the final expression for integra-tion is removed by effecting a mapping of variables as: k → x A ; k → x B ; k → − x C ; k → − x D + G . Thismapping brings the momentum Dirac delta in a com-pletely symmetric form as δ ( x A + x B + x C + x D ) . Nowwe invert the energy Dirac delta with respect to the firsttwo variables. Note that we can use the same integrationroutine for inverting the energy Dirac delta with respectto any couple of variables by simply altering the map-ping. Appendix B details the inversion of Dirac deltasand the derivation of the final structure of the expressionfor Monte Carlo integration which has the form as: Φ[ F [] , l Ax , h Ax , l Ay , h Ay , l Bx , h Bx , l By , h By , l Cx , h Cx , l Cy , h Cy , l Dx , h Dx , l Dy , h Dy , µ , µ A , µ A , µ A , µ A , µ A , µ B , µ B , µ B , µ B , µ B , µ C , µ C , µ C , µ C , µ C , µ D , µ D , µ D , µ D , µ D ]= (cid:90) h By l By dy B (cid:90) h Cx l Cx (cid:90) h Cy l Cy dx C dy C (cid:90) h Dx l Dx (cid:90) h Dy l Dy dx D dy D (cid:32) F [ x A + [ y B , x C , y C , x D , y D ] , y A [ y B , y C , y D ] , x B + [ y B , x C , y C , x D , y D ] , y B , x C , y C , x D , y D ] (cid:112) D [ y B , x C , y C , x D , y D ]Θ [ l Ax ,h Ax ,l Ay ,h Ay ] [ x A + [ y B , x C , y C , x D , y D ] , y A [ y B , y C , y D ]]Θ [ l Bx ,h Bx ,l By ,h By ] [ x B + [ y B , x C , y C , x D , y D ] , y B ]Θ [0 , ∞ ] [ D [ y B , x C , y C , x D , y D ]]+ F [ x A − [ y B , x C , y C , x D , y D ] , y A [ y B , y C , y D ] , x B − [ y B , x C , y C , x D , y D ] , y B , x C , y C , x D , y D ] (cid:112) D [ y B , x C , y C , x D , y D ]Θ [ l Ax ,h Ax ,l Ay ,h Ay ] [ x A − [ y B , x C , y C , x D , y D ] , y A [ y B , y C , y D ]]Θ [ l Bx ,h Bx ,l By ,h By ] [ x B − [ y B , x C , y C , x D , y D ] , y B ]Θ [0 , ∞ ] [ D [ y B , x C , y C , x D , y D ]] (cid:33) (8)where Θ [ ...,... ] [ ... ] represents the Heaviside function be- tween the edges along the x and y directions (i.e l Ax , h Ax or l Ay , h Ay or l Bx , h Bx ), of the elements corresponding tothe variables reduced (in this case x A or y A or x B ), Dis the discriminant of the quadratic equation obtained inthe reduction of energy Dirac delta and the basis func-tions and w e − e , in eq.(6), are grouped in F[...]. For amore detailed description of eq.8 refer to Appendix B.Eq.(8) is rid of Dirac deltas and hence is a smooth func-tion which has been converted from a 8 dimensional in-tegral ( eq.6) to a 5 dimensional integral. It can be in-tegrated, in principle, with the standard Monte Carlotechnique. V. TIME PROPAGATION
Once the scattering tensor is calculated, the popula-tion can be easily propagated in time, using eq.5, bycontracting the scattering tensor with the instantaneouspopulations. The discretisation in time for eq.5, how-ever, would require an intelligent choice of a suitable timepropagation scheme given that multiple time scales arenormally involved in the far-from-equilibrium thermali-sation dynamics. Following the conclusions of our previ-ous work[2], we choose to employ adaptive time steppingthrough DP853 numerical algorithm.
VI. NUMERICAL TESTS FOR THE CODE
In this section we outline the numerical tests conductedto confirm the consistency of the code.
A. Scaling of computational cost and storage costof the scattering tensor
In ref.[1] we show that with our numerical method, ow-ing to our choice of basis functions, the scattering tensorsize scales as ∼ N d − d , where d is the dimensionality ofthe system and N is the total number of elements in themesh (referred to as mesh resolution for convenience).Hence, for the 2D system considered in this study, thescattering tensor size should scale as ∼ N . . Similarlythe actual Wall time required for the calculation of thescattering tensor is another critical quantity of interestand we expect it to follow a similar scaling law. To verifythis, in fig.3 we plot the number of entries in the scatter-ing tensor and the computational time for the scatteringtensor against an increasing mesh resolution,N. As seenfrom the figure, the size of the scattering tensor and thecomputational time indeed scale as ∼ N . . B. Representation of the band structure
We had limited the order of polynomials to 2 whilechoosing our basis set. As a check to see if it ade-quately represents the band structure of doped graphene,
FIG. 3. Size of the scattering tensor and Computational timefor the scattering tensor vs number of mesh elementsFIG. 4. Band Structure of Graphene :Plotted from dispersionrelation(left) and Representation generated by the numericalcode (right) in the chosen basis functions. we compare the projected dispersion from eq.(4)as gener-ated by the numerical code with the actual dispersion ofdoped graphene. Fig.4 shows that the band structure ofdoped graphene is very closely represented including theDirac point (which is one of the advantages of using thechosen basis functions). The chosen cartesian mesh al-lowed for the representation of the band structure alongwith the dirac point by adjusting the mesh such that theDirac point lies exactly on a mesh node.
C. Conservation of particles, momentum andEnergy
The finite stochastic error present in the Monte CarloIntegration method breaks the inherent symmetries ofthe scattering tensor which are equivalent to particlenumber, momentum and energy conservation. Failing toconserve these quantities could lead to spurious resultsand hence this is an extremely critical issue. However,when an apposite construction of the Monte Carlo pointsis done, these conservations are automatically obeyed. Assuch, we analyzed the particle number, momentum andenergy conservation for our scattering tensor before ac-
TimeStep ∆ P (∆ K ) x (∆ K ) y ∆ E ∆ P ), Change in mo-mentum along X direction ( (∆ K ) x ), Change in momentumalong Y direction ( (∆ K ) y ) and Change in the total energy( ∆ E ) with time steps tually using it for time propagation of the population.Table I shows the change in particle number, momentumand energy over 4 time steps. Indeed, the numerical codepreserves these critical quantities to machine precision asseen from Table I. D. Is Fermi Dirac distribution a steady statesolution?
The Fermi Dirac distribution is the equilibrium dis-tribution for fermions and hence should be ideally thesteady state solution. Yet, since the numerical code dis-cretises the dispersion relation and the population, thesteady state solution for the numerical code would corre-spond to this discretised band structure and hence wouldbe slightly different than the input Fermi Dirac distri-bution. However, in principle, the thermalisation pro-cess is an exponential decay. So it is not obvious howlong we have to propagate to be sufficiently close to thethermal equilibrium. We start with the bandstructuredepicted in fig.4 and an initial Fermi Dirac distributionwith µ = 0 . eV and T=700 K and let this populationpropagate in time for 300 time steps with a time stepof dt=0.0001 to obtain a numerical steady state solu-tion. At this point we propagate the distribution for onemore time step and note the change in population. Nowwe introduce a small number of particles in the system,which corresponds to an excitation, and note the changein the population after 1 time step. Figure 5 depicts thissequence well.The change in population post excitation is more thanfive orders of magnitude greater than the change in pop-ulation after the numerically obtained distribution after300 time steps. Thus, it can be safely concluded that thenumerically obtained distribution after 300 time steps issufficiently close to the steady state solution. VII. RESULTS AND DISCUSSION
In this section we apply the code to a selected 2D sys-tem and study the thermalisations of the added out-of-equilibrium excitations. Since it is not our intention tostudy real systems, we only demonstrate the capabilitiesof the code considering excitations in doped graphene
FIG. 5. Check for confirming that Fermi Dirac distribution isthe steady state solution. a) Initial input Fermi-Dirac distri-bution b) Numerical Steady State solution c) Change in thepopulation from the numerical steady state solution after 1time step of dt=1e-4 d) Small excitation added to the numer-ical steady state solution e) Change in the population fromthe numerical steady state solution post addition of the exci-tation after 1 time step of dt=1e-4. Note the orders of magni-tude change in the population post addition of excitation ascompared to the change in the population after reaching thenumerical steady state solution. with a bandstructure shown in fig.4. We follow a method-ology as below:1.
Choose a mesh resolution:
In this study themesh resolution was kept to 20 elements in the xand y directions each.2.
Calculate the list of element combinationswith possible scatterings:
We traverse throughall the possible combinations of elements, corre-sponding to the scattering channel for which wewish to calculate the scattering tensor, to find alist of combinations where momentum and energyconservation can be satisfied.3.
Calculate the scattering tensor using thislist:
The list obtained in the previous step is usedto calculate the scattering tensor corresponding tothe scattering channel considered. This scatteringtensor will be used to time propagate the excita-tions and study the thermalisation characteristics.4.
Introduce an initial Fermi Dirac and letit stabilise to an equilibrium distributionas described in section VI.
The initial FermiDirac distribution is introduced at a chemicalpotential, µ = 0 . eV and Temperature, T = FIG. 6. Numerical steady state distribution after stabilizingthe input Fermi-Dirac distribution for 300 timesteps. K .We then let it stabilize to a numerical steadystate distribution as described in section VI. Nowwe have a numerical steady state solution as shownin figure 6.5. Introduce an excitation in this numericalsteady state solution and propagate in time .The introduced excitation is propagated in time us-ing DP853 time propagation scheme. We monitorthe change in population as it decays to a new equi-librium distribution and analyze the thermalizationdynamics.Initially we study small excitations over the equilib-rium Fermi-Dirac distribution before moving on to exci-tations which resemble the ones in real experiments, saylaser interactions with doped graphene.In principle,the scattering amplitude, w e − e , is a func-tion of the momenta, (cid:126) k n , of all the states involved in thescattering. However, following the conclusions of [3] thatthe impact of w e − e on the scattering integral is muchweaker than that of the phase space factor, denoted as P in eq.1, we approximate the value of w e − e to bea constant equal to 1 in all of the k-space. We do notintend to study real systems here and so the exact valueof the constant is irrelevant. However, as shown in Ap-pendix C, when the transferred momentum is small, inreal systems the change in population vanishes. In ournumerical code instead, the analytical treatment of Diracdeltas leads to a divergence of the joint density of statesand thereby a diverging scattering integral in such sce-narios. To rectify this, when the transferred momentumis less than or equal to the mesh width we choose a spe-cial functional form for the scattering amplitude, w e − e ,as given in Appendix C. FIG. 7. Equilibrium Scattering Rates for an electron in theband depicted in fig.4. Equal energy contours are also shownin red. Higher energies correspond to higher scattering rates.
A. Scattering Rates
Before studying the decay of excitations, we analyzethe scattering rates for the band structure depicted infig.4. Fig.7 shows the scattering rates juxtaposed withequal energy lines in red. Particles added at higher en-ergies should decay faster thereby giving higher scatter-ing rates. For the inverted cone band structure shown infig.4, energy increases radially from the center of the Bril-louin zone. Accordingly fig.7 depicts radially increasingscattering rates as expected.
B. Decay of small Excitations
As a first step towards the application of the numericalcode to 2D materials we run a few test cases with smallexcitations introduced to the numerical steady state dis-tribution in fig. 6 and then analyse their thermalisations.The band structure and the initial numerical steady statedistribution are both radially symmetric. Hence it wouldsuffice to study different excitations only in one quadrantalong the diagonal and along a line perpendicular to oneside of the domain. Fig. 8 shows the locations whereexcitations are introduced while the exact co-ordinateson the mesh for the introduced excitations are listed intable II.Figure 9 shows the decay of one of the excitations atthe locations depicted in Fig.8(right part). As the ther-malization progresses, the population distribution func-tion varies smoothly and relaxes to a new equilibriumdistribution. Fig.10 shows the change of the populationdistribution function with respect to energy on a sectionof the Brillouin zone passing through the introduced ex-citation,as the thermalization progresses.
FIG. 8. Introduced excitations with the initial FermiDiracdistribution (left) and Locations of the introduced excitationsin the Brillouin zone given by full circles (right). The dottedcircle at the center highlights the location of the initial Fermi-Dirac distribution
Since we do not consider Umklapp scatterings here,the final equilibrium distribution is not guaranteed to bea Fermi-Dirac distribution. Nevertheless, from fig.10, weextract two quantities of interest as below: β = − ∂f∂E (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) f =0 . µ = E | f =0 . The first term, β , is analogous to the K B T term in theFermi-Dirac equilibrium distribution and we use it hereas an indicator of the temperature. The second term µ is analogous to the chemical potential term in the Fermi-Dirac equilibrium distribution and we use it here as anindicator of the particle number. Together, the variationof β and µ can provide interesting insights into how thepopulation distribution function changes its shape.Fig.11 shows the variation of β and µ with time forthe same excitation considered in fig.9. At time t=0, thevalue of β corresponds to that of the numerical steadystate distribution, which is a Fermi-Dirac distribution,and then as the thermalization progresses β decreasescontinuously indicating an increase in the overall tem-perature. This was expected since the added excitationwas at higher energies and energy conservation dictatesthat the thermalised distribution should therefore be athigher energy or higher temperature than the numericalsteady state distribution. In combination with the varia-tion of β , the variation of µ suggests that at the start ofthe thermalization, the particles added at very high en-ergy gradually fall to relatively lower energy states. Theoccupancy in the proximity of the Dirac point is still rel-atively low. As the thermalisation progresses, the statesin the proximity of the Dirac point start getting com-pletely filled gradually as seen from the recovery of µ . Adecrease in the Gibbs free energy during the thermalisa-tion process in conjunction with the increased number ofparticles should lead to a decreased value of µ at the endof the thermalisation as corroborated by fig.11.Fig.12 shows the average change in population withtime at the peak of this introduced excitation. As seenfrom the figure, the thermalisation is, in fact, a non- trivial function of time and not a simple exponential.However, for comparison we show a double exponentialfit. Although it is immediately evident that a doubleexponential fit is still not sufficient and does not com-pletely represent the thermalisation as seen from figure12, nonetheless, it can still help to generate a quantitativeestimation of the time scales involved in the thermalisa-tion. τ represents the initial time scale of decay, which isgenerally smaller, while τ (cid:48) represents the thermalisationtime scale at later stages, which is generally longer.To visualise the effect of the initial energy of the in-troduced excitation on the thermalisation dynamics, werepeated the analysis for all the introduced excitationslisted in table II. Each excitation is introduced at a dif-ferent point in the Brillouin zone and hence over a dif-ferent initial population distribution. As a result, eachexcitation starts to decay from a different value of theinitial population distribution function and thermalisesto its own unique value of population distribution func-tion. To compare the thermalisation of all the intro-duced excitations on a common footing we normalize thevalue of the population distribution function at the co-ordinates of each excitation by its initial value just afterthe excitation is introduced and compare these ’normal-ized thermalisation profiles’ in fig.13. The excitationsaway from the center are at higher energies and so thescattering tensor elements should be larger in magnitudeleading to a faster decay of excitation. Indeed, fig.13shows that the decay is rapid as the excitations moveaway from the center i.e. the decay rates are in the or-der: (3,3)>(5,5)>(3,10)>(5,10)>(7,7)>(8,8) . Thisis further corroborated from the initial scattering rates( /τ ) in Table II. The initial scattering rates increase forexcitations further away from the center of the Brillouinzone. Comparing these initial scattering rates with thosein fig.2 we see a reasonable agreement. We also note thatthe radial symmetry of the problem is perfectly visible inthe thermalisation curves and the scattering rates.Finally, it is immediately evident from Fig.13 that mul-tiple time scales are involved in the thermalization pro-cess. Specifying a single time step value at the start ofRK4 routine would have been impractical and so our useof adaptive time stepping routine instead of the RK4 rou-tine is justified. C. Decay of real excitations
In reality, for doped graphene, when light excites elec-trons from the valence band to the conduction band, theexcitations will not manifest as an isolated bump like weintroduced to the numerical steady state distribution inthe last section. Rather, a better approximation for suchexcitations would be a Gaussian distribution spread overthe periphery of the cone at the energy level proportionalto the frequency of the light causing the excitation. Toanalyze the thermalization of such cases we simulated 2cases: One with a small excitation, introduced at a ra-
FIG. 9. Decay of an introduced excitation with time. Weshow here the change in the population distribution functiononly for the excitation at (3,10).FIG. 10. Variation of population distribution function for theexcitation at (3,10) with respect to energy as the thermalisa-tion progresses. dius of 6(1/nm) and 7 (1/nm) from the center of theBrillouin zone, and another with a heavy excitation at aradius of 6 (1/nm) from the center of the Brillouin zone.For both cases the population distribution function variessmoothly to a new equilibrium distribution as seen fromfig.14 and fig.18 respectively. As mentioned before thereis no UmKlapp scattering and hence the new equilibriumdistribution is not guaranteed to be a Fermi-Dirac distri-bution. However, similar to the case of isolated excita-
FIG. 11. Variation of β and µ for the excitation at (3,10) withtime.FIG. 12. Decay at the peak of the introduced excitationfrom fig.9 with time. A double exponential fit is also shownfor comparison. The constants are: a=0.1683; τ =1/822.6;b=0.009836; τ (cid:48) =1/159.Excitation No. Co-ordinates ofthe Excitation onmesh (X,Y) Scattering Rates (1 /τ ) /τ ) for each excitation. Excitations at higher energiesare expected to decay faster and accordingly the scatteringrates increase from excitation no.1 to no.6. FIG. 13. Thermalization of introduced excitations with time.For a direct comparison, we plot the normalised populationdistribution function at the peak of the introduced excitationvs time. The initial value of the population distribution func-tion, f i , is used for the normalisation. Excitations introducedfurther away from the center of Brillouin zone are at higherenergies and hence they thermalise relatively quickly.FIG. 14. Thermalization of realistic small excitation, at aradius of 6 (1/nm), with time. The total population distribu-tion with the introduced gaussian excitation thermalises to anew equilibrium distribution with time. tions considered previously, we can analyse the variationof the population distribution function with energy alongany section of the Brillouin zone (the excitation is radi-ally symmetric) as shown in fig.15.The variation of respective β and µ is shown in fig-ure 16 and fig.20 respectively. As in the previous case ofisolated excitations, beta decreases monotonically indi-cating an increase in temperature which was expected asbefore. The variation of µ is much steeper as seen fromthe figures. The higher number of particles in the intro-duced gaussian excitations as compared to the isolated FIG. 15. Energy resolved change in the distribution functionof the introduced small excitation in fig. ?? with timeFIG. 16. Variation of β and µ with time for the introducedsmall gaussian excitations at a radius of 6 (1/nm) and a radiusof 7 (1/nm). excitations cause an immediate occupation of the statesin the proximity of the Dirac point and a steeper decreaseof µ . Fig.16 also juxtaposes the difference in variation of β and µ for the same gaussian excitation introduced atdifferent initial energies. For the lower energy excita-tion at a radius of 6 (1/nm) the rise in temperature isrelatively slower and the final temperature is relativelyhigher in comparison to the higher energy excitation atradius of 7 (1/nm). Following the higher scattering ratesfor the higher energy excitation the decrease in µ is alsohigher compared to the lower energy excitation.In fig.17 and fig.21 we also present the decay with timeat the peak of the introduced gaussian excitations. Forcomparison we again fit a double exponential functionand generate a quantitative estimate of the time scalesinvolved in the thermalisation. We remind that the valueof the scattering amplitude, w e − e was chosen to be a con-stant equal to 1 in these simulations. Fitting the decay ofthe introduced excitation from our simulations (eg. fig.17or fig.21) to an actual decay data from experiments, wecan conclude on an appropriate value of the scatteringamplitude, w e − e thereby fixing the only free parameterin our method.1 FIG. 17. Decay at the peak of the introduced small excitation( radius=6 (1/nm)) with time. A double exponential fit is alsoshown for comparison. The decay constant can be inferredfrom the time constants, τ and τ (cid:48) in the exponential fit. Theconstants are: a=0.04162;b=0.06997; τ =1/1518; τ (cid:48) =1/1.823FIG. 18. Thermalization of realistic heavy excitation withtime VIII. CONCLUSIONS
In our previous work [1] we had proposed a solver forthe solution of the time dependent Boltzmann scatter-ing integral, with no approximations and with realisticband structures and matrix elements. We had furtherimproved our proposed numerical solver to include sec-ond order momentum discretisation and adaptive timestepping [2]. In this work, we showcased the applications
FIG. 19. Energy resolved change in the distribution functionof the introduced heavy excitation in fig. ?? with timeFIG. 20. Variation of β and µ with time for realistic heavyexcitationFIG. 21. Decay at the peak of the introduced heavy exci-tation with time. A double exponential fit is also shown forcomparison. The decay constant can be inferred from thetime constants, τ and τ (cid:48) in the exponential fit. The constantsare: a=0.04953;b=0.1551; τ =1/4525; τ (cid:48) =1/241.1 w e − e = 1 ),comparing ourresults to the actual experimental data we can arrive ata relevant value of scattering amplitude for the materialthereby fixing the only free parameter in our code andmaking it truly applicable to a host of novel physicalphenomenon concerning 2D materials. Appendix A: Basis functions
We need an orthonormal set of polynomials as basisfunctions. In 1D, we use normalized Lagrange Polynomi-als as the required orthonormal basis set. They have theform, Ψ Bbn ( k ) = P Bbn ( k ) , k ∈ E Bn = 0 , otherwise (A.1)where, E Bn denotes the element labelled B in band n and P Bbn ( k ) has the form, P B n = γ B n P B n = β B n k + γ B n P B n = α B n k + β B n k + γ B n (A.2)The six unknown coefficients are determined by imposingthe orthonormality condition.In 2D, we generate our basis set using the 1D LagrangePolynomials along the X and Y direction as follows:The 1D Lagrange polynomials along the X and the Y direction are: Ψ i = 1 √ ∆ i Ψ i = 2 √ √ ∆ i ( k i − C i )∆ i Ψ i = √ √ ∆ i (cid:32) (cid:16) k i − C i ∆ i (cid:17) − (cid:33) Where, i=x or y depending on if the Lagrange poly-nomials are taken along the X or the Y direction re-spectively. Accordingly, { ∆ x , ∆ y } are the widths of themesh elements in the X and Y directions respectively and{ C x , C y } are the centers of the mesh elements in the Xand the Y directions respectively. Taking the cartesianproduct of these 1D Lagrange polynomials and limitingthe highest order of the resulting polynomial to 2 sincewe want up to second order basis functions, we get anorthomormal basis set in 2D as, Ψ =Ψ x ∗ Ψ y = 1 (cid:112) ∆ x ∆ y Ψ =Ψ x ∗ Ψ y = 2 √ (cid:112) ∆ x ∆ y ( k x − C x )∆ x Ψ =Ψ x ∗ Ψ y = 2 √ (cid:112) ∆ x ∆ y ( k y − C y )∆ y Ψ =Ψ x ∗ Ψ y = 12∆ / x ∆ / y ( k x − C x )( k y − C y )Ψ =Ψ x ∗ Ψ y = √ (cid:112) ∆ x ∆ y (cid:32) (cid:16) k x − C x ∆ x (cid:17) − (cid:33) Ψ =Ψ x ∗ Ψ y = √ (cid:112) ∆ x ∆ y (cid:32) (cid:16) k y − C y ∆ y (cid:17) − (cid:33) Appendix B: Final form for Integration of theScattering tensor elements
The expressions inside the dirac deltas for momen-tum and energy conservation reduce to simple polyno-mials, as shown in eq.6, following our choice of basisfunctions. Hence both Dirac deltas can be analyticallyreduced to obtain a smooth final expression for integra-tion by Monte Carlo method. However, the dissymme-try in the expression of the momentum Dirac delta ineq.6 would result in different final expressions dependingupon the variables on which the deltas are reduced andthe order in which the the deltas are inverted. To re-move this disparity we effect a mapping of variables as : k → x A ; k → x B ; k → − x C ; k → − x D + G . Now eq.6takes the form,3 Φ[ F [] , l Ax , h Ax , l Ay , h Ay , l Bx , h Bx , l By , h By , l Cx , h Cx , l Cy , h Cy , l Dx , h Dx , l Dy , h Dy ,µ , µ A , µ A , µ A , µ A , µ A , µ B , µ B , µ B , µ B , µ B , µ C , µ C , µ C , µ C , µ C , µ D , µ D , µ D , µ D , µ D , ]= (cid:90) h Ax l Ax n (cid:90) h Ay l Ay n dx A dy A (cid:90) h Bx l Bx n (cid:90) h By l By n dx B dy B (cid:90) h Cx l Cx n (cid:90) h Cy l Cy n dx C dy C (cid:90) h Dx l Dx n (cid:90) h Dy l Dy n dx D dy D F [ x A , y A , x B , y B , x C , y c , x D , y D ] δ ( x A + x B + x C + x D ) δ ( y A + y B + y C + y D ) δ [ µ + µ A x A + µ A y A + µ A x A y A + µ A x A + µ A y A + µ B x B + µ B y B + µ B x B y B + µ B x B + µ B y B + µ C x C + µ C y C + µ C x C y C + µ C x C + µ C y C + µ D x D + µ D y D + µ D x D y D + µ D x D + µ D y D ] (B.1)where F[...] incorporates w e − e and the basis functions(Eg. Ψ Aan ( k ) ) in eq.6.Now that we have a symmetric expression in the mo-mentum Dirac delta,we can choose any variables for thereduction of the Dirac deltas. We choose to invert thethree Dirac deltas in eq.B.1 on the variables x A , y A and x B . Reducing x A and y A gives, x A [ x B , x C , x D ] = − x B − x C − x D , and y A [ y B , y C , y D ] = − y B − y C − y D .Substituting for x A and y A in the Dirac delta for en-ergy we get a quadratic equation in x B as, δ ( energy ) = αx B + β [ y B , x c , y C , x D , y D ] x B + γ [ y B , x c , y C , x D , y D ] (B.2)where, α = µ A + µ B β [ y B , x c , y C , x D , y D ] = µ B − µ A + 2 µ A ( x C + x D ) − µ A y A [ y B , y C , y D ] + µ B y B γ [ y B , x c , y C , x D , y D ] = µ + µ A y A [ y B , y C , y D ] + µ A ( y A [ y B , y C , y D ]) − µ A ( x C + x D ) y A [ y B , y C , y D ]+ µ B y B + µ B y B + 2 µ A x c x D + ( µ C − µ A ) x C + ( µ A + µ C ) x C + µ C y C + µ C y C + µ C x C y C + ( µ D − µ A ) x D + ( µ A + µ D ) x D + µ D y D + µ D y D + µ D x D y D Solving eq.B.2 gives two values of x B ( x B + and x B − )which will give two corresponding values of x A ( x A + and x A − ). Finally we substitute for x A , y A and x B in eq.B.1to get, Φ[ ... ] = (cid:90) h By l By n dy B (cid:90) h Cx l Cx n (cid:90) h Cy l Cy n dx C dy C (cid:90) h Dx l Dx n (cid:90) h Dy l Dy n dx D dy D (cid:32) F [ x A + [ y B , x C , y C , x D , y D ] , y A [ y B , y C , y D ] , x B + [ y B , x C , y C , x D , y D ] , y B , x C , y c , x D , y D ] (cid:112) D [ y B , x C , y C , x D , y D ]Θ [ l Ax ,h Ax ,l Ay ,h Ay ] [ x A + [ y B , x C , y C , x D , y D ] , y A [ y B , y C , y D ]]Θ [ l Bx ,h Bx ,l By ,h By ] [ x B + [ y B , x C , y C , x D , y D ] , y B ]Θ [0 , ∞ ] [ D [ y B , x C , y C , x D , y D ]]+ F [ x A − [ y B , x C , y C , x D , y D ] , y A [ y B , y C , y D ] , x B − [ y B , x C , y C , x D , y D ] , y B , x C , y c , x D , y D ] (cid:112) D [ y B , x C , y C , x D , y D ]Θ [ l Ax ,h Ax ,l Ay ,h Ay ] [ x A − [ y B , x C , y C , x D , y D ] , y A [ y B , y C , y D ]]Θ [ l Bx ,h Bx ,l By ,h By ] [ x B − [ y B , x C , y C , x D , y D ] , y B ]Θ [0 , ∞ ] [ D [ y B , x C , y C , x D , y D ]] (cid:33) (B.3)Here, Θ[ ... ] is the heaviside function between the limitsof the reduced variables (In this case, x A , y A and x B ) and D is the determinant of the quadratic equation in eq.B.24as, D = ( β [ y B , x c , y C , x D , y D ]) − αγ [ y B , x c , y C , x D , y D ] The expression in eq.B.3 is rid of Dirac deltas and henceis a smooth function which has been reduced from a 8dimensional integral to a 5 dimensional integral. We canuse standard Monte Carlo integration method to obtainthe integration value.
Appendix C: Inputs for the Functional form of theScattering Amplitude
Consider a four leg electron-electron scattering similarto the one shown in fig.1. We assume here that all theinvolved electron states in the scattering belong to thesame band. For a more general case for all the involvedelectron states belonging to different bands, the deriva-tion preserves the mathematical structure but requires aband index to be specified where appropriate. As suchwe drop the band index from now on.Let ˆ N k denote the number operator for the instanta-neous occupation number of a state labelled by k. We re- mind that the expectation value of this number operator, < ˆ N k > , becomes the population distribution function, f ( t, k ) , in the context of the Boltzmann Equation.From the Quantum Fokker-Planck equation for thisscattering it can be shown that the change of < ˆ N k > with time goes as: d < ˆ N k >dt ∼ [ ˆ N k , V int ] (C.1)where , V int is the number conserving interaction termfor the scattering and it is given as, V int = 12 (cid:88) k ,k ,k U k ,k ,k ,k a † k a † k a k a k [ ˆ N k , V int ] is the commutator of the number operator withthe interaction term. So eq.C.1 implies that when thenumber operator commutes with the interaction termthere is no change in the occupation of the state.For simplicity of calculation, let us consider the changein the occupation number of the state labelled by k . Thecorresponding commutator [ ˆ N k , V int ] , can be derived as: ˆ N k a † k a † k a k a k =( a † k + a † k ˆ N k ) a † k a k a k = a † k a † k a k a k + a † k ( a † k δ k k + a † k ˆ N k ) a k a k = a † k a † k a k a k + a † k a † k δ k k a k a k + a † k a † k ( − a k δ k k + a k ˆ N k ) a k = a † k a † k a k a k + a † k a † k δ k k a k a k − a † k a † k a k δ k k a k + a † k a † k a k ( − a k δ k k + a k ˆ N k )= a † k a † k a k a k + a † k a † k δ k k a k a k − a † k a † k a k δ k k a k − a † k a † k a k a k δ k k + a † k a † k a k a k ˆ N k So, [ ˆ N k , ˆ V int ] = a † k a † k a k a k + a † k a † k δ k k a k a k − a † k a † k a k δ k k a k − a † k a † k a k a k δ k k (C.2)Eq.C.2 presents two scenarios when the commutator iszero or in other words d< ˆ N k >dt = 0 .i) k = k Or k = k ii) k = k and k = k and k = k i.e. all the involvedstates are the same. This is a trivial case and we do notpursue it further.In case of scenario (i) the momentum Dirac delta dic-tates that the remaining two ’k’s must also be equal toeach other and hence the transferred momentum is zero.In real systems, the change in the occupation numberfor such scenarios must be zero. But in our numericalcode when we analytically invert the Dirac deltas, such situations could result in the divergence of the scatter-ing integral. We rectify this by introducing a functionalform of the scattering amplitude w e − e such that it varieswith the square root of the distance between the involved’k’s ( say, | (cid:126) k − (cid:126) k | Or | (cid:126) k − (cid:126) k | ) as long as the distanceis less than or equal to the width of the mesh elementand takes a constant value otherwise (The value of thisconstant for a given material can be benchmarked fromthe actual experimental data on thermalisation in thatmaterial). With this functional form of the scatteringamplitude, the scattering integral vanishes ( leading toan unchanged population) as the transferred momentumtends to 0. [1] ’Deterministic solver for the time-dependent far-from-equilibrium quantum Boltzmann equation’, Wais, Michael and Held, Karsten and Battiato, Marco, arXiv preprint5