# On discrete gradient vector fields and Laplacians of simplicial complexes

Ivan Contreras  and  Andrew R. Tawfeek<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

Discrete Morse theory, a cell complex-analog to smooth Morse theory allowing
homotopic tools in the discrete realm, has been developed over the past few
decades since its original formulation by Robin Forman in 1998. In particular,
discrete gradient vector fields on simplicial complexes capture important
topological features of the structure. We prove that the characteristic
polynomials of the Laplacian matrices of a simplicial complex are generating
functions for discrete gradient vector fields if the complex is a
triangulation of an orientable manifold. Furthermore, we provide a full
characterization of the correspondence between rooted forests in higher
dimensions and discrete gradient vector fields.

##  1\. Introduction

Since the appearance of discrete Morse theory due to Forman in [15], it has
brought a wide range of applications of this combinatorial model of Morse
theory to mathematics and the sciences alike. Some of the recent applications
include cellular sheaf cohomology [11], topology of random points [19], and
studying human speech from statistical data [8].

Here, we further demonstrate the importance of discrete Morse theory by
showing that discrete gradient fields have a novel intimate connection with
the topological structure of a simplicial complex as characterized by its
Laplacian $\Delta_{d}:=\partial_{d}\partial_{d}^{T}$. The discrete Laplacian
plays an important role in various subjects such as chip-firing on graphs [5],
critical groups in higher dimensions [12], and has even been generalized to
matroids [21].

In particular, we prove (Theorem 4.5) that when $K$ is a triangulation of an
orientable $d$-manifold, we have that

$det(\Delta_{d}+\lambda I)=\sum_{i=0}^{n}|\mathcal{M}_{m-i}(K)|\lambda^{n-i},$

where $m=|K_{d}|$ and $n=|K_{d-1}|$. Here, $\mathcal{M}_{\ell}(K)$ denote the
set of gradients on $K$ having $\ell$-many critical cells in dimension $d$ and
all cells critical in dimension $n$ for $n<d-1$. In the graph-theoretic case
the above statement holds for simple graphs (Theorem 3.5).

This paper is organized as follows. In Section 2 we review the background and
notation concerning simplicial complexes, higher dimensional rooted forests,
and discrete Morse theory. In Section 3 we prove a special case of our main
statement exclusively in the $1$-dimensional case, providing a graph theoretic
interpretation of the main result, and motivating our methods in the general
case. In Section 4, we conclude by proving the general statement, and then we
present some future directions in Section 5.

### Keywords

discrete Laplacian, simplicial complexes, discrete Morse theory, discrete
gradient vector fields, matchings, rooted forests, spectral graph theory.

## Acknowledgements

We particularly thank Alejandro Morales for his valuable suggestions, comments
and guidance during the project. We would also extend our thanks to David A.
Cox, Michael Ching, Nathan Pflueger and Nicholas A. Scoville, for very helpful
discussions and suggestions. The second author was partially supported by the
Gregory S. Call Research Fellowship.

##  2\. Background and Notation

###  2.1. Simplicial complexes

A simplicial complex111All simplicial complexes considered are finite-
dimensional and connected, unless stated otherwise. on $n$ vertices is a
collection $K$ of subsets of $[n]:=\\{1,2,\dots,n\\}$ which is closed under
taking subsets, i.e. if $f\in K$ and $g\subseteq f$, then $g\in K$. An element
$f\in K$ is called a $d$-dimensional face or cell (or $d$-face and $d$-cell)
if $|f|=d+1$. We use $K_{d}$ to denote the subset of $K$ consisting of only
$d$-cells and $K_{\leq d}$ to denote the collection of all cells of dimension
$\leq d$, called the $d$-skeleton. When $K$ consists of all subsets of $[n]$,
we call it a $n$-simplex, and denote it by $\Delta^{n}$. If $g\subsetneq f$ in
$K$, we say $g$ is a face of $f$ and write $g<f$ (with the notation chosen due
to the later discussed poset structure on cells). Cells $f\in K$ that are
contained in no other cell of $K$ are said to be maximal, and we say $K$ is an
$n$-dimensional simplicial complex if its maximal cells are at most
$n$-dimensional. For instance, a graph is a $1$-dimensional simplicial
complex.

Each $n$-dimensional simplicial complex has a geometric realization as a
subspace of $\mathbb{R}^{n+1}$, in which each $d$-face is represented by the
convex hull of $d+1$ affinely independent points.

Figure 1. An example of a $2$-dimensional simplicial complex $K$, with the
standard lexicographic orientations indicated, and the the corresponding
incidence matrix $\partial^{K}_{2}$.

To each simplicial complex $K$ we associate a chain complex $C_{\bullet}(K)$
of abelian groups222We will always consider $\mathbb{Z}$-modules.

${\cdots}$${{C_{d+1}(K)}}$${{C_{d}(K)}}$${{C_{d-1}(K)}}$${\cdots}$$\scriptstyle{\partial^{K}_{d}}$$\scriptstyle{\partial^{K}_{d+1}}$

where we let $C_{d}(K)$, for $d\geq 0$, denote the free abelian group
generated by all orientations of all simplices of $K$ of dimension $d$, under
the quotient by the relation that the sum of two opposite orientations is
zero. The map $\partial^{K}_{d}$ (or just $\partial_{d}$ if no ambiguity will
arise) is called the $d$th-boundary map, which we routinely identify with its
corresponding matrix, called the $d$th-incidence matrix of the simplicial
complex $K$, defined as follows:

  * •

the rows of $\partial_{d}$ are indexed by the $(d-1)$-faces and the columns
are indexed by the $d$-faces,

  * •

the entry $\partial_{d}[r,f]$ corresponding to a $(d-1)$-face $r$ and a
$d$-face $f$ is $(-1)^{j}$ if $f=\\{v_{0},v_{1},\dots,v_{d}\\}$ with
$v_{0}<v_{1}<\cdots<v_{d}$ and $r=f\setminus\\{v_{j}\\}$ for some $j$, and $0$
otherwise.

It is a standard exercise to show that $\partial_{d-1}\circ\partial_{d}=0$.
Elements of $C_{n}(K)$ are said to be $n$-chains. A $n$-cycle of $K$ is a non-
zero element of $C_{n}(K)$ that is in $\ker\partial_{n}$. Since for any $n\geq
1$ the image of an element $\tau\in C_{n}(K)$ is necessarily sent to the
kernel of $\partial_{n-1}$, we say a $(n-1)$-cycle is a $K$-boundary if it is
inside $\text{im }\partial_{n}$.

We define $H_{d}(K)$ to be the $d$th integral homology group of the simplicial
complex $K$ (or equivalently, of the chain complex $C_{\bullet}(K)$), given by

$H_{d}(K)=\frac{\ker\partial_{d}}{\text{im }\partial_{d+1}}.$

This group is in-fact finitely generated, so we have by the fundamental
theorem of finitely generated abelian groups that there exists primes
$p_{1},\dots,p_{s}$ and positive integers $n_{1},\dots,n_{s}$ such that

$H_{d}(K)\cong\mathbb{Z}^{r}\times\mathbb{Z}/{p_{1}^{n_{1}}}\mathbb{Z}\times\cdots\times\mathbb{Z}/{p_{s}^{n_{s}}}\mathbb{Z},$

and that this decomposition is uniquely determined by the group $H_{d}(K)$ (up
to permutation of the factors). The value $r$ is called the rank or Betti
number, denoted by $\beta_{d}(K)$ (or just $\beta_{d}$, if no ambiguity would
arise). The product of the corresponding finite groups is called the torsion
of $H_{d}(K)$. If $H_{d}(K)$ has no finite subgroups, we say $K$ is torsion-
free in dimension $d$, or that $H_{d}(K)$ is torsion-free.

Similarly, we have that the $d$th homology group of $K$ relative to a
subcomplex $A$ is as follows: given a subcomplex $A\subseteq K$, as the
boundary maps $\partial^{K}_{n}:C_{n}(K)\to C_{n-1}(K)$ of $C_{\bullet}(K)$
leave $C_{\bullet}(A)$ invariant333That is to say,
$\partial_{n}^{K}|_{C_{n}(A)}$ maps $C_{n}(A)$ to $C_{n-1}(A)$, and hence is
precisely (after a restriction of the codomain of
$\partial_{n}^{K}|_{C_{n}(A)}$ to $C_{n-1}(A)\subseteq C_{n-1}(K)$)
$\partial_{n}^{A}$ of $C_{\bullet}(A)$., this induces the a short exact
sequence

$0\rightarrow C_{\bullet}(A)\rightarrow C_{\bullet}(K)\rightarrow
C_{\bullet}(K)/C_{\bullet}(A)\rightarrow 0$

between the chain complexes, where the boundary maps of
$C_{\bullet}(K)/C_{\bullet}(A)$,

$\partial^{*}_{n}:C_{n}(K)/C_{n}(A)\to C_{n-1}(K)/C_{n-1}(A),$

are given by
$\partial^{*}_{n}(\sigma+C_{n}(A))=\partial_{n}(\sigma)+C_{n-1}(A)$. We
classically denote $C_{\bullet}(K)/C_{\bullet}(A)$ by $C_{n}(K,A)$, and write
$H_{d}(K,A)$ to denote the $d$th homology group of $K$ relative to $A$,
defined444For the readers familiar with CW complexes, a perhaps more
geometrically intuitive way to visualize this group is to consider $K$ to be a
CW complex, in which case due to the equivalence of cellular and simplicial
homology (see [16]), we have that $H_{d}(K,A)\cong H_{d}(K/A)$ for $d\geq 1$
and $H_{0}(K,A)\cong H_{0}(K/A)\oplus\mathbb{Z}$. to be the $d$th homology
group of the chain complex $C_{\bullet}(K,A)$.

Let us now introduce the object of our main interest: the $d$th Laplacian of a
simplicial complex $K$ of dimension $d$, defined to be the
$|K_{d-1}|\times|K_{d-1}|$ -matrix

$\Delta^{K}_{d}=\partial_{d}\cdot\partial_{d}^{T}.$

The rows and columns of $\Delta^{K}_{d}$ (or just $\Delta_{d}$, if no
ambiguity would arise) are indexed by $(d-1)$-faces. The $(r,s)$th entry is
given as follows.

$\Delta_{d}^{K}[r,s]=\begin{cases}\left|\coprod_{\sigma\in
K_{d}}\sigma\cap\\{r\\}\right|&\text{if $r=s$}\\\ (-1)^{i+j}&\text{if $r\neq
s$ and $\exists f\in K_{d}$ such that $r=f\setminus\\{v_{i}\\}$ and
$s=f\setminus\\{v_{j}\\}$}\\\ 0&\text{otherwise ($r,s$ not incident to a
common $d$-face).}\end{cases}$

In the case that $K$ is a $1$-dimensional simplicial complex, i.e. a graph,
one often writes $\Delta_{1}$ as just $\Delta$ and explicitly calls it the
graph Laplacian.

###  2.2. Rooted forests

We now overview some of the basic notions on rooted forests, as introduced in
[6], which will be important for obtaining our higher dimensional results.

######  Definition 2.1.

Let $K$ be a $d$-dimensional complex, and let $F$ be a subset of the
$d$-faces. We call $F$ a forest if the corresponding columns of
$\partial_{d}^{K}$ are independent.

Equivalently, note that a collection of $d$-face $F$ of a $d$-dimensional
simplicial complex $K$ is a forest if and only if $\partial^{F}_{d}$ is
injective.

######  Remark 1.

As pointed out in [6], there is a connection between forests and matroids. For
a $d$-dimensional simplicial complex $K$, the linear matroid of
$\partial^{K}_{d}$, called the simplicial matroid $M_{K}$ of $K$, is the
matroid having the ground set of $d$-faces and independent sets corresponding
to forests. The bases of $M_{K}$ are the spanning forests of $K$: the
$d$-faces corresponding to columns of $\partial^{K}_{d}$ that form a basis of
the columns.

######  Definition 2.2.

Let $K$ be a $d$-dimensional complex. Let $R$ be a subset of the
$(d-1)$-faces, and let $\overline{R}=K_{d-1}\setminus R$. We say that $R$ is
relatively-free if the rows of $\partial_{d}^{K}$ corresponding to the faces
in $\overline{R}$ are of maximal rank. We say that $R$ is relatively-
generating if the rows of $\partial^{K}_{d}$ corresponding to the faces of
$\overline{R}$ are independent. We say $R$ is a root of $K$ if it is both
relatively-free and relatively-generating, that is, the rows of
$\partial_{d}^{K}$ corresponding to the faces in $\overline{R}$ form a basis
of the rows.

A more geometric view is as follows: given a subset of the $(d-1)$-faces of a
$d$-dimensional complex $K$, we have that $R$ is relatively-free if and only
if $R$ contains no $K$-boundary (i.e. $C_{d}(R)$ contains no $K$-boundary). In
a similar vein, $R$ is relatively-generating if and only if for all cells
$\sigma\in\overline{R}$, we have that $R\cup\\{\sigma\\}$ contains a
$K$-boundary.

The conditions of being relatively-free and relatively-generating are natural
generalizations of the properties of roots in the graph-theoretical case. A
root $R$ being relatively-free generalizes the condition of containing at most
one vertex per connected component of a forest, and $R$ being relatively-
generating generalizes the condition of having at least one vertex per
connected component in a forest.

######  Definition 2.3.

Let $K$ be a $d$-dimensional simplicial complex. A rooted forest of $K$ is a
pair $(F,R)$ where $F$ is a $d$-dimensional forest of $K$, and $R$ is a root
of $F$.

Lastly, given a rooted forest $(F,R)$, a fitting orientation of $(F,R)$ is a
bijection $\psi$ between $\overline{R}:=K_{d-1}\setminus R$ and $F$, such that
each face $r\in\overline{R}$ is mapped to a face $f\in F$ containing $r$. This
generalizes the (unique) fitting orientation on a rooted forest of a graph,
where each edge is oriented towards the root of its connected component.

