(I post this note proj_ELLIPTIC, proj_NS_ELLIPTIC, proj_NS_IVP and Cliff)

Hi all, 

This is an update report of my post to proj_ELLIPTIC last week.

1. Robin BC in Bernd's thorn_BAM_Elliptic
   
   Now, works fine. I tested 67**3, 131**3 with upto 8 processors. 
   
   I checked using Poisson equation for LaneEmden exact stars. 
   This Robin BC option gives close solution to the exact one. 
   Closer than that of my previous routine which fix BC as
   phi = - totalmass / r, and this is an expected result.    
   

2. Nakamura et al's BC for solving momentum constraint
   To say short, Nakamura et al uses a trial function which falls off
   at order O(1/r), and solves the next order O(1/r**2) or O(1/r**3)
   numerically. Their method is only applicable for conformally flat 
   case, but can apply for all situation with the same routine.   
 
   I coded up their method and compared with the previous ones, that is
   to solve momentum constraints for head-on neutron stars and compare the
   extrinsic curvature between 
     (i)   solving momentum constraints with Nakamura et al's BC
     (ii)  solving momentum constraints with simply set BC equals zero
     (iii) solving Poisson and time-derivative of Poisson equations,
           and use weak field approximation.  
   The list of equations are attached in the end of this file in LaTeX
   form 4 pages (but you may need to read their original reference.)
      
   All three methods give roughly the same order, the same shape of the
   extrinsic curvature. However, the tail part (anti-direction to the 
   heading on) of the extrinsic curvature are different, and of cource
   (i) and (iii) give the similar solution. 
   
   
3. thorn_pninit (Programs for Neutron star INITial data)
   The above routines are part of my CACTUS thorn pninit. 
   This thorn is available upon your request.  The recent status was
   posted in the project page proj_PN and proj_NS_IVP. 
   

4. 
HS> Progress of Theoretical Physics, Supplement 128 (1997 Dec)
HS> It has 400 pages, and $2.5 per page. 

It was $0.25 per page. 


PS. I will visit Japan from March 8 to April 13 in order to attend 
at least two conferences, one workshop, one meeting, and give more
than two seminars.  I will be able to contact by email, but not much
for coding.  If you have any technical requests or questions, please
send them before this Wednesday. 


Hisaaki SHINKAI                               
-----------------------------------------------------------
E-mail: shinkai@wurel.wustl.edu
URL:    http://wugrav.wustl.edu/People/HISAAKI/shinkai.html
Office: Dept. of Physics, Washington University
        Campus Box 1105, One Brookings Dr. 
        St.Louis, MO 63130-4899, USA
phone:  +1-314-935-4617  fax:  +1-314-935-6219  
-----------------------------------------------------------

% 19980228  HS pninit_44.tex
% to make complete, you need pninit_sa.tex and pninit_eq.tex
%%**********************************************************************
%23456789012345678901234567890123456789012345678901234567890123456789012
%00000000111111111122222222223333333333444444444455555555556666666666777
\documentstyle[11pt]{article}

\topmargin       0.0in
\oddsidemargin   0.0cm
\evensidemargin  0.0cm
\headsep         0.0in
\textwidth       6.5in
\textheight      9.0in

%-----------------------------------------------------------------------
%------------------------- macro for begin-eqs
\def\be{\begin{equation}}
\def\en{\end{equation}}
\def\bear{\begin{eqnarray}}
\def\enar{\end{eqnarray}}
%--------------->>>>>>>>>>>> commands for number of eqs. (1.1a)(1.1b)...
% following definition is for book or report style
%\renewcommand{\theequation}{\thechapter.\theenumi\alph{equation}}
\def\bera{ \setcounter{enumi}{\value{equation}}
           \addtocounter{enumi}{1}
           \setcounter{equation}{0}
           \renewcommand{\theequation}{\theenumi\alph{equation}}
           \begin{eqnarray} }
\def\enra{ \end{eqnarray}
           \setcounter{equation}{\value{enumi}}
           \renewcommand{\theequation}{\arabic{equation}}  }
%------------------------- macro for lists
\def\been{\begin{enumerate}}
\def\enen{\end{enumerate}}
\def\beit{\begin{itemize}}
\def\enit{\end{itemize}}
\def\dsp{\displaystyle}
\def\mova{\left( M \over a\right)}

