Python 3, f_{PTO(Δ¹2-CA+BI+Π¹2-CA^-)}(2), 90 bytes

a=[0,1,2,3]*2
while a:v=a.pop();a+=[i+(len(a)-v)*(i>=v-v%4)for i in a[v:]];print(v,end="")

Try it online!
Ungolfed:

n = 4
a = [*range(n)]*2
while a:
    v = a.pop()
    l = len(a) - v
    for i in a[v:]:
        if i < v - v%n:
            a.append(i)
        else:
            a.append(i + l)
    print(v,end="")

Human-readable output:

n = 4
a = [*range(n)]*2
while a:
    v = a.pop()
    l = len(a) - v
    for i in a[v:]:
        if i < v - v%n:
            a.append(i)
        else:
            a.append(i + l)
    b = []
    for i in a:
        if i>=len(b):
            b.append(0)
        else:
            b.append(b[i] + 1)
    c = [tuple(b[n*i:n*(i+1)]) for i in range((len(b) - 1)//n + 1)]
    print(*c, sep="")

This simulates BMS, but using a more nontrivial input would require two additional bytes (replace [0,1,2,3] with [*range(9)] in the first line and %4 with %9 in the second line). In particular, it stores a concatenation of all the columns of the matrix, with each element being stored as the index of its parent (and zeroes being stored as their own index).

The bound for the runtime is not formally proven, but the same is true about the PTO(Z₂) bound in the other answers involving BMS, so I'm guessing that's ok. The matrix (0,0,0)(1,1,1)(2,2,1)(3,2,1)(2,1,0)(3,2,0) should correspond to ψ(ε_(α+1)) where α is the smallest (next recursively inaccessible)-stable ordinal (i.e. smallest β that is γ-stable with γ being a recursively inaccessible ordinal above β). According to Michael Rathjen's "An ordinal analysis of parameter free Π¹₂-comprehension", this should be the proof-theoretic ordinal of KPi+"there exists a stable ordinal", or in terms of arithmetic, Δ¹₂-CA+BI+Π¹₂-CA^-. The code reaches (0,0,0)(1,1,1)(2,2,1)(3,2,1)(3,2,0), which is quite a bit bigger (in fact it should be ψ(α) where α is stable up to the next fixed point of enumerating rec. inaccessibles), and there should be plenty of points between PTO(Δ¹₂-CA+BI+Π¹₂-CA^-) and this where SGH catches up to FGH, so it's safe to use the fast-growing hierarchy. Of course, all this depends on the exact details of what OCF is being used and which system of fundamental sequences we are choosing for the FGH, but with any reasonable choice, the bound will hold.

The number that is printed is simply the concatenation of all popped elements, so it's lower bounded by roughly 10^runtime, and since the first popped element is nonzero, it is indeed a number. Almost 20% of the code was the print statement, so if you just want long runtime, you can do this in only 74 bytes by simply removing the print. If you want to beat all the other entries below Loader's number, you can do that in 98 bytes:

a=[*range(n:=9**9)]*2
while a:v=a.pop();a+=[i+(len(a)-v)*(i>=v-v%n)for i in a[v:]];print(v,end="")

JavaScript (Node.js), H_PTO(Z2)(1000), 480 bytes

lim(BMS)[1000]>lim(BMS)[32]

let b=(e)=>{return e[e.length-1].every(it=>!it)};let a=(s,i)=>{let lnz=s.at(-1).findLastIndex(it=>!!it);let r = s.length-s.at(-1)[lnz]-1;let j=s.at(-1);s[s.length-1]=s.at(-1).map((it,id)=>id<lnz?it:!s[r][id]?0:s[r][id]+s.at(-1)[lnz]);let b=s.slice(r+1);for(var h=0;h<i;h++){b=b.map((it,id)=>it.map(it2=>it2<=id+1?it2:it2+j[lnz]));s.push(...b);}return s};s=[Array(1000).fill(0),Array(1000).fill(1)];i=1000;while(s.length){b(s)?(()=>{s=s.slice(0,-1);i++})():s=a(s,i)};console.log(i)

Beautified:

let b = (e) => {
    return e[e.length - 1].every(it => !it)
};
let a = (s, i) => {
    let lnz = s.at(-1).findLastIndex(it => !!it);
    let r = s.length - s.at(-1)[lnz] - 1;
    let j = s.at(-1);
    s[s.length - 1] = s.at(-1).map((it, id) => id < lnz ? it : !s[r][id] ? 0 : s[r][id] + s.at(-1)[lnz]);
    let b = s.slice(r + 1);
    for (var h = 0; h < i; h++) {
        b = b.map((it, id) => it.map(it2 => it2 <= id + 1 ? it2 : it2 + j[lnz]));
        s.push(...b);
    }
    return s
};
s = [Array(1000).fill(0), Array(1000).fill(1)];
i = 1000;
while (s.length) {
    b(s) ? (() => {
        s = s.slice(0, -1);
        i++
    })() : s = a(s, i)
};
console.log(i)

Clear version of pointer bms code:

function point(s,i){
    var lnz = s[s.length-1].findLastIndex(item=>!!item);
    var r = s.length-s[s.length-1][lnz]-1;
    var j=s[s.length-1]
    s[s.length-1]=s[s.length-1].map((item,id)=>id<lnz?item:!s[r][id]?0:s[r][id]+s[s.length-1][lnz]);
    var b=s.slice(r+1,s.length);
    for(var h=0;h<i;h++){
        b=b.map((item,id)=>item.map(item2=>item2<=id+1?item2:item2+j[lnz]))
        s.push(...b);
    }
    return s;
}

This is the largest entry below Loader's number we currently have so far. Try it online!

Japt, 39 Bytes, f_PTO(Z2) (32) or BMS[32]

This is a port of this answer.

T°?UÊ?ßUcUsX=Uo)Ë<XfH ?D:D+UÊ-X:T:ßLÇ%H

I made a few golfing edits compared to my previous ;T°?UÊ?ßUcUsX=Uo)m_<X-X%BÊ?Z:Z+UÊ-X:T:ßBÊo cBÊo. Since I am able to use constants, I used L=100 instead of importing the alphabet and getting its length. I also used H=32, which allowed me to shave a few bytes in creating the BMS array. I also used a few other Japt substitutions.

