On a Fractional Schr\"odinger equation in the presence of Harmonic potential

In this paper, we establish the existence of ground state solutions for a fractional Schr\"odinger equation in the presence of a harmonic trapping potential. We also address the orbital stability of standing waves. Additionally, we provide interesting numerical results about the dynamics and compare them with other types of Schr\"odinger equations. Our results explain the effect of each term of the Schr\"odinger equation : The fractional power, the power of the nonlinearity and the harmonic potential.

In recent years, a great attention has been focused on the study of problems involving the fractional Laplacian, which naturally appears in obstacle problems, phase transition, conservation laws, financial market. Nonlinear fractional Schrödinger equations have been proposed by Laskin [16,17] in order to expand the Feynman path integral, from the Brownian like to the Lévy like quantum mechanical paths. The stationary solutions of fractional nonlinear Schrödinger equations have also been intensively studied due to their huge importance in nonlinear optics and quantum mechanics [16,17,12,10]. The most interesting solutions have the special form: ψ(x, t) = e −iλt u(x), λ ∈ R, u(x) ∈ C.
They are called the standing waves. These solutions reduce (1) to a semilinear elliptic equation.
In fact, after plugging (5) into (1), we need to solve the following equation The case s = 1 has been intensively studied by many authors (See [20]). There also exist a considerable amount of results concerning the standing waves of fractional Nonlinear Schrödinger equations without the harmonic potential, we refer the readers to [3,6,5,8,11,13,18,19] and the references therein.
In this paper, we mainly focus on the solutions to (6). To the best of our knowledge, our results are new and will open the way to solve other class of fractional Schrödinger equations. This paper has two main parts: In the first part, we address the existence of standing waves through a particular variational form, whose solutions are called ground state solutions. We prove the existence of ground state solutions (Theorem 2.1), and show some qualitative properties like monotonicity and radiality (Lemma 3.2). We also proved that the ground state solutions are orbitally stable (Def 4.1, Theorem 2.2) if we have the uniqueness of the solutions for the Cauchy problem (1) (Theorem 4.1). We have also addressed the critical case σ = 2s N , which is consistent with the case s = 1 in [9]. The second part of this article deals with the numerical method to solve (1) and to establish the existence of ground state solutions as well as to establish the optimality of our conditions. In this part, we were not only able to show the existence of ground state solutions for 0 < σ < 2s N but we also gave a constrained variational problem ((62)-(63)), which was crucial to find the standing waves for the subcritical 2s N ≤ σ < 2s N −2s . The numerical results provided a good explanation of the effect of s on the ground state solution. To reach this goal, we showed the ground state solution is continuous and decreasing with respect to s in L 2 and L ∞ norm (Figure 2), which is a similar phenomenon to [15]. Besides, like GrossPitaevskii Equation [9], we examined the convergence property of λ c . It turned out this convergence property also holds true in our case. Second, we checked the stability of ground state solutions for different s. If we add a small perturbation to the initial condition, for different s, the absolute value of the solution will always have periodic behavior, which shows the orbital stability ( Figure 6). Furthermore, surprisingly, when s becomes smaller, the stability is worse, which means the oscillation amplitude in the periodic phenomenon becomes larger (Figure 8,9(a)). We then address the case where the harmonic potential is not radial, and we obtained non radial symmetrical ground state solution ( Figure  5(a),5(b)). Finally, we provided interesting numerical results for the time dynamics of FNLS.
The main difficulty of constructing ground state solutions comes from the lack of compactness of the Sobolev embeddings for the unbounded domain R N . However, by defining an appropriate function space, in which the norm of the potential is involved, we "recuperate" the compactness (see Lemma 3.1). This fact, combined with rearrangement inequalities are the key points to prove the existence of ground state solutions. In the numerical part, the presence of the harmonic potential term is challenging. In fact, one can't take Fourier transfrom directly on both sides of the equation like [15] because we have nonlinear term. Different from [9], we also can't use finite difference directly since fractional Laplacian is not a local term. Consequtently, we opted idea from [7] and use time splitting method. By our splitting, we can obtain specific solutions in each small step and also preserve the mass (3). For the ground state solutions, the classical Newton's method [9] is too slow because we have to deal with fractional Laplacian. To overcome this, we borrow idea from [2] and use normalized gradient flow (NGF) to find the ground state solutions. Moreover, for the case 2s N ≤ σ ≤ 2s N −2s , we have noticed that the energy in the original variational problem can not be bounded from below, therefore, we present a new constrained variational problem ((62)-(63)) to establish the existence of ground state solutions.
The paper is organized as follows. In section 2, we give our main results about the existence of ground state solutions and orbital stability of standing waves. In section 3, we provide the proof of the existence. Then, in section 4, we discuss the orbital stability. In section 5, we use Split-Step Fourier Spectral method to solve (1) numerically. In section 6, instead of using common iterative Newton's method, we use the NGF method to find ground states when 0 < σ ≤ 2s N −2s . Finally, in section 7, we present our numerical results for the dynamics (1) and compare them with other kinds of Schrödinger equations ( [9,15]).

