Spatio-temporal linear instability analysis for arbitrary dispersion relations on the Lefschetz thimble in multidimensional spacetime

In linear stability analysis of field quantities described by partial differential equations, the well-established classical theory is all but impossible to apply to concrete problems in its entirety even for uniform backgrounds when the spatial dimension is more than 1. In this study, using the Lefschetz thimble method, we develop a new formalism to give an explicit expression to the asymptotic behavior of linear perturbations. It is not only more mathematically rigorous than the previous theory but also useful practically in its applications to realistic problems, and, as such, has an impact on broad subjects in physics.


I. INTRODUCTION
In many dynamical problems of complex nonlinear systems, it is often important to investigate the perturbative dynamics near fixed points, or time-independent solutions, by linearization, since the original nonlinear equations are too difficult to treat directly. For systems with a finite degree of freedom, it suffices to obtain normalmode solutions with an exponential time-dependence, exp(−iωt), and see if there is a complex ω with a positive imaginary part. In the case of infinite degrees of freedom, however, this may not be the case. In particular, the linear analysis of the dynamics of field quantities defined in an unlimited spatial domain is much more involved, since the single normal mode has a zero-measure and we need to study the behavior of wave packets; the concept of "instability" itself needs to be extended a bit: solutions that grow exponentially in time at every spatial point, i.e., unstable in the ordinary sense, are called "absolutely unstable" whereas those solutions that are damped in time at each fixed point in space but are amplified exponentially if seen from an appropriate moving observer are referred to as "convectively unstable." It is then insufficient to see just the sign of the imaginary parts of frequency ω(k), this time a function of wave number k, for each normal mode.
All these things are well known, e.g., in plasma physics and a mathematically rigorous theory was established by Briggs [1] and others already in 1960's. It is given in some standard text books [2]. In this theory, the asymptotic behavior of a perturbation given as a wave packet is derived by finding a coalescence of some of the complex roots of the dispersion relation (DR) ∆(ω, k) = 0. Since they treated only the case with one spatial dimension, the wave number is not a vector but a number actually. They employed the Laplace transform and its inverse to obtain a solution in the form of complex integrals over ω and k.
The dominant contribution to these integrals at asymptotically large times comes from a pole in the upperhemisphere of the complex ω-plane, which is in turn produced by the coalescence of two roots in k of the DR and the resultant pinch of the integration contour in k at this complex value of ω. The point, at which two complex roots in k of the DR merge, is referred to as the critical point and in addition to ∆(ω, k) = 0, ∂∆(ω, k)/∂k = 0 is also satisfied there. Note that not all the critical points give the pinch of the integration contour. The theory has been applied to shear flows [3,4], jets [5], solitons [6] and even neutrino oscillations recently [7,8].
Brevdo extended this theory to 2 spatial dimensions in 1991 [9]. Although it is indeed a straightforward extension of the Briggs' theory, which just looks for root mergers repeatedly, carrying it out is actually much more difficult than in the one spatial dimension and it is essentially impossible to apply unless the DR is very simple. It is not difficult to imagine then that it is hopeless to extend the theory to higher dimensions. Interestingly, Bers et al. [10] treated the absolute/convective instabilities of electron beams in 3 spatial dimensions (i.e., in 4-dimensional spacetime,) relativistically in 1984 with the same method, seemingly an impossible task. Indeed they fell short of full application of the theory: the authors found a critical point but did not show that the integration contour is really pinched by the merging roots of the DR. This is the core of the problem. Finding critical points is an easy part of the Briggs' theory. What is most difficult is to show that a particular critical point satisfies the pinch criterion and contributes to the asymptotic behavior of the perturbation. This is because it requires knowledge on the behavior of the DR not only in the vicinity of the critical point but also in much wider a region. In the case of one spatial dimension, some people change the order of integrations and conduct the ω-integral first and then apply the formula of the steepest descent method to the remaining k-integral rather blindly, i.e., not showing that the integral contour can be modified so that it could coincide with the steepest-decent path from the particular critical point [11]. A disguise notwithstanding, it is apparent that this alternative formalism has actually the same problem as the partial application of Briggs' theory does. The crucial fact is that the critical point does not always contribute to the asymptotic behavior.
In this paper, we present a new formalism that cir-arXiv:1912.11177v1 [math-ph] 24 Dec 2019 cumvents the difficulty just mentioned and is applicable to arbitrary dimensions and to complicated DR's by utilizing the Lefschetz thimble method, which is based on algebraic topology and geometry and gives a generalization of the steepest descent method to complex manifolds [12]. It was Witten who advocated the Lefschetz thimble method in the analysis of Chern-Simons theory [13]. The method has also been applied to the sign problem in quantum chromodynamics [14]. Our new method, also being based on the Lefschetz thimble method, can not only tell which critical points contribute to the asymptotic behavior of perturbation unambiguously but also be applied easily even when ω(k) is multi-valued, since we no longer need to solve DR ∆(ω, k) = 0 as ω = ω(k).

