Murat Gunestas

Murad Mehmet

Duminda Wijesekera

Detecting Ponzi and Pyramid Business Schemes in Choreographed Web Services

Keywords: Web services, choreographies, Ponzi schemes, pyramid schemes

de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Businesses are increasingly invoking dynamic services and generating content-specific operations among choreographed web services, thereby creating service interdependencies between web services. These dynamic service interdependencies can be exploited to create new misuse activities. Some exploit the infrastructural dependencies of the services themselves while others use the infrastructural dependencies to create illegal business schemes. This paper focuses on detecting Ponzi and pyramid investment schemes created to defraud unsuspecting investors.

Illegal business schemes are difficult to detect because they are similar to legal business schemes at a microscopic level and are apparent only at the macroscopic level. Thus, they can elude local monitoring of web Table 1. Pyramid scheme.

Level

Payments of $400 (#) transactions. Ponzi and pyramid schemes [START_REF]Security and Exchange Commission, Ponzi schemes, Washington[END_REF][START_REF] Valentine | Pyramid schemes[END_REF] are difficult to differentiate from multilevel marketing schemes that either run their own businesses or invest in others. The basic dynamic of these two schemes is to rob Peter to pay Paul [START_REF] Zuckoff | Ponzi's Scheme: The True Story of a Financial Legend[END_REF]. The Ponzi scheme is named after Charles Ponzi [START_REF] Zuckoff | Ponzi's Scheme: The True Story of a Financial Legend[END_REF]. For many years, he collected money and promised returns within 90 days to investors who enrolled other investors in the scheme. Ponzi paid out the early investors using funds invested by late joiners. Other than running this scheme, Ponzi neither ran a business nor invested in other businesses and, consequently, neither incurred a profit nor a loss.

1 ($100 × 3 = $300) 1 × # 1 × # 1 × # 2 ($30 × 9 = $270) 3 × # 3 × # 3 × # 3 (
Numerous Internet-based Ponzi schemes are currently being investigated [START_REF] Savvysugar | Help! My friend is part of a pyramid scheme[END_REF][START_REF] Stewart | South Africa investigates an alleged Ponzi scheme[END_REF]. In 2006 alone, 25,000 web sites suspected of running Ponzi schemes were shut down by the U.S. Securities and Exchange Commission [7].

A classic pyramid scheme, shown in Table 1, also uses the same principle. The orchestrator who originates the scheme promises top-level investors large returns on their investments, as do the recruited investors to their potential recruits. The example in Table 1 is organized as a four-level payment scheme with a span of three, i.e., only up to four levels of ancestral recruiters profit from investments, and each recruit at every level recruits three others. A Level 1 investor recruits three others and receives $100 per recruit. Each recruit is expected to recruit three others, thereby building a recruit tree. The Level 1 investor is paid $30 (not $100) per recruit at Levels 2, 3 and 4; and does not profit from investors beyond Level 4.

The orchestrator considers the recruiting activity to be complete when he receives $400 from an investor. Therefore, the orchestrator pays his recruiters $100 for the first parent or $30 for the three immediate ancestral recruiters, paying at most $190 on an investment of $400. Consequently, an investor makes $3,810 ($2,430 + $810 + $270 + $300) on his investment of $400 to his promoter if he successfully recruits three other investors and they successfully recruit three other investors all the way down to Level 4. The scheme grows quickly, but is unsustainable in the long term and many promoters lose their investments because of the difficulty of attracting new recruits [START_REF] Valentine | Pyramid schemes[END_REF]. Financial institutions and their business partners are increasingly engaging service-oriented architectures such as dynamic brokering, creating dynamic service interdependencies that can be exploited to orchestrate illegal business practices. In order to detect such illegal schemes, it is necessary to have a comprehensive non-repudiable perspective of complex multi-party transaction models. A multi-party communication can arise in two possible ways among dynamic web services. In the first, a static communication pattern is specified and all participants follow this previously-known pattern. In the second, web services discover and transact with other services, dynamically creating choreographies that were unknown a priori. Consequently, in order to discover illegal activity, it is necessary to identify illegal business transactions using both types of choreographies.

Evidence Generation Framework

This section briefly reviews the evidence generation framework (EGF) of Gunestas, et al. [START_REF] Gunestas | An evidence generation model for web services[END_REF], which is used to track business-level choreographic misuse. The EGF has three layers (Figure 1). The pairwise evidence generation service (Level 1) generates evidence for pairwise interactions between web services. The evidence derivation service (Layer 2) derives facts from available pairwise evidence in order to refute or justify claims of agreement violations between communicating partner services. The comprehensive evidence generation service (Layer 3) generates instances of requested choreographies from Layer 2 and Layer 3 data.

The EGF provides online evidence generation and management capabilities to other web services as a web service itself. In order to use the EGF, other web services (called member services) should integrate the EGF with themselves using a centralized service access point. Thereafter, the EGF acts as a trusted third party. As a service, the EGF receives and retains service requests and responds in a cryptographicallysecure manner, retains the correspondence in secure repositories, and provides the correspondence for dispute resolution and forensic investigations. The EGF provides "evidence adapters" for all requests.

Gunestas, et al. [START_REF] Gunestas | An evidence generation model for web services[END_REF] have constructed a prototype implementation of Layer 1 and several protocols based on one-way and request-response message exchange patterns. Currently, the EGF provides evidence for non-repudiation, fairness and timeliness using digital signatures to provide proofs of receipt and delivery, to link a message to its creator/sender and to provide message integrity. For accountability, the EGF uses fair non-repudiation mechanisms that utilize trusted third parties (TTPs). Some fair exchange protocols (e.g., [START_REF] Kremer | An intensive survey of nonrepudiation protocols[END_REF]) do not use TTPs because they assume that the participants have prior knowledge of the message contents; but they are not used in the EGF because web services may not always know the expected message content. Timeliness is required because of the time-sensitive nature of most business transactions. Evidence records are based on the time observed at TTPs. EGF servers gather pairwise transactional evidence that flows between sender and receiver web services, employing inline TTPs that use the Simple Evidence Layer Protocol (SELP) or offline TTPs that use the Optimistic Evidence Layer Protocol (OELP) [START_REF] Herzberg | The Delivery and Evidence Layer[END_REF]. SELP and OELP are used by end-points to obtain non-repudiable evidence by engaging a specific message format and digital signatures.

Evidence of Observed Interactions

Web services use many kinds of messages (e.g., one-way messages and request-response messages) in order to choreograph business processes between themselves that correspond to the four proposed WSDL operation types (in-only, out-only, in-out and out-in). An external observer who is ignorant of the business processes can observe only the one-way and request-response message exchange patterns that are formalized in Definition 1.

Definition 1 (Message):

A web services message consists of three components: (i) mandatory fields corresponding to the sender, receiver and time, where the first two are URLs and the last is chosen from a set T ; (ii) optional fields corresponding to a finite set of attributes from a set A; and (iii) message content consisting of strings from an alphanumeric set C.

In this paper, the symbol | denotes string concatenation and encA(r) denotes the string obtained by encrypting r using A's key. Furthermore, if m is a message, then m.a denotes the value of its attribute a. For example, m.time is the value of the timestamp of m.

Definition 1 establishes the notation for describing the messages used to extract knowledge about externally-observable facts pertaining to choreographies. Because different choreographic specifications may select different labels for their identifier fields, in order to address naming convention problems, we use XPATH expressions to specify ID values. Furthermore, because any fabricator can produce messages, we rely on cryptographically-secure messages to ascertain reliable evidence. The messages are collected to derive "evidence objectives" -claims that are to be substantiated or refuted using the collected evidence, e.g., message origin, message properties or the intended recipient. This evidence is generated from cryptographically-secure messages. However, certain objectives such as evidence of delivery and evidence of non-availability may require messages to be signed by a TTP.

Definition 2 (Primitive Evidence Objectives):

The three primitive evidence objectives are: (i) Evidence of Origin: message m with origin A and content r|sigA(r) from A to B is said to provide evidence of origin; (ii) Evidence of Delivery: message m with content ack|sigT T P (ack|m), where TTP is a trusted third party or content ack|sigB(ack|m) and B is m.recipient is said to provide evidence of delivery; and (iii) Message Evidence: evidence of a message m is a pair (m 1 , m 2 ), where m 1 is the evidence of origin and m 2 is the evidence of delivery of m.

According to Definition 2, cryptographic evidence is required from a web service of a trusted third party for claims of origin and delivery. Interested readers are referred to [START_REF] Gunestas | Forensic web services[END_REF][START_REF] Herzberg | The Delivery and Evidence Layer[END_REF] for details about how the evidence of origin and evidence of delivery can be collected at TTPs using WS-Evidence messages that are generated via non-repudiation protocols.

The message evidence (ME) is stored in the form of log records in the EGF as described in [START_REF] Gunestas | Forensic web services[END_REF]. Because log records may contain large volumes of data, message evidence indices (MEIs) are used to refer to messages. Table 2 shows a sample index table, where the first column is an index for stored packets that have the attributes of time, sender, message string

Evidence of Global Misuse

Mining choreographies that are created due to message content from external observations require linkage parameters, which can be derived from externally-invisible message content; this renders them not very helpful to external monitors and auditors. One opportunity to obtain this information is when a victim makes a complaint. Thus, the following method can be used to detect a Ponzi-like scheme: Accept a victim complaint.

Examine the content of specimen records involving promotion or investment messages provided by the victim.

Determine the parameters in the evidence that can be linked.

Detect choreographies and design them in the dynamics of the algorithms.

Create the algorithms.

Run the algorithms in the appropriate order, possibly running more than one algorithm when required.

Collect a comprehensive set of evidence to determine if the scheme is illegal and its effect on the network.

Broadcast an alert to current and potential victims. The method described above works for hierarchical schemes. Different methods would be applied to other schemes. We apply our heuristic method to specify several algorithms for Ponzi-like schemes.

Ponzi Schemes over Web Services

Typical Ponzi-like schemes have three types of actors: a malicious investment service acting as the orchestrator and investor services acting as promoters or victims (depending upon their investments and return rates). Figure 2 illustrates how these actors collaborate to spread the financial scheme over a network of web services. The investment company, InvComp (Orchestrator), promotes A (Promoter) to recruit B into its scheme by promising a quick return on investment and encourages B to promote the scheme to other potential investors. This promoting activity (using promote messages) may not necessarily be observed in the records because promoters may choose other means to convince investors. If A invests (using an invest message) in InvComp, then we say that A has been recruited. A now starts promoting InvComp to other investors in order to get a quick return on his investment and in the process recruits B. We recognize that B has been promoted by A because of the reference value in the content field of the invest message sent by B to InvComp. In accordance with the return policy of the scheme, InvComp makes a payment to A (pay message). The choreographies between the investment company, recruiters and recruits spread in an investor web service network. When a recruit cannot attract enough investors, he loses his money (Victim). Figure 2 shows the complete and incomplete recruit choreographies.

We do not assume that we know the global scheme when mining. Instead, we assume that either a promoter web service identity or an invest message is submitted to a law enforcement agency by a victim. Following the heuristic method presented above, it is possible to find invest, pay and promote messages that contain attributes that refer to each other, demonstrating the collaboration in a pervasive manner.

Pattern Discovery

This section shows how to discover the patterns that help create comprehensive evidence of illegal business schemes. Following the heuristic method described above, we assume that a victim brings an invest message that contains a promoter web service. An algorithm executed on the evidence repository could show that the promoter web service has received pay messages from other web services and has sent an invest message to the same web service of the alleged investment company. Table 3 shows sample records corresponding to the misuse scheme.

The records in Table 3 show the pattern that keeps the fraudulent activity alive -invest messages are linked by sender header fields and promoter content fields. In other words, every promoter web service in an invest message is paid right after the invest message. The message evidence can be correlated to conclude that the promoter, victim and orchestrator may be involved in a hidden recruit choreography.

In cases where the promoting activity does not involve a web service message, the message or choreography patterns can be included as part of the recruiting activity. We call this content-based choreography "recruit." Pattern I in Table 4 is the signature for the "rob Peter to pay Paul" activity. Pattern II is used as the link between the recruiter and recruit, enabling the mining of recruit paths to create recruit trees from MEI records. For simplicity, we define patterns succinctly; however, more complex patterns may include pay messages, which increases the complexity of queries.

Table 3 shows that investors assume a promoter role within a subsequent recruit choreography pattern, thus creating a recursive investment scheme. For example, note that Investor B sends an invest message to InvComp at Time 45. At Time 67, Investor C makes a reference to B through his investment message and B receives a payment from InvComp some time later (pay message at Time 76). The same choreography can be observed between C and other subsequent investors, revealing that they recruited other investors shown in later records. Consequently, in order to detect such a scheme, it is necessary to recognize the recursive investment and payback schemes. The recursive scheme creates a recruit tree that joins a recruiter to all his recruits. By traversing a path in a recruit tree from a chosen (victim) node to the root of the tree (say recruit paths), it is possible to identify the orchestrator. The path Victim-C-B-A in Figure 2 is such a recruit path.

In this paper, we use the following notation to specify recruit trees formally. Let k be an integer that denotes the fanout of the recruit tree. Then, all finite sequences of {0, . . . , k -1} are used as identifiers for web service nodes. We use the notation 2 k<w to denote the set of finite subsequences of {0, . . . , k -1}. For example, all binary sequences can be used to index trees with fanout 2, where the left child of node x p is x p0 and the right child is x p1 , where p is a finite sequence of integers {0,1}, i.e. 2 <w . The notation p<q denotes that p is a subsequence of q where p, q ∈ k <w . The length of p ∈ k <w is denoted by |p|. Now suppose p ∈ k <w where |p|= m and p =<p o , . . . , p m-1 >. Then, the i th ancestors of p for i ≥ 1 are given by ancestor(i)=<p o , . . . , p m-i >. (iv) Define a recruit tree to be active* = LFP(f, m, E) where function f is defined in items (i) and (ii); the message m is defined in item (i); and the set of message equations E is defined in items (i) and (ii).

We denote the class of Ponzi schemes of fanout k and depth l and attribute equations E as Ponzi(k, l, E), where E is the collection of equations in Definition 3. Definition 3 provides a generic definition for Ponzi-like schemes where the number of recruits employed by any recruiter is limited to an integer k and the number of ancestors deriving a payback from the recruitment is at most an integer l.

The web service nodes are numbered by strings from {0, . . . , k -1}, resulting in trees where every node has at most n levels. Thus, the parameter p in Definition 3 denotes a path with |p| elements in such a tree. Step 0 with the empty string φ represents the recruiter in item (i) in Definition 3. Item (ii) assumes that the tree is defined up to a path p of length n and finds its next level. This step is the sequential composition of two steps. First, P p sends messages to each of its children to invest. Next, each of these children invest in I, followed by I paying the ancestors of these children. The ancestors who are paid back are limited to at most l generations. Item (iii) in Definition 3 collects all possible paths that extend the tree to the next level n + 1. Item (iv) collects all the sub-trees with depth n + 1.

Detecting Global Misuse

This section discusses how to mine global misuse instances of given patterns from log records of observed web transactions using Stream-SQL [12] and a StreamBase platform [11]. StreamSQL is an event pattern language that can be used to define queries over streams of data. StreamBase is an event processing platform that can run StreamSQL queries over input source file or database and produce outputs.

StreamSQL has several commands. CREATE INPUT STREAM creates data streams from a named file that is pre-configured according to a known schema. CREATE OUTPUT STREAM creates an output stream that is pre-configured according to a schema. A PATTERN phrase defines the search criteria from multiple input streams. A WITHIN phrase creates the maximum size of a window that moves along a collection of aligned streams when searching for a pattern.

Algorithm 1, which is based on StreamSQL, discovers Pattern I defined in Table 4. The algorithm accepts MEI records in ascending order of timestamps. The algorithm processes the pattern by filtering the records into two groups, invest and pay, using the predicates defined in Lines 6 and 7. This enables the pattern to employ the appropriate template (THEN phrase) in Line 9, i.e., invest messages are expected before pay messages. The predicates defined in Line 11 stipulate that the receiver of the invest message should be the sender of the following pay message; and the promoter value in the content of the invest message should be the receiver of the following pay message. The window size is set to 3 in Line 10. The SELECT phrase gathers the re- Detecting a few Ponzi-like recruits may not be sufficient to declare that a Ponzi scheme exists. To increase the confidence, a minimum support value is defined as a threshold; and an alert is sent only when the threshold is exceeded. Algorithm 2 can be employed to strengthen Algorithm 1 by incorporating a predefined minimum support value.

Algorithm 2 sends an alert when at least six Ponzi-like recruits are detected in the output of Algorithm 1. Note that the addition of Algorithm 2 to Algorithm 1 decreases the number of false positives.

Generating Comprehensive Evidence

This section discusses how to detect the orchestrator or the earliest known recruiter of a Ponzi scheme by tracing a potential recruit path to its root. Next, all possible paths that originate at the detected recruiter are examined to identify others who have invested in the Ponzi scheme. We begin by defining a choreography.

Suppose that the complaint includes an invest message msg and that the investment company I is used by the participants in the scheme. We define an ancestorChain(n) as:

(i) ancestorChain(0) = msg (ii) ancestorChain(n + 1) = [msg l;p (msg As a special case, we show how to compute the ancestor chain corresponding to Pattern II in Table 4 using StreamSQL in Algorithm 3. Given a recruiter, the algorithm traces ancestor recruiters to find the a window of size 6, limiting it to finding patterns within the specified period. Queries with narrower windows may cause the algorithm to miss more correlations than queries with wider windows.

Experimental Validation

The three-layered evidence generation framework provides evidence of web service misuse. The services in the bottom and middle layers have been studied by others [START_REF] Bilal | Fair BPEL processes transaction using non-repudiation protocols[END_REF][START_REF] Keller | The WSLA framework: Specifying and monitoring service level agreements for web services[END_REF][START_REF] Sahai | Automated SLA monitoring for web services[END_REF]. Therefore, we validate the top layer of the framework for which have introduced novel queries. 5 summarizes the queries and their associated patterns along with the names we have used throughout this paper. In order to test their accuracy and performance rates, we generated synthetic data and used a special simulation platform described below.

MEI-I MEI-X MEI-XX

Data Characteristics

To our knowledge, our illegal business scheme is novel and we were unable to obtain real data featuring the scheme. Consequently, we generated synthetic data in the MEI format called MEI-I that conforms with the Ponzi/pyramid schemes described above. MEI-I contains a total of 193,827 records (one record/second) over a period of three days (January 19, 2006 to January 21, 2006). The MEI-I file has a comma separated values (CSV) format. We focused on the capacity rate of data that embodies the illegal business scheme in its records.

To evaluate the performance, we created three more data sets using our seed data. The MEI-X data set contains the same Ponzi malicious activity, but is ten times larger than MEI-I in terms of the number of records. MEI-XX and MEI-L are two other sets that are 20 and 50 times larger than MEI-I, respectively (Figure 3).

Test Environment

As mentioned above, we defined use/misuse patterns using Stream-SQL and employed the StreamBase platform to detect the patterns. We used the specially-generated data described above in the MEI structure. The feed simulation platform of StreamBase was adjusted to accept this data in the form of CSV files. The platform empowers users to run StreamSQLs over any user-defined file satisfying the data schema ex-

Test Results

The performance of the algorithms was tested by executing three major queries (Table 5) over data sets of different sizes (Figure 4). As mentioned earlier, we observed that maximum proximity rates give the best decisions in the tests. Therefore, we built our performance tests using these rates in the algorithms.

Figure 8 presents the test results. Note that queries for detecting recruits and generating recruit trees take more time than queries for locating the orchestrator. The algorithm execution times for detecting Ponzi-like recruits are greater than the execution times for the other algorithms over the largest data set MEI-L. The recruit path climbing algorithm, which traverses records in the backward direction, exhibits the best performance and there is no need to use a larger window size.

We used the StreamBase Studio Feed Simulation platform in our validation tests. Other higher performance platforms are available, but we obtained reasonable performance despite using an integrated development environment instead of an enterprise environment. We also observed that the queries yield results with reasonable accuracy for the synthetic data sets. These results were obtained because we determined a reasonable window size for window-based queries for each data set. Also, we converged the evidence outcomes by tuning the time, property and key-based patterns associated with the queries.

MEI-I MEI-X MEI-XX

Conclusions

Web service choreographies can be used as the foundation for detecting illegal business schemes. Our implemented system engages Stream-SQL to define the associated use/misuse patterns and the StreamBase platform to detect the patterns in large streams of data. Although our choreographies only specify Ponzi-like schemes, the method can be used to specify (and detect) other illegal business activities [START_REF] Westphal | Data Mining for Intelligence, Fraud and Criminal Detection: Advanced Analytics and Information Sharing Technologies[END_REF] that can be mined from financial transaction repositories. Our future work will attempt to develop an online warning system that detects business schemes that appear legal from a microscopic view, but are illegal from a macroscopic perspective.

Figure 2 .

 2 Figure 2. Ponzi-like recruits over web services.

  Finally, φ denotes the empty string in k <w . Definition 3 (Recruit trees of fanout k and depth m): Suppose I is an investment company web service. Inductively define active(n) for every integer n as follows: (i) Active(0)={m} where m is defined as a message where [m.sender=p o , m.receiver=I]. (ii) Suppose active(n) has been defined and p ∈ k <w with |p|= n, then for each i ∈{0, . . . , k -1} define [active(pˆi)=msg ;p PayBack(pˆi)] where msg satisfies [msg.sender=P pˆi , msg.reciever=I and msg.content= invest] and PayBack(pˆi) satisfies PayBack(pˆi)=msg 1 ∩ p . . . ∩ p msg m where every msg i is of the form [msg i .content=pay, msg i .sender=I, msg i .receiver=ancestor(p,i)] for i = l. (iii) Let active(n + 1)=active(pˆ0) ∪ p . . . ∪ p active(pˆ(n -1)).

Algorithm 2

 2 Enhancing Ponzi Detection INPUT: PonziDetectOut from detectRecruits OUTPUT: Ponzi alerts DESCRIPTION: Counts detected Ponzi-like recruits using a window size of 6 as the minimum support. Emits Ponzi alerts when the minimum support is reached 1 SELECT "Ponzi Alerts," count() AS minSup 2 FROM PonziDetectOut [SIZE 6 TUPLES] 3 INTO PonziAlerts ; quired information about the detected pattern and emits the result to the PonziDetectOut table.

Figure 3 .

 3 Figure 3. Capacity values of test data.

Table 2 .

 2 Sample MEI table.

	ID	Time Sender Receiver Msg Content
	63..	21	A	B	r	<..invID..>
	67..	22	B	C	m	".."
	68..	23	C	B	k	<..payID..>
	and content. The first reference in the table is to a message sent by A
	to B with content <invID>.

Table 3 .

 3 MEI tuples featuring a misuse scheme.

	ID Time	Sender	Receiver	Msg	Content
	. . .
	45	B	InvComp invest Promoter=A
	55	InvComp	A	pay	150
	. . .
	67	C	InvComp invest Promoter=B
	76	InvComp	B	pay	150
	78	InvComp	A	pay	30
	. . .
	87	Victim	InvComp invest Promoter=C
	89	InvComp	C	pay	150
	92	InvComp	B	pay	30
	104	InvComp	A	pay	30

Table 4 .

 4 Ponzi scheme with fanout 1 and depth 1.

	I	invest;p pay where
		invest.sender=A and invest.reciever=B and pay.sender=B and
		pay.reciever=C and invest.prometer=pay.reciever
	II invest1;p invest2 where
		invest1.sender=C and invest1.receiver=B and invest2.sender=A
		and invest2.reciever=B and invest2.promoter= invest1.sender

  1;p k 1 ) ∪ p . . . ∪ p (msg l-1;p k l-1 ) ;p ancestor(ancestorChain(n),n)] satisfying the equations E: msg 1 .content=pay and k 1 .content=invest, . . . , msg l-1 .content=pay and k l-1 .content=invest (msg l .time < msg l .time < k 1 .time), . . . , (msg 1 .time < msg l-1 .time > k l-1 .time) and msg l .reciever=I. Also, we define Earliest(msg) as LFP(f, msg, E).

Table 5 .

 5 Query-pattern mapping.

	Query Name	Pattern	Pattern Name
	DetectRecruits	invest;p pay	invest-pay
	ClimbRecruitPath	invest1;p invest2 invest-invest
	GenerateRecruitTree invest1;p invest2 invest-invest
	Table

Table 6 .

 6 Test environmentpected by the algorithm. The platform also provides observable outputs of algorithm execution. Using the StreamBase Manager, it is also possible to observe CPU and memory usage during algorithm execution. Although StreamBase encourages the use of enterprise servers for benchmarking and improved performance, we observed that the feed simulation platform was adequate to test our queries over vast amounts of data with a reasonable resource allocation rate. Table 6 lists the hardware and software components of the test environment.

	Hardware	Software
	CPU: Intel Core2 T7400 2.16 GHz, OS: Windows XP SP2
	4 MB L2 Cache, 667 MHz FSB	JVM : SUN JDK 1.5.0.15
	Physical Memory: 2 GB, 995 MHz	StreamBase Studio: Version 6.4
	Hard Disk: 250 GB, 7200 rpm	Max Heap Size: 1024 MB

orchestrator. It identifies when the scheme began by traversing records in descending order of timestamps looking for the senders of invest messages. This is done by declaring the dynamic variable pointerPromoter in Line 4 to which the suspected promoter is passed by default (see DE-FAULT) as the orchestrator. Each time the output emits a hop (that meets the criterion in Line 7: invest message sent to the orchestrator) by the node to the receiver, the promoter value is extracted from the content of the invest message using a XPATH function (SELECT phrase in Line 9). It is then written to the output stream (Line 2) and assigned to the dynamic variable pointerPromoter in Line 5. The newly-assigned value is used as the predicate in Line 7 for locating the next message if it matches the sender value.

After the orchestrator or earliest recruits have been identified, a traceforward algorithm is used to generate the evidence. Since the algorithm above locates one of the oldest message records for the trace-forward algorithm, it helps gather the most comprehensive evidence of the scheme.

Algorithm 4 compiles the evidence using Pattern II in


Xavier Antoine
email: xavier.antoine@univ-lorraine.fr

Qinglin Tang
email: qinglintang@scu.edu.cn

Jiwei Zhang
email: jwzhang@csrc.ac.cn

Contact Qinglin

On the numerical solution and dynamical laws of nonlinear fractional Schrödinger/Gross-Pitaevskii equations

Keywords: fractional Schrödinger equation, fractional Gross-Pitaevskii equations, dynamical laws, conservation properties, discretization schemes, accuracy order, fast implementation

The purpose of this paper is to discuss some recent developments concerning the numerical simulation of space and time fractional Schrödinger and Gross-Pitaevskii equations. In particular, we address some questions related to the discretization of the models (order of accuracy and fast implementation) and clarify some of their dynamical properties. Some numerical simulations illustrate these points.

Introduction

The development of fractional partial differential equations (PDEs) has grown impressively during the last few years, because of the huge potential of emerging applications in science. In particular, some important impacts concern fractional quantum dynamics based on the space or/and time Fractional Schrödinger equation and its nonlinear version (SFNLSE or TFNLSE) [START_REF] Achar | Time fractional Schrödinger equation revisited[END_REF][START_REF] Bao | Numerical methods for computing ground state and dynamics of nonlinear relativistic Hartree equation for boson stars[END_REF][START_REF] Bhatti | Fractional Schrödinger wave equation and fractional uncertainty principle[END_REF][START_REF] Carusotto | Quantum fluids of light[END_REF][START_REF] Dong | Space-time fractional Schrödinger equation with time-independent potentials[END_REF][START_REF] Elgart | Mean field dynamics of boson stars[END_REF][START_REF] Fröhlich | Effective dynamics for boson stars[END_REF][START_REF] Lomin | Fractional-time quantum dynamics[END_REF][START_REF] Naber | Time fractional Schrödinger equation[END_REF][START_REF] Pinsker | Fractional quantum mechanics in polariton condensates with velocity-dependent mass[END_REF][START_REF] Wang | Generalized fractional Schrödinger equation with space-time fractional derivatives[END_REF]. The fractional nonlinear Schrödinger equation is used to describe the nonlocal phenomena in quantum physics and to explore the quantum behaviors of either long-range interactions or time-dependent processes with many scales [1, 45, 46, 48-50, 58, 64, 74, 81, 82]. For example, FNLSE arise in the modeling of quantum fluids of light [START_REF] Carusotto | Quantum fluids of light[END_REF], boson stars [START_REF] Bao | Numerical methods for computing ground state and dynamics of nonlinear relativistic Hartree equation for boson stars[END_REF][START_REF] Elgart | Mean field dynamics of boson stars[END_REF][START_REF] Fröhlich | Effective dynamics for boson stars[END_REF] and polariton condensates [START_REF] Pinsker | Fractional quantum mechanics in polariton condensates with velocity-dependent mass[END_REF]. Some analytical and approximate solutions have been considered e.g. for the TFSE [START_REF] Khan | Approximate solutions to time-fractional Schrödinger equation via homotopy analysis method[END_REF][START_REF] Odibat | Analytic study on time-fractional Schrödinger equations: exact solutions by GDTM[END_REF]. The complexity of FNLSE with different potentials and nonlinearities requires the development of efficient and accurate numerical simulations [START_REF] Bhrawy | A fully spectral collocation approximation for multidimensional fractional Schrödinger equations[END_REF][START_REF] Duo | Mass conservative method for solving the fractional nonlinear Schrödinger equation[END_REF][START_REF] Garrappa | Solving the time-fractional Schrödinger equation by Krylov projection methods[END_REF][START_REF] Hicdurmaz | A stable numerical method for multidimensional time fractional Schrödinger equation[END_REF][START_REF] Khaliq | A fourth-order implicit-explicit scheme for the space fractional nonlinear Schrödinger equations[END_REF][START_REF] Kirkpatrick | Fractional Schrödinger dynamics and decoherence[END_REF][START_REF] Li | The numerical computation of the time-fractional linear Schrödinger equation on unbounded domain[END_REF][START_REF] Liang | The locally extrapolated exponential splitting scheme for multi-dimensional nonlinear space-fractional Schrödinger equations[END_REF][START_REF] Mohebbi | The use of a meshless technique based on collocation and radial basis functions for solving the time fractional nonlinear Schrödinger equation arising in quantum mechanics[END_REF][START_REF] Wei | Analysis of an implicit fully discrete local discontinuous Galerkin method for the time-fractional Schrödinger equation[END_REF] to go further.

The aim of the paper is not to review and compare all the numerical methods that have been designed during the last years for the SFNLSE and TFNLSE but rather to discuss some developments of efficient, accurate, stable and physically relevant numerical schemes through examples, most particularly regarding the dynamical laws for FNLSE where many questions remain open. For the integer order NLSE, the picture is now pretty clear [START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF][START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF] and the development of high-order schemes follow a welldefined route, in particular concerning the discrete dynamical laws that a numerical scheme should mimic regarding the continuous physical equations. Concerning the Gross-Pitaevskii equation (GPE) arising in Bose-Einstein Condensation (BEC) [START_REF] Pitaevskii | Bose-Einstein Condensation[END_REF] as a generalization of the NLSE, numerous extensions to relevant physical situations are now clarified [START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF][START_REF] Antoine | Modeling and Computation of Bose-Einstein Condensates: Stationary States, Nucleation, Dynamics, Stochasticity[END_REF][START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF] (multi-components, nonlocal nonlinear interactions,...). For the fractional case, the situation is more complicated and still needs to be analyzed deeply. The aim of this paper is to contribute to the topic, by discussing the question of the dynamical laws that the space or time fractional NLSE/GPE should fulfill to yield correct physical solutions, while being efficient, accurate and stable. The schemes used here have mainly been designed by the authors and serve as examples to illustrate the purpose. Finally, the paper can be seen as a complement of [START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF] while being much more prospective and does not pretend to be a review paper.

The plan of the paper is the following. In Section 2, we recall the main results concerning the features of the most popular schemes for the standard NLSE/GPE. In Section 3, we develop the space fractional NLSE/GPE, some of its properties, the numerical schemes leading to a fast implementation with coherent physical features, and finally provide some numerical examples. In Section 4, we consider the TFNLSE and discuss the dynamical laws. This case is the most difficult to analyze since only a few theoretical results are available. Then, we explain some fast evaluation schemes for the TFNLSE and report some illustrative examples to address the difficult issue of the dynamical laws that can be expected at the discrete level, according to the underlying fractional equation. We finally end the paper with a conclusion and draw a few perspectives in Section 5.

The standard NLSE/GPE

Integer-order models and dynamical properties

We start by considering the following time-dependent cubic NLSE [START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF][START_REF] Sulem | The Nonlinear Schrödinger Equation: Self-Focusing and Wave Collapse[END_REF] i∂ t ψ(t, x) = -1 2 ∆ + V (x) + β|ψ(t, x)| 2 ψ(t, x), x ∈ R d , t > 0, (2.1)

ψ(t = 0, x) = ψ 0 (x), x ∈ R d . (2.2)
In the above equations, i = √ -1 is the complex unit, t is the time variable, x ∈ R d is the spatial variable (for d = 1, 2, 3), ψ := ψ(t, x) is the unknown complexvalued wave function, ∆ = ∇ 2 is the usual Laplace operator in dimension d (and ∇ the gradient operator) and ψ 0 := ψ 0 (x) is a given complex-valued initial data. The real-valued external potential function V := V (x) is given and its definition is related to the underlying application. For instance, for BEC, it can be chosen as a harmonic confining trap, an optical lattice potential or even a quadratic-plus-quartic function [START_REF] Huang | A Bloch decomposition based numerical method for quantum dynamics in periodic potentials[END_REF][START_REF] Huang | Numerical simulation of the nonlinear Schrödinger equation with multidimensional periodic potentials[END_REF][START_REF] Pitaevskii | Bose-Einstein Condensation[END_REF]; for nonlinear optics applications, it could be set as an attractive potential. The potential function V may also depend on the time variable t or can even be stochastic. The standard cubic nonlinearity involves the density function ρ := |ψ| 2 := ψ ψ, the parameter β being the nonlinearity strength [START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF][START_REF] Pitaevskii | Bose-Einstein Condensation[END_REF] (positive for a repulsive/defocusing interaction and negative for an attractive/focusing interaction) and ψ the complex conjugate of the wave function ψ. Other local nonlinearities are used in nonlinear optics, i.e. cubic-quintic nonlinearity or saturation of the intensity nonlinearity [START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF][START_REF] Bao | Numerical methods and comparison for computing dark and bright solitons in the nonlinear Schrödinger equation[END_REF], or even nonlocal nonlinearities like for dipolar gazes [START_REF] Bao | Dimension reduction of the Schrödinger equation with Coulomb and anisotropic confining potentials[END_REF]. The cubic NLSE (2.1) is also called the Gross-Pitaevskii Equation (GPE) in the framework of BEC [START_REF] Antoine | Modeling and Computation of Bose-Einstein Condensates: Stationary States, Nucleation, Dynamics, Stochasticity[END_REF][START_REF] Pitaevskii | Bose-Einstein Condensation[END_REF].

Important dynamical properties are verified for the solution ψ of (2.1). Here we mention the most important ones that are also hopefully expected to be fulfilled at the discrete level to get some robust numerical methods. The integer-order NLSE (2.1) is a dispersive PDE that is time reversible, i.e. it is unchanged under the change of time variable t → -t, and next taking the conjugate into the equation. A second crucial property is gauge invariance, that is, if V → V + C (C ∈ R), then the solution ψ → ψe -iCt implies that the density ρ = |ψ| 2 is unchanged. The NLSE (2.1) conserves some quantities over the time, such as the mass and energy [START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF][START_REF] Cazenave | Semilinear Schrödinger Equations[END_REF][START_REF] Sulem | The Nonlinear Schrödinger Equation: Self-Focusing and Wave Collapse[END_REF], for t ≥ 0,

N (t) := ψ(t, •) 2 2 = R d |ψ(t, x)| 2 dx ≡ R d |ψ(0, x)| 2 dx := N (0), (2.3)
E 2 (t) := R d 1 2 |∇ψ(t, x)| 2 + V (x)|ψ(t, x)| 2 + β 2 |ψ(t, x)| 4 dx ≡ E 2 (0). (2.4)
In the free-potentiel situation, i.e. V (x) ≡ 0, the momentum and angular momentum are conserved [START_REF] Sulem | The Nonlinear Schrödinger Equation: Self-Focusing and Wave Collapse[END_REF]. The NLSE (2.1) admits the plane wave solution as

ψ(t, x) = Ae i(k•x-ωt)
, where the time frequency ω, amplitude A and spatial wave number k satisfy the dispersion relation

ω = |k| 2 2 + β|A| 2 . (2.5)
More dynamical properties of the NLSE (2.1) such as the well-posedness and finite time blow-up and the properties of its solitary solutions in 1d can be found in [START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF][START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF][START_REF] Bao | Numerical methods and comparison for computing dark and bright solitons in the nonlinear Schrödinger equation[END_REF][START_REF] Cazenave | Semilinear Schrödinger Equations[END_REF][START_REF] Sulem | The Nonlinear Schrödinger Equation: Self-Focusing and Wave Collapse[END_REF] and references therein.

Overview of popular numerical schemes

To numerically simulate the dynamical system (2.1), the full-space problem usually is truncated suitably so that the perturbation of the solution at the boundary is negligible. This results therefore in an initial boundary-value problem. In this paper, we restrict ourselves to consider the case where the solution is confined within domain

D :=]a; b[ d .
This is a standard situation, particularly for BEC, where the potential is strongly confining. Therefore, the choice of the boundary condition (BC) does not essentially modify the solution if D is large enough. The specific BCs are chosen according to the spatial discretization scheme one uses. More precisely, when a second-order finite-difference (FD) scheme is applied (see subsection 4.2), we impose a homogeneous Dirichlet BC. For the Fourier pseudospectral (SP) scheme (see subsection 3.2), we will use a periodic BC. In addition, if the solution ψ is not confined in D, i.e., ψ can strike the boundary, then much more complicated BCs are required, such as the transparent, artificial and absorbing BCs as well as perfectly matched layers (see references [START_REF] Antoine | A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schrödinger equations[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF][START_REF] Zheng | A perfectly matched layer approach to the nonlinear Schrödinger wave equation[END_REF] for the integer order Schrödinger equations and [START_REF] Zhang | Efficient numerical computation of time-fractional nonlinear Schrödinger equations in unbounded domain[END_REF] for fractional order problems). However, this is out of the scope of the current paper.

Considering now the time approximation schemes for the Schrödinger equation, it is well-known that building schemes that preserve at the discrete level the physical properties discussed above in subsection 2.1 is nontrivial. We do not develop here a full discussion and rather choose to report in Table 1 the conclusion from [START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF] for four specific schemes (where the discretization of the physical quantities must be suitably defined thanks to the approximation scheme). TSSP (respectively TSFD) refers to the second-order Strang Time Splitting scheme with pseudospectral SP [START_REF] Bao | On time-splitting spectral approximation for the Schrödinger equation in the semiclassical regime[END_REF][START_REF] Bao | Numerical study of time-splitting spectral discretization of nonlinear Schödinger equations in the semi-classical regimes[END_REF][START_REF] Besse | Order Estimates in Time of Splitting Methods for the Nonlinear Schrödinger Equation[END_REF][START_REF] Lubich | On splitting methods for Schrödinger-Poisson and cubic nonlinear Schrödinger equations[END_REF] (respectively second-order FD) approximation in space. CNFD and ReFD are related to the second-order Crank-Nicolson (CN) and Relaxation schemes [START_REF] Besse | A relaxation scheme for nonlinear Schrödinger equation[END_REF] with FD.

Let us now discuss the convergence properties of the various schemes as well as some of the implementation issues, taking the one-dimensional case to explain the technical details. We assume here that we discretize in space with J -1 equally spaced interior points, i.e. x j = jh x , 0 ≤ j ≤ J (with x 0 = a and x J = b) in ]a; b[ (the extension to the d-dimensional situation is direct.) The spatial mesh size is then h = h x = (b -a)/J. In practice, we wish to compute the numerical solution from t = 0 to a maximal computational time t = T . Let us introduce the N + 1 uniformly distributed discrete times (t n ) 0≤n≤N such that: t n = n∆t, for n = 0, ..., N , with ∆t := T /N . The four schemes are unconditionally stable, i.e. there is no Courant-Friedrichs-Lewy (CFL) condition involving h and ∆t. In addition, the schemes are all second-order in time, and in space second-order for FD or spectral for SP.

The only fully explicit scheme is TSSP, where the nonlinearity and potential terms can be directly integrated in time, leading then to diagonalize the remaining freespace linear Schrödinger operator through the FFT (see also subsection 3.2) at a computational cost O(J d log J). In terms of memory storage, we need the values of the solution (ψ n j ) 0≤j≤J (with ψ n j ≈ ψ(x j , t n ), and ψ n ≈ ψ(•, t n )) only at time t n since no time memory effect is involved into the equation (indeed, i∂ t is a local operator). TSFD and ReFD are linearly implicit schemes since the nonlinearity is made explicit. Hence, each time step requires to store and solve the associated linear system, i.e. corresponding to O(J d ) coefficients to store and to O(J d log J) operations with an adapted fast linear algebra solver for structured systems. Of course, the solution must also be stored at time t n . Finally, since the CNFD scheme is fully implicit, each time step requires the solution to a nonlinear system that solved by iterative methods. As a consequence, the computational cost is clearly much larger than O(J d ), while the memory storage remains O(J d ). Table 1

(J d ) O(J d ) O(J d ) O(J d ) Computational Cost O(J d log J) O(J d ) 2 O(J d log J) 3 O(J d log J) 4
Table 1. Integer order NLSE : physical/numerical properties of various numerical methods in the ddimensional case.

Let us remark that high-order TSSP schemes [START_REF] Bao | A fourth-order time-splitting Laguerre-Hermite pseudospectral method for Bose-Einstein condensates[END_REF][START_REF] Shen | Spectral Methods. Algorithms, Analysis and Applications[END_REF][START_REF] Thalhammer | Convergence analysis of high-order time-splitting pseudo-spectral methods for nonlinear Schrödinger equations[END_REF] can also be built for the NLSE/GPE. Nevertheless, these approaches do not always apply for some equations, e.g. when the system is non-autonomous. Recent high-order schemes (exponential integrators [START_REF] Besse | High-order exponential integrators for nonlinear Schrödinger equations with application to rotating Bose-Einstein condensates[END_REF][START_REF] Thalhammer | High-order exponential operator splitting methods for time-dependent Schrödinger equations[END_REF][START_REF] Thalhammer | A numerical study of adaptive space and time discretisations for Gross-Pitaevskii equations[END_REF], IMEXSP [START_REF] Antoine | High-Order IMEX-spectral schemes for computing the dynamics of systems of nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF]...) with time-stepping techniques have been designed to get new efficient solvers that could be combined with pseudospectral approximation schemes. We also note that ReFD can be extended to include the pseudospectral approximation, resulting in the ReSP scheme given in [START_REF] Antoine | Perfectly matched layers and pseudospectral schemes for the dynamics of Gross-Pitaevskii equations with rotating term[END_REF]. When developing the ReSP method, each time step t n requires the solution to a n-dependent linear system that can be easily solved by combining an iterative Krylov subspace solver (GMRES) and an adapted preconditioner. From the analysis above, we see that TSSP provides a suitable way to solve the NLSE.