Javascript, 511 bytes

I am new to code golf (this is my second time), but certainly not to implementing or inventing large numbers. However this programs takes quite a lot of code to beat TREE(3) :( The program results in \$f_{TFBO+1}(3) = f_{TFBO}(f_{TFBO}(f_{TFBO}(3)))\$ in the fast-growing hierarchy, where \${TFBO}\$ is the Takeuti-Feferman-Buchholz Ordinal \$\psi(\Omega_{\omega+1})\$. It uses Buchholz hydras and implements the BH function. Here is my compressed code, which was done with Naruyoko's Javascript compressor:

_="69{%t.length-1D*6m#{%B]}*6s#{t.splice(m#,14*@n=1;*6e#{7==0/s#;if(9>1/$i=0;i<n;i:t'B-1]4}}7>0&&typeof t==\"number\"/$i=98t[i]<m#8i--/}@s;$j5i8j <= 98j:s't[j]4@?5s;?[j]=?[i]-1;?[l(?)]=0;s#;$k508k <= l(?)8k:t'?[k]4}7==A/B]=n+1Dn += 1;%tD*6b(k/i5[\"+\",0];$l508l < k-28l:i'A)}while(l(i)>=0/i=e(i4%nD*Console.log(b(b(b(3))));#(t)$for(let %return '.push(*\n    /){4)D5 = 6function 7if(m#8; 9l#:++/?sa@var A\"w\"Bt[9D;}";for($ of"DBA@?:987654/*'%$#")with(_.split($))_=join(pop());eval(_)

This is basically incomprehensible, with only fragments of the original code revealing themselves, so here is the original code (still somewhat golfed):

function l(t){return t.length-1;}
function m(t){return t[l(t)]}
function s(t){t.splice(m(t),1);}
var n=1;
function e(t){if(m(t)==0){s(t);if(l(t)>1){for(let i=0;i<n;i++){t.push(t[l(t)-1]);}}}if(m(t)>0&&typeof t=="number"){for(let i=l(t); t[i]<m(t); i--){}var s;for(let j = i; j <= l(t); j++){s.push(t[j]);}var sa = s;sa[j]=sa[i]-1;sa[l(sa)]=0;s(t);for(let k = 0; k <= l(sa); k++){t.push(sa[k]);}}if(m(t)=="w"){t[l(t)]=n+1;}n += 1;return t;}
function b(k){i = ["+",0];for(let l = 0; l < k-2; l++){i.push("w")}while(l(i)>=0){i=e(i);}return n;}
Console.log(b(b(b(3))));

Explanation

This program, as I mentioned earlier implements Buchholz hydras, which are stored as arrays. This is what the code does:

It defines shorthands l(t), the index of the last element of the array, m(t), the last element of the array, and s(t), which chops off the last element of the array (the rightmost head). I also define a counter n which keeps track of how many times the hydra has been expanded.

The biggest block of code is e(t), which expands the hydra. Here is how it works: it checks if the rightmost head (which it will chop off) has label zero. If it does, then it calls s(t), splicing the last item in the array, and checks if l(t)>1, i.e. the rightmost head's parent is not the root. If it is, then it uses a for loop to push n copies of the rightmost head's parent (which is not the root).

Next, if it is not zero, it checks if it is a positive integer. If it is, it uses a for loop to find i such that t[i] has a lower label than the rightmost head. Then it creates a tree s such that everything after (and including) t[i] is in s, i.e. the subtree rooted at t[i]. It also creates sa where the rightmost head is relabeled with sa[i]-1 and the rightmost head of sa is relabeled with 0. It then deletes t's rightmost head and uses a for loop to insert sa into t.

If it is not zero, then it can only be the least limit ordinal omega (represented by the string "w"). If it is, it simply relabels it with n+1. After everything, it increments n by 1 and returns the new version of t.

Next, we have the function b(n), which represents the function BH(n). It sets i to the tree [+,0] and then pushes n-2 copies of omega "w". It then repeatedly applies e to i (expands i) and returns the final output. My large number is then b(b(b(n)))

Size

Why does this produce a large enough number? Well, the original inventor, Wilfried Buchholz, proved that this function eventually dominates (overtakes) any function which is provably total in \$\Pi^1_1\$-comprehension with bar induction (\$\Pi^1_1-CA+BI\$), so the \$BH(n)\$ function is commonly approximated to \$\psi_0(\varepsilon_{\Omega_\omega +1})\$ using the fast-growing hierarchy, so my number is approximately \$f_{ψ(\varepsilon_{\Omega_\omega+1})+1}(3)\$.

The goal is \$\psi_0(\Omega^{\Omega^\omega \omega})\$, which is much less than \$\psi_0(\varepsilon_{\Omega+1})\$, let alone \$\psi_0(\varepsilon_{\Omega_\omega+1})+1\$. In fact, not only does this number defeat TREE(3), it also defeats SCG(13), which is around \$\psi_0(\Omega_\omega)\$, although it still falls massively short of Loader's number.

In terms of size, my entry would come second place, while in terms of length, it would come second-to-last. The latter is what matters more, so this is probably a really bad entry.

Python, 359 (184 without whitespace) 224 212 bytes, ~$$H_{ψ(ψ_Ι(0))[7]}(8)$$

*S,n=0,2,9
while S:
 i=r=d=0
 for a in S[::-1]:
  i-=1
  if d<S[-1]-a:
   r=r or i;d+=S[-1]-a-d>1
   if[k+d for k in S[i:]]<S[r:]:break
   p=i+1
 S.pop();n*=n
 if r:S=S[:p]+[d*i+k for i in range(n)for k in S[p:]]

I'm new to codegolf so I'm not sure if all the whitespace here counts. If not, this entry isn't the shortest entry, but vastly dominates patcail's entry in size. The code here implements Sudden Sequence System (definition here), which is known to have a perfect correspondence to ordinals below $$\psi(\psi_I(0))$$ represented in Extended Weak Buchholz function. The original idea uses sequences lexicographically less than (0,2) to represent ordinals, but to save characters the sequence itself is reversed. There are actually a few open problems here that we don't know the answer to, but strongly suspect they are true:

1. Does the limit of Extended Weak Buchholz function equal to limit of Extended Buchholz function?
2. Is Sudden Sequence System terminating?
3. Does the correspondence between the two systems hold?

Edit 1: It's broken due to a silly mistake. Now S is correctly initialized. Thanks to Naruyoko for a few great optimizations.

Edit 2: Optimized a bit

JavaScript, 190 bytes, H_ψ(εΩ+1)(9)_{Based off of this analysis}

A=[0,1,2];B=[0,1,2];for(F=C=9;F;F--){for(C++,G=0;G<=F;G++)(A[F]||A[F-G]<A[F]-H)&&(B[F]?(I=G,G=F):(H=A[F]-A[F-G],B[F-G]<B[F]&&(I=G,G=F)));for(J=0;J<C*I;J++)A[F]=A[F-I]+H,B[F]=B[F-I],F++;H=0}C

This program is a modified version of this 225B Pair-sequence number translation in JavaScript. For Pair-sequence number and their original code, see here.

The modifications done:

• It is in JavaScript instead of BASIC.
• No iteration(f_ψ(Ωω+1)->f_ψ(Ωω))
• The sequence is (0,0)(1,1)(2,2), which corresponds to ordinal ψ(ε_Ω+1)._{This is in Hardy-hierarchy ordinal}

Javascript, 111 bytes, ~ \$f_{\psi(\Omega_\omega)}(6)\$

\$f\$ is the Fast-growing Hierarchy. \$ψ\$ is Buchholz's Psi. This entry, despite being 111 bytes, dominates all of the previous entries in both size and the amount of bytes (except for Loader's number).

