\section{Introduction}

\paragraph{}
When formulating physics in a discrete geometric frame work, one is challenged by the absence
of a continuous diffeomorphism group. What is the analogue of a vector field on a
finite abstract simplicial complex $G$? Discrete theory approaches like
\cite{Forman1999,MarsdenDesbrun,RomeroSergeraert} have put forward some notions.
What about a combinatorial discrete frame work \cite{AmazingWorld}?

\paragraph{}
In order to get a continuum motion, one always can just look at paths in the
unitary group on the Hilbert space $l^2(G)$
like for example isospectral Lax deformations $D'=[B,D]$ of the exterior derivative $d$
defining $D=d+d^*$ \cite{IsospectralDirac,IsospectralDirac2}. But we would like to have a definition of
vector field which is formally identical to the classical case and which agrees with
the classical case if the differential complex comes from a Riemannian manifold.

\paragraph{}
As in the continuum, like on a Riemannian manifold \cite{AMR,FrankelGeometryPhysics}, we would like to
see vector fields related to $1$-forms, possibly moderated by a Riemannian metric,
The presence of an exterior derivative $d$
then produces potential fields $F=dV$ which then can be used to generate dynamics like
$x''=-\nabla V(x)$. In that case, the Riemannian structure allows to transfer the $1$-form $dV$
into a vector field $\nabla V$. In this note, we look at a
vector field notion in the discrete which works on any finite abstract simplicial complex $G$,
a finite set of non-empty sets $x$ invariant under the process of taking non-empty finite subsets.

\paragraph{}
We see that there is a finite dimensional Lie algebra of vector fields $X$
for which the Lie derivative $L_X$ defines a coordinate change
which commutes with exterior differentiation $d$. The coordinate changes allow to have
a basic general covariance principle. The motion can be extended to a Hamiltonian frame work so
that we could look also at analogues the Kepler problem, where a mass point moves
in a central field given by a potential associated to the geometry. Given a vector field $X$,
then the solution $\psi(t) = e^{i D_X t} \psi(0)$ of the second order and so
a Newton type equation $f_{tt} = -L_X f$ with $L_X=D_X^2$  and
$\psi(0)=f(0)-i D_X^{-1} f_t(0)$ resembles then the solution $e^{i D t} \psi(0)$ of the wave equation
$f_{tt} = -L f$ with $L=D^2$. It is the Cartan magic formula \cite{CartanRiemannianGeometry}, which is
now part of any differential topology book like \cite{AMR})
which produces the analogy between the wave equation with Hodge Laplacian $L=d d^* + d^* d = (d+d^*)^2=D^2$
and the Lie derivative $L_X = d i_X + i_X d = (d+i_X)^2 = D_X^2$. The {\bf Cartan formula}
can therefore be seen as a key to port notions of ordinary differential equations on manifolds to
discrete spaces like simplicial complexes.

\paragraph{}
The Cartan formula has been used in the past in discrete frame works (it appears in \cite{MMPDTKMD}).
Usually, in the Newton case as well as in the wave case, one does not write the dynamics
using complex coordinates. Here in the discrete, it is convenient as the equations become just
Schr\"odinger equations giving paths in the unitary group.
The wave equation case with Laplacian $L$ is the most
symmetric case, where the propagation of information happens in all directions or
(if the momentum $f_t$ is chosen accordingly) allows to force the propagation into a specific direction.
The analogue of ordinary differential equations are obtained when replacing $L$ with $L_X$,
in which  $X$ is a deterministic field, which assigns to a simplex, just one smaller
dimensional simplex.

\paragraph{}
The note draws also from insight gained in \cite{AnnieThesis,AnnieProject} and belongs
to the theme of looking for finite dimensional analogues of partial differential equations in the continuum.
The set-up for \cite{AnnieThesis} builds on work like \cite{Chapman} and is
an {\bf advection model} for a {\bf directed graph} $G$ is $u' = -Lu$, where
$L={\rm div} (V ({\rm grad}_0(u)))$ is the directed Laplacian on the graph. As in the usual (scalar) Laplacian
$L=d^*d$ for undirected graphs which leads to the {\bf heat equation} $u' = -Lu$, the
{\bf advection Laplacian} uses difference operators:
the modification of the gradient ${\rm grad}_0$ which is the maximum of ${\rm grad}$ and $0$.
The divergence ${\rm div}=d^*$ as well as ${\rm grad}=d$ are defined by the usual
exterior derivative $d$ on the graph.
The {\bf consensus model} is the situation, where the graph $G$ is replaced by its
{\bf reverse graph} $G^T$, where all directions are reversed. One can therefore focus
on advection. A central part of \cite{AnnieThesis} relates this to Markov chains. If $L=D-A$, then
$M=A D^{-1}$ is a {\bf left stochastic matrix}, a Markov operator, which maps
probability vectors to probability vectors. The kernel of $L$ is related to the
fixed points of $M$. Assume $L D^{-1} u=0$, then $(D-A) D^{-1} u=0$ and $u=A D^{-1} u = Mu$.
Perron-Frobenius allows to study the structure of the equilibria which are given by the
kernel of $L$. This concludes the diversion into advection.

\paragraph{}
What is the connection?
While the topic is related, we look here at differential equations on all differential
forms and not only on $0$-forms. Also, we don't yet really study the dynamics much and
just establish the linear algebra set-up showing that
there is an elementary way to define a Lie algebra of vector
fields in a discrete set-up. The affinity to the continuum is that the formalism is not only
similar but identical to the continuum. Whatever is done works both for Riemannian manifolds
or finite abstract simplicial complexes or more generally for a differential complex.
There are many open questions as seen at the end of this note:
integer-valued deterministic $X$ often produce integer eigenvalues of $L_X$ for example. We
would like to know when this is case appears.

