This post presents an elementary but highly general demonstration that avoids any considerations of Bessel functions. (You need to know some basic derivatives of real functions of one variable, as presented in the first section, and only the simplest properties of integrals -- linearity and positivity--; no integrals will need to be evaluated.)
Calculus preliminaries
Observe that the function defined on $(0,\infty)$ (and extended to $0$ via continuity) as
$$u(x) = x\log x$$
is strictly convex because $u^{\prime\prime}(x) = 1/x$ is everywhere positive. Because whenever $x\gt 0$ also $1+x\gt 0,$ Jensen's Inequality implies
$$\frac{1+x}{2} \log\left(\frac{1+x}{2}\right) = u\left(\frac{1+x}{2}\right) \lt \frac{1}{2}u(1) + \frac{1}{2}u(x) = \frac{1}{2}u(x).$$
Next, $u(x) \ge x-1$ for all $x \ge 0$ because
$$\frac{\mathrm d}{\mathrm d x}(u(x) - (x-1)) = \log x$$
has a unique zero at $x=1,$ which must be the point of the global minimum, and
$$u(1) - (1-1) = 1\log(1) - 0 = 0$$
is its minimal value, QED.
A general result
Let $\Omega$ be a measure space and $\mu$ be a finite measure on $\Omega.$ Divide $\mu$ by $\mu(\Omega)$ to normalize it.
Suppose that $f$ is a density function with respect to $\mu$ for which $\int_\Omega f\,\mathrm d\mu = 1:$ that is, $f$ is a probability density. We are interested in the functional
$$H[f] = \int_\Omega f \log f \,\mathrm d\mu,$$
the negative of the entropy. The preceding inequality shows that $H$ is nonnegative since
$$H[f] = H[f] + 1 - 1 = H[f] + \int_\Omega (f\log f - f + 1)\,\mathrm d\mu\ge \int_\Omega 0 \,\mathrm d\mu = 0.$$
Moreover, the constant function $1$ attains this minimum because
$$H[1] = \int_\Omega 1 \log 1\,\mathrm d\mu = \int_\Omega 0 \,\mathrm d\mu = 0.$$
What was the effect of the initial normalization of $\mu$? It did two things:
All integrals were multiplied by $1/\mu(\Omega).$
All density functions had to be multiplied by $\mu(\Omega)$ to compensate for (1).
The first effect will not change the fact that $0$ is the lower bound of $H,$ while the second effect changes $H$ to
$$\begin{aligned}
H[\mu(\Omega)f] &= \int_\Omega \mu(\Omega) f \log(\mu(\Omega) f)\,\mathrm d\mu\\
&= \mu(\Omega) \int_\Omega f \log f\,\mathrm d\mu + \mu(\Omega)\log(\mu(\Omega))\int_\Omega f\,\mathrm d\mu\\
&= \mu(\Omega) H[f] + \mu(\Omega)\log(\mu(\Omega))
\end{aligned}$$
because $\int_\Omega f\,\mathrm d\mu = 1.$ This constant additive change will not alter the order relationship between $H[1]$ and $H[f]$ -- it only changes the value of the lower bound.
(It is instructive to compare this to the analysis at https://stats.stackexchange.com/a/415436/919 which concerns probability distributions on a non-finite measure space. There, the analog of changing the measure is to rescale the variable.)
We conclude
On a space $\Omega$ of finite measure, the entropy $-H[f]$ is maximized among probability densities $f$ for the constant probability density $f = 1/\mu(\Omega)$ (attaining the minimum value $\mu(\Omega)\log\mu(\Omega)$).
Solution to the problem
A von-Mises Fisher distribution in $D$ (Euclidean) dimensions is an isotropic Normal distribution conditioned on the unit sphere. (The parameter $1/\kappa$ uniformly scales that distribution.) Applying the preceding result to the induced Lebesgue measure $\mu$ on the sphere shows that the entropy is maximized by the uniform distribution on the sphere. That is the distribution corresponding to $\kappa = 0,$ QED.