%-------------------------
\newcommand{\square}
{\kern1pt\vbox{\hrule height 1.2pt\hbox{\vrule width1.2pt\hskip 3pt
\vbox{\vskip 6pt}\hskip 3pt\vrule width 0.6pt}\hrule height 0.6pt}
\kern1pt}
\newcommand{\gsim}{\stackrel{>}{\sim}}
\def\mn{_{\mu\nu}}
\def\non{\nonumber \\}
\def\nonn{\nonumber \\ &&}
\def\pl{\partial}
\def\ptl{\partial}
\def\half{{1\over 2}}
\def\disp{\displaystyle}
%-----------------------------------------------------------------------
\begin{document}

\begin{center}
{\LARGE\bf Metric and Extrinsic curvature extraction
\footnote{This is a part of the document 
{\it Cactus thorn pninit, version 1.4}.}
}
\end{center}
\begin{flushright}
19980228 Hisaaki Shinkai\footnote{
shinkai@wurel.wustl.edu, 
Physics Dept., Washington Univ., 
St. Louis, MO 63130-4899, USA.
}
\end{flushright}

%-----------------------------------------------------------------------
%\section{}
%-----------------------------------------------------------------------
Preparing metric is essential to construct an initial data for GR
evolution. CACTUS thorn pninit supports the following two methods. 

\subsection{Method 1: Use weak field approximation}
Suppose we have the solution of Poisson equations:
\bera
\Delta \psi  &=&4 \pi G \rho, \label{poissonNR1} \\
\Delta {\psi_t}  &=& -4 \pi G \pl_j (\rho  v^j), \label{poissonNR2}.
\enra
By the weak field approximation, metric is given by following relation.
\bera
g_{00}&=&-1-{2 \over c^2}\psi -{2 \over c^2}\psi^2  
+O({1\over c^4}) \\
g_{0i}&=&-{2 \over c^2}\psi^2  +O({1\over c^3}) \\
g_{ij}&=&\delta_{ij}(1-{2 \over c^2}\psi) 
-{4 G\over 5 c^5}Q_{ij}^{[3]}
+O({1\over c^7})
\enra
The lapse function $N$ is 
\be
N=\sqrt{-g_{00}}=1+{1 \over c^2}\psi +{1 \over c^2}\psi^2 
\en
The extrinsic curvature, then,
\be
K_{ij}= - {1\over 2 N} \pl_t g_{ij} = -
\delta_{ij}({1\over c^2}\psi_t) (1+ {2\over c^2} \psi
+ {2\over c^2} \psi^2)^{1/2}
\en
where $\psi_t$ is given by solving (\ref{poissonNR2}).

%-----------------------------------------------------------------------
\subsection{Method 2: Solve constraint equations}
Another method for preparing metric is to solve the 
constraint equations by inputing a density profile and velocity. 
The Hamiltonian constraint equation
\be
R+K^2-K_{ij}K^{ij}=2 \kappa\rho_H, \label{cham}
\en
where $\kappa=8\pi$, $R$ is the 3-dimensional Ricci scalar curvature and 
$\rho_H$ is the energy density  measured by the normal line
observer $\rho_H\equiv n^\mu n^\nu T_{\mu\nu}$, where $n^\mu$ is the
unit timelike 4-vector. 

The momentum constraint equation 
\be
D_j(K^j_{~i}-\delta^j_{~i} K)=\kappa J_i \label{cmom}
\en
where $D_j$ is the 3-d covariant derivative, and $J_i$ is the 
momentum density $J_i\equiv - h_i^{~\mu} n^\nu T_{\mu\nu}$, where
$h_{\mu\nu}$ is the projection tensor $h_{\mu\nu}=g_{\mu\nu}+n_\mu n_\nu$.

We assume the perfect fluid stress-energy tensor,
\be
T_{\mu\nu}=(\rho+\rho\varepsilon+p) u_\mu u_\nu + p g_{\mu\nu}
\en
where $\rho, \varepsilon$ and $p$ are the proper mass density, the
specific internal energy and the pressure, respectively, and $u_\mu$
is the 4-velocity of the fluid. 

The actual relations with input quantities we take are
\bear
\rho_H&=&\rho+ \rho\varepsilon, \\
v^i&=&{u^i \over u^0}={\alpha J^i \over p+\rho_H}-\beta^i.
\enar
%---------- ---------- ---------- ---------- ---------- ---------- ----- 
\subsubsection{Solve the Hamiltonian constraint }

\paragraph{York-O'Murchadha's conformal approach }  %-----------

