We study dynamical phenomena in a harmonic graphene (honeycomb) lattice, consisting of equal particles connected by linear and angular springs. Equations of in-plane motion for the lattice are derived. Initial conditions typical for molecular dynamic modelling are considered. Particles have random initial velocities and zero displacements. In this case, the lattice is far from thermal equilibrium. In particular, initial kinetic and potential energies are not equal. Moreover, initial kinetic energies (and temperatures), corresponding to degrees of freedom of the unit cell, are generally different. The motion of particles leads to equilibration of kinetic and potential energies and redistribution of kinetic energy among degrees of freedom. During equilibration, the kinetic energy performs decaying high-frequency oscillations. We show that these oscillations are accurately described by an integral depending on dispersion relation and polarization matrix of the lattice. At large times, kinetic and potential energies tend to equal values. Kinetic energy is partially redistributed among degrees of freedom of the unit cell. Equilibrium distribution of the kinetic energies is accurately predicted by the non-equipartition theorem. Presented results may serve for better understanding of the approach to thermal equilibrium in graphene.

This article is part of the theme issue ‘Modelling of dynamic phenomena and localization in structured media (part 2)’.

1. Introduction

Two-dimensional materials attract significant interest from researchers due to their extraordinary physical properties [1,2]. The best-known two-dimensional material is graphene, a monolayer of carbon atoms forming a perfect hexagonal (honeycomb) lattice [3,4]. Recent advances in nanoscale manufacturing allow the fabrication of single-layer graphene of submicrometre size [5]. The perfect structure of graphene yields unique physical properties. A review of theoretical and experimental studies of graphene's mechanical and thermal properties can be found in [6,7]. In particular, in paper [8], it is shown that the strength of graphene is close to the so-called theoretical strength. Graphene also shows extraordinary heat-conducting properties, including high thermal conductivity (above 3000 W m K⁻¹, [9]). Also, although in-plane elastic properties of graphene are isotropic [10], heat transport is anisotropic [11,12].

The low concentration of defects in graphene makes it a promising candidate for investigation of the ballistic heat transport [13–16]. In contrast with conventional diffusive heat transport, in the ballistic regime, the heat is carried by waves, travelling with velocities close to the speed of sound. This phenomenon is of great interest from both theoretical and applied points of view [13]. However, a theory of ballistic heat transport is still under development. In ballistic heat transport experiments, an initial temperature distribution can be generated by ultrashort laser excitation [17–21]. Laser excitation creates a strongly non-equilibrium state in the material. Therefore, it is important to develop mathematical models describing the transition from the non-equilibrium state to thermal equilibrium.

A convenient model for the description of thermal properties of solids is a harmonic crystal, i.e. a set of particles forming a perfect lattice and interacting via linearized forces. An approach to thermal equilibrium in harmonic crystals is considered in many works [22–32]. In particular, in a pioneering paper by Klein & Prigogine [26], the behaviour of kinetic temperature, proportional to kinetic energy of thermal motion, in a one-dimensional harmonic chain with random initial conditions. An expression for the temperature is obtained using an exact analytical solution of equations of motion. Another approach based on analysis of covariances of particle velocities is proposed by Krivtsov [27]. In papers [28–30,33–36], this approach is applied to crystals with monoatomic lattice. It is shown that the approach to thermal equilibrium is accompanied by two physical processes: equilibration of kinetic and potential energies and redistribution of energy among degrees of freedom of the unit cell. Generalization for the case of complex lattices with an arbitrary number of particles per unit cell is carried out in paper [37]. In the present paper, an approach to thermal equilibrium in a two-dimensional graphene lattice is considered.

The main goal of the present paper is to investigate peculiarities of thermal equilibration in graphene. The paper is organized as follows. In §2, the geometry of the lattice is described. In §3, equations of motion of the graphene lattice are derived using the Lagrange–Euler formalism. Equations of motion are obtained using a two-parametric model. In the framework of this model, atoms interact via longitudinal and angular springs. It is the simplest model allowing a correct description of in-plane elastic properties of graphene [38] as well as other honeycomb structures [39]. Random initial conditions, typical for molecular dynamics modelling of thermal processes, are discussed in §4. In §5, a dispersion relation is derived. Exact formulae, describing the time evolution of kinetic and potential energies, are presented in §6. In §7, analytical predictions are compared with results of numerical solution of equations of motion. Redistribution of kinetic energy between degrees of freedom of the unit cell¹ is investigated in §8.

2. Geometry of the lattice

Graphene has a regular honeycomb lattice shown in figure 1a. Its unit cell contains two carbon atoms (figure 1b). Unit cells are numbered by pairs of indices {n, m}. Displacements of atoms from unit cell {n, m} are denoted as U_n,m, V_n,m. Position vectors of unit cells {n, m} and {s, p} in the reference configuration are related as

r_{n, m} = r_{s, p} + (n - s) e_{1} + (m - p) e_{2} .

2.1

Here, primitive vectors of the lattice, e₁, e₂, are represented in Cartesian basis as

e_{1} = a \frac{\sqrt{3}}{2} (i + \sqrt{3} j) and e_{2} = a \frac{\sqrt{3}}{2} (- i + \sqrt{3} j),

2.2

where a is an equilibrium bond length in graphene; in figure 1b, vector i is horizontal, vector j is vertical.

Vectors a₁, a₂, a₃ and a₄, a₅, a₆ connect atoms with the displacements U_n,m and V_n,m with their nearest neighbours (figure 1). The vectors satisfy relations a₄ = − a₁, a₅ = − a₂, a₆ = − a₃.

3. Equations of motion

In this section, the equations of motion for two atoms in a unit cell are derived using the Lagrange–Euler formalism. Each atom has 2 d.f., i.e. in-plane motions are considered. Since in harmonic approximation in-plane and out-of-plane motions are decoupled, the latter can be considered separately (see [37,40]).

The following model of interactions in graphene is considered. Each atom is connected with three nearest neighbours by linear springs (bonds) with stiffness c. Additionally, the nearest bonds between particles are connected by angular springs with stiffness g.

In the framework of the Lagrange–Euler formalism, equations of motion of the unit cell have the form:

\frac{d}{d t} \frac{\partial K}{\partial {\dot{q}}_{n, m}^{i}} - \frac{\partial K}{\partial q_{n, m}^{i}} = - \frac{\partial W}{\partial q_{n, m}^{i}}, q_{n, m} = {[\begin{matrix} U_{n, m}^{x} & U_{n, m}^{y} & V_{n, m}^{x} & V_{n, m}^{y} \end{matrix}]}^{⊤},

