# Getting rid of sub-exponential terms in an asymptotic expansion for a modified Bessel function

I am trying to get Mathematica to produce suitable asymptotic expansions for some modified Bessel functions at large argument (more specifically, the expansion in the DLMF's eq. (10.40.1)), and I'm struggling with some subexponential terms which I would like to eliminate in a systematic fashion.

More to the point, suppose I try something of the form

Series[BesselI[n, z], {z, Infinity, 1}]


which is perfectly consistent with the documentation's description of the use of Series for asymptotic expansions, and which normally works perfectly well, but which in this occasion returns

$$e^{-z} \left(e^{2 z} \left(\frac{\sqrt{\frac{1}{z}}}{\sqrt{2 \pi }}+O\left(\left(\frac{1}{z}\right)^{3/2}\right)\right) +\left(\frac{i e^{i n \pi } \sqrt{\frac{1}{z}}}{\sqrt{2 \pi }}+O\left(\left(\frac{1}{z}\right)^{3/2}\right)\right)\right) \tag{1}.$$

This is in contrast with the asymptotic series in the DLMF,

$$\mathop{I_{\nu}}\nolimits\!\left(z\right)\sim\frac{e^{z}}{(2\pi z)^{\frac{1}% {2}}}\sum_{k=0}^{\infty}(-1)^{k}\frac{a_{k}(\nu)}{z^{k}}, \qquad\qquad |\mathrm{arg}(z)|<\pi/2 -\delta<\pi/2.$$

which contains no exponentially-decreasing term. I understand where this comes from, since, for example,

BesselI[1/2, z]


evaluates to $$\frac{\sqrt{\frac{2}{\pi }} \sinh (z)}{\sqrt{z}},$$ which obviously has both exponential and subexponential terms, with the terms in $e^{-z}$ playing a role for negative $z$, which Mathematica doesn't quite (yet) have a way to know isn't the case. (I'm also baffled as to why it's factorized it in the bizarre grouping of $(1)$, but that's relatively unimportant.)

I would like to get rid of these exponentially-decaying terms, in as systematic and general a way as possible. I have tried providing assumptions of various forms (e.g. $\mathrm{Re}(z)>1$, or $\arg(z)<\pi/4$, or $z>1$, and variations on that theme) with little success. I am currently using a scheme using Delete, but it is (i) a horrible hack, and (ii) liable to fail if Series returns its output in a different order than what the Delete construct is expecting.

I could also implement the DLMF series directly, but I was hoping that Mathematica is good enough at symbolic calculus that such a step shouldn't be necessary; and in any case I feel the problem is interesting and general enough to consider without recourse to that.

Is there an in-built, or at least a cleaner, way to get this expansion?

• Have you tried to use ReplaceAll with the rule Exp[a_*z] /; a<0 -> 0? – Marius Ladegård Meyer Mar 17 '16 at 15:26
• @Marius That just multiplies everything by zero. – Emilio Pisanty Mar 17 '16 at 15:58
• @Emilio Pisanty : Mathematica's expansion is correct. Compare with formula 10.40.5 of your cited reference in which this formula is said to be the "more general (and more accurate) version of (10.40.1)" – Dr. Wolfgang Hintze Mar 17 '16 at 16:46
• @Dr.WolfgangHintze I'm not saying it's incorrect, I'm saying it's not what I want. For the stated domain of DLMF (10.40.1), ($|\arg(z)|<\pi/2-\delta<\pi/2$), the additional terms do very little. In my application, $z$ is on or near the positive real axis ($\mathrm{Re}(z)\approx10$, $|\mathrm{Im}(z)|\lesssim1$), so each additional term in the $e^z$ series is ~10 times smaller than the previous one, but the leading term in the $e^{-z}$ series is ${\sim}e^{-20}{\sim}10^{-9}$ times smaller than the leading term in the $e^z$ series. Hence them being negligible and more a hindrance than help. – Emilio Pisanty Mar 17 '16 at 16:57
• Does this help: ComplexExpand[Series[BesselI[n, 2 z], {z, \[Infinity], 2}]] /. Exp[a_ z] /; a < 0 :> 0? It might be too "aggressive" in simply treating z as real. – chuy Mar 17 '16 at 22:02

Here's a more complete answer than was possible in the comment.

The series expansion in question is

s = Series[BesselI[n, 2 z], {z, \[Infinity], 2}] // Normal // Expand


$$-\frac{n^2 e^{2 z} \left(\frac{1}{z}\right)^{3/2}}{8 \sqrt{\pi }}+\frac{i n^2 \left(\frac{1}{z}\right)^{3/2} e^{-2 z+i \pi n}}{8 \sqrt{\pi }}-\frac{i \left(\frac{1}{z}\right)^{3/2} e^{-2 z+i \pi n}}{32 \sqrt{\pi }}+\frac{i \sqrt{\frac{1}{z}} e^{-2 z+i \pi n}}{2 \sqrt{\pi }}+\frac{e^{2 z} \left(\frac{1}{z}\right)^{3/2}}{32 \sqrt{\pi }}+\frac{e^{2 z} \sqrt{\frac{1}{z}}}{2 \sqrt{\pi }}$$

This Expansion is valid for all z in the compex plane. If, however, we are interested in the expansion for Re[z] going to infinity it simplifies matters if we remove all terms containing decaying exponents.

The appropriate method is pattern identification..

A fairly general pattern identifying negative exponentials is this

pat = c_.  Exp[b_. + (a_?Negative) z ];


You can select the terms with negative exponentials in s:

Cases[s, pat]


$$\left\{\frac{i \sqrt{\frac{1}{z}} e^{-2 z+i \pi n}}{2 \sqrt{\pi }},-\frac{i \left(\frac{1}{z}\right)^{3/2} e^{-2 z+i \pi n}}{32 \sqrt{\pi }},\frac{i n^2 \left(\frac{1}{z}\right)^{3/2} e^{-2 z+i \pi n}}{8 \sqrt{\pi }}\right\}$$

or those which are not negative eponentials

Cases[s, Except[pat] ]


$$\left\{\frac{e^{2 z} \sqrt{\frac{1}{z}}}{2 \sqrt{\pi }},\frac{e^{2 z} \left(\frac{1}{z}\right)^{3/2}}{32 \sqrt{\pi }},-\frac{n^2 e^{2 z} \left(\frac{1}{z}\right)^{3/2}}{8 \sqrt{\pi }}\right\}$$

The brutal method to remove the unwanted terms from s would be

sr = s /. pat -> 0


$$-\frac{n^2 e^{2 z} \left(\frac{1}{z}\right)^{3/2}}{8 \sqrt{\pi }}+\frac{e^{2 z} \left(\frac{1}{z}\right)^{3/2}}{32 \sqrt{\pi }}+\frac{e^{2 z} \sqrt{\frac{1}{z}}}{2 \sqrt{\pi }}$$

However, this obviously works only for sums. The Cases approach is more general.

Observation: the pattern pat contains default pattern which are composed of an underline followed by a dot. These are useful to designate a term which can be missing.

Try this

t = {Exp[-z/3], Exp[3 z]};

Cases[t, pat]

(* Out[151]= {E^(-z/3)} *)


I tried to solve this problem by explicitly adding the order of terms you want to drop, and simplifying the result, but ended up being more confused by the behavior of FullSimplify. Nevertheless,

Factor[Series[BesselI[n, z], {z, Infinity, 1}] + O[z, Infinity]^(1/2) Exp[-z]]


works, and you can get rid of O[_] with Normal if desired. Factor is needed to factor out Exp[-z], and add O[z, Infinity]^(1/2) to the unwanted series. Simplify can work, in place of Factor, but may not be robust for reasons below.

Technically, the answer is complete, but regarding why I was confused by FullSimplify, you will be dismayed to find it doesn't work in place of Factor, and simply returns

E^z Sqrt[O[1/z]]

Further investigation shows FullSimplify[O[z, Infinity]^(1/2) Exp[-z]] returns the same, suggesting Mathematica believes O[z, Infinity]^(1/2) Exp[-z] and O[z, Infinity]^(1/2) Exp[z] are of the same order, which I believe is false in general, if I go by this definition of O. The reason is $|f(z)|<\frac{A e^z}{z}$ does not imply there exists $A_-$ such that $|f(z)|<\frac{A_- e^{-z}}{z}$ (e.g. $f(z>1)=A$).

At first, I thought this was a quirk of FullSimplify, which professes to use some transformations which are only generically correct, and could have neglected to consider functions growing faster than O[z, Infinity]^(1/2) Exp[-z]. But I found Limit[O[z, Infinity]^(1/2) Exp[z], z -> Infinity] gives 0, so the lack of generality already exists in O.

Nevertheless, we can avoid this issue by using only Factor, and not any other transformations, which may convert O[z, Infinity]^(1/2) Exp[-z] to O[z, Infinity]^(1/2) Exp[z], from FullSimplify or even Simplify.

• That's an interesting trick for sure. I'm working on a separate part of the code at the moment but I'll need to return there in a few weeks and I'll give this a bit more testing then. – Emilio Pisanty Jul 21 '16 at 13:30