Virial theorem, boundary conditions, and pressure for massless Dirac electrons

TL;DR

This paper derives and analyzes the virial and Hellmann--Feynman theorems for massless Dirac electrons, revealing how boundary conditions affect the Hamiltonian's Hermiticity and the resulting pressure components in solid-state systems.

Contribution

It introduces generalized theorems for massless Dirac electrons, accounting for boundary effects and boundary-induced pressures, with applications to graphene and comparison to massive electrons.

Findings

Boundary conditions break Hamiltonian Hermiticity, leading to boundary integral terms.

Electron reflections contribute to the kinetic pressure component.

Distinct anomalous pressure arises in solid-state Dirac electron systems.

Abstract

The virial and the Hellmann--Feynman theorems for massless Dirac electrons in a solid are derived and analyzed using generalized continuity equations and scaling transformations. Boundary conditions imposed on the wave function in a finite sample are shown to break the Hermiticity of the Hamiltonian resulting in additional terms in the theorems in the forms of boundary integrals. The thermodynamic pressure of the electron gas is shown to be composed of the kinetic pressure, which is related to the boundary integral in the virial theorem and arises due to electron reflections from the boundary, and the anomalous pressure, which is specific for electrons in solids. Connections between the kinetic pressure and the properties of the wave function on the boundary are drawn. The general theorems are illustrated by examples of uniform electron gas, and electrons in rectangular and circular graphene samples. The analogous consideration for ordinary massive electrons is presented for comparison.

Equations

$$-{\psi }^{\ast }[H,A]\psi +{\psi }^{\ast }(H-H^{+})A\psi =0.$$

$${\psi }^{\ast }(H-H^{+})A\psi =-i\hslash \mathrm{div}\mathbf{J}[A],$$

$$\mathbf{J}[A]=\frac{\hslash }{2mi}{\psi }^{\ast }(\nabla -{\nabla }^{+})A\psi$$

$$-\frac{i}{\hslash }⟨\psi |[H,A]|\psi ⟩+\underset{\partial \Omega }{\oint }d\mathbf{s}\cdot \mathbf{J}[A]=0,$$

$$G=\frac{\mathbf{r}\cdot \mathbf{p}+\mathbf{p}\cdot \mathbf{r}}{2}=-i\hslash \left(\mathbf{r}\cdot \nabla +\frac{D}{2}\right),$$

$$⟨\psi |-2H_{\mathrm{kin}}+\mathbf{r}\cdot \nabla U_{\mathrm{ext}}|\psi ⟩+\underset{\partial \Omega }{\oint }d\mathbf{s}\cdot \mathbf{J}[G]=0.$$

$$\delta E=⟨\psi |\delta H|\psi ⟩+⟨\psi |H-H^{+}|\delta \psi ⟩.$$

$$\frac{\partial E}{\partial \lambda }=⟨\psi \left|\frac{\partial H}{\partial \lambda }\right|\psi ⟩+⟨\psi \left|H-H^{+}\right|\frac{\partial \psi }{\partial \lambda }⟩.$$

$$\frac{\partial E}{\partial \lambda }=⟨\psi \left|\frac{\partial H}{\partial \lambda }\right|\psi ⟩-i\hslash \underset{\partial \Omega }{\oint }d\mathbf{s}\cdot \mathbf{J}\left[\frac{\partial }{\partial \lambda }\right].$$

$$D\mathcal{P}\Omega =-R\frac{\partial E}{\partial R}=i\hslash \underset{\partial \Omega }{\oint }d\mathbf{s}\cdot \mathbf{J}\left[R\frac{\partial }{\partial R}\right].$$

$$⟨\psi |-2H_{\mathrm{kin}}+\mathbf{r}\cdot \nabla U_{\mathrm{ext}}|\psi ⟩+D\mathcal{P}\Omega =0,$$

$$D\mathcal{P}\Omega =⟨\psi |\mathbf{r}\cdot \nabla U_{\mathrm{b}}|\psi ⟩.$$

$$\delta \psi ({\mathbf{r}}_0)=-\delta {\mathbf{r}}_0\cdot \nabla \psi ({\mathbf{r}}_0)+\mathcal{O}([\delta {\mathbf{r}}_0{]}^2).$$

$$R\frac{\partial \psi }{\partial R}=-\mathbf{r}\cdot \nabla \psi .$$

$$\mathbf{J}[G]=i\hslash \mathbf{J}\left[\frac{\partial }{\partial R}\right].$$

$$⟨\psi |-2H_{\mathrm{kin}}+\mathbf{r}\cdot \nabla U_{\mathrm{ext}}|\psi ⟩+D\mathcal{P}\Omega =0,$$

$$D\mathcal{P}\Omega =\underset{\partial \Omega }{\oint }d\mathbf{s}\cdot \mathbf{J}[G]=\frac{{\hslash }^2}{2m}\underset{\partial \Omega }{\oint }d\mathbf{s}\cdot \mathbf{r}|\nabla \psi {|}^2$$

$$\delta E=-\frac{{\hslash }^2}{2m}\underset{\partial \Omega }{\oint }d\mathbf{s}\cdot {\psi }^{\ast }{\nabla }^{+}\delta \psi .$$

$$\delta E=-\frac{{\hslash }^2}{2m}\underset{\partial \Omega }{\oint }d{s}_{\nu }(\delta {\mathbf{r}}_0{)}_{\mu }{\psi }^{\ast }{\nabla }_{\nu }^{+}{\nabla }_{\mu }\psi .$$

$$T_{\mu \nu }\equiv J_{\nu }[p_{\mu }]=\frac{{\hslash }^2}{2m}{\psi }^{\ast }{\nabla }_{\mu }({\nabla }_{\nu }^{+}-{\nabla }_{\nu })\psi ,$$

$$\delta E=-\underset{\partial \Omega }{\oint }d{s}_{\nu }(\delta {\mathbf{r}}_0{)}_{\mu }{T}_{\mu \nu }.$$

$$\delta E=-\underset{\partial \Omega }{\oint }ds\delta {\mathbf{r}}_0\cdot \mathbf{P},$$

$$P_{\mu }=T_{\mu \nu }{n}_{\nu },$$

$$D\mathcal{P}V=\underset{\partial \Omega }{\oint }ds\mathbf{r}\cdot \mathbf{P}=\underset{\partial \Omega }{\oint }d{s}_{\nu }{r}_{\mu }{T}_{\mu \nu }.$$

$${\rho }_N=\sum _n{w}_n|{\Psi }_n⟩⟨{\Psi }_n|,$$

$${\rho }_1(\mathbf{r},{\mathbf{r}}^{\prime })=\sum _n{w}_n\int \prod _{j}d{\mathbf{r}}_{j}d{\mathbf{r}}_j^{\prime }\sum _i\delta (\mathbf{r}-{\mathbf{r}}_i)$$

$$\times \delta ({\mathbf{r}}^{\prime }-{\mathbf{r}}_i^{\prime }){\Psi }_n({\mathbf{r}}_1\dots {\mathbf{r}}_N){\Psi }_n^{\ast }({\mathbf{r}}_1^{\prime }\dots {\mathbf{r}}_N^{\prime }),$$

$$\mathbf{J}[A]=\frac{\hslash }{2mi}{({\nabla }_{\mathbf{r}}-{\nabla }_{{\mathbf{r}}^{\prime }})A_{\mathbf{r}}{\rho }_1(\mathbf{r},{\mathbf{r}}^{\prime })|}_{{\mathbf{r}}^{\prime }=\mathbf{r}}.$$

$$⟨-2H_{\mathrm{kin}}-V_{\mathrm{int}}+\mathbf{r}\cdot \nabla U_{\mathrm{ext}}⟩+\underset{\partial \Omega }{\oint }d\mathbf{s}\cdot \mathbf{J}[G]=0$$

$$\frac{\partial E}{\partial \lambda }=\sum _n⟨{\Psi }_n\left|\frac{\partial (w_{n}H)}{\partial \lambda }\right|{\Psi }_n⟩-i\hslash \underset{\partial \Omega }{\oint }d\mathbf{s}\cdot \mathbf{J}\left[\frac{\partial }{\partial \lambda }\right].$$

$$D\mathcal{P}\Omega =i\hslash \underset{\partial \Omega }{\oint }d\mathbf{s}\cdot \mathbf{J}\left[R\frac{\partial }{\partial R}\right],$$

