Plunging in the Dirac sea using graphene quantum dots

The dynamics of low energy charge carriers in a graphene quantum dot subjected to a time-dependent local field is investigated numerically. In particular, we study a configuration where a Coulomb electric field is provided by an ion traversing the graphene sample. A Galerkin-like numerical scheme is introduced to solve the massless Dirac equation describing charge carriers subjected to space- and time-dependent electromagnetic potentials and is used to evaluate the field induced interband transitions. It is demonstrated that as the ion goes through graphene, electron-hole pairs are generated dynamically via the adiabatic pair creation mechanism around avoided crossings, similar to electron-positron pair generation in low energy heavy ion collisions.


I. INTRODUCTION
In the last two decades, Dirac materials have received an unprecedented amount of attention because they have very interesting electrical, mechanical and thermal properties [1].Graphene in particular, a 2-D array of carbon atoms arranged on a honeycomb lattice, is the quintessential Dirac material.Close to its Dirac points at a momentum |K ± | ≈ 1.7 Å−1 ≈ 12.3 eV/v F , the dispersion relation is linear and relativistic-like, being given by E = v F |p|+O(|p|/|K ± |), where v F ≈ 1.093×10 6 m/s is the Fermi velocity and p is the momentum deviation from K ± [2].In addition, the dynamics of charge carriers are described by a massless 2-D Dirac equation, as long as |p| |K ± |, confirming the Dirac material nature of graphene [3].
It was realized early that owing to this relativistic-like quantum behavior, the dynamics of charge carriers in graphene would be similar to the one of relativistic electrons and thus, could be used to simulate or investigate quantum electrodynamics (QED) processes [4,5].The Klein paradox [6,7], the Zitterbewegung [8], atomic collapse [9,10], the Schwinger process [11][12][13][14][15][16][17][18] and the Breit-Wheeler process [19] have been considered from this point of view.Therefore, whilst there may exists some qualitative difference between graphene and QED [5], this material provides a link between condensed matter and high energy physics.
In this article, electron-hole pair production in ion bombarded graphene is proposed as an analogue to QED electron-positron pair production in low energy heavy ion collisions (HIC), in the same spirit as these previous studies.The configuration considered is shown in Fig. 1.As the ion traverses the finite graphene sample (graphene quantum dot), it interacts with charge carriers in the valence band and has a certain probability to excite them in the conduction band, thus generating an electron-hole pair.In the following, we argue that under certain conditions, pair production occurs via the same mechanism as in HIC whereby some positive energy states "plunge" into the negative energy continuum (the Dirac sea).
The dynamics of charge carriers in graphene subjected to electromagnetic fields has gained momentum in the last few years, with the potential of controlling the electron dynamics in graphene-based devices [20][21][22][23].Many theoretical and experimental investigations have focused on the interaction of graphene with a (homogeneous) laser field [24][25][26][27][28].The regime of strong fields, whereby multi-photonic effects, the holy grail of atomic physics, start to be important and may lead to new physical phenomena, have also been investigated [14][15][16][17][18][29][30][31][32].For inhomogeneous fields, some work has also been performed.For example, some experiments [33] and numerical simulations using either time-dependent density functional theory (TDDFT) simulations [34][35][36][37], nonequilibrium Green functions methods [38] and molecular dynamics [39], have investigated ion bombardment of graphene samples.The main objectives of these studies was to quantify the amount of energy transferred from arXiv:2009.07648v1 [cond-mat.mes-hall]16 Sep 2020 the projectile to the target, the stopping power, and understand the formation of defects in the atomic structure as the charge traverses the honeycomb lattice.Another example is Ref. [40], where the current generated by highly charged ions (with Z = 20−30, where Z is the ion electric charge) have been estimated experimentally and using DFT calculations.The spatio-temporal dynamics of charge carriers under local optical excitation has been considered in Ref. [41].Nonetheless, the dynamics of electrons and holes subjected to localized inhomogeneous fields is not as well known as in the homogeneous case.One of the goals of this article is to fill this gap by investigating electron-hole pair production from a dynamical and localized Coulomb-like potential.
To study the main features of this dynamical and nonperturbative phenomenon, extensive numerical simulations are performed to evaluate the pair production rate and the electron density generated by the passing ion.Assuming that certain conditions (given below) are fulfilled, the charge carriers are modeled by a simple Dirac equation, allowing for a connection with QED processes and HIC.To solve this space-and time-dependent Dirac equation, a Galerkin-like numerical scheme is introduced, similar to the ones developed in Refs.[42][43][44].Finally, we discuss the challenges for observing this phenomenon experimentally.
This article is separated as follows.In section II, the main process for particle-hole production in graphene and its analogy with electron-positron pair creation are presented.Section III is devoted to the physical model for charge carriers.The numerical method is introduced in Section IV while numerical results can be found in Section V. We conclude in Section VI.Natural units where = c = 1 are used throughout the article.

