Spatio-temporal dynamics in graphene

Abstract

Temporally and spectrally resolved dynamics of optically excited carriers in graphene has been intensively studied theoretically and experimentally, whereas carrier diffusion in space has attracted much less attention. Understanding the spatio-temporal carrier dynamics is of key importance for optoelectronic applications, where carrier transport phenomena play an important role. In this work, we provide a microscopic access to the time-, momentum-, and space-resolved dynamics of carriers in graphene. We determine the diffusion coefficient to be D ≈ 360 cm² s⁻¹ and reveal the impact of carrier–phonon and carrier–carrier scattering on the diffusion process. In particular, we show that phonon-induced scattering across the Dirac cone gives rise to back-diffusion counteracting the spatial broadening of the carrier distribution.

---

The time- and momentum-resolved carrier dynamics in graphene is meanwhile well understood,^1–7 but there have been only a few studies on spatio-temporal dynamics and diffusion in graphene^8–10 and other low dimensional materials, such as carbon nanotubes¹¹ and transition metal dichalcogenides.^12–15 Kulig et al. studied¹⁵ the exciton diffusion in WS₂ and determined that the diffusion coefficient varies over two orders of magnitude with respect to the pump fluence. In graphene, pump–probe experiments performed at relatively high pump fluences^16,17 demonstrated a diffusion coefficient of D = 250 ± 140 cm² s⁻¹ on a picosecond timescale after optical excitation. The diffusion of photoexcited carriers has been studied theoretically¹⁸ with an effective Boltzmann approach, where many-particle scattering has been only considered with relaxation rates.

However, a full microscopic view on the spatio-temporal dynamics revealing the interplay between diffusion and momentum- and time-dependent scattering processes is still missing.

Exploiting the density matrix formalism^19,20 and the Wigner representation,²¹ we provide microscopic insights into the temporally, spectrally, and spatially resolved dynamics of optically excited carriers in graphene including carrier diffusion, carrier–light, carrier–phonon, and carrier–carrier scattering processes on the same microscopic footing, cf. Fig. 1. In particular, we determine the diffusion coefficient and show that the diffusion process can be tuned with experimentally accessible knobs, such as pump fluence, substrate and temperature. Furthermore, we reveal how carrier–phonon scattering counteracts the diffusion through efficient scattering across the Dirac cone resulting in an efficient back-diffusion, cf. Fig. 1(b).

Fig. 1 Spatio-temporal carrier dynamics in graphene. (a) Optically excited carriers in spatial region II diffuse to regions I and III. (b) Many-particle scattering leads to relaxation in momentum space in each spatial region. Back-scattering via carrier–phonon processes and the resulting back-diffusion counteracts the spatial distribution of carriers.

Theoretical approach

We consider a graphene sheet under local optical excitation (red arrows in Fig. 1). The optically excited carriers relax to lower energies via Coulomb- (orange arrows) and phonon-induced scattering (green arrows). The inhomogeneous optical excitation creates spatial gradients in the carrier density, giving rise to diffusion of carriers (purple arrows). To obtain microscopic access to the spatio-temporal dynamics, we derive a set of coupled equations of motion for the electron occupation probability , the microscopic polarization , and the phonon number . Here, the creation and annihilation operators and a_kv/c with momentum k are used for electrons in the valence or conduction band (v,c), while the hole occupation probability is given by ρ^h_k = 1 − ρ^v_k. The corresponding phonon operators are , b_qj with the phonon mode j and the phonon momentum q.

To introduce spatial effects, we transform the occupation probability into the Wigner formalism.^22,23 Here, we consider fluctuations of the occupation probability and perform the Fourier transformation with respect to the momentum difference q resulting in the Wigner function

with λ = e,h denoting electrons in the conduction band and holes in the valence band. Note that the Wigner function f ^λ_k(r) is a quasi-probability function, i.e. f ^λ_k(r) can be negative. Nevertheless, integration over r or k gives the actual distribution in momentum space or the carrier density in real space, i.e. or with A as the area of the graphene sheet.

