High order approximation scheme for a fractional order coupled system describing the dynamics of rotating two-component Bose-Einstein condensates

1.   Introduction

Of concern is proposing a numerical scheme based on a high-order finite difference/Legendre-Galerkin spectral method for solving the coupled Gross¢Pitaevskii equations in the dimensionless form with time and space fractional derivatives:

with the initial conditions

and the homogeneous boundary conditions

such that $\Omega = (a, b) \subset \mathbb {R}$ and $I = (0, T] \subset \mathbb {R}$. The parameters $\delta, \, \beta_{11},$ $\beta_{12},$ $\beta_{21},$ ${\beta_{22}}$ and $\lambda$ are constants to be given and $\phi_0(x)$ and ${\psi_0(x)}$ are given smooth functions.

The temporal fractional derivative is defined in Caputo sense ^[26], which means

The spatial fractional operator of Riesz type of order $\alpha$ with respect to $a \leq x \leq b$ ^[26], namely

where ${}_{a} D_{x}^{\alpha}\Psi(x, t)$ and ${}_x D_{b}^{\alpha}\Psi(x, t)$ are the left- and right-Riemann–Liouville derivatives of order $\alpha$ with respect to $x \in (a, b)$, and are defined as

and

Bose and Einstein predicted theoretically Bose-Einstein condensations (BEC) which were detected experimentally by Anderson et al. in 1995 ^[4]. The link between the spin angular momentum of the electron spin and the orbital angular momentum was performed by the spin-orbit (SO) coupling such as Rasha type and Dresselhaus type. The SO coupling had been extensively discussed experimentally. Bosons and fermions SO coupling was achieved by Jacob et al. in 2008 ^[20]. A non-dimensionalization and dimension reduction were applied in ^[12,40]. Then, a two-component BEC with an internal atomic Josephson junction (or an external driving field) can be well modeled by coupled Gross-Pitaevskii equations in dimensionless form.

Lately, Laskin extended the Feynman path integral approach over the more general Lvy-like quantum paths and derived a fractional Schrödinger equation, which modifies the integer equation by invoking the fractional Laplacian. It has been proposed to study BEC of which the particles obey a non-Gaussian distribution law ^[10,31,32], where fractional Schrödinger was named as fractional Gross¢Pitaevskii Equation (FGPE) and BEC as fractional BEC.

Due to the nonlocality of fractional differential operators, the numerical solutions of the fractional models are more complicated than the classical models. There are several analytical methods to solve fractional differential equations. However, analytical methods do not work well on most of the fractional differential equations, e.g. with nonlinearities or linear equations with time-dependent coefficients. From the numerical implementation point of view, the time-dependent Gross¢Pitaevskii equation describing the dynamics of rotating Bose-Einstein condensates and its discretization with the finite element method were considered in ^[19]. The approach in ^[14] exerted some efforts to propose a finite-difference method based on weighted-shifted Grünwald differences for solving the multi-dimensional Gross¢Pitaevskii equation, which considers fractional derivatives of the Riesz type in space, a generalized potential function and angular momentum rotation. An analysis based on a compact finite difference scheme was proposed in ^[33] for the integer-order coupled Gross¢Pitaevskii equations in one space dimension. That scheme can conserve the total mass and energy at the discrete level. In ^[21], a sine pseudo-spectral/difference scheme that preserves the discrete mass and energy was produced and analyzed for the integer-order coupled Gross¢Pitaevskii equations with Dirichlet boundary conditions in several spatial dimensions. The approach in ^[23] was devoted to analyzing the convergence of explicit finite difference schemes for computing the integer-order coupled Gross¢Pitaevskii equations in high space dimensions.