Figure 2. Examples of rooted forests. The red $1$-cells in figures (a) – (c)
and red vertex in (d) denote the roots of the corresponding forest. The arrows
represent fitting orientations.

A classic result of algebraic graph theory is that the Laplacian and rooted
forests are intimately connected. Let $\Gamma=(V_{\Gamma},E_{\Gamma})$ be a
graph on $n$ vertices having $m$ edges. Recall that an edge-subgraph of
$\Gamma$ is constructed by taking a subset $S\subseteq E_{\Gamma}$ together
with all the vertices of $\Gamma$ incident with some edge belonging to $S$.

######  Theorem 2.4 (Biggs 1974).

The coefficients $q_{i}$ of the characteristic polynomial $\det(\Delta-\lambda
I)=\sum_{i=0}^{n}(-1)^{i}q_{i}\lambda^{n-i}$, where $\Delta$ is the graph
Laplacian, are given by the formula

$q_{i}=\sum_{F\in\Omega_{i}}|V_{F}|\qquad(1\leq i\leq n),$

where $\Omega_{i}$ consists of all the edge-subgraphs of $\Gamma$ which have
$i$ edges and are forests.

The above theorem also implies Kirchoff’s Theorem, i.e. the coefficient of
$\lambda^{1}$ in $\det(\Delta-\lambda I)$ is subsequently

$\lambda_{1}\cdots\lambda_{n}=|V_{\Gamma}|\ \tau(\Gamma),$

where the $\lambda_{i}$ above are the non-zero eigenvalues of the graph
Laplacian and $\tau(\Gamma)$ is the number of spanning trees of $\Gamma$. The
following theorem from [6] generalizes Theorem 2.4 to higher dimensions.

######  Theorem 2.5 (Bernardi-Klivans 2016).

For a $d$-dimensional simplicial complex $K$, the characteristic polynomial of
the Laplacian matrix $\Delta_{d}$ gives a generating function for the rooted
forests of $K$. More precisely,

$\sum_{\text{$(F,R)$ rooted forest of
$K$}}|H_{d-1}(F,R)|^{2}\lambda^{|R|}=\det(\Delta_{d}+\lambda I),$

where $I$ is the identity matrix of dimension $|K_{d-1}|$.

###  2.3. Discrete Morse theory

Lastly, we overview the foundations of discrete Morse theory – with a
particular emphasis placed on concepts relevant to discrete gradient vector
fields.

######  Definition 2.6.

Let $K$ be a simplicial complex. A discrete Morse function on $K$ is a
function $f:K\to\mathbb{R}$ satisfying the following two inequalities for all
$\sigma\in K_{p}$:

$|\\{\tau\in K_{p-1}\ |\ \tau<\sigma\text{ and }f(\sigma)\leq f(\tau)\\}|\leq
1;$ $|\\{\tau\in K_{p+1}\ |\ \sigma<\tau\text{ and }f(\tau)\leq
f(\sigma)\\}|\leq 1.$

That is, for any given cell, there is at most one coface with a smaller or
equal value and at most one face with a larger or equal value. It is possible
that both of the above sets have zero cardinality for some $\sigma\in K_{p}$,
in which case we say the cell $\sigma$ is a critical cell of index $p$. In
other words, we say $f$ is critical at a cell $\sigma$ if and only if, locally
at $\sigma$, the function $f$ is increasing with the dimension of the cells.
We denote by $c_{i}(f)$ the number of critical cells of index $i$ for a
discrete Morse function $f$.

######  Theorem 2.7 (Forman 2002).

Let $f$ be a discrete Morse function on a $n$-dimensional simplicial complex
$K$. We then have that555Yes, the second statement is indeed called an
“inequality” in the literature.

  1. (1)

First Weak Morse Inequality: $\beta_{i}(K)\leq c_{i}(f)$ for all $i$, and

  2. (2)

Second Weak Morse Inequality:

$\sum_{i=0}^{n}(-1)^{i}c_{i}(f)=\chi(K).$

where $\chi(K)$ denotes the Euler characteristic of $K$, i.e.
$\sum_{i=0}^{n}(-1)^{i}|K_{i}|$.

Note that there always exists a discrete Morse function on a complex, e.g. the
trivial function $f:K\to\mathbb{R}$ defined by

$f(\sigma)=\dim\sigma,$

where $\dim\sigma$ is the dimension of the cell $\sigma$, is a discrete Morse
function on $X$. Under this function, every cell is critical. The constant
function is not a Morse function unless one is merely considering disjoint
unions of $0$-cells.

Historically, discrete Morse functions were studied on regular CW complexes,
though today research focuses on the case of simplicial complexes. This is
largely due to their equivalence: a regular CW complex is isomorphic to a
simplicial complex in a natural definition of the category of combinatorial
cell complexes. For an explicit discussion of this, we direct the reader to
[1].

We can now build towards discrete gradient vector fields. The following lemma
will make sure the definition we provide will be well-defined.

######  Lemma 2.8.

Suppose $\sigma$ is a $p$-cell that is not critical for a discrete Morse
function $f$. Then exactly one of the following conditions holds:

  1. (1)

There is a $(p+1)$-cell $\tau>\sigma$ with $f(\tau)\leq f(\sigma)$;

  2. (2)

There is a $(p-1)$-cell $\alpha<\sigma$ with $f(\alpha)\geq f(\sigma)$.

###### Proof.

Assume for contradiction both conditions hold for $\sigma$. We then have that
$\tau$ has another $(p-1)$-cell, say $\sigma^{\prime}$, where
$\alpha<\sigma^{\prime}<\tau$. Since $\alpha<\sigma$ and $f(\alpha)\geq
f(\sigma)$, we have that $f(\alpha)<f(\sigma^{\prime})$. Similarly,
$f(\sigma^{\prime})<f(\tau)$. Hence we get that

$f(\tau)\leq f(\sigma)\leq f(\alpha)<f(\sigma^{\prime})<f(\tau),$

which is a contradiction. ∎

Figure 3. An example of a discrete Morse function on a triangulation of a
torus.

######  Definition 2.9.

A discrete vector field $V$ on a simplicial complex $K$ is a collection of
pairs of faces where each pair $(\sigma,\tau)\in V$ satisfies:

  1. (1)

$\dim\tau=\dim\sigma+1$,

  2. (2)

$\sigma<\tau$, and

  3. (3)

every face of the complex is in at most one pair.

If $f$ is a discrete Morse function on $K$, it induces the discrete vector
field

$V_{f}:=\\{(\sigma,\tau)\in K\times K\ |\ \dim\tau=\dim\sigma+1\text{ and
}f(\tau)\leq f(\sigma)\\}.$

We say that a discrete vector field $V$ is a discrete gradient vector field if
it is induced by a discrete Morse function on $K$. If $(\sigma,\tau)\in V$, we
call $(\sigma,\tau)$ an arrow, matching, or vector, with $\tau$ the head and
$\sigma$ the tail.

Observe that fitting orientations are a special case of discrete vector
fields, where pairings occur only in the two top-most dimensions. In general,
fitting orientations do not coincide with discrete gradient vector fields, as
the latter is in-fact characterized by not having closed paths of faces.

Figure 4. Examples of fitting orientations having closed cycles on rooted
forests of triangulations of the real projective plane and the Möbius strip.

We conclude with presenting an important equivalent way of viewing discrete
gradient vector fields. Recall that given a simplicial complex $K$, we may
associate to it a Hasse diagram666We sometimes abuse notation and treat the
Hasse diagram $\mathcal{H}_{K}$ alternately as both a poset and a graph. It
will be clear from context., which we denoted by $\mathcal{H}_{K}$, defined as
follows: if $\sigma\in K$, we write $\sigma\in\mathcal{H}_{K}$ for the
corresponding node. There is an edge between two vertices
$\sigma,\tau\in\mathcal{H}_{K}$ if and only if $\tau$ is a face of co-
dimension-$1$ of $\sigma$. It is standard to draw Hasse diagrams with levels
increasing from the base (see Figure $5$).

Recall that a matching on a graph $\Gamma$ is a subset $M\subseteq E_{\Gamma}$
of the edges such that no two edges in $M$ have common vertices. Importantly
observe that if $V$ a discrete vector field on $K$, then $M_{V}:=V$ is a(n
edge) matching on $\mathcal{H}_{K}$. If $V_{f}$ is a discrete gradient vector
field, then $M_{f}:=M_{V_{f}}$ is called an acyclic matching on
$\mathcal{H}_{K}$. In this language, fitting orientations on a rooted forest
$(F,R)$ may be seen as matchings on the subposet $\mathcal{H}_{(F,R)}$ of
$\mathcal{H}_{K}$ induced by $F\cup\overline{R}$.

######  Remark 2.

Acyclic matchings on $\mathcal{H}_{K}$ are also sometimes called Morse
matchings and the computational complexiy of finding them is something of a
great deal of interest, as it allows us to find a smaller complex which is
homotopy equivalent to the original complex. For more on this, [18] is a great
place to start.

A matching $M_{V}$ on $\mathcal{H}_{K}$ from a discrete vector field $V$ are
commonly illustrated by orienting edges on $\mathcal{H}_{K}$: if $e\in M_{V}$,
it is oriented upwards, otherwise it is oriented downwards. The matching is
acyclic if and only if the corresponding oriented Hasse diagram has no
directed cycles777In fact, more is known: if an oriented Hasse diagram arising
from a matching contains a cycle, it is necessarily contained within two
levels. See §2.2 of [26] for details..

Figure 5. An example of a discrete gradient vector field $V$ on the
$2$-simplex $K=\Delta^{2}$ and the associated orientation on the Hasse diagram
$\mathcal{H}_{K}$.

We conclude by mentioning the specific sets we aim to enumerate.

######  Definition 2.10.

For a simplicial complex $K$, let $\mathcal{M}(K)$ denote the set of
equivalence classes of discrete Morse functions on $K$, where $f\sim g$ if and
only if $V_{f}=V_{g}$. Furthermore, for any $\ell\in\mathbb{Z}_{\geq 0}$,
define the subset $\mathcal{M}_{\ell}(K)\subseteq\mathcal{M}(K)$ by

$\mathcal{M}_{\ell}(K):=\\{f\in\mathcal{M}(K)\ |\ c_{d}(f)=\ell\text{ and
}c_{n}(f)=|K_{n}|\text{ for }n<d-1\\}.$

Observe it follows from the Second Weak Morse Inequality that as we have
$c_{d}(f)=\ell$ and $c_{n}(f)=|K_{n}|\text{ for }n<d-1$, then from canceling
terms, we obtain

$c_{d}(f)-c_{d-1}(f)=|K_{d}|-|K_{d-1}|,$

so as $c_{d}(f)=\ell$ for any $f\in\mathcal{M}_{\ell}(K)$, it necessarily
follows that $c_{d-1}(f)=\ell-|K_{d}|+|K_{d-1}|$ for all discrete Morse
functions $f\in\mathcal{M}_{\ell}(K)$.

An equivalent formulation of this is having $\mathcal{M}(K)$ denote the set of
acyclic matchings on $\mathcal{H}_{K}$ and $\mathcal{M}_{\ell}(K)$ denote the
set of acyclic matchings on the subposets $\mathcal{H}_{(F,R)}$ such that

$|F|=|K_{d}|-\ell\qquad\text{and}\qquad|\overline{R}|=\ell-|K_{d}|+|K_{d-1}|.$

##  3\. Gradients of $1$-dimensional Simplicial Complexes

In this section, we will discuss our results towards enumerating all the
discrete gradient vector fields that may be imposed on a $1$-dimensional
complex, i.e. a graph. We approach this by finding a generating polynomial,
and furthermore, we show that this generating polynomial is equal to the
characteristic polynomial of the graph Laplacian.

Throughout, assume $\Gamma$ is a finite, simple, and not neccesarily connected
graph and $f$ is a discrete Morse function on $\Gamma$. We commonly denote a
subgraph $\Gamma^{\prime}$ of $\Gamma$ with the notation
$\Gamma^{\prime}\subseteq\Gamma$, and write $f|_{\Gamma^{\prime}}$ for the
discrete Morse function induced on the subgraph by $f$. All graphs are
considered labeled, so their automorphism group is always trivial.

We will identify the discrete gradient vector field $V_{f}$ with a directed
graph $\Gamma_{f}$ in the following natural way: the underlying undirected
graph of $\Gamma_{f}$ is $\Gamma$ and an edge $(v,w)$ is directed in
$\Gamma_{f}$ if and only if $(\\{v\\},\\{v,w\\})\in V_{f}$, otherwise it is
considered undirected. This provides a helpful characterization of critical
cells in the $1$-dimensional case:

  * •

an edge is a critical cell of $f$ on $\Gamma$ if and only if it is undirected
in $\Gamma_{f}$,

  * •

a vertex is a critical cell of $f$ on $\Gamma$ if and only if it is a sink in
$\Gamma_{f}$.

###  3.1. Characterizing gradients

We now can begin descending towards the heart of the enumeration problem. The
counting of gradients on a graph becomes clear once we begin removing critical
edges, which is precisely what the next definition gives us the ability to
focus more closely on.

######  Definition 3.1.

Let $f$ be a discrete Morse function on the graph $\Gamma$ and $X\subseteq
E_{\Gamma}$ be the set of critical edges under $f$. Then the graph
$R(\Gamma,f)$, which we call the remainder, is the edge-subgraph of
$E_{\Gamma_{f}}\setminus X$ in the directed graph $\Gamma_{f}$.

We can now construct a natural bijection between discrete gradient vector
fields on graphs and rooted forests. Note that this was obtained through
different means in [9], relying instead on the acyclicity property of
gradients.