3.1

where U^x_n,m, U^y_n,m, V^x_n,m, V^y_n,m are components of displacement vectors U_n,m, V_n,m in Cartesian basis; qⁱ_n,m, i = 1, 2, 3, 4 is the i-th component of column q_n,m; K is the total kinetic energy of the unit cell; W is contribution of particles with displacements U_n,m, V_n,m to the total potential energy of the lattice (see formulae (3.3), (3.4), (3.5)); → p stands for the transpose sign. Here and below matrices are denoted by bold italic symbols, while invariant vectors, e.g. displacements, are denoted by bold symbols.

The total kinetic energy of the unit cell is calculated as

K = \frac{1}{2} M ({\dot{U}}_{n, m}^{2} + {\dot{V}}_{n, m}^{2}),

3.2

where M is particle mass.

Potential energy, W, is equal to the sum of energies of all linear and angular springs connected with particles from the unit cell {n, m}:

W = W^{s} + W^{b} .

3.3

Here, W^s and W^b stand for contributions of linear and angular springs respectively. The expression for W^s in harmonic approximation has the form

\begin{aligned} W^{s} & = \frac{c}{2 a^{2}} [(a_{1} \cdot (V_{n + 1, m} - U_{n, m}))^{2} + (a_{2} \cdot (V_{n, m + 1} - U_{n, m}))^{2} \\ + (a_{3} \cdot (V_{n, m} - U_{n, m}))^{2} + (a_{4} \cdot (U_{n - 1, m} - V_{n, m}))^{2} + (a_{5} \cdot (U_{n, m - 1} - V_{n, m}))^{2}] . \end{aligned}

3.4

Here, dot stands for the scalar product; a is an equilibrium length of linear springs defined by formula (2.2).

The angular part of the potential energy, W^b, depends on the change of 14 angles shown in figure 2. The following expression for W^b is used

W^{b} = \sum_{i = 1}^{14} Π_{i}^{b} and Π_{i}^{b} = \frac{1}{2} g a^{2} (Θ_{i} - Θ_{i}^{0})^{2} .

3.5

Since calculation of forces in the case of multibody interactions is not trivial [41], we show detailed derivation for a single angular spring.

Figure 2. Numbering of angles for the given unit cell.

Consider, for example, angular spring number 2 connecting bonds with vectors a₁ and a₂ (figures 1 and 2). Change of angle Θ₂ between the bonds is caused by relative rotation of vectors a₁ and a₂ due to motion of particles. In a deformed state, these vectors turn into vectors a₁ + V_n+1,m − U_n,m and a₂ + V_n,m+1 − U_n,m. Corresponding rotation angles are denoted ϕ₁, ϕ₂. The angles are calculated as

\begin{aligned} \sin ϕ_{1} & = \frac{(a_{1} \times (V_{n + 1, m} - U_{n, m})) \cdot n}{| a_{1} | | a_{1} + V_{n + 1, m} - U_{n, m} |}, \sin ϕ_{2} = \frac{(a_{2} \times (V_{n, m + 1} - U_{n, m})) \cdot n}{| a_{2} | | a_{2} + V_{n, m + 1} - U_{n, m} |} \Rightarrow \\ \Rightarrow ϕ_{1} \approx \frac{1}{a^{2}} (n \times a_{1}) \cdot (V_{n + 1, m} - U_{n, m}), ϕ_{2} \approx \frac{1}{a^{2}} (n \times a_{2}) \cdot (V_{n, m + 1} - U_{n, m}), \end{aligned}

3.6

where n = i × j is a unit vector, normal to graphene plane; × stands for the cross product. Then a change of angle Θ₂ is calculated as the angle of relative rotation ϕ₁ − ϕ₂. Therefore, the expression for the potential energy of the angular spring number 2 takes the form

Π_{2}^{b} = \frac{g}{2 a^{2}} {[(n \times a_{1}) \cdot (V_{n + 1, m} - U_{n, m}) - (n \times a_{2}) \cdot (V_{n, m + 1} - U_{n, m})]}^{2} .

3.7

Energies of other angular springs are calculated similarly.

Substituting expressions for kinetic and potential energies into the Lagrange–Euler equation (3.1), yields equations of motion for the unit cell

\begin{aligned} \begin{aligned} M {\ddot{U}}_{n, m}^{x} & = \frac{c \sqrt{3}}{4} [\sqrt{3} V_{n + 1, m}^{x} + V_{n + 1, m}^{y} - 2 \sqrt{3} U_{n, m}^{x} + \sqrt{3} V_{n, m + 1}^{x} - V_{n, m + 1}^{y}] \\ + \frac{g}{4} [2 (- U_{n - 1, m}^{x} + \sqrt{3} U_{n - 1, m}^{y} - U_{n, m - 1}^{x} - \sqrt{3} U_{n, m - 1}^{y}) + U_{n + 1, m - 1}^{x} \\ + \sqrt{3} U_{n + 1, m - 1}^{y} + U_{n - 1, m + 1}^{x} - \sqrt{3} U_{n - 1, m + 1}^{y} - 2 (U_{n + 1, m}^{x} + U_{n, m + 1}^{x}) \\ + 6 (V_{n + 1, m}^{x} - \sqrt{3} V_{n + 1, m}^{y} - 5 U_{n, m}^{x} + 4 V_{n, m}^{x} + V_{n, m + 1}^{x} + \sqrt{3} V_{n, m + 1}^{y})], \\ M {\ddot{U}}_{n, m}^{y} & = \frac{c}{4} [\sqrt{3} V_{n + 1, m}^{x} + V_{n + 1, m}^{y} - 6 U_{n, m}^{y} + 4 V_{n, m}^{y} - \sqrt{3} V_{n, m + 1}^{x} + V_{n, m + 1}^{y}] \\ + \frac{\sqrt{3} g}{4} [U_{n - 1, m + 1}^{x} - \sqrt{3} U_{n - 1, m + 1}^{y} - U_{n + 1, m - 1}^{x} - \sqrt{3} U_{n + 1, m - 1}^{y} + 2 U_{n + 1, m}^{x} \\ - 2 U_{n, m + 1}^{x} - 10 \sqrt{3} U_{n, m}^{y} + 6 (V_{n, m + 1}^{x} + \sqrt{3} V_{n, m + 1}^{y} - V_{n + 1, m}^{x} + \sqrt{3} V_{n + 1, m}^{y})] \end{aligned}\} \end{aligned}

