Consider Mathieu equation:
$$\frac{d^2}{d\xi^2}R(\xi)+(a-2q\cos(2\xi))R(\xi)=0.$$
It's a second order ODE, so it should have two linearly independent solutions. One of the choices is to denote one as Mathieu sine $\operatorname{S}(a,q,\xi)$ and another as Mathieu cosine $\operatorname{C}(a,q,\xi)$.
But many books instead take some strange analogy to Bessel functions. Namely, they say that two functions, cosine-elliptic $\operatorname{ce}_r(\xi,q)$ and sine-elliptic $\operatorname{se}_r(\xi,q)$ are solutions of the first kind analogous to one Bessel $\operatorname{J}_\nu(x)$ function and two another functions $\operatorname{fe}_r(\xi,q)$ and $\operatorname{ge}_r(\xi,q)$ are solutions of the second kind analogous to one Neumann $\operatorname{Y}_\nu(x)$ function. Actually, the analogy with Bessel functions seems to go only for modified aka radial Mathieu equation, but the number of functions is the same for angular equation.
But sine-elliptic $\operatorname{se}$ and cosine-elliptic $\operatorname{ce}$ are already linearly independent, why are they both in the first kind? And why is there another pair of linearly independent functions $\operatorname{fe}$ and $\operatorname{ge}$? What is so much different between Mathieu and Bessel equations that the former has twice the number of solutions than the latter?
My guess is that the pairs $\operatorname{se}$ and $\operatorname{fe}$ and similarly $\operatorname{ce}$ and $\operatorname{ge}$ correspond to $\operatorname{S}$ and $\operatorname{C}$ correspondingly, but have differing sets of characteristic exponents: the functions of the first kind have real characteristic exponents (corresponding to bands in solid-state physics) while those of the second kind have complex ones (corresponding to band gaps). Is this right?
Let's denote two most general linearly independent solutions of Mathieu equation with characteristic $a$ as $\operatorname{C}(a,q,x)$ and $\operatorname{S}(a,q,x)$. Now if we want the solution to be $2\pi n$-periodic, where $n\in\mathbb{N}$, we'll have to consider special forms of these solutions.
Even solution will require characteristic $a_r(q)$ where $r=k/n,\;k\in\mathbb{N}$ to satisfy boundary conditions and is usually denoted
$$\operatorname{ce}_r(x,q)=\operatorname{C}(a_r,q,x).$$
Odd solution will require characteristic $b_r(q)$ with the same constraints on $r$ with addition that $r>0$, and is usually denoted
$$\operatorname{se}_r(x,q)=\operatorname{S}(b_r,q,x).$$
These are usually normalized so that
$$\int_0^{2\pi n}(\operatorname{ce}_r(x))^2dx=\pi n$$
and the same for $\operatorname{se}$.
Now for each $a_r$ and $b_r$ there are also the "second" solutions, i.e. in addition to $\operatorname{ce}_r(x,q)$ we have $\operatorname{S}(a_r,q,x)$, and in addition to $\operatorname{se}_r(x,q)$ we also have $\operatorname{C}(b_r,q,x)$ as solutions to exact same corresponding equations. Actually, $a_r=b_r$ for any $r\not\in\mathbb{N}$. The problem is that as $r$ approaches some integer, these functions vanish due to normalization. Thus to recover these functions (which are non-periodic at integer $r$), the following functions are defined:
$$\operatorname{fe}_n(x,q)=C_n(q)(x \operatorname{ce}_n(x,q)+f_n(x,q)),$$ $$\operatorname{ge}_n(x,q)=S_n(q)(x \operatorname{se}_n(x,q)+g_n(x,q)),$$
where $f_n$ and $g_n$ are $2\pi$-periodic.
These functions can be thought of as
$$\operatorname{fe}_n(x,q)\propto\lim_{r\to n-} \frac{\operatorname{se}_r(x,q)}{r-n},$$
$$\operatorname{ge}_n(x,q)\propto\lim_{r\to n+} \frac{\operatorname{ce}_r(x,q)}{r-n}.$$
Thus, these additional functions $\operatorname{fe}$ and $\operatorname{ge}$ are nothing more than limit cases of correspondingly $\operatorname{se}$ and $\operatorname{ce}$.
I don't really see much of analogy with Bessel functions here since the second solutions are unbounded only in the case of integer $r$, and for any characteristic internal to the spectrum of Mathieu operator the equation has only bounded solutions. So there're actually no singular Neumann-like solutions in this case.