# The scalar perturbation equations # For the line of sight calculation see lineofsight.txt # For initial conditions see inits.txt # Assume maple 6 or higher (maple 4 is buggy) # Use conformal time. # Includes cdm, baryons, lambda, massive neutrinos and a scalar field # We set the CDM (accel=0) frame from the start restart; no_quint:={psi(t)=0,V(t)=0,V1(t)=0,V2(t)=0,rhopsi(t)=0,ppsi(t)=0}; no_numassive:={pinu(t)=0,qnu(t)=0,rhonu(t)=0,pnu(t)=0}; H_t:=diff(S(t),t)/S(t); dS:=S(t)*H(t); dH:=-1/6*S(t)^2*kappa*(rho(t)+3*p(t)); Friedmann:=H(t)^2=1/3*S(t)^2*kappa*rho(t)-K; #Eliminate phi(t) [Weyl tensor variable] using constraint equation cons1:=k^3*Kfac*phi(t)+1/2*kappa*S(t)^2*k*(rhox(t)+rhopi(t)*Kfac)+3/2*kappa*S(t)^2*H(t)*rhoq(t); solve(cons1,phi(t)); phi_sub:=simplify(%); #Eliminate z [perturbation in expansion rate] in favour of the curvature perturbation eta z_sub:=solve(k^2*eta(t)=kappa*S(t)^2*rhox(t) -2*k*H(t)*z(t),z(t)); #Solve for sigma(t) [shear] using constraint cons2:=2/3*(k/S(t))^2*(z(t)-Kfac*sigma(t))+kappa*rhoq(t); solve(cons2,sigma(t)); sigma_sub:=subs(z(t)=z_sub,%); solve(cons2,rhoq(t)); rhoq_sub:=%; varsubs:={z(t)=z_sub,sigma(t)=sigma_sub,phi(t)=phi_sub}; drag_t:=opac(t)*(4/3*v(t)-qg(t)); photbar_t:=rhog(t)/rhob(t); dgrho_t:=kappa*S(t)^2*rhox(t); drhob:=-3*H(t)*(rhob(t)+pb(t)); drhoc:=-3*H(t)*rhoc(t); drhog:=-4*H(t)*rhog(t); drhor:=-4*H(t)*rhor(t); drhonu:=-3*H(t)*(rhonu(t)+pnu(t)); drhopsi:=-3*H(t)*(rhopsi(t)+ppsi(t)); dclxr:=-4/3*k*z(t)-k*qr(t); dclxg:=-4/3*k*z(t)-k*qg(t); dclxb:=-(1+pb(t)/rhob(t))*k*(z(t)+v(t)) +(pb(t)/rhob(t)-c2(t))*3*H(t)*clxb(t); dclxc:=-k*z(t); dclxnu:=-k*(1+pnu(t)/rhonu(t))*z(t) - k*qnu(t) +3*H(t)*(-clxpnu(t) + clxnu(t)*pnu(t)/rhonu(t)); dclvdot:= - 2*H(t)*diff(clv(t),t) - k*z(t)*diff(psi(t),t) - k^2*clv(t) - clv(t)*S(t)^2*V2(t); dpsidot:=-2*H(t)*diff(psi(t),t)-S(t)^2*V1(t); rhopsi_t:=1/2*diff(psi(t),t)^2/S(t)^2 + V(t); ppsi_t:=1/2*diff(psi(t),t)^2/S(t)^2 - V(t); dz:=-H(t)*z(t)-1/2*kappa*S(t)^2/k*(2*rhog(t)*clxg(t)+2*rhor(t)*clxr(t)+(1+3*c2(t))*rhob(t)*clxb(t)+rhoc(t)*clxc(t)+rhonu(t)*(clxpnu(t)*3+ clxnu(t)) + 1/S(t)^2*(4*diff(psi(t),t)*diff(clv(t),t)-2*clv(t)*S(t)^2*V1(t))); Kfac_sub:= (1-3*K/k^2); dqr:=-2/3*k*pir(t)*Kfac+1/3*k*clxr(t); dqg:=-2/3*k*pig(t)*Kfac+1/3*k*clxg(t)+drag(t); dv:=-(1-3*c2(t))*H(t)*v(t)+1/(1+pb(t)/rhob(t))*(c2(t)*k*clxb(t)-photbar(t)*drag(t)); dqnu:=-H(t)*(1-3*pnu(t)/rhonu(t))*qnu(t) -k/3*(2*Kfac*pinu(t) - 3*clxpnu(t)); dsigma:=simplify(-H(t)*sigma(t)+k*phi(t)-1/2*kappa*S(t)^2/k*(rhopi(t))); simplify(subs(rhoq(t)=rhoq_sub,dsigma)); dsigma:=%; rho_t:=rhob(t)+rhoc(t)+rhor(t)+rhog(t)+rhonu(t)+rhopsi(t)+rhov(t); p_t:=1/3*(rhor(t)+rhog(t))+pb(t)+pnu(t)+ppsi(t)-rhov(t); dpb:=c2(t)*drhob; dpig:=-opac(t)*(pig(t)-polter(t)) - 3/5*k*Kf[2]*J_3(t) + 2/5*k*qg(t) + 8/15*k*sigma(t); dJ_3:=k*(3/7*pig(t)-4/7*Kf[3]*J_4(t))-opac(t)*J_3(t); dpir:=- 3/5*k*Kf[2]*G_3(t) + 2/5*k*qr(t) + 8/15*k*sigma(t); dG_3:=k*(3/7*pir(t)-4/7*Kf[3]*G_4(t)); dG_4:=k*(4/9*G_3(t)-5/9*Kf[4]*G_5(t)); G_eq:=proc(l) local Gl,Eq; Gl:=cat('G_',l)(t); Eq:=diff(Gl,t)+k/(2*l+1)*( (l+1)*Kf[l]*cat('G_',l+1)(t) - l*cat('G_',l-1)(t)); if (l < 3) then print("ERROR: only for l>2"); fi; simplify(subs(G_2(t)=pir(t),Eq)); end; rhopi_t:=rhog(t)*pig(t)+rhor(t)*pir(t) + rhonu(t)*pinu(t); rhox_t:=rhog(t)*clxg(t)+rhor(t)*clxr(t)+rhob(t)*clxb(t)+rhoc(t)*clxc(t) +rhonu(t)*clxnu(t) + 1/S(t)^2*(diff(psi(t),t)*diff(clv(t),t)+clv(t)*S(t)^2*V1(t)); rhoq_t:=rhog(t)*qg(t)+rhor(t)*qr(t)+(rhob(t)+pb(t))*v(t) + rhonu(t)*qnu(t) +k*diff(psi(t),t)*clv(t)/S(t)^2; #pressure perturbation rhoxp_t:=rhonu(t)*clxpnu(t) + 1/3*(rhog(t)*clxg(t)+rhor(t)*clxr(t)) + 1/S(t)^2*(diff(psi(t),t)*diff(clv(t),t)-clv(t)*S(t)^2*V1(t))+c2(t)*rhob(t)*clxb(t); subtots:={rhopi(t)=rhopi_t,rhoq(t)=rhoq_t,rhox(t)=rhox_t,rho(t)=rho_t,p(t)=p_t,rhoxp(t)=rhoxp_t}; subtotsquint := {rhopsi(t)=rhopsi_t, ppsi(t)=ppsi_t}; deta:=simplify(-1/k*(2*K*z(t) + kappa*S(t)^2*rhoq(t))); drhoq:=k*rhoxp(t)-4*H(t)*rhoq(t)-2/3*k*Kfac*rhopi(t); drhox:=-k*rhoq(t) - 3*H(t)*(rhox(t) + rhoxp(t)) - k*z(t)*(rho(t)+p(t)); polter_t:=2/15*(3*pig(t)/4 + 9*E2(t)/2); dE2:=-opac(t)*(E2(t) - polter(t)) - k/3*Kf[2]*E3(t); E_eq:=proc(l) local El,Eq; El:=cat('E',l)(t); Eq:=diff(El,t)+k/(2*l+1)*( (l+3)*(l-1)/(l+1)*Kf[l]*cat('E',l+1)(t) - l*cat('E',l-1)(t))+opac(t)*El; if (l < 2) then print("ERROR: only for l>=2"); fi; if (l=2) then Eq:=Eq-polter(t)*opac(t); fi; simplify(subs(E1(t)=0,Eq)); end; sublist:={diff(S(t),t)=dS,diff(rhob(t),t)=drhob,diff(rhoc(t),t)=drhoc,diff(rhor(t),t)=drhor,diff(rhog(t),t)=drhog,diff(clxr(t),t)=dclxr,diff(clxg(t),t)=dclxg,diff(clxb(t),t)=dclxb,diff(clxc(t),t)=dclxc,diff(z(t),t)=dz,diff(qr(t),t)=dqr,diff(qg(t),t)=dqg,diff(v(t),t)=dv,diff(sigma(t),t)=dsigma,diff(H(t),t)=dH,diff(exptau(t),t)=g(t),diff(pig(t),t)=dpig,diff(rhonu(t),t)=drhonu,diff(clxnu(t),t)=dclxnu,diff(eta(t),t)=deta,diff(diff(clv(t),t),t)=dclvdot,diff(diff(psi(t),t),t)=dpsidot,diff(V1(t),t)=diff(psi(t),t)*V2(t)}; subtotderivs:={diff(rhox(t),t)=drhox,diff(rhoq(t),t)=drhoq,diff(S(t),t)=dS,diff(H(t),t)=dH,diff(eta(t),t)=deta}; #End of main definitions and equations #########################################