Add Brans-Dicke in Jordan frame

1 files changed, 267 insertions, 0 deletions

diff --git a/03-brans-dicke-jor.mac b/03-brans-dicke-jor.mac
new file mode 100644
index 0000000..4e3599b
--- /dev/null
+++ b/03-brans-dicke-jor.mac
@@ -0,0 +1,267 @@
+load("setup.mac")$
+load("tensors.mac")$
+load("series.mac")$
+
+
+/* Equations */
+
+depends(λ, [u,r,θ,ϕ]);
+
+/* Einstein tensor */
+for i : 1 thru 4 do for j : 1 thru 4 do (
+  G_ab[i,j] : expand(ric[i,j] - ric_s*lg[i,j]/2));
+
+/* Stress-energy tensor of the theory */
+for i : 1 thru 4 do for j : 1 thru 4 do (
+  T_ab[i,j] : (ω/8/%pi/λ)*(
+    cdiff(λ,i)*cdiff(λ,j)
+    - lg[i,j]*(contr(ug[ρ,σ]*cdiff(λ,ρ)*cdiff(λ,σ),[ρ,σ])/2)));
+
+/* Einstein equations tensor */
+for i : 1 thru 4 do for j : 1 thru 4 do (
+  E_ab[i,j] : expand(G_ab[i,j] - (1/λ)*(
+      8*%pi*T_ab[i,j]
+      + cdiff(cdiff(λ, j), i)
+      - contr(mcs[i,j,k]*cdiff(λ,k),[k])
+      - lg[i,j]*box_op(λ))))$
+
+/* Einstein equations tensor with mixed indices */
+for i : 1 thru 4 do for j : 1 thru 4 do (
+  uE_ab[i,j] : expand(contr(ug[i,k]*E_ab[k,j], [k])))$
+
+/* solve(uE_ab[1,2]=0, diff(β,r))[1], */
+
+/* Hypersurface equations */
+hs_eqs : [
+/* Determine β */
+uE_ab[1,2]=0,
+/* Determine A */
+uE_ab[1,3]=0,
+/* uE_ab[1,4]=0 is trivial, would determine A' otherwise */
+/* Determine B */
+uE_ab[1,1]=0
+]$
+
+/* Evolution equations */
+ev_eqs : [
+/* Determine γ */
+E_ab[3,3]=0,
+/* Should be fulfilled automatically (?) */
+E_ab[4,4]=0
+/* Fulfilled trivially */
+/* E_ab[3,4]=0 */
+/* E_ab[4,3]=0 */
+]$
+
+/* Suplementary conditions */
+sc_eqs : [
+uE_ab[2,1]=0,
+uE_ab[2,3]=0
+/* Fulfilled trivially */
+/* uE_ab[2,4]=0 */
+]$
+
+hs_eqs : trigreduce(hs_eqs)$
+sc_eqs: trigreduce(sc_eqs)$
+
+hs_eqs : subst([diff(λ,ϕ)=0], hs_eqs)$
+ev_eqs : subst([diff(λ,ϕ)=0], ev_eqs)$
+sc_eqs : subst([diff(λ,ϕ)=0], sc_eqs)$
+
+/* Scalar wave equation */
+sw_eq : trigreduce(box_op(λ)=0)$
+
+
+/* Series ansatz */
+
+expansion_order : 6;
+
+/* Sequence elements inherit dependance on u and θ */
+depends([c, H, M, N, γc, βc, Ac, Bc, λc], [u, θ])$
+
+/* Main series ansatz */
+build_main_ansatz(n) := block([], [
+  γ = c/r + neg_pow_series(γc, r, 2, n),
+  β = H + neg_pow_series(βc, r, 1, n),
+  A = neg_pow_series(Ac, r, 1, n),
+  B = neg_pow_series(Bc, r, -1, n),
+  λ = neg_pow_series(λc, r, 0, n)])$
+
+main_ansatz : build_main_ansatz(expansion_order)$
+
+funcs : [A, B, β, γ, β_expr, γ_expr, λ];
+
+
+/**/
+
+/* Simplify ahead to speed up calculation */
+γβ_facts : [H=0, γc[2]=0];
+sim_facts : ev([
+  /* From Equation 4 (1/r term) */
+  diff(Φc[0],u)=0,
+  diff(Φc[0],u,2)=0,
+  /* From Equation 4 (1/r^2 term) */
+  diff(Φc[0],θ)=0,
+  diff(Φc[0],θ, 2)=0,
+  diff(diff(Φc[0],θ),u)=0,
+  diff(diff(Φc[0],u),θ)=0],
+  eval,diff);
+
+coeffs : append(γβ_facts, sim_facts);
+
+/* Equation 1 : determine β */
+e_hs_eqs[1] : expand_equation(expand(exp(2*β)*λ^2*hs_eqs[1]), coeffs)$
+ce_hs_eqs[1] : coeff_equations(e_hs_eqs[1])$
+β_coeffs : block(
+  [vals : []],
+  for i : 1 thru 5 do (
+    vals : append(vals, solve(subst(vals, ce_hs_eqs[1][i][2]), βc[i]))),
+  vals)$
+
+coeffs : append(coeffs, β_coeffs)$
+
+/* Equation 2 : determine A */
+e_hs_eqs[2] : expand_equation(expand(exp(2*β)*λ^2*hs_eqs[2]), coeffs)$
+/* BUG: diff(λc[0],θ)=0 does not get substituted */
+e_hs_eqs[2] : subst(coeffs, e_hs_eqs[2])$
+ce_hs_eqs[2] : coeff_equations(e_hs_eqs[2])$
+
+A_coeffs : block(
+  [vals : []],
+  for i : 1 thru expansion_order do (
+    if not(i=3) then (
+      vals : append(vals, solve(subst(vals, ce_hs_eqs[2][i][2]), Ac[i])))),
+  ev(vals, eval, diff, expand))$
+
+coeffs : append(coeffs, A_coeffs)$
+
+/* Equation 3 : determine B */
+print_progress : true;
+e_hs_eqs[3] : expand_equation(expand(exp(2*β)*λ^2*hs_eqs[3]), coeffs)$
+/* TODO: 17 31 34 slow due to squares? */
+/* BUG: diff(λc[0],θ)=0 does not get substituted; due to missing ev? */
+e_hs_eqs[3] : subst(coeffs, e_hs_eqs[3])$
+ce_hs_eqs[3] : coeff_equations(e_hs_eqs[3])$
+B_coeffs : block(
+  [vals : []],
+  for i : 1 thru expansion_order do (
+    if i=2 then (
+      vals : append(vals, [Bc[0]=-2*M]))
+    else (
+      vals : append(vals, solve(subst(vals, ce_hs_eqs[3][i][2]), Bc[i-2])))),
+  ev(vals, eval, diff, expand))$
+
+coeffs : append(coeffs, B_coeffs)$
+
+
+/* Evolution equations */
+
+depends(D, [u, θ]);
+
+/* Equation 4 : determine λ */
+e_sw_eq : expand_equation(sw_eq, coeffs)$
+ce_sw_eq : coeff_equations(e_sw_eq)$
+λu_coeffs : block(
+  [vals : append(sim_facts, [λc[1]=D])],
+  for i : 3 thru 6 do (
+    vals : append(vals, solve(subst(vals, ce_sw_eq[i][2]), diff(λc[i-1],u)))),
+  vals)$
+
+coeffs : append(coeffs, λu_coeffs)$
+
+depends(C, [u, θ]);
+extra_coeffs : [γc[3]=C-c^3/6, Ac[3]=2*N + 3*c*diff(c,θ) + 4*c^2*cot(θ)]$
+coeffs : append(coeffs, extra_coeffs)$
+
+/* Equation 5 : determine γ */
+e_ev_eqs[2] : expand_equation(expand(exp(2*β)*exp(2*γ)*λ^2*ev_eqs[2]/sin(θ)^2), coeffs)$
+/* BUG: diff(λc[0],θ)=0 does not get substituted; due to missing ev? */
+e_ev_eqs[2] : subst(coeffs, e_ev_eqs[2])$
+ce_ev_eqs[2] : coeff_equations(e_ev_eqs[2])$
+
+eq_quad : trigreduce(expand(trigsimp(subst(coeffs,
+        4*solve(ce_ev_eqs[2][2][2], diff(C, u))[1]))));
+
+/* Relation of c to c_AB */
+ce_ev_eqs[2][1]; /* TODO: D_u * c = c_u * D ??? */
+
+
+/* Supplementary conditions */
+
+/* Supplementary Condition 1 */
+e_sc_eqs[1] : expand_equation(expand(λ^2*sc_eqs[1]), coeffs)$
+/* BUG */
+e_sc_eqs[1] : subst(coeffs, e_sc_eqs[1])$
+eq_mono : block(
+  [eq : coeff(e_sc_eqs[1], r, -2)=0],
+  trigreduce(expand(solve(eq, diff(M,u)))));
+
+/* Supplementary Condition 2 */
+e_sc_eqs[2] : expand_equation(sc_eqs[2], coeffs);
+eq_dipo : block(
+  [eq : coeff(e_sc_eqs[2], r, -2)=0],
+  trigreduce(expand(-3*solve(eq, diff(N,u)))));
+
+
+/* Cache */
+
+path(name) := printf(false, "cache/03-brans-dicke-jor/~a.lisp", name);
+
+save(path("e_hs_eqs"), e_hs_eqs);
+save(path("ce_hs_eqs"), ce_hs_eqs);
+save(path("e_sw_eq"), e_sw_eq);
+save(path("ce_sw_eq"), ce_sw_eq);
+save(path("e_ev_eqs"), e_ev_eqs);
+save(path("ce_ev_eqs"), ce_ev_eqs);
+
+save(path("beta_coeffs"), β_coeffs);
+save(path("a_coeffs"), A_coeffs);
+save(path("b_coeffs"), B_coeffs);
+save(path("lambdau_coeffs"), λu_coeffs);
+
+save(path("e_sc_eqs"), e_sc_eqs);
+save(path("coeffs"), coeffs)
+
+save(path("eq_quad"), eq_quad);
+save(path("eq_mono"), eq_mono);
+save(path("eq_dipo"), eq_dipo);
+
+load_cache() := block(
+  [cache : [
+  "e_hs_eqs", "ce_hs_eqs",
+  "e_sw_eq", "ce_sw_eq",
+  "e_ev_eqs", "ce_ev_eqs",
+  "beta_coeffs", "a_coeffs", "b_coeffs", "lambdau_coeffs",
+  "e_sc_eqs",
+  "coeffs",
+  "eq_quad", "eq_mono", "eq_dipo"]],
+  for file in cache do load(path(file)));
+
+
+/* Misc. checks */
+
+p1 : (4/r + 2*diff(λ,r)/λ)*diff(β,r);
+p2 : - diff(λ,r,2)/λ;
+p3 : - ω*diff(λ,r)^2/λ^2;
+p4 : - sum(sum(sum(sum(uh_AB[A,B]*uh_AB[C,D]*diff(h_AB[A,C],r)*diff(h_AB[B,D],r), A, 3, 4), B, 3, 4), C, 3, 4), D, 3, 4)/4;
+eq1 : p1 + p2 + p3 + p4;
+
+expand(λ^2 * (p1 + p2 + p3 + p4))
+
+expand(λ^2 * (p1))
+
+expand(exp(2*β)*λ*hs_eqs[1])
+
+expand(λ^2 * (p4))
+
+uE_ab[1,2]=0;
+
+tmp : subst([diff(diff(λc[0], u), θ)=0, diff(λc[0], u, 2)=0],eq_mono)
+
+tmp : trigreduce(expand(trigsimp(tmp)))
+
+λu_coeffs[5]
+
+subst(diff(λu_coeffs[5], u), ev(tmp, eval))
+
+/* define eth and compare with paper */