Define quantities and useful functions

Use the “v” dimension in VEST as the k dimension in the calculation

<<Notation`
<<VEST`
DefObject[u0,1,{True,True},"(Subscript[u, 0])"]
DefObject[b0,1,{True,True},"(Subscript[b, 0])"]
DefObject[u10,1,{True,True},"(u^((0)))"](* This is just a placeholder*)
DefObject[b10,1,{True,True},"(b^((0)))"](* This is just a placeholder*)
DefObject[Bc,1,{False,False},"B"]
DefObject[Bt,2,{False,False},"B"]
DefObject[Ut,2,{False,False},"U"]
DefObject[Ω,1,{False,False},"Ω"]
DefObject[km,0,{False,True},"k"](* magnitude of k (or v) *)
DefObject[gh,1,{False,False},"(OverHat[g])"](* Gradient direction for density and velocity stratification*)
DefObject[gρ,0,{False,False},"∇ln(ρ)"](* ∇ln(ρ), same direction as ∇ln((Subscript[u, rms])) *)
DefObject[gt,0,{False,False},"∇ln((Subscript[u, rms]))"](* ∇ln((Subscript[u, rms])), no need for a separate (Subscript[b, rms]) direction since contributions all separate completely*)
DefObject[ρ0m1,0,{False,False},"((Subscript[ρ, 0])^(-1))"] (* ((Subscript[ρ, 0])^(-1)), for use in the induction equation, just a constant *)
DefObject[K,1,{False,False},"K"](* Large scale K *)
DefObject[v0stat,2,{True,True},"((OverTilde[v])^((0)))"]
DefObject[Wu,0,{True,True},"(Subscript[(OverTilde[W]), u])"](* Isotropic nonhelical u fluctuations *)
DefObject[Hu,0,{True,True},"(Subscript[(OverTilde[H]), u])"](* Isotropic helical u fluctuations *)
DefObject[b0stat,2,{True,True},"((OverTilde[b])^((0)))"]
DefObject[Wb,0,{True,True},"(Subscript[(OverTilde[W]), b])"](* Isotropic nonhelical b fluctuations *)
DefObject[Hb,0,{True,True},"(Subscript[(OverTilde[H]), b])"](* Isotropic helical b fluctuations *)
DefObject[h0stat,2,{True,True},"((OverTilde[h])^((0)))"]
DefObject[Hc,0,{True,True},"(Subscript[(OverTilde[H]), C])"](* Isotropic ub correlation *)
DefConst[{ν,η}];(* Viscosity and resitivity *)
DefIndices[{i,j,k,l,m,n,p,q,r,s,t}];