3.8

and

\begin{aligned} \begin{aligned} M {\ddot{V}}_{n, m}^{x} & = \frac{\sqrt{3} c}{4} [\sqrt{3} U_{n - 1, m}^{x} + U_{n - 1, m}^{y} - 2 \sqrt{3} V_{n, m}^{x} + \sqrt{3} U_{n, m - 1}^{x} - U_{n, m - 1}^{y}] \\ + \frac{g}{4} [6 (U_{n - 1, m}^{x} - \sqrt{3} U_{n - 1, m}^{y} + U_{n, m - 1}^{x} + \sqrt{3} U_{n, m - 1}^{y} + 4 U_{n, m}^{x}) \\ - 2 (V_{n + 1, m}^{x} - \sqrt{3} V_{n + 1, m}^{y} + V_{n, m + 1}^{x} + \sqrt{3} V_{n, m + 1}^{y}) + V_{n + 1, m - 1}^{x} - \sqrt{3} V_{n + 1, m - 1}^{y} \\ + V_{n - 1, m + 1}^{x} + \sqrt{3} V_{n - 1, m + 1}^{y} - 30 V_{n, m}^{x} - 2 (V_{n, m - 1}^{x} + V_{n - 1, m}^{x})], \\ M {\ddot{V}}_{n, m}^{y} & = \frac{c}{4} [\sqrt{3} U_{n - 1, m}^{x} + U_{n - 1, m}^{y} + 4 U_{n, m}^{y} - 6 V_{n, m}^{y} - \sqrt{3} U_{n, m - 1}^{x} + U_{n, m - 1}^{y}] \\ + \frac{\sqrt{3} g}{4} [2 \sqrt{3} (- \sqrt{3} U_{n - 1, m}^{x} + 3 U_{n - 1, m}^{y} - 5 V_{n, m}^{y} + \sqrt{3} U_{n, m - 1}^{x} + 3 U_{n, m - 1}^{y}) \\ + V_{n + 1, m - 1}^{x} - \sqrt{3} V_{n + 1, m - 1}^{y} - V_{n - 1, m + 1}^{x} - \sqrt{3} V_{n - 1, m + 1}^{y} - 2 (V_{n, m - 1}^{x} - V_{n - 1, m}^{x})] . \end{aligned}\} \end{aligned}

3.9

In further calculations, we use the following values of parameters for graphene [42]:

c = 730.2 {N m}^{- 1} g = 66.9 {N m}^{- 1}, a = 0.142 n m .

3.10

Thus equations of motion of the unit cell are given by the formulae (3.8), (3.9), (3.10). In §5, equations (3.8), (3.9) are employed for derivation of the dispersion relation.

4. Initial conditions

In this section, we introduce initial conditions typical for molecular dynamics modelling of thermal processes in lattices [27,28,37]. Initially, particles have random velocities and zero displacements. Under these initial conditions, the system is far from thermal equilibrium.

The initial conditions have the forms

q_{n, m} = 0 and {\dot{q}}_{n, m} = {\dot{q}}_{n, m}^{0} .

4.1

Here,

{\dot{q}}_{n, m}^{0}

is a column of random initial velocities of particles from unit cell {n, m} such that

⟨ {\dot{q}}_{n, m}^{0} ⟩ = 0 and ⟨ {\dot{q}}_{n, m}^{0} {\dot{q}}_{s, p}^{0} ⟩ = A δ_{D} (n - m) δ_{D} (s - p),

4.2

where

⟨ \dots ⟩

stands for mathematical expectation (in computer simulations, mathematical expectation is approximated by average over realizations with different random initial conditions); δ_D(0) = 1; δ_D(n − m) = 0 for n≠m; A is a constant matrix independent of n and m. In other words, components of

{\dot{q}}_{n, m}^{0}

are random numbers with zero mean and generally different variances given by diagonal elements of matrix A. Off-diagonal elements of matrix A are equal to covariances of initial velocities, corresponding to different degrees of freedom of the unit cell. Initial velocities of particles from different unit cells are statistically independent. An example of initial conditions (4.1), (4.2) is given by formula (7.1).

From a macroscopic point of view, initial conditions (4.1), (4.2) correspond to uniform spatial distribution of kinetic temperature, proportional to the kinetic energy of thermal motion. Under these initial conditions, initial kinetic and potential energies are not equal (potential energy is equal to zero). Moreover, kinetic energies, corresponding to different degrees of freedom of the unit cell, are generally different. Time evolution of the energy is considered in §6, 7, 8.

5. Bloch's theorem and dispersion relation

In the present section, we derive the dispersion relation for graphene lattice described by equations of motion (3.8), (3.9).

We seek for the solution of equations of motion (3.8), (3.9) in a form of propagating plane waves:

U_{n, m} = U (k) e^{i (k \cdot r_{n, m} - ω t)} and V_{n, m} = V (k) e^{i (k \cdot r_{n, m} - ω t)},

5.1

where U(k), V (k) are amplitudes of the waves; ω is a frequency; k is a wave vector; r_n,m stands for position vector of the unit cell. The Bloch (or the Bloch-Floquet) theorem states that the displacements of the atoms in a unit cell identified by the integer pair {n + s, m + p} (figure 1) satisfy the relation [43]

U_{n + s, m + p} = U_{n, m} e^{i k \cdot (s e_{1} + p e_{2})} and V_{n + s, m + p} = V_{n, m} e^{i k \cdot (s e_{1} + p e_{2})} .

5.2

As a consequence of the theorem, the proportionate change in the wave amplitude across a unit cell due to a wave propagation does not depend on the location of the unit cell within the lattice. The theorem is widely used in solid-state physics [44] and in engineering analysis of periodic structures [45,46]. Hence, we study wave propagation in the entire lattice by considering motion within a single unit cell.

Using a formula (5.2), we represent equations of motion (3.8), (3.9) in the form

M {\ddot{q}}_{n, m} + C q_{n, m} = 0, M = M E,

5.3

where C and M are 4 × 4 stiffness and mass matrices; E is 4 × 4 identity matrix; vector q_n,m is defined by formula (3.1).

We represent wave vector, k, in reciprocal basis ${\tilde{e}}_{i}$ such that ${\tilde{e}}_{i} \cdot e_{j} = 2 π δ_{i j}$ :

