
\subsection{Optimization of the stability condition by choosing $\mu$}
\section{Proof for stability: Theorem \ref{thm:stab} }
\label{sec:stab:proofs:opt}

When $k=0$, the optimization is straightforward, and the detail is omitted. The remaining of this section will be devoted to the case when $k\geq 1$, for which the analysis is more technically involved.  
From here on, we assume $1\leq k\leq 9$. With this,  we have ${\mathcal K}>1$ and $\widehat{C}_{inv}>0$. We also assume $\omega>1/2$, though not all preliminary results next depend on this assumption. One can refer to Table \ref{tab:notation} for a summary of notation.



\subsection{Preliminary lemmas}

We first state and prove some preparatory lemmas.  Lemma \ref{lem:muS:mono} and Lemma \ref{lem:tau12} can be directly verified and the proofs are skipped. 


\begin{lem} 
\label{lem:muS:mono}
\begin{itemize}
\item[(i)]  With $\omega>1/2$, there always holds   $\mu_\star\in (\frac{1}{2\omega}, 1)$.
\item[(ii)]  With $\mu_S(\lambda)$ defined
 in \eqref{eq:thm:stab5:o:3}, let its inverse be  $\lambda_S(\mu): =2(\mu-\frac{1}{2\omega})\frac{C_{inv}||v||_\infty}{\widehat{C}_{inv}||v^2||_\infty}$. 
\begin{itemize}
\item Both $\mu_S(\lambda)$ and $\lambda_S(\mu)$ are monotonically increasing. And $\mu_S(\lambda)>\frac{1}{2\omega}, \forall \lambda>0$.
\item With $\widehat{\lambda}_\star=\lambda_S(1)$, we have $\mu_S(\widehat{\lambda}_\star) = 1$. In addition, $ \mu_S(\lambda)<1\Leftrightarrow \lambda <\widehat{\lambda}_\star$. 
\item  $\mu_S(\lambda_\star)= \mu_\star$ and $\lambda_S(\mu_\star)=\lambda_\star$.
\end{itemize}
\end{itemize}
\end{lem}


\begin{lem}\label{lem:mu_star}
Consider $\mu\in (\frac{1}{2\omega}, 1]$, then 
\begin{itemize}
\item[(i)] 
\begin{equation}
 \lambda_{1} (\mu)  \leq  \lambda_{2} (\mu)\Longleftrightarrow \mu\leq \mu_\star \left(\Longleftrightarrow \frac{1}{2\omega}<\mu\leq \mu_\star<1\right),
\end{equation}
and $ \lambda_{1} (\mu_\star) =  \lambda_{2} (\mu_\star)=\lambda_\star$. In addition, $ \lambda_{1} (\mu)$ is monotonically increasing on
$(\frac{1}{2\omega}, \mu_\star]$, and   $ \lambda_{2} (\mu)$ is monotonically decreasing.
\item[(ii)] 
\begin{equation}
\lambda_S(\mu) \leq  \lambda_{1} (\mu)\Longleftrightarrow \mu\leq \mu_\star \left(\Longleftrightarrow \frac{1}{2\omega}<\mu\leq \mu_\star<1\right).
\end{equation}
\item[(iii)] 
\begin{equation}
\widehat{\lambda}_\star> \lambda_{1} (\mu), \quad \widehat{\lambda}_\star>  \lambda_{2} (\mu), \quad \forall \mu\in (\frac{1}{2\omega}, 1].
\end{equation}
\end{itemize}
\end{lem}
\begin{proof} For $\mu\in (\frac{1}{2\omega}, 1]$,  to prove (i), 
\begin{align*}
 \lambda_{1} (\mu)  \leq  \lambda_{2} (\mu)&\Longleftrightarrow\sqrt{\frac{(1-\mu)(\mu-\frac{1}{2\omega}) }{2\widehat{C}_{inv} ||v^2||_\infty} } \leq \frac{1-\mu}{4C_{inv}||v||_\infty}\\
&\Longleftrightarrow
\frac{\mu-\frac{1}{2\omega} }{\widehat{C}_{inv}||v^2||_\infty} \leq \frac{1-\mu}{8(C_{inv}||v||_\infty)^2}
\Longleftrightarrow  \mu  \leq  \mu_\star.
\end{align*}
The equality is achieved at $\mu=\mu_\star$, with the value being  $\lambda_\star$.  The monotonicity of $ \lambda_{2} (\mu)$ is straightforward.
For $ \lambda_{1} (\mu)$, note that with ${\mathcal K}>1$, we have
 $\mu_\star<\frac{1}{2}\left(1+\frac{1}{2\omega}\right)$, with $\frac{1}{2}\left(1+\frac{1}{2\omega}\right)$ being  where  $ \lambda_{1} (\mu)$ achieves its maximum. This implies that $ \lambda_{1} (\mu)$, whose square is a downward-facing parabola,  is monotonically increasing on $(\frac{1}{2\omega}, \mu_\star]$.


To prove (ii),  we proceed as below.
\begin{align*}
\lambda_S(\mu) \leq  \lambda_{1} (\mu) \Longleftrightarrow &2(\mu-\frac{1}{2\omega})\frac{C_{inv}||v||_\infty}{\widehat{C}_{inv}||v^2||_\infty}
\leq 
\sqrt{\frac{(1-\mu)(\mu-\frac{1}{2\omega} ) } {2\widehat{C}_{inv} ||v^2||_\infty} }
\notag\\
\Longleftrightarrow &  (\mu-\frac{1}{2\omega})\frac{8(C_{inv}||v||_\infty)^2}{\widehat{C}_{inv}||v^2||_\infty} \leq 1-\mu 
\Longleftrightarrow  \mu \leq \mu_\star.
\end{align*}

To prove (iii), related to  $ \lambda_{2} (\mu)$, given its being  monotonically decreasing, we only need to show 
$\widehat{\lambda}_\star> \lambda_{2} (\frac{1}{2\omega})$, which is ensured by ${\mathcal K}>1$ as below.
\begin{equation}
\widehat{\lambda}_\star> \lambda_{2} (\frac{1}{2\omega})\Longleftrightarrow
2(1-\frac{1}{2\omega})\frac{C_{inv}||v||_\infty}{\widehat{C}_{inv}||v^2||_\infty}>\frac{1-\frac{1}{2\omega}}{4C_{inv}||v||_\infty}
\Longleftrightarrow {\mathcal K}>1.
\end{equation}
Related to $ \lambda_{1} (\mu)$, from the proof of (i) of this lemma, we only need to verify 
$\widehat{\lambda}_\star> \lambda_{1} (\mu)|_{\mu= \frac{1}{2}(1+\frac{1}{2\omega})}.$ This can be argued as follows.
\begin{equation}
\widehat{\lambda}_\star> \lambda_{1} (\mu)|_{\mu= \frac{1}{2}(1+\frac{1}{2\omega})}\Longleftrightarrow
2(1-\frac{1}{2\omega})\frac{C_{inv}||v||_\infty}{\widehat{C}_{inv}||v^2||_\infty}>\frac{1-\frac{1}{2\omega}}{2\sqrt{2\widehat{C}_{inv}||v^2||_\infty }}
\Longleftrightarrow 4{\mathcal K}>1.
\end{equation}
This holds due to that ${\mathcal K}>1$. 
\end{proof}








\begin{rem}\label{rem:geometry}
Lemmas  \ref{lem:muS:mono}-\ref{lem:mu_star} tell the properties  and the relative locations of the curves $\lambda=\lambda_S(\mu)$, $\lambda= \lambda_{1} (\mu)$ and $\lambda= \lambda_{2} (\mu)$. Particularly, 
\begin{itemize}
\item
According to Lemmas  \ref{lem:muS:mono}-\ref{lem:mu_star},   the curves $\lambda=\lambda_S(\mu)$, $\lambda= \lambda_{1} (\mu)$ and $\lambda= \lambda_{2} (\mu)$ intersect at $(\mu_\star,\lambda_\star)$.
\item 
According to Lemma \ref{lem:mu_star}, to the left of  $\mu=\mu_\star$, the graph of $\lambda= \lambda_{2} (\mu)$ is above  that of $\lambda= \lambda_{1} (\mu)$, which is above the graph of $\lambda=\lambda_S(\mu)$; to the right of $\mu=\mu_\star$,  the ordering is reversed. 
\end{itemize}
\end{rem}

It is important to know the relative locations of various curves to optimize the time step condition. For general weight function $\omega$, it is nontrivial to visualize these curves, yet their  relative locations and some special points are captured in Figure \ref{fig:proof}, which is for the constant weight function $\omega\equiv 1$.  The figure can also facilitate the readers to follow and understand the analysis in this section, which is given algebraically  for general $\omega$ and has a geometric interpretation  for the special case of $\omega\equiv 1$.




\begin{figure}[!h]
  \centering
  {\includegraphics[width=0.5\textwidth]{plot4proof/proof.png}}
  \caption{ Plots with constant $\omega\equiv1$ to facilitate the understanding of  Lemmas \ref{lem:muS:mono}-\ref{lem:mu_star}. The scheme is: i) unconditionally stable when $\lambda=\varepsilon/(\sigma_m h)$ and $\mu$ fall into the gray region, ii) $\mu$-stable under $\Delta t\leq \tau_{\varepsilon,h,1}(\mu)$ in the stripped region,  iii) $\mu$-stable under $\Delta t\leq \tau_{\varepsilon,h,2}(\mu)$ in the latticed region, and iv) $\mu$-stable under $\Delta t\leq \min(\tau_{\varepsilon,h,1}(\mu), \tau_{\varepsilon,h,2}(\mu))$ in the  blank (white) region.}
  \label{fig:proof}
  \end{figure}

\begin{lem}\label{lem:tau12}
When $\frac{\varepsilon}{\sigma_m h}> \max ( \lambda_{1} (\mu),  \lambda_{2} (\mu))$, both
$ \widehat{\tau}^\omega_{\eps,h,1}  (\mu)$ and $ \widehat{\tau}^\omega_{\eps,h,2}  (\mu)$ are finite, and they satisfy 
\begin{equation}
 \widehat{\tau}^\omega_{\eps,h,1}  (\mu)=\tau_{\varepsilon, h, 1}(\mu)\leq  \widehat{\tau}^\omega_{\eps,h,2}  (\mu)=\tau_{\varepsilon, h, 2}(\mu)\Longleftrightarrow \mu\leq \mu_S(\frac{\varepsilon}{\sigma_m h}) \Longleftrightarrow 
\lambda_S(\mu)\leq \frac{\varepsilon}{\sigma_m h}.
\end{equation}
Moreover, $\tau_{\varepsilon,h,1}(\mu_S(\frac{\varepsilon}{\sigma_m h}) )= \tau_{\varepsilon,h,2}(\mu_S(\frac{\varepsilon}{\sigma_m h}) ).$
\end{lem}


\begin{lem}\label{lem:monotone}
When restricted to $\{\mu: \frac{\varepsilon}{\sigma_m h}> \lambda_{2} (\mu)\}$,  $\tau_{\varepsilon,h,2}(\mu)$ is positive and  monotonically decreasing. 
When restricted to $\{\mu\in (\frac{1}{2\omega}, \min(\mu_S(\frac{\varepsilon}{\sigma_m h}), 1)]: \frac{\varepsilon}{\sigma_m h}> \lambda_{1} (\mu) \}$,  $\tau_{\varepsilon,h,1}(\mu)$ is positive and monotonically increasing. 
\end{lem}
\begin{proof}
The definitions of $\lambda_{j}(\mu)$ ensures  $\tau_{\varepsilon,h,j}(\mu)$  is positive with $j=1, 2$ for the considered $\mu$. The monotonicity of $\tau_{\varepsilon,h,2}(\mu)$ directly comes from its being linear, and what remained will be devoted to showing the monotonicity of  $\tau_{\varepsilon,h,1}(\mu)$. 

Based on the definition of $\tau_{\varepsilon,h,1}(\mu)$ in \eqref{eqn:tau1},  we know that when $\frac{\varepsilon}{\sigma_m h}> \lambda_{1} (\mu)$, we have $2\varepsilon^2 \widehat{C}_{inv}||v^2||_\infty-(1-\mu)(\mu-\frac{1}{2\omega} )\sigma_m^2h^2>0$, and 
\begin{align*}
\tau'_{\varepsilon,h,1}(\mu)=\frac{2\varepsilon^2h^2\sigma_m\left(2\varepsilon^2 \widehat{C}_{inv}|| v^2||_\infty-(\mu-\frac{1}{2\omega})^2\sigma_m^2h^2\right)
}{(2\varepsilon^2 \widehat{C}_{inv}||v^2||_\infty-(1-\mu)(\mu-\frac{1}{2\omega} )\sigma_m^2h^2)^2}.
\end{align*}
As a result, the sign of $\tau'_{\varepsilon,h,1}(\mu)$,  same as that of
$
q(\mu):=
2\varepsilon^2 \widehat{C}_{inv}|| v^2||_\infty-(\mu-\frac{1}{2\omega})^2\sigma_m^2h^2,$
will inform  about  the monotonicity of $\tau_{\varepsilon,h,1}(\mu)$.

Consider the two roots of $q(\mu)$, which are  $\tilde{\mu}_{1,2}=\tilde{\mu}_{1,2}(\frac{\varepsilon}{\sigma_m h})=\frac{1}{2\omega}\mp \dfrac{\varepsilon}{\sigma_m h}\sqrt{2\widehat{C}_{inv}||v^2||_\infty }.$ And $q(\mu)>0$ when  $\mu \in(\tilde{\mu}_1,\tilde{\mu}_2)$.  Note that $\tilde{\mu}_1<\frac{1}{2\omega}$. One can further show that $\tilde{\mu}_2(\lambda)> \mu_S(\lambda), \; \forall \lambda>0$ as below.
\begin{align*}
\mu_S(\lambda)<\tilde{\mu}_2(\lambda)&\Longleftrightarrow\frac{1}{2\omega}+\frac{1}{2}\lambda\frac{\widehat{C}_{inv}||v^2||_\infty }{C_{inv}||v||_\infty}< \frac{1}{2\omega}+\lambda\sqrt{2\widehat{C}_{inv}||v^2||_\infty}
\notag\\
&\Longleftrightarrow \frac{\widehat{C}_{inv}||v^2||_\infty }{2C_{inv}||v||_\infty}<\sqrt{2\widehat{C}_{inv}||v^2||_\infty }\Longleftrightarrow {\mathcal K}> 1.
\end{align*}
Hence $\left(\frac{1}{2\omega}, \min(\mu_S(\frac{\varepsilon}{\sigma_m h}), 1)\right]\subset (\tilde{\mu}_1,\tilde{\mu}_2)$. And the monotonicity of $\tau_{\varepsilon,h,1}(\mu)$ will follow.





\end{proof}


 
\begin{lem}\label{lem:lambda_star}
Assume $\lambda>0$. 
\begin{itemize}
\item[(i)] $\lambda>\lambda_\star \Longleftrightarrow \lambda > \lambda_{2} (\mu_S(\lambda))$
\item[(ii)] When $\lambda\leq \widehat{\lambda}_\star$,  then $\lambda>\lambda_\star \Longleftrightarrow  \lambda> \lambda_{1} (\mu_S(\lambda))$.
\item[(iii)] When $\lambda_\star<\frac{\varepsilon}{\sigma_m h}\leq  \widehat{\lambda}_\star$,   we have $\frac{\varepsilon}{\sigma_m h}> \max ( \lambda_{1} (\mu),  \lambda_{2} (\mu))|_{\mu=\mu_S(\frac{\varepsilon}{\sigma_m h})}$.
\end{itemize}
\end{lem}
\begin{proof}
To prove (i), we proceed from the definitions of  $ \lambda_{2} (\mu)$ and $\mu_S(\lambda)$, and get
\begin{align}
\lambda > \lambda_{2} (\mu_S(\lambda)) \Longleftrightarrow & \lambda >\frac{1-\frac{1}{2\omega}-\frac{1}{2}\lambda \frac{\widehat{C}_{inv}||v^2||_\infty}{C_{inv}||v||_\infty} }{4C_{inv}||v||_\infty} \\
\Longleftrightarrow &
\left(1+ \frac{\widehat{C}_{inv}||v^2||_\infty}{8(C_{inv}||v||_\infty)^2} \right)\lambda > \frac{ 1-\frac{1}{2\omega} }{4C_{inv}||v||_\infty}
\Longleftrightarrow\lambda>  \lambda_\star.
\end{align}
To prove (ii),  we first notice  $\mu_S(\lambda)>\frac{1}{2\omega}$ holds  when $\lambda>0$. With $\lambda\leq \widehat{\lambda}_\star$, equivalently $\mu_S(\lambda)\leq 1$,  we then have
\begin{align}
\lambda> \lambda_{1} (\mu_S(\lambda)) \Longleftrightarrow & \lambda > \sqrt{\frac{(1-\frac{1}{2\omega}-\frac{1}{2}\lambda \frac{\widehat{C}_{inv}||v^2||_\infty}{C_{inv}||v||_\infty})\frac{1}{2}\lambda \frac{\widehat{C}_{inv}||v^2||_\infty}{C_{inv}||v||_\infty} }{2\widehat{C}_{inv}||v^2||_\infty} } \notag\\
\Longleftrightarrow&
\lambda > \left(1-\frac{1}{2\omega} - \frac{1}{2} \lambda \frac{\widehat{C}_{inv}||v^2||_\infty}{C_{inv}||v||_\infty} \right) \frac{1}{4C_{inv}||v||_\infty}
\Longleftrightarrow \lambda> \lambda_\star.
\end{align}
(iii) is a direct result of (i) and (ii) of this lemma.
\end{proof}
 

  

\subsection{Proof of Theorem \ref{thm:stab}: unconditionally stable region, $k\geq 1$}\label{sec:un_k1}

Based on Theorem \ref{thm:mu} and the definition of (unconditional) stability, the IMEX1-LDG method is unconditionally stable if and only if $ \Dt_{\textrm{stab}} (\varepsilon, h)=\infty$, which is equivalent to 
\begin{align}
\label{eq:uncon-pk}
\frac{\varepsilon}{\sigma_m h}\leq  \max_{\mu\in(\frac{1}{2\omega},1]}\big(\min\big( \lambda_{1} (\mu), \lambda_{2} (\mu) \big)\big).
\end{align}
Using Lemma \ref{lem:muS:mono}-(i) and  Lemma \ref{lem:mu_star}-(i), one has
\begin{equation}
 \min\big( \lambda_{1} (\mu), \lambda_{2} (\mu)\big)
