# Electrical Conductivity of Hydrogen Plasmas: Low-density Benchmarks and
Virial Expansion Including $e-e$ Collisions

G. Röpke<EMAIL_ADDRESS>Institut für Physik, Universität Rostock,
D-18051 Rostock, Germany

###### Abstract

An improved virial expansion for the low-density limit of the electrical
conductivity $\sigma(T,n)$ of hydrogen as the simplest ionic plasma is
presented. Quantum statistical methods provide exact values for the lowest
virial coefficients, which serve as a benchmark for analytical approaches to
electrical conductivity as well as for numerical results from density
functional theory based molecular dynamics simulations (DFT-MD) or path-
integral Monte Carlo (PIMC) simulations. The correction factor introduced by
Reinholz et al., Phys. Rev. E 91, 043105 (2015) is applied to describe the
inclusion of electron-electron collisions in DFT based calculations of
transport coefficients. As a benchmark, the first virial coefficient is
correctly described with this approach. The value of the second virial
coefficient is discussed, questions about its value according to DFT-MD
simulations are addressed.

##  I Introduction

Coulomb plasmas such as the hydrogen plasma are strongly interacting many-
particle systems in which the formation of bound states (atoms and their
ionisation stages) is possible. Depending on the parameter values temperature
$T$ and particle number density $n_{e}=n_{i}=n$ (we consider a charge-neutral
hydrogen plasma consisting of electrons $e$, mass $m_{e}$, and protons $i$,
mass $m_{i}$), the electrons can be bound or degenerate, so that a quantum
statistical approach is required. Other ionic plasmas with atomic nuclei,
charge number $Z$, can be treated in a similar way, but are not considered
here. Instead of the variables $T$ and $n$, which characterise the state of
the plasma, the dimensionless plasma parameter

$\Gamma=\frac{e^{2}}{4\pi\epsilon_{0}k_{\rm
B}T}\left(\frac{4\pi}{3}n\right)^{1/3}$ (1)

is usually introduced as the ratio of potential to kinetic energy in the non-
degenerate case, and the electron degeneracy parameter

$\Theta=\frac{2m_{e}k_{\rm B}T}{\hbar^{2}}(3\pi^{2}n)^{-2/3}.$ (2)

The dc conductivity $\sigma_{\rm dc}(T,n)$ of the hydrogen plasma was
calculated in the low-density limit in the seminal work by Spitzer and Härm
Spitzer53 using kinetic theory. The Boltzmann equation valid in the low-
density limit and the corresponding Fokker-Planck equation were solved, with
the result

$\sigma^{\rm Spitzer}(T,n)=\frac{(k_{\rm
B}T)^{3/2}(4\pi\epsilon_{0})^{2}}{m_{e}^{1/2}e^{2}}\frac{1}{0.84624\ln(\Theta/\Gamma)}.$
(3)

In order to treat arbitrary densities, the linear response theory was worked
out, which considers the plasma under the influence of a weak external field.
Properties of the plasma are determined by correlation functions. According to
linear response theory, the dc conductivity is related to equilibrium
fluctuations of the charge current (which is mainly the electron current due
to the large mass ratio $m_{i}/m_{e}\gg 1$), resulting in the Kubo formula
Kubo66

$\sigma(T,n)=\frac{e^{2}\beta}{3m_{e}^{2}\Omega}\int_{-\infty}^{0}dte^{-izt}\int_{0}^{1}d\lambda\,{\rm
Tr}[{\bf P\cdot P}(t+i\hbar\beta\lambda)\rho_{0}]$ (4)

where $\beta=1/k_{\rm B}T$, $\Omega$ is the volume, $\bf P$ the total momentum
of electrons, the $t$ dependence is according to the Heisenberg picture

$O(t)=e^{iHt/\hbar}Oe^{-iHt/\hbar},$ (5)

and

$\rho_{0}=e^{-\beta(H-\sum_{c}\mu_{c}N_{c})}/{\rm
Tr}e^{-\beta(H-\sum_{c}\mu_{c}N_{c})}$ (6)

is the equilibrium statistical operator. Finally the limit $z=i\eta\to 0$ has
to be taken.

Equilibrium correlation functions can be evaluated using the methods of
quantum statistics such as the Green’s function method FW ; GRBuch . This is a
perturbative approach in which an expansion with respect to the Coulomb
interaction is performed Roep88 ; RR89 ; Heidi . By applying partial
summations describing, e.g., quasiparticles, screening and the formation of
bound states, useful results are obtained for a wide range of plasma
parameters. The application to electrical conductivity was also shown in
Rcond2018 . For the relation to kinetic theory see Reinholz12 . However, exact
results are only obtained in special limiting cases. In the following, we
discuss the virial expansion of the inverse conductivity of the hydrogen
plasma, which contains these exact results. They can be employed as benchmarks
for any approach.

In order to avoid perturbation theory for the electron-ion interaction, the
method of DFT-MD simulations was worked out Desjarlais . The ion dynamics is
treated classically using molecular dynamics simulations, while the electrons
are treated as quantum particles, with the wave equations being solved using
the DFT formalism. The evaluation of the correlation function (4) is performed
numerically using the Kubo-Greenwood formula. However, the treatment of the
electron-electron ($e-e$) interaction is approximated by a suitable choice of
the exchange-correlation functional $E_{xc}$ Perdew1996 . This problem is also
a subject of the present work.

The rigorous treatment of the $e-e$ interaction is possible with the help of
path-integral Monte Carlo (PIMC) simulations, see Mil01 ; Bo20 ; Boehme22 and
references therein. However, the calculations are very complex and the current
shortcomings of this approach include the relatively small number of particles
and the sign problem for fermions.

In this work, we focus on two aspects: the virial expansion of the resistivity
of hydrogen plasmas in the low-density region, and how DFT-MD simulations can
be improved to describe the contribution of electron-electron collisions to
plasma conductivity. In Sec. II we discuss the virial expansion of the
resistivity. The DFT-MD simulations of plasma conductivity are presented in
Sec. III, and ways to account for $e-e$ collisions are discussed in Sec. IV.
Finally, we draw conclusions in Sec. V.

##  II Benchmarks from the virial expansion of the resistivity

.

###  II.1 The standard virial expansion

For the sake of simplicity, we consider hydrogen plasmas. The results can be
transferred to other plasmas with ionic charge $Ze$, since the low-density
limit considered here is dominated by the Coulomb force and short-range
interactions contribute to higher orders of the virial expansion of the
resistivity.

To investigate the conductivity of the hydrogen plasma, we introduce the
dimensionless conductivity $\sigma^{*}(T,n)$ according to

$\displaystyle\sigma(T,n)$ $\displaystyle=$ $\displaystyle\frac{(k_{\rm
B}T)^{3/2}(4\pi\epsilon_{0})^{2}}{m_{e}^{1/2}e^{2}}\;\sigma^{*}(T,n)=\frac{32405.4}{\Omega{\rm
m}}T_{\rm eV}^{3/2}\sigma^{*}(T,n)\,,$ (7)

where $T_{\rm eV}=k_{\rm B}T/{\rm eV}=11604.6\,T/{\rm K}$ denotes the
temperature measured in units eV.

Using the generalised linear response theory, a virial expansion was proposed
for the dimensionless resistivity $\rho^{*}(T,n)=1/\sigma^{*}(T,n)$ Roep88 ;
RR89 ; Redmer97

$\displaystyle\rho^{*}(T,n)=\rho_{1}(T)\ln\left(\frac{\Theta}{\Gamma}\right)+\rho_{2}(T)+{\cal
O}\left[n^{1/2}\ln\left(\frac{\Theta}{\Gamma}\right)\right]$ (8)

with ($T_{\rm Ha}=27.21137\,T_{\rm eV},n_{\rm Bohr}=n\,a_{\rm Bohr}^{3}$)

$\frac{\Theta}{\Gamma}=\frac{1}{(36\pi^{5})^{1/3}}\frac{2m_{e}}{\hbar^{2}}\frac{4\pi\epsilon_{0}}{e^{2}}\frac{(k_{\rm
B}T)^{2}}{n}=\frac{T_{\rm Ha}^{2}}{4n_{\rm
Bohr}}(96\pi^{5})^{-1/3}\propto\frac{T^{2}}{n}\,.$ (9)

It describes the well-known logarithmic density dependence, the Coulomb
logarithm, which is caused by the long-range Coulomb interaction.

According to Spitzer and Härm Spitzer53 , who calculated the value
$\gamma_{E}=0.5816$, the first virial coefficient

$\rho_{1}(T)=\rho^{\rm
Spitzer}_{1}=\frac{\pi^{3/2}}{2^{7/2}\gamma_{E}}=0.84624$ (10)

is not dependent on $T$. For the second virial coefficient, the high-
temperature limit follows from the Quantum Lenard-Balescu (QLB) approach RR89
; RRT89 ; Redmer97 ; Karachtanov16

$\lim_{T\to\infty}\rho_{2}(T)=\rho_{2}^{\rm QLB}=0.4917~{}.$ (11)

In this limit the conductivity of hydrogen plasmas is approximated as

$\sigma(T,n)=\frac{32405.4\,T_{\rm
eV}^{3/2}}{0.84624\ln(\Theta/\Gamma)+0.4917}\frac{1}{\Omega{\rm m}}.$ (12)

The temperature dependence of $\rho_{2}(T)$ is not exactly known, an
approximation can be obtained from the interpolation formula RRT89 ; Esser03

$\rho_{2}(T)\approx 0.4917+0.846\ln\left[\frac{1+8.492/T_{\rm
eV}}{1+25.83/T_{\rm eV}+167.2/T_{\rm eV}^{2}}\right]\,.$ (13)

It is a challenge for numerical simulations such as the PIMC simulation to
improve this approximation for the temperature-dependent second virial
coefficient.

The virial expansion serves as a benchmark for the theoretical and
experimental determination of the plasma conductivity in the low-density
range. To extract the virial coefficients from data for $\sigma(T,n)$, a
virial plot was proposed in RSRB21 . To determine $\rho_{1}(T)$ and
$\rho_{2}(T)$, we draw

$\rho_{1}^{\rm
eff}(T,n)=\rho^{*}(T,n)\frac{1}{\ln(\Theta/\Gamma)}=\frac{1}{\sigma^{*}(T,n)}\frac{1}{\ln(\Theta/\Gamma)}$
(14)

as a function of $x=1/\ln(\Theta/\Gamma)$. Since

$\rho_{1}^{\rm eff}(T,n)=\rho_{1}(T)+\rho_{2}(T)x+{\cal O}[n^{1/2}]$ (15)

we expect a linear relationship in the range of low densities, where the
ordinate at $x=0$ is $\rho_{1}$ and the slope is $\rho_{2}$. To perform the
extrapolation to $x=0$, the variable $x$ must be sufficiently small for linear
behavior to be observed.

This method was applied to several approaches for determining the electrical
conductivity of hydrogen plasmas, in particular to the data presented in
Grabowski . As shown in RSRB21 , several semi-empirical expressions for plasma
conductivity fail to fulfill this benchmark.

The next term in the virial expansion (8), which is of the order
$n^{1/2}\ln(\Theta/\Gamma)$, is called the Debye-Onsager relaxation effect. An
approximate expression for this third virial term was given in Ref. Roep88 .
This higher order virial coefficient can be extracted by analysing the
effective second virial coefficient