II. PLUNGING IN THE DIRAC SEA: GRAPHENE AND HIC
In the Dirac sea interpretation, pair production is generated by an external field that induces transitions between the negative and positive energy states [45,46].In condensed matter systems, this corresponds to interband transitions between the valence and conduction bands.In the second-quantization formalism, these transitions between the negative and the positive continua are actually responsible for the generation of electron-hole or electron-positron pairs [45].
More specifically, the mechanism considered in this article is adiabatic pair creation (APC) [47] (sometimes called spontaneous pair creation [46]), a dynamical QED process where the external electromagnetic field potential has some specific characteristics: 1. Time-dependent potential well that can localize the electron (form bound states or resonances) at some given time, 2. Supercritical for some time interval, where the pos- itive bound states energies dive into the negative energy continuum, 3. Adiabatic and non-perturbative time evolution.
The single-particle state of electrons subjected to external electromagnetic fields of this type is depicted in Fig. 2.This figure also includes an energy gap, which is induced by the finite size of the graphene sample (see the end of Sec.IV A 1).The lowest energy states of the positive energy continuum (valence band) are shifted towards the negative energy states (conduction band), or Dirac sea.When electronic states dive into the Dirac sea, they cross with negative energy states and interband transitions can occur, resulting in the production of electron-hole pairs.Because initially, the negative energy states are filled (Dirac sea or valence band), the transitions always occur from the negative to the positive energy states.This physical interpretation is supported by the second-quantized formulation, as shown in Sec.

III.
In QED, the Coulomb field becomes supercritical when Z ∼ 137.Ions with such high charges do not exist in nature, so the strategy used in HIC is to collide two highly charged ions (typically fully stripped Uranium atoms with Z = 92).As the two ions approach each other, the bound state energies of the combined system decrease, up to a critical distance (at approximately R cr ∼ 35 fm [46]) where the lower energy state reaches E = −mc 2 and where the field becomes supercritical.At this point, the lower bound state starts "plunging" in the Dirac sea, i.e. it enters the negative energy continuum.This allows for transitions between negative and positive energy states and thus, allows for the generation of electron-positron pairs.After the collision, the bound states return to the positive energy continuum.Although this pair creation mechanism has been predicted in the 70's [48,49], it still eludes an experimental observation [46,50,51], in part because the ion "sticking" time is very short yielding a small number of pairs.
In the considered configuration (see Fig. 1), the field is provided by an ion passing through a graphene sample.As the ion passes, the lower energy states of the conduction band are pulled down and fall within the valence band, allowing for interband transitions.This occurs when β > 1 2 , i.e. for charges of value Z ∼ 1 [9,52,53].Finally, as the ion gets further, we recover the free dynamics but electron-hole pairs have been created.This phenomenon is the dynamical analog of the atomic collapse resonances observed experimentally from artificial nuclei on graphene [10].
The main differences between pair creation in HIC and electron-hole creation in ion bombarded graphene is the presence of the mass gap in HIC, which effectively reduces the pair production rate because it entails higher ion charges for the occurrence of pair production, and the absence of Coulombic bound states before the ion interaction.In other words, in graphene, the transitions occurs between scattering states.Nevertheless, the physical principle underlying electron-hole production in ion bombarded graphene and electron-positron pair production are the same: they are produced via APC.To verify if this phenomenon can be measured experimentally, we study the spectrum of the system and perform a numerical study where the electron-hole pair density is evaluated.

