Pensieve header: The “Speedy” engine, with t=b-ϵa and with w/ω.

The “Speedy” Engine

Internal Utilities

Canonical Form:

CCF[ℰ_] := (Subscript[PP, CCF])@ExpandDenominator@ExpandNumerator@(Subscript[PP, Together])@Together[(Subscript[PP, Exp])[
Expand[ℰ] //. (E^(x_))(E^(y_)):>(E^(x+y)) /. (E^(x_)):>(E^(CCF[x]))]];
CF[ℰ_List]:=CF/@ℰ;
CF[sd_SeriesData] := MapAt[CF,sd,3];
CF[ℰ_] := (Subscript[PP, CF])@Module[
{vs=Cases[ℰ,(Subscript[(y|b|t|w|a|x|η|β|τ|ω|α|ξ), _]),∞]⋃{y,b,t,w,a,x,η,β,τ,ω,α,ξ}},
Total[CoefficientRules[Expand[ℰ], vs] /. (ps_->c_) :> CCF[c](Times@@(vs^ps))]
];
CF[ℰ_𝔼]:=CF/@ℰ; CF[(Subscript[𝔼, sp___])[ℰs___]]:=CF/@(Subscript[𝔼, sp])[ℰs];

The Kronecker δ:

Kδ/:(Subscript[Kδ, i_,j_]):=If[i===j,1,0];

Equality, multiplication, and degree-adjustment of perturbed Gaussians; 𝔼[L,Q,P] stands for (e^(L+Q))P:

𝔼/:𝔼[L1_,Q1_,P1_]≡𝔼[L2_,Q2_,P2_]:=CF[L1==L2]∧CF[Q1==Q2]∧CF[Normal[P1-P2]==0];
𝔼/:𝔼[L1_,Q1_,P1_]𝔼[L2_,Q2_,P2_]:=𝔼[L1+L2,Q1+Q2,P1*P2];
(Subscript[𝔼[L_,Q_,P_], $k_]):=𝔼[L,Q,Series[Normal@P,{ϵ,0,$k}]];

Zip and Bind

Variables and their duals:

{(SuperStar[t]),(SuperStar[b]),(SuperStar[y]),(SuperStar[w]),(SuperStar[a]),(SuperStar[x]),(SuperStar[z])}={τ,β,η,ω,α,ξ,ζ}; 
{(SuperStar[τ]),(SuperStar[β]),(SuperStar[η]),(SuperStar[ω]),(SuperStar[α]),(SuperStar[ξ]),(SuperStar[ζ])}={t,b,y,w,a,x,z}; (SuperStar[((Subscript[u_, i_]))]):=(Subscript[((SuperStar[u])), i]);

Upper to lower and lower to Upper:

U2l={(Subscript[B, i_]^(p_.)):>(E^(-p ℏ γ (Subscript[b, i]))),(B^(p_.)):>(E^(-p ℏ γ b)), (Subscript[T, i_]^(p_.)):>(E^(-p ℏ (Subscript[t, i]))),(T^(p_.)):>(E^(-p ℏ t)),(Subscript[𝒜, i_]^(p_.)):>(E^(p γ (Subscript[α, i]))),(𝒜^(p_.)):>(E^(p γ α))}; 
l2U={(E^(c_.(Subscript[b, i_])+d_.)):>(Subscript[B, i]^(-c/(ℏ γ)))(E^d),(E^(c_.b+d_.)):>(B^(-c/(ℏ γ)))(E^d),
(E^(c_.(Subscript[t, i_])+d_.)):>(Subscript[T, i]^(-c/ℏ))(E^d),(E^(c_.t+d_.)):>(T^(-c/ℏ))(E^d),
(E^(c_.(Subscript[α, i_])+d_.)):>(Subscript[𝒜, i]^(c/γ))(E^d),(E^(c_.α+d_.)):>(𝒜^(c/γ))(E^d),
(E^(ℰ_)):>(E^(Expand@ℰ))};

Derivatives in the presence of exponentiated variables:

(Subscript[D, b])[f_]:=(∂_b)f-ℏ γ B(∂_B)f; (Subscript[D, (Subscript[b, i_])])[f_]:=(∂_(Subscript[b, i]))f-ℏ γ (Subscript[B, i])(∂_(Subscript[B, i]))f;
(Subscript[D, t])[f_]:=(∂_t)f-ℏ T(∂_T)f; (Subscript[D, (Subscript[t, i_])])[f_]:=(∂_(Subscript[t, i]))f-ℏ (Subscript[T, i])(∂_(Subscript[T, i]))f; 
(Subscript[D, α])[f_]:=(∂_α)f+γ 𝒜 (∂_𝒜)f; (Subscript[D, (Subscript[α, i_])])[f_]:=(∂_(Subscript[α, i]))f+γ (Subscript[𝒜, i]) (∂_(Subscript[𝒜, i]))f;
(Subscript[D, v_])[f_]:=(∂_v)f; (Subscript[D, {v_,0}])[f_]:=f; (Subscript[D, {}])[f_]:=f; (Subscript[D, {v_,n_Integer}])[f_]:=(Subscript[D, v])[(Subscript[D, {v,n-1}])[f]];
(Subscript[D, {l_List,ls___}])[f_]:=(Subscript[D, {ls}])[(Subscript[D, l])[f]];

Finite Zips:

collect[sd_SeriesData, ζ_] := MapAt[collect[#,ζ]&,sd,3];
collect[ℰ_,ζ_] := (Subscript[PP, Collect])@Collect[ℰ,ζ];
(Subscript[Zip, {}])[P_]:=P;
(Subscript[Zip, ζs_])[Ps_List]:=(Subscript[Zip, ζs])/@Ps;
(Subscript[Zip, {ζ_,ζs___}])[P_] := (Subscript[PP, Zip])[
(collect[P//(Subscript[Zip, {ζs}]),ζ] /. f_.(ζ^(d_.)):>((Subscript[D, {(SuperStar[ζ]),d}])[f]))/.(SuperStar[ζ])->0/.(((SuperStar[ζ])/.{b->B,t->T,α->𝒜})->1)]

QZip implements the “Q-level zips” on 𝔼(L,Q,P)=(e^(L+Q))P(ϵ). Such zips regard the L variables as scalars.

(Subscript[QZip, ζs_List])@𝔼[L_,Q_,P_] := (Subscript[PP, QZip])@Module[{ζ,z,zs,c,ys,ηs,qt,zrule,ζrule,out},
zs=Table[(SuperStar[ζ]),{ζ,ζs}];
c=CF[Q/.Alternatives@@(ζs⋃zs)->0];
ys=CF@Table[(∂_ζ)(Q/.Alternatives@@zs->0),{ζ,ζs}];
ηs=CF@Table[(∂_z)(Q/.Alternatives@@ζs->0),{z,zs}];
qt=CF@Inverse@Table[(Subscript[Kδ, z,(SuperStar[ζ])])-(∂_(z,ζ))Q,{ζ,ζs},{z,zs}];
zrule=Thread[zs->CF[qt.(zs+ys)]];
ζrule=Thread[ζs->ζs+ηs.qt];
CF /@ 𝔼[L,c+ηs.qt.ys,Det[qt](Subscript[Zip, ζs])[P/.(zrule⋃ζrule)]] ];

LZip implements the “L-level zips” on 𝔼(L,Q,P)=((Pe)^(L+Q)). Such zips regard all of ((Pe)^Q) as a single”P”. Here the z’s are b and α and the ζ’s are β and a.

(Subscript[LZip, ζs_List])@𝔼[L_,Q_,P_] := (Subscript[PP, LZip])@Module[{ζ,z,zs,Zs,c,ys,ηs,lt,zrule,Zrule,ζrule,Q1,EEQ,EQ},
zs=Table[(SuperStar[ζ]),{ζ,ζs}];
Zs=zs/.{b->B,t->T,α->𝒜};
c=L/.Alternatives@@(ζs⋃zs)->0/.Alternatives@@Zs->1;
ys=Table[(∂_ζ)(L/.Alternatives@@zs->0),{ζ,ζs}];
ηs=Table[(∂_z)(L/.Alternatives@@ζs->0),{z,zs}];
lt=Inverse@Table[(Subscript[Kδ, z,(SuperStar[ζ])])-(∂_(z,ζ))L,{ζ,ζs},{z,zs}];
zrule=Thread[zs->lt.(zs+ys)];
Zrule=Join[zrule,zrule/.r_Rule:>((U=r[[1]]/.{b->B,t->T,α->𝒜})->(U/.U2l/.r//.l2U))];
ζrule=Thread[ζs->ζs+ηs.lt];
Q1=Q/.(Zrule⋃ζrule);
EEQ[ps___]:=EEQ[ps]=(Subscript[PP, "EEQ"])@(CF[(E^(-Q1))(Subscript[D, Thread[{zs,{ps}}]])[(E^Q1)]]/.{Alternatives@@zs->0,Alternatives@@Zs->1});
CF@𝔼[c+ηs.lt.ys,Q1/.{Alternatives@@zs->0,Alternatives@@Zs->1},
Det[lt]((Subscript[Zip, ζs])[(EQ@@zs)(P/.(Zrule⋃ζrule))] /. Derivative[ps___][EQ][___] :> EEQ[ps] /. _EQ->1)  ] ];

(Subscript[B, {}])[L_,R_]:=L R;
(Subscript[B, {is__}])[L_𝔼,R_𝔼] := (Subscript[PP, B])@Module[{n},
Times[
L/.Table[(Subscript[(v:b|B|t|T|w|a|x|y), i])->(Subscript[v, n@i]),{i,{is}}],
R/.Table[(Subscript[(v:β|τ|ω|α|𝒜|ξ|η), i])->(Subscript[v, n@i]),{i,{is}}]
] // (Subscript[LZip, Join@@Table[{(Subscript[β, n@i]),(Subscript[τ, n@i]),(Subscript[w, n@i]),(Subscript[a, n@i])},{i,{is}}]]) // (Subscript[QZip, Join@@Table[{(Subscript[ξ, n@i]),(Subscript[y, n@i])},{i,{is}}]]) ];
(Subscript[B, is___])[L_,R_]:=(Subscript[B, {is}])[L,R];

𝔼 morphisms with domain and range.

(Subscript[B, is_List])[(Subscript[𝔼, d1_->r1_])[L1_,Q1_,P1_],(Subscript[𝔼, d2_->r2_])[L2_,Q2_,P2_]]:=(Subscript[𝔼, (d1⋃Complement[d2,is])->(r2⋃Complement[r1,is])])@@(Subscript[B, is])[𝔼[L1,Q1,P1],𝔼[L2,Q2,P2]];
(Subscript[𝔼, d1_->r1_])[L1_,Q1_,P1_]//(Subscript[𝔼, d2_->r2_])[L2_,Q2_,P2_] := (Subscript[B, r1⋂d2])[(Subscript[𝔼, d1->r1])[L1,Q1,P1],(Subscript[𝔼, d2->r2])[L2,Q2,P2]];
(Subscript[𝔼, d1_->r1_])[L1_,Q1_,P1_]≡(Subscript[𝔼, d2_->r2_])[L2_,Q2_,P2_] ^:= (d1==d2)∧(r1==r2)∧(𝔼[L1,Q1,P1]≡𝔼[L2,Q2,P2]);
(Subscript[𝔼, d1_->r1_])[L1_,Q1_,P1_](Subscript[𝔼, d2_->r2_])[L2_,Q2_,P2_] ^:= (Subscript[𝔼, (d1⋃d2)->(r1⋃r2)])@@(𝔼[L1,Q1,P1]𝔼[L2,Q2,P2]);
(Subscript[(Subscript[𝔼, dr_])[L_,Q_,P_], $k_]):=(Subscript[𝔼, dr])@@(Subscript[𝔼[L,Q,P], $k]);
(Subscript[𝔼, _])[ℰ___][i_]:={ℰ}[[i]];

𝔼[Λ]

(Subscript[𝔼, dr_])[Λ_]:=CF@Module[{L,Λ0=Limit[Λ,ϵ->0]},(Subscript[(Subscript[𝔼, dr])[L=Λ0 /. (Subscript[(η|y|ξ|x), _])->0,Λ0-L,(E^(Λ-Λ0))], $k])/.l2U]

Exponentials as needed.

Task. Define  (Subscript[Exp, m,i,k])[P] to compute (e^(𝕆(P))) to (ϵ^k) in the using the (Subscript[m, i,i->i]) multiplication, where P is an ϵ-dependent near-docile element, giving the answer in 𝔼-form.
Methodology. If (Subscript[P, 0]):=(Subscript[P, ϵ=0]) and (e^(λ 𝕆(P)))=𝕆((e^(λ (Subscript[P, 0])))F(λ)), then F(λ=0)=1 and we have:
𝕆((e^(λ (Subscript[P, 0])))((Subscript[P, 0])F(λ)+(∂_λ)F))=𝕆((∂_λ)(e^(λ (Subscript[P, 0])))F(λ))=(∂_λ)𝕆((e^(λ (Subscript[P, 0])))F(λ))=(∂_λ)(e^(λ 𝕆(P)))=(e^(λ 𝕆(P)))𝕆(P)=𝕆((e^(λ (Subscript[P, 0])))F(λ))𝕆(P).
This is a linear ODE for F. Setting inductively (Subscript[F, k])=(Subscript[F, k-1])+(ϵ^k)φ we find that (Subscript[F, 0])=1 and solve for φ.

(* Bug: The first line is valid only if 𝕆((E^(Subscript[P, 0])))==(E^(𝕆((Subscript[P, 0])))). *)
(Subscript[Exp, m_,i_,0])[P_]:=Module[{LQ=Normal@P /. ϵ->0},
𝔼[LQ/.(Subscript[(x|y), i])->0, LQ/.(Subscript[(b|a|t), i])->0,1] ];

(Subscript[Exp, m_,i_,k_])[P_]:=Block[{$k=k},
Module[{P0,λ,φ,φs,F,j,rhs,eqn,pows,at0,atλ},
P0=Normal@P /. ϵ->0;
F=Normal@Last@(Subscript[Exp, m,i,k-1])[λ P];
While[
rhs=(Subscript[m, i,j->i])[(Subscript[(Subscript[𝔼, {}->{i}])[λ P0/.(Subscript[(x|y), i])->0,λ P0/.(Subscript[(b|a|t), i])->0, F], k])(Subscript[sσ, i->j])@(Subscript[(Subscript[𝔼, {}->{i}])[0,0,P], k])]//Last//Normal;
eqn=CF[((∂_λ)F)+P0 F-rhs];
eqn=!=0,(*do*)
pows=First/@CoefficientRules[eqn,{(Subscript[y, i]),(Subscript[b, i]),(Subscript[a, i]),(Subscript[x, i])}];
F+=Sum[(ϵ^k)(Subscript[φ, js])[λ]Times@@(({(Subscript[y, i]),(Subscript[b, i]),(Subscript[a, i]),(Subscript[x, i])})^js),{js,pows}];
rhs=(Subscript[m, i,j->i])[(Subscript[(Subscript[𝔼, {}->{i}])[λ P0/.(Subscript[(x|y), i])->0,λ P0/.(Subscript[(b|a|t), i])->0, F], k])(Subscript[sσ, i->j])@(Subscript[(Subscript[𝔼, {}->{i}])[0,0,P], k])]//Last//Normal;
eqn=CF[((∂_λ)F)+P0 F-rhs];
φs=Table[(Subscript[φ, js])[λ],{js,pows}];
at0=Table[(Subscript[φ, js])[0]==0,{js,pows}];
atλ=(#==0)&/@(pows/.CoefficientRules[eqn,{(Subscript[y, i]),(Subscript[b, i]),(Subscript[a, i]),(Subscript[x, i])}]);
F=F/. DSolve[And@@(at0⋃atλ),φs,λ][[1]]
];
(Subscript[𝔼, {}->{i}])[P0/.(Subscript[(x|y), i])->0,P0/.(Subscript[(b|a|t), i])->0, F+((O[ϵ])^(k+1))/.λ->1]] ]

“Define” Code

Define[lhs=rhs, ...] defines the lhs to be rhs, except that rhs is computed only once for each value of $k. Fancy Mathematica not for the faint of heart. Most readers should ignore.

SetAttributes[Define, HoldAll];
Define[def_,defs__] := (Define[def]; Define[defs];);
Define[(Subscript[op_, is__])=ℰ_] := Module[{SD,ii,jj,kk,isp,nis,nisp,sis},Block[{i,j,k},
ReleaseHold[Hold[
SD[(Subscript[op, nisp,$k_Integer]),(Subscript[PP, Boot])@Block[{i,j,k},(Subscript[op, isp,$k])=ℰ;(Subscript[op, nis,$k])]];
SD[(Subscript[op, isp]),(Subscript[op, {is},$k])]; SD[(Subscript[op, sis__]),(Subscript[op, {sis}])];
]/. {SD->SetDelayed,
isp->{is}/.{i->i_,j->j_,k->k_},
nis->{is}/.{i->ii,j->jj,k->kk},
nisp->{is}/.{i->ii_,j->jj_,k->kk_}
}]  ]]