(*List of all E and N and their derivatives and conjugates *)ENlist ={ Eη,Nν,dEη,dNν,ddEη,ddNν,dddEη,dddNν,Eηs,Nνs,dEηs,dNνs,ddEηs,ddNνs,dddEηs,dddNνs};
ENprint={"(Subscript[E, η])","(Subscript[N, ν])","(Subscript[E, η])'","(Subscript[N, ν])'","(Subscript[E, η])''","(Subscript[N, ν])''","(Subscript[E, η])'''","(Subscript[N, ν])'''",
"(SuperStar[(Subscript[E, η])])","(SuperStar[(Subscript[N, ν])])","(SuperStar[(Subscript[E, η])])'","(SuperStar[(Subscript[N, ν])])'","(SuperStar[(Subscript[E, η])])''","(SuperStar[(Subscript[N, ν])])''","(SuperStar[(Subscript[E, η])])'''","(SuperStar[(Subscript[N, ν])])'''"};
ENMember[_]=False;
MapThread[(ENMember[#1[]]=True;DefObject[#1,0,{False,True},#2])&,{ENlist,ENprint}];

(* Some ghetto functions to take derivatives *)
ClearAll[dEN]
SetAttributes[dEN,HoldFirst]
dEN[var_]:=Symbol["d"<>SymbolName[Head@var]][]
ClearAll[ENstar]
SetAttributes[ENstar,HoldFirst]
ENstar[var_]:=Symbol[SymbolName[Head@var]<>"s"][]

Defining (Subscript[[.], +/-]) operators, see RS06 
Be very careful with these operators!!! They do not function very well with ToCanonical and should be removed before doing this

ClearAll[kplus,kminus,u0orb0]
u0orb0[u0[ind_]]=True;
u0orb0[b0[ind_]]=True;
Notation[(NotationTemplateTag[(SubPlus[[a_]])]) ⟸ (NotationTemplateTag[kplus[a_]])]
Notation[(NotationTemplateTag[(SubMinus[[a_]])]) ⟸ (NotationTemplateTag[kminus[a_]])]
kplus[a_+b_]:=kplus[a]+kplus[b]
kminus[a_+b_]:=kminus[a]+kminus[b]
kplus[rest_ num_/;(NumericQ[num] || PD[num,Xd[],Vd[i]]==0) ]:=num kplus[rest]
kminus[rest_ num_/;(NumericQ[num] || PD[num,Xd[],Vd[i]]==0) ]:=num kminus[rest]
kplus[f_ H_/;u0orb0[H]]:=kplus[f]kplus[H]
kminus[f_ H_/;u0orb0[H]]:=kminus[f]kminus[H]
kplus[f_ PD[H_/;u0orb0[H],Xd[],Vd[ind__]]]:=kplus[f]kplus[PD[H,Xd[],Vd[ind]]]
kminus[f_ PD[H_/;u0orb0[H],Xd[],Vd[ind__]]]:=kminus[f]kminus[PD[H,Xd[],Vd[ind]]]
kpmProductRule={kplus[F_]kminus[G_]kminus[PD[H_/;u0orb0[H],Xd[],Vd[ind_]]]:>-kplus[F]kminus[PD[G,Xd[],Vd[ind]]]kminus[H],
kplus[F_]kminus[G_]kminus[PD[H_/;u0orb0[H],Xd[],Vd[i1_,i2_]]]:>kplus[F]kminus[PD[G,Xd[],Vd[i1,i2]]]kminus[H],
kminus[F_]kplus[G_]kplus[PD[H_/;u0orb0[H],Xd[],Vd[ind_]]]:>-kminus[F]kplus[PD[G,Xd[],Vd[ind]]]kplus[H],
kminus[F_]kplus[G_]kplus[PD[H_/;u0orb0[H],Xd[],Vd[i1_,i2_]]]:>kminus[F]kplus[PD[G,Xd[],Vd[i1,i2]]]kplus[H]};

Function to expand functions of k as k+K/2, figuring out how much of this expansion needs to be kept (based on whether it involves (Subscript[B, i]) or (Subscript[B, ij])), actually

(* Rules and turnStatistic... are now a bit different *)
tov0b0h0Rules = {kplus[u0[i1_]]kminus[u0[i2_]]:>v0stat[i1,i2],kplus[b0[i1_]]kminus[b0[i2_]]:>b0stat[i1,i2],kplus[u0[i1_]]kminus[b0[i2_]]:>0,
kplus[b0[i1_]]kminus[u0[i2_]]:>0(* THIS MUST BE CHANGED IF NONZERO*)};

ClearAll[findBiBij](* Finds if term contains B or ∂B (or netiher)*)
findBiBij[expr_]:=If[!FreeQ[expr,Bt[__]],"Bgrad",
If[!FreeQ[expr,Bc[_]],"Bcon",
"noB"]]

ClearAll[turnIntoStatisticsIntegral]
turnIntoStatisticsIntegral::noB="No B found, setting term to zero";
turnIntoStatisticsIntegral[expr_]:=Module[{type},
type = findBiBij[expr];
Switch[type,
"Bgrad",
expandInK[(* Don't keep ∇(Subscript[v, ij ])if (Subscript[B, ij]) term *)expr/.tov0b0h0Rules/.fullStandardRuleSet//.{kplus[a_ b_]:>kplus[a]kplus[b],kminus[a_ b_]:>kminus[a]kminus[b]},0],
"Bcon",
expandInK[(* Keep ∇(Subscript[v, ij ])if (Subscript[B, i]) term *)
expr/.tov0b0h0Rules/.fullStandardRuleSet//.{kplus[a_ b_]:>kplus[a]kplus[b],kminus[a_ b_]:>kminus[a]kminus[b]},1],
_,
Message[turnIntoStatisticsIntegral::noB];Print[expr];0 (* neglect any contributions without (Subscript[B, i]). These are all proportional to <(Subscript[u, 0])(Subscript[b, 0])>*)]

]

(* Function to get to first order in K *)
expandInKRules0 = {kminus[v[i_]]:>-v[i],kplus[v[i_]]:>v[i],
kminus[km[]^n_.]:>km[]^n,kplus[km[]^n_.]:>km[]^n,
kminus[EN_^n_.]/;ENMember[EN]:>ENstar[EN]^n,
kplus[EN_^n_.]/;ENMember[EN]:>EN^n};
expandInKRules1 = {kplus[v[i_]]:>v[i]+K[i]/2,kminus[v[i_]]:>-v[i]+K[i]/2,kminus[km[]^n_.]:>km[]^n-n/2 dot[K,v]km[]^(n-2),kplus[km[]^n_.]:>km[]^n+n/2 dot[K,v]km[]^(n-2),kminus[EN_^n_.]/;ENMember[EN]:>ENstar[EN]^n-n/2dot[K,v]/km[]ENstar[EN]^(n-1) ENstar[dEN[EN]],
kplus[EN_^n_.]/;ENMember[EN]:>EN^n+n/2dot[K,v]/km[]EN^(n-1)dEN[EN]};
ClearAll[expandInK];
expandInK::noSecond = "Cannot go beyond first order expansion!";
expandInK[expr_,order_]:=
Switch[order,
0,expr/.expandInKRules0,
1,ExpandAll[expr/.expandInKRules1]/.K[i1_]K[i2_]->0,
_,Message[expandInK::noSecond]
]

ClearAll[calculateEN](* Given some product of E and N finds the symmetric (real) part in ω *)
(* Checked (k^(-υ)) rather than (k^(-2)) here. Did not give a positive (Subscript[((Subscript[η, yx])), u]) so τ-approximation must be more complex than this (or perhaps it can be very sensitive to differences in the exact perturbation expansion?)*)
calculateEN::notAntiSym="Imaginary part not antisymmetric in ω!";
calculateEN[expr_]:=Module[{tmp,EE=(1/(η (k^2)-I ω)),NN=(1/(ν (k^2)-I ω)),EEs=(1/(η (k^2)+I ω)),NNs=(1/(ν (k^2)+I ω)),impart},
tmp=expr/.{km[]->k,Eη[]->EE,Nν[]->NN,Eηs[]->EEs,Nνs[]->NNs,
dEη[]->D[#,k]&@EE,dNν[]->D[#,k]&@NN,dEηs[]->D[#,k]&@EEs,dNνs[]->D[#,k]&@NNs,
ddEη[]->D[#,k,k]&@EE,ddNν[]->D[#,k,k]&@NN,ddEηs[]->D[#,k,k]&@EEs,ddNνs[]->D[#,k,k]&@NNs};
impart=ComplexExpand[Im[tmp]];
If[Simplify[impart≠(impart/.ω->-ω)],Message[calculateEN::notAntiSym]];
ComplexExpand[Re[tmp]]//Together
]

Rules to use in Canonicalizing

(* Now don't store these as rules in ToCanonical, apply separately (because of kplus) *)
(*Occasionally need these rules separately *)
UsqΩsqBsqRules = {Bt[i1_,i1_]->0,Ut[i1_,i1_]->0,
Bt[i1_,i2_]Bt[i3_,i4_]->0,Bc[i1_]Bt[i3_,i4_]->0,Bc[i1_]Bc[i2_]->0,
Ut[i1_,i2_]Ut[i3_,i4_]->0,Ut[i1_,i2_]Ω[i3_]->0,Ω[i1_]Ω[i2_]->0,gρ[]gρ[]->0,gρ[]Ω[i2_]->0,Ut[i1_,i2_]gρ[]->0,
(*We neglect contributions to resistivity from stratification *)
gρ[]Bt[i1_,i2_]->0};
fullStandardRuleSet=Join[{PD[EN_,Xd[],Vd[i1_]]/;ENMember[EN]:>v[i1]/km[] dEN[EN],
PD[EN_,Xd[],Vd[i1_,i2_]]/;ENMember[EN]:>(delta[i1,i2]/km[]-v[i1]v[i2]/km[]^3) dEN[EN]+v[i1]v[i2]/km[]^2  dEN[dEN[EN]],
PD[km[],Xd[],Vd[i1_]]:>v[i1]/km[],
v[i1_]**v[i1_]:>km[]^2,
((v[i1_])^2):>km[]^2,
K[i1_]K[i2_]->0,((K[i1_])^2)->0},UsqΩsqBsqRules];

(* Function to apply divergence constraint*)
VijBijDivRules = {v[i1_]v0stat[i1_,i2_]:>-1/2 K[i1]v0stat[i1,i2],
v[i1_]v0stat[i2_,i1_]:>1/2 K[i1]v0stat[i2,i1],
v[i1_]b0stat[i1_,i2_]:>-1/2 K[i1]b0stat[i1,i2],
v[i1_]b0stat[i2_,i1_]:>1/2 K[i1]b0stat[i2,i1]};
ClearAll@applyVijBijDivRules
applyVijBijDivRules[expr_]:=expr/.VijBijDivRules/.Bt[i1_,i2_]K[i3_]->0

(* Rules to change to homogenous statistics *)(* Follows convention from Radler et al 2003 Eq.(71) *)
v0b0statToWH = {v0stat[i1_,i2_]:>1/2(delta[i1,i2]-v[i1]v[i2]/km[]^2-1/2/km[]^2(v[i1]K[i2]-v[i2]K[i1]))Wu[]-I/km[]^2Hu[]leviC[i1,i2,t]v[t],b0stat[i1_,i2_]:>1/2(delta[i1,i2]-v[i1]v[i2]/km[]^2-1/2/km[]^2(v[i1]K[i2]-v[i2]K[i1]))Wb[]-I/km[]^2 Hb[]leviC[i1,i2,t]v[t],
h0stat[i1_,i2_]->0};
(* NB: This neglects anything arising from ∇H *)
turnV0B0statToWH=ExpandAll[#/.v0b0statToWH]/.{K[i1_]Hu[]->0,K[i1_]Hb[]->0,K[i1_]K[i2_]->0}&;
integrationOverK::found6="Found set of 6 (Subscript[k, i]). Need to upgrade integration function!!!";
integrationOverK = #/.v[i1_]v[i2_]v[i3_]v[i4_]v[i5_]v[i6_]:>(Message[integrationOverK::found6];0)/.
v[i1_]v[i2_]v[i3_]v[i4_]v[i5_]:>0/.
v[i1_]v[i2_]v[i3_]v[i4_]:>1/15 km[]^4(delta[i1,i2]delta[i3,i4]+delta[i1,i3]delta[i2,i4]+delta[i1,i4]delta[i2,i3])/.
v[i1_]v[i2_]v[i3_]->0/.
v[i1_]v[i2_]:>1/3 km[]^2 delta[i1,i2]/.
v[i1_]->0&;

(* Wu is now (ρ^2)q, where q is velocity specturm (see Rudiger and Kitchanov). Thus ∇W gets a contribution from both of these bits. Of course, Wb is only turbulent part. *)
GWtoGuGρRules = {K[i1_]Wu[]:>-I(2 gρ[] gh[i1]Wu[]+2gt[] gh[i1]Wu[]),
K[i1_]Wb[]:>-I(2gt[] gh[i1]Wb[])};

Separating (Subscript[B, ij]) and (Subscript[U, ij]) into parts as in Radler

(* Just convenient way to set these rules when needed*)
setUBDecompositionRules:=Module[{},
DefObject[cB,1,{False,False},"∇\[Cross]B"];
DefObject[Bts,2,{False,False},"(B^((s)))",Symmetric[{1,2}]];
DefObject[Dts,2,{False,False},"(D^((s)))",Symmetric[{1,2}]];
DefObject[Wr,1,{False,False},"W"];
defineUserRules[{Ut[i1_,i2_]:>Dts[i1,i2]-1/2leviC[i1,i2,l]Wr[l],
Bt[i1_,i2_]:>Bts[i1,i2]-1/2leviC[i1,i2,l]cB[l],
Dts[i1_,i1_]->0,Bts[i1_,i1_]->0}]
]

(* Set rules to convert to actual, D ∇\[Cross]B, W\[Cross]∇\[Cross]B, κ∇B etc. Seemed easier to do in two steps.*)
setUBDecompositionRules2:=Module[{},
(* α coefficients *)
DefObject[DB,1,{False,False},"D∘B"];
DefObject[WcB,1,{False,False},"W\[Cross]B"];
DefObject[ΩcB,1,{False,False},"Ω\[Cross]B"];
DefObject[gdΩB,1,{False,False},"(g.Ω)B"];
DefObject[gdWB,1,{False,False},"(g.W)B"];
DefObject[ΩdBgPgdBΩ,1,{False,False},"[(Ω.B)g+(g.B)Ω]"];
DefObject[WdBgPgdBW,1,{False,False},"[(W.B)g+(g.B)W]"];
DefObject[αgDB,1,{False,False},"(OverHat[α])(g,D)∘B"];
DefObject[gcB,1,{False,False},"g\[Cross]B"];
DefObject[gcΩcB,1,{False,False},"(g\[Cross]Ω)\[Cross]B"];
DefObject[gcWcB,1,{False,False},"(g\[Cross]W)\[Cross]B"];
DefObject[gdDcB,1,{False,False},"(g∘D)\[Cross]B"];
(* β coefficients *)
DefObject[DcurlB,1,{False,False},"D∘∇\[Cross]B"];
DefObject[ΩccurlB,1,{False,False},"Ω\[Cross]∇\[Cross]B"];
DefObject[WccurlB,1,{False,False},"W\[Cross]∇\[Cross]B"];
DefObject[ΩgradB,1,{False,False},"Ω∘(∇B)^((s))"];
DefObject[WgradB,1,{False,False},"W∘(∇B)^((s))"];
DefObject[κgradB,1,{False,False},"κ(D)∘(∇B)^((s))"];
defineUserRules[{(* α coefficients *)
gh[i1_]Ω[i1_]Bc[i2_]:>gdΩB[i2],
gh[i1_]Wr[i1_]Bc[i2_]:>gdWB[i2],
gh[i2_]Ω[i1_]Bc[i1_]:>1/2ΩdBgPgdBΩ[i2]-1/2gcΩcB[i2],
gh[i2_]Wr[i1_]Bc[i1_]:>1/2WdBgPgdBW[i2]-1/2gcWcB[i2],
gh[i1_]Ω[i2_]Bc[i1_]:>1/2gcΩcB[i2]+1/2ΩdBgPgdBΩ[i2],
gh[i1_]Wr[i2_]Bc[i1_]:>1/2gcWcB[i2]+1/2WdBgPgdBW[i2],
(*NB: this uses the identity 
(Subscript[B, j]) (Subscript[(D^((s))), ik]) (Subscript[g, l]) (Subscript[ε, ljk]) - (Subscript[B, j]) (Subscript[(D^((s))), kl]) (Subscript[g, l]) (Subscript[ε, ijk]) + (Subscript[B, j]) (Subscript[(D^((s))), kj]) (Subscript[g, l]) (Subscript[ε, ilk])=0 *)
leviC[i1_,i2_,i3_]gh[i2_]Bc[i3_]:>gcB[i1],
leviC[i1_,i2_,i3_]gh[i3_]Bc[i2_]:>-gcB[i1],
leviC[i1_,i2_,i3_]Dts[i2_,i4_]gh[i4_]Bc[i3_]:>gdDcB[i1],
leviC[i1_,i2_,i3_]Dts[i4_,i2_]gh[i4_]Bc[i3_]:>gdDcB[i1],
leviC[i1_,i3_,i2_]Dts[i2_,i4_]gh[i4_]Bc[i3_]:>-gdDcB[i1],
leviC[i1_,i3_,i2_]Dts[i4_,i2_]gh[i4_]Bc[i3_]:>-gdDcB[i1],
leviC[i1_,i2_,i3_]Dts[i3_,i4_]gh[i2_]Bc[i4_]:>-1/2gdDcB[i1]-1/2αgDB[i1],
leviC[i4_,i2_,i3_]Dts[i1_,i2_]gh[i3_]Bc[i4_]:>-1/2gdDcB[i1]+1/2αgDB[i1],
(* β coefficients *)
Dts[i1_,i2_]cB[i2_]:>DcurlB[i1],Dts[i1_,i2_]cB[i1_]:>DcurlB[i2],
leviC[i1_,i2_,i3_]Ω[i2_]cB[i3_]:>ΩccurlB[i1],leviC[i1_,i2_,i3_]Ω[i3_]cB[i2_]:>-ΩccurlB[i1],
leviC[i1_,i2_,i3_]Wr[i2_]cB[i3_]:>WccurlB[i1],leviC[i1_,i2_,i3_]Wr[i3_]cB[i2_]:>-WccurlB[i1],
Ω[i1_]Bts[i1_,i2_]:>ΩgradB[i2],Ω[i2_]Bts[i1_,i2_]:>ΩgradB[i1],
Wr[i1_]Bts[i1_,i2_]:>WgradB[i2],Wr[i2_]Bts[i1_,i2_]:>WgradB[i1],
leviC[i1_,i2_,i3_]Dts[i3_,i4_]Bts[i2_,i4_]:>κgradB[i1]/2,
leviC[i1_,i2_,i3_]Dts[i4_,i3_]Bts[i2_,i4_]:>κgradB[i1]/2,
leviC[i1_,i2_,i3_]Dts[i4_,i3_]Bts[i4_,i2_]:>1/2κgradB[i1]}]
]
setUBDecompositionRulesHelical:=defineUserRules[(* α coefficients *)
{Dts[i1_,i2_]Bc[i2_]:>DB[i1],Dts[i2_,i1_]Bc[i2_]:>DB[i1],
leviC[i1_,i2_,i3_]Wr[i2_]Bc[i3_]:>WcB[i1],leviC[i1_,i2_,i3_]Wr[i3_]Bc[i2_]:>-WcB[i1],
leviC[i1_,i2_,i3_]Ω[i2_]Bc[i3_]:>ΩcB[i1],leviC[i1_,i2_,i3_]Ω[i3_]Bc[i2_]:>-ΩcB[i1]}]
  
ClearAll[transportCoeff]
transportCoeff[expr_,term_]:=calculateEN[Plus@@Cases[expr,a_ term:>a,∞]]

Integrate a specific form of W to obtain idea of signs

assumpList=ν>0&&η>0&&τc>0&&λc>0&&q>0&&Pm>0&&Rm>0&&St>0&&u2>0&&b2>0&&ϵ>0&&kq>0;
Notation[(NotationTemplateTag[(Subscript[τ, c])]) ⟸ (NotationTemplateTag[τc])]
Notation[(NotationTemplateTag[(Subscript[λ, c])]) ⟸ (NotationTemplateTag[λc])]
Notation[(NotationTemplateTag[⟨(u^2)⟩]) ⟸ (NotationTemplateTag[u2])]
Notation[(NotationTemplateTag[⟨(b^2)⟩]) ⟸ (NotationTemplateTag[b2])]
Wexp = u2((2((λc)^3)τc)/(3 (((2π))^(5/2))))(((((k λc))^2)Exp[-(((k λc))^2)/2])/((1+(((ω τc))^2))(*(1+(((k λc))^(kq+4)))*)));(* Is this correct? Where does (k^2) come from?*)
Hexp = b2((2((λc)^3)τc)/(3 (((2π))^(5/2))))(((((k λc))^2)Exp[-(((k λc))^2)/2])/((1+(((ω τc))^2))(*(1+(((k λc))^(kq+4)))*)));
ClearAll[integrateOverω,integrateOverk,simpFunc]
simpFunc=Simplify[#,Assumptions->assumpList]&;
integrateOverω[expr_]:=Integrate[expr  4π (k^2),{ω,-∞,∞},Assumptions->assumpList&&k>0]
integrateOverk[expr_]:=Integrate[expr ,{k,0,∞},Assumptions->assumpList]
(* Rules to non-dimensionalize *)
nonDimensionalizePrelim=Expand[Expand[#]/.λc->(Sqrt[q η τc])/.ν->Pm η/.τc->Rm St η/u2/.b2->ϵ u2]&;
nonDimensionalizeFinal = simpFunc[#]/.(St^2)(η^2)/u2->((λc)^2)/(q^2)/.u2 ϵ->b2&;

Calculation

Radler calculation to check everything is working

(u^((1))) etc. definitions

(* This is Radler (u^((1))),(b^((0))), and (b^((1))) *)
subrules = {u10->u1sub,b10->radb0};
u1sub (Overscript[=, ind])Nν[](-Ut[i,l]u0[l]+Ut[l,k]v[l]PD[u0[i],Xd[],Vd[k]]+2 v[i]v[j]/km[]^2u0[l]Ut[j,l]+2v(OverVector[•])Ω/km[]^2 leviC[i,j,k]u0[j]v[k])

radb0(Overscript[=, ind])Eη[](I v(OverVector[•])Bc u0[i]-Bt[i,j]u0[j]-Bt[j,k]v[j]PD[u0[i],Xd[],Vd[k]])
radb1(Overscript[=, ind])Evaluate[Eη[](I v(OverVector[•])Bc u10[i]-Bt[i,j]u10[j]-Bt[j,k]v[j]PD[u10[i],Xd[],Vd[k]]-(-Ut[i,j]b10[j]-Ut[j,k]v[j]PD[b10[i],Xd[],Vd[k]]))/.subrules//ToCanonical]

Main calculation

Expand ℰ = (ℰ^((00)))+(ℰ^((10)))+(ℰ^((01)))

defineUserRules[{}];
u0transport=ToCanonical[
turnIntoStatisticsIntegral/@
ExpandAll[PDExpand[ExpandAll[PDExpand[
Expand[
leviC[i,j,k]
(kplus[u0[j]]kminus[radb0[k]]+kplus[u0[j]]kminus[radb1[k]]+kplus[u1sub[j]]kminus[radb0[k]])
]/.UsqΩsqBsqRules/.kpmProductRule]]/.kpmProductRule]]]/.VijBijDivRules/.fullStandardRuleSet//ToCanonical;

Double check that it matches Radler Eq.(44), could do others too

(*setStandardUserRules;*)
DefObject[kdOok2,0,{False,True},"((k . Ω)/(k^2))"]
(*DefObject[v0statS,2,{False,False},"((OverTilde[v])^((s)))",Symmetric[{1,2}]]
DefObject[v0statA,2,{False,False},"((OverTilde[v])^((A)))",Antisymmetric[{1,2}]]*)
ToCanonical[ToCanonical[u0transport/.Ut[_,_]->0]/.{v[i1_]Ω[i1_]:>kdOok2 km[]^2}];(*/.{v0stat[i1_,i2_]:>v0statS[i1,i2]+v0statA[i1,i2]}]/.v0statA[i1_,i2_]->0*)
(*v0statA->0 Because always occurs with an even number of ks -> disappears upon integration *)

Final answer - plug in (Subscript[((v^((0)))), ij ])= (Subscript[P, ij])W(k,ω) then integrate over k space using Radler Eq.(48)

(*Have to plug in the rules, applying ToCanonical,  a couple of times. Once for integration over k, then again to get (Subscript[U, jj])=0, (Subscript[K, i])(Subscript[K, i])=0 etc.*)
defineUserRules[{}];
finalans=ToCanonical[
ToCanonical[
integrationOverK[ToCanonical[ToCanonical[
turnV0B0statToWH/@u0transport]/.fullStandardRuleSet]]
]/.fullStandardRuleSet]

Plug in Radler’s decomposition of (Subscript[U, ij]) and (Subscript[B, ij]) tensors to look at shear-current effect and Ω×J effect

(* This is basic set in terms of ∇×B, Ω, W, D, (B^((s))) etc. *)
setUBDecompositionRules;
finalansTrans=finalans/.K[i1_]:>-I g[i1]//ToCanonical//ToCanonical;
(* This is set as written in Eq. (10), D∘∇×B, W×∇×B, etc. *)
setUBDecompositionRules2;
finalansTrans=finalansTrans//ToCanonical//ToCanonical;
setUBDecompositionRulesHelical;
finalansTrans=finalansTrans//ToCanonical

Find all the transport coefficients directly as given in Radler (mainly in the appendix)
NB: Negative sign in his definitions in Eq.(10)

All are the same as Radler except for (β^((D))), which has a differnent denominator.

coefficientSimplify = Simplify[#,Assumptions->k∈Reals&&k>0]&
Print["β coefficients"]
β0r=transportCoeff[-finalansTrans,Wu[]cB[i]] (* (β^((0))) *)
δΩr=transportCoeff[-finalansTrans,Wu[]ΩccurlB[i]]//coefficientSimplify (* (δ^((Ω))) *)
δWr=transportCoeff[-finalansTrans,Wu[]WccurlB[i]]//coefficientSimplify (* (δ^((W))) *)
κΩr=transportCoeff[-finalansTrans,Wu[]ΩgradB[i]]//coefficientSimplify(* (κ^((Ω))) *)
κWr=transportCoeff[-finalansTrans,Wu[]WgradB[i]] //coefficientSimplify(* (κ^((W))) *)
βDr=transportCoeff[-finalansTrans,Wu[]DcurlB[i]]//coefficientSimplify (* (β^((D))) *)
κDr=transportCoeff[-finalansTrans,Wu[]κgradB[i]]//coefficientSimplify (* (κ^((D))) *)
Print["α coefficients"]
γ0r = transportCoeff[-finalansTrans,Wu[]gcB[i]]//coefficientSimplify (*(γ^((0)))*)
αΩ1r = transportCoeff[-finalansTrans,Wu[]gdΩB[i]]//coefficientSimplify (* (α^(1(Ω)))*)
αΩ2r = transportCoeff[-finalansTrans,Wu[]ΩdBgPgdBΩ[i]]//coefficientSimplify (* (α^(2(Ω)))*)
αW1r = transportCoeff[-finalansTrans,Wu[]gdWB[i]]//coefficientSimplify (* (α^(1(W)))*)
αW2r = transportCoeff[-finalansTrans,Wu[]WdBgPgdBW[i]]//coefficientSimplify (* (α^(2(W)))*)
γΩr = transportCoeff[-finalansTrans,Wu[]gcΩcB[i]]//coefficientSimplify (*(γ^((Ω)))*)
γWr = transportCoeff[-finalansTrans,Wu[]gcWcB[i]]//coefficientSimplify (*(γ^((W)))*)
αDr = transportCoeff[-finalansTrans,Wu[]αgDB[i]]//coefficientSimplify (*(α^((D)))*)
γDr = transportCoeff[-finalansTrans,Wu[]gdDcB[i]]//coefficientSimplify (*(γ^((D)))*)

SCfull=δWr-1/2(κWr-βDr+κDr)/.{η->ηt/k^2,ν->νt/k^2}//Expand//Together
SCyfull=δWr-1/2(κWr+βDr-κDr)/.{η->ηt/k^2,ν->νt/k^2}//Expand//Together
numerator=-4 ((ηt)^2)(ω^2)(((ηt)^2)+(ω^2))+(3((ηt)^4)+6((ηt)^2)(ω^2)-5(ω^4))(((νt)^2)+(ω^2));(* Corrected version of D5 *)denomenator=15 (((((ηt)^2)+(ω^2)))^3) (((νt)^2)+(ω^2));
SCyfull/.Expand[numerator]->numerator/.(a_+b_)/denomenator:>a/denomenator+b/denomenator
-δΩr+1/2(κΩr)/.{η->ηt/k^2,ν->νt/k^2}//FullSimplify

δWr-1/2(κWr-βDr/2+κDr)/.{η->ηt/k^2,ν->νt/k^2}//Expand//Together
Expand[8(4 ((ηt)^2)(ω^2)(((ηt)^2)+(ω^2))(((νt)^2)+(ω^2))+(((ηt)^4)+12((ηt)^2)(ω^2)-5(ω^4))(((((νt)^2)+(ω^2)))^2)/4)](* NOT THE SAME!!! *)

Adding (Subscript[b, 0]) also

Full calculation, including small-scale (Subscript[b, rms]), stratification of density, and turbulence. 
Note that we use the notation (Subscript[u, i]) here, but this is infact (Subscript[m, i]) in the paper (i.e., ρ (Subscript[u, i]))

(u^((1))) etc. definitions

subrules = {u10->u0sub,b10->b0sub};
(* Momentum equation same as the incompressible case (with u->m), aside from a single viscosity contribution *)
u0sub (Overscript[=, ind])Nν[](
(-I ν v(OverVector[•])gh gρ u0[i]+I ν v(OverVector[•])gh gρ v[i]v[j]/km[]^2u0[j])+
(-Ut[i,l]u0[l]+Ut[l,k]v[l]PD[u0[i],Xd[],Vd[k]]+2 v[i]v[j]/km[]^2u0[l]Ut[j,l])+
(I v(OverVector[•])Bc b0[i]-I v(OverVector[•])Bc v[i]v[j]/km[]^2b0[j])-(-Bt[i,l]b0[l]+Bt[l,k]v[l]PD[b0[i],Xd[],Vd[k]]+2 v[i]v[j]/km[]^2b0[l]Bt[j,l])+2v(OverVector[•])Ω/km[]^2 leviC[i,j,k]u0[j]v[k]
)//ToCanonical
u1sub (Overscript[=, ind])Evaluate[ToCanonical[
Nν[](
(-I ν v(OverVector[•])gh gρ u10[i]+I ν v(OverVector[•])gh gρ v[i]v[j]/km[]^2u10[j])+
(-Ut[i,l]u10[l]+Ut[l,k]v[l]PD[u10[i],Xd[],Vd[k]]+2 v[i]v[j]/km[]^2u10[l]Ut[j,l])+
(I v(OverVector[•])Bc b10[i]-I v(OverVector[•])Bc v[i]v[j]/km[]^2b10[j])-(-Bt[i,l]b10[l]+Bt[l,k]v[l]PD[b10[i],Xd[],Vd[k]]+2 v[i]v[j]/km[]^2b10[l]Bt[j,l])+2v(OverVector[•])Ω/km[]^2 leviC[i,j,k]u10[j]v[k]
)/.subrules]/.UsqΩsqBsqRules]

b0sub(Overscript[=, ind])Eη[](
I  ρ0m1[] v(OverVector[•])Bc u0[i]-ρ0m1[]Bt[i,j]u0[j]-ρ0m1[]Bt[j,k]v[j]PD[u0[i],Xd[],Vd[k]]-(-Ut[i,j]b0[j]-Ut[j,k]v[j]PD[b0[i],Xd[],Vd[k]])+
ρ0m1[]gρ[] Bc[i]gh[r]u0[r]+ρ0m1[]gρ[] gh[r]Bc[j]v[j]PD[u0[i],Xd[],Vd[r]]
)//ToCanonical
b1sub(Overscript[=, ind])Evaluate[ToCanonical[
Eη[](
I ρ0m1[] v(OverVector[•])Bc u10[i]-ρ0m1[]Bt[i,j]u10[j]-ρ0m1[]Bt[j,k]v[j]PD[u10[i],Xd[],Vd[k]]-(-Ut[i,j]b10[j]-Ut[j,k]v[j]PD[b10[i],Xd[],Vd[k]])+
ρ0m1[]gρ Bc[i]gh[r]u10[r]+ρ0m1[]gρ[] gh[r]Bc[j]v[j]PD[u10[i],Xd[],Vd[r]]
)/.subrules]/.UsqΩsqBsqRules]

Main calculation

Expand ℰ = (ℰ^((00)))+(ℰ^((10)))+(ℰ^((01)))

defineUserRules[{}];
u0transport=ToCanonical[
turnIntoStatisticsIntegral/@
ExpandAll[PDExpand[ExpandAll[PDExpand[
Expand[
leviC[i,j,k]
(kplus[u0[j]]kminus[b0[k]]+kplus[u0[j]]kminus[b0sub[k]+b1sub[k]]+kplus[u0sub[j]+u1sub[j]]kminus[b0[k]]+kplus[u0sub[j]]kminus[b0sub[k]])
]/.UsqΩsqBsqRules/.kpmProductRule]]/.kpmProductRule]]]/.VijBijDivRules/.fullStandardRuleSet//ToCanonical;

Final answer - plug in (Subscript[((v^((0)))), ij ])= (Subscript[P, ij])W(k,ω) then integrate over k space using Radler Eq.(48)

(*Have to plug in the rules, applying ToCanonical,  a couple of times. Once for integration over k, then again to get (Subscript[U, jj])=0, (Subscript[K, i])(Subscript[K, i])=0 etc.*)
defineUserRules[{}];
finalans=ToCanonical[
ToCanonical[
integrationOverK[ToCanonical[ToCanonical[
turnV0B0statToWH/@u0transport]/.fullStandardRuleSet]]
]/.fullStandardRuleSet];

Plug in Radler’s decomposition of (Subscript[U, ij]) and (Subscript[B, ij]) tensors to look at shear-current effect and Ω×J effect

(* This is basic set in terms of ∇×B, Ω, W, D, (B^((s))) etc. *)
setUBDecompositionRules;
finalansTrans=finalans/.GWtoGuGρRules//ToCanonical//ToCanonical;
(* This is set as written in Eq. (10), D∘∇×B, W×∇×B, etc. *)
setUBDecompositionRules2;
finalansTrans=finalansTrans//ToCanonical//ToCanonical;
setUBDecompositionRulesHelical;
finalansTrans=finalansTrans//ToCanonical

Transport coefficients

Find all the transport coefficients directly as given in Radler (mainly in the appendix)
NB: Negative sign in his definitions in Eq.(10)

coefficientSimplify = Simplify[#,Assumptions->k∈Reals&&k>0]&;
Print["β coefficients"]
Print["((β^((0))))_u"];β0u=transportCoeff[-finalansTrans,Wu[]cB[i]]/((ρ0m1[])^2)//coefficientSimplify (* (β^((0))) *)
Print["((β^((0))))_b"];β0b=transportCoeff[-finalansTrans,Wb[]cB[i]]//coefficientSimplify
Print["((δ^((Ω))))_u"];δΩu=transportCoeff[-finalansTrans,Wu[]ΩccurlB[i]]/((ρ0m1[])^2) //coefficientSimplify(* (δ^((Ω))) *)
Print["((δ^((Ω))))_b"];δΩb=transportCoeff[-finalansTrans,Wb[]ΩccurlB[i]]//coefficientSimplify
Print["((δ^((W))))_u"];δWu=transportCoeff[-finalansTrans,Wu[]WccurlB[i]]/((ρ0m1[])^2) //coefficientSimplify(* (δ^((W))) *)
Print["((δ^((W))))_b"];δWb=transportCoeff[-finalansTrans,Wb[]WccurlB[i]]//coefficientSimplify
Print["((κ^((Ω))))_u"];κΩu=transportCoeff[-finalansTrans,Wu[]ΩgradB[i]]/((ρ0m1[])^2)//coefficientSimplify(* (κ^((Ω))) *)
Print["((κ^((Ω))))_b"];κΩb=transportCoeff[-finalansTrans,Wb[]ΩgradB[i]]//coefficientSimplify
Print["((κ^((W))))_u"];κWu=transportCoeff[-finalansTrans,Wu[]WgradB[i]]/((ρ0m1[])^2) //coefficientSimplify(* (κ^((W))) *)
Print["((κ^((W))))_b"];κWb=transportCoeff[-finalansTrans,Wb[]WgradB[i]]//coefficientSimplify
Print["((β^((D))))_u"];βDu=transportCoeff[-finalansTrans,Wu[]DcurlB[i]]/((ρ0m1[])^2)//coefficientSimplify (* (β^((D))) *)
Print["((β^((D))))_b"];βDb=transportCoeff[-finalansTrans,Wb[]DcurlB[i]]//coefficientSimplify
Print["((κ^((D))))_u"];κDu=transportCoeff[-finalansTrans,Wu[]κgradB[i]]/((ρ0m1[])^2)//coefficientSimplify (* (κ^((D))) *)
Print["((κ^((D))))_b"];κDb=transportCoeff[-finalansTrans,Wb[]κgradB[i]]//coefficientSimplify
Print["α coefficients"]
Print["((γ^((0))))_u"];
γ0u = transportCoeff[-finalansTrans,Wu[]gcB[i]]/((ρ0m1[])^2)//coefficientSimplify (*(γ^((0)))*)
Print["((γ^((0))))_b"];
γ0b = transportCoeff[-finalansTrans,Wb[]gcB[i]]//coefficientSimplify (*(γ^((0)))*)
Print["((α^((Ω) 1)))_u"]
αΩ1u = transportCoeff[-finalansTrans,Wu[]gdΩB[i]]/((ρ0m1[])^2)//coefficientSimplify (* (α^(1(Ω)))*)
Print["((α^((Ω) 1)))_b"]
αΩ1b = transportCoeff[-finalansTrans,Wb[]gdΩB[i]]//coefficientSimplify (* (α^(1(Ω)))*)
Print["((α^((Ω) 2)))_u"]
αΩ2u= transportCoeff[-finalansTrans,Wu[]ΩdBgPgdBΩ[i]]/((ρ0m1[])^2)//coefficientSimplify (* (α^(2(Ω)))*)
Print["((α^((Ω) 2)))_b"]
αΩ2b= transportCoeff[-finalansTrans,Wb[]ΩdBgPgdBΩ[i]]//coefficientSimplify (* (α^(2(Ω)))*)
Print["((α^((W) 1)))_u"]
αW1u = transportCoeff[-finalansTrans,Wu[]gdWB[i]]/((ρ0m1[])^2)//coefficientSimplify (* (α^(1(W)))*)
Print["((α^((W) 1)))_b"]
αW1b = transportCoeff[-finalansTrans,Wb[]gdWB[i]]//coefficientSimplify (* (α^(1(W)))*)
Print["((α^((W) 2)))_u"]
αW2u = transportCoeff[-finalansTrans,Wu[]WdBgPgdBW[i]]/((ρ0m1[])^2)//coefficientSimplify (* (α^(2(W)))*)
Print["((α^((W) 2)))_b"]
αW2b = transportCoeff[-finalansTrans,Wb[]WdBgPgdBW[i]]//coefficientSimplify (* (α^(2(W)))*)
Print["((γ^((Ω))))_u"]
γΩu = transportCoeff[-finalansTrans,Wu[]gcΩcB[i]]/((ρ0m1[])^2)//coefficientSimplify (*(γ^((Ω)))*)
Print["((γ^((Ω))))_b"]
γΩb = transportCoeff[-finalansTrans,Wb[]gcΩcB[i]]//coefficientSimplify (*(γ^((Ω)))*)
Print["((γ^((W))))_u"];
γWu = transportCoeff[-finalansTrans,Wu[]gcWcB[i]]/((ρ0m1[])^2)//coefficientSimplify (*(γ^((W)))*)
Print["((γ^((W))))_b"];
γWb = transportCoeff[-finalansTrans,Wb[]gcWcB[i]]//coefficientSimplify (*(γ^((W)))*)
Print["((α^((D))))_u"]
αDu = transportCoeff[-finalansTrans,Wu[]αgDB[i]]/((ρ0m1[])^2)//coefficientSimplify (*(α^((D)))*)
Print["((α^((D))))_b"]
αDb = transportCoeff[-finalansTrans,Wb[]αgDB[i]]//coefficientSimplify (*(α^((D)))*)
Print["((γ^((D))))_u"];
γDu= transportCoeff[-finalansTrans,Wu[]gdDcB[i]]/((ρ0m1[])^2)//coefficientSimplify (*(γ^((D)))*)
Print["((γ^((D))))_b"];
γDb= transportCoeff[-finalansTrans,Wb[]gdDcB[i]]//coefficientSimplify (*(γ^((D)))*)
Print["Helical α coefficients"]
Print["((α^(H (0))))_u"];αH0u=transportCoeff[-finalansTrans,Hu[]Bc[i]]/((ρ0m1[])^2)//coefficientSimplify (* (α^(H(0))) *)
Print["((α^(H (0))))_b"];αH0b=transportCoeff[-finalansTrans,Hb[]Bc[i]]//coefficientSimplify (* (α^(H(0))) *)
Print["((γ^(H (Ω)))\!)_u"];αHΩu=transportCoeff[-finalansTrans,Hu[]ΩcB[i]]/((ρ0m1[])^2)//coefficientSimplify (* (γ^(H(Ω))) *)
Print["((γ^(H (Ω)))\!)_b"];αHΩb=transportCoeff[-finalansTrans,Hb[]ΩcB[i]]//coefficientSimplify (* (γ^(H(Ω))) *)
Print["((γ^(H (W))))_u"];αHWu=transportCoeff[-finalansTrans,Hu[]WcB[i]]/((ρ0m1[])^2)//coefficientSimplify (* (γ^(H(W))) *)
Print["((γ^(H (W))))_b"];αHWb=transportCoeff[-finalansTrans,Hb[]WcB[i]]//coefficientSimplify (* (γ^(H(W))) *)
Print["((α^(H (D))))_u"];αHDu=transportCoeff[-finalansTrans,Hu[]DB[i]]/((ρ0m1[])^2)//coefficientSimplify (* (α^(H(D))) *)
Print["((α^(H (D))))_b"];αHDb=transportCoeff[-finalansTrans,Hb[]DB[i]]//coefficientSimplify (* (α^(H(D))) *)

ClearAll[printInTex]
Options[printInTex]={simplify->Simplify};
printInTex[expr_,OptionsPattern[]]:=Module[{wett},
wett=OptionValue[simplify][ expr/.{η->(OverTilde[η])/(k^2),ν->(OverTilde[ν])/(k^2)}];
TeXForm[wett]
]

γWu/.{η->(OverTilde[η])/(k^2),ν->(OverTilde[ν])/(k^2)}//Simplify

printInTex[γWu]

Cartesian shear dynamo from the shear-current/Ω×J effect

Shear-current and Ω×J in Cartesian geometry
SWtotalu is always negative (negative definite integrand), but δWu is always positive (use integration by parts). So W×J effect itself is positive, other bits cause the sign change. δWb is very similar to δWu and positive.

SWtotalu=δWu-1/2(κWu-βDu+κDu)//FullSimplify
SWtotalb=δWb-1/2(κWb-βDb+κDb)//FullSimplify
SΩtotalu=-δΩu+1/2κΩu//FullSimplify
SΩtotalb=-δΩb+1/2κΩb//FullSimplify

nonDimensionalizeFinal[nonDimensionalizePrelim[integrateOverk[integrateOverω[δΩu Wexp/.ν->η ]]]]//Expand//simpFunc//nonDimensionalizeFinal
nonDimensionalizeFinal[nonDimensionalizePrelim[integrateOverk[integrateOverω[δΩb Wexp/.ν->η ]]]]

Plot[{-(1/72) (Sqrt[q]) (-2 (Sqrt[q])+(E^(q/2)) (Sqrt[2 π]) (1+q) Erfc[((Sqrt[q])/(Sqrt[2]))]),-(1/72) (-2 (2+q)+(E^(q/2)) (Sqrt[2 π]) (Sqrt[q]) (3+q) Erfc[((Sqrt[q])/(Sqrt[2]))])},{q,0.01,10}]

Carry out integrals over Gaussian (Subscript[W, u]) and (Subscript[W, b])

SΩtotaluInt=nonDimensionalizeFinal[nonDimensionalizePrelim[integrateOverk[integrateOverω[SΩtotalu Wexp]]]]
SΩtotalbInt=nonDimensionalizeFinal[nonDimensionalizePrelim[integrateOverk[integrateOverω[SΩtotalb Hexp]]]]

SWtotaluInt=nonDimensionalizeFinal[nonDimensionalizePrelim[integrateOverk[integrateOverω[SWtotalu Wexp]]]]
SWtotalbInt=nonDimensionalizeFinal[nonDimensionalizePrelim[integrateOverk[integrateOverω[SWtotalb Hexp]]]]//Expand//simpFunc//nonDimensionalizeFinal

prec=30;
listatPm[Pmp_]:=(
sΩfu = Function[q$,Evaluate[SΩtotaluInt/((Rm^2) ((λc)^2) )/.q->q$/.Pm->Pmp//N[#,prec]&//Expand]];
sΩfb = Function[q$,Evaluate[SΩtotalbInt/((Rm^2) ((λc)^2) ϵ)/.q->q$/.Pm->Pmp//N[#,prec]&//Expand]];
qplt = 10^N[#,prec]&/@Range[-2,2,1/20];
)

Show[(listatPm[101/100];ListLogLinearPlot[{{qplt,sΩfu/@qplt}\[Transpose],{qplt,sΩfb/@qplt}\[Transpose]},Joined->True,PlotStyle->Thick(*,PlotRange->{{(10^(-2)),(10^2)},{(10^(-5)),(10^1)}}*),AxesLabel->{q,(Subscript[η, 21])}]),(listatPm[2];ListLogLinearPlot[{{qplt,sΩfu/@qplt}\[Transpose],{qplt,sΩfb/@qplt}\[Transpose]},Joined->True,PlotStyle->Dashed]),
(listatPm[1/2];ListLogLinearPlot[{{qplt,sΩfu/@qplt}\[Transpose],{qplt,sΩfb/@qplt}\[Transpose]},Joined->True,PlotStyle->Dotted])]

Simplifying the (Subscript[b, rms]) shear-current dynamo and proving positivity at Pm=1

ans=SWtotalb/.{η->ηt/(k^2),ν->νt/(k^2)}//Together
{{νη,η2,ν2},rem}=PolynomialReduce[Numerator@ans,{(((ηt)^2)+(ω^2))(((νt)^2)+(ω^2)),(((((ηt)^2)+(ω^2)))^2),(((((νt)^2)+(ω^2)))^2)},{ηt,νt,ω}]
simpSWtotalb=η2 /(Simplify[Denominator[ans]/(((((ηt)^2)+(ω^2)))^2)])+ν2 /(Simplify[Denominator[ans]/(((((νt)^2)+(ω^2)))^2)])+νη /(Simplify[Denominator[ans]/(((ηt)^2)+(ω^2))/(((νt)^2)+(ω^2))])
FullSimplify[ans==%]

simpSWtotalb=-((8 (ω^4))/(15 (((((νt)^2)+(ω^2)))^3)))-((8 ηt νt (ω^2))/(15 (((((ηt)^2)+(ω^2)))^2) (((νt)^2)+(ω^2))))+((2 ((ηt)^2) ((νt)^2)+4 ηt ((νt)^3)+4 ((ηt)^2) (ω^2)+6 (ω^4))/(15 (((ηt)^2)+(ω^2)) (((((νt)^2)+(ω^2)))^2)));
simpSWtotalbPm1=simpSWtotalb/.νt->ηt//FullSimplify
-Integrate[%,ω](* Negative with integration by parts, so ω dW/dω≤0 & ArcTan[ω]dW/dω≤0 implies this is positive *)

Series[simpSWtotalb/.νt->Pm ηt,{Pm,∞,2}]
-Integrate[%,ω]//Expand
Series[simpSWtotalb/.ηt->oPm νt,{oPm,∞,1}]
-Integrate[%,ω]//Expand

Integrate[((2 ((k^8) (ν^4)-12 (k^4) (ν^2) (ω^2)+3 (ω^4)))/(15 ((((k^4) (ν^2)+(ω^2)))^3))),ω]

Full growth rate

Last@Eigenvalues[(({{-ηyy (k^2),  η21 (k^2)}, {-S + η12 (k^2), -ηxx (k^2)}}))]//Simplify

(Sqrt[(k^2) S η21+(k^4) η12 η21])-(k^2) ηt/.k->((S η21)/(4 ηt))/.η12->0//FullSimplify

Cartesian shear dynamo from the α effect due to vertical stratification

αyyWtotalu=αDu-αW1u(*/.η->ν*)//FullSimplify(*<0 probably, is at Pm=1*)
αyyWtotalb=αDb-αW1b(*/.η->ν*)//FullSimplify(*>0 probably, is at Pm=1*)
αyyΩtotalu=αΩ1u (*/.η->ν*)//FullSimplify(*>0*)
αyyΩtotalb=αΩ1b(*/.η->ν*)//FullSimplify(*<0*)

(* Negativity with integration by parts impies integral is positive so long as ω dW/dω≤0 & ArcTan[ω]dW/dω≤0 *)
-Integrate[αyyWtotalb/.ν->η,ω]//FullSimplify
-Integrate[αyyΩtotalb,ω]

nonDimensionalizeFinal[nonDimensionalizePrelim[integrateOverk[integrateOverω[αyyWtotalu Wexp/.ν->η ]]]]

nonDimensionalizeFinal[nonDimensionalizePrelim[integrateOverk[integrateOverω[δΩu Wexp/.ν->η ]]]]//Expand//simpFunc//nonDimensionalizeFinal
nonDimensionalizeFinal[nonDimensionalizePrelim[integrateOverk[integrateOverω[δΩb Wexp/.ν->η ]]]]

Coefficients for the shear dynamo problems

Mean-fields depend only on z, U is just a basic linear shear, Ω is in the (OverHat[z]) direction, and there is vertical stratification (decreasing turbulence intensity away from the midplane). We look above the midplane here so (1/(u^2))((∂(u^2))/(∂z))=g (OverHat[z]) with g<0, and for simplicity we will assume (Subscript[g, u])=(Subscript[g, b]).

(* Define the various vectors and tensors for an accretion disk style dynamo *) 
xv={x,y,z};
Bv=Through[{Bx,By,0&}@z];
Uv={0,-S x,0};
Bij=Grad[Bv,xv];BSij = Symmetrize[Bij,Symmetric[{1,2}]]//Normal;curlB=Curl[Bv,xv];
Uij = Grad[Uv,xv];Dij=Symmetrize[Uij,Symmetric[{1,2}]]//Normal;Wi = Curl[Uv,xv];Ωi={0,0,Ω};
gi = {0,0,-gt};

Resistivity terms, there is no stratification in here

(βContribution=Curl[
-β0 curlB -βD Sum[Dij[[All,j]]curlB[[j]],{j,1,3}]-(δΩ Ωi+δW Wi)\[Cross]curlB-Sum[(κΩ Ωi+κW Wi)[[j]] BSij[[j,All]],{j,3}]-κD (Sum[LeviCivitaTensor[3][[All,j,l]]Dij[[l,k]]BSij[[j,k]]+LeviCivitaTensor[3][[All,k,l]]Dij[[l,j]]BSij[[j,k]],{j,3},{l,3},{k,3}]),xv]//Collect[#,{Bx''[z],By''[z]},Collect[#,{S,Ω}]&]&)//ColumnForm

α and γ terms, these rely on the turbulence stratification

(αContribution=Curl[
-αΩ1(gi.Ωi)Bv-αW1(gi.Wi)Bv-αΩ2((Ωi.Bv)gi+(gi.Bv)Ωi)-αW2((Wi.Bv)gi+(gi.Bv)Wi)-αD Sum[LeviCivitaTensor[3][[All,l,m]]Dij[[l,j]]gi[[m]]Bv[[j]]+LeviCivitaTensor[3][[j,l,m]]Dij[[l,All]]gi[[m]]Bv[[j]],{j,3},{l,3},{m,3}]-
(γ0 gi+γΩ gi\[Cross]Ωi+γW gi\[Cross]Wi+γD Sum[Dij[[All,j]]gi[[j]],{j,1,3}])\[Cross]Bv,
xv]//Collect[#,{Bx'[z],By'[z]},Collect[#,{S,Ω}]&]&)//ColumnForm
αContribution/.Ω->2/3 S//Collect[#,{Bx'[z],By'[z]},Simplify]&//ColumnForm

Curl[Table[(Subscript[a, i,j]),{i,3},{j,3}].Bv,xv]

Eigenvalues of the α effect dynamo, include diagonal resistivity

(({{-I k gt γ -(k^2)ηt, -I k αyy}, {I k αxx-S, -I k gt γ -(k^2)ηt}}))//Eigenvalues
(* To put this in a more easily readable form *)
Γ==((Sqrt[S k \[LeftBracketingBar](Subscript[α, yy])\[RightBracketingBar]/2])-(k^2)(Subscript[η, t]))+I(-Sign[(Subscript[α, yy])](Sqrt[S k \[LeftBracketingBar](Subscript[α, yy])\[RightBracketingBar]/2])- k (Subscript[g, t]) (Subscript[γ, (0)]))

DensityPlot[Re[Exp[-I t]Exp[I 2 π x]],{x,0,2},{t,0,20}]