# Formal first integrals and higher variational equations

Sergi Simon

► To cite this version:

Sergi Simon. Formal first integrals and higher variational equations. 2023. [https://hal.science/hal-03929379v2](https://hal.science/hal-03929379v2)

HAL Id: hal-03929379

[https://hal.science/hal-03929379v2](https://hal.science/hal-03929379v2)

Preprint submitted on 6 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Formal first integrals and higher variational equations

Sergi Simon

► To cite this version:

Sergi Simon. Formal first integrals and higher variational equations. 2024. hal-03929379v3

HAL Id: hal-03929379

https://hal.science/hal-03929379v3

Preprint submitted on 3 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Formal first integrals and higher variational equations

Sergi Simon

School of Mathematics and Physics, University of Portsmouth
Lion Gate Bldg, Lion Terrace Portsmouth PO1 3HF, UK
sergi.simon@port.ac.uk

August 10, 2023

## Abstract

The question of how Algebra can be used to solve dynamical systems and characterize chaos was first posed in a fertile mathematical context by Ziglin, Morales, Ramis and Simó using differential Galois theory. Their study was aimed at first-order, later higher-order, variational equations of Hamiltonian systems. Recent work by this author formalized a compact yet comprehensive expression of higher-order variationals as one infinite linear system, thereby simplifying the approach. More importantly, the dual of this linear system contains all information relevant to first integrals, regardless of whether the original system is Hamiltonian. This applicability to formal calculation of conserved quantities is the centerpiece of this paper, following an introduction to the requisite context. Three important examples, namely particular cases of Dixon's system, the SIR epidemiological model with vital dynamics and the Van der Pol oscillator, are tackled, and explicit convergent first integrals are provided for the first two.

**Keywords.** Integrability, Ziglin-Morales-Ramis-Simó theory, formal calculus, chaos, Dixon's system, SIR epidemiological model.

**2000 Mathematics Subject Classification:** 34A05, 37C79, 37J99, 37J30, 34M15, 34C28, 37C10, 15A69, 16W60, 13F25, 37N25.

# 1 Introduction

Given an arbitrary dynamical system, the formulation of its higher variational equations as a linear (infinite) system (LVE*) has shown potential to make strong inroads in the study of integrability [34, 35]. The study, part of which we will explain more in detail in §1.1, is twofold:

(i) on one hand, the original set of equations (LVE) is amenable to the Ziglin-Morales-Ramis-Simó non-integrability framework whenever the system is *Hamiltonian* and the first integrals whose existence is obstructed are *meromorphic* ([27, 28, 29, 30] and a long assorted array of references derived therefrom, including [34, 35]); there is further, as of yet unfinished study in the direction of situations where the system is not Hamiltonian ([16, 19, 20]).

(ii) on the other, the *dual* system (LVE)\* has jets of formal first integrals among its solutions ([3, 34]), which entails the possibility to furnish first integrals instead of finding obstructions to them; the advantage of this second approach is that the original system need not be Hamiltonian, and the formal first integrals, if convergent need not be meromorphic. The only difficulties in this case are computational in nature, namely in the context of resummation techniques.
# Formal first integrals and higher variational equations

Sergi Simon

School of Mathematics and Physics, University of Portsmouth
Lion Gate Bldg, Lion Terrace Portsmouth PO1 3HF, UK
sergi.simon@port.ac.uk

December 3, 2024

## Abstract

The question of how Algebra can be used to solve dynamical systems and characterize chaos was first posed in a fertile mathematical context by Ziglin, Morales, Ramis and Simó using differential Galois theory. Their study was aimed at first-order, later higher-order, variational equations of Hamiltonian systems. Recent work by this author formalized a compact yet comprehensive expression of higher-order variationals as one infinite linear system, thereby simplifying the approach. More importantly, the dual of this linear system contains all information relevant to first integrals, regardless of whether the original system is Hamiltonian. This applicability to formal calculation of conserved quantities is the centerpiece of this paper, following an introduction to the requisite context. Three important examples, namely particular cases of Dixon's system, the SIR epidemiological model with vital dynamics and the Van der Pol oscillator, are tackled, and explicit convergent first integrals are provided for the first two.

**Keywords.** Integrability, Ziglin-Morales-Ramis-Simó theory, formal calculus, chaos, Dixon's system, SIR epidemiological model.

**2000 Mathematics Subject Classification:** 34A05, 37C79, 37J99, 37J30, 34M15, 34C28, 37C10, 15A69, 16W60, 13F25, 37N25.

## Contents

1 Introduction 2

1.1 The algebraic study of integrability 2

1.1.1 Basics 2

1.2 Symmetric products, powers and exponentials 3

1.3 Application to power series 6

1.4 Higher-order variational equations 7

1.5 First integrals and higher-order variational equations 8

2 Jet filtering methods and computation of formal first integrals 10

2.1 Preliminaries 10

2.2 Progressive filtering and computation of formal first integrals 11

2.3 Filtering with a single infinite matrix product 12

2.3.1 Statement of the main result on infinite matrix filtering 12

2.3.2 Proof of Theorem 2.3.1 13

2.3.3 Further comments 16
In the present work we will exploit the second item (ii), namely by describing a method to automatically produce Taylor terms of formal first integrals by way of an automatic, easily recursified sequence of quadratures. First we recount the minimal background exposition necessary, then we present the main results in §2 and finally we apply it to two simple examples to test its accuracy and usefulness as well as make novel statements about the integrability of the examples themselves.

# 1.1 The algebraic study of integrability

## 1.1.1 Basics

Let $\psi(t, \cdot)$ be the flow and $\phi(t) = \psi(t, x)$ a particular solution of a given autonomous system

$$ \dot{z} = X(z), \quad X: \mathbb{C}^m \to \mathbb{C}^m, \qquad (1) $$

respectively. The **variational system** of (1) along $\phi$ has $\frac{\partial}{\partial z} \psi(t, \phi)$ as a fundamental matrix:

$$ \dot{Y} = A_1 Y, \quad A_1(t) := \left. \frac{\partial X}{\partial z} \right|_{z=\phi(t)} \in \mathrm{Mat}_n(K), \qquad (\mathrm{VE}_\phi) $$

$K = \mathbb{C}(\phi)$ being the smallest differential field ([36, Def. 1.1]) containing all entries of $\phi(t)$. $\frac{\partial^k}{\partial z^k} \psi(t, \phi)$ are multilinear $k$-forms appearing in the Taylor expansion of the flow along $\phi$:

$$ \psi(t, z) = \psi(t, \phi) + \sum_{k=1}^{\infty} \frac{1}{k!} \frac{\partial^k \psi(t, \phi)}{\partial z^k} \{z - \phi\}^k; \qquad (2) $$

$\partial_z^k \psi(t, \phi)$ also satisfy an echeloned set of systems, depending on the previous $k-1$ partial derivatives and usually called **order-k variational equations** $\mathrm{VE}_\phi^k$ ([30, p. 859-861], [34, Corollary 3]). Thus we have, (1) given, a linear system $\mathrm{VE}_\phi =: \mathrm{VE}_\phi^1 =: \mathrm{LVE}_\phi^1$ and a family of non-linear systems $\{\mathrm{VE}_\phi^k\}_{k \ge 2}$.

[34] presented an explicit linearized version $\mathrm{LVE}_\phi^k$, $k \ge 1$, by means of symmetric products $\odot$ of finite and infinite matrices based on already-existing definitions by Bekbaev, e.g. [4]. This was done in preparation for the Ziglin-Morales-Ramis-Simó (ZMRS) theoretical framework based on monodromy and differential Galois groups [27, 36, 38], but has other consequences as well. More specifically, our outcomes in [34] have two applications for system (1), Hamiltonian or not:

* full structure of $\mathrm{VE}_\phi^k$ and $\mathrm{LVE}_\phi^k$, i.e. recovering the flow, which underlies the ZMRS theoretical corpus in practicality, albeit as a tool rather than as a goal;

* a byproduct is the full structure of dual systems $(\mathrm{LVE}_\phi^k)^*$, i.e. recovering formal first integrals of (1) in ways which simplified earlier results in [3] significantly.

As said in the introduction, results in the present paper are based on the second of these applications. For examples of the first application, see [34, §6] or the bulk of [35]. See also [22, 23] for examples where the non-linearized $\mathrm{VE}_k$ were used.

# 1.2 Symmetric products, powers and exponentials

**Notation 1.2.1.** The conventions listed below were already introduced in [3, 34, 35]:

1. Multi-index modulus, arithmetic, order and lexicographic order: for $\mathbf{i} = (i_1, \dots, i_n) \in \mathbb{Z}_{\ge 0}^n$, $\mathbf{i} = |\mathbf{i}| := \sum_k i_k$; addition and subtraction are defined entrywise as usual; $\mathbf{i} \le \mathbf{j}$ means $i_k \le j_k$ for every $k \ge 1$; $\mathbf{i} <_{\mathrm{lex}} \mathbf{j}$ if $i_1 = j_1, \dots, i_{k-1} = j_{k-1}$ and $i_k < j_k$ for some $k \ge 1$.
2. Whenever such derivation is possible, we define the lexicographically sifted differential of $F(z_1, \dots, z_m)$ of order $m$ as the row vector

$$F^{(m)}(z) := \mathrm{lex} \left( \frac{\partial^m F}{\partial z_1^{i_1} \dots \partial z_n^{i_n}}(z) \right), \qquad (3)$$

where $|i| = m$ and entries are ordered as per $<_{lex}$ on multi-indices. For instance, for $n = 2$ the first two differentials would be

$$F^{(1)} = \left( \frac{\partial F}{\partial z_1}, \frac{\partial F}{\partial z_2}, \frac{\partial F}{\partial z_3} \right), \quad F^{(2)} = \left( \frac{\partial^2 F}{\partial z_1^2}, \frac{\partial^2 F}{\partial z_1 z_2}, \frac{\partial^2 F}{\partial z_1 z_3}, \frac{\partial^2 F}{\partial z_2^2}, \frac{\partial^2 F}{\partial z_2 z_3}, \frac{\partial^2 F}{\partial z_3^2} \right).$$

3. We define $d_{n,k} := \binom{n+k-1}{n-1}$, $D_{n,k} := \sum_{i=1}^k d_{n,i}$. It is easy to check there are $d_{n,k}$ $k$-plexs of integers in $\{1, \dots, n\}$, and just as many homogeneous monomials of degree $n$ in $k$ variables.

4. Given integers $k_1, \dots, k_n \ge 0$, we define the usual multinomial coefficient as

$$\binom{k_1 + \dots + k_n}{k_1, \dots, k_n} := \binom{k_1 + \dots + k_n}{k} := \frac{(k_1 + \dots + k_n)!}{k_1! k_2! \dots k_n!}.$$

For a multi-index $k \in \mathbb{Z}_{\ge 0}^n$, define $k! := k_1! \cdots k_n!$. For any two such $k, j$, we define

$$\binom{k}{p} := \frac{k_1! k_2! \cdots k_n!}{p_1! p_2! \cdots p_n! (k_1 - p_1)! (k_2 - p_2)! \cdots (k_n - p_n)!} = \binom{k_1}{p_1} \binom{k_2}{p_2} \cdots \binom{k_n}{p_n}, \quad (4)$$

and the multi-index counterpart to the multinomial, $\binom{k_1+\dots+k_m}{k_1, \dots, k_m} := \frac{(k_1+\dots+k_m)!}{k_1! k_2! \dots k_m!}$.

The compact formulation called for by Notation 1.2.1 (3) was achieved in [34] through an operation $\odot$ that had already been defined by Bekbaev (e.g. [4, 5]) and was systematized using basic categorical properties of the tensor product. Let $K$ be a field and $V$ a $K$-vector space. Let $\mathrm{Sym}^r V$ be the symmetric power of $V$. We write $w_1 \odot w_2$ for equivalence classes of tensor products of these vectors. Hence, product $\odot$ operates exactly like products of homogeneous polynomials in several variables.

**Notation 1.2.2.** When dealing with matrix sets, we will use super-indices and subindices:

1. The space of $(i, j)$-matrices $\mathrm{Mat}_{m,n}^{i,j}(K)$ is defined equivalently as the set of $d_{m,i} \times d_{n,j}$ matrices with entries in $K$, or vector space $\mathrm{Hom}_K(\mathrm{Sym}^j K^m; \mathrm{Sym}^i K^n)$.

2. It is clear from the above that $\mathrm{Mat}_n^{0,0}(K)$ is the set of all scalars $\alpha \in K$ and $\mathrm{Mat}_n^{0,k}(K)$ (resp. $\mathrm{Mat}_n^{k,0}(K)$) is made up of all row (resp. column) vectors whose entries are indexed by $d_{n,k}$ lexicographically ordered $k$-tuples.

3. Reference to $K$ may be dropped and notation may be abridged if dimensions are repeated or trivial, e.g. $\mathrm{Mat}_n^{i,j} := \mathrm{Mat}_{n,n}^{i,j}$, $\mathrm{Mat}_{m,n}^{i,j} := \mathrm{Mat}_{m,n}^{i,j}$, $\mathrm{Mat}_n := \mathrm{Mat}_n^1$, etcetera.

4. $\mathrm{Mat}^{n,m}(K)$ denotes the set of block matrices $A = (A_{i,j})_{i,j \ge 0}$ with $A_{i,j} : \mathrm{Sym}^i K^m \to \mathrm{Sym}^j K^n$, hence $A_{i,j} \in \mathrm{M}_{d_{n,i} \times d_{m,j}}(K) = \mathrm{Mat}_{n,m}^{i,j}(K)$:

$$A = \begin{pmatrix}
\cdots & \vdots & \vdots & \vdots \\
\cdots & A_{2,2} & A_{2,1} & \leftarrow A_{2,0} \\
\cdots & A_{1,2} & A_{1,1} & \leftarrow A_{1,0} \\
\cdots & A_{0,2} & A_{0,1} & \leftarrow A_{0,0}
\end{pmatrix}$$
Corollary 3]). Thus we have, (1) given, a linear system $VE_{\phi} =: VE_{\phi}^{1} =: LVE_{\phi}^{1}$ and a family of non-linear systems $\{VE_{\phi}^{k}\}_{k \ge 2}$.

[34] presented an explicit linearized version $LVE_{\phi}^{k}$, $k \ge 1$, by means of symmetric products $\odot$ of finite and infinite matrices based on already-existing definitions by Bekbaev, e.g. [4]. This was done in preparation for the Ziglin-Morales-Ramis-Simó (ZMRS) theoretical framework based on monodromy and differential Galois groups [27, 37, 39], but has other consequences as well. More specifically, our outcomes in [34] have two applications for system (1), Hamiltonian or not:

* full structure of $VE_{\phi}^{k}$ and $LVE_{\phi}^{k}$, i.e. recovering the flow, which underlies the ZMRS theoretical corpus in practicality, albeit as a tool rather than as a goal;

* a byproduct is the full structure of dual systems $(LVE_{\phi}^{k})^{*}$, i.e. recovering formal first integrals of (1) in ways which simplified earlier results in [3] significantly.

As said in the introduction, results in the present paper are based on the second of these applications. For examples of the first application, see [34, §6] or the bulk of [35]. See also [22, 23] for examples where the non-linearized $VE_{k}$ were used.

## 1.2 Symmetric products, powers and exponentials

**Notation 1.2.1.** The conventions listed below were already introduced in [3, 34, 35]:

1. Multi-index modulus, arithmetic, order and lexicographic order: for $\mathbf{i} = (i_1, \dots, i_n) \in \mathbb{Z}_{\ge 0}^n$, $i = |\mathbf{i}| := \sum_k i_k$; addition and subtraction are defined entrywise as usual; $\mathbf{i} \le \mathbf{j}$ means $i_k \le j_k$ for every $k \ge 1$; $\mathbf{i} <_{\text{lex}} \mathbf{j}$ if $i_1 = j_1, \dots, i_{k-1} = j_{k-1}$ and $i_k < j_k$ for some $k \ge 1$.

2. Whenever such derivation is possible, we define the lexicographically sifted differential of $F(z_1, \dots, z_m)$ of order $m$ as the row vector

$$F^{(m)}(z) := \mathrm{lex} \left( \frac{\partial^m F}{\partial z_1^{i_1} \dots \partial z_n^{i_n}}(z) \right), \quad (3)$$

where $|\mathbf{i}| = m$ and entries are ordered as per $<_{\text{lex}}$ on multi-indices. For instance, for $n = 2$ the first two differentials would be

$$F^{(1)} = \left( \frac{\partial F}{\partial z_1}, \frac{\partial F}{\partial z_2}, \frac{\partial F}{\partial z_3} \right), \quad F^{(2)} = \left( \frac{\partial^2 F}{\partial z_1^2}, \frac{\partial^2 F}{\partial z_1 z_2}, \frac{\partial^2 F}{\partial z_1 z_3}, \frac{\partial^2 F}{\partial z_2^2}, \frac{\partial^2 F}{\partial z_2 z_3}, \frac{\partial^2 F}{\partial z_3^2} \right).$$

3. We define $d_{n,k} := \binom{n+k-1}{n-1}$, $D_{n,k} := \sum_{i=1}^k d_{n,i}$. It is easy to check there are $d_{n,k}$ k-ples of integers in $\{1, \dots, n\}$, and just as many homogeneous monomials of degree $n$ in $k$ variables.

4. Given integers $k_1, \dots, k_n \ge 0$, we define the usual multinomial coefficient as

$$\binom{k_1 + \cdots + k_n}{k_1, \dots, k_n} := \binom{k_1 + \cdots + k_n}{\mathbf{k}} := \frac{(k_1 + \cdots + k_n)!}{k_1! k_2! \cdots k_n!}.$$

For a multi-index $\mathbf{k} \in \mathbb{Z}_{\ge 0}^n$, define $\mathbf{k}! := k_1! \cdots k_n!$. For any two such $\mathbf{k}, \mathbf{j}$, we define

$$\binom{\mathbf{k}}{\mathbf{p}} := \frac{k_1! k_2! \cdots k_n!}{p_1! p_2! \cdots p_n! (k_1 - p_1)! (k_2 - p_2)! \cdots (k_n - p_n)!} = \binom{k_1}{p_1} \binom{k_2}{p_2} \cdots \binom{k_n}{p_n}, \quad (4)$$

and the multi-index counterpart to the multinomial, $\binom{k_1+\cdots+k_m}{k_1, \dots, k_m} := \frac{(k_1+\cdots+k_m)!}{k_1! k_2! \cdots k_m!}$.
We write Mat := Mat<sup>n,n</sup> if n is unambiguous. Conversely, Mat<sup>i,j</sup><sub>n,m</sub> is embedded in Mat<sup>n,m</sup> by identifying every matrix A<sub>i,j</sub> with an element of Mat<sup>n,m</sup> equal to 0 save for block A<sub>i,j</sub>.

**Definition 1.2.3** (Symmetric products of finite and infinite matrices). [4, 34]

(i) Let $A \in \mathrm{Mat}_{m,n}^{i_1,j_1}(K)$, $B \in \mathrm{Mat}_{m,n}^{i_2,j_2}(K)$. Given any multi-index $\mathbf{k} = (k_1, \dots, k_n) \in \mathbb{Z}_{\ge 0}^n$ such that $|\mathbf{k}| = k_1 + \dots + k_n = j_1 + j_2$, define $C := A \odot B \in \mathrm{Mat}_{m,n}^{i_1+i_2, j_1+j_2}$ by

$$C(e_1^{\odot k_1} \cdots e_n^{\odot k_n}) = \frac{1}{(j_1+j_2)} \sum_{\mathbf{p}} \binom{\mathbf{k}}{j_1} A(e_1^{\odot p_1} \cdots e_n^{\odot p_n}) \odot B(e_1^{\odot k_1-p_1} \cdots e_n^{\odot k_n-p_n}), \quad (5)$$

notation abused by removing $\odot$ to reduce space, binomials as in (4) and summation taking place for specific multi-indices $\mathbf{p}$, namely those such that

$$|\mathbf{p}| = j_1 \quad \text{and} \quad 0 \le p_i \le k_i, \quad i = 1, \dots, n.$$

(ii) For any $A, B \in \mathrm{Mat}^{n,m}(K)$, define $A \odot B = C \in \mathrm{Mat}^{n,m}(K)$ by

$$C = (C_{i,j})_{i,j \ge 0}, \quad C_{i,j} = \sum_{0 \le i_1 \le i, 0 \le j_1 \le j} \binom{j}{j_1} A_{i_1, j_1} \odot B_{i-i_1, j-j_1}. \quad (6)$$

Same as always, $\odot^k$ will stand for powers built with this product.

The following is a mere exercise in induction:

**Lemma 1.2.4.** Defining $\odot_{i=1}^r A_i$ recursively by $(\odot_{i=1}^{r-1} A_i) \odot A_r$ with $A_i \in \mathrm{Mat}_{m,n}^{k_i, j_i}$,

$$(A_1 \odot \cdots \odot A_r) e^{\odot \mathbf{k}} = \frac{1}{\binom{j_1+\cdots+j_r}{j_1, j_2, \ldots, j_r}} \sum_{\mathbf{p}_1, \ldots, \mathbf{p}_r} \binom{\mathbf{k}}{\mathbf{p}_1, \ldots, \mathbf{p}_r} \odot_{i=1}^r A_i e^{\odot \mathbf{p}_i}, \quad (7)$$

if $|\mathbf{k}| = j_1 + \cdots + j_r$, sums obviously taken for $\mathbf{p}_1 + \cdots + \mathbf{p}_r = \mathbf{k}$ and $|\mathbf{p}_i| = j_i$, for every $i = 1, \dots, r$. $\square$

The following is straightforward and has already been seen e.g. in [4, 5, 34], or can be easily derived therefrom:

**Proposition 1.2.5.** For any $A, B, C$, and whenever products make sense,

$$\text{a)} A \odot B = B \odot A.$$

$$\text{b)} (A+B) \odot C = A \odot C + B \odot C.$$

$$\text{c)} (A \odot B) \odot C = A \odot (B \odot C).$$

$$\text{d)} (\alpha A) \odot B = \alpha (A \odot B) \text{ for every } \alpha \in K.$$

$$\text{e)} \text{ If } A \text{ is square and invertible, then } (A^{-1})^{\odot k} = (A^{\odot k})^{-1}.$$

$$\text{f)} A \odot B = 0 \text{ if and only if } A = 0 \text{ or } B = 0.$$

$$\text{g)} \text{ If } A \text{ is a square (1,1)-matrix, then } A\mathbf{v}_1 \odot A\mathbf{v}_2 \odot \cdots \odot A\mathbf{v}_m = A^{\odot m}\mathbf{v}_1 \odot \cdots \odot \mathbf{v}_m.$$

$$\text{h)} \text{ If } \mathbf{v} \text{ is a column vector, then } (A \odot B) \mathbf{v}^{\odot(p+q)} = (A\mathbf{v}^{\odot p}) \odot (B\mathbf{v}^{\odot q}), p, q \in \mathbb{Z}_{\ge 0}.$$
The compact formulation called for by Notation 1.2.1 (3) was achieved in [34] through an operation $\odot$ that had already been defined by Bekbaev (e.g. [4, 5]) and was systematized using basic categorical properties of the tensor product. Let $K$ be a field and $V$ a $K$-vector space. Let $\mathrm{Sym}^r V$ be the symmetric power of $V$. We write $w_1 \odot w_2$ for equivalence classes of tensor products of these vectors. Hence, product $\odot$ operates exactly like products of homogeneous polynomials in several variables.

**Notation 1.2.2.** When dealing with matrix sets, we will use super-indices and subindices:

1. The space of $(i, j)$-matrices $\mathrm{Mat}_{m,n}^{i,j}(K)$ is defined equivalently as the set of $d_{m,i} \times d_{n,j}$ matrices with entries in $K$, or vector space $\mathrm{Hom}_K(\mathrm{Sym}^j K^m; \mathrm{Sym}^i K^n)$.

2. It is clear from the above that $\mathrm{Mat}_n^{0,0}(K)$ is the set of all scalars $\alpha \in K$ and $\mathrm{Mat}_n^{0,k}(K)$ (resp. $\mathrm{Mat}_n^{k,0}(K)$) is made up of all row (resp. column) vectors whose entries are indexed by $d_{n,k}$ lexicographically ordered $k$-tuples.

3. Reference to $K$ may be dropped and notation may be abridged if dimensions are repeated or trivial, e.g. $\mathrm{Mat}_n^{i,j} := \mathrm{Mat}_{n,n}^{i,j}$, $\mathrm{Mat}_{m,n}^i := \mathrm{Mat}_{m,n}^{i,i}$, $\mathrm{Mat}_n := \mathrm{Mat}_n^1$, etcetera.

4. $\mathrm{Mat}^{n,m}(K)$ denotes the set of block matrices $A = (A_{i,j})_{i,j \ge 0}$ with $A_{i,j} : \mathrm{Sym}^i K^m \to \mathrm{Sym}^j K^n$, hence $A_{i,j} \in \mathrm{M}_{d_{n,i} \times d_{m,j}}(K) = \mathrm{Mat}_{n,m}^{i,j}(K)$:

$$A = \begin{pmatrix}
\cdots & \vdots & \vdots & \vdots \\
\cdots & A_{2,2} & A_{2,1} & \leftarrow A_{2,0} \\
\cdots & A_{1,2} & A_{1,1} & \leftarrow A_{1,0} \\
\cdots & A_{0,2} & A_{0,1} & \leftarrow A_{0,0}
\end{pmatrix}$$

We write $\mathrm{Mat} := \mathrm{Mat}^{n,n}$ if $n$ is unambiguous. Conversely, $\mathrm{Mat}_{n,m}^{i,j}$ is embedded in $\mathrm{Mat}^{n,m}$ by identifying every matrix $A_{i,j}$ with an element of $\mathrm{Mat}^{n,m}$ equal to 0 save for block $A_{i,j}$.

**Definition 1.2.3** (Symmetric products of finite and infinite matrices). [4, 34]

(i) Let $A \in \mathrm{Mat}_{m,n}^{i_1,j_1}(K)$, $B \in \mathrm{Mat}_{m,n}^{i_2,j_2}(K)$. Given any multi-index $\mathbf{k} = (k_1, \dots, k_n) \in \mathbb{Z}_{\ge 0}^n$ such that $|\mathbf{k}| = k_1 + \dots + k_n = j_1 + j_2$, define $C := A \odot B \in \mathrm{Mat}_{m,n}^{i_1+i_2, j_1+j_2}$ by

