K. C. Au's Discovering hypergeometric series with harmonic numbers via Wilf-Zeilberger seeds (arXiv:2602.08721v3) conjectures dimension formulas for the coefficient spaces of WZ-seed summands. For a summand S(a,b,c,d,e;n), V_N(S) is the \mathbb Q-span of the sequences n\mapsto[a^{i_1}b^{i_2}c^{i_3}d^{i_4}e^{i_5}]S with \sum i_j=N, taken before summation.

Claim (Conjecture 3.4(a)). For the first summand of Example II,

\sum_{N\ge0}\dim_{\mathbb Q}V_N(S_1)\,t^N=\frac{1}{(1-t)^3(1-t^2)^2}-\frac{t^2}{1-t},

so the dimensions are 1,3,7,15,29,49,79,119,174,244,\dots

Full proof (8 pages) and replayable programs: doi:10.5281/zenodo.22733350

Proof outline

  1. Sign symmetries. With T=(a-b+e)/2, u=(a-c+2e)/2, v=(2a-b-c-2d+e)/2, x=(c-a)/2, y=(b+c+e)/2, the whole summand (prefactor included) is even in x and in y, so its coefficients live in \mathbb Q[T,u,v,X,Y] with X=x^2, Y=y^2 of weight 2.
  2. Harmonic form. S_1/s_n=C\exp\big(\sum_j p_j^+H_{j,1}(n)+p_j^-H_{j,1/2}(n)\big), with p_j^\pm power sums of signed slopes and C a rational prefactor. Schneider's criterion (arXiv:0808.2596, Thm 6.1) makes the H_{j,\alpha} algebraically independent over \mathbb Q(n), so the coefficient row space M is a module over the algebra generated by the p_j^\pm.
  3. Weighted coordinates h=T-u+v, k=2T-u-v, z=u-v and two weight-2 coordinates r,s. Then p_1^\pm=\pm h+k, p_2^+=r-hk-h^2/2, p_2^-=s-hk+h^2/2, and three generators U=q_3^+, V=q_3^-, W=q_4^+ with L_0=\mathbb Q[h,k,r,s][U,V,W] give, via a 3\times3 certificate of determinant 35,
    L_0+zL_0=\mathcal K:=\{f:[z^2]f(h,0,z,0,0)=0\}.
  4. The missing monomial is real. On k=r=s=0 every p_j^\pm has no z or z^2 term, and the full prefactor has [z]R=-10/w^3 but [z^2]R=0 (w=n+h). So M\subseteq\mathcal K. The degree-1 part of C(1) supplies z, so \mathcal K\subseteq M.

Hence M=\mathcal K: every weighted-degree-N coefficient except h^{N-2}z^2, which gives the series.

Not claimed. Nothing about Conjecture 3.4(b) or the summed numerical spaces, and nothing about the independence of \zeta(3)^2 and \pi^6. Not independently refereed.

Checks run here

The ranks, computed straight from the Pochhammer product in the raw parameters a,\dots,e, with no invariant coordinates used:

python
# Coefficient ranks of Au's Example II summand, straight from the Pochhammer product in a,b,c,d,e.
# Predicted dim V_N: coefficients of 1/((1-t)^3(1-t^2)^2) - t^2/(1-t), i.e. 1,3,7,15,29,49,79.
from itertools import product
from math import comb
p, D, ROWS = 1000003, 6, 230
mons = [m for d in range(D + 1) for m in product(range(d + 1), repeat=5) if sum(m) == d]
idx = {m: i for i, m in enumerate(mons)}
E = [tuple(int(i == j) for j in range(5)) for i in range(5)]
up = [[idx.get(tuple(x + y for x, y in zip(m, E[i]))) for m in mons] for i in range(5)]
deg = [sum(m) for m in mons]
def mul(v, c0, cs):
    out = [x * c0 % p for x in v]
    for i, ci in enumerate(cs):
        if ci:
            for j, t in enumerate(up[i]):
                if t is not None and v[j]:
                    out[t] = (out[t] + ci * v[j]) % p
    return out
def div(v, c0, cs):
    inv, out = pow(c0 % p, -1, p), [0] * len(mons)
    for j in range(len(mons)):  # monomials are sorted by degree, so predecessors are done
        acc = v[j]
        m = mons[j]
        for i, ci in enumerate(cs):
            if ci and m[i]:
                acc -= ci * out[idx[tuple(x - (k == i) for k, x in enumerate(m))]]
        out[j] = acc % p * inv % p
    return out
