R Émi Carles

AND Christian Klein

Christof Sparber

ON GROUND STATE (IN-)STABILITY IN MULTI-DIMENSIONAL CUBIC-QUINTIC SCHR ÖDINGER EQUATIONS

Keywords: September 9, 2021. 1991 Mathematics Subject Classification. Primary: 35Q55. Secondary: 35C08, 65M70 Nonlinear Schrödinger equation, solitary waves, orbital stability, timesplitting method

We consider the nonlinear Schrödinger equation with a focusing cubic term and a defocusing quintic nonlinearity in dimensions two and three. The main interest of this article is the problem of orbital (in-)stability of ground state solitary waves. We recall the notions of energy minimizing versus action minimizing ground states and prove that, in general, the two must be considered as nonequivalent. We numerically investigate the orbital stability of least action ground states in the radially symmetric case, confirming existing conjectures or leading to new ones.

u |t=0 = u 0 ∈ H 1 (R d ).
The quintic modification of the cubic Schrödinger equation is a model which was introduced in the one-dimensional case in [START_REF] Pushkarov | Self-action of light beams in nonlinear media: soliton solutions[END_REF], as an approximate model in the framework of nonlinear optics. Equation (1.1) appeared more recently in the context of Bose-Einstein condensation, with d = 2 or 3: see e.g. [START_REF] Abdullaev | Stability of trapped Bose-Einstein condensates[END_REF][START_REF] Gammal | Atomic Bose-Einstein condensation with three-body intercations and collective excitations[END_REF][START_REF] Michinel | Liquid light condensates[END_REF], and [START_REF] Malomed | Vortex solitons: Old results and new perspectives[END_REF] for a review. In space dimensions d = 2 or 3, the impact of the quintic term on the dynamical properties of the solution u is stronger than in d = 1, as we shall discuss below.

Depending on the space dimension, which we always assume at most three to simplify the discussion, the nonlinearity in this model is seen to be: focusing L 2subcritical plus defocusing L 2 -critical (d = 1), focusing L 2 -critical plus defocusing H 1 -subcritical (d = 2), or focusing L 2 -supercritical plus defocusing H 1 -critical (d = 3). Recall that for the purely focusing cubic NLS, solitons exist in every dimension and finite time blow-up is possible provided d 2 (see e.g. [START_REF] Cazenave | Semilinear Schrödinger equations[END_REF]). The presence of the quintic nonlinearity prevents finite time blow-up in d = 2 or 3 (see Proposition 1.1 below), and also affects the stability of solitary waves: understanding this latter aspect more precisely is the main motivation for this paper. For a more precise discussion on the role of criticality in combined power nonlinearities see [START_REF] Tao | The nonlinear Schrödinger equation with combined power-type nonlinearities[END_REF].

The NLS (1.1) formally enjoys the following basic conservation laws:

(1) Mass:

M (u) = u(t, •) 2 L 2 (R d ) , (2)
Momentum:

P (u) = Im R d ū(t, x)∇u(t, x)dx,
(3) Energy:

E(u) = 1 2 ∇u(t, •) 2 L 2 (R d ) - 1 2 u(t, •) 4 L 4 (R d ) + 1 3 u(t, •) 6 L 6 (R d ) .
As evoked above, one important effect of the defocusing, quintic term is to prevent finite time blow-up which may occur in the purely cubic case. Indeed, the conservation of the energy, combined with Hölder's inequality, (1.2)

u 4 L 4 (R d ) u L 2 (R d ) u 3 L 6 (R d ) ,
shows that the focusing, cubic part cannot be an obstruction to the existence of a global in-time solution. More precisely, we have, in view of [START_REF] Cazenave | Semilinear Schrödinger equations[END_REF] for d = 2 and [START_REF] Zhang | On the Cauchy problem of 3-D energy-critical Schrödinger equations with subcritical perturbations[END_REF] for d = 3: Proposition 1.1 (Global well-posedness). Let d = 2, 3. For any initial data u 0 ∈ H 1 (R d ), the equation (1.1) has a unique solution u ∈ C(R; H 1 (R d )), such that u |t=0 = u 0 . This solution obeys the conservation of mass, energy, and momentum.

We note that in [START_REF] Lemesurier | Focusing and multifocusing solutions of the nonlinear Schrödinger equation[END_REF], numerical simulations are presented, in which the influence of a small defocusing quintic term on the time-evolution of a focusing cubic NLS is studied. In d = 2 and 3, and for initial data consisting of Gaussians, one obtains a time-periodic (multi-focusing) solution, similar to the one depicted in Fig. 11.

1.2. Orbital stability of action minimizing ground states. A particular class of global solutions are time-periodic solitary waves of the form u(t, x) = e iωt φ(x), with ω ∈ R and φ satisfying

(1.3) - 1 2 ∆φ + ωφ -|φ| 2 φ + |φ| 4 φ = 0, φ ∈ H 1 (R d ) \ {0}.
For d 3, solitary waves φ exist provided that the frequency ω satisfies the (necessary and sufficient) condition: 0 < ω < 3 16 , see [START_REF] Carles | Orbital stability vs. scattering in the cubic-quintic Schrödinger equation[END_REF]. Given an admissible ω ∈ (0, 3 16 ), we may then look for least action ground states, i.e. solutions φ ω (x) which minimize the action S ω (φ) = E(φ) + ωM (φ) among all nontrivial stationary solutions φ ∈ H 1 (R d ). Indeed, it is known from [START_REF] Byeon | Symmetry and monotonicity of least energy solutions[END_REF][START_REF] Cingolani | Multi-peak solutions for magnetic NLS equations without non-degeneracy conditions[END_REF] that every minimizer of the action S ω is of the form

(1.4) φ ω (x) = e iθ Q ω (x -x 0 ),
for some constant θ ∈ R, x 0 ∈ R d , and with Q ω the unique positive, radial solution to (1.3). In the following, we are mainly interested in the orbital stability of these specific solutions. To this end, we note that, as in the case of more standard, homogeneous nonlinearities (e.g. the cubic case), the NLS (1.1) enjoys three important invariances: (i) Spatial translation: if u(t, x) solves (1.1), then so does u(t, x -x 0 ), for any given x 0 ∈ R d .

(ii) Gauge: if u(t, x) solves (1.1), then so does e iθ u(t, x), for any given constant θ ∈ R. (iii) Galilean: if u(t, x) solves (1.1), then so does u(t, x -vt)e iv•x-i|v| 2 t/2 , for any given v ∈ R d . The first two invariants are seen to be present in formula (1.4). In combination with the third one, these invariants motivate the following standard notion of stability (see e.g. [START_REF] Cazenave | Semilinear Schrödinger equations[END_REF]): Definition 1.2. Let φ be a solution of (1.3). The standing wave e iωt φ(x) is orbitally stable in H 1 (R d ), if for all ε > 0, there exists δ > 0 such that if

u 0 ∈ H 1 (R d ) satisfies u 0 -φ H 1 (R d ) δ, then the solution to (1.1) with u |t=0 = u 0 satisfies sup t∈R inf θ∈R y∈R d u(t, •) -e iθ φ(• -y) H 1 (R d ) ε.
Otherwise, the standing wave is said to be unstable.

Following the breakthrough due to M. Weinstein, Grillakis, Shatah and Strauss introduced a general stability/instability criterion in [START_REF] Grillakis | Stability theory of solitary waves in the presence of symmetry. I[END_REF] (see also [START_REF] De Bièvre | Orbital stability: analysis meets geometry[END_REF]). Assuming certain spectral properties of the linearization of (1.3) about Q ω (which are satisfied in the present cubic-quintic case, see e.g. [START_REF] Carles | Orbital stability vs. scattering in the cubic-quintic Schrödinger equation[END_REF][START_REF] Killip | Solitons and scattering for the cubic-quintic nonlinear Schrödinger equation on R 3[END_REF][START_REF] Lewin | Uniqueness and non-degeneracy for a nuclear nonlinear Schrödinger equation[END_REF]), one has the following dichotomy:

(i) If ∂ ∂ω M (Q ω ) > 0, then e iωt Q ω (x) is orbitally stable, (ii) If ∂ ∂ω M (Q ω ) < 0, then e iωt Q ω (x)
is unstable. This criterion has proven extremely useful in the case of homogeneous nonlinearities, as well as in the case of mixed nonlinearities in d = 1 thanks to an explicit formula, cf. [START_REF] Iliev | Stability and instability of solitary waves for one-dimensional singular Schrödinger equations[END_REF][START_REF] Ohta | Stability and instability of standing waves for one-dimensional nonlinear Schrödinger equations with double power nonlinearity[END_REF]. In particular, when d = 1, all ground states solitary waves for (1.1) are orbitally stable. However, in the case d = 2 or 3, only partial results are currently available by using the criterion above, see Sections 2.2.2 and 2.2.3, respectively. In our numerical simulations, we will only consider radial perturbation of ground states, and thus remain in the radial framework. The notion of orbital stability then coincides with asymptotic stability up to a phase.

Remark 1.3. The dependence of δ upon ε in Definition 1.2 is unknown, in general. The proof of stability via the above criterion provides a rather explicit dependence of (a possible) δ as a function of ε, when ∂ ∂ω M (Q ω ) > 0 is known. On the other hand, the stability of the set of constrained energy minimizers (cf. Definition 2.1 below) is obtained by a non-constructive argument. In numerical simulations, tuning initial perturbations of a solitary wave which are sufficiently large to be visible, but not too large (to still adhere to the notion of stability), requires a subtle balance. 1.3. Constrained energy minimizers. Since solitary waves may be obtained by other means than minimizing the action S ω , one may want to look for alternative approaches to orbital stability. An important such alternative is obtained if for ρ > 0, we denote Γ(ρ) = u ∈ H 1 (R d ), M (u) = ρ , and assume that the constrained minimization problem

(1.5) u ∈ Γ(ρ), E(u) = inf{E(v) ; v ∈ Γ(ρ)} =: E min (ρ)
has a solution. We call such minimizers energy ground states, in order to make the distinction with least action ground states as clear as possible. Denote by E(ρ) the set of such solutions, i.e.

E(ρ)

:= {u ∈ H 1 (R d ), M (u) = ρ, E(u) = E min (ρ)}.
Now, let φ ∈ E(ρ): Then there exists a Lagrange multiplier Λ such that

dE(φ) = ΛdM (φ),
and thus, φ solves the stationary Schrödinger equation (1.3) for some (unknown) ω ∈ (0, 3 16 ). Observe that if φ ∈ E(ρ), then

