How should we think about Spherical Harmonics? Studying Quantum Mechanics I only thought about Spherical Harmonics $Y_{l,m}(\theta , \phi)$:
$$Y_{l,m}(\theta , \phi)=N_{l,m}P_{l,m}(\theta)e^{im\phi}$$
as the simultaneous eigenfunctions of $L_z$ and $L^2$.
But then I stumbled on these two statements:

The Spherical Harmonics are a complete orthonormal basis for the space of the functions defined on the sphere.


$P_{l,0}$, called Legendre's polynomials, are a complete orthonormal basis on the circle.

This two statements have been stated to me without any further explanation.
Given my current physical and mathematical background, as an undergraduate student, I struggle to understand these two; but they seem to me really crucial to comprehend the concepts of angular momentum in QM. I have a feel on what these two statements are trying to tell me, but I don't understand them precisely at all.
Are this two statements really crucial as I think or can I neglect them and continue thinking about $Y$ as eigenfunctions and nothing more without any repercussions?
Is there a simple way to explain what these two statements mean? (By simple way I mean an undergraduate level explanation, without assuming knowledge about Group Theory and such)
 A: 
The Spherical Harmonics are a complete orthonormal basis
for the space of the functions defined on the sphere.

This statement means:
You can decompose any arbitrary function $f(\theta,\phi)$
defined on the sphere (i.e. for $0\le\theta\le\pi$ and $0\le\phi\le 2\pi$)
into a series of the spherical harmonics $Y_{l,m}(\theta,\phi)$:
$$f(\theta,\phi)=\sum_{l=0}^{\infty}\sum_{m=-l}^{+l}f_{l,m}Y_{l,m}(\theta,\phi)$$
The word "complete" in the statement above means,
the function $f(\theta,\phi)$ can completely (i.e.
without any remainder) decomposed into spherical harmonics,
so that you actually have an equality in the formula above.
The word "orthonormal" in the statement above means,
the functions $Y_{l,m}(\theta,\phi)$ satisfy
the following orthonormality relations:
$$\int_0^{2\pi} d\phi \int_0^\pi d\theta\ \sin\theta
 \ Y_{l,m}(\theta,\phi)\ Y_{l',m'}^*(\theta,\phi)
 =\delta_{l,l'}\delta_{m,m'}
 =\begin{cases}
 1, &\text{ if } l=l' \text{ and } m=m' \\
 0  &\text{ else }
 \end{cases}$$
Using these orthonormality relations you can prove
the following formula to find all the coefficients
$f_{l,m}$ from the function $f(\theta,\phi)$:
$$f_{l,m}=\int_0^{2\pi}d\phi \int_0^{\pi}d\theta\ \sin\theta\ f(\theta,\phi)\ Y_{l,m}^*(\theta,\phi)$$
You can find more information about this in Wikipedia section "Spherical harmonics expansion".


$P_{l,0}$, called Legendre's polynomials, are a complete orthonormal basis on the circle.

This statement means a similar thing:
You can decompose any function $g(\theta)$ defined on the circle
(i.e. for $0\le\theta\le\pi$)
into a series of the Legendre polynomials $P_{l,0}(\theta)$:
$$g(\theta)=\sum_{l=0}^\infty g_l P_{l,0}(\theta)$$
"Complete" and "orthonormal" have similar meanings as above.
And you can find the coefficients $g_l$ from the function
$g(\theta)$ by
$$g_l=\int_0^\pi d\theta\ \sin\theta\ g(\theta)\ P_{l,0}^*(\theta)$$
A: As knowing your current knowledge I would like to a very basic interpretation of the following:

$P_{l,0}$ called Legendre's polynomials, are a complete orthonormal basis on the circle.

