A Convergent Three-Step Numerical Method to Solve a Double-Fractional Two-Component Bose–Einstein Condensate

: This manuscript introduces a discrete technique to estimate the solution of a double-fractional two-component Bose–Einstein condensate. The system consists of two coupled nonlinear parabolic partial differential equations whose solutions are two complex functions, and the spatial fractional derivatives are interpreted in the Riesz sense. Initial and homogeneous Dirichlet boundary data are imposed on a multidimensional spatial domain. To approximate the solutions, we employ a ﬁnite difference methodology. We rigorously establish the existence of numerical solutions along with the main numerical properties. Concretely, we show that the scheme is consistent in both space and time as well as stable and convergent. Numerical simulations in the one-dimensional scenario are presented in order to show the performance of the scheme. For the sake of convenience, A MATLAB code of the numerical model is provided in the appendix at the end of this work.


Introduction
In this work, p ∈ N represents the number of spatial dimensions, and T ∈ R + defines a period of time. We agree that I n = {1, 2, ..., n} and I n = I n ∪ {0}, for each n ∈ N, and assume that a i and b i are real numbers with the property that a i < b i , for all i ∈ I p . For convenience, we use the notations Ω = Π p i=1 (a i , b i ) ⊆ R p and Ω T = Ω × (0, T) and agree that Ω and Ω T denote the closures of Ω and Ω T , respectively, in the Euclidean metric topology of R p . We assume that ψ 1 : Ω T → C and ψ 2 : Ω T → C are sufficiently smooth functions and employ the conventions x = (x 1 , x 2 , . . . , x p ) ∈ Ω and i = √ −1. Throughout this manuscript, Γ will represent the usual Gamma function that extends factorials.
Definition 1 (Podlubny [1]). Suppose that f is a real function defined in all of R and assume that n ∈ N ∪ {0} and α ∈ R are such that n − 1 < α < n. We define the fractional derivative in the Riesz sense of order α of f at the point x ∈ R as Definition 2. Let ψ : Ω T → C be a function and let i ∈ I p . Suppose that α and n are a real and a natural number, respectively, satisfying the properties in the previous definition. Whenever it exists, the space-fractional partial derivative in the sense of Riesz of order α of ψ with respect to the variable x i at (x, t) ∈ Ω T is defined as Here, we extend ψ by allowing ψ ≡ 0 outside of Ω. If these fractional derivatives exist for all i ∈ I p , then we define In this manuscript, we will use D, β 11 , β 12 ,β 22 , and λ to represent constant real numbers, and we employ P : Ω → R to denote a real function. Throughout, we will assume that α 1 , α 2 ∈ R satisfy the properties 1 < α 1 ≤ 2 and 1 < α 2 ≤ 2. We let φ 1 : Ω → C and φ 2 : Ω → C be sufficiently regular functions. In this work, we tackle the numerical solution of the following double-fractional initial-boundary value problem ruled by two nonlinear coupled Gross-Pitaevskii-type equations: This system is a two-component Bose-Einstein condensate that considers the presence of an internal atomic Josephson junction. In this physical context, ψ 1 and ψ 2 denote the two-component condensate wave functions.
There are systems that describe single-component Bose-Einstein condensates. As examples, there are works that investigate round, symmetric, and central vortex states in rotating Bose-Einstein condensates [2], considering only a parameter that characterizes the interaction between particles and a constant that corresponds to the angular speed of the laser beam. As in the present work, a real-valued external trapping potential is employed in those models. From the numerical point of view, there are already works that tackle the numerical solution of single-component fractional Nose-Einstein condensates.
As an example, there are papers that computationally studied the ground states of space-fractinoal nonlinear Schrödinger-Gross-Pitaevskii equations with a rotation term and non-local nonlinear interactions [3]. However, we must clarify that the singlecomponent non-fractional case has been more studied that the double-component fully fractional scenario [4,5], possibly due to the mathematical difficulties in the double-component and fractional scenario. Clearly, this case is physically more relevant as it describes the dynamics of mixtures Bose-Einstein condensates that have been experimentally observed in the laboratory [6][7][8][9][10].
The following are crucial definitions and properties to discretize the fractional derivatives in the initial-boundary value problem in Equation (4).
(a) If k = 0, then g Lemma 2 (Wang et al. [12]). Let f : R → R be a function in the set C 5 (R), which has integrable derivatives to the fifth order on R.
for almost all x ∈ R.
The field of fractional calculus has witnessed a vertiginous development in recent years, partly due to the vast amount of potential applications in the physical sciences [13]. Today, many different fractional derivatives and integrals have been proposed. For example, some of the first fractional derivatives introduced historically in mathematics were the Riemann-Liouville fractional derivatives [14], which generalized the classical integer-order derivatives with respect to specific analytical properties [1].
The fractional derivatives in the senses of Caputo, Riesz, and Grünwald-Letnikov are also extensions of the traditional derivatives of integer order. These fractional operators are nonequivalent in general, and various applications of all of them have been proposed in science and engineering [15,16]. For instance, some reports have provided theoretical foundations for the application of fractional calculus to the theory of viscoelasticity [17], while others have proposed possible applications of fractional calculus to dynamic problems of solid mechanics [18], continuous-time financial economics [19,20], Earth system dynamics [21], mathematical modeling of biological phenomena [22], and the modeling of two-phase gas/liquid flow systems [23], to mention some potential applications.
However, it is important to recall that Riesz-type derivatives may be the only fractional derivatives that have real physical applications. This is due to a well-known result by Tarasov, which establishes that Riesz fractional derivatives result from systems with longrange interactions in a continuum-limit case [24]. In turn, systems consisting of particles with long-range interactions are useful in statistical mechanics, thermostatics [25,26], and the theory of biological oscillator networks [27].
From the mathematical point of view, many interesting avenues of investigation have been opened by the progress in fractional calculus. Indeed, the different fractional derivatives have found discrete analogues, which have been used extensively in the literature. As examples, Riesz fractional derivatives have been discretized consistently in various fashions using fractional-order centered differences [11,28] and weighted-shifted Grünwald differences [29,30].
Clearly, those discrete approaches have been studied to determine their analytical properties, and they have been used extensively to provide discrete models to solve Riesz space-fractional conservative/dissipative space-fractional wave equations [31], a Hamiltonian fractional nonlinear elastic string equation [32], an energy-preserving double fractional Klein-Gordon-Zakharov system [33], and even a Riesz space-fractional generalization with generalized time-dependent diffusion coefficient and potential of the Higgs boson equation in the de Sitter space-time [32], among other complex systems.
On the other hand, Caputo fractional derivatives have been discretized consistently using various criteria. For instance, some high-order L2-compact difference approaches have used to that end [34], as well as L1 formulas [35,36] and L1-2 methodologies [37]. Using those approaches, various numerical schemes have been proposed to efficiently solve Caputo time-fractional diffusive and wave differential equations [38][39][40]. Various potential applications of these systems have been reported in the sciences [25,41].
From the analytical point of view, the literature offers a wide range of reports that focus on the extension of integer-order methods and results for the fractional case. For example, there are various articles that tackle the existence, uniqueness, regularity, and asymptotic behavior of the solution for the fractional porous medium equation [42], nonlinear fractional diffusion equations [43], nonlinear fractional heat equations [44], the Fisher-Kolmogorov-Petrovskii-Piscounov equation with nonlinear fractional diffusion [45], fractional thinfilm equations [46], and the fractional Schrödinger equation with general nonnegative potentials [47].
From a more particular point of view, the fractional generalization of the classical vector calculus operators (that is, the gradient, divergence, curl, and Laplacian operators) has also been an active topic of research developed from different approaches. Some of the first attempts to extend those operators to the fractional scenario were proposed in [48,49] using the Nishimoto fractional derivative.
These operators were used later on in [50] to provide a physical interpretation for the fractional advection-dispersion equation for flow in heterogeneous porous media (see [51] and the references therein for a historic account of the efforts to formulate a fractional form of vector calculus). More recently, a new generalization of the Helmholtz decomposition theorem for both fractional time and space was proposed in [52,53] using the discrete Grünwald-Letnikov fractional derivative. As in the present manuscript, the authors of [52,53] considered different derivative orders assuming non-homogeneous models and non-isotropic spaces.
This manuscript is organized as follows. Section 2 is devoted to present the discrete nomenclature and the numerical model to solve the continuous problem under investigation in this work. We introduce suitable discrete norms and recall a useful theorem by Desplanques. With that result at hand, we rigorously establish the properties of existence and uniqueness of numerical solutions. Section 3 is devoted to establishing the most important computational features of our numerical scheme. More concretely, we show that the discrete model proposed in this work has second order consistency under suitable conditions on the smoothness of the solutions.
We also rigorously prove the stability property of the scheme along with its quadratic order of convergence in both time and space. In turn, Section 4 is devoted to show some illustrative simulations that exhibit the stability and the convergence of the method proposed in this work. We will also provide numerical evidence that the finite-difference scheme proposed here has a quadratic order of convergence in both space and time. Finally, we close this manuscript with some concluding remarks. For the sake of convenience to the reader, an appendix is provided at the end of this paper, in which we present the numerical scheme used to produce the simulations.

