Here is an algebraic proof, without a ton of insight. To make life easier, let's put our curve into reduced Weierstrass form, $y^2 = x^3+a_4 x + a_6$. We put $f(x) = x^3+a_4 x + a_6$. The $3$-torsion points are the flexes, meaning the points where
$\tfrac{d^2 y}{(dx)^2}=0$. We compute
$$\frac{d^2 y}{(dx)^2} = \frac{d^2 }{(dx)^2} f(x)^{1/2} = (1/2) f'' f^{-1/2} - (1/4) (f')^2 f^{-3/2}=\frac{2 f'' f - (f')^2}{4 f^{3/2}}.$$
So the $x$-coordinates of the $3$-torsion points are the roots of 
$$2 f'' f - (f')^2 = 3 x^4 + 6 a_4 x^2 + 12 a_6 x - a_4^2.$$
We deduce that the elementary symmetric polynomials in $(x_1, x_2, x_3, x_4)$ take the values
$$e_1(x)=0,\ e_2(x)=2 a_4,\ e_3(x) = -4 a_6,\ e_4(x) = -a_4^2/3 . (\ast)$$

Expanding $(y-x_1 x_2 - x_3 x_4)(y-x_1 x_3 - x_2 x_4)(y-x_1 x_4 - x_2 x_3)$ gives a polynomial in $y$ whose coefficients are elementary symmetric polynomials in $(x_1, x_2, x_3, x_4)$. By the fundamental theorem of symmetric polynomials, we can write the coefficients of this cubic as polynomials in the $e_j(x)$, and then plug in the formulas from $(\ast)$. (If you use Mathematica, the `SymmetricReduction` command will do this for you.) I get that this cubic is
$$y^3 - 2 a_4 y^2 + \tfrac{4}{3} a_4^2 y - \tfrac{8}{3} a_4^3 - 16 a_6^2$$
$$=y^3 - b_4 y^2 + \tfrac{1}{3} b_4^2 y - \tfrac{1}{3} b_4^3 - 16 a_6^2 = (y-b_4/3)^3 - \tfrac{8}{27} b_4^3 - 16 a_6^2.$$

So the values of $x_i x_j + x_k x_{\ell}$ are
$$\tfrac{b_4}{3} + \sqrt[3]{\tfrac{8}{27} b_4^3 - 16 a_6^2} = \frac{b_4 + \sqrt[3]{64 a_4^3 - 432 a_6^2}}{3} = \frac{b_4 + \sqrt[3]{16 \Delta}}{3}$$
or
$$\sqrt[3]{16 \Delta} = - b_4 + 3(x_i x_j + x_k x_{\ell}).$$
I'm not sure why I have a $16$ that you don't, but the formulas match otherwise.