$$C(e_1^{\odot k_1} \cdots e_n^{\odot k_n}) = \frac{1}{(j_1+j_2)} \sum_{\mathbf{p}} \binom{\mathbf{k}}{\mathbf{p}} A(e_1^{\odot p_1} \cdots e_n^{\odot p_n}) \odot B(e_1^{\odot k_1-p_1} \cdots e_n^{\odot k_n-p_n}), \quad (5)$$

notation abused by removing $\odot$ to reduce space, binomials as in (4) and summation taking place for specific multi-indices $\mathbf{p}$, namely those such that

$$|\mathbf{p}| = j_1 \quad \text{and} \quad 0 \le p_i \le k_i, \quad i = 1, \dots, n.$$

(ii) For any $A, B \in \mathrm{Mat}^{n,m}(K)$, define $A \odot B = C \in \mathrm{Mat}^{n,m}(K)$ by

$$C = (C_{i,j})_{i,j \ge 0}, \quad C_{i,j} = \sum_{0 \le i_1 \le i, \ 0 \le j_1 \le j} \binom{j}{j_1} A_{i_1, j_1} \odot B_{i-i_1, j-j_1}. \quad (6)$$

Same as always, $\odot^k$ will stand for powers built with this product.
**Lemma 1.2.6 ([34]).** (i) Given square A, B ∈ Mat<sub>n</sub><sup>k,k</sup> and matrices X<sub>i</sub> ∈ Mat<sub>n</sub><sup>k,j<sub>i</sub></sup>, i = 1, 2,

$$ (A \odot B)(X_1 \odot X_2) = \frac{1}{2} (AX_1 \odot BX_2 + BX_1 \odot AX_2), \qquad (8) $$

and in general for any square A<sub>1</sub>, ..., A<sub>m</sub> ∈ Mat<sub>n</sub><sup>k,k</sup> and any X<sub>i</sub> ∈ Mat<sub>n</sub><sup>k,j<sub>i</sub></sup>, i = 1, ..., m,

$$ \left( \bigodot_{i=1}^{m} A_i \right) \left( \bigodot_{i=1}^{m} X_i \right) = \frac{1}{m!} \sum_{\sigma \in \mathfrak{S}_k} \bigodot_{i=1}^{m} A_{\sigma(i)} X_i. \quad \square \qquad (9) $$

(ii) Given A ∈ Mat<sub>n</sub><sup>1,j</sup> and X<sub>1</sub>, ..., X<sub>m</sub> such that X<sub>i</sub> ∈ Mat<sub>n</sub><sup>1,q<sub>i</sub></sup>, 1 ≤ j ≤ m,

$$ \binom{m}{j} (A \odot \mathrm{Id}_{n}^{\odot m-j}) \bigodot_{i=1}^{m} X_i = \sum_{1 \le i_1 < \cdots < i_j \le m} [A (X_{i_1} \odot \cdots \odot X_{i_j})] \odot \bigodot_{s \ne i_1, \ldots, i_j} X_s. \quad \square \qquad (10) $$

(iii) Given a square matrix A ∈ Mat<sub>n</sub><sup>1,1</sup> and X<sub>1</sub>, ..., X<sub>m</sub> such that X<sub>i</sub> ∈ Mat<sub>n</sub><sup>1,j<sub>i</sub></sup>,

$$ (A \odot \mathrm{Id}_{n}^{\odot m-1}) \left( \bigodot_{i=1}^{m} X_i \right) = \frac{1}{m} \sum_{i=1}^{m} (AX_i) \odot (X_1 \odot \cdots \odot \widehat{X_i} \odot \cdots \odot X_m). \qquad (11) $$

(iv) Given square matrix X ∈ Mat<sub>n</sub><sup>k-1,k-1</sup> and vector v ∈ K<sup>n</sup>, we have, for each i = 1, ..., n,

$$ (X \odot e_i^T) (v \odot \mathrm{Id}_n^{\odot k-1}) = \frac{k-1}{k} X (v \odot \mathrm{Id}_n^{\odot k-2}) \odot e_i^T + \frac{v_i}{k} X \qquad (12) $$

**Lemma 1.2.7.** Let (K, ∂) be a differential field.

(i) For any given X ∈ Mat<sub>n</sub><sup>k<sub>1</sub>,j<sub>1</sub></sup> (K) and Y ∈ Mat<sub>n</sub><sup>k<sub>2</sub>,j<sub>2</sub></sup> (K),

$$ \partial(X \odot Y) = \partial(X) \odot Y + X \odot \partial(Y). \quad \square \qquad (13) $$

(ii) If Y is a square n × n matrix having entries in K and ∂Y = AY, then

$$ \partial \mathrm{Sym}^k Y = k (A \odot \mathrm{Sym}^{k-1} (\mathrm{Id}_n)) \mathrm{Sym}^k Y. \qquad (14) $$

(iii) If X ∈ Mat<sub>n</sub><sup>1,j<sub>1</sub></sup> and Y ∈ Mat<sub>n</sub><sup>1,j<sub>2</sub></sup> satisfy systems ∂X = AX + B<sub>1</sub> and ∂Y = AY + B<sub>2</sub> with A ∈ Mat<sub>n</sub><sup>1,1</sup>, B<sub>i</sub> ∈ Mat<sub>n</sub><sup>1,j<sub>i</sub></sup>, then symmetric product X ⊙ Y satisfies linear system

$$ \partial(X \odot Y) = 2 (A \odot \mathrm{Id}_{d_{n,k}}) (X \odot Y) + (B_1 \odot Y + B_2 \odot X). \qquad (15) $$

(iv) If ∂X<sub>i</sub> = AX<sub>i</sub> + B<sub>i</sub>, i = 1, ..., m, with X<sub>i</sub>, B<sub>i</sub> ∈ Mat<sub>n</sub><sup>1,j<sub>i</sub></sup>, A ∈ Mat<sub>n</sub><sup>1,1</sup> then

$$ \partial \bigodot_{i=1}^{m} X_i = m \left( A \odot \mathrm{Id}_{d_{n,k}}^{\odot m-1} \right) \bigodot_{i=1}^{m} X_i + \sum_{i=1}^{m} B_i \odot \bigodot_{j \ne i} X_j. \quad \square \qquad (16) $$

**Definition 1.2.8.** (See also [5]) for every matrix A ∈ Mat<sub>n,m</sub> we define the formal power series

$$ \exp_{\odot} A := 1 + A^{\odot 1} + \frac{1}{2} A^{\odot 2} + \cdots = \sum_{i=0}^{\infty} \frac{1}{i!} A^{\odot i}. $$

Whenever A = 0 save for a finite distinguished submatrix A<sub>j,k</sub>, the abuse of notation exp<sub>⊗</sub> A<sub>j,k</sub> = exp<sub>⊗</sub> A will be customary.
The following is a mere exercise in induction:

**Lemma 1.2.4.** Defining $\odot_{i=1}^r A_i$ recursively by $(\odot_{i=1}^{r-1} A_i) \odot A_r$ with $A_i \in \mathrm{Mat}_{m,n}^{k_i,j_i}$,

$$ (A_1 \odot \cdots \odot A_r) e^{\odot k} = \frac{1}{\binom{j_1+\cdots+j_r}{j_1, j_2, \ldots, j_r}} \sum_{\mathbf{p}_1, \ldots, \mathbf{p}_r} \left( \mathbf{k} \cdot \mathbf{p}_1, \ldots, \mathbf{p}_r \right) \odot \sum_{i=1}^r A_i e^{\odot \mathbf{p}_i}, \quad (7) $$

if $|\mathbf{k}| = j_1 + \cdots + j_r$, sums obviously taken for $\mathbf{p}_1 + \cdots + \mathbf{p}_r = \mathbf{k}$ and $|\mathbf{p}_i| = j_i$, for every $i = 1, \ldots, r$. $\square$

The following is straightforward and has already been seen e.g. in [4, 5, 34], or can be easily derived therefrom:

**Proposition 1.2.5.** For any A, B, C, and whenever products make sense,

a) $A \odot B = B \odot A$.

b) $(A+B) \odot C = A \odot C + B \odot C$.

c) $(A \odot B) \odot C = A \odot (B \odot C)$.

d) $(\alpha A) \odot B = \alpha (A \odot B)$ for every $\alpha \in K$.

e) If A is square and invertible, then $(A^{-1})^{\odot k} = (A^{\odot k})^{-1}$.

f) $A \odot B = 0$ if and only if $A = 0$ or $B = 0$.

g) If A is a square (1,1)-matrix, then $A\mathbf{v}_1 \odot A\mathbf{v}_2 \odot \cdots \odot A\mathbf{v}_m = A^{\odot m}\mathbf{v}_1 \odot \cdots \odot \mathbf{v}_m$.

h) If $\mathbf{v}$ is a column vector, then $(A \odot B) \mathbf{v}^{\odot(p+q)} = (A\mathbf{v}^{\odot p}) \odot (B\mathbf{v}^{\odot q})$, $p, q \in \mathbb{Z}_{\ge 0}$.

**Lemma 1.2.6 ([34]).** (i) Given square A, B $\in \mathrm{Mat}_n^{k,k}$ and matrices $X_i \in \mathrm{Mat}_n^{k,j_i}$, $i = 1, 2$,

$$ (A \odot B)(X_1 \odot X_2) = \frac{1}{2}(AX_1 \odot BX_2 + BX_1 \odot AX_2), \quad (8) $$

and in general for any square $A_1, \ldots, A_m \in \mathrm{Mat}_n^{k,k}$ and any $X_i \in \mathrm{Mat}_n^{k,j_i}$, $i = 1, \ldots, m$,

$$ \left( \bigodot_{i=1}^{m} A_i \right) \left( \bigodot_{i=1}^{m} X_i \right) = \frac{1}{m!} \sum_{\sigma \in \mathfrak{S}_k} \bigodot_{i=1}^{m} A_{\sigma(i)} X_i. \quad \square \quad (9) $$

(ii) Given $A \in \mathrm{Mat}_n^{1,j}$ and $X_1, \ldots, X_m$ such that $X_i \in \mathrm{Mat}_n^{1,q_i}$, $1 \le j \le m$,

$$ \binom{m}{j} (A \odot \mathrm{Id}_n^{\odot m-j}) \bigodot_{i=1}^{m} X_i = \sum_{1 \le i_1 < \cdots < i_j \le m} [A (X_{i_1} \odot \cdots \odot X_{i_j})] \odot \bigodot_{s \ne i_1, \ldots, i_j} X_s. \quad \square \quad (10) $$

(iii) Given a square matrix $A \in \mathrm{Mat}_n^{1,1}$ and $X_1, \ldots, X_m$ such that $X_i \in \mathrm{Mat}_n^{1,j_i}$,

$$ (A \odot \mathrm{Id}_n^{\odot m-1}) \left( \bigodot_{i=1}^{m} X_i \right) = \frac{1}{m} \sum_{i=1}^{m} (AX_i) \odot (X_1 \odot \cdots \odot \widehat{X_i} \odot \cdots \odot X_m). \quad (11) $$

(iv) Given square matrix $X \in \mathrm{Mat}_n^{k-1,k-1}$ and vector $\mathbf{v} \in K^n$, we have, for each $i = 1, \ldots, n$,

$$ (X \odot e_i^T) (\mathbf{v} \odot \mathrm{Id}_n^{\odot k-1}) = \frac{k-1}{k} X (\mathbf{v} \odot \mathrm{Id}_n^{\odot k-2}) \odot e_i^T + \frac{v_i}{k} X \quad (12) $$
**Lemma 1.2.9.**

a) For every two $A, B \in \mathrm{Mat}^{n,m}$, $\exp_{\odot}(A+B) = \exp_{\odot} A \odot \exp_{\odot} B$.

b) For every $Y \in \mathrm{Mat}^{n,m}$ and any derivation $\partial: K \to K$, $\partial \exp_{\odot} Y = (\partial Y) \odot \exp_{\odot} Y$.

c) ([4, Corollary 3]) Given square matrices $A, B \in \mathrm{Mat}_n^{1,1}$, $\exp_{\odot} AB = \exp_{\odot} A \exp_{\odot} B$.

d) In particular, for every invertible square $A \in \mathrm{Mat}_n^{1,1}$, $\exp_{\odot} A^{-1} = (\exp_{\odot} A)^{-1}$. $\square$

For examples and properties, see [34, §3.1].

**1.3 Application to power series**

**Lemma 1.3.1** ([34, Lemma 3.7]). If $F = F_1 \times \cdots \times F_m$ is a vector power series, adequate $M_F^{1,i} \in \mathrm{Mat}_{m,n}^{1,i}(K)$ render

$$F(\boldsymbol{x}) = M_F \exp_{\odot} X \quad \text{where } M_F := \left( \begin{array}{ccc} \cdots & M_F^{1,2} & M_F^{1,1} \\ \cdots & 0 & 0 \\ \cdots & 0 & 0 \end{array} \right) \in \mathrm{Mat}^{m,n}.$$

Following Definition 1.2.8, write $F(\boldsymbol{x}) = M_F \exp_{\odot} \boldsymbol{x}$ if it poses no clarity issue. $\square$

Thus every formal power series $F \in K[[\boldsymbol{x}]]$, $\boldsymbol{x} = (x_1, \dots, x_n)$ can be expressed in the form $M_F \exp_{\odot} \boldsymbol{x}$, where

$$M_F := \left( \begin{array}{ccc} \cdots & M_F^{1,2} & M_F^{1,1} \\ \cdots & 0 & 0 \\ \cdots & 0 & 0 \end{array} \right) \in \mathrm{Mat}^{1,n}(K), \quad X := \left( \begin{array}{c|c} 0 & \boldsymbol{x} \\ \hline 0 & 0 \end{array} \right). \\ \\ M_F = J_F + M_F^{1,0} := \left( \begin{array}{c|c} \cdots & M_F^{1,2} \\ \cdots & 0 \end{array} \right) + \left( \begin{array}{c|c} 0 & M_F^{1,0} \\ \hline 0 & 0 \end{array} \right). \qquad (17)$$

See [34, §3.2] for more details.

**1.4 Higher-order variational equations**

See [34] for proofs and further elaboration on everything summarized below:

**Notation 1.4.1.** For every set of indices $1 \le i_1 \le \cdots \le i_r$ such that $\sum_{j=1}^r i_j = k$, $c_{i_1, \dots, i_r}^k$ is defined as the amount of totally ordered partitions of a set of $k$ elements among subsets of sizes $i_1, \dots, i_r$:

$$c_{i_1, \dots, i_r} = c_{i_1, \dots, i_r}^k = \frac{\binom{k}{i_1 i_2 \dots i_r}}{n_1! \dots n_m!}, \quad \left\{ \begin{array}{l} (i_1, \dots, i_r) = (k_1 \stackrel{n_1}{\cdot} k_1, \dots, k_m \stackrel{n_m}{\cdot} k_m), \\ 1 \le k_1 < k_2 < \dots < k_m. \end{array} \right. \qquad (18)$$

**Proposition 1.4.2.** Let $K = \mathbb{C}(\phi)$ and $A_k$, $Y_k$ be $\partial_k X(\phi)$, $\partial_z^k \psi(t, \phi)$ minus crossed terms; let $A = J_X^\phi$, $Y = J_\phi$ be the derivative jets for $X$ and $\psi$ at $\phi$, with $A_k$, $Y_k$ as blocks. Then, if $c_i^k = \#\{\text{ordered } i_1, \dots, i_r \text{-partitions of } k \text{ elements}\}$,

$$\dot{Y}_k = \sum_{j=1}^{k} A_j \sum_{|i|=k} c_i^k Y_{i_1} \odot Y_{i_2} \odot \cdots \odot Y_{i_r}, \quad k \ge 1. \qquad (\mathrm{VE}_\phi^k)$$

**Lemma 1.4.3.** Let $Y \in \mathrm{Mat}(K)$ equal to zero outside of block row 1,k, $k \ge 1$:

$$Y := \left( \begin{array}{cccc|c} \cdots & Y_3 & Y_2 & Y_1 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{array} \right), \quad Y_i \in \mathrm{Mat}_n^{1,i}. \qquad (19)$$

Let $Z_{r,s}$, $s, r \ge 1$, be the corresponding block in $\exp_{\odot} Y$. Then,
**Lemma 1.2.7.** Let $(K, \partial)$ be a differential field.

(i) For any given $X \in \mathrm{Mat}_n^{k_1, j_1}(K)$ and $Y \in \mathrm{Mat}_n^{k_2, j_2}(K)$,
$$ \partial(X \odot Y) = \partial(X) \odot Y + X \odot \partial(Y). \quad \square \qquad (13) $$

(ii) If $Y$ is a square $n \times n$ matrix having entries in $K$ and $\partial Y = AY$, then

$$ \partial \mathrm{Sym}^k Y = k (A \odot \mathrm{Sym}^{k-1}(\mathrm{Id}_n)) \mathrm{Sym}^k Y. \qquad (14) $$

(iii) If $X \in \mathrm{Mat}_n^{1, j_1}$ and $Y \in \mathrm{Mat}_n^{1, j_2}$ satisfy systems $\partial X = AX + B_1$ and $\partial Y = AY + B_2$ with $A \in \mathrm{Mat}_n^{1, 1}$, $B_i \in \mathrm{Mat}_n^{1, j_i}$, then symmetric product $X \odot Y$ satisfies linear system

$$ \partial(X \odot Y) = 2 (A \odot \mathrm{Id}_{d_{n,k}}) (X \odot Y) + (B_1 \odot Y + B_2 \odot X). \qquad (15) $$

(iv) If $\partial X_i = AX_i + B_i$, $i = 1, \dots, m$, with $X_i, B_i \in \mathrm{Mat}_n^{1, j_i}$, $A \in \mathrm{Mat}_n^{1, 1}$ then

$$ \partial \bigcirc_{i=1}^{m} X_i = m \left( A \odot \mathrm{Id}_{d_{n,k}}^{\odot m-1} \right) \bigcirc_{i=1}^{m} X_i + \sum_{i=1}^{m} B_i \odot \bigcirc_{j \neq i} X_j. \quad \square \qquad (16) $$

**Definition 1.2.8.** (See also [5]) for every matrix $A \in \mathrm{Mat}^{n,m}$ we define the formal power series

$$ \exp_{\odot} A := 1 + A^{\odot 1} + \frac{1}{2} A^{\odot 2} + \dots = \sum_{i=0}^{\infty} \frac{1}{i!} A^{\odot i}. $$

Whenever $A = 0$ save for a finite distinguished submatrix $A_{j,k}$, the abuse of notation $\exp_{\odot} A_{j,k} = \exp_{\odot} A$ will be customary.

**Lemma 1.2.9.**

a) For every two $A, B \in \mathrm{Mat}^{n,m}$, $\exp_{\odot}(A+B) = \exp_{\odot} A \odot \exp_{\odot} B$.

b) For every $Y \in \mathrm{Mat}^{n,m}$ and any derivation $\partial: K \to K$, $\partial \exp_{\odot} Y = (\partial Y) \odot \exp_{\odot} Y$.

c) ([4, Corollary 3]) Given square matrices $A, B \in \mathrm{Mat}_n^{1, 1}$, $\exp_{\odot} AB = \exp_{\odot} A \exp_{\odot} B$.

d) In particular, for every invertible square $A \in \mathrm{Mat}_n^{1, 1}$, $\exp_{\odot} A^{-1} = (\exp_{\odot} A)^{-1}$. $\square$
For examples and properties, see [34, §3.1].

## 1.3 Application to power series

**Lemma 1.3.1** ([34, Lemma 3.7]). If $F = F_1 \times \dots \times F_m$ is a vector power series, adequate $M_F^{1, i} \in \mathrm{Mat}_{m,n}^{1, i}(K)$ render

$$ F(\boldsymbol{x}) = M_F \exp_{\odot} X \quad \text{where} \quad M_F := \left( \begin{array}{c|c} \dots & M_F^{1,2} \\ \hline \dots & M_F^{1,1} \\ \hline 0 & 0 \end{array} \right) \in \mathrm{Mat}^{m,n}. $$

Following Definition 1.2.8, write $F(\boldsymbol{x}) = M_F \exp_{\odot} \boldsymbol{x}$ if it poses no clarity issue. $\square$

Thus every formal power series $F \in K[[\boldsymbol{x}]]$, $\boldsymbol{x} = (x_1, \dots, x_n)$ can be expressed in the form $M_F \exp_{\odot} \boldsymbol{x}$, where

$$ M_F := \left( \begin{array}{c|c} \dots & M_F^{1,2} \\ \hline \dots & M_F^{1,1} \\ \hline 0 & 0 \end{array} \right) \in \mathrm{Mat}^{1,n}(K), \quad X := \left( \begin{array}{c|c} 0 & \boldsymbol{x} \\ \hline 0 & 0 \end{array} \right). \\ M_F = J_F + M_F^{1,0} := \left( \begin{array}{c|c} \dots & M_F^{1,2} \\ \hline \dots & M_F^{1,1} \\ \hline 0 & 0 \end{array} \right) + \left( \begin{array}{c|c} 0 & M_F^{1,0} \\ \hline 0 & 0 \end{array} \right). \qquad (17) $$

See [34, §3.2] for more details.
a) Row block $r$ in $\exp_{\odot} Y$ is recursively obtained in terms of row blocks 1 and $r-1$:

$$Z_{r,s} = \frac{1}{r} \sum_{j=1}^{s-r+1} \binom{s}{j} Y_j \odot Z_{r-1,s-j}. \qquad (20)$$

In particular, $Z_{r,r} = Y_1^{\odot r}$ and $Z_{r,s} = 0_{d_{n,r},d_{n,s}}$ whenever $r > s$.

b) Using Notation 1.4.1 and (18), for every $s \ge r$

$$Z_{r,s} = \sum_{i_1+\cdots+i_r=s} c_{i_1, \dots, i_r}^s Y_{i_1} \odot Y_{i_2} \odot \cdots \odot Y_{i_r}. \qquad (21)$$

**Notation 1.4.4.** $K := \mathbb{C}(\phi)$, $A_i := X^{(i)}(\phi)$, $Y_i := \text{lex} \left( \frac{\partial^i}{\partial z^i} \varphi(t, \phi) \right)$ and, per Lemma 1.4.3,

$$\Upsilon_1 = Y_1, \qquad \Upsilon_k = \left( \begin{array}{c} Z_{k,k} \\ Z_{k-1,k} \\ \vdots \\ Z_{1,k} \end{array} \right) \quad \left. \begin{array}{c} \Upsilon_{k-1} \\ \end{array} \right\} , \quad k \ge 2, \qquad (22)$$

formed by the first $k$ block rows and columns in $\Upsilon = \exp_{\odot} Y$. Define $A, Y \in \text{Mat}(K)$ as in Lemma 1.4.3 with the above $A_i, Y_i$ as blocks. Denote the canonical basis on $K^n$ (meaning the set of columns of $\text{Id}_n$) by $\{e_1, \dots, e_n\}$.

**Proposition 1.4.5** (Explicit version of $\text{LVE}_{\phi}^k$). Still following Notation 1.4.4, the infinite system

$$\boxed{\dot{X} = A_{\text{LVE}_{\phi}} X,} \qquad A_{\text{LVE}_{\phi}} := A \odot \exp_{\odot} \text{Id}_n, \qquad (\text{LVE}_{\phi})$$

has $\Upsilon := \exp_{\odot} Y$ as a solution matrix. Hence, for every $k \ge 1$,

a) the lower-triangular recursive $D_{n,k} \times D_{n,k}$ form for $\text{LVE}_{\phi}^k$ is $\dot{Y} = A_{\text{LVE}_{\phi}^k} Y$, its system matrix being obtained from the first $k$ row and column blocks of $A_{\text{LVE}_{\phi}}$:

$$A_{\text{LVE}_{\phi}^k} = \left( \begin{array}{c} \binom{k}{k-1} A_1 \odot \text{Id}_n^{\odot k-1} \\ \binom{k}{k-2} A_2 \odot \text{Id}_n^{\odot k-2} \\ \vdots \\ \binom{k}{0} A_k \end{array} \right) \quad \left. \begin{array}{c} A_{\text{LVE}_{\phi}^{k-1}} \\ \end{array} \right\} , \qquad (23)$$

b) and the principal fundamental matrix for $\text{LVE}_{\phi}^k$ is $\Upsilon_k$ from $\exp_{\odot} Y$ in Notation 1.4.4.

The construction of an infinite matrix $\Upsilon_{\text{LVE}_{\phi}} = \exp_{\odot} Y$ follows in such a way, that the first $d_{n,k} := \binom{n+k-1}{n-1}$ rows and columns are a principal fundamental matrix for $\text{LVE}_{\phi}^k$. The definition is recursive and amenable to symbolic computation.

**Example 1.4.6.** For instance, for $k=5$ we have

$$A_{\text{LVE}_{\phi}^5} = \left( \begin{array}{ccccc} 5A_1 \odot \text{Id}_n^{\odot 4} & & & & \\ 10A_2 \odot \text{Id}_n^{\odot 3} & 4A_1 \odot \text{Id}_n^{\odot 3} & & & \\ 10A_3 \odot \text{Id}_n^{\odot 2} & 6A_2 \odot \text{Id}_n^{\odot 2} & 3A_1 \odot \text{Id}_n^{\odot 2} & & \\ 5A_4 \odot \text{Id}_n & 4A_3 \odot \text{Id}_n & 3A_2 \odot \text{Id}_n & 2A_1 \odot \text{Id}_n & \\ A_5 & A_4 & A_3 & A_2 & A_1 \end{array} \right),$$
## 1.4 Higher-order variational equations

See [34] for proofs and further elaboration on everything summarized below:

**Notation 1.4.1.** For every set of indices $1 \le i_1 \le \dots \le i_r$ such that $\sum_{j=1}^r i_j = k$, $c_{i_1, \dots, i_r}^k$ is defined as the amount of totally ordered partitions of a set of $k$ elements among subsets of sizes $i_1, \dots, i_r$:

