# Structure-preserving particle methods for the Landau collision operator
using the metriplectic framework

S. Jeyakumar  Mathematical Sciences Institute, Australian National University,
Acton ACT 2601, Australia The University of Western Australia, 35 Stirling
Highway, Crawley WA 6009, Australia  M. Kraus  Max Planck Institute for Plasma
Physics, Boltzmannstraße 2, 85748 Garching, Germany  Technical University of
Munich, Department of Mathematics, Boltzmannstraße 3, 85748 Garching, Germany
M. J. Hole  Mathematical Sciences Institute, Australian National University,
Acton ACT 2601, Australia Australian Nuclear Science and Technology
Organisation, Locked Bag 2001, Kirrawee DC, NSW 2232, Australia  D. Pfefferlé
The University of Western Australia, 35 Stirling Highway, Crawley WA 6009,
Australia

###### Abstract

We present a novel family of particle discretisation methods for the nonlinear
Landau collision operator. We exploit the metriplectic structure underlying
the Vlasov–Maxwell–Landau system in order to obtain disretisation schemes that
automatically preserve mass, momentum, and energy, warrant monotonic
dissipation of entropy, and are thus guaranteed to respect the laws of
thermodynamics. In contrast to recent works that used radial basis functions
and similar methods for regularisation, here we use an auxiliary spline or
finite element representation of the distribution function to this end.
Discrete gradient methods are employed to guarantee the aforementioned
properties in the time discrete domain as well.

###### Contents

  1. 1 Introduction
  2. 2 The Metriplectic formulation
  3. 3 Semi-discretisation
  4. 4 Time discretisation
  5. 5 Conclusion

##  1 Introduction

In this work, we present a structure-preserving particle based approach to
simulating the Landau collision operator, using the metriplectic formulation
of Morrison (1986). The subject of structure-preserving algorithms to simulate
ideal dynamics in plasmas has been well studied (Morrison, 2017). In the
context of ideal kinetics in particular, such as the Vlasov–Maxwell or
Vlasov–Poisson systems, recent works include Chen et al. (2011); Markidis &
Lapenta (2011); Squire et al. (2012); Evstatiev & Shadwick (2013); Qin et al.
(2016); Burby (2017); Kraus et al. (2017); Zhang & Gamba (2017); Campos Pinto
et al. (2022). Complementing the ideal dynamics with dissipative effects, e.g.
those caused by collisions, complicates the construction of structure-
preserving numerical methods for the coupled system and is a more recent topic
of study.

Initial works studying dissipation were primarily based on grid-based
approaches (Yoon & Chang, 2014; Taitano et al., 2015; Hirvijoki & Adams, 2017;
Kraus & Hirvijoki, 2017; Shiroto & Sentoku, 2019), however, there have been
some recent advances in studying particle-based structure-preserving methods
for this application. Hirvijoki et al. (2018) approach this by varying the
weights of the marker particles, instead of their velocities. Carrillo et al.
(2020) and Hirvijoki (2021) use finite-sized marker particles to discretise
the Landau operator, in particular by convolving the delta-function particles
with a finite-sized shape function. Further work by Bailo et al. (2024)
extends such a regularisation approach to the full Vlasov–Maxwell–Landau
system.

Alternative approaches include Monte-Carlo methods, which treat collisions as
a stochastic process and effectively model their underlying microscopic
behaviour (see for e.g. Boozer & Kuo-Petravic (1981)). More recent versions of
such an approach are structure-preserving stochastic methods for collision
operators, including the Landau operator (Tyranowski, 2021a, b). Schnake et
al. (2024) considers a sparse-grid discontinuous Galerkin approach to the
Vlasov–Poisson–Lenard–Bernstein equations, which considers both a sparse grid
for the entire phase space, as well as a hybrid grid which only uses a sparse
grid in velocity space. Finally, the work of Laidin & Pareschi (2024)
considers a constrained polynomial based approach to solving collisional
equations while conserving the required moments, such as momentum and energy.

This manuscript is structured as follows. In Section 2, we give a short
overview of the metriplectic framework, the Vlasov–Maxwell–Landau system and
how the Landau collision operator can be described in terms of a metric
bracket. We pay special attention to the conservation properties and show how
systems that fit into the metriplectic framework are automatically guaranteed
to satisfy the laws of thermodynamics. In Section 3, we describe the
discretisation of the Landau operator by exploiting its metriplectic
structure. We show how an auxiliary spline or finite element representation of
the distribution function can be used for regularisation, prove conservation
of mass, momentum, and energy, and show that entropy is dissipated
monotonically, thereby verifying adherence of the semi-discrete scheme to the
laws of thermodynamics. In Section 4, we show how the discrete gradient method
can be used to obtain a discretisation in time that retains all aforementioned
conservation properties of the semi-discrete scheme.

##  2 The Metriplectic formulation

The metriplectic framework (Morrison, 1986; Kaufman & Morrison, 1982; Kaufman,
1984; Morrison, 1984a, b; Grmela, 1984a, b, 1985; Morrison & Updike, 2023)
provides a convenient way to express dynamical systems which consist of two
parts: a conservative, Hamiltonian part and a dissipative part, which are
compatible with the laws of thermodynamics. The ideal dynamics are described
via a Hamiltonian functional and a Poisson bracket, and the dissipative
dynamics through a symmetric metric bracket and a functional representing the
dissipated quantity, usually the entropy.

In this section, we give an overview of the metriplectic framework, briefly
sketch the Vlasov–Maxwell–Landau system, and then explain the Landau collision
operator in more detail. For a slightly more comprehensive exposition of the
metriplectic framework in the context of collision operators see e.g. Kraus &
Hirvijoki (2017, Section 2.1).

###  2.1 General Framework

Let us denote by $u(t,z)=(u^{1},\,u^{2},\,...,\,u^{m})^{T}$ the field
variables, defined over the domain $\Omega$ with coordinates $z$, and let
$\mathcal{A}$ be an arbitrary functional of the field variables $u$. In the
metriplectic framework, the evolution of any functional $\mathcal{A}$ of the
dynamical variables is given by:

$\dot{\mathcal{A}}=\\{\mathcal{A},\mathcal{F}\\}+(\mathcal{A},\mathcal{F}),$
(1)

where $\mathcal{F}=\mathcal{H}-\mathcal{S}$ is the free energy, i.e., the
difference between the Hamiltonian functional $\mathcal{H}$ and the entropy
functional $\mathcal{S}$. The Poisson bracket $\\{\cdot,\cdot\\}$ describes
the conservative part of the dynamics. It is a bilinear, anti-symmetric
operation that satisfies Leibniz’ rule and the Jacobi identity. The metric
bracket $(\cdot,\cdot)$ describes the dissipative part of the dynamics. It is
a bilinear, positive semi-definite, symmetric operation that satisfies
Leibniz’ rule. With the sign convention chosen here, entropy is always
dissipated.

The two brackets have to satisfy certain compatibility conditions which
ultimately guarantee compliance with the laws of thermodynamics. The entropy
functional is usually chosen to be a Casimir invariant of the Poisson bracket,
meaning that it satisfies $\\{\mathcal{S},\mathcal{A}\\}=0$ for all
functionals $\mathcal{A}$. The metric bracket is constructed to ensure that
the Hamiltonian and other physically relevant invariants $\mathcal{I}$ of the
conservative system (Casimir invariants or momentum maps) are Casimir
invariants of the metric bracket, i.e. $(\mathcal{I},\mathcal{A})=0$ for all
functionals $\mathcal{A}$. Thus, the evolution equation of (1) becomes:

$\dot{\mathcal{A}}=\\{\mathcal{A},\mathcal{H}\\}-(\mathcal{A},\mathcal{S}).$
(2)

The compatibility conditions ensure a separation of the conservative and
dissipative dynamics, as the Hamiltonian dynamics do not affect the entropy,
and the dissipative dynamics do not affect the Hamiltonian or any other
physically relevant invariant of the ideal dynamics. It is desirable to also
preserve this property in the course of discretisation.

###  2.2 H-Theorem and Energy-Casimir principle

With the aforementioned requirements that (i) $\mathcal{H}$ is a Casimir of
the metric bracket, (ii) $\mathcal{S}$ is a Casimir of the Poisson bracket,
and that (iii) the metric bracket is positive semi-definite, the metriplectic
framework automatically reproduces the First and Second Law of Thermodynamics:

$\displaystyle\dfrac{d\mathcal{H}}{dt}$
$\displaystyle=\\{\mathcal{H},\mathcal{H}\\}-(\mathcal{H},\mathcal{S})=0,$
(3a) $\displaystyle\dfrac{d\mathcal{S}}{dt}$
$\displaystyle=\\{\mathcal{S},\mathcal{H}\\}-(\mathcal{S},\mathcal{S})=-(\mathcal{S},\mathcal{S})\leq
0.$ (3b)

When a metriplectic system possesses no Casimir invariants, its equilibrium
state, $u_{eq}$, is determined fully by an energy principle,
$\delta\mathcal{F}[u_{eq}]=0$. Additionally, this state is stable given the
second variation satisfies $\delta^{2}\mathcal{F}[u_{eq}]>0$. When there are
Casimir invariants, however, the equilibrium condition must be modified to
account for the degeneracies in the state of the system (Morrison, 1998). In
this case, the equilibrium state must satisfy an energy-Casimir principle as
follows:

$\bigg{[}\delta\mathcal{F}+\sum_{i}\lambda_{i}\delta\mathcal{C}_{i}\bigg{]}\bigg{|}_{u=u_{eq}}=0,$
(4)

where $\mathcal{C}_{i}$ are the Casimir invariants, and $\lambda_{i}$ are
corresponding Lagrange multipliers which are determined by the initial
conditions of the system.

###  2.3 The Vlasov-Maxwell-Landau system

The Vlasov-Maxwell-Landau system consists of: the Vlasov equation, which
describes the evolution of the phase-space probability distribution, Maxwell’s
equations, which describe the evolution of the electromagnetic fields, and the
Landau operator, which describes collisional interactions between charged
particles. This system is described by the following set of equations:

$\displaystyle\partial_{t}f+v\cdot\nabla_{x}f+\frac{q}{m}(E+v\times
B)\cdot\nabla_{v}f=C[f],$ (5) $\displaystyle C[f]=\frac{\partial}{\partial
v}\cdot\int U(v-v^{\prime})\left(f(v^{\prime})\frac{\partial f}{\partial
v}-f(v)\frac{\partial f}{\partial v^{\prime}}\right)\text{d}v^{\prime},$ (6)
$\displaystyle\frac{1}{c^{2}}\frac{\partial E}{\partial t}=\text{curl
}B-\mu_{0}q\int vf\text{d}v,$ (7) $\displaystyle\frac{\partial B}{\partial
t}=-\text{curl }E.$ (8)

