Why does $d^n \exp(-x-x^{-1})/(dx)^n$ only have $n$ positive real zeroes?

Question

Set $f(x) = \exp(-x-x^{-1})$. An easy induction shows that $$\frac{d^n}{(dx)^n} f(x) = \phi_n(x^{-1}) f(x)$$ for $\phi_n$ a polynomial of degree $2n$. Clearly, the roots of $\phi_n(x^{-1})$ are the same as the roots of $f^{(n)}(x)$ and, by repeated use of Rolle's theorem, $f^{(n)}(x)$ has at least $n$ roots in $(0, \infty)$. Since $\phi_n$ has degree $2n$, there is room for many more positive real roots than that.

However, for $n \leq 50$, computer computations show that $\phi_n$ has exactly $n$ positive real roots! Why?

Motivation: Nothing really, I was just thinking about this question and fooling around.

Data: Here are the first $10$ values of $\phi_n(y)$:

1, -1 + y^2, 1 - 2 y^2 - 2 y^3 + y^4, -1 + 3 y^2 + 6 y^3 + 3 y^4 - 6 y^5 + y^6, 1 - 4 y^2 - 12 y^3 - 18 y^4 + 32 y^6 - 12 y^7 + y^8, -1 + 5 y^2 + 20 y^3 + 50 y^4 + 60 y^5 - 50 y^6 - 180 y^7 + 115 y^8 - 20 y^9 + y^10, 1 - 6 y^2 - 30 y^3 - 105 y^4 - 240 y^5 - 200 y^6 + 540 y^7 + 1095 y^8 - 1080 y^9 + 294 y^10 - 30 y^11 + y^12, -1 + 7 y^2 + 42 y^3 + 189 y^4 + 630 y^5 + 1295 y^6 + 420 y^7 - 5075 y^8 - 7140 y^9 + 10521 y^10 - 3990 y^11 + 623 y^12 - 42 y^13 + y^14, 1 - 8 y^2 - 56 y^3 - 308 y^4 - 1344 y^5 - 4256 y^6 - 7560 y^7 + 3430 y^8 + 48160 y^9 + 48664 y^10 - 108360 y^11 + 53788 y^12 - 11424 y^13 + 1168 y^14 - 56 y^15 + y^16, -1 + 9 y^2 + 72 y^3 + 468 y^4 + 2520 y^5 + 10668 y^6 + 31752 y^7 + 45234 y^8 - 83160 y^9 - 478674 y^10 - 330120 y^11 + 1186836 y^12 - 742392 y^13 + 201132 y^14 - 27720 y^15 + 2007 y^16 - 72 y^17 + y^18, 1 - 10 y^2 - 90 y^3 - 675 y^4 - 4320 y^5 - 22800 y^6 - 93240 y^7 - 256830 y^8 - 246960 y^9 + 1272348 y^10 + 5033700 y^11 + 1965810 y^12 - 13829760 y^13 + 10636800 y^14 - 3530520 y^15 + 614925 y^16 - 59760 y^17 + 3230 y^18 - 90 y^19 + y^20

The other roots are shrinking towards $0$ like $1/n$ while accumulating on some sort of curve. Here is a picture for $n=50$:

Answer 1

As I said in the remarks, the old post on MathStackExchange is highly relevant here to the extent that it solves the problem. The only thing you need is to turn the argument there upside down.