PTERMS = [  # (exponents of a,b,c,d,e ; power of n ; coefficient)
    ('00000', 2, -3), ('00000', 3, 10), ('00001', 1, -4), ('00001', 2, 18), ('00010', 2, -6),
    ('00100', 1, 1), ('00100', 2, -6), ('01000', 1, 2), ('01000', 2, -12), ('10000', 1, -3),
    ('10000', 2, 18), ('00002', 0, -1), ('00002', 1, 10), ('00011', 1, -8), ('00101', 0, 1),
    ('00101', 1, -8), ('00110', 1, 2), ('01001', 0, 2), ('01001', 1, -16), ('01010', 1, 4),
    ('01100', 1, 4), ('02000', 1, 4), ('10001', 0, -2), ('10001', 1, 22), ('10010', 1, -6),
    ('10100', 0, 1), ('10100', 1, -6), ('11000', 0, 1), ('11000', 1, -14), ('20000', 0, -1),
    ('20000', 1, 10), ('00003', 0, 2), ('00012', 0, -2), ('00102', 0, -2), ('00111', 0, 2),
    ('01002', 0, -4), ('01011', 0, 4), ('01101', 0, 4), ('02001', 0, 4), ('10002', 0, 6),
    ('10011', 0, -4), ('10101', 0, -4), ('10110', 0, 2), ('11001', 0, -10), ('11010', 0, 2),
    ('11100', 0, 2), ('12000', 0, 2), ('20001', 0, 6), ('20010', 0, -2), ('20100', 0, -2),
    ('21000', 0, -4), ('30000', 0, 2),
]
def P(n):
    v = [0] * len(mons)
    for ex, k, co in PTERMS:
        v[idx[tuple(map(int, ex))]] = (v[idx[tuple(map(int, ex))]] + co * n ** k) % p
    return v
def mulv(v, w):
    out = [0] * len(mons)
    for j, x in enumerate(w):
        if x:
            for i, y in enumerate(v):
                if y and deg[i] + deg[j] <= D:
                    k = idx[tuple(s + t for s, t in zip(mons[i], mons[j]))]
                    out[k] = (out[k] + x * y) % p
    return out
G = [1] + [0] * (len(mons) - 1)
G = mul(G, 1, [0, 0, 0, 2, 2])
for c0, cs in [(1, [0, 0, 0, 0, 1]), (1, [1, 0, -1, 0, 1]), (1, [2, 0, 0, -2, 2]), (1, [2, -2, -2, -2, 0]), (1, [1, -1, 0, 0, 1]), (2, [1, -1, 0, 0, 1])]:
    G = div(G, c0, cs)
rows = []
for n in range(1, ROWS + 1):
    rows.append(mulv(G, P(n)))
    for c0, cs in [(n, [0, -1, 0, 0, 0]), (n, [1, 0, 0, 0, 1]), (n, [0, 0, 1, 0, 1]), (n, [1, -1, -1, 0, 0]), (2 * n + 1, [0, 0, 0, 2, 2]), (2 * n - 1, [2, -2, 0, -2, 0])]:
        G = mul(G, c0, cs)
    for c0, cs in [(n + 1, [0, 0, 0, 0, 1]), (n + 1, [1, 0, -1, 0, 1]), (2 * n + 1, [2, 0, 0, -2, 2]), (2 * n + 1, [2, -2, -2, -2, 0]), (2 * n + 1, [1, -1, 0, 0, 1]), (2 * n + 2, [1, -1, 0, 0, 1])]:
        G = div(G, c0, cs)
def rank(M):
    M, r = [row[:] for row in M], 0
    for col in range(len(M[0])):
        piv = next((i for i in range(r, len(M)) if M[i][col]), None)
        if piv is None: continue
        M[r], M[piv] = M[piv], M[r]
        inv = pow(M[r][col], -1, p)
        M[r] = [x * inv % p for x in M[r]]
        for i in range(len(M)):
            if i != r and M[i][col]:
                f = M[i][col]; M[i] = [(x - f * y) % p for x, y in zip(M[i], M[r])]
        r += 1
    return r