$\rho_{2}^{\rm
eff}(T,n)=\rho^{*}(T,n)-\rho_{1}\ln(\Theta/\Gamma)=\rho_{2}(T)+\rho_{3}(T)n^{1/2}\ln(\Theta/\Gamma)+{\cal
O}\left(n^{1/2}\right)$ (16)

as a function of $x_{1}=n^{1/2}\ln(\Theta/\Gamma)$. The corresponding virial
plot gives $\rho_{2}(T)$ as the value of $\rho_{2}^{\rm eff}(T,n)$ at
$x_{1}=0$, while $\rho_{3}(T)$ follows from the slope at $x_{1}=0$. The third
virial term in the approximation of Ref. Roep88 is only a small correction and
is not considered here. To extract it from a virial plot, accurate simulations
or measurements of the plasma conductivity are required.

Table 1: Electrical conductivity $\sigma(T,n)$ of hydrogen plasmas from virial expansion, Eq. (12), $\sigma_{\Lambda}$ from Eq. (22), $\sigma_{\hat{\Lambda}}$ from Eq. (24). The density $\rho=1$ g cm-3 corresponds to the particle number density $n=5.98\times 10^{23}$ cm-3, the density $\rho=1.67$ g cm-3 to $n=1\times 10^{24}$ cm-3, and $\rho=10$ g cm-3 to $n=5.98\times 10^{24}$ cm-3. $\rho$ (g cm-3) |  $T$ (eV) | $\Gamma$ | $\Theta$ | $1/\ln(\Theta/\Gamma)$ |  $\sigma$ ($\Omega^{-1}$ cm-1) |  $\sigma_{\Lambda}$ ($\Omega^{-1}$ cm-1) | $1/\Lambda$ |  $\sigma_{\hat{\Lambda}}$ ($\Omega^{-1}$ cm-1) | $1/{\hat{\Lambda}}$
---|---|---|---|---|---|---|---|---|---
1 | 2000 | 0.0097773 | 77.284 | 0.1114 | 3.58498$\times 10^{6}$ | 3.58497$\times 10^{6}$ | 0.1003 | 3.58496$\times 10^{6}$ | 0.1046
1 | 1000 | 0.0195546 | 38.642 | 0.1318 | 1.48255$\times 10^{6}$ | 1.48252$\times 10^{6}$ | 0.1164 | 1.48251$\times 10^{6}$ | 0.1224
1 | 700 | 0.0279351 | 27.049 | 0.1454 | 951335 | 951229 | 0.127 | 951281 | 0.1341
1 | 400 | 0.0488865 | 15.457 | 0.1737 | 483512 | 483446 | 0.1481 | 483412 | 0.1578
1 | 200 | 0.097773 | 7.7284 | 0.2288 | 218811 | 218659 | 0.1862 | 218582 | 0.2018
1 | 100 | 0.195546 | 3.8642 | 0.3352 | 107445 | 107035 | 0.2504 | 106826 | 0.2789
1 | 70 | 0.279351 | 2.7049 | 0.4405 | 78668.2 | 77910.1 | 0.3035 | 77531 | 0.3456
1.67 | 20000 | 0.00116 | 549.05 | 0.07653 | 7.93757$\times 10^{7}$ | 7.93757$\times 10^{7}$ | 0.0711 | 7.93757$\times 10^{7}$ | 0.07327
1.67 | 10000 | 0.00232 | 274.53 | 0.08561 | 3.12362$\times 10^{7}$ | 3.12362$\times 10^{7}$ | 0.07888 | 3.12362$\times 10^{7}$ | 0.08155
1.67 | 1000 | 0.0232 | 27.453 | 0.1413 | 1.58184$\times 10^{6}$ | 1.58179$\times 10^{6}$ | 0.1239 | 1.58176$\times 10^{6}$ | 0.1306
1.67 | 100 | 0.232 | 2.7453 | 0.4047 | 125498 | 124570 | 0.2865 | 124103 | 0.324
10 | 1000 | 0.0421291 | 8.3252 | 0.1892 | 2.06434$\times 10^{6}$ | 2.06386$\times 10^{6}$ | 0.1591 | 2.06361$\times 10^{6}$ | 0.1704
10 | 100 | 0.421291 | 0.83252 | 1.468 | 303434 | 274960 | 0.5529 | 263134 | 0.687
100 | 1000 | 0.0907644 | 1.7936 | 0.3352 | 3.39771$\times 10^{6}$ | 3.38473$\times 10^{6}$ | 0.2504 | 3.37813$\times 10^{6}$ | 0.2789

###  II.2 Generalized virial expansions

The occurrence of a logarithmic term $\propto\ln(1/n)$ in the virial expansion
of the resistivity, Eq. (8), requires some discussion. We need to define a
reference density $n^{\rm ref}(T)$ to make the logarithm dimensionless,
$\ln(1/n)\to\ln(n^{\rm ref}(T)/n)$. In Eq. (8), $n^{\rm
ref}(T)=n\Theta/\Gamma$ was taken. The choice of the reference density
influences the higher virial coefficients. In principle, $\hat{n}^{\rm
ref}(T)$ can be chosen such that the next virial coefficient disappears,
$\hat{\rho}_{2}(T)=0$, as shown below.

Another problem is that the logarithm can become zero or negative, $\ln(n^{\rm
ref}/n)\leq 0$ for $n^{\rm ref}/n\leq 1$, but the conductivity or resistivity
are positively defined quantities. Therefore, it is of interest to perform a
modified virial expansion using the Coulomb logarithm $\Lambda(T,n)$ from
kinetic theory (for references see French22 ),

$\displaystyle\rho^{*}(T,n)=\rho_{1}\Lambda(T,n)+\tilde{\rho}_{2}(T)+{\cal
O}[n^{1/2}\Lambda(T,n)],$ (17)

$\Lambda(T,n)=\ln(1+b)-\frac{b}{1+b},$ (18)

where

$b(T,n)=\frac{3(k_{\rm B}T)^{2}}{\pi
n}\,\frac{4\pi\epsilon_{0}m}{e^{2}\hbar^{2}}=\frac{n_{\Lambda}(T)}{n}=\frac{3^{4/3}\pi^{2/3}}{2^{1/3}}\frac{\Theta}{\Gamma}=\frac{3}{\pi}\frac{T_{\rm
Ha}^{2}}{n_{\rm Bohr}}$ (19)

is the Born parameter $\Theta/\Gamma$ (up to a factor). We have in the low-
density limit

$\lim_{n\to 0}\Lambda(T,n)=\ln[b(T,n)]-1.$ (20)

The first virial coefficient $\rho_{1}(T)$ remains unchanged, but the higher
virial coefficients are modified if $\ln(n^{\rm ref}/n)$ with $n^{\rm
ref}=n\Theta/\Gamma$ is replaced by $\Lambda(T,n)$ in all terms of the virial
expansion. In particular, from Eq. (17) we obtain

$\lim_{T\to\infty}\tilde{\rho}_{2}(T)=\lim_{T\to\infty}\rho_{2}(T)+\rho_{1}\left[1-\frac{1}{3}\ln(81\pi^{2}/2)\right]=0.4917+0.84624[1-\frac{1}{3}\ln(81\pi^{2}/2)]=-0.351938.$
(21)

The conductivity follows as

$\sigma_{\Lambda}(T,n)=\frac{32405.4\,T_{\rm
eV}^{3/2}}{0.84624\,\Lambda-0.351938}\frac{1}{\Omega{\rm m}}$ (22)

The value of $b(T,n)$, Eq. (19), is related to the collision integral
evaluated for a screened potential, see GRBuch , but changes when the static
screening is replaced by dynamic screening. We can scale $b(T,n)$ by a factor
$y$ without changing the analytical structure of the virial expansion. It is
possible to set the second virial coefficient to zero at $T\to\infty$ if the
parameter $b$ is replaced by $\hat{b}=b/y$ with $y=1.5158$. Then we have

$\rho^{*}(T,n)=\rho_{1}\hat{\Lambda}(T,n)+\hat{\rho}_{2}(T)+{\cal
O}[n^{1/2}\hat{\Lambda}(T,n)],$ (23)

where $\hat{\Lambda}(T,n)$ follows from Eq. (18) replacing $b$ by $\hat{b}$.
We have $\lim_{T\to\infty}\hat{\rho}_{2}(T)=0$, and in the high-temperature
limit

$\sigma_{\hat{\Lambda}}(T,n)=\frac{32405.4\,T_{\rm
eV}^{3/2}}{0.84624\,\hat{\Lambda}}\frac{1}{\Omega{\rm m}}.$ (24)

Note that a $T$-dependent $y(T)$ can be introduced to compensate for the
second virial coefficient for all $T$. These modifications of the virial
expansion lead to changes in the higher virial coefficients. It is a question
to the higher order expansions whether such modifications would result in a
better convergence.

As example, the values for the conductivity for specific values of $T,n$ are
shown in Tab. 1. The values for low densities and high $T$ are considered
where the parameter $x,\tilde{x}=1/\Lambda$ or $\hat{x}=1/{\hat{\Lambda}}$ is
small, so that the virial expansion with error bars below few percent is
valid. The largest discrepancy occurs for the relatively high density 10 g/cm3
and the low temperature 100 eV. The method of virial plots can also be used to
extract the virial coefficients $\rho_{1}(T),\rho_{2}(T),\tilde{\rho}_{2}(T)$
or $\hat{\rho}_{2}(T)$ from experimental data or numerical simulations. An
example is discussed in Sec. III below, see also RSRB21 .

The method of virial plots has also been applied to thermodynamic variables
describing the equation of state of plasmas, see RDVBB24 . There, higher-order
virial coefficients and generalized virial expansions are discussed. Higher-
order virial coefficients are of interest to describe degeneracy as known from
the Lee-More model, see, e.g., French22 .

##  III The Kubo-Greenwood formula and DFT-MD simulations

The DFT-MD approach has been successfully applied to calculate the
thermodynamic properties of complex materials in a wide range of $T$ and $n$.
In the Born-Oppenheimer approximation, the ions are treated as classical
particles whose time-dependent configuration is described by molecular
dynamics simulations. They act like an external potential on the electron
subsystem. Together with the kinetic energy and Coulomb’s $e-e$ interaction, a
many-electron Schrödinger equation is formulated, which determines the many-
body wave functions. In density functional theory, the correlated many-
electron wave function is replaced by the antisymmetrized product of optimized
single-electron states, the Kohn-Sham orbitals. For the corresponding Kohn-
Sham equations, the Hamiltonian contains the kinetic energy of a non-
interacting reference system, the mean-field Coulomb interaction energy, and
an exchange-correlation energy $E_{xc}$ that accounts to some approximation
for the effects of $e-e$ correlations and antisymmetrization. A simple
approximation uses expressions for $E_{xc}$ that are derived for the uniform
electron gas. For frequency-dependent electrical conductivity, see (4) for
$z=\omega+i\eta$, the Kubo-Greenwood formula Kubo66 ; Greenwood

