CC BY
30
udc 519.688, 004.942
A. D. Panferov, A. V. Makhankov
Simulation of the effect of short optical pulses on
graphene
Abstract. The interaction of high-frequency pulsed electric fields with graphene is currently the subject of intense research. The paper presents the results of testing a software system for modeling such processes using the example of ultrashort laser pulses of the optical range with different polarizations. The authors develop the system on a base of a new theoretical approach based on the quantum kinetic equation. The approach contains a computational model for a new system of ordinary differential equations with non-linearly dependent on time and problem parameters coefficients.
The need to analyze the behavior of solutions of this system of equations in the field of changing several parameters leads to the polynomial computational complexity. The lack of knowledge of the nature of the parametric dependence of solutions requires several iterations of the choice of covering grids. The paper describes the adaptation of this modeling system for use in massively parallel computing systems.
Key words and phrases: numerical simulation, high performance computing, graphene, kinetic equation.
2010 Mathematics Subject Classification: 81T80; 81T40, 82C20
Introduction
Currently, the features of optical radiation interaction with new two-dimensional or pseudo-two-dimensional materials realized in the form of monoatomic or monomolecular layers on substrates or in the free state are intensively studied. Graphene is one of these perspective materials [1]. Due to the extremely high mobility of charge carriers [2,3] and unusual optoelectric properties, [4-6] it can be used in the development of new
Supported by RFBR according to the research project 18-07-00778.
© A. D. Panferov, A. V. Makhankov, 2019 © Saratov State University, 2019
© Program Systems: Theory and Applications (design), 2019
ultra-fast optoelectronic devices [7-9]. Particularly interesting are its nonlinear properties in the infrared region of the spectrum and optical range, leading to the generation of high-frequency harmonics [10-13]. This direction requires complex experimental technologies and using of unique equipment. Theoretical models of these processes are essential both for understanding their physical nature and for assessing the achievable parameters and characteristics.
Even with the simplicity and transparency of the physical principles of the model, the procedure for predicting the observed values on its basis can be a computationally complex resource-intensive task. The present work considers a test problem, which is devoted to demonstrating the capabilities of the developed software system for simulating the behavior of charge carriers in graphene in the presence of external electric fields of various nature using a new kinetic quantum field approach [14]. The universality of this approach makes it possible to consider electric fields with an arbitrary dependence on time, including a variable direction of action. The only limiting condition is the requirement of spatial homogeneity of the field. We test the system on the problem with realistic parameters under the conditions of experiments presented in the works [15,16].
1. Physical model
The basic kinetic equation(KE) which describing electron-hole excitations in electric field for D = 2 + 1 graphene model is [14]:
(1)
t
f(PuP2,t) = 1A (p1,p2,t)J dt'\(pup2,t') [1 - 2 f (pi,p2,0]c os 0(t,t'),
to
where
t
(2) 0(t,t') = |y dt"e(pi,p2, t ''),
t'
and parameter
(3) A ( t) evF [E1P2 - E2Pi] (3) A (pi,p2,t) = -2?-7\-
£2(pi,p2,t)
determined through the components of a two-dimensional (k = 1, 2) quasimomentum Pk = pk — (e/c)Ak(t) and quasiparticles energy
(4) £(pt)= VF VP2 = VF V(Pl)2 + (P2)2,
taking into account the specifics of the dispersion law in the neighbourhood of Dirac points. In these definitions: ft—reduced Planck constant, e— elementary charge, c—velocity of light, vF —Fermi velocity, pk—components of a quasi-particle momentum in the absence of an electric field, Ek— strength and Ak —vector potential of electric field.
For a numerical study of the behavior of KE (1) solutions, it is convenient to convert it into the form of a system of three ordinary differential equation [17]
1 2e 2e
(5) f = 2 Ли u = Л (1 — 2f) — "ft v, V = "ft U
with appropriate initial conditions f (t0) = u(t0) = v(t0) =0 for a physical system, which occupies a vacuum state before the field switched on. The primary object to describe the characteristics of the system is the distribution function f (p1,p2,t). Functions u(p1,p2,t) and v(p1,p2,t) are auxiliary and represent the features of the evolution of polarization effects and energy of the system.
To describe the electric field of short optical pulses, we define it in the form:
(6) E1(t) = E0 cos(wi) exp(—t2/2r2) E2(t) = E0 cos(wt + ф) exp(—t2/2r2).
Here ш - cyclic radiation frequency. The characteristic of the pulse duration т can be determined through an auxiliary parameter a (usually interpreted as the number of fundamental frequency waves per pulse) using relation т = a/ш. The phase shift ф between the components of the field vector determines its polarization. In the case ф = 0 polarization is linear, ф = 0.5n corresponds to circular polarization.
We have selected the parameter values in accord with the experimental conditions [16]: ш = 2n • 5.0x1014Hz, that corresponds to the wavelength Л = 600 nm (optical range). The amplitude of electric field is Eo — 2.528x104 V/m. The selected field amplitude value corresponds to focusing into a spot with a diameter of 1.5 pm at a pulse power of 50 pW. We do not fix the ф parameter. Changing it, we can simulate the process of changing the polarization of a laser pulse.
2. Software implementation
Implementation of the simulation process should provide complete information about the distribution function of the carriers at the time of the termination of the electric field. Strictly, it is f (pi;p2,t ^ to), and the integration must be performed for —to < t < to. But, since the field (6) is very quickly suppressed by the exponent exp(— t2/2r2), real time interval t;n < t < tout for describing distribution function evolution through the system (5) is finite and amounts to several t. Its specific value is chosen based on the numerical experiments from the condition of integration
results stability.
Since in real experiments the integral characteristics of the distribution function are measured, it is necessary to obtain sufficiently detailed information about f (pi,p2 : t = tout) = fout(pi,p2) to adequately reproduce the values of these characteristics in the subsequent stages of modeling the system behavior. In general, a priori we know nothing about its behavior in the momentum space. According to the conditions of the problem, only the first Brillouin zone is considered. But whether fout(pi,p2) is localized in a relatively small area of this zone, distributed over several compact areas, or relatively uniformly covers the entire zone can be found only by the results of numerical experiments. It follows that the mandatory step of solving the problem will be the determination of range arguments values pi and p2. Next problem is a sampling step selection. If fout(pi,p2) turns out to be smooth, we can use sparse grid. With the complex behavior of the distribution function, grid with a vast number of points may be required. These problems will be demonstrated by specific examples below.
An important simplifying circumstance in this case is the refusal to take into account the dissipative processes and the reverse influence of the generated quasiparticles on the active field. It allows us to solve a system of equations (5) for each point of the momentum space independently and work with the grid points in parallel. In software, this is implemented using MPI procedures. MPI was selected because it is necessary to have the opportunity to simulate on massive parallel systems without limit of scaling
level.
We solve the system of ODE using the GSL library [18]. The choice of GSL is determined by the presence of different levels objects for the organization of procedure for solving the Cauchy problem of the ODE system and the support of various methods for the numerical solution of such problem. There are implementations of several modifications of the explicit and three modifications of the implicit Runge-Kutta methods,
special features are provided for the solution of stiff systems. Effectiveness of different methods was studied in detail [19].
3. Results of numerical simulation
Figure 1 shows the steps of determining the shape of the residual distribution function formed as a result of the action of a field pulse (6). Frequency and amplitude were defined in Section 1. For two additional parameters, values are selected: ^ = 0 and a = 4. As noted above, the choice of ^ = 0 specifies the common-mode behavior x and y components of field (6). In this case, their sum is linearly polarized with the polarization plane on the diagonal x = y. The choice of a = 4 minimizes the field duration and the time required to the simulation.
Figure 1a presents a quick overview of the Brillouin zone. It was building by 169 points (grid 13x13), and it gives only an approximate idea of the solution localization area. Figure 1b was building taking into account the obtained information about the region of the distribution function localization in the momentum space. The computation were carried out in 484 points (grid 22x22). The account the refinement of the localization area has improved the resolution in the momentum space by 5 times. The figure gives a rough representation of the distribution form and does not guarantee the correct reproduction of integral characteristics. On the next iteration, the grid step reduces in 5 times by every axes. The result presents on Figure 1c. In this case equations system (5) solves for 10404 various samples of parameters (grid 102x102), that provides high potential for parallel. Resulting data allow us to be sure both in the smoothness of the studied solution and in the exact reproduction of its integral parameters. Also, they demonstrate the reflection of the spatial symmetry of the acting electric field in the momentum space.
Choosing a more realistic value for the parameter a = 16 significantly complicates the task. First of all, the duration of field action increases proportionally. But other factors affect the computational complexity of the problem. Figure 2a presents the results of computation on a grid similar to that used for the bottom panel Figure 1.
We see the distribution function changes. Its maximum value increases approximately by an one order of magnitude, but the "thickness" of the localization area decreases significantly. As the result, a grid 102x102 does not provide exact reproducing its form. We have to increase the number of grid points on each axis of the momentum space. Figure 2 b shows a result
(a) an overview of the entire Brillouin zone
Figure 1. Three stages define the shape of distribution function for residual carriers
(a) -1.0 < pi < 1.0, -1.0 < p2 < 1.0
Figure 2. Building of the distribution function for the field momentum with the parameter value a = 16, grid 102x102 and different areas
for area 0.0 < pi < 0.5, —0.5 < p2 < 0.0 with grid step reduced in 4 times. So if you do not take additional simplifying actions, the complexity of the task is approximately proportional to the cube of parameter a. This should be taken into account, since in realistic situations this parameter has rather large values.
In the model under consideration, there is another free parameter the manipulations of which allow varying the polarization of the incident
(a) 0 = 0.3n
(b) 0 = 0.5n
Figure 3. The effect of the phase shift 0 on the shape of the distribution function
radiation. Model behavior when changing the degree of polarization demonstrates Figure 3. Both figures use the same set of parameters as the first series of figures (Figure 1) and a grid 102x102. In conjuction Figure 1c, Figure 3a and Figure 3b reflects the phase shifts 0 = 0, 0 = 0.3n and 0 = 0.5n respectively. This series of figures shows the change in the symmetry of the solution. In the case of circular polarization 0 = 0.5n the residual distribution function acquires rotational symmetry with respect to the point p1 = p2 = 0. This potentially opens the possibility to reduce the dimension of the problem from 2D to 1D with a corresponding decrease in
Table 1. Computation time to solve the task for a single point
Integration method (in designations of GSL) Time (ms)
Explicit embedded Runge-Kutta-Fehlberg (4, 5) 304.7
Explicit embedded Runge-Kutta Cash-Karp (4, 5) 194.0
Explicit embedded Runge-Kutta Prince-Dormand (8, 9) 64,3
A variable-coefficient linear multistep Adams method in 144 5 Nordsieck form
its computational complexity. But in general case of the ^ values there is only a reflective symmetry with respect to the line pi = p2. The values of ^ has no noticeable effect on the computation time.
All the above results are obtained using an explicit method in the solution of the equations system (5), referred to in the GSL documentation as "embedded Runge-Kutta Prince-Dormand (8, 9)". Its excellence has been demonstrated earlier [19] for another field model. Since the field model under consideration (6) has significant differences, a control comparison of the solution time was computed per grid point for several methods that showed the best results in [19]. Tests were performed for the value a = 256 in the area maximum values of the distribution function and in the area of background values on nodes with processors Intel® Xeon® E5405 with clock frequency 2.0 GHz. The computation times for points in the area of maximum values are given. Outside this area, the computation times do not exceed these values. From the given data it follows that "embedded Runge-Kutta Prince-Dormand (8, 9)" the method still copes with the task more than twice as fast as the nearest competitor.
Conclusion
The article presents the results of testing the system for modeling the interaction of graphene with external electric fields on the example of a realistic physical problem. The result of the action on a graphene of a short optical pulse is considered. The possibility of reproducing in detail the result of such interaction in the presence of sufficient computing resources is demonstrated. A feature of the developed system is a new physical model based on the quantum kinetic equation [14]. This allows modeling without restrictions on the nature of the dependence of the active field on time, its frequency, amplitude, and direction.
The authors thank S.A. Levenec, S.A.Smolyansky and S.O.Pirogov for useful discussion of the work.
References
[1] M. Baudisch, A. Marini, J. D. Cox, T. Zhu, F. Silva, S. Teichmann, M. Massicotte, F. Koppens, L. S. Levitov, F. J. Garcia de Abajo, J. Biegert. "Ultrafast nonlinear optical response of Dirae fermions in graphene", Nature Communications, 9 (2018), 1018. I ' 47
[2] L. Wang, I. Meric, P. Y. Huang, Q. Gao, Y. Gao, H. Tran, T. Taniguchi, K. Watanabe, L. M. Campos, D. A. Muller, J. Guo, P. Kim, J. Hone, K.L. Shepard, G. R. Dean. "One-dimensional electrical contact to a two-dimensional material", Science, 342 (2013), pp. 614-617. ' 47
[3] L. Banszerus, M. Schmitz, S. Engels, M. Goldsche, K. Watanabe, T. Taniguchi, B. Beschoten, C. Stampfer. "Ballistic transport exceeding 28 fim in CVD grown graphene", Nano Lett., 16 (2016), pp. 1387-1391. 47
[4] J. Wang, Y. Hernandez, M. Lotya, J. N. Coleman, W. J. Blau. "Broadband nonlinear optical response of graphene dispersions", Advanced Materials, 21 (2009), pp. 2430-2435. 47
[5] Z. Fei, G.O. Andreev, W. Bao, L. M. Zhang, A. S. McLeod, G. Wang, M.K. Stewart, Z. Zhao, G. Dominguez, M. Thiemens, M. M. Fogler, M.J. Tauber, A. H. Castro-Neto, C.N. Lau, F. Keilmann, D.N. Basov. "Infrared nanoscopy of Dirac plasmons at the graphene-Si02 interface", Nano Lett., 11 (2011), pp. 4701-4705. 47
[6] A. E. Nikolaenko, N. Papasimakis, E. Atmatzakis, Z. Luo, Z.X. Shen, F. Angelis, S.A. Boden, E. Fabrizio, N.I. Zheludev. "Nonlinear graphene metamaterial", Appl. Phys. Lett., 100 (2012), 181109. i ' 47
[7] M. Garg, M. Zhan, T. T. Luu, H. Lakhotia, T. Klostermann, A. Guggenmos,
E. Goulielmakis. "Multi-petahertz electronic metrology", Nature, 538 (2016), pp. 359-363. 48
[8] A. Sommer, E. M. Bothschafter, S.A. Sato, C. Jakubeit, T. Latka, O. Razskazovskaya, H. Fattahi, M. Jobst, W. Schweinberger, V. Shirvanyan, V. S. Yakovlev, R. Kienberger, K. Yabana, N. Karpowicz, M. Schultze,
F. Krausz. "Attosecond nonlinear polarization and light-matter energy transfer in solids", Nature, 534 (2016), pp. 86-90. I 148
[9] T. Higuchi, C. Heide, K. Ullmann, H. B. Weber, P. Hommelhoff. "Light-field-driven currents in graphene", Nature, 550 (2017), pp. 224-228.
^48
[10] N. Kumar, J. Kumar, C. Gerstenkorn, R. Wang, H. Chiu, A.L. Smirl, H. Zhao. "Third harmonic generation in graphene and few-layer graphite films", Phys. Rev. B, 87 (2013), 121406. I 148
[11] S. Hong, J. I. Dadap, N. Petrone, P. Yeh, J. Hone, R. M. Osgood, Phys. Rev. X, 3 (2013), 021014. i 48
[12] G.X. Ni, L. Wang, M. D. Goldflam, M. Wagner, Z. Fei, A. S. McLeod, M. K. Liu, F. Keilmann, B. Ozyilmaz, A. H. Castro Neto, J. Hone, M. M. Fogler, D. N. Basov. "Ultrafast optical switching of infrared plasmon
polaritons in high-mobility graphene", Nature Photonics, 10 (2016), pp. 244-247.
[13] P. Bowlan, E. Martinez-Moreno, K. Reimann, T. Elsaesser, M. Woerner. "Ultrafast terahertz response of multilayer graphene in the nonperturbative regime", Phys. Rev. B, 89 (2014), 041408. I ' 48
[14] S. A. Smolyansky, D. V. Churochkin, V. V. Dmitriev, A. D. Panferov, B. Kampfer. "Residual currents generated from vacuum by an electric field pulse in 2+1 dimensional QED models", EPJ Web Conf., 138 (2017), 06004. d t
48 55
[15] I. Gierz, J.C. Petersen, M. Mitrano, C. Cacho, I.C. Edmond Turcu, E. Springate, A. Stohr, A. Kohler, U. Starke, A. Cavalleri. "Snapshots of non-equilibrium Dirac carrier distributions in graphene", Nature Materials, 12 (2013), pp. 1119-1124. 48
[16] K. J. Tielrooij, L. Piatkowski, M. Massicotte, A. Woessner, Q. Ma, Y. Lee, K. S. Myhro, C. N. Lau, P. Jarillo-Herrero, N. F. van Hulst, F. H. L. Koppens. "Generation of photovoltage in graphene on a femtosecond timescale through efficient carrier heating", Nature Nanotechnology, 10 (2015), pp. 437-443.
^48,49
[17] D.B. Blaschke, A.V. Prozorkevich, G. Röpke, C.D. Roberts, S.M. Schmidt, D.S. Shkirmanov, S.A. Smolyansky. "Dynamical Schwinger effect and high-intensity lasers. Realising nonperturbative QED", Eur. Phys. J. D, 55 (2009), 341.'1 t«
[18] M. Galassi, J. Theiler et all. GNU Scientific Library—GSL 2.5 documentation, https://www.gnu.org/software/gsl/, 1996-2018. t50
[19] S. A. Levenec, T. T. Verevin, A. V. Makhankov, A. D. Panferov, S. O. Pirogov. "Modeling the dynamics of massless charge carriers in a two-dimensional system", Computer science and information technology, Proceedings of the international scientific conference (July 2 - July 3, 2018, Saratov, Russia), Publ. center "Science", 2018 (Russian), (ui^tsi 55
Received
Revised
Published
19.12.2018
12.02.2019 26.03.2019
Recommended by
p.h.d. Sergey A. Amelkin
Sample citation of this publication:
Anatolii Panferov, Alexey Makhankov. "Simulation of the effect of short optical pulses on graphene". Program, Systems: Theory and Applications, 2019, 10:1(40), pp. 47-58. 10.25209/2079-3316-2019-10-1-47-58
URL http : //psta. psiras. ru/read/psta2019_l_47-58 .pdf
The same article in Russian: 10.25209/2079-3316-2019-10-1-33-46
About the authors:
Anatolii Dmitrievich Panferov
PhD in Physics and Mathematics, Deputy Head PRCNIT SSU on research and production activities. Associate Professor of the Department Discrete Mathematics and Information Technology SSU. Research interests: high-performance computing, parallel programming, numerical solution of quantum kinetic equations, modeling of the vacuum production of particles in QED, generation of carriers in semiconductors including gapless, processes in the early stages of collision of relativistic nuclei.
[Da 0000-0003-2332-0982 e-mail: panferovad@info.sgu.ru
Alexey Vladimirovich Makhankov
Graduate student of the SSU. Research interests: modeling of physical processes on high-performance computing systems, parallel programming.
e-mail:
0000-0002-9848-9734
next642009@yandex.ru
Эта же статья по-русски:
10.25209/2079-3316-2019-10-1-33-46
