ode]haskellmathescape,linenos,numbersep=5pt,frame=lines,framesep=2mm,
epl]haskellmathescape,numbersep=5pt,frame=lines,framesep=2mm, 11institutetext:
DeepFlow, Inc., 3-16-40 Fujimi-shi Tsuruse nishi 354-0026, Japan

# Automatic Differentiation With Higher Infinitesimals, or Computational
Smooth Infinitesimal Analysis in Weil Algebra††thanks: This work was supported
by the Research Institute for Mathematical Sciences, an International Joint
Usage/Research Center located in Kyoto University.

Hiromi Ishii 11

###### Abstract

We propose an algorithm to compute the $C^{\infty}$-ring structure of
arbitrary Weil algebra. It allows us to do some analysis with _higher
infinitesimals_ numerically and symbolically. To that end, we first give a
brief description of the (Forward-mode) _automatic differentiation_ (AD) in
terms of _$C^{\infty}$ -rings_. The notion of a $C^{\infty}$-ring was
introduced by Lawvere [10] and used as the fundamental building block of
_smooth infinitesimal analysis_ and _synthetic differential geometry_ [11]. We
argue that interpreting AD in terms of $C^{\infty}$-rings gives us a unifying
theoretical framework and modular ways to express multivariate partial
derivatives. In particular, we can “package” higher-order Forward-mode AD as a
Weil algebra, and take tensor products to compose them to achieve multivariate
higher-order AD. The algorithms in the present paper can also be used for a
pedagogical purpose in learning and studying smooth infinitesimal analysis as
well.

###### Keywords:

automatic differentiation smooth infinitesimal analysis Weil algebrassmooth
algebras$C^{\infty}$-ringssymbolic-numeric algorithmssymbolic
differentiationGröbner basiszero-dimensional ideals

##  1 Introduction

_Automatic Differentiation_ (or, _AD_ for short) is a method to calculate
derivatives of (piecewise) smooth functions accurately and efficiently. AD has
a long history of research, and under the recent rise of differentiable
programming in machine learning, AD has been attracting more interests than
before recently.

_Smooth Infinitesimal Analysis_ (or, _SIA_ for short), on the other hand, is
an area of mathematics that uses _nilpotent infinitesimals_ to develop the
theory of real analysis. Its central building blocks are _Weil algebras_ ,
which can be viewed as the real line augmented with nilpotent infinitesimals.
Indeed, SIA is a subarea of _Synthetic Differential Geometry_ (SDG) initiated
by Lawvere [10], which studies smooth manifolds topos-theoretically, and
higher multivariate infinitesimals play crucial roles in building theory of,
e.g. vector fields, differential forms and tangent spaces. The key observation
of Lawvere is that manifolds can be classified solely by their smooth function
ring $C^{\infty}(M)$, and both such function rings and Weil algebras are
special cases of _$C^{\infty}$ -rings_.

It has been pointed out that AD and SIA have some connection; e.g. even
Wikipedia article [14] mentions the connection between first-order Forward-
mode AD with the ring $\mathbb{R}[X]/X^{2}$ of dual numbers. However, a
precise theoretical description of this correspondence is not well-
communicated, and further generalisation of AD in terms of SIA hasn’t been
discussed in depth.

The present paper aims at filling this gap, giving a unified description of AD
in terms of $C^{\infty}$-rings and Weil algebras. Furthermore, our main
contribution is algorithms to compute the $C^{\infty}$-ring structure of a
general Weil algebra. This enables automatic differentiation done in
_arbitrary_ Weil algebras other than dual numbers, and, together with tensor
products, lets us compute higher-order multivariate partial derivatives in a
modular and composable manner, packed as Weil algebra. Such algorithms can
also be used to learn and study the theory of SIA and SDG.

This paper is organised as follows. In Section 2, we review the basic concepts
and facts on $C^{\infty}$-rings and Weil algebras. This section provides basic
theoretical background — but the proofs of proposed algorithms are, however,
not directly dependent on the content of this section. So readers can skip
this section first and go back afterwards when necessary. Subsequently, we
discuss the connection between Forward-mode automatic differentiation and Weil
algebras in Section 3. There, we see how the notion of Weil algebra and
$C^{\infty}$-ring can be applied to treat higher-order partial ADs in a
unified and general setting. Then, in Section 4, we give algorithms to compute
the $C^{\infty}$-ring structure of an arbitrary Weil algebra. These algorithms
enable us to do _automatic differentiation with higher infinitesimals_ , or
_computational smooth infinitesimal analysis_. We give some small examples in
Section 5, using our proof-of-concept implementation [7] in Haskell. Finally,
we discuss related and possible future works and conclude in Section 6.

##  2 Preliminaries

In this section, we briefly review classical definitions and facts on Weil
algebras and $C^{\infty}$-rings without proofs, which will be used in Section
4. For theoretical detail, we refer readers to Moerdijk–Reyes [11, Chapters I
and II] or Joyce [8].

We use the following notational convention:

######  Definition 1 (Notation)

Throughout the paper, we use the following notation:

  * •

$g\circ f$ denotes the composite function from $A$ to $C$ of functions $f:A\to
B$ and $g:B\to C$, that is, the function defined by $(g\circ f)(x)=g(f(x))$
for all $x\in A$.

  * •

For functions $f_{i}:Z\to X_{i}\,(1\leq i\leq n)$,
$\braket{f_{1},\dots,f_{n}}$ denotes the product of functions $f_{i}$ given by
the universality of the product objects. That is, $\braket{f_{1},\dots,f_{n}}$
is the function of type $Z\to X_{1}\times\cdots\times X_{n}$ defined by
$\braket{f_{1},\dots,f_{n}}(z)=(f_{1}(z),\dots,f_{n}(z))\in
X_{1}\times\cdots\times X_{n}$ for all $z\in Z$

######  Definition 2 (Lawvere [10])

A _$C^{\infty}$_ -ring $A$ is a product-preserving functor from the category
$\mathord{\mathrm{CartSp}}$ of finite-dimensional Euclidean spaces and smooth
maps to the category $\mathord{\mathrm{Sets}}$ of sets.

We identify $A$ with $A(\mathbb{R})$ and $A^{n}$ with $A(\mathbb{R}^{n})$. For
a map $f:\mathbb{R}^{m}\to\mathbb{R}$, we call $A(f):A^{m}\to A$ the
_$C^{\infty}$ -lifting_ of $f$ to $A$.

Intuitively, a $C^{\infty}$-ring $A$ is an $\mathbb{R}$-algebra $A$ augmented
with $m$-ary operations $A(f):A^{m}\to A$ respecting composition, projection
and product for all smooth maps $f:\mathbb{R}^{m}\to\mathbb{R}$.

One typical example of a $C^{\infty}$-ring is a formal power series ring:

######  Theorem 2.1 (Implicitly in Lawvere [10]; See [11, 1.3 Borel’s
Theorem])

The ring $\mathbb{R}\llbracket X_{1},\dots,X_{n}\rrbracket$ of formal power
series with finitely many variables has the $C^{\infty}$-ring structure via
Taylor expansion at $\boldsymbol{0}$. In particular, lifting of a smooth map
$f:\mathbb{R}^{m}\to\mathbb{R}$ is given by:

$\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}(f)(g_{1},\dots,g_{m})=\sum_{\alpha\in\mathbb{N}^{n}}\frac{\boldsymbol{X}^{\alpha}}{\alpha!}D^{\alpha}(f\circ\braket{g_{1},\dots,g_{m}})(\boldsymbol{0}),$

where $\alpha!=\alpha_{1}!\dots\alpha_{n}!$ is the multi-index factorial and
$D^{\alpha}$ is the partial differential operator to degree $\alpha$.

The $C^{\infty}$-rings of central interest in this paper are _Weil algebras_ ,
and have a deep connection with
$\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}$:

######  Definition 3 (Weil algebra)

A _Weil algebra_ $W$ is an associative $\mathbb{R}$-algebra which can be
written as $W=\mathbb{R}[X_{1},\dots,X_{n}]/I$ for some ideal
$I\subseteq\mathbb{R}[\boldsymbol{X}]$ such that
$\braket{X_{1},\dots,X_{n}}^{k}\subseteq I$ for some $k\in\mathbb{N}$.

It follows that a Weil algebra $W$ is finite-dimensional as a
$\mathbb{R}$-linear space and hence $I$ is a _zero-dimensional_ ideal. A Weil
algebra $W$ can be regarded as a real line $\mathbb{R}$ augmented with
nilpotent infinitesimals $d_{i}={[X_{i}]}_{I}$. In what follows, we identify
an element $\boldsymbol{u}\in W$ of a $k$-dimensional Weil algebra $W$ with a
$k$-dimensional vector $\boldsymbol{u}=(u_{1},\dots,u_{k})\in\mathbb{R}^{k}$
of reals.

Although it is unclear from the definition, Weil algebras have the canonical
$C^{\infty}$-structure. First note that, if $I$ is zero-dimensional, we have
$\mathbb{R}[\boldsymbol{X}]/I\simeq\mathbb{R}\llbracket\boldsymbol{X}\rrbracket/I$.
Hence, in particular, any Weil algebra $W$ can also be regarded as a quotient
ring of the formal power series by zero-dimensional ideal. Thus, together with
Theorem 2.1, the following lemma shows that any Weil algebra $W$ has the
canonical $C^{\infty}$-ring structure:

######  Lemma 1 (Implicitly in Lawvere [10]; See [11, 1.2 Proposition])

For any $C^{\infty}$-ring $A$ and a ring-theoretical ideal $I\subseteq A$, the
quotient ring $A/I$ again has the canonical $C^{\infty}$-ring structure
induced by the canonical quotient mapping:

$(A/I)(f)([x_{1}]_{I},\dots,[x_{m}]_{I})\coloneqq\left[A(f)(x_{1},\dots,x_{m})\right]_{I},$

where $x_{i}\in A$ and $f:\mathbb{R}^{m}\xrightarrow{C^{\infty}}\mathbb{R}$.
In particular, the $C^{\infty}$-structure of Weil algebra $W$ is induced by
the canonical quotient mapping to that of
$\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}$.

##  3 Connection between Automatic Differentiation and Weil Algebras

In this section, based on the basic facts on $C^{\infty}$-rings and Weil
algebras reviewed in Section 2, we describe the connection of automatic
differentiation (AD) and Weil algebra.

_Forward-mode_ AD is a technique to compute a value and differential
coefficient of given univariate composition of smooth function efficiently. It
can be implemented by ad-hoc polymorphism (or equivalently, function
overloading). For detailed implementation, we refer readers to Elliott [2] and
Kmett’s ad package [9].

Briefly speaking, in Forward-mode AD, one stores both the value and
differential coefficient simultaneously, say in a form $f(x)+f^{\prime}(x)d$
for $d$ an indeterminate variable. Then, when evaluating composite functions,
one uses the Chain Rule for implementation:

$\frac{\mathrm{d}}{\mathrm{d}x}(g\circ f)(x)=f^{\prime}(x)g^{\prime}(f(x)).$

The following definitions of functions on dual numbers illustrate the idea:

$\displaystyle(a_{1}+b_{1}d)+(a_{2}+b_{2}d)$
$\displaystyle=(a_{1}+a_{2})+(b_{1}+b_{2})d$
$\displaystyle(a_{1}+b_{1}d)\times(a_{2}+b_{2}d)$
$\displaystyle=a_{1}a_{2}+(a_{1}b_{2}+a_{2}b_{1})d$
$\displaystyle\cos(a_{1}+b_{1}d)$
$\displaystyle=\cos(a_{1})-b_{1}\sin(a_{1})d$

The last equation for $\cos$ expresses the nontrivial part of Forward-mode AD.
As mentioned above, we regard $a_{1}+b_{1}d$ as a pair
$(a_{1},b_{1})=(f(x),f^{\prime}x)$ of value and differential coefficient of
some smooth function $f$ at some point $x$. So if
$a_{2}+b_{2}d=\cos(a_{1}+b_{1}d)$, we must have $a_{2}=\cos(f(x))=\cos a_{1}$
and $b_{2}=\frac{\mathrm{d}}{\mathrm{d}x}\cos(f(x))=-b_{1}\sin(a_{1})$ by
Chain Rule. The first two equations for addition and multiplication suggest us
to regard operations on Forward-mode AD as extending the algebraic structure
of $\mathbb{R}[d]=\mathbb{R}[X]/X^{2}$. Indeed, first-order Forward-mode AD
can be identified with the arithmetic on _dual numbers_ :

###### Definition 4

The _dual number ring_ is a Weil algebra $\mathbb{R}[X]/{X^{2}}$. We often
write $d={[X]}_{I}\in\mathbb{R}[d]$ and
$\mathbb{R}[d]\coloneq\mathbb{R}[X]/{X^{2}}$.

We use an analogous notation for multivariate versions:

$\mathbb{R}[d_{1},\dots,d_{k}]\coloneq\mathbb{R}[\boldsymbol{X}]/\braket{X_{1}^{2},\dots,X_{k}^{2}}.$

Since the dual number ring $\mathbb{R}[d]$ is a Weil algebra, one can apply
Theorem 2.1 and Lemma 1 to compute its $C^{\infty}$-structure. Letting
$f:\mathbb{R}\xrightarrow{C^{\infty}}\mathbb{R}$ be a univariate smooth
function, then we can derive the $C^{\infty}$-lifting
$\mathbb{R}[d](f):\mathbb{R}[d]\to\mathbb{R}[d]$ as follows:

$\displaystyle\mathbb{R}\llbracket{}X\rrbracket(f)(a+bX)$ $\displaystyle=\>$
$\displaystyle f(a)+\frac{\mathrm{d}}{\mathrm{d}x}(f(a+bx))(0)X+\cdots$
$\displaystyle(\text{by \lx@cref{creftypecap~refnum}{thm:series-is-smooth}})$
$\displaystyle=\>$ $\displaystyle f(a)+bf^{\prime}(a)X+\cdots$
$\displaystyle\xrightarrow{X\mapsto d}\>$ $\displaystyle
f(a)+bf^{\prime}(a)d,$ $\displaystyle\therefore$
$\displaystyle\mathbb{R}[d](f)(a+bd)=f(a)+bf^{\prime}(a)d.$
$\displaystyle(\text{by \lx@cref{creftypecap~refnum}{lem:quot-ring-ideal}})$
($\ast$)

