Is an algorithm known for ordering natural numbers in the $\epsilon_0$ ordinal?

Question

I have read that there are countable ordinals which can not be reached using algebraic operations and limits starting from $1$, but $\epsilon_0 = \omega \uparrow\uparrow \omega$ is certainly not amongst them. While I can envisage well-orderings like $\omega^2$, $\omega^3$, even $\omega^\omega$ applied to the underlying set of natural numbers, I don't know where to even start to reach $\epsilon_0$. I tried applying the same logic that got me to $\omega^\omega$ (which involved writing $n$ using prime decomposition and sorting lexicographically on powers of successive primes) recursively, but I think the highest I could reach using this technique was $\omega^{\omega^\omega}$ (or possibly $\omega^{\omega^2}$, I would have to redo the derivation as it's been many years ago).

Is any algorithm on how to order $\mathbb{N}$ so that the ordering is order isomorphic to $\epsilon_0$ known, or is it known that it can not be found? If so, what is the largest order for which it is possible?

Update

Here's an explicit construction based on Mitchell Spector's answer and eliminating the need for enumeration and elimination. In other words, an explicit bijection between $\mathbb{N}_0$ and ordinals smaller than $\epsilon_0$.

Taking an ordinal $\alpha$, write it in the Cantor normal form, expressing the exponents likewise recursively. For example,

$$α = ω^{ω^2+ω.2}+ω^2.3+ω+2.$$

Construct a $n = f(α) \in \mathbb{N}$ according to the following rules:

1. $f(k\omega^β) = p^k_{f(\beta)}, \forall k\in\mathbb{N}_0, \forall \beta\in\mathord{\rm Ord}, \beta < \epsilon_0,$ where $p_j$ denotes the $j$-th prime ($p_1$ = 2),
2. $f(β+γ) = f(β)f(γ)$.

We can show this is a bijection by showing how $f^{-1}$ can be obtained:

Taking an $n \in \mathbb{N}$, write down the prime factorization. For every prime that appears in it, write it as $p_j$ for some $j\in\mathbb{N}$. Repeat the factorization in the indices until no numbers other than $1$ appear in the chains of $p_\circ$. For example,

$$1000 = 2^3 5^3 = p_1^3 p_3^3 = p_1^3 p_{p_2}^3 = p_1^3 p_{p_{p_1}}^3.$$

Create an ordinal $f^{-1}(n)$ from this decomposition according to the rules

1. $f^{-1}(kl) = f^{-1}(k) + f^{-1}(l)$, taking the larger of the ordinals first,
2. $f^{-1}(p_k) = \omega^{f^{-1}(k)}$.

This will produce the Cantor normal form (with multiples by an integer represented by sums of a number of equal ordinals). Due to the properties of factorization it is possible for any $n$ and has a unique result for each. The common language between the two worlds is a tree structure (the very same that's used to represent a state of the Hydra game) which can uniquely be found for both ordinals $<ε_0$ via the CNF and natural numbers via the "index-inherited" factorization.

Here's the mapping for a few first $n \in \mathbb{N}$: $$\begin{aligned} 1 &&&≘ 0, \\ 2 &= p_1^1 &&≘ 1, \\ 3 &= p_2 = p_{p_1} &&≘ ω, \\ 4 &= p_1^2 &&≘ 2, \\ 5 &= p_3 = p_{p_{p_1}} &&≘ ω^ω, \\ 6 &= p_1 p_{p_1} &&≘ ω + 1, \\ 7 &= p_4 = p_{p_1^2} &&≘ ω^2, \\ 8 &= p_1^3 &&≘ 3, \\ 9 &= p_2^2 = p_{p_1}^2 &&≘ ω.2, \ldots \end{aligned}$$

In an ordering pulled back from $\rm Ord$ onto $\mathbb{N}$ via $f^{-1}$, one can observe the following:

1. The successor of $n$ is $2n$,
2. If $k\perp l$, the ($k.l$)-th prime is strictly larger any positive power of $p_k$, just as those of $p_{k^2}$, $p_{k^3}$, etc., or products thereof.

There may be some glitch in the last line, but to give an idea,

$$\begin{aligned} 1\ [0] &≤ 2\ [1] ≤ 4\ [2] ≤ 8\ ≤ 16\ ≤ \ldots \\ &\quad≤ 3\ [ω] ≤ 6\ [ω+1] ≤ 12\ ≤ \ldots \\ &\quad≤ 9\ [ω.2] ≤ 18\ ≤ 36\ ≤ \ldots \\ &\quad≤ 27\ [ω.3] ≤ \ldots \\ &\quad≤ \ldots \\ &≤ 7 [ω^2] ≤ \ldots \ldots ≤ 5\ [ω^ω] ≤ \ldots \ldots ≤ 11\ [ω^{ω^ω}] ≤ \ldots\ldots ≤ 31\ [ω^{ω^{ω^ω}}] ≤ \ldots\ldots \end{aligned}$$

Answer 1

Using Cantor normal form, every ordinal less than $\varepsilon_0$ can be written as a finite arithmetic expression starting with $1$ and applying the unary operation of exponentiation with base $\omega$ and the binary operation of addition. (The ordinal $0$ is represented by the empty expression.)

Enumerate all the expressions of the above form (they're just strings on a finite alphabet): $e_0, e_1, e_2, \dots.$

Eliminate from the list any $e_n$ whose value (as an ordinal) is equal to the value (as an ordinal) of some $e_k$ for $k\lt n,$ leaving you with $e'_0, e'_1, e'_2, \dots.$

Then the relation $R$ on $\mathbb{N}$ given by $mRn$ iff the ordinal value of $e'_m$ is less than the ordinal value of $e'_n$ is a well-ordering of $\mathbb{N}$ of order type $\varepsilon_0.$

---

If you prefer, you can note that every ordinal less than $\varepsilon_0$ can be uniquely written in Cantor normal form, with each exponent written in Cantor normal form, each exponent in those Cantor normal forms written in Cantor normal form, etc. (at all superscript levels, similarly to what one sees in Goodstein's theorem). Then you can simply enumerate these expressions and skip the step of passing from $e_n$ to $e'_n,$ because, with this approach, there aren't any duplicates.

This method also has the advantage that it's easier to see how you can look at two expressions $s$ and $t$ and tell if the value represented by $s$ is less than the value represented by $t,$ without knowing much about ordinal computations.

---

The largest ordinal that is not the order type of any recursive well-ordering of $\mathbb{N}$ is much larger than $\varepsilon_0.$ It's called the Church-Kleene $\omega_1,$ and it's written as $\omega_1^{\,CK}.$