The other one just follows from here. First I assume you know what the basis vector means. You might like to recap 3D basis vectors.
Now suppose functions which belong to $L_2(-1,1)$, that are functions that square-integrable between ranges of $-1$ to $1$, the range suggest that the functions are basically the cosine whose range is the same as these functions.
$$\int_{-1}^1 dx |f(x)|^2 < \infty$$
Let's restrict ourselves to polynomials. The simplest set I can think of is
$$|\psi_0\rangle,|\psi_1\rangle \cdots \rightarrow x^0,x^1,x^2,\cdots$$
We want to build on are an orthonormal basis. The simplest basis that we can choose from the above monomials are $f_0(x)=c_0$, $f_1(x)= a_1 x+b_1$ First  So  Let's start imposing the conditions :
$$\int_{-1}^{1} |f_0(x)|^2dx=\int_{-1}^{1}dx \ c_0^2=1 \Rightarrow c_0=\frac{1}{\sqrt{2}}$$
Now I need to impose the following condition on $f_1(x)$
$$\langle \phi_i|\phi_j\rangle=\delta_{ij}$$
$$\int_{-1}^{1}dx \ f_0(x)f_1(x)=0 \Rightarrow f_1(x)=a_1 x$$
$$\int_{-1}^{1}dx \ |f_1(x)|^2=1 \Rightarrow  a_1=\sqrt{\frac{3}{2}}$$
$$\Rightarrow  f_1(x)=\sqrt{\frac{3}{2}}x$$
Next $f_2(x)=a_2x^2+b_2x+c_2$ then you have three conditions that exactly equals to number of unknown. From you can build this series.  Now do you recognize these function, let me change the normalization constant and notation a little bit:
$$\int_{-1}^{1}dx \ P_n(x)P_m(x)=\frac{2}{2n+1}\delta_{nm}$$
These are indeed Legendre Polynomials.
The normalization is due to historical reasons to get rid of those awful factors.
$$P_0(x)=1, \ \ \ \ P_1(x)=x, \ \ \ P_2(x)=\frac{3x^2-1}{2}\cdots$$
So these objects act as unit vectors in the space of function that belongs to $L_2(-1,1)$. That's the special case so you can make your own basis (changing the range or normalization factor).
In general,
$$\int_{-\infty}^{\infty} d\mu(x) \ \phi_n^*(x)\phi_m(x)=A_n\delta_{nm}$$
where $d\mu(x)$ is a measure that insure that function goes to zero ( sufficient rapidly as $|x|\rightarrow \infty$). Now If $d\mu(x)=dx e^{-x^2}$, then these are Hermite polynomial. Sometimes you would like to work with radial coordinate in that case $d\mu(x)=dx \ e^{-x}$ that laguree polynomials.

Now Why circle? That you already reorganize by Yourself. Why sphere? That to You can reorganize yourself. That's just a matter of ranges.
A: 
How should we think about Spherical Harmonics?

In short: In the same way that you think about plane waves.

Spherical harmonics (just like plane waves) are basic, essential tools. As such, they are used over a very wide variety of contexts, and each of those contexts will end up using a different "way to think" about the tools: using different properties of the tool, putting it in different frameworks, and so on.
In that regard, "how we should think" about the tool depends on what we want to do with it and what we already know about the context. This is ultimately a personal question that is up to you, and you haven't really given us much to go on to say anything useful.
That said, you seem to be confused about the relationship between the statements

*

*

spherical harmonics are eigenfunctions of the angular momentum

and


*

spherical harmonics form a basis for the functions over the sphere.

As I said above, you should think of these in the same way that you think of plane waves, for which the two statements above have direct analogues:

*

*

plane waves are eigenfunctions of the linear momentum

and


*

plane waves form a basis for the functions over the line,

in the more specific senses that

*

*the plane wave $e^{ikx}$ is an eigenfunction of $\hat p = -i\hbar \frac{\partial}{\partial x}$ with eigenvalue $k$
and


*every function $f(x)$ over real $x$ can be expressed as a superposition of plane waves, by means of its Fourier transform: $f(x) = \int_{-\infty}^\infty \tilde f(k) e^{ikx} \mathrm dk$.


So: which of these two properties is more important? it depends! it depends on what you're doing and what you care about. Maybe you need to combine the two on equal footings, maybe you need to put more weight on one or the other.
Can you forget about the second statement and just focus on the eigenfunction properties? maybe! it depends on what you're doing and what you care about. If you continue analyzing the problem, both properties will come into play eventually, but where and when depends on the path you take.
(As a rule of thumb, if a statement is unclear and it is not being actively used in the context you're in, then: yes, it is probably safe to put a pin on it and move on, and to return to understanding what it means and why it holds only when you run into it being actively used. In many cases, actually seeing it used in practice is a massive help into understanding what it's for!)

