As a set, i.e. as a von Neumann ordinal, the $\omega$-th limit ordinal $\omega^2$ is fairly complex and not so easy to visualize (for the novice). But as an explicit well-ordering of $\mathbb{N}$, there is a chance, and even more: all limit ordinals less than $\omega^2$ come as good old natural numbers.
Let $\pi:\mathbb{N}\rightarrow \mathbb{N}$ be the function that maps each natural number to its smallest prime factor. Consider the following well-ordering of $\mathbb{N}$:
$$n \preceq m :\equiv \begin{cases} n &\leq m &\mbox{if } &\pi(n) = \pi(m) \\
\pi(n) &< \pi(m) &\mbox{if } &\pi(n) \neq \pi(m) \end{cases}
$$
This ordering is of type $\omega^2$ and it well-orders $\mathbb{N}$ like this:
$$2,4,6,8,\dots,3,9,15,\dots,5,25,35,\dots,7,49,\dots,11,121,\dots$$
Note, that the limit ordinals less than $\omega^2$, i.e. $\omega,\omega\cdot 2,\omega \cdot 3,\dots$ correspond exactly to the (odd) prime numbers.
I wonder if this specific well-ordering of order-type $\omega^2$ is by
any means distinguished – as the most simple or the most natural one –
like the natural ordering of $\mathbb{N}$ is the most natural
well-ordering of type $\omega$.
[Addendum:] For example this well-ordering looks like the limit of a sequence of "natural" well-orderings of type $\omega\cdot k$:
$\omega\cdot 2 = 2,4,6,8,\dots,3,5,7,9,\dots$
the multiples of $2$, followed by the rest
$\omega\cdot 3 = 2,4,6,8,\dots,3,9,15,\dots,5,7,11,13,\dots$
the multiples of $2$, followed by the multiples of $3$ (that are not multiples of $2$), followed by the rest
In comparsion, orderings of type $\omega^2$ based on arbitrary pairing functions (see Joel's answer) come somehow out of the blue.
And I am looking for comparably easy to understand explicit well-orderings of $\mathbb{N}$ of types significantly larger than $\omega^2$, e.g. $\omega^\omega$ or $\epsilon_0$.
Best Answer
The ordinals below $\omega^2$ are exactly those of the form $\omega\cdot n+k$ for natural numbers $n$ and $k$. Thus, these are the ordinals having two digits in base $\omega$, and counting to $\omega^2$ is much like counting to one hundred in this regard.
To order $\mathbb{N}$ in order type $\omega^2$, therefore, is essentially the same as to have a pairing function on the natural numbers, where $\langle n,k\rangle$ is the number that will be assigned to $\omega\cdot n+k$.
There are a huge variety of pairing functions, many of which are simpler than others by various measures.
My favorite is the Cantor pairing function $\langle n,k\rangle = \frac{(n+k)(n+k+1)}2 +k$, which has the advantage that it is a polynomial, easy to compute and invert, and has a simple definition.
Another simple pairing function is $\langle n,k\rangle = 2^n(2k+1)-1$, which has the advantage that it is easily seen to be bijective, and is also easy to compute, invert and define.
The pairing function that is implicit in your order, in contrast, is not so simple on the computational criteria, since it requires us to factor numbers into primes, which can be difficult.
Lastly, let me say that I don't find your question to be mathematically meaningful without a clearer criteria for simplicity or naturality. Shall we minimize the size of the Turing machine that computes the order? Shall we find a smallest defining arithmetic formula? Shall we minimize the logical complexity of the definition? Or what?