path: root/notes/FourierOnS2.tex

% vim:ts=2 sw=2 et:
\documentclass[a4paper, twocolumn]{article}

%
% Packages
%

%% styling for this document
\usepackage{tex/docstyle}
%\usepackage[firstpageonly, color={[gray]{0.9}}]{draftwatermark}

%% Theorems, TODO: move to docstyle
\usepackage{amsthm}

\newtheorem{theorem}{Theorem}
\newtheorem{lemma}{Lemma}
\newtheorem{definition}{Definition}
\newtheorem{remark}{Remark}
\newtheorem{proposition}{Proposition}

%
% Metadata
%

\title{Fourier on the Surface of the Sphere}
\author{Naoki Pross, Manuel Cattaneo}
\date{\today}

%
% Document
%

\begin{document}
\maketitle
\tableofcontents

% \section{Introduction}

%% FIXME
% Before we attack the main topic, it is good to refresh the normal Fourier theory on the simplest of surfaces, which is on a flat plane. Recall that Fourier's idea is to decompose a function \(f : \Omega \to \mathbb{C}\) into a set of spectral coefficients \(c\) which can be used to reconstruct \(f\). In the beginning we will let \(\Omega = \mathbb{R}^2\), since it is very similar to the 1D Fourier analysis.

\section{Fourier on a flat surface}

%% FIXME
First, we need to discuss about the class of functions on which the Fourier works at all. We will denote the set of such functions with \(C(\Omega; \mathbb{C})\), that means: continuous bounded functions from \(\Omega\) to the complex numbers. Informally, we will refer to these as ``nice'' functions. 

%% TODO: Check that R^2 / Z^2 = (R / Z) \times (R / Z)
To start we first pick \(\Omega = \mathbb{R}^2 / \mathbb{Z}^2\), i.e. the set of ``nice'' functions that are periodic with period 1. Then we need to define an inner product in this set of functions.

\begin{definition}[Inner product in \(C(\mathbb{R}^2 / \mathbb{Z}^2; \mathbb{C})\)]
  Let \(f(\mu, \nu), g(\mu, \nu) \in C(\mathbb{R}^2 / \mathbb{Z}^2; \mathbb{C})\). The inner product between \(f\) and \(g\) is
  \[
    \langle f, g \rangle =
    \iint_{[0, 1]^2} f g^* d\mu d\nu,
  \]
  where \(g^*\) denotes the complex conjugate of \(g\).
\end{definition}

With this construction, now we just need a suitable set of basis function for the decomposition. Again, recall that in the 1D Fourier analysis the basis functions are sines and cosines. Here it is very similar, we just have two dimensions instead of one. Therefore, we let
\[
  B_{m, n}(\mu, \nu) = \sin(\pi m \mu)\sin(\pi n \nu)
\]
where \(m, n \in \mathbb{Z}\), be our basis functions in the space of ``nice'' functions from \(\mathbb{R}^2\) to \(\mathbb{C}\). Like in the one dimensional Fourier analysis, we can now define the Fourier coefficients.

\begin{definition}[Fourier coefficients]
  Let \(f(\mu, \nu) \in C(\mathbb{R}^2/\mathbb{Z}^2; \mathbb{C})\). The numbers
  \begin{align*}
    c_{m, n} &= \langle f, B_{m, n} \rangle \\
             &= \iint_{[0, 1]^2} f(\mu, \nu) \sin(\pi m \mu)\sin(\pi n \nu) d\mu d\nu,
  \end{align*}
  are called the Fourier coefficients of \(f\).
\end{definition}

And finally by the Fourier theorem we can reconstruct the original function using a Fourier series:
\begin{equation}\label{eqn:fourier-theorem}
  f(\mu, \nu) =
    \sum_{m\in\mathbb{Z}} \sum_{n\in\mathbb{Z}} c_{m, n} B_{m, n}(\mu, \nu).
\end{equation}

% TODO: keep or remove? <=> discuss fourier theorem?
% \begin{definition}[\(L^2\) norm]
%   \[
%     \|f\|_2  = \sqrt{\langle f, f \rangle}
%   \]
% \end{definition}
% 
% \begin{theorem}[Fourier Theorem]
%   Let \(f : \mathbb{R}^2 \to \mathbb{C}\) be a ``nice'' function. 
%   \[
%     \lim_{N \to \infty} \left\|
%       f - \sum_{n=-N}^N c_n e^{in}
%     \right\|_2 = 0
%   \]
% \end{theorem}

