q-Heisenberg algebra generated by a,w including an R matrix satisfying QYBE. We work with all (Subscript[b, i ])set to b.
Note how ϵ was rescaled.

ϵ /: (ϵ^(n_))/;n>1 :=0;
PBWBasis={a,w}; (*a = u/(1-t)*)
B[U@a,U@w]=-(B[U@w,U@a]=U[]+U[a,w]ϵ );

UU[l___,((x_)^(n_)),r___]:= UU[l,Sequence@@Table[x,{n}],r];
UU[l___,1,r___]:= UU[l,r];
UU[]=U[]; UU[l_,r___]:=U[l]**UU[r];
(Subscript[U, i_])[ℰ_] := ℰ /. {t->t, u_U :> UU@@Replace[u, x_ :> (Subscript[x, i]),1]};

B[x_,x_]=0;
B[U[(Subscript[(x_), i_])],U[(Subscript[(y_), i_])]]:=(Subscript[U, i])[B[U@x,U@y]];
B[U[(Subscript[(x_), i_])],U[(Subscript[(y_), j_])]] /; i=!=j := 0;
B[x_,y_] := x**y-y**x;

x_⪯y_:=OrderedQ[{x,y}]; x_≺y_:=!OrderedQ[{y,x}];
Simp[ℰ_]:= Collect[ℰ, _U, Together];