Here,
$f(x,v,t):\mathbb{R}^{3}\times\mathbb{R}^{3}\times[0,\infty)\to[0,\infty)$ is
a non-negative scalar function representing the single-particle distribution
function, and $(x,v)\in\mathbb{R}^{3}\times\mathbb{R}^{3}$ are the position
and velocity coordinates. The electric and magnetic fields are given by
$E(x,t):\mathbb{R}^{3}\times[0,\infty)\to\mathbb{R}^{3}$ and
$B(x,t):\mathbb{R}^{3}\times[0,\infty)\to\mathbb{R}^{3}$, respectively. The
tensor $U(v-v^{\prime})$ is a scaled orthogonal projector onto the orthogonal
subspace of the vector $v-v^{\prime}$, where the components are given by:

$U_{ij}(v-v^{\prime})=\frac{1}{|v-v^{\prime}|^{3}}\left(\delta_{ij}-\frac{(v_{i}-v^{\prime}_{i})(v_{j}-v^{\prime}_{j})}{|v-v^{\prime}|^{2}}\right).$
(9)

The following two equations are automatically satisfied by solutions to the
above system for all time:

$\displaystyle\text{div }E=\frac{q}{\varepsilon_{0}}\int f\text{d}v,$ (10)
$\displaystyle\text{div }B=0,$ (11)

the latter assuming that it holds for the initial conditions.

The Vlasov-Maxwell system has a well-known Hamiltonian formulation (Morrison,
1980; Marsden & Weinstein, 1982; Morrison, 1982), and so its coupling to
dissipation in the form of the Landau operator can be formulated via the
metriplectic approach.

###  2.4 The Landau Operator

In this work, we consider purely the dissipative dynamics described by the
collision operator, that is

$\partial_{t}f=C[f],$ (12)

omitting the Hamiltonian part of the system. In the metriplectic framework,
the collision operator $C[f]$ can be described as

$C[f]=-(f,\mathcal{S}),$ (13)

where the entropy functional $\mathcal{S}$ is given by:

$\mathcal{S}=-\int f\log f\text{d}v.$ (14)

The metric bracket for the Landau collision operator is given by

$(\mathcal{A},\mathcal{B})=\frac{1}{2}\int\int\left(\frac{\partial\mathcal{A}_{f}}{\partial
v}-\frac{\partial\mathcal{A}_{f^{\prime}}}{\partial
v^{\prime}}\right)f(v)U(v,v^{\prime})f(v^{\prime})\left(\frac{\partial\mathcal{B}_{f}}{\partial
v}-\frac{\partial\mathcal{B}_{f^{\prime}}}{\partial
v^{\prime}}\right)dvdv^{\prime},$ (15)

with $\mathcal{A},\mathcal{B}$ functionals, and $\mathcal{A}_{f}$ and
$\mathcal{A}_{f^{\prime}}$ denoting the functional derivatives of
$\mathcal{A}$ with respect to $f(v)$ and $f(v^{\prime})$, respectively. By the
chosen sign conventions for the entropy $\mathcal{S}$ and the bracket, entropy
is dissipated in time:

$\frac{d}{dt}\mathcal{S}=-(\mathcal{S},\mathcal{S})\leq 0.$ (16)

There are three physically important invariants in the Vlasov system that we
wish to also conserve in the metric bracket – mass, momentum, and energy:

$\mathcal{M}=m\int f\text{d}v,\hskip 18.00003pt\mathcal{P}=m\int
vf\text{d}v,\hskip 18.00003pt\mathcal{E}=m\int|v|^{2}f\text{d}v.$ (17)

We note that in the case of the full Vlasov–Maxwell–Landau system, the
conserved quantities shown above are slightly modified to take the
electromagnetic contributions into account.

The quantities in (17) are Casimir invariants of the metric bracket, i.e. they
satisfy:

$(\mathcal{A},\mathcal{M})=0,\hskip
18.00003pt(\mathcal{A},\mathcal{P})=0,\hskip
18.00003pt(\mathcal{A},\mathcal{E})=0,$ (18)

where $\mathcal{A}$ is any functional.

The equilibrium state of the system can then be obtained by applying the
energy-Casimir principle (4). Here, the equilibrium state must satisfy:

$\bigg{[}\frac{\delta\mathcal{S}}{\delta
f}+\lambda_{\mathcal{M}}\frac{\delta\mathcal{M}}{\delta
f}+\lambda_{\mathcal{P}}\cdot\frac{\delta\mathcal{P}}{\delta
f}+\lambda_{\mathcal{E}}\frac{\delta\mathcal{E}}{\delta
f}\bigg{]}\bigg{|}_{f=f_{eq}}=0.$ (19)

Using Equations (14) and (17), the above condition becomes:

$-\left(1+\log f_{eq}\right)+\lambda_{\mathcal{M}}m+\lambda_{\mathcal{P}}\cdot
mv+\lambda_{\mathcal{E}}m|v|^{2}=0,$ (20)

and so equilibrium distribution functions are of the form:

$\displaystyle f_{eq}=\exp\big{(}-(\lambda_{M}m+2\lambda_{P}\cdot
mv+\lambda_{E}m|v|^{2})\big{)}.$ (21)

Thus, we have that equilibrium distribution functions of the Landau operator
are Maxwellians, with mean, standard deviation, and normalisation determined
by initial conditions.

##  3 Semi-discretisation

We discretise the Landau collision operator by discretising its metriplectic
formulation using particles. This is to say that we do not discretise the
dynamical equations, but instead we discretise the entropy and metric bracket
and obtain discrete equations via the metriplectic approach detailed in
Section 2. This approach consists of the following main steps: (a)
discretising the distribution function, (b) discretising functionals, and (c)
discretising functional derivatives, all of which will be detailed in the
following.

###  3.1 Approximating the distribution function

The distribution function is represented by a collection of $N$ macro-
particles (Klimontovich representation):

$f_{p}(v,t)=\sum_{\alpha=1}^{N}w_{\alpha}\delta(v-v_{\alpha}(t)),$ (22)

where $v_{\alpha}(t)\in\mathbb{R}^{d}$ is the velocity of the $\alpha$-th
particle, $w_{\alpha}$ is its weight, and $d\in\\{1,2,3\\}$ is the number of
velocity dimensions being considered.

This representation of the distribution function, however, does not allow for
an immediate discretisation of the entropy functional (14) as the expression
resulting from evaluating entropy $\mathcal{S}$ on the particle distribution
$f_{p}$ is not well-defined. To overcome this, we use the approach of
Jeyakumar et al. (2023) and adopt a projected distribution function as
follows:

$f_{s}(v,t)=\sum_{i}\varphi_{i}(v)f_{i},$ (23)

where $\\{\varphi_{i}(v)\\}_{i=1}^{M}$ are a set of differentiable basis
functions with local/compact support in $\mathbb{R}^{d}$ and
$\\{f_{i}(v_{1}(t),\ldots,v_{N}(t))\\}_{i=1}^{M}$ are the coefficients of the
projected distribution function in this basis.

To obtain the coefficients $\\{f_{i}\\}_{i=1}^{M}$, we relate the two
representations of the distribution function by $L^{2}$ projection onto the
basis $\\{\varphi_{i}\\}_{i=1}^{M}$, requiring weak equality:

$\displaystyle\left<f_{s},\varphi_{j}\right>=\left<f_{p},\varphi_{j}\right>,\quad
j=1,\ldots,M,$ (24)

where

$\displaystyle\left<f_{s},\varphi_{j}\right>$ $\displaystyle=\int
f_{s}(v,t)\varphi_{j}(v)dv=\sum_{i=1}^{M}\,\mathbb{M}_{ji}f_{i},$
$\displaystyle\left<f_{p},\varphi_{j}\right>$ $\displaystyle=\int
f_{p}(v,t)\varphi_{j}(v)dv=\sum_{\alpha:v_{\alpha}\in\Omega_{\varphi_{j}}}w_{\alpha}\varphi_{j}(v_{\alpha}),$

and $\mathbb{M}_{ij}=\int\varphi_{i}(v)\varphi_{j}(v)\text{d}v$ are the
elements of the mass matrix $\mathbb{M}$ of the basis. Here, the notation
$\sum_{\alpha:v_{\alpha}\in\Omega_{\varphi_{j}}}$ signifies that we only sum
over the particles which are in the support of the basis function
$\varphi_{j}$. Thus, the coefficients of the projected distribution function
are computed as follows:

$f_{i}(v_{1}(t),\ldots,v_{N}(t))=\sum_{k=1}^{M}\mathbb{M}_{ik}^{-1}\sum_{\alpha:v_{\alpha}\in\Omega_{\varphi_{k}}}w_{\alpha}\varphi_{k}(v_{\alpha}).$
(25)

As functions of the particle velocities, these coefficients can be
differentiated with respect to individual particle velocities, $v_{\beta}$:

$\displaystyle\frac{\partial f_{i}}{\partial
v_{\beta}}=\sum_{k:v_{\beta}\in\Omega_{\varphi_{k}}}\mathbb{M}_{ik}^{-1}w_{\beta}\nabla\varphi_{k}(v_{\beta}).$
(26)

While it is not possible to evaluate the entropy (14) on the particle
distribution function $f_{p}$, evaluation on the projected distribution
function $f_{s}$ poses no issues. This second representation of the
distribution function can thus be seen as a method of regularisation in a
similar vein as previous works that used radial basis functions and similar
methods to this end (Carrillo et al., 2020; Hirvijoki, 2021).

###  3.2 Discrete functional derivatives

Following the metriplectic formulation, we will discretise the entropy
functional and bracket separately, and then apply the discrete bracket to the
discrete entropy to obtain the semi-discrete equations of motion. Discretising
the metric bracket for the Landau operator (15) requires the discretisation of
the functional derivatives appearing therein. This is facilitated by the
observation that functionals of the distribution function, $\mathcal{A}[f]$,
can be considered as functions of the particle velocities
$V=\\{v_{\alpha}\\}_{\alpha=1}^{N}$, upon evaluation on the particle
distribution $f_{p}$:

$\mathcal{A}[f_{p}]=A(V).$ (27)

By considering variations of both sides of this expression, we obtain an
approximation of the functional derivatives of $\mathcal{A}[f_{p}]$ in terms
of the function $A$ as follows111We note that, generally speaking, this
approach to discretising the functional derivatives requires the functional
being differentiated to be linear in the distribution function $f$. While, in
theory, this may seem very restrictive, in practise, this choice has no impact
on the type of “discrete functionals” on which the resulting bracket may be
evaluated. Regardless of this ansatz, the discrete bracket will be applicable
to any function of the discrete degrees of freedom, no matter if it depends
linearly or nonlinearly on the degrees of freedom, and independently from how
it is constructed, e.g. by evaluating a continuous functional on the particle
distribution, $f_{p}$, or in some other way.. Starting with the left-hand
side, variation and integration by parts gives,

$\displaystyle\delta\mathcal{A}[f_{p}]$
$\displaystyle=\int\frac{\delta\mathcal{A}}{\delta f}\,\delta f_{p}\,dv$
$\displaystyle=-\sum\limits_{a=1}^{N_{p}}w_{a}\int\dfrac{\delta\mathcal{A}}{\delta
f}\,\nabla_{v}\delta\big{(}v-v_{a}(t)\big{)}\cdot\delta v_{a}\,dv$
$\displaystyle=\sum\limits_{a=1}^{N_{p}}w_{a}\,\nabla_{v}\frac{\delta\mathcal{A}}{\delta
f}\bigg{|}_{v=v_{a}}\cdot\delta v_{a}.$ (28)

