GRADIENT FLOW APPROACH TO THE CALCULATION OF STATIONARY STATES ON NONLINEAR QUANTUM GRAPHS

. We introduce and implement a method to compute stationary states of nonlinear Schr¨odinger equations on metric graphs. Stationary states are obtained as local minimizers of the nonlinear Schr¨odinger energy at ﬁxed mass. Our method is based on a normalized gradient ﬂow for the energy (i.e. a gradient ﬂow projected on a ﬁxed mass sphere) adapted to the context of nonlinear quantum graphs. We ﬁrst prove that, at the continuous level, the normalized gradient ﬂow is well-posed, mass-preserving, energy diminishing and converges (at least locally) towards stationary states. We then establish the link between the continuous ﬂow and its discretized version. We conclude by conducting a series of numerical experiments in model situations showing the good performance of the discrete ﬂow to compute stationary states. Further experiments as well as detailed explanation of our numerical algorithm are given in a companion paper.


Introduction
Partial differential equations on (metric) graphs have a relatively recent history.Recall that a metric graph G is a collection of vertices V and edges E with lengths l e ∈ (0, ∞] associated to each edge e ∈ E. One of the earliest account of a partial differential equation set up on metric graphs is the work of Lumer [38] in 1980 on ramification spaces.Among the early milestones in the development of the theory of partial differential equations on graphs, one finds the work of Nicaise [42] on propagation of nerves impulses.Since then, the theory has known considerable developments, due in particular to the natural appearance of graphs in the modeling of various physical situations.One may refer to the survey book [20] for a broad introduction to the study of partial differential equations on networks, with a special emphasis on control problems. Among partial differential equations problems set on metric graphs, one has become increasingly popular: quantum graphs.By quantum graphs, one usually refers to a metric graph G = (V, E) equipped with a differential operator H often referred to as the Hamiltonian.The most popular example of Hamiltonian is −∆ on the edges with Kirchhoff conditions (conservation of charge and current) at the vertices (see Section 2 for a precise definition), where ∆ is the Laplace operator.The book of Berkolaiko and Kuchment [16] provides an excellent introduction to the theory of quantum graphs.
Recently, another topic has gained an incredible momentum: nonlinear quantum graphs.By this terminology, we refer to a metric graph G = (V, E) equipped with a nonlinear evolution equation of Schrödinger type i∂ t u − Hu + g(|u| 2 )u = 0, where u = u(t, x) ∈ C is the unknown wave function, t denoting the time variable and x the position on the edges of G. Whereas the research on linear quantum graphs is mainly focused on the spectral properties of the Hamiltonian, one of the main area of investigation for nonlinear quantum graphs is the existence of ground states, i.e. minimizers of the Schrödinger energy E on fixed mass M , where Indeed, ground states are considered to be the building blocks of the dynamics for the nonlinear Schrödinger equation, and being able to obtain them by a minimization process guarantees in particular their (orbital) stability.
On the theoretical side, the literature concerning ground states on quantum graphs is already too vast to be shortly summarized.A perfect introduction to the topic is furnished by the survey paper of Noja [43] and we only present a few relevant samples.
Among the model cases for graphs, the simplest ones may be star-graphs, i.e. graphs with one vertex and a finite number of semi-infinite edges attached to the vertex (see Figure 1).For this type of graphs with an Star-graph with N = 6 edges attractive Dirac type interaction at the vertex, Adami, Cacciapuoti, Finco and Noja [4,5] established under a mass condition and for sub-critical nonlinearities the existence of a (local or global) minimizer of the energy at fixed mass, with an explicit formula for the minimizer (see Section 2 for more details and explanations).For more general nonlinear quantum graphs, Adami, Serra and Tilli [9,10,11] have focused on the case of Kirchhoff-Neumann boundary conditions for non-compact connected metric graphs with a finite number of edges and vertices.In particular, they obtained a topological condition (see Assumption 2.3 (H)) under which no ground state exists.On the other hand, in some cases, metric properties of the graph and the value of the mass constraint influence the existence or non-existence of the ground state [10,11,23,44].
Another particularly interesting study is presented in the work of Marzuola and Pelinovsky [40] for the dumbbell graph.As its name indicates, the dumbbell graph is made of two circles linked by a straight edge (see Figure 2).It is shown in [40] that for small fixed mass, the minimizer of the energy is a constant.As the mass increases, several bifurcations for the ground state occur, in particular a symmetric (main part located on the central edge) and an asymmetric one (main part located on one of the circles).Numerical experiments (based on Newton's iteration scheme) complement the theoretical study in [27,40].

• •
. Dumbbell graph Among the many other interesting recent results on nonlinear quantum graphs, we mention the flower graphs studied in [35], graphs with generals operators and nonlinearities [33], periodic graphs [45], etc.
On the numerical side, however, the literature devoted to nonlinear quantum graphs is very sparse.Finite differences on graphs have been implemented in a library developed in Matlab by R. H. Goodman, available in [28] and which has been used in particular in [27,36].The work [40] is one of the rare work containing numerical computation of nonlinear ground states on graphs.In our case, we have implemented a finite difference discretization scheme (see Section 4.2) in the framework of the Grafidi library [17], a Python library which we have developed for the numerical simulation on quantum graph and which is presented in the companion paper [18].
The integrability of the cubic nonlinear Schrödinger equation on graphs is analyzed in [50], with some numerical simulation and an appendix discussing the discretization at the vertices.The fully discrete (Ablowitz-Ladik type) integrable nonlinear Schrödinger is studied in [41].Other model equations on graphs are considered in [48,49].Extension to transparent vertices conditions is proposed in [52,53,54].
Our goal in this paper is to develop numerical tools for the calculation of local minimizers of the energy at fixed mass m > 0 in the setting of generic metric graphs with non necessarily Kirchhoff vertex boundary conditions.
The numerical method that we have implemented corresponds to a normalized gradient flow: at each step of time, we evolve in the direction of the gradient of the energy and renormalize the mass of the outcome.Such scheme is popular in the physics literature under the name "imaginary time method".One of the earliest mathematical analysis was performed by Bao and Du [15].More recently, in the specific case of the nonlinear Schrödinger equation on the line R with focusing cubic nonlinearity, Faou and Jezequel [24] performed a theoretical analysis of the various levels of discretization of the method, from the continuous one to the fully discrete scheme.
At the continuous level, by considering a function ψ(t, x) on G, the normalized gradient flow is given by and we establish in Section 3 the main properties of the flow.This is our first main result, which can be stated in the following informal way.
Main result 1.1 (see Theorem 3.2).Under Assumptions 2.1 and 3.1, the continuous normalized gradient flow is well-posed, mass preserving, energy diminishing, and converges locally towards local minimizers.
Having established the adequate properties of the flow at the continuous level, we turn to the discretization process.As is explained in Section 4, several time-discretizations are possible, but the so-called Gradient Flow with Discrete Normalization has proven to be very efficient.It consists into the following process to go from ψ n (an approximation of ψ(t n , •) at discrete time t n ) to ψ n+1 : (GFDN) The space discretization can be performed using second order finite differences inside the edges.The values at the vertices are obtained by approximating by finite differences the boundary conditions at the vertices.
In our second main result, we establish the link between the continuous normalized gradient flow and its space-time discretization.
Main result 1.2 (see Section 4).The Gradient Flow with Discrete Normalization (GFDN) is a time-discretization of the continuous normalized gradient flow (CNGF).Its space discretization can be obtained by finite differences with a special treatment at the vertices.
Finally, we illustrate by numerical experiments the efficiency of our technique.We use as test case the 2-star graph with δ and δ boundary conditions at the vertex connecting the two edges.This test case has been extensively studied from a theoretical point of view (see [25,26,37] for earlier works and [1] and the references therein for more recent achievements).A sneak peek of the results presented in Section 5 is offered in Figure 3 where the almost perfect agreement between the theoretical solution and the computed one is shown in the case of a 2-star graph with attractive δ condition at the vertex.We also consider other possible types of graphs.Further numerical experiments as well as a detailed presentation of our numerical algorithm are given in the companion paper [18].
Our main achievements in the numerical experiments are summarized in the following statement.
Main result 1.3 (see Section 5).The observed convergence of the discretized flow is of order 2 in space.In the test case of a nonlinear Schrödinger equation on a star graph with two edges and attractive δ or δ interactions at the vertex, the discretized flow converges towards the explicitly known ground state.Applicability of the method to generic graphs is illustrated on the sign-post graph and the tower of bubbles graph.
The rest of this paper is organized in the following way.In Section 2, we present in details the setting in which we work and give theoretical preliminaries.In Section 3, we prove that the continuous normalized gradient flow is well-posed, energy diminishing and converges locally towards a stationary state.In Section 4, we present the space-time discretization process of the continuous flow.Finally, numerical experiments in a test case and in more elaborate settings are presented in Section 5.

Preliminaries
We start with a few preliminaries to give the precise setting in which we would like to work.
2.1.Linear quantum graphs.Let G be a metric graph, i.e. a collection of edges E and vertices V. We assume that G connected.Two vertices might be connected by several edges and one edge can link a vertex to itself.Each of the edges e ∈ E will be identified with a segment where l e is the (finite or infinite) length of the edge.A (complex valued) function ψ : G → C is a collection of one dimensional maps defined for each edge e ∈ E: We define L p (G) and H k (G) by The corresponding norms will be given by The scalar product on L 2 (G) will be given by

Re
Ie φ e ψe dx.
To denote the duality product between H 1 (G) and its dual we will use the angle brackets: Note that it is common to include in the definition of H 1 (G) a continuity condition at the vertices.In order to consider more general situations, we do not make this restriction here and we will later instead introduce the space H 1 D (G), which corresponds to the Dirichlet part of the compatibility conditions at the vertices (see (3)).Given u ∈ H 2 (G) and a vertex v ∈ V of degree d v , define u(v) ∈ R dv as the column vector where e ∼ v denotes the edges incident to the vertex v and u e (v) is the corresponding limit value of u e .The boundary conditions at the vertex v will be described by For the sake of conciseness, we use the notation for the column vector of all values at the end of the edges and the corresponding boundary conditions matrices are given by The boundary conditions considered are local at the vertices, we refrain here from taking into account more general boundary conditions.We now define on the graph a second order unbounded operator H by where the domain of H is given by and the action of H on u ∈ D(H) is given by (Hu) e = −∂ xx u e for every edge e ∈ E. We restrict ourselves to self-adjoint operators, which is known to be equivalent for H (see e.g.[16,Theorem 1.4.4]) to request that at each vertex v the d v × 2d v matrix (A v |B v ) has maximal rank and the matrix A v B * v is symmetric.In that case, for each vertex v there exist three orthogonal and mutually orthogonal operators P D,v (Dirichlet part), P N,v (Neumann part) and P R,v = Id − P D,v − P N,v (Robin part), acting on C dv and an invertible self-adjoint operator Λ v acting on the subspace P R,v C dv such that the boundary values of u ∈ D(H) at the vertex v verify Using this expression of the boundary conditions, we can express (see e.g.[16,Theorem 1.4.11]) the quadratic form corresponding to H, which we denote by Q and is given by The domain of Q is given by all functions u ∈ H 1 (G) such that at each vertex P D,v u = 0. We denote it by We now consider two examples of boundary conditions: Kirchhoff-Neumann and δ-type.We already recalled what the classical Kirchhoff-Neumann boundary conditions (1) are.In terms of the projection operator, the Dirichlet part P D,v in the Kirchhoff-Neumann case is simply the projection on the kernel of B v , given by The Neumann part is given by I − P D,v , precisely and there is no Robin part.
We consider now a vertex with a δ-type condition of strength α v ∈ R at the vertex v, which is defined for u ∈ H 2 (G) as follows: u is continuous at v, This vertex condition is analogous to the jump condition appearing in the domain of the operator for the celebrated Schrödinger operator with Dirac potential (see e.g. the reference book [13] and Section 2.2.1).In terms of A v and B v matrices, the condition takes the form . Assuming that we have δ-type conditions on the whole graph, the domain H 1 D (G) of the quadratic form Q associated with H is the space of functions of H 1 (G) continuous at each vertex, and we thus may write u(v) for the unique scalar value of u ∈ H 1 D (G) at each vertex.The quadratic form associated with H then becomes 2.2.Nonlinear quantum graphs.Having established the necessary preliminaries on linear quantum graphs in the previous section, we now turn to nonlinear quantum graphs.Given a quantum graph (G, H), we consider the nonlinear Schrödinger equation on the graph G given by where u = u(t, •) ∈ L 2 (G) is the unknown wave function, t the time variable, and f is a nonlinearity satisfying the following requirements.
Assumption 2.1.The nonlinearity f : C → C verifies the following assumptions.
• Gauge invariance: there exists Typical examples for f are power type or double power type nonlinearities where 1 < p, q < ∞.We will use the real form of the anti-derivative of f , which is given for every z ∈ C by Using the antiderivative G of g, we may also express F (as we did in the Introduction) as Observe that f is a function defined on C. Its differential df at z ∈ C might be expressed for h ∈ C by The functions on which f will be evaluated in the next sections will mostly be real-valued and for simplicity we will use the following notation when the argument of f is real: for s ∈ R we define Formally, ( 4) is a Hamiltonian system in the form where the Hamiltonian, or the energy, E is a conserved quantity defined for any u ∈ H 1 D (G) by It is a C 2 functional on H 1 D (G) and its derivative is given by with the slight abuse of notation that H here denotes the corresponding operator from H 1 (G) to its dual.From Noether's theorem, the gauge symmetry of (4) yields another conserved quantity (see e.g.[21]), the mass, given by . We are interested in this paper in the standing waves solutions for the nonlinear Schrödinger equation set on the graph.By definition, a standing wave is a solution u of ( 4) given for all e ∈ E by where ω ∈ R and the profile φ ∈ H 1 (G) is independent of time.We refer to the profile φ as stationary state.Substituting into (4) leads to the equation of the profile φ, given by Therefore, φ is a critical point of the action functional Observe that there is a natural smoothing for φ: since, with our assumptions, (ωφ Strategies abound to find critical points of the action.One particularly interesting strategy is to minimize the energy on fixed mass, as the obtained minimizer will be (following the method established by Cazenave and Lions [19]) the profile of an orbitally stable standing wave of (4) (provided minimizing sequences are compact, which is usually a key step of the proof).More precisely, given m > 0, we will be looking for φ The theoretical existence of minimizers for the problem (7) has attracted a lot of attention in the past decade and we will not attempt to give an exhaustive overview of the existing literature.Some examples have already been shortly mentioned in Section 1.In what follows, we give a few more details on the case of star graphs with two or more edges, and on the topological assumption preventing the existence of ground states.

2.2.1.
Star graphs with two or more edges.One of the simplest nontrivial graph is given by two semi-infinite half-lines connected at a vertex, with δ type condition on the vertex.In this case, the operator H is equivalent to the second order derivative on R with point interaction at 0. In this setting, existence and stability of standing waves for a focusing power-type nonlinearity was treated by Fukuizumi and co.[25,26,37], using techniques based on Grillakis-Shatah-Strauss stability theory (see [30,31] for the original papers and [21,22] for recent developments).
Various generalizations have been obtained, e.g. for a generic point interaction [6,7,8] (δ or δ boundary conditions) or in the case of non-vanishing boundary conditions at infinity [34].In particular, the following results have been obtained in [8].
Proposition 2.2.Assume that G is formed by two semi-infinite edges {e 1 , e 2 } connected at the vertex v. Let H : D(H) ⊂ L 2 (G) → L 2 (G) be the operator −∂ xx with one of the following conditions to be satisfied at the vertex.
• Attractive δ conditions: • Dipole conditions: where 1 < p < 5. Then for any m > 0 there exists up to phase shift and translation a unique minimizer to A detailed review of these results as well as announcement of new results can be found in [1].The minimizer is in fact explicitly known, and we use its explicit form in Section 5.1 to compare the outcome of our numerical experiences with the theoretical ground states.

2.2.2.
General non-compact graphs with Kirchhoff condition.The existence of ground states with prescribed mass for the focusing nonlinear Schrödinger equation on non-compact graphs G equipped with Kirchhoff boundary conditions is linked to the topology of the graph.Actually, a topological hypothesis, usually referred to as Assumption (H) can prevent a graph from having ground states for every value of the mass (see [12] for a review).For the sake of clarity, we recall that a trail in a graph is a path made of adjacent edges, in which every edge is run through exactly once.In a trail, vertices can be run through more than once.The Assumption (H) has many formulations (see [12]) but we give here only the following one.

Assumption 2.3 (Assumption (H))
. Every x ∈ G lies in a trail that contains two half-lines.Under Assumption 2.3 (H), no global minimizer exists, unless G is (up to symmetries) isomorphic to R (note that this assumption does not prevent the existence of local minimizers).Let us consider for example a general N -edges star-graph G (see Figure 1).The N star-graph with N > 2 verifies Assumption 2.3 (H), so there are no ground states in this case without adding more constraints.Another example satisfying Assumption 2.3 (H) is the triple bridge B 3 (represented in Figure 4).When we are searching to obtain ground states, we consider graphs violating Assumption 2.3 (H), for example the signpost graph or a line with a tower of bubbles (Figure 5).
. Line with a signpost graph (left) and with a tower of bubbles (right)

Continuous normalized gradient flow
We want here to show that, when the standing wave profile φ is a strict local minimizer for the energy on fixed mass, the corresponding continuous normalized gradient flow (i.e. the gradient flow of the energy projected on the mass constraint) converges towards φ.
The continuous normalized gradient flow is defined by where ψ = ψ(t, •).It is the projection of the usual gradient flow be a standing wave profile solution of (6).We define the linearized action operator L + around φ by We will assume that the bound state φ is a strict local minimizer of the energy on fixed L 2 -norm, which translates for L + into the following assumption.
Assumption 3.1.There exists κ > 0 such that for any ϕ ∈ D(H) verifying Since the pioneering work of Weinstein [51], this assumption is well known to hold (if one removes translations and phase shifts) in the classical case of Schrödinger equations on R d with subcritical power-nonlinearities (f (ϕ) = |ϕ| p−1 ϕ, 1 < p < 1 + 4/d).It is has also been established in many different cases, for example in [34,37] in the case of the 2 branches star graph with δ conditions on the vertex (which is equivalent to the line with a Dirac potential) or in [32] in the case of a 1-loop graph with Kirchhoff conditions at the vertex (which is equivalent to an interval with periodic boundary conditions).Observe that a local minimizer is not necessarily a global minimizer (see e.g.[46]).
Our main result in this section is the following.
Theorem 3.2.Let the nonlinearity f and the bound state φ be such that Assumption 2.1 and Assumption 3.1 hold.Then for every 0 < µ < κ (where κ is the coercivity constant of Assumption 3.1) there exist ε > 0 and C > 0 such that for every 8) is global (i.e.T = ∞) and converges to φ exponentially fast: Remark 3.3.In Theorem 3.2, we may choose any µ > 0 such that µ < κ, where κ is the coercivity constant in Assumption 3.1.Hence the convergence rate towards the profile φ depends on the steepness of the energy well around φ.
The proof of the theorem is divided into three parts.This is the subject of the next three subsections.

