Although this is a standard derivation, you frequently don't see it in introductory electromagnetism courses, maybe because those courses shy away from the heavy use of vector calculus. Here's the usual approach. We'll find a wave equation from Maxwell's equations.
Start with
$\nabla \times \vec{E} = -\frac{\partial\vec{B}}{\partial t}$.
Take a partial derivative of both sides with respect to time. The curl operator has no partial with respect to time, so this becomes
$\nabla \times \frac{\partial\vec{E}}{\partial t} = -\frac{\partial^2\vec{B}}{\partial t^2}$.
There's another of Maxwell's equations that tells us about $\partial\vec{E}/\partial t$.
$\nabla \times \vec{B} = \mu_0\epsilon_0\frac{\partial \vec{E}}{\partial t}$
Solve this for $\partial\vec{E}/\partial t$ and plug into the previous expression to get
$\nabla \times \frac{(\nabla \times \vec{B})}{\mu_0\epsilon_0} = -\frac{\partial^2 \vec{B}}{\partial t^2}$
the curl of curl identity lets us rewrite this as
$\frac{1}{\mu_0 \epsilon_0}\left(\nabla(\nabla \cdot \vec{B}) - \nabla^2\vec{B}\right) = -\frac{\partial^2 \vec{B}}{\partial t^2}$
But the divergence of the magnetic field is zero, so kill that term, and rearrange to
$\frac{-1}{\mu_0 \epsilon_0}\nabla^2\vec{B} + \frac{\partial^2 \vec{B}}{\partial t^2} = 0$
This is the wave equation we're seeking. One solution is
$\vec{B} = B_0 e^{i (\vec{x}\cdot\vec{k} - \omega t) }$.
This represents a plane wave traveling in the direction of the vector $\vec{k}$ with frequency $\omega$ and phase velocity $v = \omega/|\vec{k}|$. In order to be a solution, this equation needs to have
$\frac{\omega^2}{k^2} = \frac{1}{\mu_0\epsilon_0}$.
Or, setting $v = 1/\sqrt{\mu_0\epsilon_0}$
$\frac{\omega}{k} = v$
This is called the dispersion relation. The speed that electromagnetic signals travel is given by the group velocity
$\frac{d\omega}{d k} = v$
So electromagnetic signals in a vacuum travel at speed $c = 1/\sqrt{\mu_0\epsilon_0}$.
Edit
You can follow the same steps to derive the wave equation for $\vec{E}$, but you will have to assume you're in free space, i.e. $\rho = 0$.
Edit
The curl of the curl identity was wrong, there's a negative number in there
The Maxwell stress tensor is introduced as analogue to the stress tensor in continuum mechanics, and its form is derived from the equation
$$
\frac{d}{dt} \int_V (\boldsymbol{\mathscr{p}} + \boldsymbol{\mathscr{g}}) \,\mathrm dV = \oint_{S} d\mathbf S \cdot \mathbf T
$$
where $\boldsymbol{\mathscr{p}}$ is density of momentum of matter and $\boldsymbol{\mathscr{g}}$ is density of momentum of the EM field. The right-hand side is a surface integral over $S$, the boundary of the region $V$ and it looks like total surface force in continuum mechanics. The region and thus also its surface is arbitrary, the equation is valid for any choice. If the boundary is in free space, obviously there is no "surface force" in the original sense, but the equation has the same form as if there was - as in continuum mechanics.
Alternatively, one can interpret $d\mathbf S\cdot\mathbf T$ on the right-hand side not as a surface force per unit area, but as EM momentum that enters the region $V$ through $dS$ per unit time. This is perhaps better as we do not have to talk about "force of tension" acting in free space (on what? is a good question). But both views are commonly used.
EDIT:
In case everything is static, force on regular charge distribution $\rho(\mathbf x)$ (first charged particle) that is contained within region $V$ is
$$
\mathbf F = \int_V \rho(\mathbf x)\mathbf E(\mathbf x) d^3\mathbf x
$$
This can be transformed using Maxwell's equations and vector calculus into
$$
\mathbf F = \oint_S d\mathbf S \cdot \mathbf T
$$
where $S$ is boundary surface of $V$. This only signifies that electrostatic force on a charged body can be calculated either as sum of elementary forces acting "locally" or "in bulk" as "volume forces", or one can alternatively calculate the same force from the field at distance surface that encloses the body as a sum of fictitious surface forces. The latter case reminds of continuum stress force, but the connection is only mathematical - physically, there is no force on the surface since there is no matter there.
Sidenote: All this is usually derived only for regular distributions, where $\rho$ is bounded. The derivation does not work for point charges, because for them the left-hand integral does not have sense - they are a special case that leads to different theory. However, incidentally the right-hand side is still valid for point particles. This is because there is similar theorem that can be derived for point particles, with slightly different EM stress tensor, which somewhat surprisingly gives the same surface integral. This is easily seen from the fact that the force between two charged bodies does not depend on whether the bodies are points or say, uniformly charged spheres.
Back to the OP's questions, all this equivalence between the two ways to express the force breaks down once the fields cease to be static; then the presence of EM momentum in $V$ cannot be neglected. Then the force acting on a charged body is not given solely by the right-hand type of integral, but it is believed that the original equation for the force - the first equation of the OP (integral of a local expression) is still correct for extended bodies.
Best Answer
The stress tensor is defined in such a way that taking its divergence brings you back to the force density (apart from the cross product term which corresponds to the time derivative of the Poynting vector). In the case of $\mathbf{B}=\mathbf{0}$, we take the divergence in index notation, i.e.
$$f_j=\partial_i\sigma_{ij}=\epsilon_0\left(\partial_i(E_iE_j)-\frac12\delta_{ij}\partial_i(E_kE_k)\right)=\epsilon_0\left(\partial_iE_iE_j+E_i\partial_iE_j-\frac12\partial_j(E_kE_k)\right).$$
In order to see that this expression corresponds to the first term in your expression 3., we recast this in index-free notation:
$$\mathbf{f}=\epsilon_0\left((\boldsymbol{\nabla}\cdot\mathbf{E})\mathbf{E}+(\mathbf{E}\cdot\boldsymbol{\nabla})\mathbf{E}-\frac12\boldsymbol{\nabla}(\mathbf{E}\cdot\mathbf{E})\right).$$
Making use of the identity
$$\frac{1}{2} \boldsymbol{\nabla} \left( \mathbf{A}\cdot\mathbf{A} \right) = \mathbf{A} \times (\boldsymbol{\nabla} \times \mathbf{A}) + (\mathbf{A} \cdot \boldsymbol{\nabla}) \mathbf{A}$$
and the Maxwell equation (for $\mathbf{B}=\mathbf{0}$)
$$\boldsymbol{\nabla}\times\mathbf{E}=\mathbf{0},$$
we find that
$$\mathbf{f}=\epsilon_0(\boldsymbol{\nabla}\cdot\mathbf{E})\mathbf{E}.$$
From this we can see that it does not matter at which stage of the derivation the magnetic field is omitted.