######  Lemma 3.2.

For a graph $\Gamma$, we have a bijection

$\mathcal{M}(\Gamma)\cong\\{F\text{ rooted forest}\ |\
F^{*}\subseteq\Gamma\\},$

where $F^{*}$ denotes the underlying undirected graph of $F$, given by the map
$f\mapsto R(\Gamma,f)$.

###### Proof.

Consider some $f\in\mathcal{M}(\Gamma)$. Let $X\subseteq E_{\Gamma}$ denote
the set of critical edges on $\Gamma$ under $f$. If $Z\subseteq\Gamma$ is a
cycle, then by the second weak Morse inequality,

$c_{1}(f|_{Z})\geq 1=\beta_{1}(Z),$

hence the edge-subgraph of $E_{\Gamma}\setminus X$ in $\Gamma$, say $\Lambda$,
must be a forest. To see that $\Lambda$ has a distinguished root for each
component, let $T$ be a connected component. Then by the first weak Morse
inequality,

$c_{0}(f|_{T})+c_{1}(f|_{T})=1+0=\chi(T)=\beta_{0}(T),$

where the second equality is due to the construction of $\Lambda$. Thus it
follows that each connected component has a sink, therefore $\Lambda$
corresponds to some $R(\Gamma,f)$.

We now construct the inverse map. Without loss of generality888We repeat the
same process on each component for forests having more than one connected
component., assume $F=T$ is tree rooted at a vertex $v$. Denote by
$\widetilde{T}$ the graph resulting from adding the edges
$X:=E_{\Gamma}\setminus E_{T}$ to $T$ and associated additional vertices. We
define a discrete Morse function $g$ on on $\widetilde{T}$ as follows: define
$g(v)=0$, and increment the assignment to cells by $1$ while traversing
outwards along edges from the root (“against the grain”). For any edges $e\in
X$, if $v,w$ are the endpoints of $e$, define $g(e)=\max\\{g(v),g(w)\\}+1$.
The resulting assignment defines a discrete Morse function on $\widetilde{T}$
that trivially satisfies $R(\widetilde{T},g)=T$. Since we clearly have
$\widetilde{T}=\Gamma$, the proof is complete. ∎

In particular, this above map can be refined to a collection of the following
maps.

######  Corollary 3.3.

For a graph $\Gamma$, we have a bijections

$\mathcal{M}_{k}(\Gamma)\cong\\{F\text{ rooted forest}\ |\
F^{*}\subseteq\Gamma\text{ and }|E_{F}|=|E_{\Gamma}|-k\\},$

for each $0\leq k\leq|E|$, given by the map $f\mapsto R(\Gamma,f)$.

It is worth noting that by the first weak Morse inequality, we necessarily
have that $\mathcal{M}_{k}(\Gamma)$ can be non-trivial only when
$\beta_{1}(\Gamma)\leq k\leq|E|$. We highlight the interesting minimal case.

######  Corollary 3.4.

We have that for a graph $\Gamma$ on $n$ vertices that

(*) $|\mathcal{M}_{\beta_{1}}(\Gamma)|=\lambda_{1}\cdots\lambda_{n},$

where the $\lambda_{i}$ are the non-zero eigenvalues of the graph Laplacian of
$\Gamma$.

###### Proof.

Observe that

$\displaystyle|\mathcal{M}_{\beta_{1}}(\Gamma)|$ $\displaystyle=|\\{F\text{
rooted forest}\ |\ F^{*}\subseteq\Gamma\text{ and
}|E_{F}|=|E_{\Gamma}|-\beta_{1}(\Gamma)\\}|$ $\displaystyle=|\\{T\text{ tree}\
|\ \text{$T$ rooted spanning tree of $\Gamma$}\\}|$
$\displaystyle=|V_{\Gamma}|\cdot|\\{T\ |\ \text{$T$ spanning tree of
$\Gamma$}\\}|$
$\displaystyle=|V_{\Gamma}|\cdot\left(\frac{1}{|V_{\Gamma}|}\lambda_{1}\cdots\lambda_{n}\right)$
by Kirchhoff’s Theorem $\displaystyle=\lambda_{1}\cdots\lambda_{n}$

as desired. ∎

######  Remark 3.

It is worth drawing a connection between this corollary and the Morse
(simplicial) complex, an object of current interest in discrete Morse theory.
The Morse complex of $\Gamma$, $\mathfrak{M}(\Gamma)$, is the simplicial
complex whose $n$-simplices are precisely the acyclic $n$-matchings of
$\mathcal{H}_{\Gamma}$. (One can of course talk of the Morse complex of a
general simplicial complex.) It was shown [3, Proposition 10] that the the
number of simplicies of maximal dimension in $\mathfrak{M}(\Gamma)$, for a
graph $\Gamma$ on $n$ vertices, is $($$\dagger$ ‣ 1$)$ as above. Corollary 3.4
proves this statement, as these cells are precisely in correspondence with the
elements of $\mathcal{M}_{\beta_{1}}(\Gamma)$. For more on the Morse complex,
see [26, Chapter 7].

Figure 6. All the possible 16 minimal critical edge configurations on the
wheel graph $W_{4}$. Note that the black edges represent the remainders, hence
every graph in the above figure represents $4$ different gradients; one for
every different vertex considered as a root.

Observe that when $f$ is a discrete Morse function on a graph $\Gamma$ where
$c_{1}(f)=\beta_{1}+k$ for some $0\leq k\leq|E|-\beta_{1}(\Gamma)$, we have
that the number of forests is given by

$\beta_{0}(R(\Gamma,f))=c_{0}(f|_{R(\Gamma,f)})+c_{1}(f|_{R(\Gamma,f)})+\beta_{1}(R(\Gamma,f))=c_{0}(f|_{R(\Gamma,f)}),$

where the first equality is due to the second weak Morse inequality. Note
that, in general, $c_{0}(f|_{R(\Gamma,f)})\leq c_{0}(f)$, as the LHS counts
the number of critical vertices that are not embedded into the edges that are
critical.

We can now present the main result in the $1$-dimensional case.

######  Theorem 3.5.

If $\Gamma$ is a graph, then

$\det(\Delta+\lambda
I)=\sum_{i=0}^{n}|\mathcal{M}_{m-i}(\Gamma)|\lambda^{n-i},$

where $n=|V_{\Gamma}|$ and $m=|E_{\Gamma}|$.

###### Proof.

We have by Corollary 3.3 that all the elements of
$\mathcal{M}_{\beta_{1}+k}(\Gamma)$ are edge-subgraphs of $\Gamma$ with
$(n-k-1)$-many edges and are forests. Each forest $F=R(\Gamma,f)$ for
$f\in\mathcal{M}_{\beta_{1}+k}(\Gamma)$ is repeated $|V_{F}|$-many times, one
for each choice of a root, as accounted for. Lastly, observe that since we
wish for $n-k-1=i$, this implies $k=n-i-1$, and using the formula for the
Euler characteristic of a graph to simplify this expression, we may obtain

$\beta_{1}+k=\beta_{1}+n-i-1=m-i.$

After applying Theorem 2.4, the proof is complete. ∎

It can be notably observed from the above formulation that the minimal
(nonzero) term of $\det(\Delta+\lambda I)$ is the coefficient
$|\mathcal{M}_{0}(\Gamma)|$ of $\lambda^{n-m}=\lambda^{\chi(\Gamma)}$, as
$\mathcal{M}_{k}(\Gamma)=\emptyset$ for $k<0$.

######  Corollary 3.6.

The number of gradients $|\mathcal{M}(\Gamma)|$ on a graph $\Gamma$ is
$\det(\Delta+I)$.

######  Remark 4.

Theorem 3.5 has been previously obtained though in a different context: recall
the definition of the Morse complex from Remark 3. Chari and Joswig showed in
[9, Corollary 3.3] that the $f$-vector of $\mathfrak{M}(\Gamma)$ satisfies

$\det(\Delta+\lambda I)=\sum_{i=0}^{n}f_{m-i}\lambda^{n-i},$

which is exactly the content of the above theorem, as the component $f_{k}$
counts the acyclic $k$-matchings, hence corresponds to gradients having $m-k$
critical cells.

This statement connecting the Laplacian and the $f$-vector of the Morse
complex sadly does not hold in higher dimensions: as we will soon see, the
coeffecients of the characteristic polynomial of the $d$th Laplacian
$\Delta_{d}$ for a $d$-dimensional simplicial complex $K$ do not correspond
with acyclic matchings on all of $\mathcal{H}_{K}$, but rather just the two
top-most levels.

######  Example 3.7.

Consider the cycle $C_{4}$. The graph Laplacian matrix is

$\Delta=\begin{pmatrix}2&-1&0&-1\\\ -1&2&-1&0\\\ 0&-1&2&-1\\\
-1&0&-1&2\end{pmatrix}$

and so the characteristic polynomial is
$\lambda^{4}+8\lambda^{3}+20\lambda^{2}+16\lambda$. Observe the below figure
for the collection of all gradients, illustrating the content of the above
theorem.

Figure 7. The subsets $\mathcal{M}_{k}(C_{4})$ for $\beta_{1}(C_{4})\leq
k\leq|E_{C_{4}}|$. The dashed purple edges denote critical edges, the filled
red circles denote critical vertices that are present in the the remainder,
and the empty red circles denote the other critical vertices which are not
present in the remainder.

######  Example 3.8.

Consider the star graph on $n+1$ vertices, $S_{n}$. Since
$\beta_{1}(S_{n})=0$, we are not “obligated” to put any critical edges
anywhere, we merely choose any edge, if any, to be critical. In the case of
$|\mathcal{M}_{k}(S_{n})|$ for instance: we have all $\binom{n}{k}$-many
choices for a critical edge. For every choice, the remainder requires a root,
thus there are $(n+1-k)$-many vertices left to choose from. Therefore

$|\mathcal{M}_{k}(S_{n})|=\binom{n}{k}(n+1-k),$

and so by Theorem 3.5, the characteristic polynomial of the Laplacian of
$S_{n}$ is given by

$\det(\Delta+\lambda I)=\sum_{i=0}^{n+1}\binom{n}{i}(i+1)\lambda^{n-i+1}.$

###  3.2. Matching polynomials and the adjacency matrix

A simple, yet surprising connection between the graph Laplacian’s
characteristic polynomial and the characteristic polynomial of the adjacency
matrix results from Theorem 3.5. Recall that for a graph $\Gamma$, an
$r$-matching is a matching having $r$-many elements. Let us denote the number
of $r$-matchings on a graph $\Gamma$ by $m_{r}$.

######  Definition 3.9.

The matching polynomial of a graph $\Gamma$ with $n$ edges is given by

$m_{\Gamma}(x)=\sum_{k\geq 0}(-1)^{k}m_{k}x^{n-2k}.$

As corollary to Theorem 3.5, we can cross a bridge built in [25], proving that
the matching polynomial and the characteristic polynomial of the adjacency
matrix are equal if and only if the given graph is a forest. As gradients on a
forest $F$ are equivalent to acyclic matchings on the Hasse diagram
$\mathcal{H}_{F}$, which is certainly a forest as it is isomorphic to the
subdivision of $F$, we obtain the following. Let $A_{\Gamma}$ and
$\Delta_{\Gamma}$ denote the adjacency matrix and the Laplacian matrix of a
graph $\Gamma$, respectively.

######  Theorem 3.10.

Let $F$ be a forest graph on $n$-many vertices having $m$-many edges, and
write $\det(\Delta_{F}+\lambda I)=\sum_{i=\beta_{0}}^{n}a_{i}\lambda^{i}$.
Then the characteristic polynomial of the adjacency matrix of
$\mathcal{H}_{F}$ is

$\det(A_{\mathcal{H}_{F}}+\lambda I)=\sum_{k=0}^{m}a_{n-k}\lambda^{n-2k}.$

###### Proof.

It follows from the above discussion that by [25],
$m_{\mathcal{H}_{F}}(\lambda)=\det(A_{\mathcal{H}_{F}}-\lambda I)$. Because
$\mathcal{H}_{F}$ is a forest, all matchings are acyclic, so by Theorem 3.5,
$m_{k}=|\mathcal{M}_{m-k}(\Gamma)|=a_{n-k}$. ∎

##  4\. Gradients of General Simplicial Complexes

Consider now some discrete Morse function $f\in\mathcal{M}(K)$ where $K$ is a
$d$-dimensional simplicial complex. We now associate a pair $(F,R)$, where $F$
is the set999To save time and avoid unnecessary over-precision, we sometimes
use $F$ to also denote the complex generated by the set of non-critical
$d$-cells of $K$. Respectively with the root $R$. of non-critical $d$-cells of
$K$, and $R$ denotes the set of critical $(d-1)$-cells of $F\subseteq K$. As
one expects from the notation, we will show that $(F,R)$ is a rooted forest of
$K$, and this is indeed only dependent on the gradient, so the definition is
well-defined.

###  4.1. Characterizing gradients

######  Theorem 4.1.

If $f\in\mathcal{M}_{\ell}(K)$, then the pair $(F,R)$ is a rooted forest of
$K$.

###### Proof.

Since $F$ is a subcomplex of $K$, $f|_{F}$ is a discrete Morse function on
${F}$, and by definition of $F$ we have $c_{d}\left(f|_{F}\right)=0$. As
$H_{d}(F)=\ker\partial_{d}^{F}$ since $F_{d+1}=\emptyset$, and because
$c_{d}(f|_{F})\geq\beta_{d}(F)=|H_{d}(F)|$ by the first weak Morse inequality,
we obtain that $\ker\partial_{d}^{F}=0$ so $F$ is a forest.

We now show that $R$ is a root of ${F}$. If $F_{d}=\emptyset$, then
$R=\emptyset$, and we are trivially done. Assume moving forward that $F_{d}$
is non-empty.