for N in range(D + 1):
    cols = [j for j in range(len(mons)) if deg[j] == N]
    want = sum(comb(N - 2 * j + 2, 2) * (j + 1) for j in range(N // 2 + 1)) - (N >= 2)
    got = rank([[row[j] for j in cols] for row in rows])
    print(f"N={N}: {len(cols)} monomials, rank {got}, predicted {want}", "OK" if got == want else "MISMATCH")
ran3753 ms
N=0: 1 monomials, rank 1, predicted 1 OK
N=1: 5 monomials, rank 3, predicted 3 OK
N=2: 15 monomials, rank 7, predicted 7 OK
N=3: 35 monomials, rank 15, predicted 15 OK
N=4: 70 monomials, rank 29, predicted 29 OK
N=5: 126 monomials, rank 49, predicted 49 OK
N=6: 210 monomials, rank 79, predicted 79 OK

The algebra behind steps 3 and 4: generators from the slope table, the determinant, and the prefactor jet.

python
# The algebra the proof rests on: generators U,V,W from the slope table, the det-35 certificate,
# and the vanishing z^2 jet of the full prefactor on k=r=s=0.
import sympy as S
T,u,v,x,y,h,k,z,r,s,n,w = S.symbols('T u v x y h k z r s n w'); Q = S.Rational
slopes = {'+': ([T+x-y, T-x+y, T+x+y, T-x-y], [u+x, u-x, T]), '-': ([T+u-v, T-u+v], [v+y, v-y, T])}
X = -r-2*s-Q(7,4)*z**2-Q(5,2)*k*z+Q(3,4)*k**2
Y = s+Q(5,4)*z**2+k*z/2-k**2/4
def p(j, side):
    pos, neg = slopes[side]
    f = S.Poly(S.expand(Q((-1)**(j-1), j)*(sum(t**j for t in pos)-sum(t**j for t in neg))), x, y)
    g = sum(co*X**(i//2)*Y**(l//2) for (i, l), co in f.terms() if i % 2 == 0 and l % 2 == 0)
    assert all(i % 2 == 0 and l % 2 == 0 for (i, l), _ in f.terms())
    return S.expand(g.subs({T: h+z, u: h+(3*z-k)/2, v: h+(z-k)/2}, simultaneous=True))
cp, al = Q(5,6)*k**3-k*r-2*k*s, -Q(7,2)*k**2-r+2*s
cm, be = -k**3/6+k*s, k**2/2-s
c4 = Q(11,8)*k**4-Q(3,2)*k**2*r-7*k**2*s-r**2/2+4*r*s+9*s**2
ga, de = -Q(19,2)*k**3+5*k*r+24*k*s, Q(23,2)*k**2+5*r+12*s
checks = {
 'p1+ = h+k': p(1,'+')-(h+k), 'p1- = -h+k': p(1,'-')-(-h+k),
 'p2+ = r-hk-h^2/2': p(2,'+')-(r-h*k-h**2/2), 'p2- = s-hk+h^2/2': p(2,'-')-(s-h*k+h**2/2),
 'U = q3+': p(3,'+').subs(h,0)-(cp+al*z+2*z**3),
 'V = q3-': p(3,'-').subs(h,0)-(cm+be*z+k*z**2+z**3),
 'W = q4+': p(4,'+').subs(h,0)-(c4+ga*z+de*z**2+18*k*z**3+3*z**4),
}
for name, diff in checks.items():
    print(name, S.expand(diff) == 0)
print('det =', S.Matrix([[3,5,12],[2,-1,2],[1,0,-1]]).det())
Xi, Yi = S.symbols('Xi Yi')  # X, Y as independent invariants for the prefactor
Pt = (2*T**3+10*T**2*n+2*T**2*u+2*T**2*v-T**2-2*T*(Xi+Yi)+18*T*n**2+8*T*n*(u+v)-4*T*n+4*T*u*v-2*T*u
      -2*Xi*n-2*Xi*u+2*Xi*v-Xi-2*Yi*n+2*Yi*u-2*Yi*v+Yi+10*n**3+6*n**2*(u+v)-3*n**2+4*n*u*v-2*n*u)
Dt = (((n+T)**2-Xi-Yi)**2-4*Xi*Yi)*(n-Q(1,2)+T-u+v)
surf = {T: h+z, u: h+3*z/2, v: h+z/2, Xi: -Q(7,4)*z**2, Yi: Q(5,4)*z**2}
num = 2*(w**2*(10*w-3)+(30*w**2-7*w)*z+(30*w-1)*z**2+16*z**3)
den = (2*w-1)*(w**4+4*w**3*z+7*w**2*z**2+6*w*z**3+11*z**4)
print('prefactor on k=r=s=0 matches eq. (21):', S.simplify((Pt/Dt).subs(surf, simultaneous=True) - (num/den).subs(w, n+h)) == 0)
ser = S.series(num/den, z, 0, 3).removeO()
print('[z]R =', S.simplify(ser.coeff(z, 1)), '  [z^2]R =', S.simplify(ser.coeff(z, 2)))
ran1309 ms
p1+ = h+k True
p1- = -h+k True
p2+ = r-hk-h^2/2 True
p2- = s-hk+h^2/2 True
U = q3+ True
V = q3- True
W = q4+ True
det = 35
prefactor on k=r=s=0 matches eq. (21): True
[z]R = -10/w**3   [z^2]R = 0