k = k_{1} {\tilde{e}}_{1} + k_{2} {\tilde{e}}_{2}, {\tilde{e}}_{1} = \frac{2 π}{3 a} (\sqrt{3} i + j), {\tilde{e}}_{2} = \frac{2 π}{3 a} (- \sqrt{3} i + j) .

5.4

Here and below k₁, k₂∈[0;1] are dimensionless components of the wave vector. Taking these relations into account, we substitute (5.1) into (5.3) and obtain a homogeneous system of linear equations with respect to q_n,m:

(C (k_{1}, k_{2}) - ω^{2} M) q_{n, m} = 0.

5.5

The system has non-trivial solutions under the following condition:

det (Ω (k_{1}, k_{2}) - ω^{2} E) = 0 and Ω (k_{1}, k_{2}) = \frac{1}{M} C,

5.6

where Ω is referred to as the dynamical matrix [47]. The solution of the fourth order equation (5.6) with respect to ω² yields dispersion surfaces ω_i(k₁, k₂), i = 1, 2, 3, 4. Two acoustic (ω₁, ω₂) and two optical (ω₃, ω₄) surfaces obtained by numerical solution of equation (5.6) are shown in figure 3. The dispersion surfaces are periodic with a unit periodic cell corresponding to the first Brillouin zone shown in figure 4.

Figure 3. Dispersion surfaces for graphene: (a) ω₁ (acoustic), (b) ω₂ (acoustic), (c) ω₃ (optical), (d) ω₄ (optical). Frequencies are normalized by $\sqrt{c / M} .$ (Online version in colour.)

Figure 4. First Brillouin zone of graphene lattice. (Online version in colour.)

Due to the symmetry of the lattice, the irreducible zone, limited by path γ − K − M − γ, covers all possible values of wave frequencies. In the following section, it is shown that dispersion relation and dynamic matrix determine time evolution of kinetic and potential energies during the approach to thermal equilibrium.

6. Time evolution of kinetic energies: an exact solution

In this section, we consider the time evolution of kinetic energies, corresponding to different degrees of freedom of the unit cell. Initially, particle displacements are equal to zero. Therefore, potential energy is equal to zero, while kinetic energy is equal to the total energy. Motion of particles causes two physical processes: redistribution of energy among kinetic and potential forms and redistribution of kinetic energy among degrees of freedom of the unit cell. An analytical description of these processes is presented below. Note that these processes are present in both harmonic and anharmonic systems [26–29,37]. In particular, redistribution of kinetic energy among degrees of freedom of the unit cell is caused by coupling between four equations of motion (3.8), (3.9). The effect of anharmonicity is discussed in §9.

Following paper [37], we introduce the following matrix, characterizing mathematical expectation of kinetic energies of the unit cell:

K = \frac{M}{2} ⟨ {\dot{q}}_{n, m} {\dot{q}}_{n, m}^{⊤} ⟩ and {\dot{q}}_{n, m} = {[\begin{matrix} {\dot{U}}_{n, m}^{x} & {\dot{U}}_{n, m}^{y} & {\dot{V}}_{n, m}^{x} & {\dot{V}}_{n, m}^{y} \end{matrix}]}^{⊤} .

6.1

Diagonal elements of matrix K are equal to kinetic energies, corresponding to degrees of freedom of the unit cell:

\begin{aligned} K_{11} = \frac{1}{2} M ⟨ {({\dot{U}}_{n, m}^{x})}^{2} ⟩, K_{22} = \frac{1}{2} M ⟨ {({\dot{U}}_{n, m}^{y})}^{2} ⟩, K_{33} = \frac{1}{2} M ⟨ {({\dot{V}}_{n, m}^{x})}^{2} ⟩ and K_{44} = \frac{1}{2} M ⟨ {({\dot{V}}_{n, m}^{y})}^{2} ⟩ . \end{aligned}

6.2

Off-diagonal elements of matrix K characterize the correlation between different degrees of freedom. Since initial conditions are spatially homogeneous, K is independent of n, m.

In paper [37], it is shown that the time evolution of matrix K is exactly described by the formula

K = \int_{0}^{1} \int_{0}^{1} P K^{'} P^{* ⊤} d k_{1} d k_{2} and K_{i j}^{'} = \frac{1}{2} {P^{* ⊤} K_{0} P}_{i j} [\cos ((ω_{i} - ω_{j}) t) + \cos ((ω_{i} + ω_{j}) t)] .

6.3

Here, K₀ is an initial value of matrix K; { · s}_ij is element i, j of the matrix; ω_i(k)≥0, i = 1, .., 4 are branches of dispersion relation; P(k) stands for the polarization matrix, composed of normalized eigenvectors of the dynamic matrix (5.6);^* stands for complex conjugate; k₁, k₂ are dimensionless components of the wave vector (5.4).

If initial kinetic energies, corresponding to all degrees of freedom of the unit cell, are equal then the total kinetic energy of the cell has the form

K = tr K = \frac{K_{0}}{2} [1 + \frac{1}{4} \sum_{j = 1}^{4} \int_{0}^{1} \int_{0}^{1} \cos (2 ω_{j} (k_{1}, k_{2}) t) d k_{1} d k_{2}] .

6.4

Here K₀ is an initial value of the total kinetic energy of the unit cell. Given the known kinetic energy, the potential energy per unit cell can be calculated as Π = K₀ − K. At large times, the integral in the right side of formula (6.4) tends to zero. Kinetic energy tends to half of its initial value. Therefore, at large times, mathematical expectations of kinetic and potential energies are equal. This fact is in a good agreement with the results presented in paper [48].

Oscillations of the kinetic energies, described by formulae (6.3), (6.4), decay in time. At large times, the system tends to thermal equilibrium, i.e. to a state with spatially and temporary uniform distribution of kinetic temperature [37]. Kinetic and potential energies tend to some equilibrium values and K tends to K_eq. Equilibrium kinetic energies can be calculated using the non-equipartition theorem, derived in paper [37]:

K_{e q} = \frac{1}{2} \int_{0}^{1} \int_{0}^{1} P diag (P^{* ⊤} K_{0} P) P^{* ⊤} d k_{1} d k_{2},

6.5

where diag( · s) stands for the diagonal part of the matrix. The theorem relates equilibrium distribution of kinetic energy among degrees of freedom with initial conditions.