Numerical Method
Throughout this manuscript, we let M i , N ∈ N, for each i ∈ I p , and fix the spatial and temporal step-sizes h i = (b i − a i )/M i and τ = T/N, for each i ∈ I p . We will consider regular partitions of the intervals [a i , b i ] and [0, T] represented as respectively. Let J = ∏ p i=1 I M i −1 and J = ∏ p i=1 I M i , and convey that x j = (x 1,j 1 , . . . , x p,j p ), for each j = (j 1 , . . . , j p ) ∈ J. Throughout this work, we will employ the symbol (u n j , v n j ) to denote a numerical approximation to (U n j , V n j ), where U n j = ψ 1 (x j , t n ) and V n j = ψ 2 (x j , t n )), for each (j, n) ∈ J × I N . Let us agree that P j = P(x j ), for each j ∈ J. We will employ the symbol ∂J to represent the boundary of J, which is defined as the set of all multi-indexes j ∈ J that satisfy x j ∈ ∂Ω. Definition 4. We introduce the following average operators, for each α ∈ (0, 1) ∪ (1, 2], (j, n) ∈ J × I N−1 and w = u, v: Under the same conditions, we introduce the difference operators In the present manuscript, we convey that Moreover, for the sake of convenience, we also introduce the discrete fractional gradient vector x p w n j ). (15) With this nomenclature at hand, the following numerical method will be employed in this work to approximate the solution of the initial-boundary value problem Equation (4) on the set Ω × [0, T], for each (j, n) ∈ J × I N−1 : One may readily check that the general iterative formula of this numerical model makes it a three-step implicit model. On the other hand, to determine the approximations at the time t 1 , we need to employ the initial conditions u 0 Substituting these formulas into the general iterative equations at time t 0 , we obtain and For the remainder of his manuscript, we define h = (h 1 , . . . , h p ) and h * = ∏ p i=1 h i . We will employ V h to represent the collection of all complex functions on the spatial mesh {x j : j ∈ J} ⊆ R p . Finally, for each w ∈ V h and all j ∈ J, let w j = w(x j ).
for any u, v ∈ V h . We define the function · 2 : V h → R as the Euclidean norm obtained from ·, · . If u, v ∈ V h , then we agree that uv ∈ V h represents the pointwise multiplication of u and v. Moreover, |u| ∈ V h will denote the absolute value of the complex function u.
Lemma 3 (Desplanques [54]). Let A be a square complex matrix. If A is strictly diagonally dominant, then A is invertible.
is any set of initial conditions, then the finitedifference model Equation (16) has a unique solution.
Proof. The proof will be provided for the one-dimensional scenario only in view that the case for higher dimensions is similar. Note that, beforehand, the approximations u 0 , v 0 , u 1 and v 1 are provided explicitly by the initial conditions and the formulas Equations (17) and (18). Proceeding by induction, suppose that u n , v n , u n−1 , and v n−1 were calculated, for some n ∈ I N−1 . It is important to observe now that the first equation of the numerical method can be equivalently rewritten as where for each (j, n) ∈ J × I N−1 . It follows that the first identity of the numerical scheme Equation (16) can be alternatively expressed in vector form as In this expression, the square complex matrix A = (a ij ) is defined by Notice now that the diagonal entries of A satisfy On the other hand, the numbers in the real sequence (g It follows that the matrix A is non-singular by the previous lemma, which implies then that there exists a unique solution u n+1 of the vector equation Equation (24). In similar fashion, we can readily prove the existence of the vector v n+1 after noticing that the second difference equation of the discrete system Equation (16) may be rewritten in vector form as where, for each (j, n) ∈ J × I N−1 , In this expression, B is the square complex matrix given by One may readily check that this matrix is strictly diagonally dominant; therefore, Equation (28) has a unique solution v n+1 . The conclusion of this theorem readily follows now by mathematical induction.

Proof.
To prove these results, we require Taylor's theorem, Lemma 2, and the smoothness assumptions on the functions ψ 1 and ψ 2 . There exist nonnegative constants C 1,k , C 2,k , C 3,k,i , C 4,k , and C 5,k.l , which are independent of τ and h, for each i ∈ I p and k, l ∈ I 2 , with the property that for each (j, n) ∈ J × I N−1 . The inequality in the conclusion of this proposition follows now applying the triangle inequality for an appropriate number C 1 , which depends only on C 1,k , C 2,k , C 3,k,i , C 4,k , and C 5,k.l for each i ∈ I p and k, l ∈ I 2 , and the model parameters. This implies, in particular, that C 2 is independent of h and τ, as desired. The constant C 2 is obtained applying the same arguments around the node (x j , t n ) instead of (x j , t n+ 1 2 ).
We tackle now the stability and the convergence of the numerical scheme Equation (16). We require the next technical results. [55]). Suppose that (ω n ) N n=0 and (ρ n ) N n=0 are sequences of non-negative numbers. Let C ≥ 0 satisfy the inequality
The stability of the discrete method Equation (16) will require us to consider two sets of initial conditions, namely (φ 1 , φ 2 ) and (φ 1 ,φ 2 ), where all the initial data are complex functions, which are defined on Ω. The respective solutions will be denoted by (u, v) and (ũ,ṽ). This means that (u, v) satisfies Equation (16), while (ũ,ṽ) satisfies the problem for each (j, n) ∈ J × I N−1 .
Lemma 8. Let (u, v) and (ũ,ṽ) be solutions of Equation (16), respectively, and let n = u n −ũ n and ζ n = v n −ṽ n , for each n ∈ I N . For each (j, n) ∈ J × I N−1 , let us define Then, there exists a constant C ≥ 0, which is independent of h and τ, with the property that Proof. The result is straightforward.