The space FNLSE/FGPE

Space fractional models and dynamical properties

We consider now the specific Space Fractional NLSE (SFNLSE )

i∂ t ψ(t, x) = 1 2 (-∆) s + V (x) + β|ψ(t, x)| 2 ψ(t, x). (3.6)
This equation is called SFGPE in the framework of the GPE for BECs. The real-valued parameter s > 0 is the space fractional order defining the nonlocal dispersive interaction. The fractional dispersion is called superdispersion (respectively subdispersion) for s > 1 (respectively s < 1) [START_REF] Antoine | On the ground states and dynamics of space fractional nonlinear Schrödinger/Gross-Pitaevskii equations with rotation term and nonlocal nonlinear interactions[END_REF]. The fractional kinetic operator is defined via a Fourier integral operator

(-∆) s ψ = 1 (2π) d R d ψ(k) |k| 2s e ik•x dk, (3.7)
where the Fourier transform is

ψ(k) = R d ψ(x)e -ik•x dx.
For s = 1, the SFNLSE (3.6) simplifies to the standard cubic NLSE (2.1). A more general form of (3.6) can be found in [START_REF] Antoine | On the ground states and dynamics of space fractional nonlinear Schrödinger/Gross-Pitaevskii equations with rotation term and nonlocal nonlinear interactions[END_REF] where a nonlocal nonlinear interaction term given by a convolution kernel U (corresponding to a Coulomb or dipolar interaction) is added to the cubic nonlinearity, the Laplace operator involves a mass term m and finally a rotation term -ΩL z ψ (where L z = -i(x∂ y -y∂ x ) is the z-component of the angular momentum, Ω representing the rotating frequency) is included into the equation. More precisely, the system writes (with λ ∈ R)