In the following sections, formulae (6.3), (6.4), (6.5) are employed for investigation of redistribution of energy among kinetic and potential forms and redistribution of kinetic energy among degrees of freedom of the unit cell.

7. Equilibration of kinetic and potential energies

In this section, we consider behaviour of the total kinetic energy of the unit cell. Time evolution of the energy is caused by redistribution of energy among kinetic and potential forms.

To check formulae (6.3), (6.4), we compare them with results of the numerical solution of lattice dynamics equations (3.8), (3.9) with the following initial conditions:

\begin{aligned} U_{n, m} = V_{n, m} = 0, {\dot{U}}_{n, m}^{x} = β_{n, m}^{x} \sqrt{\frac{2 K_{11}^{0}}{M}}, {\dot{U}}_{n, m}^{y} = β_{n, m}^{y} \sqrt{\frac{2 K_{22}^{0}}{M}} \\ and & {\dot{V}}_{n, m}^{x} = γ_{n, m}^{x} \sqrt{\frac{2 K_{33}^{0}}{M}}, {\dot{V}}_{n, m}^{y} = γ_{n, m}^{y} \sqrt{\frac{2 K_{44}^{0}}{M}}, \end{aligned}\}

7.1

where K⁰_ii, i = 1, 2, 3, 4 are proportional to mathematical expectations of initial kinetic energies, corresponding to degrees of freedom of the unit cell; β^x_n,m, β^y_n,m, γ^x_n,m, γ^y_n,m are uncorrelated random values with zero mean and unit variance (in computer simulations, these values are uniformly distributed in the interval [ − 1;1]). Corresponding matrix K₀ is diagonal. In simulations, the graphene sheet consists of 500 × 500 unit cells under fixed boundary conditions. Since we are interested in the behaviour of energies at relatively short times, boundary conditions are not important. Numerical integration of equations of motion is carried out using a symplectic leap-frog integrator with time-step 5 × 10⁻³τ_*, where

τ_{*} = 2 π \sqrt{M / c}

. During the simulation, kinetic energies, corresponding to degrees of freedom of the unit cell, for the whole sheet are calculated. In this case, averaging over realizations is not necessary.

Consider the case where initial kinetic energies corresponding to all degrees of freedom of the unit cell are equal (K⁰₁₁ = K⁰₂₂ = K⁰₃₃ = K⁰₄₄). It can be shown that under these initial conditions, the kinetic energies remain equal at any moment in time, i.e. matrix K remains spherical. Then behaviour of kinetic energy is described by formula (6.4). Comparison with the results of numerical simulations is shown in figure 5. It is seen that analytical solution (6.4) practically coincides with the results of the numerical integration of lattice dynamics equations. The kinetic energy of the system oscillates in time and tends to half of the initial value. Since the total energy of the system is conserved, then at large times kinetic and potential energies become equal.

The oscillations of kinetic energy, shown in figure 5, are quite complicated and contain multiple frequencies. Formula (6.4) allows us to simplify analysis of the oscillations. It shows individual contributions of branches of the dispersion relation. Rewriting formula (6.4), yields

K = \frac{K_{0}}{2} + \sum_{j = 1}^{4} I_{j} and I_{j} = \frac{K_{0}}{8} \int_{0}^{1} \int_{0}^{1} \cos (2 ω_{j} (k_{1}, k_{2}) t) d k_{1} d k_{2} .

7.2

Contributions I_j are presented in figure 6. It is seen that oscillations shown in figure 5 have six main frequencies. Contributions, I₁, I₂, of acoustic branches, ω₁, ω₂, have one main frequency each, while contributions, I₃, I₄, of optical branches, ω₃, ω₄, have a form of beats (two close main frequencies). Note that all main frequencies are larger than

\sqrt{c / M}

. Particular expressions for the frequencies can be derived from asymptotic analysis of corresponding integrals (see e.g. paper [49]). The analysis is beyond the scope of the present paper.

Figure 6. Contributions of branches of dispersion relation to oscillations of the total kinetic energy in graphene. Here, I_j, j = 1, 2, 3, 4 stands for contribution of branch ω_j (see formula (6.4)). (Online version in colour.)

Analysis of formula (6.3) and numerical results also show that behaviour of the total kinetic energy is independent on initial distribution of kinetic energy among degrees of freedom of the unit cell.

Thus equilibration of kinetic and potential energies in graphene is accurately described by formula (6.4) for any distribution of the initial kinetic energy among degrees of freedom of the unit cell.

8. Redistribution of kinetic energy among degrees of freedom

In the present section, we consider the behaviour of kinetic energies K₁₁, K₂₂, K₃₃, K₄₄, corresponding to degrees of freedom of the unit cell (see formula (6.2)).

We investigate two processes. Firstly, we consider redistribution of the kinetic energy among two sublattices. Initial conditions are such that one sublattice is motionless, i.e. $K_{11}^{0} = K_{22}^{0} = \frac{1}{2} K_{0}$ , K⁰₃₃ = K⁰₄₄ = 0. Here, K₀ is the total initial kinetic energy of the unit cell. We calculate equilibrium kinetic energies by numerical evaluation of the integral in formula (6.5). Calculation yields the following values:

K_{11}^{eq} = K_{22}^{eq} = K_{33}^{eq} = K_{44}^{eq} = 0.125 K_{0} .

8.1

Remaining elements of matrix K_eq are negligible. Dependencies of kinetic energies on time, obtained by numerical solution of equations of motion and by numerical evaluation of integral (6.3) are presented in figure 7. The figure shows that kinetic energies of sublattices oscillate in time and tend to equal equilibrium values (8.1), predicted by the non-equipartition theorem (6.5). Therefore, in this case, energy is equally distributed between degrees of freedom of the unit cell. Note that kinetic energies K₁₁, K₂₂ and K₃₃, K₄₄ remain equal at any moment in time.

Figure 7. Redistribution of kinetic energy among sublattices in graphene ( $K_{11}^{0} = K_{22}^{0} = \frac{1}{2} K_{0}, K_{33}^{0} = K_{44}^{0} = 0$ ). Subplots show kinetic energies K₁₁ =K₂₂ (a) and K₃₃ =K₄₄ (b). Analytical solution (6.3) (line), simulation results (dots) and the equilibrium value (8.1) (dashed line). (Online version in colour.)