Defining conformal factor $\psi$ and setting
\be
\gamma_{ij} = \psi^4 \hat{\gamma}_{ij},  ~~~ 
\gamma^{ij} = \psi^{-4} \hat{\gamma}^{ij},  ~~~ 
\rho_B = \psi^{-n} \hat{\rho},
\label{york_conf}
\en
where $n=6$ or 8 ($n>5$ is required from uniqueness of the solution).
The transformation of the scalar curvature is
\be
R= \phi^{-4} \hat{R} - 8 \phi^{-5} \hat{\Delta} \phi.
\en
The constraint, then,  become
\bear
~8~^{(3)\!}\hat{\Delta}
\psi&=&~^{(3)\!}\hat{R}\psi-\hat{K}_{ij}^{TF}
\hat{K}^{ij}_{TF}\psi^{-7} +2[{1 \over 3}
(\mbox{tr} \hat{K})^2 - \Lambda]\psi^5
-2\kappa \rho \psi^{5}
%-2\hat{\rho}\psi^{-1} 
\nonumber \\
&=&~^{(3)\!}\hat{R}\psi-\hat{K}_{ij}^{TF}
\hat{K}^{ij}_{TF}\psi^{-7} +2[{1 \over 3}
(\mbox{tr} \hat{K})^2 - \Lambda]\psi^5
-2\kappa \hat{\rho}\psi^{5-n}
%-2\hat{\rho}\psi^{-1}
\label{inithamilt} 
%\\
%\hat{D}_j\hat{K}^{ij}_{TF}&=&{2 \over 3}\psi^6\hat{D}^i \mbox{tr}
%\hat{K} +\kappa\hat{J}^i \label{initmoment}
\enar
where $^{(3)\!}\hat{\Delta}$ is
the 3-dimensional Laplacian defined by $\hat{\gamma}_{ij}$, and
$~^{(3)\!}\hat{R}$ is  the scalar curvature.
In this approach, $\hat{\gamma}_{ij}, \mbox{tr}\hat{K}$ and the
transverse-traceless(TT) part of the extrinsic curvature
$\hat{K}^{TT}_{ij}$ (and $\hat{\rho}$ and $\hat{J}_i$ when we
include matter) are  left to our choice.


\paragraph{technical details}  %-----------
\begin{itemize}
\item input (truncated) TOV solution $\rho(r)$ as $\hat{\rho}=\rho_H$. 
\item $n=8$ is used. 
\item conformal flat spacetime is used as a trial 3-metric. 
\item 
The hamiltonian constraint (\ref{inithamilt}) is linearized around
some solution $\psi_0$, where the true solution is given by 
 $ \psi=\psi_0 + \delta \psi$, as:
\bear
8 \Delta\psi &= &E\psi+F \psi^{-7} + G \psi^5 + H \psi^{-3}
+ I \psi^{-1}  + J
\nonumber \\
&=& E \psi + F (\psi_0 + \delta \psi)^{-7} 
+ G  (\psi_0 + \delta \psi)^5 + H (\psi_0 + \delta \psi)^{-3}
+ I (\psi_0 + \delta \psi)^{-1} + J
\nonumber \\
&=& E \psi  
+ F (\psi_0^{-8} -7 \psi_0^{-6} \delta\psi)
+ G (\psi_0^5 + 5 \psi_0^{4}  \delta \psi) 
+ H (\psi_0^{-3} -3  \psi_0^{-4} \delta\psi)
+ I (\psi_0^{-1} -   \psi_0^{-2} \delta\psi)
+ J 
\nonumber \\
&=&  
[
E -7 F \psi_0^{-8} + 5 G \psi_0^{4} - 3 H  \psi_0^{-4} - 2 I  \psi_0^{-2} 
] \psi 
\nonumber \\
&& + 
[
8 F \psi_0^{-7} -4  G  \psi_0^5 + 4 H  \psi_0^{-3} + 2 I \psi_0^{-1}  +J
]
\enar
This technique is explained in \cite{anninos89}.
For example, if we choose the conformal metric is flat 
$\hat{\gamma}_{ij}=\eta_{ij}$
(then $\hat{R}=0$), $n=8$, and $K=0$ at initial slice, then 
\be
8 \Delta \psi = F \psi^{-7} + H \psi^{-3} 
=  
[ -7 F \psi_0^{-8} - 3 H  \psi_0^{-4} 
] \psi 
 + 
[
8 F \psi_0^{-7}  + 4 H \psi_0^{-3} ]
\en
We initially start with $\psi_0=1$ and keep iteration 
by updating  $\psi_0=\psi_{solved}$
until it converges. 

\item boundary is assumed asymptotically flat solution, 
\be \psi= 1 + {M_{total} \over 2r} \en
\item after solving the Hamiltonian eq., we need to convert observer normal
$\rho_B$ into
matter density $\rho$ as:
 \be \rho_B = \rho+\rho \varepsilon. \en