Here is the code:

s=JSON.stringify;P=([y,z])=>y?JSON.parse((k=s([P(y),z])).replaceAll(s(z),k)):z;for(a=b=0;a=b++<9?[a,0]:P(a););b

Here is the same code expanded out:

function P([y,z]) {
  if (y==0) {
    return z
  } else {
    k = JSON.stringify([P(y),z])
    return JSON.parse(k.replaceAll(JSON.stringify(z),k))
  }
}

for(a=b=0;a=b++<9?[a,0]:P(a););
b;

I'm going to explain both the P function and the for loop.

The Predecessor Function

The inputs of the predecessor function are binary trees with zeroes as leaf nodes. Here are some examples of binary trees:

• 0
• [0,0]
• [[0,[0,0]],[0,0]]
• [[[[[0,0],0],[0,0]],[0,0]],[[0,[0,0]],[0,[0,0]]]]
• [[[[[[[[[[0,0],0],0],0],0],0],0],0],0],0]

The Predecessor function is defined like this:

• P([0,z])=z
• P([x,y])=[P(x),y] but with all instances of y replaced with [P(x),y]
• P(0) is left undefined

Right away, we can see 0 represents the number \$0\$, and [0,z] represents the structure \$z+1\$. Natural numbers can be represented as [0,[0,[0,...[0,0]...]]] with \$n+1\$ zeroes. For example, \$1 =\$ [0,0], \$2 =\$ [0,[0,0]], \$3 =\$ [0,[0,[0,0]]], and so on.