Consider redistribution of the kinetic energy among spatial directions. Initial conditions are such that particle velocities are directed along the x axis. In this case, $K_{11}^{0} = K_{33}^{0} = \frac{1}{2} K_{0}$ , K⁰₂₂ = 0, K⁰₄₄ = 0. We calculate equilibrium kinetic energies, K^eq_ii, by numerical evaluation of the integral in formula (6.5):

K_{11}^{eq} = K_{33}^{eq} = 0.1819 K_{0} and K_{22}^{eq} = K_{44}^{eq} = 0.0681 K_{0} .

8.2

Note that in this case some off-diagonal elements of matrix K_eq are not equal to zero, notably:

K_{13}^{eq} = K_{31}^{eq} = - 0.0073 K_{0} K_{24}^{eq} = K_{42}^{eq} = - 0.0026 K_{0} .

8.3

Therefore, velocities of two particles in a unit cell are correlated. Time dependencies of kinetic energies, obtained by numerical solution of equations of motion and by numerical evaluation of integral (6.3) are presented in figure 8. The figure shows that analytical and numerical results practically coincide.

Figure 8. Redistribution of kinetic energy among degrees of freedom for two sublattices ( $K_{11}^{0} = K_{33}^{0} = \frac{1}{2} K_{0}, K_{22}^{0} = K_{44}^{0} = 0$ ). Subplots show kinetic energies K₁₁ = K₃₃(a) and K₂₂ =K₄₄ (b). Analytical solution (6.3) (line), simulation results (circles) and equilibrium values (8.2) (dashed line). (Online version in colour.)

Note that convergence to equilibrium takes much longer than in the previous example. Figure 9 shows the behaviour of the energies on larger time intervals. It is seen that at large times kinetic energies oscillate with low frequency and tend to equilibrium values predicted by the non-equipartition theorem.

Thus kinetic energy is partially redistributed among degrees of freedom of the unit cell. Equilibrium values of the kinetic energies are generally different. Though in some particular cases they are equal (see e.g. formula (8.1)). Note that off-diagonal elements of matrix K, in general, are not equal to zero, i.e. velocities of particles inside a unit cell are correlated (see formula (8.3)).

9. Conclusion

Time evolution of kinetic energies, corresponding to different degrees of freedom of the unit cell, in graphene lattice with random initial velocities was investigated. It was shown that the evolution is exactly described by formula (6.3). The formula shows that the total kinetic and potential energies of the system oscillate in time and tend to equal equilibrium values. Characteristic frequency of these oscillations is higher than $\sqrt{c / M}$ . The oscillations are independent of initial distribution of kinetic energy among degrees of freedom.

If initial kinetic energies, corresponding to degrees of freedom of the unit cell, are different then kinetic energy is partially redistributed. At large times, the energies tend to generally different equilibrium values. Equipartition of energies was observed only in some particular cases. In general, the equilibrium values of the kinetic energies are accurately predicted by the non-equipartition theorem (6.5).

Redistribution of energy among degrees of freedom of the unit cell in harmonic crystals is usually a fast, high-frequency process [29,37,49]. The energies usually oscillate in time and practically reach equilibrium values at times of order of ten periods of atomic vibrations. Characteristic frequency of these oscillations is of order of $\sqrt{c / M}$ . It was shown that the oscillations in graphene contain very low frequency of order of $10^{- 2} \sqrt{c / M}$ .

In the present paper, spatially uniform distribution of kinetic energy (or kinetic temperature) was considered. In the case of non-uniform initial temperature profile, an additional physical process, notably the ballistic heat transport, should be considered. The heat transport is much slower than equilibration of energies considered above [30,50–54]. Therefore, at short times behaviour of kinetic energy at each spatial point can be approximately described by formula (6.5). It was shown that in a spatially uniform case, if initial kinetic energies, corresponding to degrees of freedom of the unit cell, are equal then their equilibrium values are also equal. Note that this fact does not guarantee that kinetic energies will remain equal during heat transport. In particular, in paper [11], it was shown that in heat-conducting harmonic crystals the kinetic energies are generally different even if their initial values are equal.

In the present paper, anharmonic effects were neglected. In paper [29], it is shown numerically that anharmonicity of interatomic interactions causes an additional relaxation process. This process leads, in particular, to equipartition of kinetic energies, corresponding to different degrees of freedom. However, at least in the case of small anharmonicity, the relaxation process is much slower than processes described above. Therefore, we believe that at short times harmonic approximation has reasonable accuracy. Investigation of anharmonic effects would be an interesting extension of the present work.

Thus presented results may serve for better understanding of thermal equilibration and ballistic heat transport in graphene and other two-dimensional materials.

Data accessibility

This article does not contain any additional data.

Author's contributions

Both authors contributed equally to this study.

Competing interests

We declare we have no competing interests.

Funding

Work of V.A.K. was financially supported by the Russian Science Foundation under grant no. 18-11-00201.

Acknowledgments

The authors are deeply grateful to A. M. Krivtsov and S. N. Gavrilov for useful discussions.

Footnotes

1 Exchange between energies corresponding to motion of particles in different spacial directions is considered. There is no energy exchange between the normal modes.

One contribution of 12 to a theme issue ‘Modelling of dynamic phenomena and localization in structured media (part 2)’.

References

