Integrating boson stars using the shooting method

In this notebook we integrate the differential equations for a spherically symmetric boson stars.  We divide the notebook into three different cases: scalar boson stars, vector boson stars, and charged boson stars.

Reference: Astrophysical signatures of boson stars: quasinormal modes and inspiral resonances , by Caio F. B. Macedo, Paolo Pani, Vitor Cardoso and Luís C. B. Crispino & Into the lair: gravitational-wave signatures of dark matter , by Caio F.B. Macedo, Paolo Pani, Vitor Cardoso, Luis C.B. Crispino

Scalar boson stars

Here we integrate the equations for spherically symmetric mini boson stars (potential with only the mass term). Let us start defining the equations of motion for the metric and the scalar field;

eqsbkg={((-1+B[r]+r B'[r])/(r^2))+((κ (((ωs)^2) ((ϕ0[r])^2)+A[r] (Vs[((ϕ0[r])^2)]+B[r] (((ϕ0)'[r])^2))))/(4 A[r])),((-1+B[r])/(r^2))-((κ ((ωs)^2) ((ϕ0[r])^2))/(4 A[r]))+((B[r] A'[r])/(r A[r]))+(1/4) κ (Vs[((ϕ0[r])^2)]-B[r] (((ϕ0)'[r])^2)),-ϕ0[r] Vs'[((ϕ0[r])^2)]+((2 B[r] (ϕ0)'[r])/r)+(1/2) B'[r] (ϕ0)'[r]+((((ωs)^2) ϕ0[r]+(1/2) B[r] A'[r] (ϕ0)'[r])/(A[r]))+B[r] (ϕ0)''[r]};
j0=((ωs ((ϕ0[r])^2))/(A[r]));