The combination of the efficiency of finite difference quotients based interpolation formulas of L1 or L2-$1_{\sigma}$ ^[3] and Galerkin Legendre spectral method is widely used to solve numerically different kinds of fractional order differential problems. For fractional diffusion problems, we refer to ^[13,36]. For the distributed-order weakly singular integral-partial differential model, we refer to ^[1]. For nonlinear fractional Schrödinger equations with Riesz space-and Caputo time-fractional derivatives, we refer to ^[35]. For a coupled system of time and space fractional diffusion equations, we refer to ^[17]. The propagation of solitons through a new type of quantum couplers called time-space fractional quantum couplers was presented in ^[18]. Concerning the Gross¢Pitaevskii equation arising in Bose-Einstein Condensation ^[25] as a generalization of the nonlinear fractional Schrödinger equations, numerous extensions to relevant physical situations are now clarified ^[6,5,8] (multi-components, nonlocal nonlinear interactions, etc.). For the fractional case, the situation is more complicated and still needs to be analyzed deeply. Serna-Reyes et al. ^[27] introduced and theoretically analyzed various numerical techniques for approximating the solutions of a fractional extension of a double condensate system that extends the well-known Gross¢Pitaevskii equation to the fractional scenario with two interacting condensates. Antoine et al. ^[7] proposed numerical schemes for time or space fractional nonlinear Schrödinger equations with some applications in Bose-Einstein condensation. Ainsworth and Mao ^[2] established the well-posedness of the fractional partial differential equation which arises by considering the gradient flow associated with a fractional Gross¢Pitaevskii free energy functional and some basic properties of the solution. Zhang et al. ^[39] studied the ground and first excited states of the fractional Bose¢Einstein condensates which are modeled by the fractional Gross¢Pitaevskii equation. They used the weighted shifted Grünwald¢Letnikov difference method to discretize the Gross¢Pitaevskii equation. Liang et al. ^[22] introduced efficient local extrapolation of the exponential operator splitting scheme to solve the multi-dimensional space-fractional nonlinear Schrödinger equations including the space-fractional Gross-Pitaevskii equation, which is used to model optical solitons in graded-index fibers.

In this paper, our goal is to numerically solve (1.1a)-(1.1b) by implementing a combined high-order numerical approach. This approach is based on the Alikhanov high-order interpolation scheme to be used to approximate the time Caupto fractional derivatives side by side to a Galerkin-type formulation base on Legendre orthogonal polynomials basis to approximate Riesz space fractional derivatives. We used the recently introduced discrete fractional Grönwall inequalities ^[24] in discrete energy estimates to prove the unconditional convergence of the proposed scheme.

2.   Numerical scheme

We fix the following notations.

● $(\cdot, \cdot)_{0, \Omega}$ denotes the inner product on the space $L^{2}(\Omega)$ with the $L^{2}$-norm $\| {\cdot} \|_{0, \Omega}$ and the maximum norm $\left\| \cdot \right\|_\infty$.

● $C_0^{\infty}(\Omega)$ denotes the space of non-singular functions with compact support in $\Omega$.

● $H^r(\Omega)$ and $H_0^r(\Omega)$ are Sobolev spaces with the norm $\| {\cdot} \|_{H^r}$ and semi-norm $\left|{\cdot}\right|_{H^r}$.

● $\mathbb{P}_N (\Omega)$ is the space of polynomials on $\Omega$ of degree less than or equal to $N$.

● The approximation space $V^0_N$ is defined as

● $I_N$ is the interpolation operator of Legendre-Gauss-Lobatto type, $I_N :C(\bar \Omega) \to V_N$,

We also define function spaces ^[11] which will be used in the construction of the numerical scheme.

Definition 1 (Fractional Sobolev space). The fractional Sobolev space $H^{\eta}(\Omega)$ for $\eta > 0$, is defined as

endowed with the semi-norm and norm respectively as

such that $H _0 ^\eta (\Omega)$ is the closure of ${C} _0 ^\infty (\Omega)$ with respect to $\| {\cdot} \| _{H ^\eta (\Omega)}$. Also, $\mathcal{F}(\tilde \Psi)$ is the Fourier transformation of the function $\tilde \Psi$ and the zero extension of $\Psi$ outside $\Omega$ denoted by $\tilde \Psi$.