i∂ t ψ(t, x) = 1 2 -∆ + m 2 s + V (x) + β|ψ| 2 + λΦ(t, x) -ΩL z ψ(t, x), (3.8) Φ(t, x) = U * |ψ(t, x)| 2 , x ∈ R d , t > 0, d ≥ 2. (3.9)
For (3.6) and (3.8)-(3.9), the mass conservation N given by (2.3) still holds. We can also define the conservation of a fractional energy E s as

E s (t) = R d 1 2 ψ -∇ 2 +m 2 s ψ+V |ψ| 2 + β 2 |ψ| 4 + λ 2 Φ|ψ| 2 -Ω ψL z ψ dx ≡ E s (0). (3.10)
We remark that m = 0 and the terms λ and Ω in (3.10) shall be omitted if d = 1. The equation is clearly still time reversible and gauge invariant. In addition, the dispersion relation (for V = 0) now reads as

ω = |k| 2s 2 + β|A| 2 . (3.11)
In the case of (3.8)-(3.9), extensions of some standard dynamical laws (mainly related to the angular momentum expectation and center of mass) have been stated in [START_REF] Antoine | On the ground states and dynamics of space fractional nonlinear Schrödinger/Gross-Pitaevskii equations with rotation term and nonlocal nonlinear interactions[END_REF] for some values of s. Nevertheless, some questions still remain open such as the dynamical laws of the center of mass (for s > 1) and condensate widths, the well-posedness of the SFNLSE with general s and attractive interaction (i.e. β < 0).