Now consider the string [1,n] where \$n>1\$.

P([1,n])=[0,n] but replace all instances of \$n\$ with [0,n] \$\to\$ [0,[0,n]]

Therefore, [1,n] corresponds to \$n+3\$, as P(P(P([1,n]))) = n

By this logic, [2,n] corresponds to \$n+7\$, [3,n] corresponds to \$n+15\$, and [n,n] would approximately correspond to \$2^n\$. Maybe [[0,n],n] corresponds to \$2^{n+1}\$?

Not so fast!

Consider the string [[0,n],n]. One would expect this to correspond to \$2^{n+1}\$, but it is much stronger. P([[0,n],n]) \$\to\$ [P([0,n]),n] = [n,n], but then the second step would be to replace all instances of n with the entire tree, or [n,n]. This makes P([[0,n],n])=[[n,n],[n,n]] rather than [n,[n,n]].

One would ask whether this would cause an infinite loop. Let's try P([[n,n],[n,n]]). If we let J = P([n,n]), we will get:

P([[n,n],[n,n]])=[J,[n,n]] but with all instances of [n,n] replaced with [J,[n,n]]

However, there are no instances of [n,n] within J, because J is strictly less than [n,n]. Therefore, P([[n,n],[n,n]])=[J,[J,[n,n]]]. This works for all J less than [n,n].

So this means [[0,n],n] corresponds to \$2^{2^n}\$. [[0,[0,n]],n] corresponds to \$2^{2^{2^{2^n}}}\$. And finally, [[n,n],n] corresponds to \$n \uparrow\uparrow n\$. Now it is time to bring in the Middle Growing Hierarchy.

Middle Growing Hierarchy

The Middle Growing Hierarchy is defined here: https://googology.wikia.org/wiki/Middle-growing_hierarchy

One can make an approximate distinction with the Middle Growing Hierarchy.

• [0,n] corresponds to \$m(0,n) \sim n+1\$
• [1,n] corresponds to \$m(2,n) \sim n+3\$
• [2,n] corresponds to \$m(3,n) \sim n+7\$
• [n,n] corresponds to \$m(n,n) \sim n+2^n\$ and \$m(ω,n)\$
• [n,[n,n]] corresponds to \$m(n+1,n) \sim 2^{n+1}\$ (not \$m(ω+1,n)\$)
• [[0,n],n] corresponds to \$m(ω+1,n) \sim 2^{2^n}\$
• [[0,[0,n]],n] corresponds to \$m(ω+2,n) \sim 2^{2^{2^{2^n}}}\$
• [[1,n],n] corresponds to \$m(ω+3,n) \sim 2^{2^{2^{2^{2^{2^{2^{2^n}}}}}}}\$
• [[n,n],n] corresponds to \$m(ω2,n) \sim n \uparrow\uparrow n\$

One can see a correspondence with the left-hand side of the binary tree and the inner subscript of the Middle Growing Hierarchy. Let's continue the correspondence. I will omit the right hand side of the binary tree and the base of the Middle Growing Hierarchy.

• [n,n] corresponds to \$ω_2\$
• [0,[n,n]] corresponds to \$ω_2+1\$
• [n,[n,n]] corresponds to \$ω_3\$
• [[0,n],[n,n]] corresponds to \$ω_4\$
• [[1,n],[n,n]] corresponds to \$ω_6\$
• [[0,n],n] corresponds to \$ω^2\$