$$ c_{i_1, \dots, i_r} = c_{i_1, \dots, i_r}^k = \frac{\binom{k}{i_1, i_2, \dots, i_r}}{n_1! \dots n_m!}, \quad \left\{ \begin{array}{l} (i_1, \dots, i_r) = (k_1, \dots, k_1, k_2, \dots, k_m, \dots, k_m), \\ 1 \le k_1 < k_2 < \dots < k_m. \end{array} \right. \qquad (18) $$

**Proposition 1.4.2.** Let $K = \mathbb{C}(\phi)$ and $A_k$, $Y_k$ be $\partial_k X(\phi)$, $\partial_z^k \psi(t, \phi)$ minus crossed terms; let $A = J_X^\phi$, $Y = J_\phi$ be the derivative jets for $X$ and $\psi$ at $\phi$, with $A_k$, $Y_k$ as blocks. Then, if $c_i^k = \#\{\text{ordered } i_1, \dots, i_r \text{-partitions of } k \text{ elements}\}$,

$$ \dot{Y}_k = \sum_{j=1}^{k} A_j \sum_{|i|=k} c_i^k Y_{i_1} \odot Y_{i_2} \odot \dots \odot Y_{i_r}, \quad k \ge 1. \qquad (\text{VE}_\phi^k) $$

**Lemma 1.4.3.** Let $Y \in \text{Mat}(K)$ equal to zero outside of block row 1, $k \ge 1$:

$$ Y := \left( \begin{array}{cccc} \cdots & Y_3 & Y_2 & Y_1 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{array} \right), \quad Y_i \in \text{Mat}_n^{1,i}. \qquad (19) $$

Let $Z_{r,s}$, $s, r \ge 1$, be the corresponding block in $\exp_\odot Y$. Then,

a) Row block $r$ in $\exp_\odot Y$ is recursively obtained in terms of row blocks 1 and $r-1$:

$$ Z_{r,s} = \frac{1}{r} \sum_{j=1}^{s-r+1} \binom{s}{j} Y_j \odot Z_{r-1,s-j}. \qquad (20) $$

In particular, $Z_{r,r} = Y_1^{\odot r}$ and $Z_{r,s} = 0_{d_{n,r},d_{n,s}}$ whenever $r > s$.

b) Using Notation 1.4.1 and (18), for every $s \ge r$

$$ Z_{r,s} = \sum_{i_1+\cdots+i_r=s} c_{i_1, \dots, i_r}^s Y_{i_1} \odot Y_{i_2} \odot \cdots \odot Y_{i_r}. \qquad (21) $$

**Notation 1.4.4.** $K := \mathbb{C}(\phi)$, $A_i := X^{(i)}(\phi)$, $Y_i := \text{lex} \left( \frac{\partial^i}{\partial Z^i} \varphi(t, \phi) \right)$ and, per Lemma 1.4.3,

$$ \Upsilon_1 = Y_1, \quad \Upsilon_k = \left( \begin{array}{c} Z_{k,k} \\ Z_{k-1,k} \\ \vdots \\ Z_{1,k} \end{array} \right) \quad \left( \begin{array}{c} \Upsilon_{k-1} \end{array} \right), \quad k \ge 2, \qquad (22) $$

formed by the first $k$ block rows and columns in $\Upsilon = \exp_\odot Y$. Define $A, Y \in \text{Mat}(K)$ as in Lemma 1.4.3 with the above $A_i, Y_i$ as blocks. Denote the canonical basis on $K^n$ (meaning the set of columns of $\text{Id}_n$) by $\{e_1, \dots, e_n\}$.

**Proposition 1.4.5 (Explicit version of LVE}_\phi^k).** Still following Notation 1.4.4, the infinite system

$$ \dot{X} = A_{\text{LVE}_\phi} X, \qquad A_{\text{LVE}_\phi} := A \odot \exp_\odot \text{Id}_n, \qquad (\text{LVE}_\phi) $$

has $\Upsilon := \exp_\odot Y$ as a solution matrix. Hence, for every $k \ge 1$,
and the principal fundamental matrix $\Upsilon_5$ is

$$ \left( \begin{array}{cccccc} Y_1^{\odot 5} & & & & & \\ 10Y_1^{\odot 3} \odot Y_2 & Y_1^{\odot 4} & & & & \\ 10Y_1^{\odot 2} \odot Y_3 + 15Y_1 \odot Y_2^{\odot 2} & 6Y_1^{\odot 2} \odot Y_2 & Y_1^{\odot 3} & & & \\ 10Y_2 \odot Y_3 + 5Y_1 \odot Y_4 & 4Y_1 \odot Y_3 + 3Y_2 \odot Y_2 & 3Y_1 \odot Y_2 & Y_1^{\odot 2} & & \\ & Y_5 & Y_4 & Y_3 & Y_2 & Y_1 \end{array} \right), \quad (24) $$

hence $(\mathrm{VE}_\phi^k)$ for $k=5$ is the lowest row in $A_{\mathrm{LVE}_\phi^5}$ times the leftmost column in $\Upsilon_5$.

## 1.5 First integrals and higher-order variational equations

Let $F: U \subseteq \mathbb{C}^n \to \mathbb{C}^n$ be a holomorphic function and $\phi: I \subset \mathbb{C} \to U$. Firstly, the flow $\varphi(t, z)$ of $X$ admits, at least formally, Taylor expansion (2) along $\phi$ which is expressible as

$$ \varphi(t, \phi + \xi) = \phi + Y_1 \xi + \frac{1}{2} Y_2 \xi^{\odot 2} + \dots = \phi + J_{\phi} \exp_{\odot} \xi, \quad (25) $$

where $J_\phi$ is the jet for flow $\varphi(t, \cdot)$ along $\phi$, displayed as $Y$ in (19) and defined in Notation 1.4.4 – that is, the matrix whose $\odot$-exponential $\Upsilon$ is a solution matrix for $(\mathrm{LVE}_\phi)$. Secondly, the Taylor series of $F$ along $\phi$ can be written, cfr. [3, Lemma 2] and Notation 1.2.1,

$$ F(\mathbf{y} + \phi) = F(\phi) + \sum_{m=1}^{\infty} \frac{1}{m!} \langle F^{(m)}(\phi), \mathrm{Sym}^m \mathbf{y} \rangle. \quad (26) $$

Lemma 1.3.1 and (17) trivially implies (26) can be expressed as $F(\mathbf{y} + \phi) = M_F^\phi \exp_\odot \mathbf{y}$, where

$$ M_F^\phi = J_F^\phi + F^{(0)}(\phi) := \left( \begin{array}{ccc|c} \cdots & 0 & 0 & 0 \\ \cdots & F^{(2)}(\phi) & F^{(1)}(\phi) & F^{(0)}(\phi) \\ \cdots & 0 & 0 & 0 \end{array} \right) \in \mathrm{Mat}^{1,n}(K), $$

i.e. $J_F^\phi$ is the jet or horizontal strip of lex-sifted partial derivatives of $F$ at $\phi$.

Let

$$ \dot{X} = A_{\mathrm{LVE}_{\phi}^{*}} X, \qquad A_{\mathrm{LVE}_{\phi}^{*}} := -(A \odot \exp_{\odot} \mathrm{Id}_n)^{T}, \qquad (\mathrm{LVE}_{\phi}^{*}) $$

be the **adjoint** or **dual** variational system of (1) along $\phi$.

It is immediate, upon derivation of equation $\Upsilon_k \Upsilon_k^{-1} = \mathrm{Id}_{D_{n,k}}$, that

**Lemma 1.5.1.** $(\Upsilon_k^{-1})^T$ is a principal fundamental matrix of $(\mathrm{LVE}_\phi^k)^*$, $k \ge 1$, hence $\lim_k (\Upsilon_k^{-1})^T$, is a solution to $(\mathrm{LVE}_\phi^*)$.

The following was proven in [30] and recounted in [3, Lemma 7], and was may now be expressed in a simple, compact fashion:

**Lemma 1.5.2.** Let $F$ and $\phi$ be a holomorphic first integral and a non-constant solution of (1) respectively. Let $V := J_F^T$ be the transposed jet of $F$ along $\phi$. Then, $V$ is a solution of $(\mathrm{LVE}_\phi^*)$.

Hence the issue, already considered in [3], is whether some converse holds and a solution of the dual system $(\mathrm{LVE}_\phi^k)^*$ is or can be the set of first $k$ terms of the formal (or perhaps even convergent) series form of a first integral. This was addressed in the aforementioned reference and was put in compact, explicit linearized form in [34], and will be recounted below.

Define

$$ \hat{A} := \left( \begin{array}{ccc|c} \cdots & 0 & 0 & 0 \\ \cdots & A_2 & A_1 & A_0 \\ \cdots & 0 & 0 & 0 \end{array} \right), \quad A_i := X^{(i)}(\phi) \in \mathrm{Mat}_n^{1,i}(K), \quad A_0 := X(\phi) = \dot{\phi}, \quad (27) $$
a) the lower-triangular recursive $D_{n,k} \times D_{n,k}$ form for LVE$_\phi^k$ is $\dot{Y} = A_{LVE_\phi^k}Y$, its system matrix being obtained from the first $k$ row and column blocks of $A_{LVE_\phi}$:

$$A_{\text{LVE}_{\phi}^k} = \begin{pmatrix} \binom{k}{k-1} A_1 \odot \text{Id}_n^{\odot k-1} & & \\ \binom{k}{k-2} A_2 \odot \text{Id}_n^{\odot k-2} & \boxed{A_{\text{LVE}_{\phi}^{k-1}}} & \\ & \vdots & \\ \binom{k}{0} A_k & & \end{pmatrix}, \quad (23)$$

b) and the principal fundamental matrix for LVE$_\phi^k$ is $\Upsilon_k$ from $\exp_\odot Y$ in Notation 1.4.4.

The construction of an infinite matrix $\Upsilon_{LVE_\phi} = \exp_\odot Y$ follows in such a way, that the first $d_{n,k} := \binom{n+k-1}{n-1}$ rows and columns are a principal fundamental matrix for LVE$_\phi^k$. The definition is recursive and amenable to symbolic computation.

**Example 1.4.6.** For instance, for $k=5$ we have

$$A_{\text{LVE}_{\phi}^{5}} = \begin{pmatrix} 5A_{1} \odot \text{Id}_{n}^{\odot 4} & & & & \\ 10A_{2} \odot \text{Id}_{n}^{\odot 3} & 4A_{1} \odot \text{Id}_{n}^{\odot 3} & & & \\ 10A_{3} \odot \text{Id}_{n}^{\odot 2} & 6A_{2} \odot \text{Id}_{n}^{\odot 2} & 3A_{1} \odot \text{Id}_{n}^{\odot 2} & & \\ 5A_{4} \odot \text{Id}_{n} & 4A_{3} \odot \text{Id}_{n} & 3A_{2} \odot \text{Id}_{n} & 2A_{1} \odot \text{Id}_{n} & \\ A_{5} & A_{4} & A_{3} & A_{2} & A_{1} \end{pmatrix},$$

and the principal fundamental matrix $\Upsilon_5$ is

$$\begin{pmatrix} Y_1^{\odot 5} & & & & \\ 10Y_1^{\odot 3} \odot Y_2 & Y_1^{\odot 4} & & & \\ 10Y_1^{\odot 2} \odot Y_3 + 15Y_1 \odot Y_2^{\odot 2} & 6Y_1^{\odot 2} \odot Y_2 & Y_1^{\odot 3} & & \\ 10Y_2 \odot Y_3 + 5Y_1 \odot Y_4 & 4Y_1 \odot Y_3 + 3Y_2 \odot Y_2 & 3Y_1 \odot Y_2 & Y_1^{\odot 2} & \\ Y_5 & Y_4 & Y_3 & Y_2 & Y_1 \end{pmatrix}, \quad (24)$$

hence $(\text{VE}_\phi^k)$ for $k=5$ is the lowest row in $A_{\text{LVE}_\phi^5}$ times the leftmost column in $\Upsilon_5$.

## 1.5 First integrals and higher-order variational equations

Let $F: U \subseteq \mathbb{C}^n \to \mathbb{C}^n$ be a holomorphic function and $\phi: I \subset \mathbb{C} \to U$. Firstly, the flow $\varphi(t, z)$ of $X$ admits, at least formally, Taylor expansion (2) along $\phi$ which is expressible as

$$\varphi(t, \phi + \boldsymbol{\xi}) = \phi + Y_1 \boldsymbol{\xi} + \frac{1}{2} Y_2 \boldsymbol{\xi}^{\odot 2} + \dots = \phi + J_{\phi} \exp_{\odot} \boldsymbol{\xi}, \quad (25)$$

where $J_\phi$ is the jet for flow $\varphi(t, \cdot)$ along $\phi$, displayed as $Y$ in (19) and defined in Notation 1.4.4 – that is, the matrix whose $\odot$-exponential $\Upsilon$ is a solution matrix for $(\text{LVE}_\phi)$. Secondly, the Taylor series of $F$ along $\phi$ can be written, cfr. [3, Lemma 2] and Notation 1.2.1,

$$F(\mathbf{y} + \boldsymbol{\phi}) = F(\boldsymbol{\phi}) + \sum_{m=1}^{\infty} \frac{1}{m!} \langle F^{(m)}(\boldsymbol{\phi}), \text{Sym}^{m}\mathbf{y} \rangle. \quad (26)$$

Lemma 1.3.1 and (17) trivially implies (26) can be expressed as $F(\mathbf{y} + \boldsymbol{\phi}) = M_F^\phi \exp_\odot \mathbf{y}$, where

$$M_F^\phi = J_F^\phi + F^{(0)}(\phi) := \begin{pmatrix} \cdots & 0 & 0 & 0 \\ \cdots & F^{(2)}(\phi) & F^{(1)}(\phi) & F^{(0)}(\phi) \\ \cdots & 0 & 0 & 0 \end{pmatrix} \in \text{Mat}^{1,n}(K),$$
and let

$$ \hat{A}_{\text{LVE}_{\phi}} := \hat{A} \odot \exp_{\odot} \text{Id}_n = \lim_{k} \hat{A}_{\text{LVE}_{\phi}^k} = \lim_{k} \left( \begin{array}{c} \begin{pmatrix} k \\ k \end{pmatrix} X^{(0)}(\phi) \odot \text{Id}_{n}^{\odot k} \\ \begin{pmatrix} k \\ k-1 \end{pmatrix} X^{(1)}(\phi) \odot \text{Id}_{n}^{\odot k-1} \\ \begin{pmatrix} k \\ k-2 \end{pmatrix} X^{(2)}(\phi) \odot \text{Id}_{n}^{\odot k-2} \\ \vdots \\ \begin{pmatrix} k \\ 0 \end{pmatrix} X^{(k)}(\phi) \odot \text{Id}_{n}^{\odot 0} \end{array} \begin{array}{c} \hat{A}_{\text{LVE}_{\phi}^{k-1}} \\ \vdots \\ \hat{A}_{\text{LVE}_{\phi}^{k-1}} \end{array} \right), $$

The following gave a compact form to [3, Th. 12] in terms of $\odot$ and infinite matrices:

**Proposition 1.5.3** ([34, Prop. 6]). Let $F, \phi, V = (\cdots | V_3 | V_2 | V_1)$ as in Lemma 1.5.2. Then $\hat{A}_{\text{LVE}_{\phi}^T} V = 0$. More specifically, $\hat{A}_{\text{LVE}_{\phi}^{k-1}} (V_k | \cdots | V_2 | V_1) = 0$ for every $k \ge 1$, i.e.

$$ \sum_{j=0}^{k-1} \binom{k-1}{j} (A_j \odot \text{Id}_n^{k-1-j})^T V_{k-j} = 0, \quad \text{for every } k \ge 1. \quad \square \qquad (28) $$

Hence, blocks in $V_1, (V_2, V_1)^T, (V_3, V_2, V_1)^T, \dots$ having all entries in the base field $K$ and satisfying both equations in Proposition 1.5.3 and 1.5.2 are jet blocks $F^{(1)}, F^{(2)}, \dots$ of a formal series that will be a first integral if convergent. These blocks belonging to the intersection of $\ker \hat{A}_{\text{LVE}_{\phi}^{k-1}}$ and the solution subspace $\text{Sol}_K (\text{LVE}_{\phi}^k)^*$ were called *admissible* solutions of the order-$k$ adjoint system in [3]. We will call the sum constructed from an admissible solution a **formal first integral**.

For a short introduction on how changes of variables are reflected on this scheme, as well as the explicit form for the *monodromy matrix* [39] of $\text{LVE}_{\phi}^k$ along any path $\gamma$, see of [34, §4 & 5].

# 2 Jet filtering methods

## 2.1 Preliminaries

Let $f$ be a first integral of (1). The following are immediate:

**Lemma 2.1.1.** (i) For every $i \ge 1$,

$$ (fg)^{(i)} = \sum_{j=0}^{i} \binom{i}{j} f^{(j)} \odot g^{(i-j)}, \qquad (29) $$

with the understanding that symmetric product by a scalar is identical to the usual product, hence

$$ \begin{aligned} J\{fg\} &= U \odot (U^{\odot 2})^T [J\{f\} \odot J\{g\}] = U \odot (U^{\odot 2})^T [\exp_{\odot} J\{f\} \odot \exp_{\odot} J\{g\}] \\ &= U \odot [U^T \exp_{\odot} J\{f\}] \odot [U^T \exp_{\odot} J\{g\}] \end{aligned} $$

$$ \text{where } U = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} $$

(ii) For every $k \ge 1$,

$$ (f^k)^{(i)} = k! \sum_{j=1}^{k} \frac{f^{j-1}}{(j-1)!} \sum_{i_1 + \cdots + i_{k-j+1} = i} c_{i_1, \dots, i_{k-j+1}}^i \bigodot_{u=1}^{k-j+1} f^{(iu)}, \qquad (30) $$
i.e. $J_F^\phi$ is the jet or horizontal strip of lex-sifted partial derivatives of $F$ at $\phi$.

Let

$$ \dot{X} = A_{\text{LVE}_{\phi}^{*}} X, \qquad A_{\text{LVE}_{\phi}^{*}} := -(A \odot \exp_{\odot} \text{Id}_n)^T, \qquad (\text{LVE}_{\phi}^{*}) $$

be the **adjoint** or **dual** variational system of (1) along $\phi$.

It is immediate, upon derivation of equation $\Upsilon_k \Upsilon_k^{-1} = \text{Id}_{D_{n,k}}$, that

**Lemma 1.5.1.** $(\Upsilon_k^{-1})^T$ is a principal fundamental matrix of $(\text{LVE}_\phi^k)^*$, $k \ge 1$, hence $\lim_k (\Upsilon_k^{-1})^T$, is a solution to $(\text{LVE}_\phi^*)$.

The following was proven in [30] and recounted in [3, Lemma 7], and was may now be expressed in a simple, compact fashion:

**Lemma 1.5.2.** Let $F$ and $\phi$ be a holomorphic first integral and a non-constant solution of (1) respectively. Let $V := J_F^T$ be the transposed jet of $F$ along $\phi$. Then, $V$ is a solution of $(\text{LVE}_\phi^*)$.

Hence the issue, already considered in [3], is whether some converse holds and a solution of the dual system $(\text{LVE}_\phi^k)^*$ is or can be the set of first $k$ terms of the formal (or perhaps even convergent) series form of a first integral. This was addressed in the aforementioned reference and was put in compact, explicit linearized form in [34], and will be recounted below.

Define

$$ \hat{A} := \begin{pmatrix} \cdots & 0 & 0 & 0 \\ \cdots & A_2 & A_1 & A_0 \\ \cdots & 0 & 0 & 0 \end{pmatrix}, \quad A_i := X^{(i)}(\phi) \in \text{Mat}_n^{1,i}(K), \quad A_0 := X(\phi) = \dot{\phi}, \quad (27) $$

and let

$$ \hat{A}_{\text{LVE}_{\phi}} := \hat{A} \odot \exp_{\odot} \text{Id}_n = \lim_{k} \hat{A}_{\text{LVE}_{\phi}^k} = \lim_{k} \left( \begin{array}{c} \begin{pmatrix} k \\ k-1 \end{pmatrix} X^{(0)}(\phi) \odot \text{Id}_n^{\odot k} \\ \begin{pmatrix} k \\ k-2 \end{pmatrix} X^{(1)}(\phi) \odot \text{Id}_n^{\odot k-1} \\ \begin{pmatrix} k \\ k-2 \end{pmatrix} X^{(2)}(\phi) \odot \text{Id}_n^{\odot k-2} \\ \vdots \\ \begin{pmatrix} k \\ 0 \end{pmatrix} X^{(k)}(\phi) \odot \text{Id}_n^{\odot 0} \end{array} \quad \hat{A}_{\text{LVE}_{\phi}^{k-1}} \right), \quad (28) $$

The following gave a compact form to [3, Th. 12] in terms of $\odot$ and infinite matrices:

**Proposition 1.5.3** ([34, Prop. 6]). Let $F$, $\phi$, $V = (\cdots | V_3 | V_2 | V_1)$ as in Lemma 1.5.2. Then $\hat{A}_{\text{LVE}_\phi}^T V = 0$. More specifically, $\hat{A}_{\text{LVE}_\phi^{k-1}} (V_k | \cdots | V_2 | V_1) = 0$ for every $k \ge 1$, i.e.

$$ \sum_{j=0}^{k-1} \binom{k-1}{j} (A_j \odot \text{Id}_n^{k-1-j})^T V_{k-j} = 0, \quad \text{for every } k \ge 1. \quad \square \qquad (29) $$

Hence, blocks in $V_1, (V_2, V_1)^T, (V_3, V_2, V_1)^T, \dots$ having all entries in the base field $K$ and satisfying both equations in Proposition 1.5.3 and 1.5.2 are jet blocks $F^{(1)}, F^{(2)}, \dots$ of a formal series that will be a first integral if convergent. In other words:

**Definition 1.5.4** ([3]). In the above hypotheses,

(i) The vectors $(V_k, \dots, V_1)^T$ belonging to the intersection of $\ker \hat{A}_{\text{LVE}_\phi^{k-1}}^T$ and the solution subspace $\text{Sol}_K (\text{LVE}_\phi^k)^*$ are called **admissible** solutions of the order-$k$ adjoint system.

(ii) We call the sum $\sum_{k \ge 1} \frac{1}{k!} V_k z^{\odot k}$ a **formal first integral** along $\phi$ if every finite truncation thereof $\sum_{k=1}^m \frac{1}{k!} V_k z^{\odot k}$ arises from an admissible solution $(V_m, \dots, V_1)$.

For a short introduction on how changes of variables are reflected on this scheme, as well as the explicit form for the *monodromy matrix* [40] of $\text{LVE}_\phi^k$ along any path $\gamma$, see of [34, §4 & 5].
hence defining

$$E_1 = \left( \begin{array}{c|c} e_1^T & 0 \\ \hline 0 & 0 \end{array} \right), \quad U = \left( \begin{array}{c|c} 0 & 1 \\ \hline 0 & 0 \end{array} \right),$$

we have

$$J\{f^k\} = k! \left\{ U \odot (U^{\odot k})^T [\exp_{\odot} J\{f\}] \right\} = k! \left\{ U \odot (U^{\odot k})^T [\exp_{\odot} J\{f\}] \right\} \quad (31)$$

$$= U \odot [U^T \exp_{\odot} J\{f\}]^{\odot k} = U \odot [U^T J\{f\}]^{\odot k} \quad (32)$$

(iii) For every $k \ge 1$,

$$f^k (1/f^k)^{(i)} = k \sum_{j=0}^{i} \frac{(k+j-1)!}{k!} f^{-j} (-1)^j \sum_{i_1+\cdots+i_j=i} c_{i_1, \ldots, i_j}^{i} \bigotimes_{u=1}^{j} f^{(i_u)}, \quad (33)$$

hence

$$J\{1/f\} = \frac{1}{f} U \odot \left\{ \begin{pmatrix} \vdots & \vdots \\ 0 & 3! (-1/f)^3 \\ 0 & 2! (-1/f)^2 \\ 0 & 1! (-1/f) \\ \hline 0 & 1 \end{pmatrix} \right\}^T \exp_{\odot} (J\{f\} - fU) \quad (34)$$

$$= U \odot \left( \frac{1}{f} U^T J\{f\} \right)^{\odot -1} = U \odot \left( \sum_{i=0}^{\infty} (-1)^i \left[ \frac{1}{f} U^T J\{f\} - 1 \right]^{\odot i} \right) \quad (35)$$

## 2.2 Progressive filtering

Let us describe a method to start with a jet of valuation one and then proceed to compute each jet term of that particular function immediately in terms of the previous ones.

**Lemma 2.2.1.** *Define*

$$F_1 := F_{i,1} := \mathrm{Id}_n - \frac{1}{X_i^0} (X^0 \odot e_i^T). \quad (36)$$

The $n-1$ non-zero rows of $F_1 Y_1^{-1}$ are linearly independent admissible solutions of degree one.

*Proof.* Let $F_{i,1} := \mathrm{Id}_n - \frac{1}{X_i(\phi(0))} (X(\phi(0)) \odot e_i^T)$. We have $\hat{A}_{\mathrm{LVE}_\phi^0} = A_0^T = X(\phi)^T$. And $(\Upsilon_1^{-1})^T = (Y_1^{-1})^T$, the fundamental matrix for $\mathrm{VE}_1^*$. Thus checking that (28) holds equates to checking whether $A_0^T (Y_1^{-1}) F_{i,1}^T = 0_n^T$, i.e. that $F_{i,1} Y_1^{-1} A_0$ equals a column of zeros. This is easy to prove by noticing that the basic chain rule and (1) entails $\dot{A}_0 = A_1 A_0$ on one hand, and on the other Lemma 1.5.1 for $k=1$ implies $\frac{d}{dt} Y_1^{-1} = -Y_1^{-1} A_1$. Both of these facts put together imply

$$\frac{d}{dt} (Y_1^{-1} A_0) = \overleftarrow{Y_1^{-1}} X + Y^{-1} \dot{X} = -Y^{-1} A_1 X + Y^{-1} A_1 X = 0$$

and thus $\frac{d}{dt} (F_{i,1} Y_1^{-1} A_0) = F_{i,1} \frac{d}{dt} (Y_1^{-1} A_0) = 0$ as well, which means $F_{i,1} Y_1^{-1} A_0$ is a constant vector. Setting $t = t_0$ and using the fact that $Y_1$ is a principal fundamental matrix (i.e. $Y_1(t_0) =$
# 2 Jet filtering methods and computation of formal first integrals

## 2.1 Preliminaries

Let $f$ be a first integral of (1). The following are immediate:

**Lemma 2.1.1.** (i) For every $i \ge 1$,

$$ (fg)^{(i)} = \sum_{j=0}^{i} \binom{i}{j} f^{(j)} \odot g^{(i-j)}, \qquad (30) $$

with the understanding that symmetric product by a scalar is identical to the usual product, hence

$$ \begin{aligned} J\{fg\} &= U \odot (U^{\odot 2})^T [J\{f\} \odot J\{g\}] = U \odot (U^{\odot 2})^T [\exp_{\odot} J\{f\} \odot \exp_{\odot} J\{g\}] \\ &= U \odot [U^T \exp_{\odot} J\{f\}] \odot [U^T \exp_{\odot} J\{g\}] \end{aligned} $$

$$ \text{where } U = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} $$

