Quarter-plane lattice paths with interacting boundaries: Kreweras and friends
Ruijie Xu, Nicholas R. Beaton and Aleksander L. Owczarek

Some calculations accompanying the solution to reverse Kreweras walks with general boundary weights (a,b,c). Symbols and equation numbers match the FPSAC submission where possible.

Many derivations are long and tedious so in general only final expressions are given here.

Preliminaries
It will be useful to have some series to substitute into equations to check their validity.

(* shorthand to apply a function f to the terms of a series *)
ApplyToSeries[f_,S_]:=MapAt[f/@#&,S,3]

(* the step set generator for reverse Kreweras walks and the associated kernel *)
S[x_,y_]:=x+y+1/x/y
A[x_]:=Coefficient[S[x,y],y,-1];
B[y_]:=Coefficient[S[x,y],x,-1];
K[x_,y_]:=1-t S[x,y]

(* the 8 coefficients of S, used for iteration *)
epsmm=Coefficient[Coefficient[S[x,y],x,-1],y,-1];
epsm0=Coefficient[Coefficient[S[x,y],x,-1],y,0];
epsmp=Coefficient[Coefficient[S[x,y],x,-1],y,1];
eps0m=Coefficient[Coefficient[S[x,y],x,0],y,-1];
eps0p=Coefficient[Coefficient[S[x,y],x,0],y,1];
epspm=Coefficient[Coefficient[S[x,y],x,1],y,-1];
epsp0=Coefficient[Coefficient[S[x,y],x,1],y,0];
epspp=Coefficient[Coefficient[S[x,y],x,1],y,1];
ϵ=epsmm;

(* the generating function with a given number of terms *)
QQ[0]=1+O[t];
QQ[n_]:=QQ[n]=ApplyToSeries[Expand@*Simplify,QQeqn[n-1]]
(* coefficients of specific powers of x, y or both *)
QQcx[n_,i_]:=QQcx[n,i]=ApplyToSeries[Coefficient[#,x,i]&,QQ[n]]
QQcy[n_,j_]:=QQcy[n,j]=ApplyToSeries[Coefficient[#,y,j]&,QQ[n]]
QQcxy[n_,i_,j_]:=QQcxy[n,i,j]=ApplyToSeries[Coefficient[Coefficient[#,x,i],y,j]&,QQ[n]]
(* evaluating QQ at some other values of (x,y) *)
QQeval[n_,xx_,yy_]:=QQeval[n,xx,yy]=QQ[n]/.{x->xx,y->yy}
QQcxeval[n_,i_,yy_]:=QQcxeval[n,i,yy]=QQcx[n,i]/.{y->yy}
QQcyeval[n_,j_,xx_]:=QQcyeval[n,j,xx]=QQcy[n,j]/.{x->xx}

(* the generalised diagonal will also be needed, as per (2.2) *)
(* the (x^π)+(x^(2π)) is just to accommodate the way Select works on expressions with only one term *)
QQdj[n_,j_]:=QQdj[n,j]=ApplyToSeries[Select[#+(x^π)+(x^(2π)),Exponent[#,y]==Exponent[#,x]+j&]&,QQ[n]]/.{y->1}
QQdjeval[n_,j_,xx_]:=QQdjeval[n,j,xx]=QQdj[n,j]/.{x->xx}

(* the functional equation which can be iterated *)
QQeqn[n_]:=1+t S[x,y] QQ[n]-t/y A[x] QQcy[n,0]-t/x B[y] QQcx[n,0]+t/x/y epsmm QQcxy[n,0,0]+t(a-1)epsp0 x QQcy[n,0]+t(a-1)epsm0/x(QQcy[n,0]-QQcxy[n,0,0]-x QQcxy[n,1,0])+t(a-1)epspm x QQcy[n,1]+t(a-1)eps0m(QQcy[n,1]-QQcxy[n,0,1])+t(a-1)epsmm/x(QQcy[n,1]-QQcxy[n,0,1]-x QQcxy[n,1,1])+t(b-1)eps0p y QQcx[n,0]+t(b-1)eps0m/y(QQcx[n,0]-QQcxy[n,0,0]-y QQcxy[n,0,1])+t(b-1)epsmp y QQcx[n,1]+t(b-1)epsm0(QQcx[n,1]-QQcxy[n,1,0])+t(b-1)epsmm/y(QQcx[n,1]-QQcxy[n,1,0]-y QQcxy[n,1,1])+t(c-1)epsm0 QQcxy[n,1,0]+t(c-1)eps0m QQcxy[n,0,1]+t(c-1)epsmm QQcxy[n,1,1]

(* the first few terms of Q[0,0] *)
QQcxy[12,0,0]

(* this will also be useful *)
Needs["Notation`"]

Symbolize[ParsedBoxWrapper[SubscriptBox["_","_"]]]
Symbolize[ParsedBoxWrapper[SubsuperscriptBox["_","_","_"]]]

Section 2.7

(* verifying the main functional equation (2.7) *)
-x y K[x,y]Q[x,y]+x y/c+x(y-y/a-t A[x])Q[x,0]+y(x-x/b-t B[y])Q[0,y]+(x y(1/a+1/b-1/c-1)+t ϵ)Q[0,0];
%/.{Q[x_,y_]->QQeval[12,x,y]}

Section 3.1

(* verifying the full-orbit sum (3.3) *)
y(1-b+t x b)(1-a+t y a)Q[x,y]-1/x(1-a+t y a)(t b+x y-x y b)Q[1/x/y,y]+1/x(1-b+t y b)(t a+x y-x y a)Q[y,1/x/y]-y(1-a+t x a)(1-b+t y b)Q[y,x]+1/x(1-a+t x a)(t b+x y-x y b)Q[1/x/y,x]-1/x(1-b+t x b)(t a+x y-x y a)Q[x,1/x/y];
%/.{Q[x_,y_]->QQeval[12,x,y]}

(* the constant term of the full orbit sum with respect to y *)
(* this was omitted from the paper *)
-((a b (t^2)+x-b x-2 a t (x^2)+b t (x^2)+a b t (x^2))/(c t x))+(1/(a b c t x)) ((a^2) (b^2) (t^2)-(a^2) (b^2) c (t^2)+a b x-a (b^2) x+a c x-(a^2) c x-b c x-a b c x+(a^2) b c x+(b^2) c x+(a^2) (b^2) c (t^3) x-2 (a^2) b t (x^2)+a (b^2) t (x^2)+(a^2) (b^2) t (x^2)+(a^2) c t (x^2)+a b c t (x^2)-(b^2) c t (x^2)-(a^2) (b^2) c t (x^2)-(a^2) b c (t^2) (x^3)+(a^2) (b^2) c (t^2) (x^3))Q[0,0]-(((1-a+a t x) ((b^2) (t^2)+x-b x-b t (x^2)+(b^2) t (x^2)))/(b t x))Q[0,x]+(((1-b+b t x) ((a^2) (t^2)+x-a x-a t (x^2)+(a^2) t (x^2)))/(a t x))Q[x,0]-(((2 a (t^2)-2 b (t^2)-a b (t^2)-x+b x+a b (t^3) x+t (x^2)+a t (x^2)-b t (x^2)-a b t (x^2)-a (t^2) (x^3)+a b (t^2) (x^3)))/(t x))(Subscript[Q, 0]^d)[(1/x)]-(((2 a b (t^2)+2 x-a x-b x-a t (x^2)-b t (x^2)+2 a b t (x^2)))/(x^2))(Subscript[Q, 1]^d)[(1/x)];
%/.{Q[x_,y_]->QQeval[12,x,y],(Subscript[Q, 0]^d)[x_]->QQdjeval[12,0,x],(Subscript[Q, 1]^d)[x_]->QQdjeval[12,1,x]}

(* verifying equation (3.4) *)
1/x/a/b/c(a-b)(t^2a b c-x c(1-a-b+a b)+t x^2(a b-a c-b c+a b c))Q[0,0]-t^2a(b-1)(Subscript[Q, 1,0])-1/x/a(1-b+t b x)(t^2a^2+x(1-a)(1-t x a))Q[x,0]+t^2(a-1)b (Subscript[Q, 0,1])+1/x/b(1-a+t x a)(t^2b^2+x(1-b)(1-t x b))Q[0,x]-t x(a-b)/c;
%/.{Q[x_,y_]->QQeval[12,x,y],(Subscript[Q, 0,1])->QQcxy[12,0,1],(Subscript[Q, 1,0])->QQcxy[12,1,0]}

(* verifying equation (3.5) *)
1/x/c(t^2(a b+a c-b c-a b c)+x(1-b-c+b c)+t^3x a b c-t x^2(a -a b-c-a c+b c+a b c)+t^2x^3a c(b-1))Q[0,0]-1/x(t^2(2 a-2 b-a b)+t^3x a b+t x^2(1+a-b-a b)+x(1+t^2x^2a)(b-1))(Subscript[Q, 0]^d)[1/x]-t/x^2(2t^2a b+x(2-a-b)-t x^2(a+b-2 a b))(Subscript[Q, 1]^d)[1/x]+t^2(a-1)b (Subscript[Q, 0,1])-t^2a(b-1)(Subscript[Q, 1,0])-1/x/c(t^2a b-x(1-t x a)(b-1));
%/.{Q[x_,y_]->QQeval[12,x,y],(Subscript[Q, 0]^d)[x_]->QQdjeval[12,0,x],(Subscript[Q, 1]^d)[x_]->QQdjeval[12,1,x],(Subscript[Q, 1,0])->QQcxy[12,1,0],(Subscript[Q, 0,1])->QQcxy[12,0,1]}

Section 3.2

(* the coefficients in equation (3.6) are *)
(Subscript[α, x,y])=x (1-b+b t x) (1-a+a t y);
(Subscript[α, 1/x/y,y])=(((1-a+a t y) (-b t-x y+b x y))/y);
(Subscript[α, 1/x/y,x])=(((1-a+a t x) (b t+x y-b x y))/y);
α=((a b (t^2) x-a b (t^2) y+x y-a x y-b x y+a b x y+a t (x^2) y+b t (x^2) y-2 a b t (x^2) y+a b (t^2) (x^2) (y^2))/(c y));
(Subscript[α, x,0])=-(((1-b+b t x) (1-a+a t y) (a t+x y-a x y))/(a y));
(Subscript[α, 0,x])=-(((1-a+a t x) (1-b+b t y) (b t+x y-b x y))/(b y));
(Subscript[α, 0,0])=-(1/(a b c y))(-a b c t+(a^2) b c t+a (b^2) c t-(a^2) (b^2) c t+(a^2) (b^2) (t^2) x-(a^2) b c (t^2) x-a (b^2) c (t^2) x+(a^2) (b^2) c (t^2) x-(a^2) (b^2) (t^2) y+(a^2) (b^2) c (t^2) y+a b x y-(a^2) b x y-a (b^2) x y+(a^2) (b^2) x y-a c x y+(a^2) c x y-b c x y+3 a b c x y-2 (a^2) b c x y+(b^2) c x y-2 a (b^2) c x y+(a^2) (b^2) c x y-(a^2) (b^2) c (t^3) x y+(a^2) b t (x^2) y+a (b^2) t (x^2) y-2 (a^2) (b^2) t (x^2) y-(a^2) c t (x^2) y-a b c t (x^2) y+2 (a^2) b c t (x^2) y-(b^2) c t (x^2) y+2 a (b^2) c t (x^2) y-(a^2) (b^2) c t (x^2) y-a b c t x (y^2)+(a^2) b c t x (y^2)+a (b^2) c t x (y^2)-(a^2) (b^2) c t x (y^2)+(a^2) (b^2) (t^2) (x^2) (y^2)-(a^2) b c (t^2) (x^2) (y^2)-a (b^2) c (t^2) (x^2) (y^2)+(a^2) (b^2) c (t^2) (x^2) (y^2));

(* then check equation (3.6) *)
-K[x,y]((Subscript[α, x,y])Q[x,y]+(Subscript[α, 1/x/y,y])Q[1/x/y,y]+(Subscript[α, 1/x/y,x])Q[1/x/y,x])+α+(Subscript[α, x,0])Q[x,0]+(Subscript[α, 0,x])Q[0,x]+(Subscript[α, 0,0])Q[0,0];
%/.{Q[x_,y_]->QQeval[12,x,y]}

(* the roots of K *)
Δ=(1-t x)^2-4t^2/x;
(SubPlus[γ])=(1-t x+Sqrt[Δ])/(2t);
(SubMinus[γ])=(1-t x-Sqrt[Δ])/(2t);
{K[x,(SubMinus[γ])],K[x,(SubPlus[γ])]}//FullSimplify
ApplyToSeries[Expand@*PowerExpand,Series[(SubPlus[γ]),{t,0,2}]]
ApplyToSeries[Expand@*PowerExpand,Series[(SubMinus[γ]),{t,0,2}]]

(* the partial fraction expansion of 1/K in equation (3.7) *)
(Subscript[A, 1])=-1/Sqrt[Δ];
(Subscript[A, 2])=1/Sqrt[Δ];
(Subscript[A, 3])=1/Sqrt[Δ];
1/K[x,y]-((Subscript[A, 1])+(Subscript[A, 2])/(1-y/(SubPlus[γ]))+(Subscript[A, 3])/(1-(SubMinus[γ])/y))//FullSimplify
(* and then *)
ApplyToSeries[Expand@*PowerExpand,Series[(Subscript[A, 1]),{t,0,2}]]
ApplyToSeries[Expand@*PowerExpand,Series[(Subscript[A, 2])/(1-y/(SubPlus[γ])),{t,0,2}]]
ApplyToSeries[Expand@*PowerExpand,Series[(Subscript[A, 3])/(1-(SubMinus[γ])/y),{t,0,2}]]

(* the coefficients in equation (3.9) *)
(* these are messy but we will tidy things up for equation (3.10) *)
β=-(((-1+a) x (1+b (-1+t x))+((2 (t^2) (x (-1+b-b t x)+a (x-t (x^2)+b (-x+t (x^2)+(t^2) (1+(x^3))))))/((Sqrt[x-2 t (x^2)+(t^2) (-4+(x^3))]) ((Sqrt[x])-t (x^(3/2))+(Sqrt[x-2 t (x^2)+(t^2) (-4+(x^3))])))))/(c t));
(Subscript[β, x,0])=-(((1+b (-1+t x)) (x-a x (1+t x)+(a^2) t (t+(x^2))) (x+(t^2) (-2+(x^3))+(Sqrt[x]) (Sqrt[x-2 t (x^2)+(t^2) (-4+(x^3))])-t (2 (x^2)+(x^(3/2)) (Sqrt[x-2 t (x^2)+(t^2) (-4+(x^3))]))))/(a t (Sqrt[x-2 t (x^2)+(t^2) (-4+(x^3))]) ((Sqrt[x])-t (x^(3/2))+(Sqrt[x-2 t (x^2)+(t^2) (-4+(x^3))]))));
(Subscript[β, 0,x])=-((2 t (1+a (-1+t x)) (x-b x (1+t x)+(b^2) t (t+(x^2))))/(b (Sqrt[x-2 t (x^2)+(t^2) (-4+(x^3))]) ((Sqrt[x])-t (x^(3/2))+(Sqrt[x-2 t (x^2)+(t^2) (-4+(x^3))]))));
(Subscript[β, 0]^d)=(a-b) t;
(Subscript[β, 1]^d)=1-b t x+a (-1+((b t (t+(x^2)))/x));
(Subscript[β, 0,0])=(1/(a c t))((-1+a) (a-c) x (1+b (-1+t x))+(2 (t^2) (b c x (1+b (-1+t x))-a x (1+b (-1+t x)) (-c+b (1+c+c t x))+(a^2) (c x (-1+t x)-b x (-1+t x) (1+c+c t x)+(b^2) (-x+c (t^3) x-(-1+c) t (x^2)+(t^2) (1+(x^3)+c (-1+(x^3)))))))/(b (Sqrt[x-2 t (x^2)+(t^2) (-4+(x^3))]) ((Sqrt[x])-t (x^(3/2))+(Sqrt[x-2 t (x^2)+(t^2) (-4+(x^3))]))));

(* verifying equation (3.9) *)
β+(Subscript[β, x,0])Q[x,0]+(Subscript[β, 0,x])Q[0,x]+(Subscript[β, 0]^d)(Subscript[Q, 0]^d)[1/x]+(Subscript[β, 1]^d)(Subscript[Q, 1]^d)[1/x]+(Subscript[β, 0,0])Q[0,0];
%/.{Q[x_,y_]->QQeval[12,x,y],(Subscript[Q, 0]^d)[x_]->QQdjeval[12,0,x],(Subscript[Q, 1]^d)[x_]->QQdjeval[12,1,x]}

(* the coefficients in equation (3.10) *)
(Subscript[μ, x,0])=-(((-1+t x) (1-b+b t x) ((a^2) (t^2)+x-a x-a t (x^2)+(a^2) t (x^2))(2 a b (t^2)+2 x-a x-b x-a t (x^2)-b t (x^2)+2 a b t (x^2)))/(a t));
μ=(((-1+t x) (2 a b (t^2)+2 x-a x-b x-a t (x^2)-b t (x^2)+2 a b t (x^2)) (a b (t^2)-x+a x+b x-a b x-2 b t (x^2)+a b t (x^2)+a b (t^2) (x^3)))/(2 c t));
ν=(1/(2 c t))x (-1+t x) ((a^2) b (t^2)+a (b^2) (t^2)-2 (a^2) (b^2) (t^2)+a x-(a^2) x+b x-2 a b x+(a^2) b x-(b^2) x+a (b^2) x+2 (a^2) (b^2) (t^3) x+(a^2) b t (x^2)+a (b^2) t (x^2)-2 (a^2) (b^2) t (x^2)-(a^2) b (t^2) (x^3)-a (b^2) (t^2) (x^3)+2 (a^2) (b^2) (t^2) (x^3));
(Subscript[ν, 0]^d)=-(((-1+t x) ((a^2) (t^2)+x-a x-a t (x^2)+(a^2) t (x^2)) ((b^2) (t^2)+x-b x-b t (x^2)+(b^2) t (x^2)))/t);
(Subscript[μ, 0,0])=-(1/(2 a c t))(-1+t x) (2 a b (t^2)+2 x-a x-b x-a t (x^2)-b t (x^2)+2 a b t (x^2)) ((a^2) b (t^2)+(a^2) c (t^2)-a b c (t^2)-(a^2) b c (t^2)-a x+(a^2) x+a b x-(a^2) b x+2 c x-2 a c x-2 b c x+2 a b c x+(a^2) b c (t^3) x-2 a b t (x^2)+(a^2) b t (x^2)+(a^2) c t (x^2)+2 b c t (x^2)-a b c t (x^2)-(a^2) b c t (x^2)+(a^2) b (t^2) (x^3)-(a^2) c (t^2) (x^3)-a b c (t^2) (x^3)+(a^2) b c (t^2) (x^3));
(Subscript[ν, 0,0])=(1/(2 c t))x (-1+t x) (-(a^2) b (t^2)-a (b^2) (t^2)+2 (a^2) (b^2) (t^2)+(a^2) c (t^2)-(a^2) b c (t^2)+(b^2) c (t^2)-a (b^2) c (t^2)-a x+(a^2) x-b x+2 a b x-(a^2) b x+(b^2) x-a (b^2) x+2 c x-2 a c x-2 b c x+2 a b c x-2 (a^2) (b^2) (t^3) x+(a^2) b c (t^3) x+a (b^2) c (t^3) x-(a^2) b t (x^2)-a (b^2) t (x^2)+2 (a^2) (b^2) t (x^2)+(a^2) c t (x^2)-(a^2) b c t (x^2)+(b^2) c t (x^2)-a (b^2) c t (x^2)+(a^2) b (t^2) (x^3)+a (b^2) (t^2) (x^3)-2 (a^2) (b^2) (t^2) (x^3)-(a^2) c (t^2) (x^3)+(a^2) b c (t^2) (x^3)-(b^2) c (t^2) (x^3)+a (b^2) c (t^2) (x^3));
(Subscript[μ, 0,1])=-(1/2) (-1+a) b t x (-1+t x) (2 a b (t^2)+2 x-a x-b x-a t (x^2)-b t (x^2)+2 a b t (x^2));
(Subscript[ν, 0,1])=(1/2) (-1+a) (a-b) b t (x^2) (-1+t x);
(Subscript[μ, 1,0])=(1/2) a (-1+b) t x (-1+t x) (2 a b (t^2)+2 x-a x-b x-a t (x^2)-b t (x^2)+2 a b t (x^2));
(Subscript[ν, 1,0])=-(1/2) a (a-b) (-1+b) t (x^2) (-1+t x);

(* then verify equation (3.10) *)
-(Subscript[μ, x,0])Q[x,0]+(μ+ν Sqrt[Δ])+(Subscript[ν, 0]^d)Sqrt[Δ](Subscript[Q, 0]^d)[1/x]+((Subscript[μ, 0,0])+(Subscript[ν, 0,0])Sqrt[Δ])Q[0,0]+((Subscript[μ, 0,1])+(Subscript[ν, 0,1])Sqrt[Δ])(Subscript[Q, 0,1])+((Subscript[μ, 1,0])+(Subscript[ν, 1,0])Sqrt[Δ])(Subscript[Q, 1,0]);
%/.{Q[x_,y_]->QQeval[12,x,y],(Subscript[Q, 0]^d)[x_]->QQdjeval[12,0,x],(Subscript[Q, 0,1])->QQcxy[12,0,1],(Subscript[Q, 1,0])->QQcxy[12,1,0]}

(* the factorisation of Δ *)
(* using Root instead of radicals seems to improve performance *)
(* different versions of Mathematica may take Root[...] in different orders, so let's not make any assumptions *)
Off[Root::sbr]
(Subscript[d, 1])=Root[-4 (t^2)+#-2 t ((#)^2)+(t^2) ((#)^3)&,1];
(Subscript[d, 2])=Root[-4 (t^2)+#-2 t ((#)^2)+(t^2) ((#)^3)&,2];
(Subscript[d, 3])=Root[-4 (t^2)+#-2 t ((#)^2)+(t^2) ((#)^3)&,3];
{(Subscript[d, 1]),(Subscript[d, 2]),(Subscript[d, 3])};
Normal[Series[%,{t,0,1}]];
(Subscript[δ, 1])=Select[{(Subscript[d, 1]),(Subscript[d, 2]),(Subscript[d, 3])},Normal[Series[#,{t,0,1}]]==0&][[1]]
(Subscript[δ, 2])=Select[{(Subscript[d, 1]),(Subscript[d, 2]),(Subscript[d, 3])},Normal[Series[#,{t,0,1}]]==1/t+2Sqrt[t]&][[1]]
(Subscript[δ, 3])=Select[{(Subscript[d, 1]),(Subscript[d, 2]),(Subscript[d, 3])},Normal[Series[#,{t,0,1}]]==1/t-2Sqrt[t]&][[1]]

(* verify equation (3.13) *)
Δ-4t^2/(Subscript[δ, 1])(1-(Subscript[δ, 1])/x)(1-x/(Subscript[δ, 2]))(1-x/(Subscript[δ, 3]))//FullSimplify

(* then as per equation (3.14) *)
(Subscript[Δ, 0])=4t^2/(Subscript[δ, 1]);
(SubMinus[Δ])=1-(Subscript[δ, 1])/x;
(SubPlus[Δ])=(1-x/(Subscript[δ, 2]))(1-x/(Subscript[δ, 3]));
(* and verify (3.15) and (3.16) *)
Series[1/Sqrt[(SubPlus[Δ])],{t,0,4}]
Series[Sqrt[(Subscript[Δ, 0])(SubMinus[Δ])],{t,0,5}]

(* the coefficients in equation (3.17) *)
(* these are all pretty messy, but definitely algebraic *)
(Subscript[P, x,0])=((a^2) (t^2)+x-a x-a t (x^2)+(a^2) t (x^2)) (2 a b (t^2)+2 x-a x-b x-a t (x^2)-b t (x^2)+2 a b t (x^2));
σ=(1/(2 (c+b c (-1+t x))))x (Sqrt[(1-(x/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,2]))) (1-(x/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,3])))]) (2 (-1+a) a (-1+b) x (Sqrt[((t^2)/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]))]) (2 b-2 a (-1+b+b t x)+a b t Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1])-(1/4) a x (Sqrt[((t^2)/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]))]) (8 (-(-1+b) b+(a^2) (-1+2 (b^2) t ((t^2)-x+t (x^2))+b (1+t x-(t^2) (x^2)))+a (1-2 b+(b^2) (1+t x-(t^2) (x^2))))-4 a b (-b+a (-1+2 b)) t (-1+t x) Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]+a b (a+b-2 a b) (t^2) ((Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1])^2))+a (-a b (-2+a+b) (t^2)+2 a b (-1+a+b-a b) (t^2)+((2 (a^2) (b^2) (t^4))/x)+(a^2) (b^2) (t^4) ((1/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,2]))+(1/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,3])))-(((-a x (1+t x)-x (-2+b+b t x)+2 a b t (t+(x^2))) (x (-1+b-2 b t x)+a (x+b (-x+t (x^2)+(t^2) (1+(x^3))))))/(x (Sqrt[(((-x+Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,2]) (-x+Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,3]))/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,2] Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,3]))])))));
(Subscript[σ, 0,0])=(1/(2 (c+b c (-1+t x))))x (Sqrt[(1-(x/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,2]))) (1-(x/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,3])))]) (-2 a b (t^2) (-2 (-1+b) c+a (-1+b) (1+2 c)+(a^2) (1+b (-1+c (t^3))))+((2 (a^2) b (a b (-1+c)-a c+b c) (t^4))/x)+(1/4) a x (Sqrt[((t^2)/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]))]) (8 (-(-1+b) (b-2 c)+(a^2) (-1+b+2 (b^2) (t^3)-b c (t^3))+a (1+2 c-2 b (1+c)+(b^2) (1-c (t^3))))-8 ((b^2) c-a (b^2) (1+c)+(a^2) (2 (b^2)+c-b (1+c))) t x-8 (-(b^2) c+a (b^2) (1+c)+(a^2) (b-2 (b^2)-c+b c)) (t^2) (x^2)+4 ((b^2) c-a (b^2) (1+c)+(a^2) (2 (b^2)+c-b (1+c))) t Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]+4 (-(b^2) c+a (b^2) (1+c)+(a^2) (b-2 (b^2)-c+b c)) (t^2) x Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]+(-(b^2) c+a (b^2) (1+c)+(a^2) (b-2 (b^2)-c+b c)) (t^2) ((Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1])^2))+(1/2) a x (Sqrt[((t^2)/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]))]) (8 (-1+a) a (-1+b) b-4 (-1+a) (-1+b) (2 (a+b-a b t x)+a b t Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1])+c (8 ((-1+b) (-1+b t x)+(a^2) t (x+(b^2) t (2 t+(x^2))-b ((t^2)+x+t (x^2)))-a (1+t x-b (((1+t x))^2)+(b^2) t ((t^2)+x+t (x^2))))-4 (-1+a) (-1+b) t (-b+a (-1+b t x)) Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]+a b (-1+a+b-a b) (t^2) ((Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1])^2)))-2 a c (t^2) (2 (-1+b) b-a (2-5 b+3 (b^2))-(((a^2) (-x+b (2 (t^2)+x)+2 (b^2) (t^2) (-1+t x)))/x)+(((a^2) (-1+b) b (t^2))/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,2]))+(((a^2) (-1+b) b (t^2))/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,3])))+((-a x (1+t x)-x (-2+b+b t x)+2 a b t (t+(x^2))) (2 c x (1+b (-1+t x))-a (x+2 c x+b x (-1+2 t x)+b c (-2 x+t (x^2)+(t^2) (1+(x^3))))+(a^2) (x+c t (t+(x^2)-t (x^3))+b (-x+c (t^3) x-(-1+c) t (x^2)+(t^2) (1+(x^3)+c (-1+(x^3)))))))/(x (Sqrt[(((-x+Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,2]) (-x+Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,3]))/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,2] Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,3]))]))+(a (a b (-1+c)-a c+b c) (t^2) (a b (t^2) Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,3]+Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,2] (a b (t^2)-(-2+a+b) Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,3])))/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,2] Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,3]));
(Subscript[σ, 0,1])=(1/(2 (c+b c (-1+t x))))(-1+a) a b c (t^2) x (-2 a b (t^2)-2 (a-b) x (Sqrt[((t^2)/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]))])-4 a (-1+b) x (Sqrt[((t^2)/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]))])+((-a x (1+t x)-x (-2+b+b t x)+2 a b t (t+(x^2)))/(Sqrt[(((x-Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,2]) (x-Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,3]))/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,2] Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,3]))]))) (Sqrt[(1-(x/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,2]))) (1-(x/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,3])))]);
(Subscript[σ, 1,0])=(1/(2 (c+b c (-1+t x))))(a^2) (-1+b) c (t^2) x (-2 a b (t^2)+2 (a-b) x (Sqrt[((t^2)/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]))])-4 (-1+a) b x (Sqrt[((t^2)/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]))])+((a x (1+t x)+x (-2+b+b t x)-2 a b t (t+(x^2)))/(Sqrt[(((x-Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,2]) (x-Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,3]))/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,2] Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,3]))]))) (Sqrt[(1-(x/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,2]))) (1-(x/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,3])))]);