Numerical schemes and their efficient implementation

From the conclusion of subsection 2.2, a natural scheme for solving (3.6) is to use an adapted TSSP scheme since the space fractional operator (-∆) s is naturally represented through an inverse Fourier transform as a pseudodifferential operator. Since the physics is confined within a finite computational domain, then periodic BCs can be set on the boundary. From t = t n to t = t n+1 := t n + ∆t and for the one-dimensional case (d = 1, the extension being direct to d ≥ 2), we solve the system in two steps. First, one considers

i∂ t ψ(t, x) = V (x) + β|ψ| 2 ψ(t, x), x ∈]a; b[, t n ≤ t ≤ t n+1 , (3.12)
for a time step ∆t, and next we solve

i∂ t ψ(t, x) = 1 2 (-∂ xx ) s ψ(t, x), x ∈]a; b[, t n ≤ t ≤ t n+1 , (3.13)
with periodic BCs at {a; b} for the same time step. The linear subproblem (3.13) is discretized in space by the Fourier pseudo-spectral method and integrated in time exactly in the Fourier space. Similarly to s = 1, the nonlinear subproblem (3.12) preserves the density, i.e. |ψ(t,

x)| 2 ≡ |ψ(t = t n , x)| 2 = |ψ n (x)| 2 , leading to ψ(t, x) = e -i[V (x)+β|ψ n (x)| 2 ](t-tn) ψ n (x), x ∈]a; b[, t n ≤ t ≤ t n+1 . (3.14)
Let us introduce J as an even positive integer and µ = 2π /(b -a), for -J 2 ≤ ≤ J 2 -1. We denote by ψ n the numerical solution at time t = t n , with components ψ n j , 0 ≤ j ≤ J, for the initial data ψ 0 = (ψ 0 (x j )) 0≤j≤J . A second-order TSSP scheme for solving (3.6) is

ψ (1) j = J/2-1 =-J/2 e -i∆t 4 |µ | 2s (ψ n ) e iµ (xj-a) , (3.15)
ψ (2) j = ψ (1) j e -i∆t[V (xj)+β|ψ (1) j | 2 ] , (3.16
)

ψ n+1 j = J/2-1 =-J/2 e -i∆t 4 |µ | 2s (ψ (2) ) e iµ (xj-a) , (3.17)
for 0 ≤ j ≤ J and n ≥ 0. Here, (ψ n ) and (ψ (2) ) are the discrete Fourier series coefficients of ψ n and ψ (2) , respectively. The TSSP scheme (3.15)-(3.17) for s inherits of the properties for s = 1 given in Table 1 and can be extended to high-order timesplitting schemes.

The extension to the general system (3.8)-(3.9) is developed in [START_REF] Antoine | On the ground states and dynamics of space fractional nonlinear Schrödinger/Gross-Pitaevskii equations with rotation term and nonlocal nonlinear interactions[END_REF]. It essentially needs to take care about the efficient spectral approximation of the nonlocal nonlinear interaction term (based on a Gaussian-Sum evaluation) and the elimination of the rotating term through a rotating Lagrangian coordinates transformation. This latter modification results in the appearance of a new time-dependent potential which must be integrated numerically in the TSSP scheme. Let us finally remark that ReSP could also be extended to the SFNLSE, and that, more generally, any well-designed highorder time integrator (TSSP, IMEXSP...) could be applied here.

Numerical examples

In this section, we consider the effect of the fractional order s on the dynamics of the SFNLSE for various situations. We only show the results in 1d and 2d. For d = 1, 2, the potential, computational domain and time step are chosen as

V (x) = 1 2 γ 2 x x 2 , d = 1, γ 2 x x 2 + γ 2 y y 2 , d = 2, D = - 64 d , 64
d d , ∆t = 10 d-5 . (3.18)
Moreover, we take the mesh size h and initial data ψ 0 as

h = h x = 1 64 , d = 1, h = h x = h y = 1 16 , d = 2, ψ 0 (x) = φ s0 gs (x -x 0 ), (3.19)
where φ s0 gs is the ground states (gs) of the SFNLSE (3.6) (if d = 1) and/or (3.8) (with m = λ = 0, if d = 2) with potential given by V (x) (3.18) and fractional order s 0 . The ground states can be computed either by the popular gradient flow method [START_REF] Antoine | On the ground states and dynamics of space fractional nonlinear Schrödinger/Gross-Pitaevskii equations with rotation term and nonlocal nonlinear interactions[END_REF][START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF][START_REF] Bao | Numerical methods for computing ground state and dynamics of nonlinear relativistic Hartree equation for boson stars[END_REF][START_REF] Bao | Accurate and efficient numerical methods for computing ground states and dynamics of dipolar Bose-Einstein condensates via the nonuniform FFT[END_REF] or the recently developed conjugated gradient flow approach [START_REF] Antoine | Efficient spectral computation of the stationary states of rotating Bose-Einstein condensates by the preconditioned nonlinear conjugate gradient method[END_REF]. Here, we apply the method proposed in [START_REF] Antoine | On the ground states and dynamics of space fractional nonlinear Schrödinger/Gross-Pitaevskii equations with rotation term and nonlocal nonlinear interactions[END_REF].

Example 3.1. Here, we consider a simple 1d case to show the difference of dynamics for various fractional powers s. For d = 1, we take γ x = 1 in (3.18) and β = 250 in (3.6). We consider two cases of initial data ψ 0 (x) (see Eq. (3.19)):

Case 1. Fix s 0 = 1 and x 0 = 0, i.e., always take the initial data as the ground states of classical NLSE, then only vary s. Case 2. Take s 0 = s, and x 0 = 8, i.e., we take the initial data as the ground states of the SFNLSE with initial center shifted to x 0 . Then, we consider the dynamics for different fractional orders s.

All the ground states are computed by using the preconditioned nonlinear conjugate gradient method developed in [START_REF] Antoine | Efficient spectral computation of the stationary states of rotating Bose-Einstein condensates by the preconditioned nonlinear conjugate gradient method[END_REF], eventually adapted to the SFNLSE. Figs. 1 and2 show the evolution of the density |ψ| 2 , the mass N (t) and the energy E s (t) for different fractional orders s in Cases 1 and 2. From these experiments and others not shown here for brevity, we could see that: 1) The mass and energy are conserved very well during the dynamics, for all the cases; 2) The fractional order s affects the dynamics of the density |ψ| 2 significantly and qualitatively (cf. Fig. 1). Oscillations usually arise during the dynamics, and would sometimes turn to a chaotic dynamics for long-time simulations; 3) When s is larger, one needs larger domain sizes to well-resolved the solution in phase space, while a smaller mesh size h is needed if s is smaller; 4) For a center-shifted initial profile, the larger s is, the faster the center of mass oscillates (cf. Fig. 2). Again, for a larger s, oscillation of the density occurs and chaotic dynamics arises (cf. Fig. 2 gs as Type 2 and set x 0 = (0, 0) T . We perturb the trapping frequency in the y-direction, i.e., we set γ y = 1.3. The other parameters are kept unchanged. Case 6. Same as Case 5, but now perturb the trapping frequency in both x and y-directions, i.e., we set γ x = γ y = 1.3.

Figs. 3 iii) and iv) show the evolution of the mass N (t) and the energy E s (t), while Fig. 4 shows the contour plots of the density |ψ(x, t)| 2 at different times for Case 1 to 6. From these figures and other experiments not shown here for the sake of brevity, we could see that: 1) The masses and total energies are conserved well for all cases;

