## Abstract

Using the phase-space formulation of quantum mechanics, we derive a four-component Wigner equation for a system composed of spin- fermions (typically, electrons) including the Zeeman effect and the spin–orbit coupling. This Wigner equation is coupled to the appropriate Maxwell equations to form a self-consistent mean-field model. A set of semiclassical Vlasov equations with spin effects is obtained by expanding the full quantum model to first order in the Planck constant. The corresponding hydrodynamic equations are derived by taking velocity moments of the phase-space distribution function. A simple closure relation is proposed to obtain a closed set of hydrodynamic equations.

This article is part of the themed issue ‘Theoretical and computational studies of non-equilibrium and non-statistical dynamics in the gas phase, in the condensed phase and at interfaces’.

## 1. Introduction

The formulation of quantum mechanics in the phase space was first introduced by Eugene Wigner in 1932 to study quantum corrections to classical statistical mechanics [1]. The goal was to link the wave function that appears in the Schrödinger equation to a pseudo-probability distribution function defined in the classical phase space. This pseudo-probability distribution changes in time according to an evolution equation (Wigner equation) which is somewhat similar to the classical Liouville equation. Mathematically speaking, the Wigner formulation is based on the Weyl transformation [2,3], which is a general method to transform operators defined in the Hilbert space into phase-space functions.

As it is based on the classical phase space, the Wigner formulation is often a more intuitive approach than the standard Schrödinger equation, especially for problems where semiclassical considerations are important. For these reasons, it is used in many areas of quantum physics, including quantum optics [4], semiclassical analysis [5,6], electronic transport [7], nonlinear electron dynamics [8] and quantum plasma theory [9]. It is also the starting point to construct quantum hydrodynamic equations, which are approximate models obtained by taking velocity moments of the Wigner function. Such models were used in the past to study the electron dynamics in molecular systems [10], metal clusters and nanoparticles [11–13], thin metal films [14], quantum plasmas [15,16] and semiconductors [17].

The works cited above only considered the charge dynamics and disregarded the spin degrees of freedom. However, it is well known that spin effects (particularly the Zeeman splitting and spin–orbit coupling) can play a decisive role in nanometric systems such as semiconductor quantum dots [18,19] and diluted magnetic semiconductors [20,21]. The coupling between the spin degrees of freedom and the electron orbital motion is of the utmost importance in many experimental studies involving magnetized nano-objects. A particularly interesting example is the ultrafast demagnetization induced by a femtosecond laser pulse in ferromagnetic thin films [22], an effect that is not yet completely elucidated from the theoretical viewpoint. Recent time-dependent density functional theory (TDDFT) simulations suggest that the spin–orbit coupling plays a central role in the demagnetization process [23].

Phase-space models based on the Boltzmann equation [24]—and the corresponding fluid models [25]—were derived in the past to describe the dynamics of a gas where the constituents possess internal degrees of freedom (internal angular momentum). However, in these models, the spin is not treated *ab initio* as a fundamental quantity, but is rather incorporated into the transport equations to ensure the correct conservation properties. More recently, a few theoretical models that include the spin in the Wigner formalism were also developed. One approach [26] consists of defining a scalar probability distribution that evolves in an extended phase space, where the spin is treated as a classical two-component variable (related to the two angles on a unit-radius sphere) on the same footing as the position or the momentum. This approach was used to derive a Wigner equation that incorporates spin effects through the Zeeman interaction [26]. Semiclassical [27] and hydrodynamic [28] spin equations were also derived from those models, including other relativistic effects such as the spin–orbit coupling, the Darwin term and the relativistic mass correction.

An alternative approach keeps the 2×2 matrix character of the distribution function [29], so that the orbital and spin dynamics are represented by different Wigner functions. Using this approach, the corresponding Wigner equations were derived from the full Dirac theory [30]; however, their complexity makes them unsuitable for applications to condensed matter and nanophysics. A more tractable Wigner equation was derived from the Pauli (instead of Dirac) theory, but only included the Zeeman effect [31].

Both approaches (extended phase space and matrix Wigner function) are equivalent from the mathematical point of view. However, the extended phase-space approach leads to cumbersome hydrodynamic equations that are, in practice, very hard to solve, either analytically or numerically, even in the non-relativistic limit. The matrix technique, which is the one used here, separates clearly the orbital motion from the spin dynamics and leads to a simpler and more transparent hydrodynamic model.