In any case, though, the completeness properties can indeed seem a little mysterious. They've been handled quite well by the other answers, so I won't belabour that point. Instead, I will fill you in on a secret that very few textbooks will tell you: the spherical harmonics are all polynomials!
More specifically, they are polynomials of the Cartesian coordinates when the position $\mathbf r$ is restricted to the unit sphere. Once you take away all of the trigonometric complexity and so on, the $Y_l^m(\theta,\phi)$ become simple polynomials:
\begin{align}
Y_0^0(\mathbf r) & = 
1
\\
Y_1^{\pm 1}(\mathbf r) & = 
x\pm i y
\\
Y_1^0(\mathbf r) & = 
z
\\
Y_2^{\pm 2}(\mathbf r) & = 
(x\pm i y)^2
\\
Y_2^{\pm 1}(\mathbf r) & = 
z (x\pm i y)
\\
Y_2^0(\mathbf r) & = 
x^2+y^2-2 z^2
\end{align}
(where I've removed the pesky normalization constants).
The basic concept of building this family of polynomials is quite simple:

*

*Work your way up in degree, from constant to linear to quadratic to higher powers.

*Keep the whole set linearly independent.

*Within each degree, keep things as invariant as possible under rotations about the $z$ axis.

That last part might sound mysterious, but it should be relatively easy to see why it forces a preference for the combinations $x\pm iy$: if you rotate by an angle $\alpha$ in the $x,y$ plane, $(x\pm iy)$ transforms simply, to $e^{\pm i\alpha}(x\pm iy)$. If you're wondering, this is the feature that connects to the $Y_l^m$ being eigenfunctions of $\hat L_z$.
With that in place, it is fairly easy to see how the progression goes
$$
1,x\pm iy, z, (x\pm iy)^2, z(x\pm iy),\ldots
$$
... but then what's with the $x^2+y^2-2 z^2$ combination? The answer to that is that, in general, there are six separate quadratic monomials we need to include:
$$
x^2,y^2,z^2,xy,xz,yz.
\tag{*}
$$
We have already covered some of these in $(x\pm i y)^2 = x^2-y^2 \pm 2ixy$, and we've fully covered the mixed terms $xy,xz,yz$, so now we need to include two more: to keep things symmetric, let's say $x^2+y^2$ and $z^2$.
Except, here's the thing: the pure-square monomials $x^2,y^2,z^2$ in $(*)$ are not linearly independent! Why is this? well, because we're on the unit sphere, which means that these terms satisfy the identity
$$
x^2+y^2+z^2 = 1,
$$
and the constant term $1$ is already in our basis set. So, from the combinations $x^2+y^2$ and $z^2$ we can only form a single polynomial, and the correct choice turns out to be $x^2+y^2-2z^2$, to make the set not only linearly independent but also orthogonal (with respect to a straightforward inner product).
As for the rest of the spherical harmonics $-$ the ladder keeps climbing, but the steps are all basically the same as the ones I've outlined already.

Anyways, I hope this is sufficient to explain why the spherical harmonics are

a complete orthonormal basis for the space of the functions defined on the sphere.

This is just a fancy way of saying that, if you have a function $f(\mathbf r)$ that's defined for $|\mathbf r|=1$, then you can expand $f(\mathbf r)$ as a "Taylor series" of suitable polynomials in $x$, $y$ and $z$, with a slightly reduced set of polynomials to account for the restrictions on $|\mathbf r|$.
A: How familiar are you with Fourier analysis?
For instance, if you want a complete basis of functions defined on the unit square, you would pick:
$$ |n, m\rangle \equiv e^{2\pi i nx}e^{2\pi i my}$$
with all $n, m \in \mathbb Z$
Note that this is just a tensor product of the complete orthogonal functions on a 1 dimensional line segment: $[0,1]$.
Physics-wise, they are eigenfunctions of the wave equation on a plate.
So we want to do something similar on a sphere, where the square area element
$$ \mathrm dA = \mathrm dx\:\mathrm dy $$
is replaced by:
$$ \mathrm dA = \sin{\theta}\mathrm d\theta \:\mathrm d\phi = \mathrm d(\cos{\theta})\mathrm d\phi $$
So the $\phi$ part is easy, it's the same as a straight-line, but its support is $[0,2\pi]$, so the functions are:
$$ e^{im\phi} $$
The polar angle is different, it's not uniform. You wind up with the Legendre polynomials:
\begin{align}
P_0(z) & = 1 \\
P_1(z) & = z \\
P_2(z) & = \frac 1 2 (3z^2-1) \\
& \ \ \vdots \\
P_l(z) & = \frac 1 l[(2l-1)zP_{l-1}(z) - (l-1)P_{l-2}(z)]\\
& \ \ \vdots 
\end{align}
When you combine the two, you get the spherical harmonics:
$$ |l,m\rangle = Y_l^m(\theta, \phi)$$
Physics-wise, they are the solutions to the wave equation on a sphere, and that's it: any function on a sphere can be expanded into a series of spherical harmonics, just as any function on a unit square can be expanded into plane waves in $x$ and $y$.
The reason the plane waves on a plane are easier is that these are eigenfunctions of momentum, and momentum generates translations, and translations commute.
Spherical harmonics are eigenfunctions of angular momentum, and angular momentum generates rotations, and rotations do not commute.
From there, you need to get into group theory and representation theory for a deeper understanding.
Without doing that, I think the most important thing to understand is that at fixed $l$, the spherical harmonics are closed under rotations. The means, rotations mix the $m$ functions, but not the different $l$ values.
This means each $l$ has a geometric meaning on the sphere. For $l=0$:
$$ Y_0^0(\theta, \phi) = \frac 1 {\sqrt{4 \pi}} $$
is rotationally invariant. It is uniform, a.k.a.: the 'monopole', a.k.a. a perfect sphere.
The 3 $Y_1^m$ are linear combinations of $x$, $y$, and $z$, which of course are closed under rotations...it's just $(x,y,z)$ which transform like a vector. They also quantify the dipole moment of a distribution on a sphere.
The 5 $Y_2^m$ describe the various quadrupole moments. $Y_2^0$ describes cylindrically symmetric ($m=0$) oblateness/prolateness, while $Y_2^{\pm 2}$ is a elongation/squishing versus longitude.  The $Y_2^{\pm 1}$ can be diagonalized away by the correct choice of axes.
From there you move onto higher multipole moments to describe any shape. The Earth's gravitational field, for instance, is not perfectly spherical, and the current best model of it, EGM-08, is a spherical harmonic expansion to degree $l=2190$ and order $|m|=2159$, which resolves features around 19 km.
