1. Introduction
The propagation of high-intensity light pulses as they propagate through an optical fibre or waveguide is principally affected by dispersion, the frequency-dependence of the refractive index, and the Kerr nonlinearity, the intensity dependence of the refractive index.Footnote 1 While dispersion and nonlinearity can be modelled to almost arbitrary degrees of accuracy [Reference Agrawal3], considerable physical insight can be had by taking the simplest models. These usually have two elements: (i) the dispersion is quadratic, so the inverse group velocity depends linearly on frequency with a slope of β 2; and (ii) the propagation constant of the light has a component that is linearly proportional to the light intensity, with slope γ. Under these conditions, light propagation is well described by the nonlinear Schrödinger equation
\begin{equation}
i \frac{\partial \psi}{\partial z}-\frac{\beta_{2}}{2!} \frac{\partial^{2} \psi}{\partial T^{2}}+\gamma |\psi|^2\psi=0,
\end{equation}which is well-understood and is known to be integrable [Reference Zakharov and Shabat28]. Here z is the propagation coordinate, ψ is the envelope of the electric field and T is time in a frame moving with ψ. We denote the evolution parameter as z to be consistent with the optics literature.
Many generalizations of the nonlinear Schrödinger equation have been proposed and investigated, including more general nonlinear terms [Reference Biswas and Konar7], coupled equations to account for the polarization of light and the inclusion of higher derivatives, representing a small correction to the dominant quadratic dispersion [Reference Kivshar and Agrawal17, Reference Piché, Cormier and Zhu20]. Here we are interested in generalizations in which the dispersion is not dominated by the quadratic term, so that the evolution equations take the form
\begin{equation}
i \frac{\partial \psi}{\partial z}+\sum_{k=1}^{M}(-1)^k \frac{\beta_{2k}}{(2k)!} \frac{\partial^{2k} \psi}{\partial T^{2k}}+\gamma |\psi|^2\psi=0,
\end{equation} where γ > 0, 
$\beta_{2M}\neq 0$ and which reverts to Eq. (1.1) when M = 1. Here 2M is both the highest dispersion order and the highest order derivative.
Interest in Eq. (1.2) is driven by recent experiments where M takes values up to M = 5 and all lower derivatives vanish [Reference Blanco-Redondo, de Sterke, Sipe, Krauss, Eggleton and Husko8, Reference Runge, Hudson, Tam, de Sterke and Blanco-Redondo24, Reference Runge, Qiang, Alexander, Rafat, Hudson, Blanco-Redondo and de Sterke25]. In other experiments, M takes values up to M = 8 and has lower derivative terms [Reference Lourdesamy, Runge, Alexander, Hudson, Blanco-Redondo and de Sterke19]. These experiments highlight that in a waveguide, the dispersion is a controllable parameter and can be engineered to suit a particular experimental need. Given this interest, it is pertinent to investigate the formal mathematical properties of Eq. (1.2) and of the ordinary differential equation (ODE) that derive from it (see below). In this paper, we consider the possible integrability of these equations.
We focus on solitary wave solutions of the form 
$\psi(z,T)=e^{i\mu z}v(T)$, which are governed by the ODE
\begin{equation}
-\mu\, v(T)+\sum_{k=1}^{M}(-1)^k \frac{\beta_{2k}}{(2k)!} v^{(2k)}(T)+\gamma\, v(T)^3=0,
\end{equation}where v can be taken to be real and where its superscript indicates a time derivative of the associated order. Such solitary wave solutions play an important role in physical settings [Reference Agrawal3, Reference de Sterke, Runge, Hudson and Blanco-Redondo11, Reference Kivshar and Agrawal17, Reference Lourdesamy, Runge, Alexander, Hudson, Blanco-Redondo and de Sterke19, Reference Runge, Hudson, Tam, de Sterke and Blanco-Redondo24, Reference Runge, Qiang, Alexander, Rafat, Hudson, Blanco-Redondo and de Sterke25] and the corresponding exact reduction of the partial differential equation (PDE) to the above ODE is suitable to study integrability via the Ablowitz–Ramani–Segur conjecture [Reference Ablowitz, Ramani and Segur1, Reference Ablowitz, Ramani and Segur2, Reference Hone15, Reference Weiss, Tabor and Carnevale27].
Applying a scaling
\begin{equation*}
v=\left(\frac{\mu}{\gamma}\right)^{\frac{1}{2}} u,\qquad T=d\cdot t,
\end{equation*} the ODE for 
$u=u(t)$ can be written as
\begin{equation}
u^3=u+\sum_{k=1}^M \alpha_{2k} \,u^{(2k)},
\end{equation}with coefficients
\begin{equation*}
\alpha_{2k}=(-1)^{k+1}\frac{\beta_{2k}}{\mu\, (2k)!\,d^{2k}}\qquad (1\leq k\leq M).
\end{equation*}For the purpose of studying singularities in later sections (Section 2 onwards), it is convenient to choose d such that the highest order coefficient is normalized to equal
\begin{equation}
\alpha_{2M}=\frac{(M-1)!}{(3M-1)!},
\end{equation} while for the simpler examples, d is chosen such that 
$\alpha_{2M}=1$.
The requirement of the absence of odd derivatives ensures that the solution can be made real.
Some analytic solutions to Eq. (1.4) are known for specific ratios of the coefficients. Karlsson and Höök [Reference Karlsson and Höök16], for example, showed that for M = 2 and for particular values of the 
$\alpha_{2k}$ the equation has a solution in the form of the square of a hyperbolic secant. This work was subsequently generalized to higher M and to hyperbolic secants of higher powers and linear combinations thereof [Reference Kudryashov18, Reference Qiang, Alexander and de Sterke21]. More recently, analytic solutions in the form of a rapidly converging series of functions were reported by Qiang et al. [Reference Qiang, Broderick and de Sterke22]. Nonetheless, the general mathematical properties of Eq. (1.4) and the PDE Eq. (1.2) from which they are derived remain largely unexplored.
As mentioned, the nonlinear Schrödinger equation (1.1) is integrable [Reference Zakharov and Shabat28]. However, the integrability of PDE Eq. (1.2) for M > 1, or that of the ODE Eq. (1.4) which derive from it, has not been investigated systematically. Integrability would mean that the properties and tools of integrable systems, such as the existence of sufficiently many conserved quantities and that of Lax pairs, can be used to investigate these equations.
The aim of this paper is to investigate the integrability and, closely related to that, singularities of these differential equations. We first investigate the integrability of ODE (Eq. 1.4) by applying the Painlevé test [Reference Conte9, Reference Goriely12, Reference Hone15]. From this, we draw conclusions regarding the integrability of PDEs of the form of Eq. (1.2). We then analyse singularities of solutions and find that they are poles at leading order with, generally, multi-valued subdominant contributions that obstruct the passing of the Painlevé test. Through a numerical analysis, we further show that such branch-point singularities are present even in the pure quartic solitary wave when analytically continued in the complex plane.
The paper is organized as follows. In Section 2, we apply the Painlevé test to ODE (1.4). In Section 3, we investigate singularities of solutions of ODE (1.4). This is followed by a numerical analysis of a fundamental solution and its singularities in Section 4. The paper ends with a conclusion and discussion in Sections 6 and 5. Appendix A contains a technical estimate.
2. Painlevé analysis of the ODE
In this section, we investigate for what values of the parameters, if any, ODE (1.4) has the Painlevé property. As a reminder, an ODE has the Painlevé property when it does not have any movable branch points. In that case, we also call the ODE Painlevé integrable. Furthermore, an ODE is said to have the weak Painlevé property if all movable branch points are algebraic [Reference Ramani, Dorizzi and Grammaticos23].
2.1. Derivation of the resonance polynomial
Various tests have been developed that check the necessary conditions for Painlevé integrability [Reference Conte9, Reference Hone15]. For our purposes, the basic Painlevé test will suffice. The first step of this test is to find all possible exponents 
$p\in\mathbb{C}$, different from 
$0,1,\ldots,M-1$, such that substitution of 
$u(t)\sim u_0(t-t_0)^p$ into ODE (1.4) gives a consistent dominant balance. All such exponents are required to be integers for the ODE to pass the first step of the test.
Substitution gives the dominant balance
\begin{equation}
u_0^3(t-t_0)^{3p}\sim\alpha_{2M}\,u_0\, (t-t_0)^{p-2M}\,p\cdot (p-1)\cdot\ldots\cdot(p-2M+1).
\end{equation} This balance is consistent only when 
$p=-M$, in which case, necessarily 
$u_0^2=\pm 1$, due to the normalization of 
$\alpha_{2M}$ in (1.5). Thus, the only non-trivial value that the exponent can take is 
$p=-M$, which is integer and hence ODE (1.4) passes the first step in the test.
For the second part of the test, we develop u(t) in a Laurent series around 
$t=t_0$ with a pole of order M at 
$t=t_0$. Since ODE (1.4) only involves even derivatives and odd powers of u, 
$(t-t_0)^M u(t)$ must be even around 
$t=t_0$, and thus the Laurent series of u(t) takes the form
\begin{equation}
u(t)=\sum_{k=0}^\infty u_k(t-t_0)^{-M+2k}.
\end{equation} Recall further that 
$u_0=\pm 1$, and since ODE (1.4) is invariant under 
$u\mapsto -u$, we may assume 
$u_0=1$ without loss of generality.
For the ODE to pass the second part of the test, it is required that the substitution of the Laurent series into the ODE leaves 
$2M-1$ of the coefficients undetermined, so that, taking into account the free parameter t 0, the Laurent series contains a full 2M degrees of freedom.
Substituting the Laurent series into ODE (1.4) and comparing coefficients of order 
$-3M+2n$, with 
$n\geq 1$, leads to the following relation among the coefficients,
\begin{equation*}
\sum_{i,j\geq 0, i+j\leq n}u_iu_ju_{n-i-j}=\sum_{k=0}^M (3M-2n-1)_{2k}\,
\alpha_{2k}\, u_{k+n-M},
\end{equation*} where we used the notation 
$u_k=0$ for k < 0, 
$\alpha_0:=1$ and 
$(x)_n$ is the falling factorial
\begin{equation*}
(x)_n=\begin{cases}
x\cdot(x-1)\cdot\ldots\cdot(x-n+1) & \text{if }n\geq 1,\\
1 & \text{if }n=0.
\end{cases}
\end{equation*}Solving for un gives the following recurrence relation for the coefficients,
\begin{equation}
\alpha_{2M} P_M(2n) u_n=\sum_{(i,j)\in A_n}u_iu_ju_{n-i-j}-\sum_{k=0}^{M-1}
(3M-2n-1)_{2k}\,\alpha_{2k}\, u_{k+n-M},
\end{equation}where An denotes the set of indices
\begin{equation}
A_n=\{(i,j)\in\mathbb{Z}_{\geq 0}^2:i+j\leq n\}\setminus\{(n,0),(0,n),(0,0)\},
\end{equation} and the resonance polynomial 
$P_M(r)$ is defined by
\begin{align}
P_M(r):&=(r-M)\cdot(r-M-1)\cdot\ldots\cdot (r-3M+1)-\frac{3}{\alpha_{2M}}\\
&=(r-M)_{2M}-(3M)_{2M},\nonumber
\end{align}where the second equality follows from (1.5). We note that we could have bypassed the Laurent series expansion and obtained the resonance polynomial directly [Reference Hone15], but regardless, we need it for the singularity analysis in Section 3.
2.2. Roots of the resonance polynomials
The resonance polynomial 
$P_M(r)$ is a polynomial of degree 2M with a root at 
$r=-1$. For the ODE to pass the second stage of the Painlevé test, it is required that all 
$2M-1$ remaining roots of the resonance polynomial are positive integers. Remarkably, the resonance polynomial is entirely independent of any lower-order dispersion coefficients 
$\alpha_{2k}$, 
$1\leq k\leq M-1$.
For M = 1, we have
\begin{equation*}
P_1(r)=(r-1)(r-2)-6=(r+1)(r-4),
\end{equation*}so that (3.2) passes the second test, as expected. However, setting M = 2, we get
\begin{align*}
P_2(r)&=(r-2)(r-3)(r-4)(r-5)-6\cdot 5\cdot 4\cdot 3\\
&=(r+1)(r-8)(r^2-7r+30).
\end{align*} So, apart from the positive integer root r = 8, we have two non-real roots 
$r=\tfrac{1}{2}(7\pm \sqrt{71}\,i)$. This means that ODE (1.4) with M = 2 is not Painlevé integrable for any value of the lower-order dispersion coefficient α 2.
The roots 
$r=\tfrac{1}{2}(7\pm \sqrt{71}i)$ are not real and, in particular, not rational and thus indicate the existence of movable branch points which are non-algebraic. Therefore, ODE (1.4) with M = 2 does not even have the weak Painlevé property for any value of the lower-order dispersion coefficient α 2.
For general M, we have the following result.
Lemma 1. For any 
$M \geq 1$, the resonance polynomial 
$P_M(r)$ has only two real roots, 
$r=-1$ and 
$r=4M$.
Remark 1.
We remark that the general theory guarantees that 
$r=-1$ is a root of the resonance polynomial [Reference Conte9]. The second real root, 
$r=4M$, is related to the existence of an integral of motion, as explained in Remark 2.
Proof. From the formula 
$P_M(r)=(r-M)_{2M}-(3M)_{2M}$, it is clear that 
$r=-1$ and 
$r=4M$ are roots 
$P_M(r)$. To prove that those are the only real roots, we first show that 
$P_M(r)$ is negative on the interval 
$[M,3M-1]$. To this end, take any 
$r\in [M,3M-1]$ and choose an 
$n\in\mathbb{Z}$ such that 
$n\leq r \lt n+1$. Then 
$M\leq n\leq 3M-1$ and we can estimate
\begin{align*}
|(r-M)_{2M}|=\,&|r-M|\cdot |r-M-1|\cdot\ldots \cdot |r-n+1|\cdot|r-n|\\
&\cdot |r-n-1|\cdot |r-n-2|\cdot\ldots \cdot |r-3M+2|\cdot |r-3M+1|\\
\leq\, &(n+1-M)\cdot (n-M)\cdot \ldots \cdot 2\cdot 1\\
&\cdot 1 \cdot 2\cdot\ldots\cdot (3M-2-n)\cdot (3M-1-n)\\
=\,&(n+1-M)!\cdot (3M-1-n)!
\leq\, (2M)! \lt \, (3M)_{2M},
\end{align*} and thus 
$P_M(r)$ is negative for all 
$r\in [M,3M-1]$.
On the other hand, 
$P_M(r)$ is strictly increasing for 
$r \gt 3M-1$ and strictly decreasing for r < M. Since further 
$P_M(r)\sim r^{2M}$ as 
$r\rightarrow \pm \infty$, this means that 
$P_M(r)$ has exactly two real roots, one in the interval 
$(-\infty,M)$ and another in the interval 
$(3M-1,\infty)$. As we already know that 
$r=-1$ and 
$r=4M$ are roots, these are the only two real roots and the lemma follows.
According to the above lemma, ODE (1.4) fails the second step in the Painlevé test for 
$M\geq 2$, for any values of the lower-order dispersion coefficients 
$\alpha_{2k}$, 
$1\leq k\leq M-1$. Since the resonance polynomial 
$P_M(r)$ has 
$2M-2$ non-real roots, the ODE also does not have the weak Painlevé property for 
$M\geq 2$ and any values of the dispersion coefficients.
As ODE (1.4) is an exact reduction of the PDE (1.2), this means that the original PDE does not satisfy the Painlevé property. In light of the Ablowitz–Ramani–Segur conjecture [Reference Ablowitz, Ramani and Segur1, Reference Ablowitz, Ramani and Segur2, Reference Weiss, Tabor and Carnevale27], it further suggests that the original PDE (1.2) is not solvable by inverse scattering for 
$M\geq 2$ and any values of the dispersion coefficients.
3. Analysis of singularities
In this section, we study the movable singularities of ODE (1.4). We start by introducing a polynomial integral of motion that plays an important role in the analysis in this section.
3.1. An integral of motion
ODE (1.4) has at least one integral of motion, which can be obtained by subtracting the left-hand side from the right-hand side of Eq. (1.4), multiplying the result by 
$u'(t)$ and integrating, yielding
\begin{equation}
\mathcal{I}(u):=\tfrac{1}{2}u^2-\tfrac{1}{4}u^4+\sum_{k=1}^M \alpha_{2k} \mathcal{D}_{k}(u),
\end{equation} where, for 
$k\geq 1$,
\begin{equation*}
\mathcal{D}_{k}(u):=\int u^{(2k)}(t)u'(t)\,dt=(-1)^{k+1}\tfrac{1}{2}(u^{(k)})^2+\sum_{m=1}^{k-1}(-1)^{m+1} u^{(2k-m)}u^{(m)}.
\end{equation*}For example, considering the case M = 1, as normalized in Eq. (1.5),
\begin{equation}
u^3=u+\tfrac{1}{2}u'',
\end{equation}the integral of motion is given by
\begin{equation*}
\mathcal{I}=\tfrac{1}{2}u^2-\tfrac{1}{4}u^4+\tfrac{1}{4}(u')^2.
\end{equation*}3.2. Pole singularities
Pole singularities of solutions of ODE (1.4) are characterized by the following lemma.
Lemma 2. Let 
$M\geq 1$ and u(t) be a solution of ODE (1.4). If u(t) has a pole at 
$t=t_0\in\mathbb{C}$, then this pole has order M and u(t) admits a Laurent expansion around 
$t=t_0$,
\begin{equation}
u(t)=\pm(t-t_0)^{-M}\left(1+\sum_{k=1}^\infty u_k(t-t_0)^{2k}\right),
\end{equation} for a unique sign ± and value of 
$u_{2M}$, where the remaining coefficients are determined by recursion (2.3). Conversely, for any choice of 
$t_0,u_{2M}^*\in\mathbb{C}$ and sign ±, there exist a locally unique solution u(t) of ODE (1.4), meromorphic on an open neighbourhood of t 0, admitting Laurent expansion (3.3) around 
$t=t_0$ with 
$u_{2M}=u_{2M}^*$.
Remark 2.
The proof of Lemma 2 shows that the appearance of one free coefficient, 
$u_{2M}$, in expansion (3.3), is directly tied to the polynomial integral of motion (3.1), and that their values are related through a linear equation. When M = 2, this linear equation reads
\begin{equation*}
\mathcal{I}=\tfrac{89}{540}+\tfrac{212}{45}\alpha_2^2-\tfrac{2104}{15}\alpha_2^4-\tfrac{57}{10}u_4.
\end{equation*}Proof. The forward implication follows from Eqs (2.1), (2.2) and (2.3).
To prove the backward implication, take 
$t_0,u_{2M}^*\in\mathbb{C}$ and fix the sign 
$\pm=+$ without loss of generality. By Lemma 1, 
$P_M(2n)\neq 0$ for 
$1\leq n\leq 2M-1$ and thus the recursive formulas (2.3) determine the values of 
$u_1,u_2,\ldots, u_{2M-1}$ uniquely.
When 
$n=2M$, the left-hand side of Eq. (2.3) vanishes. We will show that the right-hand side vanishes as well, so that 
$u_{2M}$ is free and we may set 
$u_{2M}=u_{2M}^*$. To show this, we use the integral of motion 
$\mathcal{I}$ defined in Eq. (3.1). The equation
\begin{equation}
\mathcal{I}(u)=\mathcal{I}_0,
\end{equation} for any fixed value of 
$\mathcal{I}_0$, is an ODE for u that is essentially equivalent to the original ODE (1.4). Namely, any solution of (1.4) satisfies (3.4), for some constant value of 
$\mathcal{I}_0$; conversely, any solution of (3.4) is either a solution of (1.4) or a constant function.
When we substitute the Laurent series into (3.4), comparing terms of order 
$-4M+2n$ leads to a polynomial relation among 
$u_1,\ldots, u_{n}$. It reads
\begin{equation}
\hat{P}_M(2n) u_n=F_{M,n}(u_1,\ldots,u_{n-1})-\mathbf{1}_{n=2M}\mathcal{I}_0,
\end{equation} where 
$\mathbf{1}_{n=2M}$ is the indicator function of the condition 
$n=2M$, 
$F_{M,n}$ is a polynomial in 
$(u_1,\ldots,u_{n-1})$ and 
$\hat{P}_M(2n)$ is a polynomial in 2n, given by
\begin{align*}
\hat{P}_M(2n)&=1+\alpha_{2M} \sum_{s=1}^{M-1}(-1)^s\left[(-M)_{2M-s}(-M+2n)_s+(-M+2n)_{2M-s}(-M)_s\right] \\
&\quad +\alpha_{2M}(2M-1)_M(-M+2n)_M.
\end{align*} The value of 
$\hat{P}_M(2n)$ at 
$n=2M$ is
\begin{equation*}
\hat{P}_M(4M)=1+3M\left(\frac{1}{M+1}+\frac{1}{M+2}+\ldots+\frac{1}{3M-1}\right),
\end{equation*} which is clearly nonzero. Recurrence (3.5) is hence well-defined at 
$n=2M$ and there is a unique value of 
$\mathcal{I}_0$ for which 
$u_{2M}=u_{2M}^*$. Thus the right-hand side in the original recurrence (2.3) must vanish at 
$n=2M$. Upon having fixed the value 
$u_{2M}=u_{2M}^*$, the coefficients un, 
$n\geq 2M+1$, are all uniquely determined by the recurrence relation (2.3). To finish the proof of the lemma, it remains to be shown that the formal Laurent series is convergent on a small enough punctured neighbourhood of 
$t=t_0$. That is, there exists an 
$R\geq 1$ large enough such that the coefficients un, 
$n\geq 0$, of the Laurent series solution are bounded in modulus as
\begin{equation}
|u_n|\leq R^n\qquad (n\geq 0).
\end{equation}This estimate is given in Appendix A. The lemma follows.
Since ODE (1.4) has order 2M, whereas the expansions around poles in the above lemma only have two free parameters, t 0 and 
$u_{2M}$, it might be expected that the remaining 
$2M-2$ free parameters enter through multi-valued, subdominant contributions. We investigate this in the next section.
3.3. Dominant balance analysis and branch-point singularities
Let 
$\hat{u}(t)$ denote a Laurent series solution as in (3.3), for some choices of 
$t_0,u_{2M}\in\mathbb{C}$ and choice of sign ±. We perturb the Laurent series solution 
$\hat{u}(t)$ around the pole t 0 by adding a lower-order contribution y(t),
\begin{equation*}
u(t) = \hat{u}(t) + y(t).
\end{equation*} By taking the difference of two copies of the ODE (1.4), one with dependent variable u and the other with 
$\hat{u}$, we arrive at the following ODE for y,
\begin{equation}
y^3 + 3\,\hat{u}\,y^2 + 3\hat{u}^2y = y + \sum_{k=1}^{M}\alpha_{2k}\,y^{(2k)}.
\end{equation} We now employ the method of dominant balance to simplify Eq. (3.7). The two dominant terms are the highest-order derivative and 
$3\hat{u}^2y$. Taking the dominant contribution of 
$\hat{u}$ close to the pole, we obtain a special case of Euler’s differential equation for the leading order of y,
\begin{equation}
(t-t_0)^{2M}y^{(2M)}(t)= (3M)_{2M}\,y(t).
\end{equation} To obtain its general solution, we use the ansatz 
$y(t) = (t-t_0)^{\lambda}$ where 
$\lambda\in\mathbb{C}$. This gives the necessary and sufficient condition
\begin{equation}
P_M(M+\lambda)= 0,
\end{equation} where 
$P_M(r)$ is the resonance polynomial defined in (2.5).
We thus see that the resonance polynomial controls the multi-valued subdominant contributions of the highest order. It thus forms the focal point of the analysis in the remainder of this section. Firstly, to ensure that solutions of (3.9) provide the general solution of (3.8), we prove in the following lemma that all the roots are simple.
Lemma 3. All the roots of the resonance polynomial 
$P_M(r)$ are simple.
Proof. The derivative of 
$P_M(r)$ is given by
\begin{align*}
P_M'(r)
&=\prod_{k=0}^{2M-1} (r-M-k) \sum_{k=0}^{2M-1} (r-M-k)^{-1}\\
&= \left( P_M(r)+(3M)_{2M} \right) \sum_{k=0}^{2M-1} (r-M-k)^{-1}.
\end{align*} Now let r be a root of PM. By Lemma 1, either 
$r=-1$, 
$r=4M$ or 
$\text{Im}(r) \neq 0$. In particular, since 
$r \notin [M,3M-1]\,\cap\,\mathbb{Z}$,
\begin{equation*}
P_M'(r) =(3M)_{2M} \sum_{k=0}^{2M-1} (r-M-k)^{-1}.
\end{equation*} If 
$r=-1$, then each term in the above sum is negative and hence 
$P_M'(-1) \lt 0$. On the other hand, if 
$r=4M$, then each term is positive and thus 
$P_M'(4M) \gt 0$. Finally, if 
$\text{Im}(r) \neq 0$, then
\begin{equation*}
\text{sgn}(\Im(r-M-k)^{-1}) = -\text{sgn}(\Im(r-M-k))= -\text{sgn}(\Im(r)),
\end{equation*} for 
$0\leq k\leq 2M-1$, so the imaginary part of each term has the same sign. It follows that 
$P_M'(r)$ is not real and thus also nonzero in this case.
According to the above lemma, the general solution of (3.8) is given by
\begin{equation}
y(t)=d_1(t-t_0)^{-M-1}+d_2(t-t_0)^{3M}+\sum_{1\leq j\leq 2M-2}{s_j(t-t_0)^{\lambda_j}},
\end{equation} where 
$d_1,d_2$ and the 
$s_j\in\mathbb{C}$ are free constants and λj, 
$1\leq j\leq 2M-2$, are the solutions of (3.9), other than 
$\lambda=-M-1$ and 
$\lambda=3M$. Note that d 2 plays the same role as the free coefficient 
$u_{2M}$ in Lemma 2. Also, since y is assumed subdominant to 
$\widehat{u}$ in the dominant balance analysis, we require 
$d_1=0$. More generally, only powers λ satisfying 
$\Re \lambda \gt -M$ are permitted. In terms of 
$r=\lambda+M$, the latter condition reads 
$\Re r \gt 0$. We have the following bounds for roots of the resonance polynomial.
Lemma. All of the roots r of the resonance polynomial 
$P_M(r)$, except 
$r=-1$ and 
$r=4M$, satisfy the bounds
\begin{equation*}
-1 \lt \Re r \lt 4M.
\end{equation*}Proof. From the defining equation of 
$P_M(r)$, (2.5), roots of the polynomial coincide with solutions of
\begin{equation*}
Q_M(r)=(3M)_{2M},\qquad Q_M(r):=(r-M)_{2M}.
\end{equation*} It will be convenient to apply a change of variables 
$r\mapsto x$, where 
$r=2M-\tfrac{1}{2}+x$, so that the equation becomes
\begin{equation}
\hat{Q}_M(x)=(3M)_{2M},
\end{equation}where
\begin{equation*}
\hat{Q}_M(x)=Q_M(2M-\tfrac{1}{2}+x)=\prod_{k=0}^{M-1}(x-(\tfrac{1}{2}+k))(x+(\tfrac{1}{2}+k)).
\end{equation*} This way, the solutions are centred around x = 0 since 
$\hat{Q}_M(x)$ is a symmetric polynomial.
All the solutions of (3.11) lie on the curve
\begin{equation}
\mathcal{C}=\{x\in\mathbb{C}:|\hat{Q}_M(x)|=(3M)_{2M}\}.
\end{equation} Geometrically, this is the curve consisting of points 
$x\in\mathbb{C}$ such that the product of the distances from x to the 2M real points
\begin{equation}
-M+\tfrac{1}{2}, -M+\tfrac{3}{2}, -M+\tfrac{5}{2},\ldots,M-\tfrac{5}{2}, M-\tfrac{3}{2}, M-\tfrac{1}{2},
\end{equation} is constant and equal to 
$(3M)_{2M}$. A plot of 
$\mathcal{C}$ and the above points is given in Figure 1 for M = 5.
Figure 1. The curve 
$\mathcal{C}$ defined in (3.12) is shown in blue in the complex x-plane, with in black the solutions of (3.11) and red the points given in Eq. (3.13) for M = 5.
Note that 
$x_0=2M+\tfrac{1}{2}$ lies on the curve 
$\mathcal{C}$ as it is a solution to (3.11). Now, any point 
$x\in\mathbb{C}$, not equal to x 0, with 
$\Re x\geq x_0$, certainly lies farther away than x 0 from each of the points in (3.13). In particular, necessarily
\begin{equation*}|\hat{Q}_M(x)| \gt |\hat{Q}_M(x_0)|=(3M)_{2M}.\end{equation*} It follows that any point x on the curve 
$\mathcal{C}$ satisfies 
$\Re x\leq 2M+\tfrac{1}{2}$, and the only point on the curve where this inequality becomes an equality, is 
$x=2M+\tfrac{1}{2}$.
Furthermore, since 
$\hat{Q}_M(z)$ is an even polynomial, 
$\mathcal{C}$ is symmetric under reflection in the imaginary axis. Thus, we also obtain that each point x on the curve 
$\mathcal{C}$ satisfies 
$\Re x\geq -2M-\tfrac{1}{2}$, and the only point on the curve where this inequality becomes an equality, is 
$x=-2M-\tfrac{1}{2}$.
As 
$\mathcal{C}$ contains all the solutions to Eq. (3.1), it follows that any solution of the latter equation, apart from 
$x=\pm(2M+\tfrac{1}{2})$, satisfies
\begin{equation*}
-(2M+\tfrac{1}{2}) \lt \Re x \lt 2M+\tfrac{1}{2}.
\end{equation*} The lemma follows by translating these bounds to the roots of the polynomial 
$P_M(r)$ via 
$r=2M-\tfrac{1}{2}+x$.
As a corollary from the above lemma, we see that solutions of Eq. (3.9), other than 
$\lambda=-M-1$ and 
$\lambda=3M$, satisfy the bounds
\begin{equation*}
-M-1 \lt \Re\, \lambda \lt 3M.
\end{equation*}Now, solutions λ with
\begin{equation}
-M-1 \lt \Re\, \lambda\leq -M,
\end{equation}violate the assumptions of our dominant balance analysis. This suggests that such exponents would only appear in conjunction with other exponents through subdominant terms of the form
\begin{equation}
(t-t_0)^\Lambda,\quad \Lambda=j_0+j_1\lambda_1+j_2\lambda_2+\ldots+j_{2M-2}\lambda_{2M-2,}
\end{equation} where 
$j_k\in\mathbb{Z}_{\geq 0}$ for 
$0\leq k\leq 2M-2$, in the full asymptotic expansion of the solution y of (3.7) around t 0, see (Reference Goriely12, § 3.8).
Remark 3.
In the large M limit, the shape of the curve 
$\mathcal{C}$, defined in (3.12), is described, after a rescaling 
$x=M\, z$, by
\begin{equation}
\left\{z\in\mathbb{C}\,:\, \int_{-1}^1\log |z-t|\,dt=3\log 3-2\right\}.
\end{equation}This follows from the fact that
\begin{equation*}
\frac{1}{M}\log\frac{|\hat{Q}_M(M\,z)|}{M^{2M}}=\frac{1}{M}\sum_{k=0}^{2M-1}\log\left|z +1-\frac{2k+1}{2M}\right|\xrightarrow{M\rightarrow \infty} \int_{-1}^1\log |z-t|\,dt,
\end{equation*}and
\begin{equation*}
\frac{1}{M}\log\frac{(3M)_{2M}}{M^{2M}}=\frac{1}{M}(\log(3M)! -\log M!-2M\log M)\xrightarrow{M\rightarrow \infty} 3\log 3-2,
\end{equation*}by Stirling’s formula.
We conclude from our dominant balance analysis that the general singularities of solutions to ODE (1.4), are at leading order poles of order M, with subdominant multi-valued contributions given by complex powers with exponents that are roots of 
$P_M(M+\lambda)$. Since pole singularities, as classified in Lemma 2, only have two free parameters, they are a rare occurrence among general singularities when M > 1. We remark that, nonetheless, for special values of the parameters, there exist global solutions of ODE (1.4) on the complex plane that admit finite polynomial expressions in terms of hyperbolic secants [Reference Qiang, Alexander and de Sterke21], and thus only have pole singularities, and infinitely many of them.
4. Numerical analysis of the pure quartic solitary wave
In this section, we discuss the pure quartic solitary wave, the shape of which is governed by the fourth-order ODE
\begin{equation}
u^3=4\,u+u^{(4)},
\end{equation} where 
$u=u(t)$. Analogous to (3.1), this ODE has an integral of motion
\begin{equation}
\mathcal{I}(u):=2\,u^2-\tfrac{1}{4}u^4-\tfrac{1}{2}(u'')^2+u' u'''.
\end{equation}By applying the scaling
\begin{equation*}
u\mapsto \tilde{u}=2\, u,\quad t\mapsto \tilde{t}=30^{-\frac{1}{4}}t,
\end{equation*}the ODE is normalized as in (1.4), namely
\begin{equation*}
\tilde{u}^3=\tilde{u}+\alpha_4 \, {\tilde{u}}^{(4)},\qquad \alpha_4=\frac{1}{120}.
\end{equation*}This scaling can be used to translate all the results in the previous sections to the standard form (4.1) that we will be using here.
Our focus in this section is on the fundamental solution u(t) of Eq. (4.1). A plot of this solution is given in Figure 2. The fundamental solution can be characterized as the smooth, real, even solution on the real line with the smallest L 2-norm, which, together with all of its derivatives, vanishes exponentially as 
$t\rightarrow \pm \infty$. This solution has been studied intensively experimentally [Reference Runge, Hudson, Tam, de Sterke and Blanco-Redondo24], numerically [Reference Tam, Alexander, Blanco-Redondo and de Sterke26] and analytically [Reference Bandara, Giraldo, Broderick and Krauskopf5, Reference Qiang, Broderick and de Sterke22] on the real line. In this section, we carry out a numerical analysis of the fundamental solution in the complex plane and investigate the singularities closest to the origin. Remarkably, we find that even the fundamental solution has branch-point singularities.
Figure 2. Fundamental solitary wave solution of Eq. (4.1).
First, we consider the Cauchy data of this solution at t = 0,
\begin{equation*}
\mu_0=u(0),\quad \mu_1=u'(0),\quad \mu_2=u''(0),\quad \mu_3=u'''(0).
\end{equation*} Since u is even, we know that 
$\mu_1=\mu_3=0$. Therefore, Eq. (4.2) at t = 0 simplifies to
\begin{equation*}
\mathcal{I}(u)=2\,\mu_0^2-\tfrac{1}{4}\mu_0^4-\tfrac{1}{2}\mu_2^2.
\end{equation*} Furthermore, since u(t), as well as its derivatives, vanish exponentially as 
$t\rightarrow \pm \infty$, it follows that 
$\mathcal{I}(u)=0$. Therefore, the initial values 
$(\mu_0,\mu_2)$ lie on the plain algebraic curveFootnote 2
\begin{equation}
2\,\mu_2^2=\mu_0^2(8-\mu_0^2).
\end{equation}Using the analytic construction of the fundamental solution as a converging infinite series [Reference Hereman, Banerjee, Korpel, Assanto, van Immerzeele and Meerpoel13, Reference Qiang, Broderick and de Sterke22], we can obtain the values of the Cauchy data at t = 0 of the fundamental solution with arbitrary numerical precision. By evaluating the convergent series, truncated at the 28th exponential term, i.e. j = 27 in [Reference Qiang, Broderick and de Sterke22, Eq. (8)], at t = 0 we obtain
\begin{align*}
\mu_0&=2.539452989387891732547416944293400,\\
\mu_2&=-2.236433917320019900927316157131181,
\end{align*}with 34 digits of significance. These values indeed satisfy (4.3) up to a numerical error of size 10−32.
Since u(t) is smooth on the real line and satisfies the ODE Eq. (4.1), it must be analytic on the real line (as follows, e.g. by application of [Reference Hille14, Thm 2.2.2]). However, as solutions can generally develop branch points at arbitrary locations in the complex plane, analytic continuation away from the real line is not guaranteed to be possible and the result is generally non-unique. By the analysis in Section 3, near a singularity 
$t=t_0$ of u(t), the solution is described by
\begin{equation*}
u(t)=\hat{u}(t)+y(t),
\end{equation*} where 
$\hat{u}(t)$ is single-valued with a double pole at 
$t=t_0$,
\begin{equation*}
\hat{u}(t)=\pm\left(\sqrt{120}(t-t_0)^{-2}+\tfrac{1}{3}\sqrt{\tfrac{2}{15}}(t-t_0)^2+u_4(t-t_0)^6\right)+\mathcal{O}\left((t-t_0)^{10}\right),
\end{equation*}and y(t) is a subdominant multi-valued contribution of the form
\begin{equation}
y(t)=s_1\,(t-t_0)^{\frac{3}{2}-\frac{\sqrt{71}}{2}i}+s_2\,(t-t_0)^{\frac{3}{2}+\frac{\sqrt{71}}{2}i}+\mathcal{O}\left((t-t_0)^{\frac{3}{2}+\delta}\right),
\end{equation} for some δ > 0, since 
$\lambda=\frac{3}{2}\pm\frac{\sqrt{71}}{2}i$ are the only solutions of Eq. (3.9) other than 
$\lambda=-2$ and λ = 6 when M = 2, see also Eq. (3.10). Here, the sign ± and the values of u 4, s 1 and s 2 are free. It can be convenient to write the leading order behaviour of y(t) as
\begin{equation*}
y(t)=(t-t_0)^{\frac{3}{2}}\Big[c_1\cos\Big(\tfrac{\sqrt{71}}{2}\log(t-t_0)\Big)+c_2\sin\Big(\tfrac{\sqrt{71}}{2}\log(t-t_0)\Big)\Big]+\mathcal{O}\left((t-t_0)^{\frac{3}{2}+\delta}\right),
\end{equation*} with 
$s_1=c_1+c_2\,i$ and 
$s_2=c_1-c_2\,i$.
We investigate the singularities of the fundamental solution u(t) closest to the origin and numerically determine the values of the free constants in the above formulas for them. To determine their location, we consider the power series expansion of u(t) around t = 0,
\begin{equation*}
u(t)=\sum_{n=0}^\infty a_{n}t^{2n},\quad a_0=\mu_0,\quad a_1=\tfrac{1}{2}\mu_2.
\end{equation*}The coefficients are determined by the recursion
\begin{equation*}
(2n+4)(2n+3)(2n+2)(2n+1)a_{n+2}=\sum_{(i,j)\in A_n}a_ia_ja_{n-i-j}+(3\mu_0-4)a_n,
\end{equation*}where An is the set defined in Eq. (2.4). Recursively computing the values of the coefficients yields that, for n = 2000,
\begin{equation}
|a_{n}|^{-\frac{1}{2n}}\approx 2.74.
\end{equation} This suggests that u(t) has one or more singularities at a distance of approximately 2.74 from the origin in the complex plane. Furthermore, since 
$a_0 \gt 0$, 
$a_1 \lt 0$ and 
$3\mu_0-4 \gt 0$, it follows from the recursion that 
$(-1)^n a_{n} \gt 0$, for all 
$n\geq 0$. By the Vivanti–Pringsheim theorem, u(t) must therefore have singularities at 
$t=\pm i s$, where s the radius of convergence of u(t).
To determine the value of s more precisely, we consider the solution u(t) on the imaginary axis. We set
\begin{equation*}
v(x)=u(i x),\quad t=ix,
\end{equation*}and consider x as a real variable. Then, v(x) is also a solution of the quartic ODE (4.1), which is real-valued for real x around x = 0, with Cauchy data
\begin{equation*}
v(0)=\mu_0,\quad v'(0)=0,\quad v''(0)=-\mu_2,\quad v'''(0)=0.
\end{equation*}Using Mathematica’s inbuilt NDSolve, we find that v(x) develops a singularity at
\begin{equation}
s=2.747328905704978057791970142220888,
\end{equation} rounded off at the 34th decimal point. This means that singularities of u(t), closest to the origin, are located at 
$t=\pm i s$. We see numerically that 
$v(x)\sim \sqrt{120}(x-s)^{-2}$ as 
$x\uparrow s$. The difference
\begin{equation}
y(x):=v(x)-\Big(\sqrt{120}(x-s)^{-2}+\tfrac{1}{3}\sqrt{\tfrac{2}{15}}(x-s)^2\Big),
\end{equation}is shown in Figure 3.
According to our analysis, see Eq. (4.4), y(x) should vanish with order 
$\tfrac{3}{2}$ as 
$x\uparrow s$, which is consistent with the numerics in Figure 3. For a more precise numerical comparison, we set
\begin{equation}
h(x):=(s-x)^{-\frac{3}{2}}y(x).
\end{equation}This function is plotted in Figure 4.
Figure 4. The functions h(x) (solid blue) and 
$h_{\text{mod}}(x)$ (dashed red), defined in equations (4.8) and (4.9), on the interval 
$(1.5,s)$, where t = ix. Here is is the location of the singularity of u(t) closest to the origin with 
$\Im t\geq 0$. The numerical value of s is given in (4.6).
Then, again according to our analysis, h(x) should be asymptotic to
\begin{equation}
h_\text{mod}(x)=c_1\cos\Big(\tfrac{\sqrt{71}}{2}\log(s-x)\Big)+c_2\sin\Big(\tfrac{\sqrt{71}}{2}\log(s-x)\Big),
\end{equation}for appropriate values of the constants c 1 and c 2. A numerical fit in Mathematica yields
\begin{equation*}
c_1\approx -0.00179,\qquad c_2\approx 0.00092,
\end{equation*} for which we get excellent asymptotic agreement between h(x) and 
$h_\text{mod}(x)$ as 
$x\uparrow s$, as shown in Figure 4.
Returning to the original fundamental solution u(t), we draw the following conclusions. The singularities closest to the origin are purely imaginary and lie at 
$t=\pm i s$, with the numerical value of s given in Eq. (4.6). Around the singularity t = is, the local behaviour of u(t) is described by
\begin{equation*}
u(ix)=\sqrt{120}(s-x)^{-2}+(s-x)^{\frac{3}{2}}h_\text{mod}(x)+\tfrac{1}{3}\sqrt{\tfrac{2}{15}}(s-x)^2+\mathcal{O}\left((s-x)^\delta\right),
\end{equation*} as 
$x\rightarrow s$ with 
$|\arg (x-s)| \lt \pi$, where 
$h_\text{mod}(x)$ is the multi-valued function defined in Eq. (4.9). Here δ > 2 and we expect the estimate to hold for 
$\delta=\tfrac{7}{2}$.
These numerical results show that the fundamental solution, corresponding to the pure quartic solitary wave, has movable branch-point singularitities in the complex plane. In particular, this demonstrates that the non-Painlevé integrability of ODE (1.4) manifests itself even in the fundamental solution. It further suggests that the fundamental solution in the pure quartic case does not admit a simple exact description like the Karlsson and Höök solution [Reference Karlsson and Höök16].
5. Discussion
The generalization of the nonlinear Schrödinger equation we have explored here is characterized by the inclusion of high-order, even derivatives, corresponding to high orders of dispersion, whereas the nonlinear term has been left unchanged. These equations describe the pulses that can be generated using a setup that was recently reported [Reference Runge, Hudson, Tam, de Sterke and Blanco-Redondo24, Reference Runge, Qiang, Alexander, Rafat, Hudson, Blanco-Redondo and de Sterke25]. We note that other hierarchies based on the nonlinear Schrödinger equation, with all members integrable, have been explored by Ankiewicz et al. [Reference Ankiewicz, Kedziora, Chowdury, Bandelow and Akhmediev4] and by Bandelow et al. [Reference Bandelow, Ankiewicz, Amiranashvili and Akhmediev6]. However, in these equations most of the terms are nonlinear and include high-order derivatives. In the context of the experiments by Blanco-Redondo et al. [Reference Blanco-Redondo, de Sterke, Sipe, Krauss, Eggleton and Husko8] and Runge et al. [Reference Runge, Hudson, Tam, de Sterke and Blanco-Redondo24, Reference Runge, Qiang, Alexander, Rafat, Hudson, Blanco-Redondo and de Sterke25], such terms have no obvious physical significance.
Our work in this paper relied on the analytic method developed earlier [Reference Qiang, Broderick and de Sterke22] to find analytic expressions, in terms of a rapidly converging sum of functions, for the solitary wave solutions of the generalized nonlinear Schrödinger equation. We showed that this method can be used to explore these solutions to very high accuracy in the complex plane. This is illustrated by the remarkable fact that we could find the position of the branch point of Eq. (4.1) to over 30 decimal places. This represents a significant advantage over purely numerical methods. Additional work done since the work discussed here, has shown that this analytic method can be extended beyond the fundamental solution, to multi-peak solutions as well. The modifications to the analytic method allowing for such solutions to be found are shown in a paper currently in preparation.
6. Conclusion
Recent experimental work has demonstrated the ability to generate the solutions to generalizations of the nonlinear Schrödinger equation to higher orders of dispersion [Reference Blanco-Redondo, de Sterke, Sipe, Krauss, Eggleton and Husko8, Reference Runge, Hudson, Tam, de Sterke and Blanco-Redondo24, Reference Runge, Qiang, Alexander, Rafat, Hudson, Blanco-Redondo and de Sterke25]. Although the properties of some of the solutions have been explored, for example their functional forms [Reference Kudryashov18, Reference Qiang, Alexander and de Sterke21, Reference Qiang, Broderick and de Sterke22] and their stability [Reference Tam, Alexander, Blanco-Redondo and de Sterke26], see also [Reference Curran and Marangell10], this work has been applied to each order individually. Therefore, until now overarching statements regarding the properties of the entire hierarchy Eq. (1.2) of the resulting equations have been missing. In this paper we present the first such general results: we demonstrate that none of them are integrable in the Painlevé sense. This conclusion applies irrespective of the order of the differential equation and irrespective of the presence of possible lower order derivatives. We have further shown that the non-integrability manifests itself in physically relevant solutions like the pure quartic solitary wave, when extended into the complex plane.
Data availability statement
For access to the data and materials produced and used in this paper, please contact the corresponding author Dr Pieter Roffelsen.
Acknowledgements
The authors thank Andrew Hone, Nalini Joshi and Yang Shi for fruitful discussions.
Author contributions
Conceptualization: P.R., C.M.dS. and N.G.R.B. Methodology: P.R. Formal analysis: P.R., P.L. and J.J.M-B. Supervision: P.R., C.M.dS., N.G.R.B. and Y.L.Q. Validation: Y.L.Q. Visualization: P.R., P.L., J.J.M-B., N.G.R.B., Y.L.Q. and C.M.dS. Writing original draft: P.R., P.L. J.J.M-B. All authors reviewed and edited the final version manuscript.
Funding statement
P.R. was supported by Australian Research Council Discovery Project #DP210100129. C.M.deS. was supported by an ARC Discovery Project (DP230102200) and by the Australian Research Council Centre of Excellence in Optical Microcombs for Breakthrough Science (project CE230100006) and funded by the Australian Government.
Competing interests
None.
Ethical standards
The research meets all ethical guidelines, including adherence to the legal requirements of the study country.
Appendix A. An estimate of coefficients
In this Appendix we derive estimate (3.6). Our starting point to obtain this estimate, is the following estimate that we obtain from recursion (2.3),
\begin{equation}
|\alpha_M| |P_M(2n)| |u_n|\leq \sum_{(i,j)\in A_n}|u_iu_ju_{n-i-j}|+\sum_{k=0}^{M-1}|u_{k+n-M}|
|\alpha_k| \prod_{j=1}^{2k}|2n-3M+j|,
\end{equation}for all 
$n\geq 0$.
Since 
$P_M(2n)\sim (2n)^{2M}$ as 
$n\rightarrow \infty$, we can find an 
$N_1\geq 2M$ such that, for all 
$n\geq N_1$,
\begin{equation*}
|P_M(2n)|\geq \tfrac{1}{2} (2n)^{2M}.
\end{equation*} Then, for 
$n\geq N_1$, we infer the following estimate from (A.1),
\begin{align*}
|u_n| & \leq \frac{2}{|\alpha_M|(2n)^{2M}}\sum_{(i,j)\in A_n}|u_iu_ju_{n-i-j}|\\
&\quad +\sum_{k=0}^{M-1}|u_{k+n-M}|\,
\left|\frac{\alpha_k}{\alpha_M}\right| \frac{2}{(2n)^{2M-2k}}\prod_{j=1}^{2k}\left|\frac{2n-3M+j}{2n}\right|.
\end{align*} Next, we choose an 
$N_2\geq N_1$ such that, for all 
$n\geq N_2$ and for all 
$0\leq k\leq M-1$,
\begin{equation*}
\prod_{j=1}^{2k}\left|\frac{2n-3M+j}{2n}\right|\leq 2.
\end{equation*} Then, for 
$n\geq N_2$, we have the estimate
\begin{equation}
|u_n|\leq \frac{2}{|\alpha_M|(2n)^{2M}}\sum_{(i,j)\in A_n}|u_iu_ju_{n-i-j}|+\sum_{k=0}^{M-1}|u_{k+n-M}|\,
\left|\frac{\alpha_k}{\alpha_M}\right|.
\end{equation} Finally, we choose 
$N\geq N_2$, such that, for all 
$n\geq N$,
\begin{equation}
\frac{2}{|\alpha_M|(2n)^{2M}} (\tfrac{1}{2}(n+1)(n+2)-3)\leq \frac{2}{3},
\end{equation}and we set
\begin{equation}
R=\max\left(R_0,1,|u_1|,|u_2|^{\frac{1}{2}},|u_3|^{\frac{1}{3}},\ldots, |u_{N}|^{\frac{1}{N}}\right), \qquad R_0:=3\sum_{k=0}^{M-1}\left|\frac{\alpha_k}{\alpha_M}\right|.
\end{equation} Note that, by definition of R, 
$|u_n|\leq R^n$ for all 
$0\leq n\leq N$. To prove the same inequality for n > N, we use estimate (A.2) and induction. So, suppose n > N is such that 
$|u_m|\leq R^m$ holds for all 
$0\leq m\leq n-1$. Then we obtain, from estimate (A.2),
\begin{align*}
|u_n|
&\leq \frac{2}{|\alpha_M|(2n)^{2M}}\sum_{(i,j)\in A_n}R^n+\sum_{k=0}^{M-1}R^{k+n-M}\,
\left|\frac{\alpha_k}{\alpha_M}\right|\\
&\leq \frac{2}{|\alpha_M|(2n)^{2M}}\#(A_n)R^n+R^{n-1}\sum_{k=0}^{M-1}
\left|\frac{\alpha_k}{\alpha_M}\right|\\
&=\frac{2}{|\alpha_M|(2n)^{2M}}(\tfrac{1}{2}(n+1)(n+2)-3)R^n+R^{n-1}(\tfrac{1}{3}R_0)\\
&\leq \tfrac{2}{3}R^n+\tfrac{1}{3}R^n =R^n,
\end{align*}where, to obtain the third line, we used that the cardinality of An is 
$\#(A_n)=\tfrac{1}{2}(n+1)(n+2)-3$ as well as the definition of R 0 in (A.4). The fourth line was obtained by application of inequality (A.3) and the fact that 
$R\geq R_0$. By induction, we obtain estimate (3.6) for all 
$n\geq 0$.
Open access