So it seems like there is a jump from [[n,n],n] to [[[0,n],n],n], similar to the jump from [n,n] to [[0,n],n]. But even this doesn't capture the power of this notation.

More Ordinal Comparison

• [n,[[0,n],n]] corresponds to \$ω^2+ω\$
• [[n,n],[[0,n],n]] corresponds to \$2\timesω^2\$
• [[n,[n,n]],[[0,n],n]] corresponds to \$ω^3\$
• [[[0,n],[n,n]],[[0,n],n]] corresponds to \$ω^4\$
• [[[0,n],n],[[0,n],n]] corresponds to \$ω^ω\$

We're not even at [[0,[0,n]],n] yet, what is going on?

• [[[0,n],n],[[[0,n],n],[[0,n],n]]] corresponds to \$2\timesω^ω\$
• [[n,[[0,n],n]],[[[0,n],n],[[0,n],n]]] corresponds to \$ω^{ω+1}\$
• [[[n,n],[[0,n],n]],[[[0,n],n],[[0,n],n]]] corresponds to \$ω^{ω_2}\$
• [[n,[[n,n],[[0,n],n]]],[[[0,n],n],[[0,n],n]]] corresponds to \$ω^{ω_2+1}\$
• [[[0,[n,n]],[[0,n],n]],[[[0,n],n],[[0,n],n]]] corresponds to \$ω^{ω_3}\$
• [[[n,[n,n]],[[0,n],n]],[[[0,n],n],[[0,n],n]]] corresponds to \$ω^{ω^2}\$
• [[[[0,n],[n,n]],[[0,n],n]],[[[0,n],n],[[0,n],n]]] corresponds to \$ω^{ω^3}\$
• [[0,[0,n]],n] corresponds to \$ω^{ω^ω}\$

Where does the strength comes from? The strength lies in the fact that the binary tree notation corresponds to performing the Middle Growing Hierarchy on ordinals! Here are some examples:

• [n,n] corresponds to \$m(ω,n)=n+2^n\$, so [n,n] as an ordinal corresponds to \$ω+2^ω=ω_2\$
• [[0,n],n] corresponds to \$m(ω+1,n)=2^{n+2^n}\$, so [[0,n],n] as an ordinal corresponds to \$2^{ω+2^ω}=2^{ω_2}=ω^2\$
• [[[0,n],n],[[0,n],n]] corresponds to \$m(ω,m(ω+1,n))=2^{2^{n+2^n}}\$, so [[[0,n],n],[[0,n],n]] as an ordinal corresponds to \$2^{2^{ω+2^ω}}=ω^ω\$
• [[0,[0,n]],n] corresponds to \$m(ω+2,n)=2^{2^{2^{n+2^n}}}\$, so [[0,[0,n]],n] as an ordinal corresponds to \$2^{2^{2^{ω+2^ω}}}=ω^{ω^ω}\$

As it turns out, this pattern continues. I'm not going to go through the full analysis, but here are some more ordinal values. Remember that these are ordinals, not functions!

• [[1,n],n] corresponds to \$ω^{ω^{ω^{ω^{ω^{ω^ω}}}}}\$
• [[n,n],n] corresponds to \$ε_0\$
• [[n,[n,n]],n] corresponds to \$ζ_0\$
• [[[0,n],n],n] corresponds to \$φ(ω,0)\$
• [[[n,n],n],n] corresponds to the BHO

Speed of Notation

Essentially, if a structure \$K\$ corresponds to \$g_a (n)\$ in the slow-growing hierarchy, then the structure [K,n] corresponds to \$m_a (n)\$ in the middle-growing hierarchy. This makes the limit [[[...,n],n],n], which corresponds to the first SGH-MGH catching point, of \$ψ(Ω_ω)\$. For comparison, \$\text{TREE}(n)\$ only corresponds to the ordinal \$ψ(Ω^{Ω^ω})\$, much much smaller. The premise of this notation is essentially nested Goodstein sequences, except it works!

The Middle Growing Hierarchy corresponds closely to the Fast Growing Hierarchy, this is why I put in \$f_{ψ(Ω_ω)}\$ as it is a catching point.

Actual Value of the program

Now that we have gone through the Predecessor function, and how it corresponds to numbers, functions, and ordinals, it is time to return to the value of this program.

To extract a value from a binary tree, such as [[[0,0],0],0], one would have to repeatedly apply the predecessor function until the value crashes down to 0. As we seen before, one would have to apply the predecessor function a massive amount of time, on the order of \$m(ψ(Ω_ω),x)\$

