I want to construct the heat kernel $k_t(x,y)$, which is the fundamental solution to the heat equation $(\partial_t + \Delta_x)u(t,x) = 0$. There is something that I don't understand about using Duhamel's principle in this construction. (My context is that $\Delta$ is the Laplace-Beltrami operator---the geometer's Laplacian, i.e., a positive operator---on a compact manifold, but I think my question also makes sense in other settings.)
I have a parametrix $p^0_t(x,y)$ that is a "zeroth-order" approximation to $k_t(x,y)$, and I want to use Duhamel's principle to improve my approximation. Let's define the zeroth error $E_t^0(x,y)$ by $$ E_t^0(x,y) := (\partial_t + \Delta_x) p_t^0(x,y). $$ Duhamel's principle gives that the actual heat kernel satisfies $$ k_t = p_t^0 - \int_0^t k_{t-s_1} E^0_{s_1} \, ds_1. ~~~~~(\ast)$$ (Henceforth I will suppress the dependence on the spatial variables; by $k_{t-s_1} E^0_{s_1}$, I really mean the composition of operators that has integral kernel $\int k_{t-s_1}(x,z) E^0_{s_1}(z,y) \, dz $.) This suggests the following improvement of $p_t^0$, which I will call $p_t^1$: $$p_t^1 := p_t^0 - \int_0^t p^0_{t-s_1} E^0_{s_1} \, ds_1. $$ One can compute that the error associated to $p_t^1$, which I will call $E_t^1$, is \begin{align*} E_t^1 :=& (\partial_t + \Delta_x) p_t^1(x,y) \\ =& - \int_0^t E^0_{t-s_1} E^0_{s_1} \, ds_1. \end{align*} Applying Duhamel's principle again, we have \begin{align} k_t &= p_t^1 - \int_0^t k_{t-s_1} E^1_{s_1} \, ds_1 \\ &= p_t^1 + \int_0^t k_{t-s_1} \int_0^{s_1} E^0_{s_1-s_2} E^0_{s_2} \, ds_2 \, ds_1 \end{align} (I think one could iterate this procedure to obtain successive approximations $p_t^k$; this seems to be similar to the approach taken in, for example, Friedman's parabolic pde book.)
On the other hand, we could have taken a different approach: beginning with $(\ast)$, we have \begin{align} k_t &= p_t^0 - \int_0^t k_{t-s_1} E^0_{s_1} \, ds_1 \\ &= p_t^0 - \int_0^t p^0_{t-s_1} E^0_{s_1} \, ds_1 - \int_0^t (k_{t-s_1} - p^0_{t-s_1}) E^0_{s_1} \, ds_1 \\ &= p_t^1 + \int_0^t \int_0^{t-s_1} k_{t-s_1 - s_2} E^0_{s_2} E^0_{s_1} \, ds_2 \, ds_1, \end{align} where for the last equality I have again used $(\ast)$.
I seem to have obtained two expressions for the correction term when I approximate $k_t$ by $p_t^1$, which are not obviously equal: $$ \int_0^t k_{t-s_1} \int_0^{s_1} E^0_{s_1-s_2} E^0_{s_2} \, ds_2 \, ds_1 ~~ \text{ and } ~~ \int_0^t \int_0^{t-s_1} k_{t-s_1 - s_2} E^0_{s_2} E^0_{s_1} \, ds_2 \, ds_1, $$
Question: Are these two expressions equal, and if so, is there an easy way to see this, other than the argument above? I have a feeling there is a simple substitution that I am missing. Or have I made a silly mistake somewhere?