{e iθ φ(• -y); θ ∈ R, y ∈ R d } ⊂ E(ρ).
When the nonlinearity is homogeneous (and L 2 -subcritical), this inclusion becomes an equality, see [START_REF] Cazenave | Orbital stability of standing waves for some nonlinear Schrödinger equations[END_REF][START_REF] Cazenave | Semilinear Schrödinger equations[END_REF]. However, for non-homogeneous nonlinearities like in (1.1), relating these two constructions of solitary waves (i.e., action minimizing ground states versus constrained energy minimizers) is not obvious at all, and the issue is possibly more complex than it may appear at a glance. First, a-priori nothing guarantees that an element of E(ρ) minimizes the action. Second, and this is more subtle: consider a least action ground state Q ω , and let ρ = M (Q ω ). It is not obvious, and not necessarily true, that Q ω ∈ E(ρ). In particular, the map ρ → ω may not be one-to-one. We will prove that, unlike in the case of homogeneous nonlinearities, we may indeed have Q ω ∈ E(ρ), cf. Theorem 2.5 below. This fact should be compared to the recent results of [START_REF] Jeanjean | On global minimizers for a mass constrained problem[END_REF]: In there, the authors establish for a large class of nonlinearities (including the cubic-quintic one in space dimension d 3), that all all energy minimizing ground states are least action ground states.

In addition, they show that if ω is obtained as the Lagrange multiplier associated to the mass constrained M (u) = ρ, then any least action solution of (1.3) at this value of ω is a constrained energy minimizer with the same mass ρ.

This paper is now organized as follows: In Section 2, we shall review several known mathematical results on the (in-)stability of ground state solitary waves in d = 1, 2, 3. In particular, we recall the fact that the set of energy minimizers is orbitally stable. We then prove that the dynamics of this set can be distinguished from dispersive behavior and that in d = 3 the sets of action and energy minimizers are not equivalent. In Section 3 we numerically construct action ground states and also collect several of their qualitative properties. Numerical evidence for the orbital stability of these action ground states in d = 2 is then given in Section 4, where we will also describe the numerical algorithm used to simulate the time evolution of (1.1). Finally, we shall turn to the question of orbital stability and instability of 3D action ground states in Section 5, where we will provide numerical evidence for several conjectures on the particular nature of the instability.

Mathematical results on orbital (in-)stability

2.1. Stability for energy ground states. As a preliminary step, we shall recall the Pohozaev identities for the cubic-quintic case (for a derivation, see e.g. [START_REF] Carles | Orbital stability vs. scattering in the cubic-quintic Schrödinger equation[END_REF])

: if φ ∈ H 1 (R d ) solves the stationary Schrödinger equation (1.3), then (2.1) 1 2 R d |∇φ| 2 dx - R d |φ| 4 dx + R d |φ| 6 dx + ω R d |φ| 2 dx = 0,
as well as

(2.2) d -2 2 R d |∇φ| 2 dx - d 2 R d |φ| 4 dx + d 3 R d |φ| 6 dx + ωd R d |φ| 2 dx = 0.
The aforementioned admissible range for ω ∈ (0, 3 16 ) is one of the consequences of these identities. In addition, if d = 2, and after multiplying (2.1) by 2 and subtracting (2.2), we find

0 = ∇φ 2 L 2 (R 2 ) -φ 4 L 4 (R 2 ) + 4 3 φ 6 L ( R 2 ) = 2E(φ) + 5 6 φ 6 L ( R 2 ) .
Therefore, any solitary wave in 2D has negative energy. In the 3D case, this is not necessarily so, as we will see below.

Next, we recall the notion of orbital stability for the set of energy minimizers, as introduced in [9]: Definition 2.1. We say that solitary waves are E(ρ)-orbitally stable, if for all ε > 0, there exists δ > 0 such that if

u 0 ∈ H 1 (R d ) satisfies inf φ∈E(ρ) u 0 -φ H 1 (R d ) δ, then the solution to (1.1) with u |t=0 = u 0 satisfies sup t∈R inf φ∈E(ρ) u(t, •) -φ H 1 (R d ) ε.
This notion is weaker than the one given in Definition 1.2, in the sense that E(ρ) may be a large set. In particular, we cannot infer orbital stability of individual members φ ∈ E(ρ), as required in Definition 1.2. In the case of homogeneous nonlinearities, however, one can prove that the set E(ρ) consists of only a single element φ (up to translation and phase conjugation), and thus one recovers Definition 1.2 from the one above. To prove E(ρ)-orbital stability via the concentrationcompactness method, the main step consists in showing that the minimal energy E min (ρ) < 0. In our case this yields: Theorem 2.2 (From [START_REF] Carles | Orbital stability vs. scattering in the cubic-quintic Schrödinger equation[END_REF]). Let d = 2 or 3.

(1) If E min (ρ) < 0, then E(ρ) is not empty, and the set of energy ground states is E(ρ)-orbitally stable.

(2) There exists ρ 0 (d) > 0 such that for ρ > ρ 0 (d), E min (ρ) < 0.

In particular, it seems reasonable to expect that no dispersion is possible near elements of E(ρ), in the sense that a solution u(t, •) within Definition 2.1 cannot satisfy

(2.3) u(t) L ∞ (R d ) -→ t→∞ 0.
We will use this criterion as a guiding principle for interpreting several of our numerical findings below. However, it is not clear a priori that for ρ > 0, such that

E(ρ) = ∅, we have inf{ φ L ∞ (R d ) , φ ∈ E(ρ)} > 0.
The property E min (ρ) < 0 makes it possible to rule out this scenario.

Proposition 2.3 (Non-dispersion of energy ground states). Let d = 2 or 3, and

ρ > 0. If E min (ρ) < 0, then m p := inf{ φ L p (R d ) , φ ∈ E(ρ)} > 0,
for any p ∈ (4, ∞]. In particular, there exists ε 0 (d) such that for 0 < ε ε 0 (d), any solution u provided by Definition 2.1 satisfies

inf t∈R u(t, •) L ∞ (R d ) > 0.
The statement of this proposition does not involve the above parameter ρ 0 (d).

For practical application, we want to emphasize that if for a given mass ρ, a stationary solution has negative energy, then solutions around E(ρ) cannot disperse.

Proof. Assume, by contradiction, that there exists a sequence (φ n

) n∈N ⊂ E(ρ) such that lim n→∞ φ n L p (R d ) = 0, for some p > 4. Since φ n L 2 (R d ) = ρ, by interpolation, this implies that lim n→∞ φ n L 4 (R d ) = 0.
In turn, this yields that

E min (ρ) = lim n→∞ E(φ n ) = lim n→∞ 1 2 ∇φ n 2 L 2 (R d ) + 1 3 φ n 6 L 6 (R d ) 0, a contradiction. Now, choose 0 < ε < m 5 . Consider initial data u 0 ∈ H 1 (R d ) such that inf φ∈E(ρ) u 0 -φ H 1 (R d ) δ,
where δ stems from Definition 2.1 and orbital stability (which is ensured since E min (ρ) < 0). Then, by the E(ρ)-orbital stability and Sobolev imbedding, we have

sup t∈R inf φ∈E(ρ) u(t, •) -φ L 5 (R d ) ε.
Since for all t and all φ ∈ E(ρ),

u(t, •) L 5 (R d ) φ L 5 (R d ) -u(t, •) -φ L 5 (R d ) m 5 -u(t, •) -φ L 5 (R d ) , this implies that inf t∈R u(t, •) L 5 (R d ) m 5 -ε > 0.
In particular, since u(t,

•) L 2 = const., an interpolation between L 2 and L ∞ proves inf t∈R u(t, •) L ∞ (R d ) > 0,
and hence (2.3) cannot hold.

Further mathematical results.

In the following we review some of the known results on orbital (in-)stability of ground states in d = 1, 2, 3. Moreover, we shall prove that 3D action ground states are not necessarily energy minimizers.

2.2.1.

The purely cubic case. For the cubic Schrödinger equation

(2.4) i∂ t u + 1 2 ∆u = -|u| 2 u, (t, x) ∈ R × R d ,
with d 3, the Cauchy problem is globally well-posed for d = 1 (both in L 2 (R) and H 1 (R), since the nonlinearity is L 2 -subcritical), while finite time blow-up is possible if d = 2 or 3, see e.g. [START_REF] Cazenave | Semilinear Schrödinger equations[END_REF]. Regarding the solitary waves, the analogue of (1.3) is

(2.5) - 1 2 ∆φ + ωφ -|φ| 2 φ = 0, φ ∈ H 1 (R d ) \ {0}.
The corresponding Pohozaev identities imply that such a non-trivial solution exists only if ω > 0. Conversely, for any given ω > 0, (2.5) has a unique positive radial solution. As a matter of fact, since the nonlinearity is homogeneous, the role of ω is explicit: consider Q cubic the unique positive radial solution in the case ω = 1,

(2.6) - 1 2 ∆Q cubic + Q cubic -Q 3 cubic = 0, x ∈ R d .
Then for any ω > 0,

φ ω (x) := √ ω Q cubic (x √ ω)
is a positive radial solution to (2.5). By uniqueness of such solutions, φ ω is an action ground state, minimizing

E cubic (φ) + ωM (φ) = 1 2 ∇φ 2 L 2 (R d ) - 1 2 φ 4 L 4 (R d ) + ω φ 2 L 2 (R d ) .
In particular, we readily compute

φ ω L 2 (R d ) = ω 1/2-d/4 Q cubic L 2 (R d ) .
Recalling the Grillakis-Shatah-Strauss stability criterion, this directly implies that cubic ground states are orbitally stable in d = 1, and unstable in d = 3 (in fact, we have strong instability by blow-up, see e.g. [START_REF] Cazenave | Semilinear Schrödinger equations[END_REF]). The 2D case is L 2 -critical and instability stems from the fact that the cubic ground state has exactly zero energy (as seen from [START_REF] Weinstein | Nonlinear Schrödinger equations and sharp interpolation estimates[END_REF]). Arbitrarily small perturbations can therefore make the energy negative, which consequently leads to finite time blow-up of the associated solution u by a standard virial argument.

2.2.2. Cubic-quintic case in 2D. For d = 2, it follows from the analysis in [START_REF] Carles | Orbital stability vs. scattering in the cubic-quintic Schrödinger equation[END_REF] that any solitary wave has a mass larger than that of the cubic ground state, i.e. for any ω ∈ (0, 3 16 ), and any solution to

(1.3) φ L 2 (R 2 ) > Q cubic L 2 (R 2 ) ,
where Q cubic is the radial, positive solution to (2.6).