Theorem 3 (Nonlinear stability).
Consider the solutions (u, v) and (ũ,ṽ) of Equation (16) corresponding to initial conditions (φ 1 , φ 2 ) and (φ 1 ,φ 2 ), respectively. For each n ∈ I N , define n = u n −ũ n and ζ n = v n −ṽ n . If (C + λ)τ < 1 2 , then there is a nonnegative constant C that is not dependent on τ and h, such that Proof. Beforehand, take C as in the last result. Note that the existence of (u, v) and (ũ,ṽ) is guaranteed by Theorem 1. It is clear that ( , η) satisfies the following discrete problem, for each (j, n) ∈ J × I N−1 : Here, P n j and Q n j are as in Lemma 8. Then, there is a nonnegative real number C such that, for each n ∈ I N−1 , max{ P n 2 2 , Q n 2 2 } ≤ C n−1 2 2 + ζ n−1 2 2 + n 2 2 + ζ n 2 2 .
The number C does not depend on h and τ. Take the first identity in Equation (52) and calculate its inner product with 2µ (1) t n . On both sides of the result, next calculate the imaginary part. Similarly, consider the second identity of Equation (52) and take the inner product with 2µ (1) t ζ n . Calculate the imaginary part on both sides. After applying the identities after Theorem 6, we readily obtain that δ t n 2 2 = 2 Im P n , µ t n + 2λ Im µ (1) t ζ n , µ t n , ∀n ∈ I N−1 , After adding those two results side by side and bounding from above, we notice that there exists a nonnegative number C that is not dependent on τ and h, with Im( P n , µ t n + Q n , µ t ζ n ) Subtract the term τ(C + λ)( m+1 2 2 + ζ m+1 2 2 ) from both ends of this chain of inequalities, group the common terms, and we notice that (C + λ)τ < 1 2 . As a consequence, for each m ∈ I M−1 . The conclusion is reached then using Lemma 4 with the constants ω n = n 2 2 + ζ n 2 2 and ρ n = (1 + (C + λ 4 )τ) 0 2 2 + ζ 0 2 2 , for each n ∈ I N−1 .
Finally, we investigate now the convergence of the finite-difference method Equation (16). As mentioned previously, (u, v) represents a solution of the numerical model, while (U, V) = (ψ 1 , ψ 2 ) denotes a solution of the continuous problem. It is clear that (U, V) satisfies the discrete initial-boundary value problem Here, ρ n j and σ n j denote local truncation errors. Under the regularity assumptions on the solutions U and V stated in Theorem 2, there is a nonnegative number C 0 , such that ρ n 2 , σ n 2 ≤ C 0 (τ 2 + h 2 2 ). This number is independent of τ and h. This fact will be employed in the proof of the following theorem. (4), such that U, V ∈ C x,t (Ω T ) and P ∈ C(Ω). If τ is sufficiently small, then the corresponding solution of the finite-difference scheme Equation (16)