$$⟨-2H_{\mathrm{kin}}-V_{\mathrm{int}}+\mathbf{r}\cdot \nabla U_{\mathrm{ext}}⟩+D\mathcal{P}\Omega =0,$$

Full text

I Introduction

Discovery of graphene CastroNeto and three-dimensional Dirac and Weyl semimetals Armitage , where the electron low-energy dynamics is described by the effective Dirac equation for massless particles, uncovered the new area of solid state physics of Dirac materials Wehling . Peculiar properties of these materials motivate researchers to reconsider the conventional notions and models of quantum electron phenomena initially developed for massive electrons. One of such notions is the virial theorem and related quantum theory of stress and pressure of electron gas Marc ; Nielsen ; Maranganti ; MartinPendas .

The virial theorem for a system of interacting particles provides the relationship between average kinetic energy, Coulomb interaction energy, and external pressure Marc . Applications of the virial theorem in classical and quantum statistical physics include estimation of the system properties, derivation of equations of state, checking accuracy of quantum chemistry and density functional calculations etc. The quantum-mechanical virial theorem for a system of ordinary massive electrons can be derived by using scaling transformations of an electron wave function Fock ; Lowdin and spatial integration of the continuity equation for momentum density Ziesche ; Godfrey . The pressure can be exerted on the system by Coulomb potentials of atomic nuclei Ziesche ; Bader , by a generic external potential, or by impenetrable walls with the Dirichlet boundary conditions imposed on the electron wave function. In the last case, the Hamiltonian of the system is Hermitian only in the subspace of wave functions satisfying these conditions, and the scaling transformations drive the functions out of the Hermiticity domain Abad . The non-Hermiticity of the Hamiltonian in the presence of the scaling transformations results in emergence of the additional term in the virial theorem, which is proportional to pressure and has a form of an integral of the virial current density through the boundary Marc ; Argyres ; Bobrov ; Srebrenik . An alternative way to calculate pressure as a response of the system energy to small volume changes relies on the Hellmann–Feynman theorem. The non-Hermiticity of a Hamiltonian of enclosed system in the space of perturbed wave functions leads to emergence of the boundary integral in the Hellmann–Feynman theorem as well Esteve ; Argyres ; Lowdin ; Konstantinou ; Srebrenik . The consistency between these two definitions of pressure (see Ref. Marc , p. 289) requires fulfilment of a specific boundary relationship for the wave function on the system boundary which relates its spatial derivative to its derivative with respect to the boundary location Cottrell ; Fernandez ; Argyres . Note that the virial theorem is a particular case of more general stress theorem Nielsen ; Maranganti ; Godfrey ; Bader .

In the Dirac materials, the linear electron dispersion causes modification of the conventional virial theorem Stokes ; Sokolik . Moreover, the momentum cutoff imposed at the bottom of the valence band in order to bound the system energy from below leads to appearance of the additional term in the resulting generalized virial theorem Sokolik . However, the proper quantum-mechanical analysis of boundary contributions to the virial theorem for massless Dirac electrons is still lacking. The electron wave function in these materials is multi-component and obeys the boundary conditions which differ from the Dirichlet ones for massive electrons. For graphene, the infinite mass Berry ; McCann , zigzag, and armchair Brey ; Akhmerov ; McCann boundary conditions are used depending on the lattice edge crystal structure. For three-dimensional Dirac and Weyl semimetals, various boundary conditions are proposed Volkov ; Hashimoto . Other possible anomalies in scaling properties of a system of massless Dirac electrons can also give rise to additional terms in the virial theorem Lin .

In this paper, we derive the generalized virial and Hellmann–Feynman theorems for massless Dirac electrons, that contain additional terms coming from the non-Hermiticity of the Hamiltonian in the presence of a system boundary and from the momentum cutoff. Associating these terms with the pressure, we show that the thermodynamic pressure of Dirac electrons is the sum of the kinetic pressure, which is caused by reflections of electrons from the boundary, and the anomalous pressure, which is caused by redistribution of electron states during changes of the system area in the presence of the momentum cutoff.

For the kinetic pressure, we show that the physically relevant boundary conditions imply the boundary relationship for the wave function, which allows us to achieve consistency between definitions of the pressure based on the virial and Hellmann–Feynman theorems. In addition to the total pressure, we calculate the local pressure on the boundary and connect it with the kinetic part of the stress tensor. To illustrate the derived theorems, we consider several particular examples: uniform electron gas in graphene, and rectangular and circular graphene flakes with the appropriate boundary conditions. The calculations of the electron pressure in these examples are in agreement with the generalized virial theorem.

The article is organized as follows. In Section II we consider the virial and Hellmann–Feynman theorems for ordinary massive electrons, and derive the boundary relationship and the quantum-mechanical expressions for the electron pressure. In Section III we provide the similar analysis for massless Dirac electrons and reveal the important differences stemming from the momentum cutoff and from different form of boundary conditions. In Section IV we consider electrons in graphene samples of different geometries and demonstrate fulfilment of the general theorems in these systems, and in Section V we summarize and discuss our results. A and B are devoted to consideration of the scaling properties of uniform gases of massive and massless electrons, which are closely connected with the virial theorems. In C we consider two-band model of massive electrons and in D we calculate corrections to graphene electron gas properties caused by electron dispersion nonlinearities.

II Massive electrons

II.1 Virial theorem

The virial and Hellmann–Feynman theorems for a bounded system of ordinary massive electrons Marc can be derived from the generalized continuity equations. We will assume that the system is surrounded by infinitely high potential walls, so the wave function obeys the Dirichlet boundary conditions. Consider first a single particle with a stationary state wave function $\psi(\mathbf{r})$ obeying the Schrodinger equation $H\psi=E\psi$ locally at the point $\mathbf{r}$. Multiplying this equation by $\psi^{*}A$ from the left (where $A$ is some operator) and subtracting the Hermitian conjugate equation $\psi^{*}H^{+}=E\psi^{*}$, which is multiplied by $A\psi$ from the right, we get

[TABLE]

Since the wave function $\psi$ can disobey the Dirichlet boundary condition after action of $A$, the Hamiltonian bracketed between $\psi^{*}$ and $A\psi$ becomes, in general case, non-Hermitian. In coordinate representation this non-Hermiticity is demonstrated only by the kinetic part $H_{\mathrm{kin}}=-\hbar^{2}\nabla^{2}/2m$ of $H$, and we can write

[TABLE]

where

[TABLE]

is the single-particle generalized current density of the quantity corresponding to the operator $A$, e. g., the probability current if $A=1$ or the momentum current if $A=p_{\mu}$ Marc . Hereafter we treat $\nabla\equiv\overrightarrow{\nabla}$ and $\nabla^{+}\equiv\overleftarrow{\nabla}$ as the operators, which act on the functions, respectively, to the right and to the left. Substituting (2)–(3) in (1), integrating over the volume of the system $\Omega$ and applying the Gauss theorem, we get the formula

[TABLE]

which equates total generation rate of the quantity $A$ in the system to the flux of this quantity out of the system in a stationary state.

As $A$, we can take the virial (virial of momentum, to be more precise) operator

[TABLE]

where $D$ is the space dimensionality. If the particle moves in the external potential $U_{\mathrm{ext}}(\mathbf{r})$, then $H=H_{\mathrm{kin}}+U_{\mathrm{ext}}$ and $[H,G]=i\hbar(-2H_{\mathrm{kin}}+\mathbf{r}\cdot\nabla U_{\mathrm{ext}})$. With $A=G$, Eq. (4) takes the form of the virial theorem with the boundary term:

[TABLE]

II.2 Generalized Hellmann–Feynman theorem

Let us return to a single-particle Schrodinger equation $H|\psi\rangle=E|\psi\rangle$ written for the whole state vector $|\psi\rangle$ and admit a small variation of the Hamiltonian $\delta H$ and/or boundary conditions, resulting in a small variation $|\delta\psi\rangle$ of $|\psi\rangle$. Taking into account conservation of the wave function normalization $\langle\delta\psi|\psi\rangle+\langle\psi|\delta\psi\rangle=0$, which implies $\langle\delta\psi|H|\psi\rangle=-\langle\psi|H^{+}|\delta\psi\rangle$, we can write the variation of energy $E=\langle\psi|H|\psi\rangle$ as

[TABLE]

