Further author information: (Send correspondence to Harshit Agrawal.)
E-mail<EMAIL_ADDRESS>
Accepted for publication in Proc. SPIE, Medical Imaging 2024: Physics of
Medical Imaging

# Utilizing U-Net Architectures with Auxiliary Information for Scatter
Correction in CBCT Across Different Field-of-View Settings

Harshit Agrawal  Aalto University, Espoo, Finand.  Planmeca Oy., Helsinki,
Finland.  Ari Hietanen  Planmeca Oy., Helsinki, Finland.  Simo Särkkä  Aalto
University, Espoo, Finand.

###### Abstract

Cone-beam computed tomography (CBCT) has become a vital imaging technique in
various medical fields but scatter artifacts are a major limitation in CBCT
scanning. This challenge is exacerbated by the use of large flat panel 2D
detectors. The scatter-to-primary ratio increases significantly with the
increase in the size of FOV being scanned. Several deep learning methods,
particularly U-Net architectures, have shown promising capabilities in
estimating the scatter directly from the CBCT projections. However, the
influence of varying FOV sizes on these deep learning models remains
unexplored. Having a single neural network for the scatter estimation of
varying FOV projections can be of significant importance towards real clinical
applications. This study aims to train and evaluate the performance of a U-Net
network on a simulated dataset with varying FOV sizes. We further propose a
new method (Aux-Net) by providing auxiliary information, such as FOV size, to
the U-Net encoder. We validate our method on 30 different FOV sizes and
compare it with the U-Net. Our study demonstrates that providing auxiliary
information to the network enhances the generalization capability of the
U-Net. Our findings suggest that this novel approach outperforms the baseline
U-Net, offering a significant step towards practical application in real
clinical settings where CBCT systems are employed to scan a wide range of
FOVs.

###### keywords:

CBCT, Scatter Correction, U-Net, Auxiliary Information, FOV.

##  1 INTRODUCTION

One major drawback of cone-beam computed tomography (CBCT) systems is the
increase in scatter-to-primary signal ratio due to the use of large flat panel
2D detector and the proximity of detector to the patient. Several deep
learning-based methods have been proposed in recent years for the scatter
estimation. A U-Net architecture was proposed for low-frequency scatter
estimation for different tube voltages, anatomies, and noise levels [1]. This
architecture estimated scatter signal directly from the scatter-affected
projections. Similarly, other versions of U-Net and its variants have been
investigated for the scatter estimation [2, 3, 4].

In CBCT systems, the scatter may vary significantly with the size of scanning
field-of-view (FOV). The scatter-to-primary ratio increases 4.9-fold for a
17x12 cm FOV in comparison to a 6x6 cm FOV [5]. In modern CBCT systems, the
FOV size can be even larger, which will increase the scatter-to-primary ratio
even more. While the applicability of deep learning networks has been
demonstrated for estimating the scatter for a fixed FOV, no study has been
done to evaluate deep learning network’s performance with changes in FOV size.
Moreover, it remains unclear if a single neural network can generalize to
different FOVs. In practice, having a single network for multiple FOVs is
important as CBCT systems are employed to scan a wide range of FOVs in real
clinical settings.

In this study, we train and evaluate a single U-Net for the scatter estimation
performance [1] on varying FOV sizes. Further, we propose a simple method to
improve the performance of the single network for multiple FOV sizes by
providing auxiliary FOV size information to the encoder (Aux-Net), which
outperforms the baseline U-Net.

##  2 Materials and Methods

###  2.1 Dataset

We utilized two different sources to generate the training dataset. We used 3
CT reconstructions from HNSCC-3DCT-RT dataset [6, 7] for generating training
data. Further, 3 reconstructions were acquired by scanning 3 anthropomorphic
phantoms using GE CT scanner, out of which we used 2 CT reconstructions for
training and 1 for testing. We used Monte Carlo based X-ray simulations [8, 9]
to generate primary (without scatter) and the corresponding scatter
projections. For each reconstruction in training, we simulated total 18
different sizes of FOVs. For test phantom, we simulated a total of 30 FOV
sizes which were different than the FOV sizes used in the training. The FOVs
used in the training and the testing are given in the table 1. We simulated
the geometry of Planmeca Viso G7® CBCT scanner (Planmeca Oy, Helsinki,
Finland) for the jaw protocol. The source-to-detector distance was 700 mm and
the source-to-object distance was 490 mm. For each phantom in the training,
the projections were simulated from 0 to 210 degrees with an spacing of 2.1
degrees. We simulated 100 projections for each phantom. We further augmented
our dataset by simulating each projection 10 times with a different random
seed to generate random direction and energy of X-ray. We averaged these
simulated projections with each other to obtain 10 different noise levels for
training. For example, the average projection of two simulated projections has
different noise level than the average of 10 simulated projections. For test
dataset, we simulated 500 views for each phantom ranging from 0 to 210
degrees. We also scanned an anthropological phantom using Planmeca Viso G7®
CBCT scanner (Planmeca Oy, Helsinki, Finland) for $170\times 170$ mm FOV to
compare the results of the scatter correction methods qualitatively.

Table 1: The FOV sizes used for the training (X) and the testing (✓). diameter (mm) | height (mm)
---|---
| 30 | 40 | 50 | 60 | 70 | 80 | 90 | 100 | 110 | 120 | 130 | 140 | 150 | 160 | 170 | 180
120 | X |  |  | X |  |  | X |  |  | X |  |  | X |  |  | X
130 |  | ✓ | ✓ |  | ✓ | ✓ |  | ✓ | ✓ |  | ✓ | ✓ |  | ✓ | ✓ |
140 | X |  |  | X |  |  | X |  |  | X |  |  | X |  |  | X
150 |  | ✓ | ✓ |  | ✓ | ✓ |  | ✓ | ✓ |  | ✓ | ✓ |  | ✓ | ✓ |
160 | X |  |  | X |  |  | X |  |  | X |  |  | X |  |  | X
170 |  | ✓ | ✓ |  | ✓ | ✓ |  | ✓ | ✓ |  | ✓ | ✓ |  | ✓ | ✓ |

###  2.2 Network Architecture

The architecture of the proposed Aux-Net is shown in the Fig.1. To incorporate
auxiliary information about the FOV shape of the CBCT projections, we create
two additional channels, both having the same dimensions as the input
features. The first channel contains width information and the second channel
contains height information. The height and width are divided by the
respective detector maximum sizes, i.e. 300 mm and 250 mm, respectively. Since
magnification factor is same for all simulation geometries, the FOV shape is
directly related with the shape of the projection. These two channels are
concatenated to each level’s input features in the encoder. For each level of
the encoder, the channels are sized exactly to the input features to that
level. Aux-Net has a slightly increased 7.275697 million parameters compared
to the baseline U-Net’s 7.257552 million trainable parameters. The baseline
U-Net has same architecture as Aux-Net except the presence of auxiliary
channels.

Figure 1: Aux-Net architecture. $x$ in the first auxiliary channel is width
$w$ normalized by maximum width of detector (250 mm) and $y$ in the second
auxiliary channel is height $h$ normalized by maximum height of detector (300
mm). The auxiliary channels are recreated at each level of the encoder to
match the dimensions of the corresponding feature maps.

###  2.3 Training and Evaluation

For the training, the simulated primary and scatter projections were added and
divided by the simulated flats. The input to the model was further linearized.
All the input and target were resized to $300\times 250$. The model was
trained with a combination of mean-square-error (MSE) and high-frequency loss
to penalize the high-frequency component in the predicted scatter. For the
evaluation of neural network predictions, the mean-absolute-percentage-error
(MAPE) [1] and MSE were calculated between the upsampled neural network
scatter estimates and Monte Carlo simulated ground truth scatter.

##  3 Results

###  3.1 Scatter estimation on simulated data

We compare our method Aux-Net with baseline U-Net and Monte Carlo simulated
ground truth scatter projections. Fig. 2 shows MAPE and MSE errors between the
ground truth and the neural network-based scatter estimates. The three figures
show that proposed Aux-Net outperforms the baseline U-Net on all FOV sizes.

(a) MAPE and MSE errors for U-Net and AuxNet for FOV diameter of 130 mm.
X-axis shows the FOV heights in mm.

(b) MAPE and MSE errors for U-Net and AuxNet for FOV diameter of 150 mm.
X-axis shows the FOV heights in mm.

(c) MAPE and MSE errors for U-Net and AuxNet for FOV diameter of 170 mm.
X-axis shows the FOV heights in mm.

Figure 2: Comparison of the scatter estimates from U-Net and AuxNet with Monte
Carlo simulated ground truth scatter.

###  3.2 Scatter estimation on phantom scan

Figure 3 depicts the result of U-Net and Aux-Net based scatter corrected
reconstructions for the two axial slices at different locations. The predicted
scatter estimates were subtracted from the measurements and reconstruction was
done using FDK algorithm for CBCT. The intensity profile is more stable for
Aux-Net corrected image slices. Furthermore, the Aux-Net corrected
reconstruction is more robust in the presence of high density object in
comparison to the U-Net corrected reconstruction.

(a) Axial image slice number 113.

(b) Mean intensity (HU) profiles.

(c) Axial image slice number 221.

(d) Mean intensity (HU) profiles.

Figure 3: Comparison of the scatter corrected reconstructions from U-Net and
Aux-Net with the uncorrected reconstruction of an anthropomorphic phantom. The
intensity window of the axial images is -1000 to 2500 HU. Mean intensity
profiles of the rectangular area are plotted on the right.

##  4 Conclusion

In this study, we developed a new approach for improving deep learning based
scatter estimation across various FOV sizes. By concatenating the size
information with the input of the encoders in the U-Net, we increased the
robustness and generalization capability of the scatter estimation model. The
proposed Aux-Net outperformed the U-Net for all FOV sizes, specifically on
bigger FOVs. Moreover, the proposed method of using auxiliary information is
not confined solely to the FOV size. Including other variables that affect the
amount of scatter, such as, kV (kilovolt), patient size, object to detector
distance, scanned anatomy, and exposure time, as auxiliary information may
further enhance the robustness of the neural network and its clinical
applicability.

###### Acknowledgements.

The work was supported in-part by a Business Finland grant under TomoHead
project.

## References

  * [1] Maier, J., Eulig, E., Vöth, T., Knaup, M., Kuntz, J., Sawall, S., and Kachelrieß, M., “Real-time scatter estimation for medical CT using the deep scatter estimation: Method and robustness analysis with respect to different anatomies, dose levels, tube voltages, and data truncation,” Med. Phys. 46(1), 238–249 (2019).
  * [2] Nomura, Y., Xu, Q., Shirato, H., Shimizu, S., and Xing, L., “Projection-domain scatter correction for cone beam computed tomography using a residual convolutional neural network,” Med. Phys. 46(7), 3142–3155 (2019).
  * [3] Lee, H. and Lee, J., “A deep learning-based scatter correction of simulated x-ray images,” Electronics 8(9), 944 (2019).
  * [4] Roser, P., Birkhold, A., Preuhs, A., Syben, C., Felsner, L., Hoppe, E., Strobel, N., Kowarschik, M., Fahrig, R., and Maier, A., “X-ray scatter estimation using deep splines,” IEEE Trans. Med. Imaging 40(9), 2272–2283 (2021).
  * [5] Pauwels, R., Pittayapat, P., Sinpitaksakul, P., and Panmekiate, S., “Scatter-to-primary ratio in dentomaxillofacial cone-beam CT: effect of field of view and beam energy,” Dentomaxillofac. Radiol. 50(8), 20200597 (2021).
  * [6] Bejarano, T., Ornelas-Couto, M., and Mihaylov, I., “Head-and-neck squamous cell carcinoma patients with CT taken during pre-treatment, mid-treatment, and post-treatment (HNSCC-3DCT-RT) [Dataset].” Online (2018).  [doi:10.7937/K9/TCIA.2018.13upr2xf].
  * [7] Clark, K., Vendt, B., Smith, K., Freymann, J., Kirby, J., Koppel, P., Moore, S., Phillips, S., Maffitt, D., Pringle, M., et al., “The Cancer Imaging Archive (TCIA): Maintaining and Operating a Public Information Repository,” J. Digit. Imaging 26(6), 1045–1057 (2013).
  * [8] Badal, A. and Badano, A., “Accelerating Monte Carlo simulations of photon transport in a voxelized geometry using a massively parallel graphics processing unit,” Med. Phys. 36(11), 4878–4880 (2009).
  * [9] Agrawal, H., Hietanen, A., and Särkkä, S., “Deep Learning Based Projection Domain Metal Segmentation for Metal Artifact Reduction in Cone Beam Computed Tomography,” IEEE Access 11, 100371–100382 (2023).


Slope limiters play an essential role in maintaining the non-oscillatory behavior of high-resolution methods for nonlinear conservation laws. The family of minmod limiters serves as the prototype example. Here, we revisit the question of non-oscillatory behavior of high-resolution central schemes in terms of the  slope limiter proposed by van Albada et. al. 1982. The  van Albada (vA) limiter is smoother near extrema, and consequently,  in many cases, it outperforms  the results obtained using the standard minmod limiter. In particular, we prove that the vA limiter ensures 1D TVD stability and  demonstrate that it yields noticeable improvement in computation of one- and two-dimensional systems.

§ INTRODUCTION

We revisit the class of high-resolution, non-oscillatory schemes for the approximate solution of nonlinear conservation laws. This class of schemes
went through an intense period of development during the `80s and `90s.
By high-resolution we refer to the class of schemes which are formally  second- or higher-order methods at all computational cells, except perhaps  for finitely many critical cells, thus breaking the first-order accuracy barrier of Godunov for monotone schemes, [17]. By non-oscillatory we refer to  classes of schemes which satisfy certain type of stability bounds — weaker than monotonicity yet strong enough to exclude spurious oscillations. In this context, we mention the canonical example of  Total Variation Diminishing (TVD) schemes of Harten [14, 49, 5, 39, 40]. Other classes of non-oscillatory stable schemes followed, including finite-volume, finite-element, Discontinuous Galerkin (DG) and spectral schemes; as examples we mention the Total-Variation  (TVB) bounded schemes [43, 20, 52, 34, 57], the One-Sided Lip Condition  (OSLC) bounded E-schemes  [37, 50], and  entropy stable schemes [54, 8]. We also mention the notable class of Essentially Non-Oscillatory (ENO) and Weighted ENO schemes [15, 32, 45, 46] and their sign property stability bound [9, 8, 10].In this paper we focus on second-order Godunov-type methods. The key step of such  methods is the second-order accurate reconstruction. To fix our notations, we consider one-dimensional,  piece-wise constant scalar numerical solution at  time level $t^n := n\Dt$, which is realized in terms  of its cell-averages, $\{\baru_\alpha^n= \baru(x_\alpha,t^n)\}$, over equi-spaced mesh $\{x_{\alpha} :  x_\alpha= \alpha{\Dx}\}$,