Upon equating this expression with the variation of the right-hand side of
(27), namely

$\displaystyle\delta A(V)$
$\displaystyle=\sum\limits_{a=1}^{N_{p}}\frac{\partial A}{\partial
v_{a}}\,\delta v_{a},$ (29)

we obtain:

$\frac{\partial A}{\partial v_{\alpha}}=w_{\alpha}\frac{\partial}{\partial
v}\left.\frac{\delta\mathcal{A}}{\delta f}\right|_{v=v_{\alpha}}.$ (30)

Evaluating the bracket of Equation (15) on the particle distribution from
Equation (22), we obtain the following bracket on functionals:

$(\mathcal{A},\mathcal{B})=\frac{1}{2}\sum_{\alpha,\beta}w_{\alpha}w_{\beta}\int\left(\frac{\partial\mathcal{A}_{f}}{\partial
v}-\frac{\partial\mathcal{A}_{f^{\prime}}}{\partial
v^{\prime}}\right)\delta(v-v_{\alpha})U(v,v^{\prime})\delta(v^{\prime}-v_{\beta})\left(\frac{\partial\mathcal{B}_{f}}{\partial
v}-\frac{\partial\mathcal{B}_{f^{\prime}}}{\partial
v^{\prime}}\right)dvdv^{\prime},$ (31)

for some functionals $\mathcal{A},\mathcal{B}$. Using relation (30) to replace
the functional derivatives $\mathcal{A}_{f}$, etc., we obtain the discrete
bracket on functions:

$(A,B)_{h}=\frac{1}{2}\sum_{\alpha,\beta}w_{\alpha}w_{\beta}\left(\frac{1}{w_{\alpha}}\frac{\partial
A}{\partial v_{\alpha}}-\frac{1}{w_{\beta}}\frac{\partial A}{\partial
v_{\beta}}\right)U(v_{\alpha},v_{\beta})\left(\frac{1}{w_{\alpha}}\frac{\partial
B}{\partial v_{\alpha}}-\frac{1}{w_{\beta}}\frac{\partial B}{\partial
v_{\beta}}\right),$ (32)

where $A$, $B$ are functions of the particle velocities
$V=\\{v_{\alpha}\\}_{\alpha=1}^{N}$.

###  3.3 Discrete entropy

The remaining ingredient required for constructing the discrete collision
operator is the entropy. The entropy functional (14) is not amenable to plain
evaluation on the particle distribution $f_{p}$ as the resulting expression is
not well-defined. Thus, we evaluate the entropy functional (14) on the
projected distribution function $f_{s}$ instead, such that the discrete
entropy becomes:

$\displaystyle S_{h}$ $\displaystyle=-\int
f_{s}\log{f_{s}}\,dv=-\int\left(\sum_{i}\varphi_{i}(v)f_{i}\right)\log\left(\sum_{j}\varphi_{j}(v)f_{j}\right)\text{d}v.$
(33)

Using this expression for entropy, we now compute its derivative with respect
to the degrees of freedom (this will be required to apply the metric bracket
to the entropy):

$\displaystyle\frac{\partial S_{h}}{\partial
v_{\alpha}}=-\int\left[\sum_{i}\frac{\partial f_{i}}{\partial
v_{\alpha}}\varphi_{i}(v)\log\left(\sum_{j}\varphi_{j}(v)f_{j}\right)+\frac{\sum_{i}\varphi_{i}(v)f_{i}}{\sum_{j}\varphi_{j}(v)f_{j}}\sum_{j}\frac{\partial
f_{j}}{\partial v_{\alpha}}\varphi_{j}(v)\right]\text{d}v,$ (34)

where the derivative $\partial f_{i}/\partial v_{\alpha}$ is given by Equation
(26). Simplifying the second term and using Equation (26) gives us the
following:

$\displaystyle\frac{\partial S_{h}}{\partial v_{\alpha}}$
$\displaystyle=-\int\left[\sum_{i}\frac{\partial f_{i}}{\partial
v_{\alpha}}\varphi_{i}(v)\log\left(\sum_{j}\varphi_{j}(v)f_{j}\right)+\sum_{j}\frac{\partial
f_{j}}{\partial v_{\alpha}}\varphi_{j}(v)\right]\text{d}v,$
$\displaystyle=-\int\left[\sum_{i,k}\mathbb{M}^{-1}_{ik}w_{\alpha}\nabla\varphi_{k}(v_{\alpha})\varphi_{i}(v)\log\left(\sum_{j,l}\varphi_{j}\mathbb{M}^{-1}_{jl}\sum_{\beta}w_{\beta}\varphi_{l}(v_{\beta})\right)\right.$
$\displaystyle\hskip
54.00009pt\left.+\sum_{j,l}\mathbb{M}^{-1}_{jl}w_{\alpha}\nabla\varphi_{l}(v_{\alpha})\varphi_{j}(v)\right]\text{d}v.$
(35)

We define the following for convenience:

$\displaystyle\mathbb{L}_{k}=\sum_{i}\mathbb{M}^{-1}_{ik}\int\varphi_{i}(v)\left[1+\log\left(\sum_{j}f_{j}\varphi_{j}(v)\right)\right]\text{d}v.$
(36)

Relabelling indices in (35) and using (36), we obtain our final expression for
the derivative of the discrete entropy with respect to the degrees of freedom:

$\frac{\partial S_{h}}{\partial
v_{\alpha}}=-\sum_{k}w_{\alpha}\mathbb{L}_{k}\nabla\varphi_{k}(v_{\alpha}).$
(37)

With this, we can now compute the collision contribution for a given particle,
$\dot{v}_{\gamma}$, through the metric bracket.

###  3.4 Discrete collision operator

The discrete collisional dynamics are given by

$\displaystyle\dot{v}_{\gamma}$ $\displaystyle=-(v_{\gamma},S_{h})_{h}.$ (38)

Upon inserting the expression of the discrete bracket (32) and the derivative
of the entropy (37), we obtain:

$\displaystyle\dot{v}_{\gamma}$
$\displaystyle=\frac{1}{2}\sum_{\alpha,\beta}w_{\alpha}w_{\beta}\left(\frac{1}{w_{\alpha}}\delta_{\alpha\gamma}-\frac{1}{w_{\beta}}\delta_{\beta\gamma}\right)U(v_{\alpha},v_{\beta})\sum_{k}\mathbb{L}_{k}(\nabla\varphi_{k}(v_{\alpha})-\nabla\varphi_{k}(v_{\beta})),$
(39)
$\displaystyle=\frac{1}{2}\left[\sum_{\beta}w_{\beta}U(v_{\gamma},v_{\beta})\sum_{k}\mathbb{L}_{k}(\nabla\varphi_{k}(v_{\gamma})-\nabla\varphi_{k}(v_{\beta}))\right.$
$\displaystyle\hskip
54.00009pt-\left.\sum_{\alpha}w_{\alpha}U(v_{\alpha},v_{\gamma})\sum_{k}\mathbb{L}_{k}(\nabla\varphi_{k}(v_{\alpha})-\nabla\varphi_{k}(v_{\gamma}))\right],$
(40)
$\displaystyle=\sum_{\alpha}w_{\alpha}U(v_{\gamma},v_{\alpha})\sum_{k}\mathbb{L}_{k}(\nabla\varphi_{k}(v_{\gamma})-\nabla\varphi_{k}(v_{\alpha})).$
(41)

Here, we obtain Equation (40) through separating the two Kronecker delta terms
in (39) and contracting the appropriate indices. Equation (41) then comes from
relabelling indices and recognising that the two sums in (40) are in fact
identical by symmetry of the kernel $U(\cdot,\cdot)$. Equation (41) is the
final semi-discretised system of equations.

For the purpose of discretising the system in time, it will be useful to
rewrite the semi-discretised system of equations as a gradient system. To
obtain this form, we rearrange Equation (39) as follows with the use of (37):

$\dot{v}_{\gamma}=\mathbb{G}_{\gamma\sigma}\frac{\partial S_{h}}{\partial
v_{\sigma}}=-\sum_{\beta}\left(\frac{w_{\beta}}{w_{\gamma}}\delta_{\gamma\sigma}-\delta_{\beta\sigma}\right)U(v_{\gamma},v_{\beta})\frac{\partial
S_{h}}{\partial v_{\sigma}}.$ (42)

Here, the operator $\mathbb{G}$ can be expressed as:

$\mathbb{G}(v_{1},\dots,v_{n})=\begin{bmatrix}-\sum\limits_{\beta\neq
1}U(v_{1},v_{\beta})\frac{w_{\beta}}{w_{1}}&&U(v_{1},v_{2})&&U(v_{1},v_{3})&&\dots\\\
U(v_{2},v_{1})&&-\sum\limits_{\beta\neq
2}U(v_{2},v_{\beta})\frac{w_{\beta}}{w_{2}}&&U(v_{2},v_{3})&&\dots\\\
\vdots&&\hskip 63.0001pt\ddots&&&&\end{bmatrix}.$ (43)

The matrix $\mathbb{G}$ is a $Nd\times Nd$ block matrix, where $d$ is the
number of velocity-space dimensions and $N$ the number of particles, and each
$d\times d$ block is given by

$\displaystyle\mathbb{G}_{\gamma\sigma}=\begin{cases}-\sum\limits_{\beta\neq\gamma}U(v_{\gamma},v_{\beta})\frac{w_{\beta}}{w_{\gamma}}\qquad&\text{if}\;\gamma=\sigma,\\\
U(v_{\gamma},v_{\sigma})&\text{else}.\end{cases}$ (44)

We observe by the symmetry of the kernel $U(\cdot,\cdot)$ that the operator
$\mathbb{G}$ is also symmetric. It is also negative semi-definite. To see
this, take a vector $p\in\mathbb{R}^{dN}$ and interpret it as structured into
$N$ blocks $p_{\gamma}$ of size $d$, in analogy to the structure of
$\mathbb{G}$, so that each block $p_{\gamma}\in\mathbb{R}^{d}$. Multiplying
the operator $\mathbb{G}$ from both sides with $p$, we have:

$\displaystyle p^{T}_{\gamma}\mathbb{G}_{\gamma\sigma}p_{\sigma}$
$\displaystyle=-\sum_{\gamma,\sigma=1}^{N}\sum_{\beta=1}^{N}p^{T}_{\gamma}\left(\frac{w_{\beta}}{w_{\gamma}}\delta_{\gamma\sigma}-\delta_{\beta\sigma}\right)U(v_{\gamma},v_{\beta})p_{\sigma},$
$\displaystyle=-\sum_{\sigma=1}^{N}\sum_{\beta=1}^{N}\frac{w_{\beta}}{w_{\sigma}}p^{T}_{\sigma}U(v_{\sigma},v_{\beta})p_{\sigma}+\sum_{\gamma=1}^{N}\sum_{\beta=1}^{N}p^{T}_{\gamma}U(v_{\gamma},v_{\beta})p_{\beta},$
$\displaystyle=-\sum_{\sigma=1}^{N}\sum_{\beta=1}^{N}\left(\frac{w_{\beta}}{w_{\sigma}}p^{T}_{\sigma}U(v_{\sigma},v_{\beta})p_{\sigma}-p^{T}_{\sigma}U(v_{\sigma},v_{\beta})p_{\beta}\right).$
(45)

