Equations of motion for charged scalar and vector boson stars

Here we describe the equations to obtain spherically symmetric solutions of the system composed by a complex charged scalar and vector. We shall consider all terms combined to obtain the differential equations, using a control parameter to turn off eventually. 
We also add the perturbations, to illustrate how the procedure works to find the quasinormal modes of such configurations.

Let us first clear all the possible variables

Quit[];

LaunchKernels[]

and we define the coordinate system to be used (we shall use cos=χ to better perform the algebraic computations)

$Assumptions={r≥0,1≥χ≥-1,Element[{𝓂,ℓ},Integers]};
x={t,r,χ,ϕ};

Metric, scalar field and vector field

Metric for spherically symmetric spacetime, plus the axial and polar perturbations

Let us define the metric to be used in the equations of motion. The metric will be represented by a matrix

gdd={{g00,g01,g02,g03},{g01,g11,g12,g13},{g02,g12,g22,g23},{g03,g13,g23,g33}};

The metric components are

(*Polar perturbation*)
h00p=H0[t,r]Y[χ];
h01p=H1[t,r]Y[χ];
h11p=H2[t,r]Y[χ];
h22p=(r^2)Kp[t,r]Y[χ]/(1-χ^2);
h33p=(r^2)Kp[t,r]Y[χ](1-χ^2);
(*Axial perturbation*)
h02a=0;
h03a=-h0[t,r]Y'[χ](1-χ^2);
h12a=0;
h13a=-h1[t,r]Y'[χ] (1-χ^2);
(*defining the metric components*)
g00=-A[t,r](1+ϵp h00p); 
g01=ϵp h01p;
g02=ϵa h02a;
g03=ϵa h03a;
g11=((B[t,r])^(-1))(1+ϵp h11p);
g12=ϵa h12a;
g13=ϵa h13a;
g22=((r^2)/(1-(χ^2)))+ϵp h22p;
g23=0;
g33=(r^2)(1-(χ^2))+ϵp h33p;
(*Function for the expansions. Note that we are throwing away second order terms*)
expan[xx_]:=Factor[Normal[Series[xx,{ϵa,0,1},{ϵp,0,1}]/.{ϵp ϵa->0}]];
(*Inverse of the metric and determinant*)
guu=expan[Inverse[gdd]];
dim=Length[gdd];
DET=expan[Det[gdd]];
(*The next are just to avoid the number of brackets in the expressions of the notebook*)
Do[g[i,j]=Simplify[expan[gdd[[i]][[j]]]];,{i,1,4},{j,1,4}];
Do[gu[i,j]=Simplify[expan[guu[[i,j]]]];,{i,1,4},{j,1,4}];

Next we jump to the computations of the scalar and vector pontentials. Remember that we have two potentials. One for the EM field and other for the complex Proca Field.

DET/.{ϵp->0 }

Complex scalar, complex proca, and the electromagnetic field expressions

Defining the scalar field. Note that because the field is complex, we are separating the perturbations into two. We also factorized the background frequency ωs for the scalar field

Φ=(ϕ0[t,r]+ϵp((δϕ1[t,r])/r) Y[χ])Exp[-I ωs t];
Φb=(ϕ0[t,r]+ϵp((δϕ2[t,r])/r) Y[χ]) Exp[I ωs t];

Defining the complex Proca expression.

