(**************************************************************************
*                                                                         *
*                   Blue Axion Isocurvature Spectrum                           *
*                A mathematica package for evaluating                     *
*        isocurvature power spectrum for resonant underdamped             *
*                       fields with (SubPlus[c]) > 9/4                              *
*                                                                         *
*                             Version 1.1                                 *
*                         Dated: 30-Sept-2021                             *
*                                                                         *
*                                                                         *
*                     By Sai Chaitanya Tadepalli                          *
*                          and Daniel J Chung                             *
*                                                                         *
*          Affiliated to University of Wisconsin-Madison, USA             *
*                                                                         *
*                                                                         *
*                    Contact: stadepalli@wisc.edu                         *
*                                                                         *
*              Copyright (c) 2021  Chung Research Group                   *
*                         All rights reserved                             *
*                                                                         *
*               Supported in part by Ray Mac Donald Fund                  *
*                                                                         *
*              Refer to _______ for details about the code,              *
*        the underlying numerical model and parameter definitions.        *
*                                                                         *
*                                                                         *
* This software package is provided without guarantee of support or       *
* maintenance. The copyright holder makes no express or implied warranty  *
* of any kind with respect to this software, including implied warranties *
* of merchantability or fitness for a particular purpose, and is not      *
* liable for any damages resulting in any way from its use.               *
*                                                                         *
*                                                                         *
* Everyone is granted permission to copy, modify and redistribute this    *
* software package, provided that:                                        *
*  1. All copies contain this notice in the main program file and in the  *
*     supporting documentation.                                           *
*  2. All modified copies carry a prominent notice stating who made the   *
*     last modification and the date of such modification.                *
*  3. No charge is made for this software or works derived from it.       *
*                                                                         *
*                                                                         *
*                                                                         *
*       Input Lagrangian Paramaters: [ (SubPlus[ϕ])(0) , (Subscript[ϵ, 0]) , F , (SubMinus[c]) , (SubPlus[c]) ]         *
*                                                                         *
*  To generate isocurvature power spectrum for resonant underdamped cases *
*                                                                         *
*  1. Load all the functions                                              *
*  2. Plot spectrum using the following command                           *
*     "ListLogLogPlot[Isospectrum[e0,phip0,cm,F,cp]]"                    *
*     using the user-specific lagrangian parameters.                      *
*                                                                         *
*  For instance consider:                                                 *
*       (SubPlus[ϕ])(0)=0.1 (Subscript[M, plank])/H, (Subscript[ϵ, 0])=0, F≡(Subscript[f, PQ])/H=20.2, (SubPlus[c])=2.34, (SubMinus[c])=0.5           *
*                                                                         *
*  ListLogLogPlot[Isospectrum[0,3.3*^8,0.5,20.2,2.34],FrameLabel->{"K/a(0)H","(Subscript[(Δ^2), s])(K)/(Subscript[(ω^2), a])"},Frame->True]          *
*                             *
**************************************************************************)

(*  Clear["Global`*"]       *)

theta=0.02; (* default theta value *)

AlphaEstimation[e0_,phip0_,cm_,F_,cp_]:=(w:=Sqrt[cp-9/4];
varphi:=ArcTan[(3/2+e0)/w];
amax=7;
j :=Ceiling[-1/2+1/Pi (-varphi+2 w/3Log[(phip0/F^2 (Sqrt[cp]) /amax)])];
Tz:=1/w ((2j+1) Pi/2+varphi);alpha:=w phip0/F^2 Sec[varphi] Exp[-3/2 Tz];
alpha)

(*         STARTING FROM PERTURBED SOLN, SOLVE FOR (Subscript[T, c]) USING CUBIC POLYNOMIAL EXPN               *)