Assuming the variation $\delta\lambda$ of some parameter as a physical origin of both $\delta H$ and $|\delta\psi\rangle$, we get the Hellmann–Feynman theorem, generalized for the case of non-Hermitian Hamiltonian of a bounded system:

[TABLE]

Using (2) and the Gauss theorem, we can rewrite (8) as

[TABLE]

The boundary integrals in the virial (6) and Hellmann–Feynman (9) theorems for massive electrons Marc can be related to the pressure, as will be shown below.

II.3 Thermodynamic pressure

The thermodynamic pressure $\mathcal{P}=-\partial E/\partial\Omega$ is defined as a response of the system energy to adiabatically slow volume change. The latter can be introduced as a uniform and isotropic dilation or contraction of the system boundary points $\mathbf{r}\rightarrow\mathbf{r}(1+\delta R/R)$, where $R$ is a linear size of the system. In this approach the Hamiltonian of the system does not change, and the wave function is affected only by the change of boundary conditions. Taking $R$ as a slowly varying parameter in (9), we obtain

[TABLE]

An alternative way to calculate the pressure is to assume a finite-height confining potential $U_{\mathrm{b}}(\mathbf{r})$ on the boundary added to the Hamiltonian, so we can discard the boundary integrals, because the Hamiltonian becomes Hermitian due to vanishing of $\psi$ and $\nabla\psi$ at $|\mathbf{r}|\rightarrow\infty$. Introducing the system size dependence $U_{\mathrm{b}}(\mathbf{r})=\tilde{U}_{\mathrm{b}}(\mathbf{r}/R)$ and applying the ordinary virial (6) and Hellmann–Feynman theorems, we get:

[TABLE]

In this paper, we do not use this method to define the system boundary, because it is inapplicable in the case of massless Dirac electrons, which cannot be confined by a scalar potential CastroNeto because their energy spectrum is unbounded from below. Instead, for both kinds of electrons we define the boundary directly through the boundary conditions for the wave function, giving rise to the boundary integral in (6). Connection of this integral with the physically measurable pressure will be drawn below.

II.4 Boundary relationship

To calculate $\mathcal{P}$ using Eq. (10), we need to know the derivative $\partial\psi/\partial R$ on the system boundary. Suppose $\psi$ and $\tilde{\psi}$ are the wave functions of the same stationary state at, respectively, initial and slightly perturbed boundaries (see Fig. 1). If any boundary point $\mathbf{r}_{0}$ moves outwards on a small vector $\delta\mathbf{r}_{0}$, then the Dirichlet boundary conditions imply $\psi(\mathbf{r}_{0})=0$ and $\tilde{\psi}(\mathbf{r}_{0}+\delta\mathbf{r}_{0})=0$, and the wave function change in a fixed point is $\delta\psi(\mathbf{r}_{0})\equiv\tilde{\psi}(\mathbf{r}_{0})-\psi(\mathbf{r}_{0})=\tilde{\psi}(\mathbf{r}_{0})$. We can consider $\psi(\mathbf{r})=\psi(\mathbf{r};\mathbf{r}_{0})$ as a function of the vector $\mathbf{r}$ and the boundary position $\mathbf{r_{0}}$, so $\tilde{\psi}(\mathbf{r})=\psi(\mathbf{r};\mathbf{r}_{0}+\delta\mathbf{r}_{0})$, and $\psi(\mathbf{r};\mathbf{r}_{0})$ vanishes when its arguments coincide. Hence $\delta\psi(\mathbf{r}_{0})=\psi(\mathbf{r}_{0};\mathbf{r}_{0}+\delta\mathbf{r}_{0})=\delta\mathbf{r}_{0}\cdot\nabla_{\mathbf{r}_{0}}\psi(\mathbf{r};\mathbf{r}_{0})|_{\mathbf{r}=\mathbf{r}_{0}}+\mathcal{O}([\delta\mathbf{r}_{0}]^{2})$. By using the property $(\nabla_{\mathbf{r}}+\nabla_{\mathbf{r}_{0}})\psi(\mathbf{r};\mathbf{r}_{0})=0$, we obtain the boundary relationship

[TABLE]

We can see from Fig. 1 that $\delta\psi$ is indeed proportional to $\nabla\psi$ because the main cause of the wave function change is just a motion of the boundary and the wave function as a whole, while the change of $\nabla\psi$ provides only a second-order contribution.

Assuming a uniform and isotropic contraction or dilation of the system boundary $\mathbf{r}_{0}\rightarrow\mathbf{r}_{0}(1+\delta R/R)$, we get from (13) another version of the boundary relationship valid on the system boundary (see also Cottrell ; Fernandez ; Argyres ):

[TABLE]

In A we show how this relationship can be obtained based on the scaling arguments. Using (3), (5) (14) and the Dirichlet condition, we see that on the boundary

[TABLE]

This formula can be applied to unify (6) and (10) as the virial theorem with the pressure term:

[TABLE]

(here we have used that Dirichlet boundary conditions imply that $\nabla\psi$ is directed parallel to the normal to the boundary). It is similar to (11)–(12), but formulated for a bounded system.

II.5 Local pressure

Now let us calculate from (7) the response of the system energy on arbitrary small perturbations $\delta\mathbf{r}_{0}$ of the boundary points $\mathbf{r}_{0}$. Using (2)–(3) with $A\psi=\delta\psi$ and applying the Gauss theorem, we get

[TABLE]

The boundary relationship (13) allows us to rewrite it in the form

[TABLE]

Introducing the kinetic stress tensor (or spatial part of the stress-energy tensor, or momentum flux density) by applying (3) to the momentum operator Maranganti ,

[TABLE]

we recast (19) as

[TABLE]

Associating $\delta E$ with the work done by external forces acting on the system, which are opposite to the vector of local pressure $\mathbf{P}$ of the system itself, we write

[TABLE]

thus, given the arbitrariness of $\delta\mathbf{r}_{0}$, the comparison of (21) and (22) results in

[TABLE]

where $\mathbf{n}$ is the unit normal to the boundary and directed outside. The physical meaning of Eq. (23) is that the system pressure $P_{\mu}$ exerted to the surroundings at some point is caused by particle collisions with the boundary, which transfer momentum at the rate proportional to the normal component $\mathbf{J}[p_{\mu}]\cdot\mathbf{n}$ of the momentum flux at that point.

In the case of uniform dilation $\delta\mathbf{r}_{0}=(\delta R/R)\mathbf{r}_{0}$ of the system boundary we recover (17) with the relationship between thermodynamic and local pressures:

[TABLE]

The integrals here are independent of the choice of the origin because $\oint_{\partial\Omega}d\mathbf{s}\cdot\mathbf{P}=0$ for a system being in mechanical equilibrium with its surroundings.

II.6 Many-body system

The many-body and thermal ensemble generalizations of all calculations presented above are rather straightforward. Assume that the system state is characterized by an $N$-particle density matrix

[TABLE]

where $|\Psi_{n}\rangle$ are the eigenstates, $H|\Psi_{n}\rangle=E_{n}|\Psi_{n}\rangle$, of the many-body Hamiltonian $H=\sum_{i}[(-\hbar^{2}\nabla_{i}^{2}/2m)+U_{\mathrm{ext}}(\mathbf{r}_{i})]+(1/2)\sum_{i\neq j}V_{\mathrm{int}}(\mathbf{r}_{i}-\mathbf{r}_{j})$ with the energies $E_{n}$, entering the ensemble with the probabilities $w_{n}$. As specific examples, we can consider the many-body ground state $|\Psi_{0}\rangle$ at $T=0$, where $w_{n}=\delta_{n0}$, or the thermal state, where $w_{n}\propto\exp(-E_{n}/T)$.

Introducing the one-body density matrix

[TABLE]

we can define the many-body counterpart of the generalized current (3):

[TABLE]

Another distinction of the many-body system is the presence of the interparticle interaction $V_{\mathrm{int}}$, leading to an additional term $\propto[H,V_{\mathrm{int}}]$ in the virial theorem. For the Coulomb interaction, we have $[H,V_{\mathrm{int}}]=-i\hbar V_{\mathrm{int}}$, and the many-body virial theorem can be obtained by taking the linear combination of the single-particle ones (6) with the coefficients $w_{n}$:

[TABLE]

(here $\langle A\rangle$ is defined as $\mathrm{Tr}\,[\rho A]$). The generalized Hellmann–Feynman theorem (8) for a many-body system is

[TABLE]