(ii) For every $k \ge 1$,

$$ (f^k)^{(i)} = k! \sum_{j=1}^{k} \frac{f^{j-1}}{(j-1)!} \sum_{i_1 + \cdots + i_{k-j+1} = i} c_{i_1, \ldots, i_{k-j+1}}^{i} \bigcirc_{u=1}^{k-j+1} f^{(i_u)}, \qquad (31) $$

hence defining

$$ E_1 = \left( \begin{array}{c|c} e_1^T & 0 \\ \hline 0 & 0 \end{array} \right), \quad U = \left( \begin{array}{c|c} 0 & 1 \\ \hline 0 & 0 \end{array} \right), $$

we have

$$ J\{f^k\} = k! \left\{ U \odot (U^{\odot k})^T [\exp_{\odot} J\{f\}] \right\} = k! \left\{ U \odot (U^{\odot k})^T [\exp_{\odot} J\{f\}] \right\} \quad (32) $$

$$ = U \odot [U^T \exp_{\odot} J\{f\}]^{\odot k} = U \odot [U^T J\{f\}]^{\odot k} \qquad (33) $$

(iii) For every $k \ge 1$,

$$ f^k (1/f^k)^{(i)} = k \sum_{j=0}^{i} \frac{(k+j-1)!}{k!} f^{-j} (-1)^j \sum_{i_1 + \cdots + i_j = i} c_{i_1, \ldots, i_j}^{i} \bigcirc_{u=1}^{j} f^{(i_u)}, \qquad (34) $$

hence

$$ J\{1/f\} = \frac{1}{f} U \odot \left\{ \begin{pmatrix} \vdots & \vdots \\ 0 & 3! (-1/f)^3 \\ 0 & 2! (-1/f)^2 \\ 0 & 1! (-1/f) \\ \hline 0 & 1 \end{pmatrix}^T \exp_{\odot} (J\{f\} - fU) \right\} \qquad (35) $$

$$ = U \odot \left( \frac{1}{f} U^T J\{f\} \right)^{\odot-1} = U \odot \left( \sum_{i=0}^{\infty} (-1)^i \left[ \frac{1}{f} U^T J\{f\} - 1 \right]^{\odot i} \right) \qquad (36) $$
Id) means that for $t = t_0$, abridging notation $X^0 = (X_1^0, \dots, X_n^0) = X(\phi(t_0))$,

$$ F_{i,1}Y_1(t_0)^{-1}A_0 = F_{i,1}A_0 = X(\phi(t_0)) - \frac{1}{X_i(\phi(t_0))} (X(\phi(t_0)) \odot e_i^T) X(\phi(t_0)) \\ = \begin{pmatrix} X_1^0 \\ X_2^0 \\ \vdots \\ X_i^0 \\ \vdots \\ X_{n-1}^0 \\ X_n^0 \end{pmatrix} - \frac{1}{X_i^0} \begin{pmatrix} \boxed{0_{n \times (i-1)}} & X_1^0 & X_2^0 & \vdots & X_i^0 & \vdots & X_{n-1}^0 & X_n^0 \\ \vdots & \boxed{0_{n \times (n-i)}} & & & \vdots & & & \\ \vdots & & \boxed{0_{n \times (n-i)}} & X_1^0 & X_2^0 & \vdots & X_i^0 & \vdots \\ & & & X_{n-1}^0 & X_n^0 & X_{n-1}^0 & X_n^0 & \vdots \\ \vdots & & & & \vdots & \boxed{0_{n \times (n-i)}} & X_{n-1}^0 & X_n^0 \end{pmatrix} \begin{pmatrix} X_1^0 \\ X_2^0 \\ \vdots \\ X_i^0 \\ \vdots \\ X_{n-1}^0 \\ X_n^0 \end{pmatrix} = \mathbf{0}_n. $$

Still for $k = 1$, we need to check that the $n-1$ non-zero columns of $(F_1 Y^{-1})^T$ are linearly independent. This is a consequence of the fact that

$$ F_1^T = \left( \begin{array}{c|c} \boxed{\begin{array}{c} \text{Id}_{i-1} \\ -\frac{X_1}{X_i} \cdots -\frac{X_{i-1}}{X_i} \end{array}} & \begin{array}{c} 0 \\ \vdots \\ 0 \\ 0 \\ \vdots \\ 0 \end{array} \\ \hline \boxed{\begin{array}{c} \text{Id}_{n-i} \\ -\frac{X_{i+1}}{X_i} \cdots -\frac{X_n}{X_i} \end{array}} & \begin{array}{c} 0 \\ \vdots \\ 0 \\ \vdots \\ 0 \end{array} \end{array} \right) $$

thus in $(Y^{-1})^T F_1^T$ column $i$ equals $\mathbf{0}_n$ and its remaining $n-1$ columns are the output of subtracting multiples of column $i$ in $(Y^{-1})^T$ from its remaining (already independent) $n-1$ columns. $\square$

Let $\phi(t)$ be a particular solution of (1) whose parametric expression can be represented with one differential equation involving one coordinate function $z_i$ (case in point: (51), (57) in §3).

**Theorem 2.2.2.** In the above hypotheses, let $f_1$ be one of the columns of $(Y_1^{-1}) F_1^T$. Define $X_k := (Y_1^{\odot k})^T$. Then each term generated by the following recursion

$$ f_k := -X_k^{-1} \int X_k \left\{ \sum_{j=2}^{k} \binom{k}{j} (A_j \odot \text{Id}_{n}^{\odot k-j})^T f_{k-j+1} \right\} dx, \quad (37) $$

is the degree-$k$ term of the expansion of a formal first integral

$$ f(z) = f_1 z + \frac{1}{2} f_2 z^{\odot 2} + \frac{1}{3!} f_3 z^{\odot 3} + \dots, $$

provided that each new integral in (2.2.2) is computed by taking independent variable $x = z_i(t)$, and the limits or constants of integration at each stage can be taken so as to ensure that

$$ f_k (A_0 \odot \text{Id}_n^{\odot k-1}) + \binom{k-1}{1} f_{k-1} (A_1 \odot \text{Id}_n^{\odot k-2}) + \dots + \binom{k-1}{k-2} f_2 (A_{k-2} \odot \text{Id}_n) + f_1 A_{k-1} = 0. \quad (38) $$

*Proof.* Most of the work has been done already by Proposition 1.5.3, Lemma 2.2.1 and the properties intrinsic to $\odot$. Firstly, $X_k^{-1} = [(Y_1^{-1})^{\odot k}]^T$. Secondly, condition (38) is nothing but
## 2.2 Progressive filtering and computation of formal first integrals

Let us describe a method to start with a jet of valuation one and then proceed to compute each jet term of that particular function immediately in terms of the previous ones.

**Lemma 2.2.1.** *Define*

$$F_1 := F_{i,1} := \mathrm{Id}_n - \frac{1}{X_i^0} (X^0 \odot e_i^T). \qquad (37)$$

The $n-1$ non-zero rows of $F_1 Y_1^{-1}$ are linearly independent admissible solutions of degree one.

*Proof.* Let $F_{i,1} := \mathrm{Id}_n - \frac{1}{X_i(\phi(0))} (X(\phi(0)) \odot e_i^T)$. We have $\hat{A}_{\mathrm{LVE}_\phi^0} = A_0^T = X(\phi)^T$. And $(\Upsilon_1^{-1})^T = (Y_1^{-1})^T$, the fundamental matrix for $\mathrm{VE}_1^*$. Thus checking that (29) holds equates to checking whether $A_0^T (Y_1^{-1}) F_{i,1}^T = \mathbf{0}_n^T$, i.e. that $F_{i,1} Y_1^{-1} A_0$ equals a column of zeros. This is easy to prove by noticing that the basic chain rule and (1) entails $\hat{A}_0 = A_1 A_0$ on one hand, and on the other Lemma 1.5.1 for $k=1$ implies $\frac{d}{dt} Y^{-1} = -Y^{-1} A_1$. Both of these facts put together imply

$$\frac{d}{dt} (Y_1^{-1} A_0) = \overleftarrow{Y_1^{-1}} X + Y^{-1} \overrightarrow{X} = -Y^{-1} A_1 X + Y^{-1} A_1 X = 0$$

and thus $\frac{d}{dt} (F_{i,1} Y_1^{-1} A_0) = F_{i,1} \frac{d}{dt} (Y_1^{-1} A_0) = 0$ as well, which means $F_{i,1} Y_1^{-1} A_0$ is a constant vector. Setting $t=t_0$ and using the fact that $Y_1$ is a principal fundamental matrix (i.e. $Y_1(t_0) = \mathrm{Id}$) means that for $t=t_0$, abridging notation $X^0 = (X_1^0, \dots, X_n^0) = X(\phi(t_0))$,

$$F_{i,1} Y_1 (t_0)^{-1} A_0 = F_{i,1} A_0 = X(\phi(t_0)) - \frac{1}{X_i(\phi(t_0))} (X(\phi(t_0)) \odot e_i^T) X(\phi(t_0)) \\ = \begin{pmatrix} X_1^0 \\ X_2^0 \\ \vdots \\ X_i^0 \\ \vdots \\ X_{n-1}^0 \\ X_n^0 \end{pmatrix} - \frac{1}{X_i^0} \begin{pmatrix} \boxed{0_{n \times (i-1)}} & X_1^0 & X_2^0 & \vdots & X_i^0 & \vdots & X_{n-1}^0 & X_n^0 \\ \boxed{0_{n \times (n-i)}} & \boxed{0_{n \times (n-i)}} & X_1^0 & X_2^0 & \vdots & X_i^0 & X_{n-1}^0 & X_n^0 \end{pmatrix} \begin{pmatrix} X_1^0 \\ X_2^0 \\ \vdots \\ X_i^0 \\ \vdots \\ X_{n-1}^0 \\ X_n^0 \end{pmatrix} = \mathbf{0}_n.$$

Still for $k=1$, we need to check that the $n-1$ non-zero columns of $(F_1 Y^{-1})^T$ are linearly independent. This is a consequence of the fact that

$$F_1^T = \left( \begin{array}{c|c} \boxed{\begin{array}{c} \mathrm{Id}_{i-1} \\ -\frac{X_1}{X_i} \cdots -\frac{X_{i-1}}{X_i} \end{array}} & \begin{array}{c} 0 \\ \vdots \\ 0 \\ 0 \\ 0 \\ \vdots \\ 0 \end{array} \\ \hline \boxed{\begin{array}{c} \mathrm{Id}_{n-i} \\ -\frac{X_{i+1}}{X_i} \cdots -\frac{X_n}{X_i} \end{array}} & \begin{array}{c} 0 \\ \vdots \\ 0 \\ \vdots \\ 0 \\ \vdots \\ 0 \end{array} \end{array} \right)$$

thus in $(Y^{-1})^T F_1^T$ column $i$ equals $\mathbf{0}_n$ and its remaining $n-1$ columns are the output of subtracting multiples of column $i$ in $(Y^{-1})^T$ from its remaining (already independent) $n-1$ columns. $\square$

Let $\phi(t)$ be a particular solution of (1) whose parametric expression can be represented with one differential equation involving one single function $x(t)$, e.g. (54), (60), (67) in §3.
a repetition of (28). In order to check the consequences of defining (37), all it takes is to realize
that the first $d_{n,k} \times d_{n,k}$ equations of $(\mathrm{LVE}_{\phi}^{*})$ read

$$
\overline{[(Y_1^{-1})^{\odot k}]^T} = -(A_1 \odot \mathrm{Id}_n^{\odot k-1})^T [(Y_1^{-1})^{\odot k}]^T,
$$

which means $(Y_1^{\odot k})^T$ is a principal fundamental matrix thereof; the homogeneous part of the equations satisfied by $f_k$ displays the same matrix, thus (37) is just plain variation of constants applied to what is said in Lemma 1.5.2. □

2.3 Filtering with a single infinite matrix product

The above result provides a recursive method to compute formal first integrals of arbitrary au-
tonomous systems, thereby summing up the aims of the current paper. The question arises,
however, as to whether a single infinite filter matrix $\Phi = \exp_{\odot} F$ exists, amenable to recur-
sive computation but expressible in compact form and computed without any further quadra-
tures, which multiplied times the fundamental matrix of the dual systems (whose terms are also
computable recursively, as seen above) yields an already-trimmed matrix comprised of linearly
independent admissible solutions in a single swipe.

Building such an infinite quadrature-free filter matrix beyond $k=1$ has a value that is more
theoretical than practical in nature, since the computational inviability of ever-growing matrices
and the effort of discarding the aforementioned cross-products $g_{j_1}^{(i_1)} \odot \cdots \odot g_{j_m}^{(i_m)}$ get in the way of
approximating jets properly. So the question arises: why consider one such matrix? The answer
resides in the recent trends in categorification of dynamical systems (e.g. [7]) and the potential
to transport the techniques known in reductions of flat meromorphic connections to arbitrary
dynamical systems.

A considerable wealth of evidence shows one such quadrature-free matrix exists, at least for
systems defined by analytical vector fields, but a deeper study would warrant techniques beyond
the purview and goals of this article and will thus be left for future work:

**Conjecture 2.3.1.** Let $\phi(t)$ be a non-trivial solution for (1) and $t_0$ such that $X_i(\phi(t_0)) \neq 0$ for some $i$. Let $(\Upsilon_k^{-1})^T$ be the principal fundamental matrix for $(\mathrm{LVE}_\phi^k)^*$ equalling $\mathrm{Id}$ at $t=t_0$. Use superscript 0 to denote values at $t=t_0$ and define $F_1 := F_{i,1}$ as in (36). Then, for every $k \ge 1$, the constant matrix $\Phi_k^T$ defined by the transpose lower left $D_{n,k}$-block of $\Phi = \exp_\odot F$, where $F = (\cdots | F_3 | F_2 | F_1)$ and

$$
F_k = - \frac{1}{X_i^0} \left[ \sum_{j=0}^{k-2} \binom{k-1}{j} F_{j+1} (A_{k-j-1} \odot \mathrm{Id}^{\odot j}) \right] U_k, \quad (39)
$$

where

$$
U_k = \left[ \sum_{j=0}^{k-1} \binom{k}{j+1} (-1)^j (\mathrm{Id} - F_1)^{\odot j} \odot \mathrm{Id}^{\odot k-1-j} \right] \odot e_i^T = \left[ \bigodot_{j=1}^{k-1} (\mathrm{Id} - \zeta_k^j F_1) \right] \odot e_i^T \quad (40)
$$

and $\zeta_k = \exp(2\pi i/k)$ is the root of unity, is such that the last $n-1$ non-zero columns of
$(\Upsilon_k^{-1})^T (\Phi_k)^T$ are linearly independent admissible solutions $g_1^{(\cdot)}, g_2^{(\cdot)}, \dots, g_{n-1}^{(\cdot)} \in \text{Sol} (\text{LVE}_{\phi}^k)^* \cap$
$\ker \hat{A}_{\text{LVE}_{\phi}^{k-1}}^T$ and jets of functionally independent formal first integrals $g_1, \dots, g_{n-1}$ of valuation
1.

**Remarks 2.3.2.**
**Theorem 2.2.2.** In the above hypotheses, let $f_1$ be one of the columns of $(Y_1^{-1}) F_1^T$. Define $X_k := (Y_1^{\odot k})^T$. Then each term generated by the following recursion

$$f_k := -X_k^{-1} \int X_k \left\{ \sum_{j=2}^{k} \binom{k}{j} (A_j \odot \mathrm{Id}_n^{\odot k-j})^T f_{k-j+1} \right\} dx, \quad (38)$$

is the degree-k term of the expansion of a formal first integral

$$f(z) = f_1 z + \frac{1}{2} f_2 z^{\odot 2} + \frac{1}{3!} f_3 z^{\odot 3} + \dots,$$

provided that the limits or constants of integration at each stage (38) are such that

$$f_k (A_0 \odot \mathrm{Id}_n^{\odot k-1}) + \binom{k-1}{1} f_{k-1} (A_1 \odot \mathrm{Id}_n^{\odot k-2}) + \dots + \binom{k-1}{k-2} f_2 (A_{k-2} \odot \mathrm{Id}_n) + f_1 A_{k-1} = 0. \quad (39)$$

*Proof.* Most of the work has been done already by Proposition 1.5.3, Lemma 2.2.1 and the properties intrinsic to $\odot$. Firstly, $X_k^{-1} = [(Y_1^{-1})^{\odot k}]^T$. Secondly, condition (39) is nothing but a repetition of (29). In order to check the consequences of defining (38), all it takes is to realize that the first $d_{n,k} \times d_{n,k}$ equations of $(\mathrm{LVE}_{\phi}^{\star})$ read

$$\overline{[(Y_1^{-1})^{\odot k}]^T} = -(A_1 \odot \mathrm{Id}_n^{\odot k-1})^T [(Y_1^{-1})^{\odot k}]^T,$$

which means $(Y_1^{\odot k})^T$ is a principal fundamental matrix thereof; the homogeneous part of the equations satisfied by $f_k$ displays the same matrix, thus (38) is just plain variation of constants applied to what is said in Lemma 1.5.2. $\square$

## 2.3 Filtering with a single infinite matrix product

The above result provides a recursive method to compute formal first integrals of arbitrary autonomous systems, thereby summing up the immediate aims of the current paper. The question arises, however, as to whether a single infinite filter matrix $\Phi = \exp_{\odot} F$ exists, amenable to recursive computation but expressible in compact form and computed without any further quadratures, which multiplied times the fundamental matrix of the dual systems (whose terms are also computable recursively, as seen above) yields an already-trimmed matrix comprised of linearly independent admissible solutions in a single swipe.

Building such an infinite quadrature-free filter matrix beyond $k=1$ has a value that is more theoretical than practical in nature at this stage, since the computational inviability of ever-growing matrices and the effort of discarding the cross-products $g_{j_1}^{(i_1)} \odot \dots \odot g_{j_m}^{(i_m)}$ get in the way of approximating jets properly. So the question arises: why consider one such matrix? The answer resides in the recent trends in categorification of dynamical systems (e.g. [7]) and the potential to transport the techniques known in reductions of linear differential systems and flat meromorphic connections to arbitrary dynamical systems.

In this Section we prove such a matrix exists for meromorphic and holomorphic systems.

### 2.3.1 Statement of the main result on infinite matrix filtering

Let $\phi(t)$ be a non-trivial solution for (1) and $t_0$ such that $X_i(\phi(t_0)) \neq 0$ for some $i$. Let $\Upsilon_k$ be the principal fundamental matrix for $(\mathrm{LVE}_{\phi}^k)$ such that $\Upsilon_k(t_0) = \mathrm{Id}$ and $F_1 := F_{i,1}$ as in (37).
1. We already know all the columns of $(\Upsilon_k^{-1})^T (\Phi_k)^T$ belong to $\text{Sol}(\text{LVE}_\phi^k)^*$, because $\frac{d}{dt} (\Upsilon^{-1})^T = A_{\text{LVE}^*} (\Upsilon^{-1})^T$ where $\Upsilon = \exp_{\odot} Y$, and the above finite-order products are nothing but truncations of the right-product of $(\Upsilon^{-1})^T$ by a constant matrix. Thus it only remains to show that these truncations are also annihilated by left-multiplication by $\hat{A}_{\text{LVE}_\phi^k}^T$, and that the non-zero columns in $(Y_1^{-1})^T F_1^T$ are linearly independent; independence of the full non-zero columns in higher orders will be a trivial consequence of the latter fact.

2. In light of Remark 1., Lemma 2.2.1 basically entails that the entire conjecture holds true at level $k=1$ and, as seen in §2.2, that was all we needed to compute independent jets of formal first integrals, regardless whether 2.3.1 loses its conjectural status:

**Proposition 2.3.3.** *The conjecture is true at $t_0$ for every $k \ge 1$.*

*Proof.* All occurrences of $A_i$ from here onward will refer to constant matrices $A_i^0 = X^{(i)}(\phi(t_0))$. Case $k=1$ has already been tackled paragraphs earlier; in general, we need to prove that

$$F_k (A_0 \odot \text{Id}^{\odot k-1}) + \binom{k-1}{1} F_{k-1} (A_1 \odot \text{Id}^{k-2}) + \dots + \binom{k-1}{k-2} F_2 (A_{k-2} \odot \text{Id}) + F_1 (A_{k-1}) = 0, (41)$$

for every $k \ge 1$.

First of all, bilinearity from Prop. 1.2.5 and a simple induction argument imply, for every set of $n \times n$ matrices $M_1, \dots, M_m$,

$$\bigcirc_{i=1}^{m} (\text{Id} + M_i) = \text{Id}_{n}^{\odot m} + \sum_{j=1}^{m} \text{Id}_{n}^{\odot m-1} \odot M_j + \sum_{1 \le j_1 < j_2 \le m} \text{Id}_{n}^{\odot m-2} \odot M_{j_1} \odot M_{j_2} + \dots + M_1 \odot \dots \odot M_m$$

and thus setting $m=k-1$ and $M_i = -\zeta_k^i F_1$ very basic properties of cyclotomic polynomials imply

$$U_k = \left[ \bigcirc_{i=1}^{k-1} (\text{Id} - \zeta_k^i F_i) \right] \odot e_i^T = \left( \sum_{j=0}^{k-1} \text{Id}_n^{\odot j} \odot F_1^{\odot k-1-j} \right) \odot e_i^T \quad (42)$$

We claim that

$$U_k (A_0 \odot \text{Id}_n^{\odot k-1}) = X_i \text{Id}_n^{\odot k-1}, \quad \text{for every } k \ge 1. \quad (43)$$

Indeed, we have, using (12) and the rest of the properties of $\odot$ (all of which can be traced back to Lemma 1.2.4 or the universal property of $\odot$, [34, §2.1]),

$$
\begin{align*}
U_k (A_0 \odot \text{Id}_n^{\odot k-1})
  &= \left[ \left( \sum_{j=0}^{k-1} \text{Id}_n^{\odot j} \odot F_1^{\odot k-1-j} \right) \odot e_i^T \right] (A_0 \odot \text{Id}_n^{\odot k-1}) \\
  &= \frac{k-1}{k} \left( \sum_{j=0}^{k-1} \text{Id}_n^{\odot j} \odot F_1^{\odot k-1-j} \right) (A_0 \odot \text{Id}_n^{\odot k-2}) \odot e_i^T + \frac{X_i}{k} \left( \sum_{j=0}^{k-1} \text{Id}_n^{\odot j} \odot F_1^{\odot k-j} \right) \\
  &= \frac{k-1}{k} \sum_{j=0}^{k-2} \frac{k-1-j}{k-1} F_1^{\odot j} \odot A_0 \odot \text{Id}_n^{\odot k-2-j} \odot e_i^T + \frac{X_i}{k} \left( \sum_{j=0}^{k-1} \text{Id}_n^{\odot j} \odot F_1^{\odot k-j} \right) \\
  &= \frac{1}{k} \sum_{j=0}^{k-2} (k-1-j) F_1^{\odot j} \odot A_0 \odot \text{Id}_n^{\odot k-2-j} \odot e_i^T + \frac{X_i}{k} \left( \sum_{j=0}^{k-1} \text{Id}_n^{\odot j} \odot F_1^{\odot k-j} \right) \\
  &= X_i \text{Id}_n^{\odot k-1}.
\end{align*}
$$
**Theorem 2.3.1.** In the above conditions, assume $X$ is holomorphic on a neighborhood $U$ of $\phi(t_0)$ and $Y_1$ and $Y_1^{-1}$ are holomorphic on a neighborhood $V$ of $t_0$. For every $k \ge 1$, define the constant matrices

$$F_k = -\frac{1}{X_i^0} \left[ \sum_{j=0}^{k-2} \binom{k-1}{j} F_{j+1} (A_{k-j-1}^0 \odot \mathrm{Id}^{\odot j}) \right] U_k, \quad (40)$$

where we use superscript 0 to denote values at $t = t_0$,

$$U_k = \left[ \sum_{j=0}^{k-1} \binom{k}{j+1} (-1)^j (\mathrm{Id} - F_1)^{\odot j} \odot \mathrm{Id}^{\odot k-1-j} \right] \odot e_i^T = \left[ \bigodot_{j=1}^{k-1} (\mathrm{Id} - \zeta_k^j F_1) \right] \odot e_i^T, \quad (41)$$

and $\zeta_k = \exp(\frac{2\pi i}{k})$ is the primitive root of unity. Let $\Phi_k$ be the lower square right $D_{n,k}$-block of

$$\Phi = \exp_{\odot} F, \quad \text{where} \quad F = (\cdots | F_3 | F_2 | F_1).$$

Then, the last $n-1$ non-zero columns of $(\Upsilon_k^{-1})^T (\Phi_k)^T$ are linearly independent admissible solutions $g_1^{(\cdot)}, g_2^{(\cdot)}, \dots, g_{n-1}^{(\cdot)} \in \mathrm{Sol} (\mathrm{LVE}_{\phi}^k)^* \cap \ker \hat{A}_{\mathrm{LVE}_{\phi}^k}^T$, i.e. truncated jets of functionally independent formal first integrals $g_1, \dots, g_{n-1}$ of valuation 1, defined on a neighborhood of $\phi(t_0)$.

**Remarks 2.3.2.**

1. We already know all the columns of $(\Upsilon_k^{-1})^T (\Phi_k)^T$ belong to $\mathrm{Sol} (\mathrm{LVE}_{\phi}^k)^*$, because $\frac{d}{dt} (\Upsilon^{-1})^T = A_{\mathrm{LVE}^*} (\Upsilon^{-1})^T$ where $\Upsilon = \exp_{\odot} Y$, and the above finite-order products are nothing but truncations of the right-product of $(\Upsilon^{-1})^T$ by a constant matrix. We also know that the non-zero columns in $(Y_1^{-1})^T F_1^T$ are linearly independent, from whence follows the independence of the full non-zero columns in higher orders.

