Blow-up dynamics and spectral property in the L 2 -critical nonlinear Schrödinger equation in high dimensions
BBLOW-UP DYNAMICS AND SPECTRAL PROPERTYIN THE L -CRITICAL NONLINEAR SCHR ¨ODINGER EQUATIONIN HIGH DIMENSIONS KAI YANG, SVETLANA ROUDENKO AND YANXIANG ZHAO
Abstract.
We study stable blow-up dynamics in the L -critical nonlinear Schr¨odinger equationin high dimensions. First, we show that in dimensions d = 4 to d = 12 generic blow-up behaviorconfirms the log-log regime in our numerical simulations under the radially symmetric assumption,and asymptotic analysis, including the log-log rate and the convergence of the blow-up profilesto the rescaled ground state; this matches the description of the stable blow-up regime in the 2 d cubic NLS equation.Next, we address the question of rigorous justification of the log-log dynamics in higher di-mensions ( d ≥ d = 5 to d = 12, anda modification of it in dimensions 2 ≤ d ≤
12. This, combined with previous results of Merle-Rapha¨el, proves the log-log stable blow-up dynamics in dimensions d ≤
10 and radially stable for d ≤ Contents
1. Introduction 22. Numerical simulations of the solutions 72.1. Dynamic rescaling method 72.2. Numerical Results 103. Spectral Properties 183.1. The radial case 203.2. The non-radial case 254. Conclusions 30Appendix A: Justification of “log-log” corrections in d = 6 , ...,
12 31Appendix B: Artificial boundary conditions 34Appendix C: Computation of the potentials V and V d = 2 , , Mathematics Subject Classification.
Primary: 35Q55, 35Q40; secondary: 35P30.
Key words and phrases. stable blow-up dynamics, NLS equation, log-log blow-up, spectral property, dynamicrescaling method. a r X i v : . [ m a t h . A P ] M a r K.YANG, S. ROUDENKO, Y. ZHAO Introduction
We consider the Cauchy problem of the L -critical nonlinear Schr¨odinger (NLS) equation: (cid:40) iu t + ∆ u + | u | d u = 0 , ( t, x ) ∈ R × R d u ( x,
0) = u ∈ H ( R d ) . (1.1)The local well-posedness of the equation (1.1) in H is given by Ginibre and Velo in [15], seealso [7]. It shows that for u ∈ H ( R d ), there exists 0 < T ≤ ∞ such that there is a uniquesolution u ( t ) ∈ C ([0 , T ) , H ( R d )). We say the solution exists globally in time if T = ∞ , and thesolution blows up in finite time if T < ∞ and lim sup t (cid:37) T (cid:107)∇ u ( t ) (cid:107) L = ∞ .During their lifespan, the solutions u ( x, t ) of the Cauchy problem (1.1) satisfy the followingconservation laws Mass : M ( u ( x, t )) def = (cid:90) | u ( x, t ) | dx = (cid:90) | u ( x ) | dx, (1.2) Energy : E ( u ( x, t )) def = 12 (cid:90) |∇ u ( x, t ) | dx −
12 + d (cid:90) | u ( x, t ) | d dx = E ( u ) , (1.3) Momentum : P ( u ( x, t )) def = Im (cid:18)(cid:90) ∇ u ( x, t )¯ u ( x, t ) dx (cid:19) = Im (cid:18)(cid:90) ∇ u ( x )¯ u ( x ) dx (cid:19) . (1.4)Substituting u ( t, x ) = e it Q ( x ) into (1.1), we obtain a special class of periodic solutions with Q solving − ∆ Q + Q − | Q | d Q = 0 . (1.5)The equation (1.5) exhibits many solutions, even when restricted to those that are smooth andvanishing-at-infinity. In this work we are only interested in real positive solutions. The existenceand uniqueness of a real, positive, vanishing-at-infinity solution of (1.5) are obtained in [4] and[24], and this solution is referred to as the ground state solution, which we also denote by Q .Note that the ground state solution is radially symmetric Q = Q ( r ), moreover, it is exponentiallydecaying at infinity, [4]. For the purpose of this paper, we will simply say Q ∈ H ( R d ). In thedimension d = 1, the ground state solution of (1.5) is explicit Q ( x ) = 3 / sech / (2 x ) . (1.6)While in dimensions d ≥
2, the ground state is not explicit, its properties are known, andvarious numerical methods (e.g., renormalization method or shooting method, see, for example,the monograph [9] and reference therein) produce the ground state Q numerically.The importance of the ground state comes into play when one wants to understand the longtime behavior of solutions. In 1983, Weinstein [53] showed that solutions exist globally in time if (cid:107) u (cid:107) L < (cid:107) Q (cid:107) L . By convexity arguments on a finite variance, it is known that solutions blow upin finite time if E [ u ] <
0, see [51], [58] or [16]. Thus, solutions with (cid:107) u (cid:107) L ≥ (cid:107) Q (cid:107) L may blowup in finite time. Moreover, the minimal mass blow-up solutions happen exactly at the threshold (cid:107) u (cid:107) L = (cid:107) Q (cid:107) L , and the classification of all minimal mass blow-up solutions was done by Merlein [33] (radial) and [34] (general case). We note that all minimal mass blow-up solutions areunstable. In this paper we are interested in stable blow-up dynamics, thus, we consider the initialdata with the mass above the mass of the ground state Q . TABLE BLOWUP IN NLS IN HIGH DIMENSIONS 3
Investigations of the stable blow-up dynamics for the NLS equation (mainly in the two dimen-sional case, i.e., cubic nonlinearity) go back to 1970’s (for example, [5], [48], [17], [47]). In 1986,McLaughlin, Papanicolaou, Sulem and Sulem in [32] (see also [46]) showed that the rate of stableblow-up should be coming from scaling and equivalent to ( T − t ) − by applying the dynamicrescaling method. Previously, Talanov and collaborators (1978), Wood (1984) and Rypdal andRasmussen (1986) suggested a form of ( | ln( T − t ) | / ( T − t )) from a different approach, see [52],[56], [43]. Then, Landman, Papanicolaou, Sulem, and Sulem in [26] (also LeMesurier, Papanico-laou, Sulem and Sulem in [28]) as well as Fraiman [14], and others (see [25], [8], [23], [30]; fora more complete reference list on this subject refer to [46]) concluded from asymptotic analysis,incorporating numerical simulations, that the rate of blow-up is of the form (cid:18) ln | ln( T − t ) | ( T − t ) (cid:19) ,which is now referred to as the “log-log” law. Later in 1990, Landman, Papanicolaou, Sulem,Sulem and Wang in [27] concluded that this log-log rate is a stable blow-up rate in 2d by simulatingthe blow-up dynamics without radial symmetry assumption. Note that numerical tracking of thelog-log correction is extremely difficult as it is reached only at exceedingly huge focusing levels (oforders ∼ iu t + u xx + | u | u = 0). As it was the firstanalytical proof of the “log-log” blow-up, it opened possibilities to study blow-up solutions rigor-ously. In the series of papers [35], [36], [37], [41], [38], Merle and Rapha¨el did a systematic studyof the generic blow-up, considering initial data slightly above the ground state mass and with neg-ative energy, and obtained a detailed description of the dynamics, in particular, the convergenceof the blow-up profiles to the self-similarly rescaled Q and the blow-up rate (cid:18) ln | ln( T − t ) | ( T − t ) (cid:19) .This description was proved for any dimension, provided the so-called “spectral property” holdstrue (see the actual statement below). This spectral property was proved rigorously in the onedimensional case in [35], since the ground state Q is explicit in 1d (1.6) and the spectral propertycould be deduced from some known properties of the second-order differential operators [49]. Forthe dimensions d = 2 , ,
4, Fibich, Merle and Rapha¨el gave a numerically-assisted proof of thespectral property, using the numerical representation of Q , see [11] and [12]. For d ≥
6, theauthors in [11] indicated that the spectral property is indecisive (via their approach), further-more, there was no direct numerical simulation of the blow-up solutions in higher dimensions,even under the radial symmetry assumption (all prior L -critical NLS numerical simulations were K.YANG, S. ROUDENKO, Y. ZHAO in 2d), as it could have been a computationally difficult task at that time. Investigating stableblow-up regime in higher dimensions is the goal of this work.In this paper, we first show that in higher dimensions, 4 ≤ d ≤
12, a generic self-similar blow-up dynamics is also described by the “log-log” law, i.e., our numerical simulations show thatthe blow-up profiles converge to the self-similar ground state, furthermore, we obtain the log-logblow-up rate via derivation as in [46, Section 8] and matching refinement adapted from [1]. As weonly numerically simulate the radial case, we conclude in this first part that this log-log blow-updynamics is radially stable. Secondly, we give a proof of the spectral property for d = 5 , ..., d ≤
12 forthe initial data with negative energy and mass slightly above the mass of the ground state.For consistency, in what follows, we adopt the notation from [11] and [46]: denote by B α aneighborhood of H functions with the L norm slightly above the mass of the ground state Q , i.e., B α = (cid:8) f ∈ H ( R d ) : (cid:82) Q < (cid:82) | f | < (cid:82) Q + α (cid:9) ; the variable x = ( x i ) ≤ i ≤ d is used asthe space variable, r = | x | stands for the radial coordinate and the radial derivative ∂ r f ( r ) wetypically write as f r ; we also denote H r as the space of radial functions in H , and (cid:104) , (cid:105) as thestandard inner product on L . We define the scaling symmetry generator acting on Q by Q = d Q + x · ∇ Q, and also define Q = d Q + x · ∇ Q . (1.7) Spectral Property 1 ([11], [31]) . Let d ≥ . Consider the two real Schr¨odinger operators L = − ∆ + V and L = − ∆ + V , (1.8) where V ( r ) = 2 d (cid:18) d + 1 (cid:19) Q d − rQ r , (1.9) V ( r ) = 2 d Q d − rQ r , (1.10) and the real-valued quadratic form for u = f + ig ∈ H ( R d ) is defined as B ( u, u ) = B ( f, f ) + B ( g, g ) (1.11)= (cid:104) L f, f (cid:105) + (cid:104) L g, g (cid:105) . (1.12) Then there exists a universal constant δ > such that for any u ∈ H ( R d ) , if the followingorthogonal conditions hold (cid:104) f, Q (cid:105) = (cid:104) f, Q (cid:105) = (cid:104) f, x i Q (cid:105) ≤ i ≤ d = 0 , (1.13) (cid:104) Q , g (cid:105) = (cid:104) Q , g (cid:105) = (cid:104) ∂ x i Q, g (cid:105) ≤ i ≤ d = 0 , (1.14) then B ( u, u ) > , or more precisely, B ( u, u ) ≥ δ (cid:18)(cid:90) |∇ u | + | u | e −| x | (cid:19) . (1.15)We note that the commutator-type formulas hold for the operators L and L , which wasobserved in [35], L f = 12 [ L + ( f ) − ( L + f ) ] and L g = 12 [ L − ( g ) − ( L − g ) ] , TABLE BLOWUP IN NLS IN HIGH DIMENSIONS 5 where L + and L − are the standard linearized operators obtained from the linearization of solutionsaround the ground state Q : L + = − ∆ + 1 − (cid:18) d (cid:19) Q d , (1.16) L − = − ∆ + 1 − Q d . (1.17)We mention that the choice of the orthogonal conditions (1.13) and (1.14) in the above SpectralProperty 1 comes from the generalized null space of the matrix operator L = (cid:20) L − − L + (cid:21) (see,e.g., Weinstein [55], Buslaev-Perelman [6]).Our first result is the following: Theorem 1.1.
The Spectral Property 1 holds true in dimensions from d = 5 to d = 10 . TheSpectral Property 1 holds true in the radial case for dimensions d = 11 and d = 12 . We next modify the statement of the Spectral Property 1.
Spectral Property 2.
Let d ≥ and assume the same set up as in the Spectral Property 1,except replace the orthogonal conditions (1.14) by (cid:104) Q, g (cid:105) = (cid:104) ∂ x i Q, g (cid:105) ≤ i ≤ d = 0 . (1.18) Then B ( u, u ) > , or (1.15) is true (i.e., the same conclusion as in the Spectral Property 1 holds). Our second result is about the modified Spectral Property.
Theorem 1.2.
The Spectral Property 2 holds true in dimensions from d = 2 to d = 10 . TheSpectral Property 2 holds true in the radial case for dimensions d = 11 and d = 12 . The Spectral Property 1 added to the work of Merle-Rapha¨el [35], [36], [37], [38], and Fibich-Merle-Rapha¨el [11], now fully completes the proof of the stable blow-up dynamics in dimensions d ≤
10, and radially stable in d ≤
12 with the description of various dynamical features (profile,rate, phase, control of the remainder, etc; we note that the bounded control in the external regionis given in [20]). For completeness of this introduction, we provide a concise statement of thestable blow-up dynamics below.
Theorem 1.3 (Stable blow-up dynamics of (1.1)) . Assume the spectral property holds true, whichis now proved for d ≤ and in the radial case for d = 11 , .There exist universal constants α > and C > such that the following holds true. For u ∈ B α , let u ( t ) be the corresponding H solution to (1.1) on [0 , T ) , the interval of the maximalin forward time existence of u . (i) Description of the singularity: Assume that u ( t ) blows up in finite time, i.e., < T < ∞ .Then there exist parameters ( λ ( t ) , x ( t ) , γ ( t )) ∈ R × R d × R and asymptotic profile u ∗ ∈ L such that u ( t ) − λ ( t ) d Q (cid:18) x − x ( t ) λ ( t ) (cid:19) e iγ ( t ) → u ∗ in L ( R d ) as t → T. Moreover, the blow-up point is finite in the sense that x ( t ) → x ( T ) ∈ R d as t → T. K.YANG, S. ROUDENKO, Y. ZHAO (ii)
Estimates on the blow-up speed: for t close enough to T, either lim t → T |∇ u ( t ) | L |∇ Q | L (cid:18) T − t ln | ln( T − t ) | (cid:19) = 1 √ π , (1.19) or |∇ u ( t ) | L ≥ C ( u ) T − t . The equation (1.19) is referred to as the “log-log” blow-up rate. (iii)
Sufficient condition for “log-log” blow-up : If E ( u ) < and (cid:82) Q < (cid:82) | u | < (cid:82) Q + α ,then u ( t ) blows up in finite time with the “log-log” rate (1.19). More generally, the set ofinitial data u ∈ B α such that the corresponding solution u ( t ) to (1.1) blows up in finitetime < t < ∞ with the “log-log” speed (1.19) is open in H . Remark 1.4.
Our numerical simulations show that for the initial data u of the Gaussian-type, u ∼ Ae − r , or of the ground state-type, u ∼ A e − r , the parameter α in Theorem 1.3(iii) can betaken very large, i.e., α = ∞ . This paper consists of two parts. In the first part, Section 2, we show the results from the directnumerical simulations of solutions in the L -critical NLS equation in dimensions 4 ≤ d ≤
12 bythe dynamic rescaling method for the radially symmetric data. This shows that the blow-up rateis ( T − t ) − , possibly with some correction terms. Applying the arguments from [26] and [28],we show that the correction term indeed exists, which can be obtained via asymptotic analysis q ( t ) ≈ (cid:16) (2 π ) / ln ln | T − t | (cid:17) , and to further confirm it, we apply the approach from [1] to fit thesolution with various functional forms. This leads us to the conclusion that at least for the radialdata in the dimensions d = 4 , ...,
12, the “log-log” blow-up dynamics is generic and stable, seeSection 2.2 for details. In the second part of the paper, Section 3, we revisit the Spectral Property 1and give a numerically-assisted proof of it for dimensions up to d = 10 (general case), and for theradially symmetric case for d = 11 and d = 12; we then also establish the Spectral Property 2.In Appendices, we provide the rest of the numerical simulations in dimensions d = 6 , ...,
12 (thesimulations for d = 4 , V , in high dimensions; and also providecomparison with previous results from [11] of the Spectral Property 1 for d = 2 , , Acknowledgements : The authors would like to thank Gideon Simpson for his visit to GWUand posing some of the questions addressed in this paper as well as for helpful discussion onthis topic. The authors are grateful for the fruitful discussions, the initial guidance for numericsand comments on the initial draft of the paper to Pavel Lushnikov, who also visited GWU inApril 2017, and whose trip was partially supported by GWIMS. S.R. would like to thank PierreRapha¨el for bringing to her attention [12] and for posing the question about the orthogonalitycondition (1.18). K.Y. and Y.Z. would like to acknowledge Yongyong Cai, who hosted their visitto the Computational Science Research Center (CSRC) in Beijing during Summer 2017. S.R.was partially supported by the NSF CAREER grant DMS-1151618 as well as part of the K.Y.’sgraduate research fellowship to work on this project came from the above grant. Y.Z. was partiallysupported by the Simons Foundation through Grant No. 357963. A more general description is E ( u ) <
12 [ P ( u )] M ( u ) , see [35], [36]. TABLE BLOWUP IN NLS IN HIGH DIMENSIONS 7 Numerical simulations of the solutions
In this section, we first review the dynamic rescaling method which is applied in our numericalsimulations. Then, we present our numerical results. These numerical findings, combined withthe analysis in [26] and [28], adapted to our setting, shows the “log-log” blow-up rate.2.1.
Dynamic rescaling method.
We use the dynamic rescaling method (see [32] and [46])for simulating the blow-up solutions of the equation (1.1). The key idea of this method is toappropriately rescale the equation (1.1) in both time and space variables. Then solutions for therescaled equation exist globally in time. We restrict our numerical simulations under the radiallysymmetric assumption. Letting σ = d , where d is the dimension, we set u ( r, t ) = 1 L ( t ) σ v ( ξ, τ ) , where ξ = rL ( t ) , τ = (cid:90) t dsL ( s ) . (2.1)Given a prescribed value L (0) (e.g., L (0) = 1), we obtain the initial condition v ( ξ ) = L (0) σ u ( r ).Then the equation (1.1) becomes iv τ + ia ( τ ) (cid:16) ξv ξ + vσ (cid:17) + ∆ ξ v + | v | σ v = 0 , (2.2)with boundary conditions v ξ (0) = 0 , v ( ∞ ) = 0 , (2.3)where a = − L dLdt = − d ln Ldτ . (2.4)There are various choices for tracking the parameter L ( t ). The first choice is the following: sincewe want to bound (cid:107)∇ u ( t ) (cid:107) L as t → T , we choose the parameter L ( t ) such that the value of (cid:107)∇ ξ v ( τ ) (cid:107) L in the rescaled equation (2.2) remains constant in time, i.e., L ( t ) = (cid:18) (cid:107)∇ ξ v (cid:107) L (cid:107)∇ u ( t ) (cid:107) L (cid:19) /p , p = 2 + 2 σ − d. (2.5)Then, a ( τ ) is, consequently, a ( τ ) = − p (cid:107) v (cid:107) L (cid:90) ∞ | v | σ Im (¯ v ∆ v ) ξ d − dξ. (2.6)An alternative choice for the L ( t ) is to fix the L ∞ norm of v such that (cid:107) v (cid:107) L ∞ remains constant(see [28], [23] for such a choice), thus, set L ( t ) = (cid:18) v (0) (cid:107) u ( t ) (cid:107) L ∞ (cid:19) σ , (2.7)and hence, a ( τ ) rewrites as a ( τ ) = − σ | v (0) | Im (¯ v ∆ ξ v ) | (0 ,τ ) . (2.8)The above two choices are typical rescalings, for example, see [32] and [46].To discuss our numerical scheme, we note that the equation (2.2) is of the form iv τ + ∆ ξ v + N ( v ) = 0 , τ ∈ [0 , ∞ ) , ξ ∈ [0 , ∞ ) , (2.9) K.YANG, S. ROUDENKO, Y. ZHAO where N ( v ) = ia ( τ ) (cid:16) ξv ξ + vσ (cid:17) + | v | σ v. Before we discretize the space variable ξ , we need to map the spatial domain [0 , ∞ ) onto [ − , ξ ∈ [0 , ∞ ). One possible way to do that is to set ξ = κ z − z , where κ is a constant indicatingthe half number of the collocation points assigned on the interval [0 , κ ], and z is the ChebyshevGauss-Lobatto collocation points on the interval [ − , v ( ∞ ) = 0 on the right, we remove the lastChebyshev point, and, consequently, delete the last row and the last column of the matrix M in(2.12). This is similar with the singular behavior at x → ∞ in [3] in which the authors impose Noboundary condition at x → ∞ . The Laplacian operator can be discretized by Chebyshev-Gauss-Lobatto differentiation (for details refer to [45] and [50]). We denote the discretized Laplacianwith N + 1 collocation points by ∆ N . Note that∆ ξ v = v ξξ + d − ξ v ξ . The singularity at ξ = 0 is eliminated by applying the L’Hˆospital’s rule:lim ξ → d − ξ v ξ = ( d − v ξξ . We use the following notation for v : let v ( m ) be the discretized v at the time τ = m · δτ , where δτ is the time step and m is the iteration. The time evolution of (2.9) can be approximated bythe second order Crank-Nicolson-Adam-Bashforth method: i v ( m +1) − v ( m ) δτ + 12 (cid:0) ∆ N v ( m +1) + ∆ N v ( m ) (cid:1) + 12 (cid:0) N ( v ( m ) ) − N ( v ( m − ) (cid:1) = 0 . (2.10)We reorganize (2.10) as (cid:18) iδτ + 12 ∆ N (cid:19) v ( m +1) = (cid:18) iδτ −
12 ∆ N (cid:19) v ( m ) − (cid:0) N ( v ( m ) ) − N ( v ( m − ) (cid:1) , (2.11)which is equivalent to M v ( m +1) = F ( v ( m ) , v ( m − ) . (2.12)Therefore, the time step is updated by v ( m +1) = M − F ( v ( m ) , v ( m − ) . (2.13)The inverse of the matrix M can be calculated only once in the beginning, since M = (cid:0) idτ + ∆ N (cid:1) ,which stays the same.The boundary conditions (2.3) are imposed similar to [32], [50] and [45] as follows: For thehomogeneous Neumann boundary condition ( v ξ (0) = 0) on the left, we substitute the first rowof the matrix M by the first row of the first order Chebyshev differential matrix, and changethe first element of the vector F to 0. Because of the homogeneous Dirichlet boundary condition v ( ∞ ) = 0 on the right, we delete the last row and column of M as well as the last element of thevector F . TABLE BLOWUP IN NLS IN HIGH DIMENSIONS 9
An alternative method is to use predictor-corrector method similar to [9], i v ( m +1)pred − v ( m ) δτ + 12 (cid:16) ∆ N v ( m +1)pred + ∆ N v ( m ) (cid:17) + 12 (cid:0) N ( v ( m ) ) − N ( v m − ) (cid:1) = 0 , (P) (2.14) i v ( m +1) − v ( m ) δτ + 12 (cid:0) ∆ N v ( m +1) + ∆ N v ( m ) (cid:1) + 12 (cid:16) N ( v ( m +1)pred ) + N ( v ( m − ) (cid:17) = 0 . (C) (2.15)The above two schemes lead to the similar results. Numerical test suggests that the scheme(2.14) and (2.15) is slightly more accurate than the scheme (2.10), though it is still second orderin time and doubles the computational time. We mainly use the “predictor-corrector” scheme(2.14) and (2.15) in our simulation.Note that it is sufficient to use N = 256 collocation points. For any sufficiently smooth function f over [ − , f ( x ) can be expanded by the Chebyshev polynomials f ( x ) = ∞ (cid:88) n =0 c n T n ( x ) , c n → ∞ as n → ∞ . (2.16)Therefore, we choose a number N such that the coefficients c n are sufficiently small (generallyreaches the machine’s accuracy 10 − ) for n ≥ N . Figure 1 shows the Chebyshev coefficients | c n | of the numerical solution v in the case d = 4 at different times τ = 0 and τ = 400. One can seethat | c n | reaches the order 10 − around N = 128, and reaches the machine’s accuracy around N = 170. Hence, it is sufficient to choose any N > N . Since the Chebyshev transform can beprocessed by the fast Fourier transform (fft) (e.g., [50], [45] and [3]), we choose N = 256, a powerof 2. Similar study for the decaying behavior of the Chebyshev coefficients c n for multidomainspectral method for Schr¨odinger equation can be found in [3]. Figure 1.
The Chebyshev coefficients for the solution v in the case d = 4 at time τ = 0 (left) and τ = 400 (right). One can see the value | c n | reaches the machine’saccuracy around N = 170.We choose κ = 256 and dτ = 2 × − . We set v (0) = 1 and we choose to fix the value (cid:107) v ( τ ) (cid:107) ∞ ≡ τ . Then the parameter L ( t ) becomes L ( t ) = (1 / (cid:107) u ( t ) (cid:107) L ∞ ) σ , see [28]. Wedecided not to fix the value of (cid:107)∇ v ( τ ) (cid:107) L to be constant in high dimensions since it involvesintegrating (cid:82) ∞ ...ξ d − dξ , which is large when the dimension d increases. These two differentscalings actually lead to the same slope of L ( t ) on the log scale, i.e., (log( L ) / log( T − t )) ≈ , and also to the same decay property of a ( τ ), when we test the cases d = 4 and d = 5. Here, theterm (cid:82) ∞ ...ξ d − dξ can be reduced to the order (cid:82) − δ − ... (1 − z ) . − d dz by calculating the qunatity a ( τ ) from the variable z , where δ = 1 − z N and z N is the second last discretized ChebyshevGauss-Lobatto collocation point, i.e., a ( τ ) = − p (cid:107) v (cid:107) L (cid:90) ∞ | v | σ Im (¯ v ∆ v ) ξ d − dξ = − p (cid:107) v (cid:107) L (cid:90) − δ − | v | σ Im( α d − (1 − z ) . − d (1 + z ) d − . v zz − α d − (1 − z ) . − d (1 + z ) d − . v z +( d − α d − (1 + z ) d − . (1 − z ) . − d v z ) 1 (cid:112) (1 + z )(1 − z ) dz. The second time step v (1) can be obtained by the standard second order explicit Runge-Kuttamethod (RK2). We stopped our simulations at the dimension d = 12, see details and discussionon this in Section 3.2.2. Numerical Results.
Initial Data.
As in [32], we use the Gaussian-type data u ( r ) = Ae − r as the initial data ,which leads to the self-similar blow-up solutions concentrated at the origin (Towerne profiles).We actually work with data u = A e − r d . The reason for not taking u = A e − r as the initialdata is that the amplitude A becomes very large in higher dimensions; while the normalizationterm r d allows us to keep A reasonably small. The typical initial data is given in Table 1. d (cid:107) Q (cid:107) (cid:107) e − r d (cid:107) ˜ A A (cid:107) u (cid:107) . . . . . .
466 765 . . . . . .
81 768 8 . . . . . e
610 7880266 . . . e
711 107056593 . . . . e
812 1586849773 . . . e Table 1.
Details about the initial data used in the dimension d , u = A e − r d , alsothe threshold amplitude for blow-up ˜ A , and the L norms of the ground state Q , e − r d and u .This table lists the mass of the ground state Q , the mass of e − r d , and we also list the thresholdof the amplitude ˜ A for the blow-up v.s. global existence. We list an example of the amplitude We also ran simulations with the ground state-type of initial data u ∼ Ae − r , and the results were similar.For brevity, we only include the discussion on the Gaussian-type data. TABLE BLOWUP IN NLS IN HIGH DIMENSIONS 11 A that we use in our simulations to be specific, though any amplitude A > ˜ A gives the sameresult.Our next Table 2 shows how the quantity (cid:107) v (cid:107) L ∞ ξ is conserved in the rescaled time τ , i.e., wetrack the error E = max τ ( (cid:107) v ( τ ) (cid:107) L ∞ ξ ) − min τ ( (cid:107) v ( τ ) (cid:107) L ∞ ξ ) . (2.17) d E e − e −
10 1 e − e − e − e − e − e −
10 9 e − Table 2.
The error E from (2.17) on the conserved quantity (cid:107) v ( τ ) (cid:107) L ∞ ξ ≡ τ by using the predictor-corrector method (2.14) and (2.15) with δτ = 2 × − .We comment that one should avoid choosing initial data too flat around the origin, for example,something like super-Gaussian data u = A e − r , since this may lead to the collapsing rings insteadof the Towerne profiles which we are trying to track here, for example, see [19], [21], [22], [9].2.2.2. Blow-up rate.
We discuss our numerical results in the case d = 4, the results for thedimensions d = 5 , ...,
12 are similar, we omit them here, and refer the reader to [57]. We firstobserve that L ( t ) ∼ √ T − t, (2.18)see Figures 2. The slope of L ( t ) in the log scale is shown on the left, where we plotted log( T − t )vs. log( L ); it gives a straight line with the slope approximately (the slope is 0 . d = 2 and d = 3 similarcomputations and graphs were obtained in [32]. Figure 2.
The 4d case: the slope of L ( t ) on the log scale (left); the behavior of a ( τ ) (right).Our next task is to obtain the correction term in the rate (2.18). For that, as it is discussed in[46], one has to study the behavior of a ( τ ), defined in (2.2)-(2.4). We show the behavior of a ( τ )in Figures 2 on the right; the graph shows a very slow decay of a ( τ ), which is also similar to the2d case in [32]. The formal asymptotic analysis in [26] (see also [46, Section 8]) investigates the decay rate of a ( τ ). For that, first, substituting v = e iτ − iaξ / w into (2.2), the equation for w is obtained: iw τ + ∆ w − w + 14 b ( τ ) ξ w + | w | σ w = 0 , (2.19)where b ( τ ) = a + a τ = − L L tt . Then, studying the asymptotic behavior of solutions to (2.19), adescription of the asymptotic form of collapsing solutions very close to singularity is given in [46,Section 8.2]. In particular, [46, Proposition 8.5] shows the following behavior of parameters: a = − L t L, τ t = L − , L L tt = − b, (2.20)where b = a + a τ ≈ a obeys b τ = − N c M ν ( √ b ) ≈ − ν M e − π/ √ b (2.21)(thus, a τ (cid:28) a ) with N c = (cid:90) ∞ Q r d − dr, M = 14 (cid:90) r Q r d − dr and ν = lim r →∞ e r r d − Q ( r ) . The above characterization of parameters was first derived by Fraiman in 1985 [14] via perturba-tion arguments for the linear stability analysis around Q , and independently by Papanicolaou etal. in [28] and [26] via a solvability condition from considering dimension as a continuous variabledecreasing down to the L -critical value σ . At the first leading order, as τ → ∞ , it was shownthat a ( τ ) ≈ b / decays at the following rate a ( τ ) ≈ b / ∼ π ln τ , i.e., slower than any polynomial rate in τ , see [46, Proposition 8.6]. When more corrective termsare retained, then a ( τ ) ≈ π ln τ + c · ln ln τ , and from a τ ≈ − ν M a − e − π/a (consequence of the solvability condition) one concludes that c = 3.Furthermore, the scaling factor L ( t ) has the asymptotic form L ( t ) ≈ (cid:32) π ( T − t )ln ln( T − t ) (cid:33) . (2.22)Following [26], we further investigate numerically the correction term in L ( t ) in (2.18) andstudy the slope of a ( τ ) as a function of 1 / (ln τ + 3 ln ln τ ). In Figure 3, we observe that thefunction a ( τ ) vs 1 / (ln τ + 3 ln ln τ ) is a straight line. It would now be very tempting to makea conclusion that a ( τ ) ∼ / (ln τ + 3 ln ln τ ), as it was also predicted in d = 2 in [26] and [28].However, more care is needed here.If we change the constant c dependence in the second term of 1 / (ln τ + c ln ln τ ) and try toinvestigate the slope of a ( τ ) as a function of 1 / (ln τ + c ln ln τ ) for different c , including zero andlarge constants (we tried c = 0 , , , , TABLE BLOWUP IN NLS IN HIGH DIMENSIONS 13
Figure 3.
The curve fitting for the a ( τ ) for 4d case. This initially suggests that a ( τ ) ∼ / (ln τ + 3 ln ln τ ).at the very high focusing levels. Also, one notes that the slope of the line is not √ π as expectedfrom asymptotic analysis from [26], where it states that the correction term for a ( τ ) is given by q ( t ) ≈ (cid:32) π ln ln | T − t | (cid:33) . This is because at the time we stop our simulations (we are forced to stop as the maximalmachine’s precision is reached), the values of a ( τ ) are still far from 0, this is also observed in [26],[28], [46]. Hence, further justification for the correction term is needed and we discuss it in thenext subsection.2.2.3. Further justification of “log-log” correction.
In [9, Chapter 18], Fibich shows that the log-logregime is reached when the focusing amplitude of the solution is extremely large (e.g., a necessarycondition for the log-log to hold in dimension 2 is A (cid:29) ), which is basically impossible toobserve numerically (although see [29] for even higher order correction terms in the blow-up rate).Instead he suggests studying the reduced equations, where the equations (2.20)-(2.21) are writtenin the following form b τ ( τ ) = − ν ( b ) , a τ ( τ ) = b − a , (2.23)where ν ( b ) = c ν e − π/ √ b with c ν = 2 ν /M . Our numerical calculations give c ν ≈ . d = 2(matching § c ν ≈ .
11 in d = 4. The advantage of working with the system(2.23) as mentioned in [9] is that it can be solved by standard numerical methods over hundredsof orders of magnitudes without a significant deterioration in the numerical accuracy. We solvethis system by RK2 and then recover L ( τ ) and T − t from numerical quadratures of L ( τ ) = L (0) e − (cid:82) τ a ( s ) ds , T − t = (cid:90) ∞ τ dtdτ dτ = (cid:90) ∞ τ L ( τ ) dτ. (2.24)We show our results for d = 4 in Figure 4. On the left subfigure in Figure 4, we plot the behaviorof b ( t ) as a function of L ( t ) for three approximations: a solid blue line is the numerical solution,a dash-dot pink (straight) line is a strict adiabatic approximation b ( t ) ≡ b (0) = 10 − and thedashed green line is an asymptotic approximation of the log-log law b ( τ ) ∼ π ln τ . One can see thateven when the amplitude reaches focusing of 10 (i.e., L ( t ) ∼ − ) the log-log law is further away from the numerical solution than the strict adiabatic law, though it is slowly decreasingdown. Figure 4.
The reduced equations (2.23)-(2.24) in d = 4, b (0) = 0 . L (0) = 1, and L t (0) = 0. Left: b ( t ) as a function of L ( t ), where solid curve is a numerical solutionof the reduced system (2.23), dashdot curve is a strict adiabaticity b ≡ b (0), anddash curve is an asymptotic approximation for the log-log law. Right: the ratio of L ( t ) /L log log .On the right subfigure in Figure 4 the ratio L ( t ) /L log log is shown, where L log log = (cid:32) π ( T − t )ln ln( T − t ) (cid:33) is the anticipated blow-up rate. It is slowly growing towards 1, however, still a bit far from it.These figures also show that there is an intermediate adiabatic regime in the blow-up dynamics,which asymptotically approaches the log-log regime. The adiabatic laws are obtained in [13] and[30], and further discussed in [9].While it is challenging to observe the log-log correction term numerically, in [1] the functionalform testing was suggested and the authors succeeded in showing that among all tested func-tional forms, the log-log form stabilized the power 1 / L ( t ) ∼ (cid:18) F ( T − t ) T − t (cid:19) , (2.25)where F ( s ) = (cid:0) ln s (cid:1) γ , γ = 1 , . , . , . , ...,
0, or F ( s ) = ln ln s . Then L ( t i ) was computed at thetime t i as well as the approximation parameter ρ i = L ( t i ) L ( t i +1 ) (cid:46) ln (cid:18) F i +1 / ( T − t i +1 ) F i / ( T − t i ) (cid:19) . (2.26)Since the power ρ is expected to be (after all it is square root decay), one can check how fast theparameter ρ i converges to and what choice of F ( s ) gives the best approximation. In [1] it wasshown that F ( s ) = ln ln s provides the fastest and best parameter ρ i stabilization , moreover, it TABLE BLOWUP IN NLS IN HIGH DIMENSIONS 15 also gives the optimal quantity in several error estimates such as the standard deviation, l -normdiscrepancy, and l -norm discrepancy, which gave an extra assurance of ln ln s selection for F ( s ).In our further numerical investigation of the correction term, we use this functional approach aswell. We observe that the “log-log” correction does stabilize the power ρ close to the best, andit minimizes the approximation errors quite good (eventually) in the functional fitting. For thispurpose, we compute L ( t ) as defined in (2.25), ρ i as defined in (2.26), and (cid:15) i , the l discrepancy,defined as (following [1]) (cid:15) i = (cid:34) i − j + 1 i (cid:88) j = j (cid:18) − ρ j (cid:19) (cid:35) . (2.27)Before proceeding with our fitting results, one other remark is needed. We need to specify theprocess of calculating the quantity ( T − t i ). We take the blow-up time T to be the time when weend our simulation. we have L ( τ i +1 ) = exp(ln L ( τ i +1 )). Denoting ∆ t i +1 := t i +1 − t i , and from thelast equation of (2.1), ∆ t i +1 = δτ L ( τ i +1 ) . (2.28)Thus, the mapping for rescaled time τ back to the real time t is calculated as t ( τ i +1 ) = t (( i + 1) δτ ) := i +1 (cid:88) j =1 ∆ t j = δτ i +1 (cid:88) j =0 L ( τ j ) . (2.29)Note that as time evolves, the time difference T − t ( τ i ) will become smaller and smaller, andeventually reach saturation level (with little change), therefore we treat the stopping time t ( τ end ) = t ( τ M ) as the blow-up time T where M is the total number of iterations when reaching the stoppingcondition ( L < − ). Then, we can take T = t ( τ end ) = t ( M δτ ) = δτ M (cid:88) j =0 L ( τ j ) . (2.30)Consequently for any t i , we can calculate T − t i as T − t i = M (cid:88) j = i +1 ∆ t j = δτ M (cid:88) j = i +1 L ( τ j ) . (2.31)This indicates that instead of recording the cumulative time t i , we just need to record the elapsedtime between the two recorded data points, i.e., ∆ t i = t i +1 − t i . By doing so, it can avoid theloss of significance when adding a small number onto a larger one.Continuing with the functional fitting for the correction term, we tried several choices for F ( s ),namely, F ( s ) = 1, F ( s ) = (cid:0) ln s (cid:1) γ , γ = 0 . , . , . , . , . , . , . F ( s ) = ln ln s .In the dimension d = 4, while studying the parameter ρ i , we observe that out of all differentpowers of γ that we tried, mentioned above, in the fitting of F ( s ) = (cid:0) ln s (cid:1) γ , the powers 0 .
25 and0 . the best, see Table 3 for values of ρ i . In fact, one can firstobserve that the functional forms (cid:0) ln s (cid:1) γ decrease down to the power quite well, for example, γ = 0 .
25 gives the best approximation on steps i = 3 , , ; γ = 0 . i = 6 , , ,
9. However,the second observation is that all such forms tend to decrease down to the power and thencontinue decreasing further down (for example, the form with γ = 0 .
25 starts decreasing on steps
The fitting power ρ i from different corrections for F ( s ) i L ( t ) range 1 (ln s ) . (ln s ) . (ln s ) . (ln s ) . ln ln s e ∼ e . . . . . . e ∼ e . . . . . . e ∼ e . . . . . . e ∼ e . . . . . . e ∼ e
10 0 . . . . . . e ∼ e
11 0 . . . . . . e ∼ e
12 0 . . . . . . e ∼ e
14 0 . . . . . . e ∼ e
15 0 . . . . . . e ∼ e
16 0 . . . . . . Table 3.
Comparison of the different functional forms F ( s ) for the correction termin d = 4. “ L ( t ) range” means the values L ( t i ) ∼ L ( t i +1 ) . The blow-up rate ρ seemsto fit close to 0.5 in the intermediate regime for F ( s ) = (ln s ) γ (for example, for i = 3 , γ = 0 .
25; for i = 7 , , γ = 0 . i = 4 , , ..., i = 9 compared to the log-log form), thus,eventually not representing the appropriate correction. Note that the functional form, whichstabilizes well and stays quite close to , is the log-log form.The first observation is due to the existence of an intermediate regime in the blow-up dynam-ics, the adiabatic regime. The second observation is indicating that the adiabatic regime goesasymptotically into the log-log regime.We also computed the l discrepancy (cid:15) i , defined by (2.27), starting from the step j = 0 andaccumulating up to the step i , and show the results in Table 4. We note that in the dimension d = 4 this cumulative error is minimized the best by the log-log correction. Since the log-logregime is an asymptotic regime, we also computed the l discrepancy error starting from j = 7and up to the last reliable step i = 9, i.e., for the window of the three last approximations (seealso the discussion about the window of approximations in [1]), and recorded it in Table 5. Whileone can observe that the form F ( s ) = (cid:0) ln s (cid:1) . minimizes the discrepancy (cid:15) i the best in Table5, it is because this specific power of γ = 0 . the closest at that specificwindow. However, as discussed above, we expect it to keep decreasing, and thus, getting furtheraway from . In general, we suspect that all functional forms (cid:0) ln s (cid:1) γ will for some period of timeapproximate the power well, but then will escape away from , and thus, destabilize away fromthe blow-up regime. This is, of course, an area for further challenging numerical investigations,as well as possible testing of adiabatic regimes given by Malkin adiabatic law [30] and Fibichadiabatic law [9]. We note that the second best approximation in Table 5 is produced by thelog-log form.For dimension 5 we do a similar investigation and list the results of the functional fittings for F ( s ) = 1 , (cid:0) ln s (cid:1) . , (cid:0) ln s (cid:1) . , (cid:0) ln s (cid:1) . , (cid:0) ln s (cid:1) . , and ln ln s , in Table 6. One can observe that the TABLE BLOWUP IN NLS IN HIGH DIMENSIONS 17
The l discrepancy (cid:15) from different corrections F ( s ) i L ( t ) range 1 (ln s ) . (ln s ) . (ln s ) . (ln s ) . ln ln s e ∼ e . . e − . . . . e −
41 3 e ∼ e . . e − . . . . e −
42 7 e ∼ e . . e − . . . . e −
43 1 e ∼ e . . e − . . . . e −
44 2 e ∼ e
10 0 . . e − . e − . . . e −
45 4 e ∼ e
11 0 . . e − . e − . . . e −
46 6 e ∼ e
12 0 . . e − . e − . . . e −
47 8 e ∼ e
14 0 . . e − . e − . . . e −
48 1 e ∼ e
15 0 . . e − . e − . . . e −
49 1 e ∼ e
16 0 . . e − . e − . . . e − Table 4.
Comparison of the l discrepancy (cid:15) i in the fitting of different correctionterms in d = 4. One can see that the log-log correction minimizes this error thebest. F ( s ) 1 (ln s ) . (ln s ) . (ln s ) . (ln s ) . ln ln s (cid:15) . . e − . e − . e − . e − . e − Table 5.
Comparison of the l discrepancy (cid:15) i computed in the window from j = 7to i = 9 in the fitting of different correction terms in d = 4.forms of type (cid:0) ln s (cid:1) γ give decreasing ρ i as the step i increases; some of them reach the value 0 . γ = 0 . , .
2) and some might reach it eventually ( γ = 0 . ρ i values does not perform as good as the stabilization seenin the log-log form.We supply the l discrepancy errors (cid:15) i for the window j = 7 to i = 9 (the last 3 steps) in Table7, where the loglog fit has the smallest l deviation.We provide computations for the functional fittings for other dimensions in Appendix A.Putting all our numerical calculations, asymptotical analysis and functional fitting results to-gether, we conclude that the blow-up rate L ( t ) (with the first term correction) is given by (2.22).2.2.4. Blow-up profile.
In this subsection we investigate profiles of the blow-up solutions, andshow our results in dimensions d = 4. Figure 5 shows how the blow-up solutions v = v ( ξ, τ ) from(2.2) converge to the rescaled ground state Q , which leads to the conclusion that the profile isgiven by the rescaled (self-similar) version of the corresponding ground state L /σ Q ( rL ).We plot three different times snapshots ( τ = 2 , , τ and t variable (it is easier to distinguish and track the profiles in the rescaled τ variable, as τ → ∞ ,than in the variable t , which converges to some finite time 0 < T < ∞ , and thus, t maybeindistinguishable very close to T ). The results for other dimensions are similar, see in [57]. The fitting power ρ i from different corrections for F ( s ) i L ( t ) range 1 (ln s ) . (ln s ) . (ln s ) . (ln s ) . ln ln s e ∼ e . . . . . . e ∼ e . . . . . . e ∼ e . . . . . . e ∼ e . . . . . . e ∼ e
10 0 . . . . . . e ∼ e
11 0 . . . . . . e ∼ e
12 0 . . . . . . e ∼ e
13 0 . . . . . . e ∼ e
14 0 . . . . . . e ∼ e
16 0 . . . . . . Table 6.
Comparison of the different functional forms F ( s ) for the correctionterms in d = 5. “ L ( t ) range” means the values L ( t i ) ∼ L ( t i +1 ) . While all other formsgive a slow decrease, the loglog form seems to stabilize to the power 1 / F ( s ) 1 (ln s ) . (ln s ) . (ln s ) . (ln s ) . ln ln s (cid:15) . . e − . e − . e − . e − . e − Table 7.
The l discrepancy (cid:15) starting from j = 7 to i = 9 in the fitting ofdifferent correction terms in d = 5. This means we only consider the behaviorsclose to blow-up. One can see the log-log correction minimizes the deviation thebest at this stage.To summarize, our numerical simulations confirm that a generic blow-up in higher dimensions(in the L -critical NLS) has the log-log regime characteristics (rate and profile), and our nextgoal is to justify these observations rigorously.The analytical proof of the log-log blow-up regime (including higher dimensions) was given inseveral works of Merle and Rapha¨el, [35, 36, 37], provided the Spectral Property 1 holds. In [11],the authors were able to check it in dimensions d = 2 , , d ≥
5, while confirming the results for the lowdimensions in Appendix E. 3.
Spectral Properties
The stable “log-log” blow-up regime for the initial data with the mass slightly above the groundstate mass, M [ Q ], negative energy and zero momentum for the 1d case was proved by Merle andRapha¨el in [35] and [36], see also the work of Galina Perelman [39] in 1d. In [35] and [36] a proofof Theorem 1.3 in higher dimensions was also given, assuming the Spectral Property 1 holds true.A major obstacle for obtaining spectral properties in higher dimensions is the lack of explicit or the energy is adjusted for the non-zero momentum TABLE BLOWUP IN NLS IN HIGH DIMENSIONS 19
Figure 5.
Blow-up profile for d = 4. The blow-up solution converges to the groundstate Q . Here, the “Rescaled Q ” means L /σ Q (cid:0) rL (cid:1) .expression for the ground state Q . In 2006, Fibich, Merle and Rapha¨el made an attempt to checkthat the Spectral Property 1 holds true with a numerically-assisted proof for the dimensions upto d = 5 in [11].Before we proceed to study the Spectral Property 1 in higher dimensions ( d ≥ u ( x, t )lives in the following space: (cid:26) u is radial : (cid:90) |∇ u | + | u | e − γ | x | < ∞ (cid:27) . This implies that u r ( r ) → r → ∞ , or equivalently, u ( r ) → C as r → ∞ . On the other hand,from the analysis of the operators L and L , we know that (1 + r d − ) u ∈ L ∞ ( R d ). In d = 2, thismeans that u ∈ L ∞ ( R ), or in other words, u ( r ) → C as r → ∞ with C not necessarily beingzero. For d ≥
3, the last condition implies | u ( r ) | < C r d − , and hence, u ( ∞ ) = 0. Thus, usingthe boundary conditions u r ( L ) = 0 in d = 3 , , d = 3 ,
4. For clarification, we include Table 12 as the comparisonbetween the application of the two different boundary conditions in dimension d = 5. We alsoincluded a comparison for dimensions d = 2 , , As far as the higher dimensions d ≥
6, we think that one of the reasons that the methods from[11] can not handle d ≥ u r ( L ) = 0 for sufficiently large L (say L = 20 or L = 30), instead of u ( ∞ ) = 0. In this paper, we use the boundary condition u ( ∞ ) = 0. We also note that the same approach was used in [31] for analyzing the 3d cubicNLS equation. Our simulations show that the spectral properties holds for d ≤
10 for generalcase (not necessarily radial), and also for d = 11 and d = 12 in the radial case. We stopped ourcalculations at d = 12, since the magnitude as well as the L norm of the ground state became toolarge, and, computationally, it was not reliable to guarantee the accuracy. Moreover, the indexof both operators L and L becomes increasingly more challenging to obtain numerically as thedimension d increases beyond 12.3.1. The radial case.
In this section, we show that the Spectral Property 1 holds true from d = 5 to d = 12, and the Spectral Property 2 holds for d = 2 ...
12. We first recall the definitionsof an index of a bilinear form B , see, for example, [11] and [31]: Definition 3.1 (index of a bilinear form) . The index of a bilinear form B with respect to a vectorspace V is ind V ( B ) = min { k ∈ N | there exists a sub-space P of codimension k such that B | P is positive . } Let B and B be the bilinear forms associated with the operators L and L , see (1.10) and(1.11) in Definition (1.1). Note that the ind H ( B , ) equals to the number of negative eigenvaluesof L , . Therefore, we often refer to ind( L , ) as the number of negative eigenvalues of L , .Since the potential term V , is smooth and decays exponentially fast, according to TheoremXIII.8 in [42], which is a generalization of the Sturm Oscillation Theorem (Section XIII.7 of[42]), the operators L , have finite number of negative eigenvalues. Moreover, the number of thenegative eigenvalues can be estimated by counting the number of zeros of the solutions to thefollowing ODE: − ∂ rr U − d − r ∂ r U + V , ( r ) U = 0 ,U (0) = 1 , U r (0) = 0 . (3.1)The ODE (3.1) is a standard IVP problem, which can be solved, for example, by matlab solver“ ode45 ”. Note that when r (cid:29)
1, the equation (3.1) is essentially free (i.e., the potential termcan be neglected), see [31], and consequently, the solution must behave as U ( r ) ≈ C + C r d − . (3.2)For the L case, we apply the following statement, which is from [12]. According to thisproposition, the numerical values in Table 8 suggest that there will be no more intersections for r ≥ Proposition 3.1 (Criterion for the positivity of u , [12]) . Let u be a radial solution to − u rr ( r ) − d − r u r − V u = 0 , on r > . TABLE BLOWUP IN NLS IN HIGH DIMENSIONS 21
Let V + = max { V, } . Assume that there holds for some r > ∂ r u ( r ) u ( r ) > and • for d ≥ , ∀ r ≥ r , V + ( r ) ≤ ( d − r ; • for d = 2 , ∀ r ≥ r , V + ( r ) ≤ r (log r ) and u (cid:48) ( r ) u ( r ) ≥ r (cid:90) ∞ r V + ( r ) rdr. Then u cannot vanish for r ≥ r (see Table 8). d ∂ r u ( r ) u ( r ) 1 . e −
05 4 . e −
05 1 . e −
05 6 . e − − V ( r ) − d − r − . − . − . − . d ∂ r u ( r ) u ( r
0) 1 . e −
06 4 . e −
07 1 . e −
07 6 . e − − V ( r ) − d − r − . − . − . − . Table 8.
Values of the quantities from Proposition 3.1 at r = 6.For the case of L , the theorem is not applicable. We adopt the argument from [31]: we firstnotice that the equation (3.1) converges to the free equation (3.2). Then, we choose a largeenough interval (say L = 100) to ensure that the solution goes to a constant. The constant canbe found from C r d − i + C = U i , C r d − i − + C = U i − , (3.3)where r i and U i are discretized points of r and U ( r ). Once the constant C starts to stabilize, weconclude that the solution enters the free region and no more “zeros” will occur.We also point out that one needs to be careful in the numerical calculation of the potential V or V , since when d ≥
5, the term V = 2 d Q d − rQ r generates the negative power of Q . This may fail to describe the exponential decay property,especially, when d ≥
8. An alternative way is needed to calculate the potential V and V . Weprovide a new approach for that and discuss the details in the Appendix C.Our numerical solutions of the equation (3.1) are given in Figure 6 as an example for the case d = 5, there U stands for the solution to L and Z for L . Solutions to other dimensions of theequation (3.1) are similar. We conclude the following statement. Figure 6.
Solutions of (3.1) in d = 5. Numerical justification of Proposition 3.2.The blue line shows the behavior of U from L U = 0 and the green line shows thebehavior of Z from L Z = 0. We also show the behavior of the tail of Z , and theconstant C it approaches to as r → ∞ , see (3.3). Proposition 3.2 (indices of B , ) . For d = 5 to d = 12 , the indices of L , in the radial case are ind( L ) = 2 , ind( L ) = 1 Therefore, ind H r ( B ) = 2 , ind H r ( B ) = 1The following property shows that the indices of the bilinear forms are stable under the pertur-bations. Thus, it is sufficient to check the terms B , ( u, u ) > B , ( u, u ) > δ (cid:82) | u | e −| x | for some sufficiently small δ , given in Definition 1. Proposition 3.3 ([11], [31]) . For the operators L , (from the Spectral Property), there exists auniversal constant δ > , sufficiently small, such that for the perturbed operators ¯ L , = L , − δ e −| x | , the associated bilinear forms are stable, i.e., ind H ( ¯ B , ) = ind H ( B , ) . We return to the discussion of the proof of the spectral property, which involves solving theBVP problem L , U = f . While the numerical calculations suggest that L , are invertible, theproof of the invertibility of the operators L , in [11] and [31] works in straightforward adaptationto our cases. Proposition 3.4 (Invertibility of L , ) . Let d ∈ { , , · · · , } and f ∈ C loc ( R d ) be radiallysymmetric with | f ( r ) | ≤ e − Cr . Then there exists a unique radial solution to L , u = f with (1 + r d − ) u ∈ L ∞ . The following definition lists the numerical values of the bilinear forms which we need in theproof of the Spectral Property. It involves the computation of the BVP problem L , U = f . Wetake u r (0) = 0 as the left boundary condition, since u ( r ) is radially symmetric. We construct TABLE BLOWUP IN NLS IN HIGH DIMENSIONS 23 the artificial boundary condition u ( L ) + Ld − u r ( L ) = 0 to approximate the boundary condition u ( ∞ ) = 0 (see details in Appendix B, also the reader can refer to [31]). Definition 3.2 (numerical representation of the bilinear form) . Let the operator L i u = f , i = 1 , solve the linear BVP (cid:40) L , u = f, i = 1 , u r (0) = 0 , u ( ∞ ) = 0 . (3.4) Define L U = Q, L U = Q ; (3.5) L Z = Q , L Z = Q ; (3.6) and denote the following constants as values of the bilinear forms K = B ( U , U ) , K = B ( U , U ) , K = B ( U , U ) , K = B ( U , U ); (3.7) J = B ( Z , Z ) , J = B ( Z , Z ) , J = B ( Z , Z ) , J = B ( Z , Z ) . (3.8) We also define the determinants of matrices K and J by KK = K K − K K and J J = J J − J J . (3.9)We list the values of K ij and J ij from Definition 3.2 for dimensions 5 ≤ d ≤
12 in Table 9 andTable 10, respectively. d K K K KK | K K − | − . . − . . e − − . . − . . e − − . . − . .
411 8 e − − . . − . e − − . . − . e − − . . − .
457 24415176279909 1 e − − . . − . e
16 5 e − − . . − e
18 1 e − Table 9.
Values of the bilinear form B from (3.7) and (3.9) via K i ’s.We use two methods to solve the equation (3.4): one is the Chebyshev collocation method; theother method is the matlab solver “ bvp4c ”. These two methods lead to basically the same results,for a comparison in dimension d = 5 see Table 11 (values of K i ’s) and Table 12 (values of J i ’s).Since L and L are self-adjoint operators, the difference | K K − | , and the corresponding one for J ’s, is one way to check the numerical consistency, we list those values in the last columns of theTables 9-10; in Tables 11-12 we list the differences | K − K | and | J − J | , correspondingly.We note that we use u ( ∞ ) = 0 boundary condition when computing the values of K ’s in Table9 as well as in Table 11; moreover, in Table 11 we provide results obtained by two methods forcomparison purposes. d J J J J J | J J − | . − . . − . e −
106 611 . − . . − . e −
107 5608 . − . . − . e −
98 60626 . − . . − . e −
99 758310 . − . . − . e −
910 10852351 . − . . − e −
811 176804567 . − . . − e
17 1 e −
712 3286852523 . − . . − e
20 1 e − Table 10.
Evaluation of the bilinear form B from (3.8) and (3.9) via J i ’s. d = 5 K K K KK | K − K | “cheby” − . . − . . e − bvp4c ” − . . − . . e − Table 11.
Comparison of values of K i ’s in d = 5 between Chebyshev-collocationmethod and “ bvp4c ”. The boundary condition used here in both methods is u ( ∞ ) = 0.In the Table 12 we show two different boundary conditions: u ( ∞ ) = 0 (first row with Chebyshevcollocation method) and u (cid:48) ( L ) = 0 with Chebyshev collocation methods (second row) and with“ bvp4c ” from matlab (third row). This is for comparison of the results with lower dimensions,since it makes a difference in dimension 4 (though it does not influence the signs, thus, theconclusion of the Spectral Property), and it completely changes the results in dimension 5 (andhigher). Therefore, starting from the dimension 5 and higher, we only use boundary condition asin the first row of Table 11.B.C. J J J J J | J − J | u ( ∞ ) = 0 by “cheby” 80 . . − . − . e − u (cid:48) ( L ) = 0 by “cheby” − . − .
538 552 .
837 165185 .
398 2 e − u (cid:48) ( L ) = 0 by “ bvp4c ” − . − . . . e − Table 12.
Comparison of values of J i ’s for d = 5. Note that the boundary condi-tions affect the final results.With these bilinear forms calculated, we reach the following proposition: Proposition 3.5.
The bilinear form B ( f, f ) is coercive on the space U = { Q, Q } ⊥ and U ⊆ H r ,where H r stands for the radial functions in H . The bilinear form B ( g, g ) is coercive on the space V = { Q , Q } ⊥ and V ⊆ H r . Therefore, the Spectral Property 1 holds in the functional space U × V ⊆ H r × H r . TABLE BLOWUP IN NLS IN HIGH DIMENSIONS 25
Proof.
We outline the key steps of the proof, we refer to [11], [31] or [44] for the details, as theyare the same. Let’s consider the form B and recall from Proposition 3.2 that the ind H r ( B ) = 2.From the Table 9 we have K = B ( U , U ) < . This suggests that U is one of the negative spansof L associated with B . Similarly, we have U is the other negative span of L associated to B ,since K <
0. Moreover, the determinant
KK > B ) = 2, the remaining spans, that are orthogonalto ( U , U ) in the sense of B , or equivalently, orthogonal to Q or Q in L sense, must generatethe positive outcome, i.e., B ( f, f ) > (cid:104) f, Q (cid:105) = (cid:104) f, Q (cid:105) = 0 . Next, we consider the bilinear form B , while J = B ( Z , Z ) and J = B ( Z , Z ) generatethe positive values, we check their linear combination ¯ Z = Z + αZ with α = − J J . This valuecomes from the fact that if we calculate the quadratic form, this constant gives the minimumvalue for the bilinear form. Moreover, B ( ¯ Z, ¯ Z ) = C · J J for some constant
C >
0, i.e., the valueof the bilinear form has the same sign as the determinant. Thus,
J J < B ( ¯ Z, ¯ Z ) < Z lies on the negative span of L as we desire. Therefore, any g , which is orthogonal to thedirection ¯ Z , i.e., B ( g, ¯ Z ) = 0, or equivalently, (cid:104) g, Q (cid:105) = (cid:104) g, Q (cid:105) = 0, is orthogonal to ( Z , Z ),since ¯ Z is their linear combination, and consequently, we have B ( g, g ) >
0, which justifies thestatement of Proposition 3.5.To make the argument rigorous, note that the functions f and g , considered above, are notin H r (from the Proposition 3.4, we have (1 + r d − ) f ∈ L ∞ ( R d ) and a similar requirement for g ). We introduce an appropriate cut-off function and then take the limit. For further details, see[11], [31] or [44]. (cid:3) Finally, we provide details on the Spectral Property 2, where the orthogonal conditions for thesecond bilinear form can be changed to just one condition. We set L Z = Q and compute thequantity (cid:104) L Z, Z (cid:105) . Table 13 contains the values for the quantity (cid:104) L Z, Z (cid:105) in dimensions 2 to 6with two different boundary conditions, and then in dimensions 7 to 12 with u ( ∞ ) = 0.Boundary conditions d u (cid:48) ( L ) = 0 (cid:104) L Z, Z (cid:105) − . − . . − . − . u ( ∞ ) = 0 (cid:104) L Z, Z (cid:105) NA − . − . − . − . d (cid:104) L Z, Z (cid:105) , u ( ∞ ) = 0 − . − . − − − − Table 13.
Value of (cid:104) L Z, Z (cid:105) for different dimensions, where L Z = Q . Remark 3.3.
Since by Proposition 3.2 the index of L is one, arguing as in Proposition 3.5 andincorporating the results from Table 13, we obtain that the spectral property holds for L only withone condition: (cid:104) Q, g (cid:105) = 0 , which proves Theorem 1.2, and thus, validates the Spectral Property 2in the radial case. For completeness, we include details in the general case in Theorem 3.6.
The non-radial case.
To study the non-radial case, we rewrite our operators in the formof spherical harmonics, i.e., L ( k )1 , = − ∂ rr − d − r ∂ r + V , ( r ) + k ( k + d − r , k = 0 , , , .... (3.10) The notation L ( k )1 , and B ( k )1 , will stand for the k th spherical harmonics. When k = 0, it is simplythe radial case, which we already discussed in Section 3.1. For k > B ( k )1 , havethe same properties as L ( k )1 , , similar to the case k = 0 (see [42], [11] and [31]). We first studyindices of the k th bilinear forms. Proposition 3.6 (indices in the non-radial case) . The index of the bilinear form ind( B ( k )1 , ) equalsthe number of zeros of the IVP problem, counting from r > : L ( k )1 , U ( k ) = 0 , lim r → U ( k ) ( r ) r k = 1 , lim r → ddr U ( k ) ( r ) r k = 0 . (3.11) Numerical calculations show that for d = 4 , , , , , , , we have ind( L (1)1 ) = 1 , ind( L (2)1 ) = 0 , ind( L (1)2 ) = 0 . Consequently, ind H ( B (1)1 ) = 1 , ind H ( B (2)1 ) = 0 , ind H ( B (1)2 ) = 0 . For d = 11 , , we get ind( L (1)1 ) = 1 , ind( L (2)1 ) = 1 , ind( L (3)1 ) = 0 , ind( L (1)2 ) = 1 , ind( L (2)2 ) = 0 , and consequently, ind H ( B (1)1 ) = 1 , ind H ( B (2)1 ) = 1 , ind H ( B (3)1 ) = 0 . ind H ( B (1)2 ) = 1 , ind H ( B (2)2 ) = 0 . This proposition is obtained from the numerical solutions of the IVP in (3.11) in the corre-sponding cases, see Figures 7 as an example in the case d = 5. Figures 8 and 9 show the solutionof (3.11) in cases d = 11 ,
12 and k = 2 ,
3, since they are different from the cases d = 5 , · · · , U stands for the solution of L ( k )1 , and Z for L ( k )2 .In this non-radial case, in our numerical simulations, we want to get rid of the limit termsin the boundary conditions. We use the approach from [31]: let U ( k ) ( r ) = r k ˜ U ( k ) ( r ), then theoperator L ( k ) i becomes ˜ L ( k )1 , = ∂ rr − d − kr ∂ r + V , . (3.12)We rewrite (3.11) as follows r k (cid:18) − ∂ rr ˜ U − d − kr ∂ r ˜ U + V , ( r ) ˜ U (cid:19) = 0˜ U (0) = 1 , ˜ U r (0) = 0 . (3.13)The IVP (3.13) can be solved by matlab solver “ ode45 ”. Then the solution U can be re-constructed by U ( r ) = r k ˜ U ( r ). From [31, Sec. 3], it follows that U ( r ) satisfies the asymptoticbehavior for r (cid:29) U ≈ C r k + C r − d − k . (3.14) TABLE BLOWUP IN NLS IN HIGH DIMENSIONS 27
Figure 7.
Solutions of (3.11) in d = 5. Numerical justification of the Proposition3.6. Left figure is the for k = 1 and right figure is for k = 2. The blue line showsthe behaviors of U from L ( k )1 U = 0 and the red line shows the behaviors of Z from L ( k )2 Z = 0. To the right of each plot are zooms of the tails for U and Z . Onecan see the tails of U and Z increase or decrease with a rate r k , which justifies theasymptotic behavior in (3.14). Figure 8.
Solutions of (3.11) in d = 11 (left) and d = 12 (right), or numericaljustification of the Proposition 3.6. The blue line shows the behaviors of U from L (2)1 U = 0 and the red line shows the behaviors of Z from L (2)2 Z = 0. To the rightof each plot are zooms of the tails for U and Z . One can see U and Z increase ordecrease with a quadratic rate. This justifies the asymptotic behavior in (3.14) for k = 2.This indicates that for r (cid:29) U ( r ) either grows or decays with a polynomial rate r k . Conse-quently, no more zeros will occur. We stop calculations once we see such polynomial increase ordecrease.The following property ([42], also see [11], [31]) shows that once we found some k such thatind( L ( k )1 , ) = 0, we can stop the calculation. This avoids checking infinitely many k ’s. Figure 9.
Solutions of (3.11) in d = 11 (left) and d = 12 (right), or numericaljustification of the Proposition 3.6. The blue line shows the behaviors of U from L (3)1 U = 0 and the red line shows the behaviors of Z from L (3)2 Z = 0. To the rightof each plot are zooms of the tails for U and Z . One can see U and Z increase ordecrease with a cubic rate. This justifies the asymptotic behavior in (3.14) for k = 3. Proposition 3.7.
The index is monotonic with respect to k , that is, ind( B ( k +1)1 , ) ≤ ind( B ( k )1 , ) . Moreover, the uniqueness and stability of the indices of the bilinear forms B , can also beextended to the non-radial case: Proposition 3.8 (Stability, [11], [31]) . For the operators L ( k )1 , defined in (3.10), there exists auniversal constant δ > , sufficiently small, such that for the perturbed operators ¯ L ( k )1 , = L ( k )1 , − δ e −| x | , the associated bilinear forms are stable: ind H ( ¯ B ( k )1 , ) = ind H ( B ( k )1 , ) . Proposition 3.9 (Invertibility, [11], [31]) . Let d ∈ [5 , and f ∈ C loc ( R d ) with radial symmetryand | f ( r ) | ≤ e − Cr . Then there exists a unique radial solution to L ( k )1 , u = f with (1 + r d − k ) u ∈ L ∞ ( R d ) . Similar to the radial case, we list the numerical values of the bilinear forms, needed to provethe spectral properties for the non-radial case.
Definition 3.4 (numerical computations for the bilinear form) . Let the operators L ( k ) i , i = 1 , ,solve the linear BVP (cid:40) L ( k )1 , u = f,u r (0) = 0 , u ( ∞ ) = 0 . (3.15) TABLE BLOWUP IN NLS IN HIGH DIMENSIONS 29
Define L (1)1 U (1)1 = rQ, L (1)2 Z (1)1 = Q r ; (3.16) L (2)1 U (2)1 = rQ. (3.17) and denote the following constants as values of the bilinear forms K (1)11 = B (1)1 ( U (1)1 , U (1)1 ) , J (1)11 = B (1)2 ( Z (1)1 , Z (1)1 ) , K (2)11 = B (2)1 ( U (2)1 , U (2)1 ) . (3.18)The values of K ( k ) ij ’s and J ( k ) ij ’s are listed in Tables 14, 15 and 16. d K (1)11 − . − . − . − . d K (1)11 − . − − . − . Table 14.
Values of K (1)11 corresponding to the form B (1)1 . d J (1)11 . . . . d J (1)11 . − − Table 15.
Values of J (1)11 corresponding to the form B (1)2 . d
11 12 K (2)11 Table 16.
Values of K (2)11 corresponding to the form B (2)1 .The values of these bilinear forms are computed by Chebyshev collocation method with N =1025 collocation points on the interval L = 100. The artificial boundary condition is constructedin a similar way as in the previous Section 3.1 to approximate u ( ∞ ) = 0, see also Appendix Cfor details on artificial boundary condition.We are now ready to establish both spectral properties. Theorem 3.5.
Let the space U = { Q, Q , x i Q } ⊥ ⊆ H , and the space V = { Q , Q , Q x i } ⊥ ⊆ H ,where i = 1 , , · · · , d . In the dimensions d ≤ , the Spectral Property 1 holds on the givensubspace U × V ⊆ H × H . In the dimension d = 11 or d = 12 , the Spectral Property 1 holdsonly on the given subspace U × V ⊆ H r × H r . It is indecisive in the non-radial case.Proof. We only outline the main idea, as it follows the proof in [11] and [31]. For d ≤ L (1)1 ) = 1 , and we see K < U (1)1 generates the negative span. Forany f ∈ { U (1)1 } ⊥ in the B (1)1 sense, B (1)1 ( f, f ) >
0. For L (2)1 and L (1)2 , ind( L (2)1 ) = ind( L (1)2 ) = 0,thus B (2)1 and B (1)2 are coercive. The coercivity of both B and B implies the spectral property in the given space U × V . We need to note that f and g may not necessarily in H r , since fromthe Proposition 3.9, (1 + r d − k ) f ∈ L ∞ ( R d ) and the same for g . Similar to the radial case wediscuss above, we need to introduce an appropriate cut-off function and then take the limit aswhat the authors did in [11], [31] or [44].For the case d = 11 or d = 12. B (1)1 and B (1)2 are coercive on the subspace of H r that satisfiesthe orthogonal conditions from K < J <
0. However, when k = 2, the index of L (2)1 is still one and we get the positive values from the bilinear forms (see Table 16). Thus, thecoercivity of B (2)1 becomes indecisive. Hence, with current computations, we can only show thespectral property under the radial assumption for d = 11 or d = 12. (cid:3) Theorem 3.6.
Let the space U = { Q, Q , x i Q } ⊥ ⊆ H , and the space V = { Q, Q x i } ⊥ ⊆ H ,where i = 1 , , · · · , d . In the dimensions d ≤ , the Spectral Property 2 holds on the givensubspace U × V ⊆ H × H . In the dimension d = 11 or d = 12 , the Spectral Property 2 holdsonly on the given subspace U × V ⊆ H r × H r . It is indecisive for the non-radial case.Proof. The argument follows the proof of Theorem 3.5. From Table 13, one can see that the spanof ˜ Z , where L ˜ Z = Q , is negative. Combining with the fact that Ind( L (0)2 ) = 1, we conclude that˜ Z is the only negative span we have. Therefore, the coercivity of B (0)2 is reached on the subspaceorthogonal to Q . For the non-radial case, the proof is the same as the orthogonal conditionsto extend to the nonradial setting are the same as in the Theorem 3.5, and thus, finishing theproof. (cid:3) Remark 3.7.
In the dimensions d = 11 and d = 12 , ind( L (1)2 ) = 1 instead of the zero index, whichis not what we obtained in the lower dimensions. We double checked this with the standard 4thorder explicit Runge-Kutta method (RK4), taking the step size h = 0 . in obtaining solutionsto (3.11) . This led to the same results as in the above calculations via the matlab solver “ ode45 ”.Table 15 shows the negative values for J (1)11 in d = 11 and d = 12 , which also suggests thatthe index of L (1)2 is not zero in those dimensions. We also numerically calculated the negativeeigenvalue of L (1)2 : λ = − . in the dimension d = 11 , and λ = − . in the dimension d = 12 . This can be compared with the situation of the index of L (0)1 in the dimensions d = 2 and d = 3 , ind( L (0)1 ) = 1 when d = 2 but ind( L (0)1 ) = 2 when d = 3 , see [11] and [12] . Conclusions
In this paper we first discussed direct numerical simulations of the generic blow-up solutionsin the L -critical NLS equation in higher dimensions ( d = 4 , ...,
12) under the radial symmetryassumption. Our results show that the “log-log” law is universal for all L -critical NLS equations(at least up to d = 12). Secondly, we investigated the Spectral Property 1 in higher dimensions,which is the essential part of the analytical proof of the “log-log” regime for the cases when themass of the negative energy initial data is slightly above the mass of Q , the corresponding groundstate Q for the given dimension and nonlinearity. We confirm that the Spectral Property 1 (aswell as a modified version of it) holds from d = 5 to d = 10 in a general case, and for d = 11 and d = 12, at least, in the radial case. Therefore, we conclude that the “log-log” blow-up regime isthe stable blow-up regime in d ≤
10 and radially stable in d ≤ TABLE BLOWUP IN NLS IN HIGH DIMENSIONS 31
Appendix A: Justification of “log-log” corrections in d = 6 , ..., F ( s ) for the correction term in the blow-up rate (2.25)in the dimensions d = 6 , ...,
12. For dimensions d = 6 and 7 we list our computations for ρ i with F ( s ) = 1 , (cid:0) ln s (cid:1) γ , γ = 0 . , . , .
15 and ln ln s in Tables 17 and 18. One can observeThe fitting power ρ i from different corrections F ( s ) i L ( t ) range 1 (ln s ) . (ln s ) . (ln s ) . ln ln s e ∼ e . . . . . e ∼ e . . . . . e ∼ e . . . . . e ∼ e . . . . . e ∼ e . . . . . e ∼ e
11 0 . . . . . e ∼ e
12 0 . . . . . e ∼ e
13 0 . . . . . e ∼ e
14 0 . . . . . e ∼ e
15 0 . . . . . Table 17.
Comparison of the different functional forms F ( s ) for the correctionterm in d = 6. “ L ( t ) range” means the values L ( t i ) ∼ L ( t i +1 ) .The fitting power ρ i from different corrections F ( s ) i L ( t ) range 1 (ln s ) . (ln s ) . (ln s ) . ln ln s e ∼ e . . . . . e ∼ e . . . . . e ∼ e . . . . . e ∼ e . . . . . e ∼ e . . . . . e ∼ e
11 0 . . . . . e ∼ e
12 0 . . . . . e ∼ e
13 0 . . . . . e ∼ e
14 0 . . . . . e ∼ e
15 0 . . . . . Table 18.
Comparison of the different functional forms F ( s ) for the correctionterm in d = 7. “ L ( t ) range” means the values L ( t i ) ∼ L ( t i +1 ) .that log-log shows better stabilization in the approximation of the power ρ i , though, the otherapproximation seem to reach 0 . We also give l discrepancy (cid:15) i results (for the last three steps from j = 7 to i = 9) in Table 19for both d = 6 ,
7. It seems that the smallest error is given by the power γ = 0 .
2, which can beexplained similarly as in the dimension d = 4 (during the studied time interval, this functionalapproximation decreases down to 0 .
5, but then it will continue decreasing). Similar commentscan be made about other powers of γ . d F ( s ) = 1 F ( s ) = (ln s ) . F ( s ) = (ln s ) . F ( s ) = (ln s ) . F ( s ) = ln ln s . . e − . e − . e − . e −
47 0 . . e − . e − . e − . e − Table 19.
The l discrepancy (cid:15) i calculated from j = 7 to i = 9 for the fittingsgiven in Tables 17 - 18 in d = 6 and d = 7.For higher dimensions, for brevity, we list F ( s ) = 1 and F ( s ) = ln ln s results, see Tables 20,21, 22, 23, 24 for dimensions d = 8 , , , ,
12, correspondingly. i L ( t ) range ρ i from F ( s ) = 1 ρ i from F ( s ) = ln ln s e ∼ e . . e ∼ e . . e ∼ e . . e ∼ e . . e ∼ e
11 0 . . e ∼ e
12 0 . . e ∼ e
13 0 . . e ∼ e
14 0 . . e ∼ e
16 0 . . Table 20.
Comparison of the functional fitting in d = 8. i L ( t ) range ρ i from F ( s ) = 1 ρ i from F ( s ) = ln ln s e ∼ e . . e ∼ e . . e ∼ e . . e ∼ e . . e ∼ e
11 0 . . e ∼ e
12 0 . . e ∼ e
13 0 . . e ∼ e
14 0 . . e ∼ e
16 0 . . Table 21.
Comparison of the functional fitting in d = 9. TABLE BLOWUP IN NLS IN HIGH DIMENSIONS 33 i L ( t ) range ρ i from F ( s ) = 1 ρ i from F ( s ) = ln ln s e ∼ e . . e ∼ e . . e ∼ e . . e ∼ e . . e ∼ e
11 0 . . e ∼ e
12 0 . . e ∼ e
13 0 . . e ∼ e
14 0 . . e ∼ e
16 0 . . Table 22.
Comparison of the functional fitting in d = 10. i L ( t ) range ρ i from F ( s ) = 1 ρ i from F ( s ) = ln ln s e ∼ e . . e ∼ e . . e ∼ e . . e ∼ e . . e ∼ e
10 0 . . e ∼ e
12 0 . . e ∼ e
13 0 . . e ∼ e
14 0 . . e ∼ e
16 0 . . Table 23.
Comparison of the functional fitting in d = 11. i L ( t ) range ρ i from F ( s ) = 1 ρ i from F ( s ) = ln ln s e ∼ e . . e ∼ e . . e ∼ e . . e ∼ e . . e ∼ e
10 0 . . e ∼ e
12 0 . . e ∼ e
13 0 . . e ∼ e
15 0 . . e ∼ e
16 0 . . Table 24.
Comparison of the functional fitting in d = 12.We conclude this section with providing the l discrepancy error (computed for the last threesteps) for dimensions d = 8 , ...,
12 in Table 25. The error (cid:15) i stays on the same order ∼ − as inall previous dimensions. d F ( s ) = 1 F ( s ) = ln ln s . . e −
49 0 . . e −
410 0 . . e −
411 0 . . e −
412 0 . . e − Table 25.
The l discrepancy (cid:15) i calculated from j = 7 to i = 9 for the fittingsgiven in Tables 20 - 24 in d = 8 , , , , Appendix B: Artificial boundary conditions
The BVP problem (3.4) and (3.15) involves the boundary condition u ( ∞ ) = 0. We use theapproach in [31] to construct the artificial boundary condition to approximate u ( ∞ ) = 0.As V , and f decay exponentially fast, as it was discussed before, for r (cid:29) − ∂ rr u − d − kr ∂ r u = 0 , r (cid:29) . (4.1)This equation has the explicit solution u = C r d − k + C := C φ + C φ , (4.2)where φ = r d − k and φ = 1. Since u ( ∞ ) = 0, only span φ remains. Consequently, at r = L ,the solution to (3.15) must be linearly dependent on φ . Therefore, their Wronskian should equalto 0, i.e., det (cid:18) u ( L ) φ ( L ) u r ( L ) φ (cid:48) ( L ) (cid:19) = 0 , (4.3)which is u ( L ) + Ld − k u r ( L ) = 0 . (4.4)Therefore, the boundary condition u ( ∞ ) = 0 is substituted by u ( L ) + Ld − k u r ( L ) = 0 at thelength of the interval r = L . Appendix C: Computation of the potentials V and V When d ≥
5, the term Q d − in V , has the negative power. We rewrite it as follows, forexample for V V = 2 d Q d r Q r Q .
Recall that the ground state Q decays with the rate Q ∼ r − d − e − r for r (cid:29) Q r ∼ − Q and only differs by a polynomial power, and V decays with anexponential rate V ∼ e − d r for r (cid:29)
1. However, if we compute the V from (1.08), instead ofkeep decreasing, the potential term V becomes oscillating at the magnitude 10 − when r ≥ TABLE BLOWUP IN NLS IN HIGH DIMENSIONS 35 in d = 8. This fails to describe the “fast decay” property of the potential terms and may causetrouble in the process of finding the index of L ( k )1 , (according to [42], L , has infinitely manynegative eigenvalues if the potential decays slower than r ). In fact, this issue comes from thenumerical calculation of the ground state Q . When we compute the ground state Q , the Q itselfstops decreasing when it reaches the magnitude of 10 − , since it is less than the tolerance weset for the fixed point iteration and it approaches the machine accuracy. Nevertheless, (10 − ) d could be a relatively large number, especially when the dimension d is large, say d ≥ Q asit is below our tolerance and approaching the machine accuracy, we come up with an alternativeway to compute the potential terms, say V . We outline it next.Consider the function P ( r ) = Q ( r ) e r , (4.5)Where P decays with a polynomial rate r − d − , since Q ∼ r − d − e − r for r (cid:29)
1. Then, a straight-forward calculation gives us V = 2 d P d r (cid:18) P r P − (cid:19) e − d . (4.6)This implies that if we find the profile of P , which decays polynomially, then P can not be toosmall for r (cid:29)
1. Moreover, the last term e − d describes an exponential decay for the potential V .Putting the expression P ( r ) = Q ( r ) e r into the ground state equation (1.5), we have P rr − P r + d − r ( P r − P ) + P d +1 e − d r = 0 , ( P > ,P r (0) − P (0) = 0 ,P ( ∞ ) = 0 . (4.7)We can construct the artificial boundary condition to approximate P ( ∞ ) = 0 the same way asfor L , u = f . As described in Section 3, the linear flow of the equation (4.7) gives us P rr − P r + d − r ( P r − P ) = 0 . (4.8)The solutions of the free equation decay with the leading order P ∼ r − d − . Therefore, we canconstruct the artificial boundary condition at r = L :det (cid:32) P ( r ) r − d − P r ( r ) − d − r − d +12 (cid:33) = 0 , (4.9)i.e., P + 2 Ld − P r = 0 . (4.10)We use the Chebyshev differential matrices to discretize the first and second derivatives P r and P rr , then the discretized equation (4.7) changes to the nonlinear system MP − f ( P ) = , (4.11)with the first and last rows substituted by the prescribed boundary conditions. This system (4.11)can be solved, for example, by the matlab solver “ fsolve ”. d 5 6 7 8 9 10 11 12 (cid:107) Q − P e − r (cid:107) ∞ e −
09 5 e −
09 7 e −
09 1 e − e −
08 4 e − e −
05 6 e − (cid:107) V − ˜ V (cid:107) ∞ e −
09 2 e −
07 8 e −
06 1 e − e − e − . . Table 26.
Comparison of the quantities obtained in two different ways, where˜ V = d P d r (cid:0) P r P − (cid:1) e − d . We set L = 100, N = 1024.Figure 4 shows the comparison of the profiles for V when d = 8 and r (cid:29)
1. Table 26 showsthe difference of the ground state Q and P e − r , as well as the difference of V calculated from Q and P , denoted by V and ˜ V correspondingly. Note how this method allows us to describe the“exponential” decay and completely avoid oscillations! Figure 10.
Upper: the value of V calculated directly from Q ; lower: V obtainedfrom P . Our method calculating V P completely avoids oscillations. Appendix D: Comparison with prior results on the spectral property in d = 2 , , d = 2 , , u (cid:48) ( L ) = 0 and u ( ∞ ) = 0, generate the different outcomes, this does not affect the proofof the Spectral Property 1 in those dimensions, as the signs of the quantities under considerationremain the same, see Tables 27, 28 and 29. We do see some differences between our results and TABLE BLOWUP IN NLS IN HIGH DIMENSIONS 37 the results in [11], even when using the same boundary condition (however, this does not influencethe outcome for the Spectral Property 1). Note that the situation is different in dimension 5 (andhigher) as we discussed in the paper.Boundary conditions K K K J J J In [11, pp7-8], u (cid:48) ( L ) = 0 − . N A N A . . − . u (cid:48) ( L ) = 0 from “cheby” − . . − . . . − . u (cid:48) ( L ) = 0 from “ bvp4c ” − . . − . . . − . Table 27.
Comparison of the bilinear forms for K ij ’s and J ij ’s by using differentboundary conditions. This is for the case d = 2.Boundary conditions K K K J J J In [11, pp7-8], u (cid:48) ( L ) = 0 − . N A N A . . − . u (cid:48) ( L ) = 0 from “cheby” − . − .
058 1 . . . − . u (cid:48) ( L ) = 0 from “ bvp4c ” − . − .
058 1 . . . − . u ( ∞ ) = 0 − . − . . . . − . Table 28.
Comparison of the bilinear forms for K ij ’s and J ij ’s by using differentboundary conditions. This is for the case d = 3.Boundary conditions K K K J J J In [11, pp7-8], u (cid:48) ( L ) = 0 − . − . . − u (cid:48) ( L ) = 0 from “cheby” − . − . . . . − . u (cid:48) ( L ) = 0 from “ bvp4c ” − . − . . . . − . u ( ∞ ) = 0 − . − . . . . − . Table 29.
Comparison of the bilinear forms for K ij ’s and J ij ’s by using differentboundary conditions. This is for the case d = 4.To further justify our numerical results, we also studied the ( L -supercritical) 3d cubic NLSequation ( d = 3 and p = 3), which was discussed by Marzuola and Simpson in [31]. From Table30, one can see that our numerical results match the results of Marzuola and Simpson [31] verywell. 3d cubic case K J Numerical results in [31, pp14] 1 . − . . − . Table 30.
Comparison of the bilinear forms for K ’s and J ’s with the resultsfrom [31] obtained by “ bvp4c ”, here, L Z = Q + rQ r and J = (cid:104) L Z , Z (cid:105) . Thisis for the 3d cubic NLS equation ( L -supercritical). References [1] G. D. Akrivis, V. A. Dougalis, O. A. Karakashian and W. R. McKinney,
Numerical approximation of blow-upof radially symmetric solutions of the nonlinear Schr¨odinger equation , SIAM J. Sci. Comput. 25 (2003), no.1, 186-212.[2] G. D. Akrivis, V. A. Dougalis, O. A. Karakashian and W. R. McKinney,
Galerkin-finite element methodsfor the nonlinear Schr¨odinger equation , in “Advances on Computer Mathematics and its Applications” (E.Lipitakis ed.), World Scientific, Singapore 1993, 85-106.[3] M. Birem and C. Klein,
Multidomain spectral method for Schr¨odinger equations , Adv. Comput. Math. (2016)42, 395–423.[4] H. Berestycki and P. Lions,
Nonlinear scalar field equations. I. Existence of a ground state , Arch. Ration.Mech. Anal. 82(4)(1983), 313–345.[5] V. Budneva, V. Synakh, V. Zakharov,
Certain modes for wave collapse , Sov. J. Plasma Phys. 1. 1975, 335–338.[6] V. S. Buslaev and G. S. Perelman,
Nonlinear scattering: the states which are close to a soliton , Zap. Nauchn.Sem. POMI, 1992, vol. 200, 38–50.[7] T. Cazenave,
Semilinear Schr¨odinger equations , American Mathematical Soc, 2003.[8] S. Dyachenko, A. Newell, V. Zakharov,
Optical turbulence: weak turbulence, condensates and collapsingfilaments in the nonlinear Schr¨odinger equation , Physica D 57, 1992, 96–160.[9] G. Fibich,
The nonlinear Schr¨odinger equation, singular solutions and optical collapse , Springer, 2015.[10] G. Fibich, N. Gavish, X. Wang,
New singular solutions of the nonlinear Schr¨odinger equation , Physica D, 211(2005), 193–220.[11] G. Fibich, F. Merle, P. Rapha¨el,
Proof of a spectral property related to the singularity formation for the L critical nonlinear Schr¨odinger equation , Physica D, 220: 2006, 1–13.[12] G. Fibich, F. Merle, P. Rapha¨el, Erratum to “Proof of a spectral property related to the singularity formationfor the L critical nonlinear Schr¨odinger equation” , personal communication.[13] G. Fibich and G. Papanicolaou, Self-focusing in the perturbed and unperturbed nonlinear Schr¨odinger equationin critical dimension , SIAM J. Appl. Math., 60 (1999), 183-240.[14] G. M. Fraiman,
Asymptotic stability of manifold of self-similar solutions in self-focusing , Soviet Phys. JETP61 (1985), no. 2, 228-233; translated from Zh. ´Eksper. Teoret. Fiz. 88 (1985), no. 2, 390–400 (Russian).[15] J. Ginibre, G. Velo,
On a class of nonlinear Schr¨odinger equations I. The Cauchy problem, general case , J.Funct. Anal. 32(1)(1979), 1–32.[16] R. Glassey,
On the blowing up of solutions to the Cauchy problem for nonlinear Schr¨odinger equation , J.Math. Phys. 18 (1977), no. 9, 1794–1797.[17] M. Goldman, K. Rypdal, B. Hafizi,
Dimensionality and dissipation in Langmuir collapse , Phys. Fluids 23,1980, 945–955.[18] C. Guevara,
Global behavior of finite energy solutions to the focusing Nonlinear Schr¨odinger Equation in d -dimensions. Appl Math Res Express (2014) vol.2017, 177–243.[19] J. Holmer, R. Platte and S. Roudenko,
Blow-up criteria for the 3d cubic NLS equation , Nonlinearity, 23(2010), 977–1030.[20] J. Holmer and S. Roudenko,
Blow-up solutions on a sphere for the 3d quintic NLS in the energy space , Analysis& PDE, 5-3 (2012), 475–512.[21] J. Holmer and S. Roudenko,
A class of solutions to the 3d cubic nonlinear Schr¨odinger equation that blow-upon a circle , AMRX Appl Math Res Express (2011), vol. 2011, 23–94. doi:10.1093/amrx/abq016[22] J. Holmer and S. Roudenko,
On blow-up solutions to the 3d cubic nonlinear Schr¨odinger equation , AMRXAppl. Math. Res. Express, 1 (2007), article ID abm004, 31 pp., doi:10.1093/amrx/abm004[23] N. E. Kosmatov, V. F. Shvets, V. E. Zakharov,
Computer simulation of wave collapses in the nonlinearSchr¨odinger equation , Physica D 52 (1991), 16–35.[24] M.K. Kwong,
Uniqueness of positive solutions of ∆ u − u + u p = 0 in R n , Arch. Ration. Mech. Anal. 105 (3)(1989), 243–266. TABLE BLOWUP IN NLS IN HIGH DIMENSIONS 39 [25] M. J. Landman, B. LeMesurier, G. Papanicolaou, C. Sulem, P.-L. Sulem,
Singular solutions of the cubicSchr¨odinger equation, “integrable systems and applications” , Lecture Notes in Physics, 342 (2005), Springer-Verlag 207–217.[26] M. J. Landman, G. Papanicolaou, C. Sulem, P.-L. Sulem,
Rate of blowup for solutions of the nonlinearSchr¨odinger equation at critical dimension , Phys. Rev. A(3) 38(8) (1988), 3837–3843.[27] M. J. Landman, G. Papanicolaou, C. Sulem, P.-L. Sulem, X. Wang,
Stability of isotropic singularities for thenonlinear Schr¨odinger equation , Phys. Rev. A(3) 38(8) (1988), 3837–3843.[28] B. LeMesurier, G. Papanicolaou, C. Sulem, P.-L. Sulem,
Local structure of the self-focusing singularity of thenonlinear Schr¨odinger equation , Physica D 32 (1988), 210–226.[29] P. M. Lushnikov, S. A. Dyachenko and N. Vladimirova,
Beyond leading-order logarithmic scaling in thecatastrophic self-focusing of a laser beam in Kerr media , Physical Review A, v. 88 (2013), 13–845.[30] V. M. Malkin,
On the analytical theory for stationary self-focusing of radiation , Physica D, 64 (1993), 251–266.[31] J. Marzuola, G. Simpson,
Spectral Analysis for Matrix Hamiltonian Operators , Nonlinearity, vol. 24 (2011),389–429.[32] D. McLaughlin, G. Papanicolaou, C. Sulem, P.-L. Sulem,
Focusing singularity of the cubic Schr¨odinger equa-tion , Physical Review A, vol. 34 (1986) 1200–1210.[33] F. Merle,
On uniqueness and continuation properties after blow-up time of self-similar solutions of nonlinearSchr¨odinger equation with critical exponent and critical mass , Comm. Pure Appl. Math. 45 (1992), no. 2,203-254.[34] F. Merle,
Determination of blow-up solutions with minimal mass for nonlinear Schr¨odinger equations withcritical power , Duke Math. J. 69 (1993), no. 2, 427-454.[35] F. Merle, P. Rapha¨el,
Blow up dynamic and upper bound on the blow up rate for critical nonlinear Schr¨odingerequation , Ann. Math. 161 (2005), no. 1, 157–222.[36] F. Merle, P. Rapha¨el,
Sharp upper bound on the blow up rate for critical nonlinear Schr¨odinger equation ,Geom. Funct. Anal. 13 (2003), 591–642.[37] F. Merle, P. Rapha¨el,
On universality of blow up profile for L critical nonlinear Schr¨odinger equation , Invent.Math. 156 (2004), 565–672.[38] F. Merle, P. Rapha¨el, Sharp lower bound on the blow up rate for critical nonlinear Schr¨odinger equation , J.Amer. Math. Soc. 19 (2006), no. 1, 37–90.[39] Galina Perelman,
On the blow up phenomenon for the critical nonlinear Schr¨odinger equation in 1D , Ann.Henri Poincar´e 2(2001), 605–673.[40] Galina Perelman,
On the blow up phenomenon for the critical nonlinear Schr¨odinger equation in 1D ,S´eminaire: ´Equations aux D´eriv´ees Partielles, 19992000, Exp. No. III, 16 pp., ´Ecole Polytech., Palaiseau,2000.[41] P. Rapha¨el,
Stability of the log-log bound for blow up solutions to the critical nonlinear Schr¨odinger equation ,Math. Ann. 331 (2005), no. 3, 577–609.[42] M. Reed, B. Simon,
Method of modern mathematical physics IV. Analysis of Operators.
Academic Press 1978.[43] K. Rypdal, J. Rasmussen,
Blow-up in nonlinear Schr¨odinger equations II , Phys. Scripta 33, 1986, 498–504.[44] G. Simpson, I. Zwiers,
Vortex collapse for the L -critical nonlinear Schr¨odinger equation , J. Math. Phys.,52(8), (2011), 83–503.[45] J. Shen, T. Tang, L. Wang, Spectral Method, algorithms, analysis and applications.
Springer 2015.[46] C. Sulem, P.-L. Sulem,
The nonlinear Schr¨odinger equation. Self-focusing and wave-collapse.
Springer 1999.[47] P.-L. Sulem, C. Sulem, A. Patera,
Numerical simulation of singular solutions of the two-dimensional cubicSchr¨odinger equation , Com. Pure and Appl. Math, v. 37 (1984), 755–778.[48] V. Synakh, V. Zakharov,
The nature of the self-focusing singularity , Sov. Phys. JETP 41 (1975), 465–468.[49] E. C. Titchmarsh,
Eigenfunction expansions associated with second-order differential equations , Oxford,Clarendon Press, 1946.[50] L. N. Trefethen,
Spectral Method in Matlab.
SIAM 2000.[51] S. N. Vlasov, V. A. Petrishchev, and V. I. Talanov,
Averaged description of wave beams in linear and nonlinearmedia (the method of moments) , Radiophysics and Quantum Electronics 14 (1971) pp. 1062–1070. Translatedfrom Izvestiya Vysshikh Uchebnykh Zavedenii, Radiofizika, 14 (1971) pp. 1353–1363. [52] S. N. Vlasov, L. V. Piskunova, V. I. Talanov,
Structure of the field near a signularity arising from self-focusingin a cubically nonlinear medium , Sov. Phys. JETP 48. 1978, 808–812.[53] M. I. Weinstein,
Nonlinear Schr¨odinger equations and sharp interpolation estimates , Comm. Math. Phys. 87(1983) 567–576.[54] M. I. Weinstein,
Modulational stability of ground states of nonlinear Schr¨odinger equations , SIAM J. Anal.16 (1985) 472–491.[55] M. I. Weinstein,
Lyapunov stability of ground states of nonlinear dispersive evolution equations , Comm. PureAppl. Math., Volume 39 (1986) 51–68.[56] D. Wood,
The self-focusing singularity in the nonlinear Schr¨odinger equation , Stud. Appl. Math. 71 (1984),103–115.[57] Kai Yang,
Formation of singularities in nonlinear dispersive PDE , Pro-Quest LLC, 2018, Thesis (PhD) – TheGeorge Washington University; arXiv:1712.07647[58] V. E. Zakharov,
Collapse of Langmuir waves , Zh. Eksp. Teor. Fiz. 62, (1972), 1745–1751, (in Russian); Sov.Phys. JETP, 35 (1972), 908–914 (English)., Zh. Eksp. Teor. Fiz. 62, (1972), 1745–1751, (in Russian); Sov.Phys. JETP, 35 (1972), 908–914 (English).