Collapse vs. blow up and global existence in the generalized Constantin-Lax-Majda equation
NNoname manuscript No. (will be inserted by the editor)
Collapse vs. blow up and global existence in thegeneralized Constantin-Lax-Majda equation
Pavel M. Lushnikov , · Denis A.Silantyev · Michael Siegel October 6, 2020
Abstract
The question of finite time singularity formation vs. global exis-tence for solutions to the generalized Constantin-Lax-Majda equation is stud-ied, with particular emphasis on the influence of a parameter a which controlsthe strength of advection. For solutions on the infinite domain we find a newcritical value a c = 0 . . . . below which there is finite timesingularity formation that has a form of self-similar collapse, with the spatialextent of blow-up shrinking to zero. We find a new exact analytical collapsingsolution at a = 1 / a in the analytical continuation of the solu-tion from the real spatial coordinate into the complex plane. This singularitycontrols the leading order behaviour of the collapsing solution. For a c < a ≤
1, we find a blow-up solution in which the spatial extent of the blow-up regionexpands infinitely fast at the singularity time. For a (cid:38) .
3, we find that thesolution exists globally with exponential-like growth of the solution amplitudein time. We also consider the case of periodic boundary conditions. We iden-tify collapsing solutions for a < a c which are similar to the real line case. For a c < a ≤ .
95, we find new blow-up solutions which are neither expanding norcollapsing. For a ≥ , we identify a global existence of solutions. Department of Mathematics and Statistics, University of New Mexico, Albuquerque,MSC01 1115, NM, 87131, USA Landau Institute for Theoretical Physics, Chernogolovka, Moscow region, 142432, Russia Courant Institute of Mathematical Sciences, New York University, 251 Mercer Street NewYork, NY 10012-1110, USA Department of Mathematical Sciences and Center for Applied Mathematics and Statistics,New Jersey Institute of Technology, Newark, NJ 07102, USA a r X i v : . [ n li n . PS ] O c t P.M. Lushnikov, D.A. Silantyev and M. Siegel
In this paper we investigate finite-time singularity formation in the generalizedConstantin-Lax-Majda (CLM) equation [8,11,34] ω t = − auω x + ωu x , ω, x ∈ R , t > ,u x = H ω, (1)which is a 1D model for the advection and stretching of vorticity in a 3Dincompressible Euler fluid. Here ω and u are a scalar vorticity and velocity,respectively, a ∈ R is a parameter, and H is the Hilbert transform, H ω ( x ) := 1 π p.v. (cid:90) + ∞−∞ ω ( x (cid:48) ) x − x (cid:48) d x (cid:48) . (2)This equation, with a = 0, was first introduced by Constantin, Lax and Ma-jda [8] as a simplified model to study the possible formation of finite-timesingularities in the 3D incompressible Euler equations. It was later general-ized by DeGregorio [11] to include an advection term uω x , and by Okamoto,Sakajo and Wensch [34], who introduced the real parameter a to give differentrelative weights to advection and vortex stretching, u x ω . In addition to its re-lationship to the 3D Euler equation, (1) has a direct connection to the surfacequasi-geostrophic (SQG) equation [17].The 3D incompressible Euler equations can be written as ∂ t ω + u · ∇ ω = ω · ∇ u , x ∈ R , t > , (3) u = ∇ × ( − ∆ ) − ω . (4)The second equation above is the Biot-Savart law, which in free-space has anequivalent representation as a convolution integral u ( x , t ) = 14 π (cid:90) R ( x − y ) × ω ( y , t ) | x − y | d y . (5)The term ω · ∇ u on the right-hand side (r.h.s.) of (3), where ∇ u = S ( ω ) isa matrix of singular integrals, is known as the vortex stretching term. Stan-dard estimates from the theory of singular integral operators [37] show that (cid:107) ω (cid:107) L p ≤ (cid:107)∇ u (cid:107) L p ≤ c p (cid:107) ω (cid:107) L p for 1 < p < ∞ , which formally implies that thevortex stretching term scales quadratically in the vorticity, i.e., S ( ω ) ω ≈ ω .This term is therefore destabilizing and has the potential to generate singularbehavior. However, analysis of the regularity of Eqs. (3), (4) is greatly compli-cated by the nonlocal and matrix structure of S , and remains an outstandingopen question.In contrast to the vortex stretching term, the advection term u · ∇ ω does not cause any growth of vorticity. As a result, it has historically beenthought to play an unimportant role in the regularity of the incompress-ible Euler and Navier-Stokes equations. Recent studies, however, show that ollapse vs. blow up 3 advection-type terms can have an unexpected smoothing effect. For exam-ple, Hou and Lei [27] present numerical evidence that a finite-time singularityforms from smooth data in solutions to a reformulated version of the Navier-Stokes equations for axisymmetric flow with swirl, when the so-called convec-tion terms u r ∂ r ( ω θ /r ) + u z ∂ z ( ω θ /r ) and u r ∂ r ( u θ /r ) + u z ∂ z ( u θ /r ) are omitted.Here ( u r , u θ , u z ) and ω θ are velocity and vorticity components in cylindri-cal coordinates ( r, θ, z ). Adding the convection back is found to suppress afinite-time singularity formation. Related work on the smoothing effect of ad-vection/convection in the Euler and Navier-Stokes equations is given in [20,21,22,23,24,33].The generalized CLM equation (1) (also called the Okamoto-Sakajo-Wunschmodels in Ref. [17]) is obtained from the 3D Euler equations by replacing theadvection term u · ∇ ω with uω x and the vortex stretching term S ( ω ) ω by its1D analogue H ( ω ) ω . The Hilbert transform H is the unique singular integraloperator in 1D that preserves certain important properties of S ( ω ) [8]. In ad-dition, the 1D vortex stretching term H ( ω ) ω preserves the quadratic scaling ofthe vortex stretching term S ( ω ) ω in the 3D problem. The resulting equation(1) provides a simplified setting to understand the competition between thestabilizing effect of advection and destabilizing effect of vortex stretching. Inthis work we focus on smooth (analytic or C ∞ ) initial data which we consideras the most physically relevant. There are also a number of results on singu-larity formation for (1) in the case of Holder continuous initial data, see Refs.[6,17] for recent reviews.We summarize some of the known results, concentrating on those whichapply to smooth (analytic or C ∞ ) initial data. In the case a = 0, Constantin,Lax and Majda [8] obtained a closed-form exact solution to the initial valueproblem for (1) which develops a self-similar finite-time singularity for a classof analytic initial data. When a (cid:54) = 0, the simplifications that enable a closed-form solution no longer hold, and various analytical and numerical methodshave been applied to investigate singularity formation. Castro and Cordoba [5]proved finite-time blow-up for a < a > (cid:15) − small values of a >
0, vortexstretching dominates and Elgindi and Jeong [17] proved the existence of self-similar finite-time singularities in the form ω = 1 τ f ( ξ ) , ξ = xτ α , τ = t c − t, (6)where t c > α depends on a , approaching α = 1 inthe limit a → . Also, f ( ξ ) is an odd function, i.e. f ( − ξ ) = − f ( ξ ) , ξ ∈ R . Theproof of [17] is based on a continuation argument in a small neighborhood ofthe exact solution at a = 0. Chen, Hou and Huang [6] proved a similar resultusing a different method.The special case of a = 1 of Eq. (1) was first considered by De Gregorio[11] and has been the subject of extensive numerical computations in the P.M. Lushnikov, D.A. Silantyev and M. Siegel periodic geometry by Okamoto, Sakajo and Wensch [34]. These suggest thatsingularities do not occur in finite-time from smooth initial data on a periodicdomain. Okamoto et al. [34] use a least squares fit to the decay of Fouriermodes to track the distance δ ( t ) from the real line to the nearest singularity inthe complex- x plane. They find that δ ( t ) decays exponentially in time, whichis consistent with global existence. Global existence for a = 1 in the specificcase of non-negative (or non-positive) initial vorticity is proven by Lei et al.[28].The above analytical and numerical results might suggest the existence ofa threshold value a = a threshold below which finite-time singularities occur forsmooth initial data, and at/above which the solution exists globally in time.Okamoto et al. conjecture that a threshold = 1. However, for this value a = 1,Chen et al. [6] recently proved the existence of an “expanding” self-similarsolution (6) for the problem on x ∈ R . In this solution f ( ξ ) is an odd functionwith finite support and α = −
1. It implies that ω ( x, t ) → f (cid:48) (0) x as t → t c for any finite value of x ∈ R while the boundary of compact support expandsinfinitely fast in the spatial coordinate x as t → t c . The form of this solutionis apparently incompatible with the periodic geometry, and thus does not ruleout the possibility of global existence of the solution in that geometry when a = 1. In addition, Ref. [6] assumes smooth initial data with finite support ω ( x ) ∈ C ∞ c so it remains an open question if analytic initial data convergesto the expanding self-similar solution.We are not aware of any theory or simulation which consider solutions to(1) over a wide range of the parameter a as well as any simulation on x ∈ R addressing even the particular case a = 1. The main goal of this paper is tofill this gap by presenting theory and highly accurate computations to assesssingularity formation for a wide range of a for both the periodic geometry and x ∈ R .We obtain three main analytical results (Theorems 1-3 below). The firstone establishes the specific form of the leading-order complex singularity of f ( ξ ) in (6) and determines its dependence on a , when that singularity is ofpower-law type. The second one finds an exact solution for f ( ξ ) at a = 1 / α = 1 / . The third one proves that this type of exactsolution is impossible beyond the particular cases a = 0 and 1 / . Our spectrally accurate numerical simulations address all real values of a .We use a variable numerical precision, beyond the standard double precision,to mitigate loss of accuracy when computing poles and branch points in thecomplex plane, and employ fully resolved spatial Fourier spectra on an adap-tive grid with 8th order adaptive time stepping. Computations are performedboth for periodic boundary conditions (BC) as well as on the real line x ∈ R with the decaying BC ω ( x, t ) → x → ±∞ . (7)For the problem on R , we reformulate Eq. (1) in a new spatial variable q using a conformal mapping from Ref. [30] between the real line x ∈ R and q ∈ ( − π, π ) . Then our spectral simulations with a uniform spatial grid for ollapse vs. blow up 5 q ∈ ( − π, π ) ensure spectral precision on the corresponding highly non-uniformgrid for x ∈ R . Our results make use of two distinct types of numerical simulation. The firsttype is time-dependent simulation which allows us to establish the convergenceof generic initial conditions to the self-similar solution (6). As a by-product ofsuch simulations, we obtain values of α and the functional form of f ( ξ ) . Thesecond type of simulation directly solves the nonlinear eigenvalue problem for α to obtain the similarity solution (6) of Eq. (1) for each value of a. We solve thatnonlinear eigenvalue problem by iteration on the real line x ∈ R using a versionof the generalized Petviashvili method (GPM) [36,29,26,35,13]. In Theorem4 we show that there exists a nonstable eigenvalue for the linearization of theoriginal Petviashvili method [36] which prevents its convergence. However, theversion of GPM employed here avoids that instability.The results of the first and the second type of simulation are in excellentagreement with Theorems 1-3. The first major result of these simulations isthe discovery of a critical value a = a c = 0 . . . . (8)below which (i.e., for a < a c ) there is finite-time singularity formation, but atwhich point (i.e., for a = a c ) the singularity transitions or changes character.For a < a c the value of α is positive with f ( ξ ) an analytic function in a strip inthe complex plane of ξ containing the real line. The positive values of α ensure,in accordance with Eq. (6), that the solution shrinks in x as t → t c while thesolution amplitude diverges in that limit. This type of shrinking self-similarsolution is compatible with both kinds of boundary conditions (i.e., periodicand decaying on R ), and our simulations reveal the same type of singularityformation at t → t c . The shrinking and divergence of amplitude is qualitativelyreminiscent of the collapse in both the nonlinear Schr¨odinger equation and thePatlak-Keller-Segel equation, see e.g. Refs. [41,7,38,3,25,31]. The terminology“collapse” or “wave collapse” was first introduced in Ref. [41] in analogy withgravitational collapse and has been widely used ever since. The singularityformation found for a < a c is therefore of collapse type. We also find that α = 0 at the critical value a = a c .The second major result of our simulations is the uncovering of a quali-tatively different type of self-similar singularity formation for a c < a ≤
1, inwhich the spatial scale of the solution does not shrink. We refer to this type ofsingularity as “blow up.” An additional finding in the aforementioned rangeof the parameter a is that the blow-up solution on the real line x ∈ R andthe blow-up solution for periodic BC are qualitatively different. In the case x ∈ R we find that − ≤ α < α = − a = 1 . Thus Eq. (6)corresponds to an expanding self-similar solution. In particular, at a = 1, wefind that α = − x = 0 results in ω ( x, t ) = τ − − α xf (cid:48) (0) + O ( τ − − α x ) . It shows that the linear slope ∝ x increases to infinity as t → t c for a c ≤ a < a = 1 . Time-dependent simulations for x ∈ R P.M. Lushnikov, D.A. Silantyev and M. Siegel with analytic initial conditions and a c ≤ a ≤ t → t c to Eq. (6) with f ( ξ ) being of finite support. It extendsthe results of Ref. [6] from a = 1 to a c ≤ a ≤
1, as well as from compactly sup-ported initial conditions ω ( x ) ∈ C ∞ c to (noncompact) analytic data. For thisanalytic data, we find the time-dependent solution converges to an expandingself-similar solution of finite support as t → t c .The third major result of our simulations concerns periodic BC. While thecollapse case a < a c is similar for both x ∈ R and periodic BC, as mentionedthe case a c < a ≤ a c ≤ a ≤ x ∈ R would contradict the periodic BCas t approaches t c . Instead, we find a new self-similar blow-up solution ω ( x, t ) = 1 t c − t f ( x ) , (9)which is valid for a c < a ≤ . . Formally, we can interpret Eq. (9) as Eq. (6)with α = 0. However, periodic BC are qualitatively different from the finitesupport solution of Eq. (6) because of the nonlocality of the Hilbert transformin Eq. (1). We find that f ( x ) in Eq. (9) has a discontinuity at the periodicboundary in a high-order (or n th-order) derivative, such that n → ∞ in thelimit a → a + c , i.e. f ( x ) approaches a C ∞ function in that limit. A complexsingularity is also present in f ( x ) on the imaginary axis away from the realline, the form of which obeys Theorem 1.In the range 0 . < a < , our simulations are inconclusive regardingwhether blow up occurs. The value a = 1 is a special case for the periodicBC, with no blow up observed in our simulations for generic initial conditions.Instead, the solution exists globally with the first spatial derivative remainingbounded, while the second derivative grows exponentially in time. This agreeswith the result on global existence for the particular case a = 1 investigatedin Ref. [34].For a ≥ , we find that the solution exists globally for all initial conditionsconsidered in the case of periodic BC, while for the solution on the real linethe situation is not conclusive. In the latter case, the maximum of | ω | initiallygrows with time but this growth saturates at larger times at least for a (cid:38) . < a (cid:46) . , our simulations catastrophically looseprecision at sufficiently large times, and a conclusive determination betweenblow up and global existence of solutions is not possible.We also find from the simulations that the kinetic energy on the infiniteline x ∈ R , E K := ∞ (cid:90) −∞ u ( x, t )d x, (10)with an initially finite value approaches a constant as t → t c when a < . ± . . ± . < a ≤
1. In the case a (cid:38) . t → ollapse vs. blow up 7 ∞ . On the periodic domain x ∈ [ − π, π ] , we find the same behaviour of thekinetic energy up to a = 0 .
95. For a ≥ E K approaches a non-zero constant as t → ∞ ( a = 1) or tends to zero ( a > ω ( x, t ) and f ( ξ ) in the complexplane of x and ξ, we use both a fitting of the Fourier spectrum similar to Ref.[34] (see also Refs. [4,14,15,39] for more detail), and more general methodsof analytical continuation by rational interpolants (see Refs. [1,15,12,32]). Astime evolves, these singularities approach the real line in agreement with Eq.(6). We have formulated a system of ordinary differential equations (ODEs) de-scribing the motion of such singularities. Fourier fitting allows us to track onlysingularities which are nearest to the real axis, while rational interpolants gobeyond this, by giving information on singularities other than the closest one.In particular, it reveals that for a (cid:54) = 0 , / a < a c , there are genericallybranch points beyond the leading order singularities, consistent with Theorem3. The exceptional cases are a = 0 , / / a = 2 /
3, the third order pole coexists with additional branch points. Forother values of a , the nearest singularities are branch points. We find that for a c < a ≤ , the singularities approach the real line as t → t c in the spatialregions near the boundary of the support of f ( ξ ) . The rest of this paper is organized as follows. Section 2 establishes Theorem1, which describes the leading-order complex singularity and determines itsdependence on a . Section 3 reinterprets the results of Ref. [8] for a = 0 interms of moving complex poles and the self-similar solution (6). In Section 4,we derive an exact blow-up solution for a = 1 / a and establishes in Theorem 3 that, except for a = 0 , / x ∈ R are developed in Sections 6and 7. In particular, Section 6 reformulates Eq.(1) as a nonlinear eigenvalueproblem for the self-similar solution (6), and Section 7 rewrites Eq. (1) in anauxiliary variable q mapping the real line into the finite interval. Section 8 thendescribes the results of time-dependent numerical simulations for x ∈ R , andSection 9 presents self-similar solutions of the type (6) via numerical solutionof the nonlinear eigenvalue problem using a generalized Petviashvili method.Section 10 addresses the analytical continuation into the complex plane of x by rational approximation and uses it to study the structure of singularities.Section 11 describes the results of both time-dependent numerical simulationsand the generalized Petviashvili method for periodic BC. Section 12 providesa summary of the results and discusses future directions. Appendix A gives aderivation for the form of the Hilbert transform over x in variable q. P.M. Lushnikov, D.A. Silantyev and M. Siegel
We assume that ω ( x, t ) is an analytic function in the open strip containing x ∈ R in the complex plane x ∈ C decaying at x → ±∞ . Then we canrepresent ω as ω = ω + + ω − , (11)where ω + ( x.t ) is analytic in the upper complex half-plane x ∈ C + and ω − ( x.t )is analytic in the lower complex half-plane x ∈ C − .The Hilbert transform (2) implies that H ω = − i( ω + − ω − ) . (12)Assume that the solution exhibits a leading order singularity of power γ > x for ω at x = ± i v c , v c > , so that ω ( x, t ) = ω − γ ( t )[ x − i v c ( t )] γ + ¯ ω − γ ( t )[ x + i v c ( t )] γ + l.s.t., (13)where l.s.t designates less singular terms at x = ± i v c , i.e.lim x →± i v c [ x ∓ i v c ( t )] γ l.s.t. = 0 . (14)If we additionally assume that ω ( − x ) = − ω ( x ) , for x ∈ R , then Eq. (13)implies that ω − γ ( t )[ x − i v c ( t )] γ + ¯ ω − γ ( t )[ x + i v c ( t )] γ = − ω − γ ( t )[ − x − i v c ( t )] γ − ¯ ω − γ ( t )[ − x + i v c ( t )] γ , (15)i.e., ¯ ω − γ ( t )( − γ +1 = ω − γ ( t ) . Then we can define ω − γ ( t ) := − i e − i πγ/ ˜ ω − γ ( t ) , ˜ ω − γ ( t ) ∈ R (16)so that Eq. (13) takes the following form ω ( x, t ) = − i˜ ω − γ ( t ) (cid:18) e − i πγ/ [ x − i v c ( t )] γ − e i πγ/ [ x + i v c ( t )] γ (cid:19) + l.s.t.. (17)Using Eqs. (1),(12) and (17) we obtain that u x = H ω = ˜ ω − γ ( t ) (cid:18) e − i πγ/ [ x − i v c ( t )] γ + e i πγ/ [ x + i v c ( t )] γ (cid:19) + l.s.t., (18)and u := u + + u − = − ˜ ω − γ ( t )( γ − (cid:18) e − i πγ/ [ x − i v c ( t )] γ − + e i πγ/ [ x + i v c ( t )] γ − (cid:19) + l.s.t., (19)where we have additionally assumed that γ (cid:54) = 1. ollapse vs. blow up 9 Plugging Eqs. (17)-(19) into Eq. (1) and collecting the most singular terms ∝ [ x − i v c ( t )] − γ at x = i v c ( t ) on the right-hand side of Eq. (1) givesi e − i πγ ˜ ω − γ ( t )[ x − i v c ( t )] γ (cid:18) aγγ − − (cid:19) = 0 . (20)By assumption ω − γ ( t ) (cid:54) = 0. Then Eq. (20) implies that γ = 11 − a . (21)Thus we have proved the following: Theorem
1. If a solution of Eq. (1) has a complex conjugate pair of powerlaw singularities located at x = ± i v c for v c > γ > , then γ is determined by Eq. (21). Remark 1 . The condition γ > γ <
0, then the leading order term in Eq. (1) at x = ± i v c is ∝ [ x − i v c ( t )] . Remark 2 . Eq. (21) is in excellent agreement with the simulations of Section8. The singularities with γ < a = 0 results in γ = 1 . Also γ → ∞ for a → − . For theparticular values a = n − n , n = 1 , , , . . . , (22)we obtain the integer values γ = n resulting in complex pole singularities oforder n in Eq. (17) while the other values of a ∈ (0 ,
1) result in the branchpoints at x = ± i v c ( t ) . a = 0 The particular value of the parameter a = 0 implies from Eq. (21) that γ = 1.This case recovers the results of Ref. [8]. The general solution of Eq. (1) isimmediately obtained by noticing that Eqs. (1), (12) result in ω t = ω + t + ω − t = − i( ω + ) + i( ω − ) (23)which decouples into two independent ODEs ω + t = − i( ω + ) , ω − t = i( ω − ) . (24)The solutions of these ODEs with the generic initial conditions ω + ( x, t ) | t =0 = ω +0 ( x ) and ω − ( x, t ) | t =0 = ω − ( x ) are given by ω + ( x, t ) = ω +0 ( x )1 + i tω +0 ( x ) and ω − ( x, t ) = ω − ( x )1 − i tω − ( x ) . (25) Eqs. (11), (12) and (25) lead to the solution of Constantin-Lax-Majda equationfound in Ref. [8] ω ( x, t ) = 4 ω ( x )[2 − t H ω ( x )] + t ω ( x ) (26)for the generic initial condition ω ( x, t ) | t =0 = ω ( x ) = ω +0 ( x ) + ω − ( x ). AlsoEqs. (12) and (25) imply that (as in Ref. [8]) H ω ( x, t ) = 2 H ω ( x )[2 − t H ω ( x )] − tω ( x )[2 − t H ω ( x )] + t ω ( x ) . (27)Assume that there exists an x ∈ R such that ω ( x ) = 0 and H ω ( x ) >
0. Then Eq. (26) implies a singularity in the solution at the time t c :=2 / H ω ( x ) > . If there are multiple points x ∈ R such that ω ( x ) = 0and H ω ( x ) > t c := 2 / sup {H ω ( x ) | ω ( x ) = 0 } > x corresponds to the singularity at the earliest time t = t c . Aparticular example is any odd function ω ( x ) with respect to x = x (implyingthat ω ( x ) = 0) which is strictly positive for x > x and decays at x → ∞ .A series expansion of Eq. (26) at x → x and t → t − c implies that ω ( x, t ) = 1 t c − t ξω (cid:48) ( x )[ H ω ( x )] ([ H ω ( x )] − ξ H ω (cid:48) ( x )) + 4 ξ [ ω (cid:48) ( x )] + O (( t c − t ) ) , (28)where ξ := x − x t c − t (29)is the self-similar variable. Eqs. (28) and (29) provide a universal profile ofthe solution at t → t − c in a spatial neighborhood of x → x after we neglectthe correction term O (( t c − t ) ). That profile has the form of a sum of twocomplex poles at complex conjugate points ξ = ξ ± as follows ω ( x, t ) = i t c − t (cid:18) ξ + ξ − ξ + − ξ − ξ − ξ − (cid:19) , (30)where ξ ± = [ H ω ( x )] H ω (cid:48) ( x ) ± i ω (cid:48) ( x )] (31)are positions of poles in the complex plane of ξ. Eqs. (30) and (31) provide the exact solution of Eq. (1) for ω (cid:48) ( x ) < ω (cid:48) ( x ) < ξ + ∈ C + . This solution is asymptoticallystable with respect to perturbations of the initial condition as follows from Eq.(28). The only trivial change due to the perturbation of the initial conditionis a shift of both x and t c . ollapse vs. blow up 11 One can also recover from the solution (30) the representation (17) with γ = 1 which gives the exact solution ω ( x, t ) = − ˜ v c (cid:18) x − x − i˜ v c ( t c − t ) + 1 x − x + i˜ v c ( t c − t ) (cid:19) = − ˜ v c t c − t (cid:18) ξ − i˜ v c + 1 ξ + i˜ v c (cid:19) (32)of Eq. (1) for any values of the real constants t c , ˜ v c > x . Here withoutloss of generality we have shifted the origin in the real direction compared withthe solution (30). a = 1 / The particular value of the parameter a = 1 / γ = 2. In this section we look for the solution to Eq. (1) in the form (17)assuming that the l.s.t. are identically zero, i.e., ω ( x, t ) = i˜ ω − ( t ) (cid:18) x − x − i v c ( t )] − x − x + i v c ( t )] (cid:19) , (33)where for generality we have also allowed a shift of the origin by introducingthe arbitrary real constant x . Eq. (19) then becomes u = ˜ ω − ( t ) (cid:18) x − x − i v c ( t ) + 1 x − x + i v c ( t ) (cid:19) = 2˜ ω − ( t )( x − x )( x − x ) + v c ( t ) . (34)Plugging Eqs. (33) and (34) into Eq. (1), we find the latter equation is iden-tically satisfied provided dv c ( t ) dt = − ˜ ω − ( t )4 v c ( t ) , (35)and d ˜ ω − ( t ) dt = ˜ ω − ( t )4 v c ( t ) . (36)Solving the system of ordinary differential equations (ODEs) (35) and (36)results in v c ( t ) = ( t c − t ) / ˜ v c , ˜ ω − ( t ) = 4˜ v c t c − t ) / , (37)where ˜ v c > t c are two arbitrary real constants. Assuming the initialcondition is given at t = 0 and that t c >
0, we obtain that t = t c is the timeof singularity formation.Section 8 below shows the convergence during the evolution in time t ofthe solution of Eq. (1) to the exact solution given by Eqs. (33) and (37). The spatial extent of the solution shrinks while the maximum amplitude increasesuntil the singularity is reached at t = t c .One can rewrite the solution (33), (37) in the self-similar form as follows ω ( x, t ) = 1 t c − t v c (cid:18) ξ − i˜ v c ] − ξ + i˜ v c ] (cid:19) = 1 t c − t v c ξ ξ + ˜ v c ) , (38)where ξ := x − x ( t c − t ) / (39)is the self-similar variable.The exact solution presented in this section, the best of our knowledge, isnew. To summarize, in this section we have proven the following theorem: Theorem
2. Eqs. (38) and (39) provide an exact solution of Eq. (1) for a = 1 / t c , ˜ v c > x . Remark 3 . The decay of u ( x, t ) in Eq. (34) as x → ±∞ ensures that thekinetic energy (10) has a finite value for t < t c . In contrast, E K for the solution(32) at a = 0 is infinite. a The explicit self-similar solutions (29)-(31) and (38), (39) (corresponding to thevalues a = 0 , /
2) represent the particular situation where the leading ordersingularity in Eqs. (17) and (21) provides the exact solution with identicallyzero l.s.t. . All other values of a are addressed in the following theorem: Theorem
3. A solution (17) and (21) of Eq. (1) which decays at x → ±∞ requires l.s.t. which are not identically zero for any a ∈ R except a = 0 and a = 1 / . Proof
The case a ≥ a = 1 corresponds to the singular valueof γ as follows from Eq. (21), while a > γ <
0, contradicting theassumption of Theorem 2 that ω at x → ±∞ . Thus below we assume that a < γ > l.s.t. in Eq. (17) are identically zero.Then we plug Eq. (17) into Eq. (1) and collect terms with different powers of x − i v c ( t ) . The most singular term ∝ [ x − i v c ( t )] − γ is identically zero by Eq.(21) as follows from the proof of Theorem 1. Collecting the next most singularterms ∝ [ x − i v c ( t )] − − γ we obtain that dv c ( t ) dt = − − γ ˜ ω − γ ( t ) v γ − c ( t ) γ , (40)which generalizes Eq. (35) to arbitrary values of γ. We note that there is nooverlap between terms of different orders in this proof except in the case γ = 1,for which − γ = − γ −
1. However, this case is fully considered in Section 3 andexcluded by assumption in the statement of Theorem 3 because it correspondsto a = 0 . ollapse vs. blow up 13 Collecting the terms ∝ [ x − i v c ( t )] − γ we obtain that d ˜ ω − γ ( t ) dt = 2 − γ ( γ − ω − γ ( t ) v γc ( t ) (41)which generalizes Eq. (36) to arbitrary values of γ. However, at the next order, collecting terms ∝ [ x − i v c ( t )] − γ +1 leads to2 − γ − ( γ − γ + 1)i e − i πγ/ ˜ ω − γ ( t ) v γ +1 c ( t ) = 0 , (42)which cannot be satisfied by any nontrivial solution ˜ ω − γ ( t ) (cid:54)≡ γ = 2,i.e. a = 1 / . This contradiction completes the proof of Theorem 3.
Remark 4 . The ODE system (40) and (41) can be immediately solved forany γ resulting in v c ( t ) = ˜ v c ( t c − t ) γ ( γ +1) , ˜ ω − γ ( t ) = 2 γ ˜ v γc γ + 1 ( t c − t ) − γγ +1 , (43)where ˜ v c and t c are arbitrary real constants. Then neglecting l.s.t. , we obtainfrom Eqs. (17) and (43) the following self-similar “solution” ω ( x, t ) = − i t c − t γ ˜ v γc γ + 1 (cid:18) e − i πγ/ [ ξ − i˜ v c ] γ − e i πγ/ [ ξ + i˜ v c ] γ (cid:19) , (44)where ξ := x − x ( t c − t ) α , α = 2 γ ( γ + 1) , (45)is the self-similar variable. For γ = 1( a = 0) and γ = 2( a = 1 / γ (cid:54) = 1 ,
2. Onemay hope that even if γ (cid:54) = 1 ,
2, the self-similar solution is well approximated byEqs. (44) and (45) because (17) is the leading order singularity of the solution.However, we find below in Section 8 (see also Fig. 1) that the numericallycomputed self-similar solution has a different power scaling for ξ = x − x ( t c − t ) α than in Eq. (45), i.e. α (cid:54) = α for γ (cid:54) = 1 ,
2. This implies that the l.s.t , neglectedin (45), lead to a non-trivial modification of α compared with α . Discussion.
To address l.s.t. in Eq. (17), one can naively write a solutionas a formal expansion in powers of x ± i v c ( t ) as follows ω ( x, t ) = ω + ( x, t ) + ω − ( x, t )= ∞ (cid:88) j = − γ − i˜ ω j ( t ) (cid:16) e i πj/ [ x − i v c ( t )] j − e − i πj/ [ x + i v c ( t )] j (cid:17) , (46)where we have assumed integer values of γ as given by Eq. (22). For example,if a = 2 / γ = 3, such a series has the form Fig. 1
Dependence of α ( a ) on a , obtained via time-dependent simulations of Section 8 andvia nonlinear eigenvalue problem of Section 9. The green curve terminates at a = a c sincethe iteration used to solve the nonlinear eigenvalue problem for x ∈ R does not converge for a > a c . Also included for comparison is an approximation to α ( a ) from Eq. (45), α ( a ) = γ ( a )( γ ( a )+1) = − a ) (2 − a ) . ω ( x, t ) = ω + ( x, t ) + ω − ( x, t )= ∞ (cid:88) j = − − i˜ ω j ( t ) (cid:16) e i πj/ [ x − i v c ( t )] j − e − i πj/ [ x + i v c ( t )] j (cid:17) . (47)For non-integer values of γ we can write a more general Puiseux series involvingmultiple powers of nγ + m, n, m ∈ N , however, this is beyond the scope of thecurrent paper.A limitation of the formal series (46) and (47) is that they simultaneouslyinvolve expansions of ω + ( x, t ) at w = − i v c and ω − ( x, t ) at w = i v c . If theradiuses of convergence r c of both series exceed v c , then there is overlap of thecorresponding disks of convergence, e.g., at x = 0 . The simplest situation iswhen both radiuses of convergence exceed 2 v c , which would turn these formalseries into proper Laurent series at x = ± i v c for integer values of γ , whileallowing a straightforward evaluation of the Hilbert transform by Eq. (12).However, we find from the simulations of Section 10 that additional singu-larities exist in the complex plane of x beyond x = ± v c , which implies that v c < r c < v c for γ (cid:54) = 1 ,
2. Thus for integer values of γ (cid:54) = 1 , , one has to use amore general form for the Laurent series at x = ± i v c which is not necessarilyequivalent to (46). For example, in the case γ = 3 one has to replace Eq. (47)by the Laurent series ω ( x, t ) = ∞ (cid:88) j = − − i˜˜ ω j ( t ) e i πj/ [ x − i v c ( t )] j , (48)at x = i v c , where the series coefficients are the same for the singular part inEq. (47), i.e. ˜˜ ω j ( t ) ≡ ˜ ω j ( t ) for j = − , − , −
3, while for j ≥ ollapse vs. blow up 15 are different than ˜ ω j ( t ). One can also prove that ˜ ω − ( t ) ≡ ˜ ω − ( t ) ≡
0. Thedisadvantage of using the Laurent series (48) compared with (47) is that onecannot easily find H ω from (48) because both ω + ( x, t ) and ω − ( x, t ) contributeto the coefficients ˜˜ ω j ( t ). A technique from Ref. [12] can be used to to obtainpartial information on H ω , however a full solution (i.e. a knowledge of allcoefficients) requires global information. The results of Sections 3-5 suggest looking for a solution of Eq. (1) in thegeneral self-similar form (6). Substitution of the ansatz (6) into Eq. (1) reducesit to M f := f + αξf ξ = − a ( ∂ − ξ H f ) f ξ + f H f, (49)where M is a linear operator. One can also rewrite Eq. (49) as the system f + αξf ξ = − agf ξ + f g ξ , g = ∂ − ξ H f, (50)where u = τ α − g ( ξ ) . (51)We note that the series (46) or (48) are reduced to self-similar form pro-vided we add the restriction that v c ( t ) = ( t c − t ) α ˜ v c and ˜ ω − δ ( t ) = ( t c − t ) − αδ ˜˜ ω − δ , where ˜˜ ω − δ is an arbitrary constant and the subscript/power δ represents a general subscript/power in the series starting from δ = γ .We can iterate Eq. (49) for different values of α to find the optimal α which realizes the dominant collapse regime. To do this we have to invert theoperator M in Eq. (49) at each iteration. The equation M f = 0 has a generalsolution f ∝ | ξ | − α (52)for α (cid:54) = 0 and f ≡ α = 0 . Depending on the sign on α , this solutionis singular either at x → x → ±∞ . Thus the operator M is invertiblefor the class of smooth solutions decaying at x → ±∞ which we use below inSection 9.The condition that the solution of Eq. (49) decays at both x → ±∞ requiresa specific choice of α for each a . It forms a version of nonlinear eigenvalueproblem for α ( a ) . Section 9 finds α ( a ) by iterating Eq. (49) numerically. Asymptotics for ξ → ±∞ . Since both terms in r.h.s. of Eq. (50) arequadratic in f , then the decay f ( ξ ) , g ( ξ ) → ξ → ±∞ implies that Eq.(52) is valid for ξ → ±∞ provided α >
0. It agrees with the exact resultsof Section 3 (Eq. (30)) and Section 4 (Eq. (38)) for α = 1 and α = 1 / α <
0, Eq. (52) is inconsistent with the condition f ( ξ ) → . Thus we conclude that f ( ξ ) ≡ ξ → ±∞ for α < . (53)It implies that f ( ξ ) has the finite support for α < α = − Eq. (49) is invariant under a stretching of the self-similar coordinate ξ , ξ → Aξ, A = const ∈ R . (54)i.e., if f ( ξ ) is a solution for Eq. (49) then f ( Aξ ) is also a solution of the sameequation. Therefore if one finds a solution of Eq. (49) then it immediatelyimplies an infinite family of solutions from the stretching (54). Despite thisnonuniqueness, we find that the version of GPM employed here converges toa solution of Eqs. (50), (53). Further details are given in Section 9. The analysis of previous sections assumes the solution exists on the real line x ∈ ( −∞ , ∞ ) with the decaying BC (7). To address this infinite domain insimulations, we use the auxiliary (computational) variable q defined by x = tan (cid:16) q (cid:17) . (55)Eq. (55) maps the segment of the real line ( − π, π ) of q onto the real line( −∞ , ∞ ) of x. Extending both x and q into the complex plane, we find thatEq. (55) maps the infinite strip − π < Re ( q ) < π onto the complex plane x ∈ C , except for the half-lines ( − i ∞ , − i) and (i , +i ∞ ), with the upper half-strip being mapped onto the upper half-plane C + and the lower half-stripbeing mapped onto the lower half-plane C − . Also the boundaries of the strip, Re ( q ) = ± π are mapped onto ( − i ∞ , − i) and (i , +i ∞ ), see e.g., Refs. [15,30]for details of this mapping. Here and below we abuse notation and use thesame symbols for functions of either x or q . For example, we assume that˜ f ( q ) := f ( x ( q )) and remove the ˜ sign.Using the Jacobian of the mapping (55), dxdq = 12 cos ( q ) = 11 + cos q , (56)and the results of Appendix A, we rewrite Eqs. (1)-(2) for independent vari-ables q and t as ω t = − a (1 + cos q ) uω q + ω [ H π ω + C πω ] , q ∈ ( − π, π ) , (1 + cos q ) u q = [ H π ω + C πω ] , (57)where the Hilbert transform H π on the interval ( − π, π ) is defined by (see alsoAppendix A) H π f ( q ) := 12 π p.v. (cid:90) π − π f ( q (cid:48) )tan( q − q (cid:48) ) d q (cid:48) , (58)and the constant C πω is determined by C πω = − π (cid:90) π − π ω ( q (cid:48) ) tan( q (cid:48) q (cid:48) . (59) ollapse vs. blow up 17 We call Eq. (57) the transformed CLM equation. Note that Eq. (58) is thereduction of Eq. (2) to the class of 2 π -periodic functions, see Appendix A. Thedecaying BC (7) allows a 2 π -periodic extension of ω ( q, t ) with ω ( q, t ) | q = π +2 πn =0 , n ∈ N . It enables us to work with ω ( q, t ) in terms of a Fourier series over q. Based on the results of Section 7, we numerically solve Eq. (57) on the real line x ∈ R with a pseudo-spectral Fourier method by representing the 2 π -periodicsolution ω ( q, t ) as a sum of 2 N Fourier modes ˆ ω k ( t ) as ω ( q, t ) = k = N − (cid:88) k = − N ˆ ω k ( t ) e i kq . (60)We use 2 N uniformly spaced grid points in q from − π to π − ∆q , where ∆q = π/N . The Fast Fourier transform (FFT) allows us to efficiently findnumerical values of ˆ ω k ( t ) from values of ω ( q, t ) on that grid. The resolution N is chosen depending on the initial condition (IC) and adaptively adjustedthroughout the computation so that the spectrum ˆ ω k is fully resolved with thedesired precision. This means that | ˆ ω k | decays by 16-17 orders of magnitude at | k | ∼ N compared to max − N ≤ k ≤ N − | ˆ ω k | , down to the round-off floor of the errorfor double precision. For the multi-precision simulations which were performed,this decay is further enhanced (or equivalently, the round-off is reduced) byany desired number of orders. Below we focus on the description of doubleprecision simulations while noting that higher precision simulations were alsoextensively performed.The decay of the Fourier spectrum ˆ ω k is checked at the end of every timestep. If | ˆ ω k | is larger than the numerical round-off at | k | ∼ N at the giventime-step, then the simulation is ‘rewound’ for one time-step backwards with N increased by factor of 2, and the time-stepping is continued. Amplitudes ofthe new extra Fourier modes are set to 0, which is equivalent to performinga spectral interpolation of the solution at the newly inserted grid points in q space. Rewinding is done to avoid accumulation of error due to the tailsof the spectrum not being fully resolved at the time-step before the grid re-finement. For time-marching we use 11-stage explicit Runge-Kutta methodof 8 th order [9] with the adaptive time-step ∆t determined by the condition ∆t = CFL · min { ∆q/ ( a max q | (1 + cos q ) u | ) , / max q | (1 + cos q ) u q |} , where thenumerical constant CFL is typically chosen as CFL = 1 / , / /
16 toachieve numerical stability in the time-stepping and ensure that the error of themethod is near round-off level. Also, the scaling of ∆t with max q | (1 + cos q ) u | and max q | (1 + cos q ) u q | ensures numerical stability of the method during possi-ble singularity formation events. We additionally enforced the real-valuednessof ω ( q ) at each time-step to avoid numerical instability, since the FFT and inverse FFT lead to accumulation of a small imaginary part at the level ofround-off, which can be amplified during time evolution.Typically, we used the following two types of initial conditions (ICs):IC1: ω ( q ) = − (sin( q ) + 0 . q )) , (61)IC2: ω ( q ) = i 4 V c T c (cid:18) q ) − i V c ) − q ) + i V c ) (cid:19) , (62)where the real-line IC1 is similar in form to the periodic IC in Ref. [34] exceptfor an opposite sign. In IC2, V c and T c are the real numbers and in most ofour simulations we used V c = 1, T c = 1, for which IC2 reduces to ω ( q ) = −
43 (sin( q ) + 0 . q )) . (63)Note the first two derivatives of (63) are zero at q = ± π , i.e., ω ( n )0 ( q = ± π ) =0 for n = 0 , ,
2. Both ICs (61) and (62) are real-valued odd functions with anegative slope at q = 0, and lead to the formation of a singularity at q = 0at some moment in time for a < a c (see Eq. (8) for the definition of a c )while ω ( q, t ) stays real-valued and odd. The function ω ( q ) in IC1 is an entirefunction, and that in IC2 has two double poles at x = tan (cid:0) q (cid:1) = ± i V c in x -space or at q = ± i q c in q -space, where q c = 2 arctanh( V c ) . Note that IC2corresponds to the exact solution for the case a = 1 / t = T c (see Eq. (38)), while for other values of the parameter a, it is not an exactsolution but qualitatively resembles one on the real interval [ − π, π ] and servesas a good IC to obtain collapsing solutions.Computation of the 2 π -periodic Hilbert transform H π (see Appendix Afor the definition of H π ) is easily done in Fourier space asˆ H πk = − i sign( k ) , (64)where sign( k ) = 1 for k >
0, sign( k ) = 0 for k = 0 and sign( k ) = − k < C πω (59) in Eq. (57) is computed from the condition that H π ω ( q = − π ) + C πω = 0, i.e. − i k = N − (cid:80) k = − N ˆ ω k ( − − k sign( k ) + C πω = 0.While computing the values of u q from the second equation in (57), onehas to take special care at the point q = − π . Expanding both the left-handside (l.h.s.) and r.h.s. of that equation in a Taylor series at the point q = − π, we obtain that u q ( q = − π ) = H πqq ω ( q = − π ), which can also be computedusing ˆ ω k . The term with H πq in the Taylor series of the r.h.s. vanishes since H πq ω ( q = − π ) = (cid:80) k | k | ˆ ω k = 0 for the real-valued odd function ω ( q ), sincefor such a function we have ˆ ω − k = − ˆ ω k .For each simulation we made a least squares fit of the Fourier spectrum | ˆ ω k | at time t to the asymptotic decay model | ˆ ω k ( t ) | ≈ C ( t ) e − δ ( t ) | k | | k | p ( t ) (65) ollapse vs. blow up 19 for | k | (cid:29) C ( t ) , δ ( t ) and p ( t ) are the fitting parameters for eachvalue of t . This allows us to obtain both δ ( t ) > p ( t ) as functions of t .The value of δ ( t ) indicates the distance of the closest singularity of ω ( q ) fromthe real line in the complex q -plane, and the value of p ( t ) is related to thetype or power of that complex singularity, see Refs. [34,14,15,39] for moredetails. In particular, if the singularity in the solution is of a power law type ω ( q ) ∼ ( q − i q c ) − γ then using complex contour integration one obtains (seee.g. Ref. [4]) that | ˆ ω k | ≈ Ce − q c | k | / | k | − γ , meaning that δ = q c and p = 1 − γ (66)which follows from Eq. (65). According to Eq. (55), the distance δ x from theclosest singularity to the real line in the complex x -plane is δ x = tanh (cid:0) δ (cid:1) . Itimplies that δ x = δ + O ( δ ) for δ (cid:28) a = 2 / V c = 1, T c = 1 (i.e. Eq. (63)) are provided in Figs. 2 and 3. The maximal valuemax q | ω ( q, t ) | of the numerical solution increases from an initial value ∼ ∼ at the final simulation time t final . Fig. 3 shows the spectrum | ˆ ω k | and its fit to the model (65). This fit provides numerically extracted values ofboth δ ( t ) and p ( t ). Then δ x ( t ) = tanh (cid:16) δ ( t )2 (cid:17) is computed from δ ( t ) and fittedto δ x ( t ) ∝ ( t c − t ) α per Eq. (6). From the fit we obtain that α ≈ . x -space. The algebraic decay rate p ( t ) appears to stabilize at the value − t approaches the singularity time t c . An initial transient is not included in thedata used for the δ x ( t ) fit, since δ ( t ) and p ( t ) cannot be determined accuratelyat these times due to the spectrum | ˆ ω k | being oscillatory. These oscillationsquickly die out as the self-similar regime is approached.In order to extract α from the fit δ x ( t ) ∝ ( t c − t ) α with maximum accu-racy, we need an accurate estimate of t c . We obtain this estimate from the fitmax x | ω ( x, t ) | ∝ t c − t ) by extrapolating the numerical solution up to t = t c .Using this t c for the fit instead of t final (that one may naively use instead ofthe correct t c ) provides considerably more accurate values of α , for example,giving 7 significant digits instead of 4 when a = 2 /
3. Having a more accurateestimate of t c is especially important for simulations with a < . δ and p while fitting | ˆ ω k | tothe model (65) if we confine the least square fit to a window of data between1/4 and 1/3 of the total effective width of the spectrum (shown on the leftpart of Fig. 3 with a green color). This is due to an increase in the relativeerror of the spectrum data at the tails, as the round-off floor is approached.Moreover, the model (65) is accurate only asymptotically as | k | → ∞ so wecannot use too small values of | k | either.For 0 ≤ a < a c (with a c given by Eq. (8)) and for both IC1 (61) andIC2 (62), we find that δ x ( t ) evolves in time toward 0 while p ( t ) approaches Fig. 2
Results of the simulation of Eqs. (57)-(58) with a = 2 / ω ( q, t ), its derivative ω q ( q, t ) and u ( q, t ) for t = 1 . α extracted from the simulations as explained in the text. Thecollapse time t c is extracted from the fit (by extrapolation) to max | ω ( x, t ) | ∝ t c − t ) . Fig. 3
Left panel: The Fourier spectrum | ˆ ω k | at a particular time t = 1 . a = 2 /
3. The red line represents a fit to the model (65)with green line showing the portion of the | ˆ ω k | used for the least-squares fit. Center andright panels: time dependence of δ x ( t ) = tanh (cid:16) δ ( t )2 (cid:17) and p ( t ) recovered from the fit of thespectrum to Eq. (65) at different times. The red solid line at the center panel represents afit to the model δ x ( t ) ∼ ( t c − t ) α . ollapse vs. blow up 21 Fig. 4
Convergence of time-dependent numerical solution of Eqs. (57)-(58) with a =1 / a = 2 / t = t c , t c ≈ . t c ≈ . x -space, where x = tan( q ). Horizontal and vertical scalesare dynamically changed in both panels to exactly match the positions and amplitudes ofthe local maximum at x = x max and minimum at x = − x max . Fig. 5
Convergence of the time-dependent numerical solution of Eqs. (57)-(58) to the self-similar profile (6) as t → t c . Here a = 2 / ω ( x ) = − (cid:18) x − x +1 ) + x − x − ) (cid:19) − i3 (cid:18) x − x +2 ) − x − x − ) (cid:19) + (cid:18) x − x +3 ) + x − x − ) (cid:19) + i96 (cid:18) x − x +4 ) − x − x − ) (cid:19) , where x ± = − ± i4 , x ± = − ± i2 , x ± = 1 ± i4 , x ± = ± i8 . The solution is shown at two different moments in time, where for each time we overlaidthe self-similar profile as in Fig. 4, matching their corresponding maximum and minimumpositions horizontally and vertically. a constant value after a quick transient phase, see Fig. 3 (right panel). Weobserve spontaneous formation of a universal self-similar solution profile ofthe form (6) during time evolution (see Fig. 4). These self-similar profiles, aswell as the value of α in δ x ( t ) and the terminal value of p ( t ) as t → t c are thesame for a wide class of ICs (e.g. one can change a power of singularity in IC2from − − V c > T c > a . Table 1 provides the universal values of α and p vs a . Fig. 1 Fig. 6
Dependence of γ num ( a ) = 1 − p ( a ) using p ( a ) obtained via time-dependent sim-ulations by the fit to Eq. (65). These data are also provided in Table 1. Also shown is γ ( a ) = − a from Eq. (21) for comparison. Here we plot 1 /γ ( a ) instead of γ ( a ) for the easiercomparison. shows the dependence of α ( a ) on a . However, one can also find particular IC inwhich finite time singularities do not form. Two such choices are -IC1 and -IC2,i.e. IC1 (61) and IC2 (62) taken with the opposite sign. In these two cases wedid not observe collapse or singularity formation in finite time, but rather analgebraic-in-time approach of a singularity to the real line, δ x ( t ) ∼ /t µ , µ > ω ( q ) were found to either decay algebraically intime or approach constant values while higher-order derivatives of ω ( q ) growalgebraically in time. Other smooth generic initial conditions that were triedwere found to produce blow up after an initial transient, as exemplified in Fig.5. These transients made the simulation considerably slower (due to the needfor more modes in the spectrum of to resolve the solution down to doubleprecision round-off). However, in a space-time neighborhood of the singularitythese solutions recover the same self-similar profile as shown in Fig. 4, see alsoFig. 5. We note that the velocity u ( x, t ) evolves toward the self-similar profile(51) with max x | u | → ∞ for 0 < a < a c . Below we focus on IC1 and IC2, butthe reader should but keep in mind that they appear generic.Using the terminal values of p extracted by fits to Eq. (65) with various a , and employing Eq. (66) to recover γ from p , we confirmed the formula γ ( a ) = − a (see Theorem 1 and Eq. (21) in Section 2) and the correspondingformula p ( a ) = − a − a within 0.5% for 0 ≤ a < a c . Fig. 6 shows the numericalapproximation, γ num ( a ) = 1 − p ( a ) using values of p ( a ) from Table 1 as wellas the theoretical value γ = − a for comparison. We note that the plot of1 /γ num ( a ) in Fig. 6 stops at a = a c , since it is difficult to obtain accuratevalues of p ( a ) (and hence γ num ( a )) from time-dependent simulations when a > a c . This is due to a transition that occurs at a = a c , in which the fittedsingularity for a < a c corresponding to collapse is no longer closest to thereal- x line when a > a c . ollapse vs. blow up 23 In addition to Fourier fitting, we also extract values of α in an alternativeway (these values are called α below), using the spatial derivative of theself-similar solution (6) given by ω x ( x, t ) = 1( t c − t ) α f (cid:48) (cid:18) x ( t c − t ) α (cid:19) . (67)Using Eq. (67) we fit max x | ω x ( x, t ) | to the model max x | ω x ( x, t ) | ∝ t c − t ) α tofind α . The same scaling is valid for max q | ω q ( q, t ) | since the maximum of | ω q | (or | ω x | ) is at x = q = 0 for an odd collapsing solution. Fig. 2 (center rightpanels) shows such a fit for a = 2 /
3, where α ≈ . . For com-parison, α extracted via δ x ( t ) is α ≈ . α and α have errors of only 1 × − and 1 . × − , respec-tively (we compared them to α e = 0 . α for various a are also gathered in Table 1 forcomparison with values of α . We confirmed that α and α obtained using theabove two methods for 0 < a (cid:46) .
689 agree within a relative error of < . a < x | ω | → ∞ as t → t c according to the self-similar profile in Eq.(6). The extracted values of α , p and α for a < a = − u ( x, t ) during the temporal evolution approaches the self-similar profile (51)near the singularity location at x = q = 0. A qualitative difference for a < < a < a c ) is that the self-similar profile (51) approacheszero because α > u ( x, t ) is generally nonzero, even at t → t c . This extendsthe result of [5], who proved that there is finite-time singularity formation for a < ω ( x, ∈ C ∞ c ( R ) with H ω (0 , >
0, to examples with analytic initial data.The values of α ( a ) were challenging to obtain with more than 3-4 digitsof accuracy using time-dependent simulations in double precision arithmetic,since it involves the fitting of several unknown coefficients, including t c . Gen-erally, fitting to max x | ω x ( t ) | ∝ / ( t c − t ) α provided more accurate resultsfor α then obtaining α from fitting to δ x ( t ) ∝ ( t c − t ) α (for a < a c ). This isespecially true for values of a near a c since δ x ( t ) ∼ ( t c − t ) α decays very slowlywith α → a → a c , while amplitudes satisfy max x | ω ( t ) | ∝ / ( t c − t ) andmax x | ω x ( t ) | ∝ / ( t c − t ) α . At a = 2 / α = 0 . . . . ) the quantity δ x ( t )decays only by one order of magnitude over the simulation time (see Fig. 3,center panel) even though max x | ω ( t ) | grows from from the value ∼ (see Fig. 2, right panel, the first row).We obtained much more accurate values of α ( a ) (up to 14 digits of preci-sion) by numerically solving the nonlinear eigenvalue problem, Eq. (50), for aself-similar solution of Eq. (1) (see Section 9). In contrast, for a c we were able Fig. 7
Results from simulations of Eqs. (57)-(58) with a = − ω ( q, t ), its derivative ω q ( q, t ), and u ( q, t ) for t = 0 . α extracted from simulations as explained in the text. Fig. 8
Left panel: The Fourier spectrum | ˆ ω k | at time t = 0 . | ˆ ω k | used for the fit. Center and right panels: Time dependence of δ x ( t ) = tanh (cid:16) δ ( t )2 (cid:17) and p ( t ) recovered from fit of the spectrum to Eq. (65) at different times. The red solid lineat the center panel represents a fit to the model δ x ( t ) ∝ ( t c − t ) α . ollapse vs. blow up 25 Fig. 9
Left panel: Convergence of the time-dependent numerical solution to Eqs. (57)-(58)with a = 1 and IC2 (63) to a self-similar profile with compact support. The solution expandshorizontally and stretches vertically until blowing up at t = t c ≈ . x -space, where x = tan( q ), and is scaled both horizontally and vertically toexactly match the positions of the local maximum and minimum. Center and right panels:The time dependencies of max x | ω ( x, t ) | and the absolute value its location x max ( t ) on t . to obtain 14 digits of accuracy using both time-dependent simulations and thenonlinear eigenvalue problem with double precision arithmetic. Another 3 dig-its of precision are obtained (for a total of 17 digits of precision) if quadrupleprecision arithmetic is used in the nonlinear eigenvalue problem.We have also performed simulations specifically with a = 1 since this spe-cial case was addressed in Chen et al. [6] who proved for a = 1 the existenceof an “expanding” self-similar solution of the type (6) for the problem on x ∈ R . In this case f ( ξ ) is an odd function with a finite support and α = − ω ( x, t ) → f (cid:48) (0) x as t → t c for any finite value of x ∈ R while the boundary of the compact support expands infinitely fast intolarge | x | as t → t c . Ref. [6] assumes smooth initial data with finite support ω ( x ) ∈ C ∞ c , so it remains an open question if analytic initial data wouldconverge to this expanding self-similar solution. Our numerical findings in-deed show an approach to this kind of expanding solution starting from ananalytic IC, see Figs. 9 and 10. The solution grows in amplitude and expandsfaster than exponentially in time, which is demonstrated by semi-log plots ofmax x | ω ( x ) | ( t ) and its location x max ( t ) in the middle and right panels of Fig. 9.The solution obeys the self-similar profile (6) and forms a finite time singular-ity at t = t c . Fig. 10 (right panels) confirms the scales max x | ω ( x ) | ∝ / ( t c − t )and | ω x ( x = 0) | ∝ / ( t c − t ) α = const with α = −
1. One can also see(from the middle panel of Fig. 10) that max x | ω xx ( x ) | → ∞ as t → t c . We areable to simulate the growth in amplitude of ω ( x ) only by about one order ofmagnitude with our spectral code, since the spectrum widens very quickly as t → t c and decays slowly, i.e., | ˆ ω k ( x ) | ∼ k − , as seen in Fig.11 (left panel).The approach to a self-similar solution with compact support is expressed inthe complex x -plane by the approach of complex singularities (identified asbranch points from our simulations) located at x = x sing to the real line nearthe boundaries of compact support. The small distances | Im ( x sing ) | of thesesingularities to the real line for t near t c means that the solution is “almost of Fig. 10
Results of the same simulation as in Fig. 9 with a = 1 showing solution ω ( q, t ), ω q ( q, t ) and u ( q, t ) in q -space (left panels) as well as the same solution in x -space (cen-ter panels) at time t = 1 . . Right panels show the time dependence of their maxi-mum values as functions of ( t c − t ), where t c is the blow-up time extracted from the fit tomax | ω ( x, t ) | ∝ t c − t ) . Fig. 11
Left panel: The Fourier spectrum | ˆ ω k | at time t = 1 . a = 1. The red line represents a fit to the model (65) with green line showinga portion of the | ˆ ω k | used for the fit. Right panel: Time dependence of δ x ( t ) = tanh (cid:16) δ ( t )2 (cid:17) recovered from the fit of the spectrum to Eq. (65). The red solid line in the right panelrepresents a fit to the model δ x ( t ) ∝ ( t c − t ) α . ollapse vs. blow up 27 compact support” with “almost a jump” in the first derivative at the boundaryof ”compact support” in x -space. The singularity locations scale like x sing (cid:39) ± ( t c − t ) α x b ± i ( t c − t ) α y b (68)(i.e. there are four symmetrically located singularities), where α = − α ≈ . . Here the real constants t c , x b and y b depend on the IC. Note that α is different from α because it characterizes the approach of the solutionto the compactly supported profile (6). In contrast, the value α = − α suggests that the “almostcompactly supported” solution turns into a truly compactly supported solu-tion at t = t c , with a jump in the first derivative. Due to oscillations in thespectrum, it is difficult to accurately extract the value of α from the fit to δ x ( t ) ∼ ( t c − t ) α . However, using rational approximation via the AAA al-gorithm (see details about AAA in Section 10) we can observe two pairs ofbranch cuts with branch points approach the real line near x = ± ( t c − t ) α x b as t → t c , similar to the case a = 0 .
8. One can see from Fig. 12 (right panel)that the structure of the singularity for a = 0 . a = 1 case.For a c < a < a = 1 case, but involves different values of α . Another differencecompared to the a = 1 case is that there is a discontinuity in a higher-orderderivative at the boundary of ”compact support,” instead of a jump in the firstderivative ω x as occurs for a = 1. Figs. 12 - 14 show the results of simulationswith the parameter a = 0 . ω xx formingat the boundary of “compact support.” In this case, we are able to simulate thegrowth of the amplitude of ω ( x ) only by a factor of ∼ . This limitation isdue to the rapid widening of the spectrum with time so that it becomes almostalgebraic, | ˆ ω k ( x ) | ∼ k − p , p ≈ t → t c . This is again due to the solutionachieving “almost compact support” with a jump in the second derivative atthe boundary of the “compact support” in x -space. Fig. 13 (right) shows thegrowth of both max x | ω ( x ) | and max x | ω x ( x ) | = | ω x ( x = 0) | as functions of t c − t confirming the scales max x | ω ( x ) | ∼ / ( t c − t ) and | ω x ( x = 0) | ∼ / ( t c − t ) α with α = − . a = 1, for a c < a < t → t c according to Eq. (68).For example, when a = 0 . α = − . α ≈ . δ ( t )and p ( t ) from a fit to Eq. (65) due to the spectrum being oscillatory, see theleft panel of Fig. 14. The right panel of Fig. 14 provides the best fit which wewere able to obtain for δ ( t ). The fitting parameter p ( t ) was more sensitive tothe oscillations and did not appear to stabilize at any particular value, so wedo not provide a plot for it here. Fig. 12
Left panel: Convergence of time dependent numerical solution to Eqs. (57)-(58)with a = 0 . t = t c ≈ . x -space, where x = tan( q ), and is scaled both horizontally and verticallyto exactly match the positions of the local maximum and minimum. Center panel: The timedependence x max ( t ) of the location of max x ω ( x ). The dashed lines shows that it scales like x max ( t ) ∝ ( t c − t ) α with α (cid:39) − . t → t c . Right panel: The structure of complexsingularities at t = 1 . ω ( x ) ≈ ω AAA ( x ) = (cid:80) m − i =1 a i x − b i . Thesimple poles are shown as dots at locations b i with a size of dot scaled with log | a i | .The branch cuts are shown as lines connecting the dots, and form ‘U-shaped’ curvesin the upper and lower complex plane. The accumulation of poles approximates twopairs of branch points near the real line. The location of these branch points scale as x sing ∼ ± ( t c − t ) α x b ± i y ( t c − t ) α y b , where x , y > α = − . α ≈ . This type of oscillation in the spectrum occurs when there are two sym-metric singularities that are equally close to the real line. In this case, a moreelaborate fitting procedure with additional parameters to account for the os-cillation can yield improved results, see e.g. Ref. [2]. However, such fits are alsomore delicate to implement, and are beyond the scope of the current work.Simulations with ICs either of type -IC1 or -IC2 and a c < a ≤ x | ω ( x, t ) | and max x | u ( x, t ) | . The maximum slopemax x | ω x ( t ) | = | ω x ( x = 0 , t ) | is found to approach a constant value for a = 1while it decays for a < . Also, max x | ω xx ( x, t ) | grows algebraically as a functionof t , while δ x ( t ) decays algebraically, δ x ( t ) ∼ /t µ , µ >
0. Since these ICs donot result in a finite-time singularity formation, we do not discuss these casesin further detail.For a (cid:38) . ω has the form an an expanding self-similar functionwhich approaches a compactly supported profile (in the scaled variable ξ ) withinfinite slope at the boundary of the compact region, so that max x | ω | → x | ω x | , max x | u | → ∞ as t → ∞ (although ω x ( x = 0) → t → ∞ ).The complex singularities approach the real line in infinite time with positionsthat scale like x sing = ± x exp ( κ t ν ) ± i y exp ( − κ t ν ), where the constants κ , κ , ν , ν > a . In particular, for a = 2 we observe an initialgrowth followed by the eventual decay of max x | ω | . For a = 1 . , the amplitudeof ω initially grows and then plateaus over the times we were able to reach in ollapse vs. blow up 29 Fig. 13
Results of the same simulation as in Fig. 12 with a = 0 . ω ( q, t ), ω q ( q, t ),and u ( q, t ) in q -space (left panels) as well as in x -space (center panels) at time t = 1 . . Right panels show the time dependence of their maximum values as functions of ( t c − t ),where t c is the blow-up time extracted from the fit to max | ω ( x, t ) | ∝ t c − t ) . Fig. 14
Left panel: The Fourier spectrum | ˆ ω k | at a particular time t = 1 . a = 0 .
8. The red line represents a fit to the model (65)with the green line showing the portion of the | ˆ ω k | used for the fit. The purple line shows a fitto the rougher model (65) with δ = 0. Right panel: Time dependence of δ x ( t ) = tanh (cid:16) δ ( t )2 (cid:17) recovered from the fit of the spectrum to Eq. (65). The red solid line at the right panelrepresents a fit to the model δ x ( t ) ∝ ( t c − t ) α . simulations. For a = 1 . x | ω | so thatit is below an exponential rate. For both -IC1 and -IC2 we again observe globalexistence of the solution with decay of ω and infinite growth of ω x ( x = 0), withan infinite slope forming at x = 0 and a singularity approaching the real linelike x sing = 0 ± i y exp ( − κ t ν ), where y , κ , ν > . For 1 < a (cid:46) . , we find from simulations that initially max x | ω | grows. Thisperiod of initial growth is long, with the spectrum widening so quickly thatit was challenging to distinguish between a finite time singularity and globalexistence when a is near 1, but we have numerical evidence of global existencefor a at least as small as 1.3, as described in the previous paragraph.Here we summarize the behaviour of solutions to Eqs. (57)-(58) on x ∈ R ,and its dependence on the parameter a , for quite generic smooth IC: – a < a c with α ( a ) >
0: Collapse in ω, i.e. max x | ω | → ∞ at the finite time t c . As t → t c , solutions with generic IC approach the shrinking universal self-similar profile (6) near the spatial location of max x | ω | . As t → t c , the profilesshrink to zero width. The self-similar solution has leading order complexsingularities in agreement with Theorem 1 and Eq. (21). The location ofthese singularities approaches the real line as x sing = x ± i δ x ( t ), where δ x ( t ) ∝ ( t c − t ) α , α = α ( a ) >
0. In particular, x = 0 for both IC1 or IC2.Also u ( x, t ) near x follows the self-similar profile (51) with max x | u | → ∞ for 0 < a < a c . – a c < a ≤ α ( a ) <
0: Blow up in both ω and u at the finite time t c . As t → t c , solutions with generic IC approach the expanding self-similarprofile Eq. (6) which has compact support. As t → t c , the rate of expansionturns infinite. The complex singularities closest to the real line correspondto the boundaries of compact support, and they approach the real line as x sing ∼ ± ( t c − t ) α x b ± i ( t c − t ) α y b , where α = α ( a ) < α ( a ) > – a (cid:38) . x | ω | →
0, max x | ω x | , max x | u | →∞ and ω x ( x = 0) → t → ∞ . The complex singularities approach thereal line exponentially in time as x sing = ± x exp ( κ t ν ) ± i y exp ( − κ t ν ),where κ , κ , ν , ν > Similar to the transformation of Eq. (1) to Eqs. (57)-(58) in Section 7, weobtain a transformed equation for self-similar solutions of Eq. (50) by mappingthe interval ( − π, π ) of the auxiliary variable q onto the real line ( −∞ , ∞ ) as ξ = tan (cid:16) q (cid:17) . (69) ollapse vs. blow up 31 With this mapping Eq. (50) turns into M f := f + α sin q f q = − a (1 + cos q ) gf q + f [ H π f + C πf ] := N f, q ∈ [ − π, π ] , (1 + cos q ) g q = H π f + C πf , (70)where the 2 π -periodic Hilbert transform H π and the constant C πf are definedin Eqs. (58), (59), and the linear operator M is now defined in q space by thel.h.s. of the first Eq. in (70). We also define the quadratically nonlinear operator N such that N f represents the r.h.s. of the first Eq. in (70) with g expressedthrough the second equation in (70) as g = ∂ − q (cid:34) H π f + C πf (1 + cos q ) (cid:35) , ∂ − q p := q (cid:90) − π p ( q (cid:48) ) dq (cid:48) . (71)Then Eq. (70) takes the following operator form M f = N f. (72)A linearization of Eq. (72) about f together with Eqs. (70) and (71) resultin L δf := −M δf − a (1 + cos q ) ∂ − q (cid:34) H π δf + C πδf (1 + cos q ) (cid:35) f q − a (1 + cos q ) ∂ − q (cid:34) H π f + C πf (1 + cos q ) (cid:35) δf q + δf [ H π f + C πf ] + f [ H π δf + C πδf ] , (73)where L is the linearization operator and δf is the deviation from f. Taking δf = f in Eq. (73) and using Eqs. (70), (72) to remove the nonlinearterms in f , we prove the following theorem: Theorem
4. The solution f of Eq. (70) satisfies the relation L f = M f. (74) Corollary 1.
The invertability of the operator M (see Section 6) and Eq.(74) imply that the operator M − L has the eigenvalue λ = 1 with eigenfunc-tion f , which is the same as the solution f of Eq. (70).Similar to Eq. (60), we approximate a solution of Eq. (70) as a truncatedFourier series f ( q ) = k = N − (cid:88) k = − N ˆ f k e i kq . (75) Then the discrete Fourier transform allows us to rewrite Eq. (70) in matrixform as M ˆ f = (cid:100) N f , M := − αk αk − αk αk . . .. . . . . . − αk N αk N − , (76)where ˆ f = ( ˆ f k , ˆ f k , . . . , ˆ f k N ) T is a column vector, the tridiagonal matrix M ∈ R N × N represents the Fourier transform of the operator M and (cid:100) N f is the column vector of Fourier coefficients of N f . Also k := − N, k := − N +1 , . . . , k N := N − . Note that the tridiagonal form of M is a consequenceof the term sin( q ) = e i q − e − i q in the definition of M in Eq. (70).We solve Eq. (74) in the truncated Fourier representation (76) by iterationusing the generalized Petviashvili method (GPM) [26] which relates the n +1thiteration ˆ f n +1 to the n th iteration ˆ f n of ˆ f as followsˆ f n +1 − ˆ f n = (cid:32) [ − ˆ f n + M − (cid:91) N f n ] − (cid:18) ∆τ (cid:19) (cid:104) ˆ f n , − M ˆ f n + (cid:91) N f n (cid:105)(cid:104) ˆ f n , M ˆ f n (cid:105) ˆ f n (cid:33) ∆τ, (77)where superscripts give the iteration number, (cid:104) a , b (cid:105) := (cid:80) k = N − k = − N ¯ a k b k is thecomplex dot product and ∆τ is a parameter that controls the convergence rateof the iterations. At each iteration we need to solve Eq. (76) for ˆ f (assuming (cid:100) N f is given) to effectively compute M − (cid:91) N f n . Since M is a tridiagonal matrix,this is easily done in O ( N ) numerical operations in Fourier space. We note thatif one tries to avoid the FFT and iterate Eq. (70) directly in q space, then thecorresponding matrix M on the l.h.s. of Eq. (70) would be a full matrix andeach iteration would require O ( N ) numerical operations.A fixed point of the iteration (77) corresponds to the solution of Eq. (74).The straightforward iteration of (74) (instead of (77)) would diverge becauseof the positive eigenvalue λ = 1 of Corollary 1 for the linearized operator M − L . In contrast, Eq. (74) ensures an approximate projection into the sub-space orthogonal to the corresponding unstable eigenvector f. The originalPetviashvili method [36] is the nonlinear version of Eq. (77) for the particularvalue ∆τ = 1 and is often successful with both partial differential equations(PDEs) (see e.g. Refs. [26,40]) and nonlocal PDEs (see e.g. Ref. [29]). How-ever, the linear operator M − L generally has extra eigenvalues preventingthe convergence of the original Petviashvili method. GPM however uses thefreedom in choice of the parameter ∆τ to achieve convergence even with suchextra eigenvalues, see Refs. [14,26,40] for more discussion.An additional complication that arises in our Eq. (70), compared with thestraightforward use of GPM in general PDEs, is that we do not know α inadvance. Instead, for each value of a there is a nonlinear eigenvalue α ( a ) toEq. (70) that we need to determine. If we use a general value of α , then the ollapse vs. blow up 33 iteration (77) would not converge because the solution of Eq. (70) does notexist for such general values of α. To address this additional complication, we make an initial guess of α = α guess for fixed a and iterate Eq. (70) for α guess . If α guess < α ( a ) then the gen-eralized Petviashvili iteration (after an initial transient) shrinks towards q = 0.If α guess > α ( a ) then the solution expands away from q = 0. We used the bi-section method to determine α ( a ) for a given a . We start from a large enoughinterval [ α L , α R ], so that α ( a ) ∈ [ α L , α R ]. Then we try α guess = ( α L + α R ) / α guess , we obtainthe updated values [ α L , α R ]. These updated values ensure a factor 2 decreaseof the length of the updated interval [ α L , α R ], completing the first step ofthe bisection method. We continue such bisection steps until convergence to α ( a ) (i.e., until the residual of Eq. (70) decreases down to near round-off val-ues and does not decrease anymore). For each updated α guess we use thesolution from the previous bisection step to speed-up the convergence. Wejudged the expansion/shrinking of the solution by tracking the movement ofits maximum point which was determined as a critical point of the function f (cid:48) ( q ) = (cid:80) k = N − k = − N i k ˆ f k e i kq using spectral interpolation and a root-finding algo-rithm. Also, in order to pass over the initial transient dynamics (that dependson the initial guess of the solution) we skip 10 /∆τ − /∆τ initial GPM it-erations before judging the expansion/shrinking of the solution to classify thecurrent α guess . The larger ∆τ we used, the less iterations were needed, but toolarge a ∆τ leads to instability of the algorithm, so we need to keep it under acertain level. For the initial guess of the solution we typically used IC2 fromEq. (62) with V c = 1 / . < a < a c , and N = 64; ∆τ was reduced from0 . a = 0 . − near a c . For a < . ∆τ = 0 . − V c (down to 2 − ) and larger N (up to 2 ) because of the slowlydecaying tails of the function f ( q ) for small a (see the next paragraph). Fig.15 illustrates the convergence of the [ α L , α R ] interval to α ( a ) and convergenceof the residual of Eq. (70) with bisection iterations for a = 0 .
2, starting withan initial condition IC2 in (62) with V c = 1 / ≈ . × − (singularityis at ξ = i V c ) and N = 2 . The converged solution is shown in Fig. 16 (leftpanel) with a closest singularity at a distance ξ c = 7 . · − from the realline in ξ -space and at a distance q c = 1 . · − in q -space.We note that the symmetry (54) implies that ξ c can be stretched by anarbitrary positive constant. The iteration (77) generally converges to differentvalues of ξ c depending on IC (i.e., the zeroth iteration). After that one canrescale any such solution in ξ by any fixed value of ξ c . This rescaling freedomcan also be seen through the existence of the free parameter ˜ v c in the exactsolutions (32) and (38), (39).We computed self-similar profiles f ( ξ ) and g ( ξ ) for various values of a < a c to obtain α ( a ) shown in Table 1 as α e ( a ) . Additionally, we make sure that the f ( ξ ) profile tails scale as in Eq. (52) at ξ → ±∞ and we also fit the g ( ξ ) profiletails to the power law g ( ξ ) ∝ ξ β . (78) Fig. 15
Convergence of the interval [ α L , α R ] to α ( a ) (left panel) and convergence of theresidual of Eq. (76) (right panel) for the iteration (77) with a = 0 . . Here we used IC2 (62)with V c = 2 − ≈ . · − and N = 2 as the zeroth iteration. Fig. 17 show examples of such scaling and fit for a = 0 .
2. Several other curveswith different powers of ξ are present on the graphs for comparison. The fittedvalues of β ( a ) are given in Table 1 and Fig. 16 (right panel). Ignoring forthe moment the Hilbert transform, the integration operator ∂ − ξ involved indetermining g ( ξ ) from f ( ξ ) in Eq. (50) suggests that g ( ξ ) ∝ ξ − α +1 at ξ → ±∞ , (79)which implies that β = − α + 1 . (80)However, the Hilbert transform in Eq. (50) can affect this scaling. We findthat (80) is valid for 0 < a (cid:46) .
4, while a transition to the constant scaling β = − a ≈ .
45 as seen in Table 1 and Fig. 16 (right panel). Inparticular, the exact analytical solution (34) for a = 1 / α = 1 / β = − M as discussed in Section 6.That scaling was confirmed with high precision in our simulations so we donot show it in Table 1. For a < g ( ξ ) has two regions with twodifferent scalings, see Fig. 18 for a = − .
1. While the tail of g ( ξ ) still decays as ξ → ±∞ , there is an intermediate scaling regime which approximately obeys(80) as seen in Fig. 18 (left panel). We are able to observe this intermediatescaling for − . ≤ a <
0. Going below a = − . f ( ξ ) and g ( ξ ) decay very slowly and it requires morethan 10 grid points to achieve good accuracy. For a < β inTable 1 and in Fig. 16 (right panel) are from this intermediate scaling.We estimate that our iteration procedure provides at least 5-8 digits ofprecision of in α ( a ) and 2-3 digits of precision in β ( a ) for a ≥ .
3, when thespectrum of f ( q ) is fully resolved. The values of α ( a ) and β ( a ) were challenging ollapse vs. blow up 35 Fig. 16
Left panel: a = 0 .
2. Functions f ( ξ ) and scaled g ( ξ ) obtained by the iteration (77).Right panel: Power-law of scaling of the tails of g ( ξ ) vs. a . Fig. 17 a = 0 .
2. Left panel: Tail of f ( ξ ) from Fig. 16 (left panel). The dashed line shows thedecay of f ( ξ ) when it is approximated by its leading order singularities alone, as obtainedfrom (17), neglecting the l.s.t. Right panel: Tail of g ( ξ ) from Fig. 16 (left panel) comparedwith different power laws. Fig. 18
Plots of g ( ξ ) for a = − .
1. Left panel: Graph of g ( ξ ) showing two extrema (onemaximum and one minimum) in each half-space of ξ . The inset gives a magnified viewshowing extrema at small ξ . Right panel: Log-log plot of g ( ξ ) for positive ξ . Here g ( ξ ) = 0at ξ ≈ .
41. Solid lines show the scaling (79) and a fit to power law (78).6 P.M. Lushnikov, D.A. Silantyev and M. Siegel to obtain with more than 3-4 and ∼ a (cid:46) . α (cid:38) .
75) and especially for a < α >
1) since we couldnot resolve the Fourier spectrum | ˆ f k | down to round-off level 10 − , even with N = 2 modes. At its root, this is due to the slow decay of f ( ξ ) ∼ | ξ | − /α for | ξ | → ∞ and relatively large α .The numerical values of β in the scaling (78) are important to distinguishbetween solutions with infinite and finite energy E K (10), which as mentionedis of interest in analogy with the question of singularity formation in the 3DEuler and Navier-Stokes equations. Assuming that the solution is close to theself-similar profile (6), changing the variable from x to ξ in (10) and using theself-similar profile (51) of the velocity u ( x, t ) we obtain that E K = E selfsimK + E restK , (81)where E selfsimK = x b (cid:90) − x b u ( x )d x ∼ τ α − ξ b (cid:90) − ξ b g ( ξ )d ξ, ξ b = x b τ α , (82)is the kinetic energy of the approximately self-similar part of the solutionlocated at x ∈ [ − x b , x b ] and E restK is the kinetic energy of the numericalsolution outside of this interval. Here we define the cutoff value x = x b asthe spatial location where the numerical solution deviates from the self-similarprofile (6) by 5%, while inside of the interval [ − x b , x b ] the relative deviation isless than 5%. We determine both the variable ξ by the same type of procedureas in Fig. 4. We find from simulations with a < a c that x b ( t ) ≈ const ∼ τ . (83)Such behaviour is typical for collapsing self-similar solutions, see e.g. Ref. [38,25,16,31]. It implies that ξ b → ∞ as t → t c .There is no qualitative difference between integrals I g,ξ b := ξ b (cid:82) − ξ b g ( ξ )d ξ and I g, ∞ = ∞ (cid:82) −∞ g ( ξ )d ξ provided I g, ∞ < ∞ . The finiteness of I g, ∞ requires that β < − for the scaling of the tails of g ( ξ ) in (78). Using equation (80) we obtainthat β = − implies α = , i.e. β < − for α < . From the interpolationof the data of Table 1 we find that α = corresponds to a = 0 . ± . I g, ∞ < ∞ for a > . ± .
001 and I g, ∞ = ∞ for a < . ± . I g,ξ b is multiplied by τ α − inequation (82). This means that in the limit t → t c and for α < , there is acompetition between the decrease of τ α − and the growth of I g,ξ b as ξ b → ∞ . The scaling (80) for Eq. (78) is valid for a (cid:46) . I g,ξ b ∝ ξ β +1 b = τ − α (2 β +1) x β +1 b for a < . ± . t → t c . Then using Eqs. (80), (82) and (83) we obtain that E selfsimK ∼ ollapse vs. blow up 37 Fig. 19
Growth of the kinetic energy E K over time. Left panel: a = 0 .
2, semi-log plot of E K vs. τ = t c − t shows that E K grows slower than log( τ ) or any power of τ as t → t c .Center panel: a = 0 .
4, verification of the scaling E K ∼ τ α − in (82) with I g, ∞ < ∞ . Rightpanel: a = 1 . E K → ∞ exponentially as t → ∞ . τ ∼ const . Also since the main dynamics is happening in x ∈ [ − x b , x b ] with x b ( t ) ∼ const, we conclude that E restK → const as t → t c , so overall the growthof E K ( t ) as t → t c is very slow (i.e., slower than any power of τ ) for such a where the scaling (80) is true. This result is in excellent agreement with ourdirect calculation of E K ( t ) from time-dependent simulations which shows thatfor a < . ± .
001 the kinetic energy grows more slowly than log( τ ) or anypower of τ as t → t c ; see Fig. 19 (left panel) for a = 0 . . ± . < a ≤ E K → ∞ as t → t c (whilebeing finite for any t < t c ), since α < E K ∼ τ α − → ∞ as t → t c with I g, ∞ < ∞ ; see Fig. 19 (center panel) for a verification of this scaling when a = 0 .
4. For a (cid:38) .
3, which corresponds to an expanding solution with infinitetime singularity, E K → ∞ as t → ∞ , while being finite for any t < ∞ ; see Fig.19 (right panel) for an example with a = 1 .
5. For a > a c , the above splittingof E K into two parts is no longer valid but we nevertheless verify the claimsabove via time-dependent numerical simulation.For some values of a we computed α ( a a = 2 / α ( a ) = 0 . . . . and to compute f ( q ) up to ∼
60 digits of precision, see Fig. 20. High preci-sion computations help validate the results from double precision calculations,and allow us to obtain a good quality analytic continuation of the solution f ( ξ ) = f ( q ( ξ )) from the real line ξ ∈ R to the complex plane ξ ∈ C via theAAA-algorithm [32], see Section 10 below.
10 Analytical continuation into the complex plane by rationalapproximation and structure of singularities
Fits of the Fourier spectrum using Eq. (65) allows us to find only the singular-ity closest to the real line. A more powerful numerical technique of analytical continuation based on rational interpolants [1,15,12,32] allows us to go deeper(further away from the real line) into the complex plane, well beyond the clos-est singularity. However, analytic continuation further from the real line oftenrequires an increase in numerical precision, even well above the standard dou-ble precision [15,12]. In this paper we use a rational interpolation based on amodified version of the AAA algorithm of Ref. [32]. AAA finds an approxima-tion f AAA ( ξ ) to a complex function f ( ξ ) in barycentric form by minimizingthe L error of the approximation on the real line.The barycentric form is given by f AAA ( ξ ) := n ( ξ ) d ( ξ ) = (cid:80) mi =1 w i f i ξ − ξ i (cid:80) mi =1 w i ξ − ξ i , (84)where m ≥ ξ i are a set of real distinct support points , f i are aset of real data values , and w i are a set of real weights determined by L errorminimization. The integer m is increased until the L error between f AAA ( ξ )and f ( ξ ) on the real line is on the level of 10 − P R , where
P R is the currentworking precision. For analytic functions the error decreases exponentially in m . The Barycentric form (84) is a quotient of two polynomials n ( ξ ) and d ( ξ ).A partial fraction expansion of this quotient results in a sum of m − f polesAAA ( ξ ) = (cid:80) m − i =1 a i ξ − b i , with locations b i and residues a i determined by the values of w i and ξ i . The pole locations b i , which are zerosof d ( ξ ), are determined by solving a generalized eigenvalue problem describedin Ref. [32]. The values of the residues a i can be computed using L’Hospital’srule a i = res ( f AAA , b i ) = n ( b i ) /d (cid:48) ( b i ). If our data for an analytic function isgiven with precision P R on the real line, AAA and subsequent computationsof b i approximate the location of single poles with maximum precision ∼ P R ,double poles with precision ∼ P R/
2, and triple poles with precision ∼ P R/ | f ( ξ ) − f polesAAA ( ξ ) | ≈ − P R on the real line in the case of higher order poles if we increase the precision ofintermediate computations in the generalized eigenvalue problem by a factorof two for double poles and a factor of three for triple poles. We additionallymodified the original AAA algorithm [32] to deal with odd and even functionsmore efficiently and output more symmetrical sets of poles.In the particular case a = 2 / , we use 68-digit precision arithemetic for thenumerical solution of f ( ξ ) described at the end of Section 9, and incorporatethis into the AAA algorithm. This method shows that the closest singularitiesto the real line are a pair of the third order poles ∝ / ( ξ ± i χ c ) , in full agree-ment with Theorem 1 (Eq. (21) of Section 2) and the Fourier spectrum analysisof Section 8. The location ξ = ± i χ c (here Re ( χ c ) > Re ( χ c ) (cid:29) | Im ( χ c ) | )and the third order type of these poles are automatically approximated by theAAA algorithm as three simple poles (cid:80) i =1 a i ξ − b i lying very close to each other( | b − b | , | b − b | < . · − ) with the sum of their residues being es-sentially zero ( | (cid:80) i =1 a i | / | a | ≈ . · − ). We define the location of the ollapse vs. blow up 39 Fig. 20
Convergence of the residual (left panel) and spectrum of the solution (right panel)to Eq. (76), computed with a = 2 / V c =1 /
16 = 0 . N = 2048 in the zeroth iteration. Fig. 21
The structure of the complex singularities of the solution from Fig. 20 approximatedby a set of simple poles, f ( ξ ) ≈ f polesAAA ( ξ ) = (cid:80) m − i =1 a i ξ − b i using the AAA algorithm (leftpanel), and the relative error on the real line between the solution f ( ξ ) and its approximation f AAA ( ξ ) (right panel). The simple poles are shown as dots at locations b i with the size ofdot scaled with log | a i | . The branch cuts are approximated as lines connecting the dots.The triple poles locations are ξ ≈ ± i0 . ξ = ξ branch ≈± . ± i0 . triple pole by the average i χ c = (cid:80) i =1 b i / D := (cid:80) i =1 ( b i − i χ c ) a i is negligible, | D | ≈ . · − .In contrast, the quadrupole moment Q := (cid:80) i =1 ( b i − i χ c ) a i is distinct fromzero, | Q | ≈ . · − , so this multipole is well approximated by Q ( ξ − i χ c ) . Thecomplex conjugate point ξ = − i χ c was treated in a similar way, i.e., by anotherset of 3 poles of AAA.We find that the rest of the singularities of f ( ξ ) are branch points withbranch cuts extending from them. AAA approximates branch cuts by sets ofpoles, and Refs. [15,12] demonstrate how to recover branch cuts from this set of poles by increasing the numerical precision. The increase of numerical precisionrequires an increase in the number of poles m in rational interpolants to matchthe precision. These poles, which are located on a branch cut, become moredense with the increase in precision and thus recover the location of the branchcut in the continuous (infinite precision) limit. The main motivation for using68-digit precision in this paper was to ensure that we robustly recover branchcuts, see Fig. 21 (left panel). In the particular case a = 2 / , double precisionallows us to robustly see ∼
30 poles, whereas 68-digit precision allows us to see ∼
150 poles. The number of poles we use for a fixed precision is determinedby the minimal number of AAA poles to match the numerical precision of thesolution on the real line. Increasing the number poles beyond this minimalnumber produces spurious poles with very small residues, which is the analogof the round-off floor in the Fourier spectrum. We note that the exact shape ofthe branch cuts is not fixed analytically – the AAA algorithm simply providesa set of poles that corresponds to the smallest L error on the real axis for thegiven number of poles. Thus, the AAA approximation of the branch cut mightmove with a change of the precision. In contrast, the branch points computedby the algorithm are fixed. One can see 4 branch points in Fig. 21 (left panel),with two branch cuts going upward and coalescing on the imaginary axis andextending further to +i ∞ . Another two branch cuts extend downwards andmerge on the imaginary axis before going off to − i ∞ .Our investigations of complex singularities via AAA approximations showthat for any a , except for a = n − n , n = 1 , , , . . . (which corresponds to theinteger values γ = n in Eq. (21)), there is another pair of vertical branch cutscoming out of ξ = ± i χ c and coalescing with the rest of the branch cuts on theimaginary axis. For a < a c the side branch points are always above the mainsingularity at ξ = ± i χ c and their locations are ξ branch = ± (cid:15) ( a ) χ c ± i(1 + (cid:15) ( a )) χ c , where roughly (cid:15) ( a ) ∼ , (cid:15) ( a ) ∼
1. In particular, Re [ ξ branch ] /χ c < . Im [ ξ branch ] /χ c > a < . Re [ ξ branch ] /χ c ≈ . Im [ ξ branch ] /χ c ≈ .
64 for a = 2 / Re [ ξ branch ] /χ c ≈ . Im [ ξ branch ] /χ c ≈ .
51 near a = a c .
11 Results of time dependent simulations and Petviashviliiterations for periodic BC
Motivated by simulations of the generalized CLM equation (1) in Ref. [34]for 2 π -periodic BC with a = 1, we performed simulations for a wide range ofvalues of the parameter a. For this we used the periodic version of the Hilberttransform H π (58) in Eq. (1) instead of H .Simulations for a < a c show collapsing solutions with α >
0, and differenttypes of IC give qualitatively similar results near the collapse time t = t c asin the real line x ∈ R case with the same α ( a ) (see Table 1). Hence we do notdescribe them here. Expanding solutions for a > a c behaved differently sincethe finite spatial interval [ − π, π ] arrested the increasing width of the solutionat large enough times. Thus we focus our discussion on a > a c and presentdetailed results of our simulations, in particular the cases of a = 0 . a = 1. ollapse vs. blow up 41 Fig. 22
Left panel: Convergence of time-dependent numerical solution of Eqs. (1) and(58) with a = 0 . t = t c = 1 . . . . . The plot is scaled vertically by max x | ω | and horizontally bythe location x max ( t ) of max x | ω | . Right panel: Time dependence of | x max ( t ) | , which showsslowdown and eventual arrest of the horizontal expansion of the solution. We performed a simulation with a = 0 . ω ( x ) = −
43 [sin( x ) + 0 . x )] (85)which is qualitatively similar to the particular case (63) of IC2 (62), with q replaced by x and V c = 1 , T c = 1. After an initial spatial expansion, thesolution is arrested by the periodic boundary conditions. This arrest resultsin the qualitative change of the dynamics, see for example the right panelof Fig. 22 for the time dependence of the location x max ( t ) of max x | ω ( x ) | . Atlater times we still find a finite time blow up of the solution with max x | ω ( x ) | and max x | u ( x ) | → ∞ as t → t c . However, instead of Eq. (6), the solutionconverges to a new universal self-similar blow-up profile given by Eq. (9), asdemonstrated in left panel of Fig. 22. A comparison of Eqs. (6) and (9) revealsthat we can formally obtain Eq. (9) by setting α = 0 in Eq. (6) (although Eq.(9) has periodic boundary conditions, vs. decaying BC of Eq. (6)). We notethat taking the limit a → a − c in Eq. (6), we also obtain α = 0 . However, itremains unknown if Eq. (9) can be obtained from the continuation of Eq. (6)across a = a c .The spectrum ˆ ω k is initially exponentially decaying but expands and be-comes mostly algebraically decaying (similar to Fig. 14). Finite precision arith-metic only “sees” algebraic decay | ˆ ω k ( x ) | ∼ k − when t is close enough to t c ,see Fig. 24. This is because of a jump in ω xx forming at x = ± π , see Fig. 23(left and middle panels). Due to the spectrum being initially oscillatory it wasdifficult to accurately extract values of δ ( t ) and p ( t ) from a fit to Eq. (65),but using a nonoscillatory spectrum which emerges later in the simulation wewere able to recover some data for δ ( t ) and p ( t ) as shown in Fig. 24. There,one can see that δ ( t ) → p ( t ) → t → t c . Fig. 23
Results of the simulation of Eqs. (1) and (58) with a = 0 . ω ( x, t ), its derivatives ω x ( x, t ), ω xx ( x, t ), and u ( x, t ) at t = 1 . . Right panels:the growth of maximum values of the corresponding quantities over time.
Fig. 24
Left panel: Log-log plot of the Fourier spectrum | ˆ ω k | from Fig. 23 at t = 1 . a = 0 .
8. The red line represents a fit to the model (65) with green line showing a portionof the | ˆ ω k | used for the fit. Center and right panels: δ ( t ) and p ( t ) obtained from the fit of | ˆ ω k | to Eq. (65) at different times. Red lines in the center panel also show a fit to the model δ ( t ) ∝ ( t c − t ) α .ollapse vs. blow up 43 Fig. 25
Results of the simulation of Eqs. (1) and (58) with a = 1 and IC (85) showing ω ( x, t ), its derivatives ω x ( x, t ), ω xx ( x, t ), and u ( x, t ) at t = 2 . For a = 1 we considered two different types of ICs. The first one is IC (85),for which we observe global existence of the solution. Initially the amplitudeof the solution ω ( x ) grows in time, similar to the infinite domain case. But thisgrowth slows down at later times and eventually reaches a plateau with thethe same behaviour in u ( x ), see Fig. 25. Also max x | ω x | = | ω x ( x = 0) | remainsnearly constant throughout the simulation. We observe unbounded growth of | ω xx | near x = ± π that appears to be exponential in time. Due to the spectrumbeing oscillatory it was difficult to accurately extract values of δ ( t ) and p ( t )from a fit to Eq. (65). However, using AAA rational approximation we wereable to observe two pairs of branch cuts approach the real line near x = ± π as t → ∞ . Replacing IC (85) by the more general IC2 (62) (with q replaced by x and V c , T c = 1) is found to only alter the transient dynamics of the expandingsolution without qualitatively changing the overall behavior.The second type of IC we used for a = 1 is given by ω ( x ) = sin( x ) + 0 . x ) , (86) Fig. 26
Results of the simulation of Eqs. (1) and (58) with a = 1 and IC (86) as in Ref.[34] showing ω ( x ), its derivatives ω x ( x ), ω xx ( x ), and u ( x ) at t ≈
12 and the growth of theirmaximum values as functions of time.
Fig. 27
Left panel: Log-log plot of the Fourier spectrum | ˆ ω k | for the solution in Fig. 26 anda fit to the model (65). Center and right panels: Time dependence of δ ( t ) and p ( t ) obtainedfrom the fit to (65). Center panel also shows a fit of δ ( t ) to the stretched exponential model δ ( t ) ∼ e − κt ν .ollapse vs. blow up 45 which is the same as in Ref. [34]. It allows us to directly compare the resultsof our simulations with Ref. [34]. We obtain exactly the same plots as in Fig.1 of Ref. [34], see Fig. 26. The difference between simulations with IC (85)and IC (86) are seen by comparing Figs. 25 and 26. For example, the spatialderivatives of ω approach discontinuities at x = 0 in Fig. 25 vs. x = ± π inFig. 26. Also, | ω xx | → ∞ grows exponentially in time in both cases althoughat different locations in x . The exponentially decaying spectrum | ˆ ω k | with IC(86) widens over time and becomes mostly algebraically decaying, see Fig. 27(left panel). In this case the only singularity is near x = 0; the AAA rationalapproximation shows an approach of two vertical branch cuts to x = 0 overtime, so the spectrum is not oscillatory and we are able to easily recover δ ( t )and p ( t ) from the fit to Eq. (65). The fits show a stretched-exponential intime approach of the singularity to the real line i.e., δ ( t ) ∼ e − κt ν , see Fig. 27(middle panel). Figure 27 (middle and right panels) showing δ ( t ) and p ( t ) canbe compared with Fig. 3(a,b) of Ref. [34]. Our values of δ ( t ) match those valuesfrom Fig. 3(a) of Ref. [34] well, while values of p ( t ) do not match precisely withFig. 3(b) of Ref. [34] because they marginally depend on the particular partof spectrum | ˆ ω k | that is used for the fitting.For a > x -space is arrested by the periodic boundary conditionswith an infinite slope forming at the boundary x = ± π so that max x | ω x | → ∞ as t → ∞ (although max x | ω | , max x | u | , | ω x ( x = 0) | → t → ∞ ). The complexsingularities approach the real line in infinite time. Their positions scale like x sing ∼ ± π ± i y exp ( − κ t ν ), where y , κ , ν >
0. When a → + , we observethat max x | ω | grows for a short time and then decays. Unlike the x ∈ R case, itis relatively easy to compute accurately for a → + and we have been able toobtain numerical evidence of global existence for a as small as 1.000001. ForIC (86), we also observe global existence of the solution with decay of max x | ω | and unbounded growth of | ω x ( x = 0) | as t → ∞ . The complex singularitiesapproach the real line like x sing ∼ ± i y exp ( − κ t ν ), where y , κ , ν > x ∈ R case described in Section 9 for a ≤ .
95, while for a = 1 we have that E K → const as t → ∞ (because max x | u | → const as t → ∞ ) and for a > E K → t → ∞ (because max x | u | → t → ∞ ). Self-similar profiles from GPM . We also numerically computed the self-similar profile f ( x ) in Eq. (9) for a c < a ≤ .
85 using GPM described inSection 9 with α = 0. In contrast to Section 9, we do not need to use thecoordinate transformation (69) because f ( x ) is now 2 π -periodic with ξ ≡ x .We used GPM to solve Eq. (49) by the iteration (77) with M f and N f fromEq. (70) replaced by M f := f = − agf x + f H π f := N f,g x = H π f. (87) Fig. 28
The Fourier spectrum | ˆ ω k | of the self-similar profile (9) for a = 0 .
71 obtained byGPM iterations (77) of Eq. (87). Two fits are shown in different ranges of k with the first fitto Eq. (65) with δ (cid:54) = 0 at intermediate k and the second a power law fit ∝ | k | − p b for larger | k | . Left panel: Log-linear plot where the first fit turns into a nearly linear function. Rightpanel: Log-Log plot where the second fit turns into a nearly linear function. Fig. 29
Left panel: The Fourier spectra | ˆ ω k | of the self-similar profile (9) for various valuesof a as in the Table 2 obtained by GPM iterations (77) of Eq. (87). Right panel: p ( a ) and p b ( a ) from the Table 2 extracted from the two fits as in Fig. 28. The matrix M used in Eq. (77) now turns into the identity matrix. We donot need to solve the nonlinear eigenvalue problem for α because now α ≡ ∆τ even more thanin Section 9 to make sure the iterations converged and also had to use moreFourier modes in the spectrum, since the spectrum decay is only algebraic forthese solutions. Due to these technical limitations we were unable to explorethe range 0 . < a <
1, but we fully expect that self-similar solutions existthere because time-dependent simulations converge to self-similar profiles, atleast over the lower range a c < a (cid:46) .
95 (see Fig. 22). It was not possibleto obtain convergence in the upper range 0 . (cid:46) a < ollapse vs. blow up 47 before the computation became prohibitively slow. The behavior of solutions(blow up vs. global existence) therefore remains unknown in this range. Weconjecture that blow up occurs for all a c < a < a = 1 (as demonstrated) and for larger values of a .The Fourier spectrum of | ˆ ω k | corresponding to the self-similar profile (9)has two distinct domains for | k | (cid:29)
1. The particular case a = 0 .
71 shown inFig. 28 depicts such domains. The first domain corresponds to complex singu-larities of Theorem 1 (Eq. (21)) located at x sing = ± i δ. This domain is wellfitted by Eq. (65). From this fit we find that δ = 1 . p = − . p = − a − a = − . . . . which agrees within an accuracyof < .
05% with the numerical fit to Eq. (65). The second domain is due tocomplex singularities located at x = ± π and results in a discontinuity of high-order derivatives of ω ( x ) at the periodic boundary. This domain has the powerlaw spectrum ∝ | k | − p b (i.e., in Eq. (65) it corresponds to δ = 0 and p = p b )which is dominant for larger | k | . In the particular case of Fig. 28, we obtain p b = 9 . . . . . This implies that the 9th and higher-order derivatives of ω ( x ) have a discontinuity at the periodic boundary. All these singularities canbe seen using the AAA algorithm described in Section 10. We also find thatas a approaches to a c from the right, i.e. a → a + c , increasingly higher orderderivatives experience discontinuities at the periodic boundary, i.e. p b → ∞ as a → a + c , see Fig. 29 (right panel). These solutions with finite smoothnessat the periodic boundary can be considered the analog of the self-similar so-lutions with compact support found in Sections 8 and 9, for solutions on thereal line with a c < a ≤ δ, p and p b for various values of parameter a obtained from the fits described above. We note that the symmetry (54) is notvalid for periodic BC. Thus, the parameter δ is now fixed for each a , contraryto the case x ∈ R where it is a free parameter, cf. Section 9.Here we summarize the solution behaviour of Eqs. (1) and (58) for x ∈ [ − π, π ] and generic smooth IC depending on the parameter a : – a < a c : Behaviour of solutions is the same at t → t c as for the x ∈ R case,with collapse as in Eq. (6). – a c < a (cid:46) .
95 : Blow up both in ω and u in finite time t c with solutionapproaching the universal self-similar profile (9) as t → t c . That profile f ( x )has discontinuities in the high-order derivatives with complex singularitiestouching the real line only at x = ± π . The number of continuous derivativesbecomes infinite in the limit a → a + c . The singularities approach the realline as x sing (cid:39) ± π ± i( t c − t ) α y b , where α ( a ) > – a = 1 : Global existence of solution with a singularity approaching thereal line exponentially in time. For both IC (85) and IC (86) we findmax x | ω | , max x | u | , max x | u x | → const , max x | ω x | = | ω x ( x = 0) | = const , andmax x | ω xx | → ∞ as t → ∞ . – a > line near x = ± π and max x | ω | , max x | u | , | ω x ( x = 0) | → | ω x | → ∞ as t → ∞ .
12 Conclusions and discussion
We have performed a systematic sweep of the parameter a in the generalizedCLM equation (1) to determine the possibility of singularity formation and,when it occurs, its type, i.e., collapse vs. blow-up. We identified a new criticalvalue a = a c = 0 . . . . such that for a < a c collapse occursboth on the real line x ∈ R and for periodic BC. Here, collapse means that notonly is there a finite time singularity in which the amplitude of the solution ω ( x, t ) tends to infinity, but there is also a catastrophic shrinking of the spatialextent of the solution to zero as t → t c , described by the self-similar form (6).In the intermediate range a c < a ≤ , we found there is finite-time singularityformation for x ∈ R , with the self-similar solution (6) experiencing an infiniterate of expansion as t → t c . This type of self-similar singularity formation, inwhich the spatial domain does not collapse, is termed ‘blow up’. The power α in Eq. (6) controls collapse (for α > , a < a c ) vs. blow up ( α ≤ , a ≥ a c ). Weelucidated the dependence of α ( a ) on a via both direct numerical simulationof Eq. (1) and the solution of a nonlinear eigenvalue problem (49) using thegeneralized Petviashvili method (77). We have also performed multiprecisionsimulations (up to 68 digits of accuracy) to demonstrate the possibility ofrecovering α ( a ) and the structure of self-similar solutions with any desiredprecision.We also have found a new analytic collapsing solution of Eq. (1) for a = 1 / α = 1 / . This solution has a finite energy E K (10) forany time t < t c . Such finite energy solutions are of interest in analogy withthe problem concerning global regularity of the 3D Euler and Navier-Stokesequations with smooth initial data, see Refs. [18,19]. We found for generalvalues of a that the self-similar solution (6) is real analytic for a < a c while ithas finite support for a c < a ≤ . We identified that the blow up for periodic BC with a c < a ≤ .
95 isqualitatively different from that for x ∈ R , because the periodic BC arrestsor blocks the unbounded spatial expansion of the solution on the real line.To our surprise, such arrest does not result in the global existence of thesolution but instead leads to a new form of self-similar blow-up (9), in whichweak singularities develop at the boundaries of the periodic domain. In thelimit a → a + c , this self-similar solution turns into an infinitely smooth ( C ∞ )solution. We believe that the qualitative difference in blow up between x ∈ R and periodic BC might serve as an interesting lesson relevant to the search forsingularities in the 3D Euler equation.Both self-similar solutions (6) and (9) are nonlinearly stable, as follows fromour simulations. Quite generic classes of IC converge to these solutions duringthe temporal evolution. In the case of Eq. (6), such convergence/stability is ollapse vs. blow up 49 understood in the sense of convergence to a family of self-similar solutions, upto a rescaling in x , because of the symmetry (54) of Eq. (49).The structure of the leading order singularities in the complex plane x (which is the analytical continuation from x ∈ R ) is determined by Theo-rem 1. That result is valid for both x ∈ R and periodic BC and is in fullagreement with simulations. For a < a c the leading order singularities are theclosest singularities to the real line in the complex x -plane. For a > a c , thesesingularities still determine the structure of self-similar solutions near x = 0,while the solution near the boundaries of finite support in x ∈ R and theperiodic boundaries for periodic BC are controlled by less singular terms. Theself-similar solution profiles for these a have been found with high accuracyby solving a nonlinear eigenvalue problem. We have also proved in Theorem3 that, except for the exact closed-form solutions for a = 0 and a = 1 /
2, theanalytical structure of singularities in the complex x -plane goes beyond theleading order singularities. In particular, we numerically identified using theAAA algorithm the existence of additional, non-leading-order branch pointsfor a (cid:54) = 0 , / a (cid:38) . x ∈ R , while for periodic BC globalexistence is ensured for a ≥ . In the remaining gaps 1 < a (cid:46) . x ∈ R and 0 . < a < a will require additionalanalysis and/or substantial efforts in simulation.We suggest that among many other issues, the following questions wouldbe interesting to address in future work:1. Analytical study of the complex singularities beyond the leading ordersingularities addressed in Theorem 1. In particular, the case a = 2 / a > a c for x ∈ R ,or use a version of the method in Ref. [6] based on cubic splines. However,splines generally loose information about the analyticity of solutions in thecomplex plane. One way to improve the performance of GPM in this range of a might be to use a coordinate transform in the form of a conformal mappingwhich would simultaneously resolve the numerical grid near x = ± x b whilekeeping the analyticity of the solution intact. This type of approach has beensuggested in Ref. [30].3. Fill the gaps in our knowledge on blow up vs. global existence of solutionsin the parameter regime 1 < a (cid:46) . x ∈ R and 0 . < a < a = a c betweenself-similar solutions (6) for the case x ∈ R and Eq. (9) for periodic BC.5. Perform an analysis of the nonlinear stability of the blow-up solutions.This could be qualitatively similar to the stability of collapse in PDEs such as the nonlinear Schr¨odinger equation and the Patlak-Keller-Segel equation, seee.g. Refs. [41,7,38,3,25,31].6. Analyze the formation of singularities at the initial time t = 0 + . Thiscan give information on the type of singularities which first form in the com-plex plane, and subsequently move toward the real line. Such an analysis hasbeen previously performed for the evolution of a vortex sheet in the Kelvin-Helmholtz problem [10], which is also governed by a nonlocal PDE. However, asignificant difference between the current problem and the vortex sheet prob-lem is that here the singularities initially form at infinity in the complex plane,whereas in the vortex sheet problem they are generated at finite locations, dueto a singularities in the kernel of the nonlocal term at these locations. Conflict of interest
The authors declare that they have no conflict of interest.
Acknowledgements
The work of P.M.L. was supported by the National Science Foun-dation, grant DMS-1814619. M.S. was supported by National Science Foundation grantDMS-1909407. P.M.L. the support of the Russian Ministry of Science and Higher EducationGrant No. 075-15-2019-1893. Simulations were performed at the Texas Advanced Comput-ing Center using the Extreme Science and Engineering Discovery Environment (XSEDE),supported by NSF Grant ACI-1053575.
A Hilbert transform for transformed variable
In this Appendix we derive the expression for the Hilbert transform in the auxiliary variable q (55) of Section 7.The change of variable (55) in Eq. (2) together with (56) results in H f ( x ) = 1 π p.v. (cid:90) ∞−∞ f ( x (cid:48) ) x − x (cid:48) d x (cid:48) = 1 π p.v. (cid:90) π − π ˜ f ( q (cid:48) )tan q − tan q (cid:48) d q (cid:48) q (cid:48) = 12 π p.v. (cid:90) π − π ˜ f ( q (cid:48) ) (cid:104) q tan q (cid:48) − tan q (cid:48) (cid:16) tan q − tan q (cid:48) (cid:17)(cid:105) tan q − tan q (cid:48) d q (cid:48) = 12 π p.v. (cid:90) π − π ˜ f ( q (cid:48) )tan (cid:16) q − q (cid:48) (cid:17) d q (cid:48) − π π (cid:90) − π ˜ f ( q (cid:48) ) tan q (cid:48) q (cid:48) = H π f ( q ) + C πf , (88)where we used the identitiestan ( a − b ) = tan a − tan b a tan b and 1cos q = tan q q →± π [ H π f ( q ) + C πf ] = 0 . Also H π f ( x ), Eq. (58), is the reduction of H f ( x ), Eq. (2), to the class of 2 π -periodicfunctions. Assuming that f ( x ) is the periodic function with the period 2 π , we obtain fromollapse vs. blow up 51Eq. (2) that H f ( x ) = 1 π ∞ (cid:88) n = −∞ p.v. (cid:90) π − π f ( x (cid:48) ) x − x (cid:48) + 2 πn d x (cid:48) = 12 π p.v. (cid:90) π − π f ( x (cid:48) )tan (cid:16) x − x (cid:48) (cid:17) d x (cid:48) =: H π f ( x ) , (89)where we used the definition (58) and the identity ∞ (cid:88) n = −∞ x + 2 πn = 12 tan x . (90) References
1. Alpert, B., Greengard, L., Hagstrom, T.: Rapid evaluation of nonreflecting boundarykernels for time-domain wave propagation. SIAM J. Num. Anal. , 1138–1164 (2000)2. Baker, G., Caflisch, R.E., Siegel, M.: Singularity formation during Rayleigh–Taylor in-stability. Journal of Fluid Mechanics , 51–78 (1993)3. Brenner, M.P., Constantin, P., Kadanoff, L.P., Schenkel, A., Venkataramani, S.C.: Dif-fusion, attraction and collapse. Nonlinearity (4), 1071–1098 (1999)4. Carrier, G.F., Krook, M., Pearson, C.E.: Functions of a Complex Variable. McGraw-Hill,New York (1966)5. Castro, A., C´ordoba, D.: Infinite energy solutions of the surface quasi-geostrophic equa-tion. Adv. Math. (4), 1820–1829 (2010)6. Chen, J., Hou, T.Y., Huang, D.: On the finite time blowup of the De Gregorio modelfor the 3D Euler equation. arXiv:1905.06387 (2019)7. Childress, S., J.K., P.: Nonlinear aspect of chemotaxis , 217–237 (1981)8. Constantin, P., Lax, P.D., Majda, A.: A simple one-dimensional model for the three-dimensional vorticity equation. Commun. Pure Appl. Math. (6), 715–724 (1985)9. Cooper, G.J., Verner, J.H.: Some Explicit Runge-Kutta Methods of High Order. SIAMJournal on Numerical Analysis (3), 389–405 (1972). DOI 10.1137/0709037. URL https://doi.org/10.1137/0709037
10. Cowley, S.J., Baker, G.R., Tanveer, S.: On the formation of Moore curvature singularitiesin vortex sheets. J. Fluid Mech. , 233–267 (1999)11. De Gregorio, S.: On a one-dimensional model for the three-dimensional vorticity equa-tion. J. Stat. Phys. (5-6), 1251–1263 (1990)12. Dyachenko, A.I., Dyachenko, S.A., Lushnikov, P.M., Zakharov, V.E.: Dynamics of Polesin 2D Hydrodynamics with Free Surface: New Constants of Motion. Journal of FluidMechanics , 891–925 (2019)13. Dyachenko, S.A., Lushnikov, P.M., Korotkevich, A.O.: The complex singularity of aStokes wave. JETP Letters (11), 767–771 (2013). DOI 10.7868/S0370274X1323007014. Dyachenko, S.A., Lushnikov, P.M., Korotkevich, A.O.: The complex singularity of aStokes wave. JETP Letters (11), 675–679 (2013). DOI 10.7868/S0370274X1323007015. Dyachenko, S.A., Lushnikov, P.M., Korotkevich, A.O.: Branch Cuts of Stokes Wave onDeep Water. Part I: Numerical Solution and Pad´e Approximation. Studies in AppliedMathematics , 419–472 (2016). DOI DOI:10.1111/sapm.1212816. Dyachenko, S.A., Lushnikov, P.M., Vladimirova, N.: Logarithmic scaling of the collapsein the critical Keller-Segel equation. Nonlinearity , 3011–3041 (2013)17. Elgindi, T.M., Jeong, I.J.: On the effects of advection and vortex stretching.Archive for Rational Mechanics and Analysis , 1763–1817 (2020). DOI 10.1007/s00205-019-01455-9. URL https://doi.org/10.1007/s00205-019-01455-9
18. Fefferman, C.L.: Existence and smoothness of the Navier-Stokes equation. The millen-nium prize problems pp. 57–67 (2006)19. Gibbon, J.D.: The three-dimensional Euler equations: Where do we stand? Physica D , 1894–1904 (2008)20. Hou, T.Y., Jin, T., Liu, P.: Potential singularity for a family of models of the axisym-metric incompressible flow. J. Nonlinear Sci. (6), 2217–2247 (2018)2 P.M. Lushnikov, D.A. Silantyev and M. Siegel21. Hou, T.Y., Lei, Z., Luo, G., Wang, S., Zou, C.: On finite time singularity and globalregularity of an axisymmetric model for the 3D Euler equations. Archive for RationalMechanics and Analysis (2), 683–706 (2014)22. Hou, T.Y., Li, C.: Dynamic stability of the three-dimensional axisymmetric Navier-Stokes equations with swirl. Commun. Pure Appl. Math. (5), 661–697 (2008)23. Hou, T.Y., Li, R.: Dynamic depletion of vortex stretching and non-blowup of the 3-Dincompressible Euler equations. J. Nonlinear Sci. (6), 639–664 (2006)24. Hou, T.Y., Shi, Z., Wang, S.: On singularity formation of a 3D model for incompressibleNavier-Stokes equations. Adv. Math. (2), 607–641 (2012)25. Kuznetsov, E.A., Zakharov, V.E.: Wave Collapse. World Scientific Publishing Company,New York (2007)26. Lakoba, T.I., Yang, J.: A generalized Petviashvili iteration method for scalar and vectorHamiltonian equations with arbitrary form of nonlinearity. J. Comput. Phys. , 1668–1692 (2007)27. Lei, Z., Hou, T.Y.: On the stabilizing effect of convection in three-dimensional incom-pressible flows. Commun. Pure Appl. Math. (4), 501–564 (2009)28. Lei, Z., Liu, J., Ren, X.: On the Constantin–Lax–Majda model with convection. Comm.Math. Phys. pp. 1–19 (2019)29. Lushnikov, P.M.: Dispersion-managed soliton in a strong dispersion map limit. Opt.Lett. , 1535 – 1537 (2001)30. Lushnikov, P.M., Dyachenko, S.A., Silantyev, D.A.: New conformal mapping for adap-tive resolving of the complex singularities of Stokes wave. Proc. Roy. Soc. A ,20170198 (2017)31. Lushnikov, P.M., Dyachenko, S.A., Vladimirova, N.: Beyond leading-order logarithmicscaling in the catastrophic self-focusing of a laser beam in Kerr media. Phys. Rev. A , 013845 (2013)32. Nakatsukasa, Y., S`ete, O., Trefethen, L.N.: The AAA algorithm for rational approxi-mation. SIAM J. Sci. Comp. (3), A1494–A1522 (2018)33. Okamoto, H., Ohkitani, K.: On the role of the convection term in the equations ofmotion of incompressible fluid. J. Phys. Soc. Japan (10), 2737–2742 (2005)34. Okamoto, H., Sakajo, T., Wunsch, M.: On a generalization of the Constantin–Lax–Majda equation. Nonlinearity (10), 2447 (2008)35. Pelinovsky, D., Stepanyants, Y.: Convergence of petviashvili’s iteration method for nu-merical approximation of stationary solutions of nonlinear wave equations. SIAM J.Numer. Anal. , 1110–1127 (2004)36. Petviashvili, V.I.: Equation for an extraordinary soliton. Sov. J. Plasma Phys. , 257–258 (1976)37. Stein, E.M.: Singular integrals and differentiability properties of functions, vol. 2. Prince-ton university press (1970)38. Sulem, C., Sulem, P.L.: Nonlinear Schr¨odinger Equations: Self-Focusing and Wave Col-lapse. World Scientific, New York (1999)39. Sulem, C., Sulem, P.L., Frisch, H.: Tracing complex singularities with spectral methods.J. Comput. Phys. , 138–161 (1983)40. Yang, J.: Nonlinear Waves in Integrable and Nonintegrable Systems. SIAM (2010)41. Zakharov, V.E.: Collapse of langmuir waves. Sov. Phys. JETP , 908 (1972)ollapse vs. blow up 53 Table 1
Table of values of α , p and α extracted via fits to δ x ( t ), | ˆ ω k | and max | ω x ( x, t ) | intime-dependent simulations of Eqs. (57)-(58) for various values of a . Also shown are valuesof α e and β obtained from eigenvalue problem simulations of Eqs. (70) and (58) describedin Section 9. Accuracy of α ( a ) (for − ≤ a ≤ . α ( a ) (for − ≤ a ≤ . α e ( a ) is about 3-4 digits of precisionfor a < . a ≥ .
3, with more precision for 0 . ≤ a ≤ . a α e β p α α -5 - - 0.855 7.495 7.517-2 - - 0.680 3.444 3.422-1 - - 0.505 2.208 2.206-0.5 - - 0.335120 1.603747 1.600222-0.25 1.296593455 - 0.200942 1.303708 1.302424-0.2 1.239824952 - 0.167139 1.243558 1.242436-0.15 1.181358555 0.133308 0.130811 1.183300 1.182701-0.1 1.121312899 0.100401 0.091110 1.122630 1.122093-0.05 1.061051829 0.060633 0.047696 1.061617 1.0613340 1 0 0.004 1.000243 1.0000190.05 0.938365701 -0.070205 -0.052759 0.938381 0.9382880.1 0.876129662 -0.136336 -0.111326 0.876329 0.8763090.15 0.813179991 -0.240380 -0.176727 0.813219 0.8132150.2 0.749369952 -0.338799 -0.250265 0.749519 0.7495490.25 0.684513621 -0.460507 -0.333582 0.684650 0.6846710.265 0.664818990 -0.500444 -0.360765 0.664827 0.6648300.3 0.618374677 -0.610349 -0.428762 0.618375 0.6183770.35 0.550648498 -0.787978 -0.538583 0.550661 0.5506550.4 0.480939257 -0.939823 -0.666732 0.4809431 0.48094290.425 0.445184823 -0.97452 -0.739156 0.4451863 0.44518600.4375 0.427049782 -0.993899 -0.777804 0.4270512 0.42705080.45 0.408728507 -1 -0.818193 0.40872820 0.408728380.5 0.333333333 -1 -1.0000007 0.33333354 0.333333400.55 0.253852136994 -1 -1.222218 0.25385226 0.253852130.6 0.169098936470 -1 -1.4999991 0.16909915 0.16909893670.65 0.077532635626630 -1 -1.857141 0.07753269 0.077532635626622/3 0.045170944220367 -1 -1.999997 0.04517096 0.045170944220350.68 0.018526534283004 -1 -2.125013 0.01852675 0.018526534282700.685 0.008351682345844 -1 -2.175083 0.00835210 0.0083516823458430.689 0.000137203824593 -1 -2.219165 0.00013724 0.0001372038246030.68905 3.409705703117e-05 -1 -2.221589 3.4145e-05 3.4097057039e-050.68906 1.347443362884e-05 -1 -2.220924 1.3418e-05 1.3474433654e-050.689066 1.10065641e-06 -1 -2.221505 1.0808e-06 1.1006564176e-060.6890665 6.950143e-08 -1 -2.223142 - 6.9501438524e-080.68906653 7.632094e-09 -1 -2.222128 - 7.6321058379e-090.689066533 1.445152e-09 -1 -2.220519 - 1.4451679770e-090.6890665335 4.13992e-10 -1 -2.205923 - 4.1401557848e-100.6890665337 1.537e-12 -1 -2.220897 - 1.5519e-120.6890665337007 9.43093e-14 -1 -2.227272 - 1.1097e-130.68906653370074 1.18169e-14 -1 -2.222533 - 2.7574e-140.689066533700745 1.505397e-15 -1 -2.221208 - 1.4711e-140.6890665337007457 6.169686e-17 -1 - - -0.7 - - - - -0.022810.75 - - - - -0.134350.8 - - - - -0.260080.85 - - - - -0.403840.9 - - - - -0.571180.95 - - - - -0.766431 - - - - -1.000000056 Table 2
Table of values of δ, p and p b extracted via a fit of spectra | ˆ ω k | to the model (65),obtained from eigenvalue problem simulations of Eq. (87) for various values of a , a c < a < δ and p are extracted from the fit | ˆ ω k | ∝ exp( − δ | k | ) / | k | p to the central part ( k ∼
0) of thespectrum and p b is extracted from the fit | ˆ ω k | ∝ / | k | p b in the tails ( k (cid:29)
1) of the spectrum.Simulations with a ≥ .
71 were performed in double precision arithmetic. To see the powerlaw tail of the spectrum and extract p b in the case of a = 0 . a c < a ≤ .
695 the power law tail was not observable even in quadrupleprecision. See Fig. 29 for the spectra and plots of p ( a ) and p b ( a ). The accuracy of δ, p and p b approximately corresponds to the number of digits provided in the table. a δ p p bb