Technical challenge: Limit of von-Mises distribution approximates normal. How to take the limit? Background:
In psychophysics or the study of ant navigation it's important to represent random variables on a circle. The most popular distribution for doing so is the von-Mises distribution (the wrapped Gaussian is too messy to work with). This is the von-Mises distribution for a random variable on a circle, $\theta \in [-\pi,\pi)$:
\begin{align} 
p(\theta;\mu,\kappa) = \frac{e^{\kappa cos(\theta - \mu)}}{\int_{-\pi}^\pi e^{\kappa cos(\theta)} \text{d}\theta} = \frac{e^{\kappa cos(\theta - \mu)}}{2 I_0(\kappa)}  ,
\end{align}
where $\kappa$ is the concentration parameter (have the inverse variance $\kappa \sim \sigma^{-2}$ in mind), $\mu$ the mean and the normalisation can be expressed in terms of the modified Bessel function of the first kind, $I_0(\kappa)$. 
Problem statement:
Wikipedia (https://en.wikipedia.org/wiki/Von_Mises_distribution) claims that
\begin{align} 
\lim_{\kappa \rightarrow \infty} p(\theta;\mu,\kappa) = \text{Normal}(\theta; \mu, \kappa^{-1/2}).
\end{align}
(The limit is not a rigorous mathematical limit but more a physicist-way of saying: expand the expression and throw higher-order terms away).
It makes sense: If the distribution is highly concentrated, we can forget about the circularity and just approximate the von-Mises as normal. Mathematically it makes intuitive sense, too. The argument of the exponential ranges between $\pm \kappa$, being 0 always at $\theta = \pm \pi/2$. At the same time, the normalisation constant grows large such that the whole distribution is essentially divided by a very large number and tends to zero everywhere except around it's mode. Around the mode, the cosine can be approximated by $\cos(\epsilon) = 1- \epsilon^2/2$, yielding a nice normal distribution. Despite this, I find it incredibly hard to show this! Can anyone help me with this technical challenge?
My (unsuccessful) approach:
I've tried to use Taylor expansion for the cosine and the following series representations (http://mhtlab.uwaterloo.ca/courses/me755/web_chap4.pdf) of the Bessel function for large $\kappa$:
\begin{align}
I_0(\kappa) = \sum_{n}^\infty \frac{(\kappa/2)^{2n}}{n!\Gamma(n+1)} \approx \frac{e^\kappa}{\sqrt{2\pi \kappa}}(1 + \frac{1}{8\kappa} + \mathcal{O}(\kappa^{-2})),
\end{align}
however, it didn't get me anywhere. Does anyone have any ideas how to attack this problem? Any help is much appreciated!

I'm not entirely sure how rigorous this is so check over my work carefully.
From the pdf $\frac{e^{\kappa\cos(x-\mu)}}{2\pi I_0(\kappa)}$
You already got that $I_0(\kappa)= \frac{e^{\kappa}}{\sqrt{2 \pi \kappa}}(1+\frac{1}{8\kappa})$
In the limit only angles close to the mean will be significant therefore you can assume that $x-\mu$ is small.
For small $z$ we have $cos(z)=1-\frac{z^2}{2}$
Substitute the approximations for $cos(x)$ and $I(\kappa)$ into the pdf to get
$ \frac{1}{2\pi I_0(\kappa)} e^{\kappa \left(1-\frac{(x-\mu)^2}{2}\right)}$
$ \frac{1}{2\pi \frac{e^{\kappa}}{\sqrt{2 \pi \kappa}}(1+\frac{1}{8\kappa})} e^{\kappa \left(1-\frac{(x-\mu)^2}{2}\right)}$
Expanding out the exponent in the numerator and canceling $e^{\kappa}$ terms gives:
$ \frac{1}{2\pi \frac{1}{\sqrt{2 \pi \kappa}}(1+\frac{1}{8\kappa})} e^{ \left(\frac{\kappa(x-\mu)^2}{2}\right)}$
Neglecting the $\frac{1}{8\kappa}$ term in the denominator and tidying up the constants gives you
$ \frac{1}{\sqrt{\frac{2 \pi}{\kappa}}} e^{ \left(\frac{\kappa(x-\mu)^2}{2}\right)}$
Which is a normal distribution