Thus all that is needed to establish Theorem 2.3.1 is to prove that these truncations are also annihilated by left-multiplication by $\hat{A}_{\mathrm{LVE}_{\phi}^k}^T$. This will be done in §2.3.2.

2. In light of Remark 1., Lemma 2.2.1 basically entails that the entire Theorem holds true at level $k=1$ and, as seen in §2.2, that was all we needed to compute independent jets of formal first integrals, independently of Theorem 2.3.1.

## 2.3.2 Proof of Theorem 2.3.1

As per Remark 1., let us now prove that matrices defined by (40) and (41) are indeed lower blocks of a filter matrix $\Phi$ reducing $(\Upsilon^{-1})^T$ to the desired form (namely one such its lower right $D_{n,k} \times D_{n,k}$ blocks belong to $\ker \hat{A}_{\mathrm{LVE}_{\phi}^k}^T$). In other words, that

$$F \cdot \Upsilon^{-1} \cdot \hat{A}_{\mathrm{LVE}_{\phi}} = F \cdot \Upsilon^{-1} \cdot (\hat{A} \odot \exp_{\odot} \mathrm{Id}_n) = 0, \quad (42)$$

Our first step is proving the following.

**Proposition 2.3.3.** *Identity (42) is true at $t_0$, i.e. $F \cdot \hat{A}_{\mathrm{LVE}_{\phi}}|_{t=t_0} = 0$.*

*Proof.* Obviously we are using the fact that $\Upsilon|_{t=t_0} = \exp_{\odot} \mathrm{Id}$, i.e. an infinite identity matrix. All occurrences of $A_i$ from here onward in this proof will refer to constant matrices $A_i^0 = X^{(i)}(\phi(t_0))$. Case $k=1$ has already been tackled paragraphs earlier; in general, we need to prove that

$$F_k (A_0 \odot \mathrm{Id}^{\odot k-1}) + \binom{k-1}{1} F_{k-1} (A_1 \odot \mathrm{Id}^{k-2}) + \dots + \binom{k-1}{k-2} F_2 (A_{k-2} \odot \mathrm{Id}) + F_1 (A_{k-1}) = 0, \quad (43)$$
The rest follows from (43):

$$
\begin{align*}
F_k (A_0 \odot \mathrm{Id}_n^{\odot k-1}) &= -\frac{1}{X_i} \left[ \sum_{j=0}^{k-2} \binom{k-1}{j} F_{j+1} (A_{k-1-j} \odot \mathrm{Id}_n^{\odot j}) \right] U_k (A_0 \odot \mathrm{Id}_n^{k-1}) \\
&= -\frac{1}{X_i} \left[ \sum_{j=0}^{k-2} \binom{k-1}{j} F_{j+1} (A_{k-1-j} \odot \mathrm{Id}_n^{\odot j}) \right] X_i \mathrm{Id}_n^{\odot k-1} \\
&= -\sum_{j=0}^{k-2} \binom{k-1}{j} F_{j+1} (A_{k-1-j} \odot \mathrm{Id}_n^{\odot j})
\end{align*}
$$

which proves (41) true.

<figure><img src="image_2.png" /></figure>

The following corresponds to single-row blocks (i.e. one first integral at a time) and is a
valid precursor to future studies on how to maintain the filtered admissible structure even after
performing changes of variables, but this is not immediately necessary to our current aims:

**Conjecture 2.3.4.** Let $g$ be a first integral of (1) and $g_i^{(k)}$ its $k^{\text{th}}$ lexicographic derivative, for every $k \ge 1$. Let $g_k := g^{(k)}(\phi(t_0))$ be its value at the original point and define the matrix (constant, just like $F_k$):

$$
M_k = \begin{pmatrix}
g_1^{\odot k} & & & \\
c_{1,...,1,2} g_1^{\odot k-1} \odot g_2 & \ddots & & \\
\vdots & & g_1^{\odot 2} & \\
g_k & \cdots & g_2 & g_1
\end{pmatrix} \in \mathrm{Mat}_{1,n}^{n,n}
$$

be the first $k \times D_{n,k}$ block of $M := \exp_{\odot} (\cdots g_k | g_{k-1} | \cdots | g_2 | g_1)$. Then

$$
(\Upsilon_k^{-1})^T \Phi_k^T M_k^T = \begin{pmatrix} \frac{1}{k!} (g^k)^{(k)} & \cdots & \frac{1}{2} (g^2)^{(k)} & g^{(k)} \\ & \ddots & & \vdots \\ \frac{1}{2} (g^2)^{(2)} & g^{(2)} & & g^{(1)} \end{pmatrix}
$$

In other words: $(\Upsilon^{-1})^T \Phi^T \Gamma^T$ is an infinite matrix whose columns are the Taylor terms of the respective powers of $g$. $\square$

**Example 2.3.5.** In both cases (the quadrature-free one posed by Conjecture **2.3.1** and the proven, quadrature-derived method in Theorem **2.2.2**) every new Taylor block of a given formal first integral can be obtained via a filter matrix; the only difference is whether the matrix has a simple explicit form or not. Let us show the first stages of any filter matrix (proven or not) does to the fundamental matrix. Let us assume, for instance, $i = 1$. Stage $k = 1$ is not even a conjecture in virtue of Lemma **2.2.1** and the first order filter matrix, namely $F_1$ in (36), transforms the transpose of the fundamental matrix of the dual system $\Upsilon_1^{-1} = Y_1^{-1}$ into a matrix of the form

$$
\left( \begin{array}{c} \mathbf{0} \\ G_1 \end{array} \right) := \left( \begin{array}{c} \begin{array}{c} \mathbf{0} \\ g_1^{(1)} \\ g_2^{(1)} \\ \vdots \\ g_{n-1}^{(1)} \end{array} \\ \hline \end{array} \right) = \left( \begin{array}{ccccc} 0 & 0 & \cdots & 0 & 0 \\ g_{1,1}^{(1)} & g_{1,2}^{(1)} & \cdots & g_{1,n-1}^{(1)} & g_{1,n}^{(1)} \\ g_{2,1}^{(1)} & g_{2,2}^{(1)} & \cdots & g_{2,n-1}^{(1)} & g_{2,n}^{(1)} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ g_{n-1,1}^{(1)} & g_{n-1,2}^{(1)} & \cdots & g_{n-1,n-1}^{(1)} & g_{n-1,n}^{(1)} \end{array} \right); \qquad (44)
$$
for every $k \ge 1$.

First of all, bilinearity from Prop. 1.2.5 and a simple induction argument imply, for every set of $n \times n$ matrices $M_1, \dots, M_m$,

$$ \bigcirc_{i=1}^{m} (\mathrm{Id} + M_i) = \mathrm{Id}_{n}^{\odot m} + \sum_{j=1}^{m} \mathrm{Id}_{n}^{\odot m-1} \odot M_j + \sum_{1 \le j_1 < j_2 \le m} \mathrm{Id}_{n}^{\odot m-2} \odot M_{j_1} \odot M_{j_2} + \dots + M_1 \odot \dots \odot M_m $$

and thus setting $m = k - 1$ and $M_i = -\zeta_k^i F_1$ very basic properties of cyclotomic polynomials imply

$$ U_k = \left[ \bigcirc_{i=1}^{k-1} (\mathrm{Id} - \zeta_k^i F_i) \right] \odot e_i^T = \left( \sum_{j=0}^{k-1} \mathrm{Id}_n^{\odot j} \odot F_1^{\odot k-1-j} \right) \odot e_i^T \quad (44) $$

We claim that

$$ U_k(A_0 \odot \mathrm{Id}_n^{\odot k-1}) = X_i \mathrm{Id}_n^{\odot k-1}, \quad \text{for every } k \ge 1. \quad (45) $$

Indeed, we have, using (12) and the rest of the properties of $\odot$ (all of which can be traced back to Lemma 1.2.4 or the universal property of $\odot$, [34, §2.1]),

$$ \begin{aligned} U_k (A_0 \odot \mathrm{Id}_n^{\odot k-1}) &= \left[ \left( \sum_{j=0}^{k-1} \mathrm{Id}_n^{\odot j} \odot F_1^{\odot k-1-j} \right) \odot e_i^T \right] (A_0 \odot \mathrm{Id}_n^{\odot k-1}) \\ &= \frac{k-1}{k} \left( \sum_{j=0}^{k-1} \mathrm{Id}_n^{\odot j} \odot F_1^{\odot k-1-j} \right) (A_0 \odot \mathrm{Id}_n^{k-2}) \odot e_i^T + \frac{X_i}{k} \left( \sum_{j=0}^{k-1} \mathrm{Id}_n^{\odot j} \odot F_1^{\odot k-j} \right) \\ &= \frac{k-1}{k} \sum_{j=0}^{k-2} \frac{k-1-j}{k-1} F_1^{\odot j} \odot A_0 \odot \mathrm{Id}_n^{\odot k-2-j} \odot e_i^T + \frac{X_i}{k} \left( \sum_{j=0}^{k-1} \mathrm{Id}_n^{\odot j} \odot F_1^{\odot k-j} \right) \\ &= \frac{1}{k} \sum_{j=0}^{k-2} (k-1-j) F_1^{\odot j} \odot A_0 \odot \mathrm{Id}_n^{\odot k-2-j} \odot e_i^T + \frac{X_i}{k} \left( \sum_{j=0}^{k-1} \mathrm{Id}_n^{\odot j} \odot F_1^{\odot k-j} \right) \\ &= X_i \mathrm{Id}_n^{\odot k-1}. \end{aligned} $$

The rest follows from (45):

$$ \begin{aligned} F_k (A_0 \odot \mathrm{Id}_n^{\odot k-1}) &= -\frac{1}{X_i} \left[ \sum_{j=0}^{k-2} \binom{k-1}{j} F_{j+1} (A_{k-1-j} \odot \mathrm{Id}_n^{\odot j}) \right] U_k (A_0 \odot \mathrm{Id}_n^{k-1}) \\ &= -\frac{1}{X_i} \left[ \sum_{j=0}^{k-2} \binom{k-1}{j} F_{j+1} (A_{k-1-j} \odot \mathrm{Id}_n^{\odot j}) \right] X_i \mathrm{Id}_n^{\odot k-1} \\ &= -\sum_{j=0}^{k-2} \binom{k-1}{j} F_{j+1} (A_{k-1-j} \odot \mathrm{Id}_n^{\odot j}) \end{aligned} $$

which proves (43) true.

<figure><img src="image_2.png" /></figure>

Double induction, along with properties of $\odot$, easily proves that

$$ A_0 \odot A_u \odot \mathrm{Id}^{\odot v} = (A_0 \odot \mathrm{Id}^{\odot v+1}) (A_u \odot \mathrm{Id}^{\odot v}), \quad u, v \ge 0, $$

and just as easily establishes the following consequence thereof:
each of these row vectors $g_1^{(1)}, \dots, g_{n-1}^{(1)}$ is the potential gradient of a formal first integral and all $n-1$ of them are linearly independent. In other words, they are the outputs of Lemma 2.2.1. The next step would be to integrate the second variational system (from here onwards only quadratures are needed):

$$Y_2(t) = \Phi_1 \int_0^t \Phi_1^{-1}(\tau) A_2(\tau) \Phi_1^{\odot 2}(\tau) d\tau$$

invert the fundamental matrix for (LVE²),

$$\Upsilon_2 = \begin{pmatrix} Y_1^{\odot 2} & 0 \\ Y_2 & Y_1 \end{pmatrix}, \quad \text{thus} \quad \Psi_2 := (\Phi_2^{-1})^T$$

and multiply it by the second-order filter matrix

$$\Phi_2 = \begin{pmatrix} F_1^{\odot 2} & 0 \\ F_2 & F_1 \end{pmatrix},$$

where $F_2$ is obtained from (39), (40) or from (37), and the result we obtain is

$$\Phi_2(\Psi_2)^T = \Phi_2 \Upsilon_2^{-1} = \begin{pmatrix} G_1^{\odot 2} & 0 \\ G_2 & G_1 \end{pmatrix} = \begin{pmatrix} \mathbf{0}_{n} & \mathbf{0}^{\odot 2} \\ \mathbf{0}_{d_{n,2}} & \mathbf{0}_{n} \\ g_1^{(2)} & g_1^{(1)} \\ \vdots & \vdots \\ g_{n-1}^{(2)} & g_{n-1}^{(1)} \end{pmatrix} \quad (45)$$

Take for instance the first non-zero row of the lower block ($G_2 | G_1$), i.e. the one highlighted in gray. Imagine we perform the computations for $(LVE_\phi^3)$ and multiply $\Upsilon_3^{-1}$ by the filter matrix $\Phi_3$; we obtain a matrix having the one in (45) as a lower right $D_{n,k} \times D_{n,k}$ block:

$$\Phi_3 \Upsilon_3^{-1} = \begin{pmatrix} F_1^{\odot 3} & 0 & 0 \\ 3F_1 \odot F_2 & F_1^{\odot} & 0 \\ F_3 & F_2 & F_1 \end{pmatrix} \Upsilon_3^{-1} = \begin{pmatrix} G_1^{\odot 3} & 0 & 0 \\ 3G_1 \odot G_2 & G_1^{\odot} & 0 \\ G_3 & G_2 & G_1 \end{pmatrix} \quad (46) \\ = \begin{pmatrix} \mathbf{0}_{n} & \mathbf{0}^{\odot 3} \\ \mathbf{0}_{d_{n,1}} & \mathbf{0}_{d_{n,2}} \\ \vdots & \vdots \\ \mathbf{0}_{d_{n,n}} & \mathbf{0}_{d_{n,n}} \end{pmatrix} \odot 2$$

thus we have an increasing row, highlighted in darker gray color $\dots | g_1^{(3)} | g_1^{(2)} | g_1^{(1)}$. These vectors are precisely the terms of a formal Taylor series of a first integral $g_1(\mathbf{z})$ of (1) along solution $\phi$:

$$g_1(\mathbf{z} + \phi) = g_1(\phi) + g_1^{(1)}\mathbf{z} + \frac{1}{2}g_1^{(2)}(\phi)\mathbf{z}^{\odot 2} + \frac{1}{6}g_1^{(3)}(\phi)\mathbf{z}^{\odot 3} + \dots \quad (47)$$
**Lemma 2.3.4.** *Using the redundant notation*

$$A_0 = \left( \begin{array}{c|c} \cdots & 0 \\ \hline \cdots & 0 \\ \hline \cdots & 0 \end{array} \right),$$

and in the hypotheses of (27) and (28),

$$ \overline{\hat{A} \odot \exp_{\odot} \mathrm{Id}} = (A \odot \exp_{\odot} \mathrm{Id}) (A_0 \odot \exp_{\odot} \mathrm{Id}) - A_0 \odot A \odot \exp_{\odot} \mathrm{Id} \\ = (A \odot \exp_{\odot} \mathrm{Id}) (A_0 \odot \exp_{\odot} \mathrm{Id}) - (A_0 \odot \exp_{\odot} \mathrm{Id}) (A \odot \exp_{\odot} \mathrm{Id}), $$

in other words,

$$ \binom{i+j}{i} \overline{A_j \odot \mathrm{Id}^i} = \binom{i+j+1}{i} (A_{j+1} \odot \mathrm{Id}^{\odot i}) (A_0 \odot \mathrm{Id}^{\odot (i+j)}) - \binom{i+j}{i-1} (A_{j+1} \odot A_0 \odot \mathrm{Id}^{\odot (i-1)}), $$

for every $i, j \ge 0$. □

Let us study the successive derivatives of $C = F \cdot \Upsilon^{-1} \cdot \hat{A}_{\mathrm{LVE}_{\phi}}$. We know $\overline{\Upsilon^{-1}} = -\Upsilon^{-1} A_{\mathrm{LVE}_{\phi}}$ from Lemma 1.5.1. Using Lemma 2.3.4,

$$ \dot{\overline{C}} = F [-\Upsilon^{-1} A_{\mathrm{LVE}_{\phi}} \hat{A}_{\mathrm{LVE}_{\phi}} + \Upsilon^{-1} A_{\mathrm{LVE}_{\phi}} (A_0 \odot \exp_{\odot} \mathrm{Id}) - \Upsilon^{-1} A_0 \odot A \odot \exp_{\odot} \mathrm{Id}] $$

$$ = F \Upsilon^{-1} [-A_{\mathrm{LVE}_{\phi}} \hat{A}_{\mathrm{LVE}_{\phi}} + A_{\mathrm{LVE}_{\phi}} (A_0 \odot \exp_{\odot} \mathrm{Id}) - A_0 \odot A \odot \exp_{\odot} \mathrm{Id}] $$

$$ = -F \Upsilon^{-1} [A_{\mathrm{LVE}_{\phi}}^2 + A_0 \odot A \odot \exp_{\odot} \mathrm{Id}] = -(\exp_{\odot} F) \Upsilon^{-1} [A_{\mathrm{LVE}_{\phi}}^2 + A_0 \odot A_{\mathrm{LVE}}] $$

$$ = -F \Upsilon^{-1} [A_{\mathrm{LVE}_{\phi}}^2 + (A_0 \odot \exp_{\odot} \mathrm{Id}) A_{\mathrm{LVE}_{\phi}}] $$

$$ = -F \Upsilon^{-1} [\hat{A}_{\mathrm{LVE}_{\phi}} A_{\mathrm{LVE}_{\phi}}], $$

$$ \ddot{\overline{C}} = -F \Upsilon^{-1} \left[ \overline{\hat{A}_{\mathrm{LVE}_{\phi}} A_{\mathrm{LVE}_{\phi}} + \hat{A}_{\mathrm{LVE}_{\phi}} \dot{A}_{\mathrm{LVE}_{\phi}}} \right] $$

$$ = -F \Upsilon^{-1} \left[ -A_{\mathrm{LVE}_{\phi}} \hat{A}_{\mathrm{LVE}_{\phi}} A_{\mathrm{LVE}_{\phi}} + \overline{\hat{A}_{\mathrm{LVE}_{\phi}} A_{\mathrm{LVE}_{\phi}} + \hat{A}_{\mathrm{LVE}_{\phi}} \dot{A}_{\mathrm{LVE}_{\phi}}} \right] $$

$$ = -F \Upsilon^{-1} \left[ -A_{\mathrm{LVE}_{\phi}} \hat{A}_{\mathrm{LVE}_{\phi}} A_{\mathrm{LVE}_{\phi}} + (A_{\mathrm{LVE}_{\phi}} (A_0 \odot \exp_{\odot} \mathrm{Id}) - (A_0 \odot \exp_{\odot} \mathrm{Id}) A_{\mathrm{LVE}_{\phi}}) A_{\mathrm{LVE}_{\phi}} + \overline{\hat{A}_{\mathrm{LVE}_{\phi}} A_{\mathrm{LVE}_{\phi}} \dot{A}_{\mathrm{LVE}_{\phi}}} \right] $$

$$ = -F \Upsilon^{-1} \left[ -A_{\mathrm{LVE}_{\phi}}^2 A_{\mathrm{LVE}_{\phi}} - (A_0 \odot \exp_{\odot} \mathrm{Id}) A_{\mathrm{LVE}_{\phi}}^2 + \overline{\hat{A}_{\mathrm{LVE}_{\phi}} A_{\mathrm{LVE}_{\phi}} \dot{A}_{\mathrm{LVE}_{\phi}}} \right] $$

$$ = F \Upsilon^{-1} \left[ \hat{A}_{\mathrm{LVE}_{\phi}} A_{\mathrm{LVE}_{\phi}}^2 - \overline{\hat{A}_{\mathrm{LVE}_{\phi}} A_{\mathrm{LVE}_{\phi}} \dot{A}_{\mathrm{LVE}_{\phi}}} \right] $$

This leads to the following statement, which can be easily proved by induction: for every $i \ge 1$, there exists such that

$$ \frac{d^i C}{dt^i} = \frac{d^i}{dt^i} (F \cdot \Upsilon^{-1} \cdot \hat{A}_{\mathrm{LVE}_{\phi}}) = F \cdot \Upsilon^{-1} \cdot \hat{A}_{\mathrm{LVE}_{\phi}} \cdot M_i, \quad \text{for some } M_i \in \mathrm{Mat}(K). \quad (46) $$

Case $i = 0$ is true with $M_0 = \exp_{\odot} \mathrm{Id}$ by definition, cases $i = 1$ and $i = 2$ are true with $M_1 = -A_{\mathrm{LVE}_{\phi}}$ and $M_2 = A_{\mathrm{LVE}_{\phi}}^2 - \overline{A_{\mathrm{LVE}_{\phi}}}$ as seen above and if we assume (46) true for an arbitrary $i$, then the proof for $i+1$ works along the lines of case $i = 2$ by way of Lemma 2.3.4:

$$ \frac{d^{i+1} C}{dt^{i+1}} = \frac{d}{dt} (F \Upsilon^{-1} \hat{A}_{\mathrm{LVE}_{\phi}} M_i) = F \Upsilon^{-1} \left[ -A_{\mathrm{LVE}_{\phi}} \hat{A}_{\mathrm{LVE}_{\phi}} M_i + \overline{\hat{A}_{\mathrm{LVE}_{\phi}} M_i + \hat{A}_{\mathrm{LVE}_{\phi}} \dot{M}_i} \right] $$

$$ = F \Upsilon^{-1} \left[ -A_{\mathrm{LVE}_{\phi}} \hat{A}_{\mathrm{LVE}_{\phi}} M_i + (A_{\mathrm{LVE}_{\phi}} (A_0 \odot \exp_{\odot} \mathrm{Id}) - (A_0 \odot \exp_{\odot} \mathrm{Id}) A_{\mathrm{LVE}_{\phi}}) M_i + \overline{\hat{A}_{\mathrm{LVE}_{\phi}} \dot{M}_i} \right] $$

$$ = F \Upsilon^{-1} \left[ -\hat{A}_{\mathrm{LVE}_{\phi}} A_{\mathrm{LVE}_{\phi}} M_i + \overline{\hat{A}_{\mathrm{LVE}_{\phi}} \dot{M}_i} \right] = F \Upsilon^{-1} \hat{A}_{\mathrm{LVE}_{\phi}} \left[ -A_{\mathrm{LVE}_{\phi}} M_i + \dot{M}_i \right], $$
where $g_1^{(i)} = g_1^{(i)}(\phi)$, i.e. the implementation of (3) in **1.2.1 (2.)** to $z = \phi$. The same applies, *mutatis mutandis*, to the rest of rows of infinite block $...G_3 | G_2 | G_1$, i.e. replacing $g_1^{(i)}$ with $g_j^{(i)}$ for some other $j = 2, ..., n-1$. Needless to say the upper terms ... in (46) and higher-order counterparts contain all the cross-products arising from derivatives of products of the first integrals, e.g. $g_1^{(1)} \odot g_3^{(2)}$ would be a row in $G_1 \odot G_2$ (it is easy to check which one by way of lexicographic ordering and we invite the reader to do this themselves, using the examples in §3).

# 3 Examples

## 3.1 Dixon's system

The following system,

$$ \dot{u} = \frac{uv}{u^2 + v^2} - \alpha u, \quad \dot{v} = \frac{v^2}{u^2 + v^2} - \beta v + \beta - 1, \qquad (48) $$

where $\alpha, \beta > 0$, arose in [8] as a decoupled two-dimensional fragment of a three-dimensional dynamical model of the magnetic field of neutron stars. Several references (e.g. [2, 11, 37]) discussed the ostensible chaotic behavior of the dynamics of (48). [33], on the other hand, purported the non-chaotic character of this system based on a version of the Poincaré-Bendixson Theorem that precludes the need for compact sets comprising finitely many sets of homoclinic/heteroclinic connections. Given the existing, seldom-discussed gap between the concepts *non-chaotic* and *integrable* (regardless of which ad-hoc definition of chaos is used), we can afford to eschew the polemic altogether and focus on finding a conserved quantity for (48) using higher variational equations.

We will focus on the case $\alpha = \beta > 0$. The case $\alpha \neq \beta$, as well as the original three-dimensional model ([11, eqs (3)-(5)]), will be tackled in future work.

$\alpha = 1$ is immediate to solve and needs no further consideration. Let us focus on $\alpha \neq 1$. Previous attempts at simplifying higher-order LVE$^k$ suggest change of variable

$$ (u, v) = \left( \frac{\sin y}{x}, \frac{\cos y}{x} \right) \qquad (49) $$

which transforms (48) for $\alpha = \beta$ into

$$ \dot{x} = \alpha x (1 - x \cos y), \quad \dot{y} = -(\alpha - 1) x \sin y. \qquad (50) $$

We will use particular solution

$$ \phi(t) = (x(t), 0) \quad \text{such that} \quad \dot{x} = -\alpha(-1+x)x. \qquad (51) $$

Along $\Gamma$, the higher variational complex and its dual (LVE$^*_\phi$) feature the following blocks for $k \ge 0$,

