[proofplan]
We use the tableau definition of the Schur function as a generating function for semistandard Young tableaux of fixed shape. First we group the terms in this generating function by content, so the coefficient of each monomial is the number of tableaux with that content. Since Schur functions are symmetric, coefficients are constant on permutation orbits of exponent vectors. Each orbit sum is exactly a monomial symmetric function, and the corresponding coefficient is the Kostka number.
[/proofplan]
[step:Expand the Schur function by tableau content]
Let $n := |\lambda|$. Let $\operatorname{SSYT}(\lambda)$ denote the set of semistandard Young tableaux of shape $\lambda$ with entries in $\mathbb{N}$. For a tableau $T \in \operatorname{SSYT}(\lambda)$, define its content to be the finite-support sequence
\begin{align*}
\operatorname{cont}(T) := (\alpha_1,\alpha_2,\dots),
\end{align*}
where $\alpha_i$ is the number of entries of $T$ equal to $i$. Since $T$ has exactly $n$ boxes, $\sum_{i \geq 1} \alpha_i = n$.
For every finite-support sequence $\alpha = (\alpha_1,\alpha_2,\dots)$ of nonnegative integers, define
\begin{align*}
x^\alpha := \prod_{i \geq 1} x_i^{\alpha_i}.
\end{align*}
By the tableau definition of the Schur function,
\begin{align*}
s_\lambda = \sum_{T \in \operatorname{SSYT}(\lambda)} x^{\operatorname{cont}(T)}.
\end{align*}
Grouping tableaux according to their content gives
\begin{align*}
s_\lambda
= \sum_{\substack{\alpha \in \mathbb{N}_0^{(\mathbb{N})} \\ \sum_i \alpha_i = n}}
\#\{T \in \operatorname{SSYT}(\lambda) : \operatorname{cont}(T)=\alpha\}\,x^\alpha.
\end{align*}
[/step]
[step:Use symmetry to identify equal coefficients on each orbit]
We use the previously proved [symmetry of Schur functions](/theorems/5195) (citing a result not yet in the wiki: Schur functions are symmetric, proved using Bender-Knuth involutions). Since $s_\lambda$ is symmetric, the coefficient of $x^\alpha$ in $s_\lambda$ depends only on the partition obtained by rearranging the nonzero entries of $\alpha$ in weakly decreasing order.
[guided]
The preceding step expresses $s_\lambda$ as a sum of monomials indexed by contents. A content is not necessarily a partition: for example, $(0,2,1,0,\dots)$ records two entries equal to $2$ and one entry equal to $3$. The monomial attached to it is $x_2^2x_3$.
Now use symmetry. Since $s_\lambda$ is a symmetric function, permuting the variables does not change $s_\lambda$. Therefore, if $\alpha$ and $\beta$ are finite-support sequences obtained from one another by permuting coordinates, then the coefficients of $x^\alpha$ and $x^\beta$ in $s_\lambda$ are equal.
Let $\mu$ be the partition obtained by rearranging the nonzero entries of $\alpha$ in weakly decreasing order. The coefficient of every monomial in the orbit of $x^\mu$ is therefore the same. By the grouped expansion above, this common coefficient is
\begin{align*}
\#\{T \in \operatorname{SSYT}(\lambda) : \operatorname{cont}(T)=\mu\}.
\end{align*}
By definition, this number is the Kostka number $K_{\lambda\mu}$.
[/guided]
[/step]
[step:Collect each orbit as a monomial symmetric function]
For a partition $\mu$ of $n$, the monomial symmetric function $m_\mu$ is defined by
\begin{align*}
m_\mu = \sum_{\alpha} x^\alpha,
\end{align*}
where the sum is over all distinct finite-support sequences $\alpha$ obtained by permuting the parts of $\mu$ and adjoining zeros.
From the previous step, every monomial $x^\alpha$ in this orbit occurs in $s_\lambda$ with coefficient $K_{\lambda\mu}$. Hence the total contribution of the orbit corresponding to $\mu$ is
\begin{align*}
K_{\lambda\mu} m_\mu.
\end{align*}
Summing over all partitions $\mu$ with $|\mu|=n=|\lambda|$ gives
\begin{align*}
s_\lambda
= \sum_{\substack{\mu \text{ a partition} \\ |\mu| = |\lambda|}} K_{\lambda\mu} m_\mu.
\end{align*}
This is the desired expansion of $s_\lambda$ in the monomial symmetric basis.
[/step]