The thermodynamics pressure $\mathcal{P}$ can be calculated using (29) with $\lambda=R$. At $T=0$ we have $\mathcal{P}=-\partial E_{0}/\partial\Omega$, but at $T>0$ we need to consider the free energy $F$ to define the pressure $\mathcal{P}=-\partial F/\partial\Omega$. In both cases the derivatives of $w_{n}$ in (29) do not appear in the resulting formula for $\mathcal{P}$, and we obtain

[TABLE]

which looks equivalent to (10), although with the many-body current operator (27).

The Dirichlet boundary conditions imposed on the many-body wave function $\Psi_{n}(\mathbf{r}_{1}\ldots\mathbf{r}_{N})$ imply $\Psi_{n}=0$ when any of $\mathbf{r}_{i}$ is located on the boundary. Therefore the same boundary relationship (13) is valid for $\Psi_{n}$ when $\mathbf{r}$ is replaced by any of its arguments $\mathbf{r}_{i}$, and Eq. (15) is valid for a many-body system as well. Combining (15), (27), (28), and (30), we obtain the many-body counterpart of the single-particle virial theorem (16)–(17) with the pressure term:

[TABLE]

Derivation of (31) in the case of power-law $U_{\mathrm{ext}}$ by using scaling properties is shown in A. The consideration (21)–(23) of the local pressure can be also repeated with the many-body kinetic stress tensor

[TABLE]

Alongside with the local pressure on the boundary (23), we can define the local kinetic pressure MartinPendas ; Ziesche ; Godfrey

[TABLE]

in the bulk related to a trace of the total stress tensor consisting of the kinetic part $T$, interaction part $\sigma^{\mathrm{int}}$ defined in an appropriate gauge MartinPendas ; Nielsen ; Godfrey and the contribution of the external body forces acting on electrons with the spatial density $f_{\mu}^{\mathrm{ext}}=-\rho_{1}(\mathbf{r},\mathbf{r})\nabla_{\mu}U_{\mathrm{ext}}(\mathbf{r})$. According to the stress theorem Nielsen ; Maranganti ; Godfrey ; Bader , the spatial average of (34) should be equal to the thermodynamic pressure:

[TABLE]

III Massless Dirac electrons

III.1 Generalized virial and Hellmann–Feynman theorems

Now we will turn to massless Dirac electrons in a solid. As a specific example, we consider two-dimensional system of electrons in graphene, but our general theorems should be applicable to any other Dirac materials. The massless electrons in graphene CastroNeto have the following distinctions from the massive ones, important for our analysis: 1) their effective (single-particle) wave function $\psi=(\psi_{A\mathbf{K}},\psi_{B\mathbf{K}},\psi_{B\mathbf{K}^{\prime}},\psi_{A\mathbf{K}^{\prime}})^{T}$ is a multi-component column with the components corresponding to the sublattices $A,B$ and valleys $\mathbf{K},\mathbf{K}^{\prime}$; 2) the kinetic part of the Hamiltonian is $H_{\mathrm{kin}}=v_{\mathrm{F}}\boldsymbol{\Sigma}\cdot\mathbf{p}$, where $\boldsymbol{\Sigma}=\mathrm{diag}(\boldsymbol{\sigma},-\boldsymbol{\sigma})$ is the $(4\times 4)$ vector matrix with the vectors composed of Pauli matrices on the diagonal; 3) the boundary conditions imposed on the wave function are not the Dirichlet condition but have more diverse forms of the system of equations $M\psi=\psi$, with the $(4\times 4)$ matrices $M$ dependent on the edge types Brey ; Akhmerov ; Berry ; McCann ; 4) to define the ground state and make the system energy bounded from below, an appropriate momentum cutoff $|\mathbf{p}|<p_{\mathrm{c}}$ deep in the valence band should be introduced; this approach allows us to approximate the tight-binding model of electrons in graphene, which have a physically bounded valence band, by a simpler effective model.

Now we can repeat the calculations of Section II with taking into account that for massless Dirac electrons the single-particle

[TABLE]

and many-body

[TABLE]

expressions for the generalized current are different from (3), (27) due to different form of the non-Hermitian $H_{\mathrm{kin}}$. The single-particle density matrix $\rho_{1}$ in (37) is assumed to be a $(4\times 4)$ matrix over sublattice and valley degrees of freedom. Taking also into account that $H_{\mathrm{kin}}$ is linear in momentum and assuming Coulomb interaction in a many-body system, that imply $[H,G]=i\hbar(-H_{\mathrm{kin}}-V_{\mathrm{int}}+\mathbf{r}\cdot\nabla U_{\mathrm{ext}})$, we get the counterpart of the virial theorem (28) for massless Dirac electrons:

[TABLE]

Eq. (38) can be compared with the generalized virial theorem, obtained in Sokolik by means of scaling transformations of a many-body wave function with the imposed momentum cutoff:

[TABLE]

(see the alternative derivation in B). Here, as in the previous sections, $\mathcal{P}=-\partial E/\partial\Omega$, and $E$ should be understood as the ground state energy $E_{0}$ at $T=0$ or as the free energy $F$ at $T>0$. Comparing (38) and (39), we obtain for graphene

[TABLE]

To draw connection between the boundary term in (40) and the physical pressure caused by electron collisions with the boundary, as in Section II.5, we again need to consider the boundary relationship for the wave function.

III.2 Boundary relationship and local pressure

Due to the specific form of boundary condition $M\psi=\psi$, the boundary relationship for massless Dirac electrons will be different from (13) for massive electrons. We can consider $M\psi-\psi$ as the four-component function satisfying the Dirichlet condition on the boundary, so Eq. (13) with the replacement $\psi\rightarrow M\psi-\psi$ can be applied in this case:

[TABLE]

The matrix $M$ should be unitary, $M^{+}=M^{-1}$, and anticommuting with the normal component $\Sigma_{n}\equiv\boldsymbol{\Sigma}\cdot\mathbf{n}$ of the probability current operator, $\{M,\Sigma_{n}\}=0$, to ensure that the particles do not cross the boundary, i.e. $\mathbf{J}[1]\cdot\mathbf{n}=v_{\mathrm{F}}\psi^{+}\Sigma_{n}\psi=0$ Akhmerov ; McCann . Using these properties of $M$, we get $\psi^{+}\Sigma_{n}\delta\psi=\psi^{+}M\Sigma_{n}\delta\psi=-\psi^{+}\Sigma_{n}M\delta\psi$. Applying (41) and again the condition $\{M,\Sigma_{n}\}=0$, we obtain

[TABLE]

This is the counterpart of the boundary relationship (13) for massless Dirac electrons.

If the single-particle stationary state $\psi$ is unaffected by the momentum cutoff both before and after the boundary perturbation, then the energy change of this state can be found by using (2), (7), (36), and Gauss theorem:

[TABLE]

Using (42) and introducing the kinetic stress tensor for massless Dirac particles,

[TABLE]

we get the formulas, which are fully analogous to (21)–(23). In this derivation we did not used any specific form of $M$, requiring only the absence of particle flux through the boundary, thus its results should be applicable to any bounded system of massless Dirac electrons.

III.3 Thermodynamic pressure: kinetic and anomalous parts

Generalization of the results of Section III.2 for a many-body system should be done with caution because of the presence of the momentum cutoff. A generic many-body wave function $\Psi_{n}$ can be presented as a sum of factorized wave functions

[TABLE]

where each $\psi_{i}^{(k)}$ satisfies the boundary conditions and does not need to be an eigenfunction of the Hamiltonian. When the boundary is perturbed, several transformation occur with this function: first, the single-particle wave functions $\psi_{i}^{(k)}$, which satisfy the momentum cutoff both before and after perturbation, are changed by the values $\delta\psi_{i}^{(k)}$ obeying (42) on the boundary. Second, some terms in (45) disappear because one or several of $\psi_{i}^{(k)}$ in these terms cease to satisfy the momentum cutoff condition, and some new terms with $\psi_{i}^{(k)}$ satisfying the momentum cutoff condition after the perturbation can appear instead of the disappeared ones. Consequently, the perturbation $\delta\Psi_{n}=\delta_{\mathrm{def}}\Psi_{n}+\delta_{\mathrm{c}}\Psi_{n}$ can be presented as a sum of

[TABLE]

which is caused by deformations of the single-particle functions $\psi_{i}^{(k)}$, and $\delta_{\mathrm{c}}\Psi_{n}$, which is caused by momentum cutoff.

