We can use the notation of the source document, except I will say we are trying to compute $ab = c$ (for $a,b,c\in\mathbb{F}_p$), so I can use $x$ as a variable when discussing polynomials.
Let $a_1,\dots, a_n$ and $b_1,\dots, b_n$ be $t$-out-of-$n$ secret sharings of $a, b$.
Recall this means that they're computed by taking random $t$-degree polynomials:
$$A(x) = \sum_{i = 0}^t \alpha_i x^i,\quad B(x) = \sum_{i = 0}^t\beta_i x^i$$
Where $A(0) = a, B(0) = b$ (so random subject to this condition).
The shares are then created by evaluating the polynomials on the points $\{1,2,\dots,n\}$.
In particular:
$$\forall i \in\{1,\dots,n\} : a_i = A(i),\quad b_i = B(i)$$
Now we want to compute the product.
The polynomial $C(x) = A(x)B(x)$ has the right constant term (as $C(0) = A(0)B(0) = ab$), but is too high of degree (as you mention).
Moreover the "shares" of $C(i)$ can be computed locally, as $C(i) = A(i)B(i) = a_ib_i$.
But $C(x)$ is degree $2t$ (as you mention). We want to find some other polynomial $\mathcal{C}(x)$ such that:
- $\mathcal{C}(0) = C(0) = ab$
- $\mathcal{C}(x)$ is degree $t$
- Computing $\mathcal{C}(i)$ isn't "too expensive" (in communication) if you already know $A(i)$ and $B(i)$
What will we do then? The idea is to utilize that there are two different ways of representing degree $t$ polynomials:
- Through their $t+1$ coefficients (this is the "obvious") way
- Through the evaulation of them on (at least) $t+1$ distinct points
Either of these is enough information to uniquely reconstruct the polynomial.
The surprising thing is that you can convert from one to the other using a linear operation.
To see how we might establish this, recall that for an $n\times n$ matrix $D$ and vector $\vec{v} = (v_1,\dots, v_n)$, we have that:
$$(D\vec{v})_i = \sum_{k = 0}^{n-1}D_{i, k} v_k
$$
Note that this is similar to the expression:
$$A(x) = \sum_{i = 0}^t \alpha_i x^i$$
If we fix evaluation points $\{1,\dots, n\}$, then we have in fact that:
$$A(i) = \sum_{k = 0}^t \alpha_k i^k$$
This suggests setting $D_{i, k} = i^k$ and $v_k = \alpha_k$.
This is precisely what we'll do, by defining the Vandermonde matrix (with respect to the aforementioned evaluation points):
$$V = \begin{pmatrix}
1^0 & 1^1 & \dots & 1^{n-1}\\
2^0 & 2^1 & \dots & 2^{n-1} \\
\vdots && \ddots & \vdots\\
n^0 & n^1 & \dots & n^{n-1}
\end{pmatrix}$$
Note that:
$$V\begin{pmatrix}\alpha_0\\\vdots\\\alpha_{n-1}\end{pmatrix}
=\begin{pmatrix}
\sum_{k = 0}^{n-1} \alpha_i 1^i\\
\sum_{k = 0}^{n-1} \alpha_i 2^i\\
\vdots\\
\sum_{k = 0}^{n-1} \alpha_i n^i\\
\end{pmatrix} = \begin{pmatrix}A(1)\\
A(2)\\
\vdots\\
A(n)
\end{pmatrix} = \begin{pmatrix}a_1\\
a_2\\
\vdots\\
a_n
\end{pmatrix}
$$
So the Vandermonde Matrix precisely maps the "coefficient representation" of a polynomial to its "evaluation representation".
This ends up being extremely closely related to the Fourier transform. The Discrete Fourier Transform can be written as a Vandermonde Matrix, and the Fast Fourier Transform can be explained that it's Vandermonde Matrix is Toeplitz (and in fact circulant), so admits especially efficient representations and matrix/vector multiplications, but this is an ahistorical aside.
So, we have an (invertible) matrix that maps an $n$-dimensional "coefficient vector" representation of a polynomial to an $n$-dimensional "evaluation vector" representation of a polynomial.
For the moment, don't worry about how all of the people move around all the shares --- just make sure you understand how to do the computation.
We start with the "evaluation vector" representation $C(i) = A(i)B(i)$ for $i\in\{1,\dots,n\}$, which we can write as $\vec c = (c_1,\dots, c_n)$.
We convert this to the "coefficient vector" representation via $V^{-1}\vec{c}$.
This gives the coefficients of the polynomial $C(x) = A(x)B(x)$ as a vector.
While there are $n$ coefficients, as discussed before this polynomial (uniquely determined from $\vec{c}$) is of degree $2t$, so the higher-order coefficients are 0.
We can convert this to a degree $t$ polynomial via truncation.
Let:
$$P = \begin{pmatrix} I_{t+1} & 0\\
0 & 0\end{pmatrix}$$
Be an $n\times n$ block matrix, where $I_{t+1}$ is the $(t+1)\times (t+1)$ identity matrix.
Then $PV^{-1}\vec{c}$ will "drop" higher order coefficients, leaving a degree $t$ polynomial.
Importantly, it doesn't touch the constant term (so $\mathcal{C}(0) = C(0) = ab$ is preserved).
All that's left is to convert back to shares, so to convert from the "coefficient representation" to the "evaluation representation", again by using $V$.
Thus $VPV^{-1}\vec{c}$ will output the shares (of a degree $t$ polynomial) that you want.
Moreover, $VPV^{-1}$ can be precomputed by all participants in the protocol (it's just some $n\times n$ matrix. I could probably even write it out here, but won't).
This reduces the problem of multiplying shares to the problem of "computing a linear equation" of shares, which your source also discusses (at this link).
As this is getting long I'll leave the answer here, but if you don't understand the linear case I encourage you to ask a new question about it.
