This algebra is defined on the permutation module of the symmetric group. It was studied by K. Harada and R. Griess in the 1970s and a proof that its automorphism group is the symmetric group can also be found in Haraada's paper (see also this one) and the appendix of this paper by Dong and Griess, which explains some connections with vertex operator algebras. See also remarks in Griess's paper here.
In the 1970s and 1980s there were many constructions of algebras of this spirit, in the search for model algebraic structures acted on in some canonical way by finite simple groups (I have in mind, for example, papers of A. Ryba, for example there is one constructing a commutative nonassociative algebra having the Harada-Norton group as automorphisms). My (inexpert) sense is that this literature is more or less subsumed into that about vertex operator algebras. More precisely, for vertex operator algebras whose $0$th graded piece is $1$-dimensional and $1$st graded piece is $0$-dimensional there is a commutative nonassociative algebra on the $2$nd graded piece that has many features in common with $E_{n}$ (here I have in mind the paper of Miyamoto on Griess algebras). My understanding of vertex operator algebras is inadequate to say more (and there are some experts who post here who can maybe say something clearer).
This algebra, call it $E_{n}$, has other interesting properties. Let $\tau$ be the Killing type form defined by $\tau(x, y) = \operatorname{tr} L(x)L(y)$ where $L(x)y = xy$ is the multiplication operator. Note that $\operatorname{tr} L(x) =0$ for all $x$ (I call such algebras exact). The symmetric bilinear form $\tau$ is nondegenerate and invariant in the sense that $\tau(xy, z) = \tau(x, yz)$ (this is somehow analogous to the properties of the Killing form of a semisimple Lie algebra). The algebra $E_{n}$ is determined up to isometric automorphism by the associated harmonic cubic polynomial $\tau(xx, x)$. This algebra is conformally associative in the sense that its associator $[x, y, z]= (xy)z - x(yz)$ satisfies $$ \tau([x, y, z], w) =c\left(\tau(x, w)\tau(y,z) - \tau(x,z)\tau(y, w)\right) $$ where $c$ is the appropriate (dimension dependent) constant that I am too lazy to calculate right now. One can prove that up to isometric automorphism there is a unique $n$-dimensional commutative nonassociative algebra that thas the properties that a. $\operatorname{tr} L(x) =0$ for all x, b. its Killing type trace form $\tau$ is nondegenerate and invariant and c. it is conformally associative. (If you contact me directly I can provide you a preprint of mine where this is proved as part of a larger project.)
Another place these algebras $E_{n}$ occur is the following. Consider the Jordan algebra of real symmetric $n\times n$ matrices. Consider its deunitalization that comprises trace-free symmetric matrices with the product $$x \star y = \tfrac{1}{2}\left(xy + yx - \tfrac{1}{n}\operatorname{tr}(xy)I\right). $$ This is a commutative nonassociative algebra such that $\operatorname{tr}L(x) =0$ for all $x$ and the Killing type trace-form $\tau(x, y) = \operatorname{tr}L(x)L(y)$ is invariant (in fact it is a constant multiple of $\operatorname{tr}(xy)$). Although this algebra is far from being conformally nonassociative in the above sense, its subalgebra of (trace-free) diagonal matrices is isomorphic to the algebra $E_{n-1}$. Moreover, every every idempotent element of the deunitalized Jordan algebra is in the orbit of an idempotent of $E_{n-1}$. This can be used, for example, to prove that the deunitalized Jordan algebras are simple. (All this is written in my aforementioned preprint; contact me if you would like details).