(* then verify equation (3.17) *)
(* this takes a long time, so using fewer terms *)
-(Subscript[P, x,0])Q[x,0]+σ+(Subscript[σ, 0,0])Q[0,0]+(Subscript[σ, 0,1])(Subscript[Q, 0,1])+(Subscript[σ, 1,0])(Subscript[Q, 1,0]);
%/.{Q[x_,y_]->QQeval[11,x,y],(Subscript[Q, 0,1])->QQcxy[11,0,1],(Subscript[Q, 1,0])->QQcxy[11,1,0]}

(* the coefficients in equation (3.18) *)
(Subscript[P, 0]^d)=(-a t+(a^2) t+x-a x+(a^2) (t^2) (x^2)) (-b t+(b^2) t+x-b x+(b^2) (t^2) (x^2));
τ=-(((x^3) (-((8 (-1+a) (-1+b) (Sqrt[((t^2)/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]))]) (2 b x+2 a (x-b (t+x))+a b t x Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]))/(x^2))+(1/(x^3))(Sqrt[((t^2)/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]))]) (8 a b (-b+a (-1+2 b)) (t^2)+8 a b (a+b-2 a b) t x+8 (a (((-1+b))^2)-(-1+b) b+(a^2) (-1+b+2 (b^2) (t^3))) (x^2)+4 a b (a+b-2 a b) (t^2) x Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]-4 a b (a+b-2 a b) t (x^2) Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]+a b (a+b-2 a b) (t^2) (x^2) ((Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1])^2)-8 (-(-1+b) b (x^2)+(a^2) (-(x^2)+2 (b^2) t (t-x+(t^2) (x^2)-t (x^3))+b (t x+(x^2)+(t^2) (-1+(x^3))))+a ((x^2)-2 b (x^2)+(b^2) (t x+(x^2)+(t^2) (-1+(x^3))))) (Sqrt[1-x Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]]))+4 a b (t^2) (-a-b+2 a b-2 a b (t^2) x-((a b (t^2))/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,2]))-((a b (t^2))/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,3])))))/(16 c (Sqrt[((t^2)/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]))]) (Sqrt[1-x Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]])));
(Subscript[τ, 0,0])=-(((x^3) ((1/(x^3))(Sqrt[((t^2)/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]))]) (-8 ((b^2) c-a (b^2) (1+c)+(a^2) (2 (b^2)+c-b (1+c))) (t^2)+8 ((b^2) c-a (b^2) (1+c)+(a^2) (2 (b^2)+c-b (1+c))) t x+8 ((-1+b) (b-2 c)+(a^2) (1-2 (b^2) (t^3)+b (-1+c (t^3)))+a (-1-2 c+2 b (1+c)+(b^2) (-1+c (t^3)))) (x^2)+4 ((b^2) c-a (b^2) (1+c)+(a^2) (2 (b^2)+c-b (1+c))) (t^2) x Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]-4 ((b^2) c-a (b^2) (1+c)+(a^2) (2 (b^2)+c-b (1+c))) t (x^2) Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]+((b^2) c-a (b^2) (1+c)+(a^2) (2 (b^2)+c-b (1+c))) (t^2) (x^2) ((Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1])^2)-8 (2 c (x^2)-b (1+2 c) (x^2)+(b^2) ((x^2)+c t (x+t (-1+(x^3))))+(a^2) ((x^2)-2 (b^2) t (t-x+(t^2) (x^2)-t (x^3))+c t (x+t (-1+(x^3)))-b ((1+c) t x+(x^2)-c (t^3) (x^2)+(1+c) (t^2) (-1+(x^3))))-a (-2 b (1+c) (x^2)+(1+2 c) (x^2)+(b^2) ((1+c) t x+(x^2)-c (t^3) (x^2)+(1+c) (t^2) (-1+(x^3))))) (Sqrt[1-x Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]]))+(1/x)2 (Sqrt[((t^2)/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]))]) (-8 (-1+a) a (-1+b) b+4 (-1+a) (-1+b) (2 (a+b-((a b t)/x))+a b t Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1])-c (((8 ((-1+b) (b t-x) x+(a^2) t (x+(b^2) t (1+2 t (x^2))-b (t+x+(t^2) (x^2)))-a (x (t+x)-b (((t+x))^2)+(b^2) t (t+x+(t^2) (x^2)))))/(x^2))-4 (-1+a) (-1+b) t (-b+a (-1+((b t)/x))) Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]+a b (-1+a+b-a b) (t^2) ((Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1])^2)))-8 c (t^2) (-2 (-1+b) b+a (2-5 b+3 (b^2))+(a^2) (-1+b+2 (b^2) (t^2) (t-x)+2 b (t^2) x)-(((a^2) (-1+b) b (t^2))/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,2]))-(((a^2) (-1+b) b (t^2))/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,3])))+4 (t^2) ((2-3 b) b c+(a^2) (-c+2 (b^2) (-1+c (t^2) (t-x)+(t^2) x)+b (1+c+2 c (t^2) x))+a (2 c-6 b c+(b^2) (1+c (5-2 (t^2) x)))+((a b (-b c+a (b+c-b c)) (t^2))/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,2]))+((a b (-b c+a (b+c-b c)) (t^2))/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,3])))))/(16 c (Sqrt[((t^2)/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]))]) (Sqrt[1-x Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]])));
(Subscript[τ, 0,1])=(((-1+a) b (t^2) (x^2) (-2 a (-1+b)-((a b (t^2) x)/(Sqrt[((t^2)/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]))]))+(a-b) (-1+(Sqrt[1-x Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]]))))/(2 (Sqrt[1-x Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]])));
(Subscript[τ, 1,0])=-((a (-1+b) (t^2) (x^3) (a b (t^2)+((2 (-1+a) b (Sqrt[((t^2)/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]))]))/x)+(((a-b) (Sqrt[((t^2)/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]))]) (-1+(Sqrt[1-x Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]])))/x)))/(2 (Sqrt[((t^2)/(Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]))]) (Sqrt[1-x Root[-4 (t^2)+#1-2 t ((#1)^2)+(t^2) ((#1)^3)&,1]])));

(* verify equation (3.18) *)
-(Subscript[P, 0]^d)(Subscript[Q, 0]^d)[x]+τ+(Subscript[τ, 0,0])Q[0,0]+(Subscript[τ, 0,1])(Subscript[Q, 0,1])+(Subscript[τ, 1,0])(Subscript[Q, 1,0]);
%/.{Q[x_,y_]->QQeval[12,x,y],(Subscript[Q, 0]^d)[x_]->QQdjeval[12,0,x],(Subscript[Q, 0,1])->QQcxy[12,0,1],(Subscript[Q, 1,0])->QQcxy[12,1,0]}

Section 3.3

(* the two power series roots of (Subscript[P, x,0]) *)
x/.Solve[t^2a^2-x(1-t x a)(a-1)==0,x];
(Subscript[ρ, 1])=Select[%,Normal[ApplyToSeries[Factor@*PowerExpand,Series[#,{t,0,2}]]]==t^2a^2/(a-1)&][[1]]
x/.Solve[2t^2a b+x(2-a-b)+t x^2(2a b-a-b)==0,x];
(Subscript[ρ, 2])=Select[%,Normal[ApplyToSeries[Factor@*PowerExpand,Series[#,{t,0,3}]]]==2t^2a b/(a+b-2)&][[1]]
(* the two power series roots of (Subscript[P, 0]^d) *)
x/.Solve[(t a-x)(a-1)+t^2x^2a^2==0,x];
(Subscript[ρ, 3])=Select[%,Normal[ApplyToSeries[Factor@*PowerExpand,Series[#,{t,0,3}]]]==t a&][[1]]
x/.Solve[(t b-x)(b-1)+t^2x^2b^2==0,x];
(Subscript[ρ, 4])=Select[%,Normal[ApplyToSeries[Factor@*PowerExpand,Series[#,{t,0,3}]]]==t b&][[1]]

(* verify some more terms of the expansions *)
ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 1]),{t,0,11}]]
ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 2]),{t,0,8}]]
ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 3]),{t,0,13}]]
ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 4]),{t,0,13}]]

(* here are the algebraic forms of the ζ's *)
{ζ1,(Subscript[ζ1, 0,0]),(Subscript[ζ1, 0,1]),(Subscript[ζ1, 1,0])}=({σ,(Subscript[σ, 0,0]),(Subscript[σ, 0,1]),(Subscript[σ, 1,0])}/.x->(Subscript[ρ, 1]));
{ζ2,(Subscript[ζ2, 0,0]),(Subscript[ζ2, 0,1]),(Subscript[ζ2, 1,0])}=({σ,(Subscript[σ, 0,0]),(Subscript[σ, 0,1]),(Subscript[σ, 1,0])}/.x->(Subscript[ρ, 2]));
{ζ3,(Subscript[ζ3, 0,0]),(Subscript[ζ3, 0,1]),(Subscript[ζ3, 1,0])}=({τ,(Subscript[τ, 0,0]),(Subscript[τ, 0,1]),(Subscript[τ, 1,0])}/.x->(Subscript[ρ, 1]));
{ζ4,(Subscript[ζ4, 0,0]),(Subscript[ζ4, 0,1]),(Subscript[ζ4, 1,0])}=({τ,(Subscript[τ, 0,0]),(Subscript[τ, 0,1]),(Subscript[τ, 1,0])}/.x->(Subscript[ρ, 2]));

(* unfortunately these things have become so complicated that Mathematica struggles to deal with them *)
(* so instead we define these series expansions which it can actually handle *)
ζs1[n_]:=σ/.x->ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 1]),{t,0,n}]]
(Subscript[ζs1, 0,0])[n_]:=(Subscript[σ, 0,0])/.x->ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 1]),{t,0,n}]]
(Subscript[ζs1, 0,1])[n_]:=(Subscript[σ, 0,1])/.x->ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 1]),{t,0,n}]]
(Subscript[ζs1, 1,0])[n_]:=(Subscript[σ, 1,0])/.x->ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 1]),{t,0,n}]]
ζs2[n_]:=σ/.x->ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 2]),{t,0,n}]]
(Subscript[ζs2, 0,0])[n_]:=(Subscript[σ, 0,0])/.x->ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 2]),{t,0,n}]]
(Subscript[ζs2, 0,1])[n_]:=(Subscript[σ, 0,1])/.x->ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 2]),{t,0,n}]]
(Subscript[ζs2, 1,0])[n_]:=(Subscript[σ, 1,0])/.x->ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 2]),{t,0,n}]]
ζs3[n_]:=τ/.x->ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 3]),{t,0,n}]]
(Subscript[ζs3, 0,0])[n_]:=(Subscript[τ, 0,0])/.x->ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 3]),{t,0,n}]]
(Subscript[ζs3, 0,1])[n_]:=(Subscript[τ, 0,1])/.x->ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 3]),{t,0,n}]]
(Subscript[ζs3, 1,0])[n_]:=(Subscript[τ, 1,0])/.x->ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 3]),{t,0,n}]]
ζs4[n_]:=τ/.x->ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 4]),{t,0,n}]]
(Subscript[ζs4, 0,0])[n_]:=(Subscript[τ, 0,0])/.x->ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 4]),{t,0,n}]]
(Subscript[ζs4, 0,1])[n_]:=(Subscript[τ, 0,1])/.x->ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 4]),{t,0,n}]]
(Subscript[ζs4, 1,0])[n_]:=(Subscript[τ, 1,0])/.x->ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 4]),{t,0,n}]]