$\operatorname{Re}\left[\sigma(\omega)\right]=\frac{2\pi
e^{2}}{3m_{e}^{2}\omega\Omega}\sum_{k}w_{k}\sum_{j=1}^{N}\sum_{i=1}^{N}\sum_{\alpha=1}^{3}\big{[}f(\epsilon_{j,k})-f(\epsilon_{i,k})\big{]}|\langle{\Psi_{j,k}|\hat{p}_{\alpha}|\Psi_{i,k}}\rangle|^{2}\delta(\epsilon_{i,k}-\epsilon_{j,k}-\hbar\omega)$
(25)

was used to calculate the frequency-dependent dynamic electrical conductivity
$\sigma(\omega)$ in the long-wavelength limit RSRB21 ; Desjarlais02 ;
Mazevet05 ; Holst11 ; French2017 ; Gajdos2006 . Kohn-Sham wave functions
$\Psi_{i,k}$ from density functional theory calculations are used to calculate
the transition matrix elements of the electron momentum operator
$\hat{p}_{\alpha}$, $\alpha=\\{x,y,z\\}$. The Fermi-Dirac distribution
$f(\epsilon)$ accounts for the average occupation at energy $\epsilon$, and
the summation over momentum space $k$ contains the $k$-point weights $w_{k}$.

Due to the finite size of the simulation box, the $\delta$-function in
equation (25) must be approximated by a Gaussian with finite width, which also
prevents the direct calculation of the DC conductivity at $\omega=0$.
Therefore, the dynamic conductivity is extrapolated to the limit
$\omega\rightarrow 0$ using a Drude fit. The extrapolation procedure to
$\omega=0$ can be improved by using a frequency-dependent collision frequency
Heidi .

One of the main shortcomings of the DFT-MD approach is that the many-particle
interaction is replaced by a mean-field potential. When using product wave
functions for the many-electron system, correlations are excluded. The
exchange-correlation energy density functional reflects the Coulomb
interaction in a certain approximation, e.g., as it exists in the homogeneous
electron gas. However, it becomes problematic in the low-density limit, where
correlations, especially $e-e$ collisions, are important.

DFT-MD simulations have been successfully used to calculate the transport
properties in the warm-dense matter region for various materials. The
transition to condensed matter, the liquid or solid state, is adequately
described, in particular for the degenerate electron system. The question
whether the contribution of $e-e$ collisions to the conductivity is correctly
described also in the low-density range has long been the subject of
controversial debate. The $e-e$ Coulomb interaction is included in the
exchange-correlation energy of the DFT functional so that it was argued that
the ”ab initio” approach also describes the Coulomb interaction of the
electrons Desjarlais ; Starrett20 . However, it was pointed out Heidi that
DFT-MD is not able to reproduce the Spitzer limiting value of plasma
conductivity, what leads to significant deviations in the low-density range.
The correct behavior of the plasma conductivity at low density was determined
using generalized linear response theory (gLRT).

Figure 1: Reduced resistivity $\rho_{1}^{\rm eff}(T,x)$, Eq. (14), for
hydrogen plasma as a function of $x=1/\ln(\Theta/\Gamma)$: DFT-MD simulations
from Ref. RSRB21 and the generalized linear response result of Redmer and
Karakhtanov RR89 ; RRT89 ; Redmer97 ; Karachtanov16 in the high-temperature
limit. $\rho^{\rm Spitzer}_{1}=0.846$, Eq. (10), and $\rho^{\rm
Lorentz}_{1}=0.492$, Eq. (29), are defined in the text. The dot-dashed lines
are linear relations connecting the Spitzer or Lorentz limit with the
respective reduced resistivity value for the smallest $x$ value.

The method of virial plots was used to show that the DFT-MD simulations fail
to describe the correct low density behavior of the electrical conductivity of
the hydrogen plasma RSRB21 ; SCCS22 , see Fig. 1. Instead of the correct
virial expansion (8) for the resistivity of the hydrogen plasma,

$\rho^{*}(T,n)\approx 0.84624\ln\left(\frac{\Theta}{\Gamma}\right)+0.4917$
(26)

at high temperatures, the DFT-MD simulations are extrapolated as

$\rho^{*,\rm DFT-MD}(T,n)\approx
0.492\ln\left(\frac{\Theta}{\Gamma}\right)+0.9886.$ (27)

This is the virial expansion for the Lorentz plasma, a model plasma where the
Hamiltonian contains only the electron-ion interaction, the electron-electron
interaction is replaced by the interaction within the homogeneously charged
jellium to compensate the Coulomb divergency of the ion system,

$H^{\rm
Lorentz}=\sum_{c=i,e,k}\frac{p_{c,k}^{2}}{2m_{c}}+\frac{1}{2}\sum_{k\neq
l}V^{\prime}_{ii}({\bf r}_{i,l}-{\bf r}_{i,k})+\sum_{j,k}V^{\prime}_{ie}({\bf
r}_{e,l}-{\bf r}_{i,k})$ (28)

($V^{\prime}$ means the exclusion of the Hartree term in the Coulombic
pseudopotential interaction (screened by electrons) due to the charge
neutrality by the negative background) which has the first virial coefficient

$\rho_{1}^{\rm Lorentz}=\frac{1}{16}(2\pi^{3})^{1/2}=0.492126\,.$ (29)

The second virial coefficient $\lim_{T\to\infty}\rho_{2}^{\rm DFT-
MD}(T)=0.9886$ is determined by the Born approximation, but includes the
screening. In SCCS22 it was discussed that in addition to the screening of the
ion subsystem, which is given by the pair distribution function resulting from
the MD simulations, the screening by the electrons according to the exchange-
correlation energy density functional is also included. However, dynamical
screening is not described by $E_{xc}$, in contrast to the quantum Lenard-
Balescu equation.

We would like to mention that in the case of thermal conductivity it has been
shown that the contribution of $e-e$ collisions is not taken into account in
DFT-MD simulations Starrett20 ; Desjarlais ; French22 and results in an
additional term. Other approaches such as generalized linear response theory
can be used to solve this problem.

The DFT-MD simulations have been used very successfully to calculate the
transport properties of warm and dense matter (WDM). Electron-electron
collisions are not treated adequately. For degenerate systems, $e-e$
collisions are suppressed by Pauli blocking, so that their contribution can be
neglected and the results of DFT-MD are applicable. However, the systematic
inclusion of $e-e$ collisions is indispensable for the description of plasmas
at $\Theta>1$.

##  IV The account of electron-electron collisions in numerical simulations of
plasma conductivity

###  IV.1 The correction factor for Lorentz plasmas

We refer to a Lorentz plasma as a model system of electrons and ions, in which
the interaction of electrons with ions is considered, see Eq. (28). The
electron-electron interaction is only considered in the mean-field
approximation, so that the electron-ion interaction is screened, but the
effect of electron-electron collisions on the transport coefficients is
neglected. For this Lorentz model plasma, the solution of the linearized
Boltzmann equation is found using the relaxation time approach. The relaxation
time ansatz is possible for elastic collisions, which is the case for
$\lim(m_{i}/m_{e})\to\infty$, but not for $e-e$ collisions.

For the Lorentz model plasma, the electrical conductivity $\sigma^{\rm
Lorentz}(T,n)$ can be calculated beyond perturbation theory. The effective
interaction of electrons with the ion subsystem is treated, for example, using
the density functional theory. Similar approaches are given by the average
atom model or the Lee-More model. These approaches are applicable in the high-
density, degenerate case ($\Theta\leq 1$), in which electron-electron
collisions are suppressed due to Pauli blocking. They fail to describe real
plasmas such as hydrogen plasmas in the low-density range, as $e-e$ collisions
are not taken into account.

We introduce a correction factor

$R_{ee}(T,n)=\frac{\sigma(T,n)}{\sigma^{\rm Lorentz}(T,n)}$ (30)

which relates the plasma conductivity $\sigma(T,n)$ to the conductivity of the
Lorentz model plasma. This correction factor $R_{ee}(T,\Theta)$ (we replace
the density parameter $n$ by $\Theta$) approaches 1 for $\Theta\leq 1$ and the
value $0.591\pi^{3/2}/2^{5/2}=0.582$ for $\Theta\gg 1$.

Exact expressions for the conductivity $\sigma(T,n)$ of plasmas are obtained
from the generalized linear response theory (gLRT) in the form of correlation
functions. The kinetic theory (KT) follows as a special case Reinholz12 .
However, the analytical evaluation of the correlation functions is only
possible in some approximations using Green function diagram techniques. For
example, the dynamic screening of the electrons is treated in the random phase
approximation (RPA), the ionic structure factor is also approximated, the Born
approximation can be improved by a T-matrix approach.

An expression for $R_{ee}(T,\Theta)$ was given in Ref. Heidi , where both
$\sigma(T,n)$ and $\sigma^{\rm Lorentz}(T,n)$ are treated by the same
approximation (gLRT with two moments, Born approximation, RPA). It is near to
1 in the region of degeneracy $\Theta\leq 1$ and approaches the limiting value
0.581 for $\Theta\to\infty$. It is calculated numerically by solving some
integrals, see Heidi . For arbitrary ionic charges $Z$ a smooth interpolation
formula was given, which reads for hydrogen ($Z=1$) in the low-density,
classical limit (see Eq. (C.20) of Ref. Heidi )

$\displaystyle\lim_{T\to\infty}R_{ee}(T,\Theta)=0.594456+0.575578\frac{1}{\ln\left(0.0033305\,\Theta^{3/2}\,T_{\rm
K}^{1/2}\right)},$ (31)

where $T_{\rm K}=11604.6\,T_{\rm eV}$ is the temperature measured in Kelvin.
In the context of the virial expansion, using $x=1/\ln(\Theta/\Gamma)$, we
obtain

$\lim_{T\to\infty}R_{ee}(T,x)=0.594456+0.575578\,x$ (32)

in lowest order of $x$.

Within the framework of a refined calculation, we must investigate the effects
of dynamical screening, ion-ion structure factor, and strong collisions, see
Refs. RR89 ; Redmer97 ; Reinholz05 ; Karachtanov13 . These effects are of
minor importance for the correction due to $e-e$ collisions, both in the high-
density limit where $\lim_{\Theta\to 0}R_{ee}(\Theta)=1$ and in the low-
density limit, where they only occur in higher orders of the virial expansion.
Although these corrections are small, we analyze them in the next section
IV.2.

###  IV.2 Example: DFT-MD simulations for hydrogen plasmas

As shown in Fig. 1, the DFT-MD simulations for the effective resistivity
$\rho_{1}^{\rm eff,DFT-MD}$ are on a straight line pointing to the Lorentz
limit, see Ref. RSRB21 . The corrected values $\rho_{1}^{\rm eff,ee-
corr}=\rho_{1}^{\rm eff,DFT-MD}/R_{ee}$ are shown in Tab. 2. The correction
factor $R_{ee}$, Eq. (31), leads to a considerable increase (up to a factor
$\approx 2$ in the low-density range) of the resistivity and solves the
drawback of the correct $\rho_{1}$ benchmark by construction. However, the
slope determining the second virial coefficient seems to remain too large
compared to the high-temperature limit (11) obtained from QLB or gLRT
approaches Redmer97 ; Karachtanov16 . While the slope of the DFT-MD
simulations is 0.9965, the slope of the $e-e$ corrected virial line is 0.732.
The slope of the virial line in Fig. 1 determines the second virial
coefficient $\rho_{2}(T)$. The $e-e$ corrected value is closer to the high-
temperature benchmark 0.4917, but remains different.

