\section{Introduction} \label{sec:intro}
Fractional partial differential equations (FPDE) have emerged as a powerful tool for modeling the multiphysics and multiscale processes in numerical simulations ranging from physics~\cite{Goychuk2017,Goychuk2020,Goychuk2021}, biology~\cite{Lai2009} to quantitative finance~\cite{Coffey2004}. For example, some of the most significant and profoundly published experimental results in random flow environments, such as the cytosol and the plasma membrane of biological cells~\cite{Rubenstein2003}, crowded complex fluids and polymer solutions~\cite{Levine2001}, dense colloidal suspensions~\cite{Kremer1990} and single-file diffusion in colloidal systems~\cite{Kou2004}, are better rationalized within the fractional calculus framework. The growing number of applications of fractional derivatives in various fields of science and engineering indicate that there is a significant demand for better numerical algorithms. However, such algorithms are predominantly designed for one-dimensional (1D) problems~\cite{Murio2008}, due to the severe memory restrictions imposed by these derivatives, a challenge which we alleviate in Section~\ref{subsec:AdaptiveTime}.

The individual physics or scale components in FPDE have very different properties that are reflected in their discretization, for example, in Fractional Advection Diffusion Reaction (FADR) systems~\cite{Jannelli2022}, the discrete advection has a relatively slow (or `nonstiff') dynamics while the diffusion is a fast (or `stiff') evolving process. Implicit-Explicit (IMEX) integrators have been proposed as an attractive alternative (compared with the fully explicit or fully implicit time integration methods) where one combines the explicit (implicit) integration for the slow (fast) scale~\cite{Crouzeix1980}. In an IMEX scheme the system of equations assume the form~\cite{Ascher1995},
\begin{equation}
\frac{\partial^\alpha \Phi}{\partial t^\alpha} = {\bf f}(\Phi) + \eta {\bf g}(\Phi),
\label{eqn:IMEX}
\end{equation}
where the superscript, $\alpha$, denotes order of the fractional derivative and $\eta$ is a nonnegative parameter. In Equation~\eqref{eqn:IMEX}, the terms collected in ${\bf f}(\Phi)$ are on a slow time scale. Because they are (possibly) nonlinear, the implementation of a fully implicit integration scheme faces performance challenges, either from a poor performance of iterative solvers or from a complex Jacobian matrix associated with the problem. Therefore, it makes sense to solve this term explicitly. The terms in ${\bf g}(\Phi)$, however, are on a fast time scale, and their explicit solution may require excessively small time steps in order to maintain numerical stability. Being usually linear in nature, they can be solved implicitly without further complications. The stability of such methods is still bounded by the Courant-Friedrichs-Lewy (CFL) condition, but because the fast terms are treated implicitly, these conditions are less strict when compared to a fully explicit scheme of similar formal order of accuracy. A class of such an asymptotically stable IMEX method is introduced in Section~\ref{sec:method}.

We have limited our focus in this work, on the development, analysis and applicability of a novel family of spatiotemporal discretization method for the numerical solution of the 1D and 2D FADR systems. In the knowledge of the present authors, a detailed analysis of the numerical methods for FADR equation, in the finite difference framework, is absent. Such analysis is available for the fraction diffusion equation~\cite{Brunner2010} and the advection-diffusion equation~\cite{Khaled2005}. Section~\ref{sec:method} describes this method. Using the 1D linear FADR equation, the time-asymptotic stability analysis and the spectral analysis for the method is outlined in Section~\ref{sec:anal}. The numerical method is benchmarked using the test cases for the 2D fraction diffusion equation, proposed by Brunner~\cite{Brunner2010}, in Section~\ref{sec:MV}. The numerical results of the subdiffusive dynamics of the viscoelastic channel flow is delineated in section~\ref{sec:NS}. Section~\ref{sec:conclude} concludes with a brief discussion of the implication of these results in future studies.

\section{Methodology} \label{sec:method}
Let $\Gamma$ be a bounded domain in $\mathbb{R}^2$ with sufficiently smooth boundary $\partial \Gamma$. We present the numerical method for an initial-boundary value problem for the time-dependent FADR equation with fractional time-derivative of order $\alpha \in (0,\,1)$, as follows,
\begin{align}
&\frac{\partial^\alpha u({\bf x}, t)}{\partial t^\alpha} + K_1({\bf x}, t) \nabla \cdot u({\bf x}, t) = K_2({\bf x}, t) \nabla^2 u({\bf x}, t) + f({\bf x}, t), \quad {\bf x} \in \Gamma, \quad t \in (0, T), \nonumber \\
& u({\bf x}, 0) = u_0({\bf x}), \quad {\bf x} \in \Gamma, \nonumber \\
& \partial u({\bf x}, 0) = g_0({\bf x}), \quad {\bf x} \in \partial \Gamma, \quad t \in (0, T),
\label{eqn:method1}
\end{align}
where $\frac{\partial^\alpha u({\bf x}, t)}{\partial t^\alpha}$ denotes the Caputo fractional derivative of order $\alpha$~\cite{Podlubny1999} with respect to $t$ defined by
\begin{equation}
\frac{\partial^\alpha u({\bf x}, t)}{\partial t^\alpha} = \frac{1}{\Gamma(1-\alpha)} \int^t_0 \frac{d t'}{(t-t')^\alpha} \frac{\partial u({\bf x}, t)}{\partial t'}, \quad 0 < \alpha < 1,
\label{eqn:method2}
\end{equation}
and the operators $\nabla(\cdot)$ and $\nabla^2(\cdot)$ in equation~\eqref{eqn:method1}, are the (integer order) gradient and the Laplacian operator in $\mathbb{R}^2$. The FADR equation is related with the non-Markovian continuous-time random walk and is a model for anomalous diffusional flow-fields such as polymer melts~\cite{Chauhan2021}, flows of liquid crystals~\cite{Sircar2010} as well as biological flows including mucus~\cite{Sircar2016} and cartilage~\cite{Sircar2015}. In the next three sections, we propose the numerical method for the spatiotemporal discretization of equation~\eqref{eqn:method1}.

\subsection{Time integration} \label{subsec:Time}
The Caputo fractional time derivative in equations~(\ref{eqn:method1}, \ref{eqn:method2}) is discretized using the difference approximation, discussed in~\cite{Podlubny1999}. Suppose the time interval [0, $T$] is discretized uniformly into $n$ subintervals, define $t_k = k \Delta t, \,\, k = 0, 1,\ldots,n$, where $\Delta t  = T/n$ is the time-step. Let $u({\bf x},\,t_k)$ be the exact value of a function $u({\bf x}, t)$ at time step $t_k$. Then, the fractional derivative can be approximated as follows,
\begin{align}
\frac{\partial^\alpha u({\bf x}, t)}{\partial t^\alpha} &\approx \frac{1}{\Gamma(1-\alpha)} \sum^k_{j=0} \frac{u({\bf x},\,t_{j+1})-u({\bf x},\,t_j)}{\Delta t} \int^{(j+1)\Delta t}_{j\Delta t} \frac{d t'}{(t_{k+1}-t')^\alpha} \nonumber \\
&=\frac{1}{\Gamma(1-\alpha)} \sum^k_{j=0} \frac{u({\bf x},\,t_{j+1})-u({\bf x},\,t_j)}{\Delta t} \int^{(k-j)\Delta t}_{(k-j+1)\Delta t} t'^{-\alpha} d t' \nonumber \\
&=\frac{1}{\Gamma(1-\alpha)} \sum^k_{j=0} \frac{u({\bf x},\,t_{k+1-j})-u({\bf x},\,t_{k-j})}{\Delta t} \int^{(j+1)\Delta t}_{j\Delta t} t'^{-\alpha} d t' \nonumber \\
&=\frac{(\Delta t)^{1-\alpha}}{\Gamma(2-\alpha)} \sum^k_{j=0} \frac{u({\bf x},\,t_{k+1-j})-u({\bf x},\,t_{k-j})}{\Delta t} r^\alpha_j,
\end{align}
where the weight, $r^\alpha_j=\left[ (j+1)^{1-\alpha} - j^{1-\alpha} \right]$. Following an earlier work by the authors which utilized a standard integer-order advection-reaction-diffusion (ADR) equations~\cite{Sircar2020}, we retain `IMEX method philosophy' by explicitly discretizing the advection and the reaction term (i.~e., `$K_1\nabla \cdot u({\bf x}, t)$' term and `$f({\bf x}, t)$' term, respectively in equation~\eqref{eqn:method1}), while the diffusive term (i.~e., `$K_2({\bf x}, t) \nabla^2 u({\bf x}, t)$' term in equation~\eqref{eqn:method1}) is treated semi-implicitly, as follows,
\begin{equation}
K_2({\bf x}, t) \nabla^2 u({\bf x}, t) \approx \theta K_2({\bf x}, t) \nabla^2 u({\bf x}, t)|_{t_k} + (1-\theta) K_2({\bf x}, t) \nabla^2 u({\bf x}, t)|_{t_{k-1}}.
\label{eqn:method3}
\end{equation}
This family of time integration method (referred to as the $\theta-$method in subsequent discussion), generalizes the computationally explicit $L1$-method by Brunner~\cite{Brunner2010} as well as the fully implicit method recently proposed by Jannelli~\cite{Jannelli2022}.

\subsection{Adaptive time stepping} \label{subsec:AdaptiveTime}
When the simulation time is long, the size of `memory' in the fractional derivative approximation becomes enormously large. However, according to the `fading memory property'~\cite{Diethelm2006}, for long times, the solution of the FADR systems change more slowly than the standard integer-order ADR processes. Hence, it makes sense to employ a large time step at longer times. Let $\widetilde{u(\cdot,\, t_k)}$ be the numerical approximation for $u(\cdot,\, t_k)$. To detect this change, we define a measure between the numerical solutions of two consecutive time steps by,
\begin{equation}
\Delta u_{t_k} = \frac{\| \widetilde{u(\cdot,\, t_k)} - \widetilde{u(\cdot,\, t_{k-1})} \|_{l^2}}{|| \widetilde{u(\cdot,\, t_{k-1}) ||_{l^2}}}, \quad k = 1, \ldots, n.
\end{equation}
For some user-defined relaxation parameter $\delta$, if $\Delta u_{t_k} < \delta$, then the time spacing is geometrically increased (i.~e., $\Delta t \rightarrow 2\Delta t$), upto some prefixed value $\Delta t_{\text{max}}$.
\subsection{Spatial approximation} \label{subsec:Space}
The gradients and Laplacian in equation~\eqref{eqn:method1} are spatially approximated using the upwind difference scheme (UDS) and the second order central difference scheme (CDS), respectively. Recall that the CDS approximation of the convective terms in equation~\eqref{eqn:method1} does not model the flow-physics accurately~\cite{Sircar2020}. Furthermore, a first order upwind scheme is too diffusive, thereby necessitating the use of higher order upwind schemes. For example,
\begin{equation}
K_1({\bf x}, t) \frac{\partial u}{\partial x} \approx {K_1}_{ij} \left( \frac{u^k_{i+1,j}-u^k_{i-1,j}}{2\Delta x} \right) + q(K^+_1 u^-_x + K^-_1 u^+_x),
\end{equation}
where
\begin{align}
K^-_1 = \text{min}({K_1}_{ij}, 0) \quad &K^+_1 = \text{max}({K_1}_{ij}, 0), \nonumber\\
u^-_x = \frac{u^k_{i-2,j} - 3u^k_{i-1,j} + 3u^k_{i,j} - u^k_{i+1,j}}{3 \Delta x} \quad & u^+_x = \frac{u^k_{i-1,j} - 3u^k_{i,j} + 3u^k_{i+1,j} - u^k_{i+2,j}}{3 \Delta x},
\label{eqn:UD3}
\end{align}
where the parameter $q=0.5$ represents the third-order accurate upwind formula (UD3) and which is utilized for the interior stencil points. The use of ghost points are avoided by setting $q = 0$ for grid points immediately adjacent to the boundary, leading to a second order accurate method (UD2) at these points. Since the focus of the present work is on the development of the time integration method, we retain same the spatial approximation, outlined above, in all the subsequent test cases.

\section{Asymptotic stability and spectral analysis: linear 1D FADR equation} \label{sec:anal}
The linear 1D FADR equation~\eqref{eqn:method1} (with $K_1=c, K_2=\gamma, f(x, t) = \lambda u(x, t)$ for $x \in (-\infty, \infty)$, where $c, \gamma$ and $\lambda$ are constants specifying the advection speed, coefficient of diffusion and coefficient of reaction, respectively) serves as a model for FPDE replicating multi-scale processes.

\subsection{Asymptotic stability analysis} \label{subsec: ASA}
First, we show that the solution of the linear 1D FADR equation, coupled with periodic boundary conditions and discretized using the numerical method outlined in Section~\ref{subsec:Time}-\ref{subsec:Space}, is asymptotically stable for a limited range of CFL and Peclet numbers. The discretization of the linear 1D FADR equation, using the $\theta$-method, is as follows,
\begin{align}
&\frac{(\Delta t)^{-\alpha}}{\Gamma{(2-\alpha)}}\left\{ \sum_{j=1}^n \left(j^{(1-\alpha)}-(j-1)^{(1-\alpha)}\right) \left(u_i^{n-j+1}-u_i^{n-j}\right)\right\} + c \left(\frac{u_{i+1}^{n-1}-u_{i-1}^{n-1}}{2\Delta x}\right) + \nonumber\\
& qc \left(\frac{u_{i-2}^{n-1}-3u_{i-1}^{n-1}+3u_{i}^{n-1}-u_{i+1}^{n-1}}{3\Delta x}\right)
\nonumber \\
& = \gamma \theta \frac{u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n}}{(\Delta x)^2} + \gamma(1-\theta)\frac{u_{i+1}^{n-1}-2u_{i}^{n-1} + u_{i-1}^{n-1}}{(\Delta x)^2} + \lambda u_{i}^{n-1},
\label{eqn:step2}
 \end{align}