2) The fractional order s affects the dynamics significantly and qualitatively. For the Type 1 initial data, when s < 1 (subdispersion), the vortex lattice is more condensed (i.e. the compact support of the density is smaller than the initial setup), and rotates with breather-like dynamics. On the other hand, when s > 1 (superdispersion), the vortex lattice is less condensed and similarly rotates with breather-like dynamics. The structure of the vortex lattice is destroyed although some symmetric properties of the density are still kept for long time dynamics (cf. Figs 4 (a)-(c)). Eventually, a chaotic dynamics and even a turbulent behaviour arises for most cases; 2) For Type 2 initial data, the ring-structure is destroyed if either the initial center of mass is shifted or the trapping potential are perturbed asymmetrically (cf. Figs 4 (d)-(e)). On the contrary, if only and symmetrically the trapping potential are perturbed, the ring-structure will be kept and the ring undergo a breather-like dynamics (cf. Fig. 4 (f)), similar as the dynamics in standard NLSE (i.e. s = 1). Based on all these examples, we conclude that the dynamics of a fractional BECs modeled by the SFNLSE can be deeply affected by the fractional power s, depending on the sub-or superdispersion situation.

The time FNLSE/FGPE

Time fractional models and dynamical properties

We consider now the general Time-Fractional NLSE (TFNLSE) given by

