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

%
% Packages
%

%% styling for this document
\usepackage{tex/docstyle}

%% 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 complex exponentials. Here it is no different, we just have two dimensions instead of one. Therefore, we let
\[
  B_{m, n}(\mu, \nu) = e^{i2\pi m\mu} e^{i2\pi n\nu} = e^{i2\pi (m\mu + 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
  \[
    c_{m, n} = \langle f, B_{m, n} \rangle
             = \iint_{[0, 1]^2} f(\mu, \nu) e^{-i2\pi (m\mu + n\nu)} d\mu d\nu,
  \]
  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:
\[
  f(\mu, \nu) =
    \sum_{m\in\mathbb{Z}} \sum_{n\in\mathbb{Z}} c_{m, n} B_{m, n}(\mu, \nu).
\]

% 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 complex exponentials?}

A important question now is: Why did we choose \(B_{m, n}\) to be complex exponentials?
The answer has to do with solving other problems. That is because originally Fourier developed its theory to solve some difficult problems in thermodynamics, where he wanted to solve (among many other equations)
\begin{equation} \label{eqn:laplace-cartesian}
  \nabla^2 f(\mu, \nu) = \frac{\partial^2 f}{\partial \mu^2}
                       + \frac{\partial^2 f}{\partial \nu^2} = 0.
\end{equation}

This PDE is known as Laplace's equation, and can be solved by separation using the ansatz
\[
  f(\mu, \nu) = M(\mu)N(\nu),
\]
which when substituted into \eqref{eqn:laplace-cartesian} yields
\[
  \frac{d^2 M}{d\mu^2} N(\nu)
  + \frac{d^2 N}{d\nu^2} M(\mu) = 0.
\]
Notice that the partial derivatives have been simplified to normal derivatives. Continuing the separation method, we divide by \(M(\mu)N(\nu)\), obtaining:
\[
  \underbrace{\frac{d^2 M}{d\mu^2} \frac{1}{M(\mu)}}_{w}
  + \underbrace{\frac{d^2 N}{d\nu^2} \frac{1}{N(\nu)}}_{-w} = 0.
\]
We let \(w\) be the separation constant, and we see that this results in two almost identical problems of the form
\[
  \frac{d^2 X}{d\xi^2} = \pm w X(\xi).
\]
% TODO: prove / discuss that w (=m or n) must be integers because of the periodicity of f
The solutions to this elementary ODE are of course complex exponentials, 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 \(f \in C(\mathbb{R}^2/\mathbb{Z}^2; \mathbb{C})\), then
  \[
    \langle \nabla^2 f, B_{m, n} \rangle
    = (2\pi i)^2 \left( m^2 + n^2 \right) \langle f, B_{m, n} \rangle.
  \]
\end{lemma}
\begin{proof}
  To start, we first expand the left side of the statement:
  \begin{gather} 
    \nonumber
    \langle \nabla^2 f, B_{m, n} \rangle
    = \iint_{[0, 1]^2} \nabla^2 f B_{m,n} d\mu d\nu \\
    = \iint_{[0, 1]^2} \left(
    \frac{\partial^2 f}{\partial \mu^2} 
    + \frac{\partial^2 f}{\partial \nu^2} 
    \right) e^{-i2\pi m\mu} e^{-i2\pi n\nu} d\mu d\nu.
    \label{eqn:laplacian-coeffs-expanded}
  \end{gather}
  Since the integrand is a sum of partial derivatives, we now have 2 integrals. Notice that inside each integral we have an expressions of the form:
  \begin{equation} \label{eqn:inner-by-parts}
    \int_{[0, 1]} \frac{\partial^2 f}{\partial \xi^2} e^{-i2\pi x\xi} d\xi,
  \end{equation}
  once with \(x = m, \xi = \mu\) and the second time with \(x = n, \xi = \nu\).
  The integral \eqref{eqn:inner-by-parts} can be integrated by parts twice resulting in this ugly expression:
  \[
    e^{-i2\pi x\xi} \left(
    \frac{\partial f}{\partial \xi} - i2\pi f
    \right)\Bigg|_0^1
    + (i2\pi x)^2 \int_{[0, 1]} f e^{-i2\pi x\xi} d\xi.
  \]
  However, actually this is not too bad. That is because once we substitute the bounds two things happen: the exponential in the front always equals 1 and what is inside of the parenthesis can be rewritten as
  \[
    \frac{\partial f}{\partial\xi}(1) - \frac{\partial f}{\partial\xi}(0)
    + i2\pi \left[ f(1) - f(0) \right],
  \]
  which equals zero, since \(f\) and its derivative are continuous and periodic. Hence, we are left with two integrals, that when substituted back into \eqref{eqn:laplacian-coeffs-expanded} give:
  \begin{align*}
    \langle \nabla^2 & f, B_{m, n} \rangle = \\
    &(i2\pi m)^2 \iint_{[0, 1]^2} f e^{-i2\pi m\mu} e^{-i2\pi n\nu} d\mu d\nu \\
    &\quad + (i2\pi n)^2 \iint_{[0, 1]^2} f e^{-i2\pi m\mu} e^{-i2\pi n\nu} d\mu d\nu.
  \end{align*}
  Finally, to complete the proof we rewrite the right side using the compact notation:
  \[
    (i2\pi m)^2 \langle f, B_{m, n} \rangle
    + (i2\pi n)^2 \langle f, B_{m, n} \rangle.
    \qedhere
  \]
\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}

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}

\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
\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}