Quit[]

<<EDCRGTCcode.m

$Assumptions={Ω∈Reals, Ω>0, pv∈Reals, nx∈Integers, nx>=0, ny∈Integers, ny>=0, pv!=0};

1. Background & Linear Perturbations

1.1. Construct background, spin coefficients and GHP operators

coords = {u,v,x,y};
met = DiagonalMatrix[{-H[u,x,y],0,-1,-1}]+{{0,1,0,0},{1,0,0,0},{0,0,0,0},{0,0,0,0}};
RGtensors[met,coords];
lU = {0,1,0,0};
ld=gdd.lU;
nU={1,((H[u,x,y])/2),0,0};
nd=gdd.nU;
mU=(1/(Sqrt[2])){0,0,1,I};
md=gdd.mU;
mcU=Conjugate[mU];
mcd=gdd.mcU;
PsiH = {ψH[u,v,x,y],-4,0};
PsiHc = {ψHc[u,v,x,y],0,-4};
PsiHB = {ψHB[u,v,x,y],-4,0};
PsiHcB = {ψHcB[u,v,x,y],0,-4};
Kerr = FuncRepRules[H[u,x,y], -(Ω^2)((x^2)-(y^2))];
κ =(covD[ld].lU).mU;
τ = (covD[ld].nU).mU;
σ = (covD[ld].mU).mU;
ρ = (covD[ld].mcU).mU;
κp =(covD[nd].nU).mcU;
τp =(covD[nd].lU).mcU;
σp =(covD[nd].mcU).mcU;
ρp =(covD[nd].mU).mcU;
ωd= (1/2)(nU.covD[ld]+mU.covD[mcd]);
þ[f_]  := {covD[f[[1]]].lU -(f[[2]] lU.ωd +f[[3]] lU.Conjugate[ωd])f[[1]],f[[2]]+1,f[[3]]+1};
þp[f_]  := {covD[f[[1]]].nU -(f[[2]] nU.ωd +f[[3]] nU.Conjugate[ωd])f[[1]],f[[2]]-1,f[[3]]-1};
ð[f_]  := {covD[f[[1]]].mU -(f[[2]] mU.ωd +f[[3]] mU.Conjugate[ωd])f[[1]],f[[2]]+1,f[[3]]-1};
ðp[f_]  :={covD[f[[1]]].mcU -(f[[2]] mcU.ωd +f[[3]] mcU.Conjugate[ωd])f[[1]],f[[2]]-1,f[[3]]+1};
Ψ0 =-((( Rdddd.lU).mU).lU).mU;
Ψ1 =-((( Rdddd.lU).nU).lU).mU;
Ψ2 =-((( Rdddd.lU).mU).mcU).nU;
Ψ3 =-((( Rdddd.lU).nU).mcU).nU;
Ψ4 =-((( Rdddd.nU).mcU).nU).mcU //Simplify;

1.2. Define mode solutions

norm[pv_,nx_?EvenQ,ny_?EvenQ] := (1/π)(2^(-nx-ny))Gamma[((1-nx)/2)]Gamma[((1-ny)/2)];
norm[pv_,nx_?OddQ,ny_?EvenQ] := (1/(π nx ))  (2^(-nx-ny))Gamma[1-(nx/2)]Gamma[((1-ny)/2)];
norm[pv_,nx_?EvenQ,ny_?OddQ] :=(1/(π  ny))  (2^(-nx-ny))Gamma[((1-nx)/2)]Gamma[1-(ny/2)];
norm[pv_,nx_?OddQ,ny_?OddQ] :=(1/(π  nx ny  )) (2^(-nx-ny))Gamma[1-(nx/2)]Gamma[1-(ny/2)];
onShell[pv_,nx_,ny_] := {pu->Sign[pv]Ω (nx+(1/2))+Ω I( ny + (1/2)), puC->Sign[pv]Ω (nx+(1/2))-I Ω( ny + (1/2))};
onShell$mode[pv_,nx_,ny_] := norm[pv,nx,ny]Exp[I pu u]Exp[I pv v]Exp[-Abs[pv] Ω ((x^2)/2)]Exp[-I pv Ω ((y^2)/2)]HermiteH[nx, Sqrt[Abs[pv] Ω]x]HermiteH[ny, Sqrt[I pv Ω]y]/.onShell[pv,nx,ny];
onShell$modeC[pv_,nx_,ny_] := norm[pv,nx,ny]Exp[-I puC u]Exp[-I pv v]Exp[-Abs[pv] Ω ((x^2)/2)]Exp[I pv Ω ((y^2)/2)]HermiteH[nx, Sqrt[Abs[pv] Ω]x]HermiteH[ny, Sqrt[-I pv Ω]y]/.onShell[pv,nx,ny];
mode[pu_,pv_,nx_,ny_]:=norm[pv,nx,ny]Exp[I pu u]Exp[I pv v]Exp[-Abs[pv] Ω ((x^2)/2)]Exp[-I pv Ω ((y^2)/2)]HermiteH[nx, Sqrt[Abs[pv] Ω]x]HermiteH[ny, Sqrt[I pv Ω]y];
modeC[pu_,pv_,nx_,ny_]:= norm[pv,nx,ny]Exp[-I Conjugate[pu] u]Exp[-I pv v]Exp[-Abs[pv] Ω ((x^2)/2)]Exp[I pv Ω ((y^2)/2)]HermiteH[nx, Sqrt[Abs[pv] Ω]x]HermiteH[ny, Sqrt[-I pv Ω]y];
Hertz[pv_,nx_,ny_] := ((-2)/(pv^4))onShell$modeC[pv,nx,ny];
HertzC[pv_,nx_,ny_] := ((-2)/(pv^4))onShell$mode[pv,nx,ny];
eigvalue[pu_,pv_,nx_,ny_] := I pv (2 I pu+Ω+2 ny Ω)+(Ω+2 nx Ω) Abs[pv];