Two contributions to $\delta\Psi_{n}$ will result, through (7), in two parts of the energy response $\delta E$ to the volume change $\delta\Omega$, and, correspondingly, in two parts of the pressure. Assuming a uniform dilation of the system boundary $\delta\mathbf{r}_{0}=(\delta R/R)\mathbf{r}_{0}$ and applying (42) for each $\delta\psi_{i}^{(k)}$ in (46), we obtain the energy change due to single-particle wave function deformations: $\delta_{\mathrm{def}}E=-(\delta R/R)\oint_{\partial\Omega}d\mathbf{s}\cdot\mathbf{J}[G]$, which is analogous to that for massive electrons. We can call the corresponding part of the thermodynamic pressure $\mathcal{P}=\mathcal{P}_{\mathrm{kin}}+\mathcal{P}_{\mathrm{anom}}$ as kinetic pressure

[TABLE]

which is caused, from the physical point of view, by a transfer of momentum to the boundary during electron reflections. The second part of $\mathcal{P}$, according to (40), is equal to

[TABLE]

and can be called anomalous pressure. It is caused by redistribution of electron states during a change of $\Omega$ due to the presence of the cutoff, which provide an additional contribution to the energy change.

Thus, in contrast to massive electrons, where the thermodynamic and kinetic pressures are equal Marc , in the case of massless Dirac electrons they differ by $\mathcal{P}_{\mathrm{anom}}$, and the formulas (23), (24) and (35) are applicable only to $\mathcal{P}_{\mathrm{kin}}$.

We should note that the breaking of the equality between $\mathcal{P}$ and $\mathcal{P}_{\mathrm{kin}}$ is not restricted to the Dirac model, but can appear in any solid with a filled valence band. The previously considered Dirac model is the important example of such systems because the valence band cutoff is not just a formal tool used to bound the system energy from below, but has observable consequences, e. g., logarithmic renormalization of the Fermi velocity of electrons in graphene in induced by Coulomb interaction CastroNeto . Also, the Dirac model is convenient to work with because it uniformly describes both valence and conduction bands using just a single parameter, $v_{\mathrm{F}}$. However we can extend our consideration to the system with quadratic electron dispersion and energy gap between conduction and valence bands, which is the simplest model of semiconductor or insulator. In C we calculate kinetic and anomalous parts of the pressure for this model.

IV Examples

IV.1 Free noninteracting electrons in thermodynamic limit

Let us consider the simplest example of free two-dimensional noninteracting electrons occupying the single-particle states $\psi_{\mathbf{p}\gamma}(\mathbf{r})=e^{i\mathbf{p}\cdot\mathbf{r}/\hbar}(1,\gamma e^{i\varphi_{\mathbf{p}}})^{T}/\sqrt{2\Omega}$ with the energies $\epsilon_{\mathbf{p}\gamma}=\gamma v_{\mathrm{F}}p$ in the conduction ($\gamma=+1$) and valence ($\gamma=-1$) bands in the $\mathbf{K}$ valley of graphene. In the case of electron doping, the states in conduction and valence bands are filled up to the Fermi $p_{\mathrm{F}}$ and cutoff $p_{\mathrm{c}}$ momenta respectively [see Fig. 2(a)]. Without attributing the exact form of boundary conditions, we can reasonably assume that in the sample of the linear size $R$ the momenta $\mathbf{p}$ are quantized in the units of $2\pi\hbar/R$. This neglect of the boundary behavior should be justified in the thermodynamic limit of a large system. The total energy $E$ and number $N$ of electrons can be calculated in the thermodynamic limit by transforming sums over momenta into integrals:

[TABLE]

where the sign $s_{\mu}$ of the chemical potential distinguishes the cases of electron ($s_{\mu}=+1$) or hole ($s_{\mu}=-1$) doping, and $g=4$ is the degeneracy over valleys and spin projections. To calculate the pressure $\mathcal{P}=-(\partial E/\partial\Omega)_{N}$, we need to consider simultaneous changes of $\Omega\propto R^{2}$ and $p_{\mathrm{F}}$ which preserve $N$, and the result is:

[TABLE]

On the other hand, we can calculate kinetic and anomalous parts of $\mathcal{P}$ separately. Using (47), we find that each electron state contributes $\epsilon_{\mathbf{p}\gamma}/2\Omega$ to $\mathcal{P}_{\mathrm{kin}}$, so

[TABLE]

Note that electrons in the valence band provide large negative contribution to $\mathcal{P}_{\mathrm{kin}}$ because they have negative group velocity that implies negative momentum transfer to the boundary on collisions. The anomalous part of pressure can by calculated from (48) with taking into account that $p_{\mathrm{c}}$ and $p_{\mathrm{F}}$ should change simultaneously to preserve $N$:

[TABLE]

The sum of (52) and (53) gives (51) in agreement with the generalized virial theorem.

The origin of $\mathcal{P}_{\mathrm{anom}}$ can be traced by looking at Fig. 2 showing the case $s_{\mu}=+1$. When $R$ is slightly increased, the momentum quantization interval decreases, that shifts the energies of the occupied states closer to the Dirac point and results in the total energy change responsible for the kinetic pressure: $\delta E=-\mathcal{P}_{\mathrm{kin}}\delta\Omega$. Since $p_{\mathrm{c}}\gg p_{\mathrm{F}}$, this $\delta E$ will be positive, hence $\mathcal{P}_{\mathrm{kin}}<0$. However, due to the same decrease of the quantization interval, new unoccupied states appear at the bottom of the valence band [Fig. 2(b)]. To maintain the ground state, the system should fill these states with electrons taken from the Fermi level [Fig. 2(c)]. The number of transferred electrons is proportional to $p_{\mathrm{c}}^{2}$ and their energy changes are $-v_{\mathrm{F}}(p_{\mathrm{F}}+s_{\mu}p_{\mathrm{c}})$, thus we obtain the additional negative contribution $\delta E=-\mathcal{P}_{\mathrm{kin}}\delta\Omega$ to the energy change, giving rise to the positive anomalous pressure (53).

The nonlinear corrections to the dispersion away from the Dirac point affect the pressure quantitatively, as is shown in more detail in D, but preserve the general picture of two physically different contributions to the thermodynamic pressure.

IV.2 Interacting uniform electron gas in graphene

In the case of spatially uniform electron gas in graphene with the Coulomb interaction, the expressions (51)–(53) acquire interaction-induced corrections. Writing the energy $E=\Omega\epsilon(n,p_{\mathrm{c}})$ and electron number $N=\Omega n$ in terms of energy $\epsilon$ and electron $n$ densities, we can rewrite the generalized virial theorem (39) in the form

[TABLE]

whereas the pressure components (47)–(48) and the total pressure become

[TABLE]

The theorem (54) is equivalent to $\mathcal{P}=\mathcal{P}_{\mathrm{kin}}+\mathcal{P}_{\mathrm{anom}}$.

In the low-energy physics of graphene CastroNeto only the properties of electron gas at low doping levels $p_{\mathrm{F}}\ll p_{\mathrm{c}}$ are observable. To describe the properties of such low-doping Dirac electron gas, we introduce the regularized energy and electron densities, obtained after subtraction of the valence band contribution: $\epsilon_{\mathrm{r}}=\epsilon-\epsilon_{0}-\mu_{0}n_{\mathrm{r}}$, $n_{\mathrm{r}}=n-n_{0}=s_{\mu}gp_{\mathrm{F}}^{2}/4\pi\hbar^{2}$, where $n_{0}=gp_{\mathrm{c}}^{2}/4\pi\hbar^{2}$ is the density of electrons in the filled valence band, $\epsilon_{0}=\epsilon(n_{0},p_{\mathrm{c}})$ is the energy density of the filled valence band, $\mu_{0}=(\partial\epsilon/\partial n)_{p_{\mathrm{c}}}|_{n=n_{0}}$ is the chemical potential of the electron gas at $n_{\mathrm{r}}=0$. Thus at $n_{\mathrm{r}}=0$, when the Fermi level is located in the Dirac point, we have $\epsilon_{\mathrm{r}}=0$ and the regularized chemical potential $(\partial\epsilon_{\mathrm{r}}/\partial n_{\mathrm{r}})_{p_{\mathrm{c}}}=0$. Switching from $\epsilon(n,p_{\mathrm{c}})$ to $\epsilon_{\mathrm{r}}(n_{\mathrm{r}},p_{\mathrm{c}})$ and taking into account that, according to dimensionality, $\varepsilon_{0}\propto p_{\mathrm{c}}^{3}$ and $\mu_{0}\propto p_{\mathrm{c}}$, we obtain the regularized version of the generalized virial theorem, previously discussed in Ref. Sokolik :