Here, we have expanded out the two terms, contracted the appropriate indices
with the delta functions, and relabelled the resulting indices in the second
sum. Factoring out the left-multiplying terms, we obtain:

$\displaystyle p^{T}_{\gamma}\mathbb{G}_{\gamma\sigma}p_{\sigma}$
$\displaystyle=-\sum_{\sigma=1}^{N}\sum_{\beta=1}^{N}\frac{p^{T}_{\sigma}}{w_{\sigma}}U(v_{\sigma},v_{\beta})\left(w_{\beta}p_{\sigma}-w_{\sigma}p_{\beta}\right),$
(46)
$\displaystyle=-\frac{1}{2}\sum_{\sigma=1}^{N}\sum_{\beta=1}^{N}\left(w_{\beta}p_{\sigma}-w_{\sigma}p_{\beta}\right)^{T}\frac{U(v_{\sigma},v_{\beta})}{w_{\sigma}w_{\beta}}\left(w_{\beta}p_{\sigma}-w_{\sigma}p_{\beta}\right)\leq
0.$ (47)

The final equation is obtained by realising that swapping indices in Equation
(46) results in the same term up to a minus sign, and further re-arranging
results in the form shown above. As the matrix $U$ is a projector and
therefore positive semi-definite, and is being multiplied from the left and
right by the same vector, we have that every term in the sum is non-positive.
Thus the entire sum is non-positive and the operator $\mathbb{G}$ is negative
semi-definite.

###  3.5 Conservation of momentum and energy

We now demonstrate that the semi-discrete equations of motion preserve total
momentum and energy. We note that mass is trivially preserved by this scheme,
due to the particle weights being held constant. We start by looking at the
evolution of the total particle momentum,
$P=\sum_{\gamma}w_{\gamma}v_{\gamma}$:

$\displaystyle\frac{\text{d}}{\text{d}t}\sum_{\gamma}w_{\gamma}v_{\gamma}$
$\displaystyle=\sum_{\gamma}w_{\gamma}\dot{v}_{\gamma}=\sum_{\gamma}w_{\gamma}\mathbb{G}_{\gamma\sigma}\frac{\partial
S_{h}}{\partial v_{\sigma}},$ (48)

which gives us that $\text{d}P/\text{d}t=0$, as the product
$w_{\gamma}\mathbb{G}_{\gamma\sigma}$ is equal to a vector which is zero
element-wise. Thus, the semi-discretisation preserves total particle momentum
over time.

To demonstrate energy conservation, it will be useful to consider the action
of the operator $U(\cdot,\cdot)$. We note that $U(v,v^{\prime})$ is an
orthogonal projector onto the subspace orthogonal to the vector
$v-v^{\prime}$. It follows that
$U(v,v^{\prime})(v-v^{\prime})=U(v^{\prime},v)(v-v^{\prime})=\vec{0}$. We now
follow the same process as before to check conservation, where the energy
density is given by $E=\sum_{\gamma}w_{\gamma}v_{\gamma}^{2}/2$:

$\displaystyle\frac{\text{d}E}{\text{d}t}$
$\displaystyle=\sum_{\gamma}w_{\gamma}v_{\gamma}\dot{v}_{\gamma}=\sum_{\gamma}w_{\gamma}v_{\gamma}\mathbb{G}_{\gamma\sigma}\frac{\partial
S_{h}}{\partial v_{\sigma}}.$ (49)

Without loss of generality, we choose $\sigma=1$. If we consider the
contraction $\sum_{\gamma}w_{\gamma}v_{\gamma}\mathbb{G}_{\gamma 1}$, we
obtain:

$\displaystyle\sum_{\gamma}w_{\gamma}v_{\gamma}\mathbb{G}_{\gamma 1}$
$\displaystyle=-\sum_{\beta\neq
1}U(v_{1},v_{\beta})w_{\beta}v_{1}+U(v_{1},v_{2})w_{2}v_{2}$
$\displaystyle\hskip
54.00009pt+U(v_{1},v_{3})w_{3}v_{3}+\dots+U(v_{1},v_{n})w_{n}v_{n}$
$\displaystyle=\sum_{\beta\neq 1}w_{\beta}U(v_{1},v_{\beta})(v_{\beta}-v_{1})$
$\displaystyle=0,$ (50)

where the last line follows by the property of $U(\cdot,\cdot)$ being an
orthogonal projector. As this argument holds for all $\sigma$, it follows that
$\text{d}E/\text{d}t=0$ and the semi-discretisation preserves the energy.

###  3.6 H-Theorem

The semi-discrete collisional dynamics (38) satisfies a H-theorem, such that
entropy is monotonically dissipated. This can be seen by considering the time
evolution of the discretised entropy:

$\displaystyle\frac{d}{dt}S_{h}$ $\displaystyle=-(S_{h},S_{h})_{h}$
$\displaystyle=-\frac{1}{2}\sum_{\alpha,\beta}w_{\alpha}w_{\beta}\left(\frac{1}{w_{\alpha}}\frac{\partial
S_{h}}{\partial v_{\alpha}}-\frac{1}{w_{\beta}}\frac{\partial S_{h}}{\partial
v_{\beta}}\right)^{T}U(v_{\alpha},v_{\beta})\left(\frac{1}{w_{\alpha}}\frac{\partial
S_{h}}{\partial v_{\alpha}}-\frac{1}{w_{\beta}}\frac{\partial S_{h}}{\partial
v_{\beta}}\right),$ (51)

where $(\cdot,\cdot)_{h}$ is the semi-discrete metric bracket given in
Equation (32).

As the matrix $U$ is positive semi-definite, and the terms multiplying it in
the above equation are identical for a given $\alpha$ and $\beta$, and
$w_{\alpha},w_{\beta}\geq 0$, each term in the sum is non-negative. As such,
we have that $dS_{h}/dt\leq 0$ and so the discretised entropy is monotonically
dissipated.

We can also compute the discrete equilibrium state by using the energy-Casimir
principle, as before. The energy-Casimir principle for the discrete system
reads:

$\delta S_{h}(f_{s,eq})+\delta\sum_{i}\lambda_{i}C_{i}(f_{s,eq})=0,$ (52)

where $S_{h}(f_{s,eq})$ is the discrete entropy (33), evaluated at the
discrete equilibrium distribution $f_{s,eq}$, that is the spline
representation of the particles with equilibrated velocities $v_{\alpha,eq}$.
The $C_{i}$’s are the discrete Casimir invariants, i.e. the discrete mass,
momentum, and energy. Thus, the above equation becomes:

$\left[\delta
S_{h}+\lambda_{M_{h}}\delta\left(\sum_{\alpha}w_{\alpha}\right)+\lambda_{P_{h}}\delta\left(\sum_{\alpha}w_{\alpha}v_{\alpha}\right)+\lambda_{E_{h}}\delta\left(\sum_{\alpha}w_{\alpha}v_{\alpha}^{2}\right)\right]_{v_{\alpha}=v_{\alpha,eq}}=0,$
(53)

where $\lambda_{M_{h}}$, $\lambda_{P_{h}}$, and $\lambda_{E_{h}}$ are the
corresponding Lagrange multipliers for the discrete quantities. Computing the
variations in this equation, we have:

$\left.\sum_{\alpha}\left[\frac{\partial S_{h}}{\partial
v_{\alpha}}+\lambda_{P_{h}}w_{\alpha}+\lambda_{E_{h}}w_{\alpha}v_{\alpha}\right]\delta
v_{\alpha}\right|_{v_{\alpha}=v_{\alpha,eq}}=0,$ (54)

where the quantity $\partial S_{h}/\partial v_{\alpha}$ is given by Equation
(37), and the term corresponding to the mass does not appear as it does not
depend on the particle velocities, thus having zero variation. Thus, we find
that an equilibrium distribution’s degrees of freedom, i.e. the particle
velocities $\\{v_{\alpha,eq}\\}$ must satisfy:

$\left[\frac{\partial S_{h}}{\partial
v_{\alpha}}+\lambda_{P_{h}}w_{\alpha}+\lambda_{E_{h}}w_{\alpha}v_{\alpha}\right]_{v_{\alpha}=v_{\alpha,eq}}=0,$
(55)

for the appropriate Lagrange multipliers based on the chosen initial
conditions.

##  4 Time discretisation

In order to discretise the system of ODEs in time, we use the discrete
gradient method (Quispel & Turner, 1996; McLachlan et al., 1999). It is an
integral preserving method, initially developed for Hamiltonian systems and
others with first integrals, which can also be adopted to dissipative systems
as shown below.

To illustrate, we begin with a system of ODEs in gradient form,
$\dot{z}=G(z)\nabla S(z)$, where $z\in\mathbb{R}^{n}$. Here, $G(z)$ can be a
skew-symmetric matrix in the case of a Hamiltonian system, or a negative semi-
definite symmetric matrix in the case of a dissipative one. The system of ODEs
is then approximately solved by the following equation:

$\frac{z_{n+1}-z_{n}}{\Delta
t}=\tilde{G}(z_{n},z_{n+1})\bar{\nabla}S(z_{n},z_{n+1}),$ (56)

where $\bar{\nabla}S(z_{n},z_{n+1})$ is known as the discrete gradient, and
$\tilde{G}(z_{n},z_{n+1})$ is a matrix such that $\tilde{G}(z_{n},z_{n+1})\to
G(z_{n})$ in the limit $z_{n+1}\to z_{n}$ as $\Delta t\to 0$. The discrete
gradient, $\bar{\nabla}S(z_{n},z_{n+1})$, satisfies the following two
properties:

$\displaystyle(z_{n+1}-z_{n})\cdot\bar{\nabla}S(z_{n},z_{n+1})$
$\displaystyle=S(z_{n+1})-S(z_{n}),$ (57a)
$\displaystyle\bar{\nabla}S(z_{n},z_{n})$ $\displaystyle=\nabla S(z_{n}).$
(57b)

One such discrete gradient is the one by Gonzalez (Gonzalez, 1996), also known
as the midpoint discrete gradient:

$\bar{\nabla}S(z_{n},z_{n+1})=\nabla S(z_{n+1/2})\\\
+(z_{n+1}-z_{n})\frac{S(z_{n+1})-S(z_{n})-(z_{n+1}-z_{n})\cdot\nabla
S(z_{n+1/2})}{|z_{n+1}-z_{n}|^{2}},$ (58)

where $z_{n+1/2}=(z_{n}+z_{n+1})/2$ is the midpoint.

