\documentclass[10pt]{beamer} \input{mystyle_beamer.tex} \usetikzlibrary{arrows.meta} \usetikzlibrary{decorations.pathmorphing} \usetikzlibrary{calc} \usetikzlibrary{positioning,fit,backgrounds} \usetikzlibrary{external} \tikzexternalize[prefix=tikzcache/] \usepackage{appendixnumberbeamer} \tikzset{zigzag/.style={decorate, decoration=zigzag}} \def \L {2.} \newcommand{\lapse}{\alpha} % \newcommand{\shiftu}[1]{\ensuremath{\tensor{\beta}{^{#1}}}} % \newcommand{\shiftd}[1]{\ensuremath{\tensor{\beta}{_{#1}}}} \newcommand{\shift}[1]{\tensor{\beta}{#1}} \newcommand{\pt}{\Partial{_t}} % \newcommand{\csmetric}[1]{\tensor{\tilde{\gamma}}{#1}} \NewDocumentCommand{\smetric}{ s m }{% s=可选星号，m=索引们 \IfBooleanTF{#1} {\tensor{\tilde{\gamma}}{#2}} {\tensor{\gamma}{#2}} } \NewDocumentCommand{\Ktf}{ s m }{% s=可选星号，m=索引们 \IfBooleanTF{#1} {\tensor{\tilde{A}}{#2}} {\tensor{A}{#2}} } \NewDocumentCommand{\Gam}{ s m }{% s=可选星号，m=索引们 \IfBooleanTF{#1} {\tensor{\tilde{\Gamma}}{#2}} {\tensor{\Gamma}{#2}} } \NewDocumentCommand{\Gamd}{ s m }{% s=可选星号，m=索引们 \IfBooleanTF{#1} {\tensor{(\tilde{\Gamma}_\text{d})}{#2}} {\tensor{(\Gamma_\text{d})}{#2}} } \NewDocumentCommand{\Rictensor}{ s m }{% s=可选星号，m=索引们 \IfBooleanTF{#1} {\tensor{\tilde{R}}{#2}} {\tensor{R}{#2}} } \RenewDocumentCommand{\spaceD}{ s m }{% s=可选星号，m=索引们 \IfBooleanTF{#1} {\tensor{\tilde{D}}{#2}} {\tensor{D}{#2}} } \NewDocumentCommand{\constraintZ}{ s m }{% s=可选星号，m=索引们 \IfBooleanTF{#1} {\tensor{\tilde{Z}}{#2}} {\tensor{Z}{#2}} } \NewDocumentCommand{\constraintM}{ s m }{% s=可选星号，m=索引们 \IfBooleanTF{#1} {\tensor{\tilde{M}}{#2}} {\tensor{M}{#2}} } \NewDocumentCommand{\dlapse}{ s m }{% s=可选星号，m=索引们 \IfBooleanTF{#1} {\tensor{\tilde{\lapse}}{#2}} {\tensor{\lapse}{#2}} } \NewDocumentCommand{\dshift}{ s m }{% s=可选星号，m=索引们 \IfBooleanTF{#1} {\tensor{\tilde{\beta}}{#2}} {\tensor{\beta}{#2}} } \NewDocumentCommand{\dchi}{ s m }{% s=可选星号，m=索引们 \IfBooleanTF{#1} {\tensor{\tilde{\phi}}{#2}} {\tensor{\phi}{#2}} } \NewDocumentCommand{\dsmetric}{ s m }{% s=可选星号，m=索引们 \IfBooleanTF{#1} {\tensor{\tilde{\psi}}{#2}} {\tensor{\psi}{#2}} } \NewDocumentCommand{\constraintA}{ s m }{% s=可选星号，m=索引们 \IfBooleanTF{#1} {\tensor{\tilde{\mathcal{A}}}{#2}} {\tensor{\mathcal{A}}{#2}} } \NewDocumentCommand{\constraintB}{ s m }{% s=可选星号，m=索引们 \IfBooleanTF{#1} {\tensor{\tilde{\mathcal{B}}}{#2}} {\tensor{\mathcal{B}}{#2}} } \NewDocumentCommand{\constraintC}{ s m }{% s=可选星号，m=索引们 \IfBooleanTF{#1} {\tensor{\tilde{\mathcal{C}}}{#2}} {\tensor{\mathcal{C}}{#2}} } \NewDocumentCommand{\constraintD}{ s m }{% s=可选星号，m=索引们 \IfBooleanTF{#1} {\tensor{\tilde{\mathcal{D}}}{#2}} {\tensor{\mathcal{D}}{#2}} } \newcommand{\oc}[1]{\mathring{#1}} \newcommand{\olapse}{\oc{\lapse}} % \newcommand{\shiftu}[1]{\ensuremath{\tensor{\beta}{^{#1}}}} % \newcommand{\shiftd}[1]{\ensuremath{\tensor{\beta}{_{#1}}}} \newcommand{\oshift}[1]{\tensor{{\oc{\beta}}}{#1}} % \newcommand{\pt}{\Partial{_t}} % % \newcommand{\csmetric}[1]{\tensor{\tilde{\gamma}}{#1}} \newcommand{\ochi}{\oc{\chi}} \NewDocumentCommand{\osmetric}{ s m }{% s=可选星号，m=索引们 \IfBooleanTF{#1} {\tensor{{\oc{\tilde{\gamma}}}}{#2}} {\tensor{{\oc{\gamma}}}{#2}} } \NewDocumentCommand{\oKtf}{ s m }{% s=可选星号，m=索引们 \IfBooleanTF{#1} {\tensor{{\oc{\tilde{A}}}}{#2}} {\tensor{{\oc{A}}}{#2}} } % \NewDocumentCommand{\Gam}{ s m }{% s=可选星号，m=索引们 % \IfBooleanTF{#1} % {\tensor{\tilde{\Gamma}}{#2}} % {\tensor{\Gamma}{#2}} % } % \NewDocumentCommand{\Gamd}{ s m }{% s=可选星号，m=索引们 % \IfBooleanTF{#1} % {\tensor{(\tilde{\Gamma}_\text{d})}{#2}} % {\tensor{(\Gamma_\text{d})}{#2}} % } % \NewDocumentCommand{\Rictensor}{ s m }{% s=可选星号，m=索引们 % \IfBooleanTF{#1} % {\tensor{\tilde{R}}{#2}} % {\tensor{R}{#2}} % } % \RenewDocumentCommand{\spaceD}{ s m }{% s=可选星号，m=索引们 % \IfBooleanTF{#1} % {\tensor{\tilde{D}}{#2}} % {\tensor{D}{#2}} % } % \NewDocumentCommand{\constraintZ}{ s m }{% s=可选星号，m=索引们 % \IfBooleanTF{#1} % {\tensor{\tilde{Z}}{#2}} % {\tensor{Z}{#2}} % } % \NewDocumentCommand{\constraintM}{ s m }{% s=可选星号，m=索引们 % \IfBooleanTF{#1} % {\tensor{\tilde{M}}{#2}} % {\tensor{M}{#2}} % } % \NewDocumentCommand{\dlapse}{ s m }{% s=可选星号，m=索引们 % \IfBooleanTF{#1} % {\tensor{\tilde{\lapse}}{#2}} % {\tensor{\lapse}{#2}} % } % \NewDocumentCommand{\dshift}{ s m }{% s=可选星号，m=索引们 % \IfBooleanTF{#1} % {\tensor{\tilde{\beta}}{#2}} % {\tensor{\beta}{#2}} % } % \NewDocumentCommand{\dchi}{ s m }{% s=可选星号，m=索引们 % \IfBooleanTF{#1} % {\tensor{\tilde{\phi}}{#2}} % {\tensor{\phi}{#2}} % } \NewDocumentCommand{\odsmetric}{ s m }{% s=可选星号，m=索引们 \IfBooleanTF{#1} {\tensor{\oc{\tilde{\psi}}{}}{#2}} {\tensor{\oc{\psi}{}}{#2}} } % \NewDocumentCommand{\constraintA}{ s m }{% s=可选星号，m=索引们 % \IfBooleanTF{#1} % {\tensor{\tilde{\mathcal{A}}}{#2}} % {\tensor{\mathcal{A}}{#2}} % } % \NewDocumentCommand{\constraintB}{ s m }{% s=可选星号，m=索引们 % \IfBooleanTF{#1} % {\tensor{\tilde{\mathcal{B}}}{#2}} % {\tensor{\mathcal{B}}{#2}} % } % \NewDocumentCommand{\constraintC}{ s m }{% s=可选星号，m=索引们 % \IfBooleanTF{#1} % {\tensor{\tilde{\mathcal{C}}}{#2}} % {\tensor{\mathcal{C}}{#2}} % } % \NewDocumentCommand{\constraintD}{ s m }{% s=可选星号，m=索引们 % \IfBooleanTF{#1} % {\tensor{\tilde{\mathcal{D}}}{#2}} % {\tensor{\mathcal{D}}{#2}} % } \title{Towards Moving Puncture Simulations in a Discontinuous Galerkin Framework} \subtitle{A Progress Report} \author{Yingjie Wang} \institute{FAU} \addbibresource{ref.bib} \begin{document} \maketitle \begin{frame}{Contents} \tableofcontents \end{frame} \begin{frame}{Outline of the Project} \section{Introduction} \begin{enumerate} \item A first order evolution system, to be named FOZ4c \item An implementation in our discontinuous Galerkin code \texttt{nmesh} \item Simulations of binary black holes with large mass ratio with the moving puncture gauge \end{enumerate} \end{frame} \begin{frame}{Final Goal: binary black holes with large mass ratio} \subsection{binary black holes with large mass ratio} Why binary black holes? \begin{itemize} \item Main sources of gravitational waves \item Simplest two-body problem in general relativity \item Connecting strong field GR with post-Newtonian approximations \end{itemize} Why large mass ratio? \begin{itemize} \item LISA is able to detect extreme mass ratio inspirals (EMRIs) \item Compare with perturbation theory \end{itemize} \end{frame} \begin{frame}{Methods on evolving black holes} \subsection{Methods on evolving black holes} There are two main methods to evolve black holes in numerical relativity: \begin{itemize} \item Excision method: excise the black hole interior from the computational domain, and impose boundary conditions on the excision surface \item Moving puncture method: evolve the black hole as a puncture, and use a suitable gauge condition to avoid the singularity. \end{itemize} \end{frame} \begin{frame}{Methods on evolving black holes: excision} \begin{figure} \centering \begin{tikzpicture}[>=Latex, line cap=round, line join=round] % causal diamond \draw[thick,red,zigzag] (-\L,\L) coordinate(stl) -- (\L,\L) coordinate (str); \draw[thick,black] (\L,-\L) coordinate (sbr) -- (0,0) coordinate (bif) -- (stl); \draw[thick,black,fill=blue, fill opacity=0.2,text opacity=1] (bif) -- (str) -- (2*\L,0) node[right] (io) {$i^0$} -- (sbr); % null labels \draw[black] (1.4*\L,0.7*\L) node[right] (scrip) {$\mathcal{I}^+$} (1.5*\L,-0.6*\L) node[right] (scrip) {$\mathcal{I}^-$} (0.2*\L,-0.6*\L) node[right] (scrip) {$\mathcal{H}^-$} (0.5*\L,0.85*\L) node[right] (scrip) {$\mathcal{H}^+$}; % singularity label \draw[thick,red,<-] (0,1.05*\L) -- (0,1.2*\L) node[above] {\color{red} singularity}; % % Scwharzschild surface % \draw[thick,blue] (bif) .. controls (1.*\L,-0.35*\L) .. (2*\L,0); % \draw[thick,blue,<-] (1.75*\L,-0.1*\L) -- (1.9*\L,-0.5*\L) % -- (2*\L,-0.5*\L) node[right,align=left] % {$t=$ constant\\in Schwarzschild\\coordinates}; % excision surface \draw[thick,dashed,red] (-0.3*\L,0.3*\L) -- (0.4*\L,\L); \draw[thick,red,<-] (-0.33*\L,0.3*\L) -- (-0.5*\L,0.26*\L) node[left,align=right] {excision\\surface}; % Kerr-Schild surface \draw[green,thick] (0.325*\L,0.325*\L) .. controls (\L,0) .. (2*\L,0); \draw[green,dashed,thick] (0.325*\L,0.325*\L) -- (-0.051*\L,0.5*\L); % Kerr-Schild label \draw[green,thick,<-] (0.95*\L,0.15*\L) -- (1.2*\L,0.5*\L) -- (2*\L,0.5*\L) node[right,align=left] {time slice}; \end{tikzpicture} \caption{excision on a single Schwarzschild black hole} \end{figure} \end{frame} \begin{frame}{Methods on evolving black holes: excision} \begin{figure} \centering \includegraphics[width=0.7\textwidth]{imgs/black_hole_excision_mesh.png} \caption{An example mesh using excision method from \cite{Hemberger_2013}.} \end{figure} \end{frame} \begin{frame}{Methods on evolving black holes: moving puncture} In the isotropic coordinate for Schwarzschild black hole, the spatial metric \begin{equation} \dd{l^2} = \left( 1+\frac{M}{2r} \right)^4 (\dd{r^2} + r^2 \dd{\Omega^2}) \end{equation} is conformally flat, we can define a conformal metric $\tensor{\tilde{\gamma}}{_i_j}$ by $\tensor{\tilde{\gamma}}{_i_j} = \psi^{-4} \tensor{\gamma}{_i_j}$, where \begin{equation} \psi = 1 + \frac{M}{2r}. \end{equation} General $N$ black holes puncture initial data: \begin{equation} \psi = \sum_{i=1}^N \frac{M_i}{2|\myvec{r}-\myvec{r}_i|} + (\text{regular part}), \end{equation} again, if we factor out the conformal factor $\psi^4$, the remaining conformal metric $\tensor{\tilde{\gamma}}{_i_j}$ is regular everywhere. \end{frame} \begin{frame}{Methods on evolving black holes: moving puncture} In the moving puncture method, we won't cut out the black hole singularity, but choose a gauge condition that makes the singularity invisible to the numerical evolution. \begin{figure} \centering \includegraphics[width=0.7\textwidth, trim=4bp 4bp 4bp 4bp, clip]{imgs/moving_puncture_penrose_1.png} \caption{The numerical time slices in the moving puncture method in the Penrose diagram, from \cite{Hannam_2008}.} \end{figure} \end{frame} \begin{frame}{Methods on evolving black holes: moving puncture} \begin{figure} \centering \includegraphics[width=0.42\textwidth, trim=5bp 4bp 4bp 4bp, clip]{imgs/moving_puncture_init.png} \includegraphics[width=0.45\textwidth, trim=4bp 4bp 4bp 4bp, clip]{imgs/moving_puncture_stable.png} \caption{The embedding of the initial slice (left) and the final slice (right) in the moving puncture method, from \cite{Hannam_2008}.} \end{figure} \end{frame} \begin{frame}{Discontinuous Galerkin method and \texttt{nmesh}} \subsection{Discontinuous Galerkin method and \texttt{nmesh}} There are two main numerical methods in numerical relativity: \begin{itemize} \item Finite Difference / Finite Volume \begin{itemize} \item Easy to implement, and mature codes available \end{itemize} \item Spectral / Discontinuous Galerkin \begin{itemize} \item High accuracy for smooth solutions \end{itemize} \end{itemize} \end{frame} \begin{frame}{Discontinuous Galerkin method and \texttt{nmesh}} Core idea of DG method: if we have a first-order PDE system: \begin{equation} \pdvt{u} + A^i \pdv{u}{x^i} = S, \end{equation} then for some test functions $\{v_a\}$, we have \begin{equation} \left(v_a, \pdvt{u}\right) + \left(v_a, A^i \pdv{u}{x^i}\right) = (v_a, S), \end{equation} which reduce to a linear system of $\{ u_a = (v_a, u)\}$ after integrating by parts. Boundary terms are replaced by numerical fluxes. In \texttt{nmesh}, we use Lagrange polynomials over Gauss-Legendre points on each element. \begin{center} \begin{tikzpicture}[scale=4] \draw[thick, blue] (0,0) -- (1,0); \draw[thick, red] (1,0) -- (2,0); \pgfmathsetmacro{\a}{0.5 - 0.5*sqrt(3/7)} \pgfmathsetmacro{\b}{0.5 + 0.5*sqrt(3/7)} \pgfmathsetmacro{\c}{1.5 - 0.5*sqrt(3/7)} \pgfmathsetmacro{\d}{1.5 + 0.5*sqrt(3/7)} \foreach \x in {0,\a,0.5,\b,1} \fill[blue] (\x,0) circle (0.018); \foreach \x in {1,\c,1.5,\d,2} \fill[red] (\x,0) circle (0.018); \fill[blue] (1,0) circle (0.022); \fill[red] (1,0) circle (0.013); \end{tikzpicture} \end{center} Information exchange only happens at the boundaries of elements through numerical fluxes, which makes it efficient for parallel computing. \end{frame} \begin{frame}{Discontinuous Galerkin method and \texttt{nmesh}} \texttt{nmesh} has two features that are useful for our project: \begin{itemize} \item Adaptive mesh refinement (AMR) \item DG + FD/FV dynamically switching \end{itemize} \begin{figure} \centering \includegraphics[width=0.65\textwidth]{imgs/MPA1_W-9sn12l5_GRHD_D_t0400.pdf} \caption{The mesh in a binary star emulation with \texttt{nmesh}} \end{figure} \end{frame} \begin{frame}{Moving puncture evolution in \texttt{nmesh}?} \subsection{Moving puncture evolution in \texttt{nmesh}?} There are many $3+1$ formulations of Einstein's equations, \begin{itemize} \item We want DG method: \begin{itemize} \item a {\color{blue}first order} formulation is needed, like the generalized harmonic (GH) formulation. \end{itemize} \item We want moving puncture evolution: \begin{itemize} \item the {\color{blue}moving puncture gauge condition} is needed, and not compatible with the GH formulation. \item People usually use BSSN, Z4c or CCZ4 formulation for moving puncture evolutions, but they are second order in space, and not directly suitable for DG method. \end{itemize} \item We want stable evolution: \begin{itemize} \item a {\color{blue}strongly hyperbolic} formulation is needed. \item {\color{blue} constraint damping} is needed. \end{itemize} \end{itemize} $\implies$ we need to find a first order, strongly hyperbolic formulation of Einstein's equations with constraint damping and compatible with the moving puncture gauge condition. \end{frame} \begin{frame}{Searching for a first order Z4c}{Why Z4c?} \section{Searching for a first order Z4c} \subsection{Why Z4c?} % Let's look at the relationship between different mainly used formulations of Einstein's equations: \begin{figure} \centering \begin{tikzpicture}[ scale=0.5, >=Latex, node distance=8mm and 10mm, every node/.style={font=\small}, sys/.style={ draw, rounded corners, align=center, minimum width=15mm, minimum height=7mm, inner sep=2pt }, std/.style={sys, fill=blue!8}, fo/.style={sys, fill=red!10}, rel/.style={->, thick}, weak/.style={->, semithick, dashed}, equv/.style={<->, thick, dotted} ] \node[std] (bssn) {BSSN}; \node[std, right=of bssn] (z4) {Z4}; \node[std, right=of z4] (z4c) {Z4c}; \node[std, right=of z4c] (ccz4) {CCZ4}; % \node[fo, below=10mm of bssn] (fobssn) {FOBSSN}; \node[fo, below=8mm of z4] (gh) {GH}; % \node[fo, below=10mm of ccz4] (foccz4) {FOCCZ4}; \draw[rel] (z4) -- (z4c); \draw[rel] (z4) to[bend left=18] (ccz4); \draw[weak] (bssn) -- (z4c); % \draw[rel] (bssn) -- (fobssn); % \draw[rel] (ccz4) -- (foccz4); \draw[equv] (z4) -- (gh); % \node[draw=none, above=1mm of gh, font=\scriptsize] {full first order}; \end{tikzpicture} \caption{Commonly used systems.\footnote{Blue: second order systems; Red: first order systems. Solid arrows: direct modifications; Dashed arrows: indirect relation. Bidirectional arrows: equivalence.}} \end{figure} \vspace{-0.5cm} \begin{itemize} \item Z4 and GH evolve the physical metric, and thus not compatible with the moving puncture method. \item The rests use the conformal metric. \begin{itemize} \item Z4c damps the constraints better than BSSN, as the constraint violations propagate at the speed of light and thus move out. \item Unstable in CCZ4 for black hole spacetime was reported with some parameters. \end{itemize} \end{itemize} \end{frame} \begin{frame}{First order Z4c?} \subsection{First order Z4c?} \begin{figure} \centering \begin{tikzpicture}[ scale=0.5, >=Latex, node distance=8mm and 10mm, every node/.style={font=\small}, sys/.style={ draw, rounded corners, align=center, minimum width=15mm, minimum height=7mm, inner sep=2pt }, std/.style={sys, fill=blue!8}, fo/.style={sys, fill=red!10}, rel/.style={->, thick}, weak/.style={->, semithick, dashed}, equv/.style={<->, thick, dotted} ] \node[std] (bssn) {BSSN}; \node[std, right=of bssn] (z4) {Z4}; \node[std, right=of z4] (z4c) {Z4c}; \node[std, right=of z4c] (ccz4) {CCZ4}; \node[fo, below=8mm of bssn] (fobssn) {FOBSSN}; \node[fo, above=8mm of z4] (gh) {GH}; \node[fo, below=8mm of z4] (foz4) {FOZ4}; \node[fo, below=8mm of ccz4] (foccz4) {FOCCZ4}; \node[below=8mm of z4c] (foz4c) {?}; \draw[rel] (z4) -- (z4c); \draw[rel] (z4) to[bend left=18] (ccz4); \draw[weak] (bssn) -- (z4c); \draw[rel] (bssn) -- (fobssn); \draw[rel] (ccz4) -- (foccz4); \draw[rel] (z4) -- (foz4); \draw[equv] (z4) -- (gh); % \node[draw=none, above=1mm of gh, font=\scriptsize] {full first order}; \end{tikzpicture} \caption{Commonly used systems, and first order reductions.} \end{figure} \end{frame} \begin{frame}{Why not existing first order formulations?} \begin{itemize} \item GH \begin{itemize} \item Not compatible with moving puncture method \end{itemize} \item FOCCZ4 \begin{itemize} \item Not damping the new constraints introduced during the reduction \end{itemize} \item FOZ4 \begin{itemize} \item It's almost a copy of FOCCZ4. Not damping the new constraints introduced during the reduction \end{itemize} \item FOBSSN \begin{itemize} \item Damps only 2 out of 3 new constraints. Probably works; not well tested in the community; WIP in \texttt{nmesh} \end{itemize} \end{itemize} Since Z4c works well, we expect that a first order reduction of Z4c with proper constraint damping will also work well, and thus we want to find such a formulation. \end{frame} \begin{frame}{Hyperbolicity of first order systems} \subsection{Hyperbolicity of first order systems} Say we have a list of variables $u^I$, and a first order PDE system \begin{equation} \pdvt{u^I} + \tensor{A}{^i^I_J} \pdv{u^J}{x^i} = S^I, \label{first-order-PDE} \end{equation} where both $A$ and $S$ can depend on $u^I$ but not their derivatives. For any unit covector $\xi_i$, the \emph{Principal symbol} matrix is defined as $\tensor{P}{^I_J}(\xi) := \tensor{A}{^i^I_J} \xi_i$. The system is said to be: \begin{itemize} \item \emph{weakly hyperbolic} if for any $\xi_i$, $\tensor{P}{^I_J}(\xi)$ has only real eigenvalues. \item \emph{strongly hyperbolic} if for any $\xi_i$, $\tensor{P}{^I_J}(\xi)$ has only real eigenvalues and a complete set of eigenvectors, i.e., it is diagonalizable. \item \emph{symmetrically hyperbolic} if there exists a positive definite symmetrizer matrix $H_{IJ}$ such that $H_{IK} \tensor{P}{^K_J}(\xi)$ is always symmetric. \end{itemize} \begin{theorem} System~\eqref{first-order-PDE} is well-posed in the $L^2$ sense if and only if it is strongly hyperbolic. \end{theorem} \end{frame} \begin{frame}{First order reduction: wave equation as an example} \subsection{First order reduction: wave equation as an example} We can introduce auxiliary variables to reduce a second order PDE system to a first order one. For example, for the wave equation on flat spacetime: \begin{equation} \pt^2 {\phi} - \tensor{\delta}{^i^j} \Partial{_i} \Partial{_j} \phi = 0, \end{equation} we can introduce $\pi := \pt \phi$ and $\tensor{\psi}{_i} := \Partial{_i} \phi$, and rewrite the wave equation as \begin{equation} \begin{cases} \pt{\phi} = \pi, \\ \pt{\pi} = \tensor{\delta}{^i^j} \Partial{_i} \tensor{\psi}{_j},\\ \pt{\tensor{\psi}{_i}} = \Partial{_i} \pi. \end{cases} \end{equation} write $u^I = (\phi, \pi, \tensor{\psi}{_1}, \tensor{\psi}{_2}, \tensor{\psi}{_3})^T$, then the system can be written in the form of~\eqref{first-order-PDE} with \begin{equation} \tensor{A}{^i^I_J} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & -\tensor{e}{_i} \\ 0 & -\tensor{e}{^i} & 0 \end{pmatrix} \qc S=0, \end{equation} where $\tensor{e}{^i}$ is the standard column vector basis of $\setR^3$, while $\tensor{e}{_i}$ is the standard row vector basis of $(\setR^3)^*$. \end{frame} \begin{frame}{First order reduction: wave equation as an example} The principal symbol matrix is \begin{equation} \tensor{P}{^I_J}(\xi) = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & -\myvec{\xi}^T \\ 0 & -\myvec{\xi} & 0 \end{pmatrix}, \end{equation} it's very clear that the principal symbol matrix is already symmetric, and thus the system is symmetrically hyperbolic, and well-posed in the $L^2$ sense. The eigenvalues are $\lambda_1 = \lambda_2 = \lambda_3 = 0$, $\lambda_4 = -1$ and $\lambda_5 = 1$. The eigenvectors are \begin{equation} \begin{gathered} \myvec{e}_1 = \mqty(1,0,0,0,0)^T \qc \myvec{e}_2 = \mqty(0,0,-\xi_3,0,\xi_1)^T \qc \myvec{e}_3 = \mqty(0,0,-\xi_2,\xi_1,0)^T,\\ \myvec{e}_4 = \mqty(0,1,\xi_1,\xi_2,\xi_3)^T \qc \myvec{e}_5 = \mqty(0,-1,\xi_1,\xi_2,\xi_3)^T. \end{gathered} \end{equation} \end{frame} \begin{frame}{First order reduction: wave equation as an example}{constraints during the reduction} During the first order reduction, we introduce new variables. A solution to the new system is a solution to the original system if and only if the new variables satisfy some constraints. For example, in the wave equation case, we have the constraints \begin{equation} \mathcal{C}_i \definedby \tensor{\psi}{_i} - \Partial{_i} \phi = 0. \end{equation} The evolution of the constraints is given by \begin{equation} \pt \mathcal{C}_i = \pt \psi_i - \Partial{_i} \pt \phi = \Partial{_i} \pi - \Partial{_i} \pi = 0, \end{equation} which means that if the constraints are satisfied initially, they will be satisfied for all time. However, in numerical simulations, there will always be some constraint violation due to numerical errors, and thus it's better to add some constraint damping terms to the evolution equations to suppress the growth of constraint violation. \end{frame} \begin{frame}{First order reduction: wave equation as an example}{constraint damping} For example, we can add a constraint damping term $-\gamma \mathcal{C}_i$ to the evolution equation of $\tensor{\psi}{_i}$, \begin{equation} \begin{cases} \pt{\phi} = \pi, \\ \pt{\pi} = \tensor{\delta}{^i^j} \Partial{_i} \tensor{\psi}{_j},\\ \pt{\tensor{\psi}{_i}} = \Partial{_i} \pi - \gamma \left( \tensor{\psi}{_i} - \Partial{_i} \phi \right), \end{cases} \end{equation} where $\gamma > 0$ is a constant. Then the evolution of the constraints becomes \begin{equation} \pt \mathcal{C}_i = -\gamma \mathcal{C}_i, \end{equation} which means that the constraint violation will decay exponentially with time, and thus the system is more stable for numerical simulations. \end{frame} \begin{frame}{First order reduction: wave equation as an example}{constraint damping and hyperbolicity recheck} We have to check if the constraint damping term will change the hyperbolicity of the system. The new principal symbol matrix is \begin{equation} \tensor{P}{^I_J}(\xi) = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & -\myvec{\xi}^T \\ -\gamma \myvec{\xi} & -\myvec{\xi} & 0 \end{pmatrix}, \end{equation} which is no longer symmetric, but luckily it's still diagonalizable with exactly the same eigenvalues, and the eigenvectors are also almost the same, except that \begin{equation} \myvec{e}_1 = \mqty(1,-\gamma,0,0,0)^T, \end{equation} and thus the system is still strongly hyperbolic, and well-posed in the $L^2$ sense. \end{frame} \begin{frame}{Searching for a first order Z4c}{old workflow: by hand} \subsection{old workflow: by hand} Here is the Z4c evolution system with the puncture gauge condition: \[ \scalebox{0.6}{$\displaystyle \begin{aligned} \pt \chi ={}& \frac{2}{3} \chi \left[ \alpha \left( \hat{K} + 2\Theta \right) - \spaceD{_i} \shift{^i} \right],\\ \pt \smetric*{_i_j} ={} & - 2 \lapse \Ktf*{_i_j} + \shift{^k} \Partial{_k} \smetric*{_i_j} + 2 \smetric*{_{k(i}} \Partial{_{j)}} \shift{^k} - \frac{2}{3} \smetric*{_i_j} \Partial{_k} \shift{^k},\\ \pt \hat{K} ={}& - {\color{red}\spaceD{^i} \spaceD{_i} \lapse} + \lapse \left[ \Ktf*{_i_j} \Ktf*{^i^j} + \frac{1}{3} \left( \hat{K} + 2 \Theta \right)^2 \right] + 4 \pi \lapse \left( S + \rho \right) + \lapse \kappa_1 (1-\kappa_2) \Theta + \shift{^i} \Partial{_i} \hat{K},\\ \pt \Ktf*{_i_j} ={} & \chi \left[ - {\color{red}\spaceD{_i} \spaceD{_j} \lapse} + \lapse \left( {\color{red}\Rictensor{_i_j}} - 8 \pi \tensor{S}{_i_j} \right) \right]^{\text{tf}} + \lapse \left[ \left( \hat{K} + 2 \Theta \right) \Ktf*{_i_j} - 2 \Ktf*{^k_i} \Ktf*{_k_j} \right] + \shift{^k} \Partial{_k} \Ktf*{_i_j} + 2 \Ktf*{_{k(i}} \Partial{_{j)}} \shift{^k} - \frac{2}{3} \Ktf*{_i_j} \Partial{_k} \shift{^k},\\ \pt \Gam*{^i} ={} & - 2 \Ktf*{^i^j} \Partial{_j} \lapse + 2 \lapse \left[ \Gam*{^i_j_k} \Ktf*{^j^k} - \frac{3}{2} \Ktf*{^i^j} \Partial{_j} \ln(\chi) - \frac{1}{3} \smetric*{^i^j} \Partial{_j} \left( 2 \hat{K} + \Theta \right) - 8 \pi \smetric*{^i^j} \tensor{S}{_j} \right] + {\color{red}\smetric*{^j^k} \Partial{_j} \Partial{_k} \shift{^i}} + {\color{red}\frac{1}{3} \smetric*{^i^j} \Partial{_j} \Partial{_k} \shift{^k}} \notag\\ & {} + \shift{^j} \Partial{_j} \Gam*{^i} - \Gamd*{^j} \Partial{_j} \shift{^i} + \frac{2}{3} \Gamd*{^i} \Partial{_j} \shift{^j} - 2 \lapse \kappa_1 \left[ \Gam*{^i} - \Gamd*{^i} \right], \\ \pt \Theta ={} & \frac{1}{2} \lapse \left[ \Rscalar - \Ktf*{_i_j} \Ktf*{^i^j} + \frac{2}{3} \left( \hat{K} + 2 \Theta \right)^2 \right] - \lapse \left[ 8 \pi \rho + \kappa_1 (2+\kappa_2) \Theta \right] + \shift{^i} \Partial{_i} \Theta,\\ \pt \lapse ={} & - \mu_L \lapse^2 \hat{K} + \Ld{\beta} \lapse,\\ \pt \shift{^i} ={} & \mu_S \lapse^2 \Gam*{^i} - \eta \shift{^i} + \shift{^j} \Partial{_j} \shift{^i}. \end{aligned} $} \] where all the second order spatial derivatives are highlighted in red. \end{frame} \begin{frame}{Searching for a first order Z4c}{old workflow: by hand} We introduce the following auxiliary variables: \[ \scalebox{0.8}{$\displaystyle \dlapse{_i} \definedby \Partial{_i} \lapse \qc \dshift{^i_j} \definedby \Partial{_j} \shift{^i} \qc \dchi*{_i} \definedby \spaceD*{_i} \ln(\chi) \qc \dsmetric*{_i_j_k} \definedby \Partial{_i} \smetric*{_j_k}, $} \] and replace all the spatial derivatives of the original variables with the new auxiliary variables, we got this \[ \scalebox{0.47}{$\displaystyle \begin{aligned} \pt \chi ={}&\frac{2}{3} \chi \left[ \alpha \left( \hat{K} + 2\Theta \right) - \dshift{^i_i} - \Gam*{^j_j_k} \shift{^k} \right] + \shift{^k} \Partial{_k} \chi,\\ \pt \smetric*{_i_j} ={} & - 2 \lapse \Ktf*{_i_j} + \shift{^k} \Partial{_k} \smetric*{_i_j} + 2 \smetric*{_{k(i}} \dshift{^k_{j)}} - \frac{2}{3} \smetric*{_i_j} \dshift{^k_k},\\ \pt \hat{K} ={}& - \spaceD{^i} \dlapse{_i} + \lapse \left[ \Ktf*{_i_j} \Ktf*{^i^j} + \frac{1}{3} \left( \hat{K} + 2 \Theta \right)^2 \right] + 4 \pi \lapse \left( S + \rho \right) + \lapse \kappa_1 (1-\kappa_2) \Theta + \shift{^i} \Partial{_i} \hat{K},\\ \pt \Ktf*{_i_j} ={} & \chi \left[ - \spaceD{_{(i}} \dlapse{_{j)}} + \lapse \left( \Rictensor{_i_j} - 8 \pi \tensor{S}{_i_j} \right) \right]^{\text{tf}} + \lapse \left[ \left( \hat{K} + 2 \Theta \right) \Ktf*{_i_j} - 2 \Ktf*{^k_i} \Ktf*{_k_j} \right] + \shift{^k} \Partial{_k} \Ktf*{_i_j} + 2 \Ktf*{_{k(i}} \dshift{^k_{j)}} - \frac{2}{3} \Ktf*{_i_j} \dshift{^k_k},\\ \pt \Gam*{^i} ={} & - 2 \Ktf*{^i^j} \dlapse{_j} + 2 \lapse \left[ \Gam*{^i_j_k} \Ktf*{^j^k} - \frac{3}{2} \Ktf*{^i^j} \dchi*{_j} - \frac{1}{3} \smetric*{^i^j} \Partial{_j} \left( 2 \hat{K} + \Theta \right) - 8 \pi \smetric*{^i^j} \tensor{S}{_j} \right] + \smetric*{^j^k} \Partial{_j} \dshift{^i_k} + \frac{1}{3} \smetric*{^i^j} \Partial{_j} \dshift{^k_k} + \shift{^j} \Partial{_j} \Gam*{^i} - \Gamd*{^j} \dshift{^i_j} + \frac{2}{3} \Gamd*{^i} \dshift{^j_j} - 2 \lapse \kappa_1 \left[ \Gam*{^i} - \Gamd*{^i} \right], \\ \pt \Theta ={} & \frac{1}{2} \lapse \left[ \Rscalar - \Ktf*{_i_j} \Ktf*{^i^j} + \frac{2}{3} \left( \hat{K} + 2 \Theta \right)^2 \right] - \lapse \left[ 8 \pi \rho + \kappa_1 (2+\kappa_2) \Theta \right] + \shift{^i} \Partial{_i} \Theta,\\ \pt \lapse ={}& - \mu_L \lapse^2 \hat{K} + \shift{^i} \Partial{_i} \lapse,\\ \pt \shift{^i} ={}& \mu_S \lapse^2 \Gam*{^i} - \eta \shift{^i} + \shift{^j} \Partial{_j} \shift{^i},\\ \pt \dlapse{_i} ={} & - \pdv{\mu_L}{\lapse} \lapse^2 \hat{K} \dlapse{_i} - 2 \mu_L \lapse \hat{K} \dlapse{_i} - \mu_L \lapse^2 \Partial{_i} \hat{K} + \dshift{^j_i} \dlapse{_j} + \shift{^j} \Partial{_j} \dlapse{_i},\\ \pt \dshift{^i_j} ={} & \pdv{\mu_S}{\lapse} \lapse^2 \Gam*{^i} \dlapse{_j} + 2 \mu_S \lapse \Gam*{^i} \dlapse{_j} + \mu_S \lapse^2 \Partial{_j} \Gam*{^i} - \pdv{\eta}{\lapse} \dlapse{_j} \shift{^i} - \eta \dshift{^i_j} + \dshift{^k_j} \dshift{^i_k} + \shift{^k} \Partial{_k} \dshift{^i_j},\\ \pt \dchi*{_i} ={} & \frac{2}{3} \left[ \dlapse{_i} \left( \hat{K} + 2\Theta \right) + \lapse \Partial{_i} \left( \hat{K} + 2\Theta \right) - \Partial{_i} \dshift{^j_j} - \frac{1}{2} \shift{^k} \dsmetric*{_k_l_j} \Partial{_i} \smetric*{^l^j} - \frac{1}{2} \smetric*{^l^j} \shift{^k} \Partial{_i} \dsmetric*{_k_l_j} - \left( \Gam*{^j_j_k} - \frac{3}{2} \dchi*{_k} \right) \dshift{^k_i} \right] + \shift{^k} \Partial{_k} \dchi*{_i},\\ \pt \dsmetric*{_{ijk}} ={} & -2 \dlapse{_i} \Ktf*{_j_k} -2 \lapse \Partial{_i} \Ktf*{_j_k} + \dshift{^l_i} \dsmetric*{_l_j_k} + \shift{^l} \Partial{_l} \smetric*{_i_j_k} + 2 \dsmetric*{_{il(j}} \dshift{^l_{k)}} + 2 \smetric*{_{l(j}} \Partial{_{|i|}} \dshift{^l_{k)}} - \frac{2}{3} \dsmetric*{_i_j_k} \dshift{^l_l} - \frac{2}{3} \smetric*{_j_k} \Partial{_i} \dshift{^l_l}, \end{aligned} $} \] \end{frame} \begin{frame}{Searching for a first order Z4c}{old workflow: by hand} We have several new constraints: \begin{gather*} \constraintA{_i} \definedby \dlapse{_i} - \Partial{_i} \lapse \qc \constraintB{^i_j} \definedby \dshift{^i_j} - \Partial{_j} \shift{^i}, \\ \constraintC*{_i} \definedby \dchi*{_i} - \spaceD*{_i} \ln(\chi) \qc \constraintD*{_i_j_k} \definedby \dsmetric*{_i_j_k} - \Partial{_i} \smetric*{_j_k} \end{gather*} constraint dynamics: \[ \scalebox{0.8}{$\displaystyle \begin{aligned} \pt \constraintA{_i} &= - \pdv{\mu_L}{\lapse} \lapse^2 \hat{K} \constraintA{_i} - 2 \mu_L \lapse \hat{K} \constraintA{_i} + \left( \Partial{_i} \shift{^j} \right) \constraintA{_j} + \dlapse{_j} \constraintB{^j_i} + \shift{^j} \Partial{_i} \constraintA{_j},\\ \pt \constraintB{^i_j} &= \pdv{\mu_S}{\lapse} \lapse^2 \Gam*{^i} \constraintA{_j} + 2 \mu_S \lapse \Gam*{^i} \constraintA{_j} - \pdv{\eta}{\lapse} \shift{^i} \constraintA{_j} - \eta \constraintB{^i_j} + \dshift{^k_j} \constraintB{^i_k} + \left( \Partial{_k} \shift{^i} \right) \constraintB{^k_j} + \shift{^k} \Partial{_k} \constraintB{^i_j},\\ \pt \constraintC*{_i} &= \frac{2}{3} \left[ \left( \hat{K} + 2\Theta \right) \constraintA{_i} - \left( \Gam*{^j_j_k} - \frac{3}{2} \dchi*{_k} \right) \constraintB{^k_i} \right] + \shift{^k} \Partial{_i} \constraintC*{_k} + \left( \Partial{_i} \shift{^k} \right) \constraintC*{_k},\\ \pt \constraintD*{_i_j_k} &= -2 \Ktf*{_j_k} \constraintA{_i} + \dsmetric*{_l_j_k} \constraintB{^l_i} + \shift{^l} \Partial{_l} \constraintD*{_i_j_k} + 2 \constraintD*{_{il(j}} \dshift{^l_{k)}} - \frac{2}{3} \dshift{^l_l} \constraintD*{_i_j_k}, \end{aligned} $} \] these are closed, but not damping. \end{frame} \begin{frame}{Searching for a first order Z4c}{old workflow: by hand} The most simple terms to damp the constraints: \[ \scalebox{0.8}{$\displaystyle \begin{aligned} \pt \dlapse{_i} &= \dots - c_1 \constraintA{_i},\\ \pt \shift{^i_j} &= \dots - c_2 \constraintB{^i_j},\\ \pt \dchi*{_i} &= \dots - c_3 \constraintC*{_i},\\ \pt \dsmetric*{_{ijk}} &= \dots - c_4 \constraintD*{_i_j_k}, \end{aligned} $} \] where $c_1$, $c_2$, $c_3$ and $c_4$ are positive constants that large enough (the lower bound can be calculated). Thus we have constraint damping. \end{frame} \begin{frame}{Searching for a first order Z4c}{old workflow: by hand} {\color{red} \textbf{Unfortunately, this naive system is not even weakly hyperbolic!}} \vspace{0.2cm} We have to add some constraint terms that contribute to the principal symbol, until the principal symbol matrix is diagonalizable with real eigenvalues. Available terms to add: \[ \scalebox{0.8}{$\displaystyle \begin{gathered} \constraintA{_i} \qc \constraintB{^i_j} \qc \constraintC*{_i} \qc \constraintD*{_i_j_k},\\ \Partial{_{[i}} \constraintA{_{j]}} = \Partial{_{[i}} \dlapse{_{j]}} \qc \Partial{_{[i}} \constraintB{^j_{k]}} = \Partial{_{[i}} \dshift{^j_{k]}} \qc \Partial{_{[i}} \constraintC*{_{j]}} = \Partial{_{[i}} \dchi*{_{j]}} \qc \Partial{_{[i}} \constraintD*{_{j]kl}} = \Partial{_{[i}} \dsmetric*{_{j]kl}},\\ \Partial{_i} \left( \smetric*{^j^k} \Ktf*{_j_k} \right) \qc \Partial{_i} \left( \smetric*{^k^l} \dsmetric*{_j_k_l} \right) \end{gathered} $} \] \end{frame} \begin{frame}{Searching for a first order Z4c}{Introducing Mathematica into workflow} \subsection{Introducing Mathematica into workflow} To speed up the process and reduce the chance of making mistakes, I wrote a Mathematica notebook using a package \texttt{MathGR}. \texttt{MathGR} can do symbolic tensor calculations, and render tensor equations in a readable way: \begin{figure} \centering \includegraphics[width=0.5\textwidth, trim=4bp 4bp 4bp 4bp, clip]{imgs/mathGR1.png} \caption{A simple example of using \texttt{MathGR} to write tensor equations.} \end{figure} \end{frame} \begin{frame}{Searching for a first order Z4c}{A new idea: start from FOCCZ4} Another way to speed up is, maybe we can start from FOCCZ4: \begin{itemize} \item FOCCZ4 is already strongly hyperbolic \item FOCCZ4 only differs from our naive first order reduction of Z4c by some constraint terms, to be specific, some terms proportional to $\tensor{\tilde{Z}}{^i}$ or it's derivatives. We can remove these differing terms and watch how the hyperbolicity changes, and find out how to recover the hyperbolicity \item We also add constraint damping terms, which will break the hyperbolicity, and then find out how to recover the hyperbolicity. \end{itemize} In short, we are looking for a series of hyperbolicity-preserving modifications to FOCCZ4, to get a first order reduction of Z4c with constraint damping. This is working in progress. \end{frame} \begin{frame}{Searching for a first order Z4c}{Current status} \section{Current status} \begin{itemize} \item What we have: \begin{itemize} \item A Mathematica notebook that takes readable tensor equations as input, and automatically computes the principal symbol matrix \end{itemize} \item What we are doing: \begin{itemize} \item In a loop of modifying the evolution equations, look at the principle symbol, checking the hyperbolicity, and modifying again. \end{itemize} \item What's hard: \begin{itemize} \item The principle symbol matrix is large, and it's hard to evaluate if it's diagonalizable or not. \item There are many available terms to add, and there is no general principle on how to choose the terms that are more likely to recover the hyperbolicity. \end{itemize} \item What we got so far: \begin{itemize} \item Found some weakly hyperbolic systems. \item A roadmap on FOCCZ4 $\rightarrow$ FOZ4c. \item A Mathematica based workflow to test hyperbolicity. \item Some experience based "feeling" on which terms should be added and where. \end{itemize} \end{itemize} \end{frame} \begin{frame}{End} \centering {\Huge Thanks!} \end{frame} \begin{frame}{References} \printbibliography \end{frame} \appendix \begin{frame}{Appendix} \section{Appendix} \begin{figure} \centering \includegraphics[width=0.7\textwidth, trim=4bp 4bp 4bp 4bp, clip]{imgs/moving_puncture_penrose.png} \caption{The numerical time slices in the moving puncture method in the Penrose diagram, from \cite{Hannam_2008}.} \end{figure} \end{frame} \begin{frame}{Appendix}{Why can't we fully automate the workflow using Mathematica?} There are coefficients before each constraint term we add, and it's very common that they have to meet some conditions to recover the hyperbolicity, for example, in GH they require that $\gamma_3 = \gamma_1 \gamma_2$, and in FOBSSN they require that $\kappa^\phi = 0$. And many coefficients cannot be parameters at all, they may have to be something specific like the lapse function $\lapse$. Then what about adding all possible terms with undetermined coefficients, and solve the coefficients by the conditions of hyperbolicity? Yes, I thought about it. We need a function that takes the principal symbol matrix with lots of undetermined coefficients as input, and outputs the conditions on the coefficients for the principal symbol matrix to be diagonalizable with real eigenvalues. I do wrote such a function, but it cannot be simple. When apply to our case, it will run forever. \end{frame} \begin{frame}{Appendix}{Why can't we fully automate the workflow using Mathematica?} More details: we cannot just compute the eigenvalues and eigenvectors, then check if the eigenvectors form a complete basis. Think about this simple matrix: \begin{equation*} \mqty(1 & a \\ 0 & 1), \end{equation*} if we ask Mathematica to compute the eigenvectors, it will give us $\myvec{v}_1 = \mqty(1 \\ 0)$ and $\myvec{v}_2 = \mqty(0 \\ 0)$, because unless $a=0$, this matrix have a Jordan block. So the geometric multiplicity is less than the algebraic multiplicity except on a zero-measure set, that's why Mathematica thinks there is no second linearly independent eigenvector. If we simply check if the eigenvectors are linearly independent, we will miss the condition $a=0$ and get false. \end{frame} \begin{frame}{Appendix}{Why can't we fully automate the workflow using Mathematica?} The correct way is to check the geometric multiplicity for each eigenvalues. Say we get eigenvalues $\lambda_1, ..., \lambda_k$ with the algebraic multiplicity $m_1, ..., m_k$. Under the assumption that all these eigenvalues are different from each other, we check $\dim\ker(A - \lambda_i I)$ for each $i$. $\dim\ker(A - \lambda_i I) \leq m_i$, so we just have to find the condition for $\dim\ker(A - \lambda_i I) \geq m_i$. $\dim\ker(A - \lambda_i I) = n - \rank(A - \lambda_i I)$, so we check $\rank(A - \lambda_i I) \leq n - m_i$. The criterion is that any $(n-m_i+1)\times(n-m_i+1)$ minor of $A - \lambda_i I$ is zero. Nice, this is a computable condition! But physically we expect most modes to have the same characteristic speed $\shift{^\xi}$, resulting a algebraic multiplicity $\sim 22$, while the size of matrix is 55, then the number of minors to check for $\lambda = \shift{^\xi}$ is \begin{equation*} \binom{55}{34}^2 = 708507400701023142362023505625 = 7.085\times 10^{29}. \end{equation*} It's impossible to compute so many minors. \end{frame} \end{document}