where, $n / i$, denote the temporal / spatial index and we set the parameter, $q=0.5$. Introducing the following non-dimensional parameters: the CFL number, $N_c$, Peclet number, $Pe$ and Damk\"{o}hler number, $Da$, where
\begin{equation}
N_c = \frac{c (\Delta t)^\alpha}{\Delta x}, \quad Pe = \frac{\gamma (\Delta t)^\alpha}{(\Delta x)^2}, \quad Da = \frac{\lambda (\Delta x)}{c},
\end{equation}
and rearranging the terms in equation~\eqref{eqn:step2}, we arrive at the following equation,
\begin{align}
&\left(1\!+\! 2 Pe \theta \Gamma{(2\!-\! \alpha)}\right) u_i^n  \! -\!   Pe \theta \Gamma{(2\!-\!\alpha)}u_{i+1}^n \! -\!
Pe\theta\Gamma{(2\!-\! \alpha)}u_{i-1}^n   =   \left(1-qN_c\Gamma{(2\!-\! \alpha)}
\right. \nonumber\\
& \left. -2 Pe(1\!-\!\theta)\Gamma{(2\!-\!\alpha)} + Da N_c\Gamma{(2-\alpha)}\right)u_i^{n-1}  +  \bigg(-N_c\frac{\Gamma{(2-\alpha)}}{2}+q N_c\frac{\Gamma{(2-\alpha)}}{3}
 \nonumber \\
& + Pe(1\!-\!\theta)\Gamma{(2\!-\!\alpha)}\bigg)u_{i+1}^{n-1} \!+\! \bigg(N_c\frac{\Gamma{(2\!-\!\alpha)}}{2} + q N_c\Gamma{(2\!-\!\alpha)} \! +\! Pe (1-\theta)\Gamma{(2-\alpha)}\!\! \bigg)
\nonumber \\
&u_{i-1}^{n-1}- q N_c\frac{\Gamma{(2-\alpha)}}{3}u_{i-2}^{n-1} + \sum_{j=2}^n \left(j^{(1-\alpha)}-(j-1)^{(1-\alpha)}\right) \left(u_i^{n-j+1}  -  u_i^{n-j}\right).
\label{eqn:step3}
\end{align}
If $\widetilde{u_i^n}$ be the (numerically) approximate solution of equation~\eqref{eqn:step3}, then the round-off error, $\varepsilon_i^n = u_i^n - \widetilde{u_i^n}\quad (i = 0, \ldots, M)$, identically satisfies the same equation. Due to periodic boundary conditions, we have,
\begin{equation}
\varepsilon_0^n=\varepsilon_M^n \quad n = 1,2,......,T.
     \label{eqn:step4}
\end{equation}
Assuming the round-off error has the following form,
\begin{align}
    \varepsilon_i^n = \xi_n e^{Ik(ih)},
    \label{eqn:step5}
\end{align}
where $\xi_n = |\varepsilon_i^n|, I=\sqrt{-1}, h=\Delta x$ and $k = \frac{2\pi l}{L}$ (l: index, and L: Spatial domain length). We substitute the above relation in the equation for the round-off error, to arrive at the following form,
\begin{equation}
\mu_1 \xi_n = \mu_2 \xi_{n-1} +\sum_{j=2}^n r^\alpha_j \left(\xi_{n-j+1}-\xi_{n-j}\right),
 \label{eqn:step6}
\end{equation}
where
\begin{align}
&\mu_1 = 1+2Pe\theta\Gamma(2-\alpha)(1-\cos(kh)), \nonumber \\
&\mu_2 = 1  -  \left(  N_cq \left(1  -  \frac{2}{3}\cos(kh)\right) + 2 Pe(1  -  \theta)(1  -  \cos(kh))  -  DaN_c  +   IN_c\sin(kh)
\right. \nonumber \\
&\left. \qquad -  \frac{N_c q}{3} \left(2e^{-Ikh}  -  e^{-2Ikh}\right)   \right)\Gamma(2  -  \alpha), \nonumber \\
&r^\alpha_j = j^{(1-\alpha)}-(j-1)^{(1-\alpha)}.
\label{eqn:step7}
\end{align}
Next we prove the result that the finite difference scheme~(\ref{eqn:step6}-\ref{eqn:step7}) is asymptotically stable, in the form of the following theorem.
\begin{theorem}
The approximate solution to the 1D linear FADR equation using the $\theta$-method~\eqref{eqn:step6} with $\alpha \in (0,\,1)$, on the finite domain $x \in [-L,\, L]$ with periodic boundary conditions is asymptotically stable for all $t \ge 0$.
\label{thm:IMEX}
\end{theorem}
\begin{proof}
It suffices for us to show that the time-amplification factor, $\xi_n$ (equation~\eqref{eqn:step6}), obey the inequality,
\begin{equation}
\xi_n \le \xi_{n-1} \le \xi_{n-2} \le \ldots \le \xi_1 \le \xi_0.
\label{eqn:step8}
\end{equation}
\begin{itemize}
\item For $n=1$ in equation~(\ref{eqn:step6},\ref{eqn:step7}), we observe the $\mu_1 \ge 1$ while $\mu_2 \le 1$ for sufficiently small $N_c$, and $Pe$, thereby leading us towards the conclusion that $\xi_1 \le \xi_0$.

\item Using induction hypothesis, we assume that $\xi_{n-1} \le \xi_{n-2} \le \ldots \le \xi_1 \le \xi_0$. Next we show that $\xi_n \le \xi_{n-1}$.

\item Since $\mu_1 \ge 1$ and $\mu_2 \le 1$, equation~\eqref{eqn:step6} can be replaced with the following inequality,
\begin{equation}
\xi_n \le \xi_{n-1} +\sum_{j=2}^n r^\alpha_j \left(\xi_{n-j+1}-\xi_{n-j}\right) \le \xi_{n-1}.
\label{eqn:step9}
\end{equation}
The second inequality in equation~\eqref{eqn:step9} is because,
\begin{enumerate}[i)]
\item $r^\alpha_j$ is positive and
\item $r^\alpha_j > r^\alpha_{j+1}$, and
\item $\xi_{n-j+1}-\xi_{n-j} \le 0$ (via the induction hypothesis),
\end{enumerate}
which completes the proof.
\end{itemize}
\end{proof}
We emphasize that while the $\theta$-method is not unconditionally stable, it is time-asymptotically stable for a restricted range of $N_c$ and $Pe$.

\subsection{Spectral analysis} \label{subsec:SA}
Although the $\theta$-method is asymptotically stable, the presence of dispersion errors (through negative group velocity and large phase speed errors) would invalidate the long time integration results~\cite{Sircar2006}. Hence we couple the result in Section~\ref{subsec: ASA} along with the toolset of spectral analysis to deduce the relevant range of the parameters, $N_c, Pe$ and $Da$, for an accurate representation of the numerical solution of the 1D FADR equation.

