\subsection{The Galerkin Difference Basis}\label{GD_basis}
We now describe the Galerkin difference (GD) basis we will use. Note that the description here differs slightly different from the original description in \cite{Gdiff}, but the basis is identical.

We first consider grid points well-separated from a boundary. The goal is to construct a polynomial basis of odd degree $p$ on an equidistant grid with grid spacing $h$. The generating basis function $\Phi_p(x)$ is centered around $x=0$ and the basis itself is simply the union of the translates of the generating basis function. Thus an element in the basis centered around $x_i = ih$ becomes, $\phi_{p,i} = \Phi_p(x-ih)$.
\graphicspath{{graphs_GDbasis/}}
\begin{figure}[h]
\begin{center}
\includegraphics[width=0.27\textwidth]{GD_generating}\ \ \ \
\includegraphics[width=0.21\textwidth]{GD_all_p7}
\includegraphics[width=0.41\textwidth]{GD_all_3D}
\caption{To the left, the generating basis function $\Phi_p(x)$ for $p = 1,3,5,7,9$. Note that $\Phi_p(x) = \Phi_p(-x)$. Note that each function is vertically offset by $q-1$. In the middle, the ``right'' four degree 7 Lagrange polynomials that make up $\Phi_7(x)$, the part of $L_j$ that is used is in bold. To the right, a three dimensional representation of the generating basis function $\Phi_p(x)$ for $p$ up to 39. \label{fig:GD_generating}}
\end{center}
\end{figure}

The generating basis function $\Phi_p(x)$ is symmetric, $\Phi_p(x) = \Phi_p(-x)$, and has compact support on \mbox{$x \in [-qh,qh]$}, where $p = 2q-1$ (recall that $p$ is odd). In Figure \ref{fig:GD_generating} we display the non-zero part of $\Phi_p(x)$  for $p = 1,3,5,7,9,$ and $x>0$. In the lower left corner, where $p=1$, we recognize the classic finite element hat function and as $p$ increases we see that $\Phi_p(x)$ becomes increasingly similar to the {\em Cardinal Sinc} function.

An explicit formula for $\Phi_p(x)$ inside each of the $q$ positive intervals $[jh,(j+1)h)$ is obtained as follows. Let $L_j(x)$ is be the Lagrange interpolating polynomial on the grid $G_j = \{-jh,\ldots,(p-j)h\}$ with the property that $L_j(0) = 1$, then
\begin{eqnarray}\label{GD_basis_lagrange}
\Phi_p(x) = L_j(x), \ \ x \in [(q-j-1)h,(q-j)h)	.
\end{eqnarray}

A continuous function $u(x,t)$ can then be approximated by a linear combination of basis functions with weights corresponding to nodal values
\begin{eqnarray}\label{GD_basis_nodalformula}
u(x,t) \approx \sum_{i=k-(q-1)}^{k+q} u_i \phi_{p,i}(x), \ \  x \in [kh,(k+1)h),
\end{eqnarray}
where $u_i = u(x_i,t)$.

\subsubsection{Modification Near Boundaries}\label{GD_basis_subsection_1}
Near boundaries the basis must be modified. In \cite{Gdiff}, three approaches for handling boundaries are described: ghost basis, extrapolation basis and use of modified equations. Here we will exclusively use the extrapolation basis, which we describe next.

The extrapolation procedure draws from the standard practice to use ghost points in finite difference methods. First, note that the $q-1$ additional ghost basis functions associated with the $q-1$ first grid points outside the computational domain are the only ghost basis with support inside the computational domain. In the ghost basis approach the degrees of freedom at the ghost points are retained as unknowns but in the extrapolation approach they are eliminated in favor of modifying the basis itself near the boundary.

As the name suggests, the elimination is done by extrapolating the nodal values inside the computational domain to the ghost points. For example consider $p=3$. Then $q-1 = 1$ and one ghost point value, $u_{-1}$, must be determined. As the basis is fourth order accurate, the ghost point value is determined by fourth order accurate extrapolation,
\begin{eqnarray*}
u_{-1} &=& 4u_0 - 6u_1 + 4u_2 -u_3.
\end{eqnarray*}
 To understand how the modified basis is constructed, consider evaluating the approximation $u(x)$ inside the computational domain where the ghost basis has support. In this case this means \mbox{$x\in(x_0,x_1)$} and the approximation is
\begin{eqnarray*}
u(x) &=& u_{-1}\phi_{-1}(x) + u_{0}\phi_{0}(x) + u_{1}\phi_1(x) + u_2\phi_2(x).
\end{eqnarray*}
To obtain a value for $u_{-1}$ we use the extrapolation condition
\begin{eqnarray*}
u(x) &=& ( 4u_0 - 6u_1 + 4u_2 -u_3)\phi_{-1}(x) + u_{0}\phi_{0}(x) + u_{1}\phi_1(x) + u_2\phi_2(x) \nonumber\\
&=& u_0 \left[ \phi_0(x)+4\phi_{-1}(x)\right] +u_1\left[ \phi_1(x)-6\phi_{-1}(x) \right] + u_2 \left[ \phi_2(x)+4\phi_{-1}(x) \right] + u_3 \left[-\phi_{-1}(x) \right] \\
& =& u_0 \left[ \phi_0(x)+4\phi_{-1}(x) \right] +u_1 \left[ \phi_1(x)-6\phi_{-1}(x) \right] + u_2 \left[ \phi_2(x)+4\phi_{-1}(x) \right] + u_3 \left[\phi_3(x)-\phi_{-1}(x ) \right] .\nonumber
\end{eqnarray*}
In the last step, we used the fact the support of $\phi_3$ vanishes in $(x_0,x_1)$. Thus the modified basis functions are
\begin{equation*}
\tilde{\phi}_0 = \phi_0 + 4\phi_{-1}, \ \ \ \
\tilde{\phi}_1 = \phi_1 - 6\phi_{-1},   \ \ \ \
\tilde{\phi}_2 = \phi_2 + 4\phi_{-1},  \ \ \ \
\tilde{\phi}_3 = \phi_3 - \phi_{-1}.
\end{equation*}
The extension to larger $p$ requires the basis to be modified in a wider band near the boundaries and the extrapolation is done at the order of accuracy that matches that of the interior scheme. The handling of the right boundary is analogous.

