Series of nested integrals

I’m trying to calculate the following series of nested integrals with $\varepsilon(t)$ being a real function.

$$\sigma = 1 + \int\nolimits_{t_0}^t\mathrm dt_1 \, \varepsilon(t_1) + \int_{t_0}^t\mathrm dt_1 \int_{t_0}^{t_1}\mathrm dt_2 \,\varepsilon(t_1)\, \varepsilon(t_2) + \int_{t_0}^t\mathrm dt_1 \int_{t_0}^{t_1}\mathrm dt_2 \int_{t_0}^{t_2}\mathrm dt_3\, \varepsilon(t_1)\, \varepsilon(t_2)\, \varepsilon(t_3) + \cdots \;.$$

The result should be

$$\sigma = \exp\left(\int_{t_0}^t\mathrm dt_1\, \varepsilon(t_1)\right) = \sum_{i=0}^\infty \frac1{i!} \left(\int_{t_0}^t\mathrm dt_1 \,\varepsilon(t_1)\right)^i \;.$$

However, comparing the series term by term I already fail to prove the equivalence for the third term. Can someone clear this up for me? Can the series be rewritten to an exponential after all? I recall from my quantum mechanics course that if $\varepsilon(t)$ was an operator, non-commuting with itself for different times, then the formal result would be

$$\sigma = T_c \exp\left(\int_{t_0}^t\mathrm dt_1\, \varepsilon(t_1)\right) \;,$$

with $T_c$ being the usual time-ordering operator. However, as I said in my case $\varepsilon(t)$ is a plain function in the real domain.

Solutions Collecting From Web of "Series of nested integrals"

Expanding on my comment:

The function $$\sigma_1(t)= \exp\left[\int_{t_0}^t dt' \varepsilon(t')\right]$$ fulfills the following differential equation
$$\sigma_1'(t)= \varepsilon(t) \sigma_1(t)$$ with the boundary condition $\sigma_1(t_0) = 1.$ We will show in a next step that $$\sigma_2 (t) = 1 + \int_{t_0}^t dt_1 \varepsilon(t_1) + \int_{t_0}^t dt_1 \int_{t_0}^{t_1} dt_2 \varepsilon(t_1) \varepsilon(t_2) + \dots \;
obeys the same differential equation. Because the solution to this linear differential equation is unique, it follows that $\sigma_1(t) = \sigma_2(t)$.

Taking derivative of $\sigma_2 (t)$, we recover (term by term)
$$\sigma_2' (t) = \varepsilon(t) + \varepsilon(t) \int_{t_0}^t dt_1 \varepsilon(t_1) + \varepsilon(t)\int_{t_0}^t dt_1 \int_{t_0}^{t_1} dt_2 \varepsilon(t_1) \varepsilon(t_2) + \dots = \varepsilon(t)\sigma_2(t).$$
The boundary condition $\sigma_2(t_0) = 1$ also follows easily…

In conclusion, you prove with this that
$$\frac1{n!} \left(\int_{t_0}^t\mathrm dt \,\varepsilon(t)\right)^n = \int_{t_0}^t\mathrm dt_1 \int_{t_0}^{t_1}\mathrm dt_2 \cdots \int_{t_0}^{t_{n-1}}\mathrm dt_n\, \varepsilon(t_1)\, \varepsilon(t_2) \cdots \varepsilon(t_n). $$

Rewriting your definition as $\sigma=1+I_1+I_2+I_3+\ldots$, the explanation you are after is that $(I_1)^n=(n!)I_n$ for every $n\ge1$, since then
But the fact that $(I_1)^n=(n!)I_n$ is obvious: $(I_1)^n$ is the integral of the symmetric function
over the cube $K_n=[t_0,t]^n$. Likewise, $I_n$ is the integral of $e_n$ over the simplex $\Delta_n\subset K_n$ made of the points $(t_1,\ldots,t_n)$ such that $t_0\le t_1\le t_2\le\cdots\le t_n\le t$.

Recall that the symmetric group $\mathfrak{S}_n$ acts on $K_n$ as follows: for every $s$ in $\mathfrak{S}_n$ and $(t_1,\ldots,t_n)$ in $K_n$,

Now, $K_n$ is the union of the $n!$ simplexes $s\cdot\Delta_n$ for $s$ in $\mathfrak{S}_n$. The function $e_n$ is symmetric hence its integral on $s\cdot\Delta_n$ is independent on $s$. The simplexes $s\cdot\Delta_n$ intersect on zero measure sets hence $(I_1)^n$ is the sum over $s$ in $\mathfrak{S}_n$ of the integrals of $e_n$ on $s\cdot\Delta_n$. Each of these integrals is equal to $I_n$ and there are $n!$ of them hence you are done.