Collect[((eqsbkg[[3]])/(B[r])),{κ,(ϕ0)'[r],ϕ0[r],(ϕ0)''[r]},FullSimplify]

We shall use the potential with only the mass term and absorb the value of k in the definition of the scalar field. We also absorb the mass scale into the radial coordinate, effectively setting μ=1.

κ=1;
μ=1;
Vs[XX_]:=(μ^2)XX;

eqsbkg

now we jump for the boundary condition at the center of the star. We require regularity. We also set A(0)=1. While this choice makes the spacetime different from the standard Schwarzschild one at infinity, we can later one rescale the time coordinate to obtain the solution in the Schwarzchild coordinates. The simplest case is to use (B(0),A(0),ϕ0(0),ϕ0’(0))=(1,1,ϕc,0). We can also consider an expansion to improve convergence.

BOUNDARY={B[ror]==1,A[ror]==1,ϕ0[ror]==ϕc,ϕ0'[ror]==0};

Let us give the expression for the scalar field at infinity.

eqf=-ϕ0[r]+((2 B[r] (ϕ0)'[r])/r)+(1/2) B'[r] (ϕ0)'[r]+((((ωs)^2) ϕ0[r]+(1/2) B[r] A'[r] (ϕ0)'[r])/(A[r]))+B[r] (ϕ0)''[r]/.{B->(1&),A->(A0&),ϕ0->(((u[#])/(#))&)}//Simplify

DSolve[(((-A0+((ωs)^2)) u[r]+A0 u''[r])/(A0 r))==0,u,r]//Simplify

The wroskian of two solutions will therefore be

Wronskian[{(E^(-((r (Sqrt[A0-((ωs)^2)]))/(Sqrt[A0]))))/r,ϕ0[r]},r]//FullSimplify

Now we can define the integrator. The result of this integration will be the scalar field in a region far from the center. How far depends on how the scalar field grows or decays. The result will later on be used to find the proper value of the scalar field frequency to

INTB[ωs0_,ϕ0c_,ror0_,rinf0_]:=(
(*Defining rules for the parameters*)
param={ϕc->ϕ0c,ror->ror0,ωs->ωs0};
(*Setting the equations and the boundary conditions*)
EQS=Thread[eqsbkg=={0,0,0}]//.param;
BOU=BOUNDARY//.param;
(*Solving the equations*)
SOL=NDSolve[Join[EQS,BOU],{A,B,ϕ0},{r,ror0,rinf0},{AccuracyGoal->13,PrecisionGoal->12,MaxSteps->∞,
Method->"StiffnessSwitching",InterpolationOrder->All}][[1]]//Quiet;
RINF=Min[rinf0,(ϕ0/.SOL)[[1,1,2]]];
(*We can additionally provide the quantities that will aid to describe the spacetime in Schwarzschild coordinates*)
CC=((B[RINF])/(A[RINF]))/.SOL;
Ar=CC A[r]/.SOL;
Br=B[r]/.SOL;
ϕ0r=ϕ0[r]/.SOL;
Mt=-(1/2) (-1+B[r]) r/.SOL/.r->RINF;
mr=-(1/2) (-1+B[r]) r/.SOL;
ωss=ωs0(Sqrt[CC]);
(*Integrate to generate the particle number*) 
j0r=4π (r^2)(Sqrt[((A[r])/(B[r]))])j0/.SOL/.param;
NN0=((NN[RINF])/((16π)))/.NDSolve[{NN'[r]==j0r,NN[ror0]==0},NN,{r,ror0,RINF},{AccuracyGoal->10,PrecisionGoal->10}][[1]];
(*N0=4π NIntegrate[j0r (r^2),{r,ror0,RINF}];*)
)

Now, with the integration at hand, we can define the findroot. We shoot for the value of ωs in which the scalar field decays exponentially (the wroskian is zero). One can also shoot for the value of the scalar field going to zero

shoot[ωs0_?NumberQ,ϕ0c_,ror0_,rinf0_]:=(
INTB[ωs0,ϕ0c,ror0,rinf0];
(*ϕ0[r]/.SOL/.r->RINF*)
(r (Sqrt[A[r]-((ωs0)^2)]) ϕ0[r]+(Sqrt[A[r]]) (ϕ0[r]+r (ϕ0)'[r]))/.SOL/.r->RINF//Re
)

Let us finally try to find the solution for the boson star case. We start with a small value of the central scalar field. We can start with small values for the numerical infinity and slowly increasing the value until (desired) numerical accuracy is achieved

Monitor[FindRoot[shoot[xx,0.1,0.0001,40],{xx,1.0342939360249177`}],xx]
LogPlot[ϕ0r,{r,0.001,RINF}]
Mt
NN0

We can plot the metric and the scalar field for the configurations

Plot[{Ar,Br},{r,0.001,RINF}]

Plot[{ϕ0r},{r,0.001,RINF}]

From the value of the mass, we see that changes are within the fourth decimal place. We can therefore use this to find other solutions using as guess the initial one. To do so, we can also define the radius of the star to serve as a guide for the numerical infinity. We define the radius to be the one that contains 0.99 of the total mass.

FindRoot[-(1/2) (-1+Br) r==0.99Mt,{r,21}]

So we can use the infinity to be twice as large as the radius. Let us define the loop to find the solutions. This can take long because the step is small.

ϕii=0.01;
ϕfi=3;
pts=5000;
ωguesi=1.0032655739312906`;
Rdi=67.80808092320711`;

Monitor[Do[
(*Defining the steps*)
ϕit=Log[ϕii]+((Log[ϕfi]-Log[ϕii])/pts)i;
(*Defining the findroot for the solution and the radius*)
sww=FindRoot[shoot[xx,Exp[ϕit],0.001,1.8 Rdi],{xx,ωguesi}];
Rdi=r/.FindRoot[-(1/2) (-1+Br) r==0.99Mt,{r,Rdi}]//Quiet;
ωguesi=xx/.sww;
(*collecting the quatities*)
ϕi[i]=Exp[ϕit];
Mti[i]=Mt;
Rd[i]=Rdi;
ωii[i]=ωguesi;
ωiis[i]=ωss;
N0i[i]=NN0;
(*Tracking the points and plotting within the monitor*)
tabϕM=Table[{ϕi[j],Mti[j]},{j,0,i}];
mplot=ListPlot[tabϕM,PlotRange->All,Joined->True,Frame->True];
tabRM=Table[{Rd[j],Mti[j]},{j,0,i}];
mplot2=ListPlot[tabRM,PlotRange->All,Joined->True,Frame->True];
,{i,0,pts}],{i,mplot,mplot2}]

Exporting the points

texp=Table[{ϕi[j],Mti[j],Rd[j],N0i[j],ωiis[j],ωii[j]},{j,0,4501,5}];

SetDirectory[NotebookDirectory[]]

Charged boson stars

Quit[]

eqs={(1/4) (((κ ((((ωs+q A0[r]))^2) ((ϕ0[r])^2)+B[r] ((A0'[r])^2)))/(A[r]))+((-4+(r^2) κ Vs[((ϕ0[r])^2)]+4 r B'[r]+B[r] (4+(r^2) κ (((ϕ0)'[r])^2)))/(r^2))),((r (-r κ (((ωs+q A0[r]))^2) ((ϕ0[r])^2)+B[r] (4 A'[r]+r κ ((A0'[r])^2)))+A[r] (-4+(r^2) κ Vs[((ϕ0[r])^2)]+B[r] (4-(r^2) κ (((ϕ0)'[r])^2))))/(4 (r^2) A[r])),((2 r ϕ0[r] (((ωs)^2)+2 q ωs A0[r]+(q^2) ((A0[r])^2)-A[r] Vs'[((ϕ0[r])^2)])+r B[r] A'[r] (ϕ0)'[r]+A[r] (4 B[r] (ϕ0)'[r]+r B'[r] (ϕ0)'[r]+2 r B[r] (ϕ0)''[r]))/(2 r A[r])),-((q (ωs+q A0[r]) ((ϕ0[r])^2))/(B[r]))+(1/2) A0'[r] ((4/r)-((A'[r])/(A[r]))+((B'[r])/(B[r])))+A0''[r]};

Simplify[Solve[eqs==0/.B->(1-((m[#])/(#))&),{A'[r],m'[r],ϕ0''[r],A0''[r]}][[1]]]

((r κ (r (((ωs+q A0[r]))^2) ((ϕ0[r])^2)+(r-m[r]) ((A0'[r])^2)+A[r] (r (μ^2)ϕ0[r])))/(4 A[r]))//FullSimplify

We can change the variables to

A0[r_]:=(((-ωs+A0l[r])/q))

So that the equations will not depend on ωs directly

D[eqs,ωs]

Hence, we have the equations

eqsi=Simplify[Solve[eqs==0/.B->(1-((m[#])/(#))&),{A'[r],m'[r],ϕ0''[r],A0l''[r]}][[1]]]/.Rule->Subtract//FullSimplify

eqsi=Simplify[Solve[eqs==0,{A'[r],B'[r],ϕ0''[r],A0l''[r]}][[1]]]/.Rule->Subtract

We shall use units such that

κ=1;
μ=1;
Vs[XX_]:=(μ^2)XX;

and the particle density

j0=(((ωs+q A0[r]) ((ϕ0[r])^2))/(A[r]))//Simplify

Let us see how to implement the boundary condition at the origin. We take the expansions up to second order

Clear[A,B,ϕ0,A0l];

A[r_]:=Sum[(Subscript[Ai, i])(r^i),{i,0,2}];
B[r_]:=Sum[(Subscript[Bi, i])(r^i),{i,0,2}];
ϕ0[r_]:=Sum[(Subscript[ϕi, i])(r^i),{i,0,2}];
A0l[r_]:=Sum[(Subscript[A0i, i])(r^i),{i,0,2}];

rucoefs={(Subscript[Bi, 0])->1,(Subscript[ϕi, 1])->0,(Subscript[A0i, 1])->0,(Subscript[Ai, 1])->0,(Subscript[Bi, 1])->0,(Subscript[ϕi, 2])->-((((Subscript[A0i, 0]^2)-(Subscript[Ai, 0])) (Subscript[ϕi, 0]))/(6 (Subscript[Ai, 0]))),(Subscript[A0i, 2])->(1/6) (q^2) (Subscript[A0i, 0]) (Subscript[ϕi, 0]^2),(Subscript[Ai, 2])->(1/12) (2 (Subscript[A0i, 0]^2)-(Subscript[Ai, 0])) (Subscript[ϕi, 0]^2),(Subscript[Bi, 2])->-((((Subscript[A0i, 0]^2)+(Subscript[Ai, 0])) (Subscript[ϕi, 0]^2))/(12 (Subscript[Ai, 0])))};

Series[eqsi/.rucoefs,{r,0,1}]

Hence, we have the

bouor={A[r],B[r],ϕ0[r],ϕ0'[r],A0l[r],A0l'[r]}//.rucoefs

So we can use the above as boundary conditions. Now we have two possible choices. Either we use coordinate system in which A[0]=1 at the origin and shoot for A0l[0], or vice-versa. We pick the former.

Clear[A,B,ϕ0,A0l];

bou=Thread[{A[r],B[r],ϕ0[r],ϕ0'[r],A0l[r],A0l'[r]}==bouor/.{(Subscript[Ai, 0])->1,(Subscript[A0i, 0])->a0}];

At infinity, we have to extract the functions in order to obtain the mass and charge.

Clear[ϕ0]

A[r_]:=Sum[(Subscript[Ai, i])(r^(-i)),{i,0,2}];
B[r_]:=Sum[(Subscript[Bi, i])(r^(-i)),{i,0,2}];
A0l[r_]:=Sum[(Subscript[A0i, i])(r^(-i)),{i,0,1}];

rucoefs={(Subscript[Bi, 0])->1,(Subscript[Ai, 0])->1,(Subscript[Ai, 1])->(Subscript[Ai, 0]) (Subscript[Bi, 1]),(Subscript[Ai, 2])->((Subscript[A0i, 1]^2)/(4 (q^2))),(Subscript[Bi, 2])->((Subscript[A0i, 1]^2)/(4 (q^2) (Subscript[Ai, 0])))};

Series[eqsi/.{ϕ0->(0&)}//.rucoefs,{r,∞,3}]//Simplify

{B[r],A0l'[r],A0l[r]}//.rucoefs/.rucoefs/.(Subscript[Bi, 1])->-2M/.(Subscript[A0i, 1])->-Q 2q

((Subscript[A0i, 1]^2)/(4 (q^2) (r^2)))/.(Subscript[A0i, 1])->Q 2q

Clear[A,B,ϕ0,A0l];

Solve[{B[r]==1+((Q^2)/(r^2))-((2 M)/r),A0l'[r]==((2 q Q)/(r^2))},{M,Q}]//FullSimplify

Now we can define the integrator. As an input, we have the central value of the scalar field and field charge.

eqs

-((q A0l[r] ((ϕ0[r])^2))/(B[r]))+((A0l'[r] ((4/r)-((A'[r])/(A[r]))+((B'[r])/(B[r]))))/(2 q))+((A0l''[r])/q)/.{A0l->(-((2 q Qt)/(#))&),A->(1-2(M/(#))+((Q^2)/((#)^2))&),B->(1-2(M/(#))+((Q^2)/((#)^2))&)}

Series[((2 (q^2) Qt ((ϕ0[r])^2))/((1+((Q^2)/(r^2))-((2 M)/r)) r)),{r,∞,2}]

Solve[Qt== (((r^2) A0l'[r])/(2 q)),A0l'[r]]

A0l[r]->Integrate[((2 q Qt)/(r^2)),r]

INTB[ωs0_,q0_,ϕ0c_,ror0_,rinf0_]:=(
(*Defining rules for the parameters*)
param={(Subscript[ϕi, 0])->ϕ0c,ror->ror0,a0->ωs0,q->q0};
(*Setting the equations and the boundary conditions*)
EQS=Thread[eqsi=={0,0,0,0}]//.param;
BOU=bou//.param/.r->ror0;
(*Solving the equations*)
SOL=NDSolve[Join[EQS,BOU],{A,B,ϕ0,A0l},{r,ror0,rinf0},{AccuracyGoal->13,PrecisionGoal->12,MaxSteps->∞,Method->"StiffnessSwitching",InterpolationOrder->All}][[1]]//Quiet;
RINF=Min[rinf0,(ϕ0/.SOL)[[1,1,2]]];
(*We can additionally provide the quantities that will aid to describe the spacetime in Schwarzschild coordinates*)
CC=((B[RINF])/(A[RINF]))/.SOL;
Ar=CC A[r]/.SOL;
Br=B[r]/.SOL;
ϕ0r=ϕ0[r]/.SOL;
Qt=(((r^2) A0l'[r])/(2 q))/.SOL/.r->RINF/.param;
mr=(1/8) r (4-4 B[r]+(((r^2) ((A0l'[r])^2))/(q^2)))/.param/.SOL;
Mt=(1/8) r (4-4 B[r]+(((r^2) ((A0l'[r])^2))/(q^2)))/.SOL/.r->RINF/.param;
ωss=ωs0(Sqrt[CC]);
(*Integrate to generate the particle number*) 
j0r=4π (r^2)(Sqrt[((A[r])/(B[r]))])j0/.A0l[r]->((A0l[r])/(Sqrt[CC]))/.SOL/.param;
NN0=((NN[RINF])/((16π)))/.NDSolve[{NN'[r]==j0r,NN[ror0]==0},NN,{r,ror0,RINF},{AccuracyGoal->10,PrecisionGoal->10}][[1]];
(*N0=4π NIntegrate[j0r (r^2),{r,ror0,RINF}];*)
)

Now, with the integration at hand, we can define the findroot. We shoot for the value of ωs in which the scalar field decays exponentially

shoot[ωs0_?NumberQ,q0_,ϕ0c_,ror0_,rinf0_]:=(
INTB[ωs0,q0,ϕ0c,ror0,rinf0];
ϕ0[r]/.SOL/.r->RINF
)

Let us now see how to find the solution

SetDirectory[NotebookDirectory[]]

t1=Import["solutions/q045_chargedbs.dat"];

t1[[-1]]

Monitor[FindRoot[shoot[xx,0.45,0.3854649769437779`,0.001,55],{xx,1.083160829019893`}],xx]
LogPlot[ϕ0r//Abs,{r,0.001,RINF}]
Mt
NN0
Qt

NN0 2 0.45

FindRoot[mr==0.99Mt,{r,8}]

ϕii=0.01;
ϕfi=2;
pts=4000;
ωguesi=1.0000816850278809`;
Rdi=350;
q00=0.4999;
Mti[-1]=0;

Monitor[Do[
(*Defining the steps*)
ϕit=Log[ϕii]+((Log[ϕfi]-Log[ϕii])/pts)i;
(*Defining the findroot for the solution and the radius*)
sww=FindRoot[shoot[xx,q00,Exp[ϕit],0.001,1.8 Rdi],{xx,ωguesi}];
Rdi=r/.FindRoot[mr==0.99Mt,{r,Rdi}]//Quiet;
ωguesi=xx/.sww;
(*collecting the quatities*)
ϕi[i]=Exp[ϕit];
Mti[i]=Mt;
Qti[i]=Qt;
Rd[i]=Rdi;
ωii[i]=ωguesi;
ωiis[i]=ωss;
N0i[i]=NN0;
(*Tracking the points and plotting within the monitor*)
tabϕM=Table[{ϕi[j],Mti[j]},{j,0,i}];
mplot=ListPlot[tabϕM,PlotRange->All,Joined->True,Frame->True];
tabRM=Table[{Rd[j],Mti[j]},{j,0,i}];
mplot2=ListPlot[tabRM,PlotRange->All,Joined->True,Frame->True];
δM=SetPrecision[Mti[i]-Mti[i-1],5];
If[δM<0,ptsf=Length@tabRM;Break[];];
,{i,0,pts}],{i,mplot,mplot2,δM}]

texp=Table[{ϕi[j],Mti[j],Qti[j],Rd[j],N0i[j],ωiis[j],ωii[j]},{j,0,ptsf,3}];

ListPlot[texp[[All,{1,2}]]]

SetDirectory[NotebookDirectory[]]

Export["solutions/q04999_chargedbs.dat",texp]

Vector boson stars

Quit[]

Here we focus on only the mass case, without self-interaction

Vv[XX_]:=((μv)^2)XX;

eqsbkg={((-1+B[r]+r B'[r])/(r^2))+((κ (A[r] Vv[-(((f[r])^2)/(A[r]))+B[r] ((gf[r])^2)]+B[r] (((-ωv gf[r]+f'[r]))^2)+2 ((f[r])^2) Vv'[-(((f[r])^2)/(A[r]))+B[r] ((gf[r])^2)]))/(4 A[r])),((((4 (A[r] (-1+B[r])+r B[r] A'[r]))/(r^2))+κ (A[r] Vv[-(((f[r])^2)/(A[r]))+B[r] ((gf[r])^2)]+B[r] ((((-ωv gf[r]+f'[r]))^2)-2 A[r] ((gf[r])^2) Vv'[-(((f[r])^2)/(A[r]))+B[r] ((gf[r])^2)])))/(4 A[r])),(1/(2 r ((A[r])^2)))(r B[r] A'[r] (ωv gf[r]-f'[r])+A[r] (r (-ωv gf[r] B'[r]+B'[r] f'[r]-2 f[r] Vv'[-(((f[r])^2)/(A[r]))+B[r] ((gf[r])^2)])+B[r] (-4 ωv gf[r]+4 f'[r]-2 r ωv gf'[r]+2 r f''[r]))),((I B[r] (ωv f'[r]+gf[r] (-((ωv)^2)+A[r] Vv'[-(((f[r])^2)/(A[r]))+B[r] ((gf[r])^2)])))/(A[r])),(1/(2 r ((A[r])^2)))I (E^(-I t ωv)) (2 r B[r] ((f[r])^2) gf[r] A'[r] Vv''[-(((f[r])^2)/(A[r]))+B[r] ((gf[r])^2)]+((A[r])^2) (r gf[r] B'[r] Vv'[-(((f[r])^2)/(A[r]))+B[r] ((gf[r])^2)]+4 r ((B[r])^2) ((gf[r])^2) gf'[r] Vv''[-(((f[r])^2)/(A[r]))+B[r] ((gf[r])^2)]+2 B[r] (2 gf[r] Vv'[-(((f[r])^2)/(A[r]))+B[r] ((gf[r])^2)]+r gf'[r] Vv'[-(((f[r])^2)/(A[r]))+B[r] ((gf[r])^2)]+r ((gf[r])^3) B'[r] Vv''[-(((f[r])^2)/(A[r]))+B[r] ((gf[r])^2)]))+r A[r] (B[r] gf[r] A'[r] Vv'[-(((f[r])^2)/(A[r]))+B[r] ((gf[r])^2)]+2 f[r] (ωv Vv'[-(((f[r])^2)/(A[r]))+B[r] ((gf[r])^2)]-2 B[r] gf[r] f'[r] Vv''[-(((f[r])^2)/(A[r]))+B[r] ((gf[r])^2)])))};
j0=((B[r] gf[r] (-ωv gf[r]+f'[r]))/(A[r]));

eqsbkg[[1]]/.{f->((Sqrt[(1/κ)])f[#]&),gf->((Sqrt[(1/κ)])gf[#]&)}//Simplify

((I B[r] (ωv f'[r]+gf[r] (-((ωv)^2)+A[r] Vv'[-(((f[r])^2)/(A[r]))+B[r] ((gf[r])^2)])))/(A[r]))//Simplify

rufaa={f'[r]->"f'",f''[r]->"f''",gf[r]->"g",gf'[r]->"g'",A[r]->"A",B[r]->"B",B'[r]->"B'",f[r]->"f",ωv->"(Subscript[ω, V])",μv->"(Subscript[μ, V])"};

FullSimplify[((A[r])/(I B[r]ωv))eqsbkg[[4]]]/.{f'[r]->"f'",f''[r]->"f''",gf[r]->"g",gf'[r]->"g'",A[r]->"A",B[r]->"B",B'[r]->"B'",f[r]->"f",ωv->"(Subscript[ω, V])",μv->"(Subscript[μ, V])"}

((κ (((μv)^2) A[r] (-(((f[r])^2)/(A[r]))+B[r] ((gf[r])^2))+B[r] (-2 ((μv)^2) A[r] ((gf[r])^2)+(((-ωv gf[r]+f'[r]))^2))))/(4 A[r]))/.{f'[r]->"f'",gf[r]->"g",A[r]->"A",B[r]->"B",f[r]->"f",ωv->"(Subscript[ω, V])",μv->"(Subscript[μ, V])"}

Collect[((κ (2 ((μv)^2) ((f[r])^2)+((μv)^2) A[r] (-(((f[r])^2)/(A[r]))+B[r] ((gf[r])^2))+B[r] (((-ωv gf[r]+f'[r]))^2)))/(4 A[r])),{κ,((μv)^2)},FullSimplify]

(((ωv)^2)) ((gf[r])^2)-2 ωv gf[r] f'[r]+((f'[r])^2)//FullSimplify

(∂_r)((κ (((μv)^2) ((f[r])^2)+B[r] ((((ωv)^2)+((μv)^2) A[r]) ((gf[r])^2)-2 ωv gf[r] f'[r]+((f'[r])^2))))/(4 A[r]))

Collect[((eqsbkg[[-1]])/(I (E^(-I t ωv)) ((μv)^2) B[r] )),gf'[r],Simplify]/.rufaa

((eqsbkg[[-1]])/(I (E^(-I t ωv)) ((μv)^2) B[r] ))/.rufaa//FullSimplify

bkg=Solve[eqsbkg==0,{A'[r],B'[r],f'[r],gf'[r],f''[r]}][[1,1;;4]]/.Rule->Subtract//Simplify

bkg[[1]]

Let us find the boundary conditions at the origin

Clear[f,gf,A,B];

rucoefis={(Subscript[Ai, 1])->0,(Subscript[Bi, 1])->0,(Subscript[fi, 1])->0,(Subscript[gfi, 1])->-((ωv (Subscript[fi, 0]))/(3 (Subscript[Ai, 0])))};

f[r_]:=Sum[(Subscript[fi, i])(r^i),{i,0,1}]/.rucoefis;
gf[r_]:=Sum[(Subscript[gfi, i])(r^i),{i,0,1}]/.(Subscript[gfi, 0])->0/.rucoefis;
A[r_]:=Sum[(Subscript[Ai, i])(r^i),{i,0,1}]/.rucoefis;
B[r_]:=Sum[(Subscript[Bi, i])(r^i),{i,0,1}]/.(Subscript[Bi, 0])->1/.rucoefis;

Series[bkg,{r,0,0}]

Solve[%==0,{(Subscript[Ai, 1]),(Subscript[Bi, 1]),(Subscript[fi, 1]),(Subscript[gfi, 1])}]

bouori={A[r],B[r],f[r],gf[r]}/.(Subscript[Ai, 0])->1

Clear[f,gf,A,B];

BOUNDARY=Thread[{A[r],B[r],f[r],gf[r]}==bouori]/.r->ror

We set units such that μv=1 and κ=1;

κ=1;
μv=1;

bkg

Now we can define the integrator. The result of this integration will be the scalar field in a region far from the center. How far depends on how the scalar field grows or decays. The result will later on be used to find the proper value of the scalar field frequency to

INTB[ωs0_,ϕ0c_,ror0_,rinf0_]:=(
(*Defining rules for the parameters*)
param={(Subscript[fi, 0])->ϕ0c,ror->ror0,ωv->ωs0};
(*Setting the equations and the boundary conditions*)
EQS=Thread[bkg=={0,0,0,0}]//.param;
BOU=BOUNDARY//.param;
(*Solving the equations*)
SOL=NDSolve[Join[EQS,BOU],{A,B,f,gf},{r,ror0,rinf0},{AccuracyGoal->13,PrecisionGoal->12,MaxSteps->∞,Method->"StiffnessSwitching",InterpolationOrder->All}][[1]]//Quiet;
RINF=Min[rinf0,(f/.SOL)[[1,1,2]]];
(*We can additionally provide the quantities that will aid to describe the spacetime in Schwarzschild coordinates*)
CC=((B[RINF])/(A[RINF]))/.SOL;
Ar=CC A[r]/.SOL;
Br=B[r]/.SOL;
fr=f0[r]/.SOL;
gfr=gf[r]/.SOL;
Mt=-(1/2) (-1+B[r]) r/.SOL/.r->RINF;
mr=-(1/2) (-1+B[r]) r/.SOL;
ωss=ωs0(Sqrt[CC]);
(*Integrate to generate the particle number*) 
j0r=4π (r^2)(Sqrt[((A[r])/(B[r]))])j0/.SOL/.param;
NN0=-((NN[RINF])/((16π)))/.NDSolve[{NN'[r]==j0r,NN[ror0]==0},NN,{r,ror0,RINF},{AccuracyGoal->10,PrecisionGoal->10}][[1]];
)

Now, with the integration at hand, we can define the findroot. We shoot for the value of ωs in which the scalar field decays exponentially

shoot[ωs0_?NumberQ,ϕ0c_,ror0_,rinf0_]:=(
INTB[ωs0,ϕ0c,ror0,rinf0];
f[r]/.SOL/.r->RINF/.param
)

FindRoot[shoot[xx,0.001,0.001,150],{xx,1.0032743561481414}]
LogPlot[f[r]/.SOL//Abs,{r,0.001,RINF},PlotRange->All]

Mt

NN0

FindRoot[-(1/2) (-1+Br) r==0.99Mt,{r,80}]

ϕii=0.001;
ϕfi=15;
pts=5000;
ωguesi=1.003274704035908`;
Rdi=85;

Monitor[Do[
(*Defin
ing the steps*)
ϕit=Log[ϕii]+((Log[ϕfi]-Log[ϕii])/pts)i;
(*Defining the findroot for the solution and the radius*)
sww=FindRoot[shoot[xx,Exp[ϕit],0.001,1.8 Rdi],{xx,ωguesi}];
Rdi=r/.FindRoot[-(1/2) (-1+Br) r==0.99Mt,{r,Rdi}]//Quiet;
ωguesi=xx/.sww;
(*collecting the quatities*)
ϕi[i]=Exp[ϕit];
Mti[i]=Mt;
Rd[i]=Rdi;
ωii[i]=ωguesi;
ωiis[i]=ωss;
N0i[i]=NN0;
(*Tracking the points and plotting within the monitor*)
tabϕM=Table[{ϕi[j],Mti[j]},{j,0,i}];
mplot=ListPlot[tabϕM,PlotRange->All,Joined->True,Frame->True];
tabRM=Table[{Rd[j],Mti[j]},{j,0,i}];
mplot2=ListPlot[tabRM,PlotRange->All,Joined->True,Frame->True];
,{i,0,pts}],{i,mplot,mplot2}]

ListLogLinearPlot[tabϕM[[1;;-1;;100]]]

mplot2

texp=Table[{ϕi[j],Mti[j],Rd[j],N0i[j],ωiis[j],ωii[j]},{j,0,pts,5}];

SetDirectory[NotebookDirectory[]]

Export["solutions\vectorminiBS_parameters_fc_M_R_N_w.dat",texp]

We can see from the above that we have a maximum mass of

Max[tabRM[[All,2]]]

Oscillatons

Quit[]

eqsbkg={(1/(32 (r^2)))(16 A03[r] B02[r]-4 (r^2) κ (ω^2) ((ϕ01[r])^2)+2 (r^2) κ (μ^2) A03[r] ((ϕ01[r])^2)+24 (r^2) κ (ω^2) ϕ01[r] ϕ02[r]+4 (r^2) κ (μ^2) A03[r] ϕ01[r] ϕ02[r]+2 (r^2) κ (μ^2) A03[r] ((ϕ02[r])^2)+16 r A03[r] B02'[r]+2 (r^2) κ A03[r] B01[r] (((ϕ01)'[r])^2)+2 (r^2) κ A03[r] B02[r] (((ϕ01)'[r])^2)+2 (r^2) κ A03[r] B03[r] (((ϕ01)'[r])^2)+4 (r^2) κ A03[r] B01[r] (ϕ01)'[r] (ϕ02)'[r]+4 (r^2) κ A03[r] B02[r] (ϕ01)'[r] (ϕ02)'[r]+4 (r^2) κ A03[r] B03[r] (ϕ01)'[r] (ϕ02)'[r]+2 (r^2) κ A03[r] B01[r] (((ϕ02)'[r])^2)+2 (r^2) κ A03[r] B02[r] (((ϕ02)'[r])^2)+(r^2) κ A03[r] B03[r] (((ϕ02)'[r])^2)+A02[r] (-32+4 (r^2) κ (μ^2) ((ϕ01[r])^2)+4 (r^2) κ (μ^2) ϕ01[r] ϕ02[r]+4 (r^2) κ (μ^2) ((ϕ02[r])^2)+32 r B01'[r]+16 r B03'[r]+3 (r^2) κ B02[r] (((ϕ01)'[r])^2)+6 (r^2) κ B02[r] (ϕ01)'[r] (ϕ02)'[r]+(r^2) κ B02[r] (((ϕ02)'[r])^2)+4 B01[r] (8+(r^2) κ (((ϕ01)'[r])^2)+(r^2) κ (ϕ01)'[r] (ϕ02)'[r]+(r^2) κ (((ϕ02)'[r])^2))+2 B03[r] (8+(r^2) κ (((ϕ01)'[r])^2)+2 (r^2) κ (ϕ01)'[r] (ϕ02)'[r]+(r^2) κ (((ϕ02)'[r])^2)))+2 A01[r] (2 B02[r] (8+(r^2) κ (((ϕ01)'[r])^2)+(r^2) κ (ϕ01)'[r] (ϕ02)'[r]+(r^2) κ (((ϕ02)'[r])^2))+r (2 r κ (μ^2) ((ϕ01[r])^2)+4 r κ (μ^2) ϕ01[r] ϕ02[r]+16 B02'[r]+r κ (B03[r] ((((ϕ01)'[r]+(ϕ02)'[r]))^2)+2 B01[r] (ϕ01)'[r] ((ϕ01)'[r]+2 (ϕ02)'[r]))))),(1/(32 (r^2)))(-24 (r^2) κ (ω^2) ϕ01[r] ϕ02[r]+A03[r] (-32+4 (r^2) κ (μ^2) ((ϕ01[r])^2)+4 (r^2) κ (μ^2) ((ϕ02[r])^2)+32 r B01'[r]+2 (r^2) κ B02[r] (((ϕ01)'[r])^2)+4 (r^2) κ B02[r] (ϕ01)'[r] (ϕ02)'[r]+6 (r^2) κ B03[r] (ϕ01)'[r] (ϕ02)'[r]+(r^2) κ B02[r] (((ϕ02)'[r])^2)+4 B01[r] (8+(r^2) κ (((ϕ01)'[r])^2)+(r^2) κ (((ϕ02)'[r])^2)))+A02[r] (2 B02[r] (8+(r^2) κ (((ϕ01)'[r])^2)+2 (r^2) κ (ϕ01)'[r] (ϕ02)'[r]+(r^2) κ (((ϕ02)'[r])^2))+r (2 r κ (μ^2) ((ϕ01[r])^2)+4 r κ (μ^2) ϕ01[r] ϕ02[r]+2 r κ (μ^2) ((ϕ02[r])^2)+16 B02'[r]+2 r κ B01[r] (((ϕ01)'[r])^2)+2 r κ B03[r] (((ϕ01)'[r])^2)+4 r κ B01[r] (ϕ01)'[r] (ϕ02)'[r]+4 r κ B03[r] (ϕ01)'[r] (ϕ02)'[r]+2 r κ B01[r] (((ϕ02)'[r])^2)+r κ B03[r] (((ϕ02)'[r])^2)))+2 A01[r] (2 B03[r] (8+(r^2) κ (((ϕ01)'[r])^2)+(r^2) κ (((ϕ02)'[r])^2))+r (4 r κ (μ^2) ϕ01[r] ϕ02[r]+16 B03'[r]+r κ (4 B01[r] (ϕ01)'[r] (ϕ02)'[r]+B02[r] ((((ϕ01)'[r]+(ϕ02)'[r]))^2))))),(1/(32 (r^2)))(16 A03[r] B03[r]+4 (r^2) κ (ω^2) ((ϕ01[r])^2)+4 (r^2) κ (μ^2) A03[r] ϕ01[r] ϕ02[r]+36 (r^2) κ (ω^2) ((ϕ02[r])^2)+16 r A03[r] B03'[r]+(r^2) κ A03[r] B02[r] (((ϕ01)'[r])^2)+2 (r^2) κ A03[r] B03[r] (((ϕ01)'[r])^2)+4 (r^2) κ A03[r] B01[r] (ϕ01)'[r] (ϕ02)'[r]+2 (r^2) κ A03[r] B02[r] (ϕ01)'[r] (ϕ02)'[r]+(r^2) κ A03[r] B02[r] (((ϕ02)'[r])^2)+2 (r^2) κ A03[r] B03[r] (((ϕ02)'[r])^2)+2 A01[r] (-16+2 (r^2) κ (μ^2) ((ϕ01[r])^2)+2 (r^2) κ (μ^2) ((ϕ02[r])^2)+16 r B01'[r]+(r^2) κ B02[r] (((ϕ01)'[r])^2)+2 (r^2) κ B02[r] (ϕ01)'[r] (ϕ02)'[r]+2 (r^2) κ B03[r] (ϕ01)'[r] (ϕ02)'[r]+2 B01[r] (8+(r^2) κ (((ϕ01)'[r])^2)+(r^2) κ (((ϕ02)'[r])^2)))+A02[r] (2 B02[r] (8+(r^2) κ (((ϕ01)'[r])^2)+(r^2) κ (ϕ01)'[r] (ϕ02)'[r]+(r^2) κ (((ϕ02)'[r])^2))+r (2 r κ (μ^2) ((ϕ01[r])^2)+4 r κ (μ^2) ϕ01[r] ϕ02[r]+16 B02'[r]+r κ (B03[r] ((((ϕ01)'[r]+(ϕ02)'[r]))^2)+2 B01[r] (ϕ01)'[r] ((ϕ01)'[r]+2 (ϕ02)'[r]))))),-(1/(32 (r^2)))(-16 A03[r] B02[r]-4 (r^2) κ (ω^2) ((ϕ01[r])^2)-2 (r^2) κ (μ^2) A03[r] ((ϕ01[r])^2)+24 (r^2) κ (ω^2) ϕ01[r] ϕ02[r]-4 (r^2) κ (μ^2) A03[r] ϕ01[r] ϕ02[r]-2 (r^2) κ (μ^2) A03[r] ((ϕ02[r])^2)-32 r B02[r] A01'[r]-32 r B01[r] A02'[r]-16 r B03[r] A02'[r]-16 r B02[r] A03'[r]+2 (r^2) κ A03[r] B01[r] (((ϕ01)'[r])^2)+2 (r^2) κ A03[r] B02[r] (((ϕ01)'[r])^2)+2 (r^2) κ A03[r] B03[r] (((ϕ01)'[r])^2)+4 (r^2) κ A03[r] B01[r] (ϕ01)'[r] (ϕ02)'[r]+4 (r^2) κ A03[r] B02[r] (ϕ01)'[r] (ϕ02)'[r]+4 (r^2) κ A03[r] B03[r] (ϕ01)'[r] (ϕ02)'[r]+2 (r^2) κ A03[r] B01[r] (((ϕ02)'[r])^2)+2 (r^2) κ A03[r] B02[r] (((ϕ02)'[r])^2)+(r^2) κ A03[r] B03[r] (((ϕ02)'[r])^2)+A02[r] (32-4 (r^2) κ (μ^2) ((ϕ01[r])^2)-4 (r^2) κ (μ^2) ϕ01[r] ϕ02[r]-4 (r^2) κ (μ^2) ((ϕ02[r])^2)+3 (r^2) κ B02[r] (((ϕ01)'[r])^2)+6 (r^2) κ B02[r] (ϕ01)'[r] (ϕ02)'[r]+(r^2) κ B02[r] (((ϕ02)'[r])^2)+4 B01[r] (-8+(r^2) κ (((ϕ01)'[r])^2)+(r^2) κ (ϕ01)'[r] (ϕ02)'[r]+(r^2) κ (((ϕ02)'[r])^2))+2 B03[r] (-8+(r^2) κ (((ϕ01)'[r])^2)+2 (r^2) κ (ϕ01)'[r] (ϕ02)'[r]+(r^2) κ (((ϕ02)'[r])^2)))+2 A01[r] (2 B02[r] (-8+(r^2) κ (((ϕ01)'[r])^2)+(r^2) κ (ϕ01)'[r] (ϕ02)'[r]+(r^2) κ (((ϕ02)'[r])^2))+(r^2) κ (-2 (μ^2) ((ϕ01[r])^2)-4 (μ^2) ϕ01[r] ϕ02[r]+B03[r] ((((ϕ01)'[r]+(ϕ02)'[r]))^2)+2 B01[r] (ϕ01)'[r] ((ϕ01)'[r]+2 (ϕ02)'[r])))),-(1/(32 (r^2)))(-8 r (3 r κ (ω^2) ϕ01[r] ϕ02[r]+4 B03[r] A01'[r]+2 B02[r] A02'[r]+4 B01[r] A03'[r])+A03[r] (32-4 (r^2) κ (μ^2) ((ϕ01[r])^2)-4 (r^2) κ (μ^2) ((ϕ02[r])^2)+2 (r^2) κ B02[r] (((ϕ01)'[r])^2)+4 (r^2) κ B02[r] (ϕ01)'[r] (ϕ02)'[r]+6 (r^2) κ B03[r] (ϕ01)'[r] (ϕ02)'[r]+(r^2) κ B02[r] (((ϕ02)'[r])^2)+4 B01[r] (-8+(r^2) κ (((ϕ01)'[r])^2)+(r^2) κ (((ϕ02)'[r])^2)))+A02[r] (2 B02[r] (-8+(r^2) κ (((ϕ01)'[r])^2)+2 (r^2) κ (ϕ01)'[r] (ϕ02)'[r]+(r^2) κ (((ϕ02)'[r])^2))+(r^2) κ (-2 (μ^2) ((ϕ01[r])^2)-4 (μ^2) ϕ01[r] ϕ02[r]-2 (μ^2) ((ϕ02[r])^2)+2 B01[r] (((ϕ01)'[r])^2)+2 B03[r] (((ϕ01)'[r])^2)+4 B01[r] (ϕ01)'[r] (ϕ02)'[r]+4 B03[r] (ϕ01)'[r] (ϕ02)'[r]+2 B01[r] (((ϕ02)'[r])^2)+B03[r] (((ϕ02)'[r])^2)))+2 A01[r] (2 B03[r] (-8+(r^2) κ (((ϕ01)'[r])^2)+(r^2) κ (((ϕ02)'[r])^2))+(r^2) κ (-4 (μ^2) ϕ01[r] ϕ02[r]+4 B01[r] (ϕ01)'[r] (ϕ02)'[r]+B02[r] ((((ϕ01)'[r]+(ϕ02)'[r]))^2)))),-(1/(32 (r^2)))(-16 A03[r] B03[r]+4 (r^2) κ (ω^2) ((ϕ01[r])^2)-4 (r^2) κ (μ^2) A03[r] ϕ01[r] ϕ02[r]+36 (r^2) κ (ω^2) ((ϕ02[r])^2)-32 r B01[r] A01'[r]-16 r B02[r] A02'[r]-16 r B03[r] A03'[r]+(r^2) κ A03[r] B02[r] (((ϕ01)'[r])^2)+2 (r^2) κ A03[r] B03[r] (((ϕ01)'[r])^2)+4 (r^2) κ A03[r] B01[r] (ϕ01)'[r] (ϕ02)'[r]+2 (r^2) κ A03[r] B02[r] (ϕ01)'[r] (ϕ02)'[r]+(r^2) κ A03[r] B02[r] (((ϕ02)'[r])^2)+2 (r^2) κ A03[r] B03[r] (((ϕ02)'[r])^2)+2 A01[r] (16-2 (r^2) κ (μ^2) ((ϕ01[r])^2)-2 (r^2) κ (μ^2) ((ϕ02[r])^2)+(r^2) κ B02[r] (((ϕ01)'[r])^2)+2 (r^2) κ B02[r] (ϕ01)'[r] (ϕ02)'[r]+2 (r^2) κ B03[r] (ϕ01)'[r] (ϕ02)'[r]+2 B01[r] (-8+(r^2) κ (((ϕ01)'[r])^2)+(r^2) κ (((ϕ02)'[r])^2)))+A02[r] (2 B02[r] (-8+(r^2) κ (((ϕ01)'[r])^2)+(r^2) κ (ϕ01)'[r] (ϕ02)'[r]+(r^2) κ (((ϕ02)'[r])^2))+(r^2) κ (-2 (μ^2) ((ϕ01[r])^2)-4 (μ^2) ϕ01[r] ϕ02[r]+B03[r] ((((ϕ01)'[r]+(ϕ02)'[r]))^2)+2 B01[r] (ϕ01)'[r] ((ϕ01)'[r]+2 (ϕ02)'[r])))),(1/(32 r))(((A02[r])^2) (8 ((B01[r])^2) (4 (ϕ01)'[r]+2 (ϕ02)'[r]+2 r (ϕ01)''[r]+r (ϕ02)''[r])+4 ((B02[r])^2) (6 (ϕ01)'[r]+4 (ϕ02)'[r]+3 r (ϕ01)''[r]+2 r (ϕ02)''[r])+B03[r] (-8 r (μ^2) ϕ01[r]-8 r (μ^2) ϕ02[r]+16 B03[r] (ϕ01)'[r]+4 r B01'[r] (ϕ01)'[r]+4 r B02'[r] (ϕ01)'[r]+4 r B03'[r] (ϕ01)'[r]+12 B03[r] (ϕ02)'[r]+4 r B01'[r] (ϕ02)'[r]+4 r B02'[r] (ϕ02)'[r]+3 r B03'[r] (ϕ02)'[r]+8 r B03[r] (ϕ01)''[r]+6 r B03[r] (ϕ02)''[r])+2 B02[r] (-6 r (μ^2) ϕ01[r]-6 r (μ^2) ϕ02[r]+16 B03[r] (ϕ01)'[r]+3 r B01'[r] (ϕ01)'[r]+3 r B02'[r] (ϕ01)'[r]+2 r B03'[r] (ϕ01)'[r]+16 B03[r] (ϕ02)'[r]+3 r B01'[r] (ϕ02)'[r]+2 r B02'[r] (ϕ02)'[r]+2 r B03'[r] (ϕ02)'[r]+8 r B03[r] (ϕ01)''[r]+8 r B03[r] (ϕ02)''[r])+2 B01[r] (-8 r (μ^2) ϕ01[r]-4 r (μ^2) ϕ02[r]+24 B02[r] (ϕ01)'[r]+16 B03[r] (ϕ01)'[r]+4 r B01'[r] (ϕ01)'[r]+3 r B02'[r] (ϕ01)'[r]+2 r B03'[r] (ϕ01)'[r]+24 B02[r] (ϕ02)'[r]+16 B03[r] (ϕ02)'[r]+2 r B01'[r] (ϕ02)'[r]+3 r B02'[r] (ϕ02)'[r]+2 r B03'[r] (ϕ02)'[r]+12 r B02[r] (ϕ01)''[r]+8 r B03[r] (ϕ01)''[r]+12 r B02[r] (ϕ02)''[r]+8 r B03[r] (ϕ02)''[r]))+A02[r] (4 B01[r] ((8 r (ω^2)-4 r (μ^2) A03[r]) ϕ01[r]+(48 r (ω^2)-4 r (μ^2) A03[r]) ϕ02[r]+16 A03[r] B02[r] (ϕ01)'[r]+16 A03[r] B03[r] (ϕ01)'[r]+4 r B02[r] A01'[r] (ϕ01)'[r]+2 r B03[r] A01'[r] (ϕ01)'[r]+3 r B02[r] A02'[r] (ϕ01)'[r]+2 r B03[r] A02'[r] (ϕ01)'[r]+2 r B02[r] A03'[r] (ϕ01)'[r]+2 r B03[r] A03'[r] (ϕ01)'[r]+2 r A03[r] B01'[r] (ϕ01)'[r]+2 r A03[r] B02'[r] (ϕ01)'[r]+2 r A03[r] B03'[r] (ϕ01)'[r]+16 A03[r] B02[r] (ϕ02)'[r]+16 A03[r] B03[r] (ϕ02)'[r]+2 r B02[r] A01'[r] (ϕ02)'[r]+2 r B03[r] A01'[r] (ϕ02)'[r]+3 r B02[r] A02'[r] (ϕ02)'[r]+2 r B03[r] A02'[r] (ϕ02)'[r]+2 r B02[r] A03'[r] (ϕ02)'[r]+2 r B03[r] A03'[r] (ϕ02)'[r]+2 r A03[r] B01'[r] (ϕ02)'[r]+2 r A03[r] B02'[r] (ϕ02)'[r]+2 r A03[r] B03'[r] (ϕ02)'[r]+8 r A03[r] B02[r] (ϕ01)''[r]+8 r A03[r] B03[r] (ϕ01)''[r]+8 r A03[r] B02[r] (ϕ02)''[r]+8 r A03[r] B03[r] (ϕ02)''[r])+B03[r] (16 r ((ω^2)-(μ^2) A03[r]) ϕ01[r]+16 r (6 (ω^2)-(μ^2) A03[r]) ϕ02[r]+24 A03[r] B03[r] (ϕ01)'[r]+4 r B03[r] A01'[r] (ϕ01)'[r]+4 r B03[r] A02'[r] (ϕ01)'[r]+3 r B03[r] A03'[r] (ϕ01)'[r]+8 r A03[r] B01'[r] (ϕ01)'[r]+8 r A03[r] B02'[r] (ϕ01)'[r]+6 r A03[r] B03'[r] (ϕ01)'[r]+24 A03[r] B03[r] (ϕ02)'[r]+4 r B03[r] A01'[r] (ϕ02)'[r]+3 r B03[r] A02'[r] (ϕ02)'[r]+3 r B03[r] A03'[r] (ϕ02)'[r]+8 r A03[r] B01'[r] (ϕ02)'[r]+6 r A03[r] B02'[r] (ϕ02)'[r]+6 r A03[r] B03'[r] (ϕ02)'[r]+12 r A03[r] B03[r] (ϕ01)''[r]+12 r A03[r] B03[r] (ϕ02)''[r])+2 B02[r] (8 r ((ω^2)-(μ^2) A03[r]) ϕ01[r]+(60 r (ω^2)-8 r (μ^2) A03[r]) ϕ02[r]+32 A03[r] B03[r] (ϕ01)'[r]+4 r B03[r] A01'[r] (ϕ01)'[r]+4 r B03[r] A02'[r] (ϕ01)'[r]+4 r B03[r] A03'[r] (ϕ01)'[r]+4 r A03[r] B01'[r] (ϕ01)'[r]+4 r A03[r] B02'[r] (ϕ01)'[r]+4 r A03[r] B03'[r] (ϕ01)'[r]+24 A03[r] B03[r] (ϕ02)'[r]+4 r B03[r] A01'[r] (ϕ02)'[r]+4 r B03[r] A02'[r] (ϕ02)'[r]+3 r B03[r] A03'[r] (ϕ02)'[r]+4 r A03[r] B01'[r] (ϕ02)'[r]+4 r A03[r] B02'[r] (ϕ02)'[r]+3 r A03[r] B03'[r] (ϕ02)'[r]+16 r A03[r] B03[r] (ϕ01)''[r]+12 r A03[r] B03[r] (ϕ02)''[r])+4 ((B01[r])^2) (r (2 A01'[r] ((ϕ01)'[r]+(ϕ02)'[r])+A03'[r] ((ϕ01)'[r]+(ϕ02)'[r])+A02'[r] (2 (ϕ01)'[r]+(ϕ02)'[r]))+4 A03[r] (2 (ϕ01)'[r]+2 (ϕ02)'[r]+r ((ϕ01)''[r]+(ϕ02)''[r])))+2 ((B02[r])^2) (r (3 A01'[r] ((ϕ01)'[r]+(ϕ02)'[r])+2 A03'[r] ((ϕ01)'[r]+(ϕ02)'[r])+A02'[r] (3 (ϕ01)'[r]+2 (ϕ02)'[r]))+8 A03[r] (2 (ϕ01)'[r]+2 (ϕ02)'[r]+r ((ϕ01)''[r]+(ϕ02)''[r]))))+4 ((A01[r])^2) (8 ((B01[r])^2) (2 (ϕ01)'[r]+r (ϕ01)''[r])+B03[r] (-4 r (μ^2) ϕ02[r]+r B02'[r] ((ϕ01)'[r]+(ϕ02)'[r])+2 r (B03'[r] (ϕ01)'[r]+B01'[r] (ϕ02)'[r])+4 B03[r] (2 (ϕ01)'[r]+r (ϕ01)''[r]))+B02[r] (-4 r (μ^2) ϕ01[r]-4 r (μ^2) ϕ02[r]+8 B03[r] (ϕ01)'[r]+2 r B01'[r] (ϕ01)'[r]+2 r B02'[r] (ϕ01)'[r]+r B03'[r] (ϕ01)'[r]+8 B03[r] (ϕ02)'[r]+2 r B01'[r] (ϕ02)'[r]+r B02'[r] (ϕ02)'[r]+r B03'[r] (ϕ02)'[r]+4 r B03[r] (ϕ01)''[r]+4 r B03[r] (ϕ02)''[r])+2 ((B02[r])^2) (4 (ϕ01)'[r]+2 (ϕ02)'[r]+r (2 (ϕ01)''[r]+(ϕ02)''[r]))+2 B01[r] (-4 r (μ^2) ϕ01[r]+2 r B01'[r] (ϕ01)'[r]+r B02'[r] (ϕ01)'[r]+8 B03[r] (ϕ02)'[r]+r B02'[r] (ϕ02)'[r]+r B03'[r] (ϕ02)'[r]+4 r B03[r] (ϕ02)''[r]+4 B02[r] (2 (ϕ01)'[r]+2 (ϕ02)'[r]+r ((ϕ01)''[r]+(ϕ02)''[r]))))+2 A01[r] (2 B02[r] (4 r (2 (ω^2)-2 (μ^2) A02[r]-(μ^2) A03[r]) ϕ01[r]+4 r (12 (ω^2)-(μ^2) A02[r]-(μ^2) A03[r]) ϕ02[r]+16 A02[r] B03[r] (ϕ01)'[r]+16 A03[r] B03[r] (ϕ01)'[r]+2 r B03[r] A01'[r] (ϕ01)'[r]+2 r B03[r] A02'[r] (ϕ01)'[r]+2 r B03[r] A03'[r] (ϕ01)'[r]+4 r A02[r] B01'[r] (ϕ01)'[r]+2 r A03[r] B01'[r] (ϕ01)'[r]+3 r A02[r] B02'[r] (ϕ01)'[r]+2 r A03[r] B02'[r] (ϕ01)'[r]+2 r A02[r] B03'[r] (ϕ01)'[r]+2 r A03[r] B03'[r] (ϕ01)'[r]+16 A02[r] B03[r] (ϕ02)'[r]+16 A03[r] B03[r] (ϕ02)'[r]+2 r B03[r] A01'[r] (ϕ02)'[r]+2 r B03[r] A02'[r] (ϕ02)'[r]+2 r B03[r] A03'[r] (ϕ02)'[r]+2 r A02[r] B01'[r] (ϕ02)'[r]+2 r A03[r] B01'[r] (ϕ02)'[r]+3 r A02[r] B02'[r] (ϕ02)'[r]+2 r A03[r] B02'[r] (ϕ02)'[r]+2 r A02[r] B03'[r] (ϕ02)'[r]+2 r A03[r] B03'[r] (ϕ02)'[r]+8 r A02[r] B03[r] (ϕ01)''[r]+8 r A03[r] B03[r] (ϕ01)''[r]+8 r A02[r] B03[r] (ϕ02)''[r]+8 r A03[r] B03[r] (ϕ02)''[r])+B03[r] (120 r (ω^2) ϕ02[r]+r B03[r] (4 A01'[r] (ϕ01)'[r]+3 A03'[r] (ϕ02)'[r]+2 A02'[r] ((ϕ01)'[r]+(ϕ02)'[r]))+4 A02[r] (-2 r (μ^2) ϕ01[r]-2 r (μ^2) ϕ02[r]+4 B03[r] (ϕ01)'[r]+r B01'[r] (ϕ01)'[r]+r B02'[r] (ϕ01)'[r]+r B03'[r] (ϕ01)'[r]+4 B03[r] (ϕ02)'[r]+r B01'[r] (ϕ02)'[r]+r B02'[r] (ϕ02)'[r]+r B03'[r] (ϕ02)'[r]+2 r B03[r] (ϕ01)''[r]+2 r B03[r] (ϕ02)''[r])+A03[r] (-16 r (μ^2) ϕ01[r]+8 r B01'[r] (ϕ01)'[r]+4 r B02'[r] (ϕ01)'[r]+24 B03[r] (ϕ02)'[r]+4 r B02'[r] (ϕ02)'[r]+6 r B03'[r] (ϕ02)'[r]+12 r B03[r] (ϕ02)''[r]))+4 ((B01[r])^2) (2 r A01'[r] (ϕ01)'[r]+r A02'[r] (ϕ01)'[r]+8 A03[r] (ϕ02)'[r]+r A02'[r] (ϕ02)'[r]+r A03'[r] (ϕ02)'[r]+4 r A03[r] (ϕ02)''[r]+4 A02[r] (2 (ϕ01)'[r]+2 (ϕ02)'[r]+r ((ϕ01)''[r]+(ϕ02)''[r])))+((B02[r])^2) (r (3 A02'[r]+2 A03'[r]) ((ϕ01)'[r]+(ϕ02)'[r])+2 r A01'[r] (2 (ϕ01)'[r]+(ϕ02)'[r])+12 A02[r] (2 (ϕ01)'[r]+2 (ϕ02)'[r]+r ((ϕ01)''[r]+(ϕ02)''[r]))+8 A03[r] (2 (ϕ01)'[r]+2 (ϕ02)'[r]+r ((ϕ01)''[r]+(ϕ02)''[r])))+4 B01[r] (4 r ((ω^2)-(μ^2) A02[r]) ϕ01[r]-4 r (μ^2) A03[r] ϕ02[r]+8 A03[r] B02[r] (ϕ01)'[r]+16 A03[r] B03[r] (ϕ01)'[r]+2 r B02[r] A01'[r] (ϕ01)'[r]+2 r B02[r] A02'[r] (ϕ01)'[r]+r B03[r] A02'[r] (ϕ01)'[r]+r B02[r] A03'[r] (ϕ01)'[r]+2 r B03[r] A03'[r] (ϕ01)'[r]+r A03[r] B02'[r] (ϕ01)'[r]+2 r A03[r] B03'[r] (ϕ01)'[r]+8 A03[r] B02[r] (ϕ02)'[r]+2 r B02[r] A01'[r] (ϕ02)'[r]+2 r B03[r] A01'[r] (ϕ02)'[r]+r B02[r] A02'[r] (ϕ02)'[r]+r B03[r] A02'[r] (ϕ02)'[r]+r B02[r] A03'[r] (ϕ02)'[r]+2 r A03[r] B01'[r] (ϕ02)'[r]+r A03[r] B02'[r] (ϕ02)'[r]+4 r A03[r] B02[r] (ϕ01)''[r]+8 r A03[r] B03[r] (ϕ01)''[r]+4 r A03[r] B02[r] (ϕ02)''[r]+A02[r] (-4 r (μ^2) ϕ02[r]+8 B03[r] (ϕ01)'[r]+2 r B01'[r] (ϕ01)'[r]+2 r B02'[r] (ϕ01)'[r]+r B03'[r] (ϕ01)'[r]+8 B03[r] (ϕ02)'[r]+2 r B01'[r] (ϕ02)'[r]+r B02'[r] (ϕ02)'[r]+r B03'[r] (ϕ02)'[r]+4 r B03[r] (ϕ01)''[r]+4 r B03[r] (ϕ02)''[r]+4 B02[r] (4 (ϕ01)'[r]+2 (ϕ02)'[r]+2 r (ϕ01)''[r]+r (ϕ02)''[r]))))+A03[r] (r (2 B02[r] (8 (ω^2) ϕ01[r]+48 (ω^2) ϕ02[r]+4 B01[r] A01'[r] (ϕ01)'[r]+4 B03[r] A01'[r] (ϕ01)'[r]+4 B01[r] A02'[r] (ϕ01)'[r]+4 B03[r] A02'[r] (ϕ01)'[r]+4 B01[r] A03'[r] (ϕ01)'[r]+3 B03[r] A03'[r] (ϕ01)'[r]+4 B01[r] A01'[r] (ϕ02)'[r]+4 B03[r] A01'[r] (ϕ02)'[r]+4 B01[r] A02'[r] (ϕ02)'[r]+3 B03[r] A02'[r] (ϕ02)'[r]+4 B01[r] A03'[r] (ϕ02)'[r]+3 B03[r] A03'[r] (ϕ02)'[r])+((B02[r])^2) (4 A01'[r] ((ϕ01)'[r]+(ϕ02)'[r])+4 A02'[r] ((ϕ01)'[r]+(ϕ02)'[r])+A03'[r] (4 (ϕ01)'[r]+3 (ϕ02)'[r]))+3 ((B03[r])^2) (A02'[r] ((ϕ01)'[r]+(ϕ02)'[r])+2 (A03'[r] (ϕ01)'[r]+A01'[r] (ϕ02)'[r]))+4 B03[r] (4 (ω^2) ϕ01[r]+B01[r] (4 A01'[r] (ϕ01)'[r]+3 A03'[r] (ϕ02)'[r]+2 A02'[r] ((ϕ01)'[r]+(ϕ02)'[r])))+4 B01[r] (60 (ω^2) ϕ02[r]+B01[r] (A02'[r] ((ϕ01)'[r]+(ϕ02)'[r])+2 (A03'[r] (ϕ01)'[r]+A01'[r] (ϕ02)'[r]))))+A03[r] (16 ((B01[r])^2) (2 (ϕ01)'[r]+r (ϕ01)''[r])+3 B03[r] (-4 r (μ^2) ϕ02[r]+r B02'[r] ((ϕ01)'[r]+(ϕ02)'[r])+2 r (B03'[r] (ϕ01)'[r]+B01'[r] (ϕ02)'[r])+4 B03[r] (2 (ϕ01)'[r]+r (ϕ01)''[r]))+2 ((B02[r])^2) (8 (ϕ01)'[r]+6 (ϕ02)'[r]+4 r (ϕ01)''[r]+3 r (ϕ02)''[r])+B02[r] (-8 r (μ^2) ϕ01[r]-8 r (μ^2) ϕ02[r]+24 B03[r] (ϕ01)'[r]+4 r B01'[r] (ϕ01)'[r]+4 r B02'[r] (ϕ01)'[r]+3 r B03'[r] (ϕ01)'[r]+24 B03[r] (ϕ02)'[r]+4 r B01'[r] (ϕ02)'[r]+3 r B02'[r] (ϕ02)'[r]+3 r B03'[r] (ϕ02)'[r]+12 r B03[r] (ϕ01)''[r]+12 r B03[r] (ϕ02)''[r])+2 B01[r] (-8 r (μ^2) ϕ01[r]+4 r B01'[r] (ϕ01)'[r]+2 r B02'[r] (ϕ01)'[r]+24 B03[r] (ϕ02)'[r]+2 r B02'[r] (ϕ02)'[r]+3 r B03'[r] (ϕ02)'[r]+12 r B03[r] (ϕ02)''[r]+8 B02[r] (2 (ϕ01)'[r]+2 (ϕ02)'[r]+r ((ϕ01)''[r]+(ϕ02)''[r])))))),(1/(32 r))(((A02[r])^2) (8 ((B01[r])^2) (2 (ϕ01)'[r]+4 (ϕ02)'[r]+r (ϕ01)''[r]+2 r (ϕ02)''[r])+4 ((B02[r])^2) (4 (ϕ01)'[r]+6 (ϕ02)'[r]+2 r (ϕ01)''[r]+3 r (ϕ02)''[r])+B03[r] (-8 r (μ^2) ϕ01[r]-8 r (μ^2) ϕ02[r]+12 B03[r] (ϕ01)'[r]+4 r B01'[r] (ϕ01)'[r]+4 r B02'[r] (ϕ01)'[r]+3 r B03'[r] (ϕ01)'[r]+16 B03[r] (ϕ02)'[r]+4 r B01'[r] (ϕ02)'[r]+3 r B02'[r] (ϕ02)'[r]+4 r B03'[r] (ϕ02)'[r]+6 r B03[r] (ϕ01)''[r]+8 r B03[r] (ϕ02)''[r])+2 B01[r] (-4 r (μ^2) ϕ01[r]-8 r (μ^2) ϕ02[r]+24 B02[r] (ϕ01)'[r]+16 B03[r] (ϕ01)'[r]+2 r B01'[r] (ϕ01)'[r]+3 r B02'[r] (ϕ01)'[r]+2 r B03'[r] (ϕ01)'[r]+8 B02[r] (ϕ02)'[r]+16 B03[r] (ϕ02)'[r]+4 r B01'[r] (ϕ02)'[r]+r B02'[r] (ϕ02)'[r]+2 r B03'[r] (ϕ02)'[r]+12 r B02[r] (ϕ01)''[r]+8 r B03[r] (ϕ01)''[r]+4 r B02[r] (ϕ02)''[r]+8 r B03[r] (ϕ02)''[r])+B02[r] (-12 r (μ^2) ϕ01[r]-4 r (μ^2) ϕ02[r]+32 B03[r] (ϕ01)'[r]+6 r B01'[r] (ϕ01)'[r]+4 r B02'[r] (ϕ01)'[r]+4 r B03'[r] (ϕ01)'[r]+24 B03[r] (ϕ02)'[r]+2 r B01'[r] (ϕ02)'[r]+6 r B02'[r] (ϕ02)'[r]+3 r B03'[r] (ϕ02)'[r]+16 r B03[r] (ϕ01)''[r]+12 r B03[r] (ϕ02)''[r]))+4 ((A01[r])^2) (B03[r] (-4 r (μ^2) ϕ01[r]+r (2 B01'[r] (ϕ01)'[r]+2 B03'[r] (ϕ02)'[r]+B02'[r] ((ϕ01)'[r]+(ϕ02)'[r]))+8 B01[r] (2 (ϕ01)'[r]+r (ϕ01)''[r]))+4 ((B03[r])^2) (2 (ϕ02)'[r]+r (ϕ02)''[r])+2 B01[r] (-4 r (μ^2) ϕ02[r]+r B02'[r] (ϕ01)'[r]+r B03'[r] (ϕ01)'[r]+8 B01[r] (ϕ02)'[r]+2 r B01'[r] (ϕ02)'[r]+4 r B01[r] (ϕ02)''[r])+B02[r] (-4 r (μ^2) ϕ01[r]+8 B03[r] (ϕ01)'[r]+2 r B01'[r] (ϕ01)'[r]+r B02'[r] (ϕ01)'[r]+r B03'[r] (ϕ01)'[r]+8 B03[r] (ϕ02)'[r]+2 r B02'[r] (ϕ02)'[r]+r B03'[r] (ϕ02)'[r]+4 r B03[r] (ϕ01)''[r]+8 B01[r] (2 (ϕ01)'[r]+r (ϕ01)''[r])+4 r B03[r] (ϕ02)''[r])+2 ((B02[r])^2) (2 (ϕ01)'[r]+4 (ϕ02)'[r]+r ((ϕ01)''[r]+2 (ϕ02)''[r])))+A02[r] (r (2 B02[r] (12 (ω^2) ϕ01[r]+72 (ω^2) ϕ02[r]+4 B01[r] A01'[r] (ϕ01)'[r]+4 B03[r] A01'[r] (ϕ01)'[r]+6 B01[r] A02'[r] (ϕ01)'[r]+4 B03[r] A02'[r] (ϕ01)'[r]+4 B01[r] A03'[r] (ϕ01)'[r]+3 B03[r] A03'[r] (ϕ01)'[r]+8 B01[r] A01'[r] (ϕ02)'[r]+4 B03[r] A01'[r] (ϕ02)'[r]+2 B01[r] A02'[r] (ϕ02)'[r]+3 B03[r] A02'[r] (ϕ02)'[r]+4 B01[r] A03'[r] (ϕ02)'[r]+4 B03[r] A03'[r] (ϕ02)'[r])+4 ((B01[r])^2) (2 A01'[r] (ϕ01)'[r]+A03'[r] ((ϕ01)'[r]+(ϕ02)'[r])+A02'[r] ((ϕ01)'[r]+2 (ϕ02)'[r]))+((B03[r])^2) (3 A03'[r] ((ϕ01)'[r]+(ϕ02)'[r])+2 A01'[r] (2 (ϕ01)'[r]+(ϕ02)'[r])+A02'[r] (3 (ϕ01)'[r]+4 (ϕ02)'[r]))+((B02[r])^2) (2 A01'[r] (3 (ϕ01)'[r]+(ϕ02)'[r])+A03'[r] (4 (ϕ01)'[r]+3 (ϕ02)'[r])+A02'[r] (4 (ϕ01)'[r]+6 (ϕ02)'[r]))+4 B03[r] (36 (ω^2) ϕ02[r]+B01[r] (2 A01'[r] ((ϕ01)'[r]+(ϕ02)'[r])+2 A02'[r] ((ϕ01)'[r]+(ϕ02)'[r])+A03'[r] (2 (ϕ01)'[r]+(ϕ02)'[r]))))+2 A03[r] (2 ((B02[r])^2) (8 (ϕ01)'[r]+6 (ϕ02)'[r]+4 r (ϕ01)''[r]+3 r (ϕ02)''[r])+2 B01[r] (-4 r (μ^2) ϕ01[r]-4 r (μ^2) ϕ02[r]+16 B02[r] (ϕ01)'[r]+16 B03[r] (ϕ01)'[r]+2 r B01'[r] (ϕ01)'[r]+2 r B02'[r] (ϕ01)'[r]+2 r B03'[r] (ϕ01)'[r]+16 B02[r] (ϕ02)'[r]+8 B03[r] (ϕ02)'[r]+2 r B01'[r] (ϕ02)'[r]+2 r B02'[r] (ϕ02)'[r]+r B03'[r] (ϕ02)'[r]+8 r B02[r] (ϕ01)''[r]+8 r B03[r] (ϕ01)''[r]+8 r B02[r] (ϕ02)''[r]+4 r B03[r] (ϕ02)''[r])+B03[r] (-8 r (μ^2) ϕ01[r]-4 r (μ^2) ϕ02[r]+12 B03[r] (ϕ01)'[r]+4 r B01'[r] (ϕ01)'[r]+3 r B02'[r] (ϕ01)'[r]+3 r B03'[r] (ϕ01)'[r]+12 B03[r] (ϕ02)'[r]+2 r B01'[r] (ϕ02)'[r]+4 r B02'[r] (ϕ02)'[r]+3 r B03'[r] (ϕ02)'[r]+6 r B03[r] (ϕ01)''[r]+6 r B03[r] (ϕ02)''[r])+B02[r] (-8 r (μ^2) ϕ01[r]-8 r (μ^2) ϕ02[r]+24 B03[r] (ϕ01)'[r]+4 r B01'[r] (ϕ01)'[r]+4 r B02'[r] (ϕ01)'[r]+3 r B03'[r] (ϕ01)'[r]+32 B03[r] (ϕ02)'[r]+4 r B01'[r] (ϕ02)'[r]+3 r B02'[r] (ϕ02)'[r]+4 r B03'[r] (ϕ02)'[r]+12 r B03[r] (ϕ01)''[r]+16 r B03[r] (ϕ02)''[r])+8 ((B01[r])^2) (2 (ϕ01)'[r]+2 (ϕ02)'[r]+r ((ϕ01)''[r]+(ϕ02)''[r]))))+A03[r] (r (((B02[r])^2) (4 A01'[r] ((ϕ01)'[r]+(ϕ02)'[r])+A02'[r] (4 (ϕ01)'[r]+3 (ϕ02)'[r])+A03'[r] (3 (ϕ01)'[r]+4 (ϕ02)'[r]))+3 B03[r] (48 (ω^2) ϕ02[r]+B03[r] (2 A01'[r] (ϕ01)'[r]+2 A03'[r] (ϕ02)'[r]+A02'[r] ((ϕ01)'[r]+(ϕ02)'[r])))+2 B02[r] (72 (ω^2) ϕ02[r]+B03[r] (3 A03'[r] ((ϕ01)'[r]+(ϕ02)'[r])+2 A01'[r] (2 (ϕ01)'[r]+(ϕ02)'[r])+A02'[r] (3 (ϕ01)'[r]+4 (ϕ02)'[r]))))+4 ((B01[r])^2) (2 r A01'[r] (ϕ01)'[r]+8 A03[r] (ϕ02)'[r]+2 r A03'[r] (ϕ02)'[r]+r A02'[r] ((ϕ01)'[r]+(ϕ02)'[r])+4 r A03[r] (ϕ02)''[r])+A03[r] (2 ((B02[r])^2) (6 (ϕ01)'[r]+8 (ϕ02)'[r]+3 r (ϕ01)''[r]+4 r (ϕ02)''[r])+3 B03[r] (-4 r (μ^2) ϕ01[r]+2 r B01'[r] (ϕ01)'[r]+r B02'[r] (ϕ01)'[r]+8 B03[r] (ϕ02)'[r]+r B02'[r] (ϕ02)'[r]+2 r B03'[r] (ϕ02)'[r]+4 r B03[r] (ϕ02)''[r])+B02[r] (-8 r (μ^2) ϕ01[r]-4 r (μ^2) ϕ02[r]+24 B03[r] (ϕ01)'[r]+4 r B01'[r] (ϕ01)'[r]+3 r B02'[r] (ϕ01)'[r]+3 r B03'[r] (ϕ01)'[r]+24 B03[r] (ϕ02)'[r]+2 r B01'[r] (ϕ02)'[r]+4 r B02'[r] (ϕ02)'[r]+3 r B03'[r] (ϕ02)'[r]+12 r B03[r] (ϕ01)''[r]+12 r B03[r] (ϕ02)''[r]))+2 B01[r] (24 r (ω^2) ϕ01[r]+2 r (B03[r] (3 A03'[r] (ϕ01)'[r]+4 A01'[r] (ϕ02)'[r]+A02'[r] (2 (ϕ01)'[r]+(ϕ02)'[r]))+B02[r] (2 A01'[r] ((ϕ01)'[r]+(ϕ02)'[r])+2 A02'[r] ((ϕ01)'[r]+(ϕ02)'[r])+A03'[r] (2 (ϕ01)'[r]+(ϕ02)'[r])))+A03[r] (-8 r (μ^2) ϕ02[r]+24 B03[r] (ϕ01)'[r]+2 r B02'[r] (ϕ01)'[r]+3 r B03'[r] (ϕ01)'[r]+4 r B01'[r] (ϕ02)'[r]+r B02'[r] (ϕ02)'[r]+12 r B03[r] (ϕ01)''[r]+4 B02[r] (4 (ϕ01)'[r]+2 (ϕ02)'[r]+2 r (ϕ01)''[r]+r (ϕ02)''[r]))))+2 A01[r] (r (4 ((B01[r])^2) (A02'[r] (ϕ01)'[r]+A03'[r] (ϕ01)'[r]+2 A01'[r] (ϕ02)'[r])+((B03[r])^2) (3 A03'[r] (ϕ01)'[r]+4 A01'[r] (ϕ02)'[r]+A02'[r] (2 (ϕ01)'[r]+(ϕ02)'[r]))+((B02[r])^2) (2 A03'[r] ((ϕ01)'[r]+(ϕ02)'[r])+A02'[r] (3 (ϕ01)'[r]+(ϕ02)'[r])+2 A01'[r] ((ϕ01)'[r]+2 (ϕ02)'[r]))+2 B03[r] (12 (ω^2) ϕ01[r]+2 B01[r] (2 A01'[r] (ϕ01)'[r]+2 A03'[r] (ϕ02)'[r]+A02'[r] ((ϕ01)'[r]+(ϕ02)'[r]))+B02[r] (2 A01'[r] ((ϕ01)'[r]+(ϕ02)'[r])+2 A02'[r] ((ϕ01)'[r]+(ϕ02)'[r])+A03'[r] (2 (ϕ01)'[r]+(ϕ02)'[r])))+4 B01[r] (36 (ω^2) ϕ02[r]+B02[r] (2 A01'[r] (ϕ01)'[r]+A03'[r] ((ϕ01)'[r]+(ϕ02)'[r])+A02'[r] ((ϕ01)'[r]+2 (ϕ02)'[r]))))+2 A02[r] (8 ((B01[r])^2) (2 (ϕ01)'[r]+r (ϕ01)''[r])+B03[r] (-4 r (μ^2) ϕ01[r]-4 r (μ^2) ϕ02[r]+8 B03[r] (ϕ01)'[r]+2 r B01'[r] (ϕ01)'[r]+2 r B02'[r] (ϕ01)'[r]+2 r B03'[r] (ϕ01)'[r]+4 B03[r] (ϕ02)'[r]+2 r B01'[r] (ϕ02)'[r]+2 r B02'[r] (ϕ02)'[r]+r B03'[r] (ϕ02)'[r]+4 r B03[r] (ϕ01)''[r]+2 r B03[r] (ϕ02)''[r])+B02[r] (-4 r (μ^2) ϕ01[r]-8 r (μ^2) ϕ02[r]+16 B03[r] (ϕ01)'[r]+2 r B01'[r] (ϕ01)'[r]+3 r B02'[r] (ϕ01)'[r]+2 r B03'[r] (ϕ01)'[r]+16 B03[r] (ϕ02)'[r]+4 r B01'[r] (ϕ02)'[r]+r B02'[r] (ϕ02)'[r]+2 r B03'[r] (ϕ02)'[r]+8 r B03[r] (ϕ01)''[r]+8 r B03[r] (ϕ02)''[r])+2 ((B02[r])^2) (6 (ϕ01)'[r]+2 (ϕ02)'[r]+r (3 (ϕ01)''[r]+(ϕ02)''[r]))+2 B01[r] (-4 r (μ^2) ϕ01[r]+8 B03[r] (ϕ01)'[r]+2 r B01'[r] (ϕ01)'[r]+r B02'[r] (ϕ01)'[r]+r B03'[r] (ϕ01)'[r]+8 B03[r] (ϕ02)'[r]+2 r B02'[r] (ϕ02)'[r]+r B03'[r] (ϕ02)'[r]+4 r B03[r] (ϕ01)''[r]+4 r B03[r] (ϕ02)''[r]+4 B02[r] (2 (ϕ01)'[r]+4 (ϕ02)'[r]+r (ϕ01)''[r]+2 r (ϕ02)''[r])))+2 A03[r] (8 ((B01[r])^2) (2 (ϕ01)'[r]+r (ϕ01)''[r])+B03[r] (-8 r (μ^2) ϕ02[r]+r (3 B03'[r] (ϕ01)'[r]+4 B01'[r] (ϕ02)'[r]+B02'[r] (2 (ϕ01)'[r]+(ϕ02)'[r]))+6 B03[r] (2 (ϕ01)'[r]+r (ϕ01)''[r]))+B02[r] (-4 r (μ^2) ϕ01[r]-4 r (μ^2) ϕ02[r]+16 B03[r] (ϕ01)'[r]+2 r B01'[r] (ϕ01)'[r]+2 r B02'[r] (ϕ01)'[r]+2 r B03'[r] (ϕ01)'[r]+8 B03[r] (ϕ02)'[r]+2 r B01'[r] (ϕ02)'[r]+2 r B02'[r] (ϕ02)'[r]+r B03'[r] (ϕ02)'[r]+8 r B03[r] (ϕ01)''[r]+4 r B03[r] (ϕ02)''[r])+4 ((B02[r])^2) (2 (ϕ01)'[r]+2 (ϕ02)'[r]+r ((ϕ01)''[r]+(ϕ02)''[r]))+2 B01[r] (-4 r (μ^2) ϕ01[r]+2 r B01'[r] (ϕ01)'[r]+r B02'[r] (ϕ01)'[r]+16 B03[r] (ϕ02)'[r]+r B02'[r] (ϕ02)'[r]+2 r B03'[r] (ϕ02)'[r]+8 r B03[r] (ϕ02)''[r]+4 B02[r] (2 (ϕ01)'[r]+2 (ϕ02)'[r]+r ((ϕ01)''[r]+(ϕ02)''[r]))))))};

eqsbkg2=Solve[eqsbkg==0,{A01'[r],A02'[r],A03'[r],B01'[r],B02'[r],B03'[r],ϕ01''[r],ϕ02''[r]}][[1]]/.Rule->Subtract;

μ=1;
κ=1;

Boundary conditions

bc={A01[ror]==1,A02[ror]==(10^(-12)),A03[ror]==(10^(-12)),B01[ror]==1,B02[ror]==(10^(-12)),B03[ror]==(10^(-12)),(ϕ01)'[ror]==(10^(-12)),ϕ01[ror]==ϕc,(ϕ02)'[ror]==(10^(-12)),ϕ02[ror]==(10^(-12))};
bc2={A01[ror]==A00,A02[ror]==(10^(-12)),A03[ror]==(10^(-12)),B01[ror]==1,B02[ror]==(10^(-12)),B03[ror]==(10^(-12)),(ϕ01)'[ror]==(10^(-12)),ϕ01[ror]==ϕc,(ϕ02)'[ror]==(10^(-12)),ϕ02[ror]==(10^(-12))};

INTB[ωs0_,ϕ0c_,ror0_,rinf0_]:=(
(*Defining rules for the parameters*)
param={ϕc->ϕ0c,ror->ror0,ω->ωs0};
(*Setting the equations and the boundary conditions*)
EQS=Thread[eqsbkg=={0,0,0,0,0,0,0,0}]//.param;
BOU=bc//.param;
(*Solving the equations*)
SOL=NDSolve[Join[EQS,BOU],{A01,A02,A03,B01,B02,B03,ϕ01,ϕ02},{r,ror0,rinf0},{AccuracyGoal->13,PrecisionGoal->12,MaxSteps->∞,Method->"StiffnessSwitching",InterpolationOrder->All}][[1]]//Quiet;
RINF=Min[rinf0,(ϕ01/.SOL)[[1,1,2]]];
(*We can additionally provide the quantities that will aid to describe the spacetime in Schwarzschild coordinates*)
CC=(((B01[r]))/((A01[r])))/.SOL/.r->RINF;
Ar=CC (A01[r]+A02[r]Cos[2ωs0 t]+A03[r]Cos[4ωs0 t])/.SOL;
Br=(B01[r]+B02[r]Cos[2ωs0 t]+B03[r]Cos[4ωs0 t])/.SOL;
ϕ0r=(ϕ01[r]Cos[ωs0 t]+ϕ02[r]Cos[3ωs0 t])/.SOL;
ωss=ωs0(Sqrt[CC]);
Mt=-(1/2) r (-1+B01[r])/.SOL/.r->RINF;
mr=-(1/2) r (-1+B01[r])/.SOL;
)

INTB2[ωs0_,A0c_,ϕ0c_,ror0_,rinf0_]:=(
(*Defining rules for the parameters*)
param={ϕc->ϕ0c,A00->A0c,ror->ror0,ω->ωs0};
(*Setting the equations and the boundary conditions*)
EQS=Thread[eqsbkg=={0,0,0,0,0,0,0,0}]//.param;
BOU=bc2//.param;
(*Solving the equations*)
SOL=NDSolve[Join[EQS,BOU],{A01,A02,A03,B01,B02,B03,ϕ01,ϕ02},{r,ror0,rinf0},{AccuracyGoal->13,PrecisionGoal->12,MaxSteps->∞,Method->"StiffnessSwitching",InterpolationOrder->All}][[1]]//Quiet;
RINF=Min[rinf0,(ϕ01/.SOL)[[1,1,2]]];
(*We can additionally provide the quantities that will aid to describe the spacetime in Schwarzschild coordinates*)
CC=(((B01[r]))/((A01[r])))/.SOL/.r->RINF;
ωss=ωs0(Sqrt[CC]);
Ar=CC(A01[r]+A02[r]Cos[2ωs0 t]+A03[r]Cos[4ωs0 t])/.SOL;
Br=(B01[r]+B02[r]Cos[2ωs0 t]+B03[r]Cos[4ωs0 t])/.SOL;
ϕ0r=(ϕ01[r]Cos[ωss t]+ϕ02[r]Cos[3ωss t])/.SOL;
Mt=-(1/2) r (-1+B01[r])/.SOL/.r->RINF;
mr=-(1/2) r (-1+B01[r])/.SOL;
)

INTB[.05,0.1,0.001,10]//Timing

INTB2[.05,1,0.1,0.001,10]//Timing

shoot[ωs0_?NumberQ,ϕ0c_,ror0_,rinf0_]:=(
INTB[ωs0,ϕ0c,ror0,rinf0];
ϕ01[r]/.SOL/.r->RINF
)

shoot2[ωs0_?NumberQ,A0c_?NumberQ,ϕ0c_,ror0_,rinf0_]:=(
INTB2[ωs0,A0c,ϕ0c,ror0,rinf0];
{ϕ01[r],A02[r]}/.SOL/.r->RINF
)

shoot2[0.8428305588844424`,0.48055640747723877`,0.73617,0.001,25]

LogPlot[1-B02[r]/.SOL//Abs,{r,0.001,RINF},PlotRange->All]

Monitor[FindRoot[shoot2[xx,yy,0.73617,0.001,15],{xx,0.84},{yy,0.4}],{xx,yy}]
LogPlot[ϕ01[r]/.SOL//Abs,{r,0.001,RINF},PlotRange->All]

1.2158316764661223`(Sqrt[CC/.t->0])

Monitor[FindRoot[shoot[xx,0.73617,0.001,25],{xx,1.2158316764661223`}],xx]
LogPlot[ϕ01[r]/.SOL//Abs,{r,0.001,RINF},PlotRange->All]

ωss

Mt

(1/(Br Ar))/.t->0

Series[(1/(B01[r]+ϵ B02[r])),{ϵ,0,1}]

Ar

Series[(1/(CC(A01[r]+ϵ A02[r]Cos[2ωs0 t]) CC(B01[r]+ϵ B02[r]Cos[2ωs0 t]))),{ϵ,0,1}]

Plot[{(( ((A02[r] B01[r]+A01[r] B02[r])) )/(((A01[r])^2) ((B01[r])^2)))/.SOL},{r,0.01,RINF},PlotRange->All,Frame->True]

Animate[Plot[{Ar/.t->t0,Br/.t->t0},{r,0.01,RINF},PlotRange->{All,{0,1}}] ,{t0,0,10},AnimationRunning->False]

Plot[{ϕ01[r]/.SOL,ϕ02[r]/.SOL},{r,0.01,RINF},PlotRange->All]

Animate[Plot[Ar,{r}]]

CC/.r->RINF

ωss/.r->RINF

Trying to grasp the right value of ω

sx=Monitor[FindRoot[shoot[xx,0.01,0.001,40],{xx,1.004161359982149103`}],xx]
Plot[ϕ01[r]/.SOL,{r,0.001,RINF}]
Mt

Rdi=r/.FindRoot[-(1/2) (-1+B01[r]) r==0.99Mt/.SOL,{r,9}]

{ϕc,Mt,Rdi,ω}//.param

pts={{0.01`,0.12119206514566772`,79.47843704499802`,1.002303923987936`},{0.1`,0.3569252366663833`,25.004699895490738`,1.0236656851247683`};

ListPlot[pts[[All,{3,2}]],Joined->True]

r/.FindRoot[-(1/2) (-1+B01[r]) r==0.99Mt/.SOL,{r,90}]

ϕii=0.01;
ϕfi=.14;
pts=800;
ωguesi=1.002303923987936`;
Rdi=80;

Monitor[Do[
(*Defining the steps*)
ϕit=Log[ϕii]+((Log[ϕfi]-Log[ϕii])/pts)i;
(*Defining the findroot for the solution and the radius*)
sww=FindRoot[shoot[xx,Exp[ϕit],0.001,1.8 Rdi],{xx,ωguesi}];
Rdi=r/.FindRoot[-(1/2) (-1+B01[r]) r==0.99Mt/.SOL,{r,Rdi}]//Quiet;
ωguesi=xx/.sww;
(*collecting the quatities*)
ϕi[i]=Exp[ϕit];
Mti[i]=Mt;
Rd[i]=Rdi;
ωii[i]=ωguesi;
(*Tracking the points and plotting within the monitor*)
tabϕM=Table[{ϕi[j],Mti[j]},{j,0,i}];
mplot=ListPlot[tabϕM,PlotRange->All,Joined->True,Frame->True];
tabRM=Table[{Rd[j],Mti[j]},{j,0,i}];
mplot2=ListPlot[tabRM,PlotRange->All,Joined->True,Frame->True];
,{i,0,pts}],{i,mplot,mplot2}]

tabRM//Length

texp[[1]]

Exporting the points

texp=Table[{ϕi[j],Mti[j],Rd[j],ωii[j]},{j,0,800,5}];

texp[[-1]]

SetDirectory[NotebookDirectory[]]

Export["solutions/osci_parameters_phi_M_R_N_w_2.dat",texp]

Testing

SetDirectory[NotebookDirectory[]]

t1=Import["solutions/osci_parameters_phi_M_R_N_w.dat"];

tx=t1[[All,{1,2}]];
tx2=t1[[All,{1,-1}]];

ListPlot[tx,PlotRange->All]

FindMaximum[Interpolation[tx,xx],{xx,0.5}]

Interpolation[tx2,0.73617]

Monitor[FindRoot[shoot[xx,0.73617,0.001,25],{xx,1.2158316764661223`}],xx]
LogPlot[ϕ01[r]/.SOL//Abs,{r,0.001,RINF},PlotRange->All]

Plot[{-B01[r]-B02[r]-B03[r]+1/.SOL,-B02[r]/.SOL,-B01[r]+1/.SOL},{r,0.01,RINF}]

mr=-(1/2) r (-1+B01[r])/.SOL

Plot[{A01[r]/.SOL,A02[r]/.SOL,A03[r]/.SOL},{r,0.01,RINF}]

rr=r/.FindRoot[mr==0.99Mt,{r,8}]

((7.865420476015143`)/0.633)

D[mr,r]

tt=Piecewise[{{D[mr,r],r<rr},{0,r>rr}}]

Plot[tt,{r,0.001,13}]

(rr/Mt)

pden=ContourPlot[tt/.r->Mt(Sqrt[(x^2)+(y^2)]),{x,-20,20},{y,-20,20},PlotPoints->50,PlotRange->{{-20,13},{-13,13},All},ColorFunction->(Blend[{
{0, RGBColor[1, 0, 0, 0]},
{0.1,RGBColor[1, 0, 0, 0.1]},
{0.2,RGBColor[1, 0, 0, 0.2]},
{0.3,RGBColor[1, 0, 0, 0.3]},
{0.4,RGBColor[1, 0, 0, 0.4]},
{0.5,RGBColor[1, 0, 0, 0.5]},
{0.6,RGBColor[1, 0, 0, 0.6]},
{0.7,RGBColor[1, 0, 0, 0.7]},
{0.8,RGBColor[1, 0, 0, 0.8]},
{0.9,RGBColor[1, 0, 0, 0.9]},
{1,RGBColor[1, 0, 0, 1]}}, #1] &),ContourLines->False,Contours->5]

Comparing oscillatons, boson and vector boson stars

Quit[]

SetDirectory[NotebookDirectory[]];

tsbs=Import["solutions/miniBS_parameters_phi_M_R_N_w.dat"];
tvbs=Import["solutions/vectorminiBS_parameters_fc_M_R_N_w.dat"];
tvos=Import["solutions/osci_parameters_phi_M_R_N_w.dat"];

ListLogLogPlot[{tsbs[[All,{3,2}]],tvos[[All,{3,2}]],tvbs[[All,{3,2}]]},Joined->True,PlotRange->All]

tvbs[[-1]]

tvos[[All,{3,2}]][[1]]

tvos[[All,{3,4}]][[1]]

ListPlot[{tvbs[[All,{-2,2}]],tvbs[[All,{-2,4}]]},Joined->True,PlotRange->All]

ListPlot[{tsbs[[All,{-2,2}]],tsbs[[All,{-2,4}]]},Joined->True,PlotRange->All]