\subsubsection{Extension to Higher Dimensions}\label{GD_basis_subsection_2}
The extension to higher dimensions simply amounts to using the tensor product approximation built off the one dimensional basis. For example, in two dimensions we have
\begin{eqnarray}
u(x,y,t) \approx \sum_{i=k_x-(q-1)}^{k_x+q}\sum_{j=k_y-(q-1)}^{k_y+q}u_{i,j} \phi_{p,i,j}(x,y), \ \ (x,y)\in[k_xh,(k_x+1)h)\times[k_yh,(k_y+1)h),
\end{eqnarray}
where $u_{i,j} = u(x_i,y_j,t)$ and
\begin{equation}
\phi_{p,i,j}(x,y) = \left\{
\begin{array}{lll}
 \phi_{p,i}(x)\phi_{p,j}(y), & p <   i < N-p, & p < j < N - p, \\
 {\phi}_{p,i}(x)\tilde{\phi}_{p,j}(y), & p<i<N-p, & 0 \le j \le N,\\
 \tilde{\phi}_{p,i}(x)\phi_{p,j}(y), &  0 \le i \le N, & p< j < N-p,\\
\tilde{\phi}_{p,i}(x) \tilde{\phi}_{p,j}(y), & i,j\leq p, &  i,j \geq N-p.
\end{array}
\right.
\end{equation}

Below, for notational convenience, we will not explicitly distinguish between the modified basis functions and the interior basis functions and simply drop the tilde notation. Also, we will use the notation $\mathbb{Q}_{p,N}$ to denote the one dimensional space spanned by the $(N+1)$ Galerkin difference basis functions associated with the nodal degrees of freedom.

\subsubsection{Alternative Galerkin Difference Spaces}\label{DGDspace}

Lastly we note that instead of the locally (within element) continuous Galerkin basis, one can use discontinuous even-degree polynomials
constructed using cell-centered interpolation nodes as in \cite{GDDG,sipDGGD}. In this case we would need to introduce fluxes not only at the
boundaries of the macro-elements as we do here, but also at the boundaries of the cells within each macro-element. The construction of the
energy-based discontinuous Galerkin method would follow in the same way as presented below, and we could also construct the improved basis for
elements in regions where the wave speed is constant. We note that in this case it would be possible to make different choices for the fluxes for
the cell boundaries interior to each element and those between the macro-elements. However, here we will focus solely on the use of the
locally continuous basis functions.

\section{Efficient Formulation on Cartesian Grids}\label{GD_cost}
In this section we restrict our attention to the case of constant speed of sound within an un-mapped Cartesian element. Using the
tensor products of the compactly-supported basis functions described in Section~\ref{GD_basis} we will see that the complexity of
computing the time derivatives is unacceptable for $N$ large due to the structure of the lift matrix associated with (\ref{DG1}).
However, if we use
a simple simultaneous diagonalization of the mass and stiffness matrix we can construct practical implementations of (\ref{DG1}), (\ref{DG2}) and (\ref{DG1_extra}) whose cost scales linearly with the total number of degrees of freedom.

On each element $\Omega^k$, we approximate the solution by tensor product expansions (here for the case of three dimensions)
\begin{eqnarray*}
u(x,y,z,t) &=& \sum_{l_x=0}^{N} \sum_{l_y=0}^{N} \sum_{l_z=0}^{N} u_{l_x,l_y,l_z} \phi_{p,l_x}(x) \phi_{p,l_y}(y) \phi_{p,l_z}(z),\\
v(x,y,z,t) &=& \sum_{l_x=0}^{N} \sum_{l_y=0}^{N} \sum_{l_z=0}^{N} v_{l_x,l_y,l_z} \phi_{p,l_x}(x) \phi_{p,l_y}(y)\phi_{p,l_z}(z).
\end{eqnarray*}
On element $\Omega^k$, let $U^k$, $V^k$ be column vectors containing the nodal values of $u$ and $v$, that is $u_{l_x,l_y,l_z}$ and $v_{l_x,l_y,l_z}$, respectively. Then, in $d$-dimensions we may write the nodal based version of the method as the system of ordinary differential equations
\begin{eqnarray}
\hat{S} \frac{d U^k}{dt} &=& \hat{S} V^k  \nonumber \\
&& +
\sum_{j=1}^d
(\alpha-1)\left[(D^{R,R}-D^{L,L})V^k
- D^{R,L}V^{k+1}
+ D^{L,R}V^{k-1}\right] \nonumber\\
&& - \sum_{j=1}^d  \beta \left[ (C^{R,R}-C^{L,L}) U^k
+ C^{R,L} U^{k+1}
 - C^{L,R}U^{k-1} \right], \label{DG_matrix1}
\end{eqnarray}
\begin{eqnarray}
M \frac{d V^k}{d t }  &=& - c^2 S U^k  \nonumber\\
&& + c^2\sum_{j=1}^{d}
(1-\alpha)\left( E^{R,R} -E^{L,L} \right)U^k
+ \alpha \left(E^{R,L}U^{k+1}
- E^{L,R}U^{k-1} \right) \nonumber\\
&&- c^2\sum_{j=1}^{d} \tau\left[ \left( B^{R,R} - B^{L,L} \right) V^k
+ B^{R,L}V^{k+1}
- B^{L,R}V^{k-1} \right].\label{DG_matrix2}
\end{eqnarray}
Here we abuse the notation in that for each coordinate direction in the sums we use the superscript $ k \pm 1$ denote the element ``left'' and ``right'' of element $k$ in the $j$th direction. The definitions of the mass matrix $M$, the stiffness matrix $\hat{S}$ and the lift matrices $B,C,D,E$ will be given below.

\subsection{Complexity with Galerkin Difference Basis}
Now, in order to compute the time derivatives $\frac{d U^k}{d t }$ and $\frac{d V^k}{d t}$ we must evaluate the matrix vector products on the right hand side and the action of the lift matrices on $U^k$ and $V^k$. As the matrices are sparse it is possible to do this at a cost that scales as $\sim f(p) N^\kappa$, with $f(p)$ being a low degree polynomial in $p$ and $\kappa = d$ for the volume terms and $\kappa = d-1$ for the surface terms. In other words the cost scales linearly with the number of degrees of freedom.

Further, due to the tensor product structure of the mass matrix we have that the element mass matrix $M$ can be composed as a Kronecker product of the one dimensional matrix, which we denote $M_j$,
\[
M = \otimes_{j=1}^{d}  M_j,
\]
with $M_{j,kl} = \int_{L_j}^{R_j} \phi_{j,k}\phi_{j,l}d_{x_j}$. Now, as the one dimensional mass matrices have bandwidth $p$ so will its $LU$-factors. Let $ L_j U_j \equiv M_j$. Then by the Hadamard product property we have
\begin{equation*}
M = \otimes_{j=1}^d M_j = \otimes_{j=1}^d L_j U_j = (\otimes_{j=1}^d L_j)(\otimes_{j=1}^d U_j) \equiv LU.
\end{equation*}
Thus, as the cost of each substitution is $\mathcal{O}(p(N+1))$ the cost of solving $Mx = b$ is $\mathcal{O}((p(N+1))^d)$, which again is linear in the degrees of freedom. Unfortunately,  the stiffness matrix $\hat{S}$ is a sum of Kronecker products
\begin{equation*}
\hat{S} = S_1\otimes M_2 \otimes M_3 + M_1\otimes S_2\otimes M_3 + M_1\otimes M_2\otimes S_3,
\end{equation*}
with $S_{j,kl} = \int_{L_j}^{R_j}\frac{d\phi_{j,k}}{dx_j}\frac{d\phi_{j,l}}{dx_j}$, supplied by an extra equation from (\ref{DG1_extra}).
As is well-known, the cost of solving such a system directly is generally superlinear in the number of degrees-of-freedom, although recent
advances combining nested dissection ordering with low rank approximations can reduce this to near linear cost but with a significant prefactor
\cite{SchmitzYing3D}. The structure of (\ref{DG_matrix1}) can be exploited to rewrite it in the form:
\begin{displaymath}
\frac {dU^k}{dt} = V^k + G (U^{k-1},V^{k-1},U^k,V^k,U^{k+1},V^{k+1})^T ,
\end{displaymath}
where $G$ may be viewed as a lift matrix. As $G$ maps boundary data to volume data the cost of this operation will scale as
$\mathcal{O}((N+1)^{2d-1})$ per time step, with the larger cost of inverting $\hat{S}$ restricted to a precomputation. As $N$ may be much larger than $p$, this scaling implies that the resulting method would not be competitive with an implementation using a standard continuous Galerkin difference formulation of \cite{Gdiff} or the more recent method using a SIPG
formulation \cite{sipDGGD}. We must thus seek an improved method.

\subsection{Optimal Computational Complexity by Simultaneous Diagonalization}\label{GD_new_basis}
The above mentioned complexity for evolving $U$ is not competitive for practical computations unless we limit the element sizes, e.g.
with $N = \mathcal{O} (p^{(d-1)/d})$. In this section we follow  \cite{lynch1964direct} and show that it is possible to make a simple (computational) change of basis that results in a method with linear complexity. Precisely the new basis is found by solving the generalized eigenvalue of problem for each of the one dimensional matrices, $S_j$, $j = 1,\ldots,d$. That is, the new basis vectors are solutions to,
\begin{eqnarray*}
S_j\psi_{j,k_j} = \lambda_{j,k_j} M_j\psi_{j,k_j}, \ \ k_j=0,\ldots,N_j, \ \ j=1,\ldots,d.
\end{eqnarray*}
We normalize the eigenvectors according to
\begin{eqnarray*}
\psi_{j,k_j} \leftarrow \frac{\psi_{j,k_j}}{((\psi_{j,k_j})^TM_j\psi_{j,k_j})^{1/2}}.
\end{eqnarray*}

Let $\Psi_j $ be the matrix containing the new one dimensional basis
\[
\Psi_j = \left(
\begin{array}{cccc}
\psi_{j,0} & \psi_{j,1} & \cdots & \psi_{j,N_j}
\end{array}
\right) .
\]
Then the $d$-dimensional basis is
\begin{equation*}
\Psi = \otimes_{j=1}^d \Psi_j.
\end{equation*}
Now, we define $\bar{U}$ and $\bar{V}$ by
\[
U = \Psi \bar{U}, \ \ V = \Psi \bar{V},
\]
then equation (\ref{DG1}) and (\ref{DG2}) become
\begin{eqnarray}
&&\Psi^TS\Psi \left( \frac{d \bar{U}}{dt}-\bar{V} \right) = f_u, \label{DG1_sparse} \\
&&\Psi^TM\Psi \frac{d \bar{V}}{dt} +c^2\Psi^TS\Psi\bar{U}= f_v. \label{DG2_sparse}
\end{eqnarray}
where,
\begin{eqnarray*}
f_u &=& \sum_{j=1}^d
(\alpha-1) \Psi^T \left[(D^{R,R}-D^{L,L})\Psi \bar{V}^k
- D^{R,L}\Psi \bar{V}^{k+1}
+ D^{L,R}\Psi \bar{V}^{k-1}\right] \\
&& - \sum_{j=1}^d  \beta \Psi^T \left[ (C^{R,R}-C^{L,L}) \Psi \bar{U}^k
+ C^{R,L} \Psi \bar{U}^{k+1}
 - C^{L,R}\Psi \bar{U}^{k-1} \right],\\
f_v &=&c^2 \sum_{j=1}^{d}
(1-\alpha) \Psi^T \left( E^{R,R} -E^{L,L} \right)\Psi \bar{U}^k
+ \alpha \Psi^T  \left(E^{R,L}\Psi \bar{U}^{k+1}
- E^{L,R}\Psi \bar{U}^{k-1} \right) \\
&&- c^2\sum_{j=1}^{d} \tau \Psi^T\left[ \left( B^{R,R} - B^{L,L} \right) \Psi \bar{V}^k
+ B^{R,L}\Psi \bar{V}^{k+1}
- B^{L,R}\Psi \bar{V}^{k-1} \right].
\end{eqnarray*}

In the new basis we have that the mass matrix diagonalizes
\begin{equation*}
\Psi^T M \Psi =
\otimes_{j=1}^d (\Psi_j^T M_j \Psi_j)  = I,
\end{equation*}
as does the differentiation matrix
\begin{multline*}
\Psi^T S \Psi = \sum_{j=1}^d
(\Psi_1^T M_1 \Psi_1)
\otimes \cdots \otimes
(\Psi_{j-1}^TM_{j-1}\Psi_{j-1})
\otimes (\Psi_j^TS_j\Psi_j)\otimes(\Psi_{j+1}^TM_{j+1}\Psi_{j+1})\otimes
\cdots \\
\cdots
\otimes (\Psi_d^T M_d \Psi_d) =  \sum_j \Lambda_j,
\end{multline*}
where,
\begin{equation*}
    \Lambda_j = I
    \otimes
    \cdots
    \otimes I
    \otimes
    \text{diag}(\lambda_{j,0},\ldots,\lambda_{j,N_j})
    \otimes I
    \otimes
    \cdots
    \otimes I.
\end{equation*}
In the above equations we use $I$ to denote the identity matrix of size inferred by the context.

Note that as one of the eigenvalues of $S$ is zero one of the equations in (\ref{DG1_sparse}) vanishes. Suppose we have ordered the unknowns so that this corresponds to the first entry in $\bar{U}$, then we simply enforce the additional independent equation
\[
\frac{d \bar{U}_1}{dt} = \bar{V}_1.
\]

We thus conclude that the cost of all the volume terms scales linearly with the number of degrees of freedom. We now turn to the evaluation of the surface terms in the new basis.

Consider first the surface terms in the Galerkin difference basis. In a single dimension the elements in the four different surface terms are of the form
\begin{equation*}
\tilde{B}_{kl}^{X,Y} = \phi_{k}(X)\phi_{l}(Y), \ \
\tilde{D}_{kl}^{X,Y} = \frac{d\phi_{k}}{dx}(X)\phi_{l}(Y), \ \
\tilde{E}_{kl}^{X,Y} = \phi_{k}(X) \frac{d\phi_{l}}{dx}(Y), \ \
\tilde{C}_{kl}^{X,Y} = \frac{d\phi_{k}}{dx}(X)\frac{d\phi_{l}}{dx}(Y),
\end{equation*}
where $\{X,Y\}\in\{\{L,L\}, \{L,R\}, \{R,L\}, \{R,R\}\}$.

Now due to the local support properties of the Galerkin difference basis the number of nonzero elements in the above matrices are $1$ for $\tilde{B}^{X,Y}$, $(p+1)$ for $\tilde{D}^{X,Y}$ and $\tilde{E}^{X,Y}$, and $ (p+1)^2$ for $\tilde{C}^{X,Y}$.

The $d$-dimensional version of the surface matrices can again be constructed by Kronecker products. For example we have that
\begin{eqnarray*}
D_j^{X,Y} = M_1\otimes\cdots\otimes M_{j-1}\otimes\tilde{D}_j^{X,Y}\otimes M_{j+1}\otimes\cdots\otimes M_{d}.
\end{eqnarray*}

Applying the change of basis we have that
\begin{multline*}
\Psi^T D_j^{X,Y} \Psi = (\Psi_1^T M_1 \Psi_1)\otimes\cdots\otimes(\Psi_{j-1}^T M_{j-1}\Psi_{j-1})\otimes(\Psi_j^T \tilde{D}_j^{X,Y}\Psi_j)\otimes(\Psi_{j+1}^TM_{j+1}\Psi_{j+1})\otimes \cdots \\
\cdots \otimes(\Psi_d^T M_d \Psi_d)
= I\otimes\cdots\otimes I\otimes(\Psi_j \tilde{D}_j^{X,Y} \Psi_j)\otimes I\otimes\cdots\otimes I.
\end{multline*}
Thus applying $\Psi^T D_j^{X,Y} \Psi$ to the $(N+1)^d$ dimensional vector $\bar{V}$ can be done at a cost that scales with $ (p+1) (N+1)^d$. Similarly, the cost of applying   $\Psi^T B_j^{X,Y} \Psi$, $\Psi^T C_j^{X,Y} \Psi$, and $\Psi^T E_j^{X,Y} \Psi$ can be done at a cost of $ (N+1)^d$, $ (p+1)^2 (N+1)^d$, and $ (p+1) (N+1)^d$, respectively.

\graphicspath{{graphs_computational_cost/}}
\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.45\textwidth]{time_newbasis_dof}
\includegraphics[width=0.45\textwidth]{newbasis_change_cells}
\caption{Plots of CPU time (in seconds) in the Galerkin difference basis obtained by simultaneous diagonalization as a function of degrees of freedom (DOF) in two space dimensions. The graph on the left is for one DG element $(n = 1)$ and on the right for 36 DG elements $(n=6)$. The speed of sound is $c = 1$ and the splitting parameter for the upwind flux is $\xi = c$.\label{fig:timing_2d}}
\end{center}
\end{figure}

\subsection{Numerical Verification of the Computational Complexity}
We now present timing results that illustrate above analysis. We consider a two dimensional problem in the domain $[0,1]\times[0,1]$ and use the upwind flux (the other fluxes give similar timing results). For this problem we choose forcing and boundary conditions so that the solution is
\begin{equation*}
u(x,t) = \sin(16\pi t)\sin(16\pi x) + \cos(16\pi t)\cos(16\pi y).
\end{equation*}
We present results for both one DG element ($n=1$) and 36 DG elements ($n=6$). The degrees of freedom are $N+1$ for both the $x$ and the $y$ direction, i.e, the total degrees of freedom is $ n^2 (N+1)^2$. Specifically, we choose $N = 3,4,5\cdots,42$ when $n = 1$ and $N = 3,4,5\cdots,29$ when $n = 6$ for the improved basis. We also choose $p = 3$ which results in a fourth order accurate method.