[TABLE]

Similarly to (55), we introduce the regularized kinetic, anomalous and total pressures,

[TABLE]

so the theorem (56) states $\mathcal{P}_{\mathrm{r}}=\mathcal{P}_{\mathrm{r,kin}}+\mathcal{P}_{\mathrm{r,anom}}$. Despite the similarity of (54)–(55) and (56)–(57), the quantities entering these formulas are very different in value and dependencies on system parameters.

To reveal the meaning of $\mathcal{P}_{\mathrm{r,kin}}$ and $\mathcal{P}_{\mathrm{r,anom}}$, we can use the scaling form of energy density: $\epsilon_{\mathrm{r}}=\epsilon_{\mathrm{r}}^{(0)}f(\Lambda,r_{\mathrm{s}})$, where $\epsilon_{\mathrm{r}}^{(0)}=gv_{\mathrm{F}}p_{\mathrm{F}}^{3}/6\pi\hbar^{2}$ is the regularized energy density of noninteracting gas, $\Lambda=p_{\mathrm{c}}/p_{\mathrm{F}}$ is the dimensionless cutoff momentum, and $r_{\mathrm{s}}=e^{2}/\varepsilon\hbar v_{\mathrm{F}}$ is the Coulomb interaction scale (“fine structure constant”) for graphene, $\varepsilon$ is the dielectric constant of surrounding medium. In terms of $f$, we obtain $\mathcal{P}_{\mathrm{r,anom}}=-(\epsilon_{\mathrm{r}}^{(0)}/2)\Lambda(\partial f/\partial\Lambda)$. As we see, $\mathcal{P}_{\mathrm{r,anom}}$ is caused by the cutoff dependence of the energy, which appears only in the presence of Coulomb interaction, because $f=1$ and $\mathcal{P}_{\mathrm{r,anom}}=0$ at $r_{\mathrm{s}}=0$ (in contrast to $\mathcal{P}_{\mathrm{anom}}$, which is nonzero even at $r_{\mathrm{s}}=0$, see (53)). As shown in Ref. Sokolik , $\mathcal{P}_{\mathrm{r,anom}}$ measures the extent of scale invariance breaking due to the cutoff in the theorem (56). So, even after subtraction of the valence band contributions to energy and chemical potential, the electron gas properties continue to depend on the cutoff momentum $p_{\mathrm{c}}$ in the presence of interaction.

As a specific example, we can take the Hartree-Fock approximation, in which function $f$ can be expanded at large $\Lambda$ with sufficient accuracy as Lozovik ; Peres

[TABLE]

where $\mathcal{C}\approx 0.916$ is Catalan’s constant. In more accurate random-phase approximation Lozovik ; Barlas , the coefficients at $\ln\Lambda$ and $\Lambda^{-n}$ ($n\geqslant 0$) acquire additional nonlinear $r_{\mathrm{s}}$ dependencies. If, generally, $f=1+A(r_{\mathrm{s}})\ln\Lambda+\mathcal{O}(1/\Lambda)$ ($A=r_{\mathrm{s}}/4$ in the Hartree-Fock approximation), then the leading $\ln\Lambda$ term is related to the quantity

[TABLE]

which was discussed in Ref. Sokolik and can be evaluated from experimental data on graphene electron compressibility or quantum capacitance. At $\Lambda\rightarrow\infty$, we have $K=A(r_{\mathrm{s}})+\mathcal{O}(1/\Lambda)$. The numerical calculations show that $K$ is nearly constant in the range of doping levels of graphene accessible by using the electric field effect Sokolik . Thus the regularized anomalous pressure $\mathcal{P}_{\mathrm{r,anom}}$ can be related to experimental data and to the logarithmic term in $f$. Note that, as seen from (58), $\mathcal{P}_{\mathrm{r,kin}}$ is also changed in the presence of interaction due to renormalization of electron Fermi velocity CastroNeto .

Both regularized and nonregularized thermodynamic pressures are connected with the observable quantum capacitance of graphene $C_{\mathrm{Q}}$ Lozovik :

[TABLE]

In principle, we can calculate the total, unregularized pressures $\mathcal{P}$, $\mathcal{P}_{\mathrm{kin}}$, and $\mathcal{P}_{\mathrm{anom}}$ using the known dependence $C_{\mathrm{Q}}^{-1}(n_{\mathrm{r}})$:

[TABLE]

and $\mathcal{P}_{\mathrm{anom}}=\mathcal{P}-\mathcal{P}_{\mathrm{kin}}$. Using the power series expansion of $C_{\mathrm{Q}}^{-1}$ for the interacting system, similar to (58), we can perform the integrations analytically and obtain the interaction-induced corrections to the unregularized pressures, caused by corresponding corrections to $C_{\mathrm{Q}}^{-1}$ studied in Ref. Lozovik . However the formulas (61)–(62) contain the parameters $n_{0}$, $\epsilon_{0}$, and $\mu_{0}$ of the filled valence band, which can be the sources of additional interaction-induced corrections. In estimating these parameters, we also need to take into account deviations from the linear dispersion at large electron momenta (D).

IV.3 Rectangular graphene flake

Here we consider the single-particle states of massless Dirac electrons in a rectangular graphene sample with zigzag horizontal edges and armchair vertical edges, imposing the corresponding boundary conditions: $\psi_{A\mathbf{K}}=\psi_{A\mathbf{K}^{\prime}}=0$ at the bottom edge, $\psi_{B\mathbf{K}}=\psi_{B\mathbf{K}^{\prime}}=0$ at the top edge, $\psi_{A\mathbf{K}}+\psi_{A\mathbf{K}^{\prime}}=\psi_{B\mathbf{K}}+\psi_{B\mathbf{K}^{\prime}}=0$ on the left edge, $e^{2\pi\nu i}\psi_{A\mathbf{K}}+\psi_{A\mathbf{K}^{\prime}}=e^{2\pi\nu i}\psi_{B\mathbf{K}}+\psi_{B\mathbf{K}^{\prime}}=0$ on the right edge, where $\nu=\pm 2/3$ or 0 depending on the atomic-scale details Brey ; Akhmerov . Hereafter we set $\hbar\equiv 1$, $v_{\mathrm{F}}\equiv 1$ in the formulas; in numerical calculations, we take $v_{\mathrm{F}}=10^{6}\,\mbox{m/s}$ and such cutoff $p_{\mathrm{c}}$ that the filled valence band has two electrons per unit cell of graphene, which corresponds to the cutoff energy around 7.2 eV. In a $L_{x}\times L_{y}$ rectangle, the (not normalized) eigenstates $\psi=(e^{ik_{x}x}\sin k_{n}y$, $\mp(-1)^{n}e^{ik_{x}x}\sin k_{n}(L_{y}-y)$, $\pm(-1)^{n}e^{-ik_{x}x}\sin k_{n}(L_{y}-y)$, $-e^{-ik_{x}x}\sin k_{n}y)^{T}$ with energies $\epsilon=\pm\sqrt{k_{x}^{2}+k_{n}^{2}}$ are determined by quantum numbers $m\in\mathbb{Z}$ and $n$, where $k_{x}=(-\frac{2}{3}\nu+m)\pi/L_{x}$ (we take $\nu=0$), and $k_{n}$ is the $n$th positive root of the equation $k_{n}=-k_{x}\tan k_{n}L_{y}$. The local pressure on the boundaries as defined in (23) is constant at the zigzag edges and oscillates along the armchair ones. These oscillations occur because the zigzag edges preserve the valley of the incident wave so in the direction normal to them a standing wave pattern is formed, while the armchair edges change the valley so there is no interference between incident and reflected waves. Besides the local pressure on the boundaries, we calculate the local bulk pressure (34), which equals simply $\frac{1}{2}\epsilon\psi^{+}\psi$. It is constant over the $x$ direction and oscillates over the $y$ direction. The example of a typical pressure distribution for one single-particle state is shown in Fig. 3(a). If we consider the pressure distribution of a many-body system with many different states occupied, the oscillations of the pressure disappear, but the feature of zero boundary pressure in the $x$ direction at the angles of the flake is preserved.