$R$ is relatively-free: Observe that by Lemma 2.8, we have that for all
$\gamma\in F_{d}$, there exists a unique
$\alpha\in\overline{R}={F}_{d-1}\setminus R$ such that $\alpha<\gamma$. Hence
$\alpha\in\text{supp}\ \partial_{d}^{F}(\gamma)$. Lastly, as $C_{d-1}(F)\cong
C_{d-1}(R)\oplus C_{d-1}(\overline{R})$, this implies $\text{im
}\partial_{d}^{F}\cap C_{d-1}(R)=\\{0\\}$, as every image of a non-zero
element of $C_{d}(F)$ must have a non-trivial $C_{d-1}(\overline{R})$
component.

$R$ is relatively-generating: Consider $\overline{R}$ and note that it is non-
empty since $c_{i}(f)=|K_{i}|$, and since $c_{d}(f)=|K_{d}|-|F_{d}|$, we must
necessarily have that $c_{d-1}(f)=|K_{d-1}|-|F_{d}|$ by the weak Morse
inequality. Now let $\beta\in\overline{R}$. We then have that $\beta$ is not
critical under $f$ on $K$ (and ${F}$). Since $c_{d-2}(f)=|K_{d-2}|$, by Lemma
2.8 it follows that for all $\tau\in K_{d-2}$ that there is no $(d-1)$-cell
$\alpha>\tau$ such that $f(\alpha)\leq f(\tau)$. Yet as $\beta$ is not
critical, by Lemma 2.8 applied to $\beta$, due to the above condition on the
elements of $K_{d-2}$, we have that there exists some $d$-cell $\gamma>\beta$
such that $f(\gamma)\leq f(\beta)$.

But then $\gamma$ is not a critical $d$-cell by definition, hence it must be
an element of $F_{d}$ and a cell of $F$. Since $\dim F=d$, we have that
$(\beta,\gamma)$ is an arrow of the gradient induced by $f$ on $K$ (hence
$F$). As $\beta$ is the only non-critical $(d-1)$-cell on the boundary of
$\gamma$, we must have that the remaining boundary cells are critical.
Therefore, if we write $\partial_{d}^{F}(\gamma)$ as in $(*)$, then there
exists some $0\leq j\leq d$ where
$\gamma|_{[v_{0},\dots,\widehat{v}_{j},\dots,v_{d}]}=\beta$ and

$\beta+\sum_{\begin{smallmatrix}0\leq i\leq d\\\ i\neq
j\end{smallmatrix}}\gamma|_{[v_{0},\dots,\widehat{v}_{i},\dots,v_{d}]}=\partial_{d}^{F}(\gamma)\in\text{im
}\partial_{d}^{F},$

where $\beta\in\overline{R}$ and $\sum_{\begin{smallmatrix}0\leq i\leq d\\\
i\neq
j\end{smallmatrix}}\gamma|_{[v_{0},\dots,\widehat{v}_{i},\dots,v_{d}]}\in
C_{d}(R)$, as needed. ∎

Note that the above result indeed implies much of our $1$-dimensional work, as
if $f$ is a discrete Morse function on a graph $\Gamma$, then we have that
$(F,R)$ is precisely the data of $R(\Gamma,f)$. There is a topological
viewpoint to these rooted forests as well: we necessarily have that $F\cong
R$, that is, $F$ is (simplicially) homotopy equivalent to $R$. This theorem,
along with Corollary 4.3, answer two of the open questions posed by Mukherjee
in [23] simultaneously, namely, classifying all simplicial forests which are
homotopy equivalent to their roots and classifying all rooted forests
satisfying [23, Lemma 3.5].

######  Theorem 4.2.

Let $(F,R)$ be a rooted forest of a $d$-dimensional simplicial complex $K$.
Then $(F,R)$ corresponds to a gradient if and only if $F$ simplicial collapses
to $R$.

###### Proof.

For the forward direction, this follows from one of the main results of
Discrete Morse Theory, as stated by Forman in [14, Theorem 3.3], where we may
collapse along the gradient and preserve the homotopy type of our complex.

Since $F$ simplicial collapses to $R$, there must be a sequence of elementary
collapses

$F=F_{0}\searrow F_{1}\searrow F_{2}\searrow\cdots\searrow F_{n-1}\searrow
F_{n}=R,$

i.e. $F_{i+1}=F_{i}-\\{\sigma,\tau\\}$ where
$(\sigma,\tau)\in(F_{i})_{d-1}\times(F_{i})_{d}$ are a pair of simplicies in
$F_{i}$ such that $\sigma$ is a face of $\tau$ and $\sigma$ has no other
cofaces. Reversing this, we instead obtain a sequence of elementary expansions

($\dagger$) $R=R_{0}\nearrow R_{1}\nearrow R_{2}\nearrow\cdots\nearrow
R_{n-1}\nearrow R_{n}=F$

with $R_{i+1}=R_{i}\cup\\{\sigma,\tau\\}$ with
$(\sigma,\tau)\in(R_{i+1})_{d-1}\times(R_{i+1})_{d}$ a pair of simplices not
in $R_{i}$ where $\sigma$ is a face of $\tau$ and all other faces of $\tau$
are in $R_{i}$. Using this sequence, we can construct an acylic fitting
orientation $\psi:\overline{R}\to F$ on $(F,R)$ as follows: for $0\leq i\leq
n-1$, let $\sigma_{i}$ and $\tau_{i}$ denote the cells added to $R_{i}$ in the
elementary expansion to $R_{i+1}$. Since all faces besides $\sigma_{i}$ of
$\tau_{i}$ are in $R_{i}$, $\sigma_{i}$ must not be in the root $R$. Therefore
let $\psi(\sigma_{i})=\tau_{i}$. As $R$ ultimately expands to $F$, and at each
elementary expansion $R_{i}\nearrow R_{i+1}$ we added a $d$-cell from $F$ and
a $(d-1)$-cell from $\overline{R}$ to $R_{i}$, this implies $\psi$ is a
bijection, as needed.

We now show that $\psi$ is acyclic, i.e. the matching

$M_{\psi}=\\{(\sigma,\tau)\in F_{d-1}\times F_{d}\ |\
\sigma\in\overline{R}\text{ and }\psi(\sigma)=\tau\\}$

induced by $\psi$ on $\mathcal{H}_{K}$ is acyclic. Assume for contradiction
that this is not the case. We then have a directed sequence

$\sigma_{k_{0}},\ \tau_{k_{0}},\ \sigma_{k_{1}},\ \tau_{k_{1}},\ \dots,\
\sigma_{k_{m}},\ \tau_{k_{m}},$

for some $2\leq m\leq n$, where $\psi(\sigma_{k_{i}})=\tau_{k_{i}}$ and
$\sigma_{k_{i+1}}<\tau_{k_{i}}$ for $0\leq i\leq m$ and
$\sigma_{k_{0}}<\tau_{k_{m}}$. Consider $\\{\sigma_{k_{m}},\tau_{k_{m}}\\}$,
and note by construction of $\psi$, there exists some $0\leq j\leq n$ in the
sequence ($\dagger$ ‣ 4.1) such that
$R_{j}=R_{j-1}\cup\\{\sigma_{k_{m}},\tau_{k_{m}}\\}$. Since
$\sigma_{k_{m}}<\tau_{k_{m}-1}$ and $\tau_{k_{m-1}}=\psi(\sigma_{k_{m-1}})$,
we must have that $R_{j+1}=R_{j}\cup\\{\sigma_{k_{m-1}},\tau_{k_{m-1}}\\}$.
Continuing this process, we obtain the subsequence of ($\dagger$ ‣ 4.1) of
elementary expansions

$R_{j}\nearrow R_{j+1}\nearrow R_{j+2}\nearrow\cdots\nearrow
R_{j+(m-1)}\nearrow R_{j+m},$

where $R_{j+i}=R_{j+(i-1)}\cup\\{\sigma_{k_{m-i}},\tau_{k_{m-i}}\\}$ for
$0\leq i\leq m$. But we have that $\sigma_{k_{0}}<\tau_{k_{m}}$, so we also
must have that $R_{j+m+1}=R_{j+m}\cup\\{\sigma_{k_{m}},\tau_{k_{m}}\\}$. Yet
this is not a valid elementary expansion, as these are simplicies that are
already in $R_{j+m}$. Therefore this directed sequence cannot exist, hence
$M_{\psi}$ is acyclic, and so $\psi$ corresponds to a gradient. ∎

The backwards direction in the above proof also shows that if $(F,R)$
corresponds to a gradient, then that gradient is necessarily unique, since
there was only one possible fitting orientation we could construct. We prove
this differently and explicitly in Lemma 4.4.

Figure 8. An example of a sequence of elementary expansions from a root $R$ to
a forest $F$ giving rise to a discrete gradient vector field on the simplicial
complex generated by $F$.

######  Corollary 4.3.

With the same notation as above, if $(F,R)$ is a rooted forest that
corresponds to a discrete gradient vector field, then $H_{d-1}(F,R)\cong 0$.

###### Proof.

We have that $H_{d-1}(F,R)\cong H_{d-1}(F/R)$ by [16, Theorem 2.27] and [16,
Proposition 2.22]. By Theorem 4.2, $F$ is homotopy equivalent to $R$, hence
$F/R$ is a contractible topological space. Therefore $H_{d-1}(F/R)\cong 0$. ∎

Note the converse of Corollary 4.3 does not hold in general (see Example 4.10
below). It was proven by Mukherjee in [24, Theorem 4.2] that the converse does
indeed hold when the simplicial complex is a triangulation of an orientable
$d$-manifold.

Observe now that we may rewrite the sum of Theorem 2.5 as

($*$) $\det(\Delta_{d}+\lambda
I)=\sum_{i=0}^{n}\sum_{\begin{smallmatrix}\text{$(F,R)$ rooted forest}\\\
\text{with $|R|=i$}\end{smallmatrix}}|H_{d-1}(F,R)|^{2}\lambda^{i}$

where $n=|K_{d-1}|$. We now restrict our view to when $K$ is a triangulation
of an orientable $d$-manifold, where we are aided by the work of Mukherjee,
who closely studied the topology of higher dimensional rooted forests in this
case. Mukherjee showed in [24, Theorem 4.5] that for all rooted forests
$(F,R)$ of $K$, $F$ simplicial collapses to $R$, hence by Corollary 4.3, we
may take all the homological weight terms to be $1$. Therefore $(*)$ reduces
to $\sum_{i=1}^{n}q_{i}\lambda^{n-i}$, where $q_{i}$ is the number of rooted
forests $(F,R)$ of $K$ with $|R|=i$.

What then remains to be shown is a bijection between $\mathcal{M}_{m-i}(K)$
and rooted forests of $K$ having $i$-many $(d-1)$-cells in the root. Due to
[23, Corollary 3.7], we have that this map is surjective when $K$ is a
triangulation of an orientable $d$-manifold. In particular, it is shown that
there exist at least $|H_{d-1}(F,R)|$-many gradients giving rise to the rooted
forest $(F,R)$. We show that this map is injective by proving for an arbitrary
simplicial complex $K$, if a rooted forest comes $(F,R)$ from a gradient
$\varphi$, there cannot exist a distinct gradient $\widetilde{\varphi}$ giving
rise to the same rooted forest.

######  Lemma 4.4.

If $K$ is a simplicial complex, then the map

$\mathcal{M}_{m-i}(K)\ \rightarrow\ \left\\{(F,R)\text{ rooted forest of $K$}\
|\ |R|=i\right\\},$

with $m=|K_{d}|$ and $1\leq i\leq|K_{d-1}|$, is injective.

###### Proof.

Begin by assuming for contradiction that it does not hold. Then there exists
two discrete Morse functions
$\varphi,\widetilde{\varphi}\in\mathcal{M}_{m-i}(K)$ such that they induce the
same rooted forest $(F,R)$. By the definition of rooted forests induced from
discrete Morse functions, we have that $\varphi$ and $\widetilde{\varphi}$
have the same critical cells in dimensions $d$ and $d-1$. In other words, the
corresponding subsets $M_{\varphi},M_{\widetilde{\varphi}}\subseteq
E(\mathcal{H}_{K})$ acyclically match the subsets $F\subseteq K_{d}$ and
$\overline{R}\in K_{d-1}$ two different ways. We claim this is impossible.

Let $\overline{R}=\\{\sigma_{1},\dots,\sigma_{k}\\}$,
$F=\\{\tau_{1},\dots,\tau_{k}\\}$ for some $k\in\mathbb{N}$, and have
$\psi,\widetilde{\psi}$ denote the fitting orientations corresponding to the
gradients of $\varphi,\widetilde{\varphi}$, respectively. Assume without loss
of generality that $\psi(\sigma_{j})=\tau_{j}$ for $1\leq j\leq k$. It then
follows that $\widetilde{\psi}(\sigma_{j})=\tau_{{\pi(j)}}$ for some
permutation $\pi\in S_{k}$. Note that $\pi\neq e$, as otherwise
$\varphi=\widetilde{\varphi}$. We may then rewrite $\pi$ as a product of
disjoint cycles $\gamma_{1}\gamma_{2}\cdots\gamma_{n}$, and consider a cycle
having length greater than $1$, say (without loss of generality)
$\gamma_{1}=(\alpha_{1}\ \alpha_{2}\ \alpha_{3}\cdots\alpha_{m})$. Thus we
have that

$\widetilde{\psi}(\sigma_{\alpha_{1}})=\tau_{\alpha_{2}},\
\widetilde{\psi}(\sigma_{\alpha_{2}})=\tau_{\alpha_{3}},\ \dots,\
\widetilde{\psi}(\sigma_{\alpha_{m-1}})=\tau_{\alpha_{m}},\
\widetilde{\psi}(\sigma_{\alpha_{m}})=\tau_{\alpha_{1}}.$

