Path-ordered integrals with dirac deltas When calculating integrals with delta functions, the integration domains are usually rectangular. However, recently, I encountered a set of time ordered integrals over Dirac deltas of the following form:
$$ \int_0^1 dt_1 \int_0^{t_1} dt_2 \int_0^{t_2} dt_3 \int_0^{t_3} dt_4 K(t_1,t_2,t_3,t_4) \delta(t_1-t_2)\delta(t_3-t_4) $$
$$ \int_0^1 dt_1 \int_0^{t_1} dt_2 \int_0^{t_2} dt_3 \int_0^{t_3} dt_4 K(t_1,t_2,t_3,t_4) \delta(t_1-t_3)\delta(t_2-t_4) $$
$$ \int_0^1 dt_1 \int_0^{t_1} dt_2 \int_0^{t_2} dt_3 \int_0^{t_3} dt_4 K(t_1,t_2,t_3,t_4) \delta(t_1-t_4)\delta(t_3-t_2) $$
Where $K(t_1,t_2,t_3,t_4)$ is some smooth function of the arguments. 
For convenience, I set $K=1$ and evaluated these integrals in mathematica. All three of them are zero! This answer makes no sense to me, because suppose I add up these three integrals (for $K=1$), I would get:
$$ \int_0^1 dt_1 \int_0^{t_1} dt_2 \int_0^{t_2} dt_3 \int_0^{t_3} dt_4 \bigg(\delta(t_1-t_2)\delta(t_3-t_4)+\delta(t_1-t_3)\delta(t_2-t_4)+\delta(t_1-t_4)\delta(t_2-t_3)\bigg) $$
This expression is manifestly symmetric wrt to all permutation of $t_i$'s. Therefore, the time ordered integral can be transformed into a regular integral over rectangular domains:
$$ \frac{1}{24} \int_0^1 dt_1 \int_0^{1} dt_2 \int_0^{1} dt_3 \int_0^{1} dt_4 \bigg(\delta(t_1-t_2)\delta(t_3-t_4)+\delta(t_1-t_3)\delta(t_2-t_4)+\delta(t_1-t_4)\delta(t_2-t_3)\bigg) $$
This integral evaluates to $\frac{1}{8}$ instead of zero. 
So the question is: which one is the correct answer? And more generally, how should I think about Dirac deltas rigorously in time-ordered integrals?
 A: *

*We can replace simplex integration region with a rectangular domain $[0,1]^4$ if we multiply the integrand with 
$$ \theta(t_1\!-\!t_2)\theta(t_2\!-\!t_3)\theta(t_3\!-\!t_4), \tag{A}$$
where $\theta$ denotes the Heaviside step function. Multiplying distributions with discontinuous functions is not mathematically well-defined.  

*Nevertheless, we cannot resist to heuristically calculate what would happen with OP's integral in the case $K=1$ if we assume that $\theta(0)=\frac{1}{2}$. 


*

*If we first integrate over $t_1$ & $t_4$ the first integrand becomes $$\theta(0)\theta(t_2\!-\!t_3)\theta(0).$$ The integrations over $t_2$ & $t_3$ then yield $\frac{1}{8}$.

*If we first integrate over $t_3$ & $t_4$ the second integrand becomes $$\theta(t_1\!-\!t_2)^2\theta(t_2\!-\!t_1)
~=~\frac{1}{8}\delta_{t_1,t_2}~=~0\text{ a.e.}$$ 
Here $\delta_{t_1,t_2}$ is the Kronecker delta function$^1$, which is zero almost everywhere (a.e.). The integrations over $t_1$ & $t_2$ then yield $0$.

*Similarly, if we first integrate over $t_3$ & $t_4$ the third integrand becomes 
$$\theta(t_1\!-\!t_2)\theta(0)\theta(t_2\!-\!t_1)
~=~\frac{1}{8}\delta_{t_1,t_2}~=~0\text{ a.e.}$$
The integrations over $t_1$ & $t_2$ then yield $0$.
Hence the total integral is heuristically $\frac{1}{8}+0+0=\frac{1}{8}$, which agrees with OP's symmetrization argument. 
--
$^1$The Kronecker delta function $\delta_{t_1,t_2}$ should not be confused with the Dirac delta distribution $\delta(t_1\!-\!t_2)$.
A: One trick:
Whenever you have $\delta (t_1 - t_2) \delta (t_3 - t_4)K(t_1, t_2, t_3, t_4)$ you can replace the integrals to be over just two variables, $t_1$ and $t_3$, say, and argument $K(t_1, t_1, t_3, t_3)$.