Recall that Q ω denotes the action ground state. Having in mind the Grillakis-Shatah-Strauss theory, the following asymptotic results have been proved in [START_REF] Carles | Orbital stability vs. scattering in the cubic-quintic Schrödinger equation[END_REF][START_REF]The double-power nonlinear Schrödinger equations and its generalizations: uniqueness, non-degeneracy and applications[END_REF]: for ω ≈ 0 or ω ≈ 3 16 , the map ω → M (Q ω ) is increasing. This implies orbital stability in the sense of Definition 1.2, at least for some range of the frequency ω close to the critical values. The numerical plots of M (Q ω ) given in [START_REF]The double-power nonlinear Schrödinger equations and its generalizations: uniqueness, non-degeneracy and applications[END_REF] (see also Section 3 below) suggest that ω → M (Q ω ) is indeed increasing on the whole range ω ∈ (0, 3 16 ), and hence: Conjecture 2.4. In d = 2, all cubic-quintic action ground state solutions are orbitally stable.

In Section 4 we shall give further numerical evidence for this conjecture to be true, by performing several simulations of the time-evolution of perturbed action ground states in 2D.

Cubic-quintic case in 3D.

As established in [START_REF] Killip | Solitons and scattering for the cubic-quintic nonlinear Schrödinger equation on R 3[END_REF][START_REF]The double-power nonlinear Schrödinger equations and its generalizations: uniqueness, non-degeneracy and applications[END_REF], when d = 3, it holds:

(i) On the one hand, as ω → 0, it holds:

M (Q ω ) = 1 √ ω M (Q cubic ) + √ ω 2 Q cubic 6 L 6 (R 3 ) + O ω 3/2 ,
where Q cubic is the positive, radial solution to (2.6). (ii) On the other hand:

lim ω→3/16 M (Q ω ) = lim ω→3/16 ∂M (Q ω ) ∂ω = +∞.
According to the Grillakis-Shatah-Strauss theory, this implies that cubic-quintic action ground states in 3D are unstable near ω min = 0, and orbitally stable near ω max = 3 16 . Numerical plots in [START_REF]The double-power nonlinear Schrödinger equations and its generalizations: uniqueness, non-degeneracy and applications[END_REF] show a U-shaped curve for ω → M (Q ω ). This suggests the existence of a unique unstable branch and a unique stable branch. We shall numerically investigate the nature of instability in this case in Section 5.

Recalling the fact that the set of (constrained) energy minimizers with negative energy is indeed orbitally stable, cf. Theorem 2.2, we shall now show that solutions to (1.3) may have positive energy when d = 3. Indeed, following the approach of [17, Section 2.4], we can rescale Q ω via

ψ ω (x) := 1 √ ω Q ω x √ ω .
The new unknown ψ ω then solves

- 1 2 ∆ψ ω + ψ ω -ψ 3 ω + ωψ 5 ω = 0,
and, as established in [START_REF] Killip | Solitons and scattering for the cubic-quintic nonlinear Schrödinger equation on R 3[END_REF],

ψ ω = Q cubic + O(ω) in H 1 (R 3 ), as ω → 0.
This implies

E(Q ω ) = √ ω 1 2 ∇Q cubic 2 L 2 (R 3 ) - 1 2 Q cubic 4 L 4 (R 3 ) + O(ω) = √ ω Q cubic 2 L 2 (R 3 ) + O(ω) ,
where the last simplification stems from Pohozaev identities for Q cubic (discard the L 6 norms from (2.1) and (2.2)). Therefore, there exists ω 0 > 0 such that

E(Q ω ) > 0, ∀ω ∈ (0, ω 0 ).
Recalling that M (Q ω ) → ∞ as ω → 0, this shows that there exists unstable action ground states with positive energy and arbitrarily large mass. On the other hand, Theorem 2.2 shows that there exists ρ 0 > 0 such that for all ρ > ρ 0 , the minimization problem (1.5) has a solution, and E min (ρ) < 0. In summary this yields:

Theorem 2.5. Not all action ground states in d = 3 are energy ground states.

To our knowledge, this is the first rigorous statement which shows that the two notions of action versus energy ground states, in general, need to be considered as independent (this is also a consequence of [33, Appendix E], as pointed out in [START_REF]The double-power nonlinear Schrödinger equations and its generalizations: uniqueness, non-degeneracy and applications[END_REF], since at least for some values of ω close to zero, ∂M (Q ω )/∂ω < 0). Note that our theorem is consistent with the results of [START_REF] Jeanjean | On global minimizers for a mass constrained problem[END_REF] which establishes that the converse is true, i.e. all energy ground states are action ground states. From the point of view of dynamics, the results above leave open the possibility of having orbitally stable least action ground states, which are not members of E(ρ).

2.2.4.

Previous results in the 1D case. In the case d = 1, it is fairly natural to generalize the nonlinearity in (1.1) to recover features similar to those of (1.1) when d = 2 or 3. More precisely, consider (2.7)

i∂ t u + 1 2 ∂ 2 x u = -|u| p-1 u + |u| q-1 u, with 1 < p < q.
In dimension one, all algebraic nonlinearities are energy-subcritical and explicit solution formulas for φ ω are available in some cases (see in particular [START_REF] Iliev | Stability and instability of solitary waves for one-dimensional singular Schrödinger equations[END_REF]). In addition, the focusing part is L 2 -critical for p = 4. We therefore expect (2.7) to behave similar to (1.1) in d = 2, if we choose p = 4, and similar to (1.1) in d = 3, if p > 4. Indeed, using [START_REF] Grillakis | Stability theory of solitary waves in the presence of symmetry. I[END_REF][START_REF] Iliev | Stability and instability of solitary waves for one-dimensional singular Schrödinger equations[END_REF], it is proved in [START_REF] Ohta | Stability and instability of standing waves for one-dimensional nonlinear Schrödinger equations with double power nonlinearity[END_REF] that for p = 4, all standing waves (not necessarily ground states) are orbitally stable, while for p > 4, some are orbitally stable (for ω ≈ ω max , computed analogously to the value 3 16 for (1.1)), and some are unstable (for ω ≈ 0).

Numerical simulations have addressed the case p > 4, see [START_REF] Buslaev | Simulation of instability of bright solitons for NLS with saturating nonlinearity[END_REF][START_REF] Grikurov | Soliton's rebuilding in one-dimensional Schrödinger model with polynomial nonlinearity[END_REF][START_REF] Sulem | The nonlinear Schrödinger equation, Self-focusing and wave collapse[END_REF]. In particular, [START_REF] Grikurov | Soliton's rebuilding in one-dimensional Schrödinger model with polynomial nonlinearity[END_REF] reports simulations for perturbations of unstable action ground states, showing two possible dynamics: full dispersion, or convergence to another (stable) soliton.

Remark 2.6. In the case d = 1 and p > 4, the conclusion of Theorem 2.5 remains true, using the same proof.

Numerical construction of action ground states

3.1. Numerical algorithm. In this section, we shall discuss a numerical approach for constructing least action ground state solutions to (1.3) in dimensions d = 2 and 3. To this end, we first note that since Q ω is real and radially symmetric, it solves

(3.1) 1 2 ∂ 2 Q ω ∂r 2 + d -1 r ∂Q ω ∂r -ωQ ω + Q 3 ω -Q 5 ω = 0,
where r = |x|. In order to get an equation with regular coefficients (which consequently allows for a more efficient numerical approximation), we introduce the new independent variable

(3.2) s = r 2 ,
in which (3.1) reads

(3.3) 2s ∂ 2 Q ω ∂s 2 + d ∂Q ω ∂s -ωQ ω + Q 3 ω -Q 5 ω = 0.
Since it is known that cubic-quintic ground states are exponentially decreasing (see, e.g., [START_REF] Carles | Orbital stability vs. scattering in the cubic-quintic Schrödinger equation[END_REF]), we choose an s 0 1 such that Q ω (s 0 ) vanishes within numerical precision (which is of the order of 10 -16 here since we work in double precision). Below s 0 = 10 3 , while in the next section we shall also consider examples with s 0 = 10 4 . The numerical task is thus to find a non-trivial solution to (3.3) for given ω ∈ (0, 3 16 ), such that Q ω (numerically) satisfies the homogenous Dirichlet condition Q ω (s 0 ) = 0.

The interval [0, s 0 ] is then mapped via s = s0 2 (1 + ), ∈ [-1, 1] to the interval [-1, 1]. On the latter we introduce standard Chebyshev collocation points n = cos(nπ/N ), n = 0, . . . , N , N ∈ N to discretize the problem. For any given ω > 0 in the admissible range, the function Q ≡ Q ω is consequently approximated via the Lagrange interpolation polynomial P N ( ) of degree N , coinciding with Q at the collocation points, P N ( n ) = Q( n ), n = 0, . . . , N. Similarly, the (radial) derivative of Q is approximated via the derivative of the Lagrange polynomial, i.e.

∂ ∂s Q(s( n )) ≈ P N ( n ).
At the collocation points, this implies ∂ s Q( ) ≈ DQ, since the interpolation polynomial is obviously linear in the n , n = 0, . . . , N . Here, is the vector with components n , D is the Chebyshev differentiation matrix [START_REF] Trefethen | Software, Environments, and Tools[END_REF][START_REF] Weideman | A MATLAB differentiation matrix suite[END_REF] and Q is the vector with components Q(s( n )), n = 0, . . . , N . With the above discretization, equation (3.3) is approximated by a system of nonlinear equations for the vector Q which can be formally written in the form F(Q) = 0. The homogenous Dirichlet condition for s = s 0 is thereby implemented by eliminating the column and the line corresponding to s 0 , cf. [START_REF] Trefethen | Software, Environments, and Tools[END_REF] for more details. This shows that F(Q) = 0 is an N -dimensional system of nonlinear equations for the N components Q(s( n )), n = 1, . . . , N . This system will be solved via a Newton iteration,

(3.4) Q (m+1) = Q (m) -(Jac F(Q (m) )) -1 MQ (m) ,
where Jac F denotes the Jacobian of F, and where Q (m) , m = 0, 1, . . . denotes the m-th iterate. Note, however, that Q = 0 is always a trivial solution to (3.3), which needs to be avoided during the iteration process. In order for this algorithm to converge to our desired, non-trivial solution Q, one needs to identify a suitable initial iterate Q (0) . To do so, we shall apply a tracing or continuation technique as follows: We introduce in (3.3) a deformation parameter α ∈ [0, 1], such that for α = 0 we have only the focusing cubic nonlinearity, while for α = 1 we obtain the full cubic-quintic equation, i.e. we effectively solve