For $a>0$, we have the identity $$ \frac{d^{l+1}}{dx^{l+1}}[(x+a)^{-(m-1)}e^{-1/x}]=(x+a)^{-l}\frac d{dx}\left((x+a)^{l+1}\frac{d^{l}}{dx^{l}}[(x+a)^{-m}e^{-1/x}]\right) $$ This is just the identity $$ \frac{d^{l+1}}{dx^{l+1}}[(x+a)F(x)]=(x+a)^{-l}\frac d{dx}\left((x+a)^{l+1}\frac{d^{l}}{dx^{l}}F(x)\right) $$ aplied with $F(x)=(x+a)^{-m}e^{-1/x}$. It implies by Rolle that the $l+1$-st derivative of $(x+a)^{-(m-1)}e^{-1/x}$ has at least as many positive roots as the $l$-th derivative of $(x+a)^{-m}e^{-1/x}$ for $m=1$ and at least one more positive root for $m\ge 2$ (this difference comes from the fact that the product $(x+a)^{l+1}\frac{d^{l}}{dx^{l}}[(x+a)^{-m}e^{-1/x}]$ tends to $0$ at $+\infty$ for $m\ge 2$ but not for $m=1$. Inductively, we conclude that if $\frac{d^{l}}{dx^{l}}[(x+a)^{-m}e^{-1/x}]$ has at least $\ell+1$ positive root, then $\frac{d^{l+m}}{dx^{l+m}}[e^{-1/x}]$ has at least $l+m$ positive roots, which is absurd because it is a polynomial of degree $l+m-1$ times $x^{-2(l+m)}e^{-1/x}$.

Now we choose a huge $m$ and use $a=m$ to conclude that the $l$-th derivative of $\left(1+\frac xm\right)^{-m}e^{-1/x}$ has at most $l$ positive roots.

Letting $m\to+\infty$ with fixed $l$, we derive that the $l$-th derivative of $e^{-x-1/x}$ has at most $l$ ... well, not roots yet, but sign changes. Note however that the Rolle theorem spawns more than just roots. For analytic functions, it guarantees that between any two roots there is a sign change. Thus, if the $l$-th derivative has more than $l$ positive roots counted with multiplicity and $L$ is the highest multiplicity of a root, then the $l+L$-th derivative has more than $l+L$ sign changes. The end.

Answer 2

Here is an incomplete approach, which nevertheless explains the expected number of positive roots: For another variable $\alpha$ we define the more general family of polynomials $\phi_n(x,\alpha)\in\mathbb Q[\alpha][x]$ by \begin{equation} \phi_n(\frac{1}{x},\alpha)=e^{\alpha x+\frac{1}{x}}\cdot(e^{-\alpha x-\frac{1}{x}})^{(n)}. \end{equation} Set $\psi_n(x)=\phi_n(x,0)$. Note that $\phi_n(x,1)$ is the polynomial from the question. The idea is to use a deformation of $\phi_n(x,\alpha)$, where $\alpha$ runs from $1$ to $0$ on the real line.

First note that the $n$-th derivative of the product $e^{-\alpha x-\frac{1}{x}}=e^{-\alpha x}\cdot e^{-\frac{1}{x}}$ yields \begin{equation} \phi_n(x,\alpha)=\sum_{m=0}^n\binom{n}{m}(-\alpha)^{n-m}\psi_m(x). \end{equation} The polynomials $\psi_m(x)$ fulfill $\psi_0=1$ and $\psi_m=x^2(\psi_{m-1}-\psi_{m-1}')$. From that it follows easily that $\psi_m$ is monic of degree $2m$, and for $m\ge1$, it has the form \begin{equation} \psi_m(x)=(-1)^mm!x^{m+1}+\dots+x^{2m}. \end{equation} (The other coefficients, which we don't need here, can be computed explicitly too. The coefficient of $x^{m+1+k}$ in $\psi_m(x)$ is $\binom{m-1}{k}\frac{m!}{(k+1)!}(-1)^k$. In fact, $\psi_m(x)$ is, up to minor changes, the generalized Laguerre polynomial $L^{(1)}_{m-1}(x)$.)

We prove the question under the yet unproven assumption that $\phi_n(x,\alpha)$ is separable for all $0<\alpha\le1$. Since $\phi_n(0,\alpha)=(-\alpha)^n\ne0$, the number of positive real roots of $\phi_n(x,\alpha)$ is constant for $\alpha$ in this range.

From the formulas above, we see that the coefficient of $x^k$ in $\phi_n(x,\alpha)$ for $2\le k\le n+1$ is \begin{equation} -(-1)^nn(n-1)(n-1)\cdots(n-k+2)\alpha^{n+1-k}+\text{higher order terms in $\alpha$}, \end{equation} so they all have the same sign for $\alpha>0$ sufficiently small. Fix this $\alpha$. Thus among the smallest $n+2$ coefficients of $\phi_n(x,\alpha)$ there is exactly one sign change. The coefficients of $x^{n+1},\dots,x^{2n}$ contribute at most $n-1$ sign changes. So all together $\phi_n(x,\alpha)$ has at most $n$ sign changes. According to Descartes Rule, $\phi_n(x,\alpha)$ and then also $\phi_n(x,1)$ have at most $n$ positive roots.

Remark: One can also understand the behavior of the roots geometrically. The following images for $n=10$ and $n=11$ illustrate how the roots of $\phi_n(x,1)$ (the red ends) move to the roots of $\phi_n(x,0)$ (the blue ends). Only the two smallest positive roots of $\phi_n(x,\alpha)$ are shown:

When letting $\alpha$ run from $1$ to $0$, the following happens: The complex non-real roots of $\phi_n(x,1)$ walk to $0$, and so does the smallest positive root of $\phi_n(x,1)$ and the negative root (which exists for odd $n$).

Indeed, in the neighborhood of $\alpha=0$, those roots of $\phi_n(x,\alpha)$ which converge to $0$ behave as expected. This can be read off from a local expansion of these roots as Puiseux series in $\alpha$. The Newton polygon shows that they have $\alpha$-adic valuation $\frac{n}{n+1}$. Now $\alpha^{-n}\phi_n(x\alpha^{\frac{n}{n+1}},\alpha)$ is a polynomial in $\alpha^{\frac{n}{n+1}}$ and $x$, which upon setting $\alpha=0$ has the form $n!x^{n+1}-1$ (follows from the formulas from above). This shows that these $n+1$ roots leave $0$ at angles $2\pi k/(n+1)$ for $k=0,1,\dots,n$.

Answer 3

Update: The approach in (1) probably doesn't work either. Computing the Sturm chains numerically suggests that the inequality doesn't hold for $n=317$. Raising the precision still predicts that $2\phi_n(x)^2-x^3(\phi_n(x)\phi_n''(x)-\phi_n'^2(x))$ has $2$ positive roots.

A few small remarks which do not fit the field for comments:

(1) Maybe a minor modification of Fedor Petrov's suggestion from the comments could be made to work. It seems to be that $\frac{\phi_{n+1}(x)}{x^2\phi_n(x)}$ is increasing for $x>0$ on all the open intervals between the singularities. If that is the case, the claim follows easily. Using the easy relation \begin{equation} \phi_{n+1}(x)=(x^2-1)\phi_n(x)-x^2\phi_n'(x), \end{equation} we see that the monotonicity of this function is equivalent to \begin{equation} 2\phi_n(x)^2\ge x^3(\phi_n(x)\phi_n''(x)-\phi_n'^2(x)) \end{equation} for all $x>0$. Or more compactly, \begin{equation} \left(\frac{\phi_n'(x)}{\phi_n(x)}\right)'\le\frac{2}{x^3}. \end{equation} This inequality is rigorously proven using Sturm chains for all $n\le65$.

(2) The OP's expectation about the number of positive real roots holds for all $n\le6000$ by Descartes' Rule.

(3) I'm wondering if David's polynomials $\phi_n$ appear somewhere in the literature, or if they have a natural interpretation. Note that replacing $-x-\frac{1}{x}$ in their definition with $-x^2$ gives the Hermite polynomials. It is interesting that $x^2$ and $x+\frac{1}{x}$ are sort of the smallest non-constant rational functions which are fixed under the involutory automorphisms $x\mapsto-x$ and $x\mapsto\frac{1}{x}$, respectively. While these two automorphisms are conjugate in $\text{Aut}(\mathbb C(x))=\text{PGL}_2(\mathbb C)$, this does not seem to relate the $\phi_n$'s to the Hermite polynomials.

Answer 4

Another solution, inspired by fedja's. This one also establishes Peter Mueller's stronger statement that the coefficients of $\phi_n(x)$ have $n$ sign changes (in the sense of Descartes' Rule of Signs).

For any Laurent polynomial $f(x) \in \mathbb{R}[x, x^{-1}]$, let $D(f)$ denote the number of sign changes in the coefficients of $f$.

Claim 1 For any integers $0 \leq r \leq m$, $$D \left( \left( \frac{d}{dx} \right)^r (1-x/m)^m (1-x^{-1}/m)^m \right) \leq 2m-r.$$

Proof The only monomials with nonzero coefficient are $x^{-m-r}$, $x^{-m-r+1}$, ..., $x^{-r-1}$, $1$, $x$, ..., $x^{m-r}$, a total of $2m-r+1$ monomials. $\square$

Claim 2 For any Laurent polynomial $h(x)$ and any $a>0$, we have $D((x-a) h(x)) \geq D(h(x))+1$.

Proof This is a standard Lemma established in the course of proving Descartes' Rule of Signs. See, for example, this proof. $\square$

Claim 3 $$D \left(\frac{ \left( \frac{d}{dx} \right)^r (1-x/m)^m (1-x^{-1}/m)^m }{(1-x/m)^{m-r} (1-x^{-1}/m)^{m-r}} \right) \leq r.$$

This follows by $2m-2r$ uses of Claim 2.

Now, sending $m \to \infty$, we obtain $$D\left( \frac{\left( \frac{d}{dx} \right)^r \exp(-x-x^{-1})}{\exp(-x-x^{-1})} \right)\leq r$$ as desired.

Answer 5

I decided to try a contour integration attack. The formulas got too messy for me, but maybe someone else will finish it.

By Taylor's theorem, for $|z|<|x|$, we have $$f(x+z) = \sum_{n=0}^{\infty} \frac{z^n}{n!} f^{(n)}(x) = \left( \sum_{n=0}^{\infty} \frac{z^n \phi_n(1/x)}{n!} \right) f(x).$$ So $$\frac{\phi_n(1/x)}{n! (2 \pi i)} = \oint \frac{f(x+z)}{z^{n+1} f(x)} dz = \oint \exp{\Big (}-(x+z)^{-1} + x^{-1} - z - (n+1) \log z{\Big )} \ dz.$$ The integrand has three critical points, at the roots of $$z^3 + (n+2 x+1) z^2 + (2 n x+x^2+2 x-1) z + (n x^2+x^2)=0.$$ Working at the saddlepoint expansions at these roots was too messy for me, but maybe easy for someone else.

Answer 6

Not an answer, but the formulae are too bulky for a comment.

After a lot of experimenting I found a formula for the coefficients of $\phi_{n}(y)$. (Currently without formal proof, but the evidence is really strong.)

I conjecture that $\phi_{n}(y)$ has the following form $$ \phi_{n}(y)= (-1)^{n} + \sum_{k=2}^{2 n} c_{n,k} y^{k}, $$ where $$ c_{n,k}=c^{<}_{n,k}:=\frac{(-1)^{n+1} n!}{(n-k+1)!}\ {}_{2}F_{3}(1-\frac{k}{2},\frac{3}{2}-\frac{k}{2};2,2-k,n-k+2;4) $$ for $2\leq k \leq n $ and $$ c_{n,k}=c^{>}_{n,k}:=\frac{(-1)^{k} n!}{(k-n)!}{{n-1}\choose{k-n-1}} {}_{2}F_{3}(\frac{k}{2}-n,\frac{k}{2}-n+\frac{1}{2};1-n,k-n,k-n+1;4) $$ for $n< k \leq 2 n $. The ${}_{2}F_{3}$ are generalized hypergeometric functions.

The series of the hypergeometric functions in the formula for $c^{<}_{n,k}$ terminates, since one of the nominatorial parameters of the hypergeometric is always a non-positive integer.

The sums of the hypergeometric functions in the formulae for $c^{<}_{n,k}$ and $c^{>}_{n,k}$ terminate. The non-positive denominatorial parameters, are always overruled by one of the nominatorial parameters, since, e.g. $1-n<\frac{k}{2}-n$ for the relevant values. In cases of equality of the nominatorial and denominatorial parameters they cancel and the function collapses to a well defined ${}_{1}F_{2}$. This overruling and cancellation has to be carefully implemented in math algebra systems like Mathematica, since it is not easily recognized by them.

The sign of $c_{n,k}$ changes depending on $n$ and $k$ in irregular manner as already mentioned in other answers. The oscillations appear for $k$ in a relative small range above and below $n$ starting for $n$ above 7.

The formulae were found by using $f$ with a decoration of the variable $x$ $$ f(x)=\exp( - x \ z - \frac{y}{x} ) $$ and expand the result for $\phi$ in $x$, $z$, and $y$ using Mathematica. The formulae for the coefficients of the resulting multinomials were easily guessed and after setting $z=1$ and $y=1$ only Mathematica solvable sums remained.