\subsection{Why sines?}

A important question now is: Why did we choose \(B_{m, n}\) to be a multiplication of two sines?
The answer has to do with solving other problems. That is because originally Fourier developed its theory to solve some difficult problems in thermodynamics. 

The problem we want to start with is to find the eigenfunctions of the Laplace operator on the two-dimensional square domain $\Omega \in [0,1]^2$ with zero Dirichlet boundary conditions, namely
\begin{align} 
  \nabla^2 f &= \lambda f \label{eqn:laplace-cartesian},\\
           f(\mu, \nu) &= 0, \quad \mu, \nu \in \partial \Omega. \label{eqn:laplace-cartesian-bound}
\end{align}
What we need are orthogonal functions, so that we can expand a series as the one shown in Eq.\eqref{eqn:fourier-theorem}. If the notion of ``orthogonality'' reminds you the world of linear algebra you are on the right way. 

Eq.\eqref{eqn:laplace-cartesian} is indeed very reminiscent of the problem of finding the eigenvectors of a matrix. In fact, if you replace the Laplace operator $\nabla^2$ with a matrix $\mathbf{A}$ and the function $f$ with a vector $\vec{v}$, you get the well-known matrix equation
\begin{equation*}
\mathbf{A} \vec{v} = \lambda \vec{v},
\end{equation*}
which, when solved, leads to the eigenvalues and associated eigenvectors.

This intuition makes perfect sense, since as we know from the \emph{spectral decomposition theorem}, we can reconstruct any (diagonizable) matrix $\mathbf{A}$ by only knowing its eigenvalues and their respective eigenvectors. This can be performed according to the following equation:
\begin{equation*}
\mathbf{A} \vec{v} = \mathbf{Q} \mathbf{\Lambda} \mathbf{Q}^{-1},
\end{equation*}
where the columns of $\mathbf{Q}$ are the eigenvectors and the matrix $\mathbf{\Lambda}$ consists only of zeros except for the diagonal, in which the eigenvalues are contained.

So if we reconsider our problem, once we find a set of orthogonal functions, they can be used to reconstruct any ``nice'' function, supported on the domain $\Omega$. This is analogous to the previous case, where it was the matrix $\mathbf{A}$ that could be reconstructed using a set of orthogonal vectors.

In order to get the mentioned set, we have to solve the eigenfunction equation \eqref{eqn:laplace-cartesian} for the laplacian operator. This corresponds to solve the following PDE:
\[
 \frac{\partial^2 f}{\partial \mu^2}
                       + \frac{\partial^2 f}{\partial \nu^2} = \lambda f.
\]
The following ansatz can be performed
\[
  f(\mu, \nu) = M(\mu)N(\nu),
\]
which when substituted into \eqref{eqn:laplace-cartesian} yields to
\[
  \frac{d^2 M}{d\mu^2} N(\nu)
  + \frac{d^2 N}{d\nu^2} M(\mu) = \lambda M(\mu)N(\nu).
\]
Notice that the partial derivatives have been simplified to normal derivatives. Continuing the separation method, we divide by $M(\mu)N(\nu)$ \texttt{TODO: ASSUMPTION $f \neq 0$ (?)}, obtaining:
\[
  \frac{d^2 M}{d\mu^2} \frac{1}{M(\mu)}
  + \frac{d^2 N}{d\nu^2} \frac{1}{N(\nu)} = \lambda.
\]
We now let
\begin{align}
\alpha^2 &:= \frac{d^2 N}{d\nu^2} \frac{1}{N(\nu)} -\lambda \quad \text{and} \label{eqn:def1-simple-ode}\\
\beta^2 &:= \frac{d^2 M}{d\mu^2} \frac{1}{M(\mu)} -\lambda. \label{eqn:def2-simple-ode}
\end{align}
We see that this results in two almost identical problems of the form
\begin{align*}
\frac{d^2 M}{d\mu^2} = -\alpha^2 M(\mu), \\
\frac{d^2 N}{d\nu^2} = -\beta^2 N(\nu) .
\end{align*}
% TODO: prove / discuss that w (=m or n) must be integers because of the periodicity of f
The solutions to these elementary ODEa are a linear combination of sines and cosines, which can be written as
\begin{align*}
M(\mu) &= A_1 \sin(\alpha \mu) + A_2 \cos(\alpha \mu) \quad \text{and} \\
N(\nu) &= B_1 \sin(\beta \nu) + B_2 \cos(\beta \nu).
\end{align*}
The unknown parameters of $M(\mu)$ can be estimated using the Dirichlet boundary conditions as follows:
\begin{align*}
M(0) \stackrel{!}{=} 0 \quad &\Rightarrow  A_2 = 0, \\
N(1) \stackrel{!}{=} 0 \quad &\Rightarrow  \alpha = n\pi, \quad n \in \mathbb{N}.
\end{align*}
The same method can be used to compute the coefficients for $N(\nu)$, which are
\begin{align*}
B_2 &= 0, \\
\beta &= m\pi, \quad m \in \mathbb{N}.
\end{align*}
This means that we have just found an infinite number of eigenfunctions in the form
\begin{align*}
f_{m,n}(\mu,\nu)&= M_m(\mu)N_n(\nu) \\
                &= c_{m,n} \sin(\pi m \mu)\sin(\pi n \nu).
