Your results do not appear correct. This is easy to see, without any calculation, because in your table, your $E[X_{(1)}]$ increases with sample size $n$; plainly, the expected value of the sample minimum must get smaller (i.e. become more negative) as the sample size $n$ gets larger.
The problem is conceptually quite easy.
In brief: if $X$ ~ $N(0,1)$ with pdf $f(x)$:
... then the pdf of the 1st order statistic (in a sample of size $n$) is:
... obtained here using the OrderStat function in mathStatica, with domain of support:
Then, $E[X_{(1)}]$, for $n = 1,2,3$ can be easily obtained exactly as:
The exact $n = 3$ case is approximately $-0.846284$, which is obviously different to your workings of -1.06 (line 1 of your Table), so it seems clear something is wrong with your workings (or perhaps my understanding of what you are seeking).
For $n \ge 4$, obtaining closed-form solutions is more tricky, but even if symbolic integration proves difficult, we can always use numerical integration (to arbitrary precision if desired). This is really very easy ... here, for instance, is $E[X_{(1)}]$, for sample size $n = 1$ to 14, using Mathematica:
sol = Table[NIntegrate[x g, {x, -Infinity, Infinity}], {n, 1, 14}]
{0., -0.56419, -0.846284, -1.02938, -1.16296, -1.26721, -1.35218, -1.4236, -1.48501,
-1.53875, -1.58644, -1.62923, -1.66799, -1.70338}
All done. These values are obviously very different to those in your table (right hand column).
To consider the more general case of a $N(\mu, \sigma^2)$ parent, proceed exactly as above, starting with the general Normal pdf.
The calculation of $E\left[X^{-1}\right]$ when $X$ is a Gamma random variable with order parameter $n$ and rate parameter $\lambda$
requires recognition of the density of another Gamma random variable (with order
parameter $n-1$ and rate parameter $\lambda$) in the
integral given by the law of the unconscious statistician for
$E\left[X^{-1}\right]$. We have
$$\begin{align}
E\left[X^{-1}\right]&=
\int_0^\infty \frac 1x \cdot \underbrace{\lambda \frac{(\lambda x)^{n-1}}{\Gamma(n)}e^{-\lambda x}}_{\Gamma(n,\lambda)~\text{density}}\,\mathrm dx\\
&= \lambda\frac{\Gamma(n-1)}{\Gamma(n)}\int_0^\infty \underbrace{\lambda \frac{(\lambda x)^{n-2}}{\Gamma(n-1)}e^{-\lambda x}}_{\Gamma(n-1,\lambda)~\text{density}}\,\mathrm dx\\
&= \frac{\lambda}{n-1}
\end{align}$$
since for positive integer $k$, $\Gamma(k) = (k-1)!$.
Best Answer
Notice first
\begin{eqnarray} \Gamma(1/2) &=& \int_{0}^{\infty}y^{-1/2}e^{-y}dy \\ &=& \int_{0}^{\infty} \sqrt{2}z^{-1} e^{-z^2/2} z \ d z \qquad \text{ (substitute $y=\frac{z^2}{2}$ )} \\ &=& \frac{1}{2} \int_{-\infty}^{\infty} \sqrt{2} e^{-z^2/2} \ d z \qquad \text{ (since even function)} \\ &=& \int_{-\infty}^{\infty} \frac{1}{\sqrt{2}} z^2 e^{-z^2/2} \ d z \qquad \text{ (*using integration by parts)} \\ \end{eqnarray}
The last line (integration by parts) is valid as
\begin{eqnarray} \int_{-\infty}^{+\infty}z^{2}\left(e^{-z^2/2}\right)dz &=&\int_{-\infty}^{+\infty}z\left(ze^{-z^2/2}\right)dz \\ &=& \int_{-\infty}^{+\infty}z\left(-e^{-z^2/2}\right)'dz \\ &=& \underbrace{-ze^{-z^2/2}\Bigg|^{+\infty}_{-\infty}}_{=0}+\int_{-\infty}^{+\infty}\left(e^{-z^2/2}\right)dz\\ \end{eqnarray}
Hence we can easily see
$$\Gamma(1/2) = \frac{1}{\sqrt{\pi}} \mathbb{E}[Z^2] $$
hence
$$\Gamma(1/2) = \sqrt{\pi} \text{ □ }$$
Addition: To be clear, below I am answering exactly the question user @Sarina asked. My answer is based on a result he assumes (that the variance of the standard normal equals 1). The answer would be more complete by proving that this is the case (that involves a simple change from Cartesian to polar coordinates) but I feel this is not what user @Sarina wanted to know