\section{Detailed Derivations}
\label{sec:detailed_derivations}
%
\subsection{Deriving $\sdev\kernelT\,\kernelC\,\Lv$~\eqref{eqn:rte3d_ind_diff}}
\label{ssec:kernelT_deriv}
%
According to Eqs.~\eqref{eqn:rte3d_kernelT} and \eqref{eqn:rte3d_Lins},
%
\begin{equation}
	\label{eqn:rte3d_ind}
	(\kernelT\,\kernelC\,\Lv)(\bx, \bom)
	= \tint_0^{\dist} \trans(\bx', \bx) \, \sigS(\bx') \, \Lins(\bx', \bom) \,\D\tau,
\end{equation}
%
where $\bx' := \bx - \tau\bom$.
The assumptions stated at the beginning of \S\ref{sec:diff_rte} ensure that the integrand $\trans\,\sigS\,\Lins$ is continuous in $\tau$ for $0 < \tau < \dist$.
Thus, applying \thmref{reynolds}, which reduces to Leibniz's rule for differentiation in this case, to Eq.~\eqref{eqn:rte3d_ind} gives
%
\begin{equation}
	\label{eqn:rte3d_ind_diff_0}
	\begin{split}
		\!\!\!\!\left(\sdev\kernelT\,\kernelC\,\Lv\right)(\bx, \bom)
		=\;& \tint_0^{\dist} \sdev\left[ \trans(\bx', \bx) \, \sigS(\bx') \, \Lins(\bx', \bom) \right] \D\tau\\
		& + \ddist \, \trans(\bx_0, \bx) \, \sigS(\bx_0) \,\Lins(\bx_0, \bom),
	\end{split}
\end{equation}
%
where $\bx_0 := \bx - \dist\bom$.
Further, the derivative of the transmittance $\trans(\bx', \bx)$, which we derive in Appendix~\ref{ssec:trans_deriv}, equals
%
\begin{equation}
	\label{eqn:rte3d_trans0_diff}
	\dot{\trans}(\bx', \bx)	= -\trans(\bx', \bx) \,\intSigT(\bx, \bom, \tau).
\end{equation}
%
Lastly, combining Eqs.~\eqref{eqn:rte3d_ind_diff_0} and \eqref{eqn:rte3d_trans0_diff} yields Eq.~\eqref{eqn:rte3d_ind_diff} in \S\ref{ssec:dev_trans}.


\subsection{Differentiating Transmittance}
\label{ssec:trans_deriv}
%
Multiple terms we derived in \S\ref{sec:diff_rte} involve the derivative of the transmittance between two points $\bx$ and $\bx - t\bom$ in the medium given by
%
\begin{equation}
	\trans(\bx - t\bom, \bx) = \exp\left( -\tint_0^t \sigT(\bx - \tau\bom) \,\D\tau \right).
\end{equation}
%
Applying \thmref{reynolds} to this equation leads to
%
\begin{equation}
	\label{eqn:rte3d_trans_diff}
	\begin{split}
		\dot{\trans}(\bx - t\bom, \bx)
		&= \sdev\exp\left( -\tint_0^t \sigT(\bx - \tau\bom) \,\D\tau \right)\\
		&= -T(\bx - t\bom, \bx) \, \sdev\tint_0^t \sigT(\bx - \tau\bom) \,\D\tau\\
		&= -\trans(\bx - t\bom, \bx) \left( \intSigT(\bx, \bom, t) + \dot{t} \,\sigT(\bx - t\bom) \right),
	\end{split}
\end{equation}
%
where $\intSigT$ is defined in Eq.~\eqref{eqn:rte3d_int_sigT} as
$\intSigT(\bx, \bom, t) = \tint_0^t \dot{\sigT}(\bx - \tau\bom) \,\D\tau$.


\subsection{Deriving the Change Rate of \texorpdfstring{$\bom'$}{$\omega'$}~\eqref{eqn:rte3d_sph_vel}}
\label{ssec:dbom_derive}
%
Let $\bx$ and $\by$ be two 3D locations.
Then, the change rate of $\bom' := (\by \to \bx)$ with respect to a scene parameter $\param$ can be calculated as follows (see Figure~\ref{fig:proj_velocity}): (i) project $(\dot{\by} - \dbx)$ to the surface of a unit sphere centered at $\bx$, yielding $(\dot{\by} - \dbx)/\| \by - \bx \|$; (ii)~invert the projected velocity to align it with $\bom'$, yielding $(\dbx - \dot{\by})/\| \by - \bx \|$; (iii)~take the component (illustrated in red) of the inverted velocity that is within the tangent plane at $\bom'$, yielding $\dbom'$ expressed in Eq.~\eqref{eqn:rte3d_sph_vel}.

\begin{figure}[t]
	\centering
	\includegraphics[width=0.7\columnwidth]{proj_velocity.pdf}
	\caption{\label{fig:proj_velocity}
		Calculating the change rate of $\bom'$ defined as the unit vector pointing from $\by$ to $\bx$.
	}
\end{figure}

\subsection{Evaluating $\Delta\Lv$ Defined in Eq.~\eqref{eqn:rte3d_delta_L_def}}
\label{ssec:deltaL_deriv}
%
Given a point $\bx$ interior of the medium, the discontinuities of $\Lv(\bx, \bom)$ in $\bom$, in the absence of zero-measure light sources and ideal specular surfaces, generally arise when $\bom = (\by \to \bx)$ with $\by$ located on one of the boundary, silhouette, and sharp edges (as illustrated in Figure~\ref{fig:edge_types}).
In what follows, we derive $\Delta\Lv$ based on the type of edge on which $\by$ resides.

\begin{figure}[b]
	\centering
	\includegraphics[height=1in]{delta_L.pdf}
	\caption{\label{fig:deltaL_silhouette}
		Calculating $\Delta\Lv$ resulting from a silhouette edge.
		In this 2D example, the orange triangle is translating upward with a constant velocity relative to $\bx$ (while the yellow box stays stationary), and its top vertex creates a geometric discontinuity when viewed from $\bx$.
		Let $\bom'$ denote the direction corresponding to this discontinuity (i.e., the ray $(\bx, -\bom')$ intersects the orange triangle at its top vertex).
		Based on these observations, $\dist^-$, $\dist^+$, $\bx^-$, and $\bx^+$ defined in Eq.~\eqref{eqn:rte3d_d_lims} are illustrated.
	}
\end{figure}

\paragraph{Case 1}
\emph{Boundary and silhouette edges} can cause the following discontinuities with respect to $\bom$:
%
\begin{itemize}
	\item The distance $\dist$ for the ray $(\bx, -\bom)$ to hit the medium boundary, which is involved in the transport operator~$\kernelT$~\eqref{eqn:rte3d_kernelT}.
	\item The location $\bx_0 = \bx - \dist\bom$ where the ray intersects the boundary, which is needed by the source term~$\src$~\eqref{eqn:rte3d_src}.
\end{itemize}
%
For each $\bom$ where $\dist$ and $\bx_0$ are discontinuous, let 
%
\begin{equation}
	\label{eqn:rte3d_d_lims}
	\dist^- := \lim_{\epsilon \to 0^-} \dist(\bx, \bom + \epsilon\,\gnorm),
	\quad
	\dist^+ := \lim_{\epsilon \to 0^+} \dist(\bx, \bom + \epsilon\,\gnorm),
\end{equation}
%
where $\bx^- := \bx - \dist^-\bom$ and $\bx^+ := \bx - \dist^+\bom$ denote the one-sided limits of $\dist$ and $\bx_0$ (see Figure~\ref{fig:deltaL_silhouette}).
%
Then, $\Delta\Lv(\bx, \bom)$ is given by the difference of the interior radiance~$\Lv$ evaluated with $\dist = \dist^-,\; \bx_0 = \bx^-$ and $\dist = \dist^+,\; \bx_0 = \bx^+$.
Precisely,
%
\begin{equation}
	\label{eqn:rte3d_deltaL_silhouette}
	\begin{split}	
		& \Delta\Lv(\bx, \bom)\\
		\!\!\!\!=\;& \left[\tint_0^{\dist^-} \trans(\bx', \bx) \,(\kernelC\,\Lv)(\bx', \bom) \,\D\tau + \trans(\bx^-, \bx) \,\Ls(\bx^-, \bom)\right] \;-\\
		& \left[\tint_0^{\dist^+} \trans(\bx', \bx) \,(\kernelC\,\Lv)(\bx', \bom) \,\D\tau + \trans(\bx^+, \bx) \,\Ls(\bx^-, \bom)\right]\\
		\!\!\!\!=\;& \tint_{\dist^+}^{\dist^-} \trans(\bx', \bx) \,\sigS(\bx') \, \Lins(\bx', \bom) \,\D\tau\\
		& + \trans(\bx^-, \bx) \,\Ls(\bx^-, \bom) - \trans(\bx^+, \bx) \,\Ls(\bx^+ , \bom).
	\end{split}
\end{equation}
%
Notice that it is possible for $\dist^+ > \dist^-$.
In this case, the integral from $\dist^+$ to $\dist^-$ in Eq.~\eqref{eqn:rte3d_deltaL_silhouette} takes a negative value.

\paragraph{Case 2}
Another type of $\Lv$'s discontinuity in $\bom$ occurs when $\bx_0$ lies on a \emph{sharp edge}.
In this case, the surface normal $\gnorm$ at $\bx_0$ has jump discontinuities (despite $\dist$ and $\bx_0$ being continuous themselves). 
This makes $\Delta\Lv$ to be entirely determined by the difference of the source term $\src$ evaluated with the two different surface normals, denoted as $\gnorm^-$ and $\gnorm^+$, from the two faces sharing the sharp edge.
With the volumetric emission~$\Lve$ neglected, $Q(\bx, \bom) = \trans(\bx_0, \bx)\,\Ls(\bx_0, \bom)$ and
%
\begin{equation}
	\label{eqn:rte3d_deltaL_sharp}
		\begin{split}
		& \Delta\Lv(\bx, \bom) = \trans(\bx_0, \bx) \,\Delta\Ls(\bx_0, \bom)\\
		\!\!\!\!=\;& \trans(\bx_0, \bx) \left( \tint_{\Sph} \Delta\bsdf(\bx_0, -\bom', \bom)\,\Lv(\bx_0, \bom')\,\D\bom' + \Delta\Le(\bx_0, \bom) \right),
	\end{split}
\end{equation}
%
where $\Delta\bsdf$ and $\Delta\Le$ denote the difference between the BSDFs and emitted radiances evaluated with the normals $\gnorm^-$ and $\gnorm^+$ (see Figure~\ref{fig:deltaL_sharp}).

\begin{figure}[t]
	\centering
	\includegraphics[height=0.9in]{delta_L_sharp.pdf}
	\caption{\label{fig:deltaL_sharp}
		Calculating $\Delta\Lv$ resulting from a sharp edge.
		In this case, $\Delta\Lv$ is solely determined by the difference between the BSDF and the surface emission at $\bx_0$ evaluated at the two normals $\gnorm^-$ and $\gnorm^+$, as expressed in Eq.~\eqref{eqn:rte3d_deltaL_sharp}. 
	}
\end{figure}

