# Is this function completely monotone?

### Background

Long ago I bumped into an exercise in ordinary differential equations, which asks to find a solution to the differential equation:

$$h'(x)=\frac{1}{2(1+xh(x))}$$

It turns out that $$h(x)$$ is the inverse function of

$$G(x)=2e^{x^2}\int_0^xe^{-u^2}\,du=\sqrt{\pi}e^{x^2}\hbox{erf}(x)$$

where $$\hbox{erf}(x)$$ is the standard error function. This function $$G(x)$$ (multiplied by $$1/2$$), appears in the problem set of the second chapter of the book "Special functions and their applications", by Lebedev, where the reader is asked to prove that it satisfies a certain differential equation, namely,

$$G'(x)=2(1+xG(x))$$

whence it easily follows that the inverse $$G^{-1}(x)=h(x)$$ satisfies the ODE above for $$h$$. In Lebedev's book the reader is asked to use the differential equation for $$G$$ to derive a series expansion for the error function, valid in the entire complex plain. In terms of $$G$$, the series is as follows:

$$G(z)=\sum_{k=0}^{\infty}\frac{2^kz^{2k+1}}{(2k+1)!!}\qquad (z\in\mathbb{C})$$

Lebedev leaves it at that, but some further comments may be of interest.

A simple calculation involving Stirling's formula shows that this series has an infinite radius of convergence, hence it defines an entire function in the complex plain. Moreover, restricting to positive real values only, the series expansion shows that all derivatives of $$G$$ are positive in $$(0,\infty)$$, i.e., $$G(x)$$ is an absolutely monotone function in $$(0,\infty)$$. Moreover, since $$G'(0)\neq 0$$, there is a disc around zero where the inverse function $$G^{-1}(z)=h(z)$$ is analytic, and its series expansion can be computed from that of $$G(z)$$. However, it turns out that the radius of convergence of the series for $$G^{-1}(z)$$ around zero is finite. This perhaps is connected to the fact that $$G(z)$$ is not one-to-one in the whole complex plain, i.e., does not have a global inverse.

### The problem

Looking into various properties of the inverse function $$h(x)=G^{-1}(x)$$, I came across the function $$g(x)=\exp(G^{-1}(x)^2)=\exp(h^2(x))$$, where $$x\geq 0$$.

Recall that a function $$f(x)$$ defined on $$(0,\infty)$$ is called completely monotone if it has derivatives of all order and $$(-1)^nf^{(n)}(x)\geq 0$$ for all $$n=0,1,2,\dots$$.

The problem is this:

Is the function $$g''(\sqrt{x})$$ completely monotone in $$[0,\infty)$$?

Using the differential equation for $$h$$ above, the second derivative is

$$g''(x)=\frac{\exp(h^2(x))}{2(1+xh(x))^3}$$ so $$g''(\sqrt{x})=\frac{\exp(h^2(\sqrt{x}))}{2(1+\sqrt{x}h(\sqrt{x}))^3}$$

Since this is a composition of $$g''(x)$$ and the square root function $$\sqrt{x}$$, and since the derivative of $$\sqrt{x}$$ is completely monotone, then if we knew that $$g''(x)$$ is itself completely monotone, we would be able to deduce that $$g''(\sqrt{x})$$ is completely monotone as well. But that is not the case: $$g''(x)$$ is not completely monotone, because, for example, it is not convex. Its second derivative, the fourth derivative of $$g(x)$$, is not a positive function.

### A partial result

With $$f(x)=g''(\sqrt{x})$$, the following holds. For every $$n$$, there is some $$x_n>0$$, such that for all $$x>x_n$$ we have $$(-1)^nf^{(n)}(x)>0$$.

My proof is neither elegant nor short. Moreover, I still don't know whether $$(-1)^nf^{(n)}(x)>0$$ holds for all $$x$$.

There is lots of information and literature involving completely monotone functions. For example, the property of being completely monotone is equivalent to the property of being the Laplace transform of a positive measure, and in our case, since $$g''(0)=\frac{1}{2}<\infty$$, we would have a bounded positive measure. There are lots of tricks and special cases that help determine whether a function is completely monotone, but as far as my research has reached, I could not settle the case for this specific function. It might be interesting to discover new methods of proving -- or disproving -- that a certain function is completely monotone.

Any insights, comments, remarks and proofs are welcome.

• This problem haunts me ! I am still working. May I confess that I am happy of the last edit ? Cheers. Mar 25 at 14:41

I do not know how much these preliminary results could help you.