(* now verify that nothing goes wrong on the rhs when substituting the roots *)
ApplyToSeries[Factor,#]&/@{ζs1[2],(Subscript[ζs1, 0,0])[2],(Subscript[ζs1, 0,1])[2],(Subscript[ζs1, 1,0])[2]}
ApplyToSeries[Factor,#]&/@{ζs2[2],(Subscript[ζs2, 0,0])[2],(Subscript[ζs2, 0,1])[2],(Subscript[ζs2, 1,0])[2]}
ApplyToSeries[Factor,#]&/@{ζs3[2],(Subscript[ζs3, 0,0])[2],(Subscript[ζs3, 0,1])[2],(Subscript[ζs3, 1,0])[2]}
ApplyToSeries[Factor,#]&/@{ζs4[2],(Subscript[ζs4, 0,0])[2],(Subscript[ζs4, 0,1])[2],(Subscript[ζs4, 1,0])[2]}

(* and verify equation (3.25) for each of the four substitutions *)
ζs1[7]+(Subscript[ζs1, 0,0])[7]Q[0,0]+(Subscript[ζs1, 0,1])[7](Subscript[Q, 0,1])+(Subscript[ζs1, 1,0])[7](Subscript[Q, 1,0]);
%/.{Q[x_,y_]->Normal[QQeval[9,x,y]],(Subscript[Q, 0,1])->Normal[QQcxy[9,0,1]],(Subscript[Q, 1,0])->Normal[QQcxy[9,1,0]]}//Simplify
ζs2[7]+(Subscript[ζs2, 0,0])[7]Q[0,0]+(Subscript[ζs2, 0,1])[7](Subscript[Q, 0,1])+(Subscript[ζs2, 1,0])[7](Subscript[Q, 1,0]);
%/.{Q[x_,y_]->Normal[QQeval[9,x,y]],(Subscript[Q, 0,1])->Normal[QQcxy[9,0,1]],(Subscript[Q, 1,0])->Normal[QQcxy[9,1,0]]}//Simplify
ζs3[7]+(Subscript[ζs3, 0,0])[7]Q[0,0]+(Subscript[ζs3, 0,1])[7](Subscript[Q, 0,1])+(Subscript[ζs3, 1,0])[7](Subscript[Q, 1,0]);
%/.{Q[x_,y_]->Normal[QQeval[9,x,y]],(Subscript[Q, 0,1])->Normal[QQcxy[9,0,1]],(Subscript[Q, 1,0])->Normal[QQcxy[9,1,0]]}//Simplify
ζs4[7]+(Subscript[ζs4, 0,0])[7]Q[0,0]+(Subscript[ζs4, 0,1])[7](Subscript[Q, 0,1])+(Subscript[ζs4, 1,0])[7](Subscript[Q, 1,0]);
%/.{Q[x_,y_]->Normal[QQeval[9,x,y]],(Subscript[Q, 0,1])->Normal[QQcxy[9,0,1]],(Subscript[Q, 1,0])->Normal[QQcxy[9,1,0]]}//Simplify

(* checking the determinants *)
(Subscript[DD, 1,2,3])=ApplyToSeries[Factor,(Subscript[ζs1, 0,0])[12](Subscript[ζs2, 0,1])[12](Subscript[ζs3, 1,0])[12]-(Subscript[ζs1, 0,0])[12](Subscript[ζs2, 1,0])[12](Subscript[ζs3, 0,1])[12]-(Subscript[ζs1, 0,1])[12](Subscript[ζs2, 0,0])[12](Subscript[ζs3, 1,0])[12]+(Subscript[ζs1, 0,1])[12](Subscript[ζs2, 1,0])[12](Subscript[ζs3, 0,0])[12]+(Subscript[ζs1, 1,0])[12](Subscript[ζs2, 0,0])[12](Subscript[ζs3, 0,1])[12]-(Subscript[ζs1, 1,0])[12](Subscript[ζs2, 0,1])[12](Subscript[ζs3, 0,0])[12]]
(Subscript[DD, 1,2,4])=ApplyToSeries[Factor,(Subscript[ζs1, 0,0])[2](Subscript[ζs2, 0,1])[2](Subscript[ζs4, 1,0])[2]-(Subscript[ζs1, 0,0])[2](Subscript[ζs2, 1,0])[2](Subscript[ζs4, 0,1])[2]-(Subscript[ζs1, 0,1])[2](Subscript[ζs2, 0,0])[2](Subscript[ζs4, 1,0])[2]+(Subscript[ζs1, 0,1])[2](Subscript[ζs2, 1,0])[2](Subscript[ζs4, 0,0])[2]+(Subscript[ζs1, 1,0])[2](Subscript[ζs2, 0,0])[2](Subscript[ζs4, 0,1])[2]-(Subscript[ζs1, 1,0])[2](Subscript[ζs2, 0,1])[2](Subscript[ζs4, 0,0])[2]]
(Subscript[DD, 1,3,4])=ApplyToSeries[Factor,(Subscript[ζs1, 0,0])[2](Subscript[ζs3, 0,1])[2](Subscript[ζs4, 1,0])[2]-(Subscript[ζs1, 0,0])[2](Subscript[ζs3, 1,0])[2](Subscript[ζs4, 0,1])[2]-(Subscript[ζs1, 0,1])[2](Subscript[ζs3, 0,0])[2](Subscript[ζs4, 1,0])[2]+(Subscript[ζs1, 0,1])[2](Subscript[ζs3, 1,0])[2](Subscript[ζs4, 0,0])[2]+(Subscript[ζs1, 1,0])[2](Subscript[ζs3, 0,0])[2](Subscript[ζs4, 0,1])[2]-(Subscript[ζs1, 1,0])[2](Subscript[ζs3, 0,1])[2](Subscript[ζs4, 0,0])[2]]
(Subscript[DD, 2,3,4])=ApplyToSeries[Factor,(Subscript[ζs2, 0,0])[2](Subscript[ζs3, 0,1])[2](Subscript[ζs4, 1,0])[2]-(Subscript[ζs2, 0,0])[2](Subscript[ζs3, 1,0])[2](Subscript[ζs4, 0,1])[2]-(Subscript[ζs2, 0,1])[2](Subscript[ζs3, 0,0])[2](Subscript[ζs4, 1,0])[2]+(Subscript[ζs2, 0,1])[2](Subscript[ζs3, 1,0])[2](Subscript[ζs4, 0,0])[2]+(Subscript[ζs2, 1,0])[2](Subscript[ζs3, 0,0])[2](Subscript[ζs4, 0,1])[2]-(Subscript[ζs2, 1,0])[2](Subscript[ζs3, 0,1])[2](Subscript[ζs4, 0,0])[2]]

(* solving the systems *)
Q00sol=(Q[0,0]/.Solve[{0==zi+(Subscript[zi, 0,0])Q[0,0]+(Subscript[zi, 0,1])(Subscript[Q, 0,1])+(Subscript[zi, 1,0])(Subscript[Q, 1,0]),0==zj+(Subscript[zj, 0,0])Q[0,0]+(Subscript[zj, 0,1])(Subscript[Q, 0,1])+(Subscript[zj, 1,0])(Subscript[Q, 1,0]),0==zk+(Subscript[zk, 0,0])Q[0,0]+(Subscript[zk, 0,1])(Subscript[Q, 0,1])+(Subscript[zk, 1,0])(Subscript[Q, 1,0])},{Q[0,0],(Subscript[Q, 0,1]),(Subscript[Q, 1,0])}][[1]]);
(* using 1,2,4 *)
Q00sol/.{zi->ζs1[11],(Subscript[zi, 0,0])->(Subscript[ζs1, 0,0])[11],(Subscript[zi, 0,1])->(Subscript[ζs1, 0,1])[11],(Subscript[zi, 1,0])->(Subscript[ζs1, 1,0])[11],zj->ζs2[11],(Subscript[zj, 0,0])->(Subscript[ζs2, 0,0])[11],(Subscript[zj, 0,1])->(Subscript[ζs2, 0,1])[11],(Subscript[zj, 1,0])->(Subscript[ζs2, 1,0])[11],zk->ζs4[11],(Subscript[zk, 0,0])->(Subscript[ζs4, 0,0])[11],(Subscript[zk, 0,1])->(Subscript[ζs4, 0,1])[11],(Subscript[zk, 1,0])->(Subscript[ζs4, 1,0])[11]};
ApplyToSeries[Expand@*Factor,%]
(* using 1,3,4 *)
Q00sol/.{zi->ζs1[11],(Subscript[zi, 0,0])->(Subscript[ζs1, 0,0])[11],(Subscript[zi, 0,1])->(Subscript[ζs1, 0,1])[11],(Subscript[zi, 1,0])->(Subscript[ζs1, 1,0])[11],zj->ζs3[11],(Subscript[zj, 0,0])->(Subscript[ζs3, 0,0])[11],(Subscript[zj, 0,1])->(Subscript[ζs3, 0,1])[11],(Subscript[zj, 1,0])->(Subscript[ζs3, 1,0])[11],zk->ζs4[11],(Subscript[zk, 0,0])->(Subscript[ζs4, 0,0])[11],(Subscript[zk, 0,1])->(Subscript[ζs4, 0,1])[11],(Subscript[zk, 1,0])->(Subscript[ζs4, 1,0])[11]};
ApplyToSeries[Expand@*Factor,%]
(* using 2,3,4 *)
Q00sol/.{zi->ζs2[11],(Subscript[zi, 0,0])->(Subscript[ζs2, 0,0])[11],(Subscript[zi, 0,1])->(Subscript[ζs2, 0,1])[11],(Subscript[zi, 1,0])->(Subscript[ζs2, 1,0])[11],zj->ζs3[11],(Subscript[zj, 0,0])->(Subscript[ζs3, 0,0])[11],(Subscript[zj, 0,1])->(Subscript[ζs3, 0,1])[11],(Subscript[zj, 1,0])->(Subscript[ζs3, 1,0])[11],zk->ζs4[11],(Subscript[zk, 0,0])->(Subscript[ζs4, 0,0])[11],(Subscript[zk, 0,1])->(Subscript[ζs4, 0,1])[11],(Subscript[zk, 1,0])->(Subscript[ζs4, 1,0])[11]};
ApplyToSeries[Expand@*Factor,%]

(* and indeed, *)
QQcxy[9,0,0]

(* checking the b=1 case *)
(* give the coefficients a new name *)
{ξ1,(Subscript[ξ1, 0,0]),(Subscript[ξ1, 0,1])}=({σ,(Subscript[σ, 0,0]),(Subscript[σ, 0,1])}/.x->(Subscript[ρ, 1])/.b->1);
{ξ2,(Subscript[ξ2, 0,0]),(Subscript[ξ2, 0,1])}=({σ,(Subscript[σ, 0,0]),(Subscript[σ, 0,1])}/.x->(Subscript[ρ, 2])/.b->1);

(* again, using series is easier *)
ξs1[n_]:=σ/.b->1/.x->ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 1])/.b->1,{t,0,n}]]
(Subscript[ξs1, 0,0])[n_]:=(Subscript[σ, 0,0])/.b->1/.x->ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 1])/.b->1,{t,0,n}]]
(Subscript[ξs1, 0,1])[n_]:=(Subscript[σ, 0,1])/.b->1/.x->ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 1])/.b->1,{t,0,n}]]
ξs2[n_]:=σ/.b->1/.x->ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 2])/.b->1,{t,0,n}]]
(Subscript[ξs2, 0,0])[n_]:=(Subscript[σ, 0,0])/.b->1/.x->ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 2])/.b->1,{t,0,n}]]
(Subscript[ξs2, 0,1])[n_]:=(Subscript[σ, 0,1])/.b->1/.x->ApplyToSeries[PowerExpand@*Factor@*PowerExpand,Series[(Subscript[ρ, 2])/.b->1,{t,0,n}]]