Main results
We use a variational formulation to examine the solution to (6). First, note that if λ = 0, we can find solutions u(x) to (6) from the critical points of the functional J : H s (R N ) −→ R defined as: where . L 2 is the L 2 -norm and ∇ s u 2 is defined by with some normalization constant C N,s . We can derive (7) by multiplying smooth enough test function v(x) on both sides of (6) and taking the integral over x. However, instead of directly finding the critical points of (7), we consider a reconstructed variational problem, which can help us to find solutions with different λ and any energy. Specifically, for a fixed number c > 0, we need to solve the following constrained minimization problem. with where is a Hilbert space, with corresponding natural inner product.
We claim that for each minimizer u(x) of the constrained minimization problem (8), there exists some λ such that (u(x), λ) is a solution to (6). To prove the claim, we first consider λ as a Lagrange multiplier, then we define The minimizer to problem (8) must be the critical point of (11), satisfying: and where (12) implies (6) and (13) implies (9). In this paper, we will mainly focus on the minimizers of problem (8). The following theorem discusses the existence of such minimizers.
Remark 2.1 The condition 0 < σ < 2s N is important in our proof of the existence of minimizers. For the critical case σ = 2s N , we were able to obtain interesting results (section 7). After we construct the ground state solutions, we further investigate their stability. By the definition of (5), the ground state solution moves around a circle when time changes. Therefore we consider and prove the orbital stability of ground state solution (Def 4.1).
Theorem 2.2 Suppose that 0 < σ < 2s N and (1) has a unique solution with conserved mass (3) and energy (4), then the ground state solutions constructed in Theorem 2.1 are orbitally stable.

The minimization problem
In this section, we will establish the existence of ground state solutions of (6), the main difficulty comes from the lack of compactness of the Sobolev embeddings. Usually, at least when potential in (1) is radially symmetric and radially increasing, such a difficulty is overcame by considering the appropriate function space. More precisely, we have , which implies Σ s (R N ) can be embedding into H s . On the other hand, by Sobolov embedding theorem, H s (R N ) can be compactly em- By the classical Sobolev embedding theorem, for any fixed R, H s (|x| < R) is compactly embedded in L 2 (|x| < R). Therefore, for any bounded sequence in Σ s , we choose R n > 0 → ∞ and for each n, pick out the subsequence that converges in L 2 (|x| < R n ) from former convergence sequence in L 2 (|x| < R n−1 ), finally using the diagonal method combined with (14) we find the convergence sequence in L 2 (R N ).
Then we have a lemma showing the existence of I c and boundedness of minimizing sequence.
N , then I c > −∞ and all minimizing sequences of (8) are bounded in Σ s (R N ).
Proof. First, we prove that J (u) is bounded from below. Using the fractional Gagliardo-Nirenberg inequality [12], we certainly have for some positive constant K, where θ = N s 2s(σ+1) .
On the other hand, let > 0, and p, q > 1 such that 1 p + 1 q = 1, then, using Young's inequality, one gets Combining (15) and (16), we obtain for any u ∈ S c , . Hence, from (17) we get: Then we choose small enough in (20) to make 1 2 − > 0, which implies that I c > −∞ and that for all minimizing sequences {u n }, J (u n ) is bounded from above, which implies {u n } is bounded in Σ s (R N ) by (20). Now, we can use compactness(Lemma 3.1) and boundedness(Lemma 3.2) to prove our existence Theorem 2.1.