$$ A_k = \left\{ \begin{array}{ll} X(x, y) = \begin{pmatrix} \alpha x (1 - x \cos y) \\ -(\alpha - 1) x \sin y \end{pmatrix}, & k = 0, \\ \begin{pmatrix} -\alpha (2x - 1) & 0 \\ 0 & -(\alpha - 1) x \end{pmatrix}, & k = 1, \\ \left( \begin{array}{cc} (-1)^{k/2} 2\alpha & 0 \\ 0 & (-1)^{k/2} (-1 + \alpha) \end{array} \right), & k \text{ even} \\ \left( \begin{array}{cc} (-1)^{\frac{k+1}{2}} 2\alpha x & 0 \\ 0 & (-1)^{\frac{k+1}{2}} (-1 + \alpha) x \end{array} \right), & k \text{ odd}, \end{array} \right. \qquad (52) $$
which explicitly defines the sought-after matrices recursively as

$$M_{i+1} = -A_{\text{LVE}_{\phi}} M_i + \dot{M}_i, \quad i \ge 0.$$

Thus (46) is true for every $i \ge 1$ and all we need to do to finish the proof of Theorem 2.3.1 is ascertain that all derivatives of $C$ at $t = t_0$ are identically equal to the zero infinite matrix. But this is true on account of the fact that if we take (46) at $t = t_0$, we obtain

$$\left. \frac{d^{i+1}C}{dt^{i+1}} \right|_{t=t_0} = F \cdot \hat{A}_{\text{LVE}_{\phi}} \Big|_{t=t_0} \cdot M_i(t_0) = 0 \cdot M_i(t_0) = 0,$$

in virtue of Proposition 2.3.3. Thus the matrix $C = F\Upsilon^{-1}\hat{A}_{\text{LVE}_{\phi}}$ of functions holomorphic on the open neighborhood $I := \phi^{-1}(U) \cap V$ of $t_0$ has all of its derivatives with respect to $t$ equal to zero in $t_0$. This means $C \equiv 0$ in a neighborhood $J$ of $t_0$, thus on that neighborhood $(\Upsilon^{-1})^T F^T \in \ker \hat{A}_{\text{LVE}_{\phi}}^T$ which is the last step we needed to prove the columns of $(\Upsilon^{-1})^T F^T$ are admissible solutions defined on a neighborhood of $\{\phi(t) : t \in I \cap J\}$. Remarks 2.3.2 and Lemma 2.2.1 finish the proof. □

### 2.3.3 Further comments

The following corresponds to single-row blocks (i.e. one first integral at a time) and is a valid precursor to future studies on how to maintain the filtered admissible structure even after performing changes of variables, but this is not immediately necessary to our current aims:

**Corollary 2.3.5.** Let $g$ be a first integral of (1) and $g_i^{(k)}$ its $k^{\text{th}}$ lexicographic derivative, for every $k \ge 1$. Let $g_k := g^{(k)}(\phi(t_0))$ be its value at the original point and define the matrix (constant, just like $F_k$):

$$M_k = \begin{pmatrix} g_1^{\odot k} & & & \\ c_{1,...,1,2} g_1^{\odot k-1} \odot g_2 & \ddots & & \\ & \ddots & g_1^{\odot 2} & \\ & & g_2 & g_1 \end{pmatrix} \in \text{Mat}_{1,n}^{n,n}$$

be the first $k \times D_{n,k}$ block of $M := \exp_{\odot} (\cdots g_k | g_{k-1} | \cdots | g_2 | g_1)$. Then

$$(\Upsilon_k^{-1})^T \Phi_k^T M_k^T = \begin{pmatrix} \frac{1}{k!} (g^k)^{(k)} & \cdots & \frac{1}{2} (g^2)^{(k)} & g^{(k)} \\ & \ddots & \vdots & \vdots \\ \frac{1}{2} (g^2)^{(2)} & g^{(2)} & & g^{(1)} \end{pmatrix}$$

In other words: $(\Upsilon^{-1})^T \Phi^T M^T$ is an infinite matrix whose columns are the Taylor terms of the respective powers of $g$. □

**Example 2.3.6.** We have seen two methods of obtaining every new Taylor block of a given formal first integral: the global filter provided by Theorem 2.3.1, and the progressive quadrature method in Theorem 2.2.2. Let us remain focused on the former and assume, for instance, $i = 1$ in Lemma 2.2.1; $F_1$ in (37) transforms the transpose of the fundamental matrix of the dual
for instance a principal fundamental matrix of (VE<sub>φ</sub>) is

$$
\Upsilon_1 = Y_1 = \begin{pmatrix} \frac{(x-1)x}{(x_0-1)x_0} & 0 \\ 0 & (1-x_0)^{\frac{1}{\alpha}-1}(1-x)^{1-\frac{1}{\alpha}} \end{pmatrix}
$$

and the first-order filter matrix is very simple in this case,

$$
F_1 = \mathrm{Id}_2 - \frac{1}{X_1^0} (X^0 \odot e_1^T) = \begin{pmatrix} 0 & 0 \\ -\frac{X_2(x_0, y_0)}{X_1(x_0, y_0)} & 1 \end{pmatrix} = \begin{pmatrix} 0 & 0 \\ 0 & 1 \end{pmatrix},
$$

where $x_0 = x(0)$ and $y_0 = y(0)$; so are the higher-order matrices in $\Phi = \exp_{\odot} F$ obtained from (39), some of which are shown below:

$$
\begin{align*}
F_2 &= \begin{pmatrix} 0 & 0 & 0 \\ 0 & \frac{1-\alpha}{\alpha(x_0-1)} & 0 \end{pmatrix}, &
F_3 &= \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & \frac{(\alpha-1)(2\alpha-1)}{\alpha^2(x_0-1)^2} & 0 & 0 \end{pmatrix}, \\
F_4 &= \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & -\frac{(\alpha-1)(2\alpha-1)(3\alpha-1)}{\alpha^3(x_0-1)^3} & 0 & -\frac{(\alpha-1)(2x_0+1)}{\alpha(x_0-1)^2} \end{pmatrix}, \dots
\end{align*}
$$

Returning to $k=1$, Lemma 2.2.1 ensures that the second column of the filtered matrix

$$
(Y_1^{-1})^T F_1^T = \begin{pmatrix} 0 & 0 \\ 0 & (1-x_0)^{1-\frac{1}{\alpha}}(1-x)^{\frac{1}{\alpha}-1} \end{pmatrix},
$$

is an admissible solution of degree 1. We can normalize this to eliminate appearances of $x_0$:

$$
(\Phi_1^{-1})^T F_1^T P^T := \begin{pmatrix} 0 & 0 \\ 0 & (1-x_0)^{1-\frac{1}{\alpha}}(1-x)^{\frac{1}{\alpha}-1} \end{pmatrix} \begin{pmatrix} 1 & 0 \\ 0 & (1-x_0)^{\frac{1}{\alpha}-1} \end{pmatrix}^T = \begin{pmatrix} 0 & 0 \\ 0 & (1-x)^{\frac{1}{\alpha}-1} \end{pmatrix},
$$

the second of whose columns is still the gradient $f^{(1)}$ of a formal first integral (this is exactly
what appears as a row $g_1^{(1)}$ in (44)).

The higher order trimming (in which we discard redundant symmetric products $f^{(i_1)} \odot \cdots \odot f^{(i_k)}$, i.e. derivatives of powers of $f$) is clear once we write down the dual system of equations and filter out telescopically using either Lemma 2.1.1 and the general form of that each new term $f^{(k)}$ of our formal first integral $f$ is expressible as the following (bear in mind that $x = x(t)$ here): $f^{(k)}(x) = (\frac{d}{dx}f^{(k-1)}|_{g_k(x)})$, where

$$
\sum_{j=0}^{\frac{k-1}{2}} \left( (-1)^{j+1} (\alpha - 1) \binom{k}{-1-2j+k} g_{k-2j}(x) + (-1)^{j+1} \alpha \binom{k}{-2j+k} x g'_{k-2j}(x) \right) + \alpha g'_k(x) = 0,
$$

which means $g_k(x) = 0$ for every even $k \ge 2$ and e.g.

$$
f_1 = \left( 0 \left(1-x\right)^{\frac{1}{\alpha}-1} \right),
$$

$$
f_2 = \left( 0 \frac{(\alpha-1)(1-x)^{\frac{1}{\alpha}-2}}{\alpha} 0 \right),
$$

$$
f_3 = \left( 0 \left( \frac{1}{\alpha} - 2 \right) \left( \frac{1}{\alpha} - 1 \right) (1-x)^{\frac{1}{\alpha}-3} \quad 0 \quad \frac{(1-x)^{\frac{1}{\alpha}-2} (2\alpha + (\alpha-2)x-1)}{\alpha-2} \right),
$$

$$
f_4 = \left( \frac{d}{dx} f_3 \mid 0 \right),
$$

$$
f_5 = \left( \frac{d}{dx} f_4 \mid \frac{-(1-x)^{\frac{1}{\alpha}-3} - ((\alpha-2)x((\alpha-2)(3\alpha-4)x+\alpha(39\alpha-55))+14)+\alpha(2(53-24\alpha)\alpha-73)+16)}{(\alpha-2)^2(3\alpha-4)} \right),
$$
system $\Upsilon_1^{-1} = Y_1^{-1}$ into a matrix of the form

$$ \left( \begin{array}{c} \mathbf{0} \\ G_1 \end{array} \right) := \left( \begin{array}{c} \mathbf{0} \\ g_1^{(1)} \\ g_2^{(1)} \\ \vdots \\ g_{n-1}^{(1)} \end{array} \right) = \left( \begin{array}{ccccc} 0 & 0 & \cdots & 0 & 0 \\ g_{1,1}^{(1)} & g_{1,2}^{(1)} & \cdots & g_{1,n-1}^{(1)} & g_{1,n}^{(1)} \\ g_{2,1}^{(1)} & g_{2,2}^{(1)} & \cdots & g_{2,n-1}^{(1)} & g_{2,n}^{(1)} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ g_{n-1,1}^{(1)} & g_{n-1,2}^{(1)} & \cdots & g_{n-1,n-1}^{(1)} & g_{n-1,n}^{(1)} \end{array} \right); \quad (47) $$

each of these row vectors $g_1^{(1)}, \dots, g_{n-1}^{(1)}$ is the gradient of a formal first integral and all $n-1$ of them are linearly independent. The next step would be to integrate the second variational system (from here onwards only quadratures are needed):

$$ Y_2(t) = Y_1 \int_0^t Y_1^{-1}(\tau) A_2(\tau) Y_1^{\odot 2}(\tau) d\tau $$

invert the fundamental matrix for $(\text{LVE}^2)$,

$$ \Upsilon_2 = \begin{pmatrix} Y_1^{\odot 2} & 0 \\ Y_2 & Y_1 \end{pmatrix}, \quad \text{thus} \quad \Psi_2 := (\Upsilon_2^{-1})^T $$

and multiply it by the second-order filter matrix

$$ \Phi_2 = \begin{pmatrix} F_1^{\odot 2} & 0 \\ F_2 & F_1 \end{pmatrix}, $$

where $F_2$ is obtained from (40), (41), and the result we obtain is

$$ \Phi_2 (\Psi_2)^T = \Phi_2 \Upsilon_2^{-1} = \begin{pmatrix} G_1^{\odot 2} & 0 \\ G_2 & G_1 \end{pmatrix} = \left( \begin{array}{c|c} \mathbf{0}_n & \mathbf{0}^{\odot 2} \\ \hline g_1^{(1)} & \vdots \\ \vdots & g_{n-1}^{(1)} \\ \hline g_1^{(2)} & \vdots \\ \vdots & g_{n-1}^{(2)} \end{array} \right) \quad (48) $$

Take for instance the first non-zero row of the lower block $(G_2 | G_1)$, i.e. the one highlighted in gray. Imagine we perform the computations for $(\text{LVE}_\phi^3)$ and multiply $\Upsilon_3^{-1}$ by the filter matrix $\Phi_3$; we obtain a matrix having the one in (48) as a lower right $D_{n,k} \times D_{n,k}$ block:

$$ \Phi_3 \Upsilon_3^{-1} = \begin{pmatrix} F_1^{\odot 3} & 0 & 0 \\ 3F_1 \odot F_2 & F_1^{\odot} & 0 \\ F_3 & F_2 & F_1 \end{pmatrix} \Upsilon_3^{-1} = \begin{pmatrix} G_1^{\odot 3} & 0 & 0 \\ 3G_1 \odot G_2 & G_1^{\odot} & 0 \\ G_3 & G_2 & G_1 \end{pmatrix} \quad (49) $$

$$ = \left( \begin{array}{c|c|c} \mathbf{0}_n & \mathbf{0}^{\odot 3} & \mathbf{0}^{\odot 2} \\ \hline g_1^{(1)} & g_1^{(2)} & \vdots \\ \vdots & \vdots & \vdots \\ g_{n-1}^{(1)} & g_{n-1}^{(2)} & \vdots \\ \hline g_{n-1}^{(3)} & g_1^{(2)} & g_1^{(1)} \\ \vdots & \vdots & \vdots \\ g_{n-1}^{(n-1)} & g_{n-1}^{(2)} & g_{n-1}^{(1)} \end{array} \right) $$
This can be achieved via (37) in Theorem 2.2.2 (the easiest option even for low values of $k$) or the filter matrix described in §2.3 to compute each $f_k$ in terms of the preceding $f_{k-1}, \dots, f_1$. We invite the reader to apply either method and realize that the above terms appear with only minor attention paid at each step.

A tedious scrutiny of the general form of these vectors shows that they are precisely limit cases of the higher-order Taylor terms (along $(x(t), 0) \in \Gamma$) of the following function:

$$f = \left( x \sin^{\frac{\alpha}{1-\alpha}}(y) - \left( \frac{\alpha \cos y}{\alpha - 1} \right) {}_2F_1 \left( \frac{1}{2}, \frac{\alpha}{2(\alpha-1)} + 1; \frac{3}{2}; \cos^2 y \right) \right)^{\frac{1}{\alpha}-1}$$

where ${}_2F_1$ is the **Gaussian hypergeometric function** [1]. Removing the outer power yields a first integral of (50) as well, as can be easily checked:

$$f = x \sin^{\frac{\alpha}{1-\alpha}}(y) - \left( \frac{\alpha \cos y}{\alpha - 1} \right) {}_2F_1 \left( \frac{1}{2}, \frac{\alpha}{2(\alpha-1)} + 1; \frac{3}{2}; \cos^2 y \right).$$

Undoing the transformation (49) we obtain:

**Theorem 3.1.1.** *The following function*

$$F = u^{\frac{\alpha}{1-\alpha}} v^{\frac{1}{\alpha-1}} \left( \frac{u^2}{v^2} + 1 \right)^{\frac{1}{2(\alpha-1)}} - \frac{\alpha v {}_2F_1 \left( \frac{1}{2}, \frac{\alpha}{2(\alpha-1)} + 1; \frac{3}{2}; \frac{v^2}{u^2+v^2} \right)}{(\alpha-1)\sqrt{u^2+v^2}}$$

is a first integral of the original special CDK system. □

## 3.2 The SIR model with vital dynamics

The SIR epidemiological model with vital dynamics [10, 17] is given by

$$\left\{ \begin{aligned} \dot{S} &= \mu(n-S) - \frac{\beta SI}{n}, \\ \dot{I} &= \frac{\beta SI}{n} - I(\gamma + \mu), \\ \dot{R} &= \gamma I - \mu R. \end{aligned} \right. \qquad (53)$$

The system was first introduced in [17] strictly for the modelling of infectious diseases, but has since been applied to a number of situations allowing dynamic compartmentalization, e.g. marketing [32]. The study of its integrability has been the subject of several works already ([6, 14]) and the search for a first integral is still an open problem. We will address specific examples; the general case will be addressed in the future.

First, let us start with a very simple example whose solution we already know: $\beta = \gamma$ and $\mu = 0$ and change of variable

$$S = n(1 + xy), \quad I = -ny, \quad R = \gamma n z, \qquad (54)$$

we obtain system of equations $\frac{d}{dt}(x, y, z) = (\gamma - \gamma(x-1)xy, \gamma xy^2, -y)$, and applying Conjectures 2.3.1, 2.3.4 and (if the matrix is used) Lemma 2.1.1, the general terms of the jets telegraph a very identifiable structure: that of two first integrals for the transformed system, $g^* = \frac{\log(xy+1)}{\gamma} + z$, $f^* = -xy + \log(xy+1) + y$ which transform back by (54) into $g = \frac{n \log(\frac{S}{n}) + R}{\gamma n}$, $f = -\frac{S+I}{n} + \log(\frac{S}{n}) + 1$, thus taking us to the already-known first-integral $S + I + R - n$ (thus, $S + I + R$) and the remaining one, $\frac{n \log(\frac{S}{n}) + R}{\gamma n}$ (alternatively, $\frac{S}{n} e^{\frac{R}{n}}$). This is all in keeping with what is already known about the integrability of the system for $\mu = 0$ (see e.g. [6, §1]). The reader can fill in the blanks for filter matrices or progressive filtering (Theorem 2.2.2) and formal Taylor terms here; we will add more detail about these in the less trivial case below.
thus we have an increasing row, highlighted in darker gray color $| \cdot |$ $g_1^{(3)}$ $g_1^{(2)}$ $g_1^{(1)}$. These vectors are precisely the terms of a formal first integral $g_1(z)$ of (1) along solution $\phi$:

$$g_1(z + \phi) = g_1(\phi) + g_1^{(1)}z + \frac{1}{2}g_1^{(2)}(\phi)z^{\odot 2} + \frac{1}{6}g_1^{(3)}(\phi)z^{\odot 3} + \dots \quad (50)$$

where $g_1^{(i)} = g_1^{(i)}(\phi)$, i.e. the implementation of (3) in 1.2.1 (2.) to $z = \phi$. The same applies, *mutatis mutandis*, to the rest of rows of infinite block $\dots G_3 | G_2 | G_1$, i.e. replacing $g_1^{(i)}$ with $g_j^{(i)}$ for some other $j = 2, \dots, n-1$, which will lead to another first integral $g_j$ in lieu of $g_1$ in (50). Needless to say the upper terms $3G_1 \odot G_2$, etcetera in (49) and higher-order counterparts contain all the cross-products appearing in derivatives of products of the first integrals, e.g. $g_1^{(1)} \odot g_3^{(2)}$ would be a row in $G_1 \odot G_2$ (it is easy to check which one by way of lexicographic ordering and we invite the reader to do this themselves, using the examples in §3).

# 3 Examples

## 3.1 Dixon's system

The following system,

$$\dot{u} = \frac{uv}{u^2 + v^2} - \alpha u, \qquad \dot{v} = \frac{v^2}{u^2 + v^2} - \beta v + \beta - 1, \quad (51)$$

where $\alpha, \beta > 0$, arose in [8] as a decoupled two-dimensional fragment of a three-dimensional dynamical model of the magnetic field of neutron stars. Several references (e.g. [2, 11, 38]) discussed the ostensible chaotic behavior of the dynamics of (51). [33], on the other hand, purported the non-chaotic character of this system based on a version of the Poincaré-Bendixson Theorem that precludes the need for compact sets comprising finitely many sets of homoclinic/heteroclinic connections. Given the existing, seldom-discussed gap between the concepts *non-chaotic* and *integrable* (regardless of which ad-hoc definition of chaos is used), we can afford to eschew the polemic altogether and focus on finding a conserved quantity for (51) using higher variational equations.

We will focus on the case $\alpha = \beta > 0$. The case $\alpha \neq \beta$, as well as the original three-dimensional model ([11, eqs (3)-(5)]), will be tackled in future work.

$\alpha = 1$ is immediate to solve and needs no further consideration. Let us focus on $\alpha \neq 1$. Previous attempts at simplifying higher-order LVE$^k$ suggest change of variable

$$(u, v) = \left( \frac{\sin y}{x}, \frac{\cos y}{x} \right) \qquad (52)$$

which transforms (51) for $\alpha = \beta$ into

$$\dot{x} = \alpha x (1 - x \cos y), \qquad \dot{y} = - (\alpha - 1) x \sin y. \quad (53)$$

We will use particular solution

$$\phi(t) = (x(t), 0) \quad \text{such that} \quad \dot{x} = -\alpha(-1+x)x. \quad (54)$$

Along $\Gamma$, the higher variational complex and its dual (LVE$^*_\phi$) feature the following blocks for
To wit, the case in which $\gamma = 0$ and $\beta, \mu \neq 0$. To the best of our knowledge, the integrability of this case is an open problem in earlier references. Change of variables

$$ (S, I, R) = \left( y^2 \left( x^{-\frac{\beta}{\mu}} + 1 \right) + n, -x^{-\frac{\beta}{\mu}} y^2, y^2 z \right), \qquad (55) $$

itself suggested by an earlier try at simplifying LVE$_{\phi}^{3}$, provides us with new system

$$ \dot{x} = - \frac{\mu x \left[ y^2 \left( x^{-\frac{\beta}{\mu}} + 1 \right) + n \right]}{n}, \quad \dot{y} = - \frac{1}{2} \mu y, \quad \dot{z} = 0. \qquad (56) $$

This system admits an invariant plane $y = z = 0$ containing particular solution

$$ \phi(t) = (\hat{x}(t), 0, 0), \quad \text{where} \quad \frac{d}{dt}\hat{x} = -\mu\hat{x}. \qquad (57) $$

This leads to the variational system $\dot{Y} = (A \odot \exp_{\odot} \mathrm{Id}_n) Y$ along $\phi$ where

$$ A_{\mathrm{LVE}_{\phi}} = A \odot \exp_{\odot} \mathrm{Id}_n = \begin{pmatrix} \cdots \\ \cdots 4A_1 \odot \mathrm{Id}_n^{\odot 3} \\ \cdots 6A_2 \odot \mathrm{Id}_n^{\odot 2} \\ \cdots 4A_3 \odot \mathrm{Id}_n \\ \cdots A_4 \end{pmatrix} \begin{pmatrix} 3A_1 \odot \mathrm{Id}_n^{\odot 2} \\ 3A_2 \odot \mathrm{Id}_n \\ 2A_1 \odot \mathrm{Id}_n \\ A_3 \\ A_2 \\ A_1 \end{pmatrix}, $$

and $A_i$ defined per (27) appear as follows:

$$ A_0 = \begin{pmatrix} -\mu \hat{x}(t) \\ 0 \\ 0 \end{pmatrix} $$

$$ A_1 = \begin{pmatrix} -\mu & 0 & 0 \\ 0 & -\frac{\mu}{2} & 0 \\ 0 & 0 & 0 \end{pmatrix} $$

$$ A_2 = \begin{pmatrix} 0 & 0 & 0 & \frac{2\mu\hat{x}(t)\left(-\hat{x}(t)^{-\frac{\beta}{\mu}} - 1\right)}{n} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix} $$

$$ A_3 = \begin{pmatrix} 0 & 0 & 0 & \frac{2\left(\left(\beta-\mu\right)\hat{x}(t)^{-\frac{\beta}{\mu}} - \mu\right)}{n} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix} $$

$$ A_k = \begin{pmatrix} 0 & 0 & 0 & \frac{2\left(\left(\beta-\mu\right)\hat{x}(t)^{-\frac{\beta}{\mu}} - \mu\right)}{n} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}, \quad k \ge 4 $$

Lower blocks of the quadrature-free filter matrix are built as proposed in §2.3, let us write the first few cases

$$ F_1 = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}, F_2 = \begin{pmatrix} 0 & 0 & \boxed{0_{3 \times 4}} \\ 0 & -\frac{1}{2\hat{x}_0} & 0 \\ 0 & 0 & 0 \end{pmatrix}, F_3 = \begin{pmatrix} 0 & 0 & \boxed{0_{3 \times 8}} \\ 0 & \frac{3}{4\hat{x}_0^2} & 0 \\ 0 & 0 & 0 \end{pmatrix}, F_4 = \begin{pmatrix} 0 & 0 & \boxed{0_{3 \times 4}} & 0 \\ 0 & -\frac{15}{8\hat{x}_0^3} & \frac{3\left(\hat{x}_0^{-\frac{\beta}{\mu}} + 1\right)}{n\hat{x}_0} & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix} $$
$k \ge 0,$