One can notice that the derived $C^{\infty}$-structure in ($\ast$ ‣ 3) is
exactly the same as how to implement individual smooth functions for Forward-
mode AD. This describes the connection between Forward-mode AD and dual
numbers: Forward-mode AD is just a (partial) implementation of the
$C^{\infty}$-structure of the dual number ring $\mathbb{R}[d]$.

Let us see how this extends to higher-order cases. The most naïve way to
compute higher-order derivatives of smooth function is just to successively
differentiating it. This intuition can be expressed by duplicating the number
of the basis of dual numbers:

###### Theorem 3.1

For any $f:\mathbb{R}\xrightarrow{C^{\infty}}\mathbb{R}$ and
$\boldsymbol{x}\in\mathbb{R}^{n}$, we have:

$\mathbb{R}[d_{1},\dots,d_{k}](f)(x+d_{1}+\cdots+d_{n})=\sum_{0\leq i\leq
n}f^{(i)}(x)\sigma^{i}_{n}(\vec{d}),$

where, $\sigma^{i}_{k}(x_{1},\dots,x_{k})$ denotes the $k$-variate elementary
symmetric polynomial of degree $i$.

The above can be proven by an easy induction.

However, as one can easily see, terms in
$\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}/\braket{X_{i}^{2}}_{i}$
can grow exponentially and include duplicated coefficients. How could we
reduce such duplication and save space consumption? —this is where general
Weil algebras beyond (multivariate) dual numbers can play a role. We can get
derivatives in more succinct representation with _higher infinitesimal_ beyond
dual numbers:

###### Lemma 2

Let $I=\braket{X^{n+1}},W=\mathbb{R}[X]/I$ and $\varepsilon={[X]}_{I}$ for
$n\in\mathbb{N}$. Given $f:\mathbb{R}\xrightarrow{C^{\infty}}\mathbb{R}$ and
$a\in\mathbb{R}$, we have:

$W(f)(a+\varepsilon)=\sum_{k\leq n}\frac{f^{(k)}(a)}{k!}\varepsilon^{k}.$

In this representation, we have only $(n+1)$-terms, and hence it results in
succinct and efficient representation of derivatives.

If we duplicate such higher-order infinitesimals as much as needed, one can
likewise compute _multivariate_ higher-order derivatives all at once, up to
some multidegree $\beta$:

###### Lemma 3

Let $I=\Braket{X_{i}^{\beta_{i}+1}\ }{\ i\leq m}$,
$W=\mathbb{R}[X_{1},\dots,X_{m}]/I$, and $\varepsilon_{i}={[X_{i}]}_{I}$ for
$\beta={(\beta_{i})}_{i\leq m}\in\mathbb{N}^{m}$. For
$f:\mathbb{R}^{m}\xrightarrow{C^{\infty}}\mathbb{R}$ and
$\boldsymbol{a}={(a_{i})}_{i\leq m}\in\mathbb{R}^{m}$, we have:

$W(f)(a_{1}+\varepsilon_{1},\dots,a_{m}+\varepsilon_{m})=\sum_{\delta_{i}\leq\beta_{i}}\frac{D^{\delta}f}{\delta!}(\boldsymbol{a})\
\varepsilon_{1}^{\delta_{1}}\cdots\varepsilon_{m}^{\delta_{m}}.$

Note that the formal power series ring
$\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}$ can be viewed as the
inverse limit of
$\mathbb{R}[\boldsymbol{X}]/\braket{\boldsymbol{X}^{\beta}}$’s. In other
words, if we take a limit $\beta_{i}\to\infty$, we can compute any higher
derivative up to any finite orders; this is exactly what _Tower-mode_ AD aims
at, modulo factor $\frac{1}{\beta!}$.

In this way, we can view AD as a technique to compute higher derivatives
simultaneously by partially implementing a certain $C^{\infty}$-ring111Such
implementation is inherently a partial approximation: there are
$2^{\aleph_{0}}$-many smooth functions, but there are only countably many
computable (floating) functions.. Forward-mode AD (of first-order) computes
the $C^{\infty}$-structure of the dual number ring $\mathbb{R}[d]$; Tower-mode
AD computes that of the formal power series ring
$\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}$ (modulo reciprocal
factorial).

So far, we have used a Weil algebra of form
$\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}/I$. So, do we need to
define new ideals by hand whenever one wants to treat multiple variables? The
answer is no:

######  Lemma 4 (See [11, 4.19 Corollary])

For ideals $I\subseteq\mathbb{R}\llbracket\boldsymbol{X}\rrbracket$ and
$J\subseteq\mathbb{R}\llbracket{}\boldsymbol{Y}\rrbracket$, we have:

$\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}/I\otimes_{\mathbb{R}}\mathbb{R}\llbracket{}\boldsymbol{Y}\rrbracket/J\simeq\mathbb{R}\llbracket{}\boldsymbol{X},\boldsymbol{Y}\rrbracket/(I,J),$

where $\otimes_{\mathbb{R}}$ is a tensor product of $C^{\infty}$-rings.

Thanks to this lemma, we don’t have to define $I$ by hand every time, but can
take tensor products to compose predefined Weil algebras to compute
multivariate and higher-order derivatives. Examples of such calculations will
be presented in Section 5.

##  4 Algorithms

In this section, we will present the main results of this paper: concrete
algorithms to compute the $C^{\infty}$-structure of arbitrary Weil algebra and
their tensor products. For examples of applications of the algorithm presented
here, the reader can skip to the next section 5 to comprehend the actual use
case.

###  4.1 Computing $C^{\infty}$-structure of Weil algebra

Let us start with algorithms to compute the $C^{\infty}$-structure of a
general Weil algebra. Roughly speaking, the algorithm is threefold:

  1. 1.

A procedure deciding Weil-ness of an ideal and returning data required to
compute the $C^{\infty}$-structure (WeilTest, Algorithm 1),

  2. 2.

A procedure to compute the lifting $W(f):W^{m}\to W$ to a Weil algebra $W$
from $\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}(f)$ (LiftWeil,
Algorithm 2), and

  3. 3.

A procedure to lift smooth map $f:\mathbb{R}^{m}\to\mathbb{R}$ to the
$n$-variate formal power series ring
$\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}$ (LiftSeries,
Algorithm 3).

We start with Weil-ness testing. First, we define the basic data needed to
compute the $C^{\infty}$-structure of Weil algebras:

######  Definition 5 (Weil settings)

The _Weil setting_ of a Weil algebra $W$ consists of the following data:

  1. (1).

Monomial basis $\set{\boldsymbol{b}_{1},\dots,\boldsymbol{b}_{\ell}}$ of $W$,

  2. (2).

$M$, the multiplication table of $W$ in terms of the basis,

  3. (3).

$(k_{1},\dots,k_{n})\in\mathbb{N}^{n}$ such that $k_{i}$ is the maximum
satisfying $X_{i}^{k_{i}}\notin I$ for each $i$, and

  4. (4).

$\operatorname{NV}_{W}$, a table of representations of non-vanishing monomials
in $W$; i.e. for any $\alpha=(\alpha_{1},\dots,\alpha_{n})\in\mathbb{N}^{n}$,
if $\alpha_{i}\leq k_{i}$ for all $i$, then
$\operatorname{NV}_{W}(\boldsymbol{X}^{\alpha})=(c_{1},\dots,c_{n})\in\mathbb{R}^{k}$
satisfies $[\boldsymbol{X}^{\alpha}]_{I}=\sum_{i}c_{i}\boldsymbol{b}_{i}$.

A basis and multiplication table allow us to calculate the ordinary
$\mathbb{R}$-algebra structure of Weil algebra $W$. The latter two data,
$\vec{k}$ and $\operatorname{NV}_{W}$, are essential in computing the
$C^{\infty}$-structure of $W$. In theory, (4) is unnecessary if one stores a
Gröbner basis of $I$; but since normal form calculation modulo $G$ can take
much time in some case, we don’t store $G$ itself and use the precalculated
data $\operatorname{NV}$. It might be desirable to calculate
$\operatorname{NV}_{W}$ as lazily as possible. Since it involves Gröbner basis
computation it is more desirable to delay it as much as possible and do in an
on-demand manner.

With this definition, one can decide Weilness and compute their settings:

######  Algorithm 1 (WeilTest)

Input

An ideal $I\subseteq\mathbb{R}[X_{1},\dots,X_{n}]$

Output

Returns the Weil settings of $W=\mathbb{R}[\boldsymbol{X}]/I$ if it is a Weil
algebra; otherwise `No`.

Procedure

WeilTest

⬇

1$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}I$f $I$ is
not zero-dimensional

2 Return No

3$\set{\boldsymbol{b}_{1},\dots,\boldsymbol{b}_{\ell}}$
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\leftarrow$
Monomial basis of $I$

4$M$
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\leftarrow$
the Multiplication table of $W$

5For i in 1..n

6 $p_{i}$
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\leftarrow$
the monic generator of $I\cap\mathbb{R}[X_{i}]$

7 $\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}I$f
$p_{i}$ is not a monomial

8 Return No

9 $k_{i}$
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\leftarrow$
$\deg(p_{i})-1$

10$\operatorname{NV}_{W}$
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\leftarrow$
{}

11For $\alpha$ in $\Set{\alpha\in\mathbb{N}^{n}}{\alpha_{i}\leq k_{i}\;\forall
i\leq\ell}$

12 $c_{1}\boldsymbol{b}_{1}+\cdots+c_{\ell}\boldsymbol{b}_{\ell}$
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\leftarrow$
${\overline{\boldsymbol{X}^{\alpha}}}^{G}$

13 $\operatorname{NV}_{W}(\boldsymbol{X}^{\alpha})$
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\leftarrow$
($c_{1},\dots,c_{\ell}$)

14Return ($\vec{\boldsymbol{b}},M,\vec{k},\operatorname{NV}_{W}$)

###### Theorem 4.1

Algorithm 1 terminates and returns expected values.

###### Proof

Algorithms to decide the zero-dimensionality and calculate their
multiplication table is well-known (for details, we refer readers to
Cox–Little–O’Shea [1, Chapter 2]). So the only non-trivial part is nilpotence
detection (Lines 5 to 9). But, again, this is just a variation of radical
calculation algorithm for zero-dimensional ideals. Indeed, since each
$\mathbb{R}[X_{i}]$ is a PID, we have $X_{i}^{k}\in I\cap R[X_{i}]$ iff
$p_{i}\mathrel{\mid}X_{i}^{k}$, hence $p_{i}$ is a monomial iff $X_{i}$ is
nilpotent in $W$.

Now that we have the basis and multiplication table at hand, we can calculate
the ordinary algebraic operations just by the standard means.

With upper bounds $\vec{k}$ of powers and representations
$\operatorname{NV}_{W}$ of non-vanishing monomials, we can now compute the
$C^{\infty}$-structure of an arbitrary Weil algebra, when given a lifting of
smooth mapping $f$ to
$\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}$:

######  Algorithm 2 (LiftWeil)

Input

$I\subseteq\mathbb{R}[\boldsymbol{X}]$, an ideal where
$W=\mathbb{R}[\boldsymbol{X}]/I$ is a Weil algebra,
$\mathbb{R}\llbracket\boldsymbol{X}\rrbracket(f):\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}^{m}\to\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}$,
a lifting of a smooth map $f:\mathbb{R}^{m}\to\mathbb{R}$ to
$\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}$, and
$\vec{\boldsymbol{u}}=(\boldsymbol{u}_{1},\dots,\boldsymbol{u}_{m})\in
W^{m}$,.

Output

$\boldsymbol{v}=W(f)(\vec{\boldsymbol{u}})\in W$, the value of $f$ at
$\vec{\boldsymbol{u}}$ given by $C^{\infty}$-structure.

Procedure

LiftWeil

⬇

1$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}g_{i}$
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\leftarrow$
$(\boldsymbol{b}_{1},\dots,\boldsymbol{b}_{k})\cdot\boldsymbol{u}_{i}\in\mathbb{R}[\boldsymbol{X}]$
for i $\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\leq$
m

2h = $\sum_{\alpha}c_{\alpha}\boldsymbol{X}^{\alpha}$
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\leftarrow$
$\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}(f)(\vec{g})$

3$\boldsymbol{v}$
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\leftarrow$
0

4For $\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\alpha$
with $\alpha_{i}\leq k_{i}\,\forall i$

5 $\boldsymbol{v}$
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\leftarrow$
$\boldsymbol{v}$
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}+$
$c_{\alpha}\operatorname{NV}_{W}(\boldsymbol{X}^{\alpha})$

6Return $\boldsymbol{v}$

The termination and validity of Algorithm 2 are clear. One might feel it
problematic that Algorithm 2 requires _functions_ as its input. This can be
_any_ smooth computable functions on the coefficient type. Practically, we
expect a composite function of standard smooth floating-point functions as its
argument, for example, it can be $x\mapsto\sin(x)$, $(x,y)\mapsto e^{\sin
x}y^{2}$, and so on. In the modern programming language – like Haskell, Rust,
LISP, Ruby, etc. – one don’t need to worry about their representation, as we
can already freely write _higher-order functions_ that take functions or
closures as their arguments. Even in the low-level languages such as C/C++,
one can use function pointers or whatever to pass an arbitrary function to
another function.

Now that we can compute the $\mathbb{R}$-algebraic and $C^{\infty}$-structure
of a Weil algebra solely from its Weil setting, one can hard-code pre-
calculated Weil settings for known typical Weil algebras, such as the dual
number ring or higher infinitesimal rings of the form
$\mathbb{R}[X]/(X^{n+1})$, to reduce computational overheads.

####  4.1.1 Computing the $C^{\infty}$-structure of
$\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}$

So it remains to compute the $C^{\infty}$-structure of
$\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}$. Thanks to Theorem
2.1, we know the precise definition of $C^{\infty}$-lifting to
$\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}$:

$\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}(f)(g_{1},\dots,g_{m})=\sum_{\alpha\in\mathbb{N}^{n}}\frac{\boldsymbol{X}^{\alpha}}{\alpha!}D^{\alpha}(f\circ\braket{g_{1},\dots,g_{n}})(\boldsymbol{0}).$