In the case of the Landau operator, Equation (42) is the gradient form of the
semi-discrete system of ODE’s. Applying the discrete gradient approximation of
Equation (56), we have:

$\displaystyle\frac{v_{\gamma}^{n+1}-v_{\gamma}^{n}}{\Delta
t}=\tilde{\mathbb{G}}_{\gamma\sigma}(v_{n},v_{n+1})\bar{\nabla}_{\sigma}S_{h}(v_{n},v_{n+1}),$
(59)

where $\bar{\nabla}_{\sigma}S_{h}(v_{n},v_{n+1})$ represents the discrete
gradient of $\partial S_{h}/\partial v_{\sigma}$, satisfying the properties in
(57), and $v_{n}$ represents a vector of all particle velocities at time $n$.
Here, we take the midpoint approximation for the matrix
$\tilde{\mathbb{G}}_{\gamma\sigma}$, i.e.:

$\tilde{\mathbb{G}}_{\gamma\sigma}(v_{n},v_{n+1})=\mathbb{G}_{\gamma\sigma}\left(\frac{v_{n}+v_{n+1}}{2}\right).$
(60)

In the treatment of the full Vlasov–Maxwell-Landau system, this approach for
time discretisation of the dissipative part can immediately be adopted for the
conservative part as well, possibly in conjunction with some kind of operator
splitting.

###  4.1 Conservation of momentum and energy

Under this scheme, we can show that momentum and energy are conserved in time
as well. We first consider the total particle momentum:

$\displaystyle\sum_{\gamma}w_{\gamma}v_{\gamma}^{n+1}-\sum_{\gamma}w_{\gamma}v_{\gamma}^{n}$
$\displaystyle=\Delta
t\sum_{\gamma}w_{\gamma}\tilde{\mathbb{G}}_{\gamma\sigma}(v_{n},v_{n+1})\bar{\nabla}_{\sigma}S_{h}(v_{n},v_{n+1})$
$\displaystyle=\Delta
t\sum_{\gamma}w_{\gamma}\mathbb{G}_{\gamma\sigma}\left(\frac{v_{n}+v_{n+1}}{2}\right)\bar{\nabla}_{\sigma}S_{h}(v_{n},v_{n+1}),$
(61)

where we have used Equation (60) to obtain the second line. By the structure
of the matrix $\mathbb{G}_{\gamma\sigma}$ as shown in Equation (43), we have
that the sum $\sum_{\gamma}w_{\gamma}\mathbb{G}_{\gamma\sigma}(\cdot)=0$, and
so:

$\Delta
P=\sum_{\gamma}w_{\gamma}v_{\gamma}^{n+1}-\sum_{\gamma}w_{\gamma}v_{\gamma}^{n}=0.$
(62)

Thus, total particle momentum is conserved exactly in time.

Now, we consider the total particle energy. We have:

$\displaystyle\frac{1}{2}\sum_{\gamma}w_{\gamma}\left(v_{\gamma}^{n+1}\right)^{2}$
$\displaystyle=\frac{1}{2}\sum_{\gamma}w_{\gamma}\left[v_{\gamma}^{n}+\Delta
t\tilde{\mathbb{G}}_{\gamma\sigma}(v_{n},v_{n+1})\bar{\nabla}_{\sigma}S_{h}(v_{n},v_{n+1})\right]^{2},$
$\displaystyle=\frac{1}{2}\sum_{\gamma}w_{\gamma}\left(v_{\gamma}^{n}\right)^{2}+\Delta
t\sum_{\gamma}w_{\gamma}v_{\gamma}^{n}\tilde{\mathbb{G}}_{\gamma\sigma}(v_{n},v_{n+1})\bar{\nabla}_{\sigma}S_{h}(v_{n},v_{n+1})$
$\displaystyle\hskip 27.00005pt+\frac{1}{2}\Delta
t^{2}\sum_{\gamma}w_{\gamma}\left(\tilde{\mathbb{G}}_{\gamma\sigma}(v_{n},v_{n+1})\bar{\nabla}_{\sigma}S_{h}(v_{n},v_{n+1})\right)^{2}.$
(63)

Rewriting the second term of the above equation using the discrete gradient
update rule rearranged as $v_{\gamma}^{n}=v_{\gamma}^{n+1}-\Delta
t\tilde{\mathbb{G}}_{\gamma\sigma}(v_{n},v_{n+1})\bar{\nabla}_{\sigma}S_{h}(v_{n},v_{n+1})$,
we have the following useful expression

$\displaystyle\Delta
t\sum_{\gamma}w_{\gamma}v_{\gamma}^{n}\tilde{\mathbb{G}}_{\gamma\sigma}(v_{n},v_{n+1})\bar{\nabla}_{\sigma}S_{h}(v_{n},v_{n+1})$
$\displaystyle\hskip 27.00005pt=\Delta
t\sum_{\gamma}w_{\gamma}\left(v_{\gamma}^{n+1}-\Delta
t\tilde{\mathbb{G}}_{\gamma\beta}(v_{n},v_{n+1})\bar{\nabla}_{\beta}S_{h}(v_{n},v_{n+1})\right)\tilde{\mathbb{G}}_{\gamma\sigma}(v_{n},v_{n+1})\bar{\nabla}_{\sigma}S_{h}(v_{n},v_{n+1}),$
$\displaystyle\hskip 27.00005pt=\Delta
t\sum_{\gamma}w_{\gamma}v_{\gamma}^{n+1}\tilde{\mathbb{G}}_{\gamma\sigma}(v_{n},v_{n+1})\bar{\nabla}_{\sigma}S_{h}(v_{n},v_{n+1})$
$\displaystyle\hskip 54.00009pt-\Delta
t^{2}\sum_{\gamma}w_{\gamma}\left(\tilde{\mathbb{G}}_{\gamma\sigma}(v_{n},v_{n+1})\bar{\nabla}_{\sigma}S_{h}(v_{n},v_{n+1})\right)^{2}.$
(64)

Now, we rewrite Equation (63) by splitting up the second term into two halves.
Using Equation (4.1) to replace one of those terms, we have:

$\displaystyle\Delta E$
$\displaystyle=\frac{1}{2}\sum_{\gamma}w_{\gamma}\left(v_{\gamma}^{n+1}\right)^{2}-\frac{1}{2}\sum_{\gamma}w_{\gamma}\left(v_{\gamma}^{n}\right)^{2}$
$\displaystyle=\frac{1}{2}\Delta
t\sum_{\gamma}w_{\gamma}v_{\gamma}^{n}\tilde{\mathbb{G}}_{\gamma\sigma}(v_{n},v_{n+1})\bar{\nabla}_{\sigma}S_{h}(v_{n},v_{n+1})$
$\displaystyle\hskip 27.00005pt+\frac{1}{2}\Delta
t\sum_{\gamma}w_{\gamma}v_{\gamma}^{n+1}\tilde{\mathbb{G}}_{\gamma\sigma}(v_{n},v_{n+1})\bar{\nabla}_{\sigma}S_{h}(v_{n},v_{n+1})$
$\displaystyle\hskip 27.00005pt-\frac{1}{2}\Delta
t^{2}\sum_{\gamma}w_{\gamma}\left(\tilde{\mathbb{G}}_{\gamma\sigma}(v_{n},v_{n+1})\bar{\nabla}_{\sigma}S_{h}(v_{n},v_{n+1})\right)^{2}$
$\displaystyle\hskip 27.00005pt+\frac{1}{2}\Delta
t^{2}\sum_{\gamma}w_{\gamma}\left(\tilde{\mathbb{G}}_{\gamma\sigma}(v_{n},v_{n+1})\bar{\nabla}_{\sigma}S_{h}(v_{n},v_{n+1})\right)^{2},$
(65)

where the last two terms in the above equation cancel exactly. Factorising the
first two terms and then using the definition of
$\tilde{\mathbb{G}}_{\gamma\sigma}(\cdot,\cdot)$, we have:

$\displaystyle\Delta E$ $\displaystyle=\Delta
t\sum_{\gamma}w_{\gamma}\left(\frac{v_{\gamma}^{n}+v_{\gamma}^{n+1}}{2}\right)\tilde{\mathbb{G}}_{\gamma\sigma}(v_{n},v_{n+1})\bar{\nabla}_{\sigma}S_{h}(v_{n},v_{n+1}),$
$\displaystyle=\Delta
t\sum_{\gamma}w_{\gamma}\left(\frac{v_{\gamma}^{n}+v_{\gamma}^{n+1}}{2}\right)\mathbb{G}_{\gamma\sigma}\left(\frac{v_{n}+v_{n+1}}{2}\right)\bar{\nabla}_{\sigma}S_{h}(v_{n},v_{n+1}).$
(66)

By the arguments given in the semi-discrete energy conservation proof, we have
that $\sum_{\gamma}w_{\gamma}v_{\gamma}\mathbb{G}_{\gamma\sigma}=0$, where
$v_{\gamma}$ and $\mathbb{G}_{\gamma\sigma}$ are evaluated at the same point
in time. As such, we have that
$\sum_{\gamma}w_{\gamma}[(v_{\gamma}^{n}+v_{\gamma}^{n+1})/2]\mathbb{G}_{\gamma\sigma}[(v_{\gamma}^{n}+v_{\gamma}^{n+1})/2]=0$,
and so:

$\Delta
E=\frac{1}{2}\sum_{\gamma}w_{\gamma}\left(v_{\gamma}^{n+1}\right)^{2}-\frac{1}{2}\sum_{\gamma}w_{\gamma}\left(v_{\gamma}^{n}\right)^{2}=0.$
(67)

Thus, the total particle energy is also conserved exactly in time.

###  4.2 H-theorem

We can demonstrate that entropy is also monotonically dissipated in time. From
Equation (57), we have that:

$S(v_{n+1})-S(v_{n})=(v_{\gamma}^{n+1}-v_{\gamma}^{n})\cdot\bar{\nabla}_{\gamma}S_{h}(v_{n},v_{n+1}).$
(68)

Using the discrete gradient update rule (59), this becomes:

$S(v_{n+1})-S(v_{n})=\Delta
t\bar{\nabla}_{\sigma}S_{h}(v_{n},v_{n+1})\tilde{\mathbb{G}}_{\sigma\gamma}\bar{\nabla}_{\gamma}S_{h}(v_{n},v_{n+1}).$
(69)

As the matrix $\tilde{\mathbb{G}}$ is simply $\mathbb{G}$ evaluated at the
midpoint (c.f. Equation (60)), it is also negative semi-definite. Since the
vectors multiplying $\tilde{\mathbb{G}}$ from the left and right are
identical, we have:

$S(v_{n+1})-S(v_{n})\leq 0,$ (70)

and so entropy is monotonically dissipated in time. Thus, the H-theorem is
maintained at the discrete level.

##  5 Conclusion