\end{itemize}

%---------- ---------- ---------- ---------- ---------- ---------- ----- 
\subsubsection{Solve the momentum constraint }
In the case we put velocity terms in the initial configuration, we
need to solve the momentum constraint. 

\paragraph{York-O'Murchadha's conformal approach}   %-----------
Begin from decomposing extrinsic curvature $K_{ij}$ as
\be
K_{ij}=K_{ij}^{TF}+{1\over 3} \gamma_{ij} K
= K_{ij}^{TT} + ({\boldmath l}W)_{ij} +{1\over 3} \gamma_{ij} K 
\en
where $TF, TT$ denote trace-free and transverse-traceless,
respectively. 
$({\boldmath l}~)_{ij}$ is an expression of longitudinal
part, which is defined as 
\be
 ({\boldmath l}W)_{ij} = D_i W_j + D_j W_i -{2\over 3} \gamma_{ij}
D_k W^k 
\en
and
\bear
D_j ({\boldmath l}W)^{ij}
\equiv (\Delta_{\boldmath l}  W)^i
 &=& D_k D_k W^i +{1\over 3} D^i D_j W^j
+R^i_{~j} W^j \\
 &\equiv& (\Delta  W)^i +{1\over 3} D^i D_j W^j
+R^i_{~j} W^j. 
\enar


Defining conformal factor $\psi$ as in (\ref{york_conf}), and setting
\bera
K_{ij}^{TF} &=& \psi^{-2}\hat{K}_{ij}^{TF},  ~~~
K^{ij}_{TF} = \psi^{-10}\hat{K}^{ij}_{TF}, ~~~
K^{i}_{~j~TF} = \psi^{-6}\hat{K}^{i}_{~j~TF}
\\
J^i &=& \psi^{-10} \hat{J}^i, ~~~ 
J_i = \psi^{-6} \hat{J}_i,
\label{york_conf2}
\enra
The constraint is, then, 
\bera
&&\hat{D}_j\hat{K}^{ij}_{TF}
=\hat{D}_j({\boldmath l}W)^{ij} =
{2 \over 3}\psi^6\hat{D}^i {K} +\kappa\hat{J}^i \\
{\mbox or}~~~&&
 (\Delta  W)^i +{1\over 3} D^i D_j W^j
+R^i_{~j} W^j - {2 \over 3}\psi^6\hat{D}^i {K} =  \kappa\hat{J}^i
\label{initmoment}
\enra

After solving this equation, 
together with assumption of $K$ and $K_{ij}^{TT}$, 
the extrinsic curvature will be given by 
\be
K_{ij}=[K_{ij}^{TT} + ({\boldmath l}W)_{ij}]\psi^{-2} 
+{1\over 3}\psi^4 \hat{\gamma}_{ij}K
\en
We assume $K=const.$ to decouple two constraints, and also assume 
 $K_{ij}^{TT}=0$, 

\paragraph{technical details}   %-----------
We solve (\ref{initmoment}) in the conformally flat
$\hat{\gamma}_{ij}=\eta_{ij}$ background, 
with $K=0$ slice.  In this situation,  (\ref{initmoment}) becomes 
\be
{\Delta} W_i +{1\over 3} \ptl_i \ptl_j W_j =\kappa \hat{J}_i.
\label{initmomentflat}
\en

Two ways to solve (\ref{initmomentflat}) are described in \cite{sup97}.
\beit
\item Define a vector potential $A_i$ and a scalar potential $\hat{\psi}$ such 
that they satisfy
\bera
\Delta A_i &=& 8 \pi \hat{J}_i  \label{method1-A} \\
\Delta \hat{\psi} &=& 8 \pi \hat{J}_i x^i. \label{method1-psi} 
\enra
Solve these  $A_i$ and  $\hat{\psi}$, then $W_i$ is given by 
\be
W_i=A_i+{1\over 8} \ptl_i (\hat{\psi} - A_k x^k).
\en
\item Define $\chi\equiv \ptl_k W_k$, then (\ref{initmomentflat}) 
decouples to
\bera
{\Delta} \chi &=& {3 \over 4} \kappa \ptl_k \hat{J}_k  \label{method2-x} \\
{\Delta} W_i &=& \kappa \hat{J}_i - {1\over 3} \ptl_i \chi.
\label{method2-W} 
\enra
\enit
\paragraph{outer BC}   %-----------
We follow Oohara-Nakamura-Shibata\cite{sup97} to set BC for solving 
(\ref{method1-A}) - (\ref{method2-W}). Since the falling off behavior
of $W_i$ is different at the appearance of the linear momentum of the
system, they show a technique for renormalizing background falling 
off terms. Here, only used equations in the program are shown. 

