Prove approximation given by the physicist Max Born In an old book about optics, I have found a nice approximation, that for large l one has: $$P_l(\cos(\theta))  \sim \sqrt{\frac{2}{l \pi \sin(\theta)}} \sin \left((l+\frac{1}{2}) \theta + \frac{\pi}{4} \right)$$, where $P_l$ is the l-th Legendre polynomial. 
Unfortunately, I do not have any clue how to prove this, but maybe somebody here has an idea and an error approximation for this approximation? A reference would be sufficient too, of course. 
 A: You have full asymptotic expansion here, formulas 18.15.12-18.15.13. It can be deduced from the integral representation for Legendre polynomials (formula 18.10.2 here). 
However, I think that Born himself deduced this formula by studying the asymptotics of partial wave solutions of the 2D Helmholtz equation.
A: I recommand the approximation
$$\tag{1}P_n(\cos\vartheta)\approx\sqrt{\dfrac{\vartheta}{\sin\vartheta}}\left\{J_0\left(\sqrt{n(n+1)}\,\vartheta\right)-\dfrac{\vartheta+\dfrac{1}{\vartheta}-\dfrac{\cos\vartheta}{\sin\vartheta}}{8\sqrt{n(n+1)}}Y_0\left(\sqrt{n(n+1)}\,\vartheta\right)\right\}$$
in the interval $\left[0,\dfrac{\pi}{2}\right]$ with the Bessel functions $J_0$ and $Y_0$. The advantage of this approximation is that the phase of the oscillation of the Legendre polynomial close to $\vartheta=0$ corresponds to that of the approximation (1) while
$$\sin\left(\sqrt{n(n+1)}\vartheta+\dfrac{\pi}{4}\right)\approx\sin\left(\left(n+\dfrac{1}{2}\right)\vartheta+\dfrac{\pi}{4}\right)$$
will only be in phase with the Legendre polynomial in some interval $\left[\delta,\dfrac{\pi}{2}\right]$, $\delta>0$. You will get your approximation with the approximation
$$\tag{2}J_0(z)\approx\sqrt{\dfrac{2}{\pi z}}\cos\left(z-\dfrac{\pi}{4}\right)$$
and omitting the second term with $Y_0$ in the approximation (1).
In order to substantiate my approximation I begin with the Laplace operator $\Delta=\dfrac{\partial^2}{\partial x^2}+\dfrac{\partial^2}{\partial y^2}+\dfrac{\partial^2}{\partial z^2}$ that occurs often in physics for wave propagation. Part of the solution in cylindrical coordinates are the Bessel functions. They describe the wave propagation around the cylindrical axis. The Legendre polynomial describes the wave propagation in spherical coordinates on the sphere around the spherical axis of these spherical coordinates. Therefore it does not come with much surprise that the Bessel functions and the Legendre polynomials are connected.
The differential equation of the Legendre polynomial $P_n(x)$ is
$$(1-x^2) y'' - 2x y' + n(n+1)y = 0.$$
The substitution $x=\cos\vartheta$ gives
$$\tag{3}y''+\dfrac{\cos\vartheta}{\sin\vartheta}y'+n(n+1)y=0.$$
You will find the procedure of a variable change in differential equations at this question on math.stackexchange. We establish for the function $w(\vartheta)$
$$w=\sqrt{\dfrac{\vartheta}{\sin\vartheta}}$$
$$w'=\dfrac{1}{2}\sqrt{\dfrac{\vartheta}{\sin\vartheta}} \left\{\dfrac{1}{\vartheta}-\dfrac{\cos\vartheta}{\sin\vartheta}\right\}$$
$$w''=\dfrac{1}{4}\sqrt{\dfrac{\vartheta}{\sin\vartheta}} \left\{\dfrac{2+\cos^2\vartheta}{\sin^2\vartheta}-\dfrac{1}{\vartheta^2}-2\dfrac{\cos\vartheta}{\vartheta\sin\vartheta}\right\}$$
and substitute $y=wz$ in equation (3) and obtain after dividing through w
$$z''+\dfrac{1}{\vartheta}z'+\left\{\dfrac{1}{4}\left(1+\dfrac{1}{\sin^2\vartheta}-\dfrac{1}{\vartheta^2}\right)+n(n+1)\right\}z=0.$$
The term $1+\dfrac{1}{\sin^2\vartheta}-\dfrac{1}{\vartheta^2}$ is bounded and small in comparison to the term $n(n+1)$ for large n and this will not change after the substitution $\vartheta=\dfrac{\varphi}{\sqrt{n(n+1)}}$ and the division through $n(n+1)$:
$$\tag{4}z''+\dfrac{1}{\varphi}z'+\left\{ 1+\dfrac{s(\varphi)}{4n(n+1)} \right\}z=0$$
with
$$\tag{5}s(\varphi)=1+\dfrac{1}{\sin^2\dfrac{\varphi}{\sqrt{n(n+1)}}}-\dfrac{1}{\left(\dfrac{\varphi}{\sqrt{n(n+1)}}\right)^2}.$$
The differential equation of the Bessel functions $J_0(x)$ and $Y_0(x)$ is
$$\tag{6}y''+\dfrac{1}{x}y'+y=0$$
so that the substitution $z(\varphi)=z_s(\varphi)+J_0(\varphi)$ in equation (4) should lead to a sufficiently small function $z_s(\varphi)$ for large n:
$$\tag{7}z_s''+\dfrac{1}{\varphi}z_s'+\left\{ 1+\dfrac{s(\varphi)}{4n(n+1)}\right\}z_s=-\dfrac{s(\varphi)}{4n(n+1)}J_0(\varphi).$$
The Wronskian of this differential equation is
$$W(\varphi)=const\cdot e^{\left( -\int {\dfrac{1}{\varphi} \, d\varphi  } \right)}=\dfrac{const}{\varphi}.$$
Hence a solution $z_s(\varphi)$ is
$$\tag{8}z_s(\varphi)={-\int_0^\varphi \dfrac{z_1(\varphi)z_2(t)-z_2(\varphi)z_1(t)}{W(t)}\,\cdot\,\dfrac{s(t)}{4n(n+1)}\,J_0(t)\,dt} ={\dfrac{const}{4n(n+1)}\,\int_0^\varphi  t\,\left\{z_1(\varphi)z_2(t)-z_2(\varphi)z_1(t)\right\}\,s(t)\,J_0(t)\,dt}$$
with the 2 independant solutions $z_1(\varphi)$, $z_2(\varphi)$ of the differential equation (4). The estimations of the 2 independant solution $ P_n(\cos\vartheta)$, $Q_n(\cos\vartheta)$ of the differential equations (3)
$$ \vert P_n(\cos\vartheta) \vert \le \dfrac {2}{\sqrt {n\pi \sin \vartheta}}$$
$$ \vert Q_n(\cos\vartheta) \vert \le \sqrt {\dfrac {\pi}{n \sin \vartheta}}$$
are known for the interval $0<\vartheta<\pi$ and $n>0$ (formula 8, chapter 37, Kugelfunktionen, Lense). The relationship $y=wz$ renders
$$ \vert z_{1,2} (\cos\varphi) \vert \le \dfrac {const}{\sqrt { \varphi}}.$$
With a bounded function $s(\varphi)$ and the estimation (2) this gives in equation (8) 
$$\tag{9} {\left\vert \dfrac{1}{W(t)}\,z_1(\varphi)z_2(t)\,s(t)\,J_0(t)\, \right\vert } \le 
{ \left\vert \dfrac{t}{const} z_1(\varphi) \dfrac{const}{\sqrt {t}} \cdot const \cdot \dfrac{const}{\sqrt {t}} \right\vert } =
{ const \cdot \left\vert z_1(\varphi) \right\vert }.$$
The second term leads to a similar estimation with the approximation
$$Y_0(z)\approx\sqrt{\dfrac{2}{\pi z}}\sin\left(z-\dfrac{\pi}{4}\right).$$
The integral (8) can therefore be written to
$$z_s(\varphi)=\dfrac{const}{4n(n+1)}\,\{h_1(\varphi)z_1(\varphi)+h_2(\varphi)z_2(\varphi)\}$$
with functions $\vert h_1(\varphi)\vert$, $\vert h_2(\varphi)\vert$ < $const\cdot\varphi$ so that the function $z_s(\varphi)$ is also bounded. With $0\le\varphi\le\sqrt{n(n+1)}\dfrac{\pi}{2}$ we finally obtain
$$\vert z_1(\varphi)-J_0(\varphi)\vert \le 
{\dfrac{const}{4n(n+1)}\,\{const\cdot\varphi \cdot \dfrac {const}{\sqrt { \varphi}} + const\cdot\varphi \cdot \dfrac {const}{\sqrt { \varphi}}\}} =
{const\cdot\dfrac{\sqrt{\varphi}}{4n(n+1)}} \le
{\dfrac{const}{n^{3/2}}}
.$$
A similar estimation can be reached for $z_2(\varphi)=z_{\hat s}(\varphi)+Y_0(\varphi)$.
Though this is a good result these 2 solutions can be taken to estimate the functions $z_s(\varphi)$, $z_{\hat s}(\varphi)$ once again. We get now
$$\tag{9} z_1(t)\,J_0(t)\, \approx J_0(t)\,J_0(t)$$
and we can take half the value of the amplitude in the function $J_0(t)$ because the oscillation only gives a small contribution to the integral. This leads finally to the approximation (1) that is an excellent one. The error for n=4 without the term $Y_0$ is less than $0.02$, the error with the term $Y_0$ is less than $0.003$. The following link supplies you with a picture of this error: 
difference of the approximation (1) above with $J_0$ (green) and with $J_0, Y_0$ (red) for n=4
A: A reference is e.g. Abramowitz and Stegun, 8.10.7 with $\mu=0, \nu=l,$ and $\cos \left((l+\frac{1}{2}) \theta - \frac{\pi}{4} \right)$ replaced by $\sin \left((l+\frac{1}{2}) \theta + \frac{\pi}{4} \right)$.