Table 2: Virial representation of the dc conductivity $\sigma$ and of $\rho_{1}^{\rm eff}(T,x)$ with $x=1/\ln(\Theta/\Gamma)$: the values for $\sigma^{\rm DFT-MD}$ and $\rho_{1}^{\rm eff,DFT-MD}$ result from DFT-MD simulations RSRB21 . The $e-e$ corrected values for $\sigma^{\rm ee-corr}$, Eqs. (30), (31), and the corresponding $\rho_{1}^{\rm eff,ee-corr}$, Eq. (14), are also shown. $n$ | $k_{B}T$ | $\Gamma$ | $\Theta$ | $1/\ln(\Theta/\Gamma)$ | $\sigma^{\rm DFT-MD}$ | $\rho_{1}^{\rm eff,DFT-MD}$ | $\sigma^{\rm ee-corr}$ | $\rho_{1}^{\rm eff,ee-corr}$
---|---|---|---|---|---|---|---|---
[g/cm3] | [eV] |  |  |  | [$10^{6}/\Omega$ m] |  | [$10^{6}/\Omega$ m] |
2 | 75 | 0.3285 | 1.8257 | 0.58302 | 11.44 | 1.073 | 10.034 | 1.2234
2 | 100 | 0.24637 | 2.4343 | 0.43657 | 15.26 | 0.9269 | 12.434 | 1.1375
3 | 100 | 0.28203 | 1.8577 | 0.53047 | 16.85 | 1.020 | 14.412 | 1.1925
3 | 150 | 0.18802 | 2.7866 | 0.37092 | 25.67 | 0.8603 | 20.157 | 1.0956
4 | 150 | 0.20694 | 2.3003 | 0.41522 | 27.39 | 0.9026 | 22.058 | 1.1208

The $e-e$ correction factor (31) was obtained from a two-moment approximation.
The inclusion of higher moments of the distribution function would change the
$e-e$ correction factor and the corresponding values of
$\rho_{1},\rho_{2}(T)$. These corrections are expected to be small, see the
good convergence with increasing number of moments reported in Ref. Redmer97 .

The virial coefficient $\rho_{2}(T)$ is determined by the screening of the
Coulomb potential. The Born approximation can be applied in the high
temperature range $T_{\rm Ha}\gg 1$ so that strong collisions leading to a
T-matrix approach are not relevant. A reason for the $\rho_{2}$ discrepancy
can also be found in the treatment of screening, see the analysis of the slope
parameter in SCCS22 . The electrons should be treated as dynamic screening, in
lowest approximation by the RPA expression for the polarization function, to
obtain the QLB benchmark (11). It is not clear to what extent the dynamic
screening by electrons is already implemented in the DFT-MD simulations by the
choice of the exchange-correlation part $E_{xc}$ of the energy-density
functional.

If $\sigma^{\rm Lorentz}(T,n)$ which is required for the calculation of
$R_{ee}$, Eq. (31), treats the dynamical screening of the electron-ion
interaction differently than the DFT-MD approach $\sigma^{\rm DFT-MD}(T,n)$,
the use of the correction factor is not consistent since the slope parameter
$\rho_{2}(T)$ is not corrected. With respect to the contribution of the ions
to the screening, the ion structure factor describes the static screening
which is well approximated in the DFT-MD simulations.

In conclusion, the correction factor (31) given in Heidi reproduces the
general behavior for $\Theta\leq 1$ on $R_{ee}\approx 1$ and is a reasonable
approximation to the correct Spitzer result in the low-density limit. Only a
two-moment approximation of the electron distribution function is considered
to calculate $R_{ee}$ (31). It can be improved by considering higher moments
of the electron distribution function so that the exact value of the Spitzer
benchmark is reproduced. With respect to the second virial coefficient, the
specific choice of the $e-i$ interaction, which defines the Lorentz model
plasma as a reference system, is discussed in the following Sec. IV.3.

###  IV.3 Systematic approaches

While electron-ion collisions are well understood to calculate the
conductivity of plasmas, the inclusion of electron-electron collisions is
notoriously difficult. One reason for this is that a relaxation time approach
is only valid for elastic collisions where the single-particle energy is
conserved, what holds to a good approximation for $e-i$ collisions, but not
for $e-e$ collisions. Several approaches have been worked out to implement
$e-e$ collisions. The kinetic theory can solve this problem starting from the
linearized Boltzmann equation. The perturbation of the single electron
distribution function in a weak electrical field results from the Kohler
variational principle. In dense plasmas, gLRT is used to express the
conductivity in terms of equilibrium correlation functions Spitzer53 ; Roep88
.

To calculate the conductivity of degenerate electrons in condensed matter, the
Ziman formula is known, which considers the force-force autocorrelation
function. Calculations of the collision frequency based on the solution of the
electron-ion scattering problem are improved by a renormalization factor
$r(\omega)$ to obtain the correct dynamic collision frequency of the plasma
Reinholz00 ; Reinholz05 ; Heidi . While the Ziman formula is frequently used
to describe degenerate Coulomb plasmas such as electrons in liquid metals or
solids, the renormalization factor allows to treat also non-degenerate
plasmas, especially in the low-density range. The method of including $e-e$
correlations by means of a renormalization factor $r(\omega)$ is equivalent to
using the correction factor $R_{ee}$ and requires analogous approximations.

An improved approach to treat dense plasmas is to start from the Lorentz model
plasma, in which the $e-i$ interaction is solved, and the electron
distribution function is rigorously obtained from the relaxation time ansatz
for any degree of degeneracy. Starting from the linearized Boltzmann equation,
thermoelectric transport coefficients for Lorentz model plasmas are
calculated, see the Lee-More model LeeMore . For a dense system, the averaged
atom model was worked out, which describes the average influence of the
surrounding plasma on the atom. The interaction of the electrons with the ion
system is better described by the DFT approach, and the transition to
condensed matter is well established. The electron-electron interaction is
only considered in mean-field approximation to introduce an effective
potential for the single-electron quasiparticle states.

The success of relaxation-time models for transport coefficients to treat the
electron-ion interaction has led to develop a generalization of these models
which includes $e-e$ collisions. For example, as in Refs. Starrett20 ;
Starrett20a , the relaxation-time model by Starrett Starrett17 includes Pauli
blocking and accounts for correlations using a mean- force scattering cross-
section for electron-ion collisions, but electron-electron collisions are
accounted for only through the correction formula for $R_{ee}$ proposed in
Ref. Heidi .

A possible solution for this particular problem of including $e-e$ collisions
in the determination of plasma conductivity can be the use of path-integral
Monte Carlo (PIMC) simulations to calculate $\sigma(T,n)$. In PIMC
simulations, $e-e$ correlations are treated in an exact way. This enables the
immediate calculation of $R_{ee}(T,\Theta)$ according to (30) and the
construction of consistent interpolation formulas. In particular, the
determination of virial coefficients would be of interest. Such PIMC
simulations are not yet possible. Highly accurate results are only available
for the uniform electron gas RDVBB24 . In addition, there is the unsolved sign
problem, and PIMC simulations are currently performed for too small numbers of
particles (a few tens). It is expected that significant results for the
conductivity of hydrogen plasmas from PIMC simulations will be available in
the near future. This will overcome the drawback of DFT-MD simulations where
$e-e$ correlations are not treated precisely.

We propose to use the virial expansion (8) for the determination of the
conductivity of the hydrogen plasma in order to analyze the consistency of
different theoretical approaches using the virial coefficients as a benchmark.
Having $\rho_{1},\rho_{2}(T),$ etc. available, we can find improved versions
of the correction factor $R_{ee}(T,n)$. Since these correction factors relate
the plasma conductivity $\sigma(T,n)$ to the conductivity of a Lorentz model
plasma $\sigma^{\rm Lorentz}(T,n)$, see Eq. (30), these correction factors
depend on the definition of the Lorentz model Hamiltonian. They depend in
particular on the definition of the effective electron-ion potential, which
includes dynamic or static screening by electrons, or on the particular choice
of the exchange-correlation energy $E_{xc}$ in DFT-MD calculations. If the
determination of the conductivity and virial expansion for the Lorentz model
plasma is known, a corresponding correction factor can be defined.

In particular, we consider the correction factor

$R_{ee}^{\rm DFT-MD}(T,n)=\frac{\sigma(T,n)}{\sigma^{\rm{}^{\rm DFT-
MD}}(T,n)}$ (33)

and use the virial expansions (26), (27). The virial expansion for the high-
temperature limit results in

$\lim_{T\to\infty}R_{ee}^{\rm DFT-MD}(T,x)=0.5814+0.8304\,x$ (34)

in contrast to (32) which describes the relationship between the conductivity
of a model plasma and a Lorentz model plasma with the same form of screened
interaction. The expression (34) is the correction factor required to obtain
the conductivity of hydrogen plasmas in the low-density, high-temperature
limit after performing DFT-MD simulations. Compared to (32), which is
determined in the framework of a two-moment approximation, the first number
0.5814 results from the consideration of higher moments of the distribution
function. The second number 0.8304 results from the treatment of the screening
in the special version of the DFT approach and cannot be analyzed in this
work. The advantage of the correction factor $R_{ee}(T,\Theta)$ is that its
behavior is known in both limits: for $\Theta\leq 1$ it is close to 1, and for
$\Theta\gg 1$ it is determined by the virial expansion. Note that the
correction factor $R_{ee}(T,n,Z)$ in Heidi is given for any ionic charge $Z$,
so that it can also be applied to other substances.

##  V Conclusions

We propose an exact virial expansion (8) for the plasma conductivity to
analyze the consistency of theoretical approaches. The virial coefficients
serve as a benchmark for different approaches and can be better visualized
using virial plots. For example, various analytical calculations of the dc
conductivity $\sigma(T,n)$, which were presented in Ref. Grabowski do not
fulfill this exact requirement and do not provide accurate results in the
range of small densities. The results of DFT-MD simulations, which are
currently considered to be the most reliable, are checked by virial plots, and
shortcomings are found in the region of small densities. By benchmarking with
the virial expansion (8) for $x\to 0$, future PIMC simulations can also be
tested. These ab initio simulations become a computational challenge in the
low-density region, but the virial expansion enables extrapolation into this
range. For the construction of interpolation formulas, see Esser03 , knowledge
of the virial coefficients is an important ingredient.

A particular problem concerns the treatment of the $e-e$ interaction in DFT-MD
simulations. To a certain extent, the $e-e$ interaction is included in the
exchange-correlation energy $E_{xc}$ in order to obtain an effective mean
field for optimal single-particle states. However, with the help of virial
expansion it was shown that $e-e$ collisions are not included in the DFT-MD
calculations. By introducing a correction factor $R_{ee}$, which is obtained
from a model plasma with statically screened Coulomb interaction, the first
virial coefficient $\rho_{1}$ can be reproduced. In this work, the problems
with the second virial coefficient $\rho_{2}(T)$, which is determined by the
treatment of dynamic screening, are emphasized. Since the dynamic screening is
not considered in the DFT-MD simulations, the second virial coefficient
$\rho_{2}^{\rm QLB}$ (11) is also not obtained by the correction factor
$R_{ee}$. A numerical calculation can be considered in the context of PIMC
simulations to determine the second virial coefficient $\rho_{2}(T)$.