A. Formulation
We consider the following system of linear partial differential equations D(i∂)S(x) = 0 in (d + 1)-dimensional spacetime, where (x µ ) ≡ (t, x) are the coordinates and the corresponding derivative operators with upper indices are denoted by (∂ µ ) ≡ η µν ∂ ∂x ν = ∂ ∂t , − ∂ ∂x with the Minkowski metric η ≡ diag(1, −1, · · · , −1); D is an N × N matrix with each entry being a function of the derivative operators; S is a variable to be solved and have N components in general. It is typically a linear perturbation to some physical quantities of interest. We assume that the background is uniform and the elements of D are polynomials of the derivative operators with constant coefficients. This is not a serious limitation, since it is a common practice to apply linear analysis locally under the assumption that the typical wavelength of perturbation is much shorter than the scale height in the background.
The asymptotic behavior of S is obtained by considering the retarded Green function G that satisfies D(i∂)G(x) = δ (d+1) (x)I N and the condition G(x) = 0 for t < 0. Here δ (d+1) (x) is the (d + 1)-dimensional delta function and I N is the N × N unit matrix. The Laplace-Fourier transform enables us to express it explicitly as M is the Laplace-Fourier contour (see Fig. 1(a)), the orientation of which is given by dk 0 r ∧ dk 1 r ∧ · · · ∧ dk d r with k µ r ≡ Re k µ ; the inner product of two (d + 1)-dimensional vectors a and b is denoted by a · b ≡ η µν a µ b ν . In the above expression the asymptotic behavior of G at t → 0 is considered in a frame moving at a velocity v so that we could investigate the convective instability. Note that the absolute instability corresponds simply to the case with v = 0. Incidentally v µ is the (d + 1)-dimensional velocity given by v µ = (1, v). The inverse matrix of D(k) is given as D −1 (k) = adj D(k)/∆(k) in terms of the adjugate matrix adj D(k) and the determinant ∆(k) of D(k) and hence has pole singularities at the zeros of ∆(k). Note that the real solutions of ∆(k) = 0 give the DR of stably propagating modes.
For t > 0, the integration contour M can be modified to M on the section at Im k 1 = · · · = Im k d = 0 by the Jordan's lemma (see Fig. 1(b)) so that all the poles should be enclosed. Then the integration of Eq. (1) can be reduced to the integral of the Poincaré residue of the integrand by the residue theorem [12,15]: where the integration contour C ≡ D ∩ C × R d is the section of D at Im k 1 = · · · = Im k d = 0, with D being a locus defined by D ≡ k ∈ C d+1 |∆(k) = 0 ; the orientation of C is given by dk 1 r ∧ · · · ∧ dk d r C ; here and henceforth ∂ µ are k-derivatives. We note that the integrand is actually its restriction on D but is abbreviated for notational simplicity. The integrand is holomorphic unless there exists some k ∈ D such that d∆ = 0 on that point. In the following we assume d∆ = 0 for all k ∈ D, since it is satisfied indeed for almost all ∆ and even if not, we could always modify ∆ slightly by, for instance, adding k 0 with a small ∈ R so that it should satisfy the condition.
To obtain the asymptotic behavior of Eq. (2) in the limit of t → ∞, we utilize the Lefschetz thimble method on D embedded into the (d + 1)-dimensional complex space C d+1 , g (see Fig. 1(c)), with a Kähler metric g satisfying g αβ = gᾱβ = 0, g αβ = gβ α = δ αβ /2. It expresses the integration contour C as the sum of Lefschetz thimbles {J σ }, which are nothing but the steepest descent paths: C ∼ = σ C, K σ J σ , where the coefficient ·, · is referred to as the intersection form on D. Then Eq. (2) can be rewritten as in which the sum runs over the critical points of a Morse function h on D defined as the height function: , which corresponds to the real part of the exponent of e −ik·vt . Note that the intersection form ensures that only those critical points of relevance are picked up. This is exactly what was lacking in the previous theories.
There are three steps in order to evaluate the asymptotic limit of Eq. equations: which are the stationary condition for h constrained on D. After some calculations, they are written as We note that the index of (v µ ) = (1, −v) is lowered by the Minkowski metric η.
In the second step, we construct the dual thimbles K σ by solving ordinary differential equations. Both the Lefschetz thimbles J σ and their dual thimbles K σ consist of all points on the curves K(s) that satisfies the upward flow equation for h constrained on D: where is the orthogonal projector onto the hypersurface given by f (k) = const.; the indices with tilde stand collectively for indices with and without bar; we use the following notation for contraction: aαbα ≡ g αβ a α b β + gᾱ β a α b β . Equation (6) can be simplified by straightforward calculations as where a ≡ δ αβ a α a β , and ensures that K(s) satisfies the following relations for arbitrary v and s: • the constraint on D: ∆(K(s)) = 0 (9) • the stationarity of the phase of e −ik·vt : • the monotonicity of the amplitude of e −ik·vt : The difference between J σ and K σ is the boundary condition they satisfy: K(∞) = k σ for J σ and K(−∞) = k σ for K σ . Note that the orientations of J σ and K σ are chosen so that the following orthogonal relations should hold: J σ , K σ = δ σσ . It is incidentally mentioned that since h depends on v, so do k σ , J σ and K σ . We do not show it explicitly, though, for notational simplicity unless it is important. One needs to solve (numerically in general) the flow equation (8) with the initial condition K(−∞) = k σ for each critical point to construct the dual thimbles attached to the critical point k σ and obtain the intersection form C, K σ . This is the cost we have to pay in this method to avoid the costly and almost impossible exploration of the complex DR over the entire complex k-space. As a technical tip, we suggest to set the initial condition for the integration as K(0) = k σ + κ, where a small shift κ is chosen as where J σ is given below and δ ∈ R 3 is a small real vector. This prescription ensures that K(0) sits on the dual thimble K σ . For each δ one obtains a single upward flow line and K σ is the collection of all these upward flow lines. We do not need all of them in fact. It should be noted that many flow lines seem to be attracted to the steepest flow lines as s evolves as shown in Fig. 2. Then the whole picture of a dual thimble may not be captured unless the upward flow lines are sufficiently many. This concentration of flow lines occurs at smaller scales for smaller perturbations of κ indeed. We should hence choose appropriately large (and not too large of course) perturbations; otherwise large numbers of flow lines are required to recover the entire shape of a dual thimble.
The intersection form C, K σ is easily evaluated as follows. We focus on the imaginary parts of the spatial coordinates of C and K σ , i.e., we project them onto the section at k 0 = Re k 1 = · · · = Re k d = 0. Then C is reduced to the origin Im k 1 = · · · = Im k d = 0 whereas the projection of K σ is still d-dimensional in general. The evaluation of C, K σ is hence reduced to the investigation of whether the projection of K σ includes the origin or not. In other words, C, K σ = ±1 if the section of K σ at s → ∞ (or sufficiently large s practically), which is isomorphic to S d−1 , encloses the origin in the projected space; C, K σ = 0 otherwise. Importantly from a practical point of view, this can be visualized for the spacetime dimension less than 5.
The last step is to evaluate the integral on each Lef-schetz thimble in the limit of t → ∞. Since the Lefschetz thimble is nothing but a steepest descent path from a critical point, the integral is dominated by the contribution from the vicinity of the critical point. Then the asymptotic form of G for |k σ · vt| |k σ · x| is given as In deriving this expression, we expand the exponent of e −ik·vt up to the quadratic order at the critical point and conduct the resultant Gaussian integrals, assuming that all other factors are subdominant and ignoring their variations; the matrix J σ = (ι 1 , · · · , ι d ) is defined so that it should satisfy H σ = − t J −1 σ J −1 σ for the following matrix (14) and the order of the column vectors {ι i } are chosen so that ι 1 ∧ · · · ∧ ι d should give the orientation of the Lefschetz thimble at k σ , where we introduce the (d + 1)dimensional vectors ι i ≡ (v · ι i , ι i ) associated with ι i so that they should be tangential to D. Note also that the Hessian of h at k σ on D is then given by Re (H σ ) ij dk i ⊗ dk j (see [16] for the derivation). The choice of J σ is not unique and the results are unchanged by the transformation J σ → J σ V for V ∈ SO(d, C). Moreover, since one can express det J σ as det J σ = e i arg det Jσ |det H σ | −1/2 and its argument is invariant under the transformation J σ → J σ M for M ∈ GL(d, R), the absolute value and argument of det J σ are obtained, respectively, from | det H σ | and conveniently chosen J σ that is still amenable to the orientation of J σ .
Finally, if det H σ = 0 at some k σ , it means that h is not a Morse function and the formulation given so far cannot be applied at this point. Even in this case, however, h can be modified to become a Morse function just by changing v infinitesimally, since Morse functions are densely existent. Hence our results are always valid essentially.