Using the spectral (Fourier) representation of the approximate solution of equation~\eqref{eqn:step3}, we have $\widetilde{u_i^n} = \xi'_n e^{I(ikh)}\,\,\,(I=\sqrt{-1},\, k=\frac{2\pi l}{L})$. We define the numerical time-amplification factor, $G_{num} = \frac{\widetilde{u_i^n}}{\widetilde{u_i^{n-1}}} = \frac{\xi'_n}{\xi'_{n-1}}$. Dividing equation~\eqref{eqn:step3} with $\widetilde{u_i^{n-1}}$, we arrive at the equation governing $G_{num}$,
\begin{equation}
C_0 G_{num}^n + C_1 G_{num}^{n-1}+ \sum_{j=2}^n r_j^\alpha \left(G_{num}^{n-j+1}-G_{num}^{n-j}\right) = 0,
\label{eqn:step10}
\end{equation}
where the coefficients,
\begin{align}
C_0 = &1+2Pe\theta \Gamma{(2-\alpha)}(1-\cos(kh)) \nonumber\\
C_1 = &-1+qN_c\Gamma(2-\alpha)\left(1-\frac{4}{3}\cos(kh)+\frac{1}{3}\cos(2kh)\right)\!+\!2Pe(1-\theta)\Gamma(2-\alpha)(1-
\nonumber \\
&\cos(kh)) \!-\! Da N_c\Gamma(2\!-\! \alpha) \!+\!IN_c\Gamma(2-\alpha)\left(\!\sin(kh)+\frac{2}{3}q\sin(kh)-\frac{q}{3}\sin(2kh)\right).
\label{eqn:step11}
\end{align}
We remark that the order `n' of the polynomial~\eqref{eqn:step10} is fixed at $n=75$ in subsequent analysis, since our numerical studies have shown that an increase in the polynomial order by one, shifts the contours of the spectral variables by less that $0.001\%$.

\begin{figure}[htbp]
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig1a.png}
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig1b.png}
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig1c.png}
\vskip 1pt
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig1d.png}
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig1e.png}
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig1f.png}
\vskip 1pt
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig1g.png}
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig1h.png}
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig1i.png}
\caption{Group velocity ratio contours, $V_g$, at $\theta=0.5$, $\alpha=0.9$ and (a) $Pe=0.001, Da=-0.01$, (b) $Pe=0.001, Da=0.0$, (c) $Pe=0.001, Da=0.01$, (d) $Pe=0.01, Da=-0.01$, (e) $Pe=0.01, Da=0.0$, (f) $Pe=0.01, Da=0.01$, (g) $Pe=1.0, Da=-0.01$, (h) $Pe=1.0, Da=0.0$ and (i) $Pe=1.0, Da=0.01$.}
\label{fig:Fig1}
\end{figure}

Next, transforming the exact solution of the linear 1D FADR equation, using Fourier-Laplace transform as $u(x, t) = \int \int \hat{U}(k, \omega) e^{{\it I}(kx - \omega t)}$, we arrive at the exact dispersion relation,
\begin{equation}
\omega_{exact} = \omega = I(\lambda-\gamma k^2-c I k)^{\frac{1}{\alpha}}.
\label{eqn:step12}
\end{equation}
Similarly, we have a corresponding numerical dispersion relation for the approximate equation~\eqref{eqn:step3},
\begin{equation}
\omega_{num} = I(\lambda^N-\gamma^N k^2-c^N I k)^{\frac{1}{\alpha}},
\label{eqn:step12a}
\end{equation}
where the superscript `N' denotes the corresponding numerical values of the parameters. Since
\begin{equation}
G_{num} = \frac{u(k, t+\Delta t)}{u(k, t)} = \frac{e^{I(kx-\omega_{num} (t+\Delta t))}}{e^{I(kx-\omega_{num} t)}} = e^{-I\omega_{num}\Delta t} = e^{I \beta},
\label{eqn:step13}
\end{equation}
and the numerical phase speed is given by, $c_{num} = \frac{\omega_{num}}{k}$, we find that,
\begin{equation}
c_{num} = -\frac{\beta}{k \Delta t} = -\frac{1}{k \Delta t}\tan^{-1} \left( \frac{(G_{num})_{Imag}}{(G_{num})_{Re}} \right),
\label{eqn:step14}
\end{equation}
where the subscript `Imag' / `Re' denote the imaginary / real values, respectively. Using the expression for the exact phase speed, $c_{exact} = \frac{\omega_{exact}}{k}$ and the non-dimensional parameters described in equation~\eqref{eqn:step3}, we find the ratio of the phase speeds,
\begin{equation}
\frac{c_{num}}{c_{exact}} = \frac{-I\beta}{(DaN_c-(kh)^2Pe-I(kh)N_c)^{\frac{1}{\alpha}}},
\label{eqn:step15}
\end{equation}

Finally, the expression for the numerical and the exact group velocities are given by,
\begin{align}
& {V_g}_{num} = \left[\frac{\partial}{\partial k} (\omega_{num})\right]_{Re} = \left[\frac{\partial}{\partial k} \left(\frac{\beta}{\Delta t}\right)\right]_{Re} = \left[\frac{h}{\Delta t}\frac{d\beta}{d(kh)}\right]_{Re}\nonumber\\
&{V_g}_{exact}\! =\!\! \left[\frac{\partial}{\partial k} (\omega_{exact})\right]_{Re} \!\!\! = \! \!\left[\frac{\partial}{\partial k} \left(I(\lambda\!-\!\gamma k^2 \!-\! cI k)^{\frac{1}{\alpha}}\right)\right]_{Re} \!\!\!=\! \!\left[\frac{1}{\alpha} (c\!-\!2k\gamma I)(-I\omega)^{1\!-\!\alpha}\right]_{Re}.
\label{eqn:step16}
\end{align}

Hence,
\begin{equation}
\frac{(V_g)_{num}}{(V_g)_{exact}} = \left( \frac{\alpha}{\splitfrac{N_c\left(r^{1-\alpha}(\Delta t)^{1-\alpha}\cos{(1-\alpha)\phi}\right)}{+2 kh \left(r^{1-\alpha}(\Delta t)^{1-\alpha}\sin{(1-\alpha)\phi}\right)Pe}} \right) \left(\frac{d\beta}{d(kh)}\right)\bigg|_{Re},
\label{eqn:step17}
\end{equation}
where $\omega = r e^{I (\phi+\frac{\pi}{2})}$.
\begin{figure}[htbp]
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig2a.png}
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig2b.png}
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig2c.png}
\vskip 1pt
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig2d.png}
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig2e.png}
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig2f.png}
\vskip 1pt
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig2g.png}
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig2h.png}
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig2i.png}
\caption{Group velocity ratio contours, $V_g$, at $\theta=1.0$, $\alpha=0.9$ and (a) $Pe=0.001, Da=-0.01$, (b) $Pe=0.001, Da=0.0$, (c) $Pe=0.001, Da=0.01$, (d) $Pe=0.01, Da=-0.01$, (e) $Pe=0.01, Da=0.0$, (f) $Pe=0.01, Da=0.01$, (g) $Pe=1.0, Da=-0.01$, (h) $Pe=1.0, Da=0.0$ and (i) $Pe=1.0, Da=0.01$.}
\label{fig:Fig2}
\end{figure}

The contour plots for the group velocity ratio, $V_g = \left[\frac{V_{g, num}}{V_{g, exact}}\right]_{Re}$, for Peclet numbers, $Pe = 0.001, 0.01$ and $1.0$, and Damk\"{o}hler numbers, $Da = -0.01, 0.0$ and $0.01$, are presented for two values of $\theta$, namely, $\theta=0.5$ (in figure~\ref{fig:Fig1}) and $\theta=1.0$ (in figure~\ref{fig:Fig2}) in the $(N_c, kh)$-plane. The corresponding contour plots for the absolute phase speed error, $\Delta c = |1 - \frac{c_{num}}{c_{exact}}|$ are shown in figure~\ref{fig:Fig3} and figure~\ref{fig:Fig4} for $\theta=0.5$ and $\theta=1.0$, respectively. Regions of favorable spectral properties ($V_g > 0$ and $\Delta c \le 0.1$) are highlighted with grey color.
\begin{figure}[htbp]
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig3a.png}
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig3b.png}
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig3c.png}
\vskip 1pt
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig3d.png}
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig3e.png}
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig3f.png}
\vskip 1pt
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig3g.png}
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig3h.png}
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig3i.png}
\caption{Absolute phase speed error contours, $\Delta c$ at $\theta=0.5$, $\alpha=0.9$ and (a) $Pe=0.001, Da=-0.01$, (b) $Pe=0.001, Da=0.0$, (c) $Pe=0.001, Da=0.01$, (d) $Pe=0.01, Da=-0.01$, (e) $Pe=0.01, Da=0.0$, (f) $Pe=0.01, Da=0.01$, (g) $Pe=1.0, Da=-0.01$, (h) $Pe=1.0, Da=0.0$ and (i) $Pe=1.0, Da=0.01$.}
\label{fig:Fig3}
\end{figure}