Instead of solving (\ref{method1-A}) and (\ref{method1-psi}), we
solve 
\bera
\Delta \bar{A}_i &=& 8 \pi \bar{J}_i  \label{method1-A2} \\
\Delta \bar{\psi} &=& 8 \pi \bar{J}_i x^i. \label{method1-psi2} 
\enra
with the boundary conditions as
\bera
\bar{A}_i &=& O ({1\over r^3}) \\
\bar{\psi} &=& O ({1\over r^2}) 
\enra
where 
\be
\bar{J}_i=\hat{J}_i - M_{ij} x^j S^{(2)}.
\en
($M_{ij}$ and $S^{(2)}$ are defined below.)
The solution of (\ref{method1-A}) and (\ref{method1-psi})
are given by 
\bera
\bar{A}_i &=& \bar{A}_i +{A}_i^{(0)}\\
\hat{\psi}&=& \bar{\psi} + {\psi}^{(0)}
\enra
where ${A}_i^{(0)}$ and ${\psi}^{(0)}$ are first order falling-off
solutions of this system and are given by 
\bera
%\tilde{J}_i^{(0)} &\equiv& P_i S^{(1)} + M_{ij} x^j  S^{(2)} \\
{A}_i^{(0)} &\equiv& -2 (P_i \phi^{(1)} + M_{ij} x^j  \phi^{(2)}) \\
{\psi}^{(0)} &\equiv& -2 (P_i x^i \psi^{(1)} + M_{ij} x^i x^j  \psi^{(2)} 
+ M_{ii}  \psi^{(3)})
\enra
where 
we defined moments of momentum density $\hat{J}_i$. 
\bera
P_i&=&\int \hat{J}_i d^3 x \\
M_{ij}&=&\int \hat{J}_i x^j d^3 x 
\enra
and
\bera
%S^{(1)} &=& {15 \over 8\pi a^3} (1-\xi^2) 
%~~if~~r \le a, ~~~ 0 ~~if~~r > a \\
S^{(2)} &=&  {105 \over 8\pi a^5} (1-\xi^2) 
~~if~~r \le a, ~~~ 0 ~~if~~r > a\\
\phi^{(1)} &=&{1 \over 8 a}(15-10 \xi^2+3\xi^4)  
~~if~~r \le a, ~~~ {1\over r} ~~if~~r > a \\
\phi^{(2)} &=&{1 \over 8 a^3}(35-42 \xi^2+15\xi^4) 
~~if~~r \le a, ~~~ {1\over r^3} ~~if~~r > a \\
\psi^{(1)} &=& {1 \over 56 a}(35-42\xi^2+15\xi^4) 
~~if~~r \le a, ~~~ {a^2 \over 7 r^3} ~~if~~r > a \\
\psi^{(2)} &=& {1 \over 24 a^3}(63-90\xi^2+35\xi^4)
~~if~~r \le a, ~~~ {a^2 \over 3 r^5} ~~if~~r > a \\
\psi^{(3)} &=& {1 \over 72 a}(105-63\xi^2+27\xi^4-5\xi^6)
~~if~~r \le a, ~~~ {1\over r}-{a^2 \over 9 r^3} ~~if~~r > a 
\enra
where $\xi=r/a$, and $a$ is a constant, which we set 1.5 times binary's 
separation. 




\baselineskip .2in
\begin{thebibliography}{99}
\bibitem{yorkinitial}
N. O'Murchadha and J.W. York, Jr., Phys. Rev. {\bf D 10}, 428 (1974).
\bibitem{anninos89}
P. Anninos, J. Centrella and R. Matzner, Phys. Rev. {\bf D 39}, 2155 (1989).
\bibitem{sup97}
K. Oohara, T. Nakamura and M. Shibata, {\it A way to 3D NR},
Prog. Theor. Phys. Suppl. {\bf 128}, 183 (1997). 

\end{thebibliography}



%=======================================================================
\begin{verbatim}
Hisaaki SHINKAI
-----------------------------------------------------------
E-mail: shinkai@wurel.wustl.edu
URL:    http://wugrav.wustl.edu/People/HISAAKI/shinkai.html
Office: Dept. of Physics, Washington University
        Campus Box 1105, One Brookings Dr.
        St.Louis, MO 63130-4899, USA
phone:  +1-314-935-4617  fax:  +1-314-935-6219
-----------------------------------------------------------
\end{verbatim}
\end{document}

Created for the WUGRAV CoCoBoard Page Projects Page.

Created by The CoCoBoard.