To time the code we use the built-in function \verb|CPU_TIME()| in \verb|FORTRAN| to record the elapsed CPU time which is used to evolve the solution using the classic fourth order acurate Runge-Kutta method for $10$ timesteps steps.

In Figure \ref{fig:timing_2d}, we observe that the CPU time is proportional to $n^2(N+1)^2$ for the Galerkin difference basis obtained by the simultaneous diagonalization, which is orders of magnitude smaller than what would be required for the standard Galerkin basis.

\section{Dispersion Analysis}\label{dispersion}
To investigate how well the scheme proposed in Section \ref{GD_cost_subsection_1} preserves the wave propagation properties, we use the standard Bloch wave analysis as in \cite{ainsworth2004dispersive,AinsworthMonkMuniz,hu1999analysis,moura2015linear,hesthaven2007nodal}.  Here, we consider
\begin{equation}\label{wave_1d_problem}
\frac{\partial^2 u}{\partial t^2} = c^2\frac{\partial^2 u}{\partial x^2}, \ \ a<x<b,
\end{equation}
with initial condition $u(x,0) = e^{ikx}$ and periodic boundary conditions $u(a,t) = u(b,t)$, $\frac {\partial u}{\partial x}(a,t)= \frac {\partial u}{\partial x}(b,t)$.
We then seek spatially periodic solutions of the form
\begin{equation*}\label{solution_dispersion}
u(x,t) = e^{i(\kappa x-wt)},
\end{equation*}
from which the exact dispersion relation $\omega = \pm c\kappa$ for (\ref{wave_1d_problem}) can be found.

Next, partition the computational domain into non-overlapping uniform DG elements $I^k = [x^k,x^{k+1}]$, $k = 0,\cdots,n-1$ with $H = x^{k+1}-x^k = (b-a)/n$. For each DG element $I^k$, there are $N+1$ equidistant nodal degrees of freedom with spacing $h = H/N$. The semi-discretization in $I^k$ becomes then follows from (\ref{DG_matrix1}) and (\ref{DG_matrix2}) with $d=1$. Let the vectors $U^k = (U_0^k,U_1^k,\cdots,U_{N}^k)$ and $V^k = (V_0^k,V_1^k,\cdots,V_{N}^k)$ hold the nodal approximations of $u(x,t)$ and $v(x,t)$ in $I^k$, respectively.

We seek solutions in terms of Bloch waves
\begin{equation}\label{dispersion_transfer}
U^k_l = \hat{U}^k_le^{i(\kappa (x^k+lh)-\omega t)},\quad V^k_l = \hat{V}_l^ke^{i(\kappa (x^k+lh)-\omega t)},\ \ l=0,\cdots,N,
\end{equation}
and thus assume periodicity of the solution. That is
\begin{equation}\label{periodic_dispersion}
W^{k+1}_0 = e^{i\kappa H}W^k_0, \quad W^{k-1}_N = e^{-i\kappa H}W^k_N,
\end{equation}
where $W$ represents either $U$ or $V$. To condense the notation, we have omitted the superscript $k$ for the rest of this section. Let $Z = (U,V)^T$ and combine with (\ref{DG_matrix1})-(\ref{DG_matrix2}), (\ref{dispersion_transfer})-(\ref{periodic_dispersion}) to obtain the following eigenvalue problem
\begin{equation*}
\hat{A}Z = -i\Omega Z,\ \ \hat{A} = \hat{A}_1\hat{A}_2,
\end{equation*}
where
\begin{equation*}
\Omega = \frac{\omega H}{c},\quad \hat{A}_1 = \frac{H}{c}
\begin{pmatrix}
\hat{S}^{-1}&~\\
~&M^{-1}
\end{pmatrix},\ \
\hat{A}_2 =
\begin{pmatrix}
A_2^{11}&A_2^{12}\\
A_2^{21}&A_2^{22}
\end{pmatrix},
\end{equation*}
with
\begin{eqnarray*}
	A_2^{11} &=& \tau (C^{L,L} - C^{R,R})-\tau (e^{iK}C^{R,L}-e^{-iK}C^{L,R}),\\
	A_2^{12} &=& \hat{S}+ (-\theta)(D^{R,R}-D^{L,L})+\theta(e^{iK}D^{R,L}-e^{-iK}D^{L,R}),\\
	A_2^{21} &=&-c^2S+ \theta c^2(E^{R,R}-E^{L,L})+(1-\theta) c^2(e^{iK}E^{R,L}-e^{-iK}E^{L,R}),\\
	A_2^{22}&=&-\beta c^2(B^{R,R}-B^{L,L})-\beta c^2(e^{iK}B^{R,L}-e^{-iK}B^{L,R}),
\end{eqnarray*}
and $K = \kappa H$. Note that the values of $\Omega = \Omega_r + i\Omega_i$ are, in general, complex valued. Here $\Omega_i$ is non-positive and represents the numerical damping of the corresponding scheme, and the real part $\Omega_r$ is the approximation to  $\frac{\omega H}{c}$.

For the numerical simulations in this section, the computational domain is chosen to be $x\in[0,1]$. The order of the approximation space is set to be $p = 3$. The number of degrees of freedom in each DG element is $N+1 = 10$.
\graphicspath{{graphs_dispersion/}}
\begin{figure}[htb]
	\begin{center}
		\includegraphics[width=0.49\textwidth]{central_q3_r_2dg}
		\includegraphics[width=0.49\textwidth]{alternating_q3_r_2dg}
	\end{center}
	\caption{In this simulation, we choose $n = 2$ DG elements and $N = 9$ GD cells in each DG element. On the left we show the numerical dispersion relation for the central flux. On the right we show the numerical dispersion relation for the alternating flux. Black dot lines are the exact dispersion relation, p.-mode represents physical modes, dashed lines are for spurious modes.\label{dispersion_dissipation_c_2ele}}
\end{figure}
\begin{figure}[htb]
	\begin{center}
		\includegraphics[width=0.49\textwidth]{central_q3_i_2dg}
		\includegraphics[width=0.49\textwidth]{alternating_q3_i_2dg}
	\end{center}
	\caption{In this simulation, we choose $n = 2$ DG elements and $N = 9$ GD cells in each DG element. On the left we show the numerical dissipation for the central flux for all modes. On the right we show the numerical dissipation for the alternating flux for all modes. 'p.-mode' represents the physical, dashed lines are for spurious modes.\label{dispersion_dissipation_a_2ele}}
\end{figure}

\begin{figure}[htb]
	\begin{center}
		\includegraphics[width=0.49\textwidth]{upwind_q3_r_2dg}
		\includegraphics[width=0.49\textwidth]{upwind_q3_i_2dg}
	\end{center}
	\caption{In this simulation, we choose $n = 2$ DG elements and $N = 9$ GD cells in each DG element. On the left, we show the numerical dispersion relation for the upwind flux. Black dot lines are the exact dispersion relation,  p.-mode represents the physical and dash lines are for spurious modes. On the right, we illustrate the dissipation associated with the twenty modes.\label{dispersion_dissipation_up_2ele}}
\end{figure}

Figure \ref{dispersion_dissipation_up_2ele} presents the dispersion relation of the upwind flux. When $K$ is small, the numerical phase velocity also reflects the physical wave speed. Comparing the results in Figure \ref{dispersion_dissipation_up_2ele} and Figure \ref{dispersion_dissipation_c_2ele}--\ref{dispersion_dissipation_a_2ele}, we find that both conservative (central flux and alternating flux) and dissipative schemes (upwind flux) recover the physical mode for small $K$. The conservative schemes admit more complicated phenomena: the spurious modes do not damp for small values of $K$; for the dissipative scheme, however, the unphysical modes are strongly damped.

\begin{figure}[htb]
	\begin{center}

		\includegraphics[width=0.49\textwidth]{dispersion_different_p}
		\includegraphics[width=0.49\textwidth]{dissipation_different_p}
	\end{center}
	\caption{On the left, we show the numerical dispersion relation for the upwind flux. The black dashed line is the exact dispersion relation. We present the numerical dispersion relation for the physical mode at different orders $p = (1,3,5,7,9)$ with $N = 19$. On the right, we illustrate the dissipation associated with the different orders $p = (1,3,5,7,9)$ with $N = 19$.\label{dispersion_dissipation_up_dq}}
\end{figure}
In Figure \ref{dispersion_dissipation_up_dq}, we show the dispersion relation of physical modes of the dissipative scheme (upwind flux) for a range of orders of approximation $p = (1,3,5,7,9)$ and $N = 19$ for all different approximation degrees $p$. We see that the numerical phase velocity is very close to the physical wave speed when $K$ is small and improves for a broader range of $K$ as the order of the approximation increases.

Lastly we consider direct comparisons of the proposed method with both alternating and upwind flux choices to the
interior penalty discontinuous Galerkin method (IPDG). We also consider the effect of increasing the number of nodes
within the element. The dispersion relation for the IPDG method is discussed extensively in \cite{AinsworthMonkMuniz}, and we
use the Bloch wave formulation and penalty parameter suggested there. On the left in Figure \ref{discomp} we take $p=3$ and
apply the method proposed here with $9$ and $39$ cells per element along with IPDG. In all cases we use a single
element and compute the dispersion
relation for a range of wave numbers with the largest wave number corresponding to $5$ degrees-of-freedom per wavelength. Note
that this implies a significant disparity in the wave numbers considered. Here we see that the dispersion error for the upwind
method is significantly smaller than for the conservative alternating flux scheme, but keep in mind that the graph does not
take the dissipation error of the upwind discretization into account. Both our conservative method and IPDG have oscillatory
dispersion relations, but perhaps because of the larger wave number range involved we see more oscillations for the proposed
method. Increasing the number of nodes per element improves the results. Here IPDG is the best performer at the coarsest
discretization level, though we do recall that the proposed scheme admits larger time steps than IPDG.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.49\textwidth]{disp3k}
\includegraphics[width=0.49\textwidth]{disp7k}
\end{center}
\caption{On the left we compare the dispersion errors with $p=3$
for the alternating and upwind fluxes with $10$ and $40$ points per element as well as for IPDG scheme. On the right we show
the same for $p=7$ where the alternating and upwind fluxes are considered for $20$ and $40$ points per element.\label{discomp}}
\end{figure}

On the right we display the analogous plot for $p=7$ and largest wave number corresponding to $3$ degrees-of-freedom per
wavelength. Now the results for the proposed methods are shown with $19$ and $39$ cells per element. The comparisons between
them follow the same pattern as in the case $p=3$. However in this case the dispersion errors are larger for IPDG at the
coarse resolutions.

\section*{Conflict of interest}

The authors declare that they have no conflict of interest.

\section{Summary}\label{conclusion}
In conclusion, we have demonstrated the energy-based DG method with Galerkin difference basis functions for second-order wave equations. In particular:

\begin{itemize}
\item[a)] We derived a new basis by simultaneous diagonalization of the mass and the stiffness matrices from the Galerkin difference basis functions. The new basis reduces the computational cost of evolving the solution from superlinear complexity with respect to degrees of freedom to optimal linear complexity.
\item[b)] Using Bloch wave analysis we computed the dispersive and dissipative properties of the method.
\item[c)] By numerical experiments we showed that the spectral radius of the semi-discretization of our scheme is linearly proportional to the degree of the approximation space $p$. This translates to the ability to march in time using $p$ times  larger time steps compared with traditional element based methods such as spectral-, continuous- and discontinuous finite elements.
\item[d)] Optimal convergence was observed for problems in both one and two dimensions for all numerical fluxes when the degree $p\geq3$. The results apply both when
the number of points-per-element is refined with the number of elements fixed and when the number of points-per-element is fixed and the number of elements is
increased.
\item[e)] We illustrated that the method is not dramatically slower for variable coefficient problems if the mass matrices are inverted using the preconditioned conjugated gradient method. In this case fixing the number of points-per-element while increasig the number of elements may be the most efficient refinement
strategy.
\end{itemize}