In this work, we have described a structure-preserving particle method for
computing the collisional dynamics described by the Landau operator, using the
metriplectic formulation. This method is based on a projection based approach,
which was introduced in Jeyakumar et al. (2023). The method’s structure-
preserving properties have been demonstrated analytically in the semi-discrete
setting (only spatial discretisation), as well as in the fully discrete
setting (including discretisation in time using discrete gradients). The
H-theorem has also been shown to hold in both settings. This method can be
coupled to an appropriate particle method for the Vlasov–Poisson or
Vlasov–Maxwell equations, and the coupling methodology and benefits will be
reported in future papers.

## References

  * Bailo et al. (2024) Bailo, Rafael, Carrillo, José A. & Hu, Jingwei 2024 The Collisional Particle-In-Cell Method for the Vlasov-Maxwell-Landau Equations , arXiv: 2401.01689.
  * Boozer & Kuo-Petravic (1981) Boozer, Allen H. & Kuo-Petravic, Gioietta 1981 Monte Carlo evaluation of transport coefficients. The Physics of Fluids 24 (5), 851–859.
  * Burby (2017) Burby, J W 2017 Finite-dimensional collisionless kinetic theory. Physics of Plasmas 24 (3), 32101.
  * Campos Pinto et al. (2022) Campos Pinto, Martin, Kormann, Katharina & Sonnendrücker, Eric 2022 Variational Framework for Structure-Preserving Electromagnetic Particle-in-Cell Methods. Journal of Scientific Computing 91 (2), 46.
  * Carrillo et al. (2020) Carrillo, Jose A., Hu, Jingwei, Wang, Li & Wu, Jeremy 2020 A particle method for the homogeneous Landau equation. Journal of Computational Physics: X 7, 100066.
  * Chen et al. (2011) Chen, G., Chacón, L. & Barnes, D.C. 2011 An energy- and charge-conserving, implicit, electrostatic particle-in-cell algorithm. Journal of Computational Physics 230 (18), 7018–7036.
  * Evstatiev & Shadwick (2013) Evstatiev, Evstati G. & Shadwick, Bradley A. 2013 Variational formulation of particle algorithms for kinetic plasma simulations. Journal of Computational Physics 245, 376–398.
  * Gonzalez (1996) Gonzalez, O. 1996 Time Integration and Discrete Hamiltonian Systems. Journal of Nonlinear Science 6 (5), 449–467.
  * Grmela (1984a) Grmela, Miroslav 1984a Bracket formulation of dissipative fluid mechanics equations. Physics Letters A 102 (8), 355–358.
  * Grmela (1984b) Grmela, Miroslav 1984b Particle and bracket formulations of kinetic equations. In Fluids and Plasmas: Geometry and Dynamics, pp. 125–132. American Mathematical Society.
  * Grmela (1985) Grmela, Miroslav 1985 Bracket formulation of dissipative time evolution equations. Physics Letters A 111 (1-2), 36–40.
  * Hirvijoki (2021) Hirvijoki, Eero 2021 Structure-preserving marker-particle discretizations of Coulomb collisions for particle-in-cell codes. Plasma Physics and Controlled Fusion 63 (4), 44003.
  * Hirvijoki & Adams (2017) Hirvijoki, E. & Adams, M. F. 2017 Conservative discretization of the Landau collision integral. Physics of Plasmas 24 (3).
  * Hirvijoki et al. (2018) Hirvijoki, Eero, Kraus, Michael & Burby, Joshua W. 2018 Metriplectic particle-in-cell integrators for the landau collision operator, arXiv: 1802.05263.
  * Jeyakumar et al. (2023) Jeyakumar, Sandra, Kraus, Michael, Hole, Matthew & Pfefferlé, David 2023 A structure-preserving particle discretisation for the Lenard-Bernstein operator, arXiv: 2309.16894.
  * Kaufman (1984) Kaufman, Allan N. 1984 Dissipative hamiltonian systems: A unifying principle. Physics Letters A 100 (8), 419–422.
  * Kaufman & Morrison (1982) Kaufman, Allan N. & Morrison, Philip J. 1982 Algebraic structure of the plasma quasilinear equations. Physics Letters A 88 (8), 405–406.
  * Kraus & Hirvijoki (2017) Kraus, Michael & Hirvijoki, Eero 2017 Metriplectic integrators for the Landau collision operator. Physics of Plasmas 24 (10).
  * Kraus et al. (2017) Kraus, Michael, Kormann, Katharina, Morrison, Philip J. & Sonnendrücker, Eric 2017 GEMPIC: Geometric electromagnetic particle-in-cell methods. Journal of Plasma Physics 83 (4).
  * Laidin & Pareschi (2024) Laidin, Tino & Pareschi, Lorenzo 2024 Conservative polynomial approximations and applications to Fokker-Planck equations , arXiv: 2402.06473.
  * Markidis & Lapenta (2011) Markidis, Stefano & Lapenta, Giovanni 2011 The energy conserving particle-in-cell method. Journal of Computational Physics 230 (18), 7037–7052.
  * Marsden & Weinstein (1982) Marsden, Jerrold E. & Weinstein, Alan 1982 The Hamiltonian structure of the Maxwell-Vlasov equations. Physica D: Nonlinear Phenomena 4 (3), 394–406.
  * McLachlan et al. (1999) McLachlan, Robert I., Quispel, G. R.W. & Robidoux, Nicolas 1999 Geometric integration using discrete gradients. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 357 (1754), 1021–1045.
  * Morrison (1980) Morrison, Philip J 1980 The Maxwell-Vlasov equations as a continuous Hamiltonian system. Physics Letters 80.
  * Morrison (1982) Morrison, Philip J. 1982 Poisson brackets for fluids and plasmas. In AIP Conference Proceedings, , vol. 88, pp. 13–46. AIP.
  * Morrison (1984a) Morrison, Philip J. 1984a Bracket formulation for irreversible classical fields. Physics Letters A 100 (8), 423–427.
  * Morrison (1984b) Morrison, Philip J. 1984b Some observations regarding brackets and dissipation. Tech. Rep.. Center for Pure and Applied Mathematics Report PAM–228, University of California, Berkeley.
  * Morrison (1986) Morrison, Philip J. 1986 A paradigm for joined Hamiltonian and dissipative systems. Physica D: Nonlinear Phenomena 18 (1-3), 410–419.
  * Morrison (1998) Morrison, P. J. 1998 Hamiltonian description of the ideal fluid. Reviews of Modern Physics 70 (2), 467–521.
  * Morrison (2017) Morrison, P. J. 2017 Structure and structure-preserving algorithms for plasma physics. Physics of Plasmas 24 (5), arXiv: 1612.06734.
  * Morrison & Updike (2023) Morrison, Philip J. & Updike, Michael H. 2023 An inclusive curvature-like framework for describing dissipation: metriplectic 4-bracket dynamics, arXiv: 2306.06787.
  * Qin et al. (2016) Qin, Hong, Liu, Jian, Xiao, Jianyuan, Zhang, Ruili, He, Yang, Wang, Yulei, Sun, Yajuan, Burby, Joshua W., Ellison, Leland & Zhou, Yao 2016 Canonical symplectic particle-in-cell method for long-term large-scale simulations of the Vlasov-Maxwell equations. Nuclear Fusion 56 (1).
  * Quispel & Turner (1996) Quispel, G R W & Turner, G S 1996 Discrete gradient methods for solving odes numerically while preserving a first integral. Journal of Physics A: Mathematical and General 29 (13), L341.
  * Schnake et al. (2024) Schnake, Stefan, Kendrick, Coleman, Endeve, Eirik, Stoyanov, Miroslav, Hahn, Steven, Hauck, Cory D, Green, David L, Snyder, Phil & Canik, John 2024 Sparse-grid Discontinuous Galerkin Methods for the Vlasov-Poisson-Lenard-Bernstein Model , arXiv: 2402.06493.
  * Shiroto & Sentoku (2019) Shiroto, Takashi & Sentoku, Yasuhiko 2019 Structure-preserving strategy for conservative simulation of the relativistic nonlinear Landau-Fokker-Planck equation. Physical Review E 99 (5), 1–7.
  * Squire et al. (2012) Squire, J, Qin, H & Tang, W M 2012 Geometric integration of the Vlasov-Maxwell system with a variational particle-in-cell scheme. Physics of Plasmas 19 (8), 84501.
  * Taitano et al. (2015) Taitano, W. T., Chacón, L., Simakov, A. N. & Molvig, K. 2015 A mass, momentum, and energy conserving, fully implicit, scalable algorithm for the multi-dimensional, multi-species Rosenbluth-Fokker-Planck equation. Journal of Computational Physics 297, 357–380.
  * Tyranowski (2021a) Tyranowski, Tomasz M. 2021a Stochastic variational principles for the collisional Vlasov-Maxwell and Vlasov-Poisson equations. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 477 (2252), arXiv: 2102.09611.
  * Tyranowski (2021b) Tyranowski, Tomasz M. 2021b Stochastic variational principles for the collisional Vlasov-Maxwell and Vlasov-Poisson equations. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 477 (2252).
  * Yoon & Chang (2014) Yoon, E. S. & Chang, C. S. 2014 A Fokker-Planck-Landau collision equation solver on two-dimensional velocity grid and its application to particle-in-cell simulation. Physics of Plasmas 21 (3).
  * Zhang & Gamba (2017) Zhang, Chenglong & Gamba, Irene M. 2017 A conservative scheme for Vlasov Poisson Landau modeling collisional plasmas. Journal of Computational Physics 340, 470–497.


# The chirality-dependent second-order spin current
in systems with time-reversal symmetry

Ryosuke Hirakida  Department of Physics, The University of Tokyo, Bunkyo,
Tokyo 113-0033, Japan  Junji Fujimoto  Department of Physics, The University
of Tokyo, Bunkyo, Tokyo 113-0033, Japan  Masao Ogata  Department of Physics,
The University of Tokyo, Bunkyo, Tokyo 113-0033, Japan  Trans-scale Quantum
Science Institute, The University of Tokyo, Bunkyo, Tokyo 113-0033, Japan

###### Abstract

Spin currents proportional to the first- and second-order of the electric
field are calculated in a specific tight-binding model with time-reversal
symmetry. Specifically, a tight-binding model with time-reversal symmetry is
constructed with chiral hopping and spin-orbit coupling. The spin conductivity
of the model is calculated using the Boltzmann equation. As a result, it is
clarified that the first-order spin current of the electric field vanishes,
while the second-order spin current can be finite. Furthermore, the spin
current changes its sign by reversing the chirality of the model. The present
results reveal the existence of spin currents in systems with time-reversal
symmetry depending on the chirality of the system. They may provide useful
information for understanding the chirality-dependent spin polarization
phenomena in systems with time-reversal symmetry.

††preprint: APS/123-QED

##  I Introduction

Three-dimensional structures that do not have mirror symmetry are called
chiral. This symmetry breaking gives rise to various physical properties such
as electrical magnetochiral effect [1] and chiral soliton lattice [2], which
will be useful for spintronics or magnetic devices. Among these properties,
the spin polarization phenomenon called ‘chirality-induced spin selectivity’
(CISS) has recently attracted much attention [3, 4]. CISS was first discovered
in DNA in 2011 [5, 6] and has been intensively studied in organic molecules.
CISS in DNA was reported to produce spin polarization rates of up to
$60\text{\,}\%$ [5] and is expected to be applied to efficient spin current
generation techniques in the future. In 2020, it was shown that the inorganic
crystal CrNb3S6 also shows CISS [7, 8]. More recently, CISS was confirmed also
in inorganic crystals such as NbSi2 and TaSi2 [9].