The carrier dynamics is determined by a many-particle Hamilton operator H, where we take into account the free carrier and phonon contribution H₀, the carrier–carrier H_c–c and the carrier–phonon H_c–ph interaction accounting for Coulomb- (orange arrows) and phonon-induced scattering (green arrows), and the carrier–light coupling H_c–l (red arrows) that is treated on a semi-classical level. Details on the contributions of the many-particle Hamilton operator including the calculation of the matrix elements can be found in ref. 1 and 5.

Exploiting the Heisenberg equation of motion, we derive the equation of motion for the carrier fluctuation ρ^λ_k,q. Taking into account the free-particle Hamilton operator H₀ leads to with the electronic dispersion ε^λ_k. To determine an equation for the Wigner function we perform a Fourier transformation resulting in . To simplify this integro-differential equation we expand the Wigner function to the first order . By using r'e^iqr‘ = −i∇_qe^iqr‘ and shifting the q-derivative to the electron dispersion via partial integration, the r′-integral depends only on the exponential function resulting in δ_k,q, whereby the zeroth order of the expansion of the Wigner function vanishes. Finally, the equation of motion for the Wigner function for the free Hamilton operator reads

To derive the equations of motion for the Wigner function, the polarization and the phonon number with the full Hamilton operator we make the following assumptions: (i) we consider diffusion processes in the polarization to be small, since the latter quickly decays in momentum space and vanishes directly after the optical excitation.⁵ In contrast, the relaxation of carriers occurs on a picosecond timescale which is comparable to diffusion processes, and therefore the diffusion term can not be neglected in the equation for the Wigner function. (ii) We also neglect the phonon diffusion, since it is expected to be much slower than the electronic diffusion due to the flat phonon dispersion. (iii) We expect scattering processes between different spatial positions to be small compared to the diffusion. Now, using the Heisenberg equation of motion, we derive the full spatio-temporal graphene Bloch equations in second-order Born-Markov approximation

with the abbreviations , and with the initial Bose-distribution for phonons n^j_q,B. The equations describe the time-, momentum- and space-resolved coupled dynamics of electrons/holes, phonons, and the microscopic polarization. The dynamics of electrons in the conduction band and holes in the valence band is symmetric, but has different initial conditions for doped graphene samples. The appearing Rabi frequency is defined as with the free electron mass m₀, the vector potential A(r,t), and the optical matrix element M^vc_k = 〈kv|∇_k|kc〉. Since we study the carrier dynamics close to the Dirac point, renormalization effects can be neglected. Furthermore, we have introduced with the electronic dispersion ε^λ_k and the dephasing rate γ_k(r,t). The time-, momentum- and spatial dependent dephasing γ_k(r,t) and in- and out-scattering rates Γ^in/out_kλ(r,t) include carrier–carrier and carrier–phonon scattering channels. The dynamics of the phonon number n ^j_q(t) is driven by the emission and absorption rates^1,5 Γ^em/abs_qj(r,t). The constant γ_ph is the experimentally determined phonon decay rate.²⁴ More details on the appearing many-particle scattering and dephasing rates can be found in ref. 1 and 5. In this work, we assume that graphene lies on a SiC-substrate and is surrounded by air on the other side. This is taken into account by introducing an averaged dielectric background constant²⁵, where ε_s = 9.66 is the static screening constant of the substrate, while 1 describes the dielectric constant of air. Furthermore, the internal many-particle screening is taken into account by calculating the static limit of the Lindhard equation,^19,26 which screens the Coulomb matrix elements.

The derived set of equations resemble the semiconductor Bloch equations for spatial homogeneous systems (cf.ref. 1 and 5) up to the additional term , which describes the diffusion of carriers in the direction ∇_kε^λ_k ∝ e_k = k/|k|. As a result, carriers with different sign in momentum move in opposite directions generating locally asymmetric carrier distributions in momentum space and resulting in a local current with the Fermi velocity v_F. The sum contains both electrons in the conduction band and holes in the valence band and in a spatially homogeneous system, the mean current vanishes.

Spatio-temporal dynamics