\end{align*}
The whole solution can be built up by a linear combination of these functions
\begin{equation*}
f(\mu,\nu) = \sum_{m=0}^\infty \sum_{n=0}^\infty c_{m,n} \sin(\pi m \mu)\sin(\pi n \nu)
\end{equation*}
of which, however, the multiplication coefficient $c_{m,n}$ remains to be found. 

This problem will not be adressed because it is beyond the scope of this chapter. An intuition to solve it is to use other conditions for the function $f(\mu,\nu)$, in order to turn it into a ``Fourier-style'' problem, i.e. to find the coefficients of some basis functions $B_{m,n}(\mu,\nu)$ (in this case $B_{m,n}(\mu,\nu)=\sin(\pi m \mu)\sin(\pi n \nu)$), in order to reconstruct $f(\mu,\nu)$.

We can observe that the eigenfunctions which solve Eq.\eqref{eqn:laplace-cartesian} are the same we used to build the Fourier theory. This is not a coincidence, in fact quite the opposite: the basis functions of the Fourier decomposition were chosen such that the Laplacian operator is easy in the frequency domain. In other words, such that the expression
\[
  \langle \nabla^2 f, B_{m, n} \rangle
\]
is easy to compute. This is shown in the next lemma.

\begin{lemma}
  Let the eigenfunction $f(\mu,\nu) \in C(\mathbb{R}^2/\mathbb{Z}^2; \mathbb{C})$, then
  \begin{equation}\label{eqn:lemma1-eigenfunction}
    \langle \nabla^2 f, B_{m,n} \rangle = -\pi^2 \left( m^2 + n^2 \right) \langle f, B_{m,n} \rangle
  \end{equation}
\end{lemma}
\begin{proof}
  We know from Eq.\eqref{eqn:laplace-cartesian} that
  \begin{equation*}
    \nabla^2 f = \lambda f.
  \end{equation*}
  This implies
  \begin{equation*}
    \langle \nabla^2 f, B_{m,n} \rangle = \langle \lambda f, B_{m,n} \rangle.
  \end{equation*}
  Given the linearity of the inner product we can now write
  \begin{align*}
     \langle \lambda f, B_{m,n} \rangle &= \iint_{[0, 1]^2} \lambda f B_{m,n}^* d\mu d\nu \\
                          &= \lambda \iint_{[0, 1]^2} f B_{m,n}^* d\mu d\nu \\
                          &= \lambda \langle f, B_{m,n} \rangle.
  \end{align*}
  Proving Eq.\eqref{eqn:lemma1-eigenfunction} is therefore equivalent to proving the following relationship: 
  \begin{equation*}
  \lambda = -\pi^2 \left( m^2 + n^2 \right).
  \end{equation*}
  From the two definitions in Eqs.\eqref{eqn:def1-simple-ode}\eqref{eqn:def2-simple-ode}, one can write
  \begin{align*}
  \lambda &= -\alpha^2 -\beta^2 = -\left(\alpha^2 + \beta^2\right) \\
  &= -\pi^2 \left(m^2 + n^2\right)
  \end{align*}
  
