Suppose I have
$$ S_n = \sum_{1 \leq j < k \leq n } a_j a_k $$
Goal is to simplify it. This is what I got
$$ S_n = \sum_{j=1}^n \sum_{k = j+1}^n a_j a_k = \sum_{j=1}^n a_j \sum_{k=j+1}^n a_k $$
Is this a correct way to rewrite this sum? Is there a different way to write it?
More exactly $\;\displaystyle S_n = \sum_{j=1}^{n-1} a_j \sum_{k=j+1}^n a_k\quad$ (with $\,n-1\,$ at the upperbound)
You may write too : $\;\displaystyle S_n=\frac{\left(\sum_{k=1}^n a_k\right)^2-\sum_{k=1}^n (a_k)^2}2$
(supposing the product over the $a_k$ commutative of course! :-))
What you write is fine. There are other ways to write it. The condition $1 \leq j < k \leq n$ selects an "area" in the $j$-$k$-plane and any way to count that area is fine. For example, with reverted dependency of the indices, you can make $j$ dependent on $k$:
$S_n = \sum_{k=2}^{n} a_k \sum_{j = 1}^{k-1} a_j$
You can also count diagonals in the $j$-$k$-plane:
$S_n = \sum_{k=2}^{n} \sum_{m = 1}^{n-k+1} a_{k+m-1} a_{m}$
or antidiagonals ....
The schizophrenic way:
You can express this as a single sum of $n(n-1)/2$ terms
$$\sum_{m=1}^{\frac12n(n-1)}a_{n-k(m)}a_{n+m-\frac12k(m)(k(m)+1)},$$
where $k(m):=\left\lfloor\frac{\sqrt{8m-7}+1}2\right\rfloor$.
I offer an explanation of Yves Daoust's answer, as well as some variations.
The number of pairs $(j,k)$ that satisfy $1\leq j < k \leq n$ is $1 + 2 + 3 + \cdots + (n - 1) = \frac12 n(n-1).$ In order to turn the double sum into a single sum, the idea is to let $m$ range over $1, 2, \ldots, \frac12 n(n-1)$ to get the correct number of terms, and for each $m$ somehow obtain $i$ and $j$ as a function of $m$ in such a way that each pair of values $(i,j)$ occurs only once.
Now for $n > 1,$ compare the sum for $1\leq j < k \leq n$ with the sum for $1\leq j < k \leq n - 1.$ The larger sum has exactly $n - 1$ terms that the other sum does not, namely the terms $(1,n), (2,n), \ldots, (n-1,n).$ Assuming (for the moment) that we use the same formulas to get $j$ and $k$ from $m$ in both sums, the "new" terms in the sum for $1\leq j < k \leq n$ must be the last $n-1$ terms of the sum, that is, the terms for $\frac12 n(n-1) - (n-2) \leq m \leq \frac12 n(n-1).$ To make this work, we need a function $f$ such that $f\left(\frac12 n(n-1) - (n-2)\right) = f\left(\frac12 n(n-1)\right) = n$ but $f\left(\frac12 n(n-1) - (n-1)\right) = n - 1$ and $f\left(\frac12 n(n-1) + 1\right) = n + 1.$
The fact that the function has to be constant for a range of input values and then increase by one for the next range of input values suggests that we use some kind of rounding function, such as the "floor" function $\lfloor \cdot \rfloor,$ to force the values of some other increasing function to be zero. If we use the "floor" function for that purpose, it is sufficient to find any function $g$ such that $g\left(\frac12 n(n-1) - (n-2)\right) = n$ for all positive $n,$ and then $f(m) = \lfloor g(m) \rfloor$ will give us the desired value of $k.$
To find such a function, let $m = \frac12 k(k-1) - (k-2),$ so that $g(m) = k,$ and solve for $k.$ For a given value of $m,$ a little algebraic manipulation gives us a quadratic equation in $k$ to solve, $$ k^2 - 3k + 4 - 2m = 0, $$ which has the solutions $$ k = \frac{3 \pm \sqrt{8m - 7}}{2}. $$
It turns out that if we take the $\pm$ on the right-hand side as $+,$ and round the whole thing down to an integer, we get a function $f$ that does what we wanted: $$ f(m) = \left\lfloor \frac{3 + \sqrt{8m - 7}}{2} \right\rfloor. $$
It also turns out that due to the rounding, the derivation of this formula is a bit forgiving; anything in the form $\left\lfloor \frac12 (3 + \sqrt{8m + h}) \right\rfloor$ will do the trick as long as $-7 \leq h < 1.$
At this point we have found out how to get $k = n$ in the last $n-1$ terms of the sum over $m,$ in such a way that the sum has $1$ term with $k=2$ followed by $2$ terms with $k=3,$ then $3$ terms with $k=4,$ and so forth. Now if we can get get $j$ to count up from $1$ to $k-1$ in each of those groups of terms, we are done.
Recognizing that the terms with a particular value of $k$ have $\frac12 k(k-1) - (k-2) \leq m \leq \frac12 k(k-1),$ we can see that setting $$j = m - \frac12 k(k-1) - (k-2) + 1 = m + 2 - \frac12 k(k+1)$$ will produce values of $j$ such that $1 \leq j \leq k - 1.$
Recalling that we decided to set $k = f(m),$ the sum is then $$ \sum_{m=1}^{n(n-1)/2} a_{m + 2 - f(m)(f(m)+1)/2}a_{f(m)}. $$
You may notice that this expression looks quite different from the one in Yves Daoust's answer. That's because we made a number of arbitrary choices in deciding the order in which we would add up the terms. Rather than letting $j$ count up within each sequence of terms with a given value of $k,$ we could have let it count down from $k-1$ to $1.$ Alternatively, we could have taken all the terms for $j=1$ first (there are $n-1$ such terms with $2 \leq k \leq n$), then the terms for $j=2,$ and so forth, ending with a single term for $j=n-1$ (with $k=n$).
Or we could have started with the term for $j=n-1$ and ended with the $n-1$ terms for $j=1.$ That is how Yves Daoust arranged things. The function $$ k(m) = \left\lfloor \frac{1 + \sqrt{8m - 7}}{2} \right\rfloor $$ is equal to $f(m) - 1$ (according to the definition of $f$ above), so instead of taking on values $2, \ldots, n$ for $1 \leq m \leq \frac12 n(n-1)$ it takes on the values $1,\ldots,n-1.$ That is, for $m=1,$ the subscript $n - k(m)$ is $n-1$; for the next two terms ($2 \leq m \leq 3$) the subscript $n - k(m)$ is $n-2$; and so forth. So that's how the subscript in $a_{n-k(m)}$ is chosen.
For the other subscript, we may observe that $\frac12 k(m)(k(m)+1) \geq m,$ with equality when $m=\frac12 p(p-1)$ for some integer $p.$ The quantity $\frac12 k(m)(k(m)+1) - m$ is always in the range $\{0, \ldots, k(m) - 1\},$ so $n + m - \frac12 k(m)(k(m)+1)$ is always in the range $\{n - k(m) + 1, \ldots, n\}.$
A note on the practical use of these formulas:
I have actually used a formula related to these in order to use a uniformly-distributed pseudorandom value to simulate a probability distribution on $\{1,\ldots,n\}$ such that $P(k)$ was proportional to $k.$
A similar technique can reduce the storage space required for a triangular matrix by approximately one half.