We outline the spectral properties, namely, the group velocity ratio and the absolute phase speed error, for a specific case of the fractional order, $\alpha = 0.9$. In general, we find that in the limit of vanishingly small values of $kh$, the numerical method has favorable spectral properties (i.~e., $V_g > 0$ and $\Delta c \le 0.1$ in the limit, $kh \rightarrow 0$). The region of positive group velocities ($V_g > 0$) indicate a region where the numerical solution travels in the correct direction and the numerical instabilities in the form of q-waves are avoided~\cite{Sircar2020}. Figure~\ref{fig:Fig1} and~\ref{fig:Fig2} indicate that this favorable region is restricted to smaller values of $kh$ with larger values of $Pe$ as well as with smaller values of $\theta$. While the former observation can be attributed to the fact that a stiff diffusive term (i.~e., larger `$K_2({\bf x}, t) \nabla^2 u({\bf x}, t)$' term in equation~\eqref{eqn:method1} can be correctly approximated with smaller grid-size, $h$; the latter observation can be explained through a numerical stabilization due to the implicit treatment of fast time scale term (in this case, the diffusive term). Both of these observations are in congruence with the corresponding analysis of the integer order ADR equations~\cite{Sircar2020}. Similarly, the absolute phase error contours highlight spectrally favorable region ($\Delta c \le 0.1$, equivalently the phase errors are restricted to 10\% of the exact phase speed, see figures~\ref{fig:Fig3}, \ref{fig:Fig4}) at lower values of $Pe$ and at lower values of $\theta$, indicating the fact that the phase errors in this numerical method can be reduced by introducing finer grids~\cite{Sengupta2012}.

We summarize our discussion by indicating two sources of error which are particularly perplexing in the Direct Numerical Simulations (DNS) of FADR equations. The first source of error is the existence of $q-$waves for those set of numerical parameters for which the spatiotemporal discretization is asymptotically stable. In such a situation, the $q-$waves do not attenuate and have to be eliminated by deploying an explicit filter~\cite{Visbal2002}. Another aspect of spectral error is related to the Gibbs' phenomenon which occurs as a consequence of sharp discontinuity in the solution and which causes fictitious oscillations, a problem which can be remedied using high accuracy dispersion relation preserving schemes (which has atleast shown promise in integer order partial differential equations)~\cite{Sircar2006}.
\begin{figure}[htbp]
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig4a.png}
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig4b.png}
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig4c.png}
\vskip 1pt
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig4d.png}
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig4e.png}
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig4f.png}
\vskip 1pt
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig4g.png}
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig4h.png}
\includegraphics[width=0.33\linewidth, height=0.3\linewidth]{Fig4i.png}
\caption{Absolute phase speed error contours, $\Delta c$ at $\theta=1.0$, $\alpha=0.9$ and (a) $Pe=0.001, Da=-0.01$, (b) $Pe=0.001, Da=0.0$, (c) $Pe=0.001, Da=0.01$, (d) $Pe=0.01, Da=-0.01$, (e) $Pe=0.01, Da=0.0$, (f) $Pe=0.01, Da=0.01$, (g) $Pe=1.0, Da=-0.01$, (h) $Pe=1.0, Da=0.0$ and (i) $Pe=1.0, Da=0.01$.}
\label{fig:Fig4}
\end{figure}

\section{Method validation: 2D fractional diffusion equation} \label{sec:MV}
Using the spectrally relevant parameter values discussed in section~\ref{subsec:SA}, the $\theta$-method is verified by solving the 2D fractional diffusion equation ($K_1 = f = 0$ and $K_2 = 1$ in equation~\eqref{eqn:method1}), originally proposed by Brunner~\cite{Brunner2010}, for the case of (a) zero Dirichlet boundary conditions and (b) zero Neumann conditions, over the square domain, $(x, y) \in \Gamma = [-1,\,1]^2$. The initial condition for the two test cases can be outlined as follows,
\begin{enumerate}[i)]
\item {\bf Dirichlet case:} $u_0(x, y) = \cos\left( \frac{\pi}{2} x \right) \cos\left( \frac{\pi}{2} y \right)$,
\item {\bf Neumann case:} $u_0(x, y) = \sin\left( \frac{\pi}{2} x \right) \sin \left( \frac{\pi}{2} y \right)$,
\end{enumerate}
and the exact solution for both cases is given by
\begin{equation}
u_{Exact}(x, y, t) = E_\alpha \left( -\frac{1}{4} \pi^2 t^\alpha \right) u_0(x, y),
\label{eqn:MV1}
\end{equation}
where $E_\alpha (z) = \sum^\infty_{k=0} \frac{z^k}{\Gamma(\alpha k + 1)}, \quad \alpha > 0$, is the one-parameter Mittag-Leffler function.
\begin{figure}[htbp]
\includegraphics[width=0.495\linewidth, height=0.43\linewidth]{Fig5a.png}
\includegraphics[width=0.495\linewidth, height=0.43\linewidth]{Fig5b.png}
\caption{Relative error for the solution of the 2D fractional diffusion equation at simulation time $T=0.35$, at $\alpha=1.0, \theta=1.0$ ($\star$); $\alpha=0.9, \theta=1.0$ ($\square$); $\alpha=0.67, \theta=1.0$ ($\circ$); $\alpha=0.5, \theta=1.0$ ($\triangle$); $\alpha=1.0, \theta=0.6$ ($\ast$); $\alpha=0.9, \theta=0.6$ ($\diamond$); $\alpha=0.67, \theta=0.6$ ($\bullet$); $\alpha=0.5, \theta=0.6$ ($+$); and for (a) Dirichlet boundary conditions, and (b) Neumann boundary conditions.}
\label{fig:Fig5}
\end{figure}

The domain, $\Gamma \cup \partial \Gamma$, is discretized using $51 \times 51$ points and the relative error in $l_2-norm$, $\frac{\| \widetilde{u} - u_{Exact}\|_2}{\| u_{Exact} \|_2}$, versus discretization time, $\Delta t$, is shown in figure~\ref{fig:Fig5}. The total computation time is fixed at $T=0.35$, for all cases. We note that slope of the error curves at $\alpha=1.0$ is one (highlighted with solid curves in figure~\ref{fig:Fig5}a,b). This observation can be explained since at this value of $\alpha$, the $\theta$-method reduces to the standard, integer order, Euler method, which is first order accurate. Further, note that the slope of the error curves for $\alpha < 1.0$ is sub-linear, indicating that for the fractional-order diffusion equation, the $\theta$-method has sub-linear order of accuracy.

\section{Numerical simulation: 2D fractional viscoelastic channel flow} \label{sec:NS}
Next, we describe the incompressible, subdiffusive dynamics of a planar (2D) viscoelastic channel flow for polymer melts and present the numerical solution using the $\theta$-method. In an earlier work~\cite{Chauhan2021}, the authors have derived of the mathematical model (which is a nonlinear 2D coupled FADR equation) and the linear stability of the incompressible, subdiffusive dynamics of such flows. Using the following scales for non-dimensionalizing the governing equations: the height of the channel $H$ for {\color{black}length}, the timescale $T$ corresponding to maximum base flow velocity (i.~e., $T = (H / \mathcal{U}_0)^{1/\alpha}$) for time and $\rho \mathcal{U}^2_0$ for pressure and stresses (where $\rho, \mathcal{U}_0$ is the density and the velocity scale, respectively), we summarize the model in streamfunction-vorticity formulation as follows,
\begin{subequations} \label{eqn:FullSystem}
\begin{align}
&Re \left[ \frac{\partial^\alpha \Omega}{\partial t^\alpha} + {\bf v} \cdot \nabla \Omega \right] = \nu \nabla^2 \Omega + (1-\nu) \nabla \times \nabla \cdot ({\bf A} +  {\bf F}), \label{eqn:Momentum} \\
& \nabla^2 \Psi = -\Omega, \label{eqn:Poisson} \\
&\frac{\partial^\alpha {\bf A}}{\partial t^\alpha} + {\bf v}\cdot \nabla {\bf A} - (\nabla {\bf v})^T {\bf A} -{\bf A}\nabla {\bf v} = \frac{{\bf D} - {\bf A}}{We}, \label{eqn:ExtraStress}
\end{align}
\end{subequations}
where the variables $t, \Psi, {\bf v} = (u, v) = (\frac{\partial \Psi}{\partial y}, -\frac{\partial \Psi}{\partial x}), {\bf A}, \Omega = \nabla \times {\bf v}$ denote time, streamfunction, velocity, elastic stress tensor and vorticity, respectively. ${\bf F} = \frac{\mu}{1 - \nu} ({\bf E}:{\bf E}^T)$ is the finger tensor and ${\bf E} = e^{t(\nabla {\bf v})^{\frac{1}{\alpha}}}$ is the deformation gradient tensor (the operator $(\, : \,)$ represents the tensor inner product). The term `$\left(\nabla \times \nabla \cdot {\bf F}\right)$' represents the external body force. The parameters $\eta_s, \eta_p, \eta_0 = \eta_s + \eta_p$, $\nu = \eta_s / \eta_0$ and $\mu$ are the solvent viscosity, the polymeric contribution to the shear viscosity, the total viscosity, the viscous contribution to the total viscosity of the fluid and strength of the body force, respectively. The dimensionless groups characterizing inertia and elasticity are Reynolds number, $Re = \frac{\rho \mathcal{U}_0 H}{\eta_0}$, and Weissenberg number, $We = \frac{\lambda^\alpha \mathcal{U}_0}{H}$, respectively. The parameter, $\lambda$, is the polymer relaxation time. Note that equation~\eqref{eqn:ExtraStress} represents fractional version of the regular Oldroyd-B model for viscoelastic fluids~\cite{Sircar2019}, derived in a recently reported work~\cite{Chauhan2021}. Two specific cases of the fractional derivative are considered, namely, monomer diffusion in Rouse chain melts ($\alpha = \nicefrac{1}{2}$)~\cite{Rouse1953}, and in Zimm chain solution ($\alpha = \nicefrac{2}{3}$)~\cite{Zimm1956}.

\subsection{Initial and boundary conditions}\label{subsec:IBC}
\paragraph{Initial conditions:} A rectilinear coordinate system is used with $x, y$ denoting the channel flow direction and the transverse direction, respectively. The origin of this coordinate system is chosen at the left end of the lower wall of the channel. Initially, we assume that the flow is a plane (2D) Poiseuille flow with its variation entirely in the transverse direction, namely,
\begin{equation}
{\bf U}_0 = (y - y^2) {\bf e_x},
\label{eqn:IF}
\end{equation}
where ${\bf e_x}$ is the unit vector along x-direction. The base state elastic stress tensor, ${\bf A_0} = [{A_0}_{i j}]$, is chosen as,
\begin{align}
{A_0}_{11} = 0, \nonumber \\
{A_0}_{12} = {A_0}_{21} = 1 - 2y, \nonumber \\
{A_0}_{22} = 2 We (1 - 2y)^2,
\label{eqn:IS}
\end{align}

\paragraph{Boundary conditions:} In order to imitate an infinitely long channel, periodic boundary conditions are assumed at the flow inlet and outlet. No-slip (i.~e., $u = v = 0$) and zero tangential conditions (i.~e., $\frac{\partial u}{\partial x} = \frac{\partial v}{\partial x} = 0$) are imposed on the lower wall ($y = 0$) and the upper wall ($y = 1.0$) of the channel, respectively. Further, incompressibility constraint provides an additional condition on the walls: $\frac{\partial v}{\partial y} = 0$. Since the flow is parallel to the channel walls, the walls may be treated as streamline. Thus, the streamfunction value, $\Psi$, on the wall is set as a constant. That constant (which may be different on the lower and the upper wall) is found from the no-slip condition. Zero tangential condition imply that all tangential derivatives of streamfunction vanish on the wall. Thus, the boundary condition for vorticity is found from the Poisson equation~\eqref{eqn:Poisson},
\begin{equation}
\frac{\partial^2 \Psi}{\partial y^2}|_{wall} = -\Omega_{wall}.
\label{eqn:OmegaBC}
\end{equation}
Finally, the boundary conditions for the elastic stress tensor is constructed from equation~\eqref{eqn:ExtraStress}, coupled with the no-slip and zero tangential conditions, as follows,
\begin{align}
& A_{11} = 0, \nonumber \\
& A_{12} + We \frac{\partial^\alpha A_{12}}{\partial t^\alpha} - \frac{\partial u}{\partial y} = 0, \nonumber \\
& A_{22} + We \left( \frac{\partial^\alpha A_{22}}{\partial t^\alpha} - 2 \frac{\partial u}{\partial y} A_{12}\right)= 0.
\label{eqn:ExtraStressBC}
\end{align}

\subsection{Algorithmic details}\label{subsec:AD}
The size of the domain is chosen to be $(x, y) \in [0,\,5] \times [0,\,1]$. The domain is discretized using $76 \times 51$ points, such that the discrete points are equally spaced at $\Delta x = \frac{5}{75}$ and $\Delta y = \frac{1}{50}$ (excluding the boundary points, where the Dirichlet boundary conditions are imposed) in the $x-$ and the $y-$ directions, respectively. The variable in the $\theta$-method is fixed at $\theta=1.0$. The minimum and the maximum time-step are chosen as $\Delta t_{min}=10^{-3}$ and $\Delta t_{max}=1.6 \times 10^{-2}$, respectively. The time-step in the simulation is increased adaptively using the technique outlined in Section~\ref{subsec:AdaptiveTime}. The Poisson equation~\eqref{eqn:Poisson} is iteratively solved using the Gauss-Siedal iteration technique. Assuming $(N+1)$ (or $(M+1)$) points in the $x$ (or $y$) direction, the size of the coefficient matrix is $N(M-1) \times N(M-1)$. Hence, due to size constraints, the numerical inversion of this coefficient matrix imposes severe memory restrictions. Instead, we note that the coefficient matrix has the following block diagonal structure with $(M-1)$ blocks,
\begin{equation}
\begin{pmatrix}
    D_1 & D_2 & 0 & \cdots & 0 \\
    D_2 & \ddots & \ddots & & \vdots \\
    0 & \ddots & \ddots & \ddots& 0\\
    \vdots &  & \ddots & \ddots & D_2\\
    0 & \cdots& 0 & D_2 & D_1
  \end{pmatrix}
\label{eqn:AD1}
\end{equation}
where $D_1, D_2$ are $N \times N$ tridiagonal and diagonal matrices, respectively, such that $D_1 = \text{Tridiag} (r_1, R, r_1)$ and $D_2 = \text{Diag}~(r_2)$ and $r_1 = (\Delta x)^{-2}, r_2 = (\Delta y)^{-2}$ and $R = -2(r_1 + r_2)$. Since the exact location of the non-zero entries are known, the invocation of the entire coefficient matrix in the Guass-Siedal iteration is avoided. The solution is updated by utilizing the non-zero entries of each row on a case-by-case basis. The row diagonal dominance of the matrix~\eqref{eqn:AD1} insures that the iteration converges in finite number of steps. Similarly, the vorticity equation~\eqref{eqn:Momentum} involves another Laplacian term and its discretization involves a coefficient matrix identical to the form~\eqref{eqn:AD1} (with $r_1 = -\nu \theta (dx)^{-2}, r_2 = -\nu \theta (dy)^{-2}, R = -2(r_1 + r_2) + \frac{Re (\Delta t)^{-\alpha}}{\Gamma(2-\alpha)}$). Thus, the Laplacian term in equation~\eqref{eqn:Momentum} is treated identically as above.

\subsection{Numerical results}\label{subsec:AD}
The Zimm's model~\cite{Zimm1956} predicts the (`shear rate and polymer concentration independent') viscosity of the polymer solution by calculating the hydrodynamic interaction of flexible polymers (an idea which was originally proposed by Kirkwood~\cite{Kirkwood1954}) by approximating the chains using a bead-spring setup. In contrast, the Rouse model~\cite{Rouse1953} predicts that the viscoelastic properties of the polymer chain via a generalized Maxwell model, where the elasticity is governed by a single relaxation time, which is independent of the number of Maxwell elements (or the so-called `submolecules'). Since we do not wish to study the effects of elasticity, the Weissenberg number is fixed at $We=10.0$. The value of $\mu$ is fixed at $\mu=10^{-2}$ such that the body force term in equation~\eqref{eqn:Momentum} behaves as a perturbative term for the system which is initially at steady state~(\ref{eqn:IF}, \ref{eqn:IS}). Two different values of $\nu$ are considered: $\nu=0.3$ (elastic stress dominated case) and $\nu=0.6$ (viscous stress dominated case).

\paragraph{Zimm's model:} The vorticity contours for the Zimm's model are represented via figures~\ref{fig:Fig7} and~\ref{fig:Fig9}, for the elastic stress dominated case and the viscous stress dominated case, respectively. Observe that the `flow structures' are larger in size and magnitude for the elastic stress dominated case (i.~e., compare the range of the contour values in both figures). Also, observe that the structures appearing in the Zimm's model are significantly smaller in magnitude than the Rouse model (figures~\ref{fig:Fig6} and~\ref{fig:Fig8}). Physically, the formation of these `spatiotemporal macrostructures' are associated with the entanglement of the polymer chains at microscale~\cite{Rubenstein2003}, leading to localized, non-homogeneous regions with higher viscosity. Flows with parameter values, $\nu < 0.5$ (or the elastic stress dominated case), are those associated with higher concentration of polymers per unit volume. Experiments~\cite{Fogelson2015} have shown that non-Newtonian fluids with a larger polymer concentration, have a greater tendency (for the polymer strands) to agglomerate, the so-called `over-crowding effect'~\cite{Doi1996}.
\begin{figure}[htbp]
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig7a.png}
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig7b.png}
\vskip 1pt
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig7c.png}
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig7d.png}
\vskip 1pt
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig7e.png}
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig7f.png}
\vskip 1pt
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig7g.png}
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig7h.png}
\caption{Vorticity contours for the elastic stress-dominated ($\nu=0.3$) Zimm model ($\alpha=0.67$) case, shown at parameter values, $We=10.0$ and $Re=70$ (left column) versus $Re=1000$ (right column).}
\label{fig:Fig7}
\end{figure}