\end{proof}

% TODO: closing words: this is why fourier is useful

\section{Fourier on the Sphere}

\subsection{The hard problem}

Like in the previous case, the motivation for the construction of a Fourier theory is a hard problem involving derivatives. In this case, we want to solve problems that are spherically symmetric, something that is found very often in Physics (potential around a point charge, atomic orbitals, gravitational fields of planets, etc.).

In this case the equation for which solutions are sought is
\begin{equation} \label{eqn:surface-harmonics}
  \nabla_2^2 f(\vartheta, \varphi) = 0,
\end{equation}
where \(f\) is a function on the unit sphere and \(\nabla_2^2\) is the \emph{surface Laplacian}, which is defined to be:
\[
  \nabla_2^2 =
    \frac{1}{\sin\vartheta} \frac{\partial}{\partial \vartheta} \left(
      \sin \vartheta \frac{\partial}{\partial \vartheta}
    \right) + \frac{1}{\sin^2 \vartheta} \frac{\partial^2}{\partial \varphi^2}.
\]
The subscript is there to hint that this is a derivative on the unit sphere \(S^2\). The surface Laplacian can also be defined in term of the normal Laplacian in spherical coordinates, by removing the radial component:
\[
  \nabla_2^2 = r \nabla^2 - r \frac{\partial^2}{\partial r^2} r.
\]

Like in the flat case \eqref{eqn:surface-harmonics} is solved with a product ansatz
\[
  f(\vartheta, \varphi) = \Theta(\vartheta) \Phi(\varphi).
\]
Though, unfortunately this time the separation process is more involved, and the results more complicated. The separation with the separation variable \(m\) yields the following ODEs:
\begin{subequations}
  \begin{align}
    \label{eqn:separation-phi}
    0 &= \frac{d^2\Phi}{d\varphi^2} \frac{1}{\Phi(\varphi)} \\
    \label{eqn:separation-theta}
    0 &= \frac{1}{\sin\vartheta} \frac{d}{d\vartheta} \left(
      \sin\vartheta \frac{d\Theta}{d\vartheta}
    \right) \nonumber\\ &\qquad + \left[
      n(n+1) - \frac{m}{\sin^2\theta}
    \right] \Theta(\vartheta)
  \end{align}
\end{subequations}

Equation \eqref{eqn:separation-phi} is easy, the solutions are complex exponentials \(e^{im\varphi}\), while \eqref{eqn:separation-theta} is known as the \emph{associated Legendre equation}. Though, normally the equation is written in term of \(x\) and \(y(x)\), so \eqref{eqn:separation-theta} is brought to a more familiar form by using the substitution \(x = \cos\vartheta\) and \(y = \Theta\):
\[
  \left( 1 - x^2 \right) \frac{d^2 y}{dx^2} - 2x \frac{dy}{dx}
    + \left[ n(n+1) - \frac{m^2}{1 - x^2} \right] y(x) = 0.
\]
Finding the solutions to this equation is so involved, that it deserves its own section.

\subsection{The associated Legendre polynomials}

In this section we would like to find the solutions to the \emph{associated} Legendre equation, which is actually a generalization of Legendre's equation:
\begin{equation} \label{eqn:legendre}
  \left( 1 - x^2 \right) \frac{d^2 y}{dx^2}
    - 2x \frac{dy}{dx} + n(n + 1) y(x) = 0.
\end{equation}
Thus we first need examine the solutions to this equation before constructing the more general solution.

\begin{proposition}
  The polynomials
  \begin{equation} \label{eqn:legendre-poly}
    P_n(x) = \sum_{k=0}^{\lfloor n/2 \rfloor} 
    \frac{(-1)^k (2n-2k)!}{2^n k! (n-k)!(n-2k)!} x^{n-2k},
  \end{equation}
  are solutions to Legendre's equation \eqref{eqn:legendre} when \(n > 0\).
\end{proposition}
\begin{proof} See appendix. \end{proof}

The proof for this proposition is quite algebraically involved and is thus left in the appendix. Since this is a power series \eqref{eqn:legendre-poly} can also be rewritten using Gauss' Hypergeometric function.