3.1.
Local well-posedness of the continuous normalized gradient flow.Before proving Theorem 3.2, we establish the following local well-posedness result for the continuous normalized gradient flow (8).
Proposition 3.4.Assume that the nonlinearity f verifies Assumption 2.1.Then, for any ψ 0 ∈ H 1 D (G), there exists a unique maximal solution of the continuous normalized gradient flow (8) with T max ∈ (0, +∞].Moreover, the mass of the solution is preserved and its energy is diminishing, i.e. for all t ∈ (0, T max ) we have Proof of Proposition 3.4.Let ψ 0 ∈ H 1 D (G).We first show the second part of the statement: preservation of the mass.Let ψ be a solution of (8) as in the first part of the statement of Proposition 3.4.We have The mass is therefore preserved for (8).Set We now prove the first part of the statement (existence and uniqueness of a solution).We first consider the intermediate problem The intermediate problem (10) can be written more explicitly (using the expression ( 5) of E (ψ)) as where Q is the quadratic form associated with H and was defined in (2).Recall that the operator Moreover, there exists λ > 0 such that H −λ (this might be seen from the expression of Q given in (2) and the injection of Since f verifies Assumption 2.1, the nonlinearity f : is continuous, and, as a function f : , it is Lipschitz continuous on bounded sets.Indeed, for any Therefore, for any M > 0 and for any and a similar estimate holds for f .
The existence of the desired solution then follows from classical results in the theory of semilinear parabolic problems (see e.g.[39,47]).More precisely, there exists a unique (10) with ψ(0) = ψ 0 .
Given ψ, we now go back to the continuous normalized gradient flow (8) by proving that t → ψ L 2 is constant along the evolution in time.We have by direct calculations on (10) This is a first order linear ordinary differential equation in ψ 2 L 2 which may be solved explicitly: Since ψ 0 L 2 = α, this indeed gives ψ(t) 2 L 2 = α 2 for any t ∈ [0, T max ).Therefore ψ is also a solution of (8).Uniqueness of such a solution is a direct consequence of the uniqueness for (10) and the preservation of the mass.
Finally, we establish the energy diminishing property.Using (8) to replace E (ψ), we have where we have used the conservation of the mass in the form (ψ, ∂ t ψ) L 2 = 0 to obtain the last equality.This concludes the proof.
Having established local well posedness of the continuous normalized gradient flow (8), we turn our attention to the evolution for initial data in the vicinity of the bound state φ.The map χ is smooth and has bounded derivatives.Its inverse is the data-to-coordinates map χ −1 : As χ, the map χ −1 is smooth and has bounded derivatives.The second step of the proof of Theorem 3.2 is to decompose the continuous normalized gradient flow ( 8) by projecting it on W and φ, as is done in the following proposition.Proposition 3.5.Let T > 0 and ψ ∈ C((0, T ), D(H)) ∩ C 1 ((0, T ), L 2 (G)) be a solution of (8) such that ψ L 2 = φ L 2 and decompose ψ using the data-to-coordinates map χ −1 (ψ(t)) = (r(t), w(t)) ∈ R × D(H) given by (11): , where W is the dual of W and L + was defined in (9).
The proof of Proposition 3.5 is divided into three steps.In the first step we will consider the orthogonal decomposition of the flow along φ and W .In the second step we will project this orthogonal decomposition on the L 2 -sphere.The third and last step will make the link between the projected normalized energy derivative and the linearized action operator L + .

3.2.1.
Step 1: Orthogonal Decomposition.We first consider the orthogonal decomposition of the energy.
Consider the functional E RW : R × W → R defined for (r, w) ∈ R × W by Lemma 3.6 (Orthogonal decomposition of the energy).The functional E RW is differentiable and we have the following estimates For any element h ∈ W , we will note indifferently D w E RW (r, w)h or D w E RW (r, w), h the image of h by D w E RW (r, w).
Proof.Since E and χ are differentiable, the functional E RW is also differentiable and we have We now recall that φ ∈ H 1 D (G) satisfies (6).Given r ∈ R and w ∈ W , using (5), we have E (χ(r, w)) = Hχ(r, w) − f (χ(r, w)) We have used here the fact that f ∈ C 1 for the Taylor expansion, and that φ is bounded.We may now use this estimate in (12) to obtain ), which proves the first part of the statement.
From (13), for h ∈ W we get where to get the last line we have used that h ∈ W and thus φ, h = (φ, h) L 2 = 0.This proves the second part of the statement.