Now, we numerically evaluate the spatio-temporal graphene Bloch equations and investigate the interplay of diffusion and relaxation processes after optical excitation. We excite carriers with an optical pulse with a Gaussian profile both in time and space. We chose typical values for pulse characteristics including a temporal FWHM of 115 fs, a spatial FWHM of 265 nm, an excitation energy of 1 eV and a pump fluence of 1 μJ cm⁻². The temporally and spatially dependent carrier density n(x,t) is shown in Fig. 2(a). The diffusion of carrriers is reflected in the broadening of the carrier density in space. Normalizing the density for each time step, the broadening becomes more visible (Fig. 2(b)), since phonon- and Auger-driven interband processes give rise to a reduction of carriers with increasing time. Fig. 2(c) shows the temporal and spatial distribution of the electronic temperature. The latter has been determined via the average kinetic energy per particle. Assuming a Fermi distribution, one obtains . Note that shortly after the pulse the carrier distribution deviates from an equilibrium Fermi distribution and therefore the temperature is there not well defined.

Fig. 2 Spatio-temporal dynamics. (a) Carrier density n(x,t) is shown as a function of time and position after a local optical excitation centered at x = 0. (b) Carrier density normalized to the maximal density at each time step to highlight the diffusion of carriers. (c) Carrier temperature defined via the kinetic energy per particle. (d) The square of the width w of the spatial distribution is a measure for the diffusion coefficient D. (e) Diffusion coefficient as a function of different pump fluences, substrates and temperatures. The x-axis is normalized to the maximal value of the respective parameter.

To quantify the diffusion and to estimate the diffusion coefficient for graphene, we fit the carrier density with a Gaussian exp(−x²/w²(t)) for every time step. The temporal evolution of the width w(t) is depicted in Fig. 2(d). It is connected to an effective diffusion coefficient D via²⁷w²(t) = w₀²+ 4Dt resulting in D ≈ 360 cm² s⁻¹ for the investigated graphene sample on a SiC substrate. Our results fit well to the experimentally obtained values¹⁷ for the diffusion coefficient of D = 250 ± 140 cm² s⁻¹. In Fig. 2(e) we show the influence of pump fluence, substrate and temperature on the diffusion coefficient. We find that the temperature has the largest impact. The underlying processes will be discussed later.

So far we have investigated the diffusion behaviour by taking into account the full carrier dynamics. Now, to understand the fundamental interplay of many-particle scattering and diffusion processes we study separately the impact of different scattering mechanisms on the diffusion process, cf. Fig. 3. We start with the case without any scattering channels just considering the electron–light interaction. After the optical excitation, carriers with positive/negative momenta diffuse in opposite spatial directions according to the diffusion term in eqn (3). After approximately 100 fs the carrier separation becomes visible, as the intial carrier density distribution splits into two pronounced peaks of the same width but with half of the amplitude, cf. Fig. 3(a). Including the carrier–phonon scattering, we observe a strongly reduced spatial broadening of the carrier density and no splitting appears (Fig. 3(a)). Phonon-induced relaxation processes counteract the diffusion via back-scattering across the Dirac cone and the following back-diffusion (cf.Fig. 1). The impact of carrier–phonon scattering will be further microscopically resolved in the next section. Including only the carrier–carrier scattering, the density diffuses with the same speed as in the case without any scattering channels (cf.Fig. 3(c)). This is a consequence of the symmetry of Coulomb matrix elements, which favor parallel scattering.^28,29 Scattering across the Dirac cone is relatively inefficient and back-scattering is even forbidden. In contrast to the case without scattering, the spatial region between the two peaks contains a non-zero density. This reflects the weak but not vanishing Coulomb scattering processes bringing carriers from one to the other side of the Dirac cone.

Fig. 3 Impact of scattering channels on diffusion. (a) After optical excitation the carriers on Dirac cone branches with different sign of momenta diffuse in opposite directions. (b) Carrier–phonon scattering counteracts the diffusion and the broadening of the distribution due to back-scattering between Dirac cone branches and the following back-diffusion (c). (c) Carrier–carrier scattering does not effectively counteract the diffusion, since parallel scattering is preferred by the Coulomb matrix elements. Thus, one can still clearly observe the spatial separation of carriers with opposite momentum.