\begin{proposition}
  The polynomial \eqref{eqn:legendre-poly} can we rewritten using Gauss' Hypergeometric function
  \[
    {}_2F_1 \left( \begin{matrix}
      a_1, & a_2 \\ \multicolumn{2}{c}{b}
    \end{matrix} ; \frac{1 - x}{2} \right)
    =
    \sum_{k = 0}^\infty \frac{(a_1)_k (a_2)_k}{(b)_k} \frac{x^k}{k!},
  \]
  where the notation \((a)_k\) is for the Pochhammer Symbol
  \[
    (a)_k = a (a + 1) \ldots (a + k - 1).
  \]
  Hence for \(x \in (-1, 1)\) and \(n \in \mathbb{R}\):
  \[
    P_n (x) = {}_2F_1 \left( \begin{matrix}
      n + 1, & -n \\ \multicolumn{2}{c}{1}
    \end{matrix} ; \frac{1 - x}{2} \right).
  \]
\end{proposition}

In some applications, such as in quantum mechanics, it is more common to see it written yet in another form using Rodrigues' Formula.

\begin{proposition} The expression
  \begin{equation} \label{eqn:legendre-rodrigues}
    P_n(x) = \frac{1}{n!2^n}\frac{d^n}{dx^n}(x^2-1)^n.
  \end{equation}
  is equivalent to \eqref{eqn:legendre-poly}.
\end{proposition}

\begin{proof}
  We start expanding the term \((x^2-1)^n\); According to the binomial theorem
  \begin{equation*}
    (x^2-1)^n=\sum_{k=0}^n(-1)^k\binom{n}{k}x^{2(n-k)}.
  \end{equation*}
  Substituting the above, \eqref{eqn:legendre-rodrigues} becomes
  \begin{align*}
    \frac{1}{n!2^n} \frac{d^n}{dx^n} (x^2 - 1)^n 
      &= \frac{1}{n!2^n} \sum_{k=0}^n(-1)^k \binom{n}{k} \frac{d^n}{dx^n}x^{2(n-k)} \\
      &= \frac{1}{n!2^n} \sum_{k=0}^{\lfloor n / 2 \rfloor}(-1)^k
         \binom{n}{k}\frac{d^n}{dx^n}x^{2(n-k)}.
  \end{align*}
  Recall that
  \[
    \frac{d^n}{dx^n}x^\alpha
      % = \alpha(\alpha-1)(\alpha-2)\hdots(\alpha-n+1)x^{\alpha-n}
      = \frac{\alpha!}{(\alpha-n)!}x^{\alpha-n},
  \]
  thus
  \begin{align*}
    \frac{1}{n!2^n} &\sum_{k=0}^{\lfloor n / 2 \rfloor}(-1)^k
       \binom{n}{k}\frac{d^n}{dx^n}x^{2(n-k)} \\
    &= \frac{1}{n!2^n} \sum_{k=0}^{\lfloor n/2 \rfloor} (-1)^k \binom{n}{k}
        \frac{(2n-2k)!}{(n-2k)!}x^{n-2k} \\
    &= \frac{1}{n!2^n} \sum_{k=0}^{\lfloor n/2 \rfloor} (-1)^k \frac{n!}{k!(n-k)!}
        \frac{(2n-2k)!}{(n-2k)!}x^{n-2k} \\
    &= \sum_{k=0}^{\lfloor n/2 \rfloor} \frac{(-1)^k (2n-2k)!}{2^n k! (n-k)!(n-2k)!} x^{n-2k}.
    \qedhere
  \end{align*}
\end{proof}

Now, using the solutions to the Legendre equation we can construct the solution to the more general problem:
\begin{align}
  \left( 1 - x^2 \right) \frac{d^2 y}{dx^2} &- 2x \frac{dy}{dx} \nonumber \\
    & + \left[ n(n+1) - \frac{m^2}{1 - x^2} \right] y(x) = 0.
    \label{eqn:legendre-as}
\end{align}

This equation is considerably more difficult, and again, we will just analyze the solution.

\begin{definition}[Associated Legendre Polynomials]
  Let \(m \in \mathbb{N}_0\). The polynomials
  \begin{equation} \label{eqn:legendre-poly-as}
    P_{m, n} (x) = \left( 1 - x^2 \right)^{m/2} \frac{d^m}{dx^m} P_n (x),
  \end{equation}
  are called the associated Legendre polynomials.
\end{definition}

