[proofplan]
We compute the first variation of Perelman's $\mathcal L$-length on an arbitrary compact time interval $[a,b]\subset(0,\tau_2]$ with fixed endpoints. The scalar curvature term contributes the spatial gradient of $R$, while the kinetic term is integrated by parts along the time-dependent metric $g(\tau)$. The time dependence of the metric, through $\partial_\tau g=2\operatorname{Ric}$, produces the Ricci term, and differentiating the weight $\sqrt{\tau}$ produces the term $\frac{1}{2\tau}\dot{\gamma}$. Since the resulting integral vanishes against every compactly supported variation field, the coefficient vector field must vanish pointwise.
[/proofplan]
[step:Compute the first variation on a compact time interval]
Fix $0<a<b\leq \tau_2$. Let $\varepsilon>0$, and let $\Gamma:(-\varepsilon,\varepsilon)\times [a,b]\to M$ be a smooth variation of $\gamma|_{[a,b]}$, written $\Gamma(s,\tau)=\gamma_s(\tau)$, with fixed endpoints, so $\gamma_0=\gamma|_{[a,b]}$, $\gamma_s(a)=\gamma(a)$, and $\gamma_s(b)=\gamma(b)$ for all $s$. Define the variation field $X:[a,b]\to TM$ by
\begin{align*}
X(\tau)=\left.\frac{\partial \gamma_s(\tau)}{\partial s}\right|_{s=0}\in T_{\gamma(\tau)}M.
\end{align*}
Then $X(a)=X(b)=0$.
For each $\tau\in[a,b]$, all gradients, inner products, norms, Ricci tensors, and connections below are taken with respect to $g(\tau)$. Differentiating the scalar curvature term in the spatial direction $X(\tau)$ gives
\begin{align*}
\left.\frac{\partial}{\partial s}\right|_{s=0}
R(\gamma_s(\tau),\tau)
=
g(\nabla R(\gamma(\tau),\tau),X(\tau)).
\end{align*}
For the kinetic term, fix $\tau\in[a,b]$ and use the Levi-Civita connection of $g(\tau)$. This connection is torsion-free, so along the two-parameter map $\Gamma$ one has $\nabla_{\partial_s}\partial_\tau\Gamma=\nabla_{\partial_\tau}\partial_s\Gamma$ at time $\tau$; it is also metric-compatible. Hence
\begin{align*}
\left.\frac{\partial}{\partial s}\right|_{s=0}|\dot{\gamma}_s(\tau)|_{g(\tau)}^2=2g(\nabla_{\dot{\gamma}}X,\dot{\gamma}).
\end{align*}
Therefore criticality of $\gamma$ gives
\begin{align*}
0
=
\int_a^b \sqrt{\tau}\left(
g(\nabla R,X)
+
2g(\nabla_{\dot{\gamma}}X,\dot{\gamma})
\right)\,d\mathcal L^1(\tau).
\end{align*}
[guided]
Fix a compact interval $[a,b]\subset(0,\tau_2]$ so that the factor $\sqrt{\tau}$ is smooth and positive on the whole interval. We vary $\gamma$ through curves with the same endpoints. Thus we take $\varepsilon>0$ and a smooth map $\Gamma:(-\varepsilon,\varepsilon)\times [a,b]\to M$, written $\Gamma(s,\tau)=\gamma_s(\tau)$, with $\gamma_0=\gamma|_{[a,b]}$, $\gamma_s(a)=\gamma(a)$, and $\gamma_s(b)=\gamma(b)$. The associated variation field is the vector field $X:[a,b]\to TM$ along $\gamma$ defined by
\begin{align*}
X(\tau)=\left.\frac{\partial \gamma_s(\tau)}{\partial s}\right|_{s=0}\in T_{\gamma(\tau)}M.
\end{align*}
The fixed endpoint condition implies $X(a)=X(b)=0$.
Now differentiate the two pieces of the integrand. First, for each fixed $\tau$, the scalar curvature is a smooth function
\begin{align*}
R(\cdot,\tau):M \to \mathbb{R}.
\end{align*}
The directional derivative of this function in the direction $X(\tau)$ is its metric gradient paired with $X(\tau)$:
\begin{align*}
\left.\frac{\partial}{\partial s}\right|_{s=0}
R(\gamma_s(\tau),\tau)
=
g(\nabla R(\gamma(\tau),\tau),X(\tau)).
\end{align*}
Second, consider the kinetic energy term. For a fixed value of $\tau$, the metric is $g(\tau)$. Its Levi-Civita connection is metric-compatible and torsion-free. Torsion-freeness gives $\nabla_{\partial_s}\partial_\tau\Gamma=\nabla_{\partial_\tau}\partial_s\Gamma$ at this fixed time, so differentiating the velocity field in the variation direction is the same as covariantly differentiating the variation field along $\gamma$. Therefore
\begin{align*}
\left.\frac{\partial}{\partial s}\right|_{s=0}|\dot{\gamma}_s(\tau)|_{g(\tau)}^2=2g(\nabla_{\dot{\gamma}}X,\dot{\gamma}).
\end{align*}
The weight $\sqrt{\tau}$ does not depend on $s$, so differentiating the $\mathcal L$-length under the integral sign yields
\begin{align*}
\left.\frac{d}{ds}\right|_{s=0}
\mathcal L[\gamma_s]
=
\int_a^b \sqrt{\tau}\left(
g(\nabla R,X)
+
2g(\nabla_{\dot{\gamma}}X,\dot{\gamma})
\right)\,d\mathcal L^1(\tau).
\end{align*}
Because $\gamma$ is an $\mathcal L$-geodesic, this first variation vanishes for every fixed-endpoint variation, so
\begin{align*}
0
=
\int_a^b \sqrt{\tau}\left(
g(\nabla R,X)
+
2g(\nabla_{\dot{\gamma}}X,\dot{\gamma})
\right)\,d\mathcal L^1(\tau).
\end{align*}
[/guided]
[/step]
[step:Integrate the kinetic variation by parts using the backward Ricci flow]
For the vector fields $X$ and $\dot{\gamma}$ along $\gamma$, the product rule for the time-dependent metric gives
\begin{align*}
\frac{d}{d\tau}g(X,\dot{\gamma})
=
\frac{\partial g}{\partial \tau}(X,\dot{\gamma})
+
g(\nabla_{\dot{\gamma}}X,\dot{\gamma})
+
g(X,\nabla_{\dot{\gamma}}\dot{\gamma}).
\end{align*}
Using $\frac{\partial g}{\partial \tau}=2\operatorname{Ric}$, we obtain
\begin{align*}
g(\nabla_{\dot{\gamma}}X,\dot{\gamma})
=
\frac{d}{d\tau}g(X,\dot{\gamma})
-
2\operatorname{Ric}(X,\dot{\gamma})
-
g(X,\nabla_{\dot{\gamma}}\dot{\gamma}).
\end{align*}
Hence the integral of the kinetic variation decomposes as
\begin{align*}
\int_a^b 2\sqrt{\tau}\,g(\nabla_{\dot{\gamma}}X,\dot{\gamma})\,d\mathcal L^1(\tau)=\int_a^b \left(2\sqrt{\tau}\,\frac{d}{d\tau}g(X,\dot{\gamma})-4\sqrt{\tau}\,\operatorname{Ric}(X,\dot{\gamma})-2\sqrt{\tau}\,g(X,\nabla_{\dot{\gamma}}\dot{\gamma})\right)\,d\mathcal L^1(\tau).
\end{align*}
Integrating the total derivative term by parts and using $X(a)=X(b)=0$ gives
\begin{align*}
\int_a^b 2\sqrt{\tau}\,\frac{d}{d\tau}g(X,\dot{\gamma})\,d\mathcal L^1(\tau)=-\int_a^b \frac{1}{\sqrt{\tau}}g(X,\dot{\gamma})\,d\mathcal L^1(\tau).
\end{align*}
Substituting into the first variation identity yields
\begin{align*}
0
=
\int_a^b
g\left(
\sqrt{\tau}\nabla R
-
\frac{1}{\sqrt{\tau}}\dot{\gamma}
-
4\sqrt{\tau}\operatorname{Ric}(\dot{\gamma},\cdot)^\sharp
-
2\sqrt{\tau}\nabla_{\dot{\gamma}}\dot{\gamma},
X
\right)\,d\mathcal L^1(\tau).
\end{align*}
[guided]
The first variation still contains the derivative of the variation field, namely the term $g(\nabla_{\dot{\gamma}}X,\dot{\gamma})$. To obtain an Euler-Lagrange equation, we must move this derivative off $X$ and onto the coefficient of $X$ by [integration by parts](/theorems/2098).
For the vector fields $X$ and $\dot{\gamma}$ along $\gamma$, the metric depends on $\tau$, so the ordinary product rule has an extra term. The time-dependent metric product rule gives
\begin{align*}
\frac{d}{d\tau}g(X,\dot{\gamma})
=
\frac{\partial g}{\partial \tau}(X,\dot{\gamma})
+
g(\nabla_{\dot{\gamma}}X,\dot{\gamma})
+
g(X,\nabla_{\dot{\gamma}}\dot{\gamma}).
\end{align*}
Because the metric evolves by backward Ricci flow, $\frac{\partial g}{\partial \tau}=2\operatorname{Ric}$. Solving this identity for the derivative term gives
\begin{align*}
g(\nabla_{\dot{\gamma}}X,\dot{\gamma})
=
\frac{d}{d\tau}g(X,\dot{\gamma})
-
2\operatorname{Ric}(X,\dot{\gamma})
-
g(X,\nabla_{\dot{\gamma}}\dot{\gamma}).
\end{align*}
Substituting this into the kinetic part of the first variation yields
\begin{align*}
\int_a^b 2\sqrt{\tau}\,g(\nabla_{\dot{\gamma}}X,\dot{\gamma})\,d\mathcal L^1(\tau)=\int_a^b \left(2\sqrt{\tau}\,\frac{d}{d\tau}g(X,\dot{\gamma})-4\sqrt{\tau}\,\operatorname{Ric}(X,\dot{\gamma})-2\sqrt{\tau}\,g(X,\nabla_{\dot{\gamma}}\dot{\gamma})\right)\,d\mathcal L^1(\tau).
\end{align*}
Now apply the one-dimensional [integration by parts](/theorems/210) formula on $[a,b]$ to the scalar function $\tau \mapsto g(X(\tau),\dot{\gamma}(\tau))$ and the weight $\tau \mapsto 2\sqrt{\tau}$. Since $X(a)=X(b)=0$, the boundary term vanishes. Since $\frac{d}{d\tau}(2\sqrt{\tau})=\frac{1}{\sqrt{\tau}}$, we get
\begin{align*}
\int_a^b 2\sqrt{\tau}\,\frac{d}{d\tau}g(X,\dot{\gamma})\,d\mathcal L^1(\tau)=-\int_a^b \frac{1}{\sqrt{\tau}}g(X,\dot{\gamma})\,d\mathcal L^1(\tau).
\end{align*}
Combining this with the scalar curvature term from the first variation gives
\begin{align*}
0
=
\int_a^b
g\left(
\sqrt{\tau}\nabla R
-
\frac{1}{\sqrt{\tau}}\dot{\gamma}
-
4\sqrt{\tau}\operatorname{Ric}(\dot{\gamma},\cdot)^\sharp
-
2\sqrt{\tau}\nabla_{\dot{\gamma}}\dot{\gamma},
X
\right)\,d\mathcal L^1(\tau).
\end{align*}
This is the key reduction: the first variation is now an integral pairing an arbitrary variation field $X$ against a coefficient vector field along $\gamma$.
[/guided]
[/step]
[step:Use arbitrary compactly supported variations to identify the Euler-Lagrange field]
Define the vector field $E:[a,b]\to TM$ along $\gamma|_{[a,b]}$ by
\begin{align*}
E(\tau)=\sqrt{\tau}\nabla R(\gamma(\tau),\tau)-\frac{1}{\sqrt{\tau}}\dot{\gamma}(\tau)-4\sqrt{\tau}\operatorname{Ric}(\dot{\gamma}(\tau),\cdot)^{\sharp}-2\sqrt{\tau}\nabla_{\dot{\gamma}}\dot{\gamma}(\tau).
\end{align*}
The preceding step says that
\begin{align*}
\int_a^b g(E(\tau),X(\tau))\,d\mathcal L^1(\tau)=0
\end{align*}
for every smooth variation field $X$ along $\gamma$ with $X(a)=X(b)=0$. Conversely, every smooth compactly supported vector field $X:(a,b)\to TM$ along $\gamma$ arises from a fixed-endpoint variation. Indeed, let $K=\operatorname{supp}X\subset(a,b)$. Since $K$ is compact and $\gamma(K)$ is compact in $M$, choose finitely many coordinate charts whose domains cover $\gamma(K)$, and choose a smooth [partition of unity](/page/Partition%20of%20Unity) in the parameter $\tau$ subordinate to the corresponding open cover of $K$. On each coordinate subinterval, define the local variation by adding $sX(\tau)$ to the coordinate representative of $\gamma(\tau)$; after possibly shrinking $|s|<\varepsilon$, compactness of $K$ ensures that all perturbed coordinate representatives remain in the chart domains. Summing the local coordinate perturbations with the partition of unity gives a smooth variation of $\gamma$ on $[a,b]$ whose variation field is $X$. Because $X$ has compact support in $(a,b)$, this variation equals $\gamma$ near $a$ and $b$, so it has fixed endpoints. Therefore the identity holds for every smooth compactly supported vector field along $\gamma$ on $(a,b)$.
We claim that $E(\tau)=0$ for every $\tau\in(a,b)$. If not, choose $\tau_0\in(a,b)$ with $E(\tau_0)\neq 0$. By smoothness and positive-definiteness of $g(\tau_0)$, there is an open interval $I\subset(a,b)$ containing $\tau_0$ and a smooth vector field
\begin{align*}
Y:I\to TM, \qquad Y(\tau)\in T_{\gamma(\tau)}M,
\end{align*}
along $\gamma|_I$ such that $g(E,Y)>0$ on $I$. Let
\begin{align*}
\eta:I\to[0,\infty)
\end{align*}
be a smooth function with compact support in $I$ and not identically zero, and extend $\eta Y$ by zero to a smooth compactly supported vector field $X$ along $\gamma$ on $(a,b)$. Then
\begin{align*}
\int_a^b g(E,X)\,d\mathcal L^1(\tau)
=
\int_I \eta(\tau)g(E(\tau),Y(\tau))\,d\mathcal L^1(\tau)
>
0,
\end{align*}
contradicting the first variation identity. Therefore $E=0$ on $(a,b)$.
[/step]
[step:Divide by the positive weight and obtain the geodesic equation]
Since $\tau>0$ on $(a,b)$, the identity $E(\tau)=0$ gives
\begin{align*}
0=\sqrt{\tau}\nabla R-\frac{1}{\sqrt{\tau}}\dot{\gamma}-4\sqrt{\tau}\operatorname{Ric}(\dot{\gamma},\cdot)^\sharp-2\sqrt{\tau}\nabla_{\dot{\gamma}}\dot{\gamma}.
\end{align*}
Multiplying by $-\frac{1}{2\sqrt{\tau}}$ gives
\begin{align*}
\nabla_{\dot{\gamma}}\dot{\gamma}
-
\frac{1}{2}\nabla R
+
2\operatorname{Ric}(\dot{\gamma},\cdot)^\sharp
+
\frac{1}{2\tau}\dot{\gamma}
=
0
\end{align*}
on $(a,b)$. Because $0<a<b<\tau_2$ were arbitrary, the same identity holds for every $\tau\in(0,\tau_2)$. The curve $\gamma$, the backward Ricci flow metric $g(\tau)$, the scalar curvature $R(\cdot,\tau)$, and the Ricci tensor are smooth up to $\tau=\tau_2$ by the hypotheses implicit in the definition of an $\mathcal L$-geodesic for the flow. Hence the Euler-Lagrange expression is continuous as $\tau\uparrow\tau_2$, and the identity also holds at $\tau=\tau_2$. This is the desired Euler-Lagrange equation for $\mathcal L$-geodesics, with all geometric quantities evaluated using $g(\tau)$.
[/step]