=\left\{
\begin{array}{ll}
 \lambda_{1} (\mu), & \textrm{if}\; \mu\leq \mu_\star, \\
 \lambda_{2} (\mu), &\textrm{if}\; \mu\geq \mu_\star
\end{array}
\right.
\end{equation}
where $\mu_\star\in (\frac{1}{2\omega},  1)$, and  the inequality \eqref{eq:uncon-pk} will be simplified as
  \begin{align}
\frac{\varepsilon}{\sigma_m h}\leq 
\max\left(
 \max_{\mu\in(\frac{1}{2\omega},\mu_\star]}  \lambda_{1} (\mu), 
\max_{\mu\in[\mu_\star,1]} \lambda_{2} (\mu)
\right)
=\max\left( \lambda_{1} (\mu_\star),  \lambda_{2} (\mu_\star)\right)=\lambda_\star.
\label{eq:uncon-pk:2}
\end{align}
This gives the result in Theorem \ref{thm:stab} regarding the unconditional stability  when  $k\geq 1$.





\subsection{Proof of Theorem \ref{thm:stab}:  conditionally stable region, $1\leq k\leq 9$, $ \frac{\varepsilon}{\sigma_m h}>\lambda_\star$}



In this subsection, we focus on $\varepsilon$ and $h$ that satisfy $ \frac{\varepsilon}{\sigma_m h}>\lambda_\star$. For such $\varepsilon, h$, we have $\Delta t_{\textrm{stab}}(\varepsilon, h)<\infty$, and the IMEX1-LDG method is conditionally stable. Based on the $\mu$-stability result in Theorem \ref{thm:mu}, we want to optimize the time step condition by properly choosing $\mu$ from the admissible set, hence to get $ \Dt_{\textrm{stab}} (\varepsilon, h)$ and establish the remaining result in Theorem \ref{thm:stab}. 

\subsubsection{When $ \frac{\varepsilon}{\sigma_m h}>\widehat{\lambda}_\star$}
We start with the simplest case, that is when  $ \frac{\varepsilon}{\sigma_m h}>\widehat{\lambda}_\star$. According to Lemma \ref{lem:mu_star}-(iii),  for such $\varepsilon, h$, one has $\frac{\varepsilon}{\sigma_m h}>\max( \lambda_{1} (\mu), \lambda_{2} (\mu)), \forall \mu\in(\frac{1}{2\omega}, 1]$, hence $\tau_{\varepsilon, h, j}(\mu)<\infty$, $j=1, 2$, and $$\Delta t_{\textrm{stab}}(\varepsilon, h)=\max_{ \mu\in(\frac{1}{2\omega}, 1]}\min (\tau_{\varepsilon, h, 1}(\mu), \tau_{\varepsilon, h, 2}(\mu)).$$ 

Using the property of $\mu_S(\lambda)$ in Lemma \ref{lem:muS:mono},  we get
\begin{equation}
\frac{\varepsilon}{\sigma_m h}>\widehat{\lambda}_\star\Leftrightarrow \mu_S(\frac{\varepsilon}{\sigma_m h})>\mu_S(\widehat{\lambda}_\star)=1\Rightarrow \mu_S(\frac{\varepsilon}{\sigma_m h}) \geq \mu, \quad \forall \mu\in (\frac{1}{2\omega}, 1].
\label{eq:Dt0}
\end{equation}
Now following the comparison property in Lemma \ref{lem:tau12} and the monotonicity of $\tau_{\varepsilon, h, 1}(\mu)$ in Lemma \ref{lem:monotone}, we have, when $\frac{\varepsilon}{\sigma_m h}>\widehat{\lambda}_\star$,
$$\Delta t_{\textrm{stab}}(\varepsilon, h)=\max_{ \mu\in(\frac{1}{2\omega}, 1]\cap(\frac{1}{2\omega},\mu_S(\frac{\varepsilon}{\sigma_m h})] }\tau_{\varepsilon, h, 1}(\mu)=
\tau_{\varepsilon, h, 1}\left( \min(\mu_S(\frac{\varepsilon}{\sigma_m h}),   1) \right).$$


 
 
 \subsubsection{When $ {\lambda}_\star<\frac{\varepsilon}{\sigma_m h}\leq \widehat{\lambda}_\star$} %

From here on, we assume $ \frac{\varepsilon}{\sigma_m h}\in (\lambda_\star, \widehat{\lambda}_\star]$. 
The relation in \eqref{eq:Dt0} implies  
\begin{equation}
\mu_S(\frac{\varepsilon}{\sigma_m h})\leq 1.
\label{eq:Dt4}
\end{equation}  We decompose  $(\frac{1}{2\omega}, 1]$ into three disjoint sets $S_j(\varepsilon, h)$, $j=1, 2, 3$, defined as
\begin{align*}
S_1(\varepsilon,h) &= \left\{\mu\in  (\frac{1}{2\omega}, 1]: \frac{\varepsilon}{\sigma_m h}>\max( \lambda_{1} (\mu), \lambda_{2} (\mu))\right\},\\
S_2(\varepsilon,h) &= \left\{\mu\in  (\frac{1}{2\omega}, 1]:  \lambda_{1} (\mu)<\frac{\varepsilon}{\sigma_m h}\leq  \lambda_{2} (\mu)\right\},\\
S_3(\varepsilon,h) &= \left\{\mu\in  (\frac{1}{2\omega}, 1]:  \lambda_{2} (\mu)<\frac{\varepsilon}{\sigma_m h}\leq  \lambda_{1} (\mu))\right\}.
\end{align*}
One can refer to Figure \ref{fig:proof} to visualize the decomposition for a constant weight function $\omega\equiv 1$.
And correspondingly, 
$$\Delta t_{\textrm{stab}}(\varepsilon, h)=\max_{\mu\in (\frac{1}{2\omega}, 1]}\min( \widehat{\tau}^\omega_{\eps,h,1}  (\mu),  \widehat{\tau}^\omega_{\eps,h,2}  (\mu))=\max_{j=1, 2, 3} \Delta t_{\textrm{stab}}^{(j)}(\varepsilon, h),$$
where  
$\Delta t_{\textrm{stab}}^{(j)}(\varepsilon, h):=\max_{\mu\in S_j(\varepsilon, h)}\min( \widehat{\tau}^\omega_{\eps,h,1}  (\mu),  \widehat{\tau}^\omega_{\eps,h,2}  (\mu)).$ 
Next we will calculate  $\Delta t_{\textrm{stab}}^{(1)}(\varepsilon, h)$, and then show
$\Delta t_{\textrm{stab}}^{(1)}(\varepsilon, h)\geq \Delta t_{\textrm{stab}}^{(j)}(\varepsilon, h), j=2, 3$, therefore
\begin{equation} \Delta t_{\textrm{stab}}(\varepsilon, h)= \Delta t_{\textrm{stab}}^{(1)}(\varepsilon, h).\end{equation}



\bigskip
 {\bf Step 1:  To compute $\Delta t_{\textrm{stab}}^{(1)}(\varepsilon, h)$.}  When $\mu\in S_1(\varepsilon, h)$, we have $ \widehat{\tau}^\omega_{\eps,h,1}  (\mu)=\tau_{\varepsilon,h,1}(\mu)<\infty$, $ \widehat{\tau}^\omega_{\eps,h,2}  (\mu)=\tau_{\varepsilon,h,2}(\mu)<\infty$.
 Based on the comparison result in  Lemma \ref{lem:tau12},  and the property of $\mu_S(\lambda)$ in Lemma \ref{lem:muS:mono}, there holds
 
\begin{align}
\label{eq: mtau_12}
\min\left( \widehat{\tau}^\omega_{\eps,h,1}  (\mu), \widehat{\tau}^\omega_{\eps,h,2}  (\mu) \right)
&=\begin{cases}
\tau_{\varepsilon,h,1}(\mu), \quad & \mu \in (\frac{1}{2\omega},\mu_S(\frac{\varepsilon}{\sigma_m h}) ],\\
\tau_{\varepsilon,h,2}(\mu), \quad & \mu\in  (\mu_S(\frac{\varepsilon}{\sigma_m h}),1].
\end{cases}
\end{align}
 
With  $ \lambda_\star<\frac{\varepsilon}{\sigma_m h}\leq \widehat{\lambda}_\star$, based on 
Lemma \ref{lem:lambda_star}-(iii), we will get $\mu_S(\frac{\varepsilon}{\sigma_m h})\in S_1(\varepsilon, h)$.    
By further using the monotonicity of $\tau_{\varepsilon, h, j}(\mu), j=1, 2$ in Lemma \ref{lem:monotone}, and the fact $\tau_{\varepsilon,h,1}(\mu_S(\frac{\varepsilon}{\sigma_m h}) ) = \tau_{\varepsilon,h,2}(\mu_S(\frac{\varepsilon}{\sigma_m h}) )$ in Lemma \ref{lem:tau12},  when  $ \frac{\varepsilon}{\sigma_m h}\in (\lambda_\star, \widehat{\lambda}_\star]$,
\begin{align}
\label{eq:Dt1}
\Delta t_{\textrm{stab}}^{(1)}(\varepsilon, h) &= \max_{\mu \in S_1(\varepsilon,h)}\Big( \min\left( \widehat{\tau}^\omega_{\eps,h,1}  (\mu), \widehat{\tau}^\omega_{\eps,h,2}  (\mu)\right)\Big)\notag\\
& =
\tau_{\varepsilon,h,1}( \mu_S(\frac{\varepsilon}{\sigma_m h}) ) 
= \tau_{\varepsilon,h,1} \left( \min(\mu_S(\frac{\varepsilon}{\sigma_m h}), 1) \right).
\end{align}






\bigskip
 {\bf Step 2:  To show  $\Delta t_{\textrm{stab}}^{(2)}(\varepsilon, h)\leq \Delta t_{\textrm{stab}}^{(1)}(\varepsilon, h)$.}  When $\mu\in S_2(\varepsilon, h)$, we have $ \widehat{\tau}^\omega_{\eps,h,1}  (\mu)=\tau_{\varepsilon,h,1}(\mu)<\infty$, $ \widehat{\tau}^\omega_{\eps,h,2}  (\mu)=\infty$, hence $\min( \widehat{\tau}^\omega_{\eps,h,1}  (\mu), \widehat{\tau}^\omega_{\eps,h,2}  (\mu) ) = \tau_{\varepsilon,h,1}(\mu)$.

For any $\mu\in S_2(\varepsilon, h)$, based on Lemma \ref{lem:mu_star}, we have $\mu \leq \mu_\star$. Moreover, using the fact of 
$\mu_S(\lambda_\star) = \mu_\star$ and the monotonicity of $\mu_S(\lambda)$ in Lemma \ref{lem:muS:mono}, as well as the assumption $\frac{\varepsilon}{\sigma_m h}>\lambda_\star$, we have for $\mu\in S_2(\varepsilon, h)$,
$$
\mu \leq \mu_\star=\mu_S(\lambda_\star)< \mu_S(\frac{\varepsilon}{\sigma_m h}).
$$
Finally, we can once again use  the monotonicity of  $\tau_{\varepsilon,h,1}(\mu)$ in Lemma \ref{lem:monotone}, and conclude 
\begin{align}
\label{eq:Dt2}
\Delta t_{\textrm{stab}}^{(2)}(\varepsilon, h) &= \max_{\mu \in S_2(\varepsilon,h)}\Big( \min\left( \widehat{\tau}^\omega_{\eps,h,1}  (\mu), \widehat{\tau}^\omega_{\eps,h,2}  (\mu)\right)\Big)= \max_{\mu \in S_2(\varepsilon,h)}\tau_{\varepsilon, h, 1}(\mu)\notag\\
&\leq\tau_{\varepsilon,h,1}\left(\mu_S(\frac{\varepsilon}{\sigma_m h})\right) = \Delta t_{\textrm{stab} }^{(1)}(\varepsilon, h).
\end{align}

\bigskip
 {\bf Step 3:  To show  $\Delta t_{\textrm{stab}}^{(3)}(\varepsilon, h)\leq \Delta t_{\textrm{stab}}^{(1)}(\varepsilon, h)$.}   When $\mu\in S_3(\varepsilon, h)$, we have $ \widehat{\tau}^\omega_{\eps,h,1}  (\mu)=\infty$, $ \widehat{\tau}^\omega_{\eps,h,2}  (\mu)=\tau_{\varepsilon,h,2}(\mu)<\infty$, hence $\min( \widehat{\tau}^\omega_{\eps,h,1}  (\mu), \widehat{\tau}^\omega_{\eps,h,2}  (\mu) ) = \tau_{\varepsilon,h,2}(\mu)$.
 


Given any  $\mu \in S_3(\varepsilon, h)$, we know $ \lambda_{2} (\mu)<\frac{\varepsilon}{\sigma_m h}\leq  \lambda_{1} (\mu)$. This, combined with Lemma \ref{lem:mu_star}, implies $\mu>\mu_\star$, and additionally
\begin{align}
\frac{\varepsilon}{\sigma_m h} \leq  \lambda_{1} (\mu)< \lambda_S(\mu) \Leftrightarrow \mu> \mu_S\left(\frac{\varepsilon}{\sigma_m h}\right). 
\end{align}
The equivalency is based on the monotonicity of $\mu_S(\lambda)$ in Lemma \ref{lem:muS:mono}. Finally one can use the 
 monotonicity of $\tau_{\varepsilon,h,2}(\mu)$ in Lemma \ref{lem:monotone}, and conclude
\begin{align}
\Delta t_{\textrm{stab}}^{(3)}(\varepsilon, h) &=
\max_{\mu \in S_3(\varepsilon,h)}\Big( \min\left( \widehat{\tau}^\omega_{\eps,h,1}  (\mu), \widehat{\tau}^\omega_{\eps,h,2}  (\mu)\right)\Big)= \max_{\mu \in S_3(\varepsilon,h)}\tau_{\varepsilon, h, 2}(\mu)\notag\\
& \leq \tau_{\varepsilon,h,2}(\mu_S(\frac{\varepsilon}{\sigma_m h}) )=\tau_{\varepsilon,h,1}\left(\mu_S(\frac{\varepsilon}{\sigma_m h}) \right)
=\Delta t_\textrm{stab}^{(1)}(\varepsilon, h).\notag
\end{align}




  












\subsection{Norms, inverse inequalities, and more notation}
\label{sec:prelim}

We introduce some standard norms $||\phi||=||\phi||_{L^2({\Omega_x})}$, $|||\phi|||=(\langle ||\phi||^2\rangle)^{1/2}$,  and weighted norms
$||\phi||_s= ||\sqrt{\sigma_s}\phi||$, $|||\phi|||_s =|||\sqrt{\sigma_s}\phi|||.$
For a bounded function $\psi(v)$ of $v$,  without confusion we will write $||\psi||_\infty=||\psi||_{L^\infty({\Omega_v})}$.  Even though for our specific examples with ${\Omega_v}=[-1, 1]$ or $\{-1, 1\}$, we have $||v||_\infty=||v^2||_\infty=1$, we still keep $||v||_\infty$ and $||v^2||_\infty$ in most results, to possibly inform about the case with a more general bounded velocity space ${\Omega_v}$. 



In our analysis, the following inverse inequalities will be frequently used, and they are fairly standard in finite element analysis: there exist constants $C_{inv}=C_{inv}(k)$ and $\widehat{C}_{inv}=\widehat{C}_{inv}(k)$, such that for any $\phi \in P^k([a,b])$,
\begin{subequations}
\label{eq:inv}
\begin{align}
&|\phi(y)|^2 (b-a) \leq C_{inv}\int_a^b |\phi(x)|^2 dx, \qquad\mbox{with}\; y=a \;\mbox{or}\; b \label{eq:inv:1},\\
&(b-a)^2 \int_a^b |\phi'(x)|^2 dx \leq \widehat{C}_{inv}\int_a^b |\phi (x)|^2 dx.
\label{eq:inv:2}
\end{align}
\end{subequations}
Particularly, $C_{inv}(k)|_{k=0}=1$. Next lemma states a property of the inverse constants $\widehat{C}_{inv},  C_{inv}$.

\begin{lem} With $\Omega_v=[-1,1]$ or $\Omega_v=\{-1,1\}$, and with  $\widehat{C}_{inv}$, $C_{inv}$ from \eqref{eq:inv},  we define 
\begin{equation}
\label{def:K}
{\mathcal K}={\mathcal K}(k)=\frac{8(C_{inv}||v||_\infty)^2}{\widehat{C}_{inv}||v^2||_\infty}=\frac{8(C_{inv})^2}{\widehat{C}_{inv}}.
\end{equation} 
 Then at least for $k=1, 2, \cdots, 9$, we have ${\mathcal K}>1.$
 \label{lem:mK}
\end{lem}
\begin{proof}
Based on Lemmas 1-2 in \cite{Matt_Li} and a linear scaling, one can take $C_{inv}=(k+1)^2$ and $ \widehat{C}_{inv}=12k^4$, which can be used to verify ${\mathcal K}>1$ directly for $k=1, 2, \cdots 9$.
\end{proof}
 Sharper values of $C_{inv}(k)$ and $\widehat{C}_{inv}(k)$ can be numerically obtained for {\em each} $k$ by solving an eigenvalue problem (see  Section 4.1 in \cite{Matt_Li}), hence  one can check numerically whether ${\mathcal K}>1$ holds or not for larger $k$. Given the temporal accuracy of the IMEX1-LDG method is first order, it is more than enough for us to consider $k\leq 9$ in our analysis.

  For convenient  reference, we summarize  in Table \ref{tab:notation} the definitions of some notation arising from analysis, including $\lambda_\star, \widehat{\lambda}_\star$ and $\mu_\star$, which all depend on inverse constants 
 hence on $k$.  They also depend on the  weight function $\omega$ and the velocity space $\Omega_v$. 
 The same table also includes the definitions of 
 ${\mathcal K}$ in \eqref{def:K},  a function $\mu_S(\lambda)$ and its inverse $\lambda_S(\mu)$, as well as two more functions $\lambda_j(\mu), j=1, 2$. The place where each notation appears for the first time is also included.

  
