Good to see you found the Googology Wiki!
Beyond \varepsilon_0 we have:
\varepsilon_1 = \varepsilon_0 ^ {\varepsilon_0 ^ {\varepsilon_0 ^\ddots}}
\varepsilon_2 = \varepsilon_1 ^ {\varepsilon_1 ^ {\varepsilon_1 ^\ddots}}
\varepsilon_\omega = \sup \lbrace \varepsilon_0, \varepsilon_1, \varepsilon_2, \ldots \rbrace
and so on. Eventually we get to
\zeta_0 = \varepsilon_{\varepsilon_{\varepsilon_\cdots}}
We can set \varphi(0, \alpha) = \omega^\alpha, \varphi(1, \alpha) = \varepsilon_\alpha, \varphi(2, \alpha) = \zeta_\alpha. More generally,
\varphi(\alpha+1, \beta) = the \betath fixed point of f(\gamma) = \varphi(\alpha, \gamma).
When \alpha is a limit ordinal, \varphi(\alpha, \beta) is the \betath ordinal in the intersection in the ranges of f(\delta) = \varphi(\gamma, \delta) for all \gamma < \alpha.
So, for example,
\varphi(3, 0) = \varphi(2, \varphi(2, \varphi(2, \ldots)))
\varphi(\omega, 0) = \sup (\varphi(1, 0), \varphi (2, 0), \varphi (3, 0), \ldots)
and so on. This takes us up to
\Gamma_0 = \varphi(\varphi(\varphi(\ldots, 0),0),0)
or alternatively, \Gamma_0 is the smallest ordinal \alpha such that \alpha = \varphi(\alpha, 0).
We can continue the notation with \varphi(1, 0, 0) = \Gamma_0, and \varphi(1, 0, \alpha) is the \alphath fixed point of f(\beta) = \varphi(\beta, 0).
The ordinals \varphi (1, \alpha, \beta) are defined analogously to \varphi(\alpha, \beta), i.e. each function f(\beta) = \varphi (1, \alpha+1, \beta) enumerates the fixed points of g(\beta) = \varphi(1, \alpha, \beta), and at limit ordinals you enumerate the intersection of the ranges of previous ordinals. Then \varphi(2, 0, \alpha) is the \alphath fixed point of f(\beta) = \varphi(1, \beta, 0), and you can construct the hierarchy \varphi(2, \alpha, \beta) similarly. At limit ordinals you take the intersection of the ranges again, and so we define \varphi (\alpha, \beta, \gamma) for all ordinals \alpha, \beta, \gamma. Then we can define
\varphi (1, 0, 0, \alpha) to be the \alphath fixed point of f(\beta) = \varphi(\beta, 0, 0).
We can then define \varphi(\alpha, \beta, \gamma, \delta), \varphi(\alpha, \beta, \gamma, \delta, \epsilon), and so on. The general definition for the n-ary Veblen funciton is:
\varphi (\alpha) = \omega^{\alpha}
\varphi (\alpha_1, \alpha_2, \ldots, \alpha_n + 1, 0, \ldots, 0, \beta) is the \betath fixed point of the function f(\gamma) = \varphi(\alpha_1, \alpha_2, \ldots, \alpha_n, \gamma, 0, \ldots, 0)
When \alpha_n is a limit ordinal, \varphi (\alpha_1, \alpha_2, \ldots, \alpha_n, 0, \ldots, 0, \beta) is the \betath ordinal in the intersection of the ranges of f(\delta) = \varphi(\alpha_1, \alpha_2, \ldots, \alpha_{n-1}, \gamma, \delta, 0, \ldots, 0) for all \gamma < \alpha_n.
We define the Small Veblen Ordinal as
\sup (\varphi (1, 0), \varphi(1, 0, 0), \varphi(1, 0, 0), \ldots).
TREE(n) is larger than the fast-growing hierarchy at the level of the Small Veblen Ordinal. (how much further is not known.).
I'll go one step further: we can extend the n-ary Veblen function to transfinitely many places. Obviously, we can't write out transfinitely many variables, so we need to modify our notation: instead of
\varphi(\alpha, \beta, \gamma, \delta, \epsilon)
we write
\varphi(\alpha @ 4, \beta @ 3, \gamma @ 2, \delta @ 1, \epsilon @ 0).
So we append "@ n" to every variable, where n represents the index of the variable. This allows us to skip variables that are 0, and so we can notate things like \varphi (1 @ \omega, \alpha @ 0). \varphi(1 @ \omega, \alpha @ 0) is defined as the \alphath ordinal that is a fixed point of f(\beta) = \varphi (\beta @ n) for all n < \omega.
More generally, we define
\varphi(\alpha_1 @ \beta_1, \alpha_2 @ \beta_2, \ldots, \alpha_n + 1 @ \beta_n + 1, \gamma @ 0) is the \gammath fixed point of f(\delta) = \varphi(\alpha_1 @ \beta_1, \alpha_2 @ \beta_2, \ldots, \alpha_n @ \beta_n + 1, \delta @ \beta_n)
When \alpha_n is a limit ordinal, \varphi(\alpha_1 @ \beta_1, \alpha_2 @ \beta_2, \ldots, \alpha_n @ \beta_n + 1, \gamma @ 0) is the \gammath ordinal in the intersection of the ranges of f(\epsilon) = \varphi(\alpha_1 @ \beta_1, \alpha_2 @ \beta_2, \ldots, \delta @ \beta_n + 1, \epsilon @ \beta_n) for all \delta < \alpha_n
When \beta_n is a limit ordinal, \varphi(\alpha_1 @ \beta_1, \alpha_2 @ \beta_2, \ldots, \alpha_n + 1 @ \beta_n, \gamma @ 0) is the \gammath ordinal in the intersection of the ranges of f(\epsilon) = \varphi(\alpha_1 @ \beta_1, \alpha_2 @ \beta_2, \ldots, \alpha_n @ \beta_n, \epsilon @ \delta) for all \delta < \beta_n
When \alpha_n and \beta_n are limit ordinals, \varphi(\alpha_1 @ \beta_1, \alpha_2 @ \beta_2, \ldots, \alpha_n @ \beta_n, \gamma @ 0 is the \gammath ordinal in the intersection of the ranges of f(\zeta) = \varphi(\alpha_1 @ \beta_1, \alpha_2 @ \beta_2, \ldots, \delta @ \beta_n, \zeta @ \epsilon) for all \delta < \alpha_n, \epsilon < \beta_n
This defines \varphi(\alpha_1 @ \beta_1, \ldots, \alpha_n @ \beta_n) for all \alpha_i and \beta_i. This notation is known as Schutte's Klammersymbolen.
The smallest ordinal \alpha such that \alpha = \varphi (1 @ \alpha) is known as the Large Veblen Ordinal. I would think that TREE(n) would not reach the Large Veblen Ordinal in the fast-growing hierarchy, but I don't think this is known.
Phew! I hope this was at least somewhat comprehensible.