Carrier–phonon dynamics

To get a thorough understanding of the microscopic processes governing the spatio-temporal carrier dynamics, we investigate the spectral and spatial behaviour of the Wigner function for different times. We start with the interplay of diffusion and carrier–phonon scattering processes. The optically excited carriers scatter via optical phonons to lower energies and form enhanced carrier occupations separated by the energy of optical phonons (red regions in Fig. 4(a)). Diffusion processes lead to a spatial broadening of the carrier distribution and after approx. 1 ps the carriers have relaxed to lower energies close to the Dirac point (Fig. 4(b)).

Fig. 4 Impact of carrier–phonon scattering. (a), (b) Wigner function is shown in dependence of space and momentum for two fixed times. Time 0 ps corresponds to the maximum of the optical excitation pulse. Note that we plot here the change in the Winger function with respect to the intial time t₀ before the optical excitation, i.e. f_k(x,t) = f_k(x,t) − f_k(x,t₀). (c), (d) Illustration of the impact of diffusion by showing the Wigner function minus the case without diffusion. The overall spatial carrier distribution becomes broader at larger times due to diffusion. Backscattering with acoustic phonons leads to additional peaks at lower energies. Since here the momentum of the carriers is flipped, they diffuse back resulting in multiple sign changes in (d).

To investigate the impact of diffusion in more detail we performed the same calculation twice, but in the second computation we excluded diffusion processes. Illustrating the difference of both calculations, i.e. f_k(x) − f_k(x)^nodiff, we can directly observe the impact of diffusion on carrier–phonon scattering (Fig. 4(c) and (d)). As already discussed in the theory section carriers with positive/negative momentum diffuse in opposite spatial direction. This behaviour is illustrated in Fig. 4(c), where carriers with positive momentum diffuse from x < 0 positions (orange spots) to x > 0 positions (red spots). After 1 ps the carriers have already relaxed to energies close to the Dirac cone and below the optical phonon energy. Consequently, the scattering with acoustic phonons becomes dominant. Due to the flat dispersion of acoustic phonons with respect to the Dirac cones back-scattering across the Dirac cone is preferred, such that carriers with positive momenta are scattered to negative momenta and vice versa (Fig. 1). The inversion of momenta results in a back-diffusion, such that the overall carrier distribution stays bunched in space, cf. Fig. 3(b). The back-diffusion is shown in Fig. 4(d) by the the multiple sign change in the colored regions (red to orange to red).

Carrier–carrier dynamics

Now, we investigate the impact of carrier–carrier scattering on diffusion of optically excited carriers. An important aspect here is that Auger scattering is efficient giving rise to a carrier multiplication^30–34 that increases the overall carrier density (note the scale of the color map in Fig. 5 compared to Fig. 4(a)). This also results in a quick increase of the carrier distribution close to the Dirac cone already during the optical excitation (Fig. 5(a)). Since electrons and holes diffuse in the same direction, the conditions for carrier multiplication are still satisfied after the diffusion. The directional dependence (in momentum space) for intraband carrier–carrier scattering is determined by the Coulomb matrix element that includes a form factor proportial to¹ 1 + e^iϕ with the scattering angle ϕ. This means that parallel scattering (ϕ = 0) is the preferable scattering channel, and that for the back-scattering (ϕ = π) the amplitude of the Coulomb matrix element completely vanishes. As a result, scattering processes across the Dirac cone that change the sign of the carrier momentum (and lead to a back-diffusion) are inefficient. As a result, carriers with positive/negative momenta remain separated with respect to their spatial position – similarly to the case without any scattering (Fig. 3(a)). Fig. 5(b) illustrates that carriers with positive/negative momenta are mainly distributed towards positive/negative spatial positions.

Fig. 5 Impact of carrier–carrier scattering. The same as Fig. 4, but now only including the carrier–carrier instead of carrier–phonon scattering. Coulomb interaction is most efficient for parallel scattering along the Dirac cone in the momentum space. As a result, carrier–carrier scattering does not efficiently counteract the diffusion of carriers in opposite direction as clearly observed in (b).

