1. Introduction and background
The question of whether fractional derivatives can occur goes back to the origin of calculus and letters between l’Hôpital and Leibniz. While fractional calculus and associated equations have a long history only relatively recently have fractional partial differential equations been used effectively to describe physical systems. For example, it has been found to be important in anomalous diffusion [Reference Metzler and Klafter20, Reference Shlesinger and Klafter28, Reference Wang, Cherstvy, Chechkin, Thapa, Seno, Liu and Metzler30, Reference West, Grigolini, Metzler and Nonnenmacher33], amorphous materials [Reference Gu, Qing, Grebner, Wang and Schwarz10, Reference Pfister and Scher22, Reference Scher and Montroll26], porous media [Reference Benson, Schumer, Meerschaert and Wheatcraft8, Reference Benson, Wheatcraft and Meerschaert9, Reference Meerschaert, Zhang and Baeumer19], climate science [Reference Koscielny-Bunde, Bunde, Havlin, Eduardo Roman, Goldreich and Schellnhuber12], fractional quantum mechanics and optics [Reference Laskin14, Reference Longhi16, Reference Malomed17, Reference Qiu, Malomed, Mihalache, Zhu, Peng and He23], amongst others. Fractional equations using the Riesz fractional derivatives, Riesz transforms [Reference Riesz24] or fractional Laplacian [Reference Lischke, Pang, Gulian, Song, Glusa, Zheng, Mao, Cai, Meerschaert, Ainsworth and Karniadakis15] are effective tools when describing behaviour in complex systems because the Riesz fractional derivative is closely related to non-Gaussian statistics [Reference Meerschaert and Sikorskii18]. In porous media, the fractional Laplacian plays a central role and in fractional quantum mechanics the fractional Schrödinger equation is the key equation. Fractional media is ‘rough’ or multi-scale media that is neither regular nor random. Equations in multi-scale or fractional media can have fractional derivatives in any governing term [Reference West32].
With integer derivatives, one-dimensional (1D) non-linear Schrödinger (NLS) equation [Reference Ablowitz1, Reference Ablowitz and Segur3] and the Korteweg–de Vries (KdV) equation [Reference Korteweg and de Vries11] are well-known integrable equations possessing soliton solutions and an infinite set of conservation laws (cf. [Reference Ablowitz and Clarkson2, Reference Ablowitz and Segur3]). Integrable equations arise in non-linear dynamics and waves; they provide exactly solvable equations and are also an important element of Kolmogorov–Arnold–Moser theory which underlies our understanding of chaos. While in the space of possible non-linear evolution equations, integrable cases are extremely rare, nevertheless they occur frequently in application.
Recently, Ablowitz et al. found a new class of integrable fractional integrable systems; these include the integrable fractional KdV and NLS equations [Reference Ablowitz, Been and Carr5]; integrable fractional modified KdV and sine-Gordon equations [Reference Ablowitz, Been and Carr6] and also certain integrable fractional discrete/difference equations, i.e., integrable fractional discrete NLS equation [Reference Ablowitz, Been and Carr4]. These methods can used to construct N-soliton solutions in other integrable systems [Reference Weng, Zhang, Zang and Yan31, Reference Zhang, Weng and Yan34]. The key aspects in the methodology are having a general evolution equation such as a hierarchy of equations that can be expanded to fractional powers and completeness; these aspects are found on the direct scattering side and the solution obtained from inverse scattering.
In this paper, we will be concerned with a integrable or solvable fractional Burgers equation and its associated hierarchy. The classical Burgers equation is often considered to be the most elementary integrable/solvable non-trivial non-linear evolution equation. The Cole–Hopf transformation (cf. [Reference Ablowitz1]) leads directly to a linear evolution equation; unlike KdV, NLS or modified KdV, there is no need for inverse scattering. Here we investigate integrable fractional extensions of this Burgers equation and associated hierarchy. The results are explicit and the method is considerably simpler than the integrable fractional systems analysed using inverse scattering transform (IST). These equations are the simplest integrable/solvable fractional systems we are aware of.
Below we will first discuss linear fractional equations solvable by Fourier methods, then we will outline the main ideas associated with the integrable fractional KdV equation. We include this so that we can compare with the results we obtain for the integrable fractional Burgers equation and hierarchy. In the appendix, we provide more information about the integrable fractional KdV equation.
1.1. Linear evolution equations
Linear evolution equations have been studied extensively. For example, consider an equation of the form
\begin{equation}
q_{t} + \gamma(\partial_{x}) q_{x} = 0, ~~q(x,0)=q_0(x), ~~|x| \lt \infty
\end{equation} where 
$\gamma(x)$ is a polynomial and 
$q_0(x)$ vanishes rapidly as 
$|x| \to \infty$. These equations can be solved by Fourier transforms with the solution given by
\begin{equation}
q(x,t) = \frac{1}{2 \pi} \int_{-\infty}^{\infty} dk \hat{q}(k,0) e^{i k (x - \gamma(i k) t)},
\end{equation} where the Fourier transform of 
$q(x,t)$ is given by
\begin{equation}
\hat{q}(k,t) =\int_{-\infty}^{\infty} q(x,t) e^{-ikx}dx.
\end{equation} For example, with 
$\gamma(x)=x^2$, we have the linearized KdV equation
\begin{equation}
q_{t} + \partial_{x}^3 q = 0, ~~q(x,0)=q_0(x), ~~|x| \lt \infty.
\end{equation} Importantly, this methodology can be extended to fractional equations. For example, when 
$\gamma(k) = k^2 |k^2|^{\alpha}; ~~-1 \lt \alpha \lt 1 $ (here and below we take 
$\alpha \in \mathbb{R}$), we have
\begin{equation}
q_{t} + \partial_{x}^3|\partial_x^2|^{\alpha} q = 0, ~~q(x,0)=q_0(x), ~~|x| \lt \infty.
\end{equation} This is a linearized fractional KdV equation. Such fractional equations can be solved by Fourier transforms in a similar way to that of equations with polynomial 
$\gamma(k)$. This is based upon using the identification 
$k \to -i\partial_x$ in the equation and inserting the function 
$\gamma(ik)$ in the solution.
1.2. Fractional KdV hierarchy
Since the IST is the non-linear analog of Fourier transforms, we are motivated to extend IST to integrable fractional non-linear equations. We studied this issue in recent papers [Reference Ablowitz, Been and Carr4–Reference Ablowitz, Been and Carr6]. The main underlying issues are discussed below.
It was shown in [Reference Ablowitz, Kaup, Newell and Segur7] that, associated with the time-independent Schrödinger equation,
\begin{equation}
v_{xx} + \left( k^2 + q(x,t) \right) v = 0 , ~~|x| \lt \infty
\end{equation}was a class of solvable non-linear equations given by
\begin{equation}
q_{t} + \gamma(L^A) q_{x} = 0, ~~|x| \lt \infty
\end{equation}where LA is the operator
\begin{equation*}
L^{A} =-\frac{1}{4} \partial_{x}^2 - q + \frac{1}{2} q_{x} \int_{x}^{\infty} dy
\end{equation*} and 
$\gamma(k)$ is related to the dispersion relation of the linear part of the equation. Note the linear limit is 
$L^{A} \sim -\frac{1}{4} \partial_{x}^2 $. The operator L A is the adjoint to an operator that involves certain squared eigenfunctions associated with (1.6) [Reference Ablowitz, Kaup, Newell and Segur7]; see also the appendix.
The standard KdV hierarchy is given by
\begin{equation}
q_{t} + (L^A)^n q_{x} = 0, ~~n=0,1,2,3, \ldots.
\end{equation}When n = 1, this operator formulation yields the KdV equation
\begin{equation}
q_{t} + 6qq_{x}+q_{xxx}= 0.
\end{equation}When n = 2, we find Lax 5th order equation
\begin{equation}
q_t+q_{xxxxx} + 10qq_{xxx}+20q_xq_{xx}+30q^2q_x =0
\end{equation} and so on. In Ablowitz et al. [Reference Ablowitz, Kaup, Newell and Segur7], the most general case considered was when 
$\gamma(L^A)$ was taken to be a meromorphic function.
Suppose we wish to analyse fractional non-linear equations such as a fractional KdV equation where
\begin{equation*}
\gamma(L^A) = - 4 L^A \left | - 4 L^A \right |^{\alpha} \text{~with~} -1 \lt \alpha \lt 1.
\end{equation*} However, we are confronted with the question: what is the meaning of a fractional power of the operator LA; i.e., how to express 
$\left| L^A \right|^{\alpha}$ in physical space? In Ablowitz et al. [Reference Ablowitz, Been and Carr5], it was shown that the integrable fractional KdV equation
\begin{equation}
q_{t} - 4 L^A \left| 4 L^A \right|^{\alpha} q_{x} = 0
\end{equation}and integrable fractional KdV equation hierarchy
\begin{equation}
q_{t} + (-4L^A)^n \left| 4 L^A \right|^{\alpha} q_{x} = 0, ~~n=1,2, \ldots
\end{equation}(n = 1 is the integrable fractional KdV equation) can be written in terms of concrete functions associated with the time-independent Schrödinger equation (1.6). For α = 0, equation (1.12) provides an evolution equation with integer derivatives. For fractional α, an explicit non-local evolution equation can be written in term of the eigenfunctions of Schödinger equation (1.6). The result is (see the appendix for more information)
\begin{equation}
q_{t} +\! \int_{\Gamma_{\infty}}\!\! dk (-4k^2)^n |4 k^2|^{\alpha} \frac{\tau^2(k)}{4 \pi i k}\int_{-\infty}^{\infty}\! dy G(x,y,k) q_{y} = 0, ~~n=1,2, \ldots
\end{equation}where
\begin{equation*}
G(x,y,k) = \partial_{x} (\psi^2(x,k)) \varphi^2(y,k) - \partial_{x}(\varphi^2(x,k)) \psi^2(y,k),
\end{equation*} and 
$\psi(x,k), \varphi(x,k)$ satisfy (1.6) with 
$~~~\varphi(x,k) \sim e^{-ikx}, x \to -\infty, ~~~\psi(x,k ) \sim e^{ikx}$, 
$x \to \infty$ and where
\begin{equation*}
\varphi(x,k) \tau(k) = \psi(x,-k) + \rho(k) \psi(x,k)
\end{equation*} with 
$\tau(k), \rho(k)$ being the ‘transmission, reflection’ coefficients, respectively; it can be shown that 
$\psi(x,k)e^{-ikx}, \varphi(x,k)e^{ikx}$ are analytic in the upper half plane (UHP), 
$\tau(k)$ is meromorphic in the UHP; 
$\Gamma_{\infty}=\lim_{R \to \infty}\Gamma_R$ is the semicircular contour in the UHP from 
$k = - R$ to k = R, above all poles of 
$\tau(k)$ which are simple and finite in number and k = 0 is a removable singular point. Equation (1.13) provides explicit meaning of the integrable fractional KdV hierarchy (1.12). Finding the functions in equation (1.12) requires finding the solution to the direct side of the time-independent Schrödinger equation (1.6).
At first glance, equation (1.13) may appear unusual, but in fact linear equations solvable by Fourier transforms can be put in a similar form. In this regard, consider the linear dispersive equation
\begin{equation}
q_t+\gamma(\partial_x^2)q_x=0.
\end{equation}Using the Fourier representation of qx, we can rewrite equation (1.14) as
\begin{equation}
q_t+\gamma(\partial_x^2) \frac{1}{2\pi}\int_{-\infty}^{\infty}dk e^{ikx}\int_{-\infty}^{\infty}dy e^{-iky}q_y=0
\end{equation}or
\begin{equation}
q_t+\frac{1}{2\pi}\int_{-\infty}^{\infty}dk \gamma(-k^2) \int_{-\infty}^{\infty}dy e^{ik(x-y)}q_y=0.
\end{equation}A typical example is the fractional linear KdV hierarchy where
\begin{equation*}
\gamma(\partial_x^2)= (\partial_x^2)^n|\partial_x|^{\alpha} = \gt \gamma(-k^2)= (-k^2)^n|k^2|^{\alpha}, n=1,2, \ldots, -1 \lt \alpha \lt 1.
\end{equation*} In the linear limit equation (1.13) reduces to equation (1.16) with these choices of 
$ \gamma(k)$. We remark that linear fractional equations are intrinsically non-local.
1.3. IST solution of the fractional KdV hierarchy
In order to linearize/solve the integrable factional KdV hierarchy, we employ inverse scattering. Briefly, this can be stated as follows. Solve the Gel’fand, Levitan, Marchenko (GLM) equation:
\begin{align}
K(x,y;t) &+ F(x+y;t) + \int_{x}^{\infty} \! ds ~K(x,s;t) F(s+y;t) = 0
\end{align}
\begin{align}
F(x;t) &=F_c+F_d=~\frac{1}{2\pi} \int_{-\infty}^{\infty} dk \rho(k,t)e^{i k x}+ \sum_{j = 1}^{J} c_{j}(t) e^{- \kappa_{j} x}
\end{align} for 
$K(x,y,t)$, 
$y \geq x$; the solution of the fKdV equation is obtained from:
\begin{equation*}
q(x,t) = 2 \frac{d}{dx} K(x,x;t).
\end{equation*} In the GLM equation (1.17a–1.17b), the function 
$F(x;t)$ depends on the so-called ‘scattering’ data (including their time dependence) which includes both continuous and discrete data; continuous data: 
$\rho(k,t)= \frac{b(k,t)}{a(k,t)}, ~~~~~~~~~~ a(k,t) = a(k,0), \quad ~~~~~~~~~~b(k,t) =b(k,0) e^{- 2 i k \gamma(k^2) t}$ where 
$a(k,0),b(k,0)$ are known in terms of Wronskians of 
$\phi, \varphi$ and
\begin{equation*}
\gamma(k^2) = (-4 k^2)^n \left| 4 k^2 \right|^{\alpha}; -1 \lt \alpha \lt 1, n=1,2, \ldots.
\end{equation*} Discrete data: 
$k_j=i\kappa_j$, discrete eigenvalues; 
$a(k_j) = 1/\tau(k_j)=0, ~j=1,2, \ldots J$: constant in time, and 
$ c_j(t)=i b_j/a'(k_j)$, where at an eigenvalue 
$k_j, \phi(x,k_j)= b_j\psi(x,k_j)$; 
$c_j(t)$ are often referred to as norming constants.
We see that the only difference in the solution between KdV and fractional KdV is in the time dependence: 
$ \gamma(k^2)$. This is similar to the linear case. ‘Pure’ soliton solutions are obtained when 
$\rho(k,0)=0$ in which case 
$F=F_d$. Further, as 
$t \to \infty, ~ F_c \to 0 = \gt F \to F_d$ which in the long time limit leads to soliton solutions. A one soliton solution to fKdV is obtained when 
$\rho(k,0)=0, k_1=i\kappa_1, \kappa_1 \gt 0, c_1=2\kappa_1e^{2\kappa_1x_1}$; the one soliton solution is given by
\begin{equation*}
q_{sol}(x,t) = 2 \kappa_1^{2} \text{sech}^{2}{\{\kappa_1 ((x - x_{1}) -(4 \kappa_1^2)^{1+\alpha} t)\}}.
\end{equation*}2. Integrable fractional Burgers equation
The Burgers hierarchy is given by
\begin{equation}
u_t+\sigma \partial_xL^nu=0 \text{~where~} L= (\partial_x+u(x,t)), ~\sigma \text{~const.}
\end{equation}The Burgers hierarchy and certain associated solutions were discussed in Kudryashov and Sinelshchikov [Reference Kudryashov and Sinelshchikov13].
For n = 0, (2.1) reduces to the transport equation,
\begin{equation}
u_t+\sigma u_x=0.
\end{equation}For n = 1, we find the celebrated Burgers equation,
\begin{equation}
u_t+\sigma(u_{xx}+2uu_x)=0.
\end{equation}And for n = 2, we have
\begin{equation}
u_t+\sigma\big(u_{xxx}+3u_x^2+3uu_{xx}+3u^2u_x\big)=0,
\end{equation}which was found by Tasso [Reference Tasso29] (see also Sharma and Tasso [Reference Sharma and Tasso27]) and Olver [Reference Olver21], who used the property that the equation has an infinite number of symmetries. Equation (2.4) is sometimes referred to as the Sharma-Tasso-Olver (STO) equation.
In Kudryashov and Sinelshchikov [Reference Kudryashov and Sinelshchikov13], the solution of the Burgers hierarchy is obtained by using the Cole–Hopf transformation
\begin{equation}
u= \frac{\psi_x}{\psi}
\end{equation}and it is shown that
\begin{equation}
u_t+\sigma \partial_xL^{n}u= \partial_x\left(\frac{\psi_t+\sigma \psi_{n+1,x}}{\psi}\right)=0, ~n=0,1,2, \ldots
\end{equation} where 
$\psi_{n+1,x}$ represents n + 1 derivatives in x.
In order to obtain equation (2.6), we use the following lemma [Reference Kudryashov and Sinelshchikov13].
Lemma 2.1. Under the transformation (2.5), the operator Ln from (2.1) becomes
\begin{align*}
L^n\frac{\psi_x}{\psi} &= \Big(\partial_x+\frac{\psi_x}{\psi}\Big)^n\frac{\psi_x}{\psi} \\
&=\frac{\psi_{n+1,x}}{\psi}.
\end{align*}Proof: The equation in Lemma (2.1) is obtained by induction. When n = 0, we have
\begin{equation*}
L^0\frac{\psi_x}{\psi} = \frac{\psi_x}{\psi}.
\end{equation*}If we assume the lemma holds for n = k, i.e.,
\begin{equation*}
\Big(\partial_x+\frac{\psi_x}{\psi}\Big)^k\frac{\psi_x}{\psi}= \frac{\psi_{k+1,x}}{\psi},
\end{equation*}then
\begin{align*}
\Big(\partial_x+\frac{\psi_x}{\psi}\Big)^{k+1}\frac{\psi_x}{\psi} &=
\Big(\partial_x+\frac{\psi_x}{\psi}\Big)\Big(\partial_x+\frac{\psi_x}{\psi}\Big)^{k} \frac{\psi_x}{\psi}\\
&=\Big(\partial_x+\frac{\psi_x}{\psi}\Big)\frac{\psi_{k+1,x}}{\psi}\\
&=\frac{\psi_{k+2,x}}{\psi}-\frac{\psi_{k+1,x}}{\psi^2}\psi_x +\frac{\psi_{k+1,x}}{\psi^2}\psi_x\\
&=\frac{\psi_{k+2,x}}{\psi}.
\end{align*}Thus, by induction, the lemma holds for all n.
We can now derive equation (2.6) using the Cole–Hopf transformation (2.5) and Lemma (2.1). This follows from
\begin{align*}
u_t+\sigma \partial_xL^{n}u & = \frac{\psi_{xt}}{\psi} -\frac{\psi_x\psi_t}{\psi^2} +\sigma\partial_x\left( \frac{\psi_{n+1,x}}{\psi} \right)\\
& = \partial_x\left( \frac{\psi_t}{\psi}\right)+\sigma\partial_x\left( \frac{\psi_{n+1,x}}{\psi}\right) \\
& =\partial_x\left( \frac{\psi_t+\sigma \psi_{n+1,x}}{\psi} \right)\\
&=0.
\end{align*}The issue we address here is to understand the fractional Burgers equation and fractional Burgers hierarchy
\begin{equation}
u_t+\sigma \partial_xL^{n}|L|^{\alpha}u=0 \text{~ with~} -1 \lt \alpha \lt 1.
\end{equation}Importantly, when we operate on a pure exponential or evaluate a Fourier transform we use the following formula for fractional derivatives
\begin{equation}
\partial^n_x|\partial_x|^{\alpha} e^{\mu x}= \mu^n|\mu|^{\alpha} e^{\mu x}.
\end{equation}2.1. Inverse approach
In analogy with the fractional KdV equation/hierarchy, we will consider both the ‘direct and inverse’ side of Burgers equation/hierarchy. Analogous to the fact that the time-independent Schrödinger equation is key to both the direct and inverse side of the KdV hierarchy, the Cole–Hopf equation is central to both the direct and inverse side of the Burgers hierarchy. First, we consider the inverse, which has a natural extension in terms of the linearized equation. If we formally replace n by 
$n+|\alpha|, -1 \lt \alpha \lt 1$, then equation (2.6) takes the form
\begin{equation}
\partial_x\left(\frac{\psi_t+\sigma \psi_{n+|\alpha|+1,x}}{\psi}\right) = 0.
\end{equation}We note that equation (2.9) can be evaluated in a Fourier context and can be considered as the fractional continuation of equation (2.6).
After integrating once, equation (2.9) becomes
\begin{equation}
\psi_t+\sigma \psi_{n+|\alpha|+1,x}=C(t)\psi,
\end{equation}where C(t) is an arbitrary function of time. We note that C(t) can be absorbed into ψ by rescaling. We can rewrite equation (2.10) in the form
\begin{equation}
\psi_t+\sigma \partial_x^{n+1}|\partial_x|^{\alpha} \psi=C(t) \psi.
\end{equation} The term 
$|\partial_x|^{\alpha}$ yields well-posed solutions. This is the 1D analog of the fractional Laplacian. Under the Cole–Hopf transformation (2.5), equation (2.11) can be viewed as the solution of the inverse problem associated with the fractional Burgers hierarchy. Similarly, the analog of the direct problem for the Burgers hierarchy takes the form
\begin{equation}
u_t+\sigma \partial_xL^{n}|L|^{\alpha}u=0.
\end{equation}We will refer to the case n = 1 as the fractional Burgers equation.
Taking 
$C(t) = 0$, we can use the definition (2.8) to write solutions to equation (2.10) in terms of Fourier transforms or series. For n = 0 and σ real, equation (2.10) has oscillating exponential solutions of the form
\begin{equation*}
e^{i(k x - \sigma k |k|^{\alpha} t)}
\end{equation*}and a Fourier transform solution of the form
\begin{equation}
\psi(x,t) - 1 = \frac{1}{2\pi}\int_{-\infty}^{\infty} \widehat{\psi}(k,0) e^{i(kx-\sigma k |k|^{\alpha} t)} \, dk
\end{equation}where
\begin{equation*}
\widehat{\psi}(k,0) =\displaystyle\int_{-\infty}^{\infty} \big(\psi(x,0)-1\big)e^{-i k x} \,dx
\end{equation*} is the Fourier transform of the initial value. The initial value 
$\psi(x,t=0)$ is obtained from the Cole–Hopf transformation (2.5) as
\begin{equation}
\psi(x,0)=e^{\int_{-\infty}^x u_0(x')dx'}
\end{equation} where we normalize ψ to be unity as 
$x \to -\infty$. So, assuming 
$u(x,t)$ is decaying and initially 
$\int_{-\infty}^{\infty}u_0(x)dx=0$, there exists a Fourier transform for 
$\psi-1$.
In Figure 1, we see example evolutions for the fractional extension of the transport equation that have been calculated using the solution found from the inverse method (2.13) and computationally obtained using discrete Fourier transforms. Here, we see that the fractional term is associated with a dispersive front that resembles a dispersive shock wave.
Evolution of the fractional Burgers extension of the transport equation (n = 0) with σ = 1 and various values of α. The initial condition is 
$u(x,0) = \big(x-x_0\big)e^{-(x-x_0)^2}$ with 
$x_0=5$. (a) Contour plots. (b) Profiles taken at t = 5.
For n = 1 and 
$\sigma = -1$, we have the fractional Burgers equation with associated plane wave solutions of the form
\begin{equation*}
\psi_k(x,t) =e^{ikx-k^2|k|^{\alpha} t}
\end{equation*}which decay rapidly in time and general solution
\begin{equation}
\psi(x,t)=1 + \frac{1}{2\pi}\int_{-\infty}^{\infty} \hat{\psi}(k,0)e^{ikx-k^2|k|^{\alpha} t}dk.
\end{equation}We can also analyse periodic solutions via
\begin{equation*}
\psi(x,t)=1+ \sum_{n= -\infty}^{\infty} \hat{\psi}_n(t)e^{inx-n^2|n|^{\alpha} t},
\end{equation*}which is being used for the numerical examples.
In Figure 2, we see example evolutions for the fractional Burgers equation where solutions have been calculated using the solution found from the inverse method (2.15) and computationally obtained using discrete Fourier transforms. Unlike the n = 0 case, we see less pronounced effects from increasing the fractional constant α with a ‘lip’ developing on the shock-like side of the wave fronts. Note that there is not a smooth transition from the n = 0 case to the n = 1 case; we do not expect such a transition due to the absolute value in the Fourier solution of (2.11) which has waves of the form 
$e^{ikx-(ik)^{(n+1)}|k|^{\alpha} t}$, which maintains dispersion when n is even and diffusion when n is odd.
Evolution of the fractional Burgers equation (n = 1) with 
$\sigma = -1$ and various values of α. The initial condition is 
$u(x,0) = e^{-(x+x_0)^2/\omega} -e^{-(x-x_0)^2/\omega}$ with 
$x_0=10$ and ω = 4. (a) Contour plots. (b) Profiles taken at t = 5.
For n = 2 and σ = 1, similarly to the n = 0 case, our equation admits a family of periodic exponential solutions. The general solution takes the form
\begin{equation}
\psi(x,t) = 1 + \frac{1}{2\pi}\int_{-\infty}^{\infty} \hat{\psi}(k,0)e^{i\big(kx - \sigma k^3|k|^{\alpha} t\big)}dk.
\end{equation}In Figure 3, we see example evolutions for the n = 2 case found from the inverse method (2.16) and obtained computationally using discrete Fourier transforms. Note: the fractional equations also have solutions that exhibit blow up. This occurs when the non-linear effects are stronger and there are zeroes in the transformed equation, i.e., ψ = 0.
Evolution of the fractional Olver equation (n = 2) with σ = 1 and various values of α. The initial conditions is 
$u(x,0) = (x-x_0) e^{-(x-x_0)^2}$ with 
$x_0=80$. (a) Contour plots. (b) Profiles taken at t = 1.
Thus, using the concept of direct and inverse equations, we have seen that on the inverse side the fractional equation (2.11) can be evaluated via Fourier methods.
2.2. Direct approach
A key question remaining is to understand the meaning and method of calculation of the fractional equation (2.12) on the direct side. Here the part that needs to be understood is how to calculate 
$|L|^{\alpha}$ where 
$L=(u(x,t)+\partial_x)$. For this, we once again employ the Cole–Hopf transformation (2.5). Here the Cole–Hopf transformation is the analog of the time-independent Schrödinger equation that is associated with the KdV/fractional KdV equation.
Consider the fractional continuation of the Lemma (2.1)
\begin{align}
L^n|L|^{\alpha} \frac{\psi_x}{\psi}&=\Big(\partial_x+\frac{\psi_x}{\psi}\Big)^n|\partial_x+\psi_x|^{\alpha} \frac{\psi_x}{\psi}\notag \\ &=
\frac{\partial_x^{n+1}|\partial_x|^{\alpha}\psi}{\psi}
\end{align}which we calculate in Fourier space using the derivative formula (2.8). Then we have
\begin{align}
0&=\partial_t u+\sigma \partial_xL^n|L|^{\alpha}u
\end{align}
\begin{align}
&= \partial_t u+\sigma \partial_x\left(\frac{\partial_x^{n+1}|\partial_x|^{\alpha}\psi}{\psi}\right)
\end{align}
\begin{align}
&=\partial_tu + \sigma \partial_x \frac{\partial_x^{n}|\partial_x|^{\alpha}\psi_x}{\psi},
\end{align} where the initial condition 
$u(x,t=0)=u_0(x)$ is given on 
$|x| \lt \infty$. We calculate 
$\psi, \psi_x$ from u using the Cole–Hopf transformation (2.5) via
\begin{align}
\psi(x,t) &= e^{\int_{-\infty}^x u(x',t)dx'}
\end{align}
\begin{align}
\psi_x(x,t) &= u(x,t)e^{\int_{-\infty}^x u(x',t)dx'}
\end{align} where we normalize ψ to be unity as 
$x \to -\infty$ and assume 
$u(x,t)$ is decaying and initially 
$\int_{-\infty}^{\infty}u_0(x)dx=0$. Now, ψx has a Fourier transform and we can evaluate the fractional Burgers hierarchy using equation (2.18c).
We calculate the direct method numerically using a standard 4th order Runge–Kutta solver in time and spectral methods to compute the spatial derivatives on the discretized system obtained by taking the Fourier transform of equation (2.18c). This leads us to
\begin{align}
\mathrm{\mathbf{v}}_t & = - i \sigma \mathrm{\mathbf{k}}\,\mathcal{F}\left(\frac{i \mathrm{\mathbf{k}}^{n} |\mathrm{\mathbf{k}}|^ \alpha \mathrm{\mathbf{\psi}}_x}{\mathrm{\mathbf{\psi}}}\right)
\end{align}
\begin{align}
\mathrm{\mathbf{v}}(0) & = \mathcal{F}\big[u(\mathrm{\mathbf{x}},0) \big]
\end{align}
\begin{align}
\mathrm{with} \quad \mathrm{\mathbf{\psi}} &= \mathrm{Exp}\Bigg(\int \mathcal{F}^{-1}\big[\mathrm{\mathbf{v}}\big]\Bigg) \quad \mathrm{and} \quad \mathrm{\mathbf{\psi}}_x = \mathcal{F}^{-1}\big[\mathrm{\mathbf{v}}\big]\mathrm{\mathbf{\psi}}
\end{align} where 
$\mathrm{\mathbf{x}}$ and 
$\mathrm{\mathbf{k}}$ are compatible discretizations of the real and Fourier space respectively, 
$\mathcal{F}$ and 
$\mathcal{F}^{-1}$ are the discrete Fourier transform and inverse discrete Fourier transform and 
$\int$ is an approximation of 
$\int_{-\infty}^{x} u(x',t) dx'$ using a cumulative trapezoid rule. This scheme converges to solutions obtained from the inverse method as the number of Fourier modes taken increases; this is illustrated in Figure 4. Both the inverse and direct methods are in full agreement, though we note that the direct method requires small time stepping to converge.
(Left) Comparison of the numerical solutions to the fractional Burgers equations found using the inverse method (solid blue) and the direct method (dashed red) with α = 0.9, 
$N=2^{10}$ (Fourier modes) and time step 
$\Delta t = 6.25 \times 10^{-5}$. (Right) Difference between the inverse method and the direct method with an increasing number of Fourier modes taken.
Without using the Cole–Hopf formula (2.5) and Fourier fractional derivatives of exponentials given by equation (2.8), it is not clear how to evaluate the fractional extension of the L operator 
$L^n|L|^{\alpha}$. However, for special cases we can use the binomial expansion of Ln= 
$(\partial_x+u)^n=u^n(1+\frac{\partial_x}{u})^n$ when the operator 
$\frac{\partial_x}{u}$ is smaller than unity. In certain cases, analytic continuation can lead to a result valid when 
$\frac{\partial_x}{u}$ is not small. We also note that we can use other definitions of fractional derivatives on the forward and inverse side. However, the advantage of Fourier derivatives is that it leads to explicit results.
2.3. Shock solutions
Next we turn to studying special traveling wave (TW) solutions of these equations. We begin with Burgers equation (2.3) with 
$\sigma=-1$. A TW solution is obtained from equation (2.11) with 
$n=1, \alpha=0, C(t)=C_0$ constant,
\begin{equation}
\psi_t- \partial_x^{2}\psi=C_0 \psi.
\end{equation}A TW solution is obtained by looking for solutions of the form
\begin{equation}
\psi(x,t)=\psi(\zeta) \text{~~where~} \zeta= x-Vt, ~~ V \text{~constant}.
\end{equation}This yields the ODE
\begin{equation}
\psi_{\zeta \zeta}+V\psi_{\zeta}+C_0\psi=0.
\end{equation} Looking for solutions of the form 
$e^{r \zeta}, r$ constant, yields
\begin{equation*}
r^2+Vr+C_0=0 ~= \gt ~r^2+Vr+C_0=0,
\end{equation*}which has two solutions
\begin{equation*}
r_{\pm}=-\frac{V}{2}\pm\sqrt{\left(\frac{V}{2}\right)^2-C_0}.
\end{equation*} When 
$\frac{V}{2}^2-C_0 \gt 0$, there are two real solutions; we write the solution ψ as
\begin{equation}
\psi=c_1e^{r_+\zeta}+ c_2e^{r_-\zeta}
\end{equation}and the solution of Burgers equation as
\begin{align}
u(\zeta)&= \frac{\psi_{\zeta}}{\psi}
\end{align}
\begin{align}
&= \frac{c_1 r_+ e^{r_+\zeta}+c_2r_- e^{r_-\zeta}}{c_1e^{r_+\zeta}+ c_2e^{r_-\zeta}}
\end{align}
\begin{align}
&=\frac{r_{+} +r_- e^{(r_{-}-r_+)(\zeta-\zeta_0)}}{1+e^{(r_{-}-r_+)(\zeta-\zeta_0)}},
\end{align} where 
$\frac{c_2}{c_1}= e^{-(r_{-}-r_{+}) \zeta_0}$. This shock solution increases from 
$r_-$ as 
$\zeta \to -\infty$ to 
$r_+$ as 
$\zeta \to \infty$; we also note that 
$r_-+r_+=-V$. When 
$(\frac{V}{2})^2-C_0=0$ then r has a double root; we find
\begin{equation*}
\psi=c_1e^{-\frac{V}{2}\zeta}+ c_2\zeta e^{-\frac{V}{2}\zeta},
\end{equation*}which corresponds to the following rational solution of Burgers equation
\begin{equation*} u(\zeta)= \frac{\psi_{\zeta}}{\psi}= -\frac{V}{2}+\frac{1}{\zeta-\zeta_0},\end{equation*} where 
$\frac{c_1}{c_2}=-\zeta_0$. This solution is singular when 
$\zeta=\zeta_0$. When 
$\frac{V}{2}^2+C_0 \lt 0$, there are singular oscillatory solutions; we will not go into further detail on the these solutions here.
Now we will discuss the fractional Burgers equation (2.12) with 
$n=1, -1 \lt \alpha \lt 1, \sigma=-1$ whose solution is obtained from equation (2.11). As we did for Burgers equation, we look for TW solutions 
$\psi(x,t)=\psi(\zeta)$ with 
$C(t)=C_0$ constant. The fractional TW equation is given by
\begin{equation}
\partial_{\zeta}^2|\partial_{\zeta}|^{\alpha} \psi+V\partial_{\zeta} \psi +C_0\psi=0.
\end{equation} Solutions of the form 
$\psi=e^{r \zeta}, r$ constant, lead to an equation for r of the form
\begin{equation}
r^2|r|^{\alpha}+Vr+C_0=0.
\end{equation} Graphical analysis indicates that there can be two solutions 
$r_{\pm}$ where 
$r_+ \gt r_-$. In this case, the equation for ψ has the same form as in equation (2.24)
\begin{equation}
\psi=c_1e^{r_+\zeta}+ c_2e^{r_-\zeta}.
\end{equation}Hence,
\begin{align}
u(\zeta)&= \frac{\psi_{\zeta}}{\psi}
\end{align}
\begin{align}
&\qquad\qquad\qquad\qquad\quad\!\!=\frac{c_1 r_+ e^{r_+\zeta}+c_2r_- e^{r_-\zeta}}{c_1e^{r_+\zeta}+ c_2e^{r_-\zeta}}
\end{align}
\begin{align}
&\qquad\qquad\qquad\qquad\quad=\frac{r_{+} +r_- e^{(r_{-}-r_+)(\zeta-\zeta_0)}}{1+e^{(r_{-}-r_+)(\zeta-\zeta_0)}}
\end{align} where 
$\frac{c_2}{c_1}= e^{-(r_{-}-r_{+}) \zeta_0}$. These TW shock waves are depicted in Figure 5. For 
$C_0 = 0$, we have 
$r_+ = 0$ and 
$r_-=- V ^{1/(1+\alpha)}, V \gt 0$.
Finally, let’s consider the next equation in the hierarchy (2.12) with 
$\sigma=1, n=2$. The function ψ we will consider satisfies
\begin{equation}
\partial_t\psi + \partial_x^3|\partial_x|^{\alpha}\psi=C(t)\psi.
\end{equation} TW solutions, 
$\psi(x,t)=\psi(\zeta)$ with 
$C(t)=C_0$ constant, satisfy
\begin{equation}
\partial_{\zeta}^3|\partial_{\zeta}|^{\alpha} \psi-V\partial_{\zeta} \psi -C_0\psi=0.
\end{equation} Looking for solutions of the form 
$\psi=e^{r \zeta}, r$ constant, leads to the following equation for r
\begin{equation}
r^3|r|^{\alpha}-Vr-C_0=0.
\end{equation} For 
$-1 \lt \alpha \lt 1$ depending on the sign and size of C 0 and V, graphical analysis indicates there can be three real solutions 
$r_3 \gt r_2 \gt r_1$. Then we have
\begin{equation}
\psi=c_1e^{r_1\zeta}+ c_2e^{r_2\zeta} +c_3e^{r_3\zeta}
\end{equation}and
\begin{align}
u(\zeta) &= \frac{\psi_{\zeta}}{\psi}
\end{align}
\begin{align}
&\qquad\qquad\qquad\qquad\qquad\qquad\quad\!\!=\frac{c_1 r_1 e^{r_1\zeta}+c_2r_2e^{r_2\zeta}+c_3r_3e^{r_3\zeta} }{c_1e^{r_1\zeta}+ c_2e^{r_2\zeta} +c_3e^{r_3\zeta}}
\end{align}
\begin{align}
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\qquad\!\!=\frac{r_1 e^{(r_1-r_2)(\zeta-\zeta_-)} + r_2 +r_3e^{(r_3-r_2)(\zeta-\zeta_+)} }{e^{(r_1-r_2)(\zeta-\zeta_-)}+1+e^{(r_3-r_2)(\zeta-
\zeta_+)}}
\end{align} where 
$\frac{c_1}{c_2}= e^{-(r_{1}-r_{2})\zeta_-}, \frac{c_3}{c_2}= e^{-(r_{3}-r_{2})\zeta_+}$. So, as 
$\zeta \to \infty, u \sim r_3$ and as 
$\zeta \to -\infty$, 
$u \sim r_1$; hence we have a shock-like TW solution with some additional interior structure. Typical shock waves with a middle plateau are depicted in Figure 6. For 
$C_0 = 0$, we have 
$r_2 = 0$, 
$r_1= - V ^{1/(2+\alpha)}$ and 
$r_3= V ^{1/(2+\alpha)}, V \gt 0$.
3. Conclusion
Fractional integrable Burgers equation and hierarchy are considered. Fractional derivatives are calculated via Fourier methods. Using the Cole–Hopf transformation, we formulate these Burgers equations on the direct side and discuss how to find these solutions. On the inverse side, we use the Cole–Hopf transformation to find a fractional partial differential equation. Results from fractional integrable KdV equation and hierarchy are used as a guide.
Funding Statement
This project was partially supported by NSF under grant number DMS-2306290.
Appendix
To establish equation (1.13), we use a result involving completeness of suitable squared eigenfunctions of the time-independent Schrödinger equation. Concretely, in Sachs [34], it was shown that a rapidly decreasing function h(x) can be written in the form
\begin{align}
h(x) = \int_{\Gamma_{\infty}} dk ~\frac{\tau^2(k)}{4 \pi i k} \int_{-\infty}^{\infty} dy ~G(x,y,k) h(y)
\end{align}(note that the variable t is suppressed); in the above:
\begin{equation*}G(x,y,k) = \partial_{x} (\psi^2(x,k)) \varphi^2(y,k) - \partial_{x}(\varphi^2(x,k)) \psi^2(y,k),\end{equation*} where 
$\psi(x,k), \varphi(x,k)$ satisfy (1.6) with
\begin{equation*} ~~~\varphi(x,k) \sim e^{-ikx}, x \to -\infty, ~~~\psi(x,k ) \sim e^{ikx}, x \to \infty\end{equation*}and
\begin{equation*} \varphi(x,k) \tau(k) = \psi(x,-k) + \rho(k) \psi(x,k)\end{equation*} with 
$\tau(k), \rho(k)$ the ‘transmission, reflection’ coefficients, respectively; it can be shown that 
$\psi(x,k)e^{-ikx}, \varphi(x,k)e^{ikx}$ are analytic in the UHP and 
$\tau(k)$ is meromorphic in the UHP; 
$\Gamma_{\infty}=\lim_{R \to \infty}\Gamma_R$ is the semicircular contour in the UHP from 
$k = - R$ to k = R, above all poles of 
$\tau(k)$ which are simple and finite (J) in number; k = 0 is a removable singular point. Due to analyticity we can deform the integral along 
$\Gamma_{\infty}$:
\begin{equation*}\int_{\Gamma_{\infty}} dk= \int_{-\infty}^{\infty}dk + \sum_{j=1}^{j=J}(2 \pi i) (\text{residues}).\end{equation*} It is worth noting that in the linear limit: 
$\int_{\Gamma_{\infty}} \sim \int_{-\infty}^{\infty}; ~~\varphi(x,k) \sim e^{-ikx},\\
\psi(x,k ) \sim e^{ikx}, \tau \to 1$ and there are no poles in 
$\tau(k)$. Then the completeness relation reduces to
\begin{equation*}h(x)= \frac{1}{2 \pi} \int_{-\infty}^{\infty} dk e^{ikx} \int_{-\infty}^{\infty}dy e^{-iky} h(y),\end{equation*}which is the completeness relation for Fourier transforms.
Recall the operator
\begin{equation}
L^{A} =-\frac{1}{4} \partial_{x}^2 - q + \frac{1}{2} q_{x} \int_{x}^{\infty} dy.
\end{equation}Its adjoint is
\begin{equation}
L =-\frac{1}{4} \partial_{x}^2 - q + \frac{1}{2} \int_{-\infty}^{x} q_{y}dy.
\end{equation}In Ablowitz et al. [29], the following spectral relations are derived for L A
\begin{equation}
L^{A} V^A=k^2 V^A, ~~ \text{with~~} V^A= \partial_{x} (\psi^2(x,k)) \text{~or~} \partial_{x} (\varphi^2(x,k))
\end{equation}and L:
\begin{equation}
L V=k^2 V, ~~ \text{with~~} V= \psi^2(x,k) \text{~or~} \varphi^2(x,k).
\end{equation}Recall 
$G(x,y,k) = \partial_{x} (\psi^2(x,k)) \varphi^2(y,k) - \partial_{x}(\varphi^2(x,k)) \psi^2(y,k)$. Therefore, we have 
$L^{A}G(x,y,k)=k^2G(x,y,k)$ and similarly 
$|L^{A}|^{\alpha}G(x,y,k)=|k^2|^{\alpha} G(x,y,k)$. So with 
$\gamma(L^A) = - 4 L^A \left| 4 L^A \right|^{\alpha},~~-1 \lt \alpha \lt 1,$ using completeness the factional KdV equation given by
\begin{equation*}q_{t} - 4 L^A \left| 4 L^A \right|^{\alpha}) q_{x} = 0\end{equation*}is given by (t suppressed)
\begin{align}
q_{t} -\! \int_{\Gamma_{\infty}}\!\! dk 4k^2 |4 k^2|^{\alpha} \frac{\tau^2(k)}{4 \pi i k}\int_{-\infty}^{\infty}\! dy G(x,y,k) q_{y} = 0.
\end{align}Using the properties of squared eigenfunctions, we can write the above equation as
\begin{align}
q_{t} +\! \int_{\Gamma_{\infty}}\!\! dk |4 k^2|^{\alpha} \frac{\tau^2(k)}{4 \pi i k}\int_{-\infty}^{\infty}\! dy G(x,y,k) \left( 6 q q_{y} + q_{yyy} \right) = 0.
\end{align}Next it is shown how to go from equation (A.6) to (A.7). We use
\begin{equation*} -4k^2G(x,y,k)= \partial_{x} (\psi^2(x,k))(-4k^2) \varphi^2(y,k) - \partial_{x}(\varphi^2(x,k)) (-4k^2) \psi^2(y,k). \end{equation*}Then from
\begin{equation*} k^2(\varphi^2, \psi^2)^T(y,k)=L(\varphi^2, \psi^2)^T(y,k) \text{~where~} L=-\frac{1}{4}\partial_y^2-q(y)+\frac{1}{2}\int_{-\infty}^yq_{y'}(y')dy'\end{equation*}we have
\begin{equation*} -4k^2G(x,y,k)= \partial_{x} (\psi^2(x,k))(-4L) \varphi^2(y,k) - \partial_{x}(\varphi^2(x,k)) (-4L) \psi^2(y,k). \end{equation*}Using the operator L gives
\begin{align}
-4k^2G(x,y,k) q_y &= \partial_{x} (\psi^2(x,k))\left( ( \partial_y^2\varphi^2(y))q_y+4q\varphi^2(y)q_y-2q_{y} \int_{-\infty}^ydy'q_{y'}\varphi^2(y') \right)
\end{align}
\begin{align}
& - \partial_{x} (\varphi^2(x,k))\left( ( \partial_y^2\psi^2(y))q_y+4q\psi^2(y)q_y-2q_{y} \int_{-\infty}^ydy'q_{y'}\phi^2(y') \right).
\end{align} The 
$ \partial_y^2$ terms are integrated by parts twice to yield 
$G(x,y,k)q_{yyy}$ and then we interchange integrals in the last terms to find
\begin{equation}
-2 \int_{-\infty}^{\infty}dy q_{y} \int_{-\infty}^{y}dy'q_{y'}\varphi^2(y') = -2 \int_{-\infty}^{\infty} dy'q_{y'}\varphi^2(y')
\int_{y'}^{\infty}dyq_{y}= 2\int_{-\infty}^{\infty} dy'qq_{y}\varphi^2(y).
\end{equation}Combining terms yields
\begin{align*}
-4k^2G(x,y,k) q_y = G(x,y,k)(q_{yyy}+6qq_y),
\end{align*}which is the term on the right hand side of equation (A.7).
The above two equations (A.6) and (A.7) provide explicit meaning for the fractional KdV equation in terms of eigenfunctions and scattering data. We further note that when α = 0 we recover the KdV equation: 
$q_{t} + 6 q q_{x} + q_{xxx} = 0$.
Moreover in the linear limit: 
$\psi(x,k) \to e^{ikx},~~ \phi(x,k) \to e^{-ikx}, ~\tau \to 1$, then from equation (A.6) we have
\begin{equation*}q_{t} +\! \int_{-\infty}^{\infty}\!\! dk (-i) k^3 |k^2|^{\alpha} e^{ikx} \underbrace{\int_{-\infty}^{\infty} e^{-iky}q(y)dy}_{=\hat{q}(k)}=0.\end{equation*} Then using 
$ ~i k \to \partial_x$ we find the linear fractional KdV equation
\begin{equation*}q_t+ \partial_x^3|\partial_x^2|^{\alpha} q=0. \end{equation*}
Open access
