# Systems of Discrete Differential Equations, Constructive Algebraicity of the Solutions

Hadrien Notarantonio, Sergey Yurkevich

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

## To cite this version:

Hadrien Notarantonio, Sergey Yurkevich. Systems of Discrete Differential Equations, Constructive Algebraicity of the Solutions. 2023. hal-04250469

HAL Id: hal-04250469

https://inria.hal.science/hal-04250469v1

Preprint submitted on 19 Oct 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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Systems of Discrete Differential Equations, Constructive Algebraicity of the Solutions

Hadrien Notarantonio*, Sergey Yurkevich†

## Abstract

In this article, we study systems of $n \ge 1$, not necessarily linear, discrete differential equations (DDEs) of order $k \ge 1$ with one catalytic variable. We provide a constructive and elementary proof of algebraicity of the solutions of such equations. This part of the present article can be seen as a generalization of the pioneering work by Bousquet-Mélou and Jehanne (2006) who settled down the case $n = 1$. Moreover, we obtain effective bounds for the algebraicity degrees of the solutions and provide an algorithm for computing annihilating polynomials of the algebraic series. Finally, we carry out a first analysis in the direction of effectivity for solving systems of DDEs in view of practical applications.

# 1 Introduction

## 1.1 Context and motivation

The equations that lie in the interest of this work are so-called discrete differential equations with one catalytic variable of fixed-point type. They take the form

$$F(t, u) = f(u) + t \cdot Q(F(t, u), \Delta_a F(t, u), \dots, \Delta_a^k F(t, u), t, u), \quad (1)$$

where $k \in \mathbb{N}$ is called the order of the DDE, $f$ and $Q$ are polynomials, and (for some $a \in \mathbb{Q}$) $\Delta_a^\ell$ is the $\ell$th iteration of the discrete derivative operator $\Delta_a : \mathbb{Q}[u][[t]] \to \mathbb{Q}[u][[t]]$ defined by

$$\Delta_a F(t, u) := \frac{F(t, u) - F(t, a)}{u - a}.$$

Discrete differential equations are ubiquitous in enumerative combinatorics [25, 12, 13, 10, 5]. Indeed, enumerating discrete structures usually leads to introducing the corresponding generating function, say $G(t)$. In some of such cases, the combinatorial nature of the initial problem can be transformed directly into an algebraic or analytic question on the equation satisfied by $G(t)$. However for many practical counting problems, the initial combinatorial structure is too coarse to be translated into any meaningful equation. In these cases it is often helpful to try to solve a more refined problem, introducing new structure into the initial question and consequently a new variable $u$ in the generating function: one is then lead to study a DDE satisfied by some

<footnote>*Inria, Palaiseau (France). hadrien.notarantonio@inria.fr</footnote>

<footnote>†University of Vienna (Austria). sergey.yurkevich@univie.ac.at</footnote>
bivariate series $F(t, u) \in \mathbb{Q}[u][[t]]$, even though the interest only lies in the specialization, usually $G(t) = F(t, 0)$ or $G(t) = F(t, 1)$. Although this idea is, of course, very classical and used since decades [25, 12, 13], the name “catalytic” for such a variable $u$ was introduced only relatively recently by Zeilberger [27] in the year 2000.

**Example 1.** The so-called 2-constellations are special bi-colored planar maps (see [10, §5.3] for the definition). Let the sequence $(a_n)_{n \ge 0}$ enumerate the 2-constellations with $n$ black faces; we wish to discover properties of $G(t) = \sum_{n \ge 0} a_n t^n$ (e.g. a closed-form expression for $G(t)$ or the numbers $a_n$, asymptotics of $(a_n)_{n \ge 0}$, etc). This problem is usually refined by considering the numbers $a_{n,d}$ enumerating 2-constellations with $n$ black faces and outer degree 2d. With now more constraints on the studied enumeration (black faces and outer-degree), it is possible to show by a recursive analysis of the construction of a 2-constellation (see [10, Section 5.3] for details) that the bivariate generating function $F(t, u) := \sum_{n,d \ge 0} a_{n,d} u^d t^n \in \mathbb{Q}[u][[t]]$ satisfies the DDE of order 1

$$F(t, u) = 1 + t u F(t, u)^2 + t u \Delta_1 F(t, u). \qquad (2)$$

Note that $G(t) = F(t, 1) \in \mathbb{Q}[[t]]$ is the generating function we were initially interested in. The classical way of studying the series $G(t)$ is by considering Equation (2). For instance, by studying (2), it can be shown using [14], [10, Section 2] or [8, Proposition 2.4] that $F(t, u)$ is an algebraic function over $\mathbb{Q}(t, u)$ and consequently that $G(t)$ is algebraic over $\mathbb{Q}(t)$. Explicitly,

$$16t^3 G(t)^2 - (8t^2 + 12t - 1)tG(t) + t(t^2 + 11t - 1) = 0.$$

Knowing that $R(t, z) := 16t^3z^2 - (8t^2 + 12t - 1)tz + t(t^2 + 11t - 1) \in \mathbb{Q}[t, z]$ annihilates $G(t)$, it is possible to show that for $n \ge 1$, one has the closed-form expression

$$a_n = 3 \frac{2^{n-1}}{(n+2)(n+1)} \binom{2n}{n}$$

whose asymptotic behavior is $a_n \sim 3 \cdot 8^n / \sqrt{4\pi n^5}$, by Stirling's formula.

The algebraicity of $F(t, u)$ over $\mathbb{Q}(t, u)$ in example 1 is no coincidence. In their pioneering work [10], Bousquet-Mélou and Jehanne proved (see [10, Theorem 3]) that the unique power series solution of a functional equation of the form (1) is always an algebraic function.

In a variety of different contexts throughout combinatorics also appears their natural extension systems of DDEs of fixed-point type that are, systems of the form