Proof.
Let {u n } be a minimizing sequence of (8). By Lemma 3.2, {u n } is bounded in Σ s (R N ). Up to a subsequence, there exists u such that u n converges weakly to u in Σ s (R N ).
Since 2σ N −2s , we can further prove that u n will converge strongly to u in L 2 (R N ) and L 2σ+2 (R N ) (Lemma 3.1). In particular, On the other hand, thanks to the lower semi-continuity, we have xu which yields u is a minimizer. The second step consists in constructing a nonnegative, radial and radially decreasing minimizer. First, note that: which implies J (|u|) ≤ J (u). Then we use the Schwarz symmetrization [14]. We construct a symmetrization function u * , which is a radially-decreasing function from R N into R with the property that It's well-known [14] that Besides, from [9], [1], we also have Combining (23) and (24), we obtain Remark 3.1 By (21), and weakly convergence, we can also see xu (25) By [9], (26) implies u = |u| a.e. and (27) implies |u| = |u| * a.e..

Orbital stability
In this section, we will deal with the orbital stability of the ground state solutions. Let us introduce the appropriate Hilbert space: , which is a Hilbert space. In term of the new coordinates, the energy functional reads , we can also get J (ω) remains as a constant with time t if ω(t, x, v) is a solution to (1).
Then, for all c > 0, we set a similar constrained minimization problem where S c is defined by: We also introduce the following sets Proceeding as in [3,13], we have the following lemma: N , then the following properties hold true: (i) The energy functional J and J are of class C 1 on Σ s (R N ) and Σ s (R N ) respectively. (ii) There exists a constant C > 0 such that . (iii) All minimizing sequences for I c are bounded in Σ s (R N ) and all minimizing sequences Proof. (i) We follow the steps of Proposition 2.3 [13] by choosing g(x, t) = −t σ . For any u, v ∈ Σ s (R N ), we can see the last term of functional is of class C 1 on Σ s (R N ). Then by the definition of Σ s (R N ) (see (10)), the first two terms of the functional are of class C 1 on Σ s (R N ).
For the last term, by Hölder's inequality Therefore, there exists C > 0 such that (iii) This is a direct result of Lemma (3.2).
(iv) Let c > 0 and let {c n } ⊂ (0, ∞) such that c n → c. It suffices to prove that I cn → I c . By the definition of I cn , for any n there exists u n ∈ S cn such that From (iii), there exists a constant C 1 > 0 such that for all n, we have Set v n = c cn u n , then, for all n ∈ N, we have v n ∈ S c and u n − v n Σs( We deduce by part (ii) that there exists a positive constant K := K(C 1 ) such that From (29) and (30) we obtain Then, from (28) and (31), we obtain Combining this with the fact that lim Set v n = cn c u n , then v n ∈ S cn , there exists a constant L = L(C 2 ) such that Combining this with (28), we obtain Since lim n→∞ c n = c, we have lim sup n→∞ I cn ≤ I c .
It follows from (32) and (33) that lim n→∞ I cn = I c .
(v) This is a direct result of Remark 3.1.
(vi) First, we can see Σ s (R N ) ⊂ Σ s (R N ), and any ω ∈ Σ s (R N ), we have Second, for any ω ∈ Σ s , we have from which we can easily obtain Combine (34) and (35), we finally have I c = I c .
Now, for a fixed c > 0, we use the following definition of stability (see [4]) • For all ω 0 ∈ O c and ε > 0, there exists δ > 0 such that for all ψ 0 ∈ Σ s (R N ), we have where ψ denotes the solution of (1) corresponding to the initial data ψ 0 .

Proof.
The proof is by contradiction: Suppose that O c is not stable, then there exists 0 > 0, ω 0 ∈ O c and a sequence Φ n 0 ∈ Σ s (R N ) such that ω 0 − Φ n 0 Σs(R N ) → 0 as n → ∞, but for some sequence {t n } ⊂ R, where Φ n (t, .) is the unique solution of problem (1) corresponding to the initial condition Φ n 0 . Let ω n = Φ n (t n , .) = (u n , v n ) ∈ Σ s (R N ). Then, since ω ∈ S c and J (ω) = I c , it follows from the continuity of . 2 and J in Σ s (R N ) that Φ n 0 2 → c and J (Φ n 0 ) → I c , n → ∞. Thus, we deduce from (36) that Since {ω n } ⊂ Σ s (R N ), it is easy to see that {|ω n |} ⊂ Σ s (R N ). On the other hand, Lemma 4.1 (iii) and proof of Lemma 3.2 imply that {ω n } is bounded in Σ s (R N ) and hence, by passing to a subsequence there exists ω = (u, v) ∈ Σ s (R N ) such that Now, by a straightforward computation we obtain Thus, we obtain I c = lim Hence It follows from (39), (40) and (41) that which is equivalent to say that The boundedness of ω n in Σ s (R N ) and (42) imply that |w n | is bounded in Σ s (R N ). By using a similar argument to Lemma 3.1, there exists ϕ ∈ Σ s (R N ) such that |ω n | → ϕ in Σ s (R N ) and ϕ L 2 (R N ) = c with J (ϕ) = I c .

Numerical method for Fractional NLS with harmonic potential
In this section, we consider numerical methods to solve (1) and introduce the Split-Step Fourier Spectral method. First, we truncate (1) into a finite computational domain [−L, L] N with periodic boundary conditions: Let τ > 0 be the time step, and define the time sequence t n = nτ for n ≥ 0 and the mesh size h = 2L/J, where J is a positive even integer. The spatial grid points are where j is a N -dimension integer vector with each component between 0 and J. Denote ψ n j as the numerical approximation of the solution ψ(x j , t n ). By the definition of fractional Laplacian in (2), we use the Fourier spectral method for spatial discretization. Hence, we assume the ansatz: where we introduce the Split-step Fourier Spectral method. The main idea of this method is to solve (48) in two splitting steps from t = t n to t = t n+1 : First, by multiplying ψ * on both sides of (51) and subtracting it from its conjugate, we obtain |ψ(x, t)| = |ψ(x, t n )| for any t ∈ [t n , t n+1 ), therefore, (51) can be simplified to Second, taking Fourier transform on both sides of (52), we get We use the second order Strang splitting method with (53) and (54) as follows: where j comes from (49) and n ≥ 0. For n = 0, initial condition (48) is discretized as: This method has spectral-order accuracy in space and second order in time. Similar to [7], this method preserves discrete mass corresponding to (3) defined as 6 Numerical method to solve ground state solutions To find ground state solutions, we have to solve the following equation corresponding to u(x): As discussed previously, for σ < 2s N , we can solve (8)- (10) to find a solution to (60). In order to calculate the minimizer of J (u) in S c , we use normalized gradient flow method (NGF) [2]. We first apply the steepest gradient decent method to the energy functional J (u) without constraint. Then we project the solution back onto the sphere S c to make sure that the constraint ||u|| L 2 (R N ) = c is satisfied.
Thus, for a given sequence of time 0 = t 0 < t 1 < ... < t n with fixed time step τ , we compute the approximated solution u (n) of the partial differential equation combined with the projection onto S c at each step. Specifically, .
Here, we use semi-implicity time discretization scheme: to discretize fractional laplacian, where K, µ k are defined in (50) . Therefore in each step, we solve : where M n is defined in (59), j comes from (49) and n ≥ 0. For n = 0, we guess a starting function and discretize it as (58). We need to notice that we can only solve (8) for σ < 2s N . If σ ≥ 2s N , u 2σ+2 can not be bounded by · Σs , which will cause I c = −∞ in S c . However, we can use another constrained variational form to find standing waves to (1) for 2s N ≤ σ < 2s N −2s . For σ < 2s N −2s , we define the following constrained minimization problem: with T c = u ∈ Σ s (R N ) : u 2σ+2 = c , where Σ s (R N ) is defined as (10). For u ∈ T c , we have the estimate where C R,δ depends on R, δ. Hence, for ω < 0, if we choose R large enough to make 1 2 + ω R 2 > 0, we get for u ∈ T c . Besides, if ω > 0, (64) is greater than 0. Similar to Lemma 3.2 and Theorem 2.1, there exists a local minimizer for (62) with any ω and c. Now we see λ as the Lagrange multiplier like (11) but with a different functional then we can have the critical points u * and λ * satisfy by ∂K * (u) ∂u = 0. If ω > 0, by multiplying u * on both sides of (65) and taking the integral, we can see which implies λ > 0. Therefore, when ω > 0, we can define u ω,c = λ 1/2σ u * and obtain which means u ω,c e iωt is one standing wave solution to (1). We need to mention (66),(67) actually showed that we can find a ground state solution by solving (63) if L c = K(u * ) > 0. In fact, we have K(u * ) > 0 with ω < 0 but not very small. This is related with the smallest eigenvalue of (−∆) s + |x| 2 ( [9]). Now, for 2s N ≤ σ < 2s N −2s , we use NGF method and semi-implicity time discretization scheme to solve constrained problem (62). Similar to (61), the scheme is where M n 2σ+2 is discrete L 2σ+2 norm and j comes from (49) and n ≥ 0. For n = 0, we guess a starting function and discretize it as (58).

Numerical results
In this section, we show some numerical results, which can help us understand the ground state solution and also illustrate theoretical results. We have mainly investigated: 1. Ground state solutions with different s. 2. Ground state solutions with non-symmetric potentials. 3. Stability and dynamcis.

Numerical results of ground state solutions
First, we solve (8) numerically by (61) in one dimension N = 1 for the case s = 0.8 and σ = 1 to obtain a ground state solution u 0 (x). From figure 1(a) we see the u 0 (x) is radially decreasing as Remark 3.2. Second, we put u 0 (x) into the (1) as initial condition and investigate time evolution of standing waves (Figure 1(b)-1(d)). As expected, we see |ψ(x, t)| is conserved and the real and imaginary part of solution change periodicly with time t. By Theorem 2.1, we can obtain the existence of ground state solutions with σ < 2s N . We change s but keep σ = 1 to obtain ground state solutions with different s. From figure 2(a), we can see when s approaches to 0.5, the ground state solution becomes peaked with faster spatial decay. This is a similar result to the case without potential [15]. We also check · Σ 1 (R) of ground state solutions when s → 0.5 in (2(b)), whose growth shows regularity of ground state solution becomes worse. Then, we test another two things. The first is the relation between the constrained minimal energy in (8) and s. We calculate the discrete energy by From figure 4(a), we find the energy's dependence (E(s)) on s is monotonic. When s approaches to 0.5, the energy will approach to −∞ because of focusing nonlinear term. There are two reasons. First, we keep the L 2 norm of u (we test with same mass c), but the potential term becomes small since u gathers around 0. Second, in Lemma 3.2, we need σ < 2s N to bound u 2σ+2 , whose boundedness becomes worse when s → 0.5. This is different from [15], where they didn't use the variation form and keep the L 2 norm.
Second, we test the relationship between mass c and λ c , where c, λ c are mass and Lagrange multiplier corresponding to the minimizer (8)- (10) . By [9], in the case s = 1, there exits λ 0 such that for any σ, lim c→0 λ c = λ 0 .
We also test this with s = 0.8. From figure (4(b)), we can see for different σ, when c → 0, λ c will also converge to a same value λ 0 . Up to now, we only consider the case with radial symmetrical potential. However, when potential is not radially symmetric, we can still find standing waves to (1) using (8). We tried the case where potential is |x| 2 + a sin(2πx) with a = 1 and a = 5. From figure 5(a) and 5(b), we see if we add a nonsymmetrical perturbation to potential, we won't get radially symmetrical ground state solutions.

Dynamics and stability of ground state solution
First, we consider the case s > N σ 2 , which is covered by our Theorem 2.2. From theoretic results and figure 1(b), we can see standing waves preserve |ψ(x, t)| with time t. Therefore, we use |ψ(x, t)| to draw graphs and test its stability. We consider the case s = 0.8 and σ = 1. We first test condition (36) in Theorem 2.2 with initial condition ψ 0 (x) = (1 + e) * u 0 (x), because the scheme is mass preserving, it suffices to test energy preservation, which is showed in figure  6(a),6(b). Then, we test the stability of solution, where e is a constant number. From figure 6(c),6(d), we can see when e = 0.05, the solution almost preserves |ψ(x, t)| as we desired. When e = 0.2, the solution shows large perturbation but still has periodic behavior, similar to [9] and [15]. In figure 7, we compare · Σs distance between ground state solution and perturbed solution. We can see from figure 2(b) and theoretic results that when s → σN 2 , the regularity of ground state solutions becomes worse. This inspires us to investigate its stability relationship with s. By Theorem 2.2, Def 4.1, the orbital stability means we can find ω ∈ O c such that ω − ψ Σs(R) is small when we only have small perturbation in initial condition. This definition is hard to measure. Therefore, instead of checking the exact definition of orbital stability, we test classical stability by comparing distance between perturbed solution and ground state solution using normalized · Σs(R) distance: D(s, t) = u s (x, t) − ψ s (x, t) Σs(R) u s (x, t) Σs(R) .
We test initial condition ψ 0 (x) = 0.9 * u 0 (x) with s = 0.8 and s = 0.6. From figure 8,9(a), as expected, when s is small, its stability seems worse. To be complete, we test D(s, 1) with s between 0.51 and 1 in figure 9(b). We can see D(1, s) increases when s approaches to 0.5, which implies worse stability.  Second, we try to obtain some numerical result when we touch the critical point s = σN 2 . In this case, we can't find the ground state solution through (8) because I c = −∞. However, as we discussed before, we can find a ground state solution related to another constrained minimization problem (62). Here, we try to use the NGF method to find the ground state solution with s = 0.5, σ = 1. We first tried positive ω, but the projection step dominated the process (68). Therefore, we tried ω = −0.5 and find the method does converge to a solution. From figure 10(b)-10(d), we can see |ψ(x, t)| almost preserves with time t with periodical real and imaginary part. We use it to test the finite blow up phenomenon (ψ 0 (x) = 2u 0 (x)) appeared in the case without potential [15], and this also happens with potential (Figure 7.2). We note here we still can't find a perfect ground state solution, the reason might come from when s = 0.5, the stability of (1) is very bad.
Finally, we test some simple time dynamics of FNLS, we let ψ 0 (x) = u 0 (x)e ikx , which changes its phase but not absolute value. If s = 0.8 and k = 1, 20 ( Figure 11(a), 11(b)), the maximum point of |ψ(x, t)| will move along x periodically. We also test the L ∞ norm of ψ(x, t) (Figure 11(c)). We find it decreases first and then approaches to some limits, which is similar to the case without potential [15].