To calculate the kinetic pressure of the many-body system according to (47), which in the noninteracting case reduces to $\mathcal{P}_{\mathrm{kin}}=E/2L_{x}L_{y}$, and the anomalous pressure (48), we need to consider the total energy $E$ as a function of $p_{\mathrm{c}}$ with a constant $N$. Due to discrete nature of the energy spectrum, the derivative $(\partial E/p_{\mathrm{c}})_{N}$ is the sum of Dirac delta functions, so they should be smoothened (replaced by Lorentzians in our calculations) to get sensible result for $\mathcal{P}_{\mathrm{anom}}$. Both contributions to the pressure are shown in Fig. 3(b). In the region of large momenta, where the energy levels are spaced densely enough, the picture is expectedly very similar to the case of free Dirac electron gas, described by Eqs. (52)–(53). Note that the total thermodynamic pressure is positive everywhere.

IV.4 Circular graphene flake in Coulomb potential

For a circular flake, we impose the infinite mass boundary condition $\psi_{B}=ie^{i\phi}\psi_{A}$, where $\phi$ is the polar angle of the direction normal to the boundary Berry , which decouple the valleys and allows considering only a single valley. To study the effects of the external potential in the circular geometry, we assume the presence of a Coulomb impurity in the center. Without external potential, the solutions of the Dirac equation are given by Bessel functions; for a subcritical Coulomb potential $-Ze^{2}/r\equiv\tilde{g}\hbar v_{\mathrm{F}}/r$, $|\tilde{g}|<1/2$, the solutions are given by the Coulomb wave functions Pereira . The quantum numbers defining a solution are the half-integer angular momentum $j$, the radial quantum number $n\in\mathbb{N}$, and the energy sign $\gamma$. The electron-hole symmetry of the system requires $\epsilon_{\gamma,j,n}(\tilde{g})=-\epsilon_{-\gamma,-j,n}(-\tilde{g})$. The local pressure on the boundary (23) is uniform due to the circular symmetry of the system. The distributions of the local bulk pressure (34), which now includes the contribution of the external Coulomb force, are shown in Fig. 4 for several single-particle states. As shown in Fig. 5(a), the states with the same quantum numbers have higher quantum pressure (both the total kinetic pressure and the local boundary pressure) at higher values of the Coulomb potential parameter $\tilde{g}$ in the agreement with Eq. (47): the repulsive potential increases the pressure by pushing the electrons towards the boundary, while the attractive potential decreases the pressure by pulling the electrons to the center.

Fig. 5(b) shows the kinetic, anomalous (smoothened with Lorentzians), and total pressure of Dirac electrons in a circular flake. In comparison with Fig. 3(b), here the smaller size of the flake leads to visible deviations of the anomalous pressure from the linear trend (53) near the Dirac point, but the agreement with thermodynamic limit (52)–(53) is restored at large Fermi momenta.

V Discussion

Using the generalized continuity equation and scaling transformations, we derived and analyzed the virial and Hellmann–Feynman theorems for single- and many-electron systems with taking into account the presence of system boundaries. The boundary conditions imposed on the wave function make the Hamiltonian generally non-Hermitian, which results in appearance of additional terms in (6), (9) in the form of boundary integrals. We start with the case of massive electrons and analyze the thermodynamic pressure as a response of a system energy on small volume changes and relate the pressure to the boundary term in the virial theorem (16)–(17) using the boundary relationships (13)–(14). Besides, we find the local pressure as a response (22) of the energy on local deformations of the boundary and connect it with the kinetic part of the stress tensor (20), (23). The formulas are first derived for a single-particle system and then generalized for a many-body system in Section II.6. While the most of these relationships for massive electrons can be found elsewhere Marc ; Ziesche ; Argyres ; Abad ; Esteve ; Cottrell ; Fernandez ; Bobrov ; Srebrenik , we presented them for the sake of completeness. The connection (22)–(23) between energy change and boundary perturbations is also known in the boundary perturbation theory of the boundary-value problems Henry .

For massless Dirac electrons in a solid the similar formulas are different in some aspects because of the linear dispersion, different forms of boundary conditions for a wave function and due to the presence of momentum cutoff in the valence band. The latter results in appearance of the cutoff-induced term in the generalized virial theorem (39), and the thermodynamic pressure (40) turned out to consist of two parts. The first part is the kinetic pressure (47), which is just a sum of responses of single-particle energies weighted with their occupation numbers. Using the boundary relationships (41), (42) for massless Dirac electrons, we can relate it, analogously to (23), to the kinetic stress tensor. Thus the kinetic pressure is caused by momentum transferred by electrons to the surroundings during their reflections from the boundary. The second part is the anomalous pressure (48), which is related to the momentum cutoff and caused by redistribution of electron states during a volume change, as shown in Fig. 2. Note that the problem of consistency between kinetic and thermodynamic definitions of the pressure, resolved for ordinary massive electrons Cottrell ; Fernandez ; Argyres with the help of the boundary relationships, rises again in the case of massless Dirac electrons because of the anomalous contribution.

The example of free electrons considered in Section IV.1 demonstrates that while the total kinetic pressure of an electron gas in graphene is negative due to overwhelming contribution of the valence band, it is overcompensated by the anomalous pressure, making the thermodynamic pressure positive. The examples of rectangular and circular graphene flakes with the zigzag, armchair and infinite effective mass boundary conditions for the wave functions demonstrate fulfilment of the general theorems.

The pressure $\mathcal{P}$ studied in this paper is related to the grand thermodynamic potential $\Omega=-\mathcal{P}\Omega$ of the electron gas, so taking its derivatives with respect to the parameters can provide all thermodynamic properties. The derivatives of $\mathcal{P}$ with respect to the electron density are related to such observable quantities as electron compressibility and quantum capacitance, which were studied in graphene in the context of interaction and disorder-induced effects Lozovik . In Section IV.2 we consider connection of $\mathcal{P}$, $\mathcal{P}_{\mathrm{kin}}$, $\mathcal{P}_{\mathrm{anom}}$ with quantum capacitance. Our analysis of a bounded system allows to extend these studies by including the effects of boundary conditions in small graphene flakes. The general theorems derived here can be applied not only to graphene, but also to three-dimensional Dirac and Weyl semimetals, which host massless Dirac electrons as well Armitage ; Wehling .

The changes of the volume of the system considered in this paper concern only electron subsystem and not the crystal lattice itself. We analyze what happens with the electrons described by the effective Dirac equation when we move the boundary conditions without deforming the lattice, so $p_{\mathrm{c}}=\mathrm{const}$ in these processes. However the other kind of the system volume changes, when graphene is compressed or stretched as a whole, can be considered. In this case the energy change can be related to a total mechanical stress and total mechanical compressibility of graphene. Analysis of such graphene properties should include, besides the considered response of Dirac electrons, the responses of the core electrons and atomic nuclei, which is beyond the scope of our paper.

The problem of breaking the equality between thermodynamic and kinetic pressures due to the anomalous contribution is not unique for Dirac particles and can arise in any solid state system with the filled valence band. Change of electron number in this band proportional to the change of enclosing volume requires electron transitions between valence and conduction bands, which provide anomalous contribution to the total pressure. As the recent study Nakata suggests, such transitions can proceed through the Tamm states inside the energy gap.

The momentum cutoff deep in valence band of graphene, which results in appearance of the “anomalous” terms in the generalized virial theorem and pressure, may be considered as an artificial construct, however in solids with massless Dirac electrons it has real physical grounds, because valence band is indeed bounded in energy and momentum spaces. In graphene it leads to a finite-valued logarithmic renormalization of the Fermi velocity due to Coulomb interaction CastroNeto . Nevertheless, more accurate analysis with going beyond the Dirac model and with taking into account Tamm states on the boundaries can provide more insight into the problem of electron gas pressure in solids with unusual band structure. The study of relationship between kinetic and thermodynamic pressures with taking into account other possible anomalies can be extended to a broader context of statistical physics of confined many-particle systems.

Acknowledgments

The work was supported by the grants No. 17-02-01134 and 18-52-00002 of the Russian Foundation of Basic Research. Yu.E.L. was partly supported by the Program for Basic Research of the National Research University Higher School of Economics. A.A.S. and A.D.Z. acknowledge the support from the Foundation for the Advancement of Theoretical Physics and Mathematics “BASIS”.

