Posterior distribution of Normal Normal-inverse-Gamma Conjugacy Here is the setting:
The likelihood of data is
\begin{align}
    p(\boldsymbol{x} | \mu, \sigma^2) 
    &= (\frac{1}{2\pi \sigma^2})^{\frac{n}{2}} exp\left\{ -\frac{1}{2\sigma^2} \sum\limits_{i=1}^n (x_i - \mu)^2 \right\} \nonumber \\
    &= \frac{1}{(2\pi)^{n/2}} (\sigma^2)^{-n/2} exp\left\{ -\frac{1}{2\sigma^2} \left[  \sum\limits_{i=1}^n (x_i - \overline{x})^2 + n(\overline{x} - \mu)^2 \right] \right\}, \nonumber
\end{align}
and we use the Normal-inverse-Gamma as prior
\begin{align}
    p(\mu , \sigma^2) 
    &= \mathcal{N} (\mu | \mu_0 , \sigma^2 V_0) IG(\sigma^2 | \alpha_0 , b_0 ) \nonumber \\
    &= \frac{1}{\sqrt{2\pi V_0}} \frac{b_0^{\alpha_0}}{\Gamma(\alpha_0)}\frac{1}{\sigma} (\sigma^2)^{-\alpha_0 - 1} exp\left( -\frac{1}{2\sigma^2} [V_0^{-1}(\mu - \mu_0 )^2 + 2b_0] \right). \nonumber
\end{align}
Then, the posterior can be derivated via
\begin{align}
    p(\mu , \sigma^2 | \boldsymbol{x})
    &\propto p(\boldsymbol{x} | \mu , \sigma^2 ) p(\mu , \sigma^2) \nonumber \\
    &\propto \left[ (\sigma^2)^{-n/2} exp \left( -\frac{1}{2\sigma^2} \big[\sum\limits_{i=1}^b (x_i - \overline{x})^2 + n(\overline{x} - \mu )^2\big] \right) \right]  \times \left[ \sigma^{-1} (\sigma^2)^{-\alpha_0 - 1} exp \left( -\frac{1}{2\sigma^2} \left[ V_0^{-1} (\mu - \mu_0 )^2 + 2 b_0 \right] \right) \right] \nonumber \\
    &= \sigma^{-1} (\sigma^2)^{-(\alpha_0 + \frac{n}{2}) - 1} exp \left( -\frac{1}{2\sigma^2} \big[ V_0^{-1} (\mu - m_0 )^2 + 2 b_0 + \sum\limits_{i=1}^n (x_i - \overline{x})^2 + n(\overline{x} - \mu)^2 \big] \right) \nonumber \\
    &= \sigma^{-1} (\sigma^2)^{-(\alpha_0 + \frac{n}{2}) - 1} exp \Big\{ -\frac{1}{2\sigma^2} \Big[ (V_0^{-1} + n)(\mu - \frac{V_0^{-1} m_0 + n\overline{x}}{V_0^{-1} + n})^2 + \big(b_0 + \frac{1}{2}\sum\limits_{i=1}^n (x_i - \overline{x})^2 + \frac{V_0^{-1} n}{2(V_0^{-1} + n)} (m_0 - \overline{x})^2 \big) \Big] \Big\} \nonumber
\end{align}
We recognize this is an unnormalized Normal-inverse-Gamma distribution, therefore
\begin{align}
    p(\mu , \sigma^2 | \boldsymbol{x}) = NIG(\mu , \sigma^2 | m_n , V_n , \alpha_n , b_n ), \nonumber
\end{align}
where
\begin{align}
    m_n &= \frac{V_0^{-1} m_0 + n \overline{x}}{V_0^{-1} + n} \nonumber \\
    V_n^{-1} &= V_0^{-1} + n \nonumber \\
    \alpha_n &= \alpha_0 + \frac{n}{2} \nonumber \\
    b_n &= b_0 + \frac{1}{2}\sum\limits_{i=1}^n (x_i - \overline{x})^2 + \frac{V_0^{-1} n}{2(V_0^{-1} + n)}(m_0 - \overline{x})^2. \nonumber
\end{align}
As indicated in this paper (see Eq(200)), the last term can be further expressed as
\begin{align}
    b_n &= b_0 + \frac{1}{2} \left[ m_0^2 V_0^{-1} + \sum\limits_{i=1}^n x_i^2 - m_n^2 V_n^{-1} \right]. \nonumber 
\end{align}
But I fail to prove it, i.e.,
\begin{align}
    \sum\limits_{i=1}^n (x_i - \overline{x})^2 + \frac{V_0^{-1} n}{(V_0^{-1} + n)}(m_0 - \overline{x})^2 &= \left[ m_0^2 V_0^{-1} + \sum\limits_{i=1}^n x_i^2 - m_n^2 V_n^{-1} \right]. \nonumber 
\end{align}

This is much simpler to prove compared with your earlier question.
\begin{align}
\sum\limits_{i=1}^n (x_i - \overline{x})^2 &+ \frac{V_0^{-1} n}{(V_0^{-1} + n)}(m_0 - \overline{x})^2 = \sum\limits_{i=1}^n x_i^2 - n\overline{x}^2\\
&\qquad + \frac{V_0^{-1} n}{(V_0^{-1} + n)}(m_0^2 - 2 m_0\overline{x} +  \overline{x}^2)\\
&= \sum\limits_{i=1}^n x_i^2 - n\overline{x}^2 + \frac{V_0^{-1} (n+V_0^{-1}-V_0^{-1})}{(V_0^{-1} + n)}m_0^2\\
&\quad-2\frac{V_0^{-1} nm_0\overline{x}}{(V_0^{-1} + n)}
+\frac{(V_0^{-1}+n-n) n}{(V_0^{-1} + n)}\overline{x}^2\\
&= \sum\limits_{i=1}^n x_i^2 + V_0^{-1}m_0^2 - \frac{V_0^{-2}m_0^2}{(V_0^{-1} + n)}\\
&\quad -2\frac{V_0^{-1} m_0n\overline{x}}{(V_0^{-1} + n)}
-\frac{n^2\overline{x}^2}{(V_0^{-1} + n)}\\
&= V_0^{-1}m_0^2 + \sum\limits_{i=1}^n x_i^2 - \frac{(n\overline x+V_0^{-1} m_0)^2}{V_n^{-1}}\\
&= V_0^{-1}m_0^2 + \sum\limits_{i=1}^n x_i^2 -V_n^{-1}m_n^2
\end{align}
