There are a lot of great questions here, many of which I struggled with myself when I was first learning about the physical applications of PDEs. I'll go point by point and try to give the best explanations that I can. First, some housekeeping - whenever you use spherical coordinates, it is always good to explicitly state the coordinate transformation and notation you're using for reasons I mention in this answer and also discussed on Wolfram Mathworld. So, here we go:
I will use the $\displaystyle ( r,\theta ,\phi )$ (radial, azimuthal, polar) convention for spherical coordinates, as used by Wolfram Mathworld. Explicitly, the coordinate transformation is
\begin{equation*}
\begin{bmatrix}
x\\
y\\
z
\end{bmatrix} =\begin{bmatrix}
r\cos \theta \sin \phi \\
r\sin \theta \sin \phi \\
r\cos \phi
\end{bmatrix}
\end{equation*}
Yes, I know the theta and phi are in the opposite order as to what you're using, but this is just the notation that I'm the most comfortable with and trying to switch it around will only make me more prone to errors. Let's get on with actually answering your questions then.
First question: Why do we have complex solutions to what should be a real valued PDE?
This is probably the easiest question to answer. Take for instance the ordinary differential equation
$$y''+y=0$$
We can solve this in the complex plane using the usual methods as
$$y(x)=c_1e^{-ix}+c_2e^{ix}$$
This is the most general solution and allows for complex valued $y$, and, critically, complex valued coefficients. Going a bit further, we can break $c_1,c_2$ and $e^{-ix},e^{ix}$ into their real and imaginary parts:
$$y(x)=(a_1+ib_1)(\cos(x)-i\sin(x))+(a_2+ib_2)(\cos(x)+i\sin(x))$$
Which now allows us to split $y$ into its real and imaginary parts:
$$y(x)=\operatorname{Re}y+\operatorname{Im}y=\bigg[(a_1+a_2)\cos(x)+(b_1-b_2)\sin(x)\bigg]+i\bigg[(a_2-a_1)\sin(x)+(b_1+b_2)\cos(x)\bigg]$$
If we make the artificial restriction that $y$ must be real valued, i.e $\operatorname{Im}y=0$, we get $a_1=a_2~;~b_1=-b_2\implies c_1=c_2^*\implies y(x)=A\sin(x)+B\cos(x)$.
In short, if we were to write
$$y''+y=0\implies y(x)=c_1 e^{-ix}+c_2e^{ix}$$
This is still valid even if we only want real valued $y$! It's just implied that $c_1=c_2^*$. Why do we do this? Sometimes, writing the solution in terms of exponentials can make expressions a little bit shorter and easier to manipulate, even if it means we have to introduce complex valued coefficients.
Second question: Why does everyone just assume $l$ is an integer? And what do we do if $\alpha<-1/4$?
The second part of this question is in some ways easier to answer than the first. the answer to the first question is, sort of, "because that's the simplest possible case", but there's more to it than that (which I'll get to after I answer the second part.)
In short, complex numbers come to save the day again. Suppose $\alpha=-2$. Then,
$$l(l+1)=-2$$
$$l^2+l+2=0$$
$$l=\frac{-1\pm\sqrt{1-4\cdot2}}{2}=\frac{-1}{2}\pm i\frac{\sqrt{7}}{2}$$
This is totally allowed! While this will no longer be a Legendre polynomial, we do have what is known as the Legendre function (of the first kind):
$$\mathscr{P}_{l,m}(z)=\frac{1}{\Gamma(1-m)}\left(\frac{1+z}{1-z}\right)^{m/2}{}_2F_1\left([-l,l];[1-m];\frac{1-z}{2}\right)$$
Which uses the (Gaussian) hypergeometric function. It can be shown (though not easily) that in the case of $m,l\in \mathbb{N}$ and $m<l$ this reduces to the normal associated Legendre polynomial.
On to the first part then: why does everyone just assume $l$ is an integer? To answer this, let's solve the simpler problem of the spherical harmonics where we assume there is no $\theta$ (azimuthal) angle dependency. That is, let's have a look at the PDE
$$\nabla^2 U(r,\phi)=\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial U}{\partial r}\right)+\frac{1}{r^2\sin\phi}\frac{\partial}{\partial \phi}\left(\sin\phi \frac{\partial U}{\partial \phi}\right)$$
Assuming separability, i.e $U(r,\phi)=R(r)\Phi(\phi)$, we get (also note that I've multiplied through by $r^2$ which is fine as long as we don't care about the origin)
$$\frac{\partial}{\partial r}\left(r^2\frac{\partial R}{\partial r}\right)=\alpha R$$
And
$$\frac{1}{\sin\phi}\frac{\partial}{\partial \phi}\left(\sin\phi \frac{\partial \Phi}{\partial \phi}\right)=-\alpha \Phi$$
A bit of an inspired substitution $\alpha=l(l+1)$. This substitution is particularly nice if $\alpha$ is a positive integer of the form $n^2+n$, so that $l$ is a positive integer, but in general we just have $l=(-1\pm \sqrt{1-4\alpha})/2$. Recognizing that the $R$ equation is a Cauchy-Euler ODE, we get
$$R(r)=Ar^l+B\frac{1}{r^{l+1}}$$
Now for the $\Phi$ equation. We have, after some manipulation
$$\Phi''+\cot(\phi)\Phi'+l(l+1)\Phi=0$$
Which, after a substitution $\mu=\cos\phi$, becomes (I'll let you work out the details)
$$(1-\mu^2)\frac{\mathrm{d}^2\Phi}{\mathrm{d}\mu^2}-2\mu\frac{\mathrm{d}\Phi}{\mathrm{d}\mu}+l(l+1)\Phi=0$$
The reduced form of Legendre's differential equation.
So we can now resort to the Frobenius method and try to come up with a solution $\Phi(\mu)=\sum_{n=0}^\infty a_n \mu^n$. I'll cut to the chase (you can find the details on Wolfram Mathworld) and assert that after some work we get a recurrence
$$a_{n}=\frac{\left(l+(n+1)\right)(l-n)}{(n+1)(n+2)}a_{n-2}$$
So what's the upshot here? Yes, you are correct: $l$ doesn't have to be an integer for the above to simply be a solution to our differential equation. BUT, If $l$ wasn't an integer, then the equation $a_n=0$ has no solutions since $l+n+1=0$ or $l-n=0$ is clearly impossible if $l$ is not an integer. So, as long as $a_0\neq 0$, none of the coefficients in the expansion of $\Phi(\mu)$ would zero, i.e, it would be an infinite sum with quickly growing coefficients and thus would diverge for $\mu=1$ which corresponds to $\phi=0,\pi$. Since typically we are solving for an electric potential, it doesn't make any physical sense to have a electrical potential that blows up near the poles. However if $l$ is an integer, we will at some point reach some value of $n$ such that $a_n=0$ which would lead to $a_{n+2}=a_{n+4}=...=0$. All the previous even coefficients will be nonzero and their precise form will follow from the recurrence relation. As for the odd coefficients, i.e $a_{n+2m+1}$, we can simply choose all of these to be $0$. Since the recurrence is of second order, we need only worry about every second coefficient and we can choose the odd ones to all be zero which trivially satisfies the recurrence. Which half of the coefficients we set to $0$ depends on the parity of $l$. The precise series expansions for each case are listed on Wolfram Mathworld. To illustrate the above, the expansions will look like
$$P_l(\mu)=a_0\mu^0+0\mu^1 +a_2\mu^2+0\mu^3...+0\mu^{l-1}+a_l\mu^l+0\mu^{l+1}+0\mu^{l+2}+0+...\text{ for even }l$$
$$P_l(\mu)=0\mu^0+a_1\mu^1+0\mu^2+a_3\mu^3+...0\mu^{l-1}+a_l\mu^l+0\mu^{l+1}+0\mu^{l+2}+0+...\text{ for odd }l$$
As you can see, the expansions terminate which will mean that $\Phi$ will remain bounded at the poles.
To summarize:
No, $l$ need not be an integer in order to get something that solves the ODE. But, in order to get a solution of the ODE that remains bounded, $l$ must be an integer.
Finally, onto complete function sets and orthogonal bases. It's pretty obvious that $r^l$ form an basis (but not an orthogonal basis) for $L^2[0,\infty)$ - simply consider Taylor expansions. However, the Legendre polynomials do provide an orthogonal basis for $L^2[-1,1]$. We can arrive at this by applying the Gram-Schmidt process to $\{1,x,x^2,x^3,...\}$ on the interval $[-1,1]$. (REFERENCE: https://en.wikipedia.org/wiki/Orthogonal_functions#Polynomials) If we instead want an orthogonal basis on $[0,\infty)$, we can perform a Cayley Transform which brings us to the Legendre rational functions. Finally, about the completeness of $\{r^lY_{l,m}\}$ - I don't think this is required. $\{r^l Y_{l,m}\}$ need only form a basis, not necessarily an orthogonal one.