CubicPolyfields[e0_,phip0_,cm_,F_,cp_]:=(w:=Sqrt[cp-9/4];
varphi:=ArcTan[(3/2+e0)/w];
amax=7;
j :=Ceiling[1/2+1/Pi (-varphi+2 w/3Log[(phip0/F^2 (Sqrt[cp]) /amax)])];
Tz:=1/w ((2j-1) Pi/2+varphi);alpha:=w phip0/F^2 Sec[varphi] Exp[-3/2 Tz];If[alpha>1.13,Print["For (SubPlus[c]) = ",cp,", α = ",alpha," is outside the region of applicability. Exiting: Subroutine 'CubicPolyfields'"]; Abort[];,If[alpha<0.25,Print["For (SubPlus[c]) = ",cp,", α = ",alpha," is outside the region of applicability. Exiting: Subroutine 'CubicPolyfields'"]; Abort[];,]];
a:=alpha;
Ke:=-a F^2;
n:=If[a>1,10,14];
e:=(n/2)^0.25 Sqrt[-2/Ke]/2;
zI:=e;
Ts:=Tz-2 e;
xs:=Ts-Tz;
phipTs:=Ke (xs-3/2 xs^2+(((9 w)/8-w^3/6) xs^3)/w);
dphipTs:=Ke (1-3 xs+(3 ((9 w)/8-w^3/6) xs^2)/w);
phimTs:=F^2 phipTs/(cm+phipTs^2);
dphimTs:=-dphipTs phimTs/phipTs;
(*----------Evaluate Coefficients-----------*)
p0=phipTs;
q0=phimTs;
p1=dphipTs;
q1=dphimTs;
p2=-3/2 p1-cp p0;
q2=-3/2 q1-cm q0;
X=(p0+p1 zI+p2 zI^2);
Y=(q0+q1 zI+q2 zI^2);
p3=(-(((-F^2+2 X Y) (3 q1+(2+6 zI) q2-F^2 X+(cm+X^2) Y))/((F^2-2 X Y)^2-(9+cm+6/zI^2+X^2) (9+cp+6/zI^2+Y^2)))-(-(9+cm+6/zI^2+X^2) (3 p1+(2+6 zI) p2-F^2 Y+X (cp+Y^2)))/((F^2-2 X Y)^2-(9+cm+6/zI^2+X^2) (9+cp+6/zI^2+Y^2)))/zI^3;
q3=(-(((-F^2+2 X Y) (3 p1+(2+6 zI) p2-F^2 Y+(cp+Y^2) X))/((F^2-2 X Y)^2-(9+cm+6/zI^2+X^2) (9+cp+6/zI^2+Y^2)))-(-(9+cp+6/zI^2+Y^2) (3 q1+(2+6 zI) q2-F^2 X+Y (cm+X^2)))/((F^2-2 X Y)^2-(9+cm+6/zI^2+X^2) (9+cp+6/zI^2+Y^2)))/zI^3;
phip[z_]:=p0+p1 z+p2 z^2+p3 z^3;
phim[z_]:=q0+q1 z+q2 z^2+q3 z^3;
a1:=p0-q0; b1:=p1-q1;c1:=p2-q2;d1:=p3-q3;
zc=Solve[a1+b1 z+ c1 z^2 + d1 z^3==0,z];
Tc=Ts+z/.zc[[Which[(Im[z/.zc[[1]]]==0&& Re[z/.zc[[1]]]>0),1,(Im[z/.zc[[2]]]==0&& Re[z/.zc[[2]]]>0),2,(Im[z/.zc[[3]]]==0&& Re[z/.zc[[3]]]>0),3]]];
g[z_]:=(cp+phim[z+(Tc-Ts)]^2)-(cm+phip[z+(Tc-Ts)]^2);
e1[z_]:=(g[z]-Sqrt[g[z]^2+4F^4])/2/F^2;
dg[z_]:=g'[z];
factor[z_]:=e1[z]^2/((e1[z]^2+1)^2(g[z]^2+4F^4));
e1primesq[z_]:=dg[z]^2factor[z];
(*------------------Now expand around Tc------z=0, T=Tc-------------*)
u0=phip[Tc-Ts];
v0=phim[Tc-Ts];
u1=phip'[Tc-Ts];
v1=phim'[Tc-Ts];
u2=(1/2) (-cp u0-3 u1+(F^2) v0-u0 (v0^2));
v2=(1/2) (-cm v0-3 v1+(F^2) u0-v0 (u0^2));
xI=e; 
X1=(u0+u1 xI+u2 xI^2);
Y1=(v0+v1 xI+v2 xI^2);
u3=(-(((-F^2+2 X1 Y1) (3 v1+(2+6 xI) v2-F^2 X1+(cm+X1^2) Y1))/((F^2-2 X1 Y1)^2-(9+cm+6/xI^2+X1^2) (9+cp+6/xI^2+Y1^2)))-(-(9+cm+6/xI^2+X1^2) (3 u1+(2+6 xI) u2-F^2 Y1+X1 (cp+Y1^2)))/((F^2-2 X1 Y1)^2-(9+cm+6/xI^2+X1^2) (9+cp+6/xI^2+Y1^2)))/xI^3;
v3=(-(((-F^2+2 X1 Y1) (3 u1+(2+6 xI) u2-F^2 Y1+(cp+Y1^2) X1))/((F^2-2 X1 Y1)^2-(9+cm+6/xI^2+X1^2) (9+cp+6/xI^2+Y1^2)))-(-(9+cp+6/xI^2+Y1^2) (3 v1+(2+6 xI) v2-F^2 X1+Y1 (cm+X1^2)))/((F^2-2 X1 Y1)^2-(9+cm+6/xI^2+X1^2) (9+cp+6/xI^2+Y1^2)))/xI^3;
xi[z_]:=(u0+u1 z+u2 z^2+u3 z^3)(v0+v1 z+v2 z^2+v3 z^3)-F^2;
A=FindMinimum[{xi[z],-xI<z<xI},{z}];
m11:=cp+(v0+v1 z+v2 z^2+v3 z^3)^2;
m22:=cm+(u0+u1 z+u2 z^2+u3 z^3)^2;
m1sq[z_]:=1/2(m11+m22-Sqrt[(m11-m22)^2+4F^4]);
m2sq[z_]:=1/2(m11+m22+Sqrt[(m11-m22)^2+4F^4]);
E1=FindMaximum[{e1primesq[z],-2/F<z<0},{z}];
V1=MinValue[{-e1primesq[z]+m1sq[z],-2/F<z<2/F},{z}];
x1=FindRoot[-e1primesq[z]+m1sq[z]==0,{z,-2/F}];
x2=FindRoot[-e1primesq[z]+m1sq[z]==0,{z,2/F}];
wid=NIntegrate[{-e1primesq[z]+m1sq[z]},{z,x1[[1,2]],x2[[1,2]]}]/V1;
z2=(3.11-1.05alpha)/F/2;
Pm=(v0+v1 z+v2 z^2+v3 z^3)/.z->z2;
Pp=(u0+u1 z+u2 z^2+u3 z^3)/.z->z2;
DPm=(v1 +2v2 z+3v3 z^2)/.z->z2;
DPp=(u1 +2u2 z+3u3 z^2)/.z->z2;
phimmin=Sqrt[F^2Sqrt[cp/cm]-cp];
b=DPm-xi[z2]Pp/Pm;
Avgomega=F(2.05+((0.1327 +0.0454 F)/(1+(E^(7.86 (-0.744+a+0.0008 F))))));
nn=Sqrt[1-4cm/9];
m=-Sqrt[2]A[[1]]/3/Avgomega;
c1=-1/(((2  Sin[(1/3) (Sqrt[9-4 cm]) π])/(m π  )))(Pm D[BesselJ[-nn,x],x]+(Pm+2 b/3)BesselJ[-nn,m]/m)/.x->m;
c2=1/(((2  Sin[(1/3) (Sqrt[9-4 cm]) π])/(m π  )))(Pm D[BesselJ[nn,x],x]+(Pm+2 b/3)BesselJ[nn,m]/m)/.x->m;
TL=Tc+z2-2/(3(1-nn))Log[phimmin 2^(2-nn)m^nn Gamma[1-nn]/3/c2];
;
l={{a,Tz,Tc,-E1[[1]],V1,A[[1]],TL,Pm,b,Avgomega,m,g'[0],c1,c2}};
l)

ModelParameters[e0_,phip0_,cm_,F_,cp_]:=(SolnPart1=CubicPolyfields[e0,phip0,cm,F,cp];
alpha=SolnPart1[[1,1]];
Tc=SolnPart1[[1,3]];
E1=SolnPart1[[1,4]];
V1=SolnPart1[[1,5]];
A=-SolnPart1[[1,6]];
TL=SolnPart1[[1,7]];
Pm=SolnPart1[[1,8]];
DPm=SolnPart1[[1,9]];
omega=SolnPart1[[1,10]];
dg0=SolnPart1[[1,12]];
c11=SolnPart1[[1,13]];
c12=SolnPart1[[1,14]];
;
V1width=(3.11-1.05alpha)/F;
T0=0;
T1=Tc-V1width/2;
T2=Tc+V1width/2;
;
nn=Sqrt[1-4cm/9];
m=Sqrt[2]A/3/omega;
phims[x_]:=Exp[-3/2x](c11 BesselJ[nn,m Exp[-3/2x]]+c12 BesselJ[-nn,m Exp[-3/2x]]);
g[x_]:=phims[x]^2-F^4/phims[x]^2;
e1[x_]:=(g[x]-Sqrt[g[x]^2+4F^4])/2/F^2;
dg[x_]:=g'[x](1+(A Exp[-3/2x])^2/Sqrt[(Sqrt[4 (F^4)+((g[x])^2)])]^4)+(A Exp[-3/2x])^2/Sqrt[(Sqrt[4 (F^4)+((g[x])^2)])]^3 g[x];
e1primesq[x_]:=dg[x]^2(e1[x]^2/(e1[x]^2+1)^2 /(g[x]^2+4F^4));
;
V2=A^2 ((.139- .14/(1.08+Exp[11(alpha-.72)])))/F^2 ;
;
If[m>3.8/Sqrt[2],xroot=FindRoot[phims[x]==F,{x,.5}];
V3 = FindMinimum[{-e1primesq[x],xroot[[1,2]]-3/F<x<xroot[[1,2]]+1/F},{x}];
xpos=x/.V3[[2]][[1]];
(*V3width=NIntegrate[e1primesq[x],{x,xpos-0.5/Sqrt[9/4-V3[[1]]],xpos+0.5/Sqrt[9/4-V3[[1]]]}]/(-V3[[1]]);*)
V3width=0.7/Sqrt[9/4- V3[[1]]];
T3=Tc+xpos-V3width/2;
T4=T3+V3width;
V4=V2 Exp[-3(T4-T2)];,T3=I;V3=Null;V4=Null;T4=Null;];
{{alpha,cp,V1,V2,V3[[1]],V4,T0,T1,T2,T3,T4,TL,DPm}})

Isospectrum[e0_,phip0_,cm_,F_,cp_]:=(r=Sqrt[cp/cm];
Ipearly[K_,T_]:=Exp[-T]Exp[I K Exp[-T]]/2/Sqrt[K];
y0[K_,T_]:={{Ipearly[K,T]},{Derivative[0,1][Ipearly][K,T]}};
Parameters=ModelParameters[e0,phip0,cm,F,cp];
{{alpha,V0,V1,V2,V3,V4,T0,T1,T2,T3,T4,TL,DP}}=Parameters;
;
(*------------------------defining bessel fns for transfer matrices----------------*);
y1[n_,K_,T_]:=Exp[-3/2 T]BesselJ[n,K Exp[-T]];
y2[n_,K_,T_]:=Exp[-3/2 T]BesselY[n,K Exp[-T]];
W[n_,K_,T_]:=Derivative[0,0,1][y1][n,K,T]y2[n,K,T]-Derivative[0,0,1][y2][n,K,T]y1[n,K,T];
S[n_,K_,t0_,t1_]:={{y1[n,K,t1],y2[n,K,t1]},{Derivative[0,0,1][y1][n,K,t1],Derivative[0,0,1][y2][n,K,t1]}}.{{-Derivative[0,0,1][y2][n,K,t0],y2[n,K,t0]},{Derivative[0,0,1][y1][n,K,t0],-y1[n,K,t0]}}/W[n,K,t0];
(*-----------------------------------------------*);
y15[n_,c_,m_,T_]:=Exp[-3/2 T]BesselJ[Sqrt[9-4c]/n,(2/n) Sqrt[m] Exp[-n/2T]];
y25[n_,c_,m_,T_]:=Exp[-3/2 T]BesselY[Sqrt[9-4c]/n,(2/n) Sqrt[m] Exp[-n/2T]];
W5[n_,c_,m_,T_]:=Derivative[0,0,0,1][y15][n,c,m,T]y25[n,c,m,T]-Derivative[0,0,0,1][y25][n,c,m,T]y15[n,c,m,T];
S5[n_,m_,c_,t0_,t1_]:={{y15[n,c,m,t1],y25[n,c,m,t1]},{Derivative[0,0,0,1][y15][n,c,m,t1],Derivative[0,0,0,1][y25][n,c,m,t1]}}.{{-Derivative[0,0,0,1][y25][n,c,m,t0],y25[n,c,m,t0]},{Derivative[0,0,0,1][y15][n,c,m,t0],-y15[n,c,m,t0]}}/W5[n,c,m,t0];
(*-----------------------------------------------*);
Tx[K_]:=Log[V /K^2];
Tv[K_]:=Tx[K]-2Log[2];
Tk[K_]:=Tx[K]+2Log[2];
V:=V2/2 Exp[3 T2];
(*-----------------------------------------------------------------*);
S8[K_,c_,t0_,t1_]:=If[Tv[K]>=t0,If[Tv[K]>t1,S5[3,(3 ((1/2) ((E^(-2 t0))-(E^(-2 t1))) (K^2)+(1/3) ((E^(-3 t0))-(E^(-3 t1))) V))/((E^(-3 t0))-(E^(-3 t1))),c,t0,t1],If[Tk[K]>t1,S5[5/2,(5 ((1/2) ((E^(-2 Tv[K]))-(E^(-2 t1))) (K^2)+(1/3) ((E^(-3 Tv[K]))-(E^(-3 t1))) V))/(2 ((E^(-5 Tv[K]/2))-(E^(-5 t1/2)))) ,c,Tv[K],t1],S5[2,(2 ((1/2) ((E^(-2 Tk[K]))-(E^(-2 t1))) (K^2)+(1/3) ((E^(-3 Tk[K]))-(E^(-3 t1))) V))/((E^(-2 Tk[K]))-(E^(-2 t1))) ,c,Tk[K],t1].S5[5/2,(5 ((1/2) ((E^(-2 Tv[K]))-(E^(-2 Tk[K]))) (K^2)+(1/3) ((E^(-3 Tv[K]))-(E^(-3 Tk[K]))) V))/(2 ((E^(-5 Tv[K]/2))-(E^(-5 Tk[K]/2)))) ,c,Tv[K],Tk[K]]].S5[3,(3 ((1/2) ((E^(-2 t0))-(E^(-2 Tv[K]))) (K^2)+(1/3) ((E^(-3 t0))-(E^(-3 Tv[K]))) V))/((E^(-3 t0))-(E^(-3 Tv[K]))),c,t0,Tv[K]]],If[Tk[K]>=t0,If[Tk[K]>t1,S5[5/2,(5 ((1/2) ((E^(-2 t0))-(E^(-2 t1))) (K^2)+(1/3) ((E^(-3 t0))-(E^(-3 t1))) V))/(2 ((E^(-5 t0/2))-(E^(-5 t1/2)))),c,t0,t1],S5[2,(2 ((1/2) ((E^(-2 Tk[K]))-(E^(-2 t1))) (K^2)+(1/3) ((E^(-3 Tk[K]))-(E^(-3 t1))) V))/((E^(-2 Tk[K]))-(E^(-2 t1))),c,Tk[K],t1].S5[5/2,(5 ((1/2) ((E^(-2 t0))-(E^(-2 Tk[K]))) (K^2)+(1/3) ((E^(-3 t0))-(E^(-3 Tk[K]))) V))/(2 ((E^(-5 t0/2))-(E^(-5 Tk[K]/2)))),c,t0,Tk[K]]],S5[2,(2 ((1/2) ((E^(-2 t0))-(E^(-2 t1))) (K^2)+(1/3) ((E^(-3 t0))-(E^(-3 t1))) V))/((E^(-2 t0))-(E^(-2 t1))),c,t0,t1]].IdentityMatrix[2]];
(*------------------------------------------------------------------------------------------*);
Q:={{1,0},{-Sqrt[V2],1}};
Lam:=3/2-Sqrt[9/4-4 cm cp/(cm+cp)];
mbIntegral1:=1063/3072+106793/393216 cm/cp+ cm(TL-T2)+3/4(243-72cp+81cm/cp+8cm(2cp-27))(F/DPm)^4 ;
MB1=mbIntegral1/(TL-T2);
Clear[n];
mbIntegral2:=NIntegrate[((cm cp (-1+(n^4)) (-cm+cp (n^4)))/(((cm+cp (n^4)))^2))/.n->(1-1/3Exp[-Lam x]),{x,0,Infinity}];
MB2= mbIntegral2/(2/Lam);
(*-----------------------------------------------*);
f1[K_]:=S8[K,0,TL,50].S8[K,MB1,T2,TL].Q.S[Sqrt[9/4-V1],K,T1,T2].S[Sqrt[9/4-V0],K,T0,T1];
f3[K_]:=S8[K,0,T4+2/Lam,40].S8[K,MB2,T4,T4+2/Lam].S8[K,V3,T3,T4].S8[K,cm,T2,T3].Q.S[Sqrt[9/4-V1],K,T1,T2].S[Sqrt[9/4-V0],K,T0,T1];
;
P[K_]:=If[Im[T3]==0,f3[K],f1[K]]; 
(*-----------------------------------------------*);
ysol[K_]:=(P[K].y0[K,T0])[[1,1]];
;
term=(1+r^4)/(1+r^2)^2 / (r+ 1/r);
Cterm=(2 term/ Pi^2 /F^2/theta^2);
;
spectrum[K_]:=Abs[ysol[K]]^2 ( K^3  Cterm);
;
(*           select K interval and width for power spectrum plot             *);
Kmin=1000;
Kmax=2*^5;
dK=1000;
;
;
specdata=Table[{k ,spectrum[k]},{k,Kmin,Kmax,dK}];
(*-----------------------------------------------*);
specdata)

ListLogLogPlot[Isospectrum[0,3.3*^8,0.5,20.2,2.34],FrameLabel->{"K/a(0)H","(Subscript[(Δ^2), s])(K)/(Subscript[(ω^2), a])"},Frame->True,PlotRange->{{1000,2*^5},{.0001,1}},Joined->True,PlotMarkers->{{"*",Small}},PlotStyle->{Blue}]