Further, notice (on the right hand side column in figure~\ref{fig:Fig7} and~\ref{fig:Fig9}) the disappearance of the structures at larger Reynolds number. This observation can be attributed to the fact that the macrostructures are `washed out' of the channel at higher flow velocities. In an earlier work, we had shown how the spatiotemporal phase diagram indicated regions of `temporal stability' at larger values of $Re$~\cite{Chauhan2021}. In this work, we associate these regions of temporal stability with regions devoid of flow structures.
\begin{figure}[htbp]
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig9a.png}
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig9b.png}
\vskip 1pt
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig9c.png}
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig9d.png}
\vskip 1pt
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig9e.png}
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig9f.png}
\vskip 1pt
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig9g.png}
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig9h.png}
\caption{Vorticity contours for the viscous stress-dominated ($\nu=0.6$) Zimm model ($\alpha=0.67$) case, shown at parameter values, $We=10.0$ and $Re=70$ (left column) versus $Re=1000$ (right column).}
\label{fig:Fig9}
\end{figure}

\paragraph{Rouse model:} The Rouse model represents `thicker' fluid, or fluids with slower diffusion (figures~\ref{fig:Fig6} and~\ref{fig:Fig8}). After comparing the respective range of vorticity contours at lower as well as higher values of $Re$, we find that the macrostructures are more prominent (both in size and magnitude) in the Rouse model (in comparison with the Zimm's model). Even within the Rouse model, we find that the elastic stress dominated case ($\nu=0.3$, figure~\ref{fig:Fig6}) exhibits formation of larger structures, than the viscous stress dominated case (figure~\ref{fig:Fig8}). Again, we note that the structure formation is conspicuously absent at larger values of $Re$ (i.~e., notice the plots on the right hand side columns in figures~\ref{fig:Fig6} and~\ref{fig:Fig8}).
\begin{figure}[htbp]
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig6a.png}
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig6b.png}
\vskip 1pt
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig6c.png}
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig6d.png}
\vskip 1pt
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig6e.png}
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig6f.png}
\vskip 1pt
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig6g.png}
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig6h.png}
\caption{Vorticity contours for the elastic stress-dominated ($\nu=0.3$) Rouse model ($\alpha=0.5$) case, shown at parameter values, $We=10.0$ and $Re=70$ (left column) versus $Re=1000$ (right column).}
\label{fig:Fig6}
\end{figure}

We summarize our discussion by noting that for the selected values of the parameters, $\nu, \alpha$ and $We$: (a) The Rouse model is comparatively more unstable than the Zimm's model at low fluid inertia (or low values of $Re$), (b) The elastic stress dominated case ($\nu<0.5$ case) is comparatively more unstable than the viscous stress dominated case, and (c) Temporal stability is achieved at higher values of $Re$, irrespective of the model or the polymer concentration.

A notably `abnormal' feature in our numerical simulations is the presence of temporal stability at high inertia. While the in silico studies of the classical Oldroyd-B channel flows indicate the appearance of temporal instability for Reynolds number as low as $Re \sim 50$~\cite{Khalid2021}, temporal stability at high fluid inertia for viscoelastic flows is only recognized in experimental realizations (until now). For example, Riley~\cite{Riley1988} reported an elasticity induced flow stabilization of viscoelastic fluids coated over complaint surfaces at {\color{black}a} fairly high Reynolds number ($Re \sim 4000$). In a separate study involving ethanol gel fuels, elastic stabilization at {\color{black}a} high shear rate was attributed due to an abnormally high second normal stress difference~\cite{Nandagopalan2018}. Viscoelastic flow stabilization at higher values of $Re$, in tapered microchannels, was explained due to the presence of wall effects~\cite{Zarabadi2019}. In another in vitro study, a biofilm deacidification created a non-homogeneous environment for molecular diffusion, leading to a `subdiffusive effect' with hindered flow rates~\cite{Zarabadi2018}. These in vitro studies not only corroborate our numerical outcome, especially establishing the emergence of temporally stable region at high inertia, but also highlight the potential of FPDE in effectively capturing the flow-instability transition in subdiffusive flows.
\begin{figure}[htbp]
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig8a.png}
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig8b.png}
\vskip 1pt
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig8c.png}
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig8d.png}
\vskip 1pt
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig8e.png}
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig8f.png}
\vskip 1pt
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig8g.png}
\includegraphics[width=0.49\linewidth, height=0.3\linewidth]{Fig8h.png}
\caption{Vorticity contours for the viscous stress-dominated ($\nu=0.6$) Rouse model ($\alpha=0.5$) case, shown at parameter values, $We=10.0$ and $Re=70$ (left column) versus $Re=1000$ (right column).}
\label{fig:Fig8}
\end{figure}

\section{Conclusion} \label{sec:conclude}
This investigation addresses the development as well as the asymptotic and spectral stability of a novel class of numerical method for the spatiotemporal discretization of FPDE. Section~\ref{sec:method} presented the method, including the time integration (section~\ref{subsec:Time}) and the spatial discretization (section~\ref{subsec:Space}). Using 1D linear FADR equation, the asymptotic stability and spectral analysis was outlined in section~\ref{sec:anal}. The method was validated using the test cases for the 2D fraction diffusion equation, in section~\ref{sec:MV}. Section~\ref{sec:NS} described the numerical results of the subdiffusive dynamics of the viscoelastic channel flow. We conclude our discussion with a remark that the focus of this present work was on the development of the numerical method and not on the physical description of the subdiffusive flow dynamics. Hence, a comprehensive study on the mechanics of subdiffusive channel flow, using the numerical method developed in this work, is reported elsewhere. \vskip 5pt

\noindent \textbf{Declaration of competing interest} The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

\bibliographystyle{elsarticle-harv}


\section{Introduction}

On popular social platforms such as Twitter, Reddit or Facebook, users sometimes slightly change the spelling of a text, without changing its original semantic meaning, for various reasons.
For instance, the following sentences contain real perturbed texts (underlined):
\textit{``The \underline{democRATs} responsible for their attempted race war,"} \textit{``A fake tree burned and \underline{RepubLIEcans} are calling for,"} or \textit{``Thinking about \underline{suic1de}."}
These perturbations are often users' deliberate attempts to emphasize their opinions implicitly or explicitly (e.g., democRATs and RepubLIEcans), or
to censure offensive or sensitive wording to avoid platforms' censorship (e.g., suic1de).
Moreover, we also found that religion or nationality related perturbations  (e.g., ``muslim"$\rightarrow$``mus-lim" and ``chinese"$\rightarrow$``chi-nese") often appear in the context of   cyberbullying, racism, or hatespeech on platforms.

These human-written perturbations are a useful resource for both researchers and practitioners as they provide inductively-derived (i.e., observable) attacks by humans, as opposed to deductively-derived (i.e., hypothesized) adversarial attacks by machine learning (ML) algorithms.
Although such human-written perturbations occur frequently online, however, there has been no tool that helps discover, study, and analyze these human-written perturbations. Therefore, in this demo work, we propose {{\sc CrypText}}, a system that detects and extracts human-written perturbations from social platforms in a large scale.
The contributions of {{\sc CrypText}} can be summarized as follows:
\begin{itemize}
    \setlength\itemsep{-2pt}

    \item {{\sc CrypText}} is uniquely equipped with a dictionary of over 2M human-written tokens that are categorized into over 400K unique phonetic sounds. {{\sc CrypText}} is constantly learning new perturbations from social platforms of
    Reddit and Twitter, providing users with up-to-date perturbations. {{\sc CrypText}} also enables the users to discover perturbations from input texts and normalize or de-perturb them (i.e., correct perturbed texts).
    \item {{\sc CrypText}} also facilitates the evaluation of the robustness of textual ML models on noisy user-generated contents by perturbing inputs with realistic human-written (and not machine-generated) perturbations. {{\sc CrypText}} also helps analyze the usage patterns of human-written perturbations and their associated sentiments online.
    \item {{\sc CrypText}} is also beneficial for a variety of users and audiences including NLP researchers, linguists and social scientists. {{\sc CrypText}} will be open-sourced.
\end{itemize}

\section{Literature Review}

\subsection{Definition of Text Perturbations}

Previous research on the use of text perturbations, are only limited within the context of \textit{machine-generated} adversarial texts in the literature. Particularly, they propose different  algorithms to \textit{minimally} manipulate an input text $x$ to generate an adversarial text $x'$ that can fool a target ML model $f(\cdot)$ such that $f(x') & \neq f(x)$ and $x'$ still preserves the semantic meaning of $x$ as much as possible.
Then, we can define \textit{text perturbations} of $x$ (which are different from adversarial texts), regardless of whether they are generated by machine or written by human, as a collection of manipulated texts $P_x{=}\{x'_1, x'_2,..x'_N\}$ yet \textit{\textbf{without}} the adversarial condition (i.e., $f(x')\neq f(x)$).
From now on, we will use the notation $x'\in P_x$ to denote a perturbation (and \textit{not} an adversarial example) of $x$. Moreover, we specifically focus on character-level perturbations since they are intuitive to human.

\subsection{Machine-Generated  Perturbations}
There are several automatic character-level text manipulation algorithms proposed in the adversarial NLP literature. These manipulation strategies include swapping, deleting a character in a word (e.g., ``democrats"$\rightarrow$``demorcats", TextBugger~\cite{li2018textbugger}), replacing a character by its most probable misspell (e.g., ``republicans"$\rightarrow$``rwpublicans", TextBugger~\cite{li2018textbugger}), replacing a character by another visually similar digit or symbol (e.g., ``democrats"$\rightarrow$``dem0cr@ts" TextBugger~\cite{li2018textbugger}), an \textit{accent} (e.g., ``democrats"$\rightarrow$``d\.emocr\=ats", VIPER~\cite{VIPER}), or a \textit{homoglyph} (DeepWordBug~\cite{gao2018black}). Although these strategies are shown to be very effective in flipping textual ML models, they are hypothesized based on some known vulnerabilities of a target textual model $f(\cdot)$, and \textit{not inductively} derived from actual real-life observations. Thus, their applications are limited only to adversarial NLP settings.