$$\left\{ \begin{array}{l} F_1 = f_1(u) + t \cdot Q_1(\nabla_a^k F_1, \dots, \nabla_a^k F_n, t, u), \\ \vdots \\ F_n = f_n(u) + t \cdot Q_n(\nabla_a^k F_1, \dots, \nabla_a^k F_n, t, u), \end{array} \right. \qquad (3)$$

where for $i = 1, \dots, n$ the polynomials $f_i \in \mathbb{Q}[u]$, $Q_i \in \mathbb{Q}[y_1, \dots, y_{n(k+1)}, t, u]$ are given and we write $\nabla_a^k F := (F, \Delta_a F, \dots, \Delta_a^k F)$ for $\Delta_a$ as before, and the unknowns $F_i = F_i(t, u) \in \mathbb{Q}[u][[t]]$. As previously, $k$ is called the order of (3). Such systems of functional equations appear, for instance, in the enumeration of hard particles on planar maps [10, §5.4], inhomogeneous lattice paths [15], certain orientations with $n$ edges [7, §5], or parallelogram polyominoes [4, §7.1].

For any given functional equation or system of such, the results one typically wishes to obtain are for instance: a closed-form expression for the number of objects of a given size, a grasp on the
asymptotic behavior, or a classification of the nature of the generating function(s) (e.g. algebraic, D-finite, etc.). It is often easier to obtain a closed-formula or an asymptotic estimate by studying a witness of the nature of the generating function (e.g. an annihilating polynomial, a linear ODE, etc.), so a natural aim is to solve the system, i.e. to compute such a witness. These objectives frequently yield arbitrarily difficult challenges at the intersection with enumerative combinatorics, theoretical physics [11] and computational geometry [16].

**Main goals** We want to prove in an elementary way that the components of the unique solution $(F_1, \dots, F_n) \in \mathbb{Q}[u][[t]]^n$ of (3) are always algebraic power series over $\mathbb{Q}(t, u)$. Moreover in the setting $n > 1$ in (3), we want to design, analyze and theoretically compare geometry-driven algorithms that compute an annihilating polynomial of the specialized series $F_1(t, a)$.

Before providing a state-of-the-art and stating our contributions, we introduce a combinatorial example of systems of DDEs that we shall use intensively for illustrating this paper.

**Example 2.** The following system of DDEs for the generating function of certain planar orientations was considered in [7, Eq. (27)]:

$$ \begin{cases} (\mathbf{E}_{F_1}): F_1(t, u) = 1 + t \cdot \left( u + 2uF_1(t, u)^2 + 2uF_2(t, 1) + u \frac{F_1(t, u) - uF_1(t, 1)}{u-1} \right), \\ (\mathbf{E}_{F_2}): F_2(t, u) = t \cdot \left( 2uF_1(t, u)F_2(t, u) + uF_1(t, u) + uF_2(t, 1) + u \frac{F_2(t, u) - uF_2(t, 1)}{u-1} \right). \end{cases} \quad (4) $$

We show in section 2 that $F_1(t, 1) = 1 + 2t + 10t^2 + 66t^3 + \dots$ is algebraic over $\mathbb{Q}(t)$, and that its minimal polynomial $64t^3z_0^3 + (48t^3 - 72t^2 + 2t)z_0^2 - (15t^3 - 9t^2 - 19t + 1)z_0 + t^3 + 27t^2 - 19t + 1$ can be computed using tools coming from elimination theory.

## 1.2 Previous works

In the lines below, we present an overview of the main results regarding algebraicity of solutions of systems of DDEs. We start with the study of a single DDE for $k > 1$ (the case $k = 1$ is not of our interest since, as we shall see later in this article, systems of DDEs of order $k \ge 1$ reduce to a single functional equation that "involves" at least 2 univariate series).

In the seminal work [10], Bousquet-Mélou and Jehanne completely resolve algebraicity of solutions of DDEs (1), equivalently the case $n = 1$ in (3). Moreover, Bousquet-Mélou and Jehanne provide systematic methods for computing an annihilating polynomial of the specialized series $F(t, a)$. For proving [10, Theorem 3], the authors designed a “nonlinear kernel method” which allows one to prove that the unique solution of (1) is always an algebraic function over $\mathbb{Q}(t, u)$. Significantly in practice, this approach yields an algorithm for finding an annihilating polynomial of the specialization $F(t, a)$ and of the bivariate series $F(t, u)$. The idea of their algebraicity proof is to reduce the resolution of the DDE to solving some system of polynomial equations which has a solution whose coordinates contains the involved specializations of $F$. Their proof involves a symbolic deformation argument ensuring that the polynomial system which is constructed contains enough independent equations. For efficiency considerations for the resolution of a single DDE of the form (1), a recent algorithmic work [9] by Bostan, Safey El Din and Notarantonio targeted the intensive use of effective algebraic geometry in order to efficiently solve the underlying polynomial systems.

Regarding systems of DDEs (i.e. the case $n > 1$), the usual strategy (e.g. [10, Section 11])
is to reduce a given system to a single equation and to apply the method from [10, Section 2]. Nevertheless, since the reduced equation may not be of the form (1) anymore, the ideas of [10, Section 2] may not be applicable. In the literature, there exist two methods to overcome these theoretical issues. First, a deep theorem in commutative algebra by Popescu [23], central in the so-called “nested Artin approximation” theory, guarantees that equations of the form (5) always admit an algebraic solution (see also [7, Theorem 16] for a statement of this theorem). Note that the nested condition is automatically satisfied in this case and that the uniqueness of the solution is obvious. A drawback of using Popescu’s theorem, however, is that its proof is a priori highly non-constructive and can only be applied as a “black box”, whereas in practice one is often interested in the explicit annihilating polynomials of the solutions. Secondly, the frequent case when (5) is linear in the bivariate formal power series and their specializations was effectively solved (i.e. their proof of algebraicity yields an algorithm) in the more recent article [15] by Buchacher and Kauers by using a multi-dimensional kernel method. However even if their proof yields an algorithm, its efficiency was not discussed at all. Note that the now common multi-dimensional kernel method appears as well in the article of the same year [1] by Asinowski, Bacher, Banderier and Gittenberger.

Before the present work, however, there was no systematic approach for dealing with systems of DDEs such as (4). It is in this context that the present paper takes place.

## 1.3 Contributions

This paper is the full version of the extended abstract [22] that was published in the proceedings of the conference Formal Power Series and Algebraic Combinatorics 2023.

The first contribution of this article is a generalization of the algebraicity result [10, Theorem 3] to systems of discrete differential equations of a fixed-point type. Precisely, we prove the following theorem. Here and in the following, we denote $K$ a field of characteristic 0.

**Theorem 3.** Let $n, k \ge 1$ be integers and $f_1, \dots, f_n \in K[u], Q_1, \dots, Q_n \in K[y_1, \dots, y_{n(k+1)}, t, u]$ be polynomials. For $a \in K$, set $\nabla_a^k F := (F, \Delta_a F, \dots, \Delta_a^k F)$. Then the system of equations

$$ \left\{ \begin{array}{l} (\mathbf{E}_{F_1}): \quad F_1 = f_1(u) + t \cdot Q_1(\nabla_a^k F_1, \dots, \nabla_a^k F_n, t, u), \\ \vdots \\ (\mathbf{E}_{F_n}): \quad F_n = f_n(u) + t \cdot Q_n(\nabla_a^k F_1, \dots, \nabla_a^k F_n, t, u) \end{array} \right. \qquad (5) $$

admits a unique vector of solutions $(F_1, \dots, F_n) \in K[u][[t]]^n$, and all its components are algebraic functions over $K(t, u)$.

The key idea, analogous to the one in the proof of [10, Theorem 3], for proving this theorem is to define a deformation of (5) that ensures the applicability of a multi-dimensional analog of the “nonlinear kernel method”. Stated explicitly, we show in lemma 6 that after deforming the equations as in (7), the polynomial in $u$ defined by the determinant of the Jacobian matrix associated to the equations in (5) (considered with respect to the $F_i$) has exactly $nk$ solutions in an extension of the ring $\bigcup_{d \ge 1} \overline{K}[[t^{1/d}]]$. After a process of “duplication of variables”, we construct a zero-dimensional polynomial ideal, a non-trivial element of which must be the desired annihilating polynomial. The most technical step consists in proving the invertibility of a certain Jacobian matrix (lemma 8 and lemma 9) in order to justify the zero-dimensionality.
The second contribution is the analysis of the resulting algorithm for finding annihilating
polynomials of the power series $F_i(t, u)$ in theorem 3. From our constructive proof we deduce
a theoretical upper bound on the algebraicity degree of each $F_i$. Moreover using the radicality
of the constructed 0-dimensional ideal, and when the field $\mathbb{K}$ is effective (e.g. $\mathbb{K} = \mathbb{Q}$), we also
bound the arithmetic complexity of our algorithm, that is the number of operations $(+, -, \times, \div)$
performed in $\mathbb{K}$. Denoting by $\mathrm{totdeg}(P)$ the total degree of a multivariate polynomial $P$, we obtain
the following:

**Theorem 4.** In the setting of theorem 3, let $(F_1, \dots, F_n) \in \mathbb{K}[u][[t]]^n$ be the vector of solutions and $\delta := \max(\deg(f_1), \dots, \deg(f_n), \mathrm{totdeg}(Q_1), \dots, \mathrm{totdeg}(Q_n))$. Then the algebraicity degree of each $F_i(t, u)$ over $\mathbb{K}(t, u)$ is bounded by $n^{2n^2k^2}(k+1)^{n^2k^2(n+2)+n}\delta^{n^2k^2(n+2)+n}/(nk)!^{nk}$. Moreover if $\mathbb{K}$ is effective, there exists an algorithm computing an annihilating polynomial of any $F_i(t, a)$ in $O((nk\delta)^{40(n^2k+1)})$ arithmetic operations in $\mathbb{K}$.

Let us emphasize that despite the desperately looking exponent 40 in theorem 4, one can still
solve concrete examples from time to time as we shall see in section 2 and section 4.

The third contribution is the full algorithmic investigation of two natural schemes for solving
systems of DDEs. For each of them, we analyze the conditions under which they might be applied,
and the possible links between their respective outputs. The first algorithm consists in the classical
duplication of variables argument, by following our proof of theorem 3, and then performing a
brute force elimination of all irrelevant variables (lemma 11). The second scheme consists in
reducing the initial system of DDEs to a single polynomial functional equation where the general
method of [10, Section 2] and the recent algorithmic improvements made in [9, Section 5] might
apply. In this direction, we identify in section 4.2 sufficient conditions under which [10, Section 2]
and [9, Section 5] can systematically be used. At the end of section 4.2, we show that eq. (4) can
not be solved by the state-of-the-art based on reducing a system of DDEs to a single equation (and
then in applying any systematic method), while it can be solved by the systematic method that we
introduce in section 2.

**Structure of the paper.** In section 2, we explain our method in the case of two equations of order one under the assumption that no deformation is necessary. We summarize the method in an algorithm and showcase it explicitly on example 2. In section 3, we provide proofs of theorem 3 and theorem 4 with more details than in the extended abstract [22]. In section 4, we study and compare the output of two natural strategies for solving systems of DDEs. Ultimately, we discuss some necessary future works in section 5.

**Notations.** Throughout this article, $\mathbb{K}$ denotes a field of characteristic 0, $\overline{\mathbb{K}}$ its algebraic closure, $\mathbb{K}[[t]]$ the ring of formal power series in $t$ with coefficients in $\mathbb{K}$, and $\overline{\mathbb{K}}[[t^{\frac{1}{\delta}}]]$ the ring $\bigcup_{d \ge 1} \mathbb{K}[[t^{\frac{1}{d}}]]$ of Puiseux series with rational positive exponent in the variable $t$. Also for $n, k \ge 1$ and $a \in \mathbb{K}$, we use the compact notation $A(u)$ for any given polynomial expression of the form $A(F_1(t, u), F_1(t, a), \dots, \partial_u^{k-1} F_1(t, a), \dots, F_n(t, u), F_n(t, a), \dots, \partial_u^{k-1} F_n(t, a), t, u)$. For an integer $N > 1$, we denote by $\mathbb{K}[x_1, \dots, x_N]$ the ring of polynomials in the variables $x_1, \dots, x_N$ with coefficients in $\mathbb{K}$. For $P \in \mathbb{K}[x_1, \dots, x_N]$, we denote by $\partial_{x_i} P$ the partial derivative of $P$ with respect to the variable $x_i$ (for $i \in \{1, \dots, N\}$) and by $V(P)$ its zero set in $\overline{\mathbb{K}}^N$. For $\mathcal{I}$ an ideal of $\mathbb{K}[x_1, \dots, x_N]$, we also denote by $V(\mathcal{I})$ the zero set of $\mathcal{I}$ in $\overline{\mathbb{K}}^N$. Also, for $S$ a set of polynomials
in $K[x_1, \dots, x_N]$, we denote by $V(S)$ the zero set of the ideal generated by $S$ in $K[x_1, \dots, x_N]$. For $x$, $y$ two sets of variables, we denote $\{x\} >_{lex} \{y\}$ the monomial order such that $\{x\}$ (resp. $\{y\}$) is the usual degrevlex [18, Definition 5, § 2, Chapter 2] order over $x$ (resp. over $y$), and such that any monomial in $y$ is lower than any monomial containing at least one variable in $x$ (by default $\{x\}$ will always denote the usual degrevlex monomial order on the variables $x$). We use the soft-O notation $\tilde{O}(\cdot)$ for hiding polylogarithmic factors in the argument.

## 2 General strategy and application to a first example

Before proving our main theorem in section 3, we introduce our general method in the situation of two equations of order 1. We illustrate each step with the system (4) from example 2.

Starting with (5) for $n = 2$ and $k = 1$, we first multiply $(E_{F_1})$ and $(E_{F_2})$ by $(u-a)^{m_1}$ and $(u-a)^{m_2}$ respectively (for $m_1, m_2 \in \mathbb{N}$) in order to obtain a system with polynomial coefficients in $u$. By a slight abuse of notation, we shall still write $(E_{F_1})$ and $(E_{F_2})$ for those equations. Note that this system induces polynomials $E_1, E_2$ in $K(t)[x_1, x_2, z_0, z_1, u]$ whose specializations to $x_1 = F_1(t, u), x_2 = F_2(t, u), z_0 = F_1(t, a), z_1 = F_2(t, a)$, denoted by $E_1(u), E_2(u)$, are zero.

**Example 1 (cont.).** Multiplying $(E_{F_1})$ and $(E_{F_2})$ in example 2 by $u-1$ gives

$$ \begin{cases} E_1 = (1 - x_1) \cdot (u - 1) + t \cdot (2u^2 x_1^2 - u^2 z_0 + 2u^2 z_1 - 2ux_1^2 + u^2 + ux_1 - 2uz_1 - u), \\ E_2 = x_2 \cdot (1 - u) + t \cdot (2u^2 x_1 x_2 + u^2 x_1 - 2ux_1 x_2 - ux_1 + ux_2 - uz_1). \end{cases} $$

Note that applying the specializations $x_1 = F_1(t, u), x_2 = F_2(t, u), z_0 = F_1(t, 1), z_1 = F_2(t, 1)$ to $E_1$ and $E_2$ yields the vanishing of the induced polynomial functional equations.

In the spirit of [10], we take the derivative of both equations with respect to the variable $u$:

$$ \left( \begin{pmatrix} \partial_{x_1} E_1 \\ \partial_{x_1} E_2 \end{pmatrix}(u) \quad \begin{pmatrix} \partial_{x_2} E_1 \\ \partial_{x_2} E_2 \end{pmatrix}(u) \right) \cdot \begin{pmatrix} \partial_u F_1 \\ \partial_u F_2 \end{pmatrix} + \begin{pmatrix} (\partial_u E_1)(u) \\ (\partial_u E_2)(u) \end{pmatrix} = 0, \quad (6) $$

Define $\text{Det} := \partial_{x_1} E_1 \cdot \partial_{x_2} E_2 - \partial_{x_1} E_2 \cdot \partial_{x_2} E_1 \in K(t)[x_1, x_2, z_0, z_1, u]$. One can show that the specialization $\text{Det}(F_1(t, u), F_2(t, u), F_1(t, a), F_2(t, a), t, u) \in K[[t]][[u]]$ admits either 0, 1 or 2 distinct non-constant solutions in $u$ in $\overline{K}[[t^\star]]$. We assume that there exist 2 such solutions $U_1, U_2 \in \overline{K}[[t^\star]]$; we prove in section 3 that it is always the case up to the deformation (7).

Exploiting the common idea to [1, Proof of Theorem 3.2] and [15, Section 3], we define the vector $v := (\partial_{x_1} E_2, -\partial_{x_1} E_1) \in K(t)[x_1, x_2, z_0, z_1, u]^2$ and plug $U_1$ for $u$ into $v$ and (6). Note that $v$ is an element of the left-kernel of the square matrix in (6) mod $\text{Det}(x_1, x_2, z_0, z_1, u)$. After multiplication of (6) by $v$ on the left, we find a new polynomial relation relating the series $F_1(t, U_i), F_2(t, U_i), F_1(t, a), F_2(t, a), t$ and $U_i$, namely $\partial_{x_1} E_1 \cdot \partial_u E_2 - \partial_{x_1} E_2 \cdot \partial_u E_1 = 0$ when evaluated at $x_1 = F_1(t, U_i), x_2 = F_2(t, U_i), z_0 = F_1(t, a), z_1 = F_2(t, a), u = U_i$. We denote $P := \partial_{x_1} E_1 \cdot \partial_u E_2 - \partial_{x_1} E_2 \cdot \partial_u E_1 \in K(t)[x_1, x_2, z_0, z_1, u]$ this new polynomial.

**Remark 5.** Note that $P$ is the determinant of the matrix $\begin{pmatrix} \partial_{x_1} E_1 & \partial_{x_1} E_2 \\ \partial_u E_1 & \partial_u E_2 \end{pmatrix}$, which is not a coincidence, as we will see in the next section.

We define the polynomial system $S := (E_1, E_2, \text{Det}, P) \in K(t)[x_1, x_2, z_0, z_1, u]^4$. It admits the relevant solutions $(F_1(t, U_i), F_2(t, U_i), F_1(t, a), F_2(t, a), U_i) \in \overline{K}[[t^\star]]^5$, for $i \in \{1, 2\}$.
**Example 1 (cont.).** Continuing example 2, we find

$$ egin{cases} \operatorname{Det} = (4tu^2x_1 - 4tux_1 + tu - u + 1)(2tu^2x_1 - 2tux_1 + tu - u + 1), \\ P = -2tx_1x_2 - tx_1 + tx_2 - tz_1 - x_2 + P_1 \cdot u + P_2 \cdot u^2 + P_3 \cdot u^3, \end{cases} $$

where $P_1, P_2, P_3$ are explicit (but relatively big) polynomials in $\mathbb{Q}[x_1, x_2, z_0, z_1, t]$.

Applying in spirit the steps of [10, Section 2], we define for $i \in \{0, 1\}$ the polynomial systems $S_i := S(x_{2i+1}, x_{2i+2}, z_0, z_1, u_{i+1})$ by “duplicating” variables. In the case where the ideal $S_{\text{dup}} := \langle S_0, S_1, m \cdot (u_1 - u_2) - 1 \rangle$ has dimension 0 over $\mathbb{K}(t)$, finding an annihilating polynomial of $F_1(t, a)$ is done by computing a nonzero element of $\langle S_0, S_1, m \cdot (u_1 - u_2) - 1 \rangle \cap \mathbb{K}[z_0, t]$.

**Example 1 (cont.).** Continuing Example 2, we compute<fnref n="1" /> a generator of the polynomial ideal $\langle S_0, S_1, m \cdot (u_1 - u_2) - 1 \rangle \cap \mathbb{Q}[z_0, t]$. It has degree 13 in $z_0$ and 14 in $t$. In particular, it contains in its factors the minimal polynomial of $F_1(t, 1)$ given by $64t^3z_0^3 + (48t^3 - 72t^2 + 2t)z_0^2 - (15t^3 - 9t^2 - 19t + 1)z_0 + t^3 + 27t^2 - 19t + 1$.

We summarize the algorithm described above in the compact form given by Algorithm 1.

**Algorithm 1:** Solving systems of two discrete differential equations of order 1.

**Input:** A system of two DDEs $(E_{F_1}), (E_{F_2})$ of order 1, with $S_{\text{dup}}$ of dimension 0 over $\mathbb{K}(t)$.

**Output:** A nonzero $R \in \mathbb{K}[z_0, t]$ annihilating $F_1(t, a)$.

1. Replace $(E_{F_1})$ and $(E_{F_2})$ by their respective numerators and denote by $E_1$ and $E_2$ the associated polynomials in $\mathbb{K}(t)[x_1, x_2, z_0, z_1, u]$.

2. Compute $\operatorname{Det} := \partial_{x_1}E_1 \cdot \partial_{x_2}E_2 - \partial_{x_1}E_2 \cdot \partial_{x_2}E_1$ and $P := \partial_{x_1}E_1 \cdot \partial_{u}E_2 - \partial_{x_1}E_2 \cdot \partial_{u}E_1$.

3. Set $S := (E_1, E_2, \operatorname{Det}, P) \subset \mathbb{K}(t)[x_1, x_2, z_0, z_1, u]$.

4. For $0 \le i \le 1$, define $S_i := S(x_{2i+1}, x_{2i+2}, z_0, z_1, u_{i+1})$.

5. Define $S_{\text{dup}} := \langle S_0, S_1, m \cdot (u_1 - u_2) - 1 \rangle \subset \mathbb{K}(t)[m, x_1, x_2, x_3, x_4, z_0, z_1, u_1, u_2]$,

6. Return a nonzero element of $S_{\text{dup}} \cap \mathbb{K}[z_0, t]$.

We remark that if the same strategy as above is applied in the case of a single equation of first order of the form $F_1 = f(u) + t \cdot Q_1(F_1, \Delta_a F_1, t, u)$, the presented method simplifies to the classical method in [10] of Bousquet-Mélou and Jehanne which relies on studying the ideal $\langle E_1, \partial_{x_1}E_1, \partial_{u}E_1 \rangle$. Stated explicitly, $\partial_{x_1}E_1$ plays the role of $\operatorname{Det}$ and $\partial_{u}E_1$ plays the role of $P$.

# 3 Proofs of theorem 3 and theorem 4

## 3.1 Proof of theorem 3

As explained before, the statement and proof can be seen as a generalization of [10, Theorem 3] and [15, Theorem 2], so several steps are done analogously. Without loss of generality we assume that $a = 0$ and set $\Delta := \Delta_0$ and $\nabla := \nabla_0$.

Denote by $m_1, \dots, m_n$ the least positive integers greater than or equal to $k$ such that multiplying $(E_{F_i})$ in (5) by $u^{m_i}$ gives a polynomial equation; in other words, the multiplication by $u^{m_i}$ clears the denominators introduced by the application of $\Delta$. Set $\beta := [2M/k]$ and $\alpha := 3n^2k \cdot (\beta+1) + 3nM$,

<footnote n="1">¹All computations in this paper have been performed in Maple using the C library msolve [6].</footnote>
where $M := m_1 + \dots + m_n$. Let $\epsilon$ be a new variable, $\mathbb{L} := \mathbb{K}(\epsilon)$, and let $(\gamma_{i,j})_{1\le i,j\le n}$ be defined by $\gamma_{i,i} = i^k$ and $\gamma_{i,j} = t^\beta$ for $i \ne j$. Then, consider the following system which is a symbolic deformation of (5) with respect to the deformation parameter $\epsilon$:

$$ \left\{ \begin{array}{l} (\mathbf{E}_{G_1}): \quad G_1 = f_1(u) + t^\alpha \cdot Q_1(\nabla^k G_1, \nabla^k G_2, \dots, \nabla^k G_n, t^\alpha, u) + t \cdot \epsilon^k \cdot \sum_{i=1}^n \gamma_{1,i} \cdot \Delta^k G_i, \\ \vdots \\ (\mathbf{E}_{G_n}): \quad G_n = f_n(u) + t^\alpha \cdot Q_n(\nabla^k G_1, \nabla^k G_2, \dots, \nabla^k G_n, t^\alpha, u) + t \cdot \epsilon^k \cdot \sum_{i=1}^n \gamma_{n,i} \cdot \Delta^k G_i. \end{array} \right. \quad (7) $$

The fixed point nature of the above equations still implies the existence of a unique solution $(G_1, \dots, G_n) \in \mathbb{L}[u][[t]]^n$ (it can be seen by extracting the coefficient of $t^m$ of each $G_1, \dots, G_n$). Remark that the equalities $F_i(t^\alpha, u) = G_i(t, u, 0)$ relate the formal power series solutions of (5) and of (7). Hence, showing that each $G_i$ is algebraic over $\mathbb{L}(t, u)$ is enough to prove theorem 3. Moreover, as we will see later, the algebraicity of each $G_i$ follows from the algebraicity of the series $G_1(0), \dots, \partial_u^{k-1}G_1(0), \dots, G_n(0), \dots, \partial_u^{k-1}G_n(0)$. Here, and in what follows, we shall use the short notations

$$ G_i(u) \equiv G_i(t, u, \epsilon), \partial_0 G_i(u) \equiv G_i(u), G_i(0), \partial_u G_i(0), \dots, \partial_u^{k-1} G_i(0) $$

and $A(u) \equiv A(\partial_0 G_1, \dots, \partial_0 G_n, t, u)$ for any polynomial $A \in \mathbb{L}[X_1, \dots, X_n, t, u]$ with the notation $X_j := x_j, z_{k(j-1)}, z_{k(j-1)+1}, \dots, z_{k_{j-1}}$. Note that in the case $n=1$, this notation implies that for any $0 \le i \le k-1$, the variable $z_i$ stands for $\partial_u^i F_1(t, a)$.

Let us define $Y_{i,0} := x_i$ and $Y_{i,j} := (x_i - z_{k(i-1)} - \dots - \frac{u^{j-1}}{(j-1)!}z_{k(i-1)+j-1})/u^j$ for $1 \le i \le n$ and $1 \le j \le k$. With these definitions, multiplying each $(\mathbf{E}_{F_i})$ in (5) by $u^{m_i}$ and substituting the series $G_i$'s and their specializations by their associated variables yields the following system of polynomial equations

$$ \left\{ \begin{array}{l} E_1 := u^{m_1} \cdot (f_1(u) - x_1 + t^\alpha \cdot Q_1(Y_{1,0}, \dots, Y_{1,k}, Y_{2,0}, \dots, Y_{n,k}, t^\alpha, u) + t \cdot \epsilon^k \cdot \sum_{i=1}^n \gamma_{1,i} \cdot Y_{i,k}) = 0, \\ \vdots \\ E_n := u^{m_n} \cdot (f_n(u) - x_n + t^\alpha \cdot Q_n(Y_{1,0}, \dots, Y_{1,k}, Y_{2,0}, \dots, Y_{n,k}, t^\alpha, u) + t \cdot \epsilon^k \cdot \sum_{i=1}^n \gamma_{n,i} \cdot Y_{i,k}) = 0. \end{array} \right. \quad (8) $$

Like in (6), we take the derivative with respect to the variable $u$ of these equations and find

$$ \left( \begin{array}{ccc} (\partial_{x_1} E_1)(u) & \dots & (\partial_{x_n} E_1)(u) \\ \vdots & \ddots & \vdots \\ (\partial_{x_1} E_n)(u) & \dots & (\partial_{x_n} E_n)(u) \end{array} \right) \cdot \left( \begin{array}{c} \partial_u G_1 \\ \vdots \\ \partial_u G_n \end{array} \right) + \left( \begin{array}{c} (\partial_u E_1)(u) \\ \vdots \\ (\partial_u E_n)(u) \end{array} \right) = 0. \quad (9) $$

Let $\text{Det} \in \mathbb{L}[X_1, \dots, X_n, t][u]$ be the determinant of the square matrix $(\partial_{x_j} E_i)_{1\le i,j\le n}$. The following lemma gives the number of distinct relevant solutions in $u$ to the equation $\text{Det}(u) = 0$.

**Lemma 6.** $\text{Det}(u) = 0$ admits exactly $nk$ distinct nonzero solutions $U_1, \dots, U_{nk} \in \overline{\mathbb{L}}[[t^\frac{1}{k}]]$.

*Proof.* Note that we have

$$ \text{Det}(u) = \det \left( \begin{array}{ccc} -u^{m_1} + t\epsilon^k \gamma_{1,1} u^{m_1-k} & \dots & t\epsilon^k \gamma_{1,n} u^{m_1-k} \\ \vdots & \ddots & \vdots \\ t\epsilon^k \gamma_{n,1} u^{m_n-k} & \dots & -u^{m_n} + t\epsilon^k \gamma_{n,n} u^{m_n-k} \end{array} \right) + O(t^\alpha u^{M-nk}). $$

For every $i$ we first divide the $i$th row by $u^{m_i-k}$. Then, using the definition of $\gamma_{i,j}$ and $\alpha, \beta \ge n$, we see that the matrix above becomes diagonal mod $t^{n+1}$ and its determinant mod $t^{n+1}$ simplifies to
$u^{M-nk} \cdot \prod_{j=1}^{n} (-u^k + t\epsilon^k j^k) \mod t^{n+1}$. Hence, computing the first terms of a solution in $u$ by using Newton polygons, we find $nk$ distinct nonzero solutions in $u$ to the equation $\text{Det}(u) = 0$. We denote these solutions by $U_1, \dots, U_{nk} \in \overline{\mathbb{L}}[[t^\frac{1}{k}]]$. Their first terms are given by $\zeta^\ell \cdot t^\frac{1}{k} \cdot \epsilon, \dots, \zeta^\ell \cdot n \cdot t^\frac{1}{k} \cdot \epsilon \in \overline{\mathbb{L}}[[t^\frac{1}{k}]]$, for $\zeta$ a $k$-primitive root of unity and for all $1 \le \ell \le k$. Finally, note that the constant coefficient in $t$ of $\prod_{j=1}^{n} (-u^k + t\epsilon^k j^k)$ has degree $nk$ so by [10, Theorem 2] there cannot be more than $nk$ solutions to $\text{Det}(u) = 0$ in $\overline{\mathbb{L}}[[t^\frac{1}{k}]] \setminus \{0\}$. $\square$

Now, let $P$ be the determinant of the square matrix $(\partial_{x_j} E_i)_{1\le i,j \le n}$ where the last column $(\partial_{x_n} E_1, \dots, \partial_{x_n} E_n)$ is replaced by $(\partial_u E_1, \dots, \partial_u E_n)$, that is

$$ P := \det \begin{pmatrix} \partial_{x_1} E_1 & \dots & \partial_{x_{n-1}} E_1 & \partial_u E_1 \\ \vdots & \ddots & \vdots & \vdots \\ \partial_{x_1} E_n & \dots & \partial_{x_{n-1}} E_n & \partial_u E_n \end{pmatrix}. $$

1. Clearly, if $\text{Det}(u) = 0$ then (9) implies $P(u) = 0$, thus $P(u)$ vanishes at the roots $U_1, \dots, U_{nk}$ in lemma 6. Hence, defining the polynomial system $S$ in $\mathbb{L}[t][X_1, \dots, X_n, u]$ given by the vanishing of the set of polynomials $(E_1, \dots, E_n, \text{Det}, P)$, we see that $S$ is a system with exactly $n+2$ equations in the $nk+n+1$ variables given by $z_0, \dots, z_{nk-1}, x_1, \dots, x_n, u$ (here $t$ and $\epsilon$ are parameters). We introduce the duplicated system $S_{\text{dup}}^{\epsilon} := (S_1, \dots, S_{nk})$ in $\mathbb{L}(t)[x_1, \dots, x_{n^2k}, z_0, \dots, z_{nk-1}, u_1, \dots, u_{nk}]$ after duplicating $nk$ times each of the variables $x_i$'s, $u_i$'s and after duplicating $nk$ times the initial polynomial system $S$: all in all, we perform in spirit step 4 of Algorithm 1. This system is built from $nk(n+2)$ equations and $nk(n+2)$ variables.

The following lemma is proven in [8, Lemma 2.10] as a consequence of Hilbert's Nullstellensatz and of the Jacobian criterion [19, Theorem 16.19]. Recall that for an integer $N > 1$ and for some polynomial $g \in \mathbb{K}[x_1, \dots, x_N]$ and an ideal $\mathcal{I} \subset \mathbb{K}[x_1, \dots, x_N]$, the saturation of $\mathcal{I}$ by $g$ (also called saturated ideal) is defined by $\mathcal{I} : g^{\infty} := \{f \in \mathbb{K}[x_1, \dots, x_N] | \exists s \in \mathbb{N} \text{ s.t. } g^s \cdot f \in \mathcal{I}\}^2$.

**Lemma 7.** Assume that the Jacobian matrix $Jac_{S_{\text{dup}}^{\epsilon}}$ of $S_{\text{dup}}^{\epsilon}$ considered with respect to the variables $x_1, \dots, x_n, u_1, \dots, x_{n^2k-n+1}, \dots, x_{n^2k}, u_{nk}, z_0, \dots, z_{nk-1}$, is invertible at the point

$$ \mathcal{P} = (G_1(U_1), \dots, G_n(U_1), U_1, \dots, G_1(U_{nk}), \dots, G_n(U_{nk}), U_{nk}, G_1(0), \dots, \partial_u^{k-1} G_1(0), \dots, G_n(0), \dots, \partial_u^{k-1} G_n(0)) \in \overline{\mathbb{L}}[[t^\frac{1}{k}]]^{nk(n+1)} \times \mathbb{L}[[t]]^{nk}. $$

Denote $\mathcal{I}_{\text{dup}}^{\epsilon}$ the ideal of $\mathbb{L}(t)[x_1, \dots, x_n, u_1, \dots, x_{n^2k-n+1}, \dots, x_{n^2k}, u_{nk}, z_0, \dots, z_{nk-1}]$ generated by $S_{\text{dup}}^{\epsilon}$. Then the saturated ideal $\mathcal{I}_{\text{dup}}^{\epsilon} : \det(Jac_{S_{\text{dup}}^{\epsilon}})^{\infty}$ is zero-dimensional and radical over $\mathbb{L}(t)$. Moreover, $\mathcal{P}$ lies in the zero set of $\mathcal{I}_{\text{dup}}^{\epsilon} : \det(Jac_{S_{\text{dup}}^{\epsilon}})^{\infty}$.

Therefore, in order to conclude the algebraicity of $G_i(0), \dots, \partial_u^{k-1} G_i(0)$ over $\mathbb{L}(t)$ for all $1 \le i \le n$, it is enough to justify that $Jac_{S_{\text{dup}}^{\epsilon}}$ is invertible at $\mathcal{P}$. The idea for proving $\det(Jac_{S_{\text{dup}}^{\epsilon}})(\mathcal{P}) \ne 0$, analogous to the proof of [10, Theorem 3], is to show first that $Jac_{S_{\text{dup}}^{\epsilon}}(\mathcal{P})$ can be rewritten as a block triangular matrix. We will then show that the diagonal blocks are invertible by carefully analyzing the lowest valuation in $t$ of their associated determinants.

If $A \in \mathbb{L}[t][X_1, \dots, X_n, u]$, we shall define its “$i$th duplicated polynomial” by

$$ A^{(i)} := A(X_{n(i-1)+1}, \dots, X_{ni}, u_i). $$

<footnote n="2">²In practice, if $\{h_1, \dots, h_r\}$ is a generating set of $\mathcal{I}$, then a generating set of $\mathcal{I} : g^{\infty}$ is obtained by computing a generating set of $\langle h_1, \dots, h_r, m \cdot g - 1 \rangle \cap \mathbb{K}[x_1, \dots, x_N]$, where $m$ is an extra variable.</footnote>
Then the Jacobian matrix $Jac_{S_{\text{dup}}^{\epsilon}}(\mathcal{P})$ has the shape

$$Jac_{S_{\text{dup}}} (\mathcal{P}) = \begin{pmatrix} A_1 & & 0 & B_1 \\ & \ddots & & \vdots \\ 0 & & A_{nk} & B_{nk} \end{pmatrix} \in \overline{\mathbb{L}}[[t^\frac{1}{\ell}]]^{nk(n+2)\times nk(n+2)},$$

where for $i = 1, \dots, nk$ the matrices$^3 A_i \in \overline{\mathbb{L}}[[t^\frac{1}{\ell}]]^{(n+2)\times(n+1)}$ and $B_i \in \overline{\mathbb{L}}[[t^\frac{1}{\ell}]]^{(n+2)\times nk}$ are:

$$A_i := \begin{pmatrix} \partial_{x_1} E_1^{(i)}(U_i) & \dots & \partial_{x_n} E_1^{(i)}(U_i) & \partial_{u_i} E_1^{(i)}(U_i) \\ \vdots & \ddots & \vdots & \vdots \\ \partial_{x_1} E_n^{(i)}(U_i) & \dots & \partial_{x_n} E_n^{(i)}(U_i) & \partial_{u_i} E_n^{(i)}(U_i) \\ \partial_{x_1} \text{Det}^{(i)}(U_i) & \dots & \partial_{x_n} \text{Det}^{(i)}(U_i) & \partial_{u_i} \text{Det}^{(i)}(U_i) \\ \partial_{x_1} P^{(i)}(U_i) & \dots & \partial_{x_n} P^{(i)}(U_i) & \partial_{u_i} P^{(i)}(U_i) \end{pmatrix},$$ $$B_i := \begin{pmatrix} \partial_{z_0} E_1^{(i)}(U_i) & \dots & \partial_{z_{nk-1}} E_1^{(i)}(U_i) \\ \vdots & \ddots & \vdots \\ \partial_{z_0} E_n^{(i)}(U_i) & \dots & \partial_{z_{nk-1}} E_n^{(i)}(U_i) \\ \partial_{z_0} \text{Det}^{(i)}(U_i) & \dots & \partial_{z_{nk-1}} \text{Det}^{(i)}(U_i) \\ \partial_{z_0} P^{(i)}(U_i) & \dots & \partial_{z_{nk-1}} P^{(i)}(U_i) \end{pmatrix}.$$

Using $\text{Det}(U_i) = 0$ and (9), we see that the first $n \times (n+1)$ minor of each $A_i$ has rank at most $n-1$. Hence, after performing operations on the first $n$ rows, we can transform the $n$ th row of $A_i$ into the zero vector. It follows that after the suitable transformation and a permutation of rows, $Jac_{S_{\text{dup}}^{\epsilon}}(\mathcal{P})$ can be rewritten as a block triangular matrix. To give the precise form of the determinant of $Jac_{S_{\text{dup}}^{\epsilon}}(\mathcal{P})$, we first define

$$R := \det \begin{pmatrix} \partial_{x_1} E_1^{(i)}(U_i) & \dots & \partial_{x_{n-1}} E_1^{(i)}(U_i) & y_1 \\ \vdots & \ddots & \vdots & \vdots \\ \partial_{x_1} E_n^{(i)}(U_i) & \dots & \partial_{x_{n-1}} E_n^{(i)}(U_i) & y_n \end{pmatrix} \in \mathbb{K}[\{\partial_{x_i} E_j^{(i)}(U_i)\}_{1 \le j \le n, 1 \le \ell \le n-1}] [y_1, \dots, y_n]. \quad (10)$$

Then it follows that $\det(Jac_{S_{\text{dup}}^{\epsilon}})(\mathcal{P}) = \pm (\prod_{i=1}^{nk} \det(Jac_i(U_i))) \cdot \det(\Lambda)$, where

$$Jac_i(u) := \begin{pmatrix} \partial_{x_1} E_1^{(i)}(u) & \dots & \partial_{x_n} E_1^{(i)}(u) & \partial_{u_i} E_1^{(i)}(u) \\ \vdots & \ddots & \vdots & \vdots \\ \partial_{x_1} E_{n-1}^{(i)}(u) & \dots & \partial_{x_n} E_{n-1}^{(i)}(u) & \partial_{u_i} E_{n-1}^{(i)}(u) \\ \partial_{x_1} \text{Det}^{(i)}(u) & \dots & \partial_{x_n} \text{Det}^{(i)}(u) & \partial_{u_i} \text{Det}^{(i)}(u) \\ \partial_{x_1} P^{(i)}(u) & \dots & \partial_{x_n} P^{(i)}(u) & \partial_{u_i} P^{(i)}(u) \end{pmatrix} \in \overline{\mathbb{L}}[u][[t]]^{(n+1)\times(n+1)}, \text{ and}$$ $$\Lambda := (R(\partial_{z_j} E_1^{(i)}(U_i), \dots, \partial_{z_j} E_n^{(i)}(U_i)))_{1 \le i, j+1 \le nk} \in \overline{\mathbb{L}}[[t^\frac{1}{\ell}]]^{nk \times nk}. \quad (11)$$

The proof that this product is nonzero is the content of lemma 8 and lemma 9.

**Lemma 8.** For each $i = 1, \dots, nk$, the determinant of $Jac_i(U_i)$ is nonzero.

<footnote>$^3$In these matrices, we emphasize that notations like $\partial_{x_1} E_1^{(i)}(U_i)$ are compact forms for the specializations of the duplicated polynomial $\partial_{x_1} E_1^{(i)}$ to the values $x_{(i-1)n+1} = F_1(t, U_i(t)), \dots, x_{in} = F_n(t, U_i(t)), u_i = U_i(t), z_{(j-1)k+\ell} = (\partial_{u_i} F_j)(t, a)$.</footnote>
Proof. To prove that $\det(\mathrm{Jac}_i(U_i)) \neq 0$ we will show that $\mathrm{val}_t(\det(\mathrm{Jac}_i(U_i))) < \infty$, where $\mathrm{val}_t$ denotes the valuation in $t$. The main idea here is to expand $\det(\mathrm{Jac}_i(U_i))$ with respect to the last column and show that the least valuation comes from the product of $\partial_{u_i} \det^{(i)}(U_i)$ by the determinant of its associated submatrix$^4$, denoted by $\mathcal{M}$. For some matrix $\mathcal{A}$ whose entries are series in $t$, we shall denote by $\mathrm{val}_t(\mathcal{A})$ the matrix of the valuations in $t$ of the entries of $\mathcal{A}$.

We shall justify that the term with lowest exponent in $t$ in $\det(\mathrm{Jac}_i(U_i))$ comes from the first term in the product $\det(\mathcal{M}) \cdot \partial_{u_i} \det^{(i)}(U_i)$. By construction, the monomials in $\{E_i\}_{1\le i\le n}$ that are quadratic in the variables $\{x_j\}_{1\le j\le n}$ all carry a $t^\alpha$, so that $\partial_{x_j} \det^{(i)}(U_i) = O(t^\alpha)$ for all $1 \le j \le n$ and all $1 \le i \le nk$. Thus when expending $\det(\mathrm{Jac}_i(U_i))$ with respect to the last column of $\mathrm{Jac}_i(U_i)$, the minors associated with $\partial_{u_i} E_j^{(i)}(U_i)$ and $\partial_{u_i} P^{(i)}(U_i)$ are in $O(t^\alpha)$ (for $1 \le i \le nk$).

Thus, it remains to prove that $\mathrm{val}_t(\partial_{u_i} \det^{(i)}(U_i) \cdot \det(\mathcal{M})) < \alpha$; this is done in the remaining part of the proof. Note that the expression of $\mathrm{Det} \mod t^{n+1}$ in the proof of lemma 6 implies that

$$ \mathrm{val}_t(\partial_{u_i} \det^{(i)}(U_i)) = (M-nk)/k + n-1 + (k-1)/k = (M-1)/k. \quad (12) $$

For $\det(\mathcal{M})$ there are two cases to treat separately: Either **Case 1** we have $U_i^k = n^k t \epsilon^k + (\text{higher powers of } t)$, or **Case 2** we have $U_i^k = m^k t \epsilon^k + (\text{higher powers of } t)$ for some $m < n$. The reason for this distinction is that

$$ \partial_{x_j} E_{\ell}(u) = \begin{cases} -u^{m_{\ell}} + t \epsilon^k u^{m_{\ell}-k} \ell^k, & \text{if } \ell = j, \\ t^{\beta+1} \epsilon^k u^{m_{\ell}-k}, & \text{else,} \end{cases} \quad (13) $$

and, therefore, in **Case 1** the $(\ell, \ell)$ entry of $\mathcal{M}$ always has valuation in $t$ given by $m_\ell/k$ for each $\ell = 1, \dots, n-1$, while in **Case 2** the entry of $\mathcal{M}$ on row and column $m$ has a valuation in $t$ that depends on $\beta$.

**Case 1:** Assume that the $U_i$ of interest satisfies

$$ U_i^k = n^k t \epsilon^k + (\text{higher powers of } t). \quad (14) $$

By definition, $\mathcal{M}_{\ell,j} = \partial_{x_j} E_\ell(U_i)$ for $j = 1, \dots, n$ and $\ell = 1, \dots, n-1$ and thus from (13) we find

$$ \mathrm{val}_t(\mathcal{M})_{\ell,j} = \begin{cases} \frac{m_\ell}{k} & \text{if } \ell = j, \\ \beta + \frac{m_\ell}{k} & \text{else.} \end{cases} \quad (15) $$

Moreover, we claim that the bottom right entry of $\mathcal{M}$ has valuation $(M-1)/k$ in $t$. To prove this we compute

$$ \begin{aligned} \mathcal{M}_{n,n} = \partial_{x_n} P^{(i)}(U_i) &= \partial_{x_n} \det \begin{pmatrix} \partial_{x_1} E_1 & \dots & \partial_{x_{n-1}} E_1 & \partial_{u} E_1 \\ \vdots & \ddots & \vdots & \vdots \\ \partial_{x_1} E_n & \dots & \partial_{x_{n-1}} E_n & \partial_{u} E_n \end{pmatrix}^{(i)} (U_i) \\ &= \det \begin{pmatrix} \partial_{x_1} E_1 & \dots & \partial_{x_{n-1}} E_1 & \partial_{x_n, u}^2 E_1 \\ \vdots & \ddots & \vdots & \vdots \\ \partial_{x_1} E_n & \dots & \partial_{x_{n-1}} E_n & \partial_{x_n, u}^2 E_n \end{pmatrix}^{(i)} (U_i) \mod t^\alpha, \end{aligned} \quad (16) $$

<footnote>$^4$For $A = (a_{i,j})_{1\le i,j\le n}$ some matrix of size $n \times n$, the submatrix associated to an element $a_{i_0,j_0}$ is the $(n-1) \times (n-1)$ submatrix of $A$ obtained by deletion of the $i_0$th row and $j_0$th column from $A$.</footnote>
since $\partial_{x_n,x_j}^2 E_\ell \in O(t^\alpha)$ for each $j = 1, \dots, n$ and $\ell = 1, \dots, n$. Regarding the last column of the above matrix, a straightforward computations yields

$$ \partial_{x_n, u}^2 E_n(u) = -m_n u^{m_n-1} + t(m_n-k)n^k u^{m_n-k-1} + O(t^\alpha). \quad (17) $$

From this and (14) it follows that $\text{val}_t((\partial_{x_n,u}^2 E_\ell)^{(i)}(U_i)) = (m_n - 1)/k$. Moreover, using (13) and (15), we see that the only way to obtain in (16) a term with exponent in $t$ independent of $\beta$ is to take the product of entries of the main diagonal. The sufficiently large choice of $\beta = [2M/k]$ implies that $\text{val}_t(\mathcal{M}_{n,n}) = \text{val}_t((\partial_{x_n} P)^{(i)}(U_i)) = \frac{m_1}{k} + \dots + \frac{m_{n-1}}{k} + \frac{m_{n-1}}{k} = (M-1)/k$, where, as before, $M$ denotes $\sum_{j=1}^n m_j$.

Altogether, in this case we obtain that $\text{val}_t(\mathcal{M})$ has the shape

$$ \text{val}_t(\mathcal{M}) = \begin{pmatrix} \frac{m_1}{k} & \beta + \frac{m_1}{k} & \dots & \beta + \frac{m_1}{k} & \beta + \frac{m_1}{k} \\ \beta + \frac{m_2}{k} & \frac{m_2}{k} & \dots & \beta + \frac{m_2}{k} & \beta + \frac{m_2}{k} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ \beta + \frac{m_{n-1}}{k} & \beta + \frac{m_{n-1}}{k} & \dots & \frac{m_{n-1}}{k} & \beta + \frac{m_{n-1}}{k} \\ * & * & \dots & * & (M-1)/k \end{pmatrix} $$

It follows that the only term in the determinant of $\mathcal{M}$ whose exponent in $t$ has no dependency on $\beta$ comes from the product of the entries located on the main diagonal of $\mathcal{M}$, and that the choice $\beta = [2M/k]$ ensures that $\text{val}_t(\det(\mathcal{M})) < \beta$. Thus (12) implies $\text{val}_t(\partial_{u_i} \det^{(i)}(U_i)) \cdot \det(\mathcal{M})) = (M-1)/k + (2M - m_n - 1)/k$. Finally, from $\alpha = 3n^2k \cdot ([2M/k] + 1) + 3nM$, it follows that $\text{val}_t(\partial_{u_i} \det^{(i)}(U_i)) \cdot \det(\mathcal{M})) < \alpha$, as wanted.

**Case 2:** The second case is similar in spirit to **Case 1**, but because of (13) when $U_i$ is not of the form (14), computing the valuation in $t$ of $\det(\mathcal{M})$ is slightly more delicate. In this case

$$ U_i^k = m^k t \epsilon^k + (\text{higher powers of } t), \quad (18) $$

for some $m < n$. For the sake of better readability, we shall assume without loss of generality that $m = 1$, since the argument works equally well in the general case.

It follows from the definition of $\mathcal{M}$ and (13) that

$$ \mathcal{M} = \begin{pmatrix} -U_i^{m_1} + t\epsilon^k U_i^{m_1-k} & \dots & \epsilon^k t^{\beta+1} U_i^{m_1-k} & \epsilon^k t^{\beta+1} U_i^{m_1-k} \\ \vdots & \ddots & \vdots & \vdots \\ \epsilon^k t^{\beta+1} U_i^{m_{n-1}-k} & \dots & -U_i^{m_{n-1}} + (n-1)^k t \epsilon^k U_i^{m_{n-1}-k} & \epsilon^k t^{\beta+1} U_i^{m_{n-1}-k} \\ \partial_{x_1} P^{(i)}(U_i) & \dots & \partial_{x_{n-1}} P^{(i)}(U_i) & \partial_{x_n} P^{(i)}(U_i) \end{pmatrix} \mod t^\alpha. $$

As before in (15), for $1 \le j \le n$, we have $\text{val}_t(\epsilon^k t^{\beta+1} U_i^{m_j-k}) = \beta + \frac{m_j}{k}$. Also, assuming that $m = 1$, we have for all $2 \le j \le n-1$ that

$$ -U_i^{m_j} + j^k t \epsilon^k U_i^{m_j-k} = \lambda_{i,j} \cdot t^{\frac{m_j}{k}} + (\text{higher powers of } t), $$

for some nonzero $\lambda_{i,j} \in \overline{\mathbb{K}}(\epsilon)$. However, for $j = 1$ the term $\lambda_{i,j}$ vanishes and so we shall compute the valuation in $t$ of $-U_i^{m_1} + t \epsilon^k U_i^{m_1-k}$ in this case. From the expansion of $\det(u)$ to higher order terms

$$ \det(u) = u^{M-nk} \cdot \left[ \prod_{1 \le j \le n} (-u^k + j^k t \epsilon^k) - t^{2(\beta+1)} \epsilon^{2k} \sum_{1 \le \ell < j \le n} \prod_{b=1, b \ne \ell, b \ne j}^{n} (-u^k + b^k t \epsilon^k) \right] + O(t^{3\beta}), $$
we use Newton's method to find the second lowest term of $U_i$ with respect to the exponent in $t$:

$$U_i = \zeta^{\ell_i} t^{\frac{1}{k}} \epsilon + \lambda_i t^{2\beta + \frac{1}{k}} + (\text{higher powers of } t),$$

for some $1 \le \ell_i \le k$ and $\lambda_i \ne 0$. This implies that

$$\begin{aligned} \mathrm{val}_t(\mathcal{M}_{1,1}) &= \mathrm{val}_t(-U_i^{m_1} + t\epsilon^k U_i^{m_1-k}) = \mathrm{val}_t(U_i^{m_1-k}) + \mathrm{val}_t(-U_i^k + t\epsilon^k) \\ &= \frac{m_1-k}{k} + \frac{(k-1)}{k} + 2\beta + \frac{1}{k} = 2\beta + \frac{m_1}{k}. \end{aligned}$$

It remains to understand the valuation in $t$ of the last row of $\mathcal{M}$, i.e. of $\{\partial_{x_j} P^{(i)}(U_i)\}_{1 \le j \le n}$. The same argument as in (16) implies that for any $j = 1, \dots, n$ we have

$$(\partial_{x_j} P)^{(i)}(U_i) = \det \left( \begin{array}{cccc} -U_i^{m_1} + t\epsilon^k U_i^{m_1-k} & \dots & \epsilon^k t^{\beta+1} U_i^{m_1-k} & \partial_{x_j, u}^2 E_1^{(i)}(U_i) \\ \vdots & \ddots & \vdots & \vdots \\ \epsilon^k t^{\beta+1} U_i^{m_n-k} & \dots & \epsilon^k t^{\beta+1} U_i^{m_n-k} & \partial_{x_j, u}^2 E_n^{(i)}(U_i) \end{array} \right) \mod t^\alpha.$$

Moreover, in spirit of (17) it holds because of (13) and (18) that

$$\partial_{u,x_j}^2 E_\ell(U_i) = \begin{cases} (m_\ell - 1)/k & \text{if } j = \ell, \\ O(t^{\beta + \frac{m_\ell - 1}{k}}) & \text{else.} \end{cases}$$

Putting everything together, it follows that

$$\mathrm{val}_t(\partial_{x_j} P^{(i)}(U_i)) = \mathrm{val}_t \det \begin{pmatrix} t^{2\beta + \frac{m_1}{k}} & \dots & \dots & t^{\beta + \frac{m_1}{k}} & t^{\beta + \frac{m_1-1}{k}} \\ t^{\beta + \frac{m_2}{k}} & t^{\frac{m_2}{k}} & t^{\beta + \frac{m_1}{k}} & \dots & t^{\beta + \frac{m_2-1}{k}} \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ t^{\beta + \frac{m_j}{k}} & \dots & t^{\frac{m_j}{k}} & t^{\beta + \frac{m_j}{k}} & t^{(m_j-1)/k} \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ t^{\beta + \frac{m_{n-1}}{k}} & \dots & t^{\beta + \frac{m_{n-1}}{k}} & t^{\frac{m_{n-1}}{k}} & t^{\beta + \frac{m_{n-1}-1}{k}} \\ t^{\beta + \frac{m_n}{k}} & \dots & t^{\beta + \frac{m_n}{k}} & t^{\beta + \frac{m_n}{k}} & t^{\beta + \frac{m_n-1}{k}} \end{pmatrix},$$

where in the matrix above the $(i_0, j_0)$ entry is $t^{\beta+m_{i_0}/k}$ if $1 \le j_0 < n$ and $1 \le i_0 \ne j_0 \le n$, $t^{m_{i_0}/k}$ if $1 < i_0 = j_0 < n$, $t^{2\beta+m_{i_0}/k}$ if $i_0 = j_0 = 1$, $t^{\beta+\frac{m_{i_0}-1}{k}}$ if $j_0 = n, i_0 \ne j$, and $t^{(m_{j-1})/k}$ if $j_0 = n, i_0 = j$. From this we obtain that $\mathrm{val}_t(\partial_{x_j} P^{(i)}(U_i))$ is at least $2\beta + (M-1)/k$, except if $j=1$, since then it is given by the valuation of the product of the lower left entry, the upper right entry and the remaining entries on the diagonal of the matrix above. In other words:

$$\mathrm{val}_t(\partial_{x_j} P^{(i)}(U_i)) = \begin{cases} \beta + (M-1)/k & \text{if } j=1, \\ \ge 2\beta + (M-1)/k & \text{else.} \end{cases}$$

This means that $\mathrm{val}_t(\mathcal{M})$ has the shape

$$\left( \begin{array}{cccc} 2\beta + \frac{m_1}{k} & \beta + \frac{m_1}{k} & \cdots & \cdots & \beta + \frac{m_1}{k} \\ \beta + \frac{m_2}{k} & \frac{m_2}{k} & \beta + \frac{m_2}{k} & \cdots & \beta + \frac{m_2}{k} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \beta + \frac{m_{n-2}}{k} & \cdots & \frac{m_{n-2}}{k} & \beta + \frac{m_{n-2}}{k} & \beta + \frac{m_{n-2}}{k} \\ \beta + \frac{m_{n-1}}{k} & \cdots & \beta + \frac{m_{n-1}}{k} & \frac{m_{n-1}}{k} & \beta + \frac{m_{n-1}}{k} \\ \beta + \frac{M-1}{k} & \ge 2\beta + \frac{M-1}{k} & \cdots & \ge 2\beta + \frac{M-1}{k} & \ge 2\beta + \frac{M-1}{k} \end{array} \right).$$
We see that the valuation of $\det(\mathcal{M})$ is $2\beta + (2M - m_n - 1)/k$, since it is given by the sum of the lower left entry, the upper right entry and the remaining entries on the diagonal of the matrix above. Note that each other combination will sum to at least $4\beta$ and the choice $\beta = \lfloor 2M/k \rfloor$ ensures that this term is strictly larger.

It remains to check that also in this case our choice for $\alpha$ was large enough to guarantee that $\mathrm{val}_t(\partial_{u_i} \mathrm{Det}^{(i)}(U_i) \cdot \det(\mathcal{M})) < \alpha$. Using that $\mathrm{val}_t(\det(\mathcal{M})) = 2\beta + (2M - m_n - 1)/k$ and (12), we obtain

$$ \mathrm{val}_t(\partial_{u_i} \mathrm{Det}^{(i)}(U_i)) + \mathrm{val}_t(\det(\mathcal{M})) = 2\beta + \frac{3M - m_n - 2}{k} < 2\beta + \frac{3M}{k}. $$

As $2\beta + \frac{3M}{k} < \alpha$, this concludes the proof of lemma 8. □

**Lemma 9.** The determinant of $\Lambda$ is nonzero.

Proving that $\det(\Lambda) \neq 0$ is again done by analyzing the first terms in $t$ of $\det(\Lambda)$. More precisely, we are going to show that $\det(\Lambda)$ factors modulo $t^\alpha$ as: a product of $U_i$, the Vandermonde determinant $\prod_{i<j}(U_i - U_j)$, and a nonzero polynomial $H(t) \in \mathbb{K}[t]$.

Before starting with the proof we shall prove the following simple but useful fact:

**Lemma 10.** Let $\mathbb{F}$ be a field and $m \in \mathbb{N}$. Consider polynomials $A_1, \dots, A_m \in \mathbb{F}[x]$ and define $M := (A_i(x_j))_{1 \leq i,j \leq m}$. Then $\det(M) = 0$ if and only if there exist $\lambda_1, \dots, \lambda_m \in \mathbb{F}$ not all equal to 0 such that we have the linear combination $\sum_{i=1}^m \lambda_i \cdot A_i(x) = 0$.

*Proof.* Clearly, if the linear combination $\sum_{i=1}^m \lambda_i \cdot A_i(x) = 0$ exists, $M$ is singular and its determinant vanishes. For the other direction, we can assume without loss of generality that, up to permutation of columns of $M$ and linear operations on them, $\deg(A_1(x)) > \dots > \deg(A_n(x))$. Then, since $\det(M) = 0$, we must have that $A_i(x) = 0$ for some $i = 1, \dots, n$, because the product of the diagonal elements of $M$ cannot cancel otherwise. □

*Proof of lemma 9.* Recall from (10) and (11) that

$$ \Lambda := (R(\partial_{z_j} E_1^{(i)}(U_i), \dots, \partial_{z_j} E_n^{(i)}(U_i)))_{1 \le i,j+1 \le nk} \in \overline{\mathbb{L}}[[t^\frac{1}{\alpha}]]^{nk \times nk}, $$

where

$$ R := \det \begin{pmatrix} \partial_{x_1} E_1^{(i)}(U_i) & \dots & \partial_{x_{n-1}} E_1^{(i)}(U_i) & y_1 \\ \vdots & \ddots & \vdots & \vdots \\ \partial_{x_1} E_n^{(i)}(U_i) & \dots & \partial_{x_{n-1}} E_n^{(i)}(U_i) & y_n \end{pmatrix} \in \mathbb{K}[\{\partial_{x_\ell} E_j^{(i)}(U_i)\}_{1 \le \ell \le n-1, 1 \le j \le n}\][y_1, \dots, y_n]. $$

We shall first define and analyze the symbolic matrix $\tilde{\Lambda}$:

$$ \tilde{\Lambda}(u_1, \dots, u_n) = \tilde{\Lambda} := (\tilde{R}_j^{(i)})_{1 \le i,j+1 \le nk} \in \mathbb{L}[t][u_1, \dots, u_{nk}]^{nk \times nk}, \text{ where} $$

$$ \tilde{R}_j := \det \begin{pmatrix} \partial_{x_1} E_1 & \dots & \partial_{x_{n-1}} E_1 & \partial_{z_j} E_1 \\ \vdots & \ddots & \vdots & \vdots \\ \partial_{x_1} E_n & \dots & \partial_{x_{n-1}} E_n & \partial_{z_j} E_n \end{pmatrix} \mod t^\alpha \in \mathbb{L}[t][u]. $$

Note that the polynomial matrix $\tilde{\Lambda}$ is symbolic and it holds that

$$ \tilde{\Lambda}(U_1, \dots, U_{nk}) \equiv \Lambda \mod t^{\alpha}. \qquad (19) $$

We will prove lemma 9 in 3 steps corresponding to the 3 claims:

<footnote n="5">⁵Here $\deg(0) = -\infty$ and we allow us to write $-\infty > -\infty$.</footnote>
• Claim 1: For the symbolic matrix $\tilde{\Lambda}$ it holds that $\det(\tilde{\Lambda}) \neq 0$.

• Claim 2: We have that $\det(\tilde{\Lambda}) = \prod_{i=1}^{nk} u_i^{M-nk} \cdot \prod_{i<j} (u_i - u_j) \cdot H(t)$, for some nonzero polynomial $H(t) \in \mathbb{L}[t]$ of degree bounded by $n^2k(\beta + 1)$.

• Claim 3: It holds that $\det(\tilde{\Lambda}(U_1, \dots, U_{nk})) \mod t^\alpha \neq 0$, in particular, we have $\Lambda \neq 0$.

Before we start with the proofs of the three claims, we mention that it follows from the definition of $E_\ell$ and our deformation that for $0 \le j \le nk-1$ and $1 \le \ell \le n$ we have

$$ \partial_{z_j} E_\ell \mod t^\alpha = \begin{cases} -\frac{t \epsilon^{k \ell k}}{(j-k(\ell-1))!} u^{m_\ell + j - k(\ell-1) - k}, & \text{if } j \in \{k\ell - k, \dots, k\ell - 1\}, \\ -\frac{t^{\beta+1} \epsilon^k}{h!} u^{m_\ell + h - k} & \text{else, with } h = j \mod k \text{ and } 0 \le h < k. \end{cases} \quad (20) $$

We will crucially use the fact that $\partial_{z_j} E_\ell \mod t^\alpha$ is divisible by $u^{m_\ell-k}$ and that the quotient is a monomial in $u$ of degree at most $k-1$.

**Proof of Claim 1:** Assume that $\det(\tilde{\Lambda}) = 0$. By lemma 10 applied to $\tilde{\Lambda}$, there exist polynomials $\lambda_0, \dots, \lambda_{nk-1} \in \mathbb{L}[t]$ not all equal to 0 such that $\sum_{i=0}^{nk-1} \lambda_i \cdot \tilde{R}_i = 0$. For all $0 \le i \le nk-1$ the matrices associated to the $R_i$'s share the same $n-1$ first columns, while the last column is equal to $(\partial_{z_i} E_\ell)_{1\le\ell\le n}$. From (20) it follows that for all $1 \le i+1, \ell \le n$ and for all $0 \le j \le k-1$, we have $\partial_{z_{ik+j}} E_\ell = \frac{u^j}{j!} \partial_{z_{ik}} E_\ell$. Thus, using the multi-linearity of the determinant, it follows that for all $0 \le i \le n-1$ and $0 \le j \le k-1$, we have $\tilde{R}_{ik+j} = \frac{u^j}{j!} \cdot \tilde{R}_{ik}$. This implies that the linear combination $\sum_{i=0}^{nk-1} \lambda_i \cdot \tilde{R}_i(u) = 0$ rewrites into the form

$$ \sum_{i=0}^{n-1} \sum_{j=0}^{k-1} \lambda_{ik+j} \cdot \frac{u^j}{j!} \cdot \tilde{R}_{ik} = 0. \quad (21) $$

Moreover, the combination of (13) and (20) implies that all entries the $i$th row of the matrix defining $\tilde{R}_j$ are divisible by $u^{m_i-k}$. Define

$$ P_{\ell}(u) := u^{-m_{\ell}+k} \sum_{i=0}^{n-1} \sum_{j=0}^{k-1} \lambda_{ik+j} \cdot \frac{u^j}{j!} \cdot \partial_{z_{ik}} E_{\ell} \mod t^{\alpha} \quad \text{for } \ell = 1, \dots, n, $$

so that equations (13) and (21), as well as the multi-linearity of the determinant imply that

$$ \mathcal{B}(u) := \begin{pmatrix} -u^k + t\epsilon^k & \epsilon^k t^{\beta+1} & \dots & \epsilon^k t^{\beta+1} & P_1 \\ \epsilon^k t^{\beta+1} & -u^k + t\epsilon^k 2^k & \dots & \epsilon^k t^{\beta+1} & P_2 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ \epsilon^k t^{\beta+1} & \epsilon^k t^{\beta+1} & \dots & -u^k + t\epsilon^k (n-1)^k & P_{n-1} \\ \epsilon^k t^{\beta+1} & \epsilon^k t^{\beta+1} & \dots & \epsilon^k t^{\beta+1} & P_n \end{pmatrix} $$

is a singular polynomial matrix. Observe from (20) that $\deg_u P_\ell < k$ for $\ell = 1, \dots, n$. Using this and $\det(\mathcal{B}) = 0$ we shall now show that $P_1 = \dots = P_n = 0$.

First note that when computing the determinant of $\mathcal{B}$ using the Leibniz rule, we would get the product of all entries on the main diagonal and products of other terms which will be polynomials in $u$ of degree at most $(n-1)k-1$. The contribution of the main diagonal contains $u^{(n-1)k}P_n$. Clearly, this term cannot cancel with others, unless $P_n$ vanishes.
Fix some $1 \le \ell \le n-1$ and define $\omega_1, \dots, \omega_k \in \overline{\mathbb{K}}[\epsilon][[t^{1/\ell}]]$ to be the distinct solutions of $u^k = t\epsilon^k \ell^k - t^{\beta+1}\epsilon^k$. Clearly, $\det(\mathcal{B}(\omega_i)) = 0$ for $i = 1, \dots, k$, and by definition of $\omega_i$, the first $n-1$ entries of the $\ell$th row of $\mathcal{B}(\omega_i)$ agree with the first $n-1$ entries of its last row. Therefore, expanding along the last column, we find that

$$0 = \det(\mathcal{B}(\omega_i)) = \pm P_{\ell}(\omega_i) \cdot \det(\mathcal{B}_{\ell,n}(\omega_i)),$$

where $\mathcal{B}_{\ell,n}$ denotes the matrix $\mathcal{B}$ with the $\ell$th row and last column removed. By subtracting the $\ell$th column from all others one can easily compute $\det(\mathcal{B}_{\ell,n}(\omega_i))$, in particular it follows that it is nonzero. So we conclude that $P_{\ell}(\omega_i) = 0$ for $i = 1, \dots, k$ and, then, since $\deg_u P_{\ell}(u) < k$, it finally follows that $P_{\ell}$ vanishes identically.

By (20) we have that $\deg_u \partial_{z_{ik}} E_\ell = m_\ell - k$, therefore, looking at the coefficient of $u^j$ of $P_\ell(u) = 0$ for $j = 0, \dots, k-1$, we obtain with (20) the following linear relations for the $\lambda_i$'s:

$$ \begin{pmatrix} 1 & t^\beta & \cdots & t^\beta \\ t^\beta & 2^k & \cdots & t^\beta \\ \vdots & \vdots & \ddots & \vdots \\ t^\beta & \cdots & \cdots & n^k \end{pmatrix} \times \begin{pmatrix} \lambda_j \\ \lambda_{k+j} \\ \vdots \\ \lambda_{(n-1)k+j} \end{pmatrix} = 0, \quad \text{for all } 0 \le j \le k-1. \qquad (22) $$

As the square matrix in (22) is invertible, we obtain that $\lambda_i = 0$ for all $i = 0, \dots, nk-1$ which contradicts our assumption. Hence we have proved that $\det(\tilde{\Lambda}) \neq 0$.

**Proof of Claim 2:** We prove the following explicit factorization of $\det(\tilde{\Lambda})$:

$$ \det(\tilde{\Lambda}) = \prod_{i=1}^{nk} u_i^{M-nk} \cdot \prod_{i<j} (u_i - u_j) \cdot H(t), \qquad (23) $$

where $H(t) \in \mathbb{K}[\epsilon, t]$ is nonzero and of degree in $t$ at most $n^2(\beta + 1)$.

From (13) we know that for $0 \le j \le nk-1$, the polynomial $\tilde{R}_j$ is given by

$$ \det \begin{pmatrix} u^{m_1} + t\epsilon^k u^{m_1-k} & t^{b+1} u^{m_1-k} & \cdots & \epsilon^k t^{b+1} u^{m_1-k} & \partial_{z_j} E_1 \\ \epsilon^k t^{b+1} u^{m_2-k} & u^{m_2} + t\epsilon^k 2^k u^{m_2-k} & \cdots & \epsilon^k t^{b+1} u^{m_2-k} & \partial_{z_j} E_2 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ \epsilon^k t^{b+1} u^{m_{n-1}-k} & \epsilon^k t^{b+1} u^{m_{n-1}-k} & \cdots & u^{m_{n-1}-k} (u^k + t\epsilon^k (n-1)^k) & \partial_{z_j} E_{n-1} \\ \epsilon^k t^{b+1} u^{m_n-k} & \epsilon^k t^{b+1} u^{m_n-k} & \cdots & \epsilon^k t^{b+1} u^{m_n-k} & \partial_{z_j} E_n \end{pmatrix} \mod t^\alpha. $$

It follows from (20) that $u^{-m_\ell+k} \deg_{z_j} E_\ell$ is a monomial of degree at most $k-1$ in $u$. Since all other terms in the $\ell$th row of $\tilde{R}_j$ are trivially divisible by $u^{m_\ell-k}$, and $\alpha > n(\beta+1)$, it follows that $\tilde{R}_j$ factors as $u^{M-nk}$ times the determinant of the polynomial matrix

$$ \begin{pmatrix} u^k + t\epsilon^k & \epsilon^k t^{b+1} & \cdots & \epsilon^k t^{b+1} & u^{-m_1+k} \partial_{z_j} E_1 \\ \epsilon^k t^{b+1} & u^k + t\epsilon^k 2^k & \cdots & \epsilon^k t^{b+1} & u^{-m_2+k} \partial_{z_j} E_2 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ \epsilon^k t^{b+1} & \epsilon^k t^{b+1} & \cdots & u^k + t\epsilon^k (n-1)^k & u^{-m_{n-1}+k} \partial_{z_j} E_{n-1} \\ \epsilon^k t^{b+1} & \epsilon^k t^{b+1} & \cdots & \epsilon^k t^{b+1} & u^{-m_n+k} \partial_{z_j} E_n \end{pmatrix} \mod t^\alpha, \qquad (24) $$

where, as before, $M := \sum_{i=1}^n m_i$. Thus, $u^{M-nk}$ divides $\tilde{R}_j$. Moreover, the degree in $u$ of the first $n-1$ columns of (24) is upper-bounded by $k$ and the last column of this matrix has degree at most $k-1$, so it follows that $\deg_u(\tilde{R}_j) \le M-1$.
Let us consider $\tilde{\Lambda}$ which is by definition the matrix $(\tilde{R}_j^{(i)})_{1\le i,j+1\le nk} \in \mathbb{L}[t][u_1, \dots, u_{nk}]^{nk \times nk}$. Since $\deg_u(\tilde{R}_j) \le M-1$, we also have that $\det_{u_i}(\tilde{\Lambda}) \le M-1$ for each $i=1, \dots, nk$. From the considerations above it also follows that $\det(\tilde{\Lambda})$ is divisible by $\prod_{i=1}^{nk} u_i^{M-nk}$. Moreover, it is obvious that if $u_i=u_j$ for some $i \ne j$, the determinant of $\tilde{\Lambda}$ vanishes. Hence, we can also factor out the Vandermonde determinant $\prod_{i<j} (u_i-u_j)$ from $\det(\tilde{\Lambda})$. By comparing degrees in each $u_i$ it follows that the remaining factor is a constant with respect to $u_i$ and hence (23) holds. By Claim 1, $H(t)$ is nonzero, so it remains to prove that $\deg_t H(t) \le n^2k(\beta+1)$. This follows directly from the fact that the degree in $t$ of the $\tilde{R}_j$'s is bounded by $n(\beta+1)$ as determinants of $(n \times n)$ matrices whose polynomial coefficients have degree in $t$ upper-bounded by $\beta+1$.

Finally for this step, using the definitions $\alpha = 3n^2k \cdot (\beta+1) + 3nM$ and $\beta = \lfloor 2M/k \rfloor$, we have that $\alpha > n^2k(\beta+1)$. So that $\deg_t(\det(\tilde{\Lambda})) < \alpha$: by Claim 1, we obtain $\det(\tilde{\Lambda}) \ne 0 \mod t^\alpha$.

**Proof of Claim 3:** Having established (23) it is now enough to show that

$$ \mathrm{val}_t \left( \left( \prod_{i=1}^{nk} U_i^{M-nk} \cdot \prod_{i<j} (U_i - U_j) \cdot H(t) \right) < \alpha. \right) \quad (25) $$

Recall from lemma 6 that $U_i = t^{1/k} \epsilon_j \zeta^\ell + (\text{higher powers of } t)$ for $i=1, \dots, nk, j=1, \dots, n$ and $\zeta$ a primitive $k$th root of unity. Thus, and because $\deg_t(H) \le n^2k(\beta+1)$ by Claim 2, the left-hand side of (25) is bounded by $nk \cdot \frac{M-nk}{k} + nk(nk-1) \cdot \frac{1}{k} + n^2k(\beta+1)$, which is at most $nM+n^2k(\beta+1)$. Since $\alpha = 3n^2k \cdot (\beta+1) + 3nM$, we conclude that $\det(\tilde{\Lambda}(U_1, \dots, U_{nk})) \mod t^\alpha \ne 0$. Finally, together with (19) this implies that $\det(\tilde{\Lambda})$ does not vanish as well. $\square$

Having now proved that $\det(\mathrm{Jac}_{S_{\mathrm{dup}}^\epsilon}) \ne 0$ at $\mathcal{P}$, we can apply lemma 7 and obtain that the specialized series $G_i(0), \dots, \partial_u^{k-1}G_i(0)$ are all algebraic over $\mathbb{K}(t, \epsilon)$. The algebraicity of the complete formal power series $G_1, \dots, G_n$ over $\mathbb{K}(t, u, \epsilon)$ then follows again by [8, Lemma 2.10] from the invertibility of the Jacobian matrix of $E_1, \dots, E_n$ considered with respect to the variables $x_1, \dots, x_n$ (with $t, u, z_0, \dots, z_{nk-1}$ viewed as parameters). The equalities $F_i(t^\alpha, u) = G_i(t, u, 0)$ finally imply that $F_1, \dots, F_n$ are also algebraic over $\mathbb{K}(t, u)$.

## 3.2 Proof of theorem 4

We prove the quantitative estimates announced in theorem 4. The techniques of the proof being standard in effective algebraic geometry, we will be quite brief in the arguments. For more details on the upper bound on the algebraicity degree, we refer the reader to [8, Prop. 2.8] in the case $n=k=1$, and to [9, Prop. 3] for the more general case $n=1$ ($k$ arbitrary). For the complexity proof, we refer the reader to [8, Prop. 2.9] and [9, Prop. 4] in the two respective cases. As in the proof of theorem 3, we assume without loss of generality that $a=0$.

**Algebraicity degree bound:** In order to bound the algebraicity degree of $F_i(t, u)$, we will provide an upper bound on the algebraicity degree of $G_i(t, u, \epsilon)$ over $\mathbb{K}(t, u, \epsilon)$, where, as in the proof of theorem 3, $(G_1, \dots, G_n)$ denotes the solution of the deformed system (7). Following the lines of the proof of Theorem 3, we shall first give explicit bounds on the total degrees of all equations in $S_{\mathrm{dup}}^\epsilon$. Note that in this part we are interested in the algebraicity degree of $G_i(t, 0)$ over $\mathbb{K}(t, \epsilon)$, so in the computation of the total degree we take into account all the variables
$u, u_1, \dots, u_{nk}, x_1, \dots, x_{n^2k}, z_0, \dots, z_{nk-1}$ but not $t$ and $\epsilon$. In order to distinguish this restricted notion of total degree from the usual total degree denoted by totdeg, we will write it totdeg$_{t,\epsilon}$.

Let $\delta$ be a bound on totdeg$_{t,\epsilon}$ of $f_1, \dots, f_n, Q_1, \dots, Q_n$; then the total degrees of $E_1, \dots, E_n$ are bounded by $\delta(k+1)$. Moreover, totdeg$_{t,\epsilon}$(Det), totdeg$_{t,\epsilon}$(P) $\le n\delta(k+1)$, since Det and P are determinants of polynomial matrices of size $n \times n$ and degree at most $\delta(k+1)$.

Recall from lemma 7 that the ideal $\mathcal{I}_{\text{dup}}^{\epsilon} \subseteq \mathbb{K}(t, \epsilon)[x_1, \dots, x_{n^2k}, u_1, \dots, u_{nk}, z_0, \dots, z_{nk-1}]$ is defined by $nk$ duplications of the polynomials $E_1, \dots, E_n$, Det, P. By lemma 7 the saturated ideal $\mathcal{I}_{\text{dup}}^{\epsilon} : \det(\text{Jac}_{\mathcal{S}_{\text{dup}}^{\epsilon}})^{\infty}$ is radical and of dimension 0, and this yields bounds for the algebraicity degrees of the specializations $\{\partial_u^i G_j(0)\}_{0 \le i \le k-1, 1 \le j \le n}$. More precisely, applying the Heintz-Bézout theorem [20, Theorem 1] in the same way as in [8, Prop. 3.1], one obtains that the algebraicity degree of any $\partial_u^i G_j(0)$ is bounded by

$$ \max_{i=1,\dots,n} (\text{totdeg}_{t,\epsilon}(E_i))^{n^2k} \cdot \text{totdeg}_{t,\epsilon}(\text{Det})^{nk} \cdot \text{totdeg}_{t,\epsilon}(P)^{nk} \le n^{2nk} \cdot (\delta(k+1))^{nk(n+2)}. $$

We denote this bound by $\gamma(n, k, \delta)$.

As in [9, Prop. 4] there is a group action of the symmetric group $S_{nk}$ on the zero set associated to the ideal $\langle E_1^{(1)}, \dots, E_1^{(nk)}, \dots, E_n^{(1)}, \dots, E_n^{(nk)}, \text{Det}^{(1)}, \dots, \text{Det}^{(nk)}, P^{(1)}, \dots, P^{(nk)} \rangle$ which permutates each of the $nk$ duplicated blocks of coordinates. As this action preserves the $\{t, z_0, \dots, z_{nk-1}\}$-coordinate space, one spares the cardinality of the orbits and deduces that an algebraicity upper bound on any of the specialized series $\partial_u^i G_j(0)$ is also given by $\gamma(n, k, \delta)/(nk)!$, that is we have

$$ \mathrm{algdeg}_{\mathbb{K}(t,\epsilon)}(\partial_u^i G_j(0)) \le n^{nk} \cdot (\delta(k+1))^{nk(n+2)}/(nk)! \quad (26) $$

Finally, it remains to give algebraicity bounds for the full series $G_1, \dots, G_n$. This is done in the same fashion as above, this time over the field extension

$$ \mathbb{K}(t, \epsilon, G_1(0), \dots, \partial_u^{k-1} G_1(0), \dots, G_n(0), \dots, \partial_u^{k-1} G_n(0))/\mathbb{K}(t, \epsilon). $$

Using (26) and the multiplicativity of the degrees in field extensions yields

$$ \mathrm{algdeg}_{\mathbb{K}(t,u,\epsilon)}(G_i(u)) \le \delta^n (k+1)^n \cdot \left( \frac{n^{2nk} (\delta(k+1))^{nk(n+2)}}{(nk)!} \right)^{nk} \quad \text{for } i=1, \dots, n, $$

and, by specialization, the same bound then holds for $\mathrm{algdeg}_{\mathbb{K}(t,u)}(F_i(u))$ as well.

**Complexity estimate:** Let us recall that we fix $\mathbb{K}$ to be an effective field of characteristic 0 and we count the number of elementary operations $(+, -, \times, \div)$ in $\mathbb{K}$. For proving the arithmetic complexity estimate in theorem 4, we rely on the use of the parametric geometric resolution (see [24] for details). For this purpose we need:

* the number of parameters and of variables of the input equations of interest,

* and upper bound on the degree (with respect to all the variables and parameters) of the saturated ideal $\mathcal{I}_{\text{dup}}^{\epsilon} : \det(\text{Jac}_{\mathcal{S}_{\text{dup}}^{\epsilon}})^{\infty}$.

In our setting the parameters are $t$ and $\epsilon$ so we have two of them, and the variables are

$$ x_1, \dots, x_{n^2k}, z_0, \dots, z_{nk-1}, u_1, \dots, u_{nk}, $$
so $nk(n + 2)$ in total.

For computing an upper bound $\nu(n, k, \delta)$ on the degree of the ideal $\mathcal{I}_{\text{dup}}^{\epsilon} : \det(\text{Jac}_{\mathcal{S}_{\text{dup}}^{\epsilon}})^{\infty}$, we again apply the Heintz-Bézout theorem. Note that this time we must take into account the total degree of the input system $\mathcal{S}_{\text{dup}}^{\epsilon}$ with respect to both the variables and the parameters. By the same arguments as in the first part of the proof and using the definitions of $\alpha, \beta$, we find

$$ \nu(n, k, \delta) = n^{2nk}(\delta \cdot (k+1) + 3n^2k \cdot (2n\delta + 1) + 3nM)^{nk(n+2)}. $$

Without loss of generality, we only address the computation of an annihilating polynomial of $G_1(0)$ (and consequently of $F_1(0)$). The algorithm we propose consists of two steps:

**Step 1:** For $\lambda$ a new variable which is a random $\mathbb{K}[t, \epsilon]$-linear combination of all the variables involved in $\mathcal{S}_{\text{dup}}^{\epsilon}$, we compute two polynomials $V(t, \epsilon, \lambda)$, $W(t, \epsilon, \lambda) \in \mathbb{K}[t, \epsilon][\lambda]$, such that: $G_1(0)$ is a solution<fnref n="6" /> of $\mathcal{S}_{\text{dup}}^{\epsilon}$ and not a solution of $\det(\text{Jac}_{\mathcal{S}_{\text{dup}}^{\epsilon}}) = 0$ if and only if $G_1(0) = V(t, \epsilon, \lambda_0)/\partial_{\lambda}W(t, \epsilon, \lambda_0)$ for $\lambda_0 \in \overline{\mathbb{K}(t, \epsilon)}$ a solution of $W(t, \epsilon, \lambda) = 0$.

**Step 2:** Compute the squarefree part $R \in \mathbb{K}[t, \epsilon, z_0]$ of the resultant of $z_0 \cdot \partial_{\lambda}W(t, \epsilon, \lambda) - V(t, \epsilon, \lambda)$ and $W(t, \epsilon, \lambda)$ with respect to $\lambda$.

By definition, $\lambda_0$ is a solution of $\partial_{\lambda}W(t, \epsilon, \lambda)G_1(0) - V(t, \epsilon, \lambda) = 0$ and $W(t, \epsilon, \lambda) = 0$, so $R(t, \epsilon, G_1(0)) = 0$ by the property of the resultant. Moreover, it follows from the eigenvalue theorem [17, Theorem 1] and the fact that the dimension of the radical ideal $\mathcal{I}_{\text{dup}}^{\epsilon} : \det(\text{Jac}_{\mathcal{S}_{\text{dup}}^{\epsilon}})^{\infty}$ is zero, that $R$ is a nonzero polynomial.

By [24, Theorem 2], performing Step 1 can be done using $\tilde{O}((Ln^2k + (n^2k)^4) \cdot \nu(n, k, \delta)^3)$ operations in $\mathbb{K}$, where $L$ denotes the complexity of evaluating the duplicated system $\mathcal{S}_{\text{dup}}^{\epsilon}$ at all variables and parameters. Using the Bauer-Strassen theorem [3, Theorem 1], the cost for evaluating Det and $P$ is in $O(nL')$, where $L'$ is the cost for evaluating any of the $E_1, \dots, E_n$. The total degrees of the $E_i$'s are bounded by $d' := \delta \cdot (k+1) + 3n^2k \cdot (2n\delta+1) + 3nM$, so $L' \in O((d')^{nk+n+3})$, since $nk + n + 3$ is the number of variables and parameters. Duplicating $nk$ times the initial polynomials $E_1, \dots, E_n$, Det, $P$, it follows that $L \in O(n^3kL') \subset O(n^3k(d')^{nk+n+3})$. To do Step 2, we perform evaluation-interpolation on $t$ and $\epsilon$. This requires to use in total $O(\nu(n, k, \delta)^2)$ distinct points from the base field $\mathbb{K}$. For each specialization of $z_0 \cdot \partial_{\lambda}W(t, \epsilon, \lambda) - V(t, \epsilon, \lambda)$ and $W(t, \epsilon, \lambda)$ to $t = \theta_1 \in \mathbb{K}$ and $\epsilon = \theta_2 \in \mathbb{K}$, we compute the squarefree part of the bivariate resultant of $\partial_{\lambda}W(\theta_1, \theta_2, \lambda) \cdot z_\ell - V(\theta_1, \theta_2, \lambda)$ and $W(\theta_1, \theta_2, \lambda)$ with respect to $\lambda$. Using [26], this requires $\tilde{O}((\nu(n, k, \delta)/(nk)!)^{3.63})$ operations in $\mathbb{K}$.