3.2.2.
Step 2: Projection on the L 2 -sphere.We now make the link between the orthogonal decomposition and the mass normalization constraint.We denote the L 2 sphere of radius φ L 2 by

Consider the open subset of W given by
We define the functional r W : O W → R for any w ∈ O W by the implicit relation The functional r W can be made explicit by a direct calculation on the above equality and is given for w ∈ O W by In particular, r W is well defined and smooth.Moreover, we have in the open set O W the estimate Thus, we have a local parametrization of S φ around φ given by

Introduce the functional E
This functional can be used to describe the dynamics of the projected part of the normalized flow, as is done in the following lemma.
Lemma 3.8 (Gradient flow in local variables).Let w be as in Proposition 3.5.Then w is a solution of Proof.Observe first that r W and E W are differentiable on O W . Their differentials are given, for w ∈ O W and h ∈ W such that w + h ∈ O W , by Using the derivatives of E RW (given in ( 12)-( 13)) and r, we might express D w E W (w) in the following way: Recall that ψ ∈ C([0, T ], ) is a solution of the continuous normalized gradient flow (8) such that ψ L 2 = φ L 2 and that ψ is decomposed using the data-to-coordinates map χ −1 (ψ(t)) = (r(t), w(t)) ∈ R × W given by (11) in the following way: Since r(t) = (ψ(t), φ) L 2 − 1 = r W (w(t)), the function r of t is C 1 .The regularity of w in t is the same as the regularity of ψ.We have , which, by conservation of the L 2 -norm for ψ implies that for all t ∈ [0, T ] we have We want to convert the continuous normalized gradient flow (8) in ψ into a closed equation for w (r can be directly deduced from w by preservation of the L 2 norm).Observe first that To obtain the evolution equation for w, we take h ∈ W and compute: Since ψ is a solution of the normalized gradient flow (8), we get We also have Using ψ, h = w, h , we get the following equation: Since the previous equation holds for any h ∈ W , it can be rewritten as By conservation of the L 2 -norm in the normalized gradient flow (8), r might be inferred from w and we have for w the following closed equation Using (16) to replace (D w E RW )(r W (w), w) in ( 17), we obtain where to get the last line we have used the expression of w 2 L 2 in terms of r W (w), i.e.
This concludes the proof.