A correction factor $R_{ee}^{\rm DFT-MD}(T,n)$ was introduced which takes into
account the specific choice of exchange-correlation energy in the DFT
approach. The high-temperature limit of the virial expansion was given in Eq.
(34), and the full $T$ dependence would be of interest. The DFT-MD simulations
are very accurate at higher densities since $e-e$ correlations are strongly
suppressed due to Pauli blocking in the degenerate range ($\Theta\leq 1$).
With the knowledge of $R_{ee}^{\rm DFT-MD}(T,n)$, results for the conductivity
of hydrogen plasmas can also be obtained from DFT-MD calculations for the non-
degenerate range. In addition to the analytical evaluation of correlation
functions for higher moments of the distribution function in the framework of
gLRT, the numerical calculation in the context of PIMC simulations can be
considered to determine the second virial coefficient $\rho_{2}(T)$.

The approach described here can also be applied to other transport properties
such as thermal conductivity, thermopower, viscosity and diffusion
coefficients French22 . It is also interesting to extend the virial expansion
to substances other than hydrogen, where different ions can be formed and the
electron-ion interaction is no longer purely Coulombic.

Acknowledgments

The author would like to thank M. Bethkenhagen, M. French, R. Redmer, H.
Reinholz and M. Schörner for valuable and fruitful discussions.

## References

  * (1) J. L. Spitzer and R. Härm, Phys. Rev. 89, 977 (1953).
  * (2) R. Kubo, J. Phys. Soc. Japan 12, 570 (1957); Rep. Prog. Phys. 29, 255 (1966).
  * (3) A. L. Fetter and J. D. Walecka, Quantum theory of Many-Particle Systems (McGraw-Hill, 1971).
  * (4) G. Röpke, Nonequilibrium Statistical Physics (Wiley-VCH, Weinheim, 2013).
  * (5) G. Röpke, Phys. Rev. A 38, 3001 (1988).
  * (6) G. Röpke and R. Redmer, Phys. Rev. A 39, 907 (1989).
  * (7) H. Reinholz, G. Röpke, S. Rosmej, and R. Redmer, Phys. Rev. E 91, 043105 (2015).
  * (8) G. Röpke, TMF 194, 90 (2018); Engl. transl.: Theoret. and Math. Phys. 194, 74 (2018).
  * (9) H. Reinholz and G. Röpke, Phys. Rev. E 85, 036401 (2012).
  * (10) M. P. Desjarlais, C. R. Scullard, L. X. Benedict, H. D. Whitley, and R. Redmer, Phys. Rev. E 95, 033203 (2017).
  * (11) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * (12) B. Militzer and D. M. Ceperley, Phys. Rev. E 63, 066404 (2001).
  * (13) M. Bonitz et al., Phys. Plasmas 27, 042710 (2020).
  * (14) M. Böhme, Zh. A. Moldabekov, J. Vorberger, and T. Dornheim, Phys. Rev. Lett. 129, 066402 (2022).
  * (15) R. Redmer, Physics Reports 282, 35 (1997).
  * (16) H. Reinholz, R. Redmer, and D. Tamme, Contrib. Plasma Phys. 29, 395 (1989).
  * (17) V. S. Karakhtanov, Contrib. Plasma Phys. 56, 343 (2016).
  * (18) A. Esser, R. Redmer, and G. Röpke, Contrib. Plasma Phys. 43, 33 (2003).
  * (19) G. Röpke, M. Schörner, R. Redmer, and M. Bethkenhagen, Phys. Rev. E 104, 045204 (2021).
  * (20) P. E. Grabowski et al., High Energy Dens. Phys. 37, 100905 (2020).
  * (21) M. French, G. Röpke, M. Schörner, M. Bethkenhagen, M. P. Desjarlais, and R. Redmer, Phys. Rev. E 105, 065204 (2022).
  * (22) G. Röpke et al., arXiv:2310.17583; to be published in Phys. Rev. E.
  * (23) D. A. Greenwood, Proc. Phys. Soc. London 71, 585 (1958).
  * (24) M. P. Desjarlais, J. D. Kress, and L. A. Collins, Phys. Rev. E 66, 025401(R) (2002).
  * (25) S. Mazevet, M. P. Desjarlais, L. A. Collins, J. D. Kress, and N. H. Magee, Phys. Rev. E 71, 016409 (2005).
  * (26) B. Holst, M. French, and R. Redmer, Phys. Rev. B 83, 235120 (2011).
  * (27) M. French and R. Redmer, Phys. Plasmas 24, 092306 (2017).
  * (28) M. Gajdos, K. Hummer, G. Kresse, J. Furthmüller, and F. Bechstedt, Phys. Rev. B 73, 045112 (2006).
  * (29) N. R. Shaffer and C. E. Starrett, Phys. Rev. E 101, 053204 (2020).
  * (30) G. Röpke, Contrib. Plasma Phys. 2023, e202300002. https://doi.org/10.1002/ctpp.202300002.
  * (31) V. S. Karakhtanov, R. Redmer, H. Reinholz, and G. Röpke, Contrib. Plasma Phys. 53, 639 (2013).
  * (32) H. Reinholz,Ann. Phys. Fr. 30, 1 (2005).
  * (33) H. Reinholz, R. Redmer, G. Röpke, and A. Wierling, Phys. Rev. E 62, 5648 (2000).
  * (34) Y. Lee and R. More, Phys. Fluids 27, 1273 (1983).
  * (35) N. R. Shaffer and C. E. Starrett, Phys. Rev. E 101, 013208 (2020).
  * (36) C. E. Starrett, High Energy Density Phys. 25, 8 (2017).


# Perceptive Locomotion with Controllable Pace
and Natural Gait Transitions Over Uneven Terrains

Daniel Chee Hian Tan1∗, Jenny Zhang1∗, Michael (Meng Yee) Chuah1 and Zhibin
Li2 1Authors are with Institute for Infocomm Research, A*STAR,
Singapore2Author is with the Department of Computer Science, University
College London, United Kingdom∗Equal contribution

###### Abstract

This work developed a learning framework for perceptive legged locomotion that
combines visual feedback, proprioceptive information, and active gait
regulation of foot-ground contacts. The perception requires only one forward-
facing camera to obtain the heightmap, and the active regulation of gait paces
and traveling velocity are realized through our formulation of CPG-based high-
level imitation of foot-ground contacts. Through this framework, an end-user
has the ability to command task-level inputs to control different walking
speeds and gait frequencies according to the traversal of different terrains,
which enables more reliable negotiation with encountered obstacles. The
results demonstrated that the learned perceptive locomotion policy followed
task-level control inputs with intended behaviors, and was robust in presence
of unseen terrains and external force perturbations. A video demonstration can
be found here, and the codebase here.

##  I Introduction

Legged locomotion has achieved commendable performance without visual
perception of the external environment, which can negotiate certain vertical
heights via reactive reflexes [1, 2, 3]. However, perceptive locomotion
enables proactive terrain negotiation, anticipating and avoiding obstacles
[4], augmenting traversability of various terrains. The recent DARPA
subterranean challenge [5] exemplified the need for legged robots capable of
terrain-adaptive locomotion, in order to achieve better safety and efficiency
for mission completion in many real-life applications. Overall, legged robots
are becoming commodities and more advantageous compared to their wheeled and
tracked counterparts.

The required technologies for doing so boils down to robot perception and
control. In well established robotics streams, perception involves engineering
the extraction of visual features that can be modularly connected with
designed controllers, i.e. integrating SLAM and robot state estimation with
the existing blind locomotion modules, such as optimization-based motion
planning [6], multi-skill locomotion, [7], reinforcement learning based
trajectory adaptation [8], and meta-reinforcement learning [9]. Despite this
modularity and interpretability, traditional control architectures without
learning require significant engineering effort and expertise over multiple
disciplines for such systems engineering, resulting in longer development
cycles.

In contrast, learning-based approaches can bridge sensing and closed-loop
control by directly feeding perceptive data into the control policy. The
ability to directly design an end-to-end controller offers better scalability
in increasing the dimensionality for accessing both proprioceptive and
exteroceptive data within the feedback loop – a more integral and compact
scheme to design a perception-action control system with environment-aware
behaviors.

The real-world challenges in DARPA subterranean finals in September 2021 [5]
found that apart from perceiving irregular surfaces, it is equally important
to control different gait frequencies at different target traveling
velocities, in order to increase the success rate of terrain negotiation. This
is specially observed when many robots fell due to keeping an almost constant
trotting frequency over different traveling speed on uneven terrains,
resulting in them toppling over and mission failures.

This work aims to develop computationally efficient and effective feedback
control for terrain-aware motor skills, excluding the long-term motion
planning. Our contributions are: (i) active obstacle negotiation with only one
front camera using depth images; (ii) a learning framework with decoupled
modules for perception and motor control policy; (iii) controllable gait
frequency and velocity as high-level inputs for operators to regulate robot
gaits in accordance with different terrains.

The proposed scheme offers cost-effective solutions for perceptive locomotion,
and the straightforwardness of each sub-module makes them easily
implementable, because the perception requires only one off-the-shelf depth
camera and can be re-trained separately to embed new targeted terrain
features, and the computation of a forward neural network (NN) is very low
practically. The required PD control of robot joints is well understood and
straightforward to implement.

##  II Related Work

###  II-A Perceptive Locomotion

To augment traversability of various unstructured terrains, perceptive
locomotion offers advantages of proactive terrain negotiation without
requiring physical contact or trial and error as in blind walking [4]. For
visual perception, camera feeds are commonly used as input to a neural network
to predict safe footholds [10, 11], and LiDAR sensors can be used with
proprioceptive information to produce environment-aware control [12, 13].
Hybrid methods of RL and trajectory optimization have been used to produce
adaptive locomotion from estimated terrain heightmaps [14]. Through a training
curriculum, it was shown that motor policies can be learned using sparse
visual observations to achieve terrain-aware locomotion [15].

Here, we are motivated to explore the direct use of off-the-shelf onboard
cameras (e.g., RealSense D405) without needing complex SLAM information as
part of the feedback. We have investigated effective learning-based control
framework to use only the front-facing heightmap to encode terrain information
and achieve active traversal of various terrains at different gait paces.

###  II-B Motion Imitation for Desirable Behaviors

DeepMimic [3] and the related works [1, 16] used motion-capture to generate
full and highly-detailed joint trajectories, body position and velocity, body
orientation and angular velocity, among others. We refer to such detailed use
of reference as full configuration-level of motion imitation. The policy is
trained with an off-the-shelf RL algorithm, with a reward component to
incentivize imitation of the reference trajectories [3].

Although motion imitation incorporates prior knowledge about preferred skills,
allowing policy-learning to succeed with much simpler reward functions.
However, this too has its own limitation of requiring a suitable dataset of
reference motions, which is nontrivial to obtain. Furthermore, the same
inductive bias may in turn reduce the flexibility for the policy to learn
suitable terrain-adaptive locomotion skills. For example, DeepMimic
demonstrated limited ability to navigate ground terrain obstacles as compared
to recent methods using vanilla RL [2, 4, 14].

