The map $f:\varepsilon_{1}\rightarrow\mathbb{N}$ which I am trying to define has to send $\varepsilon_{0}$ to some natural number. Since $\varepsilon_{0}=\omega\uparrow^{2}\omega$, a potential candidate is $f(\varepsilon_{0})=P_{f(f(\omega))}=P_{f(3)}=P_{8}=19$ (from the calculations for small ordinals below).
More generally, $f((\omega\uparrow^{n}\alpha)\cdot k) = P_{f^{n}(\alpha)}^{k}$, where $n,k$ are finite, $f^{n}$ is simply $f$ composed with itself $n$ times, and $P_{i}$ denotes the $i$-th prime.
Some small ordinals will be mapped as follows: $f(0)=1$, $f(1)=f(\omega\uparrow^{1}0)=P_{f(0)}=P_{1}=2$, and more generally, $f(m)=f((\omega\uparrow^{1}0)\cdot m)=2^{m}$, $m\ge 1$. Finally, $f(\omega) = f(\omega\uparrow^{1}1)=P_{f(1)}=P_{2}=3$.
Am I going in the right direction?
UPDATE: Something bothers me: $f$ may not well defined because the $\uparrow$ does not produce a unique notation. For example, $\omega=\omega\uparrow^{1} 1 = \omega\uparrow^{2} 1$. Is there a way to avoid this pitfall? Perhaps, when an ordinal $\alpha$ can be denoted by several expressions, we take the smallest of the resulting numbers after applying $f$ to each of the expressions denoting $\alpha$? For example, $f(\omega)=f(\omega\uparrow^{1}1)$. Is this a correct argument to get out of the problem?
UPDATE: I am not sure I am using $\uparrow$ correctly: for example $\omega\uparrow^{1}\omega = \omega^{\omega}$. Then $\omega\uparrow^{2}\omega = \underbrace{\omega\uparrow^{1}(\omega\uparrow^{1}(\omega\uparrow^{1}\omega\cdots}_{\omega}=\underbrace{\omega^{\omega^{\omega^{...}}}}_{\omega}=\varepsilon_{0}$. But, what is $\omega\uparrow^{3}\omega$? $\omega\uparrow^{3}\omega=\underbrace{\omega\uparrow^{2}(\omega\uparrow^{2}(\omega\uparrow^{2}\omega\cdots}_{\omega}=\underbrace{\omega\uparrow^{1}(\omega\uparrow^{1}(\omega\uparrow^{1}\omega\cdots}_{\omega\uparrow^{2}\omega}\stackrel{?}{=}\underbrace{\omega^{\omega^{\omega^{...}}}}_{\varepsilon_{0}}=\;?$
This seems bigger than $\varepsilon_{1}$...In fact I am beginning to think it depends on how $\uparrow$ is defined, and there may be several ways to do so.
I have encountered two problems when using $\uparrow$ notation.
The first is how to define the $\uparrow$ operator properly. As far as I can tell, there is not a widely accepted definition of it for transfinite ordinals. For example, $\omega\uparrow^{2}(\omega+1)$ may be equal to $\varepsilon_{0}^{\omega}$ or to $\varepsilon_{0}^{\varepsilon_{0}}$, depending on the definition (see This and this MSE questions and their answers for at least three different definitions).
The second problem is that there does not seem to be a unique way to write and ordinal using $\uparrow$. For example, $\omega^{\omega}=\omega\uparrow^{2}2 = \omega\uparrow^{1}\omega$. Therefore, this notation does not yield an analog of the Cantor Normal Form.
In summary, the $\uparrow$ operator is not the best tool to appropriately express ordinals- even those below $\varepsilon_{1}$. My initial hope was that we could express $\alpha<\varepsilon_{1}$ uniquely as a sum $(\omega\uparrow^{i_{1}}\alpha_{1})+\cdots+(\omega\uparrow^{i_{n}}\alpha_{n})$, where $i_{k}\in\{1,2\}$ for $1\le i\le n$, and $\alpha_{i}<\varepsilon_{1}$. No chance of it, it seems.
The silver lining is that the harder I tried to work out a proper definition of $\uparrow$, the more it seemed like I should be using Veblen functions (here).In fact, it is enough to use a slight modification of Veblen functions (See Lemma 8.7 here) in order to use the following theorem:
(Theorem 8.12 here): every ordinal $0<\alpha<\Gamma_{0}$ (Feferman–Schütte ordinal) can be expressed uniquely as $\alpha = \psi_{\alpha_{1}}(\beta_{1})+\cdots+\psi_{\alpha_{n}}(\beta_{n})$, $n\ge 1$ with $\alpha_{i},\beta_{i}<\psi_{\alpha_{i}}(\beta_{i})\leq\alpha$, $1\leq i\leq n$, and such that $\psi_{\alpha_{1}}(\beta_{1})\ge\cdots\ge\psi_{\alpha_{n}}(\beta_{n})$.
Therefore, we can do a bit better than simply answering the original question, and define an explicit bijection $f:\Gamma_{0}\rightarrow\mathbb{N}$. Consider the set $\mathbb{N}^{*}=\mathbb{N}\cup\{0\}$ and the Cantor pairing function $q:\mathbb{N}^{*}\times\mathbb{N}^{*}\rightarrow\mathbb{N}^{*}$ defined by $q(x,y) = \frac{1}{2}(x+y)(x+y+1)+x$. Let $P_{0}=2$, $P_{1} = 3$, etc. denote the prime numbers. Then
Note that $f(\alpha+\beta) = f(\alpha)\cdot f(\beta)$. The ordering of the natural number thus produced looks like this:
Natural number [position in the ordering]
$$\begin{aligned} 1\ [0] &≤ 2\ [1] ≤ 4\ [2] ≤ 8\ [3] ≤ 16\ [4] ≤ \ldots \\ &\quad≤ 3\ [ω] ≤ 6\ [ω+1] ≤ 12\ [\omega + 2] ≤ \ldots \\ &\quad≤ 9\ [ω\cdot 2] ≤ 18\ [\omega\cdot 2+1] ≤ 36\ [\omega\cdot 2+2] ≤ \ldots \\ &\quad≤ 17\ [\omega^{2}] ≤ \ldots \\ &\quad≤ \ldots \\ &≤ 7\ [\omega^{\omega}] ≤ \ldots ≤ 5\ [\varepsilon_{0}] ≤ \ldots\leq 25\ [\varepsilon_{0}\cdot 2]\leq\ldots ≤ 31\ [\varepsilon_{0}\cdot\omega] ≤ \ldots\leq 1993\ [\varepsilon_{0}^{2}]\leq\ldots\\ & ≤ 3313\ [\varepsilon_{0}^{\omega}] ≤ \ldots\ldots \leq 11\ [\varepsilon_{1}]\leq\ldots\ldots\leq 37\ [\varepsilon_{2}]\leq\ldots\ldots\leq 47\ [\varepsilon_{\varepsilon_{0}}]\leq\ldots\\ &\leq 29\ [\zeta_{0}]\leq\ldots\ldots\leq 13\ [\psi_{\omega}(0)]\leq\ldots\ldots\leq 71\ [\psi_{\varepsilon_{0}}(1)]\leq\ldots \end{aligned}$$
It is not difficult to write $f^{-1}$ explicitly. Also, we could use a different pairing function in the definition of $f$. For example, $q'(x,y) = 2^{x}(2y+1) - 1$. In this case $\omega$ will be mapped to $5$ and $\varepsilon_{0}$ to $3$, etc. Therefore, in this ordering $5$ will be less than 3.
Finally, I will briefly mention that if instead of defining $q$ above as the Cantor pairing function, we define it as $q:\{(1,0)\}\cup(\{0\}\times\mathbb{N})\rightarrow \mathbb{N}$ by $$q(x,y) = \begin{cases} x &\text{if $x=1$ or $x=y=0$}\\ y+1 &\text{otherwise} \end{cases}$$ then $f|_{\varepsilon_{1}}:\varepsilon_{1}\rightarrow\mathbb{N}$ is a bijection which produces the ordering
$$\begin{aligned} 1\ [0] &≤ 2\ [1] ≤ 4\ [2] ≤ \ldots \\ &\quad≤ 5\ [ω] ≤ 10\ [ω+1] ≤ 20\ [\omega + 2] ≤ \ldots \\ &\quad≤ 25\ [ω\cdot 2] ≤ 50\ [\omega\cdot 2+1] ≤ 100\ [\omega\cdot 2+2] ≤ \ldots \\ &\quad≤ 11\ [\omega^{2}] ≤ \ldots \\ &\quad≤ \ldots \\ &≤ 13\ [\omega^{\omega}] ≤ \ldots ≤ 3\ [\varepsilon_{0}] ≤ \ldots\leq 9\ [\varepsilon_{0}\cdot 2]\leq\ldots ≤ 7\ [\varepsilon_{0}\cdot\omega] ≤ \ldots\leq 29\ [\varepsilon_{0}^{2}]\leq\ldots\\ & ≤ 19\ [\varepsilon_{0}^{\omega}] ≤ \ldots\leq 109\ [\varepsilon_{0}^{\varepsilon_{0}}]\leq\ldots \end{aligned}$$