Lemma 1 (Adjoint property). By choosing $1 < \eta < 2,$ then $\forall \, \Psi \in H_{0}^{\eta}(\Omega)$ and $\nu \in H_{0}^{\eta/2}(\Omega)$, we deduce

Discretization

The temporal domain $I$ is partitioned by $t_j = j \tau,$ $j = 0, \ 1, \ldots, M$ with $\tau = T/M$. Denote $t_{j+\sigma} = (j+\sigma)\tau = \sigma t_{j+1}+(1-\sigma)t_{j},$ for $j = 0, \ 1, \ldots, M-1$. We take $\Psi^{j+\sigma} = \Psi^{j+\sigma}(\cdot) = \Psi(\cdot, t_{j+\sigma})$.

Definition 2. Let $0 < \beta < 1$ and $\sigma = 1 - \frac {\beta} {2}$. Then we define

and

Lemma 2 (see ^[3]). L2-$1 _\sigma$ interpolation formula is given as follows. Assume that $\Psi(t) \in C^{3}[0, t_{j+1}]$, $0\leq j \leq M-1$, is formulated as

where $\delta_{t}\Psi^{r} = \Psi ^{r + 1} - \Psi^r.$

It can be also given as follows

where $d _{1} ^{(0, \beta, \sigma)} = -d _{0} ^{(0, \beta, \sigma)} = \sigma^{1-\beta}$ $\forall \, j = 0$, and $\forall \, j\geq 1,$

Accordingly, L2-$1 _\sigma$ Alikhanov formula at the node $t _{j+\sigma}$ is defined as

Lemma 3. Taylor's theorem can be used directly to obtain that identity

Initiating by L2-$1 _\sigma$ formula (2.7) to discretize the time Caputo fractional derivative of (1.1a), leads to

Define the following parameters

Then (2.9) has that equivalent form:

then, the full discrete scheme is to find $\psi_N^{j+1}, \ \phi_N^{j+1} \in V_N^0,$ $j \ge 0, \ \forall \nu \in V_N^0$ such that

where $P_N$ is a projection operator.

3.   Iterative algorithm implementation

Jacobi polynomials $J^{\alpha, \beta }_{i}(x)$ by the aid of Via the hypergeometric function can be (for $\alpha, \, \beta > -1$ and $x \in (-1, 1)$) as ^[29]:

such that the notation $(\cdot)_i$ represents the symbol of Pochhammer. Then, the equivalent three-term recurrence relation can be yielded

where

The Legendre polynomial $L_i (x)$ is a special case of the Jacobi polynomia. This means

The weight function which makes the orthogonality of Jacobi polynomials valid is given as $\omega^{\alpha, \beta}(x) = (1-x)^\alpha(1+x)^\beta$, i.e.,

where $\delta _{ij}$ is the Dirac Delta symbol, and

Lemma 4 (see for example ^[34]). For $\alpha > 0$, one has

We introduce the following rescale functions:

and we write $\wedge(x)$ as $\hat x$. The basis functions selected for the spatial discretization are given by ^[37,28]:

The function space $V_N^0$ can be specified as follows:

The approximate solutions $\psi_N^{j+1}$ and $\phi_N^{j+1}$ may be expressed as

where $\hat \psi_i^{j+1}$ and $\hat \phi_i^{j+1}$ are the unknown expansion coefficients to be determined. Choosing $v = \varphi _i, \ 0 \le i \le N - 2$. Then, the matrix representation of the Alikhanov L2-$1_{\sigma}$ Legendre-Galerkin spectral scheme has the following representation:

where

Lemma 5 (see ^[37,29]). The elements of the stiffness matrix $S$ are given by

where

and $\left\{ {x _r^{ - \frac{\alpha }{2}, - \frac{\alpha }{2}}, \varpi _r^{ - \frac{\alpha }{2}, - \frac{\alpha }{2}} } \right\}_{r = 0}^N$ are Jacobi-Gauss points and their weights with respect to the weight function $\omega^{ - \frac{\alpha }{2}, - \frac{\alpha }{2}}.$ The mass matrix ${ \hat M}$ is symmetric and its nonzero elements are given as

