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{Z}{X+Y}$ (this also takes $(1:-1:0)$ to $\infty$, so it might be something like the $x$-coordinate). The $8$ non-identity flexes map to $0$ and to the $3$ cube roots of $-1$, whose cross ratio is as required.
The next part of this answer addresses the question of whether similar
formulas exist for other powers of $\Delta$. I will interpret
"similar" as "a polynomial in the $x$ and $y$-coordinates of the
$N$-torsion points".
For a modular form $\phi$ of weight $k$, and $g = \left[
\begin{smallmatrix} a&b \\ c&d \end{smallmatrix} \right]$ in
$SL_2(\mathbb{Z})$, let $(g^{\ast} \phi)(z) = (cz+d)^{-k} \phi\left(
\tfrac{az+b}{cz+d} \right)$.
Recall that $\Delta^{1/24}$ is the Dedekind $\eta$ function, which is a modular form of weight $1/2$. I find $1/2$ integer weights confusing, so I'll only look at even powers of $\eta$. Thus, I'll be asking whether $\Delta^{k/12} = \eta^{2k}$ could be a polynomial in $x$ and $y$-coordinates. Note that $\Delta^{k/12}$ has weight $k$.
For any $g \in SL_2(\mathbb{Z})$, we have $g^{\ast} \eta^2 = \chi(g)
\eta^2$, where $\chi$ is a character from $SL_2(\mathbb{Z})$ to the
twelfth roots of unity, an explicit formula for which can be found in
the Wikipedia article linked before. So $g^{\ast}
\Delta^{k/12} = \chi(g)^k \Delta^{k/12}$.
Note that $\chi$ factors through the quotient $SL_2(\mathbb{Z}/12
\mathbb{Z})$.
Lemma For $1 \leq k \leq 12$, let $\phi$ be a cusp form of weight
$k$ (and some level) obeying $g^{\ast} \phi = \chi(g)^k \phi$. Then
$\phi$ is a scalar multiple of $\Delta^{k/12}$.
Proof: Since $\Delta^{1/12}$ is nowhere vanishing on the upper
half plane, the ratio $\phi / \Delta^{k/12}$ is holomorphic, and is
invariant for $SL_2(\mathbb{Z})$. Therefore, it is a polynomial in
$j$, and we have $\phi = \Delta^{k/12} \sum_{e=0}^d c_e j^e$ for some
polynomial in $j$. But, looking at $q$ series, the leading term of the
right hand side is $c_d q^{k/12 - d}$, and the leading power of $q$ on
the left hand side is positive, so $d$ must be $0$. $\square$
Now, let $\psi$ be a modular form of weight $1 \leq k \leq 11$. Let
$N$ be the LCM of $12$ and the level of $\psi$. Define
$$R \psi = \sum_{g \in SL_2(\mathbb{Z}/N \mathbb{Z})} \chi(g)^k g^{\ast} \psi.$$
(Since $12$ divides $N$, it makes sense to evaluate $\chi$ on
$SL_2(\mathbb{Z}/N \mathbb{Z})$; since the level of $\psi$ divides
$N$, it makes sense to talk about $g^{\ast} \psi$ similarly.)
Then $R \psi$ will be a modular form of the same weight $k$, obeying
$g^{\ast} \phi = \chi(g)^k \phi$. I claim furthermore that it will be
a cusp form. Proof: Let $\Gamma'$ be the kernel of $\chi$. A set of
coset representatives $SL_2(\mathbb{Z}/N \mathbb{Z})/\Gamma'$ is given
by $\left[ \begin{smallmatrix} 1 & x \\ 0 & 1 \\ \end{smallmatrix}
\right]$ for $0 \leq x \leq 11$. So write the sum as
$$\sum_{x=0}^{11} \chi\left( \begin{smallmatrix} 1 & x \\ 0 & 1 \\
\end{smallmatrix} \right)^k \left[ \begin{smallmatrix} 1 & x \\ 0 & 1 \\ \end{smallmatrix}
\right]^{\ast} \sum_{g \in \Gamma'/\Gamma(N)} g^{\ast} \phi.$$
The inner sum is a modular form for $\Gamma'$, which has only one
cusp. The action of $ \left[ \begin{smallmatrix} 1 & x \\ 0 & 1 \\ \end{smallmatrix}
\right]$ takes that cusp to itself. So the value of the whole sum at
that cusp is $\sum_{x=0}^{11} \chi\left( \begin{smallmatrix} 1 & x \\ 0 & 1 \\
\end{smallmatrix} \right)^k$ times the value of the inner sum at the
cusp, and that inner sum is $0$.
So, if $\psi$ is any modular form of weight $1 \leq k \leq 11$, then
$R \psi$ is a scalar multiple of $\Delta^{k/12}$. Of course, that
scalar might be $0$, but we can hope!
We now want to know that the $x$-coordinates of $N$-torsion points are
modular forms of level $N$ and weight $2$, and the $y$-coordinates are
modular forms of level $N$ and weight $3$.
So we can take a polynomial in $x$'s and $y$'s of appropriate weight,
apply the $R$-operator and hope.
If we are going to have a chance, we better make sure the level is
high enough. $\Delta^{k/12}$ has level $\tfrac{12}{GCD(k,12)}$, so we
should take $N$ divisible by this. The most obvious thing to try is to
take $N = \tfrac{12}{GCD(k,12)}$.
We identify the $N$-torsion points with $(\mathbb{Z}/N
\mathbb{Z})^2$. For $(a,b) \in (\mathbb{Z}/N
\mathbb{Z})^2$, we denote the coordinates of the corresponding torsion
point as $(x(a,b), y(a,b))$. Note that the action of
$SL_2(\mathbb{Z}/N \mathbb{Z})$ on $x(a,b)$ and $y(a,b)$ is precisely
the action on the vectors $(a,b)$. (Row or column vectors? I don't
feel like working that hard.)
I'll present cases in order of complexity:
$k=6$: $\Delta^{1/2}$ has level $2$, so we work with $2$-torsion
points. We want a polynomial of weight $6$, so we try cubics in the
$x$-variables. Applying $R$ to $x(1,0)^2 x(0,1)$, we obtain
$$x(1,0)^2 x(0,1) + x(0,1)^2 x(1,1) + x(1,1)^2 x(1,0)-x(0,1)^2 x(1,0)
- x(1,1)^2 x(0,1) - x(1,0)^2 x(1,1)$$
$$=(x(1,0) - x(0,1)) (x(1,0) - x(1,1)) (x(0,1) - x(1,1)).$$
Sure enough,
$$\Delta = (x(1,0) - x(0,1))^2 (x(1,0) - x(1,1))^2 (x(0,1) -
x(1,1))^2.$$
This is far from the shortest way to obtain this identity, but it
works.
$k=4$ This is the one the OP started with. This time, $\Delta^{1/3}$
has level $3$, so we work with $3$-torsion points. We want a
polynomial of weight $4$, so we try quadratics in the $x$-variables.
We have $x(a,b) = x(-a,-b)$, so we index the $x$-variables by the
points of $\mathbb{P}^1(\mathbb{F}_3)$, written in homogenous
coordinates $x(a:b)$.
The action of $SL_2(\mathbb{Z})$ on the $4$ points of
$\mathbb{P}^1(\mathbb{F}_3)$ is by the alternating subgroup $A_4$. So
we want to take a quadratic monomial, $x(1:0) x(0:1)$ and average it
with an order $3$ character of $A_4$. Writing $\omega$ for a cube root
of unity, and not working hard enough to figure out which one I mean,
we get that $\Delta^{1/3}$ is proportional to
$$\left( x(1:0) x(0:1) + x(1:1) x(1:2) \right) + \omega \left( x(1:0)
x(1:1) + x(0:1) x(1:2) \right) + \omega^2 \left( x(1:0) x(1:2) +
x(0:1) x(1:1) \right)$$
as desired. Something else cute happens here: If we chose the other
power of $\omega$, we get $0$. So we can use this to rewrite the
formula in simpler ways.
$k=1$ No monomial in $x$'s and $y$'s can have weight $1$.
$k=2$ or $k=3$. So we want linear monomials in $x$'s or
$y$'s. However, every torsion point is stabilized by some conjugate of
$\left[ \begin{smallmatrix} 1 & \ast \\ 0 & 1 \\ \end{smallmatrix}
\right]$, and averaging over this stabilizing subgroup gives $0$, so
we just get $0$ if we apply $R$ to any $x(a,b)$ or $y(a,b)$. (François
Brunault, in comments, states that something stronger is true: No
linear combination of $x$'s or $y$'s is ever a cusp form. I think I've
reconstructed the proof, but I'll leave it to him.)
That finishes the divisors of $12$. Also, $8=4+4$ and $10=4+6$, so
$\Delta^{8/12}$ and $\Delta^{10/12}$ are products of things we already
have. There are two other cases I find interesting: $k=5$ and
$k=9$. (The case $k=7$ seems like just a messier version of $k=5$, and $11 =
5+6$, which is why I don't care so much about them.)
$k=5$: We want $12$-torsion points, and we want weight $5$, so we want
products of an $x$ and a $y$. If I didn't screw up, if $(a,b)$ and
$(c,d)$ fail to generate $(\mathbb{Z}/12 \mathbb{Z})^2$, then $R
x(a,b) y(c,d)=0$. However, if they generate, then the result sure
doesn't look like zero! I get that
$$\sum_{\left[ \begin{smallmatrix} a&b \\ c&d \\
\end{smallmatrix} \right] \in SL_2(\mathbb{Z}/12 \mathbb{Z}) }
\chi \left( \begin{smallmatrix} a&b \\ c&d \end{smallmatrix} \right)^5 x(a,b) y(c,d)$$
should be a scalar multiple of $\Delta^{5/12}$! Does anyone have the
computational chops to work out which one?
$k=9$ This time, we can try $4$-torsion points and polynomials of
weight $9$. There are a lot of choices, but I think a very natural one
to try is $y(1,0) y(0,1) y(3,3)$. Using the identity $y(a,b) = - y(-a,
-b)$, there are only $8$ terms, which I think are the following:
$$y(1,0) y(0,1) y(1,1)+i y(1,2) y(0,1) y(1,1)-i y(1,0) y(2,1) y(1,1)-y(1,2) y(2,1) y(1,1)-i
y(1,0) y(0,1) y(1,3)-y(1,2) y(0,1) y(1,3)+y(1,0) y(2,1) y(1,3)+i y(1,2) y(2,1) y(1,3).$$
So, this is supposed to be a multiple of $\Delta^{3/4}$. Which one?