(3.5) 2s ∂ 2 Q ω,α ∂s 2 + d ∂Q ω,α ∂s -ωQ ω,α + Q 3 ω,α -αQ 5 ω,α = 0, α ∈ [0, 1],
instead of only (3.3). The cubic solitons Q ω,α=0 are numerically known, see, e.g. [START_REF] Arbunich | On a class of derivative nonlinear Schrödingertype equations in two spatial dimensions[END_REF][START_REF] Klein | Numerical study of fractional nonlinear Schrödinger equations[END_REF] (no explicit ground state formula exists in dimensions d 2). We can thus solve the discretized equation (3.5) for α = 0.1 and, say, ω = 0.1 via the Newton iteration described above. The resulting solution Q ω=0.1,α=0.1 is then used as an initial iterate for the same equation for α = 0.2, and so on, until the cubic-quintic case α = 1 is reached. In a second step, we use the ground states obtained for ω = 0.1 as an initial iterate for slightly smaller or larger ω's within the admissible range 0 < ω < 3 16 = 0.1875. In this way all examples presented below can be treated.

As an example we show in Fig. 1 the ground state solutions Q ω (r) at ω = 0.1 for the cubic NLS in blue and the cubic-quintic NLS in red. It can be seen that the situation is qualitatively different depending on the spatial dimension. Whereas in 2D, the cubic-quintic ground state has a slightly greater maximum and is slightly faster decaying than Q cubic , in 3D the cubic-quintic ground state has a much smaller maximum and a considerably larger support.

Numerical ground states in 2D.

We first consider the case d = 2 with N = 400 collocation points: In Fig. 2 we show on the left a plot of the ground state function Q ω (r) for various values of ω. It is seen that the maximum of the ground states increases with ω. The solutions also become more localized with increasing ω. On the right of the same figure, we show the L ∞ (R 2 )-norm of the ground states as function of ω. For convenience, we only consider values of ω ∈ [0.005, 0.16].

In Fig. 3 we depict the ground-state mass M (Q ω ) and energy E(Q ω ) as a function of ω. These plots are based on a total library of roughly 100 numerical ground state solutions Q on the shown range of ω. The corresponding mass-and energy-integrals are thereby computed with the Clenshaw-Curtis algorithm in s = r 2 , a spectral integration method based on the same Chebyshev collocation points as before, see [START_REF] Trefethen | Software, Environments, and Tools[END_REF]. Both M (Q ω ) and E(Q ω ) appear to be monotonic in ω. In particular, the 3.3. Numerical ground states in 3D. In the case d = 3, we use the same numerical parameters as before: In Fig. 4, we show on the left the ground states for several values of ω. It can again be seen that the maximum of Q ω increases with ω, at least up to some value ω * ≈ 0.1. For larger values of ω, however, the L ∞ (R 3 )-norm of Q ω is seen to be decreasing again.

Analogously to the 2D case, the solutions become more localized with increasing ω. Note, however, that despite its exponential decay, the 3D soliton is less localized than in the case of the purely focusing, cubic NLS, see Fig. 1 on the right. The 3D ground states in Fig. 4 on the left are also found to be less peaked than the corresponding solutions in dimension 2, see Fig. 2.

In contrast to the 2D case, the ground state mass M (Q ω ) is no longer monotonically increasing as a function of ω. Looking at Fig. 5, we see that, instead, M (Q ω ) has a minimum at ω crit ≈ 0.026. We consequently expect orbital instability of ground states Q ω for ω < ω crit , a phenomenon we shall study in more detail in Section 5. In Fig. 6, the corresponding ground state energy E(Q ω ) is seen to have a maximum at the same ω crit ≈ 0.026. are different from (1.1)). Clearly, ground states φ ω corresponding to the upper branch cannot correspond to constrained energy minimizers with mass ρ = M (Q ω ).

Orbital stability of action ground states in 2D

4.1. Numerical method for the time-evolution. In this section we will numerically study the time-evolution of (1.1) resulting from initial data u 0 given by perturbations of action ground states. We will only consider perturbations which conserve the radial symmetry. To consider more general perturbations, a full 3D code would be necessary which is beyond the scope of the present paper. This allows us to use the change of variables (3.2) and effectively solve (1.1) in the (non-singular) form (4.1)

i∂ t u + 2s ∂ 2 u ∂s 2 + d ∂u ∂s + |u| 2 u -|u| 4 u = 0, d = 2, 3.
We thereby use the same discretization for s ∈ [0, ∞) in terms of Chebyshev collocation points as detailed in the previous section.

After the spatial discretization in s, equation (4.1) is then approximated via a system of ordinary differential equations. These equations are then integrated in time using a time-splitting method in which the linear step is solved numerically via an implicit fourth-order Runge-Kutta method, see [START_REF] Klein | Numerical study of the transverse stability of the Peregrine solution[END_REF] for more details. The accuracy of this time-integration algorithm is henceforth controlled via the analytically conserved quantity E(u), which in our case nevertheless depends on time due to unavoidable numerical errors. As discussed in [START_REF] Klein | Fourth order time-stepping for low dispersion Korteweg-de Vries and nonlinear Schrödinger equations[END_REF], the numerical conservation of the relative mass tends to overestimate the numerical error by one to two orders of magnitude. We shall always aim at a numerical error below O(10 -3 ), i.e. below plotting accuracy. This means that we ensure a relative energy-conservation

∆ E = E(t) E(0) -1 of order ∆ E = O(10 -5
), or better. We shall use a single computational domain Ω = [0, s 0 ] for which we impose a homogeneous Dirichlet condition u(t, s 0 ) = 0, for all t 0. We mostly choose s 0 = 10 3 , but in some unstable situations we shall also take s 0 = 10 4 . As a basic test case, we first propagate the three-dimensional ground state Q ω , numerically found at ω = 0.1. We thereby use N t = 10 4 time-steps until a final time t f = 10. We find that the hereby obtained numerical solution u, at t = t f , satisfies max

Ω |u(t f , •) -e it f ω Q ω | = O(10 -9 ),
i.e. the same order of accuracy as reached in [START_REF] Klein | Numerical study of the transverse stability of the Peregrine solution[END_REF].

Remark 4.1. In general, it is not unproblematic to work with a homogeneous Dirichlet boundary condition on a finite numerical domain, since this could lead to unwanted reflections of the emitted radiation at the boundary, see, e.g., the discussion in [START_REF] Birem | Multidomain spectral method for Schrödinger equations[END_REF]. In our case, however, only small, rapidly decreasing perturbations of ground states are considered. It is thus possible to work on sufficiently large computational domains Ω, on which the radiation can separate from the bulk before spurious reflections from the boundary lead to noticeable effects. 4.2. Time-evolution of perturbed 2D ground states. In this subsection, we shall study the time-evolution of perturbed ground states to (1.1) in dimension d = 2. To this end, we first consider the case where

(4.2) u 0 (x) = λQ ω (x), λ > 0.
Here λ > 0 is a perturbation parameter and Q ω is a numerically obtained action ground state, at a certain admissible frequency ω ∈ (0, 3 16 ). We first study the case where ω = 0.1 and λ = 0.99, and use N t = 10 4 time steps to reach the indicated final time t f = 20. As expected, the solution u to (1.1), effectively given by (4.1), is found to be close to the exact time-periodic state

φ ω (t, x) = e iωt Q ω (x).
To this end, we show on the left of Fig. 8 the L ∞ (R 2 )-norm of the solution as a function of time. It can be seen that it approaches a final state as t → 20. The latter is found to be very close (in absolute value) to the unperturbed ground state Q ω=0.1 . Note that the L ∞ -difference is of the order of 10 -4 and thus, much smaller than the initial perturbation. As a second case, we consider the same 2D initial data (4.2), but with λ = 1.001. In Fig. 9 we again show the L ∞ -norm of the solution as a function of time. Similarly as before, a final state is reached and its maximum is again found to be very close to the unperturbed ground state Q ω=0.1 . In both cases, we find that the difference is largest for r close to the origin.

In order to illustrate that the qualitative picture found before is not due to our specific choice of perturbations, we shall also consider ground states perturbed by a small Gaussian-like perturbation, i.e.

(4.3)

u 0 (x) = Q ω=0.1 (x) ± λe -|x| 2 , λ = 0.001.
Note that we only consider smooth perturbations in this paper in order to allow for spectral accuracy in the radial coordinate, i.e., an exponential decrease of the numerical error with the number of collocation points. In Fig. 10 we show the behavior in time of the respective L ∞ -norms for the two choices ±λ. In both situations the difference between |u| at t f = 20 and Q ω=0.1 is found to be of the order O(10 -4 ). Moreover, the error (not depicted here for the sake of readability) is again found to be largest close to the 3) and ω = 0.1. On the left the case with "+" sign. On the right the case with "-" sign.

If similar perturbations are applied to other ground states Q ω , the resulting solution u behaves qualitatively similarly. Our numerical tests therefore support Conjecture 2.4. However, we also find that the smaller the choice of ω ∈ (0, 3 16 ), the longer it takes for the solution u to reach its final state. In fact, for small enough ω, damping effects within the time-oscillations of |u| become almost invisible, even if one computes up to much larger times t f = 400, see Fig. 11. 3) and ω = 0.05. On the left the case with the "+" sign. On the right the case with the "-" sign.

(In-)stability of action ground states in 3D

5.1. Stable branch. In this section, we shall study the question of (in-)stability of cubic-quintic ground states in dimension d = 3. In view of Fig. 5, we expect ground states Q ω with ω > ω crit ≈ 0.026 to be orbitally stable. That this is indeed the case, is strongly suggested by our numerical results below.

To this end, we first consider multiplicative perturbations of Q ω on the stable branch: In Figures 12 and13 we study the time-evolution of (1.1) with initial data of the form (4.2). On the left of Fig. 12 In Fig. 13 we study the analogous situation with λ = 1.001: Again, the (absolute value of the) solution u seems to settle around t f = 15 on the stable unperturbed ground state Q 0.1 . In both cases, the error between |u| and Q ω=0.1 is again found to be of the order O(10 -4 ). In Fig. 14 we show the solution u to (1.1) obtained from initial data (4.2) with ω = 0.01 and λ = 0.999. Note that this implies M (u 0 ) < M (Q ω ). The solution is seen to be purely dispersive which is also confirmed by the L ∞ -norm of the solution as a function of time (depicted on the right of the same figure). In fact, we did not discover any stable structure within the time-evolution even if we let the numerical code run for longer times. If we consider the same ground state before, but instead choose λ = 1.001, we find a different kind of instability. Now the solution u shows oscillations of high amplitude, see Figure 15.

These oscillations are even more visible in the L ∞ -norm of the solution, as depicted on the left of Fig. 16. One can see that early on the norm is growing strongly but then it appears to show damped oscillations around some final state. We conjecture that the latter corresponds to another ground state on the stable branch. To this end, we compare the maximum of |u|, obtained at the final time t f = 500, with the L ∞ -norms in our library of previously computed action ground states Q ω , cf. Fig. 4. Indeed we find good agreement of |u|, when compared to Q ω with ω = 0.047 > ω crit , see the right of Fig. 16. Thus perturbations of unstable ground states where M (u 0 ) > M (Q ω ), seem to result in solutions which eventually settle on another, stable ground state as t → +∞. Note, however, that M (Q ω=0.01 ) ≈ 79.44 while M (Q ω=0.047 ) ≈ 77.05. If the final state had the same mass as the unperturbed initial state, this would correspond to an ω ≈ 0.0495. This shows that a certain amount of mass is lost through radiation.