Theorem 4 (Convergence). Let (U, V) be a solution of Equation
Proof. The proof is analogous to the previous result. Beforehand, let n = u n − U n and ζ n = v n − V n , for each n ∈ I N . It is clear then that the pair ( , ζ) satisfies where, for each (j, n) ∈ J × I N−1 , and Proceeding now as in the proof of Theorem 3, we obtain a nonnegative number C , which does not depend on τ and h, with the property that Equation (53) is satisfied. Applying the same steps used to obtain Equations (54) and (55), we may readily obtain the following identities, valid for each n ∈ I N−1 : δ t n 2 2 = 2 Im P n , µ t n + 2λ Im µ (1) t ζ n , µ t n + 2 Im ρ n , µ n , δ t ζ n 2 2 = 2 Im Q n , µ t ζ n + 2λ Im µ (1) t n , µ t ζ n + 2 Im σ n , µ t ζ n .
From here, it is possible to check that there are numbers C 3 , C 4 ≥ 0 independent of τ and h, such that (Im P n , µ t n + Im Q n , µ t ζ n + Im ρ n , µ t n + Im σ n , µ t ζ n ) Subtract the term τ(C + λ + 1 2 ) m+1 2 2 + ζ m+1 2 2 from both ends of this series of inequalities. Notice now that Equation (45) is satisfied, for each m ∈ I N−1 , letting C 0 = C 4 , and Lemma 4 guarantees now that where C 5 = C 3 e C 0 T . From this, it follows that n 2 , ζ n 2 ≤ √ C 5 (τ 2 + h 2 2 ), for each (j, n) ∈ J × I N . The conclusion of this proposition readily follows now.