\subsection{Human-Written Perturbations}\label{sec:human_written_text}
From analyzing several public text corpus, we observe that machine-generated and human-written perturbations are very different. For example, human often add another layer of meaning on words by emphasizing a part of it (e.g., ``democRATs", ``repubLIEcans"). Moreover, hyphenation are also used to perturb a word (e.g., ``mus-lim", ``vac-cine"). Different from machine, humans perturb words by using emoticons, repeating characters (e.g., ``porn"$\rightarrow$``porrrrn"), and using phonetically similar characters (e.g., ``depression"$\rightarrow$``depresxion"). Human's manipulation strategies can also be very creative and require some implicit context to understand (e.g., ``democRATS"). Often, human-written perturbations also do not manifest any fixed rules. Hence, it is very challenging for an algorithm to systematically generate all human-written perturbations. As these strategies become more prevalent online, conventional ML models that often assume clean English input texts will face difficulties with noisy texts in real-life, especially models that deal with controversial topics (e.g., politics, vaccine mandate), or abusive contents (e.g., cyberbullying, hate-speech).

\section{Demonstration of {{\sc CrypText}} System}
As the first step toward understanding human-written perturbations online, we design and propose the {{\sc CrypText}} that discovers, utilizes and monitors human-written perturbations online. {{\sc CrypText}} has four main functions: (1) looking up perturbations, (2) normalizing perturbations, (3) manipulating texts using human-written perturbations, and (4) visualizing and analyzing online contents through the lens of text perturbations. All functions are provided via a graphical-user-interface (GUI) on top of a database (DB) system
and a set of public APIs. We elaborate each functions and its use cases below.