\begin{lemma}
  The associated Legendre polynomials \eqref{eqn:legendre-poly-as} are solutions to the associated Legendre differential equation \eqref{eqn:legendre-as}.
\end{lemma}
\begin{proof} See appendix. \end{proof}

\subsection{Spherical harmonics}

\clearpage
\appendix
\section{Proofs}
\subsection{Legendre Polynomials}

\begin{lemma}
  The polynomial
  \[
    P_n(x) = \sum_{k=0}^{\lfloor n/2 \rfloor} (-1)^k
      \frac{(2n-2k)!}{2^n k! (n-k)!(n-2k)!} x^{n-2k},
  \]
  is a solution to Legendre's equation
  \[
    \left( 1 - x^2 \right) \frac{d^2 y}{dx^2}
      - 2x \frac{dy}{dx} + n(n + 1) y(x) = 0
  \]
  for \(n > 0\).
\end{lemma}
\begin{proof}
To solve \eqref{eqn:legendre} we use the power series ansatz
\begin{equation} \label{eqn:legendre-power-ansatz}
  y(x) = \sum_{k=0}^\infty a_k x^k,
\end{equation}
from which follows that
\[
  y' = \sum_{k = 0}^\infty k a_k x^{k-1}, \text{ and }
  y'' = \sum_{k = 0}^\infty k (k-1) a_k x^{k-2}.
\]
By substituting the above and \eqref{eqn:legendre-power-ansatz} into \eqref{eqn:legendre} we get the that first term
\begin{align*}
  \Big( 1 &- x^2 \Big) y''
    = \left( 1 - x^2 \right) \sum_{k = 0}^\infty k (k-1) a_k x^{k - 2} \\
    &= \sum_{k = 0}^\infty k (k-1) a_k x^{k-2} + k (k-1) a_k x^{k} \\
    &= \sum_{k = 0}^\infty \left[
         (k + 1) (k + 2) a_{k + 2} + k (k-1) a_k 
       \right] x^{k},
\end{align*}
where in the last step to factor out \(x^k\) we shifted the index in the coefficients by 2, i.e.
\[
  \sum_{k = 0}^\infty k (k - 1)a_k x^{k - 2}
  = \sum_{k = 0}^\infty (k + 2) (k + 1) a_{k + 2} x^k.
\]
Similarly, the second term:
\[
  -2xy' = -2x \sum_{k = 0}^\infty k a_k x^{k - 1}
        = \sum_{k = 0}^\infty - 2k a_k x^k.
\]
Finally, combining the above the complete substitution yields
\begin{gather*}
  \Big( 1 - x^2 \Big) y'' - 2xy' + n(n + 1) y  = 0 \\
    \implies \sum_{k = 0}^\infty \big[
         (k + 1) (k + 2) a_{k + 2} + k (k-1) a_k \\
         \qquad \qquad
         - 2k a_k + n(n + 1) a_k
       \big] x^k = 0,
\end{gather*}
from which we can extract the recurrence relation
\begin{gather*}
  \begin{aligned}
    (k + 1) (k + 2) a_{k + 2} &+ k (k-1) a_k \\
    &- 2k a_k + n(n + 1) a_k = 0
  \end{aligned} \\
  \iff a_{k + 2} = \frac{(k-n)(k+n+1)}{(k+2)(k+1)} a_k.
\end{gather*}

% TODO: finish reviewing proof
  \texttt{[TODO: finish copying proof]}
\if 0

We can derive a recursion formula for $a_{k+2}$ from Eq.\eqref{eq:condition_2}, which can be expressed as
\begin{equation}\label{eq:recursion}
a_{k+2}= \frac{k (k-1) - 2 k + n(n+1)}{(k+1)(k+2)}a_k = \frac{(k-n)(k+n+1)}{(k+2)(k+1)}a_k.
\end{equation}
All coefficients can be calculated using the latter. 