Monitoring $H_{q}^{j+1, r} = H_{q}^{j+1 } (\psi_N^{j+1, r}, \phi_N^{j+1, r}),$ $q = 1, 2$, $r\geq 0.$ Then, the linear system (3.11) can be solved by the following iteration algorithm 1:

4.   Convergence analysis

We fix ${C}$ to be a generic positive constant which may differ from one inequality to another and is independent of $\tau,$ $N$, and $n$. Firstly, the following lemma is devoted to introducing the property of the projector operator $P_N$.

Lemma 6 (see ^[38]). $\forall \, \Psi \in H_{0}^{\frac{\alpha}{2}}(\Omega)\cap H^{s}(\Omega),$ there exists $P_{N}$ such that:

where $\epsilon$ and $s$ are real numbers satisfying $s > \frac{\alpha}{2}$.

The interpolation operator $I_{N}$ achieves the following property:

Lemma 7 (see ^[29]). Suppose that $\Psi \in H^{s}(\Omega)\ (s\geq1).$ Then,

and the constant $C > 0$ is independent of $N$.

Lemma 8 (see ^[30]). For any complex functions $\Psi, \Phi, \psi$ and $\phi$, we have

Lemma 9 (see ^[3]). Assume the existence of an absolute continuous function $\Psi(t)$ in $[0, T]$. Then,

Lemma 10 (Grönwall inequality ^[26,9]). Let $\Psi(t)\geq 0$ be a non-negative function which is locally integrable on $[0, +\infty]$ such that $\phantom {|} _0 D _t ^\beta \Psi(t) \leq \lambda \Psi(t) + b$. Then, we have $\Psi(t) \leq \Psi_{0}E_{\beta}(\lambda t^{\beta}) + b t^{\; \beta} E_{\beta, 1+\beta}(\lambda t^{\beta}),$ where the Mittag-Leffler function $E_{\beta}(z)$ and the generalized Mittag-Leffler function $E_{\beta_1, \beta_2}(z)$ are defined by

Lemma 11 (see ^[26,9]). For $0 < \beta_{1} < 2,$ and $\; \beta_2 \in \mathbb{R},$ we assume that there exists $\mu$ such that $\pi\beta_{1}/2 < \mu < \min(\pi, \pi\beta_{2}).$ Then, also there exists a constant ${\mathrm{C}} = C(\beta_1, \beta_2, \mu)$ such that $|E_{\beta_1, \beta_2}(z)|\leq \frac{{\mathrm{C}}}{1+|z|},$ for $\mu \leq |arg(z)|\leq \pi.$ In addition, if $\beta_{1}\in (0, 1),$ then we have the following properties

By fixing the following notation, and if we denote

then for $1 < \alpha \leq 2$, the semi-norm and the norm are defined by

and for any $u, \ v \in H^{\alpha/2}_{0}(\Omega).$ Then, there exist positive constants $C_1, \ C_2$ such that

The orthogonal projection operator $P_N:H_{0}^{\frac{\alpha}{2}}(\Omega)\rightarrow {V_N^0 }$ satisfies

Lemma 12 (see ^[3]). Let $\Psi(t)$ be any function defined on $\Omega$ and $0 < \alpha < 1$. If $\Psi^{(\sigma)} = \sigma \Psi^{j + 1} + (1 - \sigma) \Psi^j$ then

Lemma 13 ($L2$-$1_{\sigma}$ discrete fractional form of Grönwall inequality ^[24,15]). Suppose that the non-negative sequences $\{\omega^{j}, g^j |j = 0, 1, 2, ...\}$ satisfy $\phantom {|} _0 D _\tau ^\beta\omega^{j+\sigma}\leq \lambda_{1}\omega^{j+1}+\lambda_{2}\omega^{j}+g^{j}.$ Then, there exists of a positive constant $\tau^{*}$ such that

