Motivation: I want to solve this convolution problem on the circle: find $f$, given $g$ and $$ g(\theta) = \int_{S^1} e^{\mu\cos(\theta-\phi)}g(\phi)\ d\phi. $$ To do this, I want to find the Fourier coefficients $k_n$ of $K(\theta) = e^{\mu\cos(\theta)}$. (Here $\mu\in\mathbb{C}$ is nonzero.)
So far: To find the $n^{th}$ Fourier coefficient $k_n$, I integrated: \begin{align*} k_n &= \int_{S^1}e^{\mu\cos\theta}e^{-in\theta}\ d\theta\\ &= \int_{-\pi}^\pi e^{\mu\cos\theta}e^{-in\theta}\ d\theta \\ &= \int_{-\pi}^\pi e^{\mu\cos\theta}\cos(n\theta)\ d\theta - i\int_{-\pi}^\pi e^{\mu\cos\theta}\sin(n\theta)\ d\theta. \end{align*}
The integrand of the first term is even and a table of integrals tells me that $$\int_0^\pi e^{i\beta\cos(x)}\cos(nx)\ dx = i^n\pi J_n(\beta).$$ The integrand of the second term is odd, so it vanishes.
Thus I have $k_n = 2\pi i^n J_n(-i\mu)$. ($J_n$ is a Bessel function of first type.)
So far, so good. This all makes sense. However,
Confusion: Mathematica tells me (via the command FourierCoefficient[Exp[$\mu$ Cos[t]],t,n]) that the first several Fourier coefficients of $e^{\mu\cos\theta}$ are \begin{align*} k_0 &= I_0(\mu)\\ k_1 &= I_1(\mu)\\ k_2 &= I_2(\mu)\\ k_3 &= I_3(\mu)\\ k_4 &= \frac{1}{\mu^3}\bigg(\mu(24+\mu^2)I_0(\mu) - 8(6+\mu^2)I_2(\mu)\bigg)\\ k_5 &= \frac{1}{\mu^3}\bigg(\mu(48+\mu^2)I_1(\mu) - 12(16+\mu^2)I_2(\mu)\bigg)\\ k_6 &= \frac{1}{\mu^5}\bigg(\mu(1920+144\mu^2+\mu^4)I_0(\mu) - 6(640 + 128\mu^2+3\mu^4)I_1(\mu)\bigg)\\ k_7 &= \frac{1}{\mu^5}\bigg( \mu(5760 + 240\mu^2 + \mu^4)I_1(\mu) - 24(960+80\mu^2+\mu^4)I_2(\mu) \bigg)\\ k_8 &= \frac{1}{\mu^7}\bigg( \mu(322560 + 28800\mu^2 + 480\mu^4 + \mu^6)I_0(\mu) - 32(20160 + 4320\mu^2 + 150\mu^4 + \mu^6)I_1(\mu) \bigg)\\ &\vdots \end{align*} and they keep getting worse.
So what's going on here? On the one hand, integral tables and some common sense tell me that the Fourier coefficients of $K$ are the relatively nice $k_n = 2\pi i^n J_n(-i\mu)$. On the other hand, Mathematica is telling me that they are these nasty linear combinations of modified Bessel functions.
I know that the $I$-Bessel functions are related to the $J$-Bessel functions by $I_\alpha(x) = i^{-\alpha}(x)J_\alpha(ix)$, but that just tells me that $k_n = 2\pi i^n J_n(-i\mu) = 2\pi I_n(-\mu)$.
Am I making a dumb mistake? Is Mathematica using some recursion to express higher-order $I$-Bessel functions in terms of $I_0,I_1,I_2$?