(γ R + iγ I )(i C 0 D t ) α ψ(t, x) = - 1 2 ∆ψ + V (x)ψ + β|ψ| 2 ψ, x ∈ R d , t > 0, (4.20)
with initial data ψ(0, x) = ψ 0 (x). The operator C 0 D α t denotes the Caputo fractional derivative of order α (0 < α < 1) with respect to t [START_REF] Podlubny | Fractional Differential Equations[END_REF] and defined by

C 0 D α t ψ(t, x) = 1 Γ(1 -α) t 0 1 (t -u) α ∂ψ(u, x) ∂u du, 0 < α < 1, (4.21)
where Γ(•) is the Gamma special function. To define a general setting, we introduce a complex-valued coefficient in front of the fractional derivative, i.e. (γ R + iγ I ), where γ R > 0 and γ I ≥ 0. The reason is to get a unified writing of the various models that can be met in the literature, like the following one introduced in [1]

i( C 0 D t ) α ψ(t, x) = - 1 2 ∆ψ + V (x)ψ + β|ψ| 2 ψ, (4.22)
by choosing γ R + iγ I = i 1-α , and seen as a modified version of Naber model's [START_REF] Naber | Time fractional Schrödinger equation[END_REF] (

i C 0 D t ) α ψ(t, x) = - 1 2 ∆ψ + V (x)ψ + β|ψ| 2 ψ, (4.23)
corresponding to γ R = 1 and γ I = 0. Here, we set ( C 0 D t ) α = C 0 D α t . Therefore, one can interpret (4.22) as an extension of (4.23) with a dissipative term (up to a multiplicative constant) since indeed γ R = cos((1 -α)π/2) > 0 and γ I = sin((1 -α)π/2) ≥ 0, for 0 < α < 1. This generalizes the standard dissipative form of the GPE [START_REF] Tsubota | Vortex lattice formation in a rotating Bose-Einstein condensate[END_REF] which is used in quantum turbulence.

In [START_REF] Dong | Space-time fractional Schrödinger equation with time-independent potentials[END_REF], the authors prove that there is no conservation law for the linear version of (4.23) (i.e. for β = 0). In particular, for 0 < α < 1, they obtain that lim t→+∞

N (t) = 1 α 2 > 1, (4.24)
which means that particles are created during the fractional dynamical process (extracted from the confined potential field). In addition, they show that the energy E 2 (t) goes to a limiting value when t tends towards +∞, and is given through averaging. Therefore, since (4.22) can be seen as a version of (4.23) with dissipative term, one can expect that inversely a certain amount of particles is absorbed when β = 0, leading then to a loss of mass and a decay of the energy. To the best of our knowledge, no general result is available for the nonlinear case where probably many possible situations can arise (see subsection 4.3 for some numerical illustrations depending on the models and assumptions).

Concerning the dispersion relation, if one uses a travelling plane wave ψ(t, x) = Ae i(k•x-ωt) , and since one gets

( C 0 D t ) α (Ae i(k•x-ωt) ) = Ae ik•x ( C 0 D t ) α (e -iωt )
and

( C 0 D t ) α (e -iωt ) := -iωt 1-α E 1,2-α (-iωt),
where E a,b is the two-parameter Mittag-Leffler special function [START_REF] Haubold | Mittag-Leffler functions and their applications[END_REF] with parameters (a, b) := (1, 2 -α), then the dispersion relation writes

-(γ R + iγ I )i 1+α ωt 1-α E 1,2-α (-iωt)e iωt = |k| 2 2 + β|A| 2 , (4.25)
which simplifies to the standard dispersion relation for α = 1. Equation (4.25) is timedependent and not easy to understand at first sight. Let us nevertheless remark that one gets the asymptotic formula for the Mittag-Leffler functions

E a,b (z) = 1 a z (1-b)/a e z 1/a - N * r=1 1 Γ(b -ar) 1 z r + O( 1 z N * +1 ), for |z| → ∞, | arg z| ≤ µ,
where π/2 < µ < π and N * ∈ N * and N * > 1. In our situation, a = 1, b = 2 -α and z := -iωt, which means that | arg z| = π/2. As a consequence, we obtain the following asymptotic dispersion formula

(γ R + iγ I )ω α = |k| 2 2 + β|A| 2 , (4.26)
which now clearly translates the time fractional behavior of the equation. With the substitution ψ → ψe -iCt and from the Leibniz rule for Caputo fractional derivatives [START_REF] Podlubny | Fractional Differential Equations[END_REF] given as an infinite series expansion, the gauge change property does not hold if α = 1. In addition, the equation is not a priori time reversible, in particular because of the dissipation term that may appear. From the discussion above, we see that the situation still needs to be clarified to understand which physical dynamical laws/properties a good discrete scheme should fulfill to be acceptable. This is probably one interesting point to investigate deeper in the future and that will be discussed in subsection 4.3 through numerical simulations. Before that, we introduce in subsection 4.2 some discretization schemes as well as their efficient implementation, most particularly regarding the computational cost, memory storage and order of accuracy.

Numerical schemes and their efficient implementation

As in the previous sections, we assume that we can restrict the computational domain to a finite box D :=]a; b[ with suitable boundary conditions (homogeneous Dirichlet/Neumann or periodic boundary conditions).

Before developing some discretization schemes for the full fractional system (4.20), we first address the problem of designing fast evaluation schemes for the fractional Caputo derivative. This operator depends on the history information and various discretization schemes can be developed, such as the L1-approximation [START_REF] Langlands | The accuracy and stability of an implicit solution method for the fractional diffusion equation[END_REF][START_REF] Li | Numerical analysis and physical simulations for the time fractional radial diffusion equation[END_REF][START_REF] Lin | Finite difference/spectral approximations for time-fractional diffusion equation[END_REF][START_REF] Sun | A fully discrete difference scheme for a diffusion-wave system[END_REF] given at time t n by

C 0 D α tn ψ n = ∆t -α Γ(2 -α) n k=1 a n-k (ψ k -ψ k-1 ), (4.27)
where

a k = (k + 1) 1-α -k 1-α (k ≥ 0). If ψ ∈ C 2 ([0, T ]
), the truncation error satisfies [START_REF] Lin | Finite difference/spectral approximations for time-fractional diffusion equation[END_REF][START_REF] Sun | A fully discrete difference scheme for a diffusion-wave system[END_REF]]

| C 0 D α tn ψ n -C 0 D α tn ψ n | ≤ C∆t 2-α . ( 4

.28)

A high-order scheme for the Caputo derivative is the L2-1 σ formula proposed in [START_REF] Alikhanov | A new difference scheme for the time fractional diffusion equation[END_REF].

The scheme is approximated at the time point t n+σ , with σ = 1 -α/2, by

H D α t ψ n+σ = 1 Γ(1 -α) n k=1 tk tk-1 (Π 2,k ψ(s)) (t n+σ -s) α ds + tn+σ tn (Π 1,n ψ(s)) (t n+σ -s) α ds , (4.29)
where Π 2,n ψ(t) and Π 1,n ψ(t) are the quadratic and linear interpolations given by

Π 2,n ψ(t) = ψ n-1 (t -t n )(t -t n+1 ) 2∆t 2 -ψ n (t -t n-1 )(t -t n+1 ) ∆t 2 + ψ n+1 (t -t n-1 )(t -t n ) 2∆t 2 , Π 1,n ψ(t) = ψ n t n+1 -t ∆t + ψ n+1 t -t n ∆t .
As shown in [START_REF] Alikhanov | A new difference scheme for the time fractional diffusion equation[END_REF] under the assumption that ψ ∈ C 3 ([0, T ]), the truncation error satisfies

C 0 D α t ψ(t n+σ ) = H D α t ψ n+σ + O(∆t 3-α ), n = 0, 1, • • • , N -1.
From the approximation (4.27), one can see that evaluating the Caputo fractional derivative requires the storage of all the discrete past values of the unknown wave function ψ, i.e. ψ 0 , ψ 1 , • • • , ψ n , and O(n) flops at the n-th time step and for each spatial grid point x j , where j = 0, • • • , J. Thus, in the one-dimensional case (i.e. d = 1), the average storage is O(N J) and the total computational cost is O(N 2 J), which is clearly much larger than for the standard case (see Table 1). As a consequence, this is a severe limitation for the long-time simulation of the TFNLSE, most particularly for higher dimensional problems since the number of grid points then grows as J d . To reduce both the storage and computational cost, fast evaluation algorithms of the Caputo derivative can be built [START_REF] Jiang | Fast evaluation of the Caputo fractional derivative and its applications to fractional diffusion equations[END_REF]. The main idea of the acceleration method is to split the Caputo derivative into the sum of a local time part and a history part

