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

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?

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$:

• First one has same behavior for $n=10$, $20$ and $50$. (I feel like it should be equivalent to mine by a simple change of variables, but maybe not.) Second one has roughly $n$ roots, but not exactly: for $n=10$, $20$, $30$, $40$ and $50$, I get $10$, $22$, $34$, $48$, $60$ roots respectively. The corresponding polynomial now is degree $4n$. – David E Speyer Aug 19 '15 at 17:15
• The expected number of roots holds even up to $n=6000$ via Decartes' Rule. This is a little surprising as Descartes only gives an upper bound. While Faa di Bruno allows to compute messy expressions for the coefficients of $\phi_n$, I don't believe that one can decide the signs of them in order to apply Descartes. Experimentally, the first somewhat less than $n$ coefficients have constant sign, the somewhat less than the last $n$ coefficients alternate, and around the mid the behavior is somewhat irregular (I didn't see a good pattern). – Peter Mueller Aug 19 '15 at 22:25
• Differentiating the identity $x^2 f' = (1-x^2) f$ a total of $n$ times gives the recurrence $\phi_{n+1} = (x^2 - 2nx - 1) \phi_n + nx (-(n-1)x-2) \phi_{n-1} - n(n-1) x^2 \phi_{n-2}$ (I think). This can be combined with the identity $\phi_{n+1} = (x^2-1) \phi_n - x^2 \phi'_n$ to get a third order linear ODE for $\phi_n$, so the curve $(\phi_n(x), \phi'_n(x), \phi''_n(x))$ follows some explicit vector field in ${\bf R}^3$, which can be projected down to ${\bf RP}^2$. Perhaps the topology of this vector field is giving the claim? – Terry Tao Aug 21 '15 at 22:46
• @TerryTao That sounds like a really cool approach, which I'll think about when I get some time (unless someone else does it first). Just to check that I understand your strategy, that third order ODE has coefficients which depend on $x$, right? So we don't get a fixed vector field in $\mathbb{RP}^2$, but one that varies according to the "time" $x$. Or did I miss something? – David E Speyer Aug 21 '15 at 23:41
• Oh, good point, so it's still a 3D flow and so potentially rather complicated, unless one can use the asymptotic regime to simplify things by focusing on the top order terms. This approach, by the way, should also help explain the cardiod shape you observed for the other zeroes. – Terry Tao Aug 22 '15 at 5:44

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.

• Wow! That is a very cool argument! – David E Speyer Dec 22 '15 at 12:38

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 $$\phi_n(\frac{1}{x},\alpha)=e^{\alpha x+\frac{1}{x}}\cdot(e^{-\alpha x-\frac{1}{x}})^{(n)}.$$ 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 $$\phi_n(x,\alpha)=\sum_{m=0}^n\binom{n}{m}(-\alpha)^{n-m}\psi_m(x).$$ 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 $$\psi_m(x)=(-1)^mm!x^{m+1}+\dots+x^{2m}.$$ (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 $$-(-1)^nn(n-1)(n-1)\cdots(n-k+2)\alpha^{n+1-k}+\text{higher order terms in \alpha},$$ 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$.

• The discriminant of $\phi_n(x,\alpha)$, with respect to $x$, seems to be a polynomial in $\alpha$ with all coefficients of the same sign (or $0$). If one can prove that, this would be the missing piece. – Peter Mueller Aug 26 '15 at 10:40
• +1 fro the explanation and also for the figures. – Johannes Trost Dec 3 '15 at 8:36

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 $$\phi_{n+1}(x)=(x^2-1)\phi_n(x)-x^2\phi_n'(x),$$ we see that the monotonicity of this function is equivalent to $$2\phi_n(x)^2\ge x^3(\phi_n(x)\phi_n''(x)-\phi_n'^2(x))$$ for all $x>0$. Or more compactly, $$\left(\frac{\phi_n'(x)}{\phi_n(x)}\right)'\le\frac{2}{x^3}.$$ 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.

• For $n=317$. That's awesome. Of course, it suffices for each $n$ to find and exponent $a_n$ such that the function $t^{a_n}\phi_{n+1}(t)/\phi_n(t)$ increases... Though this already looks less perspective than initial question. – Fedor Petrov Aug 21 '15 at 19:35
• RE: (3) The $\phi_n$ appear - I don't think very helpfully, but do not read German - in an old paper of Raabe here. See II on page 95. – Benjamin Dickman Aug 23 '15 at 6:31
• Indeed, in order to approximately compute $\int_0^1\text{exp}(-x-1/x)dx$, Raabe (in 1836) computes approximately the four roots of $\phi_6(x)$ with $x>1$, and uses a theorem of Fourier (a refinement of Descartes for intervals) to show that there are no more roots. – Peter Mueller Aug 23 '15 at 8:42

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.

• In Claim 2: did you mean something like this: $D((x-a)h) /geq D(h)+1$ for nonzero $h$? Otherwise I don't understand how claim 2 implies claim 3. – Paata Ivanishvili Dec 23 '15 at 20:21
• Now fedja's solution and this argument eventually produced a proof from the book ... – Peter Mueller Dec 23 '15 at 20:37

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.

• Also by Taylor's theorem, the polynomials $\phi_n(y)$ have the generating function $$\sum_{n=0}^\infty \phi_n(y) \frac{z^n}{n!} = \exp\left(-z\left(1-\frac{y^2}{1+yz}\right)\right).$$ – Ira Gessel Dec 4 '15 at 15:58

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.

• Any chance you could use this to analyze the signs of the coefficients? Peter Mueller presented some excellent computational evidence that the Descartes bound is always $n$ (which would be an even stronger statement than the one I asked for.) – David E Speyer Dec 4 '15 at 15:49
• Peter wrote "Experimentally, the first somewhat less than n coefficients have constant sign, the somewhat less than the last n coefficients alternate, and around the mid the behavior is somewhat irregular (I didn't see a good pattern)." But your pattern above seems very simple. – David E Speyer Dec 4 '15 at 21:17
• @David Speyer You are right. Some of the hypergeometric functions for $k$ values around $n$ have negative sign, which results in a rather irregular pattern. I corrected my answer and deleted the incorrect comment. Thanks for the remark. – Johannes Trost Dec 5 '15 at 11:16