Check normalization

((((E^(I pu u+I pv v))))^(-1))mode[pu,pv,0,0]/.x->0/.y->0
((((E^(I pu u+I pv v))))^(-1))D[mode[pu,pv,3,2],x]/.x->0/.y->0
((((E^(I pu u+I pv v))))^(-1))D[mode[pu,pv,2,5],y]/.x->0/.y->0
((((E^(I pu u+I pv v))))^(-1))D[mode[pu,pv,3,5],x,y]/.x->0/.y->0

1.3. Utilities related to Hermite polynomials

MultRule[arg1_,arg2_,dummy_:r] := HermiteH[a_,arg1] :> ((((arg1/arg2)))^(a-2dummy))(((((((arg1/arg2)))^2) -1))^(dummy ))Binomial[a,2dummy] ((Factorial[2dummy])/(Factorial[dummy]))HermiteH[dummy, arg2];
xH = {x HermiteH[a_,x b_]:> (1/(2 b))(HermiteH[a+1, x b]+2a HermiteH[a-1, x b]),x ((HermiteH[a_,x b_])^2):> ((HermiteH[a, x b])/(2 b))(HermiteH[a+1, x b]+2a HermiteH[a-1, x b]), 
(x^2)HermiteH[a_, x b_] :> ((4 (-1+a) a HermiteH[-2+a,b x]+(2+4 a) HermiteH[a,b x]+HermiteH[2+a,b x])/(4 (b^2))), (x^2)((HermiteH[a_, x b_])^2):> (1/(4 (b^2)))(4 (-1+a) a HermiteH[-2+a,b x] HermiteH[a,b x]+2 a ((HermiteH[a,b x])^2)+HermiteH[1+a,b x] (2 a HermiteH[-1+a,b x]+HermiteH[1+a,b x])),
y HermiteH[a_,y b_]:> (1/(2 b))(HermiteH[a+1, y b]+2a HermiteH[a-1, y b]),y ((HermiteH[a_,y b_])^2):> ((HermiteH[a, y b])/(2 b))(HermiteH[a+1, y b]+2a HermiteH[a-1, y b]) ,
(y^2)HermiteH[a_, y b_] :> ((4 (-1+a) a HermiteH[-2+a,b y]+(2+4 a) HermiteH[a,b y]+HermiteH[2+a,b y])/(4 (b^2))), (y^2)((HermiteH[a_, y b_])^2):> (1/(4 (b^2)))(4 (-1+a) a HermiteH[-2+a,b y] HermiteH[a,b y]+2 a ((HermiteH[a,b y])^2)+HermiteH[1+a,b y] (2 a HermiteH[-1+a,b y]+HermiteH[1+a,b y]))};
SumRule[arg_,dummy_:𝒿] := {HermiteH[a_,arg]HermiteH[b_,arg] :> Binomial[a,dummy]Binomial[b,dummy](2^dummy)Factorial[dummy]HermiteH[a+b-2dummy,arg],
((HermiteH[a_,arg])^2):>((Binomial[a,dummy])^2)(2^dummy)Factorial[dummy]HermiteH[2a-2dummy,arg] };