\section{Convergence}\label{convergence}
In this section, we present numerical results to investigate the convergence of our method in both the $L^2$ norm, $\|u-u^h\|_{L^2}$, and the energy norm, $\left(\|\nabla(u - u^h)\|_{L^2}^2 + \|v-v^h\|^2\right)^{1/2}$. We consider both one dimensional and two dimensional problems. As we have discussed in earlier papers \cite{DATH_UP,appelo2020energy,zhang2019energy}, the energy norm is the starting point for theoretically establishing stability and rates of convergence of the energy DG method, and the results there apply directly to the proposed scheme. In all experiments below we use a nodal formulation associated with the basis functions in Section \ref{GD_basis} and march in time by the classic fourth order accurate Runge-Kutta method. We choose the speed of sound to be $c = 1$ and the flux splitting parameter $\xi$ in the upwind scheme to be $c$.

\begin{figure}[htb!]
	\begin{center}
		\includegraphics[width=0.48\textwidth]{l2_u_u1d}
		\includegraphics[width=0.48\textwidth]{en_u_1d}\\
		\includegraphics[width=0.48\textwidth]{l2_a_u1d}
		\includegraphics[width=0.48\textwidth]{en_a_1d}\\
		\includegraphics[width=0.48\textwidth]{l2_c_u1d}
		\includegraphics[width=0.48\textwidth]{en_c_1d}
		\caption{Plots of the $L^2$ norm error $\|u-u^h\|_{L^2}$ and the energy norm error $\left(\|\nabla(u-u^h)\|_{L^2}^2+\|v-v^h\|_{L^2}^2\right)^{1/2}$ as functions of $h$ in one dimension with upwind flux, central flux and alternating flux for periodic boundary conditions. In the legend, $p$ is the degree of the approximation space of $u$ and $v$. The filled circles are the results for the fixed degrees of freedom $N+1 = 11$ and the hollow circles are the results for the fixed number of DG elements $n = 10$.\label{fig:uv_1d_periodic}}
	\end{center}
\end{figure}
\subsection{Periodic Boundary Conditions in One Dimension} \label{numerical_1d_periodic}
To investigate the order of accuracy of our method, we solve
\begin{equation*}
\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}, \quad x\in(0,1),\quad t\geq0,
\end{equation*}
with periodic boundary conditions $u(0,t) = u(1,t)$ and with initial data so that the solution is the exact traveling wave
\begin{equation}\label{1d_sol_conv}
u(x,t) = \sin(8\pi(x-ct)).
\end{equation}
The discretization is performed on a uniform mesh with DG element vertices $x_i = iH$, $i = 0,\cdots,n$, $H = 1/n$. The vertices of the subcells for the $i$-th DG element are $x_{ij} = x_i + jh$ with $i = 0,\cdots,n-1$, $j = 0,\cdots,N$ and $h = H/N$. The problem is evolved until the final time $T = 1.075$ with the time step size $\Delta t = \mbox{CFL}\times h$, $\mbox{CFL} = 0.075/(2\pi)$ to guarantee the error is dominated by the spatial error.

We present results for the degree of the approximation space of $u^h, v^h$ being $p = (1,3,5,7)$. The mesh size $h$ is refined by either increasing the number of DG elements $n$ and fixing the degrees of freedom $N+1$ in each DG element or increasing the degrees of freedom $N+1$ in each DG element and fixing the number of DG elements $n$.

The $L^2$ norm error $\|u-u^h\|_{L^2}$ and the energy norm error $\left(\|\nabla(u-u^h)\|_{L^2}^2+\|v-v^h\|_{L^2}^2\right)^{1/2}$ are plotted against the grid spacing $h$ in Figure \ref{fig:uv_1d_periodic} with the upwind flux, the alternating flux and the central flux, respectively. Linear regression estimates of the rate of convergence can be found in Table \ref{convergence_rate_1d_Nchange} for fixed number of DG elements $n = 10$ and in Table \ref{convergence_rate_1d_Nfixed} for fixed degrees of freedom $N+1 = 11$ in each DG element. Note that we use the same (on element) mesh size $h = 1/n/N$ for these two cases when we generate the results in Figure \ref{fig:uv_1d_periodic}. From Figure \ref{fig:uv_1d_periodic}, we observe that the $L^2$ error $\|u-u^h\|_{L^2}$ becomes oscillatory when it reaches $10^{-13}$ for $p = 7$, thus we only use the first $31$ data points to estimate the convergence rate for this case. For the other cases we use all data to estimate the convergence rate. Generally speaking, the errors are comparable for the two modes of refinement ($N$ fixed and $n$ fixed), with slightly smaller errors if we fix $n$ and refine within each element. From the comparisons of the dispersion errors we do see that refinement within an element is likely to be preferable in most cases. Note that the existing theory for the method in one space dimension only establishes convergence in the energy norm and only proves optimal convergence for the upwind flux. That is, we only have proofs of convergence at the optimal order $p$ in the energy norm for the upwind flux; for the other fluxes the existing theory only guarantees a rate of $p-1$ and in particular it does not guarantee convergence for the conservative fluxes when $p=1$. Specifically we observe the following, which in most cases is better than what we can prove.

\begin{itemize}
\item[a).] When the number of DG elements $n = 10$ is fixed, from Table \ref{convergence_rate_1d_Nchange}, we observe convergence at rates exceeding $(p+1)$ in the $L^2$ norm for $p = (3,5,7)$ with all three choices for the flux; $p$-th order convergence for the alternating flux and optimal convergence $(p+1)$ for the upwind flux and the central flux when $p = 1$. In the energy norm, we again observe convergence at rates exceeding $p$ for both the upwind flux and the central flux with $p = (1,3,5,7)$, but suboptimal convergence for the alternating flux with $p = (3,5,7)$. We note that we do not expect the asymptotic convergence rates to exceed $p+1$ in $L^2$ and $p$ in the energy norm. However, as we refine within each element, the effect of the $2p$th
order interior formulas is felt, and it is possible to observe convergence at higher rates for some range of
resolutions. This effect is observed for the continuous Galerkin difference methods in \cite{Gdiff}.
\item[b).] When the degrees of freedom per element $N+1 = 11$ is fixed, from Table \ref{convergence_rate_1d_Nfixed}, we note optimal convergence $(p+1)$ in the $L^2$ norm for $p = (3,5,7)$ with all three fluxes; there is order reduction for $p = 1$. From the energy norm, we have $p$-th order convergence for $p = (1,3,5,7)$ with both the upwind flux and the central flux; for the alternating flux, we observe the suboptimal $(p-1)$ convergence rate for all $p$, and in particular no
convergence when $p=1$.
\end{itemize}

\input{manuscript_table/periodic_1d.tex}

\subsection{Dirichlet Boundary Conditions in Two Dimensions}
\graphicspath{{graphs_convergence_2d/}}
\begin{figure}[htb!]
	\begin{center}
		\includegraphics[width=0.48\textwidth]{l2_u_2d_u}
		\includegraphics[width=0.48\textwidth]{energy_uv_2d_u}\\
		\includegraphics[width=0.48\textwidth]{l2_u_2d_a}
		\includegraphics[width=0.48\textwidth]{energy_uv_2d_a}\\
		\includegraphics[width=0.48\textwidth]{l2_u_2d_c}
		\includegraphics[width=0.48\textwidth]{energy_uv_2d_c}
		\caption{Plots of the $L^2$ norm error $\|u-u^h\|_{L^2}$ and the energy norm error $\left(\|\nabla(u-u^h)\|_{L^2}^2+\|v-v^h\|_{L^2}^2\right)^{1/2}$ as functions of $h$ in two dimensions with upwind flux, central flux and alternating flux for the Dirichlet boundary condition. In the legend, $p$ is the degree of the approximation space of $u$ and $v$. The filled circles are the results for the fixed degrees of freedom $(N+1)^2 = 11^2$ and the hollow circles are the results for the fixed number of DG elements $n^2 = 10^2$.}\label{fig:uv_2d_dirichlet}
	\end{center}
\end{figure}
In this section, we examine the rate of convergence for our scheme applied to the following two dimensional problem
\begin{equation*}
\frac{\partial^2 u}{\partial t^2} = \left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right), \ \ \ (x,y)\in(0,1)\times(0,1),\ \ t\geq 0,
\end{equation*}
with initial data and boundary conditions chosen such that the exact solution is
\begin{equation*}
u(x,y,t) = \cos(15\sqrt{2}\pi ct)\sin(15\pi x)\sin(15\pi y).
\end{equation*}

The computational domain is discretized into Cartesian DG elements whose vertices are $x_i = iH$, $y_j = jH$, $i,j = 0,1,\cdots n$, $H = 1/n$. The vertices of subcells for the $ij$-th DG element are $x_{ik} = x_i + kh$ and $y_{jl} = y_j + lh$ with $i,j = 0,\cdots,n-1$ and $k,l = 0,\cdots,N$,  $h = H/N$. The problem is evolved until final time $T = 0.2$ and the time step size is $\Delta t = \mbox{CFL} h$, with $\mbox{CFL} = 0.075/(2\pi)$ to guarantee the temporal error is dominated by the spatial error.

As in the one dimensional test in Section \ref{numerical_1d_periodic}, we consider three different numerical fluxes: the upwind flux, the alternating flux, and the central flux. We again use two different ways to refine the spatial mesh size: one is to fix the degrees of freedom in each element, $(N+1)^2 = 11^2$ and change the number of DG elements $n^2$. The other is to fix the number of DG elements $n^2 = 10^2$ and change the degrees of freedom in each DG element. Again, we have the same (on element) mesh size $h = 1/n/N$ for these two cases.

\input{manuscript_table/Dirichlet_2d.tex}

The $L^2$ norm error $\|u-u^h\|_{L^2}$ and the energy norm error $\left(\|\nabla(u-u^h)\|_{L^2}^2+\|v-v^h\|_{L^2}^2\right)^{1/2}$ are presented in Figure \ref{fig:uv_2d_dirichlet} for the upwind flux, the central flux and the alternating flux, respectively. The corresponding convergence rates from linear regression are shown in Table \ref{convergence_rate_2d_Nchange} for a fixed number of DG elements $n^2 = 10^2$ and in Table \ref{convergence_rate_2d_Nfixed} for fixed degrees of freedom $(N+1)^2 = 11^2$ on each DG element. We use all data to generate the convergence rate here, but for the alternating flux with fixed degrees of freedom $(N+1)^2 = 11^2$, we use the data from the $25$ coarsest grids which excludes the outliers where the error is very small. Generally speaking, the results are similar to the one dimensional results in Section \ref{numerical_1d_periodic}. We observe optimal convergence when $p = (3,5,7)$ for all cases. However, we observe a rate of convergence $p+1$ in the $L^2$ error norm for the upwind flux when $p = 1$ with fixed degrees of freedom $(N+1)^2 = 11^2$ on each DG element. This is slightly better than the corresponding one dimensional result.

\section{Introduction}