5.3.

Other kinds of perturbations. The results described above are not due to our specific choice of perturbations. To show this, we consider initial data (5.1) u 0,± (x) = Q ω (x) ± λe -(|x|-|x0|) 2 , λ = 0.001.

For both |x 0 | = 0 and |x 0 | = 0, the solution in the case with the "+" sign looks very similar to the one depicted in Fig. 16. This fact becomes particularly clear when one compares the time-evolution of the L ∞ -norm of u depicted in Fig. 17, with the one from Fig. 16. By comparing the maximum of |u| found at the final time t f = 600 with the L ∞norm of a stable ground state, we find good agreement with Q ω=0.048 . The latter has mass M (Q ω=0.048 ) ≈ 77.95. Unfortunately, we are unable to clearly decide whether the final state is closer to Q ω=0.048 than to Q ω=0.047 .

In the case of initial data (5.1) with the "-" sign, we again find that the solution is completely dispersed, see Fig. 18. This is consistent with our earlier findings above which indicate that perturbation with M (u 0 ) < M (Q ω ) lead to purely dispersive solutions. We finally note that the situation is qualitatively similar for other values of ω on the unstable branch. In our last example, we choose ω = 0.007 < ω crit within the initial u 0 given by (5.1). We again find that perturbations with an initial mass smaller than M (Q ω=0.007 ) are purely dispersive. Perturbations with mass larger than M (Q ω=0.007 ) lead to damped oscillations around some final state, see Fig. 19. The asymptotic final state appears to be close to Q ω=0.044 on the stable branch. The mass of the unperturbed initial data M (Q ω=0.007 ) ≈ 90.57 is seen to be bigger than M (Q ω=0.044 ) ≈ 75.37, showing again that a non-negligible part of the initial mass has been radiated away.

The numerical results within this subsection can then be summarized as follows: (i) If M (u 0 ) < M (Q ω ), then the solution u to (1.1) is purely dispersive; (ii) If M (u 0 ) > M (Q ω ), then the solution to (1.1) converges, as t → +∞, to a solitary wave φ ω (t, x) = e iωt Q ω (x) plus radiation, where Q ω is a stable ground state with mass smaller that the unstable one, M (Q ω ) < M (Q ω ).

Remark 5.2. The same instability scenario was found (numerically) for perturbed solitary wave solutions to the generalized BBM equation in [START_REF] Bona | Stable and unstable solitary-wave solutions of the generalized regularized long-wave equation[END_REF], and for a version of NLS with derivative nonlinearity in [START_REF] Arbunich | On a class of derivative nonlinear Schrödingertype equations in two spatial dimensions[END_REF]. In particular, analogously to our situation, a perturbation which lowered the mass of the initial data below the one of the (unstable) solitary wave always resulted in purely dispersive solutions. In all of these cases, it remains an interesting open question to find a possible selection criterion for the specific value ω which describes the (stable) asymptotic state φ ω .

1 . Introduction 1 . 1 .

 111 Basic setting. This work is concerned with the time-evolution corresponding to the cubic-quintic nonlinear Schrödinger equation (NLS)(1.1) i∂ t u + 1 2 ∆u = -|u| 2 u + |u| 4 u, (t, x) ∈ R × R d ,in dimensions d = 2, or d = 3, and subject to initial data

Figure 1 .

 1 Figure 1. Ground state solutions Q ω=0.1 to the cubic NLS in blue and the cubic-quintic NLS in red: on the left for d = 2 and on the right for d = 3.

Figure 2 .

 2 Figure 2. Left: Ground state solutions to (3.1) in dimension d = 2 for several values of ω. Right: The L ∞ -norm of these states as a function of ω.

Figure 3 .

 3 Figure 3. M (Q ω ) and E(Q ω ) as a functions of ω in dimension d = 2.

Figure 4 .Figure 5 .Figure 6 . 8 EFigure 7 .

 45687 Figure 4. Left: Ground state solutions to (3.1) in dimension d = 3 for several values of ω. Right: The L ∞ -norm of these states as a function of ω.

Figure 8 .

 8 Figure 8. The 2D solution u to (1.1) for initial data (4.2) with ω = 0.1 and λ = 0.99. On the left the L ∞ -norm as a function of time. On the right the difference between |u| and Q ω=0.1 at the final time t f = 20.

Figure 9 .

 9 Figure 9. The 2D solution u to (1.1) for initial data (4.2) with ω = 0.1 and λ = 1.001. On the left the L ∞ -norm as a function of time. On the right the difference between |u| and Q ω=0.1 at the final time t f = 20.

Figure 10 .

 10 Figure 10. L ∞ -norm as a function of time for the 2D solution u to (1.1) with initial data (4.3) and ω = 0.1. On the left the case with "+" sign. On the right the case with "-" sign.

Figure 11 .

 11 Figure 11. L ∞ -norm as a function of time for the 2D solution u to (1.1) with initial data (4.3) and ω = 0.05. On the left the case with the "+" sign. On the right the case with the "-" sign.

  we show the L ∞ -norm of the solution u obtained in the case ω = 0.1 and λ = 0.99. On the right of the same figure, we show the difference between the unperturbed ground state Q ω=0.1 and |u| at the final time t f = 15. It can be seen that the L ∞ norm settles on a nearly constant value as t → 15.

Figure 12 .

 12 Figure 12. The 3D solution u to (1.1) for initial data (4.2) with ω = 0.1 and λ = 0.99. On the left the L ∞ -norm as a function of time. On the right the difference between |u| and Q ω=0.1 at the final time t f = 15.

Figure 13 .

 13 Figure 13. The 3D solution u to (1.1) for initial data (4.2) with ω = 0.1 and λ = 1.001. On the left the L ∞ -norm as a function of time. On the right the difference between |u| and Q ω=0.1 at the final time t f = 15.

Figure 14 .

 14 Figure 14. Left: The absolute value of the solution u to (1.1) in dimension d = 3 obtained from initial data (4.2) with ω = 0.01 and λ = 0.999. Right: The L ∞ -norm of u as a function of time.

Figure 15 .

 15 Figure 15. Solution to (1.1) in dimension d = 3 for initial data (4.2) with ω = 0.01 and λ = 1.001.

Figure 16 .

 16 Figure 16. Solution to (1.1) in dimension d = 3 for initial data (4.2) with ω = 0.01 and λ = 1.001: On the left the L ∞ -norm of u as a function of time. On the right |u| at the final time (in blue) together with the ground state Q ω=0.047 (in green).

Figure 17 .

 17 Figure 17. Solution to (1.1) in dimension d = 3 for initial data of the form (5.1) with |x 0 | = 1, and ω = 0.01: On the left the L ∞ -norm of u as a function of time. On the right |u| at the final time (in blue) together with the ground state Q ω=0.048 (in green).

Figure 18 .

 18 Figure 18. The L ∞ -norm of the solution u in dimension d = 3 obtained from initial data (5.1) with the "-" sign. On the left the case with x 0 = 0 and on the right the one with x 0 = 1.

Figure 19 .

 19 Figure 19. Left: Solution to (1.1) in dimension d = 3 for initial data (5.1) with x 0 = 0 and ω = 0.007. Right: The L ∞ -norm of u as a function of time.

Conjecture 5 . 1 .

 51 For ω < ω c , consider initial data of the formu 0 (x) = Q ω (x) + (|x|), with

This publication is based on work supported by the NSF through grant no. DMS-1348092. RC acknowledges support from Rennes Métropole, through its AIS program. CK is partially supported by the ANR-FWF project ANuI -ANR-17-CE40-0035, the isite BFC project NAANoD, the EIPHI Graduate School (contract ANR-17-EURE-0002), by the European Union Horizon 2020 research and innovation program under the Marie Sklodowska-Curie RISE 2017 grant agreement no. 778010 IPaDEGAN and the EITAG project funded by the FEDER de Bourgogne, the region Bourgogne-Franche-Comté and the EUR EIPHI.


Tatiana Duigou
email: tatiana.duigou@cea.fr

Fabien Chabuel

Gilles Dennler

Pascal Francescato

Julien Gaume

Georges Habchi

M Lagache

P Saré

Lionel Tenchine

VIPV: Thermocompression process development and simulation to integrate

Keywords: PV Applications, Integration and Storage;Sub-topic: 6.2. PV in Infrastructure, PV cell integration, process modelling, ground transportation

come   L'archive ouverte pluridisciplinaire

Introduction

Greenhouse gas emissions have almost doubled since 1970 and the transport sector emits 14% of these emissions. In this context, the automotive sector is encouraged to propose alternative solutions to gaspowered vehicles. The current trend to tackle this issue is to electrify vehicles, usually using high capacity batteries. However, these popular technologies of Electric Vehicles (EVs) require a large amount of lithium for their heavy batteries. Thus, the integration of PhotoVoltaic (PV) modules as body parts for hybrid or electric vehicles could allow to extend the vehicle's range while reducing the battery size [START_REF] Araki | To do list for research and development and international standardization to achieve the goal of running a majority of electric vehicles on solar energy[END_REF].

Topics and investigations

The automotive industry is changing to lightweight materials and more energy-ecient vehicles. PV modules have a potential market in this transition. Rigid PV panels integrating silicon cells built directly into the car body outperform thin exible solar lms in terms of eciency and reliability. This rigidity is usually brought by glass. To increase the integration of PV modules into the car, the use of composite materials is an interesting option, as they oer an ecient combination of mechanical performances and lightweight. However, the brittle behaviour of the cells and the necessary achievement of certain optical characteristics are incompatible with the usual manufacturing processes of composite panels yet. As well, lamination restricts the possibilities of integration in automotive parts. These parts may have complex shapes like double curvature and stieners which are poorly compatible with lamination. Thus, this study focuses on the development of an alternative manufacturing process for the production of composite double-curved PV modules and the study of thermomechanical stresses resulting from this process. A previous part of this research work enabled to identify a process which oers an alternative to the lamination process: the thermocompression. It allows to produce curved photovoltaic modules made of composites, and adapts to the automotive industry constraints. This study tackles both the adaptation of this process to integrate photovoltaic cells, and the simulation of the thermocompression process with a nite element calculation tool. The numerical model aims at optimizing the process parameters and module stack to limit residual stresses while minimizing the weight of the structure.

Conclusions and perspectives