3.2.3.
Step 3: Link with the linearized action.To conclude the proof of Proposition 3.5, it remains to make the link between D w E W (w) and L + .Lemma 3.9.The differential D w E W (w) can be approximated in the following way: where L + has been defined in (9).
Proof.We already have obtained in (16) the identity w.
From the estimates on (D w E RW )(r W (w), w) and (D r E RW )(r W (w), w) given in Lemma 3.6, we have where we have used the classical power series expansion ).This concludes the proof.

3.3.
Convergence of the normal part of the continuous normalized gradient flow.The second step of the proof of Theorem 3.2 is to prove convergence to 0 of the projected part w of the solution ψ of the continuous normalized gradient flow (8), provided w 0 H 1 is small (i.e.ψ 0 is close enough to the bound state φ).This is the object of the following proposition.
Proof of Proposition 3.10.Denote We remark that, for any t ∈ (0, T ), Thus, by denoting P W the orthogonal projector on W , since L + is self-adjoint and ∂ t w ∈ W verifies ( 15), we have where R is given by (18).By the coercivity estimate in Assumption 3.1 and Cauchy-Schwartz inequality, we have which implies in particular that Coming back to (19), we get Assume that w 0 H 1 < ε where ε > 0 is chosen such that (recall that 0 < µ < κ)).Since w is continuous, there exists T 0 > 0 such that for any t ∈ [0, T 0 ], we have For t ∈ [0, T 0 ] we integrate (20) in time from 0 to t and use (21) to obtain Defining the constant C 0 = L + w 0 , w 0 /κ and using again the coercivity estimate of Assumption 3.1, we get which by Gronwall inequality gives and therefore Note that there exists a constant C > 0 independent of w 0 such that thanks to (2) and Sobolev embeddings.This concludes the proof.