Just to let you know, [[[...,0],0],0] is not degenerate, unlike stuff like \$2 \uparrow\uparrow...\uparrow\uparrow2 = 4\$ in arrow notation. [[[...,0],0],0] will produce a massive number.

Here is the code again:

for(a=b=0;a=b++<9?[a,0]:P(a););b;

First, it sets a and b equal to 0. Then, it starts incrementing b. If b is less than 9, then it sets a to [a,0]. This means at b=9, a would had been already [[[[[[[[[0,0],0],0],0],0],0],0],0],0], which corresponds to a massive number. Then, the predecessor function gets repeatedly applied to a, increasing b by 1 for each application. Eventually, a is going to crash down to 0, but b will be some value far, far greater than \$\text{TREE}(3)\$, or \$\text{TREE}(\text{TREE}(...\text{TREE}(3)...))\$ with \$\text{TREE}(3)\$ nests. Finally, the program returns b.

So what?

One of the best thing about this program is how the notation enumerates the catching points between the SGH and the MGH. This program only reaches the very first catching point, but by a few simple extension, this program is able to formalize a meameamealokkapoowa oompa, surpass Strong Array Notation, and beat every single Ordinal Collapsing Function ever devised. \$ψ(Ω_ω)\$ is still a pathetically small value...

Haskell, 252 Bytes, TREE(3)+1