In this paper, we go beyond our previous work [31] by including both the Zeeman effect and the spin–orbit coupling, the latter being a relativistic effect to second order in 1/*c*. In terms of a semiclassical expansion, these terms are respectively first and second order in . Other relativistic corrections such as the Darwin term or the mass correction are neglected here, although they could be included with relative ease in our model.

First, we use a gauge-invariant formulation of the Weyl transformation [32] and the Moyal product [33] to derive a set of Wigner equations describing a system of spin- fermions. A self-consistent mean-field model can further be obtained by coupling these Wigner equations to the set of Maxwell equations for the electromagnetic fields, where the sources (charge and current densities) are related to velocity moments of the Wigner function. A related mean-field model was obtained recently by Dixit *et al.* [34] in the framework of the Schrödinger–Pauli equation.

Subsequently, we derive the corresponding semiclassical limit and obtain the Vlasov equations describing the evolution of an electron gas with spin and semirelativistic effects. In this model, the orbital dynamics is treated classically, whereas the spin is represented as a fully quantum variable. Finally, the Vlasov equations will be used to derive a hierarchy of hydrodynamic equations by taking velocity moments of the probability distribution function. This is an infinite hierarchy that needs to be closed by using some additional physical hypotheses. Although this is relatively easy for spinless systems (where the closure can be obtained by assuming a suitable equation of state), things are subtler when the spin degrees of freedom are included. Here, we shall use an intuitive closure to obtain a closed set of fluid equations with spin effects.

## 2. Quantum mechanics in the phase space

The basic idea of the phase-space formulation of quantum mechanics is to associate at each operator , depending on the position and momentum operators and , a function of the classical phase-space variables ** r** and

**. This correspondence is provided by the Weyl transformation [3,35], and is given by 2.1where is the Wigner operator defined as 2.2The inverse of the Weyl transformation can be deduced**

*p*^{1}from the above definition: 2.3where tr denotes the trace.

Considering a system in a statistical distribution described by the density operator , equation (2.1) can be used to determine the mean value of an arbitrary operator:
The Wigner function is then defined as the phase-space function associated with the density operator
2.4The Wigner function obeys the following equation of motion (Wigner equation):
2.5where the last term is referred to as the Moyal bracket
2.6The superscripts *L* and *R* mean that the derivative acts only on the left or on the right term in the parenthesis (*i*=*x*,*y*,*z* and here and in the following we use Einstein’s summation convention). The Wigner equation is the analogue of the density matrix evolution equation in the operator representation of quantum mechanics: , sometimes called the von Neumann equation. The Moyal brackets can be easily developed as a power series of , which makes the Wigner formulation particularly interesting to study the semiclassical limit. The lowest-order term leads to the standard Poisson brackets and to the equations of classical mechanics.

In equation (2.5), is the phase-space function associated with the Hamiltonian operator of the system, and they are related by equation (2.1). In order to determine the phase-space function of any arbitrary operator , one should apply the Weyl correspondence rule [3,32], defined as follows: (i) first symmetrize the operator with respect to the position and the momentum operators and , (ii) then replace by their associated classical variables. For instance, for the operator one finds 2.7where use has been made of the commutator . We note that the Weyl correspondence defined above is not unique, and one could have defined other rules leading to a different phase-space function, such as the Husimi representation [36].

In the case of a spinless particle moving in a scalar potential *V* (** r**), the Wigner evolution (equation (2.5)) reads as
2.8However, complications arise when we want to include magnetic interactions. It is well known that in the presence of magnetic fields one should use the kinetic momentum operator instead of (

*q*=−

*e*,

*e*>0 for an electron). This situation cannot be addressed by simply replacing with in the Weyl transformation. Indeed, it can be easily proven that with such substitution the Wigner function, equation (2.4), is not gauge invariant. As spin effects, such as the Zemann interaction or the spin–orbit coupling, strongly depend on the magnetic field, it is of paramount importance to work with a gauge-invariant formulation of the Weyl transformation. A gauge-independent definition of the Wigner function was first introduced by Stratonovich [37]: 2.9where the momentum

**was replaced by .**

*p*To be consistent with this new definition of the Wigner function, one should also modify the Weyl correspondence rule [32]. The procedure is identical, except that one should use instead of . The main difference is that one must also symmetrize operators with respect to the different components of , because they do not commute, i.e. , where *ϵ*_{ijk} is the Levi–Civita symbol and is the magnetic field operator (). The classical phase-space variable associated with the kinetic momentum operator is the linear momentum .

The Moyal product defined in equation (2.6) is also modified in the presence of magnetic fields. A gauge-invariant Moyal product was derived by Müller [38], and reads
2.10where is the operator corresponding to the free magnetic field case:
2.11and is a new operator that depends on the magnetic field:
2.12with *g*(*n*,*p*) = [(1 − (−1)^{p})(*n* + 1) − (1 − (−1)^{n+1})*p*],(*r*_{1},*r*_{2},*r*_{3}) = (*x*,*y*,*z*) and (*π*_{1},*π*_{2},*π*_{3}) = (*π*_{x},*π*_{y},*π*_{z}). This new definition of the Moyal product makes the calculation of the evolution equation much more cumbersome than in the unmagnetized case. Its great advantage is that it ensures that the final equations of motion are gauge invariant.

## 3. Derivation of the spin Wigner model

We consider an ensemble of fermions in the presence of an electromagnetic field ** E**,

**. We denote the Schrödinger wave function of the**

*B**μ*th particle state by 3.1where and are respectively the spin-up and spin-down components of the wave function. The evolution of the system is governed by the Pauli–Schrödinger equation 3.2and 3.3Here, is the Bohr magneton,

**is the vector made of the 2×2 Pauli matrices, and**

*σ**σ*

_{0}is the 2×2 identity matrix.

*V*and

**are, respectively, the scalar and vector potential. Equation (3.3) can be derived from the Dirac equation by means of a Foldy–Wouthuysen transformation [39,40]. This semirelativistic development leads to plenty of terms that couple the spin to the charge dynamics. In this work, we keep terms only up to second order in 1/**

*A**c*, where

*c*is the speed of light in vacuum, namely the Zeeman interaction (order 1/

*c*

^{0}) and the spin–orbit coupling (order 1/

*c*

^{2}). We neglect, however, the Darwin term and the relativistic mass correction, which are also second-order effects.

Without spin, the Wigner function is a scalar function related to the density matrix *ρ* through equation (2.9). This definition can be generalized as follows to take into account the spin degrees of freedom:
3.4where, for particles with spin , *F* is a 2×2 matrix. The elements of the density matrix *ρ*^{ηη′}(** r**,

**′,**

*r**t*), where

*η*=↑,↓, are given by 3.5In order to study the macroscopic properties of the system, it is convenient to project

*F*onto the Pauli basis set [41,42] 3.6where 3.7

and tr denotes the trace. With this definition, the particle density *n* and the spin polarization ** S** of the electron gas are easily expressed by the moments of the pseudo-distribution functions

*f*

_{0}and

**: 3.8and 3.9**

*f*In this representation, the Wigner functions have a clear physical interpretation: *f*_{0} is related to the total electron density (in phase space), whereas *f*_{k} (*k*=*x*,*y*,*z*) is related to the spin polarization in the direction *k*. In other words, *f*_{0} represents the probability of finding an electron at one point of the phase space at a given time, while *f*_{k} represents the probability of that electron having a spin-polarization probability in the direction *k*.

There exist different ways to include the spin in the Wigner formalism other than the one described above. For instance, Brodin and co-workers [26] introduced an extended phase space (** r**,

**,**

*v***), where**

*s***is a unitary vector that defines the spin direction. The corresponding probability distribution is a scalar function of the extended phase-space variables. This is in contrast with our approach, where the spin is treated as a fully quantum variable (evolving in a two-dimensional Hilbert space). Nevertheless, the two approaches are equivalent, as shown in [31]. The correspondence relations between our distribution functions**

*s**f*

_{0}(

**,**

*r***,**

*v**t*) and

**(**

*f***,**

*r***,**

*v**t*) and the scalar distribution used by Zamanian

*et al.*[26]

*f*

_{Z}(

**,**

*r***,**

*v***,**

*s**t*) read as 3.10

Let us now turn to the evolution equation for the Wigner functions *f*_{0}(** r**,

**,**

*v**t*) and

*f*

_{k}(

**,**

*r***,**

*v**t*). After some tedious calculations, developed in the electronic supplementary material, equation (2.5) leads to the following Wigner equations: 3.11 3.12where depends on the magnetic field and corresponds to a quantum shift of the velocity 3.13and and are written in terms of the electric and magnetic fields 3.14The subscripts ± mean that the corresponding quantity is evaluated at a shifted position . This particularly illuminating form of the Wigner equations was proposed by Serimaa

*et al.*[32] for the case of a charged particle without spin evolving in an external electromagnetic field.

Equations (3.11) and (3.12) can be viewed as a generalization of those obtained in our previous work [31], where only the Zeeman interaction was included. The latter has two effects: the first is to couple the spin to the orbital dynamics through the gradient of the magnetic field (terms in the equations); the second effect is the spin precession around an effective magnetic field [terms (*B*_{+}+*B*_{−})×** f**]. In addition, many new terms appear owing to the spin–orbit interaction, which can be easily identified because they are proportional to 1/

*c*

^{2}. Some of these terms couple the spin to the orbital dynamics, whereas others provide corrections to the spin precession or the Lorentz force. The physical origin of all these terms will appear clearly in the next session, when we discuss the semiclassical limit of the Wigner equations.

Equations (3.11) and (3.12) can be used, in the context of a mean-field approach, to describe the self-consistent spin dynamics of an ensemble of interacting electrons. In this case, the electric and the magnetic fields are solutions of the Maxwell equations:
3.15where we introduced some relativistic corrections to the source terms, namely a spin magnetization ** M**, a spin polarization

**and a new contribution to the current density; see equation (3.17). These corrections appear when one considers an expansion in 1/**

*P**c*of the Dirac–Maxwell equations and are consistent with the Hamiltonian of equation (3.3), as was shown recently using a Lagrangian method [34,43]. Using equation (2.9), we can transpose their results to our formulation, which yields 3.16 3.17 3.18 3.19This mean-field approach could, in principle, be extended, in the spirit of density functional theory, to include exchange and correlation effects by adding suitable potentials and fields that are functionals of the electron density [44].

## 4. Semiclassical limit and spin Vlasov model

The form of equations (3.11) and (3.12) is particularly useful to study the semiclassical limit of the model. Indeed, we can easily expand , and as a power series of 4.1 4.2 4.3From these semiclassical expansions, we note that the velocity shift has a purely quantum origin, because the leading term in the expansion is of first order in . Therefore, it has no classical counterpart. In the case of and , the leading term in the expansion simply corresponds to the classical electric or magnetic field.

To zeroth order, the equations for *f*_{0} and *f*_{i} decouple, so that one can study the particle motion irrespective of the spin degrees of freedom. To first order in , equations (3.11) and (3.12) become
4.4
4.5where the factor is hidden in the definition of the Bohr magneton . The quantum corrections in equations (4.4) and (4.5) couple the orbital and the spin dynamics through the Zeeman and spin–orbit interactions. There are no quantum corrections to the orbital electronic dynamics because they would appear only at the second order in . For instance, the Darwin term would not appear in the above equations (even if we had included it in the full Wigner equations) because it is a correction of order to the orbital motion of the electron. In summary, equations (4.4) and (4.5) represent a semiclassical model where the orbital dynamics is classical (hence the familiar Lorentz force terms), whereas the spin is treated as a fully quantum variable (two-dimensional Hilbert space).

In equations (4.4) and (4.5), the Zeeman interaction gives two contributions: (i) the term *μ*_{B}**∇***B*_{i}⋅**∇**_{v}, which represents the force exerted on a magnetic dipole by an inhomogeneous magnetic field and is at the basis of Stern–Gerlach-type experiments, and (ii) the term ** f**×

**, which describes the precession of the spin around the magnetic field lines.**

*B*The spin–orbit interaction yields a correction to the magnetic field , which is the first-order correction in the non-relativistic limit of the Thomas precession [45,46]. The other terms are related to the spin–orbit correction of the velocity operator. Indeed, in the Heisenberg picture, the velocity operator is determined by the evolution equation of the position operator
4.6where we used the Hamiltonian defined in equation (3.3). The associated phase-space function is determined by the Weyl correspondence rule and reads as
4.7This is the phase-space function that can be used to calculate the average velocity or the charge current. Therefore, the particles are transported with a modified velocity. The term (** E**×

**∇**)

_{i}in equation (4.4) and (4.5) is a direct consequence of this effect, whereas the term [

**×(**

*E***×**

*B***∇**

_{v})]

_{i}corresponds to the same velocity correction in the Lorentz force

**×**

*v***.**

*B*The spin Vlasov equations (4.4) and (4.5) are correct to second order in 1/*c* and can be viewed as a generalization of those obtained in [31], where only the Zeeman interaction was included. An alternative form of these equations was obtained by Asenjo *et al.* [28] in the extended phase-space formalism.

The Maxwell equations (3.15), combined with the spin Vlasov equations (4.4) and (4.5), form a self-consistent model for the charge and the spin dynamics of a system of interacting particles. One can show that the following quantities are conserved during the time evolution:
4.8
4.9
4.10
4.11where ** D**=

*ϵ*

_{0}

**+**

*E***and**

*P***=**

*H***−**

*B**μ*

_{0}

**. The conserved quantities are the total mass**

*M**M*

_{tot}, the total linear momentum

*P*_{tot}(sum of the particles and fields momenta), the total energy

*E*

_{tot}(sum of the kinetic, Zeeman and the electromagnetic field energies), and the total angular momentum

*J*_{tot}(sum of the orbital, spin and electromagnetic field angular momenta).

For simulation purposes, the spin Vlasov equations (4.4) and (4.5) are much easier to solve numerically than the corresponding Wigner equations (3.11) and (3.12), mainly because the former are local in space while the latter are not. The Vlasov approximation is valid when quantum effects in the orbital dynamics are small. From equation (2.6), it appears that the semiclassical expansion is valid when , where *L*_{0} and *v*_{0} are typical length and velocity scales. For a degenerate electron gas with density *n*, the typical velocity is the Fermi speed . Inserting into the previous inequality, we obtain the validity condition *L*_{0}*n*^{1/3}≫1, which means that the typical length scale must be larger than the interparticle distance *d*=*n*^{−1/3}. For this reason, the semiclassical limit is also referred to as the long-wavelength approximation. All in all, the above spin Vlasov equations constitute a valuable tool to simulate the charge and spin dynamics in condensed matter systems, particularly semiconductor and metallic nano-objects.

## 5. Hydrodynamic model with spin–orbit coupling

In this section, starting from equations (4.4) and (4.5), we derive the hydrodynamic evolution equations by taking velocity moments of the phase-space distribution functions. In addition to the particle density and spin polarization (equations (3.8) and (3.9)), we define the following macroscopic quantities:
5.1
5.2
5.3
5.4
5.5where we separated the mean fluid velocity ** u** from the velocity fluctuations

**≡**

*w***−**

*v***. Here,**

*u**P*

_{ij}and

*Q*

_{ijk}are respectively the pressure and the generalized energy flux tensors. They coincide with the analogous definitions for spinless fluids with probability distribution function

*f*

_{0}. The spin–velocity tensor represents the fluid current along the

*i*th direction of the

*α*th spin polarization vector, whereas

*Π*

_{ijα}represents the corresponding spin-pressure tensor.

^{2}The evolution equations for the above fluid quantities are obtained by the straightforward integration of equations (4.4) and (4.5) with respect to the velocity variable. One obtains 5.6 5.7 5.8 5.9 5.10where we introduced a new average velocity and a new spin current 5.11

The above corrections reflect the modification of the velocity owing to the spin–orbit coupling. Indeed, the average velocity can be immediately obtained from the velocity phase-space function, equation (4.7), yielding
5.12where *F* is the 2×2 distribution function defined in equation (3.6). The same holds for the spin current operator, which is defined as follows:
5.13where we symmetrized the operator so that it is Hermitian. Then the associated phase-space function,
5.14can be used to determine the spin current
5.15

As is always the case for hydrodynamic models, some further hypotheses are needed to close the above set of equations (5.6)–(5.10). A particularly interesting strategy, based on the maximum entropy method (MEP), was used in a previous work [31] to close the set of hydrodynamic equations in the presence of the sole Zeeman interaction. Unfortunately, when one adds the spin–orbit interaction, the MEP does not provide any conclusive analytical results (the difficulty arises from the fact that the spin–orbit interaction couples all the components of the velocity). However, an intuitive closure can be found by inspecting the evolution equation (5.10) for the pressure tensor. There, most spin-dependent terms cancel if we set
5.16The physical interpretation of the above equations is that the spin of a particle is simply transported along the mean fluid velocity. This is of course an approximation that amounts to neglecting some spin–velocity correlations [27]. With this assumption, equation (5.9) and the definition of the spin-pressure *Π*_{ijα} are no longer necessary. The system of fluid equations simplifies to
5.17
5.18
5.19
5.20In order to complete the closure procedure, one can proceed in the same way as is usually done for spinless fluids, for instance by supposing that the system is isotropic and adiabatic. The isotropy condition imposes that *P*_{ij}=(*P*/3)*δ*_{ij}, where *δ*_{ij} is the Kronecker delta, whereas the adiabaticity condition requires that the heat flux *Q*_{ijk} vanish. In that case, the pressure takes the usual form of the equation of state of an adiabatic system, i.e. *P*=*const*.×*n*^{(D+2)/D} (where *D* is the dimensionality of the system), which replaces equation (5.20). In summary, equations (5.17)–(5.19), together with the preceding expression for the pressure, constitute a closed system of hydrodynamic equations with spin–orbit effects.