Space-time discretization of the normalized gradient flow
4.1.Time discretization.The discretization scheme of the continuous normalized gradient flow (8) must provide a numerical method to obtain a minimizer of (7).We first consider the semi-discretization in time.The time step δt > 0 is chosen to be fixed and the discrete times t n are defined as t n = nδt, n 0. The semi-discrete approximation of any unknown function ψ(•, x), x ∈ G, at time t n is denoted by ψ n (x).In order to present the numerical schemes, we recall that the nonlinearity verifies Assumption 2.1 and that it is of the form where g is continuous.We also introduce the variable µ m , usually referred to as the chemical potential, Dropping the dependence on (t, x) ∈ [0, +∞) × G and using ( 5), the continuous normalized gradient flow (8) can therefore be rewritten as where ψ 0 2 L 2 = m.Several numerical methods can be considered for discretizing (22).For example, if the nonlinearity is f (ψ) = |ψ| 2 ψ, a standard Crank-Nicolson scheme would consist in where the intermediate values at t n+ 1 2 are given by This method can be proved to be energy diminishing.However, in the above discretization, we need to solve a fully nonlinear system at every time step, which is time-and resource-consuming in practical computation.Bao and Du introduced in [15] a more efficient solution: the Gradient Flow with Discrete Normalization (GFDN) method, which consists into one step of classical gradient flow followed by a mass normalization step.By setting ψ 0 = ψ 0 , it is given by It is not clear at first sight that ( 23) is indeed a discretization of ( 22), but we have the following result.
Some arguments are given in [14,15] for m = 1, we provide here a proof with additional details and extend the result for any m > 0.
Proof of Proposition 4.1.The starting point is to apply a first order splitting, also known as Lie splitting, to (22).Assuming that the approximation ψ n of ψ at time t n , of mass ψ n 2 L 2 = m, is known, the steps of the splitting scheme are as follows.
Step 1 requires to solve a nonlinear parabolic type partial differential equation.Following [15], we approximate ( 24) by a semi-implicit time discretization: So, we have ϕ n+1 = v n+1 .The interest of having a semi-implicit scheme stems from its stability property.The equation involved in Step 2 is an ordinary differential equation.In [15], its solution is approximated by The coupling of ( 26) and ( 27) leads to the GFDN method.It is not totally obvious that ( 27) is actually an approximation of the solution w(t n+1 ) to (25).The normalization part ( 27) is actually equivalent to solving the ordinary differential equation where We define the piecewise function With this definition, the gradient flow with discrete normalization method ( 23) is an approximation of with ψ(t n ) 2 L 2 = m.Actually, the system (28) has to be read as the Lie splitting approximation of and ψ 0 2 L 2 = m.Thus, it remains to make the link between ( 29) and ( 22) by determining the limit of μm (t, δt) when δt goes to 0. Let us define t * = t n that remains constant when δt → 0 and n → ∞.For t * t < t * + δt, we have Consequently, for t * t < t * + δt, we have Thus, we conclude that ( 29) is an approximation of ( 22).This finishes the proof.
The complete Gradient Flow with Discrete Normalization algorithm is therefore where Id is the identity map and ε is a tolerance value.
Remark 4.2.Since we use a splitting scheme to discretize the continuous normalized gradient flow (8), it is no longer guaranteed that ( 30) is energy diminishing.We observe in numerical experiments that it is actually almost the case.
Remark 4.3.It is possible to modify the algorithm (30) to deal with mass one unknowns.Indeed, let us consider Then, the algorithm (30) becomes 4.2.Space discretization.We obtained the discretization in time of the normalized gradient flow in the previous section.To complete the discretization of the flow, we now proceed to the space discretization of the operator H.We recall that H is defined as a Laplace operator on each edge e ∈ E with boundary conditions given for each vertex v by ) r∼v are vectors, with ∂ x ψ e (v) the outgoing derivative on e at v, and (A v , B v ) are matrices (see Section 2).
For each edge e ∈ E, we consider N e ∈ N * the number of interior points and {x e,k } 1 k Ne a uniform discretization of the interval I e = [0, l e ], i.e.
x e,0 := 0 < x e,1 < . . .< x e,Ne < x e,Ne+1 := l e , with x e,k+1 − x e,k = l e /(N e + 1) := δx e for 0 k N e (see Figure 6).We denote v 1 the vertex at x e,0 , v 2 the one at x e,Ne+1 and, for any ψ ∈ H 1 D (G), for all e ∈ E and 1 k N e , ψ e,k := ψ e (x e,k ), as well as ψ e,v := ψ e (y v ) for v ∈ {v 1 , v 2 }, where y v1 = x e,0 and y v2 = x e,Ne+1 . × . Discretization mesh of an edge e ∈ E We now assume that N e 3.For any 2 k N e − 1, the second order approximation of the Laplace operator by finite differences on e is given by ∆ψ(x e,k ) ≈ ψ e,k−1 − 2ψ e,k + ψ e,k+1 δx e
For the case k = 1 and k = N e , the approximation requires ψ e,v1 and ψ e,v2 and we have to use the boundary conditions in order to evaluate them.We use second order finite differences to approximate them as well.For −2 j 0, we denote ψ e,v1,j = ψ e (x e,|j| ) and ψ e,v2,j = ψ e (x e,Ne+j ).
We have the approximation of the outgoing derivative from Assuming that δx = δx e for every edge e ∈ E to simplify the presentation, this leads to the approximation of the boundary conditions where ψ v,j = (ψ e,v,j ) e∼v .Assuming that 2δxA v + 3B v is invertible, this is equivalent to Thus, we can explicitly express the value of ψ e,v1 (resp.ψ e,v2 ) : it depends linearly on the vectors ψ v1,−1 and ψ v1,−2 (resp.ψ v2,−1 and ψ v2,−2 ).It is then possible to deduce an approximation of the Laplace operator at x e,1 and x e,Ne .That is, there exists (α e,v Since (ψ e,v,j ) −2 j 0,v∈{v1,v2} are interior mesh points from the other edges, we limit our discretization to the interior mesh points of the graph.The approximated values of ψ at each vertex are computed using (31).
We denote ψ = (ψ e,k ) 1 k Ne,e∈E the vector in R N T , with N T = e∈E N e , representing the values of ψ at each interior mesh point of each edge of G.We introduce the matrix [H] ∈ R N T ×N T corresponding to the discretization of H on the interior of each edge of the graph, which yields the approximation Hψ ≈ [H]ψ.4.3.Space-time discretization.Finally, we obtain the Backward Euler Finite Difference (BEFD) scheme approximating (23).Let ψ 0 ∈ R N T .We compute the sequence (ψ n ) n 0 ⊂ R N T given by where [g(|ψ n | 2 )] ∈ R N T ×N T is a diagonal matrix whose diagonal is the vector g(|ψ n | 2 ) and ϕ n+1 2 is the usual 2 -norm on the graph G of ϕ n+1 .This scheme has been studied on rectangular domains (with an additional potential operator) and Dirichlet boundary conditions [15] and is known to be unconditionally stable.Since it is implicit, the computation of ϕ n+1 involves the inversion of a linear system whose matrix is where [Id] is the identity matrix, and right-hand-side is ψ n .Using the matrix [H], we may also compute the energy.For instance, in the case where g(z) = z, by using the standard 2 inner product on the graph, we obtain  To illustrate our methodology, we give below an example of a star-graph with 3 edges (see Figure 7  Remark 4.4.We have implemented this space discretization in the framework of the Grafidi library [17], a Python library which we have developed for the numerical simulation on quantum graph and which is presented in [18].Note that finite differences on graphs have also been implemented in a library developed in Matlab by R. H. Goodman, available in [28] and which has been used in particular in [27,36].

Numerical experiments
We present here various numerical computations of ground states using the Backward Euler Finite Difference scheme (32).Even though the (BEFD) method was built for a general nonlinearity, for simplicity we focus in this section on the computations of the ground states of the focusing cubic nonlinear Schrödinger (NLS) equation on a graph G, that reads Explicit exact solutions are available for (NLS) on various graphs, in particular star graphs.We use the twoedges star graph in Section 5.1 to validate our implementation of the (BEFD) method and to show its efficiency to compute ground states.We present in Section 5.2 some numerical results for non compact graphs for which no explicit solutions are available.More examples are presented in a companion paper [18].
5.1.Two-edges star-graph.The two-edges star-graph is one of the simplest graph.We identify the graph G = G 2 as the collection of two-half lines connected to a central vertex A. Each edge is referred to with index i = 1, 2 (see Figure 8).The coordinate of vertex A is therefore both x 1 = 0 and x 2 = 0.The unknown ψ of (34) can be thought as the collection each function ψ i living on the edge i = 1, 2.
• A x 1 x 2  34) on the real line is known to be the soliton.To compute it on a two-edges star-graph, we identify the real line R to the graph G 2 with Kirchhoff condition at the vertex located at x = 0 (see [2]).The Kirchhoff condition on R is with ψ denoting the usual forward derivative, whereas on G 2 , it is The energy is The minimum of the functional E NLS among functions of H 1 (R) with squared L 2 -norm equal to m > 0 is given (up to phase and translation) by and In order to simulate the two semi-infinite edges originated from the central vertex, we consider two finite edges of length 40.The graph is presented on Figure 9 (left).We discretize each edge with N e = 4000 nodes and set homogeneous Dirichlet boundary conditions at the external vertices (B and C on Figure 9).The time step is δt = 10 −2 .The mass is m = 2.The exact solution is plotted on Figure 9 (right).The initial datum is chosen as a Gaussian of mass m/2 on each edge, namely We plot on Figure 10 both the exact solution φ m (35) and the numerical one φ m,num obtained after 3000 iterations (left), as well as the error |φ m − φ m,num | (right).The error is plotted for a fixed δx.We discuss the variation of the error with respect to δx in Figure 14 and observe that the scheme is of order 2. We obtain a very close numerical solution.Since the initial data is symmetric and centered on 0, our solution is also symmetric and centered on 0. The (BEFD) method allows to compute the exact energy and show that the scheme is energy diminishing.We plot in Figure 11 the evolution of the numerical energy using (33) and the comparison with the exact energy.The scheme is clearly energy diminishing and we obtain a very good agreement with the exact energy.5.1.2.δ-condition.We consider now a δ-condition at the central vertex A of the G 2 .The unknown ψ δ = (ψ δ,1 , ψ δ,2 ) T is the collection of ψ δ,i living on each edge i = 1, 2. Recall that the boundary conditions at The parameter α is interpreted as the strength of the δ potential and we focus on the attractive case (α < 0).The mass and energy are Explicit ground state solutions were provided in [29] in the cubic case (see also [3,8] for the general case).Define a by and define the function φ δ = (φ δ,1 , φ δ,2 ) by The mass of φ δ is explicitly given by m and the function φ δ has been constructed so that it is the minimizer of E δ (ψ δ ) with constrained mass m δ .The energy might be explicitly calculated : Like in the previous section, we apply the (BEFD) method to compute the ground state.We take the same numerical parameters concerning the mesh size and the approximation graph of Figure 9  in logarithmic scale (see Figure 12, right) confirms the accuracy of the numerical solution.Finally, we plot the evolution of the energy on Figure 13.We restrict ourselves to 1000 iterations on the horizontal axis since the convergence is really fast.The agreement with the exact energy is notable.Using the exact solutions when considering Kirchhoff and δ conditions, we are able to evaluate the order of the numerical scheme with respect to the spatial mesh size.As it was described in Section 4.2, the scheme should be of second order in space.To confirm this, we make various simulations for different mesh sizes δx and present the results in Figure 14 for both Kirchhoff and δ conditions.In the two cases, the order of convergence is 2, as expected.5.1.3.δ -condition.The δ -condition on star graph is usually defined by interchanging functions and their derivatives in the definition of the δ-condition (see e.g.[16]).We prefer here to use the concept of δ on the graph corresponding to the δ interaction on the line, and give a precise definition in what follows.As the previous section, the unknown ψ δ = (ψ δ ,1 , ψ δ ,2 ) is the collection of ψ δ ,i living on each edge i = 1, 2. The boundary conditions at with β > 0. The mass and energy are Explicit ground state solutions are provided in [8].Let us consider the transcendental system We are looking for real solutions such that x − < 0 < x + .
The ground state in both cases is given (up to a phase factor) by When 4/β 2 < ω 8/β 2 , the ground state of (34) with boundary condition (36) minimizing the energy E δ with fixed mass M (φ δ ) is an odd function.When 8/β 2 < ω < +∞, the ground state is asymmetric.The parameters for the numerical simulations are β = 1, δt = 10 −2 and we keep 4000 nodes per edges to discretize the two-edges graph (see Figure 9, left).The initial data are Gaussian on both edges but contrary to δ-condition, we select a different sign for the two edges (to increase the convergence speed).Namely, ψ 0 (x 1 ) = −ρe −10x 2 1 and ψ 0 (x 2 ) = ρe −10x 2 2 , with ρ > 0 such that M δ (ψ 0 ) = M δ (φ δ ).In order to simulate both odd and asymmetric ground states, we select respectively ω = 6 and ω = 16.The exact and numerical solutions are plotted in Figure 15.When looking to the comparison in logarithmic scale (see Figure 16), we see that we obtain a very 5.2.General non-compact graphs with Kirchhoff condition.We consider the computation of ground states on non-compact graphs not satisfying Assumption 2.3 (H).We focus on the signpost and tower of bubbles graphs.Beside the fact that they exist, little is known about minimizers.Our numerical algorithm is an easy to use tool to provide conjectures on the qualitative behavior of ground states on metric graphs.
5.2.1.Signpost graph.We now consider a signpost graph (see Figure 5).We wish to compute a stationary state of the NLS equation (34).The graph has the following dimensions: the line segment is equal to 2 and the perimeter of the loop is 4. The initial mass is taken to be 1.Concerning the length of the main line (supposedly very large), it is set to 100.On the discretization side, we set the total number of grid points to 5000.Furthermore, the time step is fixed to 10 −2 with a total number of iteration equal to 5000.The resulting stationary state is shown in Figures 18 and 19.We can see that it is localized in the loop and the line segment and decreases slowly along the main line.This is consistent with [9,10].Tower of bubbles graph.Finally, we have computed a stationary state for the tower of bubbles graph, with 2 bubbles.The graph is characterized by the following dimensions: the top bubble has a perimeter of 8 and the one of the bottom loop is set to 4. The main line, which is suppose to be very large, has a length of 100.The initial mass is taken to be 1.The space discretization is set by fixing a total number of grid points to 10000 and, for the time discretization, we have the time step set to 10 −2 for a total number of iterations of 10000.
We obtain the stationary state depicted in Figure 20 and 21.It is clear that, as for the signpost graph, the ground state is localized in the two bubbles and decreases slowly along the main line.Again, this is consistent with [9,10].

Figure 3 .
Figure 3.Comparison of numerical solution to ground state for δ interaction.
where A v and B v are d v × d v matrices and u (v) is formed with the derivatives along the edges in the outgoing directions.Consider for example the classical Kirchhoff-Neumann boundary conditions at the vertex v: we require the conservation of charge, i.e. for all e and e incident to the same vertex v u e (v) = u e (v), and the conservation of current, i.e. e∼v u e (v) = 0.These conditions are expressed in terms of A v and B v by we recover the classical Kirchoff-Neumann boundary conditions.When α v = 0, the Dirichlet, Neumann and Robin projectors are given as follows.The Dirichlet projector P D,v is (as when α v = 0) the projection on the kernel of B v .There is no Neumann part and the Robin part is given by I − P D,v (which was the Neumann part for α = 0).The operator Λ v = B −1 v A v on the range of P R,v is the multiplication by αv dv
Matrix representation of H.

Figure 7 .
Figure 7.An example for a star-graph.
(a)).The operator H is given with Dirichlet boundary conditions for the exterior vertices and Kirchoff-Neumann conditions for the central vertex.We can see on Figure 7 (b) the positions of the non zero coefficients of the corresponding matrix [H] when the discretization is such that N e = 10, for each e ∈ E. The coefficients accounting for the Kirchhoff boundary condition are the ones not belonging to the tridiagonal component of the matrix.

Figure 13 .
Figure 13.Evolution of the energy when computing the ground state of(34) with δ condition compared to E δ .

Figure 18 .Figure 19 .
Figure 18.The numerical ground state on the signpost graph

Figure 20 .Figure 21 .
Figure 20.The numerical ground state on the tower of bubbles graph