update: auto commit

Makefile

@@ -1,10 +1,11 @@
 main = report
-latexmk = latexmk -pdflatex
-git_add_files = *.tex imgs/ ref.bib
+latexmk = latexmk -pdf
+git_add_files = *.tex imgs/ ref.bib Makefile latexmkrc
 
 .PHONY : main clean cleanall autopush
 
 main: $(main).tex
+	mkdir -pv tikzcache
 	$(latexmk) $<
 
 clean:
@@ -15,4 +16,4 @@ cleanall:
 
 autopush: main
 	git add $(git_add_files)
-	git diff --cached --quiet || (git commit -m "update: auto commit" && git push)
+	git diff --cached --quiet || (git commit -m "update: auto commit" && git push)

latexmkrc

@@ -0,0 +1,2 @@
+$pdflatex = "pdflatex -file-line-error -interaction=nonstopmode -shell-escape -synctex=1 %O %S";
+$clean_ext = "bbl synctex.gz thm auxlock nav snm";

mystyle_beamer.tex

@@ -27,7 +27,7 @@
 %\usepackage{pifont}
 \usepackage{bm}
 % \usepackage{nomencl}
-\usepackage[backref=true]{biblatex}
+\usepackage[backref=true,style=authoryear]{biblatex}
 
 \hypersetup{
 	colorlinks=true,
@@ -165,7 +165,7 @@
 \newcommand{\AD}[1]{\tensor{\mathcal{D}}{_{#1}}}
 \newcommand{\ADd}{\mathcal{D}}
 \newcommand{\connection}[1]{\mathcal{#1}}
-\newcommand{\Partial}[1]{\tensor{\partial}{_{#1}}}%普通导数算符
+\newcommand{\Partial}[1]{\tensor{\partial}{#1}}%普通导数算符
 \newcommand{\tPartial}[1]{\tensor{\tilde{\partial}}{_{#1}}}
 \newcommand{\Fd}[2][\tau]{\ensuremath{\frac{\mathrm{D_F}#2}{\dd{#1}}}}%费米导数
 \newcommand{\Fdd}[2][\tau]{\ensuremath{\mathrm{D_F}#2/\dd{#1}}}%行内费米导数
@@ -222,6 +222,7 @@
 
 \newcommand{\form}[1]{\ensuremath{\bm{#1}}}
 \newcommand{\Ric}{\ensuremath{R}}
+\newcommand{\Rscalar}{\ensuremath{R}}
 \newcommand{\curR}{\ensuremath{{R}}}
 \newcommand{\spacecurR}{\ensuremath{\mathcal{R}}}
 \newcommand{\vol}{\ensuremath{\varepsilon}}

report.tex

@@ -5,14 +5,195 @@
 \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}}
+}
 
-\title{Towards Moving Picture Simulations in a Discontinuous Galerkin Framework}
+\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}
 
@@ -49,19 +230,25 @@
     \begin{itemize}
       \item LISA is able to detect extreme mass ratio inspirals (EMRIs)
       \item Compare with perturbation theory
+      \begin{equation*}
+        \mu = \frac{m_1 m_2}{(m_1+m_2)^2} \rightarrow 0,
+      \end{equation*}
+      expanding in $\mu$ $\implies$ gravitational self-force 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}
+    
+  %   There are two main methods to evolve black holes in numerical relativity:
+  %   \begin{itemize}
+  %     \item 
+  %     \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}
+    \subsection{Methods on evolving black holes}
+    Excision method: excise the black hole interior from the computational domain, and impose boundary conditions on the excision surface.
     \begin{figure}
       \centering
       \begin{tikzpicture}[>=Latex, line cap=round, line join=round]
@@ -98,11 +285,12 @@
             -- (2*\L,0.5*\L) node[right,align=left]
             {time slice};
         \end{tikzpicture}
-        \caption{excision on a single Schwarchild black hole}
+        \caption{excision on a single Schwarzschild black hole}
       \end{figure}
   \end{frame}
 
   \begin{frame}{Methods on evolving black holes: excision}
+    Excision method: excise the black hole interior from the computational domain, and impose boundary conditions on the excision surface
     \begin{figure}
       \centering
       \includegraphics[width=0.7\textwidth]{imgs/black_hole_excision_mesh.png}
@@ -112,25 +300,27 @@
 
 
   \begin{frame}{Methods on evolving black holes: moving puncture}
-    In the isotropic coordinate for Schwarzschild black hole, the spical 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, with the conformal factor
-    \begin{equation*}
-      \psi = 1 + \frac{M}{2r}.
-    \end{equation*}
+    Moving puncture method: use a suitable gauge condition to make sure the singularity is not on our time slices at all.
 
-    General $N$ black hole puncture initial data:
-    \begin{equation*}
-      \psi = \sum_{i=1}^N \frac{M_i}{2|\vb*{r}-\vb*{r}_i|} + (\text{regular part}),
-    \end{equation*}
-    again, if we factor out the conformal factor $\psi$, the remaining conformal metric is regular everywhere.
+    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 mothod, we won't cut out the black hole singularity, but choose a gauge condition that makes the singularity invisible to the numerical evolution.
+    % 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}
@@ -149,7 +339,7 @@
 
   \begin{frame}{Discontinuous Galerkin method and \texttt{nmesh}}
     \subsection{Discontinuous Galerkin method and \texttt{nmesh}}
-    There are two main numerical methods in numerical relativity:
+    There are two main families of numerical methods in numerical relativity:
     \begin{itemize}
       \item Finite Difference / Finite Volume
       \begin{itemize}
@@ -164,14 +354,14 @@
 
   \begin{frame}{Discontinuous Galerkin method and \texttt{nmesh}}
     Core idea of DG method: if we have a first-order PDE system:
-    \begin{equation*}
+    \begin{equation}
       \pdvt{u} + A^i \pdv{u}{x^i} = S,
-    \end{equation*}
+    \end{equation}
     then for some test functions $\{v_a\}$, we have
-    \begin{equation*}
+    \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 intergrating by parts. Boundary terms are replaced by numerical fluxes.
+    \end{equation}
+    which reduce to an ODE 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}
@@ -201,19 +391,13 @@
     \texttt{nmesh} has two features that are useful for our project:
     \begin{itemize}
       \item Adaptive mesh refinement (AMR)
-      \begin{figure}
+      \item DG + FD/FV dynamically switching
+    \end{itemize}
+    \begin{figure}
         \centering
-        \includegraphics[width=0.7\textwidth]{imgs/MPA1_W-9sn12l5_GRHD_D_t0400.pdf}
-        \caption{The mesh in a binery star emulation with \texttt{nmesh}}
+        \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{itemize}
-  \end{frame}
-
-  \begin{frame}{Discontinuous Galerkin method and \texttt{nmesh}}
-    \texttt{nmesh} has two features that are useful for our project:
-    \begin{itemize}
-      \item DG/FV dynamically switching
-    \end{itemize}
   \end{frame}
 
   \begin{frame}{Moving puncture evolution in \texttt{nmesh}?}
@@ -238,27 +422,487 @@
     $\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}{Existing first order formulations}
-    \section{Towards first order Z4c}
-    \subsection{Existing first order formulations}
+  \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 mothod
+        \item Not compatible with moving puncture method
       \end{itemize}
       \item FOCCZ4
       \begin{itemize}
-        \item No constraint damping
+        \item Not damping the new constraints introduced during the reduction
       \end{itemize}
       \item FOZ4
       \begin{itemize}
-        \item No constraint damping
+        \item It's almost a copy of FOCCZ4. Not damping the new constraints introduced during the reduction
       \end{itemize}
       \item FOBSSN
       \begin{itemize}
-        \item Probably works; not well tested in the community; WIP in \texttt{nmesh}
+        \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}
@@ -276,4 +920,30 @@
     \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}