[proofplan]
We prove the theorem degree by degree. In each homogeneous degree $n$, the monomial symmetric functions indexed by partitions of $n$ form a known $\mathbb{Z}$-basis, and the Schur functions expand in that basis through Kostka numbers. Kostka unitriangularity says that, after ordering partitions by any linear extension of dominance order, the Schur-to-monomial transition matrix is unitriangular over $\mathbb{Z}$. Such a matrix is invertible over $\mathbb{Z}$, so the Schur functions form a $\mathbb{Z}$-basis in every degree; the graded direct sum then gives the result for all of $\operatorname{Sym}$.
[/proofplan]
[step:Work in one homogeneous degree]
Fix $n \in \mathbb{N} \cup \{0\}$. Let $\mathcal{P}_n$ denote the finite set of partitions $\lambda$ with $|\lambda| = n$, and let
\begin{align*}
r_n := |\mathcal{P}_n|
\end{align*}
denote its cardinality.
For each $\mu \in \mathcal{P}_n$, let $m_\mu \in \operatorname{Sym}^n$ denote the monomial symmetric function indexed by $\mu$. We use the standard monomial basis theorem for symmetric functions: the family
\begin{align*}
\{m_\mu : \mu \in \mathcal{P}_n\}
\end{align*}
is a $\mathbb{Z}$-basis of $\operatorname{Sym}^n$ (citing a result not yet in the wiki: Monomial Symmetric Function Basis Theorem).
[/step]
[step:Order the partitions by dominance]
Let $\unrhd$ denote the dominance order on $\mathcal{P}_n$: for $\lambda, \mu \in \mathcal{P}_n$, write $\lambda \unrhd \mu$ when
\begin{align*}
\sum_{i=1}^{k} \lambda_i \geq \sum_{i=1}^{k} \mu_i
\end{align*}
for every $k \geq 1$, where missing parts are interpreted as $0$.
Choose a linear extension of this partial order, and write
\begin{align*}
\mathcal{P}_n = \{\lambda_1, \dots, \lambda_{r_n}\}
\end{align*}
so that whenever $\lambda_i \unrhd \lambda_j$, one has $i \leq j$.
[/step]
[step:Write the Schur-to-monomial transition matrix]
For each $\lambda \in \mathcal{P}_n$, the Schur-to-monomial expansion has the form
\begin{align*}
s_\lambda = \sum_{\mu \in \mathcal{P}_n} K_{\lambda\mu} m_\mu,
\end{align*}
where $K_{\lambda\mu} \in \mathbb{Z}_{\geq 0}$ is the Kostka number. We use Kostka unitriangularity: for all $\lambda,\mu \in \mathcal{P}_n$,
\begin{align*}
K_{\lambda\lambda} = 1,
\end{align*}
and
\begin{align*}
K_{\lambda\mu} = 0 \quad \text{unless} \quad \lambda \unrhd \mu.
\end{align*}
This is the standard Schur-to-monomial triangularity theorem (citing a result not yet in the wiki: Kostka Unitriangularity).
Define the integer matrix
\begin{align*}
A_n := (a_{ij})_{1 \leq i,j \leq r_n} \in \mathbb{Z}^{r_n \times r_n}
\end{align*}
by
\begin{align*}
a_{ij} := K_{\lambda_j\lambda_i}.
\end{align*}
Thus the $j$-th column of $A_n$ is the coordinate vector of $s_{\lambda_j}$ in the ordered monomial basis $(m_{\lambda_1}, \dots, m_{\lambda_{r_n}})$.
By Kostka unitriangularity, $a_{jj} = K_{\lambda_j\lambda_j} = 1$. If $a_{ij} \neq 0$, then $K_{\lambda_j\lambda_i} \neq 0$, so $\lambda_j \unrhd \lambda_i$, hence $j \leq i$ by the chosen ordering. Therefore $a_{ij} = 0$ whenever $i < j$, and $A_n$ is lower triangular with diagonal entries all equal to $1$.
[guided]
The goal in this degree is to compare two ordered families in the same free $\mathbb{Z}$-module $\operatorname{Sym}^n$: the known basis $(m_\mu)_{\mu \in \mathcal{P}_n}$ and the candidate basis $(s_\lambda)_{\lambda \in \mathcal{P}_n}$. To do this, we write each Schur function in the monomial basis.
For each partition $\lambda \in \mathcal{P}_n$, the Schur-to-monomial formula is
\begin{align*}
s_\lambda = \sum_{\mu \in \mathcal{P}_n} K_{\lambda\mu} m_\mu,
\end{align*}
where $K_{\lambda\mu} \in \mathbb{Z}_{\geq 0}$ is the Kostka number. The triangularity information is exactly what makes this useful: Kostka unitriangularity says
\begin{align*}
K_{\lambda\lambda} = 1
\end{align*}
and
\begin{align*}
K_{\lambda\mu} = 0 \quad \text{unless} \quad \lambda \unrhd \mu.
\end{align*}
This means that, when partitions are ordered compatibly with dominance, each $s_\lambda$ is equal to $m_\lambda$ plus only monomial terms indexed by partitions lower in dominance order.
We now encode these expansions in a matrix. Define
\begin{align*}
A_n := (a_{ij})_{1 \leq i,j \leq r_n} \in \mathbb{Z}^{r_n \times r_n}
\end{align*}
by
\begin{align*}
a_{ij} := K_{\lambda_j\lambda_i}.
\end{align*}
The reason for this choice of indices is that the $j$-th column records the coefficients of $s_{\lambda_j}$ in the ordered monomial basis:
\begin{align*}
s_{\lambda_j} = \sum_{i=1}^{r_n} a_{ij} m_{\lambda_i}.
\end{align*}
The diagonal entries are
\begin{align*}
a_{jj} = K_{\lambda_j\lambda_j} = 1.
\end{align*}
If $a_{ij} \neq 0$, then $K_{\lambda_j\lambda_i} \neq 0$, so Kostka triangularity gives $\lambda_j \unrhd \lambda_i$. Since our ordering is a linear extension of dominance order, this implies $j \leq i$. Therefore no nonzero entry can occur above the diagonal, which is the condition $i < j$. Hence $A_n$ is lower triangular with every diagonal entry equal to $1$.
[/guided]
[/step]
[step:Invert the transition matrix over $\mathbb{Z}$]
Since $A_n$ is lower triangular with diagonal entries all equal to $1$, its determinant is
\begin{align*}
\det A_n = \prod_{j=1}^{r_n} 1 = 1.
\end{align*}
The [adjugate identity](/theorems/397) gives
\begin{align*}
A_n^{-1} = \operatorname{adj}(A_n),
\end{align*}
because $\det A_n = 1$. Since $A_n$ has integer entries, its adjugate matrix also has integer entries. Therefore $A_n \in GL_{r_n}(\mathbb{Z})$.
It follows that the ordered family
\begin{align*}
(s_{\lambda_1}, \dots, s_{\lambda_{r_n}})
\end{align*}
is obtained from the ordered $\mathbb{Z}$-basis
\begin{align*}
(m_{\lambda_1}, \dots, m_{\lambda_{r_n}})
\end{align*}
by an invertible integer change of basis. Hence
\begin{align*}
\{s_\lambda : \lambda \in \mathcal{P}_n\}
\end{align*}
is a $\mathbb{Z}$-basis of $\operatorname{Sym}^n$.
[/step]
[step:Pass from homogeneous degrees to the whole ring of symmetric functions]
By definition of the grading,
\begin{align*}
\operatorname{Sym} = \bigoplus_{n \geq 0} \operatorname{Sym}^n.
\end{align*}
For each $n \geq 0$, the previous step proves that the Schur functions indexed by partitions of $n$ form a $\mathbb{Z}$-basis of $\operatorname{Sym}^n$. Taking the union of these homogeneous bases over all $n \geq 0$ gives a $\mathbb{Z}$-basis of the direct sum. Therefore
\begin{align*}
\{s_\lambda : \lambda \text{ is a partition}\}
\end{align*}
is a $\mathbb{Z}$-basis of $\operatorname{Sym}$.
[/step]