As noted in Section 3, as a $C^{\infty}$-ring, the formal power series ring is
isomorphic to multivariate Tower-mode AD. It can be implemented in various
ways, such as Lazy Multivariate Tower AD [13], or nested Sparse Tower AD [9,
module Numeric.AD.Rank1.Sparse]. For reference, we include a succinct and
efficient variant mixing these two techniques in Appendix 0.A.

Both Tower-mode AD and formal power series can be represented as a formal
power series. The difference is the interpretation of coefficients in a given
series. On one hand, a coefficient of $\boldsymbol{X}^{\alpha}$ in Tower AD is
interpreted as the $\alpha$th partial differential coefficient
$D^{\alpha}f(\boldsymbol{a})$, where
$\boldsymbol{a}=(g_{1}(0),\dots,g_{m}(0))$. On the other hand, in
$\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}$ it is interpreted as
$\frac{D^{\alpha}f(\boldsymbol{a})}{\alpha!}$. To avoid the confusion, we
adopt the following convention: Tower-mode AD is represented as a function
from monomials $\boldsymbol{X}^{\alpha}$ to coefficient $\mathbb{R}$ in what
follows, whilst $\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}$ as-
is. Note that this is purely for notational and descriptional distinctions,
and does not indicate any essential difference.

With this distinction, we use the following notation and transformation:

###### Definition 6

$\mathord{\mathrm{Tower}}=\Set{f}{f:\mathbb{N}^{n}\to\mathbb{R}}$ denotes the
set of all elements of Tower-mode AD algebra. We denote $C^{\infty}$-lifting
of $f:\mathbb{R}^{m}\to\mathbb{R}$ to $\mathord{\mathrm{Tower}}$ by
$\mathord{\mathrm{Tower}}(f):\mathord{\mathrm{Tower}}^{m}\to\mathord{\mathrm{Tower}}$.

A _reciprocal factorial transformation_
$\operatorname{RF}:\mathord{\mathrm{Tower}}\to\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}$
is defined as follows:

$\operatorname{RF}\left(f\right)=\sum_{\alpha\in\mathbb{N}^{n}}\frac{f(\alpha)}{\alpha!}\boldsymbol{X}^{\alpha}.$

Then, the inverse reciprocal factorial transformation is given by:

$\operatorname{RF}^{-1}\left(\sum_{\alpha\in\mathbb{N}^{n}}c_{\alpha}\boldsymbol{X}^{\alpha}\right)=\lambda(\boldsymbol{X}^{\alpha}).\
\alpha!\cdot c_{\alpha}.$

######  Algorithm 3 (LiftSeries)

Input

$f:\mathbb{R}^{m}\xrightarrow{C^{\infty}}\mathbb{R}$, a smooth function which
admits Tower AD,
$g_{1},\dots,g_{n}\in\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}$,
formal power series.

Output

$\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}(f)(g_{1},\dots,g_{m})\in\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}$,
$C^{\infty}$-lifting to the formal power series ring.

Procedure

LiftSeries

⬇

1$\hat{f}$
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\leftarrow$
$\mathord{\mathrm{Tower}}(f)(\hat{g}_{1},\dots,\hat{g}_{m})$

2Return $\operatorname{RF}(\hat{f})$

###  4.2 Tensor product of Weil algebras

As indicated by Lemma 4, tensor products enable us to compose multiple Weil
algebras into one and use them to compute higher-order multivariate
derivatives. Here, we give a simple procedure to compute Weil settings of the
tensor product.

######  Algorithm 4 (WeilTensor)

Input

Weil settings of two Weil algebras $W_{1},W_{2}$, with
$\set{\boldsymbol{b}_{1}^{i},\dots,\boldsymbol{b}_{\ell_{i}}^{i}}$ a basis,
$(k^{i}_{1},\dots,k^{i}_{n_{i}})$ an upper bounds and $M_{i}$ a multiplication
table for each $W_{i}$.

Output

Weil settings of $W_{1}\otimes_{\mathbb{R}}W_{2}$.

Procedure

WeilTensor

⬇

1M
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\leftarrow$
{}

2For ({$\boldsymbol{b}^{1}_{L},\boldsymbol{b}^{1}_{R}$},
$(c_{1},\dots,c_{\ell_{1}})$) in $M_{1}$

3 For ({$\boldsymbol{b}^{2}_{L},\boldsymbol{b}^{2}_{R}$},
$(d_{1},\dots,d_{\ell_{1}})$) in $M_{2}$

4 M({$\boldsymbol{b}^{1}_{L}\boldsymbol{b}^{2}_{L}$,
$\boldsymbol{b}^{1}_{R}\boldsymbol{b}^{2}_{R}$})
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\leftarrow$
Convol($\vec{c}$, $\vec{d}$)

5$\operatorname{NV}_{W_{1}\otimes W_{2}}$
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\leftarrow$
{}

6For ($\boldsymbol{X}^{\alpha}$, $(c_{1},\dots,c_{\ell_{1}})$) in
$\operatorname{NV}_{W_{1}}$

7 For ($\boldsymbol{Y}^{\beta}$, $(d_{1},\dots,d_{\ell_{2}})$) in
$\operatorname{NV}_{W_{2}}$

8 $\operatorname{NV}_{W_{1}\otimes
W_{2}}(\boldsymbol{X}^{\alpha}\boldsymbol{Y}^{\beta})$
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\leftarrow$
Convol($\vec{c}$, $\vec{d}$)

9Return
$(\boldsymbol{b},M,(\vec{k}^{1},\vec{k}^{2}),\operatorname{NV}_{W_{1}\otimes
W_{2}})$

Here, Convol is a convolution of two sequences:

Procedure

$\mathop{\text{Convol}}((c_{1},\dots,c_{\ell_{1}}),(d_{1},\dots,d_{\ell_{2}}))$

⬇

1 j
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\leftarrow$
$\lfloor\frac{i}{\ell_{2}}\rfloor$; k
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\leftarrow$
$i\mod{\ell_{2}}$

2 $a_{i}$
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\leftarrow$
$c_{j}d_{k}$

3Return ($a_{1},\dots,a_{\ell_{1}\ell_{2}}$)

The validity proof is routine work.

##  5 Examples

We have implemented the algorithms introduced in the previous section on top
of two libraries: computational-algebra package [4, 6] and ad package [9]. The
code is available on GitHub [7].

###  5.1 Higher-order derivatives via dual numbers and higher infinitesimals

As indicated by Theorem 3.1 and Lemma 2, to compute higher-order derivatives
of univariate functions, we can use tensor products of Dual numbers or higher-
order infinitesimals.

Let us first compute higher-order derivatives of $\sin(x)$ up to $3$. First,
Let us use a tensor product of dual numbers: d1, d2, d3 :: Floating a =¿ Weil
(D1 —*— D1 —*— D1) a [d1, d2, d3] = map di [0..] Here, Weil w a represents the
type of Weil algebra with its setting given in w, D1 the dual number ideal
$I=(X^{2})$, and
|$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\times$|
the tensor product operator. Each $d_{i}$ corresponds to $i$th infinitesimal.

Next, we calculate higher-order differential coefficients at $x=\frac{\pi}{6}$
up to the third order:

¿¿¿ (sin (pi/6), cos (pi/6), -sin (pi/6), -cos (pi/6)) ( 0.49999999999999994,
0.8660254037844387, -0.49999999999999994, -0.8660254037844387)

¿¿¿ sin (pi/6 + d0 + d1 + d2) -0.8660254037844387 d(0) d(1) d(2) -
0.49999999999999994 d(0) d(1) \- 0.49999999999999994 d(0) d(2) -
0.49999999999999994 d(1) d(2) \+ 0.8660254037844387 d(0) + 0.8660254037844387
d(1) \+ 0.8660254037844387 d(2) + 0.49999999999999994

It is easy to see that terms of degree $i$ have the coefficients
$\sin^{(i)}(\pi/6)$. Since our implementation is polymorphic, if we apply the
same function to the type for symbolic computation, say Symbolic, we can
reconstruct symbolic differentiation and check that the result is indeed
correct symbolically:

¿¿¿ x :: Weil w Symbolic ¿¿¿ x = injectCoeff (var ”x”)

