# thenumbernine/Relativity

general relativity 3+1 formalism simulation
Fetching latest commit…
Cannot retrieve the latest commit at this time.
Type Name Latest commit message Commit time
Failed to load latest commit information.
include
src
test
Makefile
buildinfo

```GR simulation

Sources:
- "Gravitation", Misner, Thorne, Wheeler. (c) 1973
- "Numerical Relativity: Solving Einstein's Equations on the Computer", Baumgarte, Shapiro. (c) 2010
- "Introduction to 3+1 Numerical Relativity", Alcubierre. (c) 2008

Starting again from scratch now that I've gone through the ADM chapter of Misner, Thorne, & Wheeler and the 3+1 decomposition chapter of Baumgarte & Shapiro
I read somewhere that the other two books are the first two books published on the topic of numerical relativity ( http://adsabs.harvard.edu/abs/2011CQGra..28v9002G ).
I've observed this as I am progressing through the two simultaneously.  It seems Baumgarte and Shapiro are better at explaining things, taking smaller steps and providing problems for the reader.
I'm a big fan of this, seeing as I need all the help I can get.  Alcubierre, on the other hand, does a better job with the finishing details: offering in certain places more complete descriptions of solutions
(or at least more information up front) and occasionally alternative means of getting to them.  I'm only on chapter 3 of both, so maybe this will change later.

TODO

algorithm:
- lapse and shift functions
- Riemann solver implementation.  Baumgarte & Shapiro p.384 and Alcubierre p.5.6 discuss this.
- currently I'm running on an explicit integrator.  need to use implicit ones.

implementation:
- general-case covariant derivative would be nice.
this depends on the ability to tell whether a certain index is upper or lower,
which is difficult now that I made templates successively nested structures (so that symmetric tensors could take advantage of size optimizations).
- general-case second partial derivative? to construct the result as a symmetric-lower.
- symmetric & antisymmetric indexes must be next to one another.  this is courtesy of the nested implementation of the templates.
tensor<lower,lower> is - from the view of the data structure - a one-form of one-forms, and tensor<upper,upper> is a vector of vectors.
I'm not even sure how the templates could be changed to accomodate anything more flexible: maybe put the information last in the template arguments
like so tensor<L,U,U,L,U, symmetric<0,1>, antisymmetric<2,3,4>>.  that might be a nightmare to implement,
and i haven't come across any situations where the current system hasn't failed me.
- antisymmetric tensors have to be the outer-most indexes.  this is because only one half of an antisymmetric pair of elements is stored.
this is made compatible with read and write operations by returning a generated accessor object which does the sign-flipping when appropriate.
i think all that may need to be done here is to give the accessor class itself a dereference operator to forward on through templates
- which means specializations of itself when it returns reals and when it doesn't, but there may be more.
- separate read- and write-iterators.  read-iterators iterating over all read indexes. write iterators iterating over all written variables.
the two are not one-to-one.  symmetric matrixes have n*n read indexes but only n*(n+1)/2 variables capable of writing. antisymmetric have n*(n-1)/2.
doing this involves using the current dimension-per-rank-based method for read iterator bounds, but write iterators would have delegate to each index
the responsibility for incrementing their own offset, size, and amount (rank) of indexes in the iteration.
- index notation expressions via expression templates.  no more for-loops.  simply "beta_l(i) = gamma_ll(i,j) * beta_u(j);"
this depends on the separate write iterators.  I've already got this in my Lua tensor and symbolic math libraries.  In Lua.
- permute and constant indexes.  like boost::bind.  too bad I'm not using boost.

Conventions & Notations:

indexes: (taken from Baumgarte & Shapiro p.26)

a-h, o-z:	4D -- most often txyz. not sure if t should be stored first (for math index notation consistency) or last (ugly casting between 3D and 4D rank-1 tensors) just yet.
i-n:		3D -- most often xyz.

operators:

partial_u 	= partial derivative wrt coordinate u
diff_u 		= covariant derivative wrt coordinate u
D_u 		= projection derivative / covariant of spatial metric wrt coordinate u
= proj diff_v
proj 		= projection operator (Baumgarte & Shapiro eqn 2.31)
such that proj T^a_b = gamma^a_c gamma_b^d T^c_d

variables:

alpha		= lapse (distance between timeslices a coordinate is transported over dt)

beta^a		= shift (distance within the 3D hypersurface a coordiate is transported over dt)
= (0, beta^i) 	(i.e. beta^t = 0), therefore beta^a beta_a = beta^k beta_k

g_ab 		= ADM representation of 4D metric tensor
= gamma_ab - n_a n_b
= ( -alpha^2 + beta^k beta_k	,	beta_j)		<- where beta_j = gamma_ji beta^i
(			beta_i				,	gamma_ij)		\ which is equivalent to beta_j = gamma_ja beta^a since beta^t = 0 after all so we just neglect the timelike component of beta and gamma

g^ab		= contravariant/inverse ADM representation of 4D metric tensor
= gamma^ab - n^a n^b
= ( -1/alpha^2		,			beta^j/alpha^2			   )
(	beta^i/alpha^2	,	gamma^ij - beta^i beta^j / alpha^2 )

n_a			= covariant form of normal vector (i.e. normal one-form) to 3D hypersurface in 4D space
= (-alpha, 0)

n^a			= contravariant form of normal vector to 3D hypersurface in 4D space
= (1/alpha, -beta^i/alpha)
Note that <n,n> = n_a n^a = -alpha/alpha - 0*beta^i/alpha = -1

a^a			= acceleration vector of 3D hypersurface normal
= diff_n n = n^b diff_b n^a
A related, useful identity is 0 = diff_a -1 = diff_a (n^b n_b) = n^b diff_a n_b + n_b diff_a n^b <=> n^b diff_a n_b = 0

gamma^ab	= 4D projection operator / contravariant/inverse metric tensor of 3D spatial hypersurface
= g^ab + n^a n^b
= (	0	,		0 	 )
( 0	,	gamma_ij )

gamma_ab 	= 4D projection operator / metric tensor of 3D spatial hypersurface
= g_ab + n_a n_b
= ( beta^k beta_k	, 	beta_j   )
( 	beta_i		,	gamma_ij )

K_ab		= 3D->4D extrinsic curvature
= -gamma_a^c gamma_b^d diff_c n_d

conn^ijk	= 3D hypersurface Christoffel symbols of the first kind
= 1/2 (partial_k gamma_ij + partial_j gamma_ik - partial_i gamma_jk)

conn^i_jk	= 3D hypersurface Christoffel symbols of the second kind (connection coefficients)
= gamma^il conn_ljk
= 1/2 gamma^il (partial_k gamma_lj + partial_j gamma_lk - partial_l gamma_jk)

R4^a_bcd	= 4D Riemann curvature tensor
= partial_c conn^a_bd - partial_d conn^a_bc + conn^a_uc conn^u_bd - conn^a_ud conn^u_bc

R^i_jkl		= 3D hypersurface Riemann curvature tensor
= partial_k conn^i_jl - partial_l conn^i_jk + conn^i_uk conn^u_jl - conn^i_ul conn^u_jk

R4_ab		= 4D Ricci curvature tensor
= R^u_aub

R_ij		= 3D hypersurface Ricci curvature tensor
= R^k_ikj
= conn^k_ij,k - conn^k_ik,j + conn^k_lk conn^l_ij - conn^k_lj conn^l_ik
= 1/2 gamma^kl (partial_i partial_l gamma_kj + partial_k partial_j gamma_il - partial_i partial_j gamma_kl - partial_k partial_l gamma_ij) + gamma^kl (conn^m_il conn_mkj - conn^m_ij conn_mkl)

R4			= 4D Gaussian curvature
= g^ab R_ab

R			= 3D hypersurface Gaussian curvature
= gamma^ij R_ij

T_ab		= stress-energy tensor

S_ab		= spatial stress energy
= proj T_ab = gamma_a^c gamma_b^d T_cd

S_a			= momentum density
= -gamma_a^b n^c T_bc

rho			= energy density
= n^a n^b T_ab = n_a n_b T^ab

scratch paper:

is the covariant derivative of a symmetric tensor also symmetric?

sum upper-ranks as follows:
diff_k x^i = partial_k x^i + conn^i_jk x^j

sum lower-ranks as follows:
diff_k x_i = partial_k x_i - conn^j_ik x_j

is symmetry maintained?:
t^ij == t^ji <=> diff_k t^ij == diff_k t^ji ?
diff_k t^ij = partial_k t^ij + t^lj conn^i_lk + t^il conn^j_lk
diff_k t^ji = partial_k t^ji + t^li conn^j_lk + t^jl conn^i_lk
= partial_k t^ij + t^lj conn^i_lk + t^il conn^j_lk
yes if both are upper
t_ij == t_ji <=> diff_k t_ij == diff_k t_ji ?
diff_k t_ij = partial_k t_ij - t_lj conn^l_ik - t_il conn^l_jk
diff_k t_ji = partial_k t_ji - t_li conn^l_jk - t_jl conn^l_ik
= partial_k t_ij - t_lj conn^l_ik - t_il conn^l_jk
yes if both are lower
t^i_j == t^j_i <=> diff_k t^i_j == diff_k t^j_i ?
diff_k t^i_j = partial_k t^i_j + t^l_j conn^i_lk - t^i_l conn^l_jk
diff_k t^j_i = partial_k t^j_i + t^l_i conn^j_lk - t^j_l conn^l_ik
= partial_k t^i_j - t^l_j conn^l_ik + t^i_l conn^j_lk
here we get a sign issue.   if conn^l_ik = -conn^i_lk then we're good.

g^a_b = delta^a_b = constant, so (g^a_b),c = 0
g^ad_,c g_db = -g^ad g_db,c		<= for g dg, switch the lowers and uppers to switch the signs

conn^l_ik = 1/2 g^lm (g_mi,k + g_mk,i - g_ik,m)
conn^i_lk = 1/2 g^im (g_ml,k + g_mk,l - g_lk,i)
= 1/2 (g^im g_ml,k + g^im g_mk,l - g^im g_lk,i)
= -1/2 (g_lm g^mi,k + g^im g_lk,i - g^im g_mk,l)

... I'm not buying it, nor do I see it recorded anywhere

when raising the extrinsic curvature tensor, do the time components of the metric matter?
(or can it be calculated solely by the spatial metric / projection operator?)

both Baumgarte & Shapiro and Alcubierre make mention that the extrinsic curvature only matters for the spatial components
Misner, Thorne, & Wheeler and Alcubierre define K_ab = -proj diff_a n_b, and Baumgarte & Shapiro as K_ab = -proj diff_(a n_b) (where T_(a,b) = 1/2 (T_ab + T_ba) ).
Alcubierre mentions how only the spatial components of diff_a n_b are symmetric, but that due to projection we don't have to worry about the others. (p.69)
Baumgarte & Shapiro mention how the tensor itself is only defined wrt spatial components (p.34)
They go on to explain that the the trace is given by K = g^ab K_ab = gamma^ab K_ab, which would hold true if the time components of K_a0 were zero (as those of gamma^a0 are)

ultimately ...
K_ab = -proj diff_a n_b
= -gamma_a^c gamma_b^d diff_c n_d
= -(g_a^c + n_a n^c) (g_b^d + n_b n^d) (diff_c n_d)
= -diff_a n_b - n_a n^c diff_c n_b - n_b n^d diff_a n_d - n_a n_b n^c n^d diff_c n_d
= -diff_a n_b - n_a a_b - n_b * 0 - n_a n_b n^c * 0
= -diff_a n_b - n_a a_b - n_b * 0 - n_a n_b n^d a_d (alternatively)
= -diff_a n_b - n_a a_b - n_b * 0 - n_a n_b * 0
= -diff_a n_b - n_a a_b

so K_ab = (diff_0 n_0, diff_0 n_j)   (n_0 a_0, n_0 a_j)
(diff_i n_0, diff_i n_j) - (n_i a_0, n_i a_j)

...but n_a = (-alpha, 0), so n_i = 0

K_ab = (-diff_0 alpha, 0)   (-alpha a_0, -alpha a_j)
(-diff_i alpha, 0) - (	0	   , 	0	   )

K_ab = (-alpha a_0 - diff_0 alpha, 	-alpha a_j)
(		-diff_i alpha	 ,		0	  )

...and from there it looks to me that the spatial components of K_ij are always zero.
What am I doing wrong?

Another note: if you use the metric to raise the spatial tensor, i.e. gamma_a^c = gamma_ab g^bc then you get a non-symmetric tensor
that may or may not have a beta in the time/space components.
Despite this the book denotes the gammas with no particular order to a & b, as if gamma_a^b == gamma^a_b.
This would imply (as I read on p.54) that the components of gamma_at are zero as well as gamma^at
And this would mean that the identity g^ab = gamma^ab + n^a n^b is only true for the contravariant case
(leaving the covariant case to not be subject to the same identity...)

momentum density term:

S_ab = gamma_a^c gamma_b^d T_cd
S^ab = gamma^ac gamma^bd T_cd
= gamma^ai gamma^bj T_ij (since gamma^at = 0)
S^ij = gamma^ik gamma^jl T_kl

Baumgarte & Shapiro p.47
S^i = -gamma^ij n^a T_aj

S = S^a_a
= g^ab S_ab
= g^ab gamma_a^c gamma_b^d T_cd
= gamma^cb gamma_b^d T_cd
= gamma^cd T_cd
= gamma^ij T_ij (since the other gamma^ab components are zero)
(Baumgarte & Shapiro p.47) gives the following:
S_ij = gamma_ia gamma_jb T^ab	<- which may or may not include the betas in gamma_ab = g_ab - n_a n_b
S = gamma^ij S_ij
= gamma_ai gamma^ij gamma_jb T^ab
= (gamma_ac gamma^cj - gamma_at gamma^tj) gamma_jb T^ab
= (gamma_a^j - gamma_at gamma^tj) gamma_jb T^ab
= (gamma_a^j gamma_jb - gamma_at gamma^tj gamma_jb) T^ab
= (gamma_a^c gamma_cb - gamma_a^t gamma_tb - gamma_at gamma^tj gamma_jb) T^ab
= (gamma_ab - gamma_a^t gamma_tb - gamma_at (gamma^tj gamma_jb)) T^ab
= (gamma_ab - gamma_a^t gamma_tb - gamma_at (gamma^td gamma_db - gamma^tt gamma_tb)) T^ab
= (gamma_ab - gamma_a^t gamma_tb - gamma_at gamma^t_b + gamma_at gamma^tt gamma_tb) T^ab
There are our extra spatial terms.

Kerr-Schild coordinates in terms of lapse and shift
(I'm sure someone else has done this already -- and done a better job)
I'm getting this from Alcubierre p.56 and Baumgarte & Shapiro p.52

Boyer-Lindquist coordinates:
ds^2 = -(1 - 2 M r / rho^2) - 4 M a r sin(theta)^2 / rho^2 dt dphi + rho^2 / delta dr^2 + rho^2 dtheta^2 + sin(theta)^2 / rho^2 ((r^2 + a^2)^2 - a^2 delta sin(theta)^2) dphi^2
delta = r^2 - 2 M r + a^2
rho^2 = r^2 + a^2 cos(theta)^2
a = J / M
J = angular momentum
M = mass

Kerr-Schild coordinates:
ds^2 = (eta_uv + 2 H l_u l_v) dx^u dx^v
H = (r M - Q^2 / 2) / (r^2 + a^2 z^2/r^2)
l_u = (1, (r x + a y) / (r^2 + a^2), (r y - a x) / (r^2 + a^2), z/r)
eta_uv l^u l^v = g_ab l^a l^b = 0	<- null for both metrics
a = J / M = angular momentum density
M = mass
J = angular momentum
Q = total electric charge

separating out into like components matching our ADM metric:

g_tt = -1 + 2 H
g_ti = g_it = H l_i
g_ij = eta_ij + (2 - eta_ij) H l_i l_j

g_tt = -alpha^2 + beta^k beta_k
g_ti = g_it = beta_i = gamma_ij beta^j
g_ij = gamma_ij

...so...

gamma_ij = eta_ij + (2 - eta_ij) H l_i l_j
beta_i = H l_i
<=> beta^i gamma_ij = H l_j
<=> beta^i = H l_j gamma^ji
-alpha^2 + beta^k beta_k = g_tt = -1 + 2 H
<=> alpha^2 = 1 - 2 H - beta^k beta_k
<=> alpha = sqrt(1 - 2 H - beta^k beta_k)

and to boot (Baumgarte & Shapiro p.51):
K_ij = 2 H alpha / r (eta_ij - (2 + H) l_i l_j)

Constraints:

(Most this is taken from Ch.2 of Baumgarte & Shapiro.
I've got all the intermediate steps in a notebook, maybe someday I'll copy it all into here.  It's pretty long.)

Einstein's equation:
G_ab = R4_ab - 1/2 R4 g_ab = 8 pi T_ab

Gauss' equation:
R_abcd + K_ac _Kbd - K_ad K_cb = proj R4_abcd
R_abcd + K_ac _Kbd - K_ad K_cb = gamma^p_a gamma^q_b gamma^r_c gamma^s_d R4_pqrs

Codazzi-Maniardi:
D_b K_ac - D_a K_bc = gamma^p_a gamma^q_b gamma^r_c n^s R4_pqrs

Einstein's equation contracted twice with n:
2 n^a n^b G_ab = gamma^ac gamma^bd R4_abcd

Gauss' equation contracted twice with n:
R4 + 2 n^a n^b R4_ab = R + K^2 - K_ab K^ab

Stress-Energy contracted twice with n:
rho = T_ab n^a n^b

... substitute ...

R + K^2 - K_ab K^ab = 2 n^a n^b G_ab
R + K^2 - K_ab K^ab = 16 pi rho

Hamiltonian:
H = n^a n^b (G_ab - 8 pi T_ab)
H = n^a n^b G_ab - 8 pi rho
H = 1/2 (R + K^2 - K_ba K^ab) - 8 pi rho
H = 1/2 (R + K^2 - K_ji K^ij) - 8 pi rho	<- if K is purely spatial
H = 1/2 (R + K^2 - K^j_i K^i_j) - 8 pi rho

Momentum:
S_a = -gamma_a^b n^c T_bc
D_b K^b_a - D_a K = 8 pi S_a

D_j (K^ij - gamma^ij K) = 8 pi S^i
M^i = D_j (K^ij - gamma^ij K) - 8 pi S^i

Metric Conformal Factor:

psi = conformal factor
gammaBar_ij = background metric
gamma_ij = psi^4 gammaBar_ij

gamma = det(gamma_ij)	<- by definition
det(gamma_ij * s) = s^3 gamma	<- 3D determinant of scaled matrix
1 = gamma (gamma^(-1/3))^3 = det(gamma_ij * gamma^(-1/3))
gammaBar_ij = gamma^(-1/3) gamma_ij	<- such that det(gammaBar_ij) = 1

psi^4 = gamma^(1/3)
psi = gamma^(1/12)

Baumgarte & Shapiro p.56 provides a convenient formula for deriving the metric (and properties) from the conformal metric:

connection coefficients:
conn^i_jk = connBar^i_jk + 2 (delta^i_j DBar_k ln(psi) + delta^i_k DBar_j ln(psi) - gammaBar_jk gammaBar^il DBar_l ln(psi))

Ricci curvature:
R_ij = RBar_ij - 2 (DBar_i DBar_j ln(psi) + gammaBar_ij gammaBar^lm DBar_l DBar_m ln(psi)) + 4 ((DBar_i ln(psi)) (DBar_j ln(psi)) - gammaBar_ij gammaBar^lm (DBar_l ln(psi)) (DBar_m ln(psi)))

Gaussian curvature:
R = psi^-4 RBar - 8 psi^-5 gammaBar^ij DBar_i DBar_j psi

Hamiltonian constraint:
8 gammaBar^ij DBar_i DBar_j psi - psi RBar - psi^5 K^2 + psi^5 K_ij K^ij = -16 pi psi^5 rho

From our evolved ln(sqrt(gamma)) we get:
gamma^(1/2) = psi^6
ln(sqrt(gamma)) = 6 ln(psi)

all in all...

gammaBar_ij = psi^-4 gamma_ij
gamma_ij = psi^4 gammaBar_ij
gammaBar^ij = psi^4 gamma^ij
gamma^ij = psi^-4 gammaBar^ij

Conformal Traceless Extrinsic Curvature Tensor

A_ij = K_ij - 1/3 gamma_ij K

The 1/3 is only for the 3rd dimension -- just as all those ^3's are in the conformal factor.
But I'm not going to go changing things to be dimension-invariant just yet. For now our 1D case will be just line-simulating a 3D case.

ABar^ij = psi^10 A^ij	<- this is where Baumgarte & Shapiro start.  Thay say the value is arbitrary then show how 10 is a useful number.  From there ...
A^ij = psi^-10 ABar^ij
ABar^ij = psi^10 A^ij

A_ij = gamma_ik A^kl gamma_lj
= gammaBar_ik psi^4 gammaBar_ik psi^-10 ABar^kl gammaBar_lj psi^4	<- note using the gammaBar to lower the ABar
A_ij = psi^-2 ABar_ij
ABar_ij = psi^2 A_ij

While we're here:
ABar_ij ABar^ij = psi^2 A_ij psi^10 A^ij
= psi^12 (K_ij - 1/3 gamma_ij K)(K^ij - 1/3 gamma^ij K)
= psi^12 (K_ij K^ij - 1/3 gamma_ij K K^ij - 1/3 gamma^ij K K_ij + 1/9 gamma_ij gamma^ij K^2)
= psi^12 (K_ij K^ij - 2/3 K^2 + 1/3 K^2)
= psi^12 (K_ij K^ij - 1/3 K^2)

Now for the conformal transverse-traceless decomposition of the extrinsic curvature
(Hemholtz / Hodge decomposition)
(Baumgarte & Shapiro p.67, Alcubierre p.94)

separate ABar^ij into two parts:
ABar^ij = ABarTT^ij + ABarL^ij

let ABarTT be the transverse part that is divergenceless:
DBar_j ABarTT^ij = 0

let ABarL be the longitudinal part:
ABarL^ij = DBar^i W^j + DBar^j W^i - 2/3 gammaBar^ij DBar_k W^K

Define L to be the longitudinal operator / vector gradient that produces a symmetric, traceless tensor, such that
ABarL^ij = (LBar W)^ij

We can now take the derivative of the ABar^ij equality to find
DBar_j ABar^ij = DBar_j (ABarTT^ij + ABarL^ij)
= DBar_j ABarTT^ij + DBar_j ABarL^ij
= 		0		  + DBar_j ABarL^ij
= DBar_j (LBar W)^ij
= DBar_j (DBar^i W^j + DBar^j W^i - 2/3 gammaBar^ij DBar_k W^k)
= DBar_j DBar^i W^j + DBar_j DBar^j W^i - 2/3 gammaBar^ij DBar_j DBar k W^k
= DBar_j DBar^i W^j + DBar_j DBar^j W^i - 2/3 DBar^i DBar_j W^j

Ricci identity (I think that's the name):
v^a_;dc - v^a_;cd = R^a_bcd v^b
Using prefix notation
D_c D_d v^a - D_d D_c v^a = R^a_bcd v^b
Contracting one pair of indexes:
D_j D_i v^j - D_i D_j v^j = R^j_kji v^k
D_j D_i v^j - D_i D_j v^j = R_ki v^k
Index gymnastics
D_j D^i v^j - D^i D_j v^j = R^i_k v^k
D_j D^i v^j = D^i D_j v^j + R^i_k v^k

So now we can say
DBar_j ABar^ij = DBar^i DBar_j W^j + RBar^i_j W^j + gammaBar_jk DBar_j DBar^k W^i - 2/3 DBar^i DBar_j W^j
= DBar^i DBar_j W^j + RBar^i_j W^j + DBar^2 W^i - 2/3 DBar^i DBar_j W^j
= DBar^2 W^i + 1/3 DBar^i DBar_j W^j + RBar^i_j W^j

Now we define a new operator: DeltaL, such that
DBar_j ABar^ij = (DeltaBarL W)^i

so we introduced two new operators on vector W^i:
(L W)^ij = D^i W^j + D^j W^i - 2/3 gamma^ij D_k W^k
...which looks strikingly similar to...
Lie_W(gamma^-1/3 gamma_ij) = D_i W_j + D_j W_i - 2/3 gamma_ij D_k W^k
...and...
(DeltaL W)^i = D_j (L W)^ij
...which is the vector Laplacian

Baumgarte & Shapiro go on to instruct
1) provide gammaBar_ij, K, and ABarTT^ij
2) solve for W^i using the following, which is the momentum constraint after all the conformal factor substitutions
(DeltaBarL W)^i = 2/3 psi^6 gammaBar^ij DBar_j K + 8 pi psi^10 S^i
note that it is in terms of psi.
Take the W^i, construct ABarLL^ij, combine with ABarTT^ij to get ABar^ij.
Lower with the gammaBar_ij metric to find the ABar_ij tensor density, and then...
3) solve for psi the following:
8 DBar^2 psi - psi RBar - 2/3 psi^5 K^2 + psi^-7 ABar_ij ABar^ij = -16 pi psi^5 rho

If we're in a vaccuum -- such that K_ij = 0, K = 0, rho = 0 -- then the momentum constraint becomes (Baumgarte & Shapiro p. 69)
(DeltaBarL W)^i = 0
and if we enforce gammaBar_ij = eta_ij then we get
partial^j partial_j W^i + 1/3 partial^i partial_j W^j = 0
Likewise (Alcubierre p. 97) the Hamiltonian constraint simplifies to
D^2 psi = 0
...and solutions to these are known as "Bowen York" solutions.

One such solution, which agrees Kerr black holes (Baumgarte & Shapiro, p.71):
ABarL^ij = (LBar W)^ij = 6/r^3 l^(i eBar^j)^kl J_k l_l
W^i = eBar^ijk X_j J_k
J^i = angular momentum density
eBar^ijk = Levi-Civita connection in gammaBar metric
l^i = x^i / r
X^i = l^i / r^2 = x^i / r^3

I'm going to skip ahead in Baumgarte & Shapiro to p.390 / ch.11 to the BSSN chapter.
Alcubierre has this back at p.94 which I skipped over.

It starts with
gamma_ij = exp(4 phi) gammaBar_ij
i.e.
gamma_ij = (exp(phi))^4 gammaBar_ij
so psi = exp(phi) and phi = ln(psi)

gamma^1/3 = psi^4
gamma = psi^12
ln(gamma) = 12 ln(psi)
ln(sqrt(gamma)) = 6 ln(psi) = 6 phi	<- converting from old to new time-derivative code
phi = ln(sqrt(gamma)) / 6 = ln(gamma) / 12

phi replaces ln(psi), check.

next we have
K_ij = exp(4phi) ATilde_ij + 1/3 gamma_ij K
which comes from
A_ij = K_ij - 1/3 gamma_ij K	<- conformal transverse-traceless
K_ij = A_ij + 1/3 gamma_ij K
so
A_ij = exp(4phi) ATilde_ij
replacing exp(phi) = psi

From there (Baumgarte & Shapiro p.64) we scale by an arbitrary power of psi:
A^ij = psi^alpha ABar^ij

A_ij = psi^4 ATilde_ij
this is the alternative of our original choice for alpha back at Baumgarte & Shapiro p.64 which said
ABar^ij = psi^10 A^ij
i.e.
ABar^ij = exp(10phi) A^ij
and subsequently by lowering via gammaBar_ij = psi^-4 gamma_ij
ABar_ij = gammaBar_ik ABar^kl gammaBar_lj = psi^10 psi^-4 gamma_ik A^kl psi^-4 gamma_lj
ABar_ij = psi^2 A_ij

the other (which Baumgarte & Shapiro introduce on page 120) is
ATilde^ij = psi^4 A^ij
i.e.
ATilde^ij = exp(4phi) A^ij
and is lowered via gammaBar, the conformal metric.  the only difference with ABar is that it is based on the solution alpha = 4 rather than alpha = 10
ATilde_ij = gammaBar_ik ATilde^kl gammaBar_lj = psi^4 psi^-4 gamma_ik A^kl psi^-4 gamma_lj
ATilde_ij = psi^-4 A_ij
ATilde_ij = exp(-4phi) A_ij
No mention of ATilde being raised or lowered by a metric other than the conformal metric anywhere, and Alcubierre p.83 affirms my lowered ATilde identity.
ironing things out (ATilde_ij in terms of the old ADM variables)
ATilde_ij = exp(-4phi) (K_ij - 1/3 gamma_ij K)
= exp(-4phi) K_ij - 1/3 exp(-4phi) gamma_ij K
= exp(-4phi) K_ij - 1/3 gammaBar_ij K

good ol ABar which I am ditching
ABar^ij = psi^10 A^ij
ABar_ij = psi^2 A_ij
ATilde^ij = psi^4 A^ij
ATilde_ij = psi^-4 A_ij
so exp(6phi) ATilde^ij = psi^6 ATilde^ij = psi^6 psi^4 A^ij = psi^10 A^ij = ABar^ij
ABar^ij = psi^10 psi^-4 ATilde^ij = psi^6 ATilde^ij

one extra for the record
traceless(A)_ij = A_ij - 1/3 gamma_ij A
traceless(A + B)_ij = (A + B)_ij - 1/3 gamma_ij trace(A + B)
traceless(A + B)_ij = A_ij + B_ij - 1/3 gamma_ij (gamma^kl (A_kl + B_kl))
traceless(A + B)_ij = A_ij + B_ij - 1/3 gamma_ij A - 1/3 gamma_ij B
traceless(A + B)_ij = traceless(A)_ij + traceless(B)_ij

and even one more
gammaBar_ij = gamma^-1/3 gamma_ij
such that
det(gammaBar_ij) = det(gamma^-1/3 gamma_ij)
= (gamma^-1/3)^3 det(gamma_ij)
= gamma^-1 gamma
= 1
makes me wonder, for our toy 2+1 and 1+1 sims,
should we be lowering the conformal / traceless fraction from 1/3 to 1/N?
and off of that, should we be lowering the EFE trace-reversed Ricci tensor
from 2/4 to 2/N?
probably not.

```