This study focused on the adaptation of a thermocompression process to integrate photovoltaic cells. It considers both the experimental process development, and its simulation with a nite element analysis software. PV modules produced with this process were used for the experimental validation of the numerical model. The experimental tests validate the thermomechanical multistep model. This predictive model allows to optimize the composite layout and process parameters to minimize residuals stresses. Complementary experimental tests evaluate the electrical performances and ageing behaviour of the PV module with regards to IEC 61215 standard. According to our results, a composite structure integrating photovoltaic cells appears to be an interesting alternative to traditional glass structures, both in terms of mass and compatibility with the constraints of the automotive industry. Thus, the process developed could be applicable to either the automotive industry or any application needing lightweight and curved PV modules. Further research will aim at proposing an industrial scheme for this process.

Explanatory pages

Aims

The main obstacles to Vehicle Integrated PV (VIPV) are:

• The need to maximize the photovoltaic modules power output.

• Like most body parts, the structure integrating PV cells may present a double curvature to which monocrystalline silicon cells are poorly adapted [START_REF] Araki | To do list for research and development and international standardization to achieve the goal of running a majority of electric vehicles on solar energy[END_REF].

• The mechanical requirements for automotive parts: the body parts need to withstand a set of loads to which PV cells are not suited. The design and process of integration should take into account the minimization of stresses for the cells.

• The weight of the PV module: as it is supposed to provide energy to the vehicle, the integration of photovoltaic cells should not penalize the weight of the vehicle to maximize the vehicle range allowed by the energy intake of the PV module. The use of an optimized composite structure can make the achievement of this aim easier.

• The cost of production of automotive body parts, even if they integrate more functionalities.

Technologies of solar-powered cars already exist but they mostly use glass encapsulation to protect the cells, thus not minimizing the weight and optimizing the mechanical performances of the module. To keep a competitive product, it is further envisaged to design the composite structure integrating PV through a "one-shot" process. It should guaranty the previous specications and minimize the strain and residual stresses due to the heterogeneity and non-symmetry of the structure [START_REF] González | Structural composites for multifunctional applications: Current challenges and future trends[END_REF]

[3][4].
This study considers a composite structure with double-curvature integrating silicon solar cells. The aim is both to adapt the thermocompression process to the integration of Si cells in composite materials and to propose a predictive model of the residual stresses implied by this process. The numerical model enables to optimize the stack and processes parameters in order to minimize residual stresses, minimize the structure weight and optimize the mechanical performance of the structure. This model is validated with mechanical experimental tests on parts produced by the thermocompression process developed in this study.

Complementary tests aim at validating the electrical performances and ageing behaviour of the PV modules with regards to IEC 61215 standard.

Methods and results

a) Process development

A preliminary study was made to select materials adequate for the thermocompression process. First experimental tests were operated with a planar mold and enabled to identify materials and process parameters allowing to produce composite PV modules presenting no electrical failure. This work fed a set of specications for the design of a double curvature mold compatible with the thermocompression process and enabling to integrate PV cells. Further work enabled to produce double curved PV modules showing no cell failure. After these encouraging results, it was necessary to optimize the process parameters and stack to minimize weight and the internal stresses resulting from the process. A Design Of Experiment on this validated simulation work enabled identifying the best stack candidates and process parameters.

b) Numerical model

The numerical model is designed with a Finite Element Analysis software -ANSYS-to describe a module integrating the solar cells in sheets of polymer encapsulant. The front sheet is either made of a polymer sheet, or of glass-bre pre-pregs. The back sheet is made of glass bre-based thermoplastic pre-pregs. These materials allow to withstand the loads and protect the cells from the environment.

Due to symmetries and periodicity of the structure, loads and boundary conditions, the representative pattern is composed of one cell. This is representative of a 4-cell module. The model uses hexahedral solid elements with adapted sizing (g. 1a). The setup is adapted for a parametric analysis exploring the inuence of materials, thickness, ply orientations, number of plies ... on the stiness and strength of the whole structure. A multi-step thermomechanical analysis modelling the whole process is performed. A rst analysis was led to identify the simplications which could be applied, passing from the physical steps of the process to their numerical implementation (g. 2). It simulates 3 steps of the process: This parametric model enables an optimization of this non-symmetric structure to full the following objectives:

• The minimization of residual thermomechanical stresses resulting from the process. This process implies residual stresses coming from the curvature of the cells, and dierential deection between the front and back sheets -especially in the case of a polymer front sheet during the cooling phase.

• A representative thermal load, representing a thermal cycle of the international standard for PV modules IEC 61215;

• Equivalence of the bending stiness to that of a car roof;

• Mechanical strength of the cell and of the interconnections.

This exploration is led with a design of experiment method. Among others, it aims at minimizing interfacial and internal stresses (g. 1b).

The parametrical aspect of this simple model enables to scan a large number of congurations and to identify the most critical parameters among the number of plies, their orientation, the nature of the resin and of the encapsulant. The double-curved modules have been produced with a thermocompression process. Various materials have been used, including:

• Glass bre thermoset and prepregs fabrics for the front and back sheets;

• Thermoplastic polymers for the front sheet;

• Dierent thermoplastic encapsulants.

Several tests have been performed.

• The rst one aims at validating the deformed shape of the module at the outlet of its processing by means of laser metrological measures. This test enables to validate the numerical model, especially the deformation of the stack after the process.

• The second ones consist in mechanical tests. They enable to evaluate the bending stiness of the overall structure and the load for the rst damage. Impact tests were also performed. These test validate the mechanical behavior of the module for some tests representative of a few automotive qualication tests.

• Complementary electrical tests were performed to assess the reliability of the PV modules produced by thermocompression, especially their electrical behaviour under bending and during ageing tests. The structure was qualied with standard photovoltaic qualication equipment, including spectrophotometry, electroluminescence and light I-V measures before and after tests in climatic chambers (norm IEC 61215). They enable to evaluate and optimize the power supply of the panel, choosing the most appropriate materials for the front sheet.

Figure 1 :

 1 Figure 1: a) Z deformation and mesh for a double curved cell. b) Stresses in the cell resulting from the cooling in the mold.

1 .

 1 Application of the double cuvature to the cells, front and back sheet materials being inactive because they are not processed at this step. The thermal condition is 230 • C. 2. Cooling of the stack in the mold. The front and back sheet plies are active and have temperature dependent behaviors. The temperature goes down from 230 • C to 50 • C. 3. Opening of the mold at 50 • C. This step enables to quantify the residual deection and stresses coming from the thermocompression process.

Figure 2 :

 2 Figure 2: Going from process steps to the numerical implementation of the multi-step simulation.

Figure 3 :

 3 Figure 3: Double-curved PV module of 40 cm x 40 cm produced by thermocompression.


Ali Bakiri

Jane Hubert
email: jane.hubert@univ-reims.fr

Romain Reynaud

Carole Lambert

Agathe Martinez

Jean- Hugues Renault

Jean-Marc Nuzillard

Reconstruction of HMBC Correlation Networks: A Novel NMR-based Contribution to Metabolite Mixture Analysis

Here is introduced a new in silico method for the dereplication of natural metabolite mixtures based on HMBC and HSQC spectra that inform about short-range and long-range H-C correlations occurring in the carbon skeleton of individual chemical entities. Starting from the HMBC spectrum of a metabolite mixture, an algorithm was developed in order to recover individualized HMBC footprints of the mixture constituents. The collected H-C correlations are represented by a network of NMR peaks connected each other when sharing either a 1 H or 13 C chemical shift value. The obtained network is then divided into clusters using a community detection algorithm, and finally each cluster is tentatively assigned to a molecular structure by mean of a NMR chemical shift database containing the theoretical HMBC and HSQC correlation data of a range of natural metabolites. The proof of principle of this method is demonstrated on a model mixture of 3 known natural compounds and then on a real-life bark extract obtained from the common spruce (Picea abies L.).

INTRODUCTION

Identification and quantification of chemical entities directly within complex mixtures is gaining a great attention in the field of natural product research and more generally in all metabolomics studies dealing with the analysis of low molecular weight metabolites in biological systems (Wist, 2016) (Gaudêncio & Pereira, 2015) (Bingol et al., 2016) (Kim, Choi, & Verpoorte, 2011). Over the last decade, a range of new analytical strategies have emerged to accelerate metabolite identification processes, mainly based on the so-called dereplication concept (Beutler, Alvarado, Schaufelberger, Andrews, & McCloud, 1990)ref revue wolfender (Gaudêncio & Pereira, 2015;Hubert, Nuzillard, & Renault, 2015). Among the panel of analytical techniques available, Nuclear Magnetic Resonance (NMR) remains by far the most efficient method to unambiguously identify small organic molecules (Smolinska, Blanchet, Buydens, & Wijmenga, 2012) (Schripsema, 2010). NMR has the unique ability to probe structural information at the atom level for NMR-active nuclei such as 1 H, 13 C, 15 N and 31 P.

Proton ( 1 H nucleus) NMR spectroscopy is mostly used due to its high sensitivity, but suffers from frequent signal overlaps which confuse data interpretation when dealing with samples of complex composition. Conversely, 13 C NMR is less sensitive than 1 H NMR, nonetheless it provides a larger chemical shift dispersion (240 ppm) that reduces signal overlaps and leads to well-resolved spectra containing individualized narrow peaks. Along with the ability to obtain information about the carbon skeleton of all detected organic molecules, these advantages have progressively prompted 13 C NMR in the core of recent dereplication procedures ( Clendinen et al, 2014;Laude et al, 1986;Bighelli et al, 1994;Hubert et al, 2014).

In addition to 1 H and 13 The second breakthrough in mixture analysis strategies is the application of mathematical treatments to facilitate the interpretation of complex analytical data (Wiklund, n.d.). Along with the multivariate statistical analyses classically used in metabolomics, new mathematical methods such as data fusion (Acar, Bro, & Smilde, 2015) and graph theory are more and more applied to rise up the relevance of complex mixture studies (Garg et al., 2015). For instance, molecular networking has recently emerged in the context of mass spectrometry (MS) as an efficient way to handle the complexity of MS/MS data. It relies on the observation that structurally similar metabolites share similar MS/MS fragmentation patterns. On the basis of this property, the correlations between MS/MS spectral data of a natural sample are measured and used to construct a molecular network in which each cluster of similar fragments represents a specific group of structurally close metabolites. (Winnikoff, Glukhov, Watrous, Dorrestein, & Gerwick, 2014) (Allard et al., 2016;Yang et al., n.d.) Here is presented a new HMBC-based dereplication method that uses a networking approach for the deconvolution of complex NMR spectra of metabolite mixtures. This method exploits the ability of HMBC experiments to provide connectivity information between 1 H and 13 C atoms located at their vicinity. Starting from the HMBC spectrum of a metabolite mixture, the strategy aims to build the network of 1 H-13 C HMBC correlations and to highlight individualized patterns representing specific molecular fragments or even whole chemical structures. The isolated patterns are then identified by database search.

