The expectation of sum of squares of the arrival times in a Poisson process Question:

Let $\{N(t)\}$ be rate $\lambda$ Poisson process, with arrival times
$\{S_n,n=0,1,...\}$. Evaluate the expected sum of squares of the
arrival times occuring before $t$,
$$E(t)=\mathbb{E}\left[\sum_{n=1}^{N(t)}S_n^2\right]$$
where we define $\sum_{n=1}^{0} S_n=0$.

There's a similar question here. However, in this question, we are looking for the expected sum of squares. I tried to understand the answer to the similar question in order to solve the problem on my own. However, I cannot even understand that. I am confused. I do not even know where to start. Is there a formula for $S_i$'s?
 A: Let us compute the posterior of $S_i/t|N(t)$ with $i$ an integer between $1$ and $N(t)$ (as mentioned in the other post, it is a beta distribution).
The three parts are $$\begin{split}N(t)-i|S_i&\sim\text{Poisson}\left(\lambda(t-S_i)\right)\\
S_i&\sim \text{Gamma}\left(i, \lambda\right)\\
N(t)&\sim \text{Poisson}\left(\lambda t\right)\end{split}$$
$$\begin{split}f(S_i=s|N(t)=n)&=\frac{f(N(t)=n|S_i=s)f(S_i=s)}{f(N(t)=n)}\\
&=\frac{\frac{e^{-\lambda(t-s)}(\lambda(t-s))^{n-i}}{(n-i)!}\frac{\lambda^i}{\Gamma(i)}s^{i-1}e^{-\lambda s}}{\frac{e^{-\lambda t}(\lambda t)^n}{n!}}\\
&=\frac{e^{-\lambda(t-s)}\lambda^{n-i}(t-s)^{n-i}\lambda^is^{i-1}e^{-\lambda s}n!}{(n-i)!\Gamma(i)e^{-\lambda t}(\lambda t)^n}\\
&=\frac{(t-s)^{n-i}s^{i-1}n!}{(n-i)!\Gamma(i)t^n}\\
&=\frac{\Gamma(n+1)}{\Gamma(n-i+1)\Gamma(i)}\frac{(t-s)^{n-i}s^{i-1}}{t^n}\end{split}$$
Now consider the random variable $Y=S_i/t|N(t)$. By the transformation theorem it has density
$$\begin{split}
f_Y(y)&=f_{S_i}(yt)|t|\\
&=\frac{\Gamma(N(t)+1)}{\Gamma(N(t)-i+1)\Gamma(i)}
(t-yt)^{N(t)-i} (yt)^{i-1} t^{-N(t)} t
\\ 
&= \frac{\Gamma(N(t)+1)}{\Gamma(N(t)-i+1)\Gamma(i)}y^{i-1}(1-y)^{N(t)-i}, 0<y<1
\end{split}$$
which is Beta distribution with parameters $\alpha = i, \beta = N(t)-i+1$.
We seek $$\begin{split}E\left(\sum_{i=1}^{N(t)} S_i^2\right)&=E\left(\sum_{i=1}^{N(t)} (S_i/t)^2t^2\right)\\
&=t^2E\left(E\left(\sum_{i=1}^{N(t)} (S_i/t)^2\bigg|N(t) \right)\right)\\
&=t^2E\left(\sum_{i=1}^{N(t)} E((S_i/t)^2|N(t)) \right)\\
&=t^2E\left(\sum_{i=1}^{N(t)} \frac{i(i+1)}{(N(t)+1)(N(t)+2)}\right)
\end{split}$$
The last equality follows from the expected value of $(S_i/t)^2|N(t)$ being the second moment which is looked up to be $\frac{i(i+1)}{(N(t)+1)(N(t)+2)}$. The
$E[X^2]$ which $X$ is beta distribution is completely discussed here.
Therefore we have: $$
\begin{split}
E\left(\sum_{i=1}^{N(t)} S_i^2\right) &=
t^2E\left(\frac{1}{(N(t)+1)(N(t)+2)} \left[ \sum_{i=1}^{N(t)}i^2+\sum_{i=1}^{N(t)} i \right] \right)\\
&= t^2E\left(\frac{\frac{N(t)(N(t)+1)(2N(t)+1)}{6}+\frac{N(t)(N(t)+1)}{2}}{(N(t)+1)(N(t)+2)} \right)\\
&=t^2E\left( \frac{2N(t)^2+4N(t)}{6(N(t)+2)} \right) \\
&= t^2 E\left( \frac{N(t)}{3}\right)\\
&=\frac{\lambda t^3}{3}
\end{split}$$