$$x'=2 (h \,x+1)\implies x=c_1\, e^{h^2}+\sqrt{\pi }\,e^{h^2}\, \text{erf}(h)$$ For $$c_1=0$$, then $$h'=\frac 1{2 \left(1+\sqrt{\pi }\,h \, e^{h^2} \, \text{erf}(h)\right) } \qquad \qquad h''=-\frac{2h+\sqrt{\pi } \,e^{h^2} \left(2 h^2+1\right) \,\text{erf}(h)}{4 \left(1+\sqrt{\pi }\,h \, e^{h^2} \, \text{erf}(h)\right)^3}$$

The second deriative is minimum around $$h=0.341938$$ which corresponds to $$x=0.739762$$. For an approximation, expanded as a series around $$h=0$$, the numerator of $$h'''$$ is $$-4+24 h^2+\frac{224 }{3}h^4+\frac{1472}{15} h^6+O\left(h^8\right)$$ So, limited to $$O\left(h^6\right)$$, the approximate solution is $$h=\sqrt{\frac{\sqrt{249}-9}{56} }=0.347946$$

Using the expansion $$x=\sqrt \pi \, \sum_{n=0}^\infty \frac{h^{2 n+1}}{\Gamma \left(n+\frac{3}{2}\right)}$$ its reversion gives $$h=\frac x 2 \Bigg[1+ \sum_{n=1}^\infty (-1)^n\, a_n\, x^{2n} \Bigg]$$ where the first coefficients are $$\left\{\frac{1}{6},\frac{1}{15},\frac{29}{840},\frac{23}{1134},\frac{6409}{498960}, \frac{138661}{16216200},\frac{68868937}{11675664000},\frac{725772679}{173675502000},\cdots\right\}$$ which is decent up to $$x=2$$ (at this point $$h \sim 49$$).

Edit

For "large" values of $$x$$, consider the function $$f(h)=\sqrt{\pi }\,e^{h^2}\, \text{erf}(h)-x$$ Perform one single iteration of Halley method $$h_0=\sqrt{\log \left(\frac{x}{\sqrt{\pi }}\right)}\quad\implies\quad h_1= h_0 - \frac {2 f(h_0)\,f'(h_0)} {2 {[f'(h_0)]}^2 - f(h_0) \,f''(h_0)}$$

Some results $$\left( \begin{array}{ccc} x & \text{estimate} & \text{solution} \\ 2 & 0.721842253 & 0.708140254 \\ 3 & 0.876173756 & 0.876600650 \\ 4 & 0.993781108 & 0.994005821 \\ 5 & 1.082283625 & 1.082386926 \\ 6 & 1.152439394 & 1.152491989 \\ 7 & 1.210160261 & 1.210189706 \\ 8 & 1.258964788 & 1.258982566 \\ 9 & 1.301097077 & 1.301108472 \\ 10 & 1.338067681 & 1.338075343 \\ 15 & 1.474266175 & 1.474267862 \\ 20 & 1.565430010 & 1.565430596 \\ 25 & 1.633288003 & 1.633288263 \\ 30 & 1.687015375 & 1.687015510 \\ 35 & 1.731310000 & 1.731310078 \\ 40 & 1.768885355 & 1.768885404 \\ 45 & 1.801444560 & 1.801444592 \\ 50 & 1.830123798 & 1.830123820 \\ 55 & 1.855717315 & 1.855717330 \\ 60 & 1.878801252 & 1.878801264 \\ 65 & 1.899806289 & 1.899806298 \\ 70 & 1.919062481 & 1.919062488 \\ 75 & 1.936828124 & 1.936828130 \\ 80 & 1.953309008 & 1.953309013 \\ 85 & 1.968671643 & 1.968671646 \\ 90 & 1.983052585 & 1.983052585 \end{array} \right)$$

• Thank you for your answer; in fact, $h(x)$ is a function of $x$, not a constant. Mar 17 at 14:00
• @uniquesolution. I never considered that $h$ could be a constant. As soon as we write $h'$ we can analyze all its derivatives as a function of itself. Mar 18 at 12:09
• -- Oh, I see. Yes, I have a recursive formula for the coefficients of the Taylor expansion of $h$ around zero, namely: $h_{n+1}=-\sum_{k=0}^{n-1}{n\choose k}(n-k)h_{k+1}h_{n-k-1}$, where $h_j$ is the $j$'th derivative of $h$ at zero. There is also a recursive formula at non-zero points, but it is more complicated. In both cases, you do get some hold of the derivatives of $\exp(h^2(x))$, which is the function of interest, but the formulas are very complicated, and in particular, did not lead me to a solution of the problem. Mar 18 at 20:26
• By the way, by locating a singularity of the complex valued derivative, $h'(z)$, it is possible to compute the precise radius of convergence of the series for $h$, and it turns out to be the zero of the derivative of the Dawson integral. The number is approximately $0.924$ Mar 18 at 20:48
• just for the record -- this answer, albeit interesting, helps me not in solving the original problem. Mar 24 at 8:58