###  II-C Controllable Gaits using Central Pattern Generators

Central pattern generators (CPGs) are used here to achieve high-level motion
imitation. CPGs have been studied as a simple and effective biomimetic
mechanism for locomotion [17]. Past use of CPGs for locomotion largely focus
on the direct generation of a fixed motion pattern in the configuration space
[18, 19, 20, 21]. The resulted motion has limited flexibility to adapt to
diverse terrains.

Improving on this, we propose to use CPGs to generate foot-ground references
to guide a learning policy, such that it is incentivized to imitate the
reference gait and also has more flexibility to optimize around the foot-
ground references to traverse terrain obstacles more successfully. It is clear
that our CPG-guided framework is inherently different from the traditional way
of combining CPG joint angle trajectories with neural network residual joint
offsets. Our proposed method has the advantage that the learning has more
flexibility to learn richer terrain-adaptive motions in the joint space.

Recently, [22] demonstrated the advantages and feasbility of CPGs in
generating reference foot position trajectories, which were subsequently used
as a motion imitation primitive in combination with a learning-based
locomotion framework. However, the analysis was limited to a small number of
discrete, fixed gaits. We extend this further to demonstrate feasibility
across a continuous range of gait parameters. Furthermore, we reduce the bias
towards a particular fixed locomotion behavior by replacing the reference foot
position trajectories with reference foot contact states. Our new proposal
allows the policy more leeway to explore and learn terrain-adaptive
locomotion.

##  III Methods

###  III-A Proposed Architecture of Perceptive Locomotion

Figure 1: Architecture of learning perceptive locomotion by encoding both
visual and proprioceptive feedback, with CPG-based gait regulation. Additional
task-control inputs and CPG-based foot-ground contact regulation allow to
adjust trotting velocity and gait frequencies on-the-fly. Terrain perception
is represented by a latent vector from a heightmap encoder pretrained via the
scanning of various terrain surfaces.

The control architecture for our learning-based control policy is presented in
Fig. 1. The architecture is identical during training and testing. At every
time step $t$, the policy receives a state observation
$s_{t}=(p_{t},\hat{h_{t}},\hat{f_{t}},\hat{v_{t}})$ containing proprioceptive
state vector $p_{t}$, heightmap embedding $\hat{h_{t}}$, reference foot
contacts $\hat{f_{t}}$, and reference velocity $\hat{v_{t}}$. The
proprioceptive state $p_{t}\in\mathbb{R}^{54}$ includes the agent’s linear
velocity $(v^{t}_{x},v^{t}_{y},v^{t}_{z})$, and a 3-step history of body
orientation $(\theta^{t-2:t}_{roll},\theta^{t-2:t}_{pitch})$, body angular
velocity $(\omega^{t-2:t}_{roll},\omega^{t-2:t}_{pitch},\omega^{t}_{yaw})$,
and joint angles $(q^{t-2:t})$.

The heightmap embedding $\hat{h_{t}}\in\mathbb{R}^{32}$ consists of the
encoded output of the distance measurements obtained from the $12\times
16=192$ rays used in the heightmap setup. The heightmap measurement provides
distance to physical obstacle, and is clipped between 0.1 m and 8.0 m to
replicate the specifications of commercially available RGB-D cameras. More
details on the heightmap setup and encoder can be found in Section III-C.

The reference foot contacts $\hat{f_{t}}\in\mathbb{R}^{20}$ includes
$(\hat{f_{t}},\hat{f_{t+1}},\hat{f_{t+2}},\hat{f_{t+10}},\hat{f_{t+50}})$.
More details on the reference foot contacts and central pattern generator can
be found in Section III-D. From our experiments, allowing the policy to
observe future reference states enables better tracking of the reference foot
contact trajectory. In this paper, the indices are chosen to achieve a balance
between short-term and long-term future observations.

The target velocity vector $\hat{v_{t}}\in\mathbb{R}^{2}$ specifies the speed
at which the agent should move along the $x$ and $y$ axes. The policy $\pi$
outputs action $a_{t}$, which are the desired joint positions. These are fed
into joint PD control to produce joint torques for simulating dynamics in the
physics engine. The corresponding state measurements are then fed back into
the policy to close the feedback loop.

###  III-B Learning Locomotion via Reinforcement Learning

We formulate perceptive locomotion as an RL problem. Given states $s_{t}$ and
Markovian dynamics $p(s_{t+1},r_{t}|s_{t},a_{t})$ where $a_{t}$ are actions
and $r_{t}$ is the reward, the objective is to learn a control policy
$\pi(a_{t}|s_{t})$ that enables an agent to maximize its expected return for a
locomotion task.

At each time step $t$, the agent observes a set of states $s_{t}$ from the
environment and samples a set of actions $a_{t}\sim\pi(a|s_{t})$ from its
policy $\pi$. The agent then applies $a_{t}$, which results in a new state
$s_{t+1}$ and a scalar reward $r_{t}=r(s_{t},a_{t},s_{t+1})$. Such a repeated
process generates a trajectory
$\sigma=\\{(s_{k},a_{k},r_{k})\\}_{t\in\mathbb{N}}$. The objective is to
maximize the expected return $J(\pi)$:

$J(\pi)=\mathbb{E}_{\sigma\sim
p(\sigma|\pi)}\bigg{[}\sum_{t=0}^{N-1}\gamma^{t}r_{t}\bigg{]}$ (1)

where $N$ denotes the duration of the episode, and $\gamma\in[0,1]$ a time-
discount factor. Given that $r_{t}$ is bounded, the sum on the R.H.S converges
to a finite limit. $p(\sigma|\pi)$ denotes the likelihood of a trajectory
$\sigma$ under policy $\pi$:

$p(\sigma|\pi)=p(s_{0})\prod_{t=0}^{N-1}p(s_{t+1}|s_{t},a_{t})\pi(a_{t}|s_{t})$
(2)

with $p(s_{0})$ being the initial state distribution, and
$p(s_{t+1}|s_{t},a_{t})$ representing the dynamics of the system, which
determines the effects of the agent’s actions.

####  III-B1 Reward Design for Perceptive Locomotion

The reward function encourages the robot to move at the given reference
velocity $\hat{v_{t}}\in\mathbb{R}^{2}$ and follow sequences of reference foot
contacts, and penalizes it for unstable and inefficient motions. At every time
step $t$, we denote the agent’s linear velocity as $\boldsymbol{v}_{t}$, its
four foot contacts as vector $\boldsymbol{f}_{t}$, orientation as
$\boldsymbol{\theta}_{t}$, joint torques as $\boldsymbol{\tau}_{t}$, joint
velocities as $\boldsymbol{\dot{q}}_{t}$.

The reward term to encourage following the reference velocity $\hat{v_{t}}$ is
given by:

$r^{t}_{v}=\begin{cases}e^{-c_{\alpha}(d_{tgt})^{2}}&\text{if $d_{tgt}<0$}\\\
e^{1-\frac{1}{1-c_{\beta}(d_{tgt})^{2}}}&\text{if $d_{tgt}\geq 0$ and
$d_{tgt}<\frac{1}{c_{\beta}}$}\\\ 0&\text{otherwise}\end{cases}$ (3)

where
$d_{tgt}=\frac{\hat{\boldsymbol{v}_{t}}\cdot\boldsymbol{v}_{t}}{||\hat{\boldsymbol{v}_{t}}||^{2}}$
is the normalized distance from base and target velocity vectors, and
$c_{\alpha}$ and $c_{\beta}$ are pre-determined reward parameters as
$c_{\alpha}=16.0$ and $c_{\beta}=25.0$ for training the policy.

A Gaussian function is used when the agent’s velocity is less than the target
velocity in the target velocity’s direction. This is to help the agent to
learn with less restricted penalty that allows it to slow down and engage with
the terrain when necessary. On the other hand, a bump function is used when
the agent’s velocity is more than the desired one in the targeted direction.
This is to enforce the policy to follow the given reference velocity.

The reward term to penalize instability is given by:
$r^{t}_{s}=-||\begin{bmatrix}0,0,0\\\ 0,0,0\\\ 1,1,0\\\
\end{bmatrix}\cdot\boldsymbol{R}_{t}||$ , where $\boldsymbol{R}_{t}$ is the
rotation matrix derived from the quaternion of the agent’s orientation
$\boldsymbol{\theta}_{t}$. The reward term to encourage following the
reference foot contacts from control input is given by:
$r^{t}_{f}=\frac{1}{4}(\hat{\boldsymbol{f}_{t,0}}\cdot\boldsymbol{f}_{t})$.
The reward term to penalize energy usage is given by:
$r^{t}_{e}=-|\boldsymbol{\tau}_{t}\cdot\boldsymbol{\dot{q}}_{t}|$. For
training the policy, the total reward at any time $t$ is defined as the sum of
the weighted reward terms:
$r_{t}=w_{v}*r^{t}_{v}+w_{f}*r^{t}_{f}+w_{s}*r^{t}_{s}+w_{e}*r^{t}_{e}$, where
$w_{v}=1.0$, $w_{f}=2.0$, $w_{s}=1.5$, $w_{e}=0.0001$.

####  III-B2 Training Curriculum

This work focuses on quadruped trotting as it is the most common and useful
gait in real-world applications. The gait frequency is sampled uniformly
within $[1.5,2.5]$ Hz, and the reference velocity is sampled uniformly within
$[0.4,0.8]$ m/s. Gait frequency and target velocity are sampled at the start
of each episode and are held constant throughout the episode.

Three platforms and stairs of different heights are used in training - 2 cm, 5
cm, and 7 cm. We use a fixed-order curriculum for the terrains during
training: 2 cm $\rightarrow{}$level 1; 5cm $\rightarrow{}$level 2; 7 cm
$\rightarrow{}$level 3. At the start of the training, the agent begins at
level 1. When the agent has successfully traversed up and down the terrain at
the current level, then its level is increased by one. Then at the next reset,
the agent will start on a randomly selected terrain at its new level or lower.
The next selected terrain is proportionally sampled according to its
difficulty.

####  III-B3 Terrain Descriptions

Other than (i) platforms, we also tested the trained neural network controller
on other terrains: (ii) Hurdles: narrow obstacles designed to evade perception
and trip the robot. They have a step run of 0.06 m. In Figures 9 and 10, the
step height is specified between 0.03 and 0.06 m. (iii) Heightfield: a rugged
terrain with obstacles protruding from the ground. The specified heightfield
terrain area is split into 0.1 m $\times$ 0.1 m grids, where each grid unit is
set to a random height uniformly selected between 0 and the specified
perturbation height. In Figures 9 and 10, the perturbation height is set to
0.05 m. (iv) Stairs: staircases with 5 steps. Each step has a run of 0.3 m. In
Figures 9 and 10, the range of the height of steps is specified between 0.03
and 0.07 m.

###  III-C Exteroceptive Perception of Terrains

The robot updates the heightmap data $h_{t}$ in front of it at each time step
$t$. As seen in Fig. 2, the heightmap is designed with a dimension of 12
$\times$ 16 rays, each grid in the map is spaced horizontally and vertically
0.03 m apart and contains the height coordinate $z$ of the scanned area. With
a robot base height of 0.35 m, the field of view of is about 70°$\times$ 55°,
which can be easily covered by commercial-off-the-shelf RGB-D cameras.