Tuning the diffusion

Now, we can explain the dependence of the diffusion coefficient on pump fluence, substrate and temperature shown in Fig. 2(e). We find that the diffusion becomes less efficient with the increasing pump fluence (blue curve). Here, more carriers are excited resulting also in an increased number of emitted phonons. Thus, hot-phonon effects become important, i.e. an increasing number of phonons can be reabsorbed in back-scattering processes giving rise to additional channels for back-diffusion. Furthermore, we find that the diffusion coefficient is nearly independent of the substrate (red curve) entering in our calculations through the screening of the Coulomb potential. This is not surprising, since Coulomb-induced scattering processes have been shown to only play a minor role for the diffusion of carriers, cf. Fig. 3. Finally, we observe that the diffusion can be most efficiently tuned by varying the temperatures (orange curve). The lower the temperature, the weaker the carrier–phonon scattering, the less efficient is back-scattering and back-diffusion resulting in a considerably increased diffusion coefficient.

In summary, we provide a microscopic view on the spatio-temporal carrier dynamics in graphene based on the density matrix formalism in Wigner representation. We investigate the interplay of diffusion and many-particle scattering processes after a local optical excitation. In particular, we determine a diffusion coefficient of D ≈ 360 cm² s⁻¹ that agrees well with recent experimental values. Furthermore, we reveal that carrier–phonon scattering across the Dirac cone and the resulting back-diffusion are crucial ingredients to understand the spatial broadening of the carrier distribution. The gained insights are important e.g. for graphene-based photodetectors,^35–39 that are governed by the thermoelectric effect, which relies on spatial temperature gradients.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No. 785219 (Graphene Flagship). Furthermore, we acknowledge support from the Swedish Research Council (VR). The computations were performed on resources at Chalmers Centre for Computational Science and Engineering (C3SE) provided by the Swedish National Infrastructure for Computing (SNIC).

References