whenever $\tau \leq (\tau^{*})^{\beta} = 1/ (2 \Gamma (2 - \beta) \lambda _1)$ and

4.1.   Semi-discrete form convergence analysis

Theorem 1. Let $\{0 < \beta < 1, \; 1 < \alpha < 2, \; s\geq 1\}$. Assume that the solutions $\{\psi, \phi\}$ and $\{\psi_{N}, \phi_{N}\}$ of (1.1) and (2.9), respectively, are bounded. Thus, satisfying $\{\psi, \phi\} \in H^{1}(I; H_{0}^{\frac{\alpha}{2}}(\Lambda)\cap H^{s}(\Lambda))$, such that the external potential function $V = V(x)$ satisfies $V \in C(\Omega).$ Then, we get

Proof. The variational formulation is derived by taking the inner product of (1.1a) with ${\bf{v}}_1$,

Let $e = \psi-\psi_{N},$ $\zeta_e = \psi-P_{N}\psi$ and $\eta_e = P_{N}\psi-\psi_{N}.$ Then, we get $e = \zeta_e+\eta_e.$ Also, let $E = \phi-\phi_{N},$ $\zeta_E = \phi-P_{N}\phi$ and $\eta_E = P_{N}\phi-\phi_{N}.$ Hence, we get $E = \zeta_E+\eta_E.$ Now, making use of Lemma 6 and in the case of $\alpha\ne \frac{3}{2}$, then the following two estimates are obtained:

Subtracting (4.10) from (2.11), then we get

The orthogonality of $P_{N},$ yields

Taking the inner product of (4.13) with $\eta_e$ and noting (4.6), and choosing the imaginary part of the resulting equation, we get

Similarly, the imaginary part of the error difference equation concerned with (1.1b) and its semi discrete approximation in (2.11) has the following form

Adding (4.15) and (4.16), and noticing that

yields

Define

then assuming the boundness of the exact solutions $\{\psi, \phi\}$ and the approximate solutions $\{\psi_N, \phi_N\}$ for the system (1.1a)-(1.1d) and invoking Lemma 8 give

Similar analysis leads to

Define $G(\varepsilon) = |\varepsilon|^{2}\varepsilon,$ and making use of Cauchy-Schwarz inequality, we deduce

By Lemma 7, we have

Also,

accordingly and after some rather manipulations and invoking Lemma 9, we finally get

An implementation of Lemmas 10 and 11, yields that $\left(\| {\eta_e} \|+ \| {\eta_E} \| \right)^2\leq c_7 N ^{-2 s}$. The other inequality can be achieved in a similar fashion if $\alpha = \frac {3} {2}$ and $0 < \mu < \frac {1}{2}$. □

4.2.   Full-discrete form convergence analysis