Unprotect[NonCommutativeMultiply];
NonCommutativeMultiply[x_] := x;
0**_=_**0=0;
x_**U[]:=x; U[]**x_:=x;
(a_*x_U)**(b_*y_U) := If[a b === 0, 0, Simp[a b (x**y)]];
(a_*x_U)**y_ := Simp[a(x**y)]; x_**(a_*y_U) := Simp[a(x**y)];
(x_Plus)**y_ := (#**y)& /@ x; x_**(y_Plus) := (x**#)& /@ y;
U[x_]**U[y_] := (*U[x]**U[y] =*) If[x≺y, U[x,y], U[y,x]+B[U@x,U@y]];
U[x_] ** U[y1_,yy__] := (*U[x] ** U[y1,yy] =*)If[x⪯y1,U[x,y1,yy],(U@x**U@y1)**U@yy];
U[xx__,xn_]**U[yy__] := (*U[xx,xn]**U[yy] =*) U@xx**(U@xn**U@yy);

ToDegree[n_][ℰ_] := (Simp[ℰ] /. {ϵ->ℏ ϵ, (Subscript[b, i_]):>ℏ (Subscript[b, i]), (Subscript[t, i_]):>(E^(ℏ (Subscript[b, i]))),b->ℏ b,  t:>(E^(ℏ b)),x_U:>(ℏ^(Count[x,a|(Subscript[a, _])]))x} /. a_.x_U :>Normal[Series[a,{ℏ,0,n}]]*x)/.ℏ->1

(*First check that the element CC behaves like c did in 1-co.*)

CCC=UU[a,w]-(ϵ/2) UU[(a^2),(w^2)];
B[U[w],CCC]//Simp
B[U[a],CCC]//Simp

Testing Yang-Baxter

(*We need to deal with powers (((aw))^m), would be good to automate this. *)
MixedPow[r_,s_,m_]:=UU@@Flatten[Table[{r,s},{x,1,m}]]

(*Powers of (Subscript[c, j]) are now this, compare with CCC*)
CC[j_,m_] := If[m==0,U[],MixedPow[(Subscript[a, j]),(Subscript[w, j]),m-1]**(UU[(Subscript[a, j]),(Subscript[w, j])]-m(ϵ/2) UU[((Subscript[a, j])^2),((Subscript[w, j])^2)])]//Simp

(*Set c = U[a,w]-(ϵ/2) UU[(a^2),(w^2)] *)

CC[2,1]

(*The rescaling of ϵ turns up in the R-matrix only. This is where the strange rational functions in t come from. Previously they appeared in the Logos.
Note also that it is not yet in canonical order since CC contains mixed powers awawawawaw.*)

RAW[i_,j_,d_]:= Sum[(((((1-t))^n)UU[((Subscript[a, i])^n)]**((b^m)U[]+ϵ ((t-1)/(3t-1)) m (b^(m-1))U[(Subscript[a, i]),(Subscript[w, i])])**CC[j,m]**UU[((Subscript[w, j])^n)])/(m!n!))(1-(ϵ/4)((t-1)/(3t-1))(n-1)n),{m,0,d},{n,0,d-m}]//ToDegree[d]

RAW[1,2,1]
RAW[1,2,2]
RAW[1,2,3]
RAW[1,2,4]

(*Now let's check the QYBE*)
R3AW[d_] := (ToDegree[d][RAW[1,2,d]**RAW[1,3,d]]**RAW[2,3,d])-(ToDegree[d][RAW[2,3,d]**RAW[1,3,d]]**RAW[1,2,d])

R3AW[3]//ToDegree[3]

R3AW[4]//ToDegree[4]

R3AW[5]//ToDegree[5]

R3AW[6]//ToDegree[6]

(*Is this sufficient evidence?*)

R3AW[7]//ToDegree[7]

Timing[R3AW[8]//ToDegree[8]]

Timing[R3AW[9]//ToDegree[9]]

Timing[R3AW[10]//ToDegree[10]]

Logos

q = 1+ϵ;
qI[k_]:=k(1+ϵ(k-1)/2)//Expand
qFac[n_]:=n!(1+ϵ (n-1)n/4)//Expand
InvqFac[n_]:=(1-ϵ (n-1)n/4)/n!//Expand
qBin[n_,k_]:= Binomial[n,k](1+ϵ k (n-k) /2)//Expand

(*Checking the commutation relation for powers of u,w*)
WmUn[m_,n_]:=Sum[(q^((m-j)(n-j)))qBin[m,j]qBin[n,j]qFac[j]UU[(a^(n-j)),(w^(m-j))],{j,0,Min[m,n]}]
TestWmUn[m_,n_]:=-UU[(w^m),(a^n)]+WmUn[m,n]
TestWmUn[6,5]//Simp

(*Guess qLogos first at q=1*)
ToDegh[F_,x_]:=Series[F/.{α->α h,δ->δ h,β->β h},{h,0,x}]
d=8;
LHS=Sum[(α^m)(δ^k)(β^n)InvqFac[m]InvqFac[k]InvqFac[n]WmUn[m+k,n+k],{m,0,d},{k,0,d-m},{n,0,d-m-k}]//Simp;
ν=Sum[(δ^j),{j,0,d}];
RHS = Sum[(ν^(m+k+n+j+1))(α^(m+j))(δ^k)(β^(n+j))InvqFac[m]InvqFac[k]InvqFac[n]InvqFac[j]UU[(a^(n+k)),(w^(m+k))],{m,0,d},{k,0,d-m},{n,0,d-m-k},{j,0,d-m-n-k}]//Simp; 
ToDegh[(LHS-RHS)/.ϵ->0//Simp,6]

(*powers of nu as power series*)
nuA[z_]:=Sum[Binomial[z-1+x,x](δ^x),{x,0,d}]

(*Now set up the LHS and RHS to find the Logos relating them.*)
d=14;
LHS=Sum[(α^m)(δ^k)(β^n)InvqFac[m]InvqFac[k]InvqFac[n]WmUn[m+k,n+k],{m,0,d},{k,0,d-m},{n,0,d-m-k}]//Simp;

cRHS = Sum[nuA[m+k+n+j+1](α^(m+j))(δ^k)(β^(n+j))InvqFac[m]InvqFac[k]InvqFac[n]InvqFac[j]UU[(a^(n+k)),(w^(m+k))],{m,0,d},{k,0,d-m},{n,0,d-m-k},{j,0,d-m-n-k}]//Simp;

μ = 1-δ;

(*Here we guess and verify the Logos coefficient by coefficient*)
ToDegh[Coefficient[(-LHS+cRHS+
nuA[4]ϵ (
α β δ μ(μ+(1/2)α β)cRHS +
 β δ μ(μ+α β ) U[a]**cRHS+
 α δ μ(μ+α β ) cRHS**U[w]+
 (δ (μ^2)+α β (1-(δ^2)) ) U[a]**cRHS**U[w]+
(1/2) (β^2)δ μ   U[a,a]**cRHS+
(1/2) (α^2)δ μ   cRHS**U[w,w]+
β  δ μ U[a,a]**cRHS**U[w]+
α  δ μ U[a]**cRHS**U[w,w]+
(1/2)(δ^2)μ(μ+1) U[a,a]**cRHS**U[w,w]
)
),U[]],14]

(*Here we check the Logos is correct up to degree 14*)
ToDegh[Simp[(-LHS+cRHS+
nuA[4]ϵ (
α β δ μ(μ+(1/2)α β)cRHS +
 β δ μ(μ+α β ) U[a]**cRHS+
 α δ μ(μ+α β ) cRHS**U[w]+
 (δ (μ^2)+α β (1-(δ^2)) ) U[a]**cRHS**U[w]+
(1/2) (β^2)δ μ   U[a,a]**cRHS+
(1/2) (α^2)δ μ   cRHS**U[w,w]+
β  δ μ U[a,a]**cRHS**U[w]+
α  δ μ U[a]**cRHS**U[w,w]+
(1/2)(δ^2)μ(μ+1) U[a,a]**cRHS**U[w,w]
)
)],14]

Logos without q-factorials

(*Repeat with usual factorials in the denominator.*)
d=14;
LHS=Sum[(α^m)(δ^k)((β^n)/(m!k!n!))WmUn[m+k,n+k],{m,0,d},{k,0,d-m},{n,0,d-m-k}]//Simp;

cRHS = Sum[nuA[m+k+n+j+1](α^(m+j))(δ^k)((β^(n+j))/(m!k!n!j!))UU[(a^(n+k)),(w^(m+k))],{m,0,d},{k,0,d-m},{n,0,d-m-k},{j,0,d-m-n-k}]//Simp;

(*Here we guess and verify the Logos  coefficient by coefficient*)
ToDegh[Coefficient[(-LHS+cRHS+
nuA[4]ϵ (
((1/2)(δ^2)(μ^2)+α β δ μ+(1/4)(α^2)(β^2))cRHS +
β  (δ μ+(1/2)α β) U[a]**cRHS+
α  (δ μ+(1/2)α β) cRHS**U[w]+
 (δ +α β-(δ^2)) U[a]**cRHS**U[w]+
(1/2) (β^2)δ (μ+(1/2)δ)   U[a,a]**cRHS+
(1/2) (α^2)δ (μ+(1/2)δ)   cRHS**U[w,w]+
β  δ( μ+(1/2)δ) U[a,a]**cRHS**U[w]+
α  δ( μ+(1/2)δ) U[a]**cRHS**U[w,w]+
(δ^2)(μ+(1/4)(δ^2)) U[a,a]**cRHS**U[w,w]
)
),U[a,a,w,w]],10]