1. E. Malic and A. Knorr, Graphene and Carbon Nanotubes: Ultrafast Optics and Relaxation Dynamics, Wiley-VCH, 2013.
2. S. Butscher, F. Milde, M. Hirtschulz, E. Malic and A. Knorr, Appl. Phys. Lett., 2007, 91, 203103.
3. J. M. Dawlaty, S. Shivaraman, M. Chandrashekhar, F. Rana and M. G. Spencer, Appl. Phys. Lett., 2008, 92, 042116.
4. P. Plochocka, P. Kossacki, A. Golnik, T. Kazimierczuk, C. Berger, W. A. de Heer and M. Potemski, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 80, 245415.
5. E. Malic, T. Winzer, E. Bobkin and A. Knorr, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 205406.
6. D. Sun, C. Divin, M. Mihnev, T. Winzer, E. Malic, A. Knorr, J. E. Sipe, C. Berger, W. A. de Heer, P. N. First and T. B. Norris, New J. Phys., 2012, 14, 105012.
7. E. Malic, T. Winzer, F. Wendler, S. Brem, R. Jago, A. Knorr, M. Mittendorff, J. C. König-Otto, T. Plötzing, D. Neumaier, H. Schneider, M. Helm and S. Winnerl, Ann. Phys., 2018, 1700038.
8. L. Huang, G. V. Hartland, L.-Q. Chu, Luxmi, R. M. Feenstra, C. Lian, K. Tahy and H. Xing, Nano Lett., 2010, 10, 1308.
9. R. Rengel, E. Pascual and M. J. Martin, Appl. Phys. Lett., 2014, 104, 233107.
10. Y. Ishida, H. Masuda, H. Sakai, S. Ishiwata and S. Shin, Phys. Rev. B, 2016, 93, 100302.
11. J. H. Grönqvist, M. Hirtschulz, A. Knorr and M. Lindberg, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 035414.
12. T. Kato and T. Kaneko, ACS Nano, 2016, 10, 9687.
13. L. Yuan, T. Wang, T. Zhu, M. Zhou and L. Huang, J. Phys. Chem. Lett., 2017, 8, 3371.
14. R. Rosati, F. Lengers, D. E. Reiter and T. Kuhn, Phys. Rev. B, 2018, 98, 195411.
15. M. Kulig, J. Zipfel, P. Nagler, S. Blanter, C. Schüller, T. Korn, N. Paradiso, M. M. Glazov and A. Chernikov, Phys. Rev. Lett., 2018, 120, 207401.
16. B. A. Ruzicka, S. Wang, L. K. Werake, B. Weintrub, K. P. Loh and H. Zhao, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 195414.
17. B. A. Ruzicka, S. Wang, J. Liu, K.-P. Loh, J. Z. Wu and H. Zhao, Opt. Mater. Express, 2012, 2, 708.
18. F. T. Vasko and V. V. Mitin, Appl. Phys. Lett., 2012, 101, 151115.
19. H. Haug and S. W. Koch, Quantum Theory of the Optical and Electronic Properties of Semiconductors, World Scientific, 2009.
20. M. Kira, F. Jahnke, W. Hoyer and S. Koch, Prog. Quantum Electron., 1999, 23, 189.
21. E. Wigner, Phys. Rev., 1932, 40, 749.
22. O. Hess and T. Kuhn, Phys. Rev. A, 1996, 54, 3347.
23. F. Rossi and T. Kuhn, Rev. Mod. Phys., 2002, 74, 895.
24. K. Kang, D. Abdula, D. G. Cahill and M. Shim, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 165405.
25. L. Patrick and W. J. Choyke, Phys. Rev. B: Solid State, 1970, 2, 2255.
26. M. Kira and S. Koch, Prog. Quantum Electron., 2006, 30, 155.
27. R. Pathria and P. Beale, in Statistical Mechanics, Academic Press, Boston, 3rd edn, 2011, pp. 583–635.
28. E. Malic, T. Winzer and A. Knorr, Appl. Phys. Lett., 2012, 101, 213110.
29. M. Mittendorff, T. Winzer, E. Malic, A. Knorr, C. Berger, W. A. de Heer, H. Schneider, M. Helm and S. Winnerl, Nano Lett., 2014, 14, 1504.
30. T. Winzer, A. Knorr and E. Malic, Nano Lett., 2010, 10, 4839.
31. D. Brida, A. Tomadin, C. Manzoni, Y. J. Kim, A. Lombardo, S. Milana, R. R. Nair, K. S. Novoselov, A. C. Ferrari, G. Cerullo and M. Polini, Nat. Commun., 2013, 4, 1987.
32. T. Plötzing, T. Winzer, E. Malic, D. Neumaier, A. Knorr and H. Kurz, Nano Lett., 2014, 14, 5371.
33. M. Mittendorff, F. Wendler, E. Malic, A. Knorr, M. Orlita, M. Potemski, C. Berger, W. A. de Heer, H. Schneider, M. Helm and S. Winnerl, Nat. Phys., 2015, 11, 75.
34. I. Gierz, M. Mitrano, J. C. Petersen, C. Cacho, I. C. E. Turcu, E. Springate, A. Stöhr, A. Köhler, U. Starke and A. Cavalleri, J. Phys.: Condens. Matter, 2015, 27, 164204.
35. F. H. L. Koppens, T. Mueller, P. Avouris, A. C. Ferrari, M. S. Vitiello and M. Polini, Nat. Nanotechnol., 2014, 9, 780.
36. Z. Sun and H. Chang, ACS Nano, 2014, 8, 4133.
37. M. Buscema, J. O. Island, D. J. Groenendijk, S. I. Blanter, G. A. Steele, H. S. J. van der Zant and A. Castellanos-Gomez, Chem. Soc. Rev., 2015, 44, 3691.
38. S. Schuler, D. Schall, D. Neumaier, L. Dobusch, O. Bethge, B. Schwarz, M. Krall and T. Mueller, Nano Lett., 2016, 16, 7107.
39. J. C. W. Song, M. S. Rudner, C. M. Marcus and L. S. Levitov, Nano Lett., 2011, 11, 4688.

---