But note that we also have that $\sigma_{\alpha_{j}}<\tau_{\alpha_{j}}$ for
$1\leq j\leq m$ under the poset relation in $\mathcal{H}_{K}$. Therefore there
is in-fact a directed cycle in the fitting orientation $\widetilde{\psi}$:

$\sigma_{\alpha_{1}},\
\widetilde{\psi}(\sigma_{\alpha_{1}})=\tau_{\alpha_{2}},\
\sigma_{\alpha_{2}},\widetilde{\psi}(\sigma_{\alpha_{2}})=\tau_{\alpha_{3}},\
\dots,\ \sigma_{\alpha_{m}},\
\widetilde{\psi}(\sigma_{\alpha_{m}})=\tau_{\alpha_{1}},\
\sigma_{\alpha_{1}}.$

But this contradicts that $M_{\widetilde{\varphi}}$ is an acyclic matching,
therefore the map must be injective. ∎

With the above work and Lemma 4.4, we have proven the following result.

######  Theorem 4.5.

If $K$ is a triangulation of an orientable $d$-manifold, then

($\dagger$) $det(\Delta_{d}+\lambda
I)=\sum_{i=0}^{n}|\mathcal{M}_{m-i}(K)|\lambda^{n-i},$

where $m=|K_{d}|$ and $n=|K_{d-1}|$.

Note that we’ve also thus shown the following more general statement as well.

######  Theorem 4.6.

If $K$ is a simplicial complex such that for every rooted forest $(F,R)$, the
forest $F$ simplicial collapses to the root $R$, then equation $($$\dagger$ ‣
1$)$ is satisfied.

We also have a good understanding of how the gradient generating function in
$($$\dagger$ ‣ 1$)$ behaves when it is not equal to $\det(\Delta_{d}+\lambda
I)$. As the elements of $\mathcal{M}_{i}(K)$ correspond to a unique rooted
forest by Lemma 4.4 for arbitrary complexes, the main issue is that not every
rooted forest may be induced by a gradient, i.e. surjectivity of the map in
the lemma. Therefore the below statement follows from our above work and
Theorem 2.5.

######  Porism 4.7.

Let $K$ be an arbitrary $n$-dimensional simplicial complex. Then, for $0\leq
i\leq n$, the coefficient of $\lambda^{i}$ in
$\sum_{i=0}^{n}|\mathcal{M}_{m-i}(K)|\lambda^{n-i}$ is less than or equal to
the coefficient of $\lambda^{i}$ in $det(\Delta_{d}+\lambda I)$. The
difference $\varepsilon_{i}$ between the coefficients is

$\varepsilon_{i}=\sum_{\Lambda_{i}(K)}|H_{d-1}(F,R)|^{2}$

where the set $\Lambda_{i}(K)$ consists of the rooted forests $(F,R)$ of $K$
with $|R|=i$ where $F$ does not simplicial collapse to $R$.

######  Corollary 4.8.

Let $K$ be a simplicial triangulation of a non-orientable $d$-manifold. Then
$\varepsilon_{s}>0$, where $s=|K_{d-1}|-|K_{d}|$.

###### Proof.

By [23, Proposition 3.8], if $K$ is the triangulation of a non-orientable
$d$-manifold, we then have that all of $F=K_{d}$ is a forest that contains
torsion in $H_{d-1}(F)$. Let $R$ be a root of $F$. Then $H_{d-1}(R)$ does not
contain torsion due to being the top homology group of $R$, so $F$ cannot
simplicial collapse to $R$. Thus by Theorem 4.2, $(F,R)$ does not correspond
to a gradient. By [6, Lemma 14], $(F,R)$ is a rooted forest a simplicial
complex $K$ if and only if $|\overline{R}|=|F|$ and $H_{d}(F,R)=0$. Since
$F=K_{d}$, we have that $|R|=|K_{d-1}|-|K_{d}|$, thus $\varepsilon_{s}>0$. ∎

######  Example 4.9.

Let $K$ be the simplicial complex that arises from the triangulation of the
Möbius strip that appears in Figure 4. We have that

$\displaystyle det(\Delta_{d}+\lambda I)\quad$ $\displaystyle=\quad$
$\displaystyle\lambda^{10}+15\lambda^{9}+85\lambda^{8}+225\lambda^{7}+275\lambda^{6}+125\lambda^{5},$
$\displaystyle\sum_{i=0}^{n}|\mathcal{M}_{m-i}(K)|\lambda^{n-i}\quad$
$\displaystyle=\quad$
$\displaystyle\lambda^{10}+15\lambda^{9}+85\lambda^{8}+225\lambda^{7}+275\lambda^{6}+121\lambda^{5}.$

In this case $\Lambda_{5}(K)$ is non-empty, and has precisely the one rooted
forest $(F,R)$: the case where $F=K_{2}$ are all the $2$-cells and $R$ is the
boundary of the strip, as highlighted in red in Figure 4. Hence we obtain that

$\varepsilon_{5}=|H_{1}(F,R)|^{2}=|\mathbb{Z}/2\mathbb{Z}|^{2}=2^{2}=4.$

Let us illustrate a strange example where we do not have equality in
$($$\dagger$ ‣ 1$)$. The below example is not the triangulation of any
manifold. Though, it is curious as a counterexample, as it is
$\mathbb{Z}$-acyclic in positive codimension ($\mathbb{Z}$-APC), i.e.
$H_{i}(K;\mathbb{Z})=0$ for all $i<\dim K$, which is a property satisfied by
many complexes of combinatorial interest [13].

######  Example 4.10.

Consider the case where $K$ is the equatorial bipyramid, which is the
$2$-dimensional simplicial complex having seven $2$-cells and is the image
under $\partial_{3}$ of the complex generated by $\\{0,1,2,3\\}$ and
$\\{1,2,3,4\\}$. Here, we have that

$\displaystyle det(\Delta_{d}+\lambda I)\quad$ $\displaystyle=\quad$
$\displaystyle\lambda^{9}+21\lambda^{8}+174\lambda^{7}+710\lambda^{6}+1425\lambda^{5}+1125\lambda^{4},$
$\displaystyle\sum_{i=0}^{n}|\mathcal{M}_{m-i}(K)|\lambda^{n-i}\quad$
$\displaystyle=\quad$
$\displaystyle\lambda^{9}+21\lambda^{8}+174\lambda^{7}+710\lambda^{6}+1425\lambda^{5}+1119\lambda^{4}.$

We have that $\Lambda_{4}(K)$ consists of six rooted forests which are all
symmetric to one another, as indicated in Figure 9 below. As each of these
rooted forests $(F,R)$ have trivial $H_{1}(F,R)$, it follows that
$\varepsilon_{4}=|\Lambda_{4}(K)|=6$.

Figure 9. Two of the rooted forests of $\Lambda_{4}(K)$. The two darkened
$2$-cells in both cases indicate simplicies that are not in the forest. The
red $1$-cells are the roots. Cyclically permuting the labeling of the vertices
$1$, $2$, and $3$ around the equator in both shown rooted forests gives rise
to all six elements of $\Lambda_{4}(K)$.

##  5\. Future Directions

We have come to notice through close observations that gradients within
discrete Morse theory behave particularly nicely, e.g. satisfy Theorem 4.5,
when the complex is indeed a triangulation of manifold that may also have had
singularities. These complexes go by the name of simplicial pseudomanifolds,
satisfying the three following properties: for a $d$-dimensional simplicial
complex $K$,

  * •

The closure of $|K|$ equals a union of $d$-cells.

  * •

Each $(d-1)$-cell of $K$ is a face of exactly two $d$-cells of $K$.

  * •

Given two $d$-cells $\sigma,\sigma^{\prime}$ of $K$, there is a sequence of
$d$-cells of $K$

$\sigma=\sigma_{0},\sigma_{1},\dots,\sigma_{k}=\sigma^{\prime}$

such that $\sigma_{i}\cap\sigma_{i+1}$ is an $(d-1)$-cell for each $i$.

We have evidence then to believe the following conjecture is true and intend
on releasing a future work addressing this connection in further detail.

######  Conjecture 1.

Let $K$ be a $d$-dimensional simplicial complex. Then $K$ is an orientable
pseudomanifold if and only if

($\dagger$) $det(\Delta_{d}+\lambda
I)=\sum_{i=0}^{n}|\mathcal{M}_{m-i}(K)|\lambda^{n-i},$

where $m=|K_{d}|$ and $n=|K_{d-1}|$.

For instance, one can observe that the failure of the Möbius strip (Example
4.9) to be an issue with orientability, whereas the failure of the bipyramid
(Example 4.10) as an issue with branching, e.g. the existence of subcomplexes
of the forest(s) pictured in the figure above having three $d$-cells share a
single $(d-1)$-cell. One may also wonder: if the characteristic polynomials
play nicely with gradients in these cases, what else holds true from a
discrete Morse theory perspective of pseudomanifolds?

We have also had many other interesting questions raised to us, which we
invite the reader to explore. For instance, there are many variants of
discrete Morse theory that have similar concepts of gradients. One may instead
focus on general posets [4], consider instead acyclic matchings on free chain
complexes (dubbed algebraic Morse theory) [22], piecewise linear Morse theory
(see [20, Chapter 5]), among others.

######  Question 1.

Does a similar statement to $($$\dagger$ ‣ 1$)$ exist and appear in these
other variants of discrete Morse theory?

Discrete Morse theory has been used in various software relating to
topological data analysis, such as the Perseus Software Project for Rapid
Computation of Persistent Homology and the Computation Homology Project
(CHomP), and one interesting future suggestion is to use this work to pave the
way for a deeper understanding of the structure of gradients and allow for
improvements in the efficiency of these computations, e.g. developing more
optimized algorithms for finding acyclic matchings (see [18]).

######  Question 2.

Given the indicated structure of subsets of gradients as rooted forests, could
this assist in creating new optimized algorithms for determining gradients on
simplicial complexes satisfying particular constraints?

On the note of Remark 1, we are also curious about the following.

######  Question 3.

Is there a matroid interpretation of $($$\dagger$ ‣ 1$)$?

Further, recall that the set of independent sets $\mathcal{I}$ of a matroid
$M$ form a simplicial complex (the independence complex of the matroid).

######  Question 4.

If $M=(E,\mathcal{I})$ is a matroid, what do gradients on $\mathcal{I}$ reveal
about $M$? In the case of a simplicial matroid $M_{K}$ for a complex $K$, how
does $\det(\Delta_{K}+\lambda I)$ relate to $\det(\Delta_{\mathcal{I}}+\lambda
I)$?

The Laplacian we discussed here is only “half” of what is known as the
combinatorial Laplacian, defined as
$\mathcal{L}_{n}:=\partial_{n+1}\partial_{n+1}^{T}+\partial_{n}^{T}\partial_{n}$
for a $d$-dimensional simplicial complex, where $n\leq d-1$. Though we were
unable to find a connection between discrete gradient vector fields and
$\det(\mathcal{L}_{n}+\lambda I)$, it is suspected one may exist. This
formulation of the Laplacian is particularly interesting, as for instance, one
has that $H_{n}(K)\cong\ker(\Delta_{n})$ due to the combinatorial Hodge
decomposition theorem.

Lastly, there exists a connection between discrete Morse theory and Adinkras:
a decorated graph originating from the study of supersymmetric algebras in
particle physics (see [17]). In particular, the height function on an Adinkra
induces discrete Morse functions for surfaces that are supported on Adinkra
graphs.

######  Question 5.

What can be said about Adrinkas with the aid of discrete gradient vector
fields and their connection to the Laplacian?