$$
A_k = \left\{
\begin{array}{ll}
X(x,y) = \begin{pmatrix} \alpha x (1 - x \cos y) \\ -(\alpha - 1) x \sin y \end{pmatrix}, & k = 0, \\[2ex]
\begin{pmatrix} -\alpha (2x - 1) & 0 \\ 0 & -(\alpha - 1) x \end{pmatrix}, & k = 1, \\[2ex]
\left( \begin{array}{cc} (-1)^{k/2} 2\alpha & 0 \\ 0 & (-1)^{k/2} (-1 + \alpha) \end{array} \right), & k \text{ even} \\[2ex]
\left( \begin{array}{cc} (-1)^{\frac{k+1}{2}} 2\alpha x & 0 \\ 0 & (-1)^{\frac{k+1}{2}} (-1 + \alpha) x \end{array} \right), & k \text{ odd},
\end{array}
\right.
\quad (55)
$$

for instance a principal fundamental matrix of $(VE_\phi)$ is

$$
\Upsilon_1 = Y_1 = \begin{pmatrix} \frac{(x-1)x}{(x_0-1)x_0} & 0 \\ 0 & (1-x_0)^{\frac{1}{\alpha}-1}(1-x)^{1-\frac{1}{\alpha}} \end{pmatrix}
$$

and the first-order filter matrix is very simple in this case,

$$
F_1 = \mathrm{Id}_2 - \frac{1}{X_1^0} (X^0 \odot e_1^T) = \begin{pmatrix} 0 & 0 \\ -\frac{X_2(x_0, y_0)}{X_1(x_0, y_0)} & 1 \end{pmatrix} = \begin{pmatrix} 0 & 0 \\ 0 & 1 \end{pmatrix},
$$

where $x_0 = x(0)$ and $y_0 = y(0)$; so are the higher-order matrices in $\Phi = \exp_{\odot} F$ obtained from (40), some of which are shown below:

$$
F_2 = \begin{pmatrix} 0 & 0 & 0 \\ 0 & \frac{1-\alpha}{\alpha(x_0-1)} & 0 \end{pmatrix}, \quad F_3 = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & \frac{(\alpha-1)(2\alpha-1)}{\alpha^2(x_0-1)^2} & 0 & 0 \end{pmatrix}, \\
F_4 = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & -\frac{(\alpha-1)(2\alpha-1)(3\alpha-1)}{\alpha^3(x_0-1)^3} & 0 & -\frac{(\alpha-1)(2x_0+1)}{\alpha(x_0-1)^2} \end{pmatrix}, \dots
$$

Returning to $k=1$, Lemma 2.2.1 ensures that the second column of the filtered matrix

$$
(Y_1^{-1})^T F_1^T = \begin{pmatrix} 0 & 0 \\ 0 & (1-x_0)^{1-\frac{1}{\alpha}}(1-x)^{\frac{1}{\alpha}-1} \end{pmatrix},
$$

is an admissible solution of degree 1. We can normalize this to eliminate appearances of $x_0$:

$$
(\Phi_1^{-1})^T F_1^T P^T := \begin{pmatrix} 0 & 0 \\ 0 & (1-x_0)^{1-\frac{1}{\alpha}}(1-x)^{\frac{1}{\alpha}-1} \end{pmatrix} \begin{pmatrix} 1 & 0 \\ 0 & (1-x_0)^{\frac{1}{\alpha}-1} \end{pmatrix}^T = \begin{pmatrix} 0 & 0 \\ 0 & (1-x)^{\frac{1}{\alpha}-1} \end{pmatrix},
$$

the second of whose columns is still the gradient $f^{(1)}$ of a formal first integral (this is exactly
what appears as a row $g_1^{(1)}$ in (47)).

The higher order trimming (in which we discard redundant symmetric products $f^{(i_1)} \odot \dots \odot f^{(i_k)}$, i.e. derivatives of powers of $f$) is clear once we write down the dual system of equations and filter out telescopically using either Lemma 2.1.1 and the general form of that each new term $f^{(k)}$ of our formal first integral $f$ is expressible as the following (bear in mind that $x = x(t)$ here): $f^{(k)}(x) = (\frac{d}{dx}f^{(k-1)}|_{g_k(x)})$, where

$$
\sum_{j=0}^{\frac{k-1}{2}} \left( (-1)^{j+1} (\alpha - 1) \binom{k}{-1-2j+k} g_{k-2j}(x) + (-1)^{j+1} \alpha \binom{k}{-2j+k} x g'_{k-2j}(x) \right) + \alpha g'_k(x) = 0,
$$
which means $g_k(x) = 0$ for every even $k \ge 2$ and e.g.

$$f_1 = \left( 0 \ (1-x)^{\frac{1}{\alpha}-1} \right),$$

$$f_2 = \left( 0 \ \frac{(\alpha-1)(1-x)^{\frac{1}{\alpha}-2}}{\alpha} \ 0 \right),$$

$$f_3 = \left( 0 \ (\frac{1}{\alpha}-2) (\frac{1}{\alpha}-1) (1-x)^{\frac{1}{\alpha}-3} \ 0 \ \frac{(1-x)^{\frac{1}{\alpha}-2}(2\alpha+(\alpha-2)x-1)}{\alpha-2} \right),$$

$$f_4 = \left( \frac{d}{dx} f_3 \mid 0 \right),$$

$$f_5 = \left( \frac{d}{dx} f_4 \mid \frac{-(1-x)^{\frac{1}{\alpha}-3} - (-(\alpha-2)x((\alpha-2)(3\alpha-4)x+\alpha(39\alpha-55)+14)+\alpha(2(53-24\alpha)\alpha-73)+16))}{(\alpha-2)^2(3\alpha-4)} \right),$$

This can be achieved via (38) in Theorem 2.2.2 (the easiest option even for low values of $k$) or the filter matrix described in §2.3 to compute each $f_k$ in terms of the preceding $f_{k-1}, \dots, f_1$. We invite the reader to apply either method and realize that the above terms appear with only minor attention paid at each step.

A tedious scrutiny of the general form of these vectors shows that they are precisely limit cases of the higher-order Taylor terms (along $(x(t), 0) \in \Gamma$) of the following function:

$$f = \left( x \sin^{\frac{\alpha}{1-\alpha}}(y) - \left( \frac{\alpha \cos y}{\alpha - 1} \right) {}_2F_1 \left( \frac{1}{2}, \frac{\alpha}{2(\alpha-1)} + 1; \frac{3}{2}; \cos^2 y \right) \right)^{\frac{1}{\alpha}-1}$$

where ${}_2F_1$ is the **Gaussian hypergeometric function** [1]. Removing the outer power yields a first integral of (53) as well, as can be easily checked:

$$f = x \sin^{\frac{\alpha}{1-\alpha}}(y) - \left( \frac{\alpha \cos y}{\alpha - 1} \right) {}_2F_1 \left( \frac{1}{2}, \frac{\alpha}{2(\alpha-1)} + 1; \frac{3}{2}; \cos^2 y \right).$$

Undoing the transformation (52) we obtain:

**Theorem 3.1.1.** *The following function*

$$F = u^{-\frac{\alpha}{\alpha-1}} (u^2 + v^2)^{\frac{1}{2(\alpha-1)}} - \frac{\alpha v {}_2F_1\left(\frac{1}{2}, \frac{\alpha}{2(\alpha-1)} + 1; \frac{3}{2}; \frac{v^2}{u^2 + v^2}\right)}{(\alpha-1)\sqrt{u^2 + v^2}}$$

is a first integral of the original special CDK system. □

## 3.2 The SIR model with vital dynamics

The SIR epidemiological model with vital dynamics [10, 17] is given by

$$\left\{ \begin{array}{l} \dot{S} = \mu(n-S) - \frac{\beta SI}{n}, \\ \dot{I} = \frac{\beta SI}{n} - I(\gamma + \mu), \\ \dot{R} = \gamma I - \mu R. \end{array} \right. \qquad (56)$$

The system was first introduced in [17] strictly for the modelling of infectious diseases, but has since been applied to a number of situations allowing dynamic compartmentalization, e.g. marketing [32]. The study of its integrability has been the subject of several works already ([6, 14]) and the search for a first integral is still an open problem. We will address specific examples; the general case will be addressed in the future.

First, let us start with a very simple example whose solution we already know: $\beta = \gamma$ and $\mu = 0$ and change of variable

$$S = n(1 + xy), \quad I = -ny, \quad R = \gamma n z, \qquad (57)$$
where $P_1$ is as in (58) and $P_2, P_3$ are successively chosen so as vanish $\hat{x}_0$ from the matrix (in this case, all terms in $P_2, P_3$ equal to zero except for $(P_3)_{2,7} = \frac{3\left(\mu\left(-\hat{x}_0^{-\frac{\beta}{\mu}}-1\right)+\beta\right)}{n\sqrt{\hat{x}_0}(\beta-\mu)})$. This equates to weeding out copies of successive symmetric products of $f^{(i)}$, and/or $g^{(j)}$ and harks back to Lemma 2.1.1 and Conjecture 2.3.4. The resulting matrix retains its neat echeloned structure

$$ \Pi_3 \Phi_3 \Upsilon_3^{-1} = \begin{pmatrix} G_1^{\odot 3} & 0 & 0 \\ 3G_1 \odot G_2 & G_1^{\odot 2} & 0 \\ G_3 & G_2 & G_1 \end{pmatrix} $$

where lower block $G_3 | G_2 | G_1$ has independent formal first integral jets as its non-zero rows:

$$ \left( \begin{array}{ccccccccc} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{3}{4\hat{x}^{5/2}} & 0 & 0 & 0 & 0 & \frac{3\left(-\mu\hat{x}^{-\frac{\beta}{\mu}}+\beta-\mu\right)}{(\beta-\mu)n\sqrt{\hat{x}}} & 0 & 0 & 0 & 0 & -\frac{1}{2\hat{x}^{3/2}} & 0 & 0 & 0 & 0 & 0 & \frac{1}{\sqrt{\hat{x}}} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{array} \right) \tag{59} $$

in other words, the following are degree-three asymptotic expressions of formal first integrals ($z = (x, y, z)$ and $\phi = (\hat{x}, 0, 0)$ is the particular solution as usual):

$$ \begin{aligned} f(z + \phi) &= f^{(1)}(\phi)z + \frac{1}{2}f^{(2)}(\phi)z^{\odot 2} + \frac{1}{3!}f^{(3)}z^{\odot 3} + O(z^{\odot 4}) \\ &= \frac{y}{\sqrt{\hat{x}}} - \frac{xy}{2\hat{x}^{3/2}} + \frac{1}{6}\left(\frac{3y^3\left(\frac{\mu\hat{x}^{1-\frac{\beta}{\mu}}}{\mu-\beta} + \hat{x}\right)}{n\hat{x}^{3/2}} + \frac{9x^2y}{4\hat{x}^{5/2}}\right) + O(z^{\odot 4}), \\ g(z + \phi) &= g^{(1)}(\phi)z + \frac{1}{2}g^{(2)}(\phi)z^{\odot 2} + \frac{1}{3!}g^{(3)}z^{\odot 3} + O(z^{\odot 4}) = z + O(z^{\odot 4}). \end{aligned} $$

Again, Theorem 2.2.2 does the exact same job faster and with the guarantee that it works, e.g.

$$ \begin{aligned} f^{(2)} &= -\left[ (Y_1^{-1})^{\odot 2} \right]^T \left[ \int_{\hat{x}_0}^{\hat{x}} (Y_1^{\odot 2})^T f^{(1)} A_2 dx - \left( 0, \frac{1}{2\hat{x}_0^{3/2}}, 0, 0, 0, 0 \right) \right], \\ f^{(3)} &= -\left[ (Y_1^{-1})^{\odot 3} \right]^T \left[ \int_{\hat{x}_0}^{\hat{x}} (Y_1^{\odot 3})^T \left[ 3f^{(2)}(A_2 \odot \mathrm{Id}_3) + f^{(1)}A_3 \right] dx - \left( 0, -\frac{3}{4\hat{x}_0^{5/2}}, 0, 0, 0, 0, 0, 0, 0 \right) \right], \end{aligned} $$

and the reader can check that the jet terms are exactly the same as those in (59). Either method may sound like a cumbersome process but it only requires special attention in its first stages if savvy simplifications are chosen, as was the case for Dixon's system. Later stages of the bottom row remain zero and accordingly $g = z$ is easily checked to be a first integral (formal and convergent) for (56); for $f$, and much like in §3.1, pattern checks along with simple resummation yield a first integral for system (56),

$$ f = \frac{ye^{\frac{y^2}{2n}}}{\left[ e^{\frac{uy^2}{n}} \left( \frac{uy^2}{n} \right)^u \Gamma \left( 1 - u, \frac{uy^2}{n} \right) + x^u \right]^{\frac{1}{2u}}}, $$

where $u = \frac{\beta}{\mu}$ and $\Gamma(\cdot, \cdot)$ is the **incomplete Gamma function** [1]; undoing transformation (55) we obtain something that is now immediate to prove on its own:
we obtain system of equations $\frac{d}{dt}(x, y, z) = (\gamma - \gamma(x - 1)xy, \gamma xy^2, -y)$, and applying Theorem 2.3.1 and Corollary 2.3.5 (and if the matrix is used) Lemma 2.1.1, the general terms of the jets telegraph a very identifiable structure: that of two first integrals for the transformed system, $g^* = \frac{\log(xy+1)}{\gamma} + z$, $f^* = -xy + \log(xy+1) + y$ which transform back by (57) into $g = \frac{n \log(\frac{S}{n}) + R}{\gamma n}$, $f = -\frac{S+I}{n} + \log(\frac{S}{n}) + 1$, thus taking us to the already-known first-integral $S+I+R-n$ (thus, $S+I+R$) and the remaining one, $\frac{n \log(\frac{S}{n}) + R}{\gamma n}$ (alternatively, $\frac{S}{n}e^{\frac{R}{n}}$). This is all in keeping with what is already known about the integrability of the system for $\mu = 0$ (see e.g. [6, §1]). The reader can fill in the blanks for filter matrices or progressive filtering (Theorem 2.2.2) and formal Taylor terms here; we will add more detail about these in the less trivial case below.

To wit, the case in which $\gamma = 0$ and $\beta, \mu \neq 0$. To the best of our knowledge, the integrability of this case is an open problem in earlier references. Change of variables

$$ (S, I, R) = \left( y^2 \left( x^{-\frac{\beta}{\mu}} + 1 \right) + n, -x^{-\frac{\beta}{\mu}} y^2, y^2 z \right), \qquad (58) $$

itself suggested by an earlier try at simplifying LVE$_\phi^3$, provides us with new system

$$ \dot{x} = -\frac{\mu x \left[ y^2 \left( x^{-\frac{\beta}{\mu}} + 1 \right) + n \right]}{n}, \quad \dot{y} = -\frac{1}{2} \mu y, \quad \dot{z} = 0. \qquad (59) $$

This system admits an invariant plane $y = z = 0$ containing particular solution

$$ \phi(t) = (\hat{x}(t), 0, 0), \quad \text{where} \quad \frac{d}{dt}\hat{x} = -\mu\hat{x}. \qquad (60) $$

This leads to the variational system $\dot{Y} = (A \odot \exp_{\odot} \mathrm{Id}_n) Y$ along $\phi$ where

$$ A_{\mathrm{LVE}_{\phi}} = A \odot \exp_{\odot} \mathrm{Id}_n = \begin{pmatrix} \cdots \\ \cdots 4A_1 \odot \mathrm{Id}_n^{\odot 3} \\ \cdots 6A_2 \odot \mathrm{Id}_n^{\odot 2} & 3A_1 \odot \mathrm{Id}_n^{\odot 2} \\ \cdots 4A_3 \odot \mathrm{Id}_n & 3A_2 \odot \mathrm{Id}_n & 2A_1 \odot \mathrm{Id}_n \\ \cdots A_4 & A_3 & A_2 & A_1 \end{pmatrix}, $$

and $A_i$ defined per (27) appear as follows:

$$ \begin{align*} A_0 &= \begin{pmatrix} -\mu \hat{x}(t) \\ 0 \\ 0 \end{pmatrix} \\ A_1 &= \begin{pmatrix} -\mu & 0 & 0 \\ 0 & -\frac{\mu}{2} & 0 \\ 0 & 0 & 0 \end{pmatrix} \\ A_2 &= \begin{pmatrix} 0 & 0 & 0 & \frac{2\mu\hat{x}(t)\left(-\hat{x}(t)^{-\frac{\beta}{\mu}} - 1\right)}{n} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix} \\ A_3 &= \begin{pmatrix} 0 & 0 & 0 & \frac{2\left(\left(\beta-\mu\right)\hat{x}(t)^{-\frac{\beta}{\mu}} - \mu\right)}{n} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix} \\ A_k &= \begin{pmatrix} 0 & 0 & 0 & \frac{2\left(\left(\beta-\mu\right)\hat{x}(t)^{-\frac{\beta}{\mu}} - \mu\right)}{n} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}, \quad k \ge 4 \end{align*} $$
**Theorem 3.2.1.** The following two functions,

$$F = \frac{e^{-\frac{n+S+I}{2n}} \sqrt{-n+S+I}}{\left(\frac{-n+S+I}{I} - \left(\frac{\beta}{\mu n}\right)^{\beta/\mu} e^{\frac{\beta(-n+S+I)}{\mu n}} (-n+S+I)^{\beta/\mu} \Gamma\left(1-\frac{\beta}{\mu}, \frac{\beta(-n+S+I)}{\mu n}\right)\right)^{\frac{\mu}{2\beta}}},$$
$$G = \frac{R}{-n+S+I}.$$

are two functionally independent first integrals of SIR with $\gamma = 0$ (53). $\square$

It is worth noting, however, that the only easy goal we can aspire to in most cases will be that of formal first integrals with a compact identifiable general term but without a straightforward convergent expression, and resummation and deep pattern recognition techniques will need to be adapted to each problem. This will also be the subject of future research in the short term but for the time being let us see a short example below.

## 3.3 The van der Pol oscillator

This well-known system is usually presented as an epitome of the proverbial “simple yet chaotic” dynamical system [37]

$$\ddot{u} = \mu(1 - u^2)\dot{u} - u \quad (60)$$

Value $\mu = 2$ appears to have a significant, yet still unclarified role in potential qualitative simplifications in variable transformations. Thus let us fix this value for the present paper:

$$\dot{u} = v, \quad \dot{v} = 2(1 - u^2)v - u \quad (61)$$

Perform the change of variables

$$(u(t), v(t)) = \left( \frac{x(t)y(t)}{\sqrt{2}}, \left( \frac{x(t)}{\sqrt{2}} - \frac{1}{\sqrt{2}} \right) y(t) \right) \quad (62)$$

that will simplify the corresponding jet, and our transformed system is

$$\dot{x} = -(x-1)x^3 y^2 - 1, \quad \dot{y} = (x-1)x^2 y^3 + y. \quad (63)$$

The dual infinite variational system $(VE_\Gamma)^*$ along particular solution

$$\Gamma = \{(x(t), 0) : t \in \mathbb{C}\}, \quad \text{where } \dot{x}(t) = -1, \quad (64)$$

has the following very simple matrices:

$$\begin{align*}
A_1 &= \begin{pmatrix} 0 & 0 \\ 0 & 1 \end{pmatrix}, \\
A_2 &= \begin{pmatrix} 0 & 0 & -2(x(t) - 1)x(t)^3 \\ 0 & 0 & 0 \end{pmatrix}, \\
A_3 &= \begin{pmatrix} 0 & 0 & 2(3 - 4x(t))x(t)^2 & 0 \\ 0 & 0 & 0 & 6(x(t) - 1)x(t)^2 \end{pmatrix}, \\
A_4 &= \begin{pmatrix} 0 & 0 & 12(1 - 2x(t))x(t) & 0 & 0 \\ 0 & 0 & 0 & 6x(t)(3x(t) - 2) & 0 \end{pmatrix}, \\
A_5 &= \begin{pmatrix} 0 & 0 & 12 - 48x(t) & 0 & 0 & 0 \\ 0 & 0 & 0 & 36x(t) - 12 & 0 & 0 \end{pmatrix}, \\
A_6 &= \begin{pmatrix} 0 & 0 & -48 & 0 & 0 & 0 \\ 0 & 0 & 0 & 36 & 0 & 0 \end{pmatrix}, \\
A_k &= 0_{2 \times d_{2,k}}, \quad k \ge 7.
\end{align*}$$
The fact that the ODE in (64) is so trivial to solve ($x(t) = -t + C$) adds no simplicity to what follows below. The fundamental matrix of $VE_1$ and the first filter matrix are equally simple:

$$Y_1 = \begin{pmatrix} 1 & 0 \\ 0 & e^{x(0)-x(t)} \end{pmatrix}, \quad F_1 = \begin{pmatrix} 0 & 0 \\ 0 & 1 \end{pmatrix}$$

and the generic form of the admissible solutions after applying Theorem 2.2.2 is the following: $f_k = (\frac{d}{dx} f_{k-1} | g_k(x))$ where $g_k$ is a solution to

$$-(k-1)k(x-1)x^3 g'_{k-2}(x) - g'_k(x) + (k-2)(k-1)k(x-1)x^2 g_{k-2}(x) + kg_k(x) = 0. \quad (65)$$

and moreover $g_k(x) = 0$ if $k$ is even, and $g_1(x) = e^x$. Equation (65) has a solution (which of course only merits being written if $k$ is odd)