B. Maximum growth rate
The asymptotic limit of G(t, x + vt) for t → ∞, Eq. (13), is valid for arbitrary v. Then h(k σ , v) = Im(k σ · v) can be regarded as the growth rate of the contribution from a critical point k σ . The natural question next should be in which frame the growth rate of instability becomes maximum. The answer can be obtained simply as follows if the Im k 0 has a maximum in C.
We first define k m as where the last equality is satisfied because h(k, v) = Im(k · v) = Im(k 0 − k · v) and k ∈ R d for k ∈ C. It is actually one of the critical points for the frame moving at the velocity is since Eq. (5) is satisfied at k m indeed. In this frame, the dual thimble K m (v m ) associated with this critical point intersects with C at k = k m as the critical point is sitting on C. We hence have C, K m (v m ) = ±1. The corresponding growth rate is h(k m , v m ) = Im k 0 m by definition. Now we see that C has no intersection with any dual thimble K σ (v) associated with those critical points k σ (v) that satisfy This is because, if otherwise, the above inequality contradicts the relation h(k, v) ≥ h(k σ , v), which should hold for all points on K σ (v). Therefore, h(k σ , v) ≤ Im k 0 m is satisfied for all critical points that contribute to the sum in Eq. (13), which implies that the growth rate takes its maximum value Im k 0 m in the frame with the velocity v m . Although the maximum growth of instability is certainly one of the most relevant quantities, we stress that we are equally interested in the entire picture of the instabilities, convective and absolute alike, i.e., we want to know in which frame the instability grows at what rate. Such complete information can be also obtained in our formulation by evaluating Eq. (13) for all critical points.