## References

  * [1] Michael Aschbacher “Combinatorial cell complexes” In Progress in algebraic combinatorics (Fukuoka, 1993) 24, Adv. Stud. Pure Math. Math. Soc. Japan, Tokyo, 1996, pp. 1–80 DOI: 10.2969/aspm/02410001
  * [2] Rafael Ayala, Luis M. Fernández and Josè A. Vilches “Characterizing equivalent discrete Morse functions” In Bull. Braz. Math. Soc. (N.S.) 40.2, 2009, pp. 225–235 DOI: 10.1007/s00574-009-0010-3
  * [3] Rafael Ayala, Luis M. Fernández, Antonio Quintero and Josè A. Vilches “A note on the pure Morse complex of a graph” In Topology Appl. 155.17-18, 2008, pp. 2084–2089 DOI: 10.1016/j.topol.2007.04.023
  * [4] Eric Babson and Patricia Hersh “Discrete Morse functions from lexicographic orders” In Trans. Amer. Math. Soc. 357.2, 2005, pp. 509–534 DOI: 10.1090/S0002-9947-04-03495-6
  * [5] Matthew Baker and Farbod Shokrieh “Chip-firing games, potential theory on graphs, and spanning trees” In J. Combin. Theory Ser. A 120.1, 2013, pp. 164–182 DOI: 10.1016/j.jcta.2012.07.011
  * [6] Olivier Bernardi and Caroline J. Klivans “Directed rooted forests in higher dimension” In Electron. J. Combin. 23.4, 2016, pp. Paper 4.3520 DOI: 10.37236/5819
  * [7] Norman Biggs “Algebraic graph theory”, Cambridge Mathematical Library Cambridge University Press, Cambridge, 1993, pp. viii+205
  * [8] Kenneth A. Brown and Kevin P. Knudson “Nonlinear statistics of human speech data” In Internat. J. Bifur. Chaos Appl. Sci. Engrg. 19.7, 2009, pp. 2307–2319 DOI: 10.1142/S0218127409024086
  * [9] Manoj K. Chari and Michael Joswig “Complexes of discrete Morse functions” In Discrete Math. 302.1-3, 2005, pp. 39–51 DOI: 10.1016/j.disc.2004.07.027
  * [10] Ivan Contreras and Boyan Xu “The graph Laplacian and Morse inequalities” In Pacific J. Math. 300.2, 2019, pp. 331–345 DOI: 10.2140/pjm.2019.300.331
  * [11] Justin Curry, Robert Ghrist and Vidit Nanda “Discrete Morse theory for computing cellular sheaf cohomology” In Found. Comput. Math. 16.4, 2016, pp. 875–897 DOI: 10.1007/s10208-015-9266-8
  * [12] Art M. Duval, Caroline J. Klivans and Jeremy L. Martin “Critical groups of simplicial complexes” In Ann. Comb. 17.1, 2013, pp. 53–70 DOI: 10.1007/s00026-012-0168-z
  * [13] Art M. Duval, Caroline J. Klivans and Jeremy L. Martin “Simplicial and cellular trees” In Recent trends in combinatorics 159, IMA Vol. Math. Appl. Springer, [Cham], 2016, pp. 713–752 DOI: 10.1007/978-3-319-24298-9“˙28
  * [14] Robin Forman “Morse theory for cell complexes” In Adv. Math. 134.1, 1998, pp. 90–145 DOI: 10.1006/aima.1997.1650
  * [15] Robin Forman “A user’s guide to discrete Morse theory” In Sém. Lothar. Combin. 48, 2002, pp. Art. B48c35
  * [16] Allen Hatcher “Algebraic topology” Cambridge University Press, Cambridge, 2002
  * [17] Kevin Iga and Yan X. Zhang “Structural theory and classification of 2D Adinkras” In Adv. High Energy Phys., 2016, pp. Art. ID 398061312 DOI: 10.1155/2016/3980613
  * [18] Michael Joswig and Marc E. Pfetsch “Computing optimal Morse matchings” In SIAM J. Discrete Math. 20.1, 2006, pp. 11–25 DOI: 10.1137/S0895480104445885
  * [19] Matthew Kahle “Random geometric complexes” In Discrete Comput. Geom. 45.3, 2011, pp. 553–573 DOI: 10.1007/s00454-010-9319-3
  * [20] Kevin P. Knudson “Morse theory” Smooth and discrete World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2015, pp. xiv+181 DOI: 10.1142/9360
  * [21] W. Kook, V. Reiner and D. Stanton “Combinatorial Laplacians of matroid complexes” In J. Amer. Math. Soc. 13.1, 2000, pp. 129–148 DOI: 10.1090/S0894-0347-99-00316-1
  * [22] Dmitry N. Kozlov “Discrete Morse theory for free chain complexes” In C. R. Math. Acad. Sci. Paris 340.12, 2005, pp. 867–872 DOI: 10.1016/j.crma.2005.04.036
  * [23] Sajal Kumar Mukherjee “On the topology of rooted forests in higher dimensions” In Topology Appl. 247, 2018, pp. 50–56 DOI: 10.1016/j.topol.2018.07.013
  * [24] Sajal Kumar Mukherjee “On the rooted forests in triangulated closed manifolds” In Linear Multilinear Algebra 68.10, 2020, pp. 2034–2043 DOI: 10.1080/03081087.2019.1570066
  * [25] Horst Sachs “Beziehungen zwischen den in einem Graphen enthaltenen Kreisen und seinem charakteristischen Polynomialen” In Publ. Math. Debrecen 11, 1964, pp. 119–134
  * [26] Nicholas A. Scoville “Discrete Morse theory” 90, Student Mathematical Library American Mathematical Society, Providence, RI, 2019, pp. xiv+273 DOI: 10.1090/stml/090


# A new microscopic representation of the spin dynamics in quantum systems
with biquadratic exchange interactions

Mariya Iv. Trukhanova1,2<EMAIL_ADDRESS>Pavel Andreev1
<EMAIL_ADDRESS>1 Faculty of Physics, Lomonosov Moscow State
University, Leninskie Gory-1, 119991 Moscow, Russia
2 Theoretical Physics Laboratory, Nuclear Safety Institute, Russian Academy of
Sciences, B. Tulskaya 52, 115191 Moscow, Russia

###### Abstract

There is a version of the Landau-Lifshitz equation that takes into account the
Coulomb exchange interactions between atoms, expressed by the term
$\sim\bm{s}\times\triangle\bm{s}$. On the other hand, ions in the magnetic
materials have several valence electrons on the $d$-shell, and therefore the
Hamiltonian of many-electron atoms with spins $S>1$ should include a
biquadratic exchange interaction. We first propose a new fundamental
microscopic derivation of the spin density evolution equation with an explicit
form of biquadratic exchange interaction using the method of many-particle
quantum hydrodynamics. The equation for the evolution of the spin density is
obtained from the many-particle Schrodinger-Pauli equation and contains the
contributions of the usual Coulomb exchange interaction and the biquadratic
exchange. Furthermore, the derived biquadratic exchange torque in the spin
density evolution equation is proportional to the nematic tensor for the
medium of atoms with spin S = 1. Our method may be very attractive for further
studies of the magnetoelectric effect in multiferroics.

##  I Introduction

Interactions between the identical particles lead to the spin dependence of
their energy [1]. The quantum mechanical nature of the Coulomb exchange
interactions was demonstrated by Heisenberg [2], [3]. The description of
ferromagnetic and antiferromagnetic ordering is in good agreement with the
model based on the effective spin Hamiltonian of the exchange interactions of
atoms (ions) [3]

$H=J^{(1)}(\bm{r}_{ij})\bm{S}_{i}\cdot\bm{S}_{j}+J^{(2)}(\bm{r}_{ij})(\bm{S}_{i}\cdot\bm{S}_{j})^{2},$
(1)

where $J^{(1)}$, $J^{(2)}$ are the exchange integrals, and
$\bm{S}_{i},\bm{S}_{j}$ are the spins in the nodes of the lattice. The first
term is an antisymmetric exchange interaction and corresponds to the
Heisenberg Hamiltonian. The second term, called the biquadratic exchange,
needs to be considered for many-electron ions with spins $S>1/2$, and vanishes
in the lowest order in the overlap parameter of the ions wave functions. The
Ising model, see [4], was developed for two- and three-dimensional spin
lattices at low temperatures.

The description of the magnetoelectric effect in multiferroic materials, is
based on exchange interactions [5, 6]. The microscopic explanation of the
ferroelectricity of spin origin is based on the following mechanisms:
symmetric exchange interactions between the neighboring ions [7],
antisymmetric exchange interactions between noncollinear spins [8, 9], and
spin-dependent metal-ligand hybridization [10]. The symmetric exchange
interactions lead to the striction along a crystallographic direction
$\bm{\Pi}_{ij}$ [7], so that the polarization appears as a result of the
inversion symmetry breaking:
$\bm{P}_{ij}\sim\bm{\Pi}_{ij}(\bm{S}_{i}\cdot\bm{S}_{j})$. In the case of
antisymmetric exchange interactions, the polarization is induced under the
influence of spin-orbital coupling [11]. The relativistic corrections to
Anderson’s superexchange between magnetic ions in the presence of spin-orbital
coupling lead to an antisymmetric exchange, also named as Dzyaloshinskii-
Moriya interaction [11, 12], with Hamiltonian
$H_{DM}=\bm{D}_{ij}\cdot\bm{S}_{i}\times\bm{S}_{j}$, where $\bm{D}_{ij}$ is
the Dzyaloshinskii vector. The polarization that appears due to this mechanism
is determined by the expression:
$\bm{P}_{ij}\sim\bm{e}_{ij}\times(\bm{S}_{i}\times\bm{S}_{j})$, where
$\bm{e}_{ij}$ is a unit vector connecting two spins [8].

A phenomenological theory of inhomogeneous ferroelectric magnets was proposed
by M. Mostovoy and describes their thermodynamics and the behavior of
polarization in the magnetic field [13]. In this research, the Ginzburg-Landau
approach was applied to derive the expression for the electric polarization
induced by spiral spin-density wave states. A few years later, Masahito
Mochizuki suggested a microscopic model of the multiferroic perovskite
manganites [14]. In this model, the calculation is carried out for each
individual ion at the site of the crystal lattice. The Landau-Lifshitz-Gilbert
equation was derived for the Mn $3d$-spin atoms in the local magnetic field
from the Heisenberg model with additional interactions and magnetic
anisotropies.

The theory of the ferroelectric polarization which is induced by the spin-
spiral order was developed in Refs.[15, 16]. In particular, a method for
describing the electronic and ion-displacement contributions to the spin-
order-induced ferroelectric polarization was proposed in Ref. [16]. For the
evolution of the electric polarization with time, the nonlinear Duffing
oscillator equation was obtained using the Euler-Lagrange equation [17]. A
phenomenological framework for studying the switching dynamics of coupled
polarization and magnetization in single-phase multiferroic materials was
developed in Ref. [18], where the oscillatory mode of magnetization can act as
a microwave signals.

Understanding the microscopic foundations of the magnetoelectric effect in
magnetically ordered media is an important step in the development of
polarization control [19]. The magnetically induced electric dipole moment of
the cell of the crystal lattice, the collective dynamics of these dipoles, and
the formation of wave patterns are realized in three-dimensional physical
space. The collective effects require representation of the dynamics in terms
of field variables defined in three-dimensional physical space [20]. A model
that describes the dynamics of non-equilibrium systems can be developed using
the method of many-particle quantum hydrodynamics [21] \- [25].

Exchange interactions underlie spin ordering in magnetic structures.
Magnetocrystalline anisotropy, external magnetic field, biquadratic exchange,
and demagnetization field can also influence the formation of magnetic
structures. The basic dynamical equation for magnetization was derived in 1935
by Landau and Lifshitz [26]. Later, using the Lagrangian approach, Hilbert
supplemented the damping term with this equation [27, 28, 29]. The exchange
interaction is electrostatic in nature and is directly related to the Pauli
principle. While the symmetric Heisenberg exchange interaction is well studied
and its contribution to the Landau-Lifshitz equation has long been derived [3,
30], biquadratic exchange is rarely seen in theoretical articles regarding
magnetic phenomena. At the same time, magnetic atoms and ions have several
valence electrons on the $d$-shell, and therefore the Hamiltonian of many-
electron atoms has to include a biquadratic exchange. Since most of the
magnetic ions in type-II multiferroics have spin $S>1/2$ [9, 7], it is
necessary to consider the biquadratic exchange interaction in the equation of
magnetization evolution (spin density).

In this article, we propose a fundamental microscopic derivation of the
magnetization evolution equation, or the Landau-Lifshitz-like equation with
biquadratic exchange interaction. As a starting point, the many-particle
Schrodinger-Pauli equation for a system of particles (ions or atoms) is used.
Further, the exchange interactions between individual electrons are considered
in the form of the Heisenberg-like Hamiltonian of interactions of the magnetic
moments of ions. We derive the spin density evolution equation, which includes
the contribution of the Coulomb exchange interaction between valence electrons
and the biquadratic exchange interactions. For this purpose, we use the method
of many-particle quantum hydrodynamics.

##  II The theoretical model

Ferromagnetic materials have been studied using different theoretical methods
for a long time. In this article, we propose a new approach to learning
collective phenomena in magnetically ordered materials. The many-particle
quantum hydrodynamics method has been developed for several physical systems,
such as quantum plasmas [31], polarized Bose-Einstein condensates [32] and
spin-separated quantum systems [33]. The first step in our approach is the
introduction of the many-particle interaction Hamiltonian of interactions. The
second step is to determine the microscopic representation of the spin density
of the particle system. In order to obtain the equation defining the dynamics
of spin density, it is imperative to identify the macroscopic form of the spin
current density and the torque density.

We consider a system of magnetic ions localized in the nodes of the crystal
lattice and denote by $\bm{S}_{j}$ the spins of the valence electrons of ions.
We construct the model for two kinds of ions: those with spin $S=1/2$ and
$S=1$. But, this model can also approximately describe the behavior of spins
$S>1/2$. Most of the magnetic ions in ferroelectric materials have spins
greater than one. For instance, the perovskite oxide with an orthorhombic
structure, TbMnO3, has magnetic ions Mn3+ with $S$=2, or the cubic spinel
ZnCr2Se4, where Cr3+ has spins $S$=$3/2$, or the mixed-valence manganese oxide
family, RMn2O5, R is Y, Tb, Ho, Er, Tm, with Mn4+ ($S$=3/2). This means that
our theoretical model has to include higher-order exchange interactions.

The first step of our approach is the introduction of the many-particle
interaction Hamiltonian in the form

$\hat{H}=\sum^{N}_{j=1}\left(\frac{\hat{\bm{p}}_{j}^{2}}{2m_{j}}-\gamma_{j}\hat{S}^{\beta}_{j}B^{\beta}_{j}\right)$
$+\frac{1}{2}\sum_{l\neq
j}^{N}\left(U^{(1)}(\bm{r}_{lj})\hat{S}^{\beta}_{l}\hat{S}^{\beta}_{j}+U^{(2)}(\bm{r}_{lj})(\hat{S}^{\beta}_{l}\hat{S}^{\beta}_{j})^{2}\right),$
(2)

where $\hat{p}^{\alpha}_{j}=-i\hbar\partial^{\alpha}_{j}$ is the momentum
operator, $m_{j}$, $\gamma_{j}$ are the mass and gyromagnetic ratio of $j$-th
particle. The first term in the above Hamiltonian (2) is the sum of the
kinetic energies of all particles. The second term is responsible for the
action of the external magnetic field $\bm{B}_{j}$ on the particle’s magnetic
moment. The third and fourth terms characterize Coulomb (symmetric and
biquadratic) exchange interactions. In real magnetic materials and
multiferroics, the interactions described above occur between two neighbouring
spins of magnetic ions located at the nodes of the crystal lattice. For
instance, the positively charged transition metal ions Ni2+ in the case of
Ni3V2O8. In the simplest case of one-electron ions, the exchange interaction
is proportional to the scalar product of spins. On the other hand, the
quadratic correction is to be taken into consideration for many-electron
systems and for the spin-ordering of low-dimensional nanostructures. The
biquadratic exchange determines the magnetic features in multilayer materials
[34], perovskites and some oxides [35, 36], and are important for the magnetic
properties of 2D materials [36].