1.4. Check that modes are indeed solutions

We first define the wave operator, and check its expression in coordinates.

box0[f_] :=2(þp[þ[PsiH]]-ðp[ð[PsiH]])[[1]] //.Kerr//. FuncRepRules[ψH[u,v,x,y], f];
box[f_] := - D[f,{y,2}]-D[f,{x,2}]- ((x^2)-(y^2)) (Ω^2) D[f,{v,2}]+2D[f,u,v];
box0[f[u,v,x,y]]==box[f[u,v,x,y]]//Simplify

Check that off-shell modes are eigenfunctions of box, and on-shell modes have eigenvalue zero.

box[mode[pu,pv,nx,ny]]== eigvalue[pu,pv,nx,ny]*mode[pu,pv,nx,ny]//.xH //Simplify
eigvalue[pu,pv,nx,ny]==0/.onShell[pv,nx,ny]//Simplify
(box[onShell$mode[pv,nx,ny]]==0//.xH)//Simplify

2. Computation of Second Order Solutions

2.1. Define source terms

Source[lin_]:= (2*(1/2) (-6 (ðp@Nest[þ,PsiHc,3])[[1]]^2-8(ðp@þ@þ@PsiHc * ðp@þ@þ@þ@þ@PsiHc)[[1]]-2(ðp@þ@PsiHc * ðp@þ@þ@þ@þ@þ@PsiHc)[[1]] + (ðp@ðp@PsiHc * Nest[þ,PsiHc,6])[[1]] + 4(þ@þ@þ@þ@þ@PsiHc * ðp@ðp@þ@PsiHc)[[1]] + 6(þ@þ@þ@þ@PsiHc * ðp@ðp@þ@þ@PsiHc)[[1]] + 4(þ@þ@þ@PsiHc * ðp@ðp@þ@þ@þ@PsiHc)[[1]] + (þ@þ@PsiHc * ðp@ðp@þ@þ@þ@þ@PsiHc)[[1]]
)//.Kerr)//. FuncRepRules[ψHc[u,v,x,y], lin,6];
SourceMin[lin_,linc_] :=(2* (1/2)(
( þ@þ@þ@þ@þ@þ@PsiHc * ð@ð@PsiH)[[1]] + (þ@þ@PsiH * ð@ð@þ@þ@þ@þ@PsiHc)[[1]] -2  (ð@þ@PsiH * ð@þ@þ@þ@þ@þ@PsiHc)[[1]]
)//.Kerr)//.FuncRepRules[ψH[u,v,x,y], lin,6]//. FuncRepRules[ψHc[u,v,x,y], linc,6];

2.2. Routines to compute source term projections

First, extract the exponential factor

Unprotect[Conjugate];
Conjugate[Sign[a_]]:=Sign[a];
Conjugate[Ω]:=Ω;
Protect[Conjugate];

CleanExp[expr_, PU_,PV_,Δ_] := Module[
{prefactor,clean0,clean},
prefactor = Exp[I PU u] Exp[I PV v]Exp[-Abs[PV] Ω (x^2)/2]Exp[-I PV Ω (y^2)/2]Exp[-Δ Ω (x^2)/2];
clean0 = ((expr/prefactor) //ExpandAll);
clean = Total@Select[List@@clean0, FreeQ[#, Exp[_.*u+_.]]&];
Return[clean]
];

(*Coefficient of (x^k) in an expansion in terms of Hermite polynomials*)
coef[k_,pv_,n_]:=((Sqrt[Abs[pv]Ω])/((2^(n+1)) Factorial[n]Sqrt[π])) Sum[CoefficientList[HermiteH[n, Sqrt[Abs[pv]Ω]x],x][[j+1]] *( 1+(((-1))^(j+k)) )* (((Abs[pv]Ω))^(-((1+j+k)/2))) Gamma[((1+j+k)/2)],{j,0,n}];
coefY[k_,pv_,n_]:=((Sqrt[I pv Ω])/((2^(n+1)) Factorial[n]Sqrt[π])) Sum[CoefficientList[HermiteH[n, Sqrt[I pv Ω]x],x][[j+1]] *( 1+(((-1))^(j+k)) )* (((I pv Ω))^(-((1+j+k)/2))) Gamma[((1+j+k)/2)],{j,0,n}];
(*Coefficient of (x^k) Exp[-Δ (x^2)] in an expansion in terms of Hermite polynomials*)
coefG[k_,Δ_,pv_,n_]:=((Sqrt[Abs[pv]Ω])/((2^(n+1)) Factorial[n]Sqrt[π])) Sum[CoefficientList[HermiteH[n, Sqrt[Abs[pv]Ω]x],x][[j+1]] *( 1+(((-1))^(j+k)) )* (((Δ+Abs[pv]Ω))^(-((1+j+k)/2))) Gamma[((1+j+k)/2)],{j,0,n}];

RatioP[{pv1_,nx1_,ny1_},{pv2_,nx2_,ny2_},{NX_,NY_}] := Module[
{pu1,pu2,PU,PV,Δ,src,srcClean, nmaxX,nmaxY, polynomial,scoef,coefList},
pu1 = pu /. onShell[pv1, nx1, ny1];
pu2 = pu /. onShell[pv2, nx2, ny2];
PU = pu1+pu2;
PV=pv1+pv2;
Δ = Abs[pv1]+Abs[pv2]-Abs[PV];
src =Simplify[
(1/2)(Source[HertzC[pv1,nx1,ny1]+HertzC[pv2,nx2,ny2]]-Source[HertzC[pv1,nx1,ny1]]-Source[HertzC[pv2,nx2,ny2]])
];
srcClean=CleanExp[src, PU,PV,Δ];
nmaxX = Exponent[srcClean, x];
nmaxY=Exponent[srcClean,y];
coefList = CoefficientList[srcClean,{x,y}]//Simplify;
scoef = Sum[
coefList[[i+1,j+1]] coefG[i, Δ (Ω/2), PV, NX] coefY[j, PV, NY]
,
{i,0,nmaxX},{j,0,nmaxY}
];
Return[(scoef/(eigvalue[PU,PV,NX,NY]norm[PV,NX,NY]))//Simplify];
]

RatioM[{pv1_,nx1_,ny1_},{pv2_,nx2_,ny2_},{NX_,NY_}] := Module[
{pu1,pu2,PU,PV,Δ,src,srcClean, nmaxX,nmaxY, coefList,scoef},
pu1 = pu /. onShell[pv1, nx1, ny1];
pu2 = Conjugate[pu /. onShell[pv2, nx2, ny2]] ;
PU = pu1-pu2;
PV=pv1-pv2;
Δ = Abs[pv1]+Abs[pv2]-Abs[PV];
src = Simplify[
SourceMin[Hertz[pv2,nx2,ny2],HertzC[pv1,nx1,ny1]]
];
srcClean=CleanExp[src, PU,PV,Δ];
nmaxX = Exponent[srcClean, x];
nmaxY=Exponent[srcClean,y];
coefList = CoefficientList[srcClean,{x,y}]//Simplify;
scoef = Sum[
coefList[[i+1,j+1]] coefG[i, Δ (Ω/2), PV, NX] coefY[j, PV, NY]
,
{i,0,nmaxX},{j,0,nmaxY}
];
Return[(scoef/(eigvalue[PU,PV,NX,NY]norm[PV,NX,NY]))//Simplify];
]

2.3. Self-tests

First, check the projection integrals

((Sqrt[Abs[pv]Ω])/((2^n) Factorial[n]Sqrt[π]))Integrate[ (x^2)HermiteH[0, Sqrt[Abs[pv]Ω]x] Exp[-Abs[pv]Ω (x^2)],{x,-∞,∞}] /. n->0;
coef[2,pv,0]==%

((Sqrt[Abs[pv]Ω])/((2^n) Factorial[n]Sqrt[π]))Integrate[ (x^2)HermiteH[0, Sqrt[Abs[pv]Ω]x] Exp[-Abs[pv]Ω (x^2)]Exp[-Δ (x^2)]/. n->0 /. Δ->3,{x,-∞,∞}] /.n->0;
coefG[2,3,pv,0]==%

Tests the simplest case for the Plus source

r00=Simplify[RatioP[{l,0,0},{l,0,0},{0,0}],{l>0}];
PU = 2pu/.onShell[l,0,0];
PV = 2l;
ψ2 =  r00 * mode[PU,PV,0,0];
lhs=box[ψ2];
rhs=Source[HertzC[l,0,0]];
Simplify[lhs==rhs//.xH,{l>0}]

Now for the Minus source. Notice that the simplest case in this case is different

r00=Simplify[RatioM[{l,0,0},{-l,0,0},{0,0}],{l>0}];
PU = (pu/.onShell[l,0,0])-Conjugate[pu/.onShell[-l,0,0]];
PV = 2l;
ψ2 =  r00 * mode[PU,PV,0,0];
lhs=box[ψ2];
rhs=SourceMin[Hertz[-l,0,0],HertzC[l,0,0]];
Simplify[lhs==rhs//.xH,{l>0}]

3. Compute Ratios

3.1. Coupling of Fundamental Modes

First pa > pb > 0

Simplify[RatioP[{pa,0,0},{pb,0,0},{0,0}],{pa>0,pb>0}]
Simplify[RatioP[{pa,0,0},{pb,0,0},{0,0}],{pa>0,pb<0, pa+pb>=0}]
Simplify[RatioP[{pa,0,0},{pb,0,0},{0,0}],{pa>0,pb<0, pa+pb<0}]

funToOvertones=Table[{2n,Simplify[ (((pa^(n-1)) (pb^(3-n)))/(Sqrt[1-(pb/pa)]))RatioM[{pa,0,0},{pb,0,0},{2n,0}],{pa>0,pb>0,pa>pb}]},{n,0,6}];
Table[funToOvertones[[1+n, 2]]==(2^(1-2n))*Binomial[2n,n]*(((I+n(2+2I)))/((1+2(1+I)n))) ,{n,0,6}]

Next, pa > 0 > pb, pa+pb>0

Simplify[RatioM[{pa,0,0},{pb,0,0},{0,0}],{pa>0,pb<0,pa+pb>0}]
Simplify[RatioM[{pa,0,0},{pb,0,0},{0,0}],{pa>0,pb>0, pa<pb}]
Simplify[RatioM[{pa,0,0},{pb,0,0},{0,0}],{pa>0,pb>0, pa>pb}]

funToOvertones2=Table[{2n,Simplify[ (((pa^(n-1)) (pb^(3-n)))/(((1+(pb/pa)))^(9/2)))RatioP[{pa,0,0},{pb,0,0},{2n,0}],{pa>0,pb<0,pa+pb>0}]},{n,0,6}];
Table[funToOvertones2[[1+n, 2]]==- Binomial[-1/2, n]*(((I+n(2+2I)))/((1+2(1+I)n)))  ,{n,0,6}]

Scaling limits

((((((pa^(n-1)) (pb^(3-n)))/(Sqrt[1-(pb/pa)])) ))^(-1))/. pb-> pa-δ;
Simplify[Series[%,{pa,∞,1}],{pa>0, pa>δ,δ>0}]

((((((pa^(n-1)) (pb^(3-n)))/(((1+(pb/pa)))^(9/2))) ))^(-1))/. pb-> δ-pa;
Simplify[Series[%,{pa,∞,5}],{pa>0, pa>δ,δ>0}]

3.2. Self-coupling of overtones

SelfCoup[{nx_,ny_},{NX_,NY_}] := Simplify[RatioP[{p,nx,ny},{p,nx,ny},{NX,NY}],{p>0}];

selfcoupmax=Monitor[Table[{n, SelfCoup[{n,0},{2n,0}]},{n,0,10}],n];
selfcoupself=Monitor[Table[{n, SelfCoup[{n,0},{n,0}]},{n,0,10}],n];
selfcoupzero=Monitor[Table[{n, SelfCoup[{n,0},{0,0}]},{n,0,10}],n];

3.3. Symmetry (n_x, n_y)

Monitor[Table[{SelfCoup[{k,0},{2k,0}]==Simplify[- Conjugate[SelfCoup[{0,k},{0,2k}]],p∈Reals]}, {k,0,4}],k]

Monitor[Table[{SelfCoup[{k,1},{2k,2}]==Simplify[- Conjugate[SelfCoup[{1,k},{2,2k}]],p∈Reals]}, {k,0,4}],k]

Monitor[Table[{SelfCoup[{k,2},{2k,4}]==Simplify[- Conjugate[SelfCoup[{2,k},{4,2k}]],p∈Reals]}, {k,0,4}],k]

3.4. Self-coupling of y overtones

selfcoupmax=Monitor[Table[{n, SelfCoup[{0,n},{0,2n}]},{n,0,10}],n];
selfcoupself=Monitor[Table[{n, SelfCoup[{0,n},{0,n}]},{n,0,10}],n];
selfcoupzero=Monitor[Table[{n, SelfCoup[{0,n},{0,0}]},{n,0,10}],n];