In total we obtain that the arithmetic complexity is bounded by $(2nk\delta)^{O(n^2k)}$ operations in $\mathbb{K}$. Deducing an annihilating polynomial of $F_1(t, 0)$ from the one annihilating $G_1(t, \epsilon = 0, u = 0)$ is included in the previous complexities. This proves theorem 4. □

# 4 Elementary strategies for solving a system of DDEs

The present section aims at studying two strategies available at this time for solving a system of DDEs in practice. Both approaches have the advantage to be much more efficient than the algorithm resulting from theorem 3, however they are not guaranteed to always terminate. As before, the goal is, given a system of DDEs of the form (5), to compute a nonzero polynomial $R \in \mathbb{K}[t, z_0]$ such that $R(t, F_1(t, a)) = 0$.

<footnote n="6">⁶More precisely: the $z_0$th coordinate of a solution of $\mathcal{S}_{\text{dup}}^{\epsilon}$ which is not a solution of $\det(\text{Jac}_{\mathcal{S}_{\text{dup}}^{\epsilon}}) = 0$.</footnote>
First, in section 4.1, we summarize the algorithm underlying the duplication strategy that was used in the proof of theorem 3 and we analyze its arithmetic complexity in the case where no symbolic deformation is necessary (the arithmetic cost for solving a system of DDEs when such a symbolic deformation is needed was the content of theorem 4). In section 4.2, we study the practical approach introduced and used in [10, Section 11]: this approach consists in reducing the algorithmic study of a system of DDEs to the study of a single functional equation (which may not be of a fixed-point type). Once we are left with this single equation, we apply the geometry-driven algorithm from [9, Section 5] which avoids the duplication of the variables.