A. Spatially 2-dimensional dispersion relation
Here we demonstrate how our new method works, employing the following DR given in [9]: This toy model has a merit that it can be treated completely analytically. In fact the Morse function h(k, v) is shown in Fig. 3 as a function of v at the critical point k + , which gives the highest value to h. It is found that which is hence a necessary condition for instability for this DR. To study whether this is also the sufficient condition or not in the classical approach, we have to accomplish the following tasks for each v in the range given above. We first solve the simultaneous equations {∆(k) = 0, ∂ 1 ∆(k) = 0} for k 0 = {k 0 σ + is|s ≥ 0} with an arbitrary positive s for every critical point so that the trajectories (k 1 (k 0 ), k 2 (k 0 )) could be obtained. We then investigate for all points k =k on the trajectories whether there occurs a coalescence between two of the roots k 1 (k 0 ,k 2 ) of ∆(k 0 , k 1 ,k 2 ) = 0 at k 0 =k 0 . If it is true, those points on the trajectories are called coalescence roots. Finally, we have to survey further whether two of such trajectories (k 1 (k 0 ), k 2 (k 0 )) of the coalescence roots meet at one of the critical points or not. Although all these tasks can be accomplished for the current toy model, in which the DR is fairly simple and the spatial dimension is just 2, it is not difficult to imagine how tough that would be for more complicated DR's and/or higher spatial dimensions.
In our new formulation, the procedure is simpler. We have only to solve Eq. (8) with the initial condition K(−∞) = k σ to construct the dual thimble K σ for each critical point to find the intersection form C, K σ by examining whether K σ includes the origin in the Im k 1 − Im k 2 plane.
We can then verify that the condition (v 1 − 1) 2 + (v 2 − 1) 2 < 1 is indeed the sufficient condition. We show in Figs. 4(a), (b), (c) and (d) the critical point k + as well as the corresponding K + obtained by solving Eq. (8) numerically. In Figs. 4(a) and (b), K + intersects with C only once, implying C, K + = ±1 and that this critical point contributes to the asymptotic behavior. On the other hand, in Figs. 4(c) and (d), the intersection of K + with C occurs twice with opposite orientations and, as a result, C, K + = 0 in these cases. We can also confirm that the dual thimble K − associated with the other critical point k − has C, K − = 0. It turns out that the Green function G(t, x+vt) is exactly 0 at x = 0 for these values of velocity v. Note that although it is plotted in Figs. 4(a)-(d), Im k 0 is not necessary in fact to obtain the intersection form C, K σ .