While the theory of CISS in organic materials has been intensively
investigated [10], that in inorganic materials has not been explored. Because
the CISS is observed in the paramagnetic phase [7, 9], the spin polarization
occurs in the states with time-reversal symmetry. However, it is known that a
spin current linear in the electric field does not arise in single-channel
systems with time-reversal symmetry [11]. This theorem makes the construction
of the CISS theory difficult.

In this study, we construct a chiral tight-binding model with time-reversal
symmetry and evaluate the spin conductivity in the first- and second-order
with respect to the electric field. Specifically, we study a tight-binding
model with chiral hopping and spin-orbit coupling to calculate the spin
conductivity using the Boltzmann transport equation. We show that the first-
order spin conductivity vanishes in the tight-binding model, while the second-
order spin conductivity can be finite. Furthermore, the spin conductivity
changes its sign by reversing the chirality of the model. The present results
reveal the existence of spin currents in systems with time-reversal symmetry
and the existence of chirality-dependent spin currents. They provide useful
information for understanding the chirality-dependent spin polarization
phenomena in systems with time-reversal symmetry.

This paper is organized as follows. In Sect. II, we start with symmetry
analysis of spin conductivities in the first- and second-order with respect to
the electric field in the Hamiltonians with time-reversal symmetry. In Sect.
III, we introduce a specific chiral tight-binding model and calculate the
chirality dependence of the second-order spin conductivity. We also discuss
the parameters required for the chirality-dependent spin conductivity to
arise. Section IV is devoted for discussions and conclusions.

##  II Symmetry Analysis

We first discuss whether the spin current flows when time-reversal symmetry is
imposed on the Hamiltonian in general. We consider a model that is diagonal
with respect to the spin but includes sublattice degrees of freedom, i.e., we
assume that the spin $s_{z}$ ($s_{z}=\uparrow,\downarrow$) along the $z$-axis
is a good quantum number. In this case, the Hamiltonian in the wavenumber
representation can be written as

$\displaystyle\mathcal{H}({\bf\it k})=\begin{pmatrix}\mathcal{H}({\bf\it
k},\uparrow)&O\\\ O&\mathcal{H}({\bf\it k},\downarrow)\\\ \end{pmatrix},$ (1)

where $\mathcal{H}({\bf\it k},\uparrow)$ and $\mathcal{H}({\bf\it
k},\downarrow)$ are matrices whose matrix elements are between the band index
$n$.

By solving the Boltzmann transport equation within the approximation of
constant relaxation time, the spin-resolved first-order conductivity
$\sigma_{zz}(s_{z})$ and spin-resolved second-order conductivity
$\sigma_{zzz}(s_{z})$ with respect to the electric field $E_{z}$, i.e.,
$j_{z}(s_{z})=\sigma_{zz}(s_{z})E_{z}+\sigma_{zzz}(s_{z}){E_{z}}^{2}+\cdots$,
can be written as [12]

$\displaystyle\sigma_{zz}(s_{z})$
$\displaystyle=\frac{e^{2}\tau}{\hbar^{2}}\sum_{n}\int_{\mathrm{BZ}}\frac{\differential^{3}k}{(2\pi)^{3}}\
\frac{\partial^{2}\varepsilon_{n}({\bf\it k},s_{z})}{{\partial
k_{z}}^{2}}f(\varepsilon_{n}({\bf\it k},s_{z})),$ (2)
$\displaystyle\sigma_{zzz}(s_{z})$
$\displaystyle=\frac{e^{3}\tau^{2}}{\hbar^{3}}\sum_{n}\int_{\mathrm{BZ}}\frac{\differential^{3}k}{(2\pi)^{3}}\
\frac{\partial^{3}\varepsilon_{n}({\bf\it k},s_{z})}{{\partial
k_{z}}^{3}}f(\varepsilon_{n}({\bf\it k},s_{z})),$ (3)

where $e(<0)$ is the electric charge, $\tau$ represents the relaxation time,
$\varepsilon_{n}({\bf\it k},s_{z})$ represents the eigenenergy of the band
index $n$ with wave number ${\bf\it k}$ and spin $s_{z}$, and $f(\varepsilon)$
represents the Fermi distribution function.
$\int_{\mathrm{BZ}}\differential^{3}k$ means the integration over the entire
Brillouin zone. The first-order spin conductivity $\sigma^{\mathrm{s}}_{zz}$
and the second-order spin conductivity $\sigma^{\mathrm{s}}_{zzz}$ are defined
as

$\displaystyle\sigma^{\mathrm{s}}_{zz}$
$\displaystyle=\sigma_{zz}(\uparrow)-\sigma_{zz}(\downarrow),$ (4)
$\displaystyle\sigma^{\mathrm{s}}_{zzz}$
$\displaystyle=\sigma_{zzz}(\uparrow)-\sigma_{zzz}(\downarrow).$ (5)

We further suppose that the Hamiltonian in Eq. (1) has the time-reversal
symmetry that can be written as $\Theta=-is_{y}K$, where $K$ is the complex
conjugate operator. In this case, the eigenenergy of the Hamiltonian in Eq.
(1) satisfies

$\displaystyle\varepsilon_{n}({\bf\it k},s_{z})=\varepsilon_{n}(-{\bf\it
k},-s_{z}).$ (6)

Using this equation in Eqs. (2) and (3), we can readily see that

$\displaystyle\sigma_{zz}(s_{z})$ $\displaystyle=\sigma_{zz}(-s_{z}),$ (7)
$\displaystyle\sigma_{zzz}(s_{z})$ $\displaystyle=-\sigma_{zzz}(-s_{z}),$ (8)

which means that the first-order spin conductivity $\sigma^{\mathrm{s}}_{zz}$
vanishes in systems satisfying the above time reversal symmetry, while the
second-order spin conductivity $\sigma^{\mathrm{s}}_{zzz}$ can be finite even
under such conditions. Note that this argument does not hold if the degrees of
freedom of the time-reversal operator except for spin, for example, the
sublattice degrees of freedom are not an identity matrix.

##  III Model Analysis

Next, we discuss the chirality dependence of the second-order spin
conductivity in a specific tight-binding model. We introduce a three-
dimensional tight-binding model with chiral hopping shown in Fig. 1 [13]. This
model consists of two-dimensional honeycomb lattices stacked along the
$z$-axis. The Hamiltonian can be written as

$\displaystyle H$
$\displaystyle=H_{\mathrm{hop}}+H_{\mathrm{chiral}}+H_{\mathrm{SOC}},$ (9)
$\displaystyle H_{\mathrm{hop}}$ $\displaystyle=t_{1}\sum_{\langle
i,j\rangle,l,\alpha}c_{i,l,\alpha}^{\dagger}c_{j,l,\alpha}+\mathrm{h.c.},$
(10) $\displaystyle H_{\mathrm{chiral}}$
$\displaystyle=t_{2}\left(\sum_{i\in\mathrm{A},l,\alpha}\sum_{a=1}^{3}c_{i+\chi{\bf\it
b}_{a},l+1,\alpha}^{\dagger}c_{i,l,\alpha}+\sum_{i\in\mathrm{B},l,\alpha}\sum_{a=1}^{3}c_{i-\chi{\bf\it
b}_{a},l+1,\alpha}^{\dagger}c_{i,l,\alpha}\right)+\mathrm{h.c.},$ (11)
$\displaystyle H_{\mathrm{SOC}}$
$\displaystyle=i\lambda\left(\sum_{i\in\mathrm{A},l,\alpha,\beta}\sum_{a=1}^{3}(c_{i+{\bf\it
b}_{a},l,\alpha}^{\dagger}(s_{z})_{\alpha\beta}c_{i,l,\beta}-c_{i-{\bf\it
b}_{a},l,\alpha}^{\dagger}(s_{z})_{\alpha\beta}c_{i,l,\beta})\right.$
$\displaystyle\left.+\sum_{i\in\mathrm{B},l,\alpha,\beta}\sum_{a=1}^{3}(c_{i-{\bf\it
b}_{a},l,\alpha}^{\dagger}(s_{z})_{\alpha\beta}c_{i,l,\beta}-c_{i+{\bf\it
b}_{a},l,\alpha}^{\dagger}(s_{z})_{\alpha\beta}c_{i,l,\beta})\right)+\mathrm{h.c.},$
(12)

where $c_{i,l,\alpha}$ ($c_{i,l,\alpha}^{\dagger}$) is the annihilation
(creation) operator of an electron with spin $\alpha$ at the site $i$ in the
$xy$-plane of the $l$-th layer. $H_{\mathrm{hop}}$ denotes the nearest-
neighbor hopping in the $xy$-plane with the transfer integral
$t_{1}\in\mathbb{R}$. The $\langle i,j\rangle$ spans the nearest-neighbor
sites of the honeycomb lattices as shown in Fig. 11. $H_{\mathrm{chiral}}$
denotes the chiral hopping between the nearest-neighbor layers with the
transfer integral $t_{2}\in\mathbb{R}$. ${\bf\it a}_{a}\ (a=1,2,3)$ are the
vectors pointing the nearest-neighbor sites in the same plane, and ${\bf\it
b}_{a}\ (a=1,2,3)$ are the vectors pointing the next-nearest-neighbor sites in
the same plane. $\sum_{i\in\mathrm{A}}$ ($\sum_{i\in\mathrm{B}}$) means the
summation over the sites in the A (B) sublattice. $\chi=\pm 1$ is a parameter
which represents the chirality of the model. As we can see from Fig. 11, when
we change the sign of $\chi$, the chirality of the hopping is reversed.
Figures 11 and 1 show the right-handed and left-handed models, respectively.
$H_{\mathrm{SOC}}$ represents the spin-orbit coupling associated with the
next-nearest-neighbor hopping in the layer propotional to $s_{z}$ with the
coefficient $\lambda\in\mathbb{R}$. $(s_{z})_{\alpha\beta}$ is the
$(\alpha,\beta)$ component of the spin operator along the $z$-axis.

Figure 1: Tight-binding model described by the Hamiltonian in Eq. (9) [13]. 1
$xy$-plane of the present model. The vectors ${\bf\it a}_{a}$ and ${\bf\it
b}_{a}$ ($a=1,2,3$) are shown in the figure. 1 Nearest-neighbor hopping
(green) and next nearest-neighbor spin-orbit coupling (yellow) in the $xy$
plane. 11 Chiral hopping between the nearest-negihbor layers: 1 the right-
handed model ($\chi=+1$) and 1 the left-handed model ($\chi=-1$).

The Hamiltonian in the wavenumber representation can be written as