1.
Bhimanapati GRet al.2015Recent advances in two-dimensional materials beyond graphene. ACS Nano 9, 11509–11539. (10.1021/acsnano.5b05556) Crossref, PubMed, ISI, Google Scholar
2.
Xu M, Liang T, Shi M, Chen H. 2013Graphene-like two-dimensional materials. Chem. Rev. 113, 3766–3798. (10.1021/cr300263a) Crossref, PubMed, ISI, Google Scholar
3.
Allen MJ, Tung VC, Kaner RB. 2009Honeycomb carbon: a review of graphene. Chem. Rev. 110, 132–145. (10.1021/cr900070d) Crossref, ISI, Google Scholar
4.
Novoselov KS, Fal VI, Colombo L, Gellert PR, Schwab MG, Kim K. 2012A roadmap for graphene. Nature 490, 192–200. (10.1038/nature11458) Crossref, PubMed, ISI, Google Scholar
5.
Bhuyan MSA, Uddin MN, Islam MM, Bipasha FA, Hossain SS. 2016Synthesis of graphene. Int. Nano Lett. 6, 65–83. (10.1007/s40089-015-0176-1) Crossref, ISI, Google Scholar
6.
Akinwande Det al.2017A review on mechanics and mechanical properties of 2D materials—graphene and beyond. Extrem. Mech. Lett. 13, 42–77. (10.1016/j.eml.2017.01.008) Crossref, ISI, Google Scholar
7.
Costas G, Frank O, Koukaras EN, Sfyris D. 2015Graphene mechanics: current status and perspectives. Annu. Rev. Chem. Biomol. Eng. 6, 121–140. (10.1146/annurev-chembioeng-061114-123216) Crossref, PubMed, ISI, Google Scholar
8.
Liu F, Ming PM, Li J. 2007Ab initio calculation of ideal strength and phonon instability of graphene under tension. Phys. Rev. B 76, 064120. (10.1103/PhysRevB.76.064120) Crossref, ISI, Google Scholar
9.
Balandin AA. 2011Thermal properties of graphene and nanostructured carbon materials. Nature Mater. 10, 569–581. (10.1038/nmat3064) Crossref, PubMed, ISI, Google Scholar
10.
Berinskii IE, Borodich FM. 2013On the isotropic elastic properties of graphene crystal lattice. In Surface Effects in Solid Mechanics pp. 33–42. Berlin, Heidelberg: Springer. Google Scholar
11.
Kuzkin VA. 2019Unsteady ballistic heat transport in harmonic crystals with polyatomic unit cell. Cont. Mech. Thermodyn.1–27. (10.1007/s00161-019-00802-1) ISI, Google Scholar
12.
Renteria JD, Ramirez S, Malekpour H, Alonso B, Centeno A, Zurutuza A, Cocemasov AI, Nika DL, Balandin AA. 2015Strongly anisotropic thermal conductivity of free–standing reduced graphene oxide films annealed at high temperature. Adv. Funct. Mater. 25, 4664–4672. (10.1002/adfm.201501429) Crossref, ISI, Google Scholar
13.
Cahill DG, Ford WK, Goodson KE, Mahan GD, Majumdar A, Maris HJ, Merlin R, Phillpot SR. 2003Nanoscale thermal transport. J. Appl. Phys. 93, 793. (10.1063/1.1524305) Crossref, ISI, Google Scholar
14.
Chang C-W. 2016Thermal transport in low dimensions. Lect. Notes Phys. 921, 305–338. Crossref, Google Scholar
15.
Hsiao TK, Chang HK, Liou S-C, Chu M-W, Lee S-C, Chang C-W. 2013Observation of room-temperature ballistic thermal conduction persisting over 8.3 μm SiGe nanowires. Nat. Nanotech. 8, 534–538. (10.1038/nnano.2013.121) Crossref, PubMed, ISI, Google Scholar
16.
Pumarol ME, Rosamond MC, Tovee P, Petty MC, Zeze DA, Falko V, Kolosov OV. 2012Direct nanoscale imaging of ballistic and diffusive thermal transport in graphene nanostructures. Nano Lett. 12, 2906. (10.1021/nl3004946) Crossref, PubMed, ISI, Google Scholar
17.
Huberman S, Duncan RA, Chen K, Song B, Chiloyan V, Ding Z, Maznev AA, Chen G, Nelson KA. 2019Observation of second sound in graphite at temperatures above 100 K. Science 364, 375–379. (10.1126/science.aav3548) Crossref, PubMed, ISI, Google Scholar
18.
Indeitsev DA, Osipova EV. 2017A two-temperature model of optical excitation of acousticwaves in conductors. Dokl. Phys. 62, 136–140. (10.1134/S1028335817030065) Crossref, ISI, Google Scholar
19.
Johnson JA, Maznev AA, Cuffe J, Eliason JK, Minnich AJ, Kehoe T, Sotomayor Torres CM, Chen G, Nelson KA. 2013Direct measurement of room-temperature nondiffusive thermal transport over micron distances in a silicon membrane. Phys. Rev. Lett. 110, 025901. (10.1103/PhysRevLett.110.025901) Crossref, PubMed, ISI, Google Scholar
20.
Poletkin KV, Gurzadyan GG, Shang J, Kulish V. 2012Ultrafast heat transfer on nanoscale in thin gold films. Appl. Phys. B 107, 137–143. (10.1007/s00340-011-4862-z) Crossref, ISI, Google Scholar
21.
Rogers JA, Maznev AA, Banet MJ, Nelson KA. 2000Optical generation and characterization of acousticwaves in thin films: fundamentals and applications. Annu. Rev. Mater. Sci. 30, 117–157. (10.1146/annurev.matsci.30.1.117) Crossref, ISI, Google Scholar
22.
Boldrighini C, Pellegrinotti A, Triolo L. 1983Convergence to stationary states for infinite harmonic systems. J. Stat. Phys. 30, 123–155. (10.1007/BF01010871) Crossref, ISI, Google Scholar
23.
Dobrushin RL, Pellegrinotti A, Suhov Yu.M, Triolo L. 1986One-dimensional harmonic lattice caricature of hydrodynamics. J. Stat. Phys. 43, 571–607. (10.1007/BF01020654) Crossref, ISI, Google Scholar
24.
Dudnikova TV, Komech AI, Spohn H. 2003On the convergence to statistical equilibrium for harmonic crystals. J. Math. Phys. 44, 2596. (10.1063/1.1571658) Crossref, ISI, Google Scholar
25.
Hemmer PC. 1959Dynamic and stochastic types of motion in the linear chain. Norges tekniske hoiskole. Google Scholar
26.
Klein G, Prigogine I. 1953Sur la mecanique statistique des phenomenes irreversibles III. Physica 19, 1053. (10.1016/S0031-8914(53)80120-5) Crossref, Google Scholar
27.
Krivtsov AM. 2014Energy oscillations in a one-dimensional crystal. Dokl. Phys. 59, 427–430. (10.1134/S1028335814090080) Crossref, ISI, Google Scholar
28.
Kuzkin VA, Krivtsov AM. 2017High-frequency thermal processes in harmonic crystals. Dokl. Phys. 62, 85–89. (10.1134/S1028335817020070) Crossref, ISI, Google Scholar
29.
Kuzkin VA, Krivtsov AM. 2017An analytical description of transient thermal processes in harmonic crystals. Phys. Solid State 59, 1051. (10.1134/S1063783417050201) Crossref, ISI, Google Scholar
30.
Kuzkin VA, Krivtsov AM. 2017Fast and slow thermal processes in harmonic scalar lattices. J. Phys.: Condens. Matter 29, 505401. (10.1088/1361-648X/aa98eb) Crossref, PubMed, ISI, Google Scholar
31.
Lanford OE, Lebowitz JL. 1975Time evolution and ergodic properties of harmonic systems. In Lecture Notes in Physics, vol. 38, pp. 144–177. Berlin, Heidelberg; New York, NY: Springer. Google Scholar
32.
Spohn H, Lebowitz JL. 1977Stationary non-equilibrium states of infinite harmonic systems. Commun. Math. Phys. 54, 97–120. (10.1007/BF01614132) Crossref, ISI, Google Scholar
33.
Gavrilov SN, Krivtsov AM, Tsvetkov DV. 2019Heat transfer in a one-dimensional harmonic crystal in a viscous environment subjected to an external heat supply. Cont. Mech. Thermodyn. 31, 255–272. (10.1007/s00161-018-0681-3) Crossref, ISI, Google Scholar
34.
Gavrilov SN, Krivtsov AM. 2019Steady-state kinetic temperature distribution in a two-dimensional square harmonic scalar lattice lying in a viscous environment and subjected to a point heat source. Cont. Mech. Thermodyn.1–21. (10.1007/s00161-019-00782-2) ISI, Google Scholar
35.
Babenkov MB, Krivtsov AM, Tsvetkov DV. 2019Energy oscillations in 1D harmonic crystal on elastic foundation. Phys. Mesomech. 19, 60–67. Google Scholar
36.
Murachev AS, Krivtsov AM, Tsvetkov DV. 2019Thermal echo in a finite one-dimensional harmonic crystal. J. Phys. Condensed Matter 31, 095702. (10.1088/1361-648X/aaf3c6) Crossref, PubMed, ISI, Google Scholar
37.
Kuzkin VA. 2019Thermal equilibration in infinite harmonic crytals. Cont. Mech. Thermodyn. 31, 1401–1423. (10.1007/s00161-019-00758-2) Crossref, ISI, Google Scholar
38.
Berinskii IE, Borodich FM. 2013Elastic in-plane properties of 2D linearized models of graphene. Mech. Mater. 62, 60–68. (10.1016/j.mechmat.2013.03.004) Crossref, ISI, Google Scholar
39.
Berinskii IE. 2016Elastic networks to model auxetic properties of cellular materials. Int. J. Mech. Sci. 115, 481–488. (10.1016/j.ijmecsci.2016.07.038) Crossref, ISI, Google Scholar
40.
Barani E, Lobzenko IP, Korznikova EA, Soboleva EG, Dmitriev SV, Zhou K, Marjaneh AM. 2017Transverse discrete breathers in unstrained graphene. Eur. Phys. J. B 90, 1. (10.1140/epjb/e2017-70751-2) Crossref, ISI, Google Scholar
41.
Kuzkin VA. 2010Interatomic force in systems with multibody interactions. Phys. Rev. E 82, 016704. (10.1103/PhysRevE.82.016704) Crossref, ISI, Google Scholar
42.
Berinskiy IE, Krivtsov AM, Kudarova AM. 2008Determination of macroscopic characteristics for graphene layer using angle-depending atomic interactions. In Proc. of XXXVI Int. Summer School-Conference APM, p. 122. Google Scholar
43.
Kittel C. 1996Introduction to solid state physics. New York, NY: John Wiley and Sons Inc. Google Scholar
44.
Ashcroft NW, Merin ND. 1976Solid state physics Brooks. Belmont, CA: Cole. Google Scholar
45.
Mishuris GS, Movchan AB, Bercial JP. 2007Asymptotic analysis of Bloch-Floquet waves in a thin bi-material strip with a periodic array of finite-length cracks. Waves Random Complex Media 17, 511–533. (10.1080/17455030701288137) Crossref, ISI, Google Scholar
46.
Mishuris GS, Movchan AB, Slepyan LI. 2009Localised knife waves in a structured interface. J. Mech. Phys. Solids 57, 1958. (10.1016/j.jmps.2009.08.004) Crossref, ISI, Google Scholar
47.
Dove MT. 1993Introduction to lattice dynamics. London, UK: Cambridge University Press. Crossref, Google Scholar
48.
Slepyan LI. 2015On the energy partition in oscillations and waves. Proc. R. Soc. A 471, 20140838. (10.1098/rspa.2014.0838) Link, Google Scholar
49.
Tsaplin VA, Kuzkin VA. 2018Temperature oscillations in harmonic triangular lattice with random initial velocities. Lett. Mater. 8, 16–20. (10.22226/2410-3535-2018-1-16-20) Crossref, ISI, Google Scholar
50.
Gendelman OV, Savin AV. 2010Nonstationary heat conduction in one-dimensional chains with conserved momentum. Phys. Rev. E 81, 020103. (10.1103/PhysRevE.81.020103) Crossref, ISI, Google Scholar
51.
Gendelman OV, Shvartsman R, Madar B, Savin AV. 2012Nonstationary heat conduction in one-dimensional models with substrate potential. Phys. Rev. E 85, 011105. (10.1103/PhysRevE.85.011105) Crossref, ISI, Google Scholar
52.
Krivtsov AM. 2015Heat transfer in infinite harmonic one dimensional crystals. Dokl. Phys. 60, 407–411. (10.1134/S1028335815090062) Crossref, ISI, Google Scholar
53.
Krivtsov AM. 2019The ballistic heat equation for a one-dimensional harmonic crystal. In Dynamical processes in generalized continua and structures, Springer Nature. Google Scholar
54.
Sokolov AA, Krivtsov AM, Muller WH, Vilchevskaya EN. 2019Change of entropy for the one-dimensional ballistic heat equation: sinusoidal initial perturbation. Phys. Rev. E 99, 042107. (10.1103/PhysRevE.99.042107) Crossref, PubMed, ISI, Google Scholar

This Issue

Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences cover image

10 January 2020

Volume 378Issue 2162

Theme issue 'Modelling of dynamic phenomena and localization in structured media (part 2)' compiled and edited by G. Mishuris, A. Movchan and L. Truskinovsky

Article Information

DOI:https://doi.org/10.1098/rsta.2019.0114
PubMed:31760904
Published by:Royal Society
Print ISSN:1364-503X
Online ISSN:1471-2962

History:

Manuscript accepted09/07/2019
Published online25/11/2019
Published in print10/01/2020

License:

Citations and impact

Keywords

Subjects