Following Eq.\eqref{eq:recursion}, if we want to compute $a_6$ we would have
\begin{align*}
a_{6}= -\frac{(n-4)(n+5)}{6\cdot 5}a_4 &= -\frac{(n-4)(5+n)}{6 \cdot 5} -\frac{(n-2)(n+3)}{4 \cdot 3}  a_2 \\
&= -\frac{(n-4)(n+5)}{6 \cdot 5} -\frac{(n-2)(n+3)}{4 \cdot 3} -\frac{n(n+1)}{2 \cdot 1}  a_0 \\
&= -\frac{(n+5)(n+3)(n+1)n(n-2)(n-4)}{6!} a_0.
\end{align*}
One can generalize this relation for the $i^\text{th}$ even coefficient as
\begin{equation*}
a_{2k} = (-1)^k \frac{(n+(2k-1))(n+(2k-1)-2)\hdots (n-(2k-2)+2)(n-(2k-2))}{(2k)!}a_0
\end{equation*}
where $i=2k$.

A similar expression can be written for the odd coefficients $a_{2k-1}$. In this case, the equation starts from $a_1$ and to find the pattern we can write the recursion for an odd coefficient, $a_7$ for example
\begin{align*}
a_{7}= -\frac{(n-5)(n+6)}{7\cdot 6}a_5 &= - \frac{(n-5)(n+6)}{7\cdot 6} -\frac{(n-3)(n+4)}{5 \cdot 4}  a_3 \\
&= - \frac{(n-5)(n+6)}{7\cdot 6} -\frac{(n-3)(n+4)}{5 \cdot 4}  -\frac{(n-1)(n+2)}{3 \cdot 2}  a_1 \\
&= -\frac{(n+6)(n+4)(n+2)(n-1)(n-3)(n-5)}{7!} a_1.
\end{align*}
As before, we can generalize this equation for the $i^\text{th}$ odd coefficient
\begin{equation*}
a_{2k+1} = (-1)^k \frac{(n + 2k)(n+2k-2)\hdots(n-(2k-1)+2)(n-(2k-1))}{(2k+1)!}a_1
\end{equation*}
where $i=2k+1$.

Let be 
\begin{align*}
y_\text{e}^K(x) &:= \sum_{k=0}^K(-1)^k \frac{(n+(2k-1))(n+(2k-1)-2)\hdots \color{red}(n-(2k-2)+2)(n-(2k-2))}{(2k)!} x^{2k}, \\
y_\text{o}^K(x) &:= \sum_{k=0}^K(-1)^k \frac{(n + 2k)(n+2k-2)\hdots \color{blue} (n-(2k-1)+2)(n-(2k-1))}{(2k+1)!} x^{2k+1}.
\end{align*}
The solution to the Eq.\eqref{eq:legendre} can be written as
\begin{equation}\label{eq:solution}
y(x) = \lim_{K \to \infty} \left[ a_0 y_\text{e}^K(x) + a_1 y_\text{o}^K(x) \right].
\end{equation}

The colored parts can be analyzed separately:
\begin{itemize}
  \item[\textcolor{red}{\textbullet}] Suppose that $n=n_0$ is an even number. Then the red part, for a specific value of $k=k_0$, will follow the following relation:
\begin{equation*}
n_0-(2k_0-2)=0. 
\end{equation*}
From that point on, given the recursive nature of Eq.\eqref{eq:recursion}, all the subsequent coefficients will also be 0, making the sum finite.
\begin{equation*}
a_{2k}=0 \iff y_{\text{o}}^{2k}(x)=y_{\text{o}}^{2k_0}(x), \quad \forall k>k_0
\end{equation*} 
  \item[\textcolor{blue}{\textbullet}] Suppose that $n=n_0$ is an odd number. Then the blue part, for a specific value of $k=k_0$, will follow the following relation 
\begin{equation*}
n_0-(2k_0-1)=0.  
\end{equation*}
From that point on, for the same reason as before, all the subsequent coefficients will also be 0, making the sum finite.
\begin{equation*}
a_{2k+1}=0 \iff y_{\text{o}}^{2k+1}(x)=y_{\text{o}}^{2k_0+1}(x), \quad \forall k>k_0
\end{equation*}
\end{itemize} 

There is the possibility of expressing the solution in Eq.\eqref{eq:solution} in a more compact form, combining the two solutions $y_\text{o}^K(x)$ and $y_\text{e}^K(x)$. They are both a polynomial of maximum degree $n$, assuming $n \in \mathbb{N}$. In the case where $n$ is even, the polynomial solution
\begin{equation*}
\lim_{K\to \infty} y_\text{e}^K(x)
\end{equation*}
will be a finite sum. If instead $n$ is odd, will be 
\begin{equation*}
\lim_{K\to \infty} y_\text{o}^K(x)
\end{equation*}
to be a finite sum. 