$$g_k(x) = e^{kx} \int (k-1)k(\xi-1)\xi^2 e^{-k\xi} \left( (k-2)g_{k-2}(\xi) - \xi g'_{k-2}(\xi) \right) d\xi.$$

The formal first integral of (63) takes the form

$$F = \sum_{k=1}^{\infty} \frac{1}{k!} g_k(x) y^k = \sum_{i=0}^{\infty} (-1)^i e^{(2i+1)x} y^{2i+1} G_{2i+1};$$

the first terms of $G_\star$ are easy to discern,

$$G_1 = 1,$$

$$G_3 = \int_{x_0}^{x} e^{-2\xi}(\xi-1)^2 \xi^2 d\xi,$$

$$G_5 = \int_{x_0}^{x} (\xi - 1)^3 e^{-4\xi} \xi^5 d\xi + 3 \int_{\tau_0}^{\tau} (\xi - 1)^2 e^{-2\xi} \xi^2 \left( \int (\eta - 1)^2 e^{-2\eta} \eta^2 d\eta \right) d\xi$$

and the general term is

$$G_{2i+1} = \sum_{m=1}^{i} \sum_{j=1}^{\binom{i-1}{m-1}} D_{i-1, m-1, j} \Phi \left( (i-m) \mathbf{C}_{m, i-m-1, j} + \mathbf{1}_m, 2(i-m) \mathbf{C}_{m, i-m-1, j} + \mathbf{1}_m, 3(i-m) \mathbf{C}_{m, i-m-1, j} + \mathbf{1}_m \right)$$

where we define the following constant vectors: $\mathbf{1}_n := (1, \dots, 1) \in \mathbb{Z}^n$, and $C_{\star, \star, \star}$ are the columns of the symmetric product

$$\mathrm{Id}_k \odot \mathbf{1}_{d_{k,l}}^T = (C_{k,l,1} | C_{k,l,2} | \dots | C_{k,l,d_{k,l+1}}) \in \mathrm{Mat}_{k,d_{k,l+1}}(\mathbb{Z});$$

(obviously $C_{k,l,j} = 0$ for $l < 0$), $\Phi$ is defined recursively as follows: for every $\mathbf{a}, \mathbf{b}, \mathbf{c} \in \mathbb{Z}^m$,

$$\Phi(\mathbf{a}, \mathbf{b}, \mathbf{c}) := \begin{cases} \int_{x_0}^{x} (\tau - 1)^{a_1+1} e^{-(b_1+1)\tau} \tau^{c_1} d\tau, & m=1, \\ \int_{x_0}^{x} (\tau - 1)^{a_1+1} e^{-(b_1+1)\tau} \tau^{c_1} \Phi((a_2, \dots, a_m), (b_2, \dots, b_m), (c_2, \dots, c_m)) d\tau, & m>1 \end{cases}$$

and $D_{n,m,j}$ is the $j^{th}$ term in the vector $\mathbf{D}_{n,m}$ obtained from discarding terms containing repeating factors from the reverse-order version of $\mathbf{i}_n^T \odot \cdots \odot \mathbf{i}_n^T$ and $\mathbf{i}_n$ is the reverse-order sequence of the first $n$ odd integers $> 1$, e.g.

$$\mathbf{D}_{3,2} = (3 \cdot 5, 3 \cdot 7, 5 \cdot 7) = (D_{3,2,j})_{j=1,2,3,4}, \quad \mathbf{D}_{4,3} = (3 \cdot 5 \cdot 7, 3 \cdot 5 \cdot 9, 3 \cdot 7 \cdot 9, 5 \cdot 7 \cdot 9)$$
The bottom two rows are already admissible jets $\star^{(3)}$ $| \star^{(2)}$ $| \star^{(1)}$, but we can trim them further by multiplying this matrix by

$$ \Pi_3 = \begin{pmatrix} P_1^{\odot 3} & 0 & 0 \\ 3P_1 \odot P_2 & P_1^{\odot 2} & 0 \\ P_3 & P_2 & P_1 \end{pmatrix}^T $$

where $P_1$ is as in (61) and $P_2, P_3$ are successively chosen so as vanish $\hat{x}_0$ from the matrix (in this case, all terms in $P_2, P_3$ equal to zero except for $(P_3)_{2,7} = \frac{3\left(\mu\left(-\hat{x}_0^{-\frac{\beta}{\mu}} - 1\right) + \beta\right)}{n\sqrt{\hat{x}_0}(\beta - \mu)}$). This equates to weeding out copies of successive symmetric products of $f^{(i)}$, and/or $g^{(j)}$ and harks back to Lemma 2.1.1 and Conjecture 2.3.5. The resulting matrix retains its neat echeloned structure

$$ \Pi_3 \Phi_3 \Upsilon_3^{-1} = \begin{pmatrix} G_1^{\odot 3} & 0 & 0 \\ 3G_1 \odot G_2 & G_1^{\odot 2} & 0 \\ G_3 & G_2 & G_1 \end{pmatrix} $$

where lower block $G_3 | G_2 | G_1$ has independent formal first integral jets as its non-zero rows:

$$ \left( \begin{array}{ccccccccc} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{3}{4\hat{x}^{5/2}} & 0 & 0 & 0 & 0 & \frac{3\left(-\mu\hat{x}^{-\frac{\beta}{\mu}} + \beta - \mu\right)}{(\beta - \mu)n\sqrt{\hat{x}}} & 0 & 0 & 0 & 0 & -\frac{1}{2\hat{x}^{3/2}} & 0 & 0 & 0 & 0 & 0 & \frac{1}{\sqrt{\hat{x}}} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{array} \right) \quad (62) $$

in other words, the following are degree-three asymptotic expressions of formal first integrals ($z = (x, y, z)$ and $\phi = (\hat{x}, 0, 0)$ is the particular solution as usual):

$$ \begin{aligned} f(z + \phi) &= f^{(1)}(\phi)z + \frac{1}{2}f^{(2)}(\phi)z^{\odot 2} + \frac{1}{3!}f^{(3)}z^{\odot 3} + O(z^{\odot 4}) \\ &= \frac{y}{\sqrt{\hat{x}}} - \frac{xy}{2\hat{x}^{3/2}} + \frac{1}{6}\left(\frac{3y^3\left(\frac{\mu\hat{x}^{1-\frac{\beta}{\mu}}}{\mu-\beta} + \hat{x}\right)}{n\hat{x}^{3/2}} + \frac{9x^2y}{4\hat{x}^{5/2}}\right) + O(z^{\odot 4}), \\ g(z + \phi) &= g^{(1)}(\phi)z + \frac{1}{2}g^{(2)}(\phi)z^{\odot 2} + \frac{1}{3!}g^{(3)}z^{\odot 3} + O(z^{\odot 4}) = z + O(z^{\odot 4}). \end{aligned} $$

Again, Theorem 2.2.2 does the exact same job faster and with the guarantee that it works, e.g.

$$ \begin{aligned} f^{(2)} &= -\left[ (Y_1^{-1})^{\odot 2} \right]^T \left[ \int_{\hat{x}_0}^{\hat{x}} (Y_1^{\odot 2})^T f^{(1)} A_2 dx - \left( 0, \frac{1}{2\hat{x}_0^{3/2}}, 0, 0, 0, 0 \right) \right], \\ f^{(3)} &= -\left[ (Y_1^{-1})^{\odot 3} \right]^T \left[ \int_{\hat{x}_0}^{\hat{x}} (Y_1^{\odot 3})^T \left[ 3f^{(2)}(A_2 \odot \mathrm{Id}_3) + f^{(1)}A_3 \right] dx - \left( 0, -\frac{3}{4\hat{x}_0^{5/2}}, 0, 0, 0, 0, 0, 0, 0 \right) \right], \end{aligned} $$

and the reader can check that the jet terms are exactly the same as those in (62). Either method may sound like a cumbersome process but it only requires special attention in its first stages if savvy simplifications are chosen, as was the case for Dixon's system. Later stages of the bottom row remain zero and accordingly $g = z$ is easily checked to be a first integral (formal and convergent) for (59); for $f$, and much like in §3.1, pattern checks along with simple resummation yield a first integral for system (59),

$$ f = \frac{ye^{\frac{y^2}{2n}}}{\left[ e^{\frac{uy^2}{n}} \left( \frac{uy^2}{n} \right)^u \Gamma \left( 1 - u, \frac{uy^2}{n} \right) + x^u \right]^{\frac{1}{2u}}}, $$
# 4 Further comments

Ever since it first took form in [3], the use of linearized higher variational equations to characterize, obstruct or redefine integrability is proving itself to be a welcome addition to the domain of dynamical systems. This contribution arrives at a time in which other attempts (e.g. [24, 25]) at finding methods to integrate dynamical systems rather than proving them non-integrable are gradually gaining traction.

A number of facts are germane to further studies:

(i) the potential shown by higher variational equations in finding adequate changes of variables that will simplify the system into tractability and, potentially, integrability (this is most exemplified in the sample systems tackled in §3);

(ii) linked to the above, the possible application of a *Baker-Campbell-Hausdorff* [12] sorts of formula upon variable transformations, using the decomposition of a transformed fundamental matrix into a product of matrices to our advantage;

(iii) the possible amenability of results in §2 above to machine learning [26];

(iv) the applicability of this automatic quadrature-algebraic method to the application of generalized, multivariate *Padé approximants* [9, 13] to first integrals, which will be tackled in an immediate future;

(v) the fact that the autonomous system (1) need not be Hamiltonian, thus allowing any arbitrary system to be considered and losing the serious constraint posed by symplectic transformations, since variable transformations can now be chosen freely;

(vi) also related to (v): the fact that this overrides the obstacles posed by the direct study of variational systems, which only has an implementation to Hamiltonian systems and requires first integrals to be meromorphic.

# References

[1] M. Abramowitz and I. A. Stegun (eds.), *Handbook of mathematical functions with formulas, graphs, and mathematical tables*, A Wiley-Interscience Publication, John Wiley & Sons Inc., New York, 1984, Reprint of the 1972 edition, Selected Government Publications.

[2] J. Alvarez-Ramirez, J. Delgado-Fernandez and G. Espinosa-Paredes *The origin of a continuous two-dimensional "chaotic" dynamics*, Int. J. Bif. Chaos 15, 3023-3029, 2005.

[3] A. Aparicio-Monforte, M. Barkatou, S. Simon, J.-A. Weil, *Formal first integrals along solutions of differential systems I*, ISSAC 2011 Proceedings of the 36th International Symposium on Symbolic and Algebraic Computation, 19-26, ACM, New York, 2011.

[4] U. Bekbaev, *A matrix representation of composition of polynomial maps*, arXiv:0901.3179v3 [math.AC] 22 Sep 2009.

[5] U. Bekbaev, *Matrix representations for symmetric and antisymmetric multi-linear maps*, arXiv:1010.2579v1 [math.RA] 13 Oct 2010

[6] Chen, Ding, *On the Integrability of the SIR Epidemic Model with Vital Dynamics*, Hindawi Advances in Mathematical Physics, Volume 2020, Article ID 5869275, https://doi.org/10.1155/2020/5869275
where $u = \frac{\beta}{\mu}$ and $\Gamma(\cdot, \cdot)$ is the **incomplete Gamma function** [1]; undoing transformation (58) we obtain something that is now immediate to prove on its own:

**Theorem 3.2.1.** The following two functions,

$$F = \frac{e^{-\frac{n+S+I}{2n}} \sqrt{-n+S+I}}{\left(\frac{-n+S+I}{I} - \left(\frac{\beta}{\mu n}\right)^{\beta/\mu} e^{\frac{\beta(-n+S+I)}{\mu n}} (-n+S+I)^{\beta/\mu} \Gamma\left(1-\frac{\beta}{\mu}, \frac{\beta(-n+S+I)}{\mu n}\right)\right)^{\frac{\mu}{2\beta}}},$$
$$G = \frac{R}{-n+S+I}.$$

are two functionally independent first integrals of SIR with $\gamma = 0$ (56). $\square$

It is worth noting, however, that the only easy goal we can aspire to in most cases will be that of formal first integrals with a compact identifiable general term but without a straightforward convergent expression, and resummation and deep pattern recognition techniques will need to be adapted to each problem. This will also be the subject of future research in the short term; for the time being let us see a short example below.

## 3.3 The van der Pol oscillator

This well-known system is usually presented as an epitome of the proverbial “simple yet chaotic” dynamical system [38]

$$\ddot{u} = \mu(1 - u^2)\dot{u} - u \quad (63)$$

Value $\mu = 2$ appears to have a significant, yet still unclarified role in potential qualitative simplifications in variable transformations. Thus let us fix this value for the present paper:

$$\dot{u} = v, \quad \dot{v} = 2(1 - u^2)v - u \quad (64)$$

Perform the change of variables

$$(u(t), v(t)) = \left( \frac{x(t)y(t)}{\sqrt{2}}, \left( \frac{x(t)}{\sqrt{2}} - \frac{1}{\sqrt{2}} \right) y(t) \right) \quad (65)$$

that will simplify the corresponding jet, and our transformed system is

$$\dot{x} = -(x-1)x^3 y^2 - 1, \quad \dot{y} = (x-1)x^2 y^3 + y. \quad (66)$$

The dual infinite variational system $(VE_\Gamma)^*$ along particular solution

$$\Gamma = \{(x(t), 0) : t \in \mathbb{C}\}, \quad \text{where} \quad \dot{x}(t) = -1, \quad (67)$$

has the following very simple matrices:

$$A_1 = \begin{pmatrix} 0 & 0 \\ 0 & 1 \end{pmatrix},$$

$$A_2 = \begin{pmatrix} 0 & 0 & -2(x(t) - 1)x(t)^3 \\ 0 & 0 & 0 \end{pmatrix},$$

$$A_3 = \begin{pmatrix} 0 & 0 & 2(3 - 4x(t))x(t)^2 & 0 \\ 0 & 0 & 0 & 6(x(t) - 1)x(t)^2 \end{pmatrix},$$

$$A_4 = \begin{pmatrix} 0 & 0 & 12(1 - 2x(t))x(t) & 0 & 0 \\ 0 & 0 & 0 & 6x(t)(3x(t) - 2) & 0 \end{pmatrix},$$

$$A_5 = \begin{pmatrix} 0 & 0 & 12 - 48x(t) & 0 & 0 & 0 \\ 0 & 0 & 0 & 36x(t) - 12 & 0 & 0 \end{pmatrix},$$

$$A_6 = \begin{pmatrix} 0 & 0 & -48 & 0 & 0 & 0 \\ 0 & 0 & 0 & 36 & 0 & 0 \end{pmatrix},$$

$$A_k = 0_{2 \times d_{2,k}}, \quad k \ge 7.$$
[7] Chen, Hank and Girelli, Florian, *Integrability from categorification*, Preprint, arXiv:2307.03831v1.

[8] F. Cummings, D. Dixon, and P. Kaus, [1992] *Dynamical model of the magnetic field of neutron stars*, Astrophys. J. 386, 215-221.

[9] A. A. M. Cuyt, *A review of multivariate Padé approximation theory*, Journal of Computational and Applied Mathematics Volumes **12–13**, May 1985, Pages 221-232.

[10] D. J. Daley and J. Gani, *Epidemic Modeling: An Introduction*, Cambridge University Press, Cambridge, 2005.

[11] D. D. Dixon, F. W. Cummings and P. E. Kaus, *Continuous "chaotic" dynamics in two dimensions*, Phys. D 65 (1993), no. 1-2, 109–116. MR1221249

[12] E. B. Dynkin, *Evaluation of the coefficients of the Campbell-Hausdorff formula*, Dokl. Akad. Nauk SSSR **57**, 323, 1947.

[13] P. Guillaume; A. Huard, *Multivariate Padé approximation*, Journal of Computational and Applied Mathematics **121**, Issues 1-2, 1 September 2000, Pages 197-219.

[14] T. Harko, F. S. N. Lobo, M. K. Mak, *Exact analytical solutions of the Susceptible-Infected-Recovered (SIR) epidemic model and of the SIR model with equal death and birth rates*, Applied Mathematics and Computation Volume 236, 1 June 2014, Pages 184-194.

[15] A. Helmi and H. Vucetich, *Non-integrability and chaos in classical cosmology*, Phys. Lett. A **230** (1997), no. 3-4, 153–156.

[16] Huang, Kaiyin; Shi, Shaoyun; Li, Wenlei. Meromorphic non-integrability of several 3D dynamical systems. Entropy 19 (2017), no. 5, Paper No. 211, 16 pp.

[17] W. O. Kermack and A. G. McKendrick, *A contribution to the mathematical theory of epidemics*, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, vol. 115, no. 772, pp. 700–721, 1997.

[18] M. Lakshmanan and R. Sahadevan, *Painlevé analysis, Lie symmetries, and integrability of coupled nonlinear oscillators of polynomial type*, Phys. Rep. 224 (1993), no. 1-2, 93.

[19] Li, W.L. and Shi, S.Y, Galoisian obstruction to the integrability of general dynamical systems. J. Differ. Equ. 2012, 252, 5518–5534.

[20] Li, W.L.; Shi, S.Y, Corrigendum to “Galoisian obstruction to the integrability of general dynamical systems”. J. Differ. Equ. 2017, 262, 1253–1256.

[21] K. Makino and M. Berz, *Suppression of the wrapping effect by Taylor model-based verified integrators: long-term stabilization by preconditioning*, Int. J. Differ. Equ. Appl. **10** (2005), no. 4, 353–384 (2006).

[22] R. Martínez and C. Simó, *Non-integrability of the degenerate cases of the swinging Atwood's machine using higher order variational equations*, Discrete Contin. Dyn. Syst. **29** (2011), no. 1, 1–24.

[23] R. Martínez and C. Simó. *Non-integrability of Hamiltonian systems through high order variational equations: summary of results and examples*, Regul. Chaotic Dyn. **14** (2009), no. 3, 323–348.
The fact that the ODE in (67) is so trivial to solve ($x(t) = -t + C$) adds no simplicity to what follows below. The fundamental matrix of $VE_1$ and the first filter matrix are equally simple:

$$Y_1 = \begin{pmatrix} 1 & 0 \\ 0 & e^{x(0)-x(t)} \end{pmatrix}, \quad F_1 = \begin{pmatrix} 0 & 0 \\ 0 & 1 \end{pmatrix}$$

and the generic form of the admissible solutions after applying Theorem 2.2.2 is the following: $f_k = (\frac{d}{dx} f_{k-1} | g_k(x))$ where $g_k$ is a solution to

$$-(k-1)k(x-1)x^3 g'_{k-2}(x) - g'_k(x) + (k-2)(k-1)k(x-1)x^2 g_{k-2}(x) + kg_k(x) = 0. \quad (68)$$

and moreover $g_k(x) = 0$ if $k$ is even, and $g_1(x) = e^x$. Equation (68) has a solution (which of course only merits being written if $k$ is odd)

$$g_k(x) = e^{kx} \int (k-1)k(\xi-1)\xi^2 e^{-k\xi} \left( (k-2)g_{k-2}(\xi) - \xi g'_{k-2}(\xi) \right) d\xi.$$

The formal first integral of (66) takes the form

$$F = \sum_{k=1}^{\infty} \frac{1}{k!} g_k(x) y^k = \sum_{i=0}^{\infty} (-1)^i e^{(2i+1)x} y^{2i+1} G_{2i+1};$$

the first terms of $G_\star$ are easy to discern,

$$G_1 = 1,$$

$$G_3 = \int_{x_0}^{x} e^{-2\xi}(\xi-1)^2 \xi^2 d\xi,$$

$$G_5 = \int_{x_0}^{x} (\xi - 1)^3 e^{-4\xi} \xi^5 d\xi + 3 \int_{\tau_0}^{\tau} (\xi - 1)^2 e^{-2\xi} \xi^2 \left( \int (\eta - 1)^2 e^{-2\eta} \eta^2 d\eta \right) d\xi$$

and the general term is

$$G_{2i+1} = \sum_{m=1}^{i} \sum_{j=1}^{\binom{i-1}{m-1}} D_{i-1, m-1, j} \Phi \left( (i-m) \mathbf{C}_{m, i-m-1, j} + \mathbf{1}_m, 2(i-m) \mathbf{C}_{m, i-m-1, j} + \mathbf{1}_m, 3(i-m) \mathbf{C}_{m, i-m-1, j} + \mathbf{1}_m \right)$$

where we define the following constant vectors: $\mathbf{1}_n := (1, \dots, 1) \in \mathbb{Z}^n$, and $C_{\star, \star, \star}$ are the columns of the symmetric product

$$\mathrm{Id}_k \odot \mathbf{1}_{d_{k,l}}^T = (C_{k,l,1} | C_{k,l,2} | \dots | C_{k,l,d_{k,l+1}}) \in \mathrm{Mat}_{k,d_{k,l+1}}(\mathbb{Z});$$

(obviously $C_{k,l,j} = 0$ for $l < 0$), $\Phi$ is defined recursively as follows: for every $\mathbf{a}, \mathbf{b}, \mathbf{c} \in \mathbb{Z}^m$,

$$\Phi(\mathbf{a}, \mathbf{b}, \mathbf{c}) := \begin{cases} \int_{x_0}^{x} (\tau - 1)^{a_1+1} e^{-(b_1+1)\tau} \tau^{c_1} d\tau, & m=1, \\ \int_{x_0}^{x} (\tau - 1)^{a_1+1} e^{-(b_1+1)\tau} \tau^{c_1} \Phi((a_2, \dots, a_m), (b_2, \dots, b_m), (c_2, \dots, c_m)) d\tau, & m>1 \end{cases}$$

and $D_{n,m,j}$ is the $j$th term in the vector $\mathbf{D}_{n,m}$ obtained from discarding terms containing repeating factors from the reverse-order version of $\mathbf{i}_n^T \odot \cdots \odot \mathbf{i}_n^T$ and $\mathbf{i}_n$ is the reverse-order sequence of the first $n$ odd integers $> 1$, e.g.

$$\mathbf{D}_{3,2} = (3 \cdot 5, 3 \cdot 7, 5 \cdot 7) = (D_{3,2,j})_{j=1,2,3,4}, \quad \mathbf{D}_{4,3} = (3 \cdot 5 \cdot 7, 3 \cdot 5 \cdot 9, 3 \cdot 7 \cdot 9, 5 \cdot 7 \cdot 9)$$
[24] A. Mitsopoulos and M. Tsamparlis, *Higher order first integrals of autonomous dynamical systems*, J. Geom. Phys. **170** (2021).

[25] A. Mitsopoulos, M. Tsamparlis and A. Ukpong, *Higher-Order First Integrals of Autonomous Non-Riemannian Dynamical Systems*, Symmetry **15**(1):222, DOI: 10.3390/sym15010222

[26] G. Montavon, S. Lapuschkin, A. Binder, W. Samek, K.-R. Müller, *Explaining nonlinear classification decisions with deep Taylor decomposition*, Pattern Recognition, Volume 65, May 2017, Pages 211-222.

[27] J. J. Morales-Ruiz, *Differential Galois theory and non-integrability of Hamiltonian systems*, Progress in Mathematics, Birkhäuser Verlag, Basel, 1999.

[28] J. J. Morales-Ruiz and J.-P. Ramis, *Galoisian obstructions to integrability of Hamiltonian systems*. I, Methods Appl. Anal. **8** (2001), no. 1, 33-96.

[29] J. J. Morales-Ruiz and J.-P. Ramis, *A note on the non-integrability of some Hamiltonian systems with a homogeneous potential*, Methods Appl. Anal. **8** (2001), no. 1, 113-120.

[30] J. J. Morales-Ruiz, J.-P. Ramis, and Carles Simó, *Integrability of Hamiltonian systems and differential Galois groups of higher variational equations*, Ann. Sci. École Norm. Sup. (4) **40** (2007), no. 6, 845-884.

[31] J. J. Morales-Ruiz and S. Simon, *On the meromorphic non-integrability of some N-body problems*. Discrete Contin. Dyn. Syst. **24** (2009), no. 4, 1225-1273.

[32] H. S. Rodrigues and M. J. Fonseca, *Can information be spread as a virus? Viral marketing as epidemiological model*, Mathematical Methods in the Applied Sciences, vol. 39, no. 16, pp. 4780-4786, 2016.

[33] W. M. Seiler, M. Seiß, *No chaos in Dixon's system*, Internat. J. Bifur. Chaos Appl. Sci. Engrgr. 31 (2021), no. 3, Paper No. 2150044, 20 pp.

[34] S. Simon, *Linearised higher variational equations*, Discrete Contin. Dyn. Syst. **34** (2014), no. 11, 4827-4854.

[35] S. Simon, *Conditions and evidence for non-integrability in the Friedmann-Robertson-Walker Hamiltonian*, J. Nonlinear Math. Phys. **21** (2014), no. 1, 1-16.

[36] M. van der Put and M. F. Singer, *Galois theory of linear differential equations*, Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences], vol. 328, Springer-Verlag, Berlin, 2003.

[37] Sprott, J. C., *Elegant chaos. Algebraically simple chaotic flows*, World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2010. xvi+285 pp. ISBN: 978-981-283-881-0; 981-283-881-3.

[38] S. L. Ziglin, *Bifurcation of solutions and the nonexistence of first integrals in Hamiltonian mechanics*. I, Funktsional. Anal. i Prilozhen. **16** (1982), no. 3, 30-41, 96.

[39] H. Żoładek, *The monodromy group*. Mathematics Institute of the Polish Academy of Sciences. Mathematical Monographs (New Series) **67**, Birkhäuser Verlag, Basel, 2006.
# 4 Future work and further comments

Proving that Theorem 2.3.1 is still true if we drop the local analyticity assumptions on $X$ and $\Upsilon^{-1}$ can be considered to be a minor-scope open problem although, as said above, Theorem 2.2.2 renders such task unnecessary when it comes to tackling particular dynamical systems, especially if computational ease is a factor and redundant terms arising from products and powers of formal first integrals are to be discarded.

Ever since it first took form in [3], the use of linearized higher variational equations to characterize, obstruct or redefine integrability is proving itself to be a welcome addition to the domain of dynamical systems. This contribution arrives at a time in which other attempts (e.g. [24, 25]) at finding methods to *integrate* dynamical systems rather than proving them non-integrable are gradually gaining traction.

A number of facts are germane to further studies:

(i) the potential shown by higher variational equations in finding adequate changes of variables that will simplify the system into tractability and, potentially, integrability (this is most exemplified in the sample systems tackled in §3);

(ii) linked to the above, the possible application of a *Baker-Campbell-Hausdorff* [12] sorts of formula upon variable transformations, using the decomposition of a transformed fundamental matrix into a product of matrices to our advantage;

(iii) the possible amenability of results in §2 above to machine learning [26];

(iv) the applicability of this automatic quadrature-algebraic method to the application of generalized, multivariate *Padé approximants* [9, 13] to first integrals, which will be tackled in an immediate future;

(v) the fact that the autonomous system (1) need not be Hamiltonian, thus allowing any arbitrary system to be considered and losing the serious constraint posed by symplectic transformations, since variable transformations can now be chosen freely;

(vi) also related to (v): the fact that this overrides the obstacles posed by the direct study of variational systems, which only has an implementation to Hamiltonian systems and requires first integrals to be meromorphic;

(vii) and finally, the potential (mentioned at the start of §2.3 and already opened in ??) for categorical functoriality translating the study of general dynamical systems to that of linear operators.

# References

[1] M. Abramowitz and I. A. Stegun (eds.), *Handbook of mathematical functions with formulas, graphs, and mathematical tables*, A Wiley-Interscience Publication, John Wiley & Sons Inc., New York, 1984, Reprint of the 1972 edition, Selected Government Publications.

[2] J. Alvarez-Ramirez, J. Delgado-Fernandez and G. Espinosa-Paredes *The origin of a continuous two-dimensional "chaotic" dynamics*, Int. J. Bif. Chaos 15, 3023-3029, 2005.

[3] A. Aparicio-Monforte, M. Barkatou, S. Simon, J.-A. Weil, *Formal first integrals along solutions of differential systems I*, ISSAC 2011 Proceedings of the 36th International Symposium on Symbolic and Algebraic Computation, 19-26, ACM, New York, 2011.
[4] U. Bekbaev, *A matrix representation of composition of polynomial maps*, arXiv:0901.3179v3 [math.AC] 22 Sep 2009.

[5] U. Bekbaev, *Matrix representations for symmetric and antisymmetric multi-linear maps*, arXiv:1010.2579v1 [math.RA] 13 Oct 2010

[6] Chen, Ding, *On the Integrability of the SIR Epidemic Model with Vital Dynamics*, Hindawi Advances in Mathematical Physics, Volume 2020, Article ID 5869275, https://doi.org/10.1155/2020/5869275

[7] Chen, Hank and Girelli, Florian, *Integrability from categorification*, Preprint, arXiv:2307.03831v1.

[8] F. Cummings, D. Dixon, and P. Kaus, [1992] *Dynamical model of the magnetic field of neutron stars*, Astrophys. J. 386, 215-221.

[9] A. A. M. Cuyt, *A review of multivariate Padé approximation theory*, Journal of Computational and Applied Mathematics Volumes **12–13**, May 1985, Pages 221-232.

[10] D. J. Daley and J. Gani, *Epidemic Modeling: An Introduction*, Cambridge University Press, Cambridge, 2005.

[11] D. D. Dixon, F. W. Cummings and P. E. Kaus, *Continuous "chaotic" dynamics in two dimensions*, Phys. D 65 (1993), no. 1-2, 109–116. MR1221249

[12] E. B. Dynkin, *Evaluation of the coefficients of the Campbell-Hausdorff formula*, Dokl. Akad. Nauk SSSR **57**, 323, 1947.

[13] P. Guillaume; A. Huard, *Multivariate Padé approximation*, Journal of Computational and Applied Mathematics **121**, Issues 1-2, 1 September 2000, Pages 197-219.

[14] T. Harko, F. S. N. Lobo, M. K. Mak, *Exact analytical solutions of the Susceptible-Infected-Recovered (SIR) epidemic model and of the SIR model with equal death and birth rates*, Applied Mathematics and Computation Volume 236, 1 June 2014, Pages 184-194.

[15] A. Helmi and H. Vucetich, *Non-integrability and chaos in classical cosmology*, Phys. Lett. A **230** (1997), no. 3-4, 153–156.

[16] Huang, Kaiyin; Shi, Shaoyun; Li, Wenlei. Meromorphic non-integrability of several 3D dynamical systems. Entropy 19 (2017), no. 5, Paper No. 211, 16 pp.

[17] W. O. Kermack and A. G. McKendrick, *A contribution to the mathematical theory of epidemics*, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, vol. 115, no. 772, pp. 700-721, 1997.

[18] M. Lakshmanan and R. Sahadevan, *Painlevé analysis, Lie symmetries, and integrability of coupled nonlinear oscillators of polynomial type*, Phys. Rep. 224 (1993), no. 1-2, 93.

[19] Li, W.L. and Shi, S.Y, Galoisian obstruction to the integrability of general dynamical systems. J. Differ. Equ. 2012, 252, 5518–5534.

[20] Li, W.L.; Shi, S.Y, Corrigendum to “Galoisian obstruction to the integrability of general dynamical systems”. J. Differ. Equ. 2017, 262, 1253–1256.

[21] K. Makino and M. Berz, *Suppression of the wrapping effect by Taylor model-based verified integrators: long-term stabilization by preconditioning*, Int. J. Differ. Equ. Appl. **10** (2005), no. 4, 353–384 (2006).
[22] R. Martínez and C. Simó, *Non-integrability of the degenerate cases of the swinging Atwood's machine using higher order variational equations*, Discrete Contin. Dyn. Syst. **29** (2011), no. 1, 1-24.

[23] R. Martínez and C. Simó, *Non-integrability of Hamiltonian systems through high order variational equations: summary of results and examples*, Regul. Chaotic Dyn. **14** (2009), no. 3, 323-348.

[24] A. Mitsopoulos and M. Tsamparlis, *Higher order first integrals of autonomous dynamical systems*, J. Geom. Phys. **170** (2021).

[25] A. Mitsopoulos, M. Tsamparlis and A. Ukpong, *Higher-Order First Integrals of Autonomous Non-Riemannian Dynamical Systems*, Symmetry **15**(1):222, DOI: 10.3390/sym15010222

[26] G. Montavon, S. Lapuschkin, A. Binder, W. Samek, K.-R. Müller, *Explaining nonlinear classification decisions with deep Taylor decomposition*, Pattern Recognition, Volume 65, May 2017, Pages 211-222.

[27] J. J. Morales-Ruiz, *Differential Galois theory and non-integrability of Hamiltonian systems*, Progress in Mathematics, Birkhäuser Verlag, Basel, 1999.

[28] J. J. Morales-Ruiz and J.-P. Ramis, *Galoisian obstructions to integrability of Hamiltonian systems. I*, Methods Appl. Anal. **8** (2001), no. 1, 33-96.

[29] J. J. Morales-Ruiz and J.-P. Ramis, *A note on the non-integrability of some Hamiltonian systems with a homogeneous potential*, Methods Appl. Anal. **8** (2001), no. 1, 113-120.

[30] J. J. Morales-Ruiz, J.-P. Ramis, and Carles Simó, *Integrability of Hamiltonian systems and differential Galois groups of higher variational equations*, Ann. Sci. École Norm. Sup. (4) **40** (2007), no. 6, 845-884.

[31] J. J. Morales-Ruiz and S. Simon, *On the meromorphic non-integrability of some N-body problems*. Discrete Contin. Dyn. Syst. **24** (2009), no. 4, 1225-1273.

[32] H. S. Rodrigues and M. J. Fonseca, *Can information be spread as a virus? Viral marketing as epidemiological model*, Mathematical Methods in the Applied Sciences, vol. 39, no. 16, pp. 4780-4786, 2016.

[33] W. M. Seiler, M. Seiß, *No chaos in Dixon's system*, Internat. J. Bifur. Chaos Appl. Sci. Engrg. 31 (2021), no. 3, Paper No. 2150044, 20 pp.

[34] S. Simon, *Linearised higher variational equations*, Discrete Contin. Dyn. Syst. **34** (2014), no. 11, 4827-4854.

[35] S. Simon, *Conditions and evidence for non-integrability in the Friedmann-Robertson-Walker Hamiltonian*, J. Nonlinear Math. Phys. **21** (2014), no. 1, 1-16.

[36] S. Simon, *Tools for symmetric products of finite and infinite matrices*, GitHub repository https://github.com/Sergi-Simon/Tools-for-symmetric-products-of-finite-and-infinite-matrices

[37] M. van der Put and M. F. Singer, *Galois theory of linear differential equations*, Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences], vol. 328, Springer-Verlag, Berlin, 2003.
[38] Sprott, J. C., *Elegant chaos. Algebraically simple chaotic flows*, World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2010. xvi+285 pp. ISBN: 978-981-283-881-0; 981-283-881-3.

[39] S. L. Ziglin, *Bifurcation of solutions and the nonexistence of first integrals in Hamiltonian mechanics*. I, Funktsional. Anal. i Prilozhen. **16** (1982), no. 3, 30–41, 96.

[40] H. Żoładek, *The monodromy group*. Mathematics Institute of the Polish Academy of Sciences. Mathematical Monographs (New Series) **67**, Birkhäuser Verlag, Basel, 2006.