[proofplan]
We verify the algebra axioms directly from the entry formula for matrix multiplication. First we note that the product of two $n\times n$ matrices is again an $n\times n$ matrix over $k$. Then bilinearity, associativity, and the identity laws are proved by comparing individual entries and using the field axioms in $k$. These entrywise identities are exactly the defining properties of a unital associative $k$-algebra.
[/proofplan]
[step:Define the multiplication map and prove closure in $M_n(k)$]
Let
\begin{align*}
\mu:M_n(k)\times M_n(k)&\to M_n(k)
\end{align*}
denote the map defined by $\mu(A,B)=AB$, where for $A,B\in M_n(k)$ and $1\le i,j\le n$,
\begin{align*}
(AB)_{ij}=\sum_{\ell=1}^n A_{i\ell}B_{\ell j}.
\end{align*}
For each fixed pair $(i,j)$, every factor $A_{i\ell}$ and $B_{\ell j}$ lies in $k$, and since $k$ is closed under multiplication and finite addition, the finite sum $\sum_{\ell=1}^n A_{i\ell}B_{\ell j}$ lies in $k$. Thus $AB$ is an $n\times n$ matrix with entries in $k$, so $\mu$ is a well-defined map from $M_n(k)\times M_n(k)$ to $M_n(k)$.
[/step]
[step:Verify linearity in each matrix variable]
Let $A,A',B,B'\in M_n(k)$ and let $\lambda\in k$. For every $1\le i,j\le n$, distributivity in $k$ gives
\begin{align*}
(((A+A')B)_{ij})=\sum_{\ell=1}^n (A_{i\ell}+A'_{i\ell})B_{\ell j}.
\end{align*}
Hence
\begin{align*}
((A+A')B)_{ij}=\sum_{\ell=1}^n A_{i\ell}B_{\ell j}+\sum_{\ell=1}^n A'_{i\ell}B_{\ell j}.
\end{align*}
Therefore
\begin{align*}
((A+A')B)_{ij}=(AB)_{ij}+(A'B)_{ij}.
\end{align*}
Since this holds for every entry, $(A+A')B=AB+A'B$.
Similarly, for scalar multiplication in the first variable,
\begin{align*}
((\lambda A)B)_{ij}=\sum_{\ell=1}^n (\lambda A_{i\ell})B_{\ell j}.
\end{align*}
Using associativity and distributivity in $k$,
\begin{align*}
((\lambda A)B)_{ij}=\lambda\sum_{\ell=1}^n A_{i\ell}B_{\ell j}.
\end{align*}
Thus $((\lambda A)B)_{ij}=(\lambda AB)_{ij}$ for every $i,j$, so $(\lambda A)B=\lambda(AB)$.
For the second variable, the same entrywise calculation gives
\begin{align*}
(A(B+B'))_{ij}=\sum_{\ell=1}^n A_{i\ell}(B_{\ell j}+B'_{\ell j}).
\end{align*}
Hence
\begin{align*}
(A(B+B'))_{ij}=(AB)_{ij}+(AB')_{ij}.
\end{align*}
Thus $A(B+B')=AB+AB'$. Finally,
\begin{align*}
(A(\lambda B))_{ij}=\sum_{\ell=1}^n A_{i\ell}(\lambda B_{\ell j}).
\end{align*}
Since multiplication in a field is commutative and associative,
\begin{align*}
(A(\lambda B))_{ij}=\lambda\sum_{\ell=1}^n A_{i\ell}B_{\ell j}.
\end{align*}
Therefore $A(\lambda B)=\lambda(AB)$. These four identities prove that $\mu$ is $k$-bilinear.
[/step]
[step:Compare entries to prove associativity]
Let $A,B,C\in M_n(k)$. We prove $(AB)C=A(BC)$ by comparing the $(i,j)$ entries for arbitrary $1\le i,j\le n$.
By the definition of matrix multiplication,
\begin{align*}
(((AB)C)_{ij})=\sum_{r=1}^n (AB)_{ir}C_{rj}.
\end{align*}
Expanding $(AB)_{ir}$ gives
\begin{align*}
((AB)C)_{ij}=\sum_{r=1}^n \left(\sum_{\ell=1}^n A_{i\ell}B_{\ell r}\right)C_{rj}.
\end{align*}
Using distributivity and associativity in $k$ for finite sums,
\begin{align*}
((AB)C)_{ij}=\sum_{r=1}^n\sum_{\ell=1}^n A_{i\ell}B_{\ell r}C_{rj}.
\end{align*}
Because the sums are finite, we may interchange the order of summation:
\begin{align*}
((AB)C)_{ij}=\sum_{\ell=1}^n\sum_{r=1}^n A_{i\ell}B_{\ell r}C_{rj}.
\end{align*}
Factoring out $A_{i\ell}$ from the inner finite sum gives
\begin{align*}
((AB)C)_{ij}=\sum_{\ell=1}^n A_{i\ell}\left(\sum_{r=1}^n B_{\ell r}C_{rj}\right).
\end{align*}
The inner sum is $(BC)_{\ell j}$, so
\begin{align*}
((AB)C)_{ij}=\sum_{\ell=1}^n A_{i\ell}(BC)_{\ell j}.
\end{align*}
By the definition of matrix multiplication again,
\begin{align*}
((AB)C)_{ij}=(A(BC))_{ij}.
\end{align*}
Since this equality holds for every $1\le i,j\le n$, we have $(AB)C=A(BC)$.
[guided]
The associativity statement is an equality of matrices, so it is enough to prove equality of every entry. Fix indices $i,j\in\{1,\dots,n\}$. The left-associated product first forms $AB$ and then multiplies by $C$, so the entry formula gives
\begin{align*}
((AB)C)_{ij}=\sum_{r=1}^n (AB)_{ir}C_{rj}.
\end{align*}
Now expand the entry $(AB)_{ir}$ using the same definition:
\begin{align*}
((AB)C)_{ij}=\sum_{r=1}^n \left(\sum_{\ell=1}^n A_{i\ell}B_{\ell r}\right)C_{rj}.
\end{align*}
The expression involves only finite sums in the field $k$. Therefore distributivity lets us multiply $C_{rj}$ through the inner sum, and associativity in $k$ lets us remove ambiguity in the triple product:
\begin{align*}
((AB)C)_{ij}=\sum_{r=1}^n\sum_{\ell=1}^n A_{i\ell}B_{\ell r}C_{rj}.
\end{align*}
Because the index set $\{1,\dots,n\}\times\{1,\dots,n\}$ is finite, changing the order in which the pairs $(\ell,r)$ are summed does not change the result:
\begin{align*}
((AB)C)_{ij}=\sum_{\ell=1}^n\sum_{r=1}^n A_{i\ell}B_{\ell r}C_{rj}.
\end{align*}
Now the reason for this rearrangement becomes visible: for fixed $\ell$, the inner sum is almost the $(\ell,j)$ entry of $BC$. Factoring the scalar $A_{i\ell}$ out of the finite sum gives
\begin{align*}
((AB)C)_{ij}=\sum_{\ell=1}^n A_{i\ell}\left(\sum_{r=1}^n B_{\ell r}C_{rj}\right).
\end{align*}
By definition,
\begin{align*}
(BC)_{\ell j}=\sum_{r=1}^n B_{\ell r}C_{rj}.
\end{align*}
Substituting this into the previous expression yields
\begin{align*}
((AB)C)_{ij}=\sum_{\ell=1}^n A_{i\ell}(BC)_{\ell j}.
\end{align*}
The right-hand side is exactly the entry formula for $(A(BC))_{ij}$:
\begin{align*}
((AB)C)_{ij}=(A(BC))_{ij}.
\end{align*}
Since the indices $i$ and $j$ were arbitrary, every entry of $(AB)C$ equals the corresponding entry of $A(BC)$. Hence $(AB)C=A(BC)$.
[/guided]
[/step]
[step:Use the Kronecker delta entries of $I_n$ to prove the identity laws]
Let $I_n\in M_n(k)$ be the identity matrix, whose entries are defined by
\begin{align*}
(I_n)_{ij}=\delta_{ij},
\end{align*}
where $\delta_{ij}=1$ if $i=j$ and $\delta_{ij}=0$ if $i\ne j$.
Let $A\in M_n(k)$. For every $1\le i,j\le n$,
\begin{align*}
(I_nA)_{ij}=\sum_{\ell=1}^n (I_n)_{i\ell}A_{\ell j}.
\end{align*}
Since $(I_n)_{i\ell}=\delta_{i\ell}$, all terms with $\ell\ne i$ are zero and the term with $\ell=i$ is $A_{ij}$. Thus
\begin{align*}
(I_nA)_{ij}=A_{ij}.
\end{align*}
Therefore $I_nA=A$.
Similarly,
\begin{align*}
(AI_n)_{ij}=\sum_{\ell=1}^n A_{i\ell}(I_n)_{\ell j}.
\end{align*}
Since $(I_n)_{\ell j}=\delta_{\ell j}$, all terms with $\ell\ne j$ are zero and the term with $\ell=j$ is $A_{ij}$. Hence
\begin{align*}
(AI_n)_{ij}=A_{ij}.
\end{align*}
Therefore $AI_n=A$. Thus $I_n$ is a two-sided multiplicative identity for matrix multiplication on $M_n(k)$.
[/step]
[step:Conclude the unital associative algebra structure]
The space $M_n(k)$ is a $k$-[vector space](/page/Vector%20Space) under entrywise addition and scalar multiplication. The multiplication map $\mu:M_n(k)\times M_n(k)\to M_n(k)$ has been shown to be $k$-bilinear, associative, and to have $I_n$ as a two-sided identity. These are precisely the defining conditions for $M_n(k)$, equipped with ordinary matrix multiplication, to be a unital associative algebra over $k$.
[/step]