Depending on the coefficient we start with, $a_1$ or $a_0$, we will obtain the odd or even polynomial respectively. Starting with the last coefficient $a_n$ and, recursively, calculating all the others in descending order, we can express the two parts $y_\text{o}^K(x)$ and $y_\text{e}^K(x)$ with a single sum. Hence, because we start with the last coefficient, the choice concerning $a_1$ and $a_0$ will be at the end of the sum, and not at the beginning. To compact Eq.\eqref{eq:solution}, Eq.\eqref{eq:recursion} can be reconsidered to calculate the coefficient $a_{k-2}$, using $a_k$
\begin{equation*}
a_{k-2} = -\frac{(k+2)(k+1)}{(k-n)(k+n+1)}a_k
\end{equation*}
Now the game is to find a pattern, as before. Remember that $n$ is a fixed parameter of Eq.\eqref{eq:legendre}.  
\begin{align*}
a_{n-2} &= -\frac{n(n-1)}{2(2n-1)}a_n, \\
a_{n-4} &= -\frac{(n-2)(n-3)}{4(2n-3)}a_{n-2} \\
&= -\frac{(n-2)(n-3)}{4(2n-3)}-\frac{n(n-1)}{2(2n-1)}a_n.
\end{align*}
In general 
\begin{equation}\label{eq:general_recursion}
a_{n-2k} = (-1)^k \frac{n(n-1)(n-2)(n-3) \hdots (n-2k+1)}{2\cdot4\hdots 2k(2n-1)(2n-3)\hdots(2n-2k+1)}a_n
\end{equation}
The whole solution can now be written as
\begin{align}
y(x) &= a_n x^n + a_{n-2} x^{n-2} + a_{n-4} x^{n-4} + a_{n-6} x^{n-6} + \hdots + \begin{cases} 
a_1 x, \quad &\text{if } n \text{ odd} \\ 
a_0, \quad  &\text{if } n \text{ even} 
\end{cases} \nonumber \\
&= \sum_{k=0}^{\lfloor \frac{n}{2} \rfloor} a_{n-2k}x^{n-2k} \label{eq:solution_2}
\end{align}
By considering
\begin{align}
(2n-1)(2n-3)\hdots (2n-2k+1)&=\frac{2n(2n-1)(2n-2)(2n-3)\hdots(2n-2k+1)}
{2n(2n-2)(2n-4)(2n-6)\hdots(2n-2k+2)} \nonumber \\ 
&=\frac{\frac{(2n)!}{(2n-2k)!}}
{2^kn(n-1)(n-2)(n-3)\hdots(n-k+1)} \nonumber \\
&=\frac{\frac{(2n)!}{(2n-2k)!}}
{2^k\frac{n!}{(n-k)!}}=\frac{(n-k)!(2n)!}{n!(2n-2k)!2^k} \label{eq:1_sub_recursion}, \\
2 \cdot 4  \hdots 2k &= 2^r 1\cdot2 \hdots r = 2^r r!\label{eq:2_sub_recursion}, \\
n(n-1)(n-2)(n-3) \hdots (n-2k+1) &= \frac{n!}{(n-2k)!}\label{eq:3_sub_recursion}.
\end{align}
Eq.\eqref{eq:solution_2} can be rewritten as
\begin{equation}\label{eq:solution_3}
y(x)=a_n \sum_{k=0}^{\lfloor \frac{n}{2} \rfloor} (-1)^k \frac{n!^2(2n-2k)!}{k!(n-2k)!(n-k)!(2n)!}  x^{n-2k}.
\end{equation}
Eq.\eqref{eq:solution_3} is defined for any $a_n$. By letting $a_n$ be declared as
\begin{equation*}
a_{n} := \frac{(2n)!}{2^n n!^2},
\end{equation*}
the so called \emph{Legendre polynomial} emerges
\begin{equation}\label{eq:leg_poly}
P_n(x):=\sum_{k=0}^{\lfloor \frac{n}{2} \rfloor} (-1)^k \frac{(2n-2k)!}{2^n k! (n-k)!(n-2k)!} x^{n-2k} 
\end{equation}
\fi

\end{proof}


\end{document}