Consider the associated Legendre's ODE given by $$ (x^2-1)y''+2xy'+\Bigg(\frac{1}{4}-m^2-\frac{1}{x^2-1}\Bigg)y=0 $$ the solutions of which are $$ y=C_1P_{m-1/2}^1(x)+C_2Q_{m-1/2}^1(x) $$ which are the associated Legendre functions of the first and second kinds of order 1 and degree $m-1/2$. The normalization of the functions of the first kind, $P_{m-1/2}^1$, is usually taken to be the same as the classical Legendre polynomials, $P_n(x)$, just expressed in terms of the Gamma function to handle the half-integer degree. Using this normalization the normalization of the function of the second kind, $Q_{m-1/2}^1$, may be obtained through a variation of parameters procedure. That is, we can define it by considering $Q_{m-1/2}^1(x)=A(x)P_{m-1/2}^1(x)$ and through variation of parameters on the above ODE, we'll obtain $$ Q_{m-1/2}^1(x)=P_{m-1/2}^1(x)\int\frac{dx}{P_{m-1/2}^1(x)^2(x^2-1)} + C $$
We could similarly do the procedure starting with $Q_{m-1/2}^1$ to define $P_{m-1/2}^1$. My question is, given the normalization for $P_{m-1/2}^1$, I am interested in the integral $$ \int\frac{dx}{Q_{m-1/2}^1(x)^2(x^2-1)}=\frac{1}{\frac{1}{4}-m^2}\frac{P_{m-1/2}^1(x)}{Q_{m-1/2}^1(x)}+C $$ I thought the pre-factor would be based on this original normalization chosen for $P_{m-1/2}^1$, but I can't seem to figure out where the $\frac{1}{4}-m^2$ comes from. The only reason I have it is because I have numerically validated it when calculating this integral for integer $m$. Does anyone know where this comes from? Any help is greatly appreciated!
Edit: (additional context and an update)
I have undertaken an in depth study of these functions as they have appeared in my research. Thus far, I have closely followed Hobson's Theory of Spherical and Ellipsoidal Harmonics among other references and have defined the function $Q_{m-1/2}^1$ as follows:
Start with the well known Legendre polynomials, $P_n(x)$, with classical normalization such that $P_n(1)=1$ for all $n$. Then, utilize variation of parameters to define the function of the second kind $Q_n(x)$ as was done above. Finally, using the derivative identity to obtain the function of order 1: $$ Q_n^1(x)=\sqrt{x^2-1}\frac{dQ_n}{dx} $$ and then replacing $n$ with $m-1/2$. From applying method of Frobenius to the original Legendre's ODE, a series expansion in negative powers of $x$ can be obtained. Comparing this with the result of the above variation of parameters procedure, the correct choice of the first term (and thus normalization) of the function of the second kind can be found that agrees with that of $P_n(x)$.
However, it seems that simply replacing $n$ with $m-1/2$ in expressions for $P_n(x)$ (leaving the original series nonterminating and then using the derivative identity) results in some problems with singularities I did not face while doing the same with $Q_n(x)$, so I'm wondering if my issue partly lies here.
Furthermore, upon further inspection, I have verified that the following is also true: $$ \int\frac{dx}{(x^2-1)P_{m-1/2}^1(x)^2} = \frac{1}{m^2-\frac{1}{4}}\frac{Q_{m-1/2}^1(x)}{P_{m-1/2}^1(x)} + C $$ which contradicts my original statement above! The strange factor of $\frac{1}{4}-m^2$ is appearing here as well, only negative. It seems the variation of parameters procedure I described above without this mysterious factor is correct for functions of order $0$ and degree $n$, but fails here by a constant factor. Thus far, all of my calculations agree with series and integral representations of these functions found here, and elsewhere.
My main question then becomes this: why do all of my independent calculations of representations of these functions agree with the literature, yet the simple application of variation of parameters to the associated ODE (not the original Legendre's ODE) results in an extra factor of $\pm(m^2-\frac{1}{4})$ when calculating these integrals numerically using the aforementioned representations of these functions? Again, thank you!