Illustrative Simulations
We provide illustrative simulations to show the main properties of Equation (16). In particular, we provide a numerical study of the convergence of the scheme, and we confirm the validity of the conclusion of Theorem 4. The experiments provided in this section make use of the Algorithm provided in Appendix A, which is a MATLAB implementation of Equation (16).

Example 2.
We numerically studied the convergence of the scheme Equation (16), using the parameters in Table 1 along with the functions V, φ 1 and φ 2 defined in the paragraph before Example 1. In the following experiments, we fix T = 0.5 and α 1 = α 2 = 1.5, though similar results have been obtained for other values of these parameters (we did not include them to avoid redundancy). The exact solution of the problem is calculated using τ = 5 × 10 −5 and h = h 1 = 0.02. Let (u, v) be the solution of the scheme Equation (16) at the time t N , obtained using the computer parameters τ and h. Under these circumstances, let us define the complex constants τ,h (x j ) = ψ 1 (x j ) − u N j , for each j ∈ J. In turn, we define the temporal and spatial convergence rates, respectively, as With these conventions, Table 2a,b provides numerical studies of convergence in time and in space, respectively, for the problem described in this example. Various values of τ and h were chosen to that end. However, in any case, the results seem to support that the scheme Equation (16) is quadratically convergent, in agreement with Theorem 4.

Conclusions
We proposed a three-level discrete model to solve a double-fractional coupled Bose-Einstein system considering Riesz spatial fractional operators along with two different orders of differentiation. The system is a fractional generalization of the two-component Gross-Pitaevskii system investigated in physics. The solutions of the continuous system are complex; therefore, the computational complexity for solving the system is a shortcoming to be considered.
To alleviate this problem, a three-level scheme was designed to solve our continuous model. This scheme is linear and explicit, which implies that the computational implementation is relatively easy to propose. Indeed, a numerical algorithm to solve this system is presented in the appendix of this work. The code is provided in MATLAB for the sake of convenience, and the nomenclature used to describe it follows that of the theoretical description.
In this manuscript, we proved the existence and uniqueness of the solutions using arguments of linear algebra. We also established that the algorithm is spatially and temporally quadratically consistent. We proved that the finite-difference scheme is spatially and temporally quadratically convergent to the solution of the differential problem. Computational simulations are presented in this work to illustrate this fact.
The description of these recursive formulas will be described next. The initial steps Equations (17) and (18) will be observed also in our implementation. Moreover, the computational implementation will make use of the fact that the coefficients (g (α) k ) ∞ k=−∞ can be expressed in alternative form as g (α) 0 = Γ(α + 1) Γ(α/2 + 1) 2 (A9) and g (α) The following conventions and definitions will be observed in our code: Input: • a1= a 1 ,