Figure 2: Configuration of visual perception and its feature encoding: (a)
Setup of heightmap; (b) Heightmap autoencoder, and the green shaded area is
the encoder used in policy training.

To achieve such dimension reduction, we first train an autoencoder to compress
the heightmap data, see Fig. 2 for the overall architecture. Mean square error
is used as the loss function. Heightmap data is collected from simulation with
a camera fixed on a “hovering” robot across different terrain types. We used
only the encoder part of the trained heightmap autoencoder, and compressed the
exteroceptive heightmap observation $h_{t}\in\mathbb{R}^{12\times 16}$ to a
lower dimension state of $\hat{h_{t}}\in\mathbb{R}^{32}$. During the training
of motor policy, the weights of all layers in the encoder are fixed to produce
heightmap embedding $\hat{h_{t}}$, which is fed into the policy.

###  III-D Central Pattern Generators and Gait Representation

Our motivation of using CPGs is to have active control of the gait pattern in
the form of foot-ground contact sequences, which gives the RL agent more
flexibility to explore joint-space motions and diverse terrain-adaptive
skills. Our specifically formulated CPG approach can control both the contact
frequency and duty cycle of the trotting gait.

Traditional RL approach normally specifies target velocity instead of contact
patterns. Hence, for different velocities, the robot changes the step length
only. In contrast, our method adapts from [22] and formulates CPGs to generate
a total of 4 abstract phases (one per leg). The abstract phase vector is
decoded to a reference foot contact state, which is used as state input for
the neural network. As a result, the CPG does not impose joint angles, but
guides the policy only in terms of ground contacts. With the same foot-contact
references, the foot clearance and stride length are free to explore, enabling
a considerable range of terrain-adaptive behaviors.

We found that without the CPG input, the quality of learned locomotion
degrades significantly. For conciseness, we provide a brief recapitulation
below. A CPG is a dynamical system that produces rhythmic patterns with a
stable limit cycle. Fig. 3 shows the Hopf oscillator as an example of a CPG,
with converging limit cycles. A Hopf oscillator has parameters
$(\alpha,\beta,b,T,\mu,\gamma)$ and state $\rho=(x,y)\in\mathbb{R}^{2}$,
radius $r=\sqrt{x^{2}+y^{2}}$ with intrinsic update equations:

$\dot{x}=\alpha(\mu^{2}-r^{2})x+\gamma
y,\quad\dot{y}=\alpha(\mu^{2}-r^{2})y-\gamma x$ (4)

In our work, parameters are $\alpha=10.0,b=50.0,\mu=1.0,\gamma=50.0$. $T,b$
are the parameters to modulate the desired gait frequency and duty factor
respectively.

Figure 3: Advantages of CPG: (a) Convergence to a stable limit cycle (Equation
4); (b) Four synchronized coupled Hopf oscillators for trotting (Equation 5).

We initialize four Hopf oscillators with identical parameters, one for each
leg. It is desired that the four oscillators are synchronized, maintaining a
fixed phase offset $\phi\in\mathbb{R}^{4}$, where $\phi_{i}$ is the desired
phase of CPG $i$. To enable synchronization of the four oscillators, we modify
the intrinsic update equations in (4) with a coupling term $R_{ij}$:

$\displaystyle\dot{\rho}_{i}=f(\rho_{i})+\sum_{ij}R_{ij}\rho_{j},\quad
R_{ij}=\begin{bmatrix}\cos\theta_{ij}&-\sin\theta_{ij}\\\
\sin\theta_{ij}&\cos\theta_{ij}\\\ \end{bmatrix}$ (5)
$\displaystyle\theta_{ij}=\phi_{j}-\phi_{i},\quad i,j\in\\{1,2,3,4\\}$

Fig. 3 shows an example of how four oscillators beginning from a random
initial state can achieve convergence when they are updated according to
Equation 5.

Gaits are parameterized as binary reference trajectories
$f\in\\{-1,1\\}^{N\times 4}$. To generate $f_{ti}$, where $t$ is time and $i$
is the leg index, for each oscillator state $(x_{ti},y_{ti})$, the
corresponding leg phase is $\hat{\theta}_{ti}=\arctan(\frac{y_{ti}}{x_{ti}})$,
which is decoded to a binary foot contact state as
$f_{ti}=2\cdot\mathrm{I}[\hat{\theta}_{i}>0]-1$.

As a result, we obtain $f$ where $f_{ti}=1$ if and only if foot $i$ should be
on the ground at time $t$, and $-1$ otherwise. This method of representing
gaits as binary trajectories enables concise encapsulation of key properties:
(i) gait frequency which controls how rapidly a rhythmic gait is generated,
(ii) gait duty factor which represents the stance-to-swing duration ratio, and
(iii) phase offsets which determines the class of gait being executed (e.g.
walk, trot, canter, pace). In this work, we focus on the trotting gait which
is the mostly used in real applications, and therefore only regulate the
frequency part and keep the latter two constant. Fig. 4 shows the trajectories
of binary foot-ground contract for each gait frequency.

Figure 4: Foot-contact references with 0.5 duty factor and varying trotting
frequencies.

Our CPG-guided framework provides a balance between vanilla RL and direct
motion imitation. By incorporating an imitation term of contact reference, we
create an inductive bias towards known stable gait patterns, while keeping the
simplicity of reward formulation. Because the desired joint motions are not
hard-coded in the contact references, there is considerable flexibility for
the RL policy to learn terrain-adaptive skills.

##  IV Simulation Results

###  IV-A Setup of Simulation and Training

The NN policy is an actor-critic architecture with parameter sharing. The
shared backbone is a 2-layer network, each having 256 nodes and ReLU
activation. The linear prediction heads in the last layer are used to obtain
the value estimate and action distribution.

The NN policy is queried at a control frequency of 40 Hz, which predicts
desired local joint angles $a=\hat{\boldsymbol{q}}\in\mathbb{R}^{12}$. The
target joint positions are post-processed as: (i) clipped to a maximum of 0.5
radians from the default resting pose; (ii) velocity is limited by a maximum
change of 0.15 radians from the current position; and (ii) filtered using a
first-order Butterworth low-pass filter in case of outliers or spikes. Using
PD control, such updated actions $\hat{a}$ generate joint torques $\tau$ as:
$\tau=K_{p}(\hat{\boldsymbol{q}}-\boldsymbol{q})+K_{d}(\hat{\dot{\boldsymbol{q}}}-\dot{\boldsymbol{q}})$.
$K_{p}$ and $K_{d}$ are the proportional and derivative gains, and the target
joint velocities $\hat{\dot{\boldsymbol{q}}}$ are 0 here. For each leg,
$K_{p}$ are $(100,100,100)$ and $K_{d}$ are $(1,2,2)$ for the hip, thigh, and
knee joints respectively.

Parameter |  Value
---|---
Discount factor $\gamma$ |  0.95
Learning rate |  $3\mathrm{e}-5$
PPO clip threshold |  0.4
PPO batch size |  1024
PPO epochs |  20
TABLE I: PPO hyperparameters for training.

The learning task is implemented as an OpenAI Gym [23] environment and uses
PyBullet [24] for the physics, rigid-body and contact dynamics simulated at
1000 Hz. At each simulation time step, the latest target joint position is
used in PD control to generate motor torques. This work used the simulated
Unitree A1 robot [25] and its associated kinematic and dynamic parameters.

The policy was trained using Proximal Policy Optimization (PPO) [26]
implemented in Stable Baselines3 [27]. See Table I for the hyperparameters
used. Each training episode lasts for a maximum of 1000 discrete time steps.
To reduce unnecessary computation, we employ early termination of episodes
when the robot loses balance.

The learned policy was evaluated on the following terrains: (i) Platforms,
raised platforms of progressively increasing height; (ii) Hurdle, narrow
obstacles designed to evade perception and trip the robot; (iii) Heightfield,
rugged terrain with obstacles protruding from the ground; and (iv) Stairs,
staircases of progressively increasing height.

###  IV-B Controllable Trotting Gait Frequencies and Speed

The trained policy is designed to be able to regulate contact phases, and thus
the frequency of the trotting pace. Fig. 5 shows the executed foot contacts
closely followed different reference gait frequencies of 1.5, 2.0, and 2.5 Hz,
at a fixed target velocity of 0.6 m/s. A Hamming similarity of approximately
0.8 is achieved across all feet contacts.

At different gait frequencies, the trained policy can simultaneously follow a
given reference velocity. Fig. 6 shows the body velocities measured during 4
sets of simulation tests with gait frequencies of 1.5, 2.0, 2.5 and 3.0 Hz.
Each set of tests used reference forward velocities of 0.4, 0.6, 0.8, and 1.0
m/s. A total of 20 trials with different random generator seeds were run per
reference gait frequency. The systematic testing shows that a too low-
frequency gait of 1.5 Hz is unable to catch up with a comparatively high
velocities of 0.8 to 1.0 m/s, given the leg length of the robot. Needless to
say, higher-frequency gaits of 2.0, 2.5 and 3.0 Hz are able to do so. Results
in Fig. 6 show that with higher gait frequency, meaning a faster pace, the
robot follows 1.0 m/s reference better. This implies a natural correlation
between gait frequency and achievable traveling speed, given the size of a
physical robot, i.e. the higher the gait frequencies, the faster the robot can
trot.

Figure 5: Contact phases of the front right foot across different reference
gait frequencies at reference velocity 0.6 m/s, and Hamming similarity across
all feet.  Figure 6: Traveling velocities measured across tests at different
reference gait frequencies.

###  IV-C Terrain-Adaptive Gait with Natural Transitions

Further validations were done to demonstrate that the learned policy can
achieve seamless transitions between gaits. This allows to select gait
parameters and enforce how fast the robot navigates over different obstacles,
e.g., cautious slow pace on difficult surfaces and faster gait on flat ground
for maximizing task efficiency.

We evaluated such capability of terrain-adaptive gait transitions by using
scheduled control inputs (see Fig. 1), which actively modulate gait
transitions as the robot is crossing terrain obstacles. On flat ground, a gait
of 3.0 Hz and 0.8 m/s is used. When encountering terrain obstacles, this is
linearly reduced to 0.4 m/s. A total of 20 trials were tested per terrain.
Fig. 7 shows that the robot can smoothly follow the varying control inputs
while crossing various terrain.

Figure 7: Measured body velocities and foot contacts which were tested across
different terrains with controlled gait transitions. The reference gait
frequency and velocity start to decrease while encountering an area with rough
terrains or obstacles. Figure 8: Comparison of locomotion over a 9cm platform
with and without visual perception.

Figure 9 shows that the robot can smoothly follow the varying control inputs
while crossing other terrains than just platforms. On flat ground, a gait of
3.0 Hz and 0.8 m/s is used. When encountering terrain obstacles, this is
linearly reduced to 2.0 Hz and 0.4 m/s. A total of 20 trials with different
random generator seeds were tested per terrain.

Figure 9: Measured body velocities and foot contacts which were tested across
different terrains with controlled gait transitions. The reference gait
frequency and velocity start to decrease while encountering an area with rough
terrains and/or obstacles.

