Numerical modeling of stationary pseudospin waves on a graphene monoatomic films

Full Text

1. Introduction One of the main areas of the theoretical and experimental research is the study of the properties of graphene, an obvious candidate for the formation on its basis of the elemental base of future nanoelectronics and spintronics, which provides a gain by orders of magnitude in the field of speed, size and power consumption of devices for storing, transmitting and processing information. Significant successes were achieved in this area [1]-[7], however, a number of very unusual for carbon structures properties of graphenes, which are observed experimentally, have not yet found their adequate theoretical description. In particular, it was found [2], [3] that in graphene systems there is a ferromagnetic effect that present right up to room temperatures and above (Curie temperature exceeds 500 K). Consequently, samples of graphene films can have their own magnetization due to the presence of a nonzero spin density of valence electrons distributed in some way on the two-dimensional carbon lattice. According to researchers, these experimentally observed ferromagnetic properties of graphene structures require theoretical justification and the construction of an appropriate mathematical model. For the first time, the theoretical model of the spin-electron structure of the single-layer graphene film was proposed by Wallace. After it was published an extensive bibliography, which we will not concern. The ferromagnetic properties of graphene structures observed experimentally by different researchers, by their own admission, require proper justification and construction of the appropriate theoretical model. We believe that the source of ferromagnetism is the spontaneous breaking of spin symmetry in a graphene film. Quantum-chemical calculations show that the Wallace model can serve only as a first approximation, and subsequent approximations allow violation of the well-known Wallace symmetry. Quantum-chemical simulation of electron density in a monoatomic graphene film by the extended functional density method and the advanced Hartree- Fock method has demonstrated the possibility of the existence of unpaired electrons, which gives rise to spontaneous violation of spin symmetry in it, i.e., a nontrivial distribution of spin density. Experimentally it was found that such a non-trivial distribution does not have a traditional source: impurities, defects, boundaries [2]. Being caused by spontaneous symmetry breaking, such a distribution of spin density must satisfy the nonlinear phenomenological equation of the classical gauge field [8]-[14]. There can be many such fields and equations, and they all supposedly give results that coincide in a first approximation. Therefore, the scalar field on the two-dimensional continuum 4 was chosen as the first object of study as a mathematical model of spin density. Here we made a transition from a discrete set of nodes of a double hexagonal lattice, in which unpaired electrons and the corresponding electron and spin densities can be localized. The nonlinear field model proposed for the distribution of the spin density of valence electrons in a graphene film makes it possible to describe the experimentally observed ferromagnetic properties of such films. It is shown that these solutions (kinks, breathers) allow the formation of some spatially localized magnetization density configurations on the surface of a graphene film. The specified scalar field allows localized solutions: kinks and anti-kinks. The combination of kink and anti-kink is a quasiparticle, i.e. an approximate solution of the nonlinear equation for a 4 scalar field. Quantitative estimates are obtained for the energy and spatial dimensions of such a configuration. Their characteristic size was tens of nanometers. It is also shown that such configurations can constitute groups of discrete spectra. The problem of the interaction of kinks and anti-kinks with each other is considered. We also find it interesting to further consider the interaction of breathers with each other and with other physical fields, as well as the dynamics of spinons on graphene (fullerene, nanotube) nonplanar surfaces of various topologies. This work is devoted to a numerical study of the model 4 [15], [16]. Within the framework of the proposed model, the problems of approximate calculation of potential fields, approximate solution of the Schrödinger equations and simulation of control of the external field of population levels [8] are solved. 2. Model description In this work, we study a theoretical model that describes the properties of graphene monoatomic layers that form some two-dimensional surfaces associated with the presence of a nonzero spin density distribution function on these surfaces formed as a result of spontaneous violation of the spin symmetry of valence electrons of carbon atoms on these surfaces. Since the spin density is proportional to the magnetization density, this model allows us to describe the ferromagnetic properties of graphene structures. In the framework of the proposed model, a transition is made from the consideration of a discrete two-dimensional carbon lattice forming a graphene film to a continuous two-dimensional surface stretched over this lattice. The indicated two-dimensional surface is the configuration space of the model. Thus, we are making the transition to a continuum field model. In the model under consideration, a trivial, identically equal to zero, spin density configuration is admissible. However, it was experimentally established that this symmetric field configuration can be spontaneously disturbed to some physically observable. Within the framework of the two-dimensional field model under study, there is an analogue of the Goldstone theorem, known in quantum field theory, according to which an each broken generator of the initial symmetry of the field system must correspond to a massless scalar uncharged boson, which in our case is appropriate to call as a spinon. In this case, the spontaneous breaking of spin symmetry within the framework of the proposed model should lead to the presence on the graphene surfaces of quasiparticle spinons, which are vector bosons in 3-dimensional physical space and scalar pseudo-Goldstone bosons in the two-dimensional configuration space of the model, since the projection of the quasiparticle spin on the configuration space is equal to zero. It is significant that the presence of collective magnetic interactions of spinons, due to the influence of the total magnetic field, created by all spinons, on each spinon individually, leads to nonlinearity of the corresponding field equations and, as a consequence, the possibility of the existence of soliton configurations on graphene surfaces, which depend, inter alia, on shape and surface topology. In addition, the presence of collective interactions in the ensemble of spinons should lead to the appearance of an effective mass for a spinon, which should also affect the observed physical consequences, although due to the smallness of spin-spin interactions large values of this mass can hardly be expected. Based on the foregoing, it is clear that the equations for the desired scalar field given on a certain two-dimensional surface must be nonlinear and defined on this surface of an arbitrary, generally speaking, shape and topology. The form and topology in this case determine the boundary conditions for the field function. The indicated function determines the conditions for the existence, configuration and dynamics of quasiparticles of this field on a given two-dimensional surface. The indicated properties, in particular, are possessed by the field equations known in quantum theory, which describe, among other things, massless nonlinear scalar excitations. Thus, to describe spinon excitations on graphene surfaces, we use one of the variants of the nonlinear field model, which allows us to calculate the eigen solutions, effective masses, topological invariants, energy spectra, the dynamics of various nonlinear spinon configurations, as well as relaxation properties, Curie temperature, and other characteristics of statistical ensembles of spinons. 3. Interaction model Consider a model of a nonlinear one-component scalar field on a twodimensional surface, the surface density of the Lagrangian of which is set in the form: () = 1 () - 2 0 (2 - 2)2, (1) 4 where 0 and are the model parameters. In this case, the field equation has the form: 0 [ - 2] - 3 = 0 (2) In the case when depends only on one coordinate and does not depend on time, in other words: = (), and the equations of the form (2) have a set of static solutions, as well as kink and anti-kink solutions: 0 ±() = ±0 tanh (√2/2). (3) In the vicinity of zero, kinks (anti-kinks) have a domain wall separating regions with magnetization of different signs. The kink energy per unit length in the coordinate is calculated by the formula: () = ∞ ∫ d [ 1 ( )2 + (2 - 2)2] = √ 3 √ 3 . (4) 2 4 -∞ 0 8 0 The kink energy density along the coordinate is proportional to 0 [1 - tanh2 (√2/2)] and, therefore, is concentrated near zero on the domain wall. As in any ferromagnet, in the system under consideration there is a Curie temperature [1], [2], at which the system of interacting spins is disordered due to thermal motion. This leads to the decay of the kink and the destruction of the domain structure. In other words, this is the situation when the energy of the thermal motion of an elementary magnetic moment becomes comparable with its energy in the kink field. This allows a quantitative assessment of the model parameters. The substitution of the numerical values gives ≈ 15-30 nm, which looks quite plausible. The scatter of values for the domain wall thickness is related to the scatter of the available experimental data on measuring the Curie temperature. In any case, we see that the calculated domain wall thickness is tens of bond lengths in the cell. This confirms the correctness of the use of the proposed continuum model for the graphene lattice. Consider the case when the solutions for the field function are explicitly time-dependent. Then we obtain stable single kink and anti-kink solutions propagating along the coordinate with a constant velocity 0. The stable spatially localized field configurations (not only kinks) (in particular, which are solutions of equations of the form (2)) are of interest for many practical applications. For qualitative estimates we propose to use some approximate solutions of equations of the form (4), by using combinations of existing exact solutions. +In particular, we can consider the system of kink and +anti-kink interaction. In the simplest case, this can be a kink and an anti-kink located at some distance () from each other. Moreover, it is very important that kink and anti-kink interact with each other even at an infinitely large distance from each other. This is due precisely to the fact that their asymptotics at spatial infinity are nonzero. In addition, it should be noted that, due to the nonlinearity of the problem, the sum of exact solutions, generally speaking, is not an exact solution. Nevertheless, we choose the field function of the system of interacting kink and anti-kink in the simple form: (, ) = [+( + ) + -( - ) - 0], > 0. (5) A function of the form (5) for small values of is spatially localized near the point = 0. Kink and anti-kink at (5), spatially separated by a sufficiently large (compared to the thickness of the kink itself) distance, interact with each other, but,however, steadily maintain their own shape. We numerically look for such () to satisfy the equation (2). The obtained solutions will correspond to breathers, that is, stable kink - anti-kink configurations, known, for example, for equations like sine-Gordon, Korteweg-de Vries, and some others, both in discrete and continuous cases. Consider the Hamiltonian of a system whose field function of the form (5) satisfies the equation of type 2: () = ∞ ∫ d [()() + (()2 - 2) 2 2 0 ]. (6) -∞ This function can be considered as the total energy of the kink in the anti-kink field (or vice versa), and its dependence on the parameter may be formally investigated. Then the dependence of the Hamiltonian of the form (6) on the parameter corresponds to the dependence of the potential interaction energy of the kink and anti-kink on the distance between them. If there are minima in this function, one should expect the presence of bound states in the kink-anti-kink system. These will be the desired breathers. We will search for the bound states of the kink-anti-kink in the minima of potential energy: (, ) = 0, 2 2 (, ) > 0. (7) You can get the estimate: 0 ≈ 0.8. Thus, the ratio of the distance between the kink and the anti-kink to the size of the kink itself (anti-kink) is less than unity, which again confirms the correctness of the assumptions used. 3.1. Numerical implementation Since we consider functions (, ) of the form (5), depending on the parameter 0 , where ±() = ±0 tanh (√2/2) depends on the parameter 0 , the calculation of the functional (, ) was performed as follows. First, the partial derivative was calculated (the other partial derivatives of are identically equal to zero due to the choice of the form of the function (5)), then for given 0, and the numerical integral (6) is calculated. Thus, with the fixed parameter 0 = 1.0 selected, we obtain a numerical dependence of the potential interaction energy (, ) The calculations were performed for various = 0.01 … 100.0. Moreover, for each value of the dependence is () ≡ (, ) was calculated for the values of from = +0.0 to the value of () for which the values of (()) and (2()) coincide with up to the 5th digit. Since the dynamics of the breather in the model of spontaneous symmetry breaking [8], [9] is described by quantum nonrelativistic equations, we model the dynamics of the pulsation of stationary states of the breather by the Schrödinger equation. 4. Dynamics model For the quantum-mechanical stationary wave function of the breather (), the Schrödinger equation in the standard form is written: ℏ2 2 [- (, )]() = (), (8) 2 2 where is the effective mass of the breather, which is equal in this case to the sum of the masses of the free kink and anti-kink, and (, ) is the potential part of the total energy of the breather, depending on , is the energy of the corresponding stationary state. The movement of the breather along the generalized coordinate physically corresponds to a change in the distance between the kink and the anti-kink, while the center of the breather is fixed along the coordinate , only the effective width of the breather changes (the breather “breathes”). 4.1. Numerical implementation The dynamics of the breather is described by the Schrödinger equation: d2 - d2 + ()() = (). (9) We will seek a numerical solution of the equation (9) by using the Ritz- Galerkin method. Namely, under the assumption that the Schrödinger operator ̂ is strictly positive, we can construct the Ritz functional and the energy space in which it reaches its minimum. The local minima of the Ritz functional uniquely correspond to the eigenvectors of the operator ̂. The finite-dimensional approximation of the functional over some complete system of coordinate functions () minimizes the function of several variables. The minimum of this function corresponds to the solution of the system of linear algebraic equations with the Ritz matrix (̂, ). Te application of the Galerkin method to this problem leads to the eigenvalue and eigenvector problem, as does the Ritz method. You can use an arbitrary orthonormal basis in the space 2 of quadratically integrable functions on the axis as a complete system of coordinate functions. Moreover, convergence may turn out to be slow, which will lead to the need to solve systems of linear algebraic equations of large dimension. In order to build a well-conditioned Ritz matrix, one should choose a strongly minimal or almost orthogonal complete system of functions with respect to ̂. By the way, one should note that the potential energy in the equation (9) is a grid function, which complicates the analytical study of the known bases in 2. At the same time, the quantum operator of a nonlinear oscillator from [17] has a complete system of eigenvectors. The eigenfunctions of a nonlinear oscillator form a numerically strongly minimal system of functions for the operator ̂. These functions, however, are not orthonormal in the space 2. This fact complicates the finite-dimensional approximating problem and, more precisely, leads to the generalized eigenvalue problem and eigenvectors ̂ ⃗ = ̂ .⃗ Then, in the general case, the solution will look like this: ∞ () = ∑ ,(), =1 where () are basic functions. Therefore, in order to find the solution (9), it is necessary to specify some kind of functions () and find the values of the coefficients , and . From the obtained the numerical values of () it is clear that the form of the potential is similar to the potential of a nonlinear oscillator, so we take the eigenfunctions () for it. Consider the case = 0, 1, 2. Choose = 1. In this case, = 3, but for a smaller calculation error, you can choose a value of more than 3. Let = 10. For > 0 the condition must be satisfied: 20 ⩾ . It follows from this inequality that in our case 0 ⩾ 5. Consider 0 = 6. Subsequent are calculated using the formula = -1 - . Therefore, 1 = 5, 2 = 4. Also, 2 2 to find the form of the functions () we need the values = . In our case, 0 = 3, 1 = 5 , 2 = 2. We look for the form of the functions () using formulas: ⎧ 1 where {{0 () = (1 + 2)0 , { ⎨1 () = + (0) 0 (1) , { ⎩2 () = + (0) + (1) 0 (2) . (10) √ + = 1 (-√1 + 2 + √ ) . 2 d 1 + 2 +After substitution of the +values , , and simplification of the results, the functions (), = 0, 1, 2 will take the form: ⎧ {0 () = { { { 1 , (1 + 2)3 11 1 () = √ , ⎨ 2 (1 + 2)3 { { { 9 (102 - 1) (11) {2 () = ⎩ 2 (1 + 2)3 . The Ritz system takes form: 2 2 ∑ (, ) = ∑ (, ), = 0, 2, (12) =0 =0 where the operator = - [2] + (). Having written the scalar products (, ) and (, ), we get: ∞ 2 2 2 ∞ d - - ∫ d2 (∑ ()) () + ∑ ∫ () () () 0 =0 =0 0 2 ∞ - ∑ ∫ () () d = 0. =0 0 The same equation can be written in the matrix form: ( + - ) ⃗ = 0,⃗ (13) where = (), , = 0, 2, 2() ∞ 2() d = (- 0 d2 , ()) = - ∫ d2 () ∞ , = (), , = 0, 2, = ((), ()) = ∫ ()() d, 0 ∞ = (), , = 0, 2, = ( ()(), ()) = ∫ ()()() d, 0 ⃗ = (0, 1, 2). There are four unknowns in the resulting system of equations: and coordinates ⃗ = (0, 1, 2). In order for the system to be compatible, it is necessary that | + - | = 0. The expanded form: ∣ 00 + 00 - 00 10 + 10 - 10 20 + 20 - 20 ∣ ∣ 01 + 01 - 01 11 + 11 - 11 21 + 21 - 21 ∣ 02 + 02 - 02 12 + 12 - 12 22 + 22 - 22 We introduce the following replacements: ∣ ∣ = 0. (14) ∣ ∣ 1 = 00 + 00, 2 = 11 + 11, 3 = 22 + 22, 4 = 02 + 02, 5 = 10 + 10, 6 = 21 + 21, 7 = 20 + 20, 8 = 01 + 01, 9 = 12 + 12. (15) We get the equation: (1 - 00)(2 - 11)(3 - 22) + (4 - 02)(5 - 10)(6 - 21) + + (7 - 20)(8 - 01)(9 - 12) - (7 - 20)(2 - 11)(4 - 02) - -(1 -00)(2 -12)(6 -21) -(3 -22)(5 -10)(8 -01) = 0. Opening the brackets, we obtain the equation of the third degree with respect to . The resulting numerical values are: = (0.6362135, -10.452091, 24.544966). (16) 5. Conclusion Thus, we see that the kink plus anti-kink system can form bound states, that is, breathers which make up some groups corresponding to different minima of the interaction energy of the kink and anti-kink. These groups of breathers should noticeably differ in spatial dimensions. It is very likely that the energy spectra of breathers from various groups can overlap. This, in turn, allows us to pose the problem of tunneling breathers from one minimum to another, the lifetimes in each of these minima, the behavior in external fields, inverse populations, and other physical features of another, as well as the problem of the lifetimes in each of these minima and of the behavior in external fields, inverse populations, and some other physical features of the system under consideration. We also note that under certain conditions, the interaction of kinks and breathers with each other, as well as with external fields (note, not only electromagnetic, but also, for example, acoustic), should lead to the birthdestruction of particles, which, in principle, we can describe within the framework of the proposed model by secondary quantization of the considered system and calculating the matrix elements of the corresponding scattering matrix. Thus, in some sense, we are closing the circle, moving from a general quantum-field approach to a classical field, then to a quantum-mechanical model, and again to a quantum-field consideration through second quantization. It seems that the proposed approach will allow us to continue to make the necessary quantitative estimates even before the numerical simulation of the problem. We also find it interesting to further consider the dynamics of spinons on graphene (fullerene, nanotube) non-planar surfaces of various topologies.

About the authors

Lê Anh Nhật

Peoples’ Friendship University of Russia (RUDN University)

Author for correspondence.
Email: leanhnhat@tuyenquang.edu.vn

PhD student of Department of Applied Probability and Informatics

6, Miklukho-Maklaya St., Moscow 117198, Russian Federation

Konstantin P. Lovetskiy

Peoples’ Friendship University of Russia (RUDN University)

Email: lovetskiy-kp@rudn.ru

Docent, PhD in Physics and Mathematics, Associate Professor at the Department of Applied Probability and Informatics

6, Miklukho-Maklaya St., Moscow 117198, Russian Federation

Leonid A. Sevastianov

Peoples’ Friendship University of Russia (RUDN University); Bogoliubov Laboratory of Theoretical Physics Joint Institute for Nuclear Research

Email: sevastianov-la@rudn.ru

Professor, Doctor of Sciences in Physics and Mathematics, Professor at the Department of Applied Probability and Informatics

6, Miklukho-Maklaya St., Moscow 117198, Russian Federation; 6, Joliot-Curie St., Dubna, Moscow region 141980, Russia

Dmitry S. Kulyabov

Peoples’ Friendship University of Russia (RUDN University); Laboratory of Information Technologies Joint Institute for Nuclear Research

Email: kulyabov-ds@rudn.ru

Docent, Doctor of Sciences in Physics and Mathematics, Professor at the Department of Applied Probability and Informatics

6, Miklukho-Maklaya St., Moscow 117198, Russian Federation; 6, Joliot-Curie St., Dubna, Moscow region 141980, Russian Federation

References