The proof of concept of this HMBC-based atom networking approach is demonstrated here on a synthetic model mixture of 3 known natural compounds: (+)-catechin, taxifolin and Eisorhapontin, and then on a genuine extract obtained from the bark extract of Picea abies.

METHODS

Sample preparation.

The first analyzed sample was a model mixture composed of 3.9 mg of taxifolin, 0.4 mg of (+)-catechin, and 0.7 mg of E-isorhapontin (all purchased from Sigma-Aldrich, Steinheim, Germany). The three compounds were directly dissolved in 550 µL of deuterated dimethyl sulfoxide (DMSO-d6, Eurisotop, Saint-Aubin, France). The second sample was a crude methanol extract obtained from the barks of P. abies, prepared as described previously (Angelis et al., 2016). Briefly, 2 kg of barks were collected in the French Champagne-Ardenne territory in 2014 within the context of forestry activities (Voucher specimen JH-2014-9, Faculty of Pharmacy, University of Reims Champagne-Ardenne, France). After drying at 30 °C for 72 h, the barks were powdered in a hammer mill and extracted with methanol under magnetic stirring at ambient temperature for 20 h. After filtration, methanol was evaporated to dryness and an aliquot of the resulting extract (20 mg) was dissolved in 550 µL of DMSO-d6 for NMR analyses.

NMR analysis and peak picking. NMR analyses were performed at 298 K on an Avance III-600 spectrometer (Bruker, Karlsruhe, Germany) equipped with a cryoprobe optimized for 1 H detection and with cooled 1 H, 13 C and 2 H coils preamplifiers. HSQC and HMBC spectra were recorded using standard Bruker pulse programs. For the acquisition of HSQC spectra (pulse sequence hsqcedetgpsisp2.2), the size of FID was set at 1024 in the F2 dimension and 512 in the F1 dimension, and the number of scans was set at 4. For the acquisition of HMBC spectra (pulse sequence hmbcetgpl3nd), the size of FID was set at 2048 in the F2 dimension and 512 in the F1 dimension, and the number of scans was 4. The spectra were manually phased and baseline corrected using the TOPSPIN 3.2 software (Bruker). Once the HMBC spectra recorded and processed, a semi-automated peak picking was performed to obtain the list of correlations. The peak picking was performed on the different regions of the spectrum first automatically, then manually in order to remove possibly meaningless peaks caused by thermal noise, t1 noise, or residual 1 J peak pairs. The automatic peak picking process was based on the minimum intensity set to the lowest contour level (Bruker processing window in Topspin 3.2).

Prediction of HMBC and HSQC data.

A database of theoretical HMBC connectivity data of a range of natural metabolites was created as follows: the 1 H and13 C NMR chemical shift values of  2700 natural metabolites (October 2017) were calculated using the NMR Workbook Suite ACD/Labs software (Ontario, Canada). From these data, HMBC correlations (pairs of 1 H and molecule of the database, this program creates a list of all H-C atom pairs that might be visible in an HMBC spectrum, i.e. all H-C pairs in which atoms are placed 2 or 3 bonds away. The correlation between chemical shifts of atoms separated by 4 bonds are not observed in all molecular structures and were therefore not taken into account. Once this list is established, the corresponding predicted chemical shifts values of each atom pair are included in the database (Fig. 1). In the same way, the HSQC correlations of all metabolites of the database were determined for the 1 H-13 C atom pairs located at a distance of 1 bond. In addition to the 1 H and 13 C chemical shifts values, the parity of the number of H atoms attached to each C atom was also considered. Thus, an additional column was added to the HSQC correlation table in order to store the parity data. From experimental HMBC to 1 H-13 C chemical shift network. The 1 H-13 C chemical shift networks were generated from the peak list of the experimental HMBC spectra using an algorithm written in Python. The principle of this algorithm is to extract the HMBC correlation data from the experimental HMBC peak list and then to create a network edge list by connecting each peak (vertex) to its adjacent neighbors on the same row of the HMBC spectrum in the F1 dimension (having the same 13 C NMR chemical shift) or on the same column in the F2 dimension (having the same 1 H NMR chemical shift). The obtained edge list is then converted to a network graph using the igraph 0.7.1 package for Python 2.7. The clustering of network vertices was performed using the Louvain community detection algorithm (Blondel et al.). The Louvain method is a modularity optimization-based algorithm.

Modularity is a quality function that measures the strength of connectivity between the nodes of the same community. It measures the density of edges inside the community compared to edges outside the community. The modularity is defined as followed (J. Newman):

Q = ∑ (eii -ai 2 )
eii : The fraction of edges that connects vertices from the same group. eij : The fraction of edges that connects vertices from group i to vertices from group j and ai =

 eij

The Louvain algorithm is a bottom-up algorithm that starts with isolated vertices as communities that will be merged progressively into meta-communities according to the gain of modularity until an optimal graph partition is reached (Blondel et al.). The "RBER" method is a variant of the Louvain algorithm that uses Erdös-Rényi null-model (Newman, M et al 2004), in which each edge has the same probability of appearing. It includes a resolution parameter that allows the control of the community's size, in general, a higher resolution parameter leads to smaller communities and vice versa. The RBER method was applied for the detection of the HMBC peak clusters in HMBC networks, with a resolution parameter set at 0.2. This value was selected by trial and error and found to give the best results.

Metabolite identification.

Metabolite identification was performed by comparing the theoretical HMBC correlations of the database metabolites to the experimental HMBC correlations of the clusters extracted by network analysis. For this purpose, an algorithm was written in Python language to obtain the lists of experimental correlations corresponding to each cluster of the HMBC correlation network. Then, these lists of correlations are compared to the theoretical HMBC correlations of database metabolites. For each theoretical HMBC correlation of the database, the algorithm attempts to find an experimental correlation that matches both 1 H and 13 C chemical shift values within a defined tolerance radius, typically set at 1.5 ppm for 13 C NMR chemical shifts and 0.15 ppm 1 H NMR chemical shifts. The percentage of theoretical correlations that match experimental correlations for each metabolite of the database is then calculated to evaluate its likeliness to be present within the analyzed mixture.

Finally, the algorithm returns a list of metabolites for each cluster, ranked in the decreasing order of matching percentage (Fig. 2). As signals belonging to a putatively identified metabolite may be divided between several clusters, a second step is then performed to compare the theoretical HMBC correlations of the top ranked 20 molecules of each cluster to the complete experimental list of HBMC correlations. By this way, an improved list of molecules ranked according to their decreasing percentage of matching to the whole HMBC spectrum is obtained. The resulting list is submitted to an additional validation step which consists in the comparison of the theoretical HSQC correlations of the top ranked molecules to the HSQC correlations detected experimentally. This step is performed using another locally developed algorithm which compares the 1 H and 13 C chemical shifts values of the HSQC correlations together with their parity, therefore theoretical correlations are matched to the experimental one only if both their chemical shifts and parity match. From this step, an additional score (percentage of theoretical HSQC correlations that match experimental HSQC correlations) is obtained for each proposed metabolite. This score reinforces the identification accuracy.

RESULTS

Proof of concept on a model mixture. A mixture of 3 commercially available compounds (2 flavonoids and 1 stilbene), taxifolin, (+)-catechin and E-isorhapontin, was used as a model in order to establish a first proof of concept of the method. The proportions of the mixture were chosen so that a molar ratio of about 10 was set between the most and the least abundant mixture components, thus giving an indication of the concentration dynamics that can possibly be handled. After NMR analysis of the mixture, a semi-automated peak picking was performed to collect all 1 H-13 C HMBC connectivity data from the experimental spectrum. The resulting list of 1 H and 13 C chemical shift pairs was exported as a text file and used as input to create the HMBC-derived correlation network. A network of 267 nodes representing 1 H-13 C HMBC correlations was obtained (Fig. 3). The nodes sharing either a 1 H or a 13 C chemical shift value were connected by an edge. As the analyzed model mixture was composed of 3 molecules, the HMBC correlation network was expected to form 3 non-connected sub-networks representing each of the three molecules. However, only a big sub-network of 222 nodes was obtained, along with a small sub-network of 9 nodes and further sub-networks of 1 or 2 nodes. Non-expected interconnections were in fact due to the presence of wide, overlapping 1 H multiplet signals and of some identical chemical shift values between the three molecules. The major sub-network of 222 nodes was supposed to correspond to the three molecules of the model mixture. In an attempt to more efficiently isolate chemical shift clusters corresponding to individual molecule, a community detection algorithm was used. After the testing of several community detection approaches (Porter, Mason A et al), the Louvain method using the RBER method, which includes a resolution parameter, was found to be the most appropriate: i) the network is unweighted and undirected and the modularity based algorithms are the most appropriate to detect communities in this type of networks and ii) the resolution parameter allows the control of the communities size. By contrast, the default modularity algorithm tends to form very small communities that do not contain enough information for the identification of the metabolites. Applying this method with a resolution parameter set at 0.2, four main clusters (C1-C4) were obtained from the initial HMBC network (Fig. 3). The experimental 1 H-13 C HMBC correlations of each cluster were then compared to the theoretical 1 H-13 C HMBC correlations of the database metabolites. As a result, four lists of molecules having the highest likelihood to correspond to each cluster were obtained (Table 1). The proposed molecules were sorted according to their decreasing score value which represent the percentage of matching between experimental and theoretical HMBC correlations of each molecule. For cluster C1, the list of proposed molecules only contained flavan-3-ols (Table 1). The highest score was obtained for (+)-catechin (83 %) which was indeed one of the three constituents of the model mixture. For cluster C2, a list of dihydroflavonols was also obtained (Table 1), the highest score being attributed to taxifolin (94 %). For cluster C3, a list of 3 stilbenes and one tri-galloyl glucose was obtained (Table 1) with the highest score attributed to Eisorhapontin (54 %). This score was low as compared to the metabolites proposed for the two previous clusters, because the chemical structure of E-isohrapontin contains a glucose moiety and a methoxy group that do not share so many correlation with the aglycone moiety, so that the community detection algorithm affects the signals of the two moiety in separate communities. In addition, some correlations that would have allowed the connection of the communities were absent from the peak list. Cluster C4 was a small cluster containing just a few correlation values characteristic of a sugar moiety (Table 1). This cluster corresponded to the glucose part of E-isorhapontin. A second comparison was performed for the 20 top ranked molecules in each of the clusters C1-C4. In this case, all theoretical HMBC correlation values of the proposed molecules were compared to the whole experimental HMBC peak list. This step allowed to confirm the presence of taxifolin, (+)-catechin and E-isorhapontin in the model mixture with score values increased up to 97 % for taxifolin, 93 % for (+)-catechin and 70 % for the E-isorhapontin.