data T=T[T]Int
l(T n _)=1+sum(l<$>n)
a@(T n c)#T m d=any(a#)m||c==d&&n!m
l@(x:t)!(y:u)=l!u||x#y&&t!u
x!_=null x
a n=do x<-[1..n];T<$>mapM(\_->a$n-1)[2..x]<*>[1..3]
s 0=[[]]
s n=[t:p|p<-s$n-1,t<-a n,(l t<=n)>any(#t)p]
main=print$[x|x<-[0..],null$s x]!!0

Thanks for help from H.PWiz, Laikoni and Ørjan Johansen for help golfing the code!

As suggested by HyperNeutrino, my program outputs TREE(3)+1, exactly (TREE is computable as it turns out).

T n c is a tree with label c and nodes n. c should be 1, 2, or 3.

l t is the number of nodes in a tree t.

t1 # t2 is true if t1 homeomorphically embeds into t2 (based on Definition 4.4 here), and false otherwise.

a n outputs a big list of trees. The exact list isn't important. The important property is that a n contains every tree up to n nodes, with nodes being labelled with 1, 2, or 3, and maybe some more trees as well (but those other trees will also be labelled with 1, 2, or 3). It is also guaranteed to output a finite list.

s n lists all sequences length n of trees, such that the reverse (since we build it backwards) of that sequence is valid. A sequence is valid if the nth element (where we start counting at 1) has at most n nodes, and no tree homeomorphically embeds into a later one.

main prints out the smallest n such that there is no valid sequences of length n.

Since TREE(3) is defined as the length of the longest valid sequence, TREE(3)+1 is the smallest n such that there are no valid sequences of length n, which is what my program outputs.

New Ruby, 135 bytes, >> H_{ψ(φ3(Ω+1))}(9)

where H is the Hardy hierarchy, ψ is an extended version of Madore's OCF (will explain below) and φ is the Veblen function.

Try it online!

f=->a,n,b=a{c,d,e=a;a==c ?a-1:e ?a==a-[0]?[[c,d,f[e,n,b]],d-1,c]:c:[n<1||c==0?n:[f[c||b,n-1]],n,n]};h=[],k=9,k;h=f[h,p(k*=k)]while h!=0

Ungolfed: (using functions, not lambdas)

def f(a,n,b)
  c,d,e = a
  if a == c
    return a-1
  elsif e
    if a == a-[0]
      return [[c,d,f(e,n,b)],d-1,c]
    else
      return c
    end
  else
    x = c || b
    if n < 1 || c == 0
      return [n,n,n]
    else
      return [f(x,n-1,x),n,n]
    end
  end
end

k = 9
h = [[],k,k]
while (h != 0) do
  k *= k
  p k
  h = f(h,k,h)
end

Madore's extended OCF:

And (crudely) Veblen's phi function:

Explanation without ordinals:

f(a,n,b) reduces an array recursively. (if no third argument given, it takes the first argument twice.)
f(k,n,b) = k-1, k is a positive int.
f([c,d,0],n,b) = f([c,0,e],n,b) = c
f([c,d,e],n,b) = [[c,d,f(e,n,b)],d-1,c], d ≠ -1 and c ≠ 0

f([a],0,b) = [0,0,0]
f([0],n,b) = [n,n,n]
f([],n,b) = f([b],n,b)
f([a],n,b) = [f[a,n-1,a],n,n]

My program initiates k = 9, h = [[],9,9]. It then applies k = k*k and h = f(h,k) until h == 0 and outputs k.

Explanation with ordinals:

Ordinals follow the following representation: n, [], [a], [a,b,c], where n,d is a natural number and a,c are all ordinals.
x = Ord(y) if y is the syntactic version of x.
a[n,b] = Ord(f(a,n))
ω = Ord([0]) = Ord(f([a],-1,b))
n = Ord(n)
Ω = Ord([])
ψ'(a) = Ord([a])
ψ'(a)[n] = Ord(f([a],n))
φ(b,c) ≈ Ord([[0],b,c])
a(↓b)c = Ord([a,b,c]) (down-arrows/backwards associative hyper operators I designed just for ordinals)

We follow the following FS for our ordinals:
k[n,b] = k-1, k < ω
ω[n,b] = n(↓n)n
(a(↓b)0)[n,b] = (a(↓0)c)[n,b] = a
(a(↓b)c)[n,b] = (a(↓b)(c[n,b]))(↓b[n,b])a, b ≥ 0 and c > 0.
ψ'(a)[0,b] = 0(↓0)0
ψ'(a)[n,b] = (ψ'(a[n-1,a]))(↓n)ω, a > 0 and n ≥ 0. (also note that we've changed from [n,b] to [n,a].)
Ω[n,b] = ψ'(b)[n,b]

ψ'(ω∙α) ≈ ψ(α), the ordinal collapsing function described in the image above.

My program more or less initiates k = 9 and h = Ω(↑9)9, then applies k ← k² and h ← h[k,h] until h = 1 and returns k.

And so if I did this right, [[],9,9] is way bigger than the Bachmann-Howard ordinal ψ(Ω^ΩΩ...), which is way bigger than ϑ(Ω^ωω)+1.

ψ(Ω(↓9)9) > ψ(Ω(↓4)3) > ψ(Ω^ΩΩ)+1 > ψ(Ω^Ωωω)+1 > ϑ(Ω^ωω)+1

And if my analysis is correct, then we should have ψ'(Ω^Ω∙x) ~= ψ*(Ω^Ω∙x), where ψ* is the normal Madore's psi function. If this holds, then my ordinal is approximately ψ*(φ₃(Ω+ω)).

Old Ruby, 309 bytes, H_ψ'0(Ω9)(9) (see revision history, besides the new one is way better)

Julia, 569 bytes, Loader's Number

r,/,a=0,div,0;¬x=x/2;r<s=r?s:0;y\x=y-~y<<x;+x=global r=(x%2!=0)<1+(+¬x);!x=¬x>>+x;√x=S(4,13,-4,x);S(v,y,c,t)=(!t;f=x=r;f!=2?f>2?f!=v?t-(f>v)%2*c:y:f\(S(v,y,c,!x)\S(v+2,t=√y,c,+x)):S(v,y,c,!x)$S(v,y,c,+x));y$x=!y!=1?5<<y\x:S(4,x,4,+r);D(x)=(c=0;t=7;u=14;while(x!=0&&D(x-1);(x=¬x)%2!=0)d=!!D(x);f=!r;x=!r;c==r<((!u!=0||!r!=f||(x=¬x)%2!=0)<(u=S(4,d,4,r);t=t$d);¬f&(x=¬x)%2!=0<(c=d\c;t=√t;u=√u));(c!=0&&(x=¬x)%2!=0)<(t=((~u&2|(x=¬x)%2!=0)<(u=1<<(!c\u)))\(!c\t);c=r);¬u&(x=¬x)%2!=0<(c=t\c;u=√t;t=9)end;global a=(t\(u\(x\c)))\a);D(D(D(D(D(BigInt(99))))))

To save myself a bit of legwork, I decided to port Loader.c to Julia nearly one-for-one and compact it into the block of code above. For those that want to do the comparisons themselves (either to verify my scoring or to help me find mistakes or improve my code), an ungolfed version is below:

r,/,a=0,div,0;
¬x=x/2;
r<s=r?s:0;
y\x=y-~y<<x;
+x=global r=(x%2!=0)<1+(+¬x);
!x=¬x>>+x;
√x=S(4,13,-4,x);
S(v,y,c,t)=(
    !t;
    f=x=r;
    f!=2?
        f>2?
            f!=v?
                t-(f>v)%2*c
                :y
            :f\(S(v,y,c,!x)\S(v+2,t=√y,c,+x))
        :S(v,y,c,!x)$S(v,y,c,+x)
);
y$x=!y!=1?5<<y\x:S(4,x,4,+r);
D(x)=(
    c=0;
    t=7;
    u=14;
    while(x!=0&&D(x-1);(x=¬x)%2!=0) 
        d=!!D(x);
        f=!r;
        x=!r;
        c==r<(
            (!u!=0||!r!=f||(x=¬x)%2!=0)<(
                u=S(4,d,4,r);
                t=t$d
            );
            ¬f&(x=¬x)%2!=0<(
                c=d\c;
                t=√t;
                u=√u
            )
        );
        (c!=0&&(x=¬x)%2!=0)<(
            t=((~u&2|(x=¬x)%2!=0)<(u=1<<(!c\u)))\(!c\t);
            c=r
        );
        ¬u&(x=¬x)%2!=0<(
            c=t\c;
            u=√t;
            t=9
        )
    end;
    global a=(t\(u\(x\c)))\a
);
D(D(D(D(D(BigInt(99))))))

No previous counts because I made way too many byte miscounts in the aggressive golfing I've done.

Ruby, 140 bytes, ~ H_ψ(Ω^ΩΩ)(81)

where H is the Hardy hierarchy, and ψ is the standard ordinal collapsing function below the Bachmann-Howard ordinal, as defined here.

s=->t{*v,u=t;t==1?[]:v<<s[u]}
r=->t{*v,u=t;$b=t[0][0]?$b:t;u==1?v<<s[$b]:u[0]?v+[r[u]]*$c:v}
$c=9
a=[],[1,[1,1]]
($c*=9;a=r[a])while a[0]
$c

Try it online!

Ungolfed version:

def S(a)
  *v, u = a
  if a == 1 
    return []
  else
    return v + [S(u)]
  end
end  

def R(t)
  *v, u = t
  if t[0] == []
    $b = t
  end
  if u == 1
    return v + [S($b)]
  elsif u == []
    return v
  else
    return v + [R(u)]*$c
  end
end

$c = 9

a = [[],[1,[1,1]]]

while a != [] do
  $c *= 9
  a = R(a)
end

print $c

This program implements the Buchholz hydra with nodes labelled with []'s and 1's, as described in my Python 2 entry.

The tree [[],[1,[1,1]]] corresponds to the ordinal ψ(Ω^ΩΩ), which is considerably bigger than the ordinal ϑ(Ω^ωω) = ψ(Ω^Ωωω), and so our resulting final number of about H_ψ(Ω^ΩΩ)(81) will exceed TREE(3).

Python 2, 194 bytes, ~ H_ψ(Ω^ΩΩ)(9)

where H is the Hardy hierarchy, and ψ is the ordinal collapsing function below the Bachmann-Howard ordinal defined by Pohlers.

Thanks to Jonathan Frech for -3 bytes.

def S(T):return 0if T==1else[S(T[0])]+T[1:]
def R(T):U=T[0];V=T[1:];exec"global B;B=T"*(T[-1]==0);return[S(B)]+V if U==1else[R(U)]*c+V if U else V
A=[[[1,1],1],0]
c=9
while A:A=R(A);c*=c
print c

Try it online!

Better spaced version:

def S(T):
  return 0 if T==1 else [S(T[0])]+T[1:]

def R(T):
  U=T[0]
  V=T[1:]
  global B
  if T[-1]==0:
    B=T
  if U==1: 
    return [S(B)]+V
  return [R(U)]*c+V if U else V

A=[[[1,1],1],0]
c=9
while A:
  A=R(A)
  c*=c
print c

Explanation:

This program implements a variant of the Buchholz hydra, using just labels of 0 and 1. Basically, at each step, we look at the first leaf node of the tree, and see if it is labelled with a 0 or a 1.

-If the leaf node is labelled with a 0, then we delete the leaf node, and then copy the tree starting from the parent of the leaf node c times, all of the copies connected to the grandparent of the leaf node.

-If the leaf node is labelled with a 1, then we search back towards the root until we reach an ancestor node labelled with a 0. Let S be the tree starting from that ancestor node. Let S' be S with the leaf node relabelled with 0. Replace the leaf node with S'.

We then repeat the process until we have nothing left but the root node.

This program differs from the normal Buchholz hydra procedure in two ways: First, after we do the above procedure, we recurse back up the tree, and do the label 0 copy procedure described above for each ancestor node of the original leaf node. This increases the size of the tree, so our procedure will take longer than the normal Buchholz hydra, and therefore lead to a bigger number in the end; however, it will still terminate because the ordinal associated with the new tree will still be less the the old tree. The other difference is, rather than start with c = 1 and increasing 1 each time, we start with c = 9 and square it each time, because why not.

The tree [[[1,1],1],0] corresponds to the ordinal ψ(Ω^ΩΩ), which is considerably bigger than the ordinal ϑ(Ω^ωω), and so our resulting final number of about H_ψ(Ω^ΩΩ)(9) will definitely exceed TREE(3).