## 6. Conclusion

Phase-space methods can be applied to condensed matter and nanophysics to model the electron dynamics in either the quantum or the semiclassical regime. Several studies were performed in the past but neglected spin effects [8,44,47]. In this paper, we showed that phase-space methods can be conveniently generalized to include the spin dynamics at different orders. In an earlier work, we had developed a phase-space model that includes the lowest order spin term (the Zeeman effect), but neglects all relativistic corrections (spin–orbit coupling, Darwin term, mass correction, …). Here, we considered the case where both the Zeeman and the spin–orbit interaction are present (other relativistic corrections could be added with relative ease). The spin–orbit interaction plays an important role, for instance, in ultrafast spectroscopy experiments on magnetic nano-objects, where the electron spin is known to interact with the incident laser field and with the self-consistent field generated by the electron gas.

We first derived a four-component Wigner equation to describe the quantum dynamics of a system of spin- particles. These equations, together with the appropriate Maxwell equations, form a fully quantum self-consistent model to study the spin and charge dynamics in the mean-field approximation. This model is not limited to the linear response, but can deal with nonlinear effects, which are often important, particularly for large incident laser powers.

Next, using a semiclassical expansion to first order in , we obtained a four-component Vlasov equation. The orbital part of the motion is classical, i.e. the particles follow the classical phase-space trajectories, while the spin degrees of freedom are treated in a fully quantum fashion (two-dimensional Hilbert space). These spin Vlasov equations constitute a good approximation of the quantum electron dynamics for wavelengths larger than the typical inter-electron distance.