Appendix A Scaling relations for uniform system of massive electrons

Consider a system of $N$ massive electrons with Coulomb interaction $V_{\mathrm{int}}(r)=e^{2}/\varepsilon r$ in the external power-law potential $U_{\mathrm{ext}}(r)=U_{0}r^{\gamma}$, confined in the volume $\Omega$ by infinitely high potential walls. This system is described by the many-body Schrodinger equation

[TABLE]

On the boundary we impose the Dirichlet condition: $\Psi(\mathbf{r}_{1}\ldots\mathbf{r}_{N})=0$ when $\vee i:\,\mathbf{r}_{i}\in\partial\Omega$.

With the characteristic size of $\Omega$ being $R$, we can switch to the dimensionless coordinates $\boldsymbol{\rho}_{i}=\mathbf{r}_{i}/R$, and the wave function is scaled as $\Psi(\mathbf{r}_{1}\ldots\mathbf{r}_{N})=R^{-ND/2}\tilde{\Psi}(\boldsymbol{\rho}_{1}\ldots\boldsymbol{\rho}_{N})$. Introducing the dimensionless energy $\tilde{E}=mR^{2}E/\hbar^{2}$, interaction strength $r_{\mathrm{s}}=e^{2}mR/\varepsilon\hbar^{2}$, and external potential strength $\varkappa=U_{0}mR^{\gamma+2}/\hbar^{2}$, we obtain the scaled equation

[TABLE]

which does not depend on $R$ explicitly. The boundary conditions for $\tilde{\Psi}$ depend only on the shape of $\Omega$ and not on its size. As a result, we obtain the scaling forms

[TABLE]

Calculating derivatives of (65) and (66) with respect to $R$, we obtain the following scaling properties of $\Psi$ and $E$:

[TABLE]

Since $e^{2}(\partial E/\partial(e^{2}))=\langle V_{\mathrm{int}}\rangle$ and $\gamma U_{0}(\partial E/\partial U_{0})={\langle\mathbf{r}\cdot\nabla U_{\mathrm{ext}}\rangle}$, we can immediately identify (68) as the virial theorem (31) for the case of power-law $U_{\mathrm{ext}}(\mathbf{r})$. Eq. (67) should be valid in all points of space, so if one of the $\mathbf{r}_{i}$ is located on the boundary, $\partial\Psi/\partial(e^{2})$ and $\partial\Psi/\partial U_{0}$ vanish due to the Dirichlet boundary condition, so we obtain $(\sum_{i}\mathbf{r}_{i}\cdot\nabla_{i})\Psi=-R(\partial\Psi/\partial R)$ and hence the many-body counterpart of the boundary relationship (14).

Appendix B Scaling relations for uniform system of massless Dirac electrons

A system of $N$ massless Dirac electrons with Coulomb interaction $V_{\mathrm{int}}(r)=e^{2}/\varepsilon r$ in the external power-law potential $U_{\mathrm{ext}}(r)=U_{0}r^{\gamma}$, confined in the volume $\Omega$, is described by the many-body Dirac equation:

[TABLE]

Some boundary conditions of the kind $M_{i}\Psi=\Psi$, not specified explicitly here, are imposed on $\Psi(\mathbf{r}_{1}\ldots\mathbf{r}_{N})$ when $\vee i:\,\mathbf{r}_{i}\in\partial\Omega$. We should also impose the momentum cutoff condition $P_{p_{\mathrm{c}}}\Psi=\Psi$, where the operator $P_{p_{\mathrm{c}}}$ of projection on the subspace $|\mathbf{p}_{i}|\leq p_{\mathrm{c}}$ of momentum space was described in Ref. Sokolik .

As in A, we use the dimensionless coordinates $\boldsymbol{\rho}_{i}=\mathbf{r}_{i}/R$ and the scaled wave function $\Psi(\mathbf{r}_{1}\ldots\mathbf{r}_{N})=R^{-ND/2}\tilde{\Psi}(\boldsymbol{\rho}_{1}\ldots\boldsymbol{\rho}_{N})$. Introducing the dimensionless energy $\tilde{E}=RE/\hbar v_{\mathrm{F}}$, interaction constant $r_{\mathrm{s}}=e^{2}/\varepsilon\hbar v_{\mathrm{F}}$, and the external potential strength $\varkappa=U_{0}R^{\gamma+1}/\hbar v_{\mathrm{F}}$, we obtain the scaled Dirac equation:

[TABLE]

The boundary conditions for $\tilde{\Psi}$ are now independent on $R$, and the cutoff condition depends only on the dimensionless parameter $\Lambda=Rp_{\mathrm{c}}$. The resulting scaling forms of $\Psi$ and $E$ are

[TABLE]

Calculating derivatives of (71) and (72) with respect to $R$, we obtain the scaling properties:

[TABLE]

Taking into account that $\gamma U_{0}(\partial E/\partial U_{0})={\langle\mathbf{r}\cdot\nabla U_{\mathrm{ext}}\rangle}$ we obtain from (74) the generalized virial theorem (39). The equation (73) can be interpreted as a counterpart of (14) for the many-body wave function subject to momentum cutoff. The scaling analysis of a system of massless Dirac electrons can be also found in Lin .

Appendix C Massive electrons in a two-band model

Consider a two-band model with nonzero effective mass:

[TABLE]

This is the simplest model description of conduction and valence bands of a semiconductor or insulator separated by the gap $2\Delta$.

If the momentum cutoff at $p=p_{\mathrm{c}}$ is imposed in the valence band, it can be shown in the same way as in the Dirac case that the anomalous pressure arises in such system, too. In particular, in the two-dimensional noninteracting many-body system with degeneracy factor $g$, the number of particles $N$ is given by Eq. (50), while the energy, the thermodynamic pressure, and the kinetic and the anomalous contributions to it are given by the following expressions (the derivatives are taken at constant $N$):

[TABLE]

Here $s_{\mu}=\pm 1$ for, respectively, electron- and hole-doped material. Similarly to the case of noninteracting Dirac electrons (Section IV.1), at $p_{\mathrm{c}}\gg p_{\mathrm{F}}$ we have $\mathcal{P}_{\mathrm{kin}}<0$, but the total pressure is positive due to anomalous contribution.

Appendix D Corrections to pressure due to dispersion nonlinearities

For noninteracting Dirac electrons in graphene with perfectly linear dispersion, the kinetic and the anomalous pressure is given by the equations (52) and (53). However far away from the Dirac point the dispersion has nonlinear corrections CastroNeto . Here we consider the corrections to the pressure from these nonlinearities.

If the nearest-neighbor tight-binding Hamiltonian for the $2p_{z}$ orbitals of carbon atoms in graphene is expanded near the $\mathbf{K}$ point CastroNeto , the first term is linear in the momentum $p$, the next term is the trigonal warping proportional to $p^{2}\sin 3\varphi_{\mathbf{p}}$ which provides no contribution to the pressure after integration over the polar angle $\varphi_{\mathbf{p}}$, and the cubic term is: $\delta\epsilon_{\mathbf{p}\gamma}^{(3)}=-(7/64)\gamma v_{\mathrm{F}}p^{3}(d/\hbar)^{2}$ where $d$ is the interatomic distance. The corrections to the kinetic and anomalous pressure caused by $\delta\epsilon_{\mathbf{p}\gamma}^{(3)}$ are:

[TABLE]

where $\kappa=(p_{\mathrm{c}}d/\hbar)^{2}$. If we take such cutoff $p_{\mathrm{c}}$ that the filled valence band has two electrons per unit cell, as in Section IV.3, then $\kappa=4\pi/3\sqrt{3}$, and the kinetic pressure nearly halves in absolute value, while the anomalous pressure drops by one quarter near the Dirac point. The next-order correction is much less significant though.

We also consider the isotropic quadratic correction coming from the next-to-nearest-neighbor hopping, $\delta\epsilon_{\mathbf{p}\gamma}^{(2)}=\nu v_{\mathrm{F}}p^{2}d/\hbar$, where $\nu=3t_{2}/2t$, $t$ and $t_{2}$ are the nearest-neighbor and the next-to-nearest-neighbor hopping integrals. The corresponding pressure corrections are:

[TABLE]

If we take the upper bound $\nu=0.3$ CastroNeto , then the correction approximately doubles the kinetic pressure and increases the anomalous pressure by half.