A final search was then performed to check the matching between the theoretical HSQC correlations of the top ranked molecules to the HSQC chemical shift values collected in the experimental spectrum. HSQC data inform about one-bond H-C correlations and therefore provide information highly complementary to HMBC data to improve the identification accuracy. Taxifolin, (+)-catechin and E-isorhapontin achieved all an HSQC score value of 100 %, thus providing an additional proof of their presence in the model mixture. et al., 2016). In view of this chemical diversity, this extract was considered to be a suitable model for the testing of the real life applicability of the method. A semi-automated peak picking was performed to collect all 1 H-13 C HMBC connectivity data from the experimental HMBC spectrum of the extract. The resulting peak list was used to create a HMBC correlation network in the same way as described for the model mixture. A network of 642 nodes was obtained, containing one large sub-network of 556 nodes and two small sub-networks of 11 and 8 nodes, respectively. The other sub-networks contained less than 5 nodes. Applying the Louvain community detection algorithm to this network with the resolution parameter set at 0.2, six clusters were obtained (C1-C6) (Fig. 4). The experimental HMBC data of these clusters were automatically compared to the theoretical HMBC data of the database metabolites. For Cluster C1, a list of highly similar stilbenes was obtained, the top ranked molecules of the list were E-resveratrol, E-piceid, E-isorhapontin, and E-astringinin with high score values of 76 %, 76 %, 73 % and 71 %, respectively. All these metabolites shared a E-resveratrol skeleton (5-[(1E)-2-(4-Hydroxyphenyl) ethenyl]-1,3-benzenediol), either linked or not to a glycosidic moiety, suggesting a highly likely presence of these molecules in the extract. For cluster C2, a list of sugars was proposed by the database (Table 2), with the highest score obtained for saccharose (63 %). A list of terpenoids and tocopherols was proposed for cluster C3, but the obtained score values were relatively poor, all below 50 % at this stage of the procedure. For cluster C4, a list of flavonoids was obtained, the top ranked molecules were (+)-catechin and taxifolin derivatives with (+)-catechin having the highest score (52 %). Cluster C5 was a small cluster of only 11 HMBC correlations for which a list of various sugars with low score values was obtained. This cluster contained in fact the HMBC correlation values of either a simple sugar or the sugar moiety of a glycosylated molecule. Cluster C6 contained only 8 correlations that did not allow prioritizing any relevant result. The search for clusters of HMBC correlations resulted in the identification of specific groups of molecules. However, as the HMBC fingerprint of the molecules contained in the extract were in some cases divided between several clusters, a second database search was performed for the 20 top ranked molecules of each cluster to compare their theoretical HMBC correlations to the total experimental HMBC correlations. As a result, a new list of molecules was obtained with improved scores, mainly containing stilbenes, flavonoids and sugars. The new scores of the top ranked molecules of cluster C1 increased, putting the glycosylated molecules in the top of the list with score values of 95 % for E-isorhapontin, 94 % for E-piceid, 91 % for E-astringin, and 92 % for E-rhaponticin. The scores of stilbene aglycones increased up to 81 % for E-resveratrol, 82 % for E-astringinin and 79 % E-rhapontigenin. The score of saccharose identified from cluster C2 also increased up to 74 %, and the score of 3,4dimethoxyphenyl-β-D-glucopyranoside, that was initially ranked at 12 th position in the cluster C2, arrived after the second step with a score value of 84 %. For cluster C3, the scores of the proposed molecules did not increased significantly after the second step of the database search, all remaining below 60 %. The score value of (+)-catechin identified in cluster C4 was increased up to 81 %, and other flavonoids including taxifolin, dihydrokaempferol, and taxifolin 3'-O-β-D-glucopyranoside arrived also in the top ranked molecules with score values of 87 %, 86 % and 85 %, respectively.

A final search was then performed to check the matching between the theoretical HSQC correlations of the top ranked molecules to the HSQC chemical shift values collected in the experimental spectrum. The final list of proposed molecules is given in Table 2, with fourteen molecules achieving a score value of 100 %.

In summary, among the top ranked molecules obtained over the three step of our identification procedure, fourteen molecules including E-isorhapontin, E-astringin, E-piceid, saccharose, (+)-catechin, taxifolin, taxifolin 3'-O-β-D-glucopyranoside, 3,4-dimethoxyphenylβ-D-glucopyranoside, E-resveratrol, E-piceatannol, E-rhapontigenin, E-rhaponticin, maltose, and dihydrokaempferol were putatively identified in the methanol bark extract of P. abies (Table 2). In order to unambiguously validate the presence of these molecules in the extract, their 1D and 2D NMR data were manually checked by going back to the experimental 1D and 2D NMR spectra. The presence of eight compounds was fully confirmed (Fig. 6), namely Eisorhapontin, E-astringin, E-piceid, saccharose, (+)-catechin, taxifolin, taxifolin 3'-O-β-Dglucopyranoside, and 3,4-dimethoxyphenyl-β-D-glucopyranoside. E-resveratrol was easily eliminated from the list of hits because no HMBC correlation corresponding to the highly characteristic symmetrical part of the 1,3-benzenediol moiety of the resveratrol carbon skeleton was detected, a correlation that is always present in all the spectra of compounds bearing this structural feature and that is due to the existence of a strong 3J coupling (J. Am. Chem. Soc., 1969, 91 (17), pp 4940-4941, pas mieux…). Similarly, E-piceatannol and Erhapontigenin were also eliminated from the list of hits because no HMBC correlation corresponding to the symmetrical part of the 1,3-benzenediol moiety was detected by NMR.

Concerning E-rhaponticin, which is a positional isomer of E-isorhapontin, the manual checking of experimental NMR data for the validation of E-isorhapontin enabled to discard Erhaponticin from the list of hits thanks to three main indices (Fig. 5): First, we confirmed via the interpretation of the HMBC spectrum of the crude extract that the methoxy group of the putatively identified E-rhaponticin was directly linked to C-1 at 147.7 ppm. Secondly, the proton at 6.99 ppm linked to C-4 at 119.9 ppm exhibited HMBC correlations with C-2 at 146.4 ppm, with C-3 at 114.8 ppm, with C-6 at 109.1 ppm, but not with C-1 at 147.7 ppm. An HMBC correlation was also observed between the proton at 7.14 linked to C-6 and C-1 at 147.7 ppm, indicating that C1 belongs to the same carbon skeleton as C-2, C-3, C-4, and C-6.

Therefore C-1 and C-4 were placed in a para position from each other. ∑Thirdly, the experimental chemical shift value of C-3 was 114.8 ppm, a value which is very close to the predicted value of the same carbon position for E-isorhapontin at 115.0 ppm, while the same carbon position of the E-rhaponticin is predicted at 110.6 ppm. Therefore the validation of the structure of E-isorhapontin lead to discard the possible presence of its very similar isomer Erhaponticin in the methanol bark extract of P. abies. For dihydrokaempferol and maltose, the experimental chemical shift values and 2D NMR connectivity data were not all confirmed after manual checking, therefore both compounds were also discarded from the list of hits. values: H4→C2 (6.99 ppm→146.4 ppm); H4→C3 (6.99 ppm→114.8 ppm), H4→C6 (6.99 ppm, 13 C 109.1 ppm), H6→C4 (7.14 ppm→119.9 ppm), H6→C1 (7.14 ppm→147.7 ppm), H6→C2 (7.14 ppm→146.4 ppm), H8→C1 (3.9 ppm→ 147.7ppm).

Conclusion

A new method exploiting HMBC and HSQC correlation data was developed in order to assist and accelerate the identification of natural metabolites directly within mixtures, while increasing the identification accuracy thanks to the richness of the information brought by HMBC and HSQC spectra. To the best of our knowledge, this is the first application of the network theory for the exploitation of heteronuclear 2D NMR data in the context of natural mixture analysis. In this paper, the proof of concept of the method was illustrated by two samples, a simplified model mixture of three commercial compounds and a real life natural extract. In the crude methanol extract of the bark of Picea abies, the method allowed the accurate identification of 8 molecules, namely E-isorhapontin, E-astringin, E-piceid, saccharose, (+)-catechin, taxifolin, taxifolin 3'-O-β-Dglucopyranoside and 3,4-dimethoxyphenyl-O-β-D-glucopyranoside. This method achieves a better performance when applied to the analysis of simplified fractions with well resolved HMBC spectra. The performance of this method is indeed highly dependent on the HMBC spectrum quality because the latter determines the reliability of the initial peak picking step. The method allows an early stage chemical profiling of the natural extracts and it could be very complementary to MS-based dereplication by combining the different advantages of each method in order to cross validate the results and have more accurate insight into the chemical composition of the extracts.

  C NMR, 2D NMR spectra such as HMBC (Heteronuclear Multiple Bond Correlation spectroscopy) and HSQC (Heteronuclear Single Quantum Correlation spectroscopy) are extremely informative. Both provide important information about atom connectivity. Such spectra have also attracted an increasing interest in recent dereplication and metabolomics studies (Bingol & Bruschweiler, 2013 + ref pauli hsqc mask actae triterpenes)

Figure

  Figure 1| Construction of the theoretical 2D NMR database from the predicted 1 H and 13 C

Figure 2 |

 2 Figure 2 | HMBC atom networking workflow for metabolite mixture identification.

Figure 3 |

 3 Figure 3 | HMBC-network obtained after NMR analysis of the model mixture and

Figure

  Figure 4| HMBC-network obtained after NMR analysis of the bark extract of P. abies. Colors

Figure 5 |

 5 Figure 5 | Experimental HMBC correlations observed for E-isorhapontin. Chemical shift

Figure 6 |

 6 Figure 6 | Chemical structures of the eight confirmed metabolites in the P. abies bark extract.

Table 2 .

 2 Proposal list for the P. abies bark extract.

	Putatively identified metabolites	Score 1 (HMBC	Score 2 (HMBC	Score	Validation
		cluster only)	full spectrum)	HSQC
	C1 E-resveratrol E-piceid	76 76	81 94	100 100	eliminated* confirmed
	E-isorhapontin	73	95	100	confirmed
	E-piceatannol (astringinin)	71	82	100	eliminated*
	E-rhapontigenin	69	79	100	eliminated*
	E-pterostilbene	68	73	91	eliminated
	E-astringin	66	91	100	confirmed
	7,3'-dihydroxy-4'-methoxy-	66	66	75	eliminated
	isoflavone
	rhaponticin	65	92	100	eliminated*
	E-resveratrol-3-O-D-glucuronide	64	79	93	eliminated

C NMR chemical shifts) were generated using a locally written Python program. For each