AssumpCS={CS->L*(L+hd)}
AssumpCS2={CS2->CS-(hd+1)}
AssumpCS3={CS3->CS2-(hd+1)}
AssumpMs={Ms[r]->CS CS2 (hd^2)-2 hd (1+hd) CS2 (r^2) (ω^2)+(((1+hd))^2) (r^4) (ω^4)}
α=L+hd/2
Print["Bessel functions to next-to-next-to-leading-order"]
besselj[x_]=Simplify[1-((x^2)/(2*(2α+2)))+((x^4)/(4*2*(2α+2)(2α+4)))]
AssumpBesselEq={j''[ω r]->FullSimplify[-(((2 α+1)/(ω r))j'[ω r]+j[ω r])]}
AssumpBesselSub={j[ω r]->besselj[ω r],j'[ω r]->besselj'[ω r]}

AssumpTS={Ts[r]->Taa[r]-Trr[r]}
TS=(1/(Nld*Omd (hd+1)(r^2)))∫ (Subscript[N, L])[Om]Ts[r]\[DifferentialD]Om
Print["pref:"]
pref=((Nld*Omd*(hd+1)*(r^(L+hd+3)))/(2 CS))
AssumpXY={X[r]->(((ω^2)(r^(hd+7)))/(Ms[r]))TS,Y[r]->(((hd+1)CS2 (r^2) (CS hd -(hd+1)(ω^2)(r^2)))/(Ms[r])) TS}
Print["J:"]
J=Y[r]+((((hd+1))^3)/(r^(hd+2)))(∂_r)X[r]
Print["𝒯:"]
𝒯=pref*((hd/r)(∂_r)J+(ω^2)J-((hd CS3)/((hd+1)(r^2)))J);
Print["original integrand 𝒯 j[ω r]:"]
integrand=FullSimplify[𝒯 j[ω r]]
fullderivative1=(∂_r)(pref * (hd/r) * j[ω r]*J);
fullderivative2=(∂_r)(-pref*(r^(-(hd+4))) ((hd (CS3+(1+hd) (2+hd+L))-(1+hd) (r^2) (ω^2)) j[r ω]+hd (1+hd) r ω j'[r ω]) (((1+hd))^2) X[r]);
Print["After integrations by parts:"]
integrand2=FullSimplify[integrand-fullderivative1-fullderivative2]
Print["After Bessel Equation:"]
integrand3=FullSimplify[integrand2/.AssumpBesselEq]
Print["After plugging X,Y:"]
integrand4=FullSimplify[integrand3/.AssumpXY]
commonfactor=(1/(2 CS ))(1+hd) (r^(1+hd+L)) (∫(Subscript[N, L])[Om]\[DifferentialD]Om) Ts[r];
Print["without common prefactor:"]
expr=FullSimplify[integrand4/commonfactor]

Print["work on Ts[r] j'[r ω]:"]
Factorjderv=FullSimplify[FullSimplify[expr/.{j[ω r]->0}]/.AssumpCS3]
Factorjderv2=FullSimplify[FullSimplify[Factorjderv/.AssumpMs]/.AssumpCS2]
Print["work on Ts[r] j[r ω]:"]
Factorj=FullSimplify[FullSimplify[expr/.{j'[r ω]->0}]/.AssumpCS3]
Factorj2=FullSimplify[FullSimplify[FullSimplify[Factorj/.AssumpMs]/.AssumpCS2]/.AssumpCS]*((CS/(L(L+hd))))
Print["Final integrand:"]
integrand5=FullSimplify[commonfactor*(Factorjderv2+Factorj2)]
Print["test:"]
FullSimplify[FullSimplify[FullSimplify[FullSimplify[integrand4-integrand5/.AssumpCS3]/.AssumpMs]/.AssumpCS2]/.AssumpCS]
integrand6=FullSimplify[integrand5/.{(Subscript[N, L])[Om]->(((Subscript[x, L])[Om,r])/(r^L))}]
Q=FullSimplify[∫(integrand6) \[DifferentialD]r/.AssumpTS]