###  IV-D Robustness and Viability Tests

Due to unavailability of hardware, we tested the robustness in simulation to
see viability to implement on real robots. All our tests have validated that
the physical constraints, i.e., joint positions, velocities, and motor torques
are always within the hardware limits.

Fig. 8 demonstrates the effectiveness of our designed perception in the
traversal over uneven terrains. By perceiving the heightmap, the robot is able
to proactively step over a 9 cm platform. However, without perception, the 9
cm-step is comparatively high (body height is 25 cm during trotting); relying
on the proprioceptive data only, the robot failed to achieve sufficient foot-
ground clearance blindly and is hindered from moving forward.

Figure 10 shows that the robot can robustly recover from unseen collisions and
return to the reference path while crossing other terrains than just
platforms. On flat ground, a gait of 3.0 Hz and 0.8 m/s is used. When
encountering terrain obstacles, this is linearly reduced to 2.0 Hz and 0.4
m/s, similar to Figure 9. A total of 20 trials with different random generator
seeds were tested per terrain.

Figure 10: Robust perceptual locomotion over various terrains by rejecting
impact disturbances. Black arrows indicate the timing and direction of the
flying objects. In all test cases, the learned policy robustly recovers back
to the reference path.

Moreover, we demonstrate robustness of our learned policy against unexpected
external collisions. We initialized arbitrary objects of 3.0 kg moving at 10.0
m/s to impact the robot laterally. The impacts occur at the moment before the
robot encounters terrain obstacles, increasing the difficulty of terrain
negotiation. Despite these perturbations, our policy robustly recovers from
the impact and returns to the reference path after being pushed away. Fig. 11
shows the statistical results achieved across 20 unique seeds on platforms.
Without perturbations, the learned policy can achieve $100\%$ success rate
across all terrains. With perturbation, it achieved $85\%$, $100\%$, $90\%$
completion rate on Platforms, Hurdles, Heightfield respectively.

Figure 11: Robust perceptive locomotion over various terrains by withstanding
external force impacts. Black arrows indicate time and direction of the flying
objects. In all test cases, the learned policy robustly recovers back to the
reference path.

Though our work is simulation-only, nonetheless, we argue that sim-to-real
deployment is theoretically feasible. Sim-to-real methods for legged
locomotion have been extensively reported in the literature, such as domain
randomization [1], rapid motor adaptation [2], and fine-tuning in the real
world [16]. All of these are compatible with the CPG-guided learning framework
that we present in this paper.

To investigate whether our policy transfers well to realistic scenarios, we
tested the policy with additional levels of sensor noise over flat ground.
Results are presented in Table II. We note that velocity error and foot
contact similarity does not degrade significantly even as noise level
increases up to 10x the nominal noise level, demonstrating robustness to
random sensor noise. The nominal noise level for the heightmap sensor is
determined to be 0.005 m. This is obtained by multiplying 2% nominal depth
accuracy with the estimated 0.25 m average robot base height.

Noise | $v_{err}$ | $s$
---|---|---
1x | 0.10 | 0.88
5x | 0.16 | 0.86
10x | 0.12 | 0.83
25x | 0.37 | 0.69

Parameter | 1x Noise std. dev
---|---
Motor Position | 0.01 rad
Base RPY | 0.01 rad
Base ang. vel. | 0.1 rad/s
Heightmap rays | 0.005 m

TABLE II: Evaluation noise levels: (a) Gait metrics at varying levels of
noise. $v_{err}$ is RMSE in base velocity, and $s$ denotes Hamming similarity;
(b) Noise distributions at the 1x-noise level. Each is a zero-mean Gaussian
with the given standard deviation.

##  V Conclusion and Future Work

This work presents a learning framework to train perceptive locomotion which
combines the proprioceptive and exteroceptive data, and task-level enforcement
of gait frequency and velocity with natural transitions. This offers high-
level inputs to modulate gait paces for terrain-adaptive locomotion.
Specifically, the perceptive module uses the latent representation of
exteroceptive data effectively, and demonstrates robust tracking of reference
gait frequencies and velocities across different terrains. The results show
that during the “entry phase” of active negotiation with uneven terrains, the
front-facing visual perception is more critical, as rear legs can handle steps
blindly reasonably well using proprioceptive data. Moreover, the learned
policy is robust to unseen external perturbations which were not trained
before.

However, several limitations remain. Firstly, terrain negotiation can be
further improved. In the benchmarking test ($2.0$ Hz gait frequency, $0.6$ m/s
target velocity), given the same training parameters, our current policy has
different negotiation ability across different terrains: traversing steps and
hurdles up to 11 cm, stairs terrain up to 9 cm, and heightfield terrains with
perturbation up to 8 cm. In future work, additional training time and an
improved training curriculum may enable even better terrain negotiation
ability.

Secondly, the field of view from visual perception is only for the front legs
and not the rear ones. In the current setup, the policy has no information
about the terrain around the rear legs. The reason for such a camera
configuration is that nature shows that animals are able to negotiate terrains
with only front-seeing eyes, so we aim to keep the system with simplicity.
Within this work, we kept the exteroceptive heightmap at a straightforward
setting with low dimension, such that the rear legs rely on proprioceptive
information to blindly step over hurdles. Clearly, this simplicity trades off
the completeness of information and can be improved by adding a back camera.
In future work, including an additional rear heightmap may improve the
policy’s robustness and performance.

## References

  * [1] X. B. Peng _et al._ , “Learning Agile Robotic Locomotion Skills by Imitating Animals,” _arXiv:2004.00784 [cs]_ , Jul. 2020, arXiv: 2004.00784. [Online]. Available: http://arxiv.org/abs/2004.00784
  * [2] A. Kumar _et al._ , “RMA: Rapid Motor Adaptation for Legged Robots,” 2021\. [Online]. Available: https://arxiv.org/abs/2107.04034
  * [3] X. B. Peng _et al._ , “AMP: Adversarial Motion Priors for Stylized Physics-Based Character Control,” _ACM Transactions on Graphics_ , vol. 40, no. 4, pp. 1–20, Aug. 2021, arXiv: 2104.02180. [Online]. Available: http://arxiv.org/abs/2104.02180
  * [4] T. Miki _et al._ , “Learning robust perceptive locomotion for quadrupedal robots in the wild,” _Science Robotics_ , vol. 7, no. 62, p. eabk2822, Jan. 2022, arXiv:2201.08117 [cs]. [Online]. Available: http://arxiv.org/abs/2201.08117
  * [5] T. Rouček _et al._ , _DARPA Subterranean Challenge: Multi-robotic Exploration of Underground Environments_ , 03 2020, pp. 274–290.
  * [6] P. Fankhauser _et al._ , “Robust rough-terrain locomotion with a quadrupedal robot,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , 2018, pp. 5761–5768.
  * [7] C. Yang _et al._ , “Multi-expert learning of adaptive legged locomotion,” _Science Robotics_ , vol. 5, no. 49, p. eabb2174, 2020.
  * [8] S. Gangapurwala _et al._ , “Real-time trajectory adaptation for quadrupedal locomotion using deep reinforcement learning,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , 2021, pp. 5973–5979.
  * [9] T. Anne _et al._ , “Meta-learning for fast adaptive locomotion with uncertainties in environments and robot dynamics,” in _2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2021, pp. 4568–4575.
  * [10] S. Bazeille _et al._ , “Vision enhanced reactive locomotion control for trotting on rough terrain,” in _2013 IEEE Conference on Technologies for Practical Robot Applications (TePRA)_ , 2013, pp. 1–6.
  * [11] A. Siravuru _et al._ , “Deep visual perception for dynamic walking on discrete terrain,” _CoRR_ , vol. abs/1712.00916, 2017. [Online]. Available: http://arxiv.org/abs/1712.00916
  * [12] M. Fallon, “Accurate and robust localization for walking robots fusing kinematics, inertial, vision and lidar,” _Interface Focus_ , vol. 8, p. 20180015, 08 2018.
  * [13] Z. Liu _et al._ , “Learning vision-guided dynamic locomotion over challenging terrains,” _CoRR_ , vol. abs/2109.04322, 2021. [Online]. Available: https://arxiv.org/abs/2109.04322
  * [14] S. Gangapurwala _et al._ , “RLOC: Terrain-aware legged locomotion using reinforcement learning and optimal control,” _IEEE Transactions on Robotics_ , pp. 1–20, 2022. [Online]. Available: https://doi.org/10.1109%2Ftro.2022.3172469
  * [15] F. Acero _et al._ , _IEEE Robotics and Automation Letters_ , vol. 7, no. 4, pp. 8611–8618, 2022.
  * [16] L. Smith _et al._ , “Legged robots that keep on learning: Fine-tuning locomotion policies in the real world,” 2021. [Online]. Available: https://arxiv.org/abs/2110.05457
  * [17] A. J. Ijspeert, “Central pattern generators for locomotion control in animals and robots: A review,” _Neural Networks_ , vol. 21, no. 4, pp. 642–653, 2008, robotics and Neuroscience. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0893608008000804
  * [18] C. Liu _et al._ , “Central pattern generator inspired control for adaptive walking of biped robots,” _IEEE Transactions on Systems, Man, and Cybernetics: Systems_ , vol. 43, no. 5, pp. 1206–1215, 2013.
  * [19] S. Rutishauser _et al._ , “Passive compliant quadruped robot using central pattern generators for locomotion control,” in _2008 2nd IEEE RAS & EMBS International Conference on Biomedical Robotics and Biomechatronics_, 2008, pp. 710–715.
  * [20] B. R. Tietz _et al._ , “Tetraspine: Robust terrain handling on a tensegrity robot using central pattern generators,” in _2013 IEEE/ASME International Conference on Advanced Intelligent Mechatronics_ , 2013, pp. 261–267.
  * [21] D. Drotman _et al._ , “Electronics-free pneumatic circuits for controlling soft-legged robots,” _Science Robotics_ , vol. 6, no. 51, p. eaay2627, 2021\. [Online]. Available: https://www.science.org/doi/abs/10.1126/scirobotics.aay2627
  * [22] Y. Shao _et al._ , “Learning Free Gait Transition for Quadruped Robots via Phase-Guided Controller,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 2, pp. 1230–1237, Apr. 2022, arXiv:2201.00206 [cs]. [Online]. Available: http://arxiv.org/abs/2201.00206
  * [23] G. Brockman _et al._ , “Openai gym,” 2016. [Online]. Available: https://arxiv.org/abs/1606.01540
  * [24] E. Coumans and Y. Bai, “Pybullet, a python module for physics simulation for games, robotics and machine learning,” http://pybullet.org, 2016–2021.
  * [25] “Unitree A1 Robot,” https://www.unitree.com/products/a1/, 2021, [Online; accessed 9-June-2022].
  * [26] J. Schulman _et al._ , “Proximal policy optimization algorithms,” 2017. [Online]. Available: https://arxiv.org/abs/1707.06347
  * [27] A. Raffin _et al._ , “Stable-baselines3: Reliable reinforcement learning implementations,” _Journal of Machine Learning Research_ , vol. 22, no. 268, pp. 1–8, 2021. [Online]. Available: http://jmlr.org/papers/v22/20-1364.html