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]{\Delta}}{3}$$ or $$\sqrt[3]{\Delta} = - b_4 + 3(x_i x_j + x_k x_{\ell}).$$ In an earlier draft I said that $64 a_4^3 - 432 a_6^2 = 16 (4 a_4^3 - 27 a_6^2)$ was $16 \Delta$, but apparently the $\Delta$ that shows up in modular forms is $16$ times the classical discriminant of the cubic equation. Thanks to Sylivan JULIEN for pointing this out.
Here is a conceptual explanation for a big piece of this. For any $x_1$, $x_2$, $x_3$, $x_4$, note that $$\frac{(x_1 x_2 + x_3 x_4) - (x_1 x_3 + x_2 x_4)}{(x_1 x_2 + x_3 x_4) - (x_1 x_4 + x_2 x_3)} = \frac{(x_1 - x_4)(x_2 - x_3)}{(x_1 - x_4)(x_2 - x_3)}$$ which is the cross ratio $c(x_1, x_2 : x_3, x_4)$. We want to show that this ratio is a cube root of unity, so we want to show that the cross ratio of $x_1$, $x_2$, $x_3$, $x_4$ is a cube root of unity.
This computation turns out to be easiest when the cube is not in Weierstrass form but Hessian form: $X^3+Y^3+Z^3 = a XYZ$. The flexes of this curve are the $9$ points with homogenous coordinates $(1:-\zeta:0)$, $(0:1:-\zeta)$ and $(-\zeta:0:1)$ with $\zeta^3=1$. If we take $(1:-1:0)$ to be the origin of our curve, then negation is $(X:Y:Z) \mapsto (Y:X:Z)$ and we can take the quotient by negation to be given by the rational map $(X:Y:Z) \mapsto \tfrac{X+Y}{Z}$. The $8$ non-identity flexes map to $\infty$ and to the $3$ cube roots of $-1$, whose cross ratio is as required.
Warning: The next part uses modular form language that I am not completely comfortable with.
The computation with cross ratios shows that there are $A$ and $B$ such that the three values of $x_i x_j + x_k x_{\ell}$ are of the form $A + \zeta B$, for $\zeta$ running over the cube roots of $1$. I think it shouldn't be bad to show that $A$ and $B$ are modular forms of level $3$ and weight $4$. I don't have a conceptual explanation for why $B$ should be a cusp form but, if you believe it is, then I think it must be a multiple of $\Delta^{1/3}$.
As for other power of $\Delta$, I think that $\Delta^{1/6}$ is a cusp form for $\Gamma(6)$ of weight $2$. The $x$-coordinates of the $6$-torsion points (including the $2$-torsion and $3$-torsion) should be related to modular forms for $\Gamma(6)$ of weight $2$, so maybe we can find some clever linear combination of them which equals $\Delta^{1/6}$.