¿¿¿ normalise
¡$>sin(x+d0+d1+d2)((-1.0)*cosx)d(0)d(1)d(2)+(-(sinx))d(0)d(1)+(-(sinx))d(0)d(2)+(-(sinx))d(1)d(2)+(cosx)d(0)+(cosx)d(1)+(cosx)d(2)+(sinx)\par
Asstatedbefore,thetensor-of-
dualsapproachblowsthenumberoftermsexponentially.Letusseehowhigherinfinitesimalworks.\par\code
eps::Floatinga=>Weil(DOrder4)aeps=di0\par
Here,\definecolor{tcbcolback}{rgb}{.75,.75,.75}\definecolor{tcbcolframe}{rgb}{.5,.5,.5}\definecolor{tcbcol@origin}{rgb}{0,0,0}\definecolor{.}{rgb}{0,0,0}\definecolor{.}{rgb}{0,0,0}\leavevmode\hbox
to269.12pt{\vbox to16.4pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower
0.0pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{{}{}{}{}\pgfsys@beginscope\pgfsys@invoke{
}{}{}{}{}{}{}{}{}\definecolor[named]{pgffillcolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@fill{.5}\pgfsys@invoke{
}\pgfsys@fill@opacity{1.0}\pgfsys@invoke{
}{{}{}{{}}}{{}{}{{}}}{}{}{{}{}{{}}}{{}{}{{}}}{}{}{{}{}{{}}}{{}{}{{}}}{}{}{{}{}{{}}}{{}{}{{}}}{}{}\pgfsys@moveto{0.0pt}{3.69884pt}\pgfsys@lineto{0.0pt}{12.69884pt}\pgfsys@curveto{0.0pt}{14.74167pt}{1.656pt}{16.39767pt}{3.69884pt}{16.39767pt}\pgfsys@lineto{265.41759pt}{16.39767pt}\pgfsys@curveto{267.46042pt}{16.39767pt}{269.11642pt}{14.74167pt}{269.11642pt}{12.69884pt}\pgfsys@lineto{269.11642pt}{3.69884pt}\pgfsys@curveto{269.11642pt}{1.656pt}{267.46042pt}{0.0pt}{265.41759pt}{0.0pt}\pgfsys@lineto{3.69884pt}{0.0pt}\pgfsys@curveto{1.656pt}{0.0pt}{0.0pt}{1.656pt}{0.0pt}{3.69884pt}\pgfsys@closepath\pgfsys@fill\pgfsys@invoke{
}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\pgfsys@beginscope\pgfsys@invoke{
}{}{}{}{}{}{}{}{}\definecolor[named]{pgffillcolor}{rgb}{.75,.75,.75}\pgfsys@color@gray@fill{.75}\pgfsys@invoke{
}\pgfsys@fill@opacity{1.0}\pgfsys@invoke{
}{{}{}{{}}}{{}{}{{}}}{}{}{{}{}{{}}}{{}{}{{}}}{}{}{{}{}{{}}}{{}{}{{}}}{}{}{{}{}{{}}}{{}{}{{}}}{}{}\pgfsys@moveto{0.85358pt}{3.69884pt}\pgfsys@lineto{0.85358pt}{12.69884pt}\pgfsys@curveto{0.85358pt}{14.27025pt}{2.12743pt}{15.5441pt}{3.69884pt}{15.5441pt}\pgfsys@lineto{265.41759pt}{15.5441pt}\pgfsys@curveto{266.989pt}{15.5441pt}{268.26285pt}{14.27025pt}{268.26285pt}{12.69884pt}\pgfsys@lineto{268.26285pt}{3.69884pt}\pgfsys@curveto{268.26285pt}{2.12743pt}{266.989pt}{0.85358pt}{265.41759pt}{0.85358pt}\pgfsys@lineto{3.69884pt}{0.85358pt}\pgfsys@curveto{2.12743pt}{0.85358pt}{0.85358pt}{2.12743pt}{0.85358pt}{3.69884pt}\pgfsys@closepath\pgfsys@fill\pgfsys@invoke{
}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@fill@opacity{1.0}\pgfsys@invoke{
}{{{}}{{}}{{}}{{}}{{}}{{}}{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{6.5441pt}{5.69884pt}\pgfsys@invoke{
}\hbox{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hbox{\set@color{\ignorespaces\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor\lst@@@set@language{\@listingGroup{ltx_lst_identifier}{DOrder}}{\@listingGroup{ltx_lst_space}{
}}{\@listingGroup{ltx_lst_identifier}{n}}}}}}}}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}correspondstoanalgebra$R[X]/(X^n)$.Notethat,accordingto\lx@cref{creftypecap~refnum}{lem:higher-
infinitesimal},tocalculatean$n$\textsuperscript{th}derivativewehavetouse$R[X]/(X^n+1)$.\par\repl>>>(sin(pi/6),cos(pi/6),-sin(pi/6)/2,-cos(pi/6)/6)(0.49999999999999994,0.8660254037844387,-0.24999999999999997,-0.14433756729740646)\par>>>sin(pi/6+eps)-0.14433756729740646d(0)^{3}-0.24999999999999997d(0)^{2}+0.8660254037844387d(0)+0.49999999999999994\par>>>normalise<$>sin(x+eps)((-1.0)*cosx/6.0)d(0)^{3}+((-(sinx))/2.0)d(0)^{2}+(cosx)d(0)+(sinx)\par
Notethatby\lx@cref{creftypecap~refnum}{lem:higher-
infinitesimal},eachcoefficientisnotdirectlyadifferentialcoefficient,butdividedby$k!$,thatis$f(x
+ ε) = ∑_k ≤3 f(k)(x)k ! ε^k$.\par
LetusseehowtensorproductsofhigherWeilalgebrascanbeusedtomultivariatehigher-
orderpartialderivatives.Supposewewanttocalculatepartialderivativesof$f(x, y) =
e^2x siny$upto$(2, 1)$\textsuperscript{th}order.\repl
eps1,eps2::Weil(DOrder3|*|DOrder2)a(eps1,eps2)=(di0,di1)\par>>>f(2+eps1)(pi/6+eps2)94.5667698566742d(0)^{2}d(1)+54.59815003314423d(0)^{2}+94.5667698566742d(0)d(1)+54.59815003314423d(0)+47.2833849283371d(1)+27.299075016572115\par>>>normalise<$>f(x+eps1)(y+eps2)(4.0*exp(2.0*x)/2.0*cosy)d(0)^{2}d(1)+(4.0*exp(2.0*x)/2.0*siny)d(0)^{2}+(2.0*exp(2.0*x)*cosy)d(0)d(1)+(2.0*exp(2.0*x)*siny)d(0)+(exp(2.0*x)*cosy)d(1)+(exp(2.0*x)*siny)Onecanseethatthecoefficientof$d(0)^i
d(1)^j$correspondsexactlytothevalue$D^(i,j)f(x,y)/i!j!$.Inthisway,wecanfreelycomposemultipleWeilalgebratocalculatevariouspartialderivativesmodularly.\par$

###  5.2 Computation in General Weil Algebra

All examples so far were about the predefined, specific form of a Weil
algebra. Here, we demonstrate that we can determine whether the given ideal
defines Weil algebras with Algorithm 1, and do some actual calculation in
arbitrary Weil algebra.

First, we see that WeilTest rejects invalid ideals:

– R[X,Y]/(X^3 - Y), not zero-dimensional ¿¿¿ isWeil (toIdeal [x ^3 - y ::
Q[x,y]]) Nothing

– R[X]/(X^2 - 1), which is zero-dimensional but not Weil ¿¿¿ isWeil (toIdeal
[x ^2 - 1 :: Q[x]]) Nothing

Next, we try to calculate in arbitrary chosen Weil algebra,
$W=\mathbb{R}[x,y]/(x^{2}-y^{3},y^{4})$, whose corresponding meaning in AD is
unclear but is a Weil algebra as a matter of fact.

i :: Ideal (Rational[x,y]) i = toIdeal [x ^2 - y ^3, y ^4]

¿¿¿ isWeil i Just WeilSettings weilBasis =[[0,0],[0,1], …, [3,0]] ,
nonZeroVarMaxPowers = [3,3] , weilMonomDic = [([0,2],[0,0,0,1,0,0,0,0]), …,
([1,3],[0,0,0,0,0,0,0,1])] , table = [((0,0),1),((1,3),d(0)^2), …,
((3,4),d(0)^3)]

Let us see what will happen evaluating $\sin(a+d_{0}+d_{1})$, where
$d_{0}=[x]_{I},d_{1}=[y]_{I}$?

¿¿¿ withWeil i (sin (pi/4 + di 0 + di 1)) -2.7755575615628914e-17 d(0)^3 - … +
0.7071067811865476 d(0) \+ 0.7071067811865476 d(1) + 0.7071067811865475

¿¿¿ withWeil i (normalise
¡$>sin(x+di0+di1))((-1.0)*(-(sinx))/6.0+(-1.0)*cosx/6.0)d(0)^{3}+...+(cosx)d(0)+(cosx)d(1)+(sinx)\par
Carefullyanalysingeachoutput,onecanseethattheoutputcoincideswithwhatisgivenby\lx@cref{creftypecap~refnum}{thm:series-
is-smooth}and\lx@cref{creftypecap~refnum}{lem:quot-ring-ideal}.\par$

##  6 Discussions and Conclusions

We have illustrated the connection between automatic differentiation (AD) and
$C^{\infty}$-rings, especially Weil algebras. Methods of AD can be viewed as
techniques to calculate partial coefficients simultaneously by partially
implementing the $C^{\infty}$-lifting operator for certain $C^{\infty}$-ring.
Especially, Forward-mode AD computes the $C^{\infty}$-structure of the dual
number ring $\mathbb{R}[d]=\mathbb{R}[X]/X^{2}$ and Tower-mode computes that
of the formal power series ring
$\mathord{\mathbb{R}\llbracket\boldsymbol{X}\rrbracket}$.

The dual number ring $\mathbb{R}[d]$ is an archetypical example of Weil
algebra, which formalises the real line with nilpotent infinitesimals. We
generalised this view to arbitrary Weil algebras beyond dual numbers, enabling
us to compute higher-order derivatives efficiently and succinctly. We gave
general algorithms to compute the $C^{\infty}$-structure of Weil algebras.
With tensor products, one can easily compose (univariate) higher-order AD
corresponding to Weil algebras into multivariate ones.

In this last section, we briefly discuss the possible applications other than
AD, related works and future works.

###  6.1 Possible Applications and Related Works

Beside the reformulation of AD, we can argue that our methods can be used for
a pedagogical purpose in teaching _Smooth Infinitesimal Analysis_ (SIA) and
_Synthetic Differential Geometry_ (SDG). In those fields, arguing in the
appropriate intuitionistic topos, various infinitesimal spaces corresponding
to Weil algebra is used to build a theory, expressed by the following
_generalised Kock-Lawvere axiom_ [11]:

> For any Weil algebra $W$, the following evaluation map gives an isomorphism:
>
> $\displaystyle\mathop{\mathrm{ev}}:$ $\displaystyle
> W\to\mathbb{R}^{\mathop{\mathrm{Spec}}_{\mathbb{R}}{W}}$ $\displaystyle
> a\mapsto\lambda f.f(a)$

This is another way to state the fact that Weil algebras are
$C^{\infty}$-rings, viewed within some topoi. For dual numbers, their meaning
is clear: it just couples a value and their (first-order) differential
coefficient. However, solely from Kock-Lawvere axiom, it is unclear what the
result is in general cases. With the algorithms we have proposed, students can
use computers to calculate the map given by the axiom. In SIA and SDG, there
are plenty of uses of generalised infinitesimal spaces such as
$\mathbb{R}[x_{1},\dots,x_{n}]/\braket{x_{i}x_{j}}{i,j\leq n}$ or
$\mathbb{R}[x]/(x^{n})$. Hence, concrete examples for these Weil algebras can
help to understand the theory.

In the context of SDG, applying techniques in computer algebra to Weil
algebras has attained only little interest. One such example is
Nishimura–Osoekawa [12]: they apply zero-dimensional ideal algorithms to
compute the generating relation of limits of Weil algebras. However, their
purpose is to use computer algebra to ease heavy calculations needed to
develop the theory of SDG, and hence they are not interested in computing the
$C^{\infty}$-structure of Weil algebras.

Implementing AD in a functional setting has a long history. See, for example,
Elliott [2] for explanation and ad package by Kmett [9] for actual
implementation. In ad package, so-called _Skolem trick_ , or _RankN trick_ is
applied to distinguish multiple directional derivatives. We argue that our
method pursues other direction of formulations; we treat higher infinitesimals
as first-class citizens, enabling us to treat higher-order AD in a more
modular and composable manner.

###  6.2 Future Works

In SDG, $C^{\infty}$-ring and higher infinitesimals are used as fundamental
building blocks to formulate manifolds, vector fields, differential forms, and
so on. Hence, if one can extend our method to treat a general
$C^{\infty}$-ring $C^{\infty}(M)$ of real-valued smooth functions on $M$, it
can open up a new door to formulate differential geometry on a computer. With
such a formulation, we can define differential-geometric objects in more
synthetic manner using nilpotent infinitesimals – for example, one can define
the tangent space $T_{x}M$ at $x\in M$ on some manifold $M$ to be the
collection of $f:D\to M$ with $f(0)=x$, where $D$ is the set of nilpotents of
order two. Another virtue of such system is that we can treat infinitesimal
spaces (derived from Weil algebras), manifolds, functions spaces, and vector
spaces uniformly – they are all living in the same category. See
Moerdijk–Reyes [11] for more theoretical details. One major obstacle in this
direction is that, even if $C^{\infty}(M)$ is finitely _presented_ as a
$C^{\infty}$-ring, it is NOT finitely _generated_ as an $\mathbb{R}$-algebra,
but $2^{\aleph_{0}}$-generated, by its very nature. Hence, it seems impossible
to compute $C^{\infty}(M)$ in purely symbolic and direct way; we need some
workarounds or distinct formulations to overcome such obstacles.

As for connections with AD, there is also plenty of room for further
exploration. There are so many “modes” other than Forward- and Tower-modes in
AD: for examples, Reverse, Mixed, etc. amongst others. From the point of view
of Weil algebras, they are just implementation details. But such details
matter much when one takes the efficiency seriously. It might be desirable to
extend our formulation to handle such differences in implementation method.
For such direction, Elliot [3] proposes the categorical formulation. Exploring
how that approach fits with our algebraic framework could be interesting
future work, and perhaps also shed a light on the way to realise the
aforementioned computational SDG.

## Acknowledgments

The author is grateful to Prof. Akira Terui, for encouraging to write this
paper and many helpful feedbacks, and to anonymous reviewers for giving
constructive comments.

##  Appendix 0.A Succinct Multivariate Lazy Tower AD

For completeness, we include the referential implementation of the Tower-mode
AD in Haskell, which can be used in Algorithm 2. The method presented here is
a mixture of Lazy Multivariate Tower [13] and nested Sparse Tower [9]. For
details, we refer readers to the related paper by the author in RIMS Kôkyûroku
[5].

The idea is simple: we represent each partial derivative as a path in a tree
of finite width and infinite heights. A path goes down if the function is
differentiated by the 0th variable. It goes right if there will be no further
differentiation w.r.t. 0th variable, but differentiations w.r.t. remaining
variable can take place. This intuition is depicted by the following
illustration of the ternary case:

$f(\boldsymbol{a})$$f_{x}(\boldsymbol{a})$$f_{x^{2}}(\boldsymbol{a})$$\vdots$$\vdots$$\vdots$$f_{xy}(\boldsymbol{a})$$\vdots$$\vdots$$f_{xz}(\boldsymbol{a})$$\vdots$$f_{y}(\boldsymbol{a})$$f_{y^{2}}(\boldsymbol{a})$$\vdots$$\vdots$$f_{yz}(\boldsymbol{a})$$\vdots$$f_{z}(\boldsymbol{a})$$f_{z^{2}}(ba)$$\vdots$

This can be seen as a special kind of infinite trie (or prefix-tree) of
alphabets $\partial_{x_{i}}$, with available letter eventually decreasing.

This can be implemented by a (co-)inductive type as follows:

data STower n a where ZS :: !a -¿ STower 0 a SS :: !a -¿ STower (n + 1) a -¿
STower n a -¿ STower (n + 1) a

A tree can have an _infinite height_. Since Haskell is a lazy language, this
won’t eat up the memory and only necessary information will be gradually
allocated. Since making everything lazy can introduce a huge space leak, we
force each coefficient a when their corresponding data constructors are
reduced to weak head normal form, as expressed by field strictness annotation
!a.

Then a lifting operation for univariate function is given by:

liftSTower :: forall c n a. (KnownNat n, c a, forall x k. c x =¿ c (STower k
x) ) =¿ (forall x. c x =¿ x -¿ x) -¿ – ^Function (forall x. c x =¿ x -¿ x) -¿
– ^its first-order derivative STower n a -¿ STower n a liftSTower f df (ZS a)
= ZS (f a) liftSTower f df x@(SS a da dus) = SS (f a) (da * df x) (liftSTower
@c f df dus)

Here, we use type-level constraint c to represent to a subclass of smooth
functions, e.g.
$\definecolor{tcbcolback}{rgb}{.75,.75,.75}\definecolor{tcbcolframe}{rgb}{.5,.5,.5}\definecolor{tcbcol@origin}{rgb}{0,0,0}\definecolor{.}{rgb}{0,0,0}\definecolor{.}{rgb}{0,0,0}\leavevmode\hbox
to94.01pt{\vbox to16.4pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower
0.0pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{{}{}{}{}\pgfsys@beginscope\pgfsys@invoke{
}{}{}{}{}{}{}{}{}\definecolor[named]{pgffillcolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@fill{.5}\pgfsys@invoke{
}\pgfsys@fill@opacity{1.0}\pgfsys@invoke{
}{{}{}{{}}}{{}{}{{}}}{}{}{{}{}{{}}}{{}{}{{}}}{}{}{{}{}{{}}}{{}{}{{}}}{}{}{{}{}{{}}}{{}{}{{}}}{}{}\pgfsys@moveto{0.0pt}{3.69884pt}\pgfsys@lineto{0.0pt}{12.69884pt}\pgfsys@curveto{0.0pt}{14.74167pt}{1.656pt}{16.39767pt}{3.69884pt}{16.39767pt}\pgfsys@lineto{90.30617pt}{16.39767pt}\pgfsys@curveto{92.349pt}{16.39767pt}{94.005pt}{14.74167pt}{94.005pt}{12.69884pt}\pgfsys@lineto{94.005pt}{3.69884pt}\pgfsys@curveto{94.005pt}{1.656pt}{92.349pt}{0.0pt}{90.30617pt}{0.0pt}\pgfsys@lineto{3.69884pt}{0.0pt}\pgfsys@curveto{1.656pt}{0.0pt}{0.0pt}{1.656pt}{0.0pt}{3.69884pt}\pgfsys@closepath\pgfsys@fill\pgfsys@invoke{
}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\pgfsys@beginscope\pgfsys@invoke{
}{}{}{}{}{}{}{}{}\definecolor[named]{pgffillcolor}{rgb}{.75,.75,.75}\pgfsys@color@gray@fill{.75}\pgfsys@invoke{
}\pgfsys@fill@opacity{1.0}\pgfsys@invoke{
}{{}{}{{}}}{{}{}{{}}}{}{}{{}{}{{}}}{{}{}{{}}}{}{}{{}{}{{}}}{{}{}{{}}}{}{}{{}{}{{}}}{{}{}{{}}}{}{}\pgfsys@moveto{0.85358pt}{3.69884pt}\pgfsys@lineto{0.85358pt}{12.69884pt}\pgfsys@curveto{0.85358pt}{14.27025pt}{2.12743pt}{15.5441pt}{3.69884pt}{15.5441pt}\pgfsys@lineto{90.30617pt}{15.5441pt}\pgfsys@curveto{91.87758pt}{15.5441pt}{93.15143pt}{14.27025pt}{93.15143pt}{12.69884pt}\pgfsys@lineto{93.15143pt}{3.69884pt}\pgfsys@curveto{93.15143pt}{2.12743pt}{91.87758pt}{0.85358pt}{90.30617pt}{0.85358pt}\pgfsys@lineto{3.69884pt}{0.85358pt}\pgfsys@curveto{2.12743pt}{0.85358pt}{0.85358pt}{2.12743pt}{0.85358pt}{3.69884pt}\pgfsys@closepath\pgfsys@fill\pgfsys@invoke{
}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@fill@opacity{1.0}\pgfsys@invoke{
}{{{}}{{}}{{}}{{}}{{}}{{}}{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{6.5441pt}{5.69884pt}\pgfsys@invoke{
}\hbox{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hbox{\set@color{\ignorespaces\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor\lst@@@set@language{\@listingGroup{ltx_lst_identifier}{c}}}}}}}}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}=\definecolor{tcbcolback}{rgb}{.75,.75,.75}\definecolor{tcbcolframe}{rgb}{.5,.5,.5}\definecolor{tcbcol@origin}{rgb}{0,0,0}\definecolor{.}{rgb}{0,0,0}\definecolor{.}{rgb}{0,0,0}\leavevmode\hbox
to123.31pt{\vbox to16.4pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower
0.0pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{{}{}{}{}\pgfsys@beginscope\pgfsys@invoke{
}{}{}{}{}{}{}{}{}\definecolor[named]{pgffillcolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@fill{.5}\pgfsys@invoke{
}\pgfsys@fill@opacity{1.0}\pgfsys@invoke{
}{{}{}{{}}}{{}{}{{}}}{}{}{{}{}{{}}}{{}{}{{}}}{}{}{{}{}{{}}}{{}{}{{}}}{}{}{{}{}{{}}}{{}{}{{}}}{}{}\pgfsys@moveto{0.0pt}{3.69884pt}\pgfsys@lineto{0.0pt}{12.69884pt}\pgfsys@curveto{0.0pt}{14.74167pt}{1.656pt}{16.39767pt}{3.69884pt}{16.39767pt}\pgfsys@lineto{119.61183pt}{16.39767pt}\pgfsys@curveto{121.65466pt}{16.39767pt}{123.31067pt}{14.74167pt}{123.31067pt}{12.69884pt}\pgfsys@lineto{123.31067pt}{3.69884pt}\pgfsys@curveto{123.31067pt}{1.656pt}{121.65466pt}{0.0pt}{119.61183pt}{0.0pt}\pgfsys@lineto{3.69884pt}{0.0pt}\pgfsys@curveto{1.656pt}{0.0pt}{0.0pt}{1.656pt}{0.0pt}{3.69884pt}\pgfsys@closepath\pgfsys@fill\pgfsys@invoke{
}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\pgfsys@beginscope\pgfsys@invoke{
}{}{}{}{}{}{}{}{}\definecolor[named]{pgffillcolor}{rgb}{.75,.75,.75}\pgfsys@color@gray@fill{.75}\pgfsys@invoke{
}\pgfsys@fill@opacity{1.0}\pgfsys@invoke{
}{{}{}{{}}}{{}{}{{}}}{}{}{{}{}{{}}}{{}{}{{}}}{}{}{{}{}{{}}}{{}{}{{}}}{}{}{{}{}{{}}}{{}{}{{}}}{}{}\pgfsys@moveto{0.85358pt}{3.69884pt}\pgfsys@lineto{0.85358pt}{12.69884pt}\pgfsys@curveto{0.85358pt}{14.27025pt}{2.12743pt}{15.5441pt}{3.69884pt}{15.5441pt}\pgfsys@lineto{119.61183pt}{15.5441pt}\pgfsys@curveto{121.18324pt}{15.5441pt}{122.45709pt}{14.27025pt}{122.45709pt}{12.69884pt}\pgfsys@lineto{122.45709pt}{3.69884pt}\pgfsys@curveto{122.45709pt}{2.12743pt}{121.18324pt}{0.85358pt}{119.61183pt}{0.85358pt}\pgfsys@lineto{3.69884pt}{0.85358pt}\pgfsys@curveto{2.12743pt}{0.85358pt}{0.85358pt}{2.12743pt}{0.85358pt}{3.69884pt}\pgfsys@closepath\pgfsys@fill\pgfsys@invoke{
}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@fill@opacity{1.0}\pgfsys@invoke{
}{{{}}{{}}{{}}{{}}{{}}{{}}{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{6.5441pt}{5.69884pt}\pgfsys@invoke{
}\hbox{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\hbox{\set@color{\ignorespaces\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor\lst@@@set@language{\@listingGroup{ltx_lst_keyword}{Floating}}}}}}}}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}$ for elementary
functions. Constraint of form $\forall xk.\ \texttt{c}\ x=>\texttt{c}\
(\texttt{STower}\ k\ x)$ is an example of so-called _Quantified Constraints_.
This requires c to be implemented for any succinct Tower AD, provided that
their coefficient type, say x, is also an instance of c. This constraint is
used recursively when one implements an actual implementation of instance c
(STower n a). For example, Floating instance (for elementary floating point
operations) can be written as follows:

instance Floating a =¿ Floating (STower n a) where sin = liftSTower @Floating
sin cos cos = liftSTower @Floating cos (negate . sin) exp = liftSTower
@Floating exp exp …

In this way, we can implement Tower AD for a class of smooth function closed
under differentiation, just by specifying an original and their first
derivatives.

More general $n$-ary case of lifting operator is obtained in just the same
way:

liftNAry :: forall c n a m. ( c a, forall x k. (KnownNat k, c x) =¿ c (STower
k x) ) =¿ – — f, an m-ary smooth function (forall x. c x =¿ Vec m x -¿ x) -¿ –
— partial derivatives of f, – wrt. i-th variable in the i-th. Vec m (SmoothFun
c m) -¿ Vec m (STower n a) -¿ STower n a liftNAry f _Nil = constSS
$fNilliftNAryfdfsxss=casesing@lofZero->ZS(f$ constTerm
¡$>xss)Succ(k::SNatk)->SS(f$ constTerm ¡$>xss)(sum$ SV.zipWithSame (gi -¿
topDiffed gi * runSmooth fi xss) dfs xss ) (liftNAry @c f dfs $diffOther<$¿
xss)

diffOther :: STower (n + 1) a -¿ STower n a diffOther (SS __dus) = dus

## References

  * [1] D. Cox, J. Little and D. O’Shea  “Using Algebraic Geometry”  Springer, 2005  URL: http://www.cs.amherst.edu/~dac/uag.html
  * [2] Conal Elliott  “Beautiful differentiation”  In _International Conference on Functional Programming (ICFP)_ , 2009  URL: http://conal.net/papers/beautiful-differentiation
  * [3] Conal Elliott  “The Simple Essence of Automatic Differentiation”  In _Proceedings of the ACM on Programming Languages_ 2.ICFP New York, NY, USA: Association for Computing Machinery, 2018  DOI: 10.1145/3236765
  * [4] Hiromi Ishii  “A Purely Functional Computer Algebra System Embedded in Haskell”  In _Computer Algebra in Scientific Computing_ 11077, LNCS  Springer, Cham, 2018, pp. 288–303  DOI: 10.1007/978-3-319-99639-4˙20
  * [5] Hiromi Ishii  “A Succinct Multivariate Lazy Multivariate Tower AD for Weil Algebra Computation”  In _Computer Algebra_ 2185, RIMS Kôkyûroku  Research Institute for Mathematical Sciences, Kyoto University, Kyoto, Japan, 2021, pp. 104–112  arXiv:2103.11615
  * [6] Hiromi Ishii  “computational-algebra”, 2013  URL: https://konn.github.io/computational-algebra
  * [7] Hiromi Ishii  “smooth: Computational Smooth Infinitesimal Analysis”, 2020  URL: https://github.com/konn/smooth
  * [8] Dominic Joyce  “Algebraic Geometry over $C^{\infty}$-rings”, 2016  arXiv:1001.0023 [math.AG]
  * [9] Edward A. Kmett  “ad: Automatic Differentiation”, 2010  URL: https://hackage.haskell.org/package/ad
  * [10] F William Lawvere  “Categorical dynamics”  In _Topos theoretic methods in geometry_ 30 Mathematical Institute, Aarhus University Various Publications, 1979, pp. 1–28
  * [11] Ieke Moerdijk and Gonzalo E. Reyes  “Models for Smooth Infinitesimal Analysis”  Springer-Verlag New York, 1991  DOI: 10.1007/978-1-4757-4143-8
  * [12] Hirokazu Nishimura and Takeshi Osoekawa  “General Jacobi Identity Revisited Again”  In _International Journal of Theoretical Physics_ 46.11, 2007, pp. 2843–2862  DOI: 10.1007/s10773-007-9397-z
  * [13] Barak Pearlmutter and Jeffrey Siskind  “Lazy multivariate higher-order forward-mode AD”  In _Conference Record of the Annual ACM Symposium on Principles of Programming Languages_ 42, 2007, pp. 155–160  DOI: 10.1145/1190216.1190242
  * [14] Wikipedia  “Automatic Differentiation”, 2021  URL: https://en.wikipedia.org/w/index.php?title=Automatic˙differentiation&oldid=995938170#Automatic˙differentiation˙using˙dual˙numbers


# Inertial Force Transmission in Dense Granular Flows

Matthew Macaulay<EMAIL_ADDRESS>Pierre Rognon
<EMAIL_ADDRESS>School of Civil Engineering, The University of
Sydney, Australia

###### Abstract

Dense granular flows are well described by several continuum models, however,
their internal dynamics remain elusive. This study explores the contact force
distributions in simulated steady and homogenous shear flows. Results
demonstrate the existence of high magnitude contact forces in faster flows
with stiffer grains. A proposed physical mechanism explains this rate-
dependent force transmission. This analysis establishes a relation between
contact forces and grain velocities, providing an entry point to unify a range
of continuum models derived from either contact forces or grain velocity.

Granular flows transmit forces via a network of contacts between pairs of
grains. Knowing the magnitude of these forces is key to predicting the
dynamics of natural and industrial flows. First, contact forces are an
elementary vector of momentum transport [1]. As such, they control the
effective viscosity of the material. Second, large contact forces can lead to
grain comminution, which in turns affects the flow microstructure and its
rheological behaviour [2]. However, while several continuum models can
describe dense granular flows [3, 4], the nature and the origin of their
internal contact force distribution remain poorly understood.

Pioneering works focusing on quasi-static packings revealed that (i) the mean
contact force is controlled by the confining stress $P$ and the grain size $d$
as $Pd^{2}$ and (ii) some contacts are compressed to much higher levels [5,
6]. A piecewise force distribution was manifest, including an exponentially
decaying probability of finding contact forces with a magnitude above the
mean.

Contact forces in collisional flows are radically different as they arise from
binary collisions between grains. In shear flows, Bagnold’s scaling relates
the typical impact velocity to the shear rate $\dot{\gamma}$ and the grain
size $d$ as $\dot{\gamma}d$. Considering a simple elasto-inertial collision
between grains of mass $m$ and contact stiffness $k$ provides an estimate for
the maximum impact contact force of $\dot{\gamma}d\sqrt{km}$ [7]. This
highlights a rate-dependence of the contact forces.

Figure 1: Force and velocity distributions. (a) Plane shear under prescribed
pressure $P$ and shear rate $\dot{\gamma}$, showing a snapshot of the contact
forces (top half, represented by segments of width proportional to the force
magnitude) and of the grain velocities (bottom half) in a flow with $I=0.01$
and $K=10^{-3}$. Movies available online show the evolution of contact forces
and grain velocities within flows at different inertial number and softness.
(b,c) PDF of grain velocity (left) and contact forces (right) for different
inertial numbers (b) and different softness (c).  Figure 2: Contact force
fluctuations $\mathcal{F}$ for different inertial numbers $I$ ($K=10^{-3}$,
main) and different softness $K$ (inset). Symbols: numerical results; lines:
proposed model (4).

In dense granular flows, binary collisions are seldom as grains typically
experience multiple contacts. Their trajectories involve inertial
reorganisation events whose dynamics is independent on, and much faster than
the shear rate [3]. As a result, the grains’ relative velocity deviates from
Bagnold’s scaling [10, 8, 9]. The resulting distribution in grain velocity has
been shown to control the process of shear-induced diffusion [11, 12].
However, how and whether it is related to contact forces is unknown.

In this Letter, we seek to establish a physically-based model that captures
the force distribution in dense granular flows. In this aim, we have measured
contact force distribution in simulated steady and homogenous shear flows and
inferred a physical mechanism involving grain velocity fluctuations that is
able to capture these observations.

## Simulations.

We used a discrete element method to simulate granular flows in a plane shear
configuration (figure 1a). The shear cell is bi-dimensional. Both the normal
stress $P$ and the shear rate $\dot{\gamma}$ are prescribed using periodic
boundary conditions to avoid the flow heterogeneities that solid walls would
induce [13]. The granular material is comprised of $10^{4}$ grains of mean
diameter $d$, mass $m$ and density $\rho$. A uniform polydispersity of $d\pm
20\%$ is introduced on grain’s diameter to prevent shear-induced
crystallisation. Grains interact via pairwise contacts including a normal
elasto-dissipative force and a tangential friction. The supplemental material
further details these granular interactions and the simulation procedure.
Contact parameters include a coefficient of restitution and a coefficient of
friction, which are both set to $0.5$. In contrast, different values of the
contact normal stiffness $k$ will be considered. The effect of these contact
parameters on the granular rheology is discussed in [8].

The following analysis focuses on steady and homogenous flows obtained with
this configuration while prescribing two dimensionless numbers $I$ and $K$:

$I=\dot{\gamma}t_{i};\quad t_{i}=d\sqrt{\frac{\rho}{P}};\quad
K=\frac{t^{2}_{c}}{t^{2}_{i}}=\frac{Pd}{k}.$ (1)

The inertial number $I$ compares the shear time $\dot{\gamma}^{-1}$ to the
inertial time $t_{i}$. $t_{i}$ measures a characteristic time for a grain of
mass $m$ initially at rest to move over a distance $d$ under the action of a
force $Pd^{2}$ [3, 8]. The softness number $K$ measures the elastic
deformation of a grain experiencing a compressive contact force of magnitude
$Pd^{2}$; equivalently, it compares the inertial time $t_{i}$ to an elasto-
inertial collision time $t_{c}=d\sqrt{\frac{\rho d}{k}}$. In the following,
results cover the ranges $10^{-3}\leqslant I\leqslant.9$ and $10^{-5}\leqslant
K\leqslant 10^{-3}$, which correspond to the dense flow regime and the rigid
limit, characterised by $t_{c}<t_{i}<\dot{\gamma}^{-1}$ [3, 8].

Figure 3: Internal kinematic variables:  Velocity fluctuation $\delta v$,
acceleration fluctuation $\delta a$ and velocity memory time $\Psi$ measured
in flows with: (a) different inertial numbers ($K=10^{-3}$); and (b) different
softness $K$ for $I=10^{-2}$ (filled symbols) and $I=10^{-1}$ (open symbols).
Symbols: numerical results; lines: proposed models in (3), (5) and (6).
Numbers summarise the power-law exponents.

## Phenomenological scaling for force fluctuation.

First, we measured the component of the contact forces in the direction
transverse to the shear, $|f_{y}|$, to determine their probability density
function (PDF). Taking advantage of the homogeneity and steadiness of the
flows, these were measured considering all grains in the cell and $300$
sampled times spread evenly over $15$ shear deformations. Figures 1b,c show
that the force distributions are affected by both the inertial number $I$ and
the softness $K$. To characterise the width of these distributions, we
measured their standard deviation $\delta f$ which we refer to as force
fluctuation. Figure 2 shows the combined influence of $I$ and $K$ on the
normalised force fluctuation $\mathcal{F}=\frac{\delta f}{Pd^{2}}$. $\delta f$
converges to $Pd^{2}$ in the quasi-static limit $I\to 0$, which is consistent
with the findings in [5, 6]. In contrast, it becomes significantly larger at
high inertial numbers and for stiffer grains: $\delta f$ reaches values up to
one order of magnitude larger than $Pd^{2}$ in the range of parameters that
have been explored. However, these results do not immediately provide evidence
for a relation between $\delta f$ and the dimensionless numbers $I$ and $K$.

To identify this relation, we use the micro-inertia theory introduced in [1]
that relates contact force fluctuation to grain acceleration. The rationale
for this theory is that multiple contact forces acting on a grain may be
partially balanced. Their balanced part does not induce any acceleration and
scales with $Pd^{2}$. Their unbalanced part drives some grain acceleration
according to Newton’s law of motion. This ultimately yields a linear relation
between force and acceleration fluctuations $\delta f\approx Pd^{2}+m\delta
a$. Here, the acceleration fluctuation is defined as the standard deviation of
the distribution in grain acceleration in the $y$-direction. A dimensionless
formulation of this relation is:

$\mathcal{F}\approx 1+\mathcal{A},$ (2)

where $\mathcal{A}=\frac{\delta a}{d/t_{i}^{2}}$ is the normalised
acceleration fluctuation. The supplemental section shows that this relation is
valid for flows performed at different inertial numbers $I$ and with different
softness $K$. Thanks to this relation, establishing how the force fluctuation
varies with $I$ and $K$ may be achieved by finding how the acceleration
fluctuation varies with $I$ and $K$.

Figure 3 shows that the acceleration fluctuation increases with higher
inertial number and stiffer grains, seemingly with a power law:

$\frac{\delta a}{d\dot{\gamma}^{2}}\propto K^{n_{aK}}I^{n_{aI}}.$ (3)

Results are consistent with exponents $n_{aK}=-0.25$ and $n_{aI}=-1.25$.
Introducing this scaling into the force-acceleration relation (2) leads to the
following phenomenological scaling law for the force fluctuations:

$\mathcal{F}\approx 1+\alpha K^{-0.25}I^{0.75}.$ (4)

Figure 2 shows that this relation reasonably captures the force fluctuation in
the explored range of $I$ and $K$ with $\alpha=2.5$ as a fitting parameter.

Figure 4: Rescaling of force and velocity distribution. Lines are the PDF
presented in figure 1b,c, rescaling grain velocity and contact forces by
$\delta v$ and $\delta f$ defined in the models (5) and (4), respectively.
Filled areas show a normal distribution
$(2\pi)^{-\frac{1}{2}}e^{-\frac{1}{2}\left(|v_{y}|/\delta v\right)^{2}}$
(left) and an exponential distribution $e^{-\beta\frac{|f_{y}|}{\delta f}}$
(right) with $\beta=1.4$ found as best fitting the data.

## Physical process relating velocity and force fluctuations.

To identify a physical process at the origin of this phenomenological law, we
now seek to relate the acceleration fluctuation to other kinematic quantities
including grain velocities and their time persistence.

Figures 1b,c show that the grain vertical velocity distribution is affected by
the inertial number. However, unlike contact forces, the velocity distribution
is not dependent on the softness. Figure 3 reports the standard deviation of
these distributions, which we refer to as velocity fluctuations $\delta v$.
Results suggest the following power law:

$\frac{\delta v}{\dot{\gamma}d}\approx 0.5I^{n_{v}}.$ (5)

As previously proposed [8, 11, 12], they are consistent with a power exponent
$n_{v}=-0.5$. They further point out the negligible influence of the softness
$K$.

Figure 3 also reports the velocity auto-correlation time $\Psi$, defined as
$\Psi=\int_{t=0}^{\infty}\psi(t)dt$, where $\psi(t)=\langle
v_{y}(t_{0})v_{y}(t_{0}+t)\rangle/\langle v_{y}(t_{0})^{2}\rangle$ is the
velocity auto-correlation function and the angle brackets denote the average
operator on all grains and all reference times $t_{0}$. We refer to $\Psi$ as
velocity persistence time, as it measures a characteristic period of time
during which a grain velocity $v_{y}$ is sustained. Results suggest the
following power law:

$\Psi\dot{\gamma}\approx 0.35I^{n_{\Psi}}.$ (6)

As previously reported [9, 11, 12], they are consistent with a power exponent
$n_{\Psi}=0.5$. They further highlight the near-independence of the velocity
memory time on the softness $K$.

From these two observations on the velocity fluctuation and their persistence,
we infer a simplified scenario for a typical grain motion in dense granular
flows. Like a random walk, this scenario involves grains moving over steps of
length $\delta v\Psi$ at a constant velocity $\delta v$ before changing
direction. Accordingly, the acceleration during the step is null. As for
collisional flows, we postulate that the acceleration during the change of
direction is of the order of $\delta v/t_{c}$: this reflects a change in grain
velocity of $\pm\delta v$ in a period of time driven by the elasto-inertial
collision time $t_{c}$.

The characteristic mean square acceleration during such a cycle is
(considering than $t_{c}\ll\Psi$): $\langle
a^{2}\rangle=\frac{\big{(}\frac{\delta
v}{t_{c}}\big{)}^{2}t_{c}+0}{\Psi}=\frac{\delta v^{2}}{t_{c}\Psi}$. This
predicts the following relation between velocity fluctuation and acceleration
fluctuation $\delta a\equiv\langle a^{2}\rangle^{0.5}$:

$\delta a\propto\frac{\delta v}{\sqrt{t_{c}\Psi}}.$ (7)

Introducing this relation into the force-acceleration relation in (2) leads to
the following relation between the force and velocity fluctuations:

$\delta f-Pd^{2}\propto m\frac{\delta v}{\sqrt{t_{c}\Psi}}.$ (8)

Introducing in (8) the scalings for the velocity fluctuation and their
lifetime in (5) and (6) leads to the following scaling law for the force
fluctuation expressed in terms of the inertial number and the softness $K$:
$\mathcal{F}-1\propto K^{-\frac{1}{4}}I^{\frac{5}{4}}$, which is consistent
with the phenomenological model (4). The supplemental section further assesses
the validity of this scenario by comparing the measured $\delta a$ and $\delta
f$ to the predictions of Eqs (7) and (8), respectively.

Finally, figure 4 shows that the distribution of contact forces collapse onto
a single exponential distribution when rescaled by the force fluctuation
defined by (4). This distribution is similar to that found in quasi-static
packings [5, 6]. Here, we find how its width depends on the inertial number
and grain softness in dense flows. Similarly, figure 4 shows that grain
velocity distributions all collapse onto a single normal distribution when
rescaled by the velocity fluctuation modelled by (5). This further confirms
that, unlike contact forces, grain kinematic is independent on their softness.

## Conclusions.

This study highlights a significant rate-dependence of contact forces in dense
granular flows. It reveals the existence of forces much larger than the
characteristic scale $Pd^{2}$ with an exponential distribution, and introduces
a model (4) to capture the width of this distribution as a function of the
inertial number and grain softness. These results can readily be used to
foresee flow conditions likely to lead to grain breakage or abrasive wear
induced by extreme contact forces.

The random-walk mechanism introduced to explain this rate-dependence of the
force distribution highlights a general relationship between contact forces
and grain velocity fluctuation (8). This relationship is an entry point to
help reconcile two distinct existing approaches to understand the constitutive
behaviour of dense granular flows that are based on either contact forces [14,
15, 1] or grain velocity fluctuations [16, 17, 18]. Lastly, we anticipate that
analogous physical processes could govern the force distribution in other soft
materials such as suspensions, foams and emulsions, which would involve each
material’s specific mode of interaction working via an interstitial liquid.

## References

  * Macaulay and Rognon [2020a] M. Macaulay and P. Rognon, Two mechanisms of momentum transfer in granular flows, Physical Review E 101, 050901 (2020a).
  * Marks and Einav [2015] B. Marks and I. Einav, A mixture of crushing and segregation: the complexity of grainsize in natural granular flows, Geophysical Research Letters 42, 274 (2015).
  * MiDi [2004] G. MiDi, On dense granular flows, The European Physical Journal E 14, 341 (2004).
  * Kamrin and Koval [2012] K. Kamrin and G. Koval, Nonlocal constitutive relation for steady granular flow, Physical Review Letters 108, 178301 (2012).
  * Radjai  _et al._ [1998] F. Radjai, D. E. Wolf, M. Jean, and J.-J. Moreau, Bimodal character of stress transmission in granular packings, Physical review letters 80, 61 (1998).
  * Majmudar and Behringer [2005] T. S. Majmudar and R. P. Behringer, Contact force measurements and stress-induced anisotropy in granular materials, Nature 435, 1079 (2005).
  * Macaulay and Rognon [2020b] M. Macaulay and P. Rognon, Viscosity of cohesive granular flows, Soft Matter (2020b).
  * Da Cruz  _et al._ [2005] F. Da Cruz, S. Emam, M. Prochnow, J.-N. Roux, and F. Chevoir, Rheophysics of dense granular materials: Discrete simulation of plane shear flows, Physical Review E 72, 021309 (2005).
  * DeGiuli  _et al._ [2016] E. DeGiuli, J. McElwaine, and M. Wyart, A phase diagram unifies energy dissipation, kinetics, and rheology in inertial granular flows, Bulletin of the American Physical Society 61 (2016).
  * Radjai and Roux [2002] F. Radjai and S. Roux, Turbulentlike fluctuations in quasistatic flow of granular media, Physical review letters 89, 064302 (2002).
  * Kharel and Rognon [2017] P. Kharel and P. Rognon, Vortices enhance diffusion in dense granular flows, Physical review letters 119, 178001 (2017).
  * Macaulay and Rognon [2019] M. Macaulay and P. Rognon, Shear-induced diffusion in cohesive granular flows: effect of enduring clusters, Journal of Fluid Mechanics 858 (2019).
  * Rognon  _et al._ [2015] P. G. Rognon, T. Miller, B. Metzger, and I. Einav, Long-range wall perturbations in dense granular flows, Journal of Fluid Mechanics 764, 171 (2015).
  * Pouliquen and Forterre [2009] O. Pouliquen and Y. Forterre, A non-local rheology for dense granular flows, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 367, 5091 (2009).
  * Azéma and Radjaï [2014] E. Azéma and F. Radjaï, Internal structure of inertial granular flows, Physical Review Letters 112, 1 (2014).
  * Zhang and Kamrin [2017] Q. Zhang and K. Kamrin, Microscopic description of the granular fluidity field in nonlocal flow modeling, Physical Review Letters 118, 1 (2017).
  * Kim and Kamrin [2020] S. Kim and K. Kamrin, Power-law scaling in granular rheology across flow geometries, Phys. Rev. Lett. 125, 088002 (2020).
  * Gaume  _et al._ [2020] J. Gaume, G. Chambon, and M. Naaim, Microscopic origin of nonlocal rheology in dense granular materials, Phys. Rev. Lett. 125, 188001 (2020).


is symmetric with respect to the two indices, hence we can omit the directions
and consider non-directed trees. This gives rise to the following diagrammatic
expression of the partition function

$\displaystyle Z$
$\displaystyle=1+\sum_{(\theta_{1},x_{1})}\quad\raisebox{-12.0pt}{\includegraphics[width=6.50403pt]{1particle.pdf}}$
$\displaystyle\quad+\frac{1}{2!}\sum_{(\theta_{1},x_{1})}\sum_{(\theta_{2},x_{2})}\left(\raisebox{-15.0pt}{\includegraphics[width=108.405pt]{2particle.pdf}}\right)$
$\displaystyle\quad+\frac{1}{3!}\sum_{(\theta_{1},x_{1})}\sum_{(\theta_{2},x_{2})}\sum_{(\theta_{3},x_{3})}\left(\raisebox{-90.0pt}{\includegraphics[width=260.17464pt]{3particle.pdf}}\right)$
$\displaystyle\quad+\cdots,$ (158)

where $\sum_{(\theta,x)}=\int_{[-rL,rL]\times\mathbb{R}}{\rm d}x{\rm
d}\theta/(2\pi)$ and the ellipsis contains terms with higher numbers of
integrals. The integration on the phase space then yields

$Z=\,\raisebox{-24.2pt}{\includegraphics[width=303.53267pt]{partition_integrated}},$
(159)

where a blue box refers to a vertex integrated over the phase space. We have
now lost both directions and vertex labels, and so we consider “types” of
trees (unoriented, unlabelled trees with $N$ vertices). In general, a forest
made of $M$ different types of trees, each type $m$ appearing with
multiplicity $n_{m}$, with $\sum_{m}n_{m}=N$ vertices, has value

$\frac{1}{n_{1}!}V_{1}^{n_{1}}\frac{1}{n_{2}!}V_{2}^{n_{2}}\cdots\frac{1}{n_{M}!}V_{M}^{n_{M}},$
(160)

where $V_{m}$ is the value of the tree of type $m$. Here $1/n_{m}!$ is simply
a symmetry factor for the copies of the tree of type $m$. A tree of type $m$
may have further internal symmetry factors. A tree with $N$ vertices has $N-1$
edges, and the exact expression for a tree of type $m$ with edges $e_{j}$
(sets of two vertices) has expression

$V_{m}=\frac{1}{S_{m}}\frac{1}{(2\pi)^{N}}\int_{[-rL,rL]^{N}\times\mathbb{R}^{N}}\prod_{j=1}^{N}{\rm
d}x_{j}{\rm
d}\theta_{j}\,e^{-\sum_{i=1}^{N}\beta(x_{i}/L,\theta_{i})}\prod_{j=1}^{N-1}\psi_{x\theta}(x_{{\rm
e}_{j}},\theta_{{\rm e}_{j}}).$ (161)

such as in the last term of (159), which is explicitly

$V_{m}=\frac{1}{2!}\frac{1}{(2\pi)^{3}}\int_{[-rL,rL]^{3}\times\mathbb{R}^{3}}\prod_{j=1}^{3}{\rm
d}x_{j}{\rm
d}\theta_{j}\,e^{-\sum_{i=1}^{3}\beta(x_{i}/L,\theta_{i})}\psi_{x\theta}(x_{12},\theta_{12})\psi_{x\theta}(x_{13},\theta_{13}).$
(162)

Taking a further summation over $M$, the partition function turns into an
exponential of the sum over all the connected trees with an internal symmetry
factor:

$\displaystyle Z$
$\displaystyle=\sum_{M=0}^{\infty}\sum_{\begin{subarray}{c}n_{1},n_{2},\cdots\\\
\sum_{m}n_{m}=M\end{subarray}}\frac{1}{n_{1}!}V_{1}^{n_{1}}\frac{1}{n_{2}!}V_{2}^{n_{2}}\cdots\frac{1}{n_{M}!}V_{M}^{n_{M}}=\exp\left(\sum_{j=1}^{\infty}V_{j}\right)$
$\displaystyle=\exp\left[\int_{[-rL,rL]\times\mathbb{R}}\frac{{\rm d}x{\rm
d}\theta}{2\pi}Y(x,\theta)\right],$ (163)

where $Y(x,\theta)$ is the generating function of trees (i.e. sum of all the
values of the trees) that are rooted at $(x,\theta)$ [55].

##  Appendix B Thermodynamic Bethe ansatz from a saddle-point analysis

In this section we go beyond the GGE states discussed in Section 3.4 and study
local equilibrium states using large deviation theory. Specifically we study
measures of the form:

$\displaystyle{\rm d}{\varrho}$
$\displaystyle=\frac{1}{Z}\bigoplus_{N=0}^{\infty}{\rm d}^{N}{\bm{x}}{\rm
d}^{N}{\bm{p}}e^{-\int{\rm d}{x}{\rm
d}{\theta}\beta(x/L,\theta)q_{\theta}(\bm{x},\bm{p};x)}$ (164)

in the large-scale limit $L\to\infty$. Note that ${\rm d}^{N}{\bm{x}}{\rm
d}^{N}{\bm{p}}={\rm d}^{N}{\bm{y}}{\rm d}^{N}{\bm{\theta}}$ since they are
related by a canonical transformation. Using the expressions for the charge
densities (49) we can write this as:

$\displaystyle{\rm d}{\varrho}$ $\displaystyle=\frac{1}{Z}\bigoplus_{N}{\rm
d}^{N}{\bm{y}}{\rm
d}^{N}{\bm{\theta}}e^{-\sum_{i=1}^{N}\beta(x_{i}(\bm{y},\bm{\theta})/L,\theta_{i})}$
(165)

This is almost the distribution of non-interacting particles, however the
positions of particles are shifted due to the presence of other particles. Let
rescale $x=L\bar{x}$ and $y=L\bar{y}$ coordinates by $L$ (i.e. going to
macroscopic coordinates) and write (3) in the following way:

$\displaystyle\bar{y}_{i}(\bm{\bar{}}{x},\bm{\theta})$
$\displaystyle=\bar{x}_{i}+\frac{1}{L}\sum_{j\neq
i}\psi_{\theta}(\bar{x}_{i}-\bar{x}_{j},\theta_{i}-\theta_{j})=\bar{x}_{i}+\int{\rm
d}{\bar{x}^{\prime}}{\rm d}{\theta^{\prime}}\rho_{\rm
e}(\bar{x}^{\prime},\theta^{\prime})\psi_{\theta}(L(\bar{x}_{i}-\bar{x}^{\prime}),\theta_{i}-\theta^{\prime})$
(166) $\displaystyle=Y[\rho_{\rm e}](\bar{x}_{i},\theta_{i})$ (167)

By this we see that $\bar{y}_{i}$ is determined only by $\bar{x}_{i}$ and a
mean field coupling to the empirical quasi-particle density:

$\displaystyle\rho_{\rm e}(\bar{x},\theta)$
$\displaystyle=\frac{1}{L}\sum_{i}\delta(\bar{x}-\bar{x}_{i})\delta(\theta-\theta_{i})$
(168)

Note that $\bar{x}\to Y[\rho_{\rm e}](\bar{x},\theta)$ is invertible for fixed
$\rho_{\rm e}$ and $\theta$ since:

$\displaystyle\frac{{\rm d}}{{\rm d}\bar{x}}Y[\rho_{\rm
e}](\bar{x},\theta)=1+\int{\rm d}{\bar{x}^{\prime}}{\rm
d}{\theta^{\prime}}\rho_{\rm
e}(\bar{x}^{\prime},\theta^{\prime})\psi_{x\theta}(L(\bar{x}-\bar{x}^{\prime}),\theta-\theta^{\prime})>0.$
(169)

Denote its inverse for fixed $\rho_{\rm e}$ and $\theta$ by $X[\rho_{\rm
e}](\bar{y},\theta)$ and write:

$\displaystyle{\rm d}{\varrho}$ $\displaystyle=\frac{1}{Z}\bigoplus_{N}{\rm
d}^{N}{\bar{\bm{y}}}{\rm
d}^{N}{\bm{\theta}}\prod_{i=1}^{N}e^{-\beta(X[\rho_{\rm
e}](\bar{y}_{i},\theta_{i}),\theta_{i})}$ (170)
$\displaystyle=\frac{1}{Z}\bigoplus_{N}{\rm d}^{N}{\bar{\bm{y}}}{\rm
d}^{N}{\bm{\theta}}\exp\left(-L\int{\rm d}{\bar{y}}{\rm d}{\theta}\rho_{\rm
e,y}(\bar{y},\theta)\beta(X[\rho_{\rm e}](\bar{y},\theta),\theta)\right)$
(171) $\displaystyle=\frac{1}{Z}\mathcal{D}[\rho_{\rm
e,y}]\exp\left(LS[\rho_{\rm e,y}]-L\int{\rm d}{\bar{y}}{\rm
d}{\theta}\rho_{\rm e,y}(\bar{y},\theta)\beta(X[\rho_{\rm
e}](\bar{y},\theta),\theta)\right)$ (172)

where we defined $\rho_{\rm
e,y}(\bar{y},\theta)=\frac{1}{L}\sum_{i}\delta(\bar{y}-\bar{y}_{i})\delta(\theta-\theta_{i})$
and in the last step changed the measure from ${\rm d}^{N}{\bar{\bm{y}}}{\rm
d}^{N}{\bm{\theta}}$ to $\mathcal{D}[\rho_{\rm e,y}]e^{LS[\rho_{\rm e,y}]}$,
which introduces the entropy $S[\rho_{\rm e,y}]=-\int{\rm d}{\bar{y}}{\rm
d}{\theta}\rho_{\rm e,y}(\bar{y},\theta)\log\rho_{\rm e,y}(\bar{y},\theta)$
(i.e. the number of microstates per macrostate) in the change of measure. Note
that $\rho_{\rm e}(\bar{x},\theta)$ and $\rho_{\rm e,y}(\bar{y},\theta)$ are
related via a push-forward $\rho_{\rm e,y}=Y[\rho_{\rm e}]_{*}\rho_{\rm e}$,
i.e.:

$\displaystyle\rho_{\rm e,y}(\bar{y},\theta)$ $\displaystyle=\rho_{\rm
e}(X[\rho_{\rm e}](\bar{y},\theta),\theta)\frac{{\rm d}{X}[\rho_{\rm
e}](\bar{y},\theta)}{{\rm d}{\bar{y}}}$ (173)

By the large deviation principle we expect that as $L\to\infty$ the measure
${\rm d}{\varrho}$ will become concentrated on the minimizer of the rate
function:

$\displaystyle I[\rho_{\rm e,y}]$ $\displaystyle=\int{\rm d}{\bar{y}}{\rm
d}{\theta}\rho_{\rm e,y}(\bar{y},\theta)\log\rho_{\rm
e,y}(\bar{y},\theta)+\int{\rm d}{\bar{y}}{\rm d}{\theta}\rho_{\rm
e,y}(\bar{y},\theta)\beta(X[\rho_{\rm e}](\bar{y},\theta),\theta)$ (174)
$\displaystyle=\int{\rm d}{\bar{x}}{\rm d}{\theta}\rho_{\rm
e}(\bar{x},\theta)\log\left(\frac{\rho_{\rm e}(\bar{x},\theta)}{\frac{{\rm
d}{Y}[\rho_{\rm e}](\bar{x},\theta)}{{\rm d}{\bar{x}}}}\right)+\int{\rm
d}{\bar{x}}{\rm d}{\theta}\rho_{\rm e}(\bar{x},\theta)\beta(\bar{x},\theta)$
(175) $\displaystyle=\int{\rm d}{\bar{x}}{\rm d}{\theta}\rho_{\rm
e}(\bar{x},\theta)\log\left(\frac{\rho_{\rm
e}(\bar{x},\theta)e^{\beta(\bar{x},\theta)}}{1+\int{\rm
d}{\bar{x}^{\prime}}{\rm d}{\theta^{\prime}}\rho_{\rm
e}(\bar{x}^{\prime},\theta^{\prime})\psi_{x\theta}(L(\bar{x}-\bar{x}^{\prime}),\theta-\theta^{\prime})}\right)$
(176) $\displaystyle\to\int{\rm d}{\bar{x}}{\rm d}{\theta}\rho_{\rm
e}(\bar{x},\theta)\log\left(\frac{\rho_{\rm
e}(\bar{x},\theta)e^{\beta(\bar{x},\theta)}}{1+\int{\rm
d}{\theta^{\prime}}\rho_{\rm
e}(\bar{x},\theta^{\prime})\varphi(\theta-\theta^{\prime})}\right)$ (177)

Here we used that $\psi_{x\theta}(Lx,\theta)$ approaches
$\delta(x)\varphi(\theta)$ as $L\to\infty$. The minimizer of $I[\rho_{\rm
e,y}]$ satisfies:

$\displaystyle 0$ $\displaystyle=\log\left(\frac{\rho_{\rm
e}(\bar{x},\theta)e^{\beta(\bar{x},\theta)}}{1+\int{\rm
d}{\theta^{\prime}}\rho_{\rm
e}(\bar{x},\theta^{\prime})\varphi(\theta-\theta^{\prime})}\right)-\int{\rm
d}{\theta^{\prime}}\varphi(\theta-\theta^{\prime})\frac{\rho_{\rm
e}(\bar{x},\theta^{\prime})}{1+\int{\rm d}{\theta^{\prime\prime}}\rho_{\rm
e}(\bar{x},\theta^{\prime\prime})\varphi(\theta^{\prime}-\theta^{\prime\prime})}$
(178)

In analogy with the terminology of GHD let us define the ‘occupation function’
$\rho_{\rm n}(\bar{x},\theta)=\frac{2\pi\rho_{\rm
e}(\bar{x},\theta)}{1+\int{\rm d}{\theta^{\prime}}\rho_{\rm
e}(\bar{x},\theta^{\prime})\varphi(\theta-\theta^{\prime})}$, which satisfies:

$\displaystyle\rho_{\rm n}(\bar{x},\theta)e^{-\int\frac{{\rm
d}{\theta^{\prime}}}{2\pi}\varphi(\theta-\theta^{\prime})\rho_{\rm
n}(\bar{x},\theta^{\prime})}$ $\displaystyle=2\pi e^{-\beta(\bar{x},\theta)}$
(179)

This set of equations decouple for different $x$, which shows the local
density approximation.

## References

  * [1] Andrea Cavaglià, Stefano Negro, István M. Szécsényi and Roberto Tateo  “$T\bar{T}$-deformed 2D quantum field theories”  In _JHEP_ 2016.10 Springer ScienceBusiness Media LLC, 2016, pp. 112  DOI: 10.1007/jhep10(2016)112
  * [2] F.. Smirnov and A.. Zamolodchikov  “On space of integrable quantum field theories”  In _Nucl. Phys. B_ 915, 2017, pp. 363–383  DOI: 10.1016/j.nuclphysb.2016.12.014
  * [3] Alexander B. Zamolodchikov  “Expectation value of composite field T anti-T in two-dimensional quantum field theory”  In _preprint arXiv:hep-th/0401146_ , 2004  arXiv: https://arxiv.org/abs/hep-th/0401146
  * [4] Balázs Pozsgay, Yunfeng Jiang and Gábor Takács  “$T\bar{T}$-deformation and long range spin chains”  In _JHEP_ 2020.3, 2020, pp. 92  DOI: 10.1007/JHEP03(2020)092
  * [5] Balázs Pozsgay et al.  “Integrable spin chain with Hilbert space fragmentation and solvable real-time dynamics”  In _Phys. Rev. E_ 104 American Physical Society, 2021, pp. 044106  DOI: 10.1103/PhysRevE.104.044106
  * [6] Balázs Pozsgay, Tamás Gombor and Arthur Hutsalyuk  “Integrable hard-rod deformation of the Heisenberg spin chains”  In _Phys. Rev. E_ 104 American Physical Society, 2021, pp. 064124  DOI: 10.1103/PhysRevE.104.064124
  * [7] Márton Borsi, Levente Pristyák and Balázs Pozsgay  “Matrix Product Symmetries and Breakdown of Thermalization from Hard Rod Deformations”  In _Phys. Rev. Lett._ 131 American Physical Society, 2023, pp. 037101  DOI: 10.1103/PhysRevLett.131.037101
  * [8] J Cardy  “$T\bar{T}$ deformations of non-Lorentz invariant field theories”  In _preprint arXiv:1809.07849_ , 2018  URL: https://arxiv.org/abs/1809.07849
  * [9] Chantelle Esper and Sergey Frolov  “$T\bar{T}$ deformations of non-relativistic models”  In _JHEP_ 2021.6 Springer, 2021, pp. 1–34  URL: https://link.springer.com/article/10.1007/JHEP06(2021)101
  * [10] John Cardy and Benjamin Doyon  “$T\bar{T}$ deformations and the width of fundamental particles”  In _JHEP_ 2022.4 Springer, 2022, pp. 1–27  URL: https://link.springer.com/article/10.1007/JHEP04(2022)136
  * [11] Yunfeng Jiang  “$T\bar{T}$-deformed 1d Bose gas”  In _SciPost Physics_ 12.6, 2022, pp. 191  DOI: 10.21468/SciPostPhys.12.6.191
  * [12] Guzmán Hernández-Chifflet, Stefano Negro and Alessandro Sfondrini  “Flow equations for generalised $T\bar{T}$ deformations” arXiv: 1911.12233  In _Phys. Rev. Lett._ 124.20, 2020, pp. 200601  DOI: 10.1103/PhysRevLett.124.200601
  * [13] Benjamin Doyon, Joseph Durnin and Takato Yoshimura  “The Space of Integrable Systems from Generalised $T\bar{T}$-Deformations”  In _SciPost Phys._ 13 SciPost, 2022, pp. 072  DOI: 10.21468/SciPostPhys.13.3.072
  * [14] Olalla A. Castro-Alvaredo, Stefano Negro and Fabio Sailis  “Form factors and correlation functions of $\textrm{T}\overline{\textrm{T}}$-deformed integrable quantum field theories”  In _JHEP_ 09, 2023, pp. 048  DOI: 10.1007/JHEP09(2023)048
  * [15] Olalla A. Castro-Alvaredo, Stefano Negro and Fabio Sailis  “Completing the Bootstrap Program for $\mathrm{T}\overline{\mathrm{T}}$-Deformed Massive Integrable Quantum Field Theories”  In _preprint arXiv:2305.17068_ , 2023  arXiv: https://arxiv.org/abs/2305.17068
  * [16] O.. Castro-Alvaredo, S. Negro and I.. Szécsényi  “On the Representation of Minimal Form Factors in Integrable Quantum Field Theory”  In _preprint arXiv:2311.16955_ , 2023  arXiv: https://arxiv.org/abs/2311.16955
  * [17] Yunfeng Jiang  “A pedagogical review on solvable irrelevant deformations of 2D quantum field theory”  In _Commun. Theor. Phys._ 73.5 IOP Publishing, 2021, pp. 057201  URL: https://iopscience.iop.org/article/10.1088/1572-9494/abe4c9
  * [18] Andrea Hubacher  “Classical Scattering Theory in One Dimension”  In _Commun. Math. Phys._ 123, 1989, pp. 353–375  URL: https://link.springer.com/article/10.1007/BF01238807
  * [19] V.. Arnold  “Mathematical Methods of Classical Mechanics” 60, Graduate Texts in Mathematics  Springer, 1978  URL: https://link.springer.com/book/10.1007/978-1-4757-1693-1
  * [20] Stephen Parke  “Absence of particle production and factorization of the S-matrix in 1 + 1 dimensional models”  In _Nucl. Phys. B_ 174.1, 1980, pp. 166–182  DOI: 10.1016/0550-3213(80)90196-0
  * [21] Benjamin Doyon  “Generalised hydrodynamics of the classical Toda system”  In _J. Math. Phys._ 60.7, 2019, pp. 073302  DOI: 10.1063/1.5096892
  * [22] Benjamin Doyon  “Lecture Notes On Generalised Hydrodynamics”  In _SciPost Phys. Lect. Notes_ SciPost, 2020, pp. 18  DOI: 10.21468/SciPostPhysLectNotes.18
  * [23] C. Boldrighini, R.. Dobrushin and Yu.. Sukhov  “One-dimensional hard rod caricature of hydrodynamics”  In _J. Stat. Phys._ 31.3 Springer ScienceBusiness Media LLC, 1983, pp. 577–616  DOI: 10.1007/bf01019499
  * [24] Berislav Bu$$a, Katja Klobas and Toma$$ Prosen  “Rule 54: exactly solvable model of nonequilibrium statistical mechanics”  In _J. Stat. Mech. Theory Exp._ 2021.7 IOP PublishingSISSA, 2021, pp. 074001  DOI: 10.1088/1742-5468/ac096b
  * [25] David A. Croydon and Makiko Sasada  “Generalized Hydrodynamic Limit for the Box–Ball System”  In _Commun. Math. Phys._ 383.1 Springer ScienceBusiness Media LLC, 2020, pp. 427–463  DOI: 10.1007/s00220-020-03914-x
  * [26] C.. Yang and C.. Yang  “Thermodynamics of a One Dimensional System of Bosons with Repulsive Delta Function Interaction”  In _J. Stat. Phys._ 10.7, 1969, pp. 1115–1122  DOI: 10.1063/1.1664947
  * [27] Al.B. Zamolodchikov  “Thermodynamic Bethe ansatz in relativistic models: Scaling 3-state Potts and Lee-Yang models”  In _Nucl. Phys. B_ 342.3 Elsevier BV, 1990, pp. 695–720  DOI: 10.1016/0550-3213(90)90333-9
  * [28] Minoru Takahashi  “Thermodynamics of one-dimensional solvable models”  Cambridge University Press, 1999
  * [29] Jorn Mossel and Jean-Sébastien Caux  “Generalized TBA and generalized Gibbs”  In _J. Phys. A: Math. Theor._ 45.25 IOP Publishing, 2012, pp. 255001  URL: https://iopscience.iop.org/article/10.1088/1751-8113/45/25/255001
  * [30] Olalla A. Castro-Alvaredo, Benjamin Doyon and Takato Yoshimura  “Emergent Hydrodynamics in Integrable Quantum Systems Out of Equilibrium”  In _Phys. Rev. X_ 6 American Physical Society, 2016, pp. 041065  DOI: 10.1103/PhysRevX.6.041065
  * [31] Bruno Bertini, Mario Collura, Jacopo De Nardis and Maurizio Fagotti  “Transport in Out-of-Equilibrium $XXZ$ Chains: Exact Profiles of Charges and Currents”  In _Phys. Rev. Lett._ 117 American Physical Society, 2016, pp. 207201  DOI: 10.1103/PhysRevLett.117.207201
  * [32] Herbert Spohn  “Hydrodynamic scales of integrable many-particle systems”  In _preprint arXiv:2301.08504_ , 2023  URL: https://arxiv.org/abs/2301.08504
  * [33] Pablo A Ferrari, Chiara Franceschini, Dante GE Grevino and Herbert Spohn  “Hard rod hydrodynamics and the Levy Chentsov field”  In _preprint arXiv:2211.11117_ , 2022  URL: https://arxiv.org/abs/2211.11117
  * [34] Jürgen Moser  “Finitely many mass points on the line under the influence of an exponential potential – an integrable system”  In _Dynamical Systems, Theory and Applications: Battelle Seattle 1974 Rencontres_ Berlin, Heidelberg: Springer Berlin Heidelberg, 1975, pp. 467–497  DOI: 10.1007/3-540-07171-7˙12
  * [35] Enej Ilievski, Marko Medenjak, Tomaž Prosen and Lenart Zadnik  “Quasilocal charges in integrable lattice systems”  In _J. Stat. Mech. Theory Exp._ 2016.6 IOP Publishing, 2016, pp. 064008  DOI: 10.1088/1742-5468/2016/06/064008
  * [36] Benjamin Doyon  “Thermalization and Pseudolocality in Extended Quantum Systems”  In _Commun. Math. Phys._ 351.1 Springer ScienceBusiness Media LLC, 2017, pp. 155–200  DOI: 10.1007/s00220-017-2836-7
  * [37] Jacopo De Nardis, Benjamin Doyon, Marko Medenjak and Miłosz Panfil  “Correlation functions and transport coefficients in generalised hydrodynamics”  In _J. Stat. Mech. Theory Exp._ 2022.1 IOP Publishing, 2022, pp. 014002  URL: https://iopscience.iop.org/article/10.1088/1742-5468/ac3658
  * [38] Enej Ilievski, Eoin Quinn and Jean-Sébastien Caux  “From interacting particles to equilibrium statistical ensembles”  In _Phys. Rev. B_ 95 American Physical Society, 2017, pp. 115128  DOI: 10.1103/PhysRevB.95.115128
  * [39] N.. Khimchenko and Ya. Sinai  “On the description of classical reflectionless potentials”  In _Rep. Math. Phys._ 20, 1984, pp. 53–63  URL: https://www.sciencedirect.com/science/article/abs/pii/0034487784900715
  * [40] Neel Malvania et al.  “Generalized hydrodynamics in strongly interacting 1D Bose gases”  In _Science_ 373.6559, 2021, pp. 1129–1133  DOI: 10.1126/science.abf0147
  * [41] Benjamin Doyon and Friedrich Hübner  “Ab initio derivation of generalised hydrodynamics from a gas of interacting wave packets”  In _preprint arXiv:2307.09307_ , 2023  arXiv: https://arxiv.org/abs/2307.09307
  * [42] Herbert Spohn  “Large Scale Dynamics of Interacting Particles”  Springer Berlin Heidelberg, 1991  DOI: 10.1007/978-3-642-84371-6
  * [43] Benjamin Doyon, Herbert Spohn and Takato Yoshimura  “A geometric viewpoint on generalized hydrodynamics”  In _Nucl. Phys. B_ 926 Elsevier BV, 2018, pp. 570–583  DOI: 10.1016/j.nuclphysb.2017.12.002
  * [44] Riccardo Conti, Stefano Negro and Roberto Tateo  “The $T\bar{T}$ perturbation and its geometric interpretation”  In _JHEP_ 2019.2 Springer ScienceBusiness Media LLC, 2019  DOI: 10.1007/jhep02(2019)085
  * [45] Benjamin Doyon, Friedrich Hübner and Takato Yoshimura  “New classical integrable systems from generalized $T\bar{T}$-deformations”  In _preprint arXiv:2311.06369_ , 2023  arXiv: https://arxiv.org/abs/2311.06369
  * [46] Olivier Babelon and Denis Bernard  “The sine-Gordon solitons as an N-body problem”  In _Phys. Lett. B_ 317.3 Elsevier, 1993, pp. 363–368  URL: https://www.sciencedirect.com/science/article/abs/pii/037026939391009C
  * [47] O Babelon, Denis Bernard and FA Smirnov  “Quantization of solitons and the restricted sine-Gordon model”  In _Commun. Math. Phys._ 182 Springer, 1996, pp. 319–354  URL: https://link.springer.com/article/10.1007/BF02517893
  * [48] O Babelon, Denis Bernard and FA Smirnov  “Form factors, KdV and deformed hyperelliptic curves”  In _Nucl. Phys. B-Proceedings Supplements_ 58 Elsevier, 1997, pp. 21–33  URL: https://www.sciencedirect.com/science/article/abs/pii/S0920563297004106
  * [49] J.. Percus  “Equilibrium state of a classical fluid of hard rods in an external field”  In _J. Stat. Phys._ 15.6, 1976, pp. 505–511  DOI: 10.1007/BF01020803
  * [50] Vir B Bulchandani  “Modified Enskog equation for hard rods”  In _preprint arXiv:2309.15846_ , 2023  URL: https://arxiv.org/abs/2309.15846
  * [51] Lev Vidmar and Marcos Rigol  “Generalized Gibbs ensemble in integrable lattice models”  In _J. Stat. Mech. Theory Exp._ 2016.6 IOP Publishing, 2016, pp. 064007  URL: https://iopscience.iop.org/article/10.1088/1742-5468/2016/06/064007
  * [52] Fabian HL Essler and Maurizio Fagotti  “Quench dynamics and relaxation in isolated integrable quantum spin chains”  In _J. Stat. Mech. Theory Exp._ 2016.6 IOP Publishing, 2016, pp. 064002  URL: https://iopscience.iop.org/article/10.1088/1742-5468/2016/06/064002
  * [53] Michel Gaudin  “The Bethe wavefunction”  Cambridge University Press, 2014
  * [54] Jorrit Kruthoff and Onkar Parrikar  “On the flow of states under $T\overline{T}$”  In _SciPost Phys._ 9 SciPost, 2020, pp. 78  DOI: 10.21468/SciPostPhys.9.5.078
  * [55] Ivan Kostov, Didina Serban and Dinh-Long Vu  “TBA and Tree Expansion”  In _Quantum Theory and Symmetries with Lie Theory and Its Applications in Physics Volume 2_ Singapore: Springer Singapore, 2018, pp. 77–98
  * [56] Arthur Cayley  “A theorem on trees”  In _The Collected Mathematical Papers_ 13, Cambridge Library Collection - Mathematics  Cambridge University Press, 2009, pp. 26–28  DOI: 10.1017/CBO9780511703799.010
  * [57] Pierre Suret et al.  “Soliton Gas: Theory, Numerics and Experiments”  In _preprint arXiv:2304.06541_ , 2023  URL: https://arxiv.org/abs/2304.06541
  * [58] Seth Chaiken  “A Combinatorial Proof of the All Minors Matrix Tree Theorem”  In _SIAM Journal on Algebraic Discrete Methods_ 3.3, 1982, pp. 319–329  DOI: 10.1137/0603033