(*The background Proca field and its complex conjugated*)
Cvd={f[t,r],I gf[t,r],0,0}Exp[-I ωv t];
Cvbd={f[t,r],-I gf[t,r],0,0}Exp[I ωv t];
(*The perturbations and their complex conjugated*)
δCvda=ca[t,r]{0,0,0,-(1-(χ^2))Y'[χ]}Exp[-I ωv t];
δCvdp={cfp[t,r]Y[χ],chp[t,r]Y[χ],ckp[t,r](1-(χ^2))Y[χ],0}Exp[-I ωv t];
δCvdab=cab[t,r]{0,0,0,-(1-(χ^2))Y'[χ]}Exp[I ωv t];
δCvdpb={cfb[t,r]Y[χ],chb[t,r]Y[χ],ckb[t,r](1-(χ^2))Y[χ],0}Exp[I ωv t];

Defining background plus perturbations

Cv=Cvd+ϵa δCvda+ ϵp δCvdp;
Cvb=Cvbd+ϵa δCvdab+ ϵp δCvdpb;

Defining the electromagnetic field

A0d={A0[t,r],0,0,0};
δAd=fa[t,r]{0,0,0,-(1-(χ^2))Y'[χ]};
δAp={fp[t,r]Y[χ],hp[t,r]Y[χ],kp[t,r](1-(χ^2))Y[χ],0};

Background plus perturbations

Av=A0d+ϵa δAd+ϵp δAp;

Finally, we jump into computing the field equations

Computing the field equations

Before entering into the details for the equations, let us define one rule to simplify the angular function:

ruY={Y''[χ]->(((ℓ(ℓ+1)) Y[χ]-2 χ Y'[χ])/(-1+(χ^2))),D[Y''[χ]->(((ℓ(ℓ+1)) Y[χ]-2 χ Y'[χ])/(-1+(χ^2))),χ]};

Computing the curvature terms

(*Christofel*)
Γ=ParallelTable[expan[Sum[(1/2)gu[k,l](D[g[l,i],x[[j]]]+D[g[l,j],x[[i]]]-D[g[j,i],x[[l]]]),{l,1,dim}]],{k,1,dim},{i,1,dim},{j,1,dim}];
(*Riemann tensor (Subscript[(R^a), bcd])*)
Ruddd=ParallelTable[expan[D[Γ[[b,i,k]],x[[j]]]-D[Γ[[b,i,j]],x[[k]]]+Sum[Γ[[β,i,k]]Γ[[b,β,j]]-Γ[[β,i,j]]Γ[[b,β,k]],{β,1,dim}]],{b,1,dim},{i,1,dim},{j,1,dim},{k,1,dim}];
(*Riemman Tensor (Subscript[R, abcd])*)
Rdddd=ParallelTable[expan[Sum[g[l,b]Ruddd[[l,i,j,k]],{l,1,dim}]],{b,1,dim},{i,1,dim},{j,1,dim},{k,1,dim}];
(*Ricci tensor*)
Ricci=ParallelTable[expan[Sum[Ruddd[[k,i,k,j]],{k,1,dim}]],{i,1,dim},{j,1,dim}];
(*Ricci scalar*)
scalarR=expan[Sum[Ricci[[i,j]] gu[i,j],{i,1,dim},{j,1,dim}]];
(*Einstein Tensor*)
G=ParallelTable[expan[Ricci[[i,j]]-(1/2)g[i,j]scalarR],{i,1,dim},{j,1,dim}];
Gc=ParallelTable[expan[Sum[G[[i,αp]]gu[αp,j],{αp,1,dim}]],{i,1,dim},{j,1,dim}];

Equations for the EM field

(*The electromagnetic field tensor*)
Fd=ParallelTable[expan[D[Av[[i]],x[[j]]]-D[Av[[j]],x[[i]]]],{i,1,dim},{j,1,dim}];
Fu=ParallelTable[expan[Sum[gu[αp,i]gu[βp,j]Fd[[αp,βp]],{αp,1,dim},{βp,1,dim}]],{i,1,dim},{j,1,dim}];
(*Defining covariant derivative for the scalar fields*)
DD[XX_,i_]:=(D[XX,x[[i]]]-I q Av[[i]]XX);
DDc[XX_,i_]:=(D[XX,x[[i]]]+I q Av[[i]]XX);
(*Stress energy tensor for the electromagnetic and Proca fields*)
TemF=ParallelTable[expan[
-(1/4)g[i,j]Sum[Fd[[αp,βp]]Fu[[αp,βp]],{αp,1,dim},{βp,1,dim}]
+Sum[gu[αp,βp]Fd[[i,αp]]Fd[[j,βp]],{αp,1,dim},{βp,1,dim}]],{i,1,dim},{j,1,dim}];
(*Vector equations. First, we take a look into Maxwell's equations*)
covDF=ParallelTable[expan[Sum[(1/(((-DET))^(1/2)))D[(((-DET))^(1/2))Fu[[αp,i]],x[[αp]]],{αp,1,dim}]],{i,1,dim}];
eqA=ParallelTable[expan[covDF[[i]]-I ((q )/2)Sum[gu[αp,i](DDc[Φb,αp]Φ-Φb DD[Φ,αp]),{αp,1,dim}]],{i,1,dim}];

Equations for the self-interacting complex vector field

(*The Proca field tensor*)
Bd=ParallelTable[expan[D[Cv[[i]],x[[j]]]-D[Cv[[j]],x[[i]]]],{i,1,dim},{j,1,dim}];
Bdb=ParallelTable[expan[D[Cvb[[i]],x[[j]]]-D[Cvb[[j]],x[[i]]]],{i,1,dim},{j,1,dim}];
Bu=ParallelTable[expan[Sum[gu[αp,i]gu[βp,j]Bd[[αp,βp]],{αp,1,dim},{βp,1,dim}]],{i,1,dim},{j,1,dim}];
Bub=ParallelTable[expan[Sum[gu[αp,i]gu[βp,j]Bdb[[αp,βp]],{αp,1,dim},{βp,1,dim}]],{i,1,dim},{j,1,dim}];
CC=expan[Sum[Cv[[αp]]Cvb[[βp]]gu[αp,βp],{αp,1,dim},{βp,1,dim}]];
TemB=ParallelTable[expan[
-(1/4)g[i,j]Sum[Bd[[αp,βp]]Bub[[αp,βp]],{αp,1,dim},{βp,1,dim}]
+(1/2)Sum[gu[αp,βp](Bd[[i,αp]]Bdb[[j,βp]]+Bdb[[i,αp]]Bd[[j,βp]]),{αp,1,dim},{βp,1,dim}]
+(Vv'[CC]+ϵp δVv'[CC]Y[χ])(1/2)(Cv[[i]]Cvb[[j]]+Cvb[[i]]Cv[[j]])-(1/2)g[i,j](Vv[CC]+ϵp δVv[CC]Y[χ])],{i,1,dim},{j,1,dim}];
CovDB=ParallelTable[expan[Sum[(1/(((-DET))^(1/2)))D[(((-DET))^(1/2))Bu[[αp,i]],x[[αp]]],{αp,1,dim}]],{i,1,dim}];
CovDBb=ParallelTable[expan[Sum[(1/(((-DET))^(1/2)))D[(((-DET))^(1/2))Bub[[αp,i]],x[[αp]]],{αp,1,dim}]],{i,1,dim}];
eqB=ParallelTable[expan[CovDB[[i]]+Sum[(Vv'[CC])Cv[[αp]]gu[αp,i],{αp,1,dim}]],{i,1,dim}];
eqBb=ParallelTable[expan[CovDBb[[i]]+Sum[(Vv'[CC])Cvb[[αp]]gu[αp,i],{αp,1,dim}]],{i,1,dim}];
DivC=expan[Sum[(1/(((-DET))^(1/2)))D[(((-DET))^(1/2))Sum[gu[i,j](Vv'[CC]Cv[[j]]),{j,1,dim}],x[[i]]],{i,1,dim}]];
(*Noether current*)
jbmu=(I/2)ParallelTable[expan[Sum[Bub[[i,αp]]Cv[[αp]]-Bu[[i,αp]]Cvb[[αp]],{αp,1,dim}]],{i,1,dim}];

Quantities for the scalar field

(*Stress-energy tensor*)
TemΦ=(1/2)ParallelTable[expan[DD[Φ,j]DDc[Φb,i]+DD[Φ,i]DDc[Φb,j]
-g[i,j]Sum[gu[αp,βp](DDc[Φb,αp]DD[Φ,βp]),{αp,1,dim},{βp,1,dim}]
-g[i,j](Vs[Φ Φb])],{i,1,dim},{j,1,dim}];
(*Covariant derivative
 of the scalar field*)
covD=ParallelTable[expan[
D[Av[[i]] Φ,x[[j]]]-Sum[Γ[[βp,i,j]]Av[[βp]] Φ
,{βp,1,dim}]],{i,1,dim},{j,1,dim}];
(*Charged scalar field equation*)
scaleq=expan[1/Sqrt[-DET]Sum[D[Sqrt[-DET]gu[αp,βp]D[Φ,x[[αp]]],x[[βp]]],{αp,1,dim},{βp,1,dim}]
-I q Sum[gu[αp,βp]covD[[αp,βp]],{αp,1,dim},{βp,1,dim}]
-I q Sum[gu[αp,βp] Av[[αp]] D[Φ,x[[βp]]],{αp,1,dim},{βp,1,dim}]
-(q^2)Sum[gu[αp,βp] Av[[αp]] Av[[βp]]Φ,{αp,1,dim},{βp,1,dim}]
-(Vs'[Φ Φb])Φ];
(*Noether current over scalar charge*)
Jϕmu=((-I)/2)ParallelTable[expan[Sum[gu[αp,i](Φb D[Φ,x[[αp]]]-Φ D[Φb,x[[αp]]]-2I q Av[[αp]]Φb Φ),{αp,1,dim}]],{i,1,dim}];

Finally, let us compute the Einstein’s equation

(*Einstein's equations. The ct's parameters are for control to compute particular cases and should be equal to unity to consider the general case*)
EINEQ=ParallelTable[expan[G[[i,j]]-(κ/2)(ctF TemF+ctB TemB+ctΦ TemΦ)[[i,j]]],{i,1,dim},{j,1,dim}];
Gcc=ParallelTable[expan[Sum[(EINEQ[[i,αp]])gu[αp,j],{αp,1,dim}]],{i,1,dim},{j,1,dim}];

Now with everything computed, we are able to follow and write the background equations and the equations describing the non-axial perturbations of boson stars.

Particular cases

1- Uncharged scalar boson stars

Here we display from the equations of motion previously computed the background field equations for uncharged scalar field. To do so, we turn off the contributions from the two vector fields.  We also impose that the fields are stationary.

Let us first define the rules for this particular case

rulesS={ϵp->0,ϵa->0,q->0,ctB->0,ctF->0,ctΦ->1,ϕ0->(ϕ0[#2]&),A->(A[#2]&),B->(B[#2]&)};

Equation for the metric function B

eqtt=Collect[Gcc[[1,1]]/.rulesS,{κ,ωs,(ϕ0)'[r]},FullSimplify]

Equation for the metric function A

eqrr=Collect[Gcc[[2,2]]/.rulesS,{κ,ωs,(ϕ0)'[r]},FullSimplify]

((A[r] (-1+B[r])+r B[r] A'[r])/((r^2) A[r]))//FullSimplify

Equation from the angular components

eqthth=Collect[Gcc[[3,3]]/.rulesS,{κ,ωs,(ϕ0)'[r]},FullSimplify]

This equation is already satisfied by using the two first equations and the  equation for the scalar field, that is given by

eqϕ=(scaleq/(E^(-I t ωs)))/.rulesS

showing that the angular equation is automatically satisfied. First we use the

rulesbkg=Solve[{eqtt,eqrr,eqϕ,D[eqrr,r]}==0,{B'[r],A'[r],ϕ0''[r],A''[r]}][[1]]

Simplify[eqthth/.rulesbkg]

Therefore, the equations describing spherically symmetric uncharged BSs are given by

{eqtt,eqrr,eqϕ}//Simplify

We also have the equation for the Noether current

Jϕmu/.rulesS

The above equations should be supplemented with boundary conditions in order to be solved. We shall write the integrator for this part in a different notebook file.

3- Oscillatons

Oscillatons are composed by real scalar field that oscillate together with the metric function. They are a little more difficult to deal when compared to boson stars. Here we take only the dominant term to set the differential equations.  Note, once again, that we have to turn off the contributions from the vector fields. Also, we make ωs=0 to restrict the field to real values. We focus in the mass case for the scalar potential

rulesOS={ϵp->0,ϵa->0,q->0,ctB->0,ctF->0,ctΦ->1,ωs->0,Vs->((μ^2)#&)};

Equation for the metric function B

eqtt=Simplify[(Gcc[[1,1]]/.rulesOS)]

Collect[(1/4) (-(4/(r^2))+κ (μ^2) ((ϕ0[t,r])^2)+((4 Derivative[0,1][B][t,r])/r)+B[t,r] ((4/(r^2))+κ ((Derivative[0,1][ϕ0][t,r])^2))+((κ ((Derivative[1,0][ϕ0][t,r])^2))/(A[t,r]))),{κ},Simplify]

vars=Table[Cos[i ω t],{i,1,15}];
ordermetric=Thread[Table[Cos[i ω t],{i,5,30}]->Table[0,{i,5,30}]];
orderscalar=Thread[Table[Cos[i ω t],{i,4,30}]->Table[0,{i,4,30}]];

eq01=Collect[TrigReduce[eqtt A[t,r] /.{ϕ0->(ϕ01[#2]Cos[ω #1]+ϕ02[#2]Cos[3ω #1]&),A->(A01[#2]+A02[#2]Cos[2 ω #1]+A03[#2]Cos[4 ω #1]&),B->((B01[#2]+B02[#2]Cos[2ω #1]+B03[#2]Cos[4ω #1])&)}],vars,Simplify]/.ordermetric

eqrr=Simplify[(A[t,r]Gcc[[2,2]]/.rulesOS)]

Collect[((A[t,r] (-4+(r^2) κ (μ^2) ((ϕ0[t,r])^2)+B[t,r] (4-(r^2) κ ((Derivative[0,1][ϕ0][t,r])^2)))+r (4 B[t,r] Derivative[0,1][A][t,r]-r κ ((Derivative[1,0][ϕ0][t,r])^2)))/(4 (r^2))),κ,Simplify]

eq02=Collect[TrigReduce[eqrr /.{ϕ0->(ϕ01[#2]Cos[ω #1]+ϕ02[#2]Cos[3ω #1]&),A->(A01[#2]+A02[#2]Cos[2 ω #1]+A03[#2]Cos[4 ω #1]&),B->((B01[#2]+B02[#2]Cos[2ω #1]+B03[#2]Cos[4ω #1])&)}],vars,Simplify]/.ordermetric

This equation is already satisfied by using the two first equations and the

eqf=scaleq ((A[t,r])^2)B[t,r]/.rulesOS//Simplify

eq03=Collect[TrigReduce[eqf /.{ϕ0->(ϕ01[#2]Cos[ω #1]+ϕ02[#2]Cos[3ω #1]&),A->(A01[#2]+A02[#2]Cos[2 ω #1]+A03[#2]Cos[4 ω #1]&),B->((B01[#2]+B02[#2]Cos[2ω #1]+B03[#2]Cos[4ω #1])&)}]/.orderscalar,vars,Simplify]

1

Collecting the equations

eq02

eqm01={Coefficient[eq01,Cos[2 t ω]],Coefficient[eq01,Cos[4 t ω]],eq01/.{Cos[2 t ω]->0,Cos[4 t ω]->0}};
eqm02={Coefficient[eq02,Cos[2 t ω]],Coefficient[eq02,Cos[4 t ω]],eq02/.{Cos[2 t ω]->0,Cos[4 t ω]->0}};
eqm03={Coefficient[eq03,Cos[t ω]],Coefficient[eq03,Cos[3t ω]]};

Join[eqm01,eqm02,eqm03]

eqsbkg=Simplify[Solve[Join[eqm01,eqm02,eqm03]==0,{A01'[r],A02'[r],A03'[r],B01'[r],B02'[r],B03'[r],(ϕ01)''[r],(ϕ02)''[r]}][[1]]]/.Rule->Subtract

showing that the angular equation is automatically satisfied. First we use the

rulesbkg=Solve[{eqtt,eqrr,eqϕ,D[eqrr,r]}==0,{B'[r],A'[r],ϕ0''[r],A''[r]}][[1]]

Simplify[eqthth/.rulesbkg]

Therefore, the equations describing spherically symmetric uncharged BSs are given by

{eqtt,eqrr,eqϕ}//Simplify

The above equations should be supplemented with boundary conditions in order to be solved. We shall write the integrator for this part in a different notebook file.

4- Charged scalar boson stars

Now we proceed to the case in which we have a complex scalar field interacting with a real vector field. For this case, we turn off the proca field part of the stress energy tensor. We have the following rules

rulesSC={ϵp->0,ϵa->0,ctB->0,ctF->1,ctΦ->1,ϕ0->(ϕ0[#2]&),A->(A[#2]&),B->(B[#2]&),A0->(A0[#2]&)};

Equation for the metric function B

eqtt=Collect[Gcc[[1,1]]/.rulesSC,{κ,ωs,(ϕ0)'[r],q,A0'[r]},FullSimplify]

Equation for the metric function A

eqrr=Collect[Gcc[[2,2]]/.rulesSC,{κ,ωs,(ϕ0)'[r],q,A0'[r]},FullSimplify]

Equation from the angular components

eqthth=Collect[Gcc[[3,3]]/.rulesSC,{κ,ωs,(ϕ0)'[r],q,A0'[r]},FullSimplify]

This equation is already satisfied by using the two first equations and the equation for the scalar field, that is given by

eqϕ=(scaleq/(E^(-I t ωs)))/.rulesSC

And the equation for the vector field

eqvec=Collect[(((eqA[[1]]/.rulesSC))/((B[r])/(A[r]))),{A0''[r],A0'[r]},Simplify]

showing that the angular equation is automatically satisfied. First we use the

rulesbkg=Solve[{eqtt,eqrr,eqϕ,D[eqrr,r],eqvec}==0,{B'[r],A'[r],ϕ0''[r],A''[r],A0''[r]}][[1]]

Simplify[eqthth/.rulesbkg]

Therefore, the equations describing spherically symmetric uncharged BSs are given by

{eqtt,eqrr,eqϕ,eqvec}//Simplify

Noether current

Jϕmu/.rulesSC//Simplify

The above equations should be supplemented with boundary conditions in order to be solved. We shall write the integrator for this part in a different notebook file.

5- Vector boson stars

Finally let go through the vector boson star case.  To do so, we turn off the contributions from the scalar and electromagnetic fields.

Clear[Vv]

rulesSB={ϵp->0,ϵa->0,ctB->1,ctF->0,ctΦ->0,ϕ0->(0&),A->(A[#2]&),B->(B[#2]&),A0->(0&),f->(f[#2]&),gf->(gf[#2]&)};

Equation for the metric function B

eqtt=Collect[Gcc[[1,1]]/.rulesSB,{κ,ωs,(ϕ0)'[r],q,A0'[r]},FullSimplify]

Equation for the metric function A

eqrr=Collect[Gcc[[2,2]]/.rulesSB,{κ,ωs,(ϕ0)'[r],q,A0'[r]},FullSimplify]

Equation from the angular components

eqthth=Collect[Gcc[[3,3]]/.rulesSB,{κ,ωs,(ϕ0)'[r],q,A0'[r]},FullSimplify]

And the equation for the vector field

eqvec=(E^(I t ωv))eqB/.rulesSB//Simplify

showing that the angular equation is automatically satisfied. First we use the

DivC/.rulesSB//Simplify

Due to the complexity, let us focus on just the mass term

rulesbkg=Solve[{eqtt,eqrr,D[eqrr,r],eqvec[[1]],eqvec[[2]],DivC/.rulesSB}==0/.{Vv->(((μf)^2)#&)},{B'[r],A'[r],f''[r],f'[r],A''[r],gf'[r]}][[1]]

Simplify[eqthth/.rulesbkg/.{Vv->(((μf)^2)#&)}]

Therefore, the equations describing spherically symmetric uncharged BSs are given by

{eqtt,eqrr,eqvec[[1]],eqvec[[2]],DivC/.rulesSB}//Simplify

The above equations should be supplemented with boundary conditions in order to be solved. We shall write the integrator for this part in a different notebook file.

jbmu/.rulesSB/.i->I//Simplify

((B[r] gf[r] (-ωv gf[r]+f'[r]))/(A[r]))

6- General cases

eqtt=Collect[Gcc[[1,1]]/.{ϵa->0,ϵp->0,ctΦ->1,ctB->1,ctF->1,ωs->0,ωv->0},{κ},FullSimplify]

eqrr=Collect[Gcc[[2,2]]/.{ϵa->0,ϵp->0,ctΦ->1,ctB->1,ctF->1,ωs->0,ωv->0},{κ},FullSimplify]

eqthth=Collect[Gcc[[3,3]]/.{ϵa->0,ϵp->0,ctΦ->1,ctB->1,ctF->1,ωs->0,ωv->0},{κ},FullSimplify]

DivC/.{ϵa->0,ϵp->0,ctΦ->1,ctB->1,ctF->1,ωs->0,ωv->0}//Simplify

Trying perturbations

5- Vector boson stars

Finally let go through the vector boson star case.  To do so, we turn off the contributions from the scalar and electromagnetic fields.

rulesSB={ϵp->0,ctB->1,ctF->0,ctΦ->0,ϕ0->(0&),A->(A[#2]&),B->(B[#2]&),A0->(0&),f->(f[#2]&),gf->(gf[#2]&)};

Einstein’s equations

Gcc[[3,4]]/.rulesSB//.ruY//Simplify

Gcc[[1,4]]/.rulesSB//.ruY//Simplify

Gcc[[2,4]]/.rulesSB//.ruY//Simplify

Proca equations

Coefficient[eqB/.rulesSB//.ruY,ϵa]//Simplify

Coefficient[eqBb/.rulesSB//.ruY,ϵa]//Simplify