Discontinuous Galerkin methods have become a method of choice for solving first order hyperbolic equations in Friedrichs form \cite{hesthaven2007nodal}. They possess many desirable properties such as arbitrary order, robustness, geometric flexibility and explicit time evolution.  Analogous methods for second order hyperbolic equations are less well established, despite the fact that many governing equations arising in physics are in second order form. Even though it is often possible to rewrite second order hyperbolic equations in first order form, the first order formulation has some drawbacks. It almost always needs more
boundary conditions and it is only equivalent to the original second order equation for constrained data. In addition, not all second order hyperbolic equations can be rewritten as a first order system which is in Friedrichs form. In \cite{DATH_UP} Appel{\"o} and Hagstrom  proposed an energy-based DG method for second order wave equations. The idea in \cite{DATH_UP} is to introduce a new variable $v = \frac{\partial u}{\partial t}$ to transfer second order hyperbolic equations to first order systems in time only, and then seek approximations which satisfy a discrete energy equality. The method features a direct, mesh-independent approach to defining interelement fluxes. Both energy-conserving and upwind discretizations have been devised and their extension to elastic, advective and semi-linear wave equations can be found in \cite{appelo2018energy,appelo2020energy,zhang2019energy}.

Galerkin difference methods were introduced by Banks and Hagstrom \cite{Gdiff} to solve hyperbolic initial-boundary value problems. The idea is to use a Galerkin construction to derive energy stable finite difference methods. The basis functions are Lagrange functions associated with continuous piecewise polynomial approximation on a computational grid. Salient features of these methods are: the discrete approximations are uniform from grid-point to grid-point when away from domain boundaries; no new degrees of freedom are added when the approximation order increases; they do not require significantly smaller time steps as the order increases. Comparing them to summation-by-parts (SBP) difference schemes \cite{SBPrev}, they have the advantage of being directly constructable at arbitrary order and can be seamlessly interfaced with standard schemes on unstructured grids
\cite{GDgeo}. The relative disadvantage is that Galerkin difference operators typically require twice as many flops as SBP operators of the same order.

In this paper, we combine the energy based discontinuous Galerkin methods with the Galerkin difference methods to solve the second order wave equation. The corresponding mass and stiffness matrices are banded matrices because of the properties of the Galerkin difference basis functions. For high dimensional problems on structured grids, the Galerkin basis functions are the tensor product of the Galerkin basis functions in one dimension. Then the mass matrix is a Kronecker product of the mass matrices in each dimension and the stiffness matrix is a summation of the Kronecker product of the mass and stiffness matrices in each dimension. This fact indicates that the inversion of the mass matrix can be computed with a computational cost which is linear with respect to the total degrees of freedom, while the computational cost grows rapidly for the inversion of the stiffness matrix. A result of this paper is the application of the simultaneous diagonalization technique from \cite{lynch1964direct} to derive a new class of Galerkin difference basis functions to reduce the computational cost of the inversion of the stiffness matrix.

We note that the literature on high order methods for wave equations in second order form is extensive. We will not try to review all methods here, but rather mention a few that are representative of the state of the art. In the class of discontinuous Galerkin methods it is worth mentioning the symmetric interior penalty (SIPG) method of \cite{GSSwave},
the local discontinuous Galerkin (LDG) method of \cite{ChouShuXing2014} and the nonsymmetric interior penalty
method of \cite{IPDG_Elastic}. Finite difference methods include those using SBP operators \cite{Mattsson2004,Mattsson2012} as well as those using upwind methods \cite{BANKS20125854}.

The rest of the paper is organized as follows. In Section~\ref{preli}, we introduce the construction of the Galerkin difference basis functions and the energy based discontinuous Galerkin method for the second order scalar wave equation. We derive the new basis functions from the Galerkin difference basis functions and investigate the reduction in computational cost in Section~\ref{GD_cost}. In Section~\ref{dispersion}, we present the dispersion properties of the scheme. Section~\ref{spectral_radius} shows the spectral radii of the proposed scheme. Numerical experiments that illustrate optimal convergence in both $L^2$ and energy norms are given in Section~\ref{convergence}. In Section \ref{variable}, we apply the method to  a problem with variable sound wave speed. Here the new basis construction does not apply, and as an alternative we use preconditioned conjugate gradient iterations. Finally, our conclusions summarized in Section~\ref{conclusion}.

\subsection{The Energy Based Discontinuous Galerkin Method for the Wave Equation}\label{GD_cost_subsection_1}
We consider a mesh that discretizes $\Omega$ into non-overlapping box shaped elements $\Omega^{k}$ with $\Omega = \bigcup_{k=1}^{m} \Omega^k$. On each element $\Omega^k = \otimes_{l \in \mathcal{S}}^d[L_{l}^k,R_{l}^k]$,  where $\mathcal{S} = \{x,y\}$ if $d=2$ and $\mathcal{S} = \{x,y,z\}$ if $d=3$. Let $(\mathbb{Q}_{p,N})^d$ be the space of functions spanned by the tensor product of the one dimensional Galerkin difference basis on an element. Then a test function $\varphi$ in $(\mathbb{Q}_{p,N})^d$ can be expressed (in three dimensions) as
\begin{eqnarray*}
\varphi(x,y,z)= \prod_{\rho \in \mathcal{S}} \phi_{p,i_\rho}(\rho).
\end{eqnarray*}
Here we assume the same degree of approximation and the same number of degrees of freedom in each dimension, but remark that these can also be chosen independently.

Now following \cite{DATH_UP} we define our discretization by the element-wise variational statement
\begin{eqnarray} \label{DG1}
\int_{\Omega^k} c^2\nabla \varphi_u \cdot \left( \frac{\partial \nabla u}{\partial t}-\nabla v \right) d \Omega^k &=& \int_{\partial \Omega^k} c^2(\vec{n} \cdot \nabla \varphi_u) (v^\ast-v)
 \, d \Omega^k, \\
 \label{DG2}
\int_{ \Omega^k} \varphi_v \frac{\partial v}{\partial t}+c^2 \nabla \varphi_v \cdot\nabla u \, d \Omega^k &=&  \int_{\partial \Omega^k} c^2 \varphi_v \left( \vec{n} \cdot ( \nabla u)^\ast  \right) d \Omega^k,
\end{eqnarray}
for all $(\varphi_u,\varphi_v) \in (\mathbb{Q}_{p,N})^d \times (\mathbb{Q}_{p,N})^d$. As equation (\ref{DG1}) vanishes for constants we augment it by the independent equation
\begin{equation}
\int_{\Omega^k} \left( \frac{\partial u}{\partial t}-  v \right) d \Omega^k =0. \label{DG1_extra}
\end{equation}

\subsubsection{Numerical Fluxes}
Following the  notation in \cite{DATH_UP} we take $\beta,\tau\geq0$, $0\leq\alpha\leq 1$. Now, let the superscript $1$ represent data from inside an element and the superscript $2$ represent data from the outside of an element. Then we can write the numerical fluxes as
\begin{eqnarray*}
v^\ast &=& \alpha v^1+(1-\alpha)v^2-\beta(\nabla u^1\cdot {\bf n}^1 + \nabla u^2\cdot{\bf n}^2),\\
(\nabla u)^\ast &=& (1-\alpha)\nabla u^1+\alpha\nabla u^2-\tau(v{\bf n}^1+v{\bf n}^2).
\end{eqnarray*}
The above fluxes are energy-conserving when $\beta=\tau=0$, and upwind and dissipative when $\beta, \tau>0$. In the rest of the analysis, we focus on the three popular choices:
\begin{itemize}
\item[] {\it Central\ \ flux}: $\alpha = \frac{1}{2}$, $\beta = \tau = 0$,
\item[] {\it Alternating\ \ flux}: $\alpha = 1$, $\beta = \tau = 0$,
\item[] {\it Upwind\ \ flux}: $\alpha = \frac{1}{2}$, $\beta = \frac{\xi}{2}$, $\tau = \frac{1}{2\xi}$.
\end{itemize}
In the last flux $\xi$ is a flux splitting parameter with the same dimensional units as the speed of sound $c$.

We note that the possibility of choosing simple mesh-independent flux parameters is a feature of the energy-DG formulation.

\section{Preliminaries}\label{preli}
We consider the wave equation in first order form in time and second order form in space
\begin{eqnarray}
u_t &=& v, \label{eq:wave_1} \\
v_t &=& \nabla \cdot (c^2 \nabla u) + f, \ \ t>0, \ \ (x,y,z) \in \Omega, \label{eq:wave_2}
\end{eqnarray}
complemented with initial data and boundary conditions.

In what follows we will consider energy based discontinuous Galerkin methods implemented on $d$-dimensional tensor product elements. Each element will be mapped to the $d$-dimensional unit cube which will be discretized by an equidistant Cartesian grid. Precisely,
in each dimension $k$ and for each element we consider a grid discretizing the unit cube
\[
x_{k,i} = i h_k, \ \ i = 0,\ldots,N_k, \ \ h_k = 1/N_k.
\]
For simplicity here we will use the same number of intervals in each dimension, $N_k=N$. In the following $n$ denotes the number of elements in one dimension, $m$ denotes the number of elements in multiple dimensions and $p$ denotes the polynomial degree of the approximation. We will only consider $p$ odd, which corresponds to the {\em locally continuous} difference basis and will often use the integer $q=(p+1)/2$.

\section{Spectral Radius}\label{spectral_radius}
In this section, we study the spectral radius of the semi-discretization of our scheme. Consider the problem in one dimension
\[
\frac{\partial^2 u}{\partial t^2} = c^2\frac{\partial^2 u}{\partial x^2},\quad x\in[a,b],
\]
with either periodic boundary conditions or a homogeneous Dirichlet boundary condition at the left boundary and a homogeneous Neumann boundary condition at the right boundary. We emphasize that an expected advantage of Galerkin difference methods compared with standard DG schemes is a milder growth in the spectral radius
as the order is increased, and the results below confirm this expectation.

As above, the computational domain $[a,b]$ is divided into non-overlapping uniform DG elements $I^k = [x^k, x^{k+1}], k = 1,\cdots,n$ with element size $H = (b-a)/n$. Each DG element $I^k$ is partitioned into $N$ equidistant subcells with cell size $h = H/N$. Namely, we have $N+1$ degrees of freedom for each DG element. Then the displacement $u(x,t)$ and the velocity $v(x,t)$ in $I^k$ are approximated by the nodal values $U_l^k$ and $V_l^k$ with
\[u(x,t) = \sum_{l=0}^N U_l^k\phi_{p,l}(x),\quad v(x,t) = \sum_{l=0}^N V_l^k\phi_{p,l}(x),\]
respectively. Here, $\phi_{p,l}$ are $p$-th order Galerkin difference basis functions. In this experiment the speed of sound is $c = 1$ and we choose the splitting parameter in the upwind flux to be $\xi = c$.

Let the vectors $U$ and $V$ contain the nodal values $U_l^k$ and $V_l^k$
\[
U = [U_0^1,\cdots,U_N^1,\cdots,U_0^n,\cdots,U_N^n], \quad V = [V_0^1,\cdots,V_N^1,\cdots,V_0^n,\cdots,V_N^n].
\]
Then the semi-discretization (\ref{DG_matrix1})-(\ref{DG_matrix2}) can be written as a system of ordinary differential equations
\begin{equation*}
\frac{d{Z}}{dt} = \mathcal{L}Z, \quad Z = (U,V)^T.
\end{equation*}
In the experiments in this section, the computational domain is chosen to be $x\in[0,1]$, the number of DG elements is set to be $n = 1$, the degree of approximation space is given by $p = (1,3,5,7,9,11)$, and the number of degrees of freedom in the DG element is given by $N+1 = (31,61,121)$.

\graphicspath{{graphs_spectral_radius/}}
\begin{figure}[h]
	\begin{center}
		\includegraphics[width=0.48\textwidth]{spectral_radius_all}
		\includegraphics[width=0.48\textwidth]{spectral_radius_dbc}
	\end{center}
	\caption{Spectral radius as a function of the approximation degree $p$ for three different fluxes. The figures are for periodic boundary conditions (left) and Dirichlet and Neumann conditions (right).}\label{spectral_radii_vs_degree_new}
\end{figure}

Denote the central flux by C.-flux, the alternating flux by A.-flux and the upwind flux by U.-flux. Figure \ref{spectral_radii_vs_degree_new} displays the amplitude of the largest eigenvalue as a function of the degree $p$ for three different values of $N$. The subfigures display results for periodic and non-periodic boundary conditions. We observe that the spectral radius is linearly proportional to the degree $p$ for all three different numerical fluxes. This is in contrast to standard discontinuous finite elements where the spectral radius grows quadratically with $p$. As a consequence the method proposed here can march in time with $p$ times larger time steps.

\section{Problems with Variable Coefficients}\label{variable}
As was demonstrated in Section \ref{GD_cost} when applied to constant coefficient problems on Cartesian meshes the proposed method has the same complexity as a traditional finite difference method. Unfortunately the simultaneous diagonalization cannot be expected to work when the speed of sound varies in space or for a constant coefficient problem on non-Cartesian meshes (in the latter case the transformation from a physical element to the reference element will result in a variable coefficient problem). For such variable coefficient problems we will stay with the standard Galerkin difference basis and invert the matrices on the left of the element-wise equations (\ref{DG_matrix1}) and (\ref{DG_matrix2}) using the preconditioned conjugate gradient (pcg) method. In a first experiment we demonstrate that the number of pcg iterations needed in each time step is small. In a second experiment we demonstrate the ability of the method to compute the solution to a more complex application-type problem.

\subsection{Efficiency of PCG for Inverting Mass Matrices}\label{iterative_variable}
Here, we consider the second order wave equations with a variable coefficient $c^2(x,y)$ in two dimensions as follows
\begin{equation}\label{variable_problem}
\frac{\partial^2 u}{\partial t^2} = \nabla\cdot(c^2(x,y)\nabla u) + f(x,y), \quad (x,y)\in[0,1]\times[0,1],\quad t\geq0,
\end{equation}
where $c^2(x,y) = 1+x^2+y^2$. The initial conditions and the external forcing function are determined by the manufactured solution
\[
u(x,y,t) = \sin(8\sqrt{2}\pi t)\sin(8\pi x)\sin(8\pi y).
\]
We impose periodic boundary conditions.

The key point for the success of the new basis functions proposed in Section \ref{GD_new_basis} is the mass matrix and stiffness matrix in  d-dimensions are constructed by the tensor products of the corresponding mass and stiffness matrices in one dimension. From the scheme (\ref{DG1})-(\ref{DG2}), though the mass matrix for the problem (\ref{variable_problem}) still keeps the tensor product form, the elements of the stiffness matrix $S_{c^2(x,y)}$ are derived from
\[
\int_{\Omega_k} c^2(x,y)\nabla\varphi_u\cdot\frac{\partial \nabla u}{\partial t} d\Omega^k,
\]
which does not have a tensor product form for general $c^2(x,y)$. In \cite{chan2017weight,chan2017weight2}, the authors proposed weight-adjusted discontinuous Galerkin (WADG) method to handle the variable coefficient in the mass matrix. The idea there is to replace the weighted $L^2$ inner product with a weight-adjusted inner product. Unfortunately this approach is not applicable here.

As an alternative for the variable coefficient $c^2(x,y)$ problem, we compute the time derivatives in the scheme (\ref{DG_matrix1})--(\ref{DG_matrix2}) iteratively by the preconditioned conjugate gradient. As a preconditioner we use the zero fill-in incomplete Cholesky factorization of $S_{c^2(x,y)}$ and $M$ as the preconditioning matrix of the system (\ref{DG_matrix1}) and (\ref{DG_matrix2}), respectively. In particular we use the Fortran subroutines of Jones and Plassmann \cite{JonesPlassmannIC}. The mesh, time stepping and other parameters are the same as in the two dimensional example above. The degree of the approximation space of $u$ and $v$ is chosen to be $p = 3$.

\input{manuscript_table/number_of_iterations.tex}

\graphicspath{{graphs_sofar/}}
\begin{figure}[htb!]
        \begin{center}
                \includegraphics[width=0.9\textwidth]{sofar0_25}
                \includegraphics[width=0.90\textwidth]{sofar0_50}
                \includegraphics[width=0.90\textwidth]{sofar0_75}
                \includegraphics[width=0.90\textwidth]{sofar1_00}
                \includegraphics[width=0.90\textwidth]{sofar1_25}
        \end{center}
        \caption{Time snapshots at times 0.25s, 0.5s, 0.75s, 1.0s  and 1.25s. As the waves progress from left to right the guiding due to the variable speed of sound becomes increasingly pronounced.  The black inlay is a sketch of the sound speed profile. \label{fig:sofar}}
\end{figure}

In Table \ref{iteration_number}, we show the $L^2$ errors in $u$, the corresponding convergence rate and the average number of the iterations for solving systems of $\frac{dU}{dt}$ and $\frac{dV}{dt}$. We observe the $4$-th order convergence in the $L^2$ norm for $u$ which is the same with results for the constant variable $c^2$. In addition, the average number of iterations is comparably small for all different degrees of freedom $(N+1)^2$ on each DG element. The average number of the iterations for $\frac{dU}{dt}$ is less than $6$ and the average number of the iterations for $\frac{dV}{dt}$ is always $1$. From these results (which are largely representative for many experiments that we have conducted), we see that the iterative methods for solving the scheme (\ref{DG_matrix1})--(\ref{DG_matrix2}) with variable stiffness matrix or variable mass matrix work well and are not cost-prohibitive. However the observed increase in the number of iterations with the number of degrees-of-freedom per element, which is
expected for the sort of preconditioner we are using, suggests that this number should be held fixed and the number of DG elements increased as we
refine the mesh.

\subsection{Propagation of waves in an ocean channel}\label{nice_phenomenon}
As a final and more applied example we consider the propagation of sound from a point source 100 meters below the surface in the ocean. The point source has a $\sin (400\pi t)$ sinusodial time dependence.  The speed of sound in water is taken to be depth dependent with the formula for the speed of sound in water being
$$
c = 1450+50 \left( \frac{100}{D+40} + \tanh \left(\frac{D-300}{50}\right) \right).
$$
Here $D$ is the depth (with a positive value) below the surface. As the sound speed profile has an inflection point the waves will be guided as they propagate in the direction parallel to the surface. Depending on the ``strength'' of the guiding the effect will become visible at a few or many wavelengths. The profile we use here is relatively weak and thus it is important to have a numerical method that is able to accurately propagate waves with minimal error over long distances.

Here the computational domain is taken to be $(x,y) \in [0m, 4000m] \times [-2000m,0m]$ and is discretized by $100 \times 50$ elements each with $41 \times 41$ grid points. We take $p=3$ and use the same Runge-Kutta method as above. In Figure \ref{fig:sofar} we display the solution $u$ at different snapshots in time. As the wave fronts evolve to the right in the domain they concentrate in side the minima of the sound speed profile illustrating that the method is able to capture these phenomena well.

We note that for the purely depth-dependent sound speed profile considered here the diagonalization method used for constant coefficient problems
can be used. However, to illustrate its utility for a more complex problem we used the iterative method in this case also.


\section{Introduction}
	Since the classic review by Kogut and Wilson \cite{Wilson:1973jj} on the $\epsilon$ expansion and renormalization group (RG) flow, the general properties of RG flows have been the subject of active research. In the cases usually considered, once a theory starts flowing, it ends up at a fixed point where it is described by some conformal field theory (CFT).
	From a general point of view, the equations describing instances of RG flow form systems of autonomous differential equations, and the properties of such systems and the kinds of flows they admit are well understood \cite{guckenheimer2013nonlinear,arnold2012geometrical,Gukov:2016tnp,feigenbaum1978quantitative}. In particular, dynamical systems can exhibit flows more peculiar than that between distinct fixed points, and Kogut and Wilson speculated in 1974 on the possibility of limit cycles as well as ergodic and turbulent behaviour in RG flow.
	Since then, however, a number of monotonicity theorems have been proven that severely restrict the RG flow of unitary quantum field theories (QFTs). The first such theorem was Zamolodchikov's $c$-theorem \cite{zamolodchikov1986irreversibility}, which in two dimensions establishes a function that interpolates between central charges at CFTs and decreases monotonically along RG flow. Analogous theorems were proven in four dimensions ($a$-theorem) \cite{Komargodski:2011vj,Luty:2012ww} and  three dimensions  ($F$-theorem) \cite{Klebanov:2011gs,Jafferis:2011zi,Casini:2012ei}. The monotonicity implied by these theorems excludes the possibility of limit cycles, except for a loophole pointed out in \cite{morozov2003can,Curtright:2011qg}: multi-valued $c$ functions. This loophole had in fact been previously realized in certain deformed Wess-Zumino-Witten models \cite{Bernard:2001sc,LeClair:2003hj,Leclair:2003xj}, although these models required coupling constants to pass between infinity and minus infinity in order to realize cyclic RG flow. There are also examples of cyclics RG flow in quantum mechanics \cite{Glazek:1993qs,Glazek:2002hq,Bulycheva:2014twa,Gorsky:2013yba,LeClair:2002ux,Braaten:2004pg,Dawid:2017ahd}.

	Recently, ref. \cite{Jepsen:2020czw} put forward a QFT of interacting symmetric traceless matrices transforming under the action of the $O(N)$ group, while allowing $N$ to assume non-integer values. $O(N)$ models for non-integer $N$, an idea widely used in polymer physics \cite{de1979scaling}, had been previously given a formal definition in \cite{Binder:2019zqc}, which demonstrated the non-unitarity of these models.  Hence, the $c,a,F$-theorems are no longer valid and do not constrain the RG flow, and consequently  ref. \cite{Jepsen:2020czw} was able to show that the model studied therein possesses a closed limit cycle for $N$ slightly above $N_*\approx 4.475$. The main tool used to make this discovery was Hopf's theorem \cite{hopf1942bifurcation}, which guarantees the existence of a limit cycle in the vicinity of the codimension-one bifurcation known as the Andronov-Hopf bifurcation.

	Turning to dynamical systems parameterized by two real numbers, codimension-two bifurcations can be used to prove the occurrence of yet other kinds of flow. Specifically, R.Bogdanov \cite{bogdanov} and F.Takens \cite{takens2001forced} have established powerful theorems by which, from properties of autonomous differential equations known only to second order in the dynamical variables, one can deduce the existence of homoclinic orbits, i.e. flow curves that connect a fixed point to itself. In addition to mild genericity conditions, the conditions that must be satisfied in order for the theorems to apply can be checked merely by studying the stability of fixed points, despite the fact that homoclinic orbits signal global bifurcations \cite{guckenheimer2013nonlinear} since they arise when a limit cycle collides with a saddle point.

	An interesting fact about homoclinic orbits is that they can be used to diagnose chaos. In applications of the theory of dynamical systems to physics, chaotic behavior \cite{cvitanovic2017universality} occurs in many instances, such as in turbulence \cite{ruelle1971nature,zakharov2012kolmogorov}, meteorology \cite{lorenz1963deterministic}  and even in scattering amplitudes in string theory \cite{Gross:2021gsj}. Usually, chaotic behaviour is proven via numerical investigations of concrete systems. One of the few analytical tools that can hint at the emergence of chaos is a theorem due to Shilnikov \cite{shilnikov1965case} that, for systems possessing homoclinic orbits, stipulates conditions by which to show they are chaotic. Therefore, one important step towards uncovering chaotic RG flow is to establish the existence of homoclinic RG flow.

	Brief previous mention of homoclinic RG flow can be found in \cite{Oliynyk:2005ak,Kuipers:2018lux}, which study non-linear sigma models and $QCD_4$ in the Veneziano limit. These references, however, mention the phenomenon solely for the purpose of pointing out its impossibility in those contexts.

	In this short letter, we study a QFT with global $O(N)\times O(M)$ symmetry. Examining the RG flow of the theory as a function of $M$ and $N$, we determine the regime where the flow is non-monotonic. In this regime, we are able to establish the locations of a number of Bogdanov-Takens bifurcations, by which we are able to conclude that the theory exhibits homoclinic RG flow. In other words, the model contains fixed point with the peculiar property that a deformation by a {\it relevant} operator induces a flow that leads back to the original point: an RG flow where the IR and UV theories are one and the same. Homoclinic RG flow can be thought of as interpolating between the familiar type of RG flow (where a system flows from one fixed point to another) and the more exotic RG limit cycles (like limit cycles, homoclinic orbits are closed
	. In unitary QFTs, homoclinic RG flows are still forbidden by $c,a,F$-theorems, but a fixed point situated in a homoclinic orbit could possibly be described by a standard CFT, in contrast to fixed points that give rise to limit cycles by undergoing a Hopf bifurcation, and which require operators with complex scaling dimensions.

	The method we adopt can be applied more broadly to find homoclinic orbits in two-parametric families of theories. We expect the phenomenon to be present in many other QFTs.

	\section{The model}
	\label{sec:model}

	We consider an $\mathcal{N}=1$ supersymmetric model of interacting scalar superfields $\Phi^i_{ab}$ that is invariant under the action of an $O(N)\times O(M)$ group in $d=3-\epsilon$ dimensions. The superfields are traceless-symmetric matrices with respect to the action of an $O(N)$ group and vectors under the action of an $O(M)$ group. There are four singlet marginal operators
	\begin{gather}
		O_1 = \operatorname{tr}\left[ \Phi^i \Phi^i \Phi^j \Phi^j\right], \quad O_2= \operatorname{tr} \left[ \Phi^i \Phi^j \Phi^i \Phi^j\right], \notag\\
		O_3 = \operatorname{tr} \left[\Phi^i\Phi^i\right]^2,\quad O_4 = \operatorname{tr} \left[\Phi^i \Phi^j\right] \operatorname{tr} \left[\Phi^i \Phi^j\right]\,, \label{eq:Ops}
	\end{gather}
	and so the full action is
	\begin{gather}
		S = \int d^d x d^2 \theta\Big( \operatorname{tr} \left[\Phi^i D^2_\alpha \Phi^i\right] + \sum_i g_i O_i \Big)\,.
		\label{action}
	\end{gather}

	The RG flow of this model is gradient, meaning that there exists a function $F$ of the couplings and a four-by-four matrix $G_{ij}$ such that the beta functions of the theory satisfy the equation
	\begin{gather}
		\beta_i = \mu\frac{d g_i}{d \mu} = G_{ij} \frac{\partial F}{\partial g_j}\,.
	\end{gather}
	If $G_{ij}$ is positive or negative definite, this equation implies that $F$ changes monotonically with the RG flow, so that cyclic and homoclinic flow lines are impossible. By explicit computation to leading order in perturbation theory, we find that the metric has determinant
	\begin{gather}
		\det G  = \frac{1}{4}(M-1)^2(M+2)^2(N-3) \times\notag\\
		\times(N-2)^2(N+1)^2(N+4)^2(N+6)\,. \label{eq:det}
	\end{gather}
	We list the beta functions and the components of the metric in appendix \ref{sec:betas}. The zeroes in $\det G$ occur because of linear relations among the four operator of the theory at special values of $M$ and $N$, and their presence indicates that eigenvalues change sign as $N$ and $M$ are varied.
	Indeed one can check that the metric is sign-indefinite if $M\in\left(-2,1\right)$ or $N\in \left(-6,3\right)$, so that unusual RG flows are possible in this regime, and operators may develop complex scaling dimensions at real fixed points, which, in the terminology of \cite{Jepsen:2020czw}, are then termed "spooky". At integer values of $N$ and $M$, such operators are identically zero owing to the linear relations between the operators. The situation is closely analogous to the occurrence of evanescent operators at non-integer spacetime dimensions \cite{Collins:1984xc,Bos:1987fb,Dugan:1990df,Gracey:2008mf,Hogervorst:2014rta,Hogervorst:2015akt}.

	In the following, we will allow $M$ and $N$ to assume general real values. In consequence of this analytic extension, we are able to observe Hopf bifurcations taking place in the model along various curves in the space of $M$ and $N$. See fig. \eqref{fig:cross}. But while Hopf bifurcations are a type of codimension-one bifurcation widely found in one-parameter systems of autonomous systems of differential equations, we are dealing with a two-parameter system, and such systems are capable of exhibiting a richer variety of flows. The possible codimension-two bifurcations can be classed into five types \cite{arnold2012geometrical,guckenheimer2013nonlinear} -- Bautin, Bogdanov-Takens, cusp, double-Hopf, and zero-Hopf -- which signal different kinds of flow not present in generic one-parameter systems. As we shall now see, some of these possibilities are realized by the QFT with action \eqref{action}.

	\begin{figure}
		\begin{center}
			\includegraphics[width=0.45\textwidth]{crossPlot.png}
		\end{center}
		\caption{ Bifurcation diagram in the space of $M$ and $N$. Spooky fixed points are present in the areas coloured in \textcolor{colour1}{blue}. Curves in \textcolor{colour5}{yellow} signify Andronov-Hopf bifurcations. Points coloured in \textcolor{red}{red} mark the locations of Bogdanov-Takens bifurcations, and points in \textcolor{green}{green} indicate Zero-Hopf bifurcations.}
		\label{fig:cross}
	\end{figure}

	\section{Bogdanov-Takens Bifurcation}
	\label{sec:BT}
	A Bogdanov-Takens bifurcation occurs generically when, at a fixed point, two eigenvalues of the stability matrix $\left(\frac{\partial \beta_i}{\partial g_j}\right)$ tend to zero as two bifurcation parameters $M$ and $N$ are appropriately tuned.

	The following equations must then be satisfied:
	\begin{gather}
		\beta_i(g_i,N,M) = 0\,, \quad \det \left(\frac{\partial \beta_i}{\partial g_j}\right) \left(g_i,N,M\right)=0\,, \label{eq:BTcon}\\
		\operatorname{tr}\left[\textstyle\bigwedge^3\left( \frac{\partial \beta_i}{\partial g_j}\right)\right] \equiv \det \left(\frac{\partial \beta_i}{\partial g_j}\right) \operatorname{tr}\left[ \left(\frac{\partial \beta_i}{\partial g_j}\right)^{-1} \right] = 0\,. \nonumber

	\end{gather}
	Written in the form \eqref{eq:BTcon}, we see that the conditions for a BT bifurcation are polynomial equations in $g_i$, $M$, and $N$, and so by B\'ezout's theorem there exist at most a finite number of points that satisfy these conditions. We refer to such points as Bogdanov-Takens (BT) points. For the QFT we are studying perturbatively, it can be verified that the beta functions exhibit several such points, as shown in figure \ref{fig:cross}.
	Their existence can be checked to high numerical accuracy with the use of standard programs, e.g. PyDSTool \cite{pydstool}.
	Higher-loop contributions will provide corrections to the precise locations of these points, but as long as we take $\epsilon$ to be sufficiently small,  higher-order corrections will not alter the number or qualitative behaviour of BT points.

	While two eigenvalues tend to zero as we approach a BT point, right at the BT point itself we do not typically have a pair of eigenvectors with zero eigenvalues, the reason being that in this same limit, the two respective eigenvectors usually become linearly dependent. Rather, the stability matrix at a BT point has a Jordan block of size two with zero eigenvalue (see \eqref{eq:numstabmatrix} in appendix \ref{sec:deriv}). This means that the theory at the BT point possesses two operators $\mathcal{O}_{1,2}$ such that the generator $D$ of dilatations acts in the following way
	\begin{gather}
		D \mathcal{O}_1 = d \mathcal{O}_1\,,\quad D\mathcal{O}_2 =d\mathcal{O}_2+\mathcal{O}_1\,.
	\end{gather}
	The possibility of indecomposable representations of the conformal group was extensively studied in \cite{Gurarie:1993xq,Hogervorst:2016itc}. The upshot is that the BT theory constitutes a logarithmic CFT containing {\it generalized} marginal operators $\mathcal{O}_{1,2}$. In consequence, BT theories are non-unitary and we have
	\begin{gather*}
		\braket{\mathcal{O}_2(0) \mathcal{O}_2(x)} =  -\frac{2 k_\mathcal{O} \log \left|x\right|}{\left|x\right|^{2d}}\,,\quad   \braket{\mathcal{O}_1(0) \mathcal{O}_2(x)} =  \frac{ k_\mathcal{O}}{\left|x\right|^{2d}}\,,
	\end{gather*}
	for some constant $k_\mathcal{O}$.
	\begin{figure}
		\begin{center}
			\includegraphics[width=0.4\textwidth]{chart.pdf}
		\end{center}
		\caption{Bifurcation diagram around the Bogdanov-Takens bifurcation at $(M,N)=(M^\ast,N^\ast)\approx (0.2945,4.036)$. $\delta M = M-M^\ast$, $\delta N=N-N^\ast$. The \textcolor{blue}{blue} curve represents a saddle node bifurcation, the \textcolor{green}{green} curve represents a Hopf bifurcation, and the \textcolor{red}{red} curve represents a saddle homoclinic bifurcation. At the origin, these three codimension-one bifurcations coalesce.}
		\label{fig:bifurcationDiagram}
	\end{figure}

	\begin{figure*}

		\begin{subfigure}{0.45\textwidth}
			\begin{center}
				\includegraphics[width=0.7\textwidth]{region1.pdf}
			\end{center}
			\caption*{RG flow in region \raisebox{.5pt}{\textcircled{\raisebox{-.9pt} {1}}}. In this regime, there are no fixed points near the origin, and all flow curves swerve downwards in the IR.}
		\end{subfigure}
		\begin{subfigure}{0.45\textwidth}
			\begin{center}
				\includegraphics[width=0.7\textwidth]{region2.pdf}
			\end{center}
			\caption*{RG flow in region \raisebox{.5pt}{\textcircled{\raisebox{-.9pt} {2}}}. The fixed point marked in \textcolor{red}{red} is a saddle point. The \textcolor{colour1}{blue} curves flow outwards from this point in the IR, while the \textcolor{colour2}{orange} curves flow inward. The fixed point marked in \textcolor{colour4}{green} is IR unstable.  }
		\end{subfigure}
		\begin{subfigure}{0.45\textwidth}
			\begin{center}
				\includegraphics[width=0.7\textwidth]{region3.pdf}
			\end{center}
			\caption*{RG flow in region  \raisebox{.5pt}{\textcircled{\raisebox{-.9pt} {3}}}. In passing from region \raisebox{.5pt}{\textcircled{\raisebox{-.9pt} {2}}} to \raisebox{.5pt}{\textcircled{\raisebox{-.9pt} {3}}}, the fixed point marked in \textcolor{colour4}{green} has undergone a Hopf bifurcation and is now IR stable. An IR-repulsive limit cycle, marked in \textcolor{colour3}{cyan}, separates the two fixed points.}
		\end{subfigure}
		\begin{subfigure}{0.45\textwidth}
			\begin{center}
				\includegraphics[width=0.7\textwidth]{region4.pdf}
			\end{center}
			\caption*{RG flow in region \raisebox{.5pt}{\textcircled{\raisebox{-.9pt} {4}}}. In passing from region \raisebox{.5pt}{\textcircled{\raisebox{-.9pt} {3}}} to \raisebox{.5pt}{\textcircled{\raisebox{-.9pt} {4}}} the limit cycle collided with the \textcolor{red}{red} fixed point in a homoclinic bifurcation. }
		\end{subfigure}
		\caption{The topologically distinct types of RG flow in the vicinity of the Bogdanov-Takens bifurcation at $M=M^\ast \approx 0.2945$ and $N=N^\ast \approx 4.036$. The variables $y_3$ and $y_4$ are linear combinations of the four coupling constants $g_i$, with precise definitions given in appendix \eqref{sec:deriv}, and $\delta M = M - M^\ast$, $\delta N = N - N^\ast$.}
		\label{fig:4regimes}
	\end{figure*}
	The conditions \eqref{eq:BTcon} are not entirely sufficient to guarantee a BT bifurcation. One must also require smoothness and a set of inequalities that are generically true. Violations of the inequalities typically require fine-tuning of additional parameters and signal bifurcations of codimension higher than two. Incidentally, at the integer values $M=2$ and $N=3$, right on the boundary of the regimes with monotonic and non-monotonic RG flows, we observe a fixed point that satisfies \eqref{eq:BTcon}, but which fails to meet these genericity requirements and for this reason is not described by a logarithmic CFT.

	In appendix \ref{sec:deriv} we give the precise statement of the Bogdanov-Takens bifurcation theorem, and we explicitly check that it applies to an example of a BT point in the QFT we are studying, situated at  $M \approx 0.2945$ and $N \approx 4.036$. What this means is that we can transform the beta functions near the BT point into a particularly simple form, known as Bogdanov normal form:
	\begin{gather}
		\left\{\begin{matrix*}[l]
			\dot{\eta}_1=\eta_2\,, \\
			\dot{\eta}_2 = \delta_1+\delta_2 \eta_1+\eta_1^2 + s \eta_1 \eta_2 + \mathcal{O}\left(\left|\eta \right|^3\right)\,,\\
			\dot{\eta}_i=\lambda_i \eta_i\quad \text{for }i>2\,,
		\end{matrix*}\right. \label{eq:BTsys}
	\end{gather}
	where $s=- 1$, and $\delta_{1,2}$ are functions of $N$ and $M$ that vanish right at the BT point.

	By bringing the system into normal form, we can use the equations \eqref{eq:BTsys} to determine the behaviour of the system for small enough $\delta_1$ and $\delta_2$. In particular, we can constrain ourselves to studying the surface where only $\eta_1$ and $\eta_2$ are non-zero, noting that the dynamics in the transverse directions $\eta_3$ and $\eta_4$ are quite simple. Depending on the values of $\delta_1$ and $\delta_2$, the flow of $\eta_{1,2}$ falls into different topological types. The classification can be found in textbook \cite{kuznetsov2013elements} and amounts to the following. In the vicinity of the BT point at $\delta_1=\delta_2=0$, there are four regimes with qualitatively different flows:

	-- Regime \raisebox{.5pt}{\textcircled{\raisebox{-.9pt} {1}}}: The flow has no fixed point.

	In the other three regimes, the flow has two fixed-points, which we will label left and right. The right point is always a saddle point.

	-- Region \raisebox{.5pt}{\textcircled{\raisebox{-.9pt} {2}}}: The left point is unstable, and all flow lines starting near it terminate at the right fixed point.

	-- Region \raisebox{.5pt}{\textcircled{\raisebox{-.9pt} {3}}}: The left point is now stable, and a repulsive limit cycle separates the two fixed points.

	-- Region \raisebox{.5pt}{\textcircled{\raisebox{-.9pt} {4}}}: The left point is still stable, but the limit cycle has disappeared. Some flow lines starting near the right fixed point terminate at the left fixed point.

	In the case of the BT point at $(M,N) \approx (0.2945,4.036)$, the locations of these four adjoining regimes, as computed in appendix \eqref{eq:BTsys}, is shown in figure \ref{fig:bifurcationDiagram}. And the RG flow in each regime is depicted in figure \ref{fig:4regimes}.

	The four regimes are separated by different codimension-one bifurcations. Region \raisebox{.5pt}{\textcircled{\raisebox{-.9pt} {1}}} is demarcated from regions \raisebox{.5pt}{\textcircled{\raisebox{-.9pt} {2}}} and \raisebox{.5pt}{\textcircled{\raisebox{-.9pt} {4}}} by a saddle-node bifurcation happening at $\delta_1=\frac{1}{4}\delta_2^2$. Regions \raisebox{.5pt}{\textcircled{\raisebox{-.9pt} {2}}} and \raisebox{.5pt}{\textcircled{\raisebox{-.9pt} {3}}} are separated by an Andronov-Hopf bifurcation along the the half-curve $\delta_1=0$, $\delta_2< 0$. And regions \raisebox{.5pt}{\textcircled{\raisebox{-.9pt} {3}}} and \raisebox{.5pt}{\textcircled{\raisebox{-.9pt} {4}}} are separated by a saddle homoclinic bifurcation along $\delta_1=-\frac{6}{25}\delta_2^2+\ldots, \, \delta_2 <0$.

	\begin{figure}
		\centering
		\includegraphics[scale=0.6]{Homo.pdf}
		\caption{Flow diagram for a dynamical system containing a homoclinic orbit (marked in \textbf{black}), ie. a flow line that starts and ends at the same point. The system is described by equations \eqref{eq:A10} with parameters $\delta_1=-0.000453178$ and $\delta_2=-0.0440214$. The \textcolor{red}{red} and \textcolor{green}{green} dots indicate fixed points. The \textcolor{green}{green} dot is a "spooky" fixed point. The theory at the \textcolor{red}{red} dot is a homoclinic CFT.}
		\label{fig:homoclorbit}
	\end{figure}

	A saddle-node bifurcation corresponds to the collision and disappearance of two equilibria in dynamical systems. The phenomenon has been observed in a number of cases of RG flow, it happens for instance in in the critical $O(N)$ model \cite{fei2015three} and in prismatic models \cite{Giombi:2018qgp}, and has been proposed to occur in $QCD_4$ \cite{Gies:2005as,kaplan2009conformality,Gorbenko:2018ncu,Kuipers:2018lux}.

	An Andronov-Hopf bifurcation represents a change of stability at a fixed point that has complex eigenvalues. The flow near the fixed point changes between spiraling inwards and spiraling outwards and gives birth to a limit cycle. In the context of RG, this bifurcation was recently studied in \cite{Jepsen:2020czw}.

	The most interesting and new phenomenon associated to the model of the present paper happens along the homoclinic bifurcation line. Here the flow exhibits what is known as a homoclinic orbit.

	\section{Homoclinic RG flow}

	A homoclinic orbit is a flow line that connects a stable and an unstable direction of a saddle point. Figure \eqref{fig:homoclorbit} depicts the kind of homoclinic orbit generated by a BT bifurcation, with the saddle point marked by a \textcolor{red}{red} dot. The homoclinic orbit is seen to envelop another fixed point marked in \textcolor{green}{green}. In a QFT context, the \textcolor{green}{green} point is "spooky": the couplings are real, but the eigenvalues of the stability matrix $\left(\frac{\partial \beta_i}{\partial g_j}\right)$ have non-zero imaginary parts. In contrast to such spooky points, and to complex CFTs \cite{Gorbenko:2018ncu,Gorbenko:2018dtm}, the
	\textcolor{red}{red} saddle point is associated to real couplings and real eigenvalues of the stability matrix. These eigenvalues are small and have opposite signs: $\lambda_1,-\lambda_2 \ll 1$. The positive eigenvalue corresponds to a slightly relevant operator $\mathcal{O}_1$ with dimension $\Delta_1=d+\lambda_1 > d$, and the negative eigenvalue to a slightly irrelevant operator $\mathcal{O}_2$ with dimension $\Delta_2=d + \lambda_2 < d$. In this sense, the \textcolor{red}{red} saddle point corresponds to a real CFT.

	Standard RG lore states that if we perturb a system in the direction of a relevant operator, then we expect for the system to either lose conformality altogether or to flow to a different CFT. In the terminology of dynamical systems, standard RG trajectories are heteroclinic orbits. The classical example is the Wilson-Fischer fixed point: by carefully perturbing a Gaussian theory in $4-\epsilon$ dimension we flow to a weakly coupled interacting CFT, which in three dimensions interpolates to the Ising model. Homoclinic bifurcations provide exotic counterexamples to this general picture: if we perturb the system in the direction of a relevant operator, we come back to the original fixed point, which tentatively we can term a {\it homoclinic CFT}. Such RG behaviour obviously violates the $F$-theorem so that homoclinic fixed points must be non-unitary, as is generally the case for CFTs with symmetry groups of non-integer rank \cite{Binder:2019zqc}.

	If we tune the bifurcation parameters so as to approach the BT point along the saddle homoclinic bifurcation (the red curve in figure \ref{fig:bifurcationDiagram}), then the homoclinic orbit shrinks to a point and vanishes. In this limit, the \textcolor{red}{red} homoclinic CFT and the \textcolor{green}{green} spooky fixed point merge and become a logarithmic CFT.

	\section{Zero-Hopf Bifurcations: The Road to Chaos}

	The Bogdanov-Takens bifurcation is not the only codimension-two bifurcation that can be observed to take place in the model \eqref{action}. The theory also possesses two points in the space of $g_i$, $M$, and $N$ where the stability matrix has a pair of purely imaginary eigenvalues and one zero eigenvalue. Such fixed points indicate what is known as a Zero-Hopf (ZH) or a Fold-Hopf bifurcation. This type of bifurcation was classified in \cite{takens1974singularities} and can be divided into six sub-types. In the notation of \cite{guckenheimer2013nonlinear}, the model has a type I ZH bifurcation at $(M,N)\approx (0.8447,-1.807)$ and a type IIa ZH bifurcation at $(M,N)\approx (-3.816,1.188)$. At a type I bifurcation point, a saddle-node bifurcation is incident to a pitchfork bifurcation, and there are no nearby cyclic orbits. At a type IIa point, a saddle-node bifurcation is again incident to a pitchfork bifurcation, but additionally a Hopf bifurcation is also incident to the point, except that the stability coefficient of the associated limit cycle (what was referred to as the Hopf constant in \cite{Jepsen:2020czw}) exactly vanishes in a quadratic approximation, so that cubic fluctuations or higher decide the fate of the cyclic flow near a type IIa point.

	Generally, ZH bifurcation points are of particular interest because it is known that in their vicinity what is known as a Shilnikov homoclinic orbit may develop and render the system chaotic \cite{kuznetsov2013elements,shilnikov1965case}. Recently it was proven in \cite{baldoma2020hopf} that the presence of ZH bifurcations of type III guarantees the existence of a Shilnikov orbit and a nearby infinite set of saddle periodic orbits. This nontrivial invariant set can be embedded in an attracting domain, thus implying Shilnikov chaos.

	The ZH points of the model in the present paper are not of type III, and we cannot claim that the system is chaotic. It may be worthwhile to investigate if there exist other models that meet the simple criteria for the assured appearance of chaos.

	\section{Conclusion and Outlook}
	The approach suggested and adopted in \cite{Gukov:2016tnp,Kuipers:2018lux,Jepsen:2020czw} of studying the beta functions and renormalization of QFTs from the general perspective of dynamical systems provides a method of understanding the full range of possible RG flows. A powerful tool to this end is offered by Bogdanov's and Taken's bifurcation theorem \cite{hopf1942bifurcation}, which lists a simple set of conditions that guarantee the existence of a homoclinc RG orbit, and which can be checked already at first order in perturbation theory.

	In this short letter we have presented a QFT that satisfies these conditions, namely a supersymmetric model with global symmetry group $O(N) \times O(M)$, where $N$ and $M$ play the role of the bifurcation parameters of the system. We determined a number of parameter values where a BT bifurcation takes place and investigated the nearby RG flow to uncover the presence of homoclinic orbits, where the perturbation of a fixed point by a relevant operator induces an RG flow that returns to its starting point along an irrelevant direction.

	There are several bifurcation theorems that give simple criteria for other novel kinds of RG flows \cite{guckenheimer2013nonlinear,kuznetsov2013elements,arnold2012geometrical}. Some of these theorems allow for the determination of the onset of chaotic flow based on straightforward computations around fixed points \cite{baldoma2020hopf}. It would be interesting to find out if QFTs give birth to chaos when $N$ becomes fractional.

	\section*{Acknowledgments}
	We are grateful to Igor R. Klebanov for very insightful discussions and suggestions throughout the project. We are also grateful to Alexander Gorsky, Alexei Milekhin, Yuri Kuznetsov, Alexander Polyakov, Sergey Gukov, Slava Rychkov, and Bernardo Zan for valuable discussions and comments. We also thank Maikel Bosschaert for spotting a considerable typo in appendix B.
	This research was supported in part by the US NSF under Grant No. PHY-1914860.