Why is the catalecticant invariant under coordinate changes? Let $\mathbf{k}$ be a commutative $\mathbb{Q}$-algebra. (We could play the
same game over any commutative ring $\mathbf{k}$, but this would be a bit more
technical, so let me avoid it.)
Fix a nonnegative integer $n$. A binary form shall mean a homogeneous
polynomial $f=f\left(  x,y\right)  \in\mathbf{k}\left[  x,y\right]  $. Any
binary form $f\left(  x,y\right)  $ of degree $2n$ can be written uniquely in
the form $f\left(  x,y\right)  =\sum\limits_{i=0}^{2n}\dbinom{2n}{i}a_{i}
x^{2n-i}y^{i}$ for some scalars $a_{0},a_{1},\ldots,a_{2n}\in\mathbf{k}$
(which are simply the coefficients of $f$, rescaled by binomial coefficients).
Let $\mathcal{F}_{2n}$ denote the $\mathbf{k}$-module of all binary forms of
degree $2n$. The matrix ring $\mathbf{k}^{2\times2}$ acts on the $\mathbf{k}
$-module $\mathcal{F}_{2n}$ by the rule
\begin{equation}
\left(
\begin{array}[c]{cc}
a & b\\
c & d
\end{array}
\right)  \cdot f=f\left(  ax+cy,bx+dy\right)  .
\end{equation}
(In other words, a matrix $\left(
\begin{array}[c]{cc}
a & b\\
c & d
\end{array}
\right)  \in\mathbf{k}^{2\times2}$ acts on a binary form by substituting
$ax+cy$ and $bx+dy$ for the variables $x$ and $y$. This is called a "change of
variables", at least if the matrix is invertible.)
The catalecticant $\operatorname*{Cat}f$ of a binary form $f\left(
x,y\right)  $ of degree $2n$ is a scalar, defined as follows: Write $f\left(
x,y\right)  $ in the form $f\left(  x,y\right)  =\sum\limits_{i=0}^{2n}\dbinom{2n}
{i}a_{i}x^{2n-i}y^{i}$, then set
\begin{equation}
\operatorname*{Cat}f=\det\left(
\begin{array}[c]{cccc}
a_{0} & a_{1} & \cdots & a_{n}\\
a_{1} & a_{2} & \cdots & a_{n+1}\\
\vdots & \vdots & \ddots & \vdots\\
a_{n} & a_{n+1} & \cdots & a_{2n}
\end{array}
\right)  =\det\left(  \left(  a_{i+j-2}\right)  _{1\leq i\leq n+1,\ 1\leq
j\leq n+1}\right)  .
\end{equation}
A classical result, going back to
Sylvester (On the Principles of the Calculus of Forms) in 1852 (who, in turn, ascribes it
to Cayley), then says that $\operatorname*{Cat}f$ is a $\operatorname*{GL}
\nolimits_{2}\mathbf{k}$-invariant (in an appropriate sense of this word,
i.e., $\operatorname*{GL}\nolimits_{2}\mathbf{k}$ transforms
$\operatorname*{Cat}$ by multiplication with an appropriate power of the
determinant). More precisely, any binary form $f\left(  x,y\right)  $ of
degree $2n$ and every matrix $A\in\mathbf{k}^{2\times2}$ satisfy
\begin{equation}
\operatorname*{Cat}\left(  A\cdot f\right)  =\left(  \det A\right)  ^{n\left(
n+1\right)  }\cdot\operatorname*{Cat}f.
\label{catalecticant1}
\tag{1}
\end{equation}

Question. Is there a purely algebraic/combinatorial proof of \eqref{catalecticant1}?

These days, \eqref{catalecticant1} is usually seen as a consequence of the
fact that if $\mathbf{k}$ is an algebraically closed field of characteristic
$0$, then a binary form $f\left(  x,y\right)  $ of degree $2n$ satisfies
$\operatorname*{Cat}f=0$ if and only if $f$ can be written as $\sum_{j=1}
^{n}\left(  p_{j}x+q_{j}y\right)  ^{2n}$ for some $2n$ elements $p_{1}
,p_{2},\ldots,p_{n},q_{1},q_{2},\ldots,q_{n}\in\mathbf{k}$. This is proven,
e.g., in
J. M. Selig, Sylvester's Catalecticant. Once this fact is granted, the
equality \eqref{catalecticant1} follows using some technical Nullstellensatz
arguments (I believe so; the details are somewhat annoying to verify,
requiring e.g. to show that $\operatorname*{Cat}f$ is irreducible as a
polynomial in the $a_{0},a_{1},\ldots,a_{2n}$ for $n>0$).
I dislike this argument for the technicalities involved, including the passage
to an algebraically closed field and the use of the Nullstellensatz (at least for principal ideals). It seems
to me that \eqref{catalecticant1} should have a purely algebraic proof, using
properties of determinants. Sylvester might have had one, but I find his paper
impossible to decipher. Does anyone know such a proof?
 A: Dolgachev (2012, p. 57; pdf) observes that your matrix $\left(  a_{i+j-2}\right)  _{1\leq i\leq n+1,\ 1\leq j\leq n+1}$ (with determinant $\operatorname{Cat} f$) is the matrix of a symmetric bilinear form $\Omega_{\,f}$ on $\smash{\operatorname{Sym}^n(\mathbf k^{2})}$ in a certain basis, then states as obvious that $f\mapsto\Omega_{\,f}$ is a $\smash{\operatorname{GL}(\mathbf k^2)}$-equivariant map 
$$
\operatorname{Sym}^{2n}(\mathbf k^{2})^*\to\operatorname{Sym}^2(\operatorname{Sym}^n(\mathbf k^{2}))^*,
$$
making the first representation a direct summand in the second. (The latter’s full decomposition is in Fulton-Harris, Exercises 6.16, 11.31, 15.45.) That should help...

Edit: Here is a detailed proof. Write $\mathcal F_n=\operatorname{Sym}^n(\mathbf k^{2})^*$. Functoriality gives us morphisms
$$
\begin{array}{ccccc}
\operatorname{End}(\mathbf k^2) & \overset{\operatorname{Sym}^{n*}}\longrightarrow & \operatorname{End}(\mathcal F_n) & \overset{\operatorname{Sym}^2}\longrightarrow & \operatorname{End}(\operatorname{Sym}^2(\mathcal F_n))\\
A & \longmapsto & B_n & \longmapsto & C_n
\end{array}
$$
which satisfy $\det(B_n)= \det(A)^{n(n+1)/2}$ and $C_n(\Omega)=\Omega({}^tB_n\,\cdot,{}^tB_n\,\cdot)=B_n\,\Omega\, {}^tB_n$ if we regard quadratic forms as symmetric matrices; ${}^t$ is transposition $\operatorname{End}(V^*)\leftrightarrow\operatorname{End}(V)$. Thus,
$$
\det(C_n(\Omega))=\det(B_n)^2\det(\Omega)=\det(A)^{n(n+1)}\det(\Omega).
$$
So your $(1)$ will indeed follow 
from the claimed equivariance of $\smash{f\mapsto\Omega_f}$. To see the latter, note that (viewing as you do members of $\mathcal{F}_{2n}$ as functions of $\smash{(\begin{smallmatrix}x\\y\end{smallmatrix})}$) Dolgachev’s definition can be written 
$$
\Omega_{\,f}(\xi,\eta):=\tilde f(\xi\eta),
\quad\text{where}\quad
\tilde f=\tfrac1{(2n)!}\mathrm D^{2n}f(0)
$$
(polar form of $f$) and $\xi\eta\in\operatorname{Sym}^{2n}(\mathbf k^2)$ is the product of $\xi,\eta\in\operatorname{Sym}^n(\mathbf k^2)$. So equivariance boils down to the chain rule $\mathrm D^{2n}(\,f\circ {}^t\!A)(0)
= \mathrm D^{2n}f(0)\circ{}^tB_{2n}$ (e.g. Abraham et al. (1988, Ex. 2.3E)):
$$
\Omega_{\,f\circ\, {}^tA}(\xi,\eta)
= (\,\tilde f\circ{}^tB_{2n})(\xi\eta)
= \tilde f({}^tB_n\xi\,{}^tB_n\eta)
=C_n(\Omega_{\,f})(\xi,\eta).
$$
Remark: Defining coefficients $a_i$ of $f$ as you do, and writing $e_1=(\begin{smallmatrix}1\\0\end{smallmatrix})$, $e_2=(\begin{smallmatrix}0\\1\end{smallmatrix})$, the elementary relation
$
\mathrm D^{2n}(x^iy^{\,j})(0)(e_1^ke_2^l)=(\partial_x^{\,k}\partial_y^{\,l})(x^iy^{\,j})=i!j!\delta_{ik}\delta_{jl}
$
whenever $i+j=k+l=2n$, readily implies that the matrix of $\Omega_f$ indeed has $i,j$-th entry
$
\Omega_f(e_1^{\,n-i}e_2^{\,i},e_1^{\,n-j}e_2^{\,j})=a_{i+j}$ $(i, j = 0,\dots,n).
$
Note also that parts of this are already in Reznick (1992, 1.19, 1.24, 1.30, 1.31), (1996, 2.12).
A: Let $d=2n$ be the degree of your binary form $f$.
Let me introduce $n+1$ pairs of formal variables $\alpha^{(1)}=(\alpha^{(1)}_{1},\alpha^{(1)}_{2}),\ldots, \alpha^{(n+1)}=(\alpha^{(n+1)}_{1},\alpha^{(n+1)}_{2})$.
Let $\mathcal{S}$ be the polynomial
$$
\mathcal{S}=\prod_{1\le i<j\le n+1} (\alpha^{(i)}\alpha^{(j)})^2
$$
where I used the classical bracket notation for $2\times 2$ determinants,
namely,
$$
(\alpha^{(i)}\alpha^{(j)}):= \alpha^{(i)}_{1}\alpha^{(j)}_{2}-\alpha^{(i)}_{2}\alpha^{(j)}_{1}\ .
$$
Let $\mathcal{D}$ be the differential operator
$$
\mathcal{D}=f\left(\frac{\partial}{\partial\alpha^{(1)}_{1}},
\frac{\partial}{\partial\alpha^{(1)}_{2}}\right)\cdots
f\left(\frac{\partial}{\partial\alpha^{(n+1)}_{1}},
\frac{\partial}{\partial\alpha^{(n+1)}_{2}}\right)\ .
$$
Then the catalecticant is up to a numerical factor equal to the result of applying $\mathcal{D}$ to $\mathcal{S}$.
The latter is the symbolic expression of the catalecticant and corresponds to a directed graph on $n+1$ vertices which is basically the complete graph with doubled edges. So each vertex is of degree $d$. As a polynomial $\mathcal{S}$ is multihomogeneous of degree $d$ in each $\alpha^{(i)}$. More generally for any graph of this kind the resulting polynomial $\mathcal{D}\mathcal{S}$ will be an $SL_2$-invariant.
To prove equality with your Hankel determinant, see $\mathcal{S}$ as the square of a (homogenized $\mathbb{A}^1\rightarrow \mathbb{P}^1$) Vandermonde determinant. Take the square of the matrix and then hit its determinant with $\mathcal{D}$.
An alternate construction (essentially the same as the answer by François) is to take the determinant of the map
$$
\begin{array}{ccc}
{\rm Sym}^n(\mathbf{k}^2) & \longrightarrow & {\rm Sym}^n(\mathbf{k}^2)\\
g & \longmapsto & (f,g)_{n}
\end{array}
$$
where $(f,g)_n$ is the $n$-th transvectant, i.e., the equivariant projection
$$
{\rm Sym}^{2n}(\mathbf{k}^2)\otimes {\rm Sym}^{n}(\mathbf{k}^2)\longrightarrow {\rm Sym}^n(\mathbf{k}^2)
$$
You can learn more about symbolic forms and transvectants in Section 2 of my article "On the volume conjecture for classical spin networks". J. Knot Theory Ramifications 21 (2012), no. 3, 1250022.
See in particular Page 14 which contains a graphical proof of the $SL_2$-invariance property. The example is that of the canonisant of the quintic, but you can remove the three $x$ legs and you get the proof for the catalecticant of the binary quartic. Essentially you need the graphical rule for the action of $g\in SL_2$ on a binary form $F$, i.e.,

Then you need the basic identity for determinants:

and finally the rest of the proof:

Again, remove the three $x$'s to do the catalecticant of the quartic instead of the canonisant of the quintic.