\begin{equation*}
u^n = \sum^{}_{j} \baru^n_j {\mathds 1}_{[x_{\jmh}, x_{\jph}]}, \qquad  \baru(x_j,t^n) := \frac{1}{\Dx}\int^{x_{\jph}}_{x_{\jmh}}u(x,t^n) \dx, \qquad j\in {\mathbb Z}.
\end{equation*}

The key ingredient of Godunov-type schemes is a  piecewise linear  reconstruction, $u^{\Dx}(x,t^n) \approx u(x,t^n)$, which conserves the underlying cell averages, $\displaystyle \frac{1}{\Dx}\int^{x_{\jph}}_{x_{\jmh}}u^{\Dx}(x,t^n) \dx=\baru^n_j$. Thus, $u^{\Dx}(x,t^n)$ takes the form

\begin{equation}\label{eq:MUSCL interpolant}
\begin{split}
u^{\Dx}(x,t^n) = \sum_{j} p^n_j(x) {\mathds 1}_{[x_{\jmh}, x_{\jph}]}, \qquad\
p^n_j(x) = \baru^{n}_{j}+\frac{u'_j}{\Dx}(x-x_j),
\end{split}
\end{equation}

Here, $\displaystyle \frac{u'_j}{\Delta x}$ is an approximate slope that satisfies the second order of accuracy condition

\begin{equation}\label{eq:u slope}
\frac{u'_j}{\Dx} = \frac{\partial}{\partial x}u(x_j,t^n)+{\mathcal O}(\Dx).
\end{equation}

The approximate slopes, $\{u'_j\}$, are reconstructed from the known cell averages $\{\baru^{n}_{j}\}$, so that they maintain second-order accuracy (<ref>), while at the same time, they maintain a non-oscillatory property of the piecewise linear  reconstruction, $u^{\Dx}(x,t^n)$. A canonical  example for such reconstructed slopes is given by the so-called minmod limiter

\begin{equation}
u'_j= \minmod(\Du^n_{\jmh}, \Du^n_{\jph}), \qquad \Delta u^n_{\jph} := u^n_{j+1}-u^n_j,
\end{equation}

where $\minmod(\cdot,\cdot)$ is given by

\begin{equation}\label{eq:minmod}
\minmod(a,b) := \frac{sgn(a)+sgn(b)}{2}\cdot \min\{|a|,|b|\}.
\end{equation}

It follows that the corresponding reconstruction, $u^{\Dx}(x,t^n)$, satisfies the second-order accuracy at all by but the critical cells where $\frac{\partial}{\partial x}u(x_j,t^n)=0$, while satisfying the  TVD stability property

\begin{equation}\label{eq:TVD}
||u^{\Dx}(\cdot,t^n)||_{TV} \leq ||u^{n}||_{TV}:= \sum_j|\Delta u^n_{\jph}|.
\end{equation}

We now put this  second-order, TVD minmod-based reconstruction into action, in the context of approximate solution for one-dimensional nonlinear systems of conservation laws,

\begin{equation}\label{eq:CLS}
\frac{\partial }{\partial t} \bu(x,t) +\frac{\partial }{\partial x} {\mathbf f}(\bu(x,t)) = 0, \qquad  \bu(x,t):\Omega\times \R_+ \mapsto \R^k, \quad \Omega\subseteq\R,
\end{equation}

where  ${\mathbf f}$ is the $k$-vector of smooth flux functions. To this end, the reconstructed solution $\bu^{\Dx}(\cdot,t^n)$ is evolved to time $t^{n+1}=t^n+\Delta t$ by solving (<ref>) subject
to the piecewise linear data, $\bu(\cdot,t^n)=\bu^{\Dx}(\cdot,t^n)$, and  the resulting solution, $\bu(x,t^{n+1})$ is then realized by its cell averages, $\bu^{n+1} = \sum^{}_{\alpha} \barbu^{n+1}_\alpha {\mathds 1}_{[x_{\alpha+\hf}, x_{\alpha-\hf}]}$. This completed the cycle of Godunov-type scheme: it consists of the  the reconstruction step

\begin{equation}\tag*{(\theequation)${}_{\mathcal R}$}\label{eq:R}
{\mathcal R}: \barbu^n_\alpha= \sum^{}_{\alpha} \barbu^n_\alpha {\mathds 1}_{[x_{\alpha-\frac{1}{2}}, x_{\alpha+\frac{1}{2}}]} \mapsto \bu^{\Dx}(\cdot,t^n);
\end{equation}

for simplicity, we assume that the reconstruction step is carried out by implementing the piece-wise linear reconstruction (<ref>),(<ref>)  component-wise[In general, one makes use of local characteristic decomposition.]. This then is followed by the evolution step,

\begin{equation}\tag*{(\theequation)${}_{\mathcal E}$}\label{eq:E}
{\mathcal E}(t^n,\Delta t): \bu^{\Dx}(\cdot,t^n) \mapsto \bu^{n+1}(\cdot,t^{n+1}),
\end{equation}

and ending with averaging

\begin{equation}\tag*{(\theequation)${}_{\mathcal A}$}\label{eq:A}
{\mathcal A}: \bu^{n+1}(\cdot,t^{n+1})\mapsto \bu^{n+1} = \sum^{}_{\alpha} \barbu^{n+1}_\alpha {\mathds 1}_{[x_{\alpha-\hf},x_{\alpha+\hf}]}.
\end{equation}

It is well-known that shock discontinuities can — and in most cases of generic data —  will be formed at finite time with entropic solutions of (<ref>).   Thus, one must accept that the underlying solution is discontinuous, and a high-resolution  numerical scheme should be designed so that it prevents the spurious (Gibbs) oscillations associated with such discontinuities, particularly when high-order  (order $>1$) accuracy is involved. Since averaging operator is inherently stable, it  is the primary role of a non-oscillatory  reconstruction  step, ${\mathcal R}$, to  keep the stability  compatible with the evolution step, e.g., [53]. Back to the scalar case, $k=1$, stability is interpreted in the sense of satisfying the TVD property; then,  if the reconstruction is TVD, (<ref>) would imply the TVD stability of the overall scheme,

\[
\begin{split}
\big\|\sum^{}_{\alpha} \baru^{n+1}_\alpha {\mathds 1}_{[x_{\alpha+\hf}, x_{\alpha-\hf}]}\big\|_{TV} &\leq \|{\mathcal E}(t^n,\Delta t)u^{\Dx}(\cdot,t^n)\|_{TV} \\
& \leq \|u^{\Dx}(\cdot,t^n)\|_{TV} \\
& \leq \big\|\sum^{}_{\alpha} \baru^n_j {\mathds 1}_{[x_{\alpha-\frac{1}{2}}, x_{\alpha+\frac{1}{2}}]}\big\|_{TV}  \leq \ldots \leq \big\|\sum^{}_{\alpha} \baru^0_\alpha {\mathds 1}_{[x_{\alpha-\frac{1}{2}}, x_{\alpha+\frac{1}{2}}]}\big\|_{TV}.
\end{split}
\]

The last argument occupied  an arbitrary index $\alpha$ to emphasize that it may vary over the integers, $\alpha \leftrightarrow j\in {\mathbb Z}$, or over the half integers, $\alpha=j+\nicefrac{1}{2}$, or we can even allow staggered  grids with time-dependent shifts, $\alpha\mapsto\alpha^n=j+(\frac{1}{2})^{n-2[\nicefrac{n}{2}]}, \ j\in {\mathbb Z}$.

The TV stability encoded in the successive  inequalities above, highlights the pivotal  role of the $\minmod$ limiter in securing the TVD stability of the overall scheme.
This line of argument extends to a general  class of non-oscillatory limiters
which take the form

\begin{equation*}
u'_j = \psi(\Du^n_{\jmh}, \Du^n_{\jph}).
\end{equation*}

The corresponding Godunov-type scheme, $u^{n+1}(x)={\mathcal A}\,{\mathcal E}{\mathcal R}u^n(x)$ based on such $\psi$-limiter  amounts to a 5-cell stencil[To be precise, the resulting scheme is `essentially 3-point' stencil in the sense that the  corresponding minmod-based numerical flux  $F_{j+\frac{1}{2}}=F(u_{j-1},u_j,u_{j+1},u_{j+2})$ satisfies $F(\cdot,\bar{u},\bar{u},\cdot)=f(\bar{u})$, [39].].  The essential feature of these limiters-based reconstructions is forming difference stencils in the direction of smoothness, rather than in the direction of flow. Put differently, limiters include adaptive edge detectors   which are  necessary  in order to secure a non-oscillatory reconstruction, by avoiding `crossing' of the discontinuous data.
Higher order accuracy requires numerical derivatives with limiters which  occupy even wider stencils, e.g., [5]

\[
\psi(...,  \Du^n_{j-\frac{3}{2}},\Du^n_{\jmh}, \Du^n_{\jph},\Du^n_{j+\frac{3}{2}}, ...),
\]

which eventually led to the class of adaptive-based stencils in (W)ENO schemes, [15, 16, 45].
Back to second-order resolution, the TVD criteria has profound influence on the development of high resolution numerical schemes, starting with the second order flux-limiter  MUSCL schemes [26, 27] and following the systematic framework for TVD stability offered in [14]. A variety of limiters have been proposed in the 1970s and were studied in the context of TVD limiters since 1980s. We mention two examples from the systematic study of second-order limiters  found in [49].  First is the class of $\minmod_\theta$ limiters

\[
u'_j= \minmod_\theta(\Du^n_{\jmh}, \Du^n_{\jph}), \qquad \Delta u^n_{\jph} := u^n_{j+1}-u^n_j,
\]

where $\minmod_\theta$ is given by a one-parameter family of limiters

\begin{equation}\label{eq:theta minmod}
\minmod_\theta(a,b) := \left\{\begin{array}{ll}
\displaystyle s\cdot\min \Big(\theta |a|,\frac{|a+b|}{2},\theta |b|\Big), & \textnormal{if} \ sgn(a)=sgn(b):=s,\\
0, & \textnormal{if} \ sgn(a)+sgn(b)=0,\end{array}\right. \qquad 1\leq \theta \leq 2.
\end{equation}

The case $\theta=1$ recovers the vanilla version of minmod (<ref>).
As a second example we mention the van Albada (vA) limiter,

\begin{equation}\label{eq:VA eps limiter}
\psi^{\eps}_{\vA}(a,b) = \frac{(a^{2}+\eps^{2})b+(b^{2}+\eps^{2})a}{a^{2}+b^{2}+2\eps^{2}}, \hspace{0.5em} \eps^2 = {\mathcal O}(\Dx^{3}),
\end{equation}

This limiter was first proposed by van Albada, van Leer and Roberts [1] in 1982 and has been successfully applied in computations (see e.g. [35, 42, 47]), although a rigorous study of its TVD   stability was overshadowed by widespread applications of $\minmod$ limiter. A main aspect of this paper is to prove the TVD properties of vA limiter $\psi_{\vA}$. In fact, we demonstrate  that  when applied to the class of high resolution central schemes, the vA limiter yields comparable  and in certain cases superior performance over the $\minmod$ limiter.
This is attributed to the smoothness of the vA limiter at critical cells of local extrema[Specifically, the Lipschitz smoothness of $\psi_{\vA}(1,r)$ as a function of $r=\frac{\Du_{\jmh}}{\Du_{\jph}}$.],  in contrast to the $\minmod$ limiters, where $\minmod_\theta(\Du_{\jmh}, \Du_{\jph})=0$ whenever $\Du_{\jmh}\cdot \Du_{\jph}<0$.
Thus, the minmod-based polynomial reconstruction degenerates to first order piecewise constant interpolation at extrema cells, whereas the vA limiter still extracts a more accurate information in direction of smoothness, see (<ref>),(<ref>) below.

§ BRIEF REVIEW OF SECOND ORDER CENTRAL SCHEMES

We distinguish between two classes of Godunov type schemes, depending
on the averaging step in <ref>. In one approach, averaging at $t=t^{n+1}$ is taken at the same grid as in $t=t^n$, with $\alpha \leftrightarrow j, \ j\in{\mathbb Z}$.
This leads to a the class of  upwind schemes  which employ upstream-biased information to approximate spatial derivative.  The class of upwind scheme requires the solution of sequence of non-interacting Riemann solvers to realize the evolution step <ref>.
Wherever a second-order limiter is `turned off', the upwind scheme is reduced to the first-order Godunov scheme [11] which evolves the piecewise constant approximate solution, based on cell averages $\{\barbu^n_j\}$,

\begin{equation}\label{eq:godunov}
\begin{split}
\barbu^{n+1}_{j} &= \barbu^{n}_{j}-\frac{1}{\Dx}\int \limits^{t^{n+1}}_{t^n} \big[{\mathbf f}(\bu^{\Dx}(x_{\jph},t))-{\mathbf f}(\bu^{\Dx}(x_{\jmh},t))\big] \textnormal{d}t = \barbu^{n}_{j}-\lambda\big( {\mathbf f}(\bu^{*}_{\jph})-{\mathbf f}(\bu^{*}_{\jmh})\big),
\end{split}
\end{equation}

where $\displaystyle \lambda := \frac{\Dt}{\Dx}$ is the fixed mesh ratio. The states $\bu^{*}_{j\pm\frac{1}{2}}$ are obtained with exact or approximate Riemann solvers, see e.g. [7, 17, 41, 6, 56]. When applied to systems of equations, both left- and right-running waves may exist at interfaces. Hence the local characteristic decomposition of the flux ${\mathbf f}$ is required for the upwind constructions, based on  the “direction of the wind" identified by characteristic decomposition into local eigen-fields.

The class of central schemes
keeps the averaging over staggered grid, $\alpha\leftrightarrow j+\nicefrac{1}{2}, \ j\in {\mathbb Z}$, so that it process the information from upstream and downstream in an averaged, “central” manner, [36, 55].
The evolution step in the class of central schemes requires simple quadrature; no Riemann solvers are needed. However, wherever a second-order limiter is `turned off', it is reduced to the diffusive Lax-Friedrichs  (LxF) scheme [25]:

\begin{equation}\label{eq:LxF}
\barbu^{n+1}_{\jph} = \frac{\barbu^{n}_{j}+\barbu^{n}_{j+1}}{2}-\lambda({\mathbf f}(\bu^{n}_{j+1})-{\mathbf f}(\bu^{n}_{j})).
\end{equation}

In contrast to the Godunov scheme (<ref>), the LxF solver computes the cell averages on the staggered mesh, $\displaystyle \barbu(x_{\jph},t):= \frac{1}{\Dx}\int^{x_{j+1}}_{x_{j}}\bu(x,t)\dx$, recalling remark <ref>. The integration of equation (<ref>) over the staggered control volume $[x_{j}, x_{j+1}]\times [t^n,t^{n+1}]$ yields

\begin{equation}\label{eq:staggered cell average}
\barbu(x_{\jph},t^{n+1}) =\frac{1}{\Dx}\int_{x_{j}}^{x_{j+1}} \bu(x,t^{n}) \dx- \frac{1}{\Dx}\int^{t^{n+1}}_{t^n} [{\mathbf f}(\bu(x_{j+1},t))-{\mathbf f}(\bu(x_{j},t))] \dt.
\end{equation}

It is clear that the LxF scheme is a first order approximation to (<ref>). Compared with the upwind framework, the advantage of central differencing is that the solution is smooth in the neighborhood of points $\{x_{j}\}$. Indeed, under a sufficiently small time step, $\Dt$, the Riemann waves initiating from  $x_{\jph}$ do not affect the adjacent cell centers within the time slab, $t\in[t^n, t^{n+1}]$. Hence the costly characteristic decompositions can be avoided.

In spite of its simplicity, the main disadvantage of the LxF scheme is the excessive numerical dissipation which reduces its resolution at discontinuities and smooth extrema. The situation is dramatically improved by applying higher order piecewise polynomial interpolations, leading  to the second-order Nessyahu-Tadmor scheme[At the end, all fully-discrete second-order accurate schemes are one version of another of the Lax-Wendroff scheme.]. Still, wherever  a minmod limiter is `turned off',
the second-order central scheme is reduced to the difussive LxF scheme.
This is precisely why  the gain of a smooth vA limiter over the minmod limiter is more noticeable in the context of central schemes.
We turn to briefly review the second order fully discrete and semi-discrete central schemes.

§.§ The Nessyahu-Tadmor scheme

The second-order fully discrete central scheme was first proposed in 1990 [36]. The second-order of accuracy is obtained by combining the ${\mathcal R}\,{\mathcal E}{\mathcal A}$ steps. First, we evolve   a reconstructed piecewise linear solution $u(x,t^n)$  (<ref>): integration of  the equation (<ref>) over the rectangle $[x_{j}, x_{j+1}]\times [t^n,t^{n+1}]$ yields

\begin{equation*}
\barbu^{n+1}_{\jph} = \frac{1}{\Dx}\Big[\int^{x_{\jph}}_{x_{j}}\bp^n_{j}(x)\dx+\int^{x_{j+1}}_{x_{\jph}}\bp^n_{j+1}(x)\dx \Big]-\frac{1}{\Dx}\int^{t^{n+1}}_{t^n} \Big[{\mathbf f}(\bu^{\Dx}(x_{j+1},t))-{\mathbf f}(\bu^{\Dx}(x_{j},t))\Big] \dt.
\end{equation*}

Then, the time integral of the flux function can approximated by the midpoint rule at the expense of ${\mathcal O}(\Dt^3)$ local truncation error. This results in the Nessyahu-Tadmor (NT) scheme

\begin{equation}\label{eq:NT scheme}
\barbu^{n+1}_{\jph} = \frac{\barbu^{n}_{j}+\barbu^{n}_{j+1}}{2}+\frac{1}{8}(\bu'_j-\bu'_{j+1})-\lambda[{\mathbf f}(\bu^{\nph}_{j+1})-{\mathbf f}(\bu^{\nph}_{j})],
\end{equation}

The midpoint values, $\bu^{\nph}_{j}$, are predicted by Taylor expansions

\begin{equation*}
\bu^{\nph}_{j} = \barbu^n_{j}-\frac{\lambda}{2} {\mathbf f}'_j,
\end{equation*}

where the numerical derivative $f'_j$ satisfies the second order accuracy condition

\begin{equation}\label{eq:f slope}
\frac{{\mathbf f}'_j}{\Dx} = \frac{\partial }{\partial x} {\mathbf f}(\barbu^n_j)+{\mathcal O}(\Dx).
\end{equation}

Different options of evaluating ${\mathbf f}'_j$ were proposed in [36]. For example, the numerical derivative can be computed in terms of exact flux Jacobian,

\begin{equation}\label{eq:f slope jacob}
{\mathbf f}'_j = A(\barbu^n_{j})\bu'_j, \qquad A(\bu) := \frac{\partial {\mathbf f}(\bu)}{\partial \bu}.
\end{equation}

Alternatively, one can apply Jacobian-free approximations when the characteristic decomposition is computationally expensive (even sometimes inaccessible), e.g.

\begin{equation}\label{eq:f slope jacob-free}
{\mathbf f}'_j = \minmod(\Dbf_{\jmh}, \Dbf_{\jph}), \qquad \Dbf_{\jph} := {\mathbf f}(\barbu^n_{j+1})-{\mathbf f}(\barbu^n_j).
\end{equation}

Numerical results in [19, 33] have shown that the Jacobian-free version of the central scheme does not deteriorate the high resolution. When applied to systems of equations, the NT scheme inherits the simplicity of the LxF solver, i.e. the constructions can be extended in a component-wise manner without the use of Riemann solvers or characteristic decompositions.

The multidimensional formulation of the NT scheme was obtained with similar integration procedures, details can be found in [19, 2, 4]. The order of accuracy can be further improved by applying higher order piecewise polynomial interpolants, see e.g. [33, 28, 29, 30]. Related references can be found in [3].

§.§ The semi-discrete formulation

It can be shown that the numerical dissipation of the NT scheme has order ${\mathcal O}((\Dx)^{4} /\Dt)$. In the convective problems where $\Dt\sim \Dx$, the NT scheme achieves higher resolution than the first order LxF scheme due to the reduced numerical viscosity. However, it is noticed that the NT scheme and its higher order extensions do not admit any semi-discrete limits, hence they are not appropriate for small time step computations or steady-state calculations. This motivated the development of semi-discrete central schemes. We sketch the derivations along the lines of [24, 23].

We start with a piecewise polynomial approximation $\bu^{\Dx}(x,t^n)\approx u(x,t^n)$ of the form

\begin{equation*}
\bu^{\Dx}(x,t^n) = \sum_{j} \bp^n_{j}(x){\mathds 1}_{[x_{\jmh}, x_{\jph}]}.
\end{equation*}

The polynomials $p^n_j(x)$ should have desired order of accuracy and conserve the cell averages $\barbu^n_j := \baru(x_j,t^n)$, i.e.

\begin{equation*}
\frac{1}{\Dx} \int^{x_{\jph}}_{x_{\jmh}}\bp^n_j(x) \dx = \barbu^n_j.
\end{equation*}

We denote the reconstructed variables at $x_{\jph}$ from the left and the right by

\begin{equation}
\bu^{-}_{\jph} := \bp^n_j(x_{\jph}), \qquad \bu^{+}_{\jph} := \bp^{n}_{j+1}(x_{\jph}).
\end{equation}

In particular, when the  MUSCL interpolant (<ref>) is applied we have

\begin{equation}\label{eq:MUSCL reconstruction}
\bu^{-}_{\jph}: = \barbu^n_{j}+\frac{1}{2}\bu'_{j}, \hspace{1em} \bu^{+}_{\jph}: = \barbu^n_{j+1}-\frac{1}{2}\bu'_{j+1}.
\end{equation}

To obtain a semi-discrete formulation the key idea is to take cell averaging over narrower control volumes. Assume that the maximal forward and backward wave speeds at $x_{\jph}$ are estimated by $a^{+}_{\jph}\geq 0$ and $a^{-}_{\jph}\leq 0$. At the next time level $t^{n+1}:= t^n+\Dt$, the region influenced by the Riemann fan originating at $x_{\jph}$ is approximated with the interval, $I_{\jph} := [x^{-}_{\jph}, x^{+}_{\jph}]$, where

\begin{equation*}
x^{+}_{\jph}:= x_{\jph}+a^{+}_{\jph}\Dt, \hspace{1em}  x^{-}_{\jph}:= x_{\jph}+a^{-}_{\jph}\Dt.
\end{equation*}

It is clear that under a sufficiently small time step, the non-smooth regions $I_{\jph}$ are separated due to the finite wave speeds. In this way, the width of Riemann fan is bounded by $(a^{+}_{\jph}-a^{-}_{\jph})\Dt$, in contrast to the fixed width $\Dx$ in the fully discrete counterpart.

Now we reconstruct a non-oscillatory, conservative, piecewise polynomial interpolant at $t = t^{n+1}$,

\begin{equation}
\widetilde{\bw}^{n+1}(x) = \sum_{j}\Big[\widetilde{\bw}^n_{j}(x){\mathds 1}_{[x^{+}_{\jmh}, x^{-}_{\jph}]}+\widetilde{\bw}^n_{\jph}(x){\mathds 1}_{ [x^{-}_{\jph}, x^{+}_{\jph}]}\Big].
\end{equation}

The polynomials $\widetilde{\bw}^{n+1}_{j}$ and $\widetilde{\bw}^{n+1}_{\jph}$ conserve the cell averages, $\bar{\bw}^{n+1}_{j}$ and $\bar{\bw}^{n+1}_{\jph}$, over the smooth and non-smooth domains. The values of $\bar{\bw}^{n+1}_{j}$ and $\bar{\bw}^{n+1}_{\jph}$ are computed by integrating (<ref>) over the rectangular domains, $[x^{+}_{\jmh}, x^{-}_{\jph}]\times[t^n,t^{n+1}]$ and $[x^{-}_{\jph}, x^{+}_{\jph}]\times[t^n,t^{n+1}]$, respectively.

It can be shown that

\begin{equation}\label{eq:smooth average}
\begin{split}
\frac{x^{-}_{\jph}-x^{+}_{\jmh}}{\Dx}\bar{w}^{n+1}_{j} =  \baru^{n}_{j}+\lambda(a^{-}_{\jph} u^{-}_{\jph}-a^{+}_{\jmh} u^{+}_{\jmh})-\lambda[f(u^{-}_{\jph})-f(u^{+}_{\jmh})] +{\mathcal O}(\frac{\Dt^2}{\Dx}),
\end{split}
\end{equation}

\begin{equation}\label{eq:non-smooth average}
\bar{w}^{n+1}_{\jph}  = \frac{a^{+}_{\jph}u^{+}_{\jph}-a^{-}_{\jph}u^{-}_{\jph}}{a^{+}_{\jph}-a^{-}_{\jph}}-\frac{f(u^{+}_{\jph})-f(u^{-}_{\jph})}{a^{+}_{\jph}-a^{-}_{\jph}}+{\mathcal O}(\Dt).
\end{equation}

Notice that the width of Riemann fans is bounded with ${\mathcal O}(\Dt)$, we have

\begin{equation}\label{eq:non-smooth poly}
\tilde{w}^n_{j\pm\frac{1}{2}}(x) = \bar{w}^{n}_{j\pm\frac{1}{2}}+{\mathacal O}(\Dt), \hspace{1em}\forall x\in [x^{-}_{j\pm\frac{1}{2}}, x^{+}_{j\pm\frac{1}{2}}].
\end{equation}

Finally, the cell average at $t^{n+1}$ is computed with

\[
\barbu^{n+1}_{j} = \frac{1}{\Dx}\int^{x_{\jph}}_{x_{\jmh}}\widetilde{\bw}^{n+1}(x)\dx.
\]

Passing the limit $\Dt\rightarrow0$ yields a semi-discrete conservative scheme

\begin{equation}\label{eq:semi-discrete ode}
\frac{\textnormal{d}}{\dt}\barbu_{j}(t) := \lim_{\Dt\rightarrow 0} \frac{\barbu^{n+1}_{j}-\barbu^n_{j}}{\Dt} = -\frac{{\mathbf F}_{\jph}(t)-{\mathbf F}_{\jmh}(t)}{\Dx},
\end{equation}

where the numerical fluxes, ${\mathbf F}_{\jph}(t)$, are given by

\begin{equation}\label{eq:semi-discrete flux}
{\mathbf F}_{\jph}(t) := \frac{a^{+}_{\jph}{\mathbf f}(\bu^{-}_{\jph})-a^{-}_{\jph}{\mathbf f}(\bu^{+}_{\jph})}{a^{+}_{\jph}-a^{-}_{\jph}}+\frac{a^{+}_{\jph}a^{-}_{\jph}}{a^{+}_{\jph}-a^{-}_{\jph}}(\bu^{+}_{\jph}-\bu^{-}_{\jph}).
\end{equation}

There are different options to estimate the wave speeds $a^{\pm}_{\jph}$. The Kurganov-Tadmor (KT) scheme proposed in [24], uses the spectral radius of the Jacobian, $\displaystyle A(\bu):=\frac{\partial {\mathbf f}(\bu)}{\partial \bu}$,

\begin{equation*}
a^{+}_{\jph} = -a^{-}_{\jph} = a_{\jph} := \max\Big\{\rho(A(\bu^{-}_{\jph})), \rho(A(\bu^{+}_{\jph}))\Big\},
\end{equation*}

where $\rho(\cdot)$ represents the spectral radius of a matrix. Then the numerical flux can be expressed in terms of the Rusanov flux,

\begin{equation}\label{eq:Rus flux}
{\mathbf F}_{\jph}(t) = {\mathbf F}^{Rus}(\bu^{-}_{\jph}, \bu^{+}_{\jph}):= \frac{{\mathbf f}(\bu^{-}_{\jph})+{\mathbf f}(\bu^{+}_{\jph})}{2}-\frac{a_{\jph}}{2}(\bu^{+}_{\jph}-\bu^{-}_{\jph}).
\end{equation}

For strictly hyperbolic problems, the flux Jacobian $\displaystyle \frac{\partial {\mathbf f}}{\partial \bu}$ has $N$ distinct eigenvalues $\lambda_1<\cdots<\lambda_{N}$. The semi-discrete central-upwind scheme version of KNP, [23], employs a more accurate estimate of the wave speeds,

\begin{equation}\label{eq:hyper wave speeds}
\begin{split}
& a^{+}_{\jph} = \max\Big\{\lambda_{N}(A(\bu^{-}_{\jph})), \lambda_{N}(A(\bu^{+}_{\jph})),0\Big\},\\
& a^{-}_{\jph} = \min\Big\{\lambda_{1}(A(\bu^{-}_{\jph})), \lambda_{1}(A(\bu^{+}_{\jph})),0\Big\}.
\end{split}
\end{equation}

In this way, the scheme reduces the numerical dissipation by employing the HLL flux [18].

In the computations, the set of ODEs (<ref>) is integrated with an appropriate ODE solver. To preserve the overall high order of accuracy as well as the non-oscillatory property, one may apply a higher order Strong Stability Preserving (SSP) Runge-Kutta method, see e.g. [44, 13, 48, 12].

The semi-discrete central schemes retain the advantage of being Riemann-solver-free, hence the componentwise extension is allowed when solving systems of equations. In multidimensional problems, the second order constructions can be extended in a dimension-by-dimension manner. Third order extensions have also been derived with the help of piecewise parabolic interpolants, we refer the details to [24, 22].

§ ANALYSIS OF THE SMOOTH VAN ALBADA LIMITER

In this section, we study the analytical properties of the smooth van Albada limiter (<ref>). To get more insight, the limiter can be written as [26]

\begin{equation*}
\psi^{\eps}_{\vA}(a,b) = \frac{a+b}{2}\Big(1-\frac{(a-b)^2}{a^2+b^2+2\eps^2}\Big).
\end{equation*}

In the smooth regions where $a\approx b$, the limiter tends to recover the second order central finite differencing, $\displaystyle \frac{a+b}{2}$. Across the discontinuities, however, the averaged slope is biased to the smallest value among the two one-sided slopes. These mechanisms are expected to ensure the second order accuracy and prevent the undesirable numerical oscillations. We will provide rigorous proofs for these advantages. For convenience of discussion, we will consider the simplified version without $\eps$

\begin{equation}\label{eq:VA limiter}
\psi_{\vA}(a,b) = \frac{a^2 b+a b^2}{a^2+b^2}, \hspace{1em} a^2+b^2\neq 0.
\end{equation}

Indeed, as reported in [1], the computed solutions are not sensitive to the specific value of $\eps$.

§.§ The second order accuracy

We show that the van Albada limiter (<ref>) ensures a second order spatial accuracy. Back in the scalar framework, we first consider the reconstruction based on the point values, $u_j := u(x_j)$. Denote the left- and the right- sided slopes by

\begin{equation*}
\parm u_j:= \frac{u_{j}-u_{j-1}}{\Dx} \hspace{1em} \text{and} \hspace{1em} \parp u_j:= \frac{u_{j+1}-u_{j}}{\Dx}.
\end{equation*}

The numerical derivative constructed with the van Albada limiter (<ref>) is given by

\begin{equation*}
\frac{u'_j}{\Dx} = \frac{(\parm u_{j})(\parp u_{j})^{2}+(\parp u_{j})(\parm u_{j})^{2}}{(\parp u_{j})^{2}+(\parm u_{j})^{2}}.
\end{equation*}

Taylor expansion yields

\begin{equation}\label{eq:VA-point}
\frac{u'_j}{\Dx} = \parx u_{j} + \Big(\frac{1}{6}\parx^{3} u_{j}-\frac{1}{2}\frac{(\parx^{2} u_{j})^{2}}{\parx u_{j}}\Big)\Dx^{2}+{\mathcal O}(\Dx^{3}).
\end{equation}

The second order condition (<ref>) is satisfied in regions away from  critical points, $|\parx u_j|\gg \Dx$. In the framework of finite volume methods, the slope reconstruction is based on the cell averages, $\displaystyle \baru_j := \frac{1}{\Dx}\int^{x_{\jph}}_{x_{\jmh}}u(x)\dx$. Then the Taylor expansion is slightly modified,

\begin{equation}\label{eq:VA-cell}
\frac{u'_j}{\Dx} = \parx u_{j}+\Big(\frac{1}{6}\parx^{3}u_{j}+\frac{1}{24}\parx^{4}u_{j}-\frac{1}{2}\frac{\parx^{2}u_{j}}{\parx u_{j}}\Big)\Dx^{2}+{\mathcal O}(\Dx^{3}), \qquad |\parx u_j|\gg \Dx.
\end{equation}

For both point-value-based and cell-average-based reconstructions, the numerical derivative, $u'_j/\Dx$, approximates the exact derivative, $\parx u_j$, with order ${\mathcal O}(\Dx^2)$ in the non-critical regions. Near the critical points, the second order condition (<ref>) is satisfied up to $\parx u \approx \Dx$. Hence it is expected that the van Albada limiter introduces less dissipation than the minmod limiter at discontinuities and smooth extrema. Their numerical performances will be compared in section 4.

§.§ The non-oscillatory property

We are going to show that the van Albada limiter, as applied to the central schemes, generates non-oscillatory solutions. We begin with the following lemma.

The approximate slope $u'_j = \psi_{\vA}(\Du_{\jmh},\Du_{\jph})$ satisfies the following estimates

\begin{equation}\label{eq:VA bound}
\frac{1-\sqrt{2}}{2}\leq \frac{u_{j}'}{\Du_{j\pm\frac{1}{2}}} \leq \frac{1+\sqrt{2}}{2}, \hspace{1em} \end{equation}

\begin{equation}\label{eq:VA diff}
\frac{|u'_{j+1}-u'_{j}|}{|\Du_{\jph}|} \leq \sqrt{2}.
\end{equation}

Proof. Denote $\displaystyle r := \frac{\Du_{\jmh}}{\Du_{\jph}}$, then we can write

\begin{equation*}
u'_{j} = \frac{r^{2}+r}{r^{2}+1} \Du_{\jph} = \frac{r+1}{r^{2}+1}\Du_{\jmh}.
\end{equation*}

We will show that $\displaystyle \frac{1-\sqrt{2}}{2}\leq \frac{r^{2}+r}{r^{2}+1}, \hspace{0.2em}\frac{r+1}{r^{2}+1} \leq \frac{1+\sqrt{2}}{2}$.

Assume that $\displaystyle\frac{r^{2}+r}{r^{2}+1}$ is bounded from below and above by $m$ and $ M\in\R$ respectively,

\[
m\leq \frac{r^{2}+r}{r^{2}+1} \leq M, \hspace{0.5em} \forall r\in\R.
\]

This can be equivalently expressed as

\begin{equation}
\left\{\begin{array}{l}
(M-1)r^{2}-r+M \geq 0,\\
\\
(m-1)r^{2}-r+m \leq 0,
\end{array}\right. \qquad \forall r\in\R.
\end{equation}

The upper bound, $M$, should satisfy the conditions

\begin{equation*}
\left\{\begin{array}{l}
M>1,  \\
\\
\Delta =  1-4M(M-1)\leq 0.
\end{array}\right.
\end{equation*}

Solving these inequalities yields $\displaystyle M\geq \frac{1+\sqrt{2}}{2}$ and hence $\displaystyle \frac{1+\sqrt{2}}{2}$ gives the supremum of $\displaystyle \frac{r^2+r}{r^2+1}$. On the other hand, the lower bound, $m$, should satisfy

\begin{equation*}
\left\{\begin{array}{l}
m<1,  \\
\\
\Delta =  1-4m(m-1)\leq 0.
\end{array}\right.
\end{equation*}

These conditions yield $\displaystyle m\leq\frac{1-\sqrt{2}}{2}$ and hence $\displaystyle \frac{1-\sqrt{2}}{2}$ gives the infimum of $\displaystyle \frac{r^2+r}{r^2+1}$. Combining these results, we have

\begin{equation*}
\frac{1-\sqrt{2}}{2}\leq \frac{u'_j}{\Du_{\jph}} = \frac{r^2+r}{r^2+1} \leq \frac{1+\sqrt{2}}{2}.
\end{equation*}

With similar argument, we can show that $\displaystyle \frac{1-\sqrt{2}}{2}\leq\frac{r+1}{r^{2}+1} \leq \frac{1+\sqrt{2}}{2}$ and therefore

\begin{equation*}
\frac{1-\sqrt{2}}{2}\leq \frac{u'_j}{\Du_{\jmh}} = \frac{r+1}{r^2+1} \leq \frac{1+\sqrt{2}}{2}.
\end{equation*}

This completes the proof of (<ref>). The second estimate (<ref>) is a direct consequence of (<ref>).
$\hfill \square$

With the help of Lemma <ref>, we discuss the non-oscillatory property of the fully discrete and the semi-discrete schemes separately.

§.§.§ The fully discrete scheme

We consider the fully discrete NT scheme (<ref>). The scheme can be written in a conservative form [36][Recall that $\lambda$ is the fixed mesh-ratio, $\lambda=\frac{\Delta t}{\Delta x}$, which should be distinguished from the eigenvalues in (<ref>), $\lambda_1<\ldots <\lambda_N$.]

\begin{equation}
\label{eq:NT conserv}
\baru_{\jph}^{n+1} = \frac{\baru^{n}_{j}+\baru^{n}_{j+1}}{2}-\lamb(g_{j+1}-g_{j}), \qquad \lambda=\frac{\Delta t}{\Delta x},
\end{equation}

with the so-called modified numerical flux, $g_j$,  given by

\begin{equation}\label{eq:NT modified flux}
g_{j} = f(u^{\nph}_{j})+\frac{1}{8\lamb}u'_{j}.
\end{equation}

The NT scheme combined with the van Albada limiter (<ref>) can be shown to be TVD following similar approach in [39, 36]. The proof is based on the following lemma.

The scheme (<ref>) is TVD if the numerical flux, $g_j$, satisfies the following generalized condition

\begin{equation}\label{eq:NT TVD condition}
\lambda \Big|\frac{\Dg_{\jph}}{\Du^n_{\jph}}\Big|\leq\frac{1}{2}, \hspace{1em} \Dg_{\jph} := g_{j+1}-g_{j}.
\end{equation}

Proof. By (<ref>), the difference $\baru^{n+1}_{\jph}-\baru^{n+1}_{\jmh}$ amounts to

\begin{equation*}
\baru^{n+1}_{\jph}-\baru^{n+1}_{\jmh} = \Big(\frac{1}{2}-\lamb\frac{\Dg_{\jph}}{\Du_{\jph}}\Big)\Du_{\jph}+    \Big(\frac{1}{2}+\lamb\frac{\Dg_{\jmh}}{\Du_{\jmh}}\Big)\Du_{\jmh}.
\end{equation*}

Condition (<ref>) implies that the coefficients in the parenthesis are positive. Hence,

\begin{equation*}
\begin{split}
\sum_{j}|\baru^{n+1}_{\jph}-\baru^{n+1}_{\jmh}| & \leq \sum_{j} \Big(\frac{1}{2}-\lamb\frac{\Dg_{\jph}}{\Du_{\jph}}\Big)|\Du_{\jph}|+    \Big(\frac{1}{2}+\lamb\frac{\Dg_{\jmh}}{\Du_{\jmh}}\Big)|\Du_{\jmh}|\\
& = \sum_{j} \Big(\frac{1}{2}-\lamb\frac{\Dg_{\jph}}{\Du_{\jph}}\Big)|\Du_{\jph}|+    \Big(\frac{1}{2}+\lamb\frac{\Dg_{\jph}}{\Du_{\jmh}}\Big)|\Du_{\jph}|\\
& = \sum_{j}|\Du_{\jph}|.
\end{split}
\end{equation*}

We conclude $TV(\baru^{n+1})\leq TV(\baru^n)$. $\square$

Equipped with Lemma <ref> and <ref>, we prove the NT scheme is TVD.

Consider the scalar NT scheme  (<ref>) using numerical slope $u'_j$  constructed with the van Albada limiter (<ref>), and let the flux numerical derivative be chosen by $f'_j = a(\baru^n_j)u'_j$ with $a(u) := f'(u)$. Assume that the following condition is satisfied

\begin{equation}\label{eq:NT CFL}
\lambda\max_{u}|a(u)|\leq \frac{-1+\sqrt{\frac{3}{2}+\frac{3\sqrt{2}}{4}}}{1+\sqrt{2}} \approx 0.24.
\end{equation}

Then the NT scheme (<ref>) is TVD.

Proof. Denote $\displaystyle \CFL:= \lambda\max_{u}|a(u)|$. By (<ref>) we have

\begin{equation}\label{eq:NT decomp}
\lamb\Big|\frac{\Dg_{\jph}}{\Du_{\jph}}\Big|\leq \underbrace{\lamb\Big|\frac{f(u^{\nph}_{j+1})-f(u^{\nph}_{j})}{u^{\nph}_{j+1}-u^{\nph}_{j}}\Big|}_{I}
\underbrace{\Big|\frac{u^{\nph}_{j+1}-u^{\nph}_{j}}{\Du_{\jph}}\Big|}_{II}+
\underbrace{\frac{1}{8}\Big|\frac{\Delta u'_{\jph}}{\Du_{\jph}}\Big|}_{III}.
\end{equation}

By definition, the first term on the right of (<ref>) is bounded by

\begin{equation}\label{eq:NT est I}
\lamb\left|\frac{f(u^{\nph}_{j+1})-f(u^{\nph}_{j})}{u^{\nph}_{j+1}-u^{\nph}_{j}}\right| \leq \lambda\max_u |a(u)| := \CFL.
\end{equation}

Notice that (<ref>) implies $\displaystyle \Big|\frac{u'_j}{\Du_{j\pm\frac{1}{2}}}\Big|\leq \frac{1+\sqrt{2}}{2}$, we can estimate the second term on the right of (<ref>),

\begin{equation}\label{eq:NT est II}
\begin{split}
\Big|\frac{u^{\nph}_{j+1}-u^{\nph}_{j}}{\Du_{\jph}}\Big|& \leq 1+\frac{\lamb}{2}\Big|\frac{\Df'_{\jph}}{\Du_{\jph}}\Big|
\leq 1+\frac{\lamb}{2} \frac{|a_{j+1}u_{j+1}'|+|a_{j}u'_{j}|}{|\Du_{\jph}|}\\
& \leq 1+\frac{\lamb\max_{u}|a(u)|}{2} \cdot \frac{|u_{j+1}'|+|u'_{j}|}{|\Du_{\jph}|}  \leq 1+\frac{1+\sqrt{2}}{2}\CFL .
\end{split}
\end{equation}

As a direct consequence of (<ref>), the third term on the right of (<ref>) does not exceed

\begin{equation}\label{eq:NT est III}
\frac{1}{8}\Big|\frac{\Delta u'_{\jph}}{\Du_{\jph}}\Big|\leq\frac{\sqrt{2}}{8}.
\end{equation}

Using (<ref>), (<ref>) and (<ref>), the TVD condition (<ref>) boils down to

\begin{equation*}
\CFL\Big(1+\frac{1+\sqrt{2}}{2}\CFL\Big)+\frac{\sqrt{2}}{8}\leq \frac{1}{2},
\end{equation*}

which in turn  recovers the CFL condition (<ref>). $\square$

The TVD stability is not sharp in the sense that the condition  (<ref>) is to restrictive, and serves here only as a theoretically sufficient bound. The numerical experiments reported in section <ref>, indicate that the computed NT solutions with vA limiter remain non-oscillatory as long as $\displaystyle \CFL\leq \nicefrac{1}{2}$.

§.§.§ The semi-discrete scheme

The TVD stability of the van Albada limiter as applied to the non-staggered semi-discrete schemes (<ref>) can be shown following the lines of [38, 51]. The key ingredient is to represent the scheme in an appropriate incremental form which meets a certain positive condition. We prove the following result with the help of Lemma <ref>.

Let the states $u^{\pm}_{\jph}$ be computed from the MUSCL reconstruction (<ref>), with the numerical slope $u'_j$ constructed with the van Albada limiter (<ref>). Consider a generalized scalar MUSCL scheme

\begin{equation}\label{eq:generalized MUSCL scheme}
\frac{d}{dt}\baru_j(t) = -\frac{1}{\Dx}[F(u^{-}_{\jph}, u^{+}_{\jph})-F(u^{-}_{\jmh}, u^{+}_{\jmh})],
\end{equation}

and let the numerical flux $F(\cdot,\cdot)$ be monotone and consistent,

\begin{equation}\label{eq:flux monotone}
\frac{\partial}{\partial u} F(u,v) \geq 0 \quad \text{and} \quad \frac{\partial}{\partial v} F(u,v) \leq 0,
\end{equation}

\begin{equation}\label{eq:flux consistency}
F(u,u) = f(u).
\end{equation}

Then the scheme (<ref>) is TVD,

\begin{equation*}
\frac{d}{dt} [TV(\baru(t))] \leq 0.
\end{equation*}

Proof. We argue along the lines of [51], expressing the scheme (<ref>) in the incremental form

\begin{equation}\label{eq:semi-discrete incremental}
\frac{d}{dt}\baru_{j} = -\frac{1}{\Dx_{j}}C_{\jmh}\Du_{\jmh}+\frac{1}{\Dx_{j}}D_{\jph}\Du_{\jph},
\end{equation}

\begin{equation}\label{eq:incremental coeffs}
\begin{split}
&C_{\jmh} = \frac{1}{\Du_{\jmh}}[F(u^{-}_{\jph},u^{+}_{\jmh})-F(u^{-}_{\jmh}, u^{+}_{\jmh})], \\
&D_{\jph} = -\frac{1}{\Du_{\jph}}[F(u^{-}_{\jph}, u^{+}_{\jph})-F(u^{-}_{\jph},u^{+}_{\jmh})].
\end{split}
\end{equation}

We denote

\begin{equation*}
\begin{split}
&s_{\jph} := sign(\Du_{\jph}) = \left\{\begin{array}{cc}
1     &  \Du_{\jph}>0\\
\pm 1    & \Du_{\jph}=0\\
-1     & \Du_{\jph}<0,
\end{array}\right.  \\
&\chi_{j} := 1-s_{\jmh}s_{\jph} = \left\{\begin{array}{ll}
2,     &  \Du_{\jmh} \cdot \Du_{\jph}\leq 0\\
0,     & \Du_{\jmh} \cdot \Du_{\jph}> 0.
\end{array}\right.
\end{split}
\end{equation*}

Forward differencing of (<ref>) gives

\begin{equation}\label{eq:incremental forward diff}
\begin{split}
\frac{d}{dt}\Du_{\jph} &= \Big(\frac{1}{\Dx_{j+1}}D_{j+\frac{3}{2}}\Du_{j+\frac{3}{2}}-\frac{1}{\Dx_{j+1}}C_{\jph}\Du_{\jph}\Big)\\
&\qquad -\Big(\frac{1}{\Dx_{j}}D_{\jph}\Du_{\jph}-\frac{1}{\Dx_{j}}C_{\jmh}\Du_{\jmh}\Big).
\end{split}
\end{equation}

Multiply (<ref>) by $s_{\jph}$ and sum by parts, we have

\begin{equation*}
\begin{split}
\frac{d}{dt}TV(\baru(t)) & = \sum_{j} s_{\jph} \frac{d}{dt}\Du_{\jph}\\
& = -\sum_{j} \frac{1}{\Dx}[(s_{\jmh}-s_{\jph})C_{\jmh}\Du_{\jmh}+(s_{\jph}-s_{\jmh})D_{\jph}\Du_{\jph}]\\
& = -\sum_{j} \frac{\chi_{j}}{\Dx}[C_{\jmh}|\Du_{\jmh}|+D_{\jph}|\Du_{\jph}|].
\end{split}
\end{equation*}

Here we use the property that $s^{2}_{\jph} \equiv 1$. Since $\chi_j\geq 0$, it suffices to show that $C_{\jmh}$, $ D_{\jph}\geq 0$.
Indeed, we find that

\begin{equation}
\begin{split}
&C_{\jmh} = \partial_{u}F(\tilde{u}_1, u^{+}_{\jmh})\Big(1+\frac{u'_j-u'_{j-1}}{2\Du_{\jmh}}\Big), \\
&D_{\jph} = -\partial_{v}F(u^{-}_{\jph}, \tilde{u}_2)\Big(1-\frac{u'_{j+1}-u'_{j}}{2\Du_{\jph}}\Big),
\end{split}
\end{equation}

where $\tilde{u}_1$ is between $u^{-}_{\jmh}$ and $u^{-}_{\jph}$, $\tilde{u}_2$ is between $u^{+}_{\jmh}$ and $u^{+}_{\jph}$. From (<ref>) we know that the terms in the parentheses are positive. Hence the non-negativity of $C_{\jmh}$ and $D_{\jph}$ follows from (<ref>).

The semi-discrete central scheme (<ref>) can be considered a generalized MUSCL scheme, the associated numerical flux (<ref>) can expressed in the form $F_{\jph}(t) := F(u^{-}_{\jph}, u^{+}_{\jph})$, with

\begin{equation*}
F(u,v) = \frac{a^{+}f(u)-a^{-}f(v)}{a^{+}-a^{-}}+\frac{a^{+}a^{-}}{a^{+}-a^{-}}(v-u).
\end{equation*}

The monotonicity of $F(\cdot,\cdot)$ is clear since $a^- \leq f'(u), f'(v) \leq a^+$. Hence the TVD stability of the van Albada limiter as applied to the semi-discrete central scheme (<ref>) follows from Theorem <ref>.

The semi-discrecte scheme (<ref>) (or equivalently (<ref>)) advanced with foward Euler time-stepping reads

\begin{equation}\label{eq:MUSCL FE}
\baru^{n+1}_j = \baru^n_j-\lambda(F(u^{-}_{\jph}, u^{+}_{\jph})-F(u^{-}_{\jmh}, u^{+}_{\jmh})),
\end{equation}

where $u^{+}_{\jph} := \baru^{n}_j-\frac{1}{2}u'_{j+1}$ and $u^{-}_{\jph} := \baru^{n}_j+\frac{1}{2}u'_{j}$. With similar incremental representation in Theorem <ref>, it can be proved that the van Albada limiter as applied to the fully discrete non-staggered scheme (<ref>) generates TVD solutions. Applying a higher order SSP time integration method would achieve a higher order of accuracy in time without deteriorating the stability.

§ NUMERICAL EXPERIMENTS

In this section, we examine the performance of vA limiter (<ref>) with the small[Even smaller than the one indicated in (<ref>); this does not seem to affect the results.] bias, $\eps = (\Dx)^3$,  compared vs. the minmod limiter (<ref>) which  is taken, unless otherwise stated, in its vanilla version $\theta=1$. The advantage of the vA limiter over minmod${}_1$ is apparent  for both — the fully discrete NT scheme (<ref>) and the semi-discrete central-upwind  version of KNP, (<ref>), (<ref>).

In all the simulations of the NT scheme, the numerical slope of flux, ${\mathbf f}'_j$ is evaluated with the exact flux Jacobian (<ref>), and the CFL number is taken to be 0.45. The implementation of the semi-discrete scheme requires a high order time discretization. We use the third order, explicit SSPRK3 method for the time integration of method-of-lines ODEs (<ref>) with the CFL number chosen to be 0.7. To better ensure the robustness in the strong shock problems, the MUSCL reconstruction for systems of equations is performed on the characteristic variables, e.g., <cit.>, <cit.>[Here we use standard characteristic decomposition, $\Delta\bu_{j+\hf}=\sum_k \widehat{\alpha}^k_{j+\hf}\widehat{{\mathbf R}}^k_{j+\hf}$ where
$\{\widehat{{\mathbf R}}^k_{j+\hf}\}_k$ is the  eigensystem of an intermediate Jacobian: for convenience,  $\widehat{{\mathbf R}}^k_{j+\hf}$ are computed as the eigensystem at the  arithmetic average, $A\big(\hf(\bu_j+\bu_{j+1})\big)$, rather than the Roe average,  $\Delta {\mathbf f}_{j+\hf}=A^{\textnormal{Roe}}_{j+\hf}\Delta \bu_{j+\hf}$, as it does not seem to affect the MUSCL results in this context.]

§.§ One-dimensional linear advection equation

We consider the one-dimensional linear advection equation

\begin{equation*}
\frac{\partial u}{\partial t}+\frac{\partial u}{\partial x} = 0, \qquad 0\leq x\leq 2\pi,
\end{equation*}

subject to initial data,

\[
= \sin^4(\pi x) & 0\leq x \leq 1, \\
\equiv 1 & 2.09 \leq x \leq 3.09, \\
= \textnormal{hat function with height} = 1 & 4.18\leq x \leq 5.18,\\
\equiv 0 & \textnormal{elsewhere in} \ (1,2\pi),
\end{array}\right.
\]

and  $2\pi$ periodic boundary conditions (BCs). We compute the solution after one-period of revolution using  $N = 400$ equi-spaced mesh grids. Figure <ref> displays the solutions computed with the NT scheme. Both the vA limiter and the minmod limiter generate non-oscillatory solutions with comparable resolution at the shocks and smooth extrema. Figure <ref> compares the solutions computed with different limiters as applied to the semi-discrete central-upwind scheme. In this case, applying the smoother vA limiter  improves the noticeable `clipping' of the minmod${}_1$ that occurred at the peaks of the sinusoidal wave and the triangular wave.

It is instructive to compare the NT results in Figure <ref> which are found to be superior to the corresponding results of the central-upwind KT scheme in Figure <ref>. In particular, a direct comparison in Figure <ref> shows that the `clipping'
phenomenon in the three peaks is clearly stronger in the KT computation, whereas the NT scheme produces sharper peaks. We recall that the NT with minmod limiter is reduced at peaks to a diffusive Lax-Friedrichs scheme. Thus, the vA limiter helps  improving the NT resolution at peaks.   Of course, the truncation error of NT schemes of order ${\mathcal O}((\Delta x)^2/\lambda)$ implies that the NT resolution decreases for small(-er) CFL. But for near optimal CFL=0.45, the resolution of NT  scheme with vA limiter offers a competitive alternative to central-upwind schemes. This conclusion will be confirmed with linear contact waves studied in section <ref> below.

Linear advection. NT scheme, $N = 400$, periodic BCs, $t = 2\pi$.

Linear advection. Semi-discrete central-upwind scheme, $N = 400$, periodic BCs, $t = 2\pi$.

Linear advection. NT vs. semi-discrete KT scheme with vA limiter with $N = 400$ mesh points,  periodic BCs at $t = 2\pi$.

§.§ One-dimensional Burgers' equation

We show the results for the one-dimensional Burgers' equation

\begin{equation*}
\frac{\partial u}{\partial t}+\frac{\partial}{\partial x}\Big(\frac{u^2}{2}\Big) = 0
\end{equation*}

over the domain $[0, 2\pi]$. In this test case, the solution comprises of an expansion wave and a compression wave. Figures <ref>,  <ref> present the solutions computed with $N = 200$ mesh grids at $t = 2$. In both fully discrete and semi-discrete test cases, the solutions generated with the van Albada limiter capture the expansion wave and the shock wave with sharp resolution and without generating numerical oscillations. The results based on vA limiter are essentially comparable to those based on the minmod limiter.

1D Burger's equation, NT scheme, $N = 200$, $t = 2$.

1D Burgers' equation, semi-discrete central-upwind scheme, $N = 200$, $t = 2$.

§.§ One-dimensional Euler equations

We consider the one-dimensional Euler equations of gas dynamics

\begin{equation*}
\frac{\partial}{\partial t} \left[\begin{array}{c}
\rho  \\
\rho u\\
\end{array}\right]
+\frac{\partial}{\partial x}\left[\begin{array}{c}
\rho u  \\
\rho u^2+p\\
\end{array}\right] = 0,
\end{equation*}

where $\rho$, $u$, $p$ are the density, velocity and pressure, $E$ is the total energy per unit volume given by

\[
E = \frac{p}{\gamma-1}+\frac{\rho u^2}{2}, \quad \gamma = 1.4.
\]

As first test case,  we consider the Sod's shock tube problem: the computational domain is $[0, 1]$ with an interface at $x = 0.5$, and subject to initial states on the two sides of the interface  given by

\begin{equation*}
\left\{\begin{array}{c}
\begin{split}
& p_{L} = 1.0, \quad p_{R} = 0.1  \\
& \rho_{L} = 1.0, \quad \rho_{R} = 0.125 \\
& u_{L} = 0.0, \quad u_{R} = 0.0
\end{split}
\end{array}\right..
\end{equation*}

Figures <ref> and <ref> show the computed density at $t = 0.2$ with $N = 400$ mesh grids, computed with the NT and, respectively, semi-discrete central-upwind schemes. The numerical results, comparing the vA and minmod${}_1$ limiters show that the high resolution of the computation based on vA limiter  are on par with those based on the minmod limiter. In fact, we observe that the contact wave in Figure <ref> is even slightly better resolved by the NT+vA limiter compared with the central-upwind KNP+vA limiter in Figure <ref>. This agrees with the improved resolution of NT+vA limiter of linear advection discussed in remark <ref> and similar remark regarding the diffusion of the semi-discrete central-upwind KT schemes which goes back to <cit.>.

Sod's problem — density computed with NT scheme and vA limiter, $N = 400$, $t = 0.2$.

Sod's problem — density computed with semi-discrete central-upwind scheme and vA limiter, $N = 400$, $t = 0.2$.

Next, we turn to a second test case of the Osher-Shu problem. The computational domain is $[-5 , 5]$ with the interface placed at $x = -4$, subject to initial data

\begin{equation*}
\left\{\begin{array}{c}
\begin{split}
& p_{L} = 10.33333, \quad p_{R} = 1.0,  \\
& \rho_{L} = 3.857143, \quad \rho_{R} = 1+0.2\sin(5x), \\
& u_{L} = 2.6293690, \quad u_{R} = 0.0.
\end{split}
\end{array}\right.
\end{equation*}

In this problem, the solution of density consists of a discontinuity and a smooth harmonic wave. Figures <ref> and <ref> present the solutions of density given by the NT scheme and the semi-discrete central-upwind scheme computed with $N = 600$ mesh grids at $t = 1.8$. We observe  that the solutions computed with the vA limiter resolve the smooth extremum better than those computed with the minmod limiter; in particular, the solutions do not introduce spurious numerical oscillations at the location of the shock.

Osher-Shu problem — density, NT scheme, $N = 600$, $t = 1.8$.

Osher-Shu problem — density, semi-discrete central-upwind scheme, $N = 600$, $t = 1.8$.

§.§ Two-dimensional Euler equations

We examine the performance of the vA limiter in two-dimensional problems. We consider the double Mach reflection problem [58] for the two-dimensional Euler equations,

\begin{equation*}
\frac{\partial}{\partial t}\left[\begin{array}{c}
\rho  \\
\rho u\\
\rho v\\
\end{array}\right]+
\frac{\partial}{\partial x}\left[\begin{array}{c}
\rho u  \\
\rho u^2+p\\
\rho uv\\
\end{array}\right]+
\frac{\partial}{\partial y}\left[\begin{array}{c}
\rho v  \\
\rho uv\\
\rho v^2+p\\
\end{array}\right] = 0, \quad p = (\gamma-1)\cdot[E-\frac{\rho}{2}(u^2+v^2)].
\end{equation*}

The problem is initiated with a Mach 10 oblique shock positioned at $(1/6,0)$ which makes a $60^{\circ}$ angel with $x-$axis. The computational domain is set to $[0,4]\times[0,1]$. The bottom boundary ($y=0$) consists of a reflecting wall beginning at $x = 1/6$, the short region from $x = 0$ to $x = 1/6$ is imposed with the initial post-shock conditions. The left boundary $(x = 0)$ is also assigned with the initial post-shock values, and at the right boundary $(x=4)$ all the gradients are set to zero. The values along the
top boundary $(y = 1)$ are set to describe the exact motion of the Mach 10 shock. We refer to [58] for more detailed descriptions. We compute the solution of density at $t = 0.2$ with $480\times 120$ mesh grids. Figure <ref> shows the reference solution given by the 5th order WENO finite difference method [45].

Figures <ref>, <ref> display the solutions computed with the NT scheme. It is clear from the results that the solution generated by the vA limiter has obviously better quality. The shock waves and the slip lines underneath are characterized by much sharper transitions. The computational results of the semi-discrete scheme are presented in Figures <ref>, <ref>. Still, the numerical scheme captures more information of the near-wall flow when used in conjunction with the vA limiter. These results confirm the advantage of the vA limiter (<ref>) in terms of resolving multidimensional complex flow structures.

Double Mach reflection-density, 5th order WENO, $N_x = 480$, $N_y = 120$, $t = 0.2$.

Double Mach reflection-density, NT scheme, minmod limiter, $N_x = 480$, $N_y = 120$, $t = 0.2$.

Double Mach reflection-density, NT scheme, vA limiter, $N_x = 480$, $N_y = 120$, $t = 0.2$.

Double Mach reflection-density, semi-discrete central-upwind scheme, minmod limiter, $N_x = 480$, $N_y = 120$, $t = 0.2$.

Double Mach reflection-density, semi-discrete central-upwind scheme, van Albada limiter, $N_x = 480$, $N_y = 120$, $t = 0.2$.

§.§ Incompressible Euler equations

We consider the two-dimensional double shear layer problem for incompressible Euler equations,

\[
\left\{\begin{array}{l}
\partial_t \bu+\nabla \cdot (\bu\otimes\bu) = -\nabla p,  \\ \\
\nabla\cdot\bu = 0,
\end{array}\right.
\]

where $p$ denotes the pressure and $\bu = (u,v)$ is the velocity field. The initial data is given by

\[
u(x,y,0) = \left\{\begin{array}{ll}
\tanh{((y-\frac{\pi}{2})/\rho)} & y\leq\pi, \\ \\
\tanh{((\frac{3\pi}{2}-y)/\rho)}  & y>\pi,
\end{array}\right. \qquad
v(x,y,0) = \delta \sin{(x)}
\]

over the rectangular domain $[0,2\pi]\times [0,2\pi]$. We take $\displaystyle \rho = \frac{\pi}{15}$ and $\delta = 0.05$. The equations are solved with the finite difference projection method (see e.g. [31, 21, 60]), which consists of a predictor step and a corrector step

\[
\begin{split}
&\frac{\bu^*-\bu^n}{\Dt}+\nabla\cdot(\bu^n\otimes\bu^n) = 0,\\
& \frac{\bu^{n+1}-\bu^*}{\Dt}+\nabla p^{n+1} = 0 \quad \text{with}\quad  \nabla\cdot\bu^{n+1} = 0.
\end{split}
\]

The divergence free constraint, $\nabla\cdot\bu^{n+1} = 0$, leads to a Poisson equation for pressure

\[
-\Delta p^{n+1} = -\frac{1}{\Dt}\nabla\cdot\bu^*.
\]

In the predictor step, we reconstruct the convection fluxes, $\bu^n\otimes \bu^n$, with the Rusanov-type flux (<ref>) coupled with the second order MUSCL reconstructions (<ref>). In the corrector step, assuming the uniform mesh size is employed, $\Dx = \Dy = h$, the gradient, divergence and Laplace operators are approximated with the second order central difference operators,

\[
\nabla^h \phi_{i,j}= \left(\begin{array}{c}
\displaystyle \frac{\phi_{i+1,j}-\phi_{i-1,j}}{2\Dx}  \\
\displaystyle \frac{\phi_{i,j+1}-\phi_{i,j+1}}{2\Dy}
\end{array}\right),
\]

\[
\nabla^h\cdot\bu_{i,j} =  \frac{u_{i+1,j}-u_{i-1,j}}{2\Dx}+\frac{v_{i+1,j}-v_{i-1,j}}{2\Dy}, \quad \bu = (u,v),
\]

\[
\Delta^h \phi_{i,j} = \frac{\phi_{i+2,j}-2\phi_{i,j}+\phi_{i-2,j}}{4\Dx^2}+ \frac{\phi_{i,j+2}-2\phi_{i,j}+\phi_{i,j-2}}{4\Dy^2}.
\]

To achieve a higher order of accuracy in time, we apply the explicit SSPRK3 method for time marching. The time step is taken according to the CFL condition

\[
\Dt\max_{i,j}\Big(\frac{|u_{i,j}|}{\Dx}+\frac{|v_{i,j}|}{\Dy}\Big)\leq 0.5.
\]

Figure <ref> presents the solutions of vorticity at $t = 8$. We observe that the solution contours given by the vA limiter are more concentrated than those by the $\minmod$ limiter. Figure <ref> compares the solutions of $v$ along the line $x=\pi$. The vA limiter shows clear advantage over the $\minmod$ limiter in terms of capturing sharp extrema.

5th order WENO

MUSCL-van Albada

Double shear layer, contours of vorticity. $N_x = N_y = 128$, $CFL = 0.5$, $t = 8$.

Double shear layer, solutions of $v$ along $x = \pi$, $N_x = N_y = 128$, $\CFL = 0.5$, $t = 8$.

§.§ Concluding remarks

The vA limiter dates back to the 1982 work [1]. Its TVD stability and related non-oscillatory properties were overshadowed by extensive studies of the class of minmod limiters, [26, 14, 27, 49, 5, 39, 36]. In this work, we revisit the vA limiter in the context of high-resolution central schemes, showing that its non-oscillatory properties   improve   corresponding computations using  the minmod limiter, or  are at least  on par with the minmod-based results in other computations.
This was demonstrated in the 1D Osher-Shu problem depicted in Figures <ref>–<ref> and more noticeably, in the 2D problems reported in sections <ref> & <ref>.
We do not claim that the use of vA limiter always outperforms the class of minmod limiters: the vA limiter are found on par with the minmod results for 1D Sod problem reported in Figures <ref>–<ref>; indeed, this is further highlighted in Figure <ref> below, where we zoom on the density variation of the Sod problem and observe comparable results for  the vA limiter which is compared with the less diffusive minmod${}_\theta$ with $\theta=\hf(1+\sqrt{2})$ (corresponding to the upper bound on the right of (<ref>)). The main point here is to shed light on the smoother vA limiter, showing that it produces equal or even better results than the class of minmod limiters. The vA limiter was  added to the minmod limiter  in the Python code [59] (solver scheme `sd3').

Sod's problem. Zoom on density computed by NT scheme, $N = 400$ at $t = 0.2$ using vA limter  vs. minmod${}_\theta$ limiter with $\theta=1.2$.

Statement. On behalf of all authors, the corresponding author states  that this work is in compliance with Ethical Standards, that we did not participate in any Research involving Humans/and or Animals under Compliance with Ethical standards, and  that there is no conflict of interest.

[1]
Gale Dick van Albada, Bram van Leer, and WW Roberts Jr.
A comparative study of computational methods in cosmic gas dynamics.
Astronomy and Astrophysics, vol. 108, no. 1, Apr. 1982, p.
76-84., 108:76–84, 1982.

[2]
P Arminjon, A St-Cyr, and A Madrane.
New two-and three-dimensional non-oscillatory central finite volume
methods on staggered cartesian grids.
Applied numerical mathematics, 40(3):367–390, 2002.

[3]
Jorge Balbás and Eitan Tadmor.
Central Station – a collection of references on high-resolution
non-oscillatory central schemes.
<https://www.math.umd.edu/ tadmor/centpack/publications/>, 2006.

[4]
Jorge Balbás, Eitan Tadmor, and Cheng-Chin Wu.
Non-oscillatory central schemes for one-and two-dimensional MHD
equations: I.
Journal of Computational Physics, 201(1):261–285, 2004.

[5]
S Chakravarthy and Stanly Osher.
A new class of high accuracy TVD schemes for hyperbolic
In 23rd Aerospace Sciences Meeting, page 363, 1985.

[6]
Bernd Einfeldt.
On Godunov-type methods for gas dynamics.
SIAM Journal on numerical analysis, 25(2):294–318, 1988.

[7]
Björn Engquist and Stanley Osher.
One-sided difference approximations for nonlinear conservation laws.
Mathematics of Computation, 36(154):321–351, 1981.

[8]
Ulrik S Fjordholm, Siddhartha Mishra, and Eitan Tadmor.
Arbitrarily high-order accurate entropy stable essentially
nonoscillatory schemes for systems of conservation laws.
SIAM Journal on Numerical Analysis, 50(2):544–573, 2012.

[9]
Ulrik S Fjordholm, Siddhartha Mishra, and Eitan Tadmor.
ENO reconstruction and ENO interpolation are stable.
Foundations of Computational Mathematics, 13:139–159, 2013.

[10]
Ulrik S Fjordholm and Deep Ray.
A sign preserving WENO reconstruction method.
Journal of Scientific Computing, 68:42–63, 2016.

[11]
Sergei K Godunov.
A finite difference method for the numerical computation of
discontinuous solutions of the equations of fluid dynamics.
Mat. Sb, 47(271-290):134, 1959.

[12]
Sigal Gottlieb, David I Ketcheson, and Chi-Wang Shu.
Strong stability preserving Runge-Kutta and multistep time
World Scientific, 2011.

[13]
Sigal Gottlieb, Chi-Wang Shu, and Eitan Tadmor.
Strong stability-preserving high-order time discretization methods.
SIAM review, 43(1):89–112, 2001.

[14]
Ami Harten.
High resolution schemes for hyperbolic conservation laws.
Journal of computational physics, 49(3):357–393, 1983.

[15]
Ami Harten, Bjorn Engquist, Stanley Osher, and Sukumar R Chakravarthy.
Uniformly high order accurate essentially non-oscillatory schemes,
Journal of computational physics, 71(2):231–303, 1987.

[16]
Ami Harten and Stanley Osher.
Uniformly high-order accurate nonoscillatory schemes. I.
In Upwind and High-Resolution Schemes, pages 187–217.
Springer, 1997.

[17]
Amiram Harten, James M Hyman, Peter D Lax, and Barbara Keyfitz.
On finite-difference approximations and entropy conditions for
Communications on pure and applied mathematics, 29(3):297–322,

[18]
Amiram Harten, Peter D Lax, and Bram van Leer.
On upstream differencing and Godunov-type schemes for hyperbolic
conservation laws.
SIAM review, 25(1):35–61, 1983.

[19]
Guang-Shan Jiang and Eitan Tadmor.
Nonoscillatory central schemes for multidimensional hyperbolic
conservation laws.
SIAM Journal on Scientific Computing, 19(6):1892–1917, 1998.

[20]
Claes Johnson, Anders Szepessy, and Peter Hansbo.
On the convergence of shock-capturing streamline diffusion finite
element methods for hyperbolic conservation laws.
Mathematics of computation, 54(189):107–129, 1990.

[21]
Raz Kupferman and Eitan Tadmor.
A fast, high resolution, second-order central scheme for
incompressible flows.
Proceedings of the National Academy of Sciences,
94(10):4848–4852, 1997.

[22]
Alexander Kurganov and Doron Levy.
A third-order semidiscrete central scheme for conservation laws and
convection-diffusion equations.
SIAM Journal on Scientific Computing, 22(4):1461–1488, 2000.

[23]
Alexander Kurganov, Sebastian Noelle, and Guergana Petrova.
Semidiscrete central-upwind schemes for hyperbolic conservation laws
and hamilton–jacobi equations.
SIAM Journal on Scientific Computing, 23(3):707–740, 2001.

[24]
Alexander Kurganov and Eitan Tadmor.
New high-resolution central schemes for nonlinear conservation laws
and convection–diffusion equations.
Journal of Computational Physics, 160(1):241–282, 2000.

[25]
Peter D. Lax.
Weak solutions of nonlinear hyperbolic equations and their numerical
Communications on Pure and Applied Mathematics, 7:159–193,

[26]
Bram van Leer.
Towards the ultimate conservative difference scheme. v. a
second-order sequel to Godunov's method.
Journal of computational Physics, 32(1):101–136, 1979.

[27]
Bram van Leer.
On the relation between the upwind-differencing schemes of Godunov,
Engquist–Osher and Roe.
SIAM Journal on Scientific and statistical Computing,
5(1):1–20, 1984.

[28]
Doron Levy, Gabriella Puppo, and Giovanni Russo.
Central WENO schemes for hyperbolic systems of conservation laws.
ESAIM: Mathematical Modelling and Numerical Analysis,
33(3):547–571, 1999.

[29]
Doron Levy, Gabriella Puppo, and Giovanni Russo.
A third order central WENO scheme for 2D conservation laws.
Applied Numerical Mathematics, 33(1-4):415–421, 2000.

[30]
Doron Levy, Gabriella Puppo, and Giovanni Russo.
A fourth-order central WENO scheme for multidimensional hyperbolic
systems of conservation laws.
SIAM Journal on scientific computing, 24(2):480–506, 2002.

[31]
Doron Levy and Eitan Tadmor.
Non-oscillatory central schemes for the incompressible 2-d Euler
Mathematical Research Letters, 4(3):321–340, 1997.

[32]
Xu-Dong Liu, Stanley Osher, and Tony Chan.
Weighted essentially non-oscillatory schemes.
Journal of computational physics, 115(1):200–212, 1994.

[33]
Xu-Dong Liu and Eitan Tadmor.
Third order nonoscillatory central scheme for hyperbolic conservation
Numerische mathematik, 79(3):397–425, 1998.

[34]
Yingjie Liu, Chi-Wang Shu, Eitan Tadmor, and Mengping Zhang.
L2 stability analysis of the central discontinuous Galerkin method
and a comparison between the central and regular discontinuous Galerkin
ESAIM: Mathematical Modelling and Numerical Analysis,
42(4):593–607, 2008.

[35]
Wim A Mulder and Bram van Leer.
Experiments with implicit upwind methods for the Euler equations.
Journal of Computational Physics, 59(2):232–246, 1985.

[36]
Haim Nessyahu and Eitan Tadmor.
Non-oscillatory central differencing for hyperbolic conservation
Journal of computational physics, 87(2):408–463, 1990.

[37]
Stanley Osher.
Riemann solvers, the entropy condition, and difference.
SIAM Journal on Numerical Analysis, 21(2):217–235, 1984.

[38]
Stanley Osher.
Convergence of generalized MUSCL schemes.
SIAM Journal on Numerical Analysis, 22(5):947–961, 1985.

[39]
Stanley Osher and Eitan Tadmor.
On the convergence of difference approximations to scalar
conservation laws.
Mathematics of Computation, 50(181):19–51, 1988.

[40]
Serge Piperno and Sophie Depeyre.
Criteria for the design of limiters yielding efficient high
resolution tvd schemes.
Computers & fluids, 27(2):183–197, 1998.

[41]
Philip L Roe.
Approximate riemann solvers, parameter vectors, and difference
Journal of computational physics, 43(2):357–372, 1981.

[42]
Richard Saurel and Rémi Abgrall.
A simple method for compressible multifluid flows.
SIAM Journal on Scientific Computing, 21(3):1115–1145, 1999.

[43]
Chi-Wang Shu.
TVB uniformly high-order schemes for conservation laws.
Mathematics of Computation, 49(179):105–121, 1987.

[44]
Chi-Wang Shu.
Total-variation-diminishing time discretizations.
SIAM Journal on Scientific and Statistical Computing,
9(6):1073–1084, 1988.

[45]
Chi-Wang Shu.
Essentially non-oscillatory and weighted essentially non-oscillatory
schemes for hyperbolic conservation laws.
In Advanced numerical approximation of nonlinear hyperbolic
equations ,
pages 285–372. Springer, 1998.

[46]
Chi-Wang Shu.
Essentially non-oscillatory and weighted essentially non-oscillatory
Acta Numerica, 29:701–762, 2020.

[47]
Michael L Shur, Philippe R Spalart, and Michael Kh Strelets.
Noise prediction for increasingly complex jets. part i: Methods and
International journal of aeroacoustics, 4(3):213–245, 2005.

[48]
Raymond J Spiteri and Steven J Ruuth.
A new class of optimal high-order strong-stability-preserving time
discretization methods.
SIAM Journal on Numerical Analysis, 40(2):469–491, 2002.

[49]
Peter K Sweby.
High resolution schemes using flux limiters for hyperbolic
conservation laws.
SIAM journal on numerical analysis, 21(5):995–1011, 1984.

[50]
Eitan Tadmor.
Numerical viscosity and the entropy condition for conservative
difference schemes.
Mathematics of Computation, 43(168):369–381, 1984.

[51]
Eitan Tadmor.
Convenient total variation diminishing conditions for nonlinear
difference schemes.
SIAM journal on numerical analysis, 25(5):1002–1014, 1988.

[52]
Eitan Tadmor.
Total variation and error estimates for spectral viscosity
MATHEMATICS of computation, 60(201):245–256, 1993.

[53]
Eitan Tadmor.
Approximate solutions of nonlinear conservation laws and related
In Advanced numerical approximation of nonlinear hyperbolic
equations ,
pages 1–149. Springer, 1998.

[54]
Eitan Tadmor.
Entropy stability theory for difference approximations of nonlinear
conservation laws and related time-dependent problems.
Acta Numerica, 12:451–512, 2003.

[55]
Eitan Tadmor.
Selected topics in approximate solutions of nonlinear conservation
laws. high-resolution central schemes.
Nonlinear Conservation Laws and Applications, Springer, New
York, pages 101–122, 2011.

[56]
Eleuterio F Toro, Michael Spruce, and William Speares.
Restoration of the contact surface in the hll-riemann solver.
Shock waves, 4(1):25–34, 1994.

[57]
Sulin Wang and Zhengfu Xu.
Total variation bounded flux limiters for high order finite
difference schemes solving one-dimensional scalar conservation laws.
Mathematics of Computation, 88(316):691–716, 2019.

[58]
Paul Woodward and Phillip Colella.
The numerical simulation of two-dimensional fluid flow with strong
Journal of computational physics, 54(1):115–173, 1984.

[59]
Anil Zenginoglu.
centpy: Central schemes for conservation laws in python.

[60]
Ju Zhang and Thomas L Jackson.
A high-order incompressible flow solver with WENO.
Journal of Computational Physics, 228(7):2426–2442, 2009.


# On Projective Planes of Order 16 Associated with 1-rotational 2-$(52,4,1)$
Designs

Zazil Santizo Huerta, Vladimir Tonchev, and Melissa Keranen
Michigan Technological University
Department of Mathematical Sciences
Houghton, MI 49931, USA

###### Abstract

A maximal arc of degree $k$ in a finite projective plane $\cal{P}$ of order
$q=ks$ is a set of $(q-s+1)k$ points that meets every line of $\cal{P}$ in
either $k$ or $0$ points. The collection of the nonempty intersections of a
maximal arc with the lines of $\cal{P}$ is a resolvable Steiner
2-$((q-s+1)k,k,1)$ design. Necessary and sufficient conditions for a
resolvable Steiner 2-design to be embeddable as a maximal arc in a projective
plane were proved recently in [8]. Steiner designs associated with maximal
arcs in the known projective planes of order 16 were analyzed in [6], where it
was shown that some of the associated designs are embeddable in two
nonisomorphic planes.

Using MAGMA, we conducted an analysis to ascertain whether any of the 22 non-
isomorphic 1-rotational 2-$(52,4,1)$ designs, previously classified in [3],
could be embedded in maximal arcs of degree 4 within projective planes of
order 16. This paper presents a summary of our findings, revealing that
precisely only one out of the the twenty-two 1-rotational designs from [3] is
embeddable in a plane of order 16, being the Desarguesian plane $PG(2,16)$.

##  1 Introduction

We assume familiarity with basic facts and notions from combinatorial design
theory, and emphasize the following definitions needed to state and prove our
main theorem.

A 2-$(v,k,\lambda)$ design $D=\\{X,\cal{B}\\}$ is a pair of a point set $X$ of
size $v$ and a collection of $k$-subsets $\cal{B}$=$\\{B_{j}\\}_{j=1}^{b}$
called blocks, such that every pair of points is contained in exactly
$\lambda$ blocks. Every point of $X$ is contained in $r=\lambda(v-1)/(k-1)$
blocks, and the number of blocks is $b={\frac{v(v-1)}{k(k-1)}}\lambda$. A
design with $v=b$ is called symmetric, and a Steiner design is a design with
$\lambda=1$.

Let $D=\\{X,B\\}$ be a design where $X=\\{x_{1},\ldots,x_{v}\\}$ and
$B=\\{B_{1},\ldots,B_{b}\\}$. The incidence matrix of $D$ is the $v\times b$
0-1 matrix $M=(m_{i,j})$ defined by the rule

$m_{i,j}=\begin{cases}1&\hbox{if}\ x_{i}\in A_{j}\\\ 0&\hbox{if}\ x_{i}\notin
A_{j}\end{cases}$

The incidence matrix, $M$, of a $(v,b,r,k,\lambda)$-BIBD satisfies the
following properties:

  1. 1.

every column of $M$ contains exactly $k$ “1”s;

  2. 2.

every row of $M$ contains exactly $r$ “1”s;

  3. 3.

two distinct rows of $M$ both contain “1”s in exactly $\lambda$ columns.

Two designs are isomorphic if there is a bijection between their point sets
that maps every block of the first design to a block of the second design. An
automorphism of a design is any isomorphism of the design to itself. The set
of all automorphisms of $D$ form the automorphism group $Aut(D)$ of $D$.

A parallel class of a 2-$(v,k,\lambda)$ design $D$ with $v=nk$ is a set of
$n=v/k$ pairwise disjoint blocks. A resolution of $D$ is a partition of the
collection of blocks $\cal{B}$ into $r=\lambda(v-1)/(k-1)$ disjoint parallel
classes. A design is resolvable if it admits at least one resolution.

A finite projective plane of order $q$ is defined as a set of $q^{2}+q+1$
points, a set of lines, and a relation between points and lines, having the
following properties:

  * •

Any two points determine a line,

  * •

Any two lines determine a point,

  * •

Every point has $q+1$ lines through it, and

  * •

Every line contains $q+1$ points.

A projective plane of order $q\geq 2$ is a symmetric Steiner
2-$(q^{2}+q+1,q+1,1)$ design. An affine plane of order $q\geq 2$ is a Steiner
2-$(q^{2},q,1)$ design.

A maximal $(m,k)$-arc in a projective plane $\cal{P}$ of order $q=ks$ is a set
of $m=(q-s+1)k$ points that meets every line of $\cal{P}$ in either $k$ or $0$
points. The set of lines of $\cal{P}$ which have no points in a given maximal
$((sk-s+1)k,k)$-arc $\cal{A}$ determines a dual maximal $((sk-k+1)s,s)$-arc
${\cal{A}}^{\perp}$ in the dual plane ${\cal{P}}^{\perp}$.

The classical (or Desarguesian) plane $PG(2,p^{t})$ of order $n=p^{t}$, where
$p$ is prime and $t\geq 1$, has as points the 1-dimensional subspaces of the
3-dimensional vector space $V^{3}$ over the finite field of order $p^{t}$, and
as blocks (or lines), the 2-dimensional subspaces of $V^{3}$.

In Section 2, we provide a framework for understanding the relationships
between maximal arcs and Steiner 2-designs, resolutions, and their properties
within projective planes. It establishes theoretical groundwork and presents
specific findings within the context of planes of order 16.

In Section 3, we give the construction and analyze the 22 non isomorphic
resolvable 1-rotational 2-$(52,4,1)$ designs to determine whether it is
possible to embed them as maximal arcs of degree 4 within projective planes of
order 16.

In section 4, we provide a clear outline of the computational approach
employed to determine the embeddability of the given 1-rotational 2-$(52,4,1)$
designs into projective planes of order 16. The analysis of our computations
culminates in the presentation of the main theorem, which identifies exactly
one design is embeddable.

The following theorems concern the existence of maximal arcs within
Desarguesian planes.

######  Theorem 1.

[5] The Desarguesian plane $PG(2,2^{m})$ contains maximal arcs for every
$k=2^{i}$, $1\leq i<m$.

######  Theorem 2.

[2] Maximal $k$-arcs with $1<k<q$ do not exist in any Desarguesian plane,
$PG(2,q)$, if $q$ is odd.

##  2 Previous Results

Let $\cal{P}$ be a projective plane of order $q=sk$, $1<k<q$.

  * •

The nonempty intersections of a maximal $((sk-s+1)k,k)$-arc $\cal{A}$ with the
lines of $\cal{P}$ form a Steiner 2-$((sk-s+1)k,k,1)$ design $D$.

  * •

Respectively, the dual $((sk-s+1)k,k)$-arc $\cal{A}^{\perp}$ is the point set
of a Steiner 2-$((sk-k+1)s,s,1)$ design $D^{\perp}$.

######  Theorem 3.

  * •

The points of ${\cal{A}}^{\perp}$ determine a set of $(sk-k+1)s$ resolutions
of $D$.

  * •

The points of ${\cal{A}}$ determine a set of $(sk-s+1)k$ resolutions of the
design $D^{\perp}$ associated with the dual arc ${\cal{A}}^{\perp}$.

Every two resolutions of $D$ (resp. $D^{\perp}$) share one parallel class.

Tonchev gave the following definition in [6], which is motivated by properties
of Steiner designs associated with maximal arcs.

Two resolutions, $R_{1}$, $R_{2}$ of a 2-$((sk-s+1)k,k,1)$ design,

$R_{1}=P^{(1)}_{1}\cup P^{(1)}_{2}\cup\cdots P^{(1)}_{r},\
R_{2}=P^{(2)}_{1}\cup P^{(2)}_{2}\cup\cdots P^{(2)}_{r}$

are compatible if they share one parallel class, $P^{(1)}_{i}=P^{(2)}_{j}$,
and $|P^{(1)}_{i^{\prime}}\cap P^{(2)}_{j^{\prime}}|\leq 1$ for
$(i^{\prime},j^{\prime})\neq(i,j)$. More generally, a set of $m$ resolutions
$R_{1},\ldots,R_{m}$ is compatible if every two of these resolutions are
compatible.

A Steiner 2-$((sk-s+1)k,k,1)$ design associated with a maximal
$((sk-s+1)k,k)$-arc $\cal{A}$ in a plane $\cal{P}$ of order $q=ks$ admits a
set of $(sk-k+1)s$ pairwise compatible resolutions.

The following result tells us when a design is embeddable in a projective
plane as a maximal arc.

######  Theorem 4.

[8] Let $S=\\{R_{1},\ldots,R_{m}\\}$ be a set of $m$ pairwise compatible
resolutions of a 2-$((sk-s+1)k,k,1)$ design $D=\\{X,\cal{B}\\}$. Then

$m\leq(sk-k+1)s.$

The equality $m=(sk-k+1)s$ holds if and only if there exists a projective
plane $\cal{P}$ of order $q=sk$ such that $D$ is embeddable in $\cal{P}$ as a
maximal $((sk-s+1)k,k)$-arc.

Hamilton, et. al. have determined which of the known projective planes of
order 16 contain maximal $(52,4)$-arcs.

######  Theorem 5.

[7] All but possibly four of the 22 known projective planes of order 16,
(BBH2, BBS4, and their duals), contain maximal $(52,4)$-arcs.

The following theorem states that in the projective geometry over the Galois
field $GF(2)$ with 16 elements, there are exactly two distinct sets of
$(52,4)$-arcs. These arcs are projectively inequivalent, meaning they cannot
be transformed into each other through projective transformations while
preserving the structure. Furthermore, both of these arcs belong to a specific
type known as Denniston type.

######  Theorem 6.

[1] $PG(2,16)$ contains exactly two projectively inequivalent maximal
$(52,4)$-arcs, both of Denniston type.

In the following theorem, the listed planes contain $(52,4)$ maximal arcs,
which have associated 2-$(52,4,1)$ designs. Importantly, each design admits
two different sets of 52 compatible resolutions, leading to two non-isomorphic
planes of order 16.

######  Theorem 7.

[6] The planes $LMRH$, JOHN, BBH1, and JOWK, contain maximal $(52,4)$ arcs
whose associated 2-$(52,4,1)$ designs admit two different  sets of 52
compatible resolutions that lead to two nonisomorphic planes of order 16.

We now study a different question. Instead of asking which projective planes
contain a maximal arc, we would like to answer whether a given design can be
embedded as a maximal arc in a projective plane.

##  3 1-Rotational 2-$(52,4,1)$ Designs

A 2-$(v,k,\lambda)$ design is 1-rotational if it admits an automorphism of
order $v-1$.

Buratti and Zuanni [4] have shown that a 1-rotational $(52,4,1)$-RBIBD is
equivalent to a 4-set $F=\\{A_{1},A_{2},A_{3},A_{4}\\}$ of 4-subsets of
$\mathbb{Z}_{51}$ satisfying that any element of $\mathbb{Z}_{51}$ \-
$\\{0,17,34\\}$ is representable as a difference of two elements of a same
$A_{i}$, and $A_{1}\cup A_{2}\cup A_{3}\cup
A_{4}=\mathbb{Z}_{17}-\\{0\\}\pmod{17}$. They also provided a construction of
the 22 non-isomorphic 1-rotational 2-$(52,4,1)$ designs in [3].

######  Theorem 8.

[3] Up to isomorphism, there exist exactly 22 nonisomorphic resolvable
1-rotational 2-$(52,4,1)$ designs.

Given a $(51,3,4,1)$ resolvable difference family
$F=\\{A_{1},A_{2},A_{3},A_{4}\\}$. We give the associated 2-$(52,4,1)$-RBIBD.
This design has the point set $\mathbb{Z}_{51}\cup\\{\infty\\}$ and parallel
classes obtained by developing $\pmod{51}$ the starter parallel class given by

$\\{0,17,34,\infty\\}\cup\\{A_{i}+17j|1\leq i\leq 4;0\leq j<3\\}$

A complete set of representatives $F_{i}$, $i=1,\ldots,22$, for each of the 22
1-rotational designs are listed in [3].

We would like to determine which of the Buratti-Zuanni 2-$(52,4,1)$ designs
are embeddable as maximal $(52,4)$-arcs in a projective plane of order 16.

##  4 Main Theorem and Results

We wish to determine which of these designs can be embedded into a projective
plane of order 16 as a maximal $(54,4)$\- arc. Thus, to answer our question,
we need that the number of compatible resolutions in the design will be 52, by
Theorem 4.

Finding the compatible sets of maximum size, each containing 52 resolutions,
is essentially equivalent to identifying subsets of resolutions where every
pair is compatible. We now present some fundamental definitions of graph
theory, which are essential for establishing the graphical theoretical
equivalent of this problem.

A graph is a pair $G=(V,E)$ (sometimes called an undirected graph) consists of
a set $V$ of vertices (also called nodes) and a set $E$ of edges. If two
vertices in a graph are connected by an edge, we say the vertices are
adjacent. A clique is a subset of vertices of an undirected graph such that
every two distinct vertices in the clique are adjacent. By constructing a
graph $P$ where each vertex represents a resolution of $F_{i}$, and two
vertices are adjacent if the corresponding resolutions are compatible, we
effectively represent the compatibility relationships between resolutions as
edges in the graph. Since a clique in a graph is a subset of vertices where
each pair is adjacent, finding maximum-sized cliques in graph $P$ corresponds
to identifying sets of maximum size among resolutions that are compatible.
This is because each clique represents a compatible set of resolutions, and
its size determines the number of resolutions in the set. Therefore,
identifying the size of the largest cliques in $P$ concurrently reveals the
size of the largest set of compatible resolutions.

We now proceed to describe the outline of the algorithm created to compute
this number.

Algorithm Outline:

  * •

Construct the incidence matrices of the twenty-two 1-rotational designs from
[3].

  * •

Using Magma, compute all parallel classes and resolutions of each design
$F_{i}$.

  * •

For each $F_{i}$, define a graph $P$ whose vertices are the resolutions of
$F_{i}$, where two resolutions are adjacent if and only if they are
compatible.

  * •

Determine the largest size, $m$, of a clique in $P$.

  * •

If $m=52$, then by Theorem 4, $F_{i}$ is embbedable as a maximal $(52,4)$-arc
in the projective plane.

Table 1 summarizes the results of running the magma code provided in section
5. We see that $F_{17}$ is the only design that achieves the maximum clique
size of 52. This leads us to our main theorem.

######  Theorem 9.

Only one of the twenty-two 1-rotational 2-(52,4,1) designs, namely $F_{17}$,
is embeddable as a maximal (52,4)-arc in a plane of order 16, being the
Desarguesian plane $PG(2,16$).

Design[3] |  ${\bf|Aut(F_{i})|}$ |  Par. Cl. |  Res. |  Comp. Res.
---|---|---|---|---
$F_{1}$ |  51 |  34 |  1 |  1
$F_{2}$ |  51 |  17 |  1 |  1
$F_{3}$ |  51 |  17 |  1 |  1
$F_{4}$ |  51 |  17 |  1 |  1
$F_{5}$ |  51 |  34 |  1 |  1
$F_{6}$ |  51 |  17 |  1 |  1
$F_{7}$ |  51 |  34 |  1 |  1
$F_{8}$ |  51 |  17 |  1 |  1
$F_{9}$ |  51 |  34 |  1 |  1
$F_{10}$ |  51 |  34 |  1 |  1
$F_{11}$ |  51 |  17 |  1 |  1
$F_{12}$ |  204 |  442 |  18 |  1
$F_{13}$ |  204 |  289 |  18 |  1
$F_{14}$ |  204 |  153 |  1 |  1
$F_{15}$ |  102 |  306 |  1 |  1
$F_{16}$ |  408 |  51 |  1 |  1
$F_{17}$ |  408 |  2550 |  460 |  52
$F_{18}$ |  51 |  34 |  1 |  1
$F_{19}$ |  51 |  34 |  1 |  1
$F_{20}$ |  51 |  17 |  1 |  1
$F_{21}$ |  51 |  17 |  1 |  1
$F_{22}$ |  51 |  17 |  1 |  1
Table 1: The 22 Resolvable 1-rotational 2-$(52,4,1)$ designs

##  5 Magma Code

⬇

> DifferenceSet := func< S, n| [{(x + i) mod n + 1 : x in S } : i in [1 .. n]
> ] >;

> A := {18,33,22,46};

> B := {21,30,37,31};

> C := {6,45,25,43};

> D := {24,27,19,49};

> E:= { {1,18,35,52}, {2,19,36,52}, {3,20,37,52}, {4,21,38,52}, {5,22,39,52},
> {6,23,40,52}, {7,24,41,52}, {8,25,42,52}, {9,26,43,52}, {10,27,44,52},
> {11,28,45,52}, {12,29,46,\

52}, {13,30,47,52}, {14,31,48,52}, {15,32,49,52}, {16,33,50,52}, {17,34,51,52}
};

> F := Design< 2, 52| DifferenceSet(A, 51), DifferenceSet(B, 51),
> DifferenceSet(C, 51), DifferenceSet(D, 51), E >;

> pRank(F, 2);

41

> AutomorphismGroup(F);

Permutation group acting on a set of cardinality 52

Order = 408 = 2^3 * 3 * 17

(1, 50, 49, 23, 10, 29, 13, 5)(2, 25, 11, 4, 26, 37, 17, 7)(3, 51, 24, 36, 42,
45, 21, 9)(6, 27, 12, 30, 39, 18, 33, 15)(8, 28, 38, 43, 20, 34, 41, 19)(14,
31)(16, 32, 40,

44, 46, 47, 22, 35)

(2, 27, 14, 33, 17, 9, 5, 3)(4, 28, 40, 46, 49, 25, 13, 7)(6, 29, 15, 8, 30,
41, 21, 11)(10, 31, 16, 34, 43, 22, 37, 19)(12, 32, 42, 47, 24, 38, 45,
23)(18, 35)(20, 36, 44,

48, 50, 51, 26, 39)

> G :=AllParallelClasses(F);

> #G;

2550

> H := AllParallelisms(F);

> #H;

460

> V := { H[i]: i in [1 .. #H] };

> edgesSet := {};

> for i in [1 .. #H] do

for> for j in [1 .. #H] do

for|for> T := {};

for|for> for m in [1 .. 17] do

for|for|for> for n in [1 .. 17] do

for|for|for|for> if (i lt j) and (#(H[i] meet H[j]) eq 1) and
((#((SetToIndexedSet(H[i])[m]) meet (SetToIndexedSet(H[j])[n])) le 1)) then

for|for|for|for|if> T join:= {<m,n>};

for|for|for|for|if> end if;

for|for|for|for> end for;

for|for|for> end for;

for|for> if #T eq 288 then

for|for|if> edgesSet join:= {{i, j}};

for|for|if> end if;

for|for> end for;

for> end for;

> #edgesSet;

1326

> P := Graph< V| edgesSet >;

> M := MaximumClique(P);

> W := AllCliques(P, #M);

> #M;

52

> #W;

1

##  6 Conclusions and Future Work

Our goal is to validate our magma code by applying it to other
2-$(v,k,\lambda)$ 1-rotational designs that meet the criteria for embedding as
a maximal arc in a projective plane. Additionally, we aim to extend the
applicability of this code to test it on other types of Steiner 2-designs.

## References

  * [1] Ball, S. M., and Blokhuis, A. (2002). The classification of maximal arcs in small Desarguesian planes. Bulletin of the Belgian Mathematical Society : Simon Stevin, 9(3), 433-445.
  * [2] S. Ball, A. Blokhuis, F. Mazzocca, Maximal arcs in Desarguesian planes of odd order to not exist, Combinatorica 17 (1997) 31-41.
  * [3] M. Buratti and F. Zuanni, The 1-rotational $(52,4,1)$-RBIBD’s, JCMCC 30 (1999), 99 - 102.
  * [4] Buratti, Marco and Zuanni, Fulvio. (2000). Addendum to: ”$G$-invariantly resolvable Steiner 2-designs which are 1-rotational over $G$”. Bulletin of the Belgian Mathematical Society, Simon Stevin. 5. 10.36045/bbms/1103055695.
  * [5] R.H.F. Denniston, Some maximal arcs in finite projective planes, J. Combin. Theory 6 (1969) 317-319.
  * [6] M. Gezek, T. Wagner and V. D. Tonchev, Maximal arcs in projective planes of order 16 and related designs, Advances in Geometry, to appear.
  * [7] Gezek, Mustafa, ”COMBINATORIAL PROBLEMS RELATED TO CODES, DESIGNS AND FINITE GEOMETRIES”, Campus Access Dissertation, Michigan Technological University, 2017.
  * [8] V. D. Tonchev, On resolvable Steiner 2-designs and maximal arcs in projective planes, Designs, Codes and Cryptography, 84, July (2017), 165-172.