The generalized Hamiltonian (2) contains effective potentials
$U^{(1)}(\bm{r}_{lj})$ and $U^{(2)}(\bm{r}_{lj})$, which are functions of the
distance between the interacting particles $\bm{r}_{lj}$. These effective
potentials are the exchange integrals. Positive (or negative) values of the
function $U^{(1)}(\bm{r}_{lj})$ characterize the antiferromagnetic or
ferromagnetic equilibrium states. The dipole-dipole interactions between
magnetic moments are not considered in the framework of the many-particle
Hamiltonian (2). In this article, we aim to investigate Coulomb exchange
interactions. But, the method of many-particle quantum hydrodynamics was used
to investigate physical systems with dipole-dipole magnetic interactions (see
Refs. [37, 38]).

Collective processes in many-particle systems are described by macroscopic
physical fields. Let’s study the collective dynamics using the local spin
density of particles as a suitable collective variable. For this purpose the
dynamics of $N$-interacting particles is given by the Pauli-Schrodinger
equation $i\hbar\partial_{t}\psi_{s}(R,t)=\hat{H}\psi_{s}(R,t)$ with
Hamiltonian (2). In the approach of many-particle quantum hydrodynamics, the
many-particle wave function $\psi_{s}(R,t)$ depends on the $3N$ position
variables $R=\\{\bm{r}_{1},...\bm{r}_{N}\\}$ and time. The spin density in the
neighborhood of $\bm{r}$ can be introduced as the quantum average of the spin
density operator. The average value of the operator of the physical quantity
coincides with the eigenvalue of the operator if the quantum state is an
eigenstate of the operator of the physical quantity [24]. The operator of spin
density is chosen as the quantization of the microscopic local spin density,
and is given by
$\hat{\bm{s}}^{\alpha}=\sum_{k}\delta(\bm{r}-\bm{r}_{k})\hat{\bm{S}}_{k}^{\alpha}$.
As a result, the spin density can be introduced in the form

$\bm{s}^{\alpha}(\bm{r},t)=\sum_{s=s_{1},...,s_{N}}\int
dR\sum_{k=1}^{N}\delta(\bm{r}-\bm{r}_{k})\qquad\qquad$
$\qquad\qquad\qquad\qquad\times\psi^{\dagger}_{s}(R,t)\hat{\bm{S}}_{k}^{\alpha}\psi_{s}(R,t),$
(3)

where $dR=\Pi_{j=1}^{N}d\bm{r}_{j}$ is the element of the volume in
$3N$-dimensional configuration space, $N$ is the number of particles, and
$\hat{\bm{S}}_{k}^{\alpha}$ is the spin operator of the $k$-th particle.

For the case of $S=1/2$, the form of the cartesian components of the spin
operator in the cyclic basis is
$\hat{\bm{S}}_{k}=\frac{1}{2}\hat{\bm{\sigma}}_{k}$, where the Pauli matrices
are given by

$\qquad\hat{\sigma}^{x}_{k}=\begin{pmatrix}0&1\\\
1&0\end{pmatrix},\hat{\sigma}^{y}_{k}=\begin{pmatrix}0&-i\\\
i&0\end{pmatrix},\hat{\sigma}^{z}_{k}=\begin{pmatrix}1&0\\\
0&-1\end{pmatrix}.$ (4)

In the case of $S=1$, the spin operator can be represented in the form
$\hat{\bm{S}}_{k}=\hat{\bm{F}}_{k}$, where the spin matrices are chosen as

$\hat{F}^{x}_{k}=\frac{1}{\sqrt{2}}\begin{pmatrix}0&1&0\\\ 1&0&1\\\
0&1&0\end{pmatrix},$
$\hat{F}^{y}_{k}=\frac{1}{\sqrt{2}}\begin{pmatrix}0&-i&0\\\ i&0&-i\\\
0&i&0\end{pmatrix},\;\hat{F}^{z}_{k}=\begin{pmatrix}1&0&0\\\ 0&0&0\\\
0&0&-1\end{pmatrix}$ (5)

In the next step, we derive the dynamical equation for the evolution of the
spin density. For this purpose, we take the time derivative of the spin
density (3). This derivative acts on the many-particle wave function
$\psi_{s}(R,t)$ under the integral. The value of $\dot{\psi}_{s}$ is
determined by the Pauli-Schrodinger equation with Hamiltonian (2). The
resulting equation for the evolution of the spin density is derived in the
form

$\partial_{t}{s}^{\alpha}(\bm{r},t)+\partial_{\beta}j^{\alpha\beta}_{s}(\bm{r},t)\qquad\qquad\qquad$
$\qquad\qquad=\frac{2\gamma}{\hbar}\epsilon^{\alpha\beta\mu}s^{\beta}(\bm{r},t)B^{\mu}(\bm{r},t)+T^{\alpha}(\bm{r},t).$
(6)

Let us now discuss the relation between the equation (6) and the usual Landau-
Lifshitz equation ($LL$). The $LL$ equation describes the evolution of the
magnetization of magnetically ordered materials in the effective magnetic
field and leads to many physical effects [39]. In particular, the precession
of the magnetic moments in the external magnetic field is the simplest
consequence of this equation. This effect is characterized by the first term
on the right side of Eq. (6). The precession of the magnetic moments can be
caused by particle interactions, for instance, dipole-dipole interactions,
which were investigated in Refs. [37, 38]. The Coulomb exchange interaction
considered in this paper is taken into account by the last term on the right
side of Eq. (6). This effect is also called the energy of heterogeneity in
phenomenological consideration. On the other hand, the Landau-Lifshitz-Gilbert
($LLG$) equation takes into account the damping term [3]. But, in this
article, we focus on the microscopic mechanisms of the Coulomb exchange
interactions for $S>1/2$ systems and do not include any damping because there
are no degrees of freedom that can be used as an energy reservoir.

The second term on the left-hand side of Eq. (6) represents the spatial
derivative of the spin current density. The spin current tensor is a product
of two vectors: the charge current density and spin polarization [40]. The
spin-polarized current of electrons determines the spin transfer torque and
leads to the additional precession of the magnetization [41]. In other words,
if the spin current $j^{\alpha\beta}_{s}$ changes from point to point, the
magnetization or the spin density vector experiences spin transfer torque. The
microscopic representation of the spin current density can be derived using
the many-particle quantum hydrodynamics method in the form

$j^{\alpha\beta}_{s}(\bm{r},t)=\sum_{s=s_{1},...,s_{N}}\int
dR\sum_{k=1}^{N}\delta(\bm{r}-\bm{r}_{k})\frac{1}{2m_{k}}$
$\times\left((\hat{p}^{\beta}_{k}\psi_{s})^{\dagger}(R,t)\hat{\bm{S}}^{\alpha}_{k}\psi_{s}(R,t)+h.c.\right).$
(7)

The above expression is not postulated. Details of the derivation of the spin
current density can be found in Refs. [21], [37], [38].

###  II.1 Biquadratic exchange interactions

Let’s discuss Eq. (6) in details. The first term on the right-hand side of the
equation is responsible for the external magnetic field $\bm{B}$. The second
term characterizes the torque acting on the spin density via the Coulomb
exchange interactions $T^{\alpha}(\bm{r},t)$ and includes potentials
$U^{(1)}_{lj}$ and $U^{(2)}_{lj}$. After the calculations, this torque can be
represented in the following form

$T^{\alpha}(\bm{r},t)=\frac{1}{2\hbar}\epsilon^{\alpha\beta\gamma}\sum_{s}\int
dR\sum_{l\neq
j}^{N}\biggl{(}\delta(\bm{r}-\bm{r}_{j})-\delta(\bm{r}-\bm{r}_{l})\biggr{)}$
$\psi^{\dagger}_{s}\left(U^{(1)}_{lj}\hat{S}^{\beta}_{l}\hat{S}^{\gamma}_{j}+\frac{1}{2}U^{(2)}_{lj}\left(4\hat{\pi}^{\beta\mu}_{l}\hat{\pi}^{\gamma\mu}_{j}-\hat{S}^{\beta}_{l}\hat{S}^{\gamma}_{j}\right)\right)\psi_{s},$
(8)

where the nematic tensor operator is derived in the form

$\hat{\pi}^{\gamma\mu}_{l}=\frac{1}{2}\biggl{(}\hat{S}^{\gamma}_{l}\hat{S}^{\mu}_{l}+\hat{S}^{\mu}_{l}\hat{S}^{\gamma}_{l}\biggr{)},$
(9)

here we use the expression
$[\hat{S}^{\alpha}_{l},\hat{S}^{\beta}_{j}]=i\epsilon^{\alpha\beta\gamma}\delta_{lj}\hat{S}^{\gamma}_{l}.$
The higher spin systems support spin-nematic order in addition to the charge
and magnetic ones. The magnetic and nematic orders compete because increasing
one requires lowering the other [42].

Let us derive a more explicit form of the torque vector $\bm{T}$. We consider
the identical particles and introduce variables of the center of the mass for
each pair of particles $\bm{R}_{lj}=\frac{1}{2}(\bm{r}_{l}+\bm{r}_{j})$ and
the relative distance between them $\bm{r}_{lj}=\bm{r}_{l}-\bm{r}_{j}$. The
potentials $U^{(1)}_{lj}(\bm{r}_{lj})$ and $U^{(2)}_{lj}(\bm{r}_{lj})$
decrease rapidly as the relative distance between particles increases
$r_{ij}$. The delta-functions $\delta(\bm{r}-\bm{r}_{j})$,
$\delta(\bm{r}-\bm{r}_{l})$, the $l$-th and $j$-th arguments of the wave
function $\psi_{s}(...,\bf{r}_{l},...,\bf{r}_{j},...,t)$ are subjected to a
Taylor series expansion. For this purpose, the many-particle wave function can
be represented in the form
$\psi_{s}(R,t)=\psi_{s}(\bm{r}_{1},...,\bm{R}_{lj}+\frac{1}{2}\bm{r}_{lj},...,\bm{R}_{lj}-\frac{1}{2}\bm{r}_{lj},...,\bm{r}_{N},t),$
and can be decomposed into a Taylor series by a small parameter $\bm{r}_{lj}$

$\psi_{s}(R,t)=\psi_{s}(R^{\prime},t)+\frac{1}{2}r^{\alpha}_{lj}\partial^{\alpha}_{R,1}\psi_{s}(R^{\prime},t)\qquad\qquad$
$\qquad\qquad\qquad\qquad-\frac{1}{2}r^{\alpha}_{lj}\partial^{\alpha}_{R,2}\psi_{s}(R^{\prime},t),$
(10)

where $R^{\prime}=\bm{r}_{1},...,\bm{R}_{lj},...,\bm{R}_{lj},...,\bm{r}_{N}$,
$\partial_{R,1}^{\alpha}=\frac{\partial}{\partial\bm{R}_{lj}}$ is the partial
derivative with respect to $\bm{R}_{lj}$ at the $l$-th place and
$\partial_{R,2}^{\alpha}=\frac{\partial}{\partial\bm{R}_{lj}}$ is the partial
derivative with respect to $\bm{R}_{lj}$ at the $j$-th place. Integral over
$\bm{r}_{lj}$ and integral over $\bm{R}_{lj}$ are independent, so that we can
separate them explicitly $dR=dR_{N-2}d\bm{R}_{lj}d\bm{r}_{lj}$. As a result,
in the limit of the small relative distance between particles $\bm{r}_{lj}$,
in the second order of decomposition and after the serious of calculations, we
derive the non-trivial solution for the exchange torque

$T^{\alpha}(\bm{r},t)=(\hat{g}^{(1)}+\hat{g}^{(2)})\epsilon^{\alpha\beta\gamma}s^{\beta}(\bm{r},t)\triangle
s^{\gamma}(\bm{r},t)\qquad\qquad$
$\qquad\qquad\qquad\qquad-4\hat{g}^{(2)}\epsilon^{\alpha\beta\gamma}\pi^{\beta\mu}(\bm{r},t)\triangle\pi^{\gamma\mu}(\bm{r},t),$
(11)

where the constants of interactions have the form

$\hat{g}^{(1)}=-\frac{1}{6\hbar}\int\xi^{2}U^{(1)}(\xi)d\bm{\xi},\qquad$
$\qquad\qquad\hat{g}^{(2)}=\frac{1}{12\hbar}\int\xi^{2}U^{(2)}(\xi)d\bm{\xi}.$
(12)

Let’s discuss the physical meaning of terms in the definition (11). For
$U^{(2)}=0$, the first term in the expression (11) is proportional to
$\hat{g}^{(1)}$ and follows from the usual Coulomb exchange interaction
(symmetric exchange) [43]. While the equation for the spin density evolution
(6) in the absence of the external magnetic field and for the case of
$U^{(2)}=0$ has the form of the continuity equation

$\partial_{t}\bm{s}^{\alpha}(\bm{r},t)+\partial_{\mu}\cdot\bm{j}^{\alpha\mu}_{exc}(\bm{r},t)=0,$
(13)

where the expression
$j^{\alpha\mu}_{exc}(\bm{r},t)=-\hat{g}^{(1)}\epsilon^{\alpha\beta\gamma}s^{\beta}(\bm{r},t)\partial^{\mu}s^{\gamma}(\bm{r},t)$
represents the exchange spin current. The above equation is the spin angular-
momentum conservation law and has been derived in Refs. [3, 30] and sometimes
can be represented through the exchange magnetic field
$\partial_{t}\bm{s}=\gamma\bm{s}\times\bm{H}_{exc}$, where
$\bm{H}_{exc}=A_{exc}\triangle\bm{s}$ denotes the exchange field, $A_{exc}$ is
a material constant [44].