C 0 D α tn ψ = 1 Γ(1 -α) tn tn-1 ψ (u, x) (t n -u) α du + 1 Γ(1 -α) tn-1 0 ψ (u, x) (t n -u) α du := C loc (t n , x) + C hist (t n , x).
For the local part C loc (t n , x), we apply the L1-approximation, i.e.

C loc (t n , x) ≈ ψ(t n , x) -ψ(t n-1 , x) Γ(1 -α)∆t tn tn-1 1 (t n -u) α du = ψ(t n , x) -ψ(t n-1 , x) Γ(2 -α)∆t α , (4.30)
while for the history part C hist (t n , x), we integrate by part and get

C hist (t n , x) = 1 Γ(1 -α) ψ(t n-1 , x) ∆t α - ψ(t 0 , x) t α n -α tn-1 0 ψ(u, x) (t n -u) 1+α du . (4.31)
The sum-of-exponentials expansion can be used to approximate the kernel t -1-α . More precisely, for a given absolute error ε and for α, there exist some positive real numbers s and w , with = 1, • • • , N exp (N exp is the number of exponentials), such that

1 t 1+α - Nexp =1 ω e -s t ≤ ε, for all t ∈ [∆t, T ]. (4.32)
By replacing in (4.31) the kernel t -1-α by its sum-of-exponentials approximation in (4.32), we deduce

C hist (t n , x) ≈ 1 Γ(1 -α)   ψ(t n-1 , x) ∆t α - ψ(t 0 , x) t α n -α Nexp =1 ω U hist, (t n , x)   .
The term U hist, (t n , x) is defined by

U hist, (t n , x) = tn-1 0 e -(tn-u)s ψ(u, x)du,
has a simple recurrence relation

U hist, (t n , x) = e -s ∆t U hist, (t n-1 , x) + tn-1 tn-2 e -s (tn-u) ψ(u, x)du, (4.33)
with U hist, (t 0 , x) = 0. The integral can be calculated by

tn-1 tn-2
e -s (tn-u) ψ(u, x)du ≈ e -s ∆t s 2 ∆t (e -s ∆t -1+s ∆t)ψ n-1 + (1-e -s ∆t -e -s ∆t s ∆t)ψ n-2 .

Finally, the fast approximate evaluation of the Caputo fractional derivative is given by the formula

F C 0 D α tn ψ n = ψ(t n , x)-ψ(t n-1 , x) Γ(2 -α)∆t α + 1 Γ(1-α)   ψ(t n-1 , x) ∆t α - ψ(t 0 , x) t α n -α Nexp =1 ω U hist, (t n , x)   , (4.34)
for n > 0, and where U hist, (t n , x) can be obtained by the recurrence relation (4.33). As shown in [START_REF] Jiang | Fast evaluation of the Caputo fractional derivative and its applications to fractional diffusion equations[END_REF], if ψ(t) ∈ C 2 ([0, t n ]), then the truncated error is

| C 0 D α tn ψ -F C 0 D α tn ψ| ≤ C(∆t 2-α + ε). (4.35)
For the fast evaluation of the L2-1 σ scheme, we approximate the kernel (t n+σ -s) -α , for s ∈ (0, t n ), in (4.29) and use a linear interpolation Π 1,n ψ (s), for s ∈ (t n , t n+σ ), yielding

C 0 D α t ψ(t n+σ ) ≈ 1 Γ(1 -α) tn 0 ψ (s) Nexp =1 w e -s (tn+σ-s) ds + tn+σ tn (Π 1,n ψ(s)) (t n+σ -s) α ds = Nexp =1 w Γ(1 -α) tn 0 ψ (s)e -s (tn+σ-s) ds + 1 Γ(1 -α) (ψ n+1 -ψ n ) ∆t tn+σ tn 1 (t n+σ -s) α ds = Nexp =1 1 Γ(1 -α) w V n + λa 0 (ψ n+1 -ψ n ),
where

λ = ∆t -α Γ(2 -α) , a 0 = σ 1-α , V n = tn 0 ψ (s)e -s (tn+σ-s) ds.
The history part V n can be calculated by using a recursive relation and a quadratic interpolation, namely,

V n = tn 0 ψ (s)e -s (tn+σ-s) ds ≈ tn-1 t0 ψ (s)e -s (tn+σ-s) ds + tn tn-1 (Π 2,n ψ(t)) e -s (tn+σ-s) ds = e -s ∆t V n-1 + A (ψ n -ψ n-1 ) + B (ψ n+1 -ψ n ), (4.36)
where

A = 1 0 (3/2 -s)e -s ∆t(σ+1-s) ds, B = 1 0 (s -1/
2)e -s ∆t(σ+1-s) ds.

Overall, the FL2-1 σ formula for C 0 D α t proposed in [START_REF] Yan | Fast evaluation of the Caputo fractional derivative and its applications to fractional diffusion equations: a second-order scheme[END_REF] is given by

FH D α t ψ n+σ = Nexp =1 1 Γ(1 -α) w V n + λa 0 (ψ n+1 -ψ n ), (4.37)
where, if ψ(t) ∈ C 3 ([0, T ]), an error estimate is

| C 0 D α tn ψ -F H 0 D α tn ψ| ≤ C(∆t 3-α + ε). (4.38)
Meanwhile, the number of exponentials N exp is estimated [START_REF] Jiang | Fast evaluation of the Caputo fractional derivative and its applications to fractional diffusion equations[END_REF] by

N exp = O log 1 ε log log 1 ε + log T ∆t + log 1 ∆t log log 1 ε + log 1 ∆t .
For a fixed accuracy ε, we have N exp = O(log N ) for T 1 or N exp = O(log 2 N ) for T ≈ 1, where N = T /∆t. The resulting algorithm has a nearly optimal complexity, requiring O(N JN exp ) operations and a O(JN exp ) storage for solving the TFNLSE. We remark that a linear system must be solved at each time step, therefore at the same cost as for the standard case α = 1. Other efforts to speed-up the evaluation of weakly singular kernels can be found in [START_REF] Alpert | Rapid evaluation of nonreflecting boundary kernels for time-domain wave propagation[END_REF][START_REF] Ke | A fast direct method for block triangular Toeplitz-like with tridiagonal block systems from time-fractional partial differential equations[END_REF][START_REF] Lu | Fast approximate inversion of a block triangular Toeplitz matrix with applications to fractional subdiffusion equations[END_REF][START_REF] Lubich | Fast convolution for nonreflecting boundary conditions[END_REF][START_REF] Zheng | Approximation, stability and fast evaluation of exact artificial boundary condition for one-dimensional heat equation[END_REF].

Thus, we apply the L1 scheme to approximate the time-fractional Caputo derivative, the second-order central finite difference method to approximate the second-order spatial derivative ∂ xx by ∆ h , and use a linearized scheme to approximate the nonlinear term. We respectively have the direct scheme

(γ R + iγ I )(i C 0 D tn ) α ψ n j = - 1 2 ∆ h ψ n j + V j + β|ψ n-1 j | 2 ψ n j , (4.39)
and the corresponding fast scheme

(γ R + iγ I )(i F C 0 D tn ) α ψ n j = - 1 2 ∆ h ψ n j + V j + β|ψ n-1 j | 2 ψ n j , (4.40)
for 1 ≤ j ≤ J -1. As discussed in [START_REF] Zhang | Efficient numerical computation of time-fractional nonlinear Schrödinger equations in unbounded domain[END_REF], the convergence orders of the schemes (4.39) and (4.40) are respectively O(h 2 + ∆t) and O(h 2 + ∆t + ε), if the solution ψ is smooth enough. We point out that the choice of ε is far smaller than (h 2 +∆t) to not affect the convergence order. In addition, FFT-based pseudo-spectral schemes (SP) can also be easily adapted from Sections 2 and 3. Another possible scheme is to use the relaxation technique introduced by Li et al. in [START_REF] Li | Unconditionally convergent L1-Galerkin FEMs for nonlinear time-fractional Schrödinger equations[END_REF] for the direct method

(γ R + iγ I )(i C 0 D tn ) α ψ n j = - 1 2 ∆ h ψ n j + V j + β| ψn j | 2 ψ n j , (4.41)
and the fast scheme

(γ R + iγ I )(i F C 0 D tn ) α ψ n j = - 1 2 ∆ h ψ n j + V j + β| ψn j | 2 ψ n j , (4.42)
with ψn j = 2ψ n-1 j -ψ n-2 j . For a smooth enough solution, the truncated errors for the schemes (4.41) and (4.42) are respectively O(h 2 + ∆t 2-α ) and O(h 2 +∆t 2-α +ε). The proof of the convergence above requires to state a nontrivial discretized fractional Gronwall inequality [START_REF] Li | Analysis of L1-Galerkin FEMs for time-fractional nonlinear parabolic problems[END_REF].

We now consider the high-order scheme for the TFNLSE, written at time t n+σ , n = 0, 1, 2..., and for any ψ(t) ∈ C 2 ([0, t N ]), based on the approximation

ψ(t n+σ ) = σψ n+1 + (1 -σ)ψ n + O(∆t 2 ), n = 0, 1, • • • , N -1.
The time-fractional Caputo derivative is approximated by the FL2-1 σ formula and we adapt the ReFD method for the nonlinear interaction term. The TFNLSE equation is then approximated at (t n+σ , x j ) by For a smooth enough solution, the truncated error for (4.43) is O(h 2 + ∆t 2 + ε), fixing σ = 1 -α/2. Let us remark that time-splitting schemes are currently being developed for fractional PDEs [START_REF] Cao | Time-splitting schemes for fractional differential equations I: smooth solutions[END_REF]. This is probably an interesting future direction to investigate for designing fast TSSP for the TFNLSE, similarly to the schemes presented in Sections 2 and 3. We finally point out that for fractional subdiffusion problems, an essential feature is that the solution always lacks some smoothness near the initial time although it would be smooth away from t = 0 (see e.g. the discussions in [START_REF] Brunner | Numerical simulations of 2D fractional subdiffusion problems[END_REF][START_REF] Jin | An analysis of the L1 scheme for the subdiffusion equation with nonsmooth data[END_REF][START_REF] Mustapha | An implicit finite difference time-stepping method for a subdiffusion equation with spatial discretization by finite elements[END_REF][START_REF] Sakamoto | Initial value/boundary value problems for fractional diffusion-wave equations and applications to some inverse problems[END_REF][START_REF] Stynes | Error analysis of a finite difference method on graded meshes for a time-fractional diffusion equation[END_REF]). The expected order of convergence in time [START_REF] Liao | Sharp error estimate of nonuniform L1 formula for time-fractional reaction-subdiffusion equations[END_REF] needs the use of nonuniform mesh in the temporal direction, such as a graded mesh t n = T (n∆t) γ , for a well-chosen value of the parameter γ > 0.

     (1 -σ)u n+σ j + σu n-1+σ j = β|ψ n j | 2 , (γ R + iγ I )(i F C 0 D tn ) α ψ n+σ j = - 1 2 (σ∆ h ψ n+1 j + (1 -σ)∆ h ψ n j ) + (V j + u n+σ j )(σψ n+1 j + (1 -σ)ψ n j ).