# 4.1 The classical duplication of variables algorithm

## 4.1.1 General strategy

Consider the system of DDEs

$$ \left\{ \begin{array}{l} (\mathbf{E}_{\mathbf{F}_1}): \quad F_1 = f_1(u) + t \cdot Q_1(\nabla_a^k F_1, \dots, \nabla_a^k F_n, t, u), \\ \vdots \\ (\mathbf{E}_{\mathbf{F}_n}): \quad F_n = f_n(u) + t \cdot Q_n(\nabla_a^k F_1, \dots, \nabla_a^k F_n, t, u), \end{array} \right. \qquad (27) $$

introduced in theorem 3 and denote by $S_{\text{dup}}$ the duplicated polynomial system obtained by

* considering the polynomials $E_1, \dots, E_n \in \mathbb{K}(t)[x_1, \dots, x_n, u, z_0, \dots, z_{nk-1}]$ associated to the "numerator equations"<fnref n="7" /> of the system of DDEs (27), as well as the polynomials

$$ \text{Det} := \det \begin{pmatrix} \partial_{x_1} E_1 & \dots & \partial_{x_n} E_1 \\ \vdots & \ddots & \vdots \\ \partial_{x_1} E_n & \dots & \partial_{x_n} E_n \end{pmatrix} \quad \text{and} \quad P := \det \begin{pmatrix} \partial_{x_1} E_1 & \dots & \partial_{x_{n-1}} E_1 & \partial_u E_1 \\ \vdots & \ddots & \vdots & \vdots \\ \partial_{x_1} E_n & \dots & \partial_{x_{n-1}} E_n & \partial_u E_n \end{pmatrix}, $$

* duplicating $nk$ times the variables $x_1, \dots, x_n, u$ in the polynomials $(E_1, \dots, E_n, \text{Det}, P)$ in order to obtain $nk$ duplications of the polynomials $(E_1, \dots, E_n, \text{Det}, P)$ in the new variables $x_1, \dots, x_{n^2k}, u_1, \dots, u_{nk}$ and (unchanged variables) $z_0, \dots, z_{nk-1}$.

**Example 2 (cont.).** We consider example 2, where $k = 1$ and $n = 2$. We have the polynomials

$$ \begin{cases} E_1 := -(x_1 - 1)(u - 1) + tu(2ux_1^2 - uz_0 + 2uz_2 - 2x_1^2 + u + x_1 - 2z_2 - 1), \\ E_2 := -x_2(u - 1) + tu(2ux_1x_2 + ux_1 - 2x_1x_2 - x_1 + x_2 - z_2). \end{cases} $$

We compute

$$ \text{Det} = \det \begin{pmatrix} \partial_{x_1} E_1 & \partial_{x_2} E_1 \\ \partial_{x_1} E_2 & \partial_{x_2} E_2 \end{pmatrix} = (4tu^2x_1 - 4tux_1 + tu - u + 1)(2tu^2x_1 - 2tux_1 + tu - u + 1). $$

and the rather large polynomial $P = \det \begin{pmatrix} \partial_{x_1} E_1 & \partial_u E_1 \\ \partial_{x_1} E_2 & \partial_u E_2 \end{pmatrix}$. We obtain by the duplication strategy the set $S_{\text{dup}}$ of polynomials given by

$$ S_{\text{dup}} := \{ (E_1(x_1, x_2, u_1, z_0, z_1), E_2(x_1, x_2, u_1, z_0, z_1), \text{Det}(x_1, x_2, u_1, z_0, z_1), P(x_1, x_2, u_1, z_0, z_1), \\ E_1(x_3, x_4, u_2, z_0, z_1), E_2(x_3, x_4, u_2, z_0, z_1), \text{Det}(x_3, x_4, u_2, z_0, z_1), P(x_3, x_4, u_2, z_0, z_1) \}. $$

<footnote n="7">⁷For all $1 \le i \le n$, $E_i(u) = E_i(F_1, \dots, F_n, u, F_1(t, a), \dots, \partial_u^{k-1} F_1(t, a), \dots, F_n(t, a), \dots, \partial_u^{k-1} F_n(t, a))$.</footnote>
Observe that $S_{\text{dup}}$ is built from 8 equations in the 8 unknowns $x_1, x_2, x_3, x_4, u_1, u_2, z_0, z_1$ and the parameter $t$.

Assuming that we did not deform (27) using (7) (and thus that we did not introduce any deformation parameter $\epsilon$), it follows by construction that $S_{\text{dup}}$ is defined in the polynomial ring $\mathbb{K}(t)[x_1, \dots, x_{n^2k}, u_1, \dots, u_{nk}, z_0, \dots, z_{nk-1}]$, and is built from $nk(n+2)$ equations in $nk(n+2)$ unknowns ($t$ being considered as a parameter). We introduce the polynomial sat $:= t \prod_{j=1}^{nk} (\prod_{i=1, i \neq j}^{nk} (u_i - u_j))(u_j - a) \in \mathbb{K}(t)[u_1, \dots, u_{nk}]$ and define the saturated ideal$^8 \mathcal{I}_{\text{dup}}^{\infty} := \langle S_{\text{dup}} \rangle : \text{sat}^{\infty}$ in the ring $\mathbb{K}(t)[x_1, \dots, x_{n^2k}, u_1, \dots, u_{nk}, z_0, \dots, z_{nk-1}]$. In (H1) below we shall assume that $\mathcal{I}_{\text{dup}}^{\infty}$ has expected dimension 0 over $\mathbb{K}(t)$, without any multiplicity points, and we also assume that there exist $nk$ distinct solutions to $\text{Det}(u) = 0$ in $u$, all of them lying in $\overline{\mathbb{K}}[[t^\frac{1}{\epsilon}]] \setminus \overline{\mathbb{K}}$. All systems of DDEs that can be solved with the duplication approach satisfy this assumption, up to removing the multiplicity points (which is harmless for our objective of solving systems of DDEs).

**Hypothesis 1:** • $\text{Det}(u) = 0$ admits $nk$ distinct solutions (in $u$) in $\overline{\mathbb{K}}[[t^\frac{1}{\epsilon}]] \setminus \overline{\mathbb{K}}$. (H1)
• The ideal $\mathcal{I}_{\text{dup}}^{\infty}$ is radical of dimension 0 over $\mathbb{K}(t)$.

Under (H1), we denote by $U_1, \dots, U_{nk} \in \overline{\mathbb{K}}[[t^\frac{1}{\epsilon}]] \setminus \overline{\mathbb{K}}$ the distinct solutions to $\text{Det}(u) = 0$.

Recall that after duplication of the variables $x_1, \dots, x_n, u$ in the initial set of polynomials $(E_1, \dots, E_n, \text{Det}, P)$, we took the convention in section 3 to have for all $1 \le i \le n$, $1 \le j \le nk$ and all $0 \le \ell \le k - 1$, the correspondence “variables $\leftrightarrow$ values” as follows:

$$x_{(j-1)n+i} \leftrightarrow F_i(t, U_j(t)), \quad u_j \leftrightarrow U_j(t), \quad z_{(i-1)k+\ell} \leftrightarrow (\partial_u^\ell F_i)(t, a).$$

It follows from the second part of (H1) that the elimination ideal $\mathcal{I}_{\text{dup}}^{\infty} \cap \mathbb{K}[t, z_0]$ is not reduced to $\{0\}$. In spirit of [9, Proposition 2] and [8, Proposition 2.1], we can show that any nonzero element of $\mathcal{I}_{\text{dup}}^{\infty} \cap \mathbb{K}[t, z_0]$ annihilates the series $F_1(t, a)$.

**Lemma 11.** Under (H1), any element $R \in \mathcal{I}_{\text{dup}}^{\infty} \cap \mathbb{K}[t, z_0] \setminus \{0\}$ satisfies $R(t, F_1(t, a)) = 0$.

*Proof.* Using $R \in \mathcal{I}_{\text{dup}}^{\infty}$ and the definition of the saturation, it follows that there exists some $\ell \in \mathbb{N}$ such that $R \cdot (\prod_{i \neq j} (u_i - u_j)(u_i - a)t)^\ell$ can be written as a sum of multiples of the duplicated polynomials obtained from $E_1, \dots, E_n, \text{Det}, P$. Specializing the resulting expression to the values $x_{(j-1)n+i} = F_i(t, U_j(t))$, $u_j = U_j(t)$, $z_{(i-1)k+\ell} = (\partial_u^\ell F_i)(t, a)$ and using (H1) implies that $U_i(t) \neq U_j(t)$ and $U_i(t) \neq a$ (whenever $i \neq j$). This implies that $R(t, F_1(t, a)) = 0$. $\square$

**Example 2 (cont.).** Consider $S_{\text{dup}}$, and compute a nonzero generator $R \in \mathbb{Q}[t, z_0]$ of $\langle S_{\text{dup}}, m \cdot (u_1 - u_2) \cdot (u_1 - 1) \cdot (u_2 - 1)t - 1 \rangle \cap \mathbb{Q}[t, z_0]$. For this computation, which is not easy to perform with a naive approach, we rely on the recent work by the first author [21] which provides a Maple package with efficient implementations of algorithms for solving DDEs: one of these computes efficiently such an elimination polynomial. We obtain in a few seconds on a regular laptop that

$$R = (64t^3 z_0^3 + 2t(24t^2 - 36t + 1)z_0^2 + (-15t^3 + 9t^2 + 19t - 1)z_0 + t^3 + 27t^2 - 19t + 1) \\ \cdot (z_0 - 1) \cdot (2tz_0 + t - 1) \cdot (36 - 60z_0 + t \cdot \tilde{R}(t, z_0))$$

annihilates $F_1(t, 1)$; here we write $\tilde{R}(t, z_0) \in \mathbb{Q}[t, z_0]$ for an explicit but rather big polynomial. Moreover, as $F_1 = 1 + O(t)$ and $F_1$ is not a constant, we can refine our conclusion and identify that the first factor of $R$ is the minimal polynomial of $F_1(t, 1)$.

<footnote>$^8$Recall that in practice, a generating set of the ideal $\mathcal{I}_{\text{dup}}^{\infty}$ can be obtained by computing a Gröbner basis of $\langle S_{\text{dup}}, m \cdot \text{sat} - 1 \rangle \cap \mathbb{K}(t)[x_1, \dots, x_{n^2k}, u_1, \dots, u_{nk}, z_0, \dots, z_{nk-1}]$, where $m$ is an extra variable introduced in order to remove the solution set of the equation sat = 0.</footnote>
4.1.2 Refined complexity and size estimates

We prove a version of theorem 4 in the case where the symbolic deformation is not needed. In this context, we provide an upper bound on the algebraicity degree of $F_1(t, a)$ over $\mathbb{K}(t)$. Also, we estimate the arithmetic complexity for the computation of a nonzero element of $\mathcal{I}_{\text{dup}}^{\infty} \cap \mathbb{K}[t, z_0]$.

**Proposition 12.** Let $\mathbb{K}$ be a field of characteristic 0. Consider the system of DDEs (27), assume that (H1) holds and denote by $\delta$ an upper bound on the total degrees of $f_1, \dots, f_n, Q_1, \dots, Q_n$. Then the algebraicity degree of $F_1(t, a)$ over $\mathbb{K}(t)$ is bounded by $n^{2nk}(\delta(k+1)+1)^{nk(n+2)}/(nk)!$. Moreover when $\mathbb{K}$ is effective, there exists an algorithm computing a nonzero polynomial $R \in \mathbb{K}[t, z_0]$ such that $R(t, F_1(t, a)) = 0$, in $\tilde{O}((nk\delta)^{5.26n(nk+k+1)})$ operations in $\mathbb{K}$.

*Proof.* In spirit of the proof of theorem 4:

* Step 1: We compute an upper bound on the degree of the ideal $\mathcal{I}_{\text{dup}}^{\infty}$ and then deduce an upper bound on the algebraicity degree of $F_1(t, a)$.

* Step 2: We make explicit the number of variables and parameters in $S_{\text{dup}}$.

* Step 3: We describe and apply an algorithm which is itself based on the algorithm of [24, Theorem 2] in order to compute a nonzero annihilating polynomial of $F_1(t, a)$.

* Step 4: Finally, we use the arithmetic complexity of the algorithm on which [24, Theorem 2] relies in order to estimate the complexity of our algorithm.

In order to define the polynomials $E_1, \dots, E_n \in \mathbb{K}(t)[x_1, \dots, x_n, u, z_0, \dots, z_{nk-1}]$ associated to the numerator equations of (27), recall that it is necessary to multiply the $i$th DDE in (27) by a power $m_i \in \mathbb{N}$ of $(u-a)$. The polynomials $Q_i$ being of total degree upper-bounded by $\delta$, and the DDEs considered being of order $k$, it results that each $m_i$ is upper-bounded by $k\delta$. So each of the $E_i$'s has total degree upper-bounded by $d := \delta(k+1)+1$. Thus $\text{Det}$ and $P$ have their respective total degrees bounded by $nd$ as determinants of $(n \times n)$-matrices whose entries are polynomials of total degrees bounded by $d$. It remains to see that because of the radicality assumption in (H1), the Heintz-Bézout theorem [20, Theorem 1] applies and implies that the degree of the ideal $\mathcal{I}_{\text{dup}}^{\infty}$ is bounded by the product of the total degrees of the $nk(n+2)$ duplications $E_1^{(1)}, \dots, E_n^{(1)}, \text{Det}^{(1)}, P^{(1)}, \dots, E_1^{(nk)}, \dots, E_n^{(nk)}, \text{Det}^{(nk)}, P^{(nk)}$. Such a bound is given by $\gamma(n, k, \delta) := n^{2nk}d^{nk(n+2)}$. Similarly to the proof of theorem 4, there is a group action of the symmetric group $S_k$ over the zero set of $\mathcal{I}_{\text{dup}}^{\infty}$ in $\overline{\mathbb{K}(t)}^{nk(n+2)}$, obtained by permuting the duplicated blocks of coordinates, and which preserves the $(z_0, \dots, z_{nk-1})$-coordinate space. This group action implies that the algebraicity degree of $F_1(t, a)$ is upper-bounded by $\gamma(n, k, \delta)/(nk)! = n^{2nk}(\delta(k+1)+1)^{nk(n+2)}/(nk)!$.

Recall that we have one parameter $t$, and $nk(n+2)$ variables for the $x_i$'s, the $u_i$'s and the $z_i$'s. We will make explicit some $R \in \mathbb{K}[t, z_0] \setminus \{0\}$ annihilating $F_1(t, a)$. To do so, we first compute (here we use the second part of (H1) again) two polynomials $V, W \in \mathbb{K}[t, \lambda]$ such that: (i) $W$ is squarefree, (ii) $\lambda$ is a new variable which is a $\mathbb{K}[t]$-linear combination of the $nk(n+2)$ variables in Step 2, (iii) for all the zeros $\alpha \in \overline{\mathbb{K}(t)}^{nk(n+2)}$ of $S_{\text{dup}}$ that are not solutions of sat, there exists $\lambda_0 \in \overline{\mathbb{K}(t, \epsilon)}$ solution in $\lambda$ of $W(t, \lambda) = 0$ such that $V(t, \lambda_0)/\partial_{\lambda}W(t, \lambda_0)$ is the $z_0$-coordinate of $\alpha$. Using these polynomials, it is straightforward to see by applying Stickelberger's theorem [17, Theorem 1], using the radicality assumption in (H1) and then applying lemma 11 that the squarefree part
of the resultant of $z_0 \cdot \partial_\lambda W(t, \lambda) - V(t, \lambda)$ and $W(t, \lambda)$ with respect to $\lambda$ is a nonzero polynomial of $K[t, z_0]$ annihilating the series $F_1(t, a)$.

It remains to estimate the complexity of Step 3. Following the proof of theorem 4, the application of the algorithm underlying [24, Theorem 2] allows us to compute $V$ and $W$ in $\tilde{O}((n^2kL + (n^2k)^4)\gamma(n, k, \delta)^2)$ operations in $K$, where $L$ is the length of a straight-line program which evaluates the system $S_{\text{dup}}$ and the polynomial sat. Since the cost for evaluating $E_1, \dots, E_n$ is included in $O(d^{nk+n+2})$, then by the Baur-Strassen's theorem [3, Theorem 1], the complexity of evaluating the polynomials Det, $P$ is included in $O(n^2k\delta^{nk+n+2})$. Thus, evaluating $E_1, \dots, E_n$, Det, $P$ has an arithmetic cost which is in $O(n(\delta k)^{nk+n+2})$. As we considered $nk$ duplication of the polynomials $(E_1, \dots, E_n, \text{Det}, P)$, we find $L \in O(n^3kd^{nk+n+2}) \subset O(n^2k(\delta k)^{nk+n+2})$. It remains to compute the squarefree part $R \in K[t, z_0] \setminus \{0\}$ of the resultant of $z_0 \cdot \partial_\lambda W(t, \lambda) - V(t, \lambda)$ and $W(t, \lambda)$ with respect to $\lambda$. Note that under (H1), the quantity $\gamma(n, k, \delta)/(nk)!$ bounds the partial degrees of $R$. This resultant computation can be done by applying evaluation-interpolation with respect to $t$ with $O(n^{2nk}(\delta k)^{nk(n+2)}/(nk)!)$ points. For each specialization at say $t = \theta \in K$, we apply for instance [26] and deduce that the bivariate resultant computation of $z_0 \cdot \partial_\lambda W(\theta, \lambda) - V(\theta, \lambda)$ and $W(\theta, \lambda)$ with respect to $\lambda$ can be done in

$$ \tilde{O}(\gamma(n, k, \delta)^{2.63}/(nk)!^{2.63}) \subset \tilde{O}(n^{5.26nk}(\delta k)^{2.63nk(n+2)}/(nk)!^{2.63}) $$

operations in $K$. Replacing $L$ and $\gamma(n, k, \delta)$ by their values and summing up all arithmetic costs, one deduces that the arithmetic cost for computing $R$ is included in $\tilde{O}((nk\delta)^{5.26n(nk+k+1)})$. $\square$

**Remark 13.** (i) Note that the complexity in proposition 12 matches, as expected, the one that was proven in [9, Proposition 4] in the case $n = 1$. One shall however see that the total degree of $E_1$ is considered in [9, Proposition 4], whereas here we considered the total degree of $Q_1$ and $f_1$. Passing from the complexity in proposition 12 to the one in [9, Proposition 4] is done by replacing $\delta$ by $\delta/k$.

(ii) It follows from the inclusion $\tilde{O}((nk\delta)^{5.26n(nk+k+1)}) \subset \tilde{O}((nk\delta)^{40(n^2k+1)})$ that the arithmetic complexity stated in proposition 12 refines the arithmetic complexity stated in theorem 4.

As shown in section 2 and section 4.1, a natural way of solving (27) is to compute a nonzero element $R \in \mathcal{I}_{\text{dup}}^\infty \cap K[t, z_0]$. However, as already pointed out in [9, Section 3] in the case $n = 1$, this process of duplicating variables yields an exponential growth of the degree of the ideal $\mathcal{I}_{\text{dup}}^\infty$ with respect to the number $nk$ of duplications [20, Proposition 2]. It is known that the degree of the ideal generated by the polynomials in a given polynomial system is one of the main parameters controlling the complexity of solving the system (see, for example, [2] for a complexity analysis of Faugère's $F_5$ algorithm in the context of homogeneous polynomials). Thus, if one can avoid the strategy of variable duplication one can potentially significantly reduce the exponent $5.26n(nk + k + 1)$ in the arithmetic complexity result of proposition 12.

## 4.2 Reducing the system of DDEs to a single functional equation

### 4.2.1 Main strategy

In this subsection we elaborate on a strategy for solving systems of DDEs that avoids duplication of variables: it reduces the initial system to a single functional equation by eliminating the bivariate series $F_2, \dots, F_n$ from the system $(E_1(u) = 0, \dots, E_n(u) = 0)$. In a favourable situation, this reduction outputs a nonzero polynomial $E \in K(t)[x_1, u, z_0, \dots, z_{nk}]$ such that

$$ E(u) \equiv E(F_1(t, u), u, F_1(t, a), \dots, (\partial_u^{k-1} F_1)(t, a), \dots, F_n(t, a), \dots, (\partial_u^{k-1} F_n)(t, a)) = 0. $$
The following lemma ensures that this strategy works under the following assumption: the saturated ideal $\langle E_1, \dots, E_n \rangle : \mathrm{Det}^\infty \cap \mathbb{K}(t)[x_1, u, z_0, \dots, z_{nk-1}]$, denoted by $\mathcal{J}$, is principal<fnref n="9" /> – this natural assumption is observable on all examples of systems of DDEs we encountered so far.

**Lemma 14.** Consider (5) and denote $E_1, \dots, E_n \in \mathbb{K}(t)[x_1, \dots, x_n, u, z_0, \dots, z_{nk-1}]$ the polynomials obtained after taking the numerators of (5). Assume that the ideal $\mathcal{J}$ is principal, generated by some $E \in \mathbb{K}(t)[x_1, u, z_0, \dots, z_{nk-1}]$. Then $E \neq 0$ and $E(u) = 0$.

*Proof.* It results from the product of the diagonal elements in the matrix defining Det that $\mathrm{Det} \neq 0 \mod t$. Also, the evaluation $\mathrm{Det}(u) \in \mathbb{K}[[t]][[u]]$ is not equal to 0 (recall that $\mathrm{Det}(u)$ is the specialization of Det to the series $F_i(t, u)$ and $(\partial_u^\ell F_i)(t, a))$. It follows from the Jacobian criterion that the ideal $\langle E_1, \dots, E_n \rangle : \mathrm{Det}^\infty$ is radical and of dimension at most 0 over $\mathbb{K}(t, z_0, \dots, z_{nk-1})$. Thus $E \neq 0$. Moreover, as $E_1(u) = 0, \dots, E_n(u) = 0$ and $\mathrm{Det}(u) \neq 0$, it follows $E(u) = 0$. □

**Example 2 (cont.).** We compute a Gröbner basis of $\langle E_1, E_2, m \cdot \mathrm{Det} -1 \rangle \cap \mathbb{K}(t)[x_1, u, z_0, z_1]$ and observe that it is not reduced to $\{0\}$. Moreover, it is generated by a unique element

$$E := -(x_1 - 1)(u - 1) + tu(2ux_1^2 - uz_0 - 2x_1^2 + u + x_1 - 1) \in \mathbb{K}(t)[x_1, u, z_0, z_1].$$

For the remaining part of section 4.2 we assume, as in lemma 14, that the ideal $\mathcal{J}$ is principal. Thus the polynomial $E$ obtained after eliminating $x_2, \dots, x_n$ satisfies $E \neq 0$ and $E(u) = 0$.

A natural idea now is to use Bousquet-Mélou and Jehanne's method [10, Section 2]. This reduces our problem either to solving a polynomial system with now $3nk + 1$ unknowns and equations (where again the +1 comes from the saturation polynomial sat), or to applying the more recent algorithm from [9, Section 5]. Note that this method is not guaranteed to work in general because the functional equation given by $E(u) = 0$ may not be of a fixed point type for $F_1$. As we will see, in order to make these approaches work, one can require the following condition:

**Hypothesis 2:** $\bullet$ $\partial_{x_1} E(u) = 0$ admits $nk$ distinct solutions (in $u$) in $\overline{\mathbb{K}}[[t^\frac{1}{\star}]] \setminus \overline{\mathbb{K}}$, (H2)

$\bullet$ The polynomial system obtained after duplicating $nk$ times the variables $x_1, u$ in $(E, \partial_{x_1} E, \partial_u E)$ and adding sat = 0 induces an ideal of dimension 0 over $\mathbb{K}(t)$.

We emphasize that we are not aware of any system of DDEs inducing some $E$ that is not of a fixed-point type in $F_1$: when the system is generic (that is, for generic choices of $f_i$'s and $Q_i$'s) we could observe by generating lots of examples that this fixed-point nature is indeed satisfied.

This resolution strategy being different from the duplication of variables approach from section 4.1, we investigate below how their outputs compare.

### 4.2.2 Theoretical comparison with the strategy from section 4.1

The following proposition ensures that the first part of (H1) implies the first part of (H2). Recall that $\mathcal{J} := \langle E_1, \dots, E_n \rangle : \mathrm{Det}^\infty \cap \mathbb{K}(t)[x_1, u, z_0, \dots, z_{nk-1}]$.

<footnote n="9">⁹Recall that an ideal is called principal if it is generated by only one element.</footnote>
**Proposition 15.** Let $U(t) \in \overline{\mathbb{K}}[[t^\frac{1}{\lambda}]]$ be a solution of $\mathrm{Det}(u) = 0$ such that none of the $(n-1) \times (n-1)$ minors of $(\partial_{x_j} E_i(U(t)))_{1\le i,j\le n}$ is zero. Consider $E \in \mathcal{J}$. Then $U(t)$ is also a solution in $u$ of the equation $\partial_{x_1} E(u) = 0$.

*Proof.* Since $E \in \mathcal{J}$, there exist polynomials $V_1, \dots, V_n \in \overline{\mathbb{K}}(t)[x_1, \dots, x_n, u, z_0, \dots, z_{nk-1}]$ such that $E(U) = \sum_{i=1}^n E_i(U)V_i(U)$. Differentiating with respect to $x_j$ for $j = 1, \dots, n$ and using that $E_i(U) = 0$ and that $E$ does not depend on $x_j$ for $j \ge 2$, we find

$$ \begin{pmatrix} \partial_{x_1} E(U) \\ 0 \\ \vdots \\ 0 \end{pmatrix} = \begin{pmatrix} \partial_{x_1} E_1(U) & \dots & \partial_{x_1} E_n(U) \\ \vdots & \ddots & \vdots \\ \partial_{x_n} E_1(U) & \dots & \partial_{x_n} E_n(U) \end{pmatrix} \begin{pmatrix} V_1(U) \\ \vdots \\ V_n(U) \end{pmatrix}. \quad (28) $$

By definition of $U$, the matrix $(\partial_{x_j} E_i(U))_{i,j}$ is singular and each of its $(n-1) \times (n-1)$ minors is nonzero. It follows that we can express the first row of the matrix as a linear combination of the other rows, then (28) implies that $\partial_{x_1} E(U) = 0$. $\square$

It remains to understand the second part of (H2). We reuse the notations from the proof of proposition 15 by writing $E = \sum_{i=1}^n E_i \cdot V_i$ for $E$ a generator of $\mathcal{J}$ and for $V_1, \dots, V_n$ polynomials in $\overline{\mathbb{K}}(t)[x_1, \dots, x_n, u, z_0, \dots, z_{nk-1}]$. Also, we define the new geometric assumption (P):

**Hypothesis (P):** Every point $\beta \in \overline{\mathbb{K}}(t)^{nk+2}$ vanishing $E$ is the projection onto (P) the $\{x_1, u, z_0, \dots, z_{nk-1}\}$-coordinate space of a point $\alpha \in \overline{\mathbb{K}}(t)^{nk+n+1}$ vanishing simultaneously $E_1, \dots, E_{n-1}$ and $E_n$.

Observe that the zero set of $E$ in $\overline{\mathbb{K}}(t)^{nk+2}$ is, by the closure theorem [18, Theorem 3, Chap 3.2], the Zariski closure in $\overline{\mathbb{K}}(t)^{nk+2}$ of the projection of the solution set in $\overline{\mathbb{K}}(t)^{nk+n+1}$ of $\langle E_1, \dots, E_n \rangle : \mathrm{Det}^\infty$ onto the $\{x_1, u, z_0, \dots, z_{nk-1}\}$-coordinate space. Assumption (P) formulates that the boundary of the projection belongs to the projection itself.

From now on we will denote by $\underline{x}$ (resp. $\underline{z}$) the variables $x_1, \dots, x_n$ (resp. $z_0, \dots, z_{nk-1}$).

**Remark 16.** Once a Gröbner basis $G$ of the ideal $\langle E_1, \dots, E_n \rangle : \mathrm{Det}^\infty$ is computed for the block order $\{x_2, \dots, x_n\} >_{lex} \{x_1, u, z\}$, the set of points $\beta$ in (P) can be characterized by the extension theorem [18, Theorem 3, Chap 3.1] as the solutions of explicit polynomial equations and inequations built from $G$: namely $\beta$ are those solutions of $G \cap \overline{\mathbb{K}}(t)[x_1, u, z]$ that do not vanish simultaneously all the leading terms of $G \setminus (G \cap \overline{\mathbb{K}}(t)[x_1, u, z])$ for the degrevlex monomial order $\{x_2, \dots, x_n\}$ (the variables $t, x_1, u, z$ shall be seen as parameters once $G$ is computed). These leading coefficients are thus polynomials of $\overline{\mathbb{K}}(t)[x_1, u, z]$.

We formulate the below statement, whose purpose is to interpret geometrically the link between the algebraic sets $V(E, \partial_{x_1} E, \partial_u E) \subset \overline{\mathbb{K}}(t)^{nk+2}$ and $V(E_1, \dots, E_n, \mathrm{Det}, P) \subset \overline{\mathbb{K}}(t)^{nk+n+1}$.

**Proposition 17.** Assume that $\mathcal{J}$ is principal and radical. Let $E \in \overline{\mathbb{K}}(t)[x_1, u, z_0, \dots, z_{nk-1}]$ be its generator. Assume, moreover, that (P) holds and denote by $\pi_{x_1, u, z} : (\underline{x}, u, \underline{z}) \in \overline{\mathbb{K}}(t)^{nk+n+1} \mapsto (x_1, u, z) \in \overline{\mathbb{K}}(t)^{nk+2}$ the canonical projection onto the $(x_1, u, z)$-coordinate space. Then we have the inclusion $V(E, \partial_{x_1} E, \partial_u E) \subset \pi_{x_1, u, z}(V(E_1, \dots, E_n, \mathrm{Det}, P))$.
Proof. Writing $E = \sum_{i=1}^n E_i \cdot V_i$ for some $V_1, \dots, V_n \in \mathbb{K}(t)[\underline{x}, u, \underline{z}]$ and considering the derivatives $\partial_{x_i} E$ for $i = 1, \dots, n$, we find by the Leibniz rule

$$ \begin{pmatrix} \partial_{x_1} E \\ 0 \\ \vdots \\ 0 \end{pmatrix} = \begin{pmatrix} \partial_{x_1} E_1 & \dots & \partial_{x_1} E_n \\ \vdots & \ddots & \vdots \\ \partial_{x_n} E_1 & \dots & \partial_{x_n} E_n \end{pmatrix} \begin{pmatrix} V_1 \\ \vdots \\ V_n \end{pmatrix} + \begin{pmatrix} \partial_{x_1} V_1 & \dots & \partial_{x_1} V_n \\ \vdots & \ddots & \vdots \\ \partial_{x_n} V_1 & \dots & \partial_{x_n} V_n \end{pmatrix} \begin{pmatrix} E_1 \\ \vdots \\ E_n \end{pmatrix}. \quad (29) $$

Now for any $\beta \in V(E, \partial_{x_1} E)$ and for all $\alpha \in \pi_{x_1, u, \underline{z}}^{-1}(\beta) \cap V(E_1, \dots, E_n)$, we obtain that

$$ \begin{pmatrix} 0 \\ \vdots \\ 0 \end{pmatrix} = \begin{pmatrix} \partial_{x_1} E_1(\alpha) & \dots & \partial_{x_1} E_n(\alpha) \\ \vdots & \ddots & \vdots \\ \partial_{x_n} E_1(\alpha) & \dots & \partial_{x_n} E_n(\alpha) \end{pmatrix} \begin{pmatrix} V_1(\alpha) \\ \vdots \\ V_n(\alpha) \end{pmatrix}. $$

Thus it suffices to prove that $V_i(\alpha) \neq 0$ (for some $i \in \{1, \dots, n\}$) in order to see that $\mathrm{Det}(\alpha) = 0$. As $E$ generates the radical ideal $\mathcal{J}$, we know that $E$ is squarefree. Hence the algebraic set $V(E) \in \overline{\mathbb{K}(t)}^{nk+2}$ is smooth, hence the vector of partial derivatives of $E$ w.r.t the variables $x_1, u, z_0, \dots, z_{nk-1}$ does not vanish identically at $\beta$. As $E = \sum_{i=1}^n E_i \cdot V_i$, this enforces one of the $V_i$ to be nonzero at $\alpha$: else by writing all the partial derivatives of $E$ with the Leibniz rule applied to the sum expression $E = \sum_{i=1}^n E_i \cdot V_i$ and using $\alpha \in V(E_1, \dots, E_n, V_1, \dots, V_n)$ would imply that the vector of partial derivatives of $E$ with respect to $x_1, u, z_0, \dots, z_{nk-1}$ vanishes identically at $\beta$.

Now consider the vector of derivatives $(\partial_u E, \partial_{x_2} E, \dots, \partial_{x_n} E)$ and write it in the form (29). This allows to reuse the above argument to show that $P(\alpha) = 0$, for any $\alpha \in \overline{\mathbb{K}(t)}^{nk+n+1}$ vanishing simultaneously $\partial_u E, E_1, \dots, E_n$. It follows that under assumption (P), we have the claimed inclusion $V(E, \partial_{x_1} E, \partial_u E) \subset \pi_{x_1, u, \underline{z}}(V(E_1, \dots, E_n, \mathrm{Det}, P))$. $\square$

Let us mention that under the technical assumptions introduced in proposition 15 and proposition 17: if $\mathrm{Det}(u) = 0$ admits $\ell$ distinct solutions (in $u$) in $\overline{\mathbb{K}}[[t^\frac{1}{\ell}]]$ then so does $\partial_{x_1} E(u) = 0$.

Moreover still under the assumptions of proposition 15 and proposition 17, if the ideal $\mathcal{I}_{\mathrm{dup}}^\infty$ of section 4.1.1 has dimension 0 over $\mathbb{K}(t)$, then the application of [9, Section 5] to the functional equation $E(u) = 0$ outputs a nonzero polynomial $R \in \mathbb{K}[t, z_0]$ such that $R(t, F_1(t, a)) = 0$.

We gather below assumptions that we will use after. Also for any multivariate polynomial $p$, we denote by $\mathrm{SqFreePart}(p)$ the squarefree part of $p$.

**Hypothesis 3:** • (H1) and (P) hold, (H3)

• The ideal $\mathcal{J}$ is radical and principal,

• For all $U(t) \in \overline{\mathbb{K}}[[t^\frac{1}{\ell}]]$ solution of $\mathrm{Det}(u) = 0$, the $(n-1) \times (n-1)$ minors of $(\partial_{x_j} E_i(U(t)))_{1 \le i, j \le n}$ are nonzero.

**Proposition 18.** Assume that (H3) holds. Denote $R_1 \in \mathbb{K}[t, z_0]$ a generator of $\mathcal{I}_{\mathrm{dup}}^\infty \cap \mathbb{K}[t, z_0]$ and $R_2 \in \mathbb{K}[t, z_0]$ the polynomial computed by the duplication approach$^{10}$ applied to the functional equation $E(u) = 0$. Then $\mathrm{SqFreePart}(R_2)$ divides $\mathrm{SqFreePart}(R_1)$.

**Proof.** The statement is a direct consequence of proposition 15 and proposition 17. $\square$

<footnote>$^{10}$This duplication approach should in practice be replaced by the more efficient algorithm from [9, Section 5].</footnote>
In this section, so far, we compared two different strategies based for the first one on a duplication of variables argument, and for the second one on the reduction to a single equation. We showed that up to considering squarefree parts, the output polynomial of the second strategy divides, under technical assumptions, the output polynomial of the first strategy.

Even if we strongly believe that the introduced assumptions are reasonable for practical considerations, we shall however underline that the literature contains some degenerate examples that go against this believe. The following example illustrates this degeneracy and thus strengthen the relevancy of the strategy from section 4.1, even if the arithmetic complexity in proposition 12 can seem despairing.

**Example 2 (cont.).** In this case proposition 15 cannot be applied, since

$$ \det = \det \begin{pmatrix} 4tu^2x_1 - 4tux_1 + tu - u + 1 & 0 \\ * & 2tu^2x_1 - 2tux_1 + tu - u + 1 \end{pmatrix}, $$

so any of the two distinct solutions $U_1, U_2 \in \overline{\mathbb{K}}[[t^\frac{1}{\star}]]$ to the equation $\det(u) = 0$ annihilates at least two coefficients of the matrix $(\partial_{x_j} E_i)_{1 \le i, j \le 2}$. Moreover, it is straightforward to see that for the polynomial $E = -(x_1 - 1)(u - 1) + tu(2ux_1^2 - uz_0 - 2x_1^2 + u + x_1 - 1) \in \overline{\mathbb{K}}(t)[x_1, u, z_0, z_1]$ the equation $\partial_{x_1} E(u) = 0$ has only one solution (in $u$) in the ring $\overline{\mathbb{K}}[[t^\frac{1}{\star}]]$ while $E(u) = 0$ involves two univariate series $F_1(t, 1)$ and $F_2(t, 1)$. Thus, on this example, the strategy described in the current section fails.

# 5 Conclusion and future works

We proved in theorem 3 that solutions of systems of discrete differential equations are algebraic functions. The proof uses the observation that, up to a symbolic deformation pushing all degeneracy of a given system of the form (5) to some higher powers of $t$, all the computations become explicit. In addition, we obtained in theorem 4 quantitative estimates for the size of an annihilating polynomial for any solution $F_i$, together with arithmetic complexity estimates.

We also refined the complexity estimate in proposition 12 and also obtained an algebraicity bound for all the $F_i(t, a)$. Moreover, we started an algorithmic analysis of the problem of solving systems of DDEs. In this first work we compared two natural strategies for solving systems of DDEs. We identified rigorous assumptions ensuring their applicability illustrating them on the Example 2, where the reduction to a single equation fails, while the duplication approach works. Under these assumptions, the outputs of these two strategies are linked via proposition 18.

Our paper contains several research directions which we plan to investigate in the future:

**Algorithmic part** This article set the foundations of an algorithmic work that must be conducted in order to understand better the many subtle and degenerate situations that one can encounter with systems of DDEs from the literature. For instance, the following points should be addressed:

* **Validity of the assumptions:** It seems clear that up to deforming symbolically as in (7), assumption (H1) can be assumed to hold. Still, we do not know at this stage if (P) holds once (5) is deformed as in (7). Also, we would like to have results saying that the all introduced assumptions hold for generic choices of system of DDEs.

* **New algorithms for solving systems of DDEs:** In [9, Section 5] an algorithm was designed that computes in case of scalar DDEs the same output as the algorithm described in section 4.1
and that avoids the duplication of variables. This strategy can be naturally extended to the
context of systems of DDEs. Once this is done, a first task would be to compare this new
strategy with the ones described in section 4.2. Also, the case of systems of linear DDEs
contains many applications in which it is possible to use the linearity in order to speed up
the computations. All this should be clarified and studied in depth.

**Theoretical part:** In the case of more than one catalytic variables, the algebraicity of the solutions remains guaranteed by Popescu's theorem [23] under the additional assumption that the variables are “nested”. It would be interesting to reprove this result with elementary tools in the case of fixed-point equations and to introduce new algorithms for solving these equations.

**Acknowledgments** We would like to thank Mohab Safey El Din for many interesting discussions and Alin Bostan for his great support and valuable comments on an early version of this manuscript. We are also grateful to Mireille Bousquet-Mélou for exciting and enlightening discussions, and to Eva-Maria Hainzl and Wenjie Fang for pointing out to us the existence of the periodic Tamari intervals and their associated systems of DDEs. Also, we are grateful to Mathilde Bouvel and Nicholas Beaton for pointing us applications of our paper to [4].

Both authors were supported by the ANR-19-CE40-0018 De Rerum Natura project and the French-Austrian ANR-22-CE91-0007 & FWF I6130-N project EAGLES; the second author was funded by the DOC fellowship P-26101 of the ÖAW, the WTZ collaboration project FR-09/2021 and the project P34765 of the Austrian science fund FWF.
References

[1] Andrei Asinowski, Axel Bacher, Cyril Banderier, and Bernhard Gittenberger. "Analytic combinatorics of lattice paths with forbidden patterns, the vectorial kernel method, and generating functions for pushdown automata". In: Algorithmica 82.3 (2020), pp. 386-428. ISSN: 0178-4617. DOI: 10.1007/s00453-019-00623-3. URL: https://doi.org/10.1007/s00453-019-00623-3.

[2] Magali Bardet, Jean-Charles Faugère, and Bruno Salvy. "On the complexity of the $F_5$ Gröbner basis algorithm". In: J. Symbolic Comput. 70 (2015), pp. 49-70. ISSN: 0747-7171,1095-855X. DOI: 10.1016/j.jsc.2014.09.025. URL: https://doi.org/10.1016/j.jsc.2014.09.025.

[3] Walter Baur and Volker Strassen. "The complexity of partial derivatives". In: Theoret. Comput. Sci. 22.3 (1983), pp. 317-330. ISSN: 0304-3975,1879-2294. DOI: 10.1016/0304-3975(83)90110-X. URL: https://doi.org/10.1016/0304-3975(83)90110-X.

[4] Nicholas R. Beaton, Mathilde Bouvel, Veronica Guerrini, and Simone Rinaldi. "Slicing of parallelogram polyominoes: Catalan, Schröder, Baxter, and other sequences". In: Electron. J. Combin. 26.3 (2019), Paper No. 3.13, 36. DOI: 10.37236/7375. URL: https://doi.org/10.37236/7375.

[5] Olivier Bernardi and Mireille Bousquet-Mélou. "Counting colored planar maps: algebraicity results". In: J. Combin. Theory Ser. B 101.5 (2011), pp. 315-377. ISSN: 0095-8956. DOI: 10.1016/j.jctb.2011.02.003. URL: https://doi.org/10.1016/j.jctb.2011.02.003.

[6] Jérémy Berthomieu, Christian Eder, and Mohab Safey El Din. "msolve: A Library for Solving Polynomial Systems". In: 2021 International Symposium on Symbolic and Algebraic Computation. ISSAC '21. Saint Petersburg, Russia: ACM, July 2021, pp. 51-58. DOI: 10.1145/3452143.3465545.

[7] Nicolas Bonichon, Mireille Bousquet-Mélou, Paul Dorbec, and Claire Pennarun. "On the number of planar Eulerian orientations". In: European J. Combin. 65 (2017), pp. 59-91. ISSN: 0195-6698. DOI: 10.1016/j.ejc.2017.04.009.

[8] Alin Bostan, Frédéric Chyzak, Hadrien Notarantonio, and Mohab Safey El Din. "Algorithms for Discrete Differential Equations of Order 1". In: Proceedings of the 2022 International Symposium on Symbolic and Algebraic Computation. ISSAC '22. Villeneuve-d'Ascq, France: ACM, 2022, pp. 101-110. ISBN: 9781450386883.

[9] Alin Bostan, Hadrien Notarantonio, and Mohab Safey El Din. "Fast algorithms for discrete differential equations". In: Proceedings of the International Symposium on Symbolic & Algebraic Computation (ISSAC 2023). ACM, New York, [2023] ©2023, pp. 80-89. ISBN: 979-8-4007-0039-2. DOI: 10.1145/3597066.3597103. URL: https://doi.org/10.1145/3597066.3597103.

[10] Mireille Bousquet-Mélou and Arnaud Jehanne. "Polynomial equations with one catalytic variable, algebraic series and map enumeration". In: J. Combin. Theory Ser. B 96.5 (2006), pp. 623-672. ISSN: 0095-8956.

[11] E. Brézin, C. Itzykson, G. Parisi, and J. B. Zuber. "Planar diagrams". In: Comm. Math. Phys. 59.1 (1978), pp. 35-51. URL: http://projecteuclid.org/euclid.cmp/1103901558.

[12] W. G. Brown and W. T. Tutte. "On the Enumeration of Rooted Non-Separable Planar Maps". In: Canadian J. of Mathematics 16 (1964), pp. 572-577. DOI: 10.4153/CJM-1964-058-7.

[13] William G. Brown. "Enumeration of quadrangular dissections of the disk". In: Canadian J. Math. 17 (1965), pp. 302-317. ISSN: 0008-414X. DOI: 10.4153/CJM-1965-030-1.

[14] William G. Brown. "On the existence of square roots in certain rings of power series". In: Math. Ann. 158 (1965), pp. 82-89. ISSN: 0025-5831. DOI: 10.1007/BF01370732.

[15] Manfred Buchacher and Manuel Kauers. "Inhomogeneous restricted lattice walks". In: Sém. Lothar. Combin. 82B (2020), Art. 75, 12.
[16] Luca Castelli Aleardi, Olivier Devillers, and Gilles Schaeffer. "Optimal succinct representations of planar maps". In: *Computational geometry (SCG'06). ACM, New York, 2006, pp. 309–318. DOI: 10.1145/1137856.1137902*. URL: https://doi.org/10.1145/1137856.1137902.

[17] David A. Cox. "Stickelberger and the eigenvalue theorem". In: *Commutative algebra*. Springer, Cham, [2021] ©2021, pp. 283–298. ISBN: 978-3-030-89693-5; 978-3-030-89694-2. DOI: 10.1007/978-3-030-89694-2_8. URL: https://doi.org/10.1007/978-3-030-89694-2_8.

[18] David A. Cox, John Little, and Donal O'Shea. *Ideals, varieties, and algorithms*. Fourth. Undergraduate Texts in Mathematics. An introduction to computational algebraic geometry and commutative algebra. Springer, Cham, 2015, pp. xvi+646. DOI: 10.1007/978-3-319-16721-3. URL: https://doi.org/10.1007/978-3-319-16721-3.

[19] David Eisenbud. *Commutative algebra*. Vol. 150. Graduate Texts in Mathematics. With a view toward algebraic geometry. Springer-Verlag, New York, 1995, pp. xvi+785. ISBN: 0-387-94268-8; 0-387-94269-6.

[20] Joos Heintz. "Definability and fast quantifier elimination in algebraically closed fields". In: *Theoret. Comput. Sci.* 24.3 (1983), pp. 239–277. DOI: 10.1016/0304-3975(83)90002-6. URL: https://doi.org/10.1016/0304-3975(83)90002-6.

[21] Hadrien Notarantonio. *DDE-Solver: A Maple package for Discrete Differential Equations*. Preprint.

[22] Hadrien Notarantonio and Sergey Yurkevich. *Effective algebraicity for solutions of systems of functional equations with one catalytic variable*. In *FPSAC'23*, to appear (2023).

[23] Dorin Popescu. "General Néron desingularization and approximation". In: *Nagoya Math. J.* 104 (1986), pp. 85–115. ISSN: 0027-7630. DOI: 10.1017/S0027763000022698.

[24] Éric Schost. "Computing parametric geometric resolutions". In: *Appl. Algebra Engrg. Comm. Comput.* 13.5 (2003), pp. 349–393. ISSN: 0938-1279. DOI: 10.1007/s00200-002-0109-x.

[25] W. T. Tutte. "A Census of Planar Triangulations". In: *Canadian J. Math.* 14 (1962), pp. 21–38. DOI: 10.4153/CJM-1962-002-9.

[26] Gilles Villard. "On computing the resultant of generic bivariate polynomials". In: *ISSAC'18—Proceedings of the 2018 ACM International Symposium on Symbolic and Algebraic Computation*. ACM, New York, 2018, pp. 391–398. DOI: 10.1145/3208976.3209020. URL: https://doi.org/10.1145/3208976.3209020.

[27] Doron Zeilberger. "The umbral transfer-matrix method. I. Foundations". In: vol. 91. 1-2. In memory of Gian-Carlo Rota. 2000, pp. 451–463. DOI: 10.1006/jcta.2000.3110.