[proofplan]
We differentiate the time-dependent integral defining the pairing, keeping track of both the time derivative of the integrand and the time derivative of the Riemannian volume measure. Under Ricci flow, the volume measure evolves by $\partial_t d\mu_{g(t)} = -R_{g(t)}d\mu_{g(t)}$, and this cancels the scalar curvature term in the conjugate [heat equation](/page/Heat%20Equation). After the cancellation, the derivative is the difference between the two Laplace-Beltrami pairing terms, and those terms cancel by self-adjointness of the Laplace-Beltrami operator on a compact manifold without boundary.
[/proofplan]
[step:Define the time-dependent pairing and record the volume variation]
Define the time-dependent pairing $F:I\to\mathbb{R}$ by
\begin{align*}
F(t) := \int_M u(x,t)v(x,t)\,d\mu_{g(t)}(x).
\end{align*}
Since $M$ is compact and $u$, $v$, and $g(t)$ are smooth in space and time, $F$ is differentiable at every interior point of $I$.
We use the following volume variation identity under Ricci flow:
\begin{align*}
\partial_t d\mu_{g(t)} = -R_{g(t)}\,d\mu_{g(t)}.
\end{align*}
Indeed, in a coordinate chart, write $g_{ij}(x,t)$ for the components of $g(t)$ and write
\begin{align*}
\theta(x,t) := \sqrt{\det(g_{ij}(x,t))}
\end{align*}
for the local density of $d\mu_{g(t)}$ with respect to coordinate [Lebesgue measure](/page/Lebesgue%20Measure). For a symmetric $(0,2)$-tensor $h$ on $M$, define its $g(t)$-trace by
\begin{align*}
\operatorname{tr}_{g(t)} h := \sum_{i,j} g(t)^{ij}h_{ij},
\end{align*}
where $g(t)^{ij}$ are the components of the inverse metric in the chosen coordinate chart. By the determinant differential identity $\partial_t\det A(t)=\det A(t)\operatorname{tr}(A(t)^{-1}\partial_t A(t))$, applied to the positive definite matrix $A(t)=(g_{ij}(x,t))$, and by the chain rule for $s\mapsto \sqrt{s}$ on $(0,\infty)$, we have
\begin{align*}
\partial_t \theta = \frac{1}{2}\operatorname{tr}_{g(t)}(\partial_t g(t))\,\theta.
\end{align*}
Since the Ricci flow equation is $\partial_t g(t)=-2\operatorname{Ric}_{g(t)}$, this becomes
\begin{align*}
\partial_t \theta = \frac{1}{2}\operatorname{tr}_{g(t)}(-2\operatorname{Ric}_{g(t)})\,\theta.
\end{align*}
By definition of scalar curvature, $R_{g(t)}=\operatorname{tr}_{g(t)}\operatorname{Ric}_{g(t)}$, so
\begin{align*}
\partial_t \theta = -R_{g(t)}\theta.
\end{align*}
This coordinate identity is invariant under chart changes, hence gives the stated measure variation formula.
[/step]
[step:Differentiate the pairing and substitute the two evolution equations]
To justify differentiating the integral, choose a finite smooth atlas on the compact manifold $M$ and a smooth [partition of unity](/page/Partition%20of%20Unity) subordinate to it. In each chart, the localized integral is an integral over a fixed relatively compact coordinate domain with smooth density depending smoothly on $t$; compactness gives a uniform compact support, so the standard smooth-parameter differentiation theorem for Lebesgue integrals applies. Summing the finitely many localized identities and using the product rule for the smooth integrand together with the volume variation formula from the previous step, for each interior point $t \in I$ we obtain
\begin{align*}
F'(t)=\int_M \partial_t(uv)(x,t)\,d\mu_{g(t)}(x)+\int_M u(x,t)v(x,t)\,\partial_t d\mu_{g(t)}(x).
\end{align*}
Using $\partial_t(uv)=(\partial_t u)v+u(\partial_t v)$ and $\partial_t d\mu_{g(t)}=-R_{g(t)}d\mu_{g(t)}$, this gives
\begin{align*}
F'(t)=\int_M \left((\partial_t u)(x,t)v(x,t)+u(x,t)(\partial_t v)(x,t)\right)\,d\mu_{g(t)}(x)-\int_M R_{g(t)}(x)u(x,t)v(x,t)\,d\mu_{g(t)}(x).
\end{align*}
Substituting $\partial_t u=\Delta_{g(t)}u$ and $\partial_t v=-\Delta_{g(t)}v+R_{g(t)}v$ gives
\begin{align*}
F'(t)=\int_M v(x,t)\Delta_{g(t)}u(x,t)\,d\mu_{g(t)}(x)+\int_M u(x,t)\left(-\Delta_{g(t)}v(x,t)+R_{g(t)}(x)v(x,t)\right)\,d\mu_{g(t)}(x)-\int_M R_{g(t)}(x)u(x,t)v(x,t)\,d\mu_{g(t)}(x).
\end{align*}
The scalar-curvature integral from the conjugate heat equation cancels the scalar-curvature integral from the measure variation, so
\begin{align*}
F'(t)=\int_M v(x,t)\Delta_{g(t)}u(x,t)\,d\mu_{g(t)}(x)-\int_M u(x,t)\Delta_{g(t)}v(x,t)\,d\mu_{g(t)}(x).
\end{align*}
[guided]
We want to compute the derivative of
\begin{align*}
F(t)=\int_M u(x,t)v(x,t)\,d\mu_{g(t)}(x),
\end{align*}
but there are two sources of time dependence: the functions $u$ and $v$, and the measure $d\mu_{g(t)}$. Since $M$ is compact and all objects are smooth, differentiating under the integral sign is valid. The product rule gives
\begin{align*}
\partial_t(uv)(x,t)=(\partial_t u)(x,t)v(x,t)+u(x,t)(\partial_t v)(x,t).
\end{align*}
We also need the time derivative of the Riemannian measure. Let $(U_{\mathrm{chart}},\varphi)$ be a coordinate chart on $M$, where $U_{\mathrm{chart}}\subset M$ is the chart domain and $\varphi:U_{\mathrm{chart}}\to\varphi(U_{\mathrm{chart}})\subset\mathbb{R}^n$ is the coordinate map. In this chart, let $g_{ij}(x,t)$ be the components of $g(t)$ and define the local density $\theta: \varphi(U_{\mathrm{chart}})\times I\to (0,\infty)$ by
\begin{align*}
\theta(x,t):=\sqrt{\det(g_{ij}(x,t))},
\end{align*}
where $x\in\varphi(U_{\mathrm{chart}})$ is the coordinate variable. Let $\mathcal{L}^n$ denote $n$-dimensional Lebesgue measure on $\mathbb{R}^n$. Locally, the Riemannian measure is represented by $d\mu_{g(t)}=\theta(x,t)\,d\mathcal{L}^n(x)$. The determinant differential identity gives
\begin{align*}
\partial_t\theta=\frac{1}{2}\operatorname{tr}_{g(t)}(\partial_t g(t))\theta.
\end{align*}
Since $g(t)$ solves Ricci flow, $\partial_tg(t)=-2\operatorname{Ric}_{g(t)}$. Taking the $g(t)$-trace and using $R_{g(t)}=\operatorname{tr}_{g(t)}\operatorname{Ric}_{g(t)}$ gives
\begin{align*}
\partial_t\theta=-R_{g(t)}\theta.
\end{align*}
This local density identity is invariant under coordinate changes, hence
\begin{align*}
\partial_t d\mu_{g(t)} = -R_{g(t)}\,d\mu_{g(t)}.
\end{align*}
Therefore
\begin{align*}
F'(t)=\int_M \partial_t(uv)(x,t)\,d\mu_{g(t)}(x)+\int_M u(x,t)v(x,t)\,\partial_t d\mu_{g(t)}(x).
\end{align*}
Using the product rule and the volume variation identity, this becomes
\begin{align*}
F'(t)=\int_M \left((\partial_t u)(x,t)v(x,t)+u(x,t)(\partial_t v)(x,t)\right)\,d\mu_{g(t)}(x)-\int_M R_{g(t)}(x)u(x,t)v(x,t)\,d\mu_{g(t)}(x).
\end{align*}
Now we use the two PDEs. The forward heat equation gives
\begin{align*}
\partial_t u = \Delta_{g(t)}u,
\end{align*}
and the conjugate heat equation gives
\begin{align*}
\partial_t v = -\Delta_{g(t)}v+R_{g(t)}v.
\end{align*}
Substituting these identities into the differentiated pairing yields
\begin{align*}
F'(t)=\int_M v(x,t)\Delta_{g(t)}u(x,t)\,d\mu_{g(t)}(x)+\int_M u(x,t)\left(-\Delta_{g(t)}v(x,t)+R_{g(t)}(x)v(x,t)\right)\,d\mu_{g(t)}(x)-\int_M R_{g(t)}(x)u(x,t)v(x,t)\,d\mu_{g(t)}(x).
\end{align*}
The scalar curvature term from the conjugate heat equation is
\begin{align*}
\int_M R_{g(t)}(x)u(x,t)v(x,t)\,d\mu_{g(t)}(x),
\end{align*}
and the scalar curvature term from the evolving measure is its negative. These two terms cancel exactly. Hence
\begin{align*}
F'(t)
= \int_M v(x,t)\Delta_{g(t)}u(x,t)\,d\mu_{g(t)}(x)
- \int_M u(x,t)\Delta_{g(t)}v(x,t)\,d\mu_{g(t)}(x).
\end{align*}
This is the point of the conjugate heat equation: the additional $R_{g(t)}v$ term compensates precisely for the loss of fixed volume under Ricci flow.
[/guided]
[/step]
[step:Use compactness without boundary to make the Laplacian terms cancel]
Fix an interior point $t \in I$. For this fixed time, define the smooth functions $u_t:M\to\mathbb{R}$ and $v_t:M\to\mathbb{R}$ by $u_t(x):=u(x,t)$ and $v_t(x):=v(x,t)$.
Let $\nabla^{g(t)} f$ denote the Riemannian gradient of a smooth function $f: M \to \mathbb{R}$ with respect to the metric $g(t)$, characterized by $g(t)(\nabla^{g(t)}f,X)=Xf$ for every smooth vector field $X$ on $M$. Since $M$ is compact and has no boundary, the Laplace-Beltrami operator is symmetric with respect to $d\mu_{g(t)}$:
\begin{align*}
\int_M v_t(x)\Delta_{g(t)}u_t(x)\,d\mu_{g(t)}(x)
=
\int_M u_t(x)\Delta_{g(t)}v_t(x)\,d\mu_{g(t)}(x).
\end{align*}
This follows from [integration by parts](/theorems/2098) on a compact manifold without boundary. Applying the divergence identity to the smooth vector fields $v_t\nabla^{g(t)}u_t$ and $u_t\nabla^{g(t)}v_t$ gives
\begin{align*}
\int_M v_t(x)\Delta_{g(t)}u_t(x)\,d\mu_{g(t)}(x)=-\int_M \langle \nabla^{g(t)}v_t(x),\nabla^{g(t)}u_t(x)\rangle_{g(t)}\,d\mu_{g(t)}(x).
\end{align*}
and
\begin{align*}
\int_M u_t(x)\Delta_{g(t)}v_t(x)\,d\mu_{g(t)}(x)=-\int_M \langle \nabla^{g(t)}u_t(x),\nabla^{g(t)}v_t(x)\rangle_{g(t)}\,d\mu_{g(t)}(x).
\end{align*}
The two right-hand sides are equal by symmetry of the Riemannian metric $g(t)$.
[guided]
Fix an interior point $t\in I$ and define $u_t:M\to\mathbb{R}$ and $v_t:M\to\mathbb{R}$ by $u_t(x):=u(x,t)$ and $v_t(x):=v(x,t)$. At this fixed time, the metric $g(t)$ is a smooth Riemannian metric on the compact manifold $M$, and $M$ has no boundary. Therefore [integration by parts](/theorems/210) for the Laplace-Beltrami operator applies without a boundary term.
Let $\nabla^{g(t)}f$ denote the Riemannian gradient of a smooth function $f:M\to\mathbb{R}$ with respect to $g(t)$, characterized by $g(t)(\nabla^{g(t)}f,X)=Xf$ for every smooth vector field $X$ on $M$. Applying the [divergence theorem](/theorems/2754) to the smooth vector field $v_t\nabla^{g(t)}u_t$ gives
\begin{align*}
\int_M v_t(x)\Delta_{g(t)}u_t(x)\,d\mu_{g(t)}(x)=-\int_M \langle \nabla^{g(t)}v_t(x),\nabla^{g(t)}u_t(x)\rangle_{g(t)}\,d\mu_{g(t)}(x).
\end{align*}
Applying the same argument to $u_t\nabla^{g(t)}v_t$ gives
\begin{align*}
\int_M u_t(x)\Delta_{g(t)}v_t(x)\,d\mu_{g(t)}(x)=-\int_M \langle \nabla^{g(t)}u_t(x),\nabla^{g(t)}v_t(x)\rangle_{g(t)}\,d\mu_{g(t)}(x).
\end{align*}
The two integrands on the right are equal because the Riemannian metric is symmetric. Hence
\begin{align*}
\int_M v_t(x)\Delta_{g(t)}u_t(x)\,d\mu_{g(t)}(x)=\int_M u_t(x)\Delta_{g(t)}v_t(x)\,d\mu_{g(t)}(x).
\end{align*}
This is the cancellation needed after the scalar-curvature terms have already disappeared.
[/guided]
[/step]
[step:Conclude that the pairing is constant]
Combining the differentiated formula with the symmetry identity gives, for every interior point $t \in \operatorname{int}(I)$,
\begin{align*}
F'(t)=0.
\end{align*}
Since $M$ is compact and $u$, $v$, and $g(t)$ are smooth, $F$ is continuous on $I$. If $a,b\in \operatorname{int}(I)$ with $a<b$, then the ordinary [mean value theorem](/theorems/186) applied to $F:[a,b]\to\mathbb{R}$ gives a point $c\in(a,b)$ such that
\begin{align*}
F(b)-F(a)=F'(c)(b-a)=0.
\end{align*}
Thus $F$ is constant on $\operatorname{int}(I)$. For endpoints of $I$, if present, the same value is obtained by taking limits from the interior and using continuity of $F$ on $I$. Hence $F$ is constant on $I$, and equivalently
\begin{align*}
\frac{d}{dt}\int_M u(x,t)v(x,t)\,d\mu_{g(t)}(x)=0
\end{align*}
for every interior time $t\in\operatorname{int}(I)$.
[/step]