Numerical examples

Example 4.1. For our first 1d example, we consider the quadratic potential V (x) := x 2 /2 and a purely linear TFNLSE, i.e. β := 0. The initial data is taken as the ground state of the standard linear Schrödinger equation with harmonic trap: ψ 0 (x) := e -x 2 /2 /π 1/4 . The corresponding mass is then equal to 1 and the energy level is 1/2. The final time of computation is T = 100 (except for α = 0.25 where T = 800) for the bounded spatial domain D :=] -10, 10[. The discretization parameters are ∆t = 10 -3 and h = 10 -2 . We use the fast second-order scheme (4.43) (for ε = 10 -11 ). Figs. 5(a)-(c) show the evolution of the density function |ψ| 2 , its mass N and energy E 2 , for (γ R + iγ I ) = 1 (no dissipation case) and three values of the fractional power, i.e. α = 0.25, 0.5, 0.75. We observe that the Gaussian solution is modified according to α, the limit of the mass N being indeed α -2 as predicted by formula (4.24), and lim t→+∞ E 2 (t) = α -2 /2 from the numerical simulations. An approximation of the limit value of the mass can be obtained through an average over [0; T ] of the numerical masses, yielding the estimates 16.93, 4.04, 1.78 approximating the theoretical limits 16.38, 4.04, 1.78 for α = 0.25, 0.5, 0.75, respectively. We observe that the solution fluctuates according to t, gaining/loosing some mass alternatively. The solution also disperses slightly more as α goes to zero (this is more clearly seen in Example 4.2). A formal way to understand this property is by considering the dispersive relation (4.26) which can also be rewritten and expanded for large spatial frequencies |k| as

ω = |k| 2/α (2(γ R + iγ I )) 1/α (1 + 2β |k| 2 |A| 2 ) 1/α ∼ |k| 2/α (2(γ R + iγ I )) 1/α (1 + 2β α|k| 2 |A| 2 ). (4.44)
Therefore, we can roughly interpret α as s -1 in Eq. (3.11) up to a multiplicative constant (that may also be complex-valued). Then, superdispersion formally corresponds to α < 1 and subdispersion to α > 1. Nevertheless, the situation here shows that the modeling is more complex since a fractional laplacian also affects the nonlinearity, the corresponding operator acting as a fractional derivative operator for α < 1 and a regularizing integral operator for α > 1. For longer times, the solution tends to a stationary state for the TFNLSE at fixed α which is expected to be the ground state. To end this first example, we compute on Fig. 5(d) the same problem for α = 0.5 but with: (γ R + iγ I ) = i 1-α , which means that we add some dissipation to the model. Therefore, after a sufficiently long time the solution converges to zero as it can be directly observed. Example 4.2. For this second example, we consider the TFNLSE with a quadratic potential V (x) = x 2 /2, and first with a focusing cubic nonlinearity β < 0. The initial data is given by

ψ 0 (x) := A √ -β sech (A(x -x 0 ))e i(νx+θ0) , x ∈ D,
with A = 2, x 0 = ν = θ 0 = 0. In the standard situation α = 1, and without potential term i.e. V = 0, then the mass and energy can be computed analytically [START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF] as , respectively, for β = -1. As already noticed above, the superdispersion effect is visible for smaller values of α and some decoherence effect in the wave field seems to appear. In the case of superdispersion, then there is no longer a compensation between the laplacian and the nonlinear term as for the standard case, the fractional time derivative introducing some dominant dispersion thanks to (4.44). Furthermore, the energy globally grows with the time to pass from negative to positive values, E 2 being larger for a given fractional order than for the integer order case. The energy increases as α gets smaller. The behavior of the mass is difficult to predict, but the value is always smaller than for the standard case. Comparing Figs. 6(a) and 6(c), we see that the curves for the solutions are exactly the same, up to the scaling factor between the two values of β, which is also consistent with formulae (4.45) for the standard case (but with a time-dependency for the fractional situation). Fig. 6(d) finally reports the case where (γ R + iγ I ) = i 1-α for α = 0.5 and β = -1. Compared with Fig. 6(a), we clearly observe the effect of the dissipation term even if the solution is actually not zero for longer times, probably because of the nonlinear interactions.

N (t) := A √ -β , E 2 (t) := Aν 2 β + A 3 3β , ( 4
To complete the simulations, we consider now a defocusing nonlinearity β = 1 in Figs. 7(a)-(c) for α = 0.5, 0.75 and (γ R +iγ I ) = 1 and i 1-α . The simulation parameters are the same as previously, but the computational domain is ] -15; 15[. The behavior of the solution shows again some superdispersion effects, the evolution of the mass and energy being difficult to describe at first sight. When a dissipation term is active (see Figs. 7(c)), then we see that the solution tends to zero for large enough times t (the computational domain is ] -10; 10[ here).

Conclusion and perspectives

We proposed a few numerical methods for time or space fractional nonlinear Schrödinger equations with some applications in Bose-Einstein condensation. In particular, we focused on the physical properties that a scheme should fulfill at the discrete level. Numerical simulations illustrate the purpose of the paper, trying to highlight some issues for the time fractional case. Building accurate, efficient and stable schemes for the FNLSE remains new and open a lot of future important directions to investigate. Among them, let us mention the extension of some schemes to the FNLSE with time and space varying fractional powers [START_REF] Bhrawy | Highly accurate numerical schemes for multi-dimensional space variable-order fractional Schrödinger equations[END_REF], the simulation of coupled systems of FNLSE which can also integrate some nonlocal integral interactions like dipolar or Coulomb potentials, the development of 3d parallel solvers. Finally, one important point which still needs to be deeply understood concerns the link between the physics and the fractional quantum models, most particularly for modeling BECs.

Funding

Xavier Antoine thanks the support from the LIASFMA (University of Lorraine), ANR-

Figure 1 .

 1 Figure 1. Density |ψ(x, t)| 2 at different times t and evolution of the mass N (t) and energy Es(t) for different values of the fractional powers s for Case 1 in Example 3.1.

Figure 2 .

 2 Figure 2. Evolution of the density |ψ(x, t)| 2 , mass N (t) and energy Es(t) for different fractional powers s for Case 2 in Example 3.1.

Figs. 3 iCases 1 - 3 .

 313 Figs. 3 i) and ii) show the contour plot of the density |φ s0gs | 2 of these two types of ground states. Then, we consider the following six cases:Cases 1-3. Choose φ s0 gs as Type 1, let x 0 = (0, 0) T , and only change the fractional order in the SFNLSE (3.8) for the dynamics as s = 0.8, s = 1.2 and s = 1.4, respectively. The other parameters are kept unchanged. Case 4. Choose φ s0 gs as Type 2, and only shift the initial center of mass, i.e., we set x 0 = (5, 5) T . The other parameters are kept unchanged. Case 5. Choose φ s0 gs as Type 2 and set x 0 = (0, 0) T . We perturb the trapping frequency in the y-direction, i.e., we set γ y = 1.3. The other parameters are kept unchanged. Case 6. Same as Case 5, but now perturb the trapping frequency in both x and y-directions, i.e., we set γ x = γ y = 1.3.

Figure 3 . 3 Figure 4 .

 334 Figure 3. Initial data of Type 1 (i)) and Type 2 (ii)) and the evolution of the energy and mass for Cases 1 to 6 in Example 3.2.

Figure 5 .

 5 Figure 5. Evolution of the density |ψ(x, t)| 2 , mass N (t) and energy E 2 (t) for the fractional powers α = 0.25, 0.5 and 0.75, and for the linear case β = 0 for Example 4.1.

1 Figure 6 .

 16 Figure 6. Evolution of the density |ψ(x, t)| 2 , mass N (t) and energy E 2 (t) for the fractional powers α = 0.5 and 0.75 for Example 4.2 (focusing nonlinearity).

1 Figure 7 .

 17 Figure 7. Evolution of the density |ψ(x, t)| 2 , mass N (t) and energy E 2 (t) for the fractional powers α = 0.5 and 0.75 for Example 4.2 (defocusing nonlinearity).

resumes the main properties.

	Method	TSSP	CNFD	ReFD	TSFD
	Time Reversible	Yes	Yes	Yes	Yes
	Time Transverse Invariant	Yes	No	No	Yes
	Mass Conservation	Yes	Yes	Yes	Yes
	Energy Conservation	No	Yes	Yes 1	No
	Dispersion Relation	Yes	No	No	Yes
	Unconditional Stability	Yes	Yes	Yes	Yes
	Time Accuracy	2 nd	2 nd	2 nd	2 nd
	Spatial Accuracy	spectral	2 nd	2 nd	2 nd
	Explicit Scheme	Yes	No	No	No
	Memory Storage at tn	O

  for s = 1.4).

		s=0.2					s=0.8					s=1.3
				t=0					t=0					t=0
	0.09			t=4	0.09				t=3.6	0.09				t=4
				t=7					t=7					t=7
		×10 -5		t=9					t=10					t=10
	0.06	9			0.06	0.1044				0.06	0.1044
		6				0.1038					0.1038
	0.03	3			0.03	0.1032				0.03	0.1032
		0				0.1026					0.1026
	0	7.8 8.1 8.4 8.7		0	-0.8 -0.4 0 0.4 0.8		0	-0.8 -0.4 0 0.4 0.8
	-8	0	x	8	-8		0	x	8	-8		0	x	8
		s=1.8					Es(t)			1.001		N (t)
	0.09			t=0 t=0.88	15.8									s=0.2 s=0.8
				t=1.23										s=1.3
				t=2.17			s=0.2							s=1.8
	0.06	0.104 0.106			15.7		s=1.3 s=0.8			1
							s=1.8
	0.03	0.102
		0.1
	0	-1.2 -0.6 0 0.6 1.2		15.6					0.999
	-8	0	x	8	0	3	t	6	9	0	3	t	6	9

  12-MONU-0007-02 BECASIM and the hosting of the Beijing CSRC. Qinglin Tang thanks the support from the Zheng Ge Ru Foundation and the Singapore Ministry of Education Academic Research Fund Tier 2 R-146-000-223-112. Jiwei Zhang gratefully acknowledges the support of NSFCs (91430216 and U1530401).