Theorem 2 (Convergence of the uniform $L2$-$1_{\sigma}$ – Galerkin spectral scheme). Let $\{\psi, \phi\}$ and $\{\psi_{N}^n, \phi_{N}^n\}$ be solutions of (1.1) and (2.9), respectively, with the condition that the external potential function $V = V(x)$ satisfies $V \in C(\Omega)$ and suppose that the unique solution $\{\psi, \phi\} \in L^{\infty}([0, T;H _0 ^{\frac {\alpha} {2}} (\Omega) \cap H ^s (\Omega))$ is sufficiently regular in temporal and spatial directions, bounded and $\frac{\partial^{\beta}\psi}{\partial t^{\beta}}, \frac{\partial^{\beta}\phi}{\partial t^{\beta}}\in L^{\infty}([0, T;H _0 ^{\frac {\alpha} {2}} (\Omega) \cap H ^s (\Omega))$. Then, a positive constant $\tau^*$ is existed such that when $0 < \tau \leq \tau^*$, the Galerkin spectral scheme (2.9a)–(2.9b) admits a unique solution $\{\psi _{N}^n, \phi_{N}^n\}$ satisfying

where ${C}$ is a positive constant that has no dependence on $n$, $\tau$ and $N$.

Proof. The next variational formula is derived by taking the inner product of (1.1a) with ${\bf{v}}_1$,

Let $e = \psi-\psi_{N},$ $\zeta_e = \psi-P_{N}\psi$ and $\eta_e = P_{N}\psi-\psi_{N},$ then we get $e^{j+\sigma} = \zeta_e^{j+\sigma}+\eta_e^{j+\sigma}.$ Also, let $E = \phi-\phi_{N},$ $\zeta_E = \phi-P_{N}\phi$ and $\eta_E = P_{N}\phi-\phi_{N},$ then we get $E^{j+\sigma} = \zeta_E^{j+\sigma}+\eta_E^{j+\sigma}.$ Using Lemma 1, in case of $\alpha\ne \frac{3}{2},$ we get

Subtracting (4.23) from (2.11), then we obtain

The orthogonality of the operator $P_{N},$ enables one to write

Taking the inner product of (4.26) with $\eta_e^{j+\sigma}$ and noting (4.6) and choosing the imaginary part of the resulting equation, we get

Similarly, the imaginary part of the error difference equation related with (1.1b) and its fully discrete approximation in (2.11) has the following form

Adding (4.28) and (4.29), and noticing that

gives

Proceeding as in the proof of Theorem 1, we finally get

Applying the $L2$-$1_{\sigma}$ discrete fractional form of Grönwall inequality in Lemma 13, the final result (4.21) is achieved directly. Similarly, we can get the result (4.22) when $\alpha = \frac{3}{2}.$ This completes the proof of the theorem. □

The stability analysis can be proved similarly following ^[16].

5.   Numerical experiments

Example 1. Consider the following fractional order Gross–Pitaevskii coupled system:

with homogeneous boundaries

$f_1 (x, t)$ and $f_2(x, t)$ can be deduced by considering the exact solutions

By specifying, $\alpha = \beta+1 = 1.2,$ $1.5,$ $1.8$ and $N = 100$ for $\phi$ and $\psi$, respectively. It can be shown from Tables 1 and 2, that the resulting $L^2$-errors and the corresponding temporal convergence orders support the theoretical results with convergence order close to 2 in case of the smoothness of the solution. A spectral accuracy of convergence is also shown in Figures 1 and 2 by specifying $M = 1600$ for different $\alpha$ and $\beta$.

Example 2. Consider the following Gross¢Pitaevskii system

with the initial conditions

The fractality effects in time with respect to $0 < \beta\leq1$ and in space with respect to $1 < \alpha\leq 2$ affect directly the shapes and stability of the soliton solutions. This can clearly reflected in the numerical solutions given in Figures 3, 4, 5, 6 and 7. The performance of the solutions is represented in the following manner as can be observed from the different experiments. When $\alpha \neq 2$ and $\beta \neq 1$, the collision is not elastic and so an influence on the shape of solutions is observed. It can be also observed that the decrease in $\alpha$ values could make the shape of solitons changes faster. As observed in Figures 5, 6, 7 and 8, we conclude that different decay properties in the time direction coming from distinct selections of the fractional-order parameters $\beta$. These characteristics can be used in physics to tunable the sharpness of the space-time fractional Gross¢Pitaevskii equations by changing the space fractional order $\alpha$ and the time fractional order $\beta$ without changing the nonlinearity and dispersion effects.

6.   Conclusions

A high order (second temporal order and a spatial spectral accuracy) convergent numerical approach has been investigated for solving a system of fractional order coupled Gross Pitavskii equations. An algorithmic implementation of the scheme is given to simplify its numerical implementation. A theoretical analysis of the scheme shows unconditional convergence towards the true solution. This is also proven by giving some numerical experiments.

Conflict of interest

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

M. A. Zaky and A. Aldraiweesh extend their appreciation to the Distinguished Scientist Fellowship Program (DSFP) at King Saud University (Saudi Arabia).