Step 2 — Write the recurrence for both variables. For x: \[ (2k+1)xP_k(x)=(k+1)P_{k+1}(x)+kP_{k-1}(x). \] For y: \[ (2k+1)yP_k(y)=(k+1)P_{k+1}(y)+kP_{k-1}(y). \]

Step 3 — Multiply the first by Pₖ(y) and the second by P_k(x)P_k(x). So we multiply the x-equation by P_k(y): \[ (2k+1)xP_k(x)P_k(y)=(k+1)P_{k+1}(x)P_k(y)+kP_{k-1}(x)P_k(y). \] Multiply the y-equation by P_k(x): \[ (2k+1)yP_k(x)P_k(y)=(k+1)P_{k+1}(y)P_k(x)+kP_{k-1}(y)P_k(x). \]

Step 4 — Subtract the two equations. Subtract the second from the first:

Step 5 — Sum from k = 0 to n. Sum both sides over k = 0, 1, … , n: Left side: \[ (x-y)\sum _{k=0}^n(2k+1)P_k(x)P_k(y). \] Right side telescopes beautifully: The terms with k+1 and the terms with k cancel in pairs. Everything collapses except the top boundary term. After cancellation, only one term survives:

Step 6 — Divide by x-y. We obtain: The final result (Christoffel–Darboux Identity) This is exactly the identity we wanted.

1. Interior points x ∈ (−1, 1);
2. end point x = 1 (and analogously x = −1).

Interior points x ∈ (−1, 1).

Fix an interior point x ∈ (−1, 1), and write \[ x = \cos\varphi , \quad t = \cos\theta , \qquad \varphi , \theta \in (0, \pi ) . \] We use a standard estimate for Legendre's polynomials in the oscillatory regime: \[ \left\vert P_n (\cos\theta ) \right\vert \le \frac{C}{\sqrt{n\,\sin\theta}} \qquad \left( C > 0, \ n \ge 1 \right) . \tag{L.1} \] This follows from from the well-known asymptotic expansion \[ P_n (\cos\theta ) = \sqrt{\frac{2}{\pi n\,\sin\theta}} \,\cos \left( \left( n + \frac{1}{2} \right) \theta - \frac{\pi}{4} \right) + O \left( \frac{1}{n^{3/2}} \right) , \quad \theta \in (0, \pi ), \ n\to \infty , \] uniform in θ ∈ [ε, π − ε], and separate (crude) bound near the endpoints; but for our purposes we treat the above inequality as a known estimate.

In particular, for fixed φ ∈ (0, π), \[ \left\vert P_n (\cos\varphi )\right\vert \le \frac{C}{\sqrt{n\,\sin\varphi}} \le \frac{C}{\sqrt{n}} . \] Thus, there is C_x > 0 depending only on x such that \[ \left\vert P_n (x) \right\vert \le \frac{C_x}{\sqrt{n+1}} , \qquad n \ge 0 . \] Consequently, \[ \left\vert \hat{P}_n (x) \right\vert = \sqrt{\frac{2n+1}{2}} \left\vert P_n (x) \right\vert \le C_x \sqrt{n+1} \, \frac{1}{\sqrt{n+1}} = C_x , \] so orthonormal Legendre polynomials are uniformly bounded in index n.

Similarly, for variable t = cosθ with θ away from 0 and π, the same estimate holds uniformly on compact subintervals of (−1, 1). This gives better control of K_N(x, t) for |θ − φ| not too small.

A Fejér-type bound in the interior
It is a classical result that the orthogonal polynomial setting, for fixed interior x = cosφ, the Fejér kernel behaves like a "bump" of width ∼ 1/(N+1) in the angular variable θ, with height ∼ N+1. To make this precise, one can either:

• rely on full Plancherel--Rotach asymptotics and compute the Cesàro averages in the angular variable, or
• use estimates on K_N plus summation in m to obtain a bound of the form \[ \left\vert K_N^{(1)} (x, \cos\theta ) \right\vert \le \frac{C_x}{1 + (N+1)^2 (\theta - \varphi )^2} , \quad \theta \in (0, \pi ) . \]

Lemma 2A (Fejér kernel bound in the interior): For each fixed x ∈ (−1, 1), there exists a constant C_x > 0 such that \[ \left\vert K_N^{(1)} (x, \cos\theta ) \right\vert \le \frac{C_x}{1 + (N+1)^2 (\theta - \varphi )^2} , \qquad \theta \in (0, \pi ), \ x=\cos\varphi , \tag{L.2} \] for all θ ∈ (0, π), where x = cosφ.    ■