$\displaystyle\mathcal{H}({\bf\it k},s_{z},\chi)$ $\displaystyle=d_{0}({\bf\it
k})I+d_{x}({\bf\it k})\tau_{x}$ $\displaystyle+d_{y}({\bf\it
k})\tau_{y}+d_{z}({\bf\it k},s_{z},\chi)\tau_{z},$ (13) $\displaystyle
d_{0}({\bf\it k})$ $\displaystyle=2t_{2}\sum_{a=1}^{3}\cos({\bf\it
k}\cdot{\bf\it b}_{a})\cos(k_{z}c),$ (14) $\displaystyle d_{x}({\bf\it k})$
$\displaystyle=t_{1}\sum_{a=1}^{3}\cos({\bf\it k}\cdot{\bf\it a}_{a}),$ (15)
$\displaystyle d_{y}({\bf\it k})$
$\displaystyle=t_{1}\sum_{a=1}^{3}\sin({\bf\it k}\cdot{\bf\it a}_{a}),$ (16)
$\displaystyle d_{z}({\bf\it k},s_{z},\chi)$
$\displaystyle=-2\sum_{a=1}^{3}\sin({\bf\it k}\cdot{\bf\it b}_{a})(\chi
t_{2}\sin(k_{z}c)+2s_{z}\lambda),$ (17)

where $I$ is the identity matirx of size 2,
${\bf\it\tau}=(\tau_{x},\tau_{y},\tau_{z})$ are the Pauli matrices
representing the sublattice degrees of freedom. The eigenenergies of the
Hamiltonian in Eq. (13) are expressed as

$\displaystyle\varepsilon_{\pm}({\bf\it k},s_{z},\chi)$ $\displaystyle=$
$\displaystyle d_{0}({\bf\it k})\pm\sqrt{{d_{x}({\bf\it
k})}^{2}+{d_{y}({\bf\it k})}^{2}+{d_{z}({\bf\it k},s_{z},\chi)}^{2}}.$ (18)

The band structure of this model is shown in Fig. 2. Apparently,
$\varepsilon_{\pm}({\bf\it k},s_{z},\chi)$ satisfies the time reversal
condition in Eq. (6), and thus the relation in Eq. (7) holds, i.e., the first-
order spin conductivity $\sigma^{\mathrm{s}}_{zz}$ vanishes.

Figure 2: Band structure of the present model. The parameters are
$t_{2}=0.2t_{1}$ and $\lambda=0.06t_{1}$.

To study the chirality dependence of the second-order spin conductivity
$\sigma^{\mathrm{s}}_{zzz}(\chi)$, we study the dependence of
$\varepsilon_{\pm}({\bf\it k},s_{z},\chi)$ on $\chi$. From Eq. (17), we obtain

$\displaystyle d_{z}({\bf\it k},s_{z},\chi)$ $\displaystyle=-d_{z}({\bf\it
k},-s_{z},-\chi).$ (19)

Then we can readily see that

$\displaystyle\varepsilon_{\pm}({\bf\it k},s_{z},\chi)$
$\displaystyle=\varepsilon_{\pm}({\bf\it k},-s_{z},-\chi),$ (20)
$\displaystyle\frac{\partial^{3}\varepsilon_{n}({\bf\it
k},s_{z},\chi)}{{\partial k_{z}}^{3}}$
$\displaystyle=-\frac{\partial^{3}\varepsilon_{n}({\bf\it
k},-s_{z},-\chi)}{{\partial k_{z}}^{3}}.$ (21)

Considering these equations with Eqs. (3) and (18), we obtain

$\displaystyle\sigma_{zzz}(s_{z},\chi)$
$\displaystyle=\sigma_{zzz}(-s_{z},-\chi).$ (22)

Therefore, we can see

$\displaystyle\sigma_{zzz}^{\mathrm{s}}(\chi)$
$\displaystyle=-\sigma_{zzz}^{\mathrm{s}}(-\chi),$ (23)

which means that the spin conductivity of the right-handed model
$\sigma_{zzz}^{\mathrm{s}}(+1)$ and that of the left-handed model
$\sigma_{zzz}^{\mathrm{s}}(-1)$ are of opposite signs in the present chiral
tight-binding model.

In the following, we show that the chirality and the spin-orbit coupling of
the present model play crucial roles in the chirality-dependent second-order
spin current. Here we fix the chirality of the model $\chi$. When we assume
that $t_{2}=0$, i.e., the model has no chirality, we obtain

$\displaystyle\sigma_{zzz}(s_{z},\chi)=\sigma_{zzz}(-s_{z},\chi),$ (24)

because $d_{z}({\bf\it k},s_{z},\chi)=-d_{z}({\bf\it k},-s_{z},\chi)$ holds as
seen from Eq. (17). This means that the second-order spin conductivity
$\sigma^{\mathrm{s}}_{zzz}(\chi)$ vanishes under $t_{2}=0$. Next, when we
assume that $\lambda=0$, i.e., the model has no spin-orbit coupling, we obtain

$\displaystyle\sigma_{zzz}(s_{z},\chi)=\sigma_{zzz}(-s_{z},\chi),$ (25)

because $d_{z}({\bf\it k},s_{z},\chi)=d_{z}({\bf\it k},-s_{z},\chi)$ holds
from Eq. (17). This means that the second-order spin conductivity
$\sigma^{\mathrm{s}}_{zzz}(\chi)$ again vanishes under $\lambda=0$. From the
above discussion, we can conclude that the chirality and the spin-orbit
coupling are necessary for the chirality-dependent second-order spin current
in the present model.

Figure 3:  Calculated spin-resolved first-order conductivity. The parameters
are $t_{2}=0.2t_{1}$ and $\lambda=0.06t_{1}$. The temperature $T$ is set at
zero. 3 Spin-$\uparrow$ conductivity in the right-handed model. 3
Spin-$\downarrow$ conductivity in the right-handed model.

Figure 4:  Calculated second-order spin conductivity. The parameters and
temperature are the same as in Fig. 3. 4 Results for the right-handed model,
and 4 for the left-handed model.

The calculated spin-resolved first-order conductivities and second-order spin
conductivities are shown in Figs. 3 and 4. Figures 33 and 3 depict the
chemical potential dependences of the spin-resolved first-order
conductivities. The spin-$\uparrow$ conductivity has the same dependence as
that of the spin-$\downarrow$ conductivity, which results in the zero spin
conductivity in the first order with respect to the electric field, as
discussed in Sect. II. Figures 44 and 4 depict the chemical potential
dependences of the second-order spin conductivities. The spin conductivity in
the left-handed model has the opposite sign as compared with that in the
right-handed model, indicating that the sign of the spin conductivity can be
controlled by reversing the chirality, as discussed above. Figure 4 shows that
the second-order spin conductivity rapidly changes its sign depending on the
chemical potential in the region of $-2<\mu/t_{1}<2$. This will originate from
the difference in the energy dispersions of the spin-$\uparrow$ and
spin-$\downarrow$ in the region of $-2<\varepsilon/t_{1}<2$ as shown in Fig.
2.

##  IV Conclusion

We explored the existence of spin current and its chirality dependence in
systems with time-reversal symmetry. We evaluated the spin conductivity in the
system with time-reversal symmetry based on the Boltzmann transport equation.
We found that the first-order spin conductivity with respect to the electric
field vanishes, while the second-order spin conductivity can be finite. We
also evaluated the second-order spin conductivity using a specific chiral
tight-binding model. As a result, we revealed the existence of the second-
order spin conductivity depending on the chirality of the system. We also
found that the chirality and the spin-orbit coupling play crucial roles in the
chirality-dependent spin current.

We considered in this study the second-order spin conductivity with respect to
the electric field because the first-order spin conductivity vanishes for the
non-interacting single-channel systems with time-reversal symmetry [11].
However, spin polarization phenomena in the systems with spin hybridization
have not been considered yet. The role of weak antisymmetric spin-orbit
coupling in spin polarization phenomena is also discussed [7]. The mechanism
in which the spin-flipping components of spin-orbit couplings affect the spin-
polarized current remains to be investigated in future studies.

##  V Acknowledgement

This work was supported by Grants-in-Aid for Scientific Research from the
Japan Society for the Promotion of Science (Grant No. JP18H01162), and by JST-
Mirai Program (Grant No. JPMJMI19A1). R. H. was supported by Iwadare
Scholarship Foundation and the Program for Leading Graduate Schools (MERIT-
WINGS).

## References

  * Rikken  _et al._ [2001] G. L. Rikken, J. Fölling, and P. Wyder, Phys. Rev. Lett. 87, 236602 (2001).
  * Togawa  _et al._ [2016] Y. Togawa, Y. Kousaka, K. Inoue, and J. I. Kishine, J. Phys. Soc. Jpn. 85, 1 (2016).
  * Naaman and Waldeck [2012] R. Naaman and D. H. Waldeck, J. Phys. Chem. Lett. 3, 2178 (2012).
  * Naaman and Waldeck [2015] R. Naaman and D. H. Waldeck, Annu. Rev. Phys. Chem. 66, 263 (2015).
  * Göhler  _et al._ [2011] B. Göhler, V. Hamelbeck, T. Z. Markus, M. Kettner, G. F. Hanne, Z. Vager, R. Naaman, and H. Zacharias, Science 331, 894 (2011).
  * Xie  _et al._ [2011] Z. Xie, T. Z. Markus, S. R. Cohen, Z. Vager, R. Gutierrez, and R. Naaman, Nano Lett. 11, 4652 (2011).
  * Inui  _et al._ [2020] A. Inui, R. Aoki, Y. Nishiue, K. Shiota, Y. Kousaka, H. Shishido, D. Hirobe, M. Suda, J. I. Ohe, J. I. Kishine, H. M. Yamamoto, and Y. Togawa, Phys. Rev. Lett. 124, 166602 (2020).
  * Nabei  _et al._ [2020] Y. Nabei, D. Hirobe, Y. Shimamoto, K. Shiota, A. Inui, Y. Kousaka, Y. Togawa, and H. M. Yamamoto, Appl. Phys. Lett. 117, 052408 (2020).
  * Shishido  _et al._ [2021] H. Shishido, R. Sakai, Y. Hosaka, and Y. Togawa, Appl. Phys. Lett. 119, 182403 (2021).
  * Evers  _et al._ [2022] F. Evers, A. Aharony, N. Bar-Gill, O. Entin-Wohlman, P. Hedegård, O. Hod, P. Jelinek, G. Kamieniarz, M. Lemeshko, K. Michaeli, V. Mujica, R. Naaman, Y. Paltiel, S. Refaely-Abramson, O. Tal, J. Thijssen, M. Thoss, J. M. van Ruitenbeek, L. Venkataraman, D. H. Waldeck, B. Yan, and L. Kronik, Adv. Mater. 2106629, 1 (2022).
  * Bardarson [2008] J. H. Bardarson, J. Phys. A: Math. Theor. 41, 405203 (2008).
  * Sodemann and Fu [2015] I. Sodemann and L. Fu, Phys. Rev. Lett. 115, 216806 (2015).
  * Yoda  _et al._ [2015] T. Yoda, T. Yokoyama, and S. Murakami, Sci. Rep. 5, 12024 (2015).