Compute step response from impulse response of continuous-time LTI system
Since your question is not about how to obtain the answer, but specifically about the reason why the unit step $u(t)$ appears in it, I assume that you know that the response to a step input is computed by convolving the impulse response with a step:

$$y(t)=\int_{-\infty}^{\infty}h(\tau)u(t-\tau)d\tau=\int_{-\infty}^{t}h(\tau)d\tau\tag{1}$$

Since $h(t)$ is causal, i.e., $h(t)=0$ for $t&lt;0$, the lower integration limit can be changed to $0$. This implies that if in $(1)$ the upper integration limit $t$ is less than $0$, the result must be zero. Consequently, the response to a step input can be written as

$$y(t)=u(t)\cdot\int_0^{t}h(\tau)d\tau\tag{2}$$

Evaluating $(2)$ gives you the final answer, which includes the step function $u(t)$.

EDIT: Triggered by Dilip Sarwate's comment I add some extra explanation. Note that the step function in $(2)$ is actually redundant because for $t&lt;0$ the integral in $(2)$ is zero because for $t&lt;0$ we have

$$\int_0^{t}h(\tau)d\tau=-\int_{-|t|}^0h(\tau)d\tau=0,\quad t&lt;0$$

because $h(t)=0$ for $t&lt;0$. However, if $h(t)=f(t)u(t)$ with some function $f(t)$ that is not zero for $t&lt;0$ (in the given example we have $f(t)=e^{-4t}$), we cannot just replace the lower integration limit in $(1)$ by $0$ and leave out the step function:

$$\int_{-\infty}^th(\tau)d\tau=\int_{-\infty}^tf(\tau)u(\tau)d\tau=\int_{0}^tf(\tau)u(\tau)d\tau\neq \int_{0}^tf(\tau)d\tau$$

because now we generally have for $t&lt;0$

$$\int_{0}^tf(\tau)d\tau=-\int_{-|t|}^0f(\tau)d\tau\neq 0$$

The correct way to evaluate the integral for $h(t)=f(t)u(t)$ is

$$\int_{-\infty}^th(\tau)d\tau=u(t)\cdot\int_0^tf(\tau)d\tau$$