\begin{table}[htb]
	\caption{Some notation  (with the possible $\omega$-dependence suppressed) and the place of the first appearance}
		\begin{center}
		\begin{tabular}{lc}
					\hline	
		notation & the first appearance \\
			\hline		
			${\mathcal K}=\frac{8(C_{inv}||v||_\infty)^2}{\widehat{C}_{inv}||v^2||_\infty} $  & \eqref{def:K} \\\hline							
$\lambda_\star=\frac{2(1-1/(2\omega) )C_{inv}||v||_\infty}{\widehat{C}_{inv}||v^2||_\infty+8(C_{inv}||v||_\infty)^2}$ & \eqref{eq:thm:stab3:o}\\\hline
			$\mu_\star=\frac{1+\frac{1}{2\omega}{\mathcal K} }{1+{\mathcal K}}=\frac{\widehat{C}_{inv}||v^2||_\infty+4(C_{inv}||v||_\infty)^2/\omega}{\widehat{C}_{inv}||v^2||_\infty+8(C_{inv}||v||_\infty)^2}$ & \eqref{eq:thm:stab5:o:1}\\\hline	
			$\mu_S(\lambda)=\frac{1}{2\omega}+\frac{1}{2}\lambda\frac{\widehat{C}_{inv}||v^2||_\infty}{C_{inv}||v||_\infty}$
			&\eqref{eq:thm:stab5:o:3} \\ \hline
$\lambda_S(\mu)=\mu_S^{-1}(\mu)=2(\mu-\frac{1}{2\omega})\frac{C_{inv}||v||_\infty}{\widehat{C}_{inv}||v^2||_\infty}$ & Lemma \ref{lem:muS:mono}\\\hline	
$\widehat{\lambda}_\star=\lambda_S(1)=2(1-\frac{1}{2\omega})\frac{C_{inv}||v||_\infty}{\widehat{C}_{inv}||v^2||_\infty}$ & \eqref{eq:thm:stab5:o:3}\\\hline	
$ \lambda_{1} (\mu)=\sqrt{\frac{(1-\mu)(\mu-\frac{1}{2\omega}) }{2\widehat{C}_{inv} ||v^2||_\infty} }, \quad  \lambda_{2} (\mu)=\frac{1-\mu}{4C_{inv}||v||_\infty} $ & \eqref{eqn:lam12}\\\hline
		\end{tabular}
		\end{center}
		\label{tab:notation}
\end{table}


\section{Proof for AP property: Theorem \ref{thm:AP}}
\label{sec:ap:proofs}

We will first build some preparatory results in Lemma \ref{lem:AP}, before proving the main result on the AP property in Theorem \ref{thm:AP}. The three assumptions in Section \ref{sec:ap} still hold.  Let $\{\Psi_j\}_{j=1}^{N_k}$ be an orthonormal basis of $U_h^k$ with respect to the standard $L^2$ inner product of $L^2(\Omega_x)$. Recall the initialization is via the $L^2$ projection onto $U_h^k$, namely, $\rho^0_ {\eps,\Dt,h} =\pi_h \rho_\varepsilon^0$, $g^0_ {\eps,\Dt,h} =\pi_h g_\varepsilon^0$,  $u_ {\eps,\Dt,h} ^0=\pi_h(\sigma_s^{-1}q^0_\varepsilon)$. We also define $W_{ {\Dt,h} }^0=\pi_h W_0$ for $W=\rho, g$, and $u_ {\Dt,h} ^0=\pi_h(\sigma_s^{-1}q_0)$.

\begin{lem}\label{lem:AP} The following results hold. \\
(i) $q_\varepsilon^0\rightharpoonup q_0$  in $L^2(\Omega_x)$ as $\varepsilon\rightarrow 0$. \\
(ii) $\lim_{\varepsilon\rightarrow 0}\rho_ {\eps,\Dt,h} ^0=\rho^0_{ {\Dt,h} }$, 
$\lim_{\varepsilon\rightarrow 0}u_ {\eps,\Dt,h} ^0=u_{ {\Dt,h} }^0$ and
\begin{align}
&\lim_{\varepsilon\rightarrow 0}\langle \zeta, g_ {\eps,\Dt,h} ^0(x,\cdot)\rangle =\langle  \zeta, g_ {\Dt,h} ^0(x, \cdot)\rangle,\quad \forall  \zeta\in L^2(\Omega_v),\;\;\forall x\in \Omega_x,\\
\label{eq:g_initial_limit2}
&\lim_{\varepsilon\rightarrow 0}\langle  \zeta, (g_ {\eps,\Dt,h} ^0, \psi) \rangle =\langle  \zeta, (g_ {\Dt,h} ^0, \psi)\rangle,\quad \forall  \zeta\in L^2(\Omega_v), \;\;\forall \psi\in L^2({\Omega_x}).
\end{align}
(iii) $\sup_{\varepsilon}||W_ {\eps,\Dt,h} ^0||<\infty$, where $W=\rho, g, u$.\\
(iv) $\sup_{\{0<\varepsilon<\varepsilon_0(h)\}}||W_ {\eps,\Dt,h} ^1||=C_W(k,\Delta t, h,{\Omega_v})<\infty$, where $W=\rho,  u$.
\end{lem}
\begin{proof}[\bf{Proof}]
\par

 (i) Start with any $\phi\in C^\infty_0(\Omega_x)$, then
 \begin{equation}
 (q_0, \phi)= - (\rho_0, \phi_x)=-\lim_{\varepsilon\rightarrow0}(\rho_\varepsilon^0, \phi_x)=\lim_{\varepsilon\rightarrow0} (q_\varepsilon^0, \phi).
\end{equation}
 This result can be extended to any $\phi\in L^2(\Omega_x)$, hence $q_\varepsilon^0\rightharpoonup q_0$  in $L^2(\Omega_x)$ as $\varepsilon\rightarrow 0$, due to the uniform boundedness of $||q_\varepsilon^0||$ in $\varepsilon$ in Assumption 2 and $C^\infty_0(\Omega_x)$ being dense in  $L^2(\Omega_x)$.

\smallskip
(ii) With $W_\varepsilon^0$ weakly convergent to $W_0$ in $L^2(\Omega_x)$, for $W=\rho, q$, we have
$$\lim_{\varepsilon\rightarrow0}\rho_ {\eps,\Dt,h} ^0=\lim_{\varepsilon\rightarrow0}\pi_h\rho_\varepsilon^0=\lim_{\varepsilon\rightarrow0}\sum_{j=1}^{N_k}(\rho_\varepsilon^0,\Psi_j)\Psi_j=\sum_{j=1}^{N_k}\lim_{\varepsilon\rightarrow0}(\rho_\varepsilon^0,\Psi_j)\Psi_j=\sum_{j=1}^{N_k}(\rho_0,\Psi_j)\Psi_j=\pi_h\rho_0=\rho^0_ {\Dt,h} ,$$
$$\lim_{\varepsilon\rightarrow0}u_ {\eps,\Dt,h} ^0=\lim_{\varepsilon\rightarrow0}\pi_h(\sigma_s^{-1}q_\varepsilon^0)=\lim_{\varepsilon\rightarrow0}\sum_{j=1}^{N_k}(\sigma_s^{-1}q_\varepsilon^0,\Psi_j)\Psi_j=\sum_{j=1}^{N_k}(\sigma_s^{-1}q_0,\Psi_j)\Psi_j=\pi_h(\sigma_s^{-1}q_0)=u^0_ {\Dt,h} .$$

Now we consider any $ \zeta\in L^2(\Omega_v)$. With $\langle  \zeta g_\varepsilon^0\rangle$ weakly convergent to $\langle  \zeta g_0\rangle$ in $L^2(\Omega_x)$, we have for any $x\in\Omega_x$,
\begin{align}
\lim_{\varepsilon\rightarrow0}\langle  \zeta, g_ {\eps,\Dt,h} ^0(x,\cdot)\rangle
&=\lim_{\varepsilon\rightarrow0}\langle  \zeta, \sum_{j=1}^{N_k} (g_\varepsilon^0,\Psi_j)\Psi_j(x)\rangle=\sum_{j=1}^{N_k}\lim_{\varepsilon\rightarrow0}(\langle  \zeta g_\varepsilon^0\rangle,\Psi_j)\Psi_j(x) \notag\\
&=\sum_{j=1}^{N_k}(\langle  \zeta g_0\rangle,\Psi_j)\Psi_j(x)
=\langle  \zeta, g^0_ {\Dt,h} (x,\cdot)\rangle.
\end{align}
And \eqref{eq:g_initial_limit2} can be proved similarly.


\smallskip
(iii) Note that 
$$ |||g_ {\eps,\Dt,h} ^0|||^2=\langle ||g_ {\eps,\Dt,h} ^0||^2\rangle=\langle\sum_{j=1}^{N_k}(g_\varepsilon^0,\Psi_j)^2\rangle\leq |||g_\varepsilon^0|||^2\sum_{j=1}^{N_k}||\Psi_j||^2=N_k|||g_\varepsilon^0|||^2,$$
$$ ||u_ {\eps,\Dt,h} ^0||=||\pi_h(\sigma_s^{-1} q_\varepsilon^0)||\leq
||\sigma_s^{-1} q_\varepsilon^0||\leq 
\sigma_m^{-1}|| q_\varepsilon^0||.$$
With Assumption 2, we have $\sup_{\varepsilon}|||W_ {\eps,\Dt,h} ^0|||<\infty, W=g, u $. Similar proof goes to $\rho$.




\smallskip
(iv) 
Based on \eqref{eq:FDG:1T}, one has
\begin{align}
(\rho_ {\eps,\Dt,h} ^1,\testR)&=\Delta t \omega \langle v^2 \rangle l_h(u^1_ {\eps,\Dt,h} ,\testR)+(\rho^0_ {\eps,\Dt,h} ,\testR)
\notag\\
&-\Delta t l_h(\langle v (g^0_ {\eps,\Dt,h} +\omega v u_ {\eps,\Dt,h} ^0)\rangle,\testR)-(\sigma_a \rho_ {\eps,\Dt,h} ^1,\testR)
,\; \forall \testR\in U_h^k.
\label{eq:b:1}
\end{align}
Take $\testR=\rho_ {\eps,\Dt,h} ^1$, use $l_h(u^1_ {\eps,\Dt,h} ,\rho_ {\eps,\Dt,h} ^1)=-(\sigma_s u^1_ {\eps,\Dt,h} ,u^1_ {\eps,\Dt,h} )$ based on \eqref{eq:prop:1} and Assumption 3 for $\omega$, we get when $\varepsilon<\varepsilon_0(h)$,
 \begin{align}