The presence of the biquadratic exchange interaction, which is taken into
account by potential $U^{(2)}(\xi)$, makes the situation different. Although
the higher-order exchange interactions have already been studied in the
literature [36, 45, 46], the equation for the evolution of the spin density in
the presence of the biquadratic exchange interactions has been obtained in
this article for the first time. The biquadratic exchange leads to the non-
vanishing exchange spin current

$\bm{j}_{exc}^{{}^{\prime}\alpha\mu}(\bm{r},t)=-\hat{g}^{(2)}\epsilon^{\alpha\beta\gamma}s^{\beta}(\bm{r},t)\partial^{\mu}s^{\gamma}(\bm{r},t)\qquad\qquad$
$\qquad\qquad\qquad+4\hat{g}^{(2)}\epsilon^{\alpha\beta\gamma}\pi^{\beta\mu}(\bm{r},t)\partial^{\mu}\pi^{\gamma\mu}(\bm{r},t),$
(14)

where the density of nematic tensor is derived in the microscopic form

$\pi^{\beta\mu}(\bm{r},t)=\sum_{s}\int
dR\sum_{k=1}^{N}\delta(\bm{r}-\bm{r}_{k})\qquad\qquad$
$\qquad\qquad\qquad\qquad\times\psi^{\dagger}_{s}(R,t)\hat{\pi}^{\beta\mu}_{k}\psi_{s}(R,t).$
(15)

For ions with spin-1/2, the density of the nematic tensor becomes proportional
to the concentration $\pi^{\alpha\beta}=(1/4)\delta^{\alpha\beta}n(\bm{r},t)$.

###  II.2 The discussion of the results

The Landau-Lifshitz equation was first derived from phenomenological
considerations in 1935 and provides the basic equation for the evolution of
magnetization in bulk magnetics [26]. The effect of relativistic interactions,
the damping term, was introduced at a later date by Gilbert, and the
conclusive equation was called the Landau-Lifshitz-Gilbert equation [27]. This
equation allows us to take into account all possible fundamental interactions,
first of all the influence of the external magnetic field, but also the
Coulomb exchange interaction, which is a consequence of the overlap of the
wave functions of electrons [28]. The magnetocrystalline anisotropy field, the
demagnetization field, and biquadratic exchange interactions can also be
included into the equation. We focus our attention specifically on the
biquadratic exchange interaction, since there is no clear derivation of its
effect on the evolution of magnetization. In magnetic atoms and metal ions,
the $d$-shell is partially filled. Therefore, several valence electrons of the
atom, for which the values of spins will be $S>1/2$, participate in the
Coulomb exchange interaction. This means that for most many-electron magnetic
atoms, the Hamiltonian of the exchange interaction between two atoms must be
different from the Heisenberg Hamiltonian. However, in the first
approximation, in most studies, the exchange of only one pair of electrons is
taken into account by the overlapping wave functions of magnetic atoms. As a
result, the authors also use the Heisenberg approximation for many-electron
atoms, but, generally, it is not correct. For magnetic atoms with $S>1/2$,
such an approximation is often insufficient, and it is necessary to take into
account the biquadratic exchange, which we have included in the Hamiltonian of
the system (1). In the energy of magnetic inhomogeneity, an additional term
also appears and corresponds to biquadratic exchange.

On the one hand, the exchange interaction between valence electrons of many-
electron metal atoms has been considered in a large number of articles and has
a long history [3], but to obtain the influence of this exchange on the
evolution of magnetization is a difficult task. In this research, using the
method of many-particle quantum hydrodynamics, we first obtained an explicit
form of the biquadratic exchange interaction in the Landau-Lifshitz-like
equation. We derived the effect of this interaction from the basis of the
theoretical model and introduced a microscopic representation of the nematic
tensor (15) that arises in the model as a result of the derivation of the
equation.

##  III Conclusions

In this article, we present a new fundamental microscopic derivation of the
spin density evolution equation (6) with torque (11) for a magnetically
ordered material. The equation (6) is derived as part of the system of
equations for the evolution of the macroscopic characteristics of the medium.
Although we consider the magnetic ions in the nodes of the crystal lattice,
the structure of the equation (6) is analogous to the equation of
hydrodynamics. We aim to study the influence of Coulomb exchange interactions
on spin density dynamics. We focus on the microscopic derivation of the
exchange interactions using the method of many-particle quantum hydrodynamics.
As a result of derivation, the biquadratic exchange interaction accounting
leads to the torque in the equation for the evolution of the spin density.

The first step in describing the magnetoelectric effect in multiferroics is a
new model that takes into account the effect of exchange interactions on the
dynamics of spin density. Explicitly accounting for the impact of biquadratic
exchange is the first major step in this direction. Many-particle quantum
hydrodynamics allows us to study the dynamics of magnetization and
polarization in systems with exchange interactions. Moreover, this method
allows us, from first principles, to take into account all the necessary
interactions between particles and obtain microscopic definitions of the main
macroscopic observables, such as spin density, magnetization, and electric
polarization. This unique feature makes this method very attractive for
further studies of the magnetoelectric effect in multiferroics.

###### Acknowledgements.

The research of Mariya Iv. Trukhanova is supported by the Russian Science
Foundation under the grant No. 22-72-00036
(https://rscf.ru/en/project/22-72-00036/).

## References

  * [1] J. Stohr and H. C. Siegmann, Magnetism, From Fundamentals to Nanoscale Dynamics (Springer-Verlag Berlin Heidelberg, 2006).
  * [2] W. Heisenberg, Zur Theorie des Ferromagnetismus, Z. Phys. A 49, 619 (1928).
  * [3] V. S. L’vov, Nonlinear Spin Waves (Nauka, Moscow, 1987).
  * [4] Kei Yosida, Theory of magnetism (Springer, 1996).
  * [5] Shuai Dong, Jun-Ming Liu, Sang-Wook Cheong and Zhifeng Ren, Multiferroic materials and magnetoelectric physics: symmetry, entanglement, excitation, and topology, Advances in Physics 64(5), 519 (2015).
  * [6] L. D. Landau and E. Lifshitz, Electrodynamics of Continuous Media (Butterworth-Heinemann, Oxford, 2002).
  * [7] Y. Tokura, S. Seki and N. Nagaosa, Multiferroics of spin origin, Rep. Prog. Phys. 77(7), 076501 (2014).
  * [8] I. A. Sergienko and E. Dagotto, Role of the Dzyaloshinskii-Moriya interaction in multiferroic perovskites, Phys. Rev. B 73, 094434 (2006).
  * [9] H. Katsura, N. Nagaosa and A. V. Balatsky, Spin current and magnetoelectric effect in noncollinear magnets, Phys. Rev. Lett. 95, 057205 (2005).
  * [10] T. Arima, Ferroelectricity induced by proper-screw type magnetic order, J. Phys. Soc. Jpn. 76, 073702 (2007).
  * [11] I. Dzyaloshinsky, A thermodynamic theory of “weak” ferromagnetism of antiferromagnetics, J. Phys. Chem. Solids. 4, 241 (1958).
  * [12] T. Moriya, New mechanism of anisotropic superexchange interaction, Phys. Rev. Lett. 4, 228 (1960); T. Moriya, Anisotropic superexchange interaction and weak ferromagnetism, Phys. Rev. 120, (1960).
  * [13] M. Mostovoy, Ferroelectricity in Spiral Magnets, Phys. Rev. Lett. 96, 067601 (2006).
  * [14] M. Mochizuki, and N. Furukawa, Microscopic model and phase diagrams of the multiferroic perovskite manganites, Phys. Rev. B 80(13), 134416 (2009).
  * [15] H.J. Xiang, E.J. Kan, Y. Zhang, M.H. Whangbo and X.G. Gong, General theory for the ferroelectric polarization induced by spin-spiral order, Phys. Rev. Lett. 107, 157202 (2011).
  * [16] H.J. Xiang, P.S. Wang, M.H. Whangbo and X.G. Gong, Unified model of ferroelectricity induced by spin order, Phys. Rev. B 88, 054404 (2013).
  * [17] A. K. Bandyopadhyay, P. C. Ray and Venkatraman Gopalan, Dynamical systems analysis for polarization in ferroelectrics, J. Appl. Phys. 100, 114106 (2006).
  * [18] Kuntal Roy, Dynamical systems study in single-phase multiferroic materials, EPL 108, 67002 (2014).
  * [19] Y. J. Choi, C. L. Zhang, N. Lee and S. W. Cheong, Cross-control of magnetization and polarization by electric and magnetic fields with competing multiferroic and weak-ferromagnetic phases, Phys. Rev. Lett. 105, 097201 (2010).
  * [20] P. A. Andreev, L. S. Kuz’menkov, Problem with the single-particle description and the spectra of intrinsic modes of degenerate boson-fermion systems, Phys. Rev. A 78, 053624 (2008).
  * [21] Pavel A. Andreev, Extended hydrodynamics of degenerate partially spin polarized fermions with short-range interaction up to the third order by interaction radius approximation, Laser Phys. 31, 045501 (2021).
  * [22] Pavel A. Andreev, I. N. Mosaki, and Mariya Iv. Trukhanova, Quantum hydrodynamics of the spinor Bose–Einstein condensate at non-zero temperatures, Physics of Fluids 33, 067108 (2021).
  * [23] Pavel A. Andreev, Quantum hydrodynamic theory of quantum fluctuations in dipolar Bose–Einstein condensate, Chaos 31, 023120 (2021).
  * [24] P. A. Andreev, L. S. Kuz’menkov, Many-particle quantum hydrodynamics: Basic principles and fundamental applications, arXiv:1407.7770.
  * [25] Pavel A. Andreev, Hydrodynamic model of a Bose–Einstein condensate with anisotropic short-range interaction and bright solitons in a repulsive Bose–Einstein condensate, Laser Phys. 29(3), 035502 (2019).
  * [26] L. D. Landau, L. M. Lifshitz, On the theory of the dispersion of magnetic permeabmlity in ferromagnetic bodies. Physik. Zeits. Sowjetunion 8, 153 (1935).
  * [27] T. L. Gilbert, A phenomenological theory of damping in ferromagnetic materials. IEEE Trans. Magn. 40, 3443 (2004).
  * [28] M. D. Stiles and J. Miltat, Spin-transfer torque and dynamics. Topics in Appl. Phys. 101, 225 (2006).
  * [29] M. Lakshmanan, The fascinating world of the Landau–Lifshitz–Gilbert equation: an overview. Phil. Trans. R. Soc. A 369, 1280 (2011).
  * [30] U. Atxitia et. al., Fundamentals and applications of the Landau–Lifshitz–Bloch equation, J. Phys. D: Appl. Phys. 50, 033003 (2017).
  * [31] Pavel A.Andreev, Exchange effects in Coulomb quantum plasmas: Dispersion of waves in 2D and 3D quantum plasmas, Ann. Phys. 350, 198 (2014).
  * [32] Pavel A. Andreev and L. S. Kuz’menkov, Dispersion properties of transverse waves in electrically polarized BECs, J. Phys. B: At. Mol. Opt. Phys. 47, 225301 (2014).
  * [33] Pavel A. Andreev and L. S. Kuz’menkov, Separated spin-up and spin-down evolution of degenerated electrons in two-dimensional systems: Dispersion of longitudinal collective excitations in plane and nanotube geometry, EPL 113, 17001 (2016).
  * [34] J. C. Slonczewski, Fluctuation mechanism for biquadratic exchange coupling in magnetic multilayers, Phys. Rev. Lett. 67, 3172 (1991).
  * [35] N. S. Fedorova, C. Ederer, N. A. Spaldin and A. Scaramucci, Biquadratic and ring exchange interactions in orthorhombic perovskite manganites. Phys. Rev. B 91, 165122 (2015).
  * [36] A. Kartsev, M. Augustin, R.F.L. Evans, et al., Biquadratic exchange interactions in two-dimensional magnets. npj Comput. Mater. 6, 150 (2020).
  * [37] P. A. Andreev and L. S. Kuz’menkov, Waves of spin current in magnetized dielectrics, Int. J. Mod. Phys. 29, 1550077 (2015).
  * [38] L. S. Kuz’menkov, S. G. Maksimov, V. V. Fedoseev, Microscopic quantum hydrodynamics of systems of fermions: II, Theor. Math. Phys. 126, 212 (2001).
  * [39] M.M. Bogdan, A.S. Kovalev, Exact multisoliton solution of one-dimensional Landau-Lifshitz equations for an anisotropic ferromagnet, JETP Lett. 31(8), 424 (1980).
  * [40] Sadamichi Maekawa et. al., Spin and spin current—From fundamentals to recent progress, J. Appl. Phys. 133, 020902 (2023).
  * [41] T. Nan, C. X. Quintela, J. Irwin, et al. Controlling spin current polarization through non-collinear antiferromagnetism, Nat. Commun. 11, 4671 (2020).
  * [42] Y. Kawaguchi, M. Ueda, Spinor Bose–Einstein condensates, Physics Reports, 520, Issue 5, 253 (2012).
  * [43] A. V. Ivanov, E. V. Zipunova, S. A. Khilkov, Equations of Correlational Magnetodynamics for Ferromagnetic Materials, JETP Lett. 115, 153 (2022).
  * [44] I. Cimrák, A Survey on the numerics and computations for the Landau-Lifshitz equation of micromagnetism, ARCO 15, 1 (2007).
  * [45] S. Brinker, M. dos Santos Dias and S. Lounis, The chiral biquadratic pair interaction, N. J. Phys. 21, 083015 (2019).
  * [46] F. Mila, F. C. Zhang, On the origin of biquadratic exchange in spin 1 chains, Eur. Phys. J. 16, 7 (2000).