Let $\mathcal{D}(L)$ be as stated in the problem, meaning that $f \in \mathcal{D}(L)$ iff
- $f$ is continuously differentiable on $(-1,1)$,
- $f'$ is absolutely continuous on $(-1,1)$,
- $f$ and $Lf =-((1-x^{2})f')'$ are in $L^{2}(-1,1)$.
Aymptotics: The two classical solutions of $Lf=0$, both of which are in $L^{2}(-1,1)$:
$$
e_{1}(x)=\frac{1}{2}\ln\left(\frac{1-x}{1+x}\right),\;\;\; e_{2}(x)\equiv 1.
$$
If $f \in \mathcal{D}(L)$, then $Lf=0$ iff $f = Ae_{1}+Be_{2}$ for unique constants $A$, $B$. Let $(f,g)$ denote the inner-product of $f,g \in L^{2}(-1,1)$ and let $(f,g)_{a}^{b}$ denote the inner-product in $L^{2}(a,b)$ for $-1 \le a < b \le 1$. For $f, g \in \mathcal{D}(L)$ integration by parts gives
$$
(Lf,g)_{a}^{b}-(f,Lg)_{a}^{b} = (1-x^{2})(fg'-f'g)|_{a}^{b}.
$$
This seems innocent at first, but notice that the left side has limits as $a$, $b$ approach $-1$, $1$, which guarantees the existence of limits on the right side. For $g=e_{2}$, one has $e_{2}'=0$, which leads to
$$
(1-x^{2})f'(x) = f'(0)-(Lf,e_{2})_{0}^{x}
$$
So the following limit exists:
$$
A_{+}(f) = \lim_{x\rightarrow\pm 1}(1-x^{2})f'(x)
= f'(0)-\lim_{x\rightarrow 1}(Lf,e_{2})_{0}^{x}=f'(0)-(Lf,e_{2})_{0}^{1}.
$$
Furthermore, one has the asymptotic
$$
|(1-x^{2})f'(x)-A_{+}(f)|
= |(Lf,e_{2})_{x}^{1}| \le \|Lf\|_{x}^{1}\|e_{2}\|_{x}^{1}=o(\sqrt{1-x}).
$$
To obtain the second asymptotic, let $g=e_{1}$. Notice that $e_{1}'=-1/(1-x^{2})$. Therefore
$$
\begin{align}
(Lf,e_{2})_{0}^{x}
& =(1-x)^{2}(fe_{2}'-f'e_{2})|_{0}^{x} \\
& =-f(x)-(1-x^{2})f'(x)e_{2}(x)-f(0).
\end{align}
$$
Therefore,
$$
\begin{align}
f(x) & = -(Lf,e_{2})_{0}^{x}-(1-x^{2})f'(x)e_{2}(x)-f(0) \\
& = -(Lf,e_{2})_{0}^{x}+A_{+}(f)e_{2}(x)+o(\sqrt{1-x})e_{2}(x)-f(0).
\end{align}
$$
This gives
$$
f(x)-A_{+}(f)e_{2}(x)-B_{+}(f) = -(Lf,e_{2})_{0}^{x}+o(\sqrt{1-x})e_{2}(x)-f(0).
$$
So the stated limit $B_{\pm}(f)=\lim_{x\rightarrow +1}(f(x)-A_{+}e_{2}(x))$ exists and
$$
f(x)-A_{+}(f)e_{2}(x)-B_{+}(f) = o(\sqrt{1-x}e_{2})+O((Lf,e_{2})_{x}^{1})=o(\sqrt{1-x}e_{2}).
$$
Substituting back into the original identity,
$$
(Lf,g)-(f,Lg) = (1-x^{2})(fg'-f'g)|_{-1}^{1}.
$$
Near $x=1$, one has
$$
(1-x^{2})g'f-(1-x^{2})f'g \\
= (A_{+}(g)+o(\sqrt{1-x}))(B_{+}(f)+A_{+}(f)e_{2}+o(\sqrt{1-x}e_{2}))
- (A_{+}(f)+o(\sqrt{1-x}))(B_{+}(g)+A_{+}(g)e_{2}+o(\sqrt{1-x}e_{2})) \\
= A_{+}(g)B_{+}(f)-A_{+}(f)B_{+}(g)+o(\sqrt{1-x}e_{2}^{2}).
$$
Therefore,
$$
\begin{align}
(Lf,g)-(f,Lg) & =
-\left|\begin{array}{cc} A_{+}(f) & B_{+}(f) \\ A_{+}(g) & B_{+}(g)\end{array}\right|
+\left|\begin{array}{cc} A_{-}(f) & B_{-}(f) \\ B_{-}(g) & B_{-}(g)\end{array}\right|.
\end{align}
$$
Separated Conditions: We are told to assume that the unconstrained operator $L$ has adjoint $L^{\star} = L_{0}$, where $L_{0}$ is the restriction of $L$ to the functions $f \in \mathcal{D}(L)$ for which $A_{\pm}(f)=0$ and $B_{\pm}(f)=0$. It follows that $L_{0}$ is closed, and the graph of $L_{0}$ is of finite co-dimmension in the graph of $L$, which means that $L$ is also closed. So $L_{0}^{\star}=L$ also holds.
Let $\alpha,\beta \in [0,\pi)$ and define $L_{\alpha,\beta}$ as the restriction of $L$ to the $f \in \mathcal{D}(L)$ satisfying
$$
\cos\alpha A_{-}(f)+\sin\alpha B_{-}(f) = 0,\\
\cos\beta A_{+}(f) + \sin\beta B_{+}(f) = 0.
$$
We show that $L_{\alpha,\beta}$ is selfadjoint.
This operator $L_{\alpha,\beta}$ is symmetric on its domain because the two evaluation determinants are $0$ for the symmetric difference $(Lf,g)-(f,Lg)$, which is proved by noticing that the associated matrices have non-trivial homogenous solutions:
$$
\left[\begin{array}{cc} A_{-}(f) & B_{-}(f) \\ A_{-}(g) & B_{-}(g)\end{array}\right]
\left[\begin{array}{c} \cos\alpha \\ \sin\alpha\end{array}\right]=0,\;\;\;
\left[\begin{array}{cc} A_{+}(f) & B_{+}(f) \\ A_{+}(g) & B_{+}(g)\end{array}\right]
\left[\begin{array}{c} \cos\beta \\ \sin\beta\end{array}\right]=0.
$$
Therefore $L_{\alpha,\beta}\preceq L_{\alpha,\beta}^{\star}$.
Suppose $(L_{\alpha,\beta}f,g)=(f,h)$ for some $g,h\in L^{2}$ and all $f \in \mathcal{D}(L_{\alpha,\beta})$. Then this is true for all $f \in \mathcal{D}(L_{0})$ which implies that $h \in \mathcal{D}(L_{0}^{\star})=\mathcal{D}(L)$ and $Lg=h$. Because $(L_{\alpha,\beta}f,g)=(f,Lg)$ for all $f \in \mathcal{D}(L_{\alpha,\beta})$ then the left and right evaluation determinants must separately vanish. Define
$$
\Phi_{\alpha}(f) = +\cos\alpha A_{-}(f)+\sin\alpha B_{-}(f),\\
\Psi_{\alpha}(f) = -\sin\alpha A_{-}(f)+\cos\alpha B_{-}(f).
$$
Using muliplicative properties of determinant and the fact that $\Phi_{\alpha}(f)=0$ gives
$$
0= \left|\begin{array}{cc}
A_{-}(g) & B_{-}(g) \\
A_{-}(f) & B_{-}(f)
\end{array}\right|
= \left|\begin{array}{cc}
\Phi_{\alpha}(g) & \Psi_{\alpha}(g) \\
\Phi_{\alpha}(f) & \Psi_{\alpha}(f)
\end{array}\right| = \Phi_{\alpha}(g)\Psi_{\alpha}(f).
$$
This is true for all $f \in \mathcal{D}(L_{\alpha,\beta})$ which implies $\Phi_{\alpha}(g)=0$ for all $g \in \mathcal{D}(L_{\alpha,\beta}^{\star})$. The analogous result holds at the right endpoint as well. So $L_{\alpha,\beta}^{\star} \preceq L_{\alpha,\beta}$, which finishes the proof that $L_{\alpha,\beta}^{\star}=L_{\alpha,\beta}$.
Classical Conditions: The classical Legendre operator is $L_{0,0}$ which is equivalent to requiring $A_{\pm}(f)=0$ for $f \in \mathcal{D}(L)$. Or, because of the asymptotics, this is also equivalent to requiring the functions in the domain are bounded near $x=\pm 1$. Automatically such functions will have limits at $x=\pm 1$ because of the asymptotics derived for a general $f \in \mathcal{D}(L)$. So boundedness is equivalent:
$$
\mathcal{D}(L_{0,0})= \{ f \in \mathcal{D}(L) : f \mbox{ is bounded on } (-1,1) \}.
$$
This is why only eigenfunctions which are bounded near $\pm 1$ are considered valid. The only $f \in\mathcal{D}(L_{0,0})$ for which $Lf = \lambda f$ are the Legendre polynomials. The eigenvalues are $\lambda =n(n+1)$ for $n=0,1,2,3,\cdots$, and this set of eigenfunctions is a complete orthogonal basis of $L^{2}(-1,1)$. All kinds of other separted conditions are possible, but none of the others appears to be appropriate for applications to Physics and Engineering. The non-physical selfadjoint formulations do lead to orthogonal types of expansions.
Periodic Condition: I don't believe there is any reason that periodic types of conditions do not lead to selfadjoint restrictions. Conditions such as
$$
A_{+}(f) = A_{-}(f),\;\; B_{+}(f) = B_{-}(f)
$$
appear to be okay. However, $f(-1)=f(1)$, $f'(-1)=f'(1)$ does not make sense because $f \in \mathcal{D}(L)$ does not necessarily have limiting values for $f$ or for $f'$. And a single condition is not enough to determine a selfadjoint problem.
