Why are there four independent solutions of Mathieu equation instead of two? 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?
 A: You can define not ony two couples of independant solutions (C,S) and (ce,se), but as many couples as you want :
$$F=\alpha C + \beta S $$
$$G=\mu C + \nu S $$
where $\alpha, \beta, \mu, \nu$ are any constants, with condition $\alpha \nu - \beta \mu \neq 0$
Then $F$ and $G$ is a couple of independant solutions. But $C$, $S$, $F$, $G$ are not four independant solutions, just like $C$, $S$, $ce$, $se$.
A: 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.