B. Spatially 3-dimensional dispersion relation
Next we demonstrate our new method for spatially 3dimensional problem with the DR which is a straightforward extension of Eq. (18) to 4dimensional spacetime. The necessary condition for instability, under which h(k + , v) takes a positive value, is In Figs. 5, 6, 7 and 8, the dual thimble K + attached to the highest critical point k + is shown for some different v. The dual thimbles are now 3-dimensional manifold. We hence show their sections at some different s, which are isomorphic to S 2 , for better visibility. Note that the projection of one of these sections onto the section at k 0 = Re k 1 = · · · = Re k d can have self-intersections although it does not in the (d + 1)-dimensional complex space (just like the Klein bottle); this situation is indeed realized in Figs. 7 and 8. In Figs. 5 and 6, we can verify that the origin is enclosed in the sections of K + at sufficiently large s and hence C, K + = ±1. On the other hand, Figs. 7 and 8 show that the origin is not enclosed in the sections, implying C, K + = 0. We can then confirm that (v 1 − 1) 2 + (v 2 − 1) 2 + (v 3 − 1) 2 < 1 is also the sufficient condition for instability for the DR given by Eq. (19).

IV. CONCLUSION
Based on the Lefschetz thimble method, we have reformulated linear instability analysis for field quantities that obey partial differential equations in an unlimited, not necessarily one-dimensional, spatial domain. Instead of detecting coalescence of two complex roots of the DR as in the conventional theory, which is all but impossible in more than one spatial dimension in fact, we employ the intersection form defined between the original integration contour and the dual thimbles, which are obtained by solving the gradient flow equations. We have derived the explicit asymptotic formula for the retarded Green function that gives the evolution of perturbations in a wave packet form. Our new formulation is not only more mathematically rigorous than the classical method, properly picking up the critical points and associated Lefschetz thimble, or the steepest descent paths, that contribute to the asymptotic behavior, but is also practically useful mainly because we do not have to solve simultaneously the DR ∆(k) = 0 and some of its derivatives ∂ i ∆(k) = 0 explicitly and repeatedly for essentially all k ∈ C d+1 as in the conventional theory. Since the basic equation we assumed are quite generic, we believe that this formulation will have broad applications. As a matter of fact, we are currently applying the method to the analysis of collective neutrino oscillations, a nonlinear problem, which is notorious for its difficulties in direct numerical solutions and is our original motivation for this work. The results will be reported elsewhere soon. It may be also interesting to apply the method with some extensions possibly to the problem with oscillating sources.