The corresponding hydrodynamic equations were derived by taking velocity moments of the phase-space distribution functions. The spin–orbit interaction modifies considerably our earlier hydrodynamic equations [31], where the only spin effect was the Zeeman interaction. We proposed a simple, intuitive closure for the hydrodynamic equations whereby the spin is simply transported along by the fluid velocity of the electrons.

The present models (Vlasov and hydrodynamic) constitute two valuable tools to study the intertwined spin and charge dynamics in condensed matter systems and nano-objects. The challenge now is to implement these models in performing numerical codes and to study the electron dynamics in realistic nanoscale systems that are relevant to current experiments.

## Authors' contributions

J.H. performed the analytical calculations and drafted the manuscript. P.-A.H. and G.M. conceived the study, checked the calculations, and contributed to the writing of the manuscript.

## Competing interests

The authors declare no competing interests.

## Funding

We thank the Agence Nationale de la Recherche, project Labex ‘Nanostructures in Interaction with their Environment’, for financial support.

## Footnotes

One contribution of 11 to a theme issue ‘Theoretical and computational studies of non-equilibrium and non-statistical dynamics in the gas phase, in the condensed phase and at interfaces’.

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3685624.

↵1 For the demonstration, we use the following properties of the Wigner operator:

*δ*(−*r*′).*r*↵2 Strictly speaking, a pressure tensor should be defined in terms of the velocity fluctuations

*w*_{i}*w*_{j}, but this would unduly complicate the notation. Thus, we stick to the above definition of*Π*_{ijα}while still using the term ‘pressure’ for this quantity.

- Accepted October 26, 2016.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.