\paragraph{}
An other angle emerged while teaching the multi-variable Taylor theorem in \cite{Math22-2018}.
Already the single variable Taylor theorem can be seen as the solution $f(t,x)=e^{D t} f(0)$ of the
{\bf transport equation} $f_t = D f$ with $D=d/dx$, a partial differential equation.
Because also $f(t,x)=f(x+t)$ solves this
partial differential equation, we have $f(x+t)=e^{Dt} f(0)$ which becomes so the Taylor theorem,
provided the initial function $f(0)$ is real analytic.
For a multivariate function we can replace $L=D^2$ with a Lie derivative $L_X=D_X^2$ and in the case of a
constant field $X=v$, a Taylor expansion $f(x+tv) = f(x) + df(x) tv + d^2f (x) (tv)^2/2 + ...$
(The multivariable Taylor theorem can be formulated conveniently using directional (Fr\'echet)
derivatives which is how textbooks like \cite{Edwards1973,BlatterIII1974} treats the subject in higher
dimensions, avoiding tensor calculus.)
As both $L$ and $L_X$ can be written as a square $L=D^2, L_X=D_X^2$, the analogy between the wave and Newton equation
has appeared. The frame work shows that a``diffeomorphism Lie group" exists in any geometric structure with
an exterior derivative; Taylor links the vector field Lie algebra with translation.

\section{From Cartan to d'Alembert}

\paragraph{}
Given a smooth compact manifold $M$ or a simplicial complex $G$ with {\bf exterior derivative}
$d: \Lambda^p \to \Lambda^{p+1}$, then every vector field $X$ defines an {\bf interior derivative}
$i_X: \Lambda^{p} \to \Lambda^{p-1}$. The {\bf Cartan magic formula} writes the {\bf Lie derivative}
$L_X$ as $L_X=d i_X + i_X d$. From the identities $d^2=0$ and $i_X^2=0$ follows that $L_X$
commutes with $d$. We know already from the continuum, that without the $d i_X$ part, the naive directional
derivative $i_X d$ alone would not work, as it would be coordinate dependent.
A $L_X$ commutes with $d$ it leads to a chain homotopy between the complexes before and after
the coordinate transformation. Like the {\bf Hodge Laplacian} $L=d d^* + d^* d$, {\bf we can write
$L_X$ as a square}: define $D_X = d+i_X$ and $D=d+d^*$. Then, $L_X=D_X^2$ and $L=D^2$.
The {\bf directional Dirac operator} $D_X$ has also an adjoint $D_X^*$ but it is different from $D_X$
in general. Despite the notation used, the directional Dirac operator is not the directional derivative
used in calculus. The operators $d,i_X,D_X$ and $L_X$ work on the linear space of all differential forms.

\paragraph{}
If $L=D^2$ is the Hodge Laplacian with Dirac operator $D=d+d^*$, then the {\bf wave equation}
is $f_{tt} = -L f$. The {\bf directional wave equation} is the formal analogue
$f_{tt} = -L_X f$. Written in the d'Alembert form, it is $(\partial_{tt}+L_X) f=0$.
As $L$ and $L_X$ are both squares of simpler operators $D$ and $D_X$, we can factor
$(\partial_t + i D_X) (\partial_t - i D_t) f =0$ or $(\partial_t + i D) (\partial_t - i D) f =0$.
The solutions $e^{\pm i D_X}$ which with Euler's formula $e^{i Dt} = \cos(Dt) + i \sin(Dt)$ leads
to the explicit solution $f(t)=\cos(Dt) f(0) + \sin(Dt) D^{-1} f_t(0)$, where $D^{-1}$ is the
pseudo-inverse is defined as $D^{-1} f_t^\perp(0)$ if $f_t^\perp(0)$ is in the orthogonal
complement of the kernel of $D$.

\paragraph{}
So far, the solutions of the equations were real-valued functions in the Hilbert space $H=l^2(G,\mathbb{R})$.
If $G$ is a finite abstract simplicial complex, the Hilbert space is finite dimensional, and the frame work
is part of linear algebra.
It is convenient to build the {\bf complex valued wave} $\psi(0) = f(0) - i D^{-1} f_t(0)$,
(where again $D^{-1}$ is the pseudo inverse) and get $\psi(t) = e^{i Dt} \psi(0)$.
The dynamics is now the solution to the {\bf Schr\"odinger equation} $i\psi_t = -D \psi$,
where $\psi(0)$ encodes the initial position $f(0)$ in its real part and the initial velocity $f_t(0)$
in its imaginary part. This works in the same way for $D_X$. In both cases, the wave equation for the real
wave is equivalent to the Schr\"odinger equation for a complex wave.
In the case of $0$-forms, we have $L=d^* d$ and $L_X=i_X d$ is the {\bf directional derivative} in the
direction $X$. We summarize:

{\it Given a geometric space with an exterior derivative $d$, the second order
real wave equation $f_{tt} = -L f$ is equivalent to the first order complex
Schr\"odinger equation $\psi' = i D \psi$, leading to a d'Alembert solution
$\psi(t) = e^{i D t} \psi(0)$ which can then be computed using a Taylor expansion and for
which the real part of $\psi$ gives the wave solution $f(t)$. }

\section{The Lie algebra}

\paragraph{}
Let us assume now that we are in a finite dimensional geometric space with an exterior derivative
$d: \Lambda^p \to \Lambda^{p+1}$ leading to a differential complex on a graded vector space
$\Lambda=\oplus_{p=0}^{{\rm dim}(G)} \Lambda^p$.
A {\bf vector field} is defined by a linear operator $i_X$ on $\Lambda$ which maps
$\Lambda^p$ to $\Lambda^{p-1}$ and has the property that $i_X^2 = 0$. We can define a Lie algebra
multiplication $Z=[X,Y]$ by first forming $L_X = i_X d + d i_X$ which is a map from
$\Lambda^p$ to $\Lambda^p$ and then defining $Z$ through the inner derivative
$$   i_Z = i_{[X,Y]} = [ L_X,i_Y ] = L_X i_Y - i_Y L_X \; . $$

\paragraph{}
The field $Z$ can be read of from $i_Z$. We also have the {\bf Lie algebra relation}
\begin{eqnarray*}
 L_Z  &=&  i_Z d + d i_Z = L_X i_Y d            - i_Y L_X d  + d L_X i_Y  - d i_Y          L_X \\
      &=&                  L_X (L_Y-d i_Y)      - i_Y L_X d  + L_X d i_Y  - (L_Y  - i_Y d) L_X \\
      &=&                  L_X L_Y - L_Y L_X   \; .
\end{eqnarray*}

\paragraph{}
Now, if $i_X i_Y = i_Y i_X = 0$, then
$L_X i_Y - i_Y L_X = i_X L_Y - L_Y i_X$ because
inter changing $X$ and $Y$ produces a change of sign of $L_Z$.
As $L_Z = i_Z d + d i_Z$, also $i_Z$ changes sign
meaning $Z$ changes sign. So, $i_Z = L_X i_Y - i_Y L_X = - (L_Y i_X - i_X L_Y)$.

\paragraph{}
These elementary matrix identities prove
the following proposition which applies to any so derived Lie algebra of fields $X$ with $i_X^2=0$.

\begin{propo}
Every vector field $X$ defines
an operator $D_X = d+i_X$ which has as a square a Lie derivative $L_X=D_X^2$. The set of $L_X$
define a Lie algebra with $L_{[X,Y]}=L_X L_Y - L_Y L_X$. If $X$ is in supported on odd forms,
then $i_X^2=0$ and $L_X=D_X^2$ holds in the entire Lie algebra.
\end{propo}

\paragraph{}
Even so $L_X$ is not self-adjoint in general, it plays the
role of a Laplacian. In the case when $L_X$ is not diagonalizable, we can not form the pseudo inverse
when writing the dynamics in the complex but we can assume that the initial velocity is
in the image of $D_X$. With this definition, also the adjoint operator $d^*$ defines
a vector field. We can still see $d^*=i_X$ for some field $X$. It belongs to the class of vector fields
$X$ for which the eigenvalues of $D_X$ are real.

\paragraph{}
We could look at a subclass of ``deterministic fields", which have the property that for a $p$-form $f$,
the $(p-1)$-form $i_X f$ is supported on a single sub-simplex $y$ of $x$, if $f$ is supported on a single
simplex $x$. This would correspond to classical vector fields close to the discrete Morse theory frame-work
\cite{Forman1999}.  If $X$ is supported on one-dimensional simplices, this is close to the
discrete Morse theory set-up. Unlike in the continuum, these ``deterministic fields"
are not invariant under addition, nor under the Lie algebra multiplication $[X,Y]$.
When taking the commutator of $L_X$ and $L_Y$ for such fields,
then the corresponding inner derivative $i_{[X,Y]}$ connects simplices which are not
directly connected.

\paragraph{}
If the simplicial complex is one-dimensional, or if $X$ is restrained to $p$-forms with odd $p$ (or even if we like)
then $i_Z=i_{[X,Y]}$ satisfies again $i_Z^2=0$ and the factorization $L_Z=(i_Z+d)^2= D_Z^2$ holds in the
entire Lie algebra. In general, the new interior derivative is only nilpotent, because
$i_Z^{1+{\rm dim}(G)} =0$. The Mathematica procedures below allow to support $X$ onto any subset of forms but
we mostly use the case when $X$ is supported on odd-dimensional forms.

\paragraph{}
Let us briefly look at the $1$-dimensional (single-variable) classical case $M=\mathbb{R}$, where
$L_X=i_X d$ is what we understand to be the usual derivative $d/dx$.
Technically, the exterior derivative $d$ produces from a $0$-form $f \in \Lambda^0$
a $1$-form $df dx \in \Lambda^1$ which is in a different vector space than $f$. But for the constant
vector field $X=1$, the combination $i_X d$ produces again an element in $\Lambda^0$.
Since $i_X^*=0$ on $0$-forms, we have $i_X d = i_X d + d i_X = L_X$. Now
$e^{L_X t} f = \sum_{n=0}^{\infty} (d/dx)^n f(x)X^n t^n/n!$ is by the {\bf Taylor formula} equal
to $f(x+t X)$, illustrating that the derivative $d/dx$ is the generator of the translation.
The flow $\phi_X f = f(x+t X)$ is a solution of a {\bf transport equation} and not the wave equation.
To get an analogue of the later, we need a second derivative in time and so a symplectic or complex
structure. It is first a bit puzzling to see the Lie derivative $L_X$ as a second order operator.
But the Cartan formula shows that also in one dimensions, $L_X = i_X d = (d + i_X)^2 = D_X^2$ is second
order. In calculus, we usually think of the derivative $d$ as a map on a space of {\bf scalar functions} and
not as a map from $0$-forms to $1$-forms. The inner derivative $i_X$ which brings us back to $0$-forms
is silently assumed in calculus. This identification of 0-forms and 1-forms can not be done in the
discrete because the dimension $v_1=|E|$ of $1$-forms is different from the dimension
$v_0=|V|$ of $0$-forms on a graph $G=(V,E)$. Still, the general frame work applies and the
wave equation $f_{tt} = - L_X f$ can be written as $(\partial_t - i D_X)  (\partial_t + i D_X)=0$
with $D_X = d + i_X$.

\paragraph{}
We would like to point out that the McKean-Singer symmetry
\cite{McKeanSinger} which holds for simplicial complexes and the
operator $L=D^2$ \cite{knillmckeansinger} remains valid also for $L_X=D_X^2$:

\begin{propo}[McKean-Singer symmetry]
The non-zero spectrum of $L_X$ on even differential forms $\Lambda^{even}$
is the same than the non-zero spectrum of $L_X$ on odd
differential forms $\Lambda^{odd}$.
\end{propo}

\begin{proof}
The proof is the same than in the continuum \cite{Cycon}
or used in \cite{knillmckeansinger}: the operator $D_X$ which
exchanges $\Lambda^{even}$ and $\Lambda^{odd}$ gives
a translation between the eigenvectors belonging to non-zero
eigenvalues. The discrepancy between the kernels on odd and
even forms is by definition the Euler characteristic.
\end{proof}

\paragraph{}
Also the nonlinear {\bf Lax type deformation}
\cite{IsospectralDirac,IsospectralDirac2}
of the Dirac operator generalizes from $D$ to $D_X$. These Lax equations are
$$   D_X' = [B_X,D_X]  $$
where $D_X=d+i_X + b_X, B_X=d_X-d_X^*+i b_X$.
Unlike for $i_X=d^*$, where $L=L_X$ is the Hodge Laplacian, the deformation
is now not isospectral in general and therefor not expected to be integrable.
The expansion rate in different part of space or differential forms happens differently.
Still, these systems remain interesting non-linear differential equations
and the corresponding $d(t)$ still satisfies $d(t)^2=0$ producing an exterior
derivative after deformation. As in the case of the wave equation, the deformed
$d(t)$ keeps the same cohomology.

\section{Physics}

\paragraph{}
Any mathematical theory with some quantum gravitational ambitions should be able to be powerful enough
to solve the Kepler problem effectively in any scale: in the large, it should lead to the
classical Kepler problem, in the very large to relativistic motion in a Schwarzschild metric and
in the very small to the quantum dynamics in the Hydrogen atom. No current theory passes this
Kepler test: no theory can yet describe a point in the influence
of a central field classically, relativistically and quantum mechanically, not just in principle
or a perturbative patch work  but in an elegant manner, leading to quantitative results
which match experiments in all three scales. It should be able to describe
the motion of satellites or planets, also relativistically, predict the emission patters of
gravitational waves emitted by a binary system or the structure of the Mendeleev table in the small.

\paragraph{}
The general covariance principle in physics states that physical laws
are independent of the coordinate system. This means that the laws should not
only be invariant under a finite dimensional symmetry group like Euclidean or Lorentz symmetry
but they should be invariant under the diffeomorphism group of the manifold.
Additionally, in case of fibre bundles, additional gauge symmetries might apply, but this
is also part of the general covariance principle. An example are the Maxwell equations $dF=0,
d^*F=j$ leading in the Coulomb gauge $d^* A=0$ to the Poisson equation $L A=j$.
If we move into a new coordinate system, then the transported equations look the same.
An other example are the Einstein field equations $G = e T$, relating the
geometric Einstein tensor with the energy tensor $T$ using a proportionality factor $e$, the
Einstein constant. The covariance there is there the statement that $G$ and
$T$ are tensors. How can one port the general covariance principle to the discrete?
A naive request would be to look at laws only which are invariant after applying a deformation
through a vector field. Since $L_X$ and $d$ commute, any law which only involves the
exterior derivative does this. Examples are the wave, the heat or the Schr\"odinger equation.

\section{Examples}

\paragraph{}
The simplest case with a non-trivial Lie algebra is when
$G=\{ \{1\},\{2\},\{1,2\}\}$ is the Whitney complex of the complete graph $K_2$.
In that case,
$$ d   = \left[ \begin{array}{ccc} 0 & 0 & 0 \\ 0 & 0 & 0 \\ -1 & 1 & 0 \\ \end{array} \right],
   d^* = \left[ \begin{array}{ccc} 0 & 0 & -1 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \\ \end{array} \right] \; . $$
The general inner derivative (vector field) has the form
$$ i_X = \left[ \begin{array}{ccc} 0 & 0 & a \\ 0 & 0 & b \\ 0 & 0 & 0 \\ \end{array} \right] \; . $$
The general operator $D_X=d + i_X$ and $L_X=D_X^2= d i_X + i_X d$ then is
$$ D_X = \left[ \begin{array}{ccc} 0 & 0 & a \\ 0 & 0 & b \\ -1 & 1 & 0 \\ \end{array} \right],
   L_X = \left[ \begin{array}{ccc} -a & a & 0 \\ -b & b & 0 \\ 0 & 0 & b-a \\ \end{array} \right]  \; . $$
The eigenvalues of $L_X$ are $\{ 0, b-a, b-a \}$. Given an other vector field
$$ i_Y=\left[ \begin{array}{ccc} 0 & 0 & u \\ 0 & 0 & v \\ 0 & 0 & 0 \\ \end{array} \right]  \; , $$
one can form $i_Z  = L_X i_Y-i_Y L_X = i_X L_Y - L_Y i_X$ which is
$$ i_Z=\left[ \begin{array}{ccc} 0 & 0 & a v-b u \\ 0 & 0 & a v-b u \\ 0 & 0 & 0 \\ \end{array} \right] \; , $$
leading to
$$ D_Z=\left[ \begin{array}{ccc} 0 & 0 & a v-b u \\ 0 & 0 & a v-b u \\ -1 & 1 & 0 \\ \end{array} \right],
   L_Z=\left[ \begin{array}{ccc} b u-a v & a v-b u & 0 \\ b u-a v & a v-b u & 0 \\ 0 & 0 & 0 \\ \end{array} \right] \; .$$
This satisfies $L_Z = L_X L_Y - L_Y L_X$. The eigenvalues of $L_Z$ are all zero. Indeed $L_Z^2=0$.

\begin{comment}
d={{0,0,0}, {0,0,0},{-1,1,0}};
dt={{0,0,-1},{0,0,1}, {0,0,0}};
iX={{0,0,a},{0,0,b},{0,0,0}};
iY={{0,0,u},{0,0,v},{0,0,0}};
DX=iX+d; LX = DX.DX; DY=iY+d; LY = DY.DY;
iZ1= LX.iY-iY.LX; iZ2= iX.LY-LY.iX;
iZ=iZ1; LZ=iZ.d+d.iZ; DZ=iZ+d;
DZ.DZ==LZ
Simplify[LX.LY-LY.LX==LZ]
\end{comment}

\paragraph{}
Here is the general case if $G=\{ \{1\},\{2\}, \{3\},\{1,2\},\{1,3\} \}$ is the Whitney
complex of a linear graph of length $2$.
$$ d = \left[    \begin{array}{ccccc}
                   0 & 0 & 0 & 0 & 0 \\
                   0 & 0 & 0 & 0 & 0 \\
                   0 & 0 & 0 & 0 & 0 \\
                   -1 & 1 & 0 & 0 & 0 \\
                   -1 & 0 & 1 & 0 & 0 \\
                  \end{array} \right],
   d^* = \left[  \begin{array}{ccccc}
                   0 & 0 & 0 & -1 & -1 \\
                   0 & 0 & 0 & 1 & 0 \\
                   0 & 0 & 0 & 0 & 1 \\
                   0 & 0 & 0 & 0 & 0 \\
                   0 & 0 & 0 & 0 & 0 \\
                  \end{array} \right] \; . $$
With a vector fields
$$   i_X = \left[
                  \begin{array}{ccccc}
                   0 & 0 & 0 & a_1 & a_2 \\
                   0 & 0 & 0 & a_3 & 0   \\
                   0 & 0 & 0 &  0  & a_6 \\
                   0 & 0 & 0 & 0 & 0 \\
                   0 & 0 & 0 & 0 & 0 \\
                  \end{array} \right],
     D_X = \left[
                  \begin{array}{ccccc}
                   0 & 0 & 0 & a_1 & a_2 \\
                   0 & 0 & 0 & a_3 &  0  \\
                   0 & 0 & 0 &  0  & a_6 \\
                  -1 & 1 & 0 & 0 & 0 \\
                  -1 & 0 & 1 & 0 & 0 \\
                  \end{array} \right] \; . $$
This leads to
$$ L_X = \left[
                  \begin{array}{ccccc}
                   -a_1-a_2 & a_1 & a_2 & 0 & 0 \\
                   -a_3     & a_3 &  0  & 0 & 0 \\
                       -a_6 &   0 & a_6 & 0 & 0 \\
                   0 & 0 & 0 & a_3-a_1 &    -a_2 \\
                   0 & 0 & 0 &    -a_1 & a_6-a_2 \\
                  \end{array} \right] \; . $$
Given an other vector field
$$   i_Y = \left[
                  \begin{array}{ccccc}
                   0 & 0 & 0 & b_1 & b_2 \\
                   0 & 0 & 0 & b_3 &  0   \\
                   0 & 0 & 0 & 0   & b_6 \\
                   0 & 0 & 0 & 0 & 0 \\
                   0 & 0 & 0 & 0 & 0 \\
                  \end{array} \right],
     D_Y = \left[
                  \begin{array}{ccccc}
                   0 & 0 & 0 & b_1 & b_2 \\
                   0 & 0 & 0 & b_3 &  0  \\
                   0 & 0 & 0 &  0  & b_6 \\
                  -1 & 1 & 0 & 0 & 0 \\
                  -1 & 0 & 1 & 0 & 0 \\
                  \end{array} \right] \; , $$
we get
\begin{tiny}
$$ i_Z = \left[ \begin{array}{ccccc}
 0 & 0 & 0 & -a_2 b_1-a_3 b_1+a_1 (b_2+b_3) & a_2 (b_1+b_6)-(a_1+a_6) b_2 \\
                   0 & 0 & 0 & a_1 b_3-a_3  b_1 & a_2 b_3-a_3 b_2 \\
                   0 & 0 & 0 & a_1 b_6-a_6 b_1 & a_2 b_6-a_6 b_2 \\
                   0 & 0 & 0 & 0 & 0 \\
                   0 & 0 & 0 & 0 & 0 \\
                  \end{array} \right] \; . $$
\end{tiny}
\begin{comment}
and
\begin{tiny}
$$ L_Z = \left[ \begin{array}{ccccc}
                   -(-a_1-a_2) b_1-a_2 b_1+(a_3-a_1) b_1-a_1
                     b_2-(-a_1-a_2) b_2+(a_6-a_2) b_2-a_1
                     b_3-a_2 b_6 & (-a_1-a_2) b_1-(a_3-a_1)
                     b_1+a_1 b_2+a_1 b_3 & a_2 b_1+(-a_1-a_2)
                     b_2-(a_6-a_2) b_2+a_2 b_6 & 0 & 0 \\
                   a_3 b_1+a_3 b_2-a_2 b_3-a_3
                     b_3+(a_3-a_1) b_3 & -a_3 b_1+a_3
                     b_3-(a_3-a_1) b_3 & a_2 b_3-a_3 b_2 & 0 & 0 \\
                   a_6 b_1+a_6 b_2-a_1 b_6-a_6
                     b_6+(a_6-a_2) b_6 & a_1 b_6-a_6 b_1 &
                     -a_6 b_2+a_6 b_6-(a_6-a_2) b_6 & 0 & 0 \\
                   0 & 0 & 0 & -(-a_1-a_2) b_1-a_3 b_1+(a_3-a_1)
                     b_1-a_1 b_2-a_1 b_3+a_3 b_3-(a_3-a_1)
                     b_3 & -a_2 b_1-(-a_1-a_2) b_2-a_3
                     b_2+(a_6-a_2) b_2+a_2 b_3-a_2 b_6 \\
                   0 & 0 & 0 & -(-a_1-a_2) b_1+(a_3-a_1) b_1-a_6
                     b_1-a_1 b_2-a_1 b_3+a_1 b_6 & -a_2
                     b_1-(-a_1-a_2) b_2-a_6 b_2+(a_6-a_2)
                     b_2-a_2 b_6+a_6 b_6-(a_6-a_2) b_6 \\
                  \end{array} \right] \; . $$
\end{tiny}
\end{comment}
The eigenvalues of $L_X$ contain $0$ as well as the following two eigenvalues,
each with multiplicity two: \\
$\left(\pm \sqrt{a_1^2+2 a_1 (a_2-a_3+a_6)+(a_2+a_3-a_6)^2}-a_1-a_2+a_3+a_6\right)/2$.
We see that real eigenvalues are quite common but that imaginary eigenvalues of
$L_X$ can occur in the $4$-dimensional space of vector fields. The operator $L_Z$
does in general not have zero eigenvalues. They can even become complex.
We also see that the commutator $L_Z$ now tunnels between places which were not
directly connected in $G$.

\begin{comment}
d ={{0,0,0,0,0},  {0,0,0,0,0},{0,0,0,0,0},{-1,1,0,0,0},{-1,0,1,0,0}};
dt={{0,0,0,-1,-1},{0,0,0,1,0},{0,0,0,0,1},{0,0,0,0,0}, {0,0,0,0,0}};
iX={{0,0,0,a1,a2},{0,0,0,a3,a4},{0,0,0,a5,a6},{0,0,0,0,0},{0,0,0,0,0}};
iY={{0,0,0,b1,b2},{0,0,0,b3,b4},{0,0,0,b5,b6},{0,0,0,0,0},{0,0,0,0,0}};
a4=0; a5=0; b4=0; b5=0;

DX=iX+d;LX=DX.DX;DY=iY+d;LY=DY.DY;
iZ1=LX.iY-iY.LX;iZ2=iX.LY-LY.iX;
Simplify[iZ1==iZ2]
iZ=iZ1;LZ=iZ.d+d.iZ;
DZ=iZ+d; Simplify[DZ.DZ==LZ]
Simplify[LX.LY-LY.LX==LZ]
a1=Random[Integer,1]; a2=Random[Integer,1]; a3=Random[Integer,1]; a6=Random[Integer,1];
b1=Random[Integer,1]; b2=Random[Integer,1]; b3=Random[Integer,1]; b6=Random[Integer,1];
Eigenvalues[LX]
Eigenvalues[LY]
Eigenvalues[LZ]
\end{comment}

\paragraph{}
For the complete complex $G=\{\{1\},\{2\},\{3\},\{1,2\},\{2,3\},\{1,3\},\{1,2,3\}\}$,
we can look at the general
$$ D_X = d + iX = \left[ \begin{array}{ccccccc}
  0 & 0 & 0 & a & b & 0 & 0 \\
  0 & 0 & 0 & c & 0 & d & 0 \\
  0 & 0 & 0 & 0 & e & f & 0 \\
 -1 & 1 & 0 & 0 & 0 & 0 & g \\
 -1 & 0 & 1 & 0 & 0 & 0 & h \\
  0 & -1 & 1 & 0 & 0 & 0 & i \\
  0 & 0 & 0 & 1 & -1 & 1 & 0 \\
                 \end{array} \right] \; . $$
Now, given two general $D_X,D_Y$, we have $L_X=D_X^2, L_Y=D_Y^2$.
If $g=h=i=0$ then, $L_X i_Y-i_Y L_X = i_X L_Y-L_Y i_X$
and $i_Z=L_X i_Y-i_Y L_X$ has the property that $i_Z^2=0$.

\begin{comment}
d={{0,0,0,0,0,0,0},{0,0,0,0,0,0,0},{0,0,0,0,0,0,0},{-1,1,0,0,0,0,0},
{-1,0,1,0,0,0,0},{0,-1,1,0,0,0,0},{0,0,0,1,-1,1,0}};dt=Transpose[d];
iX={{0,0,0,a1,b1,0,0},{0,0,0,c1,0,d1,0},{0,0,0,0,e1,f1,0},{0,0,0,0,0,0,g1},
    {0,0,0,0,0,0,h1},{0,0,0,0,0,0,i1}, {0,0,0,0,0,0,0}};
iY={{0,0,0,a2,b2,0,0},{0,0,0,c2,0,d2,0},{0,0,0,0,e2,f2,0},{0,0,0,0,0,0,g2},
    {0,0,0,0,0,0,h2},{0,0,0,0,0,0,i2}, {0,0,0,0,0,0,0}};
g1=0; h1=0; i1=0; g2=0; h2=0; i2=0;
DX=iX+d;LX=DX.DX;DY=iY+d;LY=DY.DY;
iZ1=LX.iY-iY.LX;iZ2=iX.LY-LY.iX;
Simplify[iZ1==iZ2]
iZ=iZ1;LZ=iZ.d+d.iZ;
DZ=iZ+d; Simplify[DZ.DZ==LZ]
Simplify[LX.LY-LY.LX==LZ]
R:=Random[]; a1=R; b1=R; c1=R; d1=R; e1=R; f1=R; g1=0; h1=0; i1=0;
R:=Random[]; a2=R; b2=R; c2=R; d2=R; e2=R; f2=R; g2=0; h2=0; i2=0;
Chop[Eigenvalues[LX]]
Chop[Eigenvalues[LY]]
Chop[Eigenvalues[LZ]]
\end{comment}

\paragraph{}
Let $G=\{\{1\},\{2\},\{3\},\{4\},\{1,2\},\{2,3\},\{3,4\},\{4,1\}\}$
be the complex of the cyclic graph $C_4$. The exterior derivative $d$ and an example of an
interior derivative $i_X$ are:

\begin{tiny}
$$          d =  \left[
                 \begin{array}{cccccccc}
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 -1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & -1 & 1 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & -1 & 1 & 0 & 0 & 0 & 0 \\
 1 & 0 & 0 & -1 & 0 & 0 & 0 & 0 \\
                 \end{array}
                 \right],
           i_X = \left[
                 \begin{array}{cccccccc}
 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
                 \end{array}
                 \right] \; . $$
This leads to the Dirac operator $D=d+d^*$ and the directional Dirac operator
$D_X=d+i_X$
$$          D   = \left[
                  \begin{array}{cccccccc}
 0 & 0 & 0 & 0 & -1 & 0 & 0 & 1 \\
 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 1 & -1 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1 \\
 -1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & -1 & 1 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & -1 & 1 & 0 & 0 & 0 & 0 \\
 1 & 0 & 0 & -1 & 0 & 0 & 0 & 0 \\
                  \end{array}
                  \right],
            D_X = \left[
                  \begin{array}{cccccccc}
 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\
 -1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & -1 & 1 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & -1 & 1 & 0 & 0 & 0 & 0 \\
 1 & 0 & 0 & -1 & 0 & 0 & 0 & 0 \\
                  \end{array}
                  \right]  \; . $$
The Hodge Laplacian $L=D^2$ and Lie derivative $L_X=D_X^2$ are
$$  L = \left[
                  \begin{array}{cccccccc}
 2 & -1 & 0 & -1 & 0 & 0 & 0 & 0 \\
 -1 & 2 & -1 & 0 & 0 & 0 & 0 & 0 \\
 0 & -1 & 2 & -1 & 0 & 0 & 0 & 0 \\
 -1 & 0 & -1 & 2 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 2 & -1 & 0 & -1 \\
 0 & 0 & 0 & 0 & -1 & 2 & -1 & 0 \\
 0 & 0 & 0 & 0 & 0 & -1 & 2 & -1 \\
 0 & 0 & 0 & 0 & -1 & 0 & -1 & 2 \\
                  \end{array}
                  \right],
    L_X =  \left[
                  \begin{array}{cccccccc}
 1 & -1 & 0 & 0 & 0 & 0 & 0 & 0 \\
 -1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & -1 & 1 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & -1 & 1 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 2 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & -1 & 1 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & -1 & 1 & 0 \\
 0 & 0 & 0 & 0 & -1 & 0 & -1 & 0 \\
                  \end{array}
                  \right] \; . $$
\end{tiny}
The eigenvalues of $D$ are $\{-2, 2, -\sqrt{2}, -\sqrt{2}, \sqrt{2}, \sqrt{2}, 0, 0\}$, the eigenvalues
of $D_X$ are $\{ -\sqrt{2}, -\sqrt{2}, -1, -1, 1, 1, 0, 0 \}$.
The eigenvalues of $L_X$ are $\{ 4, 4, 2, 2, 2, 2, 0, 0 \}$.

\begin{figure}[!htpb]
\scalebox{0.64}{\includegraphics{figures/cyclegraph_spectrum.pdf}}
\scalebox{0.64}{\includegraphics{figures/wheelgraph_spectrum.pdf}}
\label{example}
\caption{
The figures show the spectrum of $D$ and $D_X$ for a cycle graph spectrum and a wheel graph spectrum.
In both cases, we start with $i_X=d^*$, then take away each entry $1$ or
$-1$ with probability $1/2$. $D_X$ still has the $\sigma(D_X)=-\sigma(D_X)$ symmetry
from the Dirac operator $D$ but the energies of the particles are smaller as it is
more difficult to travel.
}
\end{figure}

\section{Mathematica procedures}

\paragraph{}
Here is the code which computes the Dirac operator $D$ and the
vector field analogue $D_X$ for any simplicial complex $G$. The
code can be copy pasted when accessing the LaTeX source of this document
on the ArXiv. The first part computes the
matrices given in the above example. For the vector field, we chose
for $i_X$ just to take the first non-zero entry of $d^*$:

\begin{tiny}
\lstset{language=Mathematica} \lstset{frameround=fttt}
\begin{lstlisting}[frame=single]
G={{1},{2},{3},{4},{1,2},{2,3},{3,4},{4,1}};
n=Length[G];Dim=Map[Length,G]-1;f=Delete[BinCounts[Dim],1];
Orient[a_,b_]:=Module[{z,c,k=Length[a],l=Length[b]},
  If[SubsetQ[a,b] && (k==l+1),z=Complement[a,b][[1]];
  c=Prepend[b,z]; Signature[a]*Signature[c],0]];
d=Table[0,{n},{n}];d=Table[Orient[G[[i]],G[[j]]],{i,n},{j,n}];
dt=Transpose[d]; DD=d+dt; LL=DD.DD;
HX[x_]:=Block[{u=Flatten[Position[Abs[x],1]]},If[u=={},0,First[u]]];
iX=Table[0,{n},{n}];
Do[l=HX[dt[[k]]];If[l>0,iX[[k,l]]=dt[[k,l]]],{k,f[[1]]}];
DX=iX+d; LX = DX.DX;
\end{lstlisting}
\end{tiny}

\paragraph{}
Here is the code to generate $D_X$ and $L_X$ for a generate random
finite abstract simplicial complex $G$:

\begin{tiny}
\lstset{language=Mathematica} \lstset{frameround=fttt}
\begin{lstlisting}[frame=single]
(* Generate a random simplicial complex                               *)
Generate[A_]:=Delete[Union[Sort[Flatten[Map[Subsets,A],1]]],1]
R[n_,m_]:=Module[{A={},X=Range[n],k},Do[k:=1+Random[Integer,n-1];
  A=Append[A,Union[RandomChoice[X,k]]],{m}];Generate[A]];
G=Sort[R[10,20]];

(* Computation of exterior derivative                                 *)
n=Length[G];Dim=Map[Length,G]-1;f=Delete[BinCounts[Dim],1];
Orient[a_,b_]:=Module[{z,c,k=Length[a],l=Length[b]},
  If[SubsetQ[a,b] && (k==l+1),z=Complement[a,b][[1]];
  c=Prepend[b,z]; Signature[a]*Signature[c],0]];
d=Table[0,{n},{n}];d=Table[Orient[G[[i]],G[[j]]],{i,n},{j,n}];
dt=Transpose[d]; DD=d+dt; LL=DD.DD;

(* Build interior derivatives iX  and iY                              *)
UseIntegers=False;
e={}; Do[If[Length[G[[k]]]==2,e=Append[e,k]],{k,n}];
BuildField[P_]:=Module[{X,ee,iX=Table[0,{n},{n}]},
X=Table[If[UseIntegers,Random[Integer,1],Random[]],{l,Length[e]}];
Do[ee=G[[e[[l]]]]; Do[If[SubsetQ[G[[k]],ee],
  m=Position[G,Sort[Complement[G[[k]],Delete[ee,2]]]][[1,1]];
  iX[[m,k]]= If[MemberQ[P,Length[G[[m]]]],X[[l]],0]*
             Orient[G[[k]],G[[m]]]],{k,n}],{l,Length[e]}]; iX];

(* Build Laplacians LX,LY,LZ, plot spectrum of D and DX and matrices  *)
iX=BuildField[{1,3,5,7,9}]; iY=BuildField[{1,3,5,7,9}];
DX=iX+d; LX=DX.DX; DY=iY+d; LY=DY.DY; iZ1=LX.iY-iY.LX; iZ2=iX.LY-LY.iX;
Print[iZ1==iZ2]; iZ=iZ1; LZ=Chop[iZ.d+d.iZ]; DZ=iZ+d;
dx="\!\(\*SubscriptBox[\(D\), \(X\)]\)";
lx="\!\(\*SubscriptBox[\(L\), \(X\)]\)"; pl=PlotLabel;
GraphicsGrid[{{MatrixPlot[DX,pl->dx], MatrixPlot[LX,pl->lx]},
              {MatrixPlot[DD,pl->"D"], MatrixPlot[LL,pl->"L"]}}];
u1 = Sort[Chop[Eigenvalues[1.0 DX]]]; u2 = Sort[Eigenvalues[1.0 DD]];
u1=N[Round[u1*10^6]/10^6];            (* clear tiny imaginary parts  *)
S=ListPlot[{u1, u2},Joined ->True,PlotLegends ->Placed[{dx,"D"},Below],
 PlotRange -> All, PlotStyle -> {Red, Blue},
 PlotLabel -> "Spectrum of \!\(\*SubscriptBox[\(D\), \(X\)]\) and D"];

(* Compute Betti numbers, compare bosonic and fermionic part          *)
chi=Sum[-f[[k]](-1)^k,{k,Length[f]}]; f=Prepend[f,0]; m=Length[f]-1;
U=Table[v=f[[k+1]];
   Table[u=Sum[f[[l]],{l,k}];LL[[u+i,u+j]],{i,v},{j,v}],{k,m}];
Cohomology = Map[NullSpace, U]; Betti = Map[Length, Cohomology]
chi1=Sum[-Betti[[k]](-1)^k,{k,Length[Betti]}];
EV=Map[Eigenvalues,U];
EVFermi=Table[EV[[2k]],{k,Floor[Length[EV]/2]}];
EVBoson=Table[EV[[2k-1]],{k,Floor[(Length[EV]+1)/2]}];
extract[u_]:=Module[{v=Flatten[u],w={}},
   Do[If[Abs[v[[k]]]>10^(-8),w=Append[w,v[[k]]]],{k,Length[v]}];Sort[w]];
extract[EVFermi]==extract[EVBoson]

(* Now the same for LX                                               *)
U=Table[v=f[[k+1]];
   Table[u=Sum[f[[l]],{l,k}];LX[[u+i,u+j]],{i,v},{j,v}],{k,m}];
EV=Map[Eigenvalues,U];
EVFermi=Table[EV[[2k]],{k,Floor[Length[EV]/2]}];
EVBoson=Table[EV[[2k-1]],{k,Floor[(Length[EV]+1)/2]}];
extract[EVFermi]==extract[EVBoson]
{EVFermi,EVBoson}
Total[Abs[N[extract[EVFermi] - extract[EVBoson]]]]

Cohomology = Map[NullSpace, U]; Betti = Map[Length, Cohomology];
chi2=Sum[-Betti[[k]](-1)^k,{k,Length[Betti]}];
{chi,chi1,chi2}
\end{lstlisting} \end{tiny}

\begin{figure}[!htpb]
\scalebox{1.0}{\includegraphics{figures/matrices.pdf}}
\label{example}
\caption{
The matrices $D_X,L_X,D,L$ in the case of a random complex.
This was produced with the code above.
}
\end{figure}

\paragraph{}
And finally, here is the self-contained procedure which does the isospectral deformation of
the exterior derivative by deforming $D_X'=[B_X,D_X]$. In the case
$i_X=d^*$, this is the standard {\bf Lax isospectral deformation} we have seen
before \cite{IsospectralDirac,IsospectralDirac2}.
In an other extreme case, when $X=0$, then $d$ is not deformed at all.
We still have the inflationary initial decay of $d$ typical for that
type of integrable dynamical system. The decay of $d$ means by
the Connes formula that there is an expansion of space because if
the derivative operator become small, then the distances grow.

\begin{tiny}
\lstset{language=Mathematica} \lstset{frameround=fttt}
\begin{lstlisting}[frame=single]
Generate[A_]:=Delete[Union[Sort[Flatten[Map[Subsets,A],1]]],1]
R[n_,m_]:=Module[{A={},X=Range[n],k},Do[k:=1+Random[Integer,n-1];
  A=Append[A,Union[RandomChoice[X,k]]],{m}];Generate[A]];G=Sort[R[5,8]];
n=Length[G]; fv=Delete[BinCounts[Map[Length,G]],1];
cn=Length[fv];br={0};Do[br=Append[br,Last[br]+fv[[k]]],{k,cn}];

Orient[a_,b_]:=Module[{z,c,k=Length[a],l=Length[b]},
  If[SubsetQ[a,b] && (k==l+1),z=Complement[a,b][[1]];
  c=Prepend[b,z]; Signature[a]*Signature[c],0]];
d=Table[0,{n},{n}];d=Table[Orient[G[[i]],G[[j]]],{i,n},{j,n}];
dt=Transpose[d]; DD=d+dt; LL=DD.DD;

UseIntegers=False;e={}; Do[If[Length[G[[k]]]==2,e=Append[e,k]],{k,n}];
BuildField[P_]:=Module[{X,ee,iX=Table[0,{n},{n}]},
X=Table[If[UseIntegers,Random[Integer,1],Random[]],{l,Length[e]}];
Do[ee=G[[e[[l]]]]; Do[If[SubsetQ[G[[k]],ee],
  m=Position[G,Sort[Complement[G[[k]],Delete[ee,2]]]][[1,1]];
  iX[[m,k]]= If[MemberQ[P,Length[G[[m]]]],X[[l]],0]*
             Orient[G[[k]],G[[m]]]],{k,n}],{l,Length[e]}]; iX];
iX=BuildField[{1,3,5,7,9}]; iY=BuildField[{1,3,5,7,9}];
DX=iX+d; LX=DX.DX; DY=iY+d; LY=DY.DY; iZ=LX.iY-iY.LX; DZ=iZ+d;LZ=DZ.DZ;

T[A_]:=Module[{n=Length[A]},Table[If[i<=j,0,A[[i,j]]],{i,n},{j,n}]];
UT[{DD_,br_}]:=Module[{D1=T[DD]},     (* Lower triangular block  *)
Do[Do[Do[D1[[br[[k]]+i,br[[k]]+j]]=0,{i,br[[k+1]]-br[[k]]}],
 {j,br[[k+1]]-br[[k]]}],{k,Length[br]-1}];D1];
RuKu[f_,x_,s_]:=Module[{a,b,c,u,v,w,q},u=s*f[x]; (* Runge Kutta  *)
a=x+u/2;v=s*f[a];b=x+v/2;w=s*f[b];c=x+w; q=s*f[c];x+(u+2v+2w+q)/6];

DD=DX; d0=UT[{DD,br}]; e0=Conjugate[Transpose[d0]];
M=1000; delta=2/M; u={};  (* Deformation with Runge Kutta      *)
Do[d=UT[{DD,br}];e=Conjugate[Transpose[d]];
BB=d-e; CC=d+e; MM=CC.CC; b=DD-CC; VV=b.b;
B=BB+1.0*I*b; f[x_]:=B.x-x.B; DD=RuKu[f,1.0 DD,delta];
u=Append[u,Total[Abs[Flatten[Chop[d]]]]],{m,M}];
DDX=DD; LLX=DDX.DDX;

{Total[Abs[Flatten[Chop[d.d]]]], Total[Abs[Flatten[Chop[e.e]]]]}

F[x_]:=If[x==0,0,-Log[Abs[x]]];  (* Plot the size of d      *)
v=M*Table[F[u[[k+1]]]-F[u[[k]]],{k,Length[u]-1}];
ListPlot[v,PlotRange->All]
\end{lstlisting} \end{tiny}

\begin{figure}[!htpb]
\scalebox{1.0}{\includegraphics{figures/matrices2.pdf}}
\label{example}
\caption{
The matrices $D_Z,L_Z,D,L$ in the case of a random complex
with $f$-vector $(10, 45, 120, 192, 165, 73, 15, 1)$.
Also this was produced with the code above, where $Z=[X,Y]$
is the commutator of two random vector fields.
}
\end{figure}

\section{Questions}

\paragraph{}
The operators $D_X$ and $L_X$ are not symmetric in general so that complex eigenvalues can appear.
In that case, the solution $\psi(t) = e^{i D_X t} \psi(0)$ can grow exponentially.
$D_X$ still often has real eigenvalues, leading to quasi-periodic solutions as the orbits $e^{i D_X t} \psi(0)$ form
a subgroup of a finite dimensional torus, if the graph is finite.  Actually, if $i_X(k,l) \neq 0$ only for
one $l$, then we implement a deterministic vector field. In that case the eigenvalues of $L_X$ often
non-negative integers taking values in $\{0,1, \dots, {\rm dim}(G)+1 \}$. We would like to understand
the spectrum.

\paragraph{}
For $D=i_X + i_X^* + d + d^*$ we have $D^2 = L_X + L_X^* + L$. Now, if we average that over all
possible vector fields using a measure which is homogeneous, we expect the $L_X$ to
average out and get the wave equation governed by $L$. Can one make this more precise and see
the wave equation $f_{tt} = - L$ as an average of deterministic flows $f_{tt} = - L_X$?

\paragraph{}
We often integer eigenvalues of $L_X$ if $i_{X}$ has integer values. In small dimensional examples,
we can compute general formulas for the eigenvalues but integer eigenvalues also often appear for
large random simplicial complexes. Under which conditions does $L_X$ have integer eigenvalues?

\paragraph{}
The eigenvalues of $L_X$ are most of the time real if the entries of $I_X$ are non-negative
multiplies of $d^*$. They can become imaginary in general, if the signs are changed.
Can we find conditions which assures a real spectrum?

\bibliographystyle{plain}


\section{Introduction}
\label{sec-introduction}

The {\emph{lifted cover inequality}} (LCI) is \rev{a} well-known cutting plane for mixed integer programs.
\highlight{Given \rev{the so-called cover inequality}, in order to obtain an LCI, we may use
the {\textit{lifting} technique}. Using different lifting procedures, several types of LCIs
have been studied in the literatures, see \cite{Balas1975,Wolsey1975,Hammer1975,Balas1978a,Hartvigsen1992,Gu2000,Atamturk2005a,Letchford2019}. In this paper, we are
concentrated on the {\textit{sequential}} LCI, that is,
the variables are sequentially lifted one at a time.}
The sequential LCI was first studied in \highlight{\cite{Padberg1975a,Wolsey1975}}. Its effectiveness on using as \rev{a} cutting plane was demonstrated in \cite{Crowder1983}, see also \cite{VanRoy1987,Hoffman1991a,Gu1998,Wolter2006,Kaparis2010a}.
To lift each variable, a knapsack problem is required to be solved to compute the lifting coefficient. Under certain conditions, the \highlight{sequential} LCI can be computed in polynomial time, see \cite{Zemel1989,Nemhauser1988,Gu1995}. \rev{In general, however,} the complexity of computing \highlight{a sequential} LCI is still unknown. \rev{This was explicitly mentioned in \cite{Gu1999} as an open question.}
\vspace{0.2cm}

``We show that this dynamic programming algorithm may take exponential time to compute a sequential LCI that is not simple. It is still an open question whether an arbitrary LCI can be computed in polynomial time for a given minimal cover $ C $."
\vspace{0.2cm}

\rev{The above question was} also cited as an open question in \cite{Gu1995,VanHoesel2002,Atamturk2002c,Gupta2005}. We will give a
negative answer to the question by showing that the problem of computing a \highlight{sequential} LCI is $ \mathcal{NP} $-hard. Thus, unless $ \mathcal{P}=\mathcal{NP} $, there exists no polynomial time algorithm to compute \highlight{a sequential} LCI.

\rev{This paper is organized as follows. In Section \ref{sec-clci}, we review how to compute \highlight{a sequential} LCI.
In Section \ref{sec-example}, we \rev{describe} the elegant example by Gu \cite{Gu1995}, which provides exponential lifting coefficients. The main result is given in the last section, which shows the $\mathcal{NP} $-hardness of the problem of computings \highlight{a sequential} LCI. }
\section{Computing \rev{a} \highlight{sequentially} lifted cover inequality}
\label{sec-clci}
Consider the knapsack set $X : = \{ x \in \mathbb{B}^n \, : \, a^Tx  \leq b  \}$,
where $ b\in\mathbb{Z}_+$ and $a=(a_1,\cdots,a_n)^T\in\mathbb{Z}_{+}^{n}$ are given.
A subset $ C \subseteq N := \{ 1,\ldots,n \}$ is called a {\textit{cover}} of $ X $ if
the sum of the items $a_i$'s over $C$ exceeds the {\emph{knapsack capacity}} $b$; i.e.,
$ \sum_{i\in C}a_i > b$. A cover $C$ is a {\textit{minimal cover}} if and only if
\begin{equation*}
	\sum_{i \in C \backslash\{ j \}} a_i \leq b\quad \text{for all } j \in C.
\end{equation*}

For any subsets $ N_0 $ and $ N_1 $ of $ N$ with $ N_0 \cap N_1 = \varnothing $, denote $ X(N_0, N_1) $ be the following restriction set of $X$,
\begin{equation*}
X(N_0, N_1) = X \cap \{ x \in \mathbb{B}^n \,:\, x_i = 0\ \text{for } i \in N_0; \ x_i =1\ \text{for } i \in N_1 \}.
\end{equation*}
It is easy to see that, for each cover $C$, the {\textit{cover inequality}}
\begin{equation}
	\label{cover}
	\sum_{i \in C} x_i \leq c -1
\end{equation}
is valid for $  X(N \backslash C, \varnothing) $, where $ c := |C| $ is the cardinality  of $ C $.
The cover inequality \eqref{cover} is facet defining for $\text{conv}(X(N \backslash C, \varnothing))$, which is the convex hull of $ X(N\backslash C, \varnothing)$, if and only if $C$ is a minimal cover (see for example \hilight{\cite{Conforti2014}}).

We can consider to fix some variables to be ones as well. Assume that $(C, N_0, N_1)$ is a partition of $ N $ and denote
$\bar{b}= b - \sum_{i\in N_1} a_i$. In this case, the inequality \eqref{cover} is facet defining for $ \text{conv}(X(N_0,N_1)) $
if and only if $C$ is a minimal cover of $ X(N_0, N_1) $; i.e.,
\begin{eqnarray*}
	\sum_{i\in C}a_i > \bar{b}; \quad& &   \sum_{i \in C \backslash\{ j \}} a_i \leq \bar{b} \quad \text{for all } j \in C.
\end{eqnarray*}
Throughout this paper, we shall assume that $(C, N_0, N_1)$ be a partition of $ N $ and $C$ is a minimal cover of $ X(N_0, N_1) $.

In general, however, the inequality \eqref{cover} may not be valid for $ X $ if $ N_1 \neq \varnothing $. Furthermore, even if
$ N_1 = \varnothing $, such an inequality may not represent a facet of $ \text{conv}(X) $.
To obtain a strong valid inequality, we can lift the variables in $ N_0  \cup N_1 $ one at a time by solving an optimization problem sequentially.
More precisely, let $ l_1,\ldots,l_{n-c} $ be a lifting sequence such that $ N_0 \cup N_1 = \{ l_1,\ldots,l_{n-c} \}$ and
\begin{equation*}
	\sum_{i\in C} x_i + \sum_{i \in N'_0} \alpha_i x_i + \sum_{i\in N'_1} \beta_i x_i  \leq c -1 + \sum_{i \in N'_1}\beta_i
\end{equation*}
be the inequality obtained so far, where $ N'_0 \subseteq N_0 $, $ N'_1 \subseteq N_1 $ and $ N'_0 \cup N'_1 = \{ l_1,\ldots,l_j \} $. To lift the variable $ x_{l_{j+1}} $, we need to solve a knapsack subproblem depending on whether $l_{j+1}$ belongs to $N_0$ or $N_1$.
\highlight{We follow \cite{Gu1998} in referring to lifting the variable $ x_{l_{j+1}} $
as {\textit{up-lifting}} if $ l_{j+1} \in N_0 $ and {\textit{down-lifting}} if
		$ l_{j+1} \in N_1 $.}

\begin{itemize}
	\item [(i)] {\textit{Up-lifting}}. If $ l_{j+1} \in N_0 $, compute the lifting coefficient $ \alpha_{l_{j+1}} $ by solving the lifting problem
\begin{equation}
\label{uplift}
\begin{aligned}
\alpha_{l_{j+1}} =  & \ {\text{min}}
  & &  c - 1 + \sum_{i \in N'_1} \beta_i -  \sum_{i \in C}x_i - \sum_{i\in N'_0 }\alpha_i x_i  -\sum_{i \in N'_1} \beta_i x_i  \\
  &  \ \text{s.t.}  & &           \sum_{i \in C \cup N'_0 \cup N'_1} a_i x_i \leq \bar{b} + \sum_{i \in N'_1}a_i - a_{l_{j+1}} , \ x\in \mathbb{B}^{c+j}  .
\end{aligned}
\end{equation}
	\item [(ii)]  {\textit{Down-lifting}}. If $ l_{j+1} \in N_1 $, compute the lifting coefficient $ \hilight{\beta_{l_{j+1}}} $ by solving the lifting problem
	\begin{equation}
	\label{downlift}
	\begin{aligned}
	\beta_{l_{j+1}} =    &\  {\text{max}}
	& & \sum_{i \in C}  x_i + \sum_{i\in N'_0 }\alpha_i x_i  +\sum_{i \in N'_1} \beta_i x_i - c + 1 - \sum_{i \in N'_1} \beta_i    \\
	&  \ \text{s.t.}          & & \sum_{i \in C \cup N'_0 \cup N'_1} a_i x_i \leq \bar{b} + \sum_{i \in N'_1}a_i + a_{l_{j+1}}  , \ x\in \mathbb{B}^{c+j}.
	\end{aligned}
	\end{equation}
\end{itemize}
After having lifted all the variables, we obtain the \highlight{sequential LCI:}
\begin{equation}
	\label{lci}
	\sum_{i\in C} x_i + \sum_{i \in N_0} \alpha_i x_i + \sum_{i\in N_1} \beta_i x_i  \leq c -1 + \sum_{i \in N_1}\beta_i.
\end{equation}
\highlight{The procedure to obtain inequality of \eqref{lci} is first described implicitly in \cite{Wolsey1975}.}
See for example \highlight{\cite{Nemhauser1988,Gu1995,Kaparis2011}} for more details about \highlight{sequential} LCI. Here we just
notice that different lifting sequences may lead to different \highlight{sequential} LCIs.

The inequality \eqref{lci} is called a non-project LCI if $ N_1 = \varnothing $ and a project LCI if $ N_1 \neq \varnothing $ (see \cite{Gu1995}). Given a lifting sequence, the non-project LCI can be computed (\cite{Zemel1989}) in the complexity of $ \mathcal{O}(c n) $, where $c=|C|$ again.
For the project LCI, if $ C \cup N_1 $ is a minimal cover of $ X $ and the lifting sequence is enforced to $\{l^0_1,\ldots,l^0_r, l^1_1,\ldots,l^1_{|N_1|},l^0_{r+1},\ldots,l^0_{|N_0|}\} $, where $ l^j_i \in N_j  $ for $ i =1,\ldots,|N_j|$ and $j=0,1$ and
$r$ is a given integer between $1$ and $|N_0|$, Gu \cite{Gu1995} proved that
\highlight{\eqref{lci}} can be obtained in the complexity of $ \mathcal{O}(cn^3) $.
As mentioned before, however, the complexity of computing a \highlight{sequential} LCI with an arbitrary lifting sequence is still unknown.

We close this section by noting that the \highlight{sequential} LCI is invariant under scaling.
\begin{myobservation}
	\label{invariantproperty}
	Given the same partition and the lifting sequence, multiplying a positive integer to the knapsack constraint gives the same \highlight{sequential} LCI.
\end{myobservation}
\section{The example by Gu (1995)}
\label{sec-example}
In this section, we describe the elegant example constructed by Gu \cite{Gu1995},
which leads the lifting coefficients to be exponential. It is related to the
following vector $f\in\mathbb{Z}^{2r+1}$, where $r$ is a given positive integer.
\begin{equation}
	\label{exampleitems}
	f_1 = 1, \, f_2 =1 , \,f_3=1, \, \text{and} \, f_i = f_{i-2}+ f_{i-1} ,\ \text{for} \ i = 4,\ldots,2r + 1.
\end{equation}
\highlight{Notice that $ f $ is analogous to the Fibonacci sequence where the only difference
is that the first element in the Fibonacci sequence is $ 0 $.}
We give two facts on the vector $f$, which can easily be verified by induction.
\begin{myobservation}
	\label{inductionformula}
	For $ j=3,\dots,2r+1 $, it holds that $ f_{j} = \sum_{i=1}^{j-2} f_{i} $	.
\end{myobservation}
\begin{myobservation}
	\label{boundsf}
	For $ j=3,\dots,2r+1 $, it holds that  $ \hilight{\frac{1}{4} (\sqrt{2} -1) (\sqrt{2})^{j} \leq f_j \leq 2^j} $.
\end{myobservation}

Consider the knapsack set $ X $ with $ 2r+1 $ variables, where the coefficients of the knapsack constraint
are $ f_1,\ldots,f_{2r+1}$ in \eqref{exampleitems} and the associated knapsack capacity is $ b = \sum_{i=1}^{2r}f_i$.
Consider the partition $ (C,N_0, N_1) $ of $ \{1,\ldots,2r+1 \} $ with $ C=\{1,2\} $, $ N_0=\{ 4, 6,8,\ldots 2r \} $, and $ N_1 = \{ 3,5,7,\ldots, 2r+1 \} $. Since $f_1 =1 $, $ f_2=1 $, and
\begin{eqnarray*}
	\bar{b}&=&b-\sum_{i\in N_1}f_i  =b-\sum_{i=1}^{r}f_{2i+1} =b-f_3 - \sum_{i=2}^r (f_{2i-1} + f_{2i})\\
	&=& b-f_3 - \sum_{i=3}^{2r}f_i = 1,
\end{eqnarray*}
we know that $ C $ is a minimal cover of $ X(N_0,N_1)= \{ x\in \mathbb{B}^2 \, : \, x_1 + x_2 \leq 1\} $. Now consider the lifting sequence $ \{3,4,\ldots,2r+1\}$; i.e., $ l_j = j+2 $ for $ 1 \leq j \leq 2r-1 $. The following lemma is due to \cite{Gu1995}. For completeness, we provide a proof.
\begin{lemma}[Gu 1995]
	\label{guexample}
	(i) $ \alpha_{i} = f_{i} $ for $ i \in N_0 $; (ii) $ \beta_{i} = f_{i} $ for $ i \in N_1 $.
\end{lemma}
\begin{proof}
	We proceed by induction. At first, the lifting problem of the variable $ x_3 $ reads
	\begin{equation*}
		\begin{aligned}
			\hilight{\beta_{3}} =  & \ \hilight{\text{max}}
			& &   x_1 + x_2 -1  \\
			&  \ \text{s.t.}  & &         x_1 + x_2 \leq 2, \ x \in \mathbb{B}^2   .
		\end{aligned}
	\end{equation*}
	Hence $ \beta_{3}= 1 = f_3 $.

	Assume that $ \alpha_i = f_i $ for $ i \leq j $, $ i \in N_0 $ and $ \beta_i = f_i  $ for $ i \leq j $, $ i\in N_1 $, respectively. Now we consider lifting the variable $ x_{j+1} $. If $ {j + 1} \in N_0 $, the problem \eqref{uplift} reduces to
	\begin{equation}
		\label{exauplift}
		\begin{aligned}
			\alpha_{j+1} =  & \ {\text{min}}
			& &  1 + \sum_{i=1}^{(j-1)/2} f_{2i+1} - \sum_{i=1}^{j}f_i x_i  \\
			&  \ \text{s.t.}  & &          \sum_{i=1}^j f_i x_i  \leq 1 + \sum_{i=1}^{(j-1)/2} f_{2i+1} - f_{j+1},\ x\in \mathbb{B}^{j}  .
		\end{aligned}
	\end{equation}
	Since
	\begin{equation*}
		\sum_{i=1}^{(j-1)/2}f_{2i+1} = f_3 + \sum_{i=2}^{(j-1)/2}(f_{2i} + f_{2i-1})= f_ 3 + \sum_{i=3}^{j-1} f_i  = \sum_{i=1}^{j-1}f_i -1 = f_{j+1} -1,
	\end{equation*}
	where the last equality follows from Observation \ref{inductionformula}, the feasibility of the problem \eqref{exauplift} requires
$ x_i = 0 $, $ 1\leq i \leq j $ and hence $\alpha_{j+1} = f_{j+1} $. If $ j+1 \in N_1 $, the problem \eqref{downlift} reduces to
	\begin{equation}
		\label{exadownlift}
		\begin{aligned}
			\beta_{j+1} =    &\  {\text{max}}
			& & \sum_{i=1}^{j}f_i x_i - 1 - \sum_{i=1}^{j/2 -1} f_{2i+1}    \\
			&  \ \text{s.t.}          & & \sum_{i=1}^j f_i x_i  \leq 1 + \sum_{i=1}^{j/2-1} f_{2i+1} + f_{j+1}, \ x\in \mathbb{B}^{j}.
		\end{aligned}
	\end{equation}
 It is easy to verify that $\sum_{i=1}^j f_i  = 1 + \sum_{i=1}^{j/2-1}f_{2i+1} + f_{j+1}$. Hence the all-one vector is feasible
 and solves the problem \eqref{exadownlift}, yielding $\beta_{j+1}=f_{j+1}$. So the statement is true for $j+1$ as well. By induction, this lemma holds.

	\qed
\end{proof}

\vspace{0.2cm}
Lemma \ref{guexample} indicates that the \highlight{sequential} LCI for this specific example is
\begin{equation}
	\label{liftedinequality}
	\sum_{i=1}^{2r+1} f_i x_i \leq \sum_{i=1}^{2r}f_i.
\end{equation}
By Observation \ref{boundsf}, \hilight{the input size of this example is polynomial, but} the lifting coefficients $\{f_i\}$ are exponential with $r$.
This example by Gu will play an important role in the coming complexity analysis.
\section{$ \mathcal{NP} $-hardness of computing a \highlight{sequentially} lifted cover inequality}
\label{sec-np}
In this section, we show the $ \mathcal{NP} $-hardness of the problem of computing \highlight{a sequential} LCI.
To begin with, we give a basic property of the vector $ f $ in \eqref{exampleitems}.
\begin{lemma}
	\label{anyinteger}
	Let $ f  $ be defined as in \eqref{exampleitems}, where $r\ge 1$ is given.
    For any $ \tau \in \mathbb{Z}_+$ satisfying $ 0 \leq \tau \leq \sum_{i=1}^{2r+1} f_{i}$, there
	exists a subset $ S \subseteq \{1,\ldots,2r+1\} $ such that $ \tau = \sum_{i\in S} f_i  $.
\end{lemma}
\begin{proof}
	We proceed by induction on $r$. The result apparently holds for $r=1$.
	Assume that the result is true for some $r\geq 1 $. To verify the result for $r+1$,
    it suffices to consider the case that
    $  \sum_{i=1}^{2r+1} f_i  <  \tau \leq \sum_{i=1}^{2(r+1)+1} f_i$.
	In fact, from Observation \ref{inductionformula}, we have that
	\begin{eqnarray*}
		\sum_{i=1}^{2r+1} f_i  = f_{2r+3},\qquad&&\sum_{i=1}^{2r+3} f_i = f_{2r+3} + f_{2r+2}+ f_{2r+3}.
	\end{eqnarray*}
	So $ f_{2r+3} < \tau \leq f_{2r+3} + f_{2r+2}+ f_{2r+3}$. Let us define $ \tau_1 $ as
	\begin{equation}
	\label{tau1}
	\tau_1 =\left\{\begin{array}{ll}
	\tau- f_{2r+2}- f_{2r+3},  & \ \text{if} \ \tau > f_{2r+2}+ f_{2r+3}; \\
	\tau - f_{2r+3}, &\ \text{if} \  \tau \leq f_{2r+2}+ f_{2r+3}.
	\end{array}\right.
	\end{equation}

	Then it is easy to see that $ \tau_1 \leq f_{2r+3} = \sum_{i=1}^{2r+1} f_i $.
	By the induction assumption, there exists $S \subseteq \{ 1,2,3,...,2r+1 \} $ such that $ \tau_1 = \sum_{i \in S} f_i $. By picking
    \highlight{one more element $ f_{2r+3} $ and the possible element $ f_{2r+2} $}, we know that there exists some subset of $\{ 1,2,3,...,2r+3\}$ such that the sum of $f_i$'s over this subset is exactly $\tau$. Thus the result holds for $r+1$. By induction, this lemma is true.

	\qed
\end{proof}
\vspace{0.2cm}

\noindent
\highlight{Next, we introduce the \highlight{restricted partition problem (RPP)}, which is a variant of the partition problem \cite{GAREY1979}. Comparing to the partition problem,
	the RPP problem restricts the total sum of all element to some specific values and
	allows the sum of the elements in the subset equals to one more value.
The RPP problem is shown to be $ \mathcal{NP} $-hard in
	\cite{Klabjan1998a}.}
\begin{quote}
{\bf Problem RPP}.\ Given a nonnegative integer $ m $, a finite set $ K $ \highlight{of $ k $ elements with value $ \omega_i \in \mathbb{Z}_+ $} for the $ i $-th element \highlight{and}
	$ \sum_{i\in K} \omega_i  = 2(2^{m+1} -   1) $, does there exist a subset $ T \subseteq K $ such that $ \sum_{i \in T} \omega_i = 2^{m+1} - 1 $ or $ \sum_{i \in T} \omega_i = 2^{m+1} - 2 $ ?
\end{quote}
\highlight{For convenience, define $\lambda = 2^{m+1} -1 $ and then $ \sum_{i\in K} \omega_i = 2\lambda  $.
The RPP problem is still $ \mathcal{NP} $-hard when
$ 1 \leq \omega_i \leq \lambda - 1 $. To see this,
suppose there exists some $ j \in K $ such that $ \omega_j \geq \lambda $. We have the following two cases.
\begin{itemize}
	\item [(i)]
		$ \lambda \leq \omega_j \leq \lambda+1 $. In this case, it follows that
		$ \sum_{i \in K \backslash\{j\}} \omega_i = \lambda-1 $ or $ \sum_{i \in K \backslash\{j\}} \omega_i = \lambda  $. Thus, $ K \backslash\{j\}$ is the desired subset and the answer to the RPP problem is yes.
	\item [(ii)]
	$ \omega_j \geq \lambda+2 $. In this case, it is easy to see that the answer to the RPP
	problem is no.
\end{itemize}
As we can see, both cases can be solved in polynomial time.
Therefore, the statement follows and in the following, we assume that
 $ 1 \leq \omega_i \leq \lambda - 1 $.
}
We are now ready to present the main result of this paper; i.e., provide an $ \mathcal{NP} $-hardness proof for computing \highlight{a sequential} LCI. The basic idea is as follows. Firstly, we adopt Gu's example (see Section \ref{sec-example}) to make the lifting coefficients exponential. Secondly, some variables fixed at zero will be lifted, where the lifting coefficients can easily be obtained. Finally, we lift a variable fixed at one, \highlight{where the objective value of corresponding lifting problem is equal to some specific value if and only if the answer to the RPP problem is yes and hence is $ \mathcal{NP} $-hard}.
\begin{theorem}
	The problem of computing \highlight{a sequential} LCI is $ \mathcal{NP} $-hard.
\end{theorem}
\begin{proof}
    For an RPP instance with $m\ge 0$ and $K= \{ 2r+4,\ldots, n-1 \} $, we shall construct a problem of computing a \highlight{sequential} LCI
    in polynomial time.

	We construct a problem of computing a \highlight{sequential} LCI as follows.
Set $ r = m + 6 $,

	$ n=2r+4 +k $, and $ b = \sum_{i=1}^{2r+1} \lambda f_i  + \lambda (3\lambda + 6 )   $, where $ f$ is defined as in \eqref{exampleitems}. Set the coefficients of the knapsack constraint as
	\begin{equation*}
	a_i =\left\{\begin{array}{ll}
	\lambda f_i ,  & \ \text{for} \ i=1,\ldots,2r+1; \\
	\lambda(\lambda+3) + 1 , &\ \text{for} \  i=2r+2; \\
	\lambda(\lambda+3) - 1  , &\ \text{for} \  i=2r+3; \\
	\omega_i(\lambda + 1)  , &\  \text{for} \  i=2r+4,\ldots,n-1; \\
\lambda ( 3\lambda + 6 + f_{2r+1}) , &\ \text{for} \  i=n.
	\end{array}\right.
	\end{equation*}
 Define the partition $ (C, N_0, N_1) $ with $ C = \{1,2\} $, $ N_0 = \{ 4,6,8,\ldots, 2r,2r+2, 2r+3,\ldots,n-1  \} $, and $ N_1 = \{3,5,7,\ldots, 2r+1, n\} $. Finally, let the lifting sequence be $ \{\hilight{3,\ldots,n}\}  $.
We shall prove that the lifting coefficient $ \hilight{\beta_n} = f_{2r+1} + 3\lambda + 5  $ if and only if the answer to the RPP instance is yes.

Before doing this, we note that the input size of this \highlight{instance} is polynomial of that of the RPP instance. To see this, let $ L $ be the input size of the RPP instance. \highlight{It follows immediately that $ k = \mathcal{O}(L) $.
We next show that the number of elements $ n $ is satisfied with $ n = \mathcal{O}(L) $.
	Since the input size of a positive integer $ t $ is $  \log_2 (t+1) $, it follows that
	\begin{equation*}
		\log_2(2^{m+2} -2+1) =  \log_2 (\sum_{i \in K} \omega_i +1 ) \leq \sum_{i\in K}\log_2(\omega_i + 1) \leq L,
	\end{equation*}
	where the first inequality follows from $\omega_i +1 \geq 2  $ for all $ i \in K $.
	Thus, we have that $ m \leq L -2 $, which further implies that $ r=m+6 = \mathcal{O}(L) $.
	This, combined with the fact that $ k = \mathcal{O}(L) $, indicates
	$$
		n = 2r+ 1 + k =\mathcal{O}(L).
	$$
	Finally, it can be easily verified that $ a_i= \mathcal{O}(L^2) $ for $ i =1,\ldots, n $
	and $ b = \mathcal{O}(L^2)  $. This proves that the input size of the
	constructed instance is polynomial of that of the RPP instance.}
	For preparation \highlight{of the proof}, we also verify

\begin{equation}
\begin{aligned}
f_{2r+1} & \geq & & \dfrac{1}{4} (\sqrt{2} -1) (\sqrt{2})^{2r+1}=  \dfrac{1}{4} (2-\sqrt{2}) 2^r= \dfrac{1}{4} (2-\sqrt{2}) 2^{m+6} &  \\[8pt]
& > &&  2^{m+2} + 4 = 2\lambda + 6. &  \label{lambdaf}
\end{aligned}
\end{equation}

	In the following, we consider the lifting procedure. By construction, the knapsack inequality of this instance reads
	\begin{eqnarray*}
		&&\sum_{i=1}^{2r+1} \lambda f_i x_i  + [\lambda(\lambda+3)+1] x_{2r+2} + [\lambda(\lambda+3)-1] x_{2r+3} +\sum_{i=2r+4}^{n-1} \omega_i (\lambda+1) x_i     \\
		&&\qquad + \lambda ( 3\lambda + 6 + f_{2r+1}) x_n \leq \sum_{i=1}^{2r+1} \lambda f_i  + \lambda (3\lambda + 6 ).
	\end{eqnarray*}
	Since $ a_1 = \lambda  $, $ a_2 = \lambda $, and
	\begin{equation*}
		\begin{aligned}
			b- \sum_{i\in N_1}a_i = \sum_{i=1}^{2r+1} \lambda f_i  + \lambda (3\lambda + 6 )- \sum_{i=1}^{r} \lambda f_{2i+1} -\lambda ( 3\lambda + 6 + f_{2r+1}) &   \\
			= \sum_{i=1}^{2r} \lambda f_i - \sum_{i=1}^{r}\lambda f_{2i+1} =  \sum_{i=1}^{2r} \lambda f_i  -\lambda f_3 - \sum_{i=2}^{r} \lambda (f_{2i-1}+f_{2i})= \lambda,
		\end{aligned}
	\end{equation*}
	we know that $ C $ is a minimal cover of $ X(N_0,N_1) $ with the cover inequality
	\begin{equation}
		\label{proofcover}
		x_{1}+x_2 \leq 1.
	\end{equation}

\vspace*{0.3cm}
\noindent {\bf Step 1. Lifting the variables $ \boldsymbol{x_3,\ldots,x_{2r+1}} $}

\vspace*{0.2cm}\noindent Starting with the cover inequality \eqref{proofcover}, we know from Observation \ref{invariantproperty} and the inequality \eqref{liftedinequality} that, after lifting the variables $ x_3,\ldots, x_{2r+1} $, the inequality is
	\begin{equation*}
		\sum_{i=1}^{2r+1} f_i x_i \leq \sum_{i=1}^{2r} f_i.
	\end{equation*}

\vspace*{0.3cm}
\noindent {\bf Step 2. Lifting the variables $ \boldsymbol{x_{2r+2}} $ and $ \boldsymbol{x_{2r+3}} $}

\vspace*{0.2cm}\noindent We first consider the variable $ x_{2r+2} $. The associated lifting problem is

	\begin{equation}
		\label{lifting2r2}
		\begin{aligned}
			\alpha_{2r+2} =  \ & {\text{min}} &   & \sum_{i=1}^{2r} f_i  - \sum_{i=1}^{2r+1} f_i x_i                                                 \\

			                & \text{s.t.}  &   & \sum_{i=1}^{2r+1} \lambda f_i x_i  \leq \lambda \sum_{i=1}^{2r}  f_i -[\lambda(\lambda+3) +1] ,\ x \in \mathbb{B}^{2r+1}.
		\end{aligned}
	\end{equation}
	The problem \eqref{lifting2r2} is feasible at the zero vector since

	\begin{equation*}
\lambda \sum_{i=1}^{2r}  f_i -[\lambda(\lambda+3)+1] \geq \lambda (\sum_{i=1}^{2r}  f_i - \lambda- 4)  > \lambda (f_{2r+1}- \lambda - 4)> 0,
	\end{equation*}
where the last inequality follows from \eqref{lambdaf}. Let $ \bar{x} $ be the optimal solution of \eqref{lifting2r2}. Since $ f_{i} \in \mathbb{Z} $ for $ i =1,\ldots,2r+1 $, its feasibility requires
	\begin{equation*}
		\sum_{i=1}^{2r+1} f_i\bar{x}_i \leq \lfloor\frac{ \lambda \sum_{i=1}^{2r} f_i - \lambda(\lambda+3) - 1}{\lambda} \rfloor =  \sum_{i=1}^{2r} f_i  - \lambda - 4.
	\end{equation*}
	On the other hand,  from Lemma \ref{anyinteger}, we can always find an $ \bar{x} $ such that
	\begin{equation*}
	 	\sum_{i=1}^{2r+1} f_i \bar{x}_i = \sum_{i=1}^{2r}  f_i - \lambda - 4.
	\end{equation*}
	The optimality of $ \bar{x} $ gives that $ \alpha_{2r+2} = \sum_{i=1}^{2r}  f_i - \sum_{i=1}^{2r+1} f_i \bar{x}_i = \lambda + 4 $.

	Similarly, the lifting problem of $ x_{2r+3} $ reads
		\begin{align}
			\alpha_{2r+3} = \ & {\text{min}} &   & \sum_{i=1}^{2r} f_i  - \sum_{i=1}^{2r+1} f_i x_i - (\lambda + 4)x_{2r+2}                                                    \nonumber        \\

			                & \text{s.t.}  &   & \sum_{i=1}^{2r+1} \lambda f_i x_i + [\lambda(\lambda+3) +1] x_{2r+2}  \leq \lambda \sum_{i=1}^{2r}  f_i - [\lambda(\lambda+3) -1] , \nonumber\\
			                &              &   & x \in \mathbb{B}^{2r+2}. \label{lifting2r3}
		\end{align}
	Then $ \alpha_{2r+3} = \lambda + 2 $, which is achieved at an optimal solution $ \bar{x} $ satisfying
	$ \bar{x}_{2r+2}=1 $ and $ \sum_{i=1}^{2r+1}f_i \bar{x}_i = \sum_{i=1}^{2r} f_i  - 2\lambda-6 $.

\vspace*{0.3cm} \noindent {\bf Step 3. Lifting the variables $ \boldsymbol{x_{2r+4},\ldots,x_{n-1}} $}

\vspace*{0.2cm}\noindent	We shall show that $ \alpha_{i} = \omega_i $ for all $ i\in K$ by induction.
	At first, consider the lifting \hilight{of} the variable $ x_{2r+4} $. This requires to solve the problem
	\begin{equation}

		\label{lifting2r4}
		\begin{aligned}
			\alpha_{2r+4}= \ & {\text{min}} &   & \sum_{i=1}^{2r} f_i  - \sum_{i=1}^{2r+1} f_i x_i - (\lambda + 4)x_{2r+2}  - (\lambda + 2)x_{2r+3}        \\

			               & \text{s.t.}  &   & \sum_{i=1}^{2r+1} \lambda f_i x_i + [\lambda(\lambda+3) +1] x_{2r+2} +  [\lambda(\lambda+3)-1] x_{2r+3} \\
			               &              &   & \leq \lambda \sum_{i=1}^{2r}  f_i - \omega_{2r+4} (\lambda + 1) , \ x \in \mathbb{B}^{2r+3}.
		\end{aligned}
	\end{equation}
	If $ \hat{x} $  is an optimal solution of the problem \eqref{lifting2r4} with $ \hat{x}_{2r+3}=1 $,  by the feasibility and \eqref{lambdaf}, we have
	\begin{equation*}
		\sum_{i=1}^{2r+1} f_{i}\hat{x}_i  + \lambda + 2 <  \sum_{i=1}^{2r} f_{i} + \lambda + 2 <\sum_{i=1}^{2r+1} f_{i}.
	\end{equation*}
	This, together with Lemma \ref{anyinteger}, indicates that we can define a new feasible point $ \bar{x} $ such that $ \bar{x}_{2r+2}= \hat{x}_{2r+2} $, $  \bar{x}_{2r+3}= 0$, and $ \sum_{i=1}^{2r+1} f_i \bar{x}_i =\sum_{i=1}^{2r+1} f_i \hat{x}_i + \lambda + 2$.
	It is easy to see that $ \bar{x} $ and $ \hat{x}$ give the same objective values. Hence, we can assume that $ x_{2r+3} = 0 $ in the problem \eqref{lifting2r4}. Furthermore, since $ \omega_{2r+4} \leq  \lambda -1 $, similar to the problem \eqref{lifting2r3},
we can show that the optimal value of \eqref{lifting2r4} is $ \alpha_{2r+4} = \omega_{2r+4}$, which is achieved at an optimal solution
 $ \bar{x} $ satisfying $\bar{x}_{2r+3}=0$, $ \bar{x}_{2r+2}=1$, and $ \sum_{i=1}^{2r+1}f_i \bar{x}_i = \sum_{i=1}^{2r} f_i  - (\lambda+4)-\omega_{2r+4}$.

 Now assume that $ \alpha_{i} = \omega_{i} $ for $ 2r+4 \leq i \leq j $ and $j<n-1$ and consider the lifting problem of $ x_{j+1} $:

\begin{spacing}{0.9}
	\begin{align}
	\alpha_{j+1}= \ & {\text{min}} &   & \sum_{i=1}^{2r} f_i  - \sum_{i=1}^{2r+1} f_i x_i - (\lambda + 4)x_{2r+2}  - (\lambda + 2)x_{2r+3} - \sum_{i=2r+4}^{j} \omega_i x_i  \nonumber      \\

	& \text{s.t.}  &   & \sum_{i=1}^{2r+1} \lambda f_i x_i + [\lambda(\lambda+3) +1] x_{2r+2} +  [\lambda(\lambda+3) -1] x_{2r+3} +   \nonumber\\
	&              &   & \sum_{i=2r+4}^{j}  \omega_i(\lambda+1)x_i  \leq \lambda \sum_{i=1}^{2r}  f_i - \omega_{j+1} (\lambda + 1) , \ x \in \mathbb{B}^{j}. \label{lifting2rj1}
	\end{align}
\end{spacing}
\noindent We claim that there exists an optimal solution $ \bar{x} $ such that $ \bar{x}_{2r+3} = \bar{x}_{2r+4}=\cdots= \bar{x}_j = 0 $. To see this, suppose that an optimal solution $\hat{x}$ is such that $ \hat{x}_t = 1 $ for some $ t \subseteq [2r+3,\,j]$. Analogously, define a new point $ \bar{x} $ with $ \hilight{\bar{x}_{i}= \hat{x}_{i}} $ for $ 2r+2 \leq i \leq j$ and $i\ne t$, $ \bar{x}_{t} = 0 $, and $ \sum_{i=1}^{2r+1} f_i \bar{x}_i =\sum_{i=1}^{2r+1} f_i \hat{x}_i + \theta_t$, where
	\begin{equation*}
	\theta_t=\left\{\begin{array}{ll}
	\lambda + 2,  &\ \text{if} \ t = 2r+3; \\
	\omega_t, &\ \text{otherwise}.
	\end{array}\right.
	\end{equation*}
	By simple calculations, $ \hilight{\bar{x}} $ is feasible to the problem \eqref{lifting2rj1} and gives the same objective value as $ \bar{x} $. Similar to the problem \hilight{\eqref{lifting2r3}}, we can verify that $ \alpha_{j+1} = \omega_{j+1}$.
Thus by induction, we have that $ \alpha_{i} = \omega_i $ for all $ i\in K$.

\vspace*{0.3cm} \noindent {\bf Step 4. Lifting the variable $ \boldsymbol{x_{n}} $}

\vspace*{0.2cm}\noindent	Finally, we concentrate on lifting the variable $ x_n $.
	The lifting problem is

		{
		\begin{spacing}{0.9}
		\begin{align}
			\beta_{n} = \  & {\text{max}} &   & \sum_{i=1}^{2r+1} f_i x_i + (\lambda + 4)x_{2r+2} + (\lambda + 2)x_{2r+3}  + \sum_{i=2r+4}^{n-1}  \omega_i x_i  -\sum_{i=1}^{2r} f_i\nonumber \\

			             & \text{s.t.}  &   & \sum_{i=1}^{2r+1} \lambda f_i x_i + [\lambda(\lambda+3)+1] x_{2r+2}+  [\lambda(\lambda+3)-1] x_{2r+3} +                        \nonumber    \\
			             &              &   & \sum_{i=2r+4}^{n-1}\omega_i(\lambda +1)x_i  \leq \lambda \sum_{i=1}^{2r+1}  f_i + \lambda(3\lambda +6) , \ x \in \mathbb{B}^{n-1}.	\label{liftingn}
		\end{align}
	\end{spacing}
}
\noindent For convenience, denote $ g(x) $ to be the objective function in the above problem. Consider the point $ \check{x} $ with $ \check{x}_i =1 $ for $ 2r+2\leq  i\leq n-1  $ and $ \sum_{i=1}^{2r+1} f_i \check{x}_i = \sum_{i=1}^{2r+1} f_i  - \lambda - 2$.
By Lemma \ref{anyinteger}, such a point must exist. We can check that $ \check{x} $ is feasible to the problem \eqref{liftingn} and
$ g(\check{x}) = f_{2r+1} + 3\lambda+4$. Furthermore, for a binary vector $x$, if at least one of the two components $x_{2r+2}$ and $x_{2r+3}$ is equal to zero, we have that
$$ g(x) \leq  \sum_{i=1}^{2r+1} f_i  + \lambda + 4 + \sum_{2r+4}^{n-1} \omega_i - \sum_{i=1}^{2r}f_i  = f_{2r+1} + 3\lambda + 4. $$
Thus to seek better values for $\beta_n$, we may set $ x_{2r+2} = x_{2r+3} = 1 $. In this case, the problem \eqref{liftingn} reduces to

		\begin{align}
			\beta_{n} = \ & {\text{max}} &   & \sum_{i=1}^{2r+1} f_i x_i + \sum_{i=2r+4}^{n-1}  \omega_i x_i + 2\lambda + 6 -    \sum_{i=1}^{2r} f_i    	\label{liftingn-1}           \\

			             & \text{s.t.}  &   & \sum_{i=1}^{2r+1} \lambda f_i x_i  + \sum_{i=2r+4}^{n-1}\omega_i(\lambda +1)x_i  \leq \lambda \sum_{i=1}^{2r+1}  f_i + \lambda^2,\ x \in \mathbb{B}^{n-3}. \nonumber
		\end{align}

Now assume that $\bar{x}$ is an optimal solution of \eqref{liftingn-1}. Denote $ p = \sum_{i=2r+4}^{n-1} \omega_i \bar{x}_i $. It is easy to see
that $ p \leq 2\lambda $. Consider the following four cases.
	\begin{enumerate}
		\item [(a)] $ p  \le \lambda - 2 $. In this case, the knapsack constraint in the problem \eqref{liftingn-1} is trivially satisfied and the optimality of $ \bar{x} $ implies
		\highlight{$ \bar{x}_i = 1$ for $ i = 1,\ldots, 2r+1$}. \highlight{This further indicates} that
		$$ g(\bar{x}) = \sum_{i=1}^{2r+1} f_i + p + 2\lambda + 6 -  \sum_{i=1}^{2r} f_i = f_{2r+1} + p +  2\lambda + 6 \le f_{2r+1} + 3\lambda + 4.
		$$
		\item  [(b)]
		      $ p  = \lambda - 1 $.  Similar to the case (a), we have that 	$ g(\bar{x}) =  f_{2r+1}+3\lambda + 5  $.
		\item  [(c)]
		      $ p = \lambda $. In this case, the feasibility of $\bar{x}$ indicates that
		      $$
		      \sum_{i=1}^{2r+1} f_i \bar{x}_i  \leq  [  \lambda \sum_{i=1}^{2r+1} f_i + \lambda^2  - \lambda(\lambda+1) ]/\lambda = \sum_{i=1}^{2r+1} f_i -1.
		      $$ Furthermore, the optimality of $\bar{x}$ implies that $ \sum_{i=1}^{2r+1} f_i \bar{x}_i = \sum_{i=1}^{2r+1} f_i -1  $. Thus we
can also check that $ g(\bar{x})=f_{2r+1}+3\lambda + 5 $.
		\item  [(d)]
		      $ \lambda+1 \le p\leq 2\lambda  $. On one hand, the feasibility of $\bar{x}$ requires
		      \begin{equation}
		      	\label{ineqa}
		      	\sum_{i=1}^{2r+1} f_i \highlight{\bar{x}_i}   \leq \lfloor \frac{\lambda \sum_{i=1}^{2r+1} f_i + \lambda^2  - p(\lambda +1)}{\lambda}\rfloor = \sum_{i=1}^{2r+1} f_i  + \lambda - p - 2,
		      \end{equation}
		      where the last equality follows from \highlight{$ \lambda + 1 \leq p \leq 2 \lambda$}.
		      On the other hand, the optimality of $\bar{x}$ requires that the inequality in \eqref{ineqa} holds with equality, yielding
		      $$ g(\bar{x}) = \sum_{i=1}^{2r+1} f_i + \lambda - p - 2+ p + 2\lambda + 6 -  \sum_{i=1}^{2r} f_i =f_{2r+1} + 3\lambda + 4.
		      $$
	\end{enumerate}

To summarize, the lifting coefficient $ \beta_n = f_{2r+1} + 3\lambda + 5  $ for
the problem of computing the \highlight{sequential} LCI constructed in the above,  if any only if $ p = \lambda  $ or $ p = \lambda -1 $, or equivalently,
the answer to the RPP instance is yes.
Since the RPP problem is $\mathcal{NP}$-hard and the construction of the companion problem is in polynomial,
we conclude that computing \highlight{a sequential} LCI is $\mathcal{NP}$-hard. This completes the proof.

	\qed
\end{proof}
\def}\def\href#1#2#3#4{\ifstrequal{#2}{[link]}{}{#2\newline}{}\def\href#1#2#3#4{\ifstrequal{#2}{[link]}{}{#2\newline}}
\bibliographystyle{spmpsci}