&|| \rho_ {\eps,\Dt,h} ^1||^2+(\sigma_a\rho_ {\eps,\Dt,h} ^1,\rho_ {\eps,\Dt,h} ^1)+\frac{2\sigma_m\Delta t }{3}\langle v^2\rangle ||u_ {\eps,\Dt,h} ^1||^2\notag\\
\leq&
(\rho^0_ {\eps,\Dt,h} ,\rho^1_ {\eps,\Dt,h} )-\Delta t l_h(\langle v (g^0_ {\eps,\Dt,h} +\omega vu_ {\eps,\Dt,h} ^0)\rangle,\rho^1_ {\eps,\Dt,h} )
\label{eq:b:2}
\end{align}
Following some standard steps to  apply Cauchy-Schwarz inequality, Young inequality, inverse inequality (see, e.g. Lemma 3.9 in \cite{jang2014analysis}), based on Assumption 3, we can find a constant $C(k, \Delta t, h, \Omega_v)$ such that
\begin{align}
&|(\rho^0_ {\eps,\Dt,h} ,\rho^1_ {\eps,\Dt,h} )-\Delta t l_h(\langle v (g^0_ {\eps,\Dt,h} +\omega vu_ {\eps,\Dt,h} ^0,\rho^1_ {\eps,\Dt,h} )\rangle|\notag\\
\leq &C(k, \Delta t, h, \Omega_v)\left(||\rho_ {\eps,\Dt,h} ^0||+|||g_ {\eps,\Dt,h} ^0|||+||u_ {\eps,\Dt,h} ^0 ||\right)||\rho_ {\eps,\Dt,h} ^1||.
\label{eq:b:3}
\end{align}
Combining \eqref{eq:b:2}-\eqref{eq:b:3}, with $\sigma_a(x)\geq 0$, we obtain
\begin{align*}
&\sup_{0<\varepsilon<\varepsilon_0(h)}  || \rho_ {\eps,\Dt,h} ^1||\leq C(k, \Delta t, h, \Omega_v)\sup_\varepsilon(||\rho_ {\eps,\Dt,h} ^0||+|||g_ {\eps,\Dt,h} ^0|||+||u_ {\eps,\Dt,h} ^0 ||) <\infty,\\
&\sup_{0<\varepsilon<\varepsilon_0(h)} || u_ {\eps,\Dt,h} ^1||\leq \sqrt{\frac{3}{2\sigma_m\Delta t\langle v^2\rangle}}C(k, \Delta t, h, \Omega_v)\sup_\varepsilon(||\rho_ {\eps,\Dt,h} ^0||+|||g_ {\eps,\Dt,h} ^0|||+||u_ {\eps,\Dt,h} ^0 ||) <\infty.
\end{align*}
\end{proof}

We are ready to prove  Theorem \ref{thm:AP} on the AP property of the IMEX1-LDG method. 
\begin{proof}[\bf{Proof of Theorem \ref{thm:AP}.}] Let the mesh size $h$ be fixed.



{\sf Step 1:}   we first show that 
$\sup_{0<\varepsilon<\varepsilon_0(h)}||U_ {\eps,\Dt,h} ^n||<\infty$
for any $\Delta t$, $n\geq 1$, where  $W=\rho, g, q, u$. 
First note that when $\varepsilon<\varepsilon_0(h)$, from Assumption 3, we have $2>\omega>\frac{2}{3}$ and $\mu=\frac{3}{4}\in (\frac{1}{2\omega}, 1]$. Based on the $\mu$-stability result in Theorem \ref{thm:mu}, we have
\begin{align}
&||\rho_ {\eps,\Dt,h} ^{n+1}||^2+\varepsilon^2|||g_ {\eps,\Dt,h} ^{n}|||^2+\Delta t\sigma_m\left( \frac{1}{4}|||g_ {\eps,\Dt,h} ^{n}|||^2+\frac{2}{3 }\langle v^2 \rangle ||u_ {\eps,\Dt,h} ^{n+1}||^2\right)
\notag\\
\leq& E_{h,\mu=\frac{3}{4}}^{n+1}\leq  E_{h,\mu=\frac{3}{4}}^{n}\leq\cdots
\leq E_{h,\mu=\frac{3}{4}}^1
\notag\\
\leq& ||\rho_ {\eps,\Dt,h} ^1||^2+\varepsilon^2|||g_ {\eps,\Dt,h} ^0|||^2+\Delta t\sigma_M \left( \frac{1}{4}|||g_ {\eps,\Dt,h} ^{0}|||^2+2\langle v^2 \rangle ||u_ {\eps,\Dt,h} ^1||^2\right).
\label{eq:prove-AP} \end{align}
Moreover from \eqref{eq:FDG:1T:c}, we have $||q_ {\eps,\Dt,h} ^n||^2=(\sigma_s u_ {\eps,\Dt,h} ^n,q_ {\eps,\Dt,h} ^n)$, hence  $||q_ {\eps,\Dt,h} ^n||\leq \sigma_M ||u_ {\eps,\Dt,h} ^n||$. In combination of Lemma \ref{lem:AP},  the finiteness of $\sup_{0<\varepsilon<\varepsilon_0(h)}||W_ {\eps,\Dt,h} ^n||$, $\forall n\geq 1$ follows for  $W=\rho, g, q, u$.

\smallskip
{\sf Step 2:}
With Lemma \ref{lem:AP}, we only need to establish \eqref{eq:f:lim}  for any  $n\geq 1$. This is equivalent to show that for any given sequence $\{\varepsilon_m\}_{m=1}^{\infty}$, satisfying  $\lim_{m\rightarrow \infty} \varepsilon_m=0$ (we no longer emphasize that $\varepsilon$ considered here is bounded above by $\varepsilon_0(h)$),  we have
\begin{subequations}
\label{eq:aa:1}
\begin{align}
\label{eq:aa:1-1}
& \lim_{m\rightarrow \infty}W_{\varepsilon_m, \Delta t, h}^n  =W_{\Delta t, h}^n, \quad W=\rho, q, u,\\
& \lim_{m\rightarrow \infty}\langle \zeta, g_{\varepsilon_m, \Delta t, h}^n(x,\cdot)\rangle =\langle \zeta, g_ {\Dt,h} ^n(x, \cdot)\rangle,\quad \forall \zeta\in L^2(\Omega_v),\;\;\forall x\in \Omega_x,
\label{eq:aa:1-2}\\
& \lim_{m\rightarrow \infty}\langle \zeta, (g_{\varepsilon_m, \Delta t, h}^n, \psi) \rangle =\langle \zeta, (g_ {\Dt,h} ^n, \psi)\rangle,\quad \forall \zeta\in L^2(\Omega_v), \;\;\forall \psi\in L^2({\Omega_x}),
\label{eq:aa:1-3}
\end{align}
\end{subequations}
 for some $W_ {\Dt,h} ^n\in U_h^k$, with $W=\rho, q, u$,  and $g_ {\Dt,h} ^n\in G_h^k$, $\forall n\geq 1$.
   Let $W$ be any of $\rho, q, u$. Given that $U_h^k$ is finite dimensional,  the finiteness of $\sup_m ||W_{\varepsilon_m, \Delta t, h}^n||$ from  {\sf Step 1} implies that there is a subsequence $\{W_{\varepsilon_{m_r}, \Delta t, h}^n\}_{r=1}^{\infty}$ converging in $U_h^k$ under {\em any} norm as $r\rightarrow \infty$.  Let the limit be 
   \begin{equation}
   \label{eq:aa:2}
    W_{\Delta t, h}^n=\lim_{r\rightarrow \infty}W_{\varepsilon_{m_r}, \Delta t, h}^n, \quad W=\rho, q, u.
    \end{equation}

We now turn to $\{g_{\varepsilon_m, \Delta t, h}^n\}_{m=1}^\infty$. 
Note that each $g^n_{\varepsilon_m, \Delta t, h}$ can be written as $g^n_{\varepsilon_m, \Delta t, h}(x,v)=\sum_{j=1}^{N_k}\alpha_{\varepsilon_m}^{(j)}(v)\Psi_j(x)$, with
$||| g^n_{\varepsilon_m}|||=\left(\sum_{j=1}^{N_k} ||\alpha_{\varepsilon_m}^{(j)}||^2_{L^2(\Omega_v)}\right)^{1/2}$. This, in addition to the finiteness of
$\sup_m |||g_{\varepsilon_m, \Delta t, h}^n|||$ in {\sf Step 1},  indicates that $\sup_r ||\alpha_{\varepsilon_{m_r}}^{(j)}||^2_{L^2(\Omega_v)}$  is bounded for any $j=1,\cdots, N_k$. As a Hilbert space, $L^2(\Omega_v)$ is weakly sequentially compact, that is, $\{\alpha_{\varepsilon_{m_r}}^{(j)}\}_{r=1}^\infty$ has a subsequence which is weakly convergent in $L^2(\Omega_v)$. Without loss of generality, this subsequence is still denoted as $\{\alpha_{\varepsilon_{m_r}}^{(j)}\}_{r=1}^\infty$, and the weak limit when $r\rightarrow \infty$ is denoted as $\alpha_0^{(j)}\in L^2(\Omega_v)$, $\forall j$. We now define
$g_{\Delta t, h}^n(x,v)=\sum_{j=1}^{N_k}\alpha_0^{(j)}(v)\Psi_j(x)$.
It is clear that
{
$g_{\Delta t, h}^n\in G_h^k$.
}
For any $\zeta\in L^2(\Omega_v)$, and any $x\in\Omega_x$,
\begin{equation}
\label{eq:aa:3}
\lim_{r\rightarrow \infty}\langle \zeta, g_{\varepsilon_{m_r}, \Delta t, h}^n(x,\cdot)\rangle=\sum_{j=1}^{N_k} \left(\lim_{r\rightarrow\infty}\langle \zeta, \alpha_{\varepsilon_{m_r}}^{(j)}\rangle \right)\Psi_j(x)=\sum_{j=1}^{N_k} \langle \zeta, \alpha_0^{(j)}\rangle \Psi_j(x)=\langle \zeta, g_{\Delta t, h}^n(x,
\cdot)\rangle.
\end{equation}
Furthermore, we have $\forall \zeta\in L^2(\Omega_v),\;\forall \psi\in L^2({\Omega_x}),$
\begin{equation}
\label{eq:aa:4}
\lim_{r\rightarrow \infty}\langle \zeta, (g_{\varepsilon_{m_r}, \Delta t, h}^n, \psi)\rangle
=\sum_{j=1}^{N_k} \left(\lim_{r\rightarrow\infty}\langle \zeta, \alpha_{\varepsilon_{m_r}}^{(j)}\rangle \right)(\Psi_j, \psi)
=\langle \zeta, (g_{\Delta t, h}^n, \psi)\rangle=(\langle \zeta g_{\Delta t, h}^n\rangle, \psi).
\end{equation}


Use \eqref{eq:aa:2}-\eqref{eq:aa:4} for $n\geq 1$ as well as the similar result in Lemma \ref{lem:AP} for $n=0$, with   $\zeta$ taken when needed as
$v$,  $v {\bf 1}_{\{v>0\}}$,  $v{\bf 1}_{\{v<0\}}$, $v\zeta(v)$, $v\zeta(v){\bf 1}_{\{v>0\}}$,  $v\zeta(v){\bf 1}_{\{v<0\}}$, also use the property \eqref{w:prop:0} for $\omega$, 
 we have for any $n\geq 0$,
\begin{subequations}
\label{eq:aa:5}
\begin{align}
&
 \lim_{r\rightarrow \infty} l_h(\langle v (g_{\varepsilon_{m_r}, \Delta t, h}^n+ \omega|_{\varepsilon=\varepsilon_{m_r}}v u_{\varepsilon_{m_r},\Delta t, h}^n\rangle ,\testR)=
 l_h(\langle v (g_ {\Dt,h} ^n+v u_{ {\Dt,h} }^n)\rangle ,\testR) ,\quad \forall \testR\in U_h^k, \\
&
\lim_{r\rightarrow \infty} \langle \zeta, b_{h,v}(g_{\varepsilon_{m_r}, \Delta t, h}^n,\testG)\rangle=\langle \zeta, b_{h,v}(g_ {\Dt,h} ^n,\testG)\rangle,\quad \forall \zeta\in L^2(\Omega_v), \; \forall \testG\in U_h^k.
\end{align}
\end{subequations}

Now with \eqref{eq:aa:2}-\eqref{eq:aa:5} and  Lemma \ref{lem:AP} for the initial data, the numerical scheme  \eqref{eq:FDG:1T} as  $r\rightarrow \infty$  becomes,
$\forall \testQ,\;\testU,\;\testR\,\;\testG \in U_h^k$
\begin{subequations}
\label{eq:aa:6}
\begin{align}
&(q_ {\Dt,h} ^{n+1},\testQ)+d_h(\rho_ {\Dt,h} ^{n+1},\testQ)=0,
\label{eq:aa:6.1}\\
&(\sigma_s u_ {\Dt,h} ^{n+1},\testU)=(q_ {\Dt,h} ^{n+1},\testU),
\label{eq:aa:6.4}\\
&(\frac{\rho_ {\Dt,h} ^{n+1}-\rho_ {\Dt,h} ^n}{\Delta t},\testR)+l_h(\langle v(g_ {\Dt,h} ^n+v u_ {\Dt,h} ^n)\rangle,\testR)=\langle v^2\rangle l_h(u_ {\Dt,h} ^{n+1},\testR)-(\sigma_a \rho_ {\Dt,h} ^{n+1},\testR),
\label{eq:aa:6.2}\\
&(\langle \zeta\sigma_s g_ {\Dt,h} ^{n+1}\rangle ,\testG)=\langle \zeta v\rangle d_h(\rho_ {\Dt,h} ^{n+1},\testG),\quad\forall \zeta\in L^2(\Omega_v),
\label{eq:aa:6.3}
\end{align}
\end{subequations}
for $n\geq 0$.
Furthermore, \eqref{eq:aa:6.1} and \eqref{eq:aa:6.3} lead to
\begin{equation} 
\langle  (\pi_h(\sigma_s g_ {\Dt,h} ^n)+vq_ {\Dt,h} ^n, \zeta\psi) \rangle=0
\quad \forall \zeta\in L^2(\Omega_v),
 \;\psi \in U_h^k, \;\; n\geq 1.
\label{eq:aa:7}
\end{equation}
With $g_ {\Dt,h} ^{n}\in G_h^k$ hence $\pi_h(\sigma_s g_ {\Dt,h} ^n)+vq_ {\Dt,h} ^n\in L^2({\Omega_v})\times U_h^k$, \eqref{eq:aa:7} equivalently becomes
\begin{equation}
\pi_h(\sigma_s g_ {\Dt,h} ^{n})=- v q_ {\Dt,h} ^{n},\quad n\geq 1.
\label{eq:aa:8}
\end{equation}
Moreover, from \eqref{eq:aa:6.4} and \eqref{eq:aa:8}, one can get $g_ {\Dt,h} ^n+vu^n_ {\Dt,h} =0, n\geq 1$, as shown  below.
\begin{align*}
0\leq \sigma_m ||| g_ {\Dt,h} ^n+vu^n_ {\Dt,h}  |||^2
&\leq \left\langle \left(\sigma_s (g_ {\Dt,h} ^n+vu^n_ {\Dt,h} ), g_ {\Dt,h} ^n+vu^n_ {\Dt,h} \right)\right\rangle\\
&=\left\langle\left( -vq_ {\Dt,h} ^n+vq_ {\Dt,h} ^n,g_ {\Dt,h} ^n+vu^n_ {\Dt,h} \right)\right\rangle=0.
\end{align*}

Compare \eqref{eq:aa:6} and \eqref{eq:aa:8} with what we want in \eqref{eq:lim:s}, one also needs to have 
$g_ {\Dt,h} ^0+vu^0_ {\Dt,h} =0$. This can be argued based on the initial data being  well-prepared in Assumption 1. To see this, $\forall \zeta\in L^2({\Omega_v}), \forall \psi\in U_h^k$, we proceed as follows,
\begin{align}
0&=\lim_{\varepsilon\rightarrow 0}\Big(\langle \zeta ( g_\varepsilon^{0}+v \sigma_s^{-1}q_\varepsilon^{0})\rangle,  \psi\Big)=\lim_{\varepsilon\rightarrow 0}\Big((\langle \zeta g_\varepsilon^{0}\rangle, \psi)+\langle v\zeta\rangle (q_\varepsilon^{0}, \sigma_s^{-1}\psi)\Big)\notag \\
&=(\langle \zeta g_0\rangle, \psi)+\langle v\zeta\rangle (q_0, \sigma_s^{-1}\psi)
=(\langle  \zeta  g_ {\Dt,h} ^{0}\rangle, \psi)+ \langle \zeta v\rangle ( u_ {\Dt,h} ^{0}, \psi),
\label{eq:aa:7.a}
\end{align}
and this gives $\langle \zeta (g_ {\Dt,h} ^{0}+vu_ {\Dt,h} ^{0}, \psi)\rangle=0$. Note that $g_ {\Dt,h} ^{0}+vu_ {\Dt,h} ^{0}\in L^2({\Omega_v})\times U_h^k$, therefore \eqref{eq:aa:7.a} is indeed $g_ {\Dt,h} ^{0}+vu_ {\Dt,h} ^{0}=0$, and we can conclude the limiting scheme in   \eqref{eq:lim:s}.

It is easy to see the limiting scheme \eqref{eq:lim:s} is a consistent discretization for \eqref{eq:AP1:3.1.1}. Its stability can be obtained similarly as Lemma \ref{lem:sol}, with 
\begin{align}
&||\rho^{n+1}_ {\Dt,h} ||^2+\Delta t\langle v^2\rangle ||u^{n+1}_ {\Dt,h} ||_s^2+(\sigma_a \rho^{n+1}_ {\Dt,h} , \rho^{n+1}_ {\Dt,h} )=(\rho^{n}_ {\Dt,h} , \rho^{n+1}_ {\Dt,h} )\notag\\
\Rightarrow &\frac{1}{2}||\rho^{n+1}_ {\Dt,h} ||^2+\Delta t\langle v^2\rangle \sigma_m ||u^{n+1}_ {\Dt,h} ||^2\leq \frac{1}{2}||\rho^{n}_ {\Dt,h} ||^2\leq \cdots\leq \frac{1}{2}||\rho^{0}_ {\Dt,h} ||^2\leq \frac{1}{2}||\rho_0||^2. 
\end{align}


Finally, with a standard contradiction argument and the uniqueness of the solution to the system
\eqref{eq:lim:s} (see Lemma \ref{lem:sol}), we conclude the limiting functions $\rho^n_ {\Dt,h} ,  q^n_ {\Dt,h} , g^n_ {\Dt,h} , u_ {\Dt,h} ^n$ are unique, and  \eqref{eq:aa:1} holds for the entire sequence. In the case that the velocity space $\Omega_v$ is discrete, the analysis related to the convergence of $g_ {\eps,\Dt,h} ^n(\cdot, v)$ for each $v$ is just as simple as that for $\rho_ {\eps,\Dt,h} ^n$ and $q_ {\eps,\Dt,h} ^n$, and the convergence is in a strong sense as in \eqref{eq:aa:-4}.

\end{proof}


\section{Introduction}
\label{sec0:intro}

In this paper, we continue our efforts in devising and advancing mathematical understanding of asymptotic preserving (AP) methods to solve time-dependent multi-scale kinetic transport equations within  the discontinuous Galerkin (DG) framework  \cite{AP1, jang2014analysis, peng2018stability}. Particularly,  we focus on establishing energy-type numerical stability and the  AP property for some methods  proposed  in  \cite{peng2018stability} for the following model equation,
 \begin{equation}
	\label{eq:mmd:LJ0}
	{\mathcal P}^\varepsilon: \qquad 
	\varepsilon f_t + v \partial_x f = \frac{\sigma_s}{\varepsilon}\left(\langle f\rangle - f\right)-\varepsilon \sigma_a f,
\end{equation}
with periodic boundary conditions.  The function $f=f(x,v,t)$ is  the probability distribution  
function of the particles, with the space variable $x\in \Omega_x\subset {\mathbb R}$, velocity variable $v\in \Omega_v\subset {\mathbb R}$, and time  $t\geq 0$. 
 $\sigma_s(x)>0 $ and $\sigma_a(x)\geq 0$ are the scattering and absorption coefficients, respectively.
${\mathcal L}(f)=\langle f\rangle-f$ defines a  scattering operator, where  $\langle f\rangle: = \int_{\Omega_v} f d \nu$  and  $\nu$ is a measure of the velocity space satisfying $\int_{\Omega_v} 1 d\nu=1$.   The parameter $\varepsilon>0$ is the dimensionless Knudsen number, defined as the ratio of the mean free path of the particles over the characteristic length of the system.  The model \eqref {eq:mmd:LJ0} is in a diffusive scaling, and as $\varepsilon\rightarrow 0$, it approaches its  diffusive limit
\begin{equation}
	\label{eq;mmd:LJ0_lim}
	{\mathcal P}^0: \qquad 
\partial_t\rho=\langle v^2\rangle \partial_x (\partial_x \rho/\sigma_s)-\sigma_a \rho.
\end{equation}
Here $\rho= \langle f \rangle$ is the  macroscopic density. Though seemingly simple, the equation in \eqref{eq:mmd:LJ0} provides
a prototype model to study many realistic problems such as in neutron transport or radiative transfer theory both numerically and mathematically. 




To  simulate multi-scale models like that in   \eqref{eq:mmd:LJ0} effectively and reliably for a broad range of value for $\varepsilon$, AP methods are widely recognized by the scientific community (see e.g. review papers  \cite{jin2010asymptotic, degond2011asymptotic}).  These methods  are designed for the  governing model with $\varepsilon>0$. 
Additionally when $\varepsilon\rightarrow0$, the methods become consistent  and stable discretizations for the limiting model as in \eqref{eq;mmd:LJ0_lim} even on under-resolved meshes with $\Delta x, \Delta t\gg \varepsilon$.  Hence, AP methods provide  a natural transition of different regimes in multi-scale simulations.  AP methods usually involve some level of implicit treatment to deal with the stiffness of the model when $\varepsilon\ll 1$.
It is known that stability alone does not guarantee the scheme to capture the correct asymptotic limit \cite{caflisch1997uniformly, naldi1998numerical}. 


 In our recent work \cite{peng2018stability}, a family of high order  AP methods,  termed as IMEX-LDG methods, are designed for \eqref{eq:mmd:LJ0}.   The methods are based on the reformulation of the equation, and involve local DG (LDG) discretization in space \cite{cockburn1998local}, globally stiffly accurate implicit-explicit (IMEX) Runge-Kutta (RK) methods in time \cite{boscarino2013imex}, and a judicially  chosen IMEX strategy. The reformulation has  two steps: micro-macro decomposition \cite{liu2004boltzmann,Lemou_No_BC}, and addition/subtraction of a $\omega$-weighted diffusive term \cite{boscarino2013imex}. The latter is introduced to obtain fully implicit limiting schemes as $\varepsilon\rightarrow 0$, to achieve  unconditional stability of the methods in the diffusive regime with $\varepsilon\ll 1$,  hence to  circumvent the otherwise stringent parabolic type time step condition in this regime, namely, $\Delta t=O(\Delta x^2)$, of  many AP schemes whose limiting schemes are explicit
 \cite{jin1998diffusive, klar1998asymptotic, Lemou_No_BC,  AP1}. Using globally stiffly accurate IMEX RK methods in time, 
 and LDG methods in space with suitable numerical fluxes,
  the IMEX-LDG methods project the numerical solutions to the local equilibrium at both inner stages and full RK steps in the limit of  $\varepsilon\rightarrow 0$, and this is important for the AP property and seemingly also for  accuracy (see appendix of \cite{peng2018stability}).  
 In \cite{peng2018stability},  unconditional stability in the diffusive regime is observed numerically, and is confirmed by a Fourier-type stability analysis applied to the  two-velocity telegraph equation with ${\Omega_v}=\{-1, +1\}$, and constant material properties $\sigma_s=1$, $\sigma_a=0$. 
 
In this work we restrict our attention to the IMEX1-LDG method, the sub-family of the methods in  \cite{peng2018stability} that  is  first order accurate  in time and arbitrary order in space, and examine it systematically for the model with the general 
 material properties, namely with the spatially varying scattering  and absorption coefficients  $\sigma_s(x)$ and $\sigma_a(x)$. Our main objectives are two-fold.  
The first is to establish unconditional stability in the diffusive regime with $\varepsilon\ll1$ {\em as well as}  uniform stability with respect to $\varepsilon$.  By following an energy approach as in \cite{liu2010analysis, jang2014analysis}, one can get uniform stability yet fails to capture the unconditional stability for $\varepsilon\ll 1$. Note the methods examined in \cite{liu2010analysis, jang2014analysis} in the limit of $\varepsilon\rightarrow 0$ are explicit. We instead propose and work with a new  notion of $\mu$-stability, and get the stability we want by better exploring the contribution of the scattering operator. 
The stability results up to this point depend on  a parameter $\mu$. An intricate algebraic-based optimization  with respect to the admissible $\mu$ is  subsequently followed, to further maximize the unconditional stability region, while  also maximizing the allowable time step size in the regime when the method is  conditionally stable.  
 As our second objective, a rigorous asymptotic analysis is proved to show the AP property based on the uniform stability. 
 To our best knowledge, our analysis is the first to capture unconditional stability when $\varepsilon\ll 1$ along with uniform stability property for the model \eqref{eq:mmd:LJ0} with general material properties. A general form of  the weight function $\omega$ is also for the first time considered in such stability analysis. 
In this work, we keep the velocity variable continuous, and our analysis can be easily adapted  when the velocity variable is further discretized such as by discrete ordinates or $P_N$ methods \cite{pomraning1973equations}.  Our analysis can also be extended to AP methods with the same IMEX strategy yet with other spatial discretizations, as long as they satisfy some key properties, such as the adjoint property in \eqref{eq:lem1:1} (also see Lemma 3.5 
in \cite{peng2018stability})
and the stabilization as in \eqref{eq:upwind} due to the upwind treatment.   Though not presented here,   {\em a priori} error estimates can follow similarly as in \cite{jang2014analysis}, and they are uniform in $\varepsilon$ for  smooth enough solutions with uniform bounds in $\varepsilon$ under the relevant Sobolev norms. What seems to be more challenging and left to our  future endeavor
is to obtain  the stability analysis for IMEX-LDG methods with higher order temporal accuracy.



Finally we want to briefly review some related  literature especially in establishing numerical stability 
of AP methods for kinetic transport models in a diffusive scaling.  One commonly used approach is Fourier type analysis. 
For the telegraph equation with ${\Omega_v}=\{-1, +1\}$,  an analytical time step condition is given in \cite{Lemou_No_BC} via Fourier analysis to ensure uniform  $L^2$-stability of a first order finite difference AP method, while in \cite{peng2018stability}, necessary conditions on $\varepsilon, \Delta x, \Delta t$  are obtained numerically for the $p$-th order IMEX-LDG AP scheme ($p=1, 2, 3$) to ensure an $L^2$ energy non-increasing in time. The results seem to be uniform in $\varepsilon$, with  unconditional stability captured for $\varepsilon\ll 1$.    
Klar and Unterreiter in \cite{klar2002uniform} considered a formally first-order in time and second-order in space  AP scheme for the one-group transport equation with ${\Omega_v}=[-1, 1]$ and established uniform stability by first establishing the result in Fourier space and then transforming  it back to the physical space.  Their analysis assumes the $H^1$ smoothness of the initial data.  It is known that Fourier-type  analysis requires uniform meshes and the models being linear and constant-coefficient. Energy-based stability analysis on the other hand does not pose these restrictions, yet they are not always easy to get. In \cite{liu2010analysis}, Liu and Mieussens revisited the first order AP method in \cite{Lemou_No_BC} for a more general kinetic transport model and proved uniform stability following an energy approach. A similar analysis
is carried out in \cite{jang2014analysis} for the first order in time DG-IMEX1 method in \cite{AP1}. Based on the uniform stability analysis, error estimates and rigorous asymptotic analysis are also established in \cite{jang2014analysis}.  
In both \cite{peng2018stability} and here in this work, we want to capture the unconditional stability in the diffusive regime in addition to the uniform stability. Few other  theoretical works, among many, for AP methods include uniform consistency \cite{caflisch1997uniformly, klar1998asymptotic}, uniform convergence \cite{golse1999convergence,filbet2013analysis} based on the commuting diagram of AP schemes (see Fig 1.1 in \cite{golse1999convergence}), and a recent work on uniform accuracy with   IMEX multi-step methods \cite{hu2019uniform}.



The remaining of the paper is organized as follows.  In Section \ref{sec1:scheme}, we review and extend  the IMEX1-LDG method in \cite{peng2018stability}  to our model \eqref{eq:mmd:LJ0} with general material properties.  
Section \ref{sec3:stab} presents  main results on numerical stability.   Here several theorems, including Theorem \ref{thm:stab_nmu} and Theorem \ref{thm:mu}, are stated to obtain uniform stability, while capturing the unconditional stability in the diffusive regime.  
An optimization step is carried out in Theorem \ref{thm:stab} to find the best value of the parameter $\mu$ in the notion of $\mu$-stability in order to optimize  the stability results. Once uniform stability is available, the AP property of the method is stated in Theorem \ref{thm:AP} in Section \ref{sec:ap}. The proofs of all major theorems are presented in Sections \ref{sec:stab:proofs}-\ref{sec:ap:proofs} for better readability.


\section{Proof for stability: Theorem \ref{thm:stab_nmu} and Theorem \ref{thm:mu} }
\label{sec:stab:proofs}

In this section, we will  present the proof for Theorem \ref{thm:mu} first and then Theorem \ref{thm:stab_nmu}.

\begin{proof}[\bf{Proof of Theorem  \ref{thm:mu}.}] 
 Let $n\geq 1$. Take $\testR=\rho_h^{n+1}$ in \eqref{eq:FDG:1T:b} and use
Lemma \ref{lem:1} and Proposition \ref{prop:sol}, 
we get
\begin{align}
& 
\big(\frac{\rho_h^{n+1}-\rho_h^n}{\Delta t}, \rho_h^{n+1}\big)+l_h(\langle v g_h^n\rangle, \rho_h^{n+1})-\omega \langle v^2\rangle l_h(u_h^{n+1}-u_h^n,\rho_h^{n+1})\notag\\
=& 
\big(\frac{\rho_h^{n+1}-\rho_h^n}{\Delta t}, \rho_h^{n+1}\big)
+\langle  v d_h(  \rho_h^{n+1},  g_h^n)\rangle+\omega \langle v^2\rangle (\sigma_s (u_h^{n+1}-u_h^n),u_h^{n+1}) \notag\\
=& 
\frac{1}{2\Delta t}\left(||\rho_h^{n+1}||^2-||\rho_h^n||^2+||\rho_h^{n+1}-\rho_h^n||^2\right)
+\langle  v d_h(  \rho_h^{n+1},  g_h^n)\rangle\notag\\
&\;\;\;\;\;+
\frac{\omega \langle v^2\rangle}{2}(
||u_h^{n+1}||_s^2 -||u_h^{n}||_s^2+||u_h^{n+1}-u_h^n||_s^2 )=-(\sigma_a \rho_h^{n+1}, \rho_h^{n+1}).\label{eq:heheda:1}
\end{align}
Take $\testG=\varepsilon^2 g_h^{n+1}$ in \eqref{eq:FDG:1T:d}, integrate over $\Omega_v$ in $v$,  and shift index $n$ to $n-1$, we get
\begin{align}
&\varepsilon^2\langle\big(\frac{g_h^{n}-g_h^{n-1} }{\Delta t},g_h^n\big)\rangle+\varepsilon \langle b_{h,v}(g_h^{n-1},g_h^n)\rangle-\langle v d_h(\rho_h^{n},g_h^n)\rangle\notag\\
=&\frac{\varepsilon^2}{2\Delta t}\left(|||g_h^n|||^2-|||g_h^{n-1}|||^2+|||g_h^n-g_h^{n-1}|||^2\right)+\varepsilon \langle b_{h,v}(g_h^{n-1},g_h^n)\rangle-\langle v d_h(\rho_h^{n},g_h^n)\rangle \notag\\
=&-|||g_h^n|||_s^2 -\varepsilon^2\langle (\sigma_a g_h^n, g_h^n)\rangle.\label{eq:heheda:2}
\end{align}
Now we sum up \eqref{eq:heheda:1} and \eqref{eq:heheda:2},  with $E_h^n$ defined in \eqref{eq:energy_nmu}, and have
 \begin{align}
&\frac{1}{2\Delta t}(E_h^{n+1}-E_h^{n} ) +\frac{1}{2\Delta t}(||\rho_h^{n+1}-\rho_h^n||^2+\varepsilon^2|||g_h^n-g_h^{n-1} |||^2)
+\frac{\omega  \langle v^2\rangle}{2}||u_h^{n+1} -u_h^n||_s^2\notag\\
&+|||g_h^n|||_s^2+ \langle vd_h(\rho_h^{n+1}-\rho_h^n,g_h^n)\rangle
-\varepsilon\langle{b_{h,v}}(g_h^n-g_h^{n-1},g_h^n)\rangle+\varepsilon\langle {b_{h,v}}(g_h^n,g_h^n)\rangle \leq 0. 
\label{eq:process2}
\end{align}


\smallskip
To estimate $\langle vd_h(\rho_h^{n+1}-\rho_h^n,g_h^n)\rangle$ in \eqref{eq:process2},  
based on the scheme \eqref{eq:FDG:1T:a}-\eqref{eq:FDG:1T:c} and apply the Cauchy-Schwartz inequality, 
we get 
\begin{align}
|\langle vd_h(\rho_h^{n+1}-\rho_h^n,g_h^n)\rangle|&=| d_h(\rho_h^{n+1}-\rho_h^n,  \langle vg_h^n \rangle) | =|(q_h^{n+1}-q_h^n,\langle vg_h^n\rangle) |
\notag\\
&=|(\sigma_s (u_h^{n+1}-u_h^n), \langle v g_h^n\rangle)| \leq
\sqrt{\langle v^2\rangle} |||g_h^n|||_s\; ||u_h^{n+1}-u_h^n||_s.
\label{eq:process:dh}
\end{align}
The two terms in \eqref{eq:process2} involving the bilinear form $b_{h,v}$ can be handled  similarly as in \cite{jang2014analysis} (see its Lemma 3.2, particularly equations (3.22)-(3.24)). 
More specifically,  with $\langle g_h^m\rangle=0$ in Proposition \ref{prop:sol}, utilizing the upwind treatment in the proposed scheme for $v\partial_x g$, in addition to a few applications of inverse inequalities \eqref{eq:inv} and Young's inequality, it can be shown that
\begin{equation}
\langle {b_{h,v}}(g_h^n,g_h^n)\rangle =\left\langle\sum_i\frac{|v|}{2}[g_h^n]_{i-\frac{1}{2}}^2 \right\rangle,\label{eq:upwind}
\end{equation}
\begin{align}
&| \langle {b_{h,v}}(g_h^{n}-g_h^{n-1} ,g_h^n) \rangle |\notag\\
&\leq (\frac{\theta}{\sigma_m}+\eta) |||g_h^n-g_h^{n-1}|||^2+\frac{\sigma_m}{4\theta}\left\langle \sum_i \int_{I_i} (v\partial_x g_h^n)^2 dx\right\rangle+\frac{C_{inv} }{4\eta h}\sum_i \left\langle (v[g_h^n]_{i-\frac{1}{2}})^2\right \rangle\notag\\
&\leq(\frac{\theta}{\sigma_m}+\eta) |||g_h^n-g_h^{n-1}|||^2 +\frac{\widehat{C }_{inv}  ||v^2||_\infty}{4\theta h^2} |||g_h^n|||_s^2
+\frac{C_{inv} ||v||_{\infty} }{2\eta h}\left\langle \frac{|v|}{2}\sum_i [g_h^n]_{i-\frac{1}{2}}^2  \right\rangle.
\label{eq:process:lambda}
\end {align}
Here $\theta$ and $\eta$ are two  positive constants, which will be specified later.

One important step in this proof is to split $|||g_h^n|||_s^2$
in  \eqref{eq:process2} into two terms, each playing different roles,  according to some parameter $\mu\in[0,1]$ (additional conditions required for  $\mu$ will soon become clear),  with one term further rewritten based on the parallelogram identity,
\begin{align}
\label{eq:process:change_g}
|||g_h^n|||_s^2
=\mu|||g_h^n|||_s^2+ (1-\mu)\Big(\frac{1}{2}|||g_h^n|||_s^2-\frac{1}{2}|||g_h^{n-1}|||_s^2
+
\frac{1}{4}|||g_h^n-g_h^{n-1} |||_s^2 + \frac{1}{4}|||g_h^n+g_h^{n-1} |||_s^2\Big).
 \end{align}
We now combine \eqref{eq:process2}-\eqref{eq:process:change_g},  with the discrete energy $E_{h,\mu}^{n}$ defined in \eqref{eqn:energy-mu},
 and reach
\begin{align}
 \label{eq:process:final}
& \frac{1}{2\Delta t} (E^{n+1}_{h,\mu}-E^n_{h,\mu})+\varepsilon\left(1-\frac{C_{inv}||v||_\infty}{2\eta h} \right) \left\langle \frac{|v|}{2}\sum_i[g_h^n]_{i-\frac{1}{2}}^2 \right\rangle \\
           &+\left(\frac{\varepsilon^2}{2\Delta t}+\frac{1-\mu}{4}\sigma_m-\varepsilon(\frac{\theta}{\sigma_m}+\eta) \right)|||g_h^n-g_h^{n-1} |||^2+(1-\mu)||| \frac{g_h^n+g_h^{n-1} }{2}|||_s^2+\frac{1}{2\Delta t}||\rho_h^{n+1}-\rho_h^n||^2 \notag \\
           &+\frac{\omega \langle v^2\rangle}{2} ||u_h^{n+1}-u_h^n||_s^2-\sqrt{\langle v^2\rangle} |||g_h^n|||_s\; ||u_h^{n+1}-u_h^n||_s 
          +\left(\mu-\varepsilon\frac{\widehat{C}_{inv}||v^2||_\infty  }
           {4\theta h^2}\right) |||g_h^n|||_s^2 \leq 0.\notag
\end{align}
In order for the discrete energy to be non-increasing, namely,  $ E_{h,\mu}^{n+1}\leq E_{h,\mu}^{n}$, we require the  quadratic form in the final row of \eqref{eq:process:final} to be non-negative, and this can be ensured by a non-negative discriminant, leading to 
\begin{equation}
\mu-\varepsilon\frac{\widehat{C}_{inv}||v^2||_\infty }{4\theta h^2}\geq \frac{1}{2\omega}.
\label{eq:ineq:1}
\end{equation}
Additionally, we also require
\begin{align}
1-\frac{C_{inv} ||v||_\infty}{2\eta h}&\geq 0
 \label{eq:ineq:2} ,\\
\frac{\varepsilon^2}{2\Delta t}+\frac{1-\mu}{4}\sigma_m-\varepsilon(\frac{\theta}{\sigma_m}+\eta)&\geq 0
\label{eq:ineq:3}.
\end{align}

The inequality \eqref{eq:ineq:1} implies that  $\mu$ needs to be restricted as   $\mu>\frac{1}{2\omega}$. 
We now choose
$$\frac{\theta}{\sigma_m}=\eta=\frac{1}{2}\left(\frac{\varepsilon}{2\Delta t}+\frac{1-\mu}{4\varepsilon}\sigma_m \right), $$
and with this,  \eqref{eq:ineq:3} is satisfied automatically, while \eqref{eq:ineq:2} becomes
\begin{equation}
\frac{\varepsilon^2}{\Delta t}\geq\frac{4 C_{inv}||v||_\infty\varepsilon-(1-\mu)\sigma_m h }{2 h} ,
\label{eq:result:1}
\end{equation}
and \eqref{eq:ineq:1} is now 
\begin{equation}
\frac{\varepsilon^2}{\Delta t}\geq
\frac{2\varepsilon^2\widehat{C}_{inv}||v^2||_\infty-(1-\mu)(\mu-\frac{1}{2\omega})\sigma_m^2h^2}{2 (\mu-\frac{1}{2\omega} )\sigma_m h^2}.
\label{eq:result:2}
\end{equation}
When  $\frac{\varepsilon}{\sigma_m h}\leq\frac{1-\mu}{4C_{inv}||v||_\infty }$, the right hand side of \eqref{eq:result:1} is non-positive, hence \eqref{eq:result:1} holds for any time step $\Delta t$. Otherwise, the  time step needs to satisfy $\Delta t\leq \tau_{\varepsilon,h,2}(\mu)$ with $\tau_{\varepsilon,h,2}(\mu)$ defined in \eqref{eqn:tau2}. 
Similarly,  when $\frac{\varepsilon}{\sigma_m h}\leq\sqrt{\frac{(1-\mu)(\mu-\frac{1}{2\omega} ) } {2\widehat{C}_{inv}||v^2||_\infty} } $, the right hand side of \eqref{eq:result:2} is non-positive, hence \eqref{eq:result:2} holds for any time step $\Delta t>0$. Otherwise,
the time step needs to satisfy $\Delta t\leq \tau_{\varepsilon,h,1}(\mu)$ with $\tau_{\varepsilon,h,1}(\mu)$ defined in \eqref{eqn:tau1}. 
The discussions so far can be summarized into  the claims in Theorem \ref{thm:mu} when $k\geq 1$. 

When $k=0$, we have  $\partial_x g_h^{n}=0$,  and the estimate in \eqref {eq:process:lambda} can be replaced by
\begin{equation}
| \langle {b_{h,v}}(g_h^{n}-g_h^{n-1} ,g_h^n) \rangle \leq\eta |||g_h^n-g_h^{n-1}|||^2+\frac{C_{inv} ||v||_{\infty} }{2\eta h}\left\langle \frac{|v|}{2}\sum_i [g_h^n]_{i-\frac{1}{2}}^2  \right\rangle,
\label{eq:process:k0}
\end{equation}
and all analysis up to \eqref{eq:ineq:3}  holds without the terms containing $\theta$. Specifically, \eqref{eq:ineq:1}-\eqref{eq:ineq:3} become
\begin{equation}
\mu \geq \frac{1}{2\omega}, \quad
1-\frac{C_{inv} ||v||_\infty}{2\eta h}\geq 0,\quad
\frac{\varepsilon^2}{2\Delta t}+\frac{1-\mu}{4}\sigma_m -\varepsilon\eta\geq 0
\label{eq:ineq:4-6}.
\end{equation}
Now take
$$
\eta= \frac{\varepsilon}{2\Delta t}+\frac{1-\mu}{4\varepsilon}\sigma_m
$$
in \eqref{eq:ineq:4-6}, and follow a similar analysis as  above, one reaches the results for $k=0$.
\end{proof}

\begin{proof}[\bf{Proof of Theorem  \ref{thm:stab_nmu}.}]   The proof can be established by starting with  the equation \eqref{eq:process2}, and then following almost the identical  analysis  in \cite{jang2014analysis} (particularly, see equations (3.22), (3.26)-(3.28), (3.36)-(3.41) in \cite{jang2014analysis}), together with $|||g_h^n|||_s^2\geq \sigma_m |||g_h^n|||^2$ to deal with the general scattering coefficient $\sigma_s(x)$. The details are omitted.
\end{proof}






\section{Numerical stability}
\label{sec3:stab}


In this section, we will establish numerical stability for the IMEX1-LDG method following {\em an energy approach}. At the continuous level, one can derive an energy relation
\begin{equation}
\label{eq:stab:c}
\frac{1}{2} \frac{d}{dt}\Big(||\rho||^2+\varepsilon^2|||g|||^2\Big)=-\int_{\Omega_v}\int_{\Omega_x} \Big(\sigma_s g^2+\sigma_a (\rho+\varepsilon g)^2\Big) dxdv
\end{equation}
for the model \eqref{eq:mmd:LJ0},
implying the energy $|||f|||^2(t)=||\rho||^2(t)+\varepsilon^2|||g|||^2(t)$ does not grow in time.  Our numerical stability is a discrete analogue. Particularly,  we want to confirm that the method is unconditionally stable in the diffusive regime when $\varepsilon\ll 1$ and it is uniformly stable in $\varepsilon$, 
with a general form of the weight function $\omega=\omega(\varepsilon/(\sigma_m h))$  taken into account. 
Without loss of generality, we assume the mesh is uniform with $h=h_i, \forall i$. Our results can be extended to general meshes when $\frac{\max_i h_i}{\min_i h_i}$ is bounded uniformly during mesh refinement.  For easy  readability, we will present and discuss the main results in this section, and defer the proofs   to Sections \ref{sec:stab:proofs}-\ref{sec:stab:proofs:opt}.

The natural first attempt is to follow a similar analysis as in \cite{jang2014analysis}, and this will lead to the stability result in next theorem.

\begin{thm}\label{thm:stab_nmu} 
The following stability result holds for the IMEX1-LDG method, defined  as \eqref{eq:FDG:1T} with \eqref{eq:bilinear:def}-\eqref{eq:flux}, 
\begin{align}
E_h^{n+1}\leq E_h^{n}, \;\;\forall n\geq 1, \;\; \textrm{with} \; \; E_h^n:=||\rho_h^n||^2+\varepsilon^2|||g_h^{n-1}|||^2+\Delta t\omega\langle v^2\rangle ||u_h^n||_s^2,
\label{eq:energy_nmu}
\end{align}
under the time step condition,
 \begin{equation}
 \label{eq:cfl}
 \Delta t\leq \Delta t_{stab}=\left\{
  \begin{array}{ll}
     \frac{2h}{\alpha_2 \alpha_3}(\sigma_m h+\alpha_3\varepsilon), &\mbox{for}\; k=0,\\
  \frac{h}{\alpha_1+\alpha_2\alpha_3} (\sigma_mh+\min(\varepsilon, \frac{\alpha_2h}{\alpha_1})\alpha_3), &\mbox{for}\; k\geq 1.
  \end{array}
 \right.
 \end{equation}
Here $\alpha_i, i=1, 2, 3$ are defined in terms of the inverse constants and the velocity space, namely,
\begin{equation}
\label{eq:alpha}
\alpha_1=(||v||_\infty^2+\langle v^2\rangle) \widehat{C}_{\textrm{inv}}\;,\quad\alpha_2= 2(||v||_\infty+\langle|v|\rangle)C_{\textrm{inv}}\;,\quad
\alpha_3=2 ||v||_\infty C_{\textrm{inv}}.
\end{equation}
\end{thm}
Note that the time step condition in \eqref{eq:cfl} is essentially the same as the one for the DG-IMEX1 method defined in \cite{jang2014analysis}.  This theorem, on one hand, gives uniform stability  with respect to $\varepsilon$,  which is important for the AP property of the method, see Section \ref{sec:ap} and Section \ref{sec:ap:proofs}, also \cite{jang2014analysis}. On the other hand, the theorem {\em  fails} to capture the unconditional stability property of the method in the diffusive regime when  $\varepsilon\ll 1$.  

The main reason that Theorem \ref{thm:stab_nmu} missed the unconditional stability we observed numerically and predicted by Fourier analysis in \cite{peng2018stability} is that the damping mechanism associated with the scattering operator (see the right hand side term in \eqref{eq:stab:c})  has not been fully utilized in the analysis. By better exploring the contribution of the scattering operator, new   stability results can be established and they will  capture the unconditional stability  property of the method.  This indeed is one main contribution of this work.  The new stability analysis will be based on a new discrete energy $E_{h, \mu}^n$.



\begin{defn}
\label{def:mu-stable}
For any given  constant $\mu\in [0, 1]$, we define a discrete energy
\begin{align}
E_{h,\mu}^n= ||\rho_h^n||^2+\varepsilon^2|||g_h^{n-1}|||^2+\omega\Delta t \langle v^2\rangle ||u_h^n||_s^2+\Delta t(1-\mu)||| g_h^{n-1}|||_s^2.
\label{eqn:energy-mu}
\end{align} 
The IMEX1-LDG method is said to be {\em $\mu$-stable} if it satisfies
 \begin{align}
 E_{h,\mu}^{n+1}\leq E_{h,\mu}^n,\qquad \forall n\geq 1.
 \label{eq:mu:stab}
  \end{align}
If the method is $\mu$-stable for some $\mu\in[0,1]$, then it is said to be {\em stable}. If the scheme being $\mu$-stable (resp. stable) is independent of the time step size $\Delta t$, the method is further said to be {\em unconditionally $\mu$-stable} (resp. {\em unconditionally stable}). Note that $E_{h, 1}^n=E_n^n$.
\end{defn}

With respect to the $\mu$-stability above, a new stability result will be stated in next theorem under the assumption  $\omega>1/2$. When the weight function is $\omega\equiv 1$, this assumption always holds. In general, with the property $\omega\rightarrow 1$ as $\varepsilon\rightarrow 0$ in \eqref{w:prop:0}, the stability result can at least capture the property of the method in the diffusive regime.  
\begin{thm}{\bf ($\mu$-stability: $\omega>\frac{1}{2}$)}
\label{thm:mu} When 
$\omega>\frac{1}{2}$,
 the following $\mu$-stability results hold for the IMEX1-LDG method, defined as  \eqref{eq:FDG:1T} with \eqref{eq:bilinear:def}-\eqref{eq:flux}. 
 %
 \begin{itemize}
 %
 \item [(i)]  When $k=0$ and with any fixed  $\mu\in[\frac{1}{2\omega},1]$,  if
\begin{equation}
\label{eq:thm:stab1}
\frac{\varepsilon}{\sigma_m h}\leq   \lambda_{0}(\mu):=\frac{1-\mu}{2C_{inv}||v||_\infty }=\frac{1-\mu}{2||v||_\infty },
\end{equation}
the IMEX1-LDG method is unconditionally $\mu$-stable.
  Otherwise,
 the method is conditionally $\mu$-stable when the time step satisfies
 \begin{equation}
\label{eq:thm:stab2}
\Delta t\leq\tau_{\varepsilon, h, 0}(\mu):=\frac{2\varepsilon^2 h}{2 C_{inv}||v||_\infty \varepsilon-(1-\mu)\sigma_m h}=\frac{2\varepsilon^2 h}{2 ||v||_\infty \varepsilon-(1-\mu)\sigma_m h}.
\end{equation}
Here we have used $C_{inv}(k)|_{k=0}=1$. The result can be expressed  more compactly as $\Delta t\leq  \widehat{\tau}_{\eps,h,0}  (\mu)$, by introducing an extended real-valued function
\begin{equation}
 \widehat{\tau}_{\eps,h,0}  (\mu)=\begin{cases}
                         \infty,    \quad  &\text{if}\;\frac{\varepsilon}{\sigma_m h} \leq \lambda_0(\mu),\\
		        \tau_{\varepsilon,h,0}(\mu)   \quad &\text{otherwise}.
		        \end{cases}
\label{eqn:Gtau0}
\end{equation}
And the scheme is unconditionally $\mu$-stable if and only if $ \widehat{\tau}_{\eps,h,0}  (\mu)=\infty$. 


\item [(ii)]  When $k\ge 1$ and with any fixed  $\mu\in(\frac{1}{2\omega},1]$, if
 \begin{equation}
\label{eq:thm:stab3}
\frac{\varepsilon}{\sigma_m h}\leq \min\left( \lambda_{1} (\mu),  \lambda_{2} (\mu)\right),
\end{equation}
the IMEX1-LDG method is unconditionally $\mu$-stable. 
Otherwise, the method is conditionally $\mu$-stable when the time step satisfies
\begin{equation}
\label{eq:thm:stab4}
\Delta t\leq\left\{
\begin{array}{ll}
\tau_{\varepsilon,h,1}(\mu),&\textrm{if}\;
 \lambda_{1} (\mu)<\frac{\varepsilon}{\sigma_m h}\leq \lambda_{2} (\mu),\\
\tau_{\varepsilon,h,2}(\mu),&\textrm{if}\;
 \lambda_{2} (\mu)<\frac{\varepsilon}{\sigma_m h}\leq \lambda_{1} (\mu),\\
\min(\tau_{\varepsilon,h,1}(\mu), \tau_{\varepsilon,h,2}(\mu)),&\textrm{if}\;\frac{\varepsilon}{\sigma_m h}\geq \max\left( \lambda_{1} (\mu),  \lambda_{2} (\mu)\right).
\end{array}
\right.
\end{equation}
Here
\begin{subequations}
\label{eq:tau12}
\begin{align}
 \lambda_{1} (\mu)&:=\sqrt{\frac{(1-\mu)(\mu-\frac{1}{2\omega} ) } {2\widehat{C}_{inv}||v^2||_\infty} }, \qquad  \lambda_{2} (\mu):=\frac{1-\mu}{4C_{inv}||v||_\infty }, \label{eqn:lam12}\\
\tau_{\varepsilon,h,1}(\mu)&:=\frac{ 2\varepsilon^2(\mu-\frac{1}{2\omega} )h^2 \sigma_m} {2\varepsilon^2 \widehat{C}_{inv}||v^2||_\infty-(1-\mu)(\mu-\frac{1}{2\omega} )\sigma_m^2h^2 },
\label{eqn:tau1}\\
\tau_{\varepsilon,h,2}(\mu)&: =\frac{2\varepsilon^2 h}{4C_{inv}||v||_\infty\varepsilon-(1-\mu)\sigma_mh} .\label{eqn:tau2}
\end{align}
\end{subequations}
Again the results can be expressed more compactly as $\Delta t\leq \min\left( \widehat{\tau}^\omega_{\eps,h,1}  (\mu),  \widehat{\tau}^\omega_{\eps,h,2}  (\mu)\right)$, by introducing two  extended real-valued functions
\begin{equation}
\label{eq:Gtau12}
\widehat{\tau}_{\varepsilon, \mu, i}(\mu)=\begin{cases}
                         \infty,    \quad  &\text{if}\;\frac{\varepsilon}{\sigma_m h} \leq \lambda_i(\mu)\\
		        \tau_{\varepsilon,h,i}(\mu),
		         \quad &\text{otherwise}
		        \end{cases},
		        \qquad i=1, 2. 
\end{equation}
And the scheme is unconditionally $\mu$-stable if and only if $\min\left( \widehat{\tau}^\omega_{\eps,h,1}  (\mu),  \widehat{\tau}^\omega_{\eps,h,2}  (\mu)\right)=\infty$. 
 \end{itemize}
 \end{thm}
 
We can see now that with some choice of  $\mu$, this new stability result in Theorem \ref{thm:mu} captures the unconditional stability in the diffusive regime. This regime at the discrete level is characterized by \eqref{eq:thm:stab1} and \eqref{eq:thm:stab3} when $\varepsilon/(\sigma_m h)$ is relatively small.  It is also clear that the choice of $\mu$ matters when one  interprets  the results. For instance when $k=0$, with $\mu=1/(2\omega)$, the IMEX1-LDG method is unconditionally stable in the diffusive regime, yet with $\mu=1$, we  no longer see this property according to Theorem  \ref{thm:mu}.  This motivates us to further refine the results.   Based on the definition of the (unconditional) stability in Definition \ref{def:mu-stable},  we consider an optimization problem for any given $\varepsilon, h$, and look for  the ``best'' possible choice of $\mu$, that maximizes   the unconditionally stable region (that is, to maximize the allowable range of  $\varepsilon/(\sigma_m h)$ in  \eqref{eq:thm:stab1} and \eqref{eq:thm:stab3}), and possibly  also maximizes the allowable time step condition in  \eqref{eq:thm:stab2} and \eqref{eq:thm:stab4}  when the method is conditionally stable. 
The optimization process leads to Theoreom \ref{thm:stab} that comes next, with the underlying logic as 
$$\max\{\lambda: \lambda\leq \Theta(\mu, \lambda), \forall \mu\in [\mathcal{H}(\lambda), 1] \} =\max\{\lambda:  \lambda\leq \max_{\mu\in [\mathcal{H}(\lambda), 1] }\Theta(\mu, \lambda)\},
$$
if all maximums are assumed to exist, and $\Theta, \mathcal{H}$ are some continuous functions. The relation holds if $[\mathcal{H}(\lambda), 1]$ is replaced by $(\mathcal{H}(\lambda), 1]$. 
Note that the weight function in the stability results is in the form $\omega=\omega(\varepsilon/(\sigma_m h))$.
 


\begin{thm}{\bf (Stability: $\omega>\frac{1}{2}$)}
\label{thm:stab}  When 
$\omega>\frac{1}{2}$,
 the following stability results hold for the IMEX1-LDG method, defined as  \eqref{eq:FDG:1T} with \eqref{eq:bilinear:def}-\eqref{eq:flux}. 
 
 \begin{itemize}
 %
 \item [(i)]  When $k=0$,  the IMEX1-LDG method is stable when 
 \begin{equation}
 \Delta t\leq  \Dt_{\textrm{stab}, 0} (\varepsilon, h):= \max_{\mu\in[\frac{1}{2\omega}, 1]}  \widehat{\tau}_{\eps,h,0}  (\mu)
 = \widehat{\tau}_{\eps,h,0}  \left(\frac{1}{2\omega}\right).
 \end{equation}
 In particular, the method is unconditionally stable if $ \Dt_{\textrm{stab}, 0} (\varepsilon, h)=\infty$, that is, when 
\begin{equation}
\label{eq:thm:stab1:o}
\frac{\varepsilon}{\sigma_m h}\leq  \max_{\mu\in[\frac{1}{2\omega}, 1]}  \lambda_{0}(\mu)= \lambda_{0}\left(\frac{1}{2\omega}\right)
= \frac{1-\frac{1}{2\omega}}{2||v||_\infty }.
\end{equation}
  Otherwise,
 the method is conditionally stable under the time step condition 
 \begin{align}
\label{eq:thm:stab2:o}
\Delta t\leq  \max_{\mu\in[\frac{1}{2\omega}, 1]} \tau_{\varepsilon, h, 0}(\mu)= \tau_{\varepsilon, h, 0}\left(\frac{1}{2\omega}\right)
=\frac{2\varepsilon^2 h}{2 ||v||_\infty \varepsilon-(1-\frac{1}{2\omega})\sigma_m h}.
\end{align}

\item [(ii)]  
When $1\leq k \leq 9$,  the IMEX1-LDG method is stable when 
 \begin{equation}
 \Delta t\leq  \Dt_{\textrm{stab}} (\varepsilon, h):= \max_{\mu\in(\frac{1}{2\omega}, 1]} \min\left(  \widehat{\tau}^\omega_{\eps,h,1}  (\mu),  \widehat{\tau}^\omega_{\eps,h,2}  (\mu) \right).
 \end{equation}
 In particular, the method is unconditionally stable if $ \Dt_{\textrm{stab}} (\varepsilon, h)=\infty$, that is
when 
\begin{align}
\frac{\varepsilon}{\sigma_mh}&\leq
\max_{\mu\in (\frac{1}{2\omega}, 1]} \min\left( \lambda_{1} (\mu),  \lambda_{2} (\mu)\right)=
\min\left( \lambda_{1} (\mu),  \lambda_{2} (\mu)\right)|_{\mu=\mu_\star}\notag\\
&=
  \lambda_\star:=\frac{2(1-\frac{1}{2\omega} )C_{inv}||v||_\infty}{\widehat{C}_{inv}||v^2||_\infty+8(C_{inv}||v||_\infty)^2}.
\label{eq:thm:stab3:o}
\end{align}
Otherwise the method is conditionally stable under the time step condition
\begin{align}
\Delta t&\leq
\max_{\mu\in (\frac{1}{2\omega}, 1]}\min\left( \widehat{\tau}^\omega_{\eps,h,1}  (\mu),  \widehat{\tau}^\omega_{\eps,h,2}  (\mu)\right)\notag\\
&= \tau_{\varepsilon, h, 1}( \min(\mu_S(\frac{\varepsilon}{\sigma_m h}), 1)\notag\\
&=\left\{
\begin{array}{ll}
\tau_{\varepsilon, h, 1}(\mu_S(\frac{\varepsilon}{\sigma_m h}))=
\frac{4C_{inv}||v||_\infty\varepsilon^2 h}{(8(C_{inv}||v||_\infty)^2+\widehat{C}_{inv}||v^2||_\infty)\varepsilon-2C_{inv}||v||_\infty(1-\frac{1}{2\omega})\sigma_m h}, & \text{for}\;\lambda_\star<\frac{\varepsilon}{\sigma_m h}\leq \widehat{\lambda}_\star,\\
\tau_{\varepsilon, h, 1}(1)=\frac{(1-\frac{1}{2\omega})\sigma_m h^2}{\widehat{C}_{inv}||v^2||_\infty},
			&\text{for}\; \frac{\varepsilon}{\sigma_m h} > \widehat{\lambda}_\star.			
\end{array}
\right.
\label{eq:thm:stab4:o}
\end{align}
Here 
\begin{subequations}
\begin{align}
\mu_\star&=\frac{1+\frac{1}{2\omega}{\mathcal K} }{1+{\mathcal K}}=\frac{\widehat{C}_{inv}||v^2||_\infty+4(C_{inv}||v||_\infty)^2/\omega}{\widehat{C}_{inv}||v^2||_\infty+8(C_{inv}||v||_\infty)^2},\label{eq:thm:stab5:o:1}\\
\mu_S(\lambda)&=\frac{1}{2\omega}+\frac{1}{2}\lambda\frac{\widehat{C}_{inv}||v^2||_\infty}{C_{inv}||v||_\infty}, \quad
\widehat{\lambda}_\star=\mu_S^{-1}(1)=2(1-\frac{1}{2\omega})\frac{C_{inv}||v||_\infty}{\widehat{C}_{inv}||v^2||_\infty}.\label{eq:thm:stab5:o:3}
\end{align}
\end{subequations}
 \end{itemize}
  \end{thm}
  
  


\begin{rem} The results in Theorem \ref{thm:stab} also implies an alternative route to obtain this theorem. In fact,  one can establish  Theorem \ref{thm:stab} by following the proof of Theorem \ref{thm:mu} and taking $\mu=\frac{1}{2\omega}$ when $k=0$, and taking 
\begin{align}
&\mu=\mu(\varepsilon, h; k):
=\left\{
\begin{array}{ll}
\mu_\star, & \text{for}\;\frac{\varepsilon}{\sigma_m h}\leq \lambda_\star\\
\min\left(\mu_S(\frac{\varepsilon}{\sigma_m h}),1\right), & \text{for}\;\frac{\varepsilon}{\sigma_m h}>\lambda_\star,\\
\end{array}
\right.
\end{align}
in  defining the discrete energy $E_{h,\mu}^n$ in \eqref{eqn:energy-mu}, tailored for each given $\varepsilon$, $h$ (implicitly also for a given  weight function $\omega(\varepsilon/(\sigma_m h))$.
Note that $\mu$ is chosen according to $\varepsilon/(\sigma_m h)$ that describes the regime the model is in with respect to the 
discretization parameter $h$. The assumption $1\leq k\leq 9$ in this theorem is to ensure ${\mathcal K}>1$, see Lemma \ref{lem:mK}.
\end{rem}


Following the notion of the stability in Definition \ref{def:mu-stable} and with  $E_{h, 1}^n=E_h^n$, we can combine the results in Theorem \ref{thm:stab_nmu} and Theorem \ref{thm:stab}, and obtain our final results on numerical stability for a general weight function $\omega=\omega(\varepsilon/(\sigma_m h))$ that satisfies the property \eqref{w:prop:0}.


\begin{thm}\label{thm:stab_final}
The  following stability results hold for the IMEX1-LDG method, defined as  \eqref{eq:FDG:1T} with \eqref{eq:bilinear:def}-\eqref{eq:flux}. 

 \begin{itemize}
 \item[(i)] 
When $k=0$, the method is unconditionally stable, if 
\begin{equation}
\omega>\frac{1}{2}\quad \text{and}\quad \frac{\varepsilon}{\sigma_m h}\leq \frac{1-\frac{1}{2\omega}}{2||v||_\infty}.
\label{eq:thm:stab1:f}
\end{equation}
Otherwise, the method is conditionally stable under the time step condition
\begin{equation}
\Delta t\leq \max\left(
\frac{2||v||_\infty \varepsilon h+\sigma_m h^2}{ 2||v||_\infty (||v||_\infty+\langle |v|\rangle)},\;
 \frac{ 2\varepsilon^2 h \cdot{\bf 1}_{\{\omega>\frac{1}{2}\}}}{2 ||v||_\infty \varepsilon-(1-\frac{1}{2\omega})\sigma_m h}
\right).
\label{eq:thm:stab2:f}
\end{equation}


\item[(ii)] When $1\leq k \leq 9$,  the method is unconditionally stable, if 
\begin{equation}
\omega>\frac{1}{2}\quad \text{and}\quad
\frac{\varepsilon}{\sigma_m h}\leq  \lambda_\star.
\label{eq:thm:stab3:f}
\end{equation}
Otherwise, the method is conditionally stable under the time step condition
\begin{equation}
\Delta t\leq \max\left(  
\frac{h}{\alpha_1+\alpha_2\alpha_3} (\sigma_m h+\min(\varepsilon, \frac{\alpha_2h}{\alpha_1})\alpha_3),\;
 {\bf{1} }_{\{\omega>\frac{1}{2}\} }\cdot \tau_{\varepsilon,h,1}\big(\min(\mu_S(\frac{\varepsilon}{\sigma_m h}),1) \big)
 \right),
 \label{eq:thm:stab4:f}
\end{equation}
where $\alpha_i, i=1, 2, 3$ are given in \eqref{eq:alpha}.
\end{itemize}
\end{thm}



\begin{rem} 
\label{rem:stab:w}  When $k=0$, the IMEX1-LDG method, denoted as IMEX1-LDG1 method,  will be of first order in both space and time. We here will examine more  explicitly the stability results for this first order method   when the model is the telegraph equation (referred to as T model) and the one-group transport equation in slab geometry (referred to as OG model). Note that $\langle |v|\rangle =1$ for the former, and $\langle |v|\rangle =\frac{1}{2}$ for the latter.  Particularly, we want to give the results for  three  weight functions, including $\omega\equiv 1$ and $\omega=\exp(-\frac{\varepsilon}{\sigma_m h})$ (used in \cite{peng2018stability}), and a piecewise-defined $\omega$ takings value $1$ for ``relatively small'' $\varepsilon$ and $0$ for large $\varepsilon$ (used in \cite{boscarino2014high}).  Our analysis will provide some guidance on how to define such  piecewise constant $\omega$. All three examples of $\omega$ are monotonically non-increasing in $\varepsilon/(\sigma_m h)$.
First of all, for the IMEX1-LDG1 method, the result \eqref{eq:thm:stab2:f} is indeed   
\begin{equation}
\Delta t\leq \max\left(
\frac{2\varepsilon h+\sigma_m h^2}{\beta},\;
 \frac{ 2\varepsilon^2 h \cdot{\bf 1}_{\{\omega>\frac{1}{2}\}}}{2 \varepsilon-(1-\frac{1}{2\omega})\sigma_m h}
\right),
\quad
\beta=\left\{
\begin{array}{ll}
4 & \textrm{(T model)}\\
3 & \textrm{(OG model)}
\end{array}
\right..
\label{eq:thm:stab5:f}
\end{equation} 


\begin{itemize}
\item[i.)]  We first consider $\omega\equiv 1$. It is easy to verify that $ \frac{ 2\varepsilon^2 h}{2 \varepsilon-(1-\frac{1}{2\omega})\sigma_m h}\Big|_{\omega=1}\geq  \frac{2\varepsilon h+\sigma_m h^2}{\beta}$ always holds.  Then the stability results for the IMEX1-LDG1 method in  \eqref{eq:thm:stab1:f}-\eqref{eq:thm:stab2:f} become: the method is unconditionally stable when $\varepsilon/(\sigma_m h)\leq 1/4$, otherwise it is conditionally stable under the time step condition 
$\Delta t\leq \frac{4\varepsilon^2 h}{4\varepsilon-\sigma_m h}.$
Note that this stability condition 
is the same for both T and OG models, and is used in  \cite{peng2018stability} for numerical experiments. 




\item[ii.)] We next consider a piecewise constant $\omega$, taking value either $1$ or $0$. To have the largest possible unconditional  stability region, our analysis suggests 
$
\omega={\bf 1}_{\{\varepsilon/(\sigma_m h)\leq 1/4\}},
$
and the respective stability results for the IMEX1-LDG1 method become: the method is unconditionally stable when $\varepsilon/(\sigma_m h)\leq 1/4$, and it is conditionally stable when
\begin{equation}
\Delta t\leq \frac{2\varepsilon h+\sigma_m h^2}{\beta}.
\label{eq:thm:stab7:f}
\end{equation}
Note when $\omega=0$, our IMEX1-LDG1 method is just the DG1-IMEX1 method in \cite{AP1, jang2014analysis}, with \eqref{eq:thm:stab7:f} as the respective time step condition for stability. The results imply that, if we apply the IMEX1-LDG1 method with $\omega=1$ in the relatively diffusive regime, namely $\varepsilon/(\sigma_m h)\leq 1/4$, and apply the DG1-IMEX1 method elsewhere, the stability condition  will be inherited from the method used in each regime.
\item[iii.)]  The final case is for $\omega=\exp(-\varepsilon/(\sigma_m h))$. Note that  $\omega>1/2$ is equivalent to $\varepsilon/(\sigma_m h)<r_*$ with $r_*=\ln(2)\approx 0.69314718$, and  the second inequality in \eqref{eq:thm:stab1:f} is equivalent to $\varepsilon/(\sigma_m h)\leq r_\dag$, where $r_{\dag}\approx  0.19589899$ is the root of $x=(2-e^x)/4$.  While the stability results in \eqref{eq:thm:stab1:f}-\eqref{eq:thm:stab2:f} are straightforward when $\varepsilon/(\sigma_m h)\leq r_\dag$ and when $\varepsilon/(\sigma_m h)\geq r_*$, the results when $\varepsilon/(\sigma_m h)\in (r_\dag, r_*)$ would depend on the model.  
With some calculation, one can obtain the stability results for the 
IMEX1-LDG1 method with this weight function,
\begin{equation}
\textrm{T model}: \quad \Delta t\leq \left\{
\begin{array}{ll}
\infty & \text{when} \; \varepsilon/(\sigma_m h)\leq r_\dag\\
 \frac{ 2\varepsilon^2 h }{2 \varepsilon-\Big(1-\exp(\varepsilon/(\sigma_m h))/2\Big)\sigma_m h} &  \text{when} \; \varepsilon/(\sigma_m h)\in (r_\dag, r_*)\\
(2\varepsilon h+\sigma_m h^2)/4 & \text{when} \; \varepsilon/(\sigma_m h)\geq r_*
\end{array}
\right., 
\end{equation}
\begin{equation}
\textrm{OG model}: \quad \Delta t\leq \left\{
\begin{array}{ll}
\infty & \text{when} \; \varepsilon/(\sigma_m h)\leq r_\dag\\
 \frac{ 2\varepsilon^2 h }{2 \varepsilon-\Big(1-\exp(\varepsilon/(\sigma_m h))/2\Big)\sigma_m h} &  \text{when} \; \varepsilon/(\sigma_m h)\in (r_\dag, r_\circ)\\
(2\varepsilon h+\sigma_m h^2)/3 & \text{when} \; \varepsilon/(\sigma_m h)\geq r_\circ
\end{array}
\right. . 
\end{equation}
Here $r_\circ\approx 0.38161849$ is the root of $(2x+1)/3=2x^2/\Big(2x-1+\exp(x)/2\Big)$. 
\end{itemize}
\end{rem}








\section{Asymptotic preserving (AP) property}
\label{sec:ap}

In this section, we will state the main theorem on the AP property of the IMEX1-LDG method when the initial data is well-prepared, namely, $g+v\partial_x\rho/\sigma_s=O(\varepsilon)$ at $t=0$. The proof will be established in Section \ref{sec:ap:proofs} based on uniform stability property of the method. 
With $W=\rho, q, g, u$,  we write $W_\varepsilon|_{t=0}=W_\varepsilon^0$, $W|_{t=0}=W_0$, and denote the numerical solution at time $t^n$ as $W_ {\eps,\Dt,h} ^n$ to emphasize the dependence on $h$, $\Delta t$, $\varepsilon$.
Here
$q^0_\varepsilon=\partial_x\rho_\varepsilon^0$ and $q_0=\partial_x\rho^0$
are   weak derivatives of $\rho^0_\varepsilon$ and $\rho_0$, respectively.
The following assumptions are made in this section for the initial data and weight function $\omega$. 

\medskip
\noindent
{\bf Assumption 1} (weak convergence and being {\em well-prepared})
\begin{align}
\rho_\varepsilon^0\rightharpoonup \rho_0,& \quad\text{in}\quad L^2(\Omega_x) \quad\text{as} \quad\varepsilon\rightarrow 0, \label{weak-cvg1}\\
\langle \zeta g_\varepsilon^0\rangle\rightharpoonup \langle \zeta g_0\rangle, &\quad\text{in}\quad L^2(\Omega_x) \quad\text{as} \quad\varepsilon\rightarrow 0, \quad \forall \zeta\in L^2(\Omega_v),
\label{weak-cvg2}\\
\langle \zeta (g_\varepsilon^0+v\sigma_s^{-1}q_\varepsilon^0) \rangle\rightharpoonup 0, &\quad\text{in}\quad L^2(\Omega_x) \quad\text{as} \quad\varepsilon\rightarrow 0, \quad \forall \zeta\in L^2(\Omega_v). \label{weak-wp}
\end{align}
{\bf  Assumption 2} (boundedness of initial data)
\begin{align}
\sup_{\varepsilon}||\rho^0_\varepsilon||<\infty,\quad\sup_{\varepsilon}|||g_\varepsilon^0|||<\infty,
 \quad\text{and}
 \quad\sup_{\varepsilon}||q_\varepsilon^0||<\infty.
\end{align}
{\bf  Assumption 3} (boundedness for $\omega$)  For any $h$, there exists $\varepsilon_0(h)$, such that 
\begin{equation}
2/3 <\omega<2, \quad \forall \varepsilon<\varepsilon_0(h).
\end{equation}

The assumption for $\omega=\omega(\varepsilon/(\sigma_mh))$ is reasonable due to its property \eqref{w:prop:0}. The next theorem is our main result in terms of the AP property of the IMEX1-LDG method, defined as \eqref{eq:FDG:1T} with \eqref{eq:bilinear:def}-\eqref{eq:init0}.  



\begin{thm}
\label{thm:AP}
Let the mesh size $h$ be fixed. For any time step size $\Delta t$,
there exist unique
$\rho_ {\Dt,h} ^n, u_ {\Dt,h} ^n \in U_h^k$ and $g_ {\Dt,h} ^n\in G_h^k$ for $n\geq 0$, $q_ {\Dt,h} ^n \in U_h^k$ for $n\geq 1$, such that
\begin{subequations}
\label{eq:f:lim}
\begin{align}
&\lim_{\varepsilon\rightarrow 0}W_ {\eps,\Dt,h} ^n=W_ {\Dt,h} ^n,\quad W=\rho, q, u
\label{eq:aa:-3}\\
&\lim_{\varepsilon\rightarrow 0}\langle \zeta, g_ {\eps,\Dt,h} ^n(x,\cdot)\rangle =\langle \zeta, g_ {\Dt,h} ^n(x, \cdot)\rangle,\quad \forall \zeta\in L^2(\Omega_v),\;\;\forall x\in \Omega_x,
\label{eq:aa:-2}\\
&\lim_{\varepsilon\rightarrow 0}\langle \zeta, (g_ {\eps,\Dt,h} ^n, \psi) \rangle =\langle \zeta, (g_ {\Dt,h} ^n, \psi)\rangle,\quad \forall \zeta\in L^2(\Omega_v), \;\;\forall \psi\in L^2({\Omega_x}).
\label{eq:aa:-1}
\end{align}
\end{subequations}
Furthermore, they satisfy the following scheme
\begin{subequations}
\label{eq:lim:s}
\begin{align}
\label{eq:lim1}
(q_ {\Dt,h} ^{n+1},\testQ)+d_h(\rho_ {\Dt,h} ^{n+1},\testQ)&=0,
\quad\forall\testQ\in U_h^k,\\
\label{eq:lim4}
(\sigma_s u_ {\Dt,h} ^{n+1},\testU)&=(q_ {\Dt,h} ^{n+1},\testU)
\quad\forall\testU\in U_h^k,
\\
\label{eq:lim2}
(\frac{\rho_ {\Dt,h} ^{n+1}-\rho_ {\Dt,h} ^n}{\Delta t},\testR)&=\langle v^2 \rangle  l_h(u_ {\Dt,h} ^{n+1},\testR)-(\sigma_a \rho_ {\Dt,h} ^{n+1},\testR),\quad \forall\testR\in U_h^k,\\
\pi_h(\sigma_s g_ {\Dt,h} ^{n+1})&=- v q_ {\Dt,h} ^{n+1}, \quad  g_ {\Dt,h} ^n+vu^n_ {\Dt,h} =0,
\label{eq:lim3}
\end{align}
\end{subequations}
for $n\geq 0$, with the initial data $\rho_{ {\Dt,h} }^0=\pi_h\rho_0$.
This scheme is consistent and stable for the limiting equation \eqref{eq:AP1:3.1.1}, it involves a standard LDG method in space and backward Euler method in time. Therefore the IMEX1-LDG method is AP.
When the velocity space is discrete such as $\Omega_v=\{-1, 1\}$,
 \eqref{eq:aa:-2}-\eqref{eq:aa:-1} can be replaced by a stronger form
 \begin{equation}
\lim_{\varepsilon\rightarrow 0}g_ {\eps,\Dt,h} ^n(\cdot, v)=g_ {\Dt,h} ^n(\cdot, v), \quad \forall v\in \Omega_v.
\label{eq:aa:-4}
\end{equation}
\end{thm}


\begin{rem}\label{rem:uh}
Alternative to the modal form of the LDG discretization adopted in this work, one can instead consider its nodal form \cite{hesthaven2007nodal}. Most of our analysis in this work can be extended to the resulting nodal methods, with one main difference in how the local equilibrium being satisfied as $\varepsilon\rightarrow 0$. More specifically, using the nodal form, the equations in \eqref{eq:lim:s} containing $\sigma_s$ will be replaced by their nodal counterpart, namely,
$$\sigma_s(x_*)g^{n}_ {\Dt,h} (x_*, v)=-vq^{n}_ {\Dt,h} (x_*),
\quad \sigma_s(x_*)u_ {\Dt,h} ^n(x_*)=q_ {\Dt,h} ^n(x_*),$$
where $x_*$ is any nodal point in the discretization. Besides, the absorption terms $\sigma_a \rho$ and $\sigma_a g$ can be treated explicitly in the methods, and interested readers can refer to  \cite{peng2020thesis} for more details on the impact to stability and rigorous AP property. 
\end{rem}

\section{The IMEX1-LDG scheme }
\label{sec1:scheme}

In this section, we will review the IMEX1-LDG method proposed in \cite{peng2018stability} and extend it more systematically to the model \eqref{eq:mmd:LJ0} with general material properties $\sigma_s(x)$ and $\sigma_a(x)$, both  being in $L^\infty(\Omega_x)$ and satisfying $\sigma_M\geq \sigma_s(x)\geq\sigma_m >0, \sigma_a(x)\geq 0, \forall x\in{\Omega_x}$. The boundary conditions in space are periodic, 
and the velocity variable $v$ will not be discretized.




Two examples of the model \eqref{eq:mmd:LJ0} will be examined.  One is the one-group transport equation in slab geometry. Here $\Omega_v=[-1,1]$ and the measure $\nu$ is defined as
$
\int_{\Omega_v} f d\nu=\frac{1}{2}\int_{\Omega_v} f(x,v,t) dv,
$
with $dv$ being the standard Lebesgue measure. The other is the telegraph equation  with $\Omega_v=\{-1,1\}$, and $\nu$ is  a discrete measure, given as
$
\int_{\Omega_v}f d\nu=\frac{1}{2}\left(f(x,v=1,t)+f(x,v=-1,t)\right).
$
 There is little difference in  the formulation and analysis of the IMEX1-LDG method for both examples. 
 


\subsection{Reformulation}
\label{sec1.1:reformat}

The IMEX1-LDG method is defined based on a reformulation of \eqref{eq:mmd:LJ0}, which is obtained in several steps. As the first step, we rewrite the model into its  micro-macro decomposition  \cite{liu2004boltzmann, Lemou_No_BC}. Let $L^2(\Omega_v,\nu)$ be the square integrable space in $v$, with the inner product $\langle f , g  \rangle :=  \langle fg \rangle$.  Let $\Pi$ be the $L^2$ projection onto  $\textrm{Null}({\mathcal L})=\text{Span} \{1\}$, $\mathbf{I}$  be the identify operator, and  $\rho:=  \langle f \rangle=\Pi f$  be the macroscopic density. Then $f$ can be decomposed orthogonally into $f = \rho+\varepsilon g$, with $\rho$ and $g$ satisfying 
\begin{subequations}
\label{eq:mmd:LJ}
\begin{align}
&\partial_t\rho+\partial_x \langle vg \rangle  =-\sigma_a \rho,
\label{eq:mmd:LJ:1} \\
	&\partial_t g  +\frac{1}{\varepsilon}  (\mathbf{I}-\Pi)(v\partial_xg) +\frac{1}{\varepsilon^2 }v\partial_x\rho = -\frac{\sigma_s}{\varepsilon^2 }g-\sigma_a g. \label{eq:mmd:LJ:2}
\end{align}
\end{subequations}
This is the micro-macro decomposition.
As $\varepsilon\rightarrow 0$, the equations \eqref{eq:mmd:LJ} formally become
\begin{equation}
\label{eq:limit:diff:1stOrder}
\partial_t\rho+\partial_x \langle vg \rangle  =-\sigma_a\rho,\quad  \sigma_s g=- v\partial_x\rho,
\end{equation}
which is a first order form of the limiting diffusion equation,
\begin{equation}
\partial_t\rho=\langle v^2\rangle\partial_x\left(\partial_x\rho/\sigma_s \right)-\sigma_a\rho,
\label{eq:limit:diff}
\end{equation}
equipped with the compatible  initial condition.
The relation $\sigma_s g=-v\partial_x\rho$ in \eqref{eq:limit:diff:1stOrder} will be referred to as   the {\em local equilibrium}.  For the telegraph equation, the diffusion constant is $\langle v^2\rangle=1$, while for the one-group transport equation in slab geometry,  $\langle v^2\rangle=1/3$.

As the second step,  
a weighted diffusion term, $\omega \langle v^2\rangle  \partial_x ( \partial_x\rho/\sigma_s)$, is added  to both sides of \eqref{eq:mmd:LJ:1}, leading to 
\begin{subequations}
\label{eq:mmd:DIFF}
	\begin{align}
	&\partial_t\rho+\partial_x \langle vg \rangle+\omega \langle v^2\rangle \partial_x\left( \partial_x\rho/\sigma_s\right)  =\omega \langle v^2\rangle \partial_x\left(\partial_x\rho/\sigma_s\right)-\sigma_a\rho,\\
	&\partial_t g  +\frac{1}{\varepsilon}  (\mathbf{I}-\Pi)(v\partial_xg) +\frac{1}{\varepsilon^2 }v\partial_x\rho = -\frac{\sigma_s}{\varepsilon^2 }g-\sigma_a g.
	\end{align}
\end{subequations}
Here the weight function $\omega$ is non-negative and bounded. It is {\em independent} of $x$ and can depend on $\varepsilon$, satisfying
\begin{equation}
\omega\rightarrow1, \quad \text{as} \quad \varepsilon\rightarrow 0.
\label{w:prop:0}
\end{equation}
Additional properties desired for $\omega$ in general and considered  specifically in this work will be discussed in next subsection. 
The idea of reformulating a kinetic transport model in the diffusive scaling based on adding and subtracting a diffusive term was previously used in  \cite{boscarino2013imex} and \cite{dimarco2014implicit} to remove some parabolic stiffness in designing  AP schemes. One advancement we made in \cite{peng2018stability} and here is to improve the mathematical understanding of the  desired property and the role of the weight function  $\omega$, and such advancement can guide one to choose $\omega$  in practice.




With  the auxiliary variables $q=\partial_x \rho$ and $u=q/\sigma_s$, the  system \eqref{eq:mmd:DIFF} can also be written in its first order form
\begin{subequations}
\label{eq:mmd:Q}
	\begin{align}
	& q=\partial_x \rho, \qquad 	u=q/\sigma_s,  \label{eq:mmd:Q:0}\\
	& \partial_t\rho+\partial_x \langle v(g+\omega vu) \rangle=\omega \langle v^2\rangle \partial_x u -\sigma_a\rho,
 \label{eq:mmd:Q:1} \\
	&\partial_t g  +\frac{1}{\varepsilon}  (\mathbf{I}-\Pi)(v\partial_xg) +\frac{1}{\varepsilon^2 }v\partial_x\rho = -\frac{\sigma_s}{\varepsilon^2 }g-\sigma_a g,\label{eq:mmd:Q:2} 
	\end{align}
\end{subequations}
and correspondingly  its limiting system as  $\varepsilon\rightarrow 0$ now is 
\begin{equation}
\label{eq:AP1:3.1.1}
\partial_t\rho=\langle v^2\rangle \partial_x u-\sigma_a \rho,\quad q=\partial_x \rho=\sigma_s u, \quad g=-vq/\sigma_s=-vu.
\end{equation}
The property \eqref{w:prop:0} has been  used. The introduction of $u$ is to deal with the spatially varying scattering coefficient $\sigma_s$. 
Note that the term $v\partial_x\rho$ in \eqref{eq:mmd:Q:2} can be replaced by $vq$. 



\subsection{The IMEX1-LDG scheme}
\label{sec:num_sch}

To present the scheme,  we start with some notation.
For the computational domain $\Omega_x=[x_L, x_R]$ in space, a mesh,
 $x_L=x_{\frac{1}{2}}<x_{\frac{3}{2}}<\dots<x_{N+\frac{1}{2}}=x_R$, is introduced. Let
$I_i=[{x_{i-\frac{1}{2}}}, {x_{i+\frac{1}{2}}}]$ be an element, with $x_i$ as its center and $h_i$ as its length. Set $h=\max_i h_i$. ($\Delta x$ in the introduction is just $h$ here.) For any nonnegative integer $k$, we define a finite dimensional discrete space
\begin{equation}
U_h^k=\left\{u\in L^2(\Omega_x): u|_{I_i}\in P^k(I_i), \forall i\right\},
\label{eq:DiscreteSpace:1mesh}
\end{equation}
where the local space $P^k(I)$ consists of polynomials of degree at most $k$ on $I$.  
We also introduce 
\begin{equation}
G_h^k=\left\{u(\cdot, v) \in U_h^k: \;\; \int_{\Omega_v}\int_{\Omega_x} |u(x,v)|^2dxdv<\infty\right\}.
\label{eq:DSpace:g}
\end{equation}
For a function $\phi\in U_h^k$, we write  $\phi(x^\pm)=\lim_{\Delta x\rightarrow 0^\pm} \phi(x+\Delta x)$,
and $\phi^\pm_{i+\frac{1}{2}}=\phi(x^\pm_{i+\frac{1}{2}})$. The jump and average of $\phi$ at $x_{i+\frac{1}{2}}$ are defined as 
$[\phi]_{i+\frac{1}{2}}={\phi_{i+\frac{1}{2}}^+-\phi_{i+\frac{1}{2}}^-}$ and $\{\phi\}_{i+\frac{1}{2}}=\frac12(\phi_{i+\frac{1}{2}}^++\phi_{i+\frac{1}{2}}^-)$, respectively.

The IMEX1-LDG scheme in \cite{peng2018stability} involves a LDG  discretization in space and a first order globally stiffly accurate IMEX RK scheme in time. And an  IMEX strategy is adopted so that all the terms,  which are formally dominating in the regime $\varepsilon\ll 1$, are treated implicitly. 
The IMEX1-LDG scheme for the model with a general $\sigma_s$ is based on the system \eqref{eq:mmd:Q}, and it is defined as below.
 Given  $\rho_h^{n}, \;q_h^{n},\;u_h^{n} \in U_h^k$, $g_h^{n}\in G_h^k$  
 that approximate the solution $\rho$, $q=\partial_x\rho$, $u$,
 and $g$ at $t^n$, we look for 
 $\rho_h^{n+1}, \;q_h^{n+1},\;u_h^{n+1} \in U_h^k$, $g_h^{n+1}\in G_h^k$
 at $t^{n+1}=t^n+\Delta t$, such that $\forall\;   \testQ, \testU, \testR \in U_h^k$ and $ \testG\in G_h^k$,
\begin{subequations}
\label{eq:FDG:1T}
\begin{align}
\label{eq:FDG:1T:a}
&(q_h^{n+1},\testQ)+d_h(\rho_h^{n+1},\testQ)=0,\\
\label{eq:FDG:1T:c}
&(\sigma_s u_h^{n+1},\testU)=(q_h^{n+1},\testU),\\
\label{eq:FDG:1T:b}
&\big(\frac{\rho_h^{n+1}-\rho_h^n}{\Delta t}, \testR\big)+l_h(\langle v( g_h^n+\omega vu_h^n)\rangle,\testR)=\omega \langle v^2\rangle l_h(u_h^{n+1},\testR)
-\left(\sigma_a \rho_h^{n+1},\testR\right),\\
\label{eq:FDG:1T:d}
&\big(\frac{g_h^{n+1}-g_h^n}{\Delta t},\testG\big)+\frac{1}{\varepsilon} b_{h,v}(g_h^n,\testG)-\frac{v}{\varepsilon^2} d_h(\rho_h^{n+1},\testG)
=-\frac{1}{\varepsilon^2}(\sigma_s g_h^{n+1},\testG)
-\left(\sigma_a g_h^{n+1},\testG\right).
\end{align}
\end{subequations}
Here $(\cdot, \cdot)$ is the standard inner product for $L^2(\Omega_x)$.  
The bilinear forms $d_h, l_h$, and $b_{h, v}$ are all related to discrete spatial derivatives, and defined as
\begin{subequations}
\label{eq:bilinear:def}
\begin{align}
\label{eq:dh}
d_h(\rho_h, \testQ)&=\sum_i\int_{I_i} \rho_h \partial_x\testQ dx + \sum_i \breve{\rho}_{h,{i-\frac{1}{2}}}[\testQ]_{i-\frac{1}{2}} , \\
\label{eq:lh}
l_h(u_h,\testR)&=-\sum_i \int_{I_i} u_h \partial_x\testR dx-\sum_i \hat u_{h,{i-\frac{1}{2}}} [\testR]_{i-\frac{1}{2}},
\\
\label{eq:bh}
b_{h,v}(g_h,\testG)&=(({\bf{I}}-\Pi){\mathcal D}_h(g_h; v), \testG) =({\mathcal D}_h(g_h; v) - \langle{\mathcal D}_h(g_h; v)\rangle, \testG).
\end{align}
\end{subequations}
For a given $v\in \Omega_v$, the function ${\mathcal D}_h(g_h; v)\in U_h^k$ in \eqref{eq:bh} is an upwind DG discretization of  the  transport term $v \partial_x g$. It is determined by 
\begin{align}
({\mathcal D}_h(g_h; v), \testG)
=-\sum_i\left(\int_{I_i} vg_h \partial_x\testG dx\right) - \sum_i \widetilde{(vg_h)}_{i-\frac{1}{2}}[\testG]_{i-\frac{1}{2}}, \quad \forall \psi\in U_h^k,
\label{eq:mD}
\end{align}
where $\widetilde{vg}$ is the upwind  flux,
\begin{equation}
\label{eq:vg:upwind:L-1}
\widetilde{vg}:=
\left\{
\begin{array}{ll}
v g^-,&\mbox{if}\; v>0\\
v g^+,&\mbox{if}\; v<0
\end{array}
\right.
=v\{g\}-\frac{|v|}{2}[g].
\end{equation}

The terms $\breve{\rho}$ and $\hat{u}$ in \eqref{eq:dh}-\eqref{eq:lh}  are one of the following alternating flux pair,
\begin{equation} 
\label{eq:flux}
\textrm{right-left: } \;\;\;\; \breve{\rho} = {\rho}^+,\;\;\hat{u} = u^-;\qquad
\textrm{left-right: } \;\;\;\;
\breve{\rho}={\rho}^-,\;\;\hat{u} = {u}^+. 
\end{equation}
The choice of the numerical fluxes $\breve{\rho}$ and $\hat{u}$ is important for the numerical solution to stay close to the local equilibrium when $\varepsilon\ll 1$, and it contributes to the AP property of the scheme. Similar as in standard LDG methods, the auxiliary unknowns  $q_h$  and $u_h$  can be locally represented hence eliminated in terms of $\rho_h$. 

At $t=0$, the initialization is done via the $L^2$ projection $\pi_h$ onto $U_h^k$, namely, 
\begin{equation}
\label{eq:init0}
\rho_h^0(\cdot)=\pi_h \rho(\cdot, 0), \quad g_h^0(\cdot,v)=\pi_h g(\cdot,v,0),
\quad u_h^0(\cdot,v)=\pi_h (\sigma_s^{-1} \partial_x\rho).
\end{equation}

To complete the formulation of the scheme, one needs to specify the weight function $\omega$.  In our previous work \cite{peng2018stability}, Fourier-type stability analysis suggests 
that  $\omega$ should be chosen in the form of $\omega=\omega(\frac{\varepsilon}{h},\frac{\varepsilon^2}{\Delta t})$, to preserve the intrinsic scale of the underlying model. 
In this paper, we only consider  $\omega=\omega(\varepsilon/(\sigma_m h))$,  which is independent of $\varepsilon^2/\Delta t$. Some specific examples include $\omega=\exp\big({-\varepsilon}/({\sigma_m h}) \big)$ and $\omega\equiv 1$. One can also use a piecewise constant choice $\omega={\bf 1}_{\{ \varepsilon/(\sigma_m h) \leq \alpha\}}$, with some fixed positive constant  $\alpha$, see Remark \ref{rem:stab:w} for a specific choice of  $\alpha$ recommended by our stability analysis. (Here  ${\bf 1}_D$ is an indicator function with respect to a set $D$.) Note that all these choices are non-negative and independent of $x$, satisfying \eqref{w:prop:0}.  

The next lemma states the relation of bilinear forms $d_h$ and $l_h$, and this can be verified directly.
\begin{lem}
\label{lem:1}
With either alternating flux pair in \eqref{eq:flux},  the bilinear forms $b_h$ and $l_h$ are related, 
\begin{align}
l_h(\testQ,\testR)= d_h(\testR,\testQ),  \qquad \forall \varphi, \;  \phi\in U_h^k.
\label{eq:lem1:1}
\end{align}
\end{lem}

The unique solvability of the solution to the IMEX1-LDG method is given in next  proposition, together with some properties in \eqref{eq:prop:1} that can be easily  verified.
 The key to prove the first part of the proposition is the unique solvability of the problem examined in Lemma \ref{lem:sol}. 
\begin{prop}\label{prop:sol}
The IMEX1-LDG method is uniquely solvable for any $\varepsilon\geq 0$. In addition, the solution satisfies
\begin{equation}
\langle g_h^{n} \rangle=0, \; \forall n\geq 0,\qquad 
(\sigma_s u_h^m, \testU)=-l_h(\testU, \rho_h^m), \; \forall \testU\in U_h^k, \;\; \forall m\geq 1.
\label{eq:prop:1}
\end{equation}
\end{prop}



\begin{lem}\label{lem:sol}
Given $S\in L^2(\Omega_x)$ and $\gamma_j\geq 0, j=1, 2$. Consider the following problem: look for $\rho_h, q_h, u_h\in U_h^k$, such that $ \forall \testQ, \testU, \testR\in U_h^k$,
\begin{equation}
(q_h, \testQ)+d_h(\rho_h, \testQ)=0,  \quad (\sigma_s u_h, \testU)=(q_h, \testU), \quad (\rho_h, \testR)-\gamma_1 l_h(u_h, \testR)=-\gamma_2(\sigma_a \rho_h, \testR)+(S, \testR).
\label{eq:lem10:1}
\end{equation}
Then $\rho_h, q_h, u_h$ are uniquely solvable. 
\end{lem}
\begin{proof}
We first consider the homogeneous case with $S=0$. Take  $\testQ=\testU=u_h, \testR=\rho_h$, use the relation of $d_h$ and $l_h$, we get
$$(\rho_h, \rho_h)+\gamma_1(\sigma_s u_h, u_h)+\gamma_2(\sigma_a \rho_h, \rho_h)=0. $$
With $\gamma_1, \gamma_2, \sigma_s, \sigma_a$ being non-negative, one has $\rho_h=0$, and the equations in \eqref{eq:lem10:1} further ensure $q_h=u_h=0$. 
This, in combination with the linearity of the problem  as well as that both the solution and the test function are from the same finite dimensional space $U_h^k$, implies the unique solvability of the problem with the general  source term $S$.
\end{proof}

Following the formal asymptotic analysis as in \cite{peng2018stability}, we can show the IMEX1-LDG method is AP, namely as $\varepsilon\rightarrow 0$, its limiting scheme is a consistent and stable discretization of the limiting system \eqref{eq:AP1:3.1.1}, when the initial data is well-prepared. This will be stated in Section \ref{sec:ap} and proved in Section \ref{sec:ap:proofs} once the uniform stability is available. When the initial data is not well-prepared, our scheme can adopt a similar initial fix  \cite{peng2018stability} when $n=0$ to stay AP. There is no change to numerical stability, while the AP property can be established rigorously and the details are not presented in this paper. 