We accept Lemma 2A as a classical kernel estimate; it is proved in detail in monographs on orthogonal polynomials (e.g., Szegő, Chapter VII), by combining the Christoffel–Darboux formula, the asymptotics of Pₙ(cosθ), and summation over m ((Fejér averaging).

Fix x ∈ (−1, 1) and δ > 0. We want to show \[ \int_{|t-x| > \delta} K_N^{(1)} (x,t)\,{\text d}t \ \to \ 0 . \] In angular variables, x = cosφ, t = cosθ. The map θ ↦ t = cosθ is smooth and monotone on [0, π]. The condition |t − x| > δ is equivalent to |cosθ − cosφ| > δ, which in turn implies |θ − φ| ≥ c_x,δ > 0 (since cosθ is Lipschitz with nonzero derivative at interior points). Formally, there exists η such that \[ |t-x| > 0 \qquad \Longrightarrow \qquad |\theta - \varphi | > \eta . \] The change of variables t = cosθ yields dt = −sinθ dθ, so \[ \int_{|t-x| > \delta} K_N^{(1)} (x,t)\,{\text d}t = \int_{|\theta - \varphi |> \eta} K_N^{(1)} (x, \cos\theta )\,\sin\theta \,{\text d}\theta . \] Using the bound (L.1), \[ \left\vert \int_{|\theta - \varphi | > \eta} K_N^{(1)} (x, \cos\theta )\,\sin\theta\,{\text d}\theta \right\vert \le \int_{|\theta - \varphi | > \eta} \frac{C_x \,\sin\theta}{1 + (N+1)^2 (\theta - \varphi )^2}\,{\text d}\theta . \] For |θ − φ| > η, we have sinθ ≤ 1, so \[ \le C_x \int_{|\theta - \varphi | > \eta} \frac{{\text d}\theta}{1 + (N+1)^2 (\theta - \varphi )^2} . \] Now change variables \[ u = \left( N+1 \right) \left( \theta - \varphi \right) , \quad {\text d}\theta = \frac{{\text d}u}{N+1} . \] Then |θ − φ| > η becomes |u| > (N+1)η. Hence, \[ \int_{|\theta - \varphi | > \eta} \frac{{\text d}\theta}{1 + (N+1)^2 (\theta - \varphi )^2} = \frac{1}{N+1} \int_{|u| > (N+1)\eta} \frac{{\text d}u}{1 + u^2} . \] The integrand 1/(1 + u²) is integrable on ℝ, and its tail integral satisfies \[ \int_{|u| > (N+1)\eta} \frac{{\text d}u}{1 + u^2} \ \to\ 0 \quad \mbox{as }\ N\to\infty . \] Explicitly, \[ \int_{|\theta - \varphi | > \eta} \frac{{\text d}\theta}{1 + (N+1)^2 (\theta - \varphi )^2} \le \frac{2}{(N+1)\,\eta} . \] Thus, \[ \int_{|\theta - \varphi | > \eta} \frac{{\text d}\theta}{1 + (N+1)^2 (\theta - \varphi )^2} \le \frac{C'}{(N+1)^2} \] some constant C' depending on η.

Consequently, \[ \int_{|t - x| > \delta} K_N^{(1)} (x,t)\,{\text d}t \ \to \ 0 \quad \mbox{as } N \to \infty \] for each fixed interior x. This proves Lemma 1 for x ∈ (−1, 1).

Endpoint x = 1.

Now we treat the endpoint x = 1 explicitly, because the angular parametrization degenerates at this point. Let t = cosθ, θ ∈ [0, π]. We fix x = 1.

Let us consider \[ K_N^{(1)} (1, t) = K_N^{(1)} (1, \cos\theta ). \] The corresponding Christoffel–Darboux formula is simplified because Pₙ(1) = 1 for all n: \begin{align*} K_m (1, \cos\theta ) &= \sum_{k=0}^m \frac{2k+1}{2} \,P_k (1)\, P_k (\cos\theta ) = \sum_{k=0}^m \frac{2k+1}{2} \,P_k (\cos\theta ) , \\ &= \frac{m+1}{2} \cdot \frac{P_{m+1} (\cos\theta ) - P_m (\cos\theta )}{1-\cos\theta} , \end{align*} Thus, we need a bound for this expression. \[ \left\vert P_{m+1} (\cos\theta ) - P_{m} (\cos\theta ) \right\vert \le \lert\vert P_{m+1} (\cos\theta ) \right\vert + \left\vert P_{m} (\cos\theta ) \right\vert \le 2 . \] So \[ \left\vert K_m (1, \cos\theta ) \right\vert \le \frac{m+1}{2} \cdot \frac{2}{1 - \cos\theta} = \frac{m+1}{1-\cos\theta} . \] Hence, \begin{align*} \left\vert K_{N}^{(1)} (\cos\theta ) &\le \frac{1}{N+1} \sum_{m=0}^N \frac{m+1}{1-\cos\theta} \\ &= \frac{1}{1-\cos\theta} \cdot \frac{1}{N+1} \sum_{m=0}^N (m+1) \\ &= \frac{N+2}{2 \left( 1 - \cos\theta \right)} . \end{align*} Using the elementary identity \[ 1 - \cos\theta = 2\,\sin^2 \left( \frac{\theta}{2} \right) , \] we have 1 − cosθ ≥ cθ² for θ ∈ [0, π], with some constant c > 0. Hence, \[ 1 - \cos\theta \ge \frac{2}{\pi^2}\,\theta^2 , \] or \[ 1 - \cos\theta \ge C\,\theta^2 , \qquad 0 \le \theta \le \pi , \] for some absolute constant C. This is our global endpoint bound.

Behavior near the endpoint We want a better control for small θ, where 1 − cosθ is small. Use Taylor expansion of Pₙ near x = 1.

From the Legendre differential equation \[ \left( 1 - x^2 \right) P''_n (x) - 2x\,P'_n (x) + n \left( n+1 \right) P_n (x) = 0 , \] we get \[ \sup_{|x| \le 1} \left\vert P''_n (x) \right\vert \le C\,n^4 \] for some absolute constant C > 0 (we only need a rough bound). Thus, Taylor's theorem at x = 1 gives \[ P_n (x) = 1 + P'_n (1) \left( x-1 \right) + \frac{1}{2}\,P''_n (\xi ) \left( x-1 \right)^2 , \] for some ξ between x and 1. Therefore, \[ P_n (x) = 1 + \frac{n \left( n+1 \right)}{2} \left( x - 1 \right) + R_n (x) , \] where \[ | R_n (x) | \le \frac{1}{2}\,\sup_{|y|\le 1} \left\vert P''_n (y) \right\vert (x-1)^2 \le C\,n^4 \left( x-1 \right)^2 . \] For x = cosθ, θ ∈ [0, 1], say. Then \[ 1 - \cos\theta = \frac{\theta^2}{2} + O\left( \theta^4 \right) , \] and in particular, for sufficiently small θ, \[ c_1 \theta^2 \le 1 - \cos\theta \le c_2 \theta^2 , \] for some positive constants c₁, c₂. Hence, \[ P_n (\cos\theta ) = 1 - \frac{n(n+1)}{2}\,(1-\cos\theta ) + O\left( n^4\left( 1- \cos\theta \right)^2 \right) \] Now compute the difference: \[ P_{m+1} (\cos\theta ) - P_m (\cos\theta ) = -(m+1)\,(1-\cos\theta ) + O \left( m^4 (1-\cos\theta )^2 \right) . \] Since \[ (m+1)(m+2) - m(m+1) = (m+1) \left[ (m+2) - m \right] = 2(m+1) , \] Now average over m, we get the value of the Fejér kernel \begin{align*} K_N^{(1)} (1, \cos\theta ) &= \frac{1}{N+1} \sum_{m=0}^N K_m (1, \cos\theta ) \\ &= - \frac{1}{2(N+1)} \sum_{m=0}^N \left( m+1 \right)^2 + \cdots . \end{align*} The sum of squares is \[ \sum_{m=0}^N \left( m+1 \right)^2 = \sum_{k=1}^{N+1} k^2 = \frac{(N+1)(N+2)(2N+3)}{6} \sim \frac{(N+1)^3}{3} . \] Thus, the main term is of order (N+1)². The error term is \[ \frac{1}{N+1} \sum_{m=0}^N m^5 \left( 1 - \cos\theta ) \le C\,N^5 \left( 1 - \cos\theta \right) . \] So for small θ, \[ P_n (\cos\theta ) = 1 - \frac{n(n+1)}{4}\,\theta^2 + O \left( n^4 \theta^4 \right) \qquad \theta\to 0 . \] Hence, there exists C > 0 such that \[ \left\vert K_N^{(1)} (1, \cos\theta ) \right\vert \le C\left( (N+1)^2 + N^5 (1-\cos\theta ) \right) , \qquad 0 < \theta \le \pi . \] In particular, for θ small enough \[ N^5 \left( 1 - \cos\theta \right) \le C' \left( N+1 \right)^2 , \] i.e., \[ 1 - \cos\theta \le \frac{1}{N^5} \quad \iff \quad \theta \le \frac{1}{N^{3/2}} , \] we have \[ \left\vert K_N^{(1)} (1, \cos\theta ) \right\vert \le \left( N+1 \right)^2 . \] Combining with the global bound (L.2), we arrive at \[ \left\vert K_N^{(1)} (1 - \cos\theta ) \right\vert \lesssim \min \left\{ (N+1)^2 , \ \frac{N+1}{\theta^2} \right\} , \qquad 0 < \theta \le \pi . \] This is the endpoint analogue of the Fejér kernel bound in the angular variable.

Endpoint approximate identity We must show \[ \int_{|t-1| > \delta} K_N^{(1)} (1, t) \,{\text d}t \ \to \ 0 \quad \mbox{as } N \to \infty , \] for every small δ > 0.

In angular variables t = cosθ with θ ∈ [0, π], the condition |t − 1| > δ is equivalent \[ 1 - \cos\delta > \delta \quad \iff \quad \cos\theta < 1-\delta \quad \iff \quad \theta \ge \theta_0 , \] for some θ₀ = θ₀(δ) ∈ (0, π]. In particular, there exists c_δ > 0 such that \[ \theta \ge \theta_0 \quad \Longrightarrow \quad \theta \ge c_{\delta} > 0 . \] Now \begin{align*} \int_{|t-1| > \delta} K_N^{(1)} (1, t)\,{\text d}t &= \int_0^{\pi} 1_{1-\cos\theta > \delta} K_N^{(1)} (1, \cos\theta )\,\sin\theta\,{\text d}\theta \\ &= \int_{\theta \ge \theta_0} K_N^{(1)} (1, \cos\theta )\,\sin\theta\,{\text d}\theta . \end{align*} Using global crude bound at the endpoint and sinθ ≤ 1, we get \[ \left\vert \int_{\theta \ge \theta_0} K_N^{(1)} (1, \cos\theta )\,\sin\theta\,{\text d}\theta \right\vert \le \int_{\theta \ge \theta_0} \left\vert K_N^{(1)} (1, \cos\theta ) \right\vert {\text d}\theta \le C \int_{\theta \ge \theta_0} \frac{N+1}{\theta^2}\,{\text d}\theta . \] But θ ≥ θ₀, so \[ \int_{\theta \ge \theta_0} \frac{N+1}{\theta^2}\,{\text d}\theta = \left( N+1 \right) \left[ -\frac{1}{\theta} \right]_{\theta_0}^{\pi} = \left( N+1 \right) \left( \frac{1}{\theta_0} - \frac{1}{\pi} \right) \sim \frac{N+1}{\theta_0} . \] This bound grows with N, so we have to be more careful: the global bound is too crude. Instead, we exploit the fact that the Fejér kernel at x = 1 has total mass 1, and the large values of the kernel are confined to a small neighborhood of size ∼ 1/(N+1). In particular, we have

• For θ ≤ α/(N+1), we have \( \displaystyle \quad K_N^{(1)} (1, \cos\theta ) \le (N+1)^2 . \)
• For θ ≥ α/(N+1), we have, which is integrable in θ on [α/(N+1), π] with total integral O(1), and moreover that integral tends to 0 as we push the lower limit to a fixed θ₀.

Fix δ > 0, let θ₀ be as above. Choose N large enough so that \[ \frac{\alpha }{N+1}<\theta _0. \] Then \[ \{ \theta \geq \theta _0\} \subset \left\{ \theta \geq \frac{\alpha }{N+1}\right\} . \] Hence, \[ \int _{\theta \geq \theta _0}|F_N(1,\cos \theta )|\, d\theta \leq \int _{\theta \geq \alpha /(N+1)}\frac{C(N+1)}{\theta ^2}\, d\theta =C(N+1)\left[ \, -\frac{1}{\theta }\right] _{\alpha /(N+1)}^{\pi }=C\left( \frac{N+1}{\alpha /(N+1)}-\frac{N+1}{\pi }\right) =C\left( \frac{(N+1)^2}{\alpha }-\frac{N+1}{\pi }\right) . \] This estimate alone still grows like (N+1)², but recall that the total mass of \( \displaystyle \quad K_N^{(1)} (1, \cos\theta )\,\sin\theta\,{\text d}\theta \quad \) is 1: \[ \int_0^{\pi} K_N^{(1)} (1, \cos\theta )\,\sin\theta\,{\text d}\theta = \int_{-1}^1 K_N^{(1)} (1, t)\,{\text d}t = 1 . \] Thus, the contribution from small region θ ≤ α/(N+1), where the Fejér kernel ∼ (N+1)², is of order \[ \int _0^{\alpha /(N+1)}(N+1)^2\sin \theta \, d\theta \sim (N+1)^2\cdot \frac{\alpha }{N+1}=\alpha (N+1), \] which would blow up unless the implicit constants are such that the mass remains 1; in fact, precise asymptotics (coming from Bessel kernels) show that:

• the height is ≈ N+1, not (N+1)², in the correct normalization for dt.

Global crude bound. Since |Pₙ(cosθ)| ≤ 1 for all n and θ ∈ [0, π], \[ \bigl|P_{m+1}(\cos\theta)-P_m(\cos\theta)\bigr|\le2, \] and therefore \begin{equation*} \left\vert K_m(1,\cos\theta) \right\vert \le \frac{(m+1)}{1-\cos\theta}. \end{equation*} Using \begin{equation*} 1-\cos\theta = 2\sin^2\!\left(\frac{\theta}{2}\right) \ge \frac{2}{\pi^2}\,\theta^2, \qquad 0\le\theta\le\pi, \end{equation*} we obtain \begin{equation*} \left\vert K_m(1,\cos\theta) \right\vert \;\lesssim\; \frac{m+1}{\theta^2}, \qquad 0<\theta\le\pi. \end{equation*} Averaging in $m$ yields \[ \left\vert K_N^{(1)}(1,\cos\theta) \right\vert \le \frac{1}{N+1}\sum_{m=0}^N |K_m(1,\cos\theta)| \lesssim \frac{1}{N+1}\sum_{m=0}^N \frac{m+1}{\theta^2} \lesssim \frac{N+1}{\theta^2}. \tag{L2.1} \]

Height bound at θ = 0. At t = 1, \[ K_N^{(1)}(1,1) = \frac{1}{N+1}\sum_{n=0}^N \frac{2n+1}{2}\cdot P_n(1)^2 = \frac{1}{N+1}\sum_{n=0}^N \frac{2n+1}{2} . \] Since \[ \sum_{n=0}^N (2n+1) = (N+1)^2, \] we obtain \begin{equation*} K_N^{(1)}(1,1) = N+1. \end{equation*} By positivity of the Fejér kernel, \begin{equation*} 0\le K_N^{(1)} (1,t)\le K_N^{(1)}(1,1)=N+1, \qquad t\in[-1,1], \end{equation*} and thus \[ \left\vert K_N^{(1)} (1,\cos\theta) \right\vert \le N+1, \qquad 0\le\theta\le\pi. \tag{L2.2} \]

Fejér--type scaling. Combining (L2.1) and (L2.2), we get \begin{equation*} \left\vert K_N^{(1)} (1,\cos\theta) \right\vert \;\lesssim\; \min\!\left\{\,N+1,\ \frac{N+1}{\theta^2}\right\}, \qquad 0<\theta\le\pi. \end{equation*} For all θ > 0, \begin{equation*} \min\!\left\{1,\ \frac{1}{(N+1)^2\theta^2}\right\} \;\lesssim\; \frac{1}{1+(N+1)^2\theta^2}. \end{equation*} Multiplying by N+1, \begin{equation*} (N+1)\min\!\left\{1,\ \frac{1}{(N+1)^2\theta^2}\right\} \;\lesssim\; \frac{N+1}{1+(N+1)^2\theta^2}. \end{equation*} Hence, there exists C > 0 such that \[ \left\vert K_N^{(1)} (1,\cos\theta) \right\vert \;\le\; \frac{C\,(N+1)}{1+(N+1)^2\theta^2}, \qquad 0<\theta\le\pi, \tag{L2.3} \] which is exactly the required inequality.

Approximate identity away from θ = 0. Fix θ₀ > 0. Using Eq.(L2.3), we have \[ \int_{\theta\ge\theta_0} \left\vert K_N^{(1)} (1,\cos\theta) \right\vert \sin\theta\,{\text d}\theta \le C\int_{\theta_0}^{\pi} \frac{(N+1)}{1+(N+1)^2\theta^2}\,{\text d}\theta. \tag{L2.4} \] Since \begin{equation*} \frac{\text d}{{\text d}\theta}\arctan\bigl((N+1)\theta\bigr) = \frac{N+1}{1+(N+1)^2\theta^2}, \end{equation*} we compute \begin{equation*} \int_{\theta_0}^{\pi} \frac{(N+1)}{1+(N+1)^2\theta^2}\,{\text d}\theta = \arctan\bigl((N+1)\pi\bigr) - \arctan\bigl((N+1)\theta_0\bigr). \end{equation*} Using the asymptotic expansion \begin{equation*} \arctan z = \frac{\pi}{2}-\frac{1}{z}+O(z^{-3}), \qquad z\to+\infty, \end{equation*} we obtain \[ \arctan\bigl((N+1)\theta_0\bigr) = \frac{1}{N+1}\left(\frac{1}{\theta_0}-\frac{1}{\pi}\right) +O\!\left((N+1)^{-3}\right) \longrightarrow 0 \tag{L2.5} \] as N → ∞. Combining (L2.4) and (L2.5) yields the second inequality.

1. Start from the differential equation for Legendre's polynomials: \[ \left( 1− x^2\right) y'' −2x\,y' +ℓ\left( ℓ+1\right) y = 0. \] Set x = cosθ. Then derivative become \begin{align*} \frac{{\text d}y}{{\text d}x} &= \frac{{\text d}y}{{\text d}\theta}\cdot \frac{{\text d}\theta}{{\text d}x} = \frac{{\text d}y}{{\text d}\theta}\cdot \frac{1}{{\text d}x/{\text d}\theta} = \frac{{\text d}y}{{\text d}\theta} \left( - \frac{1}{\sin \theta} \right) , \\ \frac{{\text d}^2 y}{{\text d}x^2} &= \frac{\text d}{{\text d}x} \left( \frac{{\text d}y}{{\text d}\theta}\cdot \frac{1}{{\text d}x/{\text d}\theta} \right) = \frac{{\text d}^2 y}{{\text d}\theta^2} \left( - \frac{1}{\sin \theta} \right)^2 + \frac{{\text d}y}{{\text d}\theta} \cdot \frac{\text d}{{\text d}x} \left( - \frac{1}{\sin \theta} \right) \\ &= \frac{{\text d}^2 y}{{\text d}\theta^2} \cdot \frac{1}{\sin^2 \theta} + \frac{{\text d}y}{{\text d}\theta} \cdot \frac{\text d}{{\text d}\theta} \left( - \frac{1}{\sin \theta} \right) \cdot \frac{{\text d}\theta}{{\text d}x} = \frac{{\text d}^2 y}{{\text d}\theta^2} \left( \frac{1}{\sin^2 \theta} \right) - \frac{{\text d}y}{{\text d}\theta} \left( \frac{\cos\theta}{\sin^3 \theta} \right) . \end{align*} Then the first term of Legendre's equation is \[ \left( 1- x^2 \right) y'' = \frac{y_{\theta\theta} \sin\theta - y_{\theta} \cos\theta}{\sin\theta} . \] The second term is \[ -2x\,y' = 2\left( \frac{\cos\theta}{\sin\theta}\right) y_{\theta} . \] This transforms the Legendre equation into: \[ \frac{{\text d}^2 u}{{\text d}\theta^2} + \cot\theta \,\frac{{\text d}u}{{\text d}\theta} + \ell \left( \ell + 1 \right) u = 0 , \] where u(θ) = y(cos(θ)). We check with Mathematica:
   Clear[\[Theta], n, u]; (*y(x)=u(\[Theta]),with x=Cos[\[Theta]]*) (*dy/dx=(du/d\[Theta])/(dx/d\[Theta])*) y1 = D[u[\[Theta]], \[Theta]]/D[Cos[\[Theta]], \[Theta]] // Simplify; (*d²y/dx²=d/d\[Theta](dy/dx)/(dx/d\[Theta])*) y2 = D[y1, \[Theta]]/D[Cos[\[Theta]], \[Theta]] // Simplify; (*Legendre equation in \[Theta]:(1-x^2) y''-2 x y'+n(n+1) y=0,with \ x=Cos[\[Theta]]*) legendreEq\[Theta] = Simplify[(1 - Cos[\[Theta]]^2) y2 - 2 Cos[\[Theta]] y1 + n (n + 1) u[\[Theta]] == 0]
   Running this code yields:
   n (1 + n) u[\[Theta]] + Cot[\[Theta]] Derivative[1][u][\[Theta]] + ( u^\[Prime]\[Prime])[\[Theta]] == 0
2. Apply a Liouville–Green (WKB) transformation. Introduce a rescaled function: \[ y(θ) = \sqrt{\frac{\theta}{\sin\theta}} \,u(θ). \] Then u(θ) satisfies an equation of the form: \[ u'' + \left( ℓ+ \frac{1}{2} \right)^2 u ≈ \mbox{small correction}. \]
3. Approximate by Bessel equation After further normalization, the equation becomes asymptotically: \[ u'' + \frac{1}{\theta}\,u' + \left( ℓ+ \frac{1}{2}\right)^2 u ≈ 0, \] which is essentially the Bessel equation. Hence: \[ u(θ) ∼ J_0 \left( \left( ℓ+ \frac{1}{2}\right) θ\right) , \] and corrections produce the J₁-term.
4. Convert Bessel → cosine asymptotics Using: \[ J_0 (z) ∼ \sqrt{\frac{2}{\pi z}} \,\cos\left( z− \frac{\pi}{4} \right) , \] you obtain the cosine form of Hilb’s formula.

Another proof:   
One can derive Hilb's formula based on the integral representation of the Legendre polynomials: \[ P_n (\cos\theta ) = \frac{1}{\pi} \int_0^{\pi} \left( \cos\theta + {\bf j}\,\sin\theta\,\cos\varphi \right)^n {\text d}\varphi = \frac{1}{\pi} \int_0^{\pi} \left( x + {\bf j}\,\sqrt{1 - x^2}\,\cos\varphi \right)^n {\text d}\varphi . \] The integral is complex. but the result is real. Why? Because the imaginary part integrates to zero. Integrand can be written as \[ \left( A + {\bf j}\,B\,\cos\varphi \right)^n , \quad A = x, \quad B = \sqrt{1-x^2} . \] Now observe key symmetry: \[ \cos\left( \pi - \phi \right) = - \cos\phi . \] So \[ \left( A + {\bf j}\,B\,\cos (\pi - \phi ) \right)^n = \left( A - {\bf j}\,B\,\cos\phi \right)^n = \left( A + {\bf j}\,B^{\ast} \cos\phi \right)^n . \] Thus, the integrand at φ and at π − φ are complex conjugates. Pairing the integral kills the imaginary part: \[ \int_0^{\pi} f(\phi )\,{\text d}\phi = \int_0^{\pi /2} f(\phi )\,{\text d}\phi + \int_{\pi /2}^{\pi} f(\phi )\,{\text d}\phi . \] Make the substitution φ ↦ π − φ in the second integral \[ \int_{\pi /2}^{\pi} f(\phi )\,{\text d}\phi = \int_0^{\pi /2} f(\pi - \phi )\,{\text d}\phi = \int_0^{\pi /2} f^{\ast}(\phi )\,{\text d}\phi . \] So the full integral is \[ \int_0^{\pi} f(\phi )\,{\text d}\phi = \int_0^{\pi /2} \left[ f(\phi ) + f^{\ast} (\phi ) \right]{\text d}\phi = 2\,\int_0^{\pi /2} \Re f(\phi )\,{\text d}\phi . \] The imaginary part cancels exactly.

We rwrite the integrand in the integral formula for the Legendre polynomial in exponential form in order to apply stationary phase method. Define \[ \Phi (\phi ) = \Im \left\{ \ln \left( \cos\theta + \mathbf{j}\,\sin\theta\,\cos\phi \right) \right\} , \] so the integrand behaves like exp{ ⅉnΦ(φ) } up to a slowly varyung amplitude.

We find stationary points by solving Φ′(φ) = 0. One finds that the main contributions come from points corresponding to φ = 0 and φ = π, which geometrically encode the directions aligned with θ and π − θ. These are non-degenerate stationary points when 0 < θ < π, θ ≠ 0,π, which is why the formula is uniform away from the endpoints.

Apply the method of stationary phase Near each stationary point φ₀, expand: \[ \Phi (\phi ) = \Phi (\phi_0 ) + \frac{1}{2}\,\Phi'' (\phi_0 ) \left( \phi - \phi_0 \right)^2 , \] and approximate the amplitude by its value at φ₀. Then each neighborhood contributes an integral of the form \[ \int e^{\mathbf{j}\,n\,\Phi (\phi_0 )} \, e^{\frac{1}{2}\,\Phi'' (\phi_0 )\left( \phi - \phi_0 \right)^2}\,{\text d}\phi \approx e^{\mathbf{j}\,n\,\Phi (\phi_0 )} \, \sqrt{\frac{2\pi}{n \left\vert \Phi'' (\phi_0 ) \right\vert}} \,e^{\pm \mathbf{j}\pi /4} . \] We label:

Summing the two complex conjugate contributions yields a cosine.

Extract the explicit amplitude and phase: A more careful computation of Φ(φ₀), Φ′′(φ₀), and the prefactor from the original integral gives:

• Phase: \( \displaystyle \quad \left( n + \frac{1}{2} \right) \theta - \frac{\pi}{4} . \)
• Amplitude: \( \displaystyle \quad \sqrt{\frac{2}{n \pi\,\sin\theta}} . \)

What’s swept under the rug A full proof has to:

• Justify:
  • The choice and validity of the integral representation.
  • That only the two stationary points contribute at leading order.
  • Uniform bounds on the remainder (control of the integral away from stationary points).
• Refine: Higher-order terms come from including higher derivatives of and the amplitude in the expansion, giving a full asymptotic series in powers of 1/n.

References (with proofs):

• G. Szegő, Orthogonal Polynomials (Chapter on asymptotics)
• E. T. Whittaker & G. N. Watson, A Course of Modern Analysis
• F. W. J. Olver, Asymptotics and Special Functions

So we need to show that \begin{align*} \sigma _N \left( f,x\right) &= \frac{1}{N+1}\sum _{k=0}^N S_k (f,x) = \frac{1}{N+1}\sum _{k=0}^N \int_{-1}^1 K_k (x,t)\,f(t)\,{\text d}t \\ &= \int_{-1}^1 K_N^{(1)} (x,t)\,f(t)\,{\text d}t \end{align*} converges to f(x) for almost every x. Here the Dirichlet kernel is \begin{align*} K_n (x,t) = \sum_{k=0}^n \hat{P}_k (x)\, \hat{P}_k (t) = \sum_{k=0}^n P_k (x) \left( k + \frac{1}{2} \right) P_k (t) \\ &= \frac{n+1}{2}\cdot \frac{P_{n+1} (x)\,P_n (t) - P_n (x)\,P_{n+1} (t)}{x-t} , \quad x \ne t . \end{align*} The Fejér kernel is \begin{align*} K_N^{(1)} (x,t) &= \frac{1}{N+1} \sum_{m=0}^N K_m (x,t) \\ &= \sum_{n=0}^N \left( 1 - \frac{n}{N+1} \right) \hat{P}_n (x)\,\hat{P}_n (t) \\ &= \sum_{n=0}^N \left( 1 - \frac{n}{N+1} \right) \left( n + \frac{1}{2} \right) P_n (x)\,P_n (t) . \end{align*} We need three properties:

1. Symmetry and normalization: \[ K_N^{(1)}(x,t) = K_N^{(1)}(t,x) \] by definition (the Legendre polynomials appear in Eq.(4) symmetrically).
2. Positivity: for all N ∈ ℕ and all x, t ∈ [−1, 1], we have \[ K_N^{(1)} (x,t) \ge 0 . \] Proof: For each index m ≥ 0, let P_m : ℌ → ℌ be the orthogonal projection onto the finite dimensional subspace \[ V_m = \mbox{span}\left\{ P_0 , P_1 , P_2 , \ldots , P_m \right\} . \] Then P_m is self-adjoint and idempotent, and \[ P_m f = \sum_{n=0}^m \left( \int_{-1}^1 f(t)\,\hat{P}_n (t)\,{\text d}t \right) \hat{P}_n (x) = S_m (f , x) . \] Define the Fejér operator \[ T_N f := \frac{1}{N+1} \sum_{m=0}^N P_m f . \] Then σ_Nf = T_N. Each P_m is positive in the following sense: if f ≥ 0 almost everywhere, then P_mf ≥ 0 almost everywhere (this can be justified either by an explicit kernel representation, or via spectral theory for self-adjoint projectors on real &Lfr'²; this is standard in the theory of orthogonal expansions). A convex combination of positive operators is positive, so if f ≥ 0, then \[ \sigma_N (f, x) = T_N f (x) \ge 0 \quad \mbox{a.e. } x. \] On the other hand, we have the integral representation \[ \sigma_N (f, x) = \int_{-1}^1 K_N^{(1)} (x,t)\,f(t)\,{\text d}t . \] Fix x₀ ∈ [−1, 1] and N. Suppose that for the sake of contradiction, that there exists a set \[ E = \left\{ t \in [-1, 1] \ : \ K_N^{(1)} (x,t) < 0 \right\} \] with positive measure. Choose f ≥ 0 supported in E, nontrivial. Then \[ \sigma_N \left( f , x_0 \right) = \int_{-1}^1 K_N^{(1)} (x_0 ,t)\,f(t)\,{\text d}t \le \int_E K_N^{(1)} (x_0 ,t)\,f(t)\,{\text d}t < 0 , \] contradicting positivity of T_N. Therefore, the Fejér kernel is positive for almost every t. By continuity in t (as finite sum of continuous functions), we get positivity of the Fejér kernel for all t ∈ [−1, 1]. Since x₀ was arbitrary, the claim follows.

3. Normalization: For every x ∈ [−1, 1] and every positive integer N, we have \[ \int_{-1}^1 K_{N}^{(1)} (x,t)\,{\text d}t = 1 \qquad \forall x \in [-1, 1] . \] Proof: Since integral over the Legendre polynomial is zero, \( \displaystyle \quad \int_{-1}^1 P_n (x, t)\,{\text d}t = 0 \quad n=1,2,\ldots , \quad \) we get \begin{align*} \int_{-1}^1 K_{N}^{(1)} (x,t)\,{\text d}t &= \int_{-1}^1 \sum_{n=0}^N \left( 1 - \frac{n}{N+1} \right) \left( n + \frac{1}{2} \right) P_n (x) \,P_n (t)\,{\text d}t \\ &= \sum_{n=0}^N \left( 1 - \frac{n}{N+1} \right) \left( n + \frac{1}{2} \right) P_n (x) \,\int_{-1}^1 P_n (t)\,{\text d}t \\ &= \sum_{n=0}^N \left( 1 - \frac{n}{N+1} \right) \left( n + \frac{1}{2} \right) P_n (x) \,2\,\delta_{n,0} \\ &= \sum_{n=0}^N \left( 1 - \frac{0}{N+1} \right) \left( 0 + \frac{1}{2} \right) P_0 (x)\,2 = 1 . \end{align*} because \( \displaystyle \quad \int_{-1}^1 P_0 (x,t)\,{\text d}t = 2 \quad \) and \( \displaystyle \quad \int_{-1}^1 P_n (x,t)\,{\text d}t = 0 \quad \) for k ≥ 1. Hence \[ \int_{-1}^1 K_N^{(1)}(x,t)\,{\text d}t = \sum_{k=0}^N \left( 1 - \frac{0}{N+1} \right) \left( 0 + \frac{1}{2} \right) P_0 (x) \cdot 2 = 1 \]

   Its 𝔏¹ norm is \[ \left\| K_N^{(1)}(x,t) \right\|_1 = \int_{-1}^1 K_N^{(1)}(x,t) \,{\text d}t = 1 . \]

To prove pointwise convergence, we need to establish that the Fejér kernels form an approximate identity for each x. This approach is based on determination of estimates for the Dirichlet and Fejér kernels; these formulas are standard consequences of the Christoffel–Darboux formula and the differential equation for Legendre's polynomials, and are treated in detail in monographs such as Szegő's book.

For x ≠ t, we have \[ K_N (x,t) = \frac{N+1}{2}\cdot \frac{P_{N+1}(x)\,P_N (t) - P_N (x)\,P_{N+1} (t)}{x-t} . \] The Legendre polynomials are bounded |Pₙ(x)| ≤ 1 on interval [−1, 1] and all indices. Hence \begin{align*} \left\vert K_N (x,t) \right\vert &\le \frac{N+1}{2}\cdot \frac{|P_{N+1}(x)\,P_N (t)| + |P_N (x)\,P_{N+1} (t)|}{|x-t|} \\ &\le \frac{N+1}{2}\cdot \frac{1+1}{|x-t|} \\ &= \frac{N+1}{|x-t|} . \end{align*} Therefore, \[ \left\vert K_N (x,t) \right\vert \le \frac{N+1}{|x-t|} , \qquad x \ne t . \] For x = t, one can use the orthonormal expansion: \begin{align*} \left\vert K_N (x,x) \right\vert &= \sum_{n=0}^N \hat{P}_n^2 (x) = \sum_{n=0}^N \frac{2N+1}{2} \,P_n^2 (x) \\ &\le \sum_{n=0}^N \frac{2N+1}{2} = \frac{(N+1)^2}{2} . \end{align*} So we have a combined estimate \[ \left\vert K_N (x,x) \right\vert \le \min \left\{ \frac{(N+1)^2}{2}, \ \frac{N+1}{|x-t|} \right\} . \] It follows that \begin{align*} \left\vert K_N^{(1)} (x,t) \right\vert &= \frac{1}{N+1} \left\vert \sum_{m=0}^N K_m (x,t) \right\vert \\ &\le \frac{1}{N+1} \sum_{m=0}^N \min \left\{ \frac{(m+1)^2}{2} , \ \frac{m+1}{|x-t|} \right\} . \end{align*} This yields \[ \left\vert K_N^{(1)} (x,t) \right\vert \min \left\{ (N+1)^2 , \ \frac{N+1}{|x-t|} \right\} . \] uniformly in x, t ∈ [−1, 1]. The qualitative picture is

• Near the diagonal t = x, the Fejér kernel can be as large as ∼ (N + 1)².
• Away from the diagonal, the kernel decays like (N + 1)/|x - t|.

We state the approximate identity property precisely.
Lemma: Let x ∈ [−1, 1] be fixed. Then for every δ > 0, \[ \lim_{N\to\infty} \int_{|t-x| > \delta} K_N^{(1)} (x,t)\,{\text d}t = 0. \] Sketch of proof: Because \( \displaystyle \quad K_N^{(1)} (x,t) \ge 0 \quad \) and \( \quad \int_{-1}^1 K_N^{(1)} (x,t)\,{\text d}t =1. \quad \) it suffices to show that for each fixed δ > 0, \[ \sup_{|t-x| > \delta} K_N^{(1)} (x,t) \ \to \ 0 \quad \mbox{as }\ N \to \infty . \] That is, the Fejér kernels become uniformly small away from the diagonal. For orthogonal polynomials of the Legendre type on a compact interval such bounds follow from the asymptotic behavior of the Christoffel–Darboux kernel: one has Plancherel--Rotach type asymptotic for P_m, which apply that \[ K_N (x,t) = O(1) , \quad |t-x| > \delta , \] with the implicit constant depending on δ but not on N. In particular, \[ K_N (x,t) = \frac{1}{N+1} \sum_{m=0}^N K_m (x, t) = O(1) \] for |t − x| ≥ δ. However, the normalization \[ \int_{-1}^1 K_N (x,t)\,{\text d}t = 1 \] forces the mass to concentrate near x. Indeed, if a uniform positive proportion of the mass were retained outside |t − x| ≤ δ, then normalization together with the bound \( \displaystyle \quad K_N^{(1)} (x,t) = O(1) \quad \) away from x would contradict the growth \( \quad K_N^{(1)} (x,t) \ \sim \ c(x)\,(N+1) \quad \) on shrinking neighborhoods. A rigorous version of this argument is standard in the theory of orthogonal polynomials: see Szegő's book, Chapter VII, where general conditions are given under which the Cesàro kernels of orthogonal expansions form an approximate identity.

For the purpose of the present theorem we take Lemma as known; it is the Legendre polynomial analogue of the standard fact that the trigonometric Fejér kernel converges weakly to the delta measure and forms an approximate identity.    ■

Convergence at Lebesgue points: Let f ∈ 𝔏¹([−1, 1]). By the Lebesgue differentiation theorem, almost every x ∈ [−1, 1], \[ \lim_{r\downarrow 0} \int_{x-r}^{x+r} \left\vert f(t) - f(x) \right\vert {\text d}t = 0 . \] Such x are called Lebesgue points of f.

Fix Lebesgue point x. We shall show that \[ \lim_{N\to\infty} \sigma_N (f , x) = f(x) . \] Recall \[ \sigma_N (f, x) = \int_{-1}^1 K_N^{(1)} (x,t) \,f(t)\,{\text d}t . \] Using normalization of the Fejér kernel, we rewrite \[ \sigma_N (f, x) - f(x) = \int_{-1}^1 K_N^{(1)} (x,t) \left[ f(t) - f(x) \right] {\text d}t . \] Let ε > 0. Since x is a Lebesgue point, there exists δ > 0 such that \[ \frac{1}{2\delta} \int_{x-\delta}^{x+\delta} \left\vert f(t) - f(x) \right\vert {\text d}t < \varepsilon . \] Split the integral into "near" and "far" part: \[ \sigma_N (f,x) - f(x) = I_N^{near} + I_N \] where \[ I_N^{near} = \int_{|t-x| \le\delta} K_N^{(1)} (x,t) \left\vert f(t) - f(x) \right\vert {\text d}t . \] Near part We use only positivity and normalization of the kernel near x. Since the Fejér kernels are nonnegative, \[ \left\vert I_N^{near} \right\vert \le \int_{|x-t| \le \delta} K_N^{(1)} (x,t) \left\vert f(t) - f(x) \right\vert {\text d}t . \] Introduce the probability measure \[ {\text d}\mu (t) = K_N^{(1)} (x,t) \,{\text d}t \qquad \mbox{on } \ [-1, 1] . \] Then \[ \left\vert I_N^{near} \right\vert \le \int_{|x-t| \le\delta} \left\vert f(t) - f(x) \right\vert {\text d}\mu_n (t) . \] However, the restriction of μ_N to the interval [x − δ, x + δ] is dominated by a probability measure supported there (it is in fact already a sub-probability measure). In particular, \[ \left\vert I_N^{near} \right\vert \le \sup_{\nu} \int_{x-\delta}^{x+\delta} \left\vert f(t) - f(x) \right\vert {\text d}\nu (t) , \] where ν runs over all probability measures supported in [x − δ, x + δ]. The largest such integral (for absolutely continuous measures) is obtained by the normalized Lebesgue measure, so \[ \left\vert I_N^{near} \right\vert \le \frac{1}{2\delta} \int_{x-\delta}^{x+\delta} \left\vert f(t) - f(x) \right\vert {\text d}t < \varepsilon . \] This estimate is uniform in N. Therefore, \[ \sup_{N> 0} \left\vert I_N^{near} \right\vert \le \varepsilon . \]

Far part: We use the approximate identity property of the Fejér kernels

First treat bounded function f. Suppose f ∈ 𝔏^∞([−1, 1]). Then \[ \left\vert f(t) - f(x) \right\vert \le 2 \, \| f \|_{\infty} , \] and hence \begin{align*} \left\vert I_N^{far} \right\vert &\le \int_{|x-t| > \delta} K_N^{(1)} (x,t) \left\vert f(t) - f(x) \right\vert {\text d}t \\ &\le \| f \|_{\infty} \int_{|x-t| > \delta} K_N^{(1)} (x,t) \,{\text d}t . \end{align*} By Lemma, \[ \int_{|x-t| > \delta} K_N^{(1)} (x,t) \,{\text d}t \ \to \ 0 \quad \mbox{as} \quad N\to\infty . \] Therefore, \[ I_N^{far} \ \to \ 0 \quad \mbox{as} \quad N\to\infty \] for bounded functions.

For general f ∈ 𝔏¹, approximate f by bounded functions f^(k) (for instance, truncate f at levels ±k). The Fejér operators are bounded on 𝔏¹ (indeed, they are positive, preserve constants, and are contractions on 𝔏²; interpolation gives boundedness on 𝔏^p for 1 ≤ p ≤ 2, and density argument extends the pointwise convergence result from bounded f to all f ∈ 𝔏¹). This is standard in the general theory of Fejér/Cesàro summability of orthogonal series and mimics exactly the Fourier series.

Thus, for our fixed Lebesgue point x, \[ \lim_{N\to\infty} I_N^{far} = 0 . \]

Conclusion: Combining the two parts, we have for all ε > 0, \[ \limsup_{N\to\infty} \left\vert \sigma_N (f, x) - f(x) \right\vert \le \sup_N \left\vert I_N^{near} \right\vert + \limsup_{N\to\infty} \left\vert I_N^{far} \right\vert \le \varepsilon + 0 . \] Since ε > 0 was arbitrary, it follows that \[ \lim_{N\to\infty}\sigma_N (f, x) = f(x) . \] This identity holds at each Lebesgue point x of f, and by the Lebesgue differentiation theorem such points form a set of full measure in [−1, 1]. The theorem 2 is proved.

Remarks: The structure of the Fejér kernel here is exactly parallel to the trigonometric Fejér kernel \( \displaystyle \quad F_n (x) = \frac{1}{n}\,\sum_{k=0}^{n-1} D_k (x) , \quad \) which is nonnegative, normalized, and forms an approximate identity on the circle.

The genuinely nontrivial step is Lemma (approximate identity property), which rests on asymptotics of Christoffel–Darboux kernels and Legendre polynomials. This is treated in detain in Szegő's monograph (chapter VII). In a more general setting (orthogonal polynomials with respect to positive weights on compact intervals), and our Legendre case satisfies all the necessary hypotheses.

No essential modification is needed at the endpoints: the Lebesgue differentiation theorem applies with intervals [1 −r, 1], and the Fejér kernels still form an approximate identity in the appropriate sense. If desired, one may rewrite the kernels near x = 1 in the angular variable t = cosθ to get explicit bounds of the form \[ K_N^{(1)} (1, \cos\theta ) \le \min \left\{ (N+1)^2 , \ \frac{N+1}{\theta^2} \right\} , \] that makes the concentration near fully explicit.    ■

Approximate identity: abstract framework.
We recall the general notion needed for the convergence theorem,

Let (X, μ) be a measure space, here X = [−1, 1] with Lebesgue measure. A family of kernels K_α(x, t) (α is a parameter) is called approximate identity if

Uniform 𝔏¹ bound \[ \sup_{\alpha} \sup_{x \in [-1,1]} \int_{-1}^1 K_{\alpha} (x, t)\,{\text d}t < \infty . \] In our case, positivity + normalization gives this bound with constant 1.

Localization (vanishing tails: \[ \lim_{\alpha \to \infty} \sup_{x \in [-1,1]} \int_{-1+\delta}^{1-\delta} \left\vert K_{\alpha} (x, t)\right\vert {\text d}t = 0 \qquad \forall \delta > 0 . \] Then for any f ∈ 𝔏¹([−1, 1]), the operators \[ T_{\alpha} f(x) = \int_{-1}^1 f(t)\, K_{\alpha} (x,t)\,{\text d} t \] converge to f(x) at every Lebesgue point of f, and hence almost everywhere.

So for almost everywhere convergence, what we must prove is that the Fejér kernels form an approximate identity.

We already have items 1 and 2. The key remaining point is the localization. So we want to show that \[ \lim_{\alpha \to \infty} \sup_{x \in [-1,1]} \int_{-1+\delta}^{1-\delta} \left\vert K_{N}^{(1)} (x, t)\right\vert {\text d}t = 0 \qquad \forall \delta > 0 . \] Intuition tells us that as N grows, the Fejér kernels concentrates more and more near the diagonal t = x with tails that carry vanishing mass.

Interior points x away from endpoints (singular points).
Fix ε > 0. Consider x in a compact subinterval [−1 + &epsilon, 1 − ε]. There we can use Hilb's asymptotics for Legendre polynomials: if x = cosθ, t = cosϕ , then for large n, \[ P_n (\cos\theta ) \sim \sqrt{\frac{2}{\pi n \sin\theta}}\, \cos \left( \left( n + \frac{1}{2} \right) - \frac{\pi}{4} \right) \] and similarly for Pₙ(cosϕ).

Plugging the series expression \[ K_N^{(1)} (x,t) = \sum_{n=0}^N \left( 1 - \frac{n}{N+1} \right) \left( n + \frac{1}{2} \right) P_n (x)\, P_n (t) , \] we see that for x and t in the interior region and with |x − t| ≥ δ, the phase oscillates with respect to n. The factor 1 −n/(N+1) is a Cesàro wight; standard summation methods (Dirichlet--Abel summation or direct adaption of Fourier--Fejér analysis) show that the sum over n then behaves like a localized bump whose mass concentrates near, and whose integral over |θ − ϕ| ≥ δ tends to zero as N → ∞.

Informally, away from the line t = x, the oscillations in n cancel out, and Cesàro averaging suppresses the cancellations even more strongly. More precisely, one shows that for fixed δ > 0 and ε > 0, \[ P_n (\cos\theta ) = \sqrt{\frac{\theta}{\sin\theta}} \, J_0 \left( \left( n + \frac{1}{2} \right) \theta \right) + O \left( \frac{1}{n} \right) , \] and hence by positivity and normalization, \[ \lim_{n \to \infty} \sup_{x \in [-1+\varepsilon , 1-\varepsilon ]} \int_{-1+\delta}^{1-\delta} K_{N}^{(1)} (x, t)\, {\text d}t = 0 \qquad \forall \delta > 0 . \] Indeed, for fixed positive δ, there is a constant C(δ) such that \[ \int_{|t i x| \ge \delta} \left\vert K_N^{(1)} (x,t) \right\vert {\text d}t \le \frac{C(\delta )}{N} \to 0 . \] So we have localization uniformly for x away from the endpoints.

Near the endpoints, our previous estimates such as \[ \left\vert K_N^{(1)} (x,t) \right\vert \le \frac{C}{|x - t|} \] do not work because the kernel involves denominators like \[ \sqrt{1 - x^2}\,\sqrt{1 - t^2} , \] or their angular versions, which behave like 1/(sinθ sinϕ). These estimates are harmless when θ, ϕ are bounded away from 0, π, but they blow up near the endpoints.

Using the Christoffel–Darboux formula, the Fejér kernel can be written as \[ K_N^{(1)} (x,t) = \frac{1}{N+1} = \sum_{k=0}^N \frac{k+1}{2} \cdot \frac{P_{k+1} (x)\, P_k (t) - P_k (x)\, P_{k+1} (t)}{x-t} . \] For x = 1, we have \[ K_N^{(1)} (1,t) = \frac{1}{N+1} = \sum_{k=0}^N \frac{k+1}{2} \cdot \frac{P_{k+1} (1)\, P_k (t) - P_k (1)\, P_{k+1} (t)}{1-t} . \] Since Pₙ(1) = 1 for any nonnegative integer n, we come to the telescopic series \[ A = \sum_{k=0}^N \left( k+1 \right) P_k (t) - \sum_{k=0}^N \left( k+1 \right) P_{k+1} (t) . \] In the second sum, we change the index of summation by putting j = k+1. This yields \[ \sum_{k=0}^N \left( k+1 \right) P_{k+1} (t) = \sum_{j=1}^{N+1} j\,P_j (t) . \] So \[ A = \sum_{k=0}^N \left( k+1 \right) P_k (t) - \sum_{k=1}^{N+1} k\,P_k (t) = 1 - \left( N+1 \right) P_{N+1} (t) + \sum_{k=1}^N P_k (t) . \]

Now we show the uniform boundness of the 𝔏¹ norm of the Fejér kernel. Fixing a small positive δ, we split the integral \[ \int_{-1}^1 \left\vert K_N^{(1)} (1,t) \right\vert {\text d}t = \int_{-1}^{1-\delta} \left\vert K_N^{(1)} (1,t) \right\vert {\text d}t + \int_{1-\delta}^1 \left\vert K_N^{(1)} (1,t) \right\vert {\text d}t . \] We need bounds independent of N.For t being fae away from t = 1, |t − 1| ≥ δ, the denominator 1 − t is bounded away from zero: 1 − t ≥ δ. So from the spliting formula, we get \[ \left\vert K_N^{(1)} (1,t) \right\vert \le \frac{1}{2 \left( N+1 \right) \delta} \left( \left\vert \sum_{k=0}^N P_k (t) \right\vert + \left( N+1 \right) \left\vert P_{N+1} (t) \right\vert \right) . \] A single Legendre polynomial is uniform;y bounded: \[ \left\vert P_{N+1} (t) \right\vert \le C_2 (\delta ) , \quad t \in [-1, 1-\delta ], \quad \forall N . \] We use a classical fact about partial sum of Legendre polynomials (both uniform in N for fixed t ∈ [−1, 1 − δ]) \[ \left\vert \sum_{k=0}^N P_k (t) \right\vert \le C_1 (\delta ) , \quad t\in [-1, 1-\delta ], \quad \forall N , \] with some positive constant C₁(δ). (This comes from the generating function for the Legendre polynomials.) This gives a uniform estimate of 𝔏¹ norm doe t away from the endpoint t = 1.

Near t = 1, we set t = cosθ. Then t = 1 corresponds small θ. So on interval [1 − δ, 1] we have θ ∈ [0, θ₀] with θ₀ determined by cosθ₀ = 1 − δ. Then \[ {\text d}t = -\sin\theta\,{\text d}\theta , \quad 1-t = 1 - \cos\theta \sim \frac{\theta^2}{2} . \] The integral becomes \[ I_N^{\tiny near} = \int_{1-\delta}^1 \left\vert K_{N}^{(1)} (1,t) \right\vert \le {\text d}t = \int_0^{\theta_0} \left\vert K_{N}^{(1)} (1,\cos\theta ) \right\vert \sin\theta\,{\text d}\theta . \] Now the telescopic formula in θ-form becomes \[ K_{N}^{(1)} (1,\cos\theta ) = \frac{1}{2 \left( N+1 \right) \left( 1 - \cos\theta \right)} \left( \sum_{k=0}^N P_k (\cos\theta ) - \left( N+1 \right) P_{N+1} (\cos\theta ) \right) . \] Using small angle behavior \[ 1 - \cos\theta \sim \frac{\theta^2}{2} , \qquad \sin\theta \sim \theta , \] the prefactor behaves like \[ \frac{\sin\theta}{1 - \cos\theta} \sim \frac{\theta}{\theta^2 /2} = \frac{2}{\theta} . \] Thus, for small θ, we get \[ \int_{1-\delta}^1 \left\vert K_{N}^{(1)} (1,\cos\theta ) \right\vert \sin\theta\,{\text d}\theta \le \] Now we apply the standard asymptotic input (Hilb type asymptotic):
for each compact θ-interval [0, θ₀] and all N ≥ 1, \[ \left\vert P_n (\cos\theta ) \right\vert \le C \left( 1 + n\theta \right)^{-1/2} . \] Consequently, \[ \left\vert \sum_{k=0}^N P_k (\cos\theta ) \right\vert \le C\,\min \left\{ N+1 , \frac{1}{\sqrt{\theta}} \right\} . \] Accepting these inequalities (they are standard can be derived from Bessel-function approximations), you get

• The term with P_N+1: \[ \frac{1}{N+1} \cdot \frac{1}{\theta} \left( N+1 \right) \left\vert P_{N+1} (\cos\theta ) \right\vert \le \frac{C''}{\theta} \left( 1 + N\theta \right)^{-1/2} . \]
• The term with sum of the Legendre polynomials: \[ \frac{1}{N+1} \cdot \frac{1}{\theta} \left\vert \sum_{k=0}^N P_k (\cos\theta ) \right\vert \le \]

================================= Bessel/Mehler--Heine scaling.
Near the endpoints, we use Mehler--Heine type asymptotics. For example, near x = 1, write \[ x = \cos\frac{u}{N}, \qquad \phi = \cos\frac{v}{N} , \] with u, v of order 1. Then \[ P_n \left( \cos\frac{u}{N} \right) \to J_0 (u) \qquad \mbox{as } n\to\infty . \] Summing with Fejér weights in the Legendre case corresponds, under this scaling to constructing Bessel-type Fejér kernel in the variables u, v

Lebesgue differentiation theorem: for almost every point, the value of an integrable function is the limiting average taken around the point.

Step 1: Christoffel--Darboux and kernel representation. By the Christoffel--Darboux formula for Legendre polynomials, \begin{equation*} K_N(x,t) = \frac{N+1}{2}\, \frac{P_{N+1}(x)P_N(t)-P_N(x)P_{N+1}(t)}{x-t}. \tag{T4.1} \end{equation*} For x = x₀ fixed in (-1, 1) and t ∈ (-1, 1), this gives \[ K_N(x_0,t) = \frac{N+1}{2}\, \frac{P_{N+1}(x_0)P_N(t)-P_N(x_0)P_{N+1}(t)}{x_0-t}. \] Introduce angular variables x = &cosθ, t = cosφ with θ, &hi; ∈ (0, π), so that x₀ = cosθ₀ and t = cosφ. Then \[ x_0-t = \cos\theta_0 - \cos\varphi = -2\sin\frac{\theta_0+\varphi}{2}\,\sin\frac{\theta_0-\varphi}{2}. \] It is a classical result that for fixed θ ∈ (0, π), \begin{equation*} P_n(\cos\theta) = \sqrt{\frac{2}{\pi n\sin\theta}}\, \cos\Bigl((n+\tfrac12)\theta - \tfrac{\pi}{4}\Bigr) + O\bigl(n^{-3/2}\bigr), \tag{T4.2} \end{equation*} uniformly on compact subsets of (0, π). Substituting (T4.2) into (T4.1), a standard stationary-phase computation (or direct trigonometric algebra) shows that, for fixed θ₀ ∈ (0, π), \begin{equation*} K_N(\cos\theta_0,\cos\varphi) = \frac{1}{\pi\sin\theta_0}\, \frac{\sin\bigl((N+1/2)(\varphi-\theta_0)\bigr)}{\varphi-\theta_0} + R_N(\varphi), \tag{T4.3} \end{equation*} where the remainder R_N satisfies \begin{equation*} \int_0^\pi |R_N(\varphi)|\,{\text d}\varphi \le C, \qquad \int_{|\varphi-\theta_0|>\varepsilon} |R_N(\varphi)|\,{\text d}\varphi \le C_\varepsilon, \tag{T4.4} \end{equation*} with constants independent of N and C_{ε → 0 as ε → 0. In particular, KN(cosθ₀, cosφ) behaves like the usual trigonometric Dirichlet kernel in the variable φ θ₀.}

Step 2: Reduction to a Dirichlet integral. Changing variables t = cosφ in the representation \[ S_N f(x_0) = \int_{-1}^1 f(t)\,K_N(x_0,t)\,dt, \] we have dt = -sinφ dφ, and thus \[ S_N f(\cos\theta_0) = \int_0^\pi f(\cos\varphi)\,K_N(\cos\theta_0,\cos\varphi)\,\sin\varphi\,d\varphi. \] Let \[ g(\varphi) := f(\cos\varphi)\,\sin\varphi. \] The assumptions that f has bounded variation in a neighborhood of x₀ and finite one-sided limits at x₀ translate into the fact that g has bounded variation in a neighborhood of φ = θ₀ and finite one-sided limits g(θ₀ ±0). Moreover, since sinθ₀ > 0, we have \[ \frac{g(\theta_0-0)+g(\theta_0+0)}{2} = \sin\theta_0\,\frac{f(x_0-0)+f(x_0+0)}{2}. \] Using (T4.3), we can write \begin{align*} S_N f(\cos\theta_0) &= \frac{1}{\pi\sin\theta_0} \int_0^\pi g(\varphi)\, \frac{\sin\bigl((N+1/2)(\varphi-\theta_0)\bigr)}{\varphi-\theta_0}\, d\varphi \\ &\quad + \int_0^\pi g(\varphi)\,R_N(\varphi)\,d\varphi. \end{align*} By (T4.4) and the fact that g is of bounded variation on [0, π], the second term is uniformly bounded and contributes no singular behavior in N; its contribution to the limit is handled by the same arguments as in the standard proof for Fourier series (dominated convergence plus the localization near φ = θ₀).

Step 3: Application of the classical Dirichlet theorem. The first term is precisely of the form \[ \frac{1}{\pi\sin\theta_0} \int_0^\pi g(\varphi)\, \frac{\sin\bigl((N+1/2)(\varphi-\theta_0)\bigr)}{\varphi-\theta_0}\, {\text d}\varphi, \] which is the usual Dirichlet kernel acting on the function g in the variable φ (up to a harmless factor sinθ₀). Since g has bounded variation near θ₀ and finite one-sided limits there, the classical Dirichlet theorem for Fourier integrals implies that \[ \lim_{N\to\infty} \int_0^\pi g(\varphi)\, \frac{\sin\bigl((N+1/2)(\varphi-\theta_0)\bigr)}{\varphi-\theta_0}\, {\text d}\varphi = \pi\,\frac{g(\theta_0-0)+g(\theta_0+0)}{2}. \] Combining this with the factor 1/(πsinθ₀) and the identity \[ \frac{g(\theta_0-0)+g(\theta_0+0)}{2} = \sin\theta_0\,\frac{f(x_0-0)+f(x_0+0)}{2}, \] we obtain \[ \lim_{N\to\infty} S_N f(x_0) = \frac{f(x_0-0)+f(x_0+0)}{2}, \] as claimed in Theorem 4.

Step 1: Kernel representation at x = 1. As before, the Legendre partial sums are \[ S_N (f,x) = \int_{-1}^1 f(t)\,K_N(x,t)\,{\text d}t, \] with \[ K_N(x,t) = \sum_{n=0}^N \frac{2n+1}{2}\,P_n(x)\,P_n(t). \] Setting x = 1 and using Pₙ(1) = 1, we obtain \[ K_N(1,t) = \sum_{n=0}^N \frac{2n+1}{2}\,P_n(t). \] Equivalently, in terms of the ``physical'' normalization often used in Christoffel--Darboux formula, \[ \widetilde{K}_N(1,t) = \sum_{n=0}^N (2n+1)\,P_n(1)\,P_n(t) = \sum_{n=0}^N (2n+1)\,P_n(t), \] and K_N differs from \( \displaystyle \quad \widetilde{K}_N \quad \) only by a f and K_N differs from \( \displaystyle \quad \widetilde{K}_N \quad \) only by a factor ½.

The Christoffel--Darboux formula now yields \begin{equation*} \widetilde{K}_N(1,t) = \frac{(N+1)\bigl(P_{N+1}(t)-P_N(t)\bigr)}{1-t}. \tag{T5.1} \end{equation*} Thus, \begin{equation*} S_N f(1) = \frac{1}{2} \int_{-1}^1 f(t)\, \widetilde{K}_N(1,t)\,dt. \tag{T5.3} \end{equation*}

Step 2: Asymptotics of the kernel near t = 1. Introduce t = cosθ, θ ∈ [0, π]. Then t → 1 corresponds to θ → 0. We have \[ 1-t = 1-\cos\theta = 2\sin^2\frac{\theta}{2} \sim \frac{\theta^2}{2}, \qquad \theta\to0. \] Moreover, for small θ and large n, the Bessel-type asymptotics for Legendre polynomials gives \begin{equation*} P_n(\cos\theta) = J_0\bigl((n+\tfrac12)\theta\bigr) + O\bigl(n^{-1}\bigr), \tag{T5.2} \end{equation*} uniformly for θ ∈ [0, θ₀] with fixed θ₀ > 0. Using the differential equation or standard recurrences, one derives from Eq.(T5.2) that \begin{equation*} P_{N+1}(\cos\theta)-P_N(\cos\theta) = O(1) \quad\text{for $0\le\theta\le c/N$}, \end{equation*} and \begin{equation*} P_{N+1}(\cos\theta)-P_N(\cos\theta) = O\biggl(\frac{1}{(N+1)\theta}\biggr) \quad\text{for $c/N\le\theta\le\theta_0$}, \end{equation*} with constants independent of N. Substituting into Eq.(T5.1), and recalling 1 − cosθ ∼ θ²/2, one obtains the Fejér-type scaling \begin{equation*} \bigl|\widetilde{K}_N\bigl(1,\cos\theta\bigr)\bigr| \le C\,\frac{N+1}{1+(N+1)^2\theta^2}, \qquad 0<\theta\le\pi, \tag{T5.4} \end{equation*} for some constant C > 0. This is the same type of estimate as in the Fejér kernel for trigonometric series, but now measured in the angular variable θ around θ = 0.

Step 3: Use of the weighted bounded variation. Assume \[ g(x) = \frac{f(x)}{\sqrt{1-x^2}} \in BV[-1,1]. \] Near x = 1, we have \[ \sqrt{1-x^2} = \sqrt{(1-x)(1+x)} \sim \sqrt{2(1-x)}. \] In terms of θ, \[ 1-x = 1-\cos\theta \sim \frac{\theta^2}{2}, \quad \sqrt{1-x^2} \sim c\,\theta, \qquad\theta\to0, \] for some constant c > 0. Hence, \[ f(\cos\theta) = g(\cos\theta)\,\sqrt{1-\cos^2\theta} \sim c\,\theta\, g(\cos\theta), \qquad \theta\to0. \] Since g is of bounded variation on [-1, 1], in particular it has finite one-sided limit at x = 1; denote g(1-0). Then \[ f(1-0) = \lim_{\theta\to0} f(\cos\theta) = \lim_{\theta\to0} c\,\theta\,g(\cos\theta) = 0 \quad\text{if }g(1-0)\text{ is finite}, \] unless g has a singularity that compensates the vanishing of \( \displaystyle \quad \sqrt{1-x^2} . \quad \) More generally, the bounded variation of $g$ ensures that f(cosθ) behaves ``almost like'' a constant times θ near θ = 0, in a controlled way.

Returning to Eq.(T5.3), and changing variables t = &cosθ, dt = −sinθ dθ, we write \[ S_N f(1) = \frac12 \int_0^\pi f(\cos\theta)\, \widetilde{K}_N(1,\cos\theta)\,\sin\theta\,d\theta. \] Set \[ h(\theta) := \frac{f(\cos\theta)}{\sin\theta}, \] so that $h(\theta)$ is of bounded variation on (0, π)$ by the hypothesis on g and the relation \( \displaystyle \quad \sqrt{1-\cos^2\theta} = \sin\theta . \quad \) Then \[ f(\cos\theta)\,\sin\theta = h(\theta)\,\sin^2\theta, \] and the integral can be written as \[ S_N f(1) = \frac12 \int_0^\pi h(\theta)\,\sin^2\theta\, \widetilde{K}_N(1,\cos\theta)\,d\theta. \] By Eq.(T4.4), the kernel \( \displaystyle \quad \sin^2\theta\,\widetilde{K}_N(1,\cos\theta) \quad \) has total mass uniformly bounded and concentrates near θ = 0 with the Fejér-type scaling (T4.4). The function h has bounded variation on [0, π] and a finite one-sided limit at θ = 0, namely \[ h(0+0) = \lim_{\theta\to0^+} \frac{f(\cos\theta)}{\sin\theta} = \lim_{x\to1^-} \frac{f(x)}{\sqrt{1-x^2}} = g(1-0), \] by the definition of g. Therefore, the standard approximate identity argument (as for Fejér kernels) implies that \[ \lim_{N\to\infty} S_N f(1) = f(1-0). \] More explicitly, one writes \[ S_N f(1) - f(1-0) = \frac12 \int_0^\pi [h(\theta)-h(0+0)]\,\sin^2\theta\, \widetilde{K}_N(1,\cos\theta)\,d\theta + \frac12\,h(0+0) \int_0^\pi \sin^2\theta\, \bigl(\widetilde{K}_N(1,\cos\theta)-c_N\bigr)\,d\theta, \] where c_N is chosen so that \( \displaystyle \quad \int_0^\pi \sin^2\theta\,\widetilde{K}_N(1,\cos\theta)\,{\text d}\theta = c_N . \quad \) The second integral vanishes by orthogonality, and the first term tends to zero by the localization properties of the kernel (T5.4) and the bounded variation of h (exactly as in the classical proof of convergence at a point of continuity for Cesàro means, but here we have a Dirichlet-type kernel with Fejér scaling).

Thus, the first limit relation holds. The proof of similar relation at x = -1 is analogous.

Theorem: If segment [𝑎, b] is finite and an auxiliary function \( \displaystyle \quad \varphi_x (t) = \frac{f(x) - f(t)}{x-t} \quad \) belongs to 𝔏²([𝑎, b]) for fixed x ∈ [𝑎, b], and the sequence of orthonormal polynomials {pₙ) is bounded at point x, then the Fourier series with respect to these orthonormal polynomials converges to f(x).

Similarly, we can estimate the second integral and obtain \[ A_1 (x) \leqslant \frac{c_4}{\sqrt{\varepsilon}\left( 1- x^2 \right)^{1/4}} \leqslant \frac{c_5}{\varepsilon^{3/4}} , \] where constant c₅ does not depend on both, ε and n. A similar estimation is valid for the fifth integral.

Next, we get an estimation for the integral over second interval: \[ A_2 (x) \leqslant \frac{c_6}{\sqrt{\varepsilon}} \int_{\Delta_2} \frac{{\text d}t}{x-t} \leqslant \frac{c_7}{\sqrt{\varepsilon}} \int_{1/n}^1 \frac{{\text d}\tau}{\tau} = \frac{c_7}{\sqrt{\varepsilon}} \, \ln n . \] A similar estimation holds for the integral over Δ₄.

Finally, we estimate the integral over the third interval [x − 1/n, x + 1/n]. We have \[ A_3 (x) = \int_{\Delta_3} \left\vert \sum_{k=0}^n \hat{P}_k (x)\, \hat{P}_k (t) \right\vert {\text d}t \leqslant \frac{c_8 \left( n+1 \right)}{\sqrt{\varepsilon}} \int_{\Delta_3} {\text d}t \leqslant \frac{c_9}{\sqrt{\varepsilon}} . \] Uniting all these estimates, we obtain \[ L_n (x) = \int_{-1}^1 \left\vert \sum_{k=0}^n \hat{P}_k (x)\, \hat{P}_k (t) \right\vert {\text d}t \leqslant \frac{c_{10} \ln n}{\varepsilon^{3/4}} , \quad x \in [-1+\varepsilon , 1-\varepsilon ] . \] This gives \[ \lim_{n\to\infty} E_n (f)\,\ln n = 0 . \]