(* verifying equation (3.30) *)
(Subscript[DD, 1,2])=ApplyToSeries[Factor@*Simplify,(Subscript[ξs1, 0,0])[5](Subscript[ξs2, 0,1])[5]-(Subscript[ξs1, 0,1])[5](Subscript[ξs2, 0,0])[5]]

(* solving the system *)
Q00solb1=(Q[0,0]/.Solve[{0==xi+(Subscript[xi, 0,0])Q[0,0]+(Subscript[xi, 0,1])(Subscript[Q, 0,1]),0==xj+(Subscript[xj, 0,0])Q[0,0]+(Subscript[xj, 0,1])(Subscript[Q, 0,1])},{Q[0,0],(Subscript[Q, 0,1])}][[1]]);
Q00solb1/.{xi->ξs1[14],(Subscript[xi, 0,0])->(Subscript[ξs1, 0,0])[14],(Subscript[xi, 0,1])->(Subscript[ξs1, 0,1])[14],xj->ξs2[14],(Subscript[xj, 0,0])->(Subscript[ξs2, 0,0])[14],(Subscript[xj, 0,1])->(Subscript[ξs2, 0,1])[14]};
ApplyToSeries[Expand@*Factor,%]

(* and indeed *)
QQcxy[9,0,0]/.b->1