\begin{figure}[tb!]
    \centering
    \includegraphics[width=0.4\textwidth]{lookup_small.jpg}
    \caption{\textit{Look Up} function displays human-written perturbations for the input token ``amazon" in a form of an interactive 3D spherical word-cloud.}
    \label{fig:lookup}
    \vspace{-10pt}
\end{figure}

\subsection{Human-Written Token Database}
{{\sc CrypText}} relies on the database of human-written tokens of both correctly-spelled English words and their potential perturbations. This database stores all of the raw case-sensitive tokens found from several public datasets. Most of them are originally curated for detecting abusive materials online (e.g., rumors~\cite{kochkina2018all}, hatespeech~\cite{gomez2020exploring}, cyberbullying~\cite{wulczyn_thain_dixon_2017}) and often contain many perturbations. To curate the database, we first tokenize each sentence in the corpus, and encode their sounds using a customized version of the \textsc{Soundex} algorithm~\cite{stephenson1980methodology}.

Given a token, \textsc{Soundex} first fixes the first character, then matches the remaining characters one-by-one following a pre-defined rule set (e.g., \{b, f, p, v\}{$\rightarrow$}``1"). However, the original \textsc{Soundex} algorithm cannot recognize some visually-similar character manipulations that are frequently found in human-written perturbation texts (e.g., ``l"$\rightarrow$``1", ``a"$\rightarrow$``@", ``S"$\rightarrow$``5"). Moreover, fixing the first character can result in many mismatches (e.g., both ``losbian" and ``lesbian" share the same code ``L215" yet has distinct sounds). To overcome these, we thus customize the \textsc{Soundex} algorithm to encode visually-similar characters the same. We also introduce the \textit{phonetic level} parameter $\mathbf{k}$ and fix the first $\mathbf{k}+1 (\mathbf{k}{\geq}0)$ characters as part of the resulted encodings (e.g., ``losbian"$\rightarrow$``LO245", ``lesbian"$\rightarrow$``LE245" with $\mathbf{k}{\leftarrow}$1). All the extracted data using the modified \textsc{Soundex} algorithm are then stored as several hash-map $\mathcal{H}_\mathbf{k}$ with different values of $\mathbf{k}{\leq}2$ in the database. Each key and value of $\mathcal{H}_\mathbf{k}$ is a unique \textsc{Soundex} encoding and a set of all tokens sharing the same encoding at the phonetic level $\mathbf{k}$. Table \ref{tab:example} shows an example of $\mathcal{H}_1$ extracted from a set of three sentences.

\renewcommand{\tabcolsep}{7pt}
\begin{table}[tb]
    \centering
    \footnotesize
    \begin{tabular}{lr}
        \toprule
        \textbf{Key} & \textbf{Value} \\
        \cmidrule(lr){1-2}
        TH000 & \{the, thee\} \\
        DI630 & \{dirty, dirrrty\} \\
        RE4425 & \{republicans, repubLIEcans, republic@@ns\} \\
        \bottomrule
    \end{tabular}
    \caption{Example of extracted hash-map $\mathcal{H}_1$ from a corpus of two sentences \textit{``the dirrty republicans"}, \textit{``thee dirty repubLIEcans"}, \textit{``the dirty republic@@ns"}}
    \label{tab:example}
    \vspace{-15pt}
\end{table}

\subsection{\textit{Look Up}: Discovering Text Perturbations}\label{sec:lookup}

\noindent \textbf{Functionality.} {{\sc CrypText}} allows users to query to the DB a set of perturbations $P_x$ for a specific token $x$ in the form of an interactive 3D \textit{word-cloud} (Figure. \ref{fig:lookup}). To do this, {{\sc CrypText}} relies on what we called the \textit{SMS} property to characterize text perturbations. Specifically, we define a perturbation as a token that has phonetically similar \underline{S}ound, is perceived with the same \underline{M}eaning as its original word, yet
has different \underline{S}pelling.
For example, we consider ``demokRATs" as a perturbation of ``democrats" because they spell differently yet sound the same, and human readers can perceive both with the same meaning,
especially when placed in similar contexts (e.g., \textit{``Biden belongs to the democrats"} and \textit{``Biden belongs to the demokRATs"}). However, there are no automatic mechanisms that can measure the semantic similarity between two words of different spellings, especially when misspelled tokens are often encoded as a out-of-vocabulary (OOV). To overcome this, we utilize their Levenshtein edit-distance $\mathbf{d}$ in addition to their similarity in \textsc{Soundex} encodings as a proxy for semantic similarity. Intuitively, two tokens share the same semantic meaning if they phonetically sound the same and are separated by a sufficiently small Levenshtein distance.
Given the example in Table \ref{tab:example}, a search using the query $x{\leftarrow}\mathrm{``republicans"}$ with $\mathbf{k}{=}1, \mathbf{d}{=}1$ results in $P_x\leftarrow$\{republicans, repubLIEcans\}. Figure \ref{fig:lookup} shows an example output of the \textit{Look Up} function. We manually set $\mathbf{k}{=}1, \mathbf{d}{=}3$ by default. Advanced users will be able to adjust these parameters through a provided API.

\noindent \textbf{Use Cases.} One notable use case of \textit{Look Up} is \textit{keyword enrichment}, which provides additional queries to search for sensitive topics online. Specifically, it enables  users to search for contents that otherwise could have been censured or unreachable with correctly-spelled keywords. For example, only 67\% of the tweets found by Twitter's search API from Nov. 2021 using keyword ``democrats" has negative sentiment, while that number is much higher of 87\% if a search query also includes the perturbations of ``democrats" from the \textit{Look Up} function. We also observe the same results for ``republicans" (66\% v.s. 84\%) and ``vaccine" (46\% v.s. 61\%). We also found that most tweets that contain the perturbations of ``vaccine" focus more on either pushing back vaccine mandate or raising concerns on its safety and effectiveness, compared with the ones that are not perturbed.

\subsection{\textit{Normalization}: Detecting and De-Perturbing Texts}
\noindent \textbf{Functionality}. This function detects and corrects text perturbations from an input $x$ by leveraging the SMS property described in the \textit{Look Up} function. Specifically, for each $i$-th token $x_i$ in $x$, {{\sc CrypText}} retrieves a set of potential matching English words for $x_i$ given hyper-parameter $\mathbf{k}, \mathbf{d}$. A token $w^*$ is a valid perturbation if there exists an English word that shares the same sound encoding as $x_i$ at phonetic level $\mathbf{k}$ (i.e., belong to the same hash-map $\mathcal{H}_\mathbf{k}$), within an upper-bound edit-distance $\mathbf{d}$. However, there can be several candidate English words that associate with the token $x_i$. In this case, {{\sc CrypText}} further ranks them by approximating how they fit into their surrounding local context of $x_i$ in $x$. Specifically, we utilize a large pre-trained masked language model $G$ to calculate a \textit{coherency} score. Intuitively, this score is calculated by utilizing the pre-trained power of a complex language model such as BERT~\cite{devlin2018bert} to calculate how likely $w^*$ appears in the immediate context of $x_i$. Then, {{\sc CrypText}} outputs the most probable {\em de-perturbed} candidate for each $x_i{\in}x$ on its GUI (Figure. \ref{fig:normalization}). Advanced users can also retrieve all candidates $w^*$ and their coherency scores via a provided API.
\begin{figure}[tb!]
    \centering
    \includegraphics[width=0.4\textwidth]{normalize.png}
    \caption{\textit{Normalization} function displays the normalized version of the input with corrected tokens being highlighted. The GUI also shows a popup describing a token before and after normalization.}
    \label{fig:normalization}

\end{figure}

\noindent \textbf{Use Cases.} There are two notable uses cases, namely (1) de-noising inputs of textual ML models and (2) usage as an additional predictive signal for ML pipelines. Given that the majority of ML models are often trained only on clean English corpus, {{\sc CrypText}} can be used to correct all possible human-written perturbations in the training corpus.
Moreover, the presence of perturbations within a sentence can also inform potential adversarial behaviors from its writer, especially those  offensive or controversial perturbations (e.g., cyberbullying, political discourse, vaccine mandates, spreading disinformation), as part of a ML pipeline.

\begin{figure}[tb!]
    \centering
    \includegraphics[width=0.4\textwidth]{perturbation.png}
    \caption{{{\sc CrypText}} perturbs an tweet with highlighted perturbations.}
    \label{fig:perturbation}
    \vspace{-10pt}
\end{figure}

\subsection{{\textit{Perturbation}: Text Manipulation using Human-Written Perturbations}}

\noindent \textbf{Functionality}. This function helps manipulate an input text using either case-sensitive or case-insensitive human-written perturbations at a user-specified manipulation ratio $\mathbf{r}$ (e.g., 15\%, 25\%, 50\%). To do this, given an input text $x$, {{\sc CrypText}} first randomly samples a subset of tokens in $x$ according to the selected $r$ value. Then, it iteratively replaces each of the selected tokens by a perturbation that is randomly selected from the \textit{Look Up} function's output (Figure. \ref{fig:perturbation}).

\noindent \textbf{Use Case.} Compared with other text perturbation methods in the literature (e.g., TextBugger~\cite{li2018textbugger}, VIPER~\cite{VIPER}, DeepWordBug~\cite{gao2018black}), perturbations utilized by {{\sc CrypText}} are guaranteed to be observable in human-written texts. This makes {{\sc CrypText}} a more realistic tool to evaluate the robustness of textual ML models on noisy user-generated texts. This implication is applicable not only to NLP classification models, especially those that are built for abusive content moderation, but also on models developed for other tasks such as natural language inference and Q\&A. In fact, many of the industrial APIs are not well-equipped to deal with noisy human-written texts. Figure \ref{fig:google} shows that the \textit{Perspective} toxic detection, sentiment analysis, and text categorization API from Google Cloud all suffer from texts perturbed by {{\sc CrypText}}. Especially, the popular \textit{Perspective} API's accuracy drops nearly 10\% when only 25\% of input texts are perturbed. {{\sc CrypText}} also dedicates an \textit{ML benchmark} page that frequently updates our evaluation of publicly available NLP APIs and models on noisy human-written texts.

\subsection{\textit{Social Listening}: Monitoring Human-Written Perturbations Online}
\noindent \textbf{Functionality.} This function enables users to monitor the use of human-written perturbations online, especially on social platforms such as Reddit (via \textit{PushShift API\footnote{https://github.com/pushshift/api}}). Specifically, given a list of English words, {{\sc CrypText}} first searches on the social platforms all the contents using their perturbations as queries. Then, it aggregates and displays the usage patterns of each individual perturbation in both frequency and sentiment through interactive timeline charts. We refer the readers to the attached video for the full visualization of this function.
\begin{figure}[tb]
    \centering
    \includegraphics[width=0.4\textwidth]{google_api.png}
    \caption{Accuracy of Google's NLP APIs on texts perturbed by {{\sc CrypText}}.}
    \label{fig:google}
    \vspace{-10pt}
\end{figure}
\noindent \textbf{Use Case.} As described in the use case of \textit{Look Up} function (Sec. \ref{sec:lookup}), using perturbations as queries help capture a more comprehensive view on a specific topic, especially one that is sensitive or divided. \textit{Social Listening} helps achieve this by capturing contents that are often contrast or opposite to conventional opinions, which otherwise could be censured or unreachable on non-perturbed texts. Likewise, gatekeepers of social platforms also can utilize this function for better content moderation, especially in detecting and removing abusive texts on web (e.g., cyberbullying, racism, calling for violence and terrorism),  many of which are often intentionally written with misspellings to evade automatic detection.

\subsection{An Architectural Design}
\begin{figure}[t!b]
    \centering
    \includegraphics[width=0.35\textwidth]{archi.png}
    \caption{An implementation of {\sc CrypText}.}
    \label{fig:architecture}
    \vspace{-10pt}
\end{figure}

Figure \ref{fig:architecture} illustrates the architecture of the {{\sc CrypText}} implementation. {{\sc CrypText}} interacts with users via a GUI (i.e, website) and several function APIs. The back-end of {{\sc CrypText}} is programmed with the \textit{Django} and \textit{FastAPI} framework in \textit{python} language. Its front-end utilizes several powerful visualization tools such as \textit{chart.js}, \textit{TagCloud.js}, and \textit{dataTables.js} to provide the end-users with an interactive experience. All functions of {{\sc CrypText}} are equipped with secured public APIs, allowing users to utilize \textit{Look Up}, \textit{Normalization} and \textit{Perturbation} in bulks. Accessing such APIs requires an authorization token that will be provided upon request. All of the data utilized by {{\sc CrypText}} are stored in a MongoDB database. Since some queries might take a longer time to process, a \textit{Redis cache} is adapted to temporarily store and re-use recent queried results to enhance the user experience. Furthermore, we set up a crawler that regularly collects recent tweets (via Twitter's public stream API\footnote{\url{https://developer.twitter.com/}}) to continually enrich {{\sc CrypText}}'s database with novel perturbed tokens online.

\section{Limitation and Conclusion}
Currently, the \textit{Social Listening} function is limited to Reddit data and we plan to support other platforms in future. To the best of our knowledge, {{\sc CrypText}} is the first interactive database system that helps discover and utilize human-written text perturbations online. Database provided by {{\sc CrypText}} can also help evaluate the robustness of textual ML models on noisy human-written inputs and uncover sensitive online discussions. Especially, {{\sc CrypText}} constantly updating its DB with new perturbations via Twitter's stream API.

\bibliographystyle{IEEEtran}

\section{Introduction}

On popular social platforms such as Twitter, Reddit or Facebook, users sometimes slightly change the spelling of a text, without changing its original semantic meaning, for various reasons.
For instance, the following sentences contain real perturbed texts (underlined):
\textit{``The \underline{democRATs} responsible for their attempted race war,"} \textit{``A fake tree burned and \underline{RepubLIEcans} are calling for,"} or \textit{``Thinking about \underline{suic1de}."}
These perturbations are often users' deliberate attempts to emphasize their opinions implicitly or explicitly (e.g., democRATs and RepubLIEcans), or
to censure offensive or sensitive wording to avoid platforms' censorship (e.g., suic1de).
Moreover, we also found that religion or nationality related perturbations  (e.g., ``muslim"$\rightarrow$``mus-lim" and ``chinese"$\rightarrow$``chi-nese") often appear in the context of   cyberbullying, racism, or hatespeech on platforms.

These human-written perturbations are a useful resource for both researchers and practitioners as they provide inductively-derived (i.e., observable) attacks by humans, as opposed to deductively-derived (i.e., hypothesized) adversarial attacks by machine learning (ML) algorithms.
Although such human-written perturbations occur frequently online, however, there has been no tool that helps discover, study, and analyze these human-written perturbations. Therefore, in this demo work, we propose {{\sc CrypText}}, a system that detects and extracts human-written perturbations from social platforms in a large scale.
The contributions of {{\sc CrypText}} can be summarized as follows:
\begin{itemize}
    \setlength\itemsep{-2pt}

    \item {{\sc CrypText}} is uniquely equipped with a dictionary of over 2M human-written tokens that are categorized into over 400K unique phonetic sounds. {{\sc CrypText}} is constantly learning new perturbations from social platforms of
    Reddit and Twitter, providing users with up-to-date perturbations. {{\sc CrypText}} also enables the users to discover perturbations from input texts and normalize or de-perturb them (i.e., correct perturbed texts).
    \item {{\sc CrypText}} also facilitates the evaluation of the robustness of textual ML models on noisy user-generated contents by perturbing inputs with realistic human-written (and not machine-generated) perturbations. {{\sc CrypText}} also helps analyze the usage patterns of human-written perturbations and their associated sentiments online.
    \item {{\sc CrypText}} is also beneficial for a variety of users and audiences including NLP researchers, linguists and social scientists. {{\sc CrypText}} will be open-sourced.
\end{itemize}

\section{Literature Review}

\subsection{Definition of Text Perturbations}

Previous research on the use of text perturbations, are only limited within the context of \textit{machine-generated} adversarial texts in the literature. Particularly, they propose different  algorithms to \textit{minimally} manipulate an input text $x$ to generate an adversarial text $x'$ that can fool a target ML model $f(\cdot)$ such that $f(x'){\neq}f(x)$ and $x'$ still preserves the semantic meaning of $x$ as much as possible.
Then, we can define \textit{text perturbations} of $x$ (which are different from adversarial texts), regardless of whether they are generated by machine or written by human, as a collection of manipulated texts $P_x{=}\{x'_1, x'_2,..x'_N\}$ yet \textit{\textbf{without}} the adversarial condition (i.e., $f(x')\neq f(x)$).
From now on, we will use the notation $x'\in P_x$ to denote a perturbation (and \textit{not} an adversarial example) of $x$. Moreover, we specifically focus on character-level perturbations since they are intuitive to human.

\subsection{Machine-Generated  Perturbations}
There are several automatic character-level text manipulation algorithms proposed in the adversarial NLP literature. These manipulation strategies include swapping, deleting a character in a word (e.g., ``democrats"$\rightarrow$``demorcats", TextBugger~\cite{li2018textbugger}), replacing a character by its most probable misspell (e.g., ``republicans"$\rightarrow$``rwpublicans", TextBugger~\cite{li2018textbugger}), replacing a character by another visually similar digit or symbol (e.g., ``democrats"$\rightarrow$``dem0cr@ts" TextBugger~\cite{li2018textbugger}), an \textit{accent} (e.g., ``democrats"$\rightarrow$``d\.emocr\=ats", VIPER~\cite{VIPER}), or a \textit{homoglyph} (DeepWordBug~\cite{gao2018black}). Although these strategies are shown to be very effective in flipping textual ML models, they are hypothesized based on some known vulnerabilities of a target textual model $f(\cdot)$, and \textit{not inductively} derived from actual real-life observations. Thus, their applications are limited only to adversarial NLP settings.

\subsection{Human-Written Perturbations}\label{sec:human_written_text}
From analyzing several public text corpus, we observe that machine-generated and human-written perturbations are very different. For example, human often add another layer of meaning on words by emphasizing a part of it (e.g., ``democRATs", ``repubLIEcans"). Moreover, hyphenation are also used to perturb a word (e.g., ``mus-lim", ``vac-cine"). Different from machine, humans perturb words by using emoticons, repeating characters (e.g., ``porn"$\rightarrow$``porrrrn"), and using phonetically similar characters (e.g., ``depression"$\rightarrow$``depresxion"). Human's manipulation strategies can also be very creative and require some implicit context to understand (e.g., ``democRATS"). Often, human-written perturbations also do not manifest any fixed rules. Hence, it is very challenging for an algorithm to systematically generate all human-written perturbations. As these strategies become more prevalent online, conventional ML models that often assume clean English input texts will face difficulties with noisy texts in real-life, especially models that deal with controversial topics (e.g., politics, vaccine mandate), or abusive contents (e.g., cyberbullying, hate-speech).

\section{Demonstration of {{\sc CrypText}} System}
As the first step toward understanding human-written perturbations online, we design and propose the {{\sc CrypText}} that discovers, utilizes and monitors human-written perturbations online. {{\sc CrypText}} has four main functions: (1) looking up perturbations, (2) normalizing perturbations, (3) manipulating texts using human-written perturbations, and (4) visualizing and analyzing online contents through the lens of text perturbations. All functions are provided via a graphical-user-interface (GUI) on top of a database (DB) system
and a set of public APIs. We elaborate each functions and its use cases below.

\begin{figure}[tb!]
    \centering
    \includegraphics[width=0.4\textwidth]{ICDE'2023 - Demo CrypText/lookup_small.jpg}
    \caption{\textit{Look Up} function displays human-written perturbations for the input token ``amazon" in a form of an interactive 3D spherical word-cloud.}
    \label{fig:lookup}
    \vspace{-10pt}
\end{figure}

\subsection{Human-Written Token Database}
{{\sc CrypText}} relies on the database of human-written tokens of both correctly-spelled English words and their potential perturbations. This database stores all of the raw case-sensitive tokens found from several public datasets. Most of them are originally curated for detecting abusive materials online (e.g., rumors~\cite{kochkina2018all}, hatespeech~\cite{gomez2020exploring}, cyberbullying~\cite{wulczyn_thain_dixon_2017}) and often contain many perturbations. To curate the database, we first tokenize each sentence in the corpus, and encode their sounds using a customized version of the \textsc{Soundex} algorithm~\cite{stephenson1980methodology}.

Given a token, \textsc{Soundex} first fixes the first character, then matches the remaining characters one-by-one following a pre-defined rule set (e.g., \{b, f, p, v\}{$\rightarrow$}``1"). However, the original \textsc{Soundex} algorithm cannot recognize some visually-similar character manipulations that are frequently found in human-written perturbation texts (e.g., ``l"$\rightarrow$``1", ``a"$\rightarrow$``@", ``S"$\rightarrow$``5"). Moreover, fixing the first character can result in many mismatches (e.g., both ``losbian" and ``lesbian" share the same code ``L215" yet has distinct sounds). To overcome these, we thus customize the \textsc{Soundex} algorithm to encode visually-similar characters the same. We also introduce the \textit{phonetic level} parameter $\mathbf{k}$ and fix the first $\mathbf{k}+1 (\mathbf{k}{\geq}0)$ characters as part of the resulted encodings (e.g., ``losbian"$\rightarrow$``LO245", ``lesbian"$\rightarrow$``LE245" with $\mathbf{k}{\leftarrow}$1). All the extracted data using the modified \textsc{Soundex} algorithm are then stored as several hash-map $\mathcal{H}_\mathbf{k}$ with different values of $\mathbf{k}{\leq}2$ in the database. Each key and value of $\mathcal{H}_\mathbf{k}$ is a unique \textsc{Soundex} encoding and a set of all tokens sharing the same encoding at the phonetic level $\mathbf{k}$. Table \ref{tab:example} shows an example of $\mathcal{H}_1$ extracted from a set of three sentences.

\renewcommand{\tabcolsep}{7pt}
\begin{table}[tb]
    \centering
    \footnotesize
    \begin{tabular}{lr}
        \toprule
        \textbf{Key} & \textbf{Value} \\
        \cmidrule(lr){1-2}
        TH000 & \{the, thee\} \\
        DI630 & \{dirty, dirrrty\} \\
        RE4425 & \{republicans, repubLIEcans, republic@@ns\} \\
        \bottomrule
    \end{tabular}
    \caption{Example of extracted hash-map $\mathcal{H}_1$ from a corpus of two sentences \textit{``the dirrty republicans"}, \textit{``thee dirty repubLIEcans"}, \textit{``the dirty republic@@ns"}}
    \label{tab:example}
    \vspace{-15pt}
\end{table}

\subsection{\textit{Look Up}: Discovering Text Perturbations}\label{sec:lookup}

\noindent \textbf{Functionality.} {{\sc CrypText}} allows users to query to the DB a set of perturbations $P_x$ for a specific token $x$ in the form of an interactive 3D \textit{word-cloud} (Figure. \ref{fig:lookup}). To do this, {{\sc CrypText}} relies on what we called the \textit{SMS} property to characterize text perturbations. Specifically, we define a perturbation as a token that has phonetically similar \underline{S}ound, is perceived with the same \underline{M}eaning as its original word, yet
has different \underline{S}pelling.
For example, we consider ``demokRATs" as a perturbation of ``democrats" because they spell differently yet sound the same, and human readers can perceive both with the same meaning,
especially when placed in similar contexts (e.g., \textit{``Biden belongs to the democrats"} and \textit{``Biden belongs to the demokRATs"}). However, there are no automatic mechanisms that can measure the semantic similarity between two words of different spellings, especially when misspelled tokens are often encoded as a out-of-vocabulary (OOV). To overcome this, we utilize their Levenshtein edit-distance $\mathbf{d}$ in addition to their similarity in \textsc{Soundex} encodings as a proxy for semantic similarity. Intuitively, two tokens share the same semantic meaning if they phonetically sound the same and are separated by a sufficiently small Levenshtein distance.
Given the example in Table \ref{tab:example}, a search using the query $x{\leftarrow}\mathrm{``republicans"}$ with $\mathbf{k}{=}1, \mathbf{d}{=}1$ results in $P_x\leftarrow$\{republicans, repubLIEcans\}. Figure \ref{fig:lookup} shows an example output of the \textit{Look Up} function. We manually set $\mathbf{k}{=}1, \mathbf{d}{=}3$ by default. Advanced users will be able to adjust these parameters through a provided API.

\noindent \textbf{Use Cases.} One notable use case of \textit{Look Up} is \textit{keyword enrichment}, which provides additional queries to search for sensitive topics online. Specifically, it enables  users to search for contents that otherwise could have been censured or unreachable with correctly-spelled keywords. For example, only 67\% of the tweets found by Twitter's search API from Nov. 2021 using keyword ``democrats" has negative sentiment, while that number is much higher of 87\% if a search query also includes the perturbations of ``democrats" from the \textit{Look Up} function. We also observe the same results for ``republicans" (66\% v.s. 84\%) and ``vaccine" (46\% v.s. 61\%). We also found that most tweets that contain the perturbations of ``vaccine" focus more on either pushing back vaccine mandate or raising concerns on its safety and effectiveness, compared with the ones that are not perturbed.

\subsection{\textit{Normalization}: Detecting and De-Perturbing Texts}
\noindent \textbf{Functionality}. This function detects and corrects text perturbations from an input $x$ by leveraging the SMS property described in the \textit{Look Up} function. Specifically, for each $i$-th token $x_i$ in $x$, {{\sc CrypText}} retrieves a set of potential matching English words for $x_i$ given hyper-parameter $\mathbf{k}, \mathbf{d}$. A token $w^*$ is a valid perturbation if there exists an English word that shares the same sound encoding as $x_i$ at phonetic level $\mathbf{k}$ (i.e., belong to the same hash-map $\mathcal{H}_\mathbf{k}$), within an upper-bound edit-distance $\mathbf{d}$. However, there can be several candidate English words that associate with the token $x_i$. In this case, {{\sc CrypText}} further ranks them by approximating how they fit into their surrounding local context of $x_i$ in $x$. Specifically, we utilize a large pre-trained masked language model $G$ to calculate a \textit{coherency} score. Intuitively, this score is calculated by utilizing the pre-trained power of a complex language model such as BERT~\cite{devlin2018bert} to calculate how likely $w^*$ appears in the immediate context of $x_i$. Then, {{\sc CrypText}} outputs the most probable {\em de-perturbed} candidate for each $x_i{\in}x$ on its GUI (Figure. \ref{fig:normalization}). Advanced users can also retrieve all candidates $w^*$ and their coherency scores via a provided API.
\begin{figure}[tb!]
    \centering
    \includegraphics[width=0.4\textwidth]{ICDE'2023 - Demo CrypText/normalize.png}
    \caption{\textit{Normalization} function displays the normalized version of the input with corrected tokens being highlighted. The GUI also shows a popup describing a token before and after normalization.}
    \label{fig:normalization}

\end{figure}

\noindent \textbf{Use Cases.} There are two notable uses cases, namely (1) de-noising inputs of textual ML models and (2) usage as an additional predictive signal for ML pipelines. Given that the majority of ML models are often trained only on clean English corpus, {{\sc CrypText}} can be used to correct all possible human-written perturbations in the training corpus.
Moreover, the presence of perturbations within a sentence can also inform potential adversarial behaviors from its writer, especially those  offensive or controversial perturbations (e.g., cyberbullying, political discourse, vaccine mandates, spreading disinformation), as part of a ML pipeline.

\begin{figure}[tb!]
    \centering
    \includegraphics[width=0.4\textwidth]{ICDE'2023 - Demo CrypText/perturbation.png}
    \caption{{{\sc CrypText}} perturbs an tweet with highlighted perturbations.}
    \label{fig:perturbation}
    \vspace{-10pt}
\end{figure}

\subsection{{\textit{Perturbation}: Text Manipulation using Human-Written Perturbations}}

\noindent \textbf{Functionality}. This function helps manipulate an input text using either case-sensitive or case-insensitive human-written perturbations at a user-specified manipulation ratio $\mathbf{r}$ (e.g., 15\%, 25\%, 50\%). To do this, given an input text $x$, {{\sc CrypText}} first randomly samples a subset of tokens in $x$ according to the selected $r$ value. Then, it iteratively replaces each of the selected tokens by a perturbation that is randomly selected from the \textit{Look Up} function's output (Figure. \ref{fig:perturbation}).

\noindent \textbf{Use Case.} Compared with other text perturbation methods in the literature (e.g., TextBugger~\cite{li2018textbugger}, VIPER~\cite{VIPER}, DeepWordBug~\cite{gao2018black}), perturbations utilized by {{\sc CrypText}} are guaranteed to be observable in human-written texts. This makes {{\sc CrypText}} a more realistic tool to evaluate the robustness of textual ML models on noisy user-generated texts. This implication is applicable not only to NLP classification models, especially those that are built for abusive content moderation, but also on models developed for other tasks such as natural language inference and Q\&A. In fact, many of the industrial APIs are not well-equipped to deal with noisy human-written texts. Figure \ref{fig:google} shows that the \textit{Perspective} toxic detection, sentiment analysis, and text categorization API from Google Cloud all suffer from texts perturbed by {{\sc CrypText}}. Especially, the popular \textit{Perspective} API's accuracy drops nearly 10\% when only 25\% of input texts are perturbed. {{\sc CrypText}} also dedicates an \textit{ML benchmark} page that frequently updates our evaluation of publicly available NLP APIs and models on noisy human-written texts.

\subsection{\textit{Social Listening}: Monitoring Human-Written Perturbations Online}
\noindent \textbf{Functionality.} This function enables users to monitor the use of human-written perturbations online, especially on social platforms such as Reddit (via \textit{PushShift API\footnote{https://github.com/pushshift/api}}). Specifically, given a list of English words, {{\sc CrypText}} first searches on the social platforms all the contents using their perturbations as queries. Then, it aggregates and displays the usage patterns of each individual perturbation in both frequency and sentiment through interactive timeline charts. We refer the readers to the attached video for the full visualization of this function.
\begin{figure}[tb]
    \centering
    \includegraphics[width=0.4\textwidth]{ICDE'2023 - Demo CrypText/google_api.png}
    \caption{Accuracy of Google's NLP APIs on texts perturbed by {{\sc CrypText}}.}
    \label{fig:google}
    \vspace{-10pt}
\end{figure}
\noindent \textbf{Use Case.} As described in the use case of \textit{Look Up} function (Sec. \ref{sec:lookup}), using perturbations as queries help capture a more comprehensive view on a specific topic, especially one that is sensitive or divided. \textit{Social Listening} helps achieve this by capturing contents that are often contrast or opposite to conventional opinions, which otherwise could be censured or unreachable on non-perturbed texts. Likewise, gatekeepers of social platforms also can utilize this function for better content moderation, especially in detecting and removing abusive texts on web (e.g., cyberbullying, racism, calling for violence and terrorism),  many of which are often intentionally written with misspellings to evade automatic detection.

\subsection{An Architectural Design}
\begin{figure}[t!b]
    \centering
    \includegraphics[width=0.35\textwidth]{ICDE'2023 - Demo CrypText/archi.png}
    \caption{An implementation of {\sc CrypText}.}
    \label{fig:architecture}
    \vspace{-10pt}
\end{figure}

Figure \ref{fig:architecture} illustrates the architecture of the {{\sc CrypText}} implementation. {{\sc CrypText}} interacts with users via a GUI (i.e, website) and several function APIs. The back-end of {{\sc CrypText}} is programmed with the \textit{Django} and \textit{FastAPI} framework in \textit{python} language. Its front-end utilizes several powerful visualization tools such as \textit{chart.js}, \textit{TagCloud.js}, and \textit{dataTables.js} to provide the end-users with an interactive experience. All functions of {{\sc CrypText}} are equipped with secured public APIs, allowing users to utilize \textit{Look Up}, \textit{Normalization} and \textit{Perturbation} in bulks. Accessing such APIs requires an authorization token that will be provided upon request. All of the data utilized by {{\sc CrypText}} are stored in a MongoDB database. Since some queries might take a longer time to process, a \textit{Redis cache} is adapted to temporarily store and re-use recent queried results to enhance the user experience. Furthermore, we set up a crawler that regularly collects recent tweets (via Twitter's public stream API\footnote{\url{https://developer.twitter.com/}}) to continually enrich {{\sc CrypText}}'s database with novel perturbed tokens online.

\section{Limitation and Conclusion}
Currently, the \textit{Social Listening} function is limited to Reddit data and we plan to support other platforms in future. To the best of our knowledge, {{\sc CrypText}} is the first interactive database system that helps discover and utilize human-written text perturbations online. Database provided by {{\sc CrypText}} can also help evaluate the robustness of textual ML models on noisy human-written inputs and uncover sensitive online discussions. Especially, {{\sc CrypText}} constantly updating its DB with new perturbations via Twitter's stream API.

\bibliographystyle{IEEEtran}