III. GRAPHENE MODEL AND OBSERVABLES
In this article, we consider a regime where the charge transfer from graphene to the ion is minimal and where the integrity of the graphene sample is preserved.According to some experiments [40,54] and time-dependent density functional calculations [34,35], this occurs for ion energy of 1 keV to 2 MeV and higher, and low charge number Z ∼ 1−2.In these conditions, the number of defects and transferred electrons induced by the interaction of the projectile with graphene are negligible.
Provided that these conditions are fulfilled, we choose a theoretical description of charge carriers in terms of a massless Dirac equation where the speed of light c is replaced by the Fermi velocity v F ≈ 0.003c [3,5].This description is accurate as long as the energy of the charge carriers is E 2 eV, where corrections due to the tight-binding model can be neglected [3].Throughout this article, we are assuming that these conditions are fulfilled.It will be verified a posteriori that most electrons are generated with an energy lower than 2 eV.In addition, the interactions between electrons and the environment (through phonon-electron and electronelectron interactions) should be negligible.This can be controlled to a certain extent by using a high dielectric constant substrate and a very low temperature setting, which ensures a small electron-electron coupling constant α G ∼ e 2 /4π v F 1, where is the dielectric constant of the substrate.This condition also allows us to neglect recombination of carriers into photons, a process ∝ α G which should not modify the electron distribution significantly for the electronic densities reached in the system.Finally, the dynamics should occur on a time scale smaller than the thermalization time, on the order of a few tens of femtoseconds [25].Otherwise, the energy distribution of the charge carriers will be simply given by a thermal Fermi-Dirac-like distributions where the features of the dynamics have been washed out.
Considering an interaction of graphene with a charged particle, there is an azimuthal symmetry around the collision point.Then, it is possible and convenient to express the Dirac equation in polar coordinates: where ψ(t, r) is the two-component wave function, r is the radial distance, V is an angle-independent scalar potential, µ 1,2 = j z ∓ 1/2 and j z is the z-angular momentum quantum number taking half-integer values ( As we are interested into the collision of an ion with a graphene sample, the scalar potential represents the field of the moving ion with velocity v and charge Z.It is chosen as where α ≈ 1/137 is the fine structure constant and R(t) := (vt) 2 + r 2 is the distance from the charge.The singular Coulomb potential is regularized by a spherical constant charge distribution for R(t) ≤ R reg , where R reg is the charge radius.The ion radius will be approximately the size of the graphene lattice constant R reg ∼ a ≈ 0.246 nm, similar to Ref. [53].
The first observable considered in this article is the electron distribution on single-particle states (labeled by the energy E k and the angular momentum j z ), defined as where the subscripts in (out) means that the mathematical object is evaluated at time t → +∞ (t → −∞) and where â(out) † jz,k , â(out) jz,k are creation/annihilation operators in the electron-hole representation for the single particle free electron states u kjz (r) (states in the conduction band).We also introduced N v = 2, the number of Dirac valleys, and N spin = 2, the number of physical electron spin, to account for degeneracies.
To take into account the finite nature of the sample, the electronic states are defined on a compact support r ∈ [0, r max ] with a boxed boundary condition, where one of the spinor components is set to zero at r = r max , that is Physically, this corresponds to the zigzag boundary condition, which has been used numerous times to model circular graphene quantum dots [55,56].With these boundary conditions, the free states in both conduction and valence bands are discrete (k ∈ N) because electrons are confined.Other boundary conditions could also be considered like the armchair or the MIT bag model, but this is outside the scope of this article, which focuses on the dynamics of charge carriers.
Using the Furry picture, the electron distribution generated by the strong external field can be calculated from [45]: where t f , t i are the final and initial time, respectively (outside of this time interval, the external field is turned off).We also defined the inner product where the time-dependent wave function ψ kjz (r, t) is a solution to the Dirac equation ( 1) with an initial condition given by a negative eigenstate: ψ kjz (r, t i ) = v kjz (r) at initial time t i , where the field is turned on.Also, we are assuming a null chemical potential (µ = 0) to simulate QED-like process, but this could be relaxed if one is interested in other graphene initial conditions.The second observable considered is the spatial electron density ρ(r), evaluated from It can be verified that one recovers the electron distribution in Eq. ( 5) by integrating the density on the domain and by using the orthogonality of the wave function u kjz (r).

IV. NUMERICAL METHOD
To evaluate the observables, we need to determine the wave function ψ kjz (r, t).As long as αZ/v 1, i.e. for a ratio of the ion charge and velocity which is small enough, a perturbative treatment for finding the wave function is not accurate: perturbation theory does not hold.Therefore, one should solve the Dirac equation Eq. ( 1) nonpertubatively.For this purpose, we employ a high order Galerkin numerical scheme.The rationale for this choice is threefold: first, it can deal with the polar coordinate singularity (in 1/r) in a straightforward way, second, it is very efficient for the evaluation of time-independent states u kjz (r),v kjz (r) required in the calculation and finally, the order of convergence of the spatial discretization is high.
The numerical calculation proceeds in three distinct steps: 1) The time-independent free solutions u kjz (r) and v kjz (r) are determined, 2) The negative energy states are evolved according to the Dirac equation with the ion field and 3) The time-dependent wave function is projected onto positive energy state, using Eq. ( 6).Once the function U jz,k,k is determined, we have access to both observables.
A. Time-independent scheme: Rayleigh-Ritz method The time-independent Dirac equation is given by where E k is the k'th eigenenergy, φ kjz (r) is the k'th eigenstate and V t (r) = V (t, r) is the potential evaluated at some specific time.Two numerical schemes are introduced to solve this eigenvalue problem: one for the free case, when V t (r) = 0 and one for the interacting case, when V t (r) = 0.The former is used to initialize the time-dependent solver while the latter is used to study the spectral characteristics of the system.
1.The free case: Vt(r) = 0 In the free case with V t (r) = 0, the analytical solution of Eq. ( 9) is well-known (see Appendix A).In principle, it is possible to interpolate the solution with the B-spline basis set expansion given below by solving a linear system.However, it is more convenient, for the same computational complexity, to actually solve the eigenvalue problem.
When there is no potential, for j z > 0 (the case j z < 0 is discussed in Appendix B), Eq. ( 9) can be written for the first spinor component as where we focus on the positive energy states by setting φ kjz (r) = u kjz (r).The negative energy solution can be calculated from the positive one as u 1,kjz (r) = v 1,kjz (r) and u 2,kjz (r) = −v 2,kjz (r).The form given in ( 10)-( 11) is particularly convenient because the spectrum of Eq. ( 10) is bounded from below, allowing for a variational Rayleigh-Ritz (RR) method.The latter is equivalent to minimizing the following energy functional obtained from Eq. ( 10) (here, we dropped the index kj z for convenience): where the eigenvalue is λ = . This can be approximated by expanding the solution on a polynomial basis over the domain [0, r max ] as where are the basis functions and the prefactor r |µ1| was added to facilitate the implementation of the boundary condition at r = 0. Indeed, the behavior of the solution is given by u 1 (r) = r |µ1| F (r 2 ) [57], where the function can be Taylor expanded as In this work, a B-spline polynomial basis set {b i (r) have compact support, allowing for numerical sparse matrix linear algebra packages [58].Also, both the Dirac and the Schrodinger equation have been solved with high accuracy using these basis sets [42,[59][60][61][62].The B-spline set and the basis function set are related by The first basis function B 1,1 (r) deserves a special treatment to take the boundary condition at r = 0 into account.This boundary condition can actually be implemented by a suitable combination of B-spline functions, as B-splines are determined by the polynomial order p and the knot vector (see [59] for more details).In this work, an equidistant knot vector is used where each knot is separated by the element size h.In addition, the knot multiplicity is chosen in agreement with the implementation of the boundary conditions.In particular, the multiplicity is p at the endpoint r min = 0, allowing for the canceling of the derivative of the first basis function.On the other hand, the multiplicity is p−1 at r max insuring a boxed boundary condition where u 1 (r max ) = 0. Finally, the multiplicity is 1 at the interior points, guaranteeing the continuity of the solution.This yields a knot vector given by [r 1 , • • • , r 2p+n b −3 ], with knot points ordered as where Here, n b is the number of breakpoints, related to the number of B-spline as An approximation of the eigenenergy and eigenvector is found by minimizing the functional E over the coefficients of the basis function expansion.This is performed as usual in the RR method by reporting the basis set Eq. ( 13) into Eq.( 12).The latter becomes a generalized eigenvalue problem of the form where a is a vector with entries a = (a 1,1 , a 1,2 , • • • , a 1,N ) T while A and D are Hermitian matrices.The entries in these matrices essentially consists in integrals of basis functions (and their first derivative) over the simulation domain.
They are given explicitly in Appendix C 1.Because we are using B-splines for the basis functions, which have compact support, the matrices A and D are sparse.The integrals over basis functions are performed numerically using the Gauss-Legendre (GL) quadrature, while the generalized eigenvalue problem is solved using LAPACK [63].The solution of the eigenvalue problem then yields the first spinor component of the positive and negative energy states (u 1 and v 1 , respectively).The other component u 2 is then evaluated using Eq.(11).Therefore, the second spinor component is given by where a 2,i = i v F E a 1,i , which actually defines the basis function expansion for the second spinor component.This choice of basis expansion for u 2 guarantees that the relation between the two spinor components Eq. ( 11) is fulfilled everywhere on the domain.However, it implies that both components are expressed using different basis functions, similar to techniques for the Dirac equation based on kinetically balanced basis expansion [64].However, the latter is not required in the free and massless case because the spinor components can be decoupled.
It was verified in Appendix D 1 that this numerical scheme reproduces the eigenenergies and eigenstates of the free Dirac operator.The latter can be solved analytically, as shown in Appendix A. In both cases, one can note the presence of an energy gap between positive and negative energy states, given by ∆ = 2v F j 0,1 /r max .This effective gap is induced by the finite size of the graphene sample and vanishes in the limit r max → ∞.
2. The interacting case: V (t, r) = 0 In the presence of a potential, the two spinor components cannot be decoupled as in the free case.For this reason, a slightly different approach is used here, based on a direct discretization of Eq. ( 9).However, the spectrum of this equation is not bounded from below, in contrast to Eq. (10).In this case, the naive utilization of the Rayleigh-Ritz method can lead to spectral pollution, i.e. the appearance of spurious eigenstates [65,66].This problem can be mitigated by using balanced basis functions [64,65] or different variational principles [65,67].In this work, we are using kinetically balanced basis functions [68] using a similar strategy as given in Ref. [69].
From Eq. ( 9), we can obtain the following RR functional: In the kinetically basis set approach, the basis expansion is given by Reporting this basis expansion in Eq. ( 17) yields the generalized eigenvalue problem where S, C, P are matrices similar to the ones in the timeindependent case (their explicit expression is given in Appendix C 2).The vector is also different because it now contains both spinor components, ordered as T .This marks an important difference with the free case, where the solution of the eigenvalue problem yields only the first spinor component while the second component can be calculated from Eq. ( 16).In the interacting case, the solution to Eq. ( 20) gives both spinor components.The free case is slightly more efficient because the number of rows and columns of the matrices in the eigenvalue problem is half of the interacting case.

B. Time-dependent scheme: Galerkin method
Adapting the time-independent solver to the timedependent case is relatively straightforward by using a Galerkin method [44].The basis functions expansion has the same form as Eqs.( 18) and ( 19), but the coefficients become time-dependent as (for j z > 0) The time-dependent basis expansion is then substituted in the Dirac equation (1).As usual, the Galerkin method is then obtained by projecting this equation on the set of all basis functions {(r |µ1| B 1,i , 0), (0, , which allows us to obtain an equation of the form where S, C, P (t) are matrices similar to the ones in the time-independent case (their explicit expression is given in Appendix C 2).The vector contains both spinor components, ordered as a(t) = (a 1,1 (t), a 2,1 (t), • • • , a 1,N (t), a 2,N (t)) T .The initial condition is fixed to the solution obtained from the timeindependent solver by setting a a,i (t i ) = a a,i , for all i ∈ [1, N ] and a = 1, 2.
Eq. ( 23) is a large system of ordinary differential equation which yields the time dependence of the basis coefficients.We solve this system using the implicit Crank-Nicolson method [44] whereby each time iteration ∆t is obtained by solving the following linear system of equations: where we defined a n = a(t n ) and P n+ 1 2 = P (t n + ∆t/2), along with t n = t i + n∆t.The linear system is solved in parallel by using a generalized minimal residual Krylov method implemented in the Petsc high-performance library [70].This numerical scheme has a second order convergence in time.
It was verified that the numerical scheme reproduces the time evolution of a free state and converge towards the exact solution in Appendix D 2. This time evolution can be computed analytically and is given in Appendix A.

V. NUMERICAL RESULTS AND DISCUSSION
This section is devoted to the numerical results obtained from the numerical schemes described in the previous section for the field induced transitions in ion bombarded graphene quantum dots.The first set of numerical results is an analysis of the spectrum of the Coulomb potential, to demonstrate that the electron-hole creation proceeds via the APC mechanism for the parameters considered subsequently.Then, the actual dynamics of charge carriers subjected to the field of the passing ion is evaluated.

A. Spectral characteristics of the Coulomb potential and adiabatic evolution
To determine conditions for which electron-hole pair production occurs via the APC mechanism, we now perform a detailed numerical analysis of the Coulomb potential spectrum.As shown below, the eigenstates of the system provide significant details on the physical processes at play.These results will be important in the next section, for the interpretation of the more complex dynamical results.
The eigenvalues of the Coulomb potential are obtained numerically by solving the time-independent Dirac equation (9) using the numerical method described in Section IV A 2. Two main parameters of the model are varied: the ion charge Z and the graphene-charge distance z.The other simulation parameters are set to values given in Table I, the same values used for evaluating the charge carrier dynamics.For the case where Z is varied, we fix the charge on the graphene sample by setting z = vt = 0.The numerical results are displayed in Fig. 3.
As the ion charge increases, the positive energy states, states of the conduction band with energies E k > 0 for Z = 0, are shifted towards lower energies.Remarkably, there is an avoided crossing between the positive energy state with the lowest energy (PSLE) and the negative energy state with the highest energy (NSHE) when the ion charge reaches a critical value Z crit ≈ 0.36 (this is surrounded by a red diamond in Fig. 3-(a)), indicating that the Coulomb potential becomes supercritical and that electron-hole pairs can be created.Indeed, around this exceptional point, states involved in the avoided crossing become resonances (they actually acquire an imaginary part and cross in the energy imaginary plane [71][72][73]), suggesting that transitions between negative and positive energy states can occur via a non-perturbative tunneling process, where the transition probability in the adiabatic regime is given by the Landau-Zener formula [74,75].These transitions between negative and positive energy states result physically in the creation of electron-hole pairs, as will be shown in the next section.As the ion charge is increased further, other avoided crossings occur between different energy states.These crossings are highlighted by red circles in Fig. 3 for the PSLE and by blue and green diamonds and circles for the first and second states above PSLE, respectively.Each of these crossings corresponds to different potential tunnelling and pair creation pathways.
These first numerical results demonstrate that the potential can become supercritical for some (low) ion charge value, fulfilling one of the criteria of APC.However, in ion bombarded graphene, the charge is not sitting on the graphene sample.To verify if avoided crossings also occur as the ion moves away from the graphene sample, the spectrum is evaluated as a function of the ion-graphene distance z = vt, for an ion charge Z = 1.The results are displayed in Fig. 3-(b), using the same simulation parameters as when the ion charge is varied.These results show again avoided crossings between positive and negative energy states, highlighted by diamonds and circles in the figure.The first avoided crossing of the PSLE occurs at a critical time t crit ≈ −0.11 µm/v, where v is the ion velocity, implying that in the time interval [−t crit , t crit ], the potential becomes supercritical, in agreement with the APC mechanism.Other states also participate and have avoided crossings with negative energy states, allowing for other tunneling pathways.These numerical results demonstrate that the first two criteria of APC given in Sec.II are fulfilled by the physical system under consideration: the potential can localize the electron by forming bound states and it can become supercritical in some time interval.The last criterion for APC is an adiabatic time evolution of the system along some particular times where non-adiabatic transitions occur.From adiabatic perturbation theory, a condition for adiabaticity can be found.It is expressed in terms of the adiabatic factor as [76] for a transition between two eigenstates k and k .This condition is fulfilled, when either the time evolution of the potential (and eigenstates) is slow or when the energy gap between the eigenstates is large.Both conspire to make the system adiabatic, whereby the transition rate between eigenstates k and k is negligible.As an illustration, the adiabatic parameter is evaluated numerically and shown in Fig. 4 for the most important transition of our system, i.e. between PSLE and NSHE (labeled as 0 + and 0 − , respectively).This figure shows clearly that the adiabatic parameter is F 0 + ,0 − (t) > 1 at avoided crossings (for z ≈ −0.01 µm) while it obeys F 0 + ,0 − (t) 1 far from avoided crossings.This naturally leads to the following picture for the dynamics of interband field-induced transitions.When the ion is far from the graphene sample, the dynamics of charge carriers is adiabatic, no transition takes place and the eigenstates accumulate a phase.Then at some critical time (or distance), there is a nonadiabatic transition where the electron-hole pair creation process proceeds via tunnelling, as anticipated from the APC mechanism.The resulting electron distribution will be evaluated in the next section.We note that this process is actually very similar to field induced electronic transitions in molecules [77].
As expected, the adiabatic parameter reaches much higher values as the ion velocity is increased.However, the ion will spend less time around avoided crossings, resulting in a lower tunneling rate in the non-adiabatic transition.Therefore, the picture based on a sequence of adiabatic evolution followed by non-adiabatic transitions may start to fail at high enough v.

B. Field induced interband transitions
The results presented in the last section showed that in some regime, electron-hole pair creation proceeds via  APC.In this section, we are interested in the actual dynamics of charge carriers and distribution when they are subjected to the Coulomb potential of a passing ion.In short, we consider the dynamical production of electronhole pair in ion bombarded graphene quantum dots, the configuration of Fig. 1.The main observables are the electron distribution and the electron density, given in Eqs. ( 5) and ( 8), respectively.They are calculated using numerical solutions of the Dirac equation using the Galerkin method described in Section IV B. The value of important simulation parameters are summarized in Table I.It was tested empirically that converged results are obtained with these parameters by varying the time step and the number of basis functions.More details on the convergence of the numerical scheme can be found in Appendix D. Only the first half of all calculated eigenstates are evolved in time because the error on higher energy states is too large and yields inaccurate results (see Fig. 8 and the discussion in Appendix D 1).
For every numerical results presented in this section, a sum over angular momentum j z is performed.It was verified that for |j z | > 1/2, all the contribution can be neglected because it is at least an order of magnitude below the ones for The numerical results for the electron distribution and Electron distribution Electron density (cm density are displayed in Fig. 5 for an ion charge Z = 1 and for different velocities.The numerical results for the electron distribution in Fig. 5-(a) demonstrate that most electrons are created at low energies, lower than 2 eV where the Dirac model is still a valid approximation.It was demonstrated in numerical studies of HIC that the lower energy electron peak is one of the qualitative features hinting at a supercritical potential [43].
In addition, it was shown in the last section that most avoided crossings occur between low energy states, so the main contribution of these transitions should be found at low energies.For these reasons, the low energy electron peak suggests that electron-hole pairs are produced by the APC mechanism.Another feature of the electron distribution is that more electrons are generated at lower velocities.As argued in the last section, the adiabaticity of the system becomes dubious at higher velocities.Moreover, the faster ions stay for a shorter time in the vicinity of avoided crossings, resulting in a lower number of generated pairs through tunneling.From these considerations, one concludes that lower velocities are better suited for the investigation of APC.
The spatial density at a specific final time t = 65.8 fs is displayed in Fig. 5-(b).Qualitatively, the charge distribution is a ring centered around the ion impact point.The density peak evolves with time to larger radial distance r (not shown here for simplicity).As it evolves, the density is reduced by a simple geometrical effect whereby the ring radius becomes larger and the density covers a larger area.In short, these results demonstrate that electrons are created in the neighborhood of the impact point, when the ion is close to the graphene sample.In the first few instants, the excitation is (quasi-)bound to the ion and localized around r = 0.When the ion gets further from graphene, the Coulomb force is no longer strong enough to localize the charge carrier, so the latter is released and propagates freely in the radial direction, resulting in a radial current.The ion velocity does not modify this picture significantly, except for the fact that higher densities are reached at lower velocity, as expected from the results for the electron distribution.
The numerical results for the electron distribution and electron spatial density as a function of the ion charge are displayed in Fig. 6, for v = 0.001.The results for different charges are qualitatively similar: there is a low energy peak in the electron distribution and the spatial distribution is a ring propagating radially.The main quantitative difference is that more energy states are excited for higher Z, resulting in a higher number of charge carriers.The total number of charge carriers, obtained by summing all energy and angular momentum contributions as n = kjz n kjz is given in Table II and shows a linear increase of the production rate.This can be attributed to the fact that for higher ion charges, more energy states become supercritical, resulting in additional excited tunneling pathways.

VI. CONCLUSION
In this work, the local field interband transitions induced by an ion passing through a graphene quantum dot have been studied theoretically.It is argued by looking at the spectrum and adiabatic conditions that the main process responsible for these transitions is APC, the same mechanism proposed for electron-pair production in low energy HIC.As a consequence, ion-bombarded graphene could serve as a testbed for QED physics in a more forgiv- ing setting.Though APC eludes an experimental verification in HIC, it could be detected experimentally using our configuration.One possibility is to measure the radial induced current I. From the results of Section V B, we can roughly estimate that I = Q/∆t ≈ (2.0 mA) • F , where F is the fraction of the circle where the current is measured and where we assumed that the number of charges generated is n ≈ 100 while the time is estimated as ∆t ≈ 10 fs.According to this order of magnitude estimate, the induced current is in the several microamperes to low milliampere range, which is realistic for an experimental detection.A more elaborated and convincing approach would be to probe the electron distribution (Fig. 5-(a) and 6-(a)) using photoemission electron spectroscopy.However, this would be a challenging experiment, requiring a fine tuning of the ion spatial position and timing.More analysis are required to determine the feasibility of such an experiment.The theoretical investigations have been accomplished by introducing a Galerkin-like numerical scheme to solve the Dirac equation and by performing an extensive num-ber of numerical calculations.The convergence of the numerical scheme has been demonstrated empirically.In principle, the numerical method could be used for other field configurations, such as the one provided by a nanotip, but it is restricted to azymuthally symmetric external potentials.A 2D Cartesian coordinate version of the numerical method would be possible if this condition is not fulfilled, although in this case, other numerical schemes could also be considered [78,79].
Throughout the article, electron-electron interactions have been neglected.We expect that the electron distribution would not be modified significantly by these interactions, given that the obtained distribution are close to the thermal Fermi-Dirac distribution and the time scales involved.However, the effect on the spatial density could be more important and could induce a diffusion of the wave packet at larger times.This important topic is left for future work, along with the effect of different boundary conditions.Finally, the more complex process of carrier recombination could be investigated quantitatively because it could be used as a probe of electron/hole pair generation.and These equations can be decoupled and we get The second components u 2 , v 2 obeys the Bessel equation.Moreover, using the fact that the Bessel function derivative is given by and that µ 1 = µ 2 − 1, we get the solutions where N u,v are normalization constants to be determined and p := p 2 x + p 2 y = E/v F .Now, boundary conditions at r = 0 and r = r max are taken into account.First, it is required that solutions are not singular at the origin, omitting the contribution from Bessel functions of the second kind in the solution.Second, boxed boundary conditions are imposed by setting the first component to zero at the end of the domain, as u 1 (r max ) = v 1 (r max ) = 0.It is not possible to impose this condition to both spinor component at the same time and be consistent with the Dirac equation.Rather, the boxed boundary condition is set to one of the component (here the first) while the other component is evaluated from Eq. (11).Imposing this boundary condition, the momentum p is quantized and takes the value p k = j µ1,k /r max , for k ∈ N + the principal quantum number and j µ,k the zero of the Bessel function.For each value of k corresponds an energy eigenstate.
The normalization constant can now be evaluated by evaluating the product This can be calculated using [80]: Finally, it is verified that the positive and negative states are orthogonal on the domain by evaluating the product Again, using Eq.(A17), it can be demonstrated that u k |v k = 0.
Appendix B: Numerical scheme for jz < 0 For j z < 0, the time-independent Dirac equation is written in a slightly different form where the spinor components are exchanged: This choice for j z < 0 (and also the form for j z > 0) guarantees that singular terms of the form 1/r do not appear in the basis expansion nor in the energy functional, facilitating the numerical implementation.The energy functional for j z < 0 is given by  are known analytically.The first test consists in verifying the spatial convergence as a function of the element size h and the B-spline order.This is accomplished by evaluating the state with the lowest energy and by computing the numerical error, defined as the L 2 -norm of the difference between the exact and the approximate solution = ψ appr − ψ exact L2 .
The results are shown in Fig. 7 for the numerical error while the order of convergence is given in Table III.The results demonstrate that the numerical scheme converges very rapidly as h is reduced and the B-spline order is increased.As the computational time increases with the spline order because the matrices becomes less sparse, there is a compromise between accuracy and efficiency.In this article, we chose p = 4 which converges rapidly but does not induce an important overhead.
In the second test, the accuracy of the numerical scheme is verified as function of eigenstate energies.All the eigenstates obtained from the numerical method are compared to the analytical solution.This is shown in Fig. 8, for many element lengths.These results demonstrate that the error increases with energy, up to a point where it reaches O(1).This occurs because the period of spatial oscillations of the higher energy state wave functions approaches the grid resolution.According to these results, only the bottom half of all eigenstates are included in cal-culations in this article, allowing for accurate solutions.

Convergence of the time-dependent scheme
The convergence of the time-dependent scheme is analyzed by propagating the positive energy ground state (with a principal quantum number k = 0) from time t i = 0 s to t f = 6.58 × 10 −12 s and by the time step dt.The domain size is r max = 1.97 µm, the number of basis function is 40 and the spline order p = 4.The error is evaluated using again the L 2 -norm of the difference between the approximate solution and the exact solution.The analytical exact solution is given in Eqs.(A2) and (A14).The results are shown in Fig. 9.These results demonstrate that the numerical scheme converge with time, with an order of convergence of 3.96.

FIG. 1 .
FIG. 1. Ion bombarded graphene where an ion of charge Z passes through a graphene sample and generates electron-hole pairs.

FIG. 2 .
FIG.2.Single-particle spectrum of the Dirac operator for pair production via the APC mechanism.The blue color represents the possible free positive and negative energy states separated by an energy gap ∆, while the small black lines are bound states.In red are two specific states.The bottom one is shifted as the ion approaches each other, falls into the Dirac sea and crosses with negative energy states.As this occurs, it becomes unstable, its spectral width increases and transitions are possible.The top one is also shifted by the potential but the shift is not large enough to generate transitions.
(p) i (r)} i=1,••• ,Ns , where N s is the number of B-splines and p is the order of the spline polynomial, is chosen because B-splines b (p)

FIG. 3 .
FIG.3.Eigenstates as a function of (a) ion charge Z and (b) ion-graphene distance z = vt.The states of the conduction band (positive energy states) have energies E k > 0 when Z = 0 while the states of the valance band (negative energy states) have E k < 0 when Z = 0.The red diamond highlights the avoided crossing of the ground state with the first negative energy state while red circles are associated with avoided crossings of other negative energy states.The blue (green) diamonds and circles are highlighting avoided crossings for the first (second) positive excited states with negative energy states.

FIG. 5 .
FIG. 5. Electron distribution (a) and electron spatial density (b) as a function of the ion velocity v.The spatial density is shown at time t = 65.8 fs.

FIG. 6 .
FIG.6.Dependence of the electron distribution and spatial density with the ion charge Z.

the eigenvalue is again λ = E 2 v 2 FB 1 5 FIG. 7 .
FIG.7.Numerical error as a function of element size h and spline order p.

80 FIG. 8 .FIG. 9 .
FIG. 8. Numerical error as a function of eigenstate energies.The number on the abscissa represent the principal quantum number.

TABLE I .
Adiabatic parameters as function of ion graphene distance, for many velocities.Simulation parameters for field induced transitions.

TABLE II .
Total number of